diff --git a/basic/moon.go b/basic/moon.go index 7fd6a2e..14e49d3 100644 --- a/basic/moon.go +++ b/basic/moon.go @@ -1443,19 +1443,36 @@ func GetMoonRiseTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 { if math.Abs(now-180) > 0.5 { JD1 += (180 - now) * 4.0 / 60.0 / 24.0 } - hei := MoonHeight(JD1, Lon, Lat, TZ) + hei := HMoonHeight(JD1, Lon, Lat, TZ) if !(hei < -10 && math.Abs(Lat) < 60) { if hei > An { return -1 //拱 } - if MoonHeight(JD1+12/24+6/15/24, Lon, Lat, TZ) < An { + reJde := JD1 + 12.0/24.0 + 6.0/15.0/24.0 + mag := MoonTimeAngle(reJde, Lon, Lat, TZ) + if mag < 90 { + mag += 360 + } + reJde += (360 - mag) * 4.0 / 60.0 / 24.0 + if HMoonHeight(reJde, Lon, Lat, TZ) < An { return -2 //沉 } } dec := MoonSeeDec(JD1, Lon, Lat, TZ) - SJ := (180 - ArcCos(-Tan(Lat)*Tan(dec))) / 15 - JD1 += SJ/24.00 + SJ/33.00/15.00 - + tmp := (Sin(An) - Sin(dec)*Sin(Lat)) / (Cos(dec) * Cos(Lat)) + if math.Abs(tmp) <= 1 && Lat < 85 { + SJ := (180 - ArcCos(tmp)) / 15 + JD1 += SJ/24.00 + SJ/33.00/15.00 + } else { + i := 0 + for MoonHeight(JD1, Lon, Lat, TZ) < An { + i++ + JD1 += 15.0 / 60.0 / 24.0 + if i > 48 { + break + } + } + } for { JD0 := JD1 stDegree := HMoonHeight(JD0, Lon, Lat, TZ) - An @@ -1487,10 +1504,11 @@ func GetMoonDownTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 { } An = An - HeightDegreeByLat(HEI, Lat) moonang := MoonTimeAngle(JD, Lon, Lat, TZ) - if moonheight < 0 { //月亮在地平线上或在落下与下中天之间 + if moonheight < 0 { tms = (360 - moonang) / 15 JD1 += (tms/24 + (tms/24.0*12.0)/15.0/24.0) } + //月亮在地平线上或在落下与下中天之间 if moonheight > 0 && moonang < 180 { tms = (-moonang) / 15 JD1 += (tms/24.0 + (tms/24.0*12.0)/15.0/24.0) @@ -1505,18 +1523,34 @@ func GetMoonDownTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 { if math.Abs(now-360) > 0.5 { JD1 += (360 - now) * 4.0 / 60.0 / 24.0 } - hei := MoonHeight(JD1, Lon, Lat, TZ) + //JD1=月球中天时间 + hei := HMoonHeight(JD1, Lon, Lat, TZ) if !(hei > 10 && math.Abs(Lat) < 60) { if hei < An { return -2 //沉 } - if MoonHeight(JD1+12.0/24.0+6.0/15.0/24.0, Lon, Lat, TZ) > An { + reJde := JD1 + 12.0/24.0 + 6.0/15.0/24.0 + sub := 180 - MoonTimeAngle(reJde, Lon, Lat, TZ) + reJde += sub * 4.0 / 60.0 / 24.0 + if HMoonHeight(reJde, Lon, Lat, TZ) > An { return -1 //拱 } } dec := MoonSeeDec(JD1, Lon, Lat, TZ) - SJ := (ArcCos(-Tan(Lat) * Tan(dec))) / 15.0 - JD1 += SJ/24 + SJ/33.0/15.0 + tmp := (Sin(An) - Sin(dec)*Sin(Lat)) / (Cos(dec) * Cos(Lat)) + if math.Abs(tmp) <= 1 && Lat < 85 { + SJ := (ArcCos(tmp)) / 15.0 + JD1 += SJ/24 + SJ/33.0/15.0 + } else { + i := 0 + for MoonHeight(JD1, Lon, Lat, TZ) > An { + i++ + JD1 += 15.0 / 60.0 / 24.0 + if i > 48 { + break + } + } + } for { JD0 := JD1 stDegree := HMoonHeight(JD0, Lon, Lat, TZ) - An diff --git a/basic/moon_test.go b/basic/moon_test.go index 46c4987..0eccde2 100644 --- a/basic/moon_test.go +++ b/basic/moon_test.go @@ -14,6 +14,26 @@ func Benchmark_MoonRiseBench(b *testing.B) { } +func Test_MoonDown(t *testing.T) { + jde := GetNowJDE() + for i := 30.0; i < 90.0; i += 0.3 { + fmt.Println(i, GetMoonDownTime(jde, 115, float64(i), 8, 1, 0)) + } +} + +func Test_MoonRise(t *testing.T) { + //2.459984692085961e+06 113.58880556 87.36833333 8 0 0 + //2.459984692085961e+06 113.58880556 87.36833333 8 0 0 + //2.4599846948519214e+06 113.58880556 87.36833333 8 0 0 + + //cst := time.FixedZone("cst", 8*3600) + //jde := Date2JDE(time.Date(2023, 2, 9, 15, 59, 0, 0, cst)) + fmt.Println(GetMoonRiseTime(2.4599846948519214e+06, 113.58880556, 87.36833333, 8, 0, 0)) + for i := 30.0; i < 90.0; i += 0.3 { + fmt.Println(i, GetMoonRiseTime(2.459984692085961e+06, 113.588, float64(i), 8, 0, 0)) + } +} + func Test_MoonS(t *testing.T) { //fmt.Println(Sita(2451547)) //fmt.Println(MoonHeight(2451547, 115, 32, 8)) diff --git a/basic/star.go b/basic/star.go index fc096f7..735a8e1 100644 --- a/basic/star.go +++ b/basic/star.go @@ -2,6 +2,7 @@ package basic import ( . "b612.me/astro/tools" + "math" ) // StarHeight 星体的高度角 @@ -90,7 +91,70 @@ func StarAngle(RA, DEC, JD, Lon, Lat, TZ float64) float64 { } } -func StarRiseTime(jde, ra, dec, lon, lat, height, timezone float64, aero bool) (float64, error) { +func StarRiseTime(jde, ra, dec, lon, lat, height, timezone float64, aero bool) float64 { + return StarRiseDownTime(jde, ra, dec, lon, lat, height, timezone, aero, true) +} + +func StarDownTime(jde, ra, dec, lon, lat, height, timezone float64, aero bool) float64 { + return StarRiseDownTime(jde, ra, dec, lon, lat, height, timezone, aero, false) +} - return 0, nil +func StarRiseDownTime(jde, ra, dec, lon, lat, height, timezone float64, aero, isRise bool) float64 { + //jde 世界时，非力学时，当地时区 0时，无需转换力学时 + //ra,dec 瞬时天球座标，非J2000等时间天球坐标 + jde = math.Floor(jde) + 0.5 + var An float64 = 0 + if aero { + An = -0.566667 + } + An = An - HeightDegreeByLat(height, lat) + sct := StarCulminationTime(jde, ra, lon, timezone) + tmp := (Sin(An) - Sin(dec)*Sin(lat)) / (Cos(dec) * Cos(lat)) + if math.Abs(tmp) > 1 { + if StarHeight(sct, ra, dec, lon, lat, timezone) < 0 { + return -2 //极夜 + } else { + return -1 //极昼 + } + } + var JD1 float64 + if isRise { + JD1 = sct - ArcCos(tmp)/15.0/24.0 + } else { + JD1 = sct + ArcCos(tmp)/15.0/24.0 + } + for { + JD0 := JD1 + stDegree := StarHeight(JD0, ra, dec, lon, lat, timezone) - An + stDegreep := (StarHeight(JD0+0.000005, ra, dec, lon, lat, timezone) - StarHeight(JD0-0.000005, ra, dec, lon, lat, timezone)) / 0.00001 + JD1 = JD0 - stDegree/stDegreep + if math.Abs(JD1-JD0) <= 0.00001 { + break + } + } + return JD1 +} + +func StarCulminationTime(jde, ra, lon, timezone float64) float64 { + //jde 世界时，非力学时，当地时区 0时，无需转换力学时 + //ra,dec 瞬时天球座标，非J2000等时间天球坐标 + jde = math.Floor(jde) + 0.5 + JD1 := jde + Limit360(360-StarHourAngle(jde, ra, lon, timezone))/15.0/24.0*0.99726851851851851851 + limitStarHA := func(jde, ra, lon, timezone float64) float64 { + ha := StarHourAngle(jde, ra, lon, timezone) + if ha < 180 { + ha += 360 + } + return ha + } + for { + JD0 := JD1 + stDegree := limitStarHA(JD0, ra, lon, timezone) - 360 + stDegreep := (limitStarHA(JD0+0.000005, ra, lon, timezone) - SunHeight(JD0-0.000005, ra, lon, timezone)) / 0.00001 + JD1 = JD0 - stDegree/stDegreep + if math.Abs(JD1-JD0) <= 0.00001 { + break + } + } + return JD1 } diff --git a/basic/star_test.go b/basic/star_test.go index bbadba7..b85311f 100644 --- a/basic/star_test.go +++ b/basic/star_test.go @@ -8,7 +8,15 @@ import ( func Test_StarHeight(t *testing.T) { date := GetNowJDE() + 6.0/24.0 fmt.Println(JDE2Date(date)) - fmt.Println("Sirius Height:", StarHeight(date, 101.5, -16.8, 113.5, 22, 8.0)) - fmt.Println("Sirius Azimuth:", StarAzimuth(date, 101.5, -16.8, 113.5, 22, 8.0)) + fmt.Println("Sirius Height:", StarHeight(date, 101.5, -16.8, 115, 40, 8.0)) + fmt.Println("Sirius Azimuth:", StarAzimuth(date, 101.5, -16.8, 115, 40, 8.0)) } + +func Test_Star(t *testing.T) { + date := GetNowJDE() - 20/24.0 + fmt.Println(JDE2Date(date)) + fmt.Println("Sirius RiseTime:", JDE2Date(StarRiseTime(date, 101.529, -16.8, 113.568, 22.5, 0, 8.0, true))) + fmt.Println("Sirius CulminationTime:", JDE2Date(StarCulminationTime(date, 101.529, 113.568, 8.0))) + fmt.Println("Sirius DownTime:", JDE2Date(StarDownTime(date, 101.529, -16.8, 113.568, 22.5, 0, 8.0, true))) +} diff --git a/basic/stardat_test.go b/basic/stardat_test.go index d97b311..e827fd9 100644 --- a/basic/stardat_test.go +++ b/basic/stardat_test.go @@ -1,19 +1,20 @@ package basic import ( - "fmt" - "io/ioutil" "testing" ) func Test_ParseStar(t *testing.T) { //dat := []byte(`2491 9Alp CMaBD-16 1591 48915151881 257I 5423 064044.6-163444064508.9-164258227.22-08.88-1.46 0.00 -0.05 -0.03 A1Vm -0.553-1.205 +.375-008SBO 13 10.3 11.2AB 4*`) - LoadStarData() - out := "" + err := LoadStarData() + if err != nil { + t.Fatal(err) + } for _, v := range stardat { - data, _ := parseStarData(v) - out += fmt.Sprintln(data.HD, ",", data.HR) + _, err = parseStarData(v) + if err != nil { + t.Fatal(err) + } } - ioutil.WriteFile("./index.csv", []byte(out), 0644) } diff --git a/basic/sun.go b/basic/sun.go index 6c17de7..db5fd5b 100644 --- a/basic/sun.go +++ b/basic/sun.go @@ -1179,28 +1179,29 @@ func GetBanTime(JD, Lon, Lat, TZ, An float64) float64 { JD = math.Floor(JD) + 1.5 ntz := math.Round(Lon / 15) tztime := GetSunTZTime(JD, Lon, ntz) - dec := HSunSeeDec(tztime) - tmp := -Tan((math.Abs(Lat)+An)*(Lat/math.Abs(Lat))) * Tan(dec) - if math.Abs(tmp) > 1 { - if SunHeight(tztime, Lon, Lat, ntz) < An { - return -2 //极夜 - } - if SunHeight(tztime-0.5, Lon, Lat, ntz) > An { - return -1 //极昼 - } + if SunHeight(tztime, Lon, Lat, ntz) < An { + return -2 //极夜 } - tmp = -Tan(Lat) * Tan(dec) - rzsc := ArcCos(tmp) / 15 - sunrise := tztime + rzsc/24.0 + 35.0/24.0/60.0 - i := 0 - for LowSunHeight(sunrise, Lon, Lat, ntz) < An { - i++ - sunrise -= 15 / 60 / 24 - if i > 12 { - break + if SunHeight(tztime+0.5, Lon, Lat, ntz) > An { + return -1 //极昼 + } + tmp := (Sin(An) - Sin(HSunSeeDec(tztime))*Sin(Lat)) / (Cos(HSunSeeDec(tztime)) * 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 + } else { + sundown = tztime + i := 0 + for LowSunHeight(sundown, Lon, Lat, ntz) > An { + i++ + sundown += 15.0 / 60.0 / 24.0 + if i > 48 { + break + } } } - JD1 := sunrise - 5.00/24.00/60.00 + JD1 := sundown - 5.00/24.00/60.00 for { JD0 := JD1 stDegree := SunHeight(JD0, Lon, Lat, ntz) - An @@ -1217,25 +1218,26 @@ func GetAsaTime(JD, Lon, Lat, TZ, An float64) float64 { JD = math.Floor(JD) + 1.5 ntz := math.Round(Lon / 15) tztime := GetSunTZTime(JD, Lon, ntz) - dec := HSunSeeDec(tztime) - tmp := -Tan((math.Abs(Lat)+An)*(Lat/math.Abs(Lat))) * Tan(dec) - if math.Abs(tmp) > 1 { - if SunHeight(tztime, Lon, Lat, ntz) < An { - return -2 //极夜 - } - if SunHeight(tztime-0.5, Lon, Lat, ntz) > An { - return -1 //极昼 - } + if SunHeight(tztime, Lon, Lat, ntz) < An { + return -2 //极夜 } - tmp = -Tan(Lat) * Tan(dec) - rzsc := ArcCos(tmp) / 15 - sunrise := tztime - rzsc/24 - 25.0/24.0/60.0 - i := 0 - for LowSunHeight(sunrise, Lon, Lat, ntz) > An { - i++ - sunrise -= 15 / 60 / 24 - if i > 12 { - break + if SunHeight(tztime-0.5, Lon, Lat, ntz) > An { + return -1 //极昼 + } + tmp := (Sin(An) - Sin(HSunSeeDec(tztime))*Sin(Lat)) / (Cos(HSunSeeDec(tztime)) * 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 + } else { + sunrise = tztime + i := 0 + for LowSunHeight(sunrise, Lon, Lat, ntz) > An { + i++ + sunrise -= 15.0 / 60.0 / 24.0 + if i > 48 { + break + } } } JD1 := sunrise - 5.00/24.00/60.00 @@ -1275,18 +1277,30 @@ func GetSunRiseTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 { } An = An - HeightDegreeByLat(HEI, Lat) tztime := GetSunTZTime(JD, Lon, ntz) - dec := HSunSeeDec(tztime) - tmp := -Tan(Lat) * Tan(dec) - if math.Abs(tmp) > 1 { - if SunHeight(tztime, Lon, Lat, ntz) < 0 { - return -2 //极夜 - } - if SunHeight(tztime-0.5, Lon, Lat, ntz) > 0 { - return -1 //极昼 + if SunHeight(tztime, Lon, Lat, ntz) < An { + return -2 //极夜 + } + if SunHeight(tztime-0.5, Lon, Lat, ntz) > An { + return -1 //极昼 + } + //(sin(ho)-sin(φ)*sin(δ2))/(cos(φ)*cos(δ2)) + tmp := (Sin(An) - Sin(HSunSeeDec(tztime))*Sin(Lat)) / (Cos(HSunSeeDec(tztime)) * 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 + } else { + sunrise = tztime + i := 0 + //TODO:使用二分法计算 + for LowSunHeight(sunrise, Lon, Lat, ntz) > An { + i++ + sunrise -= 15.0 / 60.0 / 24.0 + if i > 48 { + break + } } } - rzsc := ArcCos(tmp) / 15 - sunrise := tztime - rzsc/24 - 5.0/24.0/60.0 JD1 := sunrise for { JD0 := JD1 @@ -1308,19 +1322,29 @@ func GetSunDownTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 { } An = An - HeightDegreeByLat(HEI, Lat) tztime := GetSunTZTime(JD, Lon, ntz) - dec := HSunSeeDec(tztime) - tmp := -Tan(Lat) * Tan(dec) - if math.Abs(tmp) > 1 { - if SunHeight(tztime, Lon, Lat, ntz) < 0 { - return -2 //极夜 - } - if SunHeight(tztime+0.5, Lon, Lat, ntz) > 0 { - return -1 //极昼 + if SunHeight(tztime, Lon, Lat, ntz) < An { + return -2 //极夜 + } + if SunHeight(tztime+0.5, Lon, Lat, ntz) > An { + return -1 //极昼 + } + tmp := (Sin(An) - Sin(HSunSeeDec(tztime))*Sin(Lat)) / (Cos(HSunSeeDec(tztime)) * 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 + } else { + sundown = tztime + i := 0 + for LowSunHeight(sundown, Lon, Lat, ntz) > An { + i++ + sundown += 15.0 / 60.0 / 24.0 + if i > 48 { + break + } } } - rzsc := ArcCos(tmp) / 15.0 - sunrise := tztime + rzsc/24 - 5.0/24.0/60.0 - JD1 := sunrise + JD1 := sundown for { JD0 := JD1 stDegree := SunHeight(JD0, Lon, Lat, ntz) - An diff --git a/basic/sun_test.go b/basic/sun_test.go index bd69d4d..a762ba8 100644 --- a/basic/sun_test.go +++ b/basic/sun_test.go @@ -66,3 +66,40 @@ func Test_SunRise(t *testing.T) { fmt.Println(JDE2Date((b))) fmt.Println(time.Now().UnixNano() - a) } + +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)) + + for i := 10.0; i < 90.0; i += 0.3 { + fmt.Println(i, GetAsaTime(jde, 113.5880556, float64(i), 8, -6)) + } +} + +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)) + } +} + +func Test_SunRiseRound(t *testing.T) { + jde := GetNowJDE() + for i := 10.0; i < 90.0; i += 0.3 { + fmt.Println(i, GetSunRiseTime(jde, 115, float64(i), 8, 0, 0)) + } +} +func Test_SunDown(t *testing.T) { + jde := GetNowJDE() + for i := 10.0; i < 90.0; i += 0.3 { + fmt.Println(i, GetSunDownTime(jde, 115, float64(i), 8, 0, 0)) + } +} + +func Test_SunAz(t *testing.T) { + cst := time.FixedZone("cst", 8*3600) + fmt.Println(SunAngle(Date2JDE(time.Date(2022, 5, 30, 11, 55, 0, 0, cst)), + 120, 30, 8.0)) +} diff --git a/sun/sun.go b/sun/sun.go index 113ed4a..d034e7b 100644 --- a/sun/sun.go +++ b/sun/sun.go @@ -86,6 +86,9 @@ func DownTime(date time.Time, lon, lat, height float64, aero bool) (time.Time, e // angle，朦影角度：可选-6 -12 -18(民用、航海、天文) func MorningTwilight(date time.Time, lon, lat, angle float64) (time.Time, error) { var err error + if date.Hour() > 12 { + date = date.Add(time.Hour * -12) + } jde := basic.Date2JDE(date) _, loc := date.Zone() timezone := float64(loc) / 3600.0 @@ -106,6 +109,9 @@ func MorningTwilight(date time.Time, lon, lat, angle float64) (time.Time, error) // angle，朦影角度：可选-6 -12 -18(民用、航海、天文) func EveningTwilight(date time.Time, lon, lat, angle float64) (time.Time, error) { var err error + if date.Hour() > 12 { + date = date.Add(time.Hour * -12) + } jde := basic.Date2JDE(date) _, loc := date.Zone() timezone := float64(loc) / 3600.0