own-tools/pkg/orbitviewer/coords.go

package orbitviewer

import (
	"math"
	"time"
)

const (
	EarthRadiusKm   = 6378.137 // WGS-84 equatorial radius in km
	EarthFlattening = 1.0 / 298.257223563
	Deg2Rad         = math.Pi / 180.0
	Rad2Deg         = 180.0 / math.Pi
)

// Vector3 represents a 3D vector [x, y, z] in km or km/s.
type Vector3 struct {
	X float64 `json:"x"`
	Y float64 `json:"y"`
	Z float64 `json:"z"`
}

// Geodetic represents latitude, longitude (in degrees) and altitude (in km).
type Geodetic struct {
	Latitude  float64 `json:"latitude"`  // Degrees (-90 to +90)
	Longitude float64 `json:"longitude"` // Degrees (-180 to +180)
	Altitude  float64 `json:"altitude"`  // Kilometers above WGS-84 ellipsoid
}

// LookAngles represents observer-relative topocentric angles.
type LookAngles struct {
	Azimuth   float64 `json:"azimuth"`   // Degrees (0 to 360, 0=North, 90=East)
	Elevation float64 `json:"elevation"` // Degrees (-90 to +90)
	Range     float64 `json:"range"`     // Distance in km
}

// JulianDate calculates the Julian Date for a given UTC time.
func JulianDate(t time.Time) float64 {
	t = t.UTC()
	year := t.Year()
	month := int(t.Month())
	day := t.Day()
	hour := t.Hour()
	min := t.Minute()
	sec := t.Second()
	nsec := t.Nanosecond()

	dayFraction := (float64(hour) + float64(min)/60.0 + (float64(sec)+float64(nsec)/1e9)/3600.0) / 24.0

	if month <= 2 {
		year--
		month += 12
	}

	a := math.Floor(float64(year) / 100.0)
	b := 2.0 - a + math.Floor(a/4.0)

	jd := math.Floor(365.25*float64(year+4716)) +
		math.Floor(30.6001*float64(month+1)) +
		float64(day) + dayFraction + b - 1524.5

	return jd
}

// GMST calculates Greenwich Mean Sidereal Time in radians for a given Julian Date.
func GMST(jd float64) float64 {
	d := jd - 2451545.0 // Days since J2000.0
	// Standard formula in degrees:
	gmstDeg := 280.46061837 + 360.98564736629*d
	// Normalize to 0 .. 360
	gmstDeg = math.Mod(gmstDeg, 360.0)
	if gmstDeg < 0 {
		gmstDeg += 360.0
	}
	return gmstDeg * Deg2Rad
}

// TEMEtoECEF converts True Equator Mean Equinox (TEME) coordinates to Earth-Centered Earth-Fixed (ECEF) coordinates.
func TEMEtoECEF(pos Vector3, jd float64) Vector3 {
	theta := GMST(jd)
	cosTheta := math.Cos(theta)
	sinTheta := math.Sin(theta)

	return Vector3{
		X: pos.X*cosTheta + pos.Y*sinTheta,
		Y: -pos.X*sinTheta + pos.Y*cosTheta,
		Z: pos.Z,
	}
}

// ECEFtoGeodetic converts ECEF coordinates (in km) to WGS-84 Geodetic (Lat, Lon, Alt).
func ECEFtoGeodetic(pos Vector3) Geodetic {
	x := pos.X
	y := pos.Y
	z := pos.Z

	a := EarthRadiusKm
	f := EarthFlattening
	b := a * (1.0 - f)
	e2 := (a*a - b*b) / (a * a)
	ep2 := (a*a - b*b) / (b * b)

	p := math.Sqrt(x*x + y*y)
	theta := math.Atan2(z*a, p*b)

	sinTheta := math.Sin(theta)
	cosTheta := math.Cos(theta)

	lat := math.Atan2(z+ep2*b*sinTheta*sinTheta*sinTheta, p-e2*a*cosTheta*cosTheta*cosTheta)
	lon := math.Atan2(y, x)

	sinLat := math.Sin(lat)
	n := a / math.Sqrt(1.0-e2*sinLat*sinLat)
	alt := p/math.Cos(lat) - n

	latDeg := lat * Rad2Deg
	lonDeg := lon * Rad2Deg
	// Normalize longitude to -180..180
	for lonDeg > 180 {
		lonDeg -= 360
	}
	for lonDeg < -180 {
		lonDeg += 360
	}

	return Geodetic{
		Latitude:  latDeg,
		Longitude: lonDeg,
		Altitude:  alt,
	}
}

// GeodeticToECEF converts WGS-84 Geodetic (lat, lon degrees, alt km) to ECEF (km).
func GeodeticToECEF(geo Geodetic) Vector3 {
	latRad := geo.Latitude * Deg2Rad
	lonRad := geo.Longitude * Deg2Rad

	a := EarthRadiusKm
	f := EarthFlattening
	e2 := 2*f - f*f

	sinLat := math.Sin(latRad)
	cosLat := math.Cos(latRad)
	sinLon := math.Sin(lonRad)
	cosLon := math.Cos(lonRad)

	n := a / math.Sqrt(1.0-e2*sinLat*sinLat)

	return Vector3{
		X: (n + geo.Altitude) * cosLat * cosLon,
		Y: (n + geo.Altitude) * cosLat * sinLon,
		Z: (n*(1.0-e2) + geo.Altitude) * sinLat,
	}
}

// CalculateLookAngles computes Topocentric Azimuth (deg), Elevation (deg), and Range (km)
// of a target ECEF position from an observer Geodetic location.
func CalculateLookAngles(obs Geodetic, targetECEF Vector3) LookAngles {
	obsECEF := GeodeticToECEF(obs)

	// Range vector in ECEF
	rx := targetECEF.X - obsECEF.X
	ry := targetECEF.Y - obsECEF.Y
	rz := targetECEF.Z - obsECEF.Z
	rangeDist := math.Sqrt(rx*rx + ry*ry + rz*rz)

	if rangeDist < 1e-6 {
		return LookAngles{Azimuth: 0, Elevation: 90, Range: 0}
	}

	latRad := obs.Latitude * Deg2Rad
	lonRad := obs.Longitude * Deg2Rad

	sinLat := math.Sin(latRad)
	cosLat := math.Cos(latRad)
	sinLon := math.Sin(lonRad)
	cosLon := math.Cos(lonRad)

	// Topocentric horizon coordinates (East, North, Up)
	east := -sinLon*rx + cosLon*ry
	north := -sinLat*cosLon*rx - sinLat*sinLon*ry + cosLat*rz
	up := cosLat*cosLon*rx + cosLat*sinLon*ry + sinLat*rz

	azRad := math.Atan2(east, north)
	azDeg := azRad * Rad2Deg
	if azDeg < 0 {
		azDeg += 360.0
	}

	elRad := math.Asin(up / rangeDist)
	elDeg := elRad * Rad2Deg

	return LookAngles{
		Azimuth:   azDeg,
		Elevation: elDeg,
		Range:     rangeDist,
	}
}

// FootprintRadiusKm returns the visual/radio horizon footprint radius on Earth's surface for an altitude in km.
func FootprintRadiusKm(altKm float64) float64 {
	if altKm <= 0 {
		return 0
	}
	r := EarthRadiusKm
	// Central angle theta: cos(theta) = r / (r + h)
	cosTheta := r / (r + altKm)
	if cosTheta > 1.0 {
		cosTheta = 1.0
	}
	theta := math.Acos(cosTheta)
	return r * theta // Arc distance on Earth surface
}

// CardinalDirection converts azimuth degrees to an 8 or 16-point cardinal compass string.
func CardinalDirection(az float64) string {
	az = math.Mod(az, 360.0)
	if az < 0 {
		az += 360.0
	}
	dirs := []string{
		"N", "NNE", "NE", "ENE", "E", "ESE", "SE", "SSE",
		"S", "SSW", "SW", "WSW", "W", "WNW", "NW", "NNW", "N",
	}
	idx := int(math.Round(az / 22.5))
	if idx >= len(dirs) {
		idx = 0
	}
	return dirs[idx]
}