Merge pull request #978 from linea-it/973-account-for-stars-proper-mo…

…tion-in-praia-star-search-task

Optimizes star search radius in occultation

predict_occultation/src/predict_occultation/pipeline/praia_occ_input_template.txt

@@ -7,7 +7,7 @@ DE435/JPL                                         | Ephemeris label
 {stars_parameters_of_occultation}| output file : star parameters of occultation (minimum geocentric distance,t_occ,t_initial,t_final, etc)
 {stars_parameters_of_occultation_plot}| output file : star parameters of occultation (Bruno Sicardy data plot format)
 {projected_search_circle}| radius (arcsec) of projected search circle (~ projected Earth figure + occultation body)
-01.0000                                           | minimum outreach radius (arcsec) for fast elimination of farway ephemeris points
+{minimum_outreach_radius}| minimum outreach radius (arcsec) for fast elimination of farway ephemeris points
 12d0  12d0                                        | exclusion range of Local Solar Time (day light): min max (hours)
 2005                                              | to (years - see below)
 +00.00000d0                                       | bofra

predict_occultation/src/predict_occultation/pipeline/run_praia_occ.py

@@ -184,12 +184,25 @@ def start_praia_occ(
 
     print("Gaia CSV: [%s]" % gaia_csv)
 
+    # Quando o diametro do objeto exitir no json, ele é passado para a função search_candidates
+    # que cria o arquivo praia_occ_star_search_12.dat. Sua função é reduzir o numero de calculo necessário
+    # especialmente para objetos de diametros pequenos.
+    object_diameter_upper_limit = obj_data.get("diameter_err_max", None)
+    object_diameter = obj_data.get("diameter", None)
+    if object_diameter_upper_limit is None:
+        if object_diameter is not None:
+            object_diameter *= 1.2
+    else:
+        object_diameter += object_diameter_upper_limit
+    print("Object Diameter: [%s]" % object_diameter)
+
     # Run PRAIA OCC Star Search 12
     # Criar arquivo .dat baseado no template.
     occultation_file = search_candidates(
         star_catalog=gaia_cat,
         object_ephemeris=eph_file,
         filename=occultation_table_filename,
+        object_diameter=object_diameter,
     )
 
     print("Occultation CSV Table: [%s]" % occultation_file)

predict_occultation/src/predict_occultation/pipeline/search_candidates.py

@@ -9,7 +9,7 @@
 from library import HMS2deg
 
 
-def get_best_projected_search_circle(object_ephemeris):
+def get_best_projected_search_radius(object_ephemeris, object_diameter):
     """
     Calculate the best projected search circle based on object ephemeris.
 
@@ -24,21 +24,30 @@ def get_best_projected_search_circle(object_ephemeris):
         The function calculates the best distance from the ephemeris and computes the projected search circle size.
     """
     earth_radius = 6371  # km
-    body_radius_compensation = 1000  # km
+    body_radius_compensation = (
+        0 if object_diameter is None else object_diameter / 2
+    )  # km
     distances = []
     with open(object_ephemeris, "r") as file:
         for i, line in enumerate(file, start=3):
             distances.append(line[120:160])
     distances = np.array(distances[3:], dtype=float)
     best_distance = distances.min()
-    projected_search_circle = (
-        2 * (earth_radius + body_radius_compensation) / best_distance
+    projected_search_diameter = (
+        2 * (earth_radius + body_radius_compensation * 2) / best_distance
     )
-    projected_search_circle *= 3600 * 180 / np.pi  # converts to arcsec
-    return np.around(projected_search_circle, 4)
+    projected_search_diameter *= 3600 * 180 / np.pi  # converts to arcsec
+    projected_search_radius = projected_search_diameter / 2
+    return np.around(projected_search_radius, 4)
 
 
-def praia_occ_input_file(star_catalog, object_ephemeris):
+def get_minimum_outreach_radius(projected_search_radius):
+    # The minimum outreach radius is X % greater than the projected search radius.
+    # All stars outside this region will be discounted to accelarate computations
+    return projected_search_radius * 1
+
+
+def praia_occ_input_file(star_catalog, object_ephemeris, object_diameter):
 
     # TODO: Alguns dos parametros podem vir da interface.
     try:
@@ -80,9 +89,19 @@ def praia_occ_input_file(star_catalog, object_ephemeris):
             data = data.replace(
                 "{stars_parameters_of_occultation_plot}", name.ljust(50)
             )
+            projected_search_circle = get_best_projected_search_radius(
+                object_ephemeris, object_diameter
+            )
             data = data.replace(
                 "{projected_search_circle}",
-                str(get_best_projected_search_circle(object_ephemeris)).ljust(50),
+                f"{projected_search_circle:2.7f}".ljust(50),
+            )
+            minimum_outreach_radius = get_minimum_outreach_radius(
+                projected_search_circle
+            )
+            data = data.replace(
+                "{minimum_outreach_radius}",
+                f"{minimum_outreach_radius:2.7f}".ljust(50),
             )
 
             with open(output, "w") as new_file:
@@ -206,7 +225,7 @@ def ascii_to_csv(inputFile, outputFile):
     np.savetxt(outputFile, newData, fmt="%s", header=colNames, delimiter=";")
 
 
-def search_candidates(star_catalog, object_ephemeris, filename):
+def search_candidates(star_catalog, object_ephemeris, filename, object_diameter):
 
     try:
         app_path = os.environ.get("APP_PATH").rstrip("/")
@@ -225,7 +244,9 @@ def search_candidates(star_catalog, object_ephemeris, filename):
 
         # Criar arquivo .dat baseado no template.
         search_input = praia_occ_input_file(
-            star_catalog=star_catalog, object_ephemeris=object_ephemeris
+            star_catalog=star_catalog,
+            object_ephemeris=object_ephemeris,
+            object_diameter=object_diameter,
         )
 
         print("PRAIA OCC .DAT: [%s]" % search_input)