%%  Aysajan Abidin
%% Assemble_noise.m 
%% Noise assembly code in 1D

function L = Assemble_noise(p)
global  h x_min x_max y_min y_max
%% Compute the Cholesky factor. First, compute Q_tilda = (Q*phi_i, phi_j),
%% where Q is the covariance operator. For the double integration the 
%% built-in function dblquad is used. For details type "help dblquad" in command window.
 Np = length(p);
 Q_tilda = zeros(Np);

%Just chage 1 in h and P if the spatial interval is not (0,1)
 h = 1/(Np+1);
 P = 0:h:1;
 
 for i = 1:Np     
     
     x_min = P(i); 
     x_max = p(i);    
     
     for j = 1:Np
 
        y_min = P(j); 
        y_max = p(j);
        Q_tilda(i, j) = 1/(h^2)*(dblquad(@integrand_1,x_min,x_max,y_min,y_max) +...
            dblquad(@integrand_2,x_max,x_max+h,y_max,y_max+h));
        Q_tilda(j, i) = Q_tilda(i, j);
        
     end                
 end
 U = chol(Q_tilda); 
 L = U';
