# 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.) */ 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"] ```

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 18228 days, 13 hours, 26 minutes ago.