弹道仿真C源程序.doc

弹道仿真C源程序 /*铅垂平面内飞行器弹道仿真C源程序*/#include "math.h"#include "stdio.h"#include "stdlib.h"#define RAD 57.3double ACX56,ACY56,AMZAF56,AMZWZ52,AXG112,AJZ32,AP2,AMC2,AGC2,ANDM2,ANDAF2,Y8,B3,L,S,SONIC,RHO;/* ACX56:阻力系数数组；ACY56：升力系数数组；*/* AMZAF56：某一质心位置下的静稳定性导数数组；AMZWF52：阻尼力矩系数导数数组；*/* AXG102：质心位置Xg随时间的变化规律； AJZ32：转动惯量Jz随时间的变化规律；*/* AP2:起飞，续航发动机的推力值； AMC2：起飞，续航发动机燃料质量秒流量值； */* AGC2：质心位置变化的始末值； ANDM2：气动数据插值所需要的Ma数的最大最小值； */* ANDAF2：气动数据插值所需要的攻角的最大最小值 ； Y8存放积分结果的数组； */* B3 ：存放三个时间值的数组；L：特征长度；S：特征面积；SONIC：声速；RHO：大气密度；*/double MA,ABC_ALPHA,CX,CY,MZAF,MZWZ,XG,JZ,ALPHA;/* MA:马赫数；ABC_ALPHA:攻角的绝对值；CX:阻力系数；CY：升力系数；MZAF:静稳定性导数值；*/* MZWZ：阻尼力矩导数值； XG：质心位置；JZ：转动惯量；ALPHA：攻角；*/FILE *fi,*fo;/* 主函数*/void main() void rk(); /* 调用rk函数*/void result(); /* 调用result函数*/int i,j;double h; /* h为步长*/ fi=fopen("input.dat","r");for(i=0;i<5;i+) for(j=0;j<6;j+)fscanf(fi,"%lf",&ACXij);for(i=0;i<5;i+) for(j=0;j<6;j+)fscanf(fi,"%lf",&ACYij);for(i=0;i<5;i+) for(j=0;j<6;j+)fscanf(fi,"%lf",&AMZAFij);for(i=0;i<5;i+) for(j=0;j<2;j+)fscanf(fi,"%lf",&AMZWZij);for(i=0;i<10;i+) for(j=0;j<2;j+)fscanf(fi,"%lf",&AXGij);for(i=0;i<3;i+) for(j=0;j<2;j+)fscanf(fi,"%lf",&AJZij);for(i=0;i<2;i+) fscanf(fi,"%lf",&APi);for(i=0;i<2;i+)fscanf(fi,"%lf",&AMCi);for(i=0;i<2;i+)fscanf(fi,"%lf",&AGCi);for(i=0;i<2;i+)fscanf(fi,"%lf",&ANDMi);for(i=0;i<2;i+)fscanf(fi,"%lf",&ANDAFi);for(i=0;i<8;i+)fscanf(fi,"%lf",&Yi);for(i=0;i<3;i+)fscanf(fi,"%lf",&Bi);fscanf(fi,"%lf%lf%lf%lf%lf",&L,&S,&SONIC,&RHO,&h);fclose(fi);fo=fopen("output.dat","w");fprintf(fo,"%s","TIME ALPHA V X Y ");fprintf(fo,"%s"," MASS THETA GTHETA WZ");fprintf(fo,"n");fprintf(fo,"%s"," s rad m/s m m ");fprintf(fo,"%s"," kg rad rad rad/s");fprintf(fo,"n");rk();doresult();rk();while(Y6>=0);/* */* 方程 右端子函数 */void dery(n,dy,Y)int n;double dy, Y;extern void interp();/* 调用插值子函数 */double q,a4,XF,YF,P,wzt,MC;q=RHO*Y1*Y1/2; /* q为动压 */MA=Y1/SONIC;ALPHA=Y4-Y2;ABC_ALPHA=fabs(ALPHA);interp();a0=sin(Y2);/*Y0是时间t,Y1是速度v，Y2是弹道倾角，Y3是转动角速度；*/a1=cos(Y2);/*Y4是俯仰角；Y5是x;Y6是y;Y7是质量；*/a2=sin(ALPHA);a3=cos(ALPHA);if (ALPHA<0)CY=-CY;XF=CX*q*S;YF=CY*q*S;wzt=Y3*L/Y1;MC=Y0<B1?AMC0:Y0<B2?AMC1:0;P=Y0<B1?(AP0+AP1):Y0<B2?AP1:0;/*选择发动机的推力*/dy0=1;dy1=(P*a3-XF-9.81*Y7*a0)/Y7;dy2=(P*a2+YF-9.81*Y7*a1)/(Y1*Y7);dy3=(MZAF*ALPHA+MZWZ*wzt)*q*S*L/JZ;dy4=Y3;dy5=Y1*a1;dy6=Y1*a0;dy7=-MC;return;/* * */* 四阶龙格库塔子函数*/void rk(n,h)int n;double h;extern void dery();double a4,old_y8,Y18,*dy;int i,j;dy=calloc(n,sizeof(double);a0=a1=h/2;a2=a3=h;dery(n,dy,Y);for(i=0;i<n;i+)old_yi=Yi;for(j=0;j<3;j+)for(i=0;i<n;i+)Y1i=old_yi+aj*dyi;Yi=Yi+aj+1*dyi/3;dery(n,dy,Y1);for(i=0;i<n;i+)Yi=Yi+a0*dyi/3;free(dy);return;/* * */* 插值子函数*/void interp()extern double int11(),int32();double mzaf0;if(Y0<B2)XG=int11(AXG,11,Y0);JZ=int11(AJZ,3,Y0);elseXG=AXG101;JZ=AJZ21;CX=int32(MA,ABS_ALPHA,5,6,ANDM,ANDAF,ACX);CY=int32(MA,ABS_ALPHA,5,6,ANDM,ANDAF,ACY);mzaf0=int32(MA,ABS_ALPHA,5,6,ANDM,ANDAF,AMZAF);if (ABS_ALPHA!=0.)MZAF=mzaf0*RAD+CY*(XG-AGC0)/(ABS_ALPHA*L);elseMZAF=mzaf0*RAD;MZWZ=int32(MA,XG,5,2,ANDM,AGC,AMZWZ);return;/* * */* 不等距单变元线性插值子函数*/double int11(double YY2,int n,double x)int i,j;double fy;for(j=0;j<=n-2;j+) if(x<=YYj+10)i=j;break;elsei=n-2;fy=YYi1+(YYi+11-YYi1)*(x-YYi0)/(YYi+10-YYi0);return(fy);/* * */* 等距单变元抛物线插值子函数 */double int31(x,N1,N2,I,A,Y0)double x,A2;int N1,N2,I;int I1;double H,AI,YE,x0,x1,x2,c0,c1,c2;H=(A1-A0)/(N1-1);I1=(x-A0)/H+1;if(I1-1)<0)I1=1;elseif(I1-(N1-2)>=0)I1=N1-2;x0=A0+(I1-1)*H;x1=x0+H;X2=X0+2*H ;c0=0.5*(x-x1)*(x-x2)/(H*H);c1=(x-x0)*(x-x2)/(H*H);c2=0.5*(x-x0)*(x-x1)/(H*H);YE=c0*Y0I1I+c1*Y0I1+1I+c2*Y0I1+2I;return(YE);/* * */* 等距双变元抛物线线性插值子函数 */double int32(x,qq,N1,N2,A,BB,Y0)double x,qq,A2,BB2;int N1,N2;extern double int31();int I;double H,Y22,RES;H=(BB1-BB0)/(N2-1);I=(qq-BB0)/H+1;if(I-1)<0)I=1;elseif(I-(N2-1)>=0)I=N2-1;I=I-1;Y20=int31(x,N1,N2,I,A,Y0);Y21=int31(x,N1,N2,I+1,A,Y0);RES=(qq-I*H-BB0)*Y21-(qq-(I+1)*H-BB0)*Y20)/H;return(RES);/* * */* 输出数据文件子函数 */void result()static int i=1;fprintf(fo,"%-7.3f",Y0);fprintf(fo,"%6.3f",ALPHA);fprintf(fo,"%8.3f",Y1);fprintf(fo,"%8.3f",Y5);fprintf