function gravitational_wave
close all; clear all; clc;

%Hanford reconstructed data
H_r=load('H_r.txt');
t_H=H_r(:,1);
A_H=H_r(:,2).*1.e-21;

t=5.2:0.0001:5.3648; %s

%Konstanter:
G=6.674e-11; %N m2/kg2
c=299792458; %m/s
r=410*3.08567758e22; %luminosity distance m, 410 Mpc
M_Sun=1.989e30;  %Sun mass in kg

%Masse af de to sorte huller, samt resulterende kvidremasse
M1=36*M_Sun; %kg
M2=29*M_Sun; %kg
M=((M1*M2)^(3/5))/(M1+M2)^(1/5); %chirp mass kg
k=(96/5)*(2*pi)^(8/3)*((G*M)/(c^3))^(5/3);
P0=0.2; %initial period

phi=[ ];  
phizero=pi;
P_orb=[];
for tt=t;
    P_orb1=abs((P0^(8/3)-8/3*k*tt)).^(3/8);
    P_GW=0.5*P_orb1;
    
    phi1=integral(@(tt) 2*pi.*1./(abs((P0^(8/3)-8/3*k*tt)).^(3/8)),0,tt)*1.65;
    
    phi=[phi phi1];
    P_orb=[P_orb P_orb1];
end
P_orb;
P_GW=0.5.*P_orb;
A=(2*(G*M).^(5/3))./(c^4*r).*(pi./P_GW).^(2/3);

h=A.*cos(phi+phizero);

figure()
hold on
plot(t_H,A_H,'Linewidth',1)
plot(t-4.9425,h,'Linewidth',1.5)
title('Signal ved Hanford-detektoren (H1)','FontSize',20, 'Fontweight','normal')
xlabel('Tid [s]','FontSize',20,'FontWeight','bold')
ylabel('Strain','FontSize',20,'FontWeight','bold')
legend({'H1 Data','Aproksimerede Signal'},'Fontsize',16)

end