function im=LineForce(cr,a,phi)
% Calculate the stress on a cylinder of radius a with a constant line force 
% from -phi/2 to phi/2. 
% Usage: im=LineForce(r,a,phi)
%

% revision history:
% 08/04/05 Mark D. Shattuck <mds> ipf.m  
% 01/30/06 mds added abs(cr)
% 04/30/07 mds made w a true measure of width.
% 10/26/08 mds modify for PointForce function
% 10/27/08 mds modify for LineForce function

r=abs(cr);
th=angle(cr);
ii=abs(th)<=(pi-phi/2);
im=zeros(size(r));

if(numel(phi)>1)
  im(ii)=fp(r(ii),a,phi(ii)/2-th(ii))-fp(r(ii),a,-phi(ii)/2-th(ii));
  im(~ii)=2*pi-(-fp(r(~ii),a,phi(~ii)/2-th(~ii))+fp(r(~ii),a,-phi(~ii)/2-th(~ii)));
else
  im(ii)=fp(r(ii),a,phi/2-th(ii))-fp(r(ii),a,-phi/2-th(ii));
  im(~ii)=2*pi-(-fp(r(~ii),a,phi/2-th(~ii))+fp(r(~ii),a,-phi/2-th(~ii)));
end

im=(r<a).*im;

end

function F=fp(r,a,th)

%F=(a-r).*r.*(a+r).*sin(th)./(a.*(a.^2+r.^2-2*a*r.*cos(th)))+2*acot((a-r)./(a+r).*cot(th/2));
F=2*acot((a-r)./(a+r).*cot(th/2));

end