Initial commit

1 files changed, 242 insertions, 0 deletions

diff --git a/src/slalib/aopqk.f b/src/slalib/aopqk.f
new file mode 100644
index 0000000..9f8a45c
--- /dev/null
+++ b/src/slalib/aopqk.f
@@ -0,0 +1,242 @@
+      SUBROUTINE sla_AOPQK (RAP, DAP, AOPRMS, AOB, ZOB, HOB, DOB, ROB)
+*+
+*     - - - - - -
+*      A O P Q K
+*     - - - - - -
+*
+*  Quick apparent to observed place (but see note 8, below, for
+*  remarks about speed).
+*
+*  Given:
+*     RAP    d      geocentric apparent right ascension
+*     DAP    d      geocentric apparent declination
+*     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 (T)
+*       (7)      pressure (P)
+*       (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)
+*
+*  Returned:
+*     AOB    d      observed azimuth (radians: N=0,E=90)
+*     ZOB    d      observed zenith distance (radians)
+*     HOB    d      observed Hour Angle (radians)
+*     DOB    d      observed Declination (radians)
+*     ROB    d      observed Right Ascension (radians)
+*
+*  Notes:
+*
+*   1)  This routine returns zenith distance rather than elevation
+*       in order to reflect the fact that no allowance is made for
+*       depression of the horizon.
+*
+*   2)  The accuracy of the result is limited by the corrections for
+*       refraction.  Providing the meteorological parameters are
+*       known accurately and there are no gross local effects, the
+*       observed RA,Dec predicted by this routine should be within
+*       about 0.1 arcsec for a zenith distance of less than 70 degrees.
+*       Even at a topocentric zenith distance of 90 degrees, the
+*       accuracy in elevation should be better than 1 arcmin;  useful
+*       results are available for a further 3 degrees, beyond which
+*       the slaRefro routine returns a fixed value of the refraction.
+*       The complementary routines slaAop (or slaAopqk) and slaOap
+*       (or slaOapqk) are self-consistent to better than 1 micro-
+*       arcsecond all over the celestial sphere.
+*
+*   3)  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.
+*
+*   4)  "Apparent" place means the geocentric apparent right ascension
+*       and declination, which is obtained from a catalogue mean place
+*       by allowing for space motion, parallax, precession, nutation,
+*       annual aberration, and the Sun's gravitational lens effect.  For
+*       star positions in the FK5 system (i.e. J2000), these effects can
+*       be applied by means of the sla_MAP etc routines.  Starting from
+*       other mean place systems, additional transformations will be
+*       needed;  for example, FK4 (i.e. B1950) mean places would first
+*       have to be converted to FK5, which can be done with the
+*       sla_FK425 etc routines.
+*
+*   5)  "Observed" Az,El means the position that would be seen by a
+*       perfect theodolite located at the observer.  This is obtained
+*       from the geocentric apparent RA,Dec by allowing for Earth
+*       orientation and diurnal aberration, rotating from equator
+*       to horizon coordinates, and then adjusting for refraction.
+*       The HA,Dec is obtained by rotating back into equatorial
+*       coordinates, using the geodetic latitude corrected for polar
+*       motion, and is the position that would be seen by a perfect
+*       equatorial located at the observer and with its polar axis
+*       aligned to the Earth's axis of rotation (n.b. not to the
+*       refracted pole).  Finally, the RA is obtained by subtracting
+*       the HA from the local apparent ST.
+*
+*   6)  To predict the required setting of a real telescope, the
+*       observed place produced by this routine would have to be
+*       adjusted for the tilt of the azimuth or polar axis of the
+*       mounting (with appropriate corrections for mount flexures),
+*       for non-perpendicularity between the mounting axes, for the
+*       position of the rotator axis and the pointing axis relative
+*       to it, for tube flexure, for gear and encoder errors, and
+*       finally for encoder zero points.  Some telescopes would, of
+*       course, exhibit other properties which would need to be
+*       accounted for at the appropriate point in the sequence.
+*
+*   7)  The star-independent apparent-to-observed-place parameters
+*       in AOPRMS may be computed by means of the sla_AOPPA routine.
+*       If nothing has changed significantly except the time, the
+*       sla_AOPPAT routine may be used to perform the requisite
+*       partial recomputation of AOPRMS.
+*
+*   8)  At zenith distances beyond about 76 degrees, the need for
+*       special care with the corrections for refraction causes a
+*       marked increase in execution time.  Moreover, the effect
+*       gets worse with increasing zenith distance.  Adroit
+*       programming in the calling application may allow the
+*       problem to be reduced.  Prepare an alternative AOPRMS array,
+*       computed for zero air-pressure;  this will disable the
+*       refraction corrections and cause rapid execution.  Using
+*       this AOPRMS array, a preliminary call to the present routine
+*       will, depending on the application, produce a rough position
+*       which may be enough to establish whether the full, slow
+*       calculation (using the real AOPRMS array) is worthwhile.
+*       For example, there would be no need for the full calculation
+*       if the preliminary call had already established that the
+*       source was well below the elevation limits for a particular
+*       telescope.
+*
+*  9)   The azimuths etc produced by the present routine are with
+*       respect to the celestial pole.  Corrections to the terrestrial
+*       pole can be computed using sla_POLMO.
+*
+*  Called:  sla_DCS2C, sla_REFZ, sla_REFRO, sla_DCC2S, sla_DRANRM
+*
+*  P.T.Wallace   Starlink   22 February 1996
+*
+*  Copyright (C) 1996 Rutherford Appleton Laboratory
+*-
+
+      IMPLICIT NONE
+
+      DOUBLE PRECISION RAP,DAP,AOPRMS(14),AOB,ZOB,HOB,DOB,ROB
+
+*  Breakpoint for fast/slow refraction algorithm:
+*  ZD greater than arctan(4), (see sla_REFCO routine)
+*  or vector Z less than cosine(arctan(Z)) = 1/sqrt(17)
+      DOUBLE PRECISION ZBREAK
+      PARAMETER (ZBREAK=0.242535625D0)
+
+      INTEGER I
+
+      DOUBLE PRECISION SPHI,CPHI,ST,V(3),XHD,YHD,ZHD,DIURAB,F,
+     :                 XHDT,YHDT,ZHDT,XAET,YAET,ZAET,AZOBS,
+     :                 ZDT,REFA,REFB,ZDOBS,DZD,DREF,CE,
+     :                 XAEO,YAEO,ZAEO,HMOBS,DCOBS,RAOBS
+
+      DOUBLE PRECISION sla_DRANRM
+
+
+
+*  Sin, cos of latitude
+      SPHI = AOPRMS(2)
+      CPHI = AOPRMS(3)
+
+*  Local apparent sidereal time
+      ST = AOPRMS(14)
+
+*  Apparent RA,Dec to Cartesian -HA,Dec
+      CALL sla_DCS2C(RAP-ST,DAP,V)
+      XHD = V(1)
+      YHD = V(2)
+      ZHD = V(3)
+
+*  Diurnal aberration
+      DIURAB = AOPRMS(4)
+      F = (1D0-DIURAB*YHD)
+      XHDT = F*XHD
+      YHDT = F*(YHD+DIURAB)
+      ZHDT = F*ZHD
+
+*  Cartesian -HA,Dec to Cartesian Az,El (S=0,E=90)
+      XAET = SPHI*XHDT-CPHI*ZHDT
+      YAET = YHDT
+      ZAET = CPHI*XHDT+SPHI*ZHDT
+
+*  Azimuth (N=0,E=90)
+      IF (XAET.EQ.0D0.AND.YAET.EQ.0D0) THEN
+         AZOBS = 0D0
+      ELSE
+         AZOBS = ATAN2(YAET,-XAET)
+      END IF
+
+*  Topocentric zenith distance
+      ZDT = ATAN2(SQRT(XAET*XAET+YAET*YAET),ZAET)
+
+*
+*  Refraction
+*  ----------
+
+*  Fast algorithm using two constant model
+      REFA = AOPRMS(11)
+      REFB = AOPRMS(12)
+      CALL sla_REFZ(ZDT,REFA,REFB,ZDOBS)
+
+*  Large zenith distance?
+      IF (COS(ZDOBS).LT.ZBREAK) THEN
+
+*     Yes: use rigorous algorithm
+
+*     Initialize loop (maximum of 10 iterations)
+         I = 1
+         DZD = 1D1
+         DO WHILE (ABS(DZD).GT.1D-10.AND.I.LE.10)
+
+*        Compute refraction using current estimate of observed ZD
+            CALL sla_REFRO(ZDOBS,AOPRMS(5),AOPRMS(6),AOPRMS(7),
+     :                     AOPRMS(8),AOPRMS(9),AOPRMS(1),
+     :                     AOPRMS(10),1D-8,DREF)
+
+*        Remaining discrepancy
+            DZD = ZDOBS+DREF-ZDT
+
+*        Update the estimate
+            ZDOBS = ZDOBS-DZD
+
+*        Increment the iteration counter
+            I = I+1
+         END DO
+      END IF
+
+*  To Cartesian Az/ZD
+      CE = SIN(ZDOBS)
+      XAEO = -COS(AZOBS)*CE
+      YAEO = SIN(AZOBS)*CE
+      ZAEO = COS(ZDOBS)
+
+*  Cartesian Az/ZD to Cartesian -HA,Dec
+      V(1) = SPHI*XAEO+CPHI*ZAEO
+      V(2) = YAEO
+      V(3) = -CPHI*XAEO+SPHI*ZAEO
+
+*  To spherical -HA,Dec
+      CALL sla_DCC2S(V,HMOBS,DCOBS)
+
+*  Right Ascension
+      RAOBS = sla_DRANRM(ST+HMOBS)
+
+*  Return the results
+      AOB = AZOBS
+      ZOB = ZDOBS
+      HOB = -HMOBS
+      DOB = DCOBS
+      ROB = RAOBS
+
+      END