Permalink
Switch branches/tags
Nothing to show
Find file
Fetching contributors…
Cannot retrieve contributors at this time
163 lines (151 sloc) 4.86 KB
SUBROUTINE sla_ALTAZ (HA, DEC, PHI,
: AZ, AZD, AZDD, EL, ELD, ELDD, PA, PAD, PADD)
*+
* - - - - - -
* A L T A Z
* - - - - - -
*
* Positions, velocities and accelerations for an altazimuth
* telescope mount.
*
* (double precision)
*
* Given:
* HA d hour angle
* DEC d declination
* PHI d observatory latitude
*
* Returned:
* AZ d azimuth
* AZD d " velocity
* AZDD d " acceleration
* EL d elevation
* ELD d " velocity
* ELDD d " acceleration
* PA d parallactic angle
* PAD d " " velocity
* PADD d " " acceleration
*
* Notes:
*
* 1) Natural units are used throughout. HA, DEC, PHI, AZ, EL
* and ZD are in radians. The velocities and accelerations
* assume constant declination and constant rate of change of
* hour angle (as for tracking a star); the units of AZD, ELD
* and PAD are radians per radian of HA, while the units of AZDD,
* ELDD and PADD are radians per radian of HA squared. To
* convert into practical degree- and second-based units:
*
* angles * 360/2pi -> degrees
* velocities * (2pi/86400)*(360/2pi) -> degree/sec
* accelerations * ((2pi/86400)**2)*(360/2pi) -> degree/sec/sec
*
* Note that the seconds here are sidereal rather than SI. One
* sidereal second is about 0.99727 SI seconds.
*
* The velocity and acceleration factors assume the sidereal
* tracking case. Their respective numerical values are (exactly)
* 1/240 and (approximately) 1/3300236.9.
*
* 2) Azimuth is returned in the range 0-2pi; north is zero,
* and east is +pi/2. Elevation and parallactic angle are
* returned in the range +/-pi. Parallactic angle is +ve for
* a star west of the meridian and is the angle NP-star-zenith.
*
* 3) The latitude is geodetic as opposed to geocentric. The
* hour angle and declination are topocentric. Refraction and
* deficiencies in the telescope mounting are ignored. The
* purpose of the routine is to give the general form of the
* quantities. The details of a real telescope could profoundly
* change the results, especially close to the zenith.
*
* 4) No range checking of arguments is carried out.
*
* 5) In applications which involve many such calculations, rather
* than calling the present routine it will be more efficient to
* use inline code, having previously computed fixed terms such
* as sine and cosine of latitude, and (for tracking a star)
* sine and cosine of declination.
*
* This revision: 29 October 2004
*
* Copyright P.T.Wallace. All rights reserved.
*
* 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., 59 Temple Place, Suite 330,
* Boston, MA 02111-1307 USA
*
*-
IMPLICIT NONE
DOUBLE PRECISION HA,DEC,PHI,AZ,AZD,AZDD,EL,ELD,ELDD,PA,PAD,PADD
DOUBLE PRECISION DPI,D2PI,TINY
PARAMETER (DPI=3.1415926535897932384626433832795D0,
: D2PI=6.283185307179586476925286766559D0,
: TINY=1D-30)
DOUBLE PRECISION SH,CH,SD,CD,SP,CP,CHCD,SDCP,X,Y,Z,RSQ,R,A,E,C,S,
: Q,QD,AD,ED,EDR,ADD,EDD,QDD
* Useful functions
SH=SIN(HA)
CH=COS(HA)
SD=SIN(DEC)
CD=COS(DEC)
SP=SIN(PHI)
CP=COS(PHI)
CHCD=CH*CD
SDCP=SD*CP
X=-CHCD*SP+SDCP
Y=-SH*CD
Z=CHCD*CP+SD*SP
RSQ=X*X+Y*Y
R=SQRT(RSQ)
* Azimuth and elevation
IF (RSQ.EQ.0D0) THEN
A=0D0
ELSE
A=ATAN2(Y,X)
END IF
IF (A.LT.0D0) A=A+D2PI
E=ATAN2(Z,R)
* Parallactic angle
C=CD*SP-CH*SDCP
S=SH*CP
IF (C*C+S*S.GT.0) THEN
Q=ATAN2(S,C)
ELSE
Q=DPI-HA
END IF
* Velocities and accelerations (clamped at zenith/nadir)
IF (RSQ.LT.TINY) THEN
RSQ=TINY
R=SQRT(RSQ)
END IF
QD=-X*CP/RSQ
AD=SP+Z*QD
ED=CP*Y/R
EDR=ED/R
ADD=EDR*(Z*SP+(2D0-RSQ)*QD)
EDD=-R*QD*AD
QDD=EDR*(SP+2D0*Z*QD)
* Results
AZ=A
AZD=AD
AZDD=ADD
EL=E
ELD=ED
ELDD=EDD
PA=Q
PAD=QD
PADD=QDD
END