Skip to content

Commit ecafaed

Browse files
StarainrtStarainrt
committed
更换岁差、章动算法
1 parent b763271 commit ecafaed

23 files changed

Lines changed: 1215 additions & 339 deletions

‎basic/all_test.go‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -27,13 +27,13 @@ func show() {
2727
fmt.Println("当前黄赤交角：", EclipticObliquity(jde-8.0/24.0, true))
2828
fmt.Println("当前日出：", JDE2Date(GetSunRiseTime(jde, 115, 32, 8, 1, 10)))
2929
fmt.Println("当前日落：", JDE2Date(GetSunSetTime(jde, 115, 32, 8, 1, 10)))
30-
fmt.Println("当前晨影 -6：", JDE2Date(GetAsaTime(jde, 115, 32, 8, -6)))
31-
fmt.Println("当前晨影 -12：", JDE2Date(GetAsaTime(jde, 115, 32, 8, -12)))
32-
fmt.Println("当前昏影 -6：", JDE2Date(GetBanTime(jde, 115, 32, 8, -6)))
33-
fmt.Println("当前昏影 -12：", JDE2Date(GetBanTime(jde, 115, 32, 8, -12)))
30+
fmt.Println("当前晨影 -6：", JDE2Date(MorningTwilight(jde, 115, 32, 8, -6)))
31+
fmt.Println("当前晨影 -12：", JDE2Date(MorningTwilight(jde, 115, 32, 8, -12)))
32+
fmt.Println("当前昏影 -6：", JDE2Date(EveningTwilight(jde, 115, 32, 8, -6)))
33+
fmt.Println("当前昏影 -12：", JDE2Date(EveningTwilight(jde, 115, 32, 8, -12)))
3434
fmt.Print("农历：")
3535
fmt.Println(GetLunar(2019, 10, 23, 8.0/24.0))
3636
fmt.Println("当前月出：", JDE2Date(GetMoonRiseTime(jde, 115, 32, 8, 1, 10)))
3737
fmt.Println("当前月落：", JDE2Date(GetMoonSetTime(jde, 115, 32, 8, 1, 10)))
38-
fmt.Println("月相：", MoonLight(jde-8.0/24.0))
38+
fmt.Println("月相：", MoonPhase(jde-8.0/24.0))
3939
}

‎basic/cal_test.go‎

Lines changed: 0 additions & 24 deletions
This file was deleted.

‎basic/calendar.go‎

Lines changed: 6 additions & 109 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,12 @@ import (
99

1010
var defDeltaTFn = DefaultDeltaTv2
1111

12+
// Date2JDE 日期转儒略日
13+
func Date2JDE(date time.Time) float64 {
14+
day := float64(date.Day()) + float64(date.Hour())/24.0 + float64(date.Minute())/24.0/60.0 + float64(date.Second())/24.0/3600.0 + float64(date.Nanosecond())/1000000000.0/3600.0/24.0
15+
return JDECalc(date.Year(), int(date.Month()), day)
16+
}
17+
1218
/*
1319
@name: 儒略日计算
1420
@dec: 计算给定时间的儒略日，1582年改力后为格里高利历，之前为儒略历
@@ -39,61 +45,6 @@ func GetNowJDE() (NowJDE float64) {
3945
return
4046
}
4147

42-
func dt_ext(y, jsd float64) float64 { // 二次曲线外推，用于数值外插
43-
dy := (y - 1820) / 100.00
44-
return -20 + jsd*dy*dy
45-
}
46-
func dt_cal(y float64) float64 { //传入年， 返回世界时UT与原子时（力学时 TD）之差, ΔT = TD - UT
47-
dt_at := []float64{
48-
-4000, 108371.7, -13036.80, 392.000, 0.0000,
49-
-500, 17201.0, -627.82, 16.170, -0.3413,
50-
-150, 12200.6, -346.41, 5.403, -0.1593,
51-
150, 9113.8, -328.13, -1.647, 0.0377,
52-
500, 5707.5, -391.41, 0.915, 0.3145,
53-
900, 2203.4, -283.45, 13.034, -0.1778,
54-
1300, 490.1, -57.35, 2.085, -0.0072,
55-
1600, 120.0, -9.81, -1.532, 0.1403,
56-
1700, 10.2, -0.91, 0.510, -0.0370,
57-
1800, 13.4, -0.72, 0.202, -0.0193,
58-
1830, 7.8, -1.81, 0.416, -0.0247,
59-
1860, 8.3, -0.13, -0.406, 0.0292,
60-
1880, -5.4, 0.32, -0.183, 0.0173,
61-
1900, -2.3, 2.06, 0.169, -0.0135,
62-
1920, 21.2, 1.69, -0.304, 0.0167,
63-
1940, 24.2, 1.22, -0.064, 0.0031,
64-
1960, 33.2, 0.51, 0.231, -0.0109,
65-
1980, 51.0, 1.29, -0.026, 0.0032,
66-
2000, 63.87, 0.1, 0, 0,
67-
2005, 64.7, 0.4, 0, 0, //一次项记为x,则 10x=0.4秒/年*(2015-2005),解得x=0.4
68-
2015, 69,
69-
}
70-
y0 := dt_at[len(dt_at)-2] //表中最后一年
71-
t0 := dt_at[len(dt_at)-1] //表中最后一年的 deltatT
72-
if y >= y0 {
73-
jsd := float64(31) // sjd是y1年之后的加速度估计
74-
// 瑞士星历表jsd=31, NASA网站jsd=32, skmap的jsd=29
75-
if y > y0+100.00 {
76-
return dt_ext(y, jsd)
77-
}
78-
v := dt_ext(y, jsd) //二次曲线外推
79-
dv := dt_ext(y0, jsd) - t0 // ye年的二次外推与te的差
80-
return (v - dv*(y0+100.00-y)/100.00)
81-
}
82-
d := dt_at
83-
var i int
84-
for i = 0; i < len(d); i += 5 {
85-
if float64(y) < d[i+5] {
86-
break
87-
// 判断年所在的区间
88-
}
89-
}
90-
t1 := (y - d[i]) / (d[i+5] - d[i]) * 10.00 //////// 三次插值， 保证精确性
91-
t2 := t1 * t1
92-
t3 := t2 * t1
93-
res := d[i+1] + d[i+2]*t1 + d[i+3]*t2 + d[i+4]*t3
94-
return (res)
95-
}
96-
9748
func DeltaT(date float64, isJDE bool) float64 {
9849
return defDeltaTFn(date, isJDE)
9950
}
@@ -108,50 +59,6 @@ func GetDeltaTFn() func(float64, bool) float64 {
10859
return defDeltaTFn
10960
}
11061

111-
func OldDefaultDeltaT(Date float64, IsJDE bool) (Result float64) { //传入年或儒略日，传出为秒
112-
var Year float64
113-
if IsJDE {
114-
dates := JDE2Date(Date)
115-
Year = float64(dates.Year()) + float64(dates.YearDay())/365.0
116-
} else {
117-
Year = Date
118-
}
119-
if Year < 2100 && Year >= 2010 {
120-
return dt_cal(Year)
121-
}
122-
return DefaultDeltaT(Date, IsJDE)
123-
}
124-
125-
func DefaultDeltaT(Date float64, IsJDE bool) (Result float64) { //传入年或儒略日，传出为秒
126-
var Year float64
127-
if IsJDE {
128-
dates := JDE2Date(Date)
129-
Year = float64(dates.Year()) + float64(dates.YearDay())/365.0
130-
} else {
131-
Year = Date
132-
}
133-
if Year < 2010 {
134-
Result = dt_cal(Year)
135-
return
136-
}
137-
if Year < 2100 && Year >= 2010 {
138-
var t = (Year - 2000.0)
139-
Result = 62.92 + 0.32217*t + 0.005589*t*t
140-
return
141-
}
142-
if Year >= 2100 && Year <= 2150 {
143-
Result = -20 + 32*(((Year-1820.0)/100.0)*((Year-1820.0)/100.0)) - 0.5628*(2150-Year)
144-
return
145-
}
146-
if Year > 2150 {
147-
//tmp=(Year-1820)/100;
148-
//Result= -20 + 32 * tmp*tmp;
149-
Result = dt_cal(Year)
150-
return
151-
}
152-
return
153-
}
154-
15562
func DefaultDeltaTv2(date float64, isJd bool) float64 { //传入年或儒略日，传出为秒
15663
if !isJd {
15764
date = JDECalc(int(date), int((date-math.Floor(date))*12)+1, (date-math.Floor(date))*365.25+1)
@@ -240,7 +147,6 @@ func DeltaTv2(jd float64) float64 {
240147
}
241148

242149
func TD2UT(JDE float64, UT2TD bool) float64 { // true 世界时转力学时CC，false 力学时转世界时VV
243-
244150
Deltat := DeltaT(JDE, true)
245151
if UT2TD {
246152
return JDE + Deltat/3600/24
@@ -249,9 +155,6 @@ func TD2UT(JDE float64, UT2TD bool) float64 { // true 世界时转力学时CC，
249155
}
250156
}
251157

252-
/*
253-
* @name: JDE转日期，输出为数组
254-
*/
255158
func JDE2Date(JD float64) time.Time {
256159
JD = JD + 0.5
257160
Z := float64(int(JD))
@@ -330,12 +233,6 @@ func JDE2DateByZone(JD float64, tz *time.Location, byZone bool) time.Time {
330233
Add(time.Duration(int64(1000000000 * tms))).In(tz)
331234
}
332235

333-
// Date2JDE 日期转儒略日
334-
func Date2JDE(date time.Time) float64 {
335-
day := float64(date.Day()) + float64(date.Hour())/24.0 + float64(date.Minute())/24.0/60.0 + float64(date.Second())/24.0/3600.0 + float64(date.Nanosecond())/1000000000.0/3600.0/24.0
336-
return JDECalc(date.Year(), int(date.Month()), day)
337-
}
338-
339236
func GetLunar(year, month, day int, tz float64) (lyear, lmonth, lday int, leap bool, result string) {
340237
julianDayEpoch := JDECalc(year, month, float64(day))
341238
// 确定农历年份

‎basic/coordinate.go‎

Lines changed: 0 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -60,27 +60,6 @@ func DecToBo(jde, ra, dec float64) float64 {
6060
return ArcSin(sinBo)
6161
}
6262

63-
/*
64-
* 赤道坐标岁差变换st end 为JDE时刻
65-
*/
66-
func ZuoBiaoSuiCha(ra, dec, st, end float64) (float64, float64) {
67-
t := (end - st) / 36525.0
68-
l := (2306.2181*t + 0.30188*t*t + 0.017998*t*t*t) / 3600
69-
z := (2306.2181*t + 1.09468*t*t + 0.018203*t*t*t) / 3600
70-
o := (2004.3109*t - 0.42665*t*t + 0.041833*t*t*t) / 3600
71-
A := Cos(dec) * Sin(ra+l)
72-
B := Cos(o)*Cos(dec)*Cos(ra+l) - Sin(o)*Sin(dec)
73-
C := Sin(o)*Cos(dec)*Cos(ra+l) + Cos(o)*Sin(dec)
74-
ras := math.Atan2(A, B)
75-
ras = ras * 180 / math.Pi
76-
if ras < 0 {
77-
ras += 360
78-
}
79-
ra = ras + z
80-
dec = ArcSin(C)
81-
return ra, dec
82-
}
83-
8463
/*
8564
* 地心坐标转站心坐标，参数分别为，地心赤经赤纬 纬度经度，jde，离地心位置au
8665
*/

‎basic/cst.go‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -156,7 +156,7 @@ func IsXZ(ra, dec, jde float64) string {
156156
if ra >= 360 {
157157
nra -= 360
158158
}
159-
nra, ndec = ZuoBiaoSuiCha(nra, dec, jde, 2451545.0)
159+
nra, ndec = Precess(nra, dec, jde, 2451545.0)
160160
if ra >= 360 && nra < 270 {
161161
nra += 360
162162
}

‎basic/jupiter.go‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -103,7 +103,7 @@ func JupiterApparentLo(JD float64) float64 {
103103
lo = Limit360(lo)
104104
//lo-=GXCLo(lo,bo,JD)/3600;
105105
//bo+=GXCBo(lo,bo,JD);
106-
lo += HJZD(JD)
106+
lo += Nutation2000Bi(JD)
107107
return lo
108108
}
109109

@@ -117,7 +117,7 @@ func JupiterApparentBo(JD float64) float64 {
117117
bo = bo * 180 / math.Pi
118118
//lo+=GXCLo(lo,bo,JD);
119119
//bo+=GXCBo(lo,bo,JD)/3600;
120-
//lo+=HJZD(JD);
120+
//lo+=Nutation2000Bi(JD);
121121
return bo
122122
}
123123

@@ -132,7 +132,7 @@ func JupiterApparentLoBo(JD float64) (float64, float64) {
132132
lo = Limit360(lo)
133133
//lo-=GXCLo(lo,bo,JD)/3600;
134134
//bo+=GXCBo(lo,bo,JD);
135-
lo += HJZD(JD)
135+
lo += Nutation2000Bi(JD)
136136
return lo, bo
137137
}
138138

‎basic/mars.go‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -100,7 +100,7 @@ func MarsApparentLo(JD float64) float64 {
100100
//bo := math.Atan2(z, math.Sqrt(x*x+y*y))
101101
lo = lo * 180.0 / math.Pi
102102
//bo = bo * 180 / math.Pi
103-
lo = Limit360(lo) + HJZD(JD)
103+
lo = Limit360(lo) + Nutation2000Bi(JD)
104104
//lo-=GXCLo(lo,bo,JD)/3600;
105105
//bo+=GXCBo(lo,bo,JD);
106106
return lo
@@ -116,7 +116,7 @@ func MarsApparentBo(JD float64) float64 {
116116
bo = bo * 180 / math.Pi
117117
//lo+=GXCLo(lo,bo,JD);
118118
//bo+=GXCBo(lo,bo,JD)/3600;
119-
//lo+=HJZD(JD);
119+
//lo+=Nutation2000Bi(JD);
120120
return bo
121121
}
122122

@@ -131,7 +131,7 @@ func MarsApparentLoBo(JD float64) (float64, float64) {
131131
lo = Limit360(lo)
132132
//lo -= GXCLo(lo, bo, JD) / 3600
133133
//bo += GXCBo(lo, bo, JD)
134-
lo += HJZD(JD)
134+
lo += Nutation2000Bi(JD)
135135
return lo, bo
136136
}
137137

‎basic/math.go‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
package basic
2+
3+
import "math"
4+
5+
const rad = math.Pi / 180.0
6+
const deg = 180.0 / math.Pi

‎basic/mercury.go‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -96,7 +96,7 @@ func MercuryApparentLo(JD float64) float64 {
9696
lo = Limit360(lo)
9797
//lo-=GXCLo(lo,bo,JD)/3600;
9898
//bo+=GXCBo(lo,bo,JD);
99-
lo += HJZD(JD)
99+
lo += Nutation2000Bi(JD)
100100
return lo
101101
}
102102

@@ -110,7 +110,7 @@ func MercuryApparentBo(JD float64) float64 {
110110
bo = bo * 180 / math.Pi
111111
//lo+=GXCLo(lo,bo,JD);
112112
//bo+=GXCBo(lo,bo,JD)/3600;
113-
//lo+=HJZD(JD);
113+
//lo+=Nutation2000Bi(JD);
114114
return bo
115115
}
116116

@@ -122,7 +122,7 @@ func MercuryApparentLoBo(JD float64) (float64, float64) {
122122
bo := math.Atan2(z, math.Sqrt(x*x+y*y))
123123
lo = lo * 180 / math.Pi
124124
bo = bo * 180 / math.Pi
125-
lo = Limit360(lo) + HJZD(JD)
125+
lo = Limit360(lo) + Nutation2000Bi(JD)
126126
//lo-=GXCLo(lo,bo,JD)/3600;
127127
//bo+=GXCBo(lo,bo,JD);
128128
return lo, bo

‎basic/moon.go‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1034,7 +1034,7 @@ func MoonAway(JD float64) float64 { //'月地距离
10341034
* @name 月球视黄经
10351035
*/
10361036
func MoonApparentLo(JD float64) float64 {
1037-
return MoonTrueLo(JD) + HJZD(JD)
1037+
return MoonTrueLo(JD) + Nutation2000Bi(JD)
10381038
}
10391039

10401040
/*
@@ -1083,12 +1083,12 @@ func MoonApparentDec(JD, lon, lat, tz float64) float64 {
10831083
return ndec
10841084
}
10851085

1086-
func MoonLight(JD float64) float64 {
1086+
func MoonPhase(JD float64) float64 {
10871087
MoonBo := HMoonTrueBo(JD)
10881088
SunLo := HSunApparentLo(JD)
10891089
MoonLo := HMoonApparentLo(JD)
10901090
tmp := Cos(MoonBo) * Cos(SunLo-MoonLo)
1091-
R := RDJL(JD) * 149597870.691
1091+
R := Distance(JD) * 149597870.691
10921092
i := R * Sin(ArcCos(tmp)) / (HMoonAway(JD) - R*tmp)
10931093
i = ArcTan(i)
10941094
if i < 0 {
@@ -1723,7 +1723,7 @@ func HMoonAway(JD float64) float64 { //'月地距离
17231723
* @name 月球视黄经
17241724
*/
17251725
func HMoonApparentLo(JD float64) float64 {
1726-
return HMoonTrueLo(JD) + HJZD(JD)
1726+
return HMoonTrueLo(JD) + Nutation2000Bi(JD)
17271727
}
17281728

17291729
func HMoonTrueRaDec(JD float64) (float64, float64) {

0 commit comments

Comments
 (0)