Phase-screen PSFs

We have available a more complicated kind of PSF model that tries to correctly model the wavefront as it passed through various “screens” such as the atmosphere or optics. It has a number of ancillary helper functions and classes associated with it to define things like the aperture and the effect of the various screens.

For PSFs drawn using real-space or Fourier methods, these utilities essentially evaluate the Fourier optics diffraction equation:

\[PSF(x,y) = \int \left| {\cal F}\left(A(u, v) e^{i \phi(u, v, x, y, t)} \right) \right|^2 dt\]

where \(x\), \(y\) are focal plane coordinates and \(u\), \(v\) are pupil plane coordinates, \(A\) is the aperture, \(\phi\) is the phase of the incident wavefront, and \(\cal F\) is the Fourier transform operator.

For method='phot', two possible strategies are available.

  1. The first strategy is to draw the PSF using Fourier methods into an InterpolatedImage, and then shoot photons from that profile. This strategy has good accuracy, but can be computationally expensive, particularly for atmospheric PSFs that need to be built up in small increments to simulate a finite exposure time.

  2. The second strategy, which can be significantly faster, especially for atmospheric PSFs, is to use the geometric optics approximation. This approximation has good accuracy for atmospheric PSFs, so we make it the default for PhaseScreenPSF. The accuracy is somewhat less good for purely optical PSFs though, so the default behavior for OpticalPSF is to use the first strategy. The geometric_shooting keyword can be used in both cases to override the default.

The main classes of note are:

Aperture

Class representing the illuminated region of pupil.

AtmosphericScreen

Class implementing phase(u, v, x, y, t) for von Karman type turbulence, with possibly evolving “non-frozen-flow” phases.

OpticalScreen

Class implementing optical aberrations using Zernike polynomial expansions in the wavefront.

UserScreen

Class implementing a user-defined phase screen.

PhaseScreenList

Python sequence type to hold multiple phase screens, for instance to simulate turbulence at different altitudes, or self-consistently model atmospheric and optical phase aberrations. A key method is PhaseScreenList.makePSF, which will take the list of phase screens, add them together linearly (Fraunhofer approximation), and evaluate the above diffraction equation to yield a PhaseScreenPSF object.

PhaseScreenPSF

A GSObject holding the evaluated PSF from a set of phase screens.

OpticalPSF

A GSObject for optical PSFs with potentially complicated pupils and Zernike aberrations.

Note

OpticalPSF is technically a kind of PhaseScreenPSF, but if you only want the optical model, you generally don’t need to bother with building any of the screens manually. The OpticalPSF class constructor will handle this for you.

SecondKick

A GSObject describing the high-k turbulence portion of an atmospheric PSF convolved by an Airy PSF. When using photon shooting with a PhaseScreenPSF, small scale (high-k) features are not well approximated by the geometric optics approximation, since this is where Fourier effects such as diffraction become important. The SecondKick class can compensate by simulating the appropriate behavior at high k values.

Atmosphere

Convenience function to quickly assemble multiple AtmosphericScreen instances into a PhaseScreenList.

class galsim.Aperture(diam, lam=None, circular_pupil=True, obscuration=0.0, nstruts=0, strut_thick=0.05, strut_angle=coord.Angle(0.0, coord.radians), oversampling=1.0, pad_factor=1.0, screen_list=None, pupil_plane_im=None, pupil_angle=coord.Angle(0.0, coord.radians), pupil_plane_scale=None, pupil_plane_size=None, gsparams=None)[source]

Class representing a telescope aperture embedded in a larger pupil plane array – for use with the PhaseScreenPSF class to create PSFs via Fourier or geometric optics.

The pupil plane array is completely specified by its size, sampling interval, and pattern of illuminated pixels. Pupil plane arrays can be specified either geometrically or using an image to indicate the illuminated pixels. In both cases, various options exist to control the pupil plane size and sampling interval.

Geometric pupil specification:

The first way to specify the details of the telescope aperture is through a series of keywords indicating the diameter, size of the central obscuration, and the nature of the struts holding up the secondary mirror (or prime focus cage, etc.). The struts are assumed to be rectangular obscurations extending from the outer edge of the pupil to the outer edge of the obscuration disk (or to the pupil center if obscuration = 0.). You can specify how many struts there are (evenly spaced in angle), how thick they are as a fraction of the pupil diameter, and what angle they start at relative to the positive y direction.

The size (in meters) and sampling interval (in meters) of the pupil plane array representing the aperture can be set directly using the the pupil_plane_size and pupil_plane_scale keywords. However, in most situations, it’s probably more convenient to let GalSim set these automatically based on the pupil geometry and the nature of the (potentially time-varying) phase aberrations from which a PSF is being derived.

The pupil plane array physical size is by default set to twice the pupil diameter producing a Nyquist sampled PSF image. While this would always be sufficient if using sinc interpolation over the PSF image for subsequent operations, GalSim by default uses the much faster (though approximate) quintic interpolant, which means that in some cases – in particular, for significantly aberrated optical PSFs without atmospheric aberrations – it may be useful to further increase the size of the pupil plane array, thereby increasing the sampling rate of the resulting PSF image. This can be done by increasing the oversampling keyword.

A caveat to the above occurs when using geometric_shooting=True to draw using photon-shooting. In this case, we only need an array just large enough to avoid clipping the pupil, which we can get by setting oversampling=0.5.

The pupil plane array physical sampling interval (which is directly related to the resulting PSF image physical size) is set by default to the same interval as would be used to avoid significant aliasing (image folding) for an obscured Airy profile with matching diameter and obscuration and for the value of folding_threshold in the optionally specified gsparams argument. If the phase aberrations are significant, however, the PSF image size computed this way may still not be sufficiently large to avoid aliasing. To further increase the pupil plane sampling rate (and hence the PSF image size), you can increase the value of the pad_factor keyword.

An additional way to set the pupil sampling interval for a particular set of phase screens (i.e., for a particular PhaseScreenList) is to provide the screens in the screen_list argument. Each screen in the list computes its own preferred sampling rate and the PhaseScreenList appropriately aggregates these. This last option also requires that a wavelength lam be specified, and is particularly helpful for creating PSFs derived from turbulent atmospheric screens.

Finally, when specifying the pupil geometrically, Aperture may choose to make a small adjustment to pupil_plane_scale in order to produce an array with a good size for FFTs. If your application depends on knowing the size and scale used with the Fourier optics framework, you can obtain these from the aper.pupil_plane_size and aper.pupil_plane_scale attributes.

Pupil image specification:

The second way to specify the pupil plane configuration is by passing in an image of it. This can be useful, for example, if the struts are not evenly spaced or are not radially directed, as is assumed by the simple model for struts described above. In this case, an exception is raised if keywords related to struts are also given. On the other hand, the obscuration keyword is still used to ensure that the PSF images are not aliased, though it is ignored during the actual construction of the pupil plane illumination pattern. Note that for complicated pupil configurations, it may be desireable to increase pad_factor for more fidelity at the expense of slower running time. Finally, the pupil_plane_im that is passed in can be rotated during internal calculations by specifying a pupil_angle keyword.

If you choose to pass in a pupil plane image, it must be a square array in which the image of the pupil is centered. The areas that are illuminated should have some value >0, and the other areas should have a value of precisely zero. Based on what the Aperture class determines is a good PSF sampling interval, the image of the pupil plane that is passed in might be zero-padded during internal calculations. (The pupil plane array size and scale values can be accessed via the aper.pupil_plane_size and aper.pupil_plane_scale attributes.) The pixel scale of the pupil plane can be specified in one of three ways. In descending order of priority, these are:

  1. The pupil_plane_scale keyword argument (units are meters).

  2. The pupil_plane_im.scale attribute (units are meters).

  3. If (1) and (2) are both None, then the scale will be inferred by assuming that the illuminated pixel farthest from the image center is at a physical distance of self.diam/2.

The pupil_plane_size and lam keywords are both ignored when constructing an Aperture from an image.

Parameters:
  • diam – Aperture diameter in meters.

  • lam – Wavelength in nanometers. [default: None]

  • circular_pupil – Adopt a circular pupil? [default: True]

  • obscuration – Linear dimension of central obscuration as fraction of aperture linear dimension. [0., 1.). [default: 0.0]

  • nstruts – Number of radial support struts to add to the central obscuration. [default: 0]

  • strut_thick – Thickness of support struts as a fraction of aperture diameter. [default: 0.05]

  • strut_angleAngle made between the vertical and the strut starting closest to it, defined to be positive in the counter-clockwise direction; must be an Angle instance. [default: 0. * galsim.degrees]

  • oversampling – Optional oversampling factor in the image plane for the PSF eventually constructed using this Aperture. Setting oversampling < 1 will produce aliasing in the PSF (not good). [default: 1.0]

  • pad_factor – Additional multiple by which to extend the PSF image to avoid folding. [default: 1.0]

  • screen_list – An optional PhaseScreenList object. If present, then get a good pupil sampling interval using this object. [default: None]

  • pupil_plane_im – The GalSim.Image, NumPy array, or name of file containing the pupil plane image, to be used instead of generating one based on the obscuration and strut parameters. [default: None]

  • pupil_angle – If pupil_plane_im is not None, rotation angle for the pupil plane (positive in the counter-clockwise direction). Must be an Angle instance. [default: 0. * galsim.degrees]

  • pupil_plane_scale – Sampling interval in meters to use for the pupil plane array. In most cases, it’s a good idea to leave this as None, in which case GalSim will attempt to find a good value automatically. The exception is when specifying the pupil arrangement via an image, in which case this keyword can be used to indicate the sampling of that image. See also pad_factor for adjusting the pupil sampling scale. [default: None]

  • pupil_plane_size – Size in meters to use for the pupil plane array. In most cases, it’s a good idea to leave this as None, in which case GalSim will attempt to find a good value automatically. See also oversampling for adjusting the pupil size. [default: None]

  • gsparams – An optional GSParams argument. [default: None]

property diam

Aperture diameter in meters

property gsparams

The GSParams of this object.

property illuminated

A boolean array indicating which positions in the pupil plane are exposed to the sky.

property npix

The number of pixels in each direction of the pupil-plane image.

property obscuration

Fraction linear obscuration of pupil.

property pupil_plane_scale

The scale_size of the pupil-plane image.

property pupil_plane_size

The size of the pupil-plane image.

samplePupil(photons, rng)[source]

Set the pupil_u and pupil_v values in the PhotonArray by sampling the current aperture.

property u

Pupil horizontal coordinate array in meters.

property v

Pupil vertical coordinate array in meters.

withGSParams(gsparams=None, **kwargs)[source]

Create a version of the current aperture with the given gsparams

class galsim.AtmosphericScreen(screen_size, screen_scale=None, altitude=0.0, r0_500=0.2, L0=25.0, vx=0.0, vy=0.0, alpha=1.0, time_step=None, rng=None, suppress_warning=False, mp_context=None)[source]

An atmospheric phase screen that can drift in the wind and evolves (“boils”) over time. The initial phases and fractional phase updates are drawn from a von Karman power spectrum, which is defined by a Fried parameter that effectively sets the amplitude of the turbulence, and an outer scale beyond which the turbulence power flattens.

AtmosphericScreen delays the actual instantiation of the phase screen array in memory until it is used for either drawing a PSF or querying the wavefront or wavefront gradient. This is to facilitate automatic truncation of the screen power spectrum depending on the use case. For example, when drawing a PhaseScreenPSF using Fourier methods, the entire power spectrum should generally be used. On the other hand, when drawing using photon-shooting and the geometric approximation, it’s better to truncate the high-k modes of the power spectrum here so that they can be handled instead by a SecondKick object (which also happens automatically; see the PhaseScreenPSF docstring). (See Peterson et al. 2015 for more details about the second kick). Querying the wavefront or wavefront gradient will instantiate the screen using the full power spectrum.

This class will normally attempt to sanity check that the screen has been appropriately instantiated depending on the use case, i.e., depending on whether it’s being used to draw with Fourier optics or geometric optics. If you want to turn this warning off, however, you can use the suppress_warning keyword argument.

If you wish to override the automatic truncation determination, then you can directly instantiate the phase screen array using the AtmosphericScreen.instantiate method.

Note that once a screen has been instantiated with a particular set of truncation parameters, it cannot be re-instantiated with another set of parameters.

Shared memory:

Instantiated AtmosphericScreen objects can consume a significant amount of memory. For example, an atmosphere with 6 screens, each extending 819.2 m and with resolution of 10 cm will consume 3 GB of memory. In contexts where both a realistic atmospheric PSF and high throughput via multiprocessing are required, allocating this 3 GB of memory once in a shared memory space accessible to each subprocess (as opposed to once per subprocess) is highly desireable. We provide a few functions here to enable such usage:

  • The mp_context keyword argument to AtmosphericScreen. This is used to indicate which multiprocessing process launching context will be used. This is important for setting up the shared memory correctly.

  • The initWorker and initWorkerArgs functions. These should be used in a call to multiprocessing.Pool to correctly inform the worker process where to find AtmosphericScreen shared memory.

A template example might look something like:

import galsim
import multiprocessing as mp

def work(i, atm):
    args, moreArgs = fn(i)
    psf = atm.makePSF(*args)
    return psf.drawImage(*moreArgs)

ctx = mp.get_context("spawn")  # "spawn" is generally the safest context available

atm = galsim.Atmosphere(..., mp_context=ctx)  # internally calls AtmosphericScreen ctor
nProc = 4  # Note, can set this to None to get a good default

# Note: Especially if you are using "fork" context, then you want to make sure to run
#       your Pool in a single_threaded context.  Even if not, it's probably a good idea
#       so each process isn't spawning lots of OpenMP (or other) threads.
with galsim.utilities.single_threaded():
    with ctx.Pool(
        nProc,
        initializer=galsim.phase_screens.initWorker,
        initargs=galsim.phase_screens.initWorkerArgs()
    ) as pool:
        results = []
        # First submit
        for i in range(10):
            results.append(pool.apply_async(work, (i, atm)))
        # Then wait to finish
        for r in results:
            r.wait()
# Turn future objects into actual returned images.
results = [r.get() for r in results]

It is also possible to manually instantiate each of the AtmosphericScreen objects in a PhaseScreenList in parallel using a process pool. This requires knowing what k-scale to truncate the screen at:

atm = galsim.Atmosphere(..., mp_context=ctx)
with galsim.utilities.single_threaded():
    with ctx.Pool(
        nProc,
        initializer=galsim.phase_screens.initWorker,
        initargs=galsim.phase_screens.initWorkerArgs()
    ) as pool:
        dummyPSF = atm.makePSF(...)
        kmax = dummyPSF.screen_kmax
        atm.instantiate(pool=pool, kmax=kmax)

Finally, the above multiprocessing shared memory tricks are only currently supported for non-time-evolving screens (alpha=1).

Pickling:

The shared memory portion of an atmospheric screen is not included by default in the pickle of an AtmosphericScreen instance. This means that while it is possible to pickle/unpickle an AtmosphericScreen as normal within a single launch of a potentially-multiprocess program, (as long as the shared memory is persistent), a different path is needed to say, create an AtmosphericScreen pickle, quit python, restart python, and unpickle the screen. The same holds for any object that wraps an AtmosphericScreen as an attribute as well, which could include, e.g., PhaseScreenList or PhaseScreenPSF.

To get around this limitation, the context manager galsim.utilities.pickle_shared can be used. For example:

screen = galsim.AtmosphericScreen(...)
with galsim.utilities.pickle_shared():
    pickle.dump(screen, open('myScreen.pkl', 'wb'))

will pickle both the screen object and any required shared memory used in its definition. Unpickling then proceeds exactly as normal:

screen = pickle.load(open('myScreen.pkl', 'rb'))
Parameters:
  • screen_size – Physical extent of square phase screen in meters. This should be large enough to accommodate the desired field-of-view of the telescope as well as the meta-pupil defined by the wind speed and exposure time. Note that the screen will have periodic boundary conditions, so while the code will still run with a small screen, this may introduce artifacts into PSFs or PSF correlation functions. Also note that screen_size may be tweaked by the initializer to ensure screen_size is a multiple of screen_scale.

  • screen_scale – Physical pixel scale of phase screen in meters. An order unity multiple of the Fried parameter is usually sufficiently small, but users should test the effects of varying this parameter to ensure robust results. [default: r0_500]

  • altitude – Altitude of phase screen in km. This is with respect to the telescope, not sea-level. [default: 0.0]

  • r0_500 – Fried parameter setting the amplitude of turbulence; contributes to “size” of the resulting atmospheric PSF. Specified at wavelength 500 nm, in units of meters. [default: 0.2]

  • L0 – Outer scale in meters. The turbulence power spectrum will smoothly approach a constant at scales larger than L0. Set to None or np.inf for a power spectrum without an outer scale. [default: 25.0]

  • vx – x-component wind velocity in meters/second. [default: 0.]

  • vy – y-component wind velocity in meters/second. [default: 0.]

  • alpha – Square root of fraction of phase that is “remembered” between time_steps (i.e., alpha**2 is the fraction remembered). The fraction sqrt(1-alpha**2) is then the amount of turbulence freshly generated in each step. Setting alpha=1.0 results in a frozen-flow atmosphere. Note that computing PSFs from frozen-flow atmospheres may be significantly faster than computing PSFs with non-frozen-flow atmospheres. If alpha != 1.0, then it is required that a time_step is also specified. [default: 1.0]

  • time_step – Time interval between phase boiling updates. Note that this is distinct from the time interval used to integrate the PSF over time, which is set by the time_step keyword argument to PhaseScreenPSF or PhaseScreenList.makePSF. If time_step is not None, then it is required that alpha is set to something other than 1.0. [default: None]

  • rng – Random number generator as a BaseDeviate. If None, then use the clock time or system entropy to seed a new generator. [default: None]

  • suppress_warning – Turn off instantiation sanity checking. (See above) [default: False]

  • enable (mp_context GalSim uses shared memory for phase screen allocation to better) – multiprocessing. Use this keyword to set the launch context for multiprocessing. Usually it will be sufficient to leave this at its default. [default: None]

Relevant SPIE paper: “Remembrance of phases past: An autoregressive method for generating realistic atmospheres in simulations” Srikar Srinath, Univ. of California, Santa Cruz; Lisa A. Poyneer, Lawrence Livermore National Lab.; Alexander R. Rudy, UCSC; S. Mark Ammons, LLNL Published in Proceedings Volume 9148: Adaptive Optics Systems IV September 2014

property altitude

The altitude of the screen in km.

instantiate(kmin=0.0, kmax=inf, check=None)[source]
Parameters:
  • kmin – Minimum k-mode to include when generating phase screens. Generally this will only be used when testing the geometric approximation for atmospheric PSFs. [default: 0]

  • kmax – Maximum k-mode to include when generating phase screens. This may be used in conjunction with SecondKick to complete the geometric approximation for atmospheric PSFs. [default: np.inf]

  • check – Sanity check indicator. If equal to ‘FFT’, then check that phase screen Fourier modes are not being truncated, which is appropriate for full Fourier optics. If equal to ‘phot’, then check that phase screen Fourier modes are being truncated, which is appropriate for the geometric optics approximation. If None, then don’t perform a check. Also, don’t perform a check if self.suppress_warning is True.

property kmax

The maximum k value being used.

property kmin

The minimum k value being used.

wavefront(u, v, t=None, theta=(coord.Angle(0.0, coord.radians), coord.Angle(0.0, coord.radians)))[source]

Compute wavefront due to atmospheric phase screen.

Wavefront here indicates the distance by which the physical wavefront lags or leads the ideal plane wave.

Parameters:
  • u – Horizontal pupil coordinate (in meters) at which to evaluate wavefront. Can be a scalar or an iterable. The shapes of u and v must match.

  • v – Vertical pupil coordinate (in meters) at which to evaluate wavefront. Can be a scalar or an iterable. The shapes of u and v must match.

  • t – Times (in seconds) at which to evaluate wavefront. Can be None, a scalar or an iterable. If None, then the internal time of the phase screens will be used for all u, v. If scalar, then the size will be broadcast up to match that of u and v. If iterable, then the shape must match the shapes of u and v. [default: None]

  • theta – Field angle at which to evaluate wavefront, as a 2-tuple of galsim.Angle instances. [default: (0.0*galsim.arcmin, 0.0*galsim.arcmin)] Only a single theta is permitted.

Returns:

Array of wavefront lag or lead in nanometers.

wavefront_gradient(u, v, t=None, theta=(coord.Angle(0.0, coord.radians), coord.Angle(0.0, coord.radians)))[source]

Compute gradient of wavefront due to atmospheric phase screen.

Parameters:
  • u – Horizontal pupil coordinate (in meters) at which to evaluate wavefront. Can be a scalar or an iterable. The shapes of u and v must match.

  • v – Vertical pupil coordinate (in meters) at which to evaluate wavefront. Can be a scalar or an iterable. The shapes of u and v must match.

  • t – Times (in seconds) at which to evaluate wavefront gradient. Can be None, a scalar or an iterable. If None, then the internal time of the phase screens will be used for all u, v. If scalar, then the size will be broadcast up to match that of u and v. If iterable, then the shape must match the shapes of u and v. [default: None]

  • theta – Field angle at which to evaluate wavefront, as a 2-tuple of galsim.Angle instances. [default: (0.0*galsim.arcmin, 0.0*galsim.arcmin)] Only a single theta is permitted.

Returns:

Arrays dWdu and dWdv of wavefront lag or lead gradient in nm/m.

class galsim.OpticalScreen(diam, tip=0.0, tilt=0.0, defocus=0.0, astig1=0.0, astig2=0.0, coma1=0.0, coma2=0.0, trefoil1=0.0, trefoil2=0.0, spher=0.0, aberrations=None, annular_zernike=False, obscuration=0.0, lam_0=500.0)[source]

Class to describe optical aberrations in terms of Zernike polynomial coefficients.

Input aberration coefficients are assumed to be supplied in units of wavelength, and correspond to the Zernike polynomials in the Noll convention defined in Noll, J. Opt. Soc. Am. 66, 207-211(1976). For a brief summary of the polynomials, refer to http://en.wikipedia.org/wiki/Zernike_polynomials#Zernike_polynomials.

Parameters:
  • diam – Diameter of pupil in meters.

  • tip – Tip aberration in units of reference wavelength. [default: 0]

  • tilt – Tilt aberration in units of reference wavelength. [default: 0]

  • defocus – Defocus in units of reference wavelength. [default: 0]

  • astig1 – Astigmatism (like e2) in units of reference wavelength. [default: 0]

  • astig2 – Astigmatism (like e1) in units of reference wavelength. [default: 0]

  • coma1 – Coma along y in units of reference wavelength. [default: 0]

  • coma2 – Coma along x in units of reference wavelength. [default: 0]

  • trefoil1 – Trefoil (one of the arrows along y) in units of reference wavelength. [default: 0]

  • trefoil2 – Trefoil (one of the arrows along x) in units of reference wavelength. [default: 0]

  • spher – Spherical aberration in units of reference wavelength. [default: 0]

  • aberrations – Optional keyword, to pass in a list, tuple, or NumPy array of aberrations in units of reference wavelength (ordered according to the Noll convention), rather than passing in individual values for each individual aberration. Note that aberrations[1] is piston (and not aberrations[0], which is unused.) This list can be arbitrarily long to handle Zernike polynomial aberrations of arbitrary order.

  • annular_zernike – Boolean indicating that aberrations specify the amplitudes of annular Zernike polynomials instead of circular Zernike polynomials. [default: False]

  • obscuration – Linear dimension of central obscuration as fraction of aperture linear dimension. [0., 1.). Note it is the user’s responsibility to ensure consistency of OpticalScreen obscuration and Aperture obscuration. [default: 0.0]

  • lam_0 – Reference wavelength in nanometers at which Zernike aberrations are being specified. [default: 500]

wavefront(u, v, t=None, theta=None)[source]

Compute wavefront due to optical phase screen.

Wavefront here indicates the distance by which the physical wavefront lags or leads the ideal plane wave.

Parameters:
  • u – Horizontal pupil coordinate (in meters) at which to evaluate wavefront. Can be a scalar or an iterable. The shapes of u and v must match.

  • v – Vertical pupil coordinate (in meters) at which to evaluate wavefront. Can be a scalar or an iterable. The shapes of u and v must match.

  • t – Ignored for OpticalScreen.

  • theta – Ignored for OpticalScreen.

Returns:

Array of wavefront lag or lead in nanometers.

wavefront_gradient(u, v, t=None, theta=None)[source]

Compute gradient of wavefront due to optical phase screen.

Parameters:
  • u – Horizontal pupil coordinate (in meters) at which to evaluate wavefront. Can be a scalar or an iterable. The shapes of u and v must match.

  • v – Vertical pupil coordinate (in meters) at which to evaluate wavefront. Can be a scalar or an iterable. The shapes of u and v must match.

  • t – Ignored for OpticalScreen.

  • theta – Ignored for OpticalScreen.

Returns:

Arrays dWdu and dWdv of wavefront lag or lead gradient in nm/m.

class galsim.UserScreen(table, diam=None, obscuration=0.0)[source]

Create a (static) user-defined phase screen.

Parameters:
  • table – LookupTable2D instance representing the wavefront as a function on the entrance pupil. Units are (meters, meters) -> nanometers.

  • diam – Diameter of entrance pupil in meters. If None, then use the length of the larger side of the LookupTable2D rectangle in table. This keyword is only used to compute a value for stepk, and thus has no effect on the wavefront() or wavefront_gradient() methods.

  • obscuration – Optional fractional circular obscuration of pupil. Like diam, only used for computing a value for stepk.

wavefront(u, v, t=None, theta=None)[source]

Evaluate wavefront from lookup table.

Wavefront here indicates the distance by which the physical wavefront lags or leads the ideal plane wave.

Parameters:
  • u – Horizontal pupil coordinate (in meters) at which to evaluate wavefront. Can be a scalar or an iterable. The shapes of u and v must match.

  • v – Vertical pupil coordinate (in meters) at which to evaluate wavefront. Can be a scalar or an iterable. The shapes of u and v must match.

  • t – Ignored for UserScreen.

  • theta – Ignored for UserScreen.

Returns:

Array of wavefront lag or lead in nanometers.

wavefront_gradient(u, v, t=None, theta=None)[source]

Evaluate gradient of wavefront from lookup table.

Parameters:
  • u – Horizontal pupil coordinate (in meters) at which to evaluate wavefront. Can be a scalar or an iterable. The shapes of u and v must match.

  • v – Vertical pupil coordinate (in meters) at which to evaluate wavefront. Can be a scalar or an iterable. The shapes of u and v must match.

  • t – Ignored for UserScreen.

  • theta – Ignored for UserScreen.

Returns:

Arrays dWdu and dWdv of wavefront lag or lead gradient in nm/m.

class galsim.PhaseScreenList(*layers)[source]

List of phase screens that can be turned into a PSF. Screens can be either atmospheric layers or optical phase screens. Generally, one would assemble a PhaseScreenList object using the function Atmosphere. Layers can be added, removed, appended, etc. just like items can be manipulated in a python list. For example:

# Create an atmosphere with three layers.
>>> screens = galsim.PhaseScreenList([galsim.AtmosphericScreen(...),
                                      galsim.AtmosphericScreen(...),
                                      galsim.AtmosphericScreen(...)])
# Add another layer
>>> screens.append(galsim.AtmosphericScreen(...))
# Remove the second layer
>>> del screens[1]
# Switch the first and second layer.  Silly, but works...
>>> screens[0], screens[1] = screens[1], screens[0]
Parameters:

layers – Sequence of phase screens.

instantiate(pool=None, _bar=None, **kwargs)[source]

Instantiate the screens in this PhaseScreenList.

Parameters:
  • pool – A multiprocessing.Pool object to use to instantiate screens in parallel.

  • **kwargs – Keyword arguments to forward to screen.instantiate().

makePSF(lam, **kwargs)[source]

Create a PSF from the current PhaseScreenList.

Parameters:
  • lam – Wavelength in nanometers at which to compute PSF.

  • t0 – Time at which to start exposure in seconds. [default: 0.0]

  • exptime – Time in seconds over which to accumulate evolving instantaneous PSF. [default: 0.0]

  • time_step – Time interval in seconds with which to sample phase screens when drawing using real-space or Fourier methods, or when using photon-shooting without the geometric optics approximation. Note that the default value of 0.025 is fairly arbitrary. For careful studies, we recommend checking that results are stable when decreasing time_step. Also note that when drawing using photon-shooting with the geometric optics approximation this keyword is ignored, as the phase screen can be sampled continuously in this case instead of at discrete intervals. [default: 0.025]

  • flux – Flux of output PSF. [default: 1.0]

  • theta – Field angle of PSF as a 2-tuple of Angle instances. [default: (0.0*galsim.arcmin, 0.0*galsim.arcmin)]

  • interpolant – Either an Interpolant instance or a string indicating which interpolant should be used. Options are ‘nearest’, ‘sinc’, ‘linear’, ‘cubic’, ‘quintic’, or ‘lanczosN’ where N should be the integer order to use. [default: galsim.Quintic()]

  • scale_unit – Units to use for the sky coordinates of the output profile. [default: galsim.arcsec]

  • ii_pad_factor – Zero-padding factor by which to extend the image of the PSF when creating the InterpolatedImage. See the InterpolatedImage docstring for more details. [default: 1.5]

  • suppress_warning – If pad_factor is too small, the code will emit a warning telling you its best guess about how high you might want to raise it. However, you can suppress this warning by using suppress_warning=True. [default: False]

  • 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: True]

  • aperAperture to use to compute PSF(s). [default: None]

  • second_kick – An optional second kick to also convolve by when using geometric photon-shooting. (This can technically be any GSObject, though usually it should probably be a SecondKick object). If None, then a good second kick will be chosen automatically based on screen_list. If False, then a second kick won’t be applied. [default: None]

  • kcrit – Critical Fourier scale (in units of 1/r0) at which to separate low-k and high-k turbulence. The default value was chosen based on comparisons between Fourier optics and geometric optics with a second kick correction. While most values of kcrit smaller than the default produce similar results, we caution the user to compare the affected geometric PSFs against Fourier optics PSFs carefully before changing this value. [default: 0.2]

  • fft_sign – The sign (+/-) to use in the exponent of the Fourier kernel when evaluating the Fourier optics PSF. As of version 2.3, GalSim uses a plus sign by default, which we believe to be consistent with, for example, how Zemax computes a Fourier optics PSF on DECam. Before version 2.3, the default was a negative sign. Input should be either the string ‘+’ or the string ‘-’. [default: ‘+’]

  • gsparams – An optional GSParams argument. [default: None]

The following are optional keywords to use to setup the aperture if aper is not provided.

Parameters:
  • diam – Aperture diameter in meters.

  • circular_pupil – Adopt a circular pupil? [default: True]

  • obscuration – Linear dimension of central obscuration as fraction of aperture linear dimension. [0., 1.). [default: 0.0]

  • nstruts – Number of radial support struts to add to the central obscuration. [default: 0]

  • strut_thick – Thickness of support struts as a fraction of aperture diameter. [default: 0.05]

  • strut_angleAngle made between the vertical and the strut starting closest to it, defined to be positive in the counter-clockwise direction; must be an Angle instance. [default: 0. * galsim.degrees]

  • oversampling – Optional oversampling factor in the image plane for the PSF eventually constructed using this Aperture. Setting oversampling < 1 will produce aliasing in the PSF (not good). [default: 1.0]

  • pad_factor – Additional multiple by which to extend the PSF image to avoid folding. [default: 1.0]

  • pupil_plane_im – The GalSim.Image, NumPy array, or name of file containing the pupil plane image, to be used instead of generating one based on the obscuration and strut parameters. [default: None]

  • pupil_angle – If pupil_plane_im is not None, rotation angle for the pupil plane (positive in the counter-clockwise direction). Must be an Angle instance. [default: 0. * galsim.degrees]

  • pupil_plane_scale – Sampling interval in meters to use for the pupil plane array. In most cases, it’s a good idea to leave this as None, in which case GalSim will attempt to find a good value automatically. The exception is when specifying the pupil arrangement via an image, in which case this keyword can be used to indicate the sampling of that image. See also pad_factor for adjusting the pupil sampling scale. [default: None]

  • pupil_plane_size – Size in meters to use for the pupil plane array. In most cases, it’s a good idea to leave this as None, in which case GalSim will attempt to find a good value automatically. See also oversampling for adjusting the pupil size. [default: None]

wavefront(u, v, t, theta=(coord.Angle(0.0, coord.radians), coord.Angle(0.0, coord.radians)))[source]

Compute cumulative wavefront due to all phase screens in PhaseScreenList.

Wavefront here indicates the distance by which the physical wavefront lags or leads the ideal plane wave (pre-optics) or spherical wave (post-optics).

Parameters:
  • u – Horizontal pupil coordinate (in meters) at which to evaluate wavefront. Can be a scalar or an iterable. The shapes of u and v must match.

  • v – Vertical pupil coordinate (in meters) at which to evaluate wavefront. Can be a scalar or an iterable. The shapes of u and v must match.

  • t – Times (in seconds) at which to evaluate wavefront. Can be None, a scalar or an iterable. If None, then the internal time of the phase screens will be used for all u, v. If scalar, then the size will be broadcast up to match that of u and v. If iterable, then the shape must match the shapes of u and v.

  • theta – Field angle at which to evaluate wavefront, as a 2-tuple of galsim.Angle instances. [default: (0.0*galsim.arcmin, 0.0*galsim.arcmin)] Only a single theta is permitted.

Returns:

Array of wavefront lag or lead in nanometers.

wavefront_gradient(u, v, t, theta=(coord.Angle(0.0, coord.radians), coord.Angle(0.0, coord.radians)))[source]

Compute cumulative wavefront gradient due to all phase screens in PhaseScreenList.

Parameters:
  • u – Horizontal pupil coordinate (in meters) at which to evaluate wavefront. Can be a scalar or an iterable. The shapes of u and v must match.

  • v – Vertical pupil coordinate (in meters) at which to evaluate wavefront. Can be a scalar or an iterable. The shapes of u and v must match.

  • t – Times (in seconds) at which to evaluate wavefront gradient. Can be None, a scalar or an iterable. If None, then the internal time of the phase screens will be used for all u, v. If scalar, then the size will be broadcast up to match that of u and v. If iterable, then the shape must match the shapes of u and v.

  • theta – Field angle at which to evaluate wavefront, as a 2-tuple of galsim.Angle instances. [default: (0.0*galsim.arcmin, 0.0*galsim.arcmin)] Only a single theta is permitted.

Returns:

Arrays dWdu and dWdv of wavefront lag or lead gradient in nm/m.

class galsim.PhaseScreenPSF(screen_list, lam, t0=0.0, exptime=0.0, time_step=0.025, flux=1.0, theta=(coord.Angle(0.0, coord.radians), coord.Angle(0.0, coord.radians)), interpolant=None, scale_unit=coord.arcsec, ii_pad_factor=None, suppress_warning=False, geometric_shooting=True, aper=None, second_kick=None, kcrit=0.2, fft_sign='+', gsparams=None, _force_stepk=0.0, _force_maxk=0.0, _bar=None, **kwargs)[source]

Bases: GSObject

A PSF surface brightness profile constructed by integrating over time the instantaneous PSF derived from a set of phase screens and an aperture.

There are two equivalent ways to construct a PhaseScreenPSF given a PhaseScreenList:

>>> psf = screen_list.makePSF(...)
>>> psf = PhaseScreenPSF(screen_list, ...)

Computing a PSF from a phase screen also requires an Aperture be specified. This can be done either directly via the aper keyword, or by setting a number of keywords that will be passed to the Aperture constructor. The aper keyword always takes precedence.

There are effectively three ways to draw a PhaseScreenPSF (or GSObject that includes a PhaseScreenPSF):

  1. Fourier optics

    This is the default, and is performed for all drawImage methods except method=’phot’. This is generally the most accurate option. For a PhaseScreenList that includes an AtmosphericScreen, however, this can be prohibitively slow. For OpticalPSF, though, this can sometimes be a good option.

  2. Photon-shooting from an image produced using Fourier optics.

    This is done if geometric_shooting=False when creating the PhaseScreenPSF, and method=’phot’ when calling drawImage. This actually performs the same calculations as the Fourier optics option above, but then proceeds by shooting photons from that result. This can sometimes be a good option for OpticalPSFs, especially if the same OpticalPSF can be reused for may objects, since the Fourier part of the process would only be performed once in this case.

  3. Photon-shooting using the “geometric approximation”.

    This is done if geometric_shooting=True when creating the PhaseScreenPSF, and method=’phot’ when calling drawImage. In this case, a completely different algorithm is used make an image. Photons are uniformly generated in the Aperture pupil, and then the phase gradient at that location is used to deflect each photon in the image plane. This method, which corresponds to geometric optics, is broadly accurate for phase screens that vary slowly across the aperture, and is usually several orders of magnitude or more faster than Fourier optics (depending on the flux of the object, of course, but much faster even for rather bright flux levels).

    One short-coming of this method is that it neglects interference effects, i.e. diffraction. For PhaseScreenList that include at least one AtmosphericScreen, a correction, dubbed the “second kick”, will automatically be applied to handle both the quickly varying modes of the screens and the diffraction pattern of the Aperture. For PhaseScreenLists without an AtmosphericScreen, the correction is simply an Airy function. Note that this correction can be overridden using the second_kick keyword argument, and also tuned to some extent using the kcrit keyword argument.

Note also that calling drawImage on a PhaseScreenPSF that uses a PhaseScreenList with any uninstantiated AtmosphericScreen will perform that instantiation, and that the details of the instantiation depend on the drawing method used, and also the kcrit keyword argument to PhaseScreenPSF. See the AtmosphericScreen docstring for more details.

Parameters:
  • screen_listPhaseScreenList object from which to create PSF.

  • lam – Wavelength in nanometers at which to compute PSF.

  • t0 – Time at which to start exposure in seconds. [default: 0.0]

  • exptime – Time in seconds over which to accumulate evolving instantaneous PSF. [default: 0.0]

  • time_step – Time interval in seconds with which to sample phase screens when drawing using real-space or Fourier methods, or when using photon-shooting without the geometric optics approximation. Note that the default value of 0.025 is fairly arbitrary. For careful studies, we recommend checking that results are stable when decreasing time_step. Also note that when drawing using photon-shooting with the geometric optics approximation this keyword is ignored, as the phase screen can be sampled continuously in this case instead of at discrete intervals. [default: 0.025]

  • flux – Flux of output PSF [default: 1.0]

  • theta – Field angle of PSF as a 2-tuple of Angle instances. [default: (0.0*galsim.arcmin, 0.0*galsim.arcmin)]

  • interpolant – Either an Interpolant instance or a string indicating which interpolant should be used. Options are ‘nearest’, ‘sinc’, ‘linear’, ‘cubic’, ‘quintic’, or ‘lanczosN’ where N should be the integer order to use. [default: galsim.Quintic()]

  • scale_unit – Units to use for the sky coordinates of the output profile. [default: galsim.arcsec]

  • ii_pad_factor – Zero-padding factor by which to extend the image of the PSF when creating the InterpolatedImage. See the InterpolatedImage docstring for more details. [default: 1.5]

  • suppress_warning – If pad_factor is too small, the code will emit a warning telling you its best guess about how high you might want to raise it. However, you can suppress this warning by using suppress_warning=True. [default: False]

  • 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: True]

  • aperAperture to use to compute PSF(s). [default: None]

  • second_kick – An optional second kick to also convolve by when using geometric photon-shooting. (This can technically be any GSObject, though usually it should probably be a SecondKick object). If None, then a good second kick will be chosen automatically based on screen_list. If False, then a second kick won’t be applied. [default: None]

  • kcrit – Critical Fourier scale (in units of 1/r0) at which to separate low-k and high-k turbulence. The default value was chosen based on comparisons between Fourier optics and geometric optics with a second kick correction. While most values of kcrit smaller than the default produce similar results, we caution the user to compare the affected geometric PSFs against Fourier optics PSFs carefully before changing this value. [default: 0.2]

  • fft_sign – The sign (+/-) to use in the exponent of the Fourier kernel when evaluating the Fourier optics PSF. As of version 2.3, GalSim uses a plus sign by default, which we believe to be consistent with, for example, how Zemax computes a Fourier optics PSF on DECam. Before version 2.3, the default was a negative sign. Input should be either the string ‘+’ or the string ‘-’. [default: ‘+’]

  • gsparams – An optional GSParams argument. [default: None]

The following are optional keywords to use to setup the aperture if aper is not provided:

Parameters:
  • diam – Aperture diameter in meters. [default: None]

  • circular_pupil – Adopt a circular pupil? [default: True]

  • obscuration – Linear dimension of central obscuration as fraction of aperture linear dimension. [0., 1.). [default: 0.0]

  • nstruts – Number of radial support struts to add to the central obscuration. [default: 0]

  • strut_thick – Thickness of support struts as a fraction of aperture diameter. [default: 0.05]

  • strut_angleAngle made between the vertical and the strut starting closest to it, defined to be positive in the counter-clockwise direction; must be an Angle instance. [default: 0. * galsim.degrees]

  • oversampling – Optional oversampling factor in the image plane for the PSF eventually constructed using this Aperture. Setting oversampling < 1 will produce aliasing in the PSF (not good). [default: 1.0]

  • pad_factor – Additional multiple by which to extend the PSF image to avoid folding. [default: 1.0]

  • pupil_plane_im – The GalSim.Image, NumPy array, or name of file containing the pupil plane image, to be used instead of generating one based on the obscuration and strut parameters. [default: None]

  • pupil_angle – If pupil_plane_im is not None, rotation angle for the pupil plane (positive in the counter-clockwise direction). Must be an Angle instance. [default: 0. * galsim.degrees]

  • pupil_plane_scale – Sampling interval in meters to use for the pupil plane array. In most cases, it’s a good idea to leave this as None, in which case GalSim will attempt to find a good value automatically. The exception is when specifying the pupil arrangement via an image, in which case this keyword can be used to indicate the sampling of that image. See also pad_factor for adjusting the pupil sampling scale. [default: None]

  • pupil_plane_size – Size in meters to use for the pupil plane array. In most cases, it’s a good idea to leave this as None, in which case GalSim will attempt to find a good value automatically. See also oversampling for adjusting the pupil size. [default: None]

property fft_sign

The sign (+/-) to use in the exponent of the Fourier kernel when evaluating the Fourier optics PSF.

property flux

The flux of the profile.

property kcrit

The critical Fourier scale being used for this object.

property screen_list

The PhaseScreenList being used for this object.

withFlux(flux)[source]

Create a version of the current object with a different flux.

This function is equivalent to obj.withScaledFlux(flux / obj.flux).

It creates a new object that has the same profile as the original, but with the surface brightness at every location rescaled such that the total flux will be the given value. Note that if flux is an SED, the return value will be a ChromaticObject with specified SED.

Parameters:

flux – The new flux for the object.

Returns:

the object with the new flux

withGSParams(gsparams=None, **kwargs)[source]

Create a version of the current object with the given GSParams.

You may either provide a GSParams instance:

>>> gsparams = galsim.GSParams(folding_threshold=1.e-4, maxk_threshold=1.e-4)
>>> obj = obj.withGSParams(gsparams)

or one or more named parameters as keyword arguments:

>>> obj = obj.withGSParams(folding_threshold=1.e-4, maxk_threshold=1.e-4)

Note

The latter style will leave all non-named parameters at their current values. It only updates the named parameters to the given values.

class galsim.SecondKick(lam, r0, diam, obscuration=0, kcrit=0.2, flux=1, scale_unit=coord.arcsec, gsparams=None)[source]

Bases: GSObject

Class describing the expectation value of the high-k turbulence portion of an atmospheric PSF convolved by an Airy PSF.

The power spectrum of atmospheric phase fluctuations is assumed to follow the von Karman prescription, but possibly modified by the addition of a critical scale below which the power is zero. (See the VonKarman docstring for more details).

As an expectation value, this profile is formally only exact in the infinite-exposure limit. However, at least for large apertures, we have found that this expectation value is approached rapidly, and can be applied for even fairly small exposure times.

The intended use for this profile is as a correction to applying the geometric approximation to PhaseScreenPSF objects when drawing using geometric photon shooting. In this case, the PhaseScreenPSF will simulate the effects of the low frequency turbulence modes, which can be treated purely using refraction, while the SecondKick handles the high frequency modes.

The geometric approximation is only valid for length scales larger than some critical scale where the effects of interference are unimportant. For smaller length scales, interference (diffraction) must be handled using an optical paradigm that acknowledges the wave nature of light, such as Fourier optics.

Fourier optics calculations are many orders of magnitude slower than geometric optics calculations for typical flux levels, however, so we implement a scale-splitting algorithm first described in Peterson et al. (2015) for the LSST PhoSim package. Essentially, phase fluctuations below a critical mode in Fourier space, labeled kcrit, are handled by the fast geometric optics calculations present in PhaseScreenPSF. Fluctuations for Fourier modes above kcrit are then calculated analytically by SecondKick. Because very many oscillations of these high-k modes both fit within a given telescope aperture and pass by the aperture during a moderate length exposure time, we can use the same analytic expectation value calculation for the high-k component of all PSFs across a field of view, thus incurring the somewhat expensive calculation for Fourier optics only once.

There are two limiting cases for this profile that may helpful for readers trying to understand how this class works. When kcrit = 0, then all turbulent modes are included, and this surface brightness profile becomes identical to the convolution of an Airy profile and a Von Karman profile. In contrast, when kcrit = inf, then none of the turbulent modes are included, and this surface brightness profile is just an Airy profile. In other words, the full effect of an Airy profile, and additionally some portion (which depends on kcrit) of a VonKarman profile are modeled.

For more details, we refer the reader to the original implementation described in

Peterson et al. 2015 ApJSS vol. 218

Parameters:
  • lam – Wavelength in nanometers

  • r0 – Fried parameter in meters.

  • diam – Aperture diameter in meters.

  • obscuration – Linear dimension of central obscuration as fraction of aperture linear dimension. [0., 1.). [default: 0.0]

  • kcrit – Critical Fourier mode (in units of 1/r0) below which the turbulence power spectrum will be truncated. [default: 0.2]

  • flux – The flux (in photons/cm^2/s) of the profile. [default: 1]

  • scale_unit – Units assumed when drawing this profile or evaluating xValue, kValue, etc. Should be a galsim.AngleUnit or a string that can be used to construct one (e.g., ‘arcsec’, ‘radians’, etc.). [default: galsim.arcsec]

  • gsparams – An optional GSParams argument. [default: None]

property diam

The input diam value.

property kcrit

The input kcrit value.

property lam

The input lam value.

property obscuration

The input obscuration value.

property r0

The input r0 value.

property scale_unit

The input scale_unit value.

withFlux(flux)[source]

Create a version of the current object with a different flux.

This function is equivalent to obj.withScaledFlux(flux / obj.flux).

It creates a new object that has the same profile as the original, but with the surface brightness at every location rescaled such that the total flux will be the given value. Note that if flux is an SED, the return value will be a ChromaticObject with specified SED.

Parameters:

flux – The new flux for the object.

Returns:

the object with the new flux

galsim.Atmosphere(screen_size, rng=None, _bar=None, **kwargs)[source]

Create an atmosphere as a list of turbulent phase screens at different altitudes. The atmosphere model can then be used to simulate atmospheric PSFs.

Simulating an atmospheric PSF is typically accomplished by first representing the 3-dimensional turbulence in the atmosphere as a series of discrete 2-dimensional phase screens. These screens may blow around in the wind, and may or may not also evolve in time. This function allows one to quickly assemble a list of atmospheric phase screens into a galsim.PhaseScreenList object, which can then be used to evaluate PSFs through various columns of atmosphere at different field angles.

The atmospheric screens currently available represent turbulence following a von Karman power spectrum. Specifically, the phase power spectrum in each screen can be written

\[\psi(\nu) = 0.023 r_0^{-5/3} \left(\nu^2 + \frac{1}{L_0^2}\right)^{11/6}\]

where \(\psi(\nu)\) is the power spectral density at spatial frequency \(\nu\), \(r_0\) is the Fried parameter (which has dimensions of length) and sets the amplitude of the turbulence, and \(L_0\) is the outer scale (also dimensions of length) beyond which the power asymptotically flattens.

Typical values for \(r_0\) are ~0.1 to 0.2 meters, which corresponds roughly to PSF FWHMs of ~0.5 to 1.0 arcsec for optical wavelengths. Note that \(r_0\) is a function of wavelength, scaling like \(r_0 \sim \lambda^{6/5}\). To reduce confusion, the input parameter here is named r0_500 and refers explicitly to the Fried parameter at a wavelength of 500 nm. The outer scale is typically in the 10s of meters and does not vary with wavelength.

To create multiple layers, simply specify keyword arguments as length-N lists instead of scalars (works for all arguments except rng). If, for any of these keyword arguments, you want to use the same value for each layer, then you can just specify the argument as a scalar and the function will automatically broadcast it into a list with length equal to the longest found keyword argument list. Note that it is an error to specify keywords with lists of different lengths (unless only one of them has length > 1).

The one exception to the above is the keyword r0_500. The effective Fried parameter for a set of atmospheric layers is:

r0_500_effective = (sum(r**(-5./3) for r in r0_500s))**(-3./5)

Providing r0_500 as a scalar or single-element list will result in broadcasting such that the effective Fried parameter for the whole set of layers equals the input argument. You can weight the contribution of each layer with the r0_weights keyword.

As an example, the following code approximately creates the atmosphere used by Jee+Tyson(2011) for their study of atmospheric PSFs for LSST. Note this code takes about ~2 minutes to run on a fast laptop, and will consume about (8192**2 pixels) * (8 bytes) * (6 screens) ~ 3 GB of RAM in its final state, and more at intermediate states.:

>>> altitude = [0, 2.58, 5.16, 7.73, 12.89, 15.46]  # km
>>> r0_500 = 0.16  # m
>>> weights = [0.652, 0.172, 0.055, 0.025, 0.074, 0.022]
>>> speed = np.random.uniform(0, 20, size=6)  # m/s
>>> direction = [np.random.uniform(0, 360)*galsim.degrees for i in range(6)]
>>> npix = 8192
>>> screen_scale = r0_500
>>> atm = galsim.Atmosphere(r0_500=r0_500, r0_weights=weights,
                            screen_size=screen_scale*npix,
                            altitude=altitude, L0=25.0, speed=speed,
                            direction=direction, screen_scale=screen_scale)

Once the atmosphere is constructed, a 15-sec exposure length, 5ms time step, monochromatic PSF at 700nm (using an 8.4 meter aperture, 0.6 fractional obscuration and otherwise default settings) takes about 7 minutes to draw on a fast laptop.:

>>> psf = atm.makePSF(lam=700.0, exptime=15.0, time_step=0.005, diam=8.4, obscuration=0.6)
>>> img1 = psf.drawImage()  # ~7 min

The same psf, if drawn using photon-shooting on the same laptop, will generate photons at a rate of about 1 million per second.:

>>> img2 = psf.drawImage(nx=32, ny=32, scale=0.2, method='phot', n_photons=1e6)  # ~1 sec.

Note that the Fourier-based calculation compute time will scale linearly with exposure time, while the photon-shooting calculation compute time will scale linearly with the number of photons being shot.

Many factors will affect the timing of results, of course, including aperture diameter, gsparams settings, pad_factor and oversampling options to makePSF, time_step and exposure time, frozen vs. non-frozen atmospheric layers, and so on. We recommend that users try varying these settings to find a balance of speed and fidelity.

Parameters:
  • r0_500 – Fried parameter setting the amplitude of turbulence; contributes to “size” of the resulting atmospheric PSF. Specified at wavelength 500 nm, in units of meters. [default: 0.2]

  • r0_weights – Weights for splitting up the contribution of r0_500 between different layers. Note that this keyword is only allowed if r0_500 is either a scalar or a single-element list. [default: None]

  • screen_size – Physical extent of square phase screen in meters. This should be large enough to accommodate the desired field-of-view of the telescope as well as the meta-pupil defined by the wind speed and exposure time. Note that the screen will have periodic boundary conditions, so the code will run with a smaller sized screen, though this may introduce artifacts into PSFs or PSF correlation functions. Note that screen_size may be tweaked by the initializer to ensure screen_size is a multiple of screen_scale.

  • screen_scale – Physical pixel scale of phase screen in meters. A fraction of the Fried parameter is usually sufficiently small, but users should test the effects of this parameter to ensure robust results. [default: same as each screen’s r0_500]

  • altitude – Altitude of phase screen in km. This is with respect to the telescope, not sea-level. [default: 0.0]

  • L0 – Outer scale in meters. The turbulence power spectrum will smoothly approach a constant at scales larger than L0. Set to None or np.inf for a power spectrum without an outer scale. [default: 25.0]

  • speed – Wind speed in meters/second. [default: 0.0]

  • direction – Wind direction as galsim.Angle [default: 0.0 * galsim.degrees]

  • alpha – Square root of fraction of phase that is “remembered” between time_steps (i.e., alpha**2 is the fraction remembered). The fraction sqrt(1-alpha**2) is then the amount of turbulence freshly generated in each step. Setting alpha=1.0 results in a frozen-flow atmosphere. Note that computing PSFs from frozen-flow atmospheres may be significantly faster than computing PSFs with non-frozen-flow atmospheres. [default: 1.0]

  • time_step – Time interval between phase boiling updates. Note that this is distinct from the time interval used when integrating the PSF over time, which is set by the time_step keyword argument to PhaseScreenPSF or PhaseScreenList.makePSF. If time_step is not None, then it is required that alpha is set to something other than 1.0. [default: None]

  • rng – Random number generator as a BaseDeviate. If None, then use the clock time or system entropy to seed a new generator. [default: None]

galsim.phase_screens.initWorkerArgs()[source]

Function used to generate worker arguments to pass to multiprocessing.Pool initializer.

See AtmosphericScreen docstring for more information.

galsim.phase_screens.initWorker(share)[source]

Worker initialization function to pass to multiprocessing.Pool initializer.

See AtmosphericScreen docstring for more information.

galsim.phase_screens.reset_shared_screens()[source]

Reset the global dict that contains screens being shared across multiprocessing processes.

This is almost never necessary. However, if you first use one multiprocessing context with initWorker and initWorkerArgs, and then switch to a different context for another multiprocessing action, the dict we use for storing the shared memory will not be usable in the new context. In this case, you should reset the global dict before starting the second multiprocessing action by running:

galsim.phase_screens.reset_shared_screens()

If you only ever use one multiprocessing context in your program, you should never need to call this.