文档介绍：通用弹道仿真计算程序
float rho_p=1600;
float w=;
float power=950000;
float theta=;
float alpha=;
float u1=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+*(d*+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]<z_k)
{
psi=chi_s*(*y)[0]*(1+lambda_s*(*y)[0]/z_k)/z_k;
}
else
{
psi=1;
}
p=(power*w*psi-*theta*phi*m*(*y)[2]*(*y)[2])/((*y)[1]+l_psi)/s;
if((*y)[0]<z_k)
{
(*f)[0]=u1*pow(p,index)/e1;
}
else
{
(*f)[0]=0;
}
(*f)[1]=(*y)[2];
(*f)[2]=s*p/(phi*m);
return;
}
main()
{
float y0[3];
int i,k,n;
float y[3],work0[3],work1[3],work2[3],work3[3];
float t0,t,min;
FILE*ftxtfile=NULL;
void fkt();
psi_0=(V0/w-1