classdef D1Staggered < sbp.OpSet
    properties
        % x_primal: "primal" grid with m points. Equidistant. Called Plus grid in Ossian's paper.
        % x_dual: "dual" grid with m+1 points. Called Minus grid in Ossian's paper.

        % D1_primal takes FROM dual grid TO primal grid
        % D1_dual takes FROM primal grid TO dual grid

        D1_primal % SBP operator approximating first derivative
        D1_dual % SBP operator approximating first derivative
        H_primal % Norm matrix
        H_dual % Norm matrix
        H_primalI % H^-1
        H_dualI % H^-1
        e_primal_l % Left boundary operator
        e_dual_l % Left boundary operator
        e_primal_r % Right boundary operator
        e_dual_r % Right boundary operator
        m % Number of grid points.
        m_primal % Number of grid points.
        m_dual % Number of grid points.
        h % Step size
        x_primal % grid
        x_dual % grid
        x
        borrowing % Struct with borrowing limits for different norm matrices
    end

    methods
        function obj = D1Staggered(m,lim,order)

          x_l = lim{1};
          x_r = lim{2};
          L = x_r-x_l;

          m_primal = m;
          m_dual = m+1;

          switch order
          case 2
            [x_primal, x_dual, Pp, Pm, Qp, Qm] = sbp.implementations.d1_staggered_2(m, L);
          case 4
            [x_primal, x_dual, Pp, Pm, Qp, Qm] = sbp.implementations.d1_staggered_4(m, L);
          case 6
            [x_primal, x_dual, Pp, Pm, Qp, Qm] = sbp.implementations.d1_staggered_6(m, L);
          otherwise
           error('Invalid operator order %d.',order);
          end

          obj.m = m;
          obj.m_primal = m_primal;
          obj.m_dual = m_dual;
          obj.x_primal = x_l + x_primal';
          obj.x_dual = x_l + x_dual';

          D1_primal = Pp\Qp;
          D1_dual = Pm\Qm;

          obj.D1_primal = D1_primal;
          obj.D1_dual = D1_dual;
          obj.H_primal = Pp;
          obj.H_dual = Pm; 

          obj.e_primal_l = sparse(m_primal,1);
          obj.e_primal_r = sparse(m_primal,1);
          obj.e_primal_l(1) = 1;
          obj.e_primal_r(m_primal) = 1;

          obj.e_dual_l = sparse(m_dual,1);
          obj.e_dual_r = sparse(m_dual,1);
          obj.e_dual_l(1) = 1;
          obj.e_dual_r(m_dual) = 1;

          obj.H_primalI = inv(obj.H_primal);
          obj.H_dualI = inv(obj.H_dual);

          obj.borrowing = [];
          obj.x = [];

        end

        function str = string(obj)
            str = [class(obj) '_' num2str(obj.order)];
        end
    end
end