LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
_refraction.py
Go to the documentation of this file.
1 # This file is part of afw.
2 #
3 # Developed for the LSST Data Management System.
4 # This product includes software developed by the LSST Project
5 # (https://www.lsst.org).
6 # See the COPYRIGHT file at the top-level directory of this distribution
7 # for details of code ownership.
8 #
9 # This program is free software: you can redistribute it and/or modify
10 # it under the terms of the GNU General Public License as published by
11 # the Free Software Foundation, either version 3 of the License, or
12 # (at your option) any later version.
13 #
14 # This program is distributed in the hope that it will be useful,
15 # but WITHOUT ANY WARRANTY; without even the implied warranty of
16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 # GNU General Public License for more details.
18 #
19 # You should have received a copy of the GNU General Public License
20 # along with this program. If not, see <https://www.gnu.org/licenses/>.
21 
22 from astropy import units
23 from astropy.units import cds
24 import numpy as np
25 
26 import lsst.geom
27 from ._coord import Weather
28 
29 __all__ = ["refraction", "differentialRefraction"]
30 
31 # The differential refractive index is the difference of the refractive index from 1.,
32 # multiplied by 1E8 to simplify the notation and equations.
33 deltaRefractScale = 1.0E8
34 
35 
36 def refraction(wavelength, elevation, observatory, weather=None):
37  """Calculate overall refraction under atmospheric and observing conditions.
38 
39  Parameters
40  ----------
41  wavelength : `float`
42  wavelength is in nm (valid for 230.2 < wavelength < 2058.6)
43  elevation : `lsst.geom.Angle`
44  Elevation of the observation, as an Angle.
45  observatory : `lsst.afw.coord.Observatory`
46  Class containing the longitude, latitude,
47  and altitude of the observatory.
48  weather : `lsst.afw.coord.Weather`, optional
49  Class containing the measured temperature, pressure, and humidity
50  at the observatory during an observation
51  If omitted, typical conditions for the observatory's elevation will be calculated.
52 
53  Returns
54  -------
55  refraction : `lsst.geom.Angle`
56  The angular refraction for light of the given wavelength,
57  under the given observing conditions.
58 
59  Notes
60  -----
61  The calculation is taken from [1]_.
62 
63  References
64  ----------
65  .. [1] R. C. Stone, "An Accurate Method for Computing Atmospheric
66  Refraction," Publications of the Astronomical Society of the Pacific,
67  vol. 108, p. 1051, 1996.
68  """
69  if wavelength < 230.2:
70  raise ValueError("Refraction calculation is valid for wavelengths between 230.2 and 2058.6 nm.")
71  if wavelength > 2058.6:
72  raise ValueError("Refraction calculation is valid for wavelengths between 230.2 and 2058.6 nm.")
73  latitude = observatory.getLatitude()
74  altitude = observatory.getElevation()
75  if weather is None:
76  weather = defaultWeather(altitude*units.meter)
77  reducedN = deltaN(wavelength, weather)/deltaRefractScale
78  temperature = extractTemperature(weather, useKelvin=True)
79  atmosScaleheightRatio = 4.5908E-6*temperature.to_value(units.Kelvin)
80 
81  # Account for oblate Earth
82  # This replicates equation 10 of Stone 1996
83  relativeGravity = (1. + 0.005302*np.sin(latitude.asRadians())**2.
84  - 0.00000583*np.sin(2.*latitude.asRadians())**2. - 0.000000315*altitude)
85 
86  # Calculate the tangent of the zenith angle.
87  tanZ = np.tan(np.pi/2. - elevation.asRadians())
88  atmosTerm1 = reducedN*relativeGravity*(1. - atmosScaleheightRatio)
89  atmosTerm2 = reducedN*relativeGravity*(atmosScaleheightRatio - reducedN/2.)
90  result = float(atmosTerm1*tanZ + atmosTerm2*tanZ**3.)*lsst.geom.radians
91  return result
92 
93 
94 def differentialRefraction(wavelength, wavelengthRef, elevation, observatory, weather=None):
95  """Calculate the differential refraction between two wavelengths.
96 
97  Parameters
98  ----------
99  wavelength : `float`
100  wavelength is in nm (valid for 230.2 < wavelength < 2058.6)
101  wavelengthRef : `float`
102  Reference wavelength, typically the effective wavelength of a filter.
103  elevation : `lsst.geom.Angle`
104  Elevation of the observation, as an Angle.
105  observatory : `lsst.afw.coord.Observatory`
106  Class containing the longitude, latitude,
107  and altitude of the observatory.
108  weather : `lsst.afw.coord.Weather`, optional
109  Class containing the measured temperature, pressure, and humidity
110  at the observatory during an observation
111  If omitted, typical conditions for the observatory's elevation will be calculated.
112 
113  Returns
114  -------
115  differentialRefraction : `lsst.geom.Angle`
116  The refraction at `wavelength` minus the refraction at `wavelengthRef`.
117  """
118  refractionStart = refraction(wavelength, elevation, observatory, weather=weather)
119  refractionEnd = refraction(wavelengthRef, elevation, observatory, weather=weather)
120  return refractionStart - refractionEnd
121 
122 
123 def deltaN(wavelength, weather):
124  """Calculate the differential refractive index of air.
125 
126  Parameters
127  ----------
128  wavelength : `float`
129  wavelength is in nanometers
130  weather : `lsst.afw.coord.Weather`
131  Class containing the measured temperature, pressure, and humidity
132  at the observatory during an observation
133 
134  Returns
135  -------
136  deltaN : `float`
137  The difference of the refractive index of air from 1.,
138  calculated as (n_air - 1)*10^8
139 
140  Notes
141  -----
142  The differential refractive index is the difference of
143  the refractive index from 1., multiplied by 1E8 to simplify
144  the notation and equations. Calculated as (n_air - 1)*10^8
145 
146  This replicates equation 14 of [1]_
147 
148  References
149  ----------
150  .. [1] R. C. Stone, "An Accurate Method for Computing Atmospheric
151  Refraction," Publications of the Astronomical Society of the Pacific,
152  vol. 108, p. 1051, 1996.
153  """
154  waveNum = 1E3/wavelength # want wave number in units 1/micron
155  dryAirTerm = 2371.34 + (683939.7/(130. - waveNum**2.)) + (4547.3/(38.9 - waveNum**2.))
156  wetAirTerm = 6487.31 + 58.058*waveNum**2. - 0.71150*waveNum**4. + 0.08851*waveNum**6.
157  return (dryAirTerm*densityFactorDry(weather) + wetAirTerm*densityFactorWater(weather))
158 
159 
160 def densityFactorDry(weather):
161  """Calculate dry air pressure term to refractive index calculation.
162 
163  Parameters
164  ----------
165  weather : `lsst.afw.coord.Weather`
166  Class containing the measured temperature, pressure, and humidity
167  at the observatory during an observation
168 
169  Returns
170  -------
171  densityFactor : `float`
172  Returns the relative density of dry air
173  at the given pressure and temperature.
174 
175  Notes
176  -----
177  This replicates equation 15 of [1]_
178 
179  References
180  ----------
181  .. [1] R. C. Stone, "An Accurate Method for Computing Atmospheric
182  Refraction," Publications of the Astronomical Society of the Pacific,
183  vol. 108, p. 1051, 1996.
184  """
185  temperature = extractTemperature(weather, useKelvin=True)
186  waterVaporPressure = humidityToPressure(weather)
187  airPressure = weather.getAirPressure()*units.pascal
188  dryPressure = airPressure - waterVaporPressure
189  eqn = dryPressure.to_value(cds.mbar)*(57.90E-8 - 9.3250E-4/temperature.to_value(units.Kelvin)
190  + 0.25844/temperature.to_value(units.Kelvin)**2.)
191  densityFactor = (1. + eqn)*dryPressure.to_value(cds.mbar)/temperature.to_value(units.Kelvin)
192  return densityFactor
193 
194 
195 def densityFactorWater(weather):
196  """Calculate water vapor pressure term to refractive index calculation.
197 
198  Parameters
199  ----------
200  weather : `lsst.afw.coord.Weather`
201  Class containing the measured temperature, pressure, and humidity
202  at the observatory during an observation
203 
204  Returns
205  -------
206  densityFactor : `float`
207  Returns the relative density of water vapor
208  at the given pressure and temperature.
209 
210  Notes
211  -----
212  This replicates equation 16 of [1]_
213 
214  References
215  ----------
216  .. [1] R. C. Stone, "An Accurate Method for Computing Atmospheric
217  Refraction," Publications of the Astronomical Society of the Pacific,
218  vol. 108, p. 1051, 1996.
219  """
220  temperature = extractTemperature(weather, useKelvin=True)
221  waterVaporPressure = humidityToPressure(weather)
222  densityEqn1 = (-2.37321E-3 + 2.23366/temperature.to_value(units.Kelvin)
223  - 710.792/temperature.to_value(units.Kelvin)**2.
224  + 7.75141E-4/temperature.to_value(units.Kelvin)**3.)
225  densityEqn2 = waterVaporPressure.to_value(cds.mbar)*(1. + 3.7E-4*waterVaporPressure.to_value(cds.mbar))
226  relativeDensity = waterVaporPressure.to_value(cds.mbar)/temperature.to_value(units.Kelvin)
227  densityFactor = (1 + densityEqn2*densityEqn1)*relativeDensity
228 
229  return densityFactor
230 
231 
232 def humidityToPressure(weather):
233  """Convert humidity and temperature to water vapor pressure.
234 
235  Parameters
236  ----------
237  weather : `lsst.afw.coord.Weather`
238  Class containing the measured temperature, pressure, and humidity
239  at the observatory during an observation
240 
241  Returns
242  -------
243  pressure : `astropy.units.Quantity`
244  The water vapor pressure in Pascals
245  calculated from the given humidity and temperature.
246 
247  Notes
248  -----
249  This replicates equations 18 & 20 of [1]_
250 
251  References
252  ----------
253  .. [1] R. C. Stone, "An Accurate Method for Computing Atmospheric
254  Refraction," Publications of the Astronomical Society of the Pacific,
255  vol. 108, p. 1051, 1996.
256  """
257  humidity = weather.getHumidity()
258  x = np.log(humidity/100.0)
259  temperature = extractTemperature(weather)
260  temperatureEqn1 = (temperature + 238.3*units.Celsius)*x + 17.2694*temperature
261  temperatureEqn2 = (temperature + 238.3*units.Celsius)*(17.2694 - x) - 17.2694*temperature
262  dewPoint = 238.3*temperatureEqn1/temperatureEqn2
263  waterVaporPressure = (4.50874 + 0.341724*dewPoint + 0.0106778*dewPoint**2 + 0.184889E-3*dewPoint**3
264  + 0.238294E-5*dewPoint**4 + 0.203447E-7*dewPoint**5)*133.32239*units.pascal
265 
266  return waterVaporPressure
267 
268 
269 def extractTemperature(weather, useKelvin=False):
270  """Thin wrapper to return the measured temperature from an observation.
271 
272  Parameters
273  ----------
274  weather : `lsst.afw.coord.Weather`
275  Class containing the measured temperature, pressure, and humidity
276  at the observatory during an observation
277  useKelvin : bool, optional
278  Set to True to return the temperature in Kelvin instead of Celsius
279  This is needed because Astropy can't easily convert
280  between Kelvin and Celsius.
281 
282  Returns
283  -------
284  temperature : `astropy.units.Quantity`
285  The temperature in Celsius, unless `useKelvin` is set.
286  """
287  temperature = weather.getAirTemperature()*units.Celsius
288  if useKelvin:
289  temperature = temperature.to(units.Kelvin, equivalencies=units.temperature())
290  return temperature
291 
292 
293 def defaultWeather(altitude):
294  """Set default local weather conditions if they are missing.
295 
296  Parameters
297  ----------
298  weather : `lsst.afw.coord.Weather`
299  Class containing the measured temperature, pressure, and humidity
300  at the observatory during an observation
301  altitude : `astropy.units.Quantity`
302  The altitude of the observatory, in meters.
303 
304  Returns
305  -------
306  default : `lsst.afw.coord.Weather`
307  Updated Weather class with any `nan` values replaced by defaults.
308  """
309  if isinstance(altitude, units.quantity.Quantity):
310  altitude2 = altitude
311  else:
312  altitude2 = altitude*units.meter
313  p0 = 101325.*units.pascal # sea level air pressure
314  g = 9.80665*units.meter/units.second**2 # typical gravitational acceleration at sea level
315  R0 = 8.31447*units.Joule/(units.mol*units.Kelvin) # gas constant
316  T0 = 19.*units.Celsius # Typical sea-level temperature
317  lapseRate = -6.5*units.Celsius/units.km # Typical rate of change of temperature with altitude
318  M = 0.0289644*units.kg/units.mol # molar mass of dry air
319 
320  temperature = T0 + lapseRate*altitude2
321  temperatureK = temperature.to(units.Kelvin, equivalencies=units.temperature())
322  pressure = p0*np.exp(-(g*M*altitude2)/(R0*temperatureK))
323  humidity = 40. # Typical humidity at many observatory sites.
324  weather = Weather((temperature/units.Celsius).value, (pressure/units.pascal).value, humidity)
325  return weather
def extractTemperature(weather, useKelvin=False)
Definition: _refraction.py:269
def deltaN(wavelength, weather)
Definition: _refraction.py:123
def refraction(wavelength, elevation, observatory, weather=None)
Definition: _refraction.py:36
def differentialRefraction(wavelength, wavelengthRef, elevation, observatory, weather=None)
Definition: _refraction.py:94