function [time,delta,lowerI,upperI]=tau(w0,l,r,cs,wcs);
% MATLAB HELP:
%
% function [time,delta,lowerI,upperI]=tau(w0,l,r,cs,wcs);
%
% inputs:
%  w0, frequency
%  l, angular degree
%  r, radius
%  cs, square of the sound speed as a function of radius
%  wcs, square of the acoustic cutoff of frequency as function of radius
% outputs:
%  time,   first bounce phase time
%  delta,  first bounce distane
%  lowerI, index into r of the lower turning point
%  upperI, index into r of the upper turning point
%
%  see ray_codes.ps for details
%
% aaron birch may 22 2002



N=1;

%
% solar radius in cm
%
R=6.9599e10;

%
% square of horizontal wavenumber
%
khs=l*(l+1)*r.^(-2);  
%
% square of radial wavenumber as function of radius
%
krs=(w0^2-wcs)./cs-khs;  

%
% find where radial wavenumber switches sign, i.e. turning points
%
left=sign(krs(1:(length(krs)-1)));
right=sign(krs(2:length(krs)));
I=find(left.*right==-1);

%
% approximate lower turning point
%
lower=(min(r(I))+min(r(I+1)))/2;  

%
% approximate upper turning point
%
upper=(max(r(I))+max(r(I+1)))/2; 

%
% is there are more than two turning points, then try
% to figure out what is going on
%
if (length(I)>2)
     disp('more than two turning points !');
     %
     % if two neighboring points are zero crossings then drop them
     %
     j=2;done=0;
     while (~done & j<=length(I))
     if (I(j)==I(j-1)+1)
     I=I([1:j-2 j+1:length(I)]);done=1;
     end;
     j=j+1;
     end;
     
     %
     % check is there still is a problem
     %
     if (length(I)>2)
     disp(' length(I) > 2, cant find turning points');
     w0,l,I
     end;
end;

lowerI=min(I)+1; %% really should be +1
upperI=max(I); %% and +0
range=(lowerI+N):(upperI-N);

%
%  now have to do some tricks to integrate over the singularities
%  at the turning points

%
% first treat the time at lower turning point 
%
lr=(lowerI-N):(lowerI+N);
fcn=cs(lr).^2.*krs(lr);
[P,S,mu]=polyfit(r(lr),fcn,1); % do linear fit
P(2)=P(2)-P(1)*mu(1)/mu(2);
P(1)=P(1)/mu(2);

ltime=4*w0*sqrt(P(1)*r(max(lr))+P(2))/P(1); %% and do the integration


%
% delta at the lower turning point
%
fcn=(r(lr).^2).*krs(lr)./khs(lr);
[P,S,mu]=polyfit(r(lr),fcn,1); % do linear fit
P(2)=P(2)-P(1)*mu(1)/mu(2);
P(1)=P(1)/mu(2);

ldelta=720*sqrt(P(1)*r(max(lr))+P(2))/P(1)/pi; %% and do the integration


%
% time at the upper turning point
%
lr=(upperI-N):(upperI+N);
fcn=cs(lr).^2.*krs(lr);
[P,S,mu]=polyfit(r(lr)/R,fcn,1); % do linear fit
P(2)=P(2)-P(1)*mu(1)/mu(2);
P(1)=P(1)/mu(2);

utime=-4*w0*R*sqrt(P(1)*r(min(lr))/R+P(2))/P(1); %% and do the integration


%
% delta at the upper turning point
%
fcn=(r(lr).^2).*krs(lr)./khs(lr);
[P,S,mu]=polyfit(r(lr)/R,fcn,1); % do linear fit
P(2)=P(2)-P(1)*mu(1)/mu(2);
P(1)=P(1)/mu(2);

udelta=-720*R*sqrt(P(1)*r(min(lr))/R+P(2))/P(1)/pi; %% and do the integration





%
% away from the turning point things are easier
%
time = integrate(r(range)/R,2*R*w0*(cs(range).*sqrt(krs(range))).^(-1));
delta= integrate(r(range),360*sqrt(khs(range))./sqrt(krs(range))./r(range)/pi);

%
% add in the corrections from the turning points
%
time=time+ltime+utime;
delta=delta+ldelta+udelta;


%
% all done
%