Merge branch 'main' into 273-joss-paper

.github/workflows/pythonpublish.yml

@@ -27,7 +27,7 @@ jobs:
     - name: Build a source distribution
       run: python3 -m build
     - name: Store the distribution packages
-      uses: actions/upload-artifact@v3
+      uses: actions/upload-artifact@v4
       with:
         name: python-package-distributions
         path: dist/
@@ -47,7 +47,7 @@ jobs:
 
     steps:
     - name: Download all the dists
-      uses: actions/download-artifact@v3
+      uses: actions/download-artifact@v4
       with:
         name: python-package-distributions
         path: dist/

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "astrohack"
-version = "0.5.3"
+version = "0.5.7"
 description = "Holography Antenna Commissioning Kit"
 readme = "README.md"
 requires-python = ">= 3.9, < 3.13"

src/astrohack/antenna/antenna_surface.py

@@ -5,11 +5,11 @@
 import toolviper.utils.logger as logger
 
 from astrohack.antenna.ring_panel import RingPanel
-from astrohack.utils import string_to_ascii_file, create_dataset_label
+from astrohack.utils import string_to_ascii_file, create_dataset_label, data_statistics, statistics_to_text
 from astrohack.utils.constants import *
 from astrohack.utils.conversion import to_db
 from astrohack.utils.conversion import convert_unit
-from astrohack.utils.text import add_prefix, bool_to_str, format_frequency
+from astrohack.utils.text import add_prefix, bool_to_str, format_frequency, format_value_unit
 from astrohack.visualization.plot_tools import well_positioned_colorbar, create_figure_and_axes, close_figure, \
     get_proper_color_map
 
@@ -131,6 +131,13 @@ def _read_panel_xds(self, inputxds):
         self.u_axis = inputxds.u.values
         self.v_axis = inputxds.u.values
         self.panel_distribution = inputxds['PANEL_DISTRIBUTION'].values
+        try:
+            self.amplitude_noise = inputxds['AMP_NOISE'].values
+        except KeyError:
+            logger.warning("Input panel file does not have amplitude noise information, noise statistics will be "
+                           "flawed")
+            self.amplitude_noise = np.full_like(self.amplitude, np.nan)
+
         try:
             self.resolution = inputxds.attrs['aperture_resolution']
         except KeyError:
@@ -175,19 +182,42 @@ def _read_xds(self, inputxds):
         self.label = create_dataset_label(inputxds.attrs['ant_name'], inputxds.attrs['ddi'])
 
     def _measure_ring_clip(self, clip_type, clip_level):
+        self.amplitude_noise = np.where(self.rad < self.telescope.diam / 2., np.nan, self.amplitude)
+        self.amplitude_noise = np.where(self.rad < self.telescope.inlim, self.amplitude, self.amplitude_noise)
+
         if clip_type == 'relative':
             clip = clip_level * np.nanmax(self.amplitude)
         elif clip_type == 'absolute':
             clip = clip_level
         elif clip_type == 'sigma':
-            noise = np.where(self.rad < self.telescope.diam / 2., np.nan, self.amplitude)
-            noiserms = np.sqrt(np.nanmean(noise ** 2))
-            clip = clip_level * noiserms
+            noise_stats = data_statistics(self.amplitude_noise)
+            clip = noise_stats['mean'] + clip_level * noise_stats['rms']
+        elif clip_type == 'noise_threshold':
+            clip = self._compute_noise_threshold_clip(clip_level)
         else:
             msg = f'Unrecognized clipping type: {clip_type}'
             raise Exception(msg)
         return clip
 
+    def _compute_noise_threshold_clip(self, threshold, step_multiplier=0.95):
+        noise_stats = data_statistics(self.amplitude_noise)
+
+        in_disk = np.where(self.rad < self.telescope.diam / 2., 1.0, np.nan)
+        in_disk = np.where(self.rad < self.telescope.inlim, np.nan, in_disk)
+        n_in_disk = np.nansum(in_disk)
+        in_disk_amp = in_disk * self.amplitude
+
+        fraction_in = 0
+        clip = noise_stats['max']
+        multiplier = 1.0
+        while fraction_in < threshold:
+            clip *= multiplier
+            data_in = np.where(in_disk_amp > clip, 1.0, 0.0)
+            fraction_in = np.sum(data_in) / n_in_disk
+            multiplier *= step_multiplier
+
+        return clip
+
     def _init_ringed(self, clip_type, clip_level):
         """
         Do the proper initialization and method association for the case of a ringed antenna
@@ -419,10 +449,7 @@ def gain_at_wavelength(self, corrected, wavelength):
         sinsum = np.nansum(np.sin(scaled_phase[dish_mask]))
         real_factor = np.sqrt(cossum**2 + sinsum**2)/np.sum(dish_mask)
 
-        u_fact = (self.u_axis[1] - self.u_axis[0]) / wavelength
-        v_fact = (self.v_axis[1] - self.v_axis[0]) / wavelength
-
-        theo_gain = fourpi * np.abs(u_fact * v_fact)
+        theo_gain = fourpi * self.telescope.diam/wavelength
         real_gain = theo_gain * real_factor
         return to_db(real_gain), to_db(theo_gain)
 
@@ -522,7 +549,11 @@ def plot_amplitude(self, basename, caller, parm_dict):
         else:
             parm_dict['z_lim'] = parm_dict['amplitude_limits']
 
-        title = "Amplitude, min={0:.5f}, max ={1:.5f} V".format(parm_dict['z_lim'][0], parm_dict['z_lim'][1])
+        amp_stats = data_statistics(np.where(self.mask, self.amplitude, np.nan))
+        noise_stats = data_statistics(self.amplitude_noise)
+        title = ('Amplitude, ' + statistics_to_text(amp_stats) + lnbr +
+                 'Noise, ' + statistics_to_text(noise_stats))
+
         plotname = add_prefix(basename, f'{caller}_amplitude')
         parm_dict['unit'] = self.amp_unit
         self._plot_map(plotname, self.amplitude, title, parm_dict)
@@ -723,14 +754,23 @@ def export_screws(self, filename, unit="mm", comment_char='#'):
         """
         outfile = f"# Screw adjustments for {self.telescope.name}'s {self.label}, pol. state {self.pol_state}\n"
         freq = clight/self.wavelength
-        out_freq = format_frequency(freq)
-        outfile += f"# Frequency = {out_freq}{lnbr}"
-        outfile += "# Adjustments are in " + unit + 2 * lnbr
+        rmses = self.get_rms(unit)
+        outfile += f"# Frequency = {format_frequency(freq)}{lnbr}"
+        if unit == 'mm':
+            outfile += f'# Antenna surface RMS before adjustment: {format_value_unit(rmses[0], unit)}\n'
+            outfile += f'# Antenna surface RMS after adjustment: {format_value_unit(rmses[1], unit)}\n'
+        else:
+            mmrms = self.get_rms('mm')
+            outfile += (f'# Antenna surface RMS before adjustment: {format_value_unit(rmses[0], unit)} or '
+                        f'{format_value_unit(mmrms[0], "mm")}\n')
+            outfile += (f'# Antenna surface RMS after adjustment: {format_value_unit(rmses[1], unit)} or '
+                        f'{format_value_unit(mmrms[1], "mm")}\n')
         outfile += "# Lower means away from subreflector" + lnbr
         outfile += "# Raise means toward the subreflector" + lnbr
         outfile += "# LOWER the panel if the number is POSITIVE" + lnbr
         outfile += "# RAISE the panel if the number is NEGATIVE" + lnbr
-        outfile += 2 * lnbr
+        outfile += "# Adjustments are in " + unit + lnbr
+        outfile += lnbr
         spc = ' '
         outfile += f'{comment_char} Panel{2*spc}'
         nscrews = len(self.telescope.screw_description)
@@ -776,6 +816,8 @@ def export_xds(self):
         xds['DEVIATION'] = xr.DataArray(self.deviation, dims=["u", "v"])
         xds['MASK'] = xr.DataArray(self.mask, dims=["u", "v"])
         xds['PANEL_DISTRIBUTION'] = xr.DataArray(self.panel_distribution, dims=["u", "v"])
+        xds['AMP_NOISE'] = xr.DataArray(self.amplitude_noise, dims=["u", "v"])
+        xds['RADIUS'] = xr.DataArray(self.rad, dims=["u", "v"])
 
         coords = {"u": self.u_axis,
                   "v": self.v_axis}

src/astrohack/antenna/panel_fitting.py

@@ -49,6 +49,7 @@ def _build_system(shape):
 #  Mean                           #
 ###################################
 
+
 def _solve_mean_opt(_self, samples):
     """
     Fit panel surface as a simple mean of its points deviation
@@ -82,6 +83,7 @@ def _correct_mean_opt(self, points):
 # Rigid                           #
 ###################################
 
+
 def _solve_rigid_opt(self, samples):
     """
     Fit panel surface using AIPS rigid model, an inclined plane.
@@ -129,6 +131,7 @@ def _correct_rigid_opt(self, points):
 # Flexible                        #
 ###################################
 
+
 def _solve_flexible_opt(self, samples):
     """
     Fit panel surface using AIPS flexible model, WHAT IS THIS MODEL???
@@ -243,6 +246,8 @@ def _correct_full_paraboloid_opt(self, points):
 #######################################
 # Co-rotated paraboloid least squares #
 #######################################
+
+
 def _solve_corotated_lst_sq_opt(self, samples):
     """
     Builds the designer matrix for least squares fitting, and calls the least_squares fitter for a corotated
@@ -268,6 +273,7 @@ def _solve_corotated_lst_sq_opt(self, samples):
     params, _, _, _ = least_squares_jit(matrix, vector)
     return params
 
+
 def _correct_corotated_lst_sq_opt(self, points):
     """
     Provides the correction on a point using the fitted model
@@ -293,6 +299,7 @@ def _correct_corotated_lst_sq_opt(self, points):
 # Co-rotated robust               #
 ###################################
 
+
 def _solve_corotated_robust_opt(self, samples):
     """
     Try fitting the Surface of a panel using the corotated least_squares method, if that fails fallback to scipy
@@ -634,5 +641,3 @@ def get_xcycval(self):
 
     def get_coords(self):
         return self.xc, self.yc, self.ix, self.iy
-
-

src/astrohack/config/extract_holog.param.json

@@ -17,6 +17,11 @@
         "baseline_average_nearest":{
             "type":["int", "str"]
         },
+        "exclude_antennas": {
+            "type":["ndarray", "list", "str"],
+            "struct type": ["int"],
+            "nullable": true
+        },
         "holog_name":{
             "type":["str"],
             "nullable": true

src/astrohack/config/holog.param.json

@@ -51,9 +51,6 @@
         "to_stokes":{
             "type":["boolean"]
         },
-        "apply_mask":{
-            "type":["boolean"]
-        },
         "phase_fit":{
             "type":["boolean", "list", "ndarray"],
             "minlength": 5,

src/astrohack/config/mds.param.json

@@ -954,6 +954,12 @@
             "type": ["string"],
             "check allowed with": "units.frequency"
         },
+        "rms_unit":{
+            "nullable": false,
+            "required": false,
+            "type": ["string"],
+            "check allowed with": "units.length"
+        },
         "parallel":{
             "nullable": false,
             "required": false,

src/astrohack/config/panel.param.json

@@ -10,11 +10,14 @@
         },
         "clip_type":{
             "type":["str"],
-            "allowed": ["absolute", "sigma", "relative"]
+            "allowed": ["absolute", "sigma", "relative", "noise_threshold"],
+            "nullable": false
         },
         "clip_level":{
-            "type":["float", "int"],
-            "min": 0
+            "type":["float", "int", "dict"],
+            "min": 0,
+            "nullable": false,
+            "struct type": ["dict"]
         },
         "panel_model":{
             "type":["str"],