/** 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.) It can be used to model anything from orbital decay to hypervelocity meteor impacts to orbital decay due to atmosphere. Using the defaults below, an orbit of a small iron sphere starts approximately circularly at an altitude of 180 km above earth's surface and decays over the course of several days until reeentry. The variant meteorCrater.frink models a hypervelocity impactor for Meteor Crater, Arizona and shows the huge amount of energy radiated by pushing through atmosphere. */ 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) // // https://frinklang.org/fsp/solve2.fsp?equations=area+%3D+pi+radius%5E2%0D%0Avolume+%3D+4%2F3+pi+radius%5E3%0D%0Adiameter+%3D+2+radius%0D%0Amass+%3D+volume+density%0D%0A%0D%0A&solveFor=&f=on&ev=on&sel_area=L&val_area=33+cm%5E2&sel_density=L&val_density=7+water&sel_diameter=S&val_diameter=&sel_mass=L&val_mass=1+kg&sel_pi=L&val_pi=3.141592653589793238&sel_radius=S&val_radius=&sel_volume=S&val_volume=&resultAs=cm&solved=on&showorig=on 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"]