Merge pull request #103 from punch-mission/develop

add imax effect, bump codecov, update citation

.github/workflows/CI.yml

@@ -33,7 +33,7 @@ jobs:
         pip install .
         pytest --cov=solpolpy tests/ --cov-report xml:/home/runner/coverage.xml
     - name: Upload coverage to Codecov
-      uses: codecov/codecov-action@v3
+      uses: codecov/codecov-action@v4
       env:
         CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
       with:

README.md

@@ -31,7 +31,7 @@ Please open a discussion or issue for help.
 - [ ] quantification and propagation of error
 - [ ] additional plotting utilities
 - [ ] more comprehensive support for 4-polarizer systems
-- [ ] functions to deal with the IMAX effect in wide-field imagers
+- [x] functions to deal with the IMAX effect in wide-field imagers
 
 ## Contributing
 We encourage all contributions.

docs/source/cite.rst

@@ -1,4 +1,5 @@
 Cite
 =======
 
-To cite the software please use [the citation from Zenodo](https://zenodo.org/records/10289143) for the version you used.
+To cite the software please use `the citation from Zenodo <https://zenodo.org/records/10289143>`_
+for the version you used.

pyproject.toml

@@ -1,11 +1,12 @@
 [project]
 name = "solpolpy"
-version = "0.1.1"
+version = "0.1.2"
 authors = [
     { name="J. Marcus Hughes", email="[email protected]"},
     { name="Matthew J. West", email="[email protected]"},
     { name="Ritesh Patel", email="[email protected]"},
     { name="Bryce M. Walbridge", email="[email protected]" },
+    { name="Chris Lowder", email="[email protected]"}
 ]
 description = "Solar polarization resolver for any instrument"
 readme = "README.md"

solpolpy/core.py

@@ -1,7 +1,9 @@
 """Core transformation functions for solpolpy"""
 import typing as t
 
+import astropy.units as u
 import networkx as nx
+import numpy as np
 from ndcube import NDCollection, NDCube
 
 from solpolpy.alpha import radial_north
@@ -12,7 +14,7 @@
 from solpolpy.polarizers import npol_to_mzp
 
 
-def resolve(input_data: t.Union[t.List[str], NDCollection], out_system: str) -> NDCollection:
+def resolve(input_data: t.Union[t.List[str], NDCollection], out_system: str, imax_effect: bool = False) -> NDCollection:
     """
     Apply - apply a polarization transformation to a set of input
     dataframes.
@@ -36,6 +38,11 @@ def resolve(input_data: t.Union[t.List[str], NDCollection], out_system: str) ->
         - "fourpol": For observations taken at sequence of four polarizer angles, i.e. 0°, 45°, 90° and 135°.
         - "npol": Set of images taken at than three polarizing angles other than MZP
 
+    imax_effect : Boolean
+        The 'IMAX effect' describes the change in apparent measured polarization angle as an result of foreshortening effects.
+        This effect becomes more pronounced for wide field polarized imagers - see Patel et al (2024, in preparation)
+        If True, applies the IMAX effect for wide field imagers as part of the resolution process.
+
     Raises
     ------
     AssertionError
@@ -66,6 +73,12 @@ def resolve(input_data: t.Union[t.List[str], NDCollection], out_system: str) ->
     equation = get_transform_equation(transform_path)
     requires_alpha = check_alpha_requirement(transform_path)
 
+    if imax_effect:
+        if (input_kind == 'MZP') and (out_system == 'MZP'):
+            input_data = resolve_imax_effect(input_data)
+        else:
+            raise UnsupportedTransformationError('IMAX effect applies only for MZP->MZP solpolpy transformations')
+
     if requires_alpha and "alpha" not in input_key:
         input_data = add_alpha(input_data)
     result = equation(input_data)
@@ -158,6 +171,83 @@ def check_alpha_requirement(path: t.List[str]) -> bool:
     return requires_alpha
 
 
+def generate_imax_matrix(arrayshape) -> np.ndarray:
+    """
+    Define an A matrix with which to convert MZP^ (camera coords) = A x MZP (solar coords)
+
+    Parameters
+    -------
+    arrayshape
+        Defined input WCS array shape for matrix generation
+
+    Returns
+    -------
+    ndarray
+        Output A matrix used in converting between camera coordinates and solar coordinates
+
+    """
+
+    # Ideal MZP wrt Solar North
+    thmzp = [-60, 0, 60] * u.degree
+
+    long_arr, lat_arr = np.meshgrid(np.linspace(-20, 20, arrayshape[0]), np.linspace(-20, 20, arrayshape[1]))
+
+    # Foreshortening (IMAX) effect on polarizer angle
+    phi_m = np.arctan2(np.tan(thmzp[0]) * np.cos(long_arr * u.degree), np.cos(lat_arr * u.degree)).to(u.degree)
+    phi_z = np.arctan2(np.tan(thmzp[1]) * np.cos(long_arr * u.degree), np.cos(lat_arr * u.degree)).to(u.degree)
+    phi_p = np.arctan2(np.tan(thmzp[2]) * np.cos(long_arr * u.degree), np.cos(lat_arr * u.degree)).to(u.degree)
+
+    phi = np.stack([phi_m, phi_z, phi_p])
+
+    # Define the A matrix
+    mat_a = np.empty((arrayshape[0], arrayshape[1], 3, 3))
+
+    for i in range(3):
+        for j in range(3):
+            mat_a[:, :, i, j] = (4 * np.cos(phi[i] - thmzp[j]) ** 2 - 1) / 3
+
+    return mat_a
+
+
+def resolve_imax_effect(input_data: NDCollection) -> NDCollection:
+    """
+    Resolves the IMAX effect for provided input data, correcting measured polarization angles for wide FOV imagers.
+
+    Parameters
+    -------
+    input_data : NDCollection
+        Input data on which to correct foreshortened polarization angles
+
+    Returns
+    -------
+    NDCollection
+        Output data with corrected polarization angles
+
+    """
+
+    data_shape = input_data['Bm'].data.shape
+    data_mzp_camera = np.zeros([data_shape[0], data_shape[1], 3, 1])
+
+    for i, key in enumerate(['Bm', 'Bz', 'Bp']):
+        data_mzp_camera[:, :, i, 0] = input_data[key].data
+
+    imax_matrix = generate_imax_matrix(data_shape)
+    try:
+        imax_matrix_inv = np.linalg.inv(imax_matrix)
+    except np.linalg.LinAlgError as err:
+        if 'Singular matrix' in str(err):
+            raise ValueError('Singular IMAX effect matrix is degenerate')
+        else:
+            raise err
+
+    data_mzp_solar = np.matmul(imax_matrix_inv, data_mzp_camera)
+
+    for i, key in enumerate(input_data.keys()):
+        input_data[key].data[:, :] = data_mzp_solar[:, :, i, 0]
+
+    return input_data
+
+
 def _determine_image_shape(input_collection: NDCollection) -> t.Tuple[int, int]:
     """Evaluates the shape of the image in the input NDCollection
 

tests/test_core.py

@@ -11,6 +11,7 @@
 from pytest import fixture
 
 from solpolpy.core import _determine_image_shape, add_alpha, determine_input_kind, resolve
+from solpolpy.errors import UnsupportedTransformationError
 
 wcs = astropy.wcs.WCS(naxis=3)
 wcs.ctype = 'WAVE', 'HPLT-TAN', 'HPLN-TAN'
@@ -43,6 +44,22 @@ def mzp_ones():
     data_out.append(("Bz", NDCube(np.array([1]), wcs=wcs, meta={'POLAR': 0})))
     data_out.append(("Bm", NDCube(np.array([1]), wcs=wcs, meta={'POLAR': -60})))
     return NDCollection(data_out, meta={}, aligned_axes="all")
+
+@fixture
+def mzp_data():
+    data_out = []
+    data_out.append(("Bp", NDCube(np.random.random([50,50]), wcs=wcs, meta={'POLAR': 60})))
+    data_out.append(("Bz", NDCube(np.random.random([50,50]), wcs=wcs, meta={'POLAR': 0})))
+    data_out.append(("Bm", NDCube(np.random.random([50,50]), wcs=wcs, meta={'POLAR': -60})))
+    return NDCollection(data_out, meta={}, aligned_axes="all")
+
+@fixture
+def bpb_data():
+    data_out = []
+    data_out.append(("B", NDCube(np.random.random([50,50]), wcs=wcs, meta={'POLAR': 'B'})))
+    data_out.append(("pB", NDCube(np.random.random([50,50]), wcs=wcs, meta={'POLAR': 'pB'})))
+    return NDCollection(data_out, meta={}, aligned_axes="all")
+
 @fixture
 def mzp_ones_alpha():
     data_out = []
@@ -260,3 +277,22 @@ def test_btbr_to_mzp(btbr_ones):
 def test_bp3_to_bthp(bp3_ones):
     result = resolve(bp3_ones, "Bthp")
     assert isinstance(result, NDCollection)
+
+
+def test_imax_effect(mzp_data):
+    result = resolve(mzp_data, "MZP", imax_effect=True)
+    assert isinstance(result, NDCollection)
+    for key in result.keys():
+        assert np.sum(result[key].data * mzp_data[key].data) != 0
+
+
+def test_imax_effect_unsupported_transformation_output(mzp_data):
+    with pytest.raises(UnsupportedTransformationError):
+        result = resolve(mzp_data, "BpB", imax_effect=True)
+        assert isinstance(result, NDCollection)
+
+
+def test_imax_effect_unsupported_transformation_input(bpb_data):
+    with pytest.raises(UnsupportedTransformationError):
+        result = resolve(bpb_data, "MZP", imax_effect=True)
+        assert isinstance(result, NDCollection)