LSST Applications g180d380827+6621f76652,g2079a07aa2+86d27d4dc4,g2305ad1205+f5a9e323a1,g2bbee38e9b+c6a8a0fb72,g337abbeb29+c6a8a0fb72,g33d1c0ed96+c6a8a0fb72,g3a166c0a6a+c6a8a0fb72,g3ddfee87b4+9a10e1fe7b,g48712c4677+c9a099281a,g487adcacf7+f2e03ea30b,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+aead732c78,g64a986408d+eddffb812c,g858d7b2824+eddffb812c,g864b0138d7+aa38e45daa,g974c55ee3d+f37bf00e57,g99cad8db69+119519a52d,g9c22b2923f+e2510deafe,g9ddcbc5298+9a081db1e4,ga1e77700b3+03d07e1c1f,gb0e22166c9+60f28cb32d,gb23b769143+eddffb812c,gba4ed39666+c2a2e4ac27,gbb8dafda3b+27317ec8e9,gbd998247f1+585e252eca,gc120e1dc64+5817c176a8,gc28159a63d+c6a8a0fb72,gc3e9b769f7+6707aea8b4,gcf0d15dbbd+9a10e1fe7b,gdaeeff99f8+f9a426f77a,ge6526c86ff+6a2e01d432,ge79ae78c31+c6a8a0fb72,gee10cc3b42+585e252eca,gff1a9f87cc+eddffb812c,v27.0.0.rc1
LSST Data Management Base Package
Loading...
Searching...
No Matches
_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
22from astropy import units
23from astropy.units import cds
24import numpy as np
25
26import lsst.geom
27from ._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.
33deltaRefractScale = 1.0E8
34
35
36def 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
94def 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
123def 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
160def 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
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
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
269def 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
293def 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
extractTemperature(weather, useKelvin=False)
deltaN(wavelength, weather)
differentialRefraction(wavelength, wavelengthRef, elevation, observatory, weather=None)
refraction(wavelength, elevation, observatory, weather=None)