内弹道程序

#include"stdio.h"
#include"math.h"
#include"stdlib.h"
#define roup 1600
#define w 14.72
#define f 950000
#define sita 0.2
#define alpha 0.001
#define u1 1.14*10^-8
#define n1 0.8275
#define e1 0.000893
#define x1 1.06
#define lamada -0.0566
#define miu 0.0
#define m 46
#define s 0.01905
#define V0 0.020027
#define lg 7.3206
#define fai 1.12667
#define P0 3*10^7
#define deda 735
#define L0 1.0513
#define Vj 1642.38
#define B 3.4559
#define lgpingjun 6.96
double Y[6],Z[4];
FILE *fo;
int main()//主函数
{
Y[0]=0,Y[1]=0,Y[2]=0.04296,Y[3]=0,Y[4]=0.0223,Y[5]=0.021,Z[0]=0,Z[1]=0,Z[2]=3*pow(10,7
),Z[3]=0; //依次分别给相对时间，相对速度,相对压强，炮弹相对位移,，相对燃烧量
//相对燃烧厚度,时间，速度,压强，炮弹位移。

void rk(int n,double h);
void result();
fo=fopen("output.txt","w");
fprintf(fo,"%s"," 相对时间绝对时间炮弹位移速度压强相对燃烧量相对燃烧厚度\n");
fprintf(fo," s m m/s pa \n");
do{
result();
rk(6,0.001);
}
while (Z[3]<lg);//Z[3]表示炮弹位移，截止条件
fclose(fo);
}
//龙格库塔算法
void rk(int n,double h)
{
externvoid dery(int n,double dy[],double Y[]);
double a[4],old_Y[6],*dy,Y1[6];
int i,j;
dy =(double*) calloc(n, sizeof(double));
a[0]=a[1]=h/2;
a[2]=a[3]=h;
dery(n,dy,Y);
for(i=0;i<n;i++)
old_Y[i]=Y[i] ;
for(j=0;j<3;j++)
{
for(i=0;i<n;i++)
{
Y1[i]=old_Y[i]+a[j]*dy[i];
Y[i]=Y[i]+a[j+1]*dy[i]/3;
}
dery(n,dy,Y1);
}
for(i=0;i<n;i++)
Y[i]=Y[i]+a[0]*dy[i]/3;
free(dy);
return;
}
//右端子式
void dery(int n,double dy[],double Y[])
{
dy[0]=1;
if(Y[5]<1 )
{
dy[4]=x1*(1+2*lamada*Y[5]+3*miu*pow(Y[5],2))*sqrt(sita/(2*B))*pow(Y[2],n1) ;
}
elseif(Y[5]>=1)
{
dy[4]=0;
}
if(Y[5]<1)
{
dy[5]=sqrt(sita/(2*B))*pow(Y[2],n1) ;
}
else
{
dy[5]=0;
}
dy[3]=Y[1];
dy[1]=sita/2*Y[2];
dy[2]=1/(Y[3]+0.541-0.2756*Y[4])*(1+deda*(alpha-1/1600)*Y[2])*dy[4]-(1+sita)/(Y[3]+0.5 41-0.2756*Y[4])*Y[2]*Y[1];
Z[0]=L0/Vj*Y[0];
Z[1]=Vj*Y[1];
Z[2]=f*deda*Y[2];
Z[3]=L0*Y[3];
return;
}
//显示结果
void result()
{
staticint i=1;
{
fprintf(fo,"%13.3f",Y[0]);//相对时间
fprintf(fo,"%13.6f",Z[0]);//绝对时间
fprintf(fo,"%13.6f",Z[3]);//炮弹位移
fprintf(fo,"%13.3f",Z[1]);//速度
fprintf(fo,"%20.3f",Z[2]);//压力
fprintf(fo,"%13.6f",Y[4]);//相对燃烧量
fprintf(fo,"%13.6f",Y[5]);//相对燃烧厚度
fprintf(fo,"\n");
i++;
}
return;
}。