From 1a47fd282af077e2eccc2824a42dc23c9d962103 Mon Sep 17 00:00:00 2001 From: Tim Jenness Date: Tue, 15 Jul 2014 18:21:16 -0700 Subject: [PATCH] Add palRefv and palAtmdsp --- Makefile.am | 2 + palAtmdsp.c | 180 ++++++++++++++++++++++++++++++++++++++++++++++++++++ palRefv.c | 155 ++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 337 insertions(+) create mode 100644 palAtmdsp.c create mode 100644 palRefv.c diff --git a/Makefile.am b/Makefile.am index 6105883..3641750 100644 --- a/Makefile.am +++ b/Makefile.am @@ -40,6 +40,7 @@ palAop.c \ palAoppa.c \ palAoppat.c \ palAopqk.c \ +palAtmdsp.c \ palCaldj.c \ palDafin.c \ palDe2h.c \ @@ -99,6 +100,7 @@ palPvobs.c \ palRdplan.c \ palRefco.c \ palRefro.c \ +palRefv.c \ palRefz.c \ palRverot.c \ palRvgalc.c \ diff --git a/palAtmdsp.c b/palAtmdsp.c new file mode 100644 index 0000000..0227797 --- /dev/null +++ b/palAtmdsp.c @@ -0,0 +1,180 @@ +/* +*+ +* 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 + +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; + } + } + +} diff --git a/palRefv.c b/palRefv.c new file mode 100644 index 0000000..85b3761 --- /dev/null +++ b/palRefv.c @@ -0,0 +1,155 @@ +/* +*+ +* Name: +* palRefv + +* Purpose: +* Adjust an unrefracted Cartesian vector to include the effect of atmospheric refraction + +* Language: +* Starlink ANSI C + +* Type of Module: +* Library routine + +* Invocation: +* void palRefv ( double vu[3], double refa, double refb, double vr[3] ); + +* Arguments: +* vu[3] = double (Given) +* Unrefracted position of the source (Az/El 3-vector) +* refa = double (Given) +* tan Z coefficient (radian) +* refb = double (Given) +* tan**3 Z coefficient (radian) +* vr[3] = double (Returned) +* Refracted position of the source (Az/El 3-vector) + +* Description: +* Adjust an unrefracted Cartesian vector to include the effect of +* atmospheric refraction, using the simple A tan Z + B tan**3 Z +* model. + +* Authors: +* TIMJ: Tim Jenness +* PTW: Patrick Wallace +* {enter_new_authors_here} + +* Notes: +* - This routine applies the adjustment for refraction in the +* opposite sense to the usual one - it takes an unrefracted +* (in vacuo) position and produces an observed (refracted) +* position, whereas the A tan Z + B tan**3 Z model strictly +* applies to the case where an observed position is to have the +* refraction removed. The unrefracted to refracted case is +* harder, and requires an inverted form of the text-book +* refraction models; the algorithm used here is equivalent to +* one iteration of the Newton-Raphson method applied to the above +* formula. +* +* - Though optimized for speed rather than precision, the present +* routine achieves consistency with the refracted-to-unrefracted +* A tan Z + B tan**3 Z model at better than 1 microarcsecond within +* 30 degrees of the zenith and remains within 1 milliarcsecond to +* beyond ZD 70 degrees. The inherent accuracy of the model is, of +* course, far worse than this - see the documentation for sla_REFCO +* for more information. +* +* - At low elevations (below about 3 degrees) the refraction +* correction is held back to prevent arithmetic problems and +* wildly wrong results. For optical/IR wavelengths, over a wide +* range of observer heights and corresponding temperatures and +* pressures, the following levels of accuracy (arcsec, worst case) +* are achieved, relative to numerical integration through a model +* atmosphere: +* +* ZD error +* +* 80 0.7 +* 81 1.3 +* 82 2.5 +* 83 5 +* 84 10 +* 85 20 +* 86 55 +* 87 160 +* 88 360 +* 89 640 +* 90 1100 +* 91 1700 } relevant only to +* 92 2600 } high-elevation sites +* +* The results for radio are slightly worse over most of the range, +* becoming significantly worse below ZD=88 and unusable beyond +* ZD=90. +* +* - See also the routine palRefz, which performs the adjustment to +* the zenith distance rather than in Cartesian Az/El coordinates. +* The present routine is faster than palRefz and, except very low down, +* is equally accurate for all practical purposes. However, beyond +* about ZD 84 degrees palRefz should be used, and for the utmost +* accuracy iterative use of palRefro should be considered. + +* 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) 2004 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 + +void palRefv ( double vu[3], double refa, double refb, double vr[3] ) { + + double x,y,z1,z,zsq,rsq,r,wb,wt,d,cd,f; + + /* Initial estimate = unrefracted vector */ + x = vu[0]; + y = vu[1]; + z1 = vu[2]; + + /* Keep correction approximately constant below about 3 deg elevation */ + z = DMAX(z1,0.05); + + /* One Newton-Raphson iteration */ + zsq = z*z; + rsq = x*x+y*y; + r = sqrt(rsq); + wb = refb*rsq/zsq; + wt = (refa+wb)/(1.0+(refa+3.0*wb)*(zsq+rsq)/zsq); + d = wt*r/z; + cd = 1.0-d*d/2.0; + f = cd*(1.0-wt); + + /* Post-refraction x,y,z */ + vr[0] = x*f; + vr[1] = y*f; + vr[2] = cd*(z+d*r)+(z1-z); +}