diff --git a/chunkie/+chnk/+axissymhelm2d/kern.m b/chunkie/+chnk/+axissymhelm2d/kern.m
index fc556404..213fa488 100644
--- a/chunkie/+chnk/+axissymhelm2d/kern.m
+++ b/chunkie/+chnk/+axissymhelm2d/kern.m
@@ -68,7 +68,9 @@
     zk2 = zks(2);
 end
 
-if strcmpi(type, 'd')
+switch lower(type)
+
+case {'d', 'double'}
     srcnorm = srcinfo.n(:,:);
     [~, grad] = chnk.axissymhelm2d.green(zk, src, targ, origin);
     nx = repmat(srcnorm(1,:), nt, 1);
@@ -91,9 +93,8 @@
             end
         end
     end 
-end
 
-if strcmpi(type, 'sprime')
+case {'sp', 'sprime'}
     targnorm = targinfo.n(:,:);
     [~, grad] = chnk.axissymhelm2d.green(zk, src, targ, origin);
 
@@ -118,9 +119,7 @@
         end
     end
 
-end
-
-if strcmpi(type, 's')
+case {'s', 'single'}
     submat = chnk.axissymhelm2d.green(zk, src, targ, origin);
     fker = @(x, s, t) fslp(x, zk, s, t, origin);
     for j=1:ns
@@ -140,15 +139,11 @@
         end
     end
     
-end
-
-
-if strcmpi(type, 'sdiff')
+case {'sdiff', 's_diff'}
     ifdiff = 1;
     submat = chnk.axissymhelm2d.green(zk, src, targ, origin, ifdiff);
-end
 
-if strcmpi(type, 'c')
+case {'c', 'combined'}
     srcnorm = srcinfo.n(:,:); 
     coef = ones(2,1);
     if (nargin == 6); coef = varargin{1}; end
@@ -175,10 +170,8 @@
             end
         end
     end
-end
 
-
-if strcmpi(type, 'dprime')
+case {'dp', 'dprime'}
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
   [~,~,hess] = chnk.axissymhelm2d.green(zk, src, targ, origin);
@@ -188,10 +181,8 @@
   nytarg = repmat((targnorm(2,:)).',1,ns);
   submat = hess(:,:,4).*nxsrc.*nxtarg - hess(:,:,5).*nysrc.*nxtarg ... 
               - hess(:,:,6).*nxsrc.*nytarg + hess(:,:,3).*nysrc.*nytarg;
-end
-
-if strcmpi(type, 'dprime_re_diff')
 
+case {'dp_re_diff', 'dprime_re_diff'}
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
   ifdiff = 2;
@@ -213,11 +204,7 @@
 
   submat = submat1-submat2;
 
-end
-
-
-
-if strcmpi(type, 'dprimediff')
+case {'dpdiff', 'dp_diff', 'dprimediff', 'dprime_diff'}
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
   ifdiff = 1;
@@ -246,9 +233,7 @@
         end
     end
 
-end
-
-if strcmpi(type, 'neu_rpcomb')
+case {'neu_rpcomb'}
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
   [~,gk,~,sikmat,gik,~,~,~,hessdiff] = ...
@@ -314,14 +299,8 @@
   submat(1:3:end, 3:3:end) = c2*spikmat;
   submat(2:3:end, 1:3:end) = -sikmat;
   submat(3:3:end, 1:3:end) = -spikmat;
-  
-  
-  
-
-
-
-
-    
+otherwise
+    error('Unknown axissymmetric Helmholtz kernel type ''%s''.', type); 
 end
 
 end
diff --git a/chunkie/+chnk/+elast2d/kern.m b/chunkie/+chnk/+elast2d/kern.m
index 075708c0..56388485 100644
--- a/chunkie/+chnk/+elast2d/kern.m
+++ b/chunkie/+chnk/+elast2d/kern.m
@@ -67,14 +67,15 @@
 r2 = x.^2 + y.^2;
 r4 = r2.^2;
 
-if (strcmpi(type,'s'))
+switch lower(type)
+case {'s', 'single'}
     logr = beta*log(r2)/2;
     rirjor2 = (kron(gamma*x./r2,[1 0; 1 0]) + ...
         kron(gamma*y./r2,[0 1; 0 1])).* ...
         (kron(x,[1 1; 0 0]) + kron(y,[0 0; 1 1]));
     mat = kron(logr+gamma/2,eye(2)) + rirjor2;
-end
-if (strcmpi(type,'strac'))
+
+case {'strac'}
     dirx = t.n(1,:); dirx = dirx(:);
     diry = t.n(2,:); diry = diry(:);
     rdotv = dirx(:).*x + diry(:).*y;
@@ -85,8 +86,8 @@
         kron(y.*dirx./r2,[0 -1; 1 0]) + ...
         kron(rdotv./r2,[1 0; 0 1]));
     mat = term1+term2;
-end
-if (strcmpi(type,'sgrad'))
+
+case {'sgrad', 'sg'}
     mat = zeros(4*m,2*n);
     tmp = beta*x./r2;
     mat(1:4:end,1:2:end) = tmp;
@@ -116,8 +117,7 @@
     tmp = gamma*(2*r2.*y-y.^2.*(2*y))./r4;
     mat(4:4:end,2:2:end) = mat(4:4:end,2:2:end) + tmp;
     
-end
-if (strcmpi(type,'d'))
+case {'d', 'double'}
     dirx = s.n(1,:); dirx = dirx(:).';
     diry = s.n(2,:); diry = diry(:).';
     rdotv = x.*dirx + y.*diry;
@@ -128,9 +128,8 @@
         kron(y.*dirx./r2,[0 1; -1 0]) + ...
         kron(rdotv./r2,[1 0; 0 1]));
     mat = -(term1+term2);
-end
 
-if (strcmpi(type,'dalt'))
+case {'dalt'}
     dirx = s.n(1,:); dirx = dirx(:).';
     diry = s.n(2,:); diry = diry(:).';
     rdotv = x.*dirx + y.*diry;
@@ -139,9 +138,8 @@
         (kron(x,[1 1; 0 0]) + kron(y,[0 0; 1 1]));
     term1 = eta*(kron(rdotv./r2,[2 0; 0 2]));
     mat = -(term1+term2);
-end
 
-if (strcmpi(type,'daltgrad'))
+case {'daltgrad', 'daltg'}
     dirx = s.n(1,:); dirx = dirx(:).';
     diry = s.n(2,:); diry = diry(:).';
     rdotv = x.*dirx + y.*diry;
@@ -185,9 +183,7 @@
         -2*eta*(diry./r2 - 2*rdotv.*y./r4);
     mat(4:4:end,2:2:end) = aij_xl;
 
-end
-
-if (strcmpi(type,'dalttrac'))
+case {'dalttrac'}
     dirx = s.n(1,:); dirx = dirx(:).';
     diry = s.n(2,:); diry = diry(:).';
     rdotv = x.*dirx + y.*diry;
@@ -239,7 +235,8 @@
     tmp = mu*(matg(2:4:end,:) + matg(3:4:end,:));
     mat(1:2:end,:) = mat(1:2:end,:) + tmp.*n2 + 2*mu*matg(1:4:end,:).*n1;
     mat(2:2:end,:) = mat(2:2:end,:) + tmp.*n1 + 2*mu*matg(4:4:end,:).*n2;
-
+otherwise
+        error('Unknown elasticity kernel type ''%s''.', type);
 end
 
 end
diff --git a/chunkie/+chnk/+flex2d/kern.m b/chunkie/+chnk/+flex2d/kern.m
index 7eefe85c..d67d77e4 100644
--- a/chunkie/+chnk/+flex2d/kern.m
+++ b/chunkie/+chnk/+flex2d/kern.m
@@ -98,14 +98,13 @@
 
 %%% STANDARD LAYER POTENTIALS
 
-if strcmpi(type, 's') % flexural wave single layer
+switch lower(type)
+case {'s', 'single'} % flexural wave single layer
 
    val = chnk.flex2d.hkdiffgreen(zk,src,targ);  
    submat = 1/(2*zk^2).*val;
 
-end
-
-if strcmpi(type, 'sprime') % normal derivative of flexural wave single layer
+case {'sp', 'sprime'} % normal derivative of flexural wave single layer
 
    targnorm = targinfo.n;
    nxtarg = repmat((targnorm(1,:)).',1,ns);
@@ -114,12 +113,10 @@
    [~,grad] = chnk.flex2d.hkdiffgreen(zk,src,targ);  
    submat = 1/(2*zk^2).*(grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg);
 
-end
-
 %%% CLAMPED PLATE KERNELS
 
 % boundary conditions applied to a point source
-if strcmpi(type, 'clamped_plate_bcs')
+case {'clamped_plate_bcs'}
     nxtarg = targinfo.n(1,:).'; 
     nytarg = targinfo.n(2,:).';  
     submat = zeros(2*nt,ns);
@@ -131,10 +128,9 @@
    
     submat(1:2:end,:) = firstbc;
     submat(2:2:end,:) = secondbc;
-end
 
 % kernels for the clamped plate integral equation
-if strcmpi(type, 'clamped_plate')
+case {'clamped_plate'}
    srcnorm = srcinfo.n;
    srctang = srcinfo.d;
    targnorm = targinfo.n;
@@ -195,10 +191,9 @@
     
   submat(2:2:end,1:2:end) = K21;
   submat(2:2:end,2:2:end) = K22;
-end
 
 % clamped plate kernels for plotting
-if strcmpi(type, 'clamped_plate_eval')
+case {'clamped_plate_eval'}
 
     submat = zeros(nt,2*ns);
 
@@ -226,13 +221,11 @@
     submat(:,1:2:end) = K1;
     submat(:,2:2:end) = K2;
 
-end
-
 
 %%% FREE PLATE KERNELS
 
 % boundary conditions applied to a point source
-if strcmpi(type, 'free_plate_bcs')
+case {'free_plate_bcs'}
     targnorm = targinfo.n;
     targtang = targinfo.d;
     targd2 = targinfo.d2;
@@ -274,10 +267,8 @@
     submat(1:2:end,:) = firstbc;
     submat(2:2:end,:) = secondbc;
 
-end
-
 % kernels for the free plate integral equation 
-if strcmpi(type, 'free_plate')
+case {'free_plate'}
    srcnorm = srcinfo.n;
    srctang = srcinfo.d;
    targnorm = targinfo.n;
@@ -394,10 +385,8 @@
     submat(3:4:end,1:2:end) = K11H;
     submat(4:4:end,1:2:end) = K21H;
 
-end
-
 % free plate kernels used for plotting 
-if strcmpi(type, 'free_plate_eval')           
+case {'free_plate_eval'}
    srcnorm = srcinfo.n;
    srctang = srcinfo.d;
    nu = varargin{1};
@@ -423,12 +412,10 @@
    submat(:,2:3:end) = K1H;
    submat(:,3:3:end) = K2;
 
-end
-
 %%% SUPPORTED PLATE KERNELS
 
 % boundary conditions applied to a point source
-if strcmpi(type, 'supported_plate_bcs')
+case {'supported_plate_bcs'}
     nxtarg = targinfo.n(1,:).'; 
     nytarg = targinfo.n(2,:).';  
     dx = targinfo.d(1,:).';
@@ -450,11 +437,10 @@
 
     submat(1:2:end,:) = firstbc;
     submat(2:2:end,:) = secondbc;
-end
 
 % kernels for the supported plate integral equation that have to be 
 % discretized using log quadrature
-if strcmpi(type, 'supported_plate_log')
+case {'supported_plate_log'}
     srcnorm = srcinfo.n;
     srctang = srcinfo.d;
     srcd2 = srcinfo.d2;
@@ -568,11 +554,9 @@
     submat(2:2:end,1:2:end) = K21; 
     submat(2:2:end,2:2:end) = K22; 
 
-end
-
 % the kernel for the supported plate integral equation that has to be 
 % discretized using smooth quadrature to avoid close evaluations
-if strcmpi(type, 'supported_plate_smooth')
+case {'supported_plate_smooth'}
     srcnorm = srcinfo.n;
     srctang = srcinfo.d;
     srcd2 = srcinfo.d2;
@@ -673,10 +657,8 @@
 
     submat(r2 == 0) = (nu - 1)*(12*kappa(r2 == 0).^3*(nu^2 - nu + 4) + kpp(r2 == 0)*(-5*nu^2 + 4*nu + 33))/(48*pi*(nu - 3)) ; % diagonal replacement
 
-end
-
 % supported plate kernels for plotting
-if strcmpi(type, 'supported_plate_eval')
+case {'supported_plate_eval'}
     srcnorm = srcinfo.n;
     srctang = srcinfo.d;
     srcd2 = srcinfo.d2;
diff --git a/chunkie/+chnk/+helm1d/kern.m b/chunkie/+chnk/+helm1d/kern.m
index f2960354..5fcec276 100644
--- a/chunkie/+chnk/+helm1d/kern.m
+++ b/chunkie/+chnk/+helm1d/kern.m
@@ -70,43 +70,39 @@
 [~,ns] = size(src);
 [~,nt] = size(targ);
 
+switch lower(type)
 % double layer
-if strcmpi(type,'d')
+case {'d', 'double'}
   srcnorm = srcinfo.n(:,:);
   [~,grad] = chnk.helm1d.green(zkE,src,targ);
   nx = repmat(srcnorm(1,:),nt,1);
   ny = repmat(srcnorm(2,:),nt,1);
   submat = -(grad(:,:,1).*nx + grad(:,:,2).*ny);
-end
 
 % normal derivative of single layer
-if strcmpi(type,'sprime')
+case {'sp', 'sprime'}
   targnorm = targinfo.n(:,:);
   [~,grad] = chnk.helm1d.green(zkE,src,targ);
   nx = repmat((targnorm(1,:)).',1,ns);
   ny = repmat((targnorm(2,:)).',1,ns);
 
   submat = (grad(:,:,1).*nx + grad(:,:,2).*ny);
-end
-
 
 % Tangential derivative of single layer
-if strcmpi(type,'stau')
+case {'stau', 'st'}
   targtan = targinfo.d(:,:);
   [~,grad] = chnk.helm1d.green(zkE,src,targ);
   dx = repmat((targtan(1,:)).',1,ns);
   dy = repmat((targtan(2,:)).',1,ns);
   ds = sqrt(dx.*dx+dy.*dy);
   submat = (grad(:,:,1).*dx + grad(:,:,2).*dy)./ds;
-end
 
 % single layer
-if strcmpi(type,'s')
+case {'s', 'single'}
   submat = chnk.helm1d.green(zkE,src,targ);
-end
 
 % normal derivative of double layer
-if strcmpi(type,'dprime')
+case {'dp', 'dprime'}
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
   [~,~,hess] = chnk.helm1d.green(zkE,src,targ);
@@ -116,10 +112,9 @@
   nytarg = repmat((targnorm(2,:)).',1,ns);
   submat = -(hess(:,:,1).*nxsrc.*nxtarg + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)...
       + hess(:,:,3).*nysrc.*nytarg);
-end
 
 % Combined field 
-if strcmpi(type,'c')
+case {'c', 'combined'}
   srcnorm = srcinfo.n(:,:);
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end
@@ -128,17 +123,15 @@
   ny = repmat(srcnorm(2,:),nt,1);
   submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny);
   submat = coef(1)*submatd + coef(2)*submats;
-end
 
 % normal derivative of combined field
-if strcmpi(type,'cprime')
+case {'cp', 'cprime'}
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
   
 
-
   % Get gradient and hessian info
   [~,grad,hess] = chnk.helm1d.green(zkE,src,targ);
 
@@ -155,11 +148,9 @@
   submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg);
 
   submat = coef(1)*submatdp + coef(2)*submatsp;
-end
-
 
 % Dirichlet and neumann data corresponding to combined field
-if strcmpi(type,'c2trans') 
+case {'c2trans'}
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end
   targnorm = targinfo.n(:,:);
@@ -187,12 +178,10 @@
   submat(1:2:2*nt,:) = coef(1)*submatd + coef(2)*submats;
   submat(2:2:2*nt,:) = coef(1)*submatdp + coef(2)*submatsp;
 
-end
-
 
 % all kernels, [c11 D, c12 S; c21 D', c22 S'] 
-if strcmpi(type,'all')
- 
+case {'all','trans_sys','ts'}
+
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
   cc = varargin{1};
@@ -220,11 +209,9 @@
   submat(1:2:2*nt,2:2:2*ns) = submats*cc(1,2);
   submat(2:2:2*nt,1:2:2*ns) = submatdp*cc(2,1);
   submat(2:2:2*nt,2:2:2*ns) = submatsp*cc(2,2);
-end
 
 % Dirichlet data/potential correpsonding to transmission rep
-if strcmpi(type,'trans_rep') 
-
+case {'trans_rep','trep'}
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end;
   srcnorm = srcinfo.n(:,:);
@@ -236,10 +223,9 @@
   submat = zeros(1*nt,2*ns);
   submat(1:1:1*nt,1:2:2*ns) = coef(1)*submatd;
   submat(1:1:1*nt,2:2:2*ns) = coef(2)*submats;
-end
 
 % Neumann data corresponding to transmission rep
-if strcmpi(type,'trans_rep_prime')
+case {'trans_rep_prime','trep_p', 'trans_rep_p'}
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end;
   targnorm = targinfo.n(:,:);
@@ -264,11 +250,9 @@
   submat = zeros(nt,2*ns);
   submat(1:1:1*nt,1:2:2*ns) = coef(1)*submatdp;
   submat(1:1:1*nt,2:2:2*ns) = coef(2)*submatsp;
-end
-
 
 % Gradient correpsonding to transmission rep
-if strcmpi(type,'trans_rep_grad')
+case {'trans_rep_grad','trep_g','trans_rep_g'}
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end;
   
@@ -290,6 +274,8 @@
     
   submat(2:2:2*nt,1:2:2*ns) = -coef(1)*(hess(:,:,2).*nxsrc + hess(:,:,3).*nysrc);
   submat(2:2:2*nt,2:2:2*ns) = coef(2)*grad(:,:,2);
+otherwise
+    error('Unknown 1D Helmholtz kernel type ''%s''.', type);
 end
 
 
diff --git a/chunkie/+chnk/+helm2d/fmm.m b/chunkie/+chnk/+helm2d/fmm.m
index de806a75..a71300d1 100644
--- a/chunkie/+chnk/+helm2d/fmm.m
+++ b/chunkie/+chnk/+helm2d/fmm.m
@@ -97,7 +97,7 @@
             varargout{1} = ( U.gradtarg(1,:).*targinfo.n(1,:) + ...
                              U.gradtarg(2,:).*targinfo.n(2,:) ).';
         otherwise
-            error('CHUNKIE:lap2d:fmm:pot', ...
+            error('CHUNKIE:helm2d:fmm:pot', ...
                 'Potentials not supported for Helmholtz kernel ''%s''.', type);
     end
 end
diff --git a/chunkie/+chnk/+helm2d/kern.m b/chunkie/+chnk/+helm2d/kern.m
index 0b535c40..194b71e2 100644
--- a/chunkie/+chnk/+helm2d/kern.m
+++ b/chunkie/+chnk/+helm2d/kern.m
@@ -72,77 +72,80 @@
 [~,ns] = size(src);
 [~,nt] = size(targ);
 
+switch lower(type)
 % double layer
-if strcmpi(type,'d')
+case {'d', 'double'}
   srcnorm = srcinfo.n(:,:);
   [~,grad] = chnk.helm2d.green(zk,src,targ);
   nx = repmat(srcnorm(1,:),nt,1);
   ny = repmat(srcnorm(2,:),nt,1);
   submat = -(grad(:,:,1).*nx + grad(:,:,2).*ny);
-end
 
 % double layer (difference)
-if strcmpi(type,'d_diff')
+case {'d_diff', 'double_diff'}
   srcnorm = srcinfo.n(:,:);
   [~,grad] = chnk.helm2d.helmdiffgreen(zk,src,targ);
   nx = repmat(srcnorm(1,:),nt,1);
   ny = repmat(srcnorm(2,:),nt,1);
   submat = -(grad(:,:,1).*nx + grad(:,:,2).*ny);
-end
 
 % normal derivative of single layer
-if strcmpi(type,'sprime')
+case {'sp', 'sprime'}
   targnorm = targinfo.n(:,:);
   [~,grad] = chnk.helm2d.green(zk,src,targ);
   nx = repmat((targnorm(1,:)).',1,ns);
   ny = repmat((targnorm(2,:)).',1,ns);
 
   submat = (grad(:,:,1).*nx + grad(:,:,2).*ny);
-end
 
 % normal derivative of single layer (difference)
-if strcmpi(type,'sprime_diff')
+case {'sp_diff', 'sprime_diff'}
   targnorm = targinfo.n(:,:);
   [~,grad] = chnk.helm2d.helmdiffgreen(zk,src,targ);
   nx = repmat((targnorm(1,:)).',1,ns);
   ny = repmat((targnorm(2,:)).',1,ns);
 
   submat = (grad(:,:,1).*nx + grad(:,:,2).*ny);
-end
 
 
 % Tangential derivative of single layer
-if strcmpi(type,'stau')
+case {'stau', 'st'}
   targtan = targinfo.d(:,:);
   [~,grad] = chnk.helm2d.green(zk,src,targ);
   dx = repmat((targtan(1,:)).',1,ns);
   dy = repmat((targtan(2,:)).',1,ns);
   ds = sqrt(dx.*dx+dy.*dy);
   submat = (grad(:,:,1).*dx + grad(:,:,2).*dy)./ds;
-end
 
 % Tangential derivative of single layer (difference)
-if strcmpi(type,'stau_diff')
+case {'stau_diff','st_diff'}
   targtan = targinfo.d(:,:);
   [~,grad] = chnk.helm2d.helmdiffgreen(zk,src,targ);
   dx = repmat((targtan(1,:)).',1,ns);
   dy = repmat((targtan(2,:)).',1,ns);
   ds = sqrt(dx.*dx+dy.*dy);
   submat = (grad(:,:,1).*dx + grad(:,:,2).*dy)./ds;
-end
 
 % single layer
-if strcmpi(type,'s')
+case {'s', 'single'}
   submat = chnk.helm2d.green(zk,src,targ);
-end
 
 % single layer (difference)
-if strcmpi(type,'s_diff')
+case {'s_diff', 'single_diff'}
   submat = chnk.helm2d.helmdiffgreen(zk,src,targ);
-end
+
+% gradient of single layer
+case {'sgrad','sg'}
+    [~,grad] = chnk.helm2d.green(zk,src,targ);
+    submat = reshape(permute(grad,[3,1,2]),[],ns);
+
+% gradient of single layer(difference)
+case {'sgrad_diff','sg_diff'}
+    [~,grad] = chnk.helm2d.helmdiffgreen(zk,src,targ);
+    submat = reshape(permute(grad,[3,1,2]),[],ns);
 
 % normal derivative of double layer
-if strcmpi(type,'dprime')
+case {'dp', 'dprime'}
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
   [~,~,hess] = chnk.helm2d.green(zk,src,targ);
@@ -152,10 +155,9 @@
   nytarg = repmat((targnorm(2,:)).',1,ns);
   submat = -(hess(:,:,1).*nxsrc.*nxtarg + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)...
       + hess(:,:,3).*nysrc.*nytarg);
-end
 
 % normal derivative of double layer (difference)
-if strcmpi(type,'dprime_diff')
+case {'dp_diff', 'dprime_diff'}
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
   [~,~,hess] = chnk.helm2d.helmdiffgreen(zk,src,targ);
@@ -165,10 +167,21 @@
   nytarg = repmat((targnorm(2,:)).',1,ns);
   submat = -(hess(:,:,1).*nxsrc.*nxtarg + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)...
       + hess(:,:,3).*nysrc.*nytarg);
-end
+
+% gradient of double layer
+case {'dgrad','dg'}
+    [~,~,hess] = chnk.helm2d.green(zk,src,targ);
+    submat = -(hess(:,:,1:2).*srcinfo.n(1,:)+hess(:,:,2:3).*srcinfo.n(2,:));
+    submat = reshape(permute(submat,[3,1,2]),2*nt,ns);
+
+% gradient of double layer (difference)
+case {'dgrad_diff','dg_diff'}
+    [~,~,hess] = chnk.helm2d.helmdiffgreen(zk,src,targ);
+    submat = -(hess(:,:,1:2).*srcinfo.n(1,:)+hess(:,:,2:3).*srcinfo.n(2,:));
+    submat = reshape(permute(submat,[3,1,2]),2*nt,ns);
 
 % Combined field 
-if strcmpi(type,'c')
+case {'c', 'combined'}
   srcnorm = srcinfo.n(:,:);
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end
@@ -177,10 +190,9 @@
   ny = repmat(srcnorm(2,:),nt,1);
   submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny);
   submat = coef(1)*submatd + coef(2)*submats;
-end
 
 % Combined field (difference)
-if strcmpi(type,'c_diff')
+case {'c_diff', 'combined_diff'}
   srcnorm = srcinfo.n(:,:);
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end
@@ -189,16 +201,13 @@
   ny = repmat(srcnorm(2,:),nt,1);
   submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny);
   submat = coef(1)*submatd + coef(2)*submats;
-end
 
 % normal derivative of combined field
-if strcmpi(type,'cprime')
+case {'cp', 'cprime'}
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
-  
-
 
   % Get gradient and hessian info
   [~,grad,hess] = chnk.helm2d.green(zk,src,targ);
@@ -216,16 +225,13 @@
   submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg);
 
   submat = coef(1)*submatdp + coef(2)*submatsp;
-end
 
 % normal derivative of combined field (difference)
-if strcmpi(type,'cprime_diff')
+case {'cp_diff', 'cprime_diff'}
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
-  
-
 
   % Get gradient and hessian info
   [~,grad,hess] = chnk.helm2d.helmdiffgreen(zk,src,targ);
@@ -243,11 +249,29 @@
   submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg);
 
   submat = coef(1)*submatdp + coef(2)*submatsp;
-end
 
+% gradient of combined field
+case {'cgrad','cg'}
+    coef = ones(2,1);
+    if(nargin == 5); coef = varargin{1}; end
+    [~,grad,hess] = chnk.helm2d.green(zk,src,targ);
 
-% Dirichlet and neumann data corresponding to combined field
-if strcmpi(type,'c2trans') 
+    submat = -(hess(:,:,1:2).*srcinfo.n(1,:)+hess(:,:,2:3).*srcinfo.n(2,:));
+    submat = coef(1)*reshape(permute(submat,[3,1,2]),2*nt,ns);
+    submat = submat+coef(2)*reshape(permute(grad,[3,1,2]),2*nt,ns);
+
+% gradient of combined field (difference)
+case {'cgrad_diff','cg_diff'}
+    coef = ones(2,1);
+    if(nargin == 5); coef = varargin{1}; end
+    [~,grad,hess] = chnk.helm2d.helmdiffgreen(zk,src,targ);
+
+    submat = -(hess(:,:,1:2).*srcinfo.n(1,:)+hess(:,:,2:3).*srcinfo.n(2,:));
+    submat = coef(1)*reshape(permute(submat,[3,1,2]),2*nt,ns);
+    submat = submat+coef(2)*reshape(permute(grad,[3,1,2]),2*nt,ns);
+
+% Dirichlet and Neumann data corresponding to combined field
+case {'c2trans', 'c2t', 'c2tr'}
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end
   targnorm = targinfo.n(:,:);
@@ -275,11 +299,9 @@
   submat(1:2:2*nt,:) = coef(1)*submatd + coef(2)*submats;
   submat(2:2:2*nt,:) = coef(1)*submatdp + coef(2)*submatsp;
 
-end
-
 
-% Dirichlet and neumann data corresponding to combined field (difference)
-if strcmpi(type,'c2trans_diff') 
+% Dirichlet and Neumann data corresponding to combined field (difference)
+case {'c2trans_diff', 'c2t_diff', 'c2tr_diff'}
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end
   targnorm = targinfo.n(:,:);
@@ -307,10 +329,8 @@
   submat(1:2:2*nt,:) = coef(1)*submatd + coef(2)*submats;
   submat(2:2:2*nt,:) = coef(1)*submatdp + coef(2)*submatsp;
 
-end
-
 % all kernels, [c11 D, c12 S; c21 D', c22 S'] 
-if strcmpi(type,'all')
+case {'all','trans_sys','tsys'}
  
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
@@ -334,15 +354,13 @@
   % D
   submatd  = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc);
   
-  
   submat(1:2:2*nt,1:2:2*ns) = submatd*cc(1,1);
   submat(1:2:2*nt,2:2:2*ns) = submats*cc(1,2);
   submat(2:2:2*nt,1:2:2*ns) = submatdp*cc(2,1);
   submat(2:2:2*nt,2:2:2*ns) = submatsp*cc(2,2);
-end
 
 % all kernels, [c11 D, c12 S; c21 D', c22 S'] (difference)
-if strcmpi(type,'all_diff')
+case {'all_diff','trans_sys_diff','ts_diff'}
  
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
@@ -371,10 +389,9 @@
   submat(1:2:2*nt,2:2:2*ns) = submats*cc(1,2);
   submat(2:2:2*nt,1:2:2*ns) = submatdp*cc(2,1);
   submat(2:2:2*nt,2:2:2*ns) = submatsp*cc(2,2);
-end
 
 % Dirichlet data/potential correpsonding to transmission rep
-if strcmpi(type,'trans_rep') 
+case {'trans_rep','trep'} 
 
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end;
@@ -387,10 +404,9 @@
   submat = zeros(1*nt,2*ns);
   submat(1:1:1*nt,1:2:2*ns) = coef(1)*submatd;
   submat(1:1:1*nt,2:2:2*ns) = coef(2)*submats;
-end
 
 % Dirichlet data/potential correpsonding to transmission rep (difference)
-if strcmpi(type,'trans_rep_diff') 
+case {'trans_rep_diff','trep_diff'} 
 
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end;
@@ -403,10 +419,9 @@
   submat = zeros(1*nt,2*ns);
   submat(1:1:1*nt,1:2:2*ns) = coef(1)*submatd;
   submat(1:1:1*nt,2:2:2*ns) = coef(2)*submats;
-end
 
 % Neumann data corresponding to transmission rep
-if strcmpi(type,'trans_rep_prime')
+case {'trans_rep_prime','trep_p', 'trans_rep_p'}
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end;
   targnorm = targinfo.n(:,:);
@@ -431,10 +446,9 @@
   submat = zeros(nt,2*ns);
   submat(1:1:1*nt,1:2:2*ns) = coef(1)*submatdp;
   submat(1:1:1*nt,2:2:2*ns) = coef(2)*submatsp;
-end
 
 % Neumann data corresponding to transmission rep (difference)
-if strcmpi(type,'trans_rep_prime_diff')
+case {'trans_rep_prime_diff','trep_p_diff', 'trans_rep_p_diff'}
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end;
   targnorm = targinfo.n(:,:);
@@ -459,11 +473,9 @@
   submat = zeros(nt,2*ns);
   submat(1:1:1*nt,1:2:2*ns) = coef(1)*submatdp;
   submat(1:1:1*nt,2:2:2*ns) = coef(2)*submatsp;
-end
-
 
 % Gradient correpsonding to transmission rep
-if strcmpi(type,'trans_rep_grad')
+case {'trans_rep_grad','trep_g','trans_rep_g'}
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end;
   
@@ -485,12 +497,10 @@
     
   submat(2:2:2*nt,1:2:2*ns) = -coef(1)*(hess(:,:,2).*nxsrc + hess(:,:,3).*nysrc);
   submat(2:2:2*nt,2:2:2*ns) = coef(2)*grad(:,:,2);
-end
-
 
 
 % Gradient correpsonding to transmission rep (difference)
-if strcmpi(type,'trans_rep_grad_diff')
+case {'trans_rep_grad_diff','trep_g_diff','trans_rep_g_diff'}
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end;
   
@@ -512,4 +522,6 @@
     
   submat(2:2:2*nt,1:2:2*ns) = -coef(1)*(hess(:,:,2).*nxsrc + hess(:,:,3).*nysrc);
   submat(2:2:2*nt,2:2:2*ns) = coef(2)*grad(:,:,2);
+otherwise
+    error('Unknown Helmholtz kernel type ''%s''.', type);
 end
\ No newline at end of file
diff --git a/chunkie/+chnk/+helm2dquas/kern.m b/chunkie/+chnk/+helm2dquas/kern.m
index ecb45682..9cc58384 100644
--- a/chunkie/+chnk/+helm2dquas/kern.m
+++ b/chunkie/+chnk/+helm2dquas/kern.m
@@ -90,17 +90,17 @@
 [~,nt] = size(targ);
 nkappa = length(kappa);
 
+switch lower(type)
 % double layer
-if strcmpi(type,'d')
+case {'d', 'double'}
   srcnorm = srcinfo.n(:,:);
   [~,grad] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
   nx = repmat(srcnorm(1,:),nkappa*nt,1);
   ny = repmat(srcnorm(2,:),nkappa*nt,1);
   submat = -(grad(:,:,1).*nx + grad(:,:,2).*ny);
-end
 
 % normal derivative of single layer
-if strcmpi(type,'sprime')
+case {'sp', 'sprime'}
   targnorm = targinfo.n(:,:);
   [~,grad] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
   nx = repmat(reshape(targnorm(1,:),1,nt,1),nkappa,1,ns);
@@ -109,23 +109,14 @@
   ny = reshape(ny,[],ns);
 
   submat = (grad(:,:,1).*nx + grad(:,:,2).*ny);
-end
-
-% x derivative of single layer
-if strcmpi(type,'sx')
-  [~,grad] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
-
-  submat = (grad(:,:,1));
-end
 
 % single later
-if strcmpi(type,'s')
+case {'s', 'single'}
   submat = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
   submat = reshape(submat,[],ns);
-end
 
 % normal derivative of double layer
-if strcmpi(type,'dprime')
+case {'dp', 'dprime'}
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
   [~,~,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
@@ -139,24 +130,21 @@
   nytarg = reshape(nytarg,[],ns);
   submat = -(hess(:,:,1).*nxsrc.*nxtarg + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)...
       + hess(:,:,3).*nysrc.*nytarg);
-end
 
 
 % gradient of single layer
-if strcmpi(type,'sgrad')
+case {'sgrad','sg'}
     [~,grad] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
     submat = reshape(permute(grad,[3,1,2]),[],ns);
-end
 
 % gradient of double layer
-if strcmpi(type,'dgrad')
+case {'dgrad','dg'}
     [~,~,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
     submat = -(hess(:,:,1:2).*reshape(srcinfo.n(1,:),1,[],1)+hess(:,:,2:3).*reshape(srcinfo.n(2,:),1,[],1));
     submat = reshape(permute(submat,[3,1,2]),[],ns);
-end
 
 % Combined field 
-if strcmpi(type,'c')
+case {'c', 'combined'}
   srcnorm = srcinfo.n(:,:);
   coef = ones(2,1);
   if(nargin == 6); coef = varargin{1}; end
@@ -165,10 +153,9 @@
   ny = repmat(srcnorm(2,:),nkappa*nt,1);
   submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny);
   submat = coef(1)*submatd + coef(2)*reshape(submats,[],ns);
-end
 
 % normal derivative of combined field
-if strcmpi(type,'cprime')
+case {'cp', 'cprime'}
   coef = ones(2,1);
   if(nargin == 6); coef = varargin{1}; end
   targnorm = targinfo.n(:,:);
@@ -192,10 +179,9 @@
   submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg);
 
   submat = coef(1)*submatdp + coef(2)*submatsp;
-end
 
-% Dirichlet and neumann data corresponding to combined field
-if strcmpi(type,'c2trans') 
+% Dirichlet and Neumann data corresponding to combined field
+case {'c2trans', 'c2t', 'c2tr'}
   coef = ones(2,1);
   if(nargin == 6); coef = varargin{1}; end
   targnorm = targinfo.n(:,:);
@@ -225,10 +211,8 @@
   submat(1:2:end,:) = coef(1)*submatd + coef(2)*reshape(submats,[],ns);
   submat(2:2:end,:) = coef(1)*submatdp + coef(2)*submatsp;
 
-end
-
 % gradient of combined field
-if strcmpi(type,'cgrad')
+case {'cgrad','cg'}
     coef = ones(2,1);
     if(nargin == 6); coef = varargin{1}; end
     [~,grad,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
@@ -239,11 +223,10 @@
     submatd = reshape(permute(submatd,[3,1,2]),[],ns);
     
     submat = coef(1)*submatd + coef(2)*submats;
-end
   
 
 % all kernels, [c11 D, c12 S; c21 D', c22 S'] 
-if strcmpi(type,'all')
+case {'all','trans_sys','tsys'}
  
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
@@ -274,10 +257,9 @@
   submat(1:2:end,2:2:2*ns) = submats*cc(1,2);
   submat(2:2:end,1:2:2*ns) = submatdp*cc(2,1);
   submat(2:2:end,2:2:2*ns) = submatsp*cc(2,2);
-end
-% Dirichlet data/potential correpsonding to transmission rep
-if strcmpi(type,'trans_rep') 
 
+% Dirichlet data/potential correpsonding to transmission rep
+case {'trans_rep','trep'} 
   coef = ones(2,1);
   if(nargin == 6); coef = varargin{1}; end;
   srcnorm = srcinfo.n(:,:);
@@ -291,10 +273,9 @@
   submat = zeros(nkappa*nt,2*ns);
   submat(1:1:end,1:2:2*ns) = coef(1)*submatd;
   submat(1:1:end,2:2:2*ns) = coef(2)*submats;
-end
 
 % Neumann data corresponding to transmission rep
-if strcmpi(type,'trans_rep_prime')
+case {'trans_rep_prime','trep_p', 'trans_rep_p'}
   coef = ones(2,1);
   if(nargin == 6); coef = varargin{1}; end;
   targnorm = targinfo.n(:,:);
@@ -319,10 +300,9 @@
   submat = zeros(nkappa*nt,2*ns);
   submat(1:1:end,1:2:2*ns) = coef(1)*submatdp;
   submat(1:1:end,2:2:2*ns) = coef(2)*submatsp;
-end
 
 % Gradient correpsonding to transmission rep
-if strcmpi(type,'trans_rep_grad')
+case {'trans_rep_grad','trep_g', 'trans_rep_g'}
   coef = ones(2,1);
   if(nargin == 6); coef = varargin{1}; end;
   
@@ -344,11 +324,9 @@
     
   submat(2:2:end,1:2:2*ns) = -coef(1)*(hess(:,:,2).*nxsrc + hess(:,:,3).*nysrc);
   submat(2:2:end,2:2:2*ns) = coef(2)*grad(:,:,2);
+otherwise
+    error('Unknown quasi-periodic Helmholtz kernel type ''%s''.', type);
 end
 
 
-
-
-
-
 end
\ No newline at end of file
diff --git a/chunkie/+chnk/+lap2d/kern.m b/chunkie/+chnk/+lap2d/kern.m
index 9a76c896..248e32f9 100644
--- a/chunkie/+chnk/+lap2d/kern.m
+++ b/chunkie/+chnk/+lap2d/kern.m
@@ -9,82 +9,151 @@
 [~,ns] = size(src);
 [~,nt] = size(targ);
 
-if strcmpi(type,'d')
+switch lower(type)
+% double layer
+case {'d', 'double'}
     %srcnorm = chnk.normal2d(srcinfo);
     [~,grad] = chnk.lap2d.green(src,targ,true);
     submat = -(grad(:,:,1).*srcinfo.n(1,:) + grad(:,:,2).*srcinfo.n(2,:));
-end
 
-if strcmpi(type,'sprime')
-    targnorm = chnk.normal2d(targinfo);
+% normal derivative of single layer
+case {'sp', 'sprime'}
+    if ~isfield(targinfo,'n')
+        targnorm = chnk.normal2d(targinfo);
+    else
+        targnorm = targinfo.n;
+    end
     [~,grad] = chnk.lap2d.green(src,targ,true);
     nx = repmat((targnorm(1,:)).',1,ns);
     ny = repmat((targnorm(2,:)).',1,ns);
 
     submat = (grad(:,:,1).*nx + grad(:,:,2).*ny);
-end
 
-if strcmpi(type,'stau')
-    targnorm = chnk.normal2d(targinfo);
+% Tangential derivative of single layer
+case {'stau'}
+    if ~isfield(targinfo,'n')
+        targnorm = chnk.normal2d(targinfo);
+    else
+        targnorm = targinfo.n;
+    end
     [~,grad] = chnk.lap2d.green(src,targ,true);
     nx = repmat((targnorm(1,:)).',1,ns);
     ny = repmat((targnorm(2,:)).',1,ns);
 
     submat = (-grad(:,:,1).*ny + grad(:,:,2).*nx);
-end
 
-if strcmpi(type,'hilb') % hilbert transform (two times the adjoint of stau)
-    srcnorm = chnk.normal2d(srcinfo);
+% Hilbert transform (two times the adjoint of stau)
+case {'hilb'} 
+    if ~isfield(srcinfo,'n')
+        srcnorm = chnk.normal2d(srcinfo);
+    else
+        srcnorm = srcinfo.n;
+    end
     [~,grad] = chnk.lap2d.green(src,targ,true);
     nx = repmat((srcnorm(1,:)),nt,1);
     ny = repmat((srcnorm(2,:)),nt,1);
 
     submat = 2*(grad(:,:,1).*ny - grad(:,:,2).*nx);
-end
 
-if strcmpi(type,'sgrad')
+% gradient of single layer
+case {'sgrad','sg'}
     [~,grad] = chnk.lap2d.green(src,targ,true);
     submat = reshape(permute(grad,[3,1,2]),2*nt,ns);
-end
 
-if strcmpi(type,'dgrad')
+% gradient of double layer
+case {'dgrad','dg'}
+    if ~isfield(srcinfo,'n')
+        srcnorm = chnk.normal2d(srcinfo);
+    else
+        srcnorm = srcinfo.n;
+    end
     [~,~,hess] = chnk.lap2d.green(src,targ,true);
-    submat = -(hess(:,:,1:2).*srcinfo.n(1,:)+hess(:,:,2:3).*srcinfo.n(2,:));
+    submat = -(hess(:,:,1:2).*srcnorm(1,:)+hess(:,:,2:3).*srcnorm(2,:));
     submat = reshape(permute(submat,[3,1,2]),2*nt,ns);
-end
 
-if strcmpi(type,'dprime')
-  targnorm = targinfo.n(:,:);
-  srcnorm = srcinfo.n(:,:);
-  [~,~,hess] = chnk.lap2d.green(src,targ);
-  nxsrc = repmat(srcnorm(1,:),nt,1);
-  nysrc = repmat(srcnorm(2,:),nt,1);
-  nxtarg = repmat((targnorm(1,:)).',1,ns);
-  nytarg = repmat((targnorm(2,:)).',1,ns);
-  submat = -(hess(:,:,1).*nxsrc.*nxtarg + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)...
+% normal derivative of double layer
+case {'dprime','dp'}
+    if ~isfield(targinfo,'n')
+        targnorm = chnk.normal2d(targinfo);
+    else
+        targnorm = targinfo.n;
+    end
+    if ~isfield(srcinfo,'n')
+        srcnorm = chnk.normal2d(srcinfo);
+    else
+        srcnorm = srcinfo.n;
+    end
+    [~,~,hess] = chnk.lap2d.green(src,targ);
+    nxsrc = repmat(srcnorm(1,:),nt,1);
+    nysrc = repmat(srcnorm(2,:),nt,1);
+    nxtarg = repmat((targnorm(1,:)).',1,ns);
+    nytarg = repmat((targnorm(2,:)).',1,ns);
+    submat = -(hess(:,:,1).*nxsrc.*nxtarg + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)...
       + hess(:,:,3).*nysrc.*nytarg);
-end
-
 
-if strcmpi(type,'s')
+% single layer
+case {'s', 'single'}
     submat = chnk.lap2d.green(src,targ);
-end
 
-if strcmpi(type,'c')
-    srcnorm = chnk.normal2d(srcinfo);
+% combined field
+case {'c', 'combined'}
+    if ~isfield(srcinfo,'n')
+        srcnorm = chnk.normal2d(srcinfo);
+    else
+        srcnorm = srcinfo.n;
+    end
     coef = ones(2,1);
     if(nargin == 4); coef = varargin{1}; end
     [s,grad] = chnk.lap2d.green(src,targ);
     nx = repmat(srcnorm(1,:),nt,1);
     ny = repmat(srcnorm(2,:),nt,1);
     submat = -coef(1)*(grad(:,:,1).*nx + grad(:,:,2).*ny) + coef(2)*s;
-end
 
-if strcmpi(type,'cgrad')
+% normal derivative of combined field
+case {'cprime','cp'}
+    if ~isfield(targinfo,'n')
+        targnorm = chnk.normal2d(targinfo);
+    else
+        targnorm = targinfo.n;
+    end
+    if ~isfield(srcinfo,'n')
+        srcnorm = chnk.normal2d(srcinfo);
+    else
+        srcnorm = srcinfo.n;
+    end
+    coef = ones(2,1);
+    if(nargin == 4); coef = varargin{1}; end
+    [~,grad,hess] = chnk.lap2d.green(src,targ);
+    nxsrc = repmat(srcnorm(1,:),nt,1);
+    nysrc = repmat(srcnorm(2,:),nt,1);
+    nxtarg = repmat((targnorm(1,:)).',1,ns);
+    nytarg = repmat((targnorm(2,:)).',1,ns);
+    submatdp = -(hess(:,:,1).*nxsrc.*nxtarg + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)...
+      + hess(:,:,3).*nysrc.*nytarg);
+
+    % S'
+    submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg);
+    
+    submat = coef(1)*submatdp + coef(2)*submatsp;
+
+
+% gradient of combined field
+case {'cgrad', 'cg'}
+    if ~isfield(srcinfo,'n')
+        srcnorm = chnk.normal2d(srcinfo);
+    else
+        srcnorm = srcinfo.n;
+    end
     coef = ones(2,1);
     if(nargin == 4); coef = varargin{1}; end
     [~,grad,hess] = chnk.lap2d.green(src,targ,true);
-    submat = -(hess(:,:,1:2).*srcinfo.n(1,:)+hess(:,:,2:3).*srcinfo.n(2,:));
+    submat = -(hess(:,:,1:2).*srcnorm(1,:)+hess(:,:,2:3).*srcnorm(2,:));
     submat = coef(1)*reshape(permute(submat,[3,1,2]),2*nt,ns);
     submat = submat+coef(2)*reshape(permute(grad,[3,1,2]),2*nt,ns);
+
+otherwise
+    error('Unknown Laplace kernel type ''%s''.', type);
 end
+
+end
+
diff --git a/chunkie/+chnk/+stok2d/kern.m b/chunkie/+chnk/+stok2d/kern.m
index 2f7d979a..d2ed1f3e 100644
--- a/chunkie/+chnk/+stok2d/kern.m
+++ b/chunkie/+chnk/+stok2d/kern.m
@@ -8,7 +8,7 @@
 
 switch lower(type)
 
-    case {'svel', 's'}
+    case {'svel', 's', 'single'}
 
         r = sqrt(r2);
         log1r = log(1./r);
@@ -65,7 +65,7 @@
             varargout = {Kxx, Kxy, Kyy};
         end
     
-    case 'sgrad'
+    case {'sgrad', 'sg'}
         r4 = r2.^2;
         r4inv = 1./r4./(4*pi*mu);
         Ku1x1 = rx.*(ry.^2 - rx.^2).*r4inv;
@@ -100,7 +100,7 @@
 
 
 
-    case {'dvel', 'd'}
+    case {'dvel', 'd', 'double'}
 
         r4 = r2.^2;
         nx = srcinfo.n(1,:);
@@ -182,7 +182,8 @@
         else
             varargout = {Kxx, Kxy, Kyy};
         end
-    case 'dgrad'
+
+    case {'dgrad', 'dg'}
         r4 = r2.^2;
         r6 = r2.^3;
         nx = srcinfo.n(1,:);
@@ -224,7 +225,7 @@
         end
 
 
-    case 'c'
+    case {'c', 'combined'}
 
         coefs = varargin{1};
         [Sxx, Sxy, Syy] = chnk.stok2d.kern(mu, srcinfo, targinfo, 's');
@@ -244,6 +245,7 @@
         else
             varargout = {Kxx, Kxy, Kyy};
         end
+
     case 'cpres'
 
         coefs = varargin{1};
@@ -282,7 +284,8 @@
         else
             varargout = {Kxx, Kxy, Kyy};
         end
-    case 'cgrad'
+
+    case {'cgrad', 'cg'}
         coefs = varargin{1};
         [Su1x1, Su1x2, Su1y1, Su1y2, Su2x1, Su2x2, Su2y1, Su2y2] = ...
            chnk.stok2d.kern(mu, srcinfo, targinfo, 'sgrad');
diff --git a/chunkie/@kernel/helm2d.m b/chunkie/@kernel/helm2d.m
index cc12d025..0675fa0b 100644
--- a/chunkie/@kernel/helm2d.m
+++ b/chunkie/@kernel/helm2d.m
@@ -18,6 +18,29 @@
 %   constructs the derivative of the combined-layer Helmholtz kernel with 
 %   wavenumber ZK and  parameter ETA, i.e., COEFS(1)*KERNEL.HELM2D('dp', ZK) + 
 %   COEFS(2)*KERNEL.HELM2D('sp', ZK).
+%
+%   KERNEL.HELM2D('all', ZK, COEFS) or KERNEL.HELM2D('tsys', ZK, COEFS) 
+%   constructs the (2x2) matrix of kernels [D, S; D' S'] scaled by coefs
+%   where coefs is (2x2) matrix, i.e. D part is scaled as
+%   COEFS(1,1)*KERNEL.HELM2D('d', zk), and so on.
+%
+%   KERNEL.HELM2D('trans_rep', ZK, COEFS) or KERNEL.HELM2D('trep', ZK, COEFS) 
+%   constructs the transmission repretsentation, i.e. the (1x2) matrix of
+%   kernels [D, S] scaled by coefs where coefs is (1x2) matrix, i.e. D part
+%   is scaled as COEFS(1)*KERNEL.HELM2D('d', zk), and so on.
+%
+%   KERNEL.HELM2D('trans_rep_p', ZK, COEFS) or KERNEL.HELM2D('trep_p', ZK, COEFS) 
+%   constructs the derivative of the transmission repretsentation, i.e. the
+%   (1x2) matrix of kernels [D', S'] scaled by coefs where coefs is (1x2)
+%   matrix, i.e. D part is scaled as COEFS(1)*KERNEL.HELM2D('dp', zk), and
+%   so on.
+%
+%   KERNEL.HELM2D('c2trans', ZK, COEFS) or KERNEL.HELM2D('c2t', ZK, COEFS) 
+%   evaluates the combined-layer Helmholtz kernel and its derivative, i.e.
+%   the (2x1) matrix of kernels [C; C'] scaled by coefs where coefs is
+%   (2x2) matrix, i.e. kernel returns 
+%   [COEFS(1,1)*KERNEL.HELM2D('d',zk)+COEFS(1,2)*KERNEL.HELM2D('s', zk); 
+%   COEFS(2,1)*KERNEL.HELM2D('dp',zk)+COEFS(2,2)*KERNEL.HELM2D('sp', zk)]
 % See also CHNK.HELM2D.KERN.
 
 % author: Dan Fortunato
@@ -95,6 +118,50 @@
         obj.fmm  = @(eps,s,t,sigma) chnk.helm2d.fmm(eps, zk, s, t, 'cprime', sigma, coefs);
         obj.sing = 'hs';
 
+    case {'all', 'trans_sys', 'tsys'}
+        if ( nargin < 3 )
+            warning('Missing transmission coefficients. Defaulting to [1,1;1,1].');
+            coefs = ones(2,2);
+        end
+        obj.type = 'all';
+        obj.opdims = [2,2];
+        obj.eval = @(s,t) chnk.helm2d.kern(zk, s, t, 'all', coefs);
+        obj.fmm  = []; 
+        obj.sing = 'hs'; 
+
+    case {'trans_rep','trep'} 
+        if ( nargin < 3 )
+            warning('Missing transmission coefficients. Defaulting to [1;1].');
+            coefs = ones(2,1);
+        end
+        obj.type = 'trep';
+        obj.opdims = [1,2];
+        obj.eval = @(s,t) chnk.helm2d.kern(zk, s, t, 'trans_rep', coefs);
+        obj.fmm  = []; 
+        obj.sing = 'pv'; 
+    
+    case {'trans_rep_prime','trep_p', 'trans_rep_p'}
+        if ( nargin < 3 )
+            warning('Missing transmission coefficients. Defaulting to [1;1].');
+            coefs = ones(2,1);
+        end
+        obj.type = 'trep_p';
+        obj.opdims = [1,2];
+        obj.eval = @(s,t) chnk.helm2d.kern(zk, s, t, 'trans_rep_prime', coefs);
+        obj.fmm  = []; 
+        obj.sing = 'hs'; 
+
+    case {'c2trans', 'c2t', 'c2tr'}
+        if ( nargin < 3 )
+            warning('Missing combined layer coefficients. Defaulting to [1,1i;1,1i].');
+            coefs = [1,1i;1,1i];
+        end
+        obj.type = 'c2tr';
+        obj.opdims = [2,1];
+        obj.eval = @(s,t) chnk.helm2d.kern(zk, s, t, 'c2trans', coefs);
+        obj.fmm  = []; 
+        obj.sing = 'hs'; 
+       
     otherwise
         error('Unknown Helmholtz kernel type ''%s''.', type);
 
diff --git a/chunkie/@kernel/helm2ddiff.m b/chunkie/@kernel/helm2ddiff.m
index b15861f8..72a3bd0d 100644
--- a/chunkie/@kernel/helm2ddiff.m
+++ b/chunkie/@kernel/helm2ddiff.m
@@ -12,6 +12,12 @@
 %   COEFS(1)*KERNEL.HELM2D('d', ZKS(1)) -
 %   COEFS(2)*KERNEL.HELM2D('d', ZKS(2))
 %
+%   KERNEL.HELM2DDIFF('c', ZKS) or KERNEL.HELM2DDIFF('combined', ZKS) 
+%   constructs the difference of combined-layer Helmholtz kernels with 
+%   wavenumbers ZKS(1) and ZKS(2), scaled by COEFS, i.e.
+%   KERNEL.HELM2D('c', ZKS(1), COEFS(:,1)) -
+%   KERNEL.HELM2D('c', ZKS(2), COEFS(:,2))
+%
 %   KERNEL.HELM2DDIFF('sp', ZKS) or KERNEL.HELM2DDIFF('sprime', ZKS) 
 %   constructs the difference of Neumann data for single-layer Helmholtz
 %   kernels with wavenumbers ZKS(1) and ZKS(2), scaled by COEFS, i.e.
@@ -24,12 +30,38 @@
 %   COEFS(1)*KERNEL.HELM2D('dp', ZKS(1)) -
 %   COEFS(2)*KERNEL.HELM2D('dp', ZKS(2))
 %
-%   KERNEL.HELM2DDIFF('all', ZKS, coefs) 
+%   KERNEL.HELM2DDIFF('cp', ZKS) or KERNEL.HELM2DDIFF('cprime', ZKS) 
+%   constructs the difference of Neumann for combined-layer Helmholtz 
+%   kernels with wavenumbers ZKS(1) and ZKS(2), scaled by COEFS, i.e.
+%   KERNEL.HELM2D('cp', ZKS(1), COEFS(:,1)) -
+%   KERNEL.HELM2D('cp', ZKS(2), COEFS(:,2))
+%
+%   KERNEL.HELM2DDIFF('all', ZKS, coefs) or 
+%   KERNEL.HELM2DDIFF('tsys', ZKS, coefs)
 %   constructs the (2x2) matrix of difference kernels of
 %   [D, S; D' S'] scaled by coefs where coefs is (2x2x2) tensor, i.e.
 %   D part is scaled as COEFS(1,1,1)*KERNEL.HELM2D('d', zks(1)) - 
 %   COEFS(1,1,2)*KERNEL.HELM2D('d', ZKS(2)), and so on.
 %
+%   KERNEL.HELM2DDIFF('trans_rep', ZK, COEFS) or KERNEL.HELM2DDIFF('trep', ZK, COEFS) 
+%   constructs the difference of transmission repretsentations, i.e. the (1x2) matrix of
+%   kernels [D, S] scaled by coefs where coefs is (1x2) matrix, i.e. D part
+%   is scaled as COEFS(1,1)*KERNEL.HELM2DDIFF('d', ZKS(1)) - 
+%   COEFS(1,2)*KERNEL.HELM2DDIFF('d', ZKS(2)), and so on.
+%
+%   KERNEL.HELM2DDIFF('trans_rep_p', ZK, COEFS) or KERNEL.HELM2DDIFF('trep_p', ZK, COEFS) 
+%   constructs the difference of derivative of the transmission repretsentation, i.e. the
+%   (1x2) matrix of kernels [D', S'] scaled by coefs where coefs is (1x2)
+%   matrix, i.e. D part is scaled as COEFS(1,1)*KERNEL.HELM2DDIFF('dp', ZKS(1)) - 
+%   COEFS(1,2)*KERNEL.HELM2DDIFF('dp', ZKS(2)), and so on.
+%
+%   KERNEL.HELM2DDIFF('c2trans', ZK, COEFS) or KERNEL.HELM2DDIFF('c2t', ZK, COEFS) 
+%   evaluates the difference of combined-layer Helmholtz kernel and its derivative, i.e.
+%   the (2x1) matrix of kernels [C; C'] scaled by coefs where coefs is
+%   (2x2) matrix, i.e. kernel returns 
+%   KERNEL.HELM2DDIFF('c2trans',zks(1),coefs(:,:,1)) - 
+%   KERNEL.HELM2DDIFF('c2trans',zks(2),coefs(:,:,2)) - 
+%
 %  COEFS is optional. If not given then COEFS is set to all ones.
 %
 % See also CHNK.HELM2D.KERN.
@@ -48,7 +80,7 @@
 obj.opdims = [1 1];
 if(nargin < 3)
     coefs = [1 1];
-    if strcmpi(type, 'all')
+    if strcmpi(type, 'all') || strcmpi(type, 'trans_sys') || strcmpi(type, 'tsys')
       coefs = ones(2,2,2);
     end
 end
@@ -74,6 +106,18 @@
             coefs(2)*chnk.helm2d.fmm(eps, zks(2), s, t, 'd', sigma);
         obj.sing = 'log';
 
+    case {'c', 'combined'}
+        if (nargin < 3)
+            coefs = ones(2,2);
+        end
+        obj.type = 'c';
+        obj.eval = @(s,t) chnk.helm2d.kern(zks(1), s, t, 'c_diff',coefs(:,1)) - ...
+                          chnk.helm2d.kern(zks(2), s, t, 'c_diff',coefs(:,2));
+        obj.fmm  = @(eps,s,t,sigma) ...
+            chnk.helm2d.fmm(eps, zks(1), s, t, 'c', sigma,coefs(:,1)) - ...
+            chnk.helm2d.fmm(eps, zks(2), s, t, 'c', sigma,coefs(:,2));
+        obj.sing = 'log';
+
     case {'sp', 'sprime'}
         obj.type = 'sp';
         obj.eval = @(s,t) coefs(1)*chnk.helm2d.kern(zks(1), s, t, 'sprime_diff') - ...
@@ -96,7 +140,22 @@
             obj.sing = 'hs';
         end
 
-    case {'all'}
+    case {'cp', 'cprime'}
+        if (nargin < 3)
+            coefs = ones(2,2);
+        end
+        obj.type = 'c';
+        obj.eval = @(s,t) chnk.helm2d.kern(zks(1), s, t, 'cp_diff',coefs(:,1)) - ...
+                          chnk.helm2d.kern(zks(2), s, t, 'cp_diff',coefs(:,2));
+        obj.fmm  = @(eps,s,t,sigma) ...
+            chnk.helm2d.fmm(eps, zks(1), s, t, 'cp', sigma,coefs(:,1)) - ...
+            chnk.helm2d.fmm(eps, zks(2), s, t, 'cp', sigma,coefs(:,2));
+        obj.sing = 'log';
+
+    case {'all', 'trans_sys', 'tsys'}
+        if (nargin < 3)
+            coefs = ones(2,2,2);
+        end
         obj.type = 'all';
         obj.opdims = [2,2];
         obj.eval = @(s,t) chnk.helm2d.kern(zks(1), s, t, 'all_diff', coefs(:,:,1)) - ...
@@ -107,6 +166,51 @@
         else
             obj.sing = 'hs'; 
         end
+        
+    case {'trans_rep','trep'} 
+        if ( nargin < 3 )
+            coefs = ones(2,2);
+        end
+        obj.type = 'trep';
+        obj.opdims = [1,2];
+        obj.eval = @(s,t) chnk.helm2d.kern(zks(1), s, t, 'trep_diff', coefs(:,1)) - ...
+                          chnk.helm2d.kern(zks(2), s, t, 'trep_diff', coefs(:,2));
+        obj.fmm  = []; 
+        if ( abs(coefs(2,1)-coefs(2,1)) < eps )
+            obj.sing = 'log';
+        else
+            obj.sing = 'pv'; 
+        end
+    
+    case {'trans_rep_prime','trep_p', 'trans_rep_p'}
+        if ( nargin < 3 )
+            coefs = ones(2,2);
+        end
+        obj.type = 'trep_p';
+        obj.opdims = [1,2];
+        obj.eval = @(s,t) chnk.helm2d.kern(zks(1), s, t, 'trep_p_diff', coefs(:,1)) - ...
+                          chnk.helm2d.kern(zks(2), s, t, 'trep_p_diff', coefs(:,2));
+        obj.fmm  = []; 
+        if ( abs(coefs(2,1)-coefs(2,2)) < eps )
+            obj.sing = 'log';
+        else
+            obj.sing = 'hs'; 
+        end
+    
+    case {'c2trans', 'c2t', 'c2tr'}
+        if ( nargin < 3 )
+            coefs =  ones(2,2,2);
+        end
+        obj.type = 'c2tr';
+        obj.opdims = [2,1];
+        obj.eval = @(s,t) chnk.helm2d.kern(zks(1), s, t, 'c2tr_diff', coefs(:,:,1)) - ...
+                          chnk.helm2d.kern(zks(2), s, t, 'c2tr_diff', coefs(:,:,2));
+        obj.fmm  = []; 
+        if ( abs(coefs(2,1,1)-coefs(2,1,2)) < eps )
+            obj.sing = 'log';
+        else
+            obj.sing = 'hs'; 
+        end
 
     otherwise
         error('Unknown Helmholtz difference kernel type ''%s''.', type);
diff --git a/chunkie/@kernel/helm2dquas.m b/chunkie/@kernel/helm2dquas.m
index d71db420..879c5e59 100644
--- a/chunkie/@kernel/helm2dquas.m
+++ b/chunkie/@kernel/helm2dquas.m
@@ -25,6 +25,33 @@
 %   parameter ETA, i.e., COEFS(1)*KERNEL.HELM2DQUAS('dp', ZK, KAPPA, D) + 
 %      COEFS(2)*KERNEL.HELM2DQUAS('sp', ZK, KAPPA, D).
 %
+%   KERNEL.HELM2DQUAS('all', ZK, KAPPA, D, COEFS) or 
+%   KERNEL.HELM2DQUAS('tsys', ZK, KAPPA, D, COEFS)
+%   constructs the (2x2) matrix of kernels [D, S; D' S'] scaled by coefs
+%   where coefs is (2x2) matrix, i.e. D part is scaled as
+%   COEFS(1,1)*KERNEL.HELM2DQUAS('d', zk), and so on.
+%
+%   KERNEL.HELM2DQUAS('trans_rep', ZK, KAPPA, D, COEFS) or 
+%   KERNEL.HELM2DQUAS('trep', ZK, KAPPA, D, COEFS) 
+%   constructs the transmission repretsentation, i.e. the (1x2) matrix of
+%   kernels [D, S] scaled by coefs where coefs is (1x2) matrix, i.e. D part
+%   is scaled as COEFS(1)*KERNEL.HELM2D('d', zk), and so on.
+%
+%   KERNEL.HELM2DQUAS('trans_rep_p', ZK, KAPPA, D, COEFS) or 
+%   KERNEL.HELM2DQUAS('trep_p', ZK, KAPPA, D, COEFS) 
+%   constructs the derivative of the transmission repretsentation, i.e. the
+%   (1x2) matrix of kernels [D', S'] scaled by coefs where coefs is (1x2)
+%   matrix, i.e. D part is scaled as 
+%   COEFS(1)*KERNEL.HELM2DQUAS('dp', zk, KAPPA, D), and so on.
+%
+%   KERNEL.HELM2DQUAS('c2trans', ZK, KAPPA, D, COEFS) or 
+%   KERNEL.HELM2DQUAS('c2t', ZK, KAPPA, D, COEFS) 
+%   evaluates the combined-layer Helmholtz kernel and its derivative, i.e.
+%   the (2x1) matrix of kernels [C; C'] scaled by coefs where coefs is
+%   (2x2) matrix, i.e. kernel returns 
+%   [COEFS(1,1)*KERNEL.HELM2DQUAS('d',zk, KAPPA, D)+COEFS(1,2)*KERNEL.HELM2DQUAS('s', zk, KAPPA, D); 
+%   COEFS(2,1)*KERNEL.HELM2DQUAS('dp',zk, KAPPA, D)+COEFS(2,2)*KERNEL.HELM2DQUAS('sp', zk, KAPPA, D)]
+%
 %   if kappa is an array of values of length nkappa then the kernel 
 %   becomes an (nkappa x 1) vector-valued kernel containing the scalar values
 %   for each kappa. this is much more efficient than separate kernel 
@@ -156,6 +183,50 @@
         obj.fmm = [];
         obj.sing = 'hs';
 
+    case {'all', 'trans_sys', 'tsys'}
+        if ( nargin < 3 )
+            warning('Missing transmission coefficients. Defaulting to [1,1;1,1].');
+            coefs = ones(2,2);
+        end
+        obj.type = 'all';
+        obj.opdims = [2,2];
+        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'all',quas_param, coefs);
+        obj.fmm  = []; 
+        obj.sing = 'hs'; 
+    
+    case {'trans_rep','trep'} 
+        if ( nargin < 3 )
+            warning('Missing transmission coefficients. Defaulting to [1;1].');
+            coefs = ones(2,1);
+        end
+        obj.type = 'trep';
+        obj.opdims = [1,2];
+        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'trep',quas_param, coefs);
+        obj.fmm  = []; 
+        obj.sing = 'pv'; 
+    
+    case {'trans_rep_prime','trep_p', 'trans_rep_p'}
+        if ( nargin < 3 )
+            warning('Missing transmission coefficients. Defaulting to [1;1].');
+            coefs = ones(2,1);
+        end
+        obj.type = 'trep_p';
+        obj.opdims = [1,2];
+        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'trep_p',quas_param, coefs);
+        obj.fmm  = []; 
+        obj.sing = 'hs'; 
+
+    case {'c2trans', 'c2t', 'c2tr'}
+        if ( nargin < 3 )
+            warning('Missing combined layer coefficients. Defaulting to [1,1i;1,1i].');
+            coefs =  [1,1i;1,1i];
+        end
+        obj.type = 'c2tr';
+        obj.opdims = [2,1];
+        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'c2t',quas_param, coefs);
+        obj.fmm  = []; 
+        obj.sing = 'hs'; 
+
     otherwise
         error('Unknown Helmholtz kernel type ''%s''.', type);
 
diff --git a/chunkie/@kernel/kernel.m b/chunkie/@kernel/kernel.m
index 482f6ed2..a5c9e2b9 100644
--- a/chunkie/@kernel/kernel.m
+++ b/chunkie/@kernel/kernel.m
@@ -10,7 +10,8 @@
 %      'laplace'    ('lap', 'l')                    's', 'd', 'sp', 'c'
 %      'helmholtz'  ('helm', 'h')                   's', 'd', 'sp', 'dp', 'c'
 %                                                   'cp'
-%      'helmholtz difference' ('helmdiff', 'hdiff') 's', 'd', 'sp', 'dp'
+%      'helmholtz difference'                       's', 'd', 'sp', 'dp'
+%         ('helmdiff', 'hdiff', 'helm_diff')
 %      'elasticity' ('elast', 'e')                  's', 'strac', 'd', 'dalt'
 %      'stokes'     ('stok', 's')                   'svel', 'spres', 
 %                                                   'strac', 'sgrad'
@@ -22,7 +23,7 @@
 %      'axis sym helmholtz'                         's' 'd' 'sp' 'c'
 %         ('axissymh', 'axissymhelm')
 %      'axis sym helmholtz difference'              's' 'd' 'sp' 'dp'
-%         ('axissymhdiff', 'axissymhelmdiff') 
+%         ('axissymhdiff', 'axissymhelmdiff', 'axissymhelm_diff') 
 %      'quasiperiodic helmholtz'                    's', 'd', 'sp', 'dp', 'c'
 %         ('helmquas', 'hq')                        'cp'
 %   The types may also be written in longer form, e.g. 'single', 'double',
@@ -99,7 +100,7 @@
                       obj = kernel.lap2d(varargin{:});
                   case {'helmholtz', 'helm', 'h'}
                       obj = kernel.helm2d(varargin{:});
-                  case {'helmholtz difference', 'helmdiff', 'hdiff'}
+                  case {'helmholtz difference', 'helmdiff', 'hdiff', 'helm_diff'}
                       obj = kernel.helm2ddiff(varargin{:});
                   case {'stokes', 'stok', 's'}
                       obj = kernel.stok2d(varargin{:});
@@ -112,7 +113,7 @@
                   case {'axis sym helmholtz', 'axissymh', 'axissymhelm'}
                       obj = kernel.axissymhelm2d(varargin{:});
                   case {'axis sym helmholtz difference', 'axissymhdiff' ...
-                           'axissymhelmdiff'}
+                           'axissymhelmdiff', 'axissymhelm_diff'}
                       obj = kernel.axissymhelm2ddiff(varargin{:});   
                   case {'quasiperiodic helmholtz', 'helmquas', 'hq'}
                       obj = kernel.helm2dquas(varargin{:});   
diff --git a/chunkie/@kernel/lap2d.m b/chunkie/@kernel/lap2d.m
index c4a26f0c..62a95df3 100644
--- a/chunkie/@kernel/lap2d.m
+++ b/chunkie/@kernel/lap2d.m
@@ -23,6 +23,13 @@
 %   coefs(1)*KERNEL.LAP2D('d') + coefs(2)*KERNEL.LAP2D('s'). If no
 %   value of coefs is specified the default is coefs = [1 1]  
 %
+%   KERNEL.LAP2D('cp', coefs) or KERNEL.LAP2D('cprime', coefs) constructs
+%   the normal derivative of the combined-layer Laplace kernel with
+%   parameter coefs
+%
+%   KERNEL.LAP2D('cg', coefs) or KERNEL.LAP2D('cgrad', coefs) constructs
+%   the gradient of the combined-layer Laplace kernel with parameter coefs
+%
 % See also CHNK.LAP2D.KERN.
 
 % author: Dan Fortunato
@@ -104,6 +111,16 @@
         obj.splitinfo.action = {'r','r'};
         obj.splitinfo.functions = @(s,t) lap2d_c_split(s,t,coefs);
 
+    case {'cp', 'cprime'}
+        if ( nargin < 2 )
+            warning('Missing combined layer parameter coefs. Defaulting to [1 1].');
+            coefs = ones(2,1);
+        end
+        obj.type = 'cp';
+        obj.eval = @(s,t) chnk.lap2d.kern(s, t, 'cprime', coefs);
+        obj.fmm  = @(eps,s,t,sigma) chnk.lap2d.fmm(eps, s, t, 'cprime', sigma, coefs);
+        obj.sing = 'hs';
+
     case {'cg', 'cgrad'}
         if ( nargin < 2 )
             warning('Missing combined layer parameter coefs. Defaulting to [1 1].');
diff --git a/devtools/test/KernDerInterleaveTest.m b/devtools/test/KernDerInterleaveTest.m
new file mode 100644
index 00000000..e2c59ae9
--- /dev/null
+++ b/devtools/test/KernDerInterleaveTest.m
@@ -0,0 +1,219 @@
+KernDerInterleaveTest0()
+KernDerInterleaveTest1()
+KernDerInterleaveTest2()
+
+function KernDerInterleaveTest0()
+% test the relationship between interleaved Helmholtz kernels
+
+% problem parameters
+zk = 1;
+coefs = [1;0.5];
+coefa = [1,0.3;-1i,0.2];
+
+src = []; src.r = [0;-1]; src.n = [1;-2];
+targ = []; targ.r = [[1.1;0.3], [0.3;0.8]]; targ.n = [[-1;0.3],[-1;0.4]];
+
+% kernels
+skern = kernel('h','s',zk);
+dkern = kernel('h','d',zk);
+ckern = kernel('h','c',zk,coefs);
+
+spkern = kernel('h','sp',zk);
+dpkern = kernel('h','dp',zk);
+cpkern = kernel('h','cp',zk,coefs);
+
+allkern = kernel('h','all',zk,coefa);
+trkern = kernel('h','trans_rep',zk,coefs);
+trpkern = kernel('h','trans_rep_prime',zk,coefs);
+c2trkern = kernel('h','c2tr',zk,coefs);
+
+sgkern = kernel(@(s,t) chnk.helm2d.kern(zk,s,t,'sgrad'));
+dgkern = kernel(@(s,t) chnk.helm2d.kern(zk,s,t,'dgrad'));
+cgkern = kernel(@(s,t) chnk.helm2d.kern(zk,s,t,'cgrad',coefs));
+trgkern = kernel(@(s,t) chnk.helm2d.kern(zk,s,t,'trans_rep_grad',coefs));
+
+% evaluate kernels
+sval = skern.eval(src,targ);
+dval = dkern.eval(src,targ);
+cval = ckern.eval(src,targ);
+trval = trkern.eval(src,targ);
+
+spval = spkern.eval(src,targ);
+dpval = dpkern.eval(src,targ);
+cpval = cpkern.eval(src,targ);
+trpval = trpkern.eval(src,targ);
+
+sgval = sgkern.eval(src,targ);
+dgval = dgkern.eval(src,targ);
+cgval = cgkern.eval(src,targ);
+trgval = trgkern.eval(src,targ);
+
+allval = allkern.eval(src,targ);
+
+c2trval = c2trkern.eval(src,targ);
+
+% test combined
+assert(norm((coefs(1)*dval + coefs(2)*sval)-cval) <  1e-12)
+assert(norm((coefs(1)*dpval + coefs(2)*spval)-cpval) <  1e-12)
+
+% test all
+all_assemb = zeros(size(allval));
+all_assemb(1:2:end, 1:2:end) = coefa(1,1)*dval;
+all_assemb(1:2:end, 2:2:end) = coefa(1,2)*sval;
+all_assemb(2:2:end, 1:2:end) = coefa(2,1)*dpval;
+all_assemb(2:2:end, 2:2:end) = coefa(2,2)*spval;
+assert(norm(all_assemb - allval) <  1e-12)
+
+% test transmission representiatoon
+assert(norm([coefs(1)*dval , coefs(2)*sval]-trval) <  1e-12)
+assert(norm([coefs(1)*dpval , coefs(2)*spval]-trpval) <  1e-12)
+
+% test gradients
+assert(norm((coefs(1)*dgval + coefs(2)*sgval)-cgval) <  1e-12)
+assert(norm([coefs(1)*dgval , coefs(2)*sgval]-trgval) <  1e-12)
+
+% compare normal derivate with grad dot n
+dpval2 = sum(reshape(targ.n, 2,2).* reshape(dgval,2,[]),1);
+assert(norm(dpval - dpval2(:)) <  1e-12)
+
+% test c2trans
+c2trval_assemb = zeros(size(c2trval));
+c2trval_assemb(1:2:end, :) = coefs(1)*dval + coefs(2)*sval;
+c2trval_assemb(2:2:end, :) = coefs(1)*dpval + coefs(2)*spval;
+assert(norm(c2trval_assemb-c2trval) <  1e-12)
+end
+
+function KernDerInterleaveTest1()
+% test the relationship between interleaved Helmholtz difference kernels
+
+% problem parameters
+zk = [1,2];
+coefs = [[1;0.5],[1;0.5]];
+coefa = [1,0.3;-1i,0.2]; coefa = repmat(coefa,[1,1,2]);
+coefb = reshape(coefs,2,1,2); coefb = repmat(coefb,[1,2,1]);
+
+src = []; src.r = [0;-1]; src.n = [1;-2];
+targ = []; targ.r = [[1.1;0.3], [0.3;0.8]]; targ.n = [[-1;0.3],[-1;0.4]];
+
+% kernels
+skern = kernel('helmdiff','s',zk);
+dkern = kernel('helmdiff','d',zk);
+ckern = kernel('helmdiff','c',zk,coefs);
+
+spkern = kernel('helmdiff','sp',zk);
+dpkern = kernel('helmdiff','dp',zk);
+cpkern = kernel('helmdiff','cp',zk,coefs);
+
+allkern = kernel('helmdiff','all',zk,coefa);
+trkern = kernel('helmdiff','trans_rep',zk,coefs);
+trpkern = kernel('helmdiff','trans_rep_prime',zk,coefs);
+c2trkern = kernel('helmdiff','c2tr',zk,coefb);
+
+sgkern = kernel(@(s,t) chnk.helm2d.kern(zk(1),s,t,'sgrad_diff') ...
+    - chnk.helm2d.kern(zk(2),s,t,'sgrad_diff'));
+dgkern = kernel(@(s,t) chnk.helm2d.kern(zk(1),s,t,'dgrad_diff') ...
+    - chnk.helm2d.kern(zk(2),s,t,'dgrad_diff'));
+cgkern = kernel(@(s,t) chnk.helm2d.kern(zk(1),s,t,'cgrad_diff',coefs) ...
+    - chnk.helm2d.kern(zk(2),s,t,'cgrad_diff',coefs));
+trgkern = kernel(@(s,t) chnk.helm2d.kern(zk(1),s,t,'trans_rep_grad_diff',coefb(:,:,1)) ...
+    - chnk.helm2d.kern(zk(2),s,t,'trans_rep_grad_diff',coefb(:,:,2)));
+
+% evaluate kernels
+sval = skern.eval(src,targ);
+dval = dkern.eval(src,targ);
+cval = ckern.eval(src,targ);
+trval = trkern.eval(src,targ);
+
+spval = spkern.eval(src,targ);
+dpval = dpkern.eval(src,targ);
+cpval = cpkern.eval(src,targ);
+trpval = trpkern.eval(src,targ);
+
+sgval = sgkern.eval(src,targ);
+dgval = dgkern.eval(src,targ);
+cgval = cgkern.eval(src,targ);
+trgval = trgkern.eval(src,targ);
+
+allval = allkern.eval(src,targ);
+
+c2trval = c2trkern.eval(src,targ);
+
+% test combined
+assert(norm((coefs(1)*dval + coefs(2)*sval)-cval) <  1e-12)
+assert(norm((coefs(1)*dpval + coefs(2)*spval)-cpval) <  1e-12)
+
+% test all
+all_assemb = zeros(size(allval));
+all_assemb(1:2:end, 1:2:end) = coefa(1,1)*dval;
+all_assemb(1:2:end, 2:2:end) = coefa(1,2)*sval;
+all_assemb(2:2:end, 1:2:end) = coefa(2,1)*dpval;
+all_assemb(2:2:end, 2:2:end) = coefa(2,2)*spval;
+assert(norm(all_assemb - allval) <  1e-12)
+
+% test transmission representiatoon
+assert(norm([coefs(1)*dval , coefs(2)*sval]-trval) <  1e-12)
+assert(norm([coefs(1)*dpval , coefs(2)*spval]-trpval) <  1e-12)
+
+% test gradients
+assert(norm((coefs(1)*dgval + coefs(2)*sgval)-cgval) <  1e-12)
+assert(norm([coefs(1)*dgval , coefs(2)*sgval]-trgval) <  1e-12)
+
+% compare normal derivate with grad dot n
+dpval2 = sum(reshape(targ.n, 2,2).* reshape(dgval,2,[]),1);
+assert(norm(dpval - dpval2(:)) <  1e-12)
+
+% test c2trans
+c2trval_assemb = zeros(size(c2trval));
+c2trval_assemb(1:2:end, :) = coefs(1)*dval + coefs(2)*sval;
+c2trval_assemb(2:2:end, :) = coefs(1)*dpval + coefs(2)*spval;
+assert(norm(c2trval_assemb-c2trval) <  1e-12)
+end
+
+function KernDerInterleaveTest2()
+% test the relationship between interleaved Laplace kernels
+
+% problem parameters
+coefs = [1;0.5];
+coefa = [1,0.3;-1i,0.2];
+
+src = []; src.r = [0;-1]; src.n = [1;-2];
+targ = []; targ.r = [[1.1;0.3], [0.3;0.8]]; targ.n = [[-1;0.3],[-1;0.4]];
+
+% kernels
+skern = kernel('l','s');
+dkern = kernel('l','d');
+ckern = kernel('l','c',coefs);
+
+spkern = kernel('l','sp');
+dpkern = kernel('l','dp');
+cpkern = kernel('l','cp',coefs);
+
+sgkern = kernel(@(s,t) chnk.lap2d.kern(s,t,'sgrad'));
+dgkern = kernel(@(s,t) chnk.lap2d.kern(s,t,'dgrad'));
+cgkern = kernel(@(s,t) chnk.lap2d.kern(s,t,'cgrad',coefs));
+
+% evaluate kernels
+sval = skern.eval(src,targ);
+dval = dkern.eval(src,targ);
+cval = ckern.eval(src,targ);
+
+spval = spkern.eval(src,targ);
+dpval = dpkern.eval(src,targ);
+cpval = cpkern.eval(src,targ);
+
+sgval = sgkern.eval(src,targ);
+dgval = dgkern.eval(src,targ);
+cgval = cgkern.eval(src,targ);
+
+% test combined
+assert(norm((coefs(1)*dval + coefs(2)*sval)-cval) <  1e-12)
+assert(norm((coefs(1)*dpval + coefs(2)*spval)-cpval) <  1e-12)
+
+% test gradients
+assert(norm((coefs(1)*dgval + coefs(2)*sgval)-cgval) <  1e-12)
+
+% compare normal derivate with grad dot n
+dpval2 = sum(reshape(targ.n, 2,2).* reshape(dgval,2,[]),1);
+assert(norm(dpval - dpval2(:)) <  1e-12)
+
+end
\ No newline at end of file