Init

libs/ballistics-lib/src/main/kotlin/ballistics/Ballistics.kt

@@ -0,0 +1,368 @@
+package ballistics
+
+import ballistics.flightLine.FlightLineCalculator
+import ballistics.mpl.MPLCalculator
+import ballistics.orbitalPoints.AbstractOrbPointsCalculator
+import ballistics.orbitalPoints.OrbitalPointsIntegrator
+import ballistics.orbitalPoints.OrbitalPointsTLE
+import ballistics.orbitalPoints.timeStepper.AbstractStepper
+import ballistics.orbitalPoints.timeStepper.TLEStepper
+import ballistics.types.BallisticsError
+import ballistics.types.EarthType
+import ballistics.types.FleghtLineSector
+import ballistics.types.FlightLine
+import ballistics.types.InitialConditions
+import ballistics.types.IntegrationType
+import ballistics.types.KeplerParams
+import ballistics.types.ModDVType
+import ballistics.types.OPKatObj
+import ballistics.types.OrbitalPoint
+import ballistics.types.PPI
+import ballistics.types.PointViewParams
+import ballistics.types.RevolutionParameter
+import ballistics.types.TLE
+import ballistics.types.TLEParams
+import ballistics.types.TangageType
+import ballistics.types.WorkCSType
+import ballistics.types.ZRV
+import ballistics.utils.astro.AstronomerJ2000
+import ballistics.utils.toDateTime
+import ballistics.zrv.ZRVStepperCalculator
+import ballistics.zrv.ZRVStepperCalculatorAsync
+import org.nstart.dep265.tletools.zeptomoby.core.Globals.sqr
+import reactor.core.publisher.Flux
+import kotlin.math.PI
+import kotlin.math.asin
+import kotlin.math.atan
+import kotlin.math.atan2
+import kotlin.math.sin
+import kotlin.math.sqrt
+
+class Ballistics {
+    var kaId: Int = -1
+
+    var kaUmn: Int = -1
+
+    val ic : MutableList<InitialConditions> = mutableListOf()
+
+    /**
+     * Модель движения
+     */
+    var modDVType: ModDVType = ModDVType.FOTO
+
+    /**
+     * Метод интегрирования
+     */
+    var integrationType: IntegrationType = IntegrationType.ADAMS7
+
+    /**
+     * Тип модели Земли
+     */
+    var earthType: EarthType = EarthType.PZ90d02
+
+    /**
+     * тип базовой системы координат (для расчета трасс полета, полособзора и матрицы планирования)
+     */
+    var workCoordinateSystem: WorkCSType = WorkCSType.WCSOrbit
+
+    /**
+     * Минимальный угол крена (для расчета полосы обзора)
+     */
+    var rollMin: Double = 0.0
+
+    /**
+     * Максимальный угол крена (для расчета полосы обзора и покрытия земной поверхности)
+     */
+    var rollMax: Double = PI / 6.0
+
+
+    var rollByModule : Boolean = true
+
+    /**
+     * Граничный угол Солнца для КА (для определения попадания момента наблюдения объекта в матрицу планирования)
+     */
+    var sunAngleMin: Double = 0.0
+
+    /**
+     * Признак использования ограничений на попадание в матрицу планирования из параметров объекта
+     */
+    var useObjConstraints: Boolean = false
+
+    /**
+     * Объект класса расчета точек орбиты
+     */
+    private var orbitCalculator: AbstractOrbPointsCalculator? = null
+
+    /**
+     * Объект класса расчета параметров ЗРВ
+     */
+    private var zrvCalculator: ZRVStepperCalculator? = null
+
+    /**
+     * Объект класса расчета параметров трассы полета и полосы обзора
+     */
+    private var flightLineCalculator: FlightLineCalculator? = null
+
+    /**
+     * Объект класса расчета матрицы планирования
+     */
+    private var mplCalculator: MPLCalculator? = null
+
+    /**
+     * Массив точек орбиты
+     */
+    val points: Iterable<OrbitalPoint>
+        get() = if (orbitCalculator != null) orbitCalculator!!.points else listOf<OrbitalPoint>()
+
+    val revolutions: Iterable<RevolutionParameter>
+        get() = if (orbitCalculator != null) orbitCalculator!!.revolutions else listOf<RevolutionParameter>()
+
+    /**
+     * Массив результатов расчета параметров ЗРВ
+     */
+    val zrv: Iterable<ZRV>
+        get() = if (zrvCalculator != null) zrvCalculator!!.zrv else listOf<ZRV>()
+
+    val flightLine: Iterable<FlightLine>
+        get() = if (flightLineCalculator != null) flightLineCalculator!!.flightLine else listOf<FlightLine>()
+
+    val mpl: Iterable<PointViewParams>
+        get() = if (mplCalculator != null) mplCalculator!!.mpl else listOf<PointViewParams>()
+
+    val coverings : Map<Int, MutableList<FleghtLineSector>>
+        get() = flightLineCalculator?.earthCoveringCalculator?.coverings?:mapOf()
+
+
+    /**
+     * Расчет точек орбиты методом интегрирования
+     * @param nu начальные условия движения центра масс КА
+     * @param tn время начала расчета
+     * @param tk время конца расчета
+     * @return код завершения расчета
+     */
+    fun calculateOrbPoints(
+        nu: InitialConditions,
+        tn: Double,
+        tk: Double,
+    ): BallisticsError {
+        orbitCalculator = OrbitalPointsIntegrator(nu, modDVType, integrationType, earthType)
+        return orbitCalculator!!.calculate(tn, tk)
+    }
+
+    fun setOrbitalPoints(
+        points: Iterable<OrbitalPoint>,
+        earthType: EarthType = EarthType.PZ90d02,
+    ) {
+        orbitCalculator = OrbitalPointsIntegrator(arrayOf(), ModDVType.FOTO, IntegrationType.ADAMS7, earthType)
+        orbitCalculator?.points?.addAll(points)
+    }
+
+
+    fun setEarthCoverage(covs : Map<Int, MutableList<FleghtLineSector>>){
+        if (orbitCalculator != null){
+            flightLineCalculator = FlightLineCalculator(orbitCalculator!!, rollMin, rollMax, workCoordinateSystem)
+            flightLineCalculator?.earthCoveringCalculator?.coverings?.putAll(covs)
+        }
+    }
+
+    /**
+     * Расчет точек орбиты методом интегрирования
+     * @param nu массив начальных условий движения центра масс КА
+     * @param tn время начала расчета
+     * @param tk время конца расчета
+     * @return код завершения расчета
+     */
+    fun calculateOrbPoints(
+        nu: Array<InitialConditions>,
+        tn: Double,
+        tk: Double,
+    ): BallisticsError {
+        nu.sortBy { it.point.t }
+
+        orbitCalculator = OrbitalPointsIntegrator(nu, modDVType, integrationType, earthType)
+
+        return orbitCalculator!!.calculate(tn, tk)
+    }
+
+    /**
+     * Расчет параметров ЗРВ
+     * @param ppi список ППИ
+     * @param tn время начала расчета
+     * @param tk время конца расчета
+     * @return код завершения расчета
+     */
+    fun calculateZRV(
+        ppi: List<PPI>,
+        tn: Double,
+        tk: Double,
+    ): BallisticsError {
+        return orbitCalculator?.let {
+            zrvCalculator = ZRVStepperCalculator(it, ppi)
+            zrvCalculator?.calculate(tn, tk)
+        } ?: BallisticsError.EMPTY_ORBITAL_POINTS
+    }
+
+    fun calculateZRVAsync(
+        ppi: List<PPI>,
+        tn: Double,
+        tk: Double,
+    ): Flux<ZRV> =
+        if (orbitCalculator != null) {
+            ZRVStepperCalculatorAsync(orbitCalculator!!, ppi).calculate(tn, tk)
+        } else {
+            Flux.empty()
+        }
+
+    fun calculateFlightLine(
+        tn: Double,
+        tk: Double,
+    ): BallisticsError {
+        return orbitCalculator?.let {
+            flightLineCalculator = FlightLineCalculator(it, rollMin, rollMax, workCoordinateSystem)
+            flightLineCalculator?.calculate(tn, tk)
+        } ?: BallisticsError.EMPTY_ORBITAL_POINTS
+    }
+
+    fun calculateMPL(
+        tn: Double,
+        tk: Double,
+        objs: Iterable<OPKatObj>,
+        tangage: Double = 0.0,
+        tangageType: TangageType = TangageType.TTProactive,
+    ): BallisticsError {
+        return flightLineCalculator?.let {
+            mplCalculator =
+                MPLCalculator(it, tangage, tangageType).apply {
+                    this.sunAngleMin = this@Ballistics.sunAngleMin
+                    this.useObjConstraints = this@Ballistics.useObjConstraints
+                }
+            mplCalculator?.krenMax = rollMax
+            mplCalculator?.krenMin = rollMin
+            mplCalculator?.krenByModule = rollByModule
+            mplCalculator?.calculate(tn, tk, objs)
+        } ?: BallisticsError.EMPTY_FLIGHTLINE_POINTS
+    }
+
+    /**
+     * Объект расчета выхода на заданное время
+     */
+    fun getStepper(): AbstractStepper? {
+        return orbitCalculator?.getStepper()
+    }
+
+    fun parseTLE(
+        tle1: String,
+        tle2: String,
+    ): OrbitalPoint {
+        try {
+            val stepper = TLEStepper(tle1, tle2, EarthType.PZ90d02)
+            return stepper.calculate(stepper.baseEpoch)
+        }catch (ex : Exception){
+            throw ex
+        }
+    }
+
+    fun parseTLE(tle: TLE): OrbitalPoint {
+        val stepper = TLEStepper(tle.tle1, tle.tle2, EarthType.PZ90d02)
+        return stepper.calculate(stepper.baseEpoch)
+    }
+
+    fun getTLEParams(tle : TLE, capt : String = "") : TLEParams{
+        try {
+            val stepper = TLEStepper(tle.tle1, tle.tle2, EarthType.PZ90d02)
+            val orb = stepper.calculate(stepper.baseEpoch)
+            return TLEParams(
+                capt, stepper.satellite.orbit.satId.toLong(),
+                stepper.satellite.orbit.orbitNum.toLong(),
+                toDateTime(orb.t),
+                stepper.satellite.orbit.inclination,
+                stepper.satellite.orbit.perigee,
+                stepper.satellite.orbit.apogee,
+                stepper.satellite.orbit.argPerigee,
+                stepper.satellite.orbit.eccentricity,
+                stepper.satellite.orbit.major,
+                stepper.satellite.orbit.minor,
+                stepper.satellite.orbit.meanAnomaly,
+                stepper.satellite.orbit.period,
+                stepper.satellite.orbit.semiMajor,
+                stepper.satellite.orbit.semiMinor,
+                stepper.satellite.orbit.raan,
+                stepper.satellite.orbit.meanMotion,
+                stepper.satellite.orbit.meanMotionTle,
+            )
+        }catch (ex : Exception){
+            throw ex
+        }
+    }
+
+    fun calculateOrbPoints(
+        tle: TLE,
+        tn: Double,
+        tk: Double,
+    ): BallisticsError {
+        orbitCalculator = OrbitalPointsTLE(tle.tle1, tle.tle2, earthType)
+        return orbitCalculator!!.calculate(tn, tk)
+    }
+
+    fun calculateOrbPoints(
+        tle: TLE,
+        duration: Double,
+    ): BallisticsError {
+        orbitCalculator = OrbitalPointsTLE(tle.tle1, tle.tle2, earthType)
+        val tnu = (orbitCalculator as OrbitalPointsTLE).stepper.baseEpoch
+        return orbitCalculator!!.calculate(tnu, tnu + duration)
+    }
+
+    fun calculateKeplerParams(point: OrbitalPoint): KeplerParams {
+        val res = KeplerParams()
+
+        val aL00 = 6.25648106E+7
+        val astro = AstronomerJ2000(earthType)
+        val ask = astro.grinvToASK(point)
+
+        val r = sqrt(sqr(ask.r.x) + sqr(ask.r.y) + sqr(ask.r.z))
+        val v = sqrt(sqr(ask.v.x) + sqr(ask.v.y) + sqr(ask.v.z))
+
+        val mu = astro.earth.fM
+        val k = (r * sqr(v)) / mu
+
+        val sinO = (ask.r.x * ask.v.x + ask.r.y * ask.v.y + ask.r.z * ask.v.z) / (r * v)
+        val cosO = sqrt(1 - sqr(sinO))
+
+        val c1 = ask.r.y * ask.v.z - ask.r.z * ask.v.y
+        val c2 = ask.r.z * ask.v.x - ask.r.x * ask.v.z
+        val c3 = ask.r.x * ask.v.y - ask.r.y * ask.v.x
+        val c = sqrt(sqr(c1) + sqr(c2) + sqr(c3))
+
+        res.ael = r / (2 - k)
+        val lambda = (1 / res.ael) * sqrt(mu / res.ael)
+
+        res.t = (2 * PI * res.ael * sqrt(res.ael)) / sqrt(astro.earth.middleRadius * aL00)
+        res.e = sqrt(1 - k * (2 - k) * (1 - sinO * sinO))
+        res.v = atan((k * sinO * cosO) / (k * sqr(cosO) - 1))
+        res.u = atan2((ask.r.z * c), (ask.r.y * c1 - ask.r.x * c2))
+        res.eA = atan((sqrt(1 - sqr(res.e)) * k * sinO * cosO) / (sqr(res.e) + (k * sqr(cosO) - 1)))
+        res.tau = ask.t + ((res.e * sin(res.eA) - res.eA) / lambda)
+        res.nakl = 0.5 * PI - asin(c3 / c)
+        res.omegab = atan2(c1, -c2)
+        res.omegam = res.u - res.v
+        res.o = asin(sinO)
+        res.rA = res.ael * (1 + res.e)
+        res.rP = res.ael * (1 - res.e)
+        res.dmv = point.t
+
+        if (res.omegam < 0) res.omegam += 2 * PI
+        if (res.tau < 0) res.tau += 1
+
+        return res
+    }
+
+
+    fun clear(){
+        orbitCalculator?.clear()
+        flightLineCalculator?.clear()
+        flightLineCalculator?.clear()
+        zrvCalculator?.clear()
+        mplCalculator?.clear()
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/flightLine/EarthCovering.kt

@@ -0,0 +1,105 @@
+package ballistics.flightLine
+
+import ballistics.types.FleghtLineSector
+import ballistics.types.FlightLine
+import kotlin.math.PI
+import kotlin.math.truncate
+
+internal class EarthCovering {
+    var coverings = mutableMapOf<Int, MutableList<FleghtLineSector>>()
+
+    fun getBInd(b: Double): Int {
+        var b2: Double = b * 180.0 / PI
+        b2 = 90.0 - b2
+        var k: Int = truncate(b2 / 2.0).toInt()
+        if (k > 89) k = 89
+        return k
+    }
+
+    fun getLInd(l2: Double): Int {
+        var l: Double = l2 * 180.0 / PI
+        if (l < 0) l += 360.0
+        var k: Int = truncate(l / 2.0).toInt()
+        if (k > 179) k = 179
+        return k
+    }
+
+    fun getInd(
+        b: Double,
+        l: Double,
+    ): Int {
+        return getBInd(b) * 180 + getLInd(l)
+    }
+
+    fun addCoverings(
+        frst: FlightLine,
+        sec: FlightLine,
+    ) {
+        val bmax = getBInd(minOf(frst.leftOuterSwath.lat, frst.rightOuterSwath.lat, sec.leftOuterSwath.lat, sec.rightOuterSwath.lat))
+        val bmin = getBInd(maxOf(frst.leftOuterSwath.lat, frst.rightOuterSwath.lat, sec.leftOuterSwath.lat, sec.rightOuterSwath.lat))
+        var lmin = getLInd(minOf(frst.leftOuterSwath.long, frst.rightOuterSwath.long, sec.leftOuterSwath.long, sec.rightOuterSwath.long))
+        var lmax = getLInd(maxOf(frst.leftOuterSwath.long, frst.rightOuterSwath.long, sec.leftOuterSwath.long, sec.rightOuterSwath.long))
+
+        if (lmax - lmin > 90) {
+            lmax = checkLmax(frst, sec)
+            lmin = checkLmin(frst, sec)
+            val l = lmin
+            lmin = lmax
+            lmax = l + 179
+        }
+
+        for (bi in bmin..bmax) {
+            for (li in lmin..lmax) {
+                val ll = if (li >= 180) li - 180 else li
+                val ind = bi * 180 + ll
+                if (coverings.containsKey(ind)) {
+                    coverings[ind]!!.add(FleghtLineSector(frst.t, sec.t))
+                } else {
+                    coverings.put(ind, mutableListOf(FleghtLineSector(frst.t, sec.t)))
+                }
+            }
+        }
+    }
+
+
+    fun checkLmax(frst: FlightLine,  sec: FlightLine) : Int{
+        var buf = 2 * PI
+        if (frst.leftOuterSwath.long > PI / 2) {
+            buf = frst.leftOuterSwath.long
+        }
+        if (frst.rightOuterSwath.long > PI / 2 && frst.rightOuterSwath.long < buf) {
+            buf = frst.rightOuterSwath.long
+        }
+        if (sec.rightOuterSwath.long > PI / 2 && sec.rightOuterSwath.long < buf) {
+            buf = sec.rightOuterSwath.long
+        }
+        if (sec.leftOuterSwath.long > PI / 2 && sec.leftOuterSwath.long < buf) {
+            buf = sec.leftOuterSwath.long
+        }
+        return getLInd(buf)
+    }
+
+    fun checkLmin(frst: FlightLine,  sec: FlightLine) : Int{
+        var buf = 0.0
+        if (frst.leftOuterSwath.long < PI / 2) {
+            buf = frst.leftOuterSwath.long
+        }
+        if (frst.rightOuterSwath.long < PI / 2 && frst.rightOuterSwath.long > buf) {
+            buf = frst.rightOuterSwath.long
+        }
+        if (sec.rightOuterSwath.long < PI / 2 && sec.rightOuterSwath.long > buf) {
+            buf = sec.rightOuterSwath.long
+        }
+        if (sec.leftOuterSwath.long < PI / 2 && sec.leftOuterSwath.long > buf) {
+            buf = sec.leftOuterSwath.long
+        }
+        return getLInd(buf)
+    }
+
+    fun getCoverings(
+        b: Double,
+        l: Double,
+    ): Iterable<FleghtLineSector> {
+        return coverings.getOrDefault(getInd(b, l), listOf<FleghtLineSector>())
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/flightLine/FlightLineCalculator.kt

@@ -0,0 +1,82 @@
+package ballistics.flightLine
+
+import ballistics.orbitalPoints.AbstractOrbPointsCalculator
+import ballistics.types.BallisticsError
+import ballistics.types.FlightLine
+import ballistics.types.OrbitalPoint
+import ballistics.types.Orientation
+import ballistics.types.WorkCSType
+import ballistics.utils.earth.getEarth
+
+
+internal class FlightLineCalculator(
+    var opc: AbstractOrbPointsCalculator,
+    var rollMin: Double,
+    var rollMax: Double,
+    var wcs: WorkCSType = WorkCSType.WCSOrbit,
+) {
+    var step = 60.0
+    val eart = getEarth(opc.earthType)
+    var flightLine = mutableListOf<FlightLine>()
+    var pc = PointOnEarthCalculator(opc.earthType, wcs)
+    var earthCoveringCalculator: EarthCovering = EarthCovering()
+    var needToCalcEarthCovering: Boolean = true
+
+
+    fun calculate(
+        tn: Double,
+        tk: Double,
+    ): BallisticsError {
+        flightLine.clear()
+        earthCoveringCalculator.coverings.clear()
+
+        var t = tn
+
+        val stepper = opc.getStepper()
+        var point: OrbitalPoint?
+
+        point = stepper.calculate(t)
+        if (point == null)
+            return BallisticsError.STEPPER_ERROR
+
+        flightLine.add(calcFl(point))
+        t += step
+
+        while (t <= tk) {
+            point = stepper.calculate(t)
+            if (point == null) {
+                println("Ошибка выхода на заданное время в середине расчета")
+                return BallisticsError.STEPPER_ERROR
+            }
+            val fl = calcFl(point)
+
+            if (needToCalcEarthCovering) {
+                earthCoveringCalculator.addCoverings(flightLine.last(), fl)
+            }
+
+            flightLine.add(fl)
+            t += step
+        }
+
+        return BallisticsError.OK
+    }
+
+    private fun calcFl(point: OrbitalPoint): FlightLine {
+        return FlightLine(
+            point.t,
+            point.vit,
+            if (point.v.z >= 0) 0 else 1,
+            pc.pointOnEarth(point, Orientation(0.0, -rollMax, 0.0))!!,
+            pc.pointOnEarth(point, Orientation(0.0, -rollMin, 0.0))!!,
+            pc.pointOnEarth(point, Orientation(0.0, 0.0, 0.0))!!,
+            pc.pointOnEarth(point, Orientation(0.0, rollMin, 0.0))!!,
+            pc.pointOnEarth(point, Orientation(0.0, rollMax, 0.0))!!,
+        )
+    }
+
+    fun clear(){
+        flightLine.clear()
+        earthCoveringCalculator.coverings.clear()
+        opc.clear()
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/flightLine/PointOnEarthCalculator.kt

@@ -0,0 +1,191 @@
+package ballistics.flightLine
+
+import ballistics.types.EarthType
+import ballistics.types.OrbitalPoint
+import ballistics.types.Orientation
+import ballistics.types.THBLPoint
+import ballistics.types.WorkCSType
+import ballistics.utils.astro.AstronomerJ2000
+import ballistics.utils.math.Matrix3D
+import ballistics.utils.math.Vector3D
+import kotlin.math.asin
+import kotlin.math.cos
+import kotlin.math.pow
+import kotlin.math.sin
+import kotlin.math.sqrt
+
+class PointOnEarthCalculator(val earthType: EarthType, private val wcs: WorkCSType) {
+    private var astro = AstronomerJ2000(earthType)
+
+    fun pointOnEarth(
+        point: OrbitalPoint,
+        orientation: Orientation,
+    ): THBLPoint? {
+        val krenMax = asin((astro.earth.polarRadius - 50000.0) / point.r.module())
+        if (orientation.kren > krenMax) {
+            orientation.kren = krenMax
+        }
+        if (orientation.kren < -krenMax) {
+            orientation.kren = -krenMax
+        }
+
+        val rabs: Vector3D
+        val vabs: Vector3D
+
+        if (wcs == WorkCSType.WCSPath) {
+            rabs = point.r
+            vabs = point.v
+        } else {
+            val ask = astro.grinvToASK(point)
+
+            rabs = ask.r
+            vabs = ask.v
+        }
+        val dd: Vector3D
+        val c = Matrix3D()
+        c.second = rabs.basis()
+        c.third = (vabs.rem(rabs)).basis()
+        c.first = c.second.rem(c.third)
+        var ct = c.transpose()
+        val aa =
+            Matrix3D(
+                Vector3D(
+                    cos(orientation.tang) * cos(orientation.risk) - sin(orientation.tang) * sin(orientation.kren) * sin(orientation.risk),
+                    -sin(orientation.tang) * cos(orientation.kren),
+                    cos(orientation.tang) * sin(orientation.risk) + sin(orientation.tang) * sin(orientation.kren) * cos(orientation.risk),
+                ),
+                Vector3D(
+                    sin(orientation.tang) * cos(orientation.risk) + cos(orientation.tang) * sin(orientation.kren) * sin(orientation.risk),
+                    cos(orientation.tang) * cos(orientation.kren),
+                    sin(orientation.tang) * sin(orientation.risk) - cos(orientation.tang) * sin(orientation.kren) * cos(orientation.risk),
+                ),
+                Vector3D(
+                    -cos(orientation.kren) * sin(orientation.risk),
+                    sin(orientation.kren),
+                    cos(orientation.kren) * cos(orientation.risk),
+                ),
+            )
+        if (wcs == WorkCSType.WCSOrbit) {
+            val g = Matrix3D()
+            g.makeOzMatrix(astro.si2000(point.t))
+            ct = g.transpose() * ct
+        }
+
+        val k = ct * aa
+
+        val d = Vector3D(0.0, 1.0, 0.0).basis()
+
+        dd = k * d
+
+        return earthIntersection(point.t, point.r, dd)
+    }
+
+    /**
+     * расчет точки пересечени отрезка, определенного вектором dd из точка rotn и земной поверхности
+     */
+    private fun earthIntersection(
+        t: Double,
+        rotn: Vector3D,
+        dd: Vector3D,
+    ): THBLPoint? {
+        var dz = dd.z
+        if (dz < 0) dz *= -1
+        val x1: Double
+        val y1: Double
+        val z1: Double
+        val x2: Double
+        val y2: Double
+        val z2: Double
+        val r: Vector3D
+        val h: Double
+
+        val ekvRadiusOverPolarRadiusSquared = (astro.earth.ekvRadius / astro.earth.polarRadius).pow(2)
+        if (dz < 0.1) {
+            val k1 = dd.y / dd.x
+            val k2 = rotn.y - k1 * rotn.x
+            val k3 = dd.z / dd.x
+            val k4 = rotn.z - k3 * rotn.x
+            val k5 = k1 * k1 + ekvRadiusOverPolarRadiusSquared * k3 * k3 + 1
+            val k6 = 2 * (k1 * k2 + ekvRadiusOverPolarRadiusSquared * k3 * k4)
+            val k7 = k2 * k2 + ekvRadiusOverPolarRadiusSquared * k4 * k4 - astro.earth.ekvRadius * astro.earth.ekvRadius
+
+            if ((k6 * k6 - k5 * k7 * 4) < 0) {
+                return null
+            }
+
+
+            x1 = (-k6 + sqrt(k6 * k6 - 4 * k5 * k7)) / (2 * k5)
+            y1 = k1 * x1 + k2
+            z1 = k3 * x1 + k4
+
+
+            x2 = (-k6 - sqrt(k6 * k6 - 4 * k5 * k7)) / (2 * k5)
+            y2 = k1 * x2 + k2
+            z2 = k3 * x2 + k4
+        } else {
+            val k1 = dd.x / dd.z
+            val k2 = rotn.x - k1 * rotn.z
+            val k3 = dd.y / dd.z
+            val k4 = rotn.y - k3 * rotn.z
+            val k5 = k1 * k1 + k3 * k3 + ekvRadiusOverPolarRadiusSquared
+            val k6 = 2 * (k1 * k2 + k3 * k4)
+            val k7 = k2 * k2 + k4 * k4 - astro.earth.ekvRadius * astro.earth.ekvRadius
+
+            if ((k6 * k6 - 4 * k5 * k7) < 0) {
+                return null
+            }
+
+
+            z1 = (-k6 + sqrt(k6 * k6 - 4 * k5 * k7)) / (2 * k5)
+            y1 = k3 * z1 + k4
+            x1 = k1 * z1 + k2
+
+
+            z2 = (-k6 - sqrt(k6 * k6 - 4 * k5 * k7)) / (2 * k5)
+            y2 = k3 * z2 + k4
+            x2 = k1 * z2 + k2
+        }
+        val d1 = Vector3D(rotn.x - x1, rotn.y - y1, rotn.z - z1).module()
+        val d2 = Vector3D(rotn.x - x2, rotn.y - y2, rotn.z - z2).module()
+
+
+        if (d1 <= d2) {
+            r = Vector3D(x1, y1, z1)
+            h = d1
+        } else {
+            r = Vector3D(x2, y2, z2)
+            h = d2
+        }
+
+        val blh = astro.earth.xyz2blh(r)
+
+        return THBLPoint(blh.lat, blh.long, h, astro.sunAngle(t, r))
+    }
+
+    /**
+     *  скорость компенсации (чтобы получить СДИ надо домножить на фокусное расстояние в мм)
+     */
+    fun calculateWD(
+        point: OrbitalPoint,
+        orientation: Orientation,
+        range: Double,
+    ): Double {
+        try {
+            val r = point.r.module()
+
+            val c1: Double = point.r.y * point.v.z - point.r.z * point.v.y
+            val c2: Double = point.r.z * point.v.x - point.r.x * point.v.z
+            val c3: Double = point.r.x * point.v.y - point.r.y * point.v.x
+            val cc = sqrt(c1 * c1 + c2 * c2 + c3 * c3)
+
+            // Трансверсальная скорость V*sin(Q)
+            val worb = cc / r
+
+            val wd = worb * (r - range * cos(orientation.kren)) / r / range
+
+            return wd
+        } catch (ex: Exception) {
+            return 0.0
+        }
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/mpl/MPLCalculator.kt

@@ -0,0 +1,143 @@
+package ballistics.mpl
+
+import ballistics.flightLine.FlightLineCalculator
+import ballistics.types.BallisticsError
+import ballistics.types.OPKatObj
+import ballistics.types.PointViewParams
+import ballistics.types.TangageType
+import ballistics.utils.math.Vector3D
+import ballistics.utils.math.equations.EquationCalculatorSpan
+import kotlin.math.PI
+import kotlin.math.abs
+
+internal class MPLCalculator(val flightLineCalculator: FlightLineCalculator, val tangage: Double, tangageType: TangageType) {
+    private val debugEasyCalc = false
+    private val stepper = flightLineCalculator.opc.getStepper()
+    private val orientCalculator = OrientOnPointCalculator(flightLineCalculator.opc.earthType, flightLineCalculator.wcs, tangageType)
+    private var currentObj = Vector3D()
+    var sunAngleMin: Double = 0.0
+    var useObjConstraints: Boolean = false
+
+    var krenMin = 0.0
+    var krenMax = 0.0
+    var krenByModule : Boolean = true
+
+    var mpl = mutableListOf<PointViewParams>()
+
+    fun calculate(
+        tn: Double,
+        tk: Double,
+        objs: Iterable<OPKatObj>,
+    ): BallisticsError {
+        mpl.clear()
+        val eqc = EquationCalculatorSpan()
+        eqc.delta = 0.001 * PI / 180.0
+        eqc.value = tangage
+
+        val samin = sunAngleMin
+//        krenMin = flightLineCalculator.rollMin
+//        krenMax = flightLineCalculator.rollMax
+
+
+        for (o in objs) {
+            if (useObjConstraints) {
+                krenMin = o.rollMin
+                krenMax = o.rollMax
+                sunAngleMin = o.sunAngleMin
+            }
+
+            currentObj = flightLineCalculator.eart.blh2xyz(o.lat, o.long, o.height)
+
+            val covs =
+                flightLineCalculator.earthCoveringCalculator.getCoverings(
+                    o.lat,
+                    o.long,
+                ).filter { (it.tStop >= tn && it.tStart <= tk) }
+            var tLast = 0.0
+            for (cov in covs) {
+                if (cov.tStart - tLast < 600.0) {
+                    continue
+                }
+
+                var isOk: Boolean
+                var traverz: Double?
+                if (!debugEasyCalc) {
+                    traverz = eqc.calculate(cov.tStart, cov.tStop, this::equation)
+                    isOk = traverz != null
+                } else {
+                    traverz = (cov.tStart + cov.tStop) / 2
+                    isOk = true
+                }
+
+                if (isOk &&  tryToAddViewParams(traverz!!, o)) {
+                        tLast = traverz
+                }
+            }
+        }
+        sunAngleMin = samin
+        mpl.sortBy { it.traverz }
+        return BallisticsError.OK
+    }
+
+    fun tryToAddViewParams(
+        t: Double,
+        o: OPKatObj,
+    ): Boolean {
+        val ka = stepper.calculate(t)
+        if (ka != null) {
+            val orient = orientCalculator.calculateOrientOnPoint(ka, currentObj)
+            val w = orientCalculator.pointInWCS(ka, currentObj)
+            val blh = flightLineCalculator.eart.xyz2blh(ka.r)
+
+            val vp =
+                PointViewParams(
+                    o.objON,
+                    o.objN,
+                    o.objUUID,
+                    o.pointNumb,
+                    ka.vit,
+                    t,
+                    blh.lat,
+                    blh.long,
+                    orient,
+                    w.module(),
+                    orientCalculator.astro.sunAngle(t, currentObj),
+                    orientCalculator.calculateVisirAngle(ka, o.lat, o.long, o.height),
+                    if (ka.v.z > 0) 0 else 1
+                    )
+
+            if ((vp.sunAngle >= sunAngleMin) && ( ( krenByModule &&
+                (abs(vp.orientation.kren) >= krenMin) &&
+                (abs(vp.orientation.kren) <= krenMax)
+                        ) ||
+                        ( !krenByModule &&
+                                (vp.orientation.kren >= krenMin) &&
+                                (vp.orientation.kren <= krenMax)
+                                )
+                        )
+            ) {
+                mpl.add(vp)
+            }
+            return true
+        }
+
+        return false
+    }
+
+    private fun equation(x: Double): Double {
+        var tang = -1.0
+        val ka = stepper.calculate(x)
+        if (ka != null) {
+            val orient = orientCalculator.calculateOrientOnPoint(ka, currentObj)
+            tang = orient.tang
+        }
+        return tang
+    }
+
+    fun clear(){
+        flightLineCalculator.clear()
+        stepper.clear()
+        mpl.clear()
+
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/mpl/OrientOnPointCalculator.kt

@@ -0,0 +1,113 @@
+package ballistics.mpl
+
+import ballistics.types.EarthType
+import ballistics.types.OrbitalPoint
+import ballistics.types.Orientation
+import ballistics.types.TangageType
+import ballistics.types.WorkCSType
+import ballistics.utils.astro.AstronomerJ2000
+import ballistics.utils.math.Matrix3D
+import ballistics.utils.math.Vector3D
+import kotlin.math.PI
+import kotlin.math.abs
+import kotlin.math.atan
+import kotlin.math.atan2
+import kotlin.math.cos
+import kotlin.math.sign
+import kotlin.math.sin
+import kotlin.math.sqrt
+
+internal class OrientOnPointCalculator(val earthType: EarthType, val wcs: WorkCSType, val tangType: TangageType) {
+    var astro = AstronomerJ2000(earthType)
+
+    fun pointInWCS(
+        ka: OrbitalPoint,
+        point: Vector3D,
+    ): Vector3D {
+        var w: Vector3D
+        if (wcs == WorkCSType.WCSOrbit) {
+
+            val kaA = astro.grinvToASK(ka)
+
+            val pa = astro.grinvToASK(point, ka.t)
+
+
+
+            val c = Matrix3D()
+            c.second = kaA.r.basis()
+            c.third = (kaA.v.rem(kaA.r)).basis()
+            c.first = c.second.rem(c.third)
+
+
+            val d: Vector3D = pa - kaA.r
+
+            w = c * d
+        } else {
+            val c = Matrix3D()
+            c.second = ka.r.basis()
+            c.third = (ka.v.rem(ka.r)).basis()
+            c.first = c.second.rem(c.third)
+
+
+            val d: Vector3D = point - ka.r
+
+            w = c * d
+        }
+
+        return w
+    }
+
+    fun calculateOrientOnPoint(
+        ka: OrbitalPoint,
+        point: Vector3D,
+    ): Orientation {
+        val orient = Orientation(0.0, 0.0, 0.0)
+        val w = pointInWCS(ka, point)
+
+        if (tangType == TangageType.TTProactive) {
+            orient.tang = PI - atan2(w.x, w.y)
+            orient.kren = atan2(w.z, w.y / cos(orient.tang)) - PI
+        } else {
+            orient.kren = atan2(w.z, w.y) - PI
+            orient.tang = PI - atan2(w.x, w.z / sin(orient.kren))
+        }
+
+        if (abs(orient.tang) > PI) {
+            orient.tang -= sign(orient.tang) * PI * 2
+        }
+        if (abs(orient.kren) > PI) {
+            orient.kren -= sign(orient.kren) * PI * 2
+        }
+
+        return orient
+    }
+
+    fun calculateVisirAngle(
+        ka: OrbitalPoint,
+        b: Double,
+        l: Double,
+        h: Double,
+    ): Double {
+        val rp = astro.earth.blh2xyz(b, l, h)
+        val r1 = Vector3D()
+        val r = Vector3D()
+
+        r1.x = ka.r.z - rp.z
+        r1.y = ka.r.x - rp.x
+        r1.z = ka.r.y - rp.y
+
+        val r4 = r1.y * cos(l) + r1.z * sin(l)
+
+        r.x = cos(b) * r1.x - sin(b) * r4
+        r.y = sin(b) * r1.x + cos(b) * r4
+        r.z = cos(l) * r1.z - sin(l) * r1.y
+
+        var angV: Double = sqrt(r.x * r.x + r.z * r.z)
+        if (angV != 0.0) {
+            angV = atan(r.y / sqrt(r.x * r.x + r.z * r.z))
+            angV = PI / 2 - angV
+        }
+
+        return angV
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/orbitalPoints/AbstractOrbPointsCalculator.kt

@@ -0,0 +1,40 @@
+package ballistics.orbitalPoints
+
+import ballistics.orbitalPoints.timeStepper.AbstractStepper
+import ballistics.types.BallisticsError
+import ballistics.types.EarthType
+import ballistics.types.OrbitalPoint
+import ballistics.types.RevolutionParameter
+import ballistics.utils.earth.getEarth
+
+abstract class AbstractOrbPointsCalculator(var earthType: EarthType) {
+    val points = mutableListOf<OrbitalPoint>()
+    val revolutions = mutableListOf<RevolutionParameter>()
+    val earth = getEarth(earthType)
+    var step = 60.0
+
+    abstract fun calculate(
+        tbegin: Double,
+        tend: Double,
+    ): BallisticsError
+
+    abstract fun getStepper(): AbstractStepper
+
+    protected abstract fun fastStep(t: Double): OrbitalPoint
+
+    protected fun addRevolution(t: Double) {
+        val p = fastStep(t)
+        val blh = earth.xyz2blh(p.r)
+        revolutions.add(RevolutionParameter(p, blh.long, blh.h))
+    }
+
+    protected fun equation(x: Double): Double {
+        val p = fastStep(x)
+        return p.r.z
+    }
+
+    fun clear(){
+        points.clear()
+        revolutions.clear()
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/orbitalPoints/OrbitalPointsIntegrator.kt

@@ -0,0 +1,232 @@
+package ballistics.orbitalPoints
+
+import ballistics.orbitalPoints.integrator.AbstractIntegrator
+import ballistics.orbitalPoints.integrator.IntegratorAdams
+import ballistics.orbitalPoints.integrator.IntegratorRK4
+import ballistics.orbitalPoints.timeStepper.AbstractStepper
+import ballistics.orbitalPoints.timeStepper.RungeStepper
+import ballistics.types.BallisticsError
+import ballistics.types.EarthType
+import ballistics.types.InitialConditions
+import ballistics.types.IntegrationType
+import ballistics.types.ModDVType
+import ballistics.types.OrbitalPoint
+import ballistics.utils.math.Vector3D
+import ballistics.utils.math.equations.AbstractEquationCalculator
+import ballistics.utils.math.equations.EquationCalculatorSpan
+import kotlin.math.abs
+
+internal class OrbitalPointsIntegrator(val moddv: ModDVType, integratorType: IntegrationType, earthType: EarthType) :
+    AbstractOrbPointsCalculator(earthType) {
+    private val integrator: AbstractIntegrator
+    private val fastStepper: IntegratorRK4
+    var mic = mutableListOf<InitialConditions>()
+    private var x = arrayOf<Double>()
+    private var vit = 0
+    val equationCalculator: AbstractEquationCalculator = EquationCalculatorSpan()
+    var mTn : Double = 0.0
+
+    init {
+        when (integratorType) {
+            IntegrationType.RUNG4 -> {
+                integrator = IntegratorRK4(moddv, earthType)
+                fastStepper = IntegratorRK4(moddv, earthType)
+            }
+            IntegrationType.ADAMS7 -> {
+                integrator = IntegratorAdams(moddv, earthType)
+                fastStepper = IntegratorRK4(moddv, earthType)
+            }
+        }
+        equationCalculator.delta = 0.001
+    }
+
+    constructor(
+        nu: InitialConditions,
+        modDVType: ModDVType,
+        integratorType: IntegrationType,
+        earthType: EarthType,
+    ) : this(modDVType, integratorType, earthType) {
+        mic.clear()
+        mic.add(nu)
+    }
+    constructor(
+        nu: Array<InitialConditions>,
+        modDVType: ModDVType,
+        integratorType: IntegrationType,
+        earthType: EarthType,
+    ) : this(modDVType, integratorType, earthType) {
+        mic.clear()
+        mic.addAll(nu)
+    }
+
+
+    /**
+     * Расчет точек орбиты
+     * @param nu начальные условия движения центра масс КА
+     * @param tn время начала расчета
+     * @param tk время конца расчета
+     */
+    fun calculate(
+        nu: InitialConditions,
+        tn: Double,
+        tk: Double,
+    ): BallisticsError {
+        mTn = tn
+        var t = nu.point.t
+        x = arrayOf(nu.point.v.x, nu.point.v.y, nu.point.v.z, nu.point.r.x, nu.point.r.y, nu.point.r.z, nu.point.t)
+
+        integrator.setSBall(nu.sBall)
+        fastStepper.setSBall(nu.sBall)
+        integrator.setF81(nu.f81)
+        fastStepper.setF81(nu.f81)
+        integrator.accelerate(x)
+        points.clear()
+        revolutions.clear()
+
+        val equationCalculator: AbstractEquationCalculator = EquationCalculatorSpan()
+        equationCalculator.delta = 0.001
+
+        vit = nu.point.vit
+
+
+        if (abs(nu.point.r.z) < 0.00001 && nu.point.t >= tn) {
+            addRevolution(nu.point.t)
+        }
+
+        while (t <= tk) {
+            integrator.nextStep(x)
+            nextPoint(t)
+            for (i in 0..6)
+                x[i] = integrator.y[i]
+            t += integrator.step
+        }
+
+        return BallisticsError.OK
+    }
+
+
+    fun nextPoint(t : Double){
+        if (x[5] < 0 && integrator.y[5] >= 0) {
+            ++vit
+            val tvuz = equationCalculator.calculate(t, t + integrator.step, this::equation)
+            tvuz?.let {
+                if (tvuz >= mTn) {
+                    addRevolution(it)
+                }
+            }
+        }
+        if (t >= mTn) {
+            val p = OrbitalPoint(t, vit, Vector3D(x[3], x[4], x[5]), Vector3D(x[0], x[1], x[2]))
+            points.add(p)
+        }
+    }
+
+
+    /**
+     * Расчет точек орбиты
+     * @param mIC массив начальных условий движения центра масс КА
+     * @param tn время начала расчета
+     * @param tk время конца расчета
+     */
+    fun calculate(
+        mIC: Array<InitialConditions>,
+        tn: Double,
+        tk: Double,
+    ): BallisticsError {
+        points.clear()
+        revolutions.clear()
+        mTn = tn
+
+
+        if (mIC.isEmpty())
+            return BallisticsError.EMPTY_NU
+
+
+        var t = mIC[0].point.t
+        x =
+            arrayOf(
+                mIC[0].point.v.x,
+                mIC[0].point.v.y,
+                mIC[0].point.v.z,
+                mIC[0].point.r.x,
+                mIC[0].point.r.y,
+                mIC[0].point.r.z,
+                mIC[0].point.t,
+            )
+        integrator.setSBall(mIC[0].sBall)
+        fastStepper.setSBall(mIC[0].sBall)
+        integrator.setF81(mIC[0].f81)
+        fastStepper.setF81(mIC[0].f81)
+        integrator.accelerate(x)
+        vit = mIC[0].point.vit
+
+        var tkk: Double
+        for (i in mIC.indices) {
+            tkk = if (i == mIC.size - 1) tk else mIC[i + 1].point.t
+            while (t < tkk) {
+                integrator.nextStep(x)
+
+                nextPoint(t)
+                for (j in 0..6)
+                    x[j] = integrator.y[j]
+                t += integrator.step
+            }
+
+            if (i < mIC.size - 1) {
+                if (points.isNotEmpty() && points.last().t >= mIC[i + 1].point.t) {
+                    points.removeLast()
+                }
+                arrayOf(
+                    mIC[i + 1].point.v.x,
+                    mIC[i + 1].point.v.y,
+                    mIC[i + 1].point.v.z,
+                    mIC[i + 1].point.r.x,
+                    mIC[i + 1].point.r.y,
+                    mIC[i + 1].point.r.z,
+                    mIC[i + 1].point.t,
+                ).also { x = it }
+                integrator.setSBall(mIC[i + 1].sBall)
+                fastStepper.setSBall(mIC[i + 1].sBall)
+                integrator.setF81(mIC[i + 1].f81)
+                fastStepper.setF81(mIC[i + 1].f81)
+                integrator.accelerate(x)
+
+                vit = mIC[i + 1].point.vit
+                t = mIC[i + 1].point.t
+            }
+        }
+        return BallisticsError.OK
+    }
+
+    override fun calculate(
+        tbegin: Double,
+        tend: Double,
+    ): BallisticsError {
+        if (mic.isEmpty()) {
+            return BallisticsError.EMPTY_NU
+        } else if (mic.size == 1) {
+            calculate(mic.first(), tbegin, tend)
+        } else {
+            calculate(mic.toTypedArray(), tbegin, tend)
+        }
+        return BallisticsError.OK
+    }
+
+    override fun getStepper(): AbstractStepper {
+        return RungeStepper(this)
+    }
+
+    override fun fastStep(t: Double): OrbitalPoint {
+        val dt = t - x[6]
+        fastStepper.step = dt
+        fastStepper.nextStep(x)
+        val p : OrbitalPoint = OrbitalPoint(
+                t,
+                vit,
+                Vector3D(fastStepper.y[3], fastStepper.y[4], fastStepper.y[5]),
+                Vector3D(fastStepper.y[0], fastStepper.y[1], fastStepper.y[2]),
+            )
+        return p
+    }
+
+}

libs/ballistics-lib/src/main/kotlin/ballistics/orbitalPoints/OrbitalPointsTLE.kt

@@ -0,0 +1,70 @@
+package ballistics.orbitalPoints
+
+import ballistics.orbitalPoints.timeStepper.AbstractStepper
+import ballistics.orbitalPoints.timeStepper.TLEStepper
+import ballistics.types.BallisticsError
+import ballistics.types.EarthType
+import ballistics.types.OrbitalPoint
+import ballistics.utils.math.equations.AbstractEquationCalculator
+import ballistics.utils.math.equations.EquationCalculatorSpan
+import kotlin.math.abs
+
+internal class OrbitalPointsTLE(str1: String, str2: String, earthType: EarthType) :
+    AbstractOrbPointsCalculator(earthType) {
+    val stepper = TLEStepper(str1, str2, earthType)
+    private var vit = 0
+
+    override fun calculate(
+        tbegin: Double,
+        tend: Double,
+    ): BallisticsError {
+        points.clear()
+        revolutions.clear()
+
+        val equationCalculator: AbstractEquationCalculator = EquationCalculatorSpan()
+        equationCalculator.delta = 0.001
+
+        vit = stepper.satellite.orbit.orbitNum
+        var t = stepper.baseEpoch.toDouble()
+        var p = stepper.calculate(t)
+        var pn: OrbitalPoint
+        if (abs(p.r.z) < 0.00001 && p.t >= tbegin) {
+            addRevolution(p.t)
+        }
+
+        try {
+            while (t <= tend) {
+                t = t + step
+                pn = stepper.calculate(t)
+                if (p.r.z < 0 && pn.r.z >= 0) {
+                    ++vit
+                    val tvuz = equationCalculator.calculate(t - step, t + step, this::equation)
+                    tvuz?.let {
+                        if (tvuz >= tbegin) {
+                            addRevolution(it)
+                        }
+                    }
+                }
+                if (t >= tbegin) {
+                    pn.vit = vit
+                    points.add(pn)
+                }
+                p = pn
+            }
+        }catch (ex : Exception){
+            return BallisticsError.STEPPER_ERROR
+        }
+
+        return BallisticsError.OK
+    }
+
+    override fun getStepper(): AbstractStepper {
+        return stepper
+    }
+
+    override fun fastStep(t: Double): OrbitalPoint {
+        val op = stepper.calculate(t)
+        op.vit = vit
+        return op
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/orbitalPoints/integrator/AbstractIntegrator.kt

@@ -0,0 +1,38 @@
+package ballistics.orbitalPoints.integrator
+
+import ballistics.orbitalPoints.integrator.modDv.AbstractMDV
+import ballistics.orbitalPoints.integrator.modDv.MDVBars
+import ballistics.orbitalPoints.integrator.modDv.MDVFoto
+import ballistics.orbitalPoints.integrator.modDv.MDVKondor
+import ballistics.types.AstroType
+import ballistics.types.BallisticsError
+import ballistics.types.EarthType
+import ballistics.types.ModDVType
+
+internal abstract class AbstractIntegrator(var mdType: ModDVType, earthType: EarthType) {
+    protected val a = arrayOf(1, 2, 2, 1)
+    protected val prav: AbstractMDV
+    var step = 60.0
+    val y = arrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
+
+    init {
+        when (mdType) {
+            ModDVType.FOTO -> prav = MDVFoto(earthType)
+            ModDVType.METEORM1 -> prav = MDVKondor(earthType, AstroType.ATJ2000, 8, true)
+            ModDVType.METEORM2 -> prav = MDVKondor(earthType, AstroType.ATJ2000, 8, true)
+            ModDVType.KONDOR -> prav = MDVKondor(earthType, AstroType.ATJ2000, 24, false)
+            ModDVType.KONDOR_PROGNOZ -> prav = MDVKondor(earthType, AstroType.ATJ2000, 24, true)
+            ModDVType.BARS -> prav = MDVBars(earthType)
+        }
+    }
+
+    abstract fun nextStep(x: Array<Double>): BallisticsError
+
+    abstract fun accelerate(x: Array<Double>): BallisticsError
+
+    abstract fun setSBall(s: Double)
+
+    fun setF81(f81: Double) {
+        prav.f81 = if (f81 < 25.0) 25.0 else f81
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/orbitalPoints/integrator/IntegratorAdams.kt

@@ -0,0 +1,157 @@
+package ballistics.orbitalPoints.integrator
+
+import ballistics.types.BallisticsError
+import ballistics.types.EarthType
+import ballistics.types.ModDVType
+
+internal class IntegratorAdams(mdType: ModDVType, earthType: EarthType) : AbstractIntegrator(mdType, earthType) {
+    private val rng: IntegratorRK4 = IntegratorRK4(mdType, earthType)
+
+    val a2 =
+        arrayOf(
+            0.11367394179894179894E-1,
+            -0.93840939153439153439E-1,
+            0.343080357142857142857,
+            -0.732035383597883597884,
+            0.1017964616402116402116E+1,
+            -0.1006919642857142857143E+1,
+            0.1156159060846560846560E+1,
+            0.304224537037037037037,
+        )
+
+    val a1 =
+        arrayOf(
+            -0.304224537037037037037,
+            0.2445163690476190476190E+1,
+            -0.8612127976190476190476E+1,
+            0.17379654431216931216931E+2,
+            -0.22027752976190476190476E+2,
+            0.18054538690476190476190E+2,
+            -0.9525206679894179894180E+1,
+            0.3589955357142857142857E+1,
+        )
+
+    var intpoints =
+        arrayOf(
+            arrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
+            arrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
+            arrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
+            arrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
+            arrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
+            arrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
+            arrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
+            arrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
+            arrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
+            arrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
+            arrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
+        )
+
+    override fun accelerate(x: Array<Double>): BallisticsError {
+        for (m in 0..6)
+            intpoints[0][m] = x[m]
+
+        var er = prav.calculate(intpoints[0])
+        if (er != BallisticsError.OK) {
+            return er
+        }
+        for (n in 0..6)
+            intpoints[9][n] = prav.y[n]
+
+        // кратность шага по Рунге шагу по Адамсу
+        val dStep = 2
+        val ihrng = -step / dStep
+        rng.step = ihrng
+        // текущий элемент массива
+        var j = 8
+
+        val fld1 = arrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
+        val fld2 = arrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
+        val fld3 = arrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
+        for (m in 0..6)
+            fld1[m] = x[m]
+
+        for (i in 0..6) {
+            for (q in 0..1) {
+                er = rng.nextStep(fld1)
+                if (er != BallisticsError.OK) {
+                    return er
+                }
+                for (m in 0..6)
+                    fld2[m] = rng.y[m]
+
+                for (m in 0..6)
+                    fld1[m] = fld2[m]
+            }
+
+            er = prav.calculate(fld1)
+            if (er != BallisticsError.OK) {
+                return er
+            }
+            for (m in 0..6)
+                fld3[m] = prav.y[m]
+
+            for (m in 0..6)
+                intpoints[j][m] = fld3[m]
+
+            j--
+        }
+
+        return BallisticsError.OK
+    }
+
+    override fun setSBall(s: Double) {
+        prav.sBall = s
+        rng.setSBall(s)
+    }
+
+    private fun copyPoints(from : Int, to : Int){
+        for (m in 0..6)
+            intpoints[from][m] = intpoints[to][m]
+    }
+
+    override fun nextStep(x: Array<Double>): BallisticsError {
+        for (m in 0..6)
+            intpoints[1][m] = x[m]
+
+        // Экстрaполяция
+        var r: Double
+        for (j in 1..6) {
+            r = 0.0
+            for (k in 1..8)
+                r = r + intpoints[1 + k][j - 1] * a1[k - 1]
+
+            intpoints[1][j - 1] = intpoints[1][j - 1] + r * step
+        }
+        intpoints[1][6] = intpoints[1][6] + step
+
+        var er = prav.calculate(intpoints[1])
+        if (er != BallisticsError.OK) {
+            return er
+        }
+        for (n in 0..6)
+            intpoints[10][n] = prav.y[n]
+
+        copyPoints(1,0)
+
+        // Интерполяция
+        for (j in 1..6) {
+            r = 0.0
+            for (k in 1..8)
+                r = r + intpoints[2 + k][j - 1] * a2[k - 1]
+
+            intpoints[1][j - 1] = intpoints[1][j - 1] + step * r
+        }
+        intpoints[1][6] = intpoints[1][6] + step
+
+        er = prav.calculate(intpoints[1])
+        if (er != BallisticsError.OK)
+            return er
+        for (n in 0..6)
+            intpoints[10][n] = prav.y[n]
+        for (m in 0..6)
+            y[m] = intpoints[1][m]
+        for (j in 0..9)
+            copyPoints(j, j+1)
+        return BallisticsError.OK
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/orbitalPoints/integrator/IntegratorRK4.kt

@@ -0,0 +1,42 @@
+package ballistics.orbitalPoints.integrator
+
+import ballistics.types.BallisticsError
+import ballistics.types.EarthType
+import ballistics.types.ModDVType
+
+internal class IntegratorRK4(mdType: ModDVType, earthType: EarthType) : AbstractIntegrator(mdType, earthType) {
+    override fun accelerate(x: Array<Double>): BallisticsError {
+        return BallisticsError.OK
+    }
+
+    override fun setSBall(s: Double) {
+        prav.sBall = s
+    }
+
+    override fun nextStep(x: Array<Double>): BallisticsError {
+        val step6 = step / 6
+        var r = arrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
+
+        for (i in 0..6) {
+            y[i] = x[i]
+            r[i] = x[i]
+        }
+        for (j in 0..3) {
+            var er = prav.calculate(r)
+            if (er != BallisticsError.OK) {
+                return er
+            }
+
+            for (i in 0..5)
+                y[i] = y[i] + step6 * a[j] * prav.y[i]
+            if (j == 3) {
+                break
+            }
+            for (i in 0..5)
+                r[i] = x[i] + step / a[j + 1] * prav.y[i]
+        }
+        y[6] = x[6] + step
+
+        return BallisticsError.OK
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/orbitalPoints/integrator/modDv/AbstractMDV.kt

@@ -0,0 +1,36 @@
+package ballistics.orbitalPoints.integrator.modDv
+
+import ballistics.types.BallisticsError
+import ballistics.types.EarthType
+import ballistics.utils.earth.AbstractEarth
+import ballistics.utils.earth.getEarth
+
+internal abstract class AbstractMDV(val earthType: EarthType) {
+    protected val hminimum = 160000.0
+    protected var r = 0.0
+    protected var v = 0.0
+    protected var fik = 0.0
+    protected var bca = 0.0
+    protected var oma = 0.0
+    protected var earth: AbstractEarth = getEarth(earthType)
+    var sBall = 0.0
+    var f81 = 100.0
+    var y = arrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
+
+    fun prepare(x: Array<Double>) {
+        var z2 = x[5] * x[5]
+        var r2 = Math.pow(x[3], 2.0) + Math.pow(x[4], 2.0) + z2
+        var v2 = x[2] * x[2] + x[1] * x[1] + x[0] * x[0]
+        r = Math.sqrt(r2)
+        v = Math.sqrt(v2)
+        fik = z2 / r2
+        var d = 5 * fik
+        var c = earth.ekvRadius * earth.ekvRadius / r2 * 1.5 * earth.c20
+        var b = earth.fM / r / r2
+        var a = b * (1 + c * (d - 1))
+        bca = 2 * b * c - a
+        oma = earth.wEarth * earth.wEarth - a
+    }
+
+    abstract fun calculate(x: Array<Double>): BallisticsError
+}

libs/ballistics-lib/src/main/kotlin/ballistics/orbitalPoints/integrator/modDv/MDVBars.kt

@@ -0,0 +1,51 @@
+package ballistics.orbitalPoints.integrator.modDv
+
+import ballistics.types.BallisticsError
+import ballistics.types.EarthType
+import ballistics.utils.atmosphere.getAtm81
+import ballistics.utils.geopotential.Geopotencial1990
+import ballistics.utils.math.Vector3D
+
+internal open class MDVBars(earthType: EarthType) : AbstractMDV(earthType) {
+    protected var geopotential = Geopotencial1990()
+
+    override fun calculate(x: Array<Double>): BallisticsError {
+        prepare(x)
+        y[0] = oma * x[3] + earth.wEarth * x[1] * 2
+        y[1] = oma * x[4] - earth.wEarth * x[0] * 2
+        y[2] = bca * x[5]
+
+        y[3] = x[0]
+        y[4] = x[1]
+        y[5] = x[2]
+
+        // УЧЕТ АТМОСФЕРЫ //
+
+        val h = r - earth.ekvRadius * (1.0 - earth.alphaEllips * fik)
+        if (h < hminimum) {
+            return BallisticsError.H_MINIMUM_ERROR
+        }
+
+        val ro = getAtm81(h)
+
+        val srv = v * ro * sBall
+
+        y[0] -= srv * x[0]
+        y[1] -= srv * x[1]
+        y[2] -= srv * x[2]
+
+        // УЧЕТ АНОМАЛИЙ ГЕОПОТЕНЦИАЛА //
+
+        var xyz = Vector3D(0.0, 0.0, 0.0)
+        var rez: BallisticsError = geopotential.anomkond(x, 16, 16, xyz)
+        if (rez != BallisticsError.OK) {
+            return rez
+        }
+
+        y[0] = y[0] + xyz.x / 1000000.0
+        y[1] = y[1] + xyz.y / 1000000.0
+        y[2] = y[2] + xyz.z / 1000000.0
+
+        return BallisticsError.OK
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/orbitalPoints/integrator/modDv/MDVFoto.kt

@@ -0,0 +1,33 @@
+package ballistics.orbitalPoints.integrator.modDv
+
+import ballistics.types.BallisticsError
+import ballistics.types.EarthType
+import ballistics.utils.atmosphere.getAtm62
+
+internal class MDVFoto(earthType: EarthType) : AbstractMDV(earthType) {
+    override fun calculate(x: Array<Double>): BallisticsError {
+        prepare(x)
+        y[0] = oma * x[3] + earth.wEarth * x[1] * 2
+        y[1] = oma * x[4] - earth.wEarth * x[0] * 2
+        y[2] = bca * x[5]
+
+        y[3] = x[0]
+        y[4] = x[1]
+        y[5] = x[2]
+
+        val h = r - earth.ekvRadius * (1 - earth.alphaEllips * fik)
+        if (h < hminimum) {
+            return BallisticsError.H_MINIMUM_ERROR
+        }
+
+        val ro = getAtm62(h, 6)
+
+        val srv = v * ro * sBall
+
+        y[0] -= srv * x[0]
+        y[1] -= srv * x[1]
+        y[2] -= srv * x[2]
+
+        return BallisticsError.OK
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/orbitalPoints/integrator/modDv/MDVKondor.kt

@@ -0,0 +1,99 @@
+package ballistics.orbitalPoints.integrator.modDv
+
+import ballistics.types.AstroType
+import ballistics.types.BallisticsError
+import ballistics.types.EarthType
+import ballistics.utils.astro.AstronomerJ2000
+import ballistics.utils.atmosphere.Atmosphere2004
+import ballistics.utils.math.Vector3D
+import java.time.LocalDateTime
+import java.time.ZoneOffset
+import kotlin.math.cos
+import kotlin.math.sin
+
+internal open class MDVKondor(
+    earthType: EarthType,
+    var astroType: AstroType,
+    var garmonics: Int,
+    var isProgn: Boolean,
+) : MDVBars(earthType) {
+    protected val amax = 0.6378136E+7
+    protected val alpha = 3.3528037E-3F
+
+    protected var astro = AstronomerJ2000(earthType)
+    protected var atm2004 = Atmosphere2004()
+
+    override fun calculate(x: Array<Double>): BallisticsError {
+        prepare(x)
+        y[0] = oma * x[3] + earth.wEarth * x[1] * 2
+        y[1] = oma * x[4] - earth.wEarth * x[0] * 2
+        y[2] = bca * x[5]
+
+        y[3] = x[0]
+        y[4] = x[1]
+        y[5] = x[2]
+
+        if (!isProgn) {
+            // УЧЕТ АТМОСФЕРЫ //
+
+            var rb1: Double = fik * alpha
+            rb1 = 1.0 - rb1
+            rb1 *= amax
+            val h: Double = r - rb1
+            if (h < hminimum) {
+                return BallisticsError.H_MINIMUM_ERROR
+            }
+
+            var sun = astro.sunCoordinates(x[6])
+            sun = astro.askToGrinvich(sun, x[6]).basis()
+
+            var day = (LocalDateTime.ofEpochSecond(x[6].toLong(), 0, ZoneOffset.UTC).dayOfYear - 1).toDouble()
+
+            var ro = atm2004.atm2004Kav(f81, day, h, x, sun.x, sun.y, sun.z)
+
+            var srv = v * ro * sBall
+
+            y[0] -= srv * x[0]
+            y[1] -= srv * x[1]
+            y[2] -= srv * x[2]
+        }
+        // УЧЕТ АНОМАЛИЙ ГЕОПОТЕНЦИАЛА //
+
+        var xyz = Vector3D(0.0, 0.0, 0.0)
+        var rez: BallisticsError = geopotential.anomkond(x, garmonics, garmonics, xyz)
+        if (rez != BallisticsError.OK) {
+            return rez
+        }
+
+        y[0] = y[0] + xyz.x / 1000000.0
+        y[1] = y[1] + xyz.y / 1000000.0
+        y[2] = y[2] + xyz.z / 1000000.0
+
+        if (!isProgn) {
+            // УЧЕТ Луны и Солнца
+            var tsol = astro.si2000(x[6])
+
+            val cg: Double = cos(tsol)
+            val sg: Double = sin(tsol)
+
+            // Перевод из ОГЭСК в АГЭСК
+
+            // Перевод из ОГЭСК в АГЭСК
+            val x2: Double = (x[3] * cg - x[4] * sg) / 1000.0
+            val y2: Double = (x[3] * sg + x[4] * cg) / 1000.0
+            val z2: Double = x[5] / 1000.0
+
+            var cor = astro.sunMoonCorrection(x[6], x2, y2, z2)
+            var ddx = (cor.x * cg + cor.y * sg) * 1000.0
+            var ddy = (cor.y * cg - cor.x * sg) * 1000.0
+            var ddz = cor.z * 1000.0
+
+            y[0] = y[0] + ddx // Учет влияния Солнца и Луны
+            y[1] = y[1] + ddy
+            y[2] = y[2] + ddz
+        }
+
+        return BallisticsError.OK
+    }
+}
+

libs/ballistics-lib/src/main/kotlin/ballistics/orbitalPoints/timeStepper/AbstractStepper.kt

@@ -0,0 +1,14 @@
+package ballistics.orbitalPoints.timeStepper
+
+import ballistics.types.OrbitalPoint
+
+interface AbstractStepper {
+    fun calculate(t: Double): OrbitalPoint?
+
+    fun calculate(
+        t: Double,
+        p: OrbitalPoint,
+    ): OrbitalPoint?
+
+    fun clear()
+}

libs/ballistics-lib/src/main/kotlin/ballistics/orbitalPoints/timeStepper/RungeStepper.kt

@@ -0,0 +1,119 @@
+package ballistics.orbitalPoints.timeStepper
+
+import ballistics.orbitalPoints.AbstractOrbPointsCalculator
+import ballistics.orbitalPoints.integrator.AbstractIntegrator
+import ballistics.orbitalPoints.integrator.IntegratorRK4
+import ballistics.types.BallisticsError
+import ballistics.types.EarthType
+import ballistics.types.ModDVType
+import ballistics.types.OrbitalPoint
+import ballistics.utils.math.Vector3D
+import kotlin.math.abs
+import kotlin.math.round
+
+class RungeStepper(var mdType: ModDVType, earthType: EarthType) : AbstractStepper {
+
+    private val integrator: AbstractIntegrator = IntegratorRK4(mdType, earthType) // ModDVType.FOTO//BARS
+    var points = mutableListOf<OrbitalPoint>()
+
+    constructor(source: AbstractOrbPointsCalculator) : this(ModDVType.FOTO, source.earthType) {
+        points = source.points
+    }
+
+    constructor(points: MutableList<OrbitalPoint>, earth: EarthType) : this(ModDVType.FOTO, earth) {
+        this.points = points
+    }
+
+    private fun findPoint(t: Double): OrbitalPoint? {
+        val cnt = points.size
+        val dt = t - points.first().t
+        var ind: Int = round(dt / integrator.step).toInt()
+        if (ind >= -1 && ind < cnt + 1) {
+            if (ind < 0) {
+                ind = 0
+            }
+            if (ind >= cnt) {
+                ind = cnt - 1
+            }
+            while ((t - points[ind].t > integrator.step) && (ind < points.size)) {
+                ++ind
+            }
+            if (t - points[ind].t < 0 && ind > 0) {
+                --ind
+            }
+            return points[ind]
+        }
+        return null
+    }
+
+    override fun calculate(t: Double): OrbitalPoint? {
+        val p = findPoint(t)
+        return p?.let {
+            val dt = t - p.t
+            val st = integrator.step
+            integrator.step = dt
+
+            val x = arrayOf(p.v.x, p.v.y, p.v.z, p.r.x, p.r.y, p.r.z, p.t)
+            val r = integrator.nextStep(x)
+            if (r != BallisticsError.OK) {
+                return null
+            }
+
+            var vit = p.vit
+            if (dt >= 0 && p.r.z > 0 && x[5] < 0 && p.v.z > 0) {
+                vit++
+            }
+            if (dt < 0 && p.r.z < 0 && x[5] > 0 && p.v.z > 0) {
+                vit--
+            }
+
+            integrator.step = st
+            OrbitalPoint(
+                integrator.y[6],
+                vit,
+                Vector3D(integrator.y[3], integrator.y[4], integrator.y[5]),
+                Vector3D(integrator.y[0], integrator.y[1], integrator.y[2]),
+            )
+        } ?: run {
+            null
+        }
+    }
+
+    override fun clear() {
+        points.clear()
+    }
+
+    override fun calculate(
+        t: Double,
+        p: OrbitalPoint,
+    ): OrbitalPoint? {
+        val dt = t - p.t
+        if (abs(dt) > 80.0) {
+            return null
+        }
+
+        val st = integrator.step
+        integrator.step = dt
+
+        val x = arrayOf(p.v.x, p.v.y, p.v.z, p.r.x, p.r.y, p.r.z, p.t)
+        if (integrator.nextStep(x) != BallisticsError.OK) {
+            return null
+        }
+
+        var vit = p.vit
+        if (dt >= 0 && p.r.z > 0 && x[5] < 0 && p.v.z > 0) {
+            vit++
+        }
+        if (dt < 0 && p.r.z < 0 && x[5] > 0 && p.v.z > 0) {
+            vit--
+        }
+
+        integrator.step = st
+        return OrbitalPoint(
+            integrator.y[6],
+            vit,
+            Vector3D(integrator.y[3], integrator.y[4], integrator.y[5]),
+            Vector3D(integrator.y[0], integrator.y[1], integrator.y[2]),
+        )
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/orbitalPoints/timeStepper/TLEStepper.kt

@@ -0,0 +1,96 @@
+package ballistics.orbitalPoints.timeStepper
+
+import ballistics.types.EarthType
+import ballistics.types.OrbitalPoint
+import ballistics.utils.astro.AstronomerJ2000
+import ballistics.utils.math.Vector3D
+import org.nstart.dep265.tletools.zeptomoby.core.Julian
+import org.nstart.dep265.tletools.zeptomoby.core.TLE
+import org.nstart.dep265.tletools.zeptomoby.orbit.Satellite
+import java.util.Calendar
+import java.util.GregorianCalendar
+import java.util.TimeZone
+
+class TLEStepper(str1: String, str2: String, earthType: EarthType) : AbstractStepper {
+    val astro = AstronomerJ2000(earthType)
+    val tleParser: TLE = TLE("", str1, str2)
+    val satellite: Satellite = Satellite(tleParser)
+    val baseEpoch: Double
+
+    init {
+        baseEpoch = extractUTCMillis(satellite) / 1000.0 + 10800
+    }
+
+    override fun clear() {
+
+    }
+
+    private fun extractUTCMillis(
+        year: Int,
+        month: Int,
+        dayOfMonth: Double,
+    ): Long {
+        val gc = GregorianCalendar()
+        gc[Calendar.YEAR] = year
+        gc[Calendar.MONTH] = month - 1
+        gc[Calendar.DAY_OF_MONTH] = dayOfMonth.toInt()
+
+        gc.timeZone = TimeZone.getTimeZone("UTC")
+
+        var dfrac = dayOfMonth - dayOfMonth.toLong()
+        dfrac *= 24.0
+
+        gc[Calendar.HOUR_OF_DAY] = dfrac.toInt()
+        dfrac = 60.0 * (dfrac - dfrac.toInt())
+
+        gc[Calendar.MINUTE] = dfrac.toInt()
+        dfrac = 60.0 * (dfrac - dfrac.toInt())
+
+        gc[Calendar.SECOND] = dfrac.toInt()
+        dfrac = 1000.0 * (dfrac - dfrac.toInt())
+
+        gc[Calendar.MILLISECOND] = dfrac.toInt()
+
+        return gc.timeInMillis
+    }
+
+    private fun extractUTCMillis(date: Julian.DateComponent): Long = extractUTCMillis(date.year, date.mon ?: 1, date.dom ?: 1.0)
+
+    private fun extractUTCMillis(sat: Satellite): Long =
+        extractUTCMillis(sat.orbit.epoch.getComponent())
+
+    override fun calculate(t: Double): OrbitalPoint {
+
+        try {
+            val utcMillis = baseEpoch
+            val dt = t - utcMillis
+            val pos = satellite.positionEci(dt / 60)
+            val ask =
+                OrbitalPoint(
+                    t,
+                    satellite.orbit.orbitNum + (dt / satellite.orbit.period).toInt(),
+                    Vector3D(
+                        pos.position.x * 1000.0,
+                        pos.position.y * 1000.0,
+                        pos.position.z * 1000.0,
+                    ),
+                    Vector3D(
+                        pos.velocity.x * 1000.0,
+                        pos.velocity.y * 1000.0,
+                        pos.velocity.z * 1000.0,
+                    ),
+                )
+
+            return astro.askToGrinvich(ask)
+        } catch (ex : Exception){
+            throw ex
+        }
+    }
+
+    override fun calculate(
+        t: Double,
+        p: OrbitalPoint,
+    ): OrbitalPoint? {
+        return calculate(t)
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/KotlinSpecificTools.kt

@@ -0,0 +1,57 @@
+package org.nstart.dep265.tletools.tools
+
+import java.lang.NumberFormatException
+import kotlin.math.min
+
+fun <Double> MutableList<Double>.mapInPlace(transform: (Double) -> Double) {
+    for (i in this.indices)
+        this[i] = transform(this[i])
+}
+
+fun <Double> MutableList<Double>.mapInPlaceIndexed(transform: (Int, Double) -> Double) {
+    for (i in this.indices)
+        this[i] = transform(i, this[i])
+}
+
+fun String.leftJustified(size: Int, char: Char = ' ', truncate: Boolean = false): String {
+    val padLen = size - length
+    val builder = StringBuilder(this.padEnd(size, char))
+
+    if (padLen <= 0 && truncate) {
+            builder.replace(length + padLen, length, "")
+    }
+
+    return builder.toString()
+}
+
+fun String.rightJustified(size: Int, char: Char = ' ', truncate: Boolean = false): String {
+    val padLen = size - length
+    val builder = StringBuilder(this.padStart(size, char))
+
+    if (padLen <= 0 && truncate) {
+            builder.replace(length + padLen, length, "")
+    }
+
+    return builder.toString()
+}
+
+fun String.substringByLength(startIndex: Int, size: Int): String = this.substring(startIndex, min(startIndex + size, this.length))
+
+operator fun <T, V: Enum<V>> Array<T>.get(index: Enum<V>): T = this[index.ordinal]
+operator fun <T, V: Enum<V>> Array<T>.set(index: Enum<V>, value: T) {
+    this[index.ordinal] = value
+}
+
+fun String?.isParsableToNum(): Boolean {
+
+    if (this == null) return false
+
+    try {
+        this.toDouble()
+    }
+    catch (e: NumberFormatException) {
+        return false
+    }
+
+    return true
+}

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/core/Globals.kt

@@ -0,0 +1,74 @@
+package org.nstart.dep265.tletools.zeptomoby.core
+
+import org.nstart.dep265.tletools.zeptomoby.core.iau.IAU
+import org.nstart.dep265.tletools.zeptomoby.core.iau.IAU76
+import org.nstart.dep265.tletools.zeptomoby.core.wgs.WGS
+import org.nstart.dep265.tletools.zeptomoby.core.wgs.WGS72
+import kotlin.math.atan
+import kotlin.math.floor
+import kotlin.math.pow
+import kotlin.math.sqrt
+
+@Suppress("unused")
+object Globals {
+
+    const val pi: Double = Math.PI
+    const val twoPi: Double = 2.0 * org.nstart.dep265.tletools.zeptomoby.core.Globals.pi
+    const val radsPerDeg = org.nstart.dep265.tletools.zeptomoby.core.Globals.pi / 180.0
+
+    const val gm          = 398601.2   // Earth gravitational constant, km^3/sec^2
+    const val geoSyncAlt  = 42241.892  // km
+    const val earthDia    = 12800.0    // km
+    const val daySidereal = 86164.09   // sec
+    const val day24HR     = 86400      // sec
+
+    const val ae          = 1.0
+    var iau: IAU = IAU76
+    var wgs: WGS = WGS72
+        set(value) {
+            field = value
+
+            org.nstart.dep265.tletools.zeptomoby.core.Globals.ck2 = org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.j2 * 0.5
+            org.nstart.dep265.tletools.zeptomoby.core.Globals.ck4 = -3.0 * org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.j4 / 8.0
+            org.nstart.dep265.tletools.zeptomoby.core.Globals.xj3 = org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.j3
+            org.nstart.dep265.tletools.zeptomoby.core.Globals.qo = org.nstart.dep265.tletools.zeptomoby.core.Globals.ae + 120.0 / org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.xkmPer
+            org.nstart.dep265.tletools.zeptomoby.core.Globals.s = org.nstart.dep265.tletools.zeptomoby.core.Globals.ae + 78.0  / org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.xkmPer
+            org.nstart.dep265.tletools.zeptomoby.core.Globals.xke = sqrt(3600.0 * org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.ge / (org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.xkmPer * org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.xkmPer * org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.xkmPer))
+            org.nstart.dep265.tletools.zeptomoby.core.Globals.qoms2t = (org.nstart.dep265.tletools.zeptomoby.core.Globals.qo - org.nstart.dep265.tletools.zeptomoby.core.Globals.s).pow(4)
+        }
+
+    var ck2              = org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.j2 * 0.5
+    var ck4              = -3.0 * org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.j4 * 0.125
+    var xj3              = org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.j3
+    var qo               = org.nstart.dep265.tletools.zeptomoby.core.Globals.ae + 120.0 / org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.xkmPer
+    var s                = org.nstart.dep265.tletools.zeptomoby.core.Globals.ae + 78.0  / org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.xkmPer
+    const val hrPerDay   = 24.0          // Hours per day   (solar)
+    const val minPerDay  = 1440.0        // Minutes per day (solar)
+    const val secPerDay  = 86400.0       // Seconds per day (solar)
+    const val omegaE     = 1.00273790934 // earth rotation per sidereal day
+    var xke              = sqrt(3600.0 * org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.ge /           //sqrt(ge) ER^3/min^2
+                                        (org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.xkmPer * org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.xkmPer * org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.xkmPer))
+    var qoms2t           = (org.nstart.dep265.tletools.zeptomoby.core.Globals.qo - org.nstart.dep265.tletools.zeptomoby.core.Globals.s).pow(4)            //(QO - S)^4 ER^4
+
+    fun sqr(x: Double) = x*x
+
+    fun fmod(numerator: Double, denominator: Double): Double {
+        val tquot = floor(numerator / denominator).toLong()
+        return numerator - tquot * denominator
+    }
+
+    fun fmod2p(arg: Double): Double {
+        val modU = org.nstart.dep265.tletools.zeptomoby.core.Globals.fmod(
+            arg,
+            org.nstart.dep265.tletools.zeptomoby.core.Globals.twoPi
+        )
+        return if (modU < 0.0) modU + org.nstart.dep265.tletools.zeptomoby.core.Globals.twoPi else modU
+    }
+
+    fun acTan(sinX: Double, cosX: Double): Double =
+        if (cosX == 0.0) (if (sinX > 0.0) org.nstart.dep265.tletools.zeptomoby.core.Globals.pi * 0.5 else 1.5 * org.nstart.dep265.tletools.zeptomoby.core.Globals.pi)
+        else (if (cosX > 0.0) atan(sinX/cosX) else org.nstart.dep265.tletools.zeptomoby.core.Globals.pi + atan(sinX/cosX))
+
+    fun rad2deg(r: Double): Double = Math.toDegrees(r)
+    fun deg2rad(d: Double): Double = Math.toRadians(d)
+}

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/core/Julian.kt

@@ -0,0 +1,120 @@
+package org.nstart.dep265.tletools.zeptomoby.core
+
+@Suppress("unused")
+open class Julian {
+
+    object EpochConst {
+        const val EPOCH_JAN0_12H_1900 = 2415020.0 // Dec 31.5 1899 = Dec 31 1899 12h UTC
+        const val EPOCH_JAN1_00H_1900 = 2415020.5 // Jan  1.0 1900 = Jan  1 1900 00h UTC
+        const val EPOCH_JAN1_12H_2000 = 2451545.0 // Jan  1.5 2000 = Jan  1 2000 12h UTC
+    }
+
+    object Static {
+        fun isLeapYear(y: Int): Boolean =
+            (y % 4 == 0) && (y % 100 != 0) || (y % 400 == 0)
+    }
+
+    data class DateComponent(val year: Int, val mon: Int? = null, val dom: Double? = null)
+
+    protected var date: Double = 0.0
+
+    constructor() {
+        initialize(2000, 1.0)
+    }
+
+    constructor(year: Int, day: Double) {
+        initialize(year, day)
+    }
+
+    constructor(year: Int, mon: Int, day: Int, hour: Int, min: Int, sec: Double = 0.0) {
+        val f1 = (275.0 * mon / 9.0).toInt()
+        val f2 = ((mon + 9.0) / 12.0).toInt()
+
+        val n = if (org.nstart.dep265.tletools.zeptomoby.core.Julian.Static.isLeapYear(year)) f1 - f2 + day - 30
+                else f1 - 2 * f2 + day - 30
+
+        val dblDay = n + (hour + (min + sec / 60.0) / 60.0) / 24.0
+        initialize(year, dblDay)
+    }
+
+    fun toGMST(): Double {
+        val ut = org.nstart.dep265.tletools.zeptomoby.core.Globals.fmod(date + 0.5, 1.0)
+        val tu = (fromJan1_12h_2000() - ut) / 36525.0
+        var gmst = 24110.54841 + tu * (8640184.812866 + tu * (0.093104 - tu * 6.2e-06))
+
+        gmst = org.nstart.dep265.tletools.zeptomoby.core.Globals.fmod(
+            gmst + org.nstart.dep265.tletools.zeptomoby.core.Globals.secPerDay * org.nstart.dep265.tletools.zeptomoby.core.Globals.omegaE * ut,
+            org.nstart.dep265.tletools.zeptomoby.core.Globals.secPerDay
+        )
+
+        if (gmst < 0.0) gmst += org.nstart.dep265.tletools.zeptomoby.core.Globals.secPerDay
+        return  org.nstart.dep265.tletools.zeptomoby.core.Globals.twoPi * (gmst / org.nstart.dep265.tletools.zeptomoby.core.Globals.secPerDay)
+    }
+
+    fun toLMST(lon: Double): Double = org.nstart.dep265.tletools.zeptomoby.core.Globals.fmod(
+        toGMST() + lon,
+        org.nstart.dep265.tletools.zeptomoby.core.Globals.twoPi
+    )
+
+    @Suppress("FunctionName")
+    fun fromJan0_12h_1900(): Double = date - org.nstart.dep265.tletools.zeptomoby.core.Julian.EpochConst.EPOCH_JAN0_12H_1900
+    @Suppress("FunctionName")
+    fun fromJan1_00h_1900(): Double = date - org.nstart.dep265.tletools.zeptomoby.core.Julian.EpochConst.EPOCH_JAN1_00H_1900
+    @Suppress("FunctionName")
+    fun fromJan1_12h_2000(): Double = date - org.nstart.dep265.tletools.zeptomoby.core.Julian.EpochConst.EPOCH_JAN1_12H_2000
+
+    fun getComponent(): org.nstart.dep265.tletools.zeptomoby.core.Julian.DateComponent {
+        val jdAdj = date + 0.5
+        val z = jdAdj.toInt()
+        val f = jdAdj - z
+        val alpha = ((z - 1867216.25) / 36524.25).toInt()
+        val a = z + 1 + alpha - (alpha * 0.25).toInt()
+        val b = a + 1524.0
+        val c = ((b - 122.1) / 365.25).toInt()
+        val d = (c * 365.25).toInt()
+        val e = ((b - d) / 30.6001).toInt()
+
+        val dom   = b - d - (e * 30.6001).toInt() + f
+        val month = if (e < 13.5) e - 1 else e - 13
+        val year  = if (month > 2.5) c - 4716 else c - 4715
+
+        return org.nstart.dep265.tletools.zeptomoby.core.Julian.DateComponent(year, month, dom)
+    }
+
+    fun addDay(day: Double) {
+        date += day
+    }
+
+    fun addHour(hr: Double) {
+        date += hr / org.nstart.dep265.tletools.zeptomoby.core.Globals.hrPerDay
+    }
+
+    fun addMin(min: Double) {
+        date += min / org.nstart.dep265.tletools.zeptomoby.core.Globals.minPerDay
+    }
+
+    fun addSec(sec: Double) {
+        date += sec / org.nstart.dep265.tletools.zeptomoby.core.Globals.secPerDay
+    }
+
+    fun spanDay (b: org.nstart.dep265.tletools.zeptomoby.core.Julian): Double = date - b.date
+    fun spanHour(b: org.nstart.dep265.tletools.zeptomoby.core.Julian): Double = spanDay(b) * org.nstart.dep265.tletools.zeptomoby.core.Globals.hrPerDay
+    fun spanMin (b: org.nstart.dep265.tletools.zeptomoby.core.Julian): Double = spanDay(b) * org.nstart.dep265.tletools.zeptomoby.core.Globals.minPerDay
+    fun spanSec (b: org.nstart.dep265.tletools.zeptomoby.core.Julian): Double = spanDay(b) * org.nstart.dep265.tletools.zeptomoby.core.Globals.secPerDay
+
+    protected fun initialize(year: Int, day: Double) {
+        assert((year > 1582) && (year < 3000))
+        assert((day >= 1.0) && (day < 367.0))
+
+        val iyear = year - 1
+
+        val a: Int = iyear / 100
+        val b: Int = 2 - a + a / 4
+
+        val newYears: Double = (365.25 * iyear).toInt() +
+                       (30.6001 * 14).toInt() +
+                       1720994.5 + b
+
+        date = newYears + day
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/core/Site.kt

@@ -0,0 +1,101 @@
+package org.nstart.dep265.tletools.zeptomoby.core
+
+import org.nstart.dep265.tletools.zeptomoby.core.Globals.deg2rad
+import org.nstart.dep265.tletools.zeptomoby.core.Globals.rad2deg
+import org.nstart.dep265.tletools.zeptomoby.core.Globals.sqr
+import org.nstart.dep265.tletools.zeptomoby.core.coord.Geo
+import org.nstart.dep265.tletools.zeptomoby.core.coord.Topo
+import org.nstart.dep265.tletools.zeptomoby.core.eci.EciTime
+import kotlin.math.*
+
+@Suppress("unused")
+open class Site {
+
+    constructor(degLat: Double, degLon: Double, kmAlt: Double, n: String) {
+        name = n
+        geo = Geo(deg2rad(degLat), deg2rad(degLon), kmAlt)
+    }
+
+    constructor(degLat: Double, degLon: Double, kmAlt: Double):
+            this(degLat, degLon, kmAlt, "")
+
+    constructor(g: Geo) {
+        name = ""
+        geo = g
+    }
+
+    override fun toString() = if (name.isEmpty()) "$geo" else "$name $geo"
+
+    fun positionEci(julian: org.nstart.dep265.tletools.zeptomoby.core.Julian): EciTime = EciTime(geo, julian)
+
+    @Deprecated("", ReplaceWith("positionEci(geo, julian)"))
+    fun position(julian: org.nstart.dep265.tletools.zeptomoby.core.Julian): EciTime = EciTime(geo, julian)
+
+    fun lookAngle(time: EciTime): Topo {
+        val date: org.nstart.dep265.tletools.zeptomoby.core.Julian = time.date
+        val eciSite = EciTime(geo, date)
+
+        val vecRgRate = org.nstart.dep265.tletools.zeptomoby.core.Vector(
+            time.velocity.x - eciSite.velocity.x,
+            time.velocity.y - eciSite.velocity.y,
+            time.velocity.z - eciSite.velocity.z
+        )
+
+        val x: Double = time.position.x - eciSite.position.x
+        val y: Double = time.position.y - eciSite.position.y
+        val z: Double = time.position.z - eciSite.position.z
+        val w: Double = sqrt(sqr(x) + sqr(y) + sqr(z))
+        val vecRange = org.nstart.dep265.tletools.zeptomoby.core.Vector(x, y, z, w)
+
+        val theta = date.toLMST(longitudeRad)
+        val sinLat = sin(latitudeRad)
+        val cosLat = cos(latitudeRad)
+        val sinTheta = sin(theta)
+        val cosTheta = cos(theta)
+
+        val topS = sinLat * (cosTheta * vecRange.x + sinTheta * vecRange.y) - cosLat * vecRange.z
+        val topE = -sinTheta * vecRange.x + cosTheta * vecRange.y
+        val topZ = cosLat * (cosTheta * vecRange.x + sinTheta * vecRange.y) + sinLat * vecRange.z
+        var az = atan(-topE / topS)
+
+        if (topS > 0.0) az += org.nstart.dep265.tletools.zeptomoby.core.Globals.pi
+        if (az < 0.0) az += org.nstart.dep265.tletools.zeptomoby.core.Globals.twoPi
+
+        var el = asin(topZ / vecRange.w)
+        val rate = (vecRange.x * vecRgRate.x +
+                vecRange.y * vecRgRate.y +
+                vecRange.z * vecRgRate.z) / vecRange.w
+
+        if (atmosphericCorrection) {
+            val saveEl = el
+
+            el += deg2rad(
+                (1.02 / tan(
+                    deg2rad(rad2deg(el) + 10.3 / (rad2deg(el) + 5.11))
+                )) / 60.0
+            )
+
+            if (el < 0.0) el = saveEl
+            if (el > org.nstart.dep265.tletools.zeptomoby.core.Globals.pi * 0.5) el = org.nstart.dep265.tletools.zeptomoby.core.Globals.pi * 0.5
+        }
+
+        return Topo(az,           // azimuth,   radians
+                    el,           // elevation, radians
+                    vecRange.w,   // range, km
+                    rate          // rate,  km / sec
+        )
+    }
+
+    val latitudeRad: Double
+        get() = geo.latitudeRad
+    val longitudeRad: Double
+        get() = geo.longitudeRad
+    val latitudeDeg: Double
+        get() = geo.latitudeDeg
+    val longitudeDeg: Double
+        get() = geo.longitudeDeg
+
+    val name: String
+    val geo: Geo
+    val atmosphericCorrection = false
+}

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/core/TLE.kt

@@ -0,0 +1,276 @@
+package org.nstart.dep265.tletools.zeptomoby.core
+
+import org.nstart.dep265.tletools.tools.get
+import org.nstart.dep265.tletools.tools.isParsableToNum
+import org.nstart.dep265.tletools.tools.set
+import org.nstart.dep265.tletools.tools.substringByLength
+import org.nstart.dep265.tletools.zeptomoby.core.TLE.ProtectedStatic.convertUnits
+
+@Suppress("unused")
+open class TLE {
+
+    object TLEPartLength {
+        const val TLE_LEN_LINE_DATA      = 69; const val TLE_LEN_LINE_NAME      = 24
+
+        const val TLE1_COL_SATNUM        =  2; const val TLE1_LEN_SATNUM        =  5
+        const val TLE1_COL_INTLDESC_A    =  9; const val TLE1_LEN_INTLDESC_A    =  2
+        const val TLE1_COL_INTLDESC_B    = 11; const val TLE1_LEN_INTLDESC_B    =  3
+        const val TLE1_COL_INTLDESC_C    = 14; const val TLE1_LEN_INTLDESC_C    =  3
+        const val TLE1_COL_EPOCH_A       = 18; const val TLE1_LEN_EPOCH_A       =  2
+        const val TLE1_COL_EPOCH_B       = 20; const val TLE1_LEN_EPOCH_B       = 12
+        const val TLE1_COL_MEANMOTIONDT  = 33; const val TLE1_LEN_MEANMOTIONDT  = 10
+        const val TLE1_COL_MEANMOTIONDT2 = 44; const val TLE1_LEN_MEANMOTIONDT2 =  8
+        const val TLE1_COL_BSTAR         = 53; const val TLE1_LEN_BSTAR         =  8
+        const val TLE1_COL_EPHEMTYPE     = 62; const val TLE1_LEN_EPHEMTYPE     =  1
+        const val TLE1_COL_ELNUM         = 64; const val TLE1_LEN_ELNUM         =  4
+
+        const val TLE2_COL_SATNUM        = 2;  const val TLE2_LEN_SATNUM        =  5
+        const val TLE2_COL_INCLINATION   = 8;  const val TLE2_LEN_INCLINATION   =  8
+        const val TLE2_COL_RAASCENDNODE  = 17; const val TLE2_LEN_RAASCENDNODE  =  8
+        const val TLE2_COL_ECCENTRICITY  = 26; const val TLE2_LEN_ECCENTRICITY  =  7
+        const val TLE2_COL_ARGPERIGEE    = 34; const val TLE2_LEN_ARGPERIGEE    =  8
+        const val TLE2_COL_MEANANOMALY   = 43; const val TLE2_LEN_MEANANOMALY   =  8
+        const val TLE2_COL_MEANMOTION    = 52; const val TLE2_LEN_MEANMOTION    = 11
+        const val TLE2_COL_REVATEPOCH    = 63; const val TLE2_LEN_REVATEPOCH    =  5
+    }
+
+    object Static {
+        fun isValidLine(str: String, num: org.nstart.dep265.tletools.zeptomoby.core.TLE.Lines): Boolean {
+            val testStr = str.trim()
+
+            return if (num == org.nstart.dep265.tletools.zeptomoby.core.TLE.Lines.LINE_ZERO) {
+                testStr.length <= org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE_LEN_LINE_NAME
+            } else {
+                testStr.length == org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE_LEN_LINE_DATA &&
+                str[0].digitToInt() == num.ordinal &&
+                str[1] == ' '
+            }
+        }
+    }
+
+    protected object ProtectedStatic {
+        fun expToAtof(exp: String): String {
+            val colSign = 0
+            val lenSign = 1
+
+            val colMantissa = 1
+            val lenMantissa = 5
+
+            val colExponent = 6
+            val lenExponent = 2
+
+            val strBuilder = StringBuilder()
+            strBuilder.append(exp.substringByLength(colSign, lenSign))
+            strBuilder.append("0.")
+            strBuilder.append(exp.substringByLength(colMantissa, lenMantissa))
+            strBuilder.append("e")
+            strBuilder.append(exp.substringByLength(colExponent, lenExponent).trimStart())
+
+            return strBuilder.toString()
+        }
+
+        fun convertUnits(value: Double, fld: org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields, units: org.nstart.dep265.tletools.zeptomoby.core.TLE.Units): Double =
+            when (fld) {
+                org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_I,
+                org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_RAAN,
+                org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_ARGPER,
+                org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_M -> if (units == org.nstart.dep265.tletools.zeptomoby.core.TLE.Units.U_RAD) org.nstart.dep265.tletools.zeptomoby.core.Globals.deg2rad(
+                    value
+                ) else value
+                else -> value
+            }
+
+        fun checkSum(cs: String): Int =
+            cs.dropLast(1).fold(0) {
+                acc: Int, ch: Char ->
+                when (true) {
+                    ch.isDigit() -> acc + ch.digitToInt()
+                    (ch == '-') -> acc + 1
+                    else -> acc
+                }
+            }
+    }
+
+    private object PrivateStatic {
+        const val strDegrees: String = " degrees"
+        const val strRevsPerDay: String = " revs / day"
+
+        fun getUnits(fld: org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields): String =
+            when (fld) {
+                org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_I,
+                org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_RAAN,
+                org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_ARGPER,
+                org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_M -> org.nstart.dep265.tletools.zeptomoby.core.TLE.PrivateStatic.strDegrees
+                org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_MMOTION -> org.nstart.dep265.tletools.zeptomoby.core.TLE.PrivateStatic.strRevsPerDay
+                else -> ""
+            }
+    }
+
+    enum class Lines {
+        LINE_ZERO,
+        LINE_ONE,
+        LINE_TWO
+    }
+
+    enum class Fields {
+        FLD_NORADNUM,
+        FLD_INTLDESC,
+        FLD_SET,       // TLE set number
+        FLD_EPOCHYEAR, // Epoch: Last two digits of year
+        FLD_EPOCHDAY,  // Epoch: Fractional Julian Day of year
+        FLD_ORBITNUM,  // Orbit at epoch
+        FLD_I,         // Inclination
+        FLD_RAAN,      // R.A. ascending node
+        FLD_E,         // Eccentricity
+        FLD_ARGPER,    // Argument of perigee
+        FLD_M,         // Mean anomaly
+        FLD_MMOTION,   // Mean motion
+        FLD_MMOTIONDT, // First time derivative of mean motion
+        FLD_MMOTIONDT2,// Second time derivative of mean motion
+        FLD_BSTAR,     // BSTAR Drag
+        FLD_LAST       // MUST be last
+    }
+
+    enum class Units {
+        U_RAD,            // radians
+        U_DEG,            // degrees
+        U_NATIVE,         // TLE format native units (no conversion)
+        U_LAST            // MUST be last
+    }
+
+    val header: String
+    val first: String
+    val second: String
+
+    private var field = arrayOfNulls<String>(org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_LAST.ordinal)
+    private var mapCache = mutableMapOf<Int, Double>()
+
+    fun getField(fld: org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields, units: org.nstart.dep265.tletools.zeptomoby.core.TLE.Units = org.nstart.dep265.tletools.zeptomoby.core.TLE.Units.U_NATIVE, bStrUnits: Boolean = false): Pair<Double, String> {
+        assert((0 <= fld.ordinal) && (fld < org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_LAST))
+        assert((0 <= units.ordinal) && (units < org.nstart.dep265.tletools.zeptomoby.core.TLE.Units.U_LAST))
+
+        val strBuilder = StringBuilder(field[fld])
+        if (bStrUnits) strBuilder.append(org.nstart.dep265.tletools.zeptomoby.core.TLE.PrivateStatic.getUnits(fld))
+        strBuilder.trim()
+
+        val key = key(units, fld)
+        val valueInDouble: Double
+
+        if (!mapCache.contains(key)) {
+            val valStr = field[fld.ordinal]
+
+            if (!valStr.isParsableToNum()) valueInDouble = 0.0
+            else {
+                valueInDouble = convertUnits(valStr!!.toDouble(), fld, units)
+                mapCache[key] = valueInDouble
+            }
+        }
+        else valueInDouble = mapCache[key] ?: 0.0
+
+        return Pair(valueInDouble, strBuilder.toString())
+    }
+
+    constructor(h: String, f: String, s: String) {
+        header = h.trimEnd()
+        first = f
+        second = s
+
+        initialize()
+    }
+
+    constructor(tle: org.nstart.dep265.tletools.zeptomoby.core.TLE) {
+        header = tle.header
+        first = tle.first
+        second = tle.second
+
+        field = tle.field
+        mapCache = tle.mapCache
+    }
+
+    protected fun initialize() {
+        if (field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_NORADNUM]?.isNotEmpty() == true) return
+
+        assert(first.isNotEmpty())
+        assert(second.isNotEmpty())
+
+        field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_NORADNUM] = first.substringByLength(
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE1_COL_SATNUM,
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE1_LEN_SATNUM
+        )
+        field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_INTLDESC] = first.substringByLength(
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE1_COL_INTLDESC_A,
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE1_LEN_INTLDESC_A + org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE1_LEN_INTLDESC_B + org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE1_LEN_INTLDESC_C
+        )
+        field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_EPOCHYEAR] = first.substringByLength(
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE1_COL_EPOCH_A,
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE1_LEN_EPOCH_A
+        )
+
+        field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_EPOCHDAY] = first.substringByLength(
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE1_COL_EPOCH_B,
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE1_LEN_EPOCH_B
+        )
+
+        if (first[org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE1_COL_MEANMOTIONDT] == '-') field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_MMOTIONDT] = "-0"
+        else field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_MMOTIONDT] = "0"
+
+        field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_MMOTIONDT] += first.substringByLength(
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE1_COL_MEANMOTIONDT + 1,
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE1_LEN_MEANMOTIONDT
+        )
+
+        field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_MMOTIONDT2] =
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.ProtectedStatic.expToAtof(
+                first.substringByLength(
+                    org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE1_COL_MEANMOTIONDT2,
+                    org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE1_LEN_MEANMOTIONDT2
+                )
+            )
+        field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_BSTAR] =
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.ProtectedStatic.expToAtof(
+                first.substringByLength(
+                    org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE1_COL_BSTAR,
+                    org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE1_LEN_BSTAR
+                )
+            )
+        field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_SET] = first.substringByLength(
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE1_COL_ELNUM,
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE1_LEN_ELNUM
+        ).trimStart()
+
+        field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_I] = second.substringByLength(
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE2_COL_INCLINATION,
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE2_LEN_INCLINATION
+        ).trimStart()
+        field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_RAAN] = second.substringByLength(
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE2_COL_RAASCENDNODE,
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE2_LEN_RAASCENDNODE
+        ).trimStart()
+
+        field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_E] = "0."
+        field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_E] += second.substringByLength(
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE2_COL_ECCENTRICITY,
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE2_LEN_ECCENTRICITY
+        )
+        field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_E] = field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_E.ordinal]!!.trimStart()
+
+        field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_ARGPER] = second.substringByLength(
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE2_COL_ARGPERIGEE,
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE2_LEN_ARGPERIGEE
+        ).trimStart()
+        field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_M] = second.substringByLength(
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE2_COL_MEANANOMALY,
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE2_LEN_MEANANOMALY
+        ).trimStart()
+        field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_MMOTION] = second.substringByLength(
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE2_COL_MEANMOTION,
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE2_LEN_MEANMOTION
+        ).trimStart()
+        field[org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields.FLD_ORBITNUM] = second.substringByLength(
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE2_COL_REVATEPOCH,
+            org.nstart.dep265.tletools.zeptomoby.core.TLE.TLEPartLength.TLE2_LEN_REVATEPOCH
+        ).trimStart()
+    }
+
+    private fun key(u: org.nstart.dep265.tletools.zeptomoby.core.TLE.Units, f: org.nstart.dep265.tletools.zeptomoby.core.TLE.Fields): Int = (u.ordinal * 100) + f.ordinal
+}

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/core/Vector.kt

@@ -0,0 +1,31 @@
+package org.nstart.dep265.tletools.zeptomoby.core
+
+import kotlin.math.abs
+import kotlin.math.acos
+import kotlin.math.sqrt
+
+@Suppress("unused")
+class Vector(
+    var x: Double = 0.0,
+    var y: Double = 0.0,
+    var z: Double = 0.0,
+    var w: Double = 0.0) {
+
+    fun sub(v: Vector) {
+        x -= v.x
+        y -= v.y
+        z -= v.z
+        w -= v.w
+    }
+
+    fun mul(f: Double) {
+        x *= f
+        y *= f
+        z *= f
+        w *= abs(f)
+    }
+
+    fun angle(v: Vector): Double = acos(dot(v) / (magnitude() * v.magnitude()))
+    fun magnitude(): Double = sqrt(x*x + y*y + z*z)
+    fun dot(v: Vector): Double = x*v.x + y*v.y + z*v.z
+}

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/core/coord/Geo.kt

@@ -0,0 +1,78 @@
+package org.nstart.dep265.tletools.zeptomoby.core.coord
+
+import org.nstart.dep265.tletools.zeptomoby.core.Globals
+import org.nstart.dep265.tletools.zeptomoby.core.Globals.acTan
+import org.nstart.dep265.tletools.zeptomoby.core.Globals.fmod
+import org.nstart.dep265.tletools.zeptomoby.core.Globals.rad2deg
+import org.nstart.dep265.tletools.zeptomoby.core.Globals.sqr
+import org.nstart.dep265.tletools.zeptomoby.core.Julian
+import org.nstart.dep265.tletools.zeptomoby.core.Vector
+import org.nstart.dep265.tletools.zeptomoby.core.eci.Eci
+import java.lang.StringBuilder
+import kotlin.math.abs
+import kotlin.math.cos
+import kotlin.math.sin
+import kotlin.math.sqrt
+
+open class Geo {
+
+    constructor(eci: Eci, date: Julian) {
+        val coords = construct(eci.position, fmod((acTan(eci.position.y, eci.position.x) - date.toGMST()), org.nstart.dep265.tletools.zeptomoby.core.Globals.twoPi))
+
+        latitudeRad = coords[0]
+        longitudeRad = coords[1]
+        altitudeKm = coords[2]
+    }
+
+    constructor(latRad: Double, lonRad: Double, altKm: Double) {
+        latitudeRad = latRad
+        longitudeRad = lonRad
+        altitudeKm = altKm
+    }
+
+    var altitudeKm: Double
+    val latitudeRad: Double
+    val longitudeRad: Double
+    val latitudeDeg: Double
+        get() = rad2deg(latitudeRad)
+    val longitudeDeg: Double
+        get() = rad2deg(longitudeRad)
+
+    override fun toString(): String {
+        val strBuilder = StringBuilder()
+        val isNorth = latitudeRad >= 0.0
+        val isEast  = longitudeRad >= 0.0
+
+        strBuilder.append(String.format("%04.3f%c", latitudeDeg, if (isNorth) 'N' else 'S'))
+        strBuilder.append(" ")
+        strBuilder.append(String.format("%05.3f%c", latitudeDeg, if (isEast) 'E' else 'W'))
+        strBuilder.append(" ")
+        strBuilder.append(String.format("%.1fm", altitudeKm * 1000.0))
+
+        return strBuilder.toString()
+    }
+
+    protected fun construct(posEcf: Vector, theta: Double): Array<Double> {
+        var fTheta = fmod(theta, org.nstart.dep265.tletools.zeptomoby.core.Globals.twoPi)
+        if (fTheta  < 0.0) fTheta += org.nstart.dep265.tletools.zeptomoby.core.Globals.twoPi
+
+        val kmSemiMaj = org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.xkmPer
+
+        val r = sqrt(sqr(posEcf.x) + sqr(posEcf.y))
+        val e2 = org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.f * (2.0 - org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.f)
+        var lat = acTan(posEcf.z, r)
+
+        val delta = 1.0e-7
+        var phi: Double
+        var c: Double
+
+        do {
+            phi = lat
+            c = 1.0 / sqrt(1.0 - e2 * sqr(sin(phi)))
+            lat = acTan(posEcf.z + kmSemiMaj * c * e2 * sin(phi), r)
+        }
+        while (abs(lat - phi) > delta)
+
+        return arrayOf(lat, fTheta, r / cos(lat) - kmSemiMaj * c)
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/core/coord/GeoTime.kt

@@ -0,0 +1,26 @@
+package org.nstart.dep265.tletools.zeptomoby.core.coord
+
+import org.nstart.dep265.tletools.zeptomoby.core.Julian
+import org.nstart.dep265.tletools.zeptomoby.core.eci.Eci
+import org.nstart.dep265.tletools.zeptomoby.core.eci.EciTime
+
+@Suppress("unused")
+class GeoTime: Geo {
+
+    constructor(geo: Geo, d: Julian):
+        this(geo.latitudeRad, geo.longitudeRad, geo.altitudeKm, d)
+
+    constructor(latRad: Double, lonRad: Double, altKm: Double, d: Julian): super(latRad, lonRad, altKm) {
+        date = d
+    }
+
+    constructor(eci: Eci, d: Julian): super(eci, d) {
+        date = d
+    }
+
+    constructor(eci: EciTime): super(eci, eci.date) {
+        date = eci.date
+    }
+
+    val date: Julian
+}

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/core/coord/Topo.kt

@@ -0,0 +1,13 @@
+package org.nstart.dep265.tletools.zeptomoby.core.coord
+
+import org.nstart.dep265.tletools.zeptomoby.core.Globals.rad2deg
+
+@Suppress("unused")
+open class Topo(val azimuthRad: Double, val elevationRad: Double,
+                val rangeKm: Double, val rangeRatekmSec: Double) {
+
+    val azimuthDeg: Double
+        get() = rad2deg(azimuthRad)
+    val elevationDeg: Double
+        get() = rad2deg(elevationRad)
+}

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/core/coord/TopoTime.kt

@@ -0,0 +1,17 @@
+package org.nstart.dep265.tletools.zeptomoby.core.coord
+
+import org.nstart.dep265.tletools.zeptomoby.core.Julian
+
+@Suppress("unused")
+class TopoTime: Topo {
+
+    constructor(topo: Topo, d: Julian):
+            this(topo.azimuthRad, topo.elevationRad, topo.rangeKm, topo.rangeRatekmSec, d)
+
+    constructor(azRad: Double, elRad: Double, range: Double ,rangeRate: Double, d: Julian)
+            : super(azRad, elRad, range, rangeRate) {
+        date = d
+    }
+
+    val date: Julian
+}

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/core/eci/Eci.kt

@@ -0,0 +1,57 @@
+package org.nstart.dep265.tletools.zeptomoby.core.eci
+
+import org.nstart.dep265.tletools.zeptomoby.core.Globals
+import org.nstart.dep265.tletools.zeptomoby.core.Globals.sqr
+import org.nstart.dep265.tletools.zeptomoby.core.Julian
+import org.nstart.dep265.tletools.zeptomoby.core.Vector
+import org.nstart.dep265.tletools.zeptomoby.core.coord.Geo
+import kotlin.math.cos
+import kotlin.math.sin
+import kotlin.math.sqrt
+
+open class Eci {
+
+    constructor(pos: Vector, vel: Vector) {
+        position = pos
+        velocity = vel
+    }
+
+    constructor(geo: Geo, date: Julian) {
+        val lat: Double = geo.latitudeRad
+        val lon: Double = geo.longitudeRad
+        val alt: Double = geo.altitudeKm
+
+        val theta: Double = date.toLMST(lon)
+        val c: Double = 1.0 / sqrt(1.0 + org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.f * (org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.f - 2.0) * sqr(sin(lat)))
+        val s: Double = sqr(1.0 - org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.f) * c
+        val achcp: Double = (org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.xkmPer * c + alt) * cos(lat)
+
+        val x = achcp * cos(theta) // km
+        val y = achcp * sin(theta) // km
+        val z = (org.nstart.dep265.tletools.zeptomoby.core.Globals.wgs.xkmPer * s + alt) * sin(lat) // km
+
+        position = Vector(x, y, z,
+            sqrt(
+                sqr(x) + sqr(y) + sqr(z) // range, km
+            )
+        )
+
+        val mFactor: Double = org.nstart.dep265.tletools.zeptomoby.core.Globals.twoPi * (org.nstart.dep265.tletools.zeptomoby.core.Globals.omegaE / org.nstart.dep265.tletools.zeptomoby.core.Globals.secPerDay)
+
+        val vx = -mFactor * position.y  // km / sec
+        val vy = mFactor * position.x   // km / sec
+        val vz = 0.0                    // km / sec
+
+        velocity = Vector(vx, vy, vz,
+            sqrt(
+                sqr(vx) + sqr(vy) // range rate km/sec^2
+            )
+        )
+    }
+
+    val position: Vector
+    val velocity: Vector
+
+    fun scalePosVector(f: Double) = position.mul(f)
+    fun scaleVelVector(f: Double) = velocity.mul(f)
+}

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/core/eci/EciTime.kt

@@ -0,0 +1,28 @@
+package org.nstart.dep265.tletools.zeptomoby.core.eci
+
+import org.nstart.dep265.tletools.zeptomoby.core.Julian
+import org.nstart.dep265.tletools.zeptomoby.core.Vector
+import org.nstart.dep265.tletools.zeptomoby.core.coord.Geo
+import org.nstart.dep265.tletools.zeptomoby.core.coord.GeoTime
+
+@Suppress("unused")
+class EciTime: Eci {
+
+    constructor(eci: Eci, d: Julian): super(eci.position, eci.velocity) {
+        date = d
+    }
+
+    constructor(pos: Vector, vel: Vector, d: Julian): super(pos, vel) {
+        date = d
+    }
+
+    constructor(geo: Geo, d: Julian): super(geo, d) {
+        date = d
+    }
+
+    constructor(geo: GeoTime): super(geo, geo.date) {
+        date = geo.date
+    }
+
+    val date: Julian
+}

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/core/exceptions/DecayException.kt

@@ -0,0 +1,5 @@
+package org.nstart.dep265.tletools.zeptomoby.core.exceptions
+
+import org.nstart.dep265.tletools.zeptomoby.core.Julian
+
+open class DecayException(val decayTime: Julian, val satName: String, msg: String): PropagationException(msg)

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/core/exceptions/PropagationException.kt

@@ -0,0 +1,5 @@
+package org.nstart.dep265.tletools.zeptomoby.core.exceptions
+
+import java.lang.Exception
+
+open class PropagationException(msg: String = ""): Exception(msg)

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/core/iau/IAU.kt

@@ -0,0 +1,13 @@
+package org.nstart.dep265.tletools.zeptomoby.core.iau
+
+interface IAU {
+
+    val au: Double // Astronomical unit (km)
+    val sr: Double // Solar radius (km)      (IAU 76)
+}
+
+object IAU76: IAU {
+
+    override val au: Double = 149597870.0
+    override val sr: Double = 696000.0
+}

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/core/wgs/WGS.kt

@@ -0,0 +1,31 @@
+package org.nstart.dep265.tletools.zeptomoby.core.wgs
+
+interface WGS {
+
+    val xkmPer: Double       // Earth equatorial radius - km
+    val f: Double            // Earth flattening
+    val ge: Double           // Earth gravitational constant
+    val j2: Double           // J2 harmonic
+    val j3: Double           // J3 harmonic
+    val j4: Double           // J4 harmonic
+}
+
+object WGS72: WGS {
+
+    override val xkmPer: Double = 6378.135
+    override val f: Double = 1.0 / 298.26
+    override val ge: Double = 398600.8
+    override val j2: Double = 1.0826158e-3
+    override val j3: Double = -2.53881e-6
+    override val j4: Double = -1.65597e-6
+}
+
+object WGS84: WGS {
+
+    override val xkmPer: Double = 6378.137
+    override val f: Double = 1.0 / 298.257223563
+    override val ge: Double = 398600.4418
+    override val j2: Double = 1.08262998905e-3
+    override val j3: Double = -2.53215306e-6
+    override val j4: Double = -1.61098761e-6
+}

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/orbit/Orbit.kt

@@ -0,0 +1,115 @@
+package org.nstart.dep265.tletools.zeptomoby.orbit
+
+import org.nstart.dep265.tletools.zeptomoby.core.Globals
+import org.nstart.dep265.tletools.zeptomoby.core.Globals.sqr
+import org.nstart.dep265.tletools.zeptomoby.core.Julian
+import org.nstart.dep265.tletools.zeptomoby.core.TLE
+import org.nstart.dep265.tletools.zeptomoby.core.eci.EciTime
+import org.nstart.dep265.tletools.zeptomoby.orbit.norad.NoradBase
+import org.nstart.dep265.tletools.zeptomoby.orbit.norad.NoradSDP4
+import org.nstart.dep265.tletools.zeptomoby.orbit.norad.NoradSGP4
+import kotlin.math.cos
+import kotlin.math.pow
+import kotlin.math.sqrt
+
+@Suppress("unused")
+open class Orbit(t: TLE) {
+
+    val inclination: Double
+    val eccentricity: Double
+    val raan: Double
+    val argPerigee: Double
+    val bStar: Double
+    val drag: Double
+    val meanMotionTle: Double
+    val meanAnomaly: Double
+
+    val epoch: Julian
+    val satId: String
+
+    val satName: String
+        get() = tle.header
+    val tleFirstLine: String
+        get() = tle.first
+    val tleSecondLine: String
+        get() = tle.second
+
+    val semiMajor: Double
+    val semiMinor: Double
+    val meanMotion: Double
+    val perigee: Double
+    val apogee: Double
+    val orbitNum: Int
+    val period: Double
+        get() {
+            if (secPeriod < 0.0) {
+                secPeriod = if (meanMotion == 0.0) 0.0
+                            else Globals.twoPi / meanMotion * 60.0
+            }
+            return secPeriod
+        }
+    val major: Double
+        get() = semiMajor * 2.0
+    val minor: Double
+        get() = semiMinor * 2.0
+
+    private val tle: TLE = t
+    private val noradModel: NoradBase
+    private var secPeriod: Double
+
+    init {
+        satId = tle.getField(TLE.Fields.FLD_NORADNUM).second
+        inclination   = radGet(TLE.Fields.FLD_I).first
+        eccentricity  = tle.getField(TLE.Fields.FLD_E).first
+        raan          = radGet(TLE.Fields.FLD_RAAN).first
+        argPerigee    = radGet(TLE.Fields.FLD_ARGPER).first
+        bStar         = tle.getField(TLE.Fields.FLD_BSTAR).first / Globals.ae
+        drag          = tle.getField(TLE.Fields.FLD_MMOTIONDT).first
+        meanMotionTle = tle.getField(TLE.Fields.FLD_MMOTION).first
+        meanAnomaly   = radGet(TLE.Fields.FLD_M).first
+        var epochYear = tle.getField(TLE.Fields.FLD_EPOCHYEAR).first.toInt()
+        val epochDay = tle.getField(TLE.Fields.FLD_EPOCHDAY).first
+        orbitNum = tle.getField(TLE.Fields.FLD_ORBITNUM).first.toInt()
+        epochYear += if (epochYear < 57) 2000 else 1900
+        epoch = Julian(epochYear, epochDay)
+        secPeriod = -1.0
+        val mm = meanMotionTle
+        val rpmin = mm * Globals.twoPi / Globals.minPerDay
+        val a1 = (Globals.xke / rpmin).pow(2.0 / 3.0)
+        val e = eccentricity
+        val i = inclination
+        val temp = 1.5 * Globals.ck2 * (3.0 * sqr(cos(i)) - 1.0) /
+                (1.0 - e * e).pow(1.5)
+        val delta1 = temp / (a1 * a1)
+        val a0 = a1 *
+                (1.0 - delta1 *
+                        ((1.0 / 3.0) + delta1 *
+                                (1.0 + 134.0 / 81.0 * delta1)))
+        val delta0 = temp / (a0 * a0)
+        meanMotion = rpmin / (1.0 + delta0)
+        semiMajor = a0 / (1.0 - delta0)
+        semiMinor = semiMajor * sqrt(1.0 - (e * e))
+        perigee = Globals.wgs.xkmPer * (semiMajor * (1.0 - e) - Globals.ae)
+        apogee = Globals.wgs.xkmPer * (semiMajor * (1.0 + e) - Globals.ae)
+        noradModel = if (Globals.twoPi / meanMotion >= 225) NoradSDP4(this)
+        else NoradSGP4(this)
+    }
+
+    fun tPlusEpoch(t: Julian) = t.spanSec(epoch)
+
+    @Deprecated("", ReplaceWith("positionEci(mpe)"))
+    fun position(mpe: Double) = positionEci(mpe)
+
+    fun positionEci(mpe: Double): EciTime {
+        val eci = noradModel.position(mpe)
+        val radiusAe: Double = Globals.wgs.xkmPer / Globals.ae
+
+        eci.scalePosVector(radiusAe) // km
+        eci.scaleVelVector(radiusAe * (Globals.minPerDay / 86400)) // km/sec
+
+        return eci
+    }
+
+    protected fun radGet(fld: TLE.Fields) = tle.getField(fld, TLE.Units.U_RAD)
+    protected fun degGet(fld: TLE.Fields) = tle.getField(fld, TLE.Units.U_DEG)
+}

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/orbit/Satellite.kt

@@ -0,0 +1,24 @@
+package org.nstart.dep265.tletools.zeptomoby.orbit
+
+import org.nstart.dep265.tletools.zeptomoby.core.Julian
+import org.nstart.dep265.tletools.zeptomoby.core.TLE
+
+@Suppress("unused")
+class Satellite(tle: TLE, n: String? = null) {
+
+    val name: String
+    val orbit: Orbit
+
+    init {
+        orbit = Orbit(tle)
+        name = n ?: orbit.satName
+    }
+
+    fun positionEci(time: Julian) = positionEci(time.spanMin(orbit.epoch))
+    fun positionEci(mpe: Double) = orbit.positionEci(mpe)
+
+    fun copy(src: Satellite) = Satellite(TLE(src.orbit.satName, src.orbit.tleFirstLine, src.orbit.tleSecondLine), name.ifEmpty { null })
+
+    val tleName: String
+        get() = orbit.satName
+}

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/orbit/norad/NoradBase.kt

@@ -0,0 +1,210 @@
+package org.nstart.dep265.tletools.zeptomoby.orbit.norad
+
+import org.nstart.dep265.tletools.zeptomoby.core.Globals
+import org.nstart.dep265.tletools.zeptomoby.core.Globals.sqr
+import org.nstart.dep265.tletools.zeptomoby.core.Vector
+import org.nstart.dep265.tletools.zeptomoby.core.eci.EciTime
+import org.nstart.dep265.tletools.zeptomoby.core.exceptions.DecayException
+import org.nstart.dep265.tletools.zeptomoby.core.exceptions.PropagationException
+import org.nstart.dep265.tletools.zeptomoby.orbit.Orbit
+import kotlin.math.*
+
+abstract class NoradBase(o: Orbit) {
+
+    abstract fun position(tsince: Double): EciTime
+    abstract fun clone(o: Orbit): NoradBase
+
+    protected fun finalPosition(incl: Double, omega: Double, e: Double, a: Double,
+                                xl: Double, xnode: Double, xn: Double, tsince: Double): EciTime {
+        if ((e * e) > 1.0) throw PropagationException("Error in satellite data")
+
+        val beta = sqrt(1.0 - e * e)
+        val axn  = e * cos(omega)
+        var temp = 1.0 / (a * beta * beta)
+
+        val sinip = sin(orbit.inclination)
+        val cosip = cos(orbit.inclination)
+        val aycof = 0.25 * a3ovk2 * sinip
+        val xlcof = (0.125 * a3ovk2 * sinip * (3.0 + 5.0 * cosip)) /
+                (1.0 + cosip)
+        val xll  = temp * xlcof * axn
+        val aynl = temp * aycof
+        val xlt  = xl + xll
+        val ayn  = e * sin(omega) + aynl
+
+        val e6a = 1.0e-6
+        val capu   = Globals.fmod2p(xlt - xnode)
+
+        var temp2  = capu
+        var temp3  = 0.0
+        var temp4  = 0.0
+        var temp5  = 0.0
+        var temp6  = 0.0
+        var sinepw = 0.0
+        var cosepw = 0.0
+        var fDone  = false
+
+        for (i in 1..10) {
+            if (fDone) break
+
+            sinepw = sin(temp2)
+            cosepw = cos(temp2)
+            temp3 = axn * sinepw
+            temp4 = ayn * cosepw
+            temp5 = axn * cosepw
+            temp6 = ayn * sinepw
+
+            val epw = (capu - temp4 + temp3 - temp2) / (1.0 - temp5 - temp6) + temp2
+
+            if (abs(epw - temp2) <= e6a) fDone = true
+            else temp2 = epw
+        }
+
+        val ecose = temp5 + temp6
+        val esine = temp3 - temp4
+        val elsq  = axn * axn + ayn * ayn
+        temp  = 1.0 - elsq
+        val pl = a * temp
+        val r  = a * (1.0 - ecose)
+        var temp1 = 1.0 / r
+        val rdot  = Globals.xke * sqrt(a) * esine * temp1
+        val rfdot = Globals.xke * sqrt(pl) * temp1
+        temp2 = a * temp1
+        val betal = sqrt(temp)
+        temp3 = 1.0 / (1.0 + betal)
+        val cosu = temp2 * (cosepw - axn + ayn * esine * temp3)
+        val sinu = temp2 * (sinepw - ayn - axn * esine * temp3)
+        val u = Globals.acTan(sinu, cosu)
+        val sin2u = 2.0 * sinu * cosu
+        val cos2u = 2.0 * cosu * cosu - 1.0
+
+        temp  = 1.0 / pl
+        temp1 = Globals.ck2 * temp
+        temp2 = temp1 * temp
+
+        val cosip2 = cosip * cosip
+        val x3thm1 = 3.0 * cosip2 - 1.0
+        val x1mth2 = 1.0 - cosip2
+        val x7thm1 = 7.0 * cosip2 - 1.0
+        val rk = r * (1.0 - 1.5 * temp2 * betal * x3thm1) +
+                0.5 * temp1 * x1mth2 * cos2u
+        val uk = u - 0.25 * temp2 * x7thm1 * sin2u
+        val xnodek = xnode + 1.5 * temp2 * cosio * sin2u
+        val xinck = incl + 1.5 * temp2 * cosio * sinio * cos2u
+        val rdotk = rdot - xn * temp1 * x1mth2 * sin2u
+        val rfdotk = rfdot + xn * temp1 * (x1mth2 * cos2u + 1.5 * x3thm1)
+
+        val sinuk = sin(uk)
+        val cosuk = cos(uk)
+        val sinik = sin(xinck)
+        val cosik = cos(xinck)
+        val sinnok = sin(xnodek)
+        val cosnok = cos(xnodek)
+        val xmx = -sinnok * cosik
+        val xmy = cosnok * cosik
+        val ux  = xmx * sinuk + cosnok * cosuk
+        val uy  = xmy * sinuk + sinnok * cosuk
+        val uz  = sinik * sinuk
+        val vx  = xmx * cosuk - cosnok * sinuk
+        val vy  = xmy * cosuk - sinnok * sinuk
+        val vz  = sinik * cosuk
+
+        val x = rk * ux
+        val y = rk * uy
+        val z = rk * uz
+
+        val vecPos = Vector(x, y, z)
+        val altKm = (vecPos.magnitude() * (Globals.wgs.xkmPer / Globals.ae))
+
+        if (altKm < Globals.wgs.xkmPer) {
+            val decayTime = orbit.epoch
+
+            decayTime.addMin(tsince)
+            throw DecayException(decayTime, "${orbit.satName}#${orbit.satId}", "altitude < xkmper")
+        }
+
+        val xdot = rdotk * ux + rfdotk * vx
+        val ydot = rdotk * uy + rfdotk * vy
+        val zdot = rdotk * uz + rfdotk * vz
+
+        val vecVel = Vector(xdot, ydot, zdot)
+
+        val gmt = orbit.epoch
+        gmt.addMin(tsince)
+
+        return EciTime(vecPos, vecVel, gmt)
+    }
+
+    val orbit: Orbit = o
+
+    val cosio: Double = cos(orbit.inclination); val sinio: Double = sin(orbit.inclination)
+    val betao2: Double; val betao: Double
+    var s4: Double
+        protected set
+    var qoms24: Double
+        protected set
+    val tsi: Double; val eta: Double
+    val eeta: Double; val coef: Double; val coef1: Double
+    val c1: Double; val c3: Double; val c4: Double
+    val a3ovk2: Double; val xmdot: Double; val omgdot: Double
+    val xnodot: Double; val xnodcf: Double; val t2cof: Double
+
+    init {
+        val theta2 = cosio * cosio
+        val x3thm1 = 3.0 * theta2 - 1.0
+        val eosq = sqr(orbit.eccentricity)
+        betao2 = 1.0 - eosq
+        betao  = sqrt(betao2)
+        val rp = orbit.semiMajor * (1.0 - orbit.eccentricity)
+        val perigee = (rp - 1.0) * Globals.wgs.xkmPer
+        s4      = Globals.s
+        qoms24  = Globals.qoms2t
+        if (perigee < 156.0) {
+            s4 = perigee - 78.0
+            if (perigee <= 98.0) s4 = 20.0
+
+            qoms24 = ((120.0 - s4) * Globals.ae / Globals.wgs.xkmPer).pow(4)
+            s4 = s4 / Globals.wgs.xkmPer + Globals.ae
+        }
+        val pinvsq = 1.0 / (sqr(orbit.semiMajor) * sqr(betao2))
+        tsi = 1.0 / (orbit.semiMajor - s4)
+        eta = orbit.semiMajor * orbit.eccentricity * tsi
+        eeta = orbit.eccentricity * eta
+        val etasq = eta * eta
+        val psisq = abs(1.0 - etasq)
+        coef  = qoms24 * tsi.pow(4.0)
+        coef1 = coef / psisq.pow(3.5)
+        val c2 = coef1 * orbit.meanMotion *
+                (orbit.semiMajor * (1.0 + 1.5 * etasq + eeta * (4.0 + etasq)) +
+                        0.75 * Globals.ck2 * tsi / psisq * x3thm1 *
+                        (8.0 + 3.0 * etasq * (8.0 + etasq)))
+        c1 = orbit.bStar * c2
+        a3ovk2 = -Globals.xj3 / Globals.ck2 * Globals.ae.pow(3.0)
+        c3 = coef * tsi * a3ovk2 * orbit.meanMotion * Globals.ae * sinio / orbit.eccentricity
+        val x1mth2 = 1.0 - theta2
+        c4 = 2.0 * orbit.meanMotion * coef1 * orbit.semiMajor * betao2 *
+                (eta * (2.0 + 0.5 * etasq) +
+                        orbit.eccentricity * (0.5 + 2.0 * etasq) -
+                        2.0 * Globals.ck2 * tsi / (orbit.semiMajor * psisq) *
+                        (-3.0 * x3thm1 * (1.0 - 2.0 * eeta + etasq * (1.5 - 0.5 * eeta)) +
+                                0.75 * x1mth2 *
+                                (2.0 * etasq - eeta * (1.0 + etasq)) *
+                                cos(2.0 * orbit.argPerigee)))
+        val theta4 = theta2 * theta2
+        val temp1  = 3.0 * Globals.ck2 * pinvsq * orbit.meanMotion
+        val temp2  = temp1 * Globals.ck2 * pinvsq
+        val temp3  = 1.25 * Globals.ck4 * pinvsq * pinvsq * orbit.meanMotion
+        xmdot = orbit.meanMotion + 0.5 * temp1 * betao * x3thm1 +
+                0.0625 * temp2 * betao *
+                (13.0 - 78.0 * theta2 + 137.0 * theta4)
+        val x1m5th = 1.0 - 5.0 * theta2
+        omgdot = -0.5 * temp1 * x1m5th + 0.0625 * temp2 *
+                (7.0 - 114.0 * theta2 +  395.0 * theta4) +
+                temp3 * (3.0 - 36.0 * theta2 + 49.0 * theta4)
+        val xhdot1 = -temp1 * cosio
+        xnodot = xhdot1 + (0.5 * temp2 * (4.0 - 19.0 * theta2) +
+                2.0 * temp3 * (3.0 - 7.0 * theta2)) * cosio
+        xnodcf = 3.5 * betao2 * xhdot1 * c1
+        t2cof  = 1.5 * c1
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/orbit/norad/NoradSDP4.kt

@@ -0,0 +1,616 @@
+package org.nstart.dep265.tletools.zeptomoby.orbit.norad
+
+import org.nstart.dep265.tletools.zeptomoby.core.Globals
+import org.nstart.dep265.tletools.zeptomoby.core.Globals.acTan
+import org.nstart.dep265.tletools.zeptomoby.core.Globals.fmod2p
+import org.nstart.dep265.tletools.zeptomoby.core.Globals.sqr
+import org.nstart.dep265.tletools.zeptomoby.core.eci.EciTime
+import org.nstart.dep265.tletools.zeptomoby.orbit.Orbit
+import kotlin.math.*
+
+open class NoradSDP4(o: Orbit): NoradBase(o) {
+
+    companion object Constants {
+        const val zes: Double = 0.01675
+        const val zel: Double = 0.05490
+        const val zns: Double = 1.19459e-05
+        const val znl: Double = 1.5835218e-04
+        const val thdt: Double = 4.3752691e-03
+    }
+
+    protected data class DeepSecularVars(val xmdf: Double, val omgadf: Double,
+                               val xnode: Double, val em: Double,
+                               val xinc: Double, val xn: Double)
+
+    protected data class DeepPeriodicsVars(val e: Double, val xinc: Double,
+                                 val omgadf: Double, val xnode: Double,
+                                 val xmam: Double)
+
+    protected data class DeepDotTermsVars(val xndot: Double, val xnddt: Double, val xldot: Double)
+
+    protected var dpE3 = 0.0; protected var dpEe2 = 0.0; protected var dpSe2 = 0.0
+    protected var dpSe3 = 0.0; protected var dpSgh2 = 0.0; protected var dpSgh3 = 0.0
+    protected var dpSgh4 = 0.0; protected var dpSh2 = 0.0; protected var dpSh3 = 0.0
+    protected var dpSi2 = 0.0; protected var dpSi3 = 0.0; protected var dpSl2 = 0.0
+    protected var dpSl3 = 0.0; protected var dpSl4 = 0.0; protected var dpXgh2 = 0.0
+    protected var dpXgh3 = 0.0; protected var dpXgh4 = 0.0; protected var dpXh2 = 0.0
+    protected var dpXh3 = 0.0; protected var dpXi2 = 0.0; protected var dpXi3 = 0.0
+    protected var dpXl2 = 0.0; protected var dpXl3 = 0.0; protected var dpXl4 = 0.0
+    protected val dpZmol: Double; protected val dpZmos: Double
+
+    protected var dpAtime: Double; protected var dpD2201 = 0.0; protected var dpD2211 = 0.0
+    protected var dpD3210 = 0.0; protected var dpD3222 = 0.0; protected var dpD4410 = 0.0
+    protected var dpD4422 = 0.0; protected var dpD5220 = 0.0; protected var dpD5232 = 0.0
+    protected var dpD5421 = 0.0; protected var dpD5433 = 0.0; protected var dpDel1 = 0.0
+    protected var dpDel2 = 0.0; protected var dpDel3 = 0.0; protected var dpSse: Double = 0.0
+    protected var dpSsg: Double = 0.0; protected var dpSsh: Double = 0.0; protected var dpSsi: Double = 0.0
+    protected var dpSsl: Double = 0.0; protected val dpStep2: Double; protected val dpStepn: Double
+    protected val dpStepp: Double; protected val dpThgr: Double; protected val dpXfact: Double
+    protected var dpXlamo: Double = 0.0; protected var dpXli: Double; protected var dpXni: Double
+
+    var gpReso: Boolean
+    var gpSync: Boolean
+
+    init {
+        val sinArg = sin(orbit.argPerigee)
+        val cosArg = cos(orbit.argPerigee)
+        val eqsq   = sqr(orbit.eccentricity)
+
+        val jd = orbit.epoch
+        dpThgr = jd.toGMST()
+
+        val eq = orbit.eccentricity
+        val aqnv = 1.0 / orbit.semiMajor
+        val xmao = orbit.meanAnomaly
+        val xpidot = omgdot + xnodot
+        val sinQ = sin(orbit.raan)
+        val cosQ = cos(orbit.raan)
+
+        val day = jd.fromJan0_12h_1900()
+        val dpiXnodce = 4.5236020 - 9.2422029e-4 * day
+        val dpiStem = sin(dpiXnodce)
+        val dpiCtem = cos(dpiXnodce)
+        val dpiZcosil = 0.91375164 - 0.03568096 * dpiCtem
+        val dpiZsinil = sqrt(1.0 - dpiZcosil * dpiZcosil)
+        val dpiZsinhl = 0.089683511 *dpiStem / dpiZsinil
+        val dpiZcoshl = sqrt(1.0 - dpiZsinhl * dpiZsinhl)
+        val dpiC = 4.7199672 + 0.22997150 * day
+        val dpiGam = 5.8351514 + 0.0019443680 * day
+
+        dpZmol = fmod2p(dpiC - dpiGam)
+
+        var dpiZx = 0.39785416 * dpiStem / dpiZsinil
+        val dpiZy = dpiZcoshl * dpiCtem + 0.91744867 * dpiZsinhl * dpiStem
+
+        dpiZx = acTan(dpiZx,dpiZy) + dpiGam - dpiXnodce
+
+        val dpiZcosgl = cos(dpiZx)
+        val dpiZsingl = sin(dpiZx)
+
+        dpZmos   = fmod2p(6.2565837 + 0.017201977 * day)
+
+        val zcosis = 0.91744867
+        val zsinis = 0.39785416
+        val c1ss   = 2.9864797e-06
+        val zsings = -0.98088458
+        val zcosgs =  0.1945905
+
+        var zcosg = zcosgs
+        var zsing = zsings
+        var zcosi = zcosis
+        var zsini = zsinis
+        var zcosh = cosQ
+        var zsinh = sinQ
+        var cc  = c1ss
+        var zn  = zns
+        var ze  = zes
+        val xnoi = 1.0 / orbit.meanMotion
+
+        var se  = 0.0; var si = 0.0; var sl = 0.0
+        var sgh = 0.0; var sh = 0.0
+
+        for (pass in 1..2) {
+            val a1 = zcosg * zcosh + zsing * zcosi * zsinh
+            val a3 = -zsing * zcosh + zcosg * zcosi * zsinh
+            val a7 = -zcosg * zsinh + zsing * zcosi * zcosh
+            val a8 = zsing * zsini
+            val a9 = zsing * zsinh + zcosg * zcosi * zcosh
+            val a10 = zcosg * zsini
+
+            val a2 = cosio * a7 + sinio * a8
+            val a4 = cosio * a9 + sinio * a10
+            val a5 = -sinio * a7 + cosio * a8
+            val a6 = -sinio * a9 + cosio * a10
+            val x1 = a1 * cosArg + a2 * sinArg
+            val x2 = a3 * cosArg + a4 * sinArg
+            val x3 = -a1 * sinArg + a2 * cosArg
+            val x4 = -a3 * sinArg + a4 * cosArg
+            val x5 = a5 * sinArg
+            val x6 = a6 * sinArg
+            val x7 = a5 * cosArg
+            val x8 = a6 * cosArg
+            val z31 = 12.0 * x1 * x1 - 3.0 * x3 * x3
+            val z32 = 24.0 * x1 * x2 - 6.0 * x3 * x4
+            val z33 = 12.0 * x2 * x2 - 3.0 * x4 * x4
+            var z1 = 3.0 * (a1 * a1 + a2 * a2) + z31 * eqsq
+            var z2 = 6.0 * (a1 * a3 + a2 * a4) + z32 * eqsq
+            var z3 = 3.0 * (a3 * a3 + a4 * a4) + z33 * eqsq
+            val z11 = -6.0 * a1 * a5 + eqsq*(-24.0 * x1 * x7 - 6.0 * x3 * x5)
+            val z12 = -6.0 * (a1 * a6 + a3 * a5) +
+                    eqsq * (-24.0 * (x2 * x7 + x1 * x8) - 6.0 * (x3 * x6 + x4 * x5))
+            val z13 = -6.0 * a3 * a6 + eqsq * (-24.0 * x2 * x8 - 6.0 * x4 * x6)
+            val z21 = 6.0 * a2 * a5 + eqsq * (24.0 * x1 * x5 - 6.0 * x3 * x7)
+            val z22 = 6.0*(a4 * a5 + a2 * a6) +
+                    eqsq * (24.0 * (x2 * x5 + x1 * x6) - 6.0 * (x4 * x7 + x3 * x8))
+            val z23 = 6.0 * a4 * a6 + eqsq*(24.0 * x2 * x6 - 6.0 * x4 * x8)
+            z1 += z1 + betao2 * z31
+            z2 += z2 + betao2 * z32
+            z3 += z3 + betao2 * z33
+            val s3 = cc * xnoi
+            val s2 = -0.5 * s3 / betao
+            val s4 = s3 * betao
+            val s1 = -15.0 * eq * s4
+            val s5 = x1 * x3 + x2 * x4
+            val s6 = x2 * x3 + x1 * x4
+            val s7 = x2 * x4 - x1 * x3
+            se = s1 * zn * s5
+            si = s2 * zn * (z11 + z13)
+            sl = -zn * s3 * (z1 + z3 - 14.0 - 6.0 * eqsq)
+            sgh = s4 * zn * (z31 + z33 - 6.0)
+            sh  = -zn * s2 * (z21 + z23)
+
+            if (orbit.inclination < 5.2359877e-2) sh = 0.0
+
+            dpEe2 =  2.0 * s1 * s6
+            dpE3  =  2.0 * s1 * s7
+            dpXi2 =  2.0 * s2 * z12
+            dpXi3 =  2.0 * s2 * (z13 - z11)
+            dpXl2 = -2.0 * s3 * z2
+            dpXl3 = -2.0 * s3 * (z3 - z1)
+            dpXl4 = -2.0 * s3 * (-21.0 - 9.0 * eqsq) * ze
+            dpXgh2 = 2.0 * s4 * z32
+            dpXgh3 = 2.0 * s4 * (z33 - z31)
+            dpXgh4 = -18.0 * s4 * ze
+            dpXh2 = -2.0 * s2 * z22
+            dpXh3 = -2.0 * s2 * (z23 - z21)
+
+            if (pass == 1) {
+                dpSse = se
+                dpSsi = si
+                dpSsl = sl
+                dpSsh = sh / sinio
+                dpSsg = sgh - cosio * dpSsh
+                dpSe2 = dpEe2
+                dpSi2 = dpXi2
+                dpSl2 = dpXl2
+                dpSgh2 = dpXgh2
+                dpSh2 = dpXh2
+                dpSe3 = dpE3
+                dpSi3 = dpXi3
+                dpSl3 = dpXl3
+                dpSgh3 = dpXgh3
+                dpSh3 = dpXh3
+                dpSl4 = dpXl4
+                dpSgh4 = dpXgh4
+                zcosg = dpiZcosgl
+                zsing = dpiZsingl
+                zcosi = dpiZcosil
+                zsini = dpiZsinil
+                zcosh = dpiZcoshl * cosQ + dpiZsinhl * sinQ
+                zsinh = sinQ * dpiZcoshl - cosQ * dpiZsinhl
+                zn = znl
+
+                val c1l = 4.7968065e-07
+                cc = c1l
+                ze = zel
+            }
+        }
+
+        dpSse += se
+        dpSsi += si
+        dpSsl += sl
+        dpSsg = dpSsg + sgh - cosio / sinio * sh
+        dpSsh += sh / sinio
+
+        gpReso = false
+        gpSync = false
+
+        val g310: Double
+        val f220: Double
+        var bfact = 0.0
+
+        if ((orbit.meanMotion > 0.0034906585) &&
+            (orbit.meanMotion < 0.0052359877)) {
+            gpReso = true
+            gpSync = true
+
+            val g200 = 1.0 + eqsq * (-2.5 + 0.8125 * eqsq)
+            g310 = 1.0 + 2.0 * eqsq
+
+            val g300 = 1.0 + eqsq * (-6.0 + 6.60937 * eqsq)
+            f220 = 0.75 * (1.0 + cosio) * (1.0 + cosio)
+
+            val f311 = 0.9375 * sinio * sinio * (1.0 + 3 * cosio) - 0.75 * (1.0 + cosio)
+            var f330 = 1.0 + cosio
+
+            val q22 = 1.7891679e-06
+            val q31 = 2.1460748e-06
+            val q33 = 2.2123015e-07
+
+            f330 *= 1.875 * f330 * f330
+            dpDel1 = 3.0 * orbit.meanMotion * orbit.meanMotion * aqnv * aqnv
+            dpDel2 = 2.0 * dpDel1 * f220 * g200 * q22
+            dpDel3 = 3.0 * dpDel1 * f330 * g300 * q33 * aqnv
+            dpDel1 *= f311 * g310 * q31 * aqnv
+            dpXlamo = xmao + orbit.raan + orbit.argPerigee - dpThgr
+            bfact = xmdot + xpidot - thdt
+            bfact += dpSsl + dpSsg + dpSsh
+        }
+        else if ((orbit.meanMotion >= 8.26e-03) &&
+                 (orbit.meanMotion <= 9.24e-03) &&
+                 (eq >= 0.5)) {
+            gpReso = true
+
+            val eoc  = eq * eqsq
+            val g201 = -0.306 - (eq - 0.64) * 0.440
+
+            val g211: Double; val g322: Double
+            val g410: Double; val g422: Double
+            val g520: Double
+
+            if (eq <= 0.65) {
+                g211 =    3.616 -  13.247  * eq + 16.290   * eqsq
+                g310 =  -19.302 +  117.390 * eq - 228.419  * eqsq + 156.591  * eoc
+                g322 = -18.9068 + 109.7927 * eq - 214.6334 * eqsq + 146.5816 * eoc
+                g410 =  -41.122 +  242.694 * eq - 471.094  * eqsq +  313.953 * eoc
+                g422 = -146.407 +  841.880 * eq - 1629.014 * eqsq + 1083.435 * eoc
+                g520 = -532.114 + 3017.977 * eq - 5740.0   * eqsq + 3708.276 * eoc
+            }
+            else {
+                g211 =   -72.099 +  331.819 * eq -  508.738 * eqsq +  266.724 * eoc
+                g310 =  -346.844 + 1582.851 * eq - 2415.925 * eqsq + 1246.113 * eoc
+                g322 =  -342.585 + 1554.908 * eq - 2366.899 * eqsq + 1215.972 * eoc
+                g410 = -1052.797 + 4758.686 * eq - 7193.992 * eqsq + 3651.957 * eoc
+                g422 = -3581.69  + 16178.11 * eq - 24462.77 * eqsq + 12422.52 * eoc
+
+                g520 = if (eq <= 0.715) 1464.74 - 4664.75 * eq + 3763.64 * eqsq
+                       else -5149.66 + 29936.92 * eq - 54087.36 * eqsq + 31324.56 * eoc
+
+            }
+
+            val g533: Double
+            val g521: Double
+            val g532: Double
+
+            if (eq < 0.7) {
+                g533 = -919.2277  + 4988.61   * eq - 9064.77   * eqsq + 5542.21  * eoc
+                g521 = -822.71072 + 4568.6173 * eq - 8491.4146 * eqsq + 5337.524 * eoc
+                g532 = -853.666   + 4690.25   * eq - 8624.77   * eqsq + 5341.4   * eoc
+            }
+            else {
+                g533 = -37995.78  + 161616.52 * eq - 229838.2  * eqsq + 109377.94 * eoc
+                g521 = -51752.104 + 218913.95 * eq - 309468.16 * eqsq + 146349.42 * eoc
+                g532 = -40023.88  + 170470.89 * eq - 242699.48 * eqsq + 115605.82 * eoc
+            }
+
+            val sini2  = sqr(sinio)
+            val theta2 = sqr(cosio)
+
+            f220 = 0.75 * (1.0 + 2.0 * cosio + theta2)
+
+            val root22 = 1.7891679e-06
+            val root32 = 3.7393792e-07
+            val root44 = 7.3636953e-09
+            val root52 = 1.1428639e-07
+            val root54 = 2.1765803e-09
+
+            val f221 = 1.5 * sini2
+            val f321 =  1.875 * sinio * (1.0 - 2.0 * cosio - 3.0 * theta2)
+            val f322 = -1.875 * sinio * (1.0 + 2.0 * cosio - 3.0 * theta2)
+            val f441 = 35.0 * sini2 * f220
+            val f442 = 39.3750 * sini2 * sini2
+            val f522 = 9.84375 * sinio * (sini2 * (1.0 - 2.0 * cosio - 5.0 * theta2) +
+            0.33333333*(-2.0 + 4.0 * cosio + 6.0 * theta2))
+            val f523 = sinio * (4.92187512 * sini2 * (-2.0 - 4.0 * cosio + 10.0 * theta2) +
+                    6.56250012 * (1.0 + 2.0 * cosio - 3.0 * theta2))
+            val f542 = 29.53125 * sinio * ( 2.0 - 8.0 * cosio + theta2 * (-12.0 + 8.0 * cosio + 10.0 * theta2))
+            val f543 = 29.53125 * sinio * (-2.0 - 8.0 * cosio + theta2 * ( 12.0 + 8.0 * cosio - 10.0 * theta2))
+            val xno2 = orbit.meanMotion * orbit.meanMotion
+            val ainv2 = aqnv * aqnv
+            var temp1 = 3.0 * xno2 * ainv2
+            var temp  = temp1 * root22
+
+            dpD2201 = temp * f220 * g201
+            dpD2211 = temp * f221 * g211
+            temp1 *= aqnv
+            temp = temp1 * root32
+            dpD3210 = temp * f321 * g310
+            dpD3222 = temp * f322 * g322
+            temp1 *= aqnv
+            temp = 2.0 * temp1 * root44
+            dpD4410 = temp * f441 * g410
+            dpD4422 = temp * f442 * g422
+            temp1 *= aqnv
+            temp  = temp1 * root52
+            dpD5220 = temp * f522 * g520
+            dpD5232 = temp * f523 * g532
+            temp = 2.0 * temp1 * root54
+            dpD5421 = temp * f542 * g521
+            dpD5433 = temp * f543 * g533
+            dpXlamo = xmao + orbit.raan + orbit.raan - dpThgr - dpThgr
+            bfact = xmdot + xnodot + xnodot - thdt - thdt
+            bfact += dpSsl + dpSsh + dpSsh
+        }
+
+        if (gpReso || gpSync) {
+            dpXfact = bfact - orbit.meanMotion
+            dpXli   = dpXlamo
+            dpXni   = orbit.meanMotion
+            dpAtime = 0.0
+            dpStepp = 720.0
+            dpStepn = -720.0
+            dpStep2 = 259200.0
+        }
+        else {
+            dpXfact = 0.0
+            dpXli   = 0.0
+            dpXni   = 0.0
+            dpAtime = 0.0
+            dpStepp = 0.0
+            dpStepn = 0.0
+            dpStep2 = 0.0
+        }
+    }
+
+    protected fun deepSecular(xmdf: Double, omgadf: Double, xnode: Double,
+                              xn: Double, tsince: Double): DeepSecularVars {
+        var xmdfLocal   = xmdf + dpSsl * tsince
+        var omgadfLocal = omgadf + dpSsg * tsince
+        var xnodeLocal  = xnode + dpSsh * tsince
+        val emmLocal    = orbit.eccentricity + dpSse * tsince
+        var xinccLocal  = orbit.inclination  + dpSsi * tsince
+
+        if (xinccLocal < 0.0) {
+            xinccLocal *= -1
+            xnodeLocal += Globals.pi
+            omgadfLocal -= Globals.pi
+        }
+
+        if (!gpReso) return DeepSecularVars(xmdfLocal, omgadfLocal, xnodeLocal,
+                                             emmLocal, xinccLocal, xn)
+
+        val xnddt: Double
+        val xndot: Double
+        val xldot: Double
+        var delt  = 0.0
+        var fDone = false
+
+        while (!fDone) {
+            if ((dpAtime == 0.0) ||
+                ((tsince >= 0.0) && (dpAtime <  0.0)) ||
+                ((tsince <  0.0) && (dpAtime >= 0.0))) {
+                delt = if (tsince < 0) dpStepn else dpStepp
+
+                dpAtime = 0.0
+                dpXni = orbit.meanMotion
+                dpXli = dpXlamo
+
+                fDone = true
+            }
+            else {
+                if (abs(tsince) < abs(dpAtime)) {
+                    delt = if (tsince >= 0.0) dpStepn else dpStepp
+                    deepCalcIntegrator(delt)
+                }
+                else {
+                    delt = if (tsince > 0.0) dpStepp else dpStepn
+                    fDone = true
+                }
+            }
+        }
+
+        while (abs(tsince - dpAtime) >= dpStepp) deepCalcIntegrator(delt)
+
+        val ft = tsince - dpAtime
+
+        val vars = deepCalcDotTerms()
+        xndot = vars.xndot
+        xnddt = vars.xnddt
+        xldot = vars.xldot
+
+        val xnnLocal = dpXni + xndot * ft + xnddt * ft * ft * 0.5
+        val xl = dpXli + xldot * ft + xndot * ft * ft * 0.5
+        val temp = -xnodeLocal + dpThgr + tsince * thdt
+
+        xmdfLocal = xl - omgadfLocal + temp
+        if (!gpSync) xmdfLocal = xl + temp + temp
+
+        return DeepSecularVars(xmdfLocal, omgadfLocal, xnodeLocal, emmLocal, xinccLocal, xnnLocal)
+    }
+
+    protected fun deepCalcDotTerms(): DeepDotTermsVars {
+        val fasx2 = 0.13130908
+        val fasx4 = 2.8843198
+        val fasx6 = 0.37448087
+
+        val xndotLocal: Double
+        var xnddtLocal: Double
+
+        if (gpSync) {
+            xndotLocal = dpDel1 * sin(dpXli - fasx2) +
+                dpDel2 * sin(2.0 * (dpXli - fasx4)) +
+                dpDel3 * sin(3.0 * (dpXli - fasx6))
+            xnddtLocal = dpDel1 * cos(dpXli - fasx2) +
+                2.0 * dpDel2 * cos(2.0 * (dpXli - fasx4)) +
+                3.0 * dpDel3 * cos(3.0 * (dpXli - fasx6))
+        }
+        else {
+            val g22 = 5.7686396
+            val g32 = 0.95240898
+            val g44 = 1.8014998
+            val g52 = 1.0508330
+            val g54 = 4.4108898
+
+            val xomi  = orbit.argPerigee + omgdot * dpAtime
+            val x2omi = xomi + xomi
+            val x2li  = dpXli + dpXli
+
+            xndotLocal = dpD2201 * sin(x2omi + dpXli - g22) +
+                dpD2211 * sin(dpXli - g22)         +
+                dpD3210 * sin( xomi + dpXli - g32) +
+                dpD3222 * sin(-xomi + dpXli - g32) +
+                dpD4410 * sin(x2omi + x2li - g44)   +
+                dpD4422 * sin(x2li - g44)           +
+                dpD5220 * sin( xomi + dpXli - g52) +
+                dpD5232 * sin(-xomi + dpXli - g52) +
+                dpD5421 * sin( xomi + x2li - g54)   +
+                dpD5433 * sin(-xomi + x2li - g54)
+
+            xnddtLocal = dpD2201 * cos(x2omi + dpXli - g22) +
+                dpD2211 * cos(dpXli - g22)         +
+                dpD3210 * cos( xomi + dpXli - g32) +
+                dpD3222 * cos(-xomi + dpXli - g32) +
+                dpD5220 * cos( xomi + dpXli - g52) +
+                dpD5232 * cos(-xomi + dpXli - g52) +
+                2.0 * (dpD4410 * cos(x2omi + x2li - g44) +
+                dpD4422 * cos(x2li - g44)         +
+                dpD5421 * cos( xomi + x2li - g54) +
+                dpD5433 * cos(-xomi + x2li - g54))
+        }
+
+        val xldotLocal = dpXni + dpXfact
+        xnddtLocal *= xldotLocal
+
+        return DeepDotTermsVars(xndotLocal, xnddtLocal, xldotLocal)
+    }
+
+    protected fun deepCalcIntegrator(delt: Double): DeepDotTermsVars {
+        val vars = deepCalcDotTerms()
+
+        dpXli += vars.xldot * delt + vars.xndot * dpStep2
+        dpXni += vars.xndot * delt + vars.xnddt * dpStep2
+        dpAtime += delt
+
+        return vars
+    }
+
+    protected fun deepPeriodics(e: Double, xinc: Double,
+                                omgadf: Double, xnode: Double,
+                                xmam: Double, tsince: Double): DeepPeriodicsVars{
+        val sinis = sin(xinc)
+        val cosis = cos(xinc)
+
+        val sghs: Double
+        val shs: Double
+        val sh1: Double
+        val pe: Double
+        val pinc: Double
+        val pl: Double
+        val sghl: Double
+
+        var zm = dpZmos + zns * tsince
+        var zf = zm + 2.0 * zes * sin(zm)
+        var sinzf = sin(zf)
+        var f2  = 0.5 * sinzf * sinzf - 0.25
+        var f3  = -0.5 * sinzf * cos(zf)
+        val ses = dpSe2 * f2 + dpSe3 * f3
+        val sis = dpSi2 * f2 + dpSi3 * f3
+        val sls = dpSl2 * f2 + dpSl3 * f3 + dpSl4 * sinzf
+
+        sghs = dpSgh2 * f2 + dpSgh3 * f3 + dpSgh4 * sinzf
+        shs  = dpSh2  * f2 + dpSh3  * f3
+        zm = dpZmol + znl * tsince
+        zf = zm + 2.0 * zel * sin(zm)
+        sinzf = sin(zf)
+        f2 = 0.5 * sinzf * sinzf - 0.25
+        f3 = -0.5 * sinzf * cos(zf)
+
+        val sel  = dpEe2 * f2 + dpE3  * f3
+        val sil  = dpXi2 * f2 + dpXi3 * f3
+        val sll  = dpXl2 * f2 + dpXl3 * f3 + dpXl4 * sinzf
+
+        sghl = dpXgh2 * f2 + dpXgh3 * f3 + dpXgh4 * sinzf
+        sh1  = dpXh2  * f2 + dpXh3  * f3
+        pe   = ses + sel
+        pinc = sis + sil
+        pl   = sls + sll
+
+        var pgh = sghs + sghl
+        var ph = shs + sh1
+
+        val xinccLocal = xinc + pinc
+        val eLocal  = e + pe
+        val omgadfLocal: Double
+        val xnodeLocal: Double
+        val xmamLocal: Double
+
+        if (orbit.inclination >= 0.2) {
+            ph /= sinio
+            pgh -= cosio * ph
+            omgadfLocal = omgadf + pgh
+            xnodeLocal = xnode + ph
+            xmamLocal = xmam + pl
+        }
+        else {
+            val sinok = sin(xnode)
+            val cosok = cos(xnode)
+            var alfdp = sinis * sinok
+            var betdp = sinis * cosok
+            val dalf  =  ph * cosok + pinc * cosis * sinok
+            val dbet  = -ph * sinok + pinc * cosis * cosok
+
+            alfdp += dalf
+            betdp += dbet
+
+            var xls = xmam + omgadf + cosis * xnode
+            val dls = pl + pgh - pinc * xnode * sinis
+
+            xls += dls
+            xnodeLocal  = acTan(alfdp, betdp)
+            xmamLocal   = xmam + pl
+            omgadfLocal = xls - xmamLocal - cos(xinccLocal) * xnodeLocal
+        }
+
+        return DeepPeriodicsVars(eLocal, xinccLocal, omgadfLocal, xnodeLocal, xmamLocal)
+    }
+
+    override fun position(tsince: Double): EciTime {
+        var xmdf = orbit.meanAnomaly + xmdot  * tsince
+        var omgadf = orbit.argPerigee  + omgdot * tsince
+        val xnoddf = orbit.raan + xnodot * tsince
+        val tsq = tsince * tsince
+        var xnode = xnoddf + xnodcf * tsq
+        val tempa = 1.0 - c1 * tsince
+        val tempe = orbit.bStar * c4 * tsince
+        val templ = t2cof * tsq
+        var xn = orbit.meanMotion
+        val em: Double
+        var xinc: Double
+
+        val vars1 = deepSecular(xmdf, omgadf, xnode, xn, tsince)
+        xmdf = vars1.xmdf
+        omgadf = vars1.omgadf
+        xnode = vars1.xnode
+        xn = vars1.xn
+        em = vars1.em
+        xinc = vars1.xinc
+
+        val a = (Globals.xke / xn).pow(2.0 / 3.0) * sqr(tempa)
+        var e = em - tempe
+        var xmam = xmdf + orbit.meanMotion * templ
+
+        val vars2 = deepPeriodics(e, xinc, omgadf, xnode, xmam, tsince)
+        e = vars2.e
+        xinc = vars2.xinc
+        omgadf = vars2.omgadf
+        xnode = vars2.xnode
+        xmam = vars2.xmam
+
+        val xl = xmam + omgadf + xnode
+        xn = Globals.xke / a.pow(1.5)
+
+        return finalPosition(xinc, omgadf, e, a, xl, xnode, xn, tsince)
+    }
+
+    override fun clone(o: Orbit): NoradBase = NoradSDP4(o)
+}

libs/ballistics-lib/src/main/kotlin/ballistics/tle_tools/zeptomoby/orbit/norad/NoradSGP4.kt

@@ -0,0 +1,89 @@
+package org.nstart.dep265.tletools.zeptomoby.orbit.norad
+
+import org.nstart.dep265.tletools.zeptomoby.core.Globals
+import org.nstart.dep265.tletools.zeptomoby.core.Globals.sqr
+import org.nstart.dep265.tletools.zeptomoby.core.eci.EciTime
+import org.nstart.dep265.tletools.zeptomoby.orbit.Orbit
+import kotlin.math.cos
+import kotlin.math.pow
+import kotlin.math.sin
+
+open class NoradSGP4(o: Orbit): NoradBase(o) {
+    protected val c5: Double
+    protected val omgcof: Double
+    protected val xmcof: Double
+    protected val delmo: Double
+    protected val sinmo: Double
+
+    init {
+        val etasq: Double = eta * eta
+
+        c5 = 2.0 * coef1 * orbit.semiMajor * betao2 * (1.0 + 2.75 * (etasq + eeta) + eeta * etasq)
+        omgcof = orbit.bStar * c3 * cos(orbit.argPerigee)
+        xmcof = -(2.0 / 3.0) * coef * orbit.bStar * Globals.ae / eeta
+        delmo = (1.0 + eta * cos(orbit.meanAnomaly)).pow(3.0)
+        sinmo = sin(orbit.meanAnomaly)
+    }
+
+    override fun position(tsince: Double): EciTime {
+        var isImp = false
+
+        if ((orbit.semiMajor * (1.0 - orbit.eccentricity) / Globals.ae) <
+            (220.0 / Globals.wgs.xkmPer + Globals.ae)) isImp = true
+
+        var d2 = 0.0
+        var d3 = 0.0
+        var d4 = 0.0
+
+        var t3cof = 0.0
+        var t4cof = 0.0
+        var t5cof = 0.0
+
+        if (!isImp) {
+            val c1sq: Double = c1 * c1
+            d2 = 4.0 * orbit.semiMajor * tsi * c1sq
+
+            val temp: Double = d2 * tsi * c1 / 3.0
+            d3 = (17.0 * orbit.semiMajor + s4) * temp
+            d4 = 0.5 * temp * orbit.semiMajor * tsi *
+                    (221.0 * orbit.semiMajor + 31.0 * s4) * c1
+            t3cof = d2 + 2.0 * c1sq
+            t4cof = 0.25 * (3.0 * d3 + c1 * (12.0 * d2 + 10.0 * c1sq))
+            t5cof = 0.2 * (3.0 * d4 + 12.0 * c1 * d3 + 6.0 *
+                    d2 * d2 + 15.0 * c1sq * (2.0 * d2 + c1sq))
+        }
+
+        val xmdf: Double = orbit.meanAnomaly + xmdot * tsince
+        val omgadf: Double = orbit.argPerigee + omgdot * tsince
+        val xnoddf: Double = orbit.raan + xnodot * tsince
+        var omega = omgadf
+        var xmp = xmdf
+        val tsq = tsince * tsince
+        val xnode: Double = xnoddf + xnodcf * tsq
+        var tempa: Double = 1.0 - c1 * tsince
+        var tempe: Double = orbit.bStar * c4 * tsince
+        var templ: Double = t2cof * tsq
+
+        if (!isImp) {
+            val delomg: Double = omgcof * tsince
+            val delm: Double = xmcof * ((1.0 + eta * cos(xmdf)).pow(3.0) - delmo)
+            val temp = delomg + delm
+            xmp = xmdf + temp
+            omega = omgadf - temp
+            val tcube = tsq * tsince
+            val tfour = tsince * tcube
+            tempa = tempa - d2 * tsq - d3 * tcube - d4 * tfour
+            tempe += orbit.bStar * c5 * (sin(xmp) - sinmo)
+            templ += t3cof * tcube + tfour * (t4cof + tsince * t5cof)
+        }
+
+        val a: Double = orbit.semiMajor * sqr(tempa)
+        val e: Double = orbit.eccentricity - tempe
+        val xl: Double = xmp + omega + xnode + orbit.meanMotion * templ
+        val xn: Double = Globals.xke / a.pow(1.5)
+
+        return finalPosition(orbit.inclination, omgadf, e, a, xl, xnode, xn, tsince)
+    }
+
+    override fun clone(o: Orbit): NoradBase = NoradSGP4(o)
+}

libs/ballistics-lib/src/main/kotlin/ballistics/types/Errors.kt

@@ -0,0 +1,44 @@
+package ballistics.types
+
+/**
+ * Код завершения расчета
+ */
+enum class BallisticsError(val value: Int) {
+    OK(0),
+    UNKNOWN_ERR(1),
+
+    EMPTY_NU(3),
+    H_MINIMUM_ERROR(4),
+    EMPTY_ORBITAL_POINTS(5),
+    EMPTY_FLIGHTLINE_POINTS(6),
+    ORB_POINTS_INTERVAL_ERROR(7),
+    STEPPER_ERROR(8),
+
+    EMPTY_PPI_LIST(9),
+    ;
+
+    companion object {
+        val map = BallisticsError.values().associate { it.value to it }
+
+        fun fromInt(
+            value: Int,
+            defValue: BallisticsError = BallisticsError.OK,
+        ): BallisticsError {
+            return map[value] ?: defValue
+        }
+    }
+}
+
+fun printError(err: BallisticsError): String {
+    return when (err) {
+        BallisticsError.OK -> "Успешно"
+        BallisticsError.UNKNOWN_ERR -> "нераспознанная ошибка"
+        BallisticsError.H_MINIMUM_ERROR -> "высота достигла минимального значения"
+        BallisticsError.EMPTY_NU -> "пустой массив начальных условий"
+        BallisticsError.EMPTY_ORBITAL_POINTS -> "не рассчитаны точки орбиты"
+        BallisticsError.ORB_POINTS_INTERVAL_ERROR -> "интервал расчета не пренадлежит интервалу расчета точек орбиты"
+        BallisticsError.EMPTY_PPI_LIST -> "пустой массив ППИ"
+        BallisticsError.STEPPER_ERROR -> "ошибка при выходе на заданное время"
+        else -> "неизвестный код ошибки"
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/types/Types.kt

@@ -0,0 +1,341 @@
+package ballistics.types
+
+import ballistics.utils.math.Vector3D
+import java.time.LocalDateTime
+import kotlin.math.PI
+
+/**
+ * Тип метода интегрирования
+ */
+enum class IntegrationType {
+    RUNG4,
+    ADAMS7,
+}
+
+/**
+ * Тип модели движения
+ */
+enum class ModDVType(val value: Int) {
+    FOTO(1),
+
+    KONDOR(6),
+
+    METEORM1(14),
+    METEORM2(15),
+
+    BARS(16),
+
+    KONDOR_PROGNOZ(19),
+    ;
+
+    companion object {
+        val map = values().associate { it.value to it }
+
+        fun fromInt(
+            value: Int,
+            defValue: ModDVType = FOTO,
+        ): ModDVType {
+            return map[value] ?: defValue
+        }
+    }
+}
+
+/**
+ * Тип модели Земли
+ */
+enum class EarthType {
+    PZ90d0,
+    PZ90d02,
+}
+
+enum class AstroType {
+    ATJ2000,
+}
+
+/**
+ * Тип базовой СК
+ */
+enum class WorkCSType {
+    /**
+     * орбитальная
+     */
+    WCSOrbit,
+
+    /**
+     * путевая / трассовая
+     */
+    WCSPath,
+}
+
+/**
+ * Тип учета угла тангажа в расчете матрицы планирования
+ */
+enum class TangageType {
+    /**
+     * Конструктивный
+     */
+    TTConstructive,
+
+    /**
+     * Упреждающий
+     */
+    TTProactive,
+}
+
+
+/**
+ * Параметры точки орбиты
+ */
+data class OrbitalPoint(val t: Double, var vit: Int, val r: Vector3D, val v: Vector3D) {
+    constructor() : this(0.0, 0, Vector3D(), Vector3D())
+}
+
+/**
+ * Параметры витка (ВУЗ)
+ */
+data class RevolutionParameter(var vuz: OrbitalPoint, var lVuz: Double, var hVuz: Double) {
+    constructor() : this(OrbitalPoint(), 0.0, 0.0)
+}
+
+/**
+ * Начальные условия движения
+ */
+class InitialConditions(var point: OrbitalPoint, var sBall: Double, var f81: Double) {
+    constructor() : this(OrbitalPoint(), 0.0, 100.0)
+}
+
+class TLE(val tle1: String, val tle2: String) {
+    constructor() : this("", "")
+}
+
+/**
+ * затенение по азимуму
+ * азимут в диапазоне 0..2*PI
+ * при переходе через 0 долготы нужно заводить 2 записи! до 0 и после 0!
+ */
+data class ShadowAU(
+    val azimuthStart: Double,
+    val azimuthEnd: Double,
+    val elevation: Double,
+)
+
+/**
+ * Параметры ППИ
+ */
+data class PPI(
+    val ppiNum: Int = 0,
+    val auNum: Int = 0,
+    val lat: Double = 0.0,
+    val long: Double = 0.0,
+    val height: Double = 0.0,
+    val elevMin: Double = 0.0,
+    val elevMax: Double = 0.0,
+    var shadowMin: List<ShadowAU>? = null,
+    var shadowMax: List<ShadowAU>? = null,
+) {
+    constructor(
+        ppiNum: Int,
+        auNum: Int,
+        lat: Double,
+        long: Double,
+        height: Double,
+        elevMin: Double,
+        elevMax: Double,
+    ) : this (ppiNum, auNum, lat, long, height, elevMin, elevMax, null, null)
+}
+
+/**
+ * Параметры ППИ(расширенные)
+ */
+internal class PPIExtParams() {
+    var nip = 0
+    var au = 0
+    var bn = 0.0
+    var xn = 0.0
+    var yn = 0.0
+    var zn = 0.0
+    var sl = 0.0
+    var cl = 0.0
+    var sb = 0.0
+    var cb = 0.0
+}
+
+/**
+ * Параметры ППИ(расширенные)
+ */
+internal class PPIExt() {
+    var params = PPIExtParams()
+    var elevMin = 0.0
+    var elevMax = 0.0
+    var zoneIn = VisibilityParametersZRV()
+    var zoneMax = VisibilityParametersZRV()
+    var zoneOut = VisibilityParametersZRV()
+    var vit = 0
+    var kaY = 0.0
+    var t1 = 0.0
+    var gam1 = 0.0
+    var a1 = 0.0
+    var vy1 = 0.0
+    var t2 = 0.0
+    var gam2 = 0.0
+    var a2 = 0.0
+    var vy2 = 0.0
+    var shadowMin: Iterable<ShadowAU>? = null
+    var shadowMax: Iterable<ShadowAU>? = null
+
+    fun isVisible(
+        azimuth: Double,
+        elevation: Double,
+    ): Boolean {
+        if (elevation > elevMax) {
+            return false
+        }
+        if (shadowMin != null) {
+            val a = if (azimuth >= 0) azimuth else azimuth + 2 * PI
+            val shad = shadowMin!!.find { it.azimuthStart <= a && a <= it.azimuthEnd }
+            if (shad != null && shad.elevation > elevation) {
+                return false
+            }
+        }
+        if (shadowMax != null) {
+            val a = if (azimuth >= 0) azimuth else azimuth + 2 * PI
+            val shad = shadowMax!!.find { it.azimuthStart <= a && a <= it.azimuthEnd }
+            if (shad != null && shad.elevation < elevation) {
+                return false
+            }
+        }
+        return true
+    }
+}
+
+/**
+ * Параметры видимости КА-ППИ
+ */
+data class VisibilityParametersZRV(var t: Double, var range: Double, var azimuth: Double, var elevation: Double) {
+    constructor() : this(0.0, 0.0, 0.0, 0.0)
+}
+
+/**
+ * Параметры зоны радиовидимости КА-ППИ
+ */
+data class ZRV(
+    val ppi: Int,
+    val au: Int,
+    val vit: Int,
+    val zoneIn: VisibilityParametersZRV,
+    val zoneMax: VisibilityParametersZRV,
+    val zoneOut: VisibilityParametersZRV,
+)
+
+class Orientation(
+    var tang: Double,
+    var kren: Double,
+    var risk: Double,
+)
+
+/**
+ * Параметры точки на Земной поверхности
+ */
+class BLHPoint(
+    var lat: Double,
+    var long: Double,
+    var h: Double,
+)
+
+/**
+ * Параметры точки на границе полосы обзора / трассе полета
+ */
+class THBLPoint(val lat: Double, val long: Double, val range: Double, val sunAngle: Double)
+
+/**
+ * Параметры трассы полета и полосы обзора
+ */
+class FlightLine(
+    val t: Double,
+    val vit: Int,
+    val pv: Int,
+    val rightOuterSwath: THBLPoint,
+    val rightInternalSwath: THBLPoint,
+    val flightLine: THBLPoint,
+    val leftInternalSwath: THBLPoint,
+    val leftOuterSwath: THBLPoint,
+)
+
+class FleghtLineSector(
+    val tStart: Double,
+    val tStop: Double,
+)
+
+/**
+ * Объект оперативного каталога
+ */
+class OPKatObj(
+    val objON: Int,
+    val objN: Int,
+    val objUUID : String,
+    val pointNumb: Int,
+    val lat: Double,
+    val long: Double,
+    val height: Double,
+    val rollMin: Double,
+    val rollMax: Double,
+    val sunAngleMin: Double,
+)
+
+/**
+ * Параметры видимости точки наблюдения
+ */
+class PointViewParams(
+    val objON: Int,
+    val objN: Int,
+    val objUUID : String,
+    val pointNumb: Int,
+    val vit: Int,
+    val traverz: Double,
+    val latTraverz: Double,
+    val longTraverz: Double,
+    val orientation: Orientation,
+    val range: Double,
+    val sunAngle: Double,
+    val sightAngle: Double,
+    val pv : Int,
+)
+
+data class KeplerParams(
+    var dmv: Double = 0.0,
+    var ael: Double = 0.0,
+    var e: Double = 0.0,
+    var nakl: Double = 0.0,
+    var omegab: Double = 0.0,
+    var omegam: Double = 0.0,
+    var u: Double = 0.0,
+    var v: Double = 0.0,
+    var o: Double = 0.0,
+    var eA: Double = 0.0,
+    var tau: Double = 0.0,
+    var t: Double = 0.0,
+    var rP: Double = 0.0,
+    var rA: Double = 0.0,
+)
+
+
+class TLEParams(
+    var satName : String = "",
+    var noradId : Long = 0,
+    var revolution : Long = 0,
+    var time : LocalDateTime = LocalDateTime.now(),
+    var inclination : Double = 0.0,
+    var perigee : Double = 0.0,
+    var apogee: Double = 0.0,
+    var argPerigee: Double = 0.0,
+    var eccentricity: Double = 0.0,
+    var major: Double = 0.0,
+    var minor: Double = 0.0,
+    var meanAnomaly: Double = 0.0,
+    var period: Double = 0.0,
+    var semiMajor: Double = 0.0,
+    var semiMinor: Double = 0.0,
+    var ascendingNode: Double = 0.0,
+    var meanMotion: Double = 0.0,
+    var meanMotionTLE: Double = 0.0
+)

libs/ballistics-lib/src/main/kotlin/ballistics/utils/DateTimeUtils.kt

@@ -0,0 +1,8 @@
+package ballistics.utils
+
+import java.time.LocalDateTime
+import java.time.ZoneOffset
+
+fun toDateTime(t: Double) = LocalDateTime.ofEpochSecond(t.toLong(), (t % 1 * 1e9).toInt(), ZoneOffset.UTC)
+
+fun fromDateTime(t: LocalDateTime) = t.toEpochSecond(ZoneOffset.UTC) + t.nano / 1e9

libs/ballistics-lib/src/main/kotlin/ballistics/utils/astro/AstronomerJ2000.kt

@@ -0,0 +1,406 @@
+package ballistics.utils.astro
+
+import ballistics.types.EarthType
+import ballistics.types.OrbitalPoint
+import ballistics.utils.earth.getEarth
+import ballistics.utils.math.Matrix3D
+import ballistics.utils.math.Vector3D
+import java.time.LocalDateTime
+import java.time.ZoneOffset
+import kotlin.math.PI
+import kotlin.math.atan
+import kotlin.math.atan2
+import kotlin.math.cos
+import kotlin.math.sin
+import kotlin.math.sqrt
+import kotlin.math.truncate
+
+class AstronomerJ2000(var earthType: EarthType) {
+    var ase = 149600000.0
+    var earth = getEarth(earthType)
+    var baseDate =
+        LocalDateTime.of(2000, 1, 1, 12, 0, 0, 0).toEpochSecond(ZoneOffset.UTC)
+            .toDouble() / 86400
+
+    fun si2000(tDtTm: Double): Double {
+        var d = tDtTm / 86400.0
+        var dat = truncate(d)
+        var tt = d - dat
+        if (tt < 0.125) {
+            dat--
+        }
+
+        var t = (dat - baseDate) / (36525.0)
+        var t2 = t * t
+        var t3 = t2 * t
+
+        var lt =
+            2.355555743494e0 + 8328.691425719e0 * t +
+                1.545547e-4 * t2 + 2.50e-7 * t3
+        var llt =
+            6.240060126913e0 + 628.301955171e0 * t - 2.681989e-6 * t2 -
+                6.593e-10 * t3
+        var ft =
+            1.6279050815375e0 + 8433.466156916e0 * t - 6.181956e-5 * t2 +
+                5.0275e-9 * t3
+        var dt =
+            5.1984665886602e0 + 7771.377145594e0 * t - 5.559397e-4 * t2 +
+                3.196e-8 * t3
+        var omt =
+            2.1824391966156e0 - 33.757044613e0 * t + 3.622625e-5 * t2 +
+                3.734e-8 * t3
+
+        var psi =
+            (-17.1996e0 - 0.01742 * t) * sin(omt) +
+                (0.2062e0 + 0.00002 * t) * sin(2 * omt) +
+                (-1.3187e0 - 0.00016 * t) * sin(2 * (ft - dt + omt)) +
+                (0.1426e0 - 0.00034 * t) * sin(llt) +
+                (-0.0517 + 0.00012 * t) * sin(llt + 2 * ft - 2 * dt + 2 * omt) +
+                (0.0217 - 0.00005 * t) * sin(-llt + 2 * ft - 2 * dt + 2 * omt) +
+                (0.0129 + 0.00001 * t) * sin(2 * ft - 2 * dt + omt)
+
+        var r =
+            (-0.2274e0 - 0.00002 * t) * sin(2 * ft + 2 * omt) +
+                (0.0712e0 + 0.00001 * t) * sin(lt) +
+                (-0.0386e0 - 0.00004 * t) * sin(2 * ft + omt) -
+                0.0301e0 * sin(lt + 2 * ft + 2 * omt) -
+                0.0158 * sin(lt - 2 * dt) +
+                0.0123 * sin(-lt + 2 * ft + 2 * omt)
+
+        psi += r
+
+        if (tt < 0.125) {
+            tt += 0.875
+        } else {
+            tt -= 0.125
+        }
+
+        r = 84381.448e0 - 46.84024e0 * t - 0.00059e0 * t2 + 0.001813e0 * t3
+        r = r * PI / 180.0 / 3600.0
+        psi = psi * PI / 180.0 / 3600.0
+        r = cos(r)
+        r = psi * r
+        var s =
+            (
+                24110.54841 + 8640184.812866 * t + 0.093104 * t2 -
+                    6.2e-6 * t3
+            ) * 2 * PI / 86400.0e0
+
+        s = s + r + 2 * PI * tt * 1.002737811906E0
+
+        var ds = s / (PI + PI)
+        var ts = truncate(ds)
+
+        s = (ds - ts) * (PI + PI)
+
+        return s
+    }
+
+    fun grinvToASK(gsk: OrbitalPoint): OrbitalPoint {
+        var m = Matrix3D()
+        m.makeOzMatrix(si2000(gsk.t))
+
+        var ra = m * gsk.r
+        var va = m * gsk.v
+
+        va.x -= earth.wEarth * ra.y
+        va.y += earth.wEarth * ra.x
+
+        return OrbitalPoint(gsk.t, gsk.vit, ra, va)
+    }
+
+    fun grinvToASK(
+        gsk: OrbitalPoint,
+        si: Double,
+    ): OrbitalPoint {
+        var m = Matrix3D()
+        m.makeOzMatrix(si)
+        var ra = m * gsk.r
+        var va = m * gsk.v
+
+        va.x -= earth.wEarth * ra.y
+        va.y += earth.wEarth * ra.x
+
+        return OrbitalPoint(gsk.t, gsk.vit, ra, va)
+    }
+
+    fun grinvToASK(
+        gsk: Vector3D,
+        t: Double,
+    ): Vector3D {
+        var m = Matrix3D()
+        m.makeOzMatrix(si2000(t))
+        return m * gsk
+    }
+
+    fun askToGrinvich(
+        ask: Vector3D,
+        t: Double,
+    ): Vector3D {
+        var m = Matrix3D()
+        m.makeOzMatrix(-si2000(t))
+        return m * ask
+    }
+
+    fun askToGrinvich(ask: OrbitalPoint): OrbitalPoint {
+        var m = Matrix3D()
+        m.makeOzMatrix(-si2000(ask.t))
+
+        var ra = m * ask.r
+        var va = m * ask.v
+
+        va.x += earth.wEarth * ra.y
+        va.y -= earth.wEarth * ra.x
+
+        return OrbitalPoint(ask.t, ask.vit, ra, va)
+    }
+
+    fun askToGrinvich(
+        ask: OrbitalPoint,
+        si: Double,
+    ): OrbitalPoint {
+        var m = Matrix3D()
+        m.makeOzMatrix(-si)
+        var ra = m * ask.r
+        var va = m * ask.v
+
+        va.x += earth.wEarth * ra.y
+        va.y -= earth.wEarth * ra.x
+
+        return OrbitalPoint(ask.t, ask.vit, ra, va)
+    }
+
+    fun moonCoordinates(tm: Double): Vector3D {
+        var d: Double = tm / 86400.0
+        d = d - 0.125
+        val t: Double = (d - baseDate) / 36525.0 // Юлианский период ot 12h 1.1.2000
+        val t2 = t * t
+        val t3 = t2 * t
+        // (* Средняя аномалия ЛУНЫ*)
+        val lL = 2.355555743494e0 + 8328.691425719e0 * t + 1.545547e-4 * t2 + 2.50e-7 * t3
+        // (* Средняя аномалия СОЛНЦА*)
+        val lS = 6.240060126913e0 + 628.301955171e0 * t - 2.681989e-6 * t2 - 6.593e-10 * t3
+        // (* Средний аргумент широты ЛУНЫ *)
+        val f = 1.6279050815375e0 + 8433.466156916e0 * t - 6.181956e-5 * t2 + 5.0275e-9 * t3
+        // (* Разность средних долгот ЛУНЫ и СОЛНЦА *)
+        d = 5.1984665886602e0 + 7771.377145594e0 * t - 5.559397e-4 * t2 + 3.196e-8 * t3
+        val d2 = d + d
+        // (*  Вычисление горизонтального паралакса и радиуса ЛУНЫ *)
+        val pl1 =
+            4.8481368110953599e-6 * (
+                3422.7 + 28.233 * cos(d2) + 3.086 * cos(lL + d2) + 186.539 * cos(lL) + 34.311 *
+                    cos(
+                        lL - d2,
+                    ) + 0.6 * cos(lL - 4 * d) - 0.399 * cos(lS) +
+                    1.917 * cos(lS - d2) - 0.978 * cos(d) + 10.165 * cos(2 * lL) -
+                    0.949 * cos(lL + lS) + 1.443 * cos(lL + lS - d2) + 1.152 * cos(lL - lS) + 0.621 * cos(3 * lL) - 0.713 *
+                    cos(
+                        lL - 2 * f,
+                    )
+            )
+        var rl = 6378.14 / pl1
+        rl = rl * 1.0e3
+        // (* Средняя долгота ЛУНЫ *)
+        val vL =
+            3.8103405583329555 + 8399.7091123847003500 * t -
+                0.0000281288897780 * t2 + 0.0000000921145994 * t3
+        // (* Вычисление эклиптической долготы ЛУНЫ*)
+        val vSin =
+            4.8481368110953599e-6 * (
+                2369.9 * sin(d2) + 191.95 * sin(lL + d2) + 22639.5 * sin(lL) - 4586.42 * sin(lL - d2) - 668.11 *
+                    sin(
+                        lS,
+                    ) - 165.14 * sin(lS - d2) - 125.15 * sin(d) + 769.01 * sin(2 * lL) - 211.65 * sin(2 * lL - d2) - 109.66 *
+                    sin(
+                        lL + lS,
+                    ) - 205.96 * sin(lL + lS - d2) + 147.69 * sin(lL - lS) - 411.6 * sin(2 * f)
+            )
+        val arg = vL + vSin
+        val cl = cos(arg)
+        val sl = sin(arg)
+        // (* Вычисление эклиптической широты Луны *)
+        val betta =
+            4.8481368110953599e-6 * (
+                117.26 * sin(f + d2) + 18461.48 * sin(f) -
+                    623.65 * sin(f - d2) + 1010.18 * sin(lL + f) -
+                    166.57 * sin(lL + f - d2) + 199.48 * sin(-lL + f + d2) -
+                    999.69 * sin(-lL + f)
+            )
+        val cb = cos(betta)
+        val sb = sin(betta)
+        // (* Средний наклон эклиптики  *)
+        val ep = 0.4090928042223290 - 0.0002270878917845 * t - 0.0000000028604007 * t2 + 0.0000000087896720 * t3
+        val ce = cos(ep)
+        val se = sin(ep)
+        var moon = Vector3D()
+        moon.x = rl * cl * cb
+        moon.y = rl * (sl * cb * ce - sb * se)
+        moon.z = rl * (sb * ce + sl * cb * se)
+
+        return moon
+    }
+
+    fun sunCoordinates(tm: Double): Vector3D {
+        var d = tm / 86400.0
+        var dd = d - 0.125
+        var t = (dd - baseDate) / 36525.0
+        var t2 = t * t
+        var t3 = t2 * t
+
+        //  Средняя долгота
+        var l =
+            4.8481368110953599e-6 * (
+                1009677.850 +
+                    (100 * 1296000.0 + 2771.270) * t + 1.089 * t2
+            )
+
+        //  Долгота перигелия
+        var p =
+            4.8481368110953599e-6 * (
+                1018578.046 + 6190.046 * t +
+                    1.666 * t2 + 0.012 * t3
+            )
+
+        //  Эксцентриситет земной орбиты
+        var ek = 0.0167086342 - 0.00004203654 * t - 0.00000012673 * t2
+
+        // Средний наклон эклиптики
+        var e1 =
+            4.8481368110953599e-6 * (
+                84381.448 - 46.84024 * t -
+                    0.00059 * t2 + 0.001813 * t3
+            )
+
+        // Истинный наклон эклиптики
+        var oMT =
+            2.1824391966156e0 - 33.757044613e0 * t + 3.622625e-5 * t2 +
+                3.734e-8 * t3
+
+        var psi = 4.8481368110953599e-6 * 9.2025e0 * cos(oMT)
+
+        e1 += psi
+
+        //  Средняя аномалия Земли
+        var g = l - p
+        var e = g + ek
+
+        for (i in 1..5)
+            e = g + ek * sin(e)
+
+        // Радиус Земли *)
+        var r = ase * (1.0 - ek * cos(e)) * 1.0e3
+
+        // Истинная аномалия
+        var sinv = ase / r * sqrt(1.0 - ek * ek) * sin(e)
+        var cosv = ase / r * (cos(e) - ek)
+        var v = atan2(sinv, cosv)
+
+        // Эклиптическая истинная долгота *)
+        var la = v + p
+
+        var sun = Vector3D()
+        // Собственно вычисления координат *)
+        sun.x = r * cos(la)
+        var y = r * sin(la)
+        sun.y = y * cos(e1)
+        sun.z = y * sin(e1)
+
+        return sun
+    }
+
+    fun sunCoordinatesGSK(t: Double): Vector3D {
+        var m = Matrix3D()
+        m.makeOzMatrix(-si2000(t))
+        var sun = sunCoordinates(t)
+        return m * sun
+    }
+
+    fun sunAngle(
+        t: Double,
+        currentObj: Vector3D,
+    ): Double {
+        var sun = sunCoordinatesGSK(t)
+
+        var blh = earth.xyz2blh(currentObj)
+
+        var m =
+            Matrix3D(
+                Vector3D(cos(blh.lat) * cos(blh.long), cos(blh.lat) * sin(blh.long), sin(blh.lat)),
+                Vector3D(-sin(blh.long), cos(blh.long), 0.0),
+                Vector3D(-sin(blh.lat) * cos(blh.long), -sin(blh.lat) * sin(blh.long), cos(blh.lat)),
+            )
+        var pp = m * currentObj
+        var sp = m * sun
+
+        var r = sp - pp
+
+        return atan(r.x / sqrt(r.z * r.z + r.y * r.y))
+    }
+
+    fun sunMoonCorrection(
+        t3: Double,
+        x: Double,
+        y: Double,
+        z: Double,
+    ): Vector3D {
+        var corrections = Vector3D()
+        val lm = 4902.7854
+        val mss = 132712518017.51
+        val kapa = 0.0
+
+        var sun = sunCoordinates(t3)
+        sun = sun * (1 / 1000.0)
+        var moon = moonCoordinates(t3)
+        moon = moon * (1 / 1000.0)
+
+        // Учет влияния Луны
+        var dpx: Double
+        var dpy: Double
+        var dpz: Double
+
+        dpx = moon.x - x
+        dpy = moon.y - y
+        dpz = moon.z - z
+
+        var r = dpx * dpx + dpy * dpy + dpz * dpz
+        var dr = sqrt(r)
+        var pp1 = r * dr
+
+        r = moon.x * moon.x + moon.y * moon.y + moon.z * moon.z
+        dr = sqrt(r)
+        var pp2 = r * dr
+
+        corrections.x = lm * (dpx / pp1 - moon.x / pp2)
+        corrections.y = lm * (dpy / pp1 - moon.y / pp2)
+        corrections.z = lm * (dpz / pp1 - moon.z / pp2)
+
+        dpx = sun.x - x
+        dpy = sun.y - y
+        dpz = sun.z - z
+
+        r = dpx * dpx + dpy * dpy + dpz * dpz
+        dr = sqrt(r)
+        pp1 = r * dr
+
+        dr = sun.module()
+        r = dr * dr
+        pp2 = r * dr
+
+        var kap = kapa
+        val rap = sqrt(x * x + y * y + z * z)
+        val cf = (sun.x * x + sun.y * y + sun.z * z) / (dr * rap)
+        if (cf <= 0) {
+            val rab1 = rap * sqrt(1 - cf * cf)
+            if (rab1 < 6378.388) kap = 0.0
+        }
+
+        val p1 = mss / pp1 * (1 - kap)
+        val p2 = mss / pp2
+        // Закоментарить для неучета Солнца
+        corrections.x += (dpx * p1 - sun.x * p2)
+        corrections.y += (dpy * p1 - sun.y * p2)
+        corrections.z += (dpz * p1 - sun.z * p2)
+
+        return corrections
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/utils/atmosphere/Atmosphere2004.kt

@@ -0,0 +1,237 @@
+package ballistics.utils.atmosphere
+
+import kotlin.math.cos
+import kotlin.math.exp
+import kotlin.math.floor
+import kotlin.math.sin
+import kotlin.math.sqrt
+import kotlin.math.truncate
+
+class Atmosphere2004() {
+    val aAtm =
+        listOf(
+            -2.53418e-02,
+            -2.44075e-03,
+            3.08389e-06,
+            2.90115e-06,
+            -4.99606e-08,
+            3.36327e-10,
+            -1.0966e-12,
+            1.73227e-15,
+            -1.06271e-18,
+        )
+
+    val fi1 =
+        listOf(
+            0.5411,
+            0.5515,
+            0.5585,
+            0.5585,
+            0.5585,
+            0.5585,
+            0.5585,
+        )
+
+    val nAtm =
+        listOf(
+            -1.10059e+00,
+            -1.37626e-02,
+            -4.05631e-05,
+            4.64077e-08,
+            1.31904e+01,
+            -1.00640e-01,
+            2.66200e-04,
+            -2.10328e-07,
+            -6.83628e-01,
+            8.17955e-03,
+            -5.37817e-06,
+            4.18766e+00,
+            -3.63633e-02,
+            1.28280e-04,
+            -1.18060e-07,
+            -3.77859e+00,
+            4.59058e-05,
+            -5.42317e-05,
+            4.67279e-08,
+            8.62648e+00,
+            -6.45632e-02,
+            1.67737e-04,
+            -1.23719e-07,
+            -5.32838e-01,
+            7.17357e-03,
+            -4.13071e-06,
+            3.27388e+00,
+            -2.61311e-02,
+            9.26510e-05,
+            -8.03080e-08,
+            -2.06474e+00,
+            -1.16758e-02,
+            -2.26591e-05,
+            2.19816e-08,
+            5.84701e+00,
+            -4.29259e-02,
+            1.09960e-04,
+            -7.53447e-08,
+            -2.04545e-01,
+            5.44465e-03,
+            -2.53356e-06,
+            8.09259e-01,
+            -7.46853e-03,
+            4.52020e-05,
+            -4.13791e-08,
+            2.27394e+00,
+            -3.94166e-02,
+            3.95835e-05,
+            -2.19992e-08,
+            4.16132e+00,
+            -2.99269e-02,
+            7.56965e-05,
+            -4.77865e-08,
+            2.42053e-02,
+            4.38092e-03,
+            -1.76272e-06,
+            1.24462e+00,
+            -7.76442e-03,
+            3.88966e-05,
+            -3.23585e-08,
+            1.67658e+00,
+            -3.55099e-02,
+            3.46586e-05,
+            -1.95730e-08,
+            3.14269e+00,
+            -2.20308e-02,
+            5.48269e-05,
+            -3.14805e-08,
+            1.71702e-01,
+            3.77273e-03,
+            -1.45257e-06,
+            1.34930e-01,
+            4.11962e-04,
+            1.88965e-05,
+            -1.74602e-08,
+            -1.37615e+00,
+            -1.68609e-02,
+            3.40709e-07,
+            1.72566e-09,
+            2.43907e+00,
+            -1.65847e-02,
+            4.06547e-05,
+            -2.08670e-08,
+            2.59209e-01,
+            3.45977e-03,
+            -1.38895e-06,
+            1.47530e-01,
+            1.95702e-03,
+            1.17537e-05,
+            -1.08184e-08,
+            2.38462e+00,
+            -3.87222e-02,
+            4.60938e-05,
+            -2.80992e-08,
+            1.76084e+00,
+            -1.10685e-02,
+            2.67151e-05,
+            -1.08944e-08,
+            2.60394e-01,
+            3.53130e-03,
+            -1.54729e-06,
+            -7.53476e-03,
+            3.97475e-03,
+            5.08266e-06,
+            -5.66222e-09,
+        )
+
+    private fun floorm(a: Double): Double {
+        return (if (a >= 0) floor(a) else -1 * (floor(a * -1)))
+    }
+
+    fun atm2004Kav(
+        f81: Double,
+        deny: Double,
+        h: Double,
+        fld3: Array<Double>,
+        lxg: Double,
+        lyg: Double,
+        lzg: Double,
+    ): Double {
+        var roh : Double
+
+        val ro0 = 1.58868e-08
+        val n0A = 3.03853
+        val n1A = 1.875e-03
+
+        val rg: Double = sqrt(fld3[3] * fld3[3] + fld3[4] * fld3[4] + fld3[5] * fld3[5])
+        val xort: Double = fld3[3] / rg
+        val yort: Double = fld3[4] / rg
+        val zort: Double = fld3[5] / rg
+
+        var j : Int
+        val f1: Double = floorm(f81 / 25.0 + 0.5)
+        var f0 = 25.0 * f1
+        if (f0 >= 225.0) {
+            f0 = 250.0
+            j = 7
+        } else {
+            j = truncate(f1 - 2.0).toInt()
+        }
+        if (j < 1) j = 1
+
+        val fUnderscore = (f81 - f0) / f0
+        val d = deny - 1.0
+        val aD: Double =
+            aAtm[0] + d * (
+                aAtm[1] + d * (
+                    aAtm[2] + d * (
+                        aAtm[3] + d * (
+                            aAtm[4] + d * (
+                                aAtm[5] + d * (aAtm[6] + d * (aAtm[7] + d * aAtm[8]))
+                            )
+                        )
+                    )
+                )
+            )
+        val h1 = h * 0.001
+        val n = n0A + h1 * n1A
+        val x0ro: Double = lxg * cos(fi1[j - 1]) - lyg * sin(fi1[j - 1])
+        val y0ro: Double = lxg * sin(fi1[j - 1]) + lyg * cos(fi1[j - 1])
+        val cosfi = xort * x0ro + yort * y0ro + zort * lzg
+        val k10 = 0.5 * (1 + cosfi)
+
+        val mk: Double = Math.pow(k10, 0.5 * n)
+        val nI = (15 * (j - 1)).toDouble()
+
+        val aN: Double =
+            nAtm[((nI + 0).toInt())] + h1 * (
+                nAtm[((nI + 1).toInt())] + h1 * (
+                    nAtm[((nI + 2).toInt()) ]+ h1 *
+                        nAtm[(
+                            (nI + 3).toInt()
+                        )]
+                )
+            )
+        val cN: Double =
+            nAtm[((nI + 4).toInt())] + h1 * (
+                nAtm[(nI + 5).toInt()] + h1 * (
+                    nAtm[(nI + 6).toInt()] + h1 *
+                            nAtm[(nI + 7).toInt()]
+                )
+            )
+        val dN: Double =
+            nAtm[(nI + 8).toInt()] + h1 * (nAtm[(nI + 9).toInt()] + h1 * nAtm[(nI + 10).toInt()])
+        val lN: Double =
+            nAtm[(nI + 11).toInt()] + h1 * (
+                nAtm[(nI + 12).toInt()] + h1 * (
+                    nAtm[(nI + 13).toInt()] + h1 *
+                            nAtm[(nI + 14).toInt()]
+                )
+            )
+        val k0 = 1.0 + lN * fUnderscore
+        val k1 = cN * mk
+        val k2 = dN * aD
+        roh = ro0 * exp(aN)
+
+        roh *= k0 * (1.0 + k1 + k2)
+
+        return roh
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/utils/atmosphere/Atmosphere62.kt

@@ -0,0 +1,16 @@
+package ballistics.utils.atmosphere
+
+internal fun getAtm62(
+    h: Double,
+    indAtm: Int,
+): Double {
+    val ht = arrayOf<Double>(0.0, 0.2e+5, 0.6e+5, 1.0e+5, 1.5e+5, 3.0e+5, 6e+5, 9e+5, 10e6)
+    val a = arrayOf<Double>(0.125, 0.9091e-2, 0.263e-4, 0.4141e-7, 0.21727e-9, 0.48607e-11, 0.8904e-13, 0.6497e-14, 1e-17)
+    val k1 = arrayOf<Double>(-0.2639e-8, 0.4407e-9, -0.256e-8, 0.14688e-8, 0.8004e-10, 0.7111e-11, 0.1831e-11, 0.0, 0.0)
+    val k2 = arrayOf<Double>(0.7825e-4, 0.16375e-3, 0.5905e-4, 0.1787e-3, 0.37336e-4, 0.15467e-4, 0.9275e-5, 0.954e-5, 1e-6)
+
+    val rab = k1[indAtm - 1] * ((h - ht[indAtm - 1]) * (h - ht[indAtm - 1])) - k2[indAtm - 1] * (h - ht[indAtm - 1])
+    val ro = a[indAtm - 1] * Math.exp(rab)
+
+    return ro
+}

libs/ballistics-lib/src/main/kotlin/ballistics/utils/atmosphere/Atmosphere81.kt

@@ -0,0 +1,35 @@
+package ballistics.utils.atmosphere
+
+internal fun getAtm81(h: Double): Double {
+    val ht = arrayOf<Double>(0.0, 0.2e+5, 0.6e+5, 1.0e+5, 1.5e+5, 3.0e+5, 6.0e+5, 9.0e+5)
+    val a = arrayOf<Double>(0.12522, 0.91907e-2, 0.31655e-4, 0.54733e-7, 0.20474e-9, 0.19019e-11, 0.11495e-13, 0.58038e-15)
+    val k1 = arrayOf<Double>(-0.20452e-8, 0.62669e-9, -0.86999e-9, 0.12870e-8, 0.10167e-9, 0.97266e-11, 0.15127e-10, 0.0)
+    val k2 = arrayOf<Double>(0.90764e-4, 0.16739e-3, 0.12378e-3, 0.17527e-3, 0.45825e-4, 0.19885e-4, 0.14474e-4, 0.39247e-5)
+
+    var i: Int
+    var rab : Double
+    var ro  : Double
+    if (h < 9.0e+5) {
+        i = 0
+        i += 1
+        rab = ht[0] - h
+        while (rab < 0.0) {
+            i += 1
+            rab = ht[i - 1] - h
+        }
+        if (rab == 0.0) {
+            ro = a[i - 1]
+        } else {
+            i -= 1
+            rab = k1[i - 1] * (h - ht[i - 1]) * (h - ht[i - 1]) - k2[i - 1] * (h - ht[i - 1])
+            ro = a[i - 1] * Math.exp(rab)
+        }
+    } else {
+        i = 9
+        i -= 1
+        rab = k1[i - 1] * (h - ht[i - 1]) * (h - ht[i - 1]) - k2[i - 1] * (h - ht[i - 1])
+        ro = a[i - 1] * Math.exp(rab)
+    }
+
+    return ro
+}

libs/ballistics-lib/src/main/kotlin/ballistics/utils/earth/AbstractEarth.kt

@@ -0,0 +1,27 @@
+package ballistics.utils.earth
+
+import ballistics.types.BLHPoint
+import ballistics.utils.math.Vector3D
+
+abstract class AbstractEarth {
+    var fM = 398600.44e9
+    var ekvRadius = 6378136.0
+    var wEarth = 0.7292115e-4
+    var alphaEllips = 1 / 298.257839303
+    var middleRadius = 6378100.0
+    var polarRadius = (1 - alphaEllips) * ekvRadius
+    var c20 = -1.08262741e-3
+    var eElips2 = 2 * alphaEllips - alphaEllips * alphaEllips
+
+    fun blh2xyz(blh: BLHPoint): Vector3D {
+        return blh2xyz(blh.lat, blh.long, blh.h)
+    }
+
+    abstract fun blh2xyz(
+        b: Double,
+        l: Double,
+        h: Double,
+    ): Vector3D
+
+    abstract fun xyz2blh(r: Vector3D): BLHPoint
+}

libs/ballistics-lib/src/main/kotlin/ballistics/utils/earth/EarthFactory.kt

@@ -0,0 +1,10 @@
+package ballistics.utils.earth
+
+import ballistics.types.EarthType
+
+fun getEarth(earth: EarthType): AbstractEarth {
+    return when (earth) {
+        EarthType.PZ90d0 -> EarthPZ90d0()
+        EarthType.PZ90d02 -> EarthPZ90d02()
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/utils/earth/EarthPZ90d0.kt

@@ -0,0 +1,57 @@
+package ballistics.utils.earth
+
+import ballistics.types.BLHPoint
+import ballistics.utils.math.Vector3D
+import kotlin.math.PI
+import kotlin.math.atan2
+import kotlin.math.hypot
+import kotlin.math.sqrt
+
+open class EarthPZ90d0 : AbstractEarth() {
+    init {
+        fM = 398600.44e9
+        ekvRadius = 6378136.0
+        wEarth = 0.7292115e-4
+        alphaEllips = 1 / 298.257839303
+        middleRadius = 6378100.0
+        polarRadius = (1 - alphaEllips) * ekvRadius
+        c20 = -1.08262741e-3
+        eElips2 = 2 * alphaEllips - alphaEllips * alphaEllips
+    }
+
+    override fun blh2xyz(
+        b: Double,
+        l: Double,
+        h: Double,
+    ): Vector3D {
+        var v = Vector3D()
+        val sb = Math.sin(b)
+        val cb = Math.cos(b)
+        val sl = Math.sin(l)
+        val cl = Math.cos(l)
+        val n = ekvRadius / Math.sqrt(1 - eElips2 * sb * sb)
+        v.x = (n + h) * cb * cl
+        v.y = (n + h) * cb * sl
+        v.z = ((1 - eElips2) * n + h) * sb
+        return v
+    }
+
+    override fun xyz2blh(r: Vector3D): BLHPoint {
+        var l = atan2(r.y, r.x)
+        if (l < 0.0) l = PI * 2.0 + l
+
+        var r1: Double = hypot(r.x, r.y)
+        val rr: Double = hypot(r1, r.z)
+        val sinFiC: Double = r.z / rr
+        val cosFiC = r1 / rr
+        val fic: Double = atan2(r.z, r1)
+        val ekvadrat: Double = (2.0 - alphaEllips) * alphaEllips
+        r1 = 1.0 - ekvadrat * cosFiC * cosFiC
+        val delta = ekvadrat * sinFiC * cosFiC / r1
+        val rE: Double = ekvRadius * sqrt(1.0 - ekvadrat) / sqrt(r1)
+        val eMal = rE * delta / rr
+        var b = fic + eMal
+        var h = ((rr - rE) * (1.0 - eMal * delta / 2.0))
+        return BLHPoint(b, l, h)
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/utils/earth/EarthPZ90d02.kt

@@ -0,0 +1,14 @@
+package ballistics.utils.earth
+
+class EarthPZ90d02 : EarthPZ90d0() {
+    init {
+        fM = 398600.4418e9
+        ekvRadius = 6378136.0
+        wEarth = 0.7292115e-4
+        alphaEllips = 1 / 298.25784
+        middleRadius = 6378100.0
+        polarRadius = (1 - alphaEllips) * ekvRadius
+        c20 = -1.08262575e-3
+        eElips2 = 2 * alphaEllips - alphaEllips * alphaEllips
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/utils/math/Matrix3D.kt

@@ -0,0 +1,58 @@
+package ballistics.utils.math
+
+class Matrix3D(var first: Vector3D, var second: Vector3D, var third: Vector3D) {
+    constructor() : this(Vector3D(), Vector3D(), Vector3D())
+
+    fun firstColumn(): Vector3D {
+        return Vector3D(first.x, second.x, third.x)
+    }
+
+    fun secondColumn(): Vector3D {
+        return Vector3D(first.y, second.y, third.y)
+    }
+
+    fun thirdColumn(): Vector3D {
+        return Vector3D(first.z, second.z, third.z)
+    }
+
+    fun empty() {
+        first.setData(0.0, 0.0, 0.0)
+        second.setData(0.0, 0.0, 0.0)
+        third.setData(0.0, 0.0, 0.0)
+    }
+
+    fun makeOxMatrix(angle: Double) {
+        first.setData(1.0, 0.0, 0.0)
+        second.setData(0.0, Math.cos(angle), -Math.sin(angle))
+        third.setData(0.0, Math.sin(angle), Math.cos(angle))
+    }
+
+    fun makeOyMatrix(angle: Double) {
+        first.setData(Math.cos(angle), 0.0, Math.sin(angle))
+        second.setData(0.0, 1.0, 0.0)
+        third.setData(-Math.sin(angle), 0.0, Math.cos(angle))
+    }
+
+    fun makeOzMatrix(angle: Double) {
+        first.setData(Math.cos(angle), -Math.sin(angle), 0.0)
+        second.setData(Math.sin(angle), Math.cos(angle), 0.0)
+        third.setData(0.0, 0.0, 1.0)
+    }
+
+    fun transpose(): Matrix3D {
+        return Matrix3D(
+            Vector3D(first.x, second.x, third.x),
+            Vector3D(first.y, second.y, third.y),
+            Vector3D(first.z, second.z, third.z),
+        )
+    }
+
+    operator fun times(other: Vector3D) = Vector3D(first * other, second * other, third * other)
+
+    operator fun times(other: Matrix3D) =
+        Matrix3D(
+            Vector3D(first * other.firstColumn(), first * other.secondColumn(), first * other.thirdColumn()),
+            Vector3D(second * other.firstColumn(), second * other.secondColumn(), second * other.thirdColumn()),
+            Vector3D(third * other.firstColumn(), third * other.secondColumn(), third * other.thirdColumn()),
+        )
+}

libs/ballistics-lib/src/main/kotlin/ballistics/utils/math/Quaternion3D.kt

@@ -0,0 +1,195 @@
+package ballistics.utils.math
+
+import kotlin.math.acos
+import kotlin.math.cos
+import kotlin.math.sin
+import kotlin.math.sqrt
+
+class Quaternion3D(
+    var q0: Double,
+    var q1: Double,
+    var q2: Double,
+    var q3: Double,
+) {
+    constructor() : this(0.0, 0.0, 0.0, 0.0)
+
+    fun makeQuat(
+        angle: Double,
+        v: Vector3D,
+    ) {
+        var vv = v.basis()
+
+        q0 = cos(angle / 2.0)
+        q1 = vv.x * sin(angle / 2.0)
+        q2 = vv.y * sin(angle / 2.0)
+        q3 = vv.z * sin(angle / 2.0)
+
+        makeBasis()
+    }
+
+    /*!
+    норма кватерниона
+     */
+    fun norm() = sqrt(q0 * q0 + q1 * q1 + q2 * q2 + q3 * q3)
+
+    fun makeBasis() {
+        var n = norm()
+        q0 /= n
+        q1 /= n
+        q2 /= n
+        q3 /= n
+    }
+
+    fun basis(): Quaternion3D {
+        var n = norm()
+        return Quaternion3D(q0 / n, q1 / n, q2 / n, q3 / n)
+    }
+
+    fun makeOpposite() {
+        var n = norm()
+        q0 /= n
+        q1 /= -n
+        q2 /= -n
+        q3 /= -n
+    }
+
+    fun opposite(): Quaternion3D {
+        var n = norm()
+        return Quaternion3D(q0 / n, -q1 / n, -q2 / n, -q3 / n)
+    }
+
+    fun module() = acos(basis().q0) * 2.0
+
+    fun matrix(): Matrix3D {
+        var x2 = q1 + q1
+        var y2 = q2 + q2
+        var z2 = q3 + q3
+
+        var xx = q1 * x2
+        var yy = q2 * y2
+        var wx = q0 * x2
+
+        var xy = q1 * y2
+        var yz = q2 * z2
+        var wy = q0 * y2
+
+        var xz = q1 * z2
+        var zz = q3 * z2
+        var wz = q0 * z2
+
+        return Matrix3D(
+            Vector3D(1.0 - (yy + zz), xy + wz, xz - wy),
+            Vector3D(xy - wz, 1.0 - (xx + zz), yz + wx),
+            Vector3D(xz + wy, yz - wx, 1.0 - (xx + yy)),
+        ).transpose()
+    }
+
+    fun fromMatrixStanley(m: Matrix3D) {
+        val t = m.first.x + m.second.y + m.third.z
+
+        val valueq0 = (1.0 + t) / 4.0
+        val valueq1 = (1.0 + 2.0 * m.first.x - t) / 4.0
+        val valueq2 = (1.0 + 2.0 * m.second.y - t) / 4.0
+        val valueq3 = (1.0 + 2.0 * m.third.z - t) / 4.0
+
+        if ((valueq0 >= valueq1) && (valueq0 >= valueq2) && (valueq0 >= valueq3)) {
+            q0 = sqrt(valueq0)
+            q1 = (m.second.z - m.third.y) / (4.0 * q0)
+            q2 = (m.third.x - m.first.z) / (4.0 * q0)
+            q3 = (m.first.y - m.second.x) / (4.0 * q0)
+        }
+
+        if ((valueq1 >= valueq0) && (valueq1 >= valueq2) && (valueq1 >= valueq3)) {
+            q1 = sqrt(valueq1)
+            q0 = (m.second.z - m.third.y) / (4.0 * q1)
+            q2 = (m.first.y + m.second.x) / (4.0 * q1)
+            q3 = (m.third.x + m.first.z) / (4.0 * q1)
+        }
+
+        if ((valueq2 >= valueq0) && (valueq2 >= valueq1) && (valueq2 >= valueq3)) {
+            q2 = sqrt(valueq2)
+            q0 = (m.third.x - m.first.z) / (4.0 * q2)
+            q1 = (m.first.y + m.second.x) / (4.0 * q2)
+            q3 = (m.second.z + m.third.y) / (4.0 * q2)
+        }
+
+        if ((valueq3 >= valueq1) && (valueq3 >= valueq2) && (valueq3 >= valueq0)) {
+            q3 = sqrt(valueq3)
+            q0 = (m.first.y - m.second.x) / (4.0 * q3)
+            q1 = (m.third.x + m.first.z) / (4.0 * q3)
+            q2 = (m.second.z + m.third.y) / (4.0 * q3)
+        }
+
+        makeBasis()
+        makeOpposite()
+    }
+
+    operator fun times(v: Vector3D): Vector3D {
+        var p = opposite()
+        var r = Vector3D()
+
+        var a = (q0 + q1) * (v.x)
+        var b = (q3 - q2) * (v.y - v.z)
+        var c = (q1 - q0) * (v.y + v.z)
+        var d = (q2 + q3) * (v.x)
+        var e = (q1 + q3) * (v.x + v.y)
+        var f = (q1 - q3) * (v.x - v.y)
+        var g = (q0 + q2) * (-v.z)
+        var h = (q0 - q2) * (v.z)
+
+        var qq0 = b + (-e - f + g + h) * 0.5
+        var qq1 = a - (e + f + g + h) * 0.5
+        var qq2 = -c + (e - f + g - h) * 0.5
+        var qq3 = -d + (e - f - g + h) * 0.5
+
+        a = (qq0 + qq1) * (p.q0 + p.q1)
+        c = (qq1 - qq0) * (p.q2 + p.q3)
+        d = (qq2 + qq3) * (p.q1 - p.q0)
+        e = (qq1 + qq3) * (p.q1 + p.q2)
+        f = (qq1 - qq3) * (p.q1 - p.q2)
+        g = (qq0 + qq2) * (p.q0 - p.q3)
+        h = (qq0 - qq2) * (p.q0 + p.q3)
+
+        r.x = a - (e + f + g + h) * 0.5
+        r.y = -c + (e - f + g - h) * 0.5
+        r.z = -d + (e - f - g + h) * 0.5
+
+        return r
+    }
+
+    operator fun times(other: Quaternion3D): Quaternion3D {
+        var a = (q0 + q1) * (other.q0 + other.q1)
+        var b = (q3 - q2) * (other.q2 - other.q3)
+        var c = (q1 - q0) * (other.q2 + other.q3)
+        var d = (q2 + q3) * (other.q1 - other.q0)
+        var e = (q1 + q3) * (other.q1 + other.q2)
+        var f = (q1 - q3) * (other.q1 - other.q2)
+        var g = (q0 + q2) * (other.q0 - other.q3)
+        var h = (q0 - q2) * (other.q0 + other.q3)
+
+        return Quaternion3D(
+            b + (-e - f + g + h) * 0.5,
+            a - (e + f + g + h) * 0.5,
+            -c + (e - f + g - h) * 0.5,
+            -d + (e - f - g + h) * 0.5,
+        )
+    }
+
+    operator fun plus(other: Quaternion3D): Quaternion3D {
+        return Quaternion3D(
+            q0 + other.q0,
+            q1 + other.q1,
+            q2 + other.q2,
+            q3 + other.q3,
+        )
+    }
+
+    operator fun minus(other: Quaternion3D): Quaternion3D {
+        return Quaternion3D(
+            q0 - other.q0,
+            q1 - other.q1,
+            q2 - other.q2,
+            q3 - other.q3,
+        )
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/utils/math/Vector3D.kt

@@ -0,0 +1,70 @@
+package ballistics.utils.math
+
+import kotlin.math.sqrt
+
+class Vector3D(ax: Double, ay: Double, az: Double) {
+    var x = ax
+    var y = ay
+    var z = az
+
+    constructor() : this(0.0, 0.0, 0.0)
+
+    fun setData(
+        xVal: Double,
+        yVal: Double,
+        zVal: Double,
+    ) {
+        x = xVal
+        y = yVal
+        z = zVal
+    }
+
+    fun module() = sqrt(x * x + y * y + z * z)
+
+    fun basis(): Vector3D {
+        var v = Vector3D()
+        v.x = x / module()
+        v.y = y / module()
+        v.z = z / module()
+        return v
+    }
+
+    override fun toString(): String {
+        return "x = $x, y = $y, z = $z"
+    }
+
+    operator fun plus(other: Vector3D) = Vector3D(x + other.x, y + other.y, z + other.z)
+
+    operator fun minus(other: Vector3D) = Vector3D(x - other.x, y - other.y, z - other.z)
+
+    operator fun times(value: Double) = Vector3D(x * value, y * value, z * value)
+
+    operator fun times(other: Vector3D) = x * other.x + y * other.y + z * other.z
+
+    operator fun rem(other: Vector3D) =
+        Vector3D(
+            y * other.z - z * other.y,
+            z * other.x - x * other.z,
+            x * other.y - y * other.x,
+        )
+
+    override fun equals(other: Any?): Boolean {
+        if (this === other) return true
+        if (javaClass != other?.javaClass) return false
+
+        other as Vector3D
+
+        if (x != other.x) return false
+        if (y != other.y) return false
+        if (z != other.z) return false
+
+        return true
+    }
+
+    override fun hashCode(): Int {
+        var result = x.hashCode()
+        result = 31 * result + y.hashCode()
+        result = 31 * result + z.hashCode()
+        return result
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/utils/math/equations/AbstractEquationCalculator.kt

@@ -0,0 +1,25 @@
+package ballistics.utils.math.equations
+
+abstract class AbstractEquationCalculator {
+    var iteration = 0
+    var maxIterations = 20
+    var delta = 0.001
+    var value = 0.0
+
+    abstract fun calculate(
+        x1: Double,
+        x2: Double,
+        f: (Double) -> Double,
+    ): Double?
+
+    inline fun check(
+        x1: Double,
+        x2: Double,
+        f: (Double) -> Double,
+    ): Boolean {
+        val y1 = f(x1)
+        val y2 = f(x2)
+
+        return ((y1 - value) * (y2 - value) <= 0)
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/utils/math/equations/EquationCalculator2Div.kt

@@ -0,0 +1,48 @@
+package ballistics.utils.math.equations
+
+import kotlin.math.abs
+
+class EquationCalculator2Div : AbstractEquationCalculator() {
+    override fun calculate(
+        x1: Double,
+        x2: Double,
+        f: (Double) -> Double,
+    ): Double? {
+        if (!check(x1, x2, f)) {
+            return null
+        }
+        iteration = 0
+        var xl = x1
+        var xr = x2
+
+        var xmid = (xl + xr) / 2.0
+        var ymid = f(xmid)
+        var y1 = f(xl)
+
+        if ((y1 - value) * (ymid - value) < 0) {
+            xr = xmid
+        } else {
+            xl = xmid
+        }
+
+        var d = abs(ymid - value)
+        while ((d > delta) && (iteration < maxIterations)) {
+            xmid = (xl + xr) / 2.0
+            ymid = f(xmid)
+            y1 = f(xl)
+            if ((y1 - value) * (ymid - value) < 0) {
+                xr = xmid
+            } else {
+                xl = xmid
+            }
+
+            d = abs(ymid - value)
+            iteration++
+        }
+        return if (iteration < maxIterations) {
+            xmid
+        } else {
+            null
+        }
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/utils/math/equations/EquationCalculatorSpan.kt

@@ -0,0 +1,44 @@
+package ballistics.utils.math.equations
+
+import kotlin.math.abs
+
+class EquationCalculatorSpan : AbstractEquationCalculator() {
+
+    override fun calculate(
+        x1: Double,
+        x2: Double,
+        f: (Double) -> Double,
+    ): Double? {
+        if (!check(x1, x2, f))
+            return null
+
+        iteration = 0
+        var xl = x1
+        var xr = x2
+
+        var yl = f(xl)
+        var yr = f(xr)
+        var xmid = xl - ((xr - xl) / ((yr - value) - (yl - value))) * (yl - value)
+        var ymid = f(xmid)
+
+        var d = abs(ymid - value) * 10000.0
+        while ((d > delta * 10000.0) && (iteration < maxIterations)) {
+            if ((yl - value) * (ymid - value) < 0) {
+                xr = xmid
+                yr = ymid
+            } else {
+                xl = xmid
+                yl = ymid
+            }
+            xmid = xl - ((xr - xl) / ((yr - value) - (yl - value))) * (yl - value)
+            ymid = f(xmid)
+            d = abs(ymid - value) * 10000.0
+            iteration++
+        }
+        return if (iteration < maxIterations) {
+            xmid
+        } else {
+            null
+        }
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/utils/math/minmax/AbstractMinMaxCalculator.kt

@@ -0,0 +1,13 @@
+package ballistics.utils.math.minmax
+
+abstract class AbstractMinMaxCalculator {
+    var delta = 0.00001
+    var maxIterations = 150
+    var sign = 1 // min = 1, max = -1
+
+    abstract fun calculate(
+        x1: Double,
+        x2: Double,
+        f: (Double) -> Double,
+    ): Double?
+}

libs/ballistics-lib/src/main/kotlin/ballistics/utils/math/minmax/MinMaxCrushingCalculator.kt

@@ -0,0 +1,33 @@
+package ballistics.utils.math.minmax
+
+class MinMaxCrushingCalculator : AbstractMinMaxCalculator() {
+    override fun calculate(
+        x1: Double,
+        x2: Double,
+        f: (Double) -> Double,
+    ): Double? {
+        var step = 1.0
+        var iteration = 0
+        var xk: Double = (x1 + x2) / 2
+        var yk: Double
+        var yk2: Double
+        yk = f(xk)
+        var xk2: Double = xk + step
+        yk2 = f(xk2)
+        var d = xk2 - xk
+        if (d < 0) d *= -1.0
+
+        while ((d >= delta) && (iteration < maxIterations)) {
+            if (yk2 * sign > yk * sign) step = -step / 3
+            xk = xk2
+            yk = yk2
+            xk2 = xk + step
+            yk2 = f(xk2)
+            d = xk2 - xk
+            if (d < 0) d *= -1.0
+            iteration++
+        }
+
+        return if (iteration < maxIterations) xk2 else null
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/utils/math/minmax/MinMaxGoldenSectCalculator.kt

@@ -0,0 +1,50 @@
+package ballistics.utils.math.minmax
+
+import kotlin.math.sqrt
+
+class MinMaxGoldenSectCalculator : AbstractMinMaxCalculator() {
+    val tau = (sqrt(5.0) - 1) / 2
+
+    override fun calculate(
+        x1: Double,
+        x2: Double,
+        f: (Double) -> Double,
+    ): Double? {
+        var iteration = 0
+        var x: Double?
+        var xl = x1
+        var xr = x2
+        var x11: Double = xl + (1 - tau) * (xr - xl)
+        var x22: Double = xl + tau * (xr - xl)
+
+        var d = x22 - x11
+        if (d < 0) d *= -1.0
+
+        //x = x1
+
+        var y1 = f(x11)
+        var y2 = f(x22)
+
+        while (d > delta) {
+            if (sign * y1 > sign * y2) {
+                xl = x11
+                x11 = x22
+                y1 = y2
+                x22 = xl + tau * (xr - xl)
+                y2 = f(x22)
+            } else {
+                xr = x22
+                x22 = x11
+                y2 = y1
+                x11 = xl + (1 - tau) * (xr - xl)
+                y1 = f(x11)
+            }
+            d = x22 - x11
+            if (d < 0) d *= -1.0
+            iteration++
+        }
+
+        x = (x11 + x22) / 2
+        return if (iteration < maxIterations) x else null
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/zrv/ZRVStepperCalculator.kt

@@ -0,0 +1,340 @@
+package ballistics.zrv
+
+import ballistics.orbitalPoints.AbstractOrbPointsCalculator
+import ballistics.types.*
+import ballistics.utils.earth.getEarth
+import ballistics.utils.math.Vector3D
+import ballistics.utils.math.equations.EquationCalculator2Div
+import ballistics.utils.math.minmax.MinMaxGoldenSectCalculator
+import kotlin.math.PI
+import kotlin.math.atan
+import kotlin.math.atan2
+import kotlin.math.cos
+import kotlin.math.sin
+import kotlin.math.sqrt
+
+/**
+ * Класс расчета ЗРВ КА-ППИ итерационным методом
+ * Осуществляется проход по точкам орбиты с заданным шагом и формирование участков
+ * видимости КА-ППИ с уточнением моментов входа и выхода в/из зоны видимости
+ *
+ */
+internal class ZRVStepperCalculator(opc: AbstractOrbPointsCalculator, ppi: List<PPI>) {
+    var opc: AbstractOrbPointsCalculator = opc
+    private var ppiExt = mutableListOf<PPIExt>()
+    private var indPPI = 0
+    private var step = 60.0
+    private val eart = getEarth(opc.earthType)
+    var ppi = ppi
+    var zrv = mutableListOf<ZRV>()
+
+    /**
+     * Расчет радиус-вектора КА в пунктовой СК текущего ППИ
+     * ППИ определяется по индексу _indPPI
+     */
+    private fun gskToPointSK(point: OrbitalPoint): Vector3D {
+        val r1 =
+            Vector3D(
+                point.r.z - ppiExt[indPPI].params.zn,
+                point.r.x - ppiExt[indPPI].params.xn,
+                point.r.y - ppiExt[indPPI].params.yn,
+            )
+
+        val r4 = r1.y * ppiExt[indPPI].params.cl + r1.z * ppiExt[indPPI].params.sl
+        return Vector3D(
+            ppiExt[indPPI].params.cb * r1.x - ppiExt[indPPI].params.sb * r4,
+            ppiExt[indPPI].params.sb * r1.x + ppiExt[indPPI].params.cb * r4,
+            ppiExt[indPPI].params.cl * r1.z - ppiExt[indPPI].params.sl * r1.y,
+        )
+    }
+
+    /**
+     * Расчет параметров видимости КА-ППИ для заданного положения КА
+     * ППИ определяется по индексу _indPPI
+     */
+    private fun calculateViewParams(point: OrbitalPoint) {
+        ppiExt[indPPI].t2 = point.t
+        val r1 = point.v.x * ppiExt[indPPI].params.cl + point.v.y * ppiExt[indPPI].params.sl
+        ppiExt[indPPI].vy2 = point.v.z * ppiExt[indPPI].params.sb + r1 * ppiExt[indPPI].params.cb
+        val r = gskToPointSK(point)
+        ppiExt[indPPI].gam2 = atan2(r.y, sqrt(r.x * r.x + r.z * r.z))
+        ppiExt[indPPI].kaY = r.y
+        ppiExt[indPPI].a2 = atan2(r.z, r.x)
+        if (ppiExt[indPPI].a2 < 0) {
+            ppiExt[indPPI].a2 += 2 * PI
+        }
+    }
+
+    /**
+     * Расчет зон радио-видимости на зааднном интервале времени
+     */
+    fun calculate(
+        tn: Double,
+        tk: Double,
+    ): BallisticsError {
+        zrv.clear()
+        if (ppi.isEmpty()) return BallisticsError.EMPTY_PPI_LIST
+
+        var gam1 : Double
+        var gam2 : Double
+
+        step = 60.0
+
+        val eqc = EquationCalculator2Div()
+        val minmax = MinMaxGoldenSectCalculator()
+        minmax.sign = -1
+        minmax.delta = 0.0001
+
+        eqc.delta = 0.001 * PI / 180.0
+        eqc.maxIterations = 50
+        ppiExt.clear()
+        for (p in ppi) {
+            ppiExt.add(
+                PPIExt().apply {
+                    params.nip = p.ppiNum
+                    params.au = p.auNum
+                    params.sb = sin(p.lat)
+                    params.cb = cos(p.lat)
+                    params.sl = sin(p.long)
+                    params.cl = cos(p.long)
+                    elevMax = p.elevMax
+                    elevMin = p.elevMin
+                    val r = eart.blh2xyz(p.lat, p.long, p.height)
+                    params.xn = r.x
+                    params.yn = r.y
+                    params.zn = r.z
+                    params.bn = p.lat
+                    shadowMin = p.shadowMin
+                    shadowMax = p.shadowMax
+                    // !!!
+                },
+            )
+        }
+
+        val stepper = opc.getStepper()
+
+        val cnt = ppiExt.count()
+        var ct = tn
+
+        // первая точка интервала расчета
+        var point = stepper.calculate(tn)
+        if (point == null) {
+            return BallisticsError.STEPPER_ERROR
+        }
+        for (i in 0..cnt - 1) {
+            indPPI = i
+            calculateViewParams(point)
+            if (ppiExt[i].gam2 > ppiExt[i].elevMin) {
+                ppiExt[i].zoneIn.t = point.t
+                ppiExt[i].vit = point.vit
+            }
+        }
+        ct += step
+        while (ct < tk) {
+            point = stepper.calculate(ct)
+            if (point == null) {
+                return BallisticsError.STEPPER_ERROR
+            }
+            for (i in 0..cnt - 1) {
+                indPPI = i
+                ppiExt[i].gam1 = ppiExt[i].gam2
+                ppiExt[i].t1 = ppiExt[i].t2
+                ppiExt[i].vy1 = ppiExt[i].vy2
+                ppiExt[i].a1 = ppiExt[i].a2
+
+                calculateViewParams(point)
+
+                gam1 = ppiExt[i].gam1
+                gam2 = ppiExt[i].gam2
+                // начало зоны
+                if (((gam1 < ppiExt[i].elevMin) && (gam2 > ppiExt[i].elevMin)) ||
+                    (gam1 > ppiExt[i].elevMin && gam2 > ppiExt[i].elevMin && ppiExt[i].vit == 0)
+                ) {
+                    eqc.value = ppiExt[i].elevMin
+                    var t = eqc.calculate(ppiExt[i].t1, ppiExt[i].t2, this::elevationAngle)
+
+                    if (t != null) {
+                        var buf = calcZoneParams(t)
+                        var it = 0
+                        while (buf.elevation <= ppiExt[i].elevMin && it < 15) {
+                            t += 0.005
+                            buf = calcZoneParams(t)
+                            ++it
+                        }
+                        ppiExt[i].zoneIn.t = t
+                        ppiExt[i].vit = point.vit
+                    } else {
+                        println(
+                            """
+                            Ошибка поиска момента входа в ЗРВ. 
+                            ${ppiExt[i].gam1 * 180 / PI} - ${ppiExt[i].elevMin * 180 / PI} - ${ppiExt[i].gam2 * 180 / PI}
+                            """.trimIndent(),
+                        )
+                        ppiExt[i].zoneIn.t = ppiExt[i].t2
+                    }
+                }
+                // / максимум!
+                if ((ppiExt[i].vy1 * ppiExt[i].vy2 < 0) &&
+                    (ppiExt[i].gam1 > ppiExt[i].elevMin) &&
+                    (ppiExt[i].kaY > 0)
+                ) {
+                    val tmax = minmax.calculate(ppiExt[i].t1, ppiExt[i].t2, this::elevationAngle)
+                    if (tmax == null) {
+                        println("Ошибка при поиске максимума ЗРВ")
+                        ppiExt[i].zoneMax = calcZoneParams((ppiExt[i].t1 + ppiExt[i].t2) / 2)
+                    } else {
+                        ppiExt[i].zoneMax = calcZoneParams(tmax)
+                    }
+                }
+                // конец зоны
+                if ((
+                        (ppiExt[i].gam1 > ppiExt[i].elevMin) &&
+                            (gam2 < ppiExt[i].elevMin && ppiExt[i].vit != 0)
+                    ) ||
+                    (gam1 < ppiExt[i].elevMin && gam2 < ppiExt[i].elevMin && ppiExt[i].vit != 0)
+                ) {
+                    eqc.value = ppiExt[i].elevMin
+                    var t = eqc.calculate(ppiExt[i].t1, ppiExt[i].t2, this::elevationAngle)
+
+                    if (t != null) {
+                        var buf = calcZoneParams(t)
+                        var it = 0
+                        while (buf.elevation <= ppiExt[i].elevMin && it < 15) {
+                            t -= 0.005
+                            buf = calcZoneParams(t)
+                            ++it
+                        }
+                        // ppiExt[i].zoneOut.t = if (buf.elevation<=ppiExt[i].elevMin) t!! else t!! - 0.1
+                        ppiExt[i].zoneOut.t = t
+                    } else {
+                        println("Ошибка поиска момента выхода в ЗРВ.")
+                        ppiExt[i].zoneOut.t = ppiExt[i].t1
+                    }
+                    if (ppiExt[i].zoneOut.t * ppiExt[i].zoneIn.t > 1) {
+                        addZone(i)
+                    } else {
+
+                        println(
+                            "Ошибка расчета граничной точки зоны ${ppiExt[i].params.nip}    ${ppiExt[i].zoneIn.t}   ${
+                                ppiExt[i].zoneOut.t
+                            }   ",
+                        )
+                    }
+                    ppiExt[i].vit = 0
+                }
+            }
+            ct += step
+        }
+        return BallisticsError.OK
+    }
+
+    private fun addZone(ind: Int) {
+        // разбиение зоны при наличии затенений
+        if (ppiExt[ind].elevMax < ppiExt[ind].zoneMax.elevation ||
+            (ppiExt[ind].shadowMin != null && ppiExt[ind].shadowMin!!.toList().isNotEmpty()) ||
+            (ppiExt[ind].shadowMax != null && ppiExt[ind].shadowMax!!.toList().isNotEmpty())
+        ) {
+            var tn: Double? = null
+            var inZone = false
+            var inZoneNext = false
+            var vzp = calcZoneParams(ppiExt[ind].zoneIn.t)
+
+            if (ppiExt[ind].isVisible(vzp.azimuth, vzp.elevation)) {
+                tn = ppiExt[ind].zoneIn.t
+                inZone = true
+            }
+            val step = 0.1
+            var t = ppiExt[ind].zoneIn.t + step
+            while (t <= ppiExt[ind].zoneOut.t) {
+                vzp = calcZoneParams(t)
+                inZoneNext = ppiExt[ind].isVisible(vzp.azimuth, vzp.elevation)
+
+                // конец видимости
+                if (inZone && !inZoneNext) {
+                    val z =
+                        ZRV(
+                            ppiExt[ind].params.nip,
+                            ppiExt[ind].params.au,
+                            ppiExt[ind].vit,
+                            calcZoneParams(tn!!),
+                            VisibilityParametersZRV(
+                                ppiExt[ind].zoneMax.t,
+                                ppiExt[ind].zoneMax.range,
+                                ppiExt[ind].zoneMax.azimuth,
+                                ppiExt[ind].zoneMax.elevation,
+                            ),
+                            calcZoneParams(t - step),
+                        )
+                    zrv.add(z)
+                }
+                // начало видимости
+                if (!inZone && inZoneNext) {
+                    tn = t
+                }
+
+                inZone = inZoneNext
+                t += step
+            }
+            if (inZone && inZoneNext) {
+                val z =
+                    ZRV(
+                        ppiExt[ind].params.nip,
+                        ppiExt[ind].params.au,
+                        ppiExt[ind].vit,
+                        calcZoneParams(tn!!),
+                        VisibilityParametersZRV(
+                            ppiExt[ind].zoneMax.t,
+                            ppiExt[ind].zoneMax.range,
+                            ppiExt[ind].zoneMax.azimuth,
+                            ppiExt[ind].zoneMax.elevation,
+                        ),
+                        calcZoneParams(ppiExt[ind].zoneOut.t),
+                    )
+                zrv.add(z)
+            }
+        } else {
+            val z =
+                ZRV(
+                    ppiExt[ind].params.nip,
+                    ppiExt[ind].params.au,
+                    ppiExt[ind].vit,
+                    calcZoneParams(ppiExt[ind].zoneIn.t),
+                    calcZoneParams(ppiExt[ind].zoneMax.t),
+                    calcZoneParams(ppiExt[ind].zoneOut.t),
+                )
+//            z.zoneIn.elevation = ppiExt[ind].elevMin
+//            z.zoneOut.elevation = ppiExt[ind].elevMin
+            zrv.add(z)
+        }
+    }
+
+    fun clear(){
+        opc.clear()
+        zrv.clear()
+    }
+
+    private fun calcZoneParams(t: Double): VisibilityParametersZRV {
+        val z = VisibilityParametersZRV()
+
+        val stepper = opc.getStepper()
+        val point: OrbitalPoint? = stepper.calculate(t) ?: return z
+
+        val r = gskToPointSK(point!!)
+
+        z.t = t
+        z.elevation = atan(r.y / sqrt(r.x * r.x + r.z * r.z))
+        z.range = r.module()
+        z.azimuth = atan2(r.z, r.x)
+        if (z.azimuth < 0) {
+            z.azimuth += 2 * PI
+        }
+        return z
+    }
+
+    private fun elevationAngle(t: Double): Double {
+        val stepper = opc.getStepper()
+        val point = stepper.calculate(t)
+        val r = gskToPointSK(point!!)
+        return atan(r.y / sqrt(r.z * r.z + r.x * r.x))
+    }
+}

libs/ballistics-lib/src/main/kotlin/ballistics/zrv/ZRVStepperCalculatorAsync.kt

@@ -0,0 +1,363 @@
+package ballistics.zrv
+
+import ballistics.orbitalPoints.AbstractOrbPointsCalculator
+import ballistics.types.*
+import ballistics.utils.earth.getEarth
+import ballistics.utils.math.Vector3D
+import ballistics.utils.math.equations.EquationCalculator2Div
+import ballistics.utils.math.minmax.MinMaxGoldenSectCalculator
+import reactor.core.publisher.Flux
+import reactor.core.publisher.FluxSink
+import java.io.IOException
+import kotlin.math.PI
+import kotlin.math.atan
+import kotlin.math.atan2
+import kotlin.math.cos
+import kotlin.math.sin
+import kotlin.math.sqrt
+
+internal class ZRVStepperCalculatorAsync(
+    private val opc: AbstractOrbPointsCalculator,
+    ppi: List<PPI>,
+) {
+
+    private val ppiExt: List<PPIExt>
+    private var step = 60.0
+    private val earth = getEarth(opc.earthType)
+
+    init {
+        ppiExt =
+            ppi.map { p ->
+                PPIExt().apply {
+                    params.nip = p.ppiNum
+                    params.au = p.auNum
+                    params.sb = sin(p.lat)
+                    params.cb = cos(p.lat)
+                    params.sl = sin(p.long)
+                    params.cl = cos(p.long)
+                    elevMax = p.elevMax
+                    elevMin = p.elevMin
+                    val r = earth.blh2xyz(p.lat, p.long, p.height)
+                    params.xn = r.x
+                    params.yn = r.y
+                    params.zn = r.z
+                    params.bn = p.lat
+                    shadowMin = p.shadowMin
+                    shadowMax = p.shadowMax
+                }
+            }
+    }
+
+    private fun gskToPointSK(
+        point: OrbitalPoint,
+        ppi: PPIExt,
+    ): Vector3D {
+        val r1 =
+            Vector3D(
+                point.r.z - ppi.params.zn,
+                point.r.x - ppi.params.xn,
+                point.r.y - ppi.params.yn,
+            )
+        val r4 = r1.y * ppi.params.cl + r1.z * ppi.params.sl
+
+        return Vector3D(
+            ppi.params.cb * r1.x - ppi.params.sb * r4,
+            ppi.params.sb * r1.x + ppi.params.cb * r4,
+            ppi.params.cl * r1.z - ppi.params.sl * r1.y,
+        )
+    }
+
+    private fun elevationAngle(
+        t: Double,
+        ppi: PPIExt,
+    ): Double {
+        val stepper = opc.getStepper()
+        val point = stepper.calculate(t)
+        val r = gskToPointSK(point!!, ppi)
+
+        return atan(r.y / sqrt(r.z * r.z + r.x * r.x))
+    }
+
+    private fun calcZoneParams(
+        t: Double,
+        ppi: PPIExt,
+    ): VisibilityParametersZRV {
+        val z = VisibilityParametersZRV()
+        val stepper = opc.getStepper()
+        val point: OrbitalPoint = stepper.calculate(t) ?: return z
+        val r = gskToPointSK(point, ppi)
+
+        z.t = t
+        z.elevation = atan(r.y / sqrt(r.x * r.x + r.z * r.z))
+        z.range = r.module()
+        z.azimuth = atan2(r.z, r.x)
+
+        if (z.azimuth < 0) {
+            z.azimuth += 2 * PI
+        }
+
+        return z
+    }
+
+    private fun checkZone(ppi : PPIExt) =
+        (ppi.elevMax < ppi.zoneMax.elevation) ||
+                ((ppi.shadowMin?.count() ?: 0) != 0) ||
+                ((ppi.shadowMax?.count() ?: 0) != 0)
+
+
+
+    private fun parceZone(
+        sink: FluxSink<ZRV>,
+        ppi: PPIExt
+    ){
+        var tn: Double? = null
+        var inZone = false
+        var inZoneNext = false
+        var vzp = calcZoneParams(ppi.zoneIn.t, ppi)
+
+        if (ppi.isVisible(vzp.azimuth, vzp.elevation)) {
+            tn = ppi.zoneIn.t
+            inZone = true
+        }
+
+        val step = 0.1
+        var t = ppi.zoneIn.t + step
+
+        while (t <= ppi.zoneOut.t) {
+            vzp = calcZoneParams(t, ppi)
+            inZoneNext = ppi.isVisible(vzp.azimuth, vzp.elevation)
+
+            // конец видимости
+            if (inZone && !inZoneNext) {
+                sink.next(
+                    ZRV(
+                        ppi.params.nip,
+                        ppi.params.au,
+                        ppi.vit,
+                        calcZoneParams(tn!!, ppi),
+                        VisibilityParametersZRV(ppi.zoneMax.t, ppi.zoneMax.range, ppi.zoneMax.azimuth, ppi.zoneMax.elevation),
+                        calcZoneParams(t - step, ppi),
+                    ),
+                )
+            }
+
+            // начало видимости
+            if (!inZone && inZoneNext) {
+                tn = t
+            }
+
+            inZone = inZoneNext
+            t += step
+        }
+
+        if (inZone && inZoneNext) {
+            sink.next(
+                ZRV(
+                    ppi.params.nip,
+                    ppi.params.au,
+                    ppi.vit,
+                    calcZoneParams(tn!!, ppi),
+                    VisibilityParametersZRV(ppi.zoneMax.t, ppi.zoneMax.range, ppi.zoneMax.azimuth, ppi.zoneMax.elevation),
+                    calcZoneParams(ppi.zoneOut.t, ppi),
+                ),
+            )
+        }
+    }
+
+
+    private fun addZone(
+        sink: FluxSink<ZRV>,
+        ppi: PPIExt,
+    ) {
+        // разбиение зоны при наличии затенений
+        if ( checkZone(ppi) ) {
+            parceZone(sink, ppi)
+        } else {
+            sink.next(
+                ZRV(
+                    ppi.params.nip,
+                    ppi.params.au,
+                    ppi.vit,
+                    calcZoneParams(ppi.zoneIn.t, ppi),
+                    calcZoneParams(ppi.zoneMax.t, ppi),
+                    calcZoneParams(ppi.zoneOut.t, ppi),
+                ),
+            )
+        }
+    }
+
+    private fun calculateViewParams(
+        point: OrbitalPoint,
+        ppi: PPIExt,
+    ) {
+        ppi.t2 = point.t
+
+        val r1 = point.v.x * ppi.params.cl + point.v.y * ppi.params.sl
+        ppi.vy2 = point.v.z * ppi.params.sb + r1 * ppi.params.cb
+
+        val r = gskToPointSK(point, ppi)
+        ppi.gam2 = atan2(r.y, sqrt(r.x * r.x + r.z * r.z))
+        ppi.kaY = r.y
+        ppi.a2 = atan2(r.z, r.x)
+
+        if (ppi.a2 < 0) {
+            ppi.a2 += 2 * PI
+        }
+    }
+
+    fun calculate(
+        tn: Double,
+        tk: Double,
+    ): Flux<ZRV> {
+        if (ppiExt.isEmpty()) {
+            return Flux.error(IOException()) // replace with custom exception
+        }
+
+        var gam1 : Double
+        var gam2 : Double
+
+        step = 60.0
+
+        val eqc = EquationCalculator2Div()
+        val minmax = MinMaxGoldenSectCalculator()
+
+        minmax.sign = -1
+        minmax.delta = 0.0001
+
+        eqc.delta = 0.001 * PI / 180.0
+        eqc.maxIterations = 50
+
+        val stepper = opc.getStepper()
+        var ct = tn
+
+        return Flux.create { sink ->
+            // первая точка интервала расчета
+            var point = stepper.calculate(tn) ?: return@create sink.error(IOException())
+
+            ppiExt.forEach {
+                calculateViewParams(point, it)
+
+                if (it.gam2 > it.elevMin) {
+                    it.zoneIn.t = point.t
+                    it.vit = point.vit
+                }
+            }
+
+            ct += step
+            while (ct < tk) {
+                point = stepper.calculate(ct) ?: return@create sink.error(IOException())
+
+                ppiExt.forEach { innerPpiLoop ->
+                    innerPpiLoop.gam1 = innerPpiLoop.gam2
+                    innerPpiLoop.t1 = innerPpiLoop.t2
+                    innerPpiLoop.vy1 = innerPpiLoop.vy2
+                    innerPpiLoop.a1 = innerPpiLoop.a2
+
+                    calculateViewParams(point, innerPpiLoop)
+
+                    gam1 = innerPpiLoop.gam1
+                    gam2 = innerPpiLoop.gam2
+                    // начало зоны
+                    if (((gam1 < innerPpiLoop.elevMin) && (gam2 > innerPpiLoop.elevMin)) ||
+                        (gam1 > innerPpiLoop.elevMin && gam2 > innerPpiLoop.elevMin && innerPpiLoop.vit == 0)
+                    ) {
+                        eqc.value = innerPpiLoop.elevMin
+                        var t =
+                            eqc.calculate(innerPpiLoop.t1, innerPpiLoop.t2) {
+                                elevationAngle(it, innerPpiLoop)
+                            }
+
+                        if (t != null) {
+                            var buf = calcZoneParams(t, innerPpiLoop)
+                            var it = 0
+
+                            while (buf.elevation <= innerPpiLoop.elevMin && it < 15) {
+                                t += 0.005
+                                buf = calcZoneParams(t, innerPpiLoop)
+                                ++it
+                            }
+
+                            innerPpiLoop.zoneIn.t = t
+                            innerPpiLoop.vit = point.vit
+                        } else {
+                            println(
+                                """
+                                Ошибка поиска момента входа в ЗРВ. 
+                                ${innerPpiLoop.gam1 * 180 / PI} - ${innerPpiLoop.elevMin * 180 / PI} - ${innerPpiLoop.gam2 * 180 / PI}
+                                """.trimIndent(),
+                            )
+                            innerPpiLoop.zoneIn.t = innerPpiLoop.t2
+                        }
+                    }
+
+                    // / максимум!
+                    if ((innerPpiLoop.vy1 * innerPpiLoop.vy2 < 0) &&
+                        (innerPpiLoop.gam1 > innerPpiLoop.elevMin) &&
+                        (innerPpiLoop.kaY > 0)
+                    ) {
+                        val tmax =
+                            minmax.calculate(innerPpiLoop.t1, innerPpiLoop.t2) {
+                                elevationAngle(it, innerPpiLoop)
+                            }
+
+                        if (tmax == null) {
+                            println("Ошибка при поиске максимума ЗРВ")
+                            innerPpiLoop.zoneMax = calcZoneParams((innerPpiLoop.t1 + innerPpiLoop.t2) / 2, innerPpiLoop)
+                        } else {
+                            innerPpiLoop.zoneMax = calcZoneParams(tmax, innerPpiLoop)
+                        }
+                    }
+                    // конец зоны
+                    if ((
+                            (innerPpiLoop.gam1 > innerPpiLoop.elevMin) &&
+                                (gam2 < innerPpiLoop.elevMin && innerPpiLoop.vit != 0)
+                        ) ||
+                        (
+                            gam1 < innerPpiLoop.elevMin && gam2 < innerPpiLoop.elevMin &&
+                                innerPpiLoop.vit != 0
+                        )
+                    ) {
+                        eqc.value = innerPpiLoop.elevMin
+
+                        var t =
+                            eqc.calculate(innerPpiLoop.t1, innerPpiLoop.t2) {
+                                elevationAngle(it, innerPpiLoop)
+                            }
+
+                        if (t != null) {
+                            var buf = calcZoneParams(t, innerPpiLoop)
+                            var it = 0
+
+                            while (buf.elevation <= innerPpiLoop.elevMin && it < 15) {
+                                t -= 0.005
+                                buf = calcZoneParams(t, innerPpiLoop)
+                                ++it
+                            }
+
+                            innerPpiLoop.zoneOut.t = t
+                        } else {
+                            println("Ошибка поиска момента выхода в ЗРВ.")
+                            innerPpiLoop.zoneOut.t = innerPpiLoop.t1
+                        }
+
+                        if (innerPpiLoop.zoneOut.t * innerPpiLoop.zoneIn.t > 1) {
+                            addZone(sink, innerPpiLoop)
+                        } else {
+                            println(
+                                "Ошибка расчета граничной точки зоны ${innerPpiLoop.params.nip}    ${innerPpiLoop.zoneIn.t}   ${
+                                    innerPpiLoop.zoneOut.t
+                                }   ",
+                            )
+                        }
+
+                        innerPpiLoop.vit = 0
+                    }
+                }
+
+                ct += step
+            }
+
+            sink.complete()
+        }
+    }
+}

libs/ballistics-lib/src/test/kotlin/Test.kt

@@ -0,0 +1,165 @@
+import ballistics.Ballistics
+import ballistics.types.InitialConditions
+import ballistics.types.IntegrationType
+import ballistics.types.ModDVType
+import ballistics.types.OrbitalPoint
+import ballistics.utils.fromDateTime
+import ballistics.utils.math.Vector3D
+import ballistics.utils.toDateTime
+import java.time.LocalDateTime
+import java.time.ZoneOffset
+
+
+fun main() {
+
+    var ts = System.currentTimeMillis()
+    val r = Ballistics()
+
+    r.modDVType = ModDVType.BARS // BARS
+    r.integrationType = IntegrationType.ADAMS7
+
+    val mnu = mutableListOf<InitialConditions>()
+    mnu.add(
+        InitialConditions(
+            OrbitalPoint(
+                LocalDateTime.of(2021, 12, 3, 0, 5, 45, 0).toEpochSecond(ZoneOffset.UTC).toDouble(),
+                1,
+                Vector3D(
+                    -5357933.7872,
+                    4328646.1787,
+                    0.0,
+                ),
+                Vector3D(
+                    941.8995373,
+                    1150.6919822,
+                    7544.9920840,
+                ),
+            ),
+            0.0061,
+            125.0,
+        ),
+    )
+
+    /*
+    mnu.add(InitialConditions(OrbitalPoint(1681135417.0,
+        1499,
+        6782429.138044466,
+     -1205586.994919729,
+     0.0,
+     -272.89561674108205,
+    -1460.1016031205972,
+     7539.209413654664),
+     0.0,
+     100.0)
+    )
+
+
+    mnu.add(InitialConditions(OrbitalPoint(1681203547.0,
+        1511,
+        2796871.797932365,
+        6295220.245166648,
+        0.0,
+        1351.432690993002,
+        -615.1950685524682,
+        7539.420044409451),
+        0.0,
+    100.0)
+    )
+
+    mnu.add(InitialConditions(OrbitalPoint(1681288710.0,
+        1526,
+        2219522.8066600505,
+        6520960.574145492,
+        0.0,
+        1401.0525909402313,
+        -491.90052480161194,
+        7539.654233369637),
+        0.0,
+        100.0)
+    )
+
+    mnu.add(InitialConditions(OrbitalPoint(1681373873.0,
+        1541,
+        1624190.3979707498,
+        6693858.330353258,
+        0.0,
+        1439.4621272728175,
+        -364.5959436771924,
+        7539.927000039702),
+        0.0,
+        100.0)
+    )
+
+    mnu.add(InitialConditions(OrbitalPoint(1681464714.0,
+        1557,
+        3663882.4989639986,
+        5832354.67692114,
+        0.0,
+        1249.158953470858,
+        -802.9876366520331,
+        7540.375561571017),
+        0.0,
+        100.0)
+    )
+
+    mnu.add(InitialConditions(OrbitalPoint(1681538522.0,
+        1570,
+        -2393483.7360840607,
+        6458206.838182452,
+        0.0,
+        1398.0568057512764,
+        501.055477122656,
+        7540.680958625912),
+        0.0,
+        100.0)
+    )
+
+     */
+
+    mnu.sortBy { it.point.t }
+
+    // r.calculateOrbPoints(mnu.toTypedArray(), mnu.first().point.t, mnu.last().point.t+86400.0 * 5 )
+    r.calculateOrbPoints(mnu.toTypedArray(), mnu.first().point.t, mnu.first().point.t + 86400.0 * 5)
+
+    println("время расчета = ${System.currentTimeMillis() - ts}мс")
+
+//    for (p in r.points)
+//        println("${p.vit}  ${LocalDateTime.ofEpochSecond(p.t.toLong(),0, ZoneOffset.UTC)}  ${p.r.x}  ${p.r.y}  ${p.r.z}  ${p.v.x}  ${p.v.y}  ${p.v.z}")
+//
+
+    ts = System.currentTimeMillis()
+    r.calculateFlightLine(mnu.first().point.t, mnu.first().point.t + 86400.0 * 5)
+    println("время расчета = ${System.currentTimeMillis() - ts}мс")
+
+//       var p = r.points.last()
+//      println("${p.vit}  ${LocalDateTime.ofEpochSecond(p.t.toLong(),0, ZoneOffset.UTC)}  ${p.x}  ${p.y}  ${p.z}  ${p.vx}  ${p.vy}  ${p.vz}")
+//
+//       var st = r.getStepper()
+//       st?.let{
+//           var p = it.calculate(nu.point.t + 120.0)//452.4)
+//           p?.let {
+//               println(
+//                   "${p.vit}  ${
+//                       LocalDateTime.ofEpochSecond(
+//                           p.t.toLong(),
+//                           0,
+//                           ZoneOffset.UTC
+//                       )
+//                   }  ${p.x}  ${p.y}  ${p.z}  ${p.vx}  ${p.vy}  ${p.vz}"
+//               )
+//           }
+//       } ?: println("is null")
+//
+// //
+// //    val ts2 = System.currentTimeMillis()
+//    var ppi = mutableListOf<PPI>(PPI(1000, 0, 0.36617007706841037, 1.8417344524829609,
+//        0.0, 0.08726646259971647, 1.5707963267948966))
+//    try{
+//        r.calculateZRV(ppi,  mnu.first().point.t, mnu.last().point.t+86400.0 * 5)
+//        println(r.zrv.size)
+//    }
+//    catch(e : Exception){
+//        println("err")
+//    }
+//    println("время расчета = ${System.currentTimeMillis()-ts}мс  ЗРВ : ${r.zrv.size}")
+}

libs/ballistics-lib/src/test/kotlin/ballistics/BallisticsTest.kt

@@ -0,0 +1,86 @@
+package ballistics
+
+import ballistics.orbitalPoints.timeStepper.TLEStepper
+import ballistics.types.BallisticsError
+import ballistics.types.EarthType
+import ballistics.types.InitialConditions
+import ballistics.types.OPKatObj
+import ballistics.types.TLE
+import org.junit.jupiter.api.Assertions.assertEquals
+import org.junit.jupiter.api.Assertions.assertTrue
+import org.junit.jupiter.api.Test
+import java.time.LocalDateTime
+import java.time.ZoneOffset
+import kotlin.math.PI
+
+internal class BallisticsTest {
+    @Test
+    fun parseTLE() {
+        var s1 = "1 51824U 22019A   22058.83400146  .00001503  00000-0  10000-3 0  9990"
+        var s2 = "2 51824  97.5213 135.6262 0019439 268.1769  91.7238 15.08515864   114"
+
+        var s = TLEStepper(s1, s2, EarthType.PZ90d02)
+        var r = s.baseEpoch
+        println("${ LocalDateTime.ofEpochSecond(r.toLong(),(r % 1 * 1e9).toInt(), ZoneOffset.UTC)} ")
+
+        var p = s.calculate(r)
+
+        println(
+            "${p.vit}  ${ p.let { LocalDateTime.ofEpochSecond(p.t.toLong(),(p.t % 1 * 1e9).toInt(), ZoneOffset.UTC)} }  " +
+                "${p.r.x} " +
+                " ${p.r.y}  " +
+                "${p.r.z} " +
+                " ${p.v.x}  " +
+                "${p.v.y} " +
+                " ${p.v.z}",
+        )
+    }
+
+    @Test
+    fun testTLE() {
+        val tle =
+//            TLE(
+//                "1 51824U 22019A   22058.83400146  .00001503  00000-0  10000-3 0  9990",
+//                "2 51824  97.5213 135.6262 0019439 268.1769  91.7238 15.08515864   114",
+//            )
+
+            TLE(
+                "1 45621U 11037Q   22018.21004500  .07428171  00000-0  19069+0 0  9999",
+                "2 45621  51.4607 277.9348 0211635  34.9395 326.5259 15.14655454 76661",
+            )
+
+
+        val bal1 = Ballistics()
+        val bal2 = Ballistics()
+
+        bal1.sunAngleMin = -90 * PI / 180
+        bal2.sunAngleMin = -90 * PI / 180
+
+        val op = bal1.parseTLE(tle)
+
+        val tn = op.t
+        val tk = tn + 86400 * 2
+
+        val r = bal1.calculateOrbPoints(tle, tn, tk)
+        assertTrue(r == BallisticsError.OK)
+        bal2.calculateOrbPoints(InitialConditions(op, 0.0, 0.0), tn, tk)
+
+        bal1.calculateFlightLine(tn, tk)
+        bal2.calculateFlightLine(tn, tk)
+
+        val objs = listOf<OPKatObj>(OPKatObj(1, 1, "t",1, 60 * PI / 180, 20 * PI / 180, 100.0, 0.0, 90.0, -90.0))
+
+        bal1.calculateMPL(tn, tk, objs)
+        bal2.calculateMPL(tn, tk, objs)
+
+        println("${bal1.mpl.count()} - ${bal2.mpl.count()}")
+
+        assertEquals(bal1.mpl.count(), bal2.mpl.count())
+
+        if (bal1.mpl.toList().isNotEmpty()) {
+            val v1 = bal1.mpl.first()
+            val v2 = bal2.mpl.first()
+            assertTrue(v1.traverz - v2.traverz < 60)
+        }
+    }
+}

libs/ballistics-lib/src/test/kotlin/ballistics/flightLine/FlightLineCalculatorTest.kt

@@ -0,0 +1,134 @@
+package ballistics.flightLine
+
+import ballistics.Ballistics
+import ballistics.types.BallisticsError
+import ballistics.types.EarthType
+import ballistics.types.InitialConditions
+import ballistics.types.IntegrationType
+import ballistics.types.ModDVType
+import ballistics.types.OrbitalPoint
+import ballistics.types.WorkCSType
+import ballistics.utils.math.Vector3D
+import org.junit.jupiter.api.Assertions.assertEquals
+import org.junit.jupiter.api.Test
+import java.time.LocalDateTime
+import java.time.ZoneOffset
+import kotlin.math.PI
+
+internal class FlightLineCalculatorTest {
+    @Test
+    fun testFL() {
+        var r = Ballistics()
+
+        r.earthType = EarthType.PZ90d02
+        r.modDVType = ModDVType.FOTO
+        r.integrationType = IntegrationType.ADAMS7
+        r.workCoordinateSystem = WorkCSType.WCSOrbit
+        r.rollMax = 30 * PI / 180
+
+        var mnu = mutableListOf<InitialConditions>()
+        var tt1 = LocalDateTime.of(2023, 6, 20, 4, 18, 40, 180000000)
+        var tt2 = LocalDateTime.of(2023, 6, 19, 3, 17, 49, 251000000)
+
+        mnu.add(
+            InitialConditions(
+                OrbitalPoint(
+                    tt1.toEpochSecond(ZoneOffset.UTC).toDouble(), // 1687238300.0,
+                    18721,
+                    Vector3D(
+                        -6603039.949152646,
+                        -1870023.1481177735,
+                        0.0,
+                    ),
+                    Vector3D(
+                        -401.5289935173786,
+                        1413.4309568106378,
+                        7547.69556621523,
+                    ),
+                ),
+                0.014542,
+                125.0,
+            ),
+        )
+        mnu.add(
+            InitialConditions(
+                OrbitalPoint(
+                    tt2.toEpochSecond(ZoneOffset.UTC).toDouble(), // 1687238300.0,
+                    18705,
+                    Vector3D(
+                        -5880984.961878717,
+                        -3535801.795084351,
+                        0.0,
+                    ),
+                    Vector3D(
+                        -759.5481166352387,
+                        1257.702405545223,
+                        7548.463786464257,
+                    ),
+                ),
+                0.014542,
+                125.0,
+            ),
+        )
+
+        mnu.sortBy { it.point.t }
+
+        var t = LocalDateTime.now().toEpochSecond(ZoneOffset.UTC)
+        var rez = r.calculateOrbPoints(mnu.toTypedArray(), mnu.first().point.t, mnu.first().point.t + 86400.0 * 3)
+        var t2 = LocalDateTime.now().toEpochSecond(ZoneOffset.UTC)
+        println("расчет точек орибты за ${t2 - t} с")
+
+//        for (fl in r.points)
+//            println("${fl.vit}  ${LocalDateTime.ofEpochSecond(fl.t.toLong(),0, ZoneOffset.UTC)} ${fl.r}")
+
+        if (rez == BallisticsError.OK) {
+            t = LocalDateTime.now().toEpochSecond(ZoneOffset.UTC)
+            var rezfl = r.calculateFlightLine(mnu.first().point.t, mnu.first().point.t + 86400.0 * 3)
+            t2 = LocalDateTime.now().toEpochSecond(ZoneOffset.UTC)
+            println("расчет ТП и ПО за ${t2 - t} с")
+            for (fl in r.flightLine) {
+                println(
+                    "${fl.vit}  ${LocalDateTime.ofEpochSecond(fl.t.toLong(),0, ZoneOffset.UTC)}  " +
+                        " ${fl.leftOuterSwath.lat * 180.0 / PI}  ${fl.leftOuterSwath.long * 180 / PI}   " +
+                        " ${fl.flightLine.lat * 180.0 / PI}  ${fl.flightLine.long * 180 / PI}   " +
+                        " ${fl.rightOuterSwath.lat * 180.0 / PI}  ${fl.rightOuterSwath.long * 180 / PI}",
+                )
+            }
+//            assertEquals(/* expected = */BallisticsError.OK, /* actual = */rezfl )
+        }
+
+        /*var opKatObj = mutableListOf<OPKatObj>()
+        var b = -PI/ 2
+        var l = 0.0
+        while (b <= PI /2){
+            while (l <= PI * 2){
+                opKatObj.add(OPKatObj(1,1,1,b,l,0.0,0.0,PI/2,0.0))
+                l += 0.2 * PI / 180
+            }
+            b += 0.2 * PI / 180
+            l = 0.0
+        }
+        println(opKatObj.count())
+        t = LocalDateTime.now().toEpochSecond(ZoneOffset.UTC)
+
+        rez= r.calculateMPL(mnu.first().point.t, mnu.first().point.t+86400.0 * 3, opKatObj)
+
+        t2 = LocalDateTime.now().toEpochSecond(ZoneOffset.UTC)
+
+        println(r.mpl.count())
+
+        println("расчет МПЛ за ${t2-t}")*/
+
+//        println("$rez")
+//        for (t in r.mpl)
+//            println("${t.vit}  ${LocalDateTime.ofEpochSecond(t.traverz.toLong(),(t.traverz %1 * 1e9).toInt(), ZoneOffset.UTC)}  ${t.orientation.kren * 180.0 / PI}" +
+//                    " ${t.range}  ${t.sunAngle * 180 / PI}")
+
+        assertEquals(
+            // expected =
+            BallisticsError.OK,
+            // actual =
+            rez,
+        )
+    }
+}

libs/ballistics-lib/src/test/kotlin/ballistics/orbitalPoints/OrbitalPointsIntegratorTest.kt

@@ -0,0 +1,115 @@
+package ballistics.orbitalPoints
+
+import ballistics.Ballistics
+import ballistics.types.*
+import ballistics.utils.math.Vector3D
+import org.junit.jupiter.api.Assertions.assertEquals
+import org.junit.jupiter.api.Assertions.assertTrue
+import org.junit.jupiter.api.Test
+import java.time.LocalDateTime
+import java.time.ZoneOffset
+import kotlin.math.PI
+
+internal class OrbitalPointsIntegratorTest {
+    @Test
+    fun testMDBARS() {
+        val r = Ballistics()
+
+        r.earthType = EarthType.PZ90d02
+        r.modDVType = ModDVType.BARS
+        r.integrationType = IntegrationType.ADAMS7
+
+        val mnu = mutableListOf<InitialConditions>()
+        mnu.add(
+            InitialConditions(
+                OrbitalPoint(
+                    LocalDateTime.of(2021, 12, 3, 0, 5, 45, 0).toEpochSecond(ZoneOffset.UTC).toDouble(),
+                    1,
+                    Vector3D(
+                        -5357933.7872,
+                        4328646.1787,
+                        0.0,
+                    ),
+                    Vector3D(
+                        941.8995373,
+                        1150.6919822,
+                        7544.9920840,
+                    ),
+                ),
+                0.0061,
+                125.0,
+            ),
+        )
+        val rez = r.calculateOrbPoints(mnu.toTypedArray(), mnu.first().point.t, mnu.first().point.t + 86400.0 * 3)
+
+        for (p in r.revolutions) {
+            println(
+                "${p.vuz.vit}  ${LocalDateTime.ofEpochSecond(p.vuz.t.toLong(),0, ZoneOffset.UTC)}  " +
+                    "${p.lVuz * 180.0 / PI} " +
+                    " ${p.hVuz}  " +
+                    "${p.vuz.r.x} " +
+                    " ${p.vuz.r.y}  " +
+                    "${p.vuz.r.z} " +
+                    " ${p.vuz.v.x}  " +
+                    "${p.vuz.v.y} " +
+                    " ${p.vuz.v.z}",
+            )
+        }
+
+        val p = r.points.last()
+        println(
+            "${p.vit}  ${LocalDateTime.ofEpochSecond(
+                p.t.toLong(),
+                0,
+                ZoneOffset.UTC,
+            )}  ${p.r.x}  ${p.r.y}  ${p.r.z}  ${p.v.x}  ${p.v.y}  ${p.v.z}",
+        )
+
+        val rc = Vector3D(r.points.last().r.x, r.points.last().r.y, r.points.last().r.z)
+        val vc = Vector3D(r.points.last().v.x, r.points.last().v.y, r.points.last().v.z)
+
+        val rex = Vector3D(4723743.916, -4261060.842, -2647685.695)
+        val vex = Vector3D(-3168.702, 797.721, -6953.781)
+
+        val dr = (rc - rex).module()
+        val dv = (vex - vc).module()
+
+        println(dr)
+        println(dv)
+
+
+        assertEquals(
+            BallisticsError.OK,
+            rez,
+        )
+        assertTrue(
+            dr < 0.01,
+            "Ошибка прогноза по радиус-вектору",
+        )
+        assertTrue(dv < 0.01, "ошибка прогноза по вектору скорости")
+    }
+
+
+//    fun initOrekitData() {
+//        val dataDir = File("src/main/resources/orekit-data-main") // or absolute path
+//        if (!dataDir.exists()) {
+//            throw IllegalStateException("Orekit data directory not found at ${dataDir.absolutePath}")
+//        }
+//        val crawler = DirectoryCrawler(dataDir)
+//        DataProvidersManager.getInstance().addProvider(crawler)
+//    }
+
+//    @Test
+//    fun orekitTest(){
+//
+//        initOrekitData()
+//
+//        val l1 = "1 00666U 18111A   25349.99139990  .00002772  00000+0  11883-3 0  9999"
+//        val l2 = "2 43876  97.3687 252.5788 0001999  83.2822 276.8643 15.23576679386900"
+//
+//        val tle = TLE(l1, l2)
+//        val sat = Satellite(org.nstart.dep265.tletools.zeptomoby.core.TLE("",l1, l2))
+//
+//        assertTrue { sat.orbit.satId.toInt() == tle.satelliteNumber }
+//    }
+}

libs/ballistics-lib/src/test/kotlin/ballistics/types/ModDVTypeTest.kt

@@ -0,0 +1,20 @@
+package ballistics.types
+
+import org.junit.jupiter.api.Assertions.assertEquals
+import org.junit.jupiter.api.Test
+
+internal class ModDVTypeTest {
+    /*!
+        Проветка приведения типов от Int к enum
+     */
+    @Test
+    fun testMDFoto() {
+        val md1 = ModDVType.fromInt(1)
+        val md2 = ModDVType.fromInt(16)
+        val md3 = ModDVType.fromInt(100)
+
+        assertEquals(ModDVType.FOTO, md1)
+        assertEquals(ModDVType.BARS, md2)
+        assertEquals(ModDVType.FOTO, md3)
+    }
+}

libs/ballistics-lib/src/test/kotlin/ballistics/zrv/ZRVStepperCalculatorTest.kt

@@ -0,0 +1,103 @@
+package ballistics.zrv
+
+import ballistics.Ballistics
+import ballistics.types.BallisticsError
+import ballistics.types.EarthType
+import ballistics.types.InitialConditions
+import ballistics.types.IntegrationType
+import ballistics.types.ModDVType
+import ballistics.types.OrbitalPoint
+import ballistics.types.PPI
+import ballistics.types.TLE
+import ballistics.utils.math.Vector3D
+import org.junit.jupiter.api.Assertions.assertTrue
+import org.junit.jupiter.api.Test
+import java.time.LocalDateTime
+import java.time.ZoneOffset
+import kotlin.math.PI
+
+internal class ZRVStepperCalculatorTest {
+    @Test
+    fun calculate() {
+        val r = Ballistics()
+
+
+        val p = r.getTLEParams(TLE("1 62141U 24224C   25324.24196334  .00027917  00000-0  10193-2 0  9998",
+            "2 62141  42.9988  56.0174 0001022  26.8436 333.2469 15.27609213 55582"), "test")
+
+
+
+
+
+        r.earthType = EarthType.PZ90d02
+        r.modDVType = ModDVType.BARS
+        r.integrationType = IntegrationType.ADAMS7
+
+        val mnu = mutableListOf<InitialConditions>()
+        mnu.add(
+            InitialConditions(
+                OrbitalPoint(
+                    LocalDateTime.of(2021, 12, 3, 0, 5, 45, 0).toEpochSecond(ZoneOffset.UTC).toDouble(),
+                    1,
+                    Vector3D(
+                        -5357933.7872,
+                        4328646.1787,
+                        0.0,
+                    ),
+                    Vector3D(
+                        941.8995373,
+                        1150.6919822,
+                        7544.9920840,
+                    ),
+                ),
+                0.0061,
+                125.0,
+            ),
+        )
+        var rez = r.calculateOrbPoints(mnu.toTypedArray(), mnu.first().point.t, mnu.first().point.t + 86400.0 * 3)
+        if (rez == BallisticsError.OK) {
+            rez =
+                r.calculateZRV(
+                    listOf(
+                        PPI(
+                            1,
+                            1,
+                            55.15 * PI / 180,
+                            38.39 * PI / 180,
+                            0.0,
+                            5 * PI / 180,
+                            80 * PI / 180,
+//                shadowMin = listOf(
+//                    ShadowAU(110.0 * PI / 180, 150 * PI / 180, 25 * PI / 180),
+//                    ShadowAU(175.0 * PI / 180, 183 * PI / 180, 90 * PI / 180),
+//                ),
+//                shadowMax = listOf(
+//                    ShadowAU(189.0 * PI / 180, 192 * PI / 180, 0 * PI / 180),
+//                )
+                        ),
+                    ),
+                    mnu.first().point.t,
+                    mnu.first().point.t + 86400.0 * 3,
+                )
+
+            if (rez == BallisticsError.OK) {
+                for (z in r.zrv)
+                    println(
+                        " ${z.vit} \t" +
+                            " ${z.au} \t" +
+                            " ${LocalDateTime.ofEpochSecond(z.zoneIn.t.toLong(), 0, ZoneOffset.UTC)}  \t" +
+                            " ${z.zoneIn.azimuth * 180 / PI} \t" +
+                            " ${z.zoneIn.elevation * 180 / PI} \t" +
+                            " ${LocalDateTime.ofEpochSecond(z.zoneOut.t.toLong(), 0, ZoneOffset.UTC)} \t" +
+                            " ${z.zoneOut.azimuth * 180 / PI} \t" +
+                            " ${z.zoneOut.elevation * 180 / PI} \t" +
+                            " ${LocalDateTime.ofEpochSecond(z.zoneMax.t.toLong(), 0, ZoneOffset.UTC)} \t" +
+                            " ${z.zoneMax.azimuth * 180 / PI} \t" +
+                            " ${z.zoneMax.elevation * 180 / PI} \t",
+                    )
+            }
+        }
+
+        assertTrue((rez == BallisticsError.OK))
+    }
+}