Download or view meteorCrater.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.)
This program models the entry of the Meteor Crater, Arizona impactor
through layers of atmosphere and demonstrates the huge energy radiation
due to pushing through layers of atmosphere (463 kilotons TNT).
This program can be altered to model many types of impactors. Some may
slow down to terminal velocity.
*/
use StandardAtmosphere.frink
/*for v = mach 0 to mach 5 step mach 1/10
println[(v->mach) + "\t" + Cd[v]]
exit[]*/
// h is height above ground (initial height)
h = 900 km
radius = 50 m / 2
volume = 4/3 pi radius^3
density = 7.874 g/cm^3
// m1 is the mass of the projectile
m1 = volume density
// Cd is drag coefficient of projectile
Cd = 1.4
// 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)
density = 7.874 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 = 9.19 km/s
// Initial velocity in the y direction (vertical, positive is up)
vy = -9.19 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 = -700 km
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 = .01 s
t = 0 s
// Energy lost to drag
Edrag = 0 J
g = new graphics
line = new polyline
height = 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 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]
// Calculate drag coefficient as a function of velocity
Cd = Cd[v]
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
// Power due to drag
Pdrag = fdrag v
// Energy lost to drag
// E = f * d = f * v * t
Edrag = Edrag + Pdrag timestep
x = x + dx
y = y + dy
line.addPoint[x/m, -y/m]
height.addPoint[t, -h/m]
geopotentialHeight = (r * h) / (r + h)
geopotentialEnergy = m1 gravity geopotentialHeight
kineticEnergy = 1/2 m1 (vx^2 + vy^2)
totalEnergy = geopotentialEnergy + kineticEnergy
//if t mod (1 s) == 0 s
println[formatFixed[t,"s",3] + "\t" + format[x,"m",8] + "\t" + format[y,"m",9] + "\t" + format[h,"m",9] + "\t" + format[adragx,"gee",4] + "\t" + format[adragy,"gee",4] + "\t" + format[adrag,"gee",4] + "\t" + format[v,"m/s",7] + "\t" + formatEng[Pdrag, "W", 3]]
} while h >= 5600 ft
initialEnergy = initialGeopotentialEnergy + initialKineticEnergy
println["Initial potential energy = $initialGeopotentialEnergy"]
println["Initial kinetic energy = $initialKineticEnergy"]
println["Final kinetic energy = " + (1/2 m1 v2)]
println["Initial energy = " + initialEnergy]
println["Energy lost to drag = $Edrag\t(" + formatEng[Edrag, "tons TNT", 3] + ")"]
println["Fraction of energy lost to drag = " + format[Edrag / initialEnergy, "percent", 4]]
g.add[line]
g.show[]
g2=new graphics
g2.add[height]
g2.show[]
// Drag coefficient equation for a sphere in the subsonic-to-hypersonic environment
// See ESTIMATING THE DRAG COEFFICIENTS OF METEORITES FOR ALL MACH NUMBER REGIMES.
// R. T. Carter, P. S. Jandir, and M. E. Kress,
// 40th Lunar and Planetary Science Conference (2009)
// They don't actually display the "quadratic" behavior which occurs below
// about Mach 0.8, so this is a curve-fit.
// For hypersonic velocities, this is about 0.92, while other authors state that
// a constant 0.7 is okay.
Cd[v] :=
{
M0 = v/mach
if (M0 >= 0.8)
{
// Mach 0.8 and above, use eq. from paper.
M1 = M0 + 0.35
return 2.1 e^(-1.2 M1) - 8.9 e^(-2.2 M1) + 0.92
} else
{
// Below Mach 0.8
// Use "quadratic" range picked off chart and fitted
return 0.424 + 0.472 M0^2
}
}
Download or view meteorCrater.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 19972 days, 23 hours, 59 minutes ago.