source: trunk/FACT++/pal/palAtmdsp.c@ 18501

Last change on this file since 18501 was 18347, checked in by tbretz, 9 years ago
File size: 6.0 KB
Line 
1/*
2*+
3* Name:
4* palAtmdsp
5
6* Purpose:
7* Apply atmospheric-dispersion adjustments to refraction coefficients
8
9* Language:
10* Starlink ANSI C
11
12* Type of Module:
13* Library routine
14
15* Invocation:
16* void palAtmdsp( double tdk, double pmb, double rh, double wl1,
17* double a1, double b1, double wl2, double *a2, double *b2 );
18
19
20* Arguments:
21* tdk = double (Given)
22* Ambient temperature, K
23* pmb = double (Given)
24* Ambient pressure, millibars
25* rh = double (Given)
26* Ambient relative humidity, 0-1
27* wl1 = double (Given)
28* Reference wavelength, micrometre (0.4 recommended)
29* a1 = double (Given)
30* Refraction coefficient A for wavelength wl1 (radians)
31* b1 = double (Given)
32* Refraction coefficient B for wavelength wl1 (radians)
33* wl2 = double (Given)
34* Wavelength for which adjusted A,B required
35* a2 = double * (Returned)
36* Refraction coefficient A for wavelength WL2 (radians)
37* b2 = double * (Returned)
38* Refraction coefficient B for wavelength WL2 (radians)
39
40* Description:
41* Apply atmospheric-dispersion adjustments to refraction coefficients.
42
43* Authors:
44* TIMJ: Tim Jenness
45* PTW: Patrick Wallace
46* {enter_new_authors_here}
47
48* Notes:
49* - To use this routine, first call palRefco specifying WL1 as the
50* wavelength. This yields refraction coefficients A1,B1, correct
51* for that wavelength. Subsequently, calls to palAtmdsp specifying
52* different wavelengths will produce new, slightly adjusted
53* refraction coefficients which apply to the specified wavelength.
54*
55* - Most of the atmospheric dispersion happens between 0.7 micrometre
56* and the UV atmospheric cutoff, and the effect increases strongly
57* towards the UV end. For this reason a blue reference wavelength
58* is recommended, for example 0.4 micrometres.
59*
60* - The accuracy, for this set of conditions:
61*
62* height above sea level 2000 m
63* latitude 29 deg
64* pressure 793 mb
65* temperature 17 degC
66* humidity 50%
67* lapse rate 0.0065 degC/m
68* reference wavelength 0.4 micrometre
69* star elevation 15 deg
70*
71* is about 2.5 mas RMS between 0.3 and 1.0 micrometres, and stays
72* within 4 mas for the whole range longward of 0.3 micrometres
73* (compared with a total dispersion from 0.3 to 20.0 micrometres
74* of about 11 arcsec). These errors are typical for ordinary
75* conditions and the given elevation; in extreme conditions values
76* a few times this size may occur, while at higher elevations the
77* errors become much smaller.
78*
79* - If either wavelength exceeds 100 micrometres, the radio case
80* is assumed and the returned refraction coefficients are the
81* same as the given ones. Note that radio refraction coefficients
82* cannot be turned into optical values using this routine, nor
83* vice versa.
84*
85* - The algorithm consists of calculation of the refractivity of the
86* air at the observer for the two wavelengths, using the methods
87* of the palRefro routine, and then scaling of the two refraction
88* coefficients according to classical refraction theory. This
89* amounts to scaling the A coefficient in proportion to (n-1) and
90* the B coefficient almost in the same ratio (see R.M.Green,
91* "Spherical Astronomy", Cambridge University Press, 1985).
92
93* History:
94* 2014-07-15 (TIMJ):
95* Initial version. A direct copy of the Fortran SLA implementation.
96* Adapted with permission from the Fortran SLALIB library.
97* {enter_further_changes_here}
98
99* Copyright:
100* Copyright (C) 2014 Tim Jenness
101* Copyright (C) 2005 Patrick Wallace
102* All Rights Reserved.
103
104* Licence:
105* This program is free software; you can redistribute it and/or
106* modify it under the terms of the GNU General Public License as
107* published by the Free Software Foundation; either version 3 of
108* the License, or (at your option) any later version.
109*
110* This program is distributed in the hope that it will be
111* useful, but WITHOUT ANY WARRANTY; without even the implied
112* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
113* PURPOSE. See the GNU General Public License for more details.
114*
115* You should have received a copy of the GNU General Public License
116* along with this program; if not, write to the Free Software
117* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
118* MA 02110-1301, USA.
119
120* Bugs:
121* {note_any_bugs_here}
122*-
123*/
124
125#include "pal.h"
126#include "palmac.h"
127#include <math.h>
128
129void palAtmdsp ( double tdk, double pmb, double rh, double wl1,
130 double a1, double b1, double wl2, double *a2, double *b2 ) {
131
132 double f,tdkok,pmbok,rhok;
133 double psat,pwo,w1,wlok,wlsq,w2,dn1,dn2;
134
135 /* Check for radio wavelengths */
136 if (wl1 > 100.0 || wl2 > 100.0) {
137
138 /* Radio: no dispersion */
139 *a2 = a1;
140 *b2 = b1;
141
142 } else {
143
144 /* Optical: keep arguments within safe bounds */
145 tdkok = DMIN(DMAX(tdk,100.0),500.0);
146 pmbok = DMIN(DMAX(pmb,0.0),10000.0);
147 rhok = DMIN(DMAX(rh,0.0),1.0);
148
149 /* Atmosphere parameters at the observer */
150 psat = pow(10.0, -8.7115+0.03477*tdkok);
151 pwo = rhok*psat;
152 w1 = 11.2684e-6*pwo;
153
154 /* Refractivity at the observer for first wavelength */
155 wlok = DMAX(wl1,0.1);
156 wlsq = wlok*wlok;
157 w2 = 77.5317e-6+(0.43909e-6+0.00367e-6/wlsq)/wlsq;
158 dn1 = (w2*pmbok-w1)/tdkok;
159
160 /* Refractivity at the observer for second wavelength */
161 wlok = DMAX(wl2,0.1);
162 wlsq = wlok*wlok;
163 w2 = 77.5317e-6+(0.43909e-6+0.00367e-6/wlsq)/wlsq;
164 dn2 = (w2*pmbok-w1)/tdkok;
165
166 /* Scale the refraction coefficients (see Green 4.31, p93) */
167 if (dn1 != 0.0) {
168 f = dn2/dn1;
169 *a2 = a1*f;
170 *b2 = b1*f;
171 if (dn1 != a1) {
172 *b2 *= (1.0+dn1*(dn1-dn2)/(2.0*(dn1-a1)));
173 }
174 } else {
175 *a2 = a1;
176 *b2 = b1;
177 }
178 }
179
180}
Note: See TracBrowser for help on using the repository browser.