/*
*+
* Name:
* palAtmdsp
* Purpose:
* Apply atmospheric-dispersion adjustments to refraction coefficients
* Language:
* Starlink ANSI C
* Type of Module:
* Library routine
* Invocation:
* void palAtmdsp( double tdk, double pmb, double rh, double wl1,
* double a1, double b1, double wl2, double *a2, double *b2 );
* Arguments:
* tdk = double (Given)
* Ambient temperature, K
* pmb = double (Given)
* Ambient pressure, millibars
* rh = double (Given)
* Ambient relative humidity, 0-1
* wl1 = double (Given)
* Reference wavelength, micrometre (0.4 recommended)
* a1 = double (Given)
* Refraction coefficient A for wavelength wl1 (radians)
* b1 = double (Given)
* Refraction coefficient B for wavelength wl1 (radians)
* wl2 = double (Given)
* Wavelength for which adjusted A,B required
* a2 = double * (Returned)
* Refraction coefficient A for wavelength WL2 (radians)
* b2 = double * (Returned)
* Refraction coefficient B for wavelength WL2 (radians)
* Description:
* Apply atmospheric-dispersion adjustments to refraction coefficients.
* Authors:
* TIMJ: Tim Jenness
* PTW: Patrick Wallace
* {enter_new_authors_here}
* Notes:
* - To use this routine, first call palRefco specifying WL1 as the
* wavelength. This yields refraction coefficients A1,B1, correct
* for that wavelength. Subsequently, calls to palAtmdsp specifying
* different wavelengths will produce new, slightly adjusted
* refraction coefficients which apply to the specified wavelength.
*
* - Most of the atmospheric dispersion happens between 0.7 micrometre
* and the UV atmospheric cutoff, and the effect increases strongly
* towards the UV end. For this reason a blue reference wavelength
* is recommended, for example 0.4 micrometres.
*
* - The accuracy, for this set of conditions:
*
* height above sea level 2000 m
* latitude 29 deg
* pressure 793 mb
* temperature 17 degC
* humidity 50%
* lapse rate 0.0065 degC/m
* reference wavelength 0.4 micrometre
* star elevation 15 deg
*
* is about 2.5 mas RMS between 0.3 and 1.0 micrometres, and stays
* within 4 mas for the whole range longward of 0.3 micrometres
* (compared with a total dispersion from 0.3 to 20.0 micrometres
* of about 11 arcsec). These errors are typical for ordinary
* conditions and the given elevation; in extreme conditions values
* a few times this size may occur, while at higher elevations the
* errors become much smaller.
*
* - If either wavelength exceeds 100 micrometres, the radio case
* is assumed and the returned refraction coefficients are the
* same as the given ones. Note that radio refraction coefficients
* cannot be turned into optical values using this routine, nor
* vice versa.
*
* - The algorithm consists of calculation of the refractivity of the
* air at the observer for the two wavelengths, using the methods
* of the palRefro routine, and then scaling of the two refraction
* coefficients according to classical refraction theory. This
* amounts to scaling the A coefficient in proportion to (n-1) and
* the B coefficient almost in the same ratio (see R.M.Green,
* "Spherical Astronomy", Cambridge University Press, 1985).
* History:
* 2014-07-15 (TIMJ):
* Initial version. A direct copy of the Fortran SLA implementation.
* Adapted with permission from the Fortran SLALIB library.
* {enter_further_changes_here}
* Copyright:
* Copyright (C) 2014 Tim Jenness
* Copyright (C) 2005 Patrick Wallace
* All Rights Reserved.
* Licence:
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation; either version 3 of
* the License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be
* useful, but WITHOUT ANY WARRANTY; without even the implied
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
* PURPOSE. See the GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
* MA 02110-1301, USA.
* Bugs:
* {note_any_bugs_here}
*-
*/
#include "pal.h"
#include "palmac.h"
#include <math.h>
void palAtmdsp ( double tdk, double pmb, double rh, double wl1,
double a1, double b1, double wl2, double *a2, double *b2 ) {
double f,tdkok,pmbok,rhok;
double psat,pwo,w1,wlok,wlsq,w2,dn1,dn2;
/* Check for radio wavelengths */
if (wl1 > 100.0 || wl2 > 100.0) {
/* Radio: no dispersion */
*a2 = a1;
*b2 = b1;
} else {
/* Optical: keep arguments within safe bounds */
tdkok = DMIN(DMAX(tdk,100.0),500.0);
pmbok = DMIN(DMAX(pmb,0.0),10000.0);
rhok = DMIN(DMAX(rh,0.0),1.0);
/* Atmosphere parameters at the observer */
psat = pow(10.0, -8.7115+0.03477*tdkok);
pwo = rhok*psat;
w1 = 11.2684e-6*pwo;
/* Refractivity at the observer for first wavelength */
wlok = DMAX(wl1,0.1);
wlsq = wlok*wlok;
w2 = 77.5317e-6+(0.43909e-6+0.00367e-6/wlsq)/wlsq;
dn1 = (w2*pmbok-w1)/tdkok;
/* Refractivity at the observer for second wavelength */
wlok = DMAX(wl2,0.1);
wlsq = wlok*wlok;
w2 = 77.5317e-6+(0.43909e-6+0.00367e-6/wlsq)/wlsq;
dn2 = (w2*pmbok-w1)/tdkok;
/* Scale the refraction coefficients (see Green 4.31, p93) */
if (dn1 != 0.0) {
f = dn2/dn1;
*a2 = a1*f;
*b2 = b1*f;
if (dn1 != a1) {
*b2 *= (1.0+dn1*(dn1-dn2)/(2.0*(dn1-a1)));
}
} else {
*a2 = a1;
*b2 = b1;
}
}
}