From b4b2cea7388e7241b42f0a7a242ea7f4c231b3bd Mon Sep 17 00:00:00 2001 From: starainrt Date: Tue, 4 Jan 2022 14:24:44 +0800 Subject: [PATCH] =?UTF-8?q?=E6=B7=BB=E5=8A=A0=E6=B5=B7=E6=8B=94=E5=9B=A0?= =?UTF-8?q?=E7=B4=A0?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- basic/all_test.go | 20 +++++++++++++----- basic/earth.go | 49 +++++++++++++++++++++++++++++++++++++++++++++ basic/earth_test.go | 13 ++++++++++++ basic/moon.go | 16 ++++++++------- basic/moon_test.go | 13 ++++++++++-- basic/star_test.go | 13 ++++++++++-- basic/sun.go | 18 +++++++++-------- basic/sun_test.go | 39 +++++++++++++++++++++++++++--------- moon/moon.go | 12 +++++++---- moon/moon_test.go | 2 +- sun/sun.go | 10 +++++---- 11 files changed, 162 insertions(+), 43 deletions(-) create mode 100644 basic/earth.go create mode 100644 basic/earth_test.go diff --git a/basic/all_test.go b/basic/all_test.go index 68b01fb..aceb76f 100644 --- a/basic/all_test.go +++ b/basic/all_test.go @@ -5,7 +5,17 @@ import ( "testing" ) -func Test_SunAndMoon(t *testing.T) { +func Test_All(t *testing.T) { + show() +} + +func Benchmark_All(b *testing.B) { + for i := 0; i < b.N; i++ { + show() + } +} + +func show() { jde := GetNowJDE() - 1 ra := HSunSeeRa(jde - 8.0/24.0) dec := HSunSeeDec(jde - 8.0/24.0) @@ -15,15 +25,15 @@ func Test_SunAndMoon(t *testing.T) { fmt.Println("当前太阳赤纬：", dec) fmt.Println("当前太阳星座：", WhichCst(ra, dec, jde)) fmt.Println("当前黄赤交角：", EclipticObliquity(jde-8.0/24.0, true)) - fmt.Println("当前日出：", JDE2Date(GetSunRiseTime(jde, 115, 32, 8, 1))) - fmt.Println("当前日落：", JDE2Date(GetSunDownTime(jde, 115, 32, 8, 1))) + fmt.Println("当前日出：", JDE2Date(GetSunRiseTime(jde, 115, 32, 8, 1, 10))) + fmt.Println("当前日落：", JDE2Date(GetSunDownTime(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.Print("农历：") fmt.Println(GetLunar(2019, 10, 23)) - fmt.Println("当前月出：", JDE2Date(GetMoonRiseTime(jde, 115, 32, 8, 1))) - fmt.Println("当前月落：", JDE2Date(GetMoonDownTime(jde, 115, 32, 8, 1))) + fmt.Println("当前月出：", JDE2Date(GetMoonRiseTime(jde, 115, 32, 8, 1, 10))) + fmt.Println("当前月落：", JDE2Date(GetMoonDownTime(jde, 115, 32, 8, 1, 10))) fmt.Println("月相：", MoonLight(jde-8.0/24.0)) } diff --git a/basic/earth.go b/basic/earth.go new file mode 100644 index 0000000..815945d --- /dev/null +++ b/basic/earth.go @@ -0,0 +1,49 @@ +package basic + +import ( + . "b612.me/astro/tools" + "math" +) + +//地球常数 + +const ( + EARTH_EQUATORIAL_RADIUS float64 = 6378137.0 + EARTH_POLAR_RADIUS float64 = 6356752.3 + EARTH_AVERAGE_RADIUS float64 = 6371393.0 +) + +// HeightDistance 高度与地平线距离的关系（单位：米） +func HeightDistance(height float64) float64 { + return math.Acos((EARTH_AVERAGE_RADIUS)/(EARTH_AVERAGE_RADIUS+height)) * EARTH_AVERAGE_RADIUS +} + +// HeightDistance 高度（单位：米）与地平线下角度的关系（单位：度） +func HeightDegree(height float64) float64 { + return math.Acos((EARTH_AVERAGE_RADIUS)/(EARTH_AVERAGE_RADIUS+height)) * 180 / math.Pi / 2 +} + +// HeightDistanceByLat 不同纬度下高度与地平线距离的关系（单位：米） +func HeightDistanceByLat(height, lat float64) float64 { + raduis := GeocentricRadius(lat) + return math.Acos((raduis)/(raduis+height)) * raduis +} + +// HeightDegreeByLat 不同纬度下高度（单位：米）与地平线下角度的关系（单位：度） +func HeightDegreeByLat(height, lat float64) float64 { + raduis := GeocentricRadius(lat) + return math.Acos((raduis)/(raduis+height)) * 180 / math.Pi / 2 +} + +// GeocentricRadius 地心直径与纬度的关系 +func GeocentricRadius(lat float64) float64 { + a := (EARTH_EQUATORIAL_RADIUS * EARTH_EQUATORIAL_RADIUS * Cos(lat)) + a *= a + b := (EARTH_POLAR_RADIUS * EARTH_POLAR_RADIUS * Sin(lat)) + b *= b + c := (EARTH_EQUATORIAL_RADIUS * Cos(lat)) + c *= c + d := (EARTH_POLAR_RADIUS * Sin(lat)) + d *= d + return math.Sqrt((a + b) / (c + d)) +} diff --git a/basic/earth_test.go b/basic/earth_test.go new file mode 100644 index 0000000..2bbf866 --- /dev/null +++ b/basic/earth_test.go @@ -0,0 +1,13 @@ +package basic + +import ( + "fmt" + "math" + "testing" +) + +func Test_EarthFn(t *testing.T) { + fmt.Println(HeightDistance(10000)) + //近似算法，差距在接受范围内？ + fmt.Println(math.Sqrt(((EARTH_AVERAGE_RADIUS)*2 + 10000) * 10000)) +} diff --git a/basic/moon.go b/basic/moon.go index b02c27c..4841934 100644 --- a/basic/moon.go +++ b/basic/moon.go @@ -1137,7 +1137,7 @@ func CalcMoonSH(Year float64, C int) float64 { stDegree := SunMoonSeek(JD0) - float64(C) stDegreep := (SunMoonSeek(JD0+0.000005) - SunMoonSeek(JD0-0.000005)) / 0.00001 JD1 = JD0 - stDegree/stDegreep - if math.Floor(JD1-JD0) <= 0.000001 { + if math.Abs(JD1-JD0) <= 0.00001 { break } } @@ -1220,7 +1220,7 @@ func CalcMoonXH(Year float64, C int) float64 { stDegree := SunMoonSeek(JD0) - float64(C) stDegreep := (SunMoonSeek(JD0+0.000005) - SunMoonSeek(JD0-0.000005)) / 0.00001 JD1 = JD0 - stDegree/stDegreep - if math.Floor(JD1-JD0) <= 0.000001 { + if math.Abs(JD1-JD0) <= 0.00001 { break } } @@ -1382,7 +1382,7 @@ func GetMoonTZTime(JD, Lon, Lat, TZ float64) float64 { //实际中天时间{ stDegree := MoonTimeAngle(JD0, Lon, Lat, TZ) - 359.599 stDegreep := (MoonTimeAngle(JD0+0.000005, Lon, Lat, TZ) - MoonTimeAngle(JD0-0.000005, Lon, Lat, TZ)) / 0.00001 JD1 = JD0 - stDegree/stDegreep - if math.Floor(JD1-JD0) <= 0.000001 { + if math.Abs(JD1-JD0) <= 0.00001 { break } } @@ -1398,7 +1398,7 @@ func MoonTimeAngle(JD, Lon, Lat, TZ float64) float64 { return timeangle } -func GetMoonRiseTime(JD, Lon, Lat, TZ, ZS float64) float64 { +func GetMoonRiseTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 { ntz := TZ TZ = Lon / 15 var An, tms float64 = 0, 0 @@ -1409,6 +1409,7 @@ func GetMoonRiseTime(JD, Lon, Lat, TZ, ZS float64) float64 { if ZS != 0 { An = -0.83333 //修正大气折射 } + An = An - HeightDegreeByLat(HEI, Lat) moonang := MoonTimeAngle(JD, Lon, Lat, TZ) if moonheight > 0 { //月亮在地平线上或在落下与下中天之间 if moonang > 180 { @@ -1448,7 +1449,7 @@ func GetMoonRiseTime(JD, Lon, Lat, TZ, ZS float64) float64 { stDegree := HMoonHeight(JD0, Lon, Lat, TZ) - An stDegreep := (HMoonHeight(JD0+0.000005, Lon, Lat, TZ) - HMoonHeight(JD0-0.000005, Lon, Lat, TZ)) / 0.00001 JD1 = JD0 - stDegree/stDegreep - if math.Floor(JD1-JD0) <= 0.000001 { + if math.Abs(JD1-JD0) <= 0.00002 { break } } @@ -1461,7 +1462,7 @@ func GetMoonRiseTime(JD, Lon, Lat, TZ, ZS float64) float64 { } } -func GetMoonDownTime(JD, Lon, Lat, TZ, ZS float64) float64 { +func GetMoonDownTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 { ntz := TZ TZ = Lon / 15 var An, tms float64 = 0, 0 @@ -1472,6 +1473,7 @@ func GetMoonDownTime(JD, Lon, Lat, TZ, ZS float64) float64 { if ZS != 0 { An = -0.83333 //修正大气折射 } + An = An - HeightDegreeByLat(HEI, Lat) moonang := MoonTimeAngle(JD, Lon, Lat, TZ) if moonheight < 0 { //月亮在地平线上或在落下与下中天之间 tms = (360 - moonang) / 15 @@ -1508,7 +1510,7 @@ func GetMoonDownTime(JD, Lon, Lat, TZ, ZS float64) float64 { stDegree := HMoonHeight(JD0, Lon, Lat, TZ) - An stDegreep := (HMoonHeight(JD0+0.000005, Lon, Lat, TZ) - HMoonHeight(JD0-0.000005, Lon, Lat, TZ)) / 0.00001 JD1 = JD0 - stDegree/stDegreep - if math.Floor(JD1-JD0) <= 0.000001 { + if math.Abs(JD1-JD0) <= 0.00002 { break } } diff --git a/basic/moon_test.go b/basic/moon_test.go index cef2c61..46c4987 100644 --- a/basic/moon_test.go +++ b/basic/moon_test.go @@ -6,13 +6,22 @@ import ( "time" ) +func Benchmark_MoonRiseBench(b *testing.B) { + jde := GetNowJDE() + for i := 0; i < b.N; i++ { + GetMoonRiseTime(jde, 115, 32, 8, 0, 10) + } + +} + func Test_MoonS(t *testing.T) { //fmt.Println(Sita(2451547)) //fmt.Println(MoonHeight(2451547, 115, 32, 8)) a := time.Now().UnixNano() - b := GetMoonRiseTime(GetNowJDE(), 115, 32, 8, 0) + b := GetMoonRiseTime(GetNowJDE(), 115, 32, 8, 0, 10) + fmt.Println(HMoonHeight(b, 115, 32, 8)) fmt.Println(time.Now().UnixNano() - a) fmt.Println(JDE2Date((b))) fmt.Println(time.Now().UnixNano() - a) //fmt.Printf("%.14f", GetMoonRiseTime(2451547, 115, 32, 8, 0)) -} \ No newline at end of file +} diff --git a/basic/star_test.go b/basic/star_test.go index 9264eff..4104337 100644 --- a/basic/star_test.go +++ b/basic/star_test.go @@ -7,10 +7,19 @@ import ( func Test_Isxz(t *testing.T) { now := GetNowJDE() - for i := 0.00; i <= 360.00; i+=0.5 { - for j := -90.00; j <= 90.00; j+=0.5 { + for i := 0.00; i <= 360.00; i += 0.5 { + for j := -90.00; j <= 90.00; j += 0.5 { fmt.Println(i, j) fmt.Println(WhichCst(float64(i), float64(j), now)) } } } + +func Benchmark_IsXZ(b *testing.B) { + jde := GetNowJDE() + for i := 0; i < b.N; i++ { + //GetNowJDE() + WhichCst(11.11, 12.12, jde) + } + +} \ No newline at end of file diff --git a/basic/sun.go b/basic/sun.go index 1c9ec3a..f562c18 100644 --- a/basic/sun.go +++ b/basic/sun.go @@ -1078,7 +1078,7 @@ func GetJQTime(Year, Angle int) float64 { //节气时间 stDegree := JQLospec(JD0) - float64(Angle) stDegreep := (JQLospec(JD0+0.000005) - JQLospec(JD0-0.000005)) / 0.00001 JD1 = JD0 - stDegree/stDegreep - if math.Floor(JD1-JD0) <= 0.000001 { + if math.Abs(JD1-JD0) <= 0.00001 { break } } @@ -1156,7 +1156,7 @@ func GetWHTime(Year, Angle int) float64 { stDegree := JQLospec(JD0) - float64(Angle) stDegreep := (JQLospec(JD0+0.000005) - JQLospec(JD0-0.000005)) / 0.00001 JD1 = JD0 - stDegree/stDegreep - if math.Floor(JD1-JD0) <= 0.000001 { + if math.Abs(JD1-JD0) <= 0.00001 { break } } @@ -1206,7 +1206,7 @@ func GetBanTime(JD, Lon, Lat, TZ, An float64) float64 { 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 - if math.Floor(JD1-JD0) < 0.000001 { + if math.Abs(JD1-JD0) < 0.00001 { break } } @@ -1244,7 +1244,7 @@ func GetAsaTime(JD, Lon, Lat, TZ, An float64) float64 { 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 - if math.Floor(JD1-JD0) < 0.000001 { + if math.Abs(JD1-JD0) < 0.00001 { break } } @@ -1266,13 +1266,14 @@ func SunTimeAngle(JD, Lon, Lat, TZ float64) float64 { /* * 精确计算，传入当日0时JDE */ -func GetSunRiseTime(JD, Lon, Lat, TZ, ZS float64) float64 { +func GetSunRiseTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 { var An float64 JD = math.Floor(JD) + 1.5 ntz := math.Round(Lon / 15) if ZS != 0 { An = -0.8333 } + An = An - HeightDegreeByLat(HEI, Lat) tztime := GetSunTZTime(JD, Lon, ntz) dec := HSunSeeDec(tztime) tmp := -Tan(Lat) * Tan(dec) @@ -1292,19 +1293,20 @@ func GetSunRiseTime(JD, Lon, Lat, TZ, ZS float64) float64 { 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 - if math.Floor(JD1-JD0) <= 0.000001 { + if math.Abs(JD1-JD0) <= 0.00001 { break } } return JD1 - ntz/24 + TZ/24 } -func GetSunDownTime(JD, Lon, Lat, TZ, ZS float64) float64 { +func GetSunDownTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 { var An float64 JD = math.Floor(JD) + 1.5 ntz := math.Round(Lon / 15) if ZS != 0 { An = -0.8333 } + An = An - HeightDegreeByLat(HEI, Lat) tztime := GetSunTZTime(JD, Lon, ntz) dec := HSunSeeDec(tztime) tmp := -Tan(Lat) * Tan(dec) @@ -1324,7 +1326,7 @@ func GetSunDownTime(JD, Lon, Lat, TZ, ZS float64) float64 { 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 - if math.Floor(JD1-JD0) <= 0.000001 { + if math.Abs(JD1-JD0) <= 0.00001 { break } } diff --git a/basic/sun_test.go b/basic/sun_test.go index a5c450a..e7a0d5e 100644 --- a/basic/sun_test.go +++ b/basic/sun_test.go @@ -22,6 +22,24 @@ func Test_SunLo(t *testing.T) { fmt.Printf("%.14f", HSunSeeLo(2458840.0134162)) } +func Benchmark_SunRise(b *testing.B) { + jde := GetNowJDE() + for i := 0; i < b.N; i++ { + //GetNowJDE() + GetSunRiseTime(jde, 115, 32, 8, 0, 10) + } + +} + +func Benchmark_SunLo(b *testing.B) { + jde := GetNowJDE() + for i := 0; i < b.N; i++ { + //GetNowJDE() + HSunSeeLo(jde) + } + +} + func Test_Cal(t *testing.T) { fmt.Println(JDE2Date(GetSolar(2020, 1, 1, false))) fmt.Println(JDE2Date(GetSolar(2020, 4, 1, false))) @@ -33,17 +51,18 @@ func Test_Cal(t *testing.T) { func Test_SunRise(t *testing.T) { a := time.Now().UnixNano() - b := GetSunRiseTime(GetNowJDE(), 115, 32, 8, 0) - b = GetSunRiseTime(GetNowJDE()+1, 115, 32, 8, 0) - b = GetSunRiseTime(GetNowJDE()+2, 115, 32, 8, 0) - b = GetSunRiseTime(GetNowJDE()+3, 115, 32, 8, 0) - b = GetSunRiseTime(GetNowJDE()+4, 115, 32, 8, 0) - b = GetSunRiseTime(GetNowJDE()+5, 115, 32, 8, 0) - b = GetSunRiseTime(GetNowJDE()+6, 115, 32, 8, 0) - b = GetSunRiseTime(GetNowJDE()+7, 115, 32, 8, 0) - b = GetSunRiseTime(GetNowJDE()+8, 115, 32, 8, 0) - b = GetSunRiseTime(GetNowJDE()+9, 115, 32, 8, 0) + //b := GetSunRiseTime(GetNowJDE(), 115, 32, 8, 0) + //b = GetSunRiseTime(GetNowJDE()+1, 115, 32, 8, 0) + //b = GetSunRiseTime(GetNowJDE()+2, 115, 32, 8, 0) + //b = GetSunRiseTime(GetNowJDE()+3, 115, 32, 8, 0) + //b = GetSunRiseTime(GetNowJDE()+4, 115, 32, 8, 0) + //b = GetSunRiseTime(GetNowJDE()+5, 115, 32, 8, 0) + //b = GetSunRiseTime(GetNowJDE()+6, 115, 32, 8, 0) + //b = GetSunRiseTime(GetNowJDE()+7, 115, 32, 8, 0) + //b = GetSunRiseTime(GetNowJDE()+8, 115, 32, 8, 0) + b := GetSunRiseTime(GetNowJDE()+9, 115, 32, 8, 0, 10) fmt.Println(time.Now().UnixNano() - a) + fmt.Println(SunHeight(b, 115, 32, 8)) fmt.Println(JDE2Date((b))) fmt.Println(time.Now().UnixNano() - a) } diff --git a/moon/moon.go b/moon/moon.go index 8076f60..0d72225 100644 --- a/moon/moon.go +++ b/moon/moon.go @@ -109,7 +109,9 @@ func CulminationTime(date time.Time, lon, lat float64) float64 { // date, 世界时（忽略此处时区） // lon，经度，东正西负 // lat，纬度，北正南负 -func RiseTime(date time.Time, lon, lat float64, aero bool) (time.Time, error) { +// height，高度 +// aero,是否进行大气修正 +func RiseTime(date time.Time, lon, lat, height float64, aero bool) (time.Time, error) { var err error if date.Hour() > 12 { date = date.Add(time.Hour * -12) @@ -121,7 +123,7 @@ func RiseTime(date time.Time, lon, lat float64, aero bool) (time.Time, error) { if aero { aeroFloat = 1 } - riseJde := basic.GetMoonRiseTime(jde, lon, lat, timezone, aeroFloat) + riseJde := basic.GetMoonRiseTime(jde, lon, lat, timezone, aeroFloat, height) if riseJde == -3 { err = ERR_NOT_TODAY } @@ -138,7 +140,9 @@ func RiseTime(date time.Time, lon, lat float64, aero bool) (time.Time, error) { // date, 世界时（忽略此处时区） // lon，经度，东正西负 // lat，纬度，北正南负 -func DownTime(date time.Time, lon, lat float64, aero bool) (time.Time, error) { +// height，高度 +// aero，大气修正 +func DownTime(date time.Time, lon, lat, height float64, aero bool) (time.Time, error) { var err error if date.Hour() > 12 { date = date.Add(time.Hour * -12) @@ -150,7 +154,7 @@ func DownTime(date time.Time, lon, lat float64, aero bool) (time.Time, error) { if aero { aeroFloat = 1 } - downJde := basic.GetMoonDownTime(jde, lon, lat, timezone, aeroFloat) + downJde := basic.GetMoonDownTime(jde, lon, lat, timezone, aeroFloat, height) if downJde == -3 { err = ERR_NOT_TODAY } diff --git a/moon/moon_test.go b/moon/moon_test.go index 2198af0..dbd8b6d 100644 --- a/moon/moon_test.go +++ b/moon/moon_test.go @@ -7,5 +7,5 @@ import ( ) func Test_Rise(t *testing.T) { - fmt.Println(RiseTime(time.Now(), 115, 32, true)) + fmt.Println(RiseTime(time.Now(), 120, 40, 10, true)) } diff --git a/sun/sun.go b/sun/sun.go index 8c3db9b..113ed4a 100644 --- a/sun/sun.go +++ b/sun/sun.go @@ -27,8 +27,9 @@ var ( // date，取日期，时区忽略 // lon，经度，东正西负 // lat，纬度，北正南负 +// height，高度 // aero，true时进行大气修正 -func RiseTime(date time.Time, lon, lat float64, aero bool) (time.Time, error) { +func RiseTime(date time.Time, lon, lat, height float64, aero bool) (time.Time, error) { var err error var aeroFloat float64 if aero { @@ -40,7 +41,7 @@ func RiseTime(date time.Time, lon, lat float64, aero bool) (time.Time, error) { jde := basic.Date2JDE(date) _, loc := date.Zone() timezone := float64(loc) / 3600.0 - riseJde := basic.GetSunRiseTime(jde, lon, lat, timezone, aeroFloat) + riseJde := basic.GetSunRiseTime(jde, lon, lat, timezone, aeroFloat, height) if riseJde == -2 { err = ERR_SUN_NEVER_RISE } @@ -54,8 +55,9 @@ func RiseTime(date time.Time, lon, lat float64, aero bool) (time.Time, error) { // date，当地时区日期，务必做时区修正 // lon，经度，东正西负 // lat，纬度，北正南负 +// height，高度 // aero，true时进行大气修正 -func DownTime(date time.Time, lon, lat float64, aero bool) (time.Time, error) { +func DownTime(date time.Time, lon, lat, height float64, aero bool) (time.Time, error) { var err error var aeroFloat float64 if aero { @@ -67,7 +69,7 @@ func DownTime(date time.Time, lon, lat float64, aero bool) (time.Time, error) { jde := basic.Date2JDE(date) _, loc := date.Zone() timezone := float64(loc) / 3600.0 - downJde := basic.GetSunDownTime(jde, lon, lat, timezone, aeroFloat) + downJde := basic.GetSunDownTime(jde, lon, lat, timezone, aeroFloat, height) if downJde == -2 { err = ERR_SUN_NEVER_RISE }