From d01fc7d1186cf7095225fc5a8fdac6fe6afd7f09 Mon Sep 17 00:00:00 2001
From: Aron Mirwald <aron.mirwald@sed.ethz.ch>
Date: Fri, 5 Jul 2024 16:43:38 +0200
Subject: [PATCH] fix problems and add warnings to estimate_a_pos

---
 seismostats/analysis/__init__.py   |  1 +
 seismostats/analysis/estimate_a.py | 10 +++++++---
 2 files changed, 8 insertions(+), 3 deletions(-)

diff --git a/seismostats/analysis/__init__.py b/seismostats/analysis/__init__.py
index 553b6a1..054794e 100644
--- a/seismostats/analysis/__init__.py
+++ b/seismostats/analysis/__init__.py
@@ -13,3 +13,4 @@
                                                 b_value_to_beta
                                                 )
 from seismostats.analysis.estimate_mc import mc_ks, mc_max_curvature
+from seismostats.analysis.estimate_a import estimate_a_positive
diff --git a/seismostats/analysis/estimate_a.py b/seismostats/analysis/estimate_a.py
index 0db2322..6865e61 100644
--- a/seismostats/analysis/estimate_a.py
+++ b/seismostats/analysis/estimate_a.py
@@ -108,7 +108,9 @@ def estimate_a_positive(
             warnings.warn(
                 "Completeness magnitude is higher than the lowest magnitude."
                 "Cutting the magnitudes to the completeness magnitude.")
-        magnitudes = magnitudes[magnitudes >= mc]
+        idx = magnitudes >= mc
+        magnitudes = magnitudes[idx]
+        times = times[idx]
 
     if dmc is None:
         dmc = delta_m
@@ -120,7 +122,7 @@ def estimate_a_positive(
     # differences
     mag_diffs = np.diff(magnitudes)
     time_diffs = np.diff(times)
-    if not np.all(time_diffs >= 0):
+    if not np.all(time_diffs >= 0 * time_diffs[0]):
         raise ValueError("Times are not ordered correctly.")
 
     # only consider events with magnitude difference >= dmc
@@ -131,7 +133,7 @@ def estimate_a_positive(
     # estimate the number of events within the time interval
     total_time = times[-1] - times[0] + np.mean(np.diff(times))
     total_time_pos = sum(time_diffs) + np.mean(time_diffs)
-    n_pos = len(mag_diffs) / total_time_pos * total_time
+    n_pos = total_time / total_time_pos * len(mag_diffs)
 
     # estimate a-value
     a = np.log10(n_pos)
@@ -146,3 +148,5 @@ def estimate_a_positive(
     # scale to reference time-interal or volume of interest
     if scaling is not None:
         a = a - np.log10(scaling)
+
+    return a