% Gravity simulation which begins with a single Jupiter-like planet
% orbiting a sun-like star, plus concentric circles of 'masslets' that act
% like an accretion disc or dust cloud around the star. The planet and
% masslets don't interact, and the star mass is assumed to be much larger
% than then mass of the planet, even after it has shed mass.
%
% After N planet rotations, the star loses fraction f of it's mass. The simulation
% uses Verlet integration to determine the subsequent dynamics for another
% M planet periods, before resetting.
function star2whitedwarf

%% INPUTS %%

N = 1;          %Number of orbital periods before mass shed
M = 3;          %Number of (original) orbital periods after mass shed
f = 0.3;        %Fraction of star mass shed
a = 5.202;      %Semi-major axis of planet orbit in AU
MS = 2.0;       %Mass of star in solar masses
dt = 0.05;      %Timestep in years
Rs = 0.8;       %Radius of star in AU
Rp = 0.1;       %Radius of planet in AU
nrings = 60;    %Number of rings of masslets
rmin = 2;       %Radius of minimum ring /AU
dr = 0.1;       %Difference in ring radii /AU
m0 = 50;        %Masslets in inner ring - this sets spacing for the outer rings
fsize = 18;     %Fontsize of graph grid
msize = 2;      %Marker size of masslets
xmax = 3*a;     %Max scale of plot

%

%% INITIAL CONDITIONS %%

%Time in years
t = 0;

%Planet
vy0 = 2*pi*sqrt( MS/a ); %Initial y velocity in AU per year
vx0 = 0; %Initial x velocity in AU per year
x0 = a; %Initial x position (from star) in AU
y0 = 0; %Initial y position (from star) in AU
xx = x0; yy = y0;

%Determine original planet period in years, and subsequent planet period
%after mass shed
T = sqrt( (a^3)/MS );

%Rings
for n=1:nrings
    %Ring radius
    r = rmin +(n-1)*dr;
    
    %Calculate masslets per ring
    m = round( m0*r/rmin );
    
    %Determine initial x,y coordinates
    theta = 0 : 2*pi/m : 2*pi*(m-1)/m;
    x0 = [x0,r*cos(theta)];
    y0 = [y0,r*sin(theta)];
    
    %Determine initial x,y velocities in AU/year
    vx0 = [vx0,-2*pi*sqrt( MS/r )*sin(theta)];
    vy0 = [vy0,2*pi*sqrt( MS/r )*cos(theta)];
end
x = x0; y = y0; vx = vx0; vy=vy0;

%

%% INITIALIZE ANIMATION %%

fig = figure('name','Star to white dwarf','color',[1 1 1],'renderer','painters',...
    'units','normalized',...
    'position',[0.05 0.05 0.85 0.85]);
theta = linspace(0,2*pi,200);
plot( 0,0,'r*','markersize',50 ); hold on;
star = patch( Rs*cos(theta), Rs*sin(theta),[1,0,0],'edgecolor',[1,0,0] );
masses = plot( x, y, 'go','markersize',msize );
planet_trail = plot( x(1), y(1),'b--' );
planet = patch( Rp*cos(theta) + x(1), Rp*sin(theta) + y(1),[0,0,1] );
axis equal
xlim([-xmax,xmax]); ylim([-xmax,xmax]); axis manual
set(gca,'fontsize',fsize);
grid on;
box on;
stop = 0;
tit = title(['Star mass M = ',num2str(MS),', t=0'],'fontsize',fsize);

%

%% CONSTRUCT SIMULATION %%

stop = 0;
iterations = ceil( ( N + M )*T/dt );
while stop ==0
    try
        n=1; MMS = MS;
        for n=1:iterations
            
            %Shed mass fraction f of star after N planet periods
            if t(n) > N*T
                MMS = MS*(1-f);
                RRs = Rs*( (1-f)^(1/3) );
                set( star, 'xdata',RRs*cos(theta),'ydata',RRs*sin(theta),'facecolor',[1,0,0]  );
            end
            
            %Update time /years
            t(n+1) = t(n) + dt;
            
            %Update distance of planet and masslets from star at (0,0)
            r = sqrt( x.^2 + y.^2 );
            
            %Update positions and velocities using Verlet method
            ax = -(4*pi^2)*MMS*x./r.^3 ;
            ay = -(4*pi^2)*MMS*y./r.^3 ;
            dx = vx*dt +0.5*ax*(dt^2);
            dy = vy*dt + 0.5*ay*(dt^2);
            x = x + dx; y = y + dy;
            r = sqrt( x.^2 + y.^2 );
            axx = -(4*pi^2)*MMS*x./r.^3 ;
            ayy = -(4*pi^2)*MMS*y./r.^3 ;
            vx = vx + 0.5*( axx+ax )*dt;
            vy = vy + 0.5*( ayy+ay )*dt;
            
            %Update planet trajectory
            xx = [xx,x(1)]; yy = [yy,y(1)];
            
            %Update iteration counter
            n = n+1;
            
            %Update plots
            set( masses, 'xdata',x,'ydata',y );
            set( planet, 'xdata',Rp*cos(theta) + x(1),'ydata',Rp*sin(theta) + y(1) );
            set( planet_trail, 'xdata',xx,'ydata',yy );
            set(tit,'string',['Star mass M = ',num2str(MMS,2),', t=',num2str(t(n),2)],...
                'fontsize',fsize);
            drawnow;
        end
        
        %Reset parameters
        x = x0; y = y0; vx = vx0; vy=vy0;
        MMS = MS;
        xx = x0(1); yy = y0(1);
        set( star, 'xdata',Rs*cos(theta),'ydata',Rs*sin(theta),'facecolor',[1,0,0]  );
    catch
        return;
    end
end

%End of code



