FACT++  1.0
palPolmo.c
Go to the documentation of this file.
1 /*
2 *+
3 * Name:
4 * palPolmo
5 
6 * Purpose:
7 * Correct for polar motion
8 
9 * Language:
10 * Starlink ANSI C
11 
12 * Type of Module:
13 * Library routine
14 
15 * Invocation:
16 * palPolmo ( double elongm, double phim, double xp, double yp,
17 * double *elong, double *phi, double *daz );
18 
19 * Arguments:
20 * elongm = double (Given)
21 * Mean logitude of the observer (radians, east +ve)
22 * phim = double (Given)
23 * Mean geodetic latitude of the observer (radians)
24 * xp = double (Given)
25 * Polar motion x-coordinate (radians)
26 * yp = double (Given)
27 * Polar motion y-coordinate (radians)
28 * elong = double * (Returned)
29 * True longitude of the observer (radians, east +ve)
30 * phi = double * (Returned)
31 * True geodetic latitude of the observer (radians)
32 * daz = double * (Returned)
33 * Azimuth correction (terrestrial-celestial, radians)
34 
35 * Description:
36 * Polar motion: correct site longitude and latitude for polar
37 * motion and calculate azimuth difference between celestial and
38 * terrestrial poles.
39 
40 * Authors:
41 * PTW: Patrick Wallace (STFC)
42 * TIMJ: Tim Jenness (Cornell)
43 * {enter_new_authors_here}
44 
45 * Notes:
46 * - "Mean" longitude and latitude are the (fixed) values for the
47 * site's location with respect to the IERS terrestrial reference
48 * frame; the latitude is geodetic. TAKE CARE WITH THE LONGITUDE
49 * SIGN CONVENTION. The longitudes used by the present routine
50 * are east-positive, in accordance with geographical convention
51 * (and right-handed). In particular, note that the longitudes
52 * returned by the sla_OBS routine are west-positive, following
53 * astronomical usage, and must be reversed in sign before use in
54 * the present routine.
55 *
56 * - XP and YP are the (changing) coordinates of the Celestial
57 * Ephemeris Pole with respect to the IERS Reference Pole.
58 * XP is positive along the meridian at longitude 0 degrees,
59 * and YP is positive along the meridian at longitude
60 * 270 degrees (i.e. 90 degrees west). Values for XP,YP can
61 * be obtained from IERS circulars and equivalent publications;
62 * the maximum amplitude observed so far is about 0.3 arcseconds.
63 *
64 * - "True" longitude and latitude are the (moving) values for
65 * the site's location with respect to the celestial ephemeris
66 * pole and the meridian which corresponds to the Greenwich
67 * apparent sidereal time. The true longitude and latitude
68 * link the terrestrial coordinates with the standard celestial
69 * models (for precession, nutation, sidereal time etc).
70 *
71 * - The azimuths produced by sla_AOP and sla_AOPQK are with
72 * respect to due north as defined by the Celestial Ephemeris
73 * Pole, and can therefore be called "celestial azimuths".
74 * However, a telescope fixed to the Earth measures azimuth
75 * essentially with respect to due north as defined by the
76 * IERS Reference Pole, and can therefore be called "terrestrial
77 * azimuth". Uncorrected, this would manifest itself as a
78 * changing "azimuth zero-point error". The value DAZ is the
79 * correction to be added to a celestial azimuth to produce
80 * a terrestrial azimuth.
81 *
82 * - The present routine is rigorous. For most practical
83 * purposes, the following simplified formulae provide an
84 * adequate approximation:
85 *
86 * elong = elongm+xp*cos(elongm)-yp*sin(elongm)
87 * phi = phim+(xp*sin(elongm)+yp*cos(elongm))*tan(phim)
88 * daz = -sqrt(xp*xp+yp*yp)*cos(elongm-atan2(xp,yp))/cos(phim)
89 *
90 * An alternative formulation for DAZ is:
91 *
92 * x = cos(elongm)*cos(phim)
93 * y = sin(elongm)*cos(phim)
94 * daz = atan2(-x*yp-y*xp,x*x+y*y)
95 *
96 * - Reference: Seidelmann, P.K. (ed), 1992. "Explanatory Supplement
97 * to the Astronomical Almanac", ISBN 0-935702-68-7,
98 * sections 3.27, 4.25, 4.52.
99 
100 * History:
101 * 2000-11-30 (PTW):
102 * SLALIB implementation.
103 * 2014-10-18 (TIMJ):
104 * Initial version in C.
105 * {enter_further_changes_here}
106 
107 * Copyright:
108 * Copyright (C) 2000 Rutherford Appleton Laboratory.
109 * Copyright (C) 2014 Cornell University
110 * All Rights Reserved.
111 
112 * Licence:
113 * This program is free software; you can redistribute it and/or
114 * modify it under the terms of the GNU General Public License as
115 * published by the Free Software Foundation; either version 3 of
116 * the License, or (at your option) any later version.
117 *
118 * This program is distributed in the hope that it will be
119 * useful, but WITHOUT ANY WARRANTY; without even the implied
120 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
121 * PURPOSE. See the GNU General Public License for more details.
122 *
123 * You should have received a copy of the GNU General Public License
124 * along with this program. If not, see <http://www.gnu.org/licenses/>.
125 
126 * Bugs:
127 * {note_any_bugs_here}
128 *-
129 */
130 
131 #include <math.h>
132 
133 #include "pal.h"
134 
135 void palPolmo ( double elongm, double phim, double xp, double yp,
136  double *elong, double *phi, double *daz ) {
137 
138  double sel,cel,sph,cph,xm,ym,zm,xnm,ynm,znm,
139  sxp,cxp,syp,cyp,zw,xt,yt,zt,xnt,ynt;
140 
141  /* Site mean longitude and mean geodetic latitude as a Cartesian vector */
142  sel=sin(elongm);
143  cel=cos(elongm);
144  sph=sin(phim);
145  cph=cos(phim);
146 
147  xm=cel*cph;
148  ym=sel*cph;
149  zm=sph;
150 
151  /* Rotate site vector by polar motion, Y-component then X-component */
152  sxp=sin(xp);
153  cxp=cos(xp);
154  syp=sin(yp);
155  cyp=cos(yp);
156 
157  zw=(-ym*syp+zm*cyp);
158 
159  xt=xm*cxp-zw*sxp;
160  yt=ym*cyp+zm*syp;
161  zt=xm*sxp+zw*cxp;
162 
163  /* Rotate also the geocentric direction of the terrestrial pole (0,0,1) */
164  xnm=-sxp*cyp;
165  ynm=syp;
166  znm=cxp*cyp;
167 
168  cph=sqrt(xt*xt+yt*yt);
169  if (cph == 0.0) xt=1.0;
170  sel=yt/cph;
171  cel=xt/cph;
172 
173  /* Return true longitude and true geodetic latitude of site */
174  if (xt != 0.0 || yt != 0.0) {
175  *elong=atan2(yt,xt);
176  } else {
177  *elong=0.0;
178  }
179  *phi=atan2(zt,cph);
180 
181  /* Return current azimuth of terrestrial pole seen from site position */
182  xnt=(xnm*cel+ynm*sel)*zt-znm*cph;
183  ynt=-xnm*sel+ynm*cel;
184  if (xnt != 0.0 || ynt != 0.0) {
185  *daz=atan2(-ynt,-xnt);
186  } else {
187  *daz=0.0;
188  }
189 
190 }
void palPolmo(double elongm, double phim, double xp, double yp, double *elong, double *phi, double *daz)
Definition: palPolmo.c:135