b612/astro
1
1
Fork 0
Code Issues Pull Requests Projects Releases 18 Wiki Activity

change delta T formula and add planet phenomena function

Browse Source
This commit is contained in:
兔子
2022-05-16 20:42:15 +08:00
parent bb07e23238
commit 97c017be08
38 changed files with 4700 additions and 779 deletions
Download Patch File Download Diff File Expand all files Collapse all files
basic/calendar.go
+2 -2
View File
@@ -104,8 +104,8 @@ func DeltaT(Date float64, IsJDE bool) (Result float64) { //传入年或儒略日
return
}
if Year < 2100 && Year >= 2010 {
//fmt.Println(Year)
Result = dt_cal(Year) //-3.2-(Year-2017)*0.029915;
var t = (Year - 2000.0)
Result = 62.92 + 0.32217*t + 0.005589*t*t
return
}
if Year >= 2100 && Year <= 2150 {
basic/coordinate.go
+56 -24
View File
@@ -9,8 +9,8 @@ import (
/*
* 坐标变换，黄道转赤道
*/
func LoToRa(lo, bo, jde float64) float64 {
ra := math.Atan2(Sin(lo)*Cos(Sita(jde)-Tan(bo)*Sin(Sita(jde))), Cos(lo))
func LoToRa(jde, lo, bo float64) float64 {
ra := math.Atan2(Sin(lo)*Cos(Sita(jde))-Tan(bo)*Sin(Sita(jde)), Cos(lo))
ra = ra * 180 / math.Pi
if ra < 0 {
ra += 360
@@ -18,11 +18,48 @@ func LoToRa(lo, bo, jde float64) float64 {
return ra
}
func BoToDec(lo, bo, jde float64) float64 {
func BoToDec(jde, lo, bo float64) float64 {
dec := ArcSin(Sin(bo)*Cos(Sita(jde)) + Cos(bo)*Sin(Sita(jde))*Sin(lo))
return dec
}
func LoBoToRaDec(jde, lo, bo float64) (float64, float64) {
dec := ArcSin(Sin(bo)*Cos(Sita(jde)) + Cos(bo)*Sin(Sita(jde))*Sin(lo))
ra := math.Atan2(Sin(lo)*Cos(Sita(jde))-Tan(bo)*Sin(Sita(jde)), Cos(lo))
ra = ra * 180 / math.Pi
if ra < 0 {
ra += 360
}
return ra, dec
}
func RaDecToLoBo(jde, ra, dec float64) (float64, float64) {
//tan(λ) = (sin(α)*cos(ε) + tan(δ)*sin(ε)) / cos(α)
//sin(β)=sin(δ)*cos(ε)-cos(δ)*sin(ε)*sin(α)
sita := Sita(jde)
sinBo := Sin(dec)*Cos(sita) - Cos(dec)*Sin(sita)*Sin(ra)
lo := math.Atan2((Sin(ra)*Cos(sita) + Tan(dec)*Sin(sita)), Cos(ra))
lo = Limit360(lo * 180 / math.Pi)
return lo, ArcSin(sinBo)
}
func RaToLo(jde, ra, dec float64) float64 {
//tan(λ) = (sin(α)*cos(ε) + tan(δ)*sin(ε)) / cos(α)
//sin(β)=sin(δ)*cos(ε)-cos(δ)*sin(ε)*sin(α)
sita := Sita(jde)
lo := math.Atan2((Sin(ra)*Cos(sita) + Tan(dec)*Sin(sita)), Cos(ra))
lo = Limit360(lo * 180 / math.Pi)
return lo
}
func DecToBo(jde, ra, dec float64) float64 {
//tan(λ) = (sin(α)*cos(ε) + tan(δ)*sin(ε)) / cos(α)
//sin(β)=sin(δ)*cos(ε)-cos(δ)*sin(ε)*sin(α)
sita := Sita(jde)
sinBo := Sin(dec)*Cos(sita) - Cos(dec)*Sin(sita)*Sin(ra)
return ArcSin(sinBo)
}
/*
* 赤道坐标岁差变换st end 为JDE时刻
*/
@@ -98,7 +135,7 @@ func ZhanXinLo(lo, bo, lat, lon, jd, au, h float64) float64 { //jd为格林尼
C := pcosi(lat, h)
S := psini(lat, h)
sinpi := Sin(0.0024427777777) / au
ra := LoToRa(lo, bo, jd)
ra := LoToRa(jd, lo, bo)
tH := Limit360(TD2UT(ApparentSiderealTime(jd), false)*15 + lon - ra)
N := Cos(lo)*Cos(bo) - C*sinpi*Cos(tH)
nlo := math.Atan2(Sin(lo)*Cos(bo)-sinpi*(S*Sin(Sita(jd))+C*Cos(Sita(jd))*Sin(tH)), N) * 180 / math.Pi
@@ -109,7 +146,7 @@ func ZhanXinBo(lo, bo, lat, lon, jd, au, h float64) float64 { //jd为格林尼
C := pcosi(lat, h)
S := psini(lat, h)
sinpi := Sin(0.0024427777777) / au
ra := LoToRa(lo, bo, jd)
ra := LoToRa(jd, lo, bo)
tH := Limit360(TD2UT(ApparentSiderealTime(jd), false)*15 + lon - ra)
N := Cos(lo)*Cos(bo) - C*sinpi*Cos(tH)
nlo := math.Atan2(Sin(lo)*Cos(bo)-sinpi*(S*Sin(Sita(jd))+C*Cos(Sita(jd))*Sin(tH)), N) * 180 / math.Pi
@@ -117,25 +154,20 @@ func ZhanXinBo(lo, bo, lat, lon, jd, au, h float64) float64 { //jd为格林尼
return nbo
}
/*
func GXCLo(lo,bo,jd float64) float64{ //光行差修正
k:=20.49552;
sunlo:=SunTrueLo(jd);
e:=Earthe(jd);
epi=earth->EarthPI(jd);
tmp=(-k*this->CosR(sunlo-lo)+e*k*this->CosR(epi-lo))/this->CosR(bo);
return tmp;
func GXCLo(lo, bo, jd float64) float64 { //光行差修正
k := 20.49552
sunlo := SunTrueLo(jd)
e := Earthe(jd)
epi := EarthPI(jd)
tmp := (-k*Cos(sunlo-lo) + e*k*Cos(epi-lo)) / Cos(bo)
return tmp
}
public function GXCBo(lo,bo,jd)
{
earth=new Earth();
k=20.49552;
sunlo=earth->SunTrueLo(jd);
e=earth->Earthe(jd);
epi=earth->EarthPI(jd);
tmp=-k*this->SinR(bo)*(this->SinR(sunlo-lo)-e*this->SinR(epi-lo));
return tmp;
func GXCBo(lo, bo, jd float64) float64 {
k := 20.49552
sunlo := SunTrueLo(jd)
e := Earthe(jd)
epi := EarthPI(jd)
tmp := -k * Sin(bo) * (Sin(sunlo-lo) - e*Sin(epi-lo))
return tmp
}
*/
basic/coordinate_test.go
+40
View File
@@ -8,3 +8,43 @@ import (
func Test_LoBo(t *testing.T) {
fmt.Printf("%.9f", dt_cal(2020.5))
}
func Test_LoBoRaDec(t *testing.T) {
jde := 2451545.0
lo, bo := RaDecToLoBo(jde, 10, 50)
fmt.Println("LO,BO", lo, bo)
fmt.Println(LoBoToRaDec(jde, lo, bo))
lo, bo = RaDecToLoBo(jde, 40, 80)
fmt.Println("LO,BO", lo, bo)
fmt.Println(LoBoToRaDec(jde, lo, bo))
lo, bo = RaDecToLoBo(jde, 90, 50)
fmt.Println("LO,BO", lo, bo)
fmt.Println(LoBoToRaDec(jde, lo, bo))
lo, bo = RaDecToLoBo(jde, 130, 50)
fmt.Println("LO,BO", lo, bo)
fmt.Println(LoBoToRaDec(jde, lo, bo))
lo, bo = RaDecToLoBo(jde, 160, 50)
fmt.Println("LO,BO", lo, bo)
fmt.Println(LoBoToRaDec(jde, lo, bo))
lo, bo = RaDecToLoBo(jde, 180, 50)
fmt.Println("LO,BO", lo, bo)
fmt.Println(LoBoToRaDec(jde, lo, bo))
lo, bo = RaDecToLoBo(jde, 210, 50)
fmt.Println("LO,BO", lo, bo)
fmt.Println(LoBoToRaDec(jde, lo, bo))
lo, bo = RaDecToLoBo(jde, 260, 50)
fmt.Println("LO,BO", lo, bo)
fmt.Println(LoBoToRaDec(jde, lo, bo))
lo, bo = RaDecToLoBo(jde, 270, 50)
fmt.Println("LO,BO", lo, bo)
fmt.Println(LoBoToRaDec(jde, lo, bo))
lo, bo = RaDecToLoBo(jde, 300, 50)
fmt.Println("LO,BO", lo, bo)
fmt.Println(LoBoToRaDec(jde, lo, bo))
lo, bo = RaDecToLoBo(jde, 350, 50)
fmt.Println("LO,BO", lo, bo)
fmt.Println(LoBoToRaDec(jde, lo, bo))
lo, bo = RaDecToLoBo(jde, 0, 50)
fmt.Println("LO,BO", lo, bo)
fmt.Println(LoBoToRaDec(jde, lo, bo))
}
basic/jupiter.go
+299 -9
View File
@@ -63,22 +63,22 @@ func AJupiterXYZ(JD float64) (float64, float64, float64) {
return x, y, z
}
func JupiterSeeRa(JD float64) float64 {
lo, bo := JupiterSeeLoBo(JD)
func JupiterApparentRa(JD float64) float64 {
lo, bo := JupiterApparentLoBo(JD)
sita := Sita(JD)
ra := math.Atan2((Sin(lo)*Cos(sita) - Tan(bo)*Sin(sita)), Cos(lo))
ra = ra * 180 / math.Pi
return Limit360(ra)
}
func JupiterSeeDec(JD float64) float64 {
lo, bo := JupiterSeeLoBo(JD)
func JupiterApparentDec(JD float64) float64 {
lo, bo := JupiterApparentLoBo(JD)
sita := Sita(JD)
dec := ArcSin(Sin(bo)*Cos(sita) + Cos(bo)*Sin(sita)*Sin(lo))
return dec
}
func JupiterSeeRaDec(JD float64) (float64, float64) {
lo, bo := JupiterSeeLoBo(JD)
func JupiterApparentRaDec(JD float64) (float64, float64) {
lo, bo := JupiterApparentLoBo(JD)
sita := Sita(JD)
ra := math.Atan2((Sin(lo)*Cos(sita) - Tan(bo)*Sin(sita)), Cos(lo))
ra = ra * 180 / math.Pi
@@ -92,7 +92,7 @@ func EarthJupiterAway(JD float64) float64 {
return to
}
func JupiterSeeLo(JD float64) float64 {
func JupiterApparentLo(JD float64) float64 {
x, y, z := AJupiterXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = AJupiterXYZ(JD - to)
@@ -107,7 +107,7 @@ func JupiterSeeLo(JD float64) float64 {
return lo
}
func JupiterSeeBo(JD float64) float64 {
func JupiterApparentBo(JD float64) float64 {
x, y, z := AJupiterXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = AJupiterXYZ(JD - to)
@@ -121,7 +121,7 @@ func JupiterSeeBo(JD float64) float64 {
return bo
}
func JupiterSeeLoBo(JD float64) (float64, float64) {
func JupiterApparentLoBo(JD float64) (float64, float64) {
x, y, z := AJupiterXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = AJupiterXYZ(JD - to)
@@ -145,3 +145,293 @@ func JupiterMag(JD float64) float64 {
Mag := -9.40 + 5*math.Log10(AwaySun*AwayEarth) + 0.0005*i
return FloatRound(Mag, 2)
}
func JupiterHeight(jde, lon, lat, timezone float64) float64 {
// 转换为世界时
utcJde := jde - timezone/24.0
// 计算视恒星时
ra, dec := JupiterApparentRaDec(TD2UT(utcJde, true))
st := Limit360(ApparentSiderealTime(utcJde)*15 + lon)
// 计算时角
H := Limit360(st - ra)
// 高度角、时角与天球座标三角转换公式
// sin(h)=sin(lat)*sin(dec)+cos(dec)*cos(lat)*cos(H)
sinHeight := Sin(lat)*Sin(dec) + Cos(dec)*Cos(lat)*Cos(H)
return ArcSin(sinHeight)
}
func JupiterAzimuth(jde, lon, lat, timezone float64) float64 {
// 转换为世界时
utcJde := jde - timezone/24.0
// 计算视恒星时
ra, dec := JupiterApparentRaDec(TD2UT(utcJde, true))
st := Limit360(ApparentSiderealTime(utcJde)*15 + lon)
// 计算时角
H := Limit360(st - ra)
// 三角转换公式
tanAzimuth := Sin(H) / (Cos(H)*Sin(lat) - Tan(dec)*Cos(lat))
Azimuth := ArcTan(tanAzimuth)
if Azimuth < 0 {
if H/15 < 12 {
return Azimuth + 360
}
return Azimuth + 180
}
if H/15 < 12 {
return Azimuth + 180
}
return Azimuth
}
func JupiterHourAngle(JD, Lon, TZ float64) float64 {
startime := Limit360(ApparentSiderealTime(JD-TZ/24)*15 + Lon)
timeangle := startime - JupiterApparentRa(TD2UT(JD-TZ/24.0, true))
if timeangle < 0 {
timeangle += 360
}
return timeangle
}
func JupiterCulminationTime(jde, lon, timezone float64) float64 {
//jde 世界时，非力学时，当地时区 0时，无需转换力学时
//ra,dec 瞬时天球座标，非J2000等时间天球坐标
jde = math.Floor(jde) + 0.5
JD1 := jde + Limit360(360-JupiterHourAngle(jde, lon, timezone))/15.0/24.0*0.99726851851851851851
limitHA := func(jde, lon, timezone float64) float64 {
ha := JupiterHourAngle(jde, lon, timezone)
if ha < 180 {
ha += 360
}
return ha
}
for {
JD0 := JD1
stDegree := limitHA(JD0, lon, timezone) - 360
stDegreep := (limitHA(JD0+0.000005, lon, timezone) - limitHA(JD0-0.000005, lon, timezone)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return JD1
}
func JupiterRiseTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 {
return jupiterRiseDown(JD, Lon, Lat, TZ, ZS, HEI, true)
}
func JupiterDownTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 {
return jupiterRiseDown(JD, Lon, Lat, TZ, ZS, HEI, false)
}
func jupiterRiseDown(JD, Lon, Lat, TZ, ZS, HEI float64, isRise bool) float64 {
var An float64
JD = math.Floor(JD) + 0.5
ntz := math.Round(Lon / 15)
if ZS != 0 {
An = -0.8333
}
An = An - HeightDegreeByLat(HEI, Lat)
tztime := JupiterCulminationTime(JD, Lon, ntz)
if JupiterHeight(tztime, Lon, Lat, ntz) < An {
return -2 //极夜
}
if JupiterHeight(tztime-0.5, Lon, Lat, ntz) > An {
return -1 //极昼
}
dec := HSunApparentDec(TD2UT(tztime-ntz/24, true))
//(sin(ho)-sin(φ)*sin(δ2))/(cos(φ)*cos(δ2))
tmp := (Sin(An) - Sin(dec)*Sin(Lat)) / (Cos(dec) * Cos(Lat))
var rise float64
if math.Abs(tmp) <= 1 {
rzsc := ArcCos(tmp) / 15
if isRise {
rise = tztime - rzsc/24 - 25.0/24.0/60.0
} else {
rise = tztime + rzsc/24 - 25.0/24.0/60.0
}
} else {
rise = tztime
i := 0
//TODO:使用二分法计算
for JupiterHeight(rise, Lon, Lat, ntz) > An {
i++
if isRise {
rise -= 15.0 / 60.0 / 24.0
} else {
rise += 15.0 / 60.0 / 24.0
}
if i > 48 {
break
}
}
}
JD1 := rise
for {
JD0 := JD1
stDegree := JupiterHeight(JD0, Lon, Lat, ntz) - An
stDegreep := (JupiterHeight(JD0+0.000005, Lon, Lat, ntz) - JupiterHeight(JD0-0.000005, Lon, Lat, ntz)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return JD1 - ntz/24 + TZ/24
}
// Pos
const JUPITER_S_PERIOD = 1 / ((1 / 365.256363004) - (1 / 4332.59))
func jupiterConjunction(jde, degree float64, next uint8) float64 {
//0=last 1=next
decSub := func(jde float64, degree float64, filter bool) float64 {
sub := Limit360(Limit360(JupiterApparentLo(jde)-HSunApparentLo(jde)) - degree)
if filter {
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
}
return sub
}
dayCost := JUPITER_S_PERIOD / 360
nowSub := decSub(jde, degree, false)
if next == 0 {
jde -= (360 - nowSub) * dayCost
} else {
jde += dayCost * nowSub
}
JD1 := jde
for {
JD0 := JD1
stDegree := decSub(JD0, degree, true)
stDegreep := (decSub(JD0+0.000005, degree, true) - decSub(JD0-0.000005, degree, true)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return TD2UT(JD1, false)
}
func LastJupiterConjunction(jde float64) float64 {
return jupiterConjunction(jde, 0, 0)
}
func NextJupiterConjunction(jde float64) float64 {
return jupiterConjunction(jde, 0, 1)
}
func LastJupiterOpposition(jde float64) float64 {
return jupiterConjunction(jde, 180, 0)
}
func NextJupiterOpposition(jde float64) float64 {
return jupiterConjunction(jde, 180, 1)
}
func NextJupiterEasternQuadrature(jde float64) float64 {
return jupiterConjunction(jde, 90, 1)
}
func LastJupiterEasternQuadrature(jde float64) float64 {
return jupiterConjunction(jde, 90, 0)
}
func NextJupiterWesternQuadrature(jde float64) float64 {
return jupiterConjunction(jde, 270, 1)
}
func LastJupiterWesternQuadrature(jde float64) float64 {
return jupiterConjunction(jde, 270, 0)
}
func jupiterRetrograde(jde float64, isLeft bool) float64 {
//0=last 1=next
decSub := func(jde float64, val float64) float64 {
sub := JupiterApparentRa(jde+val) - JupiterApparentRa(jde-val)
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
return sub / (2 * val)
}
jde = NextJupiterOpposition(jde)
if isLeft {
jde -= 60
} else {
jde += 60
}
for {
nowSub := decSub(jde, 1.0/86400.0)
if math.Abs(nowSub) > 0.55 {
jde += 2
continue
}
break
}
JD1 := jde
for {
JD0 := JD1
stDegree := decSub(JD0, 2.0/86400.0)
stDegreep := (decSub(JD0+15.0/86400.0, 2.0/86400.0) - decSub(JD0-15.0/86400.0, 2.0/86400.0)) / (30.0 / 86400.0)
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 30.0/86400.0 {
break
}
}
JD1 = JD1 - 15.0/86400.0
min := JD1
minRa := 100.0
for i := 0.0; i < 60.0; i++ {
tmp := decSub(JD1+i*0.5/86400.0, 0.5/86400.0)
if math.Abs(tmp) < math.Abs(minRa) {
minRa = tmp
min = JD1 + i*0.5/86400.0
}
}
return TD2UT(min, false)
}
func NextJupiterRetrogradeToPrograde(jde float64) float64 {
date := jupiterRetrograde(jde, false)
if date < jde {
op := NextJupiterOpposition(jde)
return jupiterRetrograde(op+10, false)
}
return date
}
func LastJupiterRetrogradeToPrograde(jde float64) float64 {
jde = LastJupiterOpposition(jde) - 10
date := jupiterRetrograde(jde, false)
if date > jde {
op := LastJupiterOpposition(jde)
return jupiterRetrograde(op-10, false)
}
return date
}
func NextJupiterProgradeToRetrograde(jde float64) float64 {
date := jupiterRetrograde(jde, true)
if date < jde {
op := NextJupiterOpposition(jde)
return jupiterRetrograde(op+10, true)
}
return date
}
func LastJupiterProgradeToRetrograde(jde float64) float64 {
jde = LastJupiterOpposition(jde) - 10
date := jupiterRetrograde(jde, true)
if date > jde {
op := LastJupiterOpposition(jde)
return jupiterRetrograde(op-10, true)
}
return date
}
basic/jupiter_test.go
+13
View File
@@ -0,0 +1,13 @@
package basic
import (
"fmt"
"testing"
)
func TestJupiter(t *testing.T) {
jde := GetNowJDE() - 6000
for i := 0.00; i < 20; i++ {
fmt.Println(jde+i*365, JDE2Date(jde+i*365), JDE2Date(NextJupiterRetrogradeToPrograde(jde+i*365)))
}
}
basic/mars.go
+325 -16
View File
@@ -63,22 +63,22 @@ func AMarsXYZ(JD float64) (float64, float64, float64) {
return x, y, z
}
func MarsSeeRa(JD float64) float64 {
lo, bo := MarsSeeLoBo(JD)
func MarsApparentRa(JD float64) float64 {
lo, bo := MarsApparentLoBo(JD)
sita := Sita(JD)
ra := math.Atan2((Sin(lo)*Cos(sita) - Tan(bo)*Sin(sita)), Cos(lo))
ra = ra * 180 / math.Pi
return Limit360(ra)
}
func MarsSeeDec(JD float64) float64 {
lo, bo := MarsSeeLoBo(JD)
func MarsApparentDec(JD float64) float64 {
lo, bo := MarsApparentLoBo(JD)
sita := Sita(JD)
dec := ArcSin(Sin(bo)*Cos(sita) + Cos(bo)*Sin(sita)*Sin(lo))
return dec
}
func MarsSeeRaDec(JD float64) (float64, float64) {
lo, bo := MarsSeeLoBo(JD)
func MarsApparentRaDec(JD float64) (float64, float64) {
lo, bo := MarsApparentLoBo(JD)
sita := Sita(JD)
ra := math.Atan2((Sin(lo)*Cos(sita) - Tan(bo)*Sin(sita)), Cos(lo))
ra = ra * 180 / math.Pi
@@ -92,22 +92,21 @@ func EarthMarsAway(JD float64) float64 {
return to
}
func MarsSeeLo(JD float64) float64 {
func MarsApparentLo(JD float64) float64 {
x, y, z := AMarsXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = AMarsXYZ(JD - to)
lo := math.Atan2(y, x)
bo := math.Atan2(z, math.Sqrt(x*x+y*y))
lo = lo * 180 / math.Pi
bo = bo * 180 / math.Pi
lo = Limit360(lo)
//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-=GXCLo(lo,bo,JD)/3600;
//bo+=GXCBo(lo,bo,JD);
lo += HJZD(JD)
return lo
}
func MarsSeeBo(JD float64) float64 {
func MarsApparentBo(JD float64) float64 {
x, y, z := AMarsXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = AMarsXYZ(JD - to)
@@ -121,7 +120,7 @@ func MarsSeeBo(JD float64) float64 {
return bo
}
func MarsSeeLoBo(JD float64) (float64, float64) {
func MarsApparentLoBo(JD float64) (float64, float64) {
x, y, z := AMarsXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = AMarsXYZ(JD - to)
@@ -130,12 +129,32 @@ func MarsSeeLoBo(JD float64) (float64, float64) {
lo = lo * 180 / math.Pi
bo = bo * 180 / math.Pi
lo = Limit360(lo)
//lo-=GXCLo(lo,bo,JD)/3600;
//bo+=GXCBo(lo,bo,JD);
//lo -= GXCLo(lo, bo, JD) / 3600
//bo += GXCBo(lo, bo, JD)
lo += HJZD(JD)
return lo, bo
}
func MarsTrueLoBo(JD float64) (float64, float64) {
x, y, z := AMarsXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = AMarsXYZ(JD - to)
lo := math.Atan2(y, x)
bo := math.Atan2(z, math.Sqrt(x*x+y*y))
lo = lo * 180 / math.Pi
bo = bo * 180 / math.Pi
lo = Limit360(lo)
return lo, bo
}
func MarsTrueLo(JD float64) float64 {
x, y, _ := AMarsXYZ(JD)
lo := math.Atan2(y, x)
lo = lo * 180 / math.Pi
lo = Limit360(lo)
return lo
}
func MarsMag(JD float64) float64 {
AwaySun := MarsR(JD)
AwayEarth := EarthMarsAway(JD)
@@ -145,3 +164,293 @@ func MarsMag(JD float64) float64 {
Mag := -1.52 + 5*math.Log10(AwaySun*AwayEarth) + 0.016*i
return FloatRound(Mag, 2)
}
func MarsHeight(jde, lon, lat, timezone float64) float64 {
// 转换为世界时
utcJde := jde - timezone/24.0
// 计算视恒星时
ra, dec := MarsApparentRaDec(TD2UT(utcJde, true))
st := Limit360(ApparentSiderealTime(utcJde)*15 + lon)
// 计算时角
H := Limit360(st - ra)
// 高度角、时角与天球座标三角转换公式
// sin(h)=sin(lat)*sin(dec)+cos(dec)*cos(lat)*cos(H)
sinHeight := Sin(lat)*Sin(dec) + Cos(dec)*Cos(lat)*Cos(H)
return ArcSin(sinHeight)
}
func MarsAzimuth(jde, lon, lat, timezone float64) float64 {
// 转换为世界时
utcJde := jde - timezone/24.0
// 计算视恒星时
ra, dec := MarsApparentRaDec(TD2UT(utcJde, true))
st := Limit360(ApparentSiderealTime(utcJde)*15 + lon)
// 计算时角
H := Limit360(st - ra)
// 三角转换公式
tanAzimuth := Sin(H) / (Cos(H)*Sin(lat) - Tan(dec)*Cos(lat))
Azimuth := ArcTan(tanAzimuth)
if Azimuth < 0 {
if H/15 < 12 {
return Azimuth + 360
}
return Azimuth + 180
}
if H/15 < 12 {
return Azimuth + 180
}
return Azimuth
}
func MarsHourAngle(JD, Lon, TZ float64) float64 {
startime := Limit360(ApparentSiderealTime(JD-TZ/24)*15 + Lon)
timeangle := startime - MarsApparentRa(TD2UT(JD-TZ/24.0, true))
if timeangle < 0 {
timeangle += 360
}
return timeangle
}
func MarsCulminationTime(jde, lon, timezone float64) float64 {
//jde 世界时，非力学时，当地时区 0时，无需转换力学时
//ra,dec 瞬时天球座标，非J2000等时间天球坐标
jde = math.Floor(jde) + 0.5
JD1 := jde + Limit360(360-MarsHourAngle(jde, lon, timezone))/15.0/24.0*0.99726851851851851851
limitHA := func(jde, lon, timezone float64) float64 {
ha := MarsHourAngle(jde, lon, timezone)
if ha < 180 {
ha += 360
}
return ha
}
for {
JD0 := JD1
stDegree := limitHA(JD0, lon, timezone) - 360
stDegreep := (limitHA(JD0+0.000005, lon, timezone) - limitHA(JD0-0.000005, lon, timezone)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return JD1
}
func MarsRiseTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 {
return marsRiseDown(JD, Lon, Lat, TZ, ZS, HEI, true)
}
func MarsDownTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 {
return marsRiseDown(JD, Lon, Lat, TZ, ZS, HEI, false)
}
func marsRiseDown(JD, Lon, Lat, TZ, ZS, HEI float64, isRise bool) float64 {
var An float64
JD = math.Floor(JD) + 0.5
ntz := math.Round(Lon / 15)
if ZS != 0 {
An = -0.8333
}
An = An - HeightDegreeByLat(HEI, Lat)
tztime := MarsCulminationTime(JD, Lon, ntz)
if MarsHeight(tztime, Lon, Lat, ntz) < An {
return -2 //极夜
}
if MarsHeight(tztime-0.5, Lon, Lat, ntz) > An {
return -1 //极昼
}
dec := HSunApparentDec(TD2UT(tztime-ntz/24, true))
//(sin(ho)-sin(φ)*sin(δ2))/(cos(φ)*cos(δ2))
tmp := (Sin(An) - Sin(dec)*Sin(Lat)) / (Cos(dec) * Cos(Lat))
var rise float64
if math.Abs(tmp) <= 1 {
rzsc := ArcCos(tmp) / 15
if isRise {
rise = tztime - rzsc/24 - 25.0/24.0/60.0
} else {
rise = tztime + rzsc/24 - 25.0/24.0/60.0
}
} else {
rise = tztime
i := 0
//TODO:使用二分法计算
for MarsHeight(rise, Lon, Lat, ntz) > An {
i++
if isRise {
rise -= 15.0 / 60.0 / 24.0
} else {
rise += 15.0 / 60.0 / 24.0
}
if i > 48 {
break
}
}
}
JD1 := rise
for {
JD0 := JD1
stDegree := MarsHeight(JD0, Lon, Lat, ntz) - An
stDegreep := (MarsHeight(JD0+0.000005, Lon, Lat, ntz) - MarsHeight(JD0-0.000005, Lon, Lat, ntz)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return JD1 - ntz/24 + TZ/24
}
// Pos
const MARS_S_PERIOD = 1 / ((1 / 365.256363004) - (1 / 686.98))
func marsConjunction(jde, degree float64, next uint8) float64 {
//0=last 1=next
decSub := func(jde float64, degree float64, filter bool) float64 {
sub := Limit360(Limit360(MarsApparentLo(jde)-HSunApparentLo(jde)) - degree)
if filter {
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
}
return sub
}
dayCost := MARS_S_PERIOD / 360
nowSub := decSub(jde, degree, false)
if next == 0 {
jde -= (360 - nowSub) * dayCost
} else {
jde += dayCost * nowSub
}
JD1 := jde
for {
JD0 := JD1
stDegree := decSub(JD0, degree, true)
stDegreep := (decSub(JD0+0.000005, degree, true) - decSub(JD0-0.000005, degree, true)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return TD2UT(JD1, false)
}
func LastMarsConjunction(jde float64) float64 {
return marsConjunction(jde, 0, 0)
}
func NextMarsConjunction(jde float64) float64 {
return marsConjunction(jde, 0, 1)
}
func LastMarsOpposition(jde float64) float64 {
return marsConjunction(jde, 180, 0)
}
func NextMarsOpposition(jde float64) float64 {
return marsConjunction(jde, 180, 1)
}
func NextMarsEasternQuadrature(jde float64) float64 {
return marsConjunction(jde, 90, 1)
}
func LastMarsEasternQuadrature(jde float64) float64 {
return marsConjunction(jde, 90, 0)
}
func NextMarsWesternQuadrature(jde float64) float64 {
return marsConjunction(jde, 270, 1)
}
func LastMarsWesternQuadrature(jde float64) float64 {
return marsConjunction(jde, 270, 0)
}
func marsRetrograde(jde float64, isLeft bool) float64 {
//0=last 1=next
decSub := func(jde float64, val float64) float64 {
sub := MarsApparentRa(jde+val) - MarsApparentRa(jde-val)
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
return sub / (2 * val)
}
jde = NextMarsOpposition(jde)
if isLeft {
jde -= 60
} else {
jde += 60
}
for {
nowSub := decSub(jde, 1.0/86400.0)
if math.Abs(nowSub) > 0.55 {
jde += 2
continue
}
break
}
JD1 := jde
for {
JD0 := JD1
stDegree := decSub(JD0, 2.0/86400.0)
stDegreep := (decSub(JD0+15.0/86400.0, 2.0/86400.0) - decSub(JD0-15.0/86400.0, 2.0/86400.0)) / (30.0 / 86400.0)
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 30.0/86400.0 {
break
}
}
JD1 = JD1 - 15.0/86400.0
min := JD1
minRa := 100.0
for i := 0.0; i < 60.0; i++ {
tmp := decSub(JD1+i*0.5/86400.0, 0.5/86400.0)
if math.Abs(tmp) < math.Abs(minRa) {
minRa = tmp
min = JD1 + i*0.5/86400.0
}
}
return TD2UT(min, false)
}
func NextMarsRetrogradeToPrograde(jde float64) float64 {
date := marsRetrograde(jde, false)
if date < jde {
op := NextMarsOpposition(jde)
return marsRetrograde(op+10, false)
}
return date
}
func LastMarsRetrogradeToPrograde(jde float64) float64 {
jde = LastMarsOpposition(jde) - 10
date := marsRetrograde(jde, false)
if date > jde {
op := LastMarsOpposition(jde)
return marsRetrograde(op-10, false)
}
return date
}
func NextMarsProgradeToRetrograde(jde float64) float64 {
date := marsRetrograde(jde, true)
if date < jde {
op := NextMarsOpposition(jde)
return marsRetrograde(op+10, true)
}
return date
}
func LastMarsProgradeToRetrograde(jde float64) float64 {
jde = LastMarsOpposition(jde) - 10
date := marsRetrograde(jde, true)
if date > jde {
op := LastMarsOpposition(jde)
return marsRetrograde(op-10, true)
}
return date
}
basic/mars_test.go
+21
View File
@@ -0,0 +1,21 @@
package basic
import (
"fmt"
"testing"
)
func TestMars(t *testing.T) {
jde := GetNowJDE() - 6000
/*
fmt.Println(JDE2Date(VenusCulminationTime(jde, 115, 8)))
fmt.Println(JDE2Date(VenusRiseTime(jde, 115, 23, 8, 0, 0)))
fmt.Println(JDE2Date(VenusDownTime(jde, 115, 23, 8, 0, 0)))
fmt.Println("----------------")
*/
//LastVenusConjunction(2.4596600340162036e+06)
//fmt.Println(2.4590359532407406e+06, JDE2Date(2.4590359532407406e+06), JDE2Date(NextVenusRetrograde(2.4590359532407406e+06)))
for i := 0.00; i < 1; i++ {
fmt.Println(jde+i*740, JDE2Date(jde+i*740), JDE2Date(LastMarsProgradeToRetrograde(jde+i*740)))
}
}
basic/mercury.go
+481 -20
View File
@@ -63,27 +63,20 @@ func AMercuryXYZ(JD float64) (float64, float64, float64) {
return x, y, z
}
func MercurySeeRa(JD float64) float64 {
lo, bo := MercurySeeLoBo(JD)
sita := Sita(JD)
ra := math.Atan2((Sin(lo)*Cos(sita) - Tan(bo)*Sin(sita)), Cos(lo))
ra = ra * 180 / math.Pi
return Limit360(ra)
func MercuryApparentRa(JD float64) float64 {
lo, bo := MercuryApparentLoBo(JD)
return LoToRa(JD, lo, bo)
}
func MercurySeeDec(JD float64) float64 {
lo, bo := MercurySeeLoBo(JD)
func MercuryApparentDec(JD float64) float64 {
lo, bo := MercuryApparentLoBo(JD)
sita := Sita(JD)
dec := ArcSin(Sin(bo)*Cos(sita) + Cos(bo)*Sin(sita)*Sin(lo))
return dec
}
func MercurySeeRaDec(JD float64) (float64, float64) {
lo, bo := MercurySeeLoBo(JD)
sita := Sita(JD)
ra := math.Atan2((Sin(lo)*Cos(sita) - Tan(bo)*Sin(sita)), Cos(lo))
ra = ra * 180 / math.Pi
dec := ArcSin(Sin(bo)*Cos(sita) + Cos(bo)*Sin(sita)*Sin(lo))
return Limit360(ra), dec
func MercuryApparentRaDec(JD float64) (float64, float64) {
lo, bo := MercuryApparentLoBo(JD)
return LoBoToRaDec(JD, lo, bo)
}
func EarthMercuryAway(JD float64) float64 {
@@ -92,7 +85,7 @@ func EarthMercuryAway(JD float64) float64 {
return to
}
func MercurySeeLo(JD float64) float64 {
func MercuryApparentLo(JD float64) float64 {
x, y, z := AMercuryXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = AMercuryXYZ(JD - to)
@@ -107,7 +100,7 @@ func MercurySeeLo(JD float64) float64 {
return lo
}
func MercurySeeBo(JD float64) float64 {
func MercuryApparentBo(JD float64) float64 {
x, y, z := AMercuryXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = AMercuryXYZ(JD - to)
@@ -121,7 +114,7 @@ func MercurySeeBo(JD float64) float64 {
return bo
}
func MercurySeeLoBo(JD float64) (float64, float64) {
func MercuryApparentLoBo(JD float64) (float64, float64) {
x, y, z := AMercuryXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = AMercuryXYZ(JD - to)
@@ -129,10 +122,9 @@ func MercurySeeLoBo(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)
lo = Limit360(lo) + HJZD(JD)
//lo-=GXCLo(lo,bo,JD)/3600;
//bo+=GXCBo(lo,bo,JD);
lo += HJZD(JD)
return lo, bo
}
@@ -145,3 +137,472 @@ func MercuryMag(JD float64) float64 {
Mag := -0.42 + 5*math.Log10(AwaySun*AwayEarth) + 0.0380*i - 0.000273*i*i + 0.000002*i*i*i
return FloatRound(Mag, 2)
}
func MercuryHeight(jde, lon, lat, timezone float64) float64 {
// 转换为世界时
utcJde := jde - timezone/24.0
// 计算视恒星时
ra, dec := MercuryApparentRaDec(TD2UT(utcJde, true))
st := Limit360(ApparentSiderealTime(utcJde)*15 + lon)
// 计算时角
H := Limit360(st - ra)
// 高度角、时角与天球座标三角转换公式
// sin(h)=sin(lat)*sin(dec)+cos(dec)*cos(lat)*cos(H)
sinHeight := Sin(lat)*Sin(dec) + Cos(dec)*Cos(lat)*Cos(H)
return ArcSin(sinHeight)
}
func MercuryAzimuth(jde, lon, lat, timezone float64) float64 {
// 转换为世界时
utcJde := jde - timezone/24.0
// 计算视恒星时
ra, dec := MercuryApparentRaDec(TD2UT(utcJde, true))
st := Limit360(ApparentSiderealTime(utcJde)*15 + lon)
// 计算时角
H := Limit360(st - ra)
// 三角转换公式
tanAzimuth := Sin(H) / (Cos(H)*Sin(lat) - Tan(dec)*Cos(lat))
Azimuth := ArcTan(tanAzimuth)
if Azimuth < 0 {
if H/15 < 12 {
return Azimuth + 360
}
return Azimuth + 180
}
if H/15 < 12 {
return Azimuth + 180
}
return Azimuth
}
func MercuryHourAngle(JD, Lon, TZ float64) float64 {
startime := Limit360(ApparentSiderealTime(JD-TZ/24)*15 + Lon)
timeangle := startime - MercuryApparentRa(TD2UT(JD-TZ/24.0, true))
if timeangle < 0 {
timeangle += 360
}
return timeangle
}
func MercuryCulminationTime(jde, lon, timezone float64) float64 {
//jde 世界时，非力学时，当地时区 0时，无需转换力学时
//ra,dec 瞬时天球座标，非J2000等时间天球坐标
jde = math.Floor(jde) + 0.5
JD1 := jde + Limit360(360-MercuryHourAngle(jde, lon, timezone))/15.0/24.0*0.99726851851851851851
limitHA := func(jde, lon, timezone float64) float64 {
ha := MercuryHourAngle(jde, lon, timezone)
if ha < 180 {
ha += 360
}
return ha
}
for {
JD0 := JD1
stDegree := limitHA(JD0, lon, timezone) - 360
stDegreep := (limitHA(JD0+0.000005, lon, timezone) - limitHA(JD0-0.000005, lon, timezone)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return JD1
}
func MercuryRiseTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 {
return mercuryRiseDown(JD, Lon, Lat, TZ, ZS, HEI, true)
}
func MercuryDownTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 {
return mercuryRiseDown(JD, Lon, Lat, TZ, ZS, HEI, false)
}
func mercuryRiseDown(JD, Lon, Lat, TZ, ZS, HEI float64, isRise bool) float64 {
var An float64
JD = math.Floor(JD) + 0.5
ntz := math.Round(Lon / 15)
if ZS != 0 {
An = -0.8333
}
An = An - HeightDegreeByLat(HEI, Lat)
tztime := MercuryCulminationTime(JD, Lon, ntz)
if MercuryHeight(tztime, Lon, Lat, ntz) < An {
return -2 //极夜
}
if MercuryHeight(tztime-0.5, Lon, Lat, ntz) > An {
return -1 //极昼
}
dec := HSunApparentDec(TD2UT(tztime-ntz/24, true))
//(sin(ho)-sin(φ)*sin(δ2))/(cos(φ)*cos(δ2))
tmp := (Sin(An) - Sin(dec)*Sin(Lat)) / (Cos(dec) * Cos(Lat))
var rise float64
if math.Abs(tmp) <= 1 {
rzsc := ArcCos(tmp) / 15
if isRise {
rise = tztime - rzsc/24 - 25.0/24.0/60.0
} else {
rise = tztime + rzsc/24 - 25.0/24.0/60.0
}
} else {
rise = tztime
i := 0
//TODO:使用二分法计算
for MercuryHeight(rise, Lon, Lat, ntz) > An {
i++
if isRise {
rise -= 15.0 / 60.0 / 24.0
} else {
rise += 15.0 / 60.0 / 24.0
}
if i > 48 {
break
}
}
}
JD1 := rise
for {
JD0 := JD1
stDegree := MercuryHeight(JD0, Lon, Lat, ntz) - An
stDegreep := (MercuryHeight(JD0+0.000005, Lon, Lat, ntz) - MercuryHeight(JD0-0.000005, Lon, Lat, ntz)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return JD1 - ntz/24 + TZ/24
}
// Pos
const MERCURY_S_PERIOD = 1 / ((1 / 87.9691) - (1 / 365.256363004))
func mercuryConjunction(jde float64, next uint8) float64 {
//0=last 1=next
decSub := func(jde float64) float64 {
sub := Limit360(MercuryApparentLo(jde) - HSunApparentLo(jde))
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
return sub
}
nowSub := decSub(jde)
// pos 大于0:远离太阳 小于0:靠近太阳
pos := math.Abs(decSub(jde+1/86400.0)) - math.Abs(nowSub)
if pos >= 0 && next == 1 && nowSub > 0 {
jde += MERCURY_S_PERIOD/8.0 + 2
}
if pos >= 0 && next == 1 && nowSub < 0 {
jde += MERCURY_S_PERIOD/6.0 + 2
}
if pos <= 0 && next == 0 && nowSub < 0 {
jde -= MERCURY_S_PERIOD/8.0 + 2
}
if pos <= 0 && next == 0 && nowSub > 0 {
jde -= MERCURY_S_PERIOD/6.0 + 2
}
for {
nowSub := decSub(jde)
pos := math.Abs(decSub(jde+1/86400.0)) - math.Abs(nowSub)
if math.Abs(nowSub) > 12 || (pos > 0 && next == 1) || (pos < 0 && next == 0) {
if next == 1 {
jde += 2
} else {
jde -= 2
}
continue
}
break
}
JD1 := jde
for {
JD0 := JD1
stDegree := decSub(JD0)
stDegreep := (decSub(JD0+0.000005) - decSub(JD0-0.000005)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return TD2UT(JD1, false)
}
func LastMercuryConjunction(jde float64) float64 {
return mercuryConjunction(jde, 0)
}
func NextMercuryConjunction(jde float64) float64 {
return mercuryConjunction(jde, 1)
}
func NextMercuryInferiorConjunction(jde float64) float64 {
date := NextMercuryConjunction(jde)
if EarthMercuryAway(date) > EarthAway(date) {
return NextMercuryConjunction(date + 2)
}
return date
}
func NextMercurySuperiorConjunction(jde float64) float64 {
date := NextMercuryConjunction(jde)
if EarthMercuryAway(date) < EarthAway(date) {
return NextMercuryConjunction(date + 2)
}
return date
}
func LastMercuryInferiorConjunction(jde float64) float64 {
date := LastMercuryConjunction(jde)
if EarthMercuryAway(date) > EarthAway(date) {
return LastMercuryConjunction(date - 2)
}
return date
}
func LastMercurySuperiorConjunction(jde float64) float64 {
date := LastMercuryConjunction(jde)
if EarthMercuryAway(date) < EarthAway(date) {
return LastMercuryConjunction(date - 2)
}
return date
}
func mercuryRetrograde(jde float64) float64 {
//0=last 1=next
decSunSub := func(jde float64) float64 {
sub := Limit360(MercuryApparentRa(jde) - SunApparentRa(jde))
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
return sub
}
decSub := func(jde float64, val float64) float64 {
sub := MercuryApparentRa(jde+val) - MercuryApparentRa(jde-val)
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
return sub / (2 * val)
}
lastHe := LastMercuryConjunction(jde)
nextHe := NextMercuryConjunction(jde)
nowSub := decSunSub(jde)
if nowSub > 0 {
jde = lastHe + ((nextHe - lastHe) / 5.0 * 3.5)
} else {
jde = lastHe + ((nextHe - lastHe) / 5.5)
}
for {
nowSub := decSub(jde, 1.0/86400.0)
if math.Abs(nowSub) > 0.55 {
jde += 2
continue
}
break
}
JD1 := jde
for {
JD0 := JD1
stDegree := decSub(JD0, 2.0/86400.0)
stDegreep := (decSub(JD0+15.0/86400.0, 2.0/86400.0) - decSub(JD0-15.0/86400.0, 2.0/86400.0)) / (30.0 / 86400.0)
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 30.0/86400.0 {
break
}
}
JD1 = JD1 - 15.0/86400.0
min := JD1
minRa := 100.0
for i := 0.0; i < 60.0; i++ {
tmp := decSub(JD1+i*0.5/86400.0, 0.5/86400.0)
if math.Abs(tmp) < math.Abs(minRa) {
minRa = tmp
min = JD1 + i*0.5/86400.0
}
}
//fmt.Println((min - lastHe) / (nextHe - lastHe))
return TD2UT(min, false)
}
func NextMercuryRetrograde(jde float64) float64 {
date := mercuryRetrograde(jde)
if date < jde {
nextHe := NextMercuryConjunction(jde)
return mercuryRetrograde(nextHe + 2)
}
return date
}
func LastMercuryRetrograde(jde float64) float64 {
lastHe := LastMercuryConjunction(jde)
date := mercuryRetrograde(lastHe + 2)
if date > jde {
lastLastHe := LastMercuryConjunction(lastHe - 2)
return mercuryRetrograde(lastLastHe + 2)
}
return date
}
func NextMercuryProgradeToRetrograde(jde float64) float64 {
date := NextMercuryRetrograde(jde)
sub := Limit360(MercuryApparentRa(date) - SunApparentRa(date))
if sub > 180 {
return NextMercuryRetrograde(date + MERCURY_S_PERIOD/2)
}
return date
}
func NextMercuryRetrogradeToPrograde(jde float64) float64 {
date := NextMercuryRetrograde(jde)
sub := Limit360(MercuryApparentRa(date) - SunApparentRa(date))
if sub < 180 {
return NextMercuryRetrograde(date + 12)
}
return date
}
func LastMercuryProgradeToRetrograde(jde float64) float64 {
date := LastMercuryRetrograde(jde)
sub := Limit360(MercuryApparentRa(date) - SunApparentRa(date))
if sub > 180 {
return LastMercuryRetrograde(date - 12)
}
return date
}
func LastMercuryRetrogradeToPrograde(jde float64) float64 {
date := LastMercuryRetrograde(jde)
sub := Limit360(MercuryApparentRa(date) - SunApparentRa(date))
if sub < 180 {
return LastMercuryRetrograde(date - MERCURY_S_PERIOD/2)
}
return date
}
func MercurySunElongation(jde float64) float64 {
lo1, bo1 := MercuryApparentLoBo(jde)
lo2 := SunApparentLo(jde)
bo2 := HSunTrueBo(jde)
return StarAngularSeparation(lo1, bo1, lo2, bo2)
}
func mercuryGreatestElongation(jde float64) float64 {
decSunSub := func(jde float64) float64 {
sub := Limit360(MercuryApparentRa(jde) - SunApparentRa(jde))
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
return sub
}
decSub := func(jde float64, val float64) float64 {
sub := MercurySunElongation(jde+val) - MercurySunElongation(jde-val)
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
return sub / (2 * val)
}
lastHe := LastMercuryConjunction(jde)
nextHe := NextMercuryConjunction(jde)
nowSub := decSunSub(jde)
if nowSub > 0 {
jde = lastHe + ((nextHe - lastHe) / 5.0 * 2.0)
} else {
jde = lastHe + ((nextHe - lastHe) / 6.0)
}
for {
nowSub := decSub(jde, 1.0/86400.0)
if math.Abs(nowSub) > 0.4 {
jde += 2
continue
}
break
}
JD1 := jde
for {
JD0 := JD1
stDegree := decSub(JD0, 2.0/86400.0)
stDegreep := (decSub(JD0+15.0/86400.0, 2.0/86400.0) - decSub(JD0-15.0/86400.0, 2.0/86400.0)) / (30.0 / 86400.0)
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 30.0/86400.0 {
break
}
}
JD1 = JD1 - 15.0/86400.0
min := JD1
minRa := 100.0
for i := 0.0; i < 60.0; i++ {
tmp := decSub(JD1+i*0.5/86400.0, 0.5/86400.0)
if math.Abs(tmp) < math.Abs(minRa) {
minRa = tmp
min = JD1 + i*0.5/86400.0
}
}
//fmt.Println((min - lastHe) / (nextHe - lastHe))
return TD2UT(min, false)
}
func NextMercuryGreatestElongation(jde float64) float64 {
date := mercuryGreatestElongation(jde)
if date < jde {
nextHe := NextMercuryConjunction(jde)
return mercuryGreatestElongation(nextHe + 2)
}
return date
}
func LastMercuryGreatestElongation(jde float64) float64 {
lastHe := LastMercuryConjunction(jde)
date := mercuryGreatestElongation(lastHe + 2)
if date > jde {
lastLastHe := LastMercuryConjunction(lastHe - 2)
return mercuryGreatestElongation(lastLastHe + 2)
}
return date
}
func NextMercuryGreatestElongationEast(jde float64) float64 {
date := NextMercuryGreatestElongation(jde)
sub := Limit360(MercuryApparentRa(date) - SunApparentRa(date))
if sub > 180 {
return NextMercuryGreatestElongation(date + 1)
}
return date
}
func NextMercuryGreatestElongationWest(jde float64) float64 {
date := NextMercuryGreatestElongation(jde)
sub := Limit360(MercuryApparentRa(date) - SunApparentRa(date))
if sub < 180 {
return NextMercuryGreatestElongation(date + 1)
}
return date
}
func LastMercuryGreatestElongationEast(jde float64) float64 {
date := LastMercuryGreatestElongation(jde)
sub := Limit360(MercuryApparentRa(date) - SunApparentRa(date))
if sub > 180 {
return LastMercuryGreatestElongation(date - 1)
}
return date
}
func LastMercuryGreatestElongationWest(jde float64) float64 {
date := LastMercuryGreatestElongation(jde)
sub := Limit360(MercuryApparentRa(date) - SunApparentRa(date))
if sub < 180 {
return LastMercuryGreatestElongation(date - 1)
}
return date
}
basic/mercury_test.go
+16
View File
@@ -0,0 +1,16 @@
package basic
import (
"fmt"
"testing"
)
func TestMercury(t *testing.T) {
jde := GetNowJDE()
fmt.Println(2.459941309513889e+06, JDE2Date(2.459941309513889e+06), JDE2Date(8.0/24.0+LastMercuryGreatestElongation(2.459941309513889e+06)))
fmt.Println("-------------for------------")
for i := 0.00; i < 700.0; i += 5 {
fmt.Println(jde+i, JDE2Date(jde+i), JDE2Date(8.0/24.0+LastMercuryGreatestElongationWest(jde+i)))
// fmt.Println("")
}
}
basic/moon.go
+29 -39
View File
@@ -1052,18 +1052,7 @@ func MoonTrueDec(JD float64) float64 {
* 月球真赤经
*/
func MoonTrueRa(JD float64) float64 {
MoonLo := MoonApparentLo(JD)
MoonBo := MoonTrueBo(JD)
tmp := (Sin(MoonLo)*Cos(Sita(JD)) - Tan(MoonBo)*Sin(Sita(JD))) / Cos(MoonLo)
tmp = ArcTan(tmp)
if MoonLo >= 90 && MoonLo < 180 {
tmp = 180 + tmp
} else if MoonLo >= 180 && MoonLo < 270 {
tmp = 180 + tmp
} else if MoonLo >= 270 && MoonLo <= 360 {
tmp = 360 + tmp
}
return tmp
return LoToRa(JD, MoonApparentLo(JD), MoonTrueBo(JD))
}
/**
@@ -1382,6 +1371,7 @@ func HMoonHeight(JD, Lon, Lat, TZ float64) float64 {
return ArcSin(tmp2)
}
// 废弃
func GetMoonTZTime(JD, Lon, Lat, TZ float64) float64 { //实际中天时间{
JD = math.Floor(JD) + 0.5
ttm := MoonTimeAngle(JD, Lon, Lat, TZ)
@@ -1401,9 +1391,33 @@ func GetMoonTZTime(JD, Lon, Lat, TZ float64) float64 { //实际中天时间{
return JD1
}
func MoonCulminationTime(jde, lon, lat, timezone float64) float64 {
//jde 世界时，非力学时，当地时区 0时，无需转换力学时
//ra,dec 瞬时天球座标，非J2000等时间天球坐标
jde = math.Floor(jde) + 0.5
JD1 := jde + Limit360(360-MoonTimeAngle(jde, lon, lat, timezone))/15.0/24.0/0.9
limitHA := func(jde, lon, timezone float64) float64 {
ha := MoonTimeAngle(jde, lon, lat, timezone)
if ha < 180 {
ha += 360
}
return ha
}
for {
JD0 := JD1
stDegree := limitHA(JD0, lon, timezone) - 360
stDegreep := (limitHA(JD0+0.000005, lon, timezone) - limitHA(JD0-0.000005, lon, timezone)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return JD1
}
func MoonTimeAngle(JD, Lon, Lat, TZ float64) float64 {
startime := Limit360(ApparentSiderealTime(JD-TZ/24)*15 + Lon)
timeangle := startime - HMoonApparentRa(JD-TZ/24, Lon, Lat, TZ)
timeangle := startime - HMoonApparentRa(JD, Lon, Lat, TZ)
if timeangle < 0 {
timeangle += 360
}
@@ -1667,20 +1681,7 @@ func HMoonApparentLo(JD float64) float64 {
}
func HMoonTrueRaDec(JD float64) (float64, float64) {
MoonLo := HMoonApparentLo(JD)
MoonBo := HMoonTrueBo(JD)
tmp := Sin(MoonBo)*Cos(Sita(JD)) + Cos(MoonBo)*Sin(Sita(JD))*Sin(MoonLo)
res := ArcSin(tmp)
tmp = (Sin(MoonLo)*Cos(Sita(JD)) - Tan(MoonBo)*Sin(Sita(JD))) / Cos(MoonLo)
tmp = ArcTan(tmp)
if MoonLo >= 90 && MoonLo < 180 {
tmp = 180 + tmp
} else if MoonLo >= 180 && MoonLo < 270 {
tmp = 180 + tmp
} else if MoonLo >= 270 && MoonLo <= 360 {
tmp = 360 + tmp
}
return tmp, res
return LoBoToRaDec(JD, HMoonApparentLo(JD), HMoonTrueBo(JD))
}
/*
@@ -1698,18 +1699,7 @@ func HMoonTrueDec(JD float64) float64 {
* 月球真赤经
*/
func HMoonTrueRa(JD float64) float64 {
MoonLo := HMoonApparentLo(JD)
MoonBo := HMoonTrueBo(JD)
tmp := (Sin(MoonLo)*Cos(Sita(JD)) - Tan(MoonBo)*Sin(Sita(JD))) / Cos(MoonLo)
tmp = ArcTan(tmp)
if MoonLo >= 90 && MoonLo < 180 {
tmp = 180 + tmp
} else if MoonLo >= 180 && MoonLo < 270 {
tmp = 180 + tmp
} else if MoonLo >= 270 && MoonLo <= 360 {
tmp = 360 + tmp
}
return tmp
return LoToRa(JD, HMoonApparentLo(JD), HMoonTrueBo(JD))
}
/**
basic/moon_test.go
+14
View File
@@ -2,6 +2,7 @@ package basic
import (
"fmt"
"math"
"testing"
"time"
)
@@ -45,3 +46,16 @@ func Test_MoonS(t *testing.T) {
fmt.Println(time.Now().UnixNano() - a)
//fmt.Printf("%.14f", GetMoonRiseTime(2451547, 115, 32, 8, 0))
}
func TestMoonCu(t *testing.T) {
jde := math.Floor(GetNowJDE() - 20.0/24.0)
n := MoonCulminationTime(jde, 115, 23, 8)
fmt.Println(JDE2Date(n))
fmt.Println(MoonTimeAngle(n, 115, 23, 8))
fmt.Println(MoonAngle(n, 115, 23, 8))
//fmt.Println(JDE2Date(jde))
//ra, dec := HMoonApparentRaDec(jde, 115, 23, 8)
//fmt.Println(tools.Format(ra/15, 1), tools.Format(dec, 0))
//fmt.Println(JDE2Date(GetMoonTZTime(jde, 115, 23, 8)))
//fmt.Println(JDE2Date(GetMoonDownTime(jde+1, 115, 23, 8, 1, 0)))
}
basic/neptune.go
+299 -9
View File
@@ -63,22 +63,22 @@ func ANeptuneXYZ(JD float64) (float64, float64, float64) {
return x, y, z
}
func NeptuneSeeRa(JD float64) float64 {
lo, bo := NeptuneSeeLoBo(JD)
func NeptuneApparentRa(JD float64) float64 {
lo, bo := NeptuneApparentLoBo(JD)
sita := Sita(JD)
ra := math.Atan2((Sin(lo)*Cos(sita) - Tan(bo)*Sin(sita)), Cos(lo))
ra = ra * 180 / math.Pi
return Limit360(ra)
}
func NeptuneSeeDec(JD float64) float64 {
lo, bo := NeptuneSeeLoBo(JD)
func NeptuneApparentDec(JD float64) float64 {
lo, bo := NeptuneApparentLoBo(JD)
sita := Sita(JD)
dec := ArcSin(Sin(bo)*Cos(sita) + Cos(bo)*Sin(sita)*Sin(lo))
return dec
}
func NeptuneSeeRaDec(JD float64) (float64, float64) {
lo, bo := NeptuneSeeLoBo(JD)
func NeptuneApparentRaDec(JD float64) (float64, float64) {
lo, bo := NeptuneApparentLoBo(JD)
sita := Sita(JD)
ra := math.Atan2((Sin(lo)*Cos(sita) - Tan(bo)*Sin(sita)), Cos(lo))
ra = ra * 180 / math.Pi
@@ -92,7 +92,7 @@ func EarthNeptuneAway(JD float64) float64 {
return to
}
func NeptuneSeeLo(JD float64) float64 {
func NeptuneApparentLo(JD float64) float64 {
x, y, z := ANeptuneXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = ANeptuneXYZ(JD - to)
@@ -107,7 +107,7 @@ func NeptuneSeeLo(JD float64) float64 {
return lo
}
func NeptuneSeeBo(JD float64) float64 {
func NeptuneApparentBo(JD float64) float64 {
x, y, z := ANeptuneXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = ANeptuneXYZ(JD - to)
@@ -121,7 +121,7 @@ func NeptuneSeeBo(JD float64) float64 {
return bo
}
func NeptuneSeeLoBo(JD float64) (float64, float64) {
func NeptuneApparentLoBo(JD float64) (float64, float64) {
x, y, z := ANeptuneXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = ANeptuneXYZ(JD - to)
@@ -145,3 +145,293 @@ func NeptuneMag(JD float64) float64 {
Mag := -6.87 + 5*math.Log10(AwaySun*AwayEarth)
return FloatRound(Mag, 2)
}
func NeptuneHeight(jde, lon, lat, timezone float64) float64 {
// 转换为世界时
utcJde := jde - timezone/24.0
// 计算视恒星时
ra, dec := NeptuneApparentRaDec(TD2UT(utcJde, true))
st := Limit360(ApparentSiderealTime(utcJde)*15 + lon)
// 计算时角
H := Limit360(st - ra)
// 高度角、时角与天球座标三角转换公式
// sin(h)=sin(lat)*sin(dec)+cos(dec)*cos(lat)*cos(H)
sinHeight := Sin(lat)*Sin(dec) + Cos(dec)*Cos(lat)*Cos(H)
return ArcSin(sinHeight)
}
func NeptuneAzimuth(jde, lon, lat, timezone float64) float64 {
// 转换为世界时
utcJde := jde - timezone/24.0
// 计算视恒星时
ra, dec := NeptuneApparentRaDec(TD2UT(utcJde, true))
st := Limit360(ApparentSiderealTime(utcJde)*15 + lon)
// 计算时角
H := Limit360(st - ra)
// 三角转换公式
tanAzimuth := Sin(H) / (Cos(H)*Sin(lat) - Tan(dec)*Cos(lat))
Azimuth := ArcTan(tanAzimuth)
if Azimuth < 0 {
if H/15 < 12 {
return Azimuth + 360
}
return Azimuth + 180
}
if H/15 < 12 {
return Azimuth + 180
}
return Azimuth
}
func NeptuneHourAngle(JD, Lon, TZ float64) float64 {
startime := Limit360(ApparentSiderealTime(JD-TZ/24)*15 + Lon)
timeangle := startime - NeptuneApparentRa(TD2UT(JD-TZ/24.0, true))
if timeangle < 0 {
timeangle += 360
}
return timeangle
}
func NeptuneCulminationTime(jde, lon, timezone float64) float64 {
//jde 世界时，非力学时，当地时区 0时，无需转换力学时
//ra,dec 瞬时天球座标，非J2000等时间天球坐标
jde = math.Floor(jde) + 0.5
JD1 := jde + Limit360(360-NeptuneHourAngle(jde, lon, timezone))/15.0/24.0*0.99726851851851851851
limitHA := func(jde, lon, timezone float64) float64 {
ha := NeptuneHourAngle(jde, lon, timezone)
if ha < 180 {
ha += 360
}
return ha
}
for {
JD0 := JD1
stDegree := limitHA(JD0, lon, timezone) - 360
stDegreep := (limitHA(JD0+0.000005, lon, timezone) - limitHA(JD0-0.000005, lon, timezone)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return JD1
}
func NeptuneRiseTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 {
return neptuneRiseDown(JD, Lon, Lat, TZ, ZS, HEI, true)
}
func NeptuneDownTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 {
return neptuneRiseDown(JD, Lon, Lat, TZ, ZS, HEI, false)
}
func neptuneRiseDown(JD, Lon, Lat, TZ, ZS, HEI float64, isRise bool) float64 {
var An float64
JD = math.Floor(JD) + 0.5
ntz := math.Round(Lon / 15)
if ZS != 0 {
An = -0.8333
}
An = An - HeightDegreeByLat(HEI, Lat)
tztime := NeptuneCulminationTime(JD, Lon, ntz)
if NeptuneHeight(tztime, Lon, Lat, ntz) < An {
return -2 //极夜
}
if NeptuneHeight(tztime-0.5, Lon, Lat, ntz) > An {
return -1 //极昼
}
dec := HSunApparentDec(TD2UT(tztime-ntz/24, true))
//(sin(ho)-sin(φ)*sin(δ2))/(cos(φ)*cos(δ2))
tmp := (Sin(An) - Sin(dec)*Sin(Lat)) / (Cos(dec) * Cos(Lat))
var rise float64
if math.Abs(tmp) <= 1 {
rzsc := ArcCos(tmp) / 15
if isRise {
rise = tztime - rzsc/24 - 25.0/24.0/60.0
} else {
rise = tztime + rzsc/24 - 25.0/24.0/60.0
}
} else {
rise = tztime
i := 0
//TODO:使用二分法计算
for NeptuneHeight(rise, Lon, Lat, ntz) > An {
i++
if isRise {
rise -= 15.0 / 60.0 / 24.0
} else {
rise += 15.0 / 60.0 / 24.0
}
if i > 48 {
break
}
}
}
JD1 := rise
for {
JD0 := JD1
stDegree := NeptuneHeight(JD0, Lon, Lat, ntz) - An
stDegreep := (NeptuneHeight(JD0+0.000005, Lon, Lat, ntz) - NeptuneHeight(JD0-0.000005, Lon, Lat, ntz)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return JD1 - ntz/24 + TZ/24
}
// Pos
const NEPTUNE_S_PERIOD = 1 / ((1 / 365.256363004) - (1 / 4332.59))
func neptuneConjunction(jde, degree float64, next uint8) float64 {
//0=last 1=next
decSub := func(jde float64, degree float64, filter bool) float64 {
sub := Limit360(Limit360(NeptuneApparentLo(jde)-HSunApparentLo(jde)) - degree)
if filter {
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
}
return sub
}
dayCost := NEPTUNE_S_PERIOD / 360
nowSub := decSub(jde, degree, false)
if next == 0 {
jde -= (360 - nowSub) * dayCost
} else {
jde += dayCost * nowSub
}
JD1 := jde
for {
JD0 := JD1
stDegree := decSub(JD0, degree, true)
stDegreep := (decSub(JD0+0.000005, degree, true) - decSub(JD0-0.000005, degree, true)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return TD2UT(JD1, false)
}
func LastNeptuneConjunction(jde float64) float64 {
return neptuneConjunction(jde, 0, 0)
}
func NextNeptuneConjunction(jde float64) float64 {
return neptuneConjunction(jde, 0, 1)
}
func LastNeptuneOpposition(jde float64) float64 {
return neptuneConjunction(jde, 180, 0)
}
func NextNeptuneOpposition(jde float64) float64 {
return neptuneConjunction(jde, 180, 1)
}
func NextNeptuneEasternQuadrature(jde float64) float64 {
return neptuneConjunction(jde, 90, 1)
}
func LastNeptuneEasternQuadrature(jde float64) float64 {
return neptuneConjunction(jde, 90, 0)
}
func NextNeptuneWesternQuadrature(jde float64) float64 {
return neptuneConjunction(jde, 270, 1)
}
func LastNeptuneWesternQuadrature(jde float64) float64 {
return neptuneConjunction(jde, 270, 0)
}
func neptuneRetrograde(jde float64, isLeft bool) float64 {
//0=last 1=next
decSub := func(jde float64, val float64) float64 {
sub := NeptuneApparentRa(jde+val) - NeptuneApparentRa(jde-val)
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
return sub / (2 * val)
}
jde = NextNeptuneOpposition(jde)
if isLeft {
jde -= 60
} else {
jde += 60
}
for {
nowSub := decSub(jde, 1.0/86400.0)
if math.Abs(nowSub) > 0.55 {
jde += 2
continue
}
break
}
JD1 := jde
for {
JD0 := JD1
stDegree := decSub(JD0, 2.0/86400.0)
stDegreep := (decSub(JD0+15.0/86400.0, 2.0/86400.0) - decSub(JD0-15.0/86400.0, 2.0/86400.0)) / (30.0 / 86400.0)
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 30.0/86400.0 {
break
}
}
JD1 = JD1 - 15.0/86400.0
min := JD1
minRa := 100.0
for i := 0.0; i < 60.0; i++ {
tmp := decSub(JD1+i*0.5/86400.0, 0.5/86400.0)
if math.Abs(tmp) < math.Abs(minRa) {
minRa = tmp
min = JD1 + i*0.5/86400.0
}
}
return TD2UT(min, false)
}
func NextNeptuneRetrogradeToPrograde(jde float64) float64 {
date := neptuneRetrograde(jde, false)
if date < jde {
op := NextNeptuneOpposition(jde)
return neptuneRetrograde(op+10, false)
}
return date
}
func LastNeptuneRetrogradeToPrograde(jde float64) float64 {
jde = LastNeptuneOpposition(jde) - 10
date := neptuneRetrograde(jde, false)
if date > jde {
op := LastNeptuneOpposition(jde)
return neptuneRetrograde(op-10, false)
}
return date
}
func NextNeptuneProgradeToRetrograde(jde float64) float64 {
date := neptuneRetrograde(jde, true)
if date < jde {
op := NextNeptuneOpposition(jde)
return neptuneRetrograde(op+10, true)
}
return date
}
func LastNeptuneProgradeToRetrograde(jde float64) float64 {
jde = LastNeptuneOpposition(jde) - 10
date := neptuneRetrograde(jde, true)
if date > jde {
op := LastNeptuneOpposition(jde)
return neptuneRetrograde(op-10, true)
}
return date
}
basic/planet_test.go
+2 -2
View File
@@ -6,10 +6,10 @@ import (
)
func Test_Ra(t *testing.T) {
ra, dec := UranusSeeRaDec(2456789.12345)
ra, dec := UranusApparentRaDec(2456789.12345)
fmt.Printf("%.14f\n%.14f\n", ra, dec)
fmt.Println(UranusMag(2456789.12345))
ra, dec = NeptuneSeeRaDec(2456789.12345)
ra, dec = NeptuneApparentRaDec(2456789.12345)
fmt.Printf("%.14f\n%.14f\n", ra, dec)
fmt.Println(NeptuneMag(2456789.12345))
}
basic/saturn.go
+300 -10
View File
@@ -63,22 +63,22 @@ func ASaturnXYZ(JD float64) (float64, float64, float64) {
return x, y, z
}
func SaturnSeeRa(JD float64) float64 {
lo, bo := SaturnSeeLoBo(JD)
func SaturnApparentRa(JD float64) float64 {
lo, bo := SaturnApparentLoBo(JD)
sita := Sita(JD)
ra := math.Atan2((Sin(lo)*Cos(sita) - Tan(bo)*Sin(sita)), Cos(lo))
ra = ra * 180 / math.Pi
return Limit360(ra)
}
func SaturnSeeDec(JD float64) float64 {
lo, bo := SaturnSeeLoBo(JD)
func SaturnApparentDec(JD float64) float64 {
lo, bo := SaturnApparentLoBo(JD)
sita := Sita(JD)
dec := ArcSin(Sin(bo)*Cos(sita) + Cos(bo)*Sin(sita)*Sin(lo))
return dec
}
func SaturnSeeRaDec(JD float64) (float64, float64) {
lo, bo := SaturnSeeLoBo(JD)
func SaturnApparentRaDec(JD float64) (float64, float64) {
lo, bo := SaturnApparentLoBo(JD)
sita := Sita(JD)
ra := math.Atan2((Sin(lo)*Cos(sita) - Tan(bo)*Sin(sita)), Cos(lo))
ra = ra * 180 / math.Pi
@@ -92,7 +92,7 @@ func EarthSaturnAway(JD float64) float64 {
return to
}
func SaturnSeeLo(JD float64) float64 {
func SaturnApparentLo(JD float64) float64 {
x, y, z := ASaturnXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = ASaturnXYZ(JD - to)
@@ -107,7 +107,7 @@ func SaturnSeeLo(JD float64) float64 {
return lo
}
func SaturnSeeBo(JD float64) float64 {
func SaturnApparentBo(JD float64) float64 {
x, y, z := ASaturnXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = ASaturnXYZ(JD - to)
@@ -121,7 +121,7 @@ func SaturnSeeBo(JD float64) float64 {
return bo
}
func SaturnSeeLoBo(JD float64) (float64, float64) {
func SaturnApparentLoBo(JD float64) (float64, float64) {
x, y, z := ASaturnXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = ASaturnXYZ(JD - to)
@@ -150,7 +150,297 @@ func SaturnRingB(JD float64) float64 {
T := (JD - 2451545) / 36525
i := 28.075216 - 0.012998*T + 0.000004*T*T
omi := 169.508470 + 1.394681*T + 0.000412*T*T
lo, bo := SaturnSeeLoBo(JD)
lo, bo := SaturnApparentLoBo(JD)
B := Sin(i)*Cos(bo)*Sin(lo-omi) - Cos(i)*Cos(bo)
return ArcSin(B)
}
func SaturnHeight(jde, lon, lat, timezone float64) float64 {
// 转换为世界时
utcJde := jde - timezone/24.0
// 计算视恒星时
ra, dec := SaturnApparentRaDec(TD2UT(utcJde, true))
st := Limit360(ApparentSiderealTime(utcJde)*15 + lon)
// 计算时角
H := Limit360(st - ra)
// 高度角、时角与天球座标三角转换公式
// sin(h)=sin(lat)*sin(dec)+cos(dec)*cos(lat)*cos(H)
sinHeight := Sin(lat)*Sin(dec) + Cos(dec)*Cos(lat)*Cos(H)
return ArcSin(sinHeight)
}
func SaturnAzimuth(jde, lon, lat, timezone float64) float64 {
// 转换为世界时
utcJde := jde - timezone/24.0
// 计算视恒星时
ra, dec := SaturnApparentRaDec(TD2UT(utcJde, true))
st := Limit360(ApparentSiderealTime(utcJde)*15 + lon)
// 计算时角
H := Limit360(st - ra)
// 三角转换公式
tanAzimuth := Sin(H) / (Cos(H)*Sin(lat) - Tan(dec)*Cos(lat))
Azimuth := ArcTan(tanAzimuth)
if Azimuth < 0 {
if H/15 < 12 {
return Azimuth + 360
}
return Azimuth + 180
}
if H/15 < 12 {
return Azimuth + 180
}
return Azimuth
}
func SaturnHourAngle(JD, Lon, TZ float64) float64 {
startime := Limit360(ApparentSiderealTime(JD-TZ/24)*15 + Lon)
timeangle := startime - SaturnApparentRa(TD2UT(JD-TZ/24.0, true))
if timeangle < 0 {
timeangle += 360
}
return timeangle
}
func SaturnCulminationTime(jde, lon, timezone float64) float64 {
//jde 世界时，非力学时，当地时区 0时，无需转换力学时
//ra,dec 瞬时天球座标，非J2000等时间天球坐标
jde = math.Floor(jde) + 0.5
JD1 := jde + Limit360(360-SaturnHourAngle(jde, lon, timezone))/15.0/24.0*0.99726851851851851851
limitHA := func(jde, lon, timezone float64) float64 {
ha := SaturnHourAngle(jde, lon, timezone)
if ha < 180 {
ha += 360
}
return ha
}
for {
JD0 := JD1
stDegree := limitHA(JD0, lon, timezone) - 360
stDegreep := (limitHA(JD0+0.000005, lon, timezone) - limitHA(JD0-0.000005, lon, timezone)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return JD1
}
func SaturnRiseTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 {
return saturnRiseDown(JD, Lon, Lat, TZ, ZS, HEI, true)
}
func SaturnDownTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 {
return saturnRiseDown(JD, Lon, Lat, TZ, ZS, HEI, false)
}
func saturnRiseDown(JD, Lon, Lat, TZ, ZS, HEI float64, isRise bool) float64 {
var An float64
JD = math.Floor(JD) + 0.5
ntz := math.Round(Lon / 15)
if ZS != 0 {
An = -0.8333
}
An = An - HeightDegreeByLat(HEI, Lat)
tztime := SaturnCulminationTime(JD, Lon, ntz)
if SaturnHeight(tztime, Lon, Lat, ntz) < An {
return -2 //极夜
}
if SaturnHeight(tztime-0.5, Lon, Lat, ntz) > An {
return -1 //极昼
}
dec := HSunApparentDec(TD2UT(tztime-ntz/24, true))
//(sin(ho)-sin(φ)*sin(δ2))/(cos(φ)*cos(δ2))
tmp := (Sin(An) - Sin(dec)*Sin(Lat)) / (Cos(dec) * Cos(Lat))
var rise float64
if math.Abs(tmp) <= 1 {
rzsc := ArcCos(tmp) / 15
if isRise {
rise = tztime - rzsc/24 - 25.0/24.0/60.0
} else {
rise = tztime + rzsc/24 - 25.0/24.0/60.0
}
} else {
rise = tztime
i := 0
//TODO:使用二分法计算
for SaturnHeight(rise, Lon, Lat, ntz) > An {
i++
if isRise {
rise -= 15.0 / 60.0 / 24.0
} else {
rise += 15.0 / 60.0 / 24.0
}
if i > 48 {
break
}
}
}
JD1 := rise
for {
JD0 := JD1
stDegree := SaturnHeight(JD0, Lon, Lat, ntz) - An
stDegreep := (SaturnHeight(JD0+0.000005, Lon, Lat, ntz) - SaturnHeight(JD0-0.000005, Lon, Lat, ntz)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return JD1 - ntz/24 + TZ/24
}
// Pos
const SATURN_S_PERIOD = 1 / ((1 / 365.256363004) - (1 / 10759.0))
func saturnConjunction(jde, degree float64, next uint8) float64 {
//0=last 1=next
decSub := func(jde float64, degree float64, filter bool) float64 {
sub := Limit360(Limit360(SaturnApparentLo(jde)-HSunApparentLo(jde)) - degree)
if filter {
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
}
return sub
}
dayCost := SATURN_S_PERIOD / 360
nowSub := decSub(jde, degree, false)
if next == 0 {
jde -= (360 - nowSub) * dayCost
} else {
jde += dayCost * nowSub
}
JD1 := jde
for {
JD0 := JD1
stDegree := decSub(JD0, degree, true)
stDegreep := (decSub(JD0+0.000005, degree, true) - decSub(JD0-0.000005, degree, true)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return TD2UT(JD1, false)
}
func LastSaturnConjunction(jde float64) float64 {
return saturnConjunction(jde, 0, 0)
}
func NextSaturnConjunction(jde float64) float64 {
return saturnConjunction(jde, 0, 1)
}
func LastSaturnOpposition(jde float64) float64 {
return saturnConjunction(jde, 180, 0)
}
func NextSaturnOpposition(jde float64) float64 {
return saturnConjunction(jde, 180, 1)
}
func NextSaturnEasternQuadrature(jde float64) float64 {
return saturnConjunction(jde, 90, 1)
}
func LastSaturnEasternQuadrature(jde float64) float64 {
return saturnConjunction(jde, 90, 0)
}
func NextSaturnWesternQuadrature(jde float64) float64 {
return saturnConjunction(jde, 270, 1)
}
func LastSaturnWesternQuadrature(jde float64) float64 {
return saturnConjunction(jde, 270, 0)
}
func saturnRetrograde(jde float64, isLeft bool) float64 {
//0=last 1=next
decSub := func(jde float64, val float64) float64 {
sub := SaturnApparentRa(jde+val) - SaturnApparentRa(jde-val)
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
return sub / (2 * val)
}
jde = NextSaturnOpposition(jde)
if isLeft {
jde -= 60
} else {
jde += 60
}
for {
nowSub := decSub(jde, 1.0/86400.0)
if math.Abs(nowSub) > 0.55 {
jde += 2
continue
}
break
}
JD1 := jde
for {
JD0 := JD1
stDegree := decSub(JD0, 2.0/86400.0)
stDegreep := (decSub(JD0+15.0/86400.0, 2.0/86400.0) - decSub(JD0-15.0/86400.0, 2.0/86400.0)) / (30.0 / 86400.0)
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 30.0/86400.0 {
break
}
}
JD1 = JD1 - 15.0/86400.0
min := JD1
minRa := 100.0
for i := 0.0; i < 60.0; i++ {
tmp := decSub(JD1+i*0.5/86400.0, 0.5/86400.0)
if math.Abs(tmp) < math.Abs(minRa) {
minRa = tmp
min = JD1 + i*0.5/86400.0
}
}
return TD2UT(min, false)
}
func NextSaturnRetrogradeToPrograde(jde float64) float64 {
date := saturnRetrograde(jde, false)
if date < jde {
op := NextSaturnOpposition(jde)
return saturnRetrograde(op+10, false)
}
return date
}
func LastSaturnRetrogradeToPrograde(jde float64) float64 {
jde = LastSaturnOpposition(jde) - 10
date := saturnRetrograde(jde, false)
if date > jde {
op := LastSaturnOpposition(jde)
return saturnRetrograde(op-10, false)
}
return date
}
func NextSaturnProgradeToRetrograde(jde float64) float64 {
date := saturnRetrograde(jde, true)
if date < jde {
op := NextSaturnOpposition(jde)
return saturnRetrograde(op+10, true)
}
return date
}
func LastSaturnProgradeToRetrograde(jde float64) float64 {
jde = LastSaturnOpposition(jde) - 10
date := saturnRetrograde(jde, true)
if date > jde {
op := LastSaturnOpposition(jde)
return saturnRetrograde(op-10, true)
}
return date
}
basic/saturn_test.go
+13
View File
@@ -0,0 +1,13 @@
package basic
import (
"fmt"
"testing"
)
func TestSaturn(t *testing.T) {
jde := GetNowJDE() - 6000
for i := 0.00; i < 20; i++ {
fmt.Println(jde+i*365, JDE2Date(jde+i*365), JDE2Date(LastSaturnProgradeToRetrograde(jde+i*365)))
}
}
basic/star.go
+13 -1
View File
@@ -150,7 +150,7 @@ func StarCulminationTime(jde, ra, lon, timezone float64) float64 {
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
stDegreep := (limitStarHA(JD0+0.000005, ra, lon, timezone) - limitStarHA(JD0-0.000005, ra, lon, timezone)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
@@ -158,3 +158,15 @@ func StarCulminationTime(jde, ra, lon, timezone float64) float64 {
}
return JD1
}
func StarAngularSeparation(ra1, dec1, ra2, dec2 float64) float64 {
//cos(d)=sinδ1 sinδ2 + cosδ1 cosδ2 cos(α1-α2)
d := Sin(dec1)*Sin(dec2) + Cos(dec1)*Cos(dec2)*Cos(ra1-ra2)
if math.Abs(d) >= 0.999999997 {
//d = √(Δα*cosδ)2+(Δδ)2
tmp1 := ((ra1 - ra2) * Cos((dec1+dec2)/2))
tmp2 := (dec1 - dec2)
return math.Sqrt(tmp1*tmp1 + tmp2*tmp2)
}
return ArcCos(d)
}
basic/star_test.go
+5
View File
@@ -20,3 +20,8 @@ func Test_Star(t *testing.T) {
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)))
}
func TestZB(t *testing.T) {
jde := GetNowJDE()
fmt.Println(LoBoToRaDec(jde, 156, 0))
}
basic/sun.go
+6 -40
View File
@@ -838,18 +838,7 @@ func SunApparentLo(JD float64) float64 { //'太阳视黄经
}
func SunApparentRa(JD float64) float64 { // '太阳视赤经
T := (JD - 2451545) / 36525
sitas := Sita(JD) + 0.00256*Cos(125.04-1934.136*T)
SunApparentRa := ArcTan(Cos(sitas) * Sin(SunApparentLo(JD)) / Cos(SunApparentLo(JD)))
tmp := SunApparentLo(JD)
if tmp >= 90 && tmp < 180 {
SunApparentRa = 180 + SunApparentRa
} else if tmp >= 180 && tmp < 270 {
SunApparentRa = 180 + SunApparentRa
} else if tmp >= 270 && tmp <= 360 {
SunApparentRa = 360 + SunApparentRa
}
return SunApparentRa
return LoToRa(JD, SunApparentLo(JD), 0)
}
func SunTrueRa(JD float64) float64 { //'太阳真赤经
@@ -934,35 +923,11 @@ func EarthAway(JD float64) float64 {
}
func HSunApparentRaDec(JD float64) (float64, float64) {
T := (JD - 2451545) / 36525
sitas := Sita(JD) + 0.00256*Cos(125.04-1934.136*T)
sitas2 := EclipticObliquity(JD, false) + 0.00256*Cos(125.04-1934.136*T)
tmp := HSunApparentLo(JD)
HSunApparentRa := ArcTan(Cos(sitas) * Sin(tmp) / Cos(tmp))
HSunApparentDec := ArcSin(Sin(sitas2) * Sin(tmp))
if tmp >= 90 && tmp < 180 {
HSunApparentRa = 180 + HSunApparentRa
} else if tmp >= 180 && tmp < 270 {
HSunApparentRa = 180 + HSunApparentRa
} else if tmp >= 270 && tmp <= 360 {
HSunApparentRa = 360 + HSunApparentRa
}
return HSunApparentRa, HSunApparentDec
return LoBoToRaDec(JD, HSunApparentLo(JD), HSunTrueBo(JD))
}
func HSunApparentRa(JD float64) float64 { // '太阳视赤经
T := (JD - 2451545) / 36525
sitas := Sita(JD) + 0.00256*Cos(125.04-1934.136*T)
tmp := HSunApparentLo(JD)
HSunApparentRa := ArcTan(Cos(sitas) * Sin(tmp) / Cos(tmp))
if tmp >= 90 && tmp < 180 {
HSunApparentRa = 180 + HSunApparentRa
} else if tmp >= 180 && tmp < 270 {
HSunApparentRa = 180 + HSunApparentRa
} else if tmp >= 270 && tmp <= 360 {
HSunApparentRa = 360 + HSunApparentRa
}
return HSunApparentRa
return LoToRa(JD, HSunApparentLo(JD), HSunTrueBo(JD))
}
func HSunTrueRa(JD float64) float64 { //'太阳真赤经
@@ -1030,10 +995,10 @@ func GetOneYearMoon(year float64) map[int]float64 {
}
return moon
}
func GetOneYearJQ(year int) map[int]float64 {
func GetOneYearJQ(year int) []float64 {
start := 270
var years int
jq := make(map[int]float64)
jq := make([]float64, 26)
for i := 1; i < 26; i++ {
angle := start + 15*(i-1)
if angle > 360 {
@@ -1047,6 +1012,7 @@ func GetOneYearJQ(year int) map[int]float64 {
jq[i] = GetJQTime(years, angle) + 8.0/24.0
// echo DateCalc(jq[i])."<br />";
}
jq[0] = jq[1]
return jq
}
basic/uranus.go
+299 -9
View File
@@ -63,22 +63,22 @@ func AUranusXYZ(JD float64) (float64, float64, float64) {
return x, y, z
}
func UranusSeeRa(JD float64) float64 {
lo, bo := UranusSeeLoBo(JD)
func UranusApparentRa(JD float64) float64 {
lo, bo := UranusApparentLoBo(JD)
sita := Sita(JD)
ra := math.Atan2((Sin(lo)*Cos(sita) - Tan(bo)*Sin(sita)), Cos(lo))
ra = ra * 180 / math.Pi
return Limit360(ra)
}
func UranusSeeDec(JD float64) float64 {
lo, bo := UranusSeeLoBo(JD)
func UranusApparentDec(JD float64) float64 {
lo, bo := UranusApparentLoBo(JD)
sita := Sita(JD)
dec := ArcSin(Sin(bo)*Cos(sita) + Cos(bo)*Sin(sita)*Sin(lo))
return dec
}
func UranusSeeRaDec(JD float64) (float64, float64) {
lo, bo := UranusSeeLoBo(JD)
func UranusApparentRaDec(JD float64) (float64, float64) {
lo, bo := UranusApparentLoBo(JD)
sita := Sita(JD)
ra := math.Atan2((Sin(lo)*Cos(sita) - Tan(bo)*Sin(sita)), Cos(lo))
ra = ra * 180 / math.Pi
@@ -92,7 +92,7 @@ func EarthUranusAway(JD float64) float64 {
return to
}
func UranusSeeLo(JD float64) float64 {
func UranusApparentLo(JD float64) float64 {
x, y, z := AUranusXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = AUranusXYZ(JD - to)
@@ -107,7 +107,7 @@ func UranusSeeLo(JD float64) float64 {
return lo
}
func UranusSeeBo(JD float64) float64 {
func UranusApparentBo(JD float64) float64 {
x, y, z := AUranusXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = AUranusXYZ(JD - to)
@@ -121,7 +121,7 @@ func UranusSeeBo(JD float64) float64 {
return bo
}
func UranusSeeLoBo(JD float64) (float64, float64) {
func UranusApparentLoBo(JD float64) (float64, float64) {
x, y, z := AUranusXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = AUranusXYZ(JD - to)
@@ -145,3 +145,293 @@ func UranusMag(JD float64) float64 {
Mag := -7.19 + 5*math.Log10(AwaySun*AwayEarth) + 0.016*i
return FloatRound(Mag, 2)
}
func UranusHeight(jde, lon, lat, timezone float64) float64 {
// 转换为世界时
utcJde := jde - timezone/24.0
// 计算视恒星时
ra, dec := UranusApparentRaDec(TD2UT(utcJde, true))
st := Limit360(ApparentSiderealTime(utcJde)*15 + lon)
// 计算时角
H := Limit360(st - ra)
// 高度角、时角与天球座标三角转换公式
// sin(h)=sin(lat)*sin(dec)+cos(dec)*cos(lat)*cos(H)
sinHeight := Sin(lat)*Sin(dec) + Cos(dec)*Cos(lat)*Cos(H)
return ArcSin(sinHeight)
}
func UranusAzimuth(jde, lon, lat, timezone float64) float64 {
// 转换为世界时
utcJde := jde - timezone/24.0
// 计算视恒星时
ra, dec := UranusApparentRaDec(TD2UT(utcJde, true))
st := Limit360(ApparentSiderealTime(utcJde)*15 + lon)
// 计算时角
H := Limit360(st - ra)
// 三角转换公式
tanAzimuth := Sin(H) / (Cos(H)*Sin(lat) - Tan(dec)*Cos(lat))
Azimuth := ArcTan(tanAzimuth)
if Azimuth < 0 {
if H/15 < 12 {
return Azimuth + 360
}
return Azimuth + 180
}
if H/15 < 12 {
return Azimuth + 180
}
return Azimuth
}
func UranusHourAngle(JD, Lon, TZ float64) float64 {
startime := Limit360(ApparentSiderealTime(JD-TZ/24)*15 + Lon)
timeangle := startime - UranusApparentRa(TD2UT(JD-TZ/24.0, true))
if timeangle < 0 {
timeangle += 360
}
return timeangle
}
func UranusCulminationTime(jde, lon, timezone float64) float64 {
//jde 世界时，非力学时，当地时区 0时，无需转换力学时
//ra,dec 瞬时天球座标，非J2000等时间天球坐标
jde = math.Floor(jde) + 0.5
JD1 := jde + Limit360(360-UranusHourAngle(jde, lon, timezone))/15.0/24.0*0.99726851851851851851
limitHA := func(jde, lon, timezone float64) float64 {
ha := UranusHourAngle(jde, lon, timezone)
if ha < 180 {
ha += 360
}
return ha
}
for {
JD0 := JD1
stDegree := limitHA(JD0, lon, timezone) - 360
stDegreep := (limitHA(JD0+0.000005, lon, timezone) - limitHA(JD0-0.000005, lon, timezone)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return JD1
}
func UranusRiseTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 {
return uranusRiseDown(JD, Lon, Lat, TZ, ZS, HEI, true)
}
func UranusDownTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 {
return uranusRiseDown(JD, Lon, Lat, TZ, ZS, HEI, false)
}
func uranusRiseDown(JD, Lon, Lat, TZ, ZS, HEI float64, isRise bool) float64 {
var An float64
JD = math.Floor(JD) + 0.5
ntz := math.Round(Lon / 15)
if ZS != 0 {
An = -0.8333
}
An = An - HeightDegreeByLat(HEI, Lat)
tztime := UranusCulminationTime(JD, Lon, ntz)
if UranusHeight(tztime, Lon, Lat, ntz) < An {
return -2 //极夜
}
if UranusHeight(tztime-0.5, Lon, Lat, ntz) > An {
return -1 //极昼
}
dec := HSunApparentDec(TD2UT(tztime-ntz/24, true))
//(sin(ho)-sin(φ)*sin(δ2))/(cos(φ)*cos(δ2))
tmp := (Sin(An) - Sin(dec)*Sin(Lat)) / (Cos(dec) * Cos(Lat))
var rise float64
if math.Abs(tmp) <= 1 {
rzsc := ArcCos(tmp) / 15
if isRise {
rise = tztime - rzsc/24 - 25.0/24.0/60.0
} else {
rise = tztime + rzsc/24 - 25.0/24.0/60.0
}
} else {
rise = tztime
i := 0
//TODO:使用二分法计算
for UranusHeight(rise, Lon, Lat, ntz) > An {
i++
if isRise {
rise -= 15.0 / 60.0 / 24.0
} else {
rise += 15.0 / 60.0 / 24.0
}
if i > 48 {
break
}
}
}
JD1 := rise
for {
JD0 := JD1
stDegree := UranusHeight(JD0, Lon, Lat, ntz) - An
stDegreep := (UranusHeight(JD0+0.000005, Lon, Lat, ntz) - UranusHeight(JD0-0.000005, Lon, Lat, ntz)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return JD1 - ntz/24 + TZ/24
}
// Pos
const URANUS_S_PERIOD = 1 / ((1 / 365.256363004) - (1 / 30799.095))
func uranusConjunction(jde, degree float64, next uint8) float64 {
//0=last 1=next
decSub := func(jde float64, degree float64, filter bool) float64 {
sub := Limit360(Limit360(UranusApparentLo(jde)-HSunApparentLo(jde)) - degree)
if filter {
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
}
return sub
}
dayCost := URANUS_S_PERIOD / 360
nowSub := decSub(jde, degree, false)
if next == 0 {
jde -= (360 - nowSub) * dayCost
} else {
jde += dayCost * nowSub
}
JD1 := jde
for {
JD0 := JD1
stDegree := decSub(JD0, degree, true)
stDegreep := (decSub(JD0+0.000005, degree, true) - decSub(JD0-0.000005, degree, true)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return TD2UT(JD1, false)
}
func LastUranusConjunction(jde float64) float64 {
return uranusConjunction(jde, 0, 0)
}
func NextUranusConjunction(jde float64) float64 {
return uranusConjunction(jde, 0, 1)
}
func LastUranusOpposition(jde float64) float64 {
return uranusConjunction(jde, 180, 0)
}
func NextUranusOpposition(jde float64) float64 {
return uranusConjunction(jde, 180, 1)
}
func NextUranusEasternQuadrature(jde float64) float64 {
return uranusConjunction(jde, 90, 1)
}
func LastUranusEasternQuadrature(jde float64) float64 {
return uranusConjunction(jde, 90, 0)
}
func NextUranusWesternQuadrature(jde float64) float64 {
return uranusConjunction(jde, 270, 1)
}
func LastUranusWesternQuadrature(jde float64) float64 {
return uranusConjunction(jde, 270, 0)
}
func uranusRetrograde(jde float64, isLeft bool) float64 {
//0=last 1=next
decSub := func(jde float64, val float64) float64 {
sub := UranusApparentRa(jde+val) - UranusApparentRa(jde-val)
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
return sub / (2 * val)
}
jde = NextUranusOpposition(jde)
if isLeft {
jde -= 60
} else {
jde += 60
}
for {
nowSub := decSub(jde, 1.0/86400.0)
if math.Abs(nowSub) > 0.55 {
jde += 2
continue
}
break
}
JD1 := jde
for {
JD0 := JD1
stDegree := decSub(JD0, 2.0/86400.0)
stDegreep := (decSub(JD0+15.0/86400.0, 2.0/86400.0) - decSub(JD0-15.0/86400.0, 2.0/86400.0)) / (30.0 / 86400.0)
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 30.0/86400.0 {
break
}
}
JD1 = JD1 - 15.0/86400.0
min := JD1
minRa := 100.0
for i := 0.0; i < 60.0; i++ {
tmp := decSub(JD1+i*0.5/86400.0, 0.5/86400.0)
if math.Abs(tmp) < math.Abs(minRa) {
minRa = tmp
min = JD1 + i*0.5/86400.0
}
}
return TD2UT(min, false)
}
func NextUranusRetrogradeToPrograde(jde float64) float64 {
date := uranusRetrograde(jde, false)
if date < jde {
op := NextUranusOpposition(jde)
return uranusRetrograde(op+10, false)
}
return date
}
func LastUranusRetrogradeToPrograde(jde float64) float64 {
jde = LastUranusOpposition(jde) - 10
date := uranusRetrograde(jde, false)
if date > jde {
op := LastUranusOpposition(jde)
return uranusRetrograde(op-10, false)
}
return date
}
func NextUranusProgradeToRetrograde(jde float64) float64 {
date := uranusRetrograde(jde, true)
if date < jde {
op := NextUranusOpposition(jde)
return uranusRetrograde(op+10, true)
}
return date
}
func LastUranusProgradeToRetrograde(jde float64) float64 {
jde = LastUranusOpposition(jde) - 10
date := uranusRetrograde(jde, true)
if date > jde {
op := LastUranusOpposition(jde)
return uranusRetrograde(op-10, true)
}
return date
}
basic/venus.go
+477 -9
View File
@@ -63,22 +63,22 @@ func AVenusXYZ(JD float64) (float64, float64, float64) {
return x, y, z
}
func VenusSeeRa(JD float64) float64 {
lo, bo := VenusSeeLoBo(JD)
func VenusApparentRa(JD float64) float64 {
lo, bo := VenusApparentLoBo(JD)
sita := Sita(JD)
ra := math.Atan2((Sin(lo)*Cos(sita) - Tan(bo)*Sin(sita)), Cos(lo))
ra = ra * 180 / math.Pi
return Limit360(ra)
}
func VenusSeeDec(JD float64) float64 {
lo, bo := VenusSeeLoBo(JD)
func VenusApparentDec(JD float64) float64 {
lo, bo := VenusApparentLoBo(JD)
sita := Sita(JD)
dec := ArcSin(Sin(bo)*Cos(sita) + Cos(bo)*Sin(sita)*Sin(lo))
return dec
}
func VenusSeeRaDec(JD float64) (float64, float64) {
lo, bo := VenusSeeLoBo(JD)
func VenusApparentRaDec(JD float64) (float64, float64) {
lo, bo := VenusApparentLoBo(JD)
sita := Sita(JD)
ra := math.Atan2((Sin(lo)*Cos(sita) - Tan(bo)*Sin(sita)), Cos(lo))
ra = ra * 180 / math.Pi
@@ -92,7 +92,7 @@ func EarthVenusAway(JD float64) float64 {
return to
}
func VenusSeeLo(JD float64) float64 {
func VenusApparentLo(JD float64) float64 {
x, y, z := AVenusXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = AVenusXYZ(JD - to)
@@ -107,7 +107,7 @@ func VenusSeeLo(JD float64) float64 {
return lo
}
func VenusSeeBo(JD float64) float64 {
func VenusApparentBo(JD float64) float64 {
x, y, z := AVenusXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = AVenusXYZ(JD - to)
@@ -121,7 +121,7 @@ func VenusSeeBo(JD float64) float64 {
return bo
}
func VenusSeeLoBo(JD float64) (float64, float64) {
func VenusApparentLoBo(JD float64) (float64, float64) {
x, y, z := AVenusXYZ(JD)
to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
x, y, z = AVenusXYZ(JD - to)
@@ -145,3 +145,471 @@ func VenusMag(JD float64) float64 {
Mag := -4.40 + 5*math.Log10(AwaySun*AwayEarth) + 0.0009*i + 0.000239*i*i - 0.00000065*i*i*i
return FloatRound(Mag, 2)
}
func VenusHeight(jde, lon, lat, timezone float64) float64 {
// 转换为世界时
utcJde := jde - timezone/24.0
// 计算视恒星时
ra, dec := VenusApparentRaDec(TD2UT(utcJde, true))
st := Limit360(ApparentSiderealTime(utcJde)*15 + lon)
// 计算时角
H := Limit360(st - ra)
// 高度角、时角与天球座标三角转换公式
// sin(h)=sin(lat)*sin(dec)+cos(dec)*cos(lat)*cos(H)
sinHeight := Sin(lat)*Sin(dec) + Cos(dec)*Cos(lat)*Cos(H)
return ArcSin(sinHeight)
}
func VenusAzimuth(jde, lon, lat, timezone float64) float64 {
// 转换为世界时
utcJde := jde - timezone/24.0
// 计算视恒星时
ra, dec := VenusApparentRaDec(TD2UT(utcJde, true))
st := Limit360(ApparentSiderealTime(utcJde)*15 + lon)
// 计算时角
H := Limit360(st - ra)
// 三角转换公式
tanAzimuth := Sin(H) / (Cos(H)*Sin(lat) - Tan(dec)*Cos(lat))
Azimuth := ArcTan(tanAzimuth)
if Azimuth < 0 {
if H/15 < 12 {
return Azimuth + 360
}
return Azimuth + 180
}
if H/15 < 12 {
return Azimuth + 180
}
return Azimuth
}
func VenusHourAngle(JD, Lon, TZ float64) float64 {
startime := Limit360(ApparentSiderealTime(JD-TZ/24)*15 + Lon)
timeangle := startime - VenusApparentRa(TD2UT(JD-TZ/24.0, true))
if timeangle < 0 {
timeangle += 360
}
return timeangle
}
func VenusCulminationTime(jde, lon, timezone float64) float64 {
//jde 世界时，非力学时，当地时区 0时，无需转换力学时
//ra,dec 瞬时天球座标，非J2000等时间天球坐标
jde = math.Floor(jde) + 0.5
JD1 := jde + Limit360(360-VenusHourAngle(jde, lon, timezone))/15.0/24.0*0.99726851851851851851
limitHA := func(jde, lon, timezone float64) float64 {
ha := VenusHourAngle(jde, lon, timezone)
if ha < 180 {
ha += 360
}
return ha
}
for {
JD0 := JD1
stDegree := limitHA(JD0, lon, timezone) - 360
stDegreep := (limitHA(JD0+0.000005, lon, timezone) - limitHA(JD0-0.000005, lon, timezone)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return JD1
}
func VenusRiseTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 {
return venusRiseDown(JD, Lon, Lat, TZ, ZS, HEI, true)
}
func VenusDownTime(JD, Lon, Lat, TZ, ZS, HEI float64) float64 {
return venusRiseDown(JD, Lon, Lat, TZ, ZS, HEI, false)
}
func venusRiseDown(JD, Lon, Lat, TZ, ZS, HEI float64, isRise bool) float64 {
var An float64
JD = math.Floor(JD) + 0.5
ntz := math.Round(Lon / 15)
if ZS != 0 {
An = -0.8333
}
An = An - HeightDegreeByLat(HEI, Lat)
tztime := VenusCulminationTime(JD, Lon, ntz)
if VenusHeight(tztime, Lon, Lat, ntz) < An {
return -2 //极夜
}
if VenusHeight(tztime-0.5, Lon, Lat, ntz) > An {
return -1 //极昼
}
dec := HSunApparentDec(TD2UT(tztime-ntz/24, true))
//(sin(ho)-sin(φ)*sin(δ2))/(cos(φ)*cos(δ2))
tmp := (Sin(An) - Sin(dec)*Sin(Lat)) / (Cos(dec) * Cos(Lat))
var rise float64
if math.Abs(tmp) <= 1 {
rzsc := ArcCos(tmp) / 15
if isRise {
rise = tztime - rzsc/24 - 25.0/24.0/60.0
} else {
rise = tztime + rzsc/24 - 25.0/24.0/60.0
}
} else {
rise = tztime
i := 0
//TODO:使用二分法计算
for VenusHeight(rise, Lon, Lat, ntz) > An {
i++
if isRise {
rise -= 15.0 / 60.0 / 24.0
} else {
rise += 15.0 / 60.0 / 24.0
}
if i > 48 {
break
}
}
}
JD1 := rise
for {
JD0 := JD1
stDegree := VenusHeight(JD0, Lon, Lat, ntz) - An
stDegreep := (VenusHeight(JD0+0.000005, Lon, Lat, ntz) - VenusHeight(JD0-0.000005, Lon, Lat, ntz)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return JD1 - ntz/24 + TZ/24
}
// Pos
const VENUS_S_PERIOD = 1 / ((1 / 224.701) - (1 / 365.256363004))
func venusConjunction(jde float64, next uint8) float64 {
//0=last 1=next
decSub := func(jde float64) float64 {
sub := Limit360(VenusApparentLo(jde) - HSunApparentLo(jde))
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
return sub
}
nowSub := decSub(jde)
pos := math.Abs(decSub(jde+1/86400.0)) - math.Abs(nowSub)
if pos >= 0 && next == 1 && nowSub > 0 {
jde += VENUS_S_PERIOD/8.0 + 2
}
if pos >= 0 && next == 1 && nowSub < 0 {
jde += VENUS_S_PERIOD/6.0 + 2
}
if pos <= 0 && next == 0 && nowSub < 0 {
jde -= VENUS_S_PERIOD/8.0 + 2
}
if pos <= 0 && next == 0 && nowSub > 0 {
jde -= VENUS_S_PERIOD/6.0 + 2
}
for {
nowSub := decSub(jde)
pos := math.Abs(decSub(jde+1/86400.0)) - math.Abs(nowSub)
if math.Abs(nowSub) > 24 || (pos > 0 && next == 1) || (pos < 0 && next == 0) {
if next == 1 {
jde += 8
} else {
jde -= 8
}
continue
}
break
}
JD1 := jde
for {
JD0 := JD1
stDegree := decSub(JD0)
stDegreep := (decSub(JD0+0.000005) - decSub(JD0-0.000005)) / 0.00001
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 0.00001 {
break
}
}
return TD2UT(JD1, false)
}
func LastVenusConjunction(jde float64) float64 {
return venusConjunction(jde, 0)
}
func NextVenusConjunction(jde float64) float64 {
return venusConjunction(jde, 1)
}
func NextVenusInferiorConjunction(jde float64) float64 {
date := NextVenusConjunction(jde)
if EarthVenusAway(date) > EarthAway(date) {
return NextVenusConjunction(date + 2)
}
return date
}
func NextVenusSuperiorConjunction(jde float64) float64 {
date := NextVenusConjunction(jde)
if EarthVenusAway(date) < EarthAway(date) {
return NextVenusConjunction(date + 2)
}
return date
}
func LastVenusInferiorConjunction(jde float64) float64 {
date := LastVenusConjunction(jde)
if EarthVenusAway(date) > EarthAway(date) {
return LastVenusConjunction(date - 2)
}
return date
}
func LastVenusSuperiorConjunction(jde float64) float64 {
date := LastVenusConjunction(jde)
if EarthVenusAway(date) < EarthAway(date) {
return LastVenusConjunction(date - 2)
}
return date
}
func venusRetrograde(jde float64) float64 {
//0=last 1=next
decSunSub := func(jde float64) float64 {
sub := Limit360(VenusApparentRa(jde) - SunApparentRa(jde))
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
return sub
}
decSub := func(jde float64, val float64) float64 {
sub := VenusApparentRa(jde+val) - VenusApparentRa(jde-val)
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
return sub / (2 * val)
}
lastHe := LastVenusConjunction(jde)
nextHe := NextVenusConjunction(jde)
nowSub := decSunSub(jde)
if nowSub > 0 {
jde = lastHe + ((nextHe - lastHe) / 5.0 * 3.5)
} else {
jde = lastHe + 10
}
for {
nowSub := decSub(jde, 1.0/86400.0)
if math.Abs(nowSub) > 0.5 {
jde += 5
continue
}
break
}
JD1 := jde
for {
JD0 := JD1
stDegree := decSub(JD0, 0.5/86400.0)
stDegreep := (decSub(JD0+10.0/86400.0, 0.5/86400.0) - decSub(JD0-10.0/86400.0, 0.5/86400.0)) / (20.0 / 86400.0)
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 20.0/86400.0 {
break
}
}
JD1 = JD1 - 10.0/86400.0
min := JD1
minRa := 100.0
for i := 0.0; i < 40.0; i++ {
tmp := decSub(JD1+i*0.5/86400.0, 0.5/86400.0)
if math.Abs(tmp) < math.Abs(minRa) {
minRa = tmp
min = JD1 + i*0.5/86400.0
}
}
//fmt.Println((min - lastHe) / (nextHe - lastHe))
return TD2UT(min, false)
}
func NextVenusRetrograde(jde float64) float64 {
date := venusRetrograde(jde)
if date < jde {
nextHe := NextVenusConjunction(jde)
return venusRetrograde(nextHe + 2)
}
return date
}
func LastVenusRetrograde(jde float64) float64 {
lastHe := LastVenusConjunction(jde)
date := venusRetrograde(lastHe + 2)
if date > jde {
lastLastHe := LastVenusConjunction(lastHe - 2)
return venusRetrograde(lastLastHe + 2)
}
return date
}
func NextVenusProgradeToRetrograde(jde float64) float64 {
date := NextVenusRetrograde(jde)
sub := Limit360(VenusApparentRa(date) - SunApparentRa(date))
if sub > 180 {
return NextVenusRetrograde(date + VENUS_S_PERIOD/2)
}
return date
}
func NextVenusRetrogradeToPrograde(jde float64) float64 {
date := NextVenusRetrograde(jde)
sub := Limit360(VenusApparentRa(date) - SunApparentRa(date))
if sub < 180 {
return NextVenusRetrograde(date + 12)
}
return date
}
func LastVenusProgradeToRetrograde(jde float64) float64 {
date := LastVenusRetrograde(jde)
sub := Limit360(VenusApparentRa(date) - SunApparentRa(date))
if sub > 180 {
return LastVenusRetrograde(date - 12)
}
return date
}
func LastVenusRetrogradeToPrograde(jde float64) float64 {
date := LastVenusRetrograde(jde)
sub := Limit360(VenusApparentRa(date) - SunApparentRa(date))
if sub < 180 {
return LastVenusRetrograde(date - VENUS_S_PERIOD/2)
}
return date
}
func VenusSunElongation(jde float64) float64 {
lo1, bo1 := VenusApparentLoBo(jde)
lo2 := SunApparentLo(jde)
bo2 := HSunTrueBo(jde)
return StarAngularSeparation(lo1, bo1, lo2, bo2)
}
func venusGreatestElongation(jde float64) float64 {
decSunSub := func(jde float64) float64 {
sub := Limit360(VenusApparentRa(jde) - SunApparentRa(jde))
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
return sub
}
decSub := func(jde float64, val float64) float64 {
sub := VenusSunElongation(jde+val) - VenusSunElongation(jde-val)
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
return sub / (2 * val)
}
lastHe := LastVenusConjunction(jde)
nextHe := NextVenusConjunction(jde)
nowSub := decSunSub(jde)
if nowSub > 0 {
jde = lastHe + ((nextHe - lastHe) / 5.0 * 2.5)
} else {
jde = lastHe + ((nextHe - lastHe) / 5.0)
}
for {
nowSub := decSub(jde, 1.0/86400.0)
if math.Abs(nowSub) > 0.15 {
jde += 5
continue
}
break
}
JD1 := jde
for {
JD0 := JD1
stDegree := decSub(JD0, 2.0/86400.0)
stDegreep := (decSub(JD0+15.0/86400.0, 2.0/86400.0) - decSub(JD0-15.0/86400.0, 2.0/86400.0)) / (30.0 / 86400.0)
JD1 = JD0 - stDegree/stDegreep
if math.Abs(JD1-JD0) <= 30.0/86400.0 {
break
}
}
JD1 = JD1 - 15.0/86400.0
min := JD1
minRa := 100.0
for i := 0.0; i < 60.0; i++ {
tmp := decSub(JD1+i*0.5/86400.0, 0.5/86400.0)
if math.Abs(tmp) < math.Abs(minRa) {
minRa = tmp
min = JD1 + i*0.5/86400.0
}
}
//fmt.Println((min - lastHe) / (nextHe - lastHe))
return TD2UT(min, false)
}
func NextVenusGreatestElongation(jde float64) float64 {
date := venusGreatestElongation(jde)
if date < jde {
nextHe := NextVenusConjunction(jde)
return venusGreatestElongation(nextHe + 2)
}
return date
}
func LastVenusGreatestElongation(jde float64) float64 {
lastHe := LastVenusConjunction(jde)
date := venusGreatestElongation(lastHe + 2)
if date > jde {
lastLastHe := LastVenusConjunction(lastHe - 2)
return venusGreatestElongation(lastLastHe + 2)
}
return date
}
func NextVenusGreatestElongationEast(jde float64) float64 {
date := NextVenusGreatestElongation(jde)
sub := Limit360(VenusApparentRa(date) - SunApparentRa(date))
if sub > 180 {
return NextVenusGreatestElongation(date + 1)
}
return date
}
func NextVenusGreatestElongationWest(jde float64) float64 {
date := NextVenusGreatestElongation(jde)
sub := Limit360(VenusApparentRa(date) - SunApparentRa(date))
if sub < 180 {
return NextVenusGreatestElongation(date + 1)
}
return date
}
func LastVenusGreatestElongationEast(jde float64) float64 {
date := LastVenusGreatestElongation(jde)
sub := Limit360(VenusApparentRa(date) - SunApparentRa(date))
if sub > 180 {
return LastVenusGreatestElongation(date - 1)
}
return date
}
func LastVenusGreatestElongationWest(jde float64) float64 {
date := LastVenusGreatestElongation(jde)
sub := Limit360(VenusApparentRa(date) - SunApparentRa(date))
if sub < 180 {
return LastVenusGreatestElongation(date - 1)
}
return date
}
basic/venus_test.go
+37
View File
@@ -0,0 +1,37 @@
package basic
import (
"fmt"
"testing"
)
func TestVenus(t *testing.T) {
jde := 2.4597161573032406e+06 - 720
/*
fmt.Println(JDE2Date(VenusCulminationTime(jde, 115, 8)))
fmt.Println(JDE2Date(VenusRiseTime(jde, 115, 23, 8, 0, 0)))
fmt.Println(JDE2Date(VenusDownTime(jde, 115, 23, 8, 0, 0)))
fmt.Println("----------------")
*/
//LastVenusConjunction(2.4596600340162036e+06)
//fmt.Println(2.4590359532407406e+06, JDE2Date(2.4590359532407406e+06), JDE2Date(NextVenusRetrograde(2.4590359532407406e+06)))
//fmt.Println(jde)
///fmt.Println(MarsTrueLoBo(jde))
//fmt.Println((jde-2451545)/36525, JDE2Date(0.2293425175054224*36525+2451545))
decSub := func(jde float64, val float64) float64 {
sub := VenusSunElongation(jde+val) - VenusSunElongation(jde-val)
if sub > 180 {
sub -= 360
}
if sub < -180 {
sub += 360
}
return sub / (2 * val)
}
_ = decSub
for i := 0.00; i < 1800.0; i += 50 {
fmt.Println(jde+i, JDE2Date(jde+i), JDE2Date(LastVenusGreatestElongationWest(jde+i)))
//fmt.Println(decSub(jde+i, 1.0/86400.0))
}
}