/*
*+
* Name:
* palPolmo
* Purpose:
* Correct for polar motion
* Language:
* Starlink ANSI C
* Type of Module:
* Library routine
* Invocation:
* palPolmo ( double elongm, double phim, double xp, double yp,
* double *elong, double *phi, double *daz );
* Arguments:
* elongm = double (Given)
* Mean logitude of the observer (radians, east +ve)
* phim = double (Given)
* Mean geodetic latitude of the observer (radians)
* xp = double (Given)
* Polar motion x-coordinate (radians)
* yp = double (Given)
* Polar motion y-coordinate (radians)
* elong = double * (Returned)
* True longitude of the observer (radians, east +ve)
* phi = double * (Returned)
* True geodetic latitude of the observer (radians)
* daz = double * (Returned)
* Azimuth correction (terrestrial-celestial, radians)
* Description:
* Polar motion: correct site longitude and latitude for polar
* motion and calculate azimuth difference between celestial and
* terrestrial poles.
* Authors:
* PTW: Patrick Wallace (STFC)
* TIMJ: Tim Jenness (Cornell)
* {enter_new_authors_here}
* Notes:
* - "Mean" longitude and latitude are the (fixed) values for the
* site's location with respect to the IERS terrestrial reference
* frame; the latitude is geodetic. TAKE CARE WITH THE LONGITUDE
* SIGN CONVENTION. The longitudes used by the present routine
* are 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.
*
* - XP and YP are the (changing) coordinates of the Celestial
* Ephemeris Pole with respect to the IERS Reference Pole.
* XP is positive along the meridian at longitude 0 degrees,
* and YP is positive along the meridian at longitude
* 270 degrees (i.e. 90 degrees west). Values for XP,YP can
* be obtained from IERS circulars and equivalent publications;
* the maximum amplitude observed so far is about 0.3 arcseconds.
*
* - "True" longitude and latitude are the (moving) values for
* the site's location with respect to the celestial ephemeris
* pole and the meridian which corresponds to the Greenwich
* apparent sidereal time. The true longitude and latitude
* link the terrestrial coordinates with the standard celestial
* models (for precession, nutation, sidereal time etc).
*
* - The azimuths produced by sla_AOP and sla_AOPQK are with
* respect to due north as defined by the Celestial Ephemeris
* Pole, and can therefore be called "celestial azimuths".
* However, a telescope fixed to the Earth measures azimuth
* essentially with respect to due north as defined by the
* IERS Reference Pole, and can therefore be called "terrestrial
* azimuth". Uncorrected, this would manifest itself as a
* changing "azimuth zero-point error". The value DAZ is the
* correction to be added to a celestial azimuth to produce
* a terrestrial azimuth.
*
* - The present routine is rigorous. For most practical
* purposes, the following simplified formulae provide an
* adequate approximation:
*
* elong = elongm+xp*cos(elongm)-yp*sin(elongm)
* phi = phim+(xp*sin(elongm)+yp*cos(elongm))*tan(phim)
* daz = -sqrt(xp*xp+yp*yp)*cos(elongm-atan2(xp,yp))/cos(phim)
*
* An alternative formulation for DAZ is:
*
* x = cos(elongm)*cos(phim)
* y = sin(elongm)*cos(phim)
* daz = atan2(-x*yp-y*xp,x*x+y*y)
*
* - Reference: Seidelmann, P.K. (ed), 1992. "Explanatory Supplement
* to the Astronomical Almanac", ISBN 0-935702-68-7,
* sections 3.27, 4.25, 4.52.
* History:
* 2000-11-30 (PTW):
* SLALIB implementation.
* 2014-10-18 (TIMJ):
* Initial version in C.
* {enter_further_changes_here}
* Copyright:
* Copyright (C) 2000 Rutherford Appleton Laboratory.
* Copyright (C) 2014 Cornell University
* 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, see <http://www.gnu.org/licenses/>.
* Bugs:
* {note_any_bugs_here}
*-
*/
#include <math.h>
#include "pal.h"
void palPolmo ( double elongm, double phim, double xp, double yp,
double *elong, double *phi, double *daz ) {
double sel,cel,sph,cph,xm,ym,zm,xnm,ynm,znm,
sxp,cxp,syp,cyp,zw,xt,yt,zt,xnt,ynt;
/* Site mean longitude and mean geodetic latitude as a Cartesian vector */
sel=sin(elongm);
cel=cos(elongm);
sph=sin(phim);
cph=cos(phim);
xm=cel*cph;
ym=sel*cph;
zm=sph;
/* Rotate site vector by polar motion, Y-component then X-component */
sxp=sin(xp);
cxp=cos(xp);
syp=sin(yp);
cyp=cos(yp);
zw=(-ym*syp+zm*cyp);
xt=xm*cxp-zw*sxp;
yt=ym*cyp+zm*syp;
zt=xm*sxp+zw*cxp;
/* Rotate also the geocentric direction of the terrestrial pole (0,0,1) */
xnm=-sxp*cyp;
ynm=syp;
znm=cxp*cyp;
cph=sqrt(xt*xt+yt*yt);
if (cph == 0.0) xt=1.0;
sel=yt/cph;
cel=xt/cph;
/* Return true longitude and true geodetic latitude of site */
if (xt != 0.0 || yt != 0.0) {
*elong=atan2(yt,xt);
} else {
*elong=0.0;
}
*phi=atan2(zt,cph);
/* Return current azimuth of terrestrial pole seen from site position */
xnt=(xnm*cel+ynm*sel)*zt-znm*cph;
ynt=-xnm*sel+ynm*cel;
if (xnt != 0.0 || ynt != 0.0) {
*daz=atan2(-ynt,-xnt);
} else {
*daz=0.0;
}
}