clear u w maxz kk vzero raytime

for ii=[1:1:numlayer]
	maxz(ii)=zees(ii+1); %layer boundary
	kk(ii)=(vees(ii+1)-vees(ii))/(zees(ii+1)-zees(ii));  % gradient
	vzero(ii)=(vees(ii)*zees(ii+1)-vees(ii+1)*zees(ii))/(zees(ii+1)-zees(ii));  
% intercept
end

maxz(numlayer)=100000; % arbitrarily deep value for half space

for pindex=[1:numberofrays]  % loop over all rays
	raytime(pindex)=0; % reset  ray travel time
	rayxbegin=5;      % SPECIFY: ray start x coordinate (u0)
	z0=0;	    % starting depth
	deltaangle= (endangle-startangle)/numberofrays;  % angle incriment
	angle=startangle + (pindex-1)*deltaangle;  % current ray angle
	p=sin(pi*angle/180)/vzero(1); % ray parameter
	for layer=[1:numlayer]   % loop over layers
		if layer==1 
			x0=rayxbegin; % set initial takeoff distance
		else
			x0=ulast;  % set ray start x at top of next layer
		end
		limz=((1/p)-vzero(layer))/kk(layer) ;    % limiting depth for ray in layer
		if kk(layer)<= 0   % if gradient is negative, it will go to next layer
			limz=maxz(layer); % set maxdepth to next layer boundary
		end
		if limz<0 , break, end     % If ray never starts, limz<0,  end ray loop
		theta0=asin(p*vees(layer));  % initial angle (use both + and - if z0 > 0)
%

		if layer < numlayer        % if more layers below
			if limz >= maxz(layer)              % check if penetrates next layer
				limz=maxz(layer);          % if so,  max depth in ith layer is maxz for that layer
			end	
		end			
													
							
		zcenter= -vzero(layer)/kk(layer) ;       % ray circle center depth
		xcenter=x0+cos(theta0)/(p*kk(layer)) ;  % ray circle center x value
%
		index=0;
		for dindex=[(layer-1)*pointsper+1: pointsper*layer]      % loop over valid depths
			index=index+1; % incriment depth index
			depth=(index-1)*(limz-zees(layer))/(pointsper-1)+zees(layer);% incriment depths over ray segment
			w(dindex,pindex)=depth;   % get depth of point
			% get x coordinate for that depth
			u(dindex,pindex)=abs(x0+(cos(theta0)-cos(asin(p*kk(layer)*...
				(-zcenter+depth))))/(p*kk(layer)));
			ulast=u(dindex,pindex);  % save x value of last point
		end

		raytime(pindex)=raytime(pindex)+abs(2*(1/kk(layer))*...
			log(tan(0.5*asin(p*(kk(layer)*depth+vzero(layer))))/...
			tan(0.5*theta0)));  % calculate travel time 
		if limz< maxz(layer), break, end; % ray doesn't hit next layer, jump out
	end  % loop to next layer

% reflect around max depth

	for ii=[1:pointsper*layer-1]  % =1 to 101 for two layers
		dindex=pointsper*layer+ii;  % index of new point
		idx=pointsper*layer-ii; % index of point to be reflected
		w(dindex,pindex)=w(idx,pindex);
		u(dindex,pindex)=2*xcenter-u(idx,pindex);
	end

	dist(pindex)=u(dindex,pindex);  % save ray distance
	rayp(pindex)=p;				% save ray parameter
end %ray loop