Download or view ballistics2D.frink in plain text format
/** This is a 2-dimensional ballistics / aerospace simulation that accurately
models air resistance (using the U.S. Standard Atmosphere) and the
curvature of the earth (in 2 dimensions.)
It can be used to model anything from orbital decay to hypervelocity meteor
impacts to orbital decay due to atmosphere.
Using the defaults below, an orbit of a small iron sphere starts
approximately circularly at an altitude of 180 km above earth's surface
and decays over the course of several days until reeentry.
The variant meteorCrater.frink
models a hypervelocity impactor for Meteor Crater, Arizona and shows
the huge amount of energy radiated by pushing through atmosphere.
*/
use StandardAtmosphere.frink
// h is height above ground (initial height)
h = 180 km
// m1 is the mass of the projectile
m1 = 1 kg
// Cd is drag coefficient of projectile
Cd = 1
// area is area of projectile in direction of travel.
// area of a sphere can be calculated from mass and density as:
// area = (3/4)^(2/3) m1^(2/3) pi^(1/3) density^(-2/3)
//
// https://frinklang.org/fsp/solve2.fsp?equations=area+%3D+pi+radius%5E2%0D%0Avolume+%3D+4%2F3+pi+radius%5E3%0D%0Adiameter+%3D+2+radius%0D%0Amass+%3D+volume+density%0D%0A%0D%0A&solveFor=&f=on&ev=on&sel_area=L&val_area=33+cm%5E2&sel_density=L&val_density=7+water&sel_diameter=S&val_diameter=&sel_mass=L&val_mass=1+kg&sel_pi=L&val_pi=3.141592653589793238&sel_radius=S&val_radius=&sel_volume=S&val_volume=&resultAs=cm&solved=on&showorig=on
density = 7 g/cm^3
area = (3/4)^(2/3) m1^(2/3) pi^(1/3) density^(-2/3)
//println["Area is " + (area->"cm^2")]
//area = 0.7 m^2
// Initial velocity in the x direction (horizontal)
vx = 7.82 km/s
// Initial velocity in the y direction (vertical, positive is up)
vy = 0 km/s
// mp is mass of planet
mp = earthmass
r = earthradius
// x and y are a cartesian coordinate system with the center of the planet
// at x=0 m, y=0 m. The projectile typically begins its journey at x=0 and
// at a given height-above-ground.
x = 0 m
y = r + h
initialGeopotentialHeight = (r * h) / (r + h)
//println["Geopotential height = " + (geopotentialHeight -> "ft")]
initialGeopotentialEnergy = m1 gravity initialGeopotentialHeight
initialKineticEnergy = 1/2 m1 (vx^2 + vy^2)
timestep = .1 s
t = 0 s
// Energy lost to drag
Edrag = 0 J
do
{
// l is distance from center of earth
l2 = x^2 + y^2
l = sqrt[l2]
h = l - r
// Angle with respect to center of the earth
alpha = arctan[x,y]
// Force due to gravity
fg = - G m1 mp / l2
// Acceleration due to gravity
ag = fg / m1
agx = ag sin[alpha]
agy = ag cos[alpha]
// Atmospheric drag
v2 = vx^2 + vy^2
v = sqrt[v2]
// Angle of travel (0 is in x direction, 90 degrees in y direction)
theta = arctan[vy, vx]
[temp, pressure] = StandardAtmosphere.getTemperatureAndPressure[h]
density = StandardAtmosphere.getDensity[h, temp, pressure]
fdrag = 1/2 density v2 Cd area
adrag = -fdrag / m1
adragx = adrag cos[theta]
adragy = adrag sin[theta]
t = t + timestep
// Total acceleration
axtotal = agx + adragx
aytotal = agy + adragy
atotal = sqrt[axtotal^2 + aytotal^2]
dvx = axtotal timestep
dvy = aytotal timestep
vx = vx + dvx
vy = vy + dvy
dx = vx timestep
dy = vy timestep
// Energy lost to drag
// E = f * d = f * v * t
dragpow = fdrag v
Edrag = Edrag + fdrag v timestep
x = x + dx
y = y + dy
geopotentialHeight = (r * h) / (r + h)
geopotentialEnergy = m1 gravity geopotentialHeight
kineticEnergy = 1/2 m1 (vx^2 + vy^2)
totalEnergy = geopotentialEnergy + kineticEnergy
// println["agx=$agx\tagy=$agy"]
// println["adragx=$adragx\tadragy=$adragy"]
// println[format[t,"s",3] + "\t" + format[x,"m",3] + "\t" + format[y,"m", 3] + "\t" + format[h,"m",3] + "\t" + format[agx,"gee",3] + "\t" + format[agy,"gee",3]]
if t mod (1 s) == 0 s
println[format[t,"s",3] + /*"\t" + format[x,"m",3] + "\t" + format[y,"m", 3] + */ "\t" + format[h,"m",3] + /*"\t" + format[adragx,"gee",8] + "\t" + format[adragy,"gee",8] +*/ "\t" + format[adrag,"gee",8] + "\t" + format[v,"mph",3] + "\t" + format[Edrag, "J", 3] + "\t" + format[dragpow,"W", 3]]
} while h >= 0 ft
println["Initial potential energy = $initialGeopotentialEnergy"]
println["Initial kinetic energy = $initialKineticEnergy"]
println["Final kinetic energy = " + (1/2 m1 v2)]
println["Initial energy = " + (initialGeopotentialEnergy + initialKineticEnergy)]
println["Energy lost to drag = $Edrag"]
Download or view ballistics2D.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen was born 19965 days, 22 hours, 40 minutes ago.