暂时使用星下点坐标作为图像左上角坐标

pkg/calculator/const.go

@@ -0,0 +1,10 @@
+package calculator
+
+// WGS84 ellipsoid constants
+const (
+	EarthRadius = 6378137.0         // 地球半径，单位米
+	a           = 6378137.0         // semi-major axis in meters
+	f           = 1 / 298.257223563 // flattening
+	e2          = 2*f - f*f         // square of eccentricity
+	J2000Epoch  = 2451545.0         // Julian date of J2000 epoch
+)

pkg/calculator/j2000.go

@@ -0,0 +1,78 @@
+package calculator
+
+import (
+	"math"
+	"time"
+)
+
+// Function to convert current J2000 position to WGS84
+func J2000ToWGS84(j2000X, j2000Y, j2000Z float64, utc time.Time) (float64, float64, float64) {
+	julianDate := UTCToJulianDate(utc)
+	gast := CalculateGAST(julianDate, utc)
+
+	// Calculate rotation matrix from J2000 to ITRS
+	rotationMatrix := CreateRotationMatrix(gast)
+
+	// Rotate J2000 position to ITRS
+	itrsX := rotationMatrix[0][0]*j2000X + rotationMatrix[0][1]*j2000Y + rotationMatrix[0][2]*j2000Z
+	itrsY := rotationMatrix[1][0]*j2000X + rotationMatrix[1][1]*j2000Y + rotationMatrix[1][2]*j2000Z
+	itrsZ := rotationMatrix[2][0]*j2000X + rotationMatrix[2][1]*j2000Y + rotationMatrix[2][2]*j2000Z
+
+	// Convert ITRS to geodetic coordinates (WGS84)
+	latitude, longitude, height := ECEFToGeodetic(itrsX, itrsY, itrsZ)
+	return latitude, longitude, height
+}
+
+// Function to convert UTC to Julian Date
+func UTCToJulianDate(t time.Time) float64 {
+	year, month, day := t.Date()
+	hour, minute, second := t.Clock()
+	fracOfDay := float64(hour)/24 + float64(minute)/(24*60) + float64(second)/(24*60*60)
+
+	if month <= 2 {
+		year -= 1
+		month += 12
+	}
+
+	A := year / 100
+	B := 2 - A + A/4
+	julianDate := float64(int(365.25*float64(year))) + float64(int(30.6001*float64(month+1))) + float64(day) + 1720994.5 + fracOfDay + float64(B)
+	return julianDate
+}
+
+// Function to calculate GAST (Greenwich Apparent Sidereal Time)
+func CalculateGAST(julianDate float64, utc time.Time) float64 {
+	// Calculate Julian centuries since J2000.0
+	T := (julianDate - J2000Epoch) / 36525.0
+
+	// Calculate Greenwich Mean Sidereal Time (GMST) in degrees
+	GMST := 280.46061837 + 360.98564736629*(julianDate-J2000Epoch) + 0.000387933*T*T - T*T*T/38710000
+	GMST = math.Mod(GMST, 360.0)
+	if GMST < 0 {
+		GMST += 360.0
+	}
+
+	// Calculate the equation of the equinoxes (simplified)
+	omega := 125.04 - 0.052954*T
+	L := 280.47 + 0.98565*T
+	epsilon := 23.4393 - 0.0000004*T
+	deltaPsi := -0.000319*math.Sin(math.Pi/180.0*omega) - 0.000024*math.Sin(math.Pi/180.0*2*L)
+
+	// Greenwich Apparent Sidereal Time (GAST)
+	GAST := GMST + deltaPsi*math.Cos(math.Pi/180.0*epsilon)
+	GAST = math.Mod(GAST, 360.0)
+	if GAST < 0 {
+		GAST += 360.0
+	}
+
+	return GAST * math.Pi / 180.0 // Convert to radians
+}
+
+// Function to create the rotation matrix for the given GAST
+func CreateRotationMatrix(gast float64) [3][3]float64 {
+	return [3][3]float64{
+		{math.Cos(gast), math.Sin(gast), 0},
+		{-math.Sin(gast), math.Cos(gast), 0},
+		{0, 0, 1},
+	}
+}

pkg/calculator/quaternion.go

@@ -0,0 +1,70 @@
+package calculator
+
+import (
+	"fmt"
+	"math"
+)
+
+// Quaternion represents a quaternion with scalar (w) and vector (x, y, z) parts
+type Quaternion struct {
+	w, x, y, z float64
+}
+
+// Quaternion multiplication
+func (q1 Quaternion) Mul(q2 Quaternion) Quaternion {
+	return Quaternion{
+		w: q1.w*q2.w - q1.x*q2.x - q1.y*q2.y - q1.z*q2.z,
+		x: q1.w*q2.x + q1.x*q2.w + q1.y*q2.z - q1.z*q2.y,
+		y: q1.w*q2.y - q1.x*q2.z + q1.y*q2.w + q1.z*q2.x,
+		z: q1.w*q2.z + q1.x*q2.y - q1.y*q2.x + q1.z*q2.w,
+	}
+}
+
+// Quaternion conjugate
+func (q Quaternion) Conjugate() Quaternion {
+	return Quaternion{w: q.w, x: -q.x, y: -q.y, z: -q.z}
+}
+
+// Rotate vector by quaternion
+func (q Quaternion) Rotate(v [3]float64) [3]float64 {
+	qv := Quaternion{w: 0, x: v[0], y: v[1], z: v[2]}
+	qConj := q.Conjugate()
+	qvRotated := q.Mul(qv).Mul(qConj)
+	return [3]float64{qvRotated.x, qvRotated.y, qvRotated.z}
+}
+
+func main() {
+	// 示例数据
+	qBI := Quaternion{w: 1, x: 0, y: 0, z: 0} // 本体相对惯性系四元数
+	posJ2000 := [3]float64{7000, 0, 0}        // J2000位置
+	// velJ2000 := [3]float64{0, 7.5, 0}         // J2000速度
+
+	// 相机参数
+	const numPixels = 9520
+	const fov = 10.0 * math.Pi / 180 // 假设视场角为10度
+
+	// 逐像素计算地面交点
+	for i := 0; i < numPixels; i++ {
+		// 计算像素点相对光轴的偏角
+		pixelOffset := (float64(i) - float64(numPixels)/2) / float64(numPixels)
+		angle := pixelOffset * fov
+
+		// 假设光轴在本体坐标系中指向-z方向，计算视线方向
+		dBody := [3]float64{-math.Sin(angle), 0, -math.Cos(angle)}
+
+		// 转换到惯性系
+		dInertial := qBI.Rotate(dBody)
+
+		// 计算地面交点（假设dInertial已经标准化）
+		k := -posJ2000[2] / dInertial[2] // 简化的交点计算
+		groundPoint := [3]float64{
+			posJ2000[0] + k*dInertial[0],
+			posJ2000[1] + k*dInertial[1],
+			posJ2000[2] + k*dInertial[2],
+		}
+
+		// 转换到地理坐标
+		lat, lon, _ := ECEFToGeodetic(groundPoint[0], groundPoint[1], groundPoint[2])
+		fmt.Printf("Pixel %d: Latitude: %f, Longitude: %f\n", i, lat, lon)
+	}
+}

pkg/calculator/wgs84.go

@@ -0,0 +1,77 @@
+package calculator
+
+import "math"
+
+func WGS84XYZtoLatLngH(X, Y, Z float64) (float64, float64, float64) {
+	return ECEFToGeodetic(X, Y, Z)
+}
+
+// Function to convert ECEF (ITRS) coordinates to geodetic coordinates (latitude, longitude, height)
+func ECEFToGeodetic(X, Y, Z float64) (float64, float64, float64) {
+	b := a * (1 - f)
+	e2Prime := e2 * (a * a) / (b * b)
+	p := math.Sqrt(X*X + Y*Y)
+	theta := math.Atan2(Z*a, p*b)
+
+	// Calculate Longitude
+	longitude := math.Atan2(Y, X)
+
+	// Calculate Latitude
+	latitude := math.Atan2(Z+e2Prime*b*math.Pow(math.Sin(theta), 3), p-e2*a*math.Pow(math.Cos(theta), 3))
+
+	// Calculate Height
+	N := a / math.Sqrt(1-e2*math.Pow(math.Sin(latitude), 2))
+	height := p/math.Cos(latitude) - N
+
+	// Convert radians to degrees for latitude and longitude
+	latitudeDeg := latitude * 180 / math.Pi
+	longitudeDeg := longitude * 180 / math.Pi
+
+	return latitudeDeg, longitudeDeg, height
+}
+
+// Converts degrees to radians.
+func degreesToRadians(degrees float64) float64 {
+	return degrees * math.Pi / 180
+}
+
+// Converts radians to degrees.
+func radiansToDegrees(radians float64) float64 {
+	return radians * 180 / math.Pi
+}
+
+// Calculate destination point given a start point (lat1, lng1), distance dx, dy (in meters).
+func CalculateDestination(lat1, lng1, dx, dy float64) (float64, float64) {
+	// Convert latitude and longitude from degrees to radians.
+	lat1Rad := degreesToRadians(lat1)
+	lng1Rad := degreesToRadians(lng1)
+
+	// Radius of the Earth at the given latitude.
+	radius := EarthRadius * math.Cos(lat1Rad)
+
+	// Calculate the new latitude in radians.
+	newLatRad := lat1Rad + dy/EarthRadius
+
+	// Calculate the new longitude in radians.
+	newLngRad := lng1Rad + dx/radius
+
+	// Convert the new latitude and longitude back to degrees.
+	newLat := radiansToDegrees(newLatRad)
+	newLng := radiansToDegrees(newLngRad)
+
+	return newLat, newLng
+}
+
+const (
+	MetersPerDegreeLatitude = 111320.0 // 1度纬度约等于111,320米
+)
+
+// Converts a distance in meters to degrees longitude at a specific latitude
+func MetersToDegreesLongitude(meters float64, latitude float64) float64 {
+	return meters / (111320.0 * cosLatitude(latitude))
+}
+
+// Approximates the cosine of the latitude in radians
+func cosLatitude(latitude float64) float64 {
+	return 1.0 / (1.0 + (latitude/90.0)*(latitude/90.0))
+}