b612
9ee2163cc7
- 新增月掩恒星和行星：支持搜索、掩甚点、全球掩带及固定地点轨迹计算 - 支持恒星星表坐标转换、有限盘面行星接触事件和月掩 SVG 输出 - 新增日月食及月掩全球投影图、时间标记和 GeoJSON 地理数据接口 - 扩展日食中心线、南北界及偏食足迹采样，支持极区投影 - 修正站心时角、月出月落、月球视半径、折射和恒星自行计算 - 优化内外行星事件搜索、边界选择、极端输入处理和计算稳定性
109 lines
2.4 KiB
Go
109 lines
2.4 KiB
Go
package basic
|
|
|
|
import "math"
|
|
|
|
func GetMoonLoops(year float64, loop int) []float64 {
|
|
var start float64
|
|
var newMoon, lastNewMoon float64
|
|
moonLoops := make([]float64, loop)
|
|
if year < 6000 {
|
|
start = year + 11.00/12.00 + 5.00/30.00/12.00
|
|
} else {
|
|
start = year + 9.00/12.00 + 5.00/30.00/12.00
|
|
}
|
|
i := 1
|
|
for j := 0; j < loop; j++ {
|
|
if year > 3000 {
|
|
newMoon = TD2UT(CalcMoonSH(start+float64(i-1)/12.5, 0)+8.0/24.0, false)
|
|
} else {
|
|
newMoon = TD2UT(CalcMoonS(start+float64(i-1)/12.5, 0)+8.0/24.0, false)
|
|
}
|
|
if i != 1 {
|
|
if newMoon == lastNewMoon {
|
|
j--
|
|
i++
|
|
continue
|
|
}
|
|
}
|
|
moonLoops[j] = newMoon
|
|
lastNewMoon = moonLoops[j]
|
|
i++
|
|
}
|
|
return moonLoops
|
|
}
|
|
|
|
func GetJieqiLoops(year, loop int) []float64 {
|
|
start := 270
|
|
jq := make([]float64, loop)
|
|
for i := 1; i <= loop; i++ {
|
|
angle := start + 15*(i-1)
|
|
if angle > 360 {
|
|
angle -= 360
|
|
}
|
|
jq[i-1] = GetJQTime(year+int(math.Ceil(float64(i-1)/24.000)), angle) + 8.0/24.0
|
|
}
|
|
return jq
|
|
}
|
|
|
|
func GetJQTime(year, angle int) float64 {
|
|
// Calculate initial day based on angle parity
|
|
var initialDay float64
|
|
if angle%2 == 0 {
|
|
initialDay = 18
|
|
} else {
|
|
initialDay = 3
|
|
}
|
|
|
|
// Calculate temporary factor for month offset
|
|
var tempFactor float64
|
|
if angle%10 != 0 {
|
|
tempFactor = float64(angle+15) / 30.0
|
|
} else {
|
|
tempFactor = float64(angle) / 30.0
|
|
}
|
|
|
|
// Calculate initial month, adjusting if超过 12
|
|
initialMonth := 3.0 + tempFactor
|
|
if initialMonth > 12.0 {
|
|
initialMonth -= 12.0
|
|
}
|
|
|
|
// Calculate initial Julian date
|
|
initialJD := JDECalc(year, int(initialMonth), initialDay)
|
|
|
|
// Set target angle for iteration; if angle is 0, use 360
|
|
targetAngle := float64(angle)
|
|
if angle == 0 {
|
|
targetAngle = 360.0
|
|
}
|
|
|
|
// Newton-Raphson iteration to find precise Julian date
|
|
currentJD := initialJD
|
|
var ok bool
|
|
currentJD, ok = eventNewtonRefine(currentJD, 0.00001, func(previousJD float64) float64 {
|
|
errorValue := JQLospec(previousJD, targetAngle) - targetAngle
|
|
derivative := (JQLospec(previousJD+0.000005, targetAngle) - JQLospec(previousJD-0.000005, targetAngle)) / 0.00001
|
|
return errorValue / derivative
|
|
})
|
|
if !ok {
|
|
return math.NaN()
|
|
}
|
|
|
|
// Convert to UT and return
|
|
return TD2UT(currentJD, false)
|
|
}
|
|
|
|
func JQLospec(jd float64, target float64) float64 {
|
|
sunLo := HSunApparentLo(jd)
|
|
if target >= 345 {
|
|
if sunLo <= 12 {
|
|
sunLo += 360
|
|
}
|
|
} else if target <= 15 {
|
|
if sunLo >= 350 {
|
|
sunLo -= 360
|
|
}
|
|
}
|
|
return sunLo
|
|
}
|