内弹道程序

#define m 46.0
#define fi1 1.02
#define pm_tong 307720000.0
#define vg 965.0
#define f 950000.0
#define k 1.2
#define alpha 0.001
#define roup 1600.0
#define d 0.152
#define p0 30000000.0
#define sita 0.2
#define X 1.06
#define Xk 1.50
#define lamada -0.0566
#define S 0.01905
FILE *fo;
double i,j;
double w,V0,L0;
double fi,pm,rg;
double posai0,deta0,Z0,Ik;
double B,k1,deta1,posai,x,posaiba,B1,Lposaiba,L,Zx,Lposai;
double xk,r,b,Zk,Lk,Yk,Yg,km,Lg,Lv0,Lc,Lsh,Ytg,Ytk,Vg,P1,P2,Vnt,L1,Pk,Pg,Ntj;
int main()
{
void cacuZ0();
void cacuB();
void cacuYk();
void result();
fo=fopen("output.txt","w");
fprintf(fo,"%s"," 装填密度相对装药量w/m 弹丸相对行程长身管长能量利用效率炮膛工作容积利用效率火药相对燃烧结束位置炮口压力火炮寿命装填参量B 等药室容积等膛容压力冲
量\n");
fprintf(fo,"%s"," kg/m^3 1 1 m 1 1 1 pa 发
1 m^3/kg m^3/kg \n");
i=735.0;
j=0.32;
w=j*m;
V0=w/i;
fi=fi1+1.0/3.0*j;
L0=V0/S;
pm=1.12*pm_tong *fi;
rg=vg*vg*sita/(2*f*j);
cacuZ0();
cacuB();
cacuYk();
km=(Yk+1-alpha*i)*pow((1-B*(k-1)/2*(1-Z0)*(1-Z0)),1/(k-1));
Yg=km/pow(1-fi*rg,1/(k-1))+alpha*i-1;
Lg=Yg*L0;
Lv0=L0/Xk;
Lc=2*d;
Lsh=Lg+Lv0+Lc;
Vg=V0*Yg;
Ytg=fi*m*vg*vg/(2*Vg*pm);
Ytk=Lk/Lg;
L1=L0*(1-alpha*i);
Pk=f*w/S*(1-(B*sita/2)*(1-Z0)*(1-Z0))/(L1+Lk);
Pg=Pk*pow((L1+Lk)/(L1+Lg),1+sita);
Ntj=200*(Yg+1)/j;
P1=V0/m;
Vnt=V0*(1+Yg);
P2=Vnt/m;
Ik=pow(B*f*w*m*fi/S/S,0.5);
result();
}
void cacuZ0()
{
posai0=(1/i-1/roup)/(f/p0+alpha-1/roup);
deta0=pow(1+4*lamada*posai0/X,0.5);
Z0=2*posai0/(X*(1+deta0));
return;
}
void cacuB()
{
double x1, x2, y1, x, y;
double fun(double B) ;
x1 =0.3 ; x2 = 10.0;
do {
x = (x1 + x2)/2;
y1 = fun(x1);
y = fun(x);
if( ( y < 0 && y1 < 0) || (y > 0 && y1 > 0) )
x1 = x;
else
x2 = x;
/*end if*/
}while(fabs(y) > 0.001);
B=x;
return;
}
double fun(double B)
{
double xm,B2;
k1=X*deta0;
deta1=1/(alpha-1/roup);
xm=k1/(B*(1+sita)/(1+(pm/(f*deta1)))-2*X*lamada);
posai=posai0+k1*xm+X*lamada*xm*xm;
posaiba=(posai0+posai)/2;
B1=B*sita/2-X*lamada;
Lposaiba=L0*(1-i/roup-i*(alpha-1/roup)*posaiba);
r=B1*posai0/(k1*k1);
b=pow(1+4*r,0.5);
Zx=pow((1-2*B1*xm/(b+1)/k1),(b+1)/(2*b))*pow((1+2*B1*xm/(b-1)/k1),((b-1)/(2*b)));
L=Lposaiba*(pow(Zx,(-B/B1))-1);
Lposai=L0*(1-i/roup-i*(alpha-1/roup)*posai);
B2=(posai-pm*S*(L+Lposai)/(f*w))/(xm*xm*sita/2)-B;
return B2;
}
void cacuYk()
{
xk=1-Z0;
Zk=pow((1-2*B1*xk/(b+1)/k1),(b+1)/(2*b))*pow((1+2*B1*xk/(b-1)/k1),((b-1)/(2*b)));
Lk=Lposaiba*(pow(Zk,(-B/B1))-1);
Yk=Lk/L0;
return;
}
void result()
{
staticint h=1;
{
fprintf(fo,"%13.6f",i);//装填密度
fprintf(fo,"%13.6f",j);//相对装药量w/m
fprintf(fo,"%13.6f",Yg);// 弹丸相对行程长
fprintf(fo,"%20.6f",Lsh);//身管长
fprintf(fo,"%20.6f",rg);//能量利用效率
fprintf(fo,"%20.6f",Ytg);//炮膛工作容积利用效率
fprintf(fo,"%20.6f",Ytk);//火药相对燃烧结束位置
fprintf(fo,"%20.6f",Pg);//炮口压力
fprintf(fo,"%20.6f",Ntj);// 火炮寿命
fprintf(fo,"%20.6f",B);//装填参量B
fprintf(fo,"%20.6f",P1);// 等药室容积
fprintf(fo,"%20.6f",P2);//等膛容
fprintf(fo,"%20.6f",Ik);//压力冲量
fprintf(fo,"%20.6f",L0);//压力冲量
fprintf(fo,"%20.6f",Lg);//压力冲量
fprintf(fo,"\n");
h++;
}
return;
}。