From 948a1e6829f1ba5548b898e47a41865d995f1b37 Mon Sep 17 00:00:00 2001
From: Aron Mirwald <aron.mirwald@sed.ethz.ch>
Date: Tue, 2 Jul 2024 20:26:14 +0200
Subject: [PATCH] adjusting the a-value estimation

---
 seismostats/analysis/estimate_a.py | 41 +++++++++++-------------------
 1 file changed, 15 insertions(+), 26 deletions(-)

diff --git a/seismostats/analysis/estimate_a.py b/seismostats/analysis/estimate_a.py
index 22e285e..ad2a84e 100644
--- a/seismostats/analysis/estimate_a.py
+++ b/seismostats/analysis/estimate_a.py
@@ -8,42 +8,33 @@
 
 def estimate_a(magnitudes: np.ndarray,
                mc: float | None = None,
-               m0: float | None = None,
+               m_ref: float | None = None,
                b_value: float | None = None,
-               T: dt.timedelta | None = None,
-               delta_T: dt.timedelta | None = None,
+               T: float | None = None,
                ) -> float:
     """Estimate the a-value of the Gutenberg-Richter (GR) law.
 
-    N = 10 ** (a - b * (m - mc)) (1)
+    N = 10 ** (a - b * (m_ref - mc)) (1)
 
-    where N is the number of events with magnitude greater than m0, T is the
-    time interval the magnitudes are observed, delta_T is the time interval
-    for the estimation of the a-value.
-    In order to make comparability across different time intervals possible,
-    the a-value is often estimated at a reference magnitude m0 and a reference
-    time interval delta_T. Then, instead of (1), the following equation is used:
-
-    N T/ delta_T = 10 ** (a - b * (m0 - mc)) (2)
+    where N is the number of events with magnitude greater than m_ref, which
+    occurred in the time interval T. T should be given as a float- to be 
+    precise, it should be the time interval scaled to the time-unit of interest.
+    E.g., if the number of events per year are of interest, T should be the
+    number of years in which the events occurred.
 
     If only the magnitudes are given, the a-value is estimated at the lowest
     magnitude in the sample, with mc = min(magnitudes). Eq. (1) then simplifies
-    to N = 10**a. Otherwise, the a-value is estimated at the reference
-    magnitude m0 and the reference time interval delta_T using
-    eq. (2).
+    to N = 10**a.
 
     Args:
         magnitudes: Magnitude sample
         mc:         Completeness magnitude. If None, the lowest magnitude is
                 used as completeness magnitude.
-        m0:         Reference magnitude for which the a-value is estimated. If
+        m_ref:      Reference magnitude for which the a-value is estimated. If
                 None, the a-value is estimated at mc.
         b_value:    b-value of the Gutenberg-Richter distribution
-        T:          Reference time interval for the estimation of the a-value.
-                If None, T is set to the whole interval of the times. Often, T
-                is set to 1 year.
-        delta_T:    Time interval which the magnitudes are observed in. Only
-                needed if T is not None.
+        T:          Relative length of the time interval in which the events
+                occurred (relative to the time unit of interest, e.g., years)
 
     Returns:
         a: a-value of the Gutenberg-Richter distribution
@@ -59,16 +50,14 @@ def estimate_a(magnitudes: np.ndarray,
     a = np.log10(len(magnitudes))
 
     # scale to reference magnitude
-    if m0 is not None:
+    if m_ref is not None:
         if b_value is None:
             raise ValueError(
                 "b_value must be provided if m0 is given")
-        a = a - b_value * (m0 - mc)
+        a = a - b_value * (m_ref - mc)
 
     # scale to reference time-interval
     if T is not None:
-        if delta_T is None:
-            raise ValueError("T must be provided if delta_T is given.")
-        a = a + np.log10(T / delta_T)
+        a = a - np.log10(T)
 
     return a