function [end_pt,X,iter]= int_traj(in_pt, in_v);
% int_traj.m -- integrates ballistic trajectory of particle, taking into account gravity 
% and rotation.

% Input values are the initial position (in_pt) and velocity vectors (in_v)(in cartesian 
% coordinates). The function integrates in 5s chunks (consisting of ~50 timesteps), then 
% checks to see if the particle has reimpacted the surface. If so, then integration is 
% stopped and the landing spot is one output (end_pt). If the particle does not reimpact 
% the surface within %%chunks (i.e., %%s), then end_pt = [0,0,0]. An array is also outputted
% (X), with all position and velocity information for each timestep, allowing a plot to be 
% drawn.

% Serina, last modified 21 June 2004
% traj_prime.m, checkpoint.m

x0 = in_pt(1);
y0 = in_pt(2);
z0 = in_pt(3);
vx0 = in_v(1);
vy0 = in_v(2);
vz0 = in_v(3);
cond0 = [x0 vx0 y0 vy0 z0 vz0];

% imports values found in parameters.m
global density_1 density_sur ell_x ell_y ell_z lay_depth omega

% length of "time" before model checks to see if particle has impacted with the surface
t0 = 0;
tf = 50;
% total number of times model will check before assuming particle is "lost"
T = 400;

i = 1;

global X
[t,x] = ode45('traj_prime',[t0 tf],cond0);
k = size(x);
k = k(1);
cond0 = x(k,:);
K = k;
X = x;

[pt_out, test_var] = checkpoint([cond0(1) cond0(3) cond0(5)], ell_x, ell_y, ell_z, 2);

% specify here how long before the model assumes the particle is lost
while i < T,
   % tests to see if particle has reimpacted the surface
   if test_var ~= 1,
      i = i+1;
      [t,x] = ode45('traj_prime',[t0 tf],cond0);
      k = size(x);
      k = k(1);
      cond0 = x(k,:);
      [pt_out, test_var] = checkpoint([cond0(1) cond0(3) cond0(5)], ell_x, ell_y, ell_z, 2);
      for j = 1:k,
         X(K + j,:) = x(j,:);
      end
      K = K + k;
      % tests here to see if particle is on hyperbolic (escape) orbit
      if (pt_out(1)^2 + pt_out(2)^2 + pt_out(3)^2) > 70*ell_x,
         G = 6.6726E-11; 
         M = 2500*4*pi/3*ell_x*ell_y*ell_z;
         E = 0.5*(x(k,2)^2 + x(k,4)^2 + x(k,6)^2) - G*M/(sqrt(x(k,1)^2 + x(k,3)^2 + x(k,5)^2));
         if E > 0,
            end_pt = [-1 0 0];
            iter = i;
            i = T; 
         end
      end
   % if particle has reimpacted ...
   else
      iter = i;
      i = T;
   end
end

if pt_out == [0 0 0],
   end_pt = [0 0 0];
   iter = T;
end

k = size(x);
k = k(1);
end_pt = [x(k,1) x(k,3) x(k,5)];
plot3(X(:,1),X(:,3),X(:,5),'r')