% Two long parallel cylinders
function Fij = ParalCylView(ri,rj,s)
  R=rj./ri; S=s./ri; C=1+R+S;
  Fij = (pi+sqrt(C.^2-(R+1).^2)-sqrt(C.^2-(R-1).^2)+ ...
        (R-1).*acos((R-1)./C)-(R+1).*acos((R+1)./C))/2/pi;
end
