Added tweaks to positions and passes calculation

2026-10-02 03:15:37 +00:00 · 2026-03-29 13:09:03 +01:00
1 parent 279594c425
commit 0bef2ffddc
@@ -97,40 +97,38 @@ class DeepSpaceObject(data: OrbitalData) : OrbitalObject(data) {
    }

    internal fun calculateSDP4(tSince: Double) {
-        synchronized(this) {
-            val temp = DoubleArray(12)
-            val xmdf = data.xmo + dsv.xmdot * tSince
-            val tsq = tSince * tSince
-            val templ = t2cof * tsq
-            dsv.xll = xmdf + dsv.xnodp * templ
-            dsv.omgadf = data.omegao + dsv.omgdot * tSince
-            val xnoddf = data.xnodeo + dsv.xnodot * tSince
-            dsv.xnode = xnoddf + xnodcf * tsq
-            val tempa = 1.0 - c1 * tSince
-            val tempe = data.bstar * c4 * tSince
-            dsv.xn = dsv.xnodp
-            dsv.t = tSince
-            deep.dpsec(data)
-            val a = (XKE / dsv.xn).pow(TWO_THIRDS) * tempa * tempa
-            dsv.em -= tempe
-            deep.dpper()
-            val xl = dsv.xll + dsv.omgadf + dsv.xnode
-            val beta = sqrt(1.0 - dsv.em * dsv.em)
-            dsv.xn = XKE / a.pow(1.5)
-            // Long period periodics
-            val axn = dsv.em * cos(dsv.omgadf)
-            temp[0] = invert(a * beta * beta)
-            val xll = temp[0] * xlcof * axn
-            val aynl = temp[0] * aycof
-            val xlt = xl + xll
-            val ayn = dsv.em * sin(dsv.omgadf) + aynl
-            // Solve Kepler's equation
-            val capu = mod2PI(xlt - dsv.xnode)
-            temp[2] = capu
-            converge(temp, axn, ayn, capu)
-            calculatePosAndVel(temp, a, axn, ayn)
-            calculatePhase(xlt, dsv.xnode, dsv.omgadf)
-        }
+        val temp = DoubleArray(12)
+        val xmdf = data.xmo + dsv.xmdot * tSince
+        val tsq = tSince * tSince
+        val templ = t2cof * tsq
+        dsv.xll = xmdf + dsv.xnodp * templ
+        dsv.omgadf = data.omegao + dsv.omgdot * tSince
+        val xnoddf = data.xnodeo + dsv.xnodot * tSince
+        dsv.xnode = xnoddf + xnodcf * tsq
+        val tempa = 1.0 - c1 * tSince
+        val tempe = data.bstar * c4 * tSince
+        dsv.xn = dsv.xnodp
+        dsv.t = tSince
+        deep.dpsec(data)
+        val a = (XKE / dsv.xn).pow(TWO_THIRDS) * tempa * tempa
+        dsv.em -= tempe
+        deep.dpper()
+        val xl = dsv.xll + dsv.omgadf + dsv.xnode
+        val beta = sqrt(1.0 - dsv.em * dsv.em)
+        dsv.xn = XKE / a.pow(1.5)
+        // Long period periodics
+        val axn = dsv.em * cos(dsv.omgadf)
+        temp[0] = invert(a * beta * beta)
+        val xll = temp[0] * xlcof * axn
+        val aynl = temp[0] * aycof
+        val xlt = xl + xll
+        val ayn = dsv.em * sin(dsv.omgadf) + aynl
+        // Solve Kepler's equation
+        val capu = mod2PI(xlt - dsv.xnode)
+        temp[2] = capu
+        converge(temp, axn, ayn, capu)
+        calculatePosAndVel(temp, a, axn, ayn)
+        calculatePhase(xlt, dsv.xnode, dsv.omgadf)
    }

    private fun calculatePosAndVel(temp: DoubleArray, a: Double, axn: Double, ayn: Double) {
@@ -141,53 +141,51 @@ class NearEarthObject(data: OrbitalData) : OrbitalObject(data) {
    }

    internal fun calculateSGP4(tSince: Double) {
-        synchronized(this) {
-            val temp = DoubleArray(9)
-            val xmdf = data.xmo + xmdot * tSince
-            val omgadf = data.omegao + omgdot * tSince
-            val xnoddf = data.xnodeo + xnodot * tSince
-            var omega = omgadf
-            var xmp = xmdf
-            val tsq = sqr(tSince)
-            val xnode = xnoddf + xnodcf * tsq
-            val bstar = data.bstar
-            var tempa = 1.0 - c1 * tSince
-            var tempe = bstar * c4 * tSince
-            var templ = t2cof * tsq
-            if (!sgp4Simple) {
-                val delomg = omgcof * tSince
-                val delm = xmcof * ((1.0 + eta * cos(xmdf)).pow(3.0) - delmo)
-                temp[0] = delomg + delm
-                xmp = xmdf + temp[0]
-                omega = omgadf - temp[0]
-                val tcube = tsq * tSince
-                val tfour = tSince * tcube
-                tempa = tempa - d2 * tsq - d3 * tcube - d4 * tfour
-                tempe += bstar * c5 * (sin(xmp) - sinmo)
-                templ += t3cof * tcube + tfour * (t4cof + tSince * t5cof)
-            }
-            val a = aodp * tempa.pow(2.0)
-            val eo = data.eccn
-            val e = eo - tempe
-            val xl = xmp + omega + xnode + xnodp * templ
-            val beta = sqrt(1.0 - e * e)
-            val xn = XKE / a.pow(1.5)
-
-            // Long period periodics
-            val axn = e * cos(omega)
-            temp[0] = invert(a * sqr(beta))
-            val xll = temp[0] * xlcof * axn
-            val aynl = temp[0] * aycof
-            val xlt = xl + xll
-            val ayn = e * sin(omega) + aynl
-
-            // Solve Kepler's equation
-            val capu = mod2PI(xlt - xnode)
-            temp[2] = capu
-            converge(temp, axn, ayn, capu)
-            calculatePosAndVel(temp, xnode, a, xn, axn, ayn)
-            calculatePhase(xlt, xnode, omgadf)
+        val temp = DoubleArray(9)
+        val xmdf = data.xmo + xmdot * tSince
+        val omgadf = data.omegao + omgdot * tSince
+        val xnoddf = data.xnodeo + xnodot * tSince
+        var omega = omgadf
+        var xmp = xmdf
+        val tsq = tSince * tSince
+        val xnode = xnoddf + xnodcf * tsq
+        val bstar = data.bstar
+        var tempa = 1.0 - c1 * tSince
+        var tempe = bstar * c4 * tSince
+        var templ = t2cof * tsq
+        if (!sgp4Simple) {
+            val delomg = omgcof * tSince
+            val delm = xmcof * ((1.0 + eta * cos(xmdf)).pow(3.0) - delmo)
+            temp[0] = delomg + delm
+            xmp = xmdf + temp[0]
+            omega = omgadf - temp[0]
+            val tcube = tsq * tSince
+            val tfour = tSince * tcube
+            tempa = tempa - d2 * tsq - d3 * tcube - d4 * tfour
+            tempe += bstar * c5 * (sin(xmp) - sinmo)
+            templ += t3cof * tcube + tfour * (t4cof + tSince * t5cof)
        }
+        val a = aodp * tempa * tempa
+        val eo = data.eccn
+        val e = eo - tempe
+        val xl = xmp + omega + xnode + xnodp * templ
+        val beta = sqrt(1.0 - e * e)
+        val xn = XKE / a.pow(1.5)
+
+        // Long period periodics
+        val axn = e * cos(omega)
+        temp[0] = invert(a * sqr(beta))
+        val xll = temp[0] * xlcof * axn
+        val aynl = temp[0] * aycof
+        val xlt = xl + xll
+        val ayn = e * sin(omega) + aynl
+
+        // Solve Kepler's equation
+        val capu = mod2PI(xlt - xnode)
+        temp[2] = capu
+        converge(temp, axn, ayn, capu)
+        calculatePosAndVel(temp, xnode, a, xn, axn, ayn)
+        calculatePhase(xlt, xnode, omgadf)
    }

    private fun calculatePosAndVel(
@@ -42,6 +42,27 @@ abstract class OrbitalObject(val data: OrbitalData) {
    var qoms24 = 0.0
    var s4 = 0.0

+    // Pre-allocated reusable vectors to avoid GC pressure in hot loops
+    private val obsPos = Vector4()
+    private val obsVel = Vector4()
+    private val rangeVector = Vector4()
+    private val rgvelVector = Vector4()
+    private val squintVector = Vector4()
+
+    // Cached observer position data to avoid recalculation when observer hasn't moved
+    private var cachedGsLat = Double.NaN
+    private var cachedGsLon = Double.NaN
+    private var cachedGsAlt = Double.NaN
+    private var cachedSinLat = 0.0
+    private var cachedCosLat = 0.0
+    private var cachedObsC = 0.0
+    private var cachedObsSq = 0.0
+    private var cachedObsAchFactor = 0.0
+    private var cachedObsZFactor = 0.0
+
+    // Cache for julian epoch to avoid recomputing every call
+    private val julEpoch: Double = juliandDateOfEpoch(data.epoch)
+
    fun willBeSeen(pos: GeoPos): Boolean {
        return if (data.meanmo < 1e-8) false
        else {
@@ -57,15 +78,13 @@ abstract class OrbitalObject(val data: OrbitalData) {
        orbitalPos = OrbitalPos()
        // Date/time at which the position and velocity were calculated
        julUTC = calcCurrentDaynum(time) + 2444238.5
-        // Convert satellite's epoch time to Julian and calculate time since epoch in minutes
-        val julEpoch = juliandDateOfEpoch(data.epoch)
+        // Calculate time since epoch in minutes
        val tsince = (julUTC - julEpoch) * MIN_PER_DAY
        calculateSDP4orSGP4(tsince)
        // Scale position and velocity vectors to km and km/sec
        convertSatState(position, velocity)
        // Calculate velocity of satellite
        magnitude(velocity)
-        val squintVector = Vector4()
        // Angles in rads, dist in km, vel in km/S. Calculate sat Az, El, Range and Range-rate.
        calculateObs(julUTC, position, velocity, pos, squintVector)
        calculateLatLonAlt(julUTC)
@@ -75,6 +94,38 @@ abstract class OrbitalObject(val data: OrbitalData) {
        return orbitalPos
    }

+    /**
+     * Lightweight elevation-only check for pass finding. Avoids full lat/lon/alt,
+     * eclipse, and squint calculations that aren't needed when just searching for
+     * horizon crossings.
+     */
+    fun getElevation(pos: GeoPos, time: Long): Double {
+        julUTC = calcCurrentDaynum(time) + 2444238.5
+        val tsince = (julUTC - julEpoch) * MIN_PER_DAY
+        calculateSDP4orSGP4(tsince)
+        convertSatState(position, velocity)
+        magnitude(velocity)
+        calculateObs(julUTC, position, velocity, pos, squintVector)
+        return orbitalPos.elevation
+    }
+
+    /**
+     * Full position calculation that also populates azimuth, altitude, etc.
+     * Used when we need all fields (AOS/LOS refinement, track computation).
+     */
+    fun getFullPosition(pos: GeoPos, time: Long): OrbitalPos {
+        orbitalPos = OrbitalPos()
+        julUTC = calcCurrentDaynum(time) + 2444238.5
+        val tsince = (julUTC - julEpoch) * MIN_PER_DAY
+        calculateSDP4orSGP4(tsince)
+        convertSatState(position, velocity)
+        magnitude(velocity)
+        calculateObs(julUTC, position, velocity, pos, squintVector)
+        calculateLatLonAlt(julUTC)
+        orbitalPos.time = time
+        return orbitalPos
+    }
+
    private fun calcCurrentDaynum(now: Long): Double {
        val then = 315446400000 // time in millis on 31Dec79 00:00:00 UTC (daynum 0)
        return (now - then) / 1000.0 / 60.0 / 60.0 / 24.0
@@ -117,38 +168,34 @@ abstract class OrbitalObject(val data: OrbitalData) {
        gsPos: GeoPos,
        squintVector: Vector4
    ) {
-        val obsPos = Vector4()
-        val obsVel = Vector4()
-        val range = Vector4()
-        val rgvel = Vector4()
        calculateUserPosVel(julianUTC, gsPos, obsPos, obsVel)
-        range.setXYZ(
+        rangeVector.setXYZ(
            positionVector.x - obsPos.x,
            positionVector.y - obsPos.y,
            positionVector.z - obsPos.z
        )
        // Save these values globally for calculating squint angles later
-        squintVector.setXYZ(range.x, range.y, range.z)
-        rgvel.setXYZ(
+        squintVector.setXYZ(rangeVector.x, rangeVector.y, rangeVector.z)
+        rgvelVector.setXYZ(
            velocityVector.x - obsVel.x,
            velocityVector.y - obsVel.y,
            velocityVector.z - obsVel.z
        )
-        magnitude(range)
-        val sinLat = sin(DEG2RAD * gsPos.latitude)
-        val cosLat = cos(DEG2RAD * gsPos.latitude)
+        magnitude(rangeVector)
+        val sinLat = cachedSinLat
+        val cosLat = cachedCosLat
        val sinTheta = sin(gsPosTheta)
        val cosTheta = cos(gsPosTheta)
-        val topS = sinLat * cosTheta * range.x + sinLat * sinTheta * range.y - cosLat * range.z
-        val topE = -sinTheta * range.x + cosTheta * range.y
-        val topZ = cosLat * cosTheta * range.x + cosLat * sinTheta * range.y + sinLat * range.z
+        val topS = sinLat * cosTheta * rangeVector.x + sinLat * sinTheta * rangeVector.y - cosLat * rangeVector.z
+        val topE = -sinTheta * rangeVector.x + cosTheta * rangeVector.y
+        val topZ = cosLat * cosTheta * rangeVector.x + cosLat * sinTheta * rangeVector.y + sinLat * rangeVector.z
        var azim = atan(-topE / topS)
        if (topS > 0.0) azim += PI
        if (azim < 0.0) azim += TWO_PI
        orbitalPos.azimuth = azim
-        orbitalPos.elevation = asin(topZ / range.w)
-        orbitalPos.distance = range.w
-        orbitalPos.distanceRate = dot(range, rgvel) / range.w
+        orbitalPos.elevation = asin(topZ / rangeVector.w)
+        orbitalPos.distance = rangeVector.w
+        orbitalPos.distanceRate = dot(rangeVector, rgvelVector) / rangeVector.w
        var elevation = orbitalPos.elevation / TWO_PI * 360.0
        if (elevation > 90) elevation = 180 - elevation
        orbitalPos.aboveHorizon = elevation - 0 > EPSILON
@@ -163,12 +210,24 @@ abstract class OrbitalObject(val data: OrbitalData) {
    ) {
        val mFactor = 7.292115E-5
        gsPosTheta = mod2PI(thetaGJD(time) + DEG2RAD * gsPos.longitude)
-        val c = invert(sqrt(1.0 + FLAT_FACT * (FLAT_FACT - 2) * sqr(sin(DEG2RAD * gsPos.latitude))))
-        val sq = sqr(1.0 - FLAT_FACT) * c
-        val achcp = (EARTH_RADIUS * c + gsPos.altitude / 1000.0) * cos(DEG2RAD * gsPos.latitude)
+
+        // Cache trig and position factors when observer position changes
+        if (gsPos.latitude != cachedGsLat || gsPos.longitude != cachedGsLon || gsPos.altitude != cachedGsAlt) {
+            cachedGsLat = gsPos.latitude
+            cachedGsLon = gsPos.longitude
+            cachedGsAlt = gsPos.altitude
+            cachedSinLat = sin(DEG2RAD * gsPos.latitude)
+            cachedCosLat = cos(DEG2RAD * gsPos.latitude)
+            cachedObsC = invert(sqrt(1.0 + FLAT_FACT * (FLAT_FACT - 2) * sqr(cachedSinLat)))
+            cachedObsSq = sqr(1.0 - FLAT_FACT) * cachedObsC
+            cachedObsAchFactor = (EARTH_RADIUS * cachedObsC + gsPos.altitude / 1000.0) * cachedCosLat
+            cachedObsZFactor = (EARTH_RADIUS * cachedObsSq + gsPos.altitude / 1000.0) * cachedSinLat
+        }
+
        obsPos.setXYZ(
-            achcp * cos(gsPosTheta), achcp * sin(gsPosTheta),
-            (EARTH_RADIUS * sq + gsPos.altitude / 1000.0) * sin(DEG2RAD * gsPos.latitude)
+            cachedObsAchFactor * cos(gsPosTheta),
+            cachedObsAchFactor * sin(gsPosTheta),
+            cachedObsZFactor
        )
        obsVel.setXYZ(-mFactor * obsPos.y, mFactor * obsPos.x, 0.0)
        magnitude(obsPos)
@@ -368,7 +427,7 @@ abstract class OrbitalObject(val data: OrbitalData) {
     * The function Delta_ET has been added to allow calculations on the
     * position of the sun. It provides the difference between UT (approximately
     * the same as UTC) and ET (now referred to as TDT) This function is based
-     * on a least squares fit of data from 1950 to 1991 and will need to be
+     * on the least squares fit of data from 1950 to 1991 and will need to be
     * updated periodically.
     *
     * Values determined using data from 1950-1991 in the 1990 Astronomical
@@ -21,7 +21,6 @@ import kotlin.math.acos
 import kotlin.math.asin
 import kotlin.math.atan2
 import kotlin.math.cos
-import kotlin.math.pow
 import kotlin.math.sin
 import kotlin.math.sqrt

@@ -50,23 +49,32 @@ data class OrbitalPos(
    }

    fun getOrbitalVelocity(): Double {
-        val earthG = 6.674 * 10.0.pow(-11)
-        val earthM = 5.98 * 10.0.pow(24)
-        val radius = 6.37 * 10.0.pow(6) + altitude * 10.0.pow(3)
-        return sqrt(earthG * earthM / radius) / 1000
+        val radius = EARTH_RADIUS_M + altitude * 1000.0
+        return sqrt(GM_EARTH / radius) / 1000.0
    }

    fun getRangeCircle(): List<GeoPos> {
-        val rangeCirclePoints = mutableListOf<GeoPos>()
-        val beta = acos(EARTH_RADIUS / (EARTH_RADIUS + altitude)) // * EARTH_RADIUS = radiusKm
+        val pointCount = 721
+        val rangeCirclePoints = ArrayList<GeoPos>(pointCount)
+        val beta = acos(EARTH_RADIUS / (EARTH_RADIUS + altitude))
+        val sinLat = sin(latitude)
+        val cosLat = cos(latitude)
+        val cosBeta = cos(beta)
+        val sinBeta = sin(beta)
        for (azimuth in 0..720) {
            val rads = azimuth * DEG2RAD
-            val lat = asin(sin(latitude) * cos(beta) + (cos(latitude) * sin(beta) * cos(rads)))
-            val lon = (longitude + atan2(
-                sin(rads) * sin(beta) * cos(latitude), cos(beta) - sin(latitude) * sin(lat)
-            ))
+            val sinRads = sin(rads)
+            val cosRads = cos(rads)
+            val lat = asin(sinLat * cosBeta + cosLat * sinBeta * cosRads)
+            val lon = longitude + atan2(sinRads * sinBeta * cosLat, cosBeta - sinLat * sin(lat))
            rangeCirclePoints.add(GeoPos(lat * RAD2DEG, lon * RAD2DEG))
        }
        return rangeCirclePoints
    }
+
+    companion object {
+        // Pre-computed constants for orbital velocity calculation
+        private const val GM_EARTH = 3.986004418E14 // m^3/s^2
+        private const val EARTH_RADIUS_M = 6.37E6 // meters
+    }
 }