#

---

---

---

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")

---

---

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}")

---

---

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")

---

---

---

If you divide 1 by 998,001 you get all three digit numbers from 000 to 999 in order, with the exception of 998

---

---

---

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)

---

---

Print prime numbers in a range with Python

#!/usr/bin/env python3
import sys
m = int(sys.argv[1])
n = int(sys.argv[2])
primes = [i for i in range(m,n) if all(i%j !=0 for j in range(2,int(i**0.5) + 1))]
print(primes)

---

---

---

---

Solve quadratic equations with impunity using Python (and complex numbers)

#!/usr/bin/env python3
import math,sys
a,b,c=map(float,sys.argv[1:])
d=b*b-4*a*c
if d>=0:x1=(-b+math.sqrt(d))/(2*a);x2=(-b-math.sqrt(d))/(2*a)
else:r=-b/(2*a);i=math.sqrt(-d)/(2*a);x1,x2=complex(r,i),complex(r,-i)
print("The function has","two real roots: {} and {}".format(x1,x2)if d>0 else
"one double root: {}".format(x1)if d==0 else
"two complex roots: {} and {}".format(x1,x2))

---

Python for all your linear algebra homework:

#!/usr/bin/env python3
import sys,ast,numpy as np
a=np.array(ast.literal_eval(sys.argv[1]))
print("Matrix:")
print(np.matrix(a))
print("\nTranspose:")
print(np.matrix(a.T))
print("\nTrace:",a.trace())
print("\nRank:",np.linalg.matrix_rank(a))
if a.shape[0]==a.shape[1]:
    print("\nDeterminant:",round(np.linalg.det(a),6))
    d,e=np.linalg.eig(a)
    print("\nEigenvalues:",np.round(d,6))
    print("\nEigenvectors:")
    print(np.round(e,6))
print("\nPseudo-inverse:")
print(np.matrix(np.linalg.pinv(a)))

 

 

---

Factors

#! /usr/bin/python3
from sys import*
from sympy.ntheory import factorint as f
for i in range(int(argv[1]),int(argv[2])+1):print(i,f(i,multiple=1))

---

Egyptian fractions with Python

#!/usr/bin/env python3
import sys
from sympy import Rational
from sympy.ntheory.egyptian_fraction import egyptian_fraction
r=Rational(*map(int,sys.argv[1].split('/')))
e=egyptian_fraction(r)
print(r,"="," + ".join(f"1/{i}" for i in e))

---

Days between two dates with Python

#!/usr/bin/env python3
from datetime import date;import sys
m1,d1,y1=map(int,sys.argv[1].split('/'));m2,d2,y2=map(int,sys.argv[2].split('/'))
print((date(y2,m2,d2)-date(y1,m1,d1)).days)

---

Here’s a #bash script I saved as “pi” somewhere in $PATH

#!/bin/bash
echo "scale=$1;a(1)*4" | bc -l

I demonstrate it this way:

$ for i in {1..20}; do pi $i; done

2.8 
3.12 
3.140 
3.1412 
3.14156 
3.141592 
3.1415924 
3.14159264 
3.141592652 
3.1415926532 
3.14159265356 
3.141592653588 
3.1415926535896 
3.14159265358976 
3.141592653589792 
3.1415926535897932 
3.14159265358979320 
3.141592653589793236 
3.1415926535897932384 
3.14159265358979323844

---

#! /usr/bin/python3
import math
import sys
a = float(sys.argv[1])
b = float(sys.argv[2])
h = ((a - b) / (a + b))**2
eps = sys.float_info.epsilon
def binom_half(n):
    coef = 1.0
    for k in range(n):
        coef *= (0.5 - k) / (k + 1)
    return coef
perimeter = 0.0
n = 0
term = 1.0  
while abs(term) > eps:
    coef = binom_half(n)
    term = (coef**2) * h**(2*n)
    perimeter += term
    n += 1
perimeter *= math.pi * (a + b)
print(f"Ellipse perimeter: {perimeter:.17f}")
print(f"Series converged after {n} terms")

---

Closing in on a square wave with the Fourier series under Python

#!/usr/bin/env python3
import numpy as n,matplotlib.pyplot as p
x=n.linspace(0,2*n.pi,1000)
def s(x,t):
    k=n.arange(1,2*t,2)
    return 4/n.pi*(n.sin(k[:,None]*x)/k[:,None]).sum(0)
N=12
colors=p.cm.rainbow(n.linspace(0,1,N))
for i,t in enumerate(range(1,N+1)):
    p.plot(x,s(x,t),color=colors[i],label=f'{2*t-1} harmonics')
p.title('Fourier Series Approximation of a Square Wave')
p.grid()
p.legend()
p.show()

---

---

---

---

Iteratively converge on LEO-deliverable payload mass for the Astrodyne Argo Epsilon single-booster first stage (RP-1/LOX) and Delta second stage (LH2/LOX) stack:

#!/usr/bin/env python3
import math
S1_WET = 46357.2
S1_DRY = 3708.6
S1_ISP = 280.0
S2_WET = 19322.5
S2_DRY = 2318.7
S2_ISP = 430.0
TARGET_DV = 7800.0
G0 = 9.80665
def stage_dv(wet, dry, payload, isp):
    m0 = wet + payload
    m1 = dry + payload
    return isp * G0 * math.log(m0 / m1)
def total_dv(payload):
    dv2 = stage_dv(
        S2_WET,
        S2_DRY,
        payload,
        S2_ISP,
    )
    stage2_stack = S2_WET + payload

    dv1 = stage_dv(
        S1_WET,
        S1_DRY,
        stage2_stack,
        S1_ISP,
    )
    return dv1, dv2, dv1 + dv2
low = 0.0
high = 100000.0
for _ in range(60):
    payload = (low + high) / 2
    dv1, dv2, total = total_dv(payload)
    if total > TARGET_DV:
        low = payload
    else:
        high = payload
payload = (low + high) / 2
dv1, dv2, total = total_dv(payload)
print(f"Target delta-v : {TARGET_DV:10.1f} m/s")
print(f"Payload        : {payload:10.1f} kg")
print(f"Stage 1 Δv     : {dv1:10.1f} m/s")
print(f"Stage 2 Δv     : {dv2:10.1f} m/s")
print(f"Total Δv       : {total:10.1f} m/s")

---

---

---

---

---

π calculation using a Machin-like arctangent formula

#!/usr/bin/python3
from decimal import Decimal, getcontext
import sys
digits = int(sys.argv[1]) if len(sys.argv) > 1 else 50
getcontext().prec = digits + 5  # a few extra digits to avoid rounding errors

def arctan(x):
    """Compute arctan(1/x) using the Taylor series."""
    x = Decimal(x)
    x2 = x * x
    term = Decimal(1) / x
    total = term
    n = 1
    sign = -1
    while True:
        term = term / x2
        delta = term / (2 * n + 1)
        if delta == 0:
            break
        total += sign * delta
        sign *= -1
        n += 1
    return total

pi = 16 * arctan(5) - 4 * arctan(239)
pi = +pi  # unary plus applies current context precision
print(f"Pi to {digits} digits:\n{str(pi)[:digits+2]}")

---

---

---

---

---

---

Generating the Weierstrass Monster with #Python

This is pathology of calculus that is continuous everywhere but differentiable nowhere. Henri Poincaré condemned it as “an outrage against common sense.” Charles Hermite called it a “deplorable evil.” Muahahahaha!

#!/usr/bin/env python3
import numpy as n,matplotlib.pyplot as p
a=.5;b=11;n_terms=50
x=n.linspace(-2,2,20000)
f=sum(a**k*n.cos(b**k*n.pi*x) for k in range(n_terms))
p.plot(x,f,'b',lw=.1)
p.title('Weierstrass Function')
p.grid()
p.show()

teresita:Desktop$ seq 3 | paste -sd* | bc

6

teresita:Desktop$ seq 30 | paste -sd* | bc

265252859812191058636308480000000

teresita:Desktop$ seq 300 | paste -sd* | bc

306057512216440636035370461297268629388588804173576999416776741259476533176716867465515291422477573349939147888701726368864263907759003154226842927906974559841225476930271954604008012215776252176854255965356903506788725264321896264299365204576448830388909753943489625436053225980776521270822437639449120128678675368305712293681943649956460498166450227716500185176546469340112226034729724066333258583506870150169794168850353752137554910289126407157154830282284937952636580145235233156936482233436799254594095276820608062232812387383880817049600000000000000000000000000000000000000000000000000000000000000000000000000

---

---

---

---

---

Calculate the nth prime with Python

#!/usr/bin/env python3
import sys
import time; max=int(sys.argv[1]); n=1; p=1
t0 = time.perf_counter()
while n<=max:
    f=0; j=2; s = int(p**0.5)
    while f < 1:
        if j >= s:
            f=2
        if p % j == 0:
            f=1
        j+=1
    if f != 1:
        n+=1;
        #print(n,p);
    p+=1
t1 = time.perf_counter()
print(f"n = {max}")
print(f"nth prime = {p-1}")
print(f"Elapsed time = {t1-t0:.6f} seconds")

---

---

---

---

Calculate phi with bc and bash

#!/usr/bin/env bash
a=0 b=1
for((i=2;i<=23;i++));{
    c=$((a+b))
    printf "Fibonacci(%d) / Fibonacci(%d) = %.20f\n" $i $((i-1)) "$(bc -l<<<"$c/$b")"
    a=$b b=$c
}

---

---

Solve for the intersection of two lines from four points with Python

#!/usr/bin/env python3
import sys,numpy as n
p,q,r,s=map(n.array,((float(sys.argv[i]),float(sys.argv[i+1])) for i in (1,3,5,7)))
M=n.c_[q-p,s-r]
d=n.linalg.det(M)
if abs(d)<1e-10:
    print("Coincident" if abs(n.cross(q-p,r-p))<1e-10 else "Parallel")
else:
    t,_=n.linalg.solve(M,r-p)
    print(p+t*(q-p))

---

---

---

---

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")

---

---

---

---

---

---

---

---

Plotting a trig function and its derivative with Python

#!/usr/bin/env python3
import numpy as n,matplotlib.pyplot as p
x=n.linspace(-4,4,400)
f=n.arctan(x**2);df=2*x/(x**4+1)
p.plot(x,f,'b',label='f(x)=arctan(x²)')
p.plot(x,df,'r',label="f'(x)=2x/(x⁴+1)")
p.axhline(0,c='k',lw=.5);p.axvline(0,c='k',lw=.5)
p.ylim(-2,2);p.xlabel('x');p.ylabel('y')
p.title("f(x) and f'(x)")
p.grid(ls=':',alpha=.7);p.legend();p.show()

---

---

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")

---

---

---

---

---

Calculate area and perimeter of any triangle

#!/usr/bin/python3
import sys, math
a, b, c = map(float, sys.argv[1:])
if a + b <= c or a + c <= b or b + c <= a:
    print("Error: not a valid triangle")
    sys.exit(1)
peri = a + b + c
area = math.sqrt((a + b + c) *
                 (a + b - c) *
                 (a - b + c) *
                 (-a + b + c)) / 4
print("Perimeter =", peri)
print("Area =", area)

---

---

#! /usr/bin/env python3
import matplotlib.pyplot as plt
import numpy as np
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
r = 10
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = r * np.outer(np.cos(u), np.sin(v))
y = r * np.outer(np.sin(u), np.sin(v))
z = r * np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x, y, z, color='linen', alpha=0.5)
theta = np.linspace(0, 2 * np.pi, 100)
z = np.zeros(100)
x = r * np.sin(theta)
y = r * np.cos(theta)
ax.plot(x, y, z, color='black', alpha=0.75)
ax.plot(z, x, y, color='black', alpha=0.75)
zeros = np.zeros(1000)
line = np.linspace(-10,10,1000)
ax.plot(line, zeros, zeros, color='black', alpha=0.75)
ax.plot(zeros, line, zeros, color='black', alpha=0.75)
ax.plot(zeros, zeros, line, color='black', alpha=0.75)
plt.show()

---

---

---

The Einstein field equations rendered in Libreoffice Math

R Rsub{%mu%nu}-1 over 2 g rsub{%mu%nu}R+%LAMBDA g rsub{%mu %nu}~=~{8%pi`G} over c^4 T rsub{%mu%nu}

---

---

---

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")

---

---

---

---

---

---