From 134c5ab5e5a16f9812b322868ffd5e02827123c0 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Sun, 19 Feb 2023 23:19:20 -0500
Subject: [PATCH 001/109] CrossValidate: add model to self

---
 src/mpol/crossval.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/mpol/crossval.py b/src/mpol/crossval.py
index b5aa81da..3fb76664 100644
--- a/src/mpol/crossval.py
+++ b/src/mpol/crossval.py
@@ -177,11 +177,11 @@ def run_crossval(self, dataset):
             #     train_set, test_set = train_set.to(self._device), test_set.to(self._device)
             
             # create a new model and optimizer for this k_fold
-            model = SimpleNet(coords=self._coords, nchan=self._gridder.nchan)
+            self.model = SimpleNet(coords=self._coords, nchan=self._gridder.nchan)
             # if hasattr(self._device,'type') and self._device.type == 'cuda': # TODO: confirm which objects need to be passed to gpu
-            #     model = model.to(self._device)
+            #     self.model = self.model.to(self._device)
 
-            optimizer = torch.optim.Adam(model.parameters(), lr=self._learn_rate)
+            optimizer = torch.optim.Adam(self.model.parameters(), lr=self._learn_rate)
 
             trainer = TrainTest(gridder=self._gridder, 
                                 optimizer=optimizer, 
@@ -199,10 +199,10 @@ def run_crossval(self, dataset):
                                 verbose=self._verbose
             )
 
-            loss, loss_history = trainer.train(model, train_set)
+            loss, loss_history = trainer.train(self.model, train_set)
             if self._store_cv_diagnostics:
                 self._cv_diagnostics['loss_histories'].append(loss_history)
-            all_scores.append(trainer.test(model, test_set))
+            all_scores.append(trainer.test(self.model, test_set))
 
         # average individual test scores to get the cross-val metric for chosen 
         # hyperparameters

From 95bf0e0113ab39f221919d8fc5a2829a7686babc Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Sun, 19 Feb 2023 23:19:45 -0500
Subject: [PATCH 002/109] CrossValidate: rename var

---
 src/mpol/crossval.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/mpol/crossval.py b/src/mpol/crossval.py
index 3fb76664..42a6e269 100644
--- a/src/mpol/crossval.py
+++ b/src/mpol/crossval.py
@@ -165,7 +165,7 @@ def run_crossval(self, dataset):
         """
         all_scores = []
         if self._store_cv_diagnostics:
-            self._cv_diagnostics = defaultdict(list) 
+            self._diagnostics = defaultdict(list) 
 
         split_iterator = self.split_dataset(dataset)
         for kk, (train_set, test_set) in enumerate(split_iterator):
@@ -201,7 +201,7 @@ def run_crossval(self, dataset):
 
             loss, loss_history = trainer.train(self.model, train_set)
             if self._store_cv_diagnostics:
-                self._cv_diagnostics['loss_histories'].append(loss_history)
+                self._diagnostics['loss_histories'].append(loss_history)
             all_scores.append(trainer.test(self.model, test_set))
 
         # average individual test scores to get the cross-val metric for chosen 
@@ -213,9 +213,9 @@ def run_crossval(self, dataset):
         return cv_score
     
     @property
-    def cv_diagnostics(self):
+    def diagnostics(self):
         """Dict containing diagnostics of the cross-validation loop"""
-        return self._cv_diagnostics
+        return self._diagnostics
 
 
 class RandomCellSplitGridded:

From 7e1823710076745f9ddcbd551dee22285e6b549b Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Sun, 19 Feb 2023 23:21:38 -0500
Subject: [PATCH 003/109] add onedim.py, func to get 1d vis fit

---
 src/mpol/onedim.py | 38 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 38 insertions(+)
 create mode 100644 src/mpol/onedim.py

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
new file mode 100644
index 00000000..db340e2c
--- /dev/null
+++ b/src/mpol/onedim.py
@@ -0,0 +1,38 @@
+import numpy as np 
+import matplotlib.pyplot as plt 
+
+from mpol.fourier import NuFFT
+from mpol.utils import torch2npy
+
+
+def get_1d_vis_fit(model, u, v, chan=0):
+    r"""
+    Obtain the 1D (radial) visibility model V(q) corresponding to a 2D MPoL 
+    image-domain model. 
+
+    Parameters
+    ----------
+    model : `torch.nn.Module` object
+        Instance of the `mpol.precomposed.SimpleNet` class
+    u : array, unit=:math:[`k\lambda`] 
+        u-coordinates at which to sample (e.g., those of the dataset)
+    v : array, unit=:math:[`k\lambda`]
+        v-coordinates at which to sample (e.g., those of the dataset)
+    chan : int, default=0
+        Channel of `model` to select
+
+    Returns
+    -------
+    q : array, unit=:math:[`k\lambda`]
+        Baselines corresponding to `u` and `v`
+    Vmod : array, unit=[Jy] # TODO: right unit?
+        Visibility amplitudes at `q`
+    """
+    q = np.hypot(u, v)
+
+    nufft = NuFFT(coords=model.coords, nchan=model.nchan, uu=u, vv=v)
+    # get model visibilities 
+    Vmod = nufft(model.icube()).detach()[chan] 
+    
+    return q, Vmod
+

From a5f636505fc8412eb2743f926a07041aaf7bc1aa Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Sun, 19 Feb 2023 23:28:56 -0500
Subject: [PATCH 004/109] add 1d radial profile func

---
 src/mpol/onedim.py | 43 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 43 insertions(+)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index db340e2c..525629d1 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -36,3 +36,46 @@ def get_1d_vis_fit(model, u, v, chan=0):
     
     return q, Vmod
 
+
+def get_radial_profile(model, l_center=0.0, m_center=0.0, bins=None, chan=0):     
+    r"""
+    Obtain a 1D (radial) brightness profile I(r) from an MPoL model. 
+
+    Parameters
+    ----------
+    model : `torch.nn.Module` object
+        Instance of the `mpol.precomposed.SimpleNet` class
+    l_center : float, default=0.0, unit=[arcsec]
+        Image center along l-axis. If None, the image center pixel will be used
+    m_center : float, default=0.0, unit=[arcsec]
+        Image center along m-axis. If None, the image center pixel will be used
+    bins : array, default=None, unit=[arcsec]
+        Radial bin edges to use in calculating I(r)
+    chan : int, default=0
+        Channel of `model` to select
+
+    Returns
+    -------
+    bin_centers : array, unit=[arcsec]
+        Radial coordinates of image at center of `bins`
+    Is : array, unit=[Jy / arcsec^2]
+        Azimuthally averaged pixel brightness at `rs`
+
+    """
+    skycube = torch2npy(model.icube.sky_cube)[chan]
+
+    # get image center in Cartesian [arcsec]
+    xx, yy = model.coords.sky_x_centers_2D, model.coords.sky_y_centers_2D
+    # shift center coordinate
+    xshift, yshift = xx - l_center, yy - m_center
+    rshift = np.hypot(xshift, yshift.T) 
+    rshift = rshift.flatten()
+
+    if bins is None:
+        step = np.hypot(model.coords.cell_size, model.coords.cell_size)
+        bins = np.arange(0.0, max(rshift), step)
+
+    Is, bin_edges = np.histogram(a=rshift, bins=bins, weights=skycube.flatten())
+    bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
+
+    return bin_centers, Is
\ No newline at end of file

From 0c0ffc831b8e9897c6d5332f5f97b4bf46df0bfc Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Sun, 19 Feb 2023 23:44:35 -0500
Subject: [PATCH 005/109] get_radial_profile: simplify args

---
 src/mpol/onedim.py | 13 ++++++-------
 1 file changed, 6 insertions(+), 7 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 525629d1..32f62ec2 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -37,18 +37,17 @@ def get_1d_vis_fit(model, u, v, chan=0):
     return q, Vmod
 
 
-def get_radial_profile(model, l_center=0.0, m_center=0.0, bins=None, chan=0):     
+def get_radial_profile(model, center=(0.0, 0.0), bins=None, chan=0):
     r"""
-    Obtain a 1D (radial) brightness profile I(r) from an MPoL model. 
+    Obtain a 1D (radial) brightness profile I(r) from an MPoL model.
 
     Parameters
     ----------
     model : `torch.nn.Module` object
         Instance of the `mpol.precomposed.SimpleNet` class
-    l_center : float, default=0.0, unit=[arcsec]
-        Image center along l-axis. If None, the image center pixel will be used
-    m_center : float, default=0.0, unit=[arcsec]
-        Image center along m-axis. If None, the image center pixel will be used
+    center : 2-tuple of float, default=(0.0, 0.0), unit=[arcsec]
+        Offset (RA, Dec) of source in image. Postive RA offset is west of
+        north. If None, the source is assumed to be at the image center pixel
     bins : array, default=None, unit=[arcsec]
         Radial bin edges to use in calculating I(r)
     chan : int, default=0
@@ -67,7 +66,7 @@ def get_radial_profile(model, l_center=0.0, m_center=0.0, bins=None, chan=0):
     # get image center in Cartesian [arcsec]
     xx, yy = model.coords.sky_x_centers_2D, model.coords.sky_y_centers_2D
     # shift center coordinate
-    xshift, yshift = xx - l_center, yy - m_center
+    xshift, yshift = xx - center[0], yy - center[1]
     rshift = np.hypot(xshift, yshift.T) 
     rshift = rshift.flatten()
 

From aa5dd0ca56212af641c37de82e1780db003a8f7f Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Mon, 20 Feb 2023 18:26:30 -0500
Subject: [PATCH 006/109] typos

---
 src/mpol/geometry.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/mpol/geometry.py b/src/mpol/geometry.py
index 1ce4e8de..ddca4715 100644
--- a/src/mpol/geometry.py
+++ b/src/mpol/geometry.py
@@ -70,7 +70,7 @@ def flat_to_observer(x, y, omega=None, incl=None, Omega=None):
 
 
 def observer_to_flat(X, Y, omega=None, incl=None, Omega=None):
-    """Rotate the from to convert a point in the observer frame (X,Y,Z) to the flat (x,y,z) frame.
+    """Rotate the frame to convert a point in the observer frame (X,Y,Z) to the flat (x,y,z) frame.
 
     It is assumed that the +Z axis points *towards* the observer. The rotation operations are the inverse of the :func:`~mpol.geometry.flat_to_observer` operations.
 
@@ -90,7 +90,7 @@ def observer_to_flat(X, Y, omega=None, incl=None, Omega=None):
         Omega (torch float tensor): A tensor representing the position angle of the ascending node in [radians]. Default 0.0
 
     Returns:
-        Two tensors representing ``(x, y)`` in the observer frame.
+        Two tensors representing ``(x, y)`` in the flat frame.
     """
 
     # Rotation matrices result in a *clockwise* rotation of the axes, as defined using the righthand rule.

From 1d288c263f9224335f567532f9761b246f6a2a99 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Mon, 20 Feb 2023 18:32:53 -0500
Subject: [PATCH 007/109] add deproject_vis routine

---
 src/mpol/geometry.py | 69 +++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 68 insertions(+), 1 deletion(-)

diff --git a/src/mpol/geometry.py b/src/mpol/geometry.py
index ddca4715..b45dfad5 100644
--- a/src/mpol/geometry.py
+++ b/src/mpol/geometry.py
@@ -1,4 +1,4 @@
-"""The geometry package provides routines for projecting and de-projecting sky images.
+"""The geometry package provides routines for projecting and de-projecting sky images and visibilities.
 """
 
 import numpy as np
@@ -129,3 +129,70 @@ def observer_to_flat(X, Y, omega=None, incl=None, Omega=None):
         y = y1
 
     return x, y
+
+
+# TODO
+def deproject_vis(u, v, V, weights, source_geom, inverse=False, rescale_flux=True): 
+    r"""
+    Deproject (or reproject) visibilities (and optionally rescale their 
+    amplitudes and weights, according to the source inclination) using routines 
+    in `frank <https://github.com/discsim/frank>`_
+
+    Parameters
+    ----------
+    u : array of real, size = N, unit = :math:`\lambda`
+        u-points of the visibilities
+    v : array of real, size = N, unit = :math:`\lambda`
+        v-points of the visibilities
+    V : array of real, size = N, unit = Jy
+        Complex visibilites
+    weights : array of real, size = N, unit = Jy
+        Weights on the visibilities    
+    source_geom : dict
+        Dictionary with source geometry parameters. Keys:
+            "inc" : float, unit = deg
+                Inclination
+            "PA" : float, unit = deg
+                Position angle, defined east of north.  
+            "dRA" : float, unit = arcsec
+                Phase centre offset in right ascension.
+            "dDec" : float, units = arcsec
+                Phase centre offset in declination.           
+    inverse : bool, default=False
+        If True, the uv-points are reprojected rather than deprojected
+    rescale_flux : bool, default=True
+        If True, the visibility amplitudes and weights are rescaled to account 
+        for the difference between the inclined (observed) brightness and the 
+        assumed face-on brightness, assuming the emission is optically thick. 
+        The source's integrated (2D) flux is assumed to be:
+            :math:`F = \cos(i) \int_r^{r=R}{I(r) 2 \pi r dr}`.
+        No rescaling would be appropriate in the optically thin limit.
+    Returns
+    -------
+    up : array of real, size = N, unit = :math:`\lambda`
+        Corrected u-points of the visibilities
+    vp : array of real, size = N, unit = :math:`\lambda`
+        Corrected v-points of the visibilities
+    Vp : array of real, size = N, unit = Jy
+        Corrected complex visibilites
+    weights_scaled : array of real, size = N, unit = Jy
+        Rescaled weights on the visibilities        
+
+    """
+    from frank.geometry import FixedGeometry
+
+    geom = FixedGeometry(**source_geom)
+
+    weights_scaled = weights * 1
+
+    if inverse:
+        up, vp, Vp = geom.apply_correction(u, v, V)
+        if rescale_flux:
+            Vp, weights_scaled = geom.rescale_total_flux(V, weights)
+
+    else:
+        if rescale_flux:
+            Vp, weights_scaled = geom.rescale_total_flux(V, weights)
+         up, vp, Vp = geom.undo_correction(u, v, V)
+
+    return up, vp, Vp, weights_scaled

From 9b43c9b7cb7d93c89ec13599ac080c3f59142948 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Mon, 20 Feb 2023 18:33:28 -0500
Subject: [PATCH 008/109] setup.py: add new type of EXTRA, with frank

---
 setup.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/setup.py b/setup.py
index a700ed3c..adde423e 100644
--- a/setup.py
+++ b/setup.py
@@ -52,10 +52,13 @@ def get_version(rel_path):
         "jupyter-cache",
         "Pillow",
     ],
+    "analysis": [
+        "frank>=1.2.1",
+    ],
 }
 
 EXTRA_REQUIRES["dev"] = (
-    EXTRA_REQUIRES["test"] + EXTRA_REQUIRES["docs"] + ["pylint", "black", "pre-commit"]
+    EXTRA_REQUIRES["test"] + EXTRA_REQUIRES["docs"] + EXTRA_REQUIRES["docs"] + ["pylint", "black", "pre-commit"]
 )
 
 

From 8ac919cf45ae9295a1fbc3d663f9d03a045ea0b8 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Mon, 20 Feb 2023 18:34:27 -0500
Subject: [PATCH 009/109] get_radial_profile: add comments, stage deprojec.

---
 src/mpol/onedim.py | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 32f62ec2..c971e2c4 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -43,6 +43,7 @@ def get_radial_profile(model, center=(0.0, 0.0), bins=None, chan=0):
 
     Parameters
     ----------
+    # TODO
     model : `torch.nn.Module` object
         Instance of the `mpol.precomposed.SimpleNet` class
     center : 2-tuple of float, default=(0.0, 0.0), unit=[arcsec]
@@ -61,13 +62,16 @@ def get_radial_profile(model, center=(0.0, 0.0), bins=None, chan=0):
         Azimuthally averaged pixel brightness at `rs`
 
     """
+    # observer_to_flat(X, Y, omega, incl, Omega) # TODO
+
     skycube = torch2npy(model.icube.sky_cube)[chan]
 
-    # get image center in Cartesian [arcsec]
+    # Cartesian pixel coordinates [arcsec]
     xx, yy = model.coords.sky_x_centers_2D, model.coords.sky_y_centers_2D
-    # shift center coordinate
+    # shift image center
     xshift, yshift = xx - center[0], yy - center[1]
-    rshift = np.hypot(xshift, yshift.T) 
+    # radial pixel coordinates
+    rshift = np.hypot(xshift, yshift) 
     rshift = rshift.flatten()
 
     if bins is None:
@@ -77,4 +81,4 @@ def get_radial_profile(model, center=(0.0, 0.0), bins=None, chan=0):
     Is, bin_edges = np.histogram(a=rshift, bins=bins, weights=skycube.flatten())
     bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
 
-    return bin_centers, Is
\ No newline at end of file
+    return bin_centers, Is

From 30556396fb465a4f6e2a2f36bfc31ab64a98a75a Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Mon, 20 Feb 2023 18:35:11 -0500
Subject: [PATCH 010/109] get_1d_vis_fit: stage vis deprojection

---
 src/mpol/onedim.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index c971e2c4..7e1bc65d 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -1,7 +1,7 @@
 import numpy as np 
-import matplotlib.pyplot as plt 
 
 from mpol.fourier import NuFFT
+from mpol.geometry import deproject_vis, observer_to_flat
 from mpol.utils import torch2npy
 
 
@@ -12,6 +12,7 @@ def get_1d_vis_fit(model, u, v, chan=0):
 
     Parameters
     ----------
+    # TODO 
     model : `torch.nn.Module` object
         Instance of the `mpol.precomposed.SimpleNet` class
     u : array, unit=:math:[`k\lambda`] 
@@ -25,14 +26,16 @@ def get_1d_vis_fit(model, u, v, chan=0):
     -------
     q : array, unit=:math:[`k\lambda`]
         Baselines corresponding to `u` and `v`
-    Vmod : array, unit=[Jy] # TODO: right unit?
+    Vmod : array, unit=[Jy] 
         Visibility amplitudes at `q`
     """
+    # deproject_vis(u, v, V, weights, source_geom, inverse, rescale_flux) # TODO
+
     q = np.hypot(u, v)
 
     nufft = NuFFT(coords=model.coords, nchan=model.nchan, uu=u, vv=v)
     # get model visibilities 
-    Vmod = nufft(model.icube()).detach()[chan] 
+    Vmod = nufft(model.icube()).detach()[chan]
     
     return q, Vmod
 

From 0b6da4b3dbf4caaaf7479582281a8addbf566063 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Mon, 20 Feb 2023 23:20:27 -0500
Subject: [PATCH 011/109] move frank import internal to func

---
 src/mpol/onedim.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 7e1bc65d..aa166b02 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -1,7 +1,7 @@
 import numpy as np 
 
 from mpol.fourier import NuFFT
-from mpol.geometry import deproject_vis, observer_to_flat
+from mpol.geometry import observer_to_flat
 from mpol.utils import torch2npy
 
 
@@ -28,7 +28,12 @@ def get_1d_vis_fit(model, u, v, chan=0):
         Baselines corresponding to `u` and `v`
     Vmod : array, unit=[Jy] 
         Visibility amplitudes at `q`
+
+    Notes
+    -----
+    This routine requires the `frank <https://github.com/discsim/frank>`_ package 
     """
+    from mpol.geometry import deproject_vis
     # deproject_vis(u, v, V, weights, source_geom, inverse, rescale_flux) # TODO
 
     q = np.hypot(u, v)

From ce507939d1984add10941e79266a41bc6707de2b Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Mon, 20 Feb 2023 23:47:06 -0500
Subject: [PATCH 012/109] get_radial_profile: bug fix (norm radial profile)

---
 src/mpol/geometry.py | 2 +-
 src/mpol/onedim.py   | 7 ++++++-
 2 files changed, 7 insertions(+), 2 deletions(-)

diff --git a/src/mpol/geometry.py b/src/mpol/geometry.py
index b45dfad5..685f600d 100644
--- a/src/mpol/geometry.py
+++ b/src/mpol/geometry.py
@@ -193,6 +193,6 @@ def deproject_vis(u, v, V, weights, source_geom, inverse=False, rescale_flux=Tru
     else:
         if rescale_flux:
             Vp, weights_scaled = geom.rescale_total_flux(V, weights)
-         up, vp, Vp = geom.undo_correction(u, v, V)
+        up, vp, Vp = geom.undo_correction(u, v, V)
 
     return up, vp, Vp, weights_scaled
diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index aa166b02..47cb01d6 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -86,7 +86,12 @@ def get_radial_profile(model, center=(0.0, 0.0), bins=None, chan=0):
         step = np.hypot(model.coords.cell_size, model.coords.cell_size)
         bins = np.arange(0.0, max(rshift), step)
 
-    Is, bin_edges = np.histogram(a=rshift, bins=bins, weights=skycube.flatten())
+    # get number of points in each radial bin
+    bin_counts, bin_edges = np.histogram(a=rshift, bins=bins, weights=None)
+    # get radial brightness
+    Is, _ = np.histogram(a=rshift, bins=bins, weights=skycube.flatten())
+    Is /= bin_counts
+    
     bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
 
     return bin_centers, Is

From 5437e0cbec3eae99b3d2be6b36493235c15fe30f Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Tue, 21 Feb 2023 11:21:12 -0500
Subject: [PATCH 013/109] deproject_vis: add TODO

---
 src/mpol/geometry.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpol/geometry.py b/src/mpol/geometry.py
index 685f600d..3fee3c08 100644
--- a/src/mpol/geometry.py
+++ b/src/mpol/geometry.py
@@ -131,7 +131,7 @@ def observer_to_flat(X, Y, omega=None, incl=None, Omega=None):
     return x, y
 
 
-# TODO
+# TODO: check giving correct result 
 def deproject_vis(u, v, V, weights, source_geom, inverse=False, rescale_flux=True): 
     r"""
     Deproject (or reproject) visibilities (and optionally rescale their 

From b4c318cf080c53314294acde3aacd31a88c7e230 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 22 Feb 2023 09:38:23 -0500
Subject: [PATCH 014/109] typo

---
 src/mpol/crossval.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpol/crossval.py b/src/mpol/crossval.py
index bf263fe5..cbe9e1fa 100644
--- a/src/mpol/crossval.py
+++ b/src/mpol/crossval.py
@@ -222,7 +222,7 @@ def run_crossval(self, dataset):
 
             loss, loss_history = trainer.train(self.model, train_set)
             if self._store_cv_diagnostics:
-                self._cv_diagnostics["loss_histories"].append(loss_history)
+                self._diagnostics["loss_histories"].append(loss_history)
             all_scores.append(trainer.test(self.model, test_set))
 
         # average individual test scores to get the cross-val metric for chosen

From dd7324fe0da01479f0e51421e742b44349b92b00 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 22 Feb 2023 15:12:48 -0500
Subject: [PATCH 015/109] get_1d_vis_fit: only intake `model`

---
 src/mpol/onedim.py | 18 +++++-------------
 1 file changed, 5 insertions(+), 13 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 47cb01d6..f3ac7cd4 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -5,20 +5,15 @@
 from mpol.utils import torch2npy
 
 
-def get_1d_vis_fit(model, u, v, chan=0):
+def get_1d_vis_fit(model, chan=0):
     r"""
     Obtain the 1D (radial) visibility model V(q) corresponding to a 2D MPoL 
     image-domain model. 
 
     Parameters
     ----------
-    # TODO 
     model : `torch.nn.Module` object
         Instance of the `mpol.precomposed.SimpleNet` class
-    u : array, unit=:math:[`k\lambda`] 
-        u-coordinates at which to sample (e.g., those of the dataset)
-    v : array, unit=:math:[`k\lambda`]
-        v-coordinates at which to sample (e.g., those of the dataset)
     chan : int, default=0
         Channel of `model` to select
 
@@ -31,17 +26,14 @@ def get_1d_vis_fit(model, u, v, chan=0):
 
     Notes
     -----
-    This routine requires the `frank <https://github.com/discsim/frank>`_ package 
+    This routine requires the `frank <https://github.com/discsim/frank>`_ package # TODO
     """
-    from mpol.geometry import deproject_vis
+    # from mpol.geometry import deproject_vis # TODO
     # deproject_vis(u, v, V, weights, source_geom, inverse, rescale_flux) # TODO
 
-    q = np.hypot(u, v)
+    q = model.coords.packed_q_centers_2D.ravel()
+    Vmod = model.fcube.vis.detach()[chan].ravel()
 
-    nufft = NuFFT(coords=model.coords, nchan=model.nchan, uu=u, vv=v)
-    # get model visibilities 
-    Vmod = nufft(model.icube()).detach()[chan]
-    
     return q, Vmod
 
 

From 8d93ec3142f54129f2d44b2d7775fbe0c69e6fa3 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 22 Feb 2023 15:26:55 -0500
Subject: [PATCH 016/109] setup.py: copy/paste typo

---
 setup.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/setup.py b/setup.py
index adde423e..d700b6b0 100644
--- a/setup.py
+++ b/setup.py
@@ -58,7 +58,7 @@ def get_version(rel_path):
 }
 
 EXTRA_REQUIRES["dev"] = (
-    EXTRA_REQUIRES["test"] + EXTRA_REQUIRES["docs"] + EXTRA_REQUIRES["docs"] + ["pylint", "black", "pre-commit"]
+    EXTRA_REQUIRES["test"] + EXTRA_REQUIRES["docs"] + EXTRA_REQUIRES["analysis"] + ["pylint", "black", "pre-commit"]
 )
 
 

From a98f71170778c7b2230048d56cf93190d4db1814 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 22 Feb 2023 15:38:22 -0500
Subject: [PATCH 017/109] CrossValidate: hold off on hyperpar changes

---
 src/mpol/crossval.py | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/src/mpol/crossval.py b/src/mpol/crossval.py
index 48fcdc08..f4ce77d3 100644
--- a/src/mpol/crossval.py
+++ b/src/mpol/crossval.py
@@ -181,7 +181,7 @@ def run_crossval(self, dataset):
         """
         all_scores = []
         if self._store_cv_diagnostics:
-            self._diagnostics = defaultdict(list)
+            self._cv_diagnostics = defaultdict(list)
 
         split_iterator = self.split_dataset(dataset)
         if self._split_diag_fig:
@@ -198,11 +198,11 @@ def run_crossval(self, dataset):
             #     train_set, test_set = train_set.to(self._device), test_set.to(self._device)
 
             # create a new model and optimizer for this k_fold
-            self.model = SimpleNet(coords=self._coords, nchan=self._imager.nchan)
+            model = SimpleNet(coords=self._coords, nchan=self._imager.nchan)
             # if hasattr(self._device,'type') and self._device.type == 'cuda': # TODO: confirm which objects need to be passed to gpu
-            #     self.model = self.model.to(self._device)
+            #     model = model.to(self._device)
 
-            optimizer = torch.optim.Adam(self.model.parameters(), lr=self._learn_rate)
+            optimizer = torch.optim.Adam(model.parameters(), lr=self._learn_rate)
 
             trainer = TrainTest(
                 imager=self._imager,
@@ -223,13 +223,13 @@ def run_crossval(self, dataset):
                 verbose=self._verbose,
             )
 
-            loss, loss_history = trainer.train(self.model, train_set)
+            loss, loss_history = trainer.train(model, train_set)
             if self._store_cv_diagnostics:
-                self._diagnostics["loss_histories"].append(loss_history)
-            all_scores.append(trainer.test(self.model, test_set))
+                self._cv_diagnostics["loss_histories"].append(loss_history)
+            all_scores.append(trainer.test(model, test_set))
             # store the most recent train figure for diagnostics
             self.train_figure = trainer.train_figure 
-            
+
         # average individual test scores to get the cross-val metric for chosen
         # hyperparameters
         cv_score = {
@@ -241,9 +241,9 @@ def run_crossval(self, dataset):
         return cv_score
 
     @property
-    def diagnostics(self):
+    def cv_diagnostics(self):
         """Dict containing diagnostics of the cross-validation loop"""
-        return self._diagnostics
+        return self._cv_diagnostics
 
 
 class RandomCellSplitGridded:

From e4c121b4e208826dee02f38eadb9840f77a7c7fa Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Fri, 24 Feb 2023 11:12:51 -0500
Subject: [PATCH 018/109] fold deproject_vis into get_1d_vis_fit

---
 src/mpol/geometry.py | 67 --------------------------------------------
 1 file changed, 67 deletions(-)

diff --git a/src/mpol/geometry.py b/src/mpol/geometry.py
index 3fee3c08..add68253 100644
--- a/src/mpol/geometry.py
+++ b/src/mpol/geometry.py
@@ -129,70 +129,3 @@ def observer_to_flat(X, Y, omega=None, incl=None, Omega=None):
         y = y1
 
     return x, y
-
-
-# TODO: check giving correct result 
-def deproject_vis(u, v, V, weights, source_geom, inverse=False, rescale_flux=True): 
-    r"""
-    Deproject (or reproject) visibilities (and optionally rescale their 
-    amplitudes and weights, according to the source inclination) using routines 
-    in `frank <https://github.com/discsim/frank>`_
-
-    Parameters
-    ----------
-    u : array of real, size = N, unit = :math:`\lambda`
-        u-points of the visibilities
-    v : array of real, size = N, unit = :math:`\lambda`
-        v-points of the visibilities
-    V : array of real, size = N, unit = Jy
-        Complex visibilites
-    weights : array of real, size = N, unit = Jy
-        Weights on the visibilities    
-    source_geom : dict
-        Dictionary with source geometry parameters. Keys:
-            "inc" : float, unit = deg
-                Inclination
-            "PA" : float, unit = deg
-                Position angle, defined east of north.  
-            "dRA" : float, unit = arcsec
-                Phase centre offset in right ascension.
-            "dDec" : float, units = arcsec
-                Phase centre offset in declination.           
-    inverse : bool, default=False
-        If True, the uv-points are reprojected rather than deprojected
-    rescale_flux : bool, default=True
-        If True, the visibility amplitudes and weights are rescaled to account 
-        for the difference between the inclined (observed) brightness and the 
-        assumed face-on brightness, assuming the emission is optically thick. 
-        The source's integrated (2D) flux is assumed to be:
-            :math:`F = \cos(i) \int_r^{r=R}{I(r) 2 \pi r dr}`.
-        No rescaling would be appropriate in the optically thin limit.
-    Returns
-    -------
-    up : array of real, size = N, unit = :math:`\lambda`
-        Corrected u-points of the visibilities
-    vp : array of real, size = N, unit = :math:`\lambda`
-        Corrected v-points of the visibilities
-    Vp : array of real, size = N, unit = Jy
-        Corrected complex visibilites
-    weights_scaled : array of real, size = N, unit = Jy
-        Rescaled weights on the visibilities        
-
-    """
-    from frank.geometry import FixedGeometry
-
-    geom = FixedGeometry(**source_geom)
-
-    weights_scaled = weights * 1
-
-    if inverse:
-        up, vp, Vp = geom.apply_correction(u, v, V)
-        if rescale_flux:
-            Vp, weights_scaled = geom.rescale_total_flux(V, weights)
-
-    else:
-        if rescale_flux:
-            Vp, weights_scaled = geom.rescale_total_flux(V, weights)
-        up, vp, Vp = geom.undo_correction(u, v, V)
-
-    return up, vp, Vp, weights_scaled

From e30b3db5a9ebfec304daa25e72959cdd2567e827 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Fri, 24 Feb 2023 11:46:00 -0500
Subject: [PATCH 019/109] get_radial_profile: clean up docstring

---
 src/mpol/geometry.py | 2 +-
 src/mpol/onedim.py   | 5 +++--
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/mpol/geometry.py b/src/mpol/geometry.py
index add68253..ddca4715 100644
--- a/src/mpol/geometry.py
+++ b/src/mpol/geometry.py
@@ -1,4 +1,4 @@
-"""The geometry package provides routines for projecting and de-projecting sky images and visibilities.
+"""The geometry package provides routines for projecting and de-projecting sky images.
 """
 
 import numpy as np
diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index f3ac7cd4..645061a8 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -50,9 +50,10 @@ def get_radial_profile(model, center=(0.0, 0.0), bins=None, chan=0):
         Offset (RA, Dec) of source in image. Postive RA offset is west of
         north. If None, the source is assumed to be at the image center pixel
     bins : array, default=None, unit=[arcsec]
-        Radial bin edges to use in calculating I(r)
+        Radial bin edges to use in calculating I(r). If None, bins will span 
+        the full image, with widths equal to the hypotenuse of the pixels
     chan : int, default=0
-        Channel of `model` to select
+        Channel of `model` to use
 
     Returns
     -------

From 6f2471b44bc40b0c834c0deebbe2e70406f8d6aa Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Fri, 24 Feb 2023 11:46:25 -0500
Subject: [PATCH 020/109] get_radial_profile: convert to rad

---
 src/mpol/onedim.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 645061a8..1034bb9c 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -63,7 +63,11 @@ def get_radial_profile(model, center=(0.0, 0.0), bins=None, chan=0):
         Azimuthally averaged pixel brightness at `rs`
 
     """
-    # observer_to_flat(X, Y, omega, incl, Omega) # TODO
+
+    # convert to radian
+    incl = geom["incl"] * np.pi/180
+    Omega = geom["Omega"] * np.pi/180
+    omega = geom["omega"] * np.pi/180
 
     skycube = torch2npy(model.icube.sky_cube)[chan]
 

From df1ecdc43f35dd0ae7dbc148e1dd1c9a6af5c80b Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Fri, 24 Feb 2023 11:47:56 -0500
Subject: [PATCH 021/109] get_radial_profile: add 'geom' arg

---
 src/mpol/onedim.py | 22 ++++++++++++++++------
 1 file changed, 16 insertions(+), 6 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 1034bb9c..01c8cbc1 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -37,18 +37,26 @@ def get_1d_vis_fit(model, chan=0):
     return q, Vmod
 
 
-def get_radial_profile(model, center=(0.0, 0.0), bins=None, chan=0):
+def get_radial_profile(model, geom, bins=None, chan=0):
     r"""
     Obtain a 1D (radial) brightness profile I(r) from an MPoL model.
 
     Parameters
     ----------
-    # TODO
     model : `torch.nn.Module` object
         Instance of the `mpol.precomposed.SimpleNet` class
-    center : 2-tuple of float, default=(0.0, 0.0), unit=[arcsec]
-        Offset (RA, Dec) of source in image. Postive RA offset is west of
-        north. If None, the source is assumed to be at the image center pixel
+    geom : dict 
+        Dictionary of source geometry. Keys:
+            "incl" : float, unit=[deg]
+                Inclination 
+            "Omega" : float, unit=[deg]
+                Position angle of the ascending node # TODO: convention?
+            "omega" : float, unit=[deg]
+                Argument of periastron
+            "dRA" : float, unit=[arcsec]
+                Phase center offset in right ascension. Positive is west of north. # TODO: convention?
+            "dDec" : float, unit=[arcsec]
+                Phase center offset in declination.
     bins : array, default=None, unit=[arcsec]
         Radial bin edges to use in calculating I(r). If None, bins will span 
         the full image, with widths equal to the hypotenuse of the pixels
@@ -69,12 +77,14 @@ def get_radial_profile(model, center=(0.0, 0.0), bins=None, chan=0):
     Omega = geom["Omega"] * np.pi/180
     omega = geom["omega"] * np.pi/180
 
+    # model pixel values
     skycube = torch2npy(model.icube.sky_cube)[chan]
 
     # Cartesian pixel coordinates [arcsec]
     xx, yy = model.coords.sky_x_centers_2D, model.coords.sky_y_centers_2D
     # shift image center
-    xshift, yshift = xx - center[0], yy - center[1]
+    xshift, yshift = xx - geom["dRA"], yy - geom["dDec"]
+
     # radial pixel coordinates
     rshift = np.hypot(xshift, yshift) 
     rshift = rshift.flatten()

From eab4579a4c868c6c2a628d3dc611a68e9aa1adb8 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Fri, 24 Feb 2023 11:49:00 -0500
Subject: [PATCH 022/109] get_radial_profile: incorporate image deprojection

---
 src/mpol/onedim.py | 14 +++++++++-----
 1 file changed, 9 insertions(+), 5 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 01c8cbc1..54fbead5 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -79,24 +79,28 @@ def get_radial_profile(model, geom, bins=None, chan=0):
 
     # model pixel values
     skycube = torch2npy(model.icube.sky_cube)[chan]
+    # TODO: scale (multiply) brightness by inclination? if so, add arg
 
     # Cartesian pixel coordinates [arcsec]
     xx, yy = model.coords.sky_x_centers_2D, model.coords.sky_y_centers_2D
     # shift image center
     xshift, yshift = xx - geom["dRA"], yy - geom["dDec"]
 
+    # deproject and rotate image 
+    xdep, ydep = observer_to_flat(xshift, yshift, omega=omega, incl=incl, Omega=Omega) # TODO: omega
+
     # radial pixel coordinates
-    rshift = np.hypot(xshift, yshift) 
-    rshift = rshift.flatten()
+    rr = np.hypot(xdep, ydep) 
+    rr = rr.flatten()
 
     if bins is None:
         step = np.hypot(model.coords.cell_size, model.coords.cell_size)
-        bins = np.arange(0.0, max(rshift), step)
+        bins = np.arange(0.0, max(rr), step)
 
     # get number of points in each radial bin
-    bin_counts, bin_edges = np.histogram(a=rshift, bins=bins, weights=None)
+    bin_counts, bin_edges = np.histogram(a=rr, bins=bins, weights=None)
     # get radial brightness
-    Is, _ = np.histogram(a=rshift, bins=bins, weights=skycube.flatten())
+    Is, _ = np.histogram(a=rr, bins=bins, weights=skycube.flatten())
     Is /= bin_counts
     
     bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

From 5c1f0e33b4bac3c22c6419e036eecff058a01b1d Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Fri, 24 Feb 2023 11:49:58 -0500
Subject: [PATCH 023/109] get_1d_vis_fit: update args, docstring

---
 src/mpol/onedim.py | 23 +++++++++++++++++++++--
 1 file changed, 21 insertions(+), 2 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 54fbead5..4e06e379 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -5,7 +5,7 @@
 from mpol.utils import torch2npy
 
 
-def get_1d_vis_fit(model, chan=0):
+def get_1d_vis_fit(model, geom, bins=None, rescale_flux=True, chan=0):
     r"""
     Obtain the 1D (radial) visibility model V(q) corresponding to a 2D MPoL 
     image-domain model. 
@@ -14,6 +14,25 @@ def get_1d_vis_fit(model, chan=0):
     ----------
     model : `torch.nn.Module` object
         Instance of the `mpol.precomposed.SimpleNet` class
+    geom : dict 
+        Dictionary of source geometry. Keys:
+            "incl" : float, unit=[deg]
+                Inclination 
+            "Omega" : float, unit=[deg]
+                Position angle of the ascending node # TODO: convention?
+            "omega" : float, unit=[deg]
+                Argument of periastron
+            "dRA" : float, unit=[arcsec]
+                Phase center offset in right ascension. Positive is west of north. # TODO: convention?
+            "dDec" : float, unit=[arcsec]
+                Phase center offset in declination
+    rescale_flux : bool, default=True # TODO
+        If True, the visibility amplitudes and weights are rescaled to account 
+        for the difference between the inclined (observed) brightness and the 
+        assumed face-on brightness, assuming the emission is optically thick. 
+        The source's integrated (2D) flux is assumed to be:
+            :math:`F = \cos(i) \int_r^{r=R}{I(r) 2 \pi r dr}`.
+        No rescaling would be appropriate in the optically thin limit.                
     chan : int, default=0
         Channel of `model` to select
 
@@ -26,7 +45,7 @@ def get_1d_vis_fit(model, chan=0):
 
     Notes
     -----
-    This routine requires the `frank <https://github.com/discsim/frank>`_ package # TODO
+    This routine requires the `frank <https://github.com/discsim/frank>`_ package
     """
     # from mpol.geometry import deproject_vis # TODO
     # deproject_vis(u, v, V, weights, source_geom, inverse, rescale_flux) # TODO

From 14c92391a1513f8e37fc1b088af996544363975e Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Fri, 24 Feb 2023 11:50:30 -0500
Subject: [PATCH 024/109] get_1d_vis_fit: get model u, v, V

---
 src/mpol/onedim.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 4e06e379..66f346d5 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -47,11 +47,12 @@ def get_1d_vis_fit(model, geom, bins=None, rescale_flux=True, chan=0):
     -----
     This routine requires the `frank <https://github.com/discsim/frank>`_ package
     """
-    # from mpol.geometry import deproject_vis # TODO
-    # deproject_vis(u, v, V, weights, source_geom, inverse, rescale_flux) # TODO
+    # model visibility amplitudes
+    Vmod = torch2npy(model.fcube.ground_vis)[chan] # TODO: or is it model.fcube.vis?
+
+    # model (u,v) coordinates [k\lambda]
+    uu, vv = model.coords.sky_u_centers_2D, model.coords.sky_v_centers_2D
 
-    q = model.coords.packed_q_centers_2D.ravel()
-    Vmod = model.fcube.vis.detach()[chan].ravel()
 
     return q, Vmod
 

From 287e6c4a73ae532193f9613087b10c73e6712d4d Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Fri, 24 Feb 2023 11:50:47 -0500
Subject: [PATCH 025/109] get_1d_vis_fit: deproject vis with frank

---
 src/mpol/onedim.py | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 66f346d5..4f49eb49 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -47,12 +47,22 @@ def get_1d_vis_fit(model, geom, bins=None, rescale_flux=True, chan=0):
     -----
     This routine requires the `frank <https://github.com/discsim/frank>`_ package
     """
+    from frank.geometry import FixedGeometry
+    geom_frank = FixedGeometry(geom["incl"], geom["Omega"], geom["dRA"], geom["dDec"])    
+
     # model visibility amplitudes
     Vmod = torch2npy(model.fcube.ground_vis)[chan] # TODO: or is it model.fcube.vis?
 
     # model (u,v) coordinates [k\lambda]
     uu, vv = model.coords.sky_u_centers_2D, model.coords.sky_v_centers_2D
 
+    # phase-shift the model visibilities and deproject the model (u,v) points
+    up, vp, Vp = geom_frank.apply_correction(uu.ravel() * 1e3, vv.ravel() * 1e3, Vmod.ravel())
+    # if rescale_flux: # TODO: be consistent w/ get_radial_profile
+    #     Vp, weights_scaled = geom_frank.rescale_total_flux(Vp, weights)
+    # convert back to [k\lambda]
+    up /= 1e3
+    vp /= 1e3
 
     return q, Vmod
 

From 985fd9a91c415dc75c54c8050d79615461690ba2 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Fri, 24 Feb 2023 11:51:06 -0500
Subject: [PATCH 026/109] get_1d_vis_fit: get radial vis profile

---
 src/mpol/onedim.py | 17 ++++++++++++++++-
 1 file changed, 16 insertions(+), 1 deletion(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 4f49eb49..b26b29e5 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -64,7 +64,22 @@ def get_1d_vis_fit(model, geom, bins=None, rescale_flux=True, chan=0):
     up /= 1e3
     vp /= 1e3
 
-    return q, Vmod
+    # baselines
+    qq = np.hypot(up, vp) 
+
+    if bins is None:
+        step = np.hypot(model.coords.du, model.coords.dv)
+        bins = np.arange(0.0, max(qq), step)
+
+    # get number of points in each radial bin
+    bin_counts, bin_edges = np.histogram(a=qq, bins=bins, weights=None)
+    # get radial brightness
+    Vs, _ = np.histogram(a=qq, bins=bins, weights=Vp)
+    Vs /= bin_counts
+    
+    bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
+
+    return bin_centers, Vs
 
 
 def get_radial_profile(model, geom, bins=None, chan=0):

From 6a47a0ea28cb224434610c2d8175a505a2673913 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Fri, 24 Feb 2023 11:56:24 -0500
Subject: [PATCH 027/109] get_radial_profile: simplify

---
 src/mpol/onedim.py | 16 +++++++---------
 1 file changed, 7 insertions(+), 9 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index b26b29e5..d6a42669 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -47,8 +47,6 @@ def get_1d_vis_fit(model, geom, bins=None, rescale_flux=True, chan=0):
     -----
     This routine requires the `frank <https://github.com/discsim/frank>`_ package
     """
-    from frank.geometry import FixedGeometry
-    geom_frank = FixedGeometry(geom["incl"], geom["Omega"], geom["dRA"], geom["dDec"])    
 
     # model visibility amplitudes
     Vmod = torch2npy(model.fcube.ground_vis)[chan] # TODO: or is it model.fcube.vis?
@@ -56,6 +54,8 @@ def get_1d_vis_fit(model, geom, bins=None, rescale_flux=True, chan=0):
     # model (u,v) coordinates [k\lambda]
     uu, vv = model.coords.sky_u_centers_2D, model.coords.sky_v_centers_2D
 
+    from frank.geometry import FixedGeometry
+    geom_frank = FixedGeometry(geom["incl"], geom["Omega"], geom["dRA"], geom["dDec"])    
     # phase-shift the model visibilities and deproject the model (u,v) points
     up, vp, Vp = geom_frank.apply_correction(uu.ravel() * 1e3, vv.ravel() * 1e3, Vmod.ravel())
     # if rescale_flux: # TODO: be consistent w/ get_radial_profile
@@ -73,7 +73,7 @@ def get_1d_vis_fit(model, geom, bins=None, rescale_flux=True, chan=0):
 
     # get number of points in each radial bin
     bin_counts, bin_edges = np.histogram(a=qq, bins=bins, weights=None)
-    # get radial brightness
+    # get radial vis amplitude
     Vs, _ = np.histogram(a=qq, bins=bins, weights=Vp)
     Vs /= bin_counts
     
@@ -117,11 +117,6 @@ def get_radial_profile(model, geom, bins=None, chan=0):
 
     """
 
-    # convert to radian
-    incl = geom["incl"] * np.pi/180
-    Omega = geom["Omega"] * np.pi/180
-    omega = geom["omega"] * np.pi/180
-
     # model pixel values
     skycube = torch2npy(model.icube.sky_cube)[chan]
     # TODO: scale (multiply) brightness by inclination? if so, add arg
@@ -132,7 +127,10 @@ def get_radial_profile(model, geom, bins=None, chan=0):
     xshift, yshift = xx - geom["dRA"], yy - geom["dDec"]
 
     # deproject and rotate image 
-    xdep, ydep = observer_to_flat(xshift, yshift, omega=omega, incl=incl, Omega=Omega) # TODO: omega
+    xdep, ydep = observer_to_flat(xshift, yshift,
+        omega=geom["omega"] * np.pi / 180, # TODO: omega
+        incl=geom["incl"] * np.pi / 180,
+        Omega=geom["Omega"] * np.pi / 180)
 
     # radial pixel coordinates
     rr = np.hypot(xdep, ydep) 

From ed65449eda4781d919e2b4b44e62b2f2c99202af Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Tue, 28 Feb 2023 00:23:37 -0500
Subject: [PATCH 028/109] get_radial_profile: simplify

---
 src/mpol/onedim.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index d6a42669..d126fd9c 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -133,8 +133,7 @@ def get_radial_profile(model, geom, bins=None, chan=0):
         Omega=geom["Omega"] * np.pi / 180)
 
     # radial pixel coordinates
-    rr = np.hypot(xdep, ydep) 
-    rr = rr.flatten()
+    rr = np.ravel(np.hypot(xdep, ydep))
 
     if bins is None:
         step = np.hypot(model.coords.cell_size, model.coords.cell_size)
@@ -143,7 +142,7 @@ def get_radial_profile(model, geom, bins=None, chan=0):
     # get number of points in each radial bin
     bin_counts, bin_edges = np.histogram(a=rr, bins=bins, weights=None)
     # get radial brightness
-    Is, _ = np.histogram(a=rr, bins=bins, weights=skycube.flatten())
+    Is, _ = np.histogram(a=rr, bins=bins, weights=np.ravel(skycube))
     Is /= bin_counts
     
     bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

From 8ed956349babbda19da723cd584d1ce772919fb8 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Mon, 13 Mar 2023 16:55:57 -0400
Subject: [PATCH 029/109] get_1d_vis_fit debugged

---
 src/mpol/onedim.py | 37 +++++++++++++++++++++++++------------
 1 file changed, 25 insertions(+), 12 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index d126fd9c..85f2d0fc 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -1,11 +1,10 @@
 import numpy as np 
 
-from mpol.fourier import NuFFT
 from mpol.geometry import observer_to_flat
 from mpol.utils import torch2npy
 
 
-def get_1d_vis_fit(model, geom, bins=None, rescale_flux=True, chan=0):
+def get_1d_vis_fit(model, geom, rescale_flux=True, bins=None, chan=0):
     r"""
     Obtain the 1D (radial) visibility model V(q) corresponding to a 2D MPoL 
     image-domain model. 
@@ -26,13 +25,17 @@ def get_1d_vis_fit(model, geom, bins=None, rescale_flux=True, chan=0):
                 Phase center offset in right ascension. Positive is west of north. # TODO: convention?
             "dDec" : float, unit=[arcsec]
                 Phase center offset in declination
-    rescale_flux : bool, default=True # TODO
+    rescale_flux : bool, default=True 
         If True, the visibility amplitudes and weights are rescaled to account 
         for the difference between the inclined (observed) brightness and the 
         assumed face-on brightness, assuming the emission is optically thick. 
         The source's integrated (2D) flux is assumed to be:
             :math:`F = \cos(i) \int_r^{r=R}{I(r) 2 \pi r dr}`.
-        No rescaling would be appropriate in the optically thin limit.                
+        No rescaling would be appropriate in the optically thin limit. 
+    bins : array, default=None, unit=[k\lambda]
+        Baseline bin edges to use in calculating V(q). If None, bins will span 
+        the model baseline distribution, with widths equal to the hypotenuse of 
+        the (u, v) coordinates                       
     chan : int, default=0
         Channel of `model` to select
 
@@ -47,19 +50,20 @@ def get_1d_vis_fit(model, geom, bins=None, rescale_flux=True, chan=0):
     -----
     This routine requires the `frank <https://github.com/discsim/frank>`_ package
     """
-
     # model visibility amplitudes
-    Vmod = torch2npy(model.fcube.ground_vis)[chan] # TODO: or is it model.fcube.vis?
+    Vmod = torch2npy(model.fcube.ground_vis)[chan]
 
     # model (u,v) coordinates [k\lambda]
     uu, vv = model.coords.sky_u_centers_2D, model.coords.sky_v_centers_2D
 
     from frank.geometry import FixedGeometry
-    geom_frank = FixedGeometry(geom["incl"], geom["Omega"], geom["dRA"], geom["dDec"])    
+    geom_frank = FixedGeometry(geom["incl"], geom["Omega"], geom["dRA"], geom["dDec"])   # TODO: right signs?
     # phase-shift the model visibilities and deproject the model (u,v) points
     up, vp, Vp = geom_frank.apply_correction(uu.ravel() * 1e3, vv.ravel() * 1e3, Vmod.ravel())
-    # if rescale_flux: # TODO: be consistent w/ get_radial_profile
-    #     Vp, weights_scaled = geom_frank.rescale_total_flux(Vp, weights)
+
+    # if the source is optically thick, rescale the deprojected Re(V)
+    if rescale_flux: 
+        Vp.real /= np.cos(geom["incl"] * np.pi / 180)
     # convert back to [k\lambda]
     up /= 1e3
     vp /= 1e3
@@ -82,7 +86,7 @@ def get_1d_vis_fit(model, geom, bins=None, rescale_flux=True, chan=0):
     return bin_centers, Vs
 
 
-def get_radial_profile(model, geom, bins=None, chan=0):
+def get_radial_profile(model, geom, bins=None, rescale_flux=True, chan=0): 
     r"""
     Obtain a 1D (radial) brightness profile I(r) from an MPoL model.
 
@@ -102,6 +106,13 @@ def get_radial_profile(model, geom, bins=None, chan=0):
                 Phase center offset in right ascension. Positive is west of north. # TODO: convention?
             "dDec" : float, unit=[arcsec]
                 Phase center offset in declination.
+    rescale_flux : bool, default=True 
+        If True, the brightness values are rescaled to account for the 
+        difference between the inclined (observed) brightness and the 
+        assumed face-on brightness, assuming the emission is optically thick. 
+        The source's integrated (2D) flux is assumed to be:
+            :math:`F = \cos(i) \int_r^{r=R}{I(r) 2 \pi r dr}`.
+        No rescaling would be appropriate in the optically thin limit. 
     bins : array, default=None, unit=[arcsec]
         Radial bin edges to use in calculating I(r). If None, bins will span 
         the full image, with widths equal to the hypotenuse of the pixels
@@ -130,7 +141,9 @@ def get_radial_profile(model, geom, bins=None, chan=0):
     xdep, ydep = observer_to_flat(xshift, yshift,
         omega=geom["omega"] * np.pi / 180, # TODO: omega
         incl=geom["incl"] * np.pi / 180,
-        Omega=geom["Omega"] * np.pi / 180)
+        Omega=geom["Omega"] * np.pi / 180,
+        opt_thick=rescale_flux,
+        )
 
     # radial pixel coordinates
     rr = np.ravel(np.hypot(xdep, ydep))
@@ -144,7 +157,7 @@ def get_radial_profile(model, geom, bins=None, chan=0):
     # get radial brightness
     Is, _ = np.histogram(a=rr, bins=bins, weights=np.ravel(skycube))
     Is /= bin_counts
-    
+
     bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
 
     return bin_centers, Is

From aa7d9d9d2b2ad2a253ddb70674bad5b8f859ec45 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Mon, 13 Mar 2023 16:59:04 -0400
Subject: [PATCH 030/109] flat_to_observer, observer_to_flat: opt_thick arg

---
 src/mpol/geometry.py | 10 ++++++++--
 src/mpol/onedim.py   |  3 ---
 2 files changed, 8 insertions(+), 5 deletions(-)

diff --git a/src/mpol/geometry.py b/src/mpol/geometry.py
index ddca4715..5b80250c 100644
--- a/src/mpol/geometry.py
+++ b/src/mpol/geometry.py
@@ -5,7 +5,7 @@
 import torch
 
 
-def flat_to_observer(x, y, omega=None, incl=None, Omega=None):
+def flat_to_observer(x, y, omega=None, incl=None, Omega=None, opt_thick=True):
     """Rotate the frame to convert a point in the flat (x,y,z) frame to observer frame (X,Y,Z).
 
     It is assumed that the +Z axis points *towards* the observer. It is assumed that the model is flat in the (x,y) frame (like a flat disk), and so the operations involving ``z`` are neglected. But the model lives in 3D Cartesian space.
@@ -24,10 +24,13 @@ def flat_to_observer(x, y, omega=None, incl=None, Omega=None):
         omega (torch float tensor): A tensor representing an argument of periastron [radians] Default 0.0.
         incl (torch float tensor): A tensor representing an inclination value [radians]. Default 0.0.
         Omega (torch float tensor): A tensor representing the position angle of the ascending node in [radians]. Default 0.0
+        opt_thick (bool): Whether the source is optically thick. If True, brightness in the observer frame is scaled (multiplied) by cos(incl). If False, 'incl' will be ignored. Default True
 
     Returns:
         Two tensors representing ``(X, Y)`` in the observer frame.
     """
+    if not opt_thick:
+        incl = None 
 
     # Rotation matrices result in a *clockwise* rotation of the axes, as defined using the righthand rule.
     # For example, looking down the z-axis, a positive angle will rotate the x,y axes clockwise.
@@ -69,7 +72,7 @@ def flat_to_observer(x, y, omega=None, incl=None, Omega=None):
     return X, Y
 
 
-def observer_to_flat(X, Y, omega=None, incl=None, Omega=None):
+def observer_to_flat(X, Y, omega=None, incl=None, Omega=None, opt_thick=True):
     """Rotate the frame to convert a point in the observer frame (X,Y,Z) to the flat (x,y,z) frame.
 
     It is assumed that the +Z axis points *towards* the observer. The rotation operations are the inverse of the :func:`~mpol.geometry.flat_to_observer` operations.
@@ -88,10 +91,13 @@ def observer_to_flat(X, Y, omega=None, incl=None, Omega=None):
         omega (torch float tensor): A tensor representing an argument of periastron [radians] Default 0.0.
         incl (torch float tensor): A tensor representing an inclination value [radians]. Default 0.0.
         Omega (torch float tensor): A tensor representing the position angle of the ascending node in [radians]. Default 0.0
+        opt_thick (bool): Whether the source is optically thick. If True, brightness in the flat frame is scaled (divided) by cos(incl). If False, 'incl' will be ignored. Default True
 
     Returns:
         Two tensors representing ``(x, y)`` in the flat frame.
     """
+    if not opt_thick:
+        incl = None 
 
     # Rotation matrices result in a *clockwise* rotation of the axes, as defined using the righthand rule.
     # For example, looking down the z-axis, a positive angle will rotate the x,y axes clockwise.
diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 85f2d0fc..8941d46b 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -125,12 +125,9 @@ def get_radial_profile(model, geom, bins=None, rescale_flux=True, chan=0):
         Radial coordinates of image at center of `bins`
     Is : array, unit=[Jy / arcsec^2]
         Azimuthally averaged pixel brightness at `rs`
-
     """
-
     # model pixel values
     skycube = torch2npy(model.icube.sky_cube)[chan]
-    # TODO: scale (multiply) brightness by inclination? if so, add arg
 
     # Cartesian pixel coordinates [arcsec]
     xx, yy = model.coords.sky_x_centers_2D, model.coords.sky_y_centers_2D

From 35c7b9339caba24468666c4eb8a8aed800e95e5a Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Mon, 13 Mar 2023 17:31:42 -0400
Subject: [PATCH 031/109] mypy: add frank ignore

---
 mypy.ini | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/mypy.ini b/mypy.ini
index 05045688..74c14072 100644
--- a/mypy.ini
+++ b/mypy.ini
@@ -12,4 +12,7 @@ ignore_missing_imports = True
 ignore_missing_imports = True
 
 [mypy-torchkbnufft.*]
+ignore_missing_imports = True
+
+[mypy-frank.*]
 ignore_missing_imports = True
\ No newline at end of file

From beb3de1df5be4651604c9206133360817a2d7c5c Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Mon, 13 Mar 2023 18:17:05 -0400
Subject: [PATCH 032/109] revert opt_thick geometry arg

---
 src/mpol/geometry.py | 8 --------
 1 file changed, 8 deletions(-)

diff --git a/src/mpol/geometry.py b/src/mpol/geometry.py
index 5b80250c..85b37051 100644
--- a/src/mpol/geometry.py
+++ b/src/mpol/geometry.py
@@ -24,14 +24,10 @@ def flat_to_observer(x, y, omega=None, incl=None, Omega=None, opt_thick=True):
         omega (torch float tensor): A tensor representing an argument of periastron [radians] Default 0.0.
         incl (torch float tensor): A tensor representing an inclination value [radians]. Default 0.0.
         Omega (torch float tensor): A tensor representing the position angle of the ascending node in [radians]. Default 0.0
-        opt_thick (bool): Whether the source is optically thick. If True, brightness in the observer frame is scaled (multiplied) by cos(incl). If False, 'incl' will be ignored. Default True
 
     Returns:
         Two tensors representing ``(X, Y)`` in the observer frame.
     """
-    if not opt_thick:
-        incl = None 
-
     # Rotation matrices result in a *clockwise* rotation of the axes, as defined using the righthand rule.
     # For example, looking down the z-axis, a positive angle will rotate the x,y axes clockwise.
     # A vector in the coordinate system will appear as though it has been rotated counter-clockwise.
@@ -91,14 +87,10 @@ def observer_to_flat(X, Y, omega=None, incl=None, Omega=None, opt_thick=True):
         omega (torch float tensor): A tensor representing an argument of periastron [radians] Default 0.0.
         incl (torch float tensor): A tensor representing an inclination value [radians]. Default 0.0.
         Omega (torch float tensor): A tensor representing the position angle of the ascending node in [radians]. Default 0.0
-        opt_thick (bool): Whether the source is optically thick. If True, brightness in the flat frame is scaled (divided) by cos(incl). If False, 'incl' will be ignored. Default True
 
     Returns:
         Two tensors representing ``(x, y)`` in the flat frame.
     """
-    if not opt_thick:
-        incl = None 
-
     # Rotation matrices result in a *clockwise* rotation of the axes, as defined using the righthand rule.
     # For example, looking down the z-axis, a positive angle will rotate the x,y axes clockwise.
     # A vector in the coordinate system will appear as though it has been rotated counter-clockwise.

From 9fc9c2631b9888d4975c311d64eb8394d656f539 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Mon, 13 Mar 2023 18:17:46 -0400
Subject: [PATCH 033/109] typos

---
 src/mpol/geometry.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/mpol/geometry.py b/src/mpol/geometry.py
index 85b37051..b7438971 100644
--- a/src/mpol/geometry.py
+++ b/src/mpol/geometry.py
@@ -5,7 +5,7 @@
 import torch
 
 
-def flat_to_observer(x, y, omega=None, incl=None, Omega=None, opt_thick=True):
+def flat_to_observer(x, y, omega=None, incl=None, Omega=None):
     """Rotate the frame to convert a point in the flat (x,y,z) frame to observer frame (X,Y,Z).
 
     It is assumed that the +Z axis points *towards* the observer. It is assumed that the model is flat in the (x,y) frame (like a flat disk), and so the operations involving ``z`` are neglected. But the model lives in 3D Cartesian space.
@@ -68,7 +68,7 @@ def flat_to_observer(x, y, omega=None, incl=None, Omega=None, opt_thick=True):
     return X, Y
 
 
-def observer_to_flat(X, Y, omega=None, incl=None, Omega=None, opt_thick=True):
+def observer_to_flat(X, Y, omega=None, incl=None, Omega=None):
     """Rotate the frame to convert a point in the observer frame (X,Y,Z) to the flat (x,y,z) frame.
 
     It is assumed that the +Z axis points *towards* the observer. The rotation operations are the inverse of the :func:`~mpol.geometry.flat_to_observer` operations.

From d9f21a27aa166d7027bade49bf768726d7b5ec93 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 22 Mar 2023 09:13:25 -0400
Subject: [PATCH 034/109] docstring spacing

---
 src/mpol/plot.py  | 2 +-
 src/mpol/utils.py | 5 ++---
 2 files changed, 3 insertions(+), 4 deletions(-)

diff --git a/src/mpol/plot.py b/src/mpol/plot.py
index 4c3a3d3a..6db00668 100644
--- a/src/mpol/plot.py
+++ b/src/mpol/plot.py
@@ -65,7 +65,7 @@ def plot_image(image, extent, cmap="inferno", norm=None, ax=None,
         Image x-axis label.
     ylab : str, default="Dec offset [arcsec]"
         Image y-axis label.
-
+        
     Returns
     -------
     im : Matplotlib imshow instance
diff --git a/src/mpol/utils.py b/src/mpol/utils.py
index ac8ced23..59e31683 100644
--- a/src/mpol/utils.py
+++ b/src/mpol/utils.py
@@ -286,7 +286,7 @@ def get_optimal_image_properties(image_width, u, v):
     For an image of desired width, determine the maximum pixel size that
     ensures Nyquist sampling of the provided spatial frequency points, and the 
     corresponding number of pixels to obtain the desired image width.
-
+    
     Parameters
     ----------
     image_width : float, unit = arcsec
@@ -294,7 +294,7 @@ def get_optimal_image_properties(image_width, u, v):
         `image_width` :math:`\times` `image_width`).
     u, v : float or array, unit = :math:`k\lambda`
         `u` and `v` spatial frequency points. 
-
+        
     Returns
     -------
     cell_size : float, unit = arcsec
@@ -302,7 +302,6 @@ def get_optimal_image_properties(image_width, u, v):
     npix : int
         Number of pixels of cell_size to equal (or slightly exceed) the image
         width (npix will be rounded up and enforced as even).
-
     Notes
     -----
     No assumption or correction is made concerning whether the spatial 

From c905d258c34d469ee60e37281e6dbe43b9a0ef40 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 30 Mar 2023 14:06:53 -0400
Subject: [PATCH 035/109] get_1d_vis_fit: intake vis and coords, not model

---
 src/mpol/onedim.py | 34 ++++++++++++++++------------------
 1 file changed, 16 insertions(+), 18 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 8941d46b..6f34d724 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -4,25 +4,27 @@
 from mpol.utils import torch2npy
 
 
-def get_1d_vis_fit(model, geom, rescale_flux=True, bins=None, chan=0):
+def get_1d_vis_fit(V, coords, geom, rescale_flux=True, bins=None):
     r"""
     Obtain the 1D (radial) visibility model V(q) corresponding to a 2D MPoL 
     image-domain model. 
 
     Parameters
     ----------
-    model : `torch.nn.Module` object
-        Instance of the `mpol.precomposed.SimpleNet` class
+    V : array
+        2D visibility amplitudes
+    coords : `mpol.coordinates.GridCoords` object
+        Instance of the `mpol.coordinates.GridCoords` class
     geom : dict 
         Dictionary of source geometry. Keys:
             "incl" : float, unit=[deg]
                 Inclination 
             "Omega" : float, unit=[deg]
-                Position angle of the ascending node # TODO: convention?
+                Position angle of the ascending node 
             "omega" : float, unit=[deg]
                 Argument of periastron
             "dRA" : float, unit=[arcsec]
-                Phase center offset in right ascension. Positive is west of north. # TODO: convention?
+                Phase center offset in right ascension. Positive is west of north. 
             "dDec" : float, unit=[arcsec]
                 Phase center offset in declination
     rescale_flux : bool, default=True 
@@ -35,49 +37,45 @@ def get_1d_vis_fit(model, geom, rescale_flux=True, bins=None, chan=0):
     bins : array, default=None, unit=[k\lambda]
         Baseline bin edges to use in calculating V(q). If None, bins will span 
         the model baseline distribution, with widths equal to the hypotenuse of 
-        the (u, v) coordinates                       
-    chan : int, default=0
-        Channel of `model` to select
+        the (u, v) coordinates
 
     Returns
     -------
     q : array, unit=:math:[`k\lambda`]
         Baselines corresponding to `u` and `v`
-    Vmod : array, unit=[Jy] 
+    Vs : array, unit=[Jy] 
         Visibility amplitudes at `q`
 
     Notes
     -----
     This routine requires the `frank <https://github.com/discsim/frank>`_ package
     """
-    # model visibility amplitudes
-    Vmod = torch2npy(model.fcube.ground_vis)[chan]
 
     # model (u,v) coordinates [k\lambda]
-    uu, vv = model.coords.sky_u_centers_2D, model.coords.sky_v_centers_2D
+    uu, vv = coords.sky_u_centers_2D, coords.sky_v_centers_2D
 
     from frank.geometry import FixedGeometry
-    geom_frank = FixedGeometry(geom["incl"], geom["Omega"], geom["dRA"], geom["dDec"])   # TODO: right signs?
+    geom_frank = FixedGeometry(geom["incl"], geom["Omega"], geom["dRA"], geom["dDec"])  # TODO: signs
     # phase-shift the model visibilities and deproject the model (u,v) points
-    up, vp, Vp = geom_frank.apply_correction(uu.ravel() * 1e3, vv.ravel() * 1e3, Vmod.ravel())
+    up, vp, Vp = geom_frank.apply_correction(uu.ravel() * 1e3, vv.ravel() * 1e3, V.ravel())
 
     # if the source is optically thick, rescale the deprojected Re(V)
     if rescale_flux: 
         Vp.real /= np.cos(geom["incl"] * np.pi / 180)
+
     # convert back to [k\lambda]
     up /= 1e3
     vp /= 1e3
-
-    # baselines
+    
     qq = np.hypot(up, vp) 
 
     if bins is None:
-        step = np.hypot(model.coords.du, model.coords.dv)
+        step = np.hypot(coords.du, coords.dv)
         bins = np.arange(0.0, max(qq), step)
 
     # get number of points in each radial bin
     bin_counts, bin_edges = np.histogram(a=qq, bins=bins, weights=None)
-    # get radial vis amplitude
+    # get radial vis amplitudes
     Vs, _ = np.histogram(a=qq, bins=bins, weights=Vp)
     Vs /= bin_counts
     

From 2ebd17e8476081e8728ffd9344fecd9e1f8ab137 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 30 Mar 2023 14:08:00 -0400
Subject: [PATCH 036/109] get_1d_vis_profile: rename

---
 src/mpol/onedim.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 6f34d724..03242122 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -4,7 +4,7 @@
 from mpol.utils import torch2npy
 
 
-def get_1d_vis_fit(V, coords, geom, rescale_flux=True, bins=None):
+def get_1d_vis_profile(V, coords, geom, rescale_flux=True, bins=None):
     r"""
     Obtain the 1D (radial) visibility model V(q) corresponding to a 2D MPoL 
     image-domain model. 

From 349d9e7dc5d65caa0fcccd34eff1483acac3f442 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 30 Mar 2023 14:08:21 -0400
Subject: [PATCH 037/109] get_1d_image_profile: rename

---
 src/mpol/onedim.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 03242122..b311f4e2 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -66,7 +66,7 @@ def get_1d_vis_profile(V, coords, geom, rescale_flux=True, bins=None):
     # convert back to [k\lambda]
     up /= 1e3
     vp /= 1e3
-    
+
     qq = np.hypot(up, vp) 
 
     if bins is None:
@@ -84,7 +84,7 @@ def get_1d_vis_profile(V, coords, geom, rescale_flux=True, bins=None):
     return bin_centers, Vs
 
 
-def get_radial_profile(model, geom, bins=None, rescale_flux=True, chan=0): 
+def get_1d_image_profile(image, coords, geom, bins=None, rescale_flux=True): 
     r"""
     Obtain a 1D (radial) brightness profile I(r) from an MPoL model.
 

From 9405943fb5d224e2ba06bf2a92572e9463e06016 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 30 Mar 2023 14:08:51 -0400
Subject: [PATCH 038/109]  get_1d_image_profile: intake image and coords

---
 src/mpol/onedim.py | 29 +++++++++++++----------------
 1 file changed, 13 insertions(+), 16 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index b311f4e2..fa4d25f9 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -86,22 +86,24 @@ def get_1d_vis_profile(V, coords, geom, rescale_flux=True, bins=None):
 
 def get_1d_image_profile(image, coords, geom, bins=None, rescale_flux=True): 
     r"""
-    Obtain a 1D (radial) brightness profile I(r) from an MPoL model.
+    Obtain a 1D (radial) brightness profile I(r) from an image.
 
     Parameters
     ----------
-    model : `torch.nn.Module` object
-        Instance of the `mpol.precomposed.SimpleNet` class
+    image : array
+        2D image array 
+    coords : `mpol.coordinates.GridCoords` object
+        Instance of the `mpol.coordinates.GridCoords` class
     geom : dict 
         Dictionary of source geometry. Keys:
             "incl" : float, unit=[deg]
                 Inclination 
             "Omega" : float, unit=[deg]
-                Position angle of the ascending node # TODO: convention?
+                Position angle of the ascending node 
             "omega" : float, unit=[deg]
                 Argument of periastron
             "dRA" : float, unit=[arcsec]
-                Phase center offset in right ascension. Positive is west of north. # TODO: convention?
+                Phase center offset in right ascension. Positive is west of north.
             "dDec" : float, unit=[arcsec]
                 Phase center offset in declination.
     rescale_flux : bool, default=True 
@@ -114,43 +116,38 @@ def get_1d_image_profile(image, coords, geom, bins=None, rescale_flux=True):
     bins : array, default=None, unit=[arcsec]
         Radial bin edges to use in calculating I(r). If None, bins will span 
         the full image, with widths equal to the hypotenuse of the pixels
-    chan : int, default=0
-        Channel of `model` to use
 
     Returns
     -------
     bin_centers : array, unit=[arcsec]
         Radial coordinates of image at center of `bins`
-    Is : array, unit=[Jy / arcsec^2]
+    Is : array, unit=[Jy / arcsec^2] (if `image` has these units)
         Azimuthally averaged pixel brightness at `rs`
     """
-    # model pixel values
-    skycube = torch2npy(model.icube.sky_cube)[chan]
 
     # Cartesian pixel coordinates [arcsec]
-    xx, yy = model.coords.sky_x_centers_2D, model.coords.sky_y_centers_2D
+    xx, yy = coords.sky_x_centers_2D, coords.sky_y_centers_2D
     # shift image center
     xshift, yshift = xx - geom["dRA"], yy - geom["dDec"]
 
     # deproject and rotate image 
     xdep, ydep = observer_to_flat(xshift, yshift,
-        omega=geom["omega"] * np.pi / 180, # TODO: omega
+        omega=geom["omega"] * np.pi / 180, 
         incl=geom["incl"] * np.pi / 180,
         Omega=geom["Omega"] * np.pi / 180,
-        opt_thick=rescale_flux,
         )
-
+    
     # radial pixel coordinates
     rr = np.ravel(np.hypot(xdep, ydep))
 
     if bins is None:
-        step = np.hypot(model.coords.cell_size, model.coords.cell_size)
+        step = np.hypot(coords.cell_size, coords.cell_size)
         bins = np.arange(0.0, max(rr), step)
 
     # get number of points in each radial bin
     bin_counts, bin_edges = np.histogram(a=rr, bins=bins, weights=None)
     # get radial brightness
-    Is, _ = np.histogram(a=rr, bins=bins, weights=np.ravel(skycube))
+    Is, _ = np.histogram(a=rr, bins=bins, weights=np.ravel(image))
     Is /= bin_counts
 
     bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

From 7a36f259f1aeb95a2dae2d0bf483769dc291f1da Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 30 Mar 2023 14:09:12 -0400
Subject: [PATCH 039/109] get_1d_image_profile: add flux rescaling

---
 src/mpol/onedim.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index fa4d25f9..cb609a75 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -150,6 +150,9 @@ def get_1d_image_profile(image, coords, geom, bins=None, rescale_flux=True):
     Is, _ = np.histogram(a=rr, bins=bins, weights=np.ravel(image))
     Is /= bin_counts
 
+    if rescale_flux: # TODO
+        Is *= np.cos(geom["incl"] * np.pi / 180)
+
     bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
 
     return bin_centers, Is

From ac01276edb1b9ff106865652d2099e1ed988b719 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 22 Jun 2023 13:12:39 -0400
Subject: [PATCH 040/109] get_1d_vis_profile: explicit deproject

---
 src/mpol/onedim.py | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index cb609a75..f6210896 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -54,10 +54,11 @@ def get_1d_vis_profile(V, coords, geom, rescale_flux=True, bins=None):
     # model (u,v) coordinates [k\lambda]
     uu, vv = coords.sky_u_centers_2D, coords.sky_v_centers_2D
 
-    from frank.geometry import FixedGeometry
-    geom_frank = FixedGeometry(geom["incl"], geom["Omega"], geom["dRA"], geom["dDec"])  # TODO: signs
-    # phase-shift the model visibilities and deproject the model (u,v) points
-    up, vp, Vp = geom_frank.apply_correction(uu.ravel() * 1e3, vv.ravel() * 1e3, V.ravel())
+    from frank.geometry import apply_phase_shift, deproject
+    # phase-shift the model visibilities
+    Vp = apply_phase_shift(uu.ravel() * 1e3, vv.ravel() * 1e3, V.ravel(), geom["dRA"], geom["dDec"], inverse=True)
+    # deproject the model (u,v) points
+    up, vp, _ = deproject(uu.ravel() * 1e3, vv.ravel() * 1e3, geom["incl"], geom["Omega"])
 
     # if the source is optically thick, rescale the deprojected Re(V)
     if rescale_flux: 
@@ -66,16 +67,15 @@ def get_1d_vis_profile(V, coords, geom, rescale_flux=True, bins=None):
     # convert back to [k\lambda]
     up /= 1e3
     vp /= 1e3
-
     qq = np.hypot(up, vp) 
 
     if bins is None:
         step = np.hypot(coords.du, coords.dv)
         bins = np.arange(0.0, max(qq), step)
 
-    # get number of points in each radial bin
+    # number of points in radial bins
     bin_counts, bin_edges = np.histogram(a=qq, bins=bins, weights=None)
-    # get radial vis amplitudes
+    # V amplitudes in radial bins
     Vs, _ = np.histogram(a=qq, bins=bins, weights=Vp)
     Vs /= bin_counts
     
@@ -144,9 +144,9 @@ def get_1d_image_profile(image, coords, geom, bins=None, rescale_flux=True):
         step = np.hypot(coords.cell_size, coords.cell_size)
         bins = np.arange(0.0, max(rr), step)
 
-    # get number of points in each radial bin
+    # number of points in radial bins
     bin_counts, bin_edges = np.histogram(a=rr, bins=bins, weights=None)
-    # get radial brightness
+    # brightness in radial bins
     Is, _ = np.histogram(a=rr, bins=bins, weights=np.ravel(image))
     Is /= bin_counts
 

From 6a3a24f135a683dc15832b6dfcfcf6484eecc421 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Tue, 17 Oct 2023 10:52:44 -0400
Subject: [PATCH 041/109] blank space

---
 src/mpol/plot.py  | 2 +-
 src/mpol/utils.py | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/mpol/plot.py b/src/mpol/plot.py
index 6db00668..4c3a3d3a 100644
--- a/src/mpol/plot.py
+++ b/src/mpol/plot.py
@@ -65,7 +65,7 @@ def plot_image(image, extent, cmap="inferno", norm=None, ax=None,
         Image x-axis label.
     ylab : str, default="Dec offset [arcsec]"
         Image y-axis label.
-        
+
     Returns
     -------
     im : Matplotlib imshow instance
diff --git a/src/mpol/utils.py b/src/mpol/utils.py
index 59e31683..92096b19 100644
--- a/src/mpol/utils.py
+++ b/src/mpol/utils.py
@@ -286,7 +286,7 @@ def get_optimal_image_properties(image_width, u, v):
     For an image of desired width, determine the maximum pixel size that
     ensures Nyquist sampling of the provided spatial frequency points, and the 
     corresponding number of pixels to obtain the desired image width.
-    
+
     Parameters
     ----------
     image_width : float, unit = arcsec
@@ -294,7 +294,7 @@ def get_optimal_image_properties(image_width, u, v):
         `image_width` :math:`\times` `image_width`).
     u, v : float or array, unit = :math:`k\lambda`
         `u` and `v` spatial frequency points. 
-        
+
     Returns
     -------
     cell_size : float, unit = arcsec

From 6d06d9e09584acde7e432e8470460191e4f8e20b Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Tue, 17 Oct 2023 10:58:03 -0400
Subject: [PATCH 042/109] create 1d tests file

---
 test/onedim_test.py | 4 ++++
 1 file changed, 4 insertions(+)
 create mode 100644 test/onedim_test.py

diff --git a/test/onedim_test.py b/test/onedim_test.py
new file mode 100644
index 00000000..ed5d1388
--- /dev/null
+++ b/test/onedim_test.py
@@ -0,0 +1,4 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+from mpol.onedim import radialI, radialV

From cf1c04d65a7994ebbe8feec73cae6043d7ae5e0d Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Tue, 17 Oct 2023 18:41:20 -0400
Subject: [PATCH 043/109] rename 1d vis profile func

---
 src/mpol/onedim.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index f6210896..f682c3fa 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -3,8 +3,7 @@
 from mpol.geometry import observer_to_flat
 from mpol.utils import torch2npy
 
-
-def get_1d_vis_profile(V, coords, geom, rescale_flux=True, bins=None):
+def radialV(V, coords, geom, rescale_flux, bins=None):
     r"""
     Obtain the 1D (radial) visibility model V(q) corresponding to a 2D MPoL 
     image-domain model. 

From e3aedf9bb66349f3b9a840fc6c5528b8af8d22c8 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Tue, 17 Oct 2023 18:41:53 -0400
Subject: [PATCH 044/109] radialV: force specify rescale_flux

---
 src/mpol/onedim.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index f682c3fa..b272d9d0 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -6,7 +6,7 @@
 def radialV(V, coords, geom, rescale_flux, bins=None):
     r"""
     Obtain the 1D (radial) visibility model V(q) corresponding to a 2D MPoL 
-    image-domain model. 
+    image. 
 
     Parameters
     ----------
@@ -26,7 +26,7 @@ def radialV(V, coords, geom, rescale_flux, bins=None):
                 Phase center offset in right ascension. Positive is west of north. 
             "dDec" : float, unit=[arcsec]
                 Phase center offset in declination
-    rescale_flux : bool, default=True 
+    rescale_flux : bool
         If True, the visibility amplitudes and weights are rescaled to account 
         for the difference between the inclined (observed) brightness and the 
         assumed face-on brightness, assuming the emission is optically thick. 
@@ -50,7 +50,7 @@ def radialV(V, coords, geom, rescale_flux, bins=None):
     This routine requires the `frank <https://github.com/discsim/frank>`_ package
     """
 
-    # model (u,v) coordinates [k\lambda]
+    # projected model (u,v) coordinates [k\lambda]
     uu, vv = coords.sky_u_centers_2D, coords.sky_v_centers_2D
 
     from frank.geometry import apply_phase_shift, deproject
@@ -59,7 +59,7 @@ def radialV(V, coords, geom, rescale_flux, bins=None):
     # deproject the model (u,v) points
     up, vp, _ = deproject(uu.ravel() * 1e3, vv.ravel() * 1e3, geom["incl"], geom["Omega"])
 
-    # if the source is optically thick, rescale the deprojected Re(V)
+    # if the source is optically thick, rescale the deprojected V(q)
     if rescale_flux: 
         Vp.real /= np.cos(geom["incl"] * np.pi / 180)
 

From 0440d60e5a2ba96e41274a0319d2af0c3c5a863f Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Tue, 17 Oct 2023 18:43:32 -0400
Subject: [PATCH 045/109] rename 1d brightness profile func

---
 src/mpol/onedim.py | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index b272d9d0..f3a88035 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -3,6 +3,13 @@
 from mpol.geometry import observer_to_flat
 from mpol.utils import torch2npy
 
+def radialI(image, coords, geom, rescale_flux, bins=None):
+    r"""
+    Obtain a 1D (radial) brightness profile I(r) from an image.
+
+    Parameters
+    ----------
+
 def radialV(V, coords, geom, rescale_flux, bins=None):
     r"""
     Obtain the 1D (radial) visibility model V(q) corresponding to a 2D MPoL 

From 9e58207fdf4da619c177d4acad103ddce3d889ed Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Tue, 17 Oct 2023 18:44:15 -0400
Subject: [PATCH 046/109] radialI: force specify rescale_flux

---
 src/mpol/onedim.py | 34 ++++++++++++++++++++++++++++++++++
 1 file changed, 34 insertions(+)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index f3a88035..5fe30c88 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -9,6 +9,40 @@ def radialI(image, coords, geom, rescale_flux, bins=None):
 
     Parameters
     ----------
+    image : array
+        2D image array 
+    coords : `mpol.coordinates.GridCoords` object
+        Instance of the `mpol.coordinates.GridCoords` class
+    geom : dict 
+        Dictionary of source geometry. Keys:
+            "incl" : float, unit=[deg]
+                Inclination 
+            "Omega" : float, unit=[deg]
+                Position angle of the ascending node 
+            "omega" : float, unit=[deg]
+                Argument of periastron
+            "dRA" : float, unit=[arcsec]
+                Phase center offset in right ascension. Positive is west of north.
+            "dDec" : float, unit=[arcsec]
+                Phase center offset in declination.
+    rescale_flux : bool
+        If True, the brightness values are rescaled to account for the 
+        difference between the inclined (observed) brightness and the 
+        assumed face-on brightness, assuming the emission is optically thick. 
+        The source's integrated (2D) flux is assumed to be:
+            :math:`F = \cos(i) \int_r^{r=R}{I(r) 2 \pi r dr}`.
+        No rescaling would be appropriate in the optically thin limit. 
+    bins : array, default=None, unit=[arcsec]
+        Radial bin edges to use in calculating I(r). If None, bins will span 
+        the full image, with widths equal to the hypotenuse of the pixels
+
+    Returns
+    -------
+    bin_centers : array, unit=[arcsec]
+        Radial coordinates of image at center of `bins`
+    Is : array, unit=[Jy / arcsec^2] (if `image` has these units)
+        Azimuthally averaged pixel brightness at `rs`
+    """
 
 def radialV(V, coords, geom, rescale_flux, bins=None):
     r"""

From 446530e57567dc9613e33c38913fca4ebedb8cab Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Tue, 17 Oct 2023 18:45:41 -0400
Subject: [PATCH 047/109] radialI: simplify

---
 src/mpol/onedim.py | 36 ++++++++++++++++++++++++++++++++++++
 1 file changed, 36 insertions(+)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 5fe30c88..39df97d4 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -44,6 +44,42 @@ def radialI(image, coords, geom, rescale_flux, bins=None):
         Azimuthally averaged pixel brightness at `rs`
     """
 
+    # projected Cartesian pixel coordinates [arcsec]
+    xx, yy = coords.sky_x_centers_2D, coords.sky_y_centers_2D
+
+    # deproject and rotate image 
+    xdep, ydep = observer_to_flat(xx, yy,
+        omega=geom["omega"] * np.pi / 180, 
+        incl=geom["incl"] * np.pi / 180,
+        Omega=geom["Omega"] * np.pi / 180,
+        )
+
+    # shift image center to source center
+    xc, yc = xdep - geom["dRA"], ydep - geom["dDec"]
+
+
+    # radial pixel coordinates
+    rr = np.ravel(np.hypot(xc, yc))
+
+    if bins is None:
+        step = np.hypot(coords.cell_size, coords.cell_size)
+        bins = np.arange(0.0, max(rr), step)
+
+    # number of points in radial bins
+    bin_counts, bin_edges = np.histogram(a=rr, bins=bins, weights=None)
+    # brightness in radial bins
+    Is, _ = np.histogram(a=rr, bins=bins, weights=np.ravel(image))
+    Is /= bin_counts
+
+    # if the source is optically thick, rescale the deprojected I(r)
+    if rescale_flux:
+        Is *= np.cos(geom["incl"] * np.pi / 180)
+
+    bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
+
+    return bin_centers, Is
+
+
 def radialV(V, coords, geom, rescale_flux, bins=None):
     r"""
     Obtain the 1D (radial) visibility model V(q) corresponding to a 2D MPoL 

From 1c5629c2e18073c39bea7cfe6aa0940285cf2906 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 19 Oct 2023 13:56:07 -0400
Subject: [PATCH 048/109] remove reimplemented func

---
 src/mpol/onedim.py | 76 +---------------------------------------------
 1 file changed, 1 insertion(+), 75 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 39df97d4..86b11d6c 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -157,78 +157,4 @@ def radialV(V, coords, geom, rescale_flux, bins=None):
     
     bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
 
-    return bin_centers, Vs
-
-
-def get_1d_image_profile(image, coords, geom, bins=None, rescale_flux=True): 
-    r"""
-    Obtain a 1D (radial) brightness profile I(r) from an image.
-
-    Parameters
-    ----------
-    image : array
-        2D image array 
-    coords : `mpol.coordinates.GridCoords` object
-        Instance of the `mpol.coordinates.GridCoords` class
-    geom : dict 
-        Dictionary of source geometry. Keys:
-            "incl" : float, unit=[deg]
-                Inclination 
-            "Omega" : float, unit=[deg]
-                Position angle of the ascending node 
-            "omega" : float, unit=[deg]
-                Argument of periastron
-            "dRA" : float, unit=[arcsec]
-                Phase center offset in right ascension. Positive is west of north.
-            "dDec" : float, unit=[arcsec]
-                Phase center offset in declination.
-    rescale_flux : bool, default=True 
-        If True, the brightness values are rescaled to account for the 
-        difference between the inclined (observed) brightness and the 
-        assumed face-on brightness, assuming the emission is optically thick. 
-        The source's integrated (2D) flux is assumed to be:
-            :math:`F = \cos(i) \int_r^{r=R}{I(r) 2 \pi r dr}`.
-        No rescaling would be appropriate in the optically thin limit. 
-    bins : array, default=None, unit=[arcsec]
-        Radial bin edges to use in calculating I(r). If None, bins will span 
-        the full image, with widths equal to the hypotenuse of the pixels
-
-    Returns
-    -------
-    bin_centers : array, unit=[arcsec]
-        Radial coordinates of image at center of `bins`
-    Is : array, unit=[Jy / arcsec^2] (if `image` has these units)
-        Azimuthally averaged pixel brightness at `rs`
-    """
-
-    # Cartesian pixel coordinates [arcsec]
-    xx, yy = coords.sky_x_centers_2D, coords.sky_y_centers_2D
-    # shift image center
-    xshift, yshift = xx - geom["dRA"], yy - geom["dDec"]
-
-    # deproject and rotate image 
-    xdep, ydep = observer_to_flat(xshift, yshift,
-        omega=geom["omega"] * np.pi / 180, 
-        incl=geom["incl"] * np.pi / 180,
-        Omega=geom["Omega"] * np.pi / 180,
-        )
-    
-    # radial pixel coordinates
-    rr = np.ravel(np.hypot(xdep, ydep))
-
-    if bins is None:
-        step = np.hypot(coords.cell_size, coords.cell_size)
-        bins = np.arange(0.0, max(rr), step)
-
-    # number of points in radial bins
-    bin_counts, bin_edges = np.histogram(a=rr, bins=bins, weights=None)
-    # brightness in radial bins
-    Is, _ = np.histogram(a=rr, bins=bins, weights=np.ravel(image))
-    Is /= bin_counts
-
-    if rescale_flux: # TODO
-        Is *= np.cos(geom["incl"] * np.pi / 180)
-
-    bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
-
-    return bin_centers, Is
+    return bin_centers, Vs
\ No newline at end of file

From 34c95b74a8f46524962ace0bb6e77998fb3c15a3 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 19 Oct 2023 13:56:41 -0400
Subject: [PATCH 049/109] radialI: add TODO checks

---
 src/mpol/onedim.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 86b11d6c..113000a6 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -57,6 +57,12 @@ def radialI(image, coords, geom, rescale_flux, bins=None):
     # shift image center to source center
     xc, yc = xdep - geom["dRA"], ydep - geom["dDec"]
 
+    # TODO: this block temp
+    # cos_PA = np.cos(geom["Omega"] * np.pi / 180)
+    # sin_PA = np.sin(geom["Omega"] * np.pi / 180)
+    # xdep =  xshift * cos_PA - yshift * sin_PA
+    # ydep = xshift * sin_PA + yshift * cos_PA
+    # xdep /= (geom["incl"] * np.pi / 180)
 
     # radial pixel coordinates
     rr = np.ravel(np.hypot(xc, yc))
@@ -72,7 +78,7 @@ def radialI(image, coords, geom, rescale_flux, bins=None):
     Is /= bin_counts
 
     # if the source is optically thick, rescale the deprojected I(r)
-    if rescale_flux:
+    if rescale_flux: # TODO
         Is *= np.cos(geom["incl"] * np.pi / 180)
 
     bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

From f7efc90f48060faed7528dd6f311cf4c8579da83 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 19 Oct 2023 14:00:58 -0400
Subject: [PATCH 050/109] start of onedim tests

---
 test/onedim_test.py | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

diff --git a/test/onedim_test.py b/test/onedim_test.py
index ed5d1388..dbddc495 100644
--- a/test/onedim_test.py
+++ b/test/onedim_test.py
@@ -1,4 +1,21 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
+from mpol.coordinates import GridCoords
 from mpol.onedim import radialI, radialV
+from mpol.precomposed import SimpleNet
+from mpol.utils import sky_gaussian_arcsec
+
+def test_radialV(coords, imager, dataset, generic_parameters):
+    # obtain a 1d radial visibility profile V(q) from 2d visibilities
+    coords = GridCoords(cell_size=0.005, npix=800)
+
+    g = sky_gaussian_arcsec(coords.sky_x_centers_2D, coords.sky_y_centers_2D)
+    
+    nchan = dataset.nchan
+    model = SimpleNet(coords=coords, nchan=nchan)
+
+def test_radialI(coords, imager, dataset, generic_parameters):
+    # obtain a 1d radial brightness profile I(r) from an image
+    nchan = dataset.nchan
+    model = SimpleNet(coords=coords, nchan=nchan)

From 8c7a9a07c11db8d44c7d8aa5dfb8070e0bf3d0aa Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Fri, 27 Oct 2023 12:56:48 -0400
Subject: [PATCH 051/109] radialI: remove 'rescale_flux'

---
 src/mpol/onedim.py | 13 +------------
 1 file changed, 1 insertion(+), 12 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 113000a6..ca4a4b13 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -3,7 +3,7 @@
 from mpol.geometry import observer_to_flat
 from mpol.utils import torch2npy
 
-def radialI(image, coords, geom, rescale_flux, bins=None):
+def radialI(image, coords, geom, bins=None):
     r"""
     Obtain a 1D (radial) brightness profile I(r) from an image.
 
@@ -25,13 +25,6 @@ def radialI(image, coords, geom, rescale_flux, bins=None):
                 Phase center offset in right ascension. Positive is west of north.
             "dDec" : float, unit=[arcsec]
                 Phase center offset in declination.
-    rescale_flux : bool
-        If True, the brightness values are rescaled to account for the 
-        difference between the inclined (observed) brightness and the 
-        assumed face-on brightness, assuming the emission is optically thick. 
-        The source's integrated (2D) flux is assumed to be:
-            :math:`F = \cos(i) \int_r^{r=R}{I(r) 2 \pi r dr}`.
-        No rescaling would be appropriate in the optically thin limit. 
     bins : array, default=None, unit=[arcsec]
         Radial bin edges to use in calculating I(r). If None, bins will span 
         the full image, with widths equal to the hypotenuse of the pixels
@@ -77,10 +70,6 @@ def radialI(image, coords, geom, rescale_flux, bins=None):
     Is, _ = np.histogram(a=rr, bins=bins, weights=np.ravel(image))
     Is /= bin_counts
 
-    # if the source is optically thick, rescale the deprojected I(r)
-    if rescale_flux: # TODO
-        Is *= np.cos(geom["incl"] * np.pi / 180)
-
     bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
 
     return bin_centers, Is

From 25b465881c8d3939c1c2f4933acb70215df2df8e Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Fri, 27 Oct 2023 13:01:09 -0400
Subject: [PATCH 052/109] radialI: var disambiguation

---
 src/mpol/onedim.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index ca4a4b13..d0d87ae9 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -48,7 +48,7 @@ def radialI(image, coords, geom, bins=None):
         )
 
     # shift image center to source center
-    xc, yc = xdep - geom["dRA"], ydep - geom["dDec"]
+    xc, yc = xx - geom["dRA"], yy - geom["dDec"]
 
     # TODO: this block temp
     # cos_PA = np.cos(geom["Omega"] * np.pi / 180)

From d890a706b98dca78cb4fb692a285bb6b76782c96 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Fri, 27 Oct 2023 13:01:35 -0400
Subject: [PATCH 053/109] radialI: custom deprojection

---
 src/mpol/onedim.py | 19 ++++++-------------
 1 file changed, 6 insertions(+), 13 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index d0d87ae9..348bcd7b 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -40,22 +40,15 @@ def radialI(image, coords, geom, bins=None):
     # projected Cartesian pixel coordinates [arcsec]
     xx, yy = coords.sky_x_centers_2D, coords.sky_y_centers_2D
 
-    # deproject and rotate image 
-    xdep, ydep = observer_to_flat(xx, yy,
-        omega=geom["omega"] * np.pi / 180, 
-        incl=geom["incl"] * np.pi / 180,
-        Omega=geom["Omega"] * np.pi / 180,
-        )
-
     # shift image center to source center
     xc, yc = xx - geom["dRA"], yy - geom["dDec"]
 
-    # TODO: this block temp
-    # cos_PA = np.cos(geom["Omega"] * np.pi / 180)
-    # sin_PA = np.sin(geom["Omega"] * np.pi / 180)
-    # xdep =  xshift * cos_PA - yshift * sin_PA
-    # ydep = xshift * sin_PA + yshift * cos_PA
-    # xdep /= (geom["incl"] * np.pi / 180)
+    # deproject image
+    cos_PA = np.cos(geom["Omega"] * np.pi / 180)
+    sin_PA = np.sin(geom["Omega"] * np.pi / 180)
+    xd = xc * cos_PA - yc * sin_PA
+    yd = xc * sin_PA + yc * cos_PA
+    xd /= np.cos(geom["incl"] * np.pi / 180)
 
     # radial pixel coordinates
     rr = np.ravel(np.hypot(xc, yc))

From 843b0e2dc73645376c9a9248f4d97c477ea93642 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Fri, 27 Oct 2023 13:03:05 -0400
Subject: [PATCH 054/109] radialI: more precise comments

---
 src/mpol/onedim.py | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 348bcd7b..cde32128 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -50,17 +50,20 @@ def radialI(image, coords, geom, bins=None):
     yd = xc * sin_PA + yc * cos_PA
     xd /= np.cos(geom["incl"] * np.pi / 180)
 
-    # radial pixel coordinates
-    rr = np.ravel(np.hypot(xc, yc))
+    # deprojected radial coordinates
+    rr = np.ravel(np.hypot(xd, yd))
 
     if bins is None:
         step = np.hypot(coords.cell_size, coords.cell_size)
         bins = np.arange(0.0, max(rr), step)
 
-    # number of points in radial bins
     bin_counts, bin_edges = np.histogram(a=rr, bins=bins, weights=None)
-    # brightness in radial bins
+
+    # cumulative binned brightness in each annulus
     Is, _ = np.histogram(a=rr, bins=bins, weights=np.ravel(image))
+
+
+    # average binned brightness in each annulus
     Is /= bin_counts
 
     bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

From 2679a2f4eb26914d58fb7ba1779907204ac689c8 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Fri, 27 Oct 2023 13:03:23 -0400
Subject: [PATCH 055/109] radialI: sensible default bins

---
 src/mpol/onedim.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index cde32128..54758715 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -54,8 +54,9 @@ def radialI(image, coords, geom, bins=None):
     rr = np.ravel(np.hypot(xd, yd))
 
     if bins is None:
+        # choose sensible bin size and range
         step = np.hypot(coords.cell_size, coords.cell_size)
-        bins = np.arange(0.0, max(rr), step)
+        bins = np.arange(min(abs(np.ravel(xc))), max(abs(np.ravel(xc))), step)
 
     bin_counts, bin_edges = np.histogram(a=rr, bins=bins, weights=None)
 

From 85e54a9e08b2188242bafd72c7c66a1a0172dc26 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Fri, 27 Oct 2023 13:03:31 -0400
Subject: [PATCH 056/109] radialI: mask empty bins

---
 src/mpol/onedim.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 54758715..260158d7 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -63,6 +63,9 @@ def radialI(image, coords, geom, bins=None):
     # cumulative binned brightness in each annulus
     Is, _ = np.histogram(a=rr, bins=bins, weights=np.ravel(image))
 
+    # mask empty bins
+    mask = (bin_counts == 0)
+    Is = np.ma.masked_where(mask, Is)
 
     # average binned brightness in each annulus
     Is /= bin_counts

From c009f23e86b8bc96c3ad18cdc0309e7b0b03b303 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Fri, 27 Oct 2023 13:04:51 -0400
Subject: [PATCH 057/109] radialV: replicate radialI layout

---
 src/mpol/onedim.py | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 260158d7..ba548d2a 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -138,16 +138,21 @@ def radialV(V, coords, geom, rescale_flux, bins=None):
     # convert back to [k\lambda]
     up /= 1e3
     vp /= 1e3
+
+    # deprojected baseline coordinates
     qq = np.hypot(up, vp) 
 
     if bins is None:
+        # choose sensible bin size and range
         step = np.hypot(coords.du, coords.dv)
         bins = np.arange(0.0, max(qq), step)
 
-    # number of points in radial bins
     bin_counts, bin_edges = np.histogram(a=qq, bins=bins, weights=None)
-    # V amplitudes in radial bins
+
+    # cumulative binned visibility amplitude in each annulus
     Vs, _ = np.histogram(a=qq, bins=bins, weights=Vp)
+
+    # average binned visibility amplitude in each annulus
     Vs /= bin_counts
     
     bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

From 10b7fb29233adc45398b9ea479d1bdda5022ca15 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Fri, 27 Oct 2023 13:04:59 -0400
Subject: [PATCH 058/109] radialV: mask empty bins

---
 src/mpol/onedim.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index ba548d2a..94c8b9f9 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -152,6 +152,10 @@ def radialV(V, coords, geom, rescale_flux, bins=None):
     # cumulative binned visibility amplitude in each annulus
     Vs, _ = np.histogram(a=qq, bins=bins, weights=Vp)
 
+    # mask empty bins
+    mask = (bin_counts == 0)
+    Vs = np.ma.masked_where(mask, Vs)
+
     # average binned visibility amplitude in each annulus
     Vs /= bin_counts
     

From 6934dfa88091336725ca0ecb1d7d17ea5e4137ef Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Fri, 27 Oct 2023 13:43:04 -0400
Subject: [PATCH 059/109] radialI: more general bin guess

---
 src/mpol/onedim.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 94c8b9f9..9c4dafab 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -56,7 +56,7 @@ def radialI(image, coords, geom, bins=None):
     if bins is None:
         # choose sensible bin size and range
         step = np.hypot(coords.cell_size, coords.cell_size)
-        bins = np.arange(min(abs(np.ravel(xc))), max(abs(np.ravel(xc))), step)
+        bins = np.arange(0.0, np.max((abs(xc.ravel()), abs(yc.ravel()))), step)
 
     bin_counts, bin_edges = np.histogram(a=rr, bins=bins, weights=None)
 

From c768689fea3620e486a8b87109a59983a3ef42a8 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Tue, 31 Oct 2023 15:52:43 -0400
Subject: [PATCH 060/109] radialI: set bin guess

---
 src/mpol/onedim.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 94c8b9f9..9c4dafab 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -56,7 +56,7 @@ def radialI(image, coords, geom, bins=None):
     if bins is None:
         # choose sensible bin size and range
         step = np.hypot(coords.cell_size, coords.cell_size)
-        bins = np.arange(min(abs(np.ravel(xc))), max(abs(np.ravel(xc))), step)
+        bins = np.arange(0.0, np.max((abs(xc.ravel()), abs(yc.ravel()))), step)
 
     bin_counts, bin_edges = np.histogram(a=rr, bins=bins, weights=None)
 

From 5e4b1929556f76157db0ec44f3bf292a328badbd Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Tue, 31 Oct 2023 18:02:56 -0400
Subject: [PATCH 061/109] add pytest fixtures for mock 1d dataset

---
 test/conftest.py | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/test/conftest.py b/test/conftest.py
index 2b6fc451..63d0f87f 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -117,6 +117,28 @@ def dataset_cont(mock_visibility_data_cont, coords):
     return averager.to_pytorch_dataset()
 
 
+@pytest.fixture(scope="session")
+def mock_1d_archive():
+    # use astropy routines to cache data
+    fname = download_file(
+        "https://zenodo.org/record/{:d}}/files/mock_profile_image.npz".format(zenodo_version),
+        cache=True,
+        pkgname="mpol",
+    )
+    
+    return np.load(fname)
+
+@pytest.fixture
+def mock_1d_models(mock_1d_archive):
+    m = mock_1d_archive
+    r = m['r']
+    i = m['i']
+    i2d = m['i2d']
+    xmax = ymax = m['xmax']
+    
+    return r, i, i2d, xmax, ymax
+
+
 @pytest.fixture
 def crossvalidation_products(mock_visibility_data):
     # test the crossvalidation with a smaller set of image / Fourier coordinates than normal,

From 9655277f34db5c7bbaf640ae4d56b6a46f8fe244 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Tue, 31 Oct 2023 19:30:26 -0400
Subject: [PATCH 062/109] add utility function to center a numpy image

---
 src/mpol/utils.py | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)

diff --git a/src/mpol/utils.py b/src/mpol/utils.py
index ac8ced23..c8806bdb 100644
--- a/src/mpol/utils.py
+++ b/src/mpol/utils.py
@@ -10,6 +10,30 @@ def torch2npy(tensor):
     return tensor.detach().cpu().numpy()
 
 
+def center_np_image(image):
+    """Wrap a numpy image array so that the 0 index is in the image center.
+    Adapted from frank.utilities.make_image (https://github.com/discsim/frank)
+
+    Parameters
+    ----------
+    image : 2D array
+        The numpy image with index 0 in the corner
+    
+    Returns
+    -------
+    image_wrapped : 2D array
+        The numpy image with index 0 in the center
+    """
+    tmp, image_wrapped = image.copy(), image.copy()
+
+    Nx, Ny = tmp.shape[0], tmp.shape[1]
+
+    tmp[:Nx//2,], tmp[Nx//2:] = image[Nx//2:], image[:Nx//2]
+    image_wrapped[:, :Ny//2], image_wrapped[:, Ny//2:] = tmp[:, Ny//2:], tmp[:, :Ny//2]
+
+    return image_wrapped
+
+
 def ground_cube_to_packed_cube(ground_cube):
     r"""
     Converts a Ground Cube to a Packed Visibility Cube for visibility-plane work. See Units and Conventions for more details.

From 062e2f19f4dcef8853e1317703e056ce20a31186 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Tue, 31 Oct 2023 19:30:50 -0400
Subject: [PATCH 063/109] add utility function to convert np to ImageCube

---
 src/mpol/utils.py | 38 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 38 insertions(+)

diff --git a/src/mpol/utils.py b/src/mpol/utils.py
index c8806bdb..8ecc2df1 100644
--- a/src/mpol/utils.py
+++ b/src/mpol/utils.py
@@ -5,6 +5,8 @@
 
 from .constants import arcsec, c_ms, cc, deg, kB
 
+from mpol.images import ImageCube 
+
 def torch2npy(tensor):
     """Make a copy of a PyTorch tensor on the CPU in numpy format, e.g. for plotting"""
     return tensor.detach().cpu().numpy()
@@ -34,6 +36,42 @@ def center_np_image(image):
     return image_wrapped
 
 
+def np_to_imagecube(image, coords, nchan=1, wrap=False):
+    """Convenience function for converting a numpy image into an MPoL ImageCube 
+    tensor (see mpol.images.ImageCube)
+
+    Parameters
+    ----------
+    image : array
+        An image in numpy format
+    coords : `mpol.coordinates.GridCoords` object
+        Instance of the `mpol.coordinates.GridCoords` class
+    nchan : int, default=1
+        Number of channels in the image. Default assumes a single 2D image
+    wrap : bool, default=False 
+        Whether to wrap the numpy image so that index 0 is in the image center
+        (FFT algorithms typically place index 0 in the image corner)
+
+    Returns
+    -------
+    icube : `mpol.images.ImageCube` object
+        The image cube tensor
+    """
+    if wrap: 
+        # move the 0 index to the image center 
+        image = center_np_image(image)
+
+    # broadcast image to (nchan, npix, npix)
+    img_packed_cube = np.broadcast_to(image, 
+                                      (nchan, coords.npix, coords.npix)).copy()
+
+    # convert to pytorch tensor
+    img_packed_tensor = torch.from_numpy(img_packed_cube)
+
+    # insert into ImageCube layer
+    return ImageCube(coords=coords, nchan=nchan, cube=img_packed_tensor)
+
+
 def ground_cube_to_packed_cube(ground_cube):
     r"""
     Converts a Ground Cube to a Packed Visibility Cube for visibility-plane work. See Units and Conventions for more details.

From 1b76a3242a8306e9ed65eb821dabd3801c62da93 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Tue, 31 Oct 2023 19:44:29 -0400
Subject: [PATCH 064/109] move np to ImageCube out of utils

---
 src/mpol/utils.py | 38 --------------------------------------
 1 file changed, 38 deletions(-)

diff --git a/src/mpol/utils.py b/src/mpol/utils.py
index 8ecc2df1..c8806bdb 100644
--- a/src/mpol/utils.py
+++ b/src/mpol/utils.py
@@ -5,8 +5,6 @@
 
 from .constants import arcsec, c_ms, cc, deg, kB
 
-from mpol.images import ImageCube 
-
 def torch2npy(tensor):
     """Make a copy of a PyTorch tensor on the CPU in numpy format, e.g. for plotting"""
     return tensor.detach().cpu().numpy()
@@ -36,42 +34,6 @@ def center_np_image(image):
     return image_wrapped
 
 
-def np_to_imagecube(image, coords, nchan=1, wrap=False):
-    """Convenience function for converting a numpy image into an MPoL ImageCube 
-    tensor (see mpol.images.ImageCube)
-
-    Parameters
-    ----------
-    image : array
-        An image in numpy format
-    coords : `mpol.coordinates.GridCoords` object
-        Instance of the `mpol.coordinates.GridCoords` class
-    nchan : int, default=1
-        Number of channels in the image. Default assumes a single 2D image
-    wrap : bool, default=False 
-        Whether to wrap the numpy image so that index 0 is in the image center
-        (FFT algorithms typically place index 0 in the image corner)
-
-    Returns
-    -------
-    icube : `mpol.images.ImageCube` object
-        The image cube tensor
-    """
-    if wrap: 
-        # move the 0 index to the image center 
-        image = center_np_image(image)
-
-    # broadcast image to (nchan, npix, npix)
-    img_packed_cube = np.broadcast_to(image, 
-                                      (nchan, coords.npix, coords.npix)).copy()
-
-    # convert to pytorch tensor
-    img_packed_tensor = torch.from_numpy(img_packed_cube)
-
-    # insert into ImageCube layer
-    return ImageCube(coords=coords, nchan=nchan, cube=img_packed_tensor)
-
-
 def ground_cube_to_packed_cube(ground_cube):
     r"""
     Converts a Ground Cube to a Packed Visibility Cube for visibility-plane work. See Units and Conventions for more details.

From efa68077b385f23c80e2785f3ed1eb6570e58842 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Tue, 31 Oct 2023 19:46:28 -0400
Subject: [PATCH 065/109] Move np to ImageCube into images.py

---
 src/mpol/images.py | 36 ++++++++++++++++++++++++++++++++++++
 1 file changed, 36 insertions(+)

diff --git a/src/mpol/images.py b/src/mpol/images.py
index 0185766e..9dae120c 100644
--- a/src/mpol/images.py
+++ b/src/mpol/images.py
@@ -316,3 +316,39 @@ def to_FITS(self, fname="cube.fits", overwrite=False, header_kwargs=None):
         hdul.writeto(fname, overwrite=overwrite)
 
         hdul.close()
+
+
+def np_to_imagecube(image, coords, nchan=1, wrap=False):
+    """Convenience function for converting a numpy image into an MPoL ImageCube 
+    tensor (see mpol.images.ImageCube)
+
+    Parameters
+    ----------
+    image : array
+        An image in numpy format
+    coords : `mpol.coordinates.GridCoords` object
+        Instance of the `mpol.coordinates.GridCoords` class
+    nchan : int, default=1
+        Number of channels in the image. Default assumes a single 2D image
+    wrap : bool, default=False 
+        Whether to wrap the numpy image so that index 0 is in the image center
+        (FFT algorithms typically place index 0 in the image corner)
+
+    Returns
+    -------
+    icube : `mpol.images.ImageCube` object
+        The image cube tensor
+    """
+    if wrap: 
+        # move the 0 index to the image center 
+        image = center_np_image(image)
+
+    # broadcast image to (nchan, npix, npix)
+    img_packed_cube = np.broadcast_to(image, 
+                                      (nchan, coords.npix, coords.npix)).copy()
+
+    # convert to pytorch tensor
+    img_packed_tensor = torch.from_numpy(img_packed_cube)
+
+    # insert into ImageCube layer
+    return ImageCube(coords=coords, nchan=nchan, cube=img_packed_tensor)

From ec422f659a68ea4cf18c567487f1be2929782274 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Tue, 31 Oct 2023 19:47:45 -0400
Subject: [PATCH 066/109] np__to_imagecube: add import

---
 src/mpol/images.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpol/images.py b/src/mpol/images.py
index 9dae120c..b0f9b2dc 100644
--- a/src/mpol/images.py
+++ b/src/mpol/images.py
@@ -341,7 +341,7 @@ def np_to_imagecube(image, coords, nchan=1, wrap=False):
     """
     if wrap: 
         # move the 0 index to the image center 
-        image = center_np_image(image)
+        image = utils.center_np_image(image)
 
     # broadcast image to (nchan, npix, npix)
     img_packed_cube = np.broadcast_to(image, 

From d99a025768edf15aba6921919ec84817506d0404 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Tue, 31 Oct 2023 23:13:51 -0400
Subject: [PATCH 067/109] add processing steps to mock_1d_image_models

---
 test/conftest.py | 19 ++++++++++++++++++-
 1 file changed, 18 insertions(+), 1 deletion(-)

diff --git a/test/conftest.py b/test/conftest.py
index 63d0f87f..86333736 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -126,7 +126,24 @@ def mock_1d_archive():
         pkgname="mpol",
     )
     
-    return np.load(fname)
+    return np.load(fname, allow_pickle=True)
+
+@pytest.fixture
+def mock_1d_image_models(mock_1d_archive):
+    m = mock_1d_archive
+    r = m['r']
+    i = m['i']
+    i2d = m['i2d']
+    xmax = ymax = m['xmax']
+    geom = m['geometry']
+
+    coords = coordinates.GridCoords(cell_size=xmax * 2 / i2d.shape[0], npix=i2d.shape[0])
+
+    i2d_packed_cube = np.broadcast_to(i2d, (1, coords.npix, coords.npix)).copy()
+    i2d_packed_tensor = torch.from_numpy(i2d_packed_cube)
+    i2d_icube = images.ImageCube(coords=coords, nchan=1, cube=i2d_packed_tensor)
+
+    return r, i, i2d, xmax, ymax, geom, coords, i2d_icube
 
 @pytest.fixture
 def mock_1d_models(mock_1d_archive):

From 343e3773a1c81dd27d852bb0d0133bdf7ac7b48f Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 00:31:17 -0400
Subject: [PATCH 068/109] mock 1d data: new filename

---
 test/conftest.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/conftest.py b/test/conftest.py
index 86333736..d1dddcd9 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -121,7 +121,7 @@ def dataset_cont(mock_visibility_data_cont, coords):
 def mock_1d_archive():
     # use astropy routines to cache data
     fname = download_file(
-        "https://zenodo.org/record/{:d}}/files/mock_profile_image.npz".format(zenodo_version),
+        "https://zenodo.org/record/{:d}/files/mock_disk_1d.npz".format(zenodo_version),
         cache=True,
         pkgname="mpol",
     )

From 9ef73ef79e25d9e6b87c2a9dcb195773513aacbf Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 00:31:44 -0400
Subject: [PATCH 069/109] simplify mock_1d_image_model fixture

---
 test/conftest.py | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/test/conftest.py b/test/conftest.py
index d1dddcd9..d7bb4b0c 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -128,8 +128,9 @@ def mock_1d_archive():
     
     return np.load(fname, allow_pickle=True)
 
+
 @pytest.fixture
-def mock_1d_image_models(mock_1d_archive):
+def mock_1d_image_model(mock_1d_archive):
     m = mock_1d_archive
     r = m['r']
     i = m['i']
@@ -139,11 +140,8 @@ def mock_1d_image_models(mock_1d_archive):
 
     coords = coordinates.GridCoords(cell_size=xmax * 2 / i2d.shape[0], npix=i2d.shape[0])
 
-    i2d_packed_cube = np.broadcast_to(i2d, (1, coords.npix, coords.npix)).copy()
-    i2d_packed_tensor = torch.from_numpy(i2d_packed_cube)
-    i2d_icube = images.ImageCube(coords=coords, nchan=1, cube=i2d_packed_tensor)
+    return r, i, i2d, xmax, ymax, geom, coords
 
-    return r, i, i2d, xmax, ymax, geom, coords, i2d_icube
 
 @pytest.fixture
 def mock_1d_models(mock_1d_archive):

From bc21790808169e9960be7d8c6204d9734115d709 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 00:33:13 -0400
Subject: [PATCH 070/109] simplify onedim_test design

---
 test/onedim_test.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/test/onedim_test.py b/test/onedim_test.py
index dbddc495..0f25fcbc 100644
--- a/test/onedim_test.py
+++ b/test/onedim_test.py
@@ -1,10 +1,9 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from mpol.coordinates import GridCoords
 from mpol.onedim import radialI, radialV
-from mpol.precomposed import SimpleNet
-from mpol.utils import sky_gaussian_arcsec
+from mpol.plot import plot_image
+
 
 def test_radialV(coords, imager, dataset, generic_parameters):
     # obtain a 1d radial visibility profile V(q) from 2d visibilities

From 43f3311cba9a3a439059eb968fc5d520cb3674b0 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 00:33:56 -0400
Subject: [PATCH 071/109] implement test of radialI

---
 test/onedim_test.py | 32 ++++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/test/onedim_test.py b/test/onedim_test.py
index 0f25fcbc..96e8a438 100644
--- a/test/onedim_test.py
+++ b/test/onedim_test.py
@@ -4,6 +4,38 @@
 from mpol.onedim import radialI, radialV
 from mpol.plot import plot_image
 
+def test_radialI(mock_1d_image_model):
+    # obtain a 1d radial brightness profile I(r) from an image    
+    
+    r, i, i2d, _, _, geom, coords = mock_1d_image_model
+
+    rtest, itest = radialI(i2d, coords, geom, bins=None)
+
+    expected = [
+        5.79780326e+10, 2.47990375e+10, 4.19794053e+09, 1.63165616e+10, 
+        2.56197452e+10, 1.86014523e+10, 1.39800643e+10, 1.14935415e+10, 
+        1.58898181e+10, 1.36344624e+10, 1.14059388e+10, 1.18389766e+10,
+        1.18678220e+10, 1.02746571e+10, 8.18228608e+09, 5.40044021e+09, 
+        2.63008657e+09, 5.61017562e+08, 4.11251459e+08, 4.85212055e+08, 
+        7.86240201e+08, 3.80008818e+09, 6.99254078e+09, 4.63422518e+09,
+        1.47955225e+09, 3.54437101e+08, 5.40245124e+08, 2.40733707e+08, 
+        3.75558288e+08, 3.09550836e+08, 3.89050755e+08, 2.94034065e+08, 
+        2.28420989e+08, 2.02024946e+08, 8.01676870e+08, 3.32150158e+09,
+        5.39650629e+09, 4.14278723e+09, 2.15660813e+09, 1.26767762e+09, 
+        9.44844287e+08, 8.38469945e+08, 6.67265501e+08, 7.01335008e+08, 
+        5.23378469e+08, 3.81449784e+08, 2.90773229e+08, 2.11627913e+08,
+        1.49714706e+08, 5.25051148e+07, 1.32583044e+08, 4.70658247e+07, 
+        3.52859146e+07, 8.45017542e+06, 1.55314270e+07, 9.15833896e+06, 
+        5.27578119e+06, 2.08955802e+06, 2.20079612e+06, 1.36049060e+07,
+        2.26196295e+07, 5.66055107e+06, 1.30125309e+05, 1.23963592e+04, 
+        3.63853801e+03, 2.74032116e+03, 1.21440814e+03, 1.00115299e+03, 
+        9.39440210e+02, 8.74909523e+02, 1.01296747e+03, 3.19296404e-02,
+        0.00000000e+00, 0.00000000e+00
+        ]
+
+    np.testing.assert_allclose(itest, expected, rtol=1e-6,
+                               err_msg="test_radialI")
+
 
 def test_radialV(coords, imager, dataset, generic_parameters):
     # obtain a 1d radial visibility profile V(q) from 2d visibilities

From eecfef4d66142946a84dc664c3a4bb20d8893011 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 00:34:15 -0400
Subject: [PATCH 072/109] add diag plot to test_radialI

---
 test/onedim_test.py | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/test/onedim_test.py b/test/onedim_test.py
index 96e8a438..94c83d75 100644
--- a/test/onedim_test.py
+++ b/test/onedim_test.py
@@ -11,6 +11,15 @@ def test_radialI(mock_1d_image_model):
 
     rtest, itest = radialI(i2d, coords, geom, bins=None)
 
+    _, ax = plt.subplots(ncols=2)
+
+    plot_image(i2d, extent=coords.img_ext, ax=ax[0], clab='Jy sr$^{-2}$')
+
+    ax[0].title('AS 209-like profile.\nGeometry: {:}'.format(geom))
+
+    ax[1].plot(r, i, 'k', label='truth')
+    ax[1].plot(rtest, itest, 'r.-', label='result')
+
     expected = [
         5.79780326e+10, 2.47990375e+10, 4.19794053e+09, 1.63165616e+10, 
         2.56197452e+10, 1.86014523e+10, 1.39800643e+10, 1.14935415e+10, 

From c61f8641ef277214244623d715ad035773a0635e Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 00:39:24 -0400
Subject: [PATCH 073/109] test_radialI: save diag fig

---
 test/onedim_test.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/test/onedim_test.py b/test/onedim_test.py
index 94c83d75..27389ec9 100644
--- a/test/onedim_test.py
+++ b/test/onedim_test.py
@@ -4,14 +4,14 @@
 from mpol.onedim import radialI, radialV
 from mpol.plot import plot_image
 
-def test_radialI(mock_1d_image_model):
+def test_radialI(mock_1d_image_model, tmp_path):
     # obtain a 1d radial brightness profile I(r) from an image    
     
     r, i, i2d, _, _, geom, coords = mock_1d_image_model
 
     rtest, itest = radialI(i2d, coords, geom, bins=None)
 
-    _, ax = plt.subplots(ncols=2)
+    fig, ax = plt.subplots(ncols=2)
 
     plot_image(i2d, extent=coords.img_ext, ax=ax[0], clab='Jy sr$^{-2}$')
 
@@ -20,6 +20,9 @@ def test_radialI(mock_1d_image_model):
     ax[1].plot(r, i, 'k', label='truth')
     ax[1].plot(rtest, itest, 'r.-', label='result')
 
+    fig.savefig(tmp_path / "test_radialI.png", dpi=300)
+    plt.close("all")
+
     expected = [
         5.79780326e+10, 2.47990375e+10, 4.19794053e+09, 1.63165616e+10, 
         2.56197452e+10, 1.86014523e+10, 1.39800643e+10, 1.14935415e+10, 

From 8f1a610ca37b2809534ec5116e704a296c7a1dba Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 00:42:43 -0400
Subject: [PATCH 074/109] test_radialI: add plot ax labels

---
 test/onedim_test.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/test/onedim_test.py b/test/onedim_test.py
index 27389ec9..253b633b 100644
--- a/test/onedim_test.py
+++ b/test/onedim_test.py
@@ -13,12 +13,15 @@ def test_radialI(mock_1d_image_model, tmp_path):
 
     fig, ax = plt.subplots(ncols=2)
 
-    plot_image(i2d, extent=coords.img_ext, ax=ax[0], clab='Jy sr$^{-2}$')
+    plot_image(i2d, extent=coords.img_ext, ax=ax[0], clab='Jy / sr')
 
     ax[0].title('AS 209-like profile.\nGeometry: {:}'.format(geom))
 
     ax[1].plot(r, i, 'k', label='truth')
     ax[1].plot(rtest, itest, 'r.-', label='result')
+    
+    ax[1].set_xlabel('r [arcsec]')
+    ax[1].set_ylabel('I [Jy / sr]')
 
     fig.savefig(tmp_path / "test_radialI.png", dpi=300)
     plt.close("all")

From 3cc2316d2e882c935cea860bc877aaef0274beeb Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 01:04:34 -0400
Subject: [PATCH 075/109] test_radialI: set bins

---
 test/onedim_test.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/test/onedim_test.py b/test/onedim_test.py
index 253b633b..9edc0eb4 100644
--- a/test/onedim_test.py
+++ b/test/onedim_test.py
@@ -9,7 +9,8 @@ def test_radialI(mock_1d_image_model, tmp_path):
     
     r, i, i2d, _, _, geom, coords = mock_1d_image_model
 
-    rtest, itest = radialI(i2d, coords, geom, bins=None)
+    bins = np.linspace(0, 2.0, 75)
+    rtest, itest = radialI(i2d, coords, geom, bins=bins)
 
     fig, ax = plt.subplots(ncols=2)
 

From 7380ecab5dbf144ff19e3fc6641804cc4ff3ca9a Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 10:41:32 -0400
Subject: [PATCH 076/109] radialV: less conservative step size

---
 src/mpol/onedim.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 9c4dafab..f4dee157 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -144,7 +144,7 @@ def radialV(V, coords, geom, rescale_flux, bins=None):
 
     if bins is None:
         # choose sensible bin size and range
-        step = np.hypot(coords.du, coords.dv)
+        step = np.hypot(coords.du, coords.dv) / 2
         bins = np.arange(0.0, max(qq), step)
 
     bin_counts, bin_edges = np.histogram(a=qq, bins=bins, weights=None)

From df6b2a6d2e2349db61fd2fea9d46685a4804b744 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 14:42:13 -0400
Subject: [PATCH 077/109] mock_1d_image_model test fixture: clarify naming

---
 src/mpol/onedim.py |  3 ---
 test/conftest.py   | 11 ++++++-----
 2 files changed, 6 insertions(+), 8 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index f4dee157..e8f4a91c 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -1,8 +1,5 @@
 import numpy as np 
 
-from mpol.geometry import observer_to_flat
-from mpol.utils import torch2npy
-
 def radialI(image, coords, geom, bins=None):
     r"""
     Obtain a 1D (radial) brightness profile I(r) from an image.
diff --git a/test/conftest.py b/test/conftest.py
index d7bb4b0c..937f4532 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -132,15 +132,16 @@ def mock_1d_archive():
 @pytest.fixture
 def mock_1d_image_model(mock_1d_archive):
     m = mock_1d_archive
-    r = m['r']
-    i = m['i']
-    i2d = m['i2d']
+    rtrue = m['rtrue']
+    itrue = m['itrue']
+    i2dtrue = m['i2dtrue']
     xmax = ymax = m['xmax']
     geom = m['geometry']
 
-    coords = coordinates.GridCoords(cell_size=xmax * 2 / i2d.shape[0], npix=i2d.shape[0])
+    coords = coordinates.GridCoords(cell_size=xmax * 2 / i2dtrue.shape[0], 
+                                    npix=i2dtrue.shape[0])
 
-    return r, i, i2d, xmax, ymax, geom, coords
+    return rtrue, itrue, i2dtrue, xmax, ymax, geom, coords
 
 
 @pytest.fixture

From 704d77094f430c9ec66d6ed496c88b5c23ed26fb Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 14:42:37 -0400
Subject: [PATCH 078/109] add mock_1d_vis_model test fixture

---
 test/conftest.py | 21 ++++++++++++++-------
 1 file changed, 14 insertions(+), 7 deletions(-)

diff --git a/test/conftest.py b/test/conftest.py
index 937f4532..d266b6f1 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -145,14 +145,21 @@ def mock_1d_image_model(mock_1d_archive):
 
 
 @pytest.fixture
-def mock_1d_models(mock_1d_archive):
+def mock_1d_vis_model(mock_1d_archive):
     m = mock_1d_archive
-    r = m['r']
-    i = m['i']
-    i2d = m['i2d']
-    xmax = ymax = m['xmax']
-    
-    return r, i, i2d, xmax, ymax
+    Vtrue = m['vis']
+    Vtrue_dep = m['vis_dep']
+    q_dep = m['baselines_dep']
+    geom = m['geometry']
+
+    xmax = m['xmax']
+    i2dtrue = m['i2dtrue']
+
+    coords = coordinates.GridCoords(cell_size=xmax * 2 / i2dtrue.shape[0],
+                                    npix=i2dtrue.shape[0])
+
+    return Vtrue, Vtrue_dep, q_dep, geom, coords
+
 
 
 @pytest.fixture

From 346a59887c1972d2ebd15efacc53d526a2839b2f Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 14:46:12 -0400
Subject: [PATCH 079/109] remove unneeded func

---
 src/mpol/utils.py | 24 ------------------------
 1 file changed, 24 deletions(-)

diff --git a/src/mpol/utils.py b/src/mpol/utils.py
index 65178bb9..92096b19 100644
--- a/src/mpol/utils.py
+++ b/src/mpol/utils.py
@@ -10,30 +10,6 @@ def torch2npy(tensor):
     return tensor.detach().cpu().numpy()
 
 
-def center_np_image(image):
-    """Wrap a numpy image array so that the 0 index is in the image center.
-    Adapted from frank.utilities.make_image (https://github.com/discsim/frank)
-
-    Parameters
-    ----------
-    image : 2D array
-        The numpy image with index 0 in the corner
-    
-    Returns
-    -------
-    image_wrapped : 2D array
-        The numpy image with index 0 in the center
-    """
-    tmp, image_wrapped = image.copy(), image.copy()
-
-    Nx, Ny = tmp.shape[0], tmp.shape[1]
-
-    tmp[:Nx//2,], tmp[Nx//2:] = image[Nx//2:], image[:Nx//2]
-    image_wrapped[:, :Ny//2], image_wrapped[:, Ny//2:] = tmp[:, Ny//2:], tmp[:, :Ny//2]
-
-    return image_wrapped
-
-
 def ground_cube_to_packed_cube(ground_cube):
     r"""
     Converts a Ground Cube to a Packed Visibility Cube for visibility-plane work. See Units and Conventions for more details.

From 647cee7f1671aff8a7a9b0087ae70530793c7419 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 14:46:58 -0400
Subject: [PATCH 080/109] test_radialI: update test points

---
 test/onedim_test.py | 48 +++++++++++++++++++++++++--------------------
 1 file changed, 27 insertions(+), 21 deletions(-)

diff --git a/test/onedim_test.py b/test/onedim_test.py
index 9edc0eb4..d4ed309c 100644
--- a/test/onedim_test.py
+++ b/test/onedim_test.py
@@ -9,7 +9,7 @@ def test_radialI(mock_1d_image_model, tmp_path):
     
     r, i, i2d, _, _, geom, coords = mock_1d_image_model
 
-    bins = np.linspace(0, 2.0, 75)
+    bins = np.linspace(0, 2.0, 100)
     rtest, itest = radialI(i2d, coords, geom, bins=bins)
 
     fig, ax = plt.subplots(ncols=2)
@@ -19,7 +19,7 @@ def test_radialI(mock_1d_image_model, tmp_path):
     ax[0].title('AS 209-like profile.\nGeometry: {:}'.format(geom))
 
     ax[1].plot(r, i, 'k', label='truth')
-    ax[1].plot(rtest, itest, 'r.-', label='result')
+    ax[1].plot(rtest, itest, 'r.-', label='recovery')
     
     ax[1].set_xlabel('r [arcsec]')
     ax[1].set_ylabel('I [Jy / sr]')
@@ -28,25 +28,31 @@ def test_radialI(mock_1d_image_model, tmp_path):
     plt.close("all")
 
     expected = [
-        5.79780326e+10, 2.47990375e+10, 4.19794053e+09, 1.63165616e+10, 
-        2.56197452e+10, 1.86014523e+10, 1.39800643e+10, 1.14935415e+10, 
-        1.58898181e+10, 1.36344624e+10, 1.14059388e+10, 1.18389766e+10,
-        1.18678220e+10, 1.02746571e+10, 8.18228608e+09, 5.40044021e+09, 
-        2.63008657e+09, 5.61017562e+08, 4.11251459e+08, 4.85212055e+08, 
-        7.86240201e+08, 3.80008818e+09, 6.99254078e+09, 4.63422518e+09,
-        1.47955225e+09, 3.54437101e+08, 5.40245124e+08, 2.40733707e+08, 
-        3.75558288e+08, 3.09550836e+08, 3.89050755e+08, 2.94034065e+08, 
-        2.28420989e+08, 2.02024946e+08, 8.01676870e+08, 3.32150158e+09,
-        5.39650629e+09, 4.14278723e+09, 2.15660813e+09, 1.26767762e+09, 
-        9.44844287e+08, 8.38469945e+08, 6.67265501e+08, 7.01335008e+08, 
-        5.23378469e+08, 3.81449784e+08, 2.90773229e+08, 2.11627913e+08,
-        1.49714706e+08, 5.25051148e+07, 1.32583044e+08, 4.70658247e+07, 
-        3.52859146e+07, 8.45017542e+06, 1.55314270e+07, 9.15833896e+06, 
-        5.27578119e+06, 2.08955802e+06, 2.20079612e+06, 1.36049060e+07,
-        2.26196295e+07, 5.66055107e+06, 1.30125309e+05, 1.23963592e+04, 
-        3.63853801e+03, 2.74032116e+03, 1.21440814e+03, 1.00115299e+03, 
-        9.39440210e+02, 8.74909523e+02, 1.01296747e+03, 3.19296404e-02,
-        0.00000000e+00, 0.00000000e+00
+        6.40747314e+10, 4.01920507e+10, 1.44803534e+10, 2.94238627e+09, 
+        1.28782935e+10, 2.68613199e+10, 2.26564596e+10, 1.81151845e+10, 
+        1.52128965e+10, 1.05640352e+10, 1.33411204e+10, 1.61124502e+10,
+        1.41500539e+10, 1.20121195e+10, 1.11770326e+10, 1.19676913e+10, 
+        1.20941686e+10, 1.09498286e+10, 9.74236410e+09, 7.99589196e+09, 
+        5.94787809e+09, 3.82074946e+09, 1.80823933e+09, 4.48414819e+08,
+        3.17808840e+08, 5.77317876e+08, 3.98851281e+08, 8.06459834e+08, 
+        2.88706161e+09, 6.09577814e+09, 6.98556762e+09, 4.47436415e+09, 
+        1.89511273e+09, 5.96604356e+08, 3.44571640e+08, 5.65906765e+08,
+        2.85854589e+08, 2.67589013e+08, 3.98357054e+08, 2.97052261e+08, 
+        3.82744591e+08, 3.52239791e+08, 2.74336969e+08, 2.28425747e+08, 
+        1.82290043e+08, 3.16077299e+08, 1.18465538e+09, 3.32239287e+09,
+        5.26718846e+09, 5.16458748e+09, 3.58114198e+09, 2.13431954e+09, 
+        1.40936556e+09, 1.04032244e+09, 9.24050422e+08, 8.46829316e+08, 
+        6.80909295e+08, 6.83812465e+08, 6.91856237e+08, 5.29227136e+08,
+        3.97557293e+08, 3.54893419e+08, 2.60997039e+08, 2.09306498e+08, 
+        1.93930693e+08, 6.97032407e+07, 6.66090083e+07, 1.40079594e+08, 
+        7.21775931e+07, 3.23902663e+07, 3.35932300e+07, 7.63318789e+06,
+        1.29740981e+07, 1.44300351e+07, 8.06249624e+06, 5.85567843e+06, 
+        1.42637174e+06, 3.21445075e+06, 1.83763663e+06, 1.16926652e+07, 
+        2.46918188e+07, 1.60206523e+07, 3.26596592e+06, 1.27837319e+05,
+        2.27104612e+04, 4.77267063e+03, 2.90467640e+03, 2.88482230e+03, 
+        1.43402521e+03, 1.54791996e+03, 7.23397046e+02, 1.02561351e+03, 
+        5.24845888e+02, 1.47320552e+03, 7.40419174e+02, 4.59029378e-03,
+        0.00000000e+00, 0.00000000e+00, 0.00000000e+00
         ]
 
     np.testing.assert_allclose(itest, expected, rtol=1e-6,

From d83bf2a3fcb6690b8ca0513bfd9890d39cb4bb32 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 14:47:40 -0400
Subject: [PATCH 081/109] test_radialV: add test assert

---
 test/onedim_test.py | 45 +++++++++++++++++++++++++++++++++++----------
 1 file changed, 35 insertions(+), 10 deletions(-)

diff --git a/test/onedim_test.py b/test/onedim_test.py
index d4ed309c..51ef90f0 100644
--- a/test/onedim_test.py
+++ b/test/onedim_test.py
@@ -59,16 +59,41 @@ def test_radialI(mock_1d_image_model, tmp_path):
                                err_msg="test_radialI")
 
 
-def test_radialV(coords, imager, dataset, generic_parameters):
+def test_radialV(mock_1d_vis_model, tmp_path):
     # obtain a 1d radial visibility profile V(q) from 2d visibilities
-    coords = GridCoords(cell_size=0.005, npix=800)
 
-    g = sky_gaussian_arcsec(coords.sky_x_centers_2D, coords.sky_y_centers_2D)
-    
-    nchan = dataset.nchan
-    model = SimpleNet(coords=coords, nchan=nchan)
+    Vtrue, Vtrue_dep, q_dep, geom, coords = mock_1d_vis_model
+
+    qtest, Vtest = radialV(Vtrue, coords, geom, rescale_flux=True, bins=None)
+
+
+    expected = [
+        2.53998336e-01,  1.59897580e-01,  8.59460326e-02,  7.42189236e-02,  
+        5.75440687e-02,  1.81324892e-02, -2.92922689e-03,  3.14354163e-03,  
+        6.72339399e-03, -8.54632390e-03, -1.73385166e-02, -4.03826092e-03,
+        1.45595908e-02,  1.61681713e-02,  5.93475866e-03,  2.45555912e-04, 
+        -7.05014619e-04, -6.09266430e-03, -1.02454088e-02, -2.80944776e-03,  
+        8.58212558e-03, 8.39132158e-03, -6.52523293e-04,  -4.34778158e-03,
+        1.08035980e-04,  3.40903070e-03,  2.26682041e-03,  2.42465437e-03,  
+        5.07968926e-03,  4.83377443e-03, 1.26300648e-03,  1.11930639e-03,  
+        6.45157513e-03, 1.05751150e-02,  9.14016956e-03,  5.92209210e-03,
+        5.18455986e-03,  5.88802559e-03,  5.49315770e-03, 4.96398638e-03,  
+        5.81115311e-03,  5.95304063e-03, 3.16208083e-03, -1.71765038e-04, 
+        -4.64532628e-04, 1.12448670e-03,  1.84297313e-03,  1.48094594e-03,
+        1.12953770e-03,  1.01370816e-03,  6.57047907e-04, 1.37570722e-04,  
+        3.00648884e-04,  1.04847404e-03, 1.16856102e-03,  3.08940761e-04, 
+        -5.65721897e-04, -8.38907531e-04, -8.71976125e-04, -1.09567680e-03,
+        -1.42077854e-03, -1.33702627e-03, -9.96839047e-04, -1.16400192e-03, 
+        -1.43584618e-03, -1.07454472e-03, -6.44900590e-04, -4.86165342e-04, 
+        -1.96851463e-04, 5.04190986e-05,  5.73950179e-05,  2.79905736e-04,
+        7.52685847e-04,  1.12546048e-03,  1.37149548e-03, 1.35835560e-03,  
+        1.06470794e-03,  8.81423014e-04, 8.76827161e-04,  9.03579902e-04,  
+        8.39818807e-04, 5.19936424e-04,  1.46415537e-04,  3.29054769e-05,
+        7.30096312e-05,  6.47553400e-05,  2.18817382e-05, 4.47955432e-06,  
+        7.34705616e-06,  9.06184045e-06, 9.45269846e-06,  1.00464939e-05,  
+        8.28166011e-06, 7.09361681e-06,  6.43221021e-06,  3.12425880e-06,
+        2.57495214e-07,  6.48560373e-07,  1.88421498e-07
+        ]
 
-def test_radialI(coords, imager, dataset, generic_parameters):
-    # obtain a 1d radial brightness profile I(r) from an image
-    nchan = dataset.nchan
-    model = SimpleNet(coords=coords, nchan=nchan)
+    np.testing.assert_allclose(Vtest, expected, rtol=1e-6,
+                               err_msg="test_radialV")    
\ No newline at end of file

From 611fee00a14d81a247bc29d19fa9f6e121ba5522 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 14:47:54 -0400
Subject: [PATCH 082/109] test_radialV: add plot

---
 test/onedim_test.py | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/test/onedim_test.py b/test/onedim_test.py
index 51ef90f0..a8f90cba 100644
--- a/test/onedim_test.py
+++ b/test/onedim_test.py
@@ -66,6 +66,18 @@ def test_radialV(mock_1d_vis_model, tmp_path):
 
     qtest, Vtest = radialV(Vtrue, coords, geom, rescale_flux=True, bins=None)
 
+    fig, ax = plt.subplots()
+
+    ax.plot(q_dep / 1e6, Vtrue_dep, 'k.', label='truth deprojected')
+    ax.plot(qtest / 1e6, Vtest, 'r.-', label='recovery')
+
+    ax.set_xlabel(r'Baseline [M$\lambda$]')
+    ax.set_ylabel('Re(V) [Jy]')
+
+    ax.title('Geometry: {:}'.format(geom))
+
+    fig.savefig(tmp_path / "test_radialV.png", dpi=300)
+    plt.close("all")
 
     expected = [
         2.53998336e-01,  1.59897580e-01,  8.59460326e-02,  7.42189236e-02,  

From 6a6894d5556e63b01cc3ded8405f7fccaf3f679d Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 17:39:46 -0400
Subject: [PATCH 083/109] mock_1d_archive: update zenodo record ref

---
 test/conftest.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/conftest.py b/test/conftest.py
index 0f991f13..4048b36e 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -120,7 +120,7 @@ def dataset_cont(mock_visibility_data_cont, coords):
 def mock_1d_archive():
     # use astropy routines to cache data
     fname = download_file(
-        "https://zenodo.org/record/{:d}/files/mock_disk_1d.npz".format(zenodo_version),
+        f"https://zenodo.org/record/{zenodo_record}/files/mock_disk_1d.npz",
         cache=True,
         pkgname="mpol",
     )

From 59ad2c7321f282fa427e3e5cc2cd8b34612c034e Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 18:08:28 -0400
Subject: [PATCH 084/109] test_radialI: consistent naming

---
 test/onedim_test.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/test/onedim_test.py b/test/onedim_test.py
index a8f90cba..8fc7f717 100644
--- a/test/onedim_test.py
+++ b/test/onedim_test.py
@@ -6,19 +6,19 @@
 
 def test_radialI(mock_1d_image_model, tmp_path):
     # obtain a 1d radial brightness profile I(r) from an image    
-    
-    r, i, i2d, _, _, geom, coords = mock_1d_image_model
+
+    rtrue, itrue, i2dtrue, _, _, geom, coords = mock_1d_image_model
 
     bins = np.linspace(0, 2.0, 100)
-    rtest, itest = radialI(i2d, coords, geom, bins=bins)
+    rtest, itest = radialI(i2dtrue, coords, geom, bins=bins)
 
     fig, ax = plt.subplots(ncols=2)
 
-    plot_image(i2d, extent=coords.img_ext, ax=ax[0], clab='Jy / sr')
+    plot_image(i2dtrue, extent=coords.img_ext, ax=ax[0], clab='Jy / sr')
 
     ax[0].title('AS 209-like profile.\nGeometry: {:}'.format(geom))
 
-    ax[1].plot(r, i, 'k', label='truth')
+    ax[1].plot(rtrue, itrue, 'k', label='truth')
     ax[1].plot(rtest, itest, 'r.-', label='recovery')
     
     ax[1].set_xlabel('r [arcsec]')

From f94c235353af889cdc4293e04c6c5c4f05916a9b Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 18:33:40 -0400
Subject: [PATCH 085/109] onedim_test: set geom

---
 test/onedim_test.py | 20 +++++++++++++-------
 1 file changed, 13 insertions(+), 7 deletions(-)

diff --git a/test/onedim_test.py b/test/onedim_test.py
index 8fc7f717..7c171041 100644
--- a/test/onedim_test.py
+++ b/test/onedim_test.py
@@ -7,20 +7,23 @@
 def test_radialI(mock_1d_image_model, tmp_path):
     # obtain a 1d radial brightness profile I(r) from an image    
 
-    rtrue, itrue, i2dtrue, _, _, geom, coords = mock_1d_image_model
+    rtrue, itrue, i2dtrue, _, _, _, coords = mock_1d_image_model
+
+    geom = {'incl': 40.0, 'Omega': 130.0, 'omega': 0.0, 'dRA': 0.3, 'dDec': -0.2}
 
     bins = np.linspace(0, 2.0, 100)
+
     rtest, itest = radialI(i2dtrue, coords, geom, bins=bins)
 
     fig, ax = plt.subplots(ncols=2)
 
     plot_image(i2dtrue, extent=coords.img_ext, ax=ax[0], clab='Jy / sr')
 
-    ax[0].title('AS 209-like profile.\nGeometry: {:}'.format(geom))
-
     ax[1].plot(rtrue, itrue, 'k', label='truth')
     ax[1].plot(rtest, itest, 'r.-', label='recovery')
     
+    ax[0].set_title(f"AS 209-like mock disk.\nGeometry {geom}")
+    
     ax[1].set_xlabel('r [arcsec]')
     ax[1].set_ylabel('I [Jy / sr]')
 
@@ -62,9 +65,13 @@ def test_radialI(mock_1d_image_model, tmp_path):
 def test_radialV(mock_1d_vis_model, tmp_path):
     # obtain a 1d radial visibility profile V(q) from 2d visibilities
 
-    Vtrue, Vtrue_dep, q_dep, geom, coords = mock_1d_vis_model
+    Vtrue, Vtrue_dep, q_dep, _, coords = mock_1d_vis_model
 
-    qtest, Vtest = radialV(Vtrue, coords, geom, rescale_flux=True, bins=None)
+    geom = {'incl': 40.0, 'Omega': 130.0, 'omega': 0.0, 'dRA': 0.3, 'dDec': -0.2}
+
+    bins = np.linspace(1,5e3,100)
+
+    qtest, Vtest = radialV(Vtrue, coords, geom, rescale_flux=True, bins=bins)
 
     fig, ax = plt.subplots()
 
@@ -73,8 +80,7 @@ def test_radialV(mock_1d_vis_model, tmp_path):
 
     ax.set_xlabel(r'Baseline [M$\lambda$]')
     ax.set_ylabel('Re(V) [Jy]')
-
-    ax.title('Geometry: {:}'.format(geom))
+    ax.set_title(f"Geometry {geom}")
 
     fig.savefig(tmp_path / "test_radialV.png", dpi=300)
     plt.close("all")

From 45e7b2310c43366a36dc4248beb33282d788c84c Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 18:46:24 -0400
Subject: [PATCH 086/109] add frank as test dep

---
 setup.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/setup.py b/setup.py
index bc5c2c4a..e4da0143 100644
--- a/setup.py
+++ b/setup.py
@@ -32,7 +32,8 @@ def get_version(rel_path):
         "requests",
         "astropy",
         "tensorboard",
-        "mypy",
+        "mypy"
+        "frank>=1.2.1",
     ],
     "docs": [
         "sphinx>=2.3.0",

From 3a9f6b1896bf635af4dbbcba6c4c493d4a1e036c Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 18:46:40 -0400
Subject: [PATCH 087/109] tests.yml: add frank to dep installs

---
 .github/workflows/tests.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 97a81b7c..958faa67 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -17,6 +17,7 @@ jobs:
         # in __init__. below we'll reinstall for the tests.
         run: |
           pip install astropy
+          pip install frank
           pip install .
       - name: Cache/Restore the .mpol folder cache
         uses: actions/cache@v3

From 35caea81a8264a4ab94e2ae633e8d2fe82b76b95 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 18:53:39 -0400
Subject: [PATCH 088/109] extra_requires: remove frank version restriction

---
 setup.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/setup.py b/setup.py
index e4da0143..b8dde7a2 100644
--- a/setup.py
+++ b/setup.py
@@ -33,7 +33,7 @@ def get_version(rel_path):
         "astropy",
         "tensorboard",
         "mypy"
-        "frank>=1.2.1",
+        "frank",
     ],
     "docs": [
         "sphinx>=2.3.0",
@@ -57,7 +57,7 @@ def get_version(rel_path):
         "arviz[all]"
     ],
     "analysis": [
-        "frank>=1.2.1",
+        "frank",
     ],
 }
 

From c2c6d21ba83f4c9bb652d2fb29f9c2de6994db7b Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 19:00:35 -0400
Subject: [PATCH 089/109] typo

---
 setup.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/setup.py b/setup.py
index b8dde7a2..d8c07d67 100644
--- a/setup.py
+++ b/setup.py
@@ -32,8 +32,8 @@ def get_version(rel_path):
         "requests",
         "astropy",
         "tensorboard",
-        "mypy"
-        "frank",
+        "mypy",
+        "frank>=1.2.1",
     ],
     "docs": [
         "sphinx>=2.3.0",
@@ -57,7 +57,7 @@ def get_version(rel_path):
         "arviz[all]"
     ],
     "analysis": [
-        "frank",
+        "frank>=1.2.1",
     ],
 }
 

From afad0a97c113ad395f95324de890d44e3b7cd9bf Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 20:41:55 -0400
Subject: [PATCH 090/109] test_radialI: clean up diag plot

---
 test/onedim_test.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/test/onedim_test.py b/test/onedim_test.py
index 7c171041..abef8519 100644
--- a/test/onedim_test.py
+++ b/test/onedim_test.py
@@ -15,17 +15,18 @@ def test_radialI(mock_1d_image_model, tmp_path):
 
     rtest, itest = radialI(i2dtrue, coords, geom, bins=bins)
 
-    fig, ax = plt.subplots(ncols=2)
+    fig, ax = plt.subplots(ncols=2, figsize=(10,5))
 
     plot_image(i2dtrue, extent=coords.img_ext, ax=ax[0], clab='Jy / sr')
 
     ax[1].plot(rtrue, itrue, 'k', label='truth')
     ax[1].plot(rtest, itest, 'r.-', label='recovery')
     
-    ax[0].set_title(f"AS 209-like mock disk.\nGeometry {geom}")
+    ax[0].set_title(f"Geometry:\n{geom}", fontsize=7)
     
     ax[1].set_xlabel('r [arcsec]')
     ax[1].set_ylabel('I [Jy / sr]')
+    ax[1].legend()
 
     fig.savefig(tmp_path / "test_radialI.png", dpi=300)
     plt.close("all")

From 55029c7225f5aa28e4509b26a145a45d617246ee Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 20:42:05 -0400
Subject: [PATCH 091/109] test_radialV: clean up diag plot

---
 test/onedim_test.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/test/onedim_test.py b/test/onedim_test.py
index abef8519..476c2144 100644
--- a/test/onedim_test.py
+++ b/test/onedim_test.py
@@ -77,11 +77,13 @@ def test_radialV(mock_1d_vis_model, tmp_path):
     fig, ax = plt.subplots()
 
     ax.plot(q_dep / 1e6, Vtrue_dep, 'k.', label='truth deprojected')
-    ax.plot(qtest / 1e6, Vtest, 'r.-', label='recovery')
+    ax.plot(qtest / 1e3, Vtest, 'r.-', label='recovery')
 
+    ax.set_xlim(-0.5, 6)
     ax.set_xlabel(r'Baseline [M$\lambda$]')
     ax.set_ylabel('Re(V) [Jy]')
-    ax.set_title(f"Geometry {geom}")
+    ax.set_title(f"Geometry {geom}", fontsize=10)
+    ax.legend()
 
     fig.savefig(tmp_path / "test_radialV.png", dpi=300)
     plt.close("all")

From 53da8871469446706b02399268bee949f71c2f5b Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Wed, 1 Nov 2023 23:00:54 -0400
Subject: [PATCH 092/109] remove unneeded func

---
 src/mpol/images.py | 38 +-------------------------------------
 1 file changed, 1 insertion(+), 37 deletions(-)

diff --git a/src/mpol/images.py b/src/mpol/images.py
index b0f9b2dc..5090fd6a 100644
--- a/src/mpol/images.py
+++ b/src/mpol/images.py
@@ -315,40 +315,4 @@ def to_FITS(self, fname="cube.fits", overwrite=False, header_kwargs=None):
         hdul = fits.HDUList([hdu])
         hdul.writeto(fname, overwrite=overwrite)
 
-        hdul.close()
-
-
-def np_to_imagecube(image, coords, nchan=1, wrap=False):
-    """Convenience function for converting a numpy image into an MPoL ImageCube 
-    tensor (see mpol.images.ImageCube)
-
-    Parameters
-    ----------
-    image : array
-        An image in numpy format
-    coords : `mpol.coordinates.GridCoords` object
-        Instance of the `mpol.coordinates.GridCoords` class
-    nchan : int, default=1
-        Number of channels in the image. Default assumes a single 2D image
-    wrap : bool, default=False 
-        Whether to wrap the numpy image so that index 0 is in the image center
-        (FFT algorithms typically place index 0 in the image corner)
-
-    Returns
-    -------
-    icube : `mpol.images.ImageCube` object
-        The image cube tensor
-    """
-    if wrap: 
-        # move the 0 index to the image center 
-        image = utils.center_np_image(image)
-
-    # broadcast image to (nchan, npix, npix)
-    img_packed_cube = np.broadcast_to(image, 
-                                      (nchan, coords.npix, coords.npix)).copy()
-
-    # convert to pytorch tensor
-    img_packed_tensor = torch.from_numpy(img_packed_cube)
-
-    # insert into ImageCube layer
-    return ImageCube(coords=coords, nchan=nchan, cube=img_packed_tensor)
+        hdul.close()
\ No newline at end of file

From ea9e97010f9db7565a5ebbc7f6a3485dd49c60c9 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 2 Nov 2023 22:01:54 -0400
Subject: [PATCH 093/109] conftest: recast 'geom' type

---
 test/conftest.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/test/conftest.py b/test/conftest.py
index 4048b36e..cad9d9f7 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -1,8 +1,9 @@
 import numpy as np
 import pytest
+import torch
 from astropy.utils.data import download_file
 
-from mpol import coordinates, gridding
+from mpol import coordinates, fourier, gridding, images
 from mpol.__init__ import zenodo_record
 
 # We need a fixture which provides mock visibilities of the sort we'd
@@ -136,6 +137,7 @@ def mock_1d_image_model(mock_1d_archive):
     i2dtrue = m['i2dtrue']
     xmax = ymax = m['xmax']
     geom = m['geometry']
+    geom = geom[()]
 
     coords = coordinates.GridCoords(cell_size=xmax * 2 / i2dtrue.shape[0], 
                                     npix=i2dtrue.shape[0])
@@ -150,6 +152,7 @@ def mock_1d_vis_model(mock_1d_archive):
     Vtrue_dep = m['vis_dep']
     q_dep = m['baselines_dep']
     geom = m['geometry']
+    geom = geom[()]
 
     xmax = m['xmax']
     i2dtrue = m['i2dtrue']

From e4d0499879d778090bbe18eda6071af1a2b63b7a Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 2 Nov 2023 22:02:39 -0400
Subject: [PATCH 094/109] mock_1d_image_model: undo image centering

---
 test/conftest.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/test/conftest.py b/test/conftest.py
index cad9d9f7..a39749ef 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -141,8 +141,11 @@ def mock_1d_image_model(mock_1d_archive):
 
     coords = coordinates.GridCoords(cell_size=xmax * 2 / i2dtrue.shape[0], 
                                     npix=i2dtrue.shape[0])
+    
+    # the center of the array is already at the center of the image -->
+    # undo this as expected by input to ImageCube
+    i2dtrue = np.flip(np.fft.fftshift(i2dtrue), 1)
 
-    return rtrue, itrue, i2dtrue, xmax, ymax, geom, coords
 
 
 @pytest.fixture

From 90f8fdc806145ad09d95b0d90d2707e5334a1024 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 2 Nov 2023 22:03:00 -0400
Subject: [PATCH 095/109] mock_1d_image_model: place image in cube

---
 test/conftest.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/test/conftest.py b/test/conftest.py
index a39749ef..097d42dc 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -146,6 +146,11 @@ def mock_1d_image_model(mock_1d_archive):
     # undo this as expected by input to ImageCube
     i2dtrue = np.flip(np.fft.fftshift(i2dtrue), 1)
 
+    # pack the numpy image array into an ImageCube
+    packed_cube = np.broadcast_to(i2dtrue, (1, coords.npix, coords.npix)).copy()
+    packed_tensor = torch.from_numpy(packed_cube)
+    cube_true = images.ImageCube(coords=coords, nchan=1, cube=packed_tensor)
+
 
 
 @pytest.fixture

From bdd599ddfd076593dac4991f517295b449b79c80 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 2 Nov 2023 22:03:13 -0400
Subject: [PATCH 096/109] mock_1d_image_model: update outputs

---
 test/conftest.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/conftest.py b/test/conftest.py
index 097d42dc..6eeeb865 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -151,6 +151,7 @@ def mock_1d_image_model(mock_1d_archive):
     packed_tensor = torch.from_numpy(packed_cube)
     cube_true = images.ImageCube(coords=coords, nchan=1, cube=packed_tensor)
 
+    return rtrue, itrue, cube_true, xmax, ymax, geom
 
 
 @pytest.fixture

From b50e98e1ba3f11187a9d116a70aa1ef554341a93 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 2 Nov 2023 22:03:37 -0400
Subject: [PATCH 097/109] mock_1d_vis_model: create fourier cube

---
 test/conftest.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/test/conftest.py b/test/conftest.py
index 6eeeb865..c77cdc94 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -169,7 +169,11 @@ def mock_1d_vis_model(mock_1d_archive):
     coords = coordinates.GridCoords(cell_size=xmax * 2 / i2dtrue.shape[0],
                                     npix=i2dtrue.shape[0])
 
-    return Vtrue, Vtrue_dep, q_dep, geom, coords
+    # create a FourierCube
+    packed_cube = np.broadcast_to(Vtrue, (1, len(Vtrue))).copy()
+    packed_tensor = torch.from_numpy(packed_cube)
+    cube_true = fourier.FourierCube(coords=coords)
+
 
 
 

From c313f10bdf83e97278a7d33ee45d8f33d389f8b5 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 2 Nov 2023 22:03:57 -0400
Subject: [PATCH 098/109] mock_1d_vis_model: put vis in ground_cube

---
 test/conftest.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/test/conftest.py b/test/conftest.py
index c77cdc94..416b5bee 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -174,6 +174,9 @@ def mock_1d_vis_model(mock_1d_archive):
     packed_tensor = torch.from_numpy(packed_cube)
     cube_true = fourier.FourierCube(coords=coords)
 
+    # insert the vis tensor into the FourierCube ('vis' would typically be 
+    # populated by taking the FFT of an image)
+    cube_true.ground_cube = packed_tensor
 
 
 

From 2209863198264e8fc001e53998493b3fc9b7ed8b Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 2 Nov 2023 22:04:07 -0400
Subject: [PATCH 099/109] mock_1d_vis_model: update outputs

---
 test/conftest.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/test/conftest.py b/test/conftest.py
index 416b5bee..39eee5e6 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -178,6 +178,8 @@ def mock_1d_vis_model(mock_1d_archive):
     # populated by taking the FFT of an image)
     cube_true.ground_cube = packed_tensor
 
+    return cube_true, Vtrue_dep, q_dep, geom
+
 
 
 @pytest.fixture

From 7bfb5c61026a0e719d7f9191d1e18447741131e0 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 2 Nov 2023 22:04:37 -0400
Subject: [PATCH 100/109] radialI: change input from array to icube

---
 src/mpol/onedim.py | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index e8f4a91c..b4ac71af 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -1,15 +1,14 @@
 import numpy as np 
+from mpol.utils import torch2npy
 
-def radialI(image, coords, geom, bins=None):
+def radialI(icube, geom, chan=0, bins=None):
     r"""
-    Obtain a 1D (radial) brightness profile I(r) from an image.
+    Obtain a 1D (radial) brightness profile I(r) from an image cube.
 
     Parameters
     ----------
-    image : array
-        2D image array 
-    coords : `mpol.coordinates.GridCoords` object
-        Instance of the `mpol.coordinates.GridCoords` class
+    icube : `mpol.images.ImageCube` object
+        Instance of the MPoL `images.ImageCube` class
     geom : dict 
         Dictionary of source geometry. Keys:
             "incl" : float, unit=[deg]

From 24b48cd84b129aea493b6d037b814e4e15477d4c Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 2 Nov 2023 22:04:54 -0400
Subject: [PATCH 101/109] radialI: add chan arg

---
 src/mpol/onedim.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index b4ac71af..c25a58ec 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -21,6 +21,8 @@ def radialI(icube, geom, chan=0, bins=None):
                 Phase center offset in right ascension. Positive is west of north.
             "dDec" : float, unit=[arcsec]
                 Phase center offset in declination.
+    chan : int, default=0
+        Channel of the image cube corresponding to the desired image
     bins : array, default=None, unit=[arcsec]
         Radial bin edges to use in calculating I(r). If None, bins will span 
         the full image, with widths equal to the hypotenuse of the pixels

From 845bd28742942d5db3c2af352c45eed181214f08 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 2 Nov 2023 22:05:14 -0400
Subject: [PATCH 102/109] radialI: update calls to icube

---
 src/mpol/onedim.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index c25a58ec..932fabb0 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -36,7 +36,7 @@ def radialI(icube, geom, chan=0, bins=None):
     """
 
     # projected Cartesian pixel coordinates [arcsec]
-    xx, yy = coords.sky_x_centers_2D, coords.sky_y_centers_2D
+    xx, yy = icube.coords.sky_x_centers_2D, icube.coords.sky_y_centers_2D
 
     # shift image center to source center
     xc, yc = xx - geom["dRA"], yy - geom["dDec"]
@@ -53,13 +53,15 @@ def radialI(icube, geom, chan=0, bins=None):
 
     if bins is None:
         # choose sensible bin size and range
-        step = np.hypot(coords.cell_size, coords.cell_size)
+        step = np.hypot(icube.coords.cell_size, icube.coords.cell_size)
         bins = np.arange(0.0, np.max((abs(xc.ravel()), abs(yc.ravel()))), step)
 
     bin_counts, bin_edges = np.histogram(a=rr, bins=bins, weights=None)
 
     # cumulative binned brightness in each annulus
-    Is, _ = np.histogram(a=rr, bins=bins, weights=np.ravel(image))
+    Is, _ = np.histogram(a=rr, bins=bins, 
+                         weights=torch2npy(icube.sky_cube[chan]).ravel()
+                         )
 
     # mask empty bins
     mask = (bin_counts == 0)

From cf835cfbee880b0b84db22afdd4c9a05c8c4f3cf Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 2 Nov 2023 22:05:27 -0400
Subject: [PATCH 103/109] radialV: change input to fcube

---
 src/mpol/onedim.py | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 932fabb0..c207e33c 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -75,17 +75,15 @@ def radialI(icube, geom, chan=0, bins=None):
     return bin_centers, Is
 
 
-def radialV(V, coords, geom, rescale_flux, bins=None):
+def radialV(fcube, geom, rescale_flux, chan=0, bins=None):
     r"""
     Obtain the 1D (radial) visibility model V(q) corresponding to a 2D MPoL 
     image. 
 
     Parameters
     ----------
-    V : array
-        2D visibility amplitudes
-    coords : `mpol.coordinates.GridCoords` object
-        Instance of the `mpol.coordinates.GridCoords` class
+    fcube : `~mpol.fourier.FourierCube` object
+        Instance of the MPoL `fourier.FourierCube` class
     geom : dict 
         Dictionary of source geometry. Keys:
             "incl" : float, unit=[deg]

From 943874e598b045278bfe21c35abdf71eda1f5540 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 2 Nov 2023 22:05:37 -0400
Subject: [PATCH 104/109] radialV: add chan arg

---
 src/mpol/onedim.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index c207e33c..51650c68 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -103,6 +103,8 @@ def radialV(fcube, geom, rescale_flux, chan=0, bins=None):
         The source's integrated (2D) flux is assumed to be:
             :math:`F = \cos(i) \int_r^{r=R}{I(r) 2 \pi r dr}`.
         No rescaling would be appropriate in the optically thin limit. 
+    chan : int, default=0
+        Channel of the image cube corresponding to the desired image        
     bins : array, default=None, unit=[k\lambda]
         Baseline bin edges to use in calculating V(q). If None, bins will span 
         the model baseline distribution, with widths equal to the hypotenuse of 
@@ -110,7 +112,7 @@ def radialV(fcube, geom, rescale_flux, chan=0, bins=None):
 
     Returns
     -------
-    q : array, unit=:math:[`k\lambda`]
+    bin_centers : array, unit=:math:[`k\lambda`]
         Baselines corresponding to `u` and `v`
     Vs : array, unit=[Jy] 
         Visibility amplitudes at `q`

From 052b4c820f95bb6f1456f51013685e459f1d58a1 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 2 Nov 2023 22:06:08 -0400
Subject: [PATCH 105/109] radialV: use fcube

---
 src/mpol/onedim.py | 24 +++++++++++++++---------
 1 file changed, 15 insertions(+), 9 deletions(-)

diff --git a/src/mpol/onedim.py b/src/mpol/onedim.py
index 51650c68..58e83131 100644
--- a/src/mpol/onedim.py
+++ b/src/mpol/onedim.py
@@ -121,15 +121,21 @@ def radialV(fcube, geom, rescale_flux, chan=0, bins=None):
     -----
     This routine requires the `frank <https://github.com/discsim/frank>`_ package
     """
+    from frank.geometry import apply_phase_shift, deproject
 
-    # projected model (u,v) coordinates [k\lambda]
-    uu, vv = coords.sky_u_centers_2D, coords.sky_v_centers_2D
+    # projected model (u,v) points [k\lambda]
+    uu, vv = fcube.coords.sky_u_centers_2D, fcube.coords.sky_v_centers_2D
 
-    from frank.geometry import apply_phase_shift, deproject
-    # phase-shift the model visibilities
-    Vp = apply_phase_shift(uu.ravel() * 1e3, vv.ravel() * 1e3, V.ravel(), geom["dRA"], geom["dDec"], inverse=True)
-    # deproject the model (u,v) points
-    up, vp, _ = deproject(uu.ravel() * 1e3, vv.ravel() * 1e3, geom["incl"], geom["Omega"])
+    # visibilities 
+    V = torch2npy(fcube.ground_cube[chan]).ravel()
+
+    # phase-shift the visibilities
+    Vp = apply_phase_shift(uu.ravel() * 1e3, vv.ravel() * 1e3, V, geom["dRA"], 
+                           geom["dDec"], inverse=True)
+    
+    # deproject the (u,v) points
+    up, vp, _ = deproject(uu.ravel() * 1e3, vv.ravel() * 1e3, geom["incl"], 
+                          geom["Omega"])
 
     # if the source is optically thick, rescale the deprojected V(q)
     if rescale_flux: 
@@ -139,12 +145,12 @@ def radialV(fcube, geom, rescale_flux, chan=0, bins=None):
     up /= 1e3
     vp /= 1e3
 
-    # deprojected baseline coordinates
+    # deprojected baselines
     qq = np.hypot(up, vp) 
 
     if bins is None:
         # choose sensible bin size and range
-        step = np.hypot(coords.du, coords.dv) / 2
+        step = np.hypot(fcube.coords.du, fcube.coords.dv) / 2
         bins = np.arange(0.0, max(qq), step)
 
     bin_counts, bin_edges = np.histogram(a=qq, bins=bins, weights=None)

From deeec2a4ffab84799990daeae56908a993f1ed96 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 2 Nov 2023 22:06:42 -0400
Subject: [PATCH 106/109] test_radialI: update data retrieval call

---
 test/onedim_test.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/test/onedim_test.py b/test/onedim_test.py
index 476c2144..fa25379f 100644
--- a/test/onedim_test.py
+++ b/test/onedim_test.py
@@ -7,9 +7,7 @@
 def test_radialI(mock_1d_image_model, tmp_path):
     # obtain a 1d radial brightness profile I(r) from an image    
 
-    rtrue, itrue, i2dtrue, _, _, _, coords = mock_1d_image_model
-
-    geom = {'incl': 40.0, 'Omega': 130.0, 'omega': 0.0, 'dRA': 0.3, 'dDec': -0.2}
+    rtrue, itrue, icube, _, _, geom = mock_1d_image_model
 
     bins = np.linspace(0, 2.0, 100)
 

From 34930687bb868c88328dc68182ca73e1360cc531 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 2 Nov 2023 22:06:56 -0400
Subject: [PATCH 107/109] test_radialV: update data retrieval call

---
 test/onedim_test.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/test/onedim_test.py b/test/onedim_test.py
index fa25379f..aa583626 100644
--- a/test/onedim_test.py
+++ b/test/onedim_test.py
@@ -64,9 +64,7 @@ def test_radialI(mock_1d_image_model, tmp_path):
 def test_radialV(mock_1d_vis_model, tmp_path):
     # obtain a 1d radial visibility profile V(q) from 2d visibilities
 
-    Vtrue, Vtrue_dep, q_dep, _, coords = mock_1d_vis_model
-
-    geom = {'incl': 40.0, 'Omega': 130.0, 'omega': 0.0, 'dRA': 0.3, 'dDec': -0.2}
+    fcube, Vtrue_dep, q_dep, geom = mock_1d_vis_model
 
     bins = np.linspace(1,5e3,100)
 

From 5bd168a4cc61179307ceaadd05d5efc29c8708ae Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 2 Nov 2023 22:07:28 -0400
Subject: [PATCH 108/109] test_radialI: use icube

---
 test/onedim_test.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/test/onedim_test.py b/test/onedim_test.py
index aa583626..0240ea12 100644
--- a/test/onedim_test.py
+++ b/test/onedim_test.py
@@ -3,6 +3,7 @@
 
 from mpol.onedim import radialI, radialV
 from mpol.plot import plot_image
+from mpol.utils import torch2npy
 
 def test_radialI(mock_1d_image_model, tmp_path):
     # obtain a 1d radial brightness profile I(r) from an image    
@@ -11,11 +12,12 @@ def test_radialI(mock_1d_image_model, tmp_path):
 
     bins = np.linspace(0, 2.0, 100)
 
-    rtest, itest = radialI(i2dtrue, coords, geom, bins=bins)
+    rtest, itest = radialI(icube, geom, bins=bins)
 
     fig, ax = plt.subplots(ncols=2, figsize=(10,5))
 
-    plot_image(i2dtrue, extent=coords.img_ext, ax=ax[0], clab='Jy / sr')
+    plot_image(np.squeeze(torch2npy(icube.sky_cube)), extent=icube.coords.img_ext, 
+               ax=ax[0], clab='Jy / sr')
 
     ax[1].plot(rtrue, itrue, 'k', label='truth')
     ax[1].plot(rtest, itest, 'r.-', label='recovery')

From fb7bf32911d9a7d8d937e3f97f0e188781059a23 Mon Sep 17 00:00:00 2001
From: Jeff Jennings <jjennings1519@gmail.com>
Date: Thu, 2 Nov 2023 22:07:43 -0400
Subject: [PATCH 109/109] test_radialV: use fcube

---
 test/onedim_test.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/onedim_test.py b/test/onedim_test.py
index 0240ea12..b78cc5ba 100644
--- a/test/onedim_test.py
+++ b/test/onedim_test.py
@@ -70,7 +70,7 @@ def test_radialV(mock_1d_vis_model, tmp_path):
 
     bins = np.linspace(1,5e3,100)
 
-    qtest, Vtest = radialV(Vtrue, coords, geom, rescale_flux=True, bins=bins)
+    qtest, Vtest = radialV(fcube, geom, rescale_flux=True, bins=bins)
 
     fig, ax = plt.subplots()