Space mission engineer: MQL recipes

Recipes

Orbit of the ISS
Hohmann transfer from LEO to GEO
Payload of an upper stage with the rocket equation
Sun-synchronous orbit and the geopotential
Position on a Molniya orbit with the Kepler equation
Slant range, elevation and pass duration over a ground station
Hyperbolic departure orbit to Mars
Orbital elements from a state vector
Ku-band downlink budget and bit error rate
Attitude matrix, Euler angle and inertia tensor

executed on MentDB

Orbit of the ISS

For a circular orbit at 420 km, compute the orbital velocity, the period, the number of orbits per day and the escape velocity at that altitude.

Pipeline

01 · radiusAltitude to orbit radiusunit calc→02 · orbitv, T and v_esc from mumath calc→03 · convertm/s to km/h, s to minunit convert

Script

admin

#1. Earth gravitational parameter and ISS orbit radius at 420 km altitude;
-> "[mu]" 3.986004418e14;
-> "[r]" (unit calc "6378.137 km + 420 km" "m");
#2. Circular velocity, period and escape velocity at that radius;
-> "[v]" (math calc "sqrt([mu]/[r])");
-> "[T]" (unit convert (math calc "2*pi*sqrt([r]^3/[mu])") "s" "min");
-> "[vesc]" (math calc "sqrt(2*[mu]/[r])");
#3. Orbits per day and speed in km/h;
concat "v = " (math round [v] 1) " m/s = " (math round (unit convert [v] "m/s" "km/h") 0) " km/h, T = " (math round [T] 2) " min, " (math round (math calc "1440/[T]") 2) " orbits/day, v_esc = " (math round (math calc "[vesc]/1000") 3) " km/s"

Result

mentdb

v = 7657.3 m/s = 27566 km/h, T = 92.97 min, 15.49 orbits/day, v_esc = 10.829 km/s

Go further

Add the J2 nodal precession of the orbit plane
Compute the drag decay with num ode

executed on MentDB

Hohmann transfer from LEO to GEO

Compute the two burns of a Hohmann transfer from a 300 km parking orbit to the geostationary orbit, the total delta-v and the time of flight.

Pipeline

01 · radiiBoth orbit radiiunit calc→02 · burnsdv1 and dv2math calc→03 · timeHalf period of the ellipseunit convert

Script

admin

#1. Radii of the parking orbit at 300 km and of the geostationary orbit;
-> "[mu]" 3.986004418e14;
-> "[r1]" (unit calc "6378.137 km + 300 km" "m");
-> "[r2]" (unit calc "6378.137 km + 35786 km" "m");
#2. Hohmann transfer: burn at perigee, then circularisation at apogee;
-> "[dv1]" (math calc "sqrt([mu]/[r1])*(sqrt(2*[r2]/([r1]+[r2]))-1)");
-> "[dv2]" (math calc "sqrt([mu]/[r2])*(1-sqrt(2*[r1]/([r1]+[r2])))");
#3. Time of flight: half the period of the transfer ellipse;
-> "[tof]" (unit convert (math calc "pi*sqrt((([r1]+[r2])/2)^3/[mu])") "s" "h");
concat "dv1 = " (math round [dv1] 1) " m/s, dv2 = " (math round [dv2] 1) " m/s, total = " (math round (math calc "[dv1]+[dv2]") 1) " m/s, transfer time = " (math round [tof] 2) " h"

Result

mentdb

dv1 = 2425.7 m/s, dv2 = 1466.8 m/s, total = 3892.6 m/s, transfer time = 5.28 h

Go further

Add the plane change from 28.5 deg at apogee
Compare with a bi-elliptic transfer

executed on MentDB

Payload of an upper stage with the rocket equation

With the Tsiolkovsky equation, get the delta-v of an upper stage, then solve for the heaviest payload that still leaves 6 km/s.

Pipeline

01 · dvdv = Isp g0 ln m0/mfmath calc→02 · buildEquation with the payloadconcat→03 · solveLargest payloadnum solve

Script

admin

#1. Upper stage: dry mass 4500 kg, propellant 107500 kg, Isp 348 s;
-> "[mdry]" 4500;
-> "[mprop]" 107500;
-> "[isp]" 348;
#2. Tsiolkovsky: delta-v with payload m, ve = Isp*g0;
-> "[dv0]" (math calc "[isp]*_g0*ln(([mdry]+[mprop])/[mdry])");
#3. Largest payload for a delta-v of 6000 m/s;
-> "[pl]" (num solve (concat [isp] "*9.80665*ln((" [mdry] "+" [mprop] "+m)/(" [mdry] "+m))-6000") "m" 0 50000);
concat "dv without payload = " (math round [dv0] 0) " m/s, max payload for 6 km/s = " (math round [pl] 0) " kg, payload fraction = " (math round (math calc "100*[pl]/([mdry]+[mprop]+[pl])") 2) " %"

Result

mentdb

dv without payload = 10970 m/s, max payload for 6 km/s = 17888 kg, payload fraction = 13.77 %

Go further

Chain two stages and optimise the mass split with num maximize
Propagate the Isp uncertainty with num uncertainty

executed on MentDB

Sun-synchronous orbit and the geopotential

The J2 flattening of the Earth makes the orbit plane precess. The script finds the inclination that gives a sun-synchronous orbit at 700 km, and the size of the zonal J2 and tesseral C22 terms of the geopotential with Legendre functions.

Pipeline

01 · orbitMean motionmath calc→02 · ssoInclination for 360 deg/yearnum solve→03 · zonalJ2 term with P2legendre(n, x)→04 · tesseralC22 term with P22assoc_legendre(n, m, x)

Script

admin

#1. Earth: mu, equatorial radius, zonal J2 and tesseral C22, circular orbit at 700 km;
-> "[mu]" 3.986004418e14; -> "[Re]" 6378137; -> "[J2]" 1.08263e-3; -> "[C22]" 1.5744e-6;
-> "[a]" (math calc "[Re]+700e3");
-> "[n]" (math calc "sqrt([mu]/[a]^3)");
#2. Sun-synchronous: J2 nodal precession must equal 360 deg per tropical year;
-> "[need]" (math calc "2*pi/(365.2422*86400)");
-> "[i]" (num solve (concat "-1.5*" [n] "*" [J2] "*(" [Re] "/" [a] ")^2*cos(x)-" [need]) "x" 1.6 2.0);
#3. Relative size of the J2 term of the geopotential at 45 deg latitude: J2 Re/a ^2 P2 of sin phi;
-> "[dJ2]" (math calc "-[J2]*([Re]/[a])^2*legendre(2, sin(45*pi/180))");
#4. Tesseral C22 term at the equator, longitude 0: C22 Re/a ^2 P22 of sin phi;
-> "[dC22]" (math calc "[C22]*([Re]/[a])^2*assoc_legendre(2, 2, 0)");
concat "{\"period_min\": " (math round (math calc "2*pi/[n]/60") 2) ", \"sun_sync_inclination_deg\": " (math round (math calc "[i]*180/pi") 3) ", \"J2_term_45deg\": " (math round_sig [dJ2] 4) ", \"C22_term_equator\": " (math round_sig [dC22] 4) "}"

Result

mentdb

{"period_min": 98.77, "sun_sync_inclination_deg": 98.188, "J2_term_45deg": -0.0002198, "C22_term_equator": 0.000003835}

Go further

Tabulate the sun-synchronous inclination from 400 to 1000 km with a for loop.
Add the J3 term with legendre(3, x): it creates the frozen orbits.

executed on MentDB

Position on a Molniya orbit with the Kepler equation

A Molniya orbit (a = 26600 km, e = 0.74) is used for high-latitude coverage. Three hours after perigee, solve the Kepler equation for the eccentric anomaly, then get the true anomaly, the altitude, the vis-viva speed and the position in the perifocal frame.

Pipeline

01 · convertSemi-major axis and time in SIunit convert→02 · anomalyMean anomaly n tmath calc→03 · keplerE - e sin E = Mnum solve→04 · positionPerifocal position as r e^i nucomplex from_polar→05 · splitx and y in kmcomplex re

Script

admin

#1. Molniya orbit: semi-major axis, eccentricity, Earth gravitational parameter;
-> "[mu]" 3.986004418e14;
-> "[a]" (unit convert 26600 "km" "m");
-> "[e]" 0.74;
#2. Mean anomaly 3 hours after perigee passage;
-> "[t]" (unit convert 3 "h" "s");
-> "[M]" (math calc "sqrt([mu]/[a]^3)*[t]");
#3. Kepler equation E - e sin E = M for the eccentric anomaly;
-> "[E]" (num solve (concat "x-" [e] "*sin(x)-" [M]) "x" 0 6.283185307179586);
#4. True anomaly, radius and vis-viva speed;
-> "[nu]" (math calc "2*atan(sqrt((1+[e])/(1-[e]))*tan([E]/2))");
-> "[r]" (math calc "[a]*(1-[e]*cos([E]))");
-> "[v]" (math calc "sqrt([mu]*(2/[r]-1/[a]))");
#5. Position in the perifocal frame as a complex number in km;
-> "[pos]" (complex from_polar (/ [r] 1000) [nu]);
concat "period " (math round (unit convert (math calc "2*pi*sqrt([a]^3/[mu])") "s" "h") 3) " h | M = " (math round [M] 4) " rad, E = " (math round [E] 4) " rad, true anomaly = " (math round (math calc "[nu]*180/pi") 2) " deg | altitude " (math round (math calc "([r]-6378137)/1000") 0) " km, speed " (math round [v] 0) " m/s | perifocal x = " (math round (complex re [pos]) 0) " km, y = " (math round (complex im [pos]) 0) " km"

Result

mentdb

period 11.993 h | M = 1.5717 rad, E = 2.179 rad, true anomaly = 157.17 deg | altitude 31469 km, speed 2465 m/s | perifocal x = -34883 km, y = 14683 km

Go further

Loop over t with num table to tabulate the ground track over one period.
Propagate the same orbit with num ode and the two-body equations to check the result.

executed on MentDB

Slant range, elevation and pass duration over a ground station

The tracking software gives the satellite in the east-north-up frame of the Kourou station. The script gets the slant range, the elevation and the azimuth without quadrant error, then the Earth central angle and the longest possible pass of a 550 km orbit above a 10 deg elevation mask.

Pipeline

01 · rangeSlant rangemath hypot→02 · anglesElevation and azimuthmath atan2→03 · coverageEarth central anglemath calc→04 · degreesRadians to degreesmath rad_to_deg

Script

admin

#1. Satellite seen from the Kourou station, local east-north-up vector in km;
-> "[e]" -412; -> "[n]" 830; -> "[u]" 260;
#2. Horizontal distance, slant range, elevation and azimuth;
-> "[hor]" (math hypot [e] [n]);
-> "[rho]" (math hypot [hor] [u]);
-> "[el]" (math rad_to_deg (math atan2 [u] [hor]));
-> "[az]" (math calc (concat "mod(" (math rad_to_deg (math atan2 [e] [n])) "+360, 360)"));
#3. Longest pass for a 550 km orbit and a 10 deg mask: Earth central angle;
-> "[Re]" 6378.137; -> "[h]" 550; -> "[mask]" (math deg_to_rad 10);
-> "[lambda]" (math calc "acos([Re]*cos([mask])/([Re]+[h]))-[mask]");
-> "[P]" (math calc "2*pi*sqrt(([Re]+[h])^3/398600.4418)");
concat "slant range = " (math round [rho] 1) " km, elevation = " (math round [el] 2) " deg, azimuth = " (math round [az] 2) " deg, Earth central angle = " (math round (math rad_to_deg [lambda]) 2) " deg, max pass = " (math round (math calc "2*[lambda]/(2*pi)*[P]/60") 1) " min"

Result

mentdb

slant range = 962.4 km, elevation = 15.67 deg, azimuth = 333.6 deg, Earth central angle = 14.96 deg, max pass = 7.9 min

Go further

Add the rotation of the Earth to the pass duration for a prograde orbit.
Build the visibility table of a day with num table on the central angle.

executed on MentDB

Hyperbolic departure orbit to Mars

A Mars probe leaves a 300 km parking orbit with a hyperbolic excess speed of 2.95 km/s. The script gets C3, the injection delta-v, the eccentricity and turning angle of the hyperbola, then the hyperbolic anomaly, the time and the true anomaly when the probe leaves the sphere of influence of the Earth.

Pipeline

01 · injectionPerigee speed and delta-vmath calc→02 · anomalyHyperbolic anomalymath acosh→03 · timeHyperbolic Kepler equationmath sinh→04 · trueTrue anomaly with tanhmath tanh→05 · arcInverse tangentmath atan

Script

admin

#1. Departure to Mars: v infinity 2.95 km/s from a 300 km parking orbit;
-> "[mu]" 398600.4418; -> "[rp]" (math calc "6378.137+300"); -> "[vinf]" 2.95;
#2. Speed at perigee, injection delta-v, eccentricity and turning angle;
-> "[vp]" (math calc "sqrt([vinf]^2+2*[mu]/[rp])");
-> "[dv]" (math calc "[vp]-sqrt([mu]/[rp])");
-> "[ecc]" (math calc "1+[rp]*[vinf]^2/[mu]");
-> "[a]" (math calc "[mu]/[vinf]^2");
#3. Hyperbolic anomaly at the edge of the sphere of influence, 925000 km;
-> "[F]" (math acosh (math calc "(925000/[a]+1)/[ecc]"));
-> "[t]" (math calc (concat "sqrt([a]^3/[mu])*([ecc]*" (math sinh [F]) "-[F])"));
#4. True anomaly from tan nu/2 = sqrt e+1 / e-1 tanh F/2;
-> "[nu]" (math calc (concat "2*" (math atan (math calc (concat "sqrt(([ecc]+1)/([ecc]-1))*" (math tanh (math calc "[F]/2")))))));
concat "C3 = " (math round (math calc "[vinf]^2") 2) " km2/s2, delta-v = " (math round [dv] 4) " km/s, e = " (math round [ecc] 4) ", turning angle = " (math round (math rad_to_deg (math calc "2*asin(1/[ecc])")) 2) " deg, F = " (math round [F] 4) ", time to SOI = " (math round (unit convert [t] "s" "d") 2) " d, true anomaly = " (math round (math rad_to_deg [nu]) 2) " deg"

Result

mentdb

C3 = 8.7 km2/s2, delta-v = 3.5914 km/s, e = 1.1458, turning angle = 121.56 deg, F = 3.6101, time to SOI = 3.15 d, true anomaly = 149.23 deg

Go further

Solve the hyperbolic Kepler equation e sinh F - F = M with num solve for a given date after injection.
Find with num minimize the parking orbit altitude that minimises the delta-v for this C3.

executed on MentDB

Orbital elements from a state vector

From the position and velocity of a spacecraft, the script computes the angular momentum, the node line and the eccentricity vector with vector algebra, then the six classical orbital elements, the right ascension of the node with atan2.

Pipeline

01 · momentumh = r x vvector cross→02 · eccentricityEccentricity vectorvector sub→03 · elementsInclination and nodemath atan2→04 · anglesPerigee and true anomalyvector angle→05 · normsNorms of h and evector norm

Script

admin

#1. State vector of a spacecraft in the geocentric equatorial frame, km and km/s;
-> "[r]" "[-6045, -3490, 2500]"; -> "[v]" "[-3.457, 6.618, 2.533]"; -> "[mu]" 398600;
#2. Angular momentum h = r x v and node line N = k x h;
-> "[h]" (vector cross [r] [v]);
-> "[N]" (vector cross "[0, 0, 1]" [h]);
#3. Eccentricity vector e = v x h / mu - r / norm r;
-> "[ev]" (vector sub (vector scale (vector cross [v] [h]) (math calc "1/[mu]")) (vector unit [r]));
#4. Inclination, right ascension of the node, argument of perigee and true anomaly;
-> "[i]" (math rad_to_deg (math acos (math calc (concat (vector dot [h] "[0, 0, 1]") "/" (vector norm [h])))));
-> "[raan]" (math calc (concat "mod(" (math rad_to_deg (math atan2 (vector dot [N] "[0, 1, 0]") (vector dot [N] "[1, 0, 0]"))) "+360, 360)"));
-> "[w]" (math rad_to_deg (vector angle [N] [ev]));
-> "[nu]" (math rad_to_deg (vector angle [ev] [r]));
concat "h = " (math round (vector norm [h]) 0) " km2/s, i = " (math round [i] 2) " deg, RAAN = " (math round [raan] 2) " deg, e = " (math round (vector norm [ev]) 4) ", argument of perigee = " (math round [w] 2) " deg, true anomaly = " (math round [nu] 2) " deg"

Result

mentdb

h = 58312 km2/s, i = 153.25 deg, RAAN = 255.28 deg, e = 0.1712, argument of perigee = 20.07 deg, true anomaly = 28.45 deg

Go further

The true anomaly from vector angle is in [0, 180]: use the sign of r.v to get the other half of the orbit.
Propagate the orbit with num ode and check that h and e stay constant.

executed on MentDB

Ku-band downlink budget and bit error rate

A GEO satellite sends 20 Mbit/s in Ku band to a small station. The script computes the free space loss, the C/N0 with the Boltzmann constant, the Eb/N0, the bit error rate of BPSK with the complementary error function, and the Eb/N0 required for a BER of 1e-6 with the inverse error function.

Pipeline

01 · constantsc and kconst get→02 · lossFree space path lossmath calc→03 · berBER = 0.5 erfcmath erfc→04 · requiredRequired Eb/N0math erf_inv

Script

admin

#1. Ku-band downlink from a GEO satellite: 12 GHz, slant range 38000 km, EIRP 48 dBW, station G/T 15 dB/K, 1.5 dB of rain and pointing losses;
-> "[f]" 12e9; -> "[d]" 38000e3; -> "[eirp]" 48; -> "[gt]" 15;
-> "[c]" (const get "c"); -> "[k]" (const get "k");
#2. Free space path loss and carrier to noise density ratio C/N0 in dBHz;
-> "[fspl]" (math calc "20*log10(4*pi*[d]*[f]/[c])");
-> "[cn0]" (math calc "[eirp]-[fspl]+[gt]-10*log10([k])-1.5");
#3. Eb/N0 for 20 Mbit/s and bit error rate of BPSK, 0.5 erfc sqrt Eb/N0;
-> "[ebn0]" (math calc "[cn0]-10*log10(20e6)");
-> "[ber]" (math calc (concat "0.5*" (math erfc (math calc "sqrt(10^([ebn0]/10))"))));
#4. Eb/N0 needed for a BER of 1e-6 with the inverse error function, and the margin;
-> "[x]" (math erf_inv (math calc "1-2*1e-6"));
-> "[req]" (math calc "10*log10([x]^2)");
concat "FSPL = " (math round [fspl] 2) " dB, C/N0 = " (math round [cn0] 2) " dBHz, Eb/N0 = " (math round [ebn0] 2) " dB, BER = " (math round_sig [ber] 3) ", required Eb/N0 = " (math round [req] 2) " dB, margin = " (math round (math calc "[ebn0]-[req]") 2) " dB"

Result

mentdb

FSPL = 205.63 dB, C/N0 = 84.47 dBHz, Eb/N0 = 11.46 dB, BER = 0.0000000606, required Eb/N0 = 10.53 dB, margin = 0.93 dB

Go further

Add the rain attenuation of ITU-R P.618 for the station climate.
Find with num solve the highest data rate that keeps a 3 dB margin.

executed on MentDB

Attitude matrix, Euler angle and inertia tensor

The attitude of a spacecraft is given by yaw, pitch and roll angles. The script builds the direction cosine matrix of the 3-2-1 sequence, checks it is a proper rotation, gets the Euler principal rotation angle from its trace and turns the diagonal body inertia tensor into the reference frame.

Pipeline

01 · rotationsElementary rotationsmath cos→02 · dcmC = R1 R2 R3matrix mul→03 · checkdet C = 1matrix det→04 · anglePrincipal angle from the tracematrix trace→05 · inertiaC^T I Cmatrix diag

Script

admin

#1. Attitude of a spacecraft: yaw 30 deg, pitch 20 deg, roll 10 deg, 3-2-1 sequence;
-> "[y]" (math deg_to_rad 30); -> "[p]" (math deg_to_rad 20); -> "[r]" (math deg_to_rad 10);
-> "[R3]" (concat "[[" (math cos [y]) "," (math sin [y]) ",0],[" (math calc "-sin([y])") "," (math cos [y]) ",0],[0,0,1]]");
-> "[R2]" (concat "[[" (math cos [p]) ",0," (math calc "-sin([p])") "],[0,1,0],[" (math sin [p]) ",0," (math cos [p]) "]]");
-> "[R1]" (concat "[[1,0,0],[0," (math cos [r]) "," (math sin [r]) "],[0," (math calc "-sin([r])") "," (math cos [r]) "]]");
#2. Direction cosine matrix C = R1 R2 R3, proper rotation check det C = 1;
-> "[C]" (matrix mul [R1] (matrix mul [R2] [R3]));
#3. Euler principal angle from the trace;
-> "[phi]" (math rad_to_deg (math acos (math calc (concat "(" (matrix trace [C]) "-1)/2"))));
#4. Body inertia tensor, diagonal, seen in the reference frame: C^T I C;
-> "[Ib]" (matrix diag "[1200, 1500, 900]");
-> "[Ir]" (matrix mul (matrix transpose [C]) (matrix mul [Ib] [C]));
concat "det C = " (math round (matrix det [C]) 12) ", principal angle = " (math round [phi] 3) " deg, I_xy in the reference frame = " (math round (matrix get [Ir] 1 2) 2) " kg m2, trace kept = " (math round (matrix trace [Ir]) 6) " kg m2"

Result

mentdb

det C = 1, principal angle = 35.817 deg, I_xy in the reference frame = -118.8 kg m2, trace kept = 3600 kg m2

Go further

Get the Euler axis with matrix eigen: it is the eigenvector of the eigenvalue 1.
Convert the matrix to a quaternion and compare with the principal angle.