diff --git a/chunkie/+chnk/+helm2d/kern.m b/chunkie/+chnk/+helm2d/kern.m
index 194b71e..33a3b00 100644
--- a/chunkie/+chnk/+helm2d/kern.m
+++ b/chunkie/+chnk/+helm2d/kern.m
@@ -40,7 +40,7 @@
 %                       [coef(1,1)*D coef(1,2)*S; coef(2,1)*D' coef(2,2)*S']
 %                type == 'c2trans' returns the combined field, and the 
 %                          normal derivative of the combined field
-%                        [coef(1)*D + coef(2)*S; coef(1)*D' + coef(2)*S']
+%                        [coef(1,1)*D + coef(1,2)*S; coef(2,1)*D' + coef(2,2)*S']
 %                type == 'trans_rep' returns the potential corresponding
 %                           to the transmission representation
 %                        [coef(1)*D coef(2)*S]
@@ -274,6 +274,7 @@
 case {'c2trans', 'c2t', 'c2tr'}
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end
+  if numel(coef) == 2, coef = repmat(coef(:).',2,1); end
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
 
@@ -296,8 +297,8 @@
 
   submat = zeros(2*nt,ns);
 
-  submat(1:2:2*nt,:) = coef(1)*submatd + coef(2)*submats;
-  submat(2:2:2*nt,:) = coef(1)*submatdp + coef(2)*submatsp;
+  submat(1:2:2*nt,:) = coef(1,1)*submatd + coef(1,2)*submats;
+  submat(2:2:2*nt,:) = coef(2,1)*submatdp + coef(2,2)*submatsp;
 
 
 % Dirichlet and Neumann data corresponding to combined field (difference)
diff --git a/chunkie/+chnk/+helm2dquas/green.m b/chunkie/+chnk/+helm2dquas/green.m
index 781b5ff..43e6770 100644
--- a/chunkie/+chnk/+helm2dquas/green.m
+++ b/chunkie/+chnk/+helm2dquas/green.m
@@ -1,5 +1,5 @@
-function [val,grad,hess] = green(src,targ,zk,kappa,d,sn,l)
-%CHNK.HELM2DQUAS.GREEN evaluate the quasiperiodic helmholtz green's function
+function [val,grad,hess] = green(src,targ,zk,kappa,d,sn,l,ising)
+%CHNK.HELM2DQUAS.GREEN evaluate the quasiperiodic Helmholtz Green's function
 % for the given sources and targets
 %
 % Input:
@@ -8,12 +8,14 @@
 %   d - period
 %   sn - precomputed lattice sum integrals 
 %       (see chnk.helm2dquas.latticecoefs)
+%   l - number of periodic copies computed explicitly
+%   ising - if set to 0, only include the periodic copies. If set to 1,
+%       include the free-space part
 %
 % see also CHNK.HELM2DQUAS.KERN
 [~,nsrc] = size(src);
 [~,ntarg] = size(targ);
 
-
 xs = repmat(src(1,:),ntarg,1);
 ys = repmat(src(2,:),ntarg,1);
 
@@ -23,6 +25,11 @@
 rx = xt-xs;
 ry = yt-ys;
 
+
+rx = rx(:);
+ry = ry(:);
+
+
 nx = fix(rx/d);
 rx = rx - nx*d;
 
@@ -34,14 +41,9 @@
 
 r = sqrt(r2);
 
-rx = rx(:);
-ry = ry(:);
-r = r(:);
-
 npt = size(r,1);
 
-ythresh = d/2;
-% ythresh = d/10;
+ythresh = 2*d/2;
 iclose = abs(ry) < ythresh;
 ifar = ~iclose;
 
@@ -52,6 +54,9 @@
 ryclose = ry(iclose);
 rclose = r(iclose);
 
+nxclose = nx(iclose);
+nptclose = size(rxclose, 1);
+
 nkappa = length(kappa);
 
 val = zeros(nkappa,npt,1);
@@ -76,55 +81,57 @@
 % beta = sqrt((xi_m.^2-zk^2));
 beta = sqrt(1i*(xi_m-zk)).*sqrt(-1i*(xi_m+zk));
 
-% fhat = exp(-beta.*sqrt(ryfar.^2))./beta.*exp(1i*xi_m.*rxfar)/2;
 fhat = exp(-beta.*sqrt(ryfar.^2) + 1i*xi_m.*rxfar)./(2*beta);
 val(:,ifar,:) = sum(fhat,3)/(d);
 if nargout > 1
 grad(:,ifar,1) = sum(1i*xi_m.*fhat,3)/d;
 grad(:,ifar,2) = sum(-beta.*(sqrt(ryfar.^2)./ryfar).*fhat,3)/d;
-% grad(:,ifar,:) =[gx,gy];
 end
 
 if nargout >2
 hess(:,ifar,1) = sum(-xi_m.^2.*fhat,3)/d;
 hess(:,ifar,2) = sum(-1i*xi_m.*beta.*(sqrt(ryfar.^2)./ryfar).*fhat,3)/d;
 hess(:,ifar,3) = sum((beta.*(sqrt(ryfar.^2)./ryfar)).^2.*fhat,3)/d;
-% hess(:,ifar,:) = [hxx,hxy,hyy];
 end
 
-
 end
 
 alpha = (exp(1i*kappa(:)*d));
 
-val_near=0;
-grad_near = zeros(1,1,2);
-hess_near = zeros(1,1,3);
+val_near= zeros(nkappa,nptclose);
+grad_near = zeros(nkappa,nptclose,2);
+hess_near = zeros(nkappa,nptclose,3);
+ls = -l:l;
 if ~isempty(rxclose)
-    for i = -l:l
+    for i = ls
+        if ising == 1
+            iuse = true(nptclose,1);
+        else
+            iuse = nxclose ~= -i;
+        end
+
         rxi = rxclose - i*d;
         if nargout>2
         [vali,gradi,hessi] = chnk.helm2d.green(zk,[0;0],[rxi.';ryclose.']);
         vali = reshape(vali,1,[],1);
         gradi = reshape(gradi,1,[],2);
         hessi = reshape(hessi,1,[],3);
-        val_near = val_near + vali.*alpha.^i;
-        grad_near = grad_near + gradi.*alpha.^i;
-        hess_near = hess_near + hessi.*alpha.^i;
+        val_near(:,iuse) = val_near(:,iuse) + vali(:,iuse).*alpha.^i;
+        grad_near(:,iuse,:) = grad_near(:,iuse,:) + gradi(:,iuse,:).*alpha.^i;
+        hess_near(:,iuse,:) = hess_near(:,iuse,:) + hessi(:,iuse,:).*alpha.^i;
         elseif nargout > 1
         [vali,gradi] = chnk.helm2d.green(zk,[0;0],[rxi.';ryclose.']);
         vali = reshape(vali,1,[],1);
         gradi = reshape(gradi,1,[],2);
-        val_near = val_near + vali.*alpha.^i;
-        grad_near = grad_near + gradi.*alpha.^i;
+        val_near(:,iuse) = val_near(:,iuse) + vali(:,iuse).*alpha.^i;
+        grad_near(:,iuse,:) = grad_near(:,iuse,:) + gradi(:,iuse,:).*alpha.^i;
         else
         vali = chnk.helm2d.green(zk,[0;0],[rxi.';ryclose.']);
         vali = reshape(vali,1,[],1);
-        val_near = val_near + vali.*alpha.^i;
+        val_near(:,iuse) = val_near(:,iuse) + vali(:,iuse).*alpha.^i;
         end
     end
-    
-    % sn = sn.';
+
     N = size(sn,2)-1;
     ns = (0:N);
     ns_use = (0:N+2);
@@ -139,6 +146,7 @@
         Js(:,i) = besselj(ns_use(i),zk*rclose);
         end
     end
+    % t1 = tic;
     eip = (rxclose+1i*ryclose)./rclose;
     eipn = reshape(eip.^ns,1,[], N+1);
     cs = (eipn+1./eipn)/2;
@@ -146,18 +154,18 @@
     Js = reshape(Js,1,[],N+3);
     sn = reshape(sn,nkappa, 1, N+1);
     
-    val_far = 0.25*1i*Js(:,:,1).*sn(:,:,1) + 0.5*1i*sum(sn(:,:,2:end).*Js(:,:,2:end-2).*cs(:,:,2:end),3);
+    tmp = reshape(Js(:,:,2:end-2).*cs(:,:,2:end),[],N);
+    val_far = 0.25*1i*Js(:,:,1).*sn(:,:,1) + 0.5*1i*sn(:,2:end)*tmp.';
     val(:,iclose) = val_near+val_far;
     
-    
-    
-    
     if nargout >1
         DJs = cat(3,-Js(:,:,2),.5*(Js(:,:,1:end-3)-Js(:,:,3:end-1)))*zk;
         ss = (eipn-1./eipn)/2i;
             
-        grad_far_p = 0.25*1i*DJs(:,:,1).*sn(:,:,1) + 0.5*1i*sum(sn(:,:,2:end).*DJs(:,:,2:end).*cs(:,:,2:end),3);
-        grad_far_t = (0.5*1i*sum(-reshape((1:N),1,1,[]).*sn(:,:,2:end).*Js(:,:,2:end-2).*ss(:,:,2:end),3))./rclose.';
+        tmp = reshape(DJs(:,:,2:end).*cs(:,:,2:end),[],N);
+        grad_far_p = 0.25*1i*DJs(:,:,1).*sn(:,:,1) + 0.5*1i*sn(:,2:end)*tmp.';
+        tmp = reshape(Js(:,:,2:end-2).*ss(:,:,2:end),[],N)./rclose;
+        grad_far_t = (0.5*1i*((-reshape((1:N),1,[]).*sn(:,2:end))*tmp.'));
         
         grad_far = cat(3,cs(:,:,2).*grad_far_p - ss(:,:,2).*grad_far_t, ss(:,:,2).*grad_far_p + cs(:,:,2).*grad_far_t);
         grad(:,iclose,:) = grad_near + grad_far; 
@@ -172,37 +180,85 @@
         tmp_n = rclose.^(-4).*(-ns.*ryclose.*Js(:,:,1:end-2).*(ns.*ryclose.*cs+2*rxclose.*ss)+ ...
             rclose.*ryclose.*(ryclose.*cs + 2*ns.*rxclose.*ss).*DJs+ ...
             rclose.^2.*rxclose.^2.*cs.*DDJs);
-        hess_far_xx = 0.25*1i*tmp_n(:,:,1).*sn(:,:,1)+.5*1i*sum(sn(:,:,2:end).*tmp_n(:,:,2:end),3);
+        tmp_n = reshape(tmp_n,[],N+1);
+        hess_far_xx = 0.25*1i*tmp_n(:,1).'.*sn(:,1)+.5*1i*sn(:,2:end)*tmp_n(:,2:end).';
 
         tmp_n = ns.*Js(:,:,1:end-2).*(ns.*rxclose.*ryclose.*cs+(rxclose.^2-ryclose.^2).*ss).*rclose.^(-4) ...
             +rclose.^(-3).*(-(rxclose.*ryclose.*cs+ns.*(rxclose.^2-ryclose.^2).*ss).*DJs+...
             rclose.*rxclose.*ryclose.*cs.*DDJs);
 
-
-        hess_far_xy = 0.25*1i*tmp_n(:,:,1).*sn(:,:,1)+.5*1i*sum(sn(:,:,2:end).*tmp_n(:,:,2:end),3);
+        tmp_n = reshape(tmp_n,[],N+1);
+        hess_far_xy = 0.25*1i*tmp_n(:,1).'.*sn(:,1)+.5*1i*sn(:,2:end)*tmp_n(:,2:end).';
 
         tmp_n = rclose.^(-4).*(-ns.*rxclose.*Js(:,:,1:end-2).*(ns.*rxclose.*cs-2*ryclose.*ss)+ ...
             rclose.*rxclose.*(rxclose.*cs - 2*ns.*ryclose.*ss).*DJs+ ...
             rclose.^2.*ryclose.^2.*cs.*DDJs);
-        hess_far_yy = 0.25*1i*tmp_n(:,:,1).*sn(:,:,1)+.5*1i*sum(sn(:,:,2:end).*tmp_n(:,:,2:end),3);
+        tmp_n = reshape(tmp_n,[],N+1);
+        hess_far_yy = 0.25*1i*tmp_n(:,1).'.*sn(:,1)+.5*1i*sn(:,2:end)*tmp_n(:,2:end).';
 
         hess_far = cat(3,hess_far_xx, hess_far_xy, hess_far_yy);
 
-
         hess(:,iclose,:) = hess_near + hess_far;
     end
 end
 
 quasi_phase = exp(1i*kappa(:)*nx(:).'*d);
 
-val = reshape(quasi_phase.*val,nkappa*ntarg,nsrc);
-if nargout>1
-grad = reshape(quasi_phase.*grad,nkappa*ntarg,nsrc,2);
-end
-if nargout>2
-hess = reshape(quasi_phase.*hess,nkappa*ntarg,nsrc,3);
-end
-% end
 
-end
+if nargout == 1
+    val = quasi_phase.*val;
+
+    if ising == 0
+        isub = (abs(nx(:)) > max(ls)) | ifar;
+
+        if any(isub)
+        vali = chnk.helm2d.green(zk,[0;0],[rx(isub).'+ nx(isub).'*d;ry(isub).']);
+        vali = reshape(vali,1,[],1);
+        val(:,isub,:) = val(:,isub,:) - vali;
+        end
+    end
 
+    val = reshape(val,nkappa*ntarg,nsrc);
+elseif nargout == 2
+    val = quasi_phase.*val;
+    grad = quasi_phase.*grad;
+    
+    if ising == 0
+        isub = (abs(nx(:)) > max(ls)) | ifar;
+    
+        if any(isub)
+        [vali, gradi] = chnk.helm2d.green(zk,[0;0],[rx(isub).' + nx(isub).'*d;ry(isub).']);
+        vali = reshape(vali,1,[],1);
+        gradi = reshape(gradi,1,[],2);
+        val(:,isub,:) = val(:,isub,:) - vali;
+        grad(:,isub,:) = grad(:,isub,:) - gradi;
+        end
+    end
+
+    val = reshape(val,nkappa*ntarg,nsrc);
+    grad = reshape(grad,nkappa*ntarg,nsrc,2);
+elseif nargout == 3
+    val = quasi_phase.*val;
+    grad = quasi_phase.*grad;
+    hess = quasi_phase.*hess;
+    
+    if ising == 0
+        isub = (abs(nx(:)) > max(ls)) | ifar;
+
+        if any(isub)
+        [vali, gradi, hessi] = chnk.helm2d.green(zk,[0;0],[rx(isub).' + nx(isub).'*d;ry(isub).']);
+        vali = reshape(vali,1,[],1);
+        gradi = reshape(gradi,1,[],2);
+        hessi = reshape(hessi,1,[],3);
+
+        val(:,isub,:) = val(:,isub,:) - vali;
+        grad(:,isub,:) = grad(:,isub,:) - gradi;
+        hess(:,isub,:) = hess(:,isub,:) - hessi;
+        end
+    end
+
+    val = reshape(val,nkappa*ntarg,nsrc);
+    grad = reshape(grad,nkappa*ntarg,nsrc,2);
+    hess = reshape(hess,nkappa*ntarg,nsrc,3);
+end
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+helm2dquas/kern.m b/chunkie/+chnk/+helm2dquas/kern.m
index 9cc5838..3ee7a63 100644
--- a/chunkie/+chnk/+helm2dquas/kern.m
+++ b/chunkie/+chnk/+helm2dquas/kern.m
@@ -46,7 +46,7 @@
 %                       [coef(1,1)*D coef(1,2)*S; coef(2,1)*D' coef(2,2)*S']
 %                type == 'c2trans' returns the combined field, and the 
 %                          normal derivative of the combined field
-%                        [coef(1)*D + coef(2)*S; coef(1)*D' + coef(2)*S']
+%                        [coef(1,1)*D + coef(1,2)*S; coef(2,1)*D' + coef(2,2)*S']
 %                type == 'trans_rep' returns the potential corresponding
 %                           to the transmission representation
 %                        [coef(1)*D coef(2)*S]
@@ -65,6 +65,8 @@
 %       quas_param.sn - precomputed lattice sum integrals, see 
 %               chnk.helm2dquas.latticecoefs
 %
+%   varargin{1} - ising: if set to 0, only include the periodic copies. 
+%                 If set to 1, include the free-space part
 %   varargin{1} - coef: length 2 array in the combined layer 
 %                 formula, 2x2 matrix for all kernels
 %                 otherwise does nothing
@@ -80,6 +82,12 @@
 src = srcinfo.r(:,:);
 targ = targinfo.r(:,:);
 
+
+ising = 1;
+if length(varargin) >1
+    ising = varargin{2};
+end 
+
 kappa = quas_param.kappa;
 d = quas_param.d;
 l = quas_param.l;
@@ -94,7 +102,7 @@
 % double layer
 case {'d', 'double'}
   srcnorm = srcinfo.n(:,:);
-  [~,grad] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
+  [~,grad] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l,ising);
   nx = repmat(srcnorm(1,:),nkappa*nt,1);
   ny = repmat(srcnorm(2,:),nkappa*nt,1);
   submat = -(grad(:,:,1).*nx + grad(:,:,2).*ny);
@@ -102,7 +110,7 @@
 % normal derivative of single layer
 case {'sp', 'sprime'}
   targnorm = targinfo.n(:,:);
-  [~,grad] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
+  [~,grad] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l,ising);
   nx = repmat(reshape(targnorm(1,:),1,nt,1),nkappa,1,ns);
   nx = reshape(nx,[],ns);
   ny = repmat(reshape(targnorm(2,:),1,nt,1),nkappa,1,ns);
@@ -110,16 +118,18 @@
 
   submat = (grad(:,:,1).*nx + grad(:,:,2).*ny);
 
+
+
 % single later
 case {'s', 'single'}
-  submat = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
+  submat = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l,ising);
   submat = reshape(submat,[],ns);
 
 % normal derivative of double layer
 case {'dp', 'dprime'}
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
-  [~,~,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
+  [~,~,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l,ising);
   nxsrc = repmat(srcnorm(1,:),nkappa*nt,1);
   nysrc = repmat(srcnorm(2,:),nkappa*nt,1);
   % nxtarg = repmat((targnorm(1,:)).',1,ns);
@@ -134,21 +144,21 @@
 
 % gradient of single layer
 case {'sgrad','sg'}
-    [~,grad] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
-    submat = reshape(permute(grad,[3,1,2]),[],ns);
+    [~,grad] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l,ising);
+    submat = reshape(permute(grad,[1,3,2]),[],ns);
 
 % gradient of double layer
 case {'dgrad','dg'}
-    [~,~,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
+    [~,~,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l,ising);
     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);
+    submat = reshape(permute(submat,[1,3,2]),[],ns);
 
 % Combined field 
 case {'c', 'combined'}
   srcnorm = srcinfo.n(:,:);
   coef = ones(2,1);
-  if(nargin == 6); coef = varargin{1}; end
-  [submats,grad] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
+  if(nargin >= 6); coef = varargin{1}; end
+  [submats,grad] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l,ising);
   nx = repmat(srcnorm(1,:),nkappa*nt,1);
   ny = repmat(srcnorm(2,:),nkappa*nt,1);
   submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny);
@@ -157,12 +167,12 @@
 % normal derivative of combined field
 case {'cp', 'cprime'}
   coef = ones(2,1);
-  if(nargin == 6); coef = varargin{1}; end
+  if(nargin >= 6); coef = varargin{1}; end
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
 
   % Get gradient and hessian info
-  [~,grad,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
+  [~,grad,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l,ising);
 
   nxsrc = repmat(srcnorm(1,:),nkappa*nt,1);
   nysrc = repmat(srcnorm(2,:),nkappa*nt,1);
@@ -180,15 +190,18 @@
 
   submat = coef(1)*submatdp + coef(2)*submatsp;
 
+
 % Dirichlet and Neumann data corresponding to combined field
 case {'c2trans', 'c2t', 'c2tr'}
   coef = ones(2,1);
-  if(nargin == 6); coef = varargin{1}; end
+  if(nargin >= 6); coef = varargin{1}; end
+  if numel(coef) == 2, coef = repmat(coef(:).',2,1); end
+
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
 
   % Get gradient and hessian info
-  [submats,grad,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
+  [submats,grad,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l,ising);
 
   nxsrc = repmat(srcnorm(1,:),nkappa*nt,1);
   nysrc = repmat(srcnorm(2,:),nkappa*nt,1);
@@ -206,21 +219,25 @@
   % D
   submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc);
 
-  submat = zeros(nkappa*2*nt,ns);
+  % submat = zeros(nkappa*2*nt,ns);
+  % submat(1:2:end,:) = coef(1)*submatd + coef(2)*reshape(submats,[],ns);
+  % submat(2:2:end,:) = coef(1)*submatdp + coef(2)*submatsp;
 
-  submat(1:2:end,:) = coef(1)*submatd + coef(2)*reshape(submats,[],ns);
-  submat(2:2:end,:) = coef(1)*submatdp + coef(2)*submatsp;
+  submat = zeros(nkappa,2,nt,ns);
+  submat(:,1,:,:) = reshape(coef(1,1)*submatd + coef(1,2)*reshape(submats,[],ns),nkappa,1,nt,[]);
+  submat(:,2,:,:) = reshape(coef(2,1)*submatdp + coef(2,2)*submatsp,nkappa,1,nt,[]);
+  submat = reshape(submat, [],ns);
 
 % gradient of combined field
 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);
+    if(nargin >= 6); coef = varargin{1}; end
+    [~,grad,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l,ising);
 
-    submats = reshape(permute(grad,[3,1,2]),[],ns);
+    submats = reshape(permute(grad,[1,3,2]),[],ns);
 
     submatd = -(hess(:,:,1:2).*srcinfo.n(1,:)+hess(:,:,2:3).*srcinfo.n(2,:));
-    submatd = reshape(permute(submatd,[3,1,2]),[],ns);
+    submatd = reshape(permute(submatd,[1,3,2]),[],ns);
     
     submat = coef(1)*submatd + coef(2)*submats;
   
@@ -233,7 +250,7 @@
   cc = varargin{1};
 
   % Get gradient and hessian info
-  [submats,grad,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
+  [submats,grad,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l,ising);
 
   nxsrc = repmat(srcnorm(1,:),nkappa*nt,1);
   nysrc = repmat(srcnorm(2,:),nkappa*nt,1);
@@ -252,44 +269,57 @@
   % D
   submatd  = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc);
   
-  submat = zeros(nkappa*2*nt,2*ns);
-  submat(1:2:end,1:2:2*ns) = submatd*cc(1,1);
-  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);
+  submat = zeros(nkappa,2,nt,2*ns);
+  submat(:,1,:,1:2:2*ns) = reshape(submatd,nkappa,1,nt,[])*cc(1,1);
+  submat(:,1,:,2:2:2*ns) = reshape(submats,nkappa,1,nt,[])*cc(1,2);
+  submat(:,2,:,1:2:2*ns) = reshape(submatdp,nkappa,1,nt,[])*cc(2,1);
+  submat(:,2,:,2:2:2*ns) = reshape(submatsp,nkappa,1,nt,[])*cc(2,2);
+  submat = reshape(submat, [],2*ns);
+% end
+% % Dirichlet data/potential correpsonding to transmission rep
+% if strcmpi(type,'trans_rep') 
+% =======
+%   submat = zeros(nkappa*2*nt,2*ns);
+%   submat(1:2:end,1:2:2*ns) = submatd*cc(1,1);
+%   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);
+% >>>>>>> master
 
 % Dirichlet data/potential correpsonding to transmission rep
 case {'trans_rep','trep'} 
   coef = ones(2,1);
-  if(nargin == 6); coef = varargin{1}; end;
+  if(nargin >= 6); coef = varargin{1}; end;
   srcnorm = srcinfo.n(:,:);
-  [submats,grad,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
+  [submats,grad,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l,ising);
 
   submats = reshape(submats,[],ns);
-  nx = repmat(srcnorm(1,:),nt,1);
-  ny = repmat(srcnorm(2,:),nt,1);
-  submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny);
+  nxsrc = repmat(srcnorm(1,:),nkappa*nt,1);
+  nysrc = repmat(srcnorm(2,:),nkappa*nt,1);
+  submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc);
 
-  submat = zeros(nkappa*nt,2*ns);
+  submat = zeros(nkappa*nt,2*ns,'like',1i);
   submat(1:1:end,1:2:2*ns) = coef(1)*submatd;
   submat(1:1:end,2:2:2*ns) = coef(2)*submats;
 
 % Neumann data corresponding to transmission rep
 case {'trans_rep_prime','trep_p', 'trans_rep_p'}
   coef = ones(2,1);
-  if(nargin == 6); coef = varargin{1}; end;
+  if(nargin >= 6); coef = varargin{1}; end;
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
   
   submat = zeros(nt,ns);
 
   % Get gradient and hessian info
-  [submats,grad,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
+  [submats,grad,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l,ising);
 
-  nxsrc = repmat(srcnorm(1,:),nt,1);
-  nysrc = repmat(srcnorm(2,:),nt,1);
-  nxtarg = repmat((targnorm(1,:)).',1,ns);
-  nytarg = repmat((targnorm(2,:)).',1,ns);
+  nxsrc = repmat(srcnorm(1,:),nkappa*nt,1);
+  nysrc = repmat(srcnorm(2,:),nkappa*nt,1);
+  nxtarg = repmat(reshape(targnorm(1,:),1,nt,1),nkappa,1,ns);
+  nxtarg = reshape(nxtarg,[],ns);
+  nytarg = repmat(reshape(targnorm(2,:),1,nt,1),nkappa,1,ns);
+  nytarg = reshape(nytarg,[],ns);
 
   % D'
   submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ...
@@ -304,28 +334,28 @@
 % Gradient correpsonding to transmission rep
 case {'trans_rep_grad','trep_g', 'trans_rep_g'}
   coef = ones(2,1);
-  if(nargin == 6); coef = varargin{1}; end;
+  if(nargin >= 6); coef = varargin{1}; end;
   
   srcnorm = srcinfo.n(:,:);
-  
-  % submat = zeros(nt,ns,6);
-  % S
-  [submats,grad,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
 
-  nxsrc = repmat(srcnorm(1,:),nt,1);
-  nysrc = repmat(srcnorm(2,:),nt,1);
-  % D
-  submatd  = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc);
+  [submats,grad,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l,ising);
 
-  submat = zeros(nkappa*2*nt,2*ns);
+  nxsrc = repmat(srcnorm(1,:),nkappa*nt,1);
+  nysrc = repmat(srcnorm(2,:),nkappa*nt,1);
+
+  submat = zeros(nkappa,2,nt,2*ns);
   
-  submat(1:2:end,1:2:2*ns) = -coef(1)*(hess(:,:,1).*nxsrc + hess(:,:,2).*nysrc);
-  submat(1:2:end,2:2:2*ns) = coef(2)*grad(:,:,1);
+  submat(:,1,:,1:2:2*ns) = reshape(-coef(1)*(hess(:,:,1).*nxsrc + hess(:,:,2).*nysrc),nkappa,1,nt,ns);
+  submat(:,1,:,2:2:2*ns) = reshape(coef(2)*grad(:,:,1),nkappa,1,nt,ns);
     
-  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);
+  submat(:,2,:,1:2:2*ns) = reshape(-coef(1)*(hess(:,:,2).*nxsrc + hess(:,:,3).*nysrc),nkappa,1,nt,ns);
+  submat(:,2,:,2:2:2*ns) = reshape(coef(2)*grad(:,:,2),nkappa,1,nt,ns);
+  submat = reshape(permute(submat,[1,3,2,4]),nkappa*2*nt,2*ns);
+  
+
 otherwise
     error('Unknown quasi-periodic Helmholtz kernel type ''%s''.', type);
+
 end
 
 
diff --git a/chunkie/@chunker/chunker.m b/chunkie/@chunker/chunker.m
index 3e3aa6c..9ce6fbe 100644
--- a/chunkie/@chunker/chunker.m
+++ b/chunkie/@chunker/chunker.m
@@ -56,6 +56,7 @@
 %   quiver(obj,varargin) - quiver plot of the chnkr points and normals
 %   scatter(obj,varargin) - scatter plot of the chnkr nodes
 %   tau = taus(obj) - unit tangents to curve
+%   cgrph = tochunkgraph(obj) - convert to chunkgraph
 %   obj = refine(obj,varargin) - refine the curve
 %   a = area(obj) - for a closed curve, area inside 
 %   s = arclength(obj) - get values of arclength along curve
@@ -363,6 +364,7 @@
         quiver(obj,varargin)
         scatter(obj,varargin)
         tau = taus(obj)
+        cgrph = tochunkgraph(obj)
         obj = refine(obj,varargin)
         a = area(obj)
         s = arclength(obj)
diff --git a/chunkie/@chunker/tochunkgraph.m b/chunkie/@chunker/tochunkgraph.m
new file mode 100644
index 0000000..ea49e94
--- /dev/null
+++ b/chunkie/@chunker/tochunkgraph.m
@@ -0,0 +1,50 @@
+function cgrph = tochunkgraph(chnkr)
+% TOCHUNKGRAPH converts a chunker into a chunkgraph with one edge per
+% connected component
+
+chnkr = sort(chnkr);
+[~,~,info] = sortinfo(chnkr);
+
+ncomp = info.ncomp;
+nchs = info.nchs;
+
+istart = 1;
+nverts = 0;
+verts = [];
+edge2verts = [];
+cparams = cell(1,ncomp);
+fchnks = cell(1,ncomp);
+
+for i = 1:ncomp
+    nch = nchs(i);
+    vs = chunkends(chnkr,[istart,istart+nch-1]);
+    
+    ifclosed = info.ifclosed(i);
+    if ifclosed
+        vs = vs(:,1);
+        verts = [verts, vs];
+        edge2verts = [edge2verts, [nverts+1;nverts+1]];
+        nverts = nverts + 1;
+    else
+        vs = vs(:, [1,4]);
+        verts = [verts, vs];
+        edge2verts = [edge2verts, [nverts+1;nverts+2]];
+        nverts = nverts + 2;
+    end
+
+    rchnkr = []; 
+    rchnkr.r  = chnkr.r(:,:,istart:(istart+nch-1));
+    rchnkr.d  = chnkr.d(:,:,istart:(istart+nch-1));
+    rchnkr.d2 = chnkr.d2(:,:,istart:(istart+nch-1));
+    fchnks{i} = chunkerpoints(rchnkr,struct('ifclosed',ifclosed));
+    cparams{i}.ifclosed = ifclosed;
+    istart = istart + nch;
+end
+
+pref = [];
+
+pref.dim = chnkr.dim;
+pref.k   = chnkr.k;
+cgrph = chunkgraph(verts,edge2verts,fchnks,cparams,pref);
+
+end
\ No newline at end of file
diff --git a/chunkie/@chunkgraph/chunkgraph.m b/chunkie/@chunkgraph/chunkgraph.m
index fe5809b..d719bd1 100644
--- a/chunkie/@chunkgraph/chunkgraph.m
+++ b/chunkie/@chunkgraph/chunkgraph.m
@@ -67,7 +67,6 @@
 %           using bounding rectangles
 %   flag = flagnear_rectangle_grid(obj,x,y) - find points defined via 
 %      a meshgrid of points that are close to the chunkgraph
-%   obj = obj.move(r0,r1,trotat,scale) - translate, rotate, etc
 %   rmin = min(obj) - minimum of coordinate values
 %   rmax = max(obj) - maximum of coordinate values
 %   onesmat = onesmat(obj) - matrix that corresponds to integration of a
@@ -94,14 +93,14 @@
 %   edgesendverts - (2 x nedges) array specifying starting and ending
 %     vertices of an edge. 
 % Optional input:
-%   fchnks - (nedges x 1) cell array of function handles, specifying a
-%     smooth curve to connect the given vertices. fchnk{j} should be a
-%     function in the format expected by CHUNKERFUNC but with the default
-%     that fchnk{j} is a function from [0,1] to the curve. If the specified
-%     curve does not actually connect the given vertices, the curve will be
-%     translated, rotated, and scaled to connect them. If the specified
-%     curve should be a loop, it is only translated to start at the correct
-%     vertex. 
+%   fchnks - (1 x nedges) cell array of function handles or chunkers,
+%     specifying a smooth curve to connect the given vertices. If fchnk{j}
+%     is a function handle, it should be a function in the format expected
+%     by CHUNKERFUNC but with the default that fchnk{j} is a function from
+%     [0,1] to the curve. If the specified curve does not actually connect
+%     the given vertices, the curve will be translated, rotated, and scaled
+%     to connect them. If the specified curve should be a loop, it is only
+%     translated to start at the correct vertex.
 %   cparams - struct or (nedges x 1) cell array of structs specifying curve
 %     parameters in the format expected by CHUNKERFUNC. 
 %   pref - struct specifying CHUNKER preferences. 
@@ -179,7 +178,7 @@
                 fchnks = [];
             end
             
-            if isa(fchnks,"function_handle")
+            if isa(fchnks,"function_handle") || isa(fchnks,"chunker")
                 fchnks0 = fchnks;
                 fchnks = cell(nedge,1);
                 for j = 1:nedge
@@ -242,7 +241,7 @@
                     chnkr = sort(chnkr);
                     %chnkr.vert = [v1,v2];
                     echnks(i) = chnkr;
-                elseif (~isempty(fchnks{i}) && isa(fchnks{i},'function_handle'))
+                else
                     if iscell(cparams)
                         cploc = cparams{i};
                     else
@@ -250,46 +249,71 @@
                     end
                     % set cploc.ifclosed in a way that makes sense
                     cploc.ifclosed = false;
+
                     i1 = edgesendverts(1,i);
                     i2 = edgesendverts(2,i);
-
                     if isnan(i1) || isnan(i2)
                         cploc.ifclosed = true;
                     end
 
-                    % chunkgraph edges need at least 4 chunks
-                    if isfield(cploc,'nchmin')
-                        cploc.nchmin = max(4,cploc.nchmin);
-                    else
-                        cploc.nchmin = 4;
-                    end
-                    
-                    ta = 0; tb = 1;
-                    if isfield(cploc,'ta')
-                        ta = cploc.ta; 
-                    else
-                        cploc.ta = ta;
-                    end
-                    if isfield(cploc,'tb')
-                        tb = cploc.tb; 
-                    else
-                        cploc.tb = tb;
-                    end
 
-                    vs =fchnks{i}([ta,tb]);
-                    if isfield(cploc,'maxchunklen')
-                        if ~isnan(i1) && ~isnan(i2)
-                            if i1 ~= i2 
-                                vfin0 = verts(:,i1);
-                                vfin1 = verts(:,i2);
-                                scale = norm(vfin1-vfin0,'fro')/norm(vs(:,2)-vs(:,1),'fro');
-                                cploc.maxchunklen = cploc.maxchunklen/scale;
+                    if isa(fchnks{i},'function_handle')
+
+    
+                        % chunkgraph edges need at least 4 chunks
+                        if isfield(cploc,'nchmin')
+                            cploc.nchmin = max(4,cploc.nchmin);
+                        else
+                            cploc.nchmin = 4;
+                        end
+                        
+                        ta = 0; tb = 1;
+                        if isfield(cploc,'ta')
+                            ta = cploc.ta; 
+                        else
+                            cploc.ta = ta;
+                        end
+                        if isfield(cploc,'tb')
+                            tb = cploc.tb; 
+                        else
+                            cploc.tb = tb;
+                        end
+    
+                        vs =fchnks{i}([ta,tb]);
+                        if isfield(cploc,'maxchunklen')
+                            if ~isnan(i1) && ~isnan(i2)
+                                if i1 ~= i2 
+                                    vfin0 = verts(:,i1);
+                                    vfin1 = verts(:,i2);
+                                    scale = norm(vfin1-vfin0,'fro')/norm(vs(:,2)-vs(:,1),'fro');
+                                    cploc.maxchunklen = cploc.maxchunklen/scale;
+                                end
                             end
                         end
+                        
+                        chnkr = chunkerfunc(fchnks{i}, cploc, pref);
+
+                    elseif isa(fchnks{i},'chunker')
+                        chnkr = fchnks{i};
+                        [~,~,info] = sortinfo(chnkr);
+                        if info.ncomp > 1
+                            msg = "CHUNKGRAPH:CONSTRUCTOR: invalid edge descriptor." + ...
+                                  " fchnks{" + num2str(i) + "} must be a " + ...
+                                  "single connected component.";
+                            warning(msg);
+                        end
+
+                        vs = chunkends(chnkr,[1,chnkr.nch]);
+                        vs = vs(:, [1,4]);
+
+                    else
+                        msg = "CHUNKGRAPH:CONSTRUCTOR: invalid edge descriptor." + ...
+                              " fchnks{" + num2str(i) + "} must be empty, a" + ...
+                              " function handle, or a chunker.";
+                        warning(msg);
                     end
-                    
-                    chnkr = chunkerfunc(fchnks{i}, cploc, pref);
                     chnkr = sort(chnkr);
+
                     if ~isnan(i1) && ~isnan(i2)
                         if i1 ~= i2 
                             vfin0 = verts(:,i1);
@@ -317,7 +341,6 @@
                             chnkr = move(chnkr,zeros(size(vs(:,1))),vfin-vs(:,1),0,1);
                         end
                     end
-
                     echnks(i) = chnkr;
                 end
             end
diff --git a/chunkie/@kernel/conj.m b/chunkie/@kernel/conj.m
new file mode 100644
index 0000000..8d84bbc
--- /dev/null
+++ b/chunkie/@kernel/conj.m
@@ -0,0 +1,22 @@
+function f = conj(f)
+% take a complex conjugate of a kernel class
+%
+
+if(isa(f.shifted_eval, 'function_handle'))
+    f.shifted_eval = @(varargin) conj(f.shifted_eval(varargin{:}));
+else
+    f.shifted_eval = [];
+
+end
+
+f.eval = @(varargin) conj(f.eval(varargin{:}));
+
+if (isa(f.fmm,'function_handle'))
+    f.fmm = @(varargin) conj(f.fmm(varargin{:}));
+else
+    f.fmm = [];
+end
+
+end
+
+
diff --git a/chunkie/@kernel/helm2dquas.m b/chunkie/@kernel/helm2dquas.m
index 879c5e5..0c4d7dc 100644
--- a/chunkie/@kernel/helm2dquas.m
+++ b/chunkie/@kernel/helm2dquas.m
@@ -1,4 +1,4 @@
-function obj = helm2dquas(type, zk, kappa, d, coefs, quad_opts)
+function obj = helm2dquas(type, zk, kappa, d, coefs, quad_opts, ising)
 %KERNEL.HELM2DQUAS   Construct the quasi periodic Helmholtz kernel.
 %   KERNEL.HELM2DQUAS('s', ZK, KAPPA, D) or 
 %   KERNEL.HELM2DQUAS('single', ZK, KAPPA, D) constructs the
@@ -86,8 +86,8 @@
 % lattice sum parameters
 l=2; N = 40; a = 15; M = 1e4;
 
-if nargin == 6
-    if isefield(quad_opts,'l')
+if nargin >= 6
+    if isfield(quad_opts,'l')
         l = quad_opts.l;
     end
     if isfield(quad_opts,'N')
@@ -101,6 +101,10 @@
     end
 end
 
+if nargin < 7
+    ising = 1;
+end
+
 ns = (0:N).';
 sn = chnk.helm2dquas.latticecoefs(ns,zk,d,kappa,(exp(1i*kappa*d)),a,M,l+1);
 
@@ -110,18 +114,13 @@
 quas_param.l = l;
 quas_param.sn = sn;
 
-% obj.params.kappa = kappa;
-% obj.params.d = d;
-% obj.params.l = l;
-% obj.params.Sn = Sn;
-
 obj.params.quas_param = quas_param;
-
+obj.params.ising = ising;
 switch lower(type)
 
     case {'s', 'single'}
         obj.type = 's';
-        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 's',quas_param);
+        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 's',quas_param,[],ising);
         obj.fmm = [];
         obj.sing = 'log';
         if isscalar(kappa)
@@ -133,7 +132,7 @@
 
     case {'d', 'double'}
         obj.type = 'd';
-        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'd',quas_param);
+        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'd',quas_param,[],ising);
         obj.fmm = [];
         obj.sing = 'log';
         if isscalar(kappa)
@@ -145,13 +144,13 @@
 
     case {'sp', 'sprime'}
         obj.type = 'sp';
-        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'sprime',quas_param);
+        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'sprime',quas_param,[],ising);
         obj.fmm = [];
         obj.sing = 'log';
 
     case {'dp', 'dprime'}
         obj.type = 'dp';
-        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'dprime',quas_param);
+        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'dprime',quas_param,[],ising);
         obj.fmm = [];
         obj.sing = 'hs';
 
@@ -162,7 +161,7 @@
         end
         obj.type = 'c';
         obj.params.coefs = coefs;
-        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'c',quas_param, coefs);
+        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'c',quas_param,coefs,ising);
         obj.fmm = [];
         obj.sing = 'log';
         if isscalar(kappa)
@@ -179,7 +178,7 @@
         end
         obj.type = 'cprime';
         obj.params.coefs = coefs;
-        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'cprime',quas_param, coefs);
+        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'cprime',quas_param,coefs,ising);
         obj.fmm = [];
         obj.sing = 'hs';
 
@@ -231,7 +230,10 @@
         error('Unknown Helmholtz kernel type ''%s''.', type);
 
 end
-
+if obj.params.ising == 0
+    obj.sing = 'smooth'; 
+    obj.splitinfo = [];
+end
 end
 
 function f = helm2dquas_s_split(zk,s,t,quas_param)
diff --git a/devtools/test/kernelopTest.m b/devtools/test/kernelopTest.m
index 5a56fc1..019c246 100644
--- a/devtools/test/kernelopTest.m
+++ b/devtools/test/kernelopTest.m
@@ -1,5 +1,7 @@
 %KERNELOPTEST verify kernel class operations are correct
+kernelopTest0()
 
+function kernelopTest0()
 % setup points
 src = []; 
 src.r = [[0;0],[0;1]];
@@ -36,5 +38,8 @@
 assert(norm(ckern1.eval(src,targ) - (skern.eval(src,targ)+dkern.eval(src,targ)))<1e-10)
 assert(norm(ckern2.eval(src,targ) - (skern.eval(src,targ)-dkern.eval(src,targ)))<1e-10)
 
-
+% conjugate a kernel
+conj_dkern = conj(dkern);
+assert(norm(conj_dkern.eval(src,targ) - conj_dkern.eval(src,targ))<1e-10)
+end
 
diff --git a/devtools/test/quasiperiodicTest.m b/devtools/test/quasiperiodicTest.m
index e30fedd..cb4b019 100644
--- a/devtools/test/quasiperiodicTest.m
+++ b/devtools/test/quasiperiodicTest.m
@@ -7,12 +7,16 @@ function quasiperiodicTest0()
 % problem parameters
 d= 1;
 zk = 1;
-kappa = [pi-0.1-1i,0.2];
+kappa = [pi-0.1-1i,0.2,-0.1,0.25i];
 coefs = [1;0.5];
 coefa = [1,0.3;-1i,0.2];
+nkappa = length(kappa);
+
+src = []; src.r = [[0;-1.1],[1;-1],[0.1;-0.3]]; src.n = [[1;-2],[2;-1],[1;0]];
+targ = []; targ.r = [[1.1;0.3],[2;0]]; targ.n = [[-1;0.3],[0.1;1]];
+ns = size(src.r,2);
+nt = size(targ.r,2);
 
-src = []; src.r = [0;-1]; src.n = [1;-2];
-targ = []; targ.r = [1.1;0.3]; targ.n = [-1;0.3];
 
 % kernels
 skern = kernel('hq','s',zk,kappa,d);
@@ -60,33 +64,71 @@ function quasiperiodicTest0()
 % test quasiperiodicitiy
 targ_s = targ; targ_s.r = targ.r + [d;0];
 svalshift = skern.eval(src,targ_s);
-assert(norm(svalshift - exp(1i*kappa(:)*d).*sval) <  1e-12)
+assert(norm(svalshift - repmat(exp(1i*kappa(:)*d),nt,1).*sval) <  1e-12)
 
 % 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;
+all_assemb = zeros(2,nkappa*nt,2*ns);
+
+all_assemb(1,:, 1:2:end) = coefa(1,1)*dval;
+all_assemb(1,:, 2:2:end) = coefa(1,2)*sval;
+all_assemb(2,:, 1:2:end) = coefa(2,1)*dpval;
+all_assemb(2,:, 2:2:end) = coefa(2,2)*spval;
+all_assemb = reshape(all_assemb,[2,nkappa,nt,2*ns]);
+all_assemb = permute(all_assemb,[2,1,3,4]);
+all_assemb = reshape(all_assemb, size(allval) );
+
 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 transmission representation
+tr_assemb = zeros(nkappa*nt,2,ns);
+tr_assemb(:,1,:) = coefs(1)*dval;
+tr_assemb(:,2,:) = coefs(2)*sval;
+assert(norm(tr_assemb(:,:)-trval) <  1e-12)
+trp_assemb = zeros(nkappa*nt,2,ns);
+trp_assemb(:,1,:) = coefs(1)*dpval;
+trp_assemb(:,2,:) = coefs(2)*spval;
+assert(norm(trp_assemb(:,:)-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)
+trg_assemb = zeros(2*nkappa*nt,2,ns);
+trg_assemb(:,1,:) = coefs(1)*dgval;
+trg_assemb(:,2,:) = coefs(2)*sgval;
+assert(norm(trg_assemb(:,:)-trgval) <  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;
+c2trval_assemb = zeros(nkappa,2,nt,ns);
+c2trval_assemb(:,1,:,:) = reshape(coefs(1)*dval + coefs(2)*sval,nkappa,1,nt,ns);
+c2trval_assemb(:,2,:,:) = reshape(coefs(1)*dpval + coefs(2)*spval,nkappa,1,nt,ns);
+c2trval_assemb = reshape(c2trval_assemb,[],ns);
 assert(norm(c2trval_assemb-c2trval) <  1e-12)
+
+% test ising parameter
+skern1a = kernel('hq','s',zk,kappa,d,[],[],0);
+skern1b = kernel('h','s',zk);
+sval1b = skern1b.eval(src,targ); 
+sval1b = repmat(reshape(sval1b,1,nt,ns), nkappa,1,1);
+sval1 = skern1a.eval(src,targ) + reshape(sval1b,nt*nkappa,ns);
+assert(norm(sval-sval1) <  1e-12)
+
+dkern1a = kernel('hq','d',zk,kappa,d,[],[],0);
+dkern1b = kernel('h','d',zk);
+dval1b = dkern1b.eval(src,targ); 
+dval1b = repmat(reshape(dval1b,1,nt,ns), nkappa,1,1);
+dval1 = dkern1a.eval(src,targ) + reshape(dval1b,nt*nkappa,ns);
+assert(norm(dval-dval1) <  1e-12)
+
+allkern1a = kernel(@(s,t) chnk.helm2dquas.kern(zk,s,t,'all',quas_param,coefa,0));
+allkern1b = kernel(@(s,t) chnk.helm2d.kern(zk,s,t,'all',coefa));
+allval1a = allkern1a.eval(src,targ);
+allval1b = allkern1b.eval(src,targ); 
+allval1b = reshape(allval1b,1,2*nt,2*ns); allval1b = repmat(allval1b,length(kappa),1,1);
+allval1b = reshape(allval1b,length(kappa)*2*nt,2*ns);
+assert(norm(allval1a+allval1b -allval) <  1e-12)
 end
 
 
@@ -97,7 +139,7 @@ function quasiperiodicTest1()
 % problem parameters
 d= 8;
 zk = 1;
-kappa = .05+0.1i;
+kappa = .05+.1i;
 
 % setup geometry
 nch = 2^3;
diff --git a/devtools/test/tochunkgraphTest.m b/devtools/test/tochunkgraphTest.m
new file mode 100644
index 0000000..5cb31b7
--- /dev/null
+++ b/devtools/test/tochunkgraphTest.m
@@ -0,0 +1,62 @@
+tochunkgraphTest0();
+
+function tochunkgraphTest0()
+%TOCHUNKGRAPHTEST tests the routine for converting chunkers to chunkgraphs
+
+seed = 8675309;
+rng(seed);
+
+% geometry parameters and construction
+
+
+cparams = [];
+cparams.eps = 1.0e-9;
+pref = []; 
+pref.k = 16;
+narms = 0;
+amp = 0.0;
+% make a circle
+chnkr = chunkerfunc(@(t) starfish(t,narms,amp),cparams,pref);
+
+% make an open arc
+cparams.ifclosed = 0;
+chnkr2 = chunkerfunc(@(t) cos_func(t,pi,1),cparams,pref);
+chnkr2 = move(chnkr2,[0;3]);
+
+% merge arcs
+chnkrs(3) = chunker();
+chnkrs(1) = chnkr;
+rfac = 1.1;
+chnkrs(2) = move(chnkr, [0;0], [3;0], 0, rfac);
+chnkrs(3) = chnkr2;
+chnkrtotal = merge(chnkrs);
+
+% form chunkgraph
+cgrph = tochunkgraph(chnkrtotal);
+
+assert(size(cgrph.verts,2) == 4)
+assert(length(cgrph.echnks) == 3)
+assert(cgrph.npt == chnkrtotal.npt)
+assert(norm(cgrph.echnks(1).r(:) - chnkr2.r(:))<1e-14)
+
+verts = [[0;0],[1;1],[3;0]]; edge2verts = [[1;2],[3;3]];
+cgrph = chunkgraph(verts,edge2verts,{chnkr2,chnkr});
+assert(isequal(cgrph.verts,verts))
+assert(isequal(cgrph.edgesendverts,edge2verts))
+
+% verify chunker shifted and scaled correctly 
+vend = cgrph.echnks(1).r(:,1);
+assert(norm(vend-verts(:,1)) < 1e-2);
+
+vend = cgrph.echnks(1).r(:,end);
+assert(norm(vend-verts(:,2)) < 1e-2);
+
+end
+
+function [r,d,d2] = cos_func(t,d,A)
+% parameterization of sinusoidal boundary with period d and amplitude A
+omega = 2*pi/d;
+r = [t(:), A*cos(omega*t(:))].';
+d = [ones(length(t),1), -omega*A*sin(omega*t(:))].';
+d2 = [zeros(length(t),1), -omega^2*A*cos(omega*t(:))].';
+end
\ No newline at end of file