标准弹道.cpp

在无投弹误差条件下

	S=0.07022;                                                 //参考面积
	SONIC=20.046 * sqrt(288.34-(5.86e-3)*Y[4]);                //声速
	RHO=1.225*pow((1.0-(2.0323e-5)*Y[4]),4.83);                //密度
	g=9.806*(1.0-2*Y[4]/6371000);                              //重力加速度g
	q=RHO*Y[1]*Y[1]/2;                                         //0.5*P*V*V
	MA=Y[1]/SONIC;                                             //马赫数
    
	ABC_ALPHA=fabs(ALPHA);                                     //正负攻角差值相同
	interp();                                                  //插值函数
	a[0]=sin(Y[2]);
	a[1]=cos(Y[2]);

    XF=CX*q*S;
	YF=(CY0*ALPHA+CY1*ALPHA*ALPHA*ALPHA)*q*S;
//	if(ALPHA>=0) YF=CY*q*S;
//	else YF=-CY*q*S;

	dy[0]=1;
	dy[1]=(-XF-g*220*a[0])/220;
	dy[2]=(YF-g*220*a[1])/(220*Y[1]);
	dy[3]=Y[1]*a[1];
	dy[4]=Y[1]*a[0];
  
	dX=dy[3];                                                  //水平速度
    dY=dy[4];                                                  //垂直速度
	dTheta=dy[2];                                              //弹道倾角
//  guozai=fabs(YF/225/g);                                     //计算需用过载
	return;
}
/*************************************************************************************/
/*                         调用差值函数计算气动力系数                                */
/*                                interp()                                           */
/*                       无返回值，有一维和二维差值两种                              */
/*************************************************************************************/
void interp()
{
    CX=int32(ABC_ALPHA,MA,7,5,GJ,MHS,ACX);                    //二维差值
//	CY=int32(ABC_ALPHA,MA,10,5,GJ,MHS,ACY);
	CY0=int11(ACY0,5,MA);                                     //一维差值
	CY1=int11(ACY1,5,MA);                                     //一维差值
//	CX=int11(ACY,6,MA);
	return;
}
/*************************************************************************************/
/*                              二维差值函数                                         */
/*                                 int32()                                           */
/*        x为要差得行值，qq为要差得列值，N1为行数，N2为列数，A为行向量，BB为列向量   */          
/*                                 XS为二维表                                        */
/*************************************************************************************/
double int32(double x,double qq,int N1,int N2,double A[],double BB[],double XS[][5])
{
	double Y2[2],RES;//Y2[2]对应于表中包含给定马赫数或攻角值的最小区间端点的函数值
	int I,head,end; //I用于存储表中最靠近且小于给定攻角的攻角在攻角数组中的下标
				    //head用于二分查找法时的首下标，end为尾下标
	/*利用二分法获得表中包含给定攻角的最小区间的攻角元素*/
	//先对列进行差值
	head=1;end=N2;
	if(qq>BB[N2-1]) I=end-2;
	else
	{
		while(end-head!=1)
		{
			I=(head+end)/2;
			if(qq<BB[I-1]) end=I;
			else head=I;
		}
		I=head-1;
	}

	/*这里将二元插值降维，转换成先对最小区间的两端点处的马赫数进行插值，获取其气动系数*/
	Y2[0]=int31(A,XS,N1,N2,I,x);
	Y2[1]=int31(A,XS,N1,N2,I+1,x);
	/*利用以上获得的气动系数及两点式线性插值公式，插得给定攻角下的气动系数*/
	RES=((qq-BB[I])*Y2[1]-(qq-BB[I+1])*Y2[0])/(BB[I+1]-BB[I]);
    
	return(RES);

}
/*************************************************************************************/
/*                       等距双变元抛物线插值子函数                                  */
/*                                 int31()                                           */
/*   x为要差得行值，N1为行数，N2为列数，A为行向量，XS为二维表，I为主差值函数的列值   */          
/*************************************************************************************/
double int31(double A[],double XS[][5],int N1,int N2,int I,double x)
{
	int i,j,k,m;
	double z,s;
	z=0.0;
	if(N1<1) return(z);                                 //若马赫数无插值项则返回0
	if(N1==1) {z=XS[0][I];return(z);}                   //只有一项时直接赋值
	if(N1==2)                                           //只有两项时用线性插值
	{
		z=(XS[0][I]*(x-A[1])-XS[1][I]*(x-A[0]))/(A[0]-A[1]);
		return(z);
	}
	
	/*超过两项时使用一元三点抛物线插值*/
	if(x<A[1]) {k=0;m=2;}
	else if(x>=A[N1-2]) {k=N1-3;m=N1-1;}
	else
	{
		/*二分法找区间*/
		k=1,m=N1;
		while(m-k!=1)
		{
			i=(k+m)/2;
			if(x<A[i-1])m=i;
			else k=i;
		}
		k=k-1;m=m-1;
		if(fabs(x-A[k])<fabs(x-A[m])) k=k-1;
		else m=m+1;
	}
	z=0.0;
	for(i=k;i<=m;i++)
	{
		s=1.0;
		for(j=k;j<=m;j++)
			if(j!=i) s=s*(x-A[j])/(A[i]-A[j]);
		z=z+s*XS[i][I];
	}
	return(z);
}
/*************************************************************************************/
/*                             不等距单变元差值函数                                  */
/*                                 int11()                                           */
/*                         x为要差得值，n为列数，yy为一维表                          */          
/*************************************************************************************/
double int11(double yy[2][5],int n, double x)
{
	int i,j;
	double fy;
	for(j=0;j<=n-2;j++)
	{
		if(x<=yy[0][j+1])
		{
			i=j;
			break;
		}
		else
			i=n-2;
	}
	fy=yy[1][i]+(yy[1][i+1]-yy[1][i])*(x-yy[0][i])/(yy[0][i+1]-yy[0][i]);
	return(fy);
}
/*************************************************************************************/
/*                             输出数据文件子函数                                    */
/*                                 result()                                          */
/*************************************************************************************/
void result(double ALPHA)
{
	fprintf(fo,"%5.2f   ",Y[0]);                               //时间
	fprintf(fo,"%5.4f   ",MA);                                 //马赫数
	fprintf(fo,"%5.3f   ",ALPHA);                              //攻角
	fprintf(fo,"%5.2f   ",Y[3]);                               //射程
	fprintf(fo,"%5.3f   ",dX);                                 //水平速度
	fprintf(fo,"%5.3f   ",Y[4]);                               //高度
    fprintf(fo,"%5.3f   ",dY);                                 //垂直速度
    fprintf(fo,"%5.3f   ",Y[1]);                               //合速度
	fprintf(fo,"%5.3f   ",Y[2]*57.3);                          //弹道倾角
	fprintf(fo,"%5.3f   ",ALPHA+Y[2]*57.3);                    //俯仰角
//	fprintf(fo,"%5.3f   ",dTheta);                             //弹道倾角角速率
//	fprintf(fo,"%5.3f   ",guozai);                             //需用过载
	fprintf(fo,"\n");
	return;
}