Merge pull request #38 from thomaswmorris/reparam_pointing

Reparametrize pointing
thomaswmorris · Dec 1, 2023 · 21d5e43 · 21d5e43
2 parents 16a3062 + 81d3fa6
commit 21d5e43
Show file tree

Hide file tree

Showing 61 changed files with 614 additions and 638 deletions.
diff --git a/docs/source/tutorials/Synthetic_AtLAST_observations.ipynb b/docs/source/tutorials/Synthetic_AtLAST_observations.ipynb
diff --git a/docs/source/tutorials/Synthetic_MUSTANG2_observations.ipynb b/docs/source/tutorials/Synthetic_MUSTANG2_observations.ipynb
diff --git a/docs/source/usage.rst b/docs/source/usage.rst
@@ -41,19 +41,19 @@ The arrays
         - max_el_vel (float): ...
         - max_az_acc (float): ...
         - max_el_acc (float): ...
-        - max_baseline (float): meters (only for ALMA)
+        - baseline (float): meters (only for ALMA)
 
 
 Pointings
 ^^^^^^^^^
 
-The next component of the simulation
+The next component of the simulation is the pointing, which defines the direction and scan pattern of the observation:::
 
-    import maria
+    daisy_scan = maria.get_pointing("daisy")
 
-    mustang = maria.get_array("MUSTANG-2")
+We can specify a custom pointing by passing extra arguments to the above, such as::
 
-We can specify a custom array
+    custom_daisy_scan = maria.get_pointing("daisy", integration_time=600, scan_options={"radius": 2})
 
 
 Sites
@@ -114,7 +114,7 @@ Simulation parameters
         - max_el_vel (float): ...
         - max_az_acc (float): ...
         - max_el_acc (float): ...
-        - max_baseline (float): meters (only for ALMA)
+        - baseline (float): meters (only for ALMA)
 
 - **pointing:** Scanning strategy
     - Predefined configurations: `stare`, `daisy`, `BAF`,

diff --git a/maria/array.py b/maria/array.py
@@ -23,10 +23,9 @@
 
 here, this_filename = os.path.split(__file__)
 
+all_array_params = utils.io.read_yaml(f"{here}/configs/params.yml")["array"]
+
 ARRAY_CONFIGS = utils.io.read_yaml(f"{here}/configs/arrays.yml")
-ARRAY_PARAMS = set()
-for key, config in ARRAY_CONFIGS.items():
-    ARRAY_PARAMS |= set(config.keys())
 
 DISPLAY_COLUMNS = ["description", "field_of_view", "primary_size"]
 supported_arrays_table = pd.DataFrame(ARRAY_CONFIGS).T
@@ -41,13 +40,59 @@ def __init__(self, invalid_array):
         )
 
 
-def get_array_config(array_name, **kwargs):
+def generate_array_offsets(geometry, field_of_view, n):
+    valid_array_types = ["flower", "hex", "square"]
+
+    if geometry == "flower":
+        phi = np.pi * (3.0 - np.sqrt(5.0))  # golden angle in radians
+        dzs = np.zeros(n).astype(complex)
+        for i in range(n):
+            dzs[i] = np.sqrt((i / (n - 1)) * 2) * np.exp(1j * phi * i)
+        od = np.abs(np.subtract.outer(dzs, dzs))
+        dzs *= field_of_view / od.max()
+        return np.c_[np.real(dzs), np.imag(dzs)]
+    if geometry == "hex":
+        return generate_hex_offsets(n, field_of_view)
+    if geometry == "square":
+        dxy_ = np.linspace(-field_of_view, field_of_view, int(np.ceil(np.sqrt(n)))) / (
+            2 * np.sqrt(2)
+        )
+        DX, DY = np.meshgrid(dxy_, dxy_)
+        return np.c_[DX.ravel()[:n], DY.ravel()[:n]]
+
+    raise ValueError(
+        "Please specify a valid array type. Valid array types are:\n"
+        + "\n".join(valid_array_types)
+    )
+
+
+def generate_hex_offsets(n, d):
+    angles = np.linspace(0, 2 * np.pi, 6 + 1)[1:] + np.pi / 2
+    zs = np.array([0])
+    layer = 0
+    while len(zs) < n:
+        for angle in angles:
+            for z in layer * np.exp(1j * angle) + np.arange(layer) * np.exp(
+                1j * (angle + 2 * np.pi / 3)
+            ):
+                zs = np.append(zs, z)
+        layer += 1
+    zs -= zs.mean()
+    zs *= 0.5 * d / np.abs(zs).max()
+
+    return np.c_[np.real(np.array(zs[:n])), np.imag(np.array(zs[:n]))]
+
+
+def get_array_config(array_name=None, **kwargs):
     if array_name not in ARRAY_CONFIGS.keys():
         raise InvalidArrayError(array_name)
-    ARRAY_CONFIG = ARRAY_CONFIGS[array_name].copy()
+    array_config = ARRAY_CONFIGS[array_name].copy()
     for k, v in kwargs.items():
-        ARRAY_CONFIG[k] = v
-    return ARRAY_CONFIG
+        if k in all_array_params.keys():
+            array_config[k] = v
+        else:
+            raise ValueError(f"'{k}' is not a valid argument for an array!")
+    return array_config
 
 
 def get_array(array_name, **kwargs):
@@ -64,7 +109,7 @@ def generate_dets_from_config(
     bands: Mapping,
     field_of_view: float,
     geometry: str = "hex",
-    max_baseline: float = 0,
+    baseline: float = 0,
     randomize_offsets: bool = True,
 ):
     dets = pd.DataFrame(
@@ -89,23 +134,26 @@ def generate_dets_from_config(
         band_dets.loc[:, "band_center"] = band_config["band_center"]
         band_dets.loc[:, "band_width"] = band_config["band_width"]
 
-        dets = pd.concat([dets, band_dets])
-        dets.index = np.arange(len(dets))
+        det_offsets_radians = np.radians(
+            generate_array_offsets(geometry, field_of_view, len(band_dets))
+        )
 
-    offsets_radians = np.radians(
-        utils.generate_array_offsets(geometry, field_of_view, len(dets))
-    )
+        # should we make another function for this?
+        det_baselines_meters = generate_array_offsets(
+            geometry, baseline, len(band_dets)
+        )
 
-    # should we make another function for this?
-    baseline = utils.generate_array_offsets(geometry, max_baseline, len(dets))
+        # if randomize_offsets:
+        #     np.random.shuffle(offsets_radians)  # this is a stupid function.
 
-    if randomize_offsets:
-        np.random.shuffle(offsets_radians)  # this is a stupid function.
+        band_dets.loc[:, "offset_x"] = det_offsets_radians[:, 0]
+        band_dets.loc[:, "offset_y"] = det_offsets_radians[:, 1]
+        band_dets.loc[:, "baseline_x"] = det_baselines_meters[:, 0]
+        band_dets.loc[:, "baseline_y"] = det_baselines_meters[:, 1]
+        band_dets.loc[:, "baseline_z"] = 0
 
-    dets.loc[:, "offset_x"] = offsets_radians[:, 0]
-    dets.loc[:, "offset_y"] = offsets_radians[:, 1]
-    dets.loc[:, "baseline_x"] = baseline[:, 0]
-    dets.loc[:, "baseline_y"] = baseline[:, 1]
+        dets = pd.concat([dets, band_dets])
+        dets.index = np.arange(len(dets))
 
     for key in ["offset_x", "offset_y", "baseline_x", "baseline_y"]:
         dets.loc[:, key] = dets.loc[:, key].astype(float)
@@ -123,6 +171,7 @@ class Array:
     primary_size: float = 5  # in meters
     field_of_view: float = 1  # in deg
     geometry: str = "hex"
+    baseline: float = 0
     max_az_vel: float = 0  # in deg/s
     max_el_vel: float = np.inf  # in deg/s
     max_az_acc: float = 0  # in deg/s^2
@@ -173,14 +222,19 @@ def from_config(cls, config):
         if isinstance(config["dets"], Mapping):
             field_of_view = config.get("field_of_view", 1)
             geometry = config.get("geometry", "hex")
+            baseline = config.get("baseline", 0)  # default to zero baseline
             dets = generate_dets_from_config(
-                config["dets"], field_of_view=field_of_view, geometry=geometry
+                config["dets"],
+                field_of_view=field_of_view,
+                geometry=geometry,
+                baseline=baseline,
             )
 
         return cls(
             description=config["description"],
             primary_size=config["primary_size"],
             field_of_view=field_of_view,
+            baseline=baseline,
             geometry=geometry,
             max_az_vel=config["max_az_vel"],
             max_el_vel=config["max_el_vel"],
@@ -219,7 +273,7 @@ def passbands(self, nu):
         return nu_mask.astype(float) / nu_mask.sum(axis=-1)[:, None]
 
     def angular_fwhm(self, z):
-        return utils.gaussian_beam_angular_fwhm(
+        return utils.beam.gaussian_beam_angular_fwhm(
             z=z,
             w_0=self.primary_size / np.sqrt(2 * np.log(2)),
             f=self.dets.band_center.values,
@@ -253,7 +307,7 @@ def plot_dets(self):
             band_res_arcmins = 60 * np.degrees(
                 1.22 * 2.998e8 / (1e9 * nom_freq * self.primary_size)
             )
-            # band_res_arcmins = 60 * np.degrees(utils.gaussian_beam_angular_fwhm(z=1e12, w_0=self.primary_size, f=nom_freq))
+
             offsets_arcmins = 60 * np.degrees(self.offsets[band_mask])
 
             ax.add_collection(

diff --git a/maria/atmosphere/base.py b/maria/atmosphere/base.py
@@ -83,7 +83,7 @@ def __init__(self, array, pointing, site, verbose=False, **kwargs):
 
     @property
     def EL(self):
-        return utils.xy_to_lonlat(
+        return utils.coords.xy_to_lonlat(
             self.array.offset_x[:, None],
             self.array.offset_y[:, None],
             self.pointing.az,