model.cpp

用来计算航天器再入质点弹道

#include"插值.h"
#include "AirModel.h"
//**********************************************************************************

//*******************************************************************
double Angle_function(double v)
{
	if(v>800&&v<=1700)
		 return 14.0;
	else if(v>1700&&v<=2500)
		 return 40.0-(2500-v)*26.0/1700;
	else if(v>2500)
		 return 40.0;
    else return 14.0;
}
//*******************************************************************
void Model(double *y,double *f)
{
	double  delta           =		0.0;
    double	D				=		0.0;
	double	L				=		0.0;
 	double	density			=		0.0;
	double	Cd				=		0.0;
	double  Cl				=		0.0;
	double  speedofsound	=		0.0;
	double  pressure		=		0.0;
    double	Sref			=		140;
	double	m				=		28000.0;
	double	g0				=		9.81;
	double	R0				=		6378135.0;
	double	v				=		y[3]*sqrt(g0*R0);
    double	angle			=		Angle_function(v);
    double	altitude		=		(y[0]-1.0)*R0;
	
	InitCalcAtmosCOESA();
	CalcAtmosCOESA(altitude,density,speedofsound,pressure);
    double	ma				=		v/speedofsound;
	
	Chazhi	function;
    Cd		=	function.insert("Data_Cx_Vec.txt",2,ma,angle,0);
    Cl		=   function.insert("Data_Cy_Vec.txt",2,ma,angle,0);

	D		=	density*g0*R0*y[3]*y[3]*Sref*Cd/(2*m*g0); 
	L		=	density*g0*R0*y[3]*y[3]*Sref*Cl/(2*m*g0);
	
	f[0]	=	y[3]*sin(y[4]);
	f[1]	=	y[3]*cos(y[4])*sin(y[5])/(y[0]*cos(y[2]));
	f[2]	=	y[3]*cos(y[4])*cos(y[5])/y[0];
	f[3]	=	-D-sin(y[4])/(y[0]*y[0]);
	f[4]	=	(L*cos(delta)+(y[3]*y[3]-1/y[0])*cos(y[4])/y[0])/y[3];
	f[5]	=	(L*sin(delta)/cos(y[4])+y[3]*y[3]*cos(y[4])*sin(y[5])*sin(y[2])/(cos(y[2])*y[0]))/y[3];
   
	printf("Cd=%lf  Cl=%lf \n",Cd,Cl);
	
}