function [r,r2]=ratio(a,b,s1,s2);
%RATIO compute ratio of hyperbolic functions
%      to allow extreme variations of arguments.
%
% CALL: r=ratio(a,b,s1,s2);
%
%       r = f(a,s1)./f(b,s2), ratio vector hyperbolic functions of same size as a and b
%     a,b = arguments vectors of the same size
%   s1,s2 = defines the hyperbolic function used, i.e.,
%           f(x,1)=cosh(x), f(x,-1)=sinh(x)
%
% Examples:
% ratio(2,1,1,1)   % gives r=cosh(2)/cosh(1)
% ratio(2,1,1,-1)  % gives r=cosh(2)/sinh(1)
% ratio(2,1,-1,1)  % gives r=sinh(2)/cosh(1)
% ratio(2,1,-1,-1) % gives r=sinh(2)/sinh(1)
%
% See also: tran

% Tested on: matlab 5.2
% history
% revised pab 07.11.2001
% -added comnsize + see also line
% -Fixed a bug: ratio(0,0,-1,-1) gave NaN but should return 1
% revised pab 11.01.2000
% - added sign(s1), sign(s2) to ensure correct calculation
% - updated documentation
% - fixed a bug in expression.
% by L. Borgman

[ec, a,b,s1,s2]=comnsize(a,b,sign(s1),sign(s2));
if ec~=0,
   error('a,b,s1 and s2 must be of common size or scalar!')
end	
r = zeros(size(a));

sc = (a==b) & (s1 == s2) & (s2==-1);
 
k = find(sc );
if any(k)
   r(k) = 1;
end	


k1 = find(~sc);
if any(k1),
   ak = a(k1);
   bk = b(k1);
   sa = s1(k1);
   sb = s2(k1);
   warning off    % fix to avoid warning messages about division by zero.    
   r(k1)=exp(ak-bk).*(1+sa*exp(-2*ak))./(1+sb*exp(-2*bk));
   warning on
 end	

return

% Old call
r=exp(a-b).*(1+sign(s1)*exp(-2*a))./(1+sign(s2)*exp(-2*b));

return