Merge pull request #77 from Johannes-Sahlmann/exact-idl_to_tel

Exact idl to tel

pysiaf/iando/read.py

@@ -286,6 +286,7 @@ def read_hst_fgs_amudotrep(file=None, version=None):
         'n_cones': re.compile(r'NUMBER OF CONES:   (?P<n_cones>\d)'),
         'date': re.compile(r'(?P<day>[ 123][0-9])-(?P<month>[A-Z][A-Z][A-Z])-(?P<year>[0-9][0-9])'),
         'cone_vector': re.compile(r'(CONE)*(CONE VECTOR)*(CONE ANGLE)'),
+        'cone_vector_tel': re.compile(r'(CONE)*(REVISED CONE VECTOR)*(PREVIOUS CONE VECTOR)'),
         'tvs': re.compile(r'(FGS TO ST TRANSFORMATION MATRICES)'),
     }
 
@@ -310,7 +311,12 @@ def read_hst_fgs_amudotrep(file=None, version=None):
                                    data_start=astropy_table_index+2, data_end=astropy_table_index + 2 + n_cones,
                                    guess=False, names=('CONE', 'X', 'Y', 'Z', 'CONE_ANGLE_DEG'))
                 # table.pprint()
-                data[fgs_id]['cone_parameters'] = table
+                data[fgs_id]['cone_parameters_fgs'] = table
+            elif key == 'cone_vector_tel':
+                table = Table.read(file, format='ascii.no_header', delimiter=' ',
+                                   data_start=astropy_table_index+2, data_end=astropy_table_index + 2 + n_cones,
+                                   guess=False, names=('CONE', 'V1', 'V2', 'V3', 'V1_PREV', 'V2_PREV', 'V3_PREV'))
+                data[fgs_id]['cone_parameters_tel'] = table
             elif key == 'tvs':
                 table = Table.read(file, format='ascii.no_header', delimiter=' ',
                                    data_start=astropy_table_index+2, data_end=astropy_table_index+2+3,

pysiaf/tests/test_aperture.py

@@ -24,31 +24,74 @@ def siaf_objects():
     return siafs
 
 
-def test_idl_to_tel():
-    """Test the transformations between ideal and telescope frames that do not use the planar approximation."""
+def test_idl_to_tel(verbose=False):
+    """Test the transformations between ideal and telescope frames."""
     siaf = Siaf('NIRISS')
 
     x_idl, y_idl = get_grid_coordinates(10, (0, 0), 100)
 
-    verbose = False
-
     for aper_name in siaf.apertures.keys():
-        # aperture
         aperture = siaf[aper_name]
 
-        for input_coordinates in ['spherical', 'tangent_plane']:
-            v2, v3 = aperture.idl_to_tel(x_idl, y_idl, method='spherical_transformation', input_coordinates=input_coordinates)
-            x_idl_2, y_idl_2 = aperture.tel_to_idl(v2, v3, method='spherical_transformation', output_coordinates=input_coordinates)
-            x_diff = np.abs(x_idl - x_idl_2)
-            y_diff = np.abs(y_idl - y_idl_2)
-            if verbose:
-                print('Aperture {} {} x_diff {} y_diff {}'.format(aper_name, input_coordinates, np.max(x_diff), np.max(y_diff)))
-            if input_coordinates == 'spherical':
-                threshold = 1e-12
-            elif input_coordinates == 'tangent_plane':
-                threshold = 5e-10
-            assert np.max(x_diff) < threshold
-            assert np.max(y_diff) < threshold
+        for idl_to_tel_method in ['planar_approximation', 'spherical']:
+            if idl_to_tel_method == 'spherical':
+                input_coordinate_types = ['polar', 'cartesian']
+            else:
+                input_coordinate_types = ['tangent_plane']
+
+            for input_coordinates in input_coordinate_types:
+                v2, v3 = aperture.idl_to_tel(x_idl, y_idl, method=idl_to_tel_method, input_coordinates=input_coordinates, output_coordinates=input_coordinates)
+                x_idl_2, y_idl_2 = aperture.tel_to_idl(v2, v3, method=idl_to_tel_method, input_coordinates=input_coordinates, output_coordinates=input_coordinates)
+                x_diff = np.abs(x_idl - x_idl_2)
+                y_diff = np.abs(y_idl - y_idl_2)
+                if verbose:
+                    print('{} {}: Aperture {} {} x_diff {} y_diff {}'.format(idl_to_tel_method, input_coordinates, aper_name, input_coordinates, np.max(x_diff), np.max(y_diff)))
+                if idl_to_tel_method == 'planar_approximation':
+                    threshold = 7e-14
+                elif idl_to_tel_method == 'spherical':
+                    if input_coordinates == 'polar':
+                        threshold = 6e-13
+                    elif input_coordinates == 'cartesian':
+                        threshold = 5e-8
+                assert np.max(x_diff) < threshold
+                assert np.max(y_diff) < threshold
+
+
+def test_hst_fgs_idl_to_tel(verbose=False):
+    """Test the transformations between ideal and telescope frames."""
+
+    siaf = Siaf('HST')
+
+    x_idl, y_idl = get_grid_coordinates(2, (0, -50), 1000, y_width=400)
+
+    for aper_name in 'FGS1 FGS2 FGS3'.split():
+        aperture = siaf[aper_name]
+        for idl_to_tel_method in ['planar_approximation', 'spherical']:
+            if idl_to_tel_method == 'spherical':
+                input_coordinate_types = ['polar', 'cartesian']
+            else:
+                input_coordinate_types = ['tangent_plane']
+
+            for input_coordinates in input_coordinate_types:
+                if input_coordinates == 'polar':
+                    v2, v3 = aperture.idl_to_tel(x_idl, y_idl, method=idl_to_tel_method,
+                                                 input_coordinates='cartesian',
+                                                 output_coordinates=input_coordinates)
+                    x_idl_2, y_idl_2 = aperture.tel_to_idl(v2, v3, method=idl_to_tel_method,
+                                                           input_coordinates=input_coordinates,
+                                                           output_coordinates='cartesian')
+                else:
+                    v2, v3 = aperture.idl_to_tel(x_idl, y_idl, method=idl_to_tel_method, input_coordinates=input_coordinates, output_coordinates=input_coordinates)
+                    x_idl_2, y_idl_2 = aperture.tel_to_idl(v2, v3, method=idl_to_tel_method, input_coordinates=input_coordinates, output_coordinates=input_coordinates)
+                x_diff = np.abs(x_idl - x_idl_2)
+                y_diff = np.abs(y_idl - y_idl_2)
+                if verbose:
+                    print('{} {}: Aperture {} {} x_diff {} y_diff {}'.format(idl_to_tel_method, input_coordinates, aper_name, input_coordinates, np.max(x_diff), np.max(y_diff)))
+
+                threshold = 2.5e-13
+
+                assert np.max(x_diff) < threshold
+                assert np.max(y_diff) < threshold
 
 
 def test_jwst_aperture_transforms(siaf_objects, verbose=False, threshold=None):

pysiaf/tests/test_projection.py

@@ -31,7 +31,8 @@ def make_grid_coordinates(centre_deg=(0., 0.)):
 
 def test_tangent_plane_projection_roundtrip(grid_coordinates):
     """Transform from RA/Dec to tangent plane and back. Check that input is recovered."""
-    centre_deg = (80., -70.)   # setting this to 0,0 fails because of 0,360 deg wrapping
+    # centre_deg = (80., -70.)
+    centre_deg = (0., 0.)
     ra_deg, dec_deg = grid_coordinates(centre_deg=centre_deg)
 
     # project to detector pixel coordinates
@@ -47,7 +48,8 @@ def test_tangent_plane_projection_roundtrip(grid_coordinates):
 
 def test_project_to_tangent_plane(grid_coordinates):
     """Compare projection code built with astropy functions with independent implementation."""
-    centre_deg = (80., -70.)  # setting this to 0,0 fails because of 0,360 deg wrapping
+    # centre_deg = (80., -70.)
+    centre_deg = (0., 0.)
     ra_deg, dec_deg = grid_coordinates(centre_deg=centre_deg)
 
     x_1, y_1 = projection.project_to_tangent_plane(ra_deg, dec_deg, centre_deg[0], centre_deg[1])
@@ -60,7 +62,8 @@ def test_project_to_tangent_plane(grid_coordinates):
 
 def test_deproject_from_tangent_plane(grid_coordinates):
     """Compare projection code built with astropy functions with independent implementation."""
-    centre_deg = (80., -70.)  # setting this to 0,0 fails because of 0,360 deg wrapping
+    # centre_deg = (80., -70.)
+    centre_deg = (0., 0.)
     ra_deg, dec_deg = grid_coordinates(centre_deg=centre_deg)
 
     x, y = projection.project_to_tangent_plane(ra_deg, dec_deg, centre_deg[0], centre_deg[1])
@@ -144,9 +147,9 @@ def tangent_plane_deprojection(e, n, alpha_ref, delta_ref):
                         - np.sin(rho)*np.sin(d0)*np.cos(B))
     alpha = alpha_ref + np.rad2deg(dalpha)
 
-    if np.any(alpha < 0.0):
-        index = np.where(alpha < 0.0)
-        alpha[index] += 360.0    # set to range 0 to 360 deg
+    # if np.any(alpha < 0.0):
+    #     index = np.where(alpha < 0.0)
+    #     alpha[index] += 360.0    # set to range 0 to 360 deg
 
     delta_rad = np.arcsin(np.sin(d0)*np.cos(rho) + np.cos(d0)*np.sin(rho)*np.cos(B))
     delta = np.rad2deg(delta_rad)

pysiaf/tests/test_rotations.py

@@ -4,11 +4,15 @@
 Authors
 -------
     Colin Cox
+    Johannes Sahlmann
 
 """
-import numpy as np
 from math import sin, cos, acos, pi
-from pysiaf.utils import rotations as rt
+import numpy as np
+
+import astropy.units as u
+
+from pysiaf.utils import rotations, tools
 
 # Some values used in both tests
 #  Set up some arbitrary values
@@ -24,6 +28,13 @@
 v2_array = np.array([200.0, -500.0, 300.0, -400.0])
 v3_array = np.array([300.0, -600.0, -400.0, 500.0])
 
+ra_deg = 20.0
+dec_deg = 70.0
+pa_deg = 15.0
+v2_arcsec = 200.0
+v3_arcsec = 300.0
+
+
 
 def test_attitude(verbose=False):
     """Test the properties of the attitude matrix in calculating positions and roll angles.
@@ -37,36 +48,36 @@ def test_attitude(verbose=False):
         set to true only if detailed print-out needed
 
     """
-    a = rt.attitude(v2, v3, ra, dec, roll)
+    a = rotations.attitude(v2, v3, ra, dec, roll)
     if verbose:
         print('attitude\n', a)
 
     #  Show that attitude matrix correctly connects given inputs in both directions
-    (ra_test1, dec_test1) = rt.pointing(a, v2, v3)
+    (ra_test1, dec_test1) = rotations.pointing(a, v2, v3)
     if verbose:
         print('RA and Dec %10.6f %10.6f %10.3e %10.3e ' % (ra_test1, dec_test1, ra_test1-ra, dec_test1 - dec))
     assert abs(ra_test1 - ra) < 1.0e-10 and abs(dec_test1 - dec) < 1.0e-10, 'Miscalculated RA or Dec'
 
-    (v2_test1, v3_test1) = rt.getv2v3(a, ra, dec)
+    (v2_test1, v3_test1) = rotations.getv2v3(a, ra, dec)
     if verbose:
         print('V2 V3 %10.6f %10.6f %10.3e %10.3e' % (v2_test1, v3_test1, v2_test1 - v2, v3_test1 - v3))
     assert abs(v2_test1 - v2) < 1.0e-10 and abs(v3_test1 - v3) < 1.0e-10, 'Miscalculated V2 or V3'
 
     #  Show that points away from defining position are correctly handled
-    (v2_test2, v3_test2) = rt.getv2v3(a, ra2, dec2)
-    (ra_test2, dec_test2) = rt.pointing(a, v2_test2, v3_test2)
+    (v2_test2, v3_test2) = rotations.getv2v3(a, ra2, dec2)
+    (ra_test2, dec_test2) = rotations.pointing(a, v2_test2, v3_test2)
     if verbose:
         print('Test 2 %10.6f %10.6f' % (v2_test2, v3_test2))
         print('Test2 RA and Dec %10.6f %10.6f' % (ra_test2, dec_test2))
         print('Test2 %10.3e %10.3e' % (ra_test2 - ra2, dec_test2 - dec2))
     assert abs(ra_test2 - ra2) < 1.0e10 and abs(dec_test2 - dec2) < 1.0e-10, 'Miscalculated RA or Dec'
 
     #  Position angles at reference point
-    pa1 = rt.posangle(a, v2, v3)
-    pa2 = rt.sky_posangle(a, ra, dec)
+    pa1 = rotations.posangle(a, v2, v3)
+    pa2 = rotations.sky_posangle(a, ra, dec)
     #  and at displaced point
-    pa3 = rt.posangle(a, v2_test2, v3_test2)
-    pa4 = rt.sky_posangle(a, ra_test2, dec_test2)
+    pa3 = rotations.posangle(a, v2_test2, v3_test2)
+    pa4 = rotations.sky_posangle(a, ra_test2, dec_test2)
     if verbose:
         print('PA tests')
         print('%10.6f %10.6f %10.6f %10.3e' % (roll, pa1, pa2, pa1 - pa2))
@@ -75,16 +86,16 @@ def test_attitude(verbose=False):
     assert abs(pa3 - pa4) < 1.0e-10, 'Disagreement for displaced point position angles'
 
     # Test some functions with arrays as input
-    unit_vectors = rt.unit(ra_array, dec_array)
-    rd = rt.radec(unit_vectors, positive_ra=True)
-    rd_test = rt.pointing(a, v2_array, v3_array)
+    unit_vectors = rotations.unit(ra_array, dec_array)
+    rd = rotations.radec(unit_vectors, positive_ra=True)
+    rd_test = rotations.pointing(a, v2_array, v3_array)
     if verbose:
         print(unit_vectors)
         print('RD\n', rd)
         print(rd_test[0])
         print(rd_test[1])
 
-    v_test = rt.getv2v3(a, ra_array, dec_array)
+    v_test = rotations.getv2v3(a, ra_array, dec_array)
     if verbose:
         for i in range(len(ra_array)):
             print(v_test[0][i], v_test[1][i])
@@ -94,6 +105,59 @@ def test_attitude(verbose=False):
     assert abs(v_test[1][0] - v3_array[0] < 1.0e-10), 'V3 values do not match'
 
 
+def test_attitude_matrix():
+    """Compare original and new attitude matrix generator functions."""
+
+    ra = ra_deg * u.deg
+    dec = dec_deg * u.deg
+    pa = pa_deg * u.deg
+    v2 = v2_arcsec * u.arcsec
+    v3 = v3_arcsec * u.arcsec
+
+    attitude = rotations.attitude(v2_arcsec, v3_arcsec, ra_deg, dec_deg, pa_deg)
+    attitude_matrix = rotations.attitude_matrix(v2, v3, ra, dec, pa)
+
+    assert np.all(attitude==attitude_matrix)
+
+
+def test_sky_to_tel():
+    """Test application of the attitude matrix"""
+
+    # test with quantities
+    ra = ra_deg * u.deg
+    dec = dec_deg * u.deg
+    pa = pa_deg * u.deg
+    v2 = v2_arcsec * u.arcsec
+    v3 = v3_arcsec * u.arcsec
+
+    attitude = rotations.attitude_matrix(v2, v3, ra, dec, pa)
+    ra_2, dec_2 = rotations.tel_to_sky(attitude, *rotations.sky_to_tel(attitude, ra, dec))
+    assert np.abs((ra - ra_2).to(u.milliarcsecond).value) < 1e-6
+    assert np.abs((dec - dec_2).to(u.milliarcsecond).value) < 1e-6
+
+    # test without quantities
+    attitude = rotations.attitude_matrix(v2_arcsec, v3_arcsec, ra_deg, dec_deg, pa_deg)
+    ra_2, dec_2 = rotations.tel_to_sky(attitude, *rotations.sky_to_tel(attitude, ra_deg, dec_deg))
+    assert np.abs(ra_deg - ra_2.to(u.deg).value)*u.deg.to(u.milliarcsecond) < 1e-6
+    assert np.abs(dec_deg - dec_2.to(u.deg).value)*u.deg.to(u.milliarcsecond) < 1e-6
+
+    # test array inputs
+    n_side = 3
+    span = 2 * u.arcmin
+    x_width = span.to(u.deg).value
+    centre_deg = (ra_deg, dec_deg)
+    ra_array_deg, dec_array_deg = tools.get_grid_coordinates(n_side, centre_deg, x_width)
+
+    ra_array_2, dec_array_2 = rotations.tel_to_sky(attitude, *rotations.sky_to_tel(attitude, ra_array_deg*u.deg, dec_array_deg*u.deg))
+    assert np.all(np.abs(ra_array_deg*u.deg - ra_array_2) < 1e-6 * u.milliarcsecond)
+    assert np.all(np.abs(dec_array_deg*u.deg - dec_array_2) < 1e-6 * u.milliarcsecond)
+
+    ra_array_2, dec_array_2 = rotations.tel_to_sky(attitude, *rotations.sky_to_tel(attitude, ra_array_deg, dec_array_deg))
+    assert np.all(np.abs(ra_array_deg*u.deg - ra_array_2) < 1e-6 * u.milliarcsecond)
+    assert np.all(np.abs(dec_array_deg*u.deg - dec_array_2) < 1e-6 * u.milliarcsecond)
+
+
+
 def test_axial_rotation(verbose=False):
     """Compare vector transformation using the attitude matrix with a single rotation about an axis.
 
@@ -111,13 +175,13 @@ def test_axial_rotation(verbose=False):
     values = np.random.rand(2)
     alpha = 2*pi*values[0]
     beta = acos(2*values[1] - 1.0)
-    u = np.array([cos(alpha)*cos(beta), sin(alpha)*cos(beta), sin(beta)])
+    vector = np.array([cos(alpha)*cos(beta), sin(alpha)*cos(beta), sin(beta)])
 
-    a = rt.attitude(v2, v3, ra, dec, roll)
-    va = np.dot(a, u)    # Transform using attitude matrix
+    a = rotations.attitude(v2, v3, ra, dec, roll)
+    va = np.dot(a, vector)    # Transform using attitude matrix
 
-    (axis, phi, quaternion) = rt.rodrigues(a)  # obtain single rotation parameters equivalent to attitude matrix
-    vb = rt.axial_rotation(axis, phi, u)  # Transform using axial rotation
+    (axis, phi, quaternion) = rotations.rodrigues(a)  # obtain single rotation parameters equivalent to attitude matrix
+    vb = rotations.axial_rotation(axis, phi, vector)  # Transform using axial rotation
     dot_product = np.dot(va, vb)
 
     if verbose:
@@ -132,3 +196,25 @@ def test_axial_rotation(verbose=False):
         print('Dot product', dot_product)
 
     assert abs(dot_product - 1.0) < 1.0e-12, 'Transforms do not agree'
+
+
+def test_unit_vector_from_cartesian():
+    """Test unit vector construction."""
+
+    # scalar inputs
+    x = 0.1
+    y = 0.2
+    unit_vector = rotations.unit_vector_from_cartesian(x=x, y=y)
+    assert (np.linalg.norm(unit_vector) - 1) < 1e-14
+
+    # scalar inputs with unit
+    x = 200 * u.arcsec
+    y = -300 * u.arcsec
+    unit_vector = rotations.unit_vector_from_cartesian(x=x, y=y)
+    assert (np.linalg.norm(unit_vector) - 1) < 1e-14
+
+    # array inputs with unit
+    x = np.linspace(-100, 100, 10) * u.arcsec
+    y = np.linspace(-500, -100, 10) * u.arcsec
+    unit_vector = rotations.unit_vector_from_cartesian(x=x, y=y)
+    assert np.all(np.abs(np.linalg.norm(unit_vector, axis=0) - 1)) < 1e-14