# Copyright (c) 2012-2023 by the GalSim developers team on GitHub
# https://github.com/GalSim-developers
#
# This file is part of GalSim: The modular galaxy image simulation toolkit.
# https://github.com/GalSim-developers/GalSim
#
# GalSim is free software: redistribution and use in source and binary forms,
# with or without modification, are permitted provided that the following
# conditions are met:
#
# 1. Redistributions of source code must retain the above copyright notice, this
# list of conditions, and the disclaimer given in the accompanying LICENSE
# file.
# 2. Redistributions in binary form must reproduce the above copyright notice,
# this list of conditions, and the disclaimer given in the documentation
# and/or other materials provided with the distribution.
#
__all__ = [ 'ChromaticObject', 'ChromaticAtmosphere', 'ChromaticSum',
'ChromaticConvolution', 'ChromaticDeconvolution',
'ChromaticAutoConvolution', 'ChromaticAutoCorrelation',
'ChromaticTransformation', 'SimpleChromaticTransformation',
'ChromaticFourierSqrtProfile', 'ChromaticOpticalPSF',
'ChromaticAiry', 'InterpolatedChromaticObject', ]
import math
import numpy as np
from astropy import units
import copy
from .gsobject import GSObject
from .sed import SED
from .bandpass import Bandpass
from .position import Position, _PositionD
from ._utilities import lazy_property, doc_inherit
from .gsparams import GSParams
from .phase_psf import OpticalPSF, Aperture
from .table import _LookupTable
from .sum import Add
from .errors import GalSimError, GalSimRangeError, GalSimSEDError, GalSimValueError
from .errors import GalSimIncompatibleValuesError, GalSimNotImplementedError, galsim_warn
from .photon_array import WavelengthSampler, PhotonArray, PhotonOp, ScaleFlux
from .random import BaseDeviate, UniformDeviate, BinomialDeviate, PoissonDeviate
from .shear import Shear, _Shear
from .interpolatedimage import InterpolatedImage
from .angle import Angle, _Angle, AngleUnit, arcsec, radians
from .airy import Airy
from .deltafunction import DeltaFunction
from . import utilities
from . import integ
from . import dcr
[docs]class ChromaticObject:
"""Base class for defining wavelength-dependent objects.
This class primarily serves as the base class for chromatic subclasses. See the docstrings for
subclasses for more details.
A ChromaticObject can be instantiated directly from an existing `GSObject`. In this case, the
newly created ChromaticObject will act in nearly the same way as the original `GSObject` works,
except that it has access to the ChromaticObject transformation methods described below (e.g.,
expand(), dilate(), shift(), withFlux(), ...) These can all take functions as arguments to
describe wavelength-dependent transformations. E.g.,::
>>> gsobj = galsim.Gaussian(fwhm=1)
>>> chrom_obj = galsim.ChromaticObject(gsobj).dilate(lambda wave: (wave/500.)**(-0.2))
In this and similar cases, the argument to the transformation method should be a python callable
that accepts wavelength in nanometers and returns whatever type the transformation method
normally accepts (so an int or float above).
One caveat to creating a ChromaticObject directly from a `GSObject` like this is that even
though the source `GSObject` instance has flux units in photons/s/cm^2, the newly formed
ChromaticObject will be interpreted as dimensionless, i.e., it will have a dimensionless `SED`
(and have its .dimensionless attribute set to True). See below for more discussion about the
dimensions of ChromaticObjects.
Another way to instantiate a ChromaticObject from a `GSObject` is to multiply by an `SED`.
This can be useful to consistently generate the same galaxy observed through different filters,
or, with `ChromaticSum`, to construct multi-component galaxies, each component with a different
`SED`. For example, a bulge+disk galaxy could be constructed::
>>> bulge_sed = user_function_to_get_bulge_spectrum()
>>> disk_sed = user_function_to_get_disk_spectrum()
>>> bulge_mono = galsim.DeVaucouleurs(half_light_radius=1.0)
>>> disk_mono = galsim.Exponential(half_light_radius=2.0)
>>> bulge = bulge_mono * bulge_sed
>>> disk = disk_mono * disk_sed
>>> gal = bulge + disk
The ``sed`` instances above describe the flux density in photons/nm/cm^2/s of an object, possibly
normalized with either the `SED.withFlux` or `SED.withMagnitude` methods (see their docstrings
for details about these and other normalization options). Note that for dimensional
consistency, in this case, the ``flux`` attribute of the multiplied `GSObject` is interpreted
as being dimensionless instead of in its normal units of [photons/s/cm^2]. The photons/s/cm^2
units are (optionally) carried by the `SED` instead, or even left out entirely if the `SED` is
dimensionless itself (see discussion on ChromaticObject dimensions below). The `GSObject`
``flux`` attribute *does* still contribute to the ChromaticObject normalization, though.
For example, the following are equivalent::
>>> chrom_obj = (sed * 3.0) * gsobj
>>> chrom_obj2 = sed * (gsobj * 3.0)
Subclasses that instantiate a ChromaticObject directly, such as `ChromaticAtmosphere`, also
exist. Even in this case, however, the underlying implementation always eventually wraps one
or more `GSObject` instances.
**Dimensions**:
ChromaticObjects can generally be sorted into two distinct types: those that represent galaxies
or stars and have dimensions of [photons/wavelength-interval/area/time/solid-angle], and those
that represent other types of wavelength dependence besides flux, like chromatic PSFs (these
have dimensions of [1/solid-angle]). The former category of ChromaticObjects will have their
``.spectral`` attribute set to True, while the latter category of ChromaticObjects will have
their ``.dimensionless`` attribute set to True. These two classes of ChromaticObjects have
different restrictions associated with them. For example, only spectral ChromaticObjects can
be drawn using `drawImage`, only ChromaticObjects of the same type can be added together, and
at most one spectral ChromaticObject can be part of a `ChromaticConvolution`.
Multiplying a dimensionless ChromaticObject a spectral `SED` produces a spectral ChromaticObject
(though note that the new object's `SED` may not be equal to the SED being multiplied by since
the original ChromaticObject may not have had unit normalization.)
**Methods**:
`evaluateAtWavelength` returns the monochromatic surface brightness profile (as a `GSObject`)
at a given wavelength (in nanometers).
`interpolate` can be used for non-separable ChromaticObjects to expedite the image rendering
process. See the docstring of that method for more details and discussion of when this is a
useful tool (and the interplay between interpolation, object transformations, and convolutions).
Also, ChromaticObject has most of the same methods as `GSObject` with the following exceptions:
The `GSObject` access methods (e.g. `GSObject.xValue`, `GSObject.maxk`, etc.) are not available.
Instead, you would need to evaluate the profile at a particular wavelength and access what you
want from that.
The `withFlux`, `withFluxDensity`, and `withMagnitude` methods will return a new chromatic
object with the appropriate spatially integrated flux, flux density, or magnitude.
The `drawImage` method draws the object as observed through a particular bandpass, so the
arguments are somewhat different. See the docstring of `drawImage` for more details.
"""
# ChromaticObjects should adhere to the following invariants:
# - Objects should define the attributes/properties:
# * .sed, .separable, .wave_list, .interpolated, .deinterpolated, .spectral, .dimensionless
# - obj.evaluateAtWavelength(lam).drawImage().array.sum() == obj.sed(lam)
# == obj.evaluateAtWavelength(lam).flux
# - if obj.spectral:
# obj.sed.calculateFlux(bandpass) == obj.calculateFlux(bandpass)
# == obj.drawImage(bandpass).array.sum()
# - .separable is a boolean indicating whether or not the profile can be factored into a
# spatial part and a spectral part.
# - .wave_list is a numpy array indicating wavelengths of particular interest, for instance, the
# wavelengths at which the sed is explicitly defined via a LookupTable. These are the
# wavelengths that will be used (in addition to those in bandpass.wave_list) when drawing an
# image of the chromatic profile.
# - .interpolated is a boolean indicating whether any part of the object hierarchy includes an
# InterpolatedChromaticObject.
# - .spectral indicates obj.sed.spectral
# - .dimensionless indicates obj.sed.dimensionless
def __init__(self, obj):
self._obj = obj
if not isinstance(obj, GSObject) and not isinstance(obj, ChromaticObject):
raise TypeError("Can only directly instantiate ChromaticObject with a GSObject "
"or ChromaticObject argument.")
self.separable = obj.separable
self.interpolated = obj.interpolated
self.deinterpolated = obj.deinterpolated
@property
def sed(self):
return self._obj.sed
@property
def SED(self):
from .deprecated import depr
depr('obj.SED', 2.5, 'obj.sed')
return self.sed
@property
def wave_list(self):
return self.sed.wave_list
@property
def gsparams(self):
"""The `GSParams` for this object.
"""
return self._obj.gsparams
@property
def redshift(self):
"""The redshift of the object.
"""
return self.sed.redshift
[docs] def withGSParams(self, gsparams=None, **kwargs):
"""Create a version of the current object with the given gsparams
Note: if this object wraps other objects (e.g. Convolution, Sum, Transformation, etc.)
those component objects will also have their gsparams updated to the new value.
"""
if gsparams == self.gsparams: return self
ret = copy.copy(self)
ret._obj = self._obj.withGSParams(gsparams, **kwargs)
return ret
@staticmethod
def _get_multiplier(sed, bandpass, wave_list):
"""Cached integral of product of sed and bandpass."""
wave_list = np.array(wave_list)
if len(wave_list) > 0:
bp = _LookupTable(wave_list, bandpass(wave_list), 'linear')
multiplier = bp.integrate_product(sed)
else:
multiplier = integ.int1d(lambda w: sed(w) * bandpass(w),
bandpass.blue_limit, bandpass.red_limit)
return multiplier
[docs] @staticmethod
def resize_multiplier_cache(maxsize):
"""Resize the cache (default size=10) containing the integral over the product of an `SED`
and a `Bandpass`, which is used by `ChromaticObject.drawImage`.
Parameters:
maxsize: The new number of products to cache.
"""
ChromaticObject._multiplier_cache.resize(maxsize)
def _fiducial_profile(self, bandpass):
"""
Return a fiducial achromatic profile of a chromatic object that can be used to estimate
default output image characteristics, or in the case of separable profiles, can be scaled to
give the monochromatic profile at any wavelength or the wavelength-integrated profile.
"""
bpwave = bandpass.effective_wavelength
bpwave, prof0 = self._approxWavelength(bpwave)
if prof0.flux != 0:
return bpwave, prof0
candidate_waves = np.concatenate(
[np.array([0.5 * (bandpass.blue_limit + bandpass.red_limit)]),
bandpass.wave_list,
self.wave_list])
# Prioritize wavelengths near the bandpass effective wavelength.
candidate_waves = candidate_waves[np.argsort(np.abs(candidate_waves - bpwave))]
for w in candidate_waves:
if bandpass.blue_limit <= w <= bandpass.red_limit:
prof0 = self.evaluateAtWavelength(w)
if prof0.flux != 0:
return w, prof0
raise GalSimError("Could not locate fiducial wavelength where SED * Bandpass is nonzero.")
def _approxWavelength(self, wave):
# If a class doesn't have any more appropriate choice, just use evaluateAtWavelength
# InterpolatedChromaticObject has a better choice when phot=True, so overrides this.
return wave, self.evaluateAtWavelength(wave)
def __eq__(self, other):
return (self is other or
(isinstance(other, ChromaticObject) and
hasattr(other, '_obj') and # not all ChromaticObjects have an _obj attribute.
self._obj == other._obj))
def __ne__(self, other): return not self.__eq__(other)
def __hash__(self): return hash(("galsim.ChromaticObject", self._obj))
def __repr__(self):
return 'galsim.ChromaticObject(%r)'%self._obj
def __str__(self):
return 'galsim.ChromaticObject(%s)'%self._obj
[docs] def interpolate(self, waves, **kwargs):
"""Build interpolation images to (possibly) speed up subsequent `drawImage` calls.
This method is used as a pre-processing step that can expedite image rendering using objects
that have to be built up as sums of `GSObject` instances with different parameters at each
wavelength, by interpolating between images at each wavelength instead of making a more
costly instantiation of the relevant `GSObject` at each value of wavelength at which the
bandpass is defined.
This routine does a costly initialization process to build up a grid of `Image` instances to
be used for the interpolation later on. However, the object can get reused with different
bandpasses, so there should not be any need to make many versions of this object, and there
is a significant savings each time it is drawn into an image.
As a general rule of thumb, chromatic objects that are separable do not benefit from this
particular optimization, whereas those that involve making `GSObject` instances with
wavelength-dependent keywords or transformations do benefit from it.
Note that the interpolation scheme is simple linear interpolation in wavelength, and no
extrapolation beyond the originally-provided range of wavelengths is permitted. However,
the overall flux at each wavelength will use the exact `SED` at that wavelength to give
more accurate final flux values. You can disable this feature by setting
``use_exact_sed = False``.
The speedup involved in using interpolation depends in part on the bandpass used for
rendering (since that determines how many full profile evaluations are involved in rendering
the image). However, for `ChromaticAtmosphere` with simple profiles like `Kolmogorov`, the
speedup in some simple example cases is roughly a factor of three, whereas for more
expensive to render profiles like the `ChromaticOpticalPSF`, the speedup is more typically a
factor of 10-50.
Achromatic transformations can be applied either before or after setting up interpolation,
with the best option depending on the application. For example, when rendering many times
with the same achromatic transformation applied, it is typically advantageous to apply the
transformation before setting up the interpolation. But there is no value in this when
applying different achromatic transformation to each object. Chromatic transformations
should be applied before setting up interpolation; attempts to render images of
`ChromaticObject` instances with interpolation followed by a chromatic transformation will
result in the interpolation being unset and the full calculation being done.
Because of the clever way that the `ChromaticConvolution` routine works, convolutions of
separable chromatic objects with non-separable ones that use interpolation will still
benefit from these optimizations. For example, a non-separable chromatic PSF that uses
interpolation, when convolved with a sum of two separable galaxy components each with their
own `SED`, will be able to take advantage of this optimization. In contrast, when
convolving two non-separable profiles that already have interpolation set up, there is no
way to take advantage of that interpolation optimization, so it will be ignored and the
full calculation will be done. However, interpolation can be set up for the convolution of
two non-separable profiles, after the convolution step. This could be beneficial for
example when convolving a chromatic optical PSF and chromatic atmosphere, before convolving
with multiple galaxy profiles.
For use cases requiring a high level of precision, we recommend a comparison between the
interpolated and the more accurate calculation for at least one case, to ensure that the
required precision has been reached.
The input parameter ``waves`` determines the input grid on which images are precomputed. It
is difficult to give completely general guidance as to how many wavelengths to choose or how
they should be spaced; some experimentation compared with the exact calculation is warranted
for each particular application. The best choice of settings might depend on how strongly
the parameters of the object depend on wavelength.
Parameters:
waves: The list, tuple, or NumPy array of wavelengths to be used when
building up the grid of images for interpolation. The wavelengths
should be given in nanometers, and they should span the full range
of wavelengths covered by any bandpass to be used for drawing an
`Image` (i.e., this class will not extrapolate beyond the given
range of wavelengths). They can be spaced any way the user likes,
not necessarily linearly, though interpolation will be linear in
wavelength between the specified wavelengths.
oversample_fac: Factor by which to oversample the stored profiles compared to the
default, which is to sample them at the Nyquist frequency for
whichever wavelength has the highest Nyquist frequency.
``oversample_fac``>1 results in higher accuracy but costlier
pre-computations (more memory and time). [default: 1]
use_exact_sed: If true, then rescale the interpolated image for a given wavelength
by the ratio of the exact `SED` at that wavelength to the linearly
interpolated `SED` at that wavelength. Thus, the flux of the
interpolated object should be correct, at the possible expense of
other features. [default: True]
Returns:
the version of the Chromatic object that uses interpolation
(This will be an `InterpolatedChromaticObject` instance.)
"""
if 'use_exact_SED' in kwargs:
from .deprecated import depr
depr('use_exact_SED', 2.5, 'use_exact_sed')
kwargs['use_exact_sed'] = kwargs.pop('use_exact_SED')
return InterpolatedChromaticObject(self, waves, **kwargs)
@property
def spectral(self):
"""Boolean indicating if the `ChromaticObject` has units compatible with a spectral density.
"""
return self.sed.spectral
@property
def dimensionless(self):
"""Boolean indicating if the `ChromaticObject` is dimensionless.
"""
return self.sed.dimensionless
@staticmethod
def _get_integrator(integrator, wave_list):
# Decide on integrator. If the user passed one of the integrators from galsim.integ, that's
# fine. Otherwise we decide based on the adopted integration rule and the presence/absence
# of `wave_list`.
if isinstance(integrator, str):
if integrator == 'quadratic':
rule = integ.quadRule
elif integrator == 'trapezoidal':
rule = integ.trapzRule
elif integrator == 'midpoint':
rule = integ.midptRule
else:
raise GalSimValueError("Unrecognized integration rule", integrator,
('trapezoidal', 'midpoint', 'quadratic'))
if len(wave_list) > 0:
integrator = integ.SampleIntegrator(rule)
else:
integrator = integ.ContinuousIntegrator(rule)
if not isinstance(integrator, integ.ImageIntegrator):
raise TypeError("Invalid type passed in for integrator!")
return integrator
[docs] def drawImage(self, bandpass, image=None, integrator='quadratic', **kwargs):
"""Base implementation for drawing an image of a `ChromaticObject`.
Some subclasses may choose to override this for specific efficiency gains. For instance,
most GalSim use cases will probably finish with a convolution, in which case
`ChromaticConvolution.drawImage` will be used.
The task of drawImage() in a chromatic context is to integrate a chromatic surface
brightness profile multiplied by the throughput of ``bandpass``, over the wavelength
interval indicated by ``bandpass``.
Several integrators are available in galsim.integ to do this integration when using the
first method (non-interpolated integration). By default, `galsim.integ.SampleIntegrator`
will be used if either ``bandpass.wave_list`` or ``self.wave_list`` have len() > 0.
If lengths of both are zero, which may happen if both the bandpass throughput and the SED
associated with ``self`` are analytic python functions for example, then
`galsim.integ.ContinuousIntegrator` will be used instead. This latter case by default will
evaluate the integrand at 250 equally-spaced wavelengths between ``bandpass.blue_limit``
and ``bandpass.red_limit``.
By default, the above two integrators will use the ``rule`` `galsim.integ.quadRule`
for integration. The midpoint rule for integration can be specified instead by passing an
integrator that has been initialized with the ``rule`` set to `galsim.integ.midptRule`.
When creating a `ContinuousIntegrator`, the number of samples ``N`` is also an argument.
For example::
>>> integrator = galsim.integ.ContinuousIntegrator(rule=galsim.integ.midptRule, N=100)
>>> image = chromatic_obj.drawImage(bandpass, integrator=integrator)
Finally, this method uses a cache to avoid recomputing the integral over the product of
the bandpass and object SED when possible (i.e., for separable profiles). Because the
cache size is finite, users may find that it is more efficient when drawing many images
to group images using the same SEDs and bandpasses together in order to hit the cache more
often. The default cache size is 10, but may be resized using the
`ChromaticObject.resize_multiplier_cache` method.
Parameters:
bandpass: A `Bandpass` object representing the filter against which to
integrate.
image: Optionally, the Image to draw onto. (See `GSObject.drawImage`
for details.) [default: None]
integrator: When doing the exact evaluation of the profile, this argument should
be one of the image integrators from galsim.integ, or a string
'trapezoidal', 'midpoint', or 'quadratic', in which case the routine
will use a `SampleIntegrator` or `ContinuousIntegrator` depending on
whether or not the object has a ``wave_list``. [default: 'quadratic',
which will try to select an appropriate integrator using the
quadratic integration rule automatically.]
**kwargs: For all other kwarg options, see `GSObject.drawImage`
Returns:
the drawn `Image`.
"""
# Store the last bandpass used and any extra kwargs.
self._last_bp = bandpass
if self.sed.dimensionless:
raise GalSimSEDError("Can only draw ChromaticObjects with spectral SEDs.", self.sed)
# setup output image using fiducial profile
wave0, prof0 = self._fiducial_profile(bandpass)
image = prof0.drawImage(image=image, setup_only=True, **kwargs)
_remove_setup_kwargs(kwargs)
# determine combined self.wave_list and bandpass.wave_list
wave_list, _, _ = utilities.combine_wave_list(self, bandpass)
# If there are photon ops, they'll probably need valid wavelengths, so add
# WavelengthSampler as the first op in the list (if one isn't already present).
if (kwargs.get('photon_ops', None)
and not any([isinstance(p,WavelengthSampler) for p in kwargs['photon_ops']])):
wave_sampler = WavelengthSampler(self.sed, bandpass)
kwargs['photon_ops'] = [wave_sampler] + kwargs['photon_ops']
if self.separable:
multiplier = ChromaticObject._multiplier_cache(self.sed, bandpass, tuple(wave_list))
prof0 *= multiplier/self.sed(wave0)
image = prof0.drawImage(image=image, **kwargs)
return image
# Note to developers, if we get to this point, then the implementaion is to integrate
# drawn images as a function of wavelength. This is quite slow.
# Moreover, for method=phot, it cannot implement some features like save_photons=True.
# The guards to avoid this happening for method=phot are mostly in the drawImage methods
# of subclasses. So if you find this happening, the solution is most likely to
# specialize something in the subclass's drawImage method.
assert kwargs.get('method', None) != 'phot'
integrator = self._get_integrator(integrator, wave_list)
# merge self.wave_list into bandpass.wave_list if using a sampling integrator
if isinstance(integrator, integ.SampleIntegrator):
if len(wave_list) < 2:
raise GalSimIncompatibleValuesError(
"Cannot use SampleIntegrator when Bandpass and SED are both analytic.",
integrator=integrator, bandpass=bandpass, sed=self.sed)
bandpass = Bandpass(_LookupTable(wave_list, bandpass(wave_list), 'linear'), 'nm')
add_to_image = kwargs.pop('add_to_image', False)
integral = integrator(self.evaluateAtWavelength, bandpass, image, kwargs)
# For performance profiling, store the number of evaluations used for the last integration
# performed. Note that this might not be very useful for ChromaticSum instances, which are
# drawn one profile at a time, and hence _last_n_eval will only represent the final
# component drawn.
self._last_n_eval = integrator.last_n_eval
# Apply integral to the initial image appropriately.
# Note: Don't do image = integral and return that for add_to_image==False.
# Remember that python doesn't actually do assignments, so this won't update the
# original image if the user provided one. The following procedure does work.
if not add_to_image:
image.setZero()
image += integral
self._last_wcs = image.wcs
return image
[docs] def drawKImage(self, bandpass, image=None, integrator='quadratic', **kwargs):
"""Base implementation for drawing the Fourier transform of a `ChromaticObject`.
The task of drawKImage() in a chromatic context is exactly analogous to the task of
`drawImage` in a chromatic context: to integrate the ``sed * bandpass`` weighted Fourier
profiles over wavelength.
See `drawImage` for details on integration options.
Parameters:
bandpass: A `Bandpass` object representing the filter against which to integrate.
image: If provided, this will be the complex `Image` onto which to draw the
k-space image. If ``image`` is None, then an automatically-sized image
will be created. If ``image`` is given, but its bounds are undefined,
then it will be resized appropriately based on the profile's size.
[default: None]
integrator: When doing the exact evaluation of the profile, this argument should be
one of the image integrators from galsim.integ, or a string
'trapezoidal', 'midpoint', or 'quadratic', in which case the routine will
use a `SampleIntegrator` or `ContinuousIntegrator` depending on whether or
not the object has a ``wave_list``. [default: 'quadratic', which will try
to select an appropriate integrator using the quadratic integration rule
automatically.]
**kwargs: For all other kwarg options, see `GSObject.drawKImage`.
Returns:
a complex `Image` instance (created if necessary)
"""
if self.sed.dimensionless:
raise GalSimSEDError("Can only drawK ChromaticObjects with spectral SEDs.", self.sed)
# setup output image (semi-arbitrarily using the bandpass effective wavelength)
prof0 = self.evaluateAtWavelength(bandpass.effective_wavelength)
image = prof0.drawKImage(image=image, setup_only=True, **kwargs)
_remove_setup_kwargs(kwargs)
# determine combined self.wave_list and bandpass.wave_list
wave_list, _ , _ = utilities.combine_wave_list(self, bandpass)
if self.separable:
multiplier = ChromaticObject._multiplier_cache(self.sed, bandpass, tuple(wave_list))
prof0 *= multiplier/self.sed(bandpass.effective_wavelength)
image = prof0.drawKImage(image=image, **kwargs)
return image
integrator = self._get_integrator(integrator, wave_list)
# merge self.wave_list into bandpass.wave_list if using a sampling integrator
if isinstance(integrator, integ.SampleIntegrator):
bandpass = Bandpass(_LookupTable(wave_list, bandpass(wave_list), 'linear'), 'nm')
add_to_image = kwargs.pop('add_to_image', False)
image_int = integrator(self.evaluateAtWavelength, bandpass, image, kwargs, doK=True)
# For performance profiling, store the number of evaluations used for the last integration
# performed. Note that this might not be very useful for ChromaticSum instances, which are
# drawn one profile at a time, and hence _last_n_eval will only represent the final
# component drawn.
self._last_n_eval = integrator.last_n_eval
# Apply integral to the initial image appropriately.
# Note: Don't do image = integral and return that for add_to_image==False.
# Remember that python doesn't actually do assignments, so this won't update the
# original image if the user provided one. The following procedure does work.
if add_to_image:
image += image_int
else:
image._copyFrom(image_int)
return image
[docs] def evaluateAtWavelength(self, wave):
"""Evaluate this chromatic object at a particular wavelength.
Parameters:
wave: Wavelength in nanometers.
Returns:
the monochromatic object at the given wavelength.
"""
# Subclasses all override this.
return self._obj.evaluateAtWavelength(wave)
def _shoot(self, photons, rng):
self._obj._shoot(photons, rng)
[docs] def applyTo(self, photon_array, local_wcs=None, rng=None):
"""Apply the chromatic profile as a convolution to an existing photon array.
This method allows instances of this class to duck type as a PhotonOp, so one can apply it
in a photon_ops list.
Parameters:
photon_array: A `PhotonArray` to apply the operator to.
local_wcs: A `LocalWCS` instance defining the local WCS for the current photon
bundle in case the operator needs this information. [default: None]
rng: A random number generator to use to effect the convolution.
[default: None]
"""
if not photon_array.hasAllocatedWavelengths():
raise GalSimError("Using ChromaticObject as a PhotonOp requires wavelengths be set")
p1 = PhotonArray(len(photon_array))
p1._copyFrom(photon_array, slice(None), slice(None), do_xy=False, do_flux=False)
obj = local_wcs.toImage(self) if local_wcs is not None else self
rng = BaseDeviate(rng)
obj._shoot(p1, rng)
photon_array.convolve(p1, rng)
# Make op* and op*= work to adjust the flux of the object
[docs] def __mul__(self, flux_ratio):
"""Scale the flux of the object by the given flux ratio, which may be an `SED`, a float, or
a univariate callable function (of wavelength in nanometers) that returns a float.
The normalization of a `ChromaticObject` is tracked through its ``.sed`` attribute, which
may have dimensions of either [photons/wavelength-interval/area/time/solid-angle] or
[1/solid-angle].
If ``flux_ratio`` is a spectral `SED` (i.e., ``flux_ratio.spectral==True``), then self.sed
must be dimensionless for dimensional consistency. The returned object will have a
spectral sed attribute. On the other hand, if ``flux_ratio`` is a dimensionless `SED`,
float, or univariate callable function, then the returned object will have ``.spectral``
and ``.dimensionless`` matching ``self.spectral`` and ``self.dimensionless``.
Parameters:
flux_ratio: The factor by which to scale the normalization of the object.
``flux_ratio`` may be a float, univariate callable function, in which
case the argument should be wavelength in nanometers and return value
the flux ratio for that wavelength, or an `SED`.
Returns:
a new object with scaled flux.
"""
return self.withScaledFlux(flux_ratio)
__rmul__ = __mul__
# Likewise for op/ and op/=
def __div__(self, other):
return self.__mul__(1./other)
__truediv__ = __div__
[docs] def withScaledFlux(self, flux_ratio):
"""Multiply the flux of the object by ``flux_ratio``
Parameters:
flux_ratio: The factor by which to scale the normalization of the object.
``flux_ratio`` may be a float, univariate callable function, in which
case the argument should be wavelength in nanometers and return value
the flux ratio for that wavelength, or an `SED`.
Returns:
a new object with scaled flux.
"""
return transform.Transform(self, flux_ratio=flux_ratio)
[docs] def withFlux(self, target_flux, bandpass):
"""Return a new `ChromaticObject` with flux through the `Bandpass` ``bandpass`` set to
``target_flux``.
Parameters:
target_flux: The desired flux normalization of the `ChromaticObject`.
bandpass: A `Bandpass` object defining a filter bandpass.
Returns:
the new normalized `ChromaticObject`.
"""
current_flux = self.calculateFlux(bandpass)
norm = target_flux/current_flux
return self * norm
[docs] def withMagnitude(self, target_magnitude, bandpass):
"""Return a new `ChromaticObject` with magnitude through ``bandpass`` set to
``target_magnitude``. Note that this requires ``bandpass`` to have been assigned a
zeropoint using `Bandpass.withZeropoint`.
Parameters:
target_magnitude: The desired magnitude of the `ChromaticObject`.
bandpass: A `Bandpass` object defining a filter bandpass.
Returns:
the new normalized `ChromaticObject`.
"""
if bandpass.zeropoint is None:
raise GalSimError("Cannot call ChromaticObject.withMagnitude on this bandpass, because"
" it does not have a zeropoint. See `Bandpass.withZeropoint`")
current_magnitude = self.calculateMagnitude(bandpass)
norm = 10**(-0.4*(target_magnitude - current_magnitude))
return self * norm
[docs] def withFluxDensity(self, target_flux_density, wavelength):
"""Return a new `ChromaticObject` with flux density set to ``target_flux_density`` at
wavelength ``wavelength``.
Parameters:
target_flux_density: The target normalization in photons/nm/cm^2/s.
wavelength: The wavelength, in nm, at which the flux density will be set.
Returns:
the new normalized `SED`.
"""
_photons = units.astrophys.photon/(units.s * units.cm**2 * units.nm)
if self.dimensionless:
raise GalSimSEDError("Cannot set flux density of dimensionless ChromaticObject.", self)
if isinstance(wavelength, units.Quantity):
wavelength_nm = wavelength.to(units.nm, units.spectral())
current_flux_density = self.sed(wavelength_nm.value)
else:
wavelength_nm = wavelength * units.nm
current_flux_density = self.sed(wavelength)
if isinstance(target_flux_density, units.Quantity):
target_flux_density = target_flux_density.to(
_photons, units.spectral_density(wavelength_nm)).value
factor = target_flux_density / current_flux_density
return self * factor
[docs] def atRedshift(self, redshift):
"""Create a version of the current object with a different redshift.
This will both adjust the SED to have the given redshift and set a ``redshift`` attribute
with the given value.
Returns:
the object with the new redshift
"""
from .deprecated import depr
depr('atRedshift', '2.5.3', 'SED.atRedshift',
'We now require that the SED of an object be appropriately redshifted before being '
'used in a chromatic object. So rather than `(obj * sed).atRedshift(z)`, you '
'should now use `obj * sed.atRedshift(z)`. For more complicated chromatic objects, '
'it is not always clear what parts of the wavelength dependence should be shifted '
'by this method, so we no longer allow the ambiguous syntax.')
return self._atRedshift(redshift)
def _atRedshift(self, redshift):
return ChromaticTransformation(self, _redshift=redshift)
[docs] def calculateCentroid(self, bandpass):
"""Determine the centroid of the wavelength-integrated surface brightness profile.
Parameters:
bandpass: The bandpass through which the observation is made.
Returns:
the centroid of the integrated surface brightness profile, as a PositionD.
"""
# if either the Bandpass or self maintain a wave_list, evaluate integrand only at
# those wavelengths.
if len(bandpass.wave_list) > 0 or len(self.wave_list) > 0:
w, _, _ = utilities.combine_wave_list(self, bandpass)
objs = [self.evaluateAtWavelength(ww) for ww in w]
fluxes = np.array([o.flux for o in objs])
centroids = [o.centroid for o in objs]
xcentroids = np.array([c.x for c in centroids])
ycentroids = np.array([c.y for c in centroids])
bp = _LookupTable(w, bandpass(w), 'linear')
flux = bp.integrate_product(_LookupTable(w, fluxes, 'linear'))
xcentroid = bp.integrate_product(_LookupTable(w, fluxes * xcentroids, 'linear')) / flux
ycentroid = bp.integrate_product(_LookupTable(w, fluxes * ycentroids, 'linear')) / flux
return _PositionD(xcentroid, ycentroid)
else:
flux_integrand = lambda w: self.evaluateAtWavelength(w).flux * bandpass(w)
def xcentroid_integrand(w):
mono = self.evaluateAtWavelength(w)
return mono.centroid.x * mono.flux * bandpass(w)
def ycentroid_integrand(w):
mono = self.evaluateAtWavelength(w)
return mono.centroid.y * mono.flux * bandpass(w)
flux = integ.int1d(flux_integrand, bandpass.blue_limit, bandpass.red_limit)
xcentroid = 1./flux * integ.int1d(xcentroid_integrand,
bandpass.blue_limit,
bandpass.red_limit)
ycentroid = 1./flux * integ.int1d(ycentroid_integrand,
bandpass.blue_limit,
bandpass.red_limit)
return _PositionD(xcentroid, ycentroid)
[docs] def calculateFlux(self, bandpass):
"""Return the flux (photons/cm^2/s) of the `ChromaticObject` through a `Bandpass` bandpass.
Parameters:
bandpass: A `Bandpass` object representing a filter, or None to compute the bolometric
flux. For the bolometric flux the integration limits will be set to
(0, infinity) unless overridden by non-None `SED` attributes ``blue_limit``
or ``red_limit``. Note that an `SED` defined using a `LookupTable`
automatically has ``blue_limit`` and ``red_limit`` set.
Returns:
the flux through the bandpass.
"""
if self.sed.dimensionless:
raise GalSimSEDError("Cannot calculate flux of dimensionless ChromaticObject.",
self.sed)
return self.sed.calculateFlux(bandpass)
[docs] def calculateMagnitude(self, bandpass):
"""Return the `ChromaticObject` magnitude through a `Bandpass` ``bandpass``.
Note that this requires ``bandpass`` to have been assigned a zeropoint using
`Bandpass.withZeropoint`.
Parameters:
bandpass: A `Bandpass` object representing a filter, or None to compute the
bolometric magnitude. For the bolometric magnitude the integration
limits will be set to (0, infinity) unless overridden by non-None `SED`
attributes ``blue_limit`` or ``red_limit``. Note that an `SED` defined
using a `LookupTable` automatically has ``blue_limit`` and ``red_limit``
set.
Returns:
the bandpass magnitude.
"""
if self.sed.dimensionless:
raise GalSimSEDError("Cannot calculate magnitude of dimensionless ChromaticObject.",
self.sed)
return self.sed.calculateMagnitude(bandpass)
# Add together ChromaticObjects and/or GSObjects
def __add__(self, other):
return ChromaticSum(self, other)
# Subtract ChromaticObjects and/or GSObjects
def __sub__(self, other):
return ChromaticSum(self, -other)
def __neg__(self):
return -1. * self
# Following functions work to apply affine transformations to a ChromaticObject.
#
[docs] def expand(self, scale):
"""Expand the linear size of the profile by the given (possibly wavelength-dependent)
scale factor ``scale``, while preserving surface brightness.
This doesn't correspond to either of the normal operations one would typically want to
do to a galaxy. The functions dilate() and magnify() are the more typical usage. But this
function is conceptually simple. It rescales the linear dimension of the profile, while
preserving surface brightness. As a result, the flux will necessarily change as well.
See dilate() for a version that applies a linear scale factor while preserving flux.
See magnify() for a version that applies a scale factor to the area while preserving surface
brightness.
Parameters:
scale: The factor by which to scale the linear dimension of the object. In
addition, ``scale`` may be a callable function, in which case the argument
should be wavelength in nanometers and the return value the scale for that
wavelength.
Returns:
the expanded object
"""
if hasattr(scale, '__call__'):
def buildScaleJac(w):
s = scale(w)
return np.array([[s,np.zeros_like(s)],[np.zeros_like(s),s]])
jac = buildScaleJac
else:
jac = None if scale == 1 else np.diag([scale, scale])
return transform.Transform(self, jac=jac)
[docs] def dilate(self, scale):
"""Dilate the linear size of the profile by the given (possibly wavelength-dependent)
``scale``, while preserving flux.
e.g. ``half_light_radius`` <-- ``half_light_radius * scale``
See expand() and magnify() for versions that preserve surface brightness, and thus
change the flux.
Parameters:
scale: The linear rescaling factor to apply. In addition, ``scale`` may be a
callable function, in which case the argument should be wavelength in
nanometers and the return value the scale for that wavelength.
Returns:
the dilated object.
"""
if hasattr(scale, '__call__'):
return self.expand(scale).withScaledFlux(lambda w: 1./scale(w)**2)
else:
return self.expand(scale).withScaledFlux(1./scale**2)
[docs] def magnify(self, mu):
"""Apply a lensing magnification, scaling the area and flux by ``mu`` at fixed surface
brightness.
This process applies a lensing magnification ``mu``, which scales the linear dimensions of the
image by the factor sqrt(mu), i.e., ``half_light_radius`` <-- ``half_light_radius * sqrt(mu)``
while increasing the flux by a factor of ``mu``. Thus, magnify() preserves surface
brightness.
See dilate() for a version that applies a linear scale factor while preserving flux.
Parameters:
mu: The lensing magnification to apply. In addition, ``mu`` may be a callable
function, in which case the argument should be wavelength in nanometers
and the return value the magnification for that wavelength.
Returns:
the magnified object.
"""
if hasattr(mu, '__call__'):
return self.expand(lambda w: np.sqrt(mu(w)))
else:
return self.expand(math.sqrt(mu))
[docs] def shear(self, *args, **kwargs):
"""Apply an area-preserving shear to this object, where arguments are either a `Shear`,
or arguments that will be used to initialize one.
For more details about the allowed keyword arguments, see the `Shear` docstring.
The shear() method precisely preserves the area. To include a lensing distortion with
the appropriate change in area, either use shear() with magnify(), or use lens(), which
combines both operations.
Note that, while gravitational shear is monochromatic, the shear method may be used for
many other use cases including some which may be wavelength-dependent, such as
intrinsic galaxy shape, telescope dilation, atmospheric PSF shape, etc. Thus, the
shear argument is allowed to be a function of wavelength like other transformations.
Parameters:
shear: The shear to be applied. Or, as described above, you may instead supply
parameters to construct a `Shear` directly. eg. ``obj.shear(g1=g1,g2=g2)``.
In addition, the ``shear`` parameter may be a callable function, in which
case the argument should be wavelength in nanometers and the return value
the shear for that wavelength, returned as a `galsim.Shear` instance.
Returns:
the sheared object.
"""
if len(args) == 1:
if kwargs:
raise TypeError("Gave both unnamed and named arguments!")
if not hasattr(args[0], '__call__') and not isinstance(args[0], Shear):
raise TypeError("Unnamed argument is not a Shear or function returning Shear!")
shear = args[0]
elif len(args) > 1:
raise TypeError("Too many unnamed arguments!")
elif 'shear' in kwargs:
# Need to break this out specially in case it is a function of wavelength
shear = kwargs.pop('shear')
if kwargs:
raise TypeError("Too many kwargs provided!")
else:
shear = Shear(**kwargs)
if hasattr(shear, '__call__'):
jac = lambda w: (shear(w).getMatrix()
if np.isscalar(w)
# Note: The .T switches from shape=(N,2,2) to (2,2,N)
# Technically it transposes each 2x2 matrix along the way, but
# the shear matrices are symmetric, so that doesn't matter.
else np.array([shear(ww).getMatrix() for ww in w])).T
else:
jac = shear.getMatrix()
return transform.Transform(self, jac=jac)
def _shear(self, shear):
"""Equivalent to `ChromaticObject.shear`, but only valid for a galsim.Shear object,
not any of the possible wavelength-dependent options.
Parameters:
shear: The `Shear` to be applied.
Returns:
the sheared object.
"""
return transform.Transform(self, shear.getMatrix())
[docs] def lens(self, g1, g2, mu):
"""Apply a lensing shear and magnification to this object.
This `ChromaticObject` method applies a lensing (reduced) shear and magnification.
The shear must be specified using the g1, g2 definition of shear (see `Shear` for details).
This is the same definition as the outputs of the `PowerSpectrum` and `NFWHalo` classes,
which compute shears according to some lensing power spectrum or lensing by an NFW dark
matter halo. The magnification determines the rescaling factor for the object area and
flux, preserving surface brightness.
While gravitational lensing is achromatic, we do allow the parameters ``g1``, ``g2``, and
``mu`` to be callable functions to be parallel to all the other transformations of
chromatic objects. In this case, the functions should take the wavelength in nanometers as
the argument, and the return values are the corresponding value at that wavelength.
Parameters:
g1: First component of lensing (reduced) shear to apply to the object.
g2: Second component of lensing (reduced) shear to apply to the object.
mu: Lensing magnification to apply to the object. This is the factor by which
the solid angle subtended by the object is magnified, preserving surface
brightness.
Returns:
the lensed object.
"""
if any(hasattr(g, '__call__') for g in (g1,g2)):
_g1 = g1
_g2 = g2
if not hasattr(g1, '__call__'): _g1 = lambda w: g1
if not hasattr(g2, '__call__'): _g2 = lambda w: g2
S = lambda w: Shear(g1=_g1(w), g2=_g2(w))
sheared = self.shear(S)
else:
sheared = self.shear(g1=g1,g2=g2)
return sheared.magnify(mu)
def _lens(self, g1, g2, mu):
"""Equivalent to `ChromaticObject.lens`, but without the overhead of some of the sanity
checks or any of the possible wavelength-dependent options.
Parameters:
g1: First component of lensing (reduced) shear to apply to the object.
g2: Second component of lensing (reduced) shear to apply to the object.
mu: Lensing magnification to apply to the object. This is the factor by which
the solid angle subtended by the object is magnified, preserving surface
brightness.
Returns:
the lensed object.
"""
shear = _Shear(g1 + 1j*g2)
return transform.Transform(self, shear.getMatrix() * math.sqrt(mu))
[docs] def rotate(self, theta):
"""Rotate this object by an `Angle` ``theta``.
Parameters:
theta: Rotation angle (`Angle` object, +ve anticlockwise). In addition, ``theta``
may be a callable function, in which case the argument should be wavelength
in nanometers and the return value the rotation angle for that wavelength,
returned as a `galsim.Angle` instance.
Returns:
the rotated object.
"""
if hasattr(theta, '__call__'):
def buildRMatrix(w):
sth = np.sin(theta(w))
cth = np.cos(theta(w))
R = np.array([[cth, -sth],
[sth, cth]], dtype=float)
return R
jac = buildRMatrix
else:
sth, cth = theta.sincos()
jac = np.array([[cth, -sth],
[sth, cth]], dtype=float)
return transform.Transform(self, jac=jac)
[docs] def shift(self, *args, **kwargs):
"""Apply a (possibly wavelength-dependent) (dx, dy) shift to this chromatic object.
For a wavelength-independent shift, you may supply ``dx,dy`` as either two arguments, as a
tuple, or as a PositionD or PositionI object.
For a wavelength-dependent shift, you may supply two functions of wavelength in nanometers
which will be interpreted as ``dx(wave)`` and ``dy(wave)``, or a single function of
wavelength in nanometers that returns either a 2-tuple, PositionD, or PositionI.
Parameters:
dx: Horizontal shift to apply (float or function).
dy: Vertical shift to apply (float or function).
Returns:
the shifted object.
"""
# This follows along the galsim.utilities.pos_args function, but has some
# extra bits to account for the possibility of dx,dy being functions.
# First unpack args/kwargs
if len(args) == 0:
# Then dx,dy need to be kwargs
# If not, then python will raise an appropriate error.
try:
dx = kwargs.pop('dx')
dy = kwargs.pop('dy')
except KeyError:
raise TypeError('shift() requires exactly 2 arguments (dx, dy)') from None
offset = None
elif len(args) == 1:
if hasattr(args[0], '__call__'):
try:
args[0](700.).x
# If the function returns a Position, recast it as a function returning
# a numpy array.
def offset_func(w):
d = args[0](w)
return np.asarray( (d.x, d.y) )
offset = offset_func
except AttributeError:
# Then it's a function returning a tuple or list or array.
# Just make sure it is actually an array to make our life easier later.
offset = lambda w: np.asarray(args[0](w))
elif isinstance(args[0], Position):
offset = np.asarray( (args[0].x, args[0].y) )
else:
# Let python raise the appropriate exception if this isn't valid.
dx, dy = args[0]
offset = np.asarray( (dx, dy) )
elif len(args) == 2:
dx = args[0]
dy = args[1]
offset = None
else:
raise TypeError("Too many arguments supplied!")
if kwargs:
raise TypeError("Got unexpected keyword arguments: %s",kwargs.keys())
if offset is None:
offset = utilities.functionize(lambda x,y:(x,y))(dx, dy)
return transform.Transform(self, offset=offset)
def _shift(self, dx, dy):
"""Equivalent to `ChromaticObject.shift`, but only valid for a scalar shift (dx, dy)
not any of the possible wavelength-dependent options.
Parameters:
dx: The x-component of the shift to apply
dy: The y-component of the shift to apply
Returns:
the shifted object.
"""
return transform.Transform(self, offset=_PositionD(dx,dy))
ChromaticObject._multiplier_cache = utilities.LRU_Cache(
ChromaticObject._get_multiplier, maxsize=10)
[docs]class InterpolatedChromaticObject(ChromaticObject):
"""A `ChromaticObject` that uses interpolation of predrawn images to speed up subsequent
rendering.
This class wraps another `ChromaticObject`, which is stored in the attribute ``deinterpolated``.
Any `ChromaticObject` can be used, although the interpolation procedure is most effective
for non-separable objects, which can sometimes be very slow to render.
Normally, you would not create an InterpolatedChromaticObject directly. It is the
return type from `ChromaticObject.interpolate`. See the description of that function
for more details.
Parameters:
original: The `ChromaticObject` to be interpolated.
waves: The list, tuple, or NumPy array of wavelengths to be used when
building up the grid of images for interpolation. The wavelengths
should be given in nanometers, and they should span the full range
of wavelengths covered by any bandpass to be used for drawing an `Image`
(i.e., this class will not extrapolate beyond the given range of
wavelengths). They can be spaced any way the user likes, not
necessarily linearly, though interpolation will be linear in
wavelength between the specified wavelengths.
oversample_fac: Factor by which to oversample the stored profiles compared to the
default, which is to sample them at the Nyquist frequency for
whichever wavelength has the highest Nyquist frequency.
``oversample_fac``>1 results in higher accuracy but costlier
pre-computations (more memory and time). [default: 1]
use_exact_sed: If true, then rescale the interpolated image for a given wavelength by
the ratio of the exact `SED` at that wavelength to the linearly
interpolated `SED` at that wavelength. Thus, the flux of the interpolated
object should be correct, at the possible expense of other features.
[default: True]
"""
def __init__(self, original, waves, oversample_fac=1.0, use_exact_sed=True,
use_exact_SED=None):
if use_exact_SED is not None:
from .deprecated import depr
depr('use_exact_SED', 2.5, 'use_exact_sed')
use_exact_sed = use_exact_SED
self.waves = np.sort(np.array(waves))
self.oversample = oversample_fac
self.use_exact_sed = use_exact_sed
self.separable = original.separable
self.interpolated = True
# Don't interpolate an interpolation. Go back to the original.
self.deinterpolated = original.deinterpolated
@property
def sed(self):
return self.deinterpolated.sed
# Note: We don't always need all of these, so only calculate them if needed.
# E.g. if photon shooting, pretty much only need objs.
@lazy_property
def objs(self):
return [ self.deinterpolated.evaluateAtWavelength(wave) for wave in self.waves ]
@lazy_property
def stepk_vals(self):
return np.array([ obj.stepk for obj in self.objs ])
@lazy_property
def maxk_vals(self):
return np.array([ obj.maxk for obj in self.objs ])
@lazy_property
def ims(self):
# Find the Nyquist scale for each, and to be safe, choose the minimum value to use for the
# array of images that is being stored.
nyquist_scale_vals = [ obj.nyquist_scale for obj in self.objs ]
scale = np.min(nyquist_scale_vals) / self.oversample
# Find the suggested image size for each object given the choice of scale, and use the
# maximum just to be safe.
possible_im_sizes = [ obj.getGoodImageSize(scale) for obj in self.objs ]
im_size = np.max(possible_im_sizes)
# Note that `no_pixel` is used (we want the object on its own, without a pixel response).
return [ obj.drawImage(scale=scale, nx=im_size, ny=im_size, method='no_pixel')
for obj in self.objs ]
@lazy_property
def fluxes(self):
return [ obj.flux for obj in self.objs ]
@property
def gsparams(self):
"""The `GSParams` for this object.
"""
return self.deinterpolated.gsparams
[docs] @doc_inherit
def withGSParams(self, gsparams=None, **kwargs):
if gsparams == self.gsparams: return self
ret = copy.copy(self)
ret.deinterpolated = self.deinterpolated.withGSParams(gsparams, **kwargs)
return ret
def __eq__(self, other):
return (self is other or
(isinstance(other, InterpolatedChromaticObject) and
self.deinterpolated == other.deinterpolated and
np.array_equal(self.waves, other.waves) and
self.oversample == other.oversample and
self.use_exact_sed == other.use_exact_sed))
def __hash__(self):
return hash(("galsim.InterpolatedChromaticObject", self.deinterpolated, tuple(self.waves),
self.oversample, self.use_exact_sed))
def __repr__(self):
s = 'galsim.InterpolatedChromaticObject(%r,%r'%(self.deinterpolated, self.waves)
if self.oversample != 1.0:
s += ', oversample_fac=%r'%self.oversample
if not self.use_exact_sed:
s += ', use_exact_sed=False'
s += ')'
return s
def __str__(self):
return 'galsim.InterpolatedChromaticObject(%s,%s)'%(self.deinterpolated, self.waves)
def _imageAtWavelength(self, wave):
"""
Get an image of the object at a particular wavelength, using linear interpolation between
the originally-stored images. Also returns values for step_k and max_k, to be used to
expedite the instantation of `InterpolatedImage`.
Parameters:
wave: Wavelength in nanometers.
Returns:
an `Image` of the object at the given wavelength.
"""
# First, some wavelength-related sanity checks.
if wave < self.waves[0] or wave > self.waves[-1]:
raise GalSimRangeError("Requested wavelength is outside the allowed range.",
wave, self.waves[0], self.waves[-1])
# Figure out where the supplied wavelength is compared to the list of wavelengths on which
# images were originally tabulated.
lower_idx, frac = _findWave(self.waves, wave)
# Actually do the interpolation for the image, stepk, and maxk.
im = _linearInterp(self.ims, frac, lower_idx)
stepk = _linearInterp(self.stepk_vals, frac, lower_idx)
maxk = _linearInterp(self.maxk_vals, frac, lower_idx)
# Rescale to use the exact flux or normalization if requested.
if self.use_exact_sed:
interp_norm = _linearInterp(self.fluxes, frac, lower_idx)
exact_norm = self.sed(wave)
im *= exact_norm/interp_norm
return im, stepk, maxk
def _approxWavelength(self, wave):
# More efficient to use one of the original objects, not a new InterpolatedImage.
k = np.searchsorted(self.waves, wave)
if k >= len(self.waves) or (k > 0 and wave-self.waves[k-1] < self.waves[k]-wave):
k = k - 1
return self.waves[k], self.objs[k]
[docs] def evaluateAtWavelength(self, wave):
"""
Evaluate this `ChromaticObject` at a particular wavelength using interpolation.
Parameters:
wave: Wavelength in nanometers.
Returns:
the monochromatic object at the given wavelength, as a `GSObject`.
"""
im, stepk, maxk = self._imageAtWavelength(wave)
return InterpolatedImage(im, _force_stepk=stepk, _force_maxk=maxk)
def _shoot(self, photons, rng):
w = photons.wavelength
if np.any((w < self.waves[0]) | (w > self.waves[-1])):
bad_waves = [w for w in photons.wavelength if w < self.waves[0] or w > self.waves[-1]]
raise GalSimRangeError("Shooting photons outside the interpolated wave_list",
bad_waves, self.waves[0], self.waves[-1])
k = np.searchsorted(self.waves, w)
k[k==0] = 1 # if k == 0, then w == min(waves). Using k=1 instead is fine for this.
#assert np.all(k > 0)
#assert np.all(k < len(self.waves))
# For each w, these are the wavelengthat that bracket w:
w0 = self.waves[k-1]
w1 = self.waves[k]
#assert np.all(w0 <= w)
#assert np.all(w <= w1)
# If we could get away with averaging photons shot at each wavelength,
# these would be relative fractions. So e.g. x = x0 f0 + x1 f1 would be the
# right weighted average to use.
f0 = (w1-w) / (w1-w0)
#f1 = (w-w0) / (w1-w0) (We don't need this quantity below.)
# Instead of averaging these, we can do the averaging probabilistically by selecting
# each photon with a probability equal to the relative weight we would have used
# in the average.
u = np.empty(len(photons))
UniformDeviate(rng).generate(u)
use_k = k - (u<f0).astype(int) # The second term is either 0 or 1.
# Draw photons from the saved profiles according to when we have selected to use each one.
for kk, obj in enumerate(self.objs):
use = (use_k == kk) # True for each photon where this is the object to shoot from
temp = PhotonArray(np.sum(use))
temp._copyFrom(photons, slice(None), use, do_xy=False, do_flux=False)
obj._shoot(temp, rng)
# It will have tried to shoot the right total flux. But that's not correct.
# Rescale it down by the fraction of the total flux we actually want in this set.
temp.scaleFlux(len(temp)/len(photons))
photons._copyFrom(temp, use, slice(None))
def _get_interp_image(self, bandpass, image=None, integrator='quadratic',
_flux_ratio=None, **kwargs):
if integrator == 'quadratic':
rule = integ.quadRule
elif integrator == 'trapezoidal':
rule = integ.trapzRule
elif integrator == 'midpoint':
rule = integ.midptRule
else:
raise GalSimValueError("Unrecognized integration rule", integrator,
('trapezoidal', 'midpoint', 'quadratic'))
if _flux_ratio is None:
_flux_ratio = lambda w: np.ones_like(w)
# _flux_ratio might be a constant float. For simplicity below, turn it into a function.
if not hasattr(_flux_ratio, '__call__'):
val = _flux_ratio
_flux_ratio = lambda w: np.full_like(w, val)
# setup output image (semi-arbitrarily using the bandpass effective wavelength).
# Note: we cannot just use self._imageAtWavelength, because that routine returns an image
# with whatever pixel scale was required to sample all the images properly. We want to set
# up an output image that has the requested pixel scale, which might change the image size
# and so on.
_, prof0 = self._fiducial_profile(bandpass)
image = prof0.drawImage(image=image, setup_only=True, **kwargs)
_remove_setup_kwargs(kwargs)
# determine combination of self.wave_list and bandpass.wave_list
wave_objs = [self, bandpass]
if isinstance(_flux_ratio, SED):
wave_objs += [_flux_ratio]
wave_list, _, _ = utilities.combine_wave_list(wave_objs)
if np.any((wave_list < self.waves[0]) | (wave_list > self.waves[-1])): # pragma: no cover
# MJ: I'm pretty sure it's impossible to hit this.
# But just in case I'm wrong, I'm leaving it here but with pragma: no cover.
bad_waves = [w for w in wave_list if w < self.waves[0] or w > self.waves[-1]]
raise GalSimRangeError("Requested wavelength is outside the allowed range.",
bad_waves, self.waves[0], self.waves[-1])
# weights are the weights to use at each of the given wavelengths for the integration.
weights = rule.calculateWeights(wave_list, bandpass)
# im_weights are the weights for the stored images.
im_weights = np.zeros(len(self.waves))
for w, wt in zip(wave_list, weights):
# Find where this is with respect to the wavelengths on which images are stored.
lower_idx, frac = _findWave(self.waves, w)
assert 0 <= lower_idx < len(self.waves)-1
# Rescale to use the exact flux or normalization if requested.
if self.use_exact_sed:
interp_norm = _linearInterp(self.fluxes, frac, lower_idx)
exact_norm = self.sed(w)
wt *= exact_norm/interp_norm
im_weights[lower_idx] += (1.-frac) * wt * _flux_ratio(w)
im_weights[lower_idx+1] += frac * wt * _flux_ratio(w)
# Do the integral as a weighted sum.
integral = sum(wt*im for wt,im in zip(im_weights, self.ims) if wt!=0)
# Get the stepk, maxk using the same weights
stepk = np.average(self.stepk_vals, weights=im_weights)
maxk = np.average(self.maxk_vals, weights=im_weights)
# Instantiate the InterpolatedImage, using these conservative stepk and maxk choices.
return InterpolatedImage(integral, _force_stepk=stepk, _force_maxk=maxk)
[docs] def drawImage(self, bandpass, image=None, integrator='quadratic', **kwargs):
"""Draw an image as seen through a particular bandpass using the stored interpolated
images at the specified wavelengths.
This integration will take place using interpolation between stored images that were
setup when the object was constructed. (See interpolate() for more details.)
Parameters:
bandpass: A `Bandpass` object representing the filter against which to
integrate.
image: Optionally, the `Image` to draw onto. (See `GSObject.drawImage`
for details.) [default: None]
integrator: The integration algorithm to use, given as a string. Either
'midpoint', 'trapezoidal', or 'quadratic' is allowed.
[default: 'quadratic']
**kwargs: For all other kwarg options, see `GSObject.drawImage`.
Returns:
the drawn `Image`.
"""
# Store the last bandpass used.
self._last_bp = bandpass
if self.sed.dimensionless:
raise GalSimSEDError("Can only draw ChromaticObjects with spectral SEDs.", self.sed)
int_im = self._get_interp_image(bandpass, image=image, integrator=integrator, **kwargs)
image = int_im.drawImage(image=image, **kwargs)
self._last_wcs = image.wcs
return image
[docs]class ChromaticAtmosphere(ChromaticObject):
"""A `ChromaticObject` implementing two atmospheric chromatic effects: differential
chromatic refraction (DCR) and wavelength-dependent seeing.
Due to DCR, blue photons land closer to the zenith than red photons. Kolmogorov turbulence
also predicts that blue photons get spread out more by the atmosphere than red photons,
specifically FWHM is proportional to wavelength^(-0.2). Both of these effects can be
implemented by wavelength-dependent shifts and dilations.
Since DCR depends on the zenith angle and the parallactic angle (which is the position angle of
the zenith measured from North through East) of the object being drawn, these must be specified
via keywords. There are four ways to specify these values:
1) explicitly provide ``zenith_angle`` as a keyword of type `Angle`, and
``parallactic_angle`` will be assumed to be 0 by default.
2) explicitly provide both ``zenith_angle`` and ``parallactic_angle`` as keywords of type
`Angle`.
3) provide the coordinates of the object ``obj_coord`` and the coordinates of the zenith
``zenith_coord`` as keywords of type `CelestialCoord`.
4) provide the coordinates of the object ``obj_coord`` as a `CelestialCoord`, the hour angle
of the object ``HA`` as an `Angle`, and the latitude of the observer ``latitude`` as an
`Angle`.
DCR also depends on temperature, pressure and water vapor pressure of the atmosphere. The
default values for these are expected to be appropriate for LSST at Cerro Pachon, Chile, but
they are broadly reasonable for most observatories.
Note that a ChromaticAtmosphere by itself is NOT the correct thing to use to draw an image of a
star. Stars (and galaxies too, of course) have an `SED` that is not flat. To draw a real star,
you should either multiply the ChromaticAtmosphere object by an `SED`, or convolve it with a
point source multiplied by an `SED`::
>>> psf = galsim.ChromaticAtmosphere(...)
>>> star = galsim.DeltaFunction() * psf_sed
>>> final_star = galsim.Convolve( [psf, star] )
>>> final_star.drawImage(bandpass = bp, ...)
Parameters:
base_obj: Fiducial PSF, equal to the monochromatic PSF at ``base_wavelength``
base_wavelength: Wavelength represented by the fiducial PSF, in nanometers.
scale_unit: Units used by base_obj for its linear dimensions.
[default: galsim.arcsec]
alpha: Power law index for wavelength-dependent seeing. [default: -0.2,
the prediction for Kolmogorov turbulence]
zenith_angle: `Angle` from object to zenith [default: 0]
parallactic_angle: Parallactic angle, i.e. the position angle of the zenith, measured
from North through East. [default: 0]
obj_coord: Celestial coordinates of the object being drawn as a
`CelestialCoord`. [default: None]
zenith_coord: Celestial coordinates of the zenith as a `CelestialCoord`.
[default: None]
HA: Hour angle of the object as an `Angle`. [default: None]
latitude: Latitude of the observer as an `Angle`. [default: None]
pressure: Air pressure in kiloPascals. [default: 69.328 kPa]
temperature: Temperature in Kelvins. [default: 293.15 K]
H2O_pressure: Water vapor pressure in kiloPascals. [default: 1.067 kPa]
"""
def __init__(self, base_obj, base_wavelength, scale_unit=None, **kwargs):
self.separable = False
self.interpolated = False
self.deinterpolated = self
self.base_obj = base_obj
self.base_wavelength = base_wavelength
self._gsparams = base_obj.gsparams
if scale_unit is None:
scale_unit = arcsec
elif isinstance(scale_unit, str):
scale_unit = AngleUnit.from_name(scale_unit)
self.scale_unit = scale_unit
self.alpha = kwargs.pop('alpha', -0.2)
self.zenith_angle, self.parallactic_angle, self.kw = dcr.parse_dcr_angles(**kwargs)
# Any remaining kwargs will get forwarded to galsim.dcr.get_refraction
# Check that they're valid
for kw in self.kw:
if kw not in ('temperature', 'pressure', 'H2O_pressure'):
raise TypeError("Got unexpected keyword: {0}".format(kw))
self.base_refraction = dcr.get_refraction(self.base_wavelength, self.zenith_angle,
**self.kw)
@property
def sed(self):
return self.base_obj.sed
@property
def gsparams(self):
"""The `GSParams` for this object.
"""
return self.base_obj.gsparams
[docs] @doc_inherit
def withGSParams(self, gsparams=None, **kwargs):
if gsparams == self.gsparams: return self
ret = copy.copy(self)
ret.base_obj = self.base_obj.withGSParams(gsparams, **kwargs)
return ret
def __eq__(self, other):
return (self is other or
(isinstance(other, ChromaticAtmosphere) and
self.base_obj == other.base_obj and
self.base_wavelength == other.base_wavelength and
self.alpha == other.alpha and
self.zenith_angle == other.zenith_angle and
self.parallactic_angle == other.parallactic_angle and
self.scale_unit == other.scale_unit and
self.kw == other.kw))
def __hash__(self):
return hash(("galsim.ChromaticAtmosphere", self.base_obj, self.base_wavelength,
self.alpha, self.zenith_angle, self.parallactic_angle, self.scale_unit,
frozenset(self.kw.items())))
def __repr__(self):
s = 'galsim.ChromaticAtmosphere(%r, base_wavelength=%r, alpha=%r'%(
self.base_obj, self.base_wavelength, self.alpha)
s += ', zenith_angle=%r, parallactic_angle=%r'%(self.zenith_angle, self.parallactic_angle)
s += ', scale_unit=%r'%(self.scale_unit)
for k,v in self.kw.items():
s += ', %s=%r'%(k,v)
s += ')'
return s
def __str__(self):
return 'galsim.ChromaticAtmosphere(%s, base_wavelength=%s, alpha=%s)'%(
self.base_obj, self.base_wavelength, self.alpha)
[docs] def build_obj(self):
"""Build a `ChromaticTransformation` object for this `ChromaticAtmosphere`.
We don't do this right away to help make `ChromaticAtmosphere` objects be picklable.
Building this is quite fast, so we do it on the fly in `evaluateAtWavelength` and
`ChromaticObject.drawImage`.
"""
def shift_fn(w):
shift_magnitude = dcr.get_refraction(w, self.zenith_angle, **self.kw)
shift_magnitude -= self.base_refraction
shift_magnitude = shift_magnitude * radians / self.scale_unit
sinp, cosp = self.parallactic_angle.sincos()
shift = (-shift_magnitude * sinp, shift_magnitude * cosp)
return shift
def jac_fn(w):
scale = (w/self.base_wavelength)**self.alpha
return np.diag([scale, scale])
flux_ratio = lambda w: (w/self.base_wavelength)**(-2.*self.alpha)
return ChromaticTransformation(self.base_obj, jac=jac_fn, offset=shift_fn,
flux_ratio=flux_ratio)
def _shoot(self, photons, rng):
# Start with the base PSF
self.base_obj._shoot(photons, rng)
w = photons.wavelength
# Apply the wavelength-dependent scaling
if self.alpha != 0.:
scale = (w/self.base_wavelength)**self.alpha
photons.x *= scale
photons.y *= scale
# Apply DCR
shift_magnitude = dcr.get_refraction(w, self.zenith_angle, **self.kw)
shift_magnitude -= self.base_refraction
shift_magnitude *= radians / self.scale_unit
sinp, cosp = self.parallactic_angle.sincos()
photons.x += -shift_magnitude * sinp
photons.y += shift_magnitude * cosp
[docs] def evaluateAtWavelength(self, wave):
"""Evaluate this chromatic object at a particular wavelength.
Parameters:
wave: Wavelength in nanometers.
Returns:
the monochromatic object at the given wavelength.
"""
return self.build_obj().evaluateAtWavelength(wave)
class SimpleChromaticTransformation(ChromaticTransformation):
"""A class for the simplest kind of chromatic object -- a GSObject times an SED.
This is a subclass of ChromaticTransformation, which just skips some calculations
that are unnecessary in this simple, but fairly common special case.
Parameters:
obj: The object to be transformed.
sed: An `SED` instance.
gsparams: An optional `GSParams` argument. See the docstring for `GSParams` for
details. [default: None]
propagate_gsparams: Whether to propagate gsparams to each of the components. This
is normally a good idea, but there may be use cases where one
would not want to do this. [default: True]
"""
def __init__(self, obj, sed=1., gsparams=None, propagate_gsparams=True):
self.chromatic = False
self.separable = True
self.interpolated = False
self._redshift = None
self._gsparams = GSParams.check(gsparams, obj.gsparams)
self._propagate_gsparams = propagate_gsparams
self._original = obj
self._jac = None
self._offset = (0.,0.)
self._flux_ratio = sed
if self._propagate_gsparams:
self._original = self._original.withGSParams(self._gsparams)
self.deinterpolated = self
@lazy_property
def sed(self):
return self._flux_ratio * self.original.flux
def _atRedshift(self, redshift):
return SimpleChromaticTransformation(self.original, self._flux_ratio.atRedshift(redshift),
self._gsparams, self._propagate_gsparams)
def __hash__(self):
if not hasattr(self, '_hash'):
self._hash = hash(("galsim.SimpleChromaticTransformation", self.original,
self._flux_ratio, self._gsparams, self._propagate_gsparams))
return self._hash
def __repr__(self):
return ('galsim.SimpleChromaticTransformation(%r, sed=%r, '
'gsparams=%r, propagate_gsparams=%r)')%(
self.original, self._flux_ratio,
self._gsparams, self._propagate_gsparams)
def __str__(self):
return str(self.original) + ' * ' + str(self.sed)
def __getstate__(self):
d = self.__dict__.copy()
d.pop('sed',None)
return d
def __setstate__(self, d):
self.__dict__ = d
def _getTransformations(self, wave):
flux_ratio = self._flux_ratio(wave)
return self._jac, self._offset, flux_ratio
def _shoot(self, photons, rng):
self._original._shoot(photons, rng)
wave = photons.wavelength
flux_ratio = self._flux_ratio(wave)
photons.flux *= flux_ratio
def drawImage(self, bandpass, image=None, integrator='quadratic', **kwargs):
"""
See `ChromaticObject.drawImage` for a full description.
Parameters:
bandpass: A `Bandpass` object representing the filter against which to
integrate.
image: Optionally, the `Image` to draw onto. (See `GSObject.drawImage`
for details.) [default: None]
integrator: When doing the exact evaluation of the profile, this argument should
be one of the image integrators from galsim.integ, or a string
'trapezoidal', 'midpoint', 'quadratic', in which case the routine will
use a `SampleIntegrator` or `ContinuousIntegrator` depending on whether
or not the object has a ``wave_list``. [default: 'quadratic',
which will try to select an appropriate integrator using the
quadratic integration rule automatically.]
If the object being transformed is an `InterpolatedChromaticObject`,
then ``integrator`` can only be a string, either 'midpoint',
'trapezoidal', or 'quadratic'.
**kwargs: For all other kwarg options, see `GSObject.drawImage`.
Returns:
the drawn `Image`.
"""
# Store the last bandpass used.
self._last_bp = bandpass
if self.sed.dimensionless:
raise GalSimSEDError("Can only draw ChromaticObjects with spectral SEDs.", self.sed)
image = ChromaticObject.drawImage(self, bandpass, image, integrator, **kwargs)
self._last_wcs = image.wcs
return image
[docs]class ChromaticSum(ChromaticObject):
"""A sum of several `ChromaticObject` and/or `GSObject` instances.
Any `GSObject` in the sum is assumed to have a flat `SED` with spectral density of 1
photon/s/cm**2/nm.
This is the type returned from `galsim.Add` if any of the objects are a `ChromaticObject`.
Typically, you do not need to construct a ChromaticSum object explicitly. Normally, you
would just use the + operator, which returns a ChromaticSum when used with chromatic objects::
>>> bulge = galsim.Sersic(n=3, half_light_radius=0.8) * bulge_sed
>>> disk = galsim.Exponential(half_light_radius=1.4) * disk_sed
>>> gal = bulge + disk
You can also use the `Add` factory function, which returns a ChromaticSum object if any of
the individual objects are chromatic::
>>> gal = galsim.Add([bulge,disk])
Parameters:
args: Unnamed args should be a list of objects to add.
gsparams: An optional `GSParams` argument. See the docstring for `GSParams` for
details. [default: None]
propagate_gsparams: Whether to propagate gsparams to each of the components. This
is normally a good idea, but there may be use cases where one
would not want to do this. [default: True]
"""
def __init__(self, *args, **kwargs):
# Check kwargs first:
gsparams = kwargs.pop("gsparams", None)
self._propagate_gsparams = kwargs.pop("propagate_gsparams", True)
# Make sure there is nothing left in the dict.
if kwargs:
raise TypeError("Got unexpected keyword argument(s): %s"%kwargs.keys())
if len(args) == 0:
# No arguments. Could initialize with an empty list but draw then segfaults. Raise an
# exception instead.
raise TypeError("Must provide at least one GSObject or ChromaticObject.")
elif len(args) == 1:
# 1 argument. Should be either a GSObject, ChromaticObject or a list of these.
if isinstance(args[0], (GSObject, ChromaticObject)):
args = [args[0]]
elif isinstance(args[0], list):
args = args[0]
else:
raise TypeError("Single input argument must be a GSObject, a ChromaticObject,"
" or list of them.")
# else args is already the list of objects
# Figure out what gsparams to use
if gsparams is None:
# If none is given, take the most restrictive combination from the obj_list.
self._gsparams = GSParams.combine([obj.gsparams for obj in args])
else:
# If something explicitly given, then use that.
self._gsparams = GSParams.check(gsparams)
self.interpolated = any(arg.interpolated for arg in args)
if self.interpolated:
self.deinterpolated = ChromaticSum([arg.deinterpolated for arg in args],
gsparams=self._gsparams)
else:
self.deinterpolated = self
# We can only add ChromaticObjects together if they're either all SED'd or all non-SED'd
dimensionless = all(a.dimensionless for a in args)
spectral = all(a.spectral for a in args)
if not (dimensionless or spectral):
raise GalSimIncompatibleValuesError(
"Cannot add dimensionless and spectral ChromaticObjects.", args=args)
# We always consider ChromaticSums to be inseparable.
# Note that separable groups can only be identified if the constituent objects have the
# *same* SED even though a proportional SED is mathematically sufficient for separability.
# It's basically impossible to identify if two SEDs are proportional (or even equal) unless
# they point to the same memory, and in practice any sums that are like this would
# almost certainly have different relative fluxes, so they would end up with different
# SED instances. This implies that there is no point wasting time trying to pull out
# separable groups.
self.separable = False
self._obj_list = []
for obj in args:
if self._propagate_gsparams:
obj = obj.withGSParams(self._gsparams)
self._obj_list.append(obj)
@lazy_property
def sed(self):
sed = self._obj_list[0].sed
for obj in self._obj_list[1:]:
sed += obj.sed
return sed
@property
def gsparams(self):
"""The `GSParams` for this object.
"""
return self._gsparams
@property
def obj_list(self):
"""The list of objects being added.
"""
return self._obj_list
[docs] @doc_inherit
def withGSParams(self, gsparams=None, **kwargs):
if gsparams == self.gsparams: return self
ret = copy.copy(self)
ret._gsparams = GSParams.check(gsparams, self.gsparams, **kwargs)
if self._propagate_gsparams:
ret._obj_list = [ obj.withGSParams(ret._gsparams) for obj in self.obj_list ]
return ret
def _atRedshift(self, redshift):
ret = copy.copy(self)
ret._obj_list = [ obj._atRedshift(redshift) for obj in self.obj_list ]
return ret
def __eq__(self, other):
return (self is other or
(isinstance(other, ChromaticSum) and
self.obj_list == other.obj_list and
self._gsparams == other._gsparams and
self._propagate_gsparams == other._propagate_gsparams))
def __hash__(self):
return hash(("galsim.ChromaticSum", tuple(self.obj_list), self._gsparams,
self._propagate_gsparams))
def __repr__(self):
return 'galsim.ChromaticSum(%r, gsparams=%r, propagate_gsparams=%r)'%(
self.obj_list, self._gsparams, self._propagate_gsparams)
def __str__(self):
str_list = [ str(obj) for obj in self.obj_list ]
return 'galsim.ChromaticSum([%s])'%', '.join(str_list)
def __getstate__(self):
d = self.__dict__.copy()
d.pop('sed',None)
return d
def __setstate__(self, d):
self.__dict__ = d
[docs] def evaluateAtWavelength(self, wave):
"""Evaluate this chromatic object at a particular wavelength ``wave``.
Parameters:
wave: Wavelength in nanometers.
Returns:
the monochromatic object at the given wavelength.
"""
return Add([obj.evaluateAtWavelength(wave) for obj in self.obj_list],
gsparams=self._gsparams, propagate_gsparams=self._propagate_gsparams)
def _shoot(self, photons, rng):
w = photons.wavelength
# We probabilistically choose a component for each photon based on the
# relative flux density of that component for the given wavelength.
comp_flux = np.array([obj.sed(w) for obj in self._obj_list])
total_flux = np.sum(comp_flux, axis=0)
prob = comp_flux / total_flux
cdf = np.cumsum(prob, axis=0)
u = np.empty(len(photons))
UniformDeviate(rng).generate(u)
use_k = np.sum((u >= cdf), axis=0)
n = len(self._obj_list)
use_k[use_k==n] = n-1 # This shouldn't be possible, but maybe with numerical rounding...
# Draw photons from the components.
for kk, obj in enumerate(self._obj_list):
use = (use_k == kk) # True for each photon where this is the object to shoot from
this_n = np.sum(use)
if this_n == 0:
continue
temp = PhotonArray(this_n)
temp._copyFrom(photons, slice(None), use, do_xy=False, do_flux=False)
obj._shoot(temp, rng)
photons._copyFrom(temp, use, slice(None))
[docs] def drawImage(self, bandpass, image=None, integrator='quadratic', **kwargs):
"""Slightly optimized draw method for `ChromaticSum` instances.
Draws each summand individually and add resulting images together. This might waste time if
two or more summands are separable and have the same `SED`, and another summand with a
different `SED` is also added, in which case the summands should be added together first and
the resulting `Sum` object can then be chromaticized. In general, however, drawing
individual sums independently can help with speed by identifying chromatic profiles that
are separable into spectral and spatial factors.
Parameters:
bandpass: A `Bandpass` object representing the filter against which to
integrate.
image: Optionally, the `Image` to draw onto. (See `GSObject.drawImage`
for details.) [default: None]
integrator: When doing the exact evaluation of the profile, this argument should
be one of the image integrators from galsim.integ, or a string
'trapezoidal', 'midpoint', 'quadratic', in which case the routine will
use a `SampleIntegrator` or `ContinuousIntegrator` depending on whether
or not the object has a ``wave_list``. [default: 'quadratic',
which will try to select an appropriate integrator using the
quadratic integration rule automatically.]
**kwargs: For all other kwarg options, see `GSObject.drawImage`.
Returns:
the drawn `Image`.
"""
# Store the last bandpass used.
self._last_bp = bandpass
if self.sed.dimensionless:
raise GalSimSEDError("Can only draw ChromaticObjects with spectral SEDs.", self.sed)
add_to_image = kwargs.pop('add_to_image', False)
n_photons = kwargs.pop('n_photons', None)
save_photons = kwargs.get('save_photons', False)
all_photons = []
if n_photons is None or kwargs.get('method', None) != 'phot':
# Use given add_to_image for the first one, then add_to_image=True for the rest.
image = self.obj_list[0].drawImage(
bandpass, image=image, add_to_image=add_to_image, **kwargs)
_remove_setup_kwargs(kwargs)
added_flux = image.added_flux
if save_photons:
all_photons.append(image.photons)
for obj in self.obj_list[1:]:
image = obj.drawImage(bandpass, image=image, add_to_image=True, **kwargs)
added_flux += image.added_flux
if save_photons:
all_photons.append(image.photons)
image.added_flux = added_flux
if save_photons:
image.photons = PhotonArray.concatenate(all_photons)
else:
# If n_photons is specified, need some special handling here.
rng = BaseDeviate(kwargs.get('rng', None))
total_flux = self.calculateFlux(bandpass)
remaining_nphot = n_photons
remaining_flux = total_flux
if kwargs.pop('poisson_flux',False):
pd = PoissonDeviate(rng, total_flux)
flux_ratio = pd() / total_flux
else:
flux_ratio = 1
added_flux = 0
first = True
for i, obj in enumerate(self.obj_list):
this_flux = obj.calculateFlux(bandpass)
if i == len(self.obj_list)-1:
this_nphot = remaining_nphot
else:
bd = BinomialDeviate(rng, remaining_nphot, this_flux / remaining_flux)
this_nphot = int(bd())
remaining_nphot -= this_nphot
remaining_flux -= this_flux
# Get the flux right based on the actual draw and possibly poisson_flux.
if this_nphot > 0:
obj *= (this_nphot * total_flux) / (n_photons * this_flux) * flux_ratio
image = obj.drawImage(bandpass, image=image, add_to_image=add_to_image,
n_photons=this_nphot, poisson_flux=False, **kwargs)
added_flux += image.added_flux
if save_photons:
all_photons.append(image.photons)
if first:
# Note: This might not be i==0.
# Do this after the first call we make to drawImage.
_remove_setup_kwargs(kwargs)
add_to_image = True
first = False
if not remaining_flux > 0:
break
image.added_flux = added_flux
if save_photons:
image.photons = PhotonArray.concatenate(all_photons)
self._last_wcs = image.wcs
return image
[docs] def withScaledFlux(self, flux_ratio):
"""Multiply the flux of the object by ``flux_ratio``
Parameters:
flux_ratio: The factor by which to scale the flux.
Returns:
the object with the new flux.
"""
new_obj = ChromaticSum([ obj.withScaledFlux(flux_ratio) for obj in self.obj_list ])
if hasattr(self, 'covspec'):
new_covspec = self.covspec * flux_ratio**2
new_obj.covspec = new_covspec
return new_obj
[docs]class ChromaticConvolution(ChromaticObject):
"""A convolution of several `ChromaticObject` and/or `GSObject` instances.
Any `GSObject` in the convolution is assumed to have a flat `SED` with spectral density of 1
photon/s/cm**2/nm.
This is the type returned from `galsim.Convolve` if any of the objects is a `ChromaticObject`.
The normal way to use this class is to use the `Convolve` factory function::
>>> gal = galsim.Sersic(n, half_light_radius) * galsim.SED(sed_file, 'nm', 'flambda')
>>> psf = galsim.ChromaticAtmosphere(...)
>>> final = galsim.Convolve([gal, psf])
The objects to be convolved may be provided either as multiple unnamed arguments (e.g.
``Convolve(psf, gal, pix)``) or as a list (e.g. ``Convolve([psf, gal, pix])``). Any number of
objects may be provided using either syntax. (Well, the list has to include at least 1 item.)
Parameters:
args: Unnamed args should be a list of objects to convolve.
real_space: Whether to use real space convolution. [default: None, which means
to automatically decide this according to whether the objects have hard
edges.]
gsparams: An optional `GSParams` argument. See the docstring for `GSParams` for
details. [default: None]
propagate_gsparams: Whether to propagate gsparams to each of the components. This
is normally a good idea, but there may be use cases where one
would not want to do this. [default: True]
"""
def __init__(self, *args, **kwargs):
# First check for number of arguments != 0
if len(args) == 0:
# No arguments. Could initialize with an empty list but draw then segfaults. Raise an
# exception instead.
raise TypeError("Must provide at least one GSObject or ChromaticObject")
elif len(args) == 1:
if isinstance(args[0], (GSObject, ChromaticObject)):
args = [args[0]]
elif isinstance(args[0], list):
args = args[0]
else:
raise TypeError("Single input argument must be a GSObject, or a ChromaticObject,"
" or list of them.")
# else args is already the list of objects
# Check kwargs
# real space convolution is not implemented for chromatic objects.
real_space = kwargs.pop("real_space", None)
if real_space:
raise GalSimNotImplementedError(
"Real space convolution of chromatic objects not implemented.")
gsparams = kwargs.pop("gsparams", None)
self._propagate_gsparams = kwargs.pop("propagate_gsparams", True)
# Figure out what gsparams to use
if gsparams is None:
# If none is given, take the most restrictive combination from the obj_list.
self._gsparams = GSParams.combine([obj.gsparams for obj in args])
else:
# If something explicitly given, then use that.
self._gsparams = GSParams.check(gsparams)
# Make sure there is nothing left in the dict.
if kwargs:
raise TypeError("Got unexpected keyword argument(s): %s"%kwargs.keys())
# Accumulate convolutant .seds. Check if more than one is spectral.
nspectral = sum(arg.spectral for arg in args)
if nspectral > 1:
raise GalSimIncompatibleValuesError(
"Cannot convolve more than one spectral ChromaticObject.", args=args)
self._obj_list = []
# Unfold convolution of convolution.
for obj in args:
if self._propagate_gsparams:
obj = obj.withGSParams(self._gsparams)
if isinstance(obj, ChromaticConvolution):
self._obj_list.extend(obj.obj_list)
else:
self._obj_list.append(obj)
self.separable = all(obj.separable for obj in self._obj_list)
self.interpolated = any(obj.interpolated for obj in self._obj_list)
if self.interpolated:
self.deinterpolated = ChromaticConvolution(
[obj.deinterpolated for obj in self._obj_list],
gsparams=self._gsparams, propagate_gsparams=self._propagate_gsparams)
else:
self.deinterpolated = self
# Check quickly whether we are convolving two non-separable things that aren't
# ChromaticSums, >1 of which uses interpolation. If so, emit a warning that the
# interpolation optimization is being ignored and full evaluation is necessary.
# For the case of ChromaticSums, as long as each object in the sum is separable (even if the
# entire object is not) then interpolation can still be used. So we do not warn about this
# here.
n_nonsep = 0
n_interp = 0
for obj in self._obj_list:
if not obj.separable and not isinstance(obj, ChromaticSum): n_nonsep += 1
if obj.interpolated: n_interp += 1
if n_nonsep>1 and n_interp>0:
galsim_warn(
"Image rendering for this convolution cannot take advantage of "
"interpolation-related optimization. Will use full profile evaluation.")
@lazy_property
def sed(self):
sed = self._obj_list[0].sed
for obj in self._obj_list[1:]:
sed *= obj.sed
return sed
@property
def gsparams(self):
"""The `GSParams` for this object.
"""
return self._gsparams
@property
def obj_list(self):
"""The list of objects being convolved.
"""
return self._obj_list
[docs] @doc_inherit
def withGSParams(self, gsparams=None, **kwargs):
if gsparams == self.gsparams: return self
ret = copy.copy(self)
ret._gsparams = GSParams.check(gsparams, self.gsparams, **kwargs)
if self._propagate_gsparams:
ret._obj_list = [ obj.withGSParams(ret._gsparams) for obj in self.obj_list ]
return ret
@staticmethod
def _get_effective_prof(insep_obj, bandpass, iimult, integrator, gsparams):
# Find scale at which to draw effective profile
# Use smallest nyquist scale among the fiducial profile and at the two limits of the bp.
_, prof0 = insep_obj._fiducial_profile(bandpass)
prof1 = insep_obj.evaluateAtWavelength(bandpass.red_limit)
prof2 = insep_obj.evaluateAtWavelength(bandpass.blue_limit)
iiscale = min(prof0.nyquist_scale, prof1.nyquist_scale, prof2.nyquist_scale)
iiscale /= 2 # This seems to be required to make test_monochromatic_sed to pass.
# Not sure why, since I thought straight nyquist should be good enough.
# But if it's needed there, it's probably worth always doing, rather than
# having that test use iimult=2. And definitions of Nyquist are somewhat
# confusing, so it's possible that we should expect to need a factor of
# 2 smaller than nyquist for the pixel scale. :-S
if iimult is not None:
iiscale /= iimult
# Prevent infinite recursive loop by using ChromaticObject.drawImage() on a
# ChromaticConvolution.
if isinstance(insep_obj, ChromaticConvolution):
effective_prof_image = ChromaticObject.drawImage(
insep_obj, bandpass, scale=iiscale,
integrator=integrator, method='no_pixel')
else:
effective_prof_image = insep_obj.drawImage(
bandpass, scale=iiscale, integrator=integrator,
method='no_pixel')
return InterpolatedImage(effective_prof_image, gsparams=gsparams)
[docs] @staticmethod
def resize_effective_prof_cache(maxsize):
"""Resize the cache containing effective profiles.
These are wavelength-integrated products of separable profile SEDs, inseparable profiles,
and Bandpasses) used by `ChromaticConvolution.drawImage`.
Parameters:
maxsize: The new number of effective profiles to cache.
"""
ChromaticConvolution._effective_prof_cache.resize(maxsize)
def __eq__(self, other):
return (self is other or
(isinstance(other, ChromaticConvolution) and
self.obj_list == other.obj_list and
self._gsparams == other._gsparams and
self._propagate_gsparams == other._propagate_gsparams))
def __hash__(self):
return hash(("galsim.ChromaticConvolution", tuple(self.obj_list), self._gsparams,
self._propagate_gsparams))
def __repr__(self):
return 'galsim.ChromaticConvolution(%r, gsparams=%r, propagate_gsparams=%r)'%(
self.obj_list, self._gsparams, self._propagate_gsparams)
def __str__(self):
str_list = [ str(obj) for obj in self.obj_list ]
return 'galsim.ChromaticConvolution([%s])'%', '.join(str_list)
def __getstate__(self):
d = self.__dict__.copy()
d.pop('sed',None)
return d
def __setstate__(self, d):
self.__dict__ = d
def _approxWavelength(self, wave):
# If any of the components prefer a different wavelength, use that for all.
achrom_objs = []
for k, obj in enumerate(self.obj_list):
new_wave, aobj = obj._approxWavelength(wave)
if new_wave != wave:
# Break the loop and use evaluateAtWavelength for everything else.
achrom_objs = ([o.evaluateAtWavelength(new_wave) for o in self.obj_list[:k]] +
[aobj] +
[o.evaluateAtWavelength(new_wave) for o in self.obj_list[k+1:]])
break
else:
achrom_objs.append(aobj)
return new_wave, convolve.Convolve(achrom_objs, gsparams=self._gsparams,
propagate_gsparams=self._propagate_gsparams)
[docs] def evaluateAtWavelength(self, wave):
"""Evaluate this chromatic object at a particular wavelength ``wave``.
Parameters:
wave: Wavelength in nanometers.
Returns:
the monochromatic object at the given wavelength.
"""
return convolve.Convolve([obj.evaluateAtWavelength(wave) for obj in self.obj_list],
gsparams=self._gsparams, propagate_gsparams=self._propagate_gsparams)
def _shoot(self, photons, rng):
raise GalSimNotImplementedError("ChromaticConvolution cannot be used as a PhotonOp")
[docs] def drawImage(self, bandpass, image=None, integrator='quadratic', iimult=None, **kwargs):
"""Optimized draw method for the `ChromaticConvolution` class.
Works by finding sums of profiles which include separable portions, which can then be
integrated before doing any convolutions, which are pushed to the end.
This method uses a cache to avoid recomputing 'effective' profiles, which are the
wavelength-integrated products of inseparable profiles, the spectral components of
separable profiles, and the bandpass. Because the cache size is finite, users may find
that it is more efficient when drawing many images to group images using the same
SEDs, bandpasses, and inseparable profiles (generally PSFs) together in order to hit the
cache more often. The default cache size is 10, but may be resized using the
`ChromaticConvolution.resize_effective_prof_cache` method.
Parameters:
bandpass: A `Bandpass` object representing the filter against which to
integrate.
image: Optionally, the `Image` to draw onto. (See `GSObject.drawImage`
for details.) [default: None]
integrator: When doing the exact evaluation of the profile, this argument should
be one of the image integrators from galsim.integ, or a string
'trapezoidal', 'midpoint', or 'quadratic', in which case the routine
will use a `SampleIntegrator` or `ContinuousIntegrator` depending on
whether or not the object has a ``wave_list``. [default: 'quadratic',
which will try to select an appropriate integrator using the
quadratic integration rule automatically.]
iimult: Oversample any intermediate `InterpolatedImage` created to hold
effective profiles by this amount. [default: None]
**kwargs: For all other kwarg options, see `GSObject.drawImage`.
Returns:
the drawn `Image`.
"""
# Store the last bandpass used.
self._last_bp = bandpass
if self.sed.dimensionless:
raise GalSimSEDError("Can only draw ChromaticObjects with spectral SEDs.", self.sed)
# `ChromaticObject.drawImage()` can just as efficiently handle separable cases.
if self.separable:
image = ChromaticObject.drawImage(self, bandpass, image=image, **kwargs)
self._last_wcs = image.wcs
return image
# Now split up any `ChromaticSum`s:
# This is the tricky part. Some notation first:
# int(f(x,y,lambda)) denotes the integral over wavelength of chromatic surface
# brightness profile f(x,y,lambda).
# (f1 * f2) denotes the convolution of surface brightness profiles f1 & f2.
# (f1 + f2) denotes the addition of surface brightness profiles f1 & f2.
#
# In general, chromatic s.b. profiles can be classified as either separable or inseparable,
# depending on whether they can be factored into spatial and spectral components or not.
# Write separable profiles as g(x,y) * h(lambda), and leave inseparable profiles as
# f(x,y,lambda).
# We will suppress the arguments `x`, `y`, `lambda`, hereforward, but generally an `f`
# refers to an inseparable profile, a `g` refers to the spatial part of a separable
# profile, and an `h` refers to the spectral part of a separable profile.
#
# Now, analyze a typical scenario, a bulge+disk galaxy model (each of which is separable,
# e.g., an SED times an exponential profile for the disk, and a different SED times a DeV
# profile for the bulge). Suppose the PSF is inseparable. (Chromatic PSF's will generally
# be inseparable since we usually think of the spatial part of the PSF being normalized to
# unit integral for any fixed wavelength.) Say there's also an achromatic pixel to
# convolve with.
# The formula for this might look like:
#
# img = int((bulge + disk) * PSF * pix)
# = int((g1 h1 + g2 h2) * f3 * g4) # note pix is lambda-independent
# = int(g1 h1 * f3 * g4 + g2 h2 * f3 * g4) # distribute the + over the *
# = int(g1 h1 * f3 * g4) + int(g2 h2 * f3 * g4) # distribute the + over the int
# = g1 * g4 * int(h1 f3) + g2 * g4 * int(h2 f3) # move lambda-indep terms out of int
#
# The result is that the integral is now inside the convolution, meaning we only have to
# compute two convolutions instead of a convolution for each wavelength at which we evaluate
# the integrand. This technique, making an "effective" PSF profile for each of the bulge
# and disk, is a significant time savings in most cases.
#
# In general, we make effective profiles by splitting up ChromaticSum items and collecting
# the inseparable terms on which to do integration first, and then finish with convolution
# last.
phot = kwargs.get('method', 'auto') == 'phot'
# This optimization is not actually helpful when photon shooting.
if not phot:
# Here is the logic to turn
# int((g1 h1 + g2 h2) * f3)
# -> g1 * int(h1 f3) + g2 * int(h2 f3)
for i, obj in enumerate(self.obj_list):
if isinstance(obj, ChromaticSum):
# say obj.obj_list = [A,B,C], where obj is a ChromaticSum object
# Assemble temporary list of convolutants excluding the ChromaticSum in question.
tmplist = list(self.obj_list)
del tmplist[i] # remove ChromaticSum object from obj_list
tmplist.append(obj.obj_list[0]) # Append first summand, i.e., A, to convolutants
# now draw this image
tmpobj = ChromaticConvolution(tmplist)
add_to_image = kwargs.pop('add_to_image', False)
image = tmpobj.drawImage(bandpass, image=image, integrator=integrator,
iimult=iimult, add_to_image=add_to_image, **kwargs)
# Now add in the rest of the summands in turn, i.e., B and C
for summand in obj.obj_list[1:]:
tmplist = list(self.obj_list)
del tmplist[i]
tmplist.append(summand)
tmpobj = ChromaticConvolution(tmplist)
# add to previously started image
_remove_setup_kwargs(kwargs)
image = tmpobj.drawImage(bandpass, image=image, integrator=integrator,
iimult=iimult, add_to_image=True, **kwargs)
# Return the image here, breaking the loop early. If there are two ChromaticSum
# instances in obj_list, then the above procedure will repeat in the recursion,
# effectively distributing the multiplication over both sums.
self._last_wcs = image.wcs
return image
# If program gets this far, the objects in obj_list should be atomic (non-ChromaticSum
# and non-ChromaticConvolution). (The latter case was dealt with in the constructor.)
# setup output image (semi-arbitrarily using the bandpass effective wavelength)
wave0, prof0 = self._fiducial_profile(bandpass)
image = prof0.drawImage(image=image, setup_only=True, **kwargs)
_remove_setup_kwargs(kwargs)
# If we are photon shooting, then we can move all non-spectral objects to the photon_ops
# list and deal with them that way. This both more accurate and more efficient for most
# chromatic PSFs.
if phot:
psfs = [obj for obj in self.obj_list if obj.dimensionless]
gals = [obj for obj in self.obj_list if obj.spectral]
assert len(gals) == 1 # Should have been checked by constructor.
gal = gals[0]
kwargs['photon_ops'] = psfs + kwargs.get('photon_ops', [])
# Need to calculate n_photons now using the fiducial profile, not gal, in case the
# PSF has an interpolated image (e.g. OpticalPSF) which needs more photons.
flux = self.calculateFlux(bandpass)
prof1 = prof0.withFlux(flux)
n_photons = kwargs.pop('n_photons', 0)
poisson_flux = kwargs.pop('poisson_flux', n_photons == 0.)
max_extra_noise = kwargs.pop('max_extra_noise', 0.)
rng = BaseDeviate(kwargs.get('rng', None))
n_photons, g = prof1._calculate_nphotons(n_photons, poisson_flux, max_extra_noise, rng)
gal *= g
return gal.drawImage(bandpass, image=image, integrator=integrator,
n_photons=n_photons, **kwargs)
# Separate convolutants into a Convolution of inseparable profiles multiplied by the
# wavelength-dependent normalization of separable profiles, and the achromatic part of
# separable profiles.
insep_obj = [obj for obj in self.obj_list if not obj.separable]
# Note that len(insep_obj) > 0, since purely separable ChromaticConvolutions were
# already handled above.
# Don't wrap in Convolution if not needed. Single item can draw itself better than
# Convolution can.
if len(insep_obj) == 1:
insep_obj = insep_obj[0]
else:
insep_obj = convolve.Convolve(insep_obj, gsparams=self._gsparams,
propagate_gsparams=self._propagate_gsparams)
sep_profs = []
for obj in self.obj_list:
if not obj.separable:
continue
wave0, prof0 = obj._fiducial_profile(bandpass)
sep_profs.append(prof0 / obj.sed(wave0))
insep_obj *= obj.sed
# Collapse inseparable profiles and chromatic normalizations into one effective profile
# Note that at this point, insep_obj.sed should *not* be None.
effective_prof = ChromaticConvolution._effective_prof_cache(
insep_obj, bandpass, iimult, integrator, self._gsparams)
# append effective profile to separable profiles (which should all be GSObjects)
sep_profs.append(effective_prof)
# finally, convolve and draw.
final_prof = convolve.Convolve(sep_profs, gsparams=self._gsparams,
propagate_gsparams=self._propagate_gsparams)
image = final_prof.drawImage(image=image, **kwargs)
self._last_wcs = image.wcs
return image
@lazy_property
def noise(self):
"""An estimate of the noise already in the profile.
"""
# Condition for being able to propagate noise:
# Exactly one of the convolutants has a .covspec attribute.
covspecs = [ obj.covspec for obj in self.obj_list if hasattr(obj, 'covspec') ]
if len(covspecs) != 1:
raise GalSimError("Cannot compute noise for ChromaticConvolution for which number "
"of convolutants with covspec attribute is not 1.")
if not hasattr(self, '_last_bp'):
raise GalSimError("Cannot compute noise for ChromaticConvolution until after drawImage "
"has been called.")
covspec = covspecs[0]
other = convolve.Convolve([obj for obj in self.obj_list if not hasattr(obj, 'covspec')])
return covspec.toNoise(self._last_bp, other, self._last_wcs) # rng=?
ChromaticConvolution._effective_prof_cache = utilities.LRU_Cache(
ChromaticConvolution._get_effective_prof, maxsize=10)
[docs]class ChromaticDeconvolution(ChromaticObject):
"""A class for deconvolving a `ChromaticObject`.
The ChromaticDeconvolution class represents a wavelength-dependent deconvolution kernel.
You may also specify a gsparams argument. See the docstring for `GSParams` for more
information about this option. Note: if ``gsparams`` is unspecified (or None), then the
ChromaticDeconvolution instance inherits the same `GSParams` as the object being deconvolved.
This is the type returned from `galsim.Deconvolve` if the argument is a `ChromaticObject`.
This is the normal way to construct this class.
Parameters:
obj: The object to deconvolve.
gsparams: An optional `GSParams` argument. See the docstring for `GSParams` for
details. [default: None]
propagate_gsparams: Whether to propagate gsparams to each of the components. This
is normally a good idea, but there may be use cases where one
would not want to do this. [default: True]
"""
def __init__(self, obj, gsparams=None, propagate_gsparams=True):
if not obj.sed.dimensionless:
raise GalSimSEDError("Cannot deconvolve by spectral ChromaticObject.", obj.sed)
self._gsparams = GSParams.check(gsparams, obj.gsparams)
self._propagate_gsparams = propagate_gsparams
if self._propagate_gsparams:
self._obj = obj.withGSParams(self._gsparams)
else:
self._obj = obj
self.separable = obj.separable
self.interpolated = obj.interpolated
if self.interpolated:
self.deinterpolated = ChromaticDeconvolution(self._obj.deinterpolated, self._gsparams,
self._propagate_gsparams)
else:
self.deinterpolated = self
@property
def sed(self):
return SED(lambda w: self._obj.sed(w)**-1, 'nm', '1')
@property
def gsparams(self):
"""The `GSParams` for this object.
"""
return self._gsparams
[docs] @doc_inherit
def withGSParams(self, gsparams=None, **kwargs):
if gsparams == self.gsparams: return self
ret = copy.copy(self)
ret._gsparams = GSParams.check(gsparams, self.gsparams, **kwargs)
if self._propagate_gsparams:
ret._obj = self._obj.withGSParams(ret._gsparams)
return ret
def __eq__(self, other):
return (self is other or
(isinstance(other, ChromaticDeconvolution) and
self._obj == other._obj and
self.gsparams == other.gsparams and
self._propagate_gsparams == other._propagate_gsparams))
def __hash__(self):
return hash(("galsim.ChromaticDeconvolution", self._obj, self.gsparams,
self._propagate_gsparams))
def __repr__(self):
return 'galsim.ChromaticDeconvolution(%r, gsparams=%r, propagate_gsparams=%r)'%(
self._obj, self.gsparams, self._propagate_gsparams)
def __str__(self):
return 'galsim.ChromaticDeconvolution(%s)'%self._obj
[docs] def evaluateAtWavelength(self, wave):
"""Evaluate this chromatic object at a particular wavelength ``wave``.
Parameters:
wave: Wavelength in nanometers.
Returns:
the monochromatic object at the given wavelength.
"""
return convolve.Deconvolve(self._obj.evaluateAtWavelength(wave), gsparams=self.gsparams,
propagate_gsparams=self._propagate_gsparams)
def _shoot(self, photons, rng):
raise GalSimNotImplementedError("ChromaticDeconvolution cannot use method='phot'")
[docs]class ChromaticAutoConvolution(ChromaticObject):
"""A special class for convolving a `ChromaticObject` with itself.
It is equivalent in functionality to ``galsim.Convolve([obj,obj])``, but takes advantage of
the fact that the two profiles are the same for some efficiency gains.
This is the type returned from `galsim.AutoConvolve` if the argument is a `ChromaticObject`.
This is the normal way to construct this class.
Parameters:
obj: The object to be convolved with itself.
real_space: Whether to use real space convolution. [default: None, which means
to automatically decide this according to whether the objects have hard
edges.]
gsparams: An optional `GSParams` argument. See the docstring for `GSParams` for
details. [default: None]
propagate_gsparams: Whether to propagate gsparams to each of the components. This
is normally a good idea, but there may be use cases where one
would not want to do this. [default: True]
"""
def __init__(self, obj, real_space=None, gsparams=None, propagate_gsparams=True):
if not obj.sed.dimensionless:
raise GalSimSEDError("Cannot autoconvolve spectral ChromaticObject.", obj.sed)
self._gsparams = GSParams.check(gsparams, obj.gsparams)
self._propagate_gsparams = propagate_gsparams
if self._propagate_gsparams:
self._obj = obj.withGSParams(self._gsparams)
else:
self._obj = obj
self._real_space = real_space
self.separable = obj.separable
self.interpolated = obj.interpolated
if self.interpolated:
self.deinterpolated = ChromaticAutoConvolution(self._obj.deinterpolated, real_space,
self._gsparams, self._propagate_gsparams)
else:
self.deinterpolated = self
@lazy_property
def sed(self):
return self._obj.sed * self._obj.sed
@property
def gsparams(self):
"""The `GSParams` for this object.
"""
return self._gsparams
[docs] @doc_inherit
def withGSParams(self, gsparams=None, **kwargs):
if gsparams == self.gsparams: return self
ret = copy.copy(self)
ret._gsparams = GSParams.check(gsparams, self.gsparams, **kwargs)
if self._propagate_gsparams:
ret._obj = self._obj.withGSParams(ret._gsparams)
return ret
def __eq__(self, other):
return (self is other or
(isinstance(other, ChromaticAutoConvolution) and
self._obj == other._obj and
self._real_space == other._real_space and
self.gsparams == other.gsparams and
self._propagate_gsparams == other._propagate_gsparams))
def __hash__(self):
return hash(("galsim.ChromaticAutoConvolution", self._obj, self._real_space, self.gsparams,
self._propagate_gsparams))
def __repr__(self):
return ('galsim.ChromaticAutoConvolution(%r, real_space=%r, gsparams=%r, '
'propagate_gsparams=%r)')%(
self._obj, self._real_space, self.gsparams, self._propagate_gsparams)
def __str__(self):
return 'galsim.ChromaticAutoConvolution(%s)'%self._obj
def __getstate__(self):
d = self.__dict__.copy()
d.pop('sed',None)
return d
def __setstate__(self, d):
self.__dict__ = d
[docs] def evaluateAtWavelength(self, wave):
"""Evaluate this chromatic object at a particular wavelength ``wave``.
Parameters:
wave: Wavelength in nanometers.
Returns:
the monochromatic object at the given wavelength.
"""
return convolve.AutoConvolve(self._obj.evaluateAtWavelength(wave), self._real_space,
self._gsparams, self._propagate_gsparams)
def _shoot(self, photons, rng):
raise GalSimNotImplementedError("ChromaticAutoConvolution cannot be used as a PhotonOp")
[docs]class ChromaticAutoCorrelation(ChromaticObject):
"""A special class for correlating a `ChromaticObject` with itself.
It is equivalent in functionality to::
galsim.Convolve([obj,obj.rotate(180.*galsim.degrees)])
but takes advantage of the fact that the two profiles are the same for some efficiency gains.
This is the type returned from `galsim.AutoCorrelate` if the argument is a `ChromaticObject`.
This is the normal way to construct this class.
Parameters:
obj: The object to be convolved with itself.
real_space: Whether to use real space convolution. [default: None, which means
to automatically decide this according to whether the objects have hard
edges.]
gsparams: An optional `GSParams` argument. See the docstring for `GSParams` for
details. [default: None]
propagate_gsparams: Whether to propagate gsparams to each of the components. This
is normally a good idea, but there may be use cases where one
would not want to do this. [default: True]
"""
def __init__(self, obj, real_space=None, gsparams=None, propagate_gsparams=True):
if not obj.sed.dimensionless:
raise GalSimSEDError("Cannot autocorrelate spectral ChromaticObject.", obj.sed)
self._gsparams = GSParams.check(gsparams, obj.gsparams)
self._propagate_gsparams = propagate_gsparams
if self._propagate_gsparams:
self._obj = obj.withGSParams(self._gsparams)
else:
self._obj = obj
self._real_space = real_space
self.separable = obj.separable
self.interpolated = obj.interpolated
if self.interpolated:
self.deinterpolated = ChromaticAutoCorrelation(self._obj.deinterpolated,
self._real_space, self._gsparams,
self._propagate_gsparams)
else:
self.deinterpolated = self
@lazy_property
def sed(self):
return self._obj.sed * self._obj.sed
@property
def gsparams(self):
"""The `GSParams` for this object.
"""
return self._gsparams
[docs] @doc_inherit
def withGSParams(self, gsparams=None, **kwargs):
if gsparams == self.gsparams: return self
ret = copy.copy(self)
ret._gsparams = GSParams.check(gsparams, self.gsparams, **kwargs)
if self._propagate_gsparams:
ret._obj = self._obj.withGSParams(ret._gsparams)
return ret
def __eq__(self, other):
return (self is other or
(isinstance(other, ChromaticAutoCorrelation) and
self._obj == other._obj and
self._real_space == other._real_space and
self.gsparams == other.gsparams and
self._propagate_gsparams == other._propagate_gsparams))
def __hash__(self):
return hash(("galsim.ChromaticAutoCorrelation", self._obj, self._real_space, self.gsparams,
self._propagate_gsparams))
def __repr__(self):
return ('galsim.ChromaticAutoCorrelation(%r, real_space=%r, gsparams=%r, '
'propagate_gsparams=%r)')%(
self._obj, self._real_space, self.gsparams, self._propagate_gsparams)
def __str__(self):
return 'galsim.ChromaticAutoCorrelation(%s)'%self._obj
def __getstate__(self):
d = self.__dict__.copy()
d.pop('sed',None)
return d
def __setstate__(self, d):
self.__dict__ = d
[docs] def evaluateAtWavelength(self, wave):
"""Evaluate this chromatic object at a particular wavelength ``wave``.
Parameters:
wave: Wavelength in nanometers.
Returns:
the monochromatic object at the given wavelength.
"""
return convolve.AutoCorrelate(self._obj.evaluateAtWavelength(wave), self._real_space,
self.gsparams, self._propagate_gsparams)
def _shoot(self, photons, rng):
raise GalSimNotImplementedError("ChromaticAutoCorrelation cannot be used as a PhotonOp")
[docs]class ChromaticFourierSqrtProfile(ChromaticObject):
"""A class for computing the Fourier-space square root of a `ChromaticObject`.
The ChromaticFourierSqrtProfile class represents a wavelength-dependent Fourier-space square
root of a profile.
You may also specify a gsparams argument. See the docstring for `GSParams` for more
information about this option. Note: if ``gsparams`` is unspecified (or None), then the
ChromaticFourierSqrtProfile inherits the same `GSParams` as the object being operated on.
The normal way to use this class is to use the `FourierSqrt` factory function::
>>> fourier_sqrt = galsim.FourierSqrt(chromatic_obj)
If ``chromatic_obj`` is indeed a `ChromaticObject`, then that function will create a
ChromaticFourierSqrtProfile object.
Parameters:
obj: The object to compute the Fourier-space square root of.
gsparams: An optional `GSParams` argument. See the docstring for `GSParams` for
details. [default: None]
propagate_gsparams: Whether to propagate gsparams to each of the components. This
is normally a good idea, but there may be use cases where one
would not want to do this. [default: True]
"""
def __init__(self, obj, gsparams=None, propagate_gsparams=True):
if not obj.sed.dimensionless:
raise GalSimSEDError("Cannot take Fourier sqrt of spectral ChromaticObject.", obj.sed)
self._gsparams = GSParams.check(gsparams, obj.gsparams)
self._propagate_gsparams = propagate_gsparams
if self._propagate_gsparams:
self._obj = obj.withGSParams(self._gsparams)
else:
self._obj = obj
self.separable = obj.separable
self.interpolated = obj.interpolated
if self.interpolated:
self.deinterpolated = ChromaticFourierSqrtProfile(
self._obj.deinterpolated, self._gsparams, self._propagate_gsparams)
else:
self.deinterpolated = self
@property
def sed(self):
return SED(lambda w: self._obj.sed(w)**0.5, 'nm', '1')
@property
def gsparams(self):
"""The `GSParams` for this object.
"""
return self._gsparams
[docs] @doc_inherit
def withGSParams(self, gsparams=None, **kwargs):
if gsparams == self.gsparams: return self
ret = copy.copy(self)
ret._gsparams = GSParams.check(gsparams, self.gsparams, **kwargs)
if self._propagate_gsparams:
ret._obj = self._obj.withGSParams(ret._gsparams)
return ret
def __eq__(self, other):
return (self is other or
(isinstance(other, ChromaticFourierSqrtProfile) and
self._obj == other._obj and
self.gsparams == other.gsparams and
self._propagate_gsparams == other._propagate_gsparams))
def __hash__(self):
return hash(("galsim.ChromaticFourierSqrtProfile", self._obj, self.gsparams,
self._propagate_gsparams))
def __repr__(self):
return 'galsim.ChromaticFourierSqrtProfile(%r, gsparams=%r, propagate_gsparams=%r)'%(
self._obj, self.gsparams, self._propagate_gsparams)
def __str__(self):
return 'galsim.ChromaticFourierSqrtProfile(%s)'%self._obj
[docs] def evaluateAtWavelength(self, wave):
"""Evaluate this chromatic object at a particular wavelength ``wave``.
Parameters:
wave: Wavelength in nanometers.
Returns:
the monochromatic object at the given wavelength.
"""
return fouriersqrt.FourierSqrt(self._obj.evaluateAtWavelength(wave), self.gsparams,
self._propagate_gsparams)
def _shoot(self, photons, rng):
raise GalSimNotImplementedError("ChromaticFourierSqrtProfile cannot use method='phot'")
[docs]class ChromaticOpticalPSF(ChromaticObject):
"""A subclass of ChromaticObject meant to represent chromatic optical PSFs.
Chromaticity plays two roles in optical PSFs. First, it determines the diffraction limit, via
the wavelength/diameter factor. Second, aberrations such as defocus, coma, etc. are typically
defined in physical distances, but their impact on the PSF depends on their size in units of
wavelength. Other aspects of the optical PSF do not require explicit specification of their
chromaticity, e.g., once the obscuration and struts are specified in units of the aperture
diameter, their chromatic dependence gets taken care of automatically. Note that the
ChromaticOpticalPSF implicitly defines diffraction limits in units of ``scale_units``, which by
default are arcsec, but can in principle be set to any of our GalSim angle units.
When using interpolation to speed up image rendering (see the `ChromaticObject.interpolate`
method for details), the ideal number of wavelengths to use across a given bandpass depends on
the application and accuracy requirements. In general it will be necessary to do a test in
comparison with a more exact calculation to ensure convergence. However, a typical calculation
might use ~10-15 samples across a typical optical bandpass, with ``oversample_fac`` in the range
1.5-2; for moderate accuracy, ~5 samples across the bandpass and ``oversample_fac=1`` may
suffice. All of these statements assume that aberrations are not very large (typically <~0.25
waves, which is commonly satisfied by space telescopes); if they are larger than that, then more
stringent settings are required.
Note that a ChromaticOpticalPSF by itself is NOT the correct thing to use to draw an image of a
star. Stars (and galaxies too, of course) have an `SED` that is not flat. To draw a real star,
you should either multiply the ChromaticOpticalPSF object by an `SED`, or convolve it with a
point source multiplied by an `SED`::
>>> psf = galsim.ChromaticOpticalPSF(...)
>>> star = galsim.DeltaFunction() * psf_sed
>>> final_star = galsim.Convolve( [psf, star] )
>>> final_star.drawImage(bandpass = bp, ...)
.. note::
When geometric_shooting is False (the default), the photon shooting implementation is
only approximately correct with respect to the wavelength dependence. It is also
not particularly fast, since it generates three optical screens to span the wavelength
range and shoots from these (with a subsequent adjustment to improve the accuracy
of this approximation). We expect that most users who want to use photon shooting in
conjunction with this class will prefer to make an InterpolatedChromaticObject
(by calling ``psf.interpolate(...)``), especially if it is a good approximation to
use the same optical PSF for a whole exposure or CCD image, so the setup time for
the interpolation is able to be amortized for many objects.
Parameters:
lam: Fiducial wavelength for which diffraction limit and aberrations are
initially defined, in nanometers.
diam: Telescope diameter in meters. Either ``diam`` or ``lam_over_diam`` must be
specified.
lam_over_diam: Ratio of (fiducial wavelength) / telescope diameter in units of
``scale_unit``. Either ``diam`` or ``lam_over_diam`` must be specified.
aberrations: An array of aberrations, in units of fiducial wavelength ``lam``. The
size and format of this array is described in the OpticalPSF docstring.
scale_unit: Units used to define the diffraction limit and draw images.
[default: galsim.arcsec]
gsparams: An optional `GSParams` argument. See the docstring for `GSParams` for
details. [default: None]
geometric_shooting: If True, then when drawing using photon shooting, use geometric
optics approximation where the photon angles are derived from the
phase screen gradient. If False, then first draw using Fourier
optics and then shoot from the derived InterpolatedImage. [default: False]
**kwargs: Any other keyword arguments to be passed to OpticalPSF, for example,
related to struts, obscuration, oversampling, etc. See OpticalPSF
docstring for a complete list of options.
"""
_req_params = { 'lam' : float }
_opt_params = { k:v for k,v in OpticalPSF._opt_params.items() if k != 'diam' }
_single_params = [ {'diam' : float, 'lam_over_diam' : float} ]
def __init__(self, lam, diam=None, lam_over_diam=None, aberrations=None,
scale_unit=None, gsparams=None, **kwargs):
# First, take the basic info.
if scale_unit is None:
scale_unit = arcsec
elif isinstance(scale_unit, str):
scale_unit = AngleUnit.from_name(scale_unit)
self.scale_unit = scale_unit
self._gsparams = GSParams.check(gsparams)
# We have to require either diam OR lam_over_diam:
if ( (diam is None and lam_over_diam is None) or
(diam is not None and lam_over_diam is not None) ):
raise GalSimIncompatibleValuesError(
"Need to specify telescope diameter OR wavelength/diam ratio",
diam=diam, lam_over_diam=lam_over_diam)
if diam is not None:
self.lam_over_diam = (1.e-9*lam/diam)*radians/self.scale_unit
self.diam = diam
else:
self.lam_over_diam = lam_over_diam
self.diam = (lam*1e-9/lam_over_diam)*radians/self.scale_unit
self.lam = lam
if aberrations is not None:
self.aberrations = np.asarray(aberrations)
if len(self.aberrations) < 12:
self.aberrations = np.append(self.aberrations, [0] * (12-len(self.aberrations)))
else:
self.aberrations = np.zeros(12)
# Pop named aberrations from kwargs so aberrations=[0,0,0,0,1] means the same as
# defocus=1 (w/ all other named aberrations 0).
for i, ab in enumerate(['tip', 'tilt', 'defocus', 'astig1', 'astig2', 'coma1', 'coma2',
'trefoil1', 'trefoil2', 'spher']):
if ab in kwargs:
self.aberrations[i+2] = kwargs.pop(ab)
self.fft_sign = kwargs.pop('fft_sign', '+')
if self.fft_sign not in ['+', '-']:
raise GalSimValueError("Invalid fft_sign", self.fft_sign, allowed_values=['+','-'])
self.geometric_shooting = kwargs.pop('geometric_shooting', False)
# All kwargs left should be relevant for the Aperture.
# If the pupil plane image is given, then we can make the aperture now to use at all
# wavelengths. Otherwise, we need to make it anew for each wavelength.
if 'aper' in kwargs:
# Then the aperture is given explicitly.
self._aper = kwargs.pop('aper')
if kwargs:
raise TypeError("Invalid kwargs provided in conjunction with aper: %s."%(
tuple(kwargs.keys())))
self._kwargs = {}
elif 'pupil_plane_im' in kwargs:
# Then the aperture is not wavelength dependent. Make it now.
self._aper = Aperture(self.diam, lam=self.lam, gsparams=self._gsparams, **kwargs)
self._kwargs = {}
else:
self._aper = None
self._kwargs = kwargs
# This will be the stepk and maxk values at self.lam. But wait until the first time
# we call evaluateAtWavelength to compute them.
self._stepk = self._maxk = None
# Define the necessary attributes for this ChromaticObject.
self.separable = False
self.interpolated = False
self.deinterpolated = self
@property
def sed(self):
return SED(1, 'nm', '1')
@property
def gsparams(self):
"""The `GSParams` for this object.
"""
return self._gsparams
[docs] @doc_inherit
def withGSParams(self, gsparams=None, **kwargs):
if gsparams == self.gsparams: return self
ret = copy.copy(self)
ret._gsparams = GSParams.check(gsparams, self.gsparams, **kwargs)
return ret
def __eq__(self, other):
return (self is other or
(isinstance(other, ChromaticOpticalPSF) and
self.lam == other.lam and
self.lam_over_diam == other.lam_over_diam and
np.array_equal(self.aberrations, other.aberrations) and
self.scale_unit == other.scale_unit and
self.gsparams == other.gsparams and
self.fft_sign == other.fft_sign and
self.geometric_shooting == other.geometric_shooting and
self._aper == other._aper and
self._kwargs == other._kwargs))
def __hash__(self):
return hash(("galsim.ChromaticOpticalPSF", self.lam, self.lam_over_diam,
tuple(self.aberrations), self.scale_unit, self.gsparams,
self.fft_sign, self.geometric_shooting, self._aper,
frozenset(self._kwargs.items())))
def __repr__(self):
s = 'galsim.ChromaticOpticalPSF(lam=%r, lam_over_diam=%r, aberrations=%r'%(
self.lam, self.lam_over_diam, self.aberrations.tolist())
if self.scale_unit != arcsec:
s += ', scale_unit=%r'%self.scale_unit
if self.fft_sign == '-':
s += ', fft_sign="-"'
if self.geometric_shooting:
s += ', geometric_shooting=True'
if self._aper is not None:
s += ', aper=%r'%self._aper
else:
for k,v in self._kwargs.items():
s += ', %s=%r'%(k,v)
s += ', gsparams=%r'%self.gsparams
s += ')'
return s
def __str__(self):
return 'galsim.ChromaticOpticalPSF(lam=%s, lam_over_diam=%s, aberrations=%s)'%(
self.lam, self.lam_over_diam, self.aberrations.tolist())
[docs] def evaluateAtWavelength(self, wave):
"""
Method to directly instantiate a monochromatic instance of this object.
Parameters:
wave: Wavelength in nanometers.
"""
if self._aper is not None:
# If aperture is not wavelength-dependent (i.e. given as an image) then we can use
# the same aperture for each wavelength and save making gratuitous copies of a
# large numpy array, which will be identical each time.
aper = self._aper
else:
# Otherwise we need to make the apeture anew each time.
# XXX: This is pretty slow. Maybe should provide an option to use a single
# apeture at the canonical wavelength when using geometric apertures?
# Note: oversampling and pad_factor need different defaults than Aperture defaults.
oversampling = self._kwargs.pop('oversampling', 1.5)
pad_factor = self._kwargs.pop('oversampling', 1.5)
aper = Aperture(diam=self.diam, lam=wave, gsparams=self._gsparams,
oversampling=oversampling, pad_factor=pad_factor, **self._kwargs)
# The aberrations were in units of wavelength for the fiducial wavelength, so we have to
# convert to units of waves for *this* wavelength.
wave_factor = self.lam / wave
# stepk and maxk also scale basically with this ratio, and they are fairly slow to
# calculate, so once we've done this once, store the results and just rescale all future
# versions with this factor.
if self._stepk is not None:
return OpticalPSF(
lam=wave, diam=self.diam,
aberrations=self.aberrations*wave_factor, scale_unit=self.scale_unit,
_force_stepk=self._stepk*wave_factor, _force_maxk=self._maxk*wave_factor,
gsparams=self.gsparams, aper=aper)
else:
ret = OpticalPSF(
lam=wave, diam=self.diam,
aberrations=self.aberrations*wave_factor, scale_unit=self.scale_unit,
gsparams=self.gsparams, aper=aper)
self._stepk = ret.stepk / wave_factor
self._maxk = ret.maxk / wave_factor
return ret
def _shoot(self, photons, rng):
if self.geometric_shooting:
# In the geometric shooting approximation, the lambda factors out, and this
# becomes the same kind of calculation we did for ChromaticAiry.
# Use the mean wavelength for the base profile.
mean_wave = np.mean(photons.wavelength)
obj = self.evaluateAtWavelength(mean_wave)
obj._shoot(photons, rng)
factor = photons.wavelength / mean_wave
photons.scaleXY(factor)
else:
# When not using geometric shooting, the following isn't exact.
# The exact method would involve doing the fourier transform for each wavelength
# in the photon list. Obviously, that's not tenable.
# So instead, we shoot with the same random seed for 3 different profiles:
# The minimum wavelength, the mean, and the maximum.
# Then interpolate between the results for each photon.
# This should (hopefully!) be good enough for most use cases if the bandpass
# isn't extremely wide and the wavelength dependence is modest over the range.
wave1 = np.min(photons.wavelength)
wave2 = np.mean(photons.wavelength)
wave3 = np.max(photons.wavelength)
prof1 = self.evaluateAtWavelength(wave1)
if wave1 == wave3:
# Interjection at this point -- if min=mean=max, then this is easy.
return prof1._shoot(photons, rng)
else:
# Otherwise we're ok dividing by wave2-wave1 and wave3-wave2 below.
assert wave2 != wave1
assert wave3 != wave2
prof2 = self.evaluateAtWavelength(wave2)
prof3 = self.evaluateAtWavelength(wave3)
# For each photon, shoot using one of these profiles according to the given
# wavelength.
# For wavelenghts with w1 < w < w2, select from prof1 or prof2 with probabilities
# P(use prof1) = (w2-w)/(w2-w1)
# P(use prof2) = (w-w1)/(w2-w1)
# Likewise when w2 < w < w3:
# P(use prof2) = (w3-w)/(w3-w2)
# P(use prof3) = (w-w2)/(w3-w2)
u = np.empty(len(photons))
UniformDeviate(rng).generate(u)
w = photons.wavelength
use_p1 = (wave1 <= w) & (w < wave2) & (u <= (wave2-w)/(wave2-wave1))
use_p2 = (wave1 <= w) & (w < wave2) & (u > (wave2-w)/(wave2-wave1))
use_p2 |= (wave2 <= w) & (w <= wave3) & (u <= (wave3-w)/(wave3-wave2))
use_p3 = (wave2 <= w) & (w <= wave3) & (u > (wave3-w)/(wave3-wave2))
assert np.all(use_p1 | use_p2 | use_p3)
assert not np.any(use_p1 & use_p2)
assert not np.any(use_p2 & use_p3)
assert not np.any(use_p1 & use_p3)
temp1 = PhotonArray(np.sum(use_p1))
temp2 = PhotonArray(np.sum(use_p2))
temp3 = PhotonArray(np.sum(use_p3))
temp1._copyFrom(photons, slice(None), use_p1, do_xy=False, do_flux=False)
temp2._copyFrom(photons, slice(None), use_p2, do_xy=False, do_flux=False)
temp3._copyFrom(photons, slice(None), use_p3, do_xy=False, do_flux=False)
prof1._shoot(temp1, rng)
prof2._shoot(temp2, rng)
prof3._shoot(temp3, rng)
temp1.scaleXY(w[use_p1]/wave1)
temp1.scaleFlux(len(temp1)/len(photons))
temp2.scaleXY(w[use_p2]/wave2)
temp2.scaleFlux(len(temp2)/len(photons))
temp3.scaleXY(w[use_p3]/wave3)
temp3.scaleFlux(len(temp3)/len(photons))
photons._copyFrom(temp1, use_p1, slice(None))
photons._copyFrom(temp2, use_p2, slice(None))
photons._copyFrom(temp3, use_p3, slice(None))
[docs]class ChromaticAiry(ChromaticObject):
"""A subclass of `ChromaticObject` meant to represent chromatic Airy profiles.
For more information about the basics of Airy profiles, please see `galsim.Airy`.
This class is a chromatic representation of Airy profiles, including the wavelength-dependent
diffraction limit. One can also get this functionality using the `ChromaticOpticalPSF` class,
but that class includes additional complications beyond a simple Airy profile, and thus has a
more complicated internal representation. For users who only want a (possibly obscured) Airy
profile, the ChromaticAiry class is likely to be a less computationally expensive and more
accurate option.
Parameters:
lam: Fiducial wavelength for which diffraction limit is initially defined, in
nanometers.
diam: Telescope diameter in meters. Either ``diam`` or ``lam_over_diam`` must be
specified.
lam_over_diam: Ratio of (fiducial wavelength) / telescope diameter in units of
``scale_unit``. Either ``diam`` or ``lam_over_diam`` must be specified.
scale_unit: Units used to define the diffraction limit and draw images.
[default: galsim.arcsec]
gsparams: An optional `GSParams` argument. See the docstring for `GSParams` for
details. [default: None]
**kwargs: Any other keyword arguments to be passed to `Airy`: either flux, or
gsparams. See `galsim.Airy` docstring for a complete description of these
options.
"""
_req_params = { 'lam' : float }
_opt_params = { 'scale_unit' : str }
_single_params = [ {'diam' : float, 'lam_over_diam' : float} ]
def __init__(self, lam, diam=None, lam_over_diam=None, scale_unit=None, gsparams=None,
**kwargs):
# First, take the basic info.
# We have to require either diam OR lam_over_diam:
if scale_unit is None:
scale_unit = arcsec
elif isinstance(scale_unit, str):
scale_unit = AngleUnit.from_name(scale_unit)
self.scale_unit = scale_unit
self._gsparams = GSParams.check(gsparams)
if ( (diam is None and lam_over_diam is None) or
(diam is not None and lam_over_diam is not None) ):
raise GalSimIncompatibleValuesError(
"Need to specify telescope diameter OR wavelength/diam ratio",
diam=diam, lam_over_diam=lam_over_diam)
if diam is not None:
self.lam_over_diam = (1.e-9*lam/diam)*radians/self.scale_unit
else:
self.lam_over_diam = float(lam_over_diam)
self.lam = float(lam)
self.kwargs = kwargs
# Define the necessary attributes for this ChromaticObject.
self.separable = False
self.interpolated = False
self.deinterpolated = self
@property
def sed(self):
return SED(1, 'nm', '1')
@property
def gsparams(self):
"""The `GSParams` for this object.
"""
return self._gsparams
[docs] @doc_inherit
def withGSParams(self, gsparams=None, **kwargs):
if gsparams == self.gsparams: return self
ret = copy.copy(self)
ret._gsparams = GSParams.check(gsparams, self.gsparams, **kwargs)
return ret
def __eq__(self, other):
return (self is other or
(isinstance(other, ChromaticAiry) and
self.lam == other.lam and
self.lam_over_diam == other.lam_over_diam and
self.scale_unit == other.scale_unit and
self.gsparams == other.gsparams and
self.kwargs == other.kwargs))
def __hash__(self):
return hash(("galsim.ChromaticAiry", self.lam, self.lam_over_diam, self.scale_unit,
self.gsparams, frozenset(self.kwargs.items())))
def __repr__(self):
s = 'galsim.ChromaticAiry(lam=%r, lam_over_diam=%r'%(self.lam, self.lam_over_diam)
if self.scale_unit != arcsec:
s += ', scale_unit=%r'%self.scale_unit
for k,v in self.kwargs.items():
s += ', %s=%r'%(k,v)
s += ', gsparams=%r'%self.gsparams
s += ')'
return s
def __str__(self):
return 'galsim.ChromaticAiry(lam=%s, lam_over_diam=%s)'%(self.lam, self.lam_over_diam)
[docs] def evaluateAtWavelength(self, wave):
"""
Method to directly instantiate a monochromatic instance of this object.
Parameters:
wave: Wavelength in nanometers.
"""
# We need to rescale the stored lam/diam by the ratio of input wavelength to stored fiducial
# wavelength.
ret = Airy(
lam_over_diam=self.lam_over_diam*(wave/self.lam), scale_unit=self.scale_unit,
gsparams=self.gsparams, **self.kwargs)
return ret
def _shoot(self, photons, rng):
# Start with the convolution at the reference wavelength
obj = Airy(lam_over_diam=self.lam_over_diam, scale_unit=self.scale_unit,
gsparams=self.gsparams, **self.kwargs)
obj._shoot(photons, rng)
# Now adjust the positions according to the wavelengths
factor = photons.wavelength / self.lam
photons.scaleXY(factor)
def _findWave(wave_list, wave):
# Helper routine to search a sorted NumPy array of wavelengths (not necessarily evenly spaced)
# to find where a particular wavelength ``wave`` would fit in, and return the index below along
# with the fraction of the way to the next entry in the array.
lower_idx = np.searchsorted(wave_list, wave)-1
# There can be edge issues, so watch out for that:
if lower_idx < 0: lower_idx = 0
if lower_idx > len(wave_list)-1: lower_idx = len(wave_list)-1
frac = (wave-wave_list[lower_idx]) / (wave_list[lower_idx+1]-wave_list[lower_idx])
return lower_idx, frac
def _linearInterp(list, frac, lower_idx):
# Helper routine for linear interpolation between values in lists (which could be lists of
# images, just not numbers, hence the need to avoid a LookupTable). Not really worth
# splitting out on its own now, but could be useful to have separate routines for the
# interpolation later on if we want to enable something other than linear interpolation.
return frac*list[lower_idx+1] + (1.-frac)*list[lower_idx]
def _remove_setup_kwargs(kwargs):
# Helper function to remove from kwargs anything that is only used for setting up image and that
# might otherwise interfere with drawImage.
kwargs.pop('dtype', None)
kwargs.pop('scale', None)
kwargs.pop('wcs', None)
kwargs.pop('nx', None)
kwargs.pop('ny', None)
kwargs.pop('bounds', None)
# Put these at the bottom to avoid circular import errors.
from . import convolve
from . import fouriersqrt
from . import transform