calculators/rocket/payload/Atmosphere.js

		+class Layer {
		+ static TEMPERATURES = {
		+ TROPOSPHERE: [{
		+ range: [0, 11],
		+ temperature: (altitude) => 288.15 - 6.5 * altitude
		+ }],
		+ STRATOSPHERE: [{
		+ range: [11, 20],
		+ temperature: (altitude) => 216.65
		+ }, {
		+ range: [20, 32],
		+ temperature: (altitude) => 196.65 + altitude
		+ }, {
		+ range: [32, 47],
		+ temperature: (altitude) => 139.05 + 2.8 * altitude
		+ }],
		+ MESOSPHERE: [{
		+ range: [47, 51],
		+ temperature: (altitude) => 270.65
		+ }, {
		+ range: [51, 71],
		+ temperature: (altitude) => 413.45 - 2.8 * altitude
		+ }, {
		+ range: [71, 84.852],
		+ temperature: (altitude) => 356.65 - 2.0 * altitude
		+ }],
		+ THERMOSPHERE: [{
		+ range: [84.852, 91],
		+ temperature: (altitude) => 186.8673
		+ }, {
		+ range: [91, 110],
		+ temperature: (altitude) => 263.1905 - 76.3232 * Math.sqrt(1 - Math.pow((altitude - 91) / -19.9429, 2))
		+ }, {
		+ range: [110, 120],
		+ temperature: (altitude) => 240 + 12 * (altitude - 110)
		+ }, {
		+ range: [120, 1000],
		+ temperature: (altitude) => {
		+ const ξ = (altitude - 120) * (6356.766 + 120) / (6356.766 + altitude);
		+ return 1000 - 640 * Math.exp(-0.01875 * ξ);
		+ }
		+ }, {
		+ range: [1000, Infinity],
		+ temperature: (altitude) => 1000
		+ }]
		+ };
		+
		+ constructor(name, range, pressure, density) {
		+ this.name = name;
		+ this.range = range;
		+ this.pressure = pressure;
		+ this.density = density;
		+ }
		+
		+ within(altitude) {
		+ return altitude >= this.range[0] && altitude <= this.range[1];
		+ }
		+
		+ temperature(altitude) {
		+ let result;
		+ const temperatures = Layer.TEMPERATURES[this.name];
		+
		+ for (let t of temperatures) {
		+ if (altitude >= t.range[0] && altitude <= t.range[1]) {
		+ result = t.temperature(altitude);
		+ break;
		+ }
		+ }
		+
		+ console.assert(result !== null && !isNaN(result));
		+ return result;
		+ }
		+}
		+
		+/**
		+ * See http://www.braeunig.us/space/atmmodel.htm
		+ */
		+class Atmosphere {
		+ static EQUATORIAL_RADIUS = 6378137.0; // m
		+ static FLATTENING = 1.0 / 298.25722356;
		+ static SPECIFIC_GAS_CONSTANT = 287.053; // J/kg-K
		+
		+ static LAYERS = [
		+ new Layer(
		+ "TROPOSPHERE",
		+ [0, 11],
		+ (altitude) => 101325.0 * Math.pow((288.15 / (288.15 - 6.5 * altitude)), (34.1632 / -6.5)),
		+ (pressure, gas, temperature) => pressure / (gas * temperature)
		+ ),
		+ new Layer(
		+ "STRATOSPHERE",
		+ [11, 20],
		+ (altitude) => 22632.06 * Math.exp(-34.1632 * (altitude - 11) / 216.65),
		+ (pressure, gas, temperature) => pressure / (gas * temperature)
		+ ),
		+ new Layer(
		+ "STRATOSPHERE",
		+ [20, 32],
		+ (altitude) => 5474.889 * Math.pow(216.65 / (216.65 + (altitude - 20)), 34.1632),
		+ (pressure, gas, temperature) => pressure / (gas * temperature)
		+ ),
		+ new Layer(
		+ "STRATOSPHERE",
		+ [32, 47],
		+ (altitude) => 868.0187 * Math.pow(228.65 / (228.65 + 2.8 * (altitude - 32)), 34.1632 / 2.8),
		+ (pressure, gas, temperature) => pressure / (gas * temperature)
		+ ),
		+ new Layer(
		+ "MESOSPHERE",
		+ [47, 51],
		+ (altitude) => 110.9063 * Math.exp(-34.1632 * (altitude - 47) / 270.65),
		+ (pressure, gas, temperature) => pressure / (gas * temperature)
		+ ),
		+ new Layer(
		+ "MESOSPHERE",
		+ [51, 71],
		+ (altitude) => 66.93887 * Math.pow(270.65 / (270.65 - 2.8 * (altitude - 51)), 34.1632 / -2.8),
		+ (pressure, gas, temperature) => pressure / (gas * temperature)
		+ ),
		+ new Layer(
		+ "MESOSPHERE",
		+ [71, 84.852],
		+ (altitude) => 3.956420 * Math.pow(214.65 / (214.65 - 2 * (altitude - 71)), 34.1632 / -2),
		+ (pressure, gas, temperature) => pressure / (gas * temperature)
		+ ),
		+ new Layer(
		+ "THERMOSPHERE",
		+ [84.852, Infinity],
		+ (altitude) =>
		+ Math.exp(-0.0000000422012 *
		+ Math.pow(altitude, 5) + 0.0000213489 *
		+ Math.pow(altitude, 4) - 0.00426388 *
		+ Math.pow(altitude, 3) + 0.421404 *
		+ Math.pow(altitude, 2) - 20.8270 * altitude + 416.225),
		+ (altitude) => Math.exp(
		+ 0.000000075691 *
		+ Math.pow(altitude, 5) - 0.0000376113 *
		+ Math.pow(altitude, 4) + 0.0074765 *
		+ Math.pow(altitude, 3) - 0.743012 *
		+ Math.pow(altitude, 2) + 36.7280 * altitude - 729.346)
		+ )
		+ ];
		+
		+ constructor(latitude) {
		+ this.latitude = latitude;
		+ }
		+
		+ radius(latitude) {
		+ const phi = latitude * (Math.PI / 180);
		+ return Atmosphere.EQUATORIAL_RADIUS * (1 - Atmosphere.FLATTENING * Math.sin(phi) ** 2);
		+ }
		+
		+ geopotential(altitude) {
		+ const radius = this.radius(this.latitude);
		+ return (altitude * radius) / (radius + altitude);
		+ }
		+
		+ layerValue(altitude, callback) {
		+ let result;
		+
		+ const gpAltitude = this.geopotential(altitude / 1000);
		+
		+ for (let layer of Atmosphere.LAYERS) {
		+ if (layer.within(gpAltitude)) {
		+ result = callback(layer, gpAltitude);
		+ break;
		+ }
		+ }
		+
		+ this.invariant(result);
		+
		+ return result;
		+ }
		+
		+ temperature(altitude) {
		+ return this.layerValue(altitude, (layer, gpAltitude) => layer.temperature(gpAltitude));
		+ }
		+
		+ airPressure(altitude) {
		+ return this.layerValue(altitude, (layer, gpAltitude) => layer.pressure(gpAltitude));
		+ }
		+
		+ airDensity(altitude) {
		+ return this.layerValue(altitude, (layer, gpAltitude) => {
		+ const temperature = layer.temperature(gpAltitude);
		+ const pressure = layer.pressure(gpAltitude);
		+
		+ return layer.density(pressure, this.constructor.SPECIFIC_GAS_CONSTANT, temperature);
		+ });
		+ }
		+
		+ invariant(result) {
		+ console.assert(result !== null && !isNaN(result));
		+ }
		+}
		+
		+export default Atmosphere;

calculators/rocket/payload/Earth.js

		+import Atmosphere from './Atmosphere.js';
		+
		class Earth {
		static KELVIN = 273.15;
		- static STATIC_PRESSURE = 101325;
		- static TEMPERATURE_LAPSE_RATE = -0.0065;
		- static GAS_CONSTANT = 8.31432;
		- static GRAVITY = 9.80665;
		- static MOLAR_MASS_AIR = 0.0289644;
		- static SPECIFIC_GAS_CONSTANT_DRY_AIR = 287.85;
		- static SPECIFIC_GAS_CONSTANT_WATER_VAPOUR = 461.495;
		- static RADIUS_EQUATOR = 6378137.0;
		- static RADIUS_POLAR = 6356752.0;
		static MASS = 5.9722e24;
		static SURFACE_GRAVITY = 9.780327;
		static SIDEREAL_TIME = 86164.0905;
		static SPEED_SOUND_0C = 331.3;
		static GRAVITATIONAL_CONST = 6.67430e-11;
		static STANDARD_GRAVITY = Earth.GRAVITATIONAL_CONST * Earth.MASS;
		-
		- constructor(temperature, humidity) {
		- this.surface = temperature;
		- this.humidity = humidity / 100;
		- }

		- kelvins(celsius) {
		- return celsius + Earth.KELVIN;
		+ constructor(latitude) {
		+ this.atmosphere = new Atmosphere(latitude);
		}

		radians(degrees) {
		return degrees * (Math.PI / 180);
		}

		radius(latitude) {
		- return (
		- Earth.RADIUS_EQUATOR *
		- (1 -
		- ((Earth.RADIUS_EQUATOR - Earth.RADIUS_POLAR) /
		- Earth.RADIUS_EQUATOR) *
		- Math.pow(Math.sin(this.radians(latitude)), 2))
		- );
		+ return this.atmosphere.radius(latitude);
		}



		soundSpeed(altitude) {
		- const temperature = this.surface - (6 * altitude) / 1000;
		-
		- return Earth.SPEED_SOUND_0C * Math.sqrt(1 + temperature / Earth.KELVIN);
		+ const temperature = this.atmosphere.temperature(altitude);
		+ return Earth.SPEED_SOUND_0C * Math.sqrt((1 + temperature) / Earth.KELVIN);
		}



		airPressure(altitude) {
		- const T = this.kelvins(this.surface);
		- const P = Earth.STATIC_PRESSURE;
		- const g = Earth.GRAVITY;
		- const M = Earth.MOLAR_MASS_AIR;
		- const R = Earth.GAS_CONSTANT;
		- const L = Earth.TEMPERATURE_LAPSE_RATE;
		-
		- let pressure;
		-
		- if (altitude <= 11000) {
		- pressure = P * Math.pow((1 + (L / T) * altitude), -(g * M) / (R * L));
		- } else {
		- const T1 = 216.65;
		- const P1 = P * Math.pow((1 + (L / T) * 11000), -(g * M) / (R * L));
		-
		- pressure = P1 * Math.exp(-(g * M) * (altitude - 11000) / (R * T1));
		- }
		-
		- return pressure;
		- }
		-
		- saturationVaporPressure(temperature) {
		- const BASE_PRESSURE = 6.1078;
		- const COEFFICIENTS = [
		- 0.99999683,
		- -0.90826951e-2,
		- 0.78736169e-4,
		- -0.61117958e-6,
		- 0.43884187e-8,
		- -0.29883885e-10,
		- 0.21874425e-12,
		- -0.17892321e-14,
		- 0.11112018e-16,
		- -0.30994571e-19
		- ];
		-
		- let polynomial = 0;
		-
		- for (let i = 0; i < COEFFICIENTS.length; i++) {
		- polynomial += COEFFICIENTS[i] * Math.pow(temperature, i);
		- }
		-
		- return (BASE_PRESSURE / Math.pow(polynomial, 8)) * 100;
		- }
		-
		- waterVaporPressure(humidity, saturationVaporPressure) {
		- return humidity * saturationVaporPressure;
		- }
		-
		- dryAirPressure(totalPressure, waterVaporPressure) {
		- return totalPressure - waterVaporPressure;
		+ return this.atmosphere.airPressure(altitude);
		}

		airDensity(altitude) {
		- const pressure = this.airPressure(altitude);
		- const kelvins = this.kelvins(this.surface);
		- const T_DRY_AIR = Earth.SPECIFIC_GAS_CONSTANT_DRY_AIR * kelvins;
		- const T_WET_AIR = Earth.SPECIFIC_GAS_CONSTANT_WATER_VAPOUR * kelvins;
		- const vapourSaturated = this.saturationVaporPressure(this.surface);
		- const vapourWater = this.waterVaporPressure(this.humidity, vapourSaturated);
		- const dryAir = this.dryAirPressure(pressure, vapourWater);
		- const airDensity = dryAir / T_DRY_AIR + vapourWater / T_WET_AIR;
		+ const density = this.atmosphere.airDensity(altitude);

		- return Math.max(0, airDensity);
		+ return Math.max(0, density);
		}
		}

calculators/rocket/payload/Rocket.js

		console.assert(specificImpulse > 1);

		+ this.wet = wet;
		this.dry = wet * 0.1;
		- this.mass = wet;
		this.targetMass = this.dry + payload;
		this.radius = diameter / 2.0;
		this.specificImpulse = specificImpulse;
		this.totalThrust = 0.0;
		this.targetAltitude = 0.0;

		- const crossSectionArea = Math.PI * Math.pow(this.radius, 2);
		+ const crossSectionArea = Math.PI * this.radius ** 2;
		this.ballisticCoefficient = dragCoefficient * crossSectionArea;

		this.flightRecorder = false;
		this.launched = false;
		this.relocated = false;
		this.world = false;

		- this.massFlowRate = NaN;
		+ this.massFlowRate = 5.0;
		}


		this.maxAcceleration = 5.0 * LOCAL_GRAVITY;
		this.rotationalSpeed = this.planet.rotationalSpeed(latitude, altitude);
		- this.exhaustVelocity = this.specificImpulse * LOCAL_GRAVITY;
		+ this.vExhaust = this.specificImpulse * LOCAL_GRAVITY;
		this.relocated = true;
		}

		launch(hVelocity, vVelocity) {
		console.assert(!this.launched);
		+
		+ this.launched = true;

		this.hVelocity = hVelocity + this.rotationalSpeed;
		this.vVelocity = vVelocity;
		- this.launched = true;
		+
		+ const dragForce = this.drag();
		+ const absDragForce = Math.abs(dragForce[1]);
		+
		+ this.massFlowRate = absDragForce / this.vExhaust * 1.001;
		+
		+ const vOrb = this.planet.orbitalSpeed(this.latitude, this.altitude);
		+ const vReq = vOrb - this.rotationalSpeed;
		+ const vRocket = this.vExhaust * Math.log(this.wet / this.targetMass);
		+
		+ console.info(`orbital and horizontal = ${vOrb} ${this.hVelocity}`);
		+ console.info(`dry + payload = ${this.targetMass}`);
		+ console.info(`Δv = ${vReq}; Δv rocket = ${vRocket}`);
		+ console.info(`mass flow rate = ${this.massFlowRate}`);
		+ console.info(`drag force = ${dragForce}`);
		+ console.info(`hVelocity = ${hVelocity}`);
		+ console.info(`vVelocity = ${vVelocity}`);
		}



		const rho = this.planet.airDensity(this.altitude);
		- const netHorizontalVelocity = this.hVelocity - this.rotationalSpeed;
		- const totalSpeed = Math.sqrt(Math.pow(netHorizontalVelocity, 2) + Math.pow(this.vVelocity, 2));
		- const magnitude = 0.5 * rho * this.ballisticCoefficient * Math.pow(totalSpeed, 2);
		- const forceRatio = [-netHorizontalVelocity / totalSpeed, -this.vVelocity / totalSpeed];
		+ const hVelocityNet = this.hVelocity - this.rotationalSpeed;
		+ const totalSpeed = Math.sqrt(hVelocityNet 2 + this.vVelocity 2);
		+ const magnitude = 0.5 * rho * this.ballisticCoefficient * totalSpeed ** 2;
		+ const forceRatio = [-hVelocityNet / totalSpeed, -this.vVelocity / totalSpeed];
		+
		+ console.info(`rho = ${rho}`);
		+ console.info(`hV (net) = ${hVelocityNet}`);
		+ console.info(`vV = ${this.vVelocity}`);
		+ console.info(`totalSpeed = ${totalSpeed}`);
		+ console.info(`bC = ${this.ballisticCoefficient}`);
		+ console.info(`magnitude = ${magnitude}`);
		+ console.info(`ratio = [${forceRatio}]`);

		return [forceRatio[0] * magnitude, forceRatio[1] * magnitude];
		}

		flying() {
		- console.assert(this.mass > 0);
		+ console.assert(this.wet > 0);
		console.assert(this.targetMass > 0);
		console.assert(this.targetAltitude > 0);

		return (
		- this.mass > this.targetMass &&
		- this.mass > this.dry &&
		+ this.wet > this.targetMass &&
		+ this.wet > this.dry &&
		this.altitude < this.targetAltitude
		);
		}

		fly(time) {
		+ console.assert(this.launched);
		console.assert(time > 0);

		const dragForce = this.drag();
		-
		- if (isNaN(this.massFlowRate)) {
		- this.massFlowRate = Math.abs(dragForce[1]) / this.exhaustVelocity * 1.001;
		- }
		-
		- console.debug(`mass flow rate = ${this.massFlowRate}`);
		+ const absDragForce = Math.abs(dragForce[1]);

		- this.totalThrust = this.exhaustVelocity * this.massFlowRate;
		+ this.totalThrust = this.vExhaust * this.massFlowRate;

		const thrusting = this.totalThrust > 0;
		const vThrust = thrusting ? -dragForce[1] : 0.0;
		- const hThrust = thrusting ? Math.sqrt(Math.pow(this.totalThrust, 2) - Math.pow(vThrust, 2)) : 0.0;
		+ const hThrust = thrusting ? Math.sqrt(this.totalThrust 2 - vThrust 2) : 0.0;

		- console.assert(this.mass > 0);
		+ console.assert(this.wet > 0);

		- this.hAcceleration = dragForce[0] / this.mass + hThrust / this.mass;
		- this.vAcceleration = dragForce[1] / this.mass + this.planet.gravity(this.latitude, this.altitude) + vThrust / this.mass;
		+ this.hAcceleration = dragForce[0] / this.wet + hThrust / this.wet;
		+ this.vAcceleration = dragForce[1] / this.wet + this.planet.gravity(this.latitude, this.altitude) + vThrust / this.wet;

		// Records the initial values and all subsequent deltas.
		this.record(time);

		this.hVelocity += this.hAcceleration * time;
		this.vVelocity += this.vAcceleration * time;

		- this.azimuth += this.hVelocity * time + 0.5 * this.hAcceleration * time * time;
		- this.altitude += this.vVelocity * time + 0.5 * this.vAcceleration * time * time;
		+ this.azimuth += this.hVelocity * time + 0.5 * this.hAcceleration * time ** 2;
		+ this.altitude += this.vVelocity * time + 0.5 * this.vAcceleration * time ** 2;

		- this.mass -= this.massFlowRate * time;
		+ this.wet -= this.massFlowRate * time;
		+
		+ console.log(`Time: ${time}, Altitude: ${this.altitude}, hVelocity: ${this.hVelocity}, vVelocity: ${this.vVelocity}`);
		+ console.log(`Drag: ${dragForce}, hAcc: ${this.hAcceleration}, vAcc: ${this.vAcceleration}, Mass: ${this.wet}, Mass Flow Rate: ${this.massFlowRate}`);

		if (this.hVelocity < this.targetVelocity) {
		const accelerationMagnitude = Math.sqrt(
		- Math.pow(this.hAcceleration, 2) +
		- Math.pow(this.vAcceleration, 2)
		+ this.hAcceleration 2 + this.vAcceleration 2
		);
		-
		- const excessiveAcceleration = accelerationMagnitude > this.maxAcceleration;
		- const reducibleFlowRate = this.exhaustVelocity * this.massFlowRate * 0.98 > Math.abs(dragForce[1]);

		- if (excessiveAcceleration && reducibleFlowRate) {
		+ if (accelerationMagnitude > this.maxAcceleration &&
		+ this.vExhaust * this.massFlowRate * 0.98 > absDragForce) {
		this.massFlowRate *= 0.99;
		+ console.log('Limiting Mfr :', this.massFlowRate);
		}

		- if (this.exhaustVelocity * this.massFlowRate < Math.abs(dragForce[1])) {
		- this.massFlowRate = Math.abs(dragForce[1]) / this.exhaustVelocity * 1.1;
		+ if (this.vExhaust * this.massFlowRate < absDragForce) {
		+ this.massFlowRate = absDragForce / this.vExhaust * 1.1;
		}
		} else {

		this.blackBox.verticalVelocity(time, this.vVelocity);
		this.blackBox.verticalAcceleration(time, this.vAcceleration);
		- this.blackBox.thrustAcceleration(time, this.totalThrust / this.mass);
		- this.blackBox.mass(time, this.mass);
		+ this.blackBox.thrustAcceleration(time, this.totalThrust / this.wet);
		+ this.blackBox.mass(time, this.wet);
		+
		+ const magnitude = Math.sqrt(this.hAcceleration 2 + this.vAcceleration 2) / 9.8;
		+ const percentSpeed = this.hVelocity / this.targetVelocity * 100.0;
		+
		+ console.info('time, alt, hVel :',
		+ time,
		+ this.altitude,
		+ this.hVelocity, percentSpeed, this.wet, magnitude)
		}
		}

calculators/rocket/payload/calculator.js


		// Inputs
		- const SURFACE_TEMPERATURE = 25; // Celsius
		- const RELATIVE_HUMIDITY = 86.34; // per cent
		- const INITIAL_VELOCITY = 8; // Mach
		- const INITIAL_LATITUDE = 1.469167; // meters
		- const INITIAL_ALTITUDE = 6212; // meters
		+ const INITIAL_VELOCITY = 8.0; // Mach
		+ const INITIAL_LATITUDE = 1.469167; // degrees
		+ const INITIAL_ALTITUDE = 6212.0; // meters
		const TARGET_ALTITUDE = 400.0 * 1000.0; // meters

		const ROCKET_DIAMETER = 0.6;
		- const WET_MASS = 250;
		- const PAYLOAD_MASS = 25;
		- const DRAG_COEFFICIENT = 0.219;
		- const SPECIFIC_IMPULSE = 1400.0; // seconds
		+ const WET_MASS = 250.0;
		+ const PAYLOAD_MASS = 25.0;
		+ const DRAG_COEFFICIENT = 0.29;
		+ const SPECIFIC_IMPULSE = 1700.0; // seconds

		- const earth = new Earth(SURFACE_TEMPERATURE, RELATIVE_HUMIDITY);
		+ const earth = new Earth(INITIAL_LATITUDE);

		const SPEED_OF_SOUND = earth.soundSpeed(INITIAL_ALTITUDE);
		+ console.info( `sound speed = ${SPEED_OF_SOUND}` );
		+
		const INITIAL_HORIZONTAL_VELOCITY = 0;
		const INITIAL_VERTICAL_VELOCITY = INITIAL_VELOCITY * SPEED_OF_SOUND;


		while (rocket.flying()) {
		- time += 0.01;
		- time = parseFloat(time.toFixed(2));
		+ time = parseFloat((time + 0.01).toFixed(2));
		rocket.fly(time);
		}

Separates atmosphere, fixes mach speed conversion

Author	djarvis <email>
Date	2024-12-02 01:45:14 GMT-0800
Commit	1ed186144c7d72934fb472519145403798ca1082
Parent	881e556
Delta	285 lines added, 134 lines removed, 151-line increase