%8/22/03
clear all;

c = 3e8;
L = 4000;
omega_laser = 2*pi*c/1e-6;
theta = 0.8;
phi = 0.2;
h_plus = 1e-20;
h_cross = 1e-20;
tau = L/c;

f_grav = 100;
data = load('test.dat');
data(1,2)=data(3,2);
data(2,2)=data(3,2);

t_final = 10/f_grav;
t_step = 1/f_grav/1000;

t0 = [0:t_step:t_final];

k_grav = 2*pi*f_grav/c;
Omega_grav = 2*pi*f_grav;


kx = sin(theta)*cos(phi);
ky = sin(theta)*sin(phi);

hxx = -cos(theta)*sin(2*phi)*h_cross+(cos(theta)^2*cos(phi)^2-sin(phi)^2)*h_plus;
hyy = cos(theta)*sin(2*phi)*h_cross+(cos(theta)^2*sin(phi)^2-cos(phi)^2)*h_plus;
Phi_x_exact = hxx*omega_laser./(2*c*k_grav.^2*(1-kx^2)).*(k_grav.*(1+kx)*2.*cos(2*pi*f_grav*t0+k_grav*(1-kx)*L*0.5).*sin(0.5*L*k_grav*(1-kx))+k_grav*(1-kx)*2.*cos(2*pi*f_grav*t0+2*L*k_grav-k_grav*(1+kx)*L*0.5).*sin(0.5*L*k_grav*(1+kx)));
Phi_y_exact = hyy*omega_laser./(2*c*k_grav.^2*(1-ky^2)).*(k_grav.*(1+ky)*2.*cos(2*pi*f_grav*t0+k_grav*(1-ky)*L*0.5).*sin(0.5*L*k_grav*(1-ky))+k_grav*(1-ky)*2.*cos(2*pi*f_grav*t0+2*L*k_grav-k_grav*(1+ky)*L*0.5).*sin(0.5*L*k_grav*(1+ky)));


h_xx = hxx*cos(2*pi*f_grav*t0);
h_yy = hyy*cos(2*pi*f_grav*t0);

h_xx_old=h_xx(1);
h_xx_old_old = h_xx(1);
h_yy_old=h_yy(1);
h_yy_old_old = h_yy(1);

for i =1:length(t0)
    h_xx_prime = (h_xx(i)-h_xx_old)/t_step;
    h_xx_prime_prime = (h_xx(i)-2*h_xx_old+h_xx_old_old)/t_step^2;
    Phi_x_approx3(i) = omega_laser*L/(4*c)*(4*h_xx(i)+3*tau*h_xx_prime+2.5*tau^2*h_xx_prime_prime);
%    Phi_x_approx3(i) = omega_laser/(2*c)*(L*h_xx(i)+L^2/(2*c)*h_xx_prime+L^3/(6*c^2)*h_xx_prime_prime)+omega_laser/(2*c)*(L*h_xx(i)+L^2/(2*c)*h_xx_prime+L^3/(6*c^2)*h_xx_prime_prime);
    h_xx_old_old = h_xx_old;
    h_xx_old = h_xx(i);
    
        h_yy_prime = (h_yy(i)-h_yy_old)/t_step;
    h_yy_prime_prime = (h_yy(i)-2*h_yy_old+h_yy_old_old)/t_step^2;
    Phi_y_approx3(i) = omega_laser*L/(4*c)*(4*h_yy(i)+3*tau*h_yy_prime+2.5*tau^2*h_yy_prime_prime);
%    Phi_x_approx3(i) = omega_laser/(2*c)*(L*h_xx(i)+L^2/(2*c)*h_xx_prime+L^3/(6*c^2)*h_xx_prime_prime)+omega_laser/(2*c)*(L*h_xx(i)+L^2/(2*c)*h_xx_prime+L^3/(6*c^2)*h_xx_prime_prime);
    h_yy_old_old = h_yy_old;
    h_yy_old = h_yy(i);

end

t0 = 1000*t0;
figure

amp_exact = max(Phi_y_exact)
amp_approx = max(data(:,2))
ky
ratio = amp_exact/amp_approx

plot(t0,Phi_y_exact,1000*data(:,1), data(:,2));
legend('Exact phase shift', 'e2e-calculated phase shift');
ylabel('\Delta\Phi_{laser}, round trip (rad)');
xlabel('Time (ms)');
title(['Approximating Laser Phase Shift: f_{GW} = ' num2str(f_grav) 'Hz']);