Sun position
Sun position algorithm
sky1.c
1 /*==============================================================================
2  * sky1.c - astronomical coordinate conversion routines, IAU 1980
3  *
4  * Author: David Hoadley
5  *
6  * Description: (see sky1.h)
7  *
8  * Copyright (c) 2020, David Hoadley <vcrumble@westnet.com.au>
9  * All rights reserved.
10  *
11  * Redistribution and use in source and binary forms, with or without
12  * modification, are permitted provided that the following conditions are met:
13  *
14  * * Redistributions of source code must retain the above copyright notice, this
15  * list of conditions and the following disclaimer.
16  * * Redistributions in binary form must reproduce the above copyright notice,
17  * this list of conditions and the following disclaimer in the documentation
18  * and/or other materials provided with the distribution.
19  *
20  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
21  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
22  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
23  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
24  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
25  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
26  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
27  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
28  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
29  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
30  * POSSIBILITY OF SUCH DAMAGE.
31  *
32  *==============================================================================
33  */
34 /*------------------------------------------------------------------------------
35  * Notes:
36  * Character set: UTF-8. (Non ASCII characters appear in this file)
37  *----------------------------------------------------------------------------*/
38 
39 /* ANSI includes etc. */
40 #include "instead-of-math.h"
41 #include <stdlib.h>
42 
43 /* Local and project includes */
44 #include "sky1.h"
45 
46 #include "astron.h"
47 #include "general.h"
48 
49 /*
50  * Local #defines and typedefs
51  */
52 DEFINE_THIS_FILE; // For use by REQUIRE() - assertions.
53 
54 /* Convert from units of 0.1 milliarcsec to radians */
55 #define MILLIARCSECx10_TO_RAD (PI / (180.0 * 3600.0 * 10000.0))
56 
57 #define NUM_TERMS 106
58 
59 /* Coefficients of fundamental arguments */
60 typedef struct {
61  int cd; // coefficient of d
62  int clp; // coefficient of l'
63  int cl; // coefficient of l
64  int cf; // coefficient of f
65  int com; // coefficient of omega
66 } Fcoeffs;
67 
68 /* Coefficients of dPsi and dEps accumulations */
69 typedef struct {
70  double ps1;
71  double ps2;
72  double ep1;
73  double ep2;
74 } PEcoeffs;
75 
76 /*
77  * Prototypes for local functions (not called from other modules)
78  */
79 
80 /*
81  * Global variables accessible by other modules
82  */
83 
84 
85 /*
86  * Local variables (not accessed by other modules)
87  */
88 /* Periodic Terms for the nutation in longitude and obliquity, sorted
89  roughly in descending order of PE coefficient magnitude.
90  (Terms in each of the following two tables MUST be in the same order) */
91 LOCAL const Fcoeffs coeffs[NUM_TERMS] =
92 { // d l' l F Om
93  { 0, 0, 0, 0, 1 }, /* 1 */
94  { -2, 0, 0, 2, 2 }, /* 9 */
95  { 0, 0, 0, 2, 2 }, /* 31 */
96  { 0, 0, 0, 0, 2 }, /* 2 */
97  { 0, 1, 0, 0, 0 }, /* 10 */
98  { 0, 0, 1, 0, 0 }, /* 32 */
99  { -2, 1, 0, 2, 2 }, /* 11 */
100  { 0, 0, 0, 2, 1 }, /* 33 */
101  { 0, 0, 1, 2, 2 }, /* 34 */
102  { -2, -1, 0, 2, 2 }, /* 12 */
103  { -2, 0, 1, 0, 0 }, /* 35 */
104  { -2, 0, 0, 2, 1 }, /* 13 */
105  { 0, 0, -1, 2, 2 }, /* 36 */
106  { 2, 0, 0, 0, 0 }, /* 37 */
107  { 0, 0, 1, 0, 1 }, /* 38 */
108  { 2, 0, -1, 2, 2 }, /* 40 */
109  { 0, 0, -1, 0, 1 }, /* 39 */
110  { 0, 0, 1, 2, 1 }, /* 41 */
111  { -2, 0, 2, 0, 0 }, /* 14 */
112  { 0, 0, -2, 2, 1 }, /* 3 */
113  { 2, 0, 0, 2, 2 }, /* 42 */
114  { 0, 0, 2, 2, 2 }, /* 45 */
115  { 0, 0, 2, 0, 0 }, /* 43 */
116  { -2, 0, 1, 2, 2 }, /* 44 */
117  { 0, 0, 0, 2, 0 }, /* 46 */
118  { -2, 0, 0, 2, 0 }, /* 15 */
119  { 0, 0, -1, 2, 1 }, /* 47 */
120  { 0, 2, 0, 0, 0 }, /* 16 */
121  { 2, 0, -1, 0, 1 }, /* 48 */
122  { -2, 2, 0, 2, 2 }, /* 18 */
123  { 0, 1, 0, 0, 1 }, /* 17 */
124  { -2, 0, 1, 0, 1 }, /* 49 */
125  { 0, -1, 0, 0, 1 }, /* 19 */
126  { 0, 0, 2, -2, 0 }, /* 4 */
127  { 2, 0, -1, 2, 1 }, /* 50 */
128 
129  { 2, 0, 1, 2, 2 }, /* 54 */
130  { -2, 1, 1, 0, 0 }, /* 51 */
131  { 0, 1, 0, 2, 2 }, /* 52 */
132  { 0, -1, 0, 2, 2 }, /* 53 */
133  { 2, 0, 0, 2, 1 }, /* 58 */
134 
135  { 2, 0, -2, 0, 1 }, /* 20 */
136  { 2, 0, 1, 0, 0 }, /* 55 */
137  { -2, 0, 2, 2, 2 }, /* 56 */
138  { 2, 0, 0, 0, 1 }, /* 57 */
139  { -2, 0, 1, 2, 1 }, /* 59 */
140 
141  { -2, -1, 0, 2, 1 }, /* 21 */
142  { -2, 0, 0, 0, 1 }, /* 60 */
143  { 0, -1, 1, 0, 0 }, /* 61 */
144  { 0, 0, 2, 2, 1 }, /* 62 */
145 
146  { -2, 0, 2, 0, 1 }, /* 22 */
147  { -2, 1, 0, 2, 1 }, /* 23 */
148  { -1, 0, 1, 0, 0 }, /* 24 */
149  { -2, 1, 0, 0, 0 }, /* 63 */
150  { 0, 0, 1, -2, 0 }, /* 64 */
151  { 1, 0, 0, 0, 0 }, /* 65 */
152 
153  { 0, 0, -2, 2, 2 }, /* 5 */
154  { -1, -1, 1, 0, 0 }, /* 6 */
155  { 0, 1, 1, 0, 0 }, /* 66 */
156  { 0, 0, 1, 2, 0 }, /* 67 */
157  { 0, -1, 1, 2, 2 }, /* 68 */
158  { 2, -1, -1, 2, 2 }, /* 69 */
159  { 0, 0, 3, 2, 2 }, /* 71 */
160  { 2, -1, 0, 2, 2 }, /* 72 */
161 
162  { -2, -2, 0, 2, 1 }, /* 7 */
163  { 0, 0, -2, 0, 1 }, /* 70 */
164  { 0, 1, 1, 2, 2 }, /* 73 */
165  { -2, 0, -1, 2, 1 }, /* 74 */
166  { 0, 0, 2, 0, 1 }, /* 75 */
167  { 0, 0, 1, 0, 2 }, /* 76 */
168  { 0, 0, 3, 0, 0, }, /* 77 */
169  { 1, 0, 0, 2, 2, }, /* 78 */
170  { 4, 0, -1, 2, 2 }, /* 82 */
171 
172  { 0, 0, 2, -2, 1 }, /* 8 */
173  { -2, 1, 2, 0, 0 }, /* 25 */
174  { 2, 0, 0, -2, 1 }, /* 26 */
175  { 2, 1, 0, -2, 0 }, /* 27 */
176  { 0, 1, 0, 0, 2 }, /* 28 */
177  { 1, 0, -1, 0, 1 }, /* 29 */
178  { -2, 1, 0, 2, 0 }, /* 30 */
179 
180  { 0, 0, -1, 0, 2 }, /* 79 */
181  { -4, 0, 1, 0, 0 }, /* 80 */
182  { 2, 0, -2, 2, 2 }, /* 81 */
183  { -4, 0, 2, 0, 0 }, /* 83 */
184  { -2, 1, 1, 2, 2 }, /* 84 */
185  { 2, 0, 1, 2, 1 }, /* 85 */
186  { 4, 0, -2, 2, 2 }, /* 86 */
187  { 0, 0, -1, 4, 2 }, /* 87 */
188  { -2, -1, 1, 0, 0 }, /* 88 */
189  { -2, 0, 2, 2, 1 }, /* 89 */
190  { 2, 0, 2, 2, 2, }, /* 90 */
191  { 2, 0, 1, 0, 1 }, /* 91 */
192  { -2, 0, 0, 4, 2 }, /* 92 */
193  { -2, 0, 3, 2, 2 }, /* 93 */
194  { -2, 0, 1, 2, 0 }, /* 94 */
195  { 0, 1, 0, 2, 1 }, /* 95 */
196  { 2, -1, -1, 0, 1 }, /* 96 */
197  { 0, 0, 0, -2, 1 }, /* 97 */
198  { -1, 0, 0, 2, 2 }, /* 98 */
199  { 2, 1, 0, 0, 0 }, /* 99 */
200  { -2, 0, 1, -2, 0 }, /* 100 */
201  { 0, -1, 0, 2, 1 }, /* 101 */
202  { -2, 1, 1, 0, 1 }, /* 102 */
203  { 2, 0, 1, -2, 0 }, /* 103 */
204  { 2, 0, 2, 0, 0 }, /* 104 */
205  { 4, 0, 0, 2, 2 }, /* 105 */
206  { 1, 1, 0, 0, 0 }, /* 106 */
207 };
208 
209 LOCAL const PEcoeffs pec[NUM_TERMS] = {
210  {-171996.0, -174.2, 92025.0, 8.9 }, /* 1 */
211  {-13187.0, -1.6, 5736.0, -3.1 }, /* 9 */
212  {-2274.0, -0.2, 977.0, -0.5 }, /* 31 */
213  { 2062.0, 0.2, -895.0, 0.5 }, /* 2 */
214  { 1426.0, -3.4, 54.0, -0.1 }, /* 10 */
215  { 712.0, 0.1, -7.0, 0.0 }, /* 32 */
216  { -517.0, 1.2, 224.0, -0.6 }, /* 11 */
217  { -386.0, -0.4, 200.0, 0.0 }, /* 33 */
218  { -301.0, 0.0, 129.0, -0.1 }, /* 34 */
219  { 217.0, -0.5, -95.0, 0.3 }, /* 12 */
220  { -158.0, 0.0, -1.0, 0.0 }, /* 35 */
221  { 129.0, 0.1, -70.0, 0.0 }, /* 13 */
222  { 123.0, 0.0, -53.0, 0.0 }, /* 36 */
223  { 63.0, 0.0, -2.0, 0.0 }, /* 37 */
224  { 63.0, 0.1, -33.0, 0.0 }, /* 38 */
225  { -59.0, 0.0, 26.0, 0.0 }, /* 40 */
226  { -58.0, -0.1, 32.0, 0.0 }, /* 39 */
227  { -51.0, 0.0, 27.0, 0.0 }, /* 41 */
228  { 48.0, 0.0, 1.0, 0.0 }, /* 14 */
229  { 46.0, 0.0, -24.0, 0.0 }, /* 3 */
230  { -38.0, 0.0, 16.0, 0.0 }, /* 42 */
231  { -31.0, 0.0, 13.0, 0.0 }, /* 45 */
232  { 29.0, 0.0, -1.0, 0.0 }, /* 43 */
233  { 29.0, 0.0, -12.0, 0.0 }, /* 44 */
234  { 26.0, 0.0, -1.0, 0.0 }, /* 46 */
235  { -22.0, 0.0, 0.0, 0.0 }, /* 15 */
236  { 21.0, 0.0, -10.0, 0.0 }, /* 47 */
237  { 17.0, -0.1, 0.0, 0.0 }, /* 16 */
238  { 16.0, 0.0, -8.0, 0.0 }, /* 48 */
239  { -16.0, 0.1, 7.0, 0.0 }, /* 18 */
240  { -15.0, 0.0, 9.0, 0.0 }, /* 17 */
241  { -13.0, 0.0, 7.0, 0.0 }, /* 49 */
242  { -12.0, 0.0, 6.0, 0.0 }, /* 19 */
243  { 11.0, 0.0, 0.0, 0.0 }, /* 4 */
244  { -10.0, 0.0, 5.0, 0.0 }, /* 50 */
245  /* 35 terms to here */
246  { -8.0, 0.0, 3.0, 0.0 }, /* 54 */
247  { -7.0, 0.0, 0.0, 0.0 }, /* 51 */
248  { 7.0, 0.0, -3.0, 0.0 }, /* 52 */
249  { -7.0, 0.0, 3.0, 0.0 }, /* 53 */
250  { -7.0, 0.0, 3.0, 0.0 }, /* 58 */
251  /* 40 terms to here */
252  { -6.0, 0.0, 3.0, 0.0 }, /* 20 */
253  { 6.0, 0.0, 0.0, 0.0 }, /* 55 */
254  { 6.0, 0.0, -3.0, 0.0 }, /* 56 */
255  { -6.0, 0.0, 3.0, 0.0 }, /* 57 */
256  { 6.0, 0.0, -3.0, 0.0 }, /* 59 */
257  /* 45 terms to here */
258  { -5.0, 0.0, 3.0, 0.0 }, /* 21 */
259  { -5.0, 0.0, 3.0, 0.0 }, /* 60 */
260  { 5.0, 0.0, 0.0, 0.0 }, /* 61 */
261  { -5.0, 0.0, 3.0, 0.0 }, /* 62 */
262  /* 49 terms to here */
263  { 4.0, 0.0, -2.0, 0.0 }, /* 22 */
264  { 4.0, 0.0, -2.0, 0.0 }, /* 23 */
265  { -4.0, 0.0, 0.0, 0.0 }, /* 24 */
266  { -4.0, 0.0, 0.0, 0.0 }, /* 63 */
267  { 4.0, 0.0, 0.0, 0.0 }, /* 64 */
268  { -4.0, 0.0, 0.0, 0.0 }, /* 65 */
269  /* 55 terms to here */
270  { -3.0, 0.0, 1.0, 0.0 }, /* 5 */
271  { -3.0, 0.0, 0.0, 0.0 }, /* 6 */
272  { -3.0, 0.0, 0.0, 0.0 }, /* 66 */
273  { 3.0, 0.0, 0.0, 0.0 }, /* 67 */
274  { -3.0, 0.0, 1.0, 0.0 }, /* 68 */
275  { -3.0, 0.0, 1.0, 0.0 }, /* 69 */
276  { -3.0, 0.0, 1.0, 0.0 }, /* 71 */
277  { -3.0, 0.0, 1.0, 0.0 }, /* 72 */
278  /* 63 terms to here */
279  { -2.0, 0.0, 1.0, 0.0 }, /* 7 */
280  { -2.0, 0.0, 1.0, 0.0 }, /* 70 */
281  { 2.0, 0.0, -1.0, 0.0 }, /* 73 */
282  { -2.0, 0.0, 1.0, 0.0 }, /* 74 */
283  { 2.0, 0.0, -1.0, 0.0 }, /* 75 */
284  { -2.0, 0.0, 1.0, 0.0 }, /* 76 */
285  { 2.0, 0.0, 0.0, 0.0 }, /* 77 */
286  { 2.0, 0.0, -1.0, 0.0 }, /* 78 */
287  { -2.0, 0.0, 1.0, 0.0 }, /* 82 */
288  /* 72 terms to here */
289  { 1.0, 0.0, 0.0, 0.0 }, /* 8 */
290  { 1.0, 0.0, 0.0, 0.0 }, /* 25 */
291  { 1.0, 0.0, 0.0, 0.0 }, /* 26 */
292  { -1.0, 0.0, 0.0, 0.0 }, /* 27 */
293  { 1.0, 0.0, 0.0, 0.0 }, /* 28 */
294  { 1.0, 0.0, 0.0, 0.0 }, /* 29 */
295  { -1.0, 0.0, 0.0, 0.0 }, /* 30 */
296 
297  { 1.0, 0.0, -1.0, 0.0 }, /* 79 */
298  { -1.0, 0.0, 0.0, 0.0 }, /* 80 */
299  { 1.0, 0.0, -1.0, 0.0 }, /* 81 */
300  { -1.0, 0.0, 0.0, 0.0 }, /* 83 */
301  { 1.0, 0.0, -1.0, 0.0 }, /* 84 */
302  { -1.0, 0.0, 1.0, 0.0 }, /* 85 */
303  { -1.0, 0.0, 1.0, 0.0 }, /* 86 */
304  { 1.0, 0.0, 0.0, 0.0 }, /* 87 */
305  { 1.0, 0.0, 0.0, 0.0 }, /* 88 */
306  { 1.0, 0.0, -1.0, 0.0 }, /* 89 */
307  { -1.0, 0.0, 0.0, 0.0 }, /* 90 */
308  { -1.0, 0.0, 0.0, 0.0 }, /* 91 */
309  { 1.0, 0.0, 0.0, 0.0 }, /* 92 */
310  { 1.0, 0.0, 0.0, 0.0 }, /* 93 */
311  { -1.0, 0.0, 0.0, 0.0 }, /* 94 */
312  { 1.0, 0.0, 0.0, 0.0 }, /* 95 */
313  { 1.0, 0.0, 0.0, 0.0 }, /* 96 */
314  { -1.0, 0.0, 0.0, 0.0 }, /* 97 */
315  { -1.0, 0.0, 0.0, 0.0 }, /* 98 */
316  { -1.0, 0.0, 0.0, 0.0 }, /* 99 */
317  { -1.0, 0.0, 0.0, 0.0 }, /* 100 */
318  { -1.0, 0.0, 0.0, 0.0 }, /* 101 */
319  { -1.0, 0.0, 0.0, 0.0 }, /* 102 */
320  { -1.0, 0.0, 0.0, 0.0 }, /* 103 */
321  { 1.0, 0.0, 0.0, 0.0 }, /* 104 */
322  { -1.0, 0.0, 0.0, 0.0 }, /* 105 */
323  { 1.0, 0.0, 0.0, 0.0 } /* 106 */
324 };
325 
326 LOCAL const int nTerms[] = { NUM_TERMS, 72, 63, 49, 35 };
327 
328 /*
329  *==============================================================================
330  *
331  * Implementation
332  *
333  *==============================================================================
334  *
335  * Global functions callable by other modules
336  *
337  *------------------------------------------------------------------------------
338  */
339 GLOBAL void sky1_frameBiasFK5(V3D_Matrix *biasM)
340 /*! Create the frame bias matrix from the IAU 2000 precession-nutation model.
341  This matrix is used to convert coordinates from the Geocentric Celestial
342  Reference System (GCRS) to FK5 reference system. (Although anything to do
343  with the ICRS/GCRS is related to the IAU2000+ precession-nutation theory,
344  this routine is included here to allow ICRS star catalogue positions to be
345  converted to apparent coordinates using the IAU1980 precession-nutation
346  theory, which is quite good enough for tracking. It is not good enough for
347  astrometry, but it will get us to better than within an arcsecond of our
348  desired apparent position).
349  \param[out] biasM Frame bias matrix \b B
350 
351  \par Reference:
352  _Astronomical Almanac_ 2007, page B28
353 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
354 {
355  const double eta0_as = -19.9 / 1000.0;
356  const double xi0_as = 9.1 / 1000.0;
357  const double dAlpha0_as = -22.9 / 1000.0;
358  V3D_Matrix cM, bM, aM;
359 
360  REQUIRE_NOT_NULL(biasM);
361 
362  /* First intermediate result. [cM] = [R2(xi0)] x [R3(dAlpha0)] */
363  v3d_createRotationMatrix(&aM, Zaxis, arcsecToRad(dAlpha0_as));
364  v3d_createRotationMatrix(&bM, Yaxis, arcsecToRad(xi0_as));
365  v3d_multMxM(&cM, &bM, &aM);
366 
367  /* Re-use aM for rotation by eta0, obtain
368  * [biasM] = [R1(eta0)] x [R2(xi0)] x [R3(dAlpha0)] */
369  v3d_createRotationMatrix(&aM, Xaxis, arcsecToRad(-eta0_as));
370  v3d_multMxM(biasM, &aM, &cM);
371 }
372 
373 
374 
375 GLOBAL void sky1_precessionIAU1976(double t0_cy,
376  double t1_cy,
377  Sky1_Prec1976 *terms)
378 /*! This procedure calculates the equatorial precession parameters ζ, z, and θ
379  which represent the rotation required to transform the FK5 equatorial
380  reference system of Julian epoch t0 to that of Julian epoch t1, according to
381  the IAU 1976 precession constants.
382  \param[in] t0_cy Julian centuries since J2000.0 of initial epoch, TT
383  timescale
384  \param[in] t1_cy Julian centuries since J2000.0 of final epoch, TT timescale
385  \param[out] terms The three precession angles ζ, z, θ.
386 
387  \par References:
388  Lieske,J.H., _Astron. Astrophys_, Vol 73, pp282-284, 1979\n
389  _Supplement to The Astronomical Almanac 1984_, HMNAO, ppS18-19
390 
391  \par When to call this function
392  It is quite likely that you will not need to call this function directly.
393  If you are tracking a celestial object, using star_getApparent() or
394  star_getTopocentric() or calling planet_getApparent() or
395  planet_getTopocentric(), they will call this routine for you.
396  Likewise, if you are tracking the object using the skyfast module, the call
397  to skyfast_init() will call either star_getApparent() or
398  planet_getApparent(), and therefore call this routine for you.
399  \par
400  The values calculated by this routine change only slowly. So if you are
401  calling it yourself, you can call it infrequently. Intervals of up to an
402  hour between calls will not introduce much error. This function is also
403  called by sky1_createNPmatrix().
404 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
405 {
406  double t0Sq;
407  double t, tSq, tCu;
408  double temp_as; // angle (arcseconds)
409 
410  REQUIRE_NOT_NULL(terms);
411 
412  /* Calculate interval between initial epoch T0 and final epoch T1 in Julian
413  centuries */
414  t = (t1_cy - t0_cy);
415 
416  /* Calculate squares and cubes of intervals */
417  tSq = t * t;
418  tCu = tSq * t;
419  t0Sq = t0_cy * t0_cy;
420 
421  /* Calculate equatorial precession parameters */
422  temp_as = (2306.2181 + 1.39656 * t0_cy - 0.000139 * t0Sq) * t
423  + (0.30188 - 0.000344 * t0_cy) * tSq + 0.017998 * tCu;
424  terms->zeta_rad = arcsecToRad(temp_as);
425 
426  temp_as = (2306.2181 + 1.39656 * t0_cy - 0.000139 * t0Sq) * t
427  + (1.09468 + 0.000066 * t0_cy) * tSq + 0.018203 * tCu;
428  terms->zed_rad = arcsecToRad(temp_as);
429 
430  temp_as = (2004.3109 - 0.85330 * t0_cy - 0.000217 * t0Sq) * t
431  - (0.42665 + 0.000217 * t0_cy) * tSq - 0.041833 * tCu;
432  terms->theta_rad = arcsecToRad(temp_as);
433 }
434 
435 
436 
437 GLOBAL void sky1_createPrec1976Matrix(const Sky1_Prec1976 *terms,
438  V3D_Matrix *precM)
439 /*! This routine calculates the precession matrix, based on angles ζ, z, and θ.
440  \a precM is the combined orthogonal rotation matrix \b P, required for
441  rigorous precession transformations using rectangular coordinates and matrix
442  methods:
443  V1 = P * V0
444  It appears to be calculated from\n
445  \b P = R3(-z) × R2(θ) × R3(-ζ)
446  \param[in] terms The three precession angles ζ, z, θ, as returned by
447  sky1_precessionIAU1976()
448  \param[out] precM Precession rotation matrix
449 
450  \par References:
451  Lieske,J.H., _Astron. Astrophys_, Vol 73, pp282-284, 1979\n
452  _Supplement to The Astronomical Almanac 1984_, HMNAO, ppS18-19
453 
454  \par When to call this function
455  Call this function after new values of ζ, z and θ have been calculated
456  by routine sky1_precessionIAU1976(). So, again, it need only be called
457  infrequently. But as mentioned for that function, it is quite likely that
458  you will not need to call this function directly. Any of the functions
459  sky1_createNPmatrix(), star_getApparent(), star_getTopocentric(),
460  planet_getApparent() or planet_getTopocentric() will call this
461  routine for you.
462 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
463 {
464  V3D_Matrix zedM;
465  V3D_Matrix thetaM;
466  V3D_Matrix zetaM;
467  V3D_Matrix tempM;
468 
469  REQUIRE_NOT_NULL(terms);
470  REQUIRE_NOT_NULL(precM);
471 
472  /* Calculate the precession matrix as multiplication of rotation matrices */
473  v3d_createRotationMatrix(&zetaM, Zaxis, -terms->zeta_rad);
474  v3d_createRotationMatrix(&thetaM, Yaxis, terms->theta_rad);
475  v3d_createRotationMatrix(&zedM, Zaxis, -terms->zed_rad);
476  v3d_multMxM(precM, &zedM, v3d_multMxM(&tempM, &thetaM, &zetaM));
477 }
478 
479 
480 
481 GLOBAL void sky1_nutationIAU1980(double t_cy,
482  int precision,
483  Sky1_Nut1980 *nut)
484 /*! Calculates the nutation in longitude and obliquity, according to the IAU
485  1980 Nutation Theory. Calculates first the fundamental nutation arguments,
486  and then a series of terms. There are 106 terms in the series, but fewer
487  may be selected for faster execution.
488  \param[in] t_cy Julian centuries since J2000.0, TT timescale
489  \param[in] precision How much precision do you want?
490  Valid range [0, 4]. Values outside this range will be
491  clamped to the range.
492  - 0 = full precision, use full 106-term series
493  - 1 = ignore terms < 0.2 milliarcseconds. 72-term series
494  - 2 = ignore terms < 0.3 milliarcseconds. 63-term series
495  - 3 = ignore terms < 0.5 milliarcseconds. 49-term series
496  - 4 = ignore terms < 1.0 milliarcseconds. 35-term series
497  \param[out] nut field \a nut->dPsi_rad - Nutation in longitude Δψ (radian)\n
498  field \a nut->dEps_rad - Nutation in obliquity Δε (radian)
499 
500  \par References:
501  Final report of the IAU Working Group on Nutation,
502  chairman P.K.Seidelmann, 1980,
503  published in _Celestial Mechanics_, Vol 27, pp79-106, 1982.\n
504  Kaplan,G.H. _USNO circular no. 163_, pp A3-6, 1981.\n
505  _Supplement to the Astronomical Almanac_ 1984.\n
506  The form of the data tables was inspired by, and extended from,
507  Reda, I. and Andreas, A. (2003), "Solar Position Algorithm
508  for Solar Radiation Applications", National Renewable Energy Laboratory,
509  NREL publication no. NREL/TP-560-34302.
510 
511  \note
512  - With precision = 0, exact agreement with tabulated values in the
513  _Astronomical Almanac_ 1987 pages B24 - B31 was found (to 3 decimal places,
514  as tabulated in the Almanac). Likewise for the 1990 edition pages B24 - B31.
515  But this IAU 1980 theory has been superseded by the more accurate IAU 2000
516  theory. This is evident when comparing outputs with the tabulated values in
517  the 2007 Almanac, pages B34 - B41 (to 4 decimal places). There are
518  differences of several milliarcseconds.
519 
520  - With precision = 1, differences up to 1.0 milliarcsecond were seen compared
521  to the precision = 0 values during testing.
522  - With precision = 2, differences up to 1.2 milliarcseconds were seen
523  - With precision = 3, differences up to 1.5 milliarcseconds were seen
524  - With precision = 4, differences of up to 5 milliarcseconds were seen.
525 
526  For tracking purposes, these are still very small.
527  The testing was not exhaustive, so take this as a rough guide only.
528 
529  \par When to call this function
530  It is quite likely that you will not need to call this function directly.
531  If you are tracking a celestial object, using star_getApparent() or
532  star_getTopocentric() or calling planet_getApparent() or
533  planet_getTopocentric(), they will call this routine for you.
534  Likewise, if you are tracking the object using the skyfast module, the call
535  to skyfast_init() will call either star_getApparent() or
536  planet_getApparent(), and therefore call this routine for you.
537  \par
538  The values calculated by this routine change only slowly. So if you are
539  calling it yourself, you can call it infrequently. Intervals of up to an
540  hour between calls will not introduce much error. This function is also
541  called by sky1_createNPmatrix().
542 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
543 {
544  const double Revs2Arcsec = 1296000.0; // Revolutions to arc seconds
545 
546  /* Fundamental Nutation arguments at date */
547  double l; // l - mean anomaly of the Moon (radian) (usually called M?)
548  double lp; // l' - mean anomaly of the Sun (Earth) (radian) (M'?)
549  double om; // Ω - longitude of the ascending node of the Moon's mean orbit
550  // on the ecliptic, measured from the mean equinox of date
551  // (radian)
552  double d; // D - mean elongation of the Moon from the Sun (radian)
553  double f; // F - Mean longitude of the Moon (L) minus mean longitude of
554  // the Moon's node (radian), = L - Ω
555 
556  double psiSum_masx10; // Nutation in longitude (units - 0.1 mas)
557  double epsSum_masx10; // Nutation in obliquity (units - 0.1 mas)
558  int i;
559  double a_rad; // angle - summation of args (radian)
560 
561  REQUIRE_NOT_NULL(nut);
562 
563  if (precision < 0) { precision = 0; }
564  if (precision >= ARRAY_SIZE(nTerms)) { precision = ARRAY_SIZE(nTerms) - 1; }
565 
566  // Calculate FUNDAMENTAL ARGUMENTS in the FK5 reference system
567 
568  // Mean elongation of the Moon from the Sun (called X0 in NREL SPA)
569  d = arcsecToRad(1072261.307 + (1236.0 * Revs2Arcsec + 1105601.328
570  + (-6.891 + 0.019 * t_cy) * t_cy) * t_cy);
571  //+
572  d = normalize(d, TWOPI);
573  //-
574 
575  // Solar Mean Anomaly = Mean longitude of the Sun minus mean longitude of
576  // the Sun's perigee (called X1 in NREL SPA)
577  lp = arcsecToRad(1287099.804 + (99.0 * Revs2Arcsec + 1292581.224
578  + (-0.577 - 0.012 * t_cy) * t_cy) * t_cy);
579  //+
580  lp = normalize(lp, TWOPI);
581  //-
582 
583  // Lunar Mean Anomaly = Mean longitude of the Moon minus mean longitude of
584  // the Moon 's perigee (called X2 in NREL SPA)
585  l = arcsecToRad(485866.733 + (1325.0 * Revs2Arcsec + 715922.633
586  + (31.31 + 0.064 * t_cy) * t_cy) * t_cy);
587  //+
588  l = normalize(l, TWOPI);
589  //-
590 
591  // Mean longitude of the Moon minus mean longitude of the Moon's node
592  // (called X3 in NREL SPA)
593  f = arcsecToRad(335778.877 + (1342.0 * Revs2Arcsec + 295263.137
594  + (-13.257 + 0.011 * t_cy) * t_cy) * t_cy);
595  //+
596  f = normalize(f, TWOPI);
597  //-
598 
599 
600  // Longitude of the mean ascending node of the lunar orbit on the
601  // ecliptic, measured from the mean equinox of date (X4 in NREL SPA)
602  om = arcsecToRad(450160.28 + (-5.0 * Revs2Arcsec - 482890.539
603  + (7.455 + 0.008 * t_cy) * t_cy) * t_cy);
604  //+
605  //om = normalize(om, TWOPI) - TWOPI;
606  //-
607 
608 
609  // Multiply through the table of nutation co-efficients and add up all the
610  // terms.
611  psiSum_masx10 = 0.0;
612  epsSum_masx10 = 0.0;
613  for (i = nTerms[precision] - 1; i >= 0; i--) {
614  a_rad = d * coeffs[i].cd + lp * coeffs[i].clp + l * coeffs[i].cl
615  + f * coeffs[i].cf + om * coeffs[i].com;
616  psiSum_masx10 += sin(a_rad) * (pec[i].ps1 + t_cy * pec[i].ps2);
617  // No need to calculate cos if its coefficients are zero. And if ep1
618  // is zero, so is ep2, so we don't need to test it.
619  if (pec[i].ep1 != 0.0) {
620  epsSum_masx10 += cos(a_rad) * (pec[i].ep1 + t_cy * pec[i].ep2);
621  }
622  }
623 
624  // Convert results to radians
625  nut->dPsi_rad = psiSum_masx10 * MILLIARCSECx10_TO_RAD;
626  nut->dEps_rad = epsSum_masx10 * MILLIARCSECx10_TO_RAD;
627 }
628 
629 
630 
631 GLOBAL void sky1_epsilon1980(double t_cy, Sky1_Nut1980 *nut)
632 /*! Calculate the obliquity of the ecliptic and the equation of the equinoxes
633  \param[in] t_cy Julian centuries since J2000.0, TT timescale
634  \param[in,out] nut [in] field \a nut->dPsi_rad - Nutation in longitude Δψ,
635  as returned by function sky1_nutationIAU1980()
636  (radian)\n
637  [in] field \a nut->dEps_rad - Nutation in obliquity Δε,
638  as returned by function sky1_nutationIAU1980()
639  (radian)\n
640  [out] field \a nut->eps0_rad - Mean obliquity of the
641  ecliptic ε0 (radian)\n
642  [out] field \a nut->eqEq_rad - Equation of the equinoxes
643  = Δψ * cos(ε0 + Δε) (radian) Note: not seconds
644  \par References
645  _Astronomical Almanac_ 1987 page B18\n
646  International Astronomical Union, Standards of Fundamental Astronomy
647  software collection, release 2017-04-20. http://www.iausofa.org
648  routine \c iauObl80()
649 
650  \par When to call this function
651  Call this function after new values of Δψ, and Δε have been calculated
652  by routine sky1_nutationIAU1980(). So, again, it need only be called
653  infrequently. But as mentioned for that function, it is quite likely that
654  you will not need to call this function directly. Any of the functions
655  sky1_createNPmatrix(), star_getApparent(), star_getTopocentric(),
656  planet_getApparent() or planet_getTopocentric() will call this
657  routine for you.
658 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
659 {
660  double eps0_as;
661 
662  REQUIRE_NOT_NULL(nut);
663 
664  // Calculate Mean obliquity in radian
665  eps0_as = 84381.448
666  + (-46.8150 + (-0.00059 + 0.001813 * t_cy) * t_cy) * t_cy;
667  nut->eps0_rad = arcsecToRad(eps0_as);
668 
669  // Calculate the Equation of the Equinoxes in radian (not in seconds, which
670  // is more usually done)
671  nut->eqEq_rad = nut->dPsi_rad * cos(nut->eps0_rad + nut->dEps_rad);
672 }
673 
674 
675 
676 GLOBAL void sky1_createNut1980Matrix(const Sky1_Nut1980 *nut,
677  V3D_Matrix *nutM)
678 /*! This routine calculates the Nutation matrix, using the nutation in
679  longitude, the nutation in obliquity, and the mean obliquity of the equator.
680  We calculate the Nutation matrix in a completely rigorous manner from the
681  product of three rotation matrices:\n
682  \b N = R1(-ε) × R3(-Δψ) × R1(ε0)
683  where ε = ε0 + Δε (true obliquity of the equator)
684  \param[in] nut Angles Δψ, Δε and ε0, as returned by functions
685  sky1_nutationIAU1980() and sky1_epsilon1980()
686  \param[out] nutM Nutation rotation matrix \b N
687 
688  \par References:
689  _Astronomical Almanac_ 2007, page B31\n
690  also
691  Mueller, p75
692 
693  \par When to call this function
694  Call it after you have called sky1_nutationIAU1980() and sky1_epsilon1980().
695  As mentioned for the other functions above, this need only be done
696  infrequently. But also as mentioned, it is quite likely that
697  you will not need to call this function directly. Any of the functions
698  sky1_createNPmatrix(), star_getApparent(), star_getTopocentric(),
699  planet_getApparent() or planet_getTopocentric() will call this
700  routine for you.
701 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
702 {
703  V3D_Matrix eps0M;
704  V3D_Matrix psiM;
705  V3D_Matrix epsM;
706  V3D_Matrix tempM;
707 
708  REQUIRE_NOT_NULL(nut);
709  REQUIRE_NOT_NULL(nutM);
710 
711  v3d_createRotationMatrix(&eps0M, Xaxis, nut->eps0_rad);
712  v3d_createRotationMatrix(&psiM, Zaxis, -nut->dPsi_rad);
713  v3d_createRotationMatrix(&epsM, Xaxis, -(nut->eps0_rad + nut->dEps_rad));
714 
715  // Multiply the three matrices in the order epsM * psiM * eps0M
716  v3d_multMxM(nutM, &epsM, v3d_multMxM(&tempM, &psiM, &eps0M));
717 }
718 
719 
720 
721 GLOBAL void sky1_createNPmatrix(double t0_cy,
722  double t1_cy,
723  int precision,
724  V3D_Matrix* npM)
725 /*! Call the various precession and nutation routines in this module to create
726  a combined precession and nutation rotation matrix, suitable for converting
727  a celestial object's coordinates from mean place at epoch \a t0_cy to
728  apparent coordinates at epoch \a t1_cy.
729  \param[in] t0_cy Julian centuries since J2000.0 of initial epoch,
730  TT timescale
731  \param[in] t1_cy Julian centuries since J2000.0 of final epoch,
732  TT timescale
733  \param[in] precision Precision specification to be passed to routine
734  sky1_nutationIAU1980(). See that routine for
735  explanation
736  \param[out] npM Combined precession and nutation matrix
737  \b NP = \b N × \b P
738 
739  This function calls routines sky1_precessionIAU1976(),
740  sky1_createPrec1976Matrix(), sky1_nutationIAU1980(), sky1_epsilon1980()
741  and sky1_createNut1980Matrix() to calculate the precession and nutation
742  matrices.
743 
744  \par When to call this function
745  Call this function if you are tracking non-solar system objects.
746  The resulting matrix changes only very slowly, so this function needs to be
747  called only once per hour.
748 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
749 {
750  V3D_Matrix nM, pM;
751  Sky1_Prec1976 prec;
752  Sky1_Nut1980 nut;
753 
754  REQUIRE_NOT_NULL(npM);
755 
756  sky1_precessionIAU1976(t0_cy, t1_cy, &prec);
757  sky1_createPrec1976Matrix(&prec, &pM);
758  sky1_nutationIAU1980(t1_cy, precision, &nut);
759  sky1_epsilon1980(t1_cy, &nut);
760  sky1_createNut1980Matrix(&nut, &nM);
761  v3d_multMxM(npM, &nM, &pM);
762 }
763 
764 
765 
766 GLOBAL double sky1_gmSiderealTimeIAU1982(double du)
767 /*! Calculate the Greenwich mean sidereal time.
768  \returns Greenwich Mean Sidereal Time (radian)
769  \param[in] du Days since J2000.0, UT1 timescale
770 
771  \par Reference:
772  _Astronomical Almanac_ 1987, page B6
773 
774  \par When to call this function
775  If you are tracking a celestial object, you need not call this function
776  directly. So long as you call sky1_appToTirs() every time around your
777  control loop, it will call this function for you.
778 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
779 {
780  /*
781  * The expression given on page B6 of the Almanac, which is
782  * GMST = 24110.54841 + 8640184.812866 Tu + 0.093104 Tu^2 - 6.2e-6 Tu^3
783  * is only defined for 0 h UT1 (midnight), meaning that it is only defined
784  * for Du = 0.5, 1.5, 2.5 etc., (since Du is days since noon, 1st January
785  * 2000). So we need to add a term which defines it for the rest of the day,
786  * which means adding 2Pi radian per day. But this new term will add an
787  * offset of Pi at every midnight, so we need to take that out again.
788  */
789 
790  // Time series coefficients for Greenwich Mean Sidereal Time
791  const double B0 = secToRad(24110.54841) - PI;
792  const double B1 = secToRad(8640184.812866);
793  const double B2 = secToRad(0.093104);
794  const double B3 = secToRad(-6.2e-6);
795 
796  double gmst_rad; // Greenwich Mean Sidereal Time (radian)
797  double tu; // Julian centuries since J2000.0, UT1 timescale
798  double dayfrac; // fractional part of du
799 
800  tu = du / JUL_CENT;
801  dayfrac = du - floor(du);
802 
803  gmst_rad = B0 + (tu * (B1 + (tu * (B2 + (tu * B3))))) + dayfrac * TWOPI;
804 
805  return normalize(gmst_rad, TWOPI);
806 }
807 
808 
809 
810 GLOBAL void sky1_appToTirs(const V3D_Vector *appV,
811  double j2kUT1_d,
812  double eqEq_rad,
813  V3D_Vector *terInterV)
814 /*! Convert a position in geocentric apparent coordinates to geocentric
815  coordinates in the Terrestrial Intermediate Reference System. This is the
816  first stage of converting apparent coordinates to topocentric coordinates.
817  The resulting vector depends upon the current rotational position of the
818  Earth. (For the second stage, to obtain topocentric coordinates, call
819  routine sky_siteTirsToTopo()).
820  \param[in] appV Position vector of apparent place
821  (unit vector in equatorial coordinates)
822  \param[in] j2kUT1_d days since J2000.0, UT1 timescale, as returned by function
823  sky_updateTimes() in the \a j2kUT1_d field of the
824  Sky_Times struct.
825  \param[in] eqEq_rad Equation of the equinoxes (radian), as returned by function
826  sky1_epsilon1980() in the \a eqEq_rad field of the
827  Sky1_Nut1980 struct.
828 
829  \param[out] terInterV Position vector in Terrestrial Intermediate Ref System
830 
831  \par When to call this function
832  When you have the position of a celestial object expressed in Apparent
833  coordinates, use this function to convert it to Terrestrial coordinates at
834  the relevant rotational position of the earth. (Effectively you are
835  converting from Right Ascension and Declination to the negative of the
836  Greenwich Hour Angle and Declination.)
837  If you are running a control loop to enable continuous tracking
838  of this object, you will need to call this function (once) every time around
839  your control loop.
840  \par
841  Follow this function with a call to sky_siteTirsToTopo() to obtain the
842  object's position in topocentric coordinates at the observing site.
843 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
844 {
845  double gast_rad; // Greenwich Apparent Sidereal Time (GAST)
846  V3D_Matrix earthRotM; // rotation matrix for current GAST
847 
848  REQUIRE_NOT_NULL(appV);
849  REQUIRE_NOT_NULL(terInterV);
850 
851  /* Get sidereal times */
852  gast_rad = sky1_gmSiderealTimeIAU1982(j2kUT1_d) + eqEq_rad;
853 
854  /* Create earth rotation matrix from the current Apparent Sidereal Time */
855  v3d_createRotationMatrix(&earthRotM, Zaxis, gast_rad);
856 
857  /* Convert apparent posn to posn in Terrestrial Intermediate Ref System */
858  v3d_multMxV(terInterV, &earthRotM, appV);
859 }
860 
861 
862 /*
863  *------------------------------------------------------------------------------
864  *
865  * Local functions (not called from other modules)
866  *
867  *------------------------------------------------------------------------------
868  */
869 /*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
sky1.h
sky1.h - astronomical coordinate conversion routines, IAU 1980
v3d_multMxM
V3D_Matrix * v3d_multMxM(V3D_Matrix *destM, const V3D_Matrix *srcM1, const V3D_Matrix *srcM2)
Multiply 3x3 matrix by 3x3 matrix to give a new 3x3 matrix, as per [destM] = [srcM1]*[srcM2].
Definition: vectors3d.c:517
general.h
general.h - definitions of general use to (standard) C programs
Sky1_Nut1980
Nutation angles and obliquity.
Definition: sky1.h:75
astron.h
astron.h - assorted definitions useful for astronomy
ARRAY_SIZE
#define ARRAY_SIZE(x__)
Because C passes arrays to functions by passing only a pointer to the zero'th element of the array,...
Definition: general.h:194
sky1_epsilon1980
void sky1_epsilon1980(double t_cy, Sky1_Nut1980 *nut)
Calculate the obliquity of the ecliptic and the equation of the equinoxes.
Definition: sky1.c:631
LOCAL
#define LOCAL
C has some very silly default behaviours.
Definition: general.h:133
PI
#define PI
,
Definition: general.h:158
GLOBAL
#define GLOBAL
See above.
Definition: general.h:134
REQUIRE_NOT_NULL
#define REQUIRE_NOT_NULL(pointer_)
,
Definition: general.h:308
Sky1_Prec1976::zed_rad
double zed_rad
Precession angle (z) (radian)
Definition: sky1.h:70
secToRad
static double secToRad(double angle_s)
Returns angle_s converted from seconds to radians.
Definition: astron.h:50
sky1_createNPmatrix
void sky1_createNPmatrix(double t0_cy, double t1_cy, int precision, V3D_Matrix *npM)
Call the various precession and nutation routines in this module to create a combined precession and ...
Definition: sky1.c:721
Sky1_Nut1980::dPsi_rad
double dPsi_rad
Nutation in longitude (Δψ) (radian)
Definition: sky1.h:76
sky1_frameBiasFK5
void sky1_frameBiasFK5(V3D_Matrix *biasM)
Create the frame bias matrix from the IAU 2000 precession-nutation model.
Definition: sky1.c:339
Sky1_Nut1980::dEps_rad
double dEps_rad
Nutation in obliquity (Δε) (radian)
Definition: sky1.h:77
TWOPI
#define TWOPI
,
Definition: general.h:160
sky1_createNut1980Matrix
void sky1_createNut1980Matrix(const Sky1_Nut1980 *nut, V3D_Matrix *nutM)
This routine calculates the Nutation matrix, using the nutation in longitude, the nutation in obliqui...
Definition: sky1.c:676
sky1_appToTirs
void sky1_appToTirs(const V3D_Vector *appV, double j2kUT1_d, double eqEq_rad, V3D_Vector *terInterV)
Convert a position in geocentric apparent coordinates to geocentric coordinates in the Terrestrial In...
Definition: sky1.c:810
V3D_Matrix
3x3 matrix.
Definition: vectors3d.h:51
arcsecToRad
static double arcsecToRad(double angle_as)
Returns angle_as converted from arcseconds to radians.
Definition: astron.h:43
instead-of-math.h
instead-of-math.h - header to be included instead of math.h
V3D_Vector
3x1 vector.
Definition: vectors3d.h:57
sky1_gmSiderealTimeIAU1982
double sky1_gmSiderealTimeIAU1982(double du)
Calculate the Greenwich mean sidereal time.
Definition: sky1.c:766
sky1_precessionIAU1976
void sky1_precessionIAU1976(double t0, double t1, Sky1_Prec1976 *terms)
This procedure calculates the equatorial precession parameters ζ, z, and θ which represent the rotati...
Definition: sky1.c:375
DEFINE_THIS_FILE
#define DEFINE_THIS_FILE
Leave DEFINE_THIS_FILE undefined.
Definition: general.h:252
Sky1_Nut1980::eps0_rad
double eps0_rad
Mean obliquity of ecliptic at date (ε0)(radian)
Definition: sky1.h:78
Sky1_Prec1976::theta_rad
double theta_rad
Precession angle (θ) (radian)
Definition: sky1.h:71
v3d_createRotationMatrix
V3D_Matrix * v3d_createRotationMatrix(V3D_Matrix *destM, V3D_AxisNames axis, double angle_rad)
Creates a matrix to rotate a coordinate system about an axis.
Definition: vectors3d.c:377
JUL_CENT
#define JUL_CENT
Length of Julian Century in days.
Definition: astron.h:27
sky1_createPrec1976Matrix
void sky1_createPrec1976Matrix(const Sky1_Prec1976 *terms, V3D_Matrix *precM)
This routine calculates the precession matrix, based on angles ζ, z, and θ.
Definition: sky1.c:437
Sky1_Nut1980::eqEq_rad
double eqEq_rad
Equation of the Equinoxes (radian)
Definition: sky1.h:79
sky1_nutationIAU1980
void sky1_nutationIAU1980(double t_cy, int precision, Sky1_Nut1980 *nut)
Calculates the nutation in longitude and obliquity, according to the IAU 1980 Nutation Theory.
Definition: sky1.c:481
Sky1_Prec1976
Precession angles.
Definition: sky1.h:68
v3d_multMxV
V3D_Vector * v3d_multMxV(V3D_Vector *destV, const V3D_Matrix *srcM, const V3D_Vector *srcV)
Multiply 3x3 matrix by 3x1 vector to give a new 3x1 vector, as per equation [destV] = [srcM] * [srcV]...
Definition: vectors3d.c:439
Sky1_Prec1976::zeta_rad
double zeta_rad
Precession angle (ζ) (radian)
Definition: sky1.h:69
normalize
static double normalize(double x, double range)
Normalizes a cyclic double precision floating point variable x to the interval [0,...
Definition: instead-of-math.h:144