function varargout=yprime(t,y,flag)
%This is the integration function used by ode45
%the inputs are the time and the solution vector, y.
%it returns the time derivative of y, yprime.
%It also specifies so called events if called
%with the appropriate flag. Those specify when
%to stop the  integration, apart from the
%maximum time which is specified when calling ode45.
%Written by Keli Karason, karason@mit.edu
switch flag
  case ''
    varargout{1}=f(t,y);
  case 'events'
    [varargout{1:3}]=events(t,y);
end

% -----------------------------------------------------------
function yprime=f(t,y)
x=y(1);
z=y(2);
px=y(3);
pz=y(4);

[cc,dcdz]=veldep(z);
    
xprime=cc^2*px;
zprime=cc^2*pz;
pxprime=0;
pzprime=-1/cc*dcdz;

yprime=[xprime; zprime; pxprime; pzprime];
    

% -----------------------------------------------------------
function [value,isterminal,direction]=events(t,y);

x=y(1);
z=y(2);

value=[z;z-50;x;x-100];
isterminal=[1;1;1;1];
direction=[-1;0;-1;0];