Merge pull request #162 from stscieisenhamer/rcal-853-velocity

RCAL-853: Allow ability to introduct velocity aberration scale factor

docs/romanisim/l1.rst

@@ -47,7 +47,26 @@ Gain
 ----
 L1 files are in units of DN.  We convert from photons to DN using gains from CRDS or a default values of 2 electron per DN, treating electrons and photons as equivalent.  This neglects the possibility that the quantum yield may be different from 1 at certain wavelengths, in particular in the blue, where `Givans+2022 <https://iopscience.iop.org/article/10.1088/1538-3873/ac46ba>`_ find a quantum field of up to 1.05.
 
+.. _velocity-aberration:
+
+Velocity Aberration
+-------------------
+
+Velocity (or stellar) aberration is the amount of change in an object's apparent
+location on the sky due to motion of the observer in relation to the object. For
+observers orbiting the Sun at approximately 1 AU, the absolute correction can be
+up to 26 arcsec. For Roman, this absolute difference will be handled by
+spacecraft operations. However, this is true for only the target; the aberration
+changes across the FOV, varying in both position and time. This differential
+aberration is implemented as a scale factor from the reference point of each
+detector. For Roman, this scale is approximately one-fifth of an arcsec at the
+corners.
+
+For simulation, since the Roman orbit is undefined and would be dominated by the
+Earth's motion anyways, Earth barycentric velocity is used.
+
 .. automodapi:: romanisim.l1
 .. automodapi:: romanisim.nonlinearity
 .. automodapi:: romanisim.persistence
-.. automodapi:: romanisim.cr 
+.. automodapi:: romanisim.cr
+.. automodapi:: romanisim.velocity_aberration

docs/romanisim/running.rst

@@ -15,12 +15,7 @@ The ``romanisim-make-image`` command line interface (CLI) has a number of argume
 this functionality::
 
     romanisim-make-image -h
-    usage: romanisim-make-image [-h] [--bandpass BANDPASS] [--boresight] [--catalog CATALOG] [--config CONFIG]
-				[--date DATE] [--level LEVEL] [--ma_table_number MA_TABLE_NUMBER] [--nobj NOBJ]
-				[--previous PREVIOUS] [--radec RADEC RADEC] [--rng_seed RNG_SEED] [--roll ROLL]
-				[--sca SCA] [--usecrds] [--webbpsf] [--truncate TRUNCATE]
-				[--pretend-spectral PRETEND_SPECTRAL] [--drop-extra-dq]
-				filename
+    usage: romanisim-make-image [-h] [--bandpass BANDPASS] [--boresight] [--catalog CATALOG] [--config CONFIG] [--date DATE] [--level LEVEL] [--ma_table_number MA_TABLE_NUMBER] [--nobj NOBJ] [--previous PREVIOUS] [--radec RADEC RADEC] [--rng_seed RNG_SEED] [--roll ROLL] [--sca SCA] [--usecrds] [--webbpsf] [--truncate TRUNCATE] [--pretend-spectral PRETEND_SPECTRAL] [--drop-extra-dq] [--scale-factor SCALE_FACTOR] filename
 
     Make a demo image.
 
@@ -37,21 +32,20 @@ this functionality::
       --level LEVEL         1 or 2, for L1 or L2 output (default: 2)
       --ma_table_number MA_TABLE_NUMBER
       --nobj NOBJ
-      --previous PREVIOUS   previous simulated file in chronological order used for persistence modeling. (default:
-			    None)
+      --previous PREVIOUS   previous simulated file in chronological order used for persistence modeling. (default: None)
       --radec RADEC RADEC   ra and dec (deg) (default: None)
       --rng_seed RNG_SEED
-      --roll ROLL           Position angle (North towards YIdl) measured at the V2Ref/V3Ref of the aperture used.
-			    (default: 0)
-      --sca SCA             SCA to simulate. Use -1 to generate images for all SCAs; include {} in filename for this
-			    mode to indicate where the SCA number should be filled, e.g. l1_wfi{}.asdf (default: 7)
+      --roll ROLL           Position angle (North towards YIdl) measured at the V2Ref/V3Ref of the aperture used. (default: 0)
+      --sca SCA             SCA to simulate. Use -1 to generate images for all SCAs; include {} in filename for this mode to indicate where the SCA number should be filled, e.g.
+                            l1_wfi{}.asdf (default: 7)
       --usecrds             Use CRDS for distortion map (default: False)
       --webbpsf             Use webbpsf for PSF (default: False)
       --truncate TRUNCATE   If set, truncate the MA table at given number of resultants. (default: None)
       --pretend-spectral PRETEND_SPECTRAL
-			    Pretend the image is spectral. exposure.type and instrument.element are updated to be
-			    grism / prism. (default: None)
+                            Pretend the image is spectral. exposure.type and instrument.element are updated to be grism / prism. (default: None)
       --drop-extra-dq       Do not store the optional simulated dq array. (default: False)
+      --scale-factor SCALE_FACTOR
+                            Velocity aberration-induced scale factor. If negative, use given time to calculated based on orbit ephemeris. (default: -1.0)
 
     EXAMPLE: romanisim-make-image output_image.asdf
 

docs/romanisim/wcs.rst

@@ -5,5 +5,28 @@ The simulator has two options for modeling distortion and world coordinate syste
 
 The latter system works by grabbing reference distortion maps for the appropriate detector and filter combinations from the Roman CRDS server.  These distortion maps are then wrapped into a galsim WCS object and fed to galsim's rendering pipeline.
 
+Velocity Aberration
+-------------------
+
+When using a distortion model, either provided or through CRDS, there is also
+opportunity to apply velocity aberration: See :ref:`velocity-aberration` for a
+general discussion. As discussed, the aberration can introduce an up-to
+one-fifth of an arcsec change across the FOV for any particular detector.
+
+The scale factor can be specified, or the simulator can calculate one. Though
+the exact orbit of Roman is unknown, it will be orbiting around the L2 point, as
+Webb does. The predominant velocity is very similar to that of Earth's at this
+point, so Earth is used, using a scale factor increase 1.01 to approximate the
+L2 orbit.
+
+In the simulator, there is an additional uncertainty introduced by including
+velocity aberration correction. On-orbit, the scale is calculated from the
+reference point of each detector. However, the simulator does not calculate a
+separate WCS for each detector, but computes a single WCS for the FOV of all the
+detectors, centered either on the boresight or the center of the full WFI,
+WFI_CEN. As a result, the velocity aberration can introduce ~10^-6 arcsec
+effect. Below shows a typical error vector across a detector.
+
+.. image:: ff_va_WFI02.png
 
 .. automodapi:: romanisim.wcs

pyproject.toml

@@ -25,7 +25,8 @@ dependencies = [
     # "rad >=0.19.2",
     # "roman_datamodels >=0.19.1",
     "rad @ git+https://github.com/spacetelescope/rad.git",
-    "roman_datamodels @ git+https://github.com/spacetelescope/roman_datamodels.git",
+    # "roman_datamodels @ git+https://github.com/spacetelescope/roman_datamodels.git",
+    "roman_datamodels @ git+https://github.com/stscieisenhamer/roman_datamodels.git@rcal-853-velocity",
     "gwcs >=0.19.0",
     "jsonschema >=4.8",
     "numpy >=1.22",

romanisim/image.py

@@ -879,10 +879,6 @@ def make_asdf(slope, slopevar_rn, slopevar_poisson, metadata=None,
     # ref_file: conceptually sound when we work from crds reference
     #     files
     # target: can push forward from APT file
-    # velocity_aberration:  I guess to first order,
-    #     ignore the detailed orbit around L2 and just project
-    #     the earth out to L2, and use that for the velocity
-    #     aberration?  Don't do until someone asks.
     # visit: start_time, end_time, total_exposures, ...?
     # wcsinfo: v2_ref, v3_ref, vparity, v3yangle, ra_ref, dec_ref
     #     roll_ref, s_region

romanisim/parameters.py

@@ -161,6 +161,14 @@
                    'detector': 'WFI07',
                    'optical_element': 'F184',
                    },
+    'ephemeris': {'time': Time('2026-01-01').mjd,
+                  'spatial_x': 0.,
+                  'spatial_y': 0.,
+                  'spatial_z': 0.,
+                  'velocity_x': 0.,
+                  'velocity_y': 0.,
+                  'velocity_z': 0.
+                  },
     'exposure': {'start_time': Time('2026-01-01T00:00:00'),
                  'type': 'WFI_IMAGE',
                  'ma_table_number': 4,
@@ -171,7 +179,10 @@
                  'target_dec': 66.0,
                  'target_aperture': 'WFI_CEN',
                  },
-    'wcsinfo': {'ra_ref': 270.0,
+    'velocity_aberration': {'scale_factor': 1.0},
+    'wcsinfo': {'aperture_name': 'WFI_CEN',
+                'pa_aperture': 0.,
+                'ra_ref': 270.0,
                 'dec_ref': 66.0,
                 'v2_ref': 0,
                 'v3_ref': 0,

romanisim/ris_make_utils.py

@@ -10,12 +10,14 @@
 from astropy import table
 from astropy import time
 from astropy import coordinates
+from astropy import units as u
 import galsim
 from galsim import roman
 import roman_datamodels
 from roman_datamodels import stnode
 from romanisim import catalog, image, wcs
 from romanisim import parameters, log
+from romanisim.util import calc_scale_factor
 import romanisim
 
 NMAP = {'apt': 'http://www.stsci.edu/Roman/APT'}
@@ -46,7 +48,7 @@ def merge_nested_dicts(dict1, dict2):
 
 
 def set_metadata(meta=None, date=None, bandpass='F087', sca=7,
-                 ma_table_number=4, truncate=None):
+                 ma_table_number=4, truncate=None, scale_factor=1.0):
     """
     Set / Update metadata parameters
 
@@ -62,6 +64,8 @@ def set_metadata(meta=None, date=None, bandpass='F087', sca=7,
         Integer identifier of the detector to simulate (starting at 1)
     ma_table_number : int
         Integer specifying which MA Table entry to use
+    scale_factor : float
+        Velocity aberration-induced scale factor
 
     Returns
     -------
@@ -90,6 +94,23 @@ def set_metadata(meta=None, date=None, bandpass='F087', sca=7,
     else:
         meta['exposure']['truncated'] = False
 
+    # Velocity aberration
+    if scale_factor <= 0.:
+        scale_factor = calc_scale_factor(meta['exposure']['start_time'], meta['wcsinfo']['ra_ref'], meta['wcsinfo']['dec_ref'])
+    meta['velocity_aberration']['scale_factor'] = scale_factor
+
+    # Fill out some ephemeris information, presuming all is earth.
+    position, velocity = coordinates.get_body_barycentric_posvel('earth', meta['exposure']['start_time'])
+    position = position.xyz.to(u.km)
+    velocity = velocity.xyz.to(u.km / u.s)
+    meta['ephemeris']['time'] = meta['exposure']['start_time'].mjd
+    meta['ephemeris']['spatial_x'] = position.value[0]
+    meta['ephemeris']['spatial_y'] = position.value[1]
+    meta['ephemeris']['spatial_z'] = position.value[2]
+    meta['ephemeris']['velocity_x'] = velocity.value[0]
+    meta['ephemeris']['velocity_y'] = velocity.value[1]
+    meta['ephemeris']['velocity_z'] = velocity.value[2]
+
     return meta
 
 

romanisim/tests/test_wcs.py

@@ -98,7 +98,7 @@ def test_wcs():
     # should be close to the reference v2 & v3 offset.
     assert np.abs(cc3.separation(cc4).to(u.arcsec).value
                   - np.hypot(*parameters.v2v3_wficen)) < 1
-    metadata = dict(pointing=dict(), instrument=dict(), wcsinfo=dict())
+    metadata = dict(pointing=dict(), instrument=dict(), wcsinfo=dict(), velocity_aberration=dict())
     metadata['instrument']['detector'] = 'WFI01'
     util.update_pointing_and_wcsinfo_metadata(metadata, wcs.GWCS(gwcs))
     assert metadata['wcsinfo']['aperture_name'] == 'WFI01_FULL'
@@ -149,3 +149,25 @@ def test_wcs_crds_match():
 
     # Expect a GWCS object as opposed to a dictionary
     assert type(twcs) is wcs.GWCS
+
+
+def test_scale_factor():
+    """Test that specifying a scale factor actually calculates a new wcs"""
+
+    # Truth aiming for.
+    cc = SkyCoord(ra=0 * u.deg, dec=0 * u.deg)
+    scale_factor = 1.1
+    truth = SkyCoord([(0.        , 0.        ), (0.03361111, 0.0336111 ),
+                      (0.0672222 , 0.06722216), (0.10083325, 0.10083312),
+                      (0.13444424, 0.13444393)],
+                     unit='deg')
+
+    # Make a simple grid of pixel coordinates to calculate sky for
+    grid = range(0, 4096, 1000)
+
+    # Create the wcs and generate results
+    distortion = make_fake_distortion_function()
+    gwcs = wcs.make_wcs(cc, distortion, scale_factor=scale_factor)
+    sky = gwcs.pixel_to_world(grid, grid)
+
+    assert all(truth.separation(sky).to(u.arcsec).value < 1e-3)

romanisim/util.py

@@ -2,15 +2,33 @@
 """
 
 import numpy as np
-from astropy.coordinates import SkyCoord
+from astropy.coordinates import SkyCoord, get_body_barycentric_posvel
 from astropy import units as u
 from astropy.time import Time
 import galsim
 import gwcs as gwcsmod
 
 from romanisim import parameters, wcs, bandpass
+from romanisim.velocity_aberration import compute_va_effects
 from scipy import integrate
 
+__all__ = ["skycoord",
+           "celestialcoord",
+           "scalergb",
+           "random_points_in_cap",
+           "random_points_in_king",
+           "random_points_at_radii",
+           "add_more_metadata",
+           "update_pointing_and_wcsinfo_metadata",
+           "king_profile",
+           "sample_king_distances",
+           "decode_context_times",
+           "default_image_meta",
+           "update_photom_keywords",
+           "merge_dicts",
+           "calc_scale_factor",
+]
+
 
 def skycoord(celestial):
     """Turn a CelestialCoord into a SkyCoord.
@@ -252,7 +270,7 @@ def update_pointing_and_wcsinfo_metadata(metadata, gwcs):
     v2v3 = distortion(*center)
     radec = gwcs(*center)
     t2sky = gwcs.get_transform('v2v3', 'world')
-    radecn = t2sky(v2v3[0], v2v3[1] + 1)
+    radecn = t2sky(v2v3[0], v2v3[1] + 100)
     roll_ref = (
         SkyCoord(radec[0] * u.deg, radec[1] * u.deg).position_angle(
         SkyCoord(radecn[0] * u.deg, radecn[1] * u.deg)))
@@ -282,6 +300,10 @@ def update_pointing_and_wcsinfo_metadata(metadata, gwcs):
     pa_v3 = pa_v3.to(u.deg).value
     metadata['pointing']['pa_v3'] = pa_v3
 
+    # Update velocity aberration meta for the reference point
+    metadata['velocity_aberration']['ra_reference'] = radec[0]
+    metadata['velocity_aberration']['dec_reference'] = radec[1]
+
 
 def king_profile(r, rc, rt):
     """Compute the King profile.
@@ -539,3 +561,34 @@ def merge_dicts(a, b):
         else:
             a[key] = b[key]
     return a
+
+
+def calc_scale_factor(date, ra, dec):
+    """Calculate velocity aberration scale factor
+
+    The L2 orbit is just a delta on the Earth's orbit. At the moment, there is no ephemeris for
+    Roman yet, the Earth's barycentric velocity is used to calculate velocity aberration. A scale
+    of 1.01 is applied to the velocity to scale, approximately, to what would be Roman's velocity.
+
+    Parameters
+    ----------
+    date : str or astropy.Time
+        The date at which to calculate the velocity.
+
+    ra, dec: float
+        The right ascension and declination of the target (or some other
+        point, such as the center of a detector) in the barycentric coordinate
+        system.  The equator and equinox should be the same as the coordinate
+        system for the velocity. In degrees.
+
+    Returns
+    -------
+    scale_factor : float
+        The velocity aberration scale factor
+    """
+    _, velocity = get_body_barycentric_posvel('earth', date)
+    velocity = 1.01 * velocity  # Move from earth to Roman.
+    xyz_velocity = velocity.xyz.to(u.km / u.s)
+    scale_factor, _, _ = compute_va_effects(*xyz_velocity.value, ra, dec)
+
+    return scale_factor