% Construct the relevant objects for C(k,k',theta)
function [C0_matrix C2_matrix C4_matrix lambda U Ratio] = FittedCov(k_survey, ps, N_over_2)


[dummy nbins] = size(k_survey');

%---------------------------------------------------------------
% 1- Construct off-diagonal Eigenvectors U(ell,lambda,k)

U = zeros(3,4,nbins);
lambda = zeros(3,4);


%ell = 0, 1 or 2 Eigenvectors
U(1,1,:) = FitU1_new2([0.0738 0.0505 0.6416 1.3145], k_survey);
U(1,2,:) = FitU2([0.1676 2.7257 0.5677 0.0955], k_survey);


%ell = 2, 4 Eigenvectors
U(2,1,:) = FitU1_new2([0.2731 0.8266 1.962 0.816], k_survey);
U(2,2,:) = FitU2([0.15772      2.4207     0.79153    0.032207], k_survey);
U(2,3,:) = FitU2([0.14426      4.0613     0.76611    -0.26272], k_survey);
U(2,4,:) = FitU2([0.14414       5.422     0.84826     0.31324], k_survey);

%ell = 4, 4 Eigenvectors
U(3,1,:) = FitU1_new2([0.060399     0.10344     0.64008      2.2584], k_survey);
U(3,2,:) = FitU2([0.1553    2.3370    0.9307   -0.1154], k_survey);
U(3,3,:) = FitU2([0.1569      3.6937     0.92304    0.040061], k_survey);
U(3,4,:) = FitU2([0.15233      5.1617      0.8899    -0.14503], k_survey);

%-----------------------

lambda(1,1) = 61.9323;
lambda(1,2) = 0;%2.0414;

lambda(2,1) = 35.7400;
lambda(2,2) = 4.4144;
lambda(2,3) = 1.7198;
lambda(2,4) = 0.9997;

lambda(3,1) = 22.0881;
lambda(3,2) = 4.5984;
lambda(3,3) = 2.2025;
lambda(3,4) = 1.4062;


%---------------------------------------------------------------
%2- Construct the normalized C_ell matrices

rho_0 = zeros(nbins);
rho_2 = zeros(nbins);
rho_4 = zeros(nbins);
for i = 1:nbins
    for j = 1:nbins
       for eigen = 1:4
           rho_0(i,j)  = rho_0(i,j)+ lambda(1,eigen)*U(1,eigen,i)*U(1,eigen,j); 
           rho_2(i,j)  = rho_2(i,j)+ lambda(2,eigen)*U(2,eigen,i)*U(2,eigen,j); 
           rho_4(i,j)  = rho_4(i,j)+ lambda(3,eigen)*U(3,eigen,i)*U(3,eigen,j);             
       end
        
       if(i==j)
           rho_0(i,j) = 1;
           rho_2(i,j) = 1;
           rho_4(i,j) = 1;
       end
    end 
end

size(rho_0)


%------------------------------------------------
%3- Get the diagonal elements: diag = Ratio*Gauss

Ratio=zeros(3,nbins);
diag=zeros(3,nbins);

Gauss_cov = (ps.^2./N_over_2*4*pi)';

Ratio(1,:)  = FitRatio([0.2095 1.9980],k_survey) +1;
Ratio(2,:)  = FitRatio([0.555 1.2224],k_survey) +1;
Ratio(3,:)  = FitRatio([1.6525 1.0674],k_survey) +1;

diag(1,:) = Gauss_cov.*Ratio(1,:)';
diag(2,:) = Gauss_cov.*Ratio(2,:)';
diag(3,:) = Gauss_cov.*Ratio(3,:)';

%--------------------------------------
%5-Construct the un-normalized C_ell matrices
C0_matrix = zeros(nbins);
C2_matrix = zeros(nbins);
C4_matrix = zeros(nbins);

for i = 1:nbins
    for j = 1:nbins

        C0_matrix(i,j)  = rho_0(i,j)*sqrt(diag(1,i)*diag(1,j));
        C2_matrix(i,j)  = rho_2(i,j)*sqrt(diag(2,i)*diag(2,j));
        C4_matrix(i,j)  = rho_4(i,j)*sqrt(diag(3,i)*diag(3,j));        
        
    end
end


