# meteorCrater.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.)     This program models the entry of the Meteor Crater, Arizona impactor     through layers of atmosphere and demonstrates the huge energy radiation     due to pushing through layers of atmosphere (463 kilotons TNT).     This program can be altered to model many types of impactors.  Some may     slow down to terminal velocity. */ use StandardAtmosphere.frink /*for v = mach 0 to mach 5 step mach 1/10    println[(v->mach) + "\t" + Cd[v]] exit[]*/ // h is height above ground (initial height) h = 900 km radius =  50 m / 2 volume = 4/3 pi radius^3 density = 7.874 g/cm^3 // m1 is the mass of the projectile m1 = volume density // Cd is drag coefficient of projectile Cd = 1.4 // 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.874 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 = 9.19 km/s // Initial velocity in the y direction (vertical, positive is up) vy = -9.19 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 = -700 km 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 = .01 s t = 0 s // Energy lost to drag Edrag = 0 J g = new graphics line = new polyline height = new polyline 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]    // Calculate drag coefficient as a function of velocity    Cd = Cd[v]    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    // Power due to drag    Pdrag = fdrag v        // Energy lost to drag    // E = f * d = f * v * t    Edrag = Edrag + Pdrag timestep        x = x + dx    y = y + dy    line.addPoint[x/m, -y/m]    height.addPoint[t, -h/m]    geopotentialHeight = (r * h) / (r + h)    geopotentialEnergy = m1 gravity geopotentialHeight    kineticEnergy = 1/2 m1 (vx^2 + vy^2)    totalEnergy = geopotentialEnergy + kineticEnergy //if t mod (1 s) == 0 s       println[formatFixed[t,"s",3] + "\t" + format[x,"m",8] + "\t" + format[y,"m",9] + "\t" + format[h,"m",9] + "\t" + format[adragx,"gee",4] + "\t" + format[adragy,"gee",4] + "\t" + format[adrag,"gee",4] + "\t" + format[v,"m/s",7] + "\t" + formatEng[Pdrag, "W", 3]] } while h >= 5600 ft initialEnergy = initialGeopotentialEnergy + initialKineticEnergy println["Initial potential energy = \$initialGeopotentialEnergy"] println["Initial kinetic energy   = \$initialKineticEnergy"] println["Final kinetic energy     = " + (1/2 m1 v2)] println["Initial energy           = " + initialEnergy] println["Energy lost to drag      = \$Edrag\t(" + formatEng[Edrag, "tons TNT", 3] + ")"] println["Fraction of energy lost to drag = " + format[Edrag / initialEnergy, "percent", 4]] g.add[line] g.show[] g2=new graphics g2.add[height] g2.show[] // Drag coefficient equation for a sphere in the subsonic-to-hypersonic environment // See ESTIMATING THE DRAG COEFFICIENTS OF METEORITES FOR ALL MACH NUMBER REGIMES. // R. T. Carter, P. S. Jandir, and M. E. Kress, // 40th Lunar and Planetary Science Conference (2009) // They don't actually display the "quadratic" behavior which occurs below // about Mach 0.8, so this is a curve-fit. // For hypersonic velocities, this is about 0.92, while other authors state that // a constant 0.7 is okay. Cd[v] := {    M0 = v/mach        if (M0 >= 0.8)    {       // Mach 0.8 and above, use eq. from paper.       M1 = M0 + 0.35       return 2.1 e^(-1.2 M1) - 8.9 e^(-2.2 M1) + 0.92    } else    {       // Below Mach 0.8       // Use "quadratic" range picked off chart and fitted       return 0.424 + 0.472 M0^2    } } ```