function [END_PTS, XXX] = impact(impact_loc, n)
% function impact.m -- calculates and plots ejecta motion
% [END_PTS, XXX] = impact(impact_loc, n)

% Input consists of projectile diameter [m] (proj_d), location of impact (impact_loc) [d]
% which is of form [lat long], and number of particles (n) whose motion is to be followed.

% Serina, last modified 21 June 2004
% calc_param.m (script), checkpoint.m, coor_conv.m, int_traj.m, ellipsd.m

% imports values found in parameters
% calls int_traj.m

close all
global density_1 density_sur ell_x ell_y ell_z lay_depth omega
global XXX END_PTS
calc_param

lat_d = impact_loc(1);
long_d = impact_loc(2);
lat = impact_loc(1)*2*pi/360;
long = impact_loc(2)*2*pi/360;

END_PTS = ones(n,4);

impact_z = ell_z * sin(lat);
ellipse_y = sqrt((1 - impact_z^2/ell_z^2) * ell_y^2);
ellipse_x = sqrt((1 - impact_z^2/ell_z^2) * ell_x^2);
impact_x = sqrt((ellipse_x*ellipse_y)^2 / (ellipse_y^2 + ellipse_x^2* tan(long)^2));
impact_y = impact_x * tan(long);

if impact_z == ell_z,
   impact_x = 0;
   impact_y = 0;
end

% assign initial velocity (magnitude and elevation angle)
% based on discussions with Onose-san, 21 June 2004
for j = 1:n,
   N(j,1) = 9*pi/36; %angle from surface tangent place, in asteroid frame
	N(j,2) = 2*pi *(j - 1 + 60) / 100; %angle in plane, in asteroid frame
   N(j,3) = 9.4; %velocity magnitude - r
end

M = coor_conv(N,2);

Rot_lat = Ry(M,(pi/2 - lat));
Rot_long = Rz(Rot_lat,long);

smax = [ell_x ell_y ell_z];
ellipsd(smax);
view(long_d + 90,lat_d + 25)
axis equal
hold
plot3(impact_x,impact_y,impact_z,'b*')
for j = 1:n,
   clear X end_pt iter
   in_pt = [impact_x impact_y impact_z];
   in_v = Rot_long(j,:);
   [end_pt,X,iter] = int_traj(in_pt, in_v);
   END_PTS(j,1:3) = end_pt;
   END_PTS(j,4) = iter;
   [k1 k2] = size(X);
   for m = 1:k1,
      if rem(m,50) == 0;
         n = m/50;
         XXX(j,n,:) = X(m,:);
      end
   end
end
hold

%figure(2)
%ellipsd(smax);
%axis equal
%view(long_d, lat_d)
%hold
%plot3(impact_x,impact_y,impact_z,'b*')
%for j = 1:n,
%   magic_str = ['plot3(X',int2str(j),'(:,1),X',int2str(j),'(:,3),X',int2str(j),'(:,5))'];
%   eval(magic_str);
%end
%hold