← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/gkrellmoon/oneiric

« back to all changes in this revision

Viewing changes to MoonRise.c

  • Committer: Bazaar Package Importer
  • Author(s): Martin Zobel-Helas
  • Date: 2004-11-19 22:55:21 UTC
  • Revision ID: james.westby@ubuntu.com-20041119225521-pmuy6kl327g5wevc
Tags: upstream-0.6
Import upstream version 0.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
MoonRise.c
 
1
/* (C) Mike Henderson <mghenderson@lanl.gov>. 
 
2
 *
 
3
 *  I've converted to glib types, removed unused variables, and piped the whole
 
4
 *  thing through indent.
 
5
 *
 
6
 *  josh buhl <jbuhl@users.sourceforge.net> 
 
7
 */
 
8
 
 
9
#include "Moon.h"
 
10
#include "MoonRise.h"
 
11
 
 
12
void
 
13
UTTohhmm(gdouble UT, gint * h, gint * m)
 
14
{
 
15
        if (UT < 0.0) {
 
16
                *h = -1.0;
 
17
                *m = -1.0;
 
18
        } else {
 
19
                *h = (gint) UT;
 
20
                *m = (gint) ((UT - (gdouble) (*h)) * 60.0 + 0.5);
 
21
                if (*m == 60) {
 
22
                        /* 
 
23
                         *  If it was 23:60 this should become 24:00
 
24
                         *  I prefer this designation to 00:00. So dont reset h to 0 when it goes above 24.
 
25
                         */
 
26
                        *h += 1;
 
27
                        *m = 0;
 
28
                }
 
29
        }
 
30
 
 
31
}
 
32
 
 
33
static gint
 
34
Interp(gdouble ym, gdouble y0, gdouble yp, gdouble * xe, gdouble * ye,
 
35
       gdouble * z1, gdouble * z2, gint * nz)
 
36
{
 
37
 
 
38
        gdouble a, b, c, d, dx;
 
39
 
 
40
        *nz = 0;
 
41
        a = 0.5 * (ym + yp) - y0;
 
42
        b = 0.5 * (yp - ym);
 
43
        c = y0;
 
44
        *xe = -b / (2.0 * a);
 
45
        *ye = (a * (*xe) + b) * (*xe) + c;
 
46
        d = b * b - 4.0 * a * c;
 
47
 
 
48
        if (d >= 0) {
 
49
                dx = 0.5 * sqrt(d) / fabs(a);
 
50
                *z1 = *xe - dx;
 
51
                *z2 = *xe + dx;
 
52
                if (fabs(*z1) <= 1.0)
 
53
                        *nz += 1;
 
54
                if (fabs(*z2) <= 1.0)
 
55
                        *nz += 1;
 
56
                if (*z1 < -1.0)
 
57
                        *z1 = *z2;
 
58
        }
 
59
 
 
60
        return (0);
 
61
}
 
62
 
 
63
 
 
64
static gdouble SinH(CTrans * c, gdouble UT)
 
65
{
 
66
 
 
67
        gdouble TU, TU2, TU3, LambdaMoon, BetaMoon, R, AGE, frac(), jd();
 
68
        gdouble RA_Moon, DEC_Moon, gmst, lmst, Tau, epsilon, Moon();
 
69
        gdouble angle2pi();
 
70
 
 
71
        TU = (jd(c->year, c->month, c->day, UT) - 2451545.0) / 36525.0;
 
72
 
 
73
        /* this is more accurate, but expensive. */
 
74
        TU2 = TU * TU;
 
75
        TU3 = TU2 * TU;
 
76
        Moon(TU, &LambdaMoon, &BetaMoon, &R, &AGE);
 
77
        LambdaMoon *= RadPerDeg;
 
78
        BetaMoon *= RadPerDeg;
 
79
        epsilon =
 
80
            (23.43929167 - 0.013004166 * TU - 1.6666667e-7 * TU2 -
 
81
             5.0277777778e-7 * TU3) * RadPerDeg;
 
82
        RA_Moon =
 
83
            angle2pi(atan2
 
84
                     (sin(LambdaMoon) * cos(epsilon) -
 
85
                      tan(BetaMoon) * sin(epsilon), cos(LambdaMoon)));
 
86
        DEC_Moon =
 
87
            asin(sin(BetaMoon) * cos(epsilon) +
 
88
                 cos(BetaMoon) * sin(epsilon) * sin(LambdaMoon));
 
89
 
 
90
 
 
91
        /* This is less accurate, but computationally less expensive */
 
92
/*      MiniMoon(TU, &RA_Moon, &DEC_Moon); */
 
93
/*      RA_Moon *= 15.0*RadPerDeg; */
 
94
/*      DEC_Moon *= RadPerDeg; */
 
95
 
 
96
        /*
 
97
         *  Compute Greenwich Mean Sidereal Time (gmst)
 
98
         */
 
99
        UT = 24.0 * frac(UT / 24.0);
 
100
        /* this is for the ephemeris meridian??? 
 
101
           gmst = 6.697374558 + 1.0027379093*UT + (8640184.812866+(0.093104-6.2e-6*TU)*TU)*TU/3600.0;
 
102
         */
 
103
        gmst =
 
104
            UT + 6.697374558 + (8640184.812866 +
 
105
                                (0.093104 -
 
106
                                 6.2e-6 * TU) * TU) * TU / 3600.0;
 
107
        lmst = 24.0 * frac((gmst - c->Glon / 15.0) / 24.0);
 
108
 
 
109
        Tau = 15.0 * lmst * RadPerDeg - RA_Moon;
 
110
 
 
111
        return (c->SinGlat * sin(DEC_Moon) +
 
112
                c->CosGlat * cos(DEC_Moon) * cos(Tau));
 
113
}
 
114
 
 
115
 
 
116
void
 
117
MoonRise(CTrans * c, gdouble * UTRise, gdouble * UTSet)
 
118
{
 
119
 
 
120
        gdouble UT, ym, y0, yp, SinH0;
 
121
        gdouble xe, ye, z1, z2;
 
122
        gint Rise, Set, nz, TimeZone;
 
123
 
 
124
        SinH0 = sin(8.0 / 60.0 * RadPerDeg);
 
125
 
 
126
        /* report rise and set times in LST */
 
127
        TimeZone = c->UT - c->LST;
 
128
 
 
129
        UT = 1.0 + TimeZone;
 
130
        *UTRise = -999.0;
 
131
        *UTSet = -999.0;
 
132
        Rise = Set = 0;
 
133
 
 
134
        ym = SinH(c, UT - 1.0) - SinH0;
 
135
 
 
136
        while ((UT <= 24.0 + TimeZone)) {
 
137
 
 
138
                y0 = SinH(c, UT) - SinH0;
 
139
                yp = SinH(c, UT + 1.0) - SinH0;
 
140
 
 
141
                Interp(ym, y0, yp, &xe, &ye, &z1, &z2, &nz);
 
142
 
 
143
                switch (nz) {
 
144
 
 
145
                case 0:
 
146
                        break;
 
147
                case 1:
 
148
                        if (ym < 0.0) {
 
149
                                *UTRise = UT + z1;
 
150
                                Rise = 1;
 
151
                        } else {
 
152
                                *UTSet = UT + z1;
 
153
                                Set = 1;
 
154
                        }
 
155
                        break;
 
156
                case 2:
 
157
                        if (ye < 0.0) {
 
158
                                *UTRise = UT + z2;
 
159
                                *UTSet = UT + z1;
 
160
                        } else {
 
161
                                *UTRise = UT + z1;
 
162
                                *UTSet = UT + z2;
 
163
                        }
 
164
                        Rise = 1;
 
165
                        Set = 1;
 
166
                        break;
 
167
                }
 
168
                ym = yp;
 
169
                UT += 2.0;
 
170
 
 
171
        }
 
172
 
 
173
        if (Rise) {
 
174
                *UTRise -= TimeZone;
 
175
                *UTRise = hour24(*UTRise);
 
176
        } else {
 
177
                *UTRise = -999.0;
 
178
        }
 
179
 
 
180
        if (Set) {
 
181
                *UTSet -= TimeZone;
 
182
                *UTSet = hour24(*UTSet);
 
183
        } else {
 
184
                *UTSet = -999.0;
 
185
        }
 
186
 
 
187
}