From ab8bb4f914b1ad19e83e290e7688b30d980fd17a Mon Sep 17 00:00:00 2001 From: Chunxiao Li Date: Mon, 16 Dec 2019 23:50:02 +0800 Subject: [PATCH] Delete NRLMSISE-00-checkpoint.ipynb --- .../NRLMSISE-00-checkpoint.ipynb | 1334 ----------------- 1 file changed, 1334 deletions(-) delete mode 100644 .ipynb_checkpoints/NRLMSISE-00-checkpoint.ipynb diff --git a/.ipynb_checkpoints/NRLMSISE-00-checkpoint.ipynb b/.ipynb_checkpoints/NRLMSISE-00-checkpoint.ipynb deleted file mode 100644 index a500d60..0000000 --- a/.ipynb_checkpoints/NRLMSISE-00-checkpoint.ipynb +++ /dev/null @@ -1,1334 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# -------------------------------------------------------------------- #\n", - "# -------------------- NRLMSISE-00 MODEL 2001 --------------------- #\n", - "# -------------------------------------------------------------------- #\n", - "\n", - "import numpy as np\n", - "from scipy.interpolate import CubicSpline\n", - "from pyshtools.legendre import PLegendreA,PlmIndex\n", - "from astropy.time import Time\n", - "from datetime import datetime,timedelta\n", - "from os import getenv,path,remove\n", - "from urllib.request import urlretrieve\n", - "\n", - "# ------------------------------------------------------------------- #\n", - "# ------------------------- READ DATA BLOCK ------------------------- #\n", - "# ------------------------------------------------------------------- #\n", - "\n", - "def nrlmsis00_data():\n", - " '''\n", - " 从 nrlmsis00_data.npz 文件中读取 nrlmsis00 所需的数据块 \n", - " Usage: pt,pd,ps,pdl,ptm,pdm,ptl,pma,sam,pavgm = nrlmsis00_data()\n", - " Inputs: \n", - " None\n", - " Outputs: \n", - " pt: [float array] TEMPERATURE \n", - " pd: [2d float array] DENSITY of HE, O, N2, Total mass, O2, AR, H, N, HOT O\n", - " ps: [float array] S PARAM\n", - " pdl: [2d float array] TURBO \n", - " ptm: [float array]\n", - " pdm: [2d float array]\n", - " ptl: [2d float array]\n", - " pma: [2d float array] \n", - " sam: [float array] SEMIANNUAL MULT SAM\n", - " pavgm: [float array] MIDDLE ATMOSPHERE AVERAGES \n", - " ''' \n", - " data = np.load('nrlmsis00_data.npz')\n", - " pt,pd,ps,pdl = data['pt'],data['pd'],data['ps'],data['pdl']\n", - " ptm,pdm,ptl,pma = data['ptm'],data['pdm'],data['ptl'],data['pma']\n", - " sam,pavgm = data['sam'],data['pavgm']\n", - " return pt,pd,ps,pdl,ptm,pdm,ptl,pma,sam,pavgm\n", - "\n", - "# ------------------------------------------------------------------- #\n", - "# ------------------------------ TSELEC ----------------------------- #\n", - "# ------------------------------------------------------------------- #\n", - "\n", - "def tselec(switches):\n", - " # len(switches) is equal to 23\n", - " flags = {'sw':np.zeros(23),'swc':np.zeros(23)}\n", - " for i in range(23):\n", - " if i != 8:\n", - " if switches[i] == 1:\n", - " flags['sw'][i] = 1\n", - " else:\n", - " flags['sw'][i] = 0\n", - " if switches[i] > 0:\n", - " flags['swc'][i] = 1\n", - " else:\n", - " flags['swc'][i] = 0\n", - " else:\n", - " flags['sw'][i] = switches[i]\n", - " flags['swc'][i] = switches[i]\n", - " return flags \n", - "\n", - "# ------------------------------------------------------------------- #\n", - "# ------------------------------ GLATF ------------------------------ #\n", - "# ------------------------------------------------------------------- #\n", - "\n", - "def glatf(lat):\n", - " c2 = np.cos(2*np.deg2rad(lat))\n", - " gv = 980.616*(1 - 0.0026373*c2)\n", - " reff = 2*gv/(3.085462E-6 + 2.27E-9*c2)*1E-5\n", - " return gv,reff\n", - "\n", - "# ------------------------------------------------------------------- #\n", - "# ------------------------------ CCOR ------------------------------- #\n", - "# ------------------------------------------------------------------- #\n", - "\n", - "def ccor(alt,r,h1,zh):\n", - " '''\n", - " CHEMISTRY/DISSOCIATION CORRECTION FOR MSIS MODELS\n", - " ALT - altitude\n", - " R - target ratio\n", - " H1 - transition scale length\n", - " ZH - altitude of 1/2 R\n", - " '''\n", - " e = (alt - zh)/h1\n", - " if e > 70:\n", - " return 1\n", - " elif e < -70:\n", - " return np.exp(r)\n", - " else:\n", - " return np.exp(r/(1 + np.exp(e)))\n", - " \n", - "# ------------------------------------------------------------------- #\n", - "# ------------------------------ CCOR2 ------------------------------ #\n", - "# ------------------------------------------------------------------- #\n", - "\n", - "def ccor2(alt,r,h1,zh,h2):\n", - " '''\n", - " CHEMISTRY/DISSOCIATION CORRECTION FOR MSIS MODELS\n", - " ALT - altitude\n", - " R - target ratio\n", - " H1 - transition scale length 1\n", - " ZH - altitude of 1/2 R\n", - " H2 - transition scale length 2 \n", - " ''' \n", - " e1 = (alt - zh)/h1\n", - " e2 = (alt - zh)/h2\n", - " if e1 > 70 or e2 > 70:\n", - " return 1\n", - " if e1 < -70 and e2 < -70:\n", - " return np.exp(r)\n", - " ex1,ex2 = np.exp([e1,e2])\n", - " ccor2v = r/(1 + 0.5*(ex1 + ex2))\n", - " return np.exp(ccor2v)\n", - "\n", - "# ------------------------------------------------------------------- #\n", - "# ------------------------------- SCALH ----------------------------- #\n", - "# ------------------------------------------------------------------- #\n", - "\n", - "def scalh(alt,xm,temp,gsurf,re):\n", - " rgas = 831.4\n", - " g = rgas*temp/(gsurf/(1 + alt/re)**2*xm)\n", - " return g\n", - "\n", - "# ------------------------------------------------------------------- #\n", - "# -------------------------------- DNET ----------------------------- #\n", - "# ------------------------------------------------------------------- #\n", - "\n", - "def dnet(dd,dm,zhm,xmm,xm):\n", - " '''\n", - " TURBOPAUSE CORRECTION FOR MSIS MODELS\n", - " Root mean density\n", - " DD - diffusive density\n", - " DM - full mixed density\n", - " ZHM - transition scale length\n", - " XMM - full mixed molecular weight\n", - " XM - species molecular weight\n", - " DNET - combined density\n", - " ''' \n", - " a = zhm/(xmm - xm)\n", - " if not (dm > 0 and dd > 0):\n", - " print('dnet log error {0:.1f} {1:.1f} {2:.1f}'.format(dm,dd,xm))\n", - " if dd == 0 and dm == 0: dd = 1\n", - " if dm == 0: return dd\n", - " if dd == 0: return dm\n", - " ylog = a*np.log(dm/dd)\n", - " if ylog < -10: return dd\n", - " if ylog > 10: return dm\n", - " a = dd*(1 + np.exp(ylog))**(1/a)\n", - " return a\n", - "\n", - "# ------------------------------------------------------------------- #\n", - "# -------------------------------- ZETA ----------------------------- #\n", - "# ------------------------------------------------------------------- #\n", - "\n", - "def zeta(zz,zl,re): \n", - " return (zz - zl)*(re + zl)/(re + zz)\n", - "\n", - "# ------------------------------------------------------------------- #\n", - "# ------------------------------- DENSM ----------------------------- #\n", - "# ------------------------------------------------------------------- #\n", - "\n", - "def densm(alt, d0, xm, tz, zn3, tn3, tgn3, zn2, tn2, tgn2,gsurf,re):\n", - " # Calculate Temperature and Density Profiles for lower atmos.\n", - " # call zeta\n", - " rgas = 831.4\n", - " densm_tmp = d0\n", - " tz_tmp = tz\n", - " mn3,mn2 = len(zn3),len(zn2)\n", - "\n", - " if alt > zn2[0]:\n", - " if xm == 0: \n", - " densm_tmp = tz\n", - " return densm_tmp,tz_tmp \n", - " else:\n", - " densm_tmp = d0\n", - " return densm_tmp,tz_tmp \n", - " \n", - " # STRATOSPHERE/MESOSPHERE TEMPERATURE\n", - " if alt > zn2[mn2-1]:\n", - " z = alt\n", - " else:\n", - " z = zn2[mn2-1]\n", - " mn = mn2\n", - " xs,ys = [np.zeros(mn) for i in range(2)]\n", - " z1,z2 = zn2[0],zn2[mn-1]\n", - " t1,t2=tn2[0],tn2[mn-1]\n", - " zg,zgdif = zeta(z,z1,re),zeta(z2,z1,re)\n", - " \n", - " # set up spline nodes\n", - " for k in range(mn):\n", - " xs[k] = zeta(zn2[k],z1,re)/zgdif\n", - " ys[k] = 1/tn2[k]\n", - " yd1 = -tgn2[0]/t1**2*zgdif\n", - " yd2 = -tgn2[1]/t2**2*zgdif*((re + z2)/(re + z1))**2\n", - "\n", - " # calculate spline coefficients\n", - " cs = CubicSpline(xs,ys,bc_type=((1,yd1),(1,yd2))) \n", - " x = zg/zgdif\n", - " y = cs(x)\n", - "\n", - " # temperature at altitude\n", - " tz_tmp = 1/y\n", - " if xm != 0:\n", - " # calaculate stratosphere / mesospehere density\n", - " glb = gsurf/(1 + z1/re)**2\n", - " gamm = xm*glb*zgdif/rgas\n", - " \n", - " # Integrate temperature profile\n", - " yi = cs.integrate(xs[0],x)\n", - " expl = gamm*yi\n", - " if expl > 50:\n", - " expl = 50\n", - " # Density at altitude\n", - " densm_tmp = densm_tmp*(t1/tz_tmp)*np.exp(-expl)\n", - " if alt > zn3[0]:\n", - " if xm == 0:\n", - " densm_tmp = tz_tmp\n", - " return densm_tmp,tz_tmp \n", - " else:\n", - " return densm_tmp,tz_tmp \n", - "\n", - " # troposhere / stratosphere temperature\n", - " z = alt\n", - " mn = mn3\n", - " xs,ys = [np.zeros(mn) for i in range(2)]\n", - " z1,z2 = zn3[0],zn3[mn-1]\n", - " t1,t2 = tn3[0],tn3[mn-1]\n", - " zg,zgdif = zeta(z,z1,re),zeta(z2,z1,re)\n", - "\n", - " # set up spline nodes\n", - " for k in range(mn):\n", - " xs[k] = zeta(zn3[k],z1,re)/zgdif\n", - " ys[k] = 1/tn3[k]\n", - " yd1 = -tgn3[0]/t1**2*zgdif\n", - " yd2 = -tgn3[1]/t2**2*zgdif*((re+z2)/(re+z1))**2\n", - "\n", - " # calculate spline coefficients\n", - " cs = CubicSpline(xs,ys,bc_type=((1,yd1),(1,yd2))) \n", - " x = zg/zgdif\n", - " y = cs(x)\n", - "\n", - " # temperature at altitude\n", - " tz_tmp = 1/y\n", - " if xm != 0:\n", - " # calaculate tropospheric / stratosphere density\n", - " glb = gsurf/(1 + z1/re)**2\n", - " gamm = xm*glb*zgdif/rgas\n", - " \n", - " # Integrate temperature profile\n", - " yi = cs.integrate(xs[0],x)\n", - " expl = gamm*yi\n", - " if expl > 50: expl = 50\n", - " # Density at altitude\n", - " densm_tmp = densm_tmp*(t1/tz_tmp)*np.exp(-expl)\n", - "\n", - " if xm == 0:\n", - " densm_tmp = tz_tmp\n", - " return densm_tmp,tz_tmp\n", - " else:\n", - " return densm_tmp,tz_tmp\n", - " \n", - "# ------------------------------------------------------------------- #\n", - "# ------------------------------- DENSU ----------------------------- #\n", - "# ------------------------------------------------------------------- #\n", - "\n", - "def densu (alt,dlb,tinf,tlb,xm,alpha,tz,zlb,s2,zn1,tn1,tgn1,gsurf,re):\n", - " # Calculate Temperature and Density Profiles for MSIS models\n", - " # New lower thermo polynomial\n", - " # call: zeta, \n", - "\n", - " rgas = 831.4\n", - " densu_tmp = 1\n", - " mn1 = len(zn1)\n", - " # joining altitudes of Bates and spline\n", - " za = zn1[0]\n", - " if alt > za:\n", - " z = alt\n", - " else:\n", - " z = za\n", - " # geopotential altitude difference from ZLB\n", - " zg2 = zeta(z,zlb,re)\n", - "\n", - " # Bates temperature\n", - " tt = tinf - (tinf - tlb)*np.exp(-s2*zg2)\n", - " ta = tz = tt\n", - " densu_tmp = tz_tmp = tz\n", - "\n", - " if alt < za:\n", - " # calculate temperature below ZA\n", - " # temperature gradient at ZA from Bates profile\n", - " dta = (tinf - ta)*s2*((re + zlb)/(re + za))**2\n", - " tgn1[0],tn1[0] = dta,ta\n", - " if alt > zn1[mn1-1]:\n", - " z = alt\n", - " else:\n", - " z = zn1[mn1-1]\n", - " mn = mn1\n", - " xs,ys = [np.zeros(mn) for i in range(2)]\n", - " z1,z2 = zn1[0],zn1[mn-1]\n", - " t1,t2 = tn1[0],tn1[mn-1]\n", - " # geopotental difference from z1\n", - " zg,zgdif = zeta(z,z1,re),zeta(z2,z1,re)\n", - " # set up spline nodes\n", - " for k in range(mn):\n", - " xs[k] = zeta(zn1[k],z1,re)/zgdif\n", - " ys[k] = 1/tn1[k]\n", - " # end node derivatives\n", - " yd1 = -tgn1[0]/t1**2*zgdif\n", - " yd2 = -tgn1[1]/t2**2*zgdif*((re + z2)/(re + z1))**2\n", - " # calculate spline coefficients\n", - " cs = CubicSpline(xs,ys,bc_type=((1,yd1),(1,yd2))) \n", - " x = zg/zgdif\n", - " y = cs(x)\n", - " # temperature at altitude\n", - " tz_tmp = 1/y\n", - " densu_tmp = tz_tmp\n", - " if xm == 0: return densu_tmp,tz_tmp\n", - " \n", - " # calculate density above za\n", - " glb = gsurf/(1 + zlb/re)**2\n", - " gamma = xm*glb/(s2*rgas*tinf)\n", - " expl = np.exp(-s2*gamma*zg2)\n", - " if expl > 50: expl = 50\n", - " if tt <= 0: expl = 50 \n", - "\n", - " # density at altitude\n", - " densa = dlb*(tlb/tt)**(1 + alpha + gamma)*expl\n", - " densu_tmp = densa\n", - " if alt >= za: return densu_tmp,tz_tmp\n", - " \n", - " # calculate density below za\n", - " glb = gsurf/(1 + z1/re)**2\n", - " gamm = xm*glb*zgdif/rgas\n", - "\n", - " # integrate spline temperatures\n", - " yi = cs.integrate(xs[0],x)\n", - " expl = gamm*yi\n", - " if expl > 50: expl = 50\n", - " if tz_tmp <= 0: expl = 50\n", - "\n", - " # density at altitude\n", - " densu_tmp = densu_tmp*(t1/tz_tmp)**(1 + alpha)*np.exp(-expl)\n", - " return densu_tmp,tz_tmp\n", - "\n", - "# ------------------------------------------------------------------- #\n", - "# --------------- 3hr Magnetic activity functions ------------------- #\n", - "# ------------------------------------------------------------------- #\n", - "\n", - "# Eq. A24d\n", - "def g0(a,p):\n", - " return (a - 4 + (p[25] - 1)*(a - 4 + (np.exp(-np.abs(p[24])*(a - 4)) - 1) / np.abs(p[24])))\n", - "\n", - "# Eq. A24c\n", - "def sumex(ex):\n", - " return (1 + (1 - ex**19)/(1 - ex)*ex**0.5)\n", - "\n", - "# Eq. A24a\n", - "def sg0(ex,p,ap):\n", - " # call sumex, g0\n", - " return (g0(ap[1],p) + g0(ap[2],p)*ex + g0(ap[3],p)*ex**2 + \\\n", - " g0(ap[4],p)*ex**3 + (g0(ap[5],p)*ex**4 + \\\n", - " g0(ap[6],p)*ex**12)*(1-ex**8)/(1-ex))/sumex(ex)\n", - "\n", - "# ------------------------------------------------------------------- #\n", - "# ------------------ Associated Legendre polynomials ---------------- #\n", - "# ------------------------------------------------------------------- #\n", - "\n", - "def lengendre(g_lat,lmax = 8):\n", - " # Index of PLegendreA_x can be calculated by PlmIndex(l,m)\n", - " x = np.sin(np.deg2rad(g_lat))\n", - " PLegendreA_x = PLegendreA(lmax,x)\n", - " return PLegendreA_x\n", - "\n", - "# ------------------------------------------------------------------- #\n", - "# ------------------------------- GLOBE7 ---------------------------- #\n", - "# ------------------------------------------------------------------- #\n", - "\n", - "def globe7(p,inputp,flags):\n", - " # CALCULATE G(L) FUNCTION \n", - " # Upper Thermosphere Parameters\n", - " # call: lengendre,sg0\n", - " t = np.zeros(15)\n", - " sr = 7.2722E-5\n", - " dr = 1.72142E-2\n", - " hr = 0.2618\n", - " \n", - " apdf = 0\n", - " apt = np.zeros(4)\n", - " tloc = inputp['lst']\n", - "\n", - " if not (flags['sw'][6]==0 and flags['sw'][7]==0 and flags['sw'][13]==0):\n", - " stloc,ctloc = np.sin(hr*tloc),np.cos(hr*tloc)\n", - " s2tloc,c2tloc = np.sin(2*hr*tloc),np.cos(2*hr*tloc)\n", - " s3tloc,c3tloc = np.sin(3*hr*tloc),np.cos(3*hr*tloc)\n", - " cd32 = np.cos(dr*(inputp['doy'] - p[31]))\n", - " cd18 = np.cos(2*dr*(inputp['doy'] - p[17]))\n", - " cd14 = np.cos(dr*(inputp['doy'] - p[13]))\n", - " cd39 = np.cos(2*dr*(inputp['doy'] - p[38]))\n", - "\n", - " # F10.7 EFFECT \n", - " df = inputp['f107'] - inputp['f107A']\n", - " dfa = inputp['f107A'] - 150\n", - " t[0] = p[19]*df*(1 + p[59]*dfa) + p[20]*df**2 + p[21]*dfa + p[29]*dfa**2\n", - " f1 = 1 + (p[47]*dfa + p[19]*df + p[20]*df**2)*flags['swc'][0]\n", - " f2 = 1 + (p[49]*dfa + p[19]*df + p[20]*df**2)*flags['swc'][0]\n", - " \n", - " plg = lengendre(inputp['g_lat'])\n", - "\n", - " # TIME INDEPENDENT \n", - " t[1] = p[1]*plg[3] + p[2]*plg[10] + p[22]*plg[21] + p[14]*plg[3]*dfa*flags['swc'][0] + p[26]*plg[1]\n", - " \n", - " # SYMMETRICAL ANNUAL \n", - " t[2] = p[18]*cd32\n", - "\n", - " # SYMMETRICAL SEMIANNUAL\n", - " t[3] = (p[15] + p[16]*plg[3])*cd18\n", - "\n", - " # ASYMMETRICAL ANNUAL\n", - " t[4] = f1*(p[9]*plg[1] + p[10]*plg[6])*cd14\n", - "\n", - " # ASYMMETRICAL SEMIANNUAL \n", - " t[5] = p[37]*plg[1]*cd39\n", - " \n", - " # DIURNAL \n", - " if flags['sw'][6]:\n", - " t71 = p[11]*plg[4]*cd14*flags['swc'][4]\n", - " t72 = p[12]*plg[4]*cd14*flags['swc'][4]\n", - " t[6] = f2*((p[3]*plg[2] + p[4]*plg[7] + p[27]*plg[16] + t71) * ctloc + (p[6]*plg[2] + p[7]*plg[7] + p[28]*plg[16] + t72)*stloc)\n", - " \n", - " # SEMIDIURNAL \n", - " if flags['sw'][7]:\n", - " t81 = (p[23]*plg[8] + p[35]*plg[17])*cd14*flags['swc'][4]\n", - " t82 = (p[33]*plg[8] + p[36]*plg[17])*cd14*flags['swc'][4]\n", - " t[7] = f2*((p[5]*plg[5] + p[41]*plg[12] + t81)*c2tloc +(p[8]*plg[5] + p[42]*plg[12] + t82)*s2tloc)\n", - "\n", - " # TERDIURNAL \n", - " if flags['sw'][13]:\n", - " t[13] = f2*((p[39]*plg[9] + (p[93]*plg[13] + p[46]*plg[24])*cd14*flags['swc'][4])*s3tloc + (p[40]*plg[9]+(p[94]*plg[13] + p[48]*plg[24])*cd14*flags['swc'][4])*c3tloc)\n", - " \n", - " # magnetic activity based on daily ap \n", - " if flags['sw'][8] == -1:\n", - " ap = inputp['ap_a']\n", - " if p[51]!= 0:\n", - " exp1 = np.exp(-10800*np.abs(p[51])/(1 + p[138]*(45 - np.abs(inputp['g_lat']))))\n", - " if exp1 > 0.99999: exp1 = 0.99999\n", - " if p[24] < 1E-4: p[24] = 1E-4\n", - " apt[0] = sg0(exp1,p,ap)\n", - " # apt[1] = sg2(exp1,p,ap)\n", - " # apt[2] = sg0(exp2,p,ap)\n", - " # apt[3] = sg2(exp2,p,ap)\n", - "\n", - " if flags['sw'][8]:\n", - " t[8] = apt[0]*(p[50] + p[96]*plg[3] + p[54]*plg[10] + \\\n", - " (p[125]*plg[1] + p[126]*plg[6] + p[127]*plg[15])*cd14*flags['swc'][4] + \\\n", - " (p[128]*plg[2] + p[129]*plg[7] + p[130]*plg[16])*flags['swc'][6]*np.cos(hr*(tloc - p[131])))\n", - " else:\n", - " apd = inputp['ap'] - 4\n", - " p44 = p[43]\n", - " p45 = p[44]\n", - " if p44 < 0: p44 = 1E-5\n", - " apdf = apd + (p45 - 1)*(apd + (np.exp(-p44*apd) - 1)/p44)\n", - " if flags['sw'][8]:\n", - " t[8]=apdf*(p[32] + p[45]*plg[3] + p[34]*plg[10] + \\\n", - " (p[100]*plg[1] + p[101]*plg[6] + p[102]*plg[15])*cd14*flags['swc'][4] +\n", - " (p[121]*plg[2] + p[122]*plg[7] + p[123]*plg[16])*flags['swc'][6]*np.cos(hr*(tloc - p[124])))\n", - "\n", - " if flags['sw'][9] and inputp['g_long'] > -1000:\n", - " # longitudinal\n", - " if flags['sw'][10]:\n", - " t[10] = (1 + p[80]*dfa*flags['swc'][0])*((p[64]*plg[4] + p[65]*plg[11] + p[66]*plg[22]\\\n", - " + p[103]*plg[2] + p[104]*plg[7] + p[105]*plg[16]\\\n", - " + flags['swc'][4]*(p[109]*plg[2] + p[110]*plg[7] + p[111]*plg[16])*cd14)*np.cos(np.deg2rad(inputp['g_long'])) \\\n", - " +(p[90]*plg[4]+p[91]*plg[11]+p[92]*plg[22] + p[106]*plg[2]+p[107]*plg[7]+p[108]*plg[16]\\\n", - " + flags['swc'][4]*(p[112]*plg[2] + p[113]*plg[7] + p[114]*plg[16])*cd14)*np.sin(np.deg2rad(inputp['g_long'])))\n", - "\n", - " # ut and mixed ut, longitude \n", - " if flags['sw'][11]:\n", - " t[11]=(1 + p[95]*plg[1])*(1 + p[81]*dfa*flags['swc'][0])*\\\n", - " (1 + p[119]*plg[1]*flags['swc'][4]*cd14)*\\\n", - " ((p[68]*plg[1] + p[69]*plg[6] + p[70]*plg[15])*np.cos(sr*(inputp['sec'] - p[71])))\n", - " t[11] += flags['swc'][10]*(p[76]*plg[8] + p[77]*plg[17] + p[78]*plg[30])*\\\n", - " np.cos(sr*(inputp['sec'] - p[79]) + 2*np.deg2rad(inputp['g_long']))*(1 + p[137]*dfa*flags['swc'][0])\n", - " \n", - " # ut, longitude magnetic activity \n", - " if flags['sw'][10]:\n", - " if flags['sw'][8] == -1:\n", - " if p[51]:\n", - " t[12] = apt[0]*flags['swc'][10]*(1 + p[132]*plg[1])*\\\n", - " ((p[52]*plg[4] + p[98]*plg[11] + p[67]*plg[22])* np.cos(np.deg2rad(inputp['g_long'] - p[97])))\\\n", - " + apt[0]*flags['swc'][10]*flags['swc'][4]*(p[133]*plg[2] + p[134]*plg[7] + p[135]*plg[16])*\\\n", - " cd14*np.cos(np.deg2rad(inputp['g_long'] - p[136])) + apt[0]*flags['swc'][11]* \\\n", - " (p[55]*plg[1] + p[56]*plg[6] + p[57]*plg[15])*np.cos(sr*(inputp['sec'] - p[58]))\n", - " else:\n", - " t[12] = apdf*flags['swc'][10]*(1 + p[120]*plg[1])*((p[60]*plg[4] + p[61]*plg[11] + p[62]*plg[22])*\\\n", - " np.cos(np.deg2rad(inputp['g_long']-p[63])))+apdf*flags['swc'][10]*flags['swc'][4]* \\\n", - " (p[115]*plg[2] + p[116]*plg[7] + p[117]*plg[16])* \\\n", - " cd14*np.cos(np.deg2rad(inputp['g_long'] - p[118])) \\\n", - " + apdf*flags['swc'][11]*(p[83]*plg[1] + p[84]*plg[6] + p[85]*plg[15])* np.cos(sr*(inputp['sec'] - p[75]))\n", - "\n", - " # parms not used: 82, 89, 99, 139-149 \n", - " tinf = p[30]\n", - " for i in range(14):\n", - " tinf = tinf + np.abs(flags['sw'][i])*t[i] \n", - " return tinf,[dfa,plg,ctloc,stloc,c2tloc,s2tloc,s3tloc,c3tloc,apdf,apt] \n", - "\n", - "# ------------------------------------------------------------------- #\n", - "# ------------------------------- GLOB7S ---------------------------- #\n", - "# ------------------------------------------------------------------- #\n", - "\n", - "def glob7s(p,inputp,flags,varli):\n", - " # VERSION OF GLOBE FOR LOWER ATMOSPHERE 10/26/99 \n", - " # call: lengendre,sg0\n", - " pset = 2\n", - " t = np.zeros(14)\n", - " dr = 1.72142E-2\n", - " [dfa,plg,ctloc,stloc,c2tloc,s2tloc,s3tloc,c3tloc,apdf,apt] = varli\n", - " \n", - " # confirm parameter set\n", - " if p[99] == 0: p[99] = pset\n", - " if p[99] != pset:\n", - " print(\"Wrong parameter set for glob7s\")\n", - " return -1\n", - "\n", - " for j in range(14):\n", - " t[j] = 0\n", - " cd32 = np.cos(dr*(inputp['doy'] - p[31]))\n", - " cd18 = np.cos(2*dr*(inputp['doy'] - p[17]))\n", - " cd14 = np.cos(dr*(inputp['doy'] - p[13]))\n", - " cd39 = np.cos(2*dr*(inputp['doy'] - p[38]))\n", - "\n", - " # F10.7 \n", - " t[0] = p[21]*dfa\n", - "\n", - " # time independent \n", - " t[1] = p[1]*plg[3] + p[2]*plg[10] + p[22]*plg[21] + p[26]*plg[1] + p[14]*plg[6] + p[59]*plg[15]\n", - "\n", - " # SYMMETRICAL ANNUAL \n", - " t[2] = (p[18] + p[47]*plg[3] + p[29]*plg[10])*cd32\n", - "\n", - " # SYMMETRICAL SEMIANNUAL \n", - " t[3] = (p[15] + p[16]*plg[3] + p[30]*plg[10])*cd18\n", - "\n", - " # ASYMMETRICAL ANNUAL \n", - " t[4] = (p[9]*plg[1] + p[10]*plg[6] + p[20]*plg[15])*cd14\n", - "\n", - " # ASYMMETRICAL SEMIANNUAL\n", - " t[5] = p[37]*plg[1]*cd39;\n", - "\n", - " # DIURNAL \n", - " if flags['sw'][6]:\n", - " t71 = p[11]*plg[4]*cd14*flags['swc'][4]\n", - " t72 = p[12]*plg[4]*cd14*flags['swc'][4]\n", - " t[6] = ((p[3]*plg[2] + p[4]*plg[7] + t71)*ctloc + (p[6]*plg[2] + p[7]*plg[7] + t72)*stloc) \n", - "\n", - " # SEMIDIURNAL\n", - " if flags['sw'][7]:\n", - " t81 = (p[23]*plg[8] + p[35]*plg[17])*cd14*flags['swc'][4]\n", - " t82 = (p[33]*plg[8] + p[36]*plg[17])*cd14*flags['swc'][4]\n", - " t[7] = ((p[5]*plg[5] + p[41]*plg[12] + t81)*c2tloc + (p[8]*plg[5] + p[42]*plg[12] + t82)*s2tloc)\n", - "\n", - " # TERDIURNAL \n", - " if flags['sw'][13]:\n", - " t[13] = p[39]*plg[9]*s3tloc + p[40]*plg[9]*c3tloc\n", - "\n", - " # MAGNETIC ACTIVITY\n", - " if flags['sw'][8]:\n", - " if flags['sw'][8]==1:\n", - " t[8] = apdf * (p[32] + p[45]*plg[3]*flags['swc'][1])\n", - " if flags['sw'][8]==-1:\n", - " t[8]=(p[50]*apt[0] + p[96]*plg[3]*apt[0]*flags['swc'][1])\n", - "\n", - " # LONGITUDINAL \n", - " if not (flags['sw'][9]==0 or flags['sw'][10]==0 or inputp['g_long']<=-1000):\n", - " t[10] = (1 + plg[1]*(p[80]*flags['swc'][4]*np.cos(dr*(inputp['doy'] - p[81]))\\\n", - " + p[85]*flags['swc'][5]*np.cos(2*dr*(inputp['doy'] - p[86])))\\\n", - " + p[83]*flags['swc'][2]*np.cos(dr*(inputp['doy'] - p[84]))\\\n", - " + p[87]*flags['swc'][3]*np.cos(2*dr*(inputp['doy'] - p[88])))\\\n", - " *((p[64]*plg[4] + p[65]*plg[11] + p[66]*plg[22]\\\n", - " + p[74]*plg[2] + p[75]*plg[7] + p[76]*plg[16])*np.cos(np.deg2rad(inputp['g_long']))\\\n", - " + (p[90]*plg[4] + p[91]*plg[11] + p[92]*plg[22]\\\n", - " + p[77]*plg[2] + p[78]*plg[7] + p[79]*plg[16])*np.sin(np.deg2rad(inputp['g_long'])))\n", - " \n", - " tt = 0\n", - " for i in range(14):\n", - " tt += np.abs(flags['sw'][i])*t[i]\n", - " return tt\n", - "\n", - "# ------------------------------------------------------------------- #\n", - "# ------------------------------- GTD7 ------------------------------ #\n", - "# ------------------------------------------------------------------- #\n", - " \n", - "def gtd7(inputp,switches):\n", - " tz = 0\n", - " zn3 = np.array([32.5,20.0,15.0,10.0,0.0])\n", - " zn2 = np.array([72.5,55.0,45.0,32.5])\n", - " zmix= 62.5\n", - " \n", - " output = {'d':{'He':0,'O':0,'N2':0,'O2':0,'AR':0,'RHO':0,'H':0,'N':0,'ANM O':0},\\\n", - " 't':{'TINF':0,'TG':0}}\n", - " \n", - " flags = tselec(switches)\n", - " \n", - " # Latitude variation of gravity (none for sw[1]=0) \n", - " xlat = inputp['g_lat']\n", - " if flags['sw'][1]==0: xlat = 45\n", - " gsurf,re = glatf(xlat)\n", - " pt,pd,ps,pdl,ptm,pdm,ptl,pma,sam,pavgm = nrlmsis00_data()\n", - " xmm = pdm[2,4]\n", - " \n", - " # THERMOSPHERE / MESOSPHERE (above zn2[0]) \n", - " if inputp['alt'] > zn2[0]:\n", - " altt = inputp['alt']\n", - " else:\n", - " altt = zn2[0]\n", - "\n", - " tmp = inputp['alt']\n", - " inputp['alt'] = altt\n", - " soutput,dm28,[meso_tn1,meso_tn2,meso_tn3,meso_tgn1,meso_tgn2,meso_tgn3],[dfa,plg,ctloc,stloc,c2tloc,s2tloc,s3tloc,c3tloc,apdf,apt] = gts7(inputp,flags,gsurf,re)\n", - " altt = inputp['alt']\n", - " inputp['alt'] = tmp\n", - " # metric adjustment \n", - " dm28m = dm28*1E6\n", - " output['t']['TINF'] = soutput['t']['TINF']\n", - " output['t']['TG'] = soutput['t']['TG']\n", - " if inputp['alt'] >= zn2[0]:\n", - " output['d'] = soutput['d']\n", - " return output\n", - " # LOWER MESOSPHERE/UPPER STRATOSPHERE (between zn3[0] and zn2[0])\n", - " # Temperature at nodes and gradients at end nodes\n", - " # Inverse temperature a linear function of spherical harmonics\n", - "\n", - " varli = [dfa,plg,ctloc,stloc,c2tloc,s2tloc,s3tloc,c3tloc,apdf,apt]\n", - "\n", - " meso_tgn2[0] = meso_tgn1[1]\n", - " meso_tn2[0] = meso_tn1[4]\n", - " meso_tn2[1] = pma[0,0]*pavgm[0]/(1-flags['sw'][19]*glob7s(pma[0], inputp, flags,varli))\n", - " meso_tn2[2] = pma[1,0]*pavgm[1]/(1-flags['sw'][19]*glob7s(pma[1], inputp, flags,varli))\n", - " meso_tn2[3] = pma[2,0]*pavgm[2]/(1-flags['sw'][19]*flags['sw'][21]*glob7s(pma[2], inputp, flags,varli))\n", - " meso_tgn2[1] = pavgm[8]*pma[9,0]*(1+flags['sw'][19]*flags['sw'][21]*glob7s(pma[9], inputp, flags,varli))*meso_tn2[3]*meso_tn2[3]/(pma[2,0]*pavgm[2])**2\n", - " meso_tn3[0] = meso_tn2[3]\n", - " \n", - " if inputp['alt'] <= zn3[0]:\n", - " # LOWER STRATOSPHERE AND TROPOSPHERE (below zn3[0])\n", - " # Temperature at nodes and gradients at end nodes\n", - " # Inverse temperature a linear function of spherical harmonics\n", - "\n", - " meso_tgn3[0] = meso_tgn2[1]\n", - " meso_tn3[1] = pma[3,0]*pavgm[3]/(1-flags['sw'][21]*glob7s(pma[3], inputp, flags,varli))\n", - " meso_tn3[2] = pma[4,0]*pavgm[4]/(1-flags['sw'][21]*glob7s(pma[4], inputp, flags,varli))\n", - " meso_tn3[3] = pma[5,0]*pavgm[5]/(1-flags['sw'][21]*glob7s(pma[5], inputp, flags,varli))\n", - " meso_tn3[4] = pma[6,0]*pavgm[6]/(1-flags['sw'][21]*glob7s(pma[6], inputp, flags,varli))\n", - " meso_tgn3[1] = pma[7,0]*pavgm[7]*(1+flags['sw'][21]*glob7s(pma[7], inputp, flags,varli)) *meso_tn3[4]*meso_tn3[4]/(pma[6,0]*pavgm[6])**2\n", - " # LINEAR TRANSITION TO FULL MIXING BELOW zn2[0] \n", - "\n", - " dmc = 0\n", - " if inputp['alt'] > zmix:\n", - " dmc = 1 - (zn2[0]-inputp['alt'])/(zn2[0] - zmix)\n", - " dz28 = soutput['d']['N2']\n", - " \n", - " # N2 density\n", - " dmr = soutput['d']['N2'] / dm28m - 1\n", - " output['d']['N2'],tz = densm(inputp['alt'],dm28m,xmm, tz, zn3, meso_tn3, meso_tgn3, zn2, meso_tn2, meso_tgn2,gsurf,re)\n", - " output['d']['N2'] = output['d']['N2'] * (1 + dmr*dmc)\n", - "\n", - " # HE density \n", - " dmr = soutput['d']['He'] / (dz28 * pdm[0,1]) - 1\n", - " output['d']['He'] = output['d']['N2'] * pdm[0,1] * (1 + dmr*dmc)\n", - "\n", - " # O density\n", - " output['d']['O'] = 0\n", - " output['d']['ANM O'] = 0\n", - "\n", - " # O2 density\n", - " dmr = soutput['d']['O2'] / (dz28 * pdm[3,1]) - 1\n", - " output['d']['O2'] = output['d']['N2'] * pdm[3,1] * (1 + dmr*dmc)\n", - "\n", - " # AR density \n", - " dmr = soutput['d']['AR'] / (dz28 * pdm[4,1]) - 1\n", - " output['d']['AR'] = output['d']['N2'] * pdm[4,1] * (1 + dmr*dmc)\n", - "\n", - " # Hydrogen density\n", - " output['d']['H'] = 0\n", - "\n", - " # Atomic nitrogen density \n", - " output['d']['N'] = 0\n", - "\n", - " # Total mass density \n", - " output['d']['RHO'] = 1.66E-24 * (4 * output['d']['He'] + 16 * output['d']['O'] + 28 * output['d']['N2']\\\n", - " + 32 * output['d']['O2'] + 40 * output['d']['AR'] + output['d']['H'] + 14 * output['d']['N'])\n", - "\n", - " output['d']['RHO'] = output['d']['RHO']/1000\n", - "\n", - " # temperature at altitude \n", - " dd,tz = densm(inputp['alt'], 1, 0, tz, zn3, meso_tn3, meso_tgn3, zn2, meso_tn2, meso_tgn2,gsurf,re)\n", - " output['t']['TG'] = tz\n", - " return output\n", - " \n", - "# ------------------------------------------------------------------- #\n", - "# ------------------------------- GTD7D ----------------------------- #\n", - "# ------------------------------------------------------------------- #\n", - "def gtd7d(inputp, flags):\n", - " output = gtd7(inputp, flags)\n", - " output['d']['RHO'] = 1.66E-24 * (4 * output['d']['He'] + 16 * output['d']['O'] + 28 * output['d']['N2']\\\n", - " + 32 * output['d']['O2'] + 40 * output['d']['AR'] + output['d']['H'] + 14 * output['d']['N'] + 16 * output['d']['ANM O'])\n", - "\n", - " output['d']['RHO'] = output['d']['RHO']/1000\n", - " return output\n", - "\n", - "\n", - "# ------------------------------------------------------------------- #\n", - "# ------------------------------- GTS7 ------------------------------ #\n", - "# ------------------------------------------------------------------- #\n", - " \n", - "def gts7(inputp,flags,gsurf,re):\n", - " # Thermospheric portion of NRLMSISE-00\n", - " # See GTD7 for more extensive comments\n", - " # alt > 72.5 km!\n", - " # call: nrlmsis00_data, globe7, densu\n", - " \n", - " output = {'d':{'He':0,'O':0,'N2':0,'O2':0,'AR':0,'RHO':0,'H':0,'N':0,'ANM O':0},\\\n", - " 't':{'TINF':0,'TG':0}}\n", - " tz = 0\n", - " dm28 = 0\n", - " meso_tn1,meso_tn3 = [np.zeros(5) for i in range(2)]\n", - " meso_tn2 = np.zeros(4)\n", - " meso_tgn1,meso_tgn2,meso_tgn3 = [np.zeros(2) for i in range(3)]\n", - " \n", - " zn1 = np.array([120.0, 110.0, 100.0, 90.0, 72.5])\n", - "\n", - " dr = 1.72142E-2\n", - " alpha = np.array([-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0])\n", - " altl = np.array([200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0])\n", - " pt,pd,ps,pdl,ptm,pdm,ptl,pma,sam,pavgm = nrlmsis00_data()\n", - " za = pdl[1,15]\n", - " zn1[0] = za\n", - " \n", - " # TINF VARIATIONS NOT IMPORTANT BELOW ZA OR ZN1(1)\n", - " if inputp['alt'] > zn1[0]:\n", - " tinf_tmp,varli = globe7(pt,inputp,flags)\n", - " tinf = ptm[0]*pt[0] * (1+flags['sw'][15]*tinf_tmp)\n", - " else:\n", - " tinf = ptm[0]*pt[0]\n", - " output['t']['TINF'] = tinf\n", - " \n", - " # GRADIENT VARIATIONS NOT IMPORTANT BELOW ZN1(5)\n", - " if inputp['alt'] > zn1[4]:\n", - " tinf_tmp,varli = globe7(ps,inputp,flags)\n", - " grad = ptm[3]*ps[0] * (1+flags['sw'][18]*tinf_tmp)\n", - " else:\n", - " grad = ptm[3]*ps[0]\n", - " tinf_tmp,varli = globe7(pd[3],inputp,flags) \n", - " tlb = ptm[1] * (1 + flags['sw'][16]*tinf_tmp)*pd[3,0]\n", - " s = grad/(tinf - tlb)\n", - " \n", - " # Lower thermosphere temp variations not significant for density above 300 km\n", - " if inputp['alt'] < 300:\n", - " meso_tn1[1] = ptm[6]*ptl[0,0]/(1.0-flags['sw'][17]*glob7s(ptl[0], inputp, flags,varli))\n", - " meso_tn1[2] = ptm[2]*ptl[1,0]/(1.0-flags['sw'][17]*glob7s(ptl[1], inputp, flags,varli))\n", - " meso_tn1[3] = ptm[7]*ptl[2,0]/(1.0-flags['sw'][17]*glob7s(ptl[2], inputp, flags,varli))\n", - " meso_tn1[4] = ptm[4]*ptl[3,0]/(1.0-flags['sw'][17]*flags['sw'][19]*glob7s(ptl[3], inputp, flags,varli))\n", - " meso_tgn1[1] = ptm[8]*pma[8,0]*(1.0+flags['sw'][17]*flags['sw'][19]*glob7s(pma[8], inputp, flags,varli))*meso_tn1[4]*meso_tn1[4]/(ptm[4]*ptl[3,0])**2\n", - " else:\n", - " meso_tn1[1]=ptm[6]*ptl[0,0]\n", - " meso_tn1[2]=ptm[2]*ptl[1,0]\n", - " meso_tn1[3]=ptm[7]*ptl[2,0]\n", - " meso_tn1[4]=ptm[4]*ptl[3,0]\n", - " meso_tgn1[1]=ptm[8]*pma[8,0]*meso_tn1[4]*meso_tn1[4]/(ptm[4]*ptl[3,0])**2\n", - " \n", - " # N2 variation factor at Zlb\n", - " tinf_tmp,varli = globe7(pd[2],inputp,flags)\n", - " g28 = flags['sw'][20]*tinf_tmp\n", - "\n", - " # VARIATION OF TURBOPAUSE HEIGHT\n", - " zhf = pdl[1,24]*(1+flags['sw'][4]*pdl[0,24]*np.sin(np.deg2rad(inputp['g_lat']))*np.cos(dr*(inputp['doy']-pt[13])))\n", - " output['t']['TINF'] = tinf\n", - " xmm = pdm[2,4]\n", - " z = inputp['alt']\n", - "\n", - " # N2 DENSITY\n", - " # Diffusive density at Zlb \n", - " db28 = pdm[2,0]*np.exp(g28)*pd[2,0]\n", - " # Diffusive density at Alt \n", - " output['d']['N2'],output['t']['TG'] = densu(z,db28,tinf,tlb,28,alpha[2],output['t']['TG'],ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re)\n", - " dd = output['d']['N2']\n", - " # Turbopause \n", - " zh28 = pdm[2,2]*zhf\n", - " zhm28 = pdm[2,3]*pdl[1,5] \n", - " xmd = 28 - xmm\n", - " # Mixed density at Zlb \n", - " b28,tz = densu(zh28,db28,tinf,tlb,xmd,(alpha[2]-1),tz,ptm[5],s, zn1,meso_tn1,meso_tgn1,gsurf,re)\n", - " if flags['sw'][14] and z <= altl[2]:\n", - " # Mixed density at Alt \n", - " dm28,tz = densu(z,b28,tinf,tlb,xmm,alpha[2],tz,ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re)\n", - " # Net density at Alt\n", - " output['d']['N2'] = dnet(output['d']['N2'],dm28,zhm28,xmm,28)\n", - " \n", - " # HE DENSITY\n", - " # Density variation factor at Zlb\n", - " tinf_tmp,varli = globe7(pd[0],inputp,flags)\n", - " g4 = flags['sw'][20]*tinf_tmp\n", - " # Diffusive density at Zlb \n", - " db04 = pdm[0,0]*np.exp(g4)*pd[0,0]\n", - " # Diffusive density at Alt \n", - " output['d']['He'],output['t']['TG'] = densu(z,db04,tinf,tlb, 4,alpha[0],output['t']['TG'],ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re)\n", - " dd = output['d']['He']\n", - " if flags['sw'][14] and z modified_time + timedelta(days=1):\n", - " remove(swfile)\n", - " print('Updating the space weather data',end=' ... ')\n", - " urlretrieve(url, swfile)\n", - " print('finished')\n", - " return swfile\n", - "\n", - "# ------------------------------------------------------------------- #\n", - "# ------------------------ READ SW DATA ----------------------------- #\n", - "# ------------------------------------------------------------------- #\n", - "\n", - "def read_sw(swfile):\n", - " # read the SPACE WEATHER DATA\n", - " sw_data = open(swfile,'r').readlines()\n", - " # read the SPACE WEATHER DATA\n", - " SW_OBS,SW_PRE = [],[]\n", - " flag1 = flag2 = 0\n", - " for line in sw_data:\n", - " if line.startswith('BEGIN OBSERVED'): \n", - " flag1 = 1\n", - " continue\n", - " if line.startswith('END OBSERVED'): flag1 = 0 \n", - " if flag1 == 1: \n", - " sw_p = line.split()\n", - " if len(sw_p) == 30:\n", - " del sw_p[24]\n", - " elif len(sw_p) == 31: \n", - " sw_p = np.delete(sw_p,[23,25]) \n", - " else: \n", - " sw_p = np.delete(sw_p,[23,24,25,27])\n", - " SW_OBS.append(sw_p)\n", - " \n", - " if line.startswith('BEGIN DAILY_PREDICTED'): \n", - " flag2 = 1\n", - " continue \n", - " if line.startswith('END DAILY_PREDICTED'): break \n", - " if flag2 == 1: SW_PRE.append(line.split()) \n", - " SW_OBS_PRE = np.vstack((np.array(SW_OBS),np.array(SW_PRE))) \n", - " # inverse sort\n", - " SW_OBS_PRE = np.flip(SW_OBS_PRE,0)\n", - " return SW_OBS_PRE\n", - "#---------------------------------------------------------- \n", - "def get_sw(SW_OBS_PRE,t_ymd,hour):\n", - " j = 0\n", - " for ymd in SW_OBS_PRE[:,:3]:\n", - " if np.array_equal(t_ymd,ymd): break\n", - " j+=1 \n", - " f107A,f107,ap = float(SW_OBS_PRE[j,27]),float(SW_OBS_PRE[j+1,26]),int(SW_OBS_PRE[j,22])\n", - " aph_tmp_b0 = SW_OBS_PRE[j,14:22] \n", - " i = int(np.floor_divide(hour,3))\n", - " ap_c = aph_tmp_b0[i]\n", - " aph_tmp_b1 = SW_OBS_PRE[j+1,14:22]\n", - " aph_tmp_b2 = SW_OBS_PRE[j+2,14:22]\n", - " aph_tmp_b3 = SW_OBS_PRE[j+3,14:22]\n", - " aph_tmp = np.hstack((aph_tmp_b3,aph_tmp_b2,aph_tmp_b1,aph_tmp_b0))[::-1].astype(np.float)\n", - " apc_index = 7-i\n", - " aph_c369 = aph_tmp[apc_index:apc_index+4]\n", - " aph_1233 = np.average(aph_tmp[apc_index+4:apc_index+12])\n", - " aph_3657 = np.average(aph_tmp[apc_index+12:apc_index+20])\n", - " aph = np.hstack((ap,aph_c369,aph_1233,aph_3657))\n", - " return f107A,f107,ap,aph\n", - "\n", - "# ------------------------------------------------------------------- #\n", - "# ------------------------------- OTHER ----------------------------- #\n", - "# ------------------------------------------------------------------- #\n", - "\n", - "def wraplon(lon):\n", - " if lon > 180:\n", - " lonwrap = lon - 360\n", - " else:\n", - " lonwrap = lon\n", - " return lonwrap \n", - "def hms2s(h,m,s):\n", - " return h*3.6E3 + m*60 + s\n", - "def hms2h(h,m,s):\n", - " return h + m/60 + s/3.6E3\n", - "\n", - "# ------------------------------------------------------------------- #\n", - "# ---------------------------- NRLMSISE00 --------------------------- #\n", - "# ------------------------------------------------------------------- #\n", - "\n", - "def nrlmsise00(t,lat,lon,alt,SW_OBS_PRE,o='Oxygen',s='NoAph'):\n", - " lon_wrap = wraplon(lon)\n", - " t_yday = t.yday.split(':')\n", - " t_ymd = t.iso.split()[0].split('-')\n", - " year,doy = int(t_yday[0]),int(t_yday[1])\n", - " sec = hms2s(int(t_yday[2]),int(t_yday[3]),float(t_yday[4]))\n", - " hour = hms2h(int(t_yday[2]),int(t_yday[3]),float(t_yday[4]))\n", - " lst = hour + lon/15\n", - " if alt > 80:\n", - " f107A,f107,ap,aph = get_sw(SW_OBS_PRE,t_ymd,hour)\n", - " else:\n", - " f107A,f107,ap,aph = 150,150,4,np.full(7,4)\n", - " inputp = {'doy':doy,'year':year,'sec':sec,'alt':alt,'g_lat':lat,'g_long':lon_wrap,'lst':lst,\\\n", - " 'f107A':f107A,'f107':f107,'ap':ap,'ap_a':aph}\n", - " \n", - " switches = np.ones(23)\n", - " if s is 'Aph':\n", - " switches[8] = -1 # -1 表示使用 3h 的地磁指数\n", - " \n", - " if o is 'Oxygen':\n", - " output = gtd7d(inputp,switches)\n", - " elif o is 'NoOxygen':\n", - " output = gtd7(inputp,switches)\n", - " else:\n", - " raise Exception(\"'{}' should be either 'Oxygen' or 'NoOxygen'\".format(o))\n", - " inputp['g_long'] = lon \n", - " return inputp,output " - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "The existing space weather data is already up to date\n" - ] - } - ], - "source": [ - "from pyatmos.msise import download_sw,read_sw\n", - "from pyatmos.atmosclasses import Coordinate\n", - "\n", - "# Download or update the space weather file from www.celestrak.com\n", - "swfile = download_sw() \n", - "# Read the space weather data\n", - "sw_obs_pre = read_sw(swfile) " - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "# Set a specific time and location\n", - "t = '2015-10-05 03:00:00' # time(UTC)\n", - "lat,lon = 25,102 # latitude and longitude [degree]\n", - "alt = 70 # altitude [km]\n", - "\n", - "# Initialize a coordinate instance by a space-time point\n", - "st = Coordinate(t,lat,lon,alt)" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [ - { - "ename": "FileNotFoundError", - "evalue": "[Errno 2] No such file or directory: '../data/nrlmsis00_data.npz'", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mFileNotFoundError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mpara_input\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpara_output\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mst\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnrlmsise00\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msw_obs_pre\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpara_input\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpara_output\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/Downloads/ATMOS/pyatmos/atmosclasses/coordinate.py\u001b[0m in \u001b[0;36mnrlmsise00\u001b[0;34m(self, sw_obs_pre, omode, aphmode)\u001b[0m\n\u001b[1;32m 104\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m'd'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m'He'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;36m74934329990.0412\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'O'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;36m71368139.39199762\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'N2'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;36m104.72048033793158\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'O2'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;36m0.09392848471935447\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'AR'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;36m1.3231114543012155e-07\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'RHO'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;36m8.914971667362366e-16\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'H'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;36m207405192640.34592\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'N'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;36m3785341.821909535\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'ANM O'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;36m1794317839.638502\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m't'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m'TINF'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;36m646.8157488121493\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'TG'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;36m646.8157488108872\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 105\u001b[0m '''\n\u001b[0;32m--> 106\u001b[0;31m \u001b[0mpara_input\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpara_output\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnrlmsise00\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mTime\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlat\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlon\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0malt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0msw_obs_pre\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0momode\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0maphmode\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 107\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mpara_input\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpara_output\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 108\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/Downloads/ATMOS/pyatmos/msise/nrlmsise00.py\u001b[0m in \u001b[0;36mnrlmsise00\u001b[0;34m(t, lat, lon, alt, SW_OBS_PRE, omode, aphmode)\u001b[0m\n\u001b[1;32m 1040\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1041\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0momode\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;34m'Oxygen'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1042\u001b[0;31m \u001b[0moutput\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgtd7d\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minputp\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mswitches\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1043\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0momode\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;34m'NoOxygen'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1044\u001b[0m \u001b[0moutput\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgtd7\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minputp\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mswitches\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/Downloads/ATMOS/pyatmos/msise/nrlmsise00.py\u001b[0m in \u001b[0;36mgtd7d\u001b[0;34m(inputp, flags)\u001b[0m\n\u001b[1;32m 722\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 723\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mgtd7d\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minputp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 724\u001b[0;31m \u001b[0moutput\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgtd7\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minputp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 725\u001b[0m output['d']['RHO'] = 1.66E-24 * (4 * output['d']['He'] + 16 * output['d']['O'] + 28 * output['d']['N2']\\\n\u001b[1;32m 726\u001b[0m + 32 * output['d']['O2'] + 40 * output['d']['AR'] + output['d']['H'] + 14 * output['d']['N'] + 16 * output['d']['ANM O'])\n", - "\u001b[0;32m~/Downloads/ATMOS/pyatmos/msise/nrlmsise00.py\u001b[0m in \u001b[0;36mgtd7\u001b[0;34m(inputp, switches)\u001b[0m\n\u001b[1;32m 626\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'sw'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mxlat\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m45\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 627\u001b[0m \u001b[0mgsurf\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mre\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mglatf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mxlat\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 628\u001b[0;31m \u001b[0mpt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpd\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mps\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpdl\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mptm\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpdm\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mptl\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpma\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0msam\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpavgm\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnrlmsis00_data\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 629\u001b[0m \u001b[0mxmm\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpdm\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 630\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/Downloads/ATMOS/pyatmos/msise/nrlmsise00.py\u001b[0m in \u001b[0;36mnrlmsis00_data\u001b[0;34m()\u001b[0m\n\u001b[1;32m 79\u001b[0m \u001b[0mpavgm\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mfloat\u001b[0m \u001b[0marray\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0mMIDDLE\u001b[0m \u001b[0mATMOSPHERE\u001b[0m \u001b[0mAVERAGES\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 80\u001b[0m ''' \n\u001b[0;32m---> 81\u001b[0;31m \u001b[0mdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mload\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'../data/nrlmsis00_data.npz'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 82\u001b[0m \u001b[0mpt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpd\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mps\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpdl\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'pt'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'pd'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'ps'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'pdl'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 83\u001b[0m \u001b[0mptm\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpdm\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mptl\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'ptm'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'pdm'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'ptl'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'pma'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/anaconda3/envs/py37/lib/python3.7/site-packages/numpy/lib/npyio.py\u001b[0m in \u001b[0;36mload\u001b[0;34m(file, mmap_mode, allow_pickle, fix_imports, encoding)\u001b[0m\n\u001b[1;32m 426\u001b[0m \u001b[0mown_fid\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 427\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 428\u001b[0;31m \u001b[0mfid\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mos_fspath\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfile\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"rb\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 429\u001b[0m \u001b[0mown_fid\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 430\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: '../data/nrlmsis00_data.npz'" - ] - } - ], - "source": [ - "para_input,para_output = st.nrlmsise00(sw_obs_pre)\n", - "print(para_input)\n", - "print(para_output)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "t = '2019-08-20 23:10:59' \n", - "lat,lon,alt = 3,5,900 \n", - "st = Coordinate(t,lat,lon,alt)\n", - "para_input,para_output = st.nrlmsise00(sw_obs_pre,aphmode = 'Aph')\n", - "print(para_input)\n", - "print(para_output)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "纬度范围为（-90，90），经度范围为（0，360）或（-180，180）；高度单位：km；输出单位: /m^3,密度单位：kg/m^3\n", - "温度单位：K; 72.5km 以下，O、H 、N、ANMO 均为零；空间天气数据每 12h 更新一次。" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "测试结果\n", - "无aph 2015-10-05 03:00:00 25 102 70 Oxygen\n", - "matlab 1027.318 219.965 8.2335e-05\n", - "matlab(areo): 1027.318 219.965 8.2335e-05\n", - "C: 1027.318 219.965 8.2335e-05\n", - "python: 1027.318 219.965 8.2335e-05\n", - "\n", - "无aph 2004-07-08 10:30:50 -65 -120 100 Oxygen\n", - "matlab: 1027.318 192.587 4.5846e-07\n", - "matlab(aero): 1027.318 192.587 4.5846e-07\n", - "C: 1027.318 192.587 4.5846e-07\n", - "python: 1027.318 192.587 4.5846e-07\n", - "\n", - "有aph 2019-08-20 23:10:59 3 5 900 Oxygen\n", - "matlab: 640.742 640.742 8.7480e-16 \n", - "matlab(aero): 640.741 640.741 8.7482e-16\n", - "C: 640.741 640.741 8.7482e-16\n", - "python: 640.741 640.741 8.7482e-16\n", - "\n", - "有aph 2010-02-15 12:18:37 85 170 500 NoOxygen\n", - "matlab: 774.172 774.171 1.2708e-13\n", - "matlab(areo): 774.168 774.166 1.2707e-13\n", - "C: 774.168 774.166 1.2707e-13\n", - "python: 774.168 774.166 1.2707e-13" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import msise00\n", - "from datetime import datetime\n", - "\n", - "atmos = msise00.run(time=datetime(2018, 5, 17, 21), altkm=300, glat=55, glon=120)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "atmos" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from numpy.linalg import norm" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "norm?" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "st.nrlmsise00?" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.7.4" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -}