sofars/astro/apco.rs
1use super::{IauAstrom, apcs, pvtob};
2use crate::pnp::c2ixys;
3use crate::vm::{anpm, cr, ir, rx, ry, rz, trxpv};
4
5/// Prepare for ICRS <−> observed, terrestrial, special
6///
7/// For a terrestrial observer, prepare star-independent astrometry
8/// parameters for transformations between ICRS and observed
9/// coordinates. The caller supplies the Earth ephemeris, the Earth
10/// rotation information and the refraction constants as well as the
11/// site coordinates.
12///
13/// This function is part of the International Astronomical Union's
14/// SOFA (Standards of Fundamental Astronomy) software collection.
15///
16/// Status: support function.
17///
18/// Given:
19/// ```
20/// date1 double TDB as a 2-part...
21/// date2 double ...Julian Date (Note 1)
22/// ebpv double[2][3] Earth barycentric PV (au, au/day, Note 2)
23/// ehp double[3] Earth heliocentric P (au, Note 2)
24/// x,y double CIP X,Y (components of unit vector)
25/// s double the CIO locator s (radians)
26/// theta double Earth rotation angle (radians)
27/// elong double longitude (radians, east +ve, Note 3)
28/// phi double latitude (geodetic, radians, Note 3)
29/// hm double height above ellipsoid (m, geodetic, Note 3)
30/// xp,yp double polar motion coordinates (radians, Note 4)
31/// sp double the TIO locator s' (radians, Note 4)
32/// refa double refraction constant A (radians, Note 5)
33/// refb double refraction constant B (radians, Note 5)
34/// ```
35/// Returned:
36/// ```
37/// astrom iauASTROM* star-independent astrometry parameters:
38/// pmt double PM time interval (SSB, Julian years)
39/// eb double[3] SSB to observer (vector, au)
40/// eh double[3] Sun to observer (unit vector)
41/// em double distance from Sun to observer (au)
42/// v double[3] barycentric observer velocity (vector, c)
43/// bm1 double sqrt(1-|v|^2): reciprocal of Lorenz factor
44/// bpn double[3][3] bias-precession-nutation matrix
45/// along double adjusted longitude (radians)
46/// xpl double polar motion xp wrt local meridian (radians)
47/// ypl double polar motion yp wrt local meridian (radians)
48/// sphi double sine of geodetic latitude
49/// cphi double cosine of geodetic latitude
50/// diurab double magnitude of diurnal aberration vector
51/// eral double "local" Earth rotation angle (radians)
52/// refa double refraction constant A (radians)
53/// refb double refraction constant B (radians)
54/// ```
55/// Notes:
56///
57/// 1) The TDB date date1+date2 is a Julian Date, apportioned in any
58/// convenient way between the two arguments. For example,
59/// JD(TDB)=2450123.7 could be expressed in any of these ways, among
60/// others:
61/// ```
62/// date1 date2
63///
64/// 2450123.7 0.0 (JD method)
65/// 2451545.0 -1421.3 (J2000 method)
66/// 2400000.5 50123.2 (MJD method)
67/// 2450123.5 0.2 (date & time method)
68/// ```
69/// The JD method is the most natural and convenient to use in cases
70/// where the loss of several decimal digits of resolution is
71/// acceptable. The J2000 method is best matched to the way the
72/// argument is handled internally and will deliver the optimum
73/// resolution. The MJD method and the date & time methods are both
74/// good compromises between resolution and convenience. For most
75/// applications of this function the choice will not be at all
76/// critical.
77///
78/// TT can be used instead of TDB without any significant impact on
79/// accuracy.
80///
81/// 2) The vectors eb, eh, and all the astrom vectors, are with respect
82/// to BCRS axes.
83///
84/// 3) The geographical coordinates are with respect to the WGS84
85/// reference ellipsoid. TAKE CARE WITH THE LONGITUDE SIGN
86/// CONVENTION: the longitude required by the present function is
87/// right-handed, i.e. east-positive, in accordance with geographical
88/// convention.
89///
90/// The adjusted longitude stored in the astrom array takes into
91/// account the TIO locator and polar motion.
92///
93/// 4) xp and yp are the coordinates (in radians) of the Celestial
94/// Intermediate Pole with respect to the International Terrestrial
95/// Reference System (see IERS Conventions), measured along the
96/// meridians 0 and 90 deg west respectively. sp is the TIO locator
97/// s', in radians, which positions the Terrestrial Intermediate
98/// Origin on the equator. For many applications, xp, yp and
99/// (especially) sp can be set to zero.
100///
101/// Internally, the polar motion is stored in a form rotated onto the
102/// local meridian.
103///
104/// 5) The refraction constants refa and refb are for use in a
105/// dZ = A*tan(Z)+B*tan^3(Z) model, where Z is the observed
106/// (i.e. refracted) zenith distance and dZ is the amount of
107/// refraction.
108///
109/// 6) It is advisable to take great care with units, as even unlikely
110/// values of the input parameters are accepted and processed in
111/// accordance with the models used.
112///
113/// 7) In cases where the caller does not wish to provide the Earth
114/// Ephemeris, the Earth rotation information and refraction
115/// constants, the function iauApco13 can be used instead of the
116/// present function. This starts from UTC and weather readings etc.
117/// and computes suitable values using other SOFA functions.
118///
119/// 8) This is one of several functions that inserts into the astrom
120/// structure star-independent parameters needed for the chain of
121/// astrometric transformations ICRS <-> GCRS <-> CIRS <-> observed.
122///
123/// The various functions support different classes of observer and
124/// portions of the transformation chain:
125/// ```
126/// functions observer transformation
127///
128/// iauApcg iauApcg13 geocentric ICRS <-> GCRS
129/// iauApci iauApci13 terrestrial ICRS <-> CIRS
130/// iauApco iauApco13 terrestrial ICRS <-> observed
131/// iauApcs iauApcs13 space ICRS <-> GCRS
132/// iauAper iauAper13 terrestrial update Earth rotation
133/// iauApio iauApio13 terrestrial CIRS <-> observed
134/// ```
135/// Those with names ending in "13" use contemporary SOFA models to
136/// compute the various ephemerides. The others accept ephemerides
137/// supplied by the caller.
138///
139/// The transformation from ICRS to GCRS covers space motion,
140/// parallax, light deflection, and aberration. From GCRS to CIRS
141/// comprises frame bias and precession-nutation. From CIRS to
142/// observed takes account of Earth rotation, polar motion, diurnal
143/// aberration and parallax (unless subsumed into the ICRS <-> GCRS
144/// transformation), and atmospheric refraction.
145///
146/// 9) The context structure astrom produced by this function is used by
147/// iauAtioq, iauAtoiq, iauAtciq* and iauAticq*.
148///
149/// Called:
150/// ```
151/// iauIr initialize r-matrix to identity
152/// iauRz rotate around Z-axis
153/// iauRy rotate around Y-axis
154/// iauRx rotate around X-axis
155/// iauAnpm normalize angle into range +/- pi
156/// iauC2ixys celestial-to-intermediate matrix, given X,Y and s
157/// iauPvtob position/velocity of terrestrial station
158/// iauTrxpv product of transpose of r-matrix and pv-vector
159/// iauApcs astrometry parameters, ICRS-GCRS, space observer
160/// iauCr copy r-matrix
161/// ```
162pub fn apco(
163 date1: f64,
164 date2: f64,
165 ebpv: &[[f64; 3]; 2],
166 ehp: &[f64; 3],
167 x: f64,
168 y: f64,
169 s: f64,
170 theta: f64,
171 elong: f64,
172 phi: f64,
173 hm: f64,
174 xp: f64,
175 yp: f64,
176 sp: f64,
177 refa: f64,
178 refb: f64,
179 astrom: &mut IauAstrom,
180) {
181 let r = &mut [[0.0; 3]; 3];
182 let (mut a, mut b, eral, c);
183
184 // Form the rotation matrix, CIRS to apparent [HA,Dec].
185 ir(r);
186 rz(theta + sp, r);
187 ry(-xp, r);
188 rx(-yp, r);
189 rz(elong, r);
190
191 // Solve for local Earth rotation angle.
192 a = r[0][0];
193 b = r[0][1];
194 eral = if a != 0.0 || b != 0.0 {
195 b.atan2(a)
196 } else {
197 0.0
198 };
199 astrom.eral = eral;
200
201 // Solve for polar motion [X,Y] with respect to local meridian.
202 a = r[0][0];
203 c = r[0][2];
204 astrom.xpl = c.atan2((a * a + b * b).sqrt());
205 a = r[1][2];
206 b = r[2][2];
207 astrom.ypl = if a != 0.0 || b != 0.0 {
208 -a.atan2(b)
209 } else {
210 0.0
211 };
212
213 // Adjusted longitude.
214 astrom.along = anpm(eral - theta);
215
216 // Functions of latitude.
217 astrom.sphi = phi.sin();
218 astrom.cphi = phi.cos();
219
220 // Refraction constants.
221 astrom.refa = refa;
222 astrom.refb = refb;
223
224 // Disable the (redundant) diurnal aberration step.
225 astrom.diurab = 0.0;
226
227 /* CIO based BPN matrix. */
228 c2ixys(x, y, s, r);
229
230 /* Observer's geocentric position and velocity (m, m/s, CIRS). */
231 let pvc = &mut [[0.0; 3]; 2];
232 pvtob(elong, phi, hm, xp, yp, sp, theta, pvc);
233
234 /* Rotate into GCRS. */
235 let pv = &mut [[0.0; 3]; 2];
236 trxpv(r, pvc, pv);
237
238 /* ICRS <-> GCRS parameters. */
239 apcs(date1, date2, pv, &ebpv, &ehp, astrom);
240
241 /* Store the CIO based BPN matrix. */
242 cr(r, &mut astrom.bpn);
243}