function [Cx Cy SZB Mz]=mianyuan(xn,on,al,Ma)
a=0.1;
b=0.5;
D=3.192;
L=4.2*D;
gam=1.4;
x=-0:(-L/xn):-L;
S=pi*(L/4*tan(16.5*pi/180)+D/2)^2;
for(i=1:xn+1)
if (x(i)>-a)
r(i)=sqrt(abs((1-(x(i)+a)^2/a^2)*b^2));
else if(x(i)>-L/4)
r(i)=b-(x(i)+a)*(D/2-b)/(0.25*L-a);
else if(x(i)>-(1-0.25)*L)
r(i)=D/2;
else
r(i)=D/2-(x(i)+(1-0.25)*L)*tan(16.5*pi/180);
end
end
end
end
o=0:2*pi/on:2*pi;
hold on;
for(i=1:xn+1)
for(j=1:on)
z(i,j)=r(i)*sin(o(j));
y(i,j)=r(i)*cos(o(j));
end
end
X=cell(xn,on);Y=cell(xn,on);Z=cell(xn,on);nn=cell(xn,on);
v=[-299.46*Ma*cos(al/180*pi),0,299.46*Ma*sin(al/180*pi)];
y(:,on+1)=y(:,1);z(:,on+1)=z(:,1);
for i=1:xn
for j=1:on
X{i,j}=[x(i),x(i),x(i+1),x(i+1)];
X0(i,j)=(x(i)+x(i)+x(i+1)+x(i+1))/4;
Y{i,j}=[y(i,j),y(i,j+1),y(i+1,j+1),y(i+1,j)];
Y0(i,j)=(y(i,j)+y(i,j+1)+y(i+1,j+1)+y(i+1,j));
Z{i,j}=[z(i,j),z(i,j+1),z(i+1,j+1),z(i+1,j)];
Z0(i,j)=(z(i,j)+z(i,j+1)+z(i+1,j+1)+z(i+1,j));
T1=[(X{i,j}(3)-X{i,j}(1)),(Y{i,j}(3)-Y{i,j}(1)),(Z{i,j}(3)-Z{i,j}(1))];
T2=[(X{i,j}(4)-X{i,j}(2)),(Y{i,j}(4)-Y{i,j}(2)),(Z{i,j}(4)-Z{i,j}(2))];
NN=cross(T1,T2);
Si(i,j)=sqrt(dot(NN,NN));
nn{i,j}=NN/sqrt(dot(NN,NN));
theta(i,j)=pi/2-acos(dot(nn{i,j},v)/sqrt(dot(v,v)));
end
end
for i=1:xn
for j=1:on
if theta(i,j)>=0
Cp(i,j)=2/(gam*Ma^2)*(((gam+1)/2*Ma^2)^(gam/(gam-1))*1/((2*gam/(gam+1)*Ma^2-(gam-1)/(gam+1))^(1/(gam-1)))-1)*(sin(theta(i,j)))^2;
else Cp(i,j)=0;
end
Cxti(i,j)=Cp(i,j)*nn{i,j}(1)*Si(i,j)/S; %%ij点阻力系数
Cyti(i,j)=Cp(i,j)*nn{i,j}(3)*Si(i,j)/S; %%ij点升力系数
Czti(i,j)=Cp(i,j)*nn{i,j}(2)*Si(i,j)/S; %%ij点侧向力系数
Mzi(i,j)=(-nn{i,j}(3)*X0(i,j)+nn{i,j}(1)*Z0(i,j))*Cp(i,j)*Si(i,j)/(L*S);%%ij点力矩系数
end
end
Cxt=sum(sum(Cxti));%%弹体阻力系数
Cyt=sum(sum(Cyti));%%弹体升力系数
Czt=sum(sum(Czti));%%弹体侧向力系数
Mz=sum(sum(Mzi)); %%弹体力矩系数
Cy=Cxt*sin(al/180*pi)+Cyt*cos(al/180*pi);%%升力系数
Cx=-Cxt*cos(al/180*pi)+Cyt*sin(al/180*pi);%%阻力系数
SZB=(Cy/Cx);%%升阻比
end
