coordinate frames transformation

pkg/calculator/sofa.go

@@ -0,0 +1,54 @@
+package calculator
+
+import (
+	"math"
+	"time"
+
+	"github.com/hebl/gofa"
+)
+
+func ECItoECEF(j2000X, j2000Y, j2000Z float64, utc time.Time) (float64, float64, float64) {
+	j2000Position := [3]float64{j2000X, j2000Y, j2000Z}
+	var ecef [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 坐标
+	gofa.Rxp(rc2t, j2000Position, &ecef)
+	return ecef[0], ecef[1], ecef[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)
+	return lat * 180 / math.Pi, lon * 180 / math.Pi, alt
+}