package orbitviewer

import (
	"errors"
	"fmt"
	"math"
	"strconv"
	"strings"
	"time"
)

// WGS-72 constants used in standard SGP4/SDP4
const (
	sgp4Ck2    = 5.413080e-4   // 0.5 * J2
	sgp4Ck4    = 0.62098875e-6 // -0.375 * J4
	sgp4E6a    = 1.0e-6
	sgp4Qo     = 120.0        // Parameter for atmospheric model
	sgp4So     = 78.0         // Parameter for atmospheric model
	sgp4Xke    = 0.0743669161 // sqrt(GM) in (Earth radii)^1.5 / min
	sgp4Xkmper = 6378.135     // Earth equatorial radius in km (WGS-72)
	sgp4Xj2    = 1.082616e-3  // J2 harmonic
	sgp4Xj3    = -0.253881e-5 // J3 harmonic
	sgp4Xj4    = -1.65597e-6  // J4 harmonic
	sgp4TwoPi  = 2.0 * math.Pi
)

// TLE represents a decoded Two-Line Element set.
type TLE struct {
	Name             string    `json:"name"`
	Line1            string    `json:"line1"`
	Line2            string    `json:"line2"`
	NoradID          int       `json:"norad_id"`
	Classification   string    `json:"classification"`
	IntlDesignator   string    `json:"intl_designator"`
	EpochYear        int       `json:"epoch_year"`
	EpochDay         float64   `json:"epoch_day"`
	EpochTime        time.Time `json:"epoch_time"`
	BStar            float64   `json:"bstar"`
	InclinationDeg   float64   `json:"inclination_deg"`
	RaanDeg          float64   `json:"raan_deg"`
	Eccentricity     float64   `json:"eccentricity"`
	ArgPerigeeDeg    float64   `json:"arg_perigee_deg"`
	MeanAnomalyDeg   float64   `json:"mean_anomaly_deg"`
	MeanMotionRevDay float64   `json:"mean_motion_rev_day"` // revs per day
	RevNumber        int       `json:"rev_number"`
}

// SGP4Model contains initialized constants and precomputed values for orbit propagation.
type SGP4Model struct {
	TLE TLE

	// Orbital elements at epoch
	a0, e0, i0, om0, w0, m0 float64
	n0                      float64 // rad / min
	bstar                   float64

	// Precomputed constants
	cosi0, sini0          float64
	theta2, theta4        float64
	beta0, beta02, beta04 float64
	ao, no0               float64
	isImpSimple           bool

	// Secular rate terms
	xmdot, omgdot, xnodot      float64
	c1, c2, c3, c4, c5         float64
	d2, d3, d4                 float64
	t2cof, t3cof, t4cof, t5cof float64
	eta                        float64
	omgcof, xmcof              float64
	delmo, sinmao              float64
}

// ParseTLE parses two (or three) lines of standard NORAD TLE text into a TLE struct.
func ParseTLE(name, line1, line2 string) (TLE, error) {
	line1 = strings.TrimSpace(line1)
	line2 = strings.TrimSpace(line2)

	if len(line1) < 68 || len(line2) < 68 {
		return TLE{}, errors.New("invalid TLE: line length must be at least 68 characters")
	}

	if line1[0] != '1' || line2[0] != '2' {
		return TLE{}, errors.New("invalid TLE: line 1 must start with '1' and line 2 with '2'")
	}

	noradID1, err := strconv.Atoi(strings.TrimSpace(line1[2:7]))
	if err != nil {
		return TLE{}, fmt.Errorf("invalid NORAD ID in line 1: %w", err)
	}
	noradID2, err := strconv.Atoi(strings.TrimSpace(line2[2:7]))
	if err != nil || noradID1 != noradID2 {
		return TLE{}, fmt.Errorf("NORAD ID mismatch between line 1 (%d) and line 2 (%d)", noradID1, noradID2)
	}

	class := string(line1[7])
	intlDesig := strings.TrimSpace(line1[9:17])

	epochYr2, err := strconv.Atoi(strings.TrimSpace(line1[18:20]))
	if err != nil {
		return TLE{}, fmt.Errorf("invalid epoch year: %w", err)
	}
	var fullYear int
	if epochYr2 < 57 {
		fullYear = 2000 + epochYr2
	} else {
		fullYear = 1900 + epochYr2
	}

	epochDay, err := strconv.ParseFloat(strings.TrimSpace(line1[20:32]), 64)
	if err != nil {
		return TLE{}, fmt.Errorf("invalid epoch day: %w", err)
	}

	// Calculate epoch timestamp
	startOfYear := time.Date(fullYear, 1, 1, 0, 0, 0, 0, time.UTC)
	epochDuration := time.Duration(epochDay * float64(24*time.Hour))
	epochTime := startOfYear.Add(epochDuration - 24*time.Hour)

	// BSTAR drag term: e.g. " 12345-4" -> 0.12345e-4
	bstarStr := strings.TrimSpace(line1[53:61])
	bstar, err := parseExponential(bstarStr)
	if err != nil {
		bstar = 0.0
	}

	// Line 2 elements
	inc, err := strconv.ParseFloat(strings.TrimSpace(line2[8:16]), 64)
	if err != nil {
		return TLE{}, fmt.Errorf("invalid inclination: %w", err)
	}

	raan, err := strconv.ParseFloat(strings.TrimSpace(line2[17:25]), 64)
	if err != nil {
		return TLE{}, fmt.Errorf("invalid RAAN: %w", err)
	}

	eccStr := "0." + strings.TrimSpace(line2[26:33])
	ecc, err := strconv.ParseFloat(eccStr, 64)
	if err != nil {
		return TLE{}, fmt.Errorf("invalid eccentricity: %w", err)
	}

	argPerigee, err := strconv.ParseFloat(strings.TrimSpace(line2[34:42]), 64)
	if err != nil {
		return TLE{}, fmt.Errorf("invalid arg of perigee: %w", err)
	}

	meanAnomaly, err := strconv.ParseFloat(strings.TrimSpace(line2[43:51]), 64)
	if err != nil {
		return TLE{}, fmt.Errorf("invalid mean anomaly: %w", err)
	}

	meanMotion, err := strconv.ParseFloat(strings.TrimSpace(line2[52:63]), 64)
	if err != nil {
		return TLE{}, fmt.Errorf("invalid mean motion: %w", err)
	}

	revNum, _ := strconv.Atoi(strings.TrimSpace(line2[63:68]))

	if name == "" {
		name = fmt.Sprintf("SAT-%d", noradID1)
	}

	return TLE{
		Name:             name,
		Line1:            line1,
		Line2:            line2,
		NoradID:          noradID1,
		Classification:   class,
		IntlDesignator:   intlDesig,
		EpochYear:        fullYear,
		EpochDay:         epochDay,
		EpochTime:        epochTime,
		BStar:            bstar,
		InclinationDeg:   inc,
		RaanDeg:          raan,
		Eccentricity:     ecc,
		ArgPerigeeDeg:    argPerigee,
		MeanAnomalyDeg:   meanAnomaly,
		MeanMotionRevDay: meanMotion,
		RevNumber:        revNum,
	}, nil
}

func parseExponential(s string) (float64, error) {
	s = strings.TrimSpace(s)
	if s == "" || s == "00000-0" || s == "00000+0" {
		return 0, nil
	}
	sign := 1.0
	switch s[0] {
	case '-':
		sign = -1.0
		s = s[1:]
	case '+':
		s = s[1:]
	}
	// Format is typically 5 digits followed by sign and exponent: e.g. "12345-4"
	n := len(s)
	if n >= 2 && (s[n-2] == '-' || s[n-2] == '+') {
		mantissaStr := "0." + s[:n-2]
		expStr := s[n-2:]
		m, err := strconv.ParseFloat(mantissaStr, 64)
		if err != nil {
			return 0, err
		}
		e, err := strconv.Atoi(expStr)
		if err != nil {
			return 0, err
		}
		return sign * m * math.Pow10(e), nil
	}
	return strconv.ParseFloat(s, 64)
}

// NewSGP4Model initializes the SGP4 propagation model from a TLE.
func NewSGP4Model(tle TLE) (*SGP4Model, error) {
	m := &SGP4Model{
		TLE:   tle,
		i0:    tle.InclinationDeg * Deg2Rad,
		om0:   tle.RaanDeg * Deg2Rad,
		e0:    tle.Eccentricity,
		w0:    tle.ArgPerigeeDeg * Deg2Rad,
		m0:    tle.MeanAnomalyDeg * Deg2Rad,
		n0:    tle.MeanMotionRevDay * sgp4TwoPi / 1440.0, // rad/min
		bstar: tle.BStar,
	}

	if m.e0 < 0 || m.e0 >= 0.999 {
		return nil, fmt.Errorf("eccentricity %.4f out of range for SGP4", m.e0)
	}

	// Constants
	m.cosi0 = math.Cos(m.i0)
	m.sini0 = math.Sin(m.i0)
	m.theta2 = m.cosi0 * m.cosi0
	m.theta4 = m.theta2 * m.theta2
	m.beta02 = 1.0 - m.e0*m.e0
	m.beta0 = math.Sqrt(m.beta02)
	m.beta04 = m.beta02 * m.beta02

	// Recover original mean motion (no0) and semi-major axis (ao)
	a1 := math.Pow(sgp4Xke/m.n0, 2.0/3.0)
	del1 := 1.5 * sgp4Ck2 * (3.0*m.theta2 - 1.0) / (m.beta0 * m.beta02 * a1 * a1)
	ao := a1 * (1.0 - del1*(1.0/3.0+del1*(1.0+del1*134.0/81.0)))
	delo := 1.5 * sgp4Ck2 * (3.0*m.theta2 - 1.0) / (m.beta0 * m.beta02 * ao * ao)
	m.no0 = m.n0 / (1.0 + delo)
	m.ao = ao / (1.0 - delo)
	m.a0 = m.ao

	// Perigee altitude
	perigeeKm := (m.ao*(1.0-m.e0) - 1.0) * sgp4Xkmper
	m.isImpSimple = false
	if perigeeKm < 220.0 {
		m.isImpSimple = true
	}

	s4 := sgp4So/sgp4Xkmper + 1.0
	qoms24 := math.Pow((sgp4Qo-sgp4So)/sgp4Xkmper, 4.0)
	if perigeeKm < 156.0 {
		s4 = perigeeKm - 78.0
		if perigeeKm <= 98.0 {
			s4 = 20.0
		}
		qoms24 = math.Pow((sgp4Qo-s4)/sgp4Xkmper, 4.0)
		s4 = s4/sgp4Xkmper + 1.0
	}

	pinv := 1.0 / (m.ao * m.ao * m.beta04)
	tsi := 1.0 / (m.ao - s4)
	m.eta = m.ao * m.e0 * tsi
	etasq := m.eta * m.eta
	eeta := m.e0 * m.eta
	psisq := math.Abs(1.0 - etasq)
	coef := qoms24 * math.Pow(tsi, 4.0)
	coef1 := coef / math.Pow(psisq, 3.5)

	m.c2 = coef1 * m.no0 * (m.ao*(1.0+1.5*etasq+eeta*(4.0+etasq)) +
		0.75*sgp4Ck2*tsi/psisq*(8.0+3.0*etasq*(8.0+etasq))*(3.0*m.theta2-1.0))
	m.c1 = m.bstar * m.c2
	m.c3 = coef * tsi * (-2.0 * sgp4Xj3 / sgp4Ck2) * m.no0 * m.sini0 / m.e0
	if m.e0 < 1e-4 {
		m.c3 = 0.0
	}
	m.c4 = 2.0 * m.no0 * coef1 * m.ao * m.beta02 * (m.eta*(2.0+0.5*etasq) +
		m.e0*(0.5+2.0*etasq) -
		2.0*sgp4Ck2*tsi/(m.ao*psisq)*(3.0*(1.0-3.0*m.theta2)*(1.0-2.0*eeta+etasq*(1.5-0.5*eeta))+
			0.75*(1.0-m.theta2)*(2.0*etasq-eeta*(1.0+etasq))*math.Cos(2.0*m.w0)))
	m.c5 = 2.0 * coef1 * m.ao * m.beta02 * (1.0 + 2.75*(etasq+eeta) + eeta*etasq)

	// Secular rates
	temp1 := 3.0 * sgp4Ck2 * pinv * m.no0
	temp2 := temp1 * sgp4Ck2 * pinv
	temp3 := 1.25 * sgp4Ck4 * pinv * pinv * m.no0

	m.xmdot = m.no0 + 0.5*temp1*m.beta0*(3.0*m.theta2-1.0) +
		0.0625*temp2*m.beta0*(13.0-78.0*m.theta2+137.0*m.theta4)

	m.omgdot = -0.5*temp1*(1.0-5.0*m.theta2) +
		0.0625*temp2*(7.0-114.0*m.theta2+395.0*m.theta4) +
		temp3*(3.0-36.0*m.theta2+49.0*m.theta4)

	m.xnodot = -temp1*m.cosi0 +
		(0.5*temp2*(4.0-19.0*m.theta2)+2.0*temp3*(3.0-7.0*m.theta2))*m.cosi0

	if !m.isImpSimple {
		c1sq := m.c1 * m.c1
		m.d2 = 4.0 * m.ao * tsi * c1sq
		temp := m.d2 * tsi * m.c1 / 3.0
		m.d3 = (17.0*m.ao + s4) * temp
		m.d4 = 0.5 * temp * m.ao * tsi * (221.0*m.ao + 31.0*s4) * m.c1
		m.t3cof = m.d2 + 2.0*c1sq
		m.t4cof = 0.25 * (3.0*m.d3 + m.c1*(12.0*m.d2+10.0*c1sq))
		m.t5cof = 0.2 * (3.0*m.d4 + 12.0*m.c1*m.d3 + 6.0*m.d2*m.d2 + 15.0*c1sq*(2.0*m.d2+c1sq))

		m.omgcof = m.bstar * m.c3 * math.Cos(m.w0)
		m.xmcof = 0.0
		if m.e0 > 1e-4 && eeta > 0 {
			m.xmcof = -2.0 / 3.0 * coef * m.bstar / eeta
		}
		m.delmo = math.Pow(1.0+m.eta*math.Cos(m.m0), 3.0)
		m.sinmao = math.Sin(m.m0)
	}

	return m, nil
}

// PropagateMinutes calculates position (km) and velocity (km/s) in TEME at tsince minutes from epoch.
func (m *SGP4Model) PropagateMinutes(tsince float64) (Vector3, Vector3, error) {
	// Update for secular gravity and atmospheric drag
	xmdf := m.m0 + m.xmdot*tsince
	omgadf := m.w0 + m.omgdot*tsince
	xnoddf := m.om0 + m.xnodot*tsince
	omega := omgadf
	xmp := xmdf
	tsq := tsince * tsince
	xnode := xnoddf + 0.5*m.c1*tsq
	tempa := 1.0 - m.c1*tsince
	tempe := m.bstar * m.c4 * tsince
	templ := m.t2cof * tsq

	if !m.isImpSimple {
		delomg := m.omgcof * tsince
		delm := m.xmcof * (math.Pow(1.0+m.eta*math.Cos(xmdf), 3) - m.delmo)
		temp := delomg + delm
		xmp = xmdf + temp
		omega = omgadf - temp
		tcube := tsq * tsince
		tfour := tsq * tsq
		tempa = tempa - m.d2*tsq - m.d3*tcube - m.d4*tfour
		tempe = tempe + m.bstar*m.c5*(math.Sin(xmp)-m.sinmao)
		templ = templ + m.t3cof*tcube + tfour*(m.t4cof+tsince*m.t5cof)
	}

	a := m.ao * tempa * tempa
	e := m.e0 - tempe
	xl := xmp + omega + xnode + m.no0*templ
	beta := math.Sqrt(1.0 - e*e)
	n := sgp4Xke / math.Pow(a, 1.5)

	// Long period periodic terms
	axn := e * math.Cos(omega)
	temp := 1.0 / (a * beta * beta)
	xll := temp * 0.5 * sgp4Ck2 * 1.5 * (3.0*m.theta2 - 1.0)
	ayn := e*math.Sin(omega) + temp*(-0.5*sgp4Xj3/sgp4Ck2*m.sini0)
	xlt := xl - temp*0.25*(-0.5*sgp4Xj3/sgp4Ck2*m.sini0)*axn*(3.0+5.0*m.cosi0)/(1.0+m.cosi0)

	// Solve Kepler's Equation for (U = E + w)
	u := math.Mod(xlt-xnode, sgp4TwoPi)
	if u < 0 {
		u += sgp4TwoPi
	}
	eo1 := u
	tem5 := 1.0
	for iter := 0; iter < 10 && math.Abs(tem5) >= sgp4E6a; iter++ {
		sineo1 := math.Sin(eo1)
		coseo1 := math.Cos(eo1)
		tem5 = 1.0 - coseo1*axn - sineo1*ayn
		tem5 = (u - ayn*coseo1 + axn*sineo1 - eo1) / tem5
		if math.Abs(tem5) > 1.0 {
			if tem5 > 0 {
				tem5 = 1.0
			} else {
				tem5 = -1.0
			}
		}
		eo1 += tem5
	}

	// Short period preliminary quantities
	sineo1 := math.Sin(eo1)
	coseo1 := math.Cos(eo1)
	ecose := axn*coseo1 + ayn*sineo1
	esine := axn*sineo1 - ayn*coseo1
	el2 := axn*axn + ayn*ayn
	pl := a * (1.0 - el2)
	if pl < 0 {
		return Vector3{}, Vector3{}, errors.New("satellite decayed into Earth")
	}

	r := a * (1.0 - ecose)
	rdot := sgp4Xke * math.Sqrt(a) * esine / r
	rfdot := sgp4Xke * math.Sqrt(pl) / r
	betal := math.Sqrt(1.0 - el2)
	temp = esine / (1.0 + betal)
	sinu := a / r * (sineo1 - ayn - axn*temp)
	cosu := a / r * (coseo1 - axn + ayn*temp)
	su := math.Atan2(sinu, cosu)
	sin2u := 2.0 * sinu * cosu
	cos2u := 2.0*cosu*cosu - 1.0
	temp = 1.0 / pl
	temp1 := sgp4Ck2 * temp
	temp2 := temp1 * temp

	// Update for short periodics
	rk := r*(1.0-1.5*temp2*betal*(3.0*m.theta2-1.0)) + 0.5*temp1*(1.0-m.theta2)*cos2u
	uk := su - 0.25*temp2*(7.0*m.theta2-1.0)*sin2u
	xnodek := xnode + 1.5*temp2*m.cosi0*sin2u
	xinck := m.i0 + 1.5*temp2*m.cosi0*m.sini0*cos2u
	rdotk := rdot - n*temp1*(1.0-m.theta2)*sin2u
	rfdotk := rfdot + n*temp1*((1.0-m.theta2)*cos2u+1.5*(3.0*m.theta2-1.0))

	// Orientation vectors in TEME
	sinuk := math.Sin(uk)
	cosuk := math.Cos(uk)
	sinik := math.Sin(xinck)
	cosik := math.Cos(xinck)
	sinnok := math.Sin(xnodek)
	cosnok := math.Cos(xnodek)

	mx := -sinnok * cosik
	my := cosnok * cosik
	ux := mx*sinuk + cosnok*cosuk
	uy := my*sinuk + sinnok*cosuk
	uz := sinik * sinuk
	vx := mx*cosuk - cosnok*sinuk
	vy := my*cosuk - sinnok*sinuk
	vz := sinik * cosuk

	// Position and velocity in km and km/s
	pos := Vector3{
		X: rk * ux * sgp4Xkmper,
		Y: rk * uy * sgp4Xkmper,
		Z: rk * uz * sgp4Xkmper,
	}

	vel := Vector3{
		X: (rdotk*ux + rfdotk*vx) * sgp4Xkmper / 60.0,
		Y: (rdotk*uy + rfdotk*vy) * sgp4Xkmper / 60.0,
		Z: (rdotk*uz + rfdotk*vz) * sgp4Xkmper / 60.0,
	}

	_ = xll
	return pos, vel, nil
}

// Additional helper field values
var (
	_ = sgp4Xj4
)

// PropagateTime calculates position (km) and velocity (km/s) in TEME at an absolute UTC time.
func (m *SGP4Model) PropagateTime(t time.Time) (Vector3, Vector3, error) {
	tsince := t.Sub(m.TLE.EpochTime).Minutes()
	return m.PropagateMinutes(tsince)
}

// PositionECEF calculates ECEF coordinates at time t.
func (m *SGP4Model) PositionECEF(t time.Time) (Vector3, Vector3, error) {
	posTEME, velTEME, err := m.PropagateTime(t)
	if err != nil {
		return Vector3{}, Vector3{}, err
	}
	jd := JulianDate(t)
	posECEF := TEMEtoECEF(posTEME, jd)
	// Approximate velocity in ECEF
	velECEF := TEMEtoECEF(velTEME, jd)
	return posECEF, velECEF, nil
}

// GeodeticPosition calculates WGS-84 Geodetic position (lat, lon, alt) at time t.
func (m *SGP4Model) GeodeticPosition(t time.Time) (Geodetic, Vector3, error) {
	posECEF, velECEF, err := m.PositionECEF(t)
	if err != nil {
		return Geodetic{}, Vector3{}, err
	}
	geo := ECEFtoGeodetic(posECEF)
	return geo, velECEF, nil
}