From 1cfd49919df93a0adc8cc8c45b5dd5bf98b26d7b Mon Sep 17 00:00:00 2001 From: Josh Meyers Date: Mon, 16 Sep 2024 11:35:21 -0700 Subject: [PATCH] Add lam/diam/r0 Quantity support to gsobjects --- galsim/airy.py | 15 ++++- galsim/kolmogorov.py | 23 +++++++- galsim/phase_psf.py | 9 ++- galsim/second_kick.py | 15 ++++- galsim/vonkarman.py | 24 +++++++- tests/test_config_gsobject.py | 100 +++++++++++++++++++++++++++++++++- 6 files changed, 172 insertions(+), 14 deletions(-) diff --git a/galsim/airy.py b/galsim/airy.py index ff61d3e150..301b406edf 100644 --- a/galsim/airy.py +++ b/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. diff --git a/galsim/kolmogorov.py b/galsim/kolmogorov.py index 5b7813cbfc..615a542a64 100644 --- a/galsim/kolmogorov.py +++ b/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( diff --git a/galsim/phase_psf.py b/galsim/phase_psf.py index 1284659061..c459f440ff 100644 --- a/galsim/phase_psf.py +++ b/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): diff --git a/galsim/second_kick.py b/galsim/second_kick.py index 1bb5372c70..055ac1e428 100644 --- a/galsim/second_kick.py +++ b/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: diff --git a/galsim/vonkarman.py b/galsim/vonkarman.py index 1eb2ec6b6f..33e41eb8d2 100644 --- a/galsim/vonkarman.py +++ b/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: diff --git a/tests/test_config_gsobject.py b/tests/test_config_gsobject.py index e76184a00d..68627e80c6 100644 --- a/tests/test_config_gsobject.py +++ b/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,6 +340,14 @@ 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): @@ -341,6 +355,88 @@ def test_kolmogorov(): 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"},