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