--- a/+scheme/Staggered1DAcousticsVariable.m	Mon Jun 18 16:33:53 2018 -0700
+++ b/+scheme/Staggered1DAcousticsVariable.m	Mon Jun 18 16:36:24 2018 -0700
@@ -218,26 +218,106 @@
                 
          end
           
+         % Uses the hat variable method for the interface coupling,
+         % see hypsyst_varcoeff.pdf in the hypsyst repository.
+         % Notation: u for left side of interface, v for right side. 
          function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
 
-            error('Staggered1DAcoustics, interface not implemented');
+            % Get boundary matrices
+            switch boundary
+                case 'l'
+                    A_v = obj.A_l;
+                    B_v = obj.B_l;
+                    Hi_v = obj.Hi;
+                    e_v = obj.e_l;
+
+                    A_u = neighbour_scheme.A_r;
+                    B_u = neighbour_scheme.B_r;
+                    Hi_u = neighbour_scheme.Hi;
+                    e_u = neighbour_scheme.e_r;
+                case 'r'
+                    A_u = obj.A_r;
+                    B_u = obj.B_r;
+                    Hi_u = obj.Hi;
+                    e_u = obj.e_r;
+
+                    A_v = neighbour_scheme.A_l;
+                    B_v = neighbour_scheme.B_l;
+                    Hi_v = neighbour_scheme.Hi;
+                    e_v = neighbour_scheme.e_l;
+            end
+            C_u = inv(A_u)*B_u;
+            C_v = inv(A_v)*B_v;
+
+            % Diagonalize C 
+            [T_u, ~] = eig(C_u);
+            Lambda_u = inv(T_u)*C_u*T_u;
+            lambda_u = diag(Lambda_u);
+            S_u = inv(T_u);
+
+            [T_v, ~] = eig(C_v);
+            Lambda_v = inv(T_v)*C_v*T_v;
+            lambda_v = diag(Lambda_v);
+            S_v = inv(T_v);
+
+            % Identify in- and outgoing characteristic variables
+            Im_u = lambda_u < 0;
+            Ip_u = lambda_u > 0;
+            Im_v = lambda_v < 0;
+            Ip_v = lambda_v > 0;
+
+            Tp_u = T_u(:,Ip_u);
+            Tm_u = T_u(:,Im_u);
+            Sp_u = S_u(Ip_u,:);
+            Sm_u = S_u(Im_u,:);
+
+            Tp_v = T_v(:,Ip_v);
+            Tm_v = T_v(:,Im_v);
+            Sp_v = S_v(Ip_v,:);
+            Sm_v = S_v(Im_v,:);
+
+            % Create S_tilde and T_tilde
+            Stilde = [Sp_u; Sm_v];
+            Ttilde = inv(Stilde);
+            Ttilde_p = Ttilde(:,1); 
+            Ttilde_m = Ttilde(:,2);
+
+            % Solve for penalty parameters theta_1,2
+            LHS = Ttilde_m*Sm_v*Tm_u;
+            RHS = B_u*Tm_u;
+            th_u = RHS./LHS; 
+            TH_u = diag(th_u);
+
+            LHS = Ttilde_p*Sp_u*Tp_v;
+            RHS = -B_v*Tp_v;
+            th_v = RHS./LHS; 
+            TH_v = diag(th_v);
+
+            % Construct penalty matrices
+            Z_u = TH_u*Ttilde_m*Sm_v;
+            Z_u = sparse(Z_u);
+
+            Z_v = TH_v*Ttilde_p*Sp_u;
+            Z_v = sparse(Z_v);
+
+            closure_u = Hi_u*e_u*Z_u*e_u';
+            penalty_u = -Hi_u*e_u*Z_u*e_v';
+            closure_v = Hi_v*e_v*Z_v*e_v';
+            penalty_v = -Hi_v*e_v*Z_v*e_u';
 
              switch boundary
-                 % Upwind coupling
-                 case {'l','left'}
-                     tau = -1*obj.e_l;
-                     closure = obj.Hi*tau*obj.e_l';       
-                     penalty = -obj.Hi*tau*neighbour_scheme.e_r';
-                 case {'r','right'}
-                     tau = 0*obj.e_r;
-                     closure = obj.Hi*tau*obj.e_r';       
-                     penalty = -obj.Hi*tau*neighbour_scheme.e_l';
+                 case 'l'
+                     closure = closure_v;
+                     penalty = penalty_v;
+                 case 'r'
+                     closure = closure_u;
+                     penalty = penalty_u;
              end
                  
          end
       
         function N = size(obj)
-            N = obj.m;
+            N = 2*obj.m+1;
         end
 
     end