Merge pull request #38 from astrojarred/astrojarred/issue37

Make `followup` module more robust (Addresses #37)

README.md

@@ -27,6 +27,8 @@
 - [Authors](#️authors)
 - [Instructions](#instructions)
   - [GW Observations](#instructions-gw-obs)
+    - [Example Simulation](#example-simulation)
+    - [Example followup calculation](#example-followup)
   - [Reading Results](#instructions-reading)
   - [Plotting Heatmaps](#instructions-plotting)
 
@@ -91,7 +93,7 @@ This code simulates observations of simulated gravitational wave events to deter
 
 - a dictionary object containing detection information and parameters of the event itself (extracted from the model)
 
-#### Example
+#### Example Simulation<a name ="example-simulation"></a>
 
 ```python
 from astropy import units as u
@@ -146,6 +148,54 @@ print(f"Observation time at delay={delay_time} is {res_ebl['obs_time']} with EBL
 # Obs time at delay=1800.0 s is 1292.0 s with EBL=franceschini
 ```
 
+### Example followup calculation<a name ="example-followup"></a>
+
+```python
+import astropy.units as u
+import pandas as pd
+from gravitational_wave_toy import followup
+
+lookup_talbe = "./O5_gammapy_observations_v4.parquet"
+
+# optional, but it's recommended to load the DataFrame first save time
+# otherwise you can directly pass the filepath to the get_exposure method
+lookup_df = pd.read_parquet(lookup_talbe)
+
+event_id = 1
+delay = 10 * u.s
+site = "north"
+zenith = 60
+ebl = "franceschini"
+
+
+followup.get_exposure(
+    event_id=event_id,
+    delay=delay,
+    site=site,
+    zenith=zenith,
+    extrapolation_df=lookup_df,
+    ebl=ebl,
+)
+
+# returns, e.g.
+# {
+#     'long': <Quantity 2.623 rad>,
+#     'lat': <Quantity 0.186 rad>,
+#     'eiso': <Quantity 2.67e+50 erg>,
+#     'dist': <Quantity 466000. kpc>,
+#     'obs_time': <Quantity 169. s>,
+#     'error_message': '',
+#     'angle': <Quantity 24.521 deg>,
+#     'ebl_model': 'franceschini',
+#     'min_energy': <Quantity 0.02 TeV>,
+#     'max_energy': <Quantity 10. TeV>,
+#     'seen': True,
+#     'start_time': <Quantity 10. s>,
+#     'end_time': <Quantity 179. s>,
+#     'id': 4
+# }
+```
+
 ### Reading Results<a name = "instructions-reading"></a>
 
 Note: The most recent simulations for CTA for the O5 observing run are [available on the CTA XWiki](https://cta.cloud.xwiki.com/xwiki/wiki/phys/view/Transients%20WG/Chasing%20the%20counterpart%20of%20GW%20alerts%20with%20CTA/O5%20Observation%20times%20with%20gw-toy%20package/) for CTA members.

gravitational_wave_toy/followup.py

@@ -66,14 +66,22 @@ def extrapolate_obs_time(
     if delay < min(event_dict.keys()):
         res["error_message"] = f"Minimum delay is {min(event_dict.keys())} seconds for this simulation"
         res["obs_time"] = -1
+        raise ValueError(f"Minimum delay is {min(event_dict.keys())} seconds for this simulation [{delay}s requested]")
+    elif delay > max(event_dict.keys()):
+        print(f"Warning: delay is greater than maximum delay of {max(event_dict.keys())}s for this simulation [{delay}s requested], value will be extrapolated.")
 
     # remove negative values
     pos_event_dict = {k: v for k, v in event_dict.items() if v > 0}
+
+    if not pos_event_dict:
+        res["error_message"] = f"Event is never detectable under the observation conditions {filters}"
+        res["obs_time"] = -1
+        return res
 
     # perform log interpolation
     log_event_dict = {log10(k): log10(v) for k, v in pos_event_dict.items()}
 
-    interp = interp1d(list(log_event_dict.keys()), list(log_event_dict.values()), kind="linear")
+    interp = interp1d(list(log_event_dict.keys()), list(log_event_dict.values()), kind="linear", fill_value="extrapolate")
 
     try:
         res["obs_time"] = 10**interp(log10(delay))
@@ -187,12 +195,22 @@ def get_exposure(
 
         obs_time = obs_info["obs_time"]
         if obs_time > 0:
-            obs_time = round(obs_time / target_precision.value) * target_precision
+            if obs_time > max_time.value:
+                obs_info["error_message"] = f"Exposure time of {int(obs_time)} s exceeds maximum time"
+                obs_time = -1
+            else:
+                obs_time = round(obs_time / target_precision.value) * target_precision
 
         # rename key
-        obs_info["angle"] = obs_info.pop("theta_view")
+        obs_info["angle"] = obs_info.pop("theta_view") * u.deg
         obs_info["ebl_model"] = obs_info.pop("irf_ebl_model")
 
+        # add other units
+        obs_info["long"] = obs_info["long"] * u.rad
+        obs_info["lat"] = obs_info["lat"] * u.rad
+        obs_info["eiso"] = obs_info["eiso"] * u.erg
+        obs_info["dist"] = obs_info["dist"] * u.kpc
+
         other_info = {
             "min_energy": min_energy,
             "max_energy": max_energy,