forked from neihousaigaai/own-lunar-calendar
-
Notifications
You must be signed in to change notification settings - Fork 0
/
AL_offline.py
195 lines (171 loc) · 6.66 KB
/
AL_offline.py
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
"""
Download from: https://www.informatik.uni-leipzig.de/~duc/amlich/AL.py
(c) 2006 Ho Ngoc Duc.
Astronomical algorithms
from the book "Astronomical Algorithms" by Jean Meeus, 1998
"""
import math
def jdFromDate(dd, mm, yy):
# Compute the (integral) Julian day number of day dd/mm/yyyy,
# i.e., the number of days between 1/1/4713 BC (Julian calendar) and dd/mm/yyyy.
a = int((14 - mm) / 12.)
y = yy + 4800 - a
m = mm + 12*a - 3
jd = dd + int((153*m + 2) / 5.) + 365*y + int(y/4.) - int(y/100.) + int(y/400.) - 32045
if (jd < 2299161):
jd = dd + int((153*m + 2)/5.) + 365*y + int(y/4.) - 32083
return jd
def jdToDate(jd): # Convert a Julian day number to day/month/year. jd is an integer.
if (jd > 2299160): ## After 5/10/1582, Gregorian calendar
a = jd + 32044
b = int((4*a + 3) / 146097.)
c = a - int((b*146097) / 4.)
else:
b = 0
c = jd + 32082
d = int((4*c + 3) / 1461.)
e = c - int((1461*d) / 4.)
m = int((5*e + 2) / 153.)
day = e - int((153*m + 2) / 5.) + 1
month = m + 3 - 12*int(m / 10.)
year = b*100 + d - 4800 + int(m / 10.)
return [day, month, year]
def NewMoon(k):
# Compute the time of the k-th new moon after the new moon of 1/1/1900 13:52 UCT
# (measured as the number of days since 1/1/4713 BC noon UCT, e.g., 2451545.125 is 1/1/2000 15:00 UTC.
# Returns a floating number, e.g., 2415079.9758617813 for k=2 or 2414961.935157746 for k=-2.
## Time in Julian centuries from 1900 January 0.5
T = k / 1236.85
T2 = T * T
T3 = T2 * T
dr = math.pi / 180.
Jd1 = 2415020.75933 + 29.53058868*k + 0.0001178*T2 - 0.000000155*T3
Jd1 = Jd1 + 0.00033*math.sin((166.56 + 132.87*T - 0.009173*T2)*dr)
## Mean new moon
M = 359.2242 + 29.10535608*k - 0.0000333*T2 - 0.00000347*T3
## Sun's mean anomaly
Mpr = 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
C1 = (0.1734 - 0.000393*T)*math.sin(M*dr) + 0.0021*math.sin(2*dr*M)
C1 = C1 - 0.4068*math.sin(Mpr*dr) + 0.0161*math.sin(dr*2*Mpr)
C1 = C1 - 0.0004*math.sin(dr*3*Mpr)
C1 = C1 + 0.0104*math.sin(dr*2*F) - 0.0051*math.sin(dr*(M + Mpr))
C1 = C1 - 0.0074*math.sin(dr*(M - Mpr)) + 0.0004*math.sin(dr*(2*F + M))
C1 = C1 - 0.0004*math.sin(dr*(2*F - M)) - 0.0006*math.sin(dr*(2*F + Mpr))
C1 = C1 + 0.0010*math.sin(dr*(2*F - Mpr)) + 0.0005*math.sin(dr*(2*Mpr + M))
if (T < -11):
deltat= 0.001 + 0.000839*T + 0.0002261*T2 - 0.00000845*T3 - 0.000000081*T*T3
else:
deltat= -0.000278 + 0.000265*T + 0.000262*T2
JdNew = Jd1 + C1 - deltat
return JdNew
def SunLongitude(jdn):
# Compute the longitude of the sun at any time.
# Parameter: floating number jdn, the number of days since 1/1/4713 BC noon.
T = (jdn - 2451545.0) / 36525.
## Time in Julian centuries
## from 2000-01-01 12:00:00 GMT
T2 = T * T
dr = math.pi / 180. ## degree to radian
M = 357.52910 + 35999.05030*T - 0.0001559*T2 - 0.00000048*T*T2
## mean anomaly, degree
L0 = 280.46645 + 36000.76983*T + 0.0003032*T2
## mean longitude, degree
DL = (1.914600 - 0.004817*T - 0.000014*T2) * math.sin(dr*M)
DL += (0.019993 - 0.000101*T) *math.sin(dr*2*M) + 0.000290*math.sin(dr*3*M)
L = L0 + DL ## true longitude, degree
L = L * dr
L = L - math.pi*2*(int(L / (math.pi*2)))
#### Normalize to (0, 2*math.pi)
return L
def getSunLongitude(dayNumber, timeZone):
# Compute sun position at midnight of the day with the given Julian day number.
# The time zone if the time difference between local time and UTC: 7.0 for UTC+7:00.
# The function returns a number between 0 and 11.
# From the day after March equinox and the 1st major term after March equinox, 0 is returned.
# After that, return 1, 2, 3 ...
return int(SunLongitude(dayNumber - 0.5 - timeZone/24.) / math.pi*6)
def getNewMoonDay(k, timeZone):
# Compute the day of the k-th new moon in the given time zone.
# The time zone if the time difference between local time and UTC: 7.0 for UTC+7:00.
return int(NewMoon(k) + 0.5 + timeZone / 24.)
def getLunarMonth11(yy, timeZone):
# Find the day that starts the lunar month 11of the given year for the given time zone.
# off = jdFromDate(31, 12, yy) - 2415021.076998695
off = jdFromDate(31, 12, yy) - 2415021.
k = int(off / 29.530588853)
nm = getNewMoonDay(k, timeZone)
sunLong = getSunLongitude(nm, timeZone)
#### sun longitude at local midnight
if (sunLong >= 9):
nm = getNewMoonDay(k - 1, timeZone)
return nm
def getLeapMonthOffset(a11, timeZone):
# Find the index of the leap month after the month starting on the day a11.
k = int((a11 - 2415021.076998695) / 29.530588853 + 0.5)
last = 0
i = 1 ## start with month following lunar month 11
arc = getSunLongitude(getNewMoonDay(k + i, timeZone), timeZone)
while True:
last = arc
i += 1
arc = getSunLongitude(getNewMoonDay(k + i, timeZone), timeZone)
if not (arc != last and i < 14):
break
return i - 1
def S2L(dd, mm, yy, timeZone = 7):
# Convert solar date dd/mm/yyyy to the corresponding lunar date.
dayNumber = jdFromDate(dd, mm, yy)
k = int((dayNumber - 2415021.076998695) / 29.530588853)
monthStart = getNewMoonDay(k + 1, timeZone)
if (monthStart > dayNumber):
monthStart = getNewMoonDay(k, timeZone)
# alert(dayNumber + " -> " + monthStart)
a11 = getLunarMonth11(yy, timeZone)
b11 = a11
if (a11 >= monthStart):
lunarYear = yy
a11 = getLunarMonth11(yy - 1, timeZone)
else:
lunarYear = yy + 1
b11 = getLunarMonth11(yy + 1, timeZone)
lunarDay = dayNumber - monthStart + 1
diff = int((monthStart - a11) / 29.)
lunarLeap = 0
lunarMonth = diff + 11
if (b11 - a11 > 365):
leapMonthDiff = getLeapMonthOffset(a11, timeZone)
if (diff >= leapMonthDiff):
lunarMonth = diff + 10
if (diff == leapMonthDiff):
lunarLeap = 1
if (lunarMonth > 12):
lunarMonth = lunarMonth - 12
if (lunarMonth >= 11 and diff < 4):
lunarYear -= 1
return [lunarDay, lunarMonth, lunarYear, lunarLeap]
def L2S(lunarD, lunarM, lunarY, lunarLeap, tZ = 7):
# Convert a lunar date to the corresponding solar date.
if (lunarM < 11):
a11 = getLunarMonth11(lunarY - 1, tZ)
b11 = getLunarMonth11(lunarY, tZ)
else:
a11 = getLunarMonth11(lunarY, tZ)
b11 = getLunarMonth11(lunarY + 1, tZ)
k = int(0.5 + (a11 - 2415021.076998695) / 29.530588853)
off = lunarM - 11
if (off < 0):
off += 12
if (b11 - a11 > 365):
leapOff = getLeapMonthOffset(a11, tZ)
leapM = leapOff - 2
if (leapM < 0):
leapM += 12
if (lunarLeap != 0 and lunarM != leapM):
return [0, 0, 0]
elif (lunarLeap != 0 or off >= leapOff):
off += 1
monthStart = getNewMoonDay(k + off, tZ)
return jdToDate(monthStart + lunarD - 1)