← Back to branch summary

~ubuntu-branches/debian/squeeze/stellarium/squeeze

« back to all changes in this revision

Viewing changes to src/planetsephems/sideral_time.c

  • Committer: Bazaar Package Importer
  • Author(s): Cédric Delfosse
  • Date: 2008-05-19 21:28:23 UTC
  • mfrom: (3.1.5 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080519212823-m5nfiuntxstxzxj7
Tags: 0.9.1-4
http://bugs.debian.org/479906
Add libxcursor-dev, libxfixes-dev, libxinerama-dev, libqt4-opengl-dev to
build-deps (Closes: #479906)

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/planetsephems/sideral_time.c
 
1
/*
 
2
This program is free software; you can redistribute it and/or modify
 
3
it under the terms of the GNU Library General Public License as published by
 
4
the Free Software Foundation; either version 2 of the License, or
 
5
(at your option) any later version.
 
6
 
 
7
This program 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
 
10
GNU 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 Liam Girdwood <liam@nova-ioe.org>
 
17
 
 
18
*/
 
19
 
 
20
#include <math.h>
 
21
 
 
22
#ifndef M_PI
 
23
#define M_PI           3.14159265358979323846
 
24
#endif
 
25
 
 
26
/* puts a large angle in the correct range 0 - 360 degrees */
 
27
double range_degrees(double d)
 
28
{
 
29
        d = fmod( d, 360.);
 
30
        if(d<0.) d += 360.;
 
31
        return d;
 
32
}
 
33
 
 
34
/* puts a large angle in the correct range 0 - 2PI radians */
 
35
double range_radians (double r)
 
36
{
 
37
        r = fmod( r, 2.*M_PI );
 
38
        if (r<0.) r += 2.*M_PI;
 
39
        return r;
 
40
}
 
41
 
 
42
#define TERMS 63
 
43
#define LN_NUTATION_EPOCH_THRESHOLD 0.1
 
44
 
 
45
 
 
46
/* Nutation is a period oscillation of the Earths rotational axis around it's mean position.*/
 
47
 
 
48
/*
 
49
 Contains Nutation in longitude, obliquity and ecliptic obliquity.
 
50
 Angles are expressed in degrees.
 
51
*/
 
52
struct ln_nutation
 
53
{
 
54
        double longitude;       /*!< Nutation in longitude */
 
55
        double obliquity;       /*!< Nutation in obliquity */
 
56
        double ecliptic;        /*!< Obliquity of the ecliptic */
 
57
};
 
58
 
 
59
struct nutation_arguments
 
60
{
 
61
        double D;
 
62
        double M;
 
63
        double MM;
 
64
        double F;
 
65
        double O;
 
66
};
 
67
 
 
68
struct nutation_coefficients
 
69
{
 
70
        double longitude1;
 
71
        double longitude2;
 
72
        double obliquity1;
 
73
        double obliquity2;
 
74
};
 
75
 
 
76
/* arguments and coefficients taken from table 21A on page 133 */
 
77
 
 
78
const static struct nutation_arguments arguments[TERMS] = {
 
79
        {0.0,   0.0,    0.0,    0.0,    1.0},
 
80
        {-2.0,  0.0,    0.0,    2.0,    2.0},
 
81
        {0.0,   0.0,    0.0,    2.0,    2.0},
 
82
        {0.0,   0.0,    0.0,    0.0,    2.0},
 
83
        {0.0,   1.0,    0.0,    0.0,    0.0},
 
84
        {0.0,   0.0,    1.0,    0.0,    0.0},
 
85
        {-2.0,  1.0,    0.0,    2.0,    2.0},
 
86
        {0.0,   0.0,    0.0,    2.0,    1.0},
 
87
        {0.0,   0.0,    1.0,    2.0,    2.0},
 
88
        {-2.0,  -1.0,   0.0,    2.0,    2.0},
 
89
        {-2.0,  0.0,    1.0,    0.0,    0.0},
 
90
        {-2.0,  0.0,    0.0,    2.0,    1.0},
 
91
        {0.0,   0.0,    -1.0,   2.0,    2.0},
 
92
        {2.0,   0.0,    0.0,    0.0,    0.0},
 
93
        {0.0,   0.0,    1.0,    0.0,    1.0},
 
94
        {2.0,   0.0,    -1.0,   2.0,    2.0},
 
95
        {0.0,   0.0,    -1.0,   0.0,    1.0},
 
96
        {0.0,   0.0,    1.0,    2.0,    1.0},
 
97
        {-2.0,  0.0,    2.0,    0.0,    0.0},
 
98
        {0.0,   0.0,    -2.0,   2.0,    1.0},
 
99
        {2.0,   0.0,    0.0,    2.0,    2.0},
 
100
        {0.0,   0.0,    2.0,    2.0,    2.0},
 
101
        {0.0,   0.0,    2.0,    0.0,    0.0},
 
102
        {-2.0,  0.0,    1.0,    2.0,    2.0},
 
103
        {0.0,   0.0,    0.0,    2.0,    0.0},
 
104
        {-2.0,  0.0,    0.0,    2.0,    0.0},
 
105
        {0.0,   0.0,    -1.0,   2.0,    1.0},
 
106
        {0.0,   2.0,    0.0,    0.0,    0.0},
 
107
        {2.0,   0.0,    -1.0,   0.0,    1.0},
 
108
        {-2.0,  2.0,    0.0,    2.0,    2.0},
 
109
        {0.0,   1.0,    0.0,    0.0,    1.0},
 
110
        {-2.0,  0.0,    1.0,    0.0,    1.0},
 
111
        {0.0,   -1.0,   0.0,    0.0,    1.0},
 
112
        {0.0,   0.0,    2.0,    -2.0,   0.0},
 
113
        {2.0,   0.0,    -1.0,   2.0,    1.0},
 
114
        {2.0,   0.0,    1.0,    2.0,    2.0},
 
115
        {0.0,   1.0,    0.0,    2.0,    2.0},
 
116
        {-2.0,  1.0,    1.0,    0.0,    0.0},
 
117
        {0.0,   -1.0,   0.0,    2.0,    2.0},
 
118
        {2.0,   0.0,    0.0,    2.0,    1.0},
 
119
        {2.0,   0.0,    1.0,    0.0,    0.0},
 
120
        {-2.0,  0.0,    2.0,    2.0,    2.0},
 
121
        {-2.0,  0.0,    1.0,    2.0,    1.0},
 
122
        {2.0,   0.0,    -2.0,   0.0,    1.0},
 
123
        {2.0,   0.0,    0.0,    0.0,    1.0},
 
124
        {0.0,   -1.0,   1.0,    0.0,    0.0},
 
125
        {-2.0,  -1.0,   0.0,    2.0,    1.0},
 
126
        {-2.0,  0.0,    0.0,    0.0,    1.0},
 
127
        {0.0,   0.0,    2.0,    2.0,    1.0},
 
128
        {-2.0,  0.0,    2.0,    0.0,    1.0},
 
129
        {-2.0,  1.0,    0.0,    2.0,    1.0},
 
130
        {0.0,   0.0,    1.0,    -2.0,   0.0},
 
131
        {-1.0,  0.0,    1.0,    0.0,    0.0},
 
132
        {-2.0,  1.0,    0.0,    0.0,    0.0},
 
133
        {1.0,   0.0,    0.0,    0.0,    0.0},
 
134
        {0.0,   0.0,    1.0,    2.0,    0.0},
 
135
        {0.0,   0.0,    -2.0,   2.0,    2.0},
 
136
        {-1.0,  -1.0,   1.0,    0.0,    0.0},
 
137
        {0.0,   1.0,    1.0,    0.0,    0.0},
 
138
        {0.0,   -1.0,   1.0,    2.0,    2.0},
 
139
        {2.0,   -1.0,   -1.0,   2.0,    2.0},
 
140
        {0.0,   0.0,    3.0,    2.0,    2.0},
 
141
        {2.0,   -1.0,   0.0,    2.0,    2.0}};
 
142
 
 
143
const static struct nutation_coefficients coefficients[TERMS] = {
 
144
        {-171996.0,     -174.2, 92025.0,8.9},
 
145
        {-13187.0,      -1.6,   5736.0, -3.1},
 
146
        {-2274.0,        0.2,   977.0,  -0.5},
 
147
        {2062.0,        0.2,    -895.0,    0.5},
 
148
        {1426.0,    -3.4,    54.0,    -0.1},
 
149
        {712.0,    0.1,    -7.0,    0.0},
 
150
        {-517.0,    1.2,    224.0,    -0.6},
 
151
        {-386.0,    -0.4,    200.0,    0.0},
 
152
        {-301.0,    0.0,    129.0,    -0.1},
 
153
        {217.0,    -0.5,    -95.0,    0.3},
 
154
        {-158.0,    0.0,    0.0,    0.0},
 
155
        {129.0, 0.1,    -70.0,  0.0},
 
156
        {123.0, 0.0,    -53.0,  0.0},
 
157
        {63.0,  0.0,    0.0,    0.0},
 
158
        {63.0,  1.0,    -33.0,  0.0},
 
159
        {-59.0, 0.0,    26.0,   0.0},
 
160
        {-58.0, -0.1,   32.0,   0.0},
 
161
        {-51.0, 0.0,    27.0,   0.0},
 
162
        {48.0,  0.0,    0.0,    0.0},
 
163
        {46.0,  0.0,    -24.0,  0.0},
 
164
        {-38.0, 0.0,    16.0,   0.0},
 
165
        {-31.0, 0.0,    13.0,   0.0},
 
166
        {29.0,  0.0,    0.0,    0.0},
 
167
        {29.0,  0.0,    -12.0,  0.0},
 
168
        {26.0,  0.0,    0.0,    0.0},
 
169
        {-22.0, 0.0,    0.0,    0.0},
 
170
        {21.0,  0.0,    -10.0,  0.0},
 
171
        {17.0,  -0.1,   0.0,    0.0},
 
172
        {16.0,  0.0,    -8.0,   0.0},
 
173
        {-16.0, 0.1,    7.0,    0.0},
 
174
        {-15.0, 0.0,    9.0,    0.0},
 
175
        {-13.0, 0.0,    7.0,    0.0},
 
176
        {-12.0, 0.0,    6.0,    0.0},
 
177
        {11.0,  0.0,    0.0,    0.0},
 
178
        {-10.0, 0.0,    5.0,    0.0},
 
179
        {-8.0,  0.0,    3.0,    0.0},
 
180
        {7.0,   0.0,    -3.0,   0.0},
 
181
        {-7.0,  0.0,    0.0,    0.0},
 
182
        {-7.0,  0.0,    3.0,    0.0},
 
183
        {-7.0,  0.0,    3.0,    0.0},
 
184
        {6.0,   0.0,    0.0,    0.0},
 
185
        {6.0,   0.0,    -3.0,   0.0},
 
186
        {6.0,   0.0,    -3.0,   0.0},
 
187
        {-6.0,  0.0,    3.0,    0.0},
 
188
        {-6.0,  0.0,    3.0,    0.0},
 
189
        {5.0,   0.0,    0.0,    0.0},
 
190
        {-5.0,  0.0,    3.0,    0.0},
 
191
        {-5.0,  0.0,    3.0,    0.0},
 
192
        {-5.0,  0.0,    3.0,    0.0},
 
193
        {4.0,   0.0,    0.0,    0.0},
 
194
        {4.0,   0.0,    0.0,    0.0},
 
195
        {4.0,   0.0,    0.0,    0.0},
 
196
        {-4.0,  0.0,    0.0,    0.0},
 
197
        {-4.0,  0.0,    0.0,    0.0},
 
198
        {-4.0,  0.0,    0.0,    0.0},
 
199
        {3.0,   0.0,    0.0,    0.0},
 
200
        {-3.0,  0.0,    0.0,    0.0},
 
201
        {-3.0,  0.0,    0.0,    0.0},
 
202
        {-3.0,  0.0,    0.0,    0.0},
 
203
        {-3.0,  0.0,    0.0,    0.0},
 
204
        {-3.0,  0.0,    0.0,    0.0},
 
205
        {-3.0,  0.0,    0.0,    0.0},
 
206
        {-3.0,  0.0,    0.0,    0.0}};
 
207
 
 
208
/* cache values */
 
209
static double c_JD = 0.0, c_longitude = 0.0, c_obliquity = 0.0, c_ecliptic = 0.0; 
 
210
 
 
211
 
 
212
/* Calculate nutation of longitude and obliquity in degrees from Julian Ephemeris Day
 
213
* params : JD Julian Day, nutation Pointer to store nutation.
 
214
* Chapter 21 pg 131-134 Using Table 21A */
 
215
void get_nutation (double JD, struct ln_nutation * nutation)
 
216
{
 
217
 
 
218
        double D,M,MM,F,O,T,T2,T3;
 
219
        double coeff_sine, coeff_cos;
 
220
        int i;
 
221
 
 
222
        /* should we bother recalculating nutation */
 
223
        if (fabs(JD - c_JD) > LN_NUTATION_EPOCH_THRESHOLD)
 
224
        {
 
225
                /* set the new epoch */
 
226
                c_JD = JD;
 
227
 
 
228
                /* set ecliptic */
 
229
                c_ecliptic = 23.0 + 26.0 / 60.0 + 27.407 / 3600.0;
 
230
 
 
231
                /* calc T */
 
232
                T = (JD - 2451545.0)/36525;
 
233
                T2 = T * T;
 
234
                T3 = T2 * T;
 
235
 
 
236
                /* calculate D,M,M',F and Omega */
 
237
                D = 297.85036 + 445267.111480 * T - 0.0019142 * T2 + T3 / 189474.0;
 
238
                M = 357.52772 + 35999.050340 * T - 0.0001603 * T2 - T3 / 300000.0;
 
239
                MM = 134.96298 + 477198.867398 * T + 0.0086972 * T2 + T3 / 56250.0;
 
240
                F = 93.2719100 + 483202.017538 * T - 0.0036825 * T2 + T3 / 327270.0;
 
241
                O = 125.04452 - 1934.136261 * T + 0.0020708 * T2 + T3 / 450000.0;
 
242
 
 
243
                /* convert to radians */
 
244
                D *= M_PI/180.;
 
245
                M *= M_PI/180.;
 
246
                MM *= M_PI/180.;
 
247
                F *= M_PI/180.;
 
248
                O *= M_PI/180.;
 
249
 
 
250
                /* calc sum of terms in table 21A */
 
251
                for (i=0; i< TERMS; i++)
 
252
                {
 
253
                        /* calc coefficients of sine and cosine */
 
254
                        coeff_sine = (coefficients[i].longitude1 + (coefficients[i].longitude2 * T));
 
255
                        coeff_cos = (coefficients[i].obliquity1 + (coefficients[i].obliquity2 * T));
 
256
        
 
257
                        /* sum the arguments */
 
258
                        if (arguments[i].D != 0)
 
259
                        {
 
260
                                c_longitude += coeff_sine * (sin (arguments[i].D * D));
 
261
                                c_obliquity += coeff_cos * (cos (arguments[i].D * D));
 
262
                        }
 
263
                        if (arguments[i].M != 0)
 
264
                        {
 
265
                                c_longitude += coeff_sine * (sin (arguments[i].M * M));
 
266
                                c_obliquity += coeff_cos * (cos (arguments[i].M * M));
 
267
                        }
 
268
                        if (arguments[i].MM != 0)
 
269
                        {
 
270
                                c_longitude += coeff_sine * (sin (arguments[i].MM * MM));
 
271
                                c_obliquity += coeff_cos * (cos (arguments[i].MM * MM));
 
272
                        }
 
273
                        if (arguments[i].F != 0)
 
274
                        {
 
275
                                c_longitude += coeff_sine * (sin (arguments[i].F * F));
 
276
                                c_obliquity += coeff_cos * (cos (arguments[i].F * F));
 
277
                        }
 
278
                        if (arguments[i].O != 0)
 
279
                        {
 
280
                                c_longitude += coeff_sine * (sin (arguments[i].O * O));
 
281
                                c_obliquity += coeff_cos * (cos (arguments[i].O * O));
 
282
                        }
 
283
                }    
 
284
 
 
285
                /* change to degrees */
 
286
                c_longitude /= 36000000.;
 
287
                c_obliquity /= 36000000.;
 
288
 
 
289
                c_ecliptic += c_obliquity;
 
290
        }
 
291
 
 
292
        /* return results */
 
293
        nutation->longitude = c_longitude;
 
294
        nutation->obliquity = c_obliquity;
 
295
        nutation->ecliptic = c_ecliptic;
 
296
}
 
297
 
 
298
/* Calculate the mean sidereal time at the meridian of Greenwich of a given date.
 
299
 * returns apparent sidereal time (degree).
 
300
 * Formula 11.1, 11.4 pg 83 */
 
301
double get_mean_sidereal_time (double JD)
 
302
{
 
303
    double sidereal;
 
304
    double T;
 
305
    
 
306
    T = (JD - 2451545.0) / 36525.0;
 
307
        
 
308
    /* calc mean angle */
 
309
    sidereal = 280.46061837 + (360.98564736629 * (JD - 2451545.0)) + (0.000387933 * T * T) - (T * T * T / 38710000.0);
 
310
    
 
311
    /* add a convenient multiple of 360 degrees */
 
312
    sidereal = range_degrees (sidereal);
 
313
 
 
314
    return sidereal;
 
315
} 
 
316
 
 
317
 
 
318
/* Calculate the apparent sidereal time at the meridian of Greenwich of a given date.
 
319
 * returns apparent sidereal time (degree).
 
320
 * Formula 11.1, 11.4 pg 83 */
 
321
double get_apparent_sidereal_time (double JD)
 
322
{
 
323
   double correction, sidereal;
 
324
   struct ln_nutation nutation;  
 
325
   
 
326
   /* get the mean sidereal time */
 
327
   sidereal = get_mean_sidereal_time (JD);
 
328
        
 
329
   /* add corrections for nutation in longitude and for the true obliquity of 
 
330
   the ecliptic */   
 
331
   get_nutation (JD, &nutation); 
 
332
    
 
333
   correction = (nutation.longitude * cos (nutation.obliquity*M_PI/180.));
 
334
 
 
335
   sidereal += correction;
 
336
   
 
337
   return (sidereal);
 
338
}