{ "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 }