# 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.
#
from .errors import GalSimIncompatibleValuesError
from .celestial import CelestialCoord
from .angle import degrees, Angle
[docs]def air_refractive_index_minus_one(wave, pressure=69.328, temperature=293.15, H2O_pressure=1.067):
"""Return the refractive index of air as function of wavelength.
Uses the formulae given in Filippenko (1982), which appear to come from Edlen (1953),
and Coleman, Bozman, and Meggers (1960). The units of the original formula are non-SI,
being mmHg for pressure (and water vapor pressure), and degrees C for temperature. This
function accepts SI units, however, and transforms them when plugging into the formula.
The default values for temperature, pressure and water vapor pressure are expected to be
appropriate for LSST at Cerro Pachon, Chile, but they are broadly reasonable for most
observatories.
Parameters:
wave: Wavelength array in nanometers
pressure: Air pressure in kiloPascals.
temperature: Temperature in Kelvins.
H2O_pressure: Water vapor pressure in kiloPascals.
Returns:
the refractive index minus 1.
"""
P = pressure * 7.50061683 # kPa -> mmHg
T = temperature - 273.15 # K -> C
W = H2O_pressure * 7.50061683 # kPa -> mmHg
sigma_squared = 1.0 / (wave * 1.e-3)**2.0 # inverse wavenumber squared in micron^-2
n_minus_one = (64.328 + (29498.1 / (146.0 - sigma_squared))
+ (255.4 / (41.0 - sigma_squared))) * 1.e-6
n_minus_one *= P * (1.0 + (1.049 - 0.0157 * T) * 1.e-6 * P) / (720.883 * (1.0 + 0.003661 * T))
n_minus_one -= (0.0624 - 0.000680 * sigma_squared)/(1.0 + 0.003661 * T) * W * 1.e-6
return n_minus_one
[docs]def get_refraction(wave, zenith_angle, **kwargs):
"""Compute the angle of refraction for a photon entering the atmosphere.
Photons refract when transitioning from space, where the refractive index n = 1.0 exactly, to
air, where the refractive index is slightly greater than 1.0. This function computes the
change in zenith angle for a photon with a given wavelength. Output is a positive number of
radians, even though the apparent zenith angle technically decreases due to this effect.
Parameters:
wave: Wavelength array in nanometers
zenith_angle: `Angle` from object to zenith
**kwargs: Keyword arguments to pass to air_refractive_index() to override default
pressure, temperature, and/or H2O_pressure.
Returns:
the absolute value of change in zenith angle in radians.
"""
nm1 = air_refractive_index_minus_one(wave, **kwargs)
# The following line is equivalent to:
# n_squared = (nm1 + 1)**2
# r0 = (n_squared - 1.0) / (2.0 * n_squared)
r0 = nm1 * (nm1+2) / 2.0 / (nm1**2 + 2*nm1 + 1)
return r0 * zenith_angle.tan()
[docs]def zenith_parallactic_angles(obj_coord, zenith_coord=None, HA=None, latitude=None):
"""Compute the zenith angle and parallactic angle of a celestial coordinate, given either
the celestial coordinate of the zenith, or equivalently, the hour angle of the coordinate and
the latitude of the observer. This is useful for the function `ChromaticAtmosphere`.
Parameters:
obj_coord: A `CelestialCoord` object for which the zenith and parallactic
angles will be computed.
zenith_coord: A `CelestialCoord` indicating the coordinates of the zenith.
HA: The hour angle (as an `Angle`) of the coordinate for which the
zenith and parallactic angles will be computed.
latitude: The observer's latitude, as an `Angle`.
Returns:
the tuple (zenith_angle, parallactic_angle), each of which is an `Angle`.
"""
if zenith_coord is None:
if HA is None or latitude is None:
raise GalSimIncompatibleValuesError(
"Must provide either zenith_coord or (HA, latitude).",
HA=HA, latitude=latitude, zenith_coord=zenith_coord)
zenith_coord = CelestialCoord(HA + obj_coord.ra, latitude)
else:
if HA is not None or latitude is not None:
raise GalSimIncompatibleValuesError(
"Cannot provide both zenith_coord and (HA, latitude).",
HA=HA, latitude=latitude, zenith_coord=zenith_coord)
zenith_angle = obj_coord.distanceTo(zenith_coord)
NCP = CelestialCoord(0.0*degrees, 90*degrees)
parallactic_angle = obj_coord.angleBetween(NCP, zenith_coord)
return zenith_angle, parallactic_angle
[docs]def parse_dcr_angles(**kwargs):
"""Parse the various options for specifying the zenith angle and parallactic angle
in `ChromaticAtmosphere`.
Parameters:
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]
**kwargs: For convenience, any other kwargs are returned back for further
processing.
Returns:
zenith_angle, parallactic_angle, kw, where kw is any other kwargs that aren't relevant.
"""
if 'zenith_angle' in kwargs:
zenith_angle = kwargs.pop('zenith_angle')
parallactic_angle = kwargs.pop('parallactic_angle', 0.0*degrees)
if not isinstance(zenith_angle, Angle):
raise TypeError("zenith_angle must be a galsim.Angle.")
if not isinstance(parallactic_angle, Angle):
raise TypeError("parallactic_angle must be a galsim.Angles.")
elif 'obj_coord' in kwargs:
obj_coord = kwargs.pop('obj_coord')
zenith_coord = kwargs.pop('zenith_coord', None)
HA = kwargs.pop('HA', None)
latitude = kwargs.pop('latitude', None)
zenith_angle, parallactic_angle = zenith_parallactic_angles(
obj_coord=obj_coord, zenith_coord=zenith_coord, HA=HA, latitude=latitude)
else:
raise TypeError("Need to specify zenith_angle and parallactic_angle.")
return zenith_angle, parallactic_angle, kwargs