1335 lines
69 KiB
Plaintext
1335 lines
69 KiB
Plaintext
{
|
||
"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<altl[0]:\n",
|
||
" # Turbopause \n",
|
||
" zh04 = pdm[0,2]\n",
|
||
" # Mixed density at Zlb\n",
|
||
" b04,output['t']['TG'] = densu(zh04,db04,tinf,tlb,4-xmm,alpha[0]-1,output['t']['TG'],ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re)\n",
|
||
" # Mixed density at Alt\n",
|
||
" dm04,output['t']['TG'] = densu(z,b04,tinf,tlb,xmm,0,output['t']['TG'],ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re)\n",
|
||
" zhm04 = zhm28\n",
|
||
" # Net density at Alt\n",
|
||
" output['d']['He'] = dnet(output['d']['He'],dm04,zhm04,xmm,4)\n",
|
||
" # Correction to specified mixing ratio at ground \n",
|
||
" rl = np.log(b28*pdm[0,1]/b04)\n",
|
||
" zc04 = pdm[0,4]*pdl[1,0]\n",
|
||
" hc04 = pdm[0,5]*pdl[1,1]\n",
|
||
" # Net density corrected at Alt \n",
|
||
" output['d']['He'] = output['d']['He']*ccor(z,rl,hc04,zc04) \n",
|
||
" \n",
|
||
" # O DENSITY \n",
|
||
" # Density variation factor at Zlb \n",
|
||
" tinf_tmp,varli = globe7(pd[1],inputp,flags)\n",
|
||
" g16 = flags['sw'][20]*tinf_tmp\n",
|
||
" # Diffusive density at Zlb \n",
|
||
" db16 = pdm[1,0]*np.exp(g16)*pd[1,0]\n",
|
||
" # Diffusive density at Alt \n",
|
||
" output['d']['O'],output['t']['TG'] = densu(z,db16,tinf,tlb,16,alpha[1],output['t']['TG'],ptm[5],s, zn1,meso_tn1,meso_tgn1,gsurf,re)\n",
|
||
" dd = output['d']['O']\n",
|
||
" if flags['sw'][14] and z <= altl[1]:\n",
|
||
" # Turbopause \n",
|
||
" zh16 = pdm[1,2]\n",
|
||
" # Mixed density at Zlb \n",
|
||
" b16,output['t']['TG'] = densu(zh16,db16,tinf,tlb,16-xmm,alpha[1]-1, output['t']['TG'],ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re)\n",
|
||
" # Mixed density at Alt \n",
|
||
" dm16,output['t']['TG'] = densu(z,b16,tinf,tlb,xmm,0,output['t']['TG'],ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re)\n",
|
||
" zhm16 = zhm28\n",
|
||
" # Net density at Alt \n",
|
||
" output['d']['O'] = dnet(output['d']['O'],dm16,zhm16,xmm,16)\n",
|
||
" rl = pdm[1,1]*pdl[1,16]*(1+flags['sw'][0]*pdl[0,23]*(inputp['f107A']-150))\n",
|
||
" hc16 = pdm[1,5]*pdl[1,3]\n",
|
||
" zc16 = pdm[1,4]*pdl[1,2]\n",
|
||
" hc216 = pdm[1,5]*pdl[1,4]\n",
|
||
" output['d']['O'] = output['d']['O']*ccor2(z,rl,hc16,zc16,hc216)\n",
|
||
" # Chemistry correction \n",
|
||
" hcc16 = pdm[1,7]*pdl[1,13]\n",
|
||
" zcc16 = pdm[1,6]*pdl[1,12]\n",
|
||
" rc16 = pdm[1,3]*pdl[1,14]\n",
|
||
" # Net density corrected at Alt\n",
|
||
" output['d']['O'] = output['d']['O']*ccor(z,rc16,hcc16,zcc16)\n",
|
||
" \n",
|
||
" # O2 DENSITY\n",
|
||
" # Density variation factor at Zlb \n",
|
||
" tinf_tmp,varli = globe7(pd[4],inputp,flags)\n",
|
||
" g32 = flags['sw'][20]*tinf_tmp\n",
|
||
" # Diffusive density at Zlb \n",
|
||
" db32 = pdm[3,0]*np.exp(g32)*pd[4,0]\n",
|
||
" # Diffusive density at Alt \n",
|
||
" output['d']['O2'],output['t']['TG'] = densu(z,db32,tinf,tlb, 32,alpha[3],output['t']['TG'],ptm[5],s, zn1,meso_tn1,meso_tgn1,gsurf,re)\n",
|
||
" dd = output['d']['O2'];\n",
|
||
" if flags['sw'][14]:\n",
|
||
" if z <= altl[3]:\n",
|
||
" # Turbopause \n",
|
||
" zh32 = pdm[3,2]\n",
|
||
" # Mixed density at Zlb\n",
|
||
" b32,output['t']['TG'] = densu(zh32,db32,tinf,tlb,32-xmm,alpha[3]-1, output['t']['TG'],ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re)\n",
|
||
" # Mixed density at Alt \n",
|
||
" dm32,output['t']['TG'] = densu(z,b32,tinf,tlb,xmm,0,output['t']['TG'],ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re)\n",
|
||
" zhm32 = zhm28\n",
|
||
" # Net density at Alt\n",
|
||
" output['d']['O2'] = dnet(output['d']['O2'],dm32,zhm32,xmm,32)\n",
|
||
" # Correction to specified mixing ratio at ground \n",
|
||
" rl = np.log(b28*pdm[3,1]/b32)\n",
|
||
" hc32 = pdm[3,5]*pdl[1,7]\n",
|
||
" zc32 = pdm[3,4]*pdl[1,6]\n",
|
||
" output['d']['O2'] = output['d']['O2']*ccor(z,rl,hc32,zc32)\n",
|
||
"\n",
|
||
" # Correction for general departure from diffusive equilibrium above Zlb */\n",
|
||
" hcc32 = pdm[3,7]*pdl[1,22]\n",
|
||
" hcc232 = pdm[3,7]*pdl[0,22]\n",
|
||
" zcc32 = pdm[3,6]*pdl[1,21]\n",
|
||
" rc32 = pdm[3,3]*pdl[1,23]*(1+flags['sw'][0]*pdl[0,23]*(inputp['f107A']-150))\n",
|
||
" # Net density corrected at Alt \n",
|
||
" output['d']['O2'] = output['d']['O2']*ccor2(z,rc32,hcc32,zcc32,hcc232)\n",
|
||
" # AR DENSITY\n",
|
||
" # Density variation factor at Zlb \n",
|
||
" tinf_tmp,varli = globe7(pd[5],inputp,flags)\n",
|
||
" g40 = flags['sw'][20]*tinf_tmp\n",
|
||
" # Diffusive density at Zlb \n",
|
||
" db40 = pdm[4,0]*np.exp(g40)*pd[5,0]\n",
|
||
" # Diffusive density at Alt\n",
|
||
" output['d']['AR'],output['t']['TG'] = densu(z,db40,tinf,tlb, 40,alpha[4],output['t']['TG'],ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re)\n",
|
||
" dd = output['d']['AR']\n",
|
||
" if flags['sw'][14] and z <= altl[4]:\n",
|
||
" # Turbopause\n",
|
||
" zh40 = pdm[4,2]\n",
|
||
" # Mixed density at Zlb \n",
|
||
" b40,output['t']['TG'] = densu(zh40,db40,tinf,tlb,40-xmm,alpha[4]-1,output['t']['TG'],ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re)\n",
|
||
" # Mixed density at Alt\n",
|
||
" dm40,output['t']['TG'] = densu(z,b40,tinf,tlb,xmm,0,output['t']['TG'],ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re)\n",
|
||
" zhm40 = zhm28\n",
|
||
" # Net density at Alt \n",
|
||
" output['d']['AR'] = dnet(output['d']['AR'],dm40,zhm40,xmm,40)\n",
|
||
" # Correction to specified mixing ratio at ground \n",
|
||
" rl = np.log(b28*pdm[4,1]/b40)\n",
|
||
" hc40 = pdm[4,5]*pdl[1,9]\n",
|
||
" zc40 = pdm[4,4]*pdl[1,8]\n",
|
||
" # Net density corrected at Alt\n",
|
||
" output['d']['AR'] = output['d']['AR']*ccor(z,rl,hc40,zc40)\n",
|
||
" \n",
|
||
" # HYDROGEN DENSITY \n",
|
||
" # Density variation factor at Zlb */\n",
|
||
" tinf_tmp,varli = globe7(pd[6], inputp, flags)\n",
|
||
" g1 = flags['sw'][20]*tinf_tmp\n",
|
||
" # Diffusive density at Zlb \n",
|
||
" db01 = pdm[5,0]*np.exp(g1)*pd[6,0]\n",
|
||
" # Diffusive density at Alt\n",
|
||
" output['d']['H'],output['t']['TG']=densu(z,db01,tinf,tlb,1,alpha[6],output['t']['TG'],ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re)\n",
|
||
" dd = output['d']['H']\n",
|
||
" if flags['sw'][14] and z <= altl[6]:\n",
|
||
" # Turbopause \n",
|
||
" zh01 = pdm[5,2]\n",
|
||
" # Mixed density at Zlb\n",
|
||
" b01,output['t']['TG'] = densu(zh01,db01,tinf,tlb,1-xmm,alpha[6]-1, output['t']['TG'],ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re)\n",
|
||
" # Mixed density at Alt \n",
|
||
" dm01,output['t']['TG'] = densu(z,b01,tinf,tlb,xmm,0,output['t']['TG'],ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re)\n",
|
||
" zhm01 = zhm28\n",
|
||
" # Net density at Alt\n",
|
||
" output['d']['H'] = dnet(output['d']['H'],dm01,zhm01,xmm,1)\n",
|
||
" # Correction to specified mixing ratio at ground \n",
|
||
" rl = np.log(b28*pdm[5,1]*np.abs(pdl[1,17])/b01)\n",
|
||
" hc01 = pdm[5,5]*pdl[1,11]\n",
|
||
" zc01 = pdm[5,4]*pdl[1,10]\n",
|
||
" output['d']['H'] = output['d']['H']*ccor(z,rl,hc01,zc01)\n",
|
||
" # Chemistry correction \n",
|
||
" hcc01 = pdm[5,7]*pdl[1,19]\n",
|
||
" zcc01 = pdm[5,6]*pdl[1,18]\n",
|
||
" rc01 = pdm[5,3]*pdl[1,20]\n",
|
||
" # Net density corrected at Alt\n",
|
||
" output['d']['H'] = output['d']['H']*ccor(z,rc01,hcc01,zcc01)\n",
|
||
" \n",
|
||
" # ATOMIC NITROGEN DENSITY \n",
|
||
" # Density variation factor at Zlb */\n",
|
||
" tinf_tmp,varli = globe7(pd[7],inputp,flags)\n",
|
||
" g14 = flags['sw'][20]*tinf_tmp\n",
|
||
" # Diffusive density at Zlb \n",
|
||
" db14 = pdm[6,0]*np.exp(g14)*pd[7,0]\n",
|
||
" # Diffusive density at Alt \n",
|
||
" output['d']['N'],output['t']['TG']=densu(z,db14,tinf,tlb,14,alpha[7],output['t']['TG'],ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re)\n",
|
||
" dd = output['d']['N']\n",
|
||
" if flags['sw'][14] and z <= altl[7]:\n",
|
||
" # Turbopause\n",
|
||
" zh14 = pdm[6,2]\n",
|
||
" # Mixed density at Zlb\n",
|
||
" b14,output['t']['TG'] = densu(zh14,db14,tinf,tlb,14-xmm,alpha[7]-1, output['t']['TG'],ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re)\n",
|
||
" # Mixed density at Alt \n",
|
||
" dm14,output['t']['TG'] = densu(z,b14,tinf,tlb,xmm,0,output['t']['TG'],ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re)\n",
|
||
" zhm14 = zhm28\n",
|
||
" # Net density at Alt\n",
|
||
" output['d']['N'] = dnet(output['d']['N'],dm14,zhm14,xmm,14)\n",
|
||
" # Correction to specified mixing ratio at ground \n",
|
||
" rl = np.log(b28*pdm[6,1]*np.abs(pdl[0,2])/b14)\n",
|
||
" hc14 = pdm[6,5]*pdl[0,1]\n",
|
||
" zc14 = pdm[6,4]*pdl[0,0]\n",
|
||
" output['d']['N'] = output['d']['N']*ccor(z,rl,hc14,zc14)\n",
|
||
" # Chemistry correction\n",
|
||
" hcc14 = pdm[6,7]*pdl[0,4]\n",
|
||
" zcc14 = pdm[6,6]*pdl[0,3]\n",
|
||
" rc14 = pdm[6,3]*pdl[0,5]\n",
|
||
" # Net density corrected at Alt\n",
|
||
" output['d']['N'] = output['d']['N']*ccor(z,rc14,hcc14,zcc14)\n",
|
||
" \n",
|
||
" # Anomalous OXYGEN DENSITY \n",
|
||
" tinf_tmp,varli = globe7(pd[8],inputp,flags)\n",
|
||
" g16h = flags['sw'][20]*tinf_tmp\n",
|
||
" db16h = pdm[7,0]*np.exp(g16h)*pd[8,0]\n",
|
||
" tho = pdm[7,9]*pdl[0,6]\n",
|
||
" dd,output['t']['TG'] = densu(z,db16h,tho,tho,16,alpha[8],output['t']['TG'],ptm[5],s, zn1,meso_tn1,meso_tgn1,gsurf,re)\n",
|
||
" zsht = pdm[7,5]\n",
|
||
" zmho = pdm[7,4]\n",
|
||
" zsho = scalh(zmho,16,tho,gsurf,re)\n",
|
||
" output['d']['ANM O'] = dd*np.exp(-zsht/zsho*(np.exp(-(z-zmho)/zsht)-1))\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",
|
||
"\n",
|
||
" # temperature \n",
|
||
" z = inputp['alt']\n",
|
||
" ddum,output['t']['TG'] = densu(z,1, tinf, tlb, 0, 0, output['t']['TG'], ptm[5], s, zn1, meso_tn1, meso_tgn1,gsurf,re)\n",
|
||
"\n",
|
||
" # convert to g/cm^3 \n",
|
||
" for key in output['d'].keys():\n",
|
||
" output['d'][key] = output['d'][key]*1.0E6\n",
|
||
" output['d']['RHO'] = output['d']['RHO']/1000 \n",
|
||
" return output,dm28,[meso_tn1,meso_tn2,meso_tn3,meso_tgn1,meso_tgn2,meso_tgn3],varli\n",
|
||
"\n",
|
||
"# ------------------------------------------------------------------- #\n",
|
||
"# ----------------- DOWNLOAD AND UPDATE SW DATA -------------------- #\n",
|
||
"# ------------------------------------------------------------------- #\n",
|
||
"\n",
|
||
"def update_sw(direc=None):\n",
|
||
" \n",
|
||
" if direc is None:\n",
|
||
" home = getenv('HOME')\n",
|
||
" direc = home + '/src/sw-data/'\n",
|
||
" \n",
|
||
" swfile = direc + 'SW-All.txt'\n",
|
||
" url = 'https://www.celestrak.com/SpaceData/SW-All.txt'\n",
|
||
" \n",
|
||
" if not path.exists(swfile):\n",
|
||
" print('Downloading the latest space weather data')\n",
|
||
" urlretrieve(url, swfile)\n",
|
||
" else:\n",
|
||
" modified_time = datetime.fromtimestamp(path.getmtime(swfile))\n",
|
||
" if datetime.now() > 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<ipython-input-3-89022222782d>\u001b[0m in \u001b[0;36m<module>\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
|
||
}
|