nans

Author: scivision

AltitudeProfile.py

@@ -15,11 +15,11 @@
     p.add_argument('--f107a',type=float,default=200.)
     p.add_argument('--ap',type=int,default=4)
     p = p.parse_args()
-
+# %% user parameters
     altkm = np.arange(p.alt[0], p.alt[1], p.alt[2])
     dtime = parse(p.time)
-#%%
+# %% run IRI90 across altitude
     iono =  pyiri90.runiri(dtime,altkm,p.latlon,p.f107, p.f107a, ap=p.ap)
-
+# %% plots
     pyiri90.plots.plotiono(iono,dtime,p.latlon,p.f107,p.f107a,p.ap)
     show()

TimeProfile.py

@@ -0,0 +1,27 @@
+#!/usr/bin/env python
+from dateutil.parser import parse
+from datetime import timedelta
+from matplotlib.pyplot import show
+#
+import pyiri90, pyiri90.plots
+
+if __name__ == '__main__':
+    from argparse import ArgumentParser
+    p = ArgumentParser()
+    p.add_argument('-t','--trange',help='START STOP STEP (hours) time [UTC]',nargs=3,
+                   default=('2017-08-01','2017-08-02',1.))
+    p.add_argument('--alt',help='altitude [km]',type=float, default=150.)
+    p.add_argument('-c','--latlon',help='geodetic coordinates of simulation',
+                   type=float,default=(65,-147.5))
+    p.add_argument('--f107',type=float,default=200.)
+    p.add_argument('--f107a', type=float,default=200.)
+    p.add_argument('--ap', type=int, default=4)
+    p = p.parse_args()
+# %% user parameters
+    tlim = (parse(p.trange[0]), parse(p.trange[1]))
+    dt = timedelta(hours=p.trange[2])
+# %% run IRI90 across time
+    iono =  pyiri90.timeprofile(tlim, dt, p.altkm, p.latlon, p.f107, p.f107a, ap=p.ap)
+# %% plots
+    pyiri90.plots.plottime(iono)
+    show()

fortran/iri90.f

@@ -160,13 +160,12 @@ SUBROUTINE IRI90(JF,JMAG,ALATI,ALONG,RZ12,MMDD,DHOUR,
 
       logical, intent(in) :: JF(12)
       integer,intent(in) :: JMAG,MMDD
-      real,intent(in) :: alati,along,rz12
+      real,intent(in) :: alati,along,rz12, dhour, zkm(nz)
       real,intent(out) :: outf(11,nz),oarr(30)
 
-      dimension zkm(nz)
-      character*(*) drect
-      character*50 path
-      character*10 filename
+      character(*) drect
+      character(50) path
+      character(10) filename
       INTEGER               DAYNR,DDO,DO2,SEASON,SEADAY
       REAL               LATI,LONGI,MO2,MO,MODIP,NMF2,MAGBR
       REAL                NMF1,NME,NMD,MM,MLAT,MLONG,NOBO2
@@ -432,7 +431,7 @@ SUBROUTINE IRI90(JF,JMAG,ALATI,ALONG,RZ12,MMDD,DHOUR,
 
 8448       write(monito,8449) path
 8449       format(' IRI90: File ',A50,'not found')
-       stop
+       error stop
 C
 C LINEAR INTERPOLATION IN SOLAR ACTIVITY
 C
@@ -2548,10 +2547,10 @@ SUBROUTINE LNGLSN ( N, A, B, AUS)
                      A(N,K) = ( B(K) - AMAX ) / A(K,K)
               ENDIF
 6       CONTINUE
-       RETURN
-       END
-C
-C
+
+       END SUBROUTINE LNGLSN
+
+
        SUBROUTINE LSKNM ( N, M, M0, M1, HM, SC, HX, W, X, Y, VAR, SING)
 C --------------------------------------------------------------------
 C   DETERMINES LAY-FUNCTIONS AMPLITUDES FOR A NUMBER OF CONSTRAINTS:
@@ -2999,10 +2998,10 @@ SUBROUTINE CIRA86(IDAY,SEC,GLAT,GLONG,STL,F107A,TINF,TLB,SIGMA)
       TLB = 386.0 * ( 1. + T1+T4+T5+T7+T8 ) * 0.976619
 C  Sigma      [Eq. A5]
       SIGMA = G0 / ( TINF - TLB )
-      RETURN
-      END
-C
-C
+
+      END SUBROUTINE CIRA86
+
+
       pure real FUNCTION TN(H,TINF,TLBD,S)
       implicit none
 
@@ -3016,8 +3015,8 @@ pure real FUNCTION TN(H,TINF,TLBD,S)
       TN = TINF - TLBD * EXP ( - S * ZG2 )
 
       END function TN
-C
-C
+
+
       pure real FUNCTION DTNDH(H,TLBD,S)
       implicit None
 

pyiri90/__init__.py

@@ -1,6 +1,7 @@
 """
 IRI-90 international reference ionosphere in Python
 """
+from datetime import datetime,timedelta
 from pathlib import Path
 import numpy as np
 from xarray import DataArray
@@ -10,23 +11,39 @@
 rdir = Path(__file__).parent
 Ts = ['Tn','Ti','Te']
 
-def runiri(dt,z,glatlon,f107,f107a,ap):
+
+def datetimerange(start:datetime, end:datetime, step:timedelta) -> list:
+    """like range() for datetime!"""
+    assert isinstance(start, datetime)
+    assert isinstance(end, datetime)
+    assert isinstance(step, timedelta)
+
+    return [start + i*step for i in range((end-start) // step)]
+
+
+def runiri(t, altkm:float, glatlon:tuple, f107:float, f107a:float, ap:int):
     glat,glon = glatlon
     jmag=0 # coordinates are: 0:geographic 1: magnetic
 
     JF = np.array((1,1,1,1,0,1,1,1,1,1,1,0),bool) #Solomon 1993 version of IRI
     #JF = (1,1,1) + (0,0,0) +(1,)*14 + (0,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,1,1) #for 2013 version of IRI
 
+    monthday = (t.strftime('%m%d'))
+    hourfrac = t.hour + t.minute//60+ t.second//3600
 #%% call IRI
-    outf,oarr = iri90.iri90(JF,jmag,glat,glon % 360., -f107,
-                       dt.strftime('%m%d'),
-                       dt.hour+dt.minute//60+dt.second//3600,
-                       z,str(rdir/'data')+'/')
+    outf,oarr = iri90.iri90(JF,jmag,
+                            glat,glon % 360.,
+                            -f107,
+                            monthday,
+                            hourfrac,
+                            altkm,
+                            str(rdir/'data')+'/')
 #%% arrange output
     iono = DataArray(outf[:9,:].T,
-                     coords={'alt_km':z,
+                     coords={'alt_km':altkm,
                              'sim':['ne','Tn','Ti','Te','nO+','nH+','nHe+','nO2+','nNO+']},
-                     dims=['alt_km','sim'])
+                     dims=['alt_km','sim'],
+                     attrs={'f107':f107, 'f107a':f107a, 'ap':ap, 'glatlon':glatlon})
 
 #    i=(iono['Ti']<iono['Tn']).values
 #    iono.ix[i,'Ti'] = iono.ix[i,'Tn']
@@ -41,7 +58,12 @@ def runiri(dt,z,glatlon,f107,f107a,ap):
     #iono['nN+'] = iono['ne'] * outf[10,:]/100.
 # %% negative indicates undefined
     for c in iono.sim:
-
         iono.loc[iono.loc[:,c]<= 0., c] = np.nan
 
     return iono
+
+def timeprofile(tlim:tuple, dt:timedelta, altkm:float, glatlon:tuple, f107:float, f107a:float, ap:int):
+    """compute IRI90 at a single altiude, over time range"""
+
+    t = datetimerange(tlim[0], tlim[1], dt)
+