Download or view X-15.frink in plain text format
/** This is a 2-dimensional ballistics / aerospace simulation that accurately
models air resistance (using the U.S. Standard Atmosphere as implemented
in StandardAtmosphere.frink ) 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
// h0 is height above ground (initial height)
// Highest X-15 flight, flight 149, 22 August 1963
// https://en.wikipedia.org/wiki/List_of_X-15_flights
h0 = 354200. ft
h = h0
// m1 is the mass of the projectile
// https://en.wikipedia.org/wiki/North_American_X-15
//m1 = 33500 lb // Gross weight
m1 = 14600. lb // Empty weight
// Cd is drag coefficient of projectile
//
// THERE IS LOTS OF DATA HERE FOR THE X-15
// https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=e65b16435676e88edffaaba028ee0d03abc2bd93
Cd = 0.0645
// 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)
// THERE IS LOTS OF DATA HERE FOR THE X-15
// https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=e65b16435676e88edffaaba028ee0d03abc2bd93
area = 13.60 ft^2
//println["Area is " + (area->"cm^2")]
//area = 0.7 m^2
// Initial velocity in the x direction (horizontal)
// (Orbital velocity can be found as v = sqrt[G planetmass / r] )
vx = 3794 mph
// Initial velocity in the y direction (vertical, positive is up)
vy = 0 km/s
planetmass = 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 sqrt[vx^2 + vy^2]^2
timestep = .1 s
t = 0 s
// Energy lost to drag
Edrag = 0 J
g = new graphics
p = new polyline
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 planetmass / 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
p.addPoint[x/m,-y/m]
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[vx,"mph",3] + "\t" + format[vy,"mph",3] + "\t" + format[v,"mph",3] + "\t" + format[Edrag, "J", 3] + "\t" + format[dragpow,"W", 3]]
} while h >= 0 ft
initialEnergy = initialGeopotentialEnergy + initialKineticEnergy
println["Initial potential energy = $initialGeopotentialEnergy"]
println["Initial kinetic energy = $initialKineticEnergy"]
println["Initial total energy = " + (initialGeopotentialEnergy + initialKineticEnergy)]
println["Final kinetic energy = " + (1/2 m1 v2)]
println["Energy lost to drag = $Edrag"]
println["Fraction of energy lost to drag = " + format[Edrag / initialEnergy, "percent", 4]]
[left, top, right, bottom] = getBoundingBox[p]
g.saveClip[]
g.clipRectSides[left,top,right,bottom]
g.color[0,0,1]
g.fillEllipseCenter[0,0,2 earthradius/m, 2 earthradius/m]
g.color[0,0,0]
g.add[p]
g.restoreClip[]
g.font["SansSerif", 10 km/m]
g.caption[formatFix[x,"miles", 0]]
g.caption[formatFix[h0,"miles", 0], "left", 90 deg]
g.show[]
g.write["X-15.png",400,undef]
// Frink solver:
// https://frinklang.org/fsp/solve2.fsp?equations=E+%3D+mass+gravity+height%0D%0AE+%3D+1%2F2+mass+velocity%5E2%0D%0Adistance+%3D+height+glideslope%0D%0A&solveFor=&f=m&ev=on&sel_distance=S&val_distance=&sel_E=S&val_E=&sel_glideslope=L&val_glideslope=+4&sel_gravity=L&val_gravity=196133%2F20000+m+s%5E-2&sel_height=S&val_height=+&sel_mass=S&val_mass=+&sel_velocity=L&val_velocity=4554.+mph&resultAs=miles&solved=on&showorig=on
glideslope = 4/1
maxdistance = glideslope v^2 / (2 gravity)
println["Maximum distance at glideslope of $glideslope:1 is " + format[maxdistance, "miles", 0]]
Download or view X-15.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, eliasen@mindspring.com