%Planetary Distruction by Stellar Tides
%Michael Josling Hill
%michaelhill@berkeley.edu



clear;
clf;


global mStar mPlanet rPlanet aPlanet G rDestroy
mStar    = 1.98e33; %Star mass (g)
mPlanet = 5.97e27; %Planet mass (g)
rPlanet = 6.3675e9; %Planet radius (cm)
aPlanet = 7.4e11; %Semi-major axis (cm) ~.05 AU
G       = 6.67428e-8; %cm^3/(g*s^2)

% initial positions (cm)
rx0 = aPlanet;
ry0 = aPlanet/6.5;
% initial velocities (cm/s)
vx0 = 6.0e6; 
vy0 = -1.1e7;

Y = [rx0;ry0;vx0;vy0]; % by convention, should be a column vector






%Planet's binding gravitational force vs differential stellar gravitation
rDestroy = rPlanet*power(mStar/mPlanet, 1/3); %cm 

fprintf(['Tidal Destruction Radius = ', num2str(rDestroy,'%10.2e\n'), ' cm \n']);



%Draw red tidal destruction zone
plot(0,0,'X');
fluff = (1/6)*rx0;
axis ([-(rx0+2*fluff) rx0+2*fluff -(rx0+fluff) rx0+fluff])
N=50;
t = 2*pi/N*(1:N);
x = rDestroy*sin(t);
y = rDestroy*cos(t);
hold on;
fill(x, y,'r');

%Draw star
t = 2*pi/N*(1:N);
xx = (1/10)*rDestroy*sin(t);
yy = (1/10)*rDestroy*cos(t);
hold on;
fill(xx, yy,'y');




t    =  0.0;
tmax = 4.32e5; %seconds 
dt   =  900; %seconds 
n    = tmax/dt; % number of steps

tt = zeros(floor(n+1),1);
Yt = zeros(floor(n+1),floor(size(Y,1)));
tt(1)   = t;
Yt(1,:) = Y';


title('Perturbed Orbit Wanders Into Destuction Zone','fontsize',17);
xlabel('Distance (cm)')
ylabel('Distance (cm)')

pcb = 0; %planet plot callback
for i=1:n

    f1 = kepler_ode(t     , Y                );
    f2 = kepler_ode(t+dt/2, Y+f1 .*  (dt/2)  );
    f3 = kepler_ode(t+dt/2, Y+f2 .*  (dt/2)  );
    f4 = kepler_ode(t+dt  , Y+f3 .*   dt     );
   
    Y = Y + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
    
    t = t+dt;
    tt(i+1)   = t;  % this records the time at each step
    Yt(i+1,:) = Y'; % this records the vector Y at each step --> Yt(:,1) has all x coordinates
    
    set(pcb,'Visible','off');
    hold on; %hold previous draws
    pcb = plot(Yt(i,1),Yt(i,2),'-ob','MarkerSize',20); %plot new point
    hold on;
    plot(Yt(i,1),Yt(i,2),'.b', 'linewidth', 2); %draw line
    drawnow; %draw it immediately
    
    if( abs(sqrt(Yt(i,1)^2 + Yt(i,2)^2)) <= rDestroy)
        break; %The planet wandered into the tidal destruction zone;
    end
    
end


