Download or view MarsMeteor.frink in plain text format
/** This is a 2-dimensional ballistics / aerospace simulation that attempts
to model a Mars meteorite re-entry.
This program can be altered to model many types of impactors. Some may
slow down to terminal velocity (which means they probably actually
burned up.)
*/
use MarsAtmosphere.frink
// h is height above ground (initial height)
h = 14 km
diameter = 1 cm
radius = diameter/2
volume = 4/3 pi radius^3
density = 7.874 g/cm^3
// m1 is the mass of the projectile
m1 = volume density
println["Mass is $m1"]
// 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)
//area = (3/4)^(2/3) m1^(2/3) pi^(1/3) density^(-2/3)
area = pi radius^2
println["Area is " + (area->"cm^2")]
//area = 0.7 m^2
// Initial velocity in the x direction (horizontal)
vx = 0 km/s
// Initial velocity in the y direction (vertical, positive is up)
vy = -13 km/s
// mp is mass of planet
mp = marsmass
r = marsradius
gravity = G mp / r^2
// 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 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 = .005 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 Mars
l2 = x^2 + y^2
l = sqrt[l2]
h = l - r
// Angle with respect to center of Mars
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, adensity] = tpd = MarsAtmosphere.getTPD[h]
// Calculate drag coefficient as a function of velocity
Cd = Cd[v]
fdrag = 1/2 adensity v2 Cd area
// println["tpd=$tpd, Fdrag = $fdrag"]
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 >= 0 m
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 MarsMeteor.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 19979 days, 8 hours, 50 minutes ago.