From 543abcafa5501cb67e18aa0311f15b4e9df7dee8 Mon Sep 17 00:00:00 2001 From: starainrt Date: Thu, 18 Sep 2025 13:16:04 +0800 Subject: [PATCH] =?UTF-8?q?=E6=9B=B4=E6=8D=A2=E5=B2=81=E5=B7=AE=E3=80=81?= =?UTF-8?q?=E7=AB=A0=E5=8A=A8=E7=AE=97=E6=B3=95?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- basic/all_test.go | 10 +- basic/cal_test.go | 24 -- basic/calendar.go | 115 +------- basic/coordinate.go | 21 -- basic/cst.go | 2 +- basic/jupiter.go | 6 +- basic/mars.go | 6 +- basic/math.go | 6 + basic/mercury.go | 6 +- basic/moon.go | 8 +- basic/neptune.go | 6 +- basic/nutation.go | 702 ++++++++++++++++++++++++++++++++++++++++++++ basic/precess.go | 309 +++++++++++++++++++ basic/saturn.go | 6 +- basic/star.go | 66 ++++- basic/stardat.go | 2 +- basic/sun.go | 194 ++++-------- basic/sun_test.go | 10 +- basic/uranus.go | 6 +- basic/venus.go | 6 +- moon/moon.go | 2 +- sun/sun.go | 32 +- tools/math.go | 9 + 23 files changed, 1215 insertions(+), 339 deletions(-) delete mode 100644 basic/cal_test.go create mode 100644 basic/math.go create mode 100644 basic/nutation.go create mode 100644 basic/precess.go diff --git a/basic/all_test.go b/basic/all_test.go index bfca415..1c6dae0 100644 --- a/basic/all_test.go +++ b/basic/all_test.go @@ -27,13 +27,13 @@ func show() { fmt.Println("当前黄赤交角：", EclipticObliquity(jde-8.0/24.0, true)) fmt.Println("当前日出：", JDE2Date(GetSunRiseTime(jde, 115, 32, 8, 1, 10))) fmt.Println("当前日落：", JDE2Date(GetSunSetTime(jde, 115, 32, 8, 1, 10))) - fmt.Println("当前晨影 -6：", JDE2Date(GetAsaTime(jde, 115, 32, 8, -6))) - fmt.Println("当前晨影 -12：", JDE2Date(GetAsaTime(jde, 115, 32, 8, -12))) - fmt.Println("当前昏影 -6：", JDE2Date(GetBanTime(jde, 115, 32, 8, -6))) - fmt.Println("当前昏影 -12：", JDE2Date(GetBanTime(jde, 115, 32, 8, -12))) + fmt.Println("当前晨影 -6：", JDE2Date(MorningTwilight(jde, 115, 32, 8, -6))) + fmt.Println("当前晨影 -12：", JDE2Date(MorningTwilight(jde, 115, 32, 8, -12))) + fmt.Println("当前昏影 -6：", JDE2Date(EveningTwilight(jde, 115, 32, 8, -6))) + fmt.Println("当前昏影 -12：", JDE2Date(EveningTwilight(jde, 115, 32, 8, -12))) fmt.Print("农历：") fmt.Println(GetLunar(2019, 10, 23, 8.0/24.0)) fmt.Println("当前月出：", JDE2Date(GetMoonRiseTime(jde, 115, 32, 8, 1, 10))) fmt.Println("当前月落：", JDE2Date(GetMoonSetTime(jde, 115, 32, 8, 1, 10))) - fmt.Println("月相：", MoonLight(jde-8.0/24.0)) + fmt.Println("月相：", MoonPhase(jde-8.0/24.0)) } diff --git a/basic/cal_test.go b/basic/cal_test.go deleted file mode 100644 index e66cd7c..0000000 --- a/basic/cal_test.go +++ /dev/null @@ -1,24 +0,0 @@ -package basic - -import ( - "fmt" - "testing" - "time" -) - -func Test_JDECALC(t *testing.T) { - i := 10 - for i > 0 { - fmt.Printf("%.18f\n", GetNowJDE()) - time.Sleep(time.Second) - i-- - } -} - -func Test_TDUT(t *testing.T) { - fmt.Printf("%.18f\n", DeltaT(2119.5, false)) -} - -func Test_JDE2(t *testing.T) { - fmt.Println(JDE2Date(2458868.500)) -} diff --git a/basic/calendar.go b/basic/calendar.go index e83b2dd..6f9f1f9 100644 --- a/basic/calendar.go +++ b/basic/calendar.go @@ -9,6 +9,12 @@ import ( var defDeltaTFn = DefaultDeltaTv2 +// Date2JDE 日期转儒略日 +func Date2JDE(date time.Time) float64 { + 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 + return JDECalc(date.Year(), int(date.Month()), day) +} + /* @name: 儒略日计算 @dec: 计算给定时间的儒略日，1582年改力后为格里高利历，之前为儒略历 @@ -39,61 +45,6 @@ func GetNowJDE() (NowJDE float64) { return } -func dt_ext(y, jsd float64) float64 { // 二次曲线外推，用于数值外插 - dy := (y - 1820) / 100.00 - return -20 + jsd*dy*dy -} -func dt_cal(y float64) float64 { //传入年， 返回世界时UT与原子时（力学时 TD）之差, ΔT = TD - UT - dt_at := []float64{ - -4000, 108371.7, -13036.80, 392.000, 0.0000, - -500, 17201.0, -627.82, 16.170, -0.3413, - -150, 12200.6, -346.41, 5.403, -0.1593, - 150, 9113.8, -328.13, -1.647, 0.0377, - 500, 5707.5, -391.41, 0.915, 0.3145, - 900, 2203.4, -283.45, 13.034, -0.1778, - 1300, 490.1, -57.35, 2.085, -0.0072, - 1600, 120.0, -9.81, -1.532, 0.1403, - 1700, 10.2, -0.91, 0.510, -0.0370, - 1800, 13.4, -0.72, 0.202, -0.0193, - 1830, 7.8, -1.81, 0.416, -0.0247, - 1860, 8.3, -0.13, -0.406, 0.0292, - 1880, -5.4, 0.32, -0.183, 0.0173, - 1900, -2.3, 2.06, 0.169, -0.0135, - 1920, 21.2, 1.69, -0.304, 0.0167, - 1940, 24.2, 1.22, -0.064, 0.0031, - 1960, 33.2, 0.51, 0.231, -0.0109, - 1980, 51.0, 1.29, -0.026, 0.0032, - 2000, 63.87, 0.1, 0, 0, - 2005, 64.7, 0.4, 0, 0, //一次项记为x,则 10x=0.4秒/年*(2015-2005),解得x=0.4 - 2015, 69, - } - y0 := dt_at[len(dt_at)-2] //表中最后一年 - t0 := dt_at[len(dt_at)-1] //表中最后一年的 deltatT - if y >= y0 { - jsd := float64(31) // sjd是y1年之后的加速度估计 - // 瑞士星历表jsd=31, NASA网站jsd=32, skmap的jsd=29 - if y > y0+100.00 { - return dt_ext(y, jsd) - } - v := dt_ext(y, jsd) //二次曲线外推 - dv := dt_ext(y0, jsd) - t0 // ye年的二次外推与te的差 - return (v - dv*(y0+100.00-y)/100.00) - } - d := dt_at - var i int - for i = 0; i < len(d); i += 5 { - if float64(y) < d[i+5] { - break - // 判断年所在的区间 - } - } - t1 := (y - d[i]) / (d[i+5] - d[i]) * 10.00 //////// 三次插值， 保证精确性 - t2 := t1 * t1 - t3 := t2 * t1 - res := d[i+1] + d[i+2]*t1 + d[i+3]*t2 + d[i+4]*t3 - return (res) -} - func DeltaT(date float64, isJDE bool) float64 { return defDeltaTFn(date, isJDE) } @@ -108,50 +59,6 @@ func GetDeltaTFn() func(float64, bool) float64 { return defDeltaTFn } -func OldDefaultDeltaT(Date float64, IsJDE bool) (Result float64) { //传入年或儒略日，传出为秒 - var Year float64 - if IsJDE { - dates := JDE2Date(Date) - Year = float64(dates.Year()) + float64(dates.YearDay())/365.0 - } else { - Year = Date - } - if Year < 2100 && Year >= 2010 { - return dt_cal(Year) - } - return DefaultDeltaT(Date, IsJDE) -} - -func DefaultDeltaT(Date float64, IsJDE bool) (Result float64) { //传入年或儒略日，传出为秒 - var Year float64 - if IsJDE { - dates := JDE2Date(Date) - Year = float64(dates.Year()) + float64(dates.YearDay())/365.0 - } else { - Year = Date - } - if Year < 2010 { - Result = dt_cal(Year) - return - } - if Year < 2100 && Year >= 2010 { - var t = (Year - 2000.0) - Result = 62.92 + 0.32217*t + 0.005589*t*t - return - } - if Year >= 2100 && Year <= 2150 { - Result = -20 + 32*(((Year-1820.0)/100.0)*((Year-1820.0)/100.0)) - 0.5628*(2150-Year) - return - } - if Year > 2150 { - //tmp=(Year-1820)/100; - //Result= -20 + 32 * tmp*tmp; - Result = dt_cal(Year) - return - } - return -} - func DefaultDeltaTv2(date float64, isJd bool) float64 { //传入年或儒略日，传出为秒 if !isJd { 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 { } func TD2UT(JDE float64, UT2TD bool) float64 { // true 世界时转力学时CC，false 力学时转世界时VV - Deltat := DeltaT(JDE, true) if UT2TD { return JDE + Deltat/3600/24 @@ -249,9 +155,6 @@ func TD2UT(JDE float64, UT2TD bool) float64 { // true 世界时转力学时CC， } } -/* - * @name: JDE转日期，输出为数组 - */ func JDE2Date(JD float64) time.Time { JD = JD + 0.5 Z := float64(int(JD)) @@ -330,12 +233,6 @@ func JDE2DateByZone(JD float64, tz *time.Location, byZone bool) time.Time { Add(time.Duration(int64(1000000000 * tms))).In(tz) } -// Date2JDE 日期转儒略日 -func Date2JDE(date time.Time) float64 { - 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 - return JDECalc(date.Year(), int(date.Month()), day) -} - func GetLunar(year, month, day int, tz float64) (lyear, lmonth, lday int, leap bool, result string) { julianDayEpoch := JDECalc(year, month, float64(day)) // 确定农历年份 diff --git a/basic/coordinate.go b/basic/coordinate.go index 6d6f0a5..a09f9f9 100644 --- a/basic/coordinate.go +++ b/basic/coordinate.go @@ -60,27 +60,6 @@ func DecToBo(jde, ra, dec float64) float64 { return ArcSin(sinBo) } -/* - * 赤道坐标岁差变换st end 为JDE时刻 - */ -func ZuoBiaoSuiCha(ra, dec, st, end float64) (float64, float64) { - t := (end - st) / 36525.0 - l := (2306.2181*t + 0.30188*t*t + 0.017998*t*t*t) / 3600 - z := (2306.2181*t + 1.09468*t*t + 0.018203*t*t*t) / 3600 - o := (2004.3109*t - 0.42665*t*t + 0.041833*t*t*t) / 3600 - A := Cos(dec) * Sin(ra+l) - B := Cos(o)*Cos(dec)*Cos(ra+l) - Sin(o)*Sin(dec) - C := Sin(o)*Cos(dec)*Cos(ra+l) + Cos(o)*Sin(dec) - ras := math.Atan2(A, B) - ras = ras * 180 / math.Pi - if ras < 0 { - ras += 360 - } - ra = ras + z - dec = ArcSin(C) - return ra, dec -} - /* * 地心坐标转站心坐标，参数分别为，地心赤经赤纬 纬度经度，jde，离地心位置au */ diff --git a/basic/cst.go b/basic/cst.go index 685f5f5..36d2437 100644 --- a/basic/cst.go +++ b/basic/cst.go @@ -156,7 +156,7 @@ func IsXZ(ra, dec, jde float64) string { if ra >= 360 { nra -= 360 } - nra, ndec = ZuoBiaoSuiCha(nra, dec, jde, 2451545.0) + nra, ndec = Precess(nra, dec, jde, 2451545.0) if ra >= 360 && nra < 270 { nra += 360 } diff --git a/basic/jupiter.go b/basic/jupiter.go index bc26cb8..bc0e85e 100644 --- a/basic/jupiter.go +++ b/basic/jupiter.go @@ -103,7 +103,7 @@ func JupiterApparentLo(JD float64) float64 { lo = Limit360(lo) //lo-=GXCLo(lo,bo,JD)/3600; //bo+=GXCBo(lo,bo,JD); - lo += HJZD(JD) + lo += Nutation2000Bi(JD) return lo } @@ -117,7 +117,7 @@ func JupiterApparentBo(JD float64) float64 { bo = bo * 180 / math.Pi //lo+=GXCLo(lo,bo,JD); //bo+=GXCBo(lo,bo,JD)/3600; - //lo+=HJZD(JD); + //lo+=Nutation2000Bi(JD); return bo } @@ -132,7 +132,7 @@ func JupiterApparentLoBo(JD float64) (float64, float64) { lo = Limit360(lo) //lo-=GXCLo(lo,bo,JD)/3600; //bo+=GXCBo(lo,bo,JD); - lo += HJZD(JD) + lo += Nutation2000Bi(JD) return lo, bo } diff --git a/basic/mars.go b/basic/mars.go index eecf8e6..9323cd5 100644 --- a/basic/mars.go +++ b/basic/mars.go @@ -100,7 +100,7 @@ func MarsApparentLo(JD float64) float64 { //bo := math.Atan2(z, math.Sqrt(x*x+y*y)) lo = lo * 180.0 / math.Pi //bo = bo * 180 / math.Pi - lo = Limit360(lo) + HJZD(JD) + lo = Limit360(lo) + Nutation2000Bi(JD) //lo-=GXCLo(lo,bo,JD)/3600; //bo+=GXCBo(lo,bo,JD); return lo @@ -116,7 +116,7 @@ func MarsApparentBo(JD float64) float64 { bo = bo * 180 / math.Pi //lo+=GXCLo(lo,bo,JD); //bo+=GXCBo(lo,bo,JD)/3600; - //lo+=HJZD(JD); + //lo+=Nutation2000Bi(JD); return bo } @@ -131,7 +131,7 @@ func MarsApparentLoBo(JD float64) (float64, float64) { lo = Limit360(lo) //lo -= GXCLo(lo, bo, JD) / 3600 //bo += GXCBo(lo, bo, JD) - lo += HJZD(JD) + lo += Nutation2000Bi(JD) return lo, bo } diff --git a/basic/math.go b/basic/math.go new file mode 100644 index 0000000..9001b20 --- /dev/null +++ b/basic/math.go @@ -0,0 +1,6 @@ +package basic + +import "math" + +const rad = math.Pi / 180.0 +const deg = 180.0 / math.Pi diff --git a/basic/mercury.go b/basic/mercury.go index cfcb2c2..9ad01ad 100644 --- a/basic/mercury.go +++ b/basic/mercury.go @@ -96,7 +96,7 @@ func MercuryApparentLo(JD float64) float64 { lo = Limit360(lo) //lo-=GXCLo(lo,bo,JD)/3600; //bo+=GXCBo(lo,bo,JD); - lo += HJZD(JD) + lo += Nutation2000Bi(JD) return lo } @@ -110,7 +110,7 @@ func MercuryApparentBo(JD float64) float64 { bo = bo * 180 / math.Pi //lo+=GXCLo(lo,bo,JD); //bo+=GXCBo(lo,bo,JD)/3600; - //lo+=HJZD(JD); + //lo+=Nutation2000Bi(JD); return bo } @@ -122,7 +122,7 @@ func MercuryApparentLoBo(JD float64) (float64, float64) { bo := math.Atan2(z, math.Sqrt(x*x+y*y)) lo = lo * 180 / math.Pi bo = bo * 180 / math.Pi - lo = Limit360(lo) + HJZD(JD) + lo = Limit360(lo) + Nutation2000Bi(JD) //lo-=GXCLo(lo,bo,JD)/3600; //bo+=GXCBo(lo,bo,JD); return lo, bo diff --git a/basic/moon.go b/basic/moon.go index 039d631..0c77b54 100644 --- a/basic/moon.go +++ b/basic/moon.go @@ -1034,7 +1034,7 @@ func MoonAway(JD float64) float64 { //'月地距离 * @name 月球视黄经 */ func MoonApparentLo(JD float64) float64 { - return MoonTrueLo(JD) + HJZD(JD) + return MoonTrueLo(JD) + Nutation2000Bi(JD) } /* @@ -1083,12 +1083,12 @@ func MoonApparentDec(JD, lon, lat, tz float64) float64 { return ndec } -func MoonLight(JD float64) float64 { +func MoonPhase(JD float64) float64 { MoonBo := HMoonTrueBo(JD) SunLo := HSunApparentLo(JD) MoonLo := HMoonApparentLo(JD) tmp := Cos(MoonBo) * Cos(SunLo-MoonLo) - R := RDJL(JD) * 149597870.691 + R := Distance(JD) * 149597870.691 i := R * Sin(ArcCos(tmp)) / (HMoonAway(JD) - R*tmp) i = ArcTan(i) if i < 0 { @@ -1723,7 +1723,7 @@ func HMoonAway(JD float64) float64 { //'月地距离 * @name 月球视黄经 */ func HMoonApparentLo(JD float64) float64 { - return HMoonTrueLo(JD) + HJZD(JD) + return HMoonTrueLo(JD) + Nutation2000Bi(JD) } func HMoonTrueRaDec(JD float64) (float64, float64) { diff --git a/basic/neptune.go b/basic/neptune.go index 8a6f879..5cbc6f5 100644 --- a/basic/neptune.go +++ b/basic/neptune.go @@ -103,7 +103,7 @@ func NeptuneApparentLo(JD float64) float64 { lo = Limit360(lo) //lo-=GXCLo(lo,bo,JD)/3600; //bo+=GXCBo(lo,bo,JD); - lo += HJZD(JD) + lo += Nutation2000Bi(JD) return lo } @@ -117,7 +117,7 @@ func NeptuneApparentBo(JD float64) float64 { bo = bo * 180 / math.Pi //lo+=GXCLo(lo,bo,JD); //bo+=GXCBo(lo,bo,JD)/3600; - //lo+=HJZD(JD); + //lo+=Nutation2000Bi(JD); return bo } @@ -132,7 +132,7 @@ func NeptuneApparentLoBo(JD float64) (float64, float64) { lo = Limit360(lo) //lo-=GXCLo(lo,bo,JD)/3600; //bo+=GXCBo(lo,bo,JD); - lo += HJZD(JD) + lo += Nutation2000Bi(JD) return lo, bo } diff --git a/basic/nutation.go b/basic/nutation.go new file mode 100644 index 0000000..3afa827 --- /dev/null +++ b/basic/nutation.go @@ -0,0 +1,702 @@ +package basic + +import ( + "math" +) + +/* +黄赤交角、nutation==true时，计算交角章动 +*/ +func EclipticObliquity(jde float64, nutation bool) float64 { + eps := Obliquity1980(jde) + if nutation { + eps += Nutation2000Bs(jde) + } + return eps +} + +func Sita(JD float64) float64 { + return EclipticObliquity(JD, true) +} + +// 黄经章动 1980 +func Nutation1980i(jd float64) float64 { // '黄经章动 + res, _ := Nutation1980(jd) + return res +} + +// 交角章动1980 +func Nutation1980s(jd float64) float64 { //交角章动 + _, res := Nutation1980(jd) + return res +} + +// 黄经章动 2000B +func Nutation2000Bi(jd float64) float64 { // '黄经章动 + res, _ := Nutation2000B(jd) + return res +} + +// 交角章动2000B +func Nutation2000Bs(jd float64) float64 { //交角章动 + _, res := Nutation2000B(jd) + return res +} + +// 定义 IAU 1980 章动系数 +var coefficientsNut1980 = [][9]float64{ + {0, 0, 0, 0, 1, -171996.0, -174.2, 92025.0, 8.9}, + {0, 0, 0, 0, 2, 2062.0, 0.2, -895.0, 0.5}, + {-2, 0, 2, 0, 1, 46.0, 0.0, -24.0, 0.0}, + {2, 0, -2, 0, 0, 11.0, 0.0, 0.0, 0.0}, + {-2, 0, 2, 0, 2, -3.0, 0.0, 1.0, 0.0}, + {1, -1, 0, -1, 0, -3.0, 0.0, 0.0, 0.0}, + {0, -2, 2, -2, 1, -2.0, 0.0, 1.0, 0.0}, + {2, 0, -2, 0, 1, 1.0, 0.0, 0.0, 0.0}, + {0, 0, 2, -2, 2, -13187.0, -1.6, 5736.0, -3.1}, + {0, 1, 0, 0, 0, 1426.0, -3.4, 54.0, -0.1}, + {0, 1, 2, -2, 2, -517.0, 1.2, 224.0, -0.6}, + {0, -1, 2, -2, 2, 217.0, -0.5, -95.0, 0.3}, + {0, 0, 2, -2, 1, 129.0, 0.1, -70.0, 0.0}, + {2, 0, 0, -2, 0, 48.0, 0.0, 1.0, 0.0}, + {0, 0, 2, -2, 0, -22.0, 0.0, 0.0, 0.0}, + {0, 2, 0, 0, 0, 17.0, -0.1, 0.0, 0.0}, + {0, 1, 0, 0, 1, -15.0, 0.0, 9.0, 0.0}, + {0, 2, 2, -2, 2, -16.0, 0.1, 7.0, 0.0}, + {0, -1, 0, 0, 1, -12.0, 0.0, 6.0, 0.0}, + {-2, 0, 0, 2, 1, -6.0, 0.0, 3.0, 0.0}, + {0, -1, 2, -2, 1, -5.0, 0.0, 3.0, 0.0}, + {2, 0, 0, -2, 1, 4.0, 0.0, -2.0, 0.0}, + {0, 1, 2, -2, 1, 4.0, 0.0, -2.0, 0.0}, + {1, 0, 0, -1, 0, -4.0, 0.0, 0.0, 0.0}, + {2, 1, 0, -2, 0, 1.0, 0.0, 0.0, 0.0}, + {0, 0, -2, 2, 1, 1.0, 0.0, 0.0, 0.0}, + {0, 1, -2, 2, 0, -1.0, 0.0, 0.0, 0.0}, + {0, 1, 0, 0, 2, 1.0, 0.0, 0.0, 0.0}, + {-1, 0, 0, 1, 1, 1.0, 0.0, 0.0, 0.0}, + {0, 1, 2, -2, 0, -1.0, 0.0, 0.0, 0.0}, + {0, 0, 2, 0, 2, -2274.0, -0.2, 977.0, -0.5}, + {1, 0, 0, 0, 0, 712.0, 0.1, -7.0, 0.0}, + {0, 0, 2, 0, 1, -386.0, -0.4, 200.0, 0.0}, + {1, 0, 2, 0, 2, -301.0, 0.0, 129.0, -0.1}, + {1, 0, 0, -2, 0, -158.0, 0.0, -1.0, 0.0}, + {-1, 0, 2, 0, 2, 123.0, 0.0, -53.0, 0.0}, + {0, 0, 0, 2, 0, 63.0, 0.0, -2.0, 0.0}, + {1, 0, 0, 0, 1, 63.0, 0.1, -33.0, 0.0}, + {-1, 0, 0, 0, 1, -58.0, -0.1, 32.0, 0.0}, + {-1, 0, 2, 2, 2, -59.0, 0.0, 26.0, 0.0}, + {1, 0, 2, 0, 1, -51.0, 0.0, 27.0, 0.0}, + {0, 0, 2, 2, 2, -38.0, 0.0, 16.0, 0.0}, + {2, 0, 0, 0, 0, 29.0, 0.0, -1.0, 0.0}, + {1, 0, 2, -2, 2, 29.0, 0.0, -12.0, 0.0}, + {2, 0, 2, 0, 2, -31.0, 0.0, 13.0, 0.0}, + {0, 0, 2, 0, 0, 26.0, 0.0, -1.0, 0.0}, + {-1, 0, 2, 0, 1, 21.0, 0.0, -10.0, 0.0}, + {-1, 0, 0, 2, 1, 16.0, 0.0, -8.0, 0.0}, + {1, 0, 0, -2, 1, -13.0, 0.0, 7.0, 0.0}, + {-1, 0, 2, 2, 1, -10.0, 0.0, 5.0, 0.0}, + {1, 1, 0, -2, 0, -7.0, 0.0, 0.0, 0.0}, + {0, 1, 2, 0, 2, 7.0, 0.0, -3.0, 0.0}, + {0, -1, 2, 0, 2, -7.0, 0.0, 3.0, 0.0}, + {1, 0, 2, 2, 2, -8.0, 0.0, 3.0, 0.0}, + {1, 0, 0, 2, 0, 6.0, 0.0, 0.0, 0.0}, + {2, 0, 2, -2, 2, 6.0, 0.0, -3.0, 0.0}, + {0, 0, 0, 2, 1, -6.0, 0.0, 3.0, 0.0}, + {0, 0, 2, 2, 1, -7.0, 0.0, 3.0, 0.0}, + {1, 0, 2, -2, 1, 6.0, 0.0, -3.0, 0.0}, + {0, 0, 0, -2, 1, -5.0, 0.0, 3.0, 0.0}, + {1, -1, 0, 0, 0, 5.0, 0.0, 0.0, 0.0}, + {2, 0, 2, 0, 1, -5.0, 0.0, 3.0, 0.0}, + {0, 1, 0, -2, 0, -4.0, 0.0, 0.0, 0.0}, + {1, 0, -2, 0, 0, 4.0, 0.0, 0.0, 0.0}, + {0, 0, 0, 1, 0, -4.0, 0.0, 0.0, 0.0}, + {1, 1, 0, 0, 0, -3.0, 0.0, 0.0, 0.0}, + {1, 0, 2, 0, 0, 3.0, 0.0, 0.0, 0.0}, + {1, -1, 2, 0, 2, -3.0, 0.0, 1.0, 0.0}, + {-1, -1, 2, 2, 2, -3.0, 0.0, 1.0, 0.0}, + {-2, 0, 0, 0, 1, -2.0, 0.0, 1.0, 0.0}, + {3, 0, 2, 0, 2, -3.0, 0.0, 1.0, 0.0}, + {0, -1, 2, 2, 2, -3.0, 0.0, 1.0, 0.0}, + {1, 1, 2, 0, 2, 2.0, 0.0, -1.0, 0.0}, + {-1, 0, 2, -2, 1, -2.0, 0.0, 1.0, 0.0}, + {2, 0, 0, 0, 1, 2.0, 0.0, -1.0, 0.0}, + {1, 0, 0, 0, 2, -2.0, 0.0, 1.0, 0.0}, + {3, 0, 0, 0, 0, 2.0, 0.0, 0.0, 0.0}, + {0, 0, 2, 1, 2, 2.0, 0.0, -1.0, 0.0}, + {-1, 0, 0, 0, 2, 1.0, 0.0, -1.0, 0.0}, + {1, 0, 0, -4, 0, -1.0, 0.0, 0.0, 0.0}, + {-2, 0, 2, 2, 2, 1.0, 0.0, -1.0, 0.0}, + {-1, 0, 2, 4, 2, -2.0, 0.0, 1.0, 0.0}, + {2, 0, 0, -4, 0, -1.0, 0.0, 0.0, 0.0}, + {1, 1, 2, -2, 2, 1.0, 0.0, -1.0, 0.0}, + {1, 0, 2, 2, 1, -1.0, 0.0, 1.0, 0.0}, + {-2, 0, 2, 4, 2, -1.0, 0.0, 1.0, 0.0}, + {-1, 0, 4, 0, 2, 1.0, 0.0, 0.0, 0.0}, + {1, -1, 0, -2, 0, 1.0, 0.0, 0.0, 0.0}, + {2, 0, 2, -2, 1, 1.0, 0.0, -1.0, 0.0}, + {2, 0, 2, 2, 2, -1.0, 0.0, 0.0, 0.0}, + {1, 0, 0, 2, 1, -1.0, 0.0, 0.0, 0.0}, + {0, 0, 4, -2, 2, 1.0, 0.0, 0.0, 0.0}, + {3, 0, 2, -2, 2, 1.0, 0.0, 0.0, 0.0}, + {1, 0, 2, -2, 0, -1.0, 0.0, 0.0, 0.0}, + {0, 1, 2, 0, 1, 1.0, 0.0, 0.0, 0.0}, + {-1, -1, 0, 2, 1, 1.0, 0.0, 0.0, 0.0}, + {0, 0, -2, 0, 1, -1.0, 0.0, 0.0, 0.0}, + {0, 0, 2, -1, 2, -1.0, 0.0, 0.0, 0.0}, + {0, 1, 0, 2, 0, -1.0, 0.0, 0.0, 0.0}, + {1, 0, -2, -2, 0, -1.0, 0.0, 0.0, 0.0}, + {0, -1, 2, 0, 1, -1.0, 0.0, 0.0, 0.0}, + {1, 1, 0, -2, 1, -1.0, 0.0, 0.0, 0.0}, + {1, 0, -2, 2, 0, -1.0, 0.0, 0.0, 0.0}, + {2, 0, 0, 2, 0, 1.0, 0.0, 0.0, 0.0}, + {0, 0, 2, 4, 2, -1.0, 0.0, 0.0, 0.0}, + {0, 1, 0, 1, 0, 1.0, 0.0, 0.0, 0.0}, +} + +// fmod 函数实现浮点数取模 +func fmod(f, m float64) float64 { + return f - m*math.Floor(f/m) +} + +// frac 函数获取小数部分 +func frac(f float64) float64 { + return f - math.Floor(f) +} + +// Obliquity1980 计算黄赤交角，单位为度 +// 公式 3.222-1 +func Obliquity1980(jd float64) float64 { + T := (jd - 2451545.0) / 36525.0 + as2r := ((1.0 / 3600.0) * math.Pi) / 180.0 + eps := 84381.448 - 46.8150*T - 0.00059*T*T + 0.001813*T*T*T // 84381.448 = 23d26'21.448 转换为角秒 + + return math.Mod(eps*as2r, 2*math.Pi) * deg +} + +// Nutation1980 计算 IAU 1980 章动模型 +// 返回黄经章动 (dpsi) 和交角章动 (deps)，单位为度 +func Nutation1980(jd float64) (float64, float64) { + t := (jd - 2451545.0) / 36525.0 + + twoPI := 2 * math.Pi + as2r := ((1.0 / 3600.0) * math.Pi) / 180.0 + + // 表 3.222.2 - 计算章动参数 + l := fmod((485866.733+(715922.633+(31.310+0.064*t)*t)*t)*as2r+frac(1325.0*t)*twoPI, twoPI) + lp := fmod((1287099.804+(1292581.224+(-0.577-0.012*t)*t)*t)*as2r+frac(99.0*t)*twoPI, twoPI) + F := fmod((335778.877+(295263.137+(-13.257+0.011*t)*t)*t)*as2r+frac(1342.0*t)*twoPI, twoPI) + D := fmod((1072261.307+(1105601.328+(-6.891+0.019*t)*t)*t)*as2r+frac(1236.0*t)*twoPI, twoPI) + O := fmod((450160.280+(-482890.539+(7.455+0.008*t)*t)*t)*as2r+frac(-5.0*t)*twoPI, twoPI) + + deps := 0.0 + dpsi := 0.0 + + // 公式 3.222-6 - 计算章动 + for i := len(coefficientsNut1980) - 1; i >= 0; i-- { + sumargs := coefficientsNut1980[i][0]*l + + coefficientsNut1980[i][1]*lp + + coefficientsNut1980[i][2]*F + + coefficientsNut1980[i][3]*D + + coefficientsNut1980[i][4]*O + + deps += math.Cos(sumargs) * (coefficientsNut1980[i][7] + coefficientsNut1980[i][8]*t) + dpsi += math.Sin(sumargs) * (coefficientsNut1980[i][5] + coefficientsNut1980[i][6]*t) + } + + deps = (deps * as2r) / 10000.0 + dpsi = (dpsi * as2r) / 10000.0 + + return dpsi * deg, deps * deg +} + +// Nutation2000B 计算 IAU 2000B 章动模型 +// 返回交角章动 (de) 和黄经章动 (dp)，单位为度 +func Nutation2000B(jd float64) (float64, float64) { + // 常量定义 + as2r := 4.848136811095359935899141e-6 // 角秒到弧度的转换因子 + twopi := 6.283185307179586476925287 // 2π + + T := (jd - 2451545) / 36525 + L := math.Mod(485868.249036+1717915923.2178*T, 1296000.0) * as2r + Lp := math.Mod(1287104.79305+129596581.0481*T, 1296000.0) * as2r + F := math.Mod(335779.526232+1739527262.8478*T, 1296000.0) * as2r + D := math.Mod(1072260.70369+1602961601.2090*T, 1296000.0) * as2r + Om := math.Mod(450160.398036-6962890.5431*T, 1296000.0) * as2r + + dp := 0.0 + de := 0.0 + var arg, sinarg, cosarg float64 + + // 以下是 IAU 2000B 章动模型的各项计算 + // 每行对应一个章动项 + arg = math.Mod(L+Lp+2*F+-2*D+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 1290 * sinarg + de += -556 * cosarg + + arg = math.Mod(-1*L+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 1405*sinarg + 4*cosarg + de += -610*cosarg + 2*sinarg + + arg = math.Mod(-2*L+2*F+2*D+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 1383*sinarg + -2*cosarg + de += -594*cosarg + -2*sinarg + + arg = math.Mod(L+2*F+2*D+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += -1331*sinarg + 8*cosarg + de += 663*cosarg + 4*sinarg + + arg = math.Mod(-2*Lp+2*F+-2*D+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += -1283 * sinarg + de += 672 * cosarg + + arg = math.Mod(-1*L+Lp+D+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 1314 * sinarg + de += -700 * cosarg + + arg = math.Mod(-1*L+2*F+4*D+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += -1521*sinarg + 9*cosarg + de += 647*cosarg + 4*sinarg + + arg = math.Mod(2*F+D+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 1660*sinarg + -5*cosarg + de += -710*cosarg + -2*sinarg + + arg = math.Mod(-1*L+D, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 4026*sinarg + -353*cosarg + de += -553*cosarg + -139*sinarg + + arg = math.Mod(L+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += -1981 * sinarg + de += 854 * cosarg + + arg = math.Mod(-1*L+2*F+-2*D+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += -1987*sinarg + -6*cosarg + de += 1073*cosarg + -2*sinarg + + arg = math.Mod(L+2*F, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 3339*sinarg + -13*cosarg + de += -107*cosarg + 1*sinarg + + arg = math.Mod(L+Lp, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += -3389*sinarg + 5*cosarg + de += 35*cosarg + -2*sinarg + + arg = math.Mod(-1*L+Lp+D, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 3276*sinarg + 1*cosarg + de += -9 * cosarg + + arg = math.Mod(2*L+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 2179*sinarg + -2*cosarg + de += -1129*cosarg + -2*sinarg + + arg = math.Mod(L+Lp+2*F+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 2481*sinarg + -7*cosarg + de += -1062*cosarg + -3*sinarg + + arg = math.Mod(-2*L+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += -2294*sinarg + -10*cosarg + de += 1266*cosarg + -4*sinarg + + arg = math.Mod(-1*Lp+2*F+2*D+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += -2647*sinarg + 11*cosarg + de += 1129*cosarg + 5*sinarg + + arg = math.Mod(-1*L+2*F, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += -4056*sinarg + 5*cosarg + de += 40*cosarg + -2*sinarg + + arg = math.Mod(-1*L+-1*Lp+2*F+2*D+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += -2819*sinarg + 7*cosarg + de += 1207*cosarg + 3*sinarg + + arg = math.Mod(D, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += -4230*sinarg + 5*cosarg + de += -20*cosarg + -2*sinarg + + arg = math.Mod(L+-1*Lp+2*F+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += -2878*sinarg + 8*cosarg + de += 1232*cosarg + 4*sinarg + + arg = math.Mod(-1*Lp+2*D, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 4348*sinarg + -10*cosarg + de += -81*cosarg + 2*sinarg + + arg = math.Mod(3*L+2*F+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += -2904*sinarg + 15*cosarg + de += 1233*cosarg + 7*sinarg + + arg = math.Mod(-2*L+2*F+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += -3075*sinarg + -2*cosarg + de += 1313*cosarg + -1*sinarg + + arg = math.Mod(L+-1*Lp, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 4725*sinarg + -6*cosarg + de += -41*cosarg + 3*sinarg + + arg = math.Mod(Lp+2*F+-2*D+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 3579*sinarg + 5*cosarg + de += -1900*cosarg + 1*sinarg + + arg = math.Mod(L+2*D, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 6579*sinarg + -24*cosarg + de += -199*cosarg + 2*sinarg + + arg = math.Mod(2*L+-2*D+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 4065*sinarg + 6*cosarg + de += -2206*cosarg + 1*sinarg + + arg = math.Mod(-1*L+-1*Lp+2*D, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 7350*sinarg + -8*cosarg + de += -51*cosarg + 4*sinarg + + arg = math.Mod(-2*D+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-4940+-11*T)*sinarg + -21*cosarg + de += 2720*cosarg + -9*sinarg + + arg = math.Mod(-1*Lp+2*F+-2*D+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-4752+-11*T)*sinarg + -3*cosarg + de += 2719*cosarg + -3*sinarg + + arg = math.Mod(2*L+2*F+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += -5350*sinarg + 21*cosarg + de += 2695*cosarg + 12*sinarg + + arg = math.Mod(-2*L+2*D+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-5774+-11*T)*sinarg + -15*cosarg + de += 3041*cosarg + -5*sinarg + + arg = math.Mod(2*L+2*F+-2*D+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 6443*sinarg + -7*cosarg + de += -2768*cosarg + -4*sinarg + + arg = math.Mod(L+2*F+-2*D+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (5800+10*T)*sinarg + 2*cosarg + de += -3045*cosarg + -1*sinarg + + arg = math.Mod(2*D+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-6302+-11*T)*sinarg + 2*cosarg + de += 3272*cosarg + 4*sinarg + + arg = math.Mod(-1*Lp+2*F+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-7141+21*T)*sinarg + 8*cosarg + de += 3070*cosarg + 4*sinarg + + arg = math.Mod(2*F+2*D+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-6637+-11*T)*sinarg + 25*cosarg + de += 3353*cosarg + 14*sinarg + + arg = math.Mod(Lp+2*F+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (7566+-21*T)*sinarg + -11*cosarg + de += -3250*cosarg + -5*sinarg + + arg = math.Mod(-2*L+2*F, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += -11024*sinarg + -14*cosarg + de += 104*cosarg + 2*sinarg + + arg = math.Mod(L+2*F+2*D+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += -7691*sinarg + 44*cosarg + de += 3268*cosarg + 19*sinarg + + arg = math.Mod(2*Lp, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (16707+-85*T)*sinarg + -10*cosarg + de += (168+-1*T)*cosarg + 10*sinarg + + arg = math.Mod(-1*L+2*F+2*D+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += -10204*sinarg + 25*cosarg + de += 5222*cosarg + 15*sinarg + + arg = math.Mod(-1*Lp+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-12654+11*T)*sinarg + 63*cosarg + de += 6415*cosarg + 26*sinarg + + arg = math.Mod(L+-2*D+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-12873+-10*T)*sinarg + -37*cosarg + de += 6953*cosarg + -14*sinarg + + arg = math.Mod(-2*F+2*D, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 21783*sinarg + 13*cosarg + de += -167*cosarg + 13*sinarg + + arg = math.Mod(2*Lp+2*F+-2*D+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-15794+72*T)*sinarg + -16*cosarg + de += (6850+-42*T)*cosarg + -5*sinarg + + arg = math.Mod(-1*L+2*D+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (15164+10*T)*sinarg + 11*cosarg + de += -8001*cosarg + -1*sinarg + + arg = math.Mod(Lp+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-14053+-25*T)*sinarg + 79*cosarg + de += (8551+-2*T)*cosarg + -45*sinarg + + arg = math.Mod(2*F, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 25887*sinarg + -66*cosarg + de += -550*cosarg + 11*sinarg + + arg = math.Mod(2*L, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 29243*sinarg + -74*cosarg + de += -609*cosarg + 13*sinarg + + arg = math.Mod(-1*L+2*F+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (20441+21*T)*sinarg + 10*cosarg + de += -10758*cosarg + -3*sinarg + + arg = math.Mod(L+2*F+-2*D+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 28593*sinarg + -1*cosarg + de += (-12338+10*T)*cosarg + -3*sinarg + + arg = math.Mod(2*L+2*F+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-31046+-1*T)*sinarg + 131*cosarg + de += (13238+-11*T)*cosarg + 59*sinarg + + arg = math.Mod(-2*L+2*D, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += -47722*sinarg + -18*cosarg + de += 477*cosarg + -25*sinarg + + arg = math.Mod(-2*Lp+2*F+-2*D+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += 32481 * sinarg + de += -13870 * cosarg + + arg = math.Mod(2*F+2*D+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-38571+-1*T)*sinarg + 158*cosarg + de += (16452+-11*T)*cosarg + 68*sinarg + + arg = math.Mod(2*D, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (63384+11*T)*sinarg + -150*cosarg + de += -1220*cosarg + 29*sinarg + + arg = math.Mod(-2*L+2*F+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (45893+50*T)*sinarg + 31*cosarg + de += (-24236+-10*T)*cosarg + 20*sinarg + + arg = math.Mod(L+2*F+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-51613+-42*T)*sinarg + 129*cosarg + de += 26366*cosarg + 78*sinarg + + arg = math.Mod(-1*L+2*F+2*D+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-59641+-11*T)*sinarg + 149*cosarg + de += (25543+-11*T)*cosarg + 66*sinarg + + arg = math.Mod(-1*L+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-57976+-63*T)*sinarg + -189*cosarg + de += 31429*cosarg + -75*sinarg + + arg = math.Mod(L+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (63110+63*T)*sinarg + 27*cosarg + de += -33228*cosarg + -9*sinarg + + arg = math.Mod(-1*L+2*D, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (156994+10*T)*sinarg + -168*cosarg + de += -1235*cosarg + 82*sinarg + + arg = math.Mod(-1*L+2*F+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (123457+11*T)*sinarg + 19*cosarg + de += (-53311+32*T)*cosarg + -4*sinarg + + arg = math.Mod(2*F+-2*D+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (128227+137*T)*sinarg + 181*cosarg + de += (-68982+-9*T)*cosarg + 39*sinarg + + arg = math.Mod(-1*Lp+2*F+-2*D+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (215829+-494*T)*sinarg + 111*cosarg + de += (-95929+299*T)*cosarg + 132*sinarg + + arg = math.Mod(L+2*F+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-301461+-36*T)*sinarg + 816*cosarg + de += (129025+-63*T)*cosarg + 367*sinarg + + arg = math.Mod(2*F+Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-387298+-367*T)*sinarg + 380*cosarg + de += (200728+18*T)*cosarg + 318*sinarg + + arg = math.Mod(L, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (711159+73*T)*sinarg + -872*cosarg + de += -6750*cosarg + 358*sinarg + + arg = math.Mod(Lp+2*F+-2*D+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-516821+1226*T)*sinarg + -524*cosarg + de += (224386+-677*T)*cosarg + -174*sinarg + + arg = math.Mod(Lp, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (1475877+-3633*T)*sinarg + 11817*cosarg + de += (73871+-184*T)*cosarg + -1924*sinarg + + arg = math.Mod(2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (2074554+207*T)*sinarg + -698*cosarg + de += (-897492+470*T)*cosarg + -291*sinarg + + arg = math.Mod(2*F+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-2276413+-234*T)*sinarg + 2796*cosarg + de += (978459+-485*T)*cosarg + 1374*sinarg + + arg = math.Mod(2*F+-2*D+2*Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-13170906+-1675*T)*sinarg + -13696*cosarg + de += (5730336+-3015*T)*cosarg + -4587*sinarg + + arg = math.Mod(Om, twopi) + sinarg = math.Sin(arg) + cosarg = math.Cos(arg) + dp += (-172064161+-174666*T)*sinarg + 33386*cosarg + de += (92052331+9086*T)*cosarg + 15377*sinarg + + de *= as2r / 1e7 + dp *= as2r / 1e7 + + // 行星章动修正 + dp += -0.135 * (as2r / 1e3) + de += 0.388 * (as2r / 1e3) + + return dp * deg, de * deg +} diff --git a/basic/precess.go b/basic/precess.go new file mode 100644 index 0000000..ccbae56 --- /dev/null +++ b/basic/precess.go @@ -0,0 +1,309 @@ +package basic + +import "math" + +// Precess 函数实现了基于 Explanatory Supplement 2nd ed. table 3.211.1 的岁差计算 +// ra0, dec0: 初始赤经和赤纬（弧度） +// JD0: 初始儒略日 +// JD: 目标儒略日 +// 返回: 目标历元的赤经和赤纬（弧度） +func PrecessIAU1976(ra0, dec0, JD0, JD float64) (float64, float64) { + // 定义常量 + ra0 *= rad + dec0 *= rad + t := (JD - JD0) / 36525.0 + T := (JD0 - 2451545.0) / 36525.0 + + // Explanatory Supplement 2nd ed. table 3.211.1 + zeta := ((2306.2181+1.39656*T-0.000139*T*T)*t + + (0.30188-0.000344*T)*t*t + 0.017998*t*t*t) / 3600.0 * rad + + z := ((2306.2181+1.39656*T-0.000139*T*T)*t + + (1.09468+0.000066*T)*t*t + 0.018203*t*t*t) / 3600.0 * rad + + theta := ((2004.3109-0.85330*T-0.000217*T*T)*t - + (0.42665+0.000217*T)*t*t - 0.041833*t*t*t) / 3600.0 * rad + + A := math.Cos(dec0) * math.Sin(ra0+z) + B := math.Cos(theta)*math.Cos(dec0)*math.Cos(ra0+z) - math.Sin(theta)*math.Sin(dec0) + C := math.Sin(theta)*math.Cos(dec0)*math.Cos(ra0+z) + math.Cos(theta)*math.Sin(dec0) + + raMinZ := math.Atan2(A, B) + dec := math.Asin(C) + ra := raMinZ + zeta + + return ra / rad, dec / rad +} + +// 常量定义 +const ( + AS2R = 4.848136811095359935899141e-6 // 角秒到弧度的转换因子 + D2PI = 6.283185307179586476925287 // 2π + EPS0 = 84381.406 * AS2R // J2000.0历元的黄赤交角(弧度) +) + +// 赤道坐标结构体 +type Equatorial struct { + RA float64 // 赤经(度) + Dec float64 // 赤纬(度) +} + +// 三维向量 +type Vector3 [3]float64 + +// 三维矩阵 +type Matrix3 [3][3]float64 + +// 计算长期岁差的黄极向量 +func ltpPECL(epj float64) Vector3 { + // 多项式系数 + pqPol := [2][4]float64{ + {+5851.607687, -0.1189000, -0.00028913, +0.000000101}, + {-1600.886300, +1.1689818, -0.00000020, -0.000000437}, + } + + // 周期项系数 + pqPer := [8][5]float64{ + {708.15, -5486.751211, -684.661560, 667.666730, -5523.863691}, + {2309.00, -17.127623, 2446.283880, -2354.886252, -549.747450}, + {1620.00, -617.517403, 399.671049, -428.152441, -310.998056}, + {492.20, 413.442940, -356.652376, 376.202861, 421.535876}, + {1183.00, 78.614193, -186.387003, 184.778874, -36.776172}, + {622.00, -180.732815, -316.800070, 335.321713, -145.278396}, + {882.00, -87.676083, 198.296071, -185.138669, -34.744450}, + {547.00, 46.140315, 101.135679, -120.972830, 22.885731}, + } + + // 计算从J2000.0起算的儒略世纪数 + t := (epj - 2000.0) / 100.0 + + // 初始化P和Q累加器 + p, q := 0.0, 0.0 + + // 周期项 + for i := 0; i < 8; i++ { + w := D2PI * t + a := w / pqPer[i][0] + s := math.Sin(a) + c := math.Cos(a) + p += c*pqPer[i][1] + s*pqPer[i][3] + q += c*pqPer[i][2] + s*pqPer[i][4] + } + + // 多项式项 + w := 1.0 + for i := 0; i < 4; i++ { + p += pqPol[0][i] * w + q += pqPol[1][i] * w + w *= t + } + + // 转换为弧度 + p *= AS2R + q *= AS2R + + // 形成黄极向量 + z := math.Sqrt(math.Max(1.0-p*p-q*q, 0.0)) + s := math.Sin(EPS0) + c := math.Cos(EPS0) + + return Vector3{p, -q*c - z*s, -q*s + z*c} +} + +// 计算长期岁差的赤道极向量 +func ltpPEQU(epj float64) Vector3 { + // 多项式系数 + xyPol := [2][4]float64{ + {+5453.282155, +0.4252841, -0.00037173, -0.000000152}, + {-73750.930350, -0.7675452, -0.00018725, +0.000000231}, + } + + // 周期项系数 + xyPer := [14][5]float64{ + {256.75, -819.940624, 75004.344875, 81491.287984, 1558.515853}, + {708.15, -8444.676815, 624.033993, 787.163481, 7774.939698}, + {274.20, 2600.009459, 1251.136893, 1251.296102, -2219.534038}, + {241.45, 2755.175630, -1102.212834, -1257.950837, -2523.969396}, + {2309.00, -167.659835, -2660.664980, -2966.799730, 247.850422}, + {492.20, 871.855056, 699.291817, 639.744522, -846.485643}, + {396.10, 44.769698, 153.167220, 131.600209, -1393.124055}, + {288.90, -512.313065, -950.865637, -445.040117, 368.526116}, + {231.10, -819.415595, 499.754645, 584.522874, 749.045012}, + {1610.00, -538.071099, -145.188210, -89.756563, 444.704518}, + {620.00, -189.793622, 558.116553, 524.429630, 235.934465}, + {157.87, -402.922932, -23.923029, -13.549067, 374.049623}, + {220.30, 179.516345, -165.405086, -210.157124, -171.330180}, + {1200.00, -9.814756, 9.344131, -44.919798, -22.899655}, + } + + // 计算从J2000.0起算的儒略世纪数 + t := (epj - 2000.0) / 100.0 + + // 初始化X和Y累加器 + x, y := 0.0, 0.0 + + // 周期项 + for i := 0; i < 14; i++ { + w := D2PI * t + a := w / xyPer[i][0] + s := math.Sin(a) + c := math.Cos(a) + x += c*xyPer[i][1] + s*xyPer[i][3] + y += c*xyPer[i][2] + s*xyPer[i][4] + } + + // 多项式项 + w := 1.0 + for i := 0; i < 4; i++ { + x += xyPol[0][i] * w + y += xyPol[1][i] * w + w *= t + } + + // 转换为弧度 + x *= AS2R + y *= AS2R + + // 形成赤道极向量 + w = x*x + y*y + z := 0.0 + if w < 1.0 { + z = math.Sqrt(1.0 - w) + } + + return Vector3{x, y, z} +} + +// 向量叉乘 +func pxp(a, b Vector3) Vector3 { + return Vector3{ + a[1]*b[2] - a[2]*b[1], + a[2]*b[0] - a[0]*b[2], + a[0]*b[1] - a[1]*b[0], + } +} + +// 向量归一化 +func pn(p Vector3) (Vector3, float64) { + w := math.Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]) + if w == 0.0 { + return Vector3{0, 0, 0}, 0 + } + return Vector3{p[0] / w, p[1] / w, p[2] / w}, w +} + +// 计算长期岁差矩阵 +func ltpPMAT(epj float64) Matrix3 { + // 计算赤道极向量和黄极向量 + peqr := ltpPEQU(epj) + pecl := ltpPECL(epj) + + // 计算春分点向量 + v := pxp(peqr, pecl) + eqx, _ := pn(v) + + // 计算中间行向量 + v = pxp(peqr, eqx) + + // 形成旋转矩阵 + return Matrix3{ + {eqx[0], eqx[1], eqx[2]}, + {v[0], v[1], v[2]}, + {peqr[0], peqr[1], peqr[2]}, + } +} + +// 计算包含GCRS框架偏差的长期岁差矩阵 +func ltpPBMAT(epj float64) Matrix3 { + // 框架偏差 (IERS Conventions 2010, Eqs. 5.21 and 5.33) + dx := -0.016617 * AS2R + de := -0.0068192 * AS2R + dr := -0.0146 * AS2R + + // 计算基本岁差矩阵 + rp := ltpPMAT(epj) + + // 应用偏差 + return Matrix3{ + { + rp[0][0] - rp[0][1]*dr + rp[0][2]*dx, + rp[0][0]*dr + rp[0][1] + rp[0][2]*de, + -rp[0][0]*dx - rp[0][1]*de + rp[0][2], + }, + { + rp[1][0] - rp[1][1]*dr + rp[1][2]*dx, + rp[1][0]*dr + rp[1][1] + rp[1][2]*de, + -rp[1][0]*dx - rp[1][1]*de + rp[1][2], + }, + { + rp[2][0] - rp[2][1]*dr + rp[2][2]*dx, + rp[2][0]*dr + rp[2][1] + rp[2][2]*de, + -rp[2][0]*dx - rp[2][1]*de + rp[2][2], + }, + } +} + +// 将赤道坐标转换为直角坐标向量 +func raDecToVector(ra, dec float64) Vector3 { + raRad := ra * math.Pi / 180 + decRad := dec * math.Pi / 180 + cosDec := math.Cos(decRad) + return Vector3{ + math.Cos(raRad) * cosDec, + math.Sin(raRad) * cosDec, + math.Sin(decRad), + } +} + +// 将直角坐标向量转换为赤道坐标 +func vectorToRaDec(v Vector3) (float64, float64) { + dec := math.Asin(v[2]) * 180 / math.Pi + + ra := math.Atan2(v[1], v[0]) * 180 / math.Pi + if ra < 0 { + ra += 360 + } + + return ra, dec +} + +// 应用长期岁差转换 +func Precess(ra, dec, jdFrom, jdTo float64) (float64, float64) { + // 将儒略日转换为儒略纪元 + epjFrom := 2000.0 + (jdFrom-2451545.0)/365.25 + epjTo := 2000.0 + (jdTo-2451545.0)/365.25 + + // 计算从起始历元到J2000.0的逆矩阵 + rpFrom := ltpPMAT(epjFrom) + rpFromInv := Matrix3{ + {rpFrom[0][0], rpFrom[1][0], rpFrom[2][0]}, + {rpFrom[0][1], rpFrom[1][1], rpFrom[2][1]}, + {rpFrom[0][2], rpFrom[1][2], rpFrom[2][2]}, + } + + // 计算从J2000.0到目标历元的矩阵 + rpTo := ltpPMAT(epjTo) + + // 计算复合旋转矩阵 + var rpFinal Matrix3 + for i := 0; i < 3; i++ { + for j := 0; j < 3; j++ { + for k := 0; k < 3; k++ { + rpFinal[i][j] += rpTo[i][k] * rpFromInv[k][j] + } + } + } + + // 将赤道坐标转换为直角坐标向量 + vector := raDecToVector(ra, dec) + + // 应用旋转矩阵 + rotatedVector := Vector3{ + rpFinal[0][0]*vector[0] + rpFinal[0][1]*vector[1] + rpFinal[0][2]*vector[2], + rpFinal[1][0]*vector[0] + rpFinal[1][1]*vector[1] + rpFinal[1][2]*vector[2], + rpFinal[2][0]*vector[0] + rpFinal[2][1]*vector[1] + rpFinal[2][2]*vector[2], + } + + // 将直角坐标向量转换回赤道坐标 + return vectorToRaDec(rotatedVector) +} diff --git a/basic/saturn.go b/basic/saturn.go index c804321..70f0c6b 100644 --- a/basic/saturn.go +++ b/basic/saturn.go @@ -103,7 +103,7 @@ func SaturnApparentLo(JD float64) float64 { lo = Limit360(lo) //lo-=GXCLo(lo,bo,JD)/3600; //bo+=GXCBo(lo,bo,JD); - lo += HJZD(JD) + lo += Nutation2000Bi(JD) return lo } @@ -117,7 +117,7 @@ func SaturnApparentBo(JD float64) float64 { bo = bo * 180 / math.Pi //lo+=GXCLo(lo,bo,JD); //bo+=GXCBo(lo,bo,JD)/3600; - //lo+=HJZD(JD); + //lo+=Nutation2000Bi(JD); return bo } @@ -132,7 +132,7 @@ func SaturnApparentLoBo(JD float64) (float64, float64) { lo = Limit360(lo) //lo-=GXCLo(lo,bo,JD)/3600; //bo+=GXCBo(lo,bo,JD); - lo += HJZD(JD) + lo += Nutation2000Bi(JD) return lo, bo } diff --git a/basic/star.go b/basic/star.go index eddebd6..6113056 100644 --- a/basic/star.go +++ b/basic/star.go @@ -21,7 +21,7 @@ func StarHeight(jde, ra, dec, lon, lat, timezone float64) float64 { return ArcSin(sinHeight) } -//StarAzimuth 星体的方位角 +// StarAzimuth 星体的方位角 // 传入 jde时间、瞬时赤经、瞬时赤纬、经度、纬度、时区，jde时间应为时区时间 // 返回方位角，单位为度，正北为0，度数顺时针增加，取值范围[0-360) func StarAzimuth(jde, ra, dec, lon, lat, timezone float64) float64 { @@ -46,7 +46,7 @@ func StarAzimuth(jde, ra, dec, lon, lat, timezone float64) float64 { return Azimuth } -//StarHourAngle 星体的时角 +// StarHourAngle 星体的时角 // 传入 jde时间、瞬时赤经、瞬时赤纬、经度、时区，jde时间应为时区时间 // 返回时角 func StarHourAngle(jde, ra, lon, timezone float64) float64 { @@ -58,17 +58,69 @@ func StarHourAngle(jde, ra, lon, timezone float64) float64 { return Limit360(st - ra) } -// MeanSiderealTime 不含章动下的恒星时 +// MeanSiderealTime 平恒星时 func MeanSiderealTime(JD float64) float64 { - T := (JD - 2451545) / 36525 - return (Limit360(280.46061837+360.98564736629*(JD-2451545.0)+0.000387933*T*T-T*T*T/38710000) / 15) + return MeanSiderealTime2006(JD) } // ApparentSiderealTime 视恒星时，计算章动 func ApparentSiderealTime(JD float64) float64 { - tmp := MeanSiderealTime(JD) - return tmp + HJZD(JD)*Cos(Sita(JD))/15 + return ApparentSiderealTime2006(JD) } + +// MeanSiderealTime1982 不含章动下的恒星时 +func MeanSiderealTime1982(JD float64) float64 { + T := (JD - 2451545) / 36525 + return (Limit360(280.46061837+360.98564736629*(JD-2451545.0)+0.000387933*T*T-T*T*T/38710000) / 15) +} + +// ApparentSiderealTime1982 视恒星时，计算章动 +func ApparentSiderealTime1982(JD float64) float64 { + tmp := MeanSiderealTime1982(JD) + return tmp + Nutation2000Bi(JD)*Cos(Sita(JD))/15 +} + +// EarthRotationAngle 计算地球自转角 (ERA) +// jd_ut1: UT1 时间的儒略日 +// 返回值: 地球自转角 (弧度) +func EarthRotationAngle(jd_ut1 float64) float64 { + t := jd_ut1 - 2451545.0 + frac := math.Mod(jd_ut1, 1.0) + + era := math.Mod(math.Pi*2*(0.7790572732640+0.00273781191135448*t+frac), math.Pi*2) + if era < 0 { + era += math.Pi * 2 + } + + return era +} + +// MeanSiderealTime2006 计算格林尼治平恒星时 (GMST) +// jd_ut1: UT1 时间的儒略日 +// jd_tt: TT 时间的儒略日 +// 返回值: 格林尼治平恒星时 (弧度) +func MeanSiderealTime2006(jd_ut1 float64) float64 { + jd_tt := TD2UT(jd_ut1, true) + t := (jd_tt - 2451545.0) / 36525.0 + era := EarthRotationAngle(jd_ut1) + + // 公式 2.12 + gmst := math.Mod(era+(0.014506+4612.15739966*t+1.39667721*t*t+ + -0.00009344*t*t*t+0.00001882*t*t*t*t)/60/60*math.Pi/180, math.Pi*2) + + if gmst < 0 { + gmst += math.Pi * 2 + } + + return gmst * deg / 15 +} + +// ApparentSiderealTime2006 视恒星时，计算章动 +func ApparentSiderealTime2006(JD float64) float64 { + tmp := MeanSiderealTime2006(JD) + return tmp + Nutation2000Bi(JD)*Cos(Sita(JD))/15 +} + func StarAngle(RA, DEC, JD, Lon, Lat, TZ float64) float64 { //JD=JD-8/24+TZ/24; calcjd := JD - TZ/24 diff --git a/basic/stardat.go b/basic/stardat.go index 9c201df..151f232 100644 --- a/basic/stardat.go +++ b/basic/stardat.go @@ -9344,7 +9344,7 @@ func StarDataByHR(hr int) (StarData, error) { func (s InnerStarData) RaDecByJde(jde float64) (float64, float64) { //计算自行 year := ((jde - 2451545.0) / 365.2422) - return ZuoBiaoSuiCha(s.Ra+(year*s.PmRA/3600), s.Dec+(year*s.PmDec/3600), 2451545.0, jde) + return Precess(s.Ra+(year*s.PmRA/3600), s.Dec+(year*s.PmDec/3600), 2451545.0, jde) } func (s StarData) RaDecByDate(date time.Time) (float64, float64) { diff --git a/basic/sun.go b/basic/sun.go index 179342e..ebd4158 100644 --- a/basic/sun.go +++ b/basic/sun.go @@ -7,68 +7,7 @@ import ( . "b612.me/astro/tools" ) -/* -黄赤交角、nutation==true时，计算交角章动 -*/ -func EclipticObliquity(jde float64, nutation bool) float64 { - U := (jde - 2451545) / 3652500.000 - sita := 23.000 + 26.000/60.000 + 21.448/3600.000 - ((4680.93*U - 1.55*U*U + 1999.25*U*U*U - 51.38*U*U*U*U - 249.67*U*U*U*U*U - 39.05*U*U*U*U*U*U + 7.12*U*U*U*U*U*U*U + 27.87*U*U*U*U*U*U*U*U + 5.79*U*U*U*U*U*U*U*U*U + 2.45*U*U*U*U*U*U*U*U*U*U) / 3600) - if nutation { - return sita + JJZD(jde) - } else { - return sita - } -} - -func Sita(JD float64) float64 { - return EclipticObliquity(JD, true) -} - -/* - * @name 黄经章动 - */ -func HJZD(jd float64) float64 { // '黄经章动 - t := (jd - 2451545) / 36525.000 - d := 297.8502042 + 445267.1115168*t - 0.0016300*t*t + t*t*t/545868 - t*t*t*t/113065000 - m := SunM(jd) - n := MoonM(jd) - f := MoonLonX(jd) - o := 125.04452 - 1934.136261*t + 0.0020708*t*t + t*t*t/450000 - tp := [][]float64{{0, 0, 0, 0, 1, -171996, -174.2 * t}, {-2, 0, 0, 2, 2, -13187, -1.6 * t}, {0, 0, 0, 2, 2, -2274, -0.2 * t}, {0, 0, 0, 0, 2, 2062, 0.2 * t}, {0, 1, 0, 0, 0, 1426, -3.4 * t}, {0, 0, 1, 0, 0, 712, 0.1 * t}, {-2, 1, 0, 2, 2, -517, 1.2 * t}, {0, 0, 0, 2, 1, -386, -0.4 * t}, {0, 0, 1, 2, 2, -301, 0}, {-2, -1, 0, 2, 2, 217, -0.5 * t}, {-2, 0, 1, 0, 0, -158, 0}, {-2, 0, 0, 2, 1, 129, 0.1 * t}, {0, 0, -1, 2, 2, 123, 0}, {2, 0, 0, 0, 0, 63, 0}, {0, 0, 1, 0, 1, 63, 0.1 * t}, {2, 0, -1, 2, 2, -59, 0}, {0, 0, -1, 0, 1, -58, -0.1 * t}, {0, 0, 1, 2, 1, -51, 0}, {-2, 0, 2, 0, 0, 48, 0}, {0, 0, -2, 2, 1, 46, 0}, {2, 0, 0, 2, 2, -38, 0}, {0, 0, 2, 2, 2, -31, 0}, {0, 0, 2, 0, 0, 29, 0}, {-2, 0, 1, 2, 2, 29, 0}, {0, 0, 0, 2, 0, 26, 0}, {-2, 0, 0, 2, 0, -22, 0}, {0, 0, -1, 2, 1, 21, 0}, {0, 2, 0, 0, 0, 17, -0.1 * t}, {2, 0, -1, 0, 1, 16, 0}, {-2, 2, 0, 2, 2, -16, 0.1 * t}, {0, 1, 0, 0, 1, -15, 0}, {-2, 0, 1, 0, 1, -13, 0}, {0, -1, 0, 0, 1, -12, 0}, {0, 0, 2, -2, 0, 11, 0}, {2, 0, -1, 2, 1, -10, 0}, {2, 0, 1, 2, 2, -8, 0}, {0, 1, 0, 2, 2, 7, 0}, {-2, 1, 1, 0, 0, -7, 0}, {0, -1, 0, 2, 2, -7, 0}, {2, 0, 0, 2, 1, -7, 0}, {2, 0, 1, 0, 0, 6, 0}, {-2, 0, 2, 2, 2, 6, 0}, {-2, 0, 1, 2, 1, 6, 0}, {2, 0, -2, 0, 1, -6, 0}, {2, 0, 0, 0, 1, -6, 0}, {0, -1, 1, 0, 0, 5, 0}, {-2, -1, 0, 2, 1, -5, 0}, {-2, 0, 0, 0, 1, -5, 0}, {0, 0, 2, 2, 1, -5, 0}, {-2, 0, 2, 0, 1, 4, 0}, {-2, 1, 0, 2, 1, 4, 0}, {0, 0, 1, -2, 0, 4, 0}, {-1, 0, 1, 0, 0, -4, 0}, {-2, 1, 0, 0, 0, -4, 0}, {1, 0, 0, 0, 0, -4, 0}, {0, 0, 1, 2, 0, 3, 0}, {0, 0, -2, 2, 2, -3, 0}, {-1, -1, 1, 0, 0, -3, 0}, {0, 1, 1, 0, 0, -3, 0}, {0, -1, 1, 2, 2, -3, 0}, {2, -1, -1, 2, 2, -3, 0}, {0, 0, 3, 2, 2, -3, 0}, {2, -1, 0, 2, 2, -3, 0}} - var s float64 - for i := 0; i < len(tp); i++ { - s += (tp[i][5] + tp[i][6]) * Sin(d*tp[i][0]+m*tp[i][1]+n*tp[i][2]+f*tp[i][3]+o*tp[i][4]) - } - //P=-17.20*Sin(o)-1.32*Sin(2*280.4665 + 36000.7698*t)-0.23*Sin(2*218.3165 + 481267.8813*t )+0.21*Sin(2*o); - //return P/3600; - return (s / 10000) / 3600 -} - -/* - * 交角章动 - */ -func JJZD(jd float64) float64 { //交角章动 - t := (jd - 2451545) / 36525 - //d = 297.85036 +455267.111480*t - 0.0019142*t*t+ t*t*t/189474; - //m = 357.52772 + 35999.050340*t - 0.0001603*t*t- t*t*t/300000; - //n= 134.96298 + 477198.867398*t + 0.0086972*t*t + t*t*t/56250; - //f = 93.27191 + 483202.017538*t - 0.0036825*t*t + t*t*t/327270; - d := 297.8502042 + 445267.1115168*t - 0.0016300*t*t + t*t*t/545868 - t*t*t*t/113065000 - m := SunM(jd) - n := MoonM(jd) - f := MoonLonX(jd) - o := 125.04452 - 1934.136261*t + 0.0020708*t*t + t*t*t/450000 - tp := [][]float64{{0, 0, 0, 0, 1, 92025, 8.9 * t}, {-2, 0, 0, 2, 2, 5736, -3.1 * t}, {0, 0, 0, 2, 2, 977, -0.5 * t}, {0, 0, 0, 0, 2, -895, 0.5 * t}, {0, 1, 0, 0, 0, 54, -0.1 * t}, {0, 0, 1, 0, 0, -7, 0}, {-2, 1, 0, 2, 2, 224, -0.6 * t}, {0, 0, 0, 2, 1, 200, 0}, {0, 0, 1, 2, 2, 129, -0.1 * t}, {-2, -1, 0, 2, 2, -95, 0.3 * t}, {-2, 0, 0, 2, 1, -70, 0}, {0, 0, -1, 2, 2, -53, 0}, {2, 0, 0, 0, 0, 63, 0}, {0, 0, 1, 0, 1, -33, 0}, {2, 0, -1, 2, 2, 26, 0}, {0, 0, -1, 0, 1, 32, 0}, {0, 0, 1, 2, 1, 27, 0}, {0, 0, -2, 2, 1, -24, 0}, {2, 0, 0, 2, 2, 16, 0}, {0, 0, 2, 2, 2, 13, 0}, {-2, 0, 1, 2, 2, -12, 0}, {0, 0, -1, 2, 1, -10, 0}, {2, 0, -1, 0, 1, -8, 0}, {-2, 2, 0, 2, 2, 7, 0}, {0, 1, 0, 0, 1, 9, 0}, {-2, 0, 1, 0, 1, 7, 0}, {0, -1, 0, 0, 1, 6, 0}, {2, 0, -1, 2, 1, 5, 0}, {2, 0, 1, 2, 2, 3, 0}, {0, 1, 0, 2, 2, -3, 0}, {0, -1, 0, 2, 2, 3, 0}, {2, 0, 0, 2, 1, 3, 0}, {-2, 0, 2, 2, 2, -3, 0}, {-2, 0, 1, 2, 1, -3, 0}, {2, 0, -2, 0, 1, 3, 0}, {2, 0, 0, 0, 1, 3, 0}, {-2, -1, 0, 2, 1, 3, 0}, {-2, 0, 0, 0, 1, 3, 0}, {0, 0, 2, 2, 1, 3, 0}} - var s float64 = 0 - for i := 0; i < len(tp); i++ { - s += (tp[i][5] + tp[i][6]) * Cos(d*tp[i][0]+m*tp[i][1]+n*tp[i][2]+f*tp[i][3]+o*tp[i][4]) - } - return s / 10000 / 3600 -} - -/* -@name 太阳几何黄经 -*/ +// SunLo 太阳几何黄经 func SunLo(jd float64) float64 { T := (jd - 2451545) / 365250 SunLo := 280.4664567 + 360007.6982779*T + 0.03032028*T*T + T*T*T/49931 - T*T*T*T/15299 - T*T*T*T*T/1988000 @@ -94,12 +33,14 @@ func EarthPI(JD float64) float64 { //近日點經度 T := (JD - 2451545) / 36525 return 102.93735 + 1.71953*T + 000046*T*T } + func SunMidFun(JD float64) float64 { //'太阳中间方程 T := (JD - 2451545) / 36525 M := SunM(JD) SunMidFun := (1.9146-0.004817*T-0.000014*T*T)*Sin(M) + (0.019993-0.000101*T)*Sin(2*M) + 0.00029*Sin(3*M) return SunMidFun } + func SunTrueLo(JD float64) float64 { // '太阳真黄经 SunTrueLo := SunLo(JD) + SunMidFun(JD) return SunTrueLo @@ -121,19 +62,18 @@ func SunApparentRaDec(JD float64) (float64, float64) { } func SunTrueRa(JD float64) float64 { //'太阳真赤经 - sitas := Sita(JD) - SunTrueRa := ArcTan(Cos(sitas) * Sin(SunTrueLo(JD)) / Cos(SunTrueLo(JD))) + sunTrueRa := ArcTan(Cos(sitas) * Sin(SunTrueLo(JD)) / Cos(SunTrueLo(JD))) //Select Case SunTrueLo(JD) tmp := SunTrueLo(JD) if tmp >= 90 && tmp < 180 { - SunTrueRa = 180 + SunTrueRa + sunTrueRa = 180 + sunTrueRa } else if tmp >= 180 && tmp < 270 { - SunTrueRa = 180 + SunTrueRa + sunTrueRa = 180 + sunTrueRa } else if tmp >= 270 && tmp <= 360 { - SunTrueRa = 360 + SunTrueRa + sunTrueRa = 360 + sunTrueRa } - return SunTrueRa + return sunTrueRa } func SunApparentDec(JD float64) float64 { // '太阳视赤纬 @@ -148,9 +88,9 @@ func SunTrueDec(JD float64) float64 { // '太阳真赤纬 SunTrueDec := ArcSin(Sin(sitas) * Sin(SunTrueLo(JD))) return SunTrueDec } -func SunTime(JD float64) float64 { //均时差 - tm := (SunLo(JD) - 0.0057183 - (HSunApparentRa(JD)) + (HJZD(JD))*Cos(Sita(JD))) / 15 +func SunTime(JD float64) float64 { //均时差 + tm := (SunLo(JD) - 0.0057183 - (HSunApparentRa(JD)) + (Nutation2000Bi(JD))*Cos(Sita(JD))) / 15 if tm > 23 { tm = -24 + tm } @@ -158,7 +98,6 @@ func SunTime(JD float64) float64 { //均时差 } func SunSC(Lo, JD float64) float64 { //黄道上的岁差，仅黄纬=0时 - t := (JD - 2451545) / 36525 //n := 47.0029/3600*t - 0.03302/3600*t*t + 0.000060/3600*t*t*t //m := 174.876384/3600 - 869.8089/3600*t + 0.03536/3600*t*t @@ -166,6 +105,7 @@ func SunSC(Lo, JD float64) float64 { //黄道上的岁差，仅黄纬=0时 return Lo + pk } +// 高精度，使用VSOP87 func HSunTrueLo(JD float64) float64 { L := planet.WherePlanet(0, 0, JD) return L @@ -178,15 +118,7 @@ func HSunTrueBo(JD float64) float64 { func HSunApparentLo(JD float64) float64 { L := HSunTrueLo(JD) - /* - t := (JD - 2451545) / 365250.0 - R := planet.WherePlanet(-1, 2, JD) - t2 := t * t - t3 := t2 * t //千年数的各次方 - R += (-0.0020 + 0.0044*t + 0.0213*t2 - 0.0250*t3) - L = L + HJZD(JD) - 20.4898/R/3600.00 - */ - L = L + HJZD(JD) + SunLoGXC(JD) + L = L + Nutation2000Bi(JD) + SunLoGXC(JD) return L } @@ -209,35 +141,25 @@ func HSunApparentRa(JD float64) float64 { // '太阳视赤经 return LoToRa(JD, HSunApparentLo(JD), HSunTrueBo(JD)) } -func HSunTrueRa(JD float64) float64 { //'太阳真赤经 +func HSunTrueRa(JD float64) float64 { tmp := HSunTrueLo(JD) sitas := Sita(JD) - HSunTrueRa := ArcTan(Cos(sitas) * Sin(tmp) / Cos(tmp)) - //Select Case SunTrueLo(JD) - if tmp >= 90 && tmp < 180 { - HSunTrueRa = 180 + HSunTrueRa - } else if tmp >= 180 && tmp < 270 { - HSunTrueRa = 180 + HSunTrueRa - } else if tmp >= 270 && tmp <= 360 { - HSunTrueRa = 360 + HSunTrueRa - } - return HSunTrueRa + + numerator := Cos(sitas) * Sin(tmp) + denominator := Cos(tmp) + + return ArcTan2(numerator, denominator) } func HSunApparentDec(JD float64) float64 { // '太阳视赤纬 - T := (JD - 2451545) / 36525 - sitas := EclipticObliquity(JD, false) + 0.00256*Cos(125.04-1934.136*T) - HSunApparentDec := ArcSin(Sin(sitas) * Sin(HSunApparentLo(JD))) - return HSunApparentDec + return ArcSin(Sin(EclipticObliquity(JD, true)) * Sin(HSunApparentLo(JD))) } func HSunTrueDec(JD float64) float64 { // '太阳真赤纬 - sitas := EclipticObliquity(JD, false) - HSunTrueDec := ArcSin(Sin(sitas) * Sin(HSunTrueLo(JD))) - return HSunTrueDec + return ArcSin(Sin(EclipticObliquity(JD, false)) * Sin(HSunTrueLo(JD))) } -func RDJL(jd float64) float64 { //ri di ju li +func Distance(jd float64) float64 { //ri di ju li f := SunMidFun(jd) m := SunM(jd) e := Earthe(jd) @@ -270,7 +192,6 @@ func GetMoonLoops(year float64, loop int) []float64 { moon[j] = newMoon tmp = moon[j] i++ - // echo DateCalc(moon[i])."
"; } return moon } @@ -419,10 +340,8 @@ func GetWHTime(Year, Angle int) float64 { return TD2UT(JD1, false) } -/* - * 太阳中天时刻，通过均时差计算 - */ -func GetSunTZTime(JD, Lon, TZ float64) float64 { //实际中天时间 +// 太阳中天时刻，通过均时差计算 +func CulminationTime(JD, Lon, TZ float64) float64 { //实际中天时间 JD = math.Floor(JD) tmp := (TZ*15 - Lon) * 4 / 60 return JD + tmp/24.0 - SunTime(JD)/24.0 @@ -431,23 +350,23 @@ func GetSunTZTime(JD, Lon, TZ float64) float64 { //实际中天时间 /* * 昏朦影传入 当天0时时刻 */ -func GetBanTime(JD, Lon, Lat, TZ, An float64) float64 { +func EveningTwilight(JD, Lon, Lat, TZ, An float64) float64 { JD = math.Floor(JD) + 1.5 ntz := math.Round(Lon / 15) - tztime := GetSunTZTime(JD, Lon, ntz) - if SunHeight(tztime, Lon, Lat, ntz) < An { + culminationTime := CulminationTime(JD, Lon, ntz) + if SunHeight(culminationTime, Lon, Lat, ntz) < An { return -2 //极夜 } - if SunHeight(tztime+0.5, Lon, Lat, ntz) > An { + if SunHeight(culminationTime+0.5, Lon, Lat, ntz) > An { return -1 //极昼 } - tmp := (Sin(An) - Sin(HSunApparentDec(tztime))*Sin(Lat)) / (Cos(HSunApparentDec(tztime)) * Cos(Lat)) + tmp := (Sin(An) - Sin(HSunApparentDec(culminationTime))*Sin(Lat)) / (Cos(HSunApparentDec(culminationTime)) * Cos(Lat)) var sundown float64 if math.Abs(tmp) <= 1 && Lat < 85 { rzsc := ArcCos(tmp) / 15 - sundown = tztime + rzsc/24.0 + 35.0/24.0/60.0 + sundown = culminationTime + rzsc/24.0 + 35.0/24.0/60.0 } else { - sundown = tztime + sundown = culminationTime i := 0 for LowSunHeight(sundown, Lon, Lat, ntz) > An { i++ @@ -470,42 +389,51 @@ func GetBanTime(JD, Lon, Lat, TZ, An float64) float64 { return JD1 - ntz/24 + TZ/24 } -func GetAsaTime(JD, Lon, Lat, TZ, An float64) float64 { +func MorningTwilight(JD, Lon, Lat, TZ, An float64) float64 { + // 调整到中午12点 JD = math.Floor(JD) + 1.5 + + // 计算时区 ntz := math.Round(Lon / 15) - tztime := GetSunTZTime(JD, Lon, ntz) - if SunHeight(tztime, Lon, Lat, ntz) < An { - return -2 //极夜 + + // 计算太阳上中天时间 + culminationTime := CulminationTime(JD, Lon, ntz) + + // 检查极夜和极昼条件 + if SunHeight(culminationTime, Lon, Lat, ntz) < An { + return -2 // 极夜 } - if SunHeight(tztime-0.5, Lon, Lat, ntz) > An { - return -1 //极昼 + if SunHeight(culminationTime-0.5, Lon, Lat, ntz) > An { + return -1 // 极昼 } - tmp := (Sin(An) - Sin(HSunApparentDec(tztime))*Sin(Lat)) / (Cos(HSunApparentDec(tztime)) * Cos(Lat)) + + // 计算日出时间 + sunDec := HSunApparentDec(culminationTime) + tmp := (Sin(An) - Sin(sunDec)*Sin(Lat)) / (Cos(sunDec) * Cos(Lat)) + var sunrise float64 if math.Abs(tmp) <= 1 && Lat < 85 { - rzsc := ArcCos(tmp) / 15 - sunrise = tztime - rzsc/24 - 25.0/24.0/60.0 + hourAngle := ArcCos(tmp) / 15 + sunrise = culminationTime - hourAngle/24 - 25.0/(24.0*60.0) } else { - sunrise = tztime - i := 0 - for LowSunHeight(sunrise, Lon, Lat, ntz) > An { - i++ - sunrise -= 15.0 / 60.0 / 24.0 - if i > 48 { - break - } + sunrise = culminationTime + for i := 0; i < 48 && LowSunHeight(sunrise, Lon, Lat, ntz) > An; i++ { + sunrise -= 15.0 / (60.0 * 24.0) // 每次减少15分钟 } } - JD1 := sunrise - 5.00/24.00/60.00 + + JD1 := sunrise - 5.0/(24.0*60.0) for { JD0 := JD1 - stDegree := SunHeight(JD0, Lon, Lat, ntz) - An - stDegreep := (SunHeight(JD0+0.000005, Lon, Lat, ntz) - SunHeight(JD0-0.000005, Lon, Lat, ntz)) / 0.00001 - JD1 = JD0 - stDegree/stDegreep + heightDiff := SunHeight(JD0, Lon, Lat, ntz) - An + heightDerivative := (SunHeight(JD0+0.000005, Lon, Lat, ntz) - SunHeight(JD0-0.000005, Lon, Lat, ntz)) / 0.00001 + JD1 = JD0 - heightDiff/heightDerivative + if math.Abs(JD1-JD0) < 0.00001 { break } } + return JD1 - ntz/24 + TZ/24 } @@ -544,7 +472,7 @@ func calculateSunRiseSetTime(julianDay, longitude, latitude, timeZone, zenithShi sunAngle = sunAngle - HeightDegreeByLat(height, latitude) // 获取太阳上中天时间 - solarNoonTime := GetSunTZTime(julianDay, longitude, naturalTimeZone) + solarNoonTime := CulminationTime(julianDay, longitude, naturalTimeZone) // 检查极夜极昼条件 polarCondition := checkPolarConditions(solarNoonTime, longitude, latitude, naturalTimeZone, sunAngle, isSunrise) diff --git a/basic/sun_test.go b/basic/sun_test.go index 9b5bac1..1954d56 100644 --- a/basic/sun_test.go +++ b/basic/sun_test.go @@ -22,12 +22,12 @@ func Test_Jq(t *testing.T) { func TestZD(t *testing.T) { jde := 2452982.9872345612 - zd := HJZD(jde) + zd := Nutation2000Bi(jde) fmt.Println(zd) if zd != -0.003746747950462434 { t.Fatal("not equal") } - zd = JJZD(jde) + zd = Nutation1980s(jde) fmt.Println(zd) if zd != 0.001513453926274198 { t.Fatal("not equal") @@ -121,10 +121,10 @@ func Test_SunRise(t *testing.T) { func Test_SunTwilightMo(t *testing.T) { cst := time.FixedZone("cst", 8*3600) jde := Date2JDE(time.Date(2023, 10, 3, 15, 59, 0, 0, cst)) - fmt.Println(GetAsaTime(jde, 113.58880556, 87.66833333, 8, -6)) + fmt.Println(MorningTwilight(jde, 113.58880556, 87.66833333, 8, -6)) for i := 10.0; i < 90.0; i += 0.3 { - fmt.Println(i, GetAsaTime(jde, 125.45506654, float64(i), 8, -6)) + fmt.Println(i, MorningTwilight(jde, 125.45506654, float64(i), 8, -6)) } } @@ -132,7 +132,7 @@ func Test_SunTwilightEv(t *testing.T) { cst := time.FixedZone("cst", 8*3600) jde := Date2JDE(time.Date(2023, 10, 3, 15, 59, 0, 0, cst)) for i := 10.0; i < 90.0; i += 0.3 { - fmt.Println(i, GetBanTime(jde, 115, float64(i), 8, -18)) + fmt.Println(i, EveningTwilight(jde, 115, float64(i), 8, -18)) } } diff --git a/basic/uranus.go b/basic/uranus.go index 8fe23eb..2086cec 100644 --- a/basic/uranus.go +++ b/basic/uranus.go @@ -103,7 +103,7 @@ func UranusApparentLo(JD float64) float64 { lo = Limit360(lo) //lo-=GXCLo(lo,bo,JD)/3600; //bo+=GXCBo(lo,bo,JD); - lo += HJZD(JD) + lo += Nutation2000Bi(JD) return lo } @@ -117,7 +117,7 @@ func UranusApparentBo(JD float64) float64 { bo = bo * 180 / math.Pi //lo+=GXCLo(lo,bo,JD); //bo+=GXCBo(lo,bo,JD)/3600; - //lo+=HJZD(JD); + //lo+=Nutation2000Bi(JD); return bo } @@ -132,7 +132,7 @@ func UranusApparentLoBo(JD float64) (float64, float64) { lo = Limit360(lo) //lo-=GXCLo(lo,bo,JD)/3600; //bo+=GXCBo(lo,bo,JD); - lo += HJZD(JD) + lo += Nutation2000Bi(JD) return lo, bo } diff --git a/basic/venus.go b/basic/venus.go index 71498bc..d5f957f 100644 --- a/basic/venus.go +++ b/basic/venus.go @@ -103,7 +103,7 @@ func VenusApparentLo(JD float64) float64 { lo = Limit360(lo) //lo-=GXCLo(lo,bo,JD)/3600; //bo+=GXCBo(lo,bo,JD); - lo += HJZD(JD) + lo += Nutation2000Bi(JD) return lo } @@ -117,7 +117,7 @@ func VenusApparentBo(JD float64) float64 { bo = bo * 180 / math.Pi //lo+=GXCLo(lo,bo,JD); //bo+=GXCBo(lo,bo,JD)/3600; - //lo+=HJZD(JD); + //lo+=Nutation2000Bi(JD); return bo } @@ -132,7 +132,7 @@ func VenusApparentLoBo(JD float64) (float64, float64) { lo = Limit360(lo) //lo-=GXCLo(lo,bo,JD)/3600; //bo+=GXCBo(lo,bo,JD); - lo += HJZD(JD) + lo += Nutation2000Bi(JD) return lo, bo } diff --git a/moon/moon.go b/moon/moon.go index f7f2f0b..a0fa83e 100644 --- a/moon/moon.go +++ b/moon/moon.go @@ -254,7 +254,7 @@ func PhaseDesc(date time.Time) string { // 返回Date对应UTC世界时的月相大小 func Phase(date time.Time) float64 { jde := basic.Date2JDE(date.UTC()) - return basic.MoonLight(basic.TD2UT(jde, true)) + return basic.MoonPhase(basic.TD2UT(jde, true)) } // ShuoYue 朔月 diff --git a/sun/sun.go b/sun/sun.go index b637b6b..43e46bc 100644 --- a/sun/sun.go +++ b/sun/sun.go @@ -107,7 +107,7 @@ func MorningTwilight(date time.Time, lon, lat, angle float64) (time.Time, error) jde := basic.Date2JDE(date) _, loc := date.Zone() timezone := float64(loc) / 3600.0 - calcJde := basic.GetAsaTime(jde, lon, lat, timezone, angle) + calcJde := basic.MorningTwilight(jde, lon, lat, timezone, angle) if calcJde == -2 { err = ERR_TWILIGHT_NOT_EXISTS } @@ -131,7 +131,7 @@ func EveningTwilight(date time.Time, lon, lat, angle float64) (time.Time, error) _, loc := date.Zone() timezone := float64(loc) / 3600.0 //不需要进行力学时转换，会在GetBanTime中转换 - calcJde := basic.GetBanTime(jde, lon, lat, timezone, angle) + calcJde := basic.EveningTwilight(jde, lon, lat, timezone, angle) if calcJde == -2 { err = ERR_TWILIGHT_NOT_EXISTS } @@ -152,22 +152,40 @@ func EclipticObliquity(date time.Time, nutation bool) float64 { return basic.EclipticObliquity(jde, nutation) } -// EclipticNutation 黄经章动 +// EclipticNutation 黄经章动(2000b) // 返回date对应UTC世界时的黄经章动 func EclipticNutation(date time.Time) float64 { //转换为UTC时间 jde := basic.Date2JDE(date.UTC()) //进行力学时转换与章动计算 - return basic.HJZD(basic.TD2UT(jde, true)) + return basic.Nutation2000Bi(basic.TD2UT(jde, true)) } -// AxialtiltNutation 交角章动 +// EclipticNutation1980 黄经章动(iau 1980) +// 返回date对应UTC世界时的黄经章动 +func EclipticNutation1980(date time.Time) float64 { + //转换为UTC时间 + jde := basic.Date2JDE(date.UTC()) + //进行力学时转换与章动计算 + return basic.Nutation2000Bi(basic.TD2UT(jde, true)) +} + +// AxialtiltNutation 交角章动(2000b) // 返回date对应UTC世界时的交角章动 func AxialtiltNutation(date time.Time) float64 { //转换为UTC时间 jde := basic.Date2JDE(date.UTC()) //进行力学时转换与章动计算 - return basic.JJZD(basic.TD2UT(jde, true)) + return basic.Nutation2000Bs(basic.TD2UT(jde, true)) +} + +// AxialtiltNutation1980 交角章动(1980) +// 返回date对应UTC世界时的交角章动 +func AxialtiltNutation1980(date time.Time) float64 { + //转换为UTC时间 + jde := basic.Date2JDE(date.UTC()) + //进行力学时转换与章动计算 + return basic.Nutation1980s(basic.TD2UT(jde, true)) } // GeometricLo 太阳几何黄经 @@ -274,7 +292,7 @@ func CulminationTime(date time.Time, lon float64) time.Time { jde := basic.Date2JDE(date.Add(time.Duration(-1*date.Hour())*time.Hour)) + 0.5 _, loc := date.Zone() timezone := float64(loc) / 3600.0 - calcJde := basic.GetSunTZTime(jde, lon, timezone) - timezone/24.00 + calcJde := basic.CulminationTime(jde, lon, timezone) - timezone/24.00 return basic.JDE2DateByZone(calcJde, date.Location(), false) } diff --git a/tools/math.go b/tools/math.go index adf64a0..de8c0ac 100644 --- a/tools/math.go +++ b/tools/math.go @@ -28,6 +28,15 @@ func ArcTan(x float64) float64 { return (math.Atan(x) / math.Pi * 180.00000) } +// ArcTan2 计算两个变量的反正切并转换为角度，处理所有象限 +func ArcTan2(y, x float64) float64 { + angle := math.Atan2(y, x) * 180.0 / math.Pi + if angle < 0 { + angle += 360.0 + } + return angle +} + func FloatRound(f float64, n int) float64 { p := math.Pow10(n) return math.Floor(f*p+0.5) / p