forked from janczer/goMoonPhase
-
Notifications
You must be signed in to change notification settings - Fork 0
/
package.go
337 lines (283 loc) · 9.65 KB
/
package.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
package moonphase
import (
"math"
"time"
)
type Moon struct {
phase float64
illum float64
age float64
dist float64
angdia float64
sundist float64
sunangdia float64
pdata float64
quarters [8]float64
timespace float64
}
var synmonth float64 = 29.53058868 // Synodic month (new Moon to new Moon)
func New(t time.Time) (moonP *Moon) {
moonP = new(Moon)
// Astronomical constants
var epoch float64 = 2444238.5 // 1989 January 0.0
//Constants defining the Sun's apparent orbit
var elonge float64 = 278.833540 // Ecliptic longitude of the Sun at epoch 1980.0
var elongp float64 = 282.596403 // Ecliptic longitude of the Sun at perigee
var eccent float64 = 0.016718 // Eccentricity of Earth's orbit
var sunsmax float64 = 1.495985e8 // Sun's angular size, degrees, at semi-major axis distance
var sunangsiz float64 = 0.533128
// Elements of the Moon's orbit, epoch 1980.0
var mmlong float64 = 64.975464 // Moon's mean longitude at the epoch
var mmlongp float64 = 349.383063 // Mean longitude of the perigee at the epoch
var mecc float64 = 0.054900 // Eccentricity of the Moon's orbit
var mangsiz float64 = 0.5181 // Moon's angular size at distance a from Earth
var msmax float64 = 384401 // Semi-major axis of Moon's orbit in km
moonP.timespace = float64(t.Unix())
moonP.pdata = utcToJulian(float64(t.Unix()))
// Calculation of the Sun's position
var day = moonP.pdata - epoch // Date within epoch
var n float64 = fixangle((360 / 365.2422) * day) // Mean anomaly of the Sun
var m float64 = fixangle(n + elonge - elongp) // Convert from perigee co-orginates to epoch 1980.0
var ec = kepler(m, eccent) // Solve equation of Kepler
ec = math.Sqrt((1+eccent)/(1-eccent)) * math.Tan(ec/2)
ec = 2 * rad2deg(math.Atan(ec)) // True anomaly
var lambdasun float64 = fixangle(ec + elongp) // Sun's geocentric ecliptic longitude
var f float64 = ((1 + eccent*cos(deg2rad(ec))) / (1 - eccent*eccent)) // Orbital distance factor
var sunDist float64 = sunsmax / f // Distance to Sun in km
var sunAng float64 = f * sunangsiz // Sun's angular size in degrees
// Calsulation of the Moon's position
var ml float64 = fixangle(13.1763966*day + mmlong) // Moon's mean longitude
var mm float64 = fixangle(ml - 0.1114041*day - mmlongp) // Moon's mean anomaly
var ev float64 = 1.2739 * sin(deg2rad(2*(ml-lambdasun)-mm)) // Evection
var ae float64 = 0.1858 * sin(deg2rad(m)) // Annual equation
var a3 float64 = 0.37 * sin(deg2rad(m)) // Correction term
var mmP float64 = mm + ev - ae - a3 // Corrected anomaly
var mec float64 = 6.2886 * sin(deg2rad(mmP)) // Correction for the equation of the centre
var a4 float64 = 0.214 * sin(deg2rad(2*mmP)) // Another correction term
var lP float64 = ml + ev + mec - ae + a4 // Corrected longitude
var v float64 = 0.6583 * sin(deg2rad(2*(lP-lambdasun))) // Variation
var lPP float64 = lP + v // True longitude
// Calculation of the phase of the Moon
var moonAge float64 = lPP - lambdasun // Age of the Moon in degrees
var moonPhase float64 = (1 - cos(deg2rad(moonAge))) / 2 // Phase of the Moon
// Distance of moon from the centre of the Earth
var moonDist float64 = (msmax * (1 - mecc*mecc)) / (1 + mecc*cos(deg2rad(mmP+mec)))
var moonDFrac float64 = moonDist / msmax
var moonAng float64 = mangsiz / moonDFrac // Moon's angular diameter
// store result
moonP.phase = fixangle(moonAge) / 360 // Phase (0 to 1)
moonP.illum = moonPhase // Illuminated fraction (0 to 1)
moonP.age = synmonth * moonP.phase // Age of moon (days)
moonP.dist = moonDist // Distance (kilometres)
moonP.angdia = moonAng // Angular diameter (degreees)
moonP.sundist = sunDist // Distance to Sun (kilometres)
moonP.sunangdia = sunAng // Sun's angular diameter (degrees)
moonP.phaseHunt()
return moonP
}
func sin(a float64) float64 {
return math.Sin(a)
}
func cos(a float64) float64 {
return math.Cos(a)
}
func rad2deg(r float64) float64 {
return (r * 180) / math.Pi
}
func deg2rad(d float64) float64 {
return (d * math.Pi) / 180
}
func fixangle(a float64) float64 {
return (a - 360*math.Floor(a/360))
}
func kepler(m float64, ecc float64) float64 {
epsilon := 0.000001
m = deg2rad(m)
e := m
var delta float64
delta = e - ecc*math.Sin(e) - m
e -= delta / (1 - ecc*math.Cos(e))
for math.Abs(delta) > epsilon {
delta = e - ecc*math.Sin(e) - m
e -= delta / (1 - ecc*math.Cos(e))
}
return e
}
func (m *Moon) phaseHunt() {
var sdate float64 = utcToJulian(m.timespace)
var adate float64 = sdate - 45
var ats float64 = m.timespace - 86400*45
t := time.Unix(int64(ats), 0)
var yy float64 = float64(t.Year())
var mm float64 = float64(t.Month())
var k1 float64 = math.Floor(float64(yy+((mm-1)*(1/12))-1900) * 12.3685)
var nt1 float64 = meanPhase(adate, k1)
adate = nt1
var nt2, k2 float64
for {
adate += synmonth
k2 = k1 + 1
nt2 = meanPhase(adate, k2)
if math.Abs(nt2-sdate) < 0.75 {
nt2 = truePhase(k2, 0.0)
}
if nt1 <= sdate && nt2 > sdate {
break
}
nt1 = nt2
k1 = k2
}
var data [8]float64
data[0] = truePhase(k1, 0.0)
data[1] = truePhase(k1, 0.25)
data[2] = truePhase(k1, 0.5)
data[3] = truePhase(k1, 0.75)
data[4] = truePhase(k2, 0.0)
data[5] = truePhase(k2, 0.25)
data[6] = truePhase(k2, 0.5)
data[7] = truePhase(k2, 0.75)
for i := 0; i < 8; i++ {
m.quarters[i] = (data[i] - 2440587.5) * 86400 // convert to UNIX time
}
}
func utcToJulian(t float64) float64 {
return t/86400 + 2440587.5
}
func julianToUtc(t float64) float64 {
return t*86400 + 2440587.5
}
/**
Calculates time of the mean new Moon for a given
base date. This argument K to this function is the
precomputed synodic month index, given by:
K = (year - 1900) * 12.3685
where year is expressed as a year aand fractional year
*/
func meanPhase(sdate float64, k float64) float64 {
// Time in Julian centuries from 1900 January 0.5
var t float64 = (sdate - 2415020.0) / 36525
var t2 float64 = t * t
var t3 float64 = t2 * t
var nt float64
nt = 2415020.75933 + synmonth*k +
0.0001178*t2 -
0.000000155*t3 +
0.00033*sin(deg2rad(166.56+132.87*t-0.009173*t2))
return nt
}
func truePhase(k float64, phase float64) float64 {
k += phase // Add phase to new moon time
var t float64 = k / 1236.85 // Time in Julian centures from 1900 January 0.5
var t2 float64 = t * t
var t3 float64 = t2 * t
var pt float64
pt = 2415020.75933 + synmonth*k +
0.0001178*t2 -
0.000000155*t3 +
0.00033*sin(deg2rad(166.56+132.87*t-0.009173*t2))
var m, mprime, f float64
m = 359.2242 + 29.10535608*k - 0.0000333*t2 - 0.00000347*t3 // Sun's mean anomaly
mprime = 306.0253 + 385.81691806*k + 0.0107306*t2 + 0.00001236*t3 // Moon's mean anomaly
f = 21.2964 + 390.67050646*k - 0.0016528*t2 - 0.00000239*t3 // Moon's argument of latitude
if phase < 0.01 || math.Abs(phase-0.5) < 0.01 {
// Corrections for New and Full Moon
pt += (0.1734-0.000393*t)*sin(deg2rad(m)) +
0.0021*sin(deg2rad(2*m)) -
0.4068*sin(deg2rad(mprime)) +
0.0161*sin(deg2rad(2*mprime)) -
0.0004*sin(deg2rad(3*mprime)) +
0.0104*sin(deg2rad(2*f)) -
0.0051*sin(deg2rad(m+mprime)) -
0.0074*sin(deg2rad(m-mprime)) +
0.0004*sin(deg2rad(2*f+m)) -
0.0004*sin(deg2rad(2*f-m)) -
0.0006*sin(deg2rad(2*f+mprime)) +
0.0010*sin(deg2rad(2*f-mprime)) +
0.0005*sin(deg2rad(m+2*mprime))
} else if math.Abs(phase-0.25) < 0.01 || math.Abs(phase-0.75) < 0.01 {
pt += (0.1721-0.0004*t)*sin(deg2rad(m)) +
0.0021*sin(deg2rad(2*m)) -
0.6280*sin(deg2rad(mprime)) +
0.0089*sin(deg2rad(2*mprime)) -
0.0004*sin(deg2rad(3*mprime)) +
0.0079*sin(deg2rad(2*f)) -
0.0119*sin(deg2rad(m+mprime)) -
0.0047*sin(deg2rad(m-mprime)) +
0.0003*sin(deg2rad(2*f+m)) -
0.0004*sin(deg2rad(2*f-m)) -
0.0006*sin(deg2rad(2*f+mprime)) +
0.0021*sin(deg2rad(2*f-mprime)) +
0.0003*sin(deg2rad(m+2*mprime)) +
0.0004*sin(deg2rad(m-2*mprime)) -
0.0003*sin(deg2rad(2*m+mprime))
if phase < 0.5 { // First quarter correction
pt += 0.0028 - 0.0004*cos(deg2rad(m)) + 0.0003*cos(deg2rad(mprime))
} else { // Last quarter correction
pt += -0.0028 + 0.0004*cos(deg2rad(m)) - 0.0003*cos(deg2rad(mprime))
}
}
return pt
}
//func (m *Moon) getPhase(n int8) float64 {
// return m.quarters[n]
//}
func (m *Moon) Phase() float64 {
return m.phase
}
func (m *Moon) Illumination() float64 {
return m.illum
}
func (m *Moon) Age() float64 {
return m.age
}
func (m *Moon) Distance() float64 {
return m.dist
}
func (m *Moon) Diameter() float64 {
return m.angdia
}
func (m *Moon) SunDistance() float64 {
return m.sundist
}
func (m *Moon) SunDiameter() float64 {
return m.sunangdia
}
func (m *Moon) NewMoon() float64 {
return m.quarters[0]
}
func (m *Moon) FirstQuarter() float64 {
return m.quarters[1]
}
func (m *Moon) FullMoon() float64 {
return m.quarters[2]
}
func (m *Moon) LastQuarter() float64 {
return m.quarters[3]
}
func (m *Moon) NextNewMoon() float64 {
return m.quarters[4]
}
func (m *Moon) NextFirstQuarter() float64 {
return m.quarters[1]
}
func (m *Moon) NextFullMoon() float64 {
return m.quarters[6]
}
func (m *Moon) NextLastQuarter() float64 {
return m.quarters[7]
}
func (m *Moon) PhaseName() string {
names := map[int]string{
0: "New Moon",
1: "Waxing Crescent",
2: "First Quarter",
3: "Waxing Gibbous",
4: "Full Moon",
5: "Waning Gibbous",
6: "Third Quarter",
7: "Waning Crescent",
8: "New Moon",
}
i := int(math.Floor((m.phase + 0.0625) * 8))
return names[i]
}