--- a/+scheme/Schrodinger2dCurve.m	Sat Feb 25 12:44:01 2017 +0100
+++ b/+scheme/Schrodinger2dCurve.m	Fri Mar 03 16:19:41 2017 +0100
@@ -134,8 +134,8 @@
                 ti_tau = parametrization.Ti.points(obj.p{5}(t),obj.p{6}(t),obj.p{7}(t),obj.p{8}(t));
                 
                 lcoords=points(obj.grid);
-                [obj.xm,obj.ym]= ti.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2));
-                [x_tau,y_tau]= ti_tau.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2));
+                [obj.xm,obj.ym]= ti.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_v,2));
+                [x_tau,y_tau]= ti_tau.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_v,2));
                 x = reshape(obj.xm,obj.m_tot,1);
                 y = reshape(obj.ym,obj.m_tot,1);
                 obj.x = x;
@@ -198,55 +198,59 @@
         %       neighbour_scheme    is an instance of Scheme that should be interfaced to.
         %       neighbour_boundary  is a string specifying which boundary to interface to.
         function [closure, penalty] = boundary_condition(obj, boundary,~)
-                    [e, d_n, d_t, coeff_t, coeff_n s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = obj.get_boundary_ops(boundary);
+                    [e, d_n, d_t, coeff_t, coeff_n s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g,I,Ji] = obj.get_boundary_ops(boundary);
                  
-                    a_t = spdiag(coeff_t);
-                    a_n = spdiag(coeff_n);
-                    F = (s  * a_n *d_n' + s * a_t*d_t')';
+                    a_t =  @(t) spdiag(coeff_t(t));
+                    a_n = @(t) spdiag(coeff_n(t));
+                    Ji = spdiag(Ji);
+                    
+                    F = @(t)(s * (d_n * Ji)' + s * a_t(t)  *d_t')';
                     tau1  = 1;       
-                    a = spdiag(g);
-                    tau2 =  (-1*s*a - abs(a))/4;
+                    a  = @(t)spdiag(g(t));
+                    tau2 = @(t) (-1*s*a(t) - abs(a(t)))/4;
 
-                    penalty_parameter_1 = @(t) 1*1i*halfnorm_inv_n*halfnorm_inv_t*F*e'*halfnorm_t*e;
-                    penalty_parameter_2 = @(t) halfnorm_inv_n*e*tau2;
+                    penalty_parameter_1 = @(t) 1*1i*halfnorm_inv_n*halfnorm_inv_t*F(t)*e'*halfnorm_t*e;
+                    penalty_parameter_2 = @(t) halfnorm_inv_n*e*tau2(t);
 
                     closure = @(t) obj.Ji*obj.c^2 * penalty_parameter_1(t)*e' + obj.Ji* penalty_parameter_2(t)*e';
-                    penalty = @(t) -obj.Ji*obj.c^2 * penalty_parameter_1(t)*e'+  obj.Ji*penalty_parameter_2(t)*e';
+                    penalty = @(t) -obj.Ji*obj.c^2 * penalty_parameter_1(t)*e'-  obj.Ji*penalty_parameter_2(t)*e';
                 
         end
         
         function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
             % u denotes the solution in the own domain
             % v denotes the solution in the neighbour domain
-            [e_u, d_n_u, d_t_u,  coeff_t_u, coeff_n_u,s_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t,gamm_u, I_u] = obj.get_boundary_ops(boundary);
-            [e_v, d_n_v, d_t_v,  coeff_t_v, coeff_n_v s_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t,gamm_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
+            [e_u, d_n_u, d_t_u,  coeff_t_u, coeff_n_u,s_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t,gamm_u, I_u,JI_u] = obj.get_boundary_ops(boundary);
+            [e_v, d_n_v, d_t_v,  coeff_t_v, coeff_n_v s_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t,gamm_v, I_v,JI_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
 
-            a_n_u = spdiag(coeff_n_u);
-            a_t_u = spdiag(coeff_t_u);
-            a_n_v = spdiag(coeff_n_v);
-            a_t_v = spdiag(coeff_t_v);
+            a_n_u = @(t) spdiag(coeff_n_u(t));
+            a_t_u = @(t) spdiag(coeff_t_u(t));
+            a_n_v = @(t) spdiag(coeff_n_v(t));
+            a_t_v = @(t) spdiag(coeff_t_v(t));
+            
+            Ji_u = spdiag(JI_u);
+            Ji_v =  spdiag(JI_v);
 
-            F_u = (s_u * a_n_u * d_n_u' + s_u * a_t_u*d_t_u')';
-            F_v = (s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v')';
+            F_u = @(t)((d_n_u*Ji_u)' + a_t_u(t)*d_t_u')';
+            F_v = @(t)((d_n_v*Ji_v)' + a_t_v(t)*d_t_v')';
             
-            a = spdiag(gamm_u);
-                  
-            u = obj;
-            v = neighbour_scheme;
+            a = @(t)spdiag(gamm_u(t));
+       
+            tau =  s_u*1*1i/2;            
+            sig =  -s_u*1*1i/2;
+            gamm = @(t) (-s_u*a(t))/2;
+            
+            penalty_parameter_1 = @(t) halfnorm_inv_u_n*(tau*halfnorm_inv_u_t*F_u(t)*e_u'*halfnorm_u_t*e_u);
+            penalty_parameter_2 = @(t) halfnorm_inv_u_n * e_u  * (sig );
+            penalty_parameter_3 = @(t) halfnorm_inv_u_n * e_u  * (gamm(t) );
 
-            tau = -1/2*1i;            
-            sig = 1/2*1i;
-            gamm =  (-1*s_u*a)/2;
-            
-            penalty_parameter_1 =@(t) halfnorm_inv_u_n*(e_u*tau + sig*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t*e_u);
-            penalty_parameter_2 =@(t) halfnorm_inv_u_n * e_u  * (-sig + gamm );
 
-            closure =@(t) obj.Ji*obj.c^2 * ( penalty_parameter_1(t)*e_u' + penalty_parameter_2(t)*F_u');
-            penalty =@(t) obj.Ji*obj.c^2 * (-penalty_parameter_1(t)*e_v' + penalty_parameter_2(t)*F_v');
+            closure =@(t) obj.Ji*obj.c^2 * ( penalty_parameter_1(t)*e_u' + penalty_parameter_2(t)*F_u(t)' + penalty_parameter_3(t)*e_u');
+            penalty =@(t) obj.Ji*obj.c^2 * ( -penalty_parameter_1(t)*e_v' - penalty_parameter_2(t)*F_v(t)' - penalty_parameter_3(t)*e_v');
         end
 
 
-        function [e, d_n, d_t, coeff_t,coeff_n, s,  halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g, I] = get_boundary_ops(obj, boundary)
+        function [e, d_n, d_t, coeff_t,coeff_n, s,  halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g, I,Ji] = get_boundary_ops(obj, boundary)
 
             % gridMatrix = zeros(obj.m(2),obj.m(1));
             % gridMatrix(:) = 1:numel(gridMatrix);
@@ -261,9 +265,9 @@
                     s = -1;
 
                     I = ind(1,:);
-                    coeff_t = obj.a12(I);
-                    coeff_n = obj.a11(I);
-                    g=obj.g_1(I);
+                    coeff_t = @(t)obj.a12(I);
+                    coeff_n =@(t) obj.a11(I);
+                    g = @(t)obj.g_1(I);
                 case 'e'
                     e = obj.e_e;
                     d_n = obj.du_e;
@@ -271,9 +275,9 @@
                     s = 1;
 
                     I = ind(end,:);
-                    coeff_t = obj.a12(I);
-                    coeff_n = obj.a11(I);
-                    g=obj.g_1(I);
+                    coeff_t = @(t)obj.a12(I);
+                    coeff_n = @(t)obj.a11(I);
+                    g = @(t)obj.g_1(I);
                 case 's'
                     e = obj.e_s;
                     d_n = obj.dv_s;
@@ -281,9 +285,9 @@
                     s = -1;
 
                     I = ind(:,1)';
-                    coeff_t = obj.a12(I);
-                    coeff_n = obj.a11(I);
-                    g=obj.g_2(I);
+                    coeff_t = @(t)obj.a12(I);
+                    coeff_n = @(t)obj.a11(I);
+                    g = @(t)obj.g_2(I);
                 case 'n'
                     e = obj.e_n;
                     d_n = obj.dv_n;
@@ -291,9 +295,9 @@
                     s = 1;
 
                     I = ind(:,end)';
-                    coeff_t = obj.a12(I);
-                    coeff_n = obj.a11(I);
-                    g=obj.g_2(I);
+                    coeff_t = @(t)obj.a12(I);
+                    coeff_n = @(t)obj.a11(I);
+                    g = @(t)obj.g_2(I);
                 otherwise
                     error('No such boundary: boundary = %s',boundary);
             end
@@ -309,10 +313,21 @@
                     halfnorm_inv_t = obj.Hiu;
                     halfnorm_t = obj.Hu;
             end
+            fis = diag(obj.Ji);
+            Ji = fis(I);
         end
 
         function N = size(obj)
             N = prod(obj.m);
         end
     end
+    methods(Static)
+        % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u
+        % and bound_v of scheme schm_v.
+        %   [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l')
+        function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
+            [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
+            [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
+        end
+    end
 end
\ No newline at end of file