22 from astropy
import units
23 from astropy.units
import cds
27 from ._coord
import Weather
29 __all__ = [
"refraction",
"differentialRefraction"]
33 deltaRefractScale = 1.0E8
36 def refraction(wavelength, elevation, observatory, weather=None):
37 """Calculate overall refraction under atmospheric and observing conditions.
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.
55 refraction : `lsst.geom.Angle`
56 The angular refraction for light of the given wavelength,
57 under the given observing conditions.
61 The calculation is taken from [1]_.
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.
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()
77 reducedN =
deltaN(wavelength, weather)/deltaRefractScale
79 atmosScaleheightRatio = 4.5908E-6*temperature.to_value(units.Kelvin)
83 relativeGravity = (1. + 0.005302*np.sin(latitude.asRadians())**2.
84 - 0.00000583*np.sin(2.*latitude.asRadians())**2. - 0.000000315*altitude)
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
95 """Calculate the differential refraction between two wavelengths.
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.
115 differentialRefraction : `lsst.geom.Angle`
116 The refraction at `wavelength` minus the refraction at `wavelengthRef`.
118 refractionStart =
refraction(wavelength, elevation, observatory, weather=weather)
119 refractionEnd =
refraction(wavelengthRef, elevation, observatory, weather=weather)
120 return refractionStart - refractionEnd
124 """Calculate the differential refractive index of air.
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
137 The difference of the refractive index of air from 1.,
138 calculated as (n_air - 1)*10^8
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
146 This replicates equation 14 of [1]_
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.
154 waveNum = 1E3/wavelength
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.
161 """Calculate dry air pressure term to refractive index calculation.
165 weather : `lsst.afw.coord.Weather`
166 Class containing the measured temperature, pressure, and humidity
167 at the observatory during an observation
171 densityFactor : `float`
172 Returns the relative density of dry air
173 at the given pressure and temperature.
177 This replicates equation 15 of [1]_
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.
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)
196 """Calculate water vapor pressure term to refractive index calculation.
200 weather : `lsst.afw.coord.Weather`
201 Class containing the measured temperature, pressure, and humidity
202 at the observatory during an observation
206 densityFactor : `float`
207 Returns the relative density of water vapor
208 at the given pressure and temperature.
212 This replicates equation 16 of [1]_
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.
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
233 """Convert humidity and temperature to water vapor pressure.
237 weather : `lsst.afw.coord.Weather`
238 Class containing the measured temperature, pressure, and humidity
239 at the observatory during an observation
243 pressure : `astropy.units.Quantity`
244 The water vapor pressure in Pascals
245 calculated from the given humidity and temperature.
249 This replicates equations 18 & 20 of [1]_
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.
257 humidity = weather.getHumidity()
258 x = np.log(humidity/100.0)
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
266 return waterVaporPressure
270 """Thin wrapper to return the measured temperature from an observation.
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.
284 temperature : `astropy.units.Quantity`
285 The temperature in Celsius, unless `useKelvin` is set.
287 temperature = weather.getAirTemperature()*units.Celsius
289 temperature = temperature.to(units.Kelvin, equivalencies=units.temperature())
294 """Set default local weather conditions if they are missing.
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.
306 default : `lsst.afw.coord.Weather`
307 Updated Weather class with any `nan` values replaced by defaults.
309 if isinstance(altitude, units.quantity.Quantity):
312 altitude2 = altitude*units.meter
313 p0 = 101325.*units.pascal
314 g = 9.80665*units.meter/units.second**2
315 R0 = 8.31447*units.Joule/(units.mol*units.Kelvin)
316 T0 = 19.*units.Celsius
317 lapseRate = -6.5*units.Celsius/units.km
318 M = 0.0289644*units.kg/units.mol
320 temperature = T0 + lapseRate*altitude2
321 temperatureK = temperature.to(units.Kelvin, equivalencies=units.temperature())
322 pressure = p0*np.exp(-(g*M*altitude2)/(R0*temperatureK))
324 weather = Weather((temperature/units.Celsius).value, (pressure/units.pascal).value, humidity)
def extractTemperature(weather, useKelvin=False)
def densityFactorDry(weather)
def defaultWeather(altitude)
def humidityToPressure(weather)
def deltaN(wavelength, weather)
def refraction(wavelength, elevation, observatory, weather=None)
def differentialRefraction(wavelength, wavelengthRef, elevation, observatory, weather=None)
def densityFactorWater(weather)