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
}