Merge pull request #62 from rahil-makadia/dev

program code mapping fix; gaia ades bias nan-value fix

grss/__init__.py

@@ -4,7 +4,6 @@
 from . import utils
 
 # get version from version.txt
-with open(f'{utils.grss_path}/version.txt', 'r', encoding='utf-8') as f:
-    __version__ = f.read().strip()
+__version__ = open(f'{utils.grss_path}/version.txt', 'r', encoding='utf-8').read().strip()
 
 utils.initialize()

grss/fit/fit_optical.py

@@ -343,6 +343,8 @@ def add_gaia_obs(obs_df, t_min_tdb=None, t_max_tdb=None, gaia_dr='gaiadr3', verb
         obs_df.loc[idx, 'sigRA'] = ra_sig
         obs_df.loc[idx, 'sigDec'] = dec_sig
         obs_df.loc[idx, 'sigCorr'] = corr
+        obs_df.loc[idx, 'biasRA'] = 0.0
+        obs_df.loc[idx, 'biasDec'] = 0.0
         obs_df.loc[idx, 'ctr'] = ctr
         obs_df.loc[idx, 'sys'] = sys
         obs_df.loc[idx, 'pos1'] = data['x_gaia_geocentric']
@@ -487,9 +489,8 @@ def get_packed_prog_id(code):
         raise RuntimeError ("Illegal packed prog ID " + code + " in xml")
     try:
         index = unpack_letters[code[1]]
-        if code[0] == 1:
-            index += 52
-        packed = prog_codes[index]
+        index += 62*unpack_letters[code[0]]
+        packed = prog_codes[index] if index <= 93 else ' '
     except Exception as exc:
         raise RuntimeError ("Illegal packed prog ID " + code + " in xml") from exc
     return packed
@@ -514,8 +515,8 @@ def get_unpacked_prog_id(code):
         index = prog_codes.index(code)
     except Exception as exc:
         raise RuntimeError("Illegal program code " + code) from exc
-    first = '1' if index > 51 else '0'
-    second = pack_letters[index%52]
+    first = '1' if index > 61 else '0'
+    second = pack_letters[index%62]
     return first + second
 
 def apply_station_weight_rules(group, obs_df, cols, verbose):

grss/prop/prop_parallel.py

@@ -311,6 +311,8 @@ def _reconstruct_impact_params(row):
     impact = libgrss.ImpactParameters()
     impact.t = row['t']
     impact.xRel = row['xRel']
+    impact.tMap = row['tMap']
+    impact.xRelMap = row['xRelMap']
     impact.dist = row['dist']
     impact.vel = row['vel']
     impact.vInf = row['vInf']

grss/prop/prop_utils.py

@@ -128,7 +128,10 @@ def get_scale_factor(body_id):
     units : str
         units of the scale factor
     """
-    if body_id in {1,199,"Mercury Barycenter"}:
+    if body_id in {10,"Sun"}:
+        scale_factor = 696000.0
+        units = r"R$_\odot$"
+    elif body_id in {1,199,"Mercury Barycenter"}:
         scale_factor = 2440.53
         units = "R$_{Mercury}$"
     elif body_id in {2,299,"Venus Barycenter"}:
@@ -137,6 +140,9 @@ def get_scale_factor(body_id):
     elif body_id in {399,"Earth"}:
         scale_factor = 6378.137
         units = r"R$_\oplus$"
+    elif body_id in {301,"Moon"}:
+        scale_factor = 1737.4
+        units = "R$_{Moon}$"
     elif body_id in {4,499,"Mars Barycenter"}:
         scale_factor = 3396.19
         units = "R$_{Mars}$"
@@ -606,9 +612,11 @@ def plot_bplane(ca_list, plot_offset=False, scale_coords=False, n_std=3, units_k
         # t_dev computed in partials_to_ellipse
         t_map_mean = ca_list[0].tMap
     else:
-        t_mean = np.nan
-        t_dev = np.nan
-        t_map_mean = np.nan
+        print("WARNING: No sigma points, analytic info, or enough data",
+                "to reliably compute mean time")
+        t_mean = np.mean(times)
+        t_dev = np.std(times)
+        t_map_mean = np.mean(map_times)
     if impact_any:
         t_mean_str = Time(t_mean, format='mjd', scale='tdb').utc.iso
         t_map_mean = Time(t_map_mean, format='mjd', scale='tdb').utc.iso

grss/version.txt

@@ -1 +1 @@
-4.0.0
+4.0.1