# ballistics2D.frink

```/** 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      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) // // https://frinklang.org/fsp/solve2.fsp?equations=area+%3D+pi+radius%5E2%0D%0Avolume+%3D+4%2F3+pi+radius%5E3%0D%0Adiameter+%3D+2+radius%0D%0Amass+%3D+volume+density%0D%0A%0D%0A&solveFor=&f=on&ev=on&sel_area=L&val_area=33+cm%5E2&sel_density=L&val_density=7+water&sel_diameter=S&val_diameter=&sel_mass=L&val_mass=1+kg&sel_pi=L&val_pi=3.141592653589793238&sel_radius=S&val_radius=&sel_volume=S&val_volume=&resultAs=cm&solved=on&showorig=on 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"] ```