Merge branch 'master' into #628

CHANGELOG.md

@@ -10,12 +10,29 @@ New Features
 - Added information about PSF size and shape to the data structure that is
   returned by EstimateShear(). (#612)
 - Added Spergel(2010) profile GSObject (#616).
+- Added an option to the ChromaticObject class that allows for image rendering
+  via interpolation between stored images.  This option can speed up the image
+  rendering process compared to brute force evaluation for chromatic objects
+  with basic properties that are wavelength-dependent. (#618)
+- Added new `ChromaticAiry` and `ChromaticOpticalPSF` classes for representing
+  optical PSFs.  These new classes allow the diffraction effects and (in the 
+  latter case) aberrations to be wavelength-dependent. (#618)
 - Enable initializing a DES_PSFEx object using a pyfits HDU directly instead
   of a filename. (#626)
+- Added TopHat class implementing a circular tophat profile. (#639)
+
 
 Bug Fixes and Improvements
 --------------------------
 
+- Changed the implementation of drawing Box and Pixel profiles in real space
+  (i.e. without being convolved by anything) to actually draw the surface 
+  brightness at the center of each pixel.  This is what all other profiles do,
+  but had not been what a Box or Pixel did. (#639)
+
 
 Updates to config options
 -------------------------
+
+- Added TopHat type. (#639)
+

examples/demo3.py

@@ -169,11 +169,17 @@ def main(argv):
     atmos = atmos.shear(e=atmos_e, beta=atmos_beta*galsim.radians)
     logger.debug('Made atmospheric PSF profile')
 
-    # Define the optical part of the PSF.
-    # The first argument of OpticalPSF below is lambda/diam, which needs to be in arcsec,
-    # so do the calculation:
+    # Define the optical part of the PSF:
+    # The first argument of OpticalPSF below is lambda/diam (wavelength of light / telescope
+    # diameter), which needs to be in the same units used to specify the image scale.  We are using
+    # arcsec for that, so we have to self-consistently use arcsec here, using the following
+    # calculation:
     lam_over_diam = lam * 1.e-9 / tel_diam # radians
     lam_over_diam *= 206265  # arcsec
+    # Note that we could also have made GalSim do the conversion for us if we did not know the right
+    # factor:
+    # lam_over_diam = lam * 1.e-9 / tel_diam * galsim.radians
+    # lam_over_diam = lam_over_diam / galsim.arcsec
     logger.debug('Calculated lambda over diam = %f arcsec', lam_over_diam)
     # The rest of the values should be given in units of the wavelength of the incident light.
     optics = galsim.OpticalPSF(lam_over_diam, 

examples/demo9.py

@@ -31,7 +31,7 @@
 
 New features introduced in this demo:
 
-- psf = OpticalPSF(..., trefoil1, trefoil2, nstruts, strut_thick, strut_angle)
+- psf = OpticalPSF(lam, diam, ..., trefoil1, trefoil2, nstruts, strut_thick, strut_angle)
 - im = galsim.ImageS(xsize, ysize, wcs)
 - pos = galsim.PositionD(x, y)
 - nfw = galsim.NFWHalo(mass, conc, z, omega_m, omega_lam)
@@ -82,7 +82,10 @@ def main(argv):
     image_size = 512       # pixels
     sky_level = 1.e2       # ADU / arcsec^2
 
-    psf_lam_over_D = 0.077 # (900nm / 2.4m) * 206265 arcsec/rad = 0.077 arcsec
+    psf_D = 2.4            # meters
+    psf_lam = 900.0        # nanometers; note that OpticalPSF will automatically convert units to
+                           # get lam/diam in units of arcsec, unless told otherwise.  In this case,
+                           # that is (900e-9m / 2.4m) * 206265 arcsec/rad = 0.077 arcsec.
     psf_obsc = 0.125       # (0.3m / 2.4m) = 0.125
     psf_nstruts = 4
     psf_strut_thick = 0.07
@@ -241,7 +244,7 @@ def yfunc1(u,v):
         # Make the PSF profile outside the loop to minimize the (significant) OpticalPSF 
         # construction overhead.
         psf = galsim.OpticalPSF(
-            lam_over_diam=psf_lam_over_D, obscuration=psf_obsc,
+            lam=psf_lam, diam=psf_D, obscuration=psf_obsc,
             nstruts=psf_nstruts, strut_thick=psf_strut_thick, strut_angle=psf_strut_angle,
             defocus=psf_defocus, astig1=psf_astig1, astig2=psf_astig2,
             coma1=psf_coma1, coma2=psf_coma2, trefoil1=psf_trefoil1, trefoil2=psf_trefoil2)

examples/demo9.yaml

@@ -31,7 +31,7 @@
 #
 # New features introduced in this demo:
 #
-# - obj type : OpticalPSF (..., trefoil1, trefoil2, nstruts, strut_thick, strut_angle)
+# - obj type : OpticalPSF (lam, diam, ..., trefoil1, trefoil2, nstruts, strut_thick, strut_angle)
 # - shear_value : NFWHaloShear (redshift)
 # - float_value : NFWHaloMagnification (redshift)
 # - input : nfw_halo (mass, conc, redshift)
@@ -55,12 +55,15 @@
 # Define the PSF profile
 psf : 
     # We use OpticalPSF with a few new features we haven't mentioned before:
-    # 1) We introduce the parameters trefoil1 and trefoil2, which are 3rd order Zernike terms.
-    # 2) We specify support struts that hold the secondary mirror.  You can specify how
+    # 1) We separately specify the light wavelength and telescope diameter, relying on OpticalPSF to
+    #    convert to appropriate units (arcsec by default, though the routine can be told otherwise).
+    # 2) We introduce the parameters trefoil1 and trefoil2, which are 3rd order Zernike terms.
+    # 3) We specify support struts that hold the secondary mirror.  You can specify how
     #    many struts to have, their thickness relative to the diameter, and the angle of the
     #    first one relative to the vertical axis.  Here, we use 4 supports.
-    type : OpticalPSF 
-    lam_over_diam : 0.077     # (900nm / 2.4m) * 206265 arcsec/rad = 0.077 arcsec
+    type : OpticalPSF
+    lam : 900.0               # nanometers
+    diam : 2.4                # meters: (900e-9m / 2.4m) * 206265 arcsec/rad = 0.077 arcsec
     obscuration : 0.125       # (0.3m / 2.4m) = 0.125
     nstruts : 4
     strut_thick : 0.07        # The default is 0.05

examples/json/demo9.json

@@ -14,7 +14,7 @@
 
 "#" : "New features introduced in this demo:",
 
-"#" : "- obj type : OpticalPSF (..., trefoil1, trefoil2, nstruts, strut_thick, strut_angle)",
+"#" : "- obj type : OpticalPSF (lam, diam, ..., trefoil1, trefoil2, nstruts, strut_thick, strut_angle)",
 "#" : "- shear_value : NFWHaloShear (redshift)",
 "#" : "- float_value : NFWHaloMagnification (redshift)",
 "#" : "- input : nfw_halo (mass, conc, redshift)",
@@ -35,13 +35,17 @@
 "#" : "Define the PSF profile",
 "psf" : {
     "#" : "We use OpticalPSF with a few new features we haven't mentioned before:",
-    "#" : "1) We introduce the parameters trefoil1 and trefoil2, which are 3rd order Zernike",
+    "#" : "1) We separately specify the light wavelength and telescope diameter, relying on"
+    "#" : "   OpticalPSF to convert to appropriate units (arcsec by default, though the routine"
+    "#" : "   can be told otherwise)."
+    "#" : "2) We introduce the parameters trefoil1 and trefoil2, which are 3rd order Zernike",
     "#" : "   terms.",
-    "#" : "2) We specify support struts that hold the secondary mirror.  You can specify how",
+    "#" : "3) We specify support struts that hold the secondary mirror.  You can specify how",
     "#" : "   many struts to have, their thickness relative to the diameter, and the angle of the",
     "#" : "   first one relative to the vertical axis.  Here, we use 4 supports.",
     "type" : "OpticalPSF",
-    "lam_over_diam" : 0.077,
+    "lam" : 900.0,
+    "diam" : 2.4,
     "obscuration" : 0.125,
     "nstruts" : 4,
     "strut_thick" : 0.07,

galsim/__init__.py

@@ -124,7 +124,7 @@
 from correlatednoise import CorrelatedNoise, getCOSMOSNoise, UncorrelatedNoise
 
 # GSObject
-from base import GSParams, GSObject, Gaussian, Moffat, Airy, Kolmogorov, Pixel, Box
+from base import GSParams, GSObject, Gaussian, Moffat, Airy, Kolmogorov, Pixel, Box, TopHat
 from base import Exponential, Sersic, DeVaucouleurs, Spergel
 from real import RealGalaxy, RealGalaxyCatalog, simReal
 from optics import OpticalPSF
@@ -141,6 +141,7 @@
 from chromatic import ChromaticObject, ChromaticAtmosphere, Chromatic, ChromaticSum
 from chromatic import ChromaticConvolution, ChromaticDeconvolution, ChromaticAutoConvolution
 from chromatic import ChromaticAutoCorrelation
+from chromatic import ChromaticOpticalPSF, ChromaticAiry
 from sed import SED
 from bandpass import Bandpass
 

galsim/base.py

@@ -1183,7 +1183,7 @@ def drawImage(self, image=None, nx=None, ny=None, bounds=None, scale=None, wcs=N
                     "This is probably an error.  Normally, you should let GalSim "
                     "handle the Pixel convolution for you.  If you want to handle the Pixel "
                     "convolution yourself, you can use method=no_pixel.  Or if you really meant "
-                    "for your profile to include the Pixel and also have GalSim convolve by"
+                    "for your profile to include the Pixel and also have GalSim convolve by "
                     "an _additional_ Pixel, you can suppress this warning by using method=fft.")
 
         # Check for scale if using nx, ny, or bounds
@@ -1260,7 +1260,7 @@ def draw(self, *args, **kwargs):
             return self.drawImage(*args, method='sb', **kwargs)
 
     def drawShoot(self, *args, **kwargs):
-        """An obsolete synonym for obj.drawImage(methos='phot')
+        """An obsolete synonym for obj.drawImage(method='phot')
         """
         normalization = kwargs.pop('normalization','f')
         if normalization in ['flux','f']:
@@ -1587,12 +1587,35 @@ class Airy(GSObject):
         >>> diam = 4.0    # meters
         >>> lam_over_diam = (lam * 1.e-9) / diam  # radians
         >>> lam_over_diam *= 206265  # Convert to arcsec
+        >>> airy = galsim.Airy(lam_over_diam)
+
+    Or, use separate keywords for the telescope diameter and wavelength in meters and nanometers,
+    respectively:
+
+        >>> airy = galsim.Airy(lam=lam, diam=diam)
+
+    in which case the user can also choose what units to use for internal descriptions of the light
+    profile using the `scale_unit` keyword (default: galsim.arcsec).  When drawing images, users
+    should then use units of `scale_unit` to specify the pixel scale.
+
 
     @param lam_over_diam    The parameter that governs the scale size of the profile.
                             See above for details about calculating it.
+    @param lam              Lambda (wavelength) in units of nanometers.  Must be supplied with
+                            `diam`, and in this case, image scales (`scale`) should be specified in
+                            units of `scale_unit`.
+    @param diam             Telescope diameter in units of meters.  Must be supplied with
+                            `lam`, and in this case, image scales (`scale`) should be specified in
+                            units of `scale_unit`.
     @param obscuration      The linear dimension of a central obscuration as a fraction of the
                             pupil dimension.  [default: 0]
     @param flux             The flux (in photons) of the profile. [default: 1]
+    @param scale_unit       Units used to define the diffraction limit and draw images, if the user
+                            has supplied a separate value for `lam` and `diam`.  Note that the
+                            results of calling methods like getFWHM() will be returned in units of
+                            `scale_unit`, as well.  Should be either a galsim.AngleUnit, or a string
+                            that can be used to construct one (e.g., 'arcsec', 'radians', etc.).
+                            [default: galsim.arcsec]
     @param gsparams         An optional GSParams argument.  See the docstring for GSParams for
                             details. [default: None]
 
@@ -1608,15 +1631,35 @@ class Airy(GSObject):
     The latter two are only available if the obscuration is 0.
     """
 
-    # Initialization parameters of the object, with type information
-    _req_params = { "lam_over_diam" : float }
-    _opt_params = { "flux" : float , "obscuration" : float }
-    _single_params = []
+    # Initialization parameters of the object, with type information.
+    # Note that this is not quite right; it's true that either lam_over_diam or lam should be
+    # supplied, but if lam is supplied then diam is required.  Errors in which parameters are used
+    # may be caught either by config or by the python code itself, depending on the particular
+    # error.
+    _req_params = { }
+    _opt_params = { "flux" : float , "obscuration" : float, "diam" : float,
+                    "scale_unit" : str }
+    _single_params = [{ "lam_over_diam" : float , "lam" : float } ]
     _takes_rng = False
     _takes_logger = False
 
     # --- Public Class methods ---
-    def __init__(self, lam_over_diam, obscuration=0., flux=1., gsparams=None):
+    def __init__(self, lam_over_diam=None, lam=None, diam=None, obscuration=0., flux=1.,
+                 scale_unit=galsim.arcsec, gsparams=None):
+        # Parse arguments: either lam_over_diam in arbitrary units, or lam in nm and diam in m.
+        # If the latter, then get lam_over_diam in units of `scale_unit`, as specified in
+        # docstring.
+        if lam_over_diam is not None:
+            if lam is not None or diam is not None:
+                raise TypeError("If specifying lam_over_diam, then do not specify lam or diam")
+        else:
+            if lam is None or diam is None:
+                raise TypeError("If not specifying lam_over_diam, then specify lam AND diam")
+            # In this case we're going to use scale_unit, so parse it in case of string input:
+            if isinstance(scale_unit, basestring):
+                scale_unit = galsim.angle.get_angle_unit(scale_unit)
+            lam_over_diam = (1.e-9*lam/diam)*(galsim.radians/scale_unit)
+
         GSObject.__init__(
             self, galsim._galsim.SBAiry(lam_over_diam=lam_over_diam, obscuration=obscuration,
                                         flux=flux, gsparams=gsparams))
@@ -1784,13 +1827,6 @@ class Pixel(GSObject):
 
         >>> scale = pixel.getScale()
 
-    Note: We have not implemented drawing a sheared or rotated Pixel in real space.  It's a
-          bit tricky to get right at the edges where fractional fluxes are required.
-          Fortunately, this is almost never needed.  Pixels are almost always convolved by
-          something else rather than drawn by themselves, in which case either the fourier
-          space method is used, or photon shooting.  Both of these are implemented in GalSim.
-          If you need to draw sheared or rotated Pixels in real space, please file an issue, and
-          maybe we'll implement that function.  Until then, you will get an exception if you try.
     """
 
     # Initialization parameters of the object, with type information
@@ -1831,13 +1867,6 @@ class Box(GSObject):
         >>> width = box.getWidth()
         >>> height = box.getHeight()
 
-    Note: We have not implemented drawing a sheared or rotated Box in real space.  It's a
-          bit tricky to get right at the edges where fractional fluxes are required.
-          Fortunately, this is almost never needed.  Box profiles are almost always convolved
-          by something else rather than drawn by themselves, in which case either the fourier
-          space method is used, or photon shooting.  Both of these are implemented in GalSim.
-          If you need to draw sheared or rotated Boxes in real space, please file an issue, and
-          maybe we'll implement that function.  Until then, you will get an exception if you try.
     """
 
     # Initialization parameters of the object, with type information
@@ -1864,6 +1893,45 @@ def getHeight(self):
         return self.SBProfile.getHeight()
 
 
+class TopHat(GSObject):
+    """A class describing a radial tophat profile.  This profile is a constant value within some
+    radius, and zero outside this radius.
+
+    Initialization
+    --------------
+
+    @param radius           The radius of the TopHat, where the surface brightness drops to 0.
+    @param flux             The flux (in photons) of the profile. [default: 1]
+    @param gsparams         An optional GSParams argument.  See the docstring for GSParams for
+                            details. [default: None]
+
+    Methods
+    -------
+
+    In addition to the usual GSObject methods, TopHat has the following access method:
+
+        >>> radius = tophat.getRadius()
+
+    """
+
+    # Initialization parameters of the object, with type information
+    _req_params = { "radius" : float }
+    _opt_params = { "flux" : float }
+    _single_params = []
+    _takes_rng = False
+    _takes_logger = False
+
+    # --- Public Class methods ---
+    def __init__(self, radius, flux=1., gsparams=None):
+        radius = float(radius)
+        GSObject.__init__(self, galsim._galsim.SBTopHat(radius, flux=flux, gsparams=gsparams))
+
+    def getRadius(self):
+        """Return the radius of the tophat profile.
+        """
+        return self.SBProfile.getRadius()
+
+
 class Sersic(GSObject):
     """A class describing a Sersic profile.