Skip to content

Commit 4131a89

Browse files
mihbormihbor
committed
los, solarPosition and Example10_09
1 parent 455a1e1 commit 4131a89

7 files changed

Lines changed: 263 additions & 12 deletions

File tree

‎core/src/commonMain/kotlin/ltd/mbor/sciko/orbital/LoS.kt‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
package ltd.mbor.sciko.orbital
2+
3+
import ltd.mbor.sciko.linalg.norm
4+
import org.jetbrains.kotlinx.multik.api.linalg.dot
5+
import org.jetbrains.kotlinx.multik.ndarray.data.D1
6+
import org.jetbrains.kotlinx.multik.ndarray.data.MultiArray
7+
8+
fun los(rSat: MultiArray<Double, D1>, rSun: MultiArray<Double, D1>, rPlanet: Double = rEarth): Boolean {
9+
10+
val theta = acos((rSat dot rSun) / rSat.norm() / rSun.norm())
11+
val thetaSat = acos(rPlanet / rSat.norm())
12+
val thetaSun = acos(rPlanet / rSun.norm())
13+
14+
return thetaSat + thetaSun <= theta
15+
}

‎core/src/commonMain/kotlin/ltd/mbor/sciko/orbital/RKF45.kt‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,8 @@ import kotlin.math.max
1414
import kotlin.math.min
1515
import kotlin.math.pow
1616

17+
private const val EPS = 2.2204e-16
18+
1719
private val a = mk.ndarray(mk[0.0, 1.0 / 4, 3.0 / 8, 12.0 / 13, 1.0, 1.0 / 2])
1820
private val b = mk.ndarray(mk[
1921
mk[ 0.0, 0.0, 0.0, 0.0, 0.0],
@@ -42,7 +44,7 @@ fun rkf45(
4244
tspan: ClosedRange<Double>,
4345
y0: MultiArray<Double, D1>,
4446
tolerance: Double = 1e-8,
45-
initialStep: Double = (tspan.endInclusive - tspan.start)/100,
47+
h0: Double = (tspan.endInclusive - tspan.start)/100,
4648
outerFunction: (Double, MultiArray<Double, D1>) -> MultiArray<Double, D1> = { _, y -> y },
4749
terminateFunction: TerminateFunction? = null,
4850
odeFunction: (Double, MultiArray<Double, D1>) -> MultiArray<Double, D1>,
@@ -54,7 +56,7 @@ fun rkf45(
5456
var y = y0
5557
val tOut = mutableListOf(t)
5658
val yOut = mutableListOf(y)
57-
var h = initialStep
59+
var h = h0
5860

5961
while (t < tf) {
6062
val hmin = 16 * Double.MIN_VALUE
@@ -75,7 +77,7 @@ fun rkf45(
7577
val teMax = te.map { abs(it) }.max()!!
7678
val ymax = y.map { abs(it) }.max()!!
7779
val teAllowed = tolerance * max(ymax, 1.0)
78-
val delta = (teAllowed / (teMax + Double.MIN_VALUE)).pow(1.0 / 5)
80+
val delta = (teAllowed / (teMax + EPS)).pow(0.2)
7981

8082
if (teMax <= teAllowed) {
8183
h = min(h, tf - t)

‎core/src/commonMain/kotlin/ltd/mbor/sciko/orbital/SolarPosition.kt‎

Lines changed: 49 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,49 @@
1+
package ltd.mbor.sciko.orbital
2+
import org.jetbrains.kotlinx.multik.api.mk
3+
import org.jetbrains.kotlinx.multik.api.ndarray
4+
import org.jetbrains.kotlinx.multik.ndarray.data.D1
5+
import org.jetbrains.kotlinx.multik.ndarray.data.MultiArray
6+
import org.jetbrains.kotlinx.multik.ndarray.operations.times
7+
import kotlin.math.PI
8+
import kotlin.math.cos
9+
import kotlin.math.sin
10+
11+
fun solarPosition(jd: Double): Triple<Double, Double, MultiArray<Double, D1>> {
12+
// Astronomical unit (km)
13+
val au = 149597870.691
14+
15+
// Julian days since J2000
16+
val n = jd - 2451545
17+
18+
// Julian centuries since J2000
19+
// val cy = n / 36525
20+
21+
// Mean anomaly (radians)
22+
val m = (357.528 + 0.9856003 * n).degrees.mod(2*PI)
23+
val mDeg = (357.528 + 0.9856003 * n).mod(360.0)
24+
25+
// Mean longitude (radians)
26+
val l = (280.460 + 0.98564736 * n).degrees.mod(2*PI)
27+
val lDeg = (280.460 + 0.98564736 * n).mod(360.0)
28+
29+
// Apparent ecliptic longitude (radians)
30+
val lambda = (l + 1.915.degrees * sin(m) + 0.020.degrees * sin(2 * m)).mod(2*PI)
31+
val lambdaDeg = lambda.toDegrees()
32+
33+
// Obliquity of the ecliptic (radians)
34+
val eps = (23.439 - 0.0000004 * n).degrees.mod(2*PI)
35+
val epsDeg = (23.439 - 0.0000004 * n).mod(360.0)
36+
37+
val (sl, ce) = sin(lambda) to cos(eps)
38+
39+
// Unit vector from earth to sun
40+
val u = mk.ndarray(mk[cos(lambda), sin(lambda) * cos(eps), sin(lambda) * sin(eps)])
41+
42+
// Distance from earth to sun (km)
43+
val rS = (1.00014 - 0.01671 * cos(m) - 0.000140 * cos(2 * m)) * au
44+
45+
// Geocentric position vector (km)
46+
val rSVector = rS * u
47+
48+
return Triple(lambda, eps, rSVector)
49+
}

‎examples/src/jvmMain/kotlin/ltd/mbor/sciko/orbital/examples/Example10_01.kt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -79,7 +79,7 @@ fun main() {
7979
t0..tf,
8080
y0,
8181
tolerance = tolerance,
82-
initialStep = initialStep,
82+
h0 = initialStep,
8383
terminateFunction = terminate,
8484
) { t, y -> rates(t, y, mu, RE, wE, CD, m, A) }
8585

‎examples/src/jvmMain/kotlin/ltd/mbor/sciko/orbital/examples/Example10_02.kt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -86,7 +86,7 @@ fun main() {
8686
t0s..t,
8787
dy0,
8888
tolerance = 1e-15,
89-
initialStep = delT / 10.0,
89+
h0 = delT / 10.0,
9090
) { ti, f -> ratesEnckeJ2(ti, f, R0s, V0s, t0s, mu, RE, J2c) }
9191

9292
val z = yOut.last()

‎examples/src/jvmMain/kotlin/ltd/mbor/sciko/orbital/examples/Example10_06.kt‎

Lines changed: 2 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -8,12 +8,7 @@ import org.jetbrains.kotlinx.multik.api.ndarray
88
import org.jetbrains.kotlinx.multik.ndarray.data.D1
99
import org.jetbrains.kotlinx.multik.ndarray.data.MultiArray
1010
import org.jetbrains.kotlinx.multik.ndarray.data.get
11-
import kotlin.math.PI
12-
import kotlin.math.abs
13-
import kotlin.math.cos
14-
import kotlin.math.pow
15-
import kotlin.math.sin
16-
import kotlin.math.sqrt
11+
import kotlin.math.*
1712

1813
/**
1914
* Kotlin port of MATLAB Example_10_6.m
@@ -58,7 +53,7 @@ fun main() {
5853
t0..tf,
5954
y0,
6055
tolerance = 1e-12,
61-
initialStep = T0 / 1000.0,
56+
h0 = T0 / 1000.0,
6257
odeFunction = { t, y -> ratesGaussJ2(t, y, mu, RE, J2c) }
6358
)
6459

‎examples/src/jvmMain/kotlin/ltd/mbor/sciko/orbital/examples/Example10_09.kt‎

Lines changed: 190 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,190 @@
1+
package ltd.mbor.sciko.orbital.examples
2+
3+
import ltd.mbor.sciko.linalg.norm
4+
import ltd.mbor.sciko.orbital.*
5+
import org.jetbrains.kotlinx.kandy.letsplot.export.save
6+
import org.jetbrains.kotlinx.multik.api.KEEngineType
7+
import org.jetbrains.kotlinx.multik.api.mk
8+
import org.jetbrains.kotlinx.multik.api.ndarray
9+
import org.jetbrains.kotlinx.multik.ndarray.data.D1
10+
import org.jetbrains.kotlinx.multik.ndarray.data.MultiArray
11+
import org.jetbrains.kotlinx.multik.ndarray.data.get
12+
import kotlin.math.*
13+
14+
/**
15+
* Kotlin port of MATLAB Example_10_09.m
16+
*
17+
* Solves Gauss planetary equations with solar radiation pressure (SRP)
18+
* over 3 years, including Earth eclipse (line-of-sight) shadow function.
19+
*/
20+
fun main() {
21+
mk.setEngine(KEEngineType)
22+
23+
// Conversion factors
24+
val hours = 3600.0
25+
val days = 24.0 * hours
26+
27+
// Constants
28+
val mu = muEarth // km^3/s^2
29+
val c = 2.998e8 // m/s
30+
val S = 1367.0 // W/m^2
31+
val Psr = S / c // Pa = N/m^2
32+
33+
// Spacecraft properties
34+
val CR = 2.0 // radiation pressure coefficient
35+
val m = 100.0 // kg
36+
val As = 200.0 // m^2
37+
38+
// Initial orbital parameters (given)
39+
val a0 = 10085.44 // km
40+
val e0 = 0.025422
41+
val i0 = 88.3924.degrees
42+
val RA0 = 45.38124.degrees
43+
val TA0 = 343.4268.degrees
44+
val w0 = 227.493.degrees
45+
46+
// Inferred parameters
47+
val h0 = sqrt(mu * a0 * (1 - e0.pow(2)))
48+
val T0 = 2 * PI / sqrt(mu) * a0.pow(1.5)
49+
50+
// Initial element state vector
51+
val coe0 = mk.ndarray(mk[h0, e0, RA0, i0, w0, TA0])
52+
53+
// Time span
54+
val JD0 = 2438400.5 // 6 January 1964 0 UT
55+
val t0 = 0.0
56+
val tf = 3.0 * 365.0 * days
57+
58+
val (tOut, yOut) = rkf45(
59+
t0..tf,
60+
coe0,
61+
tolerance = 1e-10,
62+
h0 = T0 / 1000.0,
63+
) { t, y ->
64+
ratesSRP(t, y, mu, JD0, days, Psr, CR, As, m)
65+
}
66+
67+
val h = yOut.map { it[0] }
68+
val e = yOut.map { it[1] }
69+
val RA = yOut.map { it[2] }
70+
val inc = yOut.map { it[3] }
71+
val w = yOut.map { it[4] }
72+
val TA = yOut.map { it[5] }
73+
74+
// a = h^2 / mu / (1 - e^2)
75+
val a = yOut.map { it[0].pow(2) / mu / (1 - it[1].pow(2)) }
76+
77+
val tDays = tOut.map { it / days }
78+
79+
// Plots similar to MATLAB example
80+
plotScalar(
81+
Triple(tDays, h.map { it - h0 }, "Δh (km^2/s)")
82+
, xLabel = "days", yLabel = "Δh (km^2/s)").save("Example10_09_h.png")
83+
84+
plotScalar(
85+
Triple(tDays, e.map { it - e0 }, "Δe")
86+
, xLabel = "days", yLabel = "Δe").save("Example10_09_e.png")
87+
88+
plotScalar(
89+
Triple(tDays, RA.map { (it - RA0).toDegrees() }, "ΔΩ (deg)")
90+
, xLabel = "days", yLabel = "deg").save("Example10_09_RA.png")
91+
92+
plotScalar(
93+
Triple(tDays, inc.map { (it - i0).toDegrees() }, "Δi (deg)")
94+
, xLabel = "days", yLabel = "deg").save("Example10_09_i.png")
95+
96+
plotScalar(
97+
Triple(tDays, w.map { (it - w0).toDegrees() }, "Δω (deg)")
98+
, xLabel = "days", yLabel = "deg").save("Example10_09_w.png")
99+
100+
plotScalar(
101+
Triple(tDays, a.map { it - a0 }, "Δa (km)")
102+
, xLabel = "days", yLabel = "km").save("Example10_09_a.png")
103+
}
104+
105+
private fun ratesSRP(
106+
t: Double,
107+
f: MultiArray<Double, D1>,
108+
mu: Double,
109+
JD0: Double,
110+
days: Double,
111+
Psr: Double,
112+
CR: Double,
113+
As: Double,
114+
m: Double,
115+
): MultiArray<Double, D1> {
116+
val h = f[0]
117+
val e = f[1]
118+
val RA = f[2]
119+
val i = f[3]
120+
val w = f[4]
121+
val TA = f[5]
122+
123+
val u = w + TA // argument of latitude
124+
125+
// State vectors from COEs at current time
126+
val coe = f
127+
val (R, V) = svFromCoe(coe, mu)
128+
val r = R.norm()
129+
130+
// Sun position at current Julian date
131+
val JD = JD0 + t / days
132+
val (lambda, eps, rSun) = solarPosition(JD) // km, earth->sun
133+
134+
// Eclipse shadow function (1 in sunlight, 0 in shadow)
135+
val nu = if (los(R, rSun)) 0.0 else 1.0
136+
137+
// Solar radiation pressure acceleration magnitude (km/s^2)
138+
val pSR = nu * (Psr * CR * As / m) / 1000.0
139+
140+
// // Compute RSW frame unit vectors
141+
// val rHat = R / r
142+
// val hVec = R cross V
143+
// val wHat = hVec / hVec.norm()
144+
// val sHat = wHat cross rHat
145+
//
146+
// // Unit vector from Earth to Sun
147+
// val uSun = rSun / rSun.norm()
148+
//
149+
// // Components of sun direction in RSW frame
150+
// val ur = (uSun dot rHat)
151+
// val us = (uSun dot sHat)
152+
// val uw = (uSun dot wHat)
153+
//
154+
// val si = sin(i)
155+
156+
val sl = sin(lambda); val cl = cos(lambda)
157+
val se = sin(eps); val ce = cos(eps)
158+
val sW = sin(RA); val cW = cos(RA)
159+
val si = sin(i); val ci = cos(i)
160+
val su = sin(u); val cu = cos(u)
161+
val sT = sin(TA); val cT = cos(TA)
162+
163+
val ur = sl*ce*cW*ci*su + sl*ce*sW*cu - cl*sW*ci*su + cl*cW*cu + sl*se*si*su
164+
165+
val us = sl*ce*cW*ci*cu - sl*ce*sW*su - cl*sW*ci*cu - cl*cW*su + sl*se*si*cu
166+
167+
val uw = - sl*ce*cW*si + cl*sW*si + sl*se*ci
168+
169+
// Gauss planetary equations with SRP (matching MATLAB lines 178-191)
170+
val hdot = -pSR * r * us
171+
172+
val edot = -pSR * (
173+
(h / mu) * sT * ur + (1.0 / (mu * h)) * ((h.pow(2) + mu * r) * cT + mu * e * r) * us
174+
)
175+
176+
val TAdot = h / r.pow(2) - pSR / (e * h) * (
177+
h.pow(2) / mu * cT * ur - (r + h.pow(2) / mu) * sT * us
178+
)
179+
180+
val RAdot = -pSR * r / h / si * su * uw
181+
182+
val idot = -pSR * r / h * cu * uw
183+
184+
val wdot = -pSR * (
185+
-1.0 / (e * h) * ((h * h / mu) * cT * ur - (r + h * h / mu) * sT * us) -
186+
r * sin(u) / h / si * cos(i) * uw
187+
)
188+
189+
return mk.ndarray(mk[hdot, edot, RAdot, idot, wdot, TAdot])
190+
}

0 commit comments

Comments
 (0)