Add lam/diam/r0 Quantity support to gsobjects

galsim/airy.py

@@ -19,6 +19,7 @@
 __all__ = [ 'Airy' ]
 
 import math
+import astropy.units as u
 
 from . import _galsim
 from .gsobject import GSObject
@@ -80,13 +81,16 @@ class Airy(GSObject):
         gsparams:       An optional `GSParams` argument. [default: None]
     """
     _req_params = { }
-    _opt_params = { "flux" : float , "obscuration" : float, "diam" : float,
-                    "scale_unit" : str }
+    _opt_params = { "flux" : float ,
+                    "obscuration" : float,
+                    "diam" : (float, u.Quantity),
+                    "scale_unit" : str
+                  }
     # 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.
-    _single_params = [{ "lam_over_diam" : float , "lam" : float } ]
+    _single_params = [{ "lam_over_diam" : float , "lam" : (float, u.Quantity) } ]
 
     # For an unobscured Airy, we have the following factor which can be derived using the
     # integral result given in the Wikipedia page (http://en.wikipedia.org/wiki/Airy_disk),
@@ -108,6 +112,11 @@ def __init__(self, lam_over_diam=None, lam=None, diam=None, obscuration=0., flux
         self._flux = float(flux)
         self._gsparams = GSParams.check(gsparams)
 
+        if isinstance(lam, u.Quantity):
+            lam = lam.to_value(u.nm)
+        if isinstance(diam, u.Quantity):
+            diam = diam.to_value(u.m)
+
         # 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.

galsim/kolmogorov.py

@@ -19,6 +19,7 @@
 __all__ = [ 'Kolmogorov' ]
 
 import numpy as np
+import astropy.units as u
 import math
 
 from . import _galsim
@@ -120,12 +121,21 @@ class Kolmogorov(GSObject):
         fwhm:               The full-width half-max size
         half_light_radius:  The half-light radius
     """
-    _opt_params = { "flux" : float, "r0" : float, "r0_500" : float, "scale_unit" : str }
+    _opt_params = {
+        "flux" : float,
+        "r0" : (float, u.Quantity),
+        "r0_500" : (float, u.Quantity),
+        "scale_unit" : str
+    }
     # Note that this is not quite right; it's true that exactly one of these 4 should be supplied,
     # but if lam is supplied then r0 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.
-    _single_params = [ { "lam_over_r0" : float, "fwhm" : float, "half_light_radius" : float,
-                         "lam" : float } ]
+    _single_params = [ {
+        "lam_over_r0" : float,
+        "fwhm" : float,
+        "half_light_radius" : float,
+        "lam" : (float, u.Quantity)
+    } ]
 
     # The FWHM of the Kolmogorov PSF is ~0.976 lambda/r0 (e.g., Racine 1996, PASP 699, 108).
     # In SBKolmogorov.cpp we refine this factor to 0.9758634299
@@ -166,6 +176,13 @@ def __init__(self, lam_over_r0=None, fwhm=None, half_light_radius=None, lam=None
         self._flux = float(flux)
         self._gsparams = GSParams.check(gsparams)
 
+        if isinstance(lam, u.Quantity):
+            lam = lam.to_value(u.nm)
+        if isinstance(r0, u.Quantity):
+            r0 = r0.to_value(u.m)
+        if isinstance(r0_500, u.Quantity):
+            r0_500 = r0_500.to_value(u.m)
+
         if fwhm is not None :
             if any(item is not None for item in (lam_over_r0, lam, r0, r0_500, half_light_radius)):
                 raise GalSimIncompatibleValuesError(

galsim/phase_psf.py

@@ -20,6 +20,7 @@
 
 from heapq import heappush, heappop
 import numpy as np
+import astropy.units as u
 import copy
 
 from . import fits
@@ -1807,7 +1808,7 @@ class OpticalPSF(GSObject):
         gsparams:           An optional `GSParams` argument. [default: None]
     """
     _opt_params = {
-        "diam": float,
+        "diam": (float, u.Quantity),
         "defocus": float,
         "astig1": float,
         "astig2": float,
@@ -1833,7 +1834,7 @@ class OpticalPSF(GSObject):
         "pupil_plane_size": float,
         "scale_unit": str,
         "fft_sign": str}
-    _single_params = [{"lam_over_diam": float, "lam": float}]
+    _single_params = [{"lam_over_diam": float, "lam": (float, u.Quantity)}]
 
     _has_hard_edges = False
     _is_axisymmetric = False
@@ -1852,6 +1853,10 @@ def __init__(self, lam_over_diam=None, lam=None, diam=None, tip=0., tilt=0., def
                  pupil_angle=0.*radians, scale_unit=arcsec, fft_sign='+', gsparams=None,
                  _force_stepk=0., _force_maxk=0.,
                  suppress_warning=False, geometric_shooting=False):
+        if isinstance(lam, u.Quantity):
+            lam = lam.to_value(u.nm)
+        if isinstance(diam, u.Quantity):
+            diam = diam.to_value(u.m)
         if fft_sign not in ['+', '-']:
             raise GalSimValueError("Invalid fft_sign", fft_sign, allowed_values=['+','-'])
         if isinstance(scale_unit, str):

galsim/second_kick.py

@@ -18,6 +18,8 @@
 
 __all__ = [ 'SecondKick' ]
 
+import astropy.units as u
+
 from . import _galsim
 from .gsobject import GSObject
 from .gsparams import GSParams
@@ -84,7 +86,11 @@ class SecondKick(GSObject):
                         construct one (e.g., 'arcsec', 'radians', etc.). [default: galsim.arcsec]
         gsparams:       An optional `GSParams` argument. [default: None]
     """
-    _req_params = { "lam" : float, "r0" : float, "diam" : float }
+    _req_params = {
+        "lam" : (float, u.Quantity),
+        "r0" : (float, u.Quantity),
+        "diam" : (float, u.Quantity),
+    }
     _opt_params = { "obscuration" : float, "kcrit" : float, "flux" : float, "scale_unit" : str }
 
     _has_hard_edges = False
@@ -94,6 +100,13 @@ class SecondKick(GSObject):
 
     def __init__(self, lam, r0, diam, obscuration=0, kcrit=0.2, flux=1,
                  scale_unit=arcsec, gsparams=None):
+        if isinstance(lam, u.Quantity):
+            lam = lam.to_value(u.nm)
+        if isinstance(r0, u.Quantity):
+            r0 = r0.to_value(u.m)
+        if isinstance(diam, u.Quantity):
+            diam = diam.to_value(u.m)
+
         if isinstance(scale_unit, str):
             self._scale_unit = AngleUnit.from_name(scale_unit)
         else:

galsim/vonkarman.py

@@ -19,6 +19,7 @@
 __all__ = [ 'VonKarman' ]
 
 import numpy as np
+import astropy.units as u
 
 from . import _galsim
 from .gsobject import GSObject
@@ -102,9 +103,17 @@ class VonKarman(GSObject):
                             keyword.  [default: False]
         gsparams:           An optional `GSParams` argument. [default: None]
     """
-    _req_params = { "lam" : float }
-    _opt_params = { "L0" : float, "flux" : float, "scale_unit" : str, "do_delta" : bool }
-    _single_params = [ { "r0" : float, "r0_500" : float } ]
+    _req_params = { "lam" : (float, u.Quantity) }
+    _opt_params = {
+        "L0" : (float, u.Quantity),
+        "flux" : float,
+        "scale_unit" : str,
+        "do_delta" : bool
+    }
+    _single_params = [ {
+        "r0" : (float, u.Quantity),
+        "r0_500" : (float, u.Quantity)
+    } ]
 
     _has_hard_edges = False
     _is_axisymmetric = True
@@ -113,6 +122,15 @@ class VonKarman(GSObject):
 
     def __init__(self, lam, r0=None, r0_500=None, L0=25.0, flux=1, scale_unit=arcsec,
                  force_stepk=0.0, do_delta=False, suppress_warning=False, gsparams=None):
+        if isinstance(lam, u.Quantity):
+            lam = lam.to_value(u.nm)
+        if isinstance(r0, u.Quantity):
+            r0 = r0.to_value(u.m)
+        if isinstance(r0_500, u.Quantity):
+            r0_500 = r0_500.to_value(u.m)
+        if isinstance(L0, u.Quantity):
+            L0 = L0.to_value(u.m)
+
         # We lose stability if L0 gets too large.  This should be close enough to infinity for
         # all practical purposes though.
         if L0 > 1e10:

tests/test_config_gsobject.py

@@ -18,6 +18,7 @@
 
 import logging
 import numpy as np
+import astropy.units as u
 import os
 import sys
 
@@ -222,6 +223,7 @@ def test_airy():
                     'gsparams' : { 'xvalue_accuracy' : 1.e-2 }
                  },
         'gal6' : { 'type' : 'Airy' , 'lam' : 400., 'diam' : 4.0, 'scale_unit' : 'arcmin' },
+        'gal7' : { 'type' : 'Airy' , 'lam' : 400*u.nm, 'diam' : '4000 mm', 'scale_unit' : 'arcmin' },
         'bad1' : { 'type' : 'Airy' , 'lam_over_diam' : 0.4, 'lam' : 400, 'diam' : 10 },
         'bad2' : { 'type' : 'Airy' , 'flux' : 1.3 },
         'bad3' : { 'type' : 'Airy' , 'lam_over_diam' : 0.4, 'obsc' : 0.3, 'flux' : 100 },
@@ -259,6 +261,8 @@ def test_airy():
     gal6a = galsim.config.BuildGSObject(config, 'gal6')[0]
     gal6b = galsim.Airy(lam=400., diam=4., scale_unit=galsim.arcmin)
     gsobject_compare(gal6a, gal6b)
+    gal7a = galsim.config.BuildGSObject(config, 'gal7')[0]
+    gsobject_compare(gal6a, gal7a)
 
     # Make sure they don't match when using the default GSParams
     gal5c = galsim.Airy(lam_over_diam=45)
@@ -298,6 +302,8 @@ def test_kolmogorov():
         'gal5' : { 'type' : 'Kolmogorov' , 'lam_over_r0' : 1, 'flux' : 50,
                    'gsparams' : { 'integration_relerr' : 1.e-2, 'integration_abserr' : 1.e-4 }
                  },
+        'gal6' : { 'type' : 'Kolmogorov' , 'lam' : '400 nm', 'r0_500' : '15 cm' },
+        'gal7' : { 'type' : 'Kolmogorov' , 'lam' : '$4000*u.Angstrom', 'r0' : 0.18 },
         'bad1' : { 'type' : 'Kolmogorov' , 'fwhm' : 2, 'lam_over_r0' : 3, 'flux' : 100 },
         'bad2' : { 'type' : 'Kolmogorov', 'flux' : 100 },
         'bad3' : { 'type' : 'Kolmogorov' , 'lam_over_r0' : 2, 'lam' : 400, 'r0' : 0.15 },
@@ -334,13 +340,103 @@ def test_kolmogorov():
     with assert_raises(AssertionError):
         gsobject_compare(gal5a, gal5c)
 
+    gal6a = galsim.config.BuildGSObject(config, 'gal6')[0]
+    gal6b = galsim.Kolmogorov(lam=400*u.nm, r0_500=15*u.cm)
+    gsobject_compare(gal6a, gal6b)
+
+    gal7a = galsim.config.BuildGSObject(config, 'gal7')[0]
+    gal7b = galsim.Kolmogorov(lam=4000*u.Angstrom, r0=0.18)
+    gsobject_compare(gal7a, gal7b)
+
     with assert_raises(galsim.GalSimConfigError):
         galsim.config.BuildGSObject(config, 'bad1')
     with assert_raises(galsim.GalSimConfigError):
         galsim.config.BuildGSObject(config, 'bad2')
     with assert_raises(galsim.GalSimConfigError):
         galsim.config.BuildGSObject(config, 'bad3')
 
+
+@timer
+def test_VonKarman():
+    """Test various ways to build a VonKarman
+    """
+    config = {
+        'gal1' : { 'type' : 'VonKarman' , 'lam' : 500, 'r0' : 0.2 },
+        'gal2' : { 'type' : 'VonKarman' , 'lam' : 760, 'r0_500' : 0.2 },
+        'gal3' : { 'type' : 'VonKarman' , 'lam' : 450*u.nm, 'r0_500' : '20 cm', 'flux' : 1.e6,
+                   'ellip' : { 'type' : 'QBeta' , 'q' : 0.6, 'beta' : 0.39 * galsim.radians }
+                 },
+        'gal4' : { 'type' : 'VonKarman' , 'lam' : 500, 'r0' : 0.2,
+                   'dilate' : 3, 'ellip' : galsim.Shear(e1=0.3),
+                   'rotate' : 12 * galsim.degrees,
+                   'lens' : {
+                       'shear' : galsim.Shear(g1=0.03, g2=-0.05),
+                       'mu' : 1.03,
+                   },
+                   'shift' : { 'type' : 'XY', 'x' : 0.7, 'y' : -1.2 },
+                 },
+        'gal5' : { 'type' : 'VonKarman' , 'lam' : 500, 'r0' : 0.2, 'do_delta' : True,
+                   'gsparams' : { 'integration_relerr' : 1.e-2, 'integration_abserr' : 1.e-4 }
+                 },
+        'bad1' : { 'type' : 'VonKarman' , 'fwhm' : 2, 'lam_over_r0' : 3, 'flux' : 100 },
+        'bad2' : { 'type' : 'VonKarman', 'flux' : 100 },
+        'bad3' : { 'type' : 'VonKarman' , 'lam' : 400, 'r0' : 0.15, 'r0_500' : 0.12 },
+    }
+
+    gal1a = galsim.config.BuildGSObject(config, 'gal1')[0]
+    gal1b = galsim.VonKarman(lam = 500, r0 = 0.2)
+    gsobject_compare(gal1a, gal1b)
+
+    gal2a = galsim.config.BuildGSObject(config, 'gal2')[0]
+    gal2b = galsim.VonKarman(lam = 760, r0_500 = 0.2)
+    gsobject_compare(gal2a, gal2b)
+
+    gal3a = galsim.config.BuildGSObject(config, 'gal3')[0]
+    gal3b = galsim.VonKarman(lam = 450*u.nm, r0_500 = 20*u.cm, flux = 1.e6)
+    gal3b = gal3b.shear(q = 0.6, beta = 0.39 * galsim.radians)
+    gsobject_compare(gal3a, gal3b)
+
+
+    gal4a = galsim.config.BuildGSObject(config, 'gal4')[0]
+    gal4b = galsim.VonKarman(lam = 500, r0 = 0.2)
+    gal4b = gal4b.dilate(3).shear(e1 = 0.3).rotate(12 * galsim.degrees)
+    gal4b = gal4b.lens(0.03, -0.05, 1.03).shift(dx = 0.7, dy = -1.2)
+    gsobject_compare(gal4a, gal4b)
+
+    gal5a = galsim.config.BuildGSObject(config, 'gal5')[0]
+    gsparams = galsim.GSParams(integration_relerr=1.e-2, integration_abserr=1.e-4)
+    gal5b = galsim.VonKarman(lam=500, r0=0.2, do_delta=True, gsparams=gsparams)
+    gsobject_compare(gal5a, gal5b)
+
+    with assert_raises(galsim.GalSimConfigError):
+        galsim.config.BuildGSObject(config, 'bad1')
+    with assert_raises(galsim.GalSimConfigError):
+        galsim.config.BuildGSObject(config, 'bad2')
+    with assert_raises(galsim.GalSimConfigError):
+        galsim.config.BuildGSObject(config, 'bad3')
+
+
+@timer
+def test_secondKick():
+    config = {
+        'gal1' : { 'type' : 'SecondKick' , 'lam' : 500, 'r0' : 0.2, 'diam': 4.0 },
+        'gal2' : { 'type' : 'SecondKick' , 'lam' : '0.5 micron', 'r0' : '10 cm', 'diam': 8.2*u.m, 'flux' : 100 },
+        'gal3' : { 'type' : 'SecondKick' , 'lam' : '$2*450*u.nm', 'r0' : 0.2, 'diam': 4.0, 'obscuration' : 0.2, 'kcrit' : 0.1 },
+    }
+
+    gal1a = galsim.config.BuildGSObject(config, 'gal1')[0]
+    gal1b = galsim.SecondKick(lam = 500, r0 = 0.2, diam=4.0)
+    gsobject_compare(gal1a, gal1b)
+
+    gal2a = galsim.config.BuildGSObject(config, 'gal2')[0]
+    gal2b = galsim.SecondKick(lam = 0.5*u.micron, r0 = 10*u.cm, diam=8.2*u.m, flux = 100)
+    gsobject_compare(gal2a, gal2b)
+
+    gal3a = galsim.config.BuildGSObject(config, 'gal3')[0]
+    gal3b = galsim.SecondKick(lam = 2*450, r0 = 0.2, diam=4.0, obscuration=0.2, kcrit=0.1)
+    gsobject_compare(gal3a, gal3b)
+
+
 @timer
 def test_opticalpsf():
     """Test various ways to build a OpticalPSF
@@ -371,10 +467,10 @@ def test_opticalpsf():
                    'pupil_plane_im' :
                        os.path.join(".","Optics_comparison_images","sample_pupil_rolled.fits"),
                    'pupil_angle' : 27.*galsim.degrees },
-        'gal6' : {'type' : 'OpticalPSF' , 'lam' : 874.0, 'diam' : 7.4, 'flux' : 70.,
+        'gal6' : {'type' : 'OpticalPSF' , 'lam' : '874 nm', 'diam' : '7.4 m', 'flux' : 70.,
                   'aberrations' : [0.06, 0.12, -0.08, 0.07, 0.04, 0.0, 0.0, -0.13],
                   'obscuration' : 0.1 },
-        'gal7' : {'type' : 'OpticalPSF' , 'lam' : 874.0, 'diam' : 7.4, 'aberrations' : []},
+        'gal7' : {'type' : 'OpticalPSF' , 'lam' : 0.874*u.micron, 'diam' : '$740*u.cm', 'aberrations' : []},
         'bad1' : {'type' : 'OpticalPSF' , 'lam' : 874.0, 'diam' : 7.4, 'lam_over_diam' : 0.2},
         'bad2' : {'type' : 'OpticalPSF' , 'lam_over_diam' : 0.2,
                   'aberrations' : "0.06, 0.12, -0.08, 0.07, 0.04, 0.0, 0.0, -0.13"},