function [V]=twosphv2(ia,ib,d,c);
%V = twosphv2(ia,ib,d,c);
%ia, Inner iadius of A
%ib, Inner iadius of B
%d, Thickness of the wall/shell
%c, Difference between centers
%
%Made by CA for SBE
%
%Corrected by CA, April 1998, gave twice the volume!!

A=0;
oa=ia+d;
ob=ib+d;

format long
format compact

pl=1;
if pl==1,
   x=[floor(-oa):((ib+c+ia)/1000):ceil(ob+c)];
   ya1=sqrt(ia.^2-x.^2);
   t=find(imag(ya1)~=0);ya1(t)=NaN.*ones(size(t));
   
   yb1=sqrt(ib.^2-(x-c).^2);
   t=find(imag(yb1)~=0);yb1(t)=NaN.*ones(size(t));
   
   ya2=sqrt(oa.^2-x.^2);
   t=find(imag(ya2)~=0);ya2(t)=NaN.*ones(size(t));
   
   yb2=sqrt(ob.^2-(x-c).^2);
   t=find(imag(yb2)~=0);yb2(t)=NaN.*ones(size(t));
   
   plot(x,abs(ya1)),hold on,
   plot(x,abs(ya2))
   plot(x,abs(yb1))
   plot(x,abs(yb2))
end;

p1=-(ob^2-oa^2-c^2)/(2*c);q1=sqrt(oa^2-p1^2);
p2=-(ib^2-oa^2-c^2)/(2*c);q2=sqrt(oa^2-p2^2);
p3=-(ob^2-ia^2-c^2)/(2*c);q3=sqrt(ia^2-p3^2);
p4=-(ib^2-ia^2-c^2)/(2*c);q4=sqrt(ia^2-p4^2);

if pl==1,
   [p1 q1],plot(real(p1),real(q1),'o');text(real(p1),real(q1),'1');
   [p2 q2],plot(real(p2),real(q2),'o');text(real(p2),real(q2),'2');
   [p3 q3],plot(real(p3),real(q3),'o');text(real(p3),real(q3),'3');
   [p4 q4],plot(real(p4),real(q4),'o');text(real(p4),real(q4),'4');
   set(gca,'DataAspectRatio',[1 1 1])
end;

if -oa<=(c-ib) & (c-ib)<=-ia,
   A1=intv1([-oa p2],oa);
   A2=intv2([c-ib p2],ib,c);
   A=A1-A2;
   if pl==1,
      plot([-oa p2],[0 0],'b+');
      plot([c-ib p2],[0 0],'ro');
      xlabel('Case 1');
   end;
end; %Tested, OK

if -oa<=(c-ob) & (c-ob)<=-ia & real(p4)<0,
   A1=intv2([c-ob p1],ob,c);
   A2=intv1([p1 p2],oa);
   A3=intv1([-ia p4],ia);
   A4=intv2([p4 p2],ib,c);
   A=A1+A2-A3-A4;
   if pl==1,
      plot([c-ob p1],[0 0],'b+');
      plot([p1 p2],[0 0],'b+');
      plot([-ia p4],[0 0],'ro');
      plot([p4 p2],[0 0],'ro');
      xlabel('Case 2');
   end;
end; %Tested, OK

if -ia<c-ob & c-ib<ia,
   A1=intv2([p3 p1],ob,c);
   A2=intv1([p1 p2],oa);
   A3=intv1([p3 p4],ia);
   A4=intv2([p4 p2],ib,c);
   A=A1+A2-A3-A4;
   if pl==1,
      plot([p3 p1],[0 0],'b+');
      plot([p1 p2],[0 0],'b+');
      plot([p3 p4],[0 0],'ro');
      plot([p4 p2],[0 0],'ro');
      xlabel('Case 3');
   end;
end; %Tested, OK 

if ia<=(c-ib) & (c-ib)<=oa,
   A1=intv2([p3 p1],ob,c);
   A2=intv1([p1 p2],oa);
   A3=intv1([p3 ia],ia);
   A4=intv2([c-ib p2],ib,c);
   A=A1+A2-A3-A4;
   if pl==1
      plot([p3 p1],[0 0],'b+');
      plot([p1 p2],[0 0],'b+');
      plot([p3 ia],[0 0],'ro');
      plot([c-ib p2],[0 0],'ro');
      xlabel('Case 4');
   end;
end; 

if oa<=(c-ib) & (c-ob)<=oa,
   A1=intv2([c-ob p1],ob,c);
   A2=intv1([p1 oa],oa);
   A=A1+A2;
   if pl==1,
      plot([c-ob p1],[0 0],'b+');
      plot([p1 oa],[0 0],'b+');
      xlabel('Case 5');
      
   end;
end; 

V=A;
if pl==1,
   axis([min(x) max(x) 0 max([oa ob])]);
   title(['ia,oa,ib,ob,c = (' num2str(ia) ',' num2str(oa) ',' num2str(ib) ',' num2str(ob) ',' num2str(c) ')' ])
   ylabel(['V= ' num2str(V)]);
   grid off
   zoom on
   drawnow
   hold off
end;

