Terminal Cruise

Interplanetary mission velocity change and time of flight using Python (example used is outbound leg to Mars):

#!/usr/bin/env python3
import sys
from math import sqrt, pi, radians, sin
muP, muD, r1, rp1, muA, r2, rp2, inc = map(float, sys.argv[1:])
a = (r1 + r2) / 2.0
tof = pi * sqrt(a**3 / muP)
v1 = sqrt(muP / r1)
v2 = sqrt(muP / r2)
vt1 = sqrt(muP * (2.0 / r1 - 1.0 / a))
vt2 = sqrt(muP * (2.0 / r2 - 1.0 / a))
vinf1 = abs(vt1 - v1)
vinf2 = abs(v2 - vt2)
vc1 = sqrt(muD / rp1)
vp1 = sqrt(vinf1**2 + 2.0 * muD / rp1)
dv1 = vp1 - vc1
vc2 = sqrt(muA / rp2)
vp2 = sqrt(vinf2**2 + 2.0 * muA / rp2)
dv2 = vp2 - vc2
plane_speed = min(vt1, vt2)
plane = 2.0 * plane_speed * sin(radians(inc) / 2.0)
print(f"Transfer time   {tof/86400:12.3f}    days")
print(f"Injection Δv       {dv1:12.6f} km/s")
print(f"Capture Δv         {dv2:12.6f} km/s")
print()
print(f"Best case total Δv {dv1+dv2:12.6f} km/s")
print(f"Plane change       {plane  :12.6f} km/s")
print(f"Worst case total Δv{dv1+dv2+plane:12.6f} km/s")

Calculate a great circle route in kilometers:

#!/usr/bin/env python3
import sys
from math import radians, degrees, sin, cos, sqrt, asin
lat1, lon1, lat2, lon2 = map(float, sys.argv[1:5])
dLat = radians(lat2 - lat1)
dLon = radians(lon2 - lon1)
lat1 = radians(lat1)
lat2 = radians(lat2)
a = sin(dLat / 2)**2 + cos(lat1) * cos(lat2) * sin(dLon / 2)**2
c = 2 * asin(sqrt(a))
print(6372.8 * c)

Nashville International Airport (BNA)

  N 36°7.2',   W 86°40.2'     (36.12,   -86.67)          

Los Angeles International Airport (LAX)

  N 33°56.4',  W 118°24.0'    (33.94,  -118.40)

Find barycentric L1 and L2 Lagrange points from normalized mass ratio

#!/usr/bin/env python3
import sys
from scipy.optimize import brentq
R = float(sys.argv[1])
m1 = R / (R + 1.0)
m2 = 1.0 / (R + 1.0)
x1 = -m2
x2 =  m1
def f(x):
    r1 = abs(x - x1)
    r2 = abs(x - x2)
    return (
        x
        - m1 * (x - x1) / r1**3
        - m2 * (x - x2) / r2**3
    )
eps = 1e-12
L1 = brentq(f, x1 + eps, x2 - eps)
L2 = brentq(f, x2 + eps, x2 + 2.0)
r1 = x2 - L1
r2 = L2 - x2
print(f"Primary   : {x1:.12f}")
print(f"Secondary : {x2:.12f}")
print(f"L1        : {L1:.12f}")
print(f"L2        : {L2:.12f}")
print(f"L1 offset : {r1:.12f}")
print(f"L2 offset : {r2:.12f}")

Galactic Cartesian coordinates from RA and DEC (Alpha Centauri used for the example, and compared to a run with Celestial coordinates)

#!/usr/bin/env python3
import sys
from astropy import units as u
from astropy.coordinates import SkyCoord
c = SkyCoord(
    " ".join(sys.argv[1:7]),
    unit=(u.hourangle, u.deg),
    distance=float(sys.argv[7]))
g = c.galactic
print(g.cartesian)

Python script for calculating the separation in AU between two stars in a binary system with known masses and period of revolution (the example uses Kruger 60 A/B):

#!/usr/bin/env python3
import math
import sys
m1 = float(sys.argv[1])
m2 = float(sys.argv[2])
P = float(sys.argv[3])
a = ((m1 + m2) * P**2) ** (1/3)
print(f"Total mass        : {m1 + m2:.4f} M☉")
print(f"Orbital period    : {P:.2f} years")
print(f"Separation    : {a:.3f} AU")

With Python calculate the “comfort zone” of a star system based on visual magnitude and distance. Example is for Ran aka Epsilon Eridani (at Sol the comfort zone happens to be at 1 AU, what a coincidence!)

#!/usr/bin/env python3
import math
import sys
SUN_ABS_MAG = 4.83
LY_PER_PC = 3.26156
m = float(sys.argv[1])      
d_ly = float(sys.argv[2])
d_pc = d_ly / LY_PER_PC
M = m - 5 * (math.log10(d_pc) - 1)
L = 10 ** ((SUN_ABS_MAG - M) / 2.5)
R = math.sqrt(L)
print(f"Distance           : {d_ly:.3f} ly ({d_pc:.3f} pc)")
print(f"Absolute magnitude : {M:.3f}")
print(f"Luminosity         : {L:.3f} L☉")
print(f"Goldilocks radius  : {R:.3f} AU")

Distance between stars with linear algebra in #Python (example used is the milk run from Rigilkent to Proxima):

#!/usr/bin/env python3
import sys
import numpy as np
p1 = np.array(sys.argv[1:4], dtype=float)
p2 = np.array(sys.argv[4:7], dtype=float)
print(np.linalg.norm(p2 - p1))

Tsiolkovsky rocket equation

#!/usr/bin/env python3
import math
import sys
wet = float(sys.argv[1])
dry = float(sys.argv[2])
isp = float(sys.argv[3])
g0 = 9.80665
delta_v = isp * g0 * math.log(wet / dry)
print(f"{delta_v:.1f} m/s")

Use Python to reverse stdin by words (rev will do this by characters)

#!/usr/bin/env python3
import sys
for line in sys.stdin:
    print(' '.join(line.split()[::-1]))

Python replacement for the uniq command.

#!/usr/bin/env python3
import sys
seen = set()
for line in sys.stdin:
    if line not in seen:
        seen.add(line)
        sys.stdout.write(line)

This is not equivalent to sort -u because it preserves the original order, and it’s not equivalent to uniq because it removes duplicates even when they’re not adjacent.