Add Dunai 2000: beta and inclination

corny/corn.py

@@ -265,6 +265,9 @@ def C_press(p, pref=1013.15, beta=136, method='Zreda et al. (2012)', inclination
     if method == 'Zreda et al. (2012)':
         r = np.exp((p - pref) / beta)
     elif method == 'Dunai et al. (2000)':
+        # Todo: 
+        # from .grains.Dunai2000 import get_beta
+        # beta = get_beta(lon, lat, alt, date)
         beta = 129.55 + 19.85 / (1+np.exp((inclination-62.05)/5.43))**3.59
         r = np.exp((p - pref) / beta)
     else:

corny/grains/Dunai2000.py

+import numpy as np
+
+def get_beta(lon, lat, alt, t):
+    """
+    Takes coordinates and date 
+    Returns barometric neutron attenuation coefficient, beta 
+    Example:
+        beta = get_beta(12.42, 51.58, 100, datetime(2014,7,1))
+    """
+    i = get_inclination(lon, lat, alt, t)
+    L = Dunai_2000(i)
+    beta = 1/L
+    return(beta)
+
+def Dunai_2000(inclination=0):
+    """
+    Takes inclination of the magnetic field 
+    Returns neutron attenuation length (L=1/beta)
+    Example:
+        data = pandas.DataFrame()
+        data["inclination"] = np.arange(0,90)
+        data["L"] = Dunai_2000(data["inclination"])
+    """
+    y = 129.55
+    x = 62.05
+    a = 19.85
+    b = -5.43
+    c = 3.59
+    L = y+a/(1+np.exp((x-inclination)/b))**c
+    return(L)
+
+# %%
+
+def get_inclination(lon, lat, alt, t):
+    """
+    Takes coordinates and date 
+    Returns inclination of the Geomagnetic field
+    Example:
+        incl = get_inclination(12.42, 51.58, 100, datetime(2014,7,1))
+    """
+    import ppigrf
+    Be, Bn, Bu = ppigrf.igrf(lon, lat, alt, t)
+    inclination = np.arctan(-Bu/np.sqrt(Be**2 + Bn**2))
+    inclination_deg = np.degrees(inclination)
+    return(float(inclination_deg))
+