通用弹道仿真计算程序
float rho_p=1600;
float w=1.16;
float power=950000;
float theta=0.25;
float alpha=0.001;
float u1=1.95e-8;
float index=0.83;
float e1=0.00055;
float d=0.00055;
float chi=0.75;
float lambda=0.12;
float mu=0.0;
float chi_s=1.696;
float lambda_s=-0.4104;
float m=2.8;
float s=0.00266;
float V0=0.00151;
float l_g=3.624;
float phi=1.168;
float l_psi;
float l_0;
float z_k;
float psi;
float psi_0;
float p;
float p1,p2,t1,t2,l1,l2;
int bijiao=0; /*以上为全局变量，内弹道参数值*/ void rk(n,t0,t,min,y0,y,work0,work1,work2,work3,fkt)/*变步长龙格库塔法子程序*/
int n;
float t0,t,min,y0[]; void fkt();/*右端函数*/
float (*y)[];
float work0[],work1[],work2[],work3[];
{
float temp[5];
int m,i,j,k,exit;
float step,t1,ttemp,tempd,c;
m=1;
exit=0;
while(exit==0)
{
step=(t-t0)/m;
for(i=1;i<=n;i++)
(*y)[i-1]=y0[i-1];
t1=t0;
temp[0]=0.5*step;
temp[1]=temp[0];
temp[2]=temp[0];
temp[3]=step;
temp[4]=step;
for(k=1;k<=m;k++)
{
ttemp=t1;
for(i=1;i<=n;i++)
{
work1[i-1]=(*y)[i-1];
work3[i-1]=(*y)[i-1];
}
for(j=1;j<=4;j++)
{
fkt(work1,work2);
t1=ttemp+temp[j];
for(i=1;i<=n;i++)
{
work1[i-1]=work3[i-1]+temp[j-1]*work2[i-1];
(*y)[i-1]=(*y)[i-1]+temp[j]*work2[i-1]/3;
}
}
}
m=m+m;
if(m!=2)
{
tempd=0;
for(i=1;i<=n;i++)
{
if((*y)[i-1]>1)
c=((*y)[i-1]-work0[i-1])/(*y)[i-1];
else
c=(*y)[i-1]-work0[i-1];
if(c>tempd)
tempd=c;
}
if(tempd<=min)exit=1;
}
if(exit==0)
{
for(i=1;i<=n;i++)
work0[i-1]=(*y)[i-1];
}
}
}
#include
#include
void fkt(float(*y)[],float(*f)[]) {
z_k=1+0.2596*(d*0.5+e1)/e1;
l_0=V0/s;
l_psi=l_0*(1-(w/V0)/rho_p-(w/V0)*(alpha-1/rho_p)*psi);
if((*y)[0]<1)
{
psi=chi*(*y)[0]*(1+lambda*(*y)[0]+mu*(*y)[0]*(*y)[0]);
}
else if((*y)[0]=l_g)break;
fprintf(ftxtfile,"%f %f %f %f\n",t,y[1],y[2],p);
for(i=1;i<=3;i++)
y0[i-1]=y[i-1];
t0=t;
}
fclose(ftxtfile);
}
浙公网安备 33021202002483