📄 module1.bas
字号:
Attribute VB_Name = "Module1"
Private Const PSI_COR As Double = 0#
Private Const EPS_COR As Double = 0#
Private Const T0 As Double = 2451545#
Private Const C_ As Double = 173.14463348
Private Const FN0 As Double = 0
Private Const TWOPI As Double = 6.28318530717959
Private Const RAD2SEC As Double = 206264.806247096
Private Const DEG2RAD As Double = 1.74532925199433E-02
Private Const RAD2DEG As Double = 57.2957795130823
Private Const PI As Double = 3.14159265358979
Private Const ecc As Double = 0.01671022
Private Const rev As Double = 1296000#
Private Const obl As Double = 23.43929111
Private Const factor As Double = 0.0000001
Private Const period As Double = 29.53058853
Private Const errlimit As Double = 1# / 12# 'Two hour 节气时差精度
Dim initialized As Boolean
Dim vterms() As Double, vmoons() As Double, vmonth() As Double
Dim lastnew As Double, lastmon As Double, nextnew As Double, pre_year
Dim c_lng(), clngx(), cobl(), coblx(), nav1(), nav2(), nav(), llng(), llngx(), lobl(), loblx(), pm(), pa(), pl(), pn(), con()
Dim CHshengxiao(), CHtiangan(), CHdizhi(), daysinmonth(), CHjieqi(), daynamesCH(), CHwuxing(), miscchar()
Private Function initial_const()
If Not initialized Then
c_lng = Array(1#, 1#, -1#, -1#, 1#, -1#, -1#, -1#, -1#, -1#, -1#, 1#, -1#, 1#, -1#, 1#, 1#, -1#, -1#, 1#, 1#, -1#, 1#, -1#, 1#, -1#, -1#, -1#, 1#, -1#, -1#, 1#, -1#, 1#, 2#, 2#, 2#, 2#, 2#, -2#, 2#, 2#, 2#, 3#, -3#, -3#, 3#, -3#, 3#, -3#, 3#, 4#, 4#, -4#, -4#, 4#, -4#, 5#, 5#, 5#, -5#, 6#, 6#, 6#, -6#, 6#, -7#, 7#, 7#, -7#, -8#, 10#, 11#, 12#, -13#, -15#, -16#, -16#, 17#, -21#, -22#, 26#, 29#, 29#, -31#, -38#, -46#, 48#, -51#, 58#, 59#, 63#, 63#, -123#, 129#, -158#, -217#, -301#, -386#, -517#, 712#, 1426#, 2062#, -2274#, -13187#, -171996#)
clngx = Array(0.1, -0.1, 0.1, 0.1, 0.1, 0.1, 0.2, -0.2, -0.4, 0.5, 1.2, -1.6, -3.4, -174.2)
cobl = Array(1#, 1#, 1#, -1#, -1#, -1#, 1#, 1#, 1#, 1#, 1#, -1#, 1#, -1#, 1#, -1#, -1#, -1#, 1#, -1#, 1#, 1#, -1#, -2#, -2#, -2#, 3#, 3#, -3#, 3#, 3#, -3#, 3#, 3#, -3#, 3#, 3#, 5#, 6#, 7#, -7#, 7#, -8#, 9#, -10#, -12#, 13#, 16#, -24#, 26#, 27#, 32#, -33#, -53#, 54#, -70#, -95#, 129#, 200#, 224#, -895#, 977#, 5736#, 92025#)
coblx = Array(-0.1, -0.1, 0.3, 0.5, -0.5, -0.6, -3.1, 8.9)
nav1 = Array(0, 0, 1, 0, 2, 1, 3, 0, 4, 0)
nav2 = Array(0, 0, 0, 5, 1, 1, 3, 3, 4, 4)
nav = Array(2, 0, 1, 1, 5, 2, 2, 0, 2, 1, 0, 3, 2, 5, 8, 1, 17, 8, 1, 18, 0, 2, 0, 8, 0, 1, 3, 2, 1, 8, 0, 17, 1, 1, 15, 1, 2, 21, 1, 1, 2, 8, 2, 0, 29, 1, 21, 2, 2, 1, 29, 2, 0, 9, 2, 5, 4, 2, 0, 4, 0, 1, 9, 2, 1, 4, 0, 2, 9, 2, 2, 4, 1, 14, 44, 2, 0, 45, 2, 5, 44, 2, 50, 0, 1, 36, 2, 2, 5, 45, 1, 37, 2, 2, 1, 45, 2, 1, 44, 2, 53, 1, 2, 8, 4, 1, 40, 3, 2, 17, 4, 2, 0, 64, 1, 39, 8, 2, 27, 4, 1, 50, 18, 1, 21, 47, 2, 44, 3, 2, 44, 8, 2, 45, 8, 1, 46, 8, 0, 67, 2, 1, 5, 74, 1, 0, 74, 2, 50, 8, 1, 5, 78, 2, 17, 53, 2, 53, 8, 2, 0, 80, 2, 0, 81, 0, 7, 79, 1, 7, 81, 2, 1, 81, 2, 24, 44, 1, 1, 79, 2, 27, 44)
llng = Array(57, 25, 82, 34, 41, 66, 33, 36, 19, 88, 18, 104, 93, 84, 47, 28, 83, 86, 69, 75, 89, 30, 58, 73, 46, 77, 23, 32, 59, 72, 31, 16, 74, 22, 98, 38, 62, 96, 37, 35, 6, 76, 85, 51, 26, 10, 13, 63, 105, 52, 102, 67, 99, 15, 24, 14, 3, 100, 65, 11, 55, 68, 20, 87, 64, 95, 27, 60, 61, 80, 91, 94, 12, 43, 71, 42, 97, 70, 7, 49, 29, 2, 5, 92, 50, 78, 56, 17, 48, 40, 90, 8, 39, 54, 81, 21, 103, 53, 45, 101, 0, 1, 9, 44, 79, 4)
llngx = Array(81, 7, 97, 0, 39, 40, 9, 44, 45, 103, 101, 79, 1, 4)
lobl = Array(51, 98, 17, 21, 5, 2, 63, 105, 38, 52, 102, 62, 96, 37, 35, 76, 36, 88, 85, 104, 93, 84, 83, 67, 99, 8, 68, 100, 60, 61, 91, 87, 64, 80, 95, 65, 55, 94, 43, 97, 0, 71, 70, 42, 49, 92, 50, 78, 56, 90, 48, 40, 39, 54, 1, 81, 103, 53, 45, 101, 9, 44, 79, 4)
loblx = Array(53, 1, 103, 9, 44, 101, 79, 4)
pm = Array(1047.349, 3497.898, 22903#, 19412.2)
pa = Array(5.203363, 9.53707, 19.191264, 30.068963)
pl = Array(0.60047, 0.871693, 5.466933, 5.32116)
pn = Array(0.001450138, 0.0005841727, 0.0002047497, 0.0001043891)
con = Array(Array(403406#, 0#, 4.721964, 1.621043), Array(195207#, -97597#, 5.937458, 62830.348067), Array(119433#, -59715#, 1.115589, 62830.821524), Array(112392#, -56188#, 5.781616, 62829.634302), _
Array(3891#, -1556#, 5.5474, 125660.5691), Array(2819#, -1126#, 1.512, 125660.9845), Array(1721#, -861#, 4.1897, 62832.4766), Array(0#, 941#, 1.163, 0.813), Array(660#, -264#, 5.415, 125659.31), _
Array(350#, -163#, 4.315, 57533.85), Array(334#, 0#, 4.553, -33.931), Array(314#, 309#, 5.198, 777137.715), Array(268#, -158#, 5.989, 78604.191), Array(242#, 0#, 2.911, 5.412), _
Array(234#, -54#, 1.423, 39302.098), Array(158#, 0#, 0.061, -34.861), Array(132#, -93#, 2.317, 115067.698), Array(129#, -20#, 3.193, 15774.337), Array(114#, 0#, 2.828, 5296.67), _
Array(99#, -47#, 0.52, 58849.27), Array(93#, 0#, 4.65, 5296.11), Array(86#, 0#, 4.35, -3980.7), Array(78#, -33#, 2.75, 52237.69), Array(72#, -32#, 4.5, 55076.47), Array(68#, 0#, 3.23, 261.08), _
Array(64#, -10#, 1.22, 15773.85), Array(46#, -16#, 0.14, 188491.03), Array(38#, 0#, 3.44, -7756.55), Array(37#, 0#, 4.37, 264.89), Array(32#, -24#, 1.14, 117906.27), Array(29#, -13#, 2.84, 55075.75), _
Array(28#, 0#, 5.96, -7961.39), Array(27#, -9#, 5.09, 188489.81), Array(27#, 0#, 1.72, 2132.19), Array(25#, -17#, 2.56, 109771.03), Array(24#, -11#, 1.92, 54868.56), Array(21#, 0#, 0.09, 25443.93), _
Array(21#, 31#, 5.98, -55731.43), Array(20#, -10#, 4.03, 60697.74), Array(18#, 0#, 4.27, 2132.79), Array(17#, -12#, 0.79, 109771.63), Array(14#, 0#, 4.24, -7752.82), Array(13#, -5#, 2.01, 188491.91), _
Array(13#, 0#, 2.65, 207.81), Array(13#, 0#, 4.98, 29424.63), Array(12#, 0#, 0.93, -7.99), Array(10#, 0#, 2.21, 46941.14), Array(10#, 0#, 3.59, -68.29), Array(10#, 0#, 1.5, 21463.25), Array(10#, -9#, 2.55, 157208.4))
' Define the time unit 'u', measured in units of 10000 Julian years from J2000.0.
CHshengxiao = Array("鼠", "牛", "虎", "兔", "龙", "蛇", "马", "羊", "猴", "鸡", "狗", "猪")
CHtiangan = Array("甲", "乙", "丙", "丁", "戊", "己", "庚", "辛", "壬", "癸")
CHdizhi = Array("子", "丑", "寅", "卯", "辰", "巳", "午", "未", "申", "酉", "戌", "亥")
daysinmonth = Array(31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31)
CHjieqi = Array("小寒", "大寒", "立春", "雨水", "惊蛰", "春分", "清明", "谷雨", "立夏", "小满", "芒种", "夏至", "小暑", "大暑", "立秋", "处暑", "白露", "秋分", "寒露", "霜降", "立冬", "小雪", "大雪", "冬至")
daynamesCH = Array("Sun 日", "Mon 一", "Tue 二", "Wed 三", "Thu 四", "Fri 五", "Sat 六")
CHwuxing = Array("甲子乙丑海中金", "丙寅丁卯炉中火", "戊辰已巳大林木", "庚午辛未路旁土", "壬申癸酉剑锋金", "甲戌乙亥山头火", "丙子丁丑涧下水", "戊寅已卯城头土", "庚辰辛巳白蜡金", "壬午癸未杨柳木", "甲申乙酉泉中水", "丙戌丁亥屋上土", "戊子已丑霹雳火", "庚寅辛卯松柏木", "壬辰癸巳长流水", "甲午乙未沙中金", "丙申丁酉山下火", "戊戌已亥平地木", "庚子辛丑壁上土", "壬寅癸卯金箔金", "甲辰乙巳覆灯火", "丙午丁未天河水", "戊申已酉大泽土", "庚戌辛亥钗钏金", "壬子癸丑桑柘木", "甲寅乙卯大溪水", "丙辰丁巳沙中土", "戊午已未天上火", "庚申辛酉石榴木", "壬戌癸亥大海水")
miscchar = Array("初", "一", "二", "三", "四", "五", "六", "七", "八", "九", "十", "廿")
initialized = True
End If
End Function
Private Function earthtilt(ByVal tjd As Double, ByRef mobl As Double, ByRef tobl As Double, ByRef eq As Double, ByRef dpsi As Double, ByRef deps As Double)
Dim tjd_last As Double, T As Double, dp As Double, de As Double
Dim d_psi As Double, d_eps As Double, mean_obliq As Double, true_obliq As Double, eq_eq As Double, args(5) As Double
T = (tjd - T0) / 36525#
If Abs(tjd - tjd_last) > 0.000001 Then Call nutation_angles(T, dp, de)
d_psi = dp + PSI_COR
d_eps = de + EPS_COR
mean_obliq = 84381.448 - 46.815 * T - 0.00059 * T ^ 2 + 0.001813 * T ^ 3
true_obliq = mean_obliq + d_eps
mean_obliq = mean_obliq / 3600#
true_obliq = true_obliq / 3600#
Call fund_args(T, args)
eq_eq = d_psi * Cos(mean_obliq * DEG2RAD) + (0.00264 * Sin(args(4)) + 0.000063 * Sin(2# * args(4)))
eq_eq = eq_eq / 15#
tjd_last = tjd
dpsi = d_psi
deps = d_eps
eq = eq_eq
mobl = mean_obliq
tobl = true_obliq
End Function
Private Function aberration(ByRef pos() As Double, ByRef ve() As Double, ByVal lighttime As Double, ByRef pos2() As Double) As Long
Dim p1mag As Double, vemag As Double, beta As Double, dot As Double, cosd As Double, gammai As Double, p As Double, q As Double, r As Double
If lighttime = 0 Then
p1mag = Sqr(pos(0) ^ 2 + pos(1) ^ 2 + pos(2) ^ 2)
lighttime = p1mag / C_
Else
p1mag = lighttime * C_
End If
vemag = Sqr(ve(0) * ve(0) + ve(1) * ve(1) + ve(2) * ve(2))
beta = vemag / C_
dot = pos(0) * ve(0) + pos(1) * ve(1) + pos(2) * ve(2)
cosd = dot / (p1mag * vemag)
gammai = Sqr(1# - beta ^ 2)
p = beta * cosd
q = (1# + p / (1# + gammai)) * lighttime
r = 1# + p
For j = 0 To 2
pos2(j) = (gammai * pos(j) + q * ve(j)) / r
Next j
End Function
Private Function precession(ByVal tjd1 As Double, ByRef pos() As Double, ByVal tjd2 As Double, ByRef pos2() As Double)
Dim xx As Double, yx As Double, zx As Double, xy As Double, yy As Double, zy As Double, xz As Double, yz As Double, zz As Double, T As Double
Dim t1 As Double, t02 As Double, t2 As Double, t3 As Double, zeta0 As Double, zee As Double, theta As Double, cz0 As Double, sz0 As Double, ct As Double, st As Double, cz As Double, sz As Double
T = (tjd1 - T0) / 36525#
t1 = (tjd2 - tjd1) / 36525#
t02 = T * T
t2 = t1 * t1
t3 = t2 * t1
zeta0 = (2306.2181 + 1.39656 * T - 0.000139 * t02) * t1 + (0.30188 - 0.000344 * T) * t2 + 0.017998 * t3
zee = (2306.2181 + 1.39656 * T - 0.000139 * t02) * t1 + (1.09468 + 0.000066 * T) * t2 + 0.018203 * t3
theta = (2004.3109 - 0.8533 * T - 0.000217 * t02) * t1 + (-0.42665 - 0.000217 * T) * t2 - 0.041833 * t3
zeta0 = zeta0 / RAD2SEC
zee = zee / RAD2SEC
theta = theta / RAD2SEC
cz0 = Cos(zeta0)
sz0 = Sin(zeta0)
ct = Cos(theta)
st = Sin(theta)
cz = Cos(zee)
sz = Sin(zee)
xx = cz0 * ct * cz - sz0 * sz
yx = -sz0 * ct * cz - cz0 * sz
zx = -st * cz
xy = cz0 * ct * sz + sz0 * cz
yy = -sz0 * ct * sz + cz0 * cz
zy = -st * sz
xz = cz0 * st
yz = -sz0 * st
zz = ct
pos2(0) = xx * pos(0) + yx * pos(1) + zx * pos(2)
pos2(1) = xy * pos(0) + yy * pos(1) + zy * pos(2)
pos2(2) = xz * pos(0) + yz * pos(1) + zz * pos(2)
End Function
Private Function nutate(ByVal tjd As Double, ByVal fn As Long, ByRef pos() As Double, ByRef pos2() As Double) As Long
Dim cobm As Double, sobm As Double, cobt As Double, sobt As Double, cpsi As Double, spsi As Double, xx As Double, yx As Double, zx As Double, xy As Double, yy As Double, zy As Double, xz As Double, yz As Double, zz As Double, oblm As Double, oblt As Double, eqeq As Double, psi As Double, eps As Double
Call earthtilt(tjd, oblm, oblt, eqeq, psi, eps)
cobm = Cos(oblm * DEG2RAD)
sobm = Sin(oblm * DEG2RAD)
cobt = Cos(oblt * DEG2RAD)
sobt = Sin(oblt * DEG2RAD)
cpsi = Cos(psi / RAD2SEC)
spsi = Sin(psi / RAD2SEC)
xx = cpsi
yx = -spsi * cobm
zx = -spsi * sobm
xy = spsi * cobt
yy = cpsi * cobm * cobt + sobm * sobt
zy = cpsi * sobm * cobt - cobm * sobt
xz = spsi * sobt
yz = cpsi * cobm * sobt - sobm * cobt
zz = cpsi * sobm * sobt + cobm * cobt
If fn = 0 Then
pos2(0) = xx * pos(0) + yx * pos(1) + zx * pos(2)
pos2(1) = xy * pos(0) + yy * pos(1) + zy * pos(2)
pos2(2) = xz * pos(0) + yz * pos(1) + zz * pos(2)
Else
pos2(0) = xx * pos(0) + xy * pos(1) + xz * pos(2)
pos2(1) = yx * pos(0) + yy * pos(1) + yz * pos(2)
pos2(2) = zx * pos(0) + zy * pos(1) + zz * pos(2)
End If
End Function
Private Function nutation_angles(ByVal T As Double, ByRef longnutation As Double, ByRef obliqnutation As Double) As Long
Dim i As Long, ii As Long, i1 As Long, i2 As Long, iop As Long
Dim A(5) As Double, angle As Double, cc As Double, ss1 As Double, cs As Double, sc As Double, c(106) As Double, s(106) As Double, lng As Double, lngx As Double, obl As Double, oblx As Double
Call fund_args(T, A)
i = 0
For ii = 0 To 9 Step 2
angle = A(nav1(ii)) * CDbl((nav1(1 + ii) + 1))
c(i) = Cos(angle)
s(i) = Sin(angle)
i = i + 1
Next ii
i = 5
For ii = 0 To 9 Step 2
i1 = nav2(ii)
i2 = nav2(1 + ii)
c(i) = c(i1) * c(i2) - s(i1) * s(i2)
s(i) = s(i1) * c(i2) + c(i1) * s(i2)
i = i + 1
Next ii
i = 10
For ii = 0 To 182 Step 3
iop = nav(ii)
i1 = nav(1 + ii)
i2 = nav(2 + ii)
Select Case iop
Case 0
c(i) = c(i1) * c(i2) - s(i1) * s(i2)
s(i) = s(i1) * c(i2) + c(i1) * s(i2)
i = i + 1
Case 1
c(i) = c(i1) * c(i2) + s(i1) * s(i2)
s(i) = s(i1) * c(i2) - c(i1) * s(i2)
i = i + 1
Case 2
cc = c(i1) * c(i2)
ss1 = s(i1) * s(i2)
sc = s(i1) * c(i2)
cs = c(i1) * s(i2)
c(i) = cc - ss1
s(i) = sc + cs
i = i + 1
c(i) = cc + ss1
s(i) = sc - cs
i = i + 1
End Select
If iop = 3 Then Exit For
Next ii
lng = 0#
For i = 0 To 105
lng = lng + c_lng(i) * s(llng(i))
Next i
lngx = 0#
For i = 0 To 13
lngx = lngx + clngx(i) * s(llngx(i))
Next i
obl = 0#
For i = 0 To 63
obl = obl + cobl(i) * c(lobl(i))
Next i
oblx = 0#
For i = 0 To 7
oblx = oblx + coblx(i) * c(loblx(i))
Next i
longnutation = (lng + T * lngx) / 10000#
obliqnutation = (obl + T * oblx) / 10000#
End Function
Private Function fund_args(ByVal T As Double, ByRef A() As Double)
A(0) = 2.35554839354394 + T * (8328.6914228839 + T * (1.51795163555396E-04 + 3.10280755910103E-07 * T))
A(1) = 6.24003593932602 + T * (628.301956024184 + T * (-2.79737494000202E-06 - 5.81776417331443E-08 * T))
A(2) = 1.62790193397196 + T * (8433.46615831845 + T * (-6.42717497046912E-05 + 5.3329504922049E-08 * T))
A(3) = 5.19846951357992 + T * (7771.37714617064 + T * (-3.34085107652581E-05 + 9.21145994108118E-08 * T))
A(4) = 2.18243862436099 + T * (-33.7570459337535 + T * (3.61428599267159E-05 + 3.87850944887629E-08 * T))
For i = 0 To 4
A(i) = fmod(A(i), TWOPI)
If (A(i) < 0#) Then A(i) = A(i) + TWOPI
Next i
End Function
Private Function radec2vector(ByVal ra As Double, ByVal dec As Double, ByVal dist As Double, ByRef vector() As Double)
vector(0) = dist * Cos(DEG2RAD * dec) * Cos(DEG2RAD * 15# * ra)
vector(1) = dist * Cos(DEG2RAD * dec) * Sin(DEG2RAD * 15# * ra)
vector(2) = dist * Sin(DEG2RAD * dec)
End Function
Private Function tdb2tdt(ByVal tdb As Double, ByRef tdtjd As Double, ByRef secdiff As Double)
Dim tdays As Double, M As Double, l As Double, lj As Double, E As Double
tdays = tdb - T0
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -