Public
Edited
Apr 4, 2023
G = math.gravitationConstant // Gravitational constant
mbody = math.unit(5.9724e24,'kg') // Mass of Earth
mu = math.multiply(G, mbody) // Standard gravitational parameter
g0 = math.gravity // Standard gravity: used for calculating propellant consumption (dmdt)
r0 = math.unit(6371, "km") // Mean radius of Earth
t0 = math.unit(0,"s") // Simulation start
dt = math.unit(0.5, "s") // Simulation timestep
tfinal = math.unit(149.5, "s") // Simulation duration
isp_sea = math.unit(282, "s") // Specific impulse (at sea level)
isp_vac = math.unit(311, "s") // Specific impulse (in vacuum)
gamma0 = math.unit(89.99970, "deg") // Initial pitch angle (90 deg is vertical)
v0 = math.unit(0.9, "m/s") // Initial velocity (must be non-zero because ODE is ill-conditioned)
phi0 = math.unit(0, "deg") // Initial orbital reference angle
m1 = math.unit(433100, "kg") // First stage mass
m2 = math.unit(111500, "kg") // Second stage mass
m3 = math.unit(1700, "kg") // Third stage / fairing mass
mp = math.unit(5000, "kg") // Payload mass
m0 = m1+m2+m3+mp // Initial mass of rocket
dm = math.unit(2750, "kg/s") // Mass flow rate
rocket_radius = math.unit(3.66, 'm')
A = math.multiply(math.pow(rocket_radius, 2), math.pi) // Area of the rocket
dragCoef = 0.2 // Drag coefficient
function odeRK4(odeFunc, y0, t0, tRange, h) {
// Define an array to hold the results
var results = [];
// Initialize the dependent variable and the independent variable
var y = y0;
var t = t0;
// Loop over the time range and compute the solution
for (var i = 0; t < tRange[1]; i++) {
// Compute the four intermediate slopes
var k1 = h * odeFunc(y, t);
var k2 = h * odeFunc(y + k1/2, t + h/2);
var k3 = h * odeFunc(y + k2/2, t + h/2);
var k4 = h * odeFunc(y + k3, t + h);
// Compute the new value of the dependent variable
y = y + 1/6 * (k1 + 2*k2 + 2*k3 + k4);
// Compute the new value of the independent variable
t = t + h;
// Add the new values to the results array
results.push({t: t, y: y});
}
// Return the results
return results;
}
//Define the equations of motion. We just thrust into current direction of motion, e.g. making a gravity turn.
gravity = (r) => math.divide(mu, math.dotPow(r,2))
angVel = (r, v, gamma) => math.multiply(math.divide(v,r), math.cos(gamma), math.rad) // Angular velocity of rocket around moon
density = (r) => math.divide(
math.multiply(math.unit(1.2250, "kg/m^3"), math.exp(-g0 * (r - r0))),
math.unit(83246.8, "m^2/s^2")) // Assume constant temperature
drag = (r, v) => math.multiply(1/2, density(r), math.square(v, 2), A , dragCoef)
# pressure ~ density for constant temperature
isp(r) = isp_vac + (isp_sea - isp_vac) * density(r)/density(r0)
thrust(r) = g0 * isp(r) * dm
# It is important to maintain the same argument order for each of these functions.
drdt(r, v, m, phi, gamma, t) = v sin(gamma)
dvdt(r, v, m, phi, gamma, t) = - gravity(r) * sin(gamma) + (thrust(r) - drag(r, v)) / m
dmdt(r, v, m, phi, gamma, t) = - dm
dphidt(r, v, m, phi, gamma, t) = angVel(r, v, gamma)
dgammadt(r, v, m, phi, gamma, t) = angVel(r, v, gamma) - gravity(r) * cos(gamma) / v * rad
dtdt(r, v, m, phi, gamma, t) = 1
// Remember to maintain the same variable order in the call to ndsolve.
result_stage1 = ndsolve([drdt, dvdt, dmdt, dphidt, dgammadt, dtdt], [r0, v0, m0, phi0, gamma0, t0], dt, tfinal)
//# Reset initial conditions for interstage flight
dm = 0 kg/s
tfinal = 10 s
x = flatten(result_stage1[end,:])
x[3] = m2+m3+mp # New mass after stage seperation
result_interstage = ndsolve([drdt, dvdt, dmdt, dphidt, dgammadt, dtdt], x, dt, tfinal)
# Reset initial conditions for stage 2 flight
dm = 270.8 kg/s
isp_vac = 348 s
tfinal = 350 s
x = flatten(result_interstage[end,:])
result_stage2 = ndsolve([drdt, dvdt, dmdt, dphidt, dgammadt, dtdt], x, dt, tfinal)
# Reset initial conditions for unpowered flight
dm = 0 kg / s
tfinal = 900 s
dt = 10 s
x = flatten(result_stage2[end,:])
result_unpowered1 = ndsolve([drdt, dvdt, dmdt, dphidt, dgammadt, dtdt], x, dt, tfinal)
# Reset initial conditions for final orbit insertion
dm = 270.8 kg / s
tfinal = 39 s
dt = 0.5 s
x = flatten(result_unpowered1[end,:])
result_insertion = ndsolve([drdt, dvdt, dmdt, dphidt, dgammadt, dtdt], x, dt, tfinal)
# Reset initial conditions for unpowered flight
dm = 0 kg / s
tfinal = 250 s
dt = 10 s
x = flatten(result_insertion[end,:])
result_unpowered2 = ndsolve([drdt, dvdt, dmdt, dphidt, dgammadt, dtdt], x, dt, tfinal)`
math = import("https://cdn.skypack.dev/mathjs@11.8")
Purpose-built for displays of data
Observable is your go-to platform for exploring data and creating expressive data visualizations. Use reactive JavaScript notebooks for prototyping and a collaborative canvas for visual data exploration and dashboard creation.
