function [r_path,theta_path,dtau_group,dtau_phase,path_length]=raypath(w0,l,r,cs,wcs);
% MATLAB HELP:
%
% function [r_path,theta_path,dtau_group,dtau_phase,path_length]=raypath(w0,l,r,cs,wcs);
%
% inputs:
%  w0, frequency of the ray
%  l, angular degree
%  r, grid in radius
%  cs, square of the sound speed as function of radius
%  wcs, square of the acoustic cutoff frequency as function of radius
% outputs:
%  r_path, radius along the ray path
%  theta_path, heliocentric angle along the path
%  dtau_group, group time along the ray
%  dtau_phase, phase time along the ray
%  path_length, total length of the ray (from upper to lower
%  turning point)
%
% NOTE: this code only computes the ray path from the upper turning
% point to the to lower turning point.  
%
%
% see ray_codes.ps for details
%
% aaron birch may 22 2002
 


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

%
% square of horizontal wavenumber
%
khs=l*(l+1)*r.^(-2);  

%
% square of radial wavenumber as a function of radius
%
krs=(w0^2-wcs)./cs-khs;  

%
% find where krs 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;  

%
% index of first point above the lower turning pt
%
first_index=min(I+1);

%
% index of last point before upper turning pt
%
last_index =max(I);    


%
% NOTE:  after this point i have essentially not documented this code
%   I give a description of the equations of ray tracing an
%   outline of the numerical approach in ray_codes.ps
%






%
% define for later
%
num=sqrt(khs)./r;  %% numerator of the integral
f=krs;             %% argument of the sqrt in the denom.

%
% do the lower turning point;
%

i=first_index;
j=2;

a=num(i-1)-( num(i)-num(i-1)) / (r(i)-r(i-1))*r(i-1) ;
b=( num(i)-num(i-1)) / (r(i)-r(i-1));

c=krs(i-1)-(krs(i)-krs(i-1))/ (r(i)-r(i-1))*r(i-1);
d=(krs(i)-krs(i-1))/(r(i)-r(i-1));

theta(1)=0;
r_path(1)=-c/d;

theta(j)=2*a*sqrt(c+d*r(i))/d+(2*r(i)/3/d-4*c/3/d^2)*b*sqrt(c+d*r(i));
r_path(j)=r(i);
j=j+1;


dtheta=zeros(1,last_index-first_index+3);



for (i=first_index:(last_index-1))

  a=num(i)-( num(i+1)-num(i)) / (r(i+1)-r(i))*r(i) ;
  b=( num(i+1)-num(i)) / (r(i+1)-r(i));
  c=krs(i)-(krs(i+1)-krs(i))/ (r(i+1)-r(i))*r(i);
  d=(krs(i+1)-krs(i))/(r(i+1)-r(i));

  dtheta(j)=2*a*sqrt(c+d*r(i+1))/d+(2*r(i+1)/3/d-4*c/3/d^2)*b*sqrt(c+d*r(i+1));
  dtheta(j)=dtheta(j)-2*a*sqrt(c+d*r(i))/d-(2*r(i)/3/d-4*c/3/d^2)*b*sqrt(c+d*r(i));
  theta(j)=theta(j-1)+dtheta(j);
  r_path(j)=r(i+1);
  
  j=j+1;

end;

i=last_index;
a=num(i)-( num(i+1)-num(i)) / (r(i+1)-r(i))*r(i) ;
b=( num(i+1)-num(i) ) / (r(i+1)-r(i)) ;
c=krs(i)-(krs(i+1)-krs(i))/ (r(i+1)-r(i))*r(i);
d=(krs(i+1)-krs(i))/(r(i+1)-r(i));


dtheta(j)=-2*a*sqrt(c+d*r(i))/d-(2*r(i)/3/d-4*c/3/d^2)*b*sqrt(c+d*r(i));
theta(j)=theta(j-1)+dtheta(j);
r_path(j)=-c/d;

theta_path=theta;

dtau_group=0;
dtau_phase=0;
path_length=0;
dtau_group=zeros(1,length(r_path));
dtau_phase=zeros(1,length(r_path));
dtau_group(1)=0;
dtau_phase(1)=0;

for (i=2:(length(r_path)-1))
  rm=(r_path(i)+r_path(i-1))/2;
  ds=sqrt(rm^2*(theta(i)-theta(i-1))^2+(r_path(i)-r_path(i-1))^2);
  c=interp1(r,sqrt(cs),rm );
  wac=interp1(r,wcs,rm);

  dtau_group(i)=dtau_group(i-1)+ds/c/sqrt(1-wac/w0^2);
  dtau_phase(i)=dtau_phase(i-1)+sqrt(1-wac/w0^2)*ds/c;
  path_length=path_length+ds;

  tds(i)=ds;
end;



dtau_phase=real(dtau_phase);
dtau_group=real(dtau_group);
r_path=real(r_path);
theta_path=real(theta_path);
path_length=real(path_length);