%%  Function that constructs artificial dissipation operators using undivided differences
function D = dissipationOperator(m, order, Hinv, scaling)
    % TBD: Add or remove D_2 and Dp/Dm?
    % d2=[1 2 1];
    % D_2=(diag(ones(m-1,1),-1)-2*diag(ones(m,1),0)+ ...
    % diag(ones(m-1,1),1));
    % D_2(1,1:3)=[d2];D_2(m,m-2:m)=[d2];
    % %Dm
    % DD_m=(diag(ones(m-1,1),+1)-diag(ones(m,1),0));
    % DD_m(m,m-1:m)=[-1 1];
    % %Dp
    % DD_p=(-diag(ones(m-1,1),-1)+diag(ones(m,1),0));
    % DD_p(1,1:2)=[-1 1];

    switch order
        case 1
            DD_1=(diag(ones(m-1,1),+1)-diag(ones(m,1),0));
            DD_1(m,m-1:m)=[-1 1];
            D = DD_2'*DD_2;
        case 2
            dd2=0*[1 -2 1];
            DD_2=(diag(ones(m-1,1),-1)-2*diag(ones(m,1),0)+ ...
            	  diag(ones(m-1,1),1));
            DD_2(1,1:3)=[dd2];DD_2(m,m-2:m)=[dd2];
            D = DD_2'*DD_2;
        case 3
            d3=0*[-1 3 -3 1];
            DD_3=(-diag(ones(m-2,1),-2)+3*diag(ones(m-1,1),-1)-3*diag(ones(m,1),0)+ ...
                  diag(ones(m-1,1),1));
            DD_3(1:2,1:4)=[d3;d3];
            DD_3(m,m-3:m)=[d3];
            D = DD_3'*DD_3;
        case 4
            default_arg('scaling', 1/12);
            d4=0*[1 -4 6 -4 1];
            DD_4=(diag(ones(m-2,1),2)-4*diag(ones(m-1,1),1)+6*diag(ones(m,1),0)-4*diag(ones(m-1,1),-1)+diag(ones(m-2,1),-2));
            DD_4(1:2,1:5)=[d4;d4];DD_4(m-1:m,m-4:m)=[d4;d4];
            D = DD_4'*DD_4;
        case 5
            d5=0*[-1 5 -10 10 -5 1];
            DD_5=(-diag(ones(m-3,1),-3)+5*diag(ones(m-2,1),-2)-10*diag(ones(m-1,1),-1)+10*diag(ones(m,1),0)-5*diag(ones(m-1,1),1)+diag(ones(m-2,1),2));
            DD_5(1:3,1:6)=[d5;d5;d5];
            DD_5(m-1:m,m-5:m)=[d5;d5];
            D = DD_5'*DD_5; 
        case 6
            default_arg('scaling', 1/60);
            d6=0*[1 -6 15 -20 15 -6 1];
            DD_6=(diag(ones(m-3,1),3)-6*diag(ones(m-2,1),2)+15*diag(ones(m-1,1),1)-20*diag(ones(m,1),0)+15*diag(ones(m-1,1),-1)-6*diag(ones(m-2,1),-2)+diag(ones(m-3,1),-3));
            DD_6(1:3,1:7)=[d6;d6;d6];DD_6(m-2:m,m-6:m)=[d6;d6;d6];
            D = DD_6'*DD_6;
        case 7
            d7=0*[-1 7 -21 35 -35 21 -7 1]; 
            DD_7=(-diag(ones(m-4,1),-4)+7*diag(ones(m-3,1),-3)-21*diag(ones(m-2,1),-2)+35*diag(ones(m-1,1),-1)-35*diag(ones(m,1),0)+21*diag(ones(m-1,1),1)-7*diag(ones(m-2,1),2)+diag(ones(m-3,1),3));
            DD_7(1:4,1:8)=[d7;d7;d7;d7];
            DD_7(m-2:m,m-7:m)=[d7;d7;d7];
            D = DD_7'*DD_7;
        case 9
            d9=0*[-1 9 -36 84 -126 126 -84 36 -9 1]; 
            DD_9=(-diag(ones(m-5,1),-5)+9*diag(ones(m-4,1),-4)-36*diag(ones(m-3,1),-3)+84*diag(ones(m-2,1),-2)-126*diag(ones(m-1,1),-1)+126*diag(ones(m,1),0)-84*diag(ones(m-1,1),1)+36*diag(ones(m-2,1),2)-9*diag(ones(m-3,1),3)+diag(ones(m-4,1),4));
            DD_9(1:5,1:10)=[d9;d9;d9;d9;d9];
            DD_9(m-3:m,m-9:m)=[d9;d9;d9;d9];
            D = DD_9'*DD_9;
        otherwise
            error('Order not yet supported', order);
    end
    D = scaling*sparse(Hinv*D);
end