diff options
Diffstat (limited to 'src/slalib/aoppa.f')
| -rw-r--r-- | src/slalib/aoppa.f | 176 |
1 files changed, 176 insertions, 0 deletions
diff --git a/src/slalib/aoppa.f b/src/slalib/aoppa.f new file mode 100644 index 0000000..b497fdb --- /dev/null +++ b/src/slalib/aoppa.f @@ -0,0 +1,176 @@ + SUBROUTINE sla_AOPPA (DATE, DUT, ELONGM, PHIM, HM, + : XP, YP, TDK, PMB, RH, WL, TLR, AOPRMS) +*+ +* - - - - - - +* A O P P A +* - - - - - - +* +* Precompute apparent to observed place parameters required by +* sla_AOPQK and sla_OAPQK. +* +* Given: +* DATE d UTC date/time (modified Julian Date, JD-2400000.5) +* DUT d delta UT: UT1-UTC (UTC seconds) +* ELONGM d mean longitude of the observer (radians, east +ve) +* PHIM d mean geodetic latitude of the observer (radians) +* HM d observer's height above sea level (metres) +* XP d polar motion x-coordinate (radians) +* YP d polar motion y-coordinate (radians) +* TDK d local ambient temperature (DegK; std=273.155D0) +* PMB d local atmospheric pressure (mB; std=1013.25D0) +* RH d local relative humidity (in the range 0D0-1D0) +* WL d effective wavelength (micron, e.g. 0.55D0) +* TLR d tropospheric lapse rate (DegK/metre, e.g. 0.0065D0) +* +* Returned: +* AOPRMS d(14) star-independent apparent-to-observed parameters: +* +* (1) geodetic latitude (radians) +* (2,3) sine and cosine of geodetic latitude +* (4) magnitude of diurnal aberration vector +* (5) height (HM) +* (6) ambient temperature (TDK) +* (7) pressure (PMB) +* (8) relative humidity (RH) +* (9) wavelength (WL) +* (10) lapse rate (TLR) +* (11,12) refraction constants A and B (radians) +* (13) longitude + eqn of equinoxes + sidereal DUT (radians) +* (14) local apparent sidereal time (radians) +* +* Notes: +* +* 1) It is advisable to take great care with units, as even +* unlikely values of the input parameters are accepted and +* processed in accordance with the models used. +* +* 2) The DATE argument is UTC expressed as an MJD. This is, +* strictly speaking, improper, because of leap seconds. However, +* as long as the delta UT and the UTC are consistent there +* are no difficulties, except during a leap second. In this +* case, the start of the 61st second of the final minute should +* begin a new MJD day and the old pre-leap delta UT should +* continue to be used. As the 61st second completes, the MJD +* should revert to the start of the day as, simultaneously, +* the delta UTC changes by one second to its post-leap new value. +* +* 3) The delta UT (UT1-UTC) is tabulated in IERS circulars and +* elsewhere. It increases by exactly one second at the end of +* each UTC leap second, introduced in order to keep delta UT +* within +/- 0.9 seconds. +* +* 4) IMPORTANT -- TAKE CARE WITH THE LONGITUDE SIGN CONVENTION. +* The longitude required by the present routine is east-positive, +* in accordance with geographical convention (and right-handed). +* In particular, note that the longitudes returned by the +* sla_OBS routine are west-positive, following astronomical +* usage, and must be reversed in sign before use in the present +* routine. +* +* 5) The polar coordinates XP,YP can be obtained from IERS +* circulars and equivalent publications. The maximum amplitude +* is about 0.3 arcseconds. If XP,YP values are unavailable, +* use XP=YP=0D0. See page B60 of the 1988 Astronomical Almanac +* for a definition of the two angles. +* +* 6) The height above sea level of the observing station, HM, +* can be obtained from the Astronomical Almanac (Section J +* in the 1988 edition), or via the routine sla_OBS. If P, +* the pressure in millibars, is available, an adequate +* estimate of HM can be obtained from the expression +* +* HM ~ -29.3D0*TSL*LOG(P/1013.25D0). +* +* where TSL is the approximate sea-level air temperature in +* deg K (see Astrophysical Quantities, C.W.Allen, 3rd edition, +* section 52). Similarly, if the pressure P is not known, +* it can be estimated from the height of the observing +* station, HM as follows: +* +* P ~ 1013.25D0*EXP(-HM/(29.3D0*TSL)). +* +* Note, however, that the refraction is proportional to the +* pressure and that an accurate P value is important for +* precise work. +* +* 7) Repeated, computationally-expensive, calls to sla_AOPPA for +* times that are very close together can be avoided by calling +* sla_AOPPA just once and then using sla_AOPPAT for the subsequent +* times. Fresh calls to sla_AOPPA will be needed only when changes +* in the precession have grown to unacceptable levels or when +* anything affecting the refraction has changed. +* +* Called: sla_GEOC, sla_REFCO, sla_EQEQX, sla_AOPPAT +* +* P.T.Wallace Starlink 6 September 1999 +* +* Copyright (C) 1999 P.T.Wallace and CCLRC +*- + + IMPLICIT NONE + + DOUBLE PRECISION DATE,DUT,ELONGM,PHIM,HM,XP,YP,TDK,PMB, + : RH,WL,TLR,AOPRMS(14) + + DOUBLE PRECISION sla_EQEQX + +* 2Pi + DOUBLE PRECISION D2PI + PARAMETER (D2PI=6.283185307179586476925287D0) + +* Seconds of time to radians + DOUBLE PRECISION S2R + PARAMETER (S2R=7.272205216643039903848712D-5) + +* Speed of light (AU per day) + DOUBLE PRECISION C + PARAMETER (C=173.14463331D0) + +* Ratio between solar and sidereal time + DOUBLE PRECISION SOLSID + PARAMETER (SOLSID=1.00273790935D0) + + DOUBLE PRECISION CPHIM,XT,YT,ZT,XC,YC,ZC,ELONG,PHI,UAU,VAU + + + +* Observer's location corrected for polar motion + CPHIM = COS(PHIM) + XT = COS(ELONGM)*CPHIM + YT = SIN(ELONGM)*CPHIM + ZT = SIN(PHIM) + XC = XT-XP*ZT + YC = YT+YP*ZT + ZC = XP*XT-YP*YT+ZT + IF (XC.EQ.0D0.AND.YC.EQ.0D0) THEN + ELONG = 0D0 + ELSE + ELONG = ATAN2(YC,XC) + END IF + PHI = ATAN2(ZC,SQRT(XC*XC+YC*YC)) + AOPRMS(1) = PHI + AOPRMS(2) = SIN(PHI) + AOPRMS(3) = COS(PHI) + +* Magnitude of the diurnal aberration vector + CALL sla_GEOC(PHI,HM,UAU,VAU) + AOPRMS(4) = D2PI*UAU*SOLSID/C + +* Copy the refraction parameters and compute the A & B constants + AOPRMS(5) = HM + AOPRMS(6) = TDK + AOPRMS(7) = PMB + AOPRMS(8) = RH + AOPRMS(9) = WL + AOPRMS(10) = TLR + CALL sla_REFCO(HM,TDK,PMB,RH,WL,PHI,TLR,1D-10, + : AOPRMS(11),AOPRMS(12)) + +* Longitude + equation of the equinoxes + sidereal equivalent of DUT +* (ignoring change in equation of the equinoxes between UTC and TDB) + AOPRMS(13) = ELONG+sla_EQEQX(DATE)+DUT*SOLSID*S2R + +* Sidereal time + CALL sla_AOPPAT(DATE,AOPRMS) + + END |