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

% Total Number of Initial Data 

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

% Initialize Arrays.. Yij, Zij as in Armstrong, Tinto et al.

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

Z21 = zeros(N1, 1); 
Z12 = zeros(N1, 1); 
Z23 = zeros(N1, 1); 
Z32 = zeros(N1,1 ); 
Z31 = zeros(N1, 1); 
Z13 = zeros(N1,1 ); 

% Initialize combination arrays: X = combination Alpha, X1 = combination X

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

%Laser Noise: Spectrum from Tinto et al.

%C1 = (1.064 * (10^(-13)))*randn(N1, 1); 
%C1s = (1.064 * (10^(-13)))*randn(N1, 1); 
%C2 = (1.064 * (10^(-13)))*randn(N1, 1); 
%C2s = (1.064 * (10^(-13)))*randn(N1, 1); 
%C3 = (1.064 * (10^(-13)))*randn(N1, 1); 
%C3s = (1.064 * (10^(-13)))*randn(N1, 1); 

C1 = LaserNoise(N1, delta); 
C1s = LaserNoise(N1, delta); 
C2 = LaserNoise(N1, delta); 
C2s = LaserNoise(N1, delta); 
C3 = LaserNoise(N1, delta); 
C3s = LaserNoise(N1, delta); 


%Shot Noise at detectors: Reference Pre Phase A report as well as Tinto & Armstrong

Y21_shot = ShotNoise(N1, delta);
Y12_shot = ShotNoise(N1, delta); 
Y31_shot = ShotNoise(N1, delta); 
Y13_shot = ShotNoise(N1, delta); 
Y23_shot = ShotNoise(N1, delta); 
Y32_shot = ShotNoise(N1, delta); 

%acceleration noise : Reference Pre Phase A report as well as Tinto & Armstrong

for I = 1:3

    v1(:,I) = AccelNoise(N1, delta); 
    v1s(:,I) = AccelNoise(N1, delta); 
    v2(:, I) = AccelNoise (N1, delta); 
    v2s(:, I) = AccelNoise (N1, delta); 
    v3(:, I) = AccelNoise (N1, delta); 
    v3s(:, I) = AccelNoise (N1, delta); 
end

%Velocity noise: Assumed to have the same spectra as that of the acceleration noise 
%except that its several (6) orders of magnitude larger. This assumption is due to 
%lack of references: 

for I = 1:3

    V1(:, I) = VelocityNoise(N1, delta); 
    V1s(:, I) = VelocityNoise(N1, delta); 
    V2 (:, I) = VelocityNoise (N1, delta); 
    V2s(:, I) = VelocityNoise (N1, delta); 
    V3(:, I) = VelocityNoise (N1, delta); 
    V3s(:, I) = VelocityNoise (N1, delta); 
end


%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))); 

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

% 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)
        C32 = 0; 
        V32 = [0 0 0];
        Y21(I)= 0; 
        Z21(I) = 0; 
    else
        C32 = C3( (t - l2)*(1/delta));
        V32 = V3( (t-l2)*(1/delta), :);
        Y21(I) = C32 - dot(n2, V32) + (2*dot(n2, v1s(I, :))) - dot(n2, V1s(I, :)) - C1s(I) + Y21_shot(I);     
        Z21(I) = C1(I) + (2*dot(n3, (v1(I, :)-V1(I, :)))) - C1s(I); 
        
    end
    
    
    %calculation of Y32

    if (t - l3 < delta)
        C13 = 0; 
        V13 = [0 0 0];
        Y32(I)= 0; 
        Z32(I) = 0; 
    else
        C13 = C1( (t - l3)*(1/delta));
        V13 = V1( (t-l3)*(1/delta), :);
        Y32(I) = C13 - dot(n3, V13) + (2*dot(n3, v2s(I, :))) - dot(n3, V2s(I, :)) - C2s(I) + Y32_shot(I);     
        Z32(I) = C2(I) + (2*dot(n1, (v2(I, :)-V2(I, :)))) - C2s(I); 
    end
    %calculation of Y13

    if (t - l1 < delta)
        C21 = 0; 
        V21 = [0 0 0];
        Y13(I) = 0; 
        Z13(I) = 0; 
    else
        C21 = C2( (t - l1)*(1/delta));
        V21 = V2( (t-l1)*(1/delta), :);
        Y13(I) = C21 - dot(n1, V21) + (2*dot(n1, v3s(I, :))) - dot(n1, V3s(I, :)) - C3s(I) + Y13_shot(I);     
        Z13(I) = C3(I) + (2*dot(n2, (v3(I, :)-V3(I, :)))) - C3s(I); 
    end
    
    
    %calculation of Y31
    
    if (t - l3 < delta) 
        C23s = 0; 
        V23s = [0 0 0]; 
        Y31(I) = 0; 
        Z31(I) = 0; 
    else
        C23s = C2s ( (t-l3)*(1/delta)); 
        V23s = V2s ( (t-l3)*(1/delta), :); 
        Y31(I) = C23s + dot(n3, V23s) - (2*dot(n3, v1(I, :))) + dot(n3, V1(I, :)) - C1(I) + Y31_shot(I);    
        Z31(I) = C1s(I) - (2*dot(n2, (v1s(I, :) - V1s(I, :)))) - C1(I); 
    end
    
    
   %calculation of Y12
    
    if (t - l1 < delta) 
        C31s = 0; 
        V31s = [0 0 0]; 
        Y12(I) = 0; 
        Z12(I) = 0; 
    else
        C31s = C3s ( (t-l1)*(1/delta)); 
        V31s = V3s ( (t-l1)*(1/delta), :); 
        Y12(I) = C31s + dot(n1, V31s) - (2*dot(n1, v2(I, :))) + dot(n1, V2(I, :)) - C2(I) + Y12_shot(I);    
        Z12(I) = C2s(I) - (2*dot(n3, (v2s(I, :) - V2s(I, :)))) - C2(I); 

    end
    
   %calculation of Y23
    
    if (t - l2 < delta) 
        C12s = 0; 
        V12s = [0 0 0]; 
        Y23(I) = 0; 
        Z23(I) = 0; 
    else
        C12s = C1s ( (t-l2)*(1/delta)); 
        V12s = V1s ( (t-l2)*(1/delta), :); 
        Y23(I) = C12s + dot(n2, V12s) - (2*dot(n2, v3(I, :))) + dot(n2, V3(I, :)) - C3(I) + Y23_shot(I);    
        Z23(I) = C3s(I) - (2*dot(n1, (v3s(I, :) - V3s(I, :)))) - C3(I); 
    end
    
end

clear C1
clear C1s
clear C2
clear C2s
clear C3
clear C3s
clear V1
clear V1s
clear V2
clear V2s
clear V3
clear V3s
clear v1
clear v1s
clear v2
clear v2s
clear v3
clear v3s
clear Y21_shot
clear Y23_shot
clear Y12_shot
clear Y13_shot
clear Y31_shot
clear Y32_shot




%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));   
        
        X(I) = X(I) - (Z13( (t-l2)*(1/delta))/2); 
        X(I) = X(I) - (Z13( (t - l1 - l3)*(1/delta))/2); 
        X(I) = X(I) - (Z21(I)/2); 
        X(I) = X(I) - (Z21( (t-l1-l2-l3)*(1/delta))/2);
        X(I) = X(I) - (Z32( (t-l3)*(1/delta))/2); 
        X(I) = X(I) - (Z32( (t- l1 -l2)*(1/delta))/2); 
        
        X(I) = X(I) + (Z23( (t-l2)*(1/delta))/2); 
        X(I) = X(I) + (Z23( (t - l1 - l3)*(1/delta))/2); 
        X(I) = X(I) + (Z31(I)/2); 
        X(I) = X(I) + (Z31( (t - l1 - l2 -l3)*(1/delta))/2); 
        X(I) = X(I) + (Z12( (t-l3)*(1/delta))/2); 
        X(I) = X(I) + (Z12( (t-l1 -l2)*(1/delta))/2); 

    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); 
        
        X1(I) = X1(I) - (Z21((t - l2 - l2 -l3 -l3)*(1/delta))/2); 
        X1(I) = X1(I) + (Z21((t - l3 - l3)*(1/delta))/2); 
        X1(I) = X1(I) + (Z21((t-l2 -l2)*(1/delta))/2); 
        X1(I) = X1(I) - (Z21(I)/2); 
        
        X1(I) = X1(I) + ( Z31((t - l2 - l2 -l3 -l3)*(1/delta))/2); 
        X1(I) = X1(I) - (Z31( (t-l3-l3)*(1/delta))/2); 
        X1(I) = X1(I) - (Z31( (t-l2-l2)*(1/delta))/2); 
        X1(I) = X1(I) + (Z31(I)/2); 
    end

end