diff --git a/pkg/calculator/sofa.go b/pkg/calculator/sofa.go index 1d7abe1..85b853a 100644 --- a/pkg/calculator/sofa.go +++ b/pkg/calculator/sofa.go @@ -8,7 +8,7 @@ import ( ) func ECItoECEF(j2000X, j2000Y, j2000Z float64, utc time.Time) (float64, float64, float64) { - j2000Position := [3]float64{j2000X, j2000Y, j2000Z} + eci := [3]float64{j2000X, j2000Y, j2000Z} var ecef [3]float64 // 当前的 TT 时间 (Terrestrial Time) @@ -43,10 +43,54 @@ func ECItoECEF(j2000X, j2000Y, j2000Z float64, utc time.Time) (float64, float64, gofa.C2t06a(tt1, tt2, ut11, ut12, xp, yp, &rc2t) // 计算 ECEF 坐标 - gofa.Rxp(rc2t, j2000Position, &ecef) + gofa.Rxp(rc2t, eci, &ecef) return ecef[0], ecef[1], ecef[2] } +func ECEFtoECI(X84, Y84, Z84 float64, utc time.Time) (float64, float64, float64) { + ecef := [3]float64{X84, Y84, Z84} + var eci [3]float64 + + // 当前的 TT 时间 (Terrestrial Time) + var tt1, tt2 float64 + gofa.Dtf2d("UTC", + utc.Year(), + int(utc.Month()), + utc.Day(), + utc.Hour(), + utc.Minute(), + float64(utc.Second()), + &tt1, &tt2) + + var djm0, mjd float64 + gofa.Cal2jd(int(utc.Year()), int(utc.Month()), int(utc.Day()), &djm0, &mjd) + var dut1 float64 = 0.0562 + xp := 0.0 + yp := 0.0 + + if e, ok := EOP.Get(mjd); ok { + xp = e.Xp + yp = e.Yp + dut1 = e.Dut1 + } + + // 计算 UT1 时间 世界时UT1与协调世界时UTC + var ut11, ut12 float64 + gofa.Utcut1(tt1, tt2, dut1, &ut11, &ut12) + + // 计算 C2T 矩阵 (ICRS -> CIRS -> ECEF) + var rc2t [3][3]float64 + gofa.C2t06a(tt1, tt2, ut11, ut12, xp, yp, &rc2t) + + // 计算 ECEF 到 ECI 的转换矩阵 + var rc2tT [3][3]float64 + gofa.Tr(rc2t, &rc2tT) + + // 将 ECEF 坐标转换为 ECI 坐标 + gofa.Rxp(rc2tT, ecef, &eci) + return eci[0], eci[1], eci[2] +} + func ECEFGeocentricToGeodetic(x, y, z float64) (float64, float64, float64) { var lat, lon, alt float64 gofa.Gc2gd(gofa.WGS84, [3]float64{x, y, z}, &lon, &lat, &alt) diff --git a/pkg/producer/aux.go b/pkg/producer/aux.go index c93fa97..171e4f7 100644 --- a/pkg/producer/aux.go +++ b/pkg/producer/aux.go @@ -114,9 +114,17 @@ func (r *Registrator) SetSceneBoundary(scene *Scene) (topLeft, bottomRight orb.P endTime := time.Unix(int64(auxilary.ReferenceTime2000)+int64(ae.UTCTimeSec), int64(ae.Microsecond)*1000).UTC() startPos84 := []float64{as.W84PosX, as.W84PosY, as.W84PosZ} - startPosECI := []float64{as.J2000PosX, as.J2000PosY, as.J2000PosZ} + startPosECI := []float64{ + as.J2000PosX + as.J2000VelX*float64(as.Microsecond)/10e6, + as.J2000PosY + as.J2000VelY*float64(as.Microsecond)/10e6, + as.J2000PosZ + as.J2000VelZ*float64(as.Microsecond)/10e6, + } endPos84 := []float64{ae.W84PosX, ae.W84PosY, ae.W84PosZ} - endPosECI := []float64{ae.J2000PosX, ae.J2000PosY, ae.J2000PosZ} + endPosECI := []float64{ + ae.J2000PosX + ae.J2000VelX*float64(ae.Microsecond)/10e6, + ae.J2000PosY + ae.J2000VelY*float64(ae.Microsecond)/10e6, + ae.J2000PosZ + ae.J2000VelZ*float64(ae.Microsecond)/10e6, + } // FIXME: GPS 拟合效果不佳 // x0 := float64(r.auxHeads[startPosInAux].TimeSec) + float64(r.auxHeads[startPosInAux].TimeSecFrac)/10e6