% Major function 
% free throw on asteroids 3D model, A rotating body
% input the size, density of small asteroids, and the number of particles
% released, R, density, and n

% the trajectories of those particles that redeposited onto the surface
% would be recorded

% the five sections in the end are used to plot different types of figures,


clear

global R density       %define the disk surface

R=245*30;       %asteroid radius, unit: m
density=1.26;                    %asteroid density, unit: g/cc

omega=4.06*10^(-4);             %angular velocity, unit s^-1
G=6.674*10^(-11);               %gravity constant                   
M=4/3*pi()*R^3*density*1000;                 %mass of asteroid, unit: kg
V_max=sqrt(8/3*G*pi()*density*1000)*R;             % escape velocity

% ejection Latitude from Cheslay et al. 2020, from -1 to 1, stepsize=0.2
sin_belta=[62/623 18/293 44/867 5/82 19/182 175/817 47/225 39/635 5/56 26/529];

n=50000;                            % the number of particles
V_terminal=rand(n,1).*(5-0.05)+0.05;   %intial velocity of the grains, unit m/s
%V_terminal=rand(n,1).*(V_max-0.05)+0.05;  % initial velocity of the grains, max=escape

sigma=rand(n,1)*pi();
V_z0=V_terminal.*sin(sigma);             % velocity in the z0 direction (not z direction)
V_xy0=abs(V_terminal.*cos(sigma));       % velocity in the xy0 plane
alpha=rand(n,1)*2*pi();                  % 
%belta=rand(n,1).*pi()/2;                 % initial position angle/Latitude
%belta=rand(n,1).*pi()-pi()/2;                 % initial position angle
%belta=zeros(n,1)+pi()/4;
belta=double.empty;

for k=1:10
    range=asin(-1+k*0.2)-asin(-1+(k-1)*0.2);
    belta_latitude=rand(round(sin_belta(k)*n),1)*range+asin(-1+(k-1)*0.2);
    belta=[belta; belta_latitude];
end

if length(belta)<n                            %belta may be smaller than n because of rounding operation
    for blank_num=1:(n-length(belta))
        belta=[belta;0];
    end
end

if length(belta)>n                            %belta may be greater than n because of rounding operation
        belta=belta(((length(belta)-n)+1):end);
end

y0=zeros(n,1)+R*cos(belta); z0=zeros(n,1)+R*sin(belta); x0=zeros(n,1);      % initial position, unit m
V_x=V_xy0.*cos(alpha);
V_rot=R*cos(belta)*omega;                                                   % rotational velocity
V_x=V_x-V_rot;                                                              % compound velocity in the x direction
V_y=V_z0.*cos(belta)+V_xy0.*sin(alpha).*cos(pi()/2-belta);
V_z=V_z0.*sin(belta)+V_xy0.*sin(alpha).*cos(pi()-belta);

tspan=1.0E+08;                           % maximum running time IMPORTANT!
l=zeros(n,1);                            % the length of trajectory
DL=nan(n,1);                             % ratio of trajectory length to asteroid radius
Te=nan(n,1);                             % flying time
displa_angle=nan(n,1);                   % the angle between the initial point and the landing site
landing_latitude=nan(n,1);
% ODE event
% when the particle fall back onto the surface, stop integration
% result export
% result=zeros(1,12);

for i=1:n
    opts=odeset('RelTol',1e-9,'AbsTol',1e-12,'events',@fallback_aster);          % adjust the relative and absolute error tolerances, and introduce fallback event, Refine 20
    zyx0=[z0(i),V_z(i),y0(i),V_y(i),x0(i),V_x(i)]; %initial velocity and position  
    [t,zyx,te,ye,ie]=ode45(@gra_aster,[0,tspan],zyx0,opts);           %zyx:z, vz, y, vy, x, vx; te: time when stop integration

    if isempty(te)==0                              % if te is not empty
        Te(i)=te;       
        for j=2:length(zyx(:,1))                                                        % calculate the length of trajectory
            deltaL=sqrt((zyx(j,1)-zyx(j-1,1))^2+(zyx(j,3)-zyx(j-1,3))^2+(zyx(j,5)-zyx(j-1,5))^2);
            l(i)=l(i)+deltaL;
        end
        l(i)=l(i)/1000;                                     %moving distance, unit: km
        % calculate the angle between the initial point and the landing site
        % noticing that the launching site also moved
        y0(i)=R*cos(belta(i))*cos(omega*te);
        x0(i)=R*cos(belta(i))*sin(omega*te);
        displa_angle(i)=acos((zyx(end,1)*z0(i)+zyx(end,3)*y0(i)+zyx(end,5)*x0(i))/R^2);
        landing_latitude(i)=atan(zyx(end,1)/sqrt(zyx(end,3)^2+zyx(end,5)^2));
    else 
        l(i)=nan;                                           %escape the gravity
    end

    DL(i)=l(i)/R*1000;                                  %ratio of trajectory length to asteroid radius
    zyx=zyx/1000;                                       %zyx, unit: km 
%    r_final(i)=radial_distance(end);                             %record the final r position

   %   plot1=plot(radial_distance,zyx(:,1),'--','linewidth',1); hold on % r vs z
   %   plot1.Color(4) = 0.1; hold on;
   
      % plot3(zyx(:,5),zyx(:,3),zyx(:,1),'--','linewidth',3); hold on   % x,y,z
      %loglog(zyx(:,3),zyx(:,1),'--'); hold on
   %  xlabel('X'); ylabel('Y'); zlabel('Z'); hold off
   % plot(radial_distance,zyx(:,1),['-',markerstyle(i)],'linewidth',2,'MarkerSize',6); hold on % r vs z
   % plot(zyx(:,3),zyx(:,5),['-',markerstyle(i)],'linewidth',2,'MarkerSize',6); hold on %plot x vs y
    
    %result=[result;V_terminal(i) z0(i)/AU apha(i) r_final(i) zyx(end,1) zyx(end,3) zyx(end,5) t(end) t_eff_co(i) t_eff_fo(i) t_eff_fa(i) t_eff_bl(i)];
end
% the number of NaN, particle escape the gravity
Num_escape=sum(isnan(Te));

%{
r=R/1000;                             % unit km
axiszoom=2;                           % zoom in or out the axis
[X Y Z]=sphere(50);         % ellipse shape
%C=zeros(length(X),length(X));
CO(:,:,1) = ones(length(X)).*0; % blue
CO(:,:,2) = ones(length(X)).*1; % blue
CO(:,:,3) = ones(length(X)).*1; % blue
%surf(r*X,r*Y,r*Z,CO,'FaceAlpha',.9,'EdgeAlpha',.2,'EdgeColor','none'); hold on;    %Alpha controls the transparency. 
surf(r*X,r*Y,r*Z,'EdgeColor','none');
%set(gca, 'YDir','reverse');       % reverse the X direction
set(gca,'linewidth',1);
set(gca,'FontSize',13);
%set(gca,'XAxisLocation','origin');
%set(gca,'YAxisLocation','origin');
%set(gca,'CameraPosition',[0, 1, 0]);
%set(gca,'visible','off');
hAxis = gca;
hAxis.XRuler.FirstCrossoverValue  = 0; % X crossover with Y axis
hAxis.YRuler.FirstCrossoverValue  = 0; % Y crossover with X axis
hAxis.ZRuler.FirstCrossoverValue  = 0;
hAxis.ZRuler.SecondCrossoverValue = 0; % Z crossover with Y axis
hAxis.XRuler.SecondCrossoverValue = 0; % X crossover with Z axis
hAxis.YRuler.SecondCrossoverValue = 0; % Y crossover with Z axis
xlim([-r.*axiszoom r.*axiszoom]);
ylim([-r.*axiszoom r.*axiszoom]);
zlim([-r.*axiszoom r.*axiszoom]);
xlabel('X','FontSize',20); ylabel('Y','FontSize',20); zlabel('Z','FontSize',20); hold off
%}

%%{
X=displa_angle;X = X(~any(isnan(X),2));X = X(~all(isnan(X),2));X = X(~isnan(X));
histogram(X,50,'Normalization','probability'); 
xlabel('Displacement angle (radian)','FontSize',20);
ylabel('Normalized frequency','FontSize',20); 
ax = gca;
ax.LineWidth=1.5;
ax.FontSize=20;
xlim([0 pi()]);
%ylim([0 0.05]);           % modify the y axis
%yticks(0:0.01:0.05)
%hold off
%%}
large_transport=sum(displa_angle>pi()/4);

%scatter3(sigma,belta,V_terminal,[],displa_angle);
%colorbar;
%scatter(belta,V_terminal,[],displa_angle);
%colorbar('FontSize',25,'Ticks',[0,pi()/4,pi()/2,pi()/4*3,pi()],'TickLabels',{'0','1/4\bf\pi','1/2\bf\pi','3/4\bf\pi','\bf\pi'});
%xlabel('Latitute','FontSize',20);
%ylabel('Initial velocity','FontSize',20);

%{
scatter(belta,landing_latitude,[],displa_angle);
cb=colorbar;
cb.Title.String ='Displacement';
cb.FontSize=20; cb.Ticks=[0,pi()/4,pi()/2,pi()/4*3,pi()];cb.TickLabels={'0','1/4\bf\pi','1/2\bf\pi','3/4\bf\pi','\bf\pi'};
xlabel('Launching Latitude','FontSize',20);
ylabel('Landing Latitude','FontSize',20);
%}

%{
X=landing_latitude;X = X(~any(isnan(X),2));X = X(~all(isnan(X),2));X = X(~isnan(X));
histogram(sin(X),10,'Normalization','probability');
xlabel('Sine Landing Latitude','FontSize',20);
ylabel('Normalized frequncy','FontSize',20); 
ax = gca;
ax.LineWidth=1.5;
ax.FontSize=20;
xlim([-1 1]);
ylim([0 0.17]);
%}

% writematrix(X)
% open('Bennu_P50000-N2171-Impact hist.fig');
% open('Bennu_P50000-N2171_disp.fig');
% hl = get(gca,'Children');


