clear
p0=7;
p=12;
w1=30/pi;
w2=60/pi;
e=w2/w1;
g=p0/p;
hold on
for r=0:0.1:p+p0
r1=r./p;
if r>=p-p0 & r<=p+p0
c0=acos((g.^2.*(1.0./2.0)+r1.^2.*(1.0./2.0)-1.0./2.0)./(g.*r1));
elseif r<p-p0
c0=pi;
end
v=@(c0) p.*w1.*sqrt(r1.^2.*(e+1.0).^2+e.^2.*g.^2-e.*g.*r1.*cos(c0).*(e+1.0).*2.0);
z1=quad(v,-1*c0,c0);
plot(r,z1,'*')
grid on
end