diff --git a/app/build.gradle b/app/build.gradle index e5a84ac6..a024cd57 100644 --- a/app/build.gradle +++ b/app/build.gradle @@ -51,6 +51,7 @@ android { } dependencies { + implementation project(':core') implementation "com.google.android.material:material:$material_version" implementation "androidx.constraintlayout:constraintlayout:$constraint_layout_version" implementation "androidx.lifecycle:lifecycle-livedata-ktx:$lifecycle_version" diff --git a/build.gradle b/build.gradle index 6a891745..a2b6403a 100644 --- a/build.gradle +++ b/build.gradle @@ -1,13 +1,14 @@ buildscript { ext { gradle_version = '4.1.3' - gradle_plugin_version = '1.4.32' + kotlin_version = '1.4.32' + coroutines_version = '1.4.1' material_version = '1.3.0' constraint_layout_version = '2.0.4' lifecycle_version = '2.3.1' navigation_version = '2.3.5' preference_version = '1.1.1' - room_version = '2.2.6' + room_version = '2.3.0' hilt_version = '2.33-beta' retrofit_version = '2.9.0' predict4java_version = '1.3.1' @@ -24,7 +25,7 @@ buildscript { dependencies { classpath "com.android.tools.build:gradle:$gradle_version" classpath "com.google.dagger:hilt-android-gradle-plugin:$hilt_version" - classpath "org.jetbrains.kotlin:kotlin-gradle-plugin:$gradle_plugin_version" + classpath "org.jetbrains.kotlin:kotlin-gradle-plugin:$kotlin_version" } } diff --git a/core/.gitignore b/core/.gitignore new file mode 100644 index 00000000..42afabfd --- /dev/null +++ b/core/.gitignore @@ -0,0 +1 @@ +/build \ No newline at end of file diff --git a/core/build.gradle b/core/build.gradle new file mode 100644 index 00000000..d217f9bb --- /dev/null +++ b/core/build.gradle @@ -0,0 +1,14 @@ +plugins { + id 'java-library' + id 'kotlin' +} + +java { + sourceCompatibility = JavaVersion.VERSION_1_8 + targetCompatibility = JavaVersion.VERSION_1_8 +} + +dependencies { + implementation "org.jetbrains.kotlin:kotlin-stdlib-jdk8:$kotlin_version" + implementation "org.jetbrains.kotlinx:kotlinx-coroutines-android:$coroutines_version" +} \ No newline at end of file diff --git a/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/DeepSpaceSat.kt b/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/DeepSpaceSat.kt new file mode 100644 index 00000000..c7f21e57 --- /dev/null +++ b/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/DeepSpaceSat.kt @@ -0,0 +1,849 @@ +/* + * Look4Sat. Amateur radio satellite tracker and pass predictor. + * Copyright (C) 2019-2021 Arty Bishop (bishop.arty@gmail.com) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ +package com.rtbishop.look4sat.predict4kotlin + +import kotlin.math.* + +class DeepSpaceSat(tle: TLE) : Satellite(tle) { + + private val c1: Double + private val c4: Double + private val x1mth2: Double + private val x3thm1: Double + private val xlcof: Double + private val xnodcf: Double + private val t2cof: Double + private val aycof: Double + private val x7thm1: Double + private val deep: DeepSpaceCalculator + private val dsv = DeepSpaceValueObject() + + init { + // Recover original mean motion (xnodp) and semimajor axis (aodp) from input elements + val a1 = (xke / super.tle.xno).pow(twoThirds) + dsv.cosio = cos(super.tle.xincl) + dsv.theta2 = dsv.cosio * dsv.cosio + x3thm1 = 3.0 * dsv.theta2 - 1 + dsv.eosq = super.tle.eccn * super.tle.eccn + dsv.betao2 = 1.0 - dsv.eosq + dsv.betao = sqrt(dsv.betao2) + val del1 = 1.5 * ck2 * x3thm1 / (a1 * a1 * dsv.betao * dsv.betao2) + val ao = a1 * (1.0 - del1 * (0.5 * twoThirds + del1 * (1.0 + 134.0 / 81.0 * del1))) + val delo = 1.5 * ck2 * x3thm1 / (ao * ao * dsv.betao * dsv.betao2) + dsv.xnodp = super.tle.xno / (1.0 + delo) + dsv.aodp = ao / (1.0 - delo) + // For perigee below 156 km, the values of S and QOMS2T are altered + setPerigee((dsv.aodp * (1.0 - super.tle.eccn) - 1.0) * earthRadius) + val pinvsq = invert(dsv.aodp * dsv.aodp * dsv.betao2 * dsv.betao2) + dsv.sing = sin(super.tle.omegao) + dsv.cosg = cos(super.tle.omegao) + val tsi = invert(dsv.aodp - s4) + val eta = dsv.aodp * super.tle.eccn * tsi + val etasq = eta * eta + val eeta = super.tle.eccn * eta + val psisq = abs(1.0 - etasq) + val coef = qoms24 * tsi.pow(4.0) + val coef1 = coef / psisq.pow(3.5) + val c2 = coef1 * dsv.xnodp * (dsv.aodp * (1.0 + 1.5 * etasq + eeta * (4.0 + etasq)) + + 0.75 * ck2 * tsi / psisq * x3thm1 * (8.0 + 3.0 * etasq * (8.0 + etasq))) + c1 = super.tle.bstar * c2 + dsv.sinio = sin(super.tle.xincl) + val a3ovk2 = -j3Harmonic / ck2 + x1mth2 = 1.0 - dsv.theta2 + c4 = + 2 * dsv.xnodp * coef1 * dsv.aodp * dsv.betao2 * (eta * (2.0 + 0.5 * etasq) + super.tle.eccn + * (0.5 + 2 * etasq) - 2 * ck2 * tsi / (dsv.aodp * psisq) + * (-3 * x3thm1 * (1.0 - 2 * eeta + etasq * (1.5 - 0.5 * eeta)) + (0.75 * x1mth2 + * (2.0 * etasq - eeta * (1.0 + etasq)) * cos(2.0 * super.tle.omegao)))) + val theta4 = dsv.theta2 * dsv.theta2 + val temp1 = 3.0 * ck2 * pinvsq * dsv.xnodp + val temp2 = temp1 * ck2 * pinvsq + val temp3 = 1.25 * ck4 * pinvsq * pinvsq * dsv.xnodp + dsv.xmdot = + dsv.xnodp + 0.5 * temp1 * dsv.betao * x3thm1 + 0.0625 * temp2 * dsv.betao * (13 - 78 * dsv.theta2 + 137 * theta4) + val x1m5th = 1.0 - 5 * dsv.theta2 + dsv.omgdot = + -0.5 * temp1 * x1m5th + 0.0625 * temp2 * (7.0 - 114 * dsv.theta2 + 395 * theta4) + temp3 * (3.0 - 36 * dsv.theta2 + 49 * theta4) + val xhdot1 = -temp1 * dsv.cosio + dsv.xnodot = + xhdot1 + (0.5 * temp2 * (4.0 - 19 * dsv.theta2) + 2 * temp3 * (3.0 - 7 * dsv.theta2)) * dsv.cosio + xnodcf = 3.5 * dsv.betao2 * xhdot1 * c1 + t2cof = 1.5 * c1 + xlcof = 0.125 * a3ovk2 * dsv.sinio * (3.0 + 5 * dsv.cosio) / (1.0 + dsv.cosio) + aycof = 0.25 * a3ovk2 * dsv.sinio + x7thm1 = 7.0 * dsv.theta2 - 1 + deep = DeepSpaceCalculator(dsv) + } + + fun calculateSDP4(tSince: Double) { + synchronized(this) { + val temp = DoubleArray(12) + val xmdf = tle.xmo + dsv.xmdot * tSince + val tsq = tSince * tSince + val templ = t2cof * tsq + dsv.xll = xmdf + dsv.xnodp * templ + dsv.omgadf = tle.omegao + dsv.omgdot * tSince + val xnoddf = tle.xnodeo + dsv.xnodot * tSince + dsv.xnode = xnoddf + xnodcf * tsq + val tempa = 1.0 - c1 * tSince + val tempe = tle.bstar * c4 * tSince + dsv.xn = dsv.xnodp + dsv.t = tSince + deep.dpsec(tle) + val a = (xke / dsv.xn).pow(twoThirds) * tempa * tempa + dsv.em = dsv.em - tempe + deep.dpper() + val xl = dsv.xll + dsv.omgadf + dsv.xnode + val beta = sqrt(1.0 - dsv.em * dsv.em) + dsv.xn = xke / a.pow(1.5) + + // Long period periodics + val axn = dsv.em * cos(dsv.omgadf) + temp[0] = invert(a * beta * beta) + val xll = temp[0] * xlcof * axn + val aynl = temp[0] * aycof + val xlt = xl + xll + val ayn = dsv.em * sin(dsv.omgadf) + aynl + + // Solve Kepler's equation + val capu = mod2PI(xlt - dsv.xnode) + temp[2] = capu + converge(temp, axn, ayn, capu) + calculatePosAndVel(temp, a, axn, ayn) + } + } + + private fun calculatePosAndVel(temp: DoubleArray, a: Double, axn: Double, ayn: Double) { + val ecose = temp[5] + temp[6] + val esine = temp[3] - temp[4] + val elsq = axn * axn + ayn * ayn + temp[0] = 1.0 - elsq + val pl = a * temp[0] + temp[9] = a * (1.0 - ecose) + temp[1] = invert(temp[9]) + temp[10] = xke * sqrt(a) * esine * temp[1] + temp[11] = xke * sqrt(pl) * temp[1] + temp[2] = a * temp[1] + val betal = sqrt(temp[0]) + temp[3] = invert(1.0 + betal) + val cosu = temp[2] * (temp[8] - axn + ayn * esine * temp[3]) + val sinu = temp[2] * (temp[7] - ayn - axn * esine * temp[3]) + val u = atan2(sinu, cosu) + val sin2u = 2.0 * sinu * cosu + val cos2u = 2.0 * cosu * cosu - 1 + temp[0] = invert(pl) + temp[1] = ck2 * temp[0] + temp[2] = temp[1] * temp[0] + + // Update for short periodics + val rk = temp[9] * (1.0 - 1.5 * temp[2] * betal * x3thm1) + 0.5 * temp[1] * x1mth2 * cos2u + val uk = u - 0.25 * temp[2] * x7thm1 * sin2u + val xnodek = dsv.xnode + 1.5 * temp[2] * dsv.cosio * sin2u + val xinck = dsv.xinc + 1.5 * temp[2] * dsv.cosio * dsv.sinio * cos2u + val rdotk = temp[10] - dsv.xn * temp[1] * x1mth2 * sin2u + val rfdotk = temp[11] + dsv.xn * temp[1] * (x1mth2 * cos2u + 1.5 * x3thm1) + super.calculatePosAndVel(rk, uk, xnodek, xinck, rdotk, rfdotk) + } + + inner class DeepSpaceValueObject { + var eosq = 0.0 + var sinio = 0.0 + var cosio = 0.0 + var betao = 0.0 + var aodp = 0.0 + var theta2 = 0.0 + var sing = 0.0 + var cosg = 0.0 + var betao2 = 0.0 + var xmdot = 0.0 + var omgdot = 0.0 + var xnodot = 0.0 + var xnodp = 0.0 + + // Used by dpsec and dpper parts of Deep() + var xll = 0.0 + var omgadf = 0.0 + var xnode = 0.0 + var em = 0.0 + var xinc = 0.0 + var xn = 0.0 + var t = 0.0 + + // Used by thetg and Deep() + var ds50 = 0.0 + } + + inner class DeepSpaceCalculator(private val dsv: DeepSpaceValueObject) { + + private val zSinis = 3.9785416E-1 + private val zSings = -9.8088458E-1 + private val zNs = 1.19459E-5 + private val c1ss = 2.9864797E-6 + private val zEs = 1.675E-2 + private val zNl = 1.5835218E-4 + private val c1l = 4.7968065E-7 + private val zEl = 5.490E-2 + private val root22 = 1.7891679E-6 + private val root32 = 3.7393792E-7 + private val root44 = 7.3636953E-9 + private val root52 = 1.1428639E-7 + private val root54 = 2.1765803E-9 + private val tHdt = 4.3752691E-3 + private val q22 = 1.7891679E-6 + private val q31 = 2.1460748E-6 + private val q33 = 2.2123015E-7 + private val g22 = 5.7686396 + private val g32 = 9.5240898E-1 + private val g44 = 1.8014998 + private val g52 = 1.0508330 + private val g54 = 4.4108898 + + private val thgr: Double + private val xnq: Double + private val xqncl: Double + private val omegaq: Double + private var zmol = 0.0 + private var zmos = 0.0 + + // Many fields below cannot be final because they are iteratively refined + private var savtsn = 0.0 + private var ee2 = 0.0 + private var e3 = 0.0 + private var xi2 = 0.0 + private var xl2 = 0.0 + private var xl3 = 0.0 + private var xl4 = 0.0 + private var xgh2 = 0.0 + private var xgh3 = 0.0 + private var xgh4 = 0.0 + private var xh2 = 0.0 + private var xh3 = 0.0 + private var sse = 0.0 + private var ssi = 0.0 + private var ssg = 0.0 + private var xi3 = 0.0 + private var se2 = 0.0 + private var si2 = 0.0 + private var sl2 = 0.0 + private var sgh2 = 0.0 + private var sh2 = 0.0 + private var se3 = 0.0 + private var si3 = 0.0 + private var sl3 = 0.0 + private var sgh3 = 0.0 + private var sh3 = 0.0 + private var sl4 = 0.0 + private var sgh4 = 0.0 + private var ssl = 0.0 + private var ssh = 0.0 + private var d3210 = 0.0 + private var d3222 = 0.0 + private var d4410 = 0.0 + private var d4422 = 0.0 + private var d5220 = 0.0 + private var d5232 = 0.0 + private var d5421 = 0.0 + private var d5433 = 0.0 + private var del1 = 0.0 + private var del2 = 0.0 + private var del3 = 0.0 + private var fasx2 = 0.0 + private var fasx4 = 0.0 + private var fasx6 = 0.0 + private var xlamo = 0.0 + private val xfact: Double + private var xni: Double + private var atime: Double + private val stepp: Double + private val stepn: Double + private val step2: Double + private var preep = 0.0 + private var pl = 0.0 + private var sghs = 0.0 + private var xli: Double + private var d2201 = 0.0 + private var d2211 = 0.0 + private var sghl = 0.0 + private var sh1 = 0.0 + private var pinc = 0.0 + private var pe = 0.0 + private var shs = 0.0 + private var zsingl = 0.0 + private var zcosgl = 0.0 + private var zsinhl = 0.0 + private var zcoshl = 0.0 + private var zsinil = 0.0 + private var zcosil = 0.0 + private var a1 = 0.0 + private var a2 = 0.0 + private var a3 = 0.0 + private var a4 = 0.0 + private var a5 = 0.0 + private var a6 = 0.0 + private var a7 = 0.0 + private var a8 = 0.0 + private var a9 = 0.0 + private var a10 = 0.0 + private var ainv2 = 0.0 + private var alfdp = 0.0 + private val aqnv: Double + private var sgh = 0.0 + private var sini2 = 0.0 + private var sinis = 0.0 + private var sinok = 0.0 + private var sh = 0.0 + private var si = 0.0 + private var sil = 0.0 + private val day: Double + private var betdp = 0.0 + private var dalf = 0.0 + private var bfact = 0.0 + private var c = 0.0 + private var cc = 0.0 + private var cosis = 0.0 + private var cosok = 0.0 + private val cosq: Double + private var ctem = 0.0 + private var f322 = 0.0 + private var zx = 0.0 + private var zy = 0.0 + private var dbet = 0.0 + private var dls = 0.0 + private var eoc = 0.0 + private val eq: Double + private var f2 = 0.0 + private var f220 = 0.0 + private var f221 = 0.0 + private var f3 = 0.0 + private var f311 = 0.0 + private var f321 = 0.0 + private var xnoh = 0.0 + private var f330 = 0.0 + private var f441 = 0.0 + private var f442 = 0.0 + private var f522 = 0.0 + private var f523 = 0.0 + private var f542 = 0.0 + private var f543 = 0.0 + private var g200 = 0.0 + private var g201 = 0.0 + private var g211 = 0.0 + private var pgh = 0.0 + private var ph = 0.0 + private var s1 = 0.0 + private var s2 = 0.0 + private var s3 = 0.0 + private var s4 = 0.0 + private var s5 = 0.0 + private var s6 = 0.0 + private var s7 = 0.0 + private var se = 0.0 + private var sel = 0.0 + private var ses = 0.0 + private var xls = 0.0 + private var g300 = 0.0 + private var g310 = 0.0 + private var g322 = 0.0 + private var g410 = 0.0 + private var g422 = 0.0 + private var g520 = 0.0 + private var g521 = 0.0 + private var g532 = 0.0 + private var g533 = 0.0 + private var gam = 0.0 + private val sinq: Double + private var sinzf = 0.0 + private var sis = 0.0 + private var sl = 0.0 + private var sll = 0.0 + private var sls = 0.0 + private var stem = 0.0 + private var temp = 0.0 + private var temp1 = 0.0 + private var x1 = 0.0 + private var x2 = 0.0 + private var x2li = 0.0 + private var x2omi = 0.0 + private var x3 = 0.0 + private var x4 = 0.0 + private var x5 = 0.0 + private var x6 = 0.0 + private var x7 = 0.0 + private var x8 = 0.0 + private var xl = 0.0 + private var xldot = 0.0 + private val xmao: Double + private var xnddt = 0.0 + private var xndot = 0.0 + private var xno2 = 0.0 + private var xnodce = 0.0 + private var xnoi = 0.0 + private var xomi = 0.0 + private val xpidot: Double + private var z1 = 0.0 + private var z11 = 0.0 + private var z12 = 0.0 + private var z13 = 0.0 + private var z2 = 0.0 + private var z21 = 0.0 + private var z22 = 0.0 + private var z23 = 0.0 + private var z3 = 0.0 + private var z31 = 0.0 + private var z32 = 0.0 + private var z33 = 0.0 + private var ze = 0.0 + private var zf = 0.0 + private var zm = 0.0 + private var zn = 0.0 + private var zsing = 0.0 + private var zsinh = 0.0 + private var zsini = 0.0 + private var zcosg = 0.0 + private var zcosh = 0.0 + private var zcosi = 0.0 + private var delt = 0.0 + private var ft = 0.0 + private var resonance: Boolean + private var synchronous: Boolean + private var doLoop = false + private var epochRestart = false + + init { + thgr = thetaG(tle.epoch) + eq = tle.eccn + xnq = dsv.xnodp + aqnv = invert(dsv.aodp) + xqncl = tle.xincl + xmao = tle.xmo + xpidot = dsv.omgdot + dsv.xnodot + sinq = sin(tle.xnodeo) + cosq = cos(tle.xnodeo) + omegaq = tle.omegao + // Initialize lunar solar terms, days since 1900 Jan 0.5 + day = dsv.ds50 + 18261.5 + if (abs(day - preep) > 1.0E-6) { + preep = day + xnodce = 4.5236020 - 9.2422029E-4 * day + stem = sin(xnodce) + ctem = cos(xnodce) + zcosil = 0.91375164 - 0.03568096 * ctem + zsinil = sqrt(1.0 - zcosil * zcosil) + zsinhl = 0.089683511 * stem / zsinil + zcoshl = sqrt(1.0 - zsinhl * zsinhl) + c = 4.7199672 + 0.22997150 * day + gam = 5.8351514 + 0.0019443680 * day + zmol = mod2PI(c - gam) + zx = 0.39785416 * stem / zsinil + zy = zcoshl * ctem + 0.91744867 * zsinhl * stem + zx = atan2(zx, zy) + zx = gam + zx - xnodce + zcosgl = cos(zx) + zsingl = sin(zx) + zmos = mod2PI(6.2565837 + 0.017201977 * day) + } else { + zmol = 0.0 + zmos = 0.0 + } + doSolarTerms() + + // Geopotential resonance initialization for 12 hour orbits + resonance = false + synchronous = false + if (!(xnq < 0.0052359877 && xnq > 0.0034906585)) { + if (xnq < 0.00826 || xnq > 0.00924) + if (eq < 0.5) + // calculateResonance + resonance = true + eoc = eq * dsv.eosq + g201 = -0.306 - (eq - 0.64) * 0.440 + if (eq <= 0.65) { + g211 = 3.616 - 13.247 * eq + 16.290 * dsv.eosq + g310 = -19.302 + 117.390 * eq - 228.419 * dsv.eosq + 156.591 * eoc + g322 = -18.9068 + 109.7927 * eq - 214.6334 * dsv.eosq + 146.5816 * eoc + g410 = -41.122 + 242.694 * eq - 471.094 * dsv.eosq + 313.953 * eoc + g422 = -146.407 + 841.880 * eq - 1629.014 * dsv.eosq + 1083.435 * eoc + g520 = -532.114 + 3017.977 * eq - 5740 * dsv.eosq + 3708.276 * eoc + } else { + g211 = -72.099 + 331.819 * eq - 508.738 * dsv.eosq + 266.724 * eoc + g310 = -346.844 + 1582.851 * eq - 2415.925 * dsv.eosq + 1246.113 * eoc + g322 = -342.585 + 1554.908 * eq - 2366.899 * dsv.eosq + 1215.972 * eoc + g410 = -1052.797 + 4758.686 * eq - 7193.992 * dsv.eosq + 3651.957 * eoc + g422 = -3581.69 + 16178.11 * eq - 24462.77 * dsv.eosq + 12422.52 * eoc + g520 = + if (eq <= 0.715) 1464.74 - 4664.75 * eq + 3763.64 * dsv.eosq + else -5149.66 + 29936.92 * eq - 54087.36 * dsv.eosq + 31324.56 * eoc + } + if (eq < 0.7) { + g533 = -919.2277 + 4988.61 * eq - 9064.77 * dsv.eosq + 5542.21 * eoc + g521 = -822.71072 + 4568.6173 * eq - 8491.4146 * dsv.eosq + 5337.524 * eoc + g532 = -853.666 + 4690.25 * eq - 8624.77 * dsv.eosq + 5341.4 * eoc + } else { + g533 = -37995.78 + 161616.52 * eq - 229838.2 * dsv.eosq + 109377.94 * eoc + g521 = -51752.104 + 218913.95 * eq - 309468.16 * dsv.eosq + 146349.42 * eoc + g532 = -40023.88 + 170470.89 * eq - 242699.48 * dsv.eosq + 115605.82 * eoc + } + sini2 = dsv.sinio * dsv.sinio + f220 = 0.75 * (1.0 + 2 * dsv.cosio + dsv.theta2) + f221 = 1.5 * sini2 + f321 = 1.875 * dsv.sinio * (1.0 - 2 * dsv.cosio - 3.0 * dsv.theta2) + f322 = -1.875 * dsv.sinio * (1.0 + 2 * dsv.cosio - 3.0 * dsv.theta2) + f441 = 35 * sini2 * f220 + f442 = 39.3750 * sini2 * sini2 + f522 = + 9.84375 * dsv.sinio * (sini2 * (1.0 - 2 * dsv.cosio - 5 * dsv.theta2) + 0.33333333 * (-2 + 4 * dsv.cosio + 6 * dsv.theta2)) + f523 = + dsv.sinio * (4.92187512 * sini2 * (-2 - 4 * dsv.cosio + 10 * dsv.theta2) + 6.56250012 * (1.0 + 2 * dsv.cosio - 3.0 * dsv.theta2)) + f542 = + 29.53125 * dsv.sinio * (2.0 - 8 * dsv.cosio + dsv.theta2 * (-12 + 8 * dsv.cosio + 10 * dsv.theta2)) + f543 = + 29.53125 * dsv.sinio * (-2 - 8 * dsv.cosio + dsv.theta2 * (12 + 8 * dsv.cosio - 10 * dsv.theta2)) + xno2 = xnq * xnq + ainv2 = aqnv * aqnv + temp1 = 3.0 * xno2 * ainv2 + temp = temp1 * root22 + d2201 = temp * f220 * g201 + d2211 = temp * f221 * g211 + temp1 *= aqnv + temp = temp1 * root32 + d3210 = temp * f321 * g310 + d3222 = temp * f322 * g322 + temp1 *= aqnv + temp = 2.0 * temp1 * root44 + d4410 = temp * f441 * g410 + d4422 = temp * f442 * g422 + temp1 *= aqnv + temp = temp1 * root52 + d5220 = temp * f522 * g520 + d5232 = temp * f523 * g532 + temp = 2.0 * temp1 * root54 + d5421 = temp * f542 * g521 + d5433 = temp * f543 * g533 + xlamo = xmao + tle.xnodeo + tle.xnodeo - thgr - thgr + bfact = dsv.xmdot + dsv.xnodot + dsv.xnodot - tHdt - tHdt + bfact += ssl + ssh + ssh + } else { + // Init synchronous resonance terms + resonance = true + synchronous = true + g200 = 1.0 + dsv.eosq * (-2.5 + 0.8125 * dsv.eosq) + g310 = 1.0 + 2 * dsv.eosq + g300 = 1.0 + dsv.eosq * (-6 + 6.60937 * dsv.eosq) + f220 = 0.75 * (1.0 + dsv.cosio) * (1.0 + dsv.cosio) + f311 = + 0.9375 * dsv.sinio * dsv.sinio * (1.0 + 3.0 * dsv.cosio) - 0.75 * (1.0 + dsv.cosio) + f330 = 1.0 + dsv.cosio + f330 *= 1.875 * f330 * f330 + del1 = 3.0 * xnq * xnq * aqnv * aqnv + del2 = 2.0 * del1 * f220 * g200 * q22 + del3 = 3.0 * del1 * f330 * g300 * q33 * aqnv + del1 *= f311 * g310 * q31 * aqnv + fasx2 = 0.13130908 + fasx4 = 2.8843198 + fasx6 = 0.37448087 + xlamo = xmao + tle.xnodeo + tle.omegao - thgr + bfact = dsv.xmdot + xpidot - tHdt + bfact += ssl + ssg + ssh + } + xfact = bfact - xnq + // Init integrator + xli = xlamo + xni = xnq + atime = 0.0 + stepp = 720.0 + stepn = -720.0 + step2 = 259200.0 + } + + // Entrance for lunar-solar periodics + fun dpper() { + sinis = sin(dsv.xinc) + cosis = cos(dsv.xinc) + if (abs(savtsn - dsv.t) >= 30) { + savtsn = dsv.t + zm = zmos + zNs * dsv.t + zf = zm + 2 * zEs * sin(zm) + sinzf = sin(zf) + f2 = 0.5 * sinzf * sinzf - 0.25 + f3 = -0.5 * sinzf * cos(zf) + ses = se2 * f2 + se3 * f3 + sis = si2 * f2 + si3 * f3 + sls = sl2 * f2 + sl3 * f3 + sl4 * sinzf + sghs = sgh2 * f2 + sgh3 * f3 + sgh4 * sinzf + shs = sh2 * f2 + sh3 * f3 + zm = zmol + zNl * dsv.t + zf = zm + 2 * zEl * sin(zm) + sinzf = sin(zf) + f2 = 0.5 * sinzf * sinzf - 0.25 + f3 = -0.5 * sinzf * cos(zf) + sel = ee2 * f2 + e3 * f3 + sil = xi2 * f2 + xi3 * f3 + sll = xl2 * f2 + xl3 * f3 + xl4 * sinzf + sghl = xgh2 * f2 + xgh3 * f3 + xgh4 * sinzf + sh1 = xh2 * f2 + xh3 * f3 + pe = ses + sel + pinc = sis + sil + pl = sls + sll + } + pgh = sghs + sghl + ph = shs + sh1 + dsv.xinc = dsv.xinc + pinc + dsv.em = dsv.em + pe + if (xqncl >= 0.2) { + /* Apply periodics directly */ + ph /= dsv.sinio + pgh -= dsv.cosio * ph + dsv.omgadf = dsv.omgadf + pgh + dsv.xnode = dsv.xnode + ph + dsv.xll = dsv.xll + pl + } else { + applyPeriodics() + // This is a patch to Lyddane modification suggested by Rob Matson + if (abs(xnoh - dsv.xnode) > Math.PI) { + if (dsv.xnode < xnoh) dsv.xnode += twoPi else dsv.xnode -= twoPi + } + dsv.xll = dsv.xll + pl + dsv.omgadf = xls - dsv.xll - cos(dsv.xinc) * dsv.xnode + } + } + + // Entrance for deep space secular effects + fun dpsec(tle: TLE) { + dsv.xll = dsv.xll + ssl * dsv.t + dsv.omgadf = dsv.omgadf + ssg * dsv.t + dsv.xnode = dsv.xnode + ssh * dsv.t + dsv.em = tle.eccn + sse * dsv.t + dsv.xinc = tle.xincl + ssi * dsv.t + if (dsv.xinc < 0) { + dsv.xinc = -dsv.xinc + dsv.xnode = dsv.xnode + Math.PI + dsv.omgadf = dsv.omgadf - Math.PI + } + if (!resonance) return + do processEpochRestartLoop() while (doLoop && epochRestart) + dsv.xn = xni + xndot * ft + xnddt * ft * ft * 0.5 + xl = xli + xldot * ft + xndot * ft * ft * 0.5 + temp = -dsv.xnode + thgr + dsv.t * tHdt + if (synchronous) dsv.xll = xl - dsv.omgadf + temp else dsv.xll = xl + temp + temp + } + + private fun doSolarTerms() { + savtsn = 1E20 + zcosg = 1.945905E-1 + zsing = zSings + zcosi = 9.1744867E-1 + zsini = zSinis + zcosh = cosq + zsinh = sinq + cc = c1ss + zn = zNs + ze = zEs + xnoi = invert(xnq) + calculateSolarTerms() + calculateLunarTerms() + calculateSolarTerms() // Solar terms done again after Lunar terms are done + sse += se + ssi += si + ssl += sl + ssg = ssg + sgh - dsv.cosio / dsv.sinio * sh + ssh += sh / dsv.sinio + } + + private fun calculateLunarTerms() { + sse = se + ssi = si + ssl = sl + ssh = sh / dsv.sinio + ssg = sgh - dsv.cosio * ssh + se2 = ee2 + si2 = xi2 + sl2 = xl2 + sgh2 = xgh2 + sh2 = xh2 + se3 = e3 + si3 = xi3 + sl3 = xl3 + sgh3 = xgh3 + sh3 = xh3 + sl4 = xl4 + sgh4 = xgh4 + zcosg = zcosgl + zsing = zsingl + zcosi = zcosil + zsini = zsinil + zcosh = zcoshl * cosq + zsinhl * sinq + zsinh = sinq * zcoshl - cosq * zsinhl + zn = zNl + cc = c1l + ze = zEl + } + + private fun calculateSolarTerms() { + a1 = zcosg * zcosh + zsing * zcosi * zsinh + a3 = -zsing * zcosh + zcosg * zcosi * zsinh + a7 = -zcosg * zsinh + zsing * zcosi * zcosh + a8 = zsing * zsini + a9 = zsing * zsinh + zcosg * zcosi * zcosh + a10 = zcosg * zsini + a2 = dsv.cosio * a7 + dsv.sinio * a8 + a4 = dsv.cosio * a9 + dsv.sinio * a10 + a5 = -dsv.sinio * a7 + dsv.cosio * a8 + a6 = -dsv.sinio * a9 + dsv.cosio * a10 + x1 = a1 * dsv.cosg + a2 * dsv.sing + x2 = a3 * dsv.cosg + a4 * dsv.sing + x3 = -a1 * dsv.sing + a2 * dsv.cosg + x4 = -a3 * dsv.sing + a4 * dsv.cosg + x5 = a5 * dsv.sing + x6 = a6 * dsv.sing + x7 = a5 * dsv.cosg + x8 = a6 * dsv.cosg + z31 = 12 * x1 * x1 - 3.0 * x3 * x3 + z32 = 24 * x1 * x2 - 6 * x3 * x4 + z33 = 12 * x2 * x2 - 3.0 * x4 * x4 + z1 = 3.0 * (a1 * a1 + a2 * a2) + z31 * dsv.eosq + z2 = 6.0 * (a1 * a3 + a2 * a4) + z32 * dsv.eosq + z3 = 3.0 * (a3 * a3 + a4 * a4) + z33 * dsv.eosq + z11 = -6 * a1 * a5 + dsv.eosq * (-24 * x1 * x7 - 6 * x3 * x5) + z12 = + -6 * (a1 * a6 + a3 * a5) + dsv.eosq * (-24 * (x2 * x7 + x1 * x8) - 6 * (x3 * x6 + x4 * x5)) + z13 = -6 * a3 * a6 + dsv.eosq * (-24 * x2 * x8 - 6 * x4 * x6) + z21 = 6.0 * a2 * a5 + dsv.eosq * (24 * x1 * x5 - 6 * x3 * x7) + z22 = + 6.0 * (a4 * a5 + a2 * a6) + dsv.eosq * (24 * (x2 * x5 + x1 * x6) - 6 * (x4 * x7 + x3 * x8)) + z23 = 6.0 * a4 * a6 + dsv.eosq * (24 * x2 * x6 - 6 * x4 * x8) + z1 += z1 + dsv.betao2 * z31 + z2 += z2 + dsv.betao2 * z32 + z3 += z3 + dsv.betao2 * z33 + s3 = cc * xnoi + s2 = -0.5 * s3 / dsv.betao + s4 = s3 * dsv.betao + s1 = -15 * eq * s4 + s5 = x1 * x3 + x2 * x4 + s6 = x2 * x3 + x1 * x4 + s7 = x2 * x4 - x1 * x3 + se = s1 * zn * s5 + si = s2 * zn * (z11 + z13) + sl = -zn * s3 * (z1 + z3 - 14 - 6 * dsv.eosq) + sgh = s4 * zn * (z31 + z33 - 6) + sh = -zn * s2 * (z21 + z23) + if (xqncl < 5.2359877E-2) sh = 0.0 + ee2 = 2.0 * s1 * s6 + e3 = 2.0 * s1 * s7 + xi2 = 2.0 * s2 * z12 + xi3 = 2.0 * s2 * (z13 - z11) + xl2 = -2 * s3 * z2 + xl3 = -2 * s3 * (z3 - z1) + xl4 = -2 * s3 * (-21 - 9 * dsv.eosq) * ze + xgh2 = 2.0 * s4 * z32 + xgh3 = 2.0 * s4 * (z33 - z31) + xgh4 = -18 * s4 * ze + xh2 = -2 * s2 * z22 + xh3 = -2 * s2 * (z23 - z21) + } + + private fun processEpochRestartLoop() { + if (atime == 0.0 || dsv.t >= 0 && atime < 0 || dsv.t < 0 && atime >= 0) { + calculateDelta() + atime = 0.0 + xni = xnq + xli = xlamo + } else if (abs(dsv.t) >= abs(atime)) calculateDelta() + processNotEpochRestartLoop() + } + + private fun calculateDelta() { + delt = if (dsv.t < 0) stepn else stepp + } + + private fun processNotEpochRestartLoop() { + do { + if (abs(dsv.t - atime) >= stepp) { + doLoop = true + epochRestart = false + } else { + ft = dsv.t - atime + doLoop = false + } + if (abs(dsv.t) < abs(atime)) { + delt = if (dsv.t >= 0) stepn else stepp + doLoop = doLoop or epochRestart + } + if (synchronous) { + xndot = del1 * sin(xli - fasx2) + del2 * sin(2.0 * (xli - fasx4)) + +del3 * sin(3.0 * (xli - fasx6)) + xnddt = del1 * cos(xli - fasx2) + 2 * del2 * cos(2.0 * (xli - fasx4)) + +3.0 * del3 * cos(3.0 * (xli - fasx6)) + } else { + xomi = omegaq + dsv.omgdot * atime + x2omi = xomi + xomi + x2li = xli + xli + xndot = + d2201 * sin(x2omi + xli - g22) + d2211 * sin(xli - g22) + (d3210 + * sin(xomi + xli - g32)) + d3222 * sin(-xomi + xli - g32) + (d4410 + * sin(x2omi + x2li - g44)) + d4422 * sin(x2li - g44) + (d5220 + * sin(xomi + xli - g52)) + d5232 * sin(-xomi + xli - g52) + (d5421 + * sin(xomi + x2li - g54)) + d5433 * sin(-xomi + x2li - g54) + xnddt = + d2201 * cos(x2omi + xli - g22) + d2211 * cos(xli - g22) + (d3210 + * cos(xomi + xli - g32)) + d3222 * cos(-xomi + xli - g32) + (d5220 + * cos(xomi + xli - g52)) + d5232 * cos(-xomi + xli - g52) + (2 + * (d4410 * cos(x2omi + x2li - g44) + d4422 * cos(x2li - g44) + (d5421 + * cos(xomi + x2li - g54)) + d5433 * cos(-xomi + x2li - g54))) + } + xldot = xni + xfact + xnddt *= xldot + if (doLoop) { + xli += xldot * delt + xndot * step2 + xni += xndot * delt + xnddt * step2 + atime += delt + } + } while (doLoop && !epochRestart) + } + + // Apply periodics with Lyddane modification + private fun applyPeriodics() { + sinok = sin(dsv.xnode) + cosok = cos(dsv.xnode) + alfdp = sinis * sinok + betdp = sinis * cosok + dalf = ph * cosok + pinc * cosis * sinok + dbet = -ph * sinok + pinc * cosis * cosok + alfdp += dalf + betdp += dbet + dsv.xnode = mod2PI(dsv.xnode) + xls = dsv.xll + dsv.omgadf + cosis * dsv.xnode + dls = pl + pgh - pinc * dsv.xnode * sinis + xls += dls + xnoh = dsv.xnode + dsv.xnode = atan2(alfdp, betdp) + } + + // Calculates the Greenwich Mean Sidereal Time for an epoch, valid 1957 through 2056 + private fun thetaG(epoch: Double): Double { + var year = floor(epoch * 1E-3) + var dayOfYear = (epoch * 1E-3 - year) * 1000.0 + year = if (year < 57) year + 2000 else year + 1900 + val dayFloor = floor(dayOfYear) + val dayFraction = dayOfYear - dayFloor + dayOfYear = dayFloor + val jd = julianDateOfYear(year) + dayOfYear + dsv.ds50 = jd - 2433281.5 + dayFraction + return mod2PI(6.3003880987 * dsv.ds50 + 1.72944494) + } + } +} diff --git a/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/GroundPos.kt b/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/GroundPos.kt new file mode 100644 index 00000000..44cfeb51 --- /dev/null +++ b/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/GroundPos.kt @@ -0,0 +1,20 @@ +/* + * Look4Sat. Amateur radio satellite tracker and pass predictor. + * Copyright (C) 2019-2021 Arty Bishop (bishop.arty@gmail.com) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ +package com.rtbishop.look4sat.predict4kotlin + +data class GroundPos(val lat: Double, val lon: Double, val alt: Double, val name: String = "base") diff --git a/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/NearEarthSat.kt b/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/NearEarthSat.kt new file mode 100644 index 00000000..3bc53450 --- /dev/null +++ b/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/NearEarthSat.kt @@ -0,0 +1,221 @@ +/* + * Look4Sat. Amateur radio satellite tracker and pass predictor. + * Copyright (C) 2019-2021 Arty Bishop (bishop.arty@gmail.com) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ +package com.rtbishop.look4sat.predict4kotlin + +import kotlin.math.* + +class NearEarthSat(tle: TLE) : Satellite(tle) { + + private val aodp: Double + private val aycof: Double + private val c1: Double + private val c4: Double + private val c5: Double + private val cosio: Double + private var d2 = 0.0 + private var d3 = 0.0 + private var d4 = 0.0 + private val delmo: Double + private val omgcof: Double + private val eta: Double + private val omgdot: Double + private val sinio: Double + private val xnodp: Double + private val sinmo: Double + private val t2cof: Double + private var t3cof = 0.0 + private var t4cof = 0.0 + private var t5cof = 0.0 + private val x1mth2: Double + private val x3thm1: Double + private val x7thm1: Double + private val xmcof: Double + private val xmdot: Double + private val xnodcf: Double + private val xnodot: Double + private val xlcof: Double + private val sgp4Simple: Boolean + + init { + // Recover original mean motion (xnodp) and semimajor axis (aodp) from input elements + val a1 = (xke / super.tle.xno).pow(twoThirds) + cosio = cos(super.tle.xincl) + val theta2 = sqr(cosio) + x3thm1 = 3.0 * theta2 - 1.0 + val eo = super.tle.eccn + val eosq = sqr(eo) + val betao2 = 1.0 - eosq + val betao = sqrt(betao2) + val del1 = 1.5 * ck2 * x3thm1 / (sqr(a1) * betao * betao2) + val ao = a1 * (1.0 - del1 * (0.5 * twoThirds + del1 * (1.0 + 134.0 / 81.0 * del1))) + val delo = 1.5 * ck2 * x3thm1 / (sqr(ao) * betao * betao2) + xnodp = super.tle.xno / (1.0 + delo) + aodp = ao / (1.0 - delo) + + // For perigee less than 220 kilometers, the "simple" flag is set + sgp4Simple = aodp * (1.0 - eo) < 220 / earthRadius + 1.0 + + // For perigees below 156 km, the values of S and QOMS2T are altered + setPerigee((aodp * (1.0 - eo) - 1.0) * earthRadius) + val pinvsq = invert(sqr(aodp) * sqr(betao2)) + val tsi = invert(aodp - s4) + eta = aodp * eo * tsi + val etasq = eta * eta + val eeta = eo * eta + val psisq = abs(1.0 - etasq) + val coef = qoms24 * tsi.pow(4.0) + val coef1 = coef / psisq.pow(3.5) + val bstar = super.tle.bstar + val c2 = coef1 * xnodp * (aodp * (1.0 + 1.5 * etasq + eeta * (4.0 + etasq)) + 0.75 + * ck2 * tsi / psisq * x3thm1 * (8.0 + 3.0 * etasq * (8.0 + etasq))) + c1 = bstar * c2 + sinio = sin(super.tle.xincl) + val a3ovk2 = -j3Harmonic / ck2 + val c3 = coef * tsi * a3ovk2 * xnodp * sinio / eo + x1mth2 = 1.0 - theta2 + val omegao = super.tle.omegao + c4 = 2 * xnodp * coef1 * aodp * betao2 * (eta * (2.0 + 0.5 * etasq) + eo * (0.5 + 2 * etasq) + - 2 * ck2 * tsi / (aodp * psisq) * (-3 * x3thm1 * (1.0 - 2 * eeta + etasq + * (1.5 - 0.5 * eeta)) + 0.75 * x1mth2 * (2.0 * etasq - eeta * (1.0 + etasq)) + * cos(2.0 * omegao))) + c5 = 2.0 * coef1 * aodp * betao2 * (1.0 + 2.75 * (etasq + eeta) + eeta * etasq) + val theta4 = sqr(theta2) + val temp1 = 3.0 * ck2 * pinvsq * xnodp + val temp2 = temp1 * ck2 * pinvsq + val temp3 = 1.25 * ck4 * pinvsq * pinvsq * xnodp + xmdot = + xnodp + 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 + omgcof = bstar * c3 * cos(omegao) + xmcof = -twoThirds * coef * bstar / eeta + xnodcf = 3.5 * betao2 * xhdot1 * c1 + t2cof = 1.5 * c1 + xlcof = 0.125 * a3ovk2 * sinio * (3.0 + 5 * cosio) / (1.0 + cosio) + aycof = 0.25 * a3ovk2 * sinio + val xmo = super.tle.xmo + delmo = (1.0 + eta * cos(xmo)).pow(3.0) + sinmo = sin(xmo) + x7thm1 = 7.0 * theta2 - 1 + if (!sgp4Simple) { + val c1sq = sqr(c1) + d2 = 4.0 * aodp * tsi * c1sq + val temp = d2 * tsi * c1 / 3.0 + d3 = (17 * aodp + s4) * temp + d4 = 0.5 * temp * aodp * tsi * (221 * aodp + 31 * s4) * c1 + t3cof = d2 + 2 * c1sq + t4cof = 0.25 * (3.0 * d3 + c1 * (12 * d2 + 10 * c1sq)) + t5cof = 0.2 * (3.0 * d4 + 12 * c1 * d3 + 6 * d2 * d2 + 15 * c1sq * (2.0 * d2 + c1sq)) + } else { + d2 = 0.0 + d3 = 0.0 + d4 = 0.0 + t3cof = 0.0 + t4cof = 0.0 + t5cof = 0.0 + } + } + + fun calculateSGP4(tSince: Double) { + synchronized(this) { + val temp = DoubleArray(9) + val xmdf = tle.xmo + xmdot * tSince + val omgadf = tle.omegao + omgdot * tSince + val xnoddf = tle.xnodeo + xnodot * tSince + var omega = omgadf + var xmp = xmdf + val tsq = sqr(tSince) + val xnode = xnoddf + xnodcf * tsq + val bstar = tle.bstar + var tempa = 1.0 - c1 * tSince + var tempe = bstar * c4 * tSince + var templ = t2cof * tsq + if (!sgp4Simple) { + val delomg = omgcof * tSince + val delm = xmcof * ((1.0 + eta * cos(xmdf)).pow(3.0) - delmo) + temp[0] = delomg + delm + xmp = xmdf + temp[0] + omega = omgadf - temp[0] + val tcube = tsq * tSince + val tfour = tSince * tcube + tempa = tempa - d2 * tsq - d3 * tcube - d4 * tfour + tempe += bstar * c5 * (sin(xmp) - sinmo) + templ += t3cof * tcube + tfour * (t4cof + tSince * t5cof) + } + val a = aodp * tempa.pow(2.0) + val eo = tle.eccn + val e = eo - tempe + val xl = xmp + omega + xnode + xnodp * templ + val beta = sqrt(1.0 - e * e) + val xn = xke / a.pow(1.5) + + // Long period periodics + val axn = e * cos(omega) + temp[0] = invert(a * sqr(beta)) + val xll = temp[0] * xlcof * axn + val aynl = temp[0] * aycof + val xlt = xl + xll + val ayn = e * sin(omega) + aynl + + // Solve Kepler's equation + val capu = mod2PI(xlt - xnode) + temp[2] = capu + converge(temp, axn, ayn, capu) + calculatePosAndVel(temp, xnode, a, xn, axn, ayn) + } + } + + private fun calculatePosAndVel( + temp: DoubleArray, xnode: Double, a: Double, + xn: Double, axn: Double, ayn: Double + ) { + val ecose = temp[5] + temp[6] + val esine = temp[3] - temp[4] + val elsq = sqr(axn) + sqr(ayn) + temp[0] = 1.0 - elsq + val pl = a * temp[0] + val r = a * (1.0 - ecose) + temp[1] = invert(r) + val rdot = xke * sqrt(a) * esine * temp[1] + val rfdot = xke * sqrt(pl) * temp[1] + temp[2] = a * temp[1] + val betal = sqrt(temp[0]) + temp[3] = invert(1.0 + betal) + val cosu = temp[2] * (temp[8] - axn + ayn * esine * temp[3]) + val sinu = temp[2] * (temp[7] - ayn - axn * esine * temp[3]) + val u = atan2(sinu, cosu) + val sin2u = 2.0 * sinu * cosu + val cos2u = 2.0 * cosu * cosu - 1 + temp[0] = invert(pl) + temp[1] = ck2 * temp[0] + temp[2] = temp[1] * temp[0] + + // Update for short periodics + val rk = r * (1.0 - 1.5 * temp[2] * betal * x3thm1) + 0.5 * temp[1] * x1mth2 * cos2u + val uk = u - 0.25 * temp[2] * x7thm1 * sin2u + val xnodek = xnode + 1.5 * temp[2] * cosio * sin2u + val xinck = tle.xincl + 1.5 * temp[2] * cosio * sinio * cos2u + val rdotk = rdot - xn * temp[1] * x1mth2 * sin2u + val rfdotk = rfdot + xn * temp[1] * (x1mth2 * cos2u + 1.5 * x3thm1) + super.calculatePosAndVel(rk, uk, xnodek, xinck, rdotk, rfdotk) + } +} diff --git a/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/PassPredictor.kt b/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/PassPredictor.kt new file mode 100644 index 00000000..423a09ae --- /dev/null +++ b/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/PassPredictor.kt @@ -0,0 +1,180 @@ +/* + * Look4Sat. Amateur radio satellite tracker and pass predictor. + * Copyright (C) 2019-2021 Arty Bishop (bishop.arty@gmail.com) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ +package com.rtbishop.look4sat.predict4kotlin + +import java.util.* + +class PassPredictor(private val satellite: Satellite, private val qth: GroundPos) { + + private val oneQuarterOrbitMin = (24.0 * 60.0 / satellite.tle.meanmo / 4.0).toInt() + private val speedOfLight = 2.99792458E8 + + fun getDownlinkFreq(freq: Long, date: Date): Long { + val rangeRate = getSatPos(date).rangeRate + return (freq.toDouble() * (speedOfLight - rangeRate * 1000.0) / speedOfLight).toLong() + } + + fun getUplinkFreq(freq: Long, date: Date): Long { + val rangeRate = getSatPos(date).rangeRate + return (freq.toDouble() * (speedOfLight + rangeRate * 1000.0) / speedOfLight).toLong() + } + + fun getSatPos(date: Date): SatPos { + return satellite.getPosition(qth, date) + } + + fun getPositions(refDate: Date, stepSec: Int, minBefore: Int, orbits: Double): List { + val positions = mutableListOf() + val orbitalPeriod = 24 * 60 / satellite.tle.meanmo + val endDate = Date(refDate.time + (orbitalPeriod * orbits * 60L * 1000L).toLong()) + val startDate = Date(refDate.time - minBefore * 60L * 1000L) + var currentDate = startDate + while (currentDate.before(endDate)) { + positions.add(getSatPos(currentDate)) + currentDate = Date(currentDate.time + stepSec * 1000) + } + return positions + } + + fun getPasses(refDate: Date, hoursAhead: Int, windBack: Boolean): List { + val passes = mutableListOf() + val endDate = Date(refDate.time + hoursAhead * 60L * 60L * 1000L) + var startDate = refDate + var shouldWindBack = windBack + var lastAosDate: Date + var count = 0 + if (satellite.willBeSeen(qth)) { + if (satellite.tle.isDeepspace) { + passes.add(nextDeepSpacePass(refDate)) + } else { + do { + if (count > 0) shouldWindBack = false + val pass = nextNearEarthPass(startDate, shouldWindBack) + lastAosDate = pass.aosDate + passes.add(pass) + startDate = + Date(pass.losDate.time + (oneQuarterOrbitMin * 3) * 60L * 1000L) + count++ + } while (lastAosDate < endDate) + } + } + return passes + } + + private fun nextDeepSpacePass(refDate: Date): SatPass { + val satPos = getSatPos(refDate) + val id = satellite.tle.catnum + val name = satellite.tle.name + val isDeep = satellite.tle.isDeepspace + val aos = Date(refDate.time - 24 * 60L * 60L * 1000L).time + val los = Date(refDate.time + 24 * 60L * 60L * 1000L).time + val tca = Date((aos + los) / 2).time + val az = Math.toDegrees(satPos.azimuth) + val elev = Math.toDegrees(satPos.elevation) + val alt = satPos.altitude + return SatPass(id, name, isDeep, aos, az, los, az, tca, az, alt, elev, this) + } + + private fun nextNearEarthPass(refDate: Date, windBack: Boolean = false): SatPass { + val calendar = Calendar.getInstance(TimeZone.getTimeZone("UTC")).apply { + clear() + timeInMillis = refDate.time + } + val id = satellite.tle.catnum + val name = satellite.tle.name + val isDeep = satellite.tle.isDeepspace + + var elevation: Double + var maxElevation = 0.0 + var alt = 0.0 + var tcaAz = 0.0 + + // wind back time 1/4 of an orbit + if (windBack) calendar.add(Calendar.MINUTE, -oneQuarterOrbitMin) + var satPos = getSatPos(calendar.time) + + if (satPos.elevation > 0.0) { + // move forward in 30 second intervals until the sat goes below the horizon + do { + calendar.add(Calendar.SECOND, 30) + satPos = getSatPos(calendar.time) + } while (satPos.elevation > 0.0) + // move forward 3/4 of an orbit + calendar.add(Calendar.MINUTE, oneQuarterOrbitMin * 3) + } + + // find the next time sat comes above the horizon + do { + calendar.add(Calendar.SECOND, 60) + satPos = getSatPos(calendar.time) + elevation = satPos.elevation + if (elevation > maxElevation) { + maxElevation = elevation + alt = satPos.altitude + tcaAz = Math.toDegrees(satPos.azimuth) + } + } while (satPos.elevation < 0.0) + + // refine to 3 seconds + calendar.add(Calendar.SECOND, -60) + do { + calendar.add(Calendar.SECOND, 3) + satPos = getSatPos(calendar.time) + elevation = satPos.elevation + if (elevation > maxElevation) { + maxElevation = elevation + alt = satPos.altitude + tcaAz = Math.toDegrees(satPos.azimuth) + } + } while (satPos.elevation < 0.0) + + val aos = satPos.time.time + val aosAz = Math.toDegrees(satPos.azimuth) + + // find when sat goes below + do { + calendar.add(Calendar.SECOND, 30) + satPos = getSatPos(calendar.time) + elevation = satPos.elevation + if (elevation > maxElevation) { + maxElevation = elevation + alt = satPos.altitude + tcaAz = Math.toDegrees(satPos.azimuth) + } + } while (satPos.elevation > 0.0) + + // refine to 3 seconds + calendar.add(Calendar.SECOND, -30) + do { + calendar.add(Calendar.SECOND, 3) + satPos = getSatPos(calendar.time) + elevation = satPos.elevation + if (elevation > maxElevation) { + maxElevation = elevation + alt = satPos.altitude + tcaAz = Math.toDegrees(satPos.azimuth) + } + } while (satPos.elevation > 0.0) + + val los = satPos.time.time + val losAz = Math.toDegrees(satPos.azimuth) + val tca = Date((aos + los) / 2).time + val elev = Math.toDegrees(maxElevation) + return SatPass(id, name, isDeep, aos, aosAz, los, losAz, tca, tcaAz, alt, elev, this) + } +} diff --git a/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/Position.kt b/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/Position.kt new file mode 100644 index 00000000..3f5af4a7 --- /dev/null +++ b/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/Position.kt @@ -0,0 +1,20 @@ +/* + * Look4Sat. Amateur radio satellite tracker and pass predictor. + * Copyright (C) 2019-2021 Arty Bishop (bishop.arty@gmail.com) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ +package com.rtbishop.look4sat.predict4kotlin + +data class Position(val lat: Double, val lon: Double) diff --git a/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/SatPass.kt b/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/SatPass.kt new file mode 100644 index 00000000..61ea418a --- /dev/null +++ b/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/SatPass.kt @@ -0,0 +1,45 @@ +/* + * Look4Sat. Amateur radio satellite tracker and pass predictor. + * Copyright (C) 2019-2021 Arty Bishop (bishop.arty@gmail.com) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ +package com.rtbishop.look4sat.predict4kotlin + +import java.util.* + +data class SatPass( + val catNum: Int, + val name: String, + val isDeepSpace: Boolean, + private val aosTime: Long, + val aosAzimuth: Double, + private val losTime: Long, + val losAzimuth: Double, + private val tcaTime: Long, + val tcaAzimuth: Double, + val altitude: Double, + val maxElevation: Double, + val predictor: PassPredictor, + var progress: Int = 0 +) { + val aosDate: Date + get() = Date(aosTime) + + val losDate: Date + get() = Date(losTime) + + val tcaDate: Date + get() = Date(tcaTime) +} diff --git a/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/SatPos.kt b/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/SatPos.kt new file mode 100644 index 00000000..00bba065 --- /dev/null +++ b/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/SatPos.kt @@ -0,0 +1,100 @@ +/* + * Look4Sat. Amateur radio satellite tracker and pass predictor. + * Copyright (C) 2019-2021 Arty Bishop (bishop.arty@gmail.com) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ +package com.rtbishop.look4sat.predict4kotlin + +import java.util.* +import kotlin.math.* + +class SatPos { + + private val earthRadiusKm = 6.378137E3 + private val r0 = 6378.16 + + // Radians + var azimuth = 0.0 + var elevation = 0.0 + var latitude = 0.0 + var longitude = 0.0 + + var time = Date() + var range = 0.0 + var rangeRate = 0.0 + var phase = 0.0 + var altitude = 0.0 + var theta = 0.0 + var eclipseDepth = 0.0 + var eclipsed = false + var aboveHorizon = false + + fun getDate(): Date { + return Date(time.time) + } + + fun getRangeCircleRadiusKm(): Double { + return 0.5 * (12756.33 * acos(earthRadiusKm / (earthRadiusKm + altitude))) + } + + fun getRangeCircle(incrementDegrees: Double = 1.0): List { + val positions = mutableListOf() + val radiusKm = this.getRangeCircleRadiusKm() + val lat = this.latitude + val lon = this.longitude + val beta = radiusKm / r0 + var tempAzimuth = 0 + while (tempAzimuth < 360) { + val azimuth = tempAzimuth / 360.0 * 2.0 * Math.PI + var rangelat = asin(sin(lat) * cos(beta) + cos(azimuth) * sin(beta) * cos(lat)) + val num = (cos(beta) - (sin(lat) * sin(rangelat))) + val den = cos(lat) * cos(rangelat) + var rangelon = if (tempAzimuth == 0 && (beta > ((Math.PI / 2.0) - lat))) { + lon + Math.PI + } else if (tempAzimuth == 180 && (beta > ((Math.PI / 2.0) - lat))) { + lon + Math.PI + } else if (abs(num / den) > 1.0) { + lon + } else { + if ((180 - tempAzimuth) >= 0) { + lon - acos(num / den) + } else { + lon + acos(num / den) + } + } + while (rangelon < 0.0) rangelon += Math.PI * 2.0 + while (rangelon > Math.PI * 2.0) rangelon -= Math.PI * 2.0 + rangelat = (rangelat / (2.0 * Math.PI)) * 360.0 + rangelon = (rangelon / (2.0 * Math.PI)) * 360.0 + + // if (rangelong < 180.0) { + // rangelong = -rangelong; + // } + // else if (rangelong > 180.0) { + // rangelong = 360.0 - rangelong; + // } + // + // if (rangelat < 90.0) { + // rangelat = -rangelat; + // } + // else if (rangelat > 90.0) { + // rangelat = 180.0 - rangelat; + // } + positions.add(Position(rangelat, rangelon)) + tempAzimuth += incrementDegrees.toInt() + } + return positions + } +} diff --git a/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/Satellite.kt b/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/Satellite.kt new file mode 100644 index 00000000..5698e39e --- /dev/null +++ b/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/Satellite.kt @@ -0,0 +1,389 @@ +/* + * Look4Sat. Amateur radio satellite tracker and pass predictor. + * Copyright (C) 2019-2021 Arty Bishop (bishop.arty@gmail.com) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ +package com.rtbishop.look4sat.predict4kotlin + +import java.io.InputStream +import java.util.* +import java.util.concurrent.atomic.AtomicReference +import kotlin.math.* + +abstract class Satellite(val tle: TLE) { + + private val flatFactor = 3.35281066474748E-3 + private val deg2Rad = 1.745329251994330E-2 + private val secPerDay = 8.6400E4 + private val minPerDay = 1.44E3 + private val epsilon = 1.0E-12 + private val position = Vector4() + private val velocity = Vector4() + private var perigee = 0.0 + val earthRadius = 6378.137 + val j3Harmonic = -2.53881E-6 + val twoPi = Math.PI * 2.0 + val twoThirds = 2.0 / 3.0 + val ck2 = 5.413079E-4 + val ck4 = 6.209887E-7 + var qoms24 = 0.0 + var s4 = 0.0 + val xke = 7.43669161E-2 + + fun willBeSeen(pos: GroundPos): Boolean { + return if (tle.meanmo < 1e-8) false else { + var lin = tle.incl + if (lin >= 90.0) lin = 180.0 - lin + val sma = 331.25 * exp(ln(1440.0 / tle.meanmo) * (2.0 / 3.0)) + val apogee = sma * (1.0 + tle.eccn) - earthRadius + acos(earthRadius / (apogee + earthRadius)) + lin * deg2Rad > abs(pos.lat * deg2Rad) + } + } + + fun getPredictor(pos: GroundPos): PassPredictor { + return PassPredictor(this, pos) + } + + fun getPosition(pos: GroundPos, time: Date): SatPos { + val satPos = SatPos() + // Date/time at which the position and velocity were calculated + val julUTC = calcCurrentDaynum(time) + 2444238.5 + // Convert satellite's epoch time to Julian and calculate time since epoch in minutes + val julEpoch = juliandDateOfEpoch(tle.epoch) + val tsince = (julUTC - julEpoch) * minPerDay + calculateSDP4orSGP4(tsince) + // Scale position and velocity vectors to km and km/sec + convertSatState(position, velocity) + // Calculate velocity of satellite + magnitude(velocity) + val squintVector = Vector4() + // Angles in rads, dist in km, vel in km/S. Calculate sat Az, El, Range and Range-rate. + calculateObs(julUTC, position, velocity, pos, squintVector, satPos) + calculateLatLonAlt(julUTC, satPos, position) + satPos.time = time + return satPos + } + + // Read the system clock and return the number of days since 31Dec79 00:00:00 UTC (daynum 0) + private fun calcCurrentDaynum(date: Date): Double { + val now = date.time + val sgp4Epoch = Calendar.getInstance(TimeZone.getTimeZone("UTC:UTC")) + sgp4Epoch.clear() + sgp4Epoch[1979, 11, 31, 0, 0] = 0 + val then = sgp4Epoch.timeInMillis + val millis = now - then + return millis / 1000.0 / 60.0 / 60.0 / 24.0 + } + + private fun juliandDateOfEpoch(epoch: Double): Double { + var year = floor(epoch * 1E-3) + val day = (epoch * 1E-3 - year) * 1000.0 + year = if (year < 57) year + 2000 else year + 1900 + return julianDateOfYear(year) + day + } + + fun julianDateOfYear(theYear: Double): Double { + val aYear = theYear - 1 + var i = floor(aYear / 100).toLong() + val a = i + i = a / 4 + val b = 2 - a + i + i = floor(365.25 * aYear).toLong() + i += (30.6001 * 14).toLong() + return i + 1720994.5 + b + } + + private fun calculateSDP4orSGP4(tsince: Double) { + if (tle.isDeepspace) (this as DeepSpaceSat).calculateSDP4(tsince) + else (this as NearEarthSat).calculateSGP4(tsince) + } + + // Converts the sat position and velocity vectors to km and km/sec + private fun convertSatState(pos: Vector4, vel: Vector4) { + scaleVector(earthRadius, pos) + scaleVector(earthRadius * minPerDay / secPerDay, vel) + } + + // Calculates the topocentric coordinates of the object with ECI pos and vel at time + private fun calculateObs( + julianUTC: Double, + positionVector: Vector4, + velocityVector: Vector4, + gsPos: GroundPos, + squintVector: Vector4, + satPos: SatPos + ) { + val obsPos = Vector4() + val obsVel = Vector4() + val range = Vector4() + val rgvel = Vector4() + val gsPosTheta = AtomicReference() + calculateUserPosVel(julianUTC, gsPos, gsPosTheta, obsPos, obsVel) + range.setXYZ( + positionVector.x - obsPos.x, + positionVector.y - obsPos.y, + positionVector.z - obsPos.z + ) + // Save these values globally for calculating squint angles later + squintVector.setXYZ(range.x, range.y, range.z) + rgvel.setXYZ( + velocityVector.x - obsVel.x, + velocityVector.y - obsVel.y, + velocityVector.z - obsVel.z + ) + magnitude(range) + val sinLat = sin(deg2Rad * gsPos.lat) + val cosLat = cos(deg2Rad * gsPos.lat) + val sinTheta = sin(gsPosTheta.get()) + val cosTheta = cos(gsPosTheta.get()) + val topS = sinLat * cosTheta * range.x + sinLat * sinTheta * range.y - cosLat * range.z + val topE = -sinTheta * range.x + cosTheta * range.y + val topZ = cosLat * cosTheta * range.x + cosLat * sinTheta * range.y + sinLat * range.z + var azim = atan(-topE / topS) + if (topS > 0.0) azim += Math.PI + if (azim < 0.0) azim += twoPi + satPos.azimuth = azim + satPos.elevation = asin(topZ / range.w) + satPos.range = range.w + satPos.rangeRate = dot(range, rgvel) / range.w + } + + // Returns the ECI position and velocity of the observer + private fun calculateUserPosVel( + time: Double, + gsPos: GroundPos, + gsPosTheta: AtomicReference, + obsPos: Vector4, + obsVel: Vector4 + ) { + val mFactor = 7.292115E-5 + gsPosTheta.set(mod2PI(thetaGJD(time) + deg2Rad * gsPos.lon)) + val c = + invert(sqrt(1.0 + flatFactor * (flatFactor - 2) * sqr(sin(deg2Rad * gsPos.lat)))) + val sq = sqr(1.0 - flatFactor) * c + val achcp = (earthRadius * c + gsPos.alt / 1000.0) * cos(deg2Rad * gsPos.lat) + obsPos.setXYZ( + achcp * cos(gsPosTheta.get()), achcp * sin(gsPosTheta.get()), + (earthRadius * sq + gsPos.alt / 1000.0) * sin(deg2Rad * gsPos.lat) + ) + obsVel.setXYZ(-mFactor * obsPos.y, mFactor * obsPos.x, 0.0) + magnitude(obsPos) + magnitude(obsVel) + } + + // Calculate the geodetic position of an object given its ECI pos and time + private fun calculateLatLonAlt( + time: Double, + satPos: SatPos, + position: Vector4 = this.position + ) { + satPos.theta = atan2(position.y, position.x) + satPos.longitude = mod2PI(satPos.theta - thetaGJD(time)) + val r = sqrt(sqr(position.x) + sqr(position.y)) + val e2 = flatFactor * (2.0 - flatFactor) + satPos.latitude = atan2(position.z, r) + var phi: Double + var c: Double + var i = 0 + var converged: Boolean + do { + phi = satPos.latitude + c = invert(sqrt(1.0 - e2 * sqr(sin(phi)))) + satPos.latitude = atan2(position.z + earthRadius * c * e2 * sin(phi), r) + converged = abs(satPos.latitude - phi) < epsilon + } while (i++ < 10 && !converged) + satPos.altitude = r / cos(satPos.latitude) - earthRadius * c + var temp = satPos.latitude + if (temp > Math.PI / 2.0) { + temp -= twoPi + satPos.latitude = temp + } + } + + fun calculatePosAndVel( + rk: Double, uk: Double, xnodek: Double, + xinck: Double, rdotk: Double, rfdotk: Double + ) { + // Orientation vectors + 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 + // Position and velocity + position.setXYZ(ux, uy, uz) + position.multiply(rk) + velocity.x = rdotk * ux + rfdotk * vx + velocity.y = rdotk * uy + rfdotk * vy + velocity.z = rdotk * uz + rfdotk * vz + } + + class Vector4 { + var w = 0.0 + var x = 0.0 + var y = 0.0 + var z = 0.0 + + fun multiply(multiplier: Double) { + x *= multiplier + y *= multiplier + z *= multiplier + } + + fun setXYZ(xValue: Double, yValue: Double, zValue: Double) { + x = xValue + y = yValue + z = zValue + } + } + + fun sqr(arg: Double): Double { + return arg * arg + } + + fun invert(value: Double): Double { + return 1.0 / value + } + + // Calculates the modulus of 2 * PI + fun mod2PI(value: Double): Double { + var retVal = value + val i = (retVal / twoPi).toInt() + retVal -= i * twoPi + if (retVal < 0.0) retVal += twoPi + return retVal + } + + // Solves Keplers' Equation + fun converge(temp: DoubleArray, axn: Double, ayn: Double, capu: Double) { + var converged = false + var i = 0 + do { + temp[7] = sin(temp[2]) + temp[8] = cos(temp[2]) + temp[3] = axn * temp[7] + temp[4] = ayn * temp[8] + temp[5] = axn * temp[8] + temp[6] = ayn * temp[7] + val epw = (capu - temp[4] + temp[3] - temp[2]) / (1.0 - temp[5] - temp[6]) + temp[2] + if (abs(epw - temp[2]) <= epsilon) converged = true else temp[2] = epw + } while (i++ < 10 && !converged) + } + + // Sets perigee and checks and adjusts the calculation if the perigee is less tan 156KM + fun setPerigee(perigee: Double) { + this.perigee = perigee + checkPerigee() + } + + // Checks and adjusts the calculation if the perigee is less tan 156KM + private fun checkPerigee() { + s4 = 1.012229 + qoms24 = 1.880279E-09 + if (perigee < 156.0) { + s4 = if (perigee <= 98.0) 20.0 else perigee - 78.0 + qoms24 = ((120 - s4) / earthRadius).pow(4.0) + s4 = s4 / earthRadius + 1.0 + } + } + + // Calculates the dot product of two vectors + private fun dot(v1: Vector4, v2: Vector4): Double { + return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z + } + + // Returns fractional part of double argument + private fun fraction(arg: Double): Double { + return arg - floor(arg) + } + + // Calculates scalar magnitude of a vector4 argument + private fun magnitude(v: Vector4) { + v.w = sqrt(sqr(v.x) + sqr(v.y) + sqr(v.z)) + } + + private fun modulus(arg1: Double): Double { + var returnValue = arg1 + val i = floor(returnValue / secPerDay).toInt() + returnValue -= i * secPerDay + if (returnValue < 0.0) returnValue += secPerDay + return returnValue + } + + // Multiplies the vector v1 by the scalar k + private fun scaleVector(k: Double, v: Vector4) { + v.multiply(k) + magnitude(v) + } + + private fun thetaGJD(theJD: Double): Double { + val earthRotPerSidDay = 1.00273790934 + val ut = fraction(theJD + 0.5) + val aJD = theJD - ut + val tu = (aJD - 2451545.0) / 36525.0 + var gmst = 24110.54841 + tu * (8640184.812866 + tu * (0.093104 - tu * 6.2E-6)) + gmst = modulus(gmst + secPerDay * earthRotPerSidDay * ut) + return twoPi * gmst / secPerDay + } + + companion object { + + fun createSat(tle: TLE): Satellite { + return if (tle.isDeepspace) DeepSpaceSat(tle) + else NearEarthSat(tle) + } + + fun importTLE(tleStream: InputStream): List { + val importedTles = mutableListOf() + val currentTLE = arrayOf(String(), String(), String()) + var line = 0 + tleStream.bufferedReader().forEachLine { + if (line != 2) { + currentTLE[line] = it + line++ + } else { + currentTLE[line] = it + importedTles.add(parseTLE(currentTLE)) + line = 0 + } + } + return importedTles + } + + private fun parseTLE(tle: Array): TLE { + val name: String = tle[0].trim() + val epoch: Double = tle[1].substring(18, 32).toDouble() + val meanmo: Double = tle[2].substring(52, 63).toDouble() + val eccn: Double = 1.0e-07 * tle[2].substring(26, 33).toDouble() + val incl: Double = tle[2].substring(8, 16).toDouble() + val raan: Double = tle[2].substring(17, 25).toDouble() + val argper: Double = tle[2].substring(34, 42).toDouble() + val meanan: Double = tle[2].substring(43, 51).toDouble() + val catnum: Int = tle[1].substring(2, 7).trim().toInt() + val bstar: Double = 1.0e-5 * tle[1].substring(53, 59).toDouble() / + 10.0.pow(tle[1].substring(60, 61).toDouble()) + return TLE(name, epoch, meanmo, eccn, incl, raan, argper, meanan, catnum, bstar) + } + } +} diff --git a/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/TLE.kt b/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/TLE.kt new file mode 100644 index 00000000..0dd2eee1 --- /dev/null +++ b/core/src/main/java/com/rtbishop/look4sat/predict4kotlin/TLE.kt @@ -0,0 +1,38 @@ +/* + * Look4Sat. Amateur radio satellite tracker and pass predictor. + * Copyright (C) 2019-2021 Arty Bishop (bishop.arty@gmail.com) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ +package com.rtbishop.look4sat.predict4kotlin + +data class TLE( + val name: String, + val epoch: Double, + val meanmo: Double, + val eccn: Double, + val incl: Double, + val raan: Double, + val argper: Double, + val meanan: Double, + val catnum: Int, + val bstar: Double +) { + val xincl: Double = Math.toRadians(incl) + val xnodeo: Double = Math.toRadians(raan) + val omegao: Double = Math.toRadians(argper) + val xmo: Double = Math.toRadians(meanan) + val xno: Double = meanmo * Math.PI * 2.0 / 1440 + val isDeepspace: Boolean = meanmo < 6.4 +} diff --git a/settings.gradle b/settings.gradle index 1c788441..51e1b291 100644 --- a/settings.gradle +++ b/settings.gradle @@ -1,2 +1,3 @@ +include ':core' include ':app' rootProject.name='Look4Sat'