三自由度质点弹道方程（vs2017 c++ ）

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

作者在以前用c++写的三自由度弹道程序的基础上模仿matlab，针对弹道微分方程组和RK4函数的写法进行了优化

示例：

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;        // 弹道系数

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

// 弹道微分方程组函数
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();

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 };

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;
}

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);
}

开始运行三自由度弹道程序：

仿真结果文件位置：

结果曲线：

（1）射程-射高

（2）时间-速度

（3）时间-弹道倾角