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 19985 days, 23 hours, 45 minutes ago.