% Sampling Frequency = 2 Hz: LISA Band : (10^-4) - 1 Hz. 
delta = 1/2; 

% Total Number of Initial Data 

N1 = (1/delta) * 3600 * 24*7;

%Wave co-ordinates
phi = 45; %azimuthal angle (0, 180)
theta = 45; %polar angle (0, 360)

%Wave amplitude
H = (10^(-21)); 
%Wave frequency (Hertz)
w = 10^(-2); 
%Polarisation Stuff
tau = 45; 
phase = 45; 

phi = phi * pi/180; 
theta = theta * pi/180; 
tau = tau * pi/180; 
phase = phase * pi/180; 
w = 2*pi*w; 

Y21 = zeros(N1, 1); 
Y12 = zeros(N1, 1); 
Y23 = zeros(N1, 1); 
Y32 = zeros(N1,1 ); 
Y31 = zeros(N1, 1); 
Y13 = zeros(N1,1 ); 

X = zeros(N1,1 ); 
X1 = zeros(N1, 1); 

%Space Craft Positions: (in light seconds):

s1 = 10*[0 1 0]'; 
s2 = 10*[-sin(2*pi/3) cos(2*pi/3) 0]';
s3 = 10*[sin(2*pi/3) cos(2*pi/3) 0]';

l3 = sqrt(dot((s1-s2), (s1-s2))); 
l1 = sqrt(dot((s2-s3), (s2-s3))); 
l2 = sqrt(dot((s3-s1), (s3-s1))); 
l = sqrt(dot(s1, s1)); 

n3 = (s1 - s2)/l3; 
n1 = (s2 - s3)/l1;
n2 = (s3 - s1)/l2; 


% Calculation of the effect of a GW wave on the LISA Spacecraft. 
% See Armstrong, Tinto, Estabrook et al. 

k = [(sin(phi)*cos(theta)) (sin(phi)*sin(theta)) cos(phi)]';

mu1 = dot(k, s1); 
mu2 = dot(k, s2); 
mu3 = dot(k, s3); 


xi = [1 0 0]; 
xj = [0 1 0];
xk = [0 0 1];

psi1p = (1/2) * (1/(1-((dot(k, n1)^2)))); 
psi1c = psi1p; 
psi1p = H * psi1p * ((dot(n1, xi)^2) - (dot(n1, xj)^2)) * sin(tau); 
psi1c = H * psi1c * (dot(n1, xi))*dot(n1, xj)*cos(tau); 

psi2p = (1/2) * (1/(1-((dot(k, n2)^2)))); 
psi2c = psi2p; 
psi2p = H * psi2p * ((dot(n2, xi)^2) - (dot(n2, xj)^2))*sin(tau); 
psi2c = H * psi2c * (dot(n2, xi))*dot(n2, xj)*cos(tau); 

psi3p = (1/2) * (1/(1-((dot(k, n3)^2)))); 
psi3c = psi3p; 
psi3p = H * psi3p * ((dot(n3, xi)^2) - (dot(n3, xj)^2))*sin(tau); 
psi3c = H * psi3c * (dot(n3, xi))*dot(n3, xj)*cos(tau); 




% Round offs due to rate of data production

%l3 = delta * floor(l3*(1/delta)); 
%l2 = delta * floor(l2*(1/delta)); 
%l1 = delta * floor(l1*(1/delta)); 

for I = 1:N1
    t = I * delta; 
    
    %calculation of Y21

    if (t - l2 < delta)
        Y21(I)= 0; 
    else
        t1 = t - (mu3*l) - l2; 
        Y21(I) = (psi2p*sin((w*t1) + phase)) + (psi2c*sin(w*t1));
        t1 = t - (mu1*l); 
        Y21(I) = Y21(I) - ((psi2p*sin((w*t1)+phase)) + (psi2c*sin(w*t1)));        
        Y21(I) = Y21(I) * (1 - ((l/l2)*(mu3 - mu1))); 
        
    end
    
    
    %calculation of Y32

    if (t - l3 < delta)
        Y32(I)= 0; 
    else
        t1 = t - (mu1*l) - l3; 
        Y32(I) = (psi3p*sin((w*t1) + phase)) + (psi3c*sin(w*t1));
        t1 = t - (mu2*l); 
        Y32(I) = Y32(I) - ((psi3p*sin((w*t1)+phase)) + (psi3c*sin(w*t1)));        
        Y32(I) = Y32(I) * (1 - ((l/l3)*(mu1 - mu2))); 
        
    end
    %calculation of Y13

    if (t - l1 < delta)
        Y13(I) = 0; 
    else
        t1 = t - (mu2*l) - l1; 
        Y13(I) = (psi1p*sin((w*t1) + phase)) + (psi1c*sin(w*t1));
        t1 = t - (mu3*l); 
        Y13(I) = Y13(I) - ((psi1p*sin((w*t1)+phase)) + (psi1c*sin(w*t1)));        
        Y13(I) = Y13(I) * (1 - ((l/l1)*(mu2 - mu3))); 
        
    end
    
    
    %calculation of Y31
    
    if (t - l3 < delta) 
        Y31(I) = 0; 
    else
        t1 = t - (mu2*l) - l3; 
        Y31(I) = (psi3p*sin((w*t1)+phase)) + (psi3c*sin(w*t1)); 
        t1 = t - (mu1*l);
        Y31(I) = Y31(I) - ((psi3p*sin((w*t1)+phase)) + (psi3c*(sin(w*t1)))); 
        Y31(I) = Y31(I) * (1 + ((l/l3)*(mu1 - mu2)));         
    end
    
    
   %calculation of Y12
    
    if (t - l1 < delta) 
        Y12(I) = 0; 
    else
        t1 = t - (mu3*l) - l1; 
        Y12(I) = (psi1p*sin((w*t1)+phase)) + (psi1c*sin(w*t1)); 
        t1 = t - (mu2*l);
        Y12(I) = Y12(I) - ((psi1p*sin((w*t1)+phase)) + (psi1c*(sin(w*t1)))); 
        Y12(I) = Y12(I) * (1 + ((l/l1)*(mu2 - mu3)));         
        
    end
    
   %calculation of Y23
    
    if (t - l2 < delta) 
        Y23(I) = 0; 
    else
        t1 = t - (mu1*l) - l2; 
        Y23(I) = (psi2p*sin((w*t1)+phase)) + (psi2c*sin(w*t1)); 
        t1 = t - (mu3*l);
        Y23(I) = Y23(I) - ((psi2p*sin((w*t1)+phase)) + (psi2c*(sin(w*t1)))); 
        Y23(I) = Y23(I) * (1 + ((l/l2)*(mu3 - mu1)));         
        
    end
    
end

% Round offs due to rate of data production

l3 = delta * floor(l3*(1/delta)); 
l2 = delta * floor(l2*(1/delta)); 
l1 = delta * floor(l1*(1/delta)); 



%calculation of combination alpha

for I = 1:N1
    t = I*delta; 
    
    if ( t - 2*(l1 + l2 + l3 ) < delta) 
        X(I) = 0; 
    else
        X(I) = Y21(I); 
        X(I) = X(I) - Y31(I);        
        X(I) = X(I) + Y13((t-l2)*(1/delta));        
        X(I) = X(I) - Y12((t-l3)*(1/delta));        
        X(I) = X(I) + Y32((t-l1-l2)*(1/delta));   
        X(I) = X(I) - Y23((t - l1-l3)*(1/delta));   
    end
end

%calculation for combination X

for I = 1:N1
    t = I*delta; 
    
    if ( (t - 2*(l2 + l2 + l3 +l3)) < delta) 
        X1(I) = 0; 
    else
        X1(I) = Y32 ( (t - (l3+l2+l2))*(1/delta)); 
        X1(I) = X1(I) - Y23 ( (t - (l2+l3+l3))*(1/delta));        
        X1(I) = X1(I) + Y31((t-l2-l2)*(1/delta));        
        X1(I) = X1(I) - Y21((t-l3-l3)*(1/delta));        
        X1(I) = X1(I) + Y23((t-l2)*(1/delta));   
        X1(I) = X1(I) - Y32((t-l3)*(1/delta));   
        X1(I) = X1(I) + Y21(I); 
        X1(I) = X1(I) - Y31(I); 
    end

end