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

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 18350 days, 1 hours, 47 minutes ago.