ballistics2D.frink

View or download 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.)
*/


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)

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"]


View or download 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 17927 days, 13 hours, 36 minutes ago.