diff options
Diffstat (limited to 'math/slalib/refz.f')
| -rw-r--r-- | math/slalib/refz.f | 170 |
1 files changed, 170 insertions, 0 deletions
diff --git a/math/slalib/refz.f b/math/slalib/refz.f new file mode 100644 index 00000000..6bdae80f --- /dev/null +++ b/math/slalib/refz.f @@ -0,0 +1,170 @@ + SUBROUTINE slREFZ (ZU, REFA, REFB, ZR) +*+ +* - - - - - +* R E F Z +* - - - - - +* +* Adjust an unrefracted zenith distance to include the effect of +* atmospheric refraction, using the simple A tan Z + B tan**3 Z +* model (plus special handling for large ZDs). +* +* Given: +* ZU dp unrefracted zenith distance of the source (radian) +* REFA dp tan Z coefficient (radian) +* REFB dp tan**3 Z coefficient (radian) +* +* Returned: +* ZR dp refracted zenith distance (radian) +* +* Notes: +* +* 1 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 formula used here is based on the +* Newton-Raphson method. For the utmost numerical consistency +* with the refracted to unrefracted model, two iterations are +* carried out, achieving agreement at the 1D-11 arcseconds level +* for a ZD of 80 degrees. The inherent accuracy of the model +* is, of course, far worse than this - see the documentation for +* slRFCO for more information. +* +* 2 At ZD 83 degrees, the rapidly-worsening A tan Z + B tan^3 Z +* model is abandoned and an empirical formula takes over. 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: +* +* ZR error +* +* 80 0.7 +* 81 1.3 +* 82 2.4 +* 83 4.7 +* 84 6.2 +* 85 6.4 +* 86 8 +* 87 10 +* 88 15 +* 89 30 +* 90 60 +* 91 150 } relevant only to +* 92 400 } high-elevation sites +* +* For radio wavelengths the errors are typically 50% larger than +* the optical figures and by ZD 85 deg are twice as bad, worsening +* rapidly below that. To maintain 1 arcsec accuracy down to ZD=85 +* at the Green Bank site, Condon (2004) has suggested amplifying +* the amount of refraction predicted by slREFZ below 10.8 deg +* elevation by the factor (1+0.00195*(10.8-E_t)), where E_t is the +* unrefracted elevation in degrees. +* +* The high-ZD model is scaled to match the normal model at the +* transition point; there is no glitch. +* +* 3 Beyond 93 deg zenith distance, the refraction is held at its +* 93 deg value. +* +* 4 See also the routine slREFV, which performs the adjustment in +* Cartesian Az/El coordinates, and with the emphasis on speed +* rather than numerical accuracy. +* +* Reference: +* +* Condon,J.J., Refraction Corrections for the GBT, PTCS/PN/35.2, +* NRAO Green Bank, 2004. +* +* P.T.Wallace Starlink 9 April 2004 +* +* Copyright (C) 2004 Rutherford Appleton Laboratory +* +* License: +* 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 2 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 (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +* Copyright (C) 1995 Association of Universities for Research in Astronomy Inc. +*- + + IMPLICIT NONE + + DOUBLE PRECISION ZU,REFA,REFB,ZR + +* Radians to degrees + DOUBLE PRECISION R2D + PARAMETER (R2D=57.29577951308232D0) + +* Largest usable ZD (deg) + DOUBLE PRECISION D93 + PARAMETER (D93=93D0) + +* Coefficients for high ZD model (used beyond ZD 83 deg) + DOUBLE PRECISION C1,C2,C3,C4,C5 + PARAMETER (C1=+0.55445D0, + : C2=-0.01133D0, + : C3=+0.00202D0, + : C4=+0.28385D0, + : C5=+0.02390D0) + +* ZD at which one model hands over to the other (radians) + DOUBLE PRECISION Z83 + PARAMETER (Z83=83D0/R2D) + +* High-ZD-model prediction (deg) for that point + DOUBLE PRECISION REF83 + PARAMETER (REF83=(C1+C2*7D0+C3*49D0)/(1D0+C4*7D0+C5*49D0)) + + DOUBLE PRECISION ZU1,ZL,S,C,T,TSQ,TCU,REF,E,E2 + + + +* Perform calculations for ZU or 83 deg, whichever is smaller + ZU1 = MIN(ZU,Z83) + +* Functions of ZD + ZL = ZU1 + S = SIN(ZL) + C = COS(ZL) + T = S/C + TSQ = T*T + TCU = T*TSQ + +* Refracted ZD (mathematically to better than 1 mas at 70 deg) + ZL = ZL-(REFA*T+REFB*TCU)/(1D0+(REFA+3D0*REFB*TSQ)/(C*C)) + +* Further iteration + S = SIN(ZL) + C = COS(ZL) + T = S/C + TSQ = T*T + TCU = T*TSQ + REF = ZU1-ZL+ + : (ZL-ZU1+REFA*T+REFB*TCU)/(1D0+(REFA+3D0*REFB*TSQ)/(C*C)) + +* Special handling for large ZU + IF (ZU.GT.ZU1) THEN + E = 90D0-MIN(D93,ZU*R2D) + E2 = E*E + REF = (REF/REF83)*(C1+C2*E+C3*E2)/(1D0+C4*E+C5*E2) + END IF + +* Return refracted ZD + ZR = ZU-REF + + END |