/** 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 } }