121 lines
3.7 KiB
Go
121 lines
3.7 KiB
Go
package orbitviewer
|
|
|
|
import (
|
|
"math"
|
|
"time"
|
|
)
|
|
|
|
// SunPositionECEF returns the position vector of the Sun in ECEF coordinates (km) for a given Julian Date.
|
|
func SunPositionECEF(jd float64) Vector3 {
|
|
d := jd - 2451545.0 // days from J2000.0
|
|
|
|
// Mean anomaly of the Sun
|
|
g := (357.529 + 0.98560028*d) * Deg2Rad
|
|
// Mean longitude of the Sun
|
|
q := 280.459 + 0.98564736*d
|
|
// Geocentric apparent ecliptic longitude
|
|
l := (q + 1.915*math.Sin(g) + 0.020*math.Sin(2*g)) * Deg2Rad
|
|
// Obliquity of ecliptic
|
|
e := (23.439 - 0.00000036*d) * Deg2Rad
|
|
|
|
// Distance to Sun in AU
|
|
rAU := 1.00014 - 0.01671*math.Cos(g) - 0.00014*math.Cos(2*g)
|
|
rKm := rAU * 149597870.7 // km
|
|
|
|
// Sun position in Equatorial coordinates (ECI)
|
|
sinL := math.Sin(l)
|
|
cosL := math.Cos(l)
|
|
cosE := math.Cos(e)
|
|
sinE := math.Sin(e)
|
|
|
|
sunECI := Vector3{
|
|
X: rKm * cosL,
|
|
Y: rKm * cosE * sinL,
|
|
Z: rKm * sinE * sinL,
|
|
}
|
|
|
|
// Rotate ECI to ECEF by GMST
|
|
return TEMEtoECEF(sunECI, jd)
|
|
}
|
|
|
|
// ObserverSunElevation calculates the Sun's elevation angle (in degrees) for an observer location at time t.
|
|
func ObserverSunElevation(obs Geodetic, t time.Time) float64 {
|
|
jd := JulianDate(t)
|
|
sunECEF := SunPositionECEF(jd)
|
|
angles := CalculateLookAngles(obs, sunECEF)
|
|
return angles.Elevation
|
|
}
|
|
|
|
// IsSatelliteIlluminated determines if a satellite at satECEF is illuminated by the Sun (not eclipsed by Earth).
|
|
// Uses a cylindrical Earth shadow model with EarthRadiusKm.
|
|
func IsSatelliteIlluminated(satECEF Vector3, jd float64) bool {
|
|
sunECEF := SunPositionECEF(jd)
|
|
|
|
// Vector from Earth center to Sun (unit vector)
|
|
sunDist := math.Sqrt(sunECEF.X*sunECEF.X + sunECEF.Y*sunECEF.Y + sunECEF.Z*sunECEF.Z)
|
|
ux := sunECEF.X / sunDist
|
|
uy := sunECEF.Y / sunDist
|
|
uz := sunECEF.Z / sunDist
|
|
|
|
// Projection of satellite vector along Sun direction
|
|
dot := satECEF.X*ux + satECEF.Y*uy + satECEF.Z*uz
|
|
|
|
// If satellite is on day side of Earth (dot >= 0), it is illuminated
|
|
if dot >= 0 {
|
|
return true
|
|
}
|
|
|
|
// Satellite is on night side of Earth (dot < 0).
|
|
// Distance squared from shadow axis:
|
|
satDistSq := satECEF.X*satECEF.X + satECEF.Y*satECEF.Y + satECEF.Z*satECEF.Z
|
|
distFromAxisSq := satDistSq - dot*dot
|
|
|
|
// If distance from Earth shadow axis is greater than Earth radius, satellite is in sunlight
|
|
return distFromAxisSq > (EarthRadiusKm * EarthRadiusKm)
|
|
}
|
|
|
|
// EstimateVisualMagnitude estimates apparent visual magnitude for a satellite.
|
|
// standardMag: intrinsic magnitude at 1000 km distance and 50% illumination (e.g. -1.3 for ISS, ~5.0 for typical cube/small sat).
|
|
func EstimateVisualMagnitude(standardMag float64, satECEF, obsECEF Vector3, jd float64) float64 {
|
|
sunECEF := SunPositionECEF(jd)
|
|
|
|
// Satellite to Observer vector
|
|
dx := obsECEF.X - satECEF.X
|
|
dy := obsECEF.Y - satECEF.Y
|
|
dz := obsECEF.Z - satECEF.Z
|
|
rangeKm := math.Sqrt(dx*dx + dy*dy + dz*dz)
|
|
if rangeKm < 1.0 {
|
|
rangeKm = 1.0
|
|
}
|
|
|
|
// Satellite to Sun vector
|
|
sx := sunECEF.X - satECEF.X
|
|
sy := sunECEF.Y - satECEF.Y
|
|
sz := sunECEF.Z - satECEF.Z
|
|
sunDist := math.Sqrt(sx*sx + sy*sy + sz*sz)
|
|
|
|
// Cosine of phase angle (angle between Sun and Observer as seen from Satellite)
|
|
cosPhase := (dx*sx + dy*sy + dz*sz) / (rangeKm * sunDist)
|
|
if cosPhase > 1.0 {
|
|
cosPhase = 1.0
|
|
} else if cosPhase < -1.0 {
|
|
cosPhase = -1.0
|
|
}
|
|
phaseAngle := math.Acos(cosPhase)
|
|
|
|
// Phase function (diffuse sphere approximation):
|
|
// phi(beta) = (sin(beta) + (pi - beta)*cos(beta)) / pi
|
|
sinP := math.Sin(phaseAngle)
|
|
cosP := math.Cos(phaseAngle)
|
|
phi := (sinP + (math.Pi-phaseAngle)*cosP) / math.Pi
|
|
if phi <= 0.001 {
|
|
phi = 0.001
|
|
}
|
|
|
|
// Apparent magnitude:
|
|
// m = M_std - 15.75 + 5 * log10(range) - 2.5 * log10(phi)
|
|
// (Normalized so standardMag is roughly apparent mag at 1000km range and optimal phase)
|
|
mag := standardMag + 5.0*math.Log10(rangeKm/1000.0) - 2.5*math.Log10(phi)
|
|
return math.Round(mag*10.0) / 10.0
|
|
}
|