added NITE and MPF refraction models

.gitignore

@@ -1,4 +1,5 @@
 env
+set_simulator*py
 test_animation.py
 NITE*py
 test3.py

CHANGELOG.md

@@ -5,8 +5,9 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 
 
 ## 3.2.0
-Added NITE model.  Peng (2023), DOI: 10.1109/TGRS.2023.3332422
-Internally the NITE model is known as refraction model 5.
+Added NITE model.  In gnssrefl this is refl_model 5.
+
+For details see Peng (2023), DOI: 10.1109/TGRS.2023.3332422
 While the NITE publication advocates using estimated zenith troposphere delays,
 we provide no such access here, and use Saastamoinen model corrections using
 standard met outputs created by Generic Mapping functions.  If someone wants
@@ -17,6 +18,8 @@ using the optional apriori_rh. The NITE model should be much more accurate
 than the default refraction model (1). It is not clear how much more precise
 it is.  I would expect it to be more precise for tall sites. 
 
+Added MPF model - name may change.  In gnssrefl this is refr_model 6.
+
 ## 3.1.5
 
 quickplt has been significantly updated to allow plotting of SNR files.

gnssrefl/gnssir_v2.py

@@ -218,13 +218,12 @@ def gnssir_guts_v2(station,year,doy, snr_type, extension,lsp):
             # I do not understand why all these extra parameters are sent to this 
             # function as they are not used. Maybe I was doing some testing.
             ele=refr.Ulich_Bending_Angle(snrD[:,1],N_ant,lsp,p,T,snrD[:,3],snrD[:,0])
-        elif irefr == 5:
+        elif (irefr == 5 ) or (irefr == 6):
             # NITE MODEL
             # remove ele < 1.5 cause it blows up
             i=snrD[:,1] > 1.5
             snrD = snrD[i,:]
             N = len(snrD[:,1])
-            print('NITE refraction correction, Peng et al. remove data < 1.5 degrees')
             # elevation angles, degrees
             ele=snrD[:,1]
             # zenith angle in radians
@@ -234,11 +233,14 @@ def gnssir_guts_v2(station,year,doy, snr_type, extension,lsp):
             [mpf, dmpf]=[refr.mpf_tot(gmf_h,gmf_w,zhd,zwd),refr.mpf_tot(dgmf_h,dgmf_w,zhd,zwd)]
 
             # inputs apriori RH, elevation angle, refractive index, zenith delay, mpf ?, dmpf?
-            dE_NITE=refr.Equivilent_Angle_Corr_NITE(RH_apriori, ele, N_ant, ztd, mpf, dmpf)
-            ele = ele + dE_NITE 
-        elif irefr == 6:
-            print('MPF refraction correction, Wiliams and Nievinski? Exiting')
-            return
+            if irefr == 5:
+                print('NITE refraction correction, Peng et al. remove data < 1.5 degrees')
+                dE_NITE=refr.Equivalent_Angle_Corr_NITE(RH_apriori, ele, N_ant, ztd, mpf, dmpf)
+                ele = ele + dE_NITE 
+            else: 
+                print('MPF refraction correction, Wiliams and Nievinski')
+                dE_MPF= refr.Equivalent_Angle_Corr_mpf(ele,mpf,N_ant,RH_apriori)
+                ele = ele + dE_MPF 
 
         elif irefr == 0:
             print('No refraction correction applied ')
@@ -465,15 +467,13 @@ def set_refraction_params(station, dmjd,lsp):
         irefr = refraction_model
         refr.readWrite_gpt2_1w(xdir, station, lsp['lat'], lsp['lon'])
         # time varying was originally set to no for now (it = 1)
-        # now allows time varying for models 2 and 4
-        if refraction_model in [2, 4, 5]:
-        #if (refraction_model == 2) or (refraction_model == 4):
+        # now allows time varying for models 2, 4 and now MPF and NITE
+        if refraction_model in [2, 4, 5, 6]:
             it = 0
             print('Using time-varying refraction parameter inputs')
         #elif (refraction_model == 5):
         else:
             it = 1
-        # I believe this uses the wrong height
         dlat = lsp['lat']*np.pi/180; dlong = lsp['lon']*np.pi/180; ht = lsp['ht']
         p,T,dT,Tm,e,ah,aw,la,undu = refr.gpt2_1w(station, dmjd,dlat,dlong,ht,it)
         # e is water vapor pressure, so humidity ??

gnssrefl/gps.py

@@ -5099,8 +5099,10 @@ def rapid_gfz_orbits(year,month,day):
         #print(littlename, ' already exists on disk')
         return littlename, fdir, True 
     try:
+        #g.replace_wget(url,littlename + '.gz')
         wget.download(url,littlename + '.gz')
-        subprocess.call(['gunzip', littlename + '.gz'])
+        if os.path.isfile(littlename + '.gz'):
+            subprocess.call(['gunzip', littlename + '.gz'])
     except:
         print('Problems downloading Rapid GFZ orbit')
         print(url)

gnssrefl/quickplt.py

@@ -27,7 +27,7 @@ def parse_arguments():
     parser.add_argument("-ymdhm", help="True/T for when columns 1-5 are year month day hour minute", type=str,default=None)
     parser.add_argument("-xlabel", type=str, help="optional x-axis label", default=None)
     parser.add_argument("-ylabel", type=str, help="optional y-axis label", default=None)
-    parser.add_argument("-symbol", help="plot symbol, e.g. * or -", type=str,default=None)
+    parser.add_argument("-symbol", help="plot symbol, e.g. o or -", type=str,default=None)
     parser.add_argument("-title", help="optional title", type=str,default=None)
     parser.add_argument("-outfile", help="optional filename for plot. Must end in png", type=str,default=None)
     parser.add_argument("-xlimits", nargs="*",type=float, help="optional x-axis limits", default=None)
@@ -416,7 +416,7 @@ def run_quickplt (filename: str, xcol: str, ycol: str, errorcol: int=None, mjd:
     # i.e. using default
     if symbol is None:
         if yerrors:
-            ax.errorbar(tval, yval, yerr=tvd[:,errorcol], fmt='.',color='blue')
+            ax.errorbar(tval, yval, yerr=tvd[:,errorcol], fmt='.')
         else:
             ax.plot(tval, yval, 'b.')
     else:
@@ -430,7 +430,8 @@ def run_quickplt (filename: str, xcol: str, ycol: str, errorcol: int=None, mjd:
         if symbol is None:
             ax.plot(tval2, yval2, 'r.')
         else:
-            ax.plot(tval2, yval2, color='red', fmt=symbol)
+            # ??
+            ax.plot(tval2, yval2, symbol )
 
     myplt.grid()
     myplt.ylabel(ylabel)

gnssrefl/refraction.py

@@ -641,7 +641,7 @@ def refrc_Rueger(drypress,vpress,temp):
 
     """
 
-    # Rüeger's "best average", for 375 ppm CO2, 77.690 for 392 ppm CO2
+    # Rueger's "best average", for 375 ppm CO2, 77.690 for 392 ppm CO2
     [K1r,K2r,K3r]=[77.689 ,71.2952 ,375463.]     
 
     # Rueger,IAG recommend, in ppm
@@ -656,26 +656,29 @@ def refrc_Rueger(drypress,vpress,temp):
     vdensity=(vpress*100.)*Mw/(Rgas*temp)   
     totaldensity=drydensity+vdensity
 
-    #divide 100 because hPa in Rueger formula, result in ppm
+    #divide by 100 because hPa in Rueger formula, result in ppm
     Nhydro=(K1r*Rgas*totaldensity/Md)/100.   
     K2rr=K2r-K1r*(Mw/Md)
     Nwet=K2rr * vpress / temp + K3r * vpress / (temp ** 2)
     # in ppm
-    ref=[round(Nrueger,4),round(Nhydro,4),round(Nwet,4)]    
+    ref=[Nrueger,Nhydro,Nwet]    
+    #ref=[round(Nrueger,4),round(Nhydro,4),round(Nwet,4)]    
 
     return ref
 
 
-def Equivilent_Angle_Corr_NITE(Hr_apr, e_T, N_ant, ztd_ant, mpf_tot, dmpf_de_tot):
+def Equivalent_Angle_Corr_NITE(Hr_apr, e_T, N_ant, ztd_ant, mpf_tot, dmpf_de_tot):
     """
 
     This function computes the "equvilent" angular correction to apply the 
     NITE formula on the true elevation angle ele_eqv = e_T + dele
 
 
     Equation (24) in Peng (2023), DOI: 10.1109/TGRS.2023.3332422
-    The variable substitude method can be found in Strandberg, J. (2020). New
-    methods and applications for interferometric GNSS reflectometry. Chalmers Tekniska Hogskola (Sweden).
+
+    The variable substitude method can be found in Strandberg, J. (2020). 
+    New methods and applications for interferometric GNSS reflectometry. 
+    Chalmers Tekniska Hogskola (Sweden).
 
     Parameters
     ----------
@@ -738,7 +741,7 @@ def gmf_deriv(dmjd,dlat,dlon,dhgt,zd):
         height in meters
     zd: float
         zenith distance in radians ??? ( is this really what you mean??
-        I suspect it is the zenith angle ... in radians
+        KL: I suspect it is the zenith angle ... in radians
 
     Returns
     -------
@@ -1112,45 +1115,6 @@ def N_layer(N_antenna, Hr):
     #return round(Nl, 4)     #in ppm
     return Nl      #in ppm
 
-def saastam_wet(p,T):
-    """
-    This was never finished
-
-    Parameters
-    ----------
-    p : float
-        total pressure, hPa
-    T : float
-        temperature in Kelvin
-    e : float
-        partial pressure of water vapor
-    h_rel : float
-        relative humidity
-    h : float
-        geodetic height above sea level (meters)
-
-    Returns
-    -------
-    WZD : float
-        wet zenith delay, meters
-
-    """
-    #TM : mean temperature of the water vapor
-    # got these values from UNB package
-    TM = T * (1 - BETA * RD/DEN)
-    MD = 28.9644
-    MW = 18.0152;
-    K1 = 77.604;
-    K2 = 64.79;
-    K3 = 3.776e5;
-    R  = 8314.34;
-    C1 = 2.2768e-03;
-    K2PRIM = K2 - K1*(MW/MD);
-    RD     = R / MD;
-
-    WZD = 1.0E-6 * (K2PRIM + K3/TM) * RD * E/DEN
-
-    return WZD
 
 def saastam2(press, lat, height):
     """
@@ -1171,9 +1135,9 @@ def saastam2(press, lat, height):
     press : float
         atmospheric total pressure in hPa
     lat : float
-        latitude of the station in degree
+        latitude of the station, degrees
     height : float
-        height of the station in meters ??? Which kind of height ????
+        ellipsoidal height of the station in meters 
 
     Returns
     -------
@@ -1192,7 +1156,9 @@ def saastam2(press, lat, height):
 def mpf_tot(gmf_h, gmf_w, zhd, zwd):  
     """
     Finds the total mapping function by weighting the hydrostatic and wet mapping 
-    function with the zenith hydrostatic and wet delay. Author: Peng
+    function with the zenith hydrostatic and wet delay. 
+
+    Author: Peng Feng
     
     Parameters
     ----------
@@ -1225,7 +1191,7 @@ def dmpf_dh(ele, dhgt):
 
     Boehm, J., Werl, B., & Schuh, H. (2006). Troposphere mapping functions for
     GPS and very long baseline interferometry from European Centre for Medium‐Range
-    Weather Forecasts operational analysis data. Journal of geophysical research: solid earth, 111(B2).
+    Weather Forecasts operational analysis data. JGR: Solid Earth, 111(B2).
 
     Parameters
     ----------
@@ -1279,7 +1245,7 @@ def Ulich_Bending_Angle_original(ele, N0):
     return np.degrees(r * f)
 
 
-def Equivilent_Angle_Corr_mpf(ele, mpf_tot, N0, Hr_apr):
+def Equivalent_Angle_Corr_mpf(ele, mpf_tot, N0, Hr_apr):
     """
     This function computes the "equvilent" angular correction to apply the 
     tropospheric delay calculated with the mapping function.
@@ -1318,7 +1284,7 @@ def Equivilent_Angle_Corr_mpf(ele, mpf_tot, N0, Hr_apr):
 
 def asknewet(e, Tm, lambda_val):
     """
-    This function determines the zenith wet delay based on the equation 22 by Askne and Nordius (1987)
+    Determines the zenith wet delay based on the equation 22 by Askne and Nordius (1987)
 
     Askne and Nordius, Estimation of tropospheric delay for microwaves from surface weather data,
     Radio Science, Vol 22(3): 379-386, 1987.