内弹道程序

#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]

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

old_Y[i]=Y[i] ;

for(j=0;j<3;j++)

{

for(i=0;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

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;

}