弹箭三自由度质点弹道计算程序

作者在以前用c++写的三自由度弹道程序的基础上模仿matlab内置的ode45函数、微分方程组写法，对c++版三自由度的弹道微分方程组和龙格库塔函数的写法进行了优化。

一、环境参数及结构特征参数

示例：在这里选用的是155榴弹

const double PI = 3.1415926;
const double G0 = 9.7803253359;        // 地球表面重力加速度(m/s2)
const double rho0 = 1.205;             // 地面大气密度（kg/m3）
const double T0 = 288.15;              // 地面温度（K）
const double tatal0 = 289.1;           // 地面虚温（K）
const double P0 = 0.1e6;               // 地面气压（MPa）
const double R1 = 287;                 // 空气常数 J/KG/K
const double k = 1.4;		
	       
double i = 0.85;                       // 榴弹弹形系数
double G = 6.328e-3;
double Mass = 45.54;                    // 弹丸重量（kg）
double D = 0.15494;                      // 弹丸口径（m）
double C = i*D*D * 1000 / Mass;        // 弹道系数

二、建立存放参数的 vector 数组及定义其它函数

2.1 存放参数的 vector 数组

std::vector<double> t_s;       // 时间
std::vector<double> vx_s;      // 速度在x轴上的分量
std::vector<double> vy_s;	   // 速度在y轴上的分量
std::vector<double> vz_s;	   // 速度在z轴上的分量
std::vector<double> y_s;       // 射高
std::vector<double> x_s;       // 射程
std::vector<double> z_s;       // 横偏
std::vector<double> v_s;       // 速度
std::vector<double> theat_s;   // 弹道倾角

2.2 声明函数

// 弹道微分方程组函数
void dery(int n, double dy[], double Y[]);
// 龙格库塔函数
void RK4(double rk_step);
// 虚温随高度变化函数
double tatal(double y);
// 空气密度函数
double rho(double P, double y);
// 声速
double Cs(double y);
// 阻力函数
double Gvcs(double y, double v, double vr);
// 相对速度
double Vr(double vx, double vy, double vz, double wx, double wz);
// 绝对速度
double Velocity(double vx, double vy, double vz);
// 仿真计算程序函数
void calculateReentry();
// 输出结果文件
void outData();

2.3 定义弹体飞行参数

double wx = 0.0;               // 纵风（m/s)
double wz = 0.0;               // 横风（m/s)
double v = 897;                // 初速（m/s）
double theat = 50 * PI / 180;  // 初始弹道倾角（deg）
double x = 0;                  // 水平距离（m）
double y = 0;                  // 垂直距离（m）
double z = 0;                  // 横偏距离（m）
double P = P0;                 // 压强（MPa）
double vx = v * cos(theat);    // 速度在x轴上的分量
double vy = v * sin(theat);	   // 速度在y轴上的分量
double vz = 0;      		   // 速度在z轴上的分量
double dt = 0.01;              // 步长(s)
double ct = 0;                 // 初始时间(s)
const int NumbVars = 9;        // 参数变量个数
// 将参数初始值放入Y数组
double Y[NumbVars] = { 0,vx,vy,vz,x,y,z,P0,theat };

三、针对微分方程组函数和RK4函数展开介绍

3.1 微分方程组函数

这里的微分方程组写法模仿了matlab，大大缩小了代码量。

/*弹丸质心运动微分方程组*/
void dery(int n, double dy[], double Y[])
{
	v = Velocity(Y[1], Y[2], Y[3]);
	vr = Vr(Y[1], Y[2], Y[3], wx, wz);
	double a = rho(Y[7], Y[5]) / rho0;
	double b = Gvcs(Y[5], v, vr);
    
    // 微分方程组
	dy[0] = 1.0;
	dy[1] = -C * a * b * (Y[1] - wx);
	dy[2] = -C * a * b * Y[2] - G0;
	dy[3] = -C * a * b * (Y[3] - wz);
	dy[4] = Y[1];
	dy[5] = Y[2];
	dy[6] = Y[3];
	dy[7] = -rho(Y[7], Y[5])*G0*Y[2];
	dy[8] = -G0 * cos(Y[8]) / v;
}

3.2 RK4函数

这里的龙格库塔函数模仿了matlab中的ode45函数的写法，同样极大地缩小了代码量。

void RK4(double rk_step)
{
    // rk_step为dt，时间步长
	double a[4], old_y[NumbVars], Y1[NumbVars], *dy;
	int i, j;
	a[0] = rk_step / 2; 
	a[1] = rk_step / 2; 
	a[2] = rk_step; 
	a[3] = rk_step;
	dy = (double*)calloc(NumbVars, sizeof(double));
	dery(NumbVars, dy, Y);
	for (i = 0; i < NumbVars; i++)  old_y[i] = Y[i];
	for (j = 0; j < 3; j++)
	{
		for (i = 0; i < NumbVars; i++)
		{
			Y1[i] = old_y[i] + a[j] * dy[i];
			Y[i] = Y[i] + a[j + 1] * dy[i] / 3;
		}
		dery(NumbVars, dy, Y1);
	}
	for (i = 0; i < NumbVars; i++)  Y[i] = Y[i] + a[0] * dy[i] / 3;
	free(dy);
}