From 0e0097c8968e53db3775397971ab8456bc9b2522 Mon Sep 17 00:00:00 2001 From: Chunxiao Li Date: Mon, 7 Jun 2021 12:59:24 +0800 Subject: [PATCH] A part separated from the source file 'nrlmsise00' --- pyatmos/msise/nrlmsise00.py | 1084 +-------------------------- pyatmos/msise/nrlmsise00_subfunc.py | 879 ++++++++++++++++++++++ 2 files changed, 914 insertions(+), 1049 deletions(-) create mode 100644 pyatmos/msise/nrlmsise00_subfunc.py diff --git a/pyatmos/msise/nrlmsise00.py b/pyatmos/msise/nrlmsise00.py index 6c9b6ef..00665ff 100644 --- a/pyatmos/msise/nrlmsise00.py +++ b/pyatmos/msise/nrlmsise00.py @@ -1,1066 +1,48 @@ -# -------------------------------------------------------------------- # -# ------------------- nrlmisise-00 model 2001 -------------------- # -# -------------------------------------------------------------------- # - -''' -pyatmos nrlmsise00 - -This submodule defines the following functions: - -nrlmsis00_data - Read the data block from nrlmsis00_data.npz - -tselec - - -glatf - - -ccor - - -ccor2 - - -scalh - - -dnet - - -zeta - - -densm - - -densu - - -g0 - -sumex - -sg0 - -lengendre - -globe7 - -glob7s - -gtd7 - -dtd7d - -gts7 - -nrlmsise00 -''' - import numpy as np -from scipy.interpolate import CubicSpline from astropy.time import Time -from pyshtools.legendre import PLegendreA,PlmIndex -import pkg_resources +from .nrlmsise00_subfunc import gtd7,gtd7d from .spaceweather import get_sw from ..utils.utils import wraplon,hms_conver +from ..class_atmos import ATMOS -# ======================== read data block ========================== # - -def nrlmsis00_data(): - ''' - Read the data block from nrlmsis00_data.npz - - Usage: pt,pd,ps,pdl,ptm,pdm,ptl,pma,sam,pavgm = nrlmsis00_data() - - Inputs: - None - - Outputs: - pt: [float array] TEMPERATURE - pd: [2d float array] DENSITY of HE, O, N2, Total mass, O2, AR, H, N, HOT O - ps: [float array] S PARAM - pdl: [2d float array] TURBO - ptm: [float array] - pdm: [2d float array] - ptl: [2d float array] - pma: [2d float array] - sam: [float array] SEMIANNUAL MULT SAM - pavgm: [float array] MIDDLE ATMOSPHERE AVERAGES - ''' - data_path = pkg_resources.resource_filename('pyatmos', 'data/') - data = np.load(data_path+'nrlmsis00_data.npz') - pt,pd,ps,pdl = data['pt'],data['pd'],data['ps'],data['pdl'] - ptm,pdm,ptl,pma = data['ptm'],data['pdm'],data['ptl'],data['pma'] - sam,pavgm = data['sam'],data['pavgm'] - return pt,pd,ps,pdl,ptm,pdm,ptl,pma,sam,pavgm - -# ============================= tselec ============================== # - -def tselec(switches): - # len(switches) is equal to 23 - flags = {'sw':np.zeros(23),'swc':np.zeros(23)} - for i in range(23): - if i != 8: - if switches[i] == 1: - flags['sw'][i] = 1 - else: - flags['sw'][i] = 0 - if switches[i] > 0: - flags['swc'][i] = 1 - else: - flags['swc'][i] = 0 - else: - flags['sw'][i] = switches[i] - flags['swc'][i] = switches[i] - return flags - -# ============================= glatf =============================== # - -def glatf(lat): - c2 = np.cos(2*np.deg2rad(lat)) - gv = 980.616*(1 - 0.0026373*c2) - reff = 2*gv/(3.085462E-6 + 2.27E-9*c2)*1E-5 - return gv,reff - -# ============================= ccor ================================ # - -def ccor(alt,r,h1,zh): - ''' - Chemistry/dissociation correction for msis models - - alt - altitude - r - target ratio - h1 - transition scale length - zh - altitude of 1/2 R - ''' - e = (alt - zh)/h1 - if e > 70: - return 1 - elif e < -70: - return np.exp(r) - else: - return np.exp(r/(1 + np.exp(e))) - -# ============================== ccor2 ============================== # - -def ccor2(alt,r,h1,zh,h2): - ''' - Chemistry/dissociation correction for msis models - - alt - altitude - r - target ratio - f1 - transition scale length 1 - zh - altitude of 1/2 R - h2 - transition scale length 2 - ''' - e1 = (alt - zh)/h1 - e2 = (alt - zh)/h2 - if e1 > 70 or e2 > 70: - return 1 - if e1 < -70 and e2 < -70: - return np.exp(r) - ex1,ex2 = np.exp([e1,e2]) - ccor2v = r/(1 + 0.5*(ex1 + ex2)) - return np.exp(ccor2v) - -# =============================== scalh ============================= # - -def scalh(alt,xm,temp,gsurf,re): - rgas = 831.4 - g = rgas*temp/(gsurf/(1 + alt/re)**2*xm) - return g - -# ================================ dnet ============================= # - -def dnet(dd,dm,zhm,xmm,xm): - ''' - Turbopause correction for msis models - Root mean density - dd - diffusive density - dm - full mixed density - zhm - transition scale length - xmm - full mixed molecular weight - xm - species molecular weight - dnet - combined density - ''' - a = zhm/(xmm - xm) - if not (dm > 0 and dd > 0): - print('dnet log error {0:.1f} {1:.1f} {2:.1f}'.format(dm,dd,xm)) - if dd == 0 and dm == 0: dd = 1 - if dm == 0: return dd - if dd == 0: return dm - ylog = a*np.log(dm/dd) - if ylog < -10: return dd - if ylog > 10: return dm - a = dd*(1 + np.exp(ylog))**(1/a) - return a - -# ================================ zeta ============================= # - -def zeta(zz,zl,re): - return (zz - zl)*(re + zl)/(re + zz) - -# =============================== densm ============================= # - -def densm(alt, d0, xm, tz, zn3, tn3, tgn3, zn2, tn2, tgn2,gsurf,re): - # Calculate Temperature and Density Profiles for lower atmos. - # call zeta - rgas = 831.4 - densm_tmp = d0 - tz_tmp = tz - mn3,mn2 = len(zn3),len(zn2) - - if alt > zn2[0]: - if xm == 0: - densm_tmp = tz - return densm_tmp,tz_tmp - else: - densm_tmp = d0 - return densm_tmp,tz_tmp - - # stratosphere/mesosphere temperature - if alt > zn2[mn2-1]: - z = alt - else: - z = zn2[mn2-1] - mn = mn2 - xs,ys = [np.zeros(mn) for i in range(2)] - z1,z2 = zn2[0],zn2[mn-1] - t1,t2=tn2[0],tn2[mn-1] - zg,zgdif = zeta(z,z1,re),zeta(z2,z1,re) - - # set up spline nodes - for k in range(mn): - xs[k] = zeta(zn2[k],z1,re)/zgdif - ys[k] = 1/tn2[k] - yd1 = -tgn2[0]/t1**2*zgdif - yd2 = -tgn2[1]/t2**2*zgdif*((re + z2)/(re + z1))**2 - - # calculate spline coefficients - cs = CubicSpline(xs,ys,bc_type=((1,yd1),(1,yd2))) - x = zg/zgdif - y = cs(x) - - # temperature at altitude - tz_tmp = 1/y - if xm != 0: - # calaculate stratosphere / mesospehere density - glb = gsurf/(1 + z1/re)**2 - gamm = xm*glb*zgdif/rgas - - # Integrate temperature profile - yi = cs.integrate(xs[0],x) - expl = gamm*yi - if expl > 50: - expl = 50 - # Density at altitude - densm_tmp = densm_tmp*(t1/tz_tmp)*np.exp(-expl) - if alt > zn3[0]: - if xm == 0: - densm_tmp = tz_tmp - return densm_tmp,tz_tmp - else: - return densm_tmp,tz_tmp - - # troposhere / stratosphere temperature - z = alt - mn = mn3 - xs,ys = [np.zeros(mn) for i in range(2)] - z1,z2 = zn3[0],zn3[mn-1] - t1,t2 = tn3[0],tn3[mn-1] - zg,zgdif = zeta(z,z1,re),zeta(z2,z1,re) - - # set up spline nodes - for k in range(mn): - xs[k] = zeta(zn3[k],z1,re)/zgdif - ys[k] = 1/tn3[k] - yd1 = -tgn3[0]/t1**2*zgdif - yd2 = -tgn3[1]/t2**2*zgdif*((re+z2)/(re+z1))**2 - - # calculate spline coefficients - cs = CubicSpline(xs,ys,bc_type=((1,yd1),(1,yd2))) - x = zg/zgdif - y = cs(x) - - # temperature at altitude - tz_tmp = 1/y - if xm != 0: - # calaculate tropospheric / stratosphere density - glb = gsurf/(1 + z1/re)**2 - gamm = xm*glb*zgdif/rgas - - # Integrate temperature profile - yi = cs.integrate(xs[0],x) - expl = gamm*yi - if expl > 50: expl = 50 - # Density at altitude - densm_tmp = densm_tmp*(t1/tz_tmp)*np.exp(-expl) - - if xm == 0: - densm_tmp = tz_tmp - return densm_tmp,tz_tmp - else: - return densm_tmp,tz_tmp - -# =============================== densu ============================= # - -def densu (alt,dlb,tinf,tlb,xm,alpha,tz,zlb,s2,zn1,tn1,tgn1,gsurf,re): - # Calculate Temperature and Density Profiles for MSIS models - # New lower thermo polynomial - # call: zeta, - - rgas = 831.4 - densu_tmp = 1 - mn1 = len(zn1) - # joining altitudes of Bates and spline - za = zn1[0] - if alt > za: - z = alt - else: - z = za - # geopotential altitude difference from ZLB - zg2 = zeta(z,zlb,re) - - # Bates temperature - tt = tinf - (tinf - tlb)*np.exp(-s2*zg2) - ta = tz = tt - densu_tmp = tz_tmp = tz - - if alt < za: - # calculate temperature below ZA - # temperature gradient at ZA from Bates profile - dta = (tinf - ta)*s2*((re + zlb)/(re + za))**2 - tgn1[0],tn1[0] = dta,ta - if alt > zn1[mn1-1]: - z = alt - else: - z = zn1[mn1-1] - mn = mn1 - xs,ys = [np.zeros(mn) for i in range(2)] - z1,z2 = zn1[0],zn1[mn-1] - t1,t2 = tn1[0],tn1[mn-1] - # geopotental difference from z1 - zg,zgdif = zeta(z,z1,re),zeta(z2,z1,re) - # set up spline nodes - for k in range(mn): - xs[k] = zeta(zn1[k],z1,re)/zgdif - ys[k] = 1/tn1[k] - # end node derivatives - yd1 = -tgn1[0]/t1**2*zgdif - yd2 = -tgn1[1]/t2**2*zgdif*((re + z2)/(re + z1))**2 - # calculate spline coefficients - cs = CubicSpline(xs,ys,bc_type=((1,yd1),(1,yd2))) - x = zg/zgdif - y = cs(x) - # temperature at altitude - tz_tmp = 1/y - densu_tmp = tz_tmp - if xm == 0: return densu_tmp,tz_tmp - - # calculate density above za - glb = gsurf/(1 + zlb/re)**2 - gamma = xm*glb/(s2*rgas*tinf) - expl = np.exp(-s2*gamma*zg2) - if expl > 50: expl = 50 - if tt <= 0: expl = 50 - - # density at altitude - densa = dlb*(tlb/tt)**(1 + alpha + gamma)*expl - densu_tmp = densa - if alt >= za: return densu_tmp,tz_tmp - - # calculate density below za - glb = gsurf/(1 + z1/re)**2 - gamm = xm*glb*zgdif/rgas - - # integrate spline temperatures - yi = cs.integrate(xs[0],x) - expl = gamm*yi - if expl > 50: expl = 50 - if tz_tmp <= 0: expl = 50 - - # density at altitude - densu_tmp = densu_tmp*(t1/tz_tmp)**(1 + alpha)*np.exp(-expl) - return densu_tmp,tz_tmp - -# =============== 3hr magnetic activity functions =================== # - -# Eq. A24d -def g0(a,p): - return (a - 4 + (p[25] - 1)*(a - 4 + (np.exp(-np.abs(p[24])*(a - 4)) - 1) / np.abs(p[24]))) - -# Eq. A24c -def sumex(ex): - return (1 + (1 - ex**19)/(1 - ex)*ex**0.5) - -# Eq. A24a -def sg0(ex,p,ap): - # call sumex, g0 - return (g0(ap[1],p) + g0(ap[2],p)*ex + g0(ap[3],p)*ex**2 + \ - g0(ap[4],p)*ex**3 + (g0(ap[5],p)*ex**4 + \ - g0(ap[6],p)*ex**12)*(1-ex**8)/(1-ex))/sumex(ex) - -# ================== Associated Legendre polynomials ================ # - -def lengendre(g_lat,lmax = 8): - # Index of PLegendreA_x can be calculated by PlmIndex(l,m) - x = np.sin(np.deg2rad(g_lat)) - PLegendreA_x = PLegendreA(lmax,x) - return PLegendreA_x - -# =============================== globe7 ============================ # - -def globe7(p,inputp,flags): - # calculate G(L) function - # Upper Thermosphere Parameters - # call: lengendre,sg0 - t = np.zeros(15) - sr = 7.2722E-5 - dr = 1.72142E-2 - hr = 0.2618 - - apdf = 0 - apt = np.zeros(4) - tloc = inputp['lst'] - - if not (flags['sw'][6]==0 and flags['sw'][7]==0 and flags['sw'][13]==0): - stloc,ctloc = np.sin(hr*tloc),np.cos(hr*tloc) - s2tloc,c2tloc = np.sin(2*hr*tloc),np.cos(2*hr*tloc) - s3tloc,c3tloc = np.sin(3*hr*tloc),np.cos(3*hr*tloc) - cd32 = np.cos(dr*(inputp['doy'] - p[31])) - cd18 = np.cos(2*dr*(inputp['doy'] - p[17])) - cd14 = np.cos(dr*(inputp['doy'] - p[13])) - cd39 = np.cos(2*dr*(inputp['doy'] - p[38])) - - # F10.7 effect - df = inputp['f107'] - inputp['f107A'] - dfa = inputp['f107A'] - 150 - t[0] = p[19]*df*(1 + p[59]*dfa) + p[20]*df**2 + p[21]*dfa + p[29]*dfa**2 - f1 = 1 + (p[47]*dfa + p[19]*df + p[20]*df**2)*flags['swc'][0] - f2 = 1 + (p[49]*dfa + p[19]*df + p[20]*df**2)*flags['swc'][0] - - plg = lengendre(inputp['g_lat']) - - # time independent - 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] - - # symmetrical annual - t[2] = p[18]*cd32 - - # symmetrical semiannual - t[3] = (p[15] + p[16]*plg[3])*cd18 - - # asymmetrical annual - t[4] = f1*(p[9]*plg[1] + p[10]*plg[6])*cd14 - - # asymmetrical semiannual - t[5] = p[37]*plg[1]*cd39 - - # diurnal - if flags['sw'][6]: - t71 = p[11]*plg[4]*cd14*flags['swc'][4] - t72 = p[12]*plg[4]*cd14*flags['swc'][4] - 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) - - # semiannual - if flags['sw'][7]: - t81 = (p[23]*plg[8] + p[35]*plg[17])*cd14*flags['swc'][4] - t82 = (p[33]*plg[8] + p[36]*plg[17])*cd14*flags['swc'][4] - t[7] = f2*((p[5]*plg[5] + p[41]*plg[12] + t81)*c2tloc +(p[8]*plg[5] + p[42]*plg[12] + t82)*s2tloc) - - # terdiurnal - if flags['sw'][13]: - 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) - - # magnetic activity based on daily ap - if flags['sw'][8] == -1: - ap = inputp['ap_a'] - if p[51]!= 0: - exp1 = np.exp(-10800*np.abs(p[51])/(1 + p[138]*(45 - np.abs(inputp['g_lat'])))) - if exp1 > 0.99999: exp1 = 0.99999 - if p[24] < 1E-4: p[24] = 1E-4 - apt[0] = sg0(exp1,p,ap) - # apt[1] = sg2(exp1,p,ap) - # apt[2] = sg0(exp2,p,ap) - # apt[3] = sg2(exp2,p,ap) - - if flags['sw'][8]: - t[8] = apt[0]*(p[50] + p[96]*plg[3] + p[54]*plg[10] + \ - (p[125]*plg[1] + p[126]*plg[6] + p[127]*plg[15])*cd14*flags['swc'][4] + \ - (p[128]*plg[2] + p[129]*plg[7] + p[130]*plg[16])*flags['swc'][6]*np.cos(hr*(tloc - p[131]))) - else: - apd = inputp['ap'] - 4 - p44 = p[43] - p45 = p[44] - if p44 < 0: p44 = 1E-5 - apdf = apd + (p45 - 1)*(apd + (np.exp(-p44*apd) - 1)/p44) - if flags['sw'][8]: - t[8]=apdf*(p[32] + p[45]*plg[3] + p[34]*plg[10] + \ - (p[100]*plg[1] + p[101]*plg[6] + p[102]*plg[15])*cd14*flags['swc'][4] + - (p[121]*plg[2] + p[122]*plg[7] + p[123]*plg[16])*flags['swc'][6]*np.cos(hr*(tloc - p[124]))) - - if flags['sw'][9] and inputp['g_lon'] > -1000: - # longitudinal - if flags['sw'][10]: - t[10] = (1 + p[80]*dfa*flags['swc'][0])*((p[64]*plg[4] + p[65]*plg[11] + p[66]*plg[22]\ - + p[103]*plg[2] + p[104]*plg[7] + p[105]*plg[16]\ - + flags['swc'][4]*(p[109]*plg[2] + p[110]*plg[7] + p[111]*plg[16])*cd14)*np.cos(np.deg2rad(inputp['g_lon'])) \ - +(p[90]*plg[4]+p[91]*plg[11]+p[92]*plg[22] + p[106]*plg[2]+p[107]*plg[7]+p[108]*plg[16]\ - + flags['swc'][4]*(p[112]*plg[2] + p[113]*plg[7] + p[114]*plg[16])*cd14)*np.sin(np.deg2rad(inputp['g_lon']))) - - # ut and mixed ut, longitude - if flags['sw'][11]: - t[11]=(1 + p[95]*plg[1])*(1 + p[81]*dfa*flags['swc'][0])*\ - (1 + p[119]*plg[1]*flags['swc'][4]*cd14)*\ - ((p[68]*plg[1] + p[69]*plg[6] + p[70]*plg[15])*np.cos(sr*(inputp['sec'] - p[71]))) - t[11] += flags['swc'][10]*(p[76]*plg[8] + p[77]*plg[17] + p[78]*plg[30])*\ - np.cos(sr*(inputp['sec'] - p[79]) + 2*np.deg2rad(inputp['g_lon']))*(1 + p[137]*dfa*flags['swc'][0]) - - # ut, longitude magnetic activity - if flags['sw'][10]: - if flags['sw'][8] == -1: - if p[51]: - t[12] = apt[0]*flags['swc'][10]*(1 + p[132]*plg[1])*\ - ((p[52]*plg[4] + p[98]*plg[11] + p[67]*plg[22])* np.cos(np.deg2rad(inputp['g_lon'] - p[97])))\ - + apt[0]*flags['swc'][10]*flags['swc'][4]*(p[133]*plg[2] + p[134]*plg[7] + p[135]*plg[16])*\ - cd14*np.cos(np.deg2rad(inputp['g_lon'] - p[136])) + apt[0]*flags['swc'][11]* \ - (p[55]*plg[1] + p[56]*plg[6] + p[57]*plg[15])*np.cos(sr*(inputp['sec'] - p[58])) - else: - t[12] = apdf*flags['swc'][10]*(1 + p[120]*plg[1])*((p[60]*plg[4] + p[61]*plg[11] + p[62]*plg[22])*\ - np.cos(np.deg2rad(inputp['g_lon']-p[63])))+apdf*flags['swc'][10]*flags['swc'][4]* \ - (p[115]*plg[2] + p[116]*plg[7] + p[117]*plg[16])* \ - cd14*np.cos(np.deg2rad(inputp['g_lon'] - p[118])) \ - + apdf*flags['swc'][11]*(p[83]*plg[1] + p[84]*plg[6] + p[85]*plg[15])* np.cos(sr*(inputp['sec'] - p[75])) - - # parms not used: 82, 89, 99, 139-149 - tinf = p[30] - for i in range(14): - tinf = tinf + np.abs(flags['sw'][i])*t[i] - return tinf,[dfa,plg,ctloc,stloc,c2tloc,s2tloc,s3tloc,c3tloc,apdf,apt] - -# =============================== glob7s ============================ # - -def glob7s(p,inputp,flags,varli): - # version of globe for lower atmosphere 10/26/99 - # call: lengendre,sg0 - pset = 2 - t = np.zeros(14) - dr = 1.72142E-2 - [dfa,plg,ctloc,stloc,c2tloc,s2tloc,s3tloc,c3tloc,apdf,apt] = varli - - # confirm parameter set - if p[99] == 0: p[99] = pset - if p[99] != pset: - print("Wrong parameter set for glob7s") - return -1 - - for j in range(14): - t[j] = 0 - cd32 = np.cos(dr*(inputp['doy'] - p[31])) - cd18 = np.cos(2*dr*(inputp['doy'] - p[17])) - cd14 = np.cos(dr*(inputp['doy'] - p[13])) - cd39 = np.cos(2*dr*(inputp['doy'] - p[38])) - - # F10.7 - t[0] = p[21]*dfa - - # time independent - 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] - - # symmetrical annual - t[2] = (p[18] + p[47]*plg[3] + p[29]*plg[10])*cd32 - - # symmetrical semiannual - t[3] = (p[15] + p[16]*plg[3] + p[30]*plg[10])*cd18 - - # asymmetrical annual - t[4] = (p[9]*plg[1] + p[10]*plg[6] + p[20]*plg[15])*cd14 - - # asymmetrical semiannual - t[5] = p[37]*plg[1]*cd39; - - # diurnal - if flags['sw'][6]: - t71 = p[11]*plg[4]*cd14*flags['swc'][4] - t72 = p[12]*plg[4]*cd14*flags['swc'][4] - t[6] = ((p[3]*plg[2] + p[4]*plg[7] + t71)*ctloc + (p[6]*plg[2] + p[7]*plg[7] + t72)*stloc) - - # semidiurnal - if flags['sw'][7]: - t81 = (p[23]*plg[8] + p[35]*plg[17])*cd14*flags['swc'][4] - t82 = (p[33]*plg[8] + p[36]*plg[17])*cd14*flags['swc'][4] - t[7] = ((p[5]*plg[5] + p[41]*plg[12] + t81)*c2tloc + (p[8]*plg[5] + p[42]*plg[12] + t82)*s2tloc) - - # terdiurnal - if flags['sw'][13]: - t[13] = p[39]*plg[9]*s3tloc + p[40]*plg[9]*c3tloc - - # magnetic activity - if flags['sw'][8]: - if flags['sw'][8]==1: - t[8] = apdf * (p[32] + p[45]*plg[3]*flags['swc'][1]) - if flags['sw'][8]==-1: - t[8]=(p[50]*apt[0] + p[96]*plg[3]*apt[0]*flags['swc'][1]) - - # longitudinal - if not (flags['sw'][9]==0 or flags['sw'][10]==0 or inputp['g_lon']<=-1000): - t[10] = (1 + plg[1]*(p[80]*flags['swc'][4]*np.cos(dr*(inputp['doy'] - p[81]))\ - + p[85]*flags['swc'][5]*np.cos(2*dr*(inputp['doy'] - p[86])))\ - + p[83]*flags['swc'][2]*np.cos(dr*(inputp['doy'] - p[84]))\ - + p[87]*flags['swc'][3]*np.cos(2*dr*(inputp['doy'] - p[88])))\ - *((p[64]*plg[4] + p[65]*plg[11] + p[66]*plg[22]\ - + p[74]*plg[2] + p[75]*plg[7] + p[76]*plg[16])*np.cos(np.deg2rad(inputp['g_lon']))\ - + (p[90]*plg[4] + p[91]*plg[11] + p[92]*plg[22]\ - + p[77]*plg[2] + p[78]*plg[7] + p[79]*plg[16])*np.sin(np.deg2rad(inputp['g_lon']))) - - tt = 0 - for i in range(14): - tt += np.abs(flags['sw'][i])*t[i] - return tt - -# =============================== gtd7 ============================== # - -def gtd7(inputp,switches): - tz = 0 - zn3 = np.array([32.5,20.0,15.0,10.0,0.0]) - zn2 = np.array([72.5,55.0,45.0,32.5]) - zmix= 62.5 - - output = {'d':{'He':0,'O':0,'N2':0,'O2':0,'AR':0,'RHO':0,'H':0,'N':0,'ANM O':0},\ - 't':{'TINF':0,'TG':0}} - - flags = tselec(switches) - - # Latitude variation of gravity (none for sw[1]=0) - xlat = inputp['g_lat'] - if flags['sw'][1]==0: xlat = 45 - gsurf,re = glatf(xlat) - pt,pd,ps,pdl,ptm,pdm,ptl,pma,sam,pavgm = nrlmsis00_data() - xmm = pdm[2,4] - - # thermosphere/mesosphere (above zn2[0]) - if inputp['alt'] > zn2[0]: - altt = inputp['alt'] - else: - altt = zn2[0] - - tmp = inputp['alt'] - inputp['alt'] = altt - 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) - altt = inputp['alt'] - inputp['alt'] = tmp - # metric adjustment - dm28m = dm28*1E6 - output['t']['TINF'] = soutput['t']['TINF'] - output['t']['TG'] = soutput['t']['TG'] - if inputp['alt'] >= zn2[0]: - output['d'] = soutput['d'] - return output - # lower mesosphere/upper stratosphere (between zn3[0] and zn2[0]) - # Temperature at nodes and gradients at end nodes - # Inverse temperature a linear function of spherical harmonics - - varli = [dfa,plg,ctloc,stloc,c2tloc,s2tloc,s3tloc,c3tloc,apdf,apt] - - meso_tgn2[0] = meso_tgn1[1] - meso_tn2[0] = meso_tn1[4] - meso_tn2[1] = pma[0,0]*pavgm[0]/(1-flags['sw'][19]*glob7s(pma[0], inputp, flags,varli)) - meso_tn2[2] = pma[1,0]*pavgm[1]/(1-flags['sw'][19]*glob7s(pma[1], inputp, flags,varli)) - meso_tn2[3] = pma[2,0]*pavgm[2]/(1-flags['sw'][19]*flags['sw'][21]*glob7s(pma[2], inputp, flags,varli)) - 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 - meso_tn3[0] = meso_tn2[3] - - if inputp['alt'] <= zn3[0]: - - # lower stratosphere and troposphere (below zn3[0]) - # Temperature at nodes and gradients at end nodes - # Inverse temperature a linear function of spherical harmonics - - meso_tgn3[0] = meso_tgn2[1] - meso_tn3[1] = pma[3,0]*pavgm[3]/(1-flags['sw'][21]*glob7s(pma[3], inputp, flags,varli)) - meso_tn3[2] = pma[4,0]*pavgm[4]/(1-flags['sw'][21]*glob7s(pma[4], inputp, flags,varli)) - meso_tn3[3] = pma[5,0]*pavgm[5]/(1-flags['sw'][21]*glob7s(pma[5], inputp, flags,varli)) - meso_tn3[4] = pma[6,0]*pavgm[6]/(1-flags['sw'][21]*glob7s(pma[6], inputp, flags,varli)) - 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 - - # linear transition to full mixing below znz[0] - - dmc = 0 - if inputp['alt'] > zmix: - dmc = 1 - (zn2[0]-inputp['alt'])/(zn2[0] - zmix) - dz28 = soutput['d']['N2'] - - # N2 density - dmr = soutput['d']['N2'] / dm28m - 1 - output['d']['N2'],tz = densm(inputp['alt'],dm28m,xmm, tz, zn3, meso_tn3, meso_tgn3, zn2, meso_tn2, meso_tgn2,gsurf,re) - output['d']['N2'] = output['d']['N2'] * (1 + dmr*dmc) - - # HE density - dmr = soutput['d']['He'] / (dz28 * pdm[0,1]) - 1 - output['d']['He'] = output['d']['N2'] * pdm[0,1] * (1 + dmr*dmc) - - # O density - output['d']['O'] = 0 - output['d']['ANM O'] = 0 - - # O2 density - dmr = soutput['d']['O2'] / (dz28 * pdm[3,1]) - 1 - output['d']['O2'] = output['d']['N2'] * pdm[3,1] * (1 + dmr*dmc) - - # AR density - dmr = soutput['d']['AR'] / (dz28 * pdm[4,1]) - 1 - output['d']['AR'] = output['d']['N2'] * pdm[4,1] * (1 + dmr*dmc) - - # Hydrogen density - output['d']['H'] = 0 - - # Atomic nitrogen density - output['d']['N'] = 0 - - # Total mass density - output['d']['RHO'] = 1.66E-24 * (4 * output['d']['He'] + 16 * output['d']['O'] + 28 * output['d']['N2']\ - + 32 * output['d']['O2'] + 40 * output['d']['AR'] + output['d']['H'] + 14 * output['d']['N']) - - output['d']['RHO'] = output['d']['RHO']/1000 - - # temperature at altitude - dd,tz = densm(inputp['alt'], 1, 0, tz, zn3, meso_tn3, meso_tgn3, zn2, meso_tn2, meso_tgn2,gsurf,re) - output['t']['TG'] = tz - return output - -# =============================== gtd7d ============================= # - -def gtd7d(inputp, flags): - output = gtd7(inputp, flags) - output['d']['RHO'] = 1.66E-24 * (4 * output['d']['He'] + 16 * output['d']['O'] + 28 * output['d']['N2']\ - + 32 * output['d']['O2'] + 40 * output['d']['AR'] + output['d']['H'] + 14 * output['d']['N'] + 16 * output['d']['ANM O']) - - output['d']['RHO'] = output['d']['RHO']/1e3 - return output - -# =============================== gts7 ------------------------------ # - -def gts7(inputp,flags,gsurf,re): - # Thermospheric portion of NRLMSISE-00 - # See gtd7 for more extensive comments - # alt > 72.5 km! - # call: nrlmsis00_data, globe7, densu - - output = {'d':{'He':0,'O':0,'N2':0,'O2':0,'AR':0,'RHO':0,'H':0,'N':0,'ANM O':0},\ - 't':{'TINF':0,'TG':0}} - tz = 0 - dm28 = 0 - meso_tn1,meso_tn3 = [np.zeros(5) for i in range(2)] - meso_tn2 = np.zeros(4) - meso_tgn1,meso_tgn2,meso_tgn3 = [np.zeros(2) for i in range(3)] - - zn1 = np.array([120.0, 110.0, 100.0, 90.0, 72.5]) - - dr = 1.72142E-2 - alpha = np.array([-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0]) - altl = np.array([200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0]) - pt,pd,ps,pdl,ptm,pdm,ptl,pma,sam,pavgm = nrlmsis00_data() - za = pdl[1,15] - zn1[0] = za - - # tinf variations not important below za or zn1[0] - if inputp['alt'] > zn1[0]: - tinf_tmp,varli = globe7(pt,inputp,flags) - tinf = ptm[0]*pt[0] * (1+flags['sw'][15]*tinf_tmp) - else: - tinf = ptm[0]*pt[0] - output['t']['TINF'] = tinf - - # gradient variations not important below zn1[4] - if inputp['alt'] > zn1[4]: - tinf_tmp,varli = globe7(ps,inputp,flags) - grad = ptm[3]*ps[0] * (1+flags['sw'][18]*tinf_tmp) - else: - grad = ptm[3]*ps[0] - tinf_tmp,varli = globe7(pd[3],inputp,flags) - tlb = ptm[1] * (1 + flags['sw'][16]*tinf_tmp)*pd[3,0] - s = grad/(tinf - tlb) - - # Lower thermosphere temp variations not significant for density above 300 km - if inputp['alt'] < 300: - meso_tn1[1] = ptm[6]*ptl[0,0]/(1.0-flags['sw'][17]*glob7s(ptl[0], inputp, flags,varli)) - meso_tn1[2] = ptm[2]*ptl[1,0]/(1.0-flags['sw'][17]*glob7s(ptl[1], inputp, flags,varli)) - meso_tn1[3] = ptm[7]*ptl[2,0]/(1.0-flags['sw'][17]*glob7s(ptl[2], inputp, flags,varli)) - meso_tn1[4] = ptm[4]*ptl[3,0]/(1.0-flags['sw'][17]*flags['sw'][19]*glob7s(ptl[3], inputp, flags,varli)) - 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 - else: - meso_tn1[1]=ptm[6]*ptl[0,0] - meso_tn1[2]=ptm[2]*ptl[1,0] - meso_tn1[3]=ptm[7]*ptl[2,0] - meso_tn1[4]=ptm[4]*ptl[3,0] - meso_tgn1[1]=ptm[8]*pma[8,0]*meso_tn1[4]*meso_tn1[4]/(ptm[4]*ptl[3,0])**2 - - # N2 variation factor at Zlb - tinf_tmp,varli = globe7(pd[2],inputp,flags) - g28 = flags['sw'][20]*tinf_tmp - - # variation of turbopause height - 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]))) - output['t']['TINF'] = tinf - xmm = pdm[2,4] - z = inputp['alt'] - - # N2 density - # Diffusive density at Zlb - db28 = pdm[2,0]*np.exp(g28)*pd[2,0] - # Diffusive density at Alt - 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) - dd = output['d']['N2'] - # Turbopause - zh28 = pdm[2,2]*zhf - zhm28 = pdm[2,3]*pdl[1,5] - xmd = 28 - xmm - # Mixed density at Zlb - b28,tz = densu(zh28,db28,tinf,tlb,xmd,(alpha[2]-1),tz,ptm[5],s, zn1,meso_tn1,meso_tgn1,gsurf,re) - if flags['sw'][14] and z <= altl[2]: - # Mixed density at Alt - dm28,tz = densu(z,b28,tinf,tlb,xmm,alpha[2],tz,ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re) - # Net density at Alt - output['d']['N2'] = dnet(output['d']['N2'],dm28,zhm28,xmm,28) - - # HE density - # Density variation factor at Zlb - tinf_tmp,varli = globe7(pd[0],inputp,flags) - g4 = flags['sw'][20]*tinf_tmp - # Diffusive density at Zlb - db04 = pdm[0,0]*np.exp(g4)*pd[0,0] - # Diffusive density at Alt - 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) - dd = output['d']['He'] - if flags['sw'][14] and z>> t = '2019-08-20 23:10:59' - >>> lat,lon,alt = 3,5,900 - >>> params,rho,T,nd = nrlmsise00(t,(lat,lon,alt),swdata) - >>> print(params) - >>> print(rho) - >>> print(T) - >>> print(nd) + NRLMSISE-00 is a semi-empirical, global reference atmospheric model of the Earth from ground level to the exosphere up to 2000 km. + A primary use of this model is to aid predictions of satellite orbital decay due to the atmospheric drag. Usage: - params,rho,T,nd = nrlmsise00(t,(lat,lon,alt),swdata,omode='Oxygen',aphmode = 'Aph') + nrl00 = nrlmsise00(t,(lat,lon,alt),swdata,aphmode = 'Aph') Inputs: t -> [str] time(UTC) - location -> [tuple/list] geodetic latitude, longitude in [degree], and altitude in [km] + location -> [tuple/list] geodetic latitude[degree], longitude[degree], and altitude[km] Parameters: - anmomode -> [bool] whether to add anomalous oxygen component to the density; defalut = True - aphmode -> [bool] whether to use the 3h geomagnetic index, default = True + aphmode -> [bool, optional, default = True] whether to use the 3h geomagnetic index Output: - params -> [dict] parameters involved in computations, including Year, DayOfYear(DOY), SecondOfDay(SOD), - Latitude[deg], Longitude[deg], Altitude[km], LocalSolarTime(LST), f107Average(f107A) [10^-22 W/m^2/Hz] - f107Daily(f107D) [10^-22 W/m^2/Hz], ApDaily(ApD), Ap3Hourly(Ap3H) - rho -> [float] total mass density, [kg/m^3] - T -> [tuple] exospheric and local temperature, [K] - nd -> [dict] number density of components, [1/m^3]; mainly include Helium(He), Oxygen(O), Oxygen(O2), Nitrogen(N), - Nitrogen(N2), Argon(Ar), Hydrogen(H), Anomalous Oxygen(ANM O) + nrl00 -> instance of class ATMOS, where its attributes include + rho -> [float] total mass density[kg/m^3] + T -> [tuple] local temperature[K] + nd -> [dict] number density of components[1/m^3], including Helium(He), Oxygen(O), Oxygen(O2), Nitrogen(N), + Nitrogen(N2), Argon(Ar), Hydrogen(H), Anomalous Oxygen(ANM O) + + Examples: + >>> from pyatmos import download_sw_nrlmsise00,read_sw_nrlmsise00 + >>> # Download or update the space weather file from www.celestrak.com + >>> swfile = download_sw_nrlmsise00() + >>> # Read the space weather data + >>> swdata = read_sw_nrlmsise00(swfile) + >>> + >>> from pyatmos import nrlmsise00 + >>> # Set a specific time and location + >>> t = '2014-07-22 22:18:45' # time(UTC) + >>> lat,lon,alt = 25,102,600 # latitude, longitude in [degree], and altitude in [km] + >>> nrl00 = nrlmsise00(t,(lat,lon,alt),swdata) + >>> print(nrl00.rho) # [kg/m^3] + >>> print(nrl00.T) # [K] + >>> print(nrl00.nd) # composition in [1/m^3] """ lat,lon,h = location @@ -1079,13 +61,14 @@ def nrlmsise00(t,location,SW_OBS_PRE,anmomode=True,aphmode=True): f107A,f107,ap,aph = get_sw(SW_OBS_PRE,t_ymd,hour) else: f107A,f107,ap,aph = 150,150,4,np.full(7,4) + inputp = {'doy':doy,'year':year,'sec':sec,'alt':alt,'g_lat':lat,'g_lon':lon_wrap,'lst':lst,\ 'f107A':f107A,'f107':f107,'ap':ap,'ap_a':aph} switches = np.ones(23) if aphmode: switches[8] = -1 # -1 indicates the use of 3h geomagnetic index - if anmomode: + if alt > 500: output = gtd7d(inputp,switches) else: output = gtd7(inputp,switches) @@ -1096,4 +79,7 @@ def nrlmsise00(t,location,SW_OBS_PRE,anmomode=True,aphmode=True): rho = output['d']['RHO'] T = (output['t']['TINF'],output['t']['TG']) nd = {'He':output['d']['He'],'O':output['d']['O'],'N2':output['d']['N2'],'O2':output['d']['O2'],'Ar':output['d']['AR'],'H':output['d']['H'],'N':output['d']['N'],'ANM O':output['d']['ANM O']} - return params,rho,T,nd + + info = {'rho':rho,'T':output['t']['TG'],'nd':nd} + + return ATMOS(info) diff --git a/pyatmos/msise/nrlmsise00_subfunc.py b/pyatmos/msise/nrlmsise00_subfunc.py new file mode 100644 index 0000000..0e9d224 --- /dev/null +++ b/pyatmos/msise/nrlmsise00_subfunc.py @@ -0,0 +1,879 @@ +import numpy as np +from scipy.interpolate import CubicSpline +from pyshtools.legendre import PLegendreA,PlmIndex +import pkg_resources + +def nrlmsis00_data(): + ''' + Read the data block from nrlmsis00_data.npz + ''' + data_path = pkg_resources.resource_filename('pyatmos', 'data/') + data = np.load(data_path+'nrlmsis00_data.npz') + pt,pd,ps,pdl = data['pt'],data['pd'],data['ps'],data['pdl'] + ptm,pdm,ptl,pma = data['ptm'],data['pdm'],data['ptl'],data['pma'] + sam,pavgm = data['sam'],data['pavgm'] + return pt,pd,ps,pdl,ptm,pdm,ptl,pma,sam,pavgm + +def tselec(switches): + flags = {'sw':np.zeros(23),'swc':np.zeros(23)} + for i in range(23): + if i != 8: + if switches[i] == 1: + flags['sw'][i] = 1 + else: + flags['sw'][i] = 0 + if switches[i] > 0: + flags['swc'][i] = 1 + else: + flags['swc'][i] = 0 + else: + flags['sw'][i] = switches[i] + flags['swc'][i] = switches[i] + return flags + +def glatf(lat): + c2 = np.cos(2*np.deg2rad(lat)) + gv = 980.616*(1 - 0.0026373*c2) + reff = 2*gv/(3.085462E-6 + 2.27E-9*c2)*1E-5 + return gv,reff + +def ccor(alt,r,h1,zh): + e = (alt - zh)/h1 + if e > 70: + return 1 + elif e < -70: + return np.exp(r) + else: + return np.exp(r/(1 + np.exp(e))) + +def ccor2(alt,r,h1,zh,h2): + e1 = (alt - zh)/h1 + e2 = (alt - zh)/h2 + if e1 > 70 or e2 > 70: + return 1 + if e1 < -70 and e2 < -70: + return np.exp(r) + ex1,ex2 = np.exp([e1,e2]) + ccor2v = r/(1 + 0.5*(ex1 + ex2)) + return np.exp(ccor2v) + +def scalh(alt,xm,temp,gsurf,re): + rgas = 831.4 + g = rgas*temp/(gsurf/(1 + alt/re)**2*xm) + return g + +def dnet(dd,dm,zhm,xmm,xm): + ''' + Turbopause correction for msis models + ''' + a = zhm/(xmm - xm) + if not (dm > 0 and dd > 0): + print('dnet log error {0:.1f} {1:.1f} {2:.1f}'.format(dm,dd,xm)) + if dd == 0 and dm == 0: dd = 1 + if dm == 0: return dd + if dd == 0: return dm + ylog = a*np.log(dm/dd) + if ylog < -10: return dd + if ylog > 10: return dm + a = dd*(1 + np.exp(ylog))**(1/a) + return a + +def zeta(zz,zl,re): + return (zz - zl)*(re + zl)/(re + zz) + +def densm(alt, d0, xm, tz, zn3, tn3, tgn3, zn2, tn2, tgn2,gsurf,re): + rgas = 831.4 + densm_tmp = d0 + tz_tmp = tz + mn3,mn2 = len(zn3),len(zn2) + + if alt > zn2[0]: + if xm == 0: + densm_tmp = tz + return densm_tmp,tz_tmp + else: + densm_tmp = d0 + return densm_tmp,tz_tmp + + # stratosphere/mesosphere temperature + if alt > zn2[mn2-1]: + z = alt + else: + z = zn2[mn2-1] + mn = mn2 + xs,ys = [np.zeros(mn) for i in range(2)] + z1,z2 = zn2[0],zn2[mn-1] + t1,t2=tn2[0],tn2[mn-1] + zg,zgdif = zeta(z,z1,re),zeta(z2,z1,re) + + # set up spline nodes + for k in range(mn): + xs[k] = zeta(zn2[k],z1,re)/zgdif + ys[k] = 1/tn2[k] + yd1 = -tgn2[0]/t1**2*zgdif + yd2 = -tgn2[1]/t2**2*zgdif*((re + z2)/(re + z1))**2 + + # calculate spline coefficients + cs = CubicSpline(xs,ys,bc_type=((1,yd1),(1,yd2))) + x = zg/zgdif + y = cs(x) + + # temperature at altitude + tz_tmp = 1/y + if xm != 0: + # calaculate stratosphere/mesospehere density + glb = gsurf/(1 + z1/re)**2 + gamm = xm*glb*zgdif/rgas + + # Integrate temperature profile + yi = cs.integrate(xs[0],x) + expl = gamm*yi + if expl > 50: + expl = 50 + # Density at altitude + densm_tmp = densm_tmp*(t1/tz_tmp)*np.exp(-expl) + if alt > zn3[0]: + if xm == 0: + densm_tmp = tz_tmp + return densm_tmp,tz_tmp + else: + return densm_tmp,tz_tmp + + # troposhere/stratosphere temperature + z = alt + mn = mn3 + xs,ys = [np.zeros(mn) for i in range(2)] + z1,z2 = zn3[0],zn3[mn-1] + t1,t2 = tn3[0],tn3[mn-1] + zg,zgdif = zeta(z,z1,re),zeta(z2,z1,re) + + # set up spline nodes + for k in range(mn): + xs[k] = zeta(zn3[k],z1,re)/zgdif + ys[k] = 1/tn3[k] + yd1 = -tgn3[0]/t1**2*zgdif + yd2 = -tgn3[1]/t2**2*zgdif*((re+z2)/(re+z1))**2 + + # calculate spline coefficients + cs = CubicSpline(xs,ys,bc_type=((1,yd1),(1,yd2))) + x = zg/zgdif + y = cs(x) + + # temperature at altitude + tz_tmp = 1/y + if xm != 0: + # calaculate tropospheric / stratosphere density + glb = gsurf/(1 + z1/re)**2 + gamm = xm*glb*zgdif/rgas + + # Integrate temperature profile + yi = cs.integrate(xs[0],x) + expl = gamm*yi + if expl > 50: expl = 50 + # Density at altitude + densm_tmp = densm_tmp*(t1/tz_tmp)*np.exp(-expl) + + if xm == 0: + densm_tmp = tz_tmp + return densm_tmp,tz_tmp + else: + return densm_tmp,tz_tmp + +def densu (alt,dlb,tinf,tlb,xm,alpha,tz,zlb,s2,zn1,tn1,tgn1,gsurf,re): + ''' + Calculate Temperature and Density Profiles for MSIS models + ''' + rgas = 831.4 + densu_tmp = 1 + mn1 = len(zn1) + # joining altitudes of Bates and spline + za = zn1[0] + if alt > za: + z = alt + else: + z = za + # geopotential altitude difference from ZLB + zg2 = zeta(z,zlb,re) + + # Bates temperature + tt = tinf - (tinf - tlb)*np.exp(-s2*zg2) + ta = tz = tt + densu_tmp = tz_tmp = tz + + if alt < za: + # calculate temperature below ZA + # temperature gradient at ZA from Bates profile + dta = (tinf - ta)*s2*((re + zlb)/(re + za))**2 + tgn1[0],tn1[0] = dta,ta + if alt > zn1[mn1-1]: + z = alt + else: + z = zn1[mn1-1] + mn = mn1 + xs,ys = [np.zeros(mn) for i in range(2)] + z1,z2 = zn1[0],zn1[mn-1] + t1,t2 = tn1[0],tn1[mn-1] + # geopotental difference from z1 + zg,zgdif = zeta(z,z1,re),zeta(z2,z1,re) + # set up spline nodes + for k in range(mn): + xs[k] = zeta(zn1[k],z1,re)/zgdif + ys[k] = 1/tn1[k] + # end node derivatives + yd1 = -tgn1[0]/t1**2*zgdif + yd2 = -tgn1[1]/t2**2*zgdif*((re + z2)/(re + z1))**2 + # calculate spline coefficients + cs = CubicSpline(xs,ys,bc_type=((1,yd1),(1,yd2))) + x = zg/zgdif + y = cs(x) + # temperature at altitude + tz_tmp = 1/y + densu_tmp = tz_tmp + if xm == 0: return densu_tmp,tz_tmp + + # calculate density above za + glb = gsurf/(1 + zlb/re)**2 + gamma = xm*glb/(s2*rgas*tinf) + expl = np.exp(-s2*gamma*zg2) + if expl > 50: expl = 50 + if tt <= 0: expl = 50 + + # density at altitude + densa = dlb*(tlb/tt)**(1 + alpha + gamma)*expl + densu_tmp = densa + if alt >= za: return densu_tmp,tz_tmp + + # calculate density below za + glb = gsurf/(1 + z1/re)**2 + gamm = xm*glb*zgdif/rgas + + # integrate spline temperatures + yi = cs.integrate(xs[0],x) + expl = gamm*yi + if expl > 50: expl = 50 + if tz_tmp <= 0: expl = 50 + + # density at altitude + densu_tmp = densu_tmp*(t1/tz_tmp)**(1 + alpha)*np.exp(-expl) + return densu_tmp,tz_tmp + +# =============== 3hr magnetic activity functions =================== # + +# Eq. A24d +def g0(a,p): + return (a - 4 + (p[25] - 1)*(a - 4 + (np.exp(-np.abs(p[24])*(a - 4)) - 1) / np.abs(p[24]))) + +# Eq. A24c +def sumex(ex): + return (1 + (1 - ex**19)/(1 - ex)*ex**0.5) + +# Eq. A24a +def sg0(ex,p,ap): + # call sumex, g0 + return (g0(ap[1],p) + g0(ap[2],p)*ex + g0(ap[3],p)*ex**2 + \ + g0(ap[4],p)*ex**3 + (g0(ap[5],p)*ex**4 + \ + g0(ap[6],p)*ex**12)*(1-ex**8)/(1-ex))/sumex(ex) + +# =============== 3hr magnetic activity functions =================== # + +def lengendre(g_lat,lmax = 8): + # Index of PLegendreA_x can be calculated by PlmIndex(l,m) + x = np.sin(np.deg2rad(g_lat)) + PLegendreA_x = PLegendreA(lmax,x) + return PLegendreA_x + +def globe7(p,inputp,flags): + ''' + Calculate G(L) function + ''' + t = np.zeros(15) + sr = 7.2722E-5 + dr = 1.72142E-2 + hr = 0.2618 + + apdf = 0 + apt = np.zeros(4) + tloc = inputp['lst'] + + if not (flags['sw'][6]==0 and flags['sw'][7]==0 and flags['sw'][13]==0): + stloc,ctloc = np.sin(hr*tloc),np.cos(hr*tloc) + s2tloc,c2tloc = np.sin(2*hr*tloc),np.cos(2*hr*tloc) + s3tloc,c3tloc = np.sin(3*hr*tloc),np.cos(3*hr*tloc) + cd32 = np.cos(dr*(inputp['doy'] - p[31])) + cd18 = np.cos(2*dr*(inputp['doy'] - p[17])) + cd14 = np.cos(dr*(inputp['doy'] - p[13])) + cd39 = np.cos(2*dr*(inputp['doy'] - p[38])) + + # F10.7 effect + df = inputp['f107'] - inputp['f107A'] + dfa = inputp['f107A'] - 150 + t[0] = p[19]*df*(1 + p[59]*dfa) + p[20]*df**2 + p[21]*dfa + p[29]*dfa**2 + f1 = 1 + (p[47]*dfa + p[19]*df + p[20]*df**2)*flags['swc'][0] + f2 = 1 + (p[49]*dfa + p[19]*df + p[20]*df**2)*flags['swc'][0] + + plg = lengendre(inputp['g_lat']) + + # time independent + 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] + + # symmetrical annual + t[2] = p[18]*cd32 + + # symmetrical semiannual + t[3] = (p[15] + p[16]*plg[3])*cd18 + + # asymmetrical annual + t[4] = f1*(p[9]*plg[1] + p[10]*plg[6])*cd14 + + # asymmetrical semiannual + t[5] = p[37]*plg[1]*cd39 + + # diurnal + if flags['sw'][6]: + t71 = p[11]*plg[4]*cd14*flags['swc'][4] + t72 = p[12]*plg[4]*cd14*flags['swc'][4] + 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) + + # semiannual + if flags['sw'][7]: + t81 = (p[23]*plg[8] + p[35]*plg[17])*cd14*flags['swc'][4] + t82 = (p[33]*plg[8] + p[36]*plg[17])*cd14*flags['swc'][4] + t[7] = f2*((p[5]*plg[5] + p[41]*plg[12] + t81)*c2tloc +(p[8]*plg[5] + p[42]*plg[12] + t82)*s2tloc) + + # terdiurnal + if flags['sw'][13]: + 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) + + # magnetic activity based on daily ap + if flags['sw'][8] == -1: + ap = inputp['ap_a'] + if p[51]!= 0: + exp1 = np.exp(-10800*np.abs(p[51])/(1 + p[138]*(45 - np.abs(inputp['g_lat'])))) + if exp1 > 0.99999: exp1 = 0.99999 + if p[24] < 1E-4: p[24] = 1E-4 + apt[0] = sg0(exp1,p,ap) + # apt[1] = sg2(exp1,p,ap) + # apt[2] = sg0(exp2,p,ap) + # apt[3] = sg2(exp2,p,ap) + + if flags['sw'][8]: + t[8] = apt[0]*(p[50] + p[96]*plg[3] + p[54]*plg[10] + \ + (p[125]*plg[1] + p[126]*plg[6] + p[127]*plg[15])*cd14*flags['swc'][4] + \ + (p[128]*plg[2] + p[129]*plg[7] + p[130]*plg[16])*flags['swc'][6]*np.cos(hr*(tloc - p[131]))) + else: + apd = inputp['ap'] - 4 + p44 = p[43] + p45 = p[44] + if p44 < 0: p44 = 1E-5 + apdf = apd + (p45 - 1)*(apd + (np.exp(-p44*apd) - 1)/p44) + if flags['sw'][8]: + t[8]=apdf*(p[32] + p[45]*plg[3] + p[34]*plg[10] + \ + (p[100]*plg[1] + p[101]*plg[6] + p[102]*plg[15])*cd14*flags['swc'][4] + + (p[121]*plg[2] + p[122]*plg[7] + p[123]*plg[16])*flags['swc'][6]*np.cos(hr*(tloc - p[124]))) + + if flags['sw'][9] and inputp['g_lon'] > -1000: + # longitudinal + if flags['sw'][10]: + t[10] = (1 + p[80]*dfa*flags['swc'][0])*((p[64]*plg[4] + p[65]*plg[11] + p[66]*plg[22]\ + + p[103]*plg[2] + p[104]*plg[7] + p[105]*plg[16]\ + + flags['swc'][4]*(p[109]*plg[2] + p[110]*plg[7] + p[111]*plg[16])*cd14)*np.cos(np.deg2rad(inputp['g_lon'])) \ + +(p[90]*plg[4]+p[91]*plg[11]+p[92]*plg[22] + p[106]*plg[2]+p[107]*plg[7]+p[108]*plg[16]\ + + flags['swc'][4]*(p[112]*plg[2] + p[113]*plg[7] + p[114]*plg[16])*cd14)*np.sin(np.deg2rad(inputp['g_lon']))) + + # ut and mixed ut, longitude + if flags['sw'][11]: + t[11]=(1 + p[95]*plg[1])*(1 + p[81]*dfa*flags['swc'][0])*\ + (1 + p[119]*plg[1]*flags['swc'][4]*cd14)*\ + ((p[68]*plg[1] + p[69]*plg[6] + p[70]*plg[15])*np.cos(sr*(inputp['sec'] - p[71]))) + t[11] += flags['swc'][10]*(p[76]*plg[8] + p[77]*plg[17] + p[78]*plg[30])*\ + np.cos(sr*(inputp['sec'] - p[79]) + 2*np.deg2rad(inputp['g_lon']))*(1 + p[137]*dfa*flags['swc'][0]) + + # ut, longitude magnetic activity + if flags['sw'][10]: + if flags['sw'][8] == -1: + if p[51]: + t[12] = apt[0]*flags['swc'][10]*(1 + p[132]*plg[1])*\ + ((p[52]*plg[4] + p[98]*plg[11] + p[67]*plg[22])* np.cos(np.deg2rad(inputp['g_lon'] - p[97])))\ + + apt[0]*flags['swc'][10]*flags['swc'][4]*(p[133]*plg[2] + p[134]*plg[7] + p[135]*plg[16])*\ + cd14*np.cos(np.deg2rad(inputp['g_lon'] - p[136])) + apt[0]*flags['swc'][11]* \ + (p[55]*plg[1] + p[56]*plg[6] + p[57]*plg[15])*np.cos(sr*(inputp['sec'] - p[58])) + else: + t[12] = apdf*flags['swc'][10]*(1 + p[120]*plg[1])*((p[60]*plg[4] + p[61]*plg[11] + p[62]*plg[22])*\ + np.cos(np.deg2rad(inputp['g_lon']-p[63])))+apdf*flags['swc'][10]*flags['swc'][4]* \ + (p[115]*plg[2] + p[116]*plg[7] + p[117]*plg[16])* \ + cd14*np.cos(np.deg2rad(inputp['g_lon'] - p[118])) \ + + apdf*flags['swc'][11]*(p[83]*plg[1] + p[84]*plg[6] + p[85]*plg[15])* np.cos(sr*(inputp['sec'] - p[75])) + + # parms not used: 82, 89, 99, 139-149 + tinf = p[30] + for i in range(14): + tinf = tinf + np.abs(flags['sw'][i])*t[i] + return tinf,[dfa,plg,ctloc,stloc,c2tloc,s2tloc,s3tloc,c3tloc,apdf,apt] + +def glob7s(p,inputp,flags,varli): + pset = 2 + t = np.zeros(14) + dr = 1.72142E-2 + [dfa,plg,ctloc,stloc,c2tloc,s2tloc,s3tloc,c3tloc,apdf,apt] = varli + + # confirm parameter set + if p[99] == 0: p[99] = pset + if p[99] != pset: + print("Wrong parameter set for glob7s") + return -1 + + for j in range(14): + t[j] = 0 + cd32 = np.cos(dr*(inputp['doy'] - p[31])) + cd18 = np.cos(2*dr*(inputp['doy'] - p[17])) + cd14 = np.cos(dr*(inputp['doy'] - p[13])) + cd39 = np.cos(2*dr*(inputp['doy'] - p[38])) + + # F10.7 + t[0] = p[21]*dfa + + # time independent + 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] + + # symmetrical annual + t[2] = (p[18] + p[47]*plg[3] + p[29]*plg[10])*cd32 + + # symmetrical semiannual + t[3] = (p[15] + p[16]*plg[3] + p[30]*plg[10])*cd18 + + # asymmetrical annual + t[4] = (p[9]*plg[1] + p[10]*plg[6] + p[20]*plg[15])*cd14 + + # asymmetrical semiannual + t[5] = p[37]*plg[1]*cd39; + + # diurnal + if flags['sw'][6]: + t71 = p[11]*plg[4]*cd14*flags['swc'][4] + t72 = p[12]*plg[4]*cd14*flags['swc'][4] + t[6] = ((p[3]*plg[2] + p[4]*plg[7] + t71)*ctloc + (p[6]*plg[2] + p[7]*plg[7] + t72)*stloc) + + # semidiurnal + if flags['sw'][7]: + t81 = (p[23]*plg[8] + p[35]*plg[17])*cd14*flags['swc'][4] + t82 = (p[33]*plg[8] + p[36]*plg[17])*cd14*flags['swc'][4] + t[7] = ((p[5]*plg[5] + p[41]*plg[12] + t81)*c2tloc + (p[8]*plg[5] + p[42]*plg[12] + t82)*s2tloc) + + # terdiurnal + if flags['sw'][13]: + t[13] = p[39]*plg[9]*s3tloc + p[40]*plg[9]*c3tloc + + # magnetic activity + if flags['sw'][8]: + if flags['sw'][8]==1: + t[8] = apdf * (p[32] + p[45]*plg[3]*flags['swc'][1]) + if flags['sw'][8]==-1: + t[8]=(p[50]*apt[0] + p[96]*plg[3]*apt[0]*flags['swc'][1]) + + # longitudinal + if not (flags['sw'][9]==0 or flags['sw'][10]==0 or inputp['g_lon']<=-1000): + t[10] = (1 + plg[1]*(p[80]*flags['swc'][4]*np.cos(dr*(inputp['doy'] - p[81]))\ + + p[85]*flags['swc'][5]*np.cos(2*dr*(inputp['doy'] - p[86])))\ + + p[83]*flags['swc'][2]*np.cos(dr*(inputp['doy'] - p[84]))\ + + p[87]*flags['swc'][3]*np.cos(2*dr*(inputp['doy'] - p[88])))\ + *((p[64]*plg[4] + p[65]*plg[11] + p[66]*plg[22]\ + + p[74]*plg[2] + p[75]*plg[7] + p[76]*plg[16])*np.cos(np.deg2rad(inputp['g_lon']))\ + + (p[90]*plg[4] + p[91]*plg[11] + p[92]*plg[22]\ + + p[77]*plg[2] + p[78]*plg[7] + p[79]*plg[16])*np.sin(np.deg2rad(inputp['g_lon']))) + + tt = 0 + for i in range(14): + tt += np.abs(flags['sw'][i])*t[i] + return tt + +def gtd7(inputp,switches): + tz = 0 + zn3 = np.array([32.5,20.0,15.0,10.0,0.0]) + zn2 = np.array([72.5,55.0,45.0,32.5]) + zmix= 62.5 + + output = {'d':{'He':0,'O':0,'N2':0,'O2':0,'AR':0,'RHO':0,'H':0,'N':0,'ANM O':0},\ + 't':{'TINF':0,'TG':0}} + + flags = tselec(switches) + + # Latitude variation of gravity (none for sw[1]=0) + xlat = inputp['g_lat'] + if flags['sw'][1]==0: xlat = 45 + gsurf,re = glatf(xlat) + pt,pd,ps,pdl,ptm,pdm,ptl,pma,sam,pavgm = nrlmsis00_data() + xmm = pdm[2,4] + + # thermosphere/mesosphere (above zn2[0]) + if inputp['alt'] > zn2[0]: + altt = inputp['alt'] + else: + altt = zn2[0] + + tmp = inputp['alt'] + inputp['alt'] = altt + 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) + altt = inputp['alt'] + inputp['alt'] = tmp + # metric adjustment + dm28m = dm28*1E6 + output['t']['TINF'] = soutput['t']['TINF'] + output['t']['TG'] = soutput['t']['TG'] + if inputp['alt'] >= zn2[0]: + output['d'] = soutput['d'] + return output + + varli = [dfa,plg,ctloc,stloc,c2tloc,s2tloc,s3tloc,c3tloc,apdf,apt] + + meso_tgn2[0] = meso_tgn1[1] + meso_tn2[0] = meso_tn1[4] + meso_tn2[1] = pma[0,0]*pavgm[0]/(1-flags['sw'][19]*glob7s(pma[0], inputp, flags,varli)) + meso_tn2[2] = pma[1,0]*pavgm[1]/(1-flags['sw'][19]*glob7s(pma[1], inputp, flags,varli)) + meso_tn2[3] = pma[2,0]*pavgm[2]/(1-flags['sw'][19]*flags['sw'][21]*glob7s(pma[2], inputp, flags,varli)) + 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 + meso_tn3[0] = meso_tn2[3] + + if inputp['alt'] <= zn3[0]: + + meso_tgn3[0] = meso_tgn2[1] + meso_tn3[1] = pma[3,0]*pavgm[3]/(1-flags['sw'][21]*glob7s(pma[3], inputp, flags,varli)) + meso_tn3[2] = pma[4,0]*pavgm[4]/(1-flags['sw'][21]*glob7s(pma[4], inputp, flags,varli)) + meso_tn3[3] = pma[5,0]*pavgm[5]/(1-flags['sw'][21]*glob7s(pma[5], inputp, flags,varli)) + meso_tn3[4] = pma[6,0]*pavgm[6]/(1-flags['sw'][21]*glob7s(pma[6], inputp, flags,varli)) + 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 + + # linear transition to full mixing below znz[0] + + dmc = 0 + if inputp['alt'] > zmix: + dmc = 1 - (zn2[0]-inputp['alt'])/(zn2[0] - zmix) + dz28 = soutput['d']['N2'] + + # N2 density + dmr = soutput['d']['N2'] / dm28m - 1 + output['d']['N2'],tz = densm(inputp['alt'],dm28m,xmm, tz, zn3, meso_tn3, meso_tgn3, zn2, meso_tn2, meso_tgn2,gsurf,re) + output['d']['N2'] = output['d']['N2'] * (1 + dmr*dmc) + + # HE density + dmr = soutput['d']['He'] / (dz28 * pdm[0,1]) - 1 + output['d']['He'] = output['d']['N2'] * pdm[0,1] * (1 + dmr*dmc) + + # O density + output['d']['O'] = 0 + output['d']['ANM O'] = 0 + + # O2 density + dmr = soutput['d']['O2'] / (dz28 * pdm[3,1]) - 1 + output['d']['O2'] = output['d']['N2'] * pdm[3,1] * (1 + dmr*dmc) + + # AR density + dmr = soutput['d']['AR'] / (dz28 * pdm[4,1]) - 1 + output['d']['AR'] = output['d']['N2'] * pdm[4,1] * (1 + dmr*dmc) + + # Hydrogen density + output['d']['H'] = 0 + + # Atomic nitrogen density + output['d']['N'] = 0 + + # Total mass density + output['d']['RHO'] = 1.66E-24 * (4 * output['d']['He'] + 16 * output['d']['O'] + 28 * output['d']['N2']\ + + 32 * output['d']['O2'] + 40 * output['d']['AR'] + output['d']['H'] + 14 * output['d']['N']) + + output['d']['RHO'] = output['d']['RHO']/1000 + + # temperature at altitude + dd,tz = densm(inputp['alt'], 1, 0, tz, zn3, meso_tn3, meso_tgn3, zn2, meso_tn2, meso_tgn2,gsurf,re) + output['t']['TG'] = tz + return output + + +def gtd7d(inputp, flags): + output = gtd7(inputp, flags) + output['d']['RHO'] = 1.66E-24 * (4 * output['d']['He'] + 16 * output['d']['O'] + 28 * output['d']['N2']\ + + 32 * output['d']['O2'] + 40 * output['d']['AR'] + output['d']['H'] + 14 * output['d']['N'] + 16 * output['d']['ANM O']) + + output['d']['RHO'] = output['d']['RHO']/1e3 + return output + +def gts7(inputp,flags,gsurf,re): + + output = {'d':{'He':0,'O':0,'N2':0,'O2':0,'AR':0,'RHO':0,'H':0,'N':0,'ANM O':0},\ + 't':{'TINF':0,'TG':0}} + tz = 0 + dm28 = 0 + meso_tn1,meso_tn3 = [np.zeros(5) for i in range(2)] + meso_tn2 = np.zeros(4) + meso_tgn1,meso_tgn2,meso_tgn3 = [np.zeros(2) for i in range(3)] + + zn1 = np.array([120.0, 110.0, 100.0, 90.0, 72.5]) + + dr = 1.72142E-2 + alpha = np.array([-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0]) + altl = np.array([200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0]) + pt,pd,ps,pdl,ptm,pdm,ptl,pma,sam,pavgm = nrlmsis00_data() + za = pdl[1,15] + zn1[0] = za + + # tinf variations not important below za or zn1[0] + if inputp['alt'] > zn1[0]: + tinf_tmp,varli = globe7(pt,inputp,flags) + tinf = ptm[0]*pt[0] * (1+flags['sw'][15]*tinf_tmp) + else: + tinf = ptm[0]*pt[0] + output['t']['TINF'] = tinf + + # gradient variations not important below zn1[4] + if inputp['alt'] > zn1[4]: + tinf_tmp,varli = globe7(ps,inputp,flags) + grad = ptm[3]*ps[0] * (1+flags['sw'][18]*tinf_tmp) + else: + grad = ptm[3]*ps[0] + tinf_tmp,varli = globe7(pd[3],inputp,flags) + tlb = ptm[1] * (1 + flags['sw'][16]*tinf_tmp)*pd[3,0] + s = grad/(tinf - tlb) + + # Lower thermosphere temp variations not significant for density above 300 km + if inputp['alt'] < 300: + meso_tn1[1] = ptm[6]*ptl[0,0]/(1.0-flags['sw'][17]*glob7s(ptl[0], inputp, flags,varli)) + meso_tn1[2] = ptm[2]*ptl[1,0]/(1.0-flags['sw'][17]*glob7s(ptl[1], inputp, flags,varli)) + meso_tn1[3] = ptm[7]*ptl[2,0]/(1.0-flags['sw'][17]*glob7s(ptl[2], inputp, flags,varli)) + meso_tn1[4] = ptm[4]*ptl[3,0]/(1.0-flags['sw'][17]*flags['sw'][19]*glob7s(ptl[3], inputp, flags,varli)) + 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 + else: + meso_tn1[1]=ptm[6]*ptl[0,0] + meso_tn1[2]=ptm[2]*ptl[1,0] + meso_tn1[3]=ptm[7]*ptl[2,0] + meso_tn1[4]=ptm[4]*ptl[3,0] + meso_tgn1[1]=ptm[8]*pma[8,0]*meso_tn1[4]*meso_tn1[4]/(ptm[4]*ptl[3,0])**2 + + # N2 variation factor at Zlb + tinf_tmp,varli = globe7(pd[2],inputp,flags) + g28 = flags['sw'][20]*tinf_tmp + + # variation of turbopause height + 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]))) + output['t']['TINF'] = tinf + xmm = pdm[2,4] + z = inputp['alt'] + + # N2 density + # Diffusive density at Zlb + db28 = pdm[2,0]*np.exp(g28)*pd[2,0] + # Diffusive density at Alt + 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) + dd = output['d']['N2'] + # Turbopause + zh28 = pdm[2,2]*zhf + zhm28 = pdm[2,3]*pdl[1,5] + xmd = 28 - xmm + # Mixed density at Zlb + b28,tz = densu(zh28,db28,tinf,tlb,xmd,(alpha[2]-1),tz,ptm[5],s, zn1,meso_tn1,meso_tgn1,gsurf,re) + if flags['sw'][14] and z <= altl[2]: + # Mixed density at Alt + dm28,tz = densu(z,b28,tinf,tlb,xmm,alpha[2],tz,ptm[5],s,zn1,meso_tn1,meso_tgn1,gsurf,re) + # Net density at Alt + output['d']['N2'] = dnet(output['d']['N2'],dm28,zhm28,xmm,28) + + # HE density + # Density variation factor at Zlb + tinf_tmp,varli = globe7(pd[0],inputp,flags) + g4 = flags['sw'][20]*tinf_tmp + # Diffusive density at Zlb + db04 = pdm[0,0]*np.exp(g4)*pd[0,0] + # Diffusive density at Alt + 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) + dd = output['d']['He'] + if flags['sw'][14] and z