mperrin / webbpsf

James Webb Space Telescope PSF simulation tool - NOTE THIS VERSION OF REPO IS SUPERCEDED BY spacetelescope/webbpsf
BSD 3-Clause "New" or "Revised" License
16 stars 15 forks source link

NIRISS GR700XD filter PSF calculation fails on Quantity comparison #148

Closed josePhoenix closed 7 years ago

josePhoenix commented 7 years ago

Kevin Volk reports being unable to calculate a GR700XD PSF with WebbPSF v0.5.1. Other pupil masks work. Seems to be the result of using Astropy units internally without some necessary conversion.

Minimal example:

import webbpsf
niriss = webbpsf.NIRISS()
niriss.pupil_mask = 'GR700XD'
niriss.calc_psf(monochromatic=1e-6, fov_pixels=2)

Output:

In [1]: %paste
import webbpsf
niriss = webbpsf.NIRISS()
niriss.pupil_mask = 'GR700XD'
niriss.calc_psf(monochromatic=1e-6, fov_pixels=2)

## -- End pasted text --
/Users/jlong/software/webbpsf/webbpsf/optics.py:341: VisibleDeprecationWarning: boolean index did not match indexed array along dimension 0; dimension is 1024 but corresponding boolean dimension is 1
  sag[self.amplitude == 0] = 0 # no OPD in opaque regions (makes no difference in propagation but improves display)
---------------------------------------------------------------------------
UnitsError                                Traceback (most recent call last)
<ipython-input-1-1ad50033a008> in <module>()
      2 niriss = webbpsf.NIRISS()
      3 niriss.pupil_mask = 'GR700XD'
----> 4 niriss.calc_psf(monochromatic=1e-6, fov_pixels=2)

/Users/jlong/software/poppy/poppy/instrument.py in calc_psf(self, outfile, source, nlambda, monochromatic, fov_arcsec, fov_pixels, oversample, detector_oversample, fft_oversample, rebin, clobber, display, save_intermediates, return_intermediates, normalize)
    263         self.optsys = self._getOpticalSystem(fov_arcsec=fov_arcsec, fov_pixels=fov_pixels,
    264                                              fft_oversample=fft_oversample, detector_oversample=detector_oversample,
--> 265                                              options=local_options)
    266         self._check_for_aliasing(wavelens)
    267         # and use it to compute the PSF (the real work happens here, in code in poppy.py)

/Users/jlong/software/webbpsf/webbpsf/webbpsf_core.py in _getOpticalSystem(self, fft_oversample, detector_oversample, fov_arcsec, fov_pixels, options)
    591         optsys = SpaceTelescopeInstrument._getOpticalSystem(self,
    592             fft_oversample=fft_oversample, detector_oversample=detector_oversample, fov_arcsec=fov_arcsec, fov_pixels=fov_pixels,
--> 593             options=options)
    594         optsys.planes[0].display_annotate = utils.annotate_ote_entrance_coords
    595         #optsys.planes[-2].display_annotate = utils.annotate_sky_pupil_coords

/Users/jlong/software/webbpsf/webbpsf/webbpsf_core.py in _getOpticalSystem(self, fft_oversample, detector_oversample, fov_arcsec, fov_pixels, options)
    444         if self.image_mask is not None or self.pupil_mask is not None or ('force_coron' in options and options['force_coron']):
    445             _log.debug("Adding coronagraph/spectrograph optics...")
--> 446             optsys, trySAM, SAM_box_size = self._addAdditionalOptics(optsys, oversample=fft_oversample)
    447         else: trySAM = False
    448

/Users/jlong/software/webbpsf/webbpsf/webbpsf_core.py in _addAdditionalOptics(self, optsys, oversample)
   1273             optsys.planes[-1].wavefront_display_hint='intensity'
   1274         elif self.pupil_mask == 'GR700XD':
-> 1275             optsys.add_pupil(optic = NIRISS_GR700XD_Grism(shift=shift))
   1276
   1277         elif (self.pupil_mask  is None and self.image_mask is not None):

/Users/jlong/software/webbpsf/webbpsf/optics.py in __init__(self, name, which, shift)
    266
    267         # perform an initial population of the OPD array for display etc.
--> 268         tmp = self.get_phasor(poppy.Wavefront(2e-6))
    269
    270     def get_opd(self, wave):

/Users/jlong/software/poppy/poppy/optics.py in get_phasor(self, wave)
    106         scale = 2. * np.pi / wavelength.to(u.meter).value
    107
--> 108         return self.get_transmission(wave) * np.exp (1.j * self.get_opd(wave) * scale)
    109
    110     def getPhasor(self,wave):

/Users/jlong/software/webbpsf/webbpsf/optics.py in get_opd(self, wave)
    345
    346         # scale for index of refraction
--> 347         index = self.ZnS_index(wavelength)
    348         opd = sag *  ( index-1)
    349         _log.debug(" Scaling for ZnS index of refraction {0} at {1:.3g} microns".format(index, wavelength*1e6))

/Users/jlong/software/webbpsf/webbpsf/optics.py in ZnS_index(self, wavelength, temperature)
    403             S_i =       S_ij[i,0]      + S_ij[i,1]     *T + S_ij[i,2]     *T**2.0 + S_ij[i,3]     *T**3.0 + S_ij[i,4]     *T**4.0
    404             lambda_i =  lambda_ij[i,0] + lambda_ij[i,1]*T + lambda_ij[i,2]*T**2.0 + lambda_ij[i,3]*T**3.0 + lambda_ij[i,4]*T**4.0
--> 405             n2minus1 += S_i*lambda_micron**2.0/(lambda_micron**2.0 - lambda_i**2.0)
    406
    407         cleartran_index = np.sqrt(1.0 + n2minus1)

/Users/jlong/software/astropy/astropy/units/quantity.py in __array_prepare__(self, obj, context)
    376                                      "argument is not a quantity (unless the "
    377                                      "latter is all zero/infinity/nan)"
--> 378                                      .format(function.__name__))
    379             except TypeError:
    380                 # _can_have_arbitrary_unit failed: arg could not be compared

UnitsError: Can only apply 'subtract' function to dimensionless quantities when other argument is not a quantity (unless the latter is all zero/infinity/nan)