Merge pull request #496 from mperrin/field_dep_wl

Add field dependence into weak lens model
spacetelescope · Nov 18, 2021 · b32536c · b32536c
2 parents de249ec + 0865072
commit b32536c
Show file tree

Hide file tree

Showing 4 changed files with 227 additions and 73 deletions.
diff --git a/webbpsf/optics.py b/webbpsf/optics.py
@@ -1172,6 +1172,30 @@ def _fix_zgrid_NaNs(xgrid, ygrid, zgrid, rot_ang=0):
     return zgrid
 
 
+def _get_initial_pupil_sampling(instrument):
+    """Utility function to retrieve the sampling of the first plane in some optical system.
+
+    Returns: npix, pixelscale
+    """
+    # Determine the pupil sampling of the first aperture in the
+    # instrument's optical system
+    if isinstance(instrument.pupil, poppy.OpticalElement):
+        # This branch needed to handle the OTE Linear Model case
+        npix = instrument.pupil.shape[0]
+        pixelscale = instrument.pupil.pixelscale
+    else:
+        # these branches to handle FITS files, by name or as an object
+        if isinstance(instrument.pupil, fits.HDUList):
+            pupilheader = instrument.pupil[0].header
+        else:
+            pupilfile = os.path.join(instrument._datapath, "OPD", instrument.pupil)
+            pupilheader = fits.getheader(pupilfile)
+
+        npix = pupilheader['NAXIS1']
+        pixelscale = pupilheader['PUPLSCAL'] * units.meter / units.pixel
+    return npix, pixelscale
+
+
 # Field dependent aberration class for JWST instruments
 class WebbFieldDependentAberration(poppy.OpticalElement):
     """ Field dependent aberration generated from Zernikes measured in ISIM CV testing
@@ -1932,3 +1956,143 @@ def display(self, *args, **kwargs):
             kwargs.update({'opd_vmax': 2.5e-7})
 
         return super().display(*args, **kwargs)
+
+
+class NIRCamFieldDependentWeakLens(poppy.OpticalElement):
+    """Higher-fidelity model of NIRCam weak lens(es), based on calibrated as-built performance
+    and field dependence.
+
+    Includes field-dependent variations in defocus power, and in astigmatism. Includes variation of the
+    +4 lens' effective OPD when used in a pair with either the +8 or -8 lens.
+
+    These are modeled as the specific values from the nearest neighbor ISIM CV calibration point,
+    with no interpolation between them included at this time.
+
+    See R. Telfer, 'NIRCam Weak Lens Characterization and Performance', JWST-REF-046515
+
+    Parameters
+    -----------
+    name : str
+        WLP8, WLM8, WLP4, WLM4, WLP12.
+
+    center_fp_only : bool
+        For debugging; override to set no field dependence and just use the average center field point power
+
+    include_power, include_astigmatism : bool
+        Can be used to selectively enable/disable parts of the optical model. Intended for debugging; should no
+        need to be set by users in general.
+
+    """
+
+    def __init__(self, name='WLP8', instrument=None, center_fp_only=False, verbose=False, include_power=True,
+                 include_astigmatism=True, **kwargs):
+        super().__init__(name=name)
+
+        self.ref_wavelength = 2.12e-6  # reference wavelength for defocus
+
+        self.verbose = verbose
+        if instrument is None:
+            self.module = 'A'
+            self.v2v3_coords = (0, -468 / 60)
+            npix = 1024
+        else:
+            self.module = instrument.module
+            self.v2v3_coords = instrument._tel_coords()
+            npix, pixelscale = _get_initial_pupil_sampling(instrument)
+
+        self.ztable_full = None
+
+        ## REFERENCE:
+        # NIRCam weak lenses, values from WSS config file, PRDOPSFLT-027
+        #                  A         B
+        # WLP4_diversity =   8.27309     8.3443         diversity in microns
+        # WLP8_diversity =  16.4554     16.5932
+        # WLM8_diversity = -16.4143    -16.5593
+        # WL_wavelength =    2.12                       Wavelength, in microns
+
+        if center_fp_only or instrument is None:
+            # use the center field point power only. No field dependence
+
+            # Power in P-V waves at center field point in optical model
+            # JWST-REF-046515, table 2      Mod A:    Mod B:
+            power_at_center_fp = {'WLM8': (-8.0188, -7.9521),
+                                  'WLM4': (-4.0285, -3.9766),
+                                  'WLP4': (3.9797, 3.9665),
+                                  'WLP8': (8.0292, 7.9675),
+                                  'WLP12': (12.0010, 11.9275)}
+
+            power_pv = power_at_center_fp[self.name][0 if self.module == 'A' else 1]
+            astig0 = 0
+            astig45 = 0
+
+        else:
+            closest_fp = self.find_closest_isim_fp_name(instrument)
+            if verbose: print(closest_fp)
+            power_pv, astig0, astig45 = self.lookup_empirical_lens_power(name, closest_fp)
+
+        self.power_pv_waves = power_pv
+        pv2rms_norm = self.ref_wavelength / (2 * np.sqrt(3))  # convert desired PV waves to RMS microns for power
+        # since the below function wants inputs in RMS
+
+        self.power_rms_microns = power_pv * pv2rms_norm
+
+        zernike_coefficients = np.zeros(6)
+        if include_power:
+            zernike_coefficients[3] = self.power_rms_microns
+        if include_astigmatism:
+            zernike_coefficients[4] = astig0
+            zernike_coefficients[5] = astig45
+        self.zernike_coefficients = zernike_coefficients
+
+        self.opd = poppy.zernike.opd_from_zernikes(
+            zernike_coefficients,
+            npix=npix,
+            outside=0
+        )
+        self.amplitude = np.ones_like(self.opd)
+
+    def find_closest_isim_fp_name(self, instr):
+        """Find the closest ISIM CV field point to a given instrument object,
+        i.e. the field point closest to the configured detector and coordinates
+        """
+
+        if self.ztable_full is None:
+            zernike_file = os.path.join(utils.get_webbpsf_data_path(), "si_zernikes_isim_cv3.fits")
+            self.ztable_full = Table.read(zernike_file)
+
+        lookup_name = f"NIRCam{instr.channel.upper()[0]}W{instr.module}"
+        ztable = self.ztable_full[self.ztable_full['instrument'] == lookup_name]
+
+        self._ztable = ztable
+        self._instr = instr
+        telcoords_am = instr._tel_coords().to(units.arcmin).value
+        if self.verbose: print(telcoords_am)
+        r = np.sqrt((telcoords_am[0] - ztable['V2']) ** 2 + (telcoords_am[1] - ztable['V3']) ** 2)
+        # Save closest ISIM CV3 WFE measured field point for reference
+        row = ztable[r == r.min()]
+        return row['field_point_name']
+
+    def lookup_empirical_lens_power(self, lens_name, field_point_name):
+        """ Lookup lens power and astigmatism versus field position, from empirical calibrations from ISIM CV testing
+
+        """
+        mypath = os.path.dirname(os.path.abspath(__file__)) + os.sep
+        wl_data_file = os.path.join(mypath, 'otelm', 'NIRCam_WL_Empirical_Power.csv')
+        wl_data = Table.read(wl_data_file, comment='#', header_start=1)
+
+        field_point_row = wl_data[wl_data['Field'] == field_point_name]
+        if self.verbose: print(field_point_row)
+
+        defocus_name = lens_name[2:]
+
+        power = field_point_row[defocus_name].data[0]
+        # Fringe zernike coefficients, from Telfer's table
+        z5 = field_point_row[defocus_name+"_Z5"].data[0]
+        z6 = field_point_row[defocus_name + "_Z6"].data[0]
+
+       # Have to convert Zernike normalization and order from fringe to noll, and nanometers to meters
+        astig0 = z6 / np.sqrt(6)*1e-9
+        astig45 = z5 / np.sqrt(6)*1e-9
+
+        if self.verbose: print(power)
+        return power, astig0, astig45
diff --git a/webbpsf/otelm/NIRCam_WL_Empirical_Power.csv b/webbpsf/otelm/NIRCam_WL_Empirical_Power.csv
@@ -0,0 +1,28 @@
+# NIRCam weak lens field dependent model, from Telfer 2021 JWST-REF-046515. Extracted from WeakLens_model.xlsx, tab 'wl_zerns_out'
+Module,Field,M8,M4,P4,P8,P12,M8_Z5,M8_Z6,P8_Z5,P8_Z6,P4_Z5,P4_Z6,P4P_Z5,P4P_Z6,P4M_Z5,P4M_Z6,P12_Z5,P12_Z6,M4_Z5,M4_Z6
+NIRCAMA,MIMF5,-7.79821,-3.88883,3.90222,7.79827,11.69333,11.8284,2.17329,-21.8102,-16.3917,2.66046,-25.383,0.459634,-22.7297,4.86129,-28.0362,-21.350566,-39.1214,16.68969,-25.86291
+NIRCAMA,MIMF6,-8.13141,-4.03775,4.08152,8.09957,12.16895,-4.49114,70.5055,-12.1629,-48.1705,31.0931,-68.4776,28.8704,-58.3407,33.3157,-78.6145,16.7075,-106.5112,28.82456,-8.109
+NIRCAMA,MIMF7,-7.9826,-4.02226,3.95172,7.94726,11.89036,-8.04144,-98.2402,-4.15686,59.3203,-25.0256,-3.45167,-24.089,1.04783,-25.9622,-7.95118,-28.24586,60.36813,-34.00364,-106.19138
+NIRCAMA,MIMF8,-7.62823,-3.84582,3.77703,7.60088,11.37253,19.0152,101.842,-41.7031,-98.7936,-44.6527,-25.5296,-47.6855,-27.6104,-41.6198,-23.4488,-89.3886,-126.404,-22.6046,78.3932
+NIRCAMA,MIMF9,-7.76541,-3.86337,3.89249,7.73848,11.62143,42.421,-110.782,-40.343,80.5466,16.0404,78.0426,21.1356,83.8357,10.9452,72.2496,-19.2074,164.3823,53.3662,-38.5324
+NIRCAMA,ISIM27,-7.87065,-3.9045,3.95512,7.85172,11.7958,-34.1561,-6.88395,15.6639,-3.72862,51.7335,-5.37579,50.1727,-3.17852,53.2942,-7.57307,65.8366,-6.90714,19.1381,-14.45702
+NIRCAMA,ISIM28,-7.99082,-3.99459,3.98663,7.97666,11.95369,51.563,-30.091,-46.2265,8.65496,-27.1178,-25.1817,-24.9901,-20.0175,-29.2456,-30.3459,-71.2166,-11.36254,22.3174,-60.4369
+NIRCAMA,ISIM29,-7.68073,-3.86401,3.80751,7.66073,11.45903,-41.5124,28.0907,10.0073,-42.3709,-14.2167,-15.0881,-14.8878,-10.5147,-13.5457,-19.6615,-4.8805,-52.8856,-55.0581,8.4292
+NIRCAMA,ISIM30,-7.62334,-3.80508,3.80806,7.61108,11.40895,74.4127,-26.7787,-76.8888,7.36615,-32.0242,42.1425,-28.8518,44.2886,-35.1967,39.9964,-105.7406,51.65475,39.216,13.2177
+NIRCAMA,ISIM26,-7.87941,-3.91516,3.94997,7.86507,11.80076,-17.1011,12.4959,6.7678,-7.32779,31.8001,-26.5136,33.2232,-18.088,30.377,-34.9392,39.991,-25.41579,13.2759,-22.4433
+NIRCAMA,ISIM40,-7.70295,-3.83391,3.85723,7.68943,11.53484,29.6077,-56.5997,-34.8163,37.0005,5.60501,45.6048,1.4279,52.8182,9.78212,38.3914,-33.3884,89.8187,39.38982,-18.2083
+NIRCAMA,ISIM41,-7.81231,-3.92045,3.88231,7.79565,11.6684,0.511318,-41.7969,-14.2855,13.1297,-20.8165,-9.27656,-21.5857,-9.18833,-20.0473,-9.36479,-35.8712,3.94137,-19.535982,-51.16169
+NIRCAMA,ISIM42,-7.62891,-3.82608,3.79297,7.61624,11.39935,21.2021,43.1701,-32.5087,-43.3742,-30.1805,-3.31541,-29.5299,4.43202,-30.8312,-11.0628,-62.0386,-38.94218,-9.6291,32.1073
+NIRCAMB,ISIM39,-7.7182,-3.86053,3.84659,7.70127,11.53678,14.0622,-13.8051,-15.347,1.08718,-60.8283,8.51465,-64.0182,6.43076,-57.6385,10.5985,-79.3652,7.51794,-43.5763,-3.2066
+NIRCAMB,ISIM6,-8.11717,-4.04257,4.05868,8.06502,12.10368,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+NIRCAMB,ISIM7,-7.98604,-3.98188,3.99419,7.93712,11.92135,-4.53499,76.7957,11.9644,-55.164,-77.3838,-60.3099,-82.2059,-63.666,-72.5618,-56.9537,-70.2415,-118.83,-77.09679,19.842
+NIRCAMB,ISIM8,-7.63938,-3.85989,3.7701,7.59747,11.35818,24.8685,-133.524,-35.7338,98.8711,-27.6675,32.1466,-31.4086,33.4579,-23.9264,30.8353,-67.1424,132.329,0.9421,-102.6887
+NIRCAMB,ISIM9,-7.76115,-3.92753,3.8203,7.71282,11.5198,46.0212,109.709,-47.6019,-86.7643,-33.4214,-7.23049,-31.7692,-4.75414,-35.0736,-9.70684,-79.3711,-91.51844,10.9476,100.00216
+NIRCAMB,ISIM31,-7.84441,-3.92669,3.90735,7.8213,11.71829,-21.8609,4.68125,21.7674,-16.5231,-37.4203,34.3587,-38.5679,37.0143,-36.2727,31.703,-16.8005,20.4912,-58.1336,36.38425
+NIRCAMB,ISIM32,-7.97726,-3.9661,4.00167,7.9517,11.94388,52.609,2.33611,-34.3161,-17.085,-106.923,-16.1214,-114.397,-14.8251,-99.4493,-17.4176,-148.7131,-31.9101,-46.8403,-15.08149
+NIRCAMB,ISIM33,-7.67305,-3.85003,3.81075,7.64228,11.44076,-41.4347,-46.7611,14.207,35.3524,-22.739,-9.28694,-27.9206,-9.11191,-17.5574,-9.46198,-13.7136,26.24049,-58.9921,-56.22308
+NIRCAMB,ISIM34,-7.61892,-3.84672,3.76126,7.59118,11.3415,84.8539,8.33055,-80.0563,-4.23418,-49.6867,6.72334,-50.2095,10.2347,-49.1639,3.21197,-130.2658,6.00052,35.69,11.54252
+NIRCAMB,MIMF10,-7.86285,-3.92203,3.92729,7.83413,11.74789,17.4975,-33.9691,-10.0449,8.30446,-73.9564,27.6646,-79.7687,28.1188,-68.1441,27.2104,-89.8136,36.42326,-50.6466,-6.7587
+NIRCAMB,ISIM43,-7.66466,-3.8582,3.79092,7.63245,11.40783,38.077,35.0516,-30.7943,-29.0379,-39.2717,9.30139,-37.8107,8.6219,-40.7326,9.98087,-68.605,-20.416,-2.6556,45.03247
+NIRCAMB,ISIM44,-7.79924,-3.89141,3.89664,7.76948,11.65492,8.96435,23.2658,-3.31207,-9.75147,-71.8903,-23.3948,-80.1137,-21.1682,-63.6669,-25.6214,-83.42577,-30.91967,-54.70255,-2.3556
+NIRCAMB,ISIM45,-7.6129,-3.82934,3.7704,7.5887,11.34594,30.8074,-58.2029,-32.4144,41.953,-37.7069,8.2259,-36.3188,3.54661,-39.0949,12.9052,-68.7332,45.49961,-8.2875,-45.2977
diff --git a/webbpsf/tests/test_nircam.py b/webbpsf/tests/test_nircam.py
@@ -297,9 +297,12 @@ def test_defocus(fov_arcsec=1, display=False):
     via either a weak lens, or via the options dict,
     and we get consistent results either way.
 
+    Note this is now an *inexact* comparison, because the weak lenses now include non-ideal effects, in particular field dependent astigmatism
+
     Test for #59 among other things
     """
     nrc = webbpsf_core.NIRCam()
+    nrc.set_position_from_aperture_name('NRCA3_FP1')
     nrc.pupilopd=None
     nrc.include_si_wfe=False
 
@@ -313,7 +316,7 @@ def test_defocus(fov_arcsec=1, display=False):
     nrc.options['defocus_wavelength']=2.12e-6
     psf_2 = nrc.calc_psf(nlambda=1, fov_arcsec=fov_arcsec, oversample=1, display=False, add_distortion=False)
 
-    assert np.allclose(psf[0].data, psf_2[0].data), "Defocused PSFs calculated two ways don't agree"
+    assert np.allclose(psf[0].data, psf_2[0].data, atol=1e-4), "Defocused PSFs calculated two ways don't agree as precisely as expected"
 
     if display:
         import webbpsf
@@ -332,21 +335,21 @@ def test_ways_to_specify_weak_lenses():
     testcases = (
         # FILTER  PUPIL   EXPECTED_DEFOCUS
         # Test methods directly specifying a single element
-        ('F212N', 'WLM8', 'Weak Lens -8'),
-        ('F200W', 'WLP8', 'Weak Lens +8'),
-        ('F187N', 'WLP8', 'Weak Lens +8'),
+        ('F212N', 'WLM8', 'WLM8'),
+        ('F200W', 'WLP8', 'WLP8'),
+        ('F187N', 'WLP8', 'WLP8'),
         # Note WLP4 can be specified as filter or pupil element or both
-        ('WLP4', 'WLP4', 'Weak Lens +4'),
-        (None, 'WLP4', 'Weak Lens +4'),
-        ('WLP4', None, 'Weak Lens +4'),
+        ('WLP4', 'WLP4', 'WLP4'),
+        (None, 'WLP4', 'WLP4'),
+        ('WLP4', None, 'WLP4'),
         # Test methods directly specifying a pair of elements stacked together
-        ('WLP4', 'WLM8', 'Weak Lens Pair -4'),
-        ('WLP4', 'WLP8', 'Weak Lens Pair +12'),
+        ('WLP4', 'WLM8', 'WLM4'),
+        ('WLP4', 'WLP8', 'WLP12'),
         # Test methods using virtual pupil elements WLM4 and WLP12
-        ('WLP4', 'WLM4', 'Weak Lens Pair -4'),
-        ('WLP4', 'WLP12', 'Weak Lens Pair +12'),
-        ('F212N', 'WLM4', 'Weak Lens Pair -4'),
-        ('F212N', 'WLP12', 'Weak Lens Pair +12'),
+        ('WLP4', 'WLM4', 'WLM4'),
+        ('WLP4', 'WLP12', 'WLP12'),
+        ('F212N', 'WLM4', 'WLM4'),
+        ('F212N', 'WLP12', 'WLP12'),
     )
 
     nrc = webbpsf_core.NIRCam()