% 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 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');