Matlab绘制箕舌线

a=input('a=?');
x=-4*a:0.01:4*a;
y=8.*a.^3./(x.^2+4.*a.^2);
b=input('b=?');
c=input('c=?');
bclcm=lcm(b,c);
th=-2*pi*bclcm:0.02:2*pi*bclcm;
rou = zeros(size(th));
for i = 1:length(rou)
    rou(i)=0.6*a+0.4*a*sin(b*th(i)/c);
end
thr = -2*pi:0.01:2*pi;
plot(x,y,'b',rou.*cos(th),a+rou.*sin(th),'y',a.*cos(thr),a+a.*sin(thr),'m',x,0,'g');%,'LineWidth','4'
set(findobj(get(gca,'Children'),'Color','b'),'LineWidth',2);
title('versoria');
axis('equal');