//This code was compiled by "g++ bolidesimul.cpp" and written in Dev-C++
//(available under the GNU General Public License at
//
//output is stored to "data.txt," which is nicely viewable through an editor like
//emacs (available under the GNU General Public License at
//
//
//Notes:
//
//The code is currently set up to compute the entry of a stony bolide of mass
//5.6*10^8 kg, initial speed 30 km/s (so Ekin0 = 60MTon) initial theta 30 degrees
//vardot denotes the time derivative of var. Other derivatives are written as
//dvar_dsomeothervar
//
//Many constants and variables are given a #define var_name vec[index] statement,
//their name followed by an underscore. This is because it was convient to pass
//all the variables or all the constants as a single vector to functions, but it
//would be inconvient to write formulas in terms of the variables' indeces inside
//the vector.
//
//All units are SI, though some conversions are made when the results are written
//to a file
#include
#include
float timeDeriv(int index, float sVec[], float cVec[]); //takes the time derivative of the variable in sVec[index], sVec is the state vector, cVec is a vector of constants
float min(float x, float y); //returns minimum of x,y
float rho_a(float h); //returns atmospheric density at altitude h
float g(float h); //returns gravitational acceleration at altitude h
#define v_ sVec[0]
#define theta_ sVec[1]
#define h_ sVec[2]
#define m_ sVec[3]
#define R_ sVec[4]
#define Rdot_ sVec[5]
#define Ekin_ sVec[6]
#define Ekindot_ sVec[7]
#define A_ sVec[8]
#define dEkin_dh_ sVec[9]
#define fragmented_ sVec[10]
#define pi_ cVec[0]
#define sigma_ cVec[1]
#define T_ cVec[2]
#define CH_ cVec[3]
#define CL_ cVec[4]
#define CD_ cVec[5]
#define Q_ cVec[6]
#define SF_ cVec[7]
#define rho_m_ cVec[8]
#define S_ cVec[9]
main()
{
const float dt = 0.001; //s, simulation time step
const int maxTimeSteps = 30000;
int stepsTaken = 0; //counts number of steps taken in simulation in case we stop before maxTimeSteps, as when h < 0 or m < 0
//put constants into handy constVector cVec
float cVec[10];
//constants
pi_ = 3.14159;
sigma_ = 5.6704*pow(10.0,-8.0); //W/(m^2 K^4) Stefan-Boltzmann constant
T_ = 25000; //K, shock wave gas temperature
CH_ = 0.1; //heat transfer coefficient
CL_ = 0.001; //lift coefficient
//constants related to bolide
CD_ = 0.5; //drag coefficient for rough sphere
Q_ = 8*1000000; //J/kg heat of ablation for stony asteroid
SF_ = 1.21; //shape factor for sphere, SF is solution to pi*R^2=SF*(M/rho_m)^(2/3) for sphere
rho_m_ = 3500; //kg/m^3 stony asteroid density
S_ = pow(10.0,7.0); //Pa yield strength of stony asteroid
//put variables into handy stateVector sVec
float sVec[11];
//initial parameters
v_ = 30000; //m/s, initial speed of bolide
theta_ = 30*pi_/180; //in degrees, measured from horizon
h_ = 100000; //m initial altitude
m_ = 5.6*pow(10.0,8.0); //kg, initial mass of bolide
A_ = SF_*pow(m_/rho_m_,2.0/3.0); //cross-sectional area, pow(x,y) = x^y
R_ = sqrt(A_/pi_);
Rdot_ = 0; //rate of change of radius, useful after t > t_crit
Ekin_ = 0.5*m_*v_*v_; //initial kinetic energy in J
fragmented_ = 0; //initially not fragmented
Ekindot_ = timeDeriv(6,sVec,cVec); //time derivative of kinetic energy in J/s
dEkin_dh_ = Ekindot_/(-v_*sin(theta_)); //J/m dEkin/dh = (dE/dt)*(dt/dh) = Ekindot/(-v*sin(theta)) derivative of energy with respect to altitude
float dataArray[maxTimeSteps][11]; //stores an instance of the state vector at each time step
//run simulation
for (int step = 0; step < maxTimeSteps; step++)
{
//store current state vector to data array
for(int j = 0; j < 11; j++) //11 variables
dataArray[step][j] = sVec[j];
//update the state vector
//Runge-Kutta 4 to integrate equations of motion for v and theta
for(int j = 0; j <= 1; j++) //these are of form y'(t) ~ f(y(t)) with no explicit t dependence
{
float k1 = timeDeriv(j,sVec,cVec); //k1 = f(y_n)
sVec[j] += (dt/2.0)*k1;
float k2 = timeDeriv(j,sVec,cVec); //k2 = f(y_n + (dt/2)*k1)
sVec[j] += -(dt/2.0)*k1 + (dt/2.0)*k2;
float k3 = timeDeriv(j,sVec,cVec); //k3 = f(y_n + (dt/2)*k2)
sVec[j] += -(dt/2.0)*k2 + dt*k3;
float k4 = timeDeriv(j,sVec,cVec); //k4 = f(y_n + dt*k3)
sVec[j] += -dt*k3; //restores sVec[k] to original value
sVec[j] += (dt/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4); //y_(n+1) = y_n + (dt/6)*(k1+2*k2+2*k3+k4)
}
for (int j = 2; j <= 3; j++) //h'(t) and m'(t) do not depend on derivatives of h(t) and m(t) or on t,
sVec[j] += dt*timeDeriv(j,sVec,cVec); //so k1=k2=k3=k4 and the RK4 method reduces to f -> f + dt*f'
//for post-fragmentation dynamics of R
//r'' = const/r, so we make this system first order by defining y1 = r, y2 = r'
//z = (y1 ; y2), z' = (y1' ; y2') = (y2 ; const/y1) and integrate this using Runge-Kutta 4
//in our case z = {r rdot}, z' = {rdot rdotdot}
if (fragmented_ > 0.5)
{
float k1[2] = {timeDeriv(4,sVec,cVec), timeDeriv(5,sVec,cVec)}; //{rdot, rdotdot}
R_ += (dt/2.0)*k1[0]; //(y1)
Rdot_ += (dt/2.0)*k1[1]; //(y2)
float k2[2] = {timeDeriv(4,sVec,cVec), timeDeriv(5,sVec,cVec)};
R_ += -(dt/2.0)*k1[0] + (dt/2.0)*k2[0];
Rdot_ += -(dt/2.0)*k1[1] + (dt/2.0)*k2[1];
float k3[2] = {timeDeriv(4,sVec,cVec), timeDeriv(5,sVec,cVec)};
R_ += -(dt/2.0)*k2[0] + dt*k3[0];
Rdot_ += -(dt/2.0)*k2[1] + dt*k3[1];
float k4[2] = {timeDeriv(4,sVec,cVec), timeDeriv(5,sVec,cVec)};
R_ += -dt*k3[0];
Rdot_ += -dt*k3[0];
R_ += (dt/6.0)*(k1[0] + 2.0*k2[0] + 2.0*k3[0] + k4[0]);
Rdot_ += (dt/6.0)*(k1[1] + 2.0*k2[1] + 2.0*k3[1] + k4[1]);
}
else //prefragmentation
{
R_ = sqrt(A_/pi_);
Rdot_ = 0; //this does not affect pre-fragmentation dynamics
}
//update other variables
Ekin_ = 0.5*m_*v_*v_; //Ekin = 0.5*m*v^2
Ekindot_ = timeDeriv(6,sVec,cVec); //Ekindot
//update cross-section area
if (fragmented_ < 0.5) //not fragmented
A_ = SF_*pow(m_/rho_m_,2.0/3.0);
else
A_ = pi_*R_*R_;
//check for fragmentation
if (CD_*rho_a(h_)*v_*v_/4.0 > S_) //equation (7)
fragmented_ = 1.0;
else
{
if (fragmented_ > 0.5) //if has already fragmented, stays fragmented
fragmented_ = 1.0;
else
fragmented_ = 0.0;
}
//this is the variable we plot in figures 2-6
dEkin_dh_ = Ekindot_/(-v_*sin(theta_)); //dEkin/dh = (dE/dt)*(dh/dt)^(-1) = Ekindot/(-v*sin(theta))
//check to see if should halt if h < 0 (impacted surface) or m < 0 (ablation has depleted mass to 0)
if ((h_ < 0) || (m_ < 0))
break;
else
stepsTaken++;
}
//writes data to file
FILE * file;
file = fopen("data.txt","w");
fprintf(file,"t\t\tv\t\ttheta\t\th\t\tm\t\tR\t\tRdot\t\tEkin\t\tEkindot\tA\t\tdEkin_dh\tfragmented\n"); //header line
for (int j = 0; j < stepsTaken; j++)
{
fprintf(file,"%.3e\t",j*dt); //prints time
for (int k = 0; k < 10; k++)
{
if (k == 2) //h
dataArray[j][k] = dataArray[j][k]/1000.0; //convert h to km
if (k == 9) //dEkin_dh
dataArray[j][k] = dataArray[j][k]/(4.184*pow(10.0,12.0)); //convert dEkin_dh to MTon/km
fprintf(file,"%.3e\t",dataArray[j][k]); //prints data
}
fprintf(file,"%i\n",int(dataArray[j][10])); //prints 0 if not fragmented, 1 otherwise
}
fclose(file);
}
float timeDeriv(int index, float sVec[], float cVec[]) //takes the time derivative of the variable in sVec[index], sVec is the state vector, cVec is a vector of constants
{ //cVec is vector of constants
switch(index)
{
case 0: //vdot, equation (2)
return (-(CD_/m_)*rho_a(h_)*A_*v_*v_ + g(h_)*sin(theta_));
break;
case 1: //thetadot, equation (5)
return (g(h_)*cos(theta_)/v_ - CL_*rho_a(h_)*A_*v_/(2*m_));
break;
case 2: //hdot, equation (6)
return (-v_*sin(theta_));
break;
case 3: //mdot, equation (3)
return (-(A_/Q_)*min((CH_*rho_a(h_)*pow(v_,3.0)),sigma_*pow(T_,4.0)));
break;
case 4: //if fragmented Rdot = Rdot (from Rdot_)
if (fragmented_ > 0.5) //if fragmented
return (Rdot_);
else
return 0; //this does not affect pre-fragmentation dynamics
break;
case 5: //Rdotdot, equation (9)
if (fragmented_ > 0.5) // if fragmented
return (CD_*rho_a(h_)*v_*v_/(16.0*R_*rho_m_/(3.0*pi_)));
else
return 0; //this does not affect pre-fragmentation dynamics
break;
case 6: //Ekindot, equation (1)
return (m_*v_*timeDeriv(0,sVec,cVec) + 0.5*v_*v_*timeDeriv(3,sVec,cVec));
break;
default:
printf("Warning: timeDeriv called with uknown index value\n");
return 0;
break;
}
}
float min(float x, float y)
{
if (x