Skip to content

Commit

Permalink
simplify
Browse files Browse the repository at this point in the history
  • Loading branch information
@scivision
scivision committed Mar 5, 2018
1 parent 803b41b commit e195488
Show file tree
Hide file tree
Showing 2 changed files with 134 additions and 167 deletions.
18 changes: 10 additions & 8 deletions examples/TimeProfile.py
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,8 @@
#!/usr/bin/env python
""" Time Profile: IRI2016 """
from pyiri2016 import IRI2016,IRI2016Profile
import pyiri2016
from numpy import arange
from datetime import timedelta
from matplotlib.pyplot import figure, show
try:
import ephem
Expand All @@ -11,11 +12,12 @@
hrlim = [0, 24] # [0,24] does whole day(s)
hrstp = 0.25 # time step [decimal hours]
#lat = -11.95; lon = -76.77
lat = 65; lon = -147.5
alts = arange(100, 200, 10)
# %% run
sim = IRI2016Profile(hrlim=hrlim, hrstp=hrstp, lat=lat, lon=lon, alt=150.,
option='time', verbose=False, time='2012-08-21')
glat = 0
glon = 0
alt_km = arange(120, 180, 10)
# %% ru
sim = pyiri2016.timeprofile(('2012-08-21','2012-08-22'),timedelta(hours=0.25),
alt_km,glat,glon)

hrbins = arange(hrlim[0], hrlim[1] + hrstp, hrstp)

Expand All @@ -34,7 +36,7 @@

pn = axs[0]
NmF2 = sim.b[0, index]
NmF1 = IRI2016()._RmNeg(sim.b[2, index])
NmF1 = pyiri2016.IRI2016()._RmNeg(sim.b[2, index])
NmE = sim.b[4, index]
pn.plot(hrbins, NmF2, label='N$_m$F$_2$')
pn.plot(hrbins, NmF1, label='N$_m$F$_1$')
Expand All @@ -48,7 +50,7 @@

pn = axs[1]
hmF2 = sim.b[1, index]
hmF1 = IRI2016()._RmNeg(sim.b[3, index])
hmF1 = pyiri2016.IRI2016()._RmNeg(sim.b[3, index])
hmE = sim.b[5, index]
pn.plot(hrbins, hmF2, label='h$_m$F$_2$')
pn.plot(hrbins, hmF1, label='h$_m$F$_1$')
Expand Down
283 changes: 124 additions & 159 deletions pyiri2016/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,12 @@

def datetimerange(start:datetime, end:datetime, step:timedelta) -> list:
"""like range() for datetime!"""
if isinstance(start,str):
start = parse(start)

if isinstance(end,str):
end = parse(end)

assert isinstance(start, datetime)
assert isinstance(end, datetime)
assert isinstance(step, timedelta)
Expand All @@ -23,122 +29,122 @@ def __init__(self):
self.iriDataFolder = Path(__file__).parent / 'data'


def Switches(self):
"""
IRI switches to turn on/off several options
"""

jf = np.ones(50,dtype=bool)

# i 1 0 standard version
# ------------------------------------------------------------------------
jf[ 1 - 1] = 1; # 1 Ne computed Ne not computed 1
jf[ 2 - 1] = 1; # 2 Te, Ti computed Te, Ti not computed 1
jf[ 3 - 1] = 1; # 3 Ne & Ni computed Ni not computed 1
jf[ 4 - 1] = 0; # 4 B0 - Table option B0 - other models jf(31) 0
jf[ 5 - 1] = 0; # 5 foF2 - CCIR foF2 - URSI 0
jf[ 6 - 1] = 0; # 6 Ni - DS-95 & DY-85 Ni - RBV-10 & TTS-03 0
jf[ 7 - 1] = 1; # 7 Ne - Tops: f10.7<188 f10.7 unlimited 1
jf[ 8 - 1] = 1; # 8 foF2 from model foF2 or NmF2 - user input 1
jf[ 9 - 1] = 1; # 9 hmF2 from model hmF2 or M3000F2 - user input 1
jf[10 - 1] = 1; # 10 Te - Standard Te - Using Te/Ne correlation 1
jf[11 - 1] = 1; # 11 Ne - Standard Profile Ne - Lay-function formalism 1
jf[12 - 1] = 1; # 12 Messages to unit 6 to meesages.text on unit 11 1
jf[13 - 1] = 1; # 13 foF1 from model foF1 or NmF1 - user input 1
jf[14 - 1] = 1; # 14 hmF1 from model hmF1 - user input (only Lay version)1
jf[15 - 1] = 1; # 15 foE from model foE or NmE - user input 1
jf[16 - 1] = 1; # 16 hmE from model hmE - user input 1
jf[17 - 1] = 1; # 17 Rz12 from file Rz12 - user input 1
jf[18 - 1] = 1; # 18 IGRF dip, magbr, modip old FIELDG using POGO68/10 for 1973 1
jf[19 - 1] = 1; # 19 F1 probability model critical solar zenith angle (old) 1
jf[20 - 1] = 1; # 20 standard F1 standard F1 plus L condition 1
jf[21 - 1] = 1; # 21 ion drift computed ion drift not computed 0
jf[22 - 1] = 0; # 22 ion densities in % ion densities in m-3 1
jf[23 - 1] = 0; # 23 Te_tops (Aeros,ISIS) Te_topside (TBT-2011) 0
jf[24 - 1] = 0; # 24 D-region: IRI-95 Special: 3 D-region models (FIRI) 1
jf[25 - 1] = 1; # 25 F107D from APF107.DAT F107D user input (oarr(41)) 1
jf[26 - 1] = 0; # 26 foF2 storm model no storm updating 1
jf[27 - 1] = 1; # 27 IG12 from file IG12 - user 1
jf[28 - 1] = 0; # 28 spread-F probability not computed 0
jf[29 - 1] = 0; # 29 IRI01-topside new options as def. by JF(30) 0
jf[30 - 1] = 0; # 30 IRI01-topside corr. NeQuick topside model 0
# (29,30) = (1,1) IRIold, (0,1) IRIcor, (0,0) NeQuick, (1,0) Gulyaeva
jf[31 - 1] = 1; # 31 B0,B1 ABT-2009 B0 Gulyaeva h0.5 1
jf[32 - 1] = 1; # 32 F10.7_81 from file PF10.7_81 - user input (oarr(46)) 1
jf[33 - 1] = 0; # 33 Auroral boundary model on/off true/false 0
jf[34 - 1] = 0; # 34 Messages on Messages off 1
jf[35 - 1] = 0; # 35 foE storm model no foE storm updating 0
# .. ....
# 50 ....
# ------------------------------------------------------------------

return jf
def Switches():
"""
IRI switches to turn on/off several options
"""

jf = np.ones(50,dtype=bool)

# i 1 0 standard version
# ------------------------------------------------------------------------
jf[ 1 - 1] = 1; # 1 Ne computed Ne not computed 1
jf[ 2 - 1] = 1; # 2 Te, Ti computed Te, Ti not computed 1
jf[ 3 - 1] = 1; # 3 Ne & Ni computed Ni not computed 1
jf[ 4 - 1] = 0; # 4 B0 - Table option B0 - other models jf(31) 0
jf[ 5 - 1] = 0; # 5 foF2 - CCIR foF2 - URSI 0
jf[ 6 - 1] = 0; # 6 Ni - DS-95 & DY-85 Ni - RBV-10 & TTS-03 0
jf[ 7 - 1] = 1; # 7 Ne - Tops: f10.7<188 f10.7 unlimited 1
jf[ 8 - 1] = 1; # 8 foF2 from model foF2 or NmF2 - user input 1
jf[ 9 - 1] = 1; # 9 hmF2 from model hmF2 or M3000F2 - user input 1
jf[10 - 1] = 1; # 10 Te - Standard Te - Using Te/Ne correlation 1
jf[11 - 1] = 1; # 11 Ne - Standard Profile Ne - Lay-function formalism 1
jf[12 - 1] = 1; # 12 Messages to unit 6 to meesages.text on unit 11 1
jf[13 - 1] = 1; # 13 foF1 from model foF1 or NmF1 - user input 1
jf[14 - 1] = 1; # 14 hmF1 from model hmF1 - user input (only Lay version)1
jf[15 - 1] = 1; # 15 foE from model foE or NmE - user input 1
jf[16 - 1] = 1; # 16 hmE from model hmE - user input 1
jf[17 - 1] = 1; # 17 Rz12 from file Rz12 - user input 1
jf[18 - 1] = 1; # 18 IGRF dip, magbr, modip old FIELDG using POGO68/10 for 1973 1
jf[19 - 1] = 1; # 19 F1 probability model critical solar zenith angle (old) 1
jf[20 - 1] = 1; # 20 standard F1 standard F1 plus L condition 1
jf[21 - 1] = 1; # 21 ion drift computed ion drift not computed 0
jf[22 - 1] = 0; # 22 ion densities in % ion densities in m-3 1
jf[23 - 1] = 0; # 23 Te_tops (Aeros,ISIS) Te_topside (TBT-2011) 0
jf[24 - 1] = 0; # 24 D-region: IRI-95 Special: 3 D-region models (FIRI) 1
jf[25 - 1] = 1; # 25 F107D from APF107.DAT F107D user input (oarr(41)) 1
jf[26 - 1] = 0; # 26 foF2 storm model no storm updating 1
jf[27 - 1] = 1; # 27 IG12 from file IG12 - user 1
jf[28 - 1] = 0; # 28 spread-F probability not computed 0
jf[29 - 1] = 0; # 29 IRI01-topside new options as def. by JF(30) 0
jf[30 - 1] = 0; # 30 IRI01-topside corr. NeQuick topside model 0
# (29,30) = (1,1) IRIold, (0,1) IRIcor, (0,0) NeQuick, (1,0) Gulyaeva
jf[31 - 1] = 1; # 31 B0,B1 ABT-2009 B0 Gulyaeva h0.5 1
jf[32 - 1] = 1; # 32 F10.7_81 from file PF10.7_81 - user input (oarr(46)) 1
jf[33 - 1] = 0; # 33 Auroral boundary model on/off true/false 0
jf[34 - 1] = 0; # 34 Messages on Messages off 1
jf[35 - 1] = 0; # 35 foE storm model no foE storm updating 0
# .. ....
# 50 ....
# ------------------------------------------------------------------

return jf

#%%

def IRI(self, time, altkm, glat, glon, ap=None, f107=None, ssn=None, var=None):
def IRI( time, altkm, glat, glon, ap=None, f107=None, ssn=None, var=None):

if isinstance(time, str):
time = parse(time)
if isinstance(time, str):
time = parse(time)

altkm = np.atleast_1d(altkm)
altkm = np.atleast_1d(altkm)

# doy = squeeze(TimeUtilities().CalcDOY(year, month, dom))

# IRI options
jf = self.Switches()

# additional "input parameters" (necessary to scale the empirical model results
# to measurements)
addinp = -np.ones(12)

#------------------------------------------------------------------------------
#
if time.year < 1958:

addinp[10 - 1] = ssn # RZ12 (This switches 'jf[17 - 1]' to '0' and
# uses correlation function to estimate IG12)

jf[25 - 1] = 1 # 25 F107D from APF107.DAT F107D user input (oarr(41)) 1
jf[27 - 1] = 1 # 27 IG12 from file IG12 - user 1
jf[32 - 1] = 1 # 32 F10.7_81 from file PF10.7_81 - user input (oarr(46)) 1
# IRI options
jf = Switches()

else: # case for solar and geomagnetic indices from files
# additional "input parameters" (necessary to scale the empirical model results
# to measurements)
# addinp = -np.ones(12)

jf[17 - 1] = 1 # 17 Rz12 from file Rz12 - user input 1
jf[26 - 1] = 1 # 26 foF2 storm model no storm updating 1
jf[35 - 1] = 1 # 35 foE storm model no foE storm updating 0
#
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#
# if time.year < 1958:
#
# addinp[10 - 1] = ssn # RZ12 (This switches 'jf[17 - 1]' to '0' and
# # uses correlation function to estimate IG12)
#
# jf[25 - 1] = 1 # 25 F107D from APF107.DAT F107D user input (oarr(41)) 1
# jf[27 - 1] = 1 # 27 IG12 from file IG12 - user 1
# jf[32 - 1] = 1 # 32 F10.7_81 from file PF10.7_81 - user input (oarr(46)) 1
#
# else: # case for solar and geomagnetic indices from files
#
# jf[17 - 1] = 1 # 17 Rz12 from file Rz12 - user input 1
# jf[26 - 1] = 1 # 26 foF2 storm model no storm updating 1
# jf[35 - 1] = 1 # 35 foE storm model no foE storm updating 0
# #
# #------------------------------------------------------------------------------

mmdd = 100*time.month + time.day # month and dom (MMDD)
mmdd = 100*time.month + time.day # month and dom (MMDD)

# %% more inputs
jmag = 0 # 0: geographic; 1: geomagnetic
# iut = 0 # 0: for LT; 1: for UT
# height = 300. # in km
# h_tec_max = 2000 # 0: no TEC; otherwise: upper boundary for integral
# ivar = var # 1: altitude; 2: latitude; 3: longitude; ...
jmag = 0 # 0: geographic; 1: geomagnetic
# iut = 0 # 0: for LT; 1: for UT
# height = 300. # in km
# h_tec_max = 2000 # 0: no TEC; otherwise: upper boundary for integral
# ivar = var # 1: altitude; 2: latitude; 3: longitude; ...

# Ionosphere (IRI)
# Ionosphere (IRI)
# a, b = iriwebg(jmag, jf, glat, glon, int(time.year), mmdd, iut, time.hour,
# height, h_tec_max, ivar, ivbeg, ivend, ivstp, addinp, self.iriDataFolder)

# hour + 25 denotes UTC time
# hour + 25 denotes UTC time

outf,oarr = iri2016.iri_sub(jf, jmag, glat, glon,
time.year, mmdd, time.hour+25, altkm,
proot/'data/')
outf,oarr = iri2016.iri_sub(jf, jmag, glat, glon,
time.year, mmdd, time.hour+25, altkm,
proot/'data/')

# %% collect output
iri = xarray.DataArray(outf[:9,:].T,
coords={'alt_km':altkm, 'sim':simout},
dims=['alt_km','sim'],
attrs={'f107':oarr[40], 'ap':oarr[50],
'glat':glat,'glon':glon,'time':time,
'NmF2':oarr[0], 'hmF2':oarr[1],
'B0':oarr[9]
})
iri = xarray.DataArray(outf[:9,:].T,
coords={'alt_km':altkm, 'sim':simout},
dims=['alt_km','sim'],
attrs={'f107':oarr[40], 'ap':oarr[50],
'glat':glat,'glon':glon,'time':time,
'NmF2':oarr[0], 'hmF2':oarr[1],
'B0':oarr[9]
})

# FIRI Ne (in m-3)
# iri_ne_firi = self._RmNeg(a[13 - 1, :])[0]
Expand All @@ -152,7 +158,27 @@ def IRI(self, time, altkm, glat, glon, ap=None, f107=None, ssn=None, var=None):
# iriadd = { 'NmF2' : b[1 - 1, :][0], 'hmF2' : b[2 - 1, :][0],
# 'B0' : b[10 - 1, :][0] }

return iri
return iri


def timeprofile(tlim:tuple, dt:timedelta,
altkm:np.ndarray, glat:float, glon:float) -> xarray.DataArray:
"""compute IRI90 at a single altiude, over time range"""

T = datetimerange(tlim[0], tlim[1], dt)

iono = xarray.DataArray(np.empty((len(T),altkm.size,9)),
coords={'time':T,'alt_km':altkm, 'sim':simout},
dims=['time','alt_km','sim'],
)

for t in T:
iri = IRI(t, altkm, glat, glon)
iono.loc[t,...] = iri

iono.attrs = iri.attrs

return iono


class IRI2016Profile(IRI2016):
Expand Down Expand Up @@ -199,25 +225,14 @@ def __init__(self, time, alt_km, altlim=[90.,150.], altstp=2., htecmax=0,
self.jmag = jmag
self.lat = lat
self.lon = lon
self.month = time.month
self.dom = time.day
self.mmdd = self.month * 100 + self.dom
self.year = time.year
self.iut = iut
self.hour = time.hour
self.alt = alt_km

self.verbose = verbose
self.numstp = int((self.vend - self.vbeg) / self.vstp) + 1

if option == 'vertical':
self.HeiProfile()
elif option == 'lat':
self.LatProfile()
elif option == 'lon':
self.LonProfile()
elif option == 'time':
self.HrProfile()



# def _CallIRI(self):
Expand Down Expand Up @@ -284,53 +299,3 @@ def HeiProfile(self):
print('%8.3f %10.3f %8.3f %8.3f %8.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %7.3f %7.3f %7.3f %7.3f' %
(varval, edens, edratio, a[2,i], a[3,i], a[4,i], a[5,i], a[10,i], a[6,i],
a[7,i], a[8,i], a[9,i], b[32,i], b[38,i], b[40,i], b[45,i], b[50,i], b[51,i]))


def LatProfile(self):

self._CallIRI()

self._GetTitle()

if self.verbose:

lat = np.arange(self.vbeg, self.numstp*self.vstp, self.vstp)

print('\tGLON\tGLAT\tNmF2\t\thmF2\tB0')
for j in range(lat.size):
print('%8.3f %8.3f %8.3e %8.3f %8.3f' % (self.lon, lat[j], self.b[0, j], self.b[1, j], self.b[9, j]))


def LonProfile(self):

self._CallIRI()

self._GetTitle()

if self.verbose:

lon = np.arange(self.vbeg, self.numstp*self.vstp, self.vstp)

print('\tGLON\tGLAT\tNmF2\t\thmF2\tB0')
for j in range(lon.size):
print('%8.3f %8.3f %8.3e %8.3f %8.3f' % (lon[j], self.lat, self.b[0, j], self.b[1, j], self.b[9, j]))


def HrProfile(self):

self._CallIRI()

self._GetTitle()

t = np.arange(self.vbeg,self.numstp*self.vstp,self.vstp)

self.out = xarray.DataArray(self.a[:9,:t.size].T,
coords={'time':t,
'sim':['ne','Tn','Ti','Te','nO+','nH+','nHe+','nO2+','nNO+']},
dims=['time','sim'])

if self.verbose:

print(' GLON GLAT\tHR\tNmF2\thmF2\tB0')
for j in range(t.size):
print('%8.3f %8.3f %8.3f %8.3e %8.3f %8.3f' % (self.lon, self.lat, t[j], self.b[0, j], self.b[1, j], self.b[9, j]))

0 comments on commit e195488

Please sign in to comment.