这是我自己写的代码,不知道哪里出问题了
syms r1 r2 rb1 rb2 a ra1 ra2 y Vy V1 V2 y02 ha1 ha2 d1 d2 da1 da2 q1 q2 vo uo u02 u02o x2 G Gs i;
m=3;z1=51;z2=52;z0=25;ha=0.8;u=20;x0=0;i=1;A=[];
v=u/180*pi;
r1=m*z1/2; r2=m*z2/2;rb1=r1*cos(v);rb2=r2*cos(v);a=(z1+z2)*m/2;ha1=ha*m; ha2=ha*m; ra2=r2-ha2; ra1=r1-ha1;
for uo=1:0.01:90
vo=uo/180*pi;
ao=a*cos(v)/cos(vo);
y=(ao-a)/m;
for u02o=1:0.01:90
v02o=u02o/180*pi;
V2=((tan(v)-v)-(tan(v02o)-v02o))*(z2-z0);
V1=((tan(v)-v)-(tan(vo)-vo))*(z2-z1)-V2;
x1=V1/(2*tan(v));
y02=((z2-z0)/2)*(cos(v)/cos(v02o)-1);
Vy=y-y02+x1;
if Vy>0;
G=ra2+ao-ra1;
else
continue;
end ;
if G>0;
d1=m*z1;d2=m*z2;da1=(z1+2*ha)*m;da2=(z2+2*ha)*m;
ua1=acos(d1*cos(u)/da1);ua2=acos(d2*cos(u)/da2);
va1=ua1*pi/180; va2=ua2*pi/180;
q1=acos((ra2^2-ra1^2-ao^2)/(2*ra1*ao));
q2=acos((ra2^2-ra1^2+ao^2)/(2*ra2*ao));
ca=(z1*(tan(va1)-tan(vo))-z2*(tan(va2)-tan(vo)))/(2*pi);%ca为重合度
else
continue;
end;
if abs(cos(q1))<=1 && abs(cos(q2))<=1 && ca>1;
Gs=z1*(q1+tan(va1)-va1)-z2*(q2+tan(va2)-va2)+(tan(uo)-uo)*(z2-z1);
else
continue;
end ;
if abs(Gs-0.05)<0.001;
A(i,1)=x1;A(i,2)=y02;A(i,3)=vo;A(i,4)=v02o;A(i,5)=a;A(i,6)=ao;A(i,7)=ca;A(i,8)=Gs;
i=i+1;
else
continue;
end ;
end
end