From a5dd01f387dc7f57e465e7e2840d913a0d25bd45 Mon Sep 17 00:00:00 2001
From: Tristan Goodwill <48371600+tristangdwl@users.noreply.github.com>
Date: Mon, 30 Jun 2025 11:39:29 -0500
Subject: [PATCH 01/12] update shape of vector valued quasiperiodic kernels

---
 chunkie/+chnk/+helm2d/kern.m      |  7 +--
 chunkie/+chnk/+helm2dquas/green.m | 22 +++-----
 chunkie/+chnk/+helm2dquas/kern.m  | 90 ++++++++++++++++---------------
 chunkie/@kernel/helm2dquas.m      |  8 +--
 devtools/test/quasiperiodicTest.m | 14 ++---
 5 files changed, 64 insertions(+), 77 deletions(-)

diff --git a/chunkie/+chnk/+helm2d/kern.m b/chunkie/+chnk/+helm2d/kern.m
index 0b535c4..401e0e4 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]
@@ -250,6 +250,7 @@
 if strcmpi(type,'c2trans') 
   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(:,:);
 
@@ -272,8 +273,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;
 
 end
 
diff --git a/chunkie/+chnk/+helm2dquas/green.m b/chunkie/+chnk/+helm2dquas/green.m
index 781b5ff..5393e11 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
+%CHNK.HELM2DQUAS.GREEN evaluate the quasiperiodic Helmholtz Green's function
 % for the given sources and targets
 %
 % Input:
@@ -13,7 +13,6 @@
 [~,nsrc] = size(src);
 [~,ntarg] = size(targ);
 
-
 xs = repmat(src(1,:),ntarg,1);
 ys = repmat(src(2,:),ntarg,1);
 
@@ -40,8 +39,7 @@
 
 npt = size(r,1);
 
-ythresh = d/2;
-% ythresh = d/10;
+ythresh = 2*d/2;
 iclose = abs(ry) < ythresh;
 ifar = ~iclose;
 
@@ -76,23 +74,19 @@
 % 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));
@@ -101,7 +95,9 @@
 grad_near = zeros(1,1,2);
 hess_near = zeros(1,1,3);
 if ~isempty(rxclose)
-    for i = -l:l
+    ls = [-l:-1, 1:l];
+    ls = -l:l;
+    for i = ls
         rxi = rxclose - i*d;
         if nargout>2
         [vali,gradi,hessi] = chnk.helm2d.green(zk,[0;0],[rxi.';ryclose.']);
@@ -123,8 +119,7 @@
         val_near = val_near + vali.*alpha.^i;
         end
     end
-    
-    % sn = sn.';
+
     N = size(sn,2)-1;
     ns = (0:N);
     ns_use = (0:N+2);
@@ -149,9 +144,6 @@
     val_far = 0.25*1i*Js(:,:,1).*sn(:,:,1) + 0.5*1i*sum(sn(:,:,2:end).*Js(:,:,2:end-2).*cs(:,:,2:end),3);
     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;
@@ -202,7 +194,5 @@
 if nargout>2
 hess = reshape(quasi_phase.*hess,nkappa*ntarg,nsrc,3);
 end
-% end
-
 end
 
diff --git a/chunkie/+chnk/+helm2dquas/kern.m b/chunkie/+chnk/+helm2dquas/kern.m
index ecb4568..bb5f387 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]
@@ -80,6 +80,7 @@
 src = srcinfo.r(:,:);
 targ = targinfo.r(:,:);
 
+
 kappa = quas_param.kappa;
 d = quas_param.d;
 l = quas_param.l;
@@ -111,13 +112,6 @@
   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')
   submat = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
@@ -145,21 +139,21 @@
 % gradient of single layer
 if strcmpi(type,'sgrad')
     [~,grad] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
-    submat = reshape(permute(grad,[3,1,2]),[],ns);
+    submat = reshape(permute(grad,[1,3,2]),[],ns);
 end
 
 % gradient of double layer
 if strcmpi(type,'dgrad')
     [~,~,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);
+    submat = reshape(permute(submat,[1,3,2]),[],ns);
 end
 
 % Combined field 
 if strcmpi(type,'c')
   srcnorm = srcinfo.n(:,:);
   coef = ones(2,1);
-  if(nargin == 6); coef = varargin{1}; end
+  if(nargin >= 6); coef = varargin{1}; end
   [submats,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);
@@ -170,7 +164,7 @@
 % normal derivative of combined field
 if strcmpi(type,'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(:,:);
 
@@ -196,8 +190,9 @@
 
 % Dirichlet and neumann data corresponding to combined field
 if strcmpi(type,'c2trans') 
-  coef = ones(2,1);
-  if(nargin == 6); coef = varargin{1}; end
+  coef = ones(1,2);
+  if(nargin >= 6); coef = varargin{1}; end
+  if numel(coef) == 2, coef = repmat(coef(:).',2,1); end
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
 
@@ -220,23 +215,26 @@
   % D
   submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc);
 
-  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 = 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 = 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);
 end
 
 % gradient of combined field
 if strcmpi(type,'cgrad')
     coef = ones(2,1);
-    if(nargin == 6); coef = varargin{1}; end
+    if(nargin >= 6); coef = varargin{1}; end
     [~,grad,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
 
-    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;
 end
@@ -269,26 +267,27 @@
   % 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') 
 
   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 = 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;
 end
@@ -296,7 +295,7 @@
 % Neumann data corresponding to transmission rep
 if strcmpi(type,'trans_rep_prime')
   coef = ones(2,1);
-  if(nargin == 6); coef = varargin{1}; end;
+  if(nargin >= 6); coef = varargin{1}; end;
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
   
@@ -305,10 +304,12 @@
   % Get gradient and hessian info
   [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);
-  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 ...
@@ -324,7 +325,7 @@
 % Gradient correpsonding to transmission rep
 if strcmpi(type,'trans_rep_grad')
   coef = ones(2,1);
-  if(nargin == 6); coef = varargin{1}; end;
+  if(nargin >= 6); coef = varargin{1}; end;
   
   srcnorm = srcinfo.n(:,:);
   
@@ -332,18 +333,19 @@
   % 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);
+  nxsrc = repmat(srcnorm(1,:),nkappa*nt,1);
+  nysrc = repmat(srcnorm(2,:),nkappa*nt,1);
   % D
   submatd  = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc);
 
-  submat = zeros(nkappa*2*nt,2*ns);
+  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) = -coef(1)*(hess(:,:,1).*nxsrc + hess(:,:,2).*nysrc);
+  submat(:,1,2:2:2*ns) = coef(2)*grad(:,:,1);
     
-  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) = -coef(1)*(hess(:,:,2).*nxsrc + hess(:,:,3).*nysrc);
+  submat(:,2,2:2:2*ns) = coef(2)*grad(:,:,2);
+  submat = reshape(submat,nkappa*2*nt,2*ns);
 end
 
 
diff --git a/chunkie/@kernel/helm2dquas.m b/chunkie/@kernel/helm2dquas.m
index d71db42..a469177 100644
--- a/chunkie/@kernel/helm2dquas.m
+++ b/chunkie/@kernel/helm2dquas.m
@@ -60,7 +60,7 @@
 l=2; N = 40; a = 15; M = 1e4;
 
 if nargin == 6
-    if isefield(quad_opts,'l')
+    if isfield(quad_opts,'l')
         l = quad_opts.l;
     end
     if isfield(quad_opts,'N')
@@ -83,11 +83,6 @@
 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;
 
 switch lower(type)
@@ -160,7 +155,6 @@
         error('Unknown Helmholtz kernel type ''%s''.', type);
 
 end
-
 end
 
 function f = helm2dquas_s_split(zk,s,t,quas_param)
diff --git a/devtools/test/quasiperiodicTest.m b/devtools/test/quasiperiodicTest.m
index e30fedd..30302aa 100644
--- a/devtools/test/quasiperiodicTest.m
+++ b/devtools/test/quasiperiodicTest.m
@@ -68,10 +68,10 @@ function quasiperiodicTest0()
 
 % 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(1:2, 1:2:end) = coefa(1,1)*dval;
+all_assemb(1:2, 2:2:end) = coefa(1,2)*sval;
+all_assemb(3:end, 1:2:end) = coefa(2,1)*dpval;
+all_assemb(3:end, 2:2:end) = coefa(2,2)*spval;
 assert(norm(all_assemb - allval) <  1e-12)
 
 % test transmission representiatoon
@@ -84,8 +84,8 @@ function quasiperiodicTest0()
 
 % 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(1:2, :) = coefs(1)*dval + coefs(2)*sval;
+c2trval_assemb(3:end, :) = coefs(1)*dpval + coefs(2)*spval;
 assert(norm(c2trval_assemb-c2trval) <  1e-12)
 end
 
@@ -97,7 +97,7 @@ function quasiperiodicTest1()
 % problem parameters
 d= 8;
 zk = 1;
-kappa = .05+0.1i;
+kappa = .05+.1i;
 
 % setup geometry
 nch = 2^3;

From 458dd2e1389414892ec4599cfe1bd7f166844293 Mon Sep 17 00:00:00 2001
From: Tristan Goodwill <48371600+tristangdwl@users.noreply.github.com>
Date: Mon, 30 Jun 2025 11:41:05 -0500
Subject: [PATCH 02/12] update quasi green ls

---
 chunkie/+chnk/+helm2dquas/green.m | 1 -
 1 file changed, 1 deletion(-)

diff --git a/chunkie/+chnk/+helm2dquas/green.m b/chunkie/+chnk/+helm2dquas/green.m
index 5393e11..2aaca02 100644
--- a/chunkie/+chnk/+helm2dquas/green.m
+++ b/chunkie/+chnk/+helm2dquas/green.m
@@ -95,7 +95,6 @@
 grad_near = zeros(1,1,2);
 hess_near = zeros(1,1,3);
 if ~isempty(rxclose)
-    ls = [-l:-1, 1:l];
     ls = -l:l;
     for i = ls
         rxi = rxclose - i*d;

From db745cf4f1806407dd619f45528f629b47b790e5 Mon Sep 17 00:00:00 2001
From: Tristan Goodwill <48371600+tristangdwl@users.noreply.github.com>
Date: Fri, 11 Jul 2025 14:35:05 -0400
Subject: [PATCH 03/12] add tochunkgraph

---
 chunkie/@chunker/chunker.m       |  13 +---
 chunkie/@chunker/tochunkgraph.m  |  50 ++++++++++++++
 chunkie/@chunkgraph/chunkgraph.m | 114 ++++++++++++++++++-------------
 devtools/test/tochunkgraphTest.m |  62 +++++++++++++++++
 4 files changed, 181 insertions(+), 58 deletions(-)
 create mode 100644 chunkie/@chunker/tochunkgraph.m
 create mode 100644 devtools/test/tochunkgraphTest.m

diff --git a/chunkie/@chunker/chunker.m b/chunkie/@chunker/chunker.m
index 3e3aa6c..466388f 100644
--- a/chunkie/@chunker/chunker.m
+++ b/chunkie/@chunker/chunker.m
@@ -29,12 +29,9 @@
 % chunker methods: 
 %   chunker(p,t,w) - construct an empty chunker with given preferences and
 %       precomputed Legendre nodes/weights (optional)
-%   obj = plus(v,obj) - provides translation via v + obj
-%   obj = mtimes(A,obj) - provides affine transformation via A*obj
-%   obj = rotate(obj,theta,r0,r1) - rotate by angle 
-%   obj = reflect(obj,theta,r0,r1) - reflect across line
 %   obj = addchunk(obj,nchadd) - add nchadd chunks to the structure 
 %       (initialized with zeros)
+%   obj = obj.move(r0,r1,trotat,scale) - translate, rotate, etc
 %   obj = makedatarows(obj,nrows) - add nrows rows to the data storage.
 %   [obj,info] = sort(obj) - sort the chunks so that adjacent chunks are
 %        stored sequentially
@@ -56,6 +53,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
@@ -68,7 +66,6 @@
 %              number of connected components, etc
 %   [re,taue] = chunkends(obj,ich) - get the endpoints of chunks
 %   flag = flagnear(obj,pts,opts) - flag points near the boundary
-%   obj = obj.move(r0,r1,trotat,scale) - translate, rotate, etc
 
 
 % author: Travis Askham (askhamwhat@gmail.com)
@@ -363,6 +360,7 @@
         quiver(obj,varargin)
         scatter(obj,varargin)
         tau = taus(obj)
+        cgrph = tochunkgraph(obj)
         obj = refine(obj,varargin)
         a = area(obj)
         s = arclength(obj)
@@ -373,11 +371,6 @@
         [re,taue] = chunkends(obj,ich)
         flag = flagnear(obj,pts,opts)
         kappa = signed_curvature(obj)
-        obj = plus(v,obj)
-        obj = mtimes(A,obj)
-        obj = rotate(obj,theta,r0,r1)
-        obj = reflect(obj,theta,r0,r1)
-        du = arclengthder(obj,u)
     end
     methods(Static)
         obj = chunkerfunc(fcurve,varargin)
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..b020b4b 100644
--- a/chunkie/@chunkgraph/chunkgraph.m
+++ b/chunkie/@chunkgraph/chunkgraph.m
@@ -50,10 +50,6 @@
 %         meets with other chunk ends in a vertex. 
 %
 % chunkgraph methods:
-%   obj = plus(v,obj) - provides translation via v + obj
-%   obj = mtimes(A,obj) - provides affine transformation via A*obj
-%   obj = rotate(obj,theta,r0,r1) - rotate by angle 
-%   obj = reflect(obj,theta,r0,r1) - reflect across line
 %   plot(obj, varargin) - plot the chunkgraph
 %   quiver(obj, varargin) - quiver plot of chunkgraph points and normals
 %   plot_regions(obj, iflabel) - plot the chunkgraph with region and 
@@ -94,14 +90,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 +175,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 +238,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 +246,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 +338,6 @@
                             chnkr = move(chnkr,zeros(size(vs(:,1))),vfin-vs(:,1),0,1);
                         end
                     end
-
                     echnks(i) = chnkr;
                 end
             end
@@ -433,10 +453,8 @@
         scatter(obj, varargin)
         rres = datares(obj, opts)
         edge_regs = find_edge_regions(obj)
-        obj = plus(v,obj)
-        obj = mtimes(A,obj)        
-        obj = rotate(obj,theta,r0,r1)
-        obj = reflect(obj,theta,r0,r1)   
+
+        
     end
 
     methods(Static)
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

From 8ee0e79128c0698a1f9b76c81d6f26454f64c83a Mon Sep 17 00:00:00 2001
From: Tristan Goodwill <48371600+tristangdwl@users.noreply.github.com>
Date: Fri, 11 Jul 2025 14:53:41 -0400
Subject: [PATCH 04/12] repair documentation

---
 chunkie/@chunker/chunker.m       | 10 +++++++++-
 chunkie/@chunkgraph/chunkgraph.m | 11 ++++++++---
 2 files changed, 17 insertions(+), 4 deletions(-)

diff --git a/chunkie/@chunker/chunker.m b/chunkie/@chunker/chunker.m
index 466388f..0a644cb 100644
--- a/chunkie/@chunker/chunker.m
+++ b/chunkie/@chunker/chunker.m
@@ -29,9 +29,12 @@
 % chunker methods: 
 %   chunker(p,t,w) - construct an empty chunker with given preferences and
 %       precomputed Legendre nodes/weights (optional)
+%   obj = plus(v,obj) - provides translation via v + obj
+%   obj = mtimes(A,obj) - provides affine transformation via A*obj
+%   obj = rotate(obj,theta,r0,r1) - rotate by angle 
+%   obj = reflect(obj,theta,r0,r1) - reflect across line
 %   obj = addchunk(obj,nchadd) - add nchadd chunks to the structure 
 %       (initialized with zeros)
-%   obj = obj.move(r0,r1,trotat,scale) - translate, rotate, etc
 %   obj = makedatarows(obj,nrows) - add nrows rows to the data storage.
 %   [obj,info] = sort(obj) - sort the chunks so that adjacent chunks are
 %        stored sequentially
@@ -66,6 +69,7 @@
 %              number of connected components, etc
 %   [re,taue] = chunkends(obj,ich) - get the endpoints of chunks
 %   flag = flagnear(obj,pts,opts) - flag points near the boundary
+%   obj = obj.move(r0,r1,trotat,scale) - translate, rotate, etc
 
 
 % author: Travis Askham (askhamwhat@gmail.com)
@@ -371,6 +375,10 @@
         [re,taue] = chunkends(obj,ich)
         flag = flagnear(obj,pts,opts)
         kappa = signed_curvature(obj)
+        obj = plus(v,obj)
+        obj = mtimes(A,obj)
+        obj = rotate(obj,theta,r0,r1)
+        obj = reflect(obj,theta,r0,r1)
     end
     methods(Static)
         obj = chunkerfunc(fcurve,varargin)
diff --git a/chunkie/@chunkgraph/chunkgraph.m b/chunkie/@chunkgraph/chunkgraph.m
index b020b4b..d719bd1 100644
--- a/chunkie/@chunkgraph/chunkgraph.m
+++ b/chunkie/@chunkgraph/chunkgraph.m
@@ -50,6 +50,10 @@
 %         meets with other chunk ends in a vertex. 
 %
 % chunkgraph methods:
+%   obj = plus(v,obj) - provides translation via v + obj
+%   obj = mtimes(A,obj) - provides affine transformation via A*obj
+%   obj = rotate(obj,theta,r0,r1) - rotate by angle 
+%   obj = reflect(obj,theta,r0,r1) - reflect across line
 %   plot(obj, varargin) - plot the chunkgraph
 %   quiver(obj, varargin) - quiver plot of chunkgraph points and normals
 %   plot_regions(obj, iflabel) - plot the chunkgraph with region and 
@@ -63,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
@@ -453,8 +456,10 @@
         scatter(obj, varargin)
         rres = datares(obj, opts)
         edge_regs = find_edge_regions(obj)
-
-        
+        obj = plus(v,obj)
+        obj = mtimes(A,obj)        
+        obj = rotate(obj,theta,r0,r1)
+        obj = reflect(obj,theta,r0,r1)   
     end
 
     methods(Static)

From d56a254622ad0338b2dd5451d8dec93d7fd54d4e Mon Sep 17 00:00:00 2001
From: Tristan Goodwill <48371600+tristangdwl@users.noreply.github.com>
Date: Tue, 15 Jul 2025 21:30:39 -0400
Subject: [PATCH 05/12] add complex conjugates of kernels

---
 chunkie/@kernel/conj.m       | 22 ++++++++++++++++++++++
 devtools/test/kernelopTest.m |  7 ++++++-
 2 files changed, 28 insertions(+), 1 deletion(-)
 create mode 100644 chunkie/@kernel/conj.m

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/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
 

From a9ac55aabe8a77e8d24bc81a58beb170cba2ee10 Mon Sep 17 00:00:00 2001
From: Tristan Goodwill <48371600+tristangdwl@users.noreply.github.com>
Date: Fri, 12 Sep 2025 14:28:19 -0500
Subject: [PATCH 06/12] add ising option to helm2dquas.green

---
 chunkie/+chnk/+helm2dquas/green.m | 92 ++++++++++++++++++++++++-------
 chunkie/+chnk/+helm2dquas/kern.m  | 35 +++++++-----
 chunkie/@kernel/helm2dquas.m      | 26 +++++----
 3 files changed, 109 insertions(+), 44 deletions(-)

diff --git a/chunkie/+chnk/+helm2dquas/green.m b/chunkie/+chnk/+helm2dquas/green.m
index 2aaca02..215c3d9 100644
--- a/chunkie/+chnk/+helm2dquas/green.m
+++ b/chunkie/+chnk/+helm2dquas/green.m
@@ -1,4 +1,4 @@
-function [val,grad,hess] = green(src,targ,zk,kappa,d,sn,l)
+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
 %
@@ -8,6 +8,9 @@
 %   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);
@@ -50,6 +53,9 @@
 ryclose = ry(iclose);
 rclose = r(iclose);
 
+nxclose = nx(iclose);
+nptclose = size(rxclose, 1);
+
 nkappa = length(kappa);
 
 val = zeros(nkappa,npt,1);
@@ -91,31 +97,37 @@
 
 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);
 if ~isempty(rxclose)
     ls = -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
 
@@ -179,19 +191,61 @@
 
         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
 
+if nargout == 1
+    val = quasi_phase.*val;
+
+    if ising == 0
+        isub = (abs(nx) > max(ls)) | ifar;
+
+        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
+
+    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;
+    
+        [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
+
+    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;
+
+        [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
+
+    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 bb5f387..26b8a06 100644
--- a/chunkie/+chnk/+helm2dquas/kern.m
+++ b/chunkie/+chnk/+helm2dquas/kern.m
@@ -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
@@ -81,6 +83,11 @@
 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 +101,7 @@
 % double layer
 if strcmpi(type,'d')
   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);
@@ -103,7 +110,7 @@
 % normal derivative of single layer
 if strcmpi(type,'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);
@@ -114,7 +121,7 @@
 
 % single later
 if strcmpi(type,'s')
-  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);
 end
 
@@ -122,7 +129,7 @@
 if strcmpi(type,'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);
@@ -138,13 +145,13 @@
 
 % gradient of single layer
 if strcmpi(type,'sgrad')
-    [~,grad] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
+    [~,grad] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l,ising);
     submat = reshape(permute(grad,[1,3,2]),[],ns);
 end
 
 % gradient of double layer
 if strcmpi(type,'dgrad')
-    [~,~,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,[1,3,2]),[],ns);
 end
@@ -154,7 +161,7 @@
   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);
+  [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);
@@ -169,7 +176,7 @@
   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);
@@ -197,7 +204,7 @@
   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);
@@ -229,7 +236,7 @@
 if strcmpi(type,'cgrad')
     coef = ones(2,1);
     if(nargin >= 6); coef = varargin{1}; end
-    [~,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);
 
     submats = reshape(permute(grad,[1,3,2]),[],ns);
 
@@ -248,7 +255,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);
@@ -280,7 +287,7 @@
   coef = ones(2,1);
   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);
   nxsrc = repmat(srcnorm(1,:),nkappa*nt,1);
@@ -302,7 +309,7 @@
   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,:),nkappa*nt,1);
   nysrc = repmat(srcnorm(2,:),nkappa*nt,1);
@@ -331,7 +338,7 @@
   
   % submat = zeros(nt,ns,6);
   % S
-  [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);
diff --git a/chunkie/@kernel/helm2dquas.m b/chunkie/@kernel/helm2dquas.m
index a469177..9747bcb 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
@@ -59,7 +59,7 @@
 % lattice sum parameters
 l=2; N = 40; a = 15; M = 1e4;
 
-if nargin == 6
+if nargin >= 6
     if isfield(quad_opts,'l')
         l = quad_opts.l;
     end
@@ -74,6 +74,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);
 
@@ -89,37 +93,37 @@
 
     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)
         obj.splitinfo = [];
         obj.splitinfo.type = {[0 0 0 0],[1 0 0 0]};
         obj.splitinfo.action = {'r','r'};
-        obj.splitinfo.functions = @(s,t) helm2dquas_s_split(zk,s,t,quas_param);
+        obj.splitinfo.functions = @(s,t) helm2dquas_s_split(zk,s,t,quas_param,ising);
         end
 
     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)
         obj.splitinfo = [];
         obj.splitinfo.type = {[0 0 0 0],[1 0 0 0],[0 0 -1 0]};
         obj.splitinfo.action = {'r','r','r'};
-        obj.splitinfo.functions = @(s,t) helm2dquas_d_split(zk,s,t,quas_param);
+        obj.splitinfo.functions = @(s,t) helm2dquas_d_split(zk,s,t,quas_param,ising);
         end
 
     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';
 
@@ -130,14 +134,14 @@
         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)
         obj.splitinfo = [];
         obj.splitinfo.type = {[0 0 0 0],[1 0 0 0],[0 0 -1 0]};
         obj.splitinfo.action = {'r','r','r'};
-        obj.splitinfo.functions = @(s,t) helm2dquas_c_split(zk,s,t,quas_param,coefs);
+        obj.splitinfo.functions = @(s,t) helm2dquas_c_split(zk,s,t,quas_param,coefs,ising);
         end
 
     case {'cp', 'cprime'}
@@ -147,7 +151,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';
 

From d1495bbb2fd38271050bfb407980ce0d9958e6fa Mon Sep 17 00:00:00 2001
From: Tristan Goodwill <48371600+tristangdwl@users.noreply.github.com>
Date: Fri, 12 Sep 2025 20:18:04 -0500
Subject: [PATCH 07/12] test ising

---
 chunkie/+chnk/+helm2dquas/green.m |  8 +++++++-
 chunkie/@kernel/helm2dquas.m      | 20 ++++++++++++--------
 devtools/test/quasiperiodicTest.m | 19 +++++++++++++++++++
 3 files changed, 38 insertions(+), 9 deletions(-)

diff --git a/chunkie/+chnk/+helm2dquas/green.m b/chunkie/+chnk/+helm2dquas/green.m
index 215c3d9..24c457f 100644
--- a/chunkie/+chnk/+helm2dquas/green.m
+++ b/chunkie/+chnk/+helm2dquas/green.m
@@ -100,8 +100,8 @@
 val_near= zeros(nkappa,nptclose);
 grad_near = zeros(nkappa,nptclose,2);
 hess_near = zeros(nkappa,nptclose,3);
+ls = -l:l;
 if ~isempty(rxclose)
-    ls = -l:l;
     for i = ls
         if ising == 1
             iuse = true(nptclose,1);
@@ -204,9 +204,11 @@
     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);
@@ -217,11 +219,13 @@
     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);
@@ -234,6 +238,7 @@
     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);
@@ -242,6 +247,7 @@
         val(:,isub,:) = val(:,isub,:) - vali;
         grad(:,isub,:) = grad(:,isub,:) - gradi;
         hess(:,isub,:) = hess(:,isub,:) - hessi;
+        end
     end
 
     val = reshape(val,nkappa*ntarg,nsrc);
diff --git a/chunkie/@kernel/helm2dquas.m b/chunkie/@kernel/helm2dquas.m
index 9747bcb..41ff3ed 100644
--- a/chunkie/@kernel/helm2dquas.m
+++ b/chunkie/@kernel/helm2dquas.m
@@ -88,42 +88,42 @@
 quas_param.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,ising);
+        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 's',quas_param,[],ising);
         obj.fmm = [];
         obj.sing = 'log';
         if isscalar(kappa)
         obj.splitinfo = [];
         obj.splitinfo.type = {[0 0 0 0],[1 0 0 0]};
         obj.splitinfo.action = {'r','r'};
-        obj.splitinfo.functions = @(s,t) helm2dquas_s_split(zk,s,t,quas_param,ising);
+        obj.splitinfo.functions = @(s,t) helm2dquas_s_split(zk,s,t,quas_param);
         end
 
     case {'d', 'double'}
         obj.type = 'd';
-        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'd',quas_param,ising);
+        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'd',quas_param,[],ising);
         obj.fmm = [];
         obj.sing = 'log';
         if isscalar(kappa)
         obj.splitinfo = [];
         obj.splitinfo.type = {[0 0 0 0],[1 0 0 0],[0 0 -1 0]};
         obj.splitinfo.action = {'r','r','r'};
-        obj.splitinfo.functions = @(s,t) helm2dquas_d_split(zk,s,t,quas_param,ising);
+        obj.splitinfo.functions = @(s,t) helm2dquas_d_split(zk,s,t,quas_param);
         end
 
     case {'sp', 'sprime'}
         obj.type = 'sp';
-        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'sprime',quas_param,ising);
+        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,ising);
+        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'dprime',quas_param,[],ising);
         obj.fmm = [];
         obj.sing = 'hs';
 
@@ -141,7 +141,7 @@
         obj.splitinfo = [];
         obj.splitinfo.type = {[0 0 0 0],[1 0 0 0],[0 0 -1 0]};
         obj.splitinfo.action = {'r','r','r'};
-        obj.splitinfo.functions = @(s,t) helm2dquas_c_split(zk,s,t,quas_param,coefs,ising);
+        obj.splitinfo.functions = @(s,t) helm2dquas_c_split(zk,s,t,quas_param,coefs);
         end
 
     case {'cp', 'cprime'}
@@ -159,6 +159,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/quasiperiodicTest.m b/devtools/test/quasiperiodicTest.m
index 30302aa..7685cdb 100644
--- a/devtools/test/quasiperiodicTest.m
+++ b/devtools/test/quasiperiodicTest.m
@@ -87,6 +87,25 @@ function quasiperiodicTest0()
 c2trval_assemb(1:2, :) = coefs(1)*dval + coefs(2)*sval;
 c2trval_assemb(3:end, :) = coefs(1)*dpval + coefs(2)*spval;
 assert(norm(c2trval_assemb-c2trval) <  1e-12)
+
+
+skern1a = kernel('hq','s',zk,kappa,d,[],[],0);
+skern1b = kernel('h','s',zk);
+sval1 = skern1a.eval(src,targ) + skern1b.eval(src,targ);
+assert(norm(sval-sval1) <  1e-12)
+
+dkern1a = kernel('hq','d',zk,kappa,d,[],[],0);
+dkern1b = kernel('h','d',zk);
+dval1 = dkern1a.eval(src,targ) + dkern1b.eval(src,targ);
+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,2); allval1b = repmat(allval1b,length(kappa),1,1);
+allval1b = reshape(allval1b,length(kappa)*2,2);
+assert(norm(allval1a+allval1b -allval) <  1e-12)
 end
 
 

From 1184bd04e18285878c7f2a08c6f26ada6eef0abd Mon Sep 17 00:00:00 2001
From: Tristan Goodwill <48371600+tristangdwl@users.noreply.github.com>
Date: Tue, 16 Sep 2025 09:05:51 -0500
Subject: [PATCH 08/12] update quasi periodic green

---
 chunkie/+chnk/+helm2dquas/green.m | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/chunkie/+chnk/+helm2dquas/green.m b/chunkie/+chnk/+helm2dquas/green.m
index 24c457f..d5e9c6f 100644
--- a/chunkie/+chnk/+helm2dquas/green.m
+++ b/chunkie/+chnk/+helm2dquas/green.m
@@ -202,7 +202,7 @@
     val = quasi_phase.*val;
 
     if ising == 0
-        isub = (abs(nx) > max(ls)) | ifar;
+        isub = (abs(nx(:)) > max(ls)) | ifar;
 
         if any(isub)
         vali = chnk.helm2d.green(zk,[0;0],[rx(isub).' + nx(isub).'*d;ry(isub).']);
@@ -217,7 +217,7 @@
     grad = quasi_phase.*grad;
     
     if ising == 0
-        isub = (abs(nx) > max(ls)) | ifar;
+        isub = (abs(nx(:)) > max(ls)) | ifar;
     
         if any(isub)
         [vali, gradi] = chnk.helm2d.green(zk,[0;0],[rx(isub).' + nx(isub).'*d;ry(isub).']);
@@ -236,7 +236,7 @@
     hess = quasi_phase.*hess;
     
     if ising == 0
-        isub = (abs(nx) > max(ls)) | ifar;
+        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).']);

From ed213942ff601da264f8cd1fdd7ed75326e7af3a Mon Sep 17 00:00:00 2001
From: Travis Askham <askhamwhat@gmail.com>
Date: Wed, 1 Oct 2025 10:21:05 -0400
Subject: [PATCH 09/12] Update chunker.m

restoring arclengthder method
---
 chunkie/@chunker/chunker.m | 1 +
 1 file changed, 1 insertion(+)

diff --git a/chunkie/@chunker/chunker.m b/chunkie/@chunker/chunker.m
index 0a644cb..9ce6fbe 100644
--- a/chunkie/@chunker/chunker.m
+++ b/chunkie/@chunker/chunker.m
@@ -379,6 +379,7 @@
         obj = mtimes(A,obj)
         obj = rotate(obj,theta,r0,r1)
         obj = reflect(obj,theta,r0,r1)
+        du = arclengthder(obj,u)
     end
     methods(Static)
         obj = chunkerfunc(fcurve,varargin)

From cd844610e4a97d5a42354bbca3f7d705178e7bd2 Mon Sep 17 00:00:00 2001
From: Tristan Goodwill <48371600+tristangdwl@users.noreply.github.com>
Date: Fri, 3 Oct 2025 16:27:57 -0500
Subject: [PATCH 10/12] fix ising bug

---
 chunkie/+chnk/+helm2dquas/green.m | 11 +++---
 chunkie/+chnk/+helm2dquas/kern.m  |  8 +---
 devtools/test/quasiperiodicTest.m | 65 +++++++++++++++++++++----------
 3 files changed, 52 insertions(+), 32 deletions(-)

diff --git a/chunkie/+chnk/+helm2dquas/green.m b/chunkie/+chnk/+helm2dquas/green.m
index d5e9c6f..acbb083 100644
--- a/chunkie/+chnk/+helm2dquas/green.m
+++ b/chunkie/+chnk/+helm2dquas/green.m
@@ -25,6 +25,11 @@
 rx = xt-xs;
 ry = yt-ys;
 
+
+rx = rx(:);
+ry = ry(:);
+
+
 nx = fix(rx/d);
 rx = rx - nx*d;
 
@@ -36,10 +41,6 @@
 
 r = sqrt(r2);
 
-rx = rx(:);
-ry = ry(:);
-r = r(:);
-
 npt = size(r,1);
 
 ythresh = 2*d/2;
@@ -205,7 +206,7 @@
         isub = (abs(nx(:)) > max(ls)) | ifar;
 
         if any(isub)
-        vali = chnk.helm2d.green(zk,[0;0],[rx(isub).' + nx(isub).'*d;ry(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
diff --git a/chunkie/+chnk/+helm2dquas/kern.m b/chunkie/+chnk/+helm2dquas/kern.m
index 746b249..3ee7a63 100644
--- a/chunkie/+chnk/+helm2dquas/kern.m
+++ b/chunkie/+chnk/+helm2dquas/kern.m
@@ -337,15 +337,11 @@
   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,ising);
 
   nxsrc = repmat(srcnorm(1,:),nkappa*nt,1);
   nysrc = repmat(srcnorm(2,:),nkappa*nt,1);
-  % D
-  % submatd  = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc);
 
   submat = zeros(nkappa,2,nt,2*ns);
   
@@ -354,7 +350,7 @@
     
   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(submat,nkappa*2*nt,2*ns);
+  submat = reshape(permute(submat,[1,3,2,4]),nkappa*2*nt,2*ns);
   
 
 otherwise
diff --git a/devtools/test/quasiperiodicTest.m b/devtools/test/quasiperiodicTest.m
index 7685cdb..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,51 +64,70 @@ 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, 1:2:end) = coefa(1,1)*dval;
-all_assemb(1:2, 2:2:end) = coefa(1,2)*sval;
-all_assemb(3:end, 1:2:end) = coefa(2,1)*dpval;
-all_assemb(3: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, :) = coefs(1)*dval + coefs(2)*sval;
-c2trval_assemb(3: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);
-sval1 = skern1a.eval(src,targ) + skern1b.eval(src,targ);
+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);
-dval1 = dkern1a.eval(src,targ) + dkern1b.eval(src,targ);
+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,2); allval1b = repmat(allval1b,length(kappa),1,1);
-allval1b = reshape(allval1b,length(kappa)*2,2);
+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
 

From 499bcf6d802be2dd8b3cae1e8240128236fef2c7 Mon Sep 17 00:00:00 2001
From: Tristan Goodwill <48371600+tristangdwl@users.noreply.github.com>
Date: Fri, 3 Oct 2025 16:50:39 -0500
Subject: [PATCH 11/12] update matvec shape

---
 chunkie/+chnk/+helm2dquas/green.m | 64 +++++++++++++++++++++++++++----
 1 file changed, 57 insertions(+), 7 deletions(-)

diff --git a/chunkie/+chnk/+helm2dquas/green.m b/chunkie/+chnk/+helm2dquas/green.m
index d5e9c6f..73f9c5d 100644
--- a/chunkie/+chnk/+helm2dquas/green.m
+++ b/chunkie/+chnk/+helm2dquas/green.m
@@ -145,6 +145,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;
@@ -152,15 +153,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; 
@@ -175,24 +179,70 @@
         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
+
+
+    % val_far = 0.25*1i*Js(:,:,1).*sn(:,:,1) + 0.5*1i*sum(sn(:,:,2:end).*Js(:,:,2:end-2).*cs(:,:,2:end),3);
+    % 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.';
+    % 
+    %     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; 
+    % end
+    % if nargout > 2
+    %     DDJs = cat(3,.5*(Js(:,:,3)-Js(:,:,1)),.25*(Js(:,:,4)-3*Js(:,:,2)),.25*(Js(:,:,1:end-4)-2*Js(:,:,3:end-2)+Js(:,:,5:end)))*zk^2;
+    %     rclose = rclose.';
+    %     rxclose = rxclose.';
+    %     ryclose = ryclose.';
+    %     ns = reshape(ns,1,1,[]);
+    % 
+    %     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 = 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 = 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);
+    % 
+    %     hess_far = cat(3,hess_far_xx, hess_far_xy, hess_far_yy);
+    % 
+    %     hess(:,iclose,:) = hess_near + hess_far;
+    % end
+    % toc(t1);
 end
 
 quasi_phase = exp(1i*kappa(:)*nx(:).'*d);

From 0e5ec8cbada2a3594b1884b1967e1507a83b9d43 Mon Sep 17 00:00:00 2001
From: Tristan Goodwill <48371600+tristangdwl@users.noreply.github.com>
Date: Fri, 3 Oct 2025 17:07:20 -0500
Subject: [PATCH 12/12] tidy up

---
 chunkie/+chnk/+helm2dquas/green.m | 44 -------------------------------
 1 file changed, 44 deletions(-)

diff --git a/chunkie/+chnk/+helm2dquas/green.m b/chunkie/+chnk/+helm2dquas/green.m
index 78c0a30..43e6770 100644
--- a/chunkie/+chnk/+helm2dquas/green.m
+++ b/chunkie/+chnk/+helm2dquas/green.m
@@ -200,50 +200,6 @@
 
         hess(:,iclose,:) = hess_near + hess_far;
     end
-
-
-    % val_far = 0.25*1i*Js(:,:,1).*sn(:,:,1) + 0.5*1i*sum(sn(:,:,2:end).*Js(:,:,2:end-2).*cs(:,:,2:end),3);
-    % 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.';
-    % 
-    %     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; 
-    % end
-    % if nargout > 2
-    %     DDJs = cat(3,.5*(Js(:,:,3)-Js(:,:,1)),.25*(Js(:,:,4)-3*Js(:,:,2)),.25*(Js(:,:,1:end-4)-2*Js(:,:,3:end-2)+Js(:,:,5:end)))*zk^2;
-    %     rclose = rclose.';
-    %     rxclose = rxclose.';
-    %     ryclose = ryclose.';
-    %     ns = reshape(ns,1,1,[]);
-    % 
-    %     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 = 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 = 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);
-    % 
-    %     hess_far = cat(3,hess_far_xx, hess_far_xy, hess_far_yy);
-    % 
-    %     hess(:,iclose,:) = hess_near + hess_far;
-    % end
-    % toc(t1);
 end
 
 quasi_phase = exp(1i*kappa(:)*nx(:).'*d);