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
   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" + (Edrag -> "J")]
} 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 17645 days, 8 hours, 32 minutes ago.