function [LL, LL0,LLDecay1,LLDecay2]=L_reg_order2(Memory,order)

%Function that makes the regularization matrix L. 
% 
% Inputs: 
% 
% M = Memory
% order = Order of the model
%Outputs:
%Size of the matrixes is dependent on the model order. 
%LL = 1st-order Tikhonov regularization matrix
% LL0 = 0-th order Tikhonov regularization matrix
% LLDecay1 = Decay regularization matrix 1st-order kernel regularization
% LLDecay2 = Decay regularization matrix 2nd-order kernel regularization
%%
LLDecay1 = NaN;
LLDecay2 = NaN;
%% 1st-order Tikhonov, 1st-order kernel
Upper = triu(-ones((Memory+1)),1)-triu(-ones((Memory+1)),2);
Lower = tril(-ones((Memory+1)),-1)-tril(-ones((Memory+1)),-2);
L1 = Upper + Lower;
DiagL1 = diag(-sum(L1,2));
L1 = L1 + DiagL1;
%% 1st-order Tikhonov, 2nd-order kernel
a = reshape( repmat(1:1:(Memory+1),Memory+1,1),1,[]);
b = repmat(1:1:Memory+1,1,Memory+1);
c = a;
d = b;

L2 = zeros([length(a),length(a)]);

for i = 1:1:length(a)
    for j = 1:1:length(a)
        if j==i
            L2(j,i)=0;
        elseif abs(a(i)-c(j))<=1 && abs(b(i)-d(j)) <=1 && (abs(a(i)-c(j))+abs(b(i)-d(j)))<=1
            L2(j,i)=-1;
        end
    end
end

DiagLL = diag(-sum(L2,2));
L2 = L2 + DiagLL;


%% 1st-order Tikhonov, 3rd-order kernel
if order==3

    aa = repmat(a,1,Memory+1);
    bb = repmat(b,1,Memory+1);
    ee = reshape( repmat(1:1:(Memory+1),(Memory+1)^2,1),1,[]);

    cc = aa;
    dd = bb;
    ff = ee;

    L3 = zeros([length(aa),length(aa)]);

    for ii = 1:1:length(aa)
            for jj = 1:1:length(aa)
                if jj==ii
                    L3(jj,ii)=0;
                elseif abs(aa(ii)-cc(jj))<=1 && abs(bb(ii)-dd(jj)) <=1 && abs(ee(ii)-ff(jj))<=1 && (abs(aa(ii)-cc(jj))+abs(bb(ii)-dd(jj))+abs(ee(ii)-ff(jj)))<=1
                    L3(jj,ii)=-1;
                end
            end
    end

    DiagLL3 = diag(-sum(L3,2));
    L3 = L3 + DiagLL3;
end

%% Decay regularization, 1st-order kernel
LLDecay1 = zeros(size(L1));    
    for i = 1:1:length(L1(:,1))
        LLDecay1(i,i) = (exp((5*i/(Memory+1))))-1;
    end    
LLDecay1 = LLDecay1.*(1/max(max(LLDecay1)));

%% Decay regularization, 2nd-order kernel
if order ==2        
    LLDecay2 = zeros([length(a),length(a)]);
    for i = 1:1:length(a)
        for j = 1:1:length(a)
            if j==i
                jj = max(a(i),b(j));
                LLDecay2(j,i)=(exp(5*jj/(Memory+1)))-1;
            end
        end
    end 
end
LLDecay2 = LLDecay2.*(1/(max(max(LLDecay2))));

%% Connect all together: First and second order kernel regularization matrized 
%is added together block diagonal to make the entire regularization matrixes.

if order ==1
    LL = L1;
elseif order ==2
    LL = blkdiag(L1,L2);

    LL10 = zeros(size(L1));
    LL20 = zeros(size(L2));

    LLDecay1 = blkdiag(LLDecay1, LL20); 
    LLDecay2 = blkdiag(LL10, LLDecay2); 
end
    
LL0 = eye(size(LL));   %0th-order Tikonov regularization, for all kernels orders
end


