1 /*
2 * This library is free software; you can redistribute it and/or
3 * modify it under the terms of the GNU Lesser General Public
4 * License as published by the Free Software Foundation; either
5 * version 2 of the License, or (at your option) any later version.
6 *
7 * This library is distributed in the hope that it will be useful,
8 * but WITHOUT ANY WARRANTY; without even the implied warranty of
9 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
10 * Lesser General Public License for more details.
11 *
12 * You should have received a copy of the GNU General Public License
13 * along with this program; if not, write to the Free Software
14 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
15 *
16 * Copyright (C) 2000 - 2005 Liam Girdwood
17 */
18
19 #include <libnova/dynamical_time.h>
20 #include <stdio.h>
21 #define TERMS 192
22
23 struct year_TD
24 {
25 int year;
26 double TD;
27 };
28
29 /* Stephenson and Houlden for years prior to 948 A.D.*/
30 static double get_dynamical_diff_sh1 (double JD);
31
32 /* Stephenson and Houlden for years between 948 A.D. and 1600 A.D.*/
33 static double get_dynamical_diff_sh2 (double JD);
34
35 /* Table 9.a pg 72 for years 1620..1992.*/
36 static double get_dynamical_diff_table (double JD);
37
38 /* get the dynamical time diff in the near past / future 1992 .. 2010 */
39 static double get_dynamical_diff_near (double JD);
40
41 /* uses equation 9.1 pg 73 to calc JDE for othe JD values */
42 static double get_dynamical_diff_other (double JD);
43
44
45 /* dynamical time in seconds for every second year from 1620 to 1992 */
46 const static double delta_t[TERMS] =
47 { 124.0, 115.0, 106.0, 98.0, 91.0,
48 85.0, 79.0, 74.0, 70.0, 65.0,
49 62.0, 58.0, 55.0, 53.0, 50.0,
50 48.0, 46.0, 44.0, 42.0, 40.0,
51 37.0, 35.0, 33.0, 31.0, 28.0,
52 26.0, 24.0, 22.0, 20.0, 18.0,
53 16.0, 14.0, 13.0, 12.0, 11.0,
54 10.0, 9.0, 9.0, 9.0, 9.0,
55 9.0, 9.0, 9.0, 9.0, 10.0,
56 10.0, 10.0, 10.0, 10.0, 11.0,
57 11.0, 11.0, 11.0, 11.0, 11.0,
58 11.0, 12.0, 12.0, 12.0, 12.0,
59 12.0, 12.0, 13.0, 13.0, 13.0,
60 13.0, 14.0, 14.0, 14.0, 15.0,
61 15.0, 15.0, 15.0, 16.0, 16.0,
62 16.0, 16.0, 16.0, 17.0, 17.0,
63 17.0, 17.0, 17.0, 17.0, 17.0,
64 17.0, 16.0, 16.0, 15.0, 14.0,
65 13.7, 13.1, 12.7, 12.5, 12.5,
66 12.5, 12.5, 12.5, 12.5, 12.3,
67 12.0, 11.4, 10.6, 9.6, 8.6,
68 7.5, 6.6, 6.0, 5.7, 5.6,
69 5.7, 5.9, 6.2, 6.5, 6.8,
70 7.1, 7.3, 7.5, 7.7, 7.8,
71 7.9, 7.5, 6.4, 5.4, 2.9,
72 1.6, -1.0, -2.7, -3.6, -4.7,
73 -5.4, -5.2, -5.5, -5.6, -5.8,
74 -5.9, -6.2, -6.4, -6.1, -4.7,
75 -2.7, 0.0, 2.6, 5.4, 7.7,
76 10.5, 13.4, 16.0, 18.2, 20.2,
77 21.2, 22.4, 23.5, 23.9, 24.3,
78 24.0, 23.9, 23.9, 23.7, 24.0,
79 24.3, 25.3, 26.2, 27.3, 28.2,
80 29.1, 30.0, 30.7, 31.4, 32.2,
81 33.1, 34.0, 35.0, 36.5, 38.3,
82 40.2, 42.2, 44.5, 46.5, 48.5,
83 50.5, 52.2, 53.8, 54.9, 55.8,
84 56.9, 58.3
85 };
86
87
88 /* Stephenson and Houlden for years prior to 948 A.D.*/
get_dynamical_diff_sh1(double JD)89 static double get_dynamical_diff_sh1 (double JD)
90 {
91 double TD,E;
92
93 /* number of centuries from 948 */
94 E = (JD - 2067314.5) / 36525.0;
95
96 TD = 1830.0 - 405.0 * E + 46.5 * E * E;
97 return (TD);
98 }
99
100 /* Stephenson and Houlden for years between 948 A.D. and 1600 A.D.*/
get_dynamical_diff_sh2(double JD)101 static double get_dynamical_diff_sh2 (double JD)
102 {
103 double TD,t;
104
105 /* number of centuries from 1850 */
106 t = (JD - 2396758.5) / 36525.0;
107
108 TD = 22.5 * t * t;
109 return TD;
110 }
111
112 /* Table 9.a pg 72 for years 1600..1992.*/
113 /* uses interpolation formula 3.3 on pg 25 */
get_dynamical_diff_table(double JD)114 static double get_dynamical_diff_table (double JD)
115 {
116 double TD = 0;
117 double a,b,c,n;
118 int i;
119
120 /* get no days since 1620 and divide by 2 years */
121 i = (int)((JD - 2312752.5) / 730.5);
122
123 /* get the base interpolation factor in the table */
124 if (i > (TERMS - 2))
125 i = TERMS - 2;
126
127 /* calc a,b,c,n */
128 a = delta_t[i+1] - delta_t[i];
129 b = delta_t[i+2] - delta_t[i+1];
130 c = a - b;
131 n = ((JD - (2312752.5 + (730.5 * i))) / 730.5);
132
133 TD = delta_t[i+1] + n / 2 * (a + b + n * c);
134
135 return TD;
136 }
137
138 /* get the dynamical time diff in the near past / future 1992 .. 2010 */
139 /* uses interpolation formula 3.3 on pg 25 */
get_dynamical_diff_near(double JD)140 static double get_dynamical_diff_near (double JD)
141 {
142 double TD = 0;
143 /* TD for 1990, 2000, 2010 */
144 double delta_T[3] = {56.86, 63.83, 70.0};
145 double a,b,c,n;
146
147 /* calculate TD by interpolating value */
148 a = delta_T[1] - delta_T[0];
149 b = delta_T[2] - delta_T[1];
150 c = b - a;
151
152 /* get number of days since 2000 and divide by 10 years */
153 n = (JD - 2451544.5) / 3652.5;
154 TD = delta_T[1] + (n / 2) * (a + b + n * c);
155
156 return TD;
157 }
158
159 /* uses equation 9.1 pg 73 to calc JDE for othe JD values */
get_dynamical_diff_other(double JD)160 static double get_dynamical_diff_other (double JD)
161 {
162 double TD;
163 double a;
164
165 a = (JD - 2382148);
166 a *= a;
167
168 TD = -15 + a / 41048480;
169
170 return (TD);
171 }
172
173 /*! \fn double ln_get_dynamical_time_diff (double JD)
174 * \param JD Julian Day
175 * \return TD
176 *
177 * Calculates the dynamical time (TD) difference in seconds (delta T) from
178 * universal time.
179 */
180 /* Equation 9.1 on pg 73.
181 */
ln_get_dynamical_time_diff(double JD)182 double ln_get_dynamical_time_diff (double JD)
183 {
184 double TD;
185
186 /* check when JD is, and use corresponding formula */
187 /* check for date < 948 A.D. */
188 if ( JD < 2067314.5 )
189 /* Stephenson and Houlden */
190 TD = get_dynamical_diff_sh1 (JD);
191 else if ( JD >= 2067314.5 && JD < 2305447.5 )
192 /* check for date 948..1600 A.D. Stephenson and Houlden */
193 TD = get_dynamical_diff_sh2 (JD);
194 else if ( JD >= 2312752.5 && JD < 2448622.5 )
195 /* check for value in table 1620..1992 interpolation of table */
196 TD = get_dynamical_diff_table (JD);
197 else if ( JD >= 2448622.5 && JD <= 2455197.5 )
198 /* check for near future 1992..2010 interpolation */
199 TD = get_dynamical_diff_near (JD);
200 else
201 /* other time period outside */
202 TD = get_dynamical_diff_other (JD);
203
204 return TD;
205 }
206
207
208 /*! \fn double ln_get_jde (double JD)
209 * \param JD Julian Day
210 * \return Julian Ephemeris day
211 *
212 * Calculates the Julian Ephemeris Day (JDE) from the given julian day
213 */
214
ln_get_jde(double JD)215 double ln_get_jde (double JD)
216 {
217 double JDE;
218 double secs_in_day = 24 * 60 * 60;
219
220 JDE = JD + ln_get_dynamical_time_diff (JD) / secs_in_day;
221
222 return JDE;
223 }
224