ballistics2D.frink

Download or view 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.)

    It can be used to model anything from orbital decay to hypervelocity meteor
    impacts to orbital decay due to atmosphere.

    Using the defaults below, an orbit of a small iron sphere starts
    approximately circularly at an altitude of 180 km above earth's surface
    and decays over the course of several days until reeentry.

    The variant meteorCrater.frink 
    ( https://frinklang.org/fsp/colorize.fsp?f=meteorCrater.frink )
    models a hypervelocity impactor for Meteor Crater, Arizona and shows
    the huge amount of energy radiated by pushing through atmosphere.
*/


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


Download or view 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 19350 days, 23 hours, 36 minutes ago.