添加海拔因素

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(GetSunRiseTime(jde, 115, 32, 8, 1, 10)))
-	fmt.Println("当前日落：", JDE2Date(GetSunDownTime(jde, 115, 32, 8, 1)))
+	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(GetMoonRiseTime(jde, 115, 32, 8, 1, 10)))
-	fmt.Println("当前月落：", JDE2Date(GetMoonDownTime(jde, 115, 32, 8, 1)))
+	fmt.Println("当前月落：", JDE2Date(GetMoonDownTime(jde, 115, 32, 8, 1, 10)))
 	fmt.Println("月相：", MoonLight(jde-8.0/24.0))
 }

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))
+}

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))
+}

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
 		}
 	}

basic/moon_test.go

@@ -6,11 +6,20 @@ 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)

basic/star_test.go

@@ -14,3 +14,12 @@ func Test_Isxz(t *testing.T) {
 		}
 	}
 }
+
+func Benchmark_IsXZ(b *testing.B) {
+	jde := GetNowJDE()
+	for i := 0; i < b.N; i++ {
+		//GetNowJDE()
+		WhichCst(11.11, 12.12, jde)
+	}
+
+}

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
 		}
 	}

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(), 115, 32, 8, 0)
-	b = GetSunRiseTime(GetNowJDE()+1, 115, 32, 8, 0)
+	//b = GetSunRiseTime(GetNowJDE()+1, 115, 32, 8, 0)
-	b = GetSunRiseTime(GetNowJDE()+2, 115, 32, 8, 0)
+	//b = GetSunRiseTime(GetNowJDE()+2, 115, 32, 8, 0)
-	b = GetSunRiseTime(GetNowJDE()+3, 115, 32, 8, 0)
+	//b = GetSunRiseTime(GetNowJDE()+3, 115, 32, 8, 0)
-	b = GetSunRiseTime(GetNowJDE()+4, 115, 32, 8, 0)
+	//b = GetSunRiseTime(GetNowJDE()+4, 115, 32, 8, 0)
-	b = GetSunRiseTime(GetNowJDE()+5, 115, 32, 8, 0)
+	//b = GetSunRiseTime(GetNowJDE()+5, 115, 32, 8, 0)
-	b = GetSunRiseTime(GetNowJDE()+6, 115, 32, 8, 0)
+	//b = GetSunRiseTime(GetNowJDE()+6, 115, 32, 8, 0)
-	b = GetSunRiseTime(GetNowJDE()+7, 115, 32, 8, 0)
+	//b = GetSunRiseTime(GetNowJDE()+7, 115, 32, 8, 0)
-	b = GetSunRiseTime(GetNowJDE()+8, 115, 32, 8, 0)
+	//b = GetSunRiseTime(GetNowJDE()+8, 115, 32, 8, 0)
-	b = GetSunRiseTime(GetNowJDE()+9, 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)
 }

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
 	}

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))
 }

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
 	}