%An Exoplanet Transit 
%Data and Code by Erin Leonard

%The goal of this code is to create two side-by-side plots, one displaying
%the data with a best fit line, the other a visualization of an exoplanet
%transit. The data I will use was taken on October 24th while observing
%Wasp-12b for a duration of about 200 minutes. 

clear
clf
pause on

%First, read in the data from a .txt file
file=fopen('data.txt')
B=fscanf(file, '%g %g', [2 inf])
fclose(file)

%Then find the best fit polynomial to fit the data (least squares)
x1=[1:192];
p=polyfit(x1,B(1,:),6)
fx=polyval(p,x1)

%With the data read in for the first plot, we now turn our attention to the
%second plot which will display a planet orbiting a picture of a sun-like
%star using ODEs

%Globalize the variables
global ms G M

ms=10^6
G=3*10^(-6)
dt=10^(-3)
t=1
tmax=2200*dt
M=1

%Initialize both the velocity vectors and the position vectors
vx=zeros(tmax/dt+1,1);
vy=zeros(tmax/dt+1,1);
vx(1)=0;
vy(1)=2;
x=zeros(tmax/dt+1,1);
y=zeros(tmax/dt+1,1);
x(1)=1;
y(1)=0;

%Iterate over a loop to create an array of the position and velocity of the
%planet over time
for i=1:tmax/dt+1;
    r=sqrt(x(i)^2+y(i)^2);
    vx(i+1)=(-G*ms/r^3)*x(i)*dt+vx(i);
    vy(i+1)=(-G*ms/r^3)*y(i)*dt+vy(i);
    x(i+1)=vx(i)*dt+x(i);
    y(i+1)=vy(i)*dt+y(i);
    
end

f=0 %for movie commands

%Now that we have all of our data we can start plotting. Since we want the
%plots to follow each other, be somewhat synchronized, we want to put them
%in the same loop so that they plot points at the same time

for i=1:10:190  %number of data points 
    
    subplot(2,2,1)
    hold all
    plot(x1,B(1,:),'+b', x1, B(2,:), '+b')
    plot(x1,fx,'-r','LineWidth',3)  
    axis([0.0 200.0 0.8 1.2])
    plot(i, fx(i), '-oy', 'LineWidth',7)
    title('Data of Exoplanet Transit')
    xlabel('Time [min]')
    ylabel('Percent Intensity')
    pause(0.1)
    
    t=10*i+200  %conversion from data point number to time in the ODE
    
    subplot (2,2,2)
    hold all
    img = imread('sun.png');
    ximage=[-1.5 1.0; -1.5 1.0];
    yimage=[0 0;0 0];
    zimage=[1.0 1.0; -1.0 -1.0];
    surf(ximage, yimage, zimage, 'CData',img, 'Facecolor', 'texturemap');
    plot(x(t),y(t), '-og', 'LineWidth', 7)
    title('Visual of Star and Exoplanet')
    axis([-1.5 1.0 -1.5 1.0])
    axis off
    view([0,1,0])
    pause(0.1)
    
    f=f+1 %increases the label of the .png by one for every loop
    saveas(gcf, ['transit_mov', sprintf('%03d',f),'.png'])
end


