updates

jurjen93 · Jul 18, 2024 · f7d3e06 · f7d3e06
1 parent daba1fe
commit f7d3e06
Show file tree

Hide file tree

Showing 7 changed files with 440 additions and 299 deletions.
diff --git a/catalogue_helpers/associate_components.py b/catalogue_helpers/associate_components.py
diff --git a/catalogue_helpers/random_forest_outlier_detection.py b/catalogue_helpers/random_forest_outlier_detection.py
diff --git a/ms_helpers/concat_with_dummies.py b/ms_helpers/concat_with_dummies.py
@@ -2,7 +2,7 @@
 Concat measurement sets with dummies
 
 Example:
-    python concat_with_dummies.py --ms *.ms --concat_name concat.ms --data_column CORRECTED_DATA
+    python concat_with_dummies.py --concat_name concat.ms --freq_avg 4 --time_avg 8 *.ms
 """
 
 __author__ = "Jurjen de Jong"

diff --git a/ms_helpers/ms_flagger.py b/ms_helpers/ms_flagger.py
@@ -28,20 +28,31 @@ def main():
     if args.timerange is not None:
         min_time, max_time = [f for f in args.timerange.split('-')]
 
+
     for ms in args.ms:
+
+        expr = ''
+
         command = ['DP3',
                    'msout.storagemanager=dysco',
                    f'msin={ms}',
                    f'msout=flagged_{ms.split("/")[-1]}',
                    'steps=[flag]',
                    'flag.type=preflagger',
-                   f'flag.baseline={args.ant}']
+                   f'flag.flag1.baseline={args.ant}']
+        if args.freqrange is not None or args.timerange is not None:
+            expr += 'flag1'
         # frequency flagging
         if args.freqrange is not None:
-            command += [f"flag.freqrange='[{min_freq} .. {max_freq} MHz]'"]
+            command += [f"flag.flag2.freqrange='[{min_freq} .. {max_freq} MHz]'"]
+            expr += ' and flag2'
         # time flagging
         if args.timerange is not None:
-            command += [f"flag.reltime='[{min_time} .. {max_time}]'"]
+            command += [f"flag.flag3.reltime='[{min_time} .. {max_time}]'"]
+            expr += ' and flag3'
+
+        command += [f'flag.expr="{expr}"']
+
 
         print(' '.join(command))
         os.system(' '.join(command))

diff --git a/ms_merger.py b/ms_merger.py
@@ -176,16 +176,14 @@ def time_resolution(resolution, fov_diam, time_smearing=0.95):
     :return: integration time in seconds
     """
 
-    #TODO: Make dependent on frequency
-
     # Convert distance from degrees to radians
     distance_from_phase_center_rad = np.deg2rad(fov_diam/2)
 
     # Calculate angular resolution (radians)
     angular_resolution_rad = resolution*4.8481*1e-6
 
-    int_time = (np.sqrt((1-time_smearing)/1.2288710615597145e-09)/
-                          distance_from_phase_center_rad*angular_resolution_rad)
+    int_time = 2.9*10**4*(angular_resolution_rad*np.sqrt(1-time_smearing)/
+                          distance_from_phase_center_rad)
 
     return int_time
 
@@ -321,7 +319,7 @@ def mjd_seconds_to_lst_seconds(mjd_seconds, longitude_deg=52.909):
     location = EarthLocation(lon=longitude_deg * u.deg, lat=0. * u.deg)
 
     # Calculate LST in hours
-    lst_hours = time_utc.sidereal_time('mean', longitude=location.lon).hour
+    lst_hours = time_utc.sidereal_time('apparent', longitude=location.lon).hour
 
     # Convert LST from hours to seconds
     lst_seconds = lst_hours * 3600.0
@@ -1269,7 +1267,7 @@ def __init__(self, msin: list = None, outname: str = 'empty.ms', chunkmem: float
         self.num_cpus = psutil.cpu_count(logical=True)
         total_memory = psutil.virtual_memory().total / (1024 ** 3)  # in GB
         target_chunk_size = total_memory / chunkmem
-        self.chunk_size = min(int(target_chunk_size * (1024 ** 3) / np.dtype(np.float128).itemsize/2/self.freq_len), 2_000_000)
+        self.chunk_size = min(int(target_chunk_size * (1024 ** 3) / np.dtype(np.float128).itemsize/2/self.freq_len), 500_000)
         print(f"\n---------------\nChunk size ==> {self.chunk_size}")