diff --git a/chunkie/+chnk/+helm1d/green.m b/chunkie/+chnk/+helm1d/green.m
new file mode 100644
index 0000000..e1fdc76
--- /dev/null
+++ b/chunkie/+chnk/+helm1d/green.m
@@ -0,0 +1,49 @@
+function [val,grad,hess] = green(zkE,src,targ)
+%CHNK.HELM1D.GREEN evaluate the 1D helmholtz green's function
+% for the given sources and targets
+
+[~,ns] = size(src);
+[~,nt] = size(targ);
+
+xs = repmat(src(1,:),nt,1);
+ys = repmat(src(2,:),nt,1);
+
+xt = repmat(targ(1,:).',1,ns);
+yt = repmat(targ(2,:).',1,ns);
+
+rx = xt-xs;
+ry = yt-ys;
+
+rx2 = rx.*rx;
+ry2 = ry.*ry;
+
+r2 = rx2+ry2;
+
+r = sqrt(r2);
+
+
+if nargout > 0
+    val = exp(1j*zkE*r);
+end
+
+[m,n] = size(xs);
+
+if nargout > 1
+    grad = zeros(m,n,2);
+    
+    grad(:,:,1) = (1/2).*(rx./r).*exp(1j*zkE*r);
+    grad(:,:,2) = (1/2).*(ry./r).*exp(1j*zkE*r);
+    grad = 2*1j*zkE*grad;
+end
+
+if nargout > 2
+
+    hess = zeros(m,n,3);
+
+    r3 = r.^3;
+    
+    hess(:,:,1) = (1./(2*r3)).*(1j*zkE*r.*rx2 - rx2 + r2).*exp(1j*zkE*r);
+    hess(:,:,2) = (1./(2*r3)).*(1j*zkE*r.*rx.*ry - rx.*ry).*exp(1j*zkE*r);
+    hess(:,:,3) = (1./(2*r3)).*(1j*zkE*r.*ry2 - ry2 + r2).*exp(1j*zkE*r);
+    hess = 2*1j*zkE*hess;
+end
diff --git a/chunkie/+chnk/+helm1d/kern.m b/chunkie/+chnk/+helm1d/kern.m
new file mode 100644
index 0000000..39dbbf9
--- /dev/null
+++ b/chunkie/+chnk/+helm1d/kern.m
@@ -0,0 +1,295 @@
+function submat= kern(zkE,srcinfo,targinfo,type,varargin)
+%CHNK.HELM1D.KERN standard Helmholtz layer potential kernels in 1D
+% 
+% Syntax: submat = chnk.heml1d.kern(zkE,srcinfo,targingo,type,varargin)
+%
+% Let x be targets and y be sources for these formulas, with
+% n_x and n_y the corresponding unit normals at those points
+% (if defined). Note that the normal information is obtained
+% by taking the perpendicular to the provided tangential deriviative
+% info and normalizing  
+%  
+% Kernels based on G(x,y) = e^{iE|x-y|}
+%
+% D(x,y) = \nabla_{n_y} G(x,y)
+% S(x,y) = G(x,y)
+% S'(x,y) = \nabla_{n_x} G(x,y)
+% D'(x,y) = \nabla_{n_x} \nabla_{n_y} G(x,y)
+%
+% Input:
+%   zkE - complex number, Helmholtz wave number
+%   srcinfo - description of sources in ptinfo struct format, i.e.
+%                ptinfo.r - positions (2,:) array
+%                ptinfo.d - first derivative in underlying
+%                     parameterization (2,:)
+%                ptinfo.d2 - second derivative in underlying
+%                     parameterization (2,:)
+%   targinfo - description of targets in ptinfo struct format,
+%                if info not relevant (d/d2) it doesn't need to
+%                be provided. sprime requires tangent info in
+%                targinfo.d
+%   type - string, determines kernel type
+%                type == 'd', double layer kernel D
+%                type == 's', single layer kernel S
+%                type == 'sprime', normal derivative of single
+%                      layer S'
+%                type == 'dprime', normal derivative of double layer D'
+%                type == 'c', combined layer kernel coef(1) D + coef(2) S
+%                type == 'stau', tangential derivative of single layer
+%                type == 'all', returns all four layer potentials, 
+%                       [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']
+%                type == 'trans_rep' returns the potential corresponding
+%                           to the transmission representation
+%                        [coef(1)*D coef(2)*S]
+%                type == 'trans_rep_prime' returns the normal derivative
+%                          corresponding to the transmission representation
+%                        [coef(1)*D' coef(2)*S']
+%                type == 'trans_rep_grad' returns the gradient corresponding
+%                         to the transmission representation
+%                        [coef(1)*d_x D coef(2)*d_x S;
+%                         coef(1)*d_y D coef(2)*d_y S]
+%
+%   varargin{1} - coef: length 2 array in the combined layer 
+%                 formula, 2x2 matrix for all kernels
+%                 otherwise does nothing
+%
+% Output:
+%   submat - the evaluation of the selected kernel for the
+%            provided sources and targets. the number of
+%            rows equals the number of targets and the
+%            number of columns equals the number of sources  
+%
+% see also CHNK.HELM1D.GREEN
+  
+src = srcinfo.r;
+targ = targinfo.r;
+
+[~,ns] = size(src);
+[~,nt] = size(targ);
+
+% double layer
+if strcmpi(type,'d')
+  srcnorm = srcinfo.n;
+  [~,grad] = chnk.helm1d.green(zkE,src,targ);
+  nx = repmat(srcnorm(1,:),nt,1);
+  ny = repmat(srcnorm(2,:),nt,1);
+  submat = -(grad(:,:,1).*nx + grad(:,:,2).*ny);
+end
+
+% normal derivative of single layer
+if strcmpi(type,'sprime')
+  targnorm = targinfo.n;
+  [~,grad] = chnk.helm1d.green(zkE,src,targ);
+  nx = repmat((targnorm(1,:)).',1,ns);
+  ny = repmat((targnorm(2,:)).',1,ns);
+
+  submat = (grad(:,:,1).*nx + grad(:,:,2).*ny);
+end
+
+
+% Tangential derivative of single layer
+if strcmpi(type,'stau')
+  targtan = targinfo.d;
+  [~,grad] = chnk.helm1d.green(zkE,src,targ);
+  dx = repmat((targtan(1,:)).',1,ns);
+  dy = repmat((targtan(2,:)).',1,ns);
+  ds = sqrt(dx.*dx+dy.*dy);
+  submat = (grad(:,:,1).*dx + grad(:,:,2).*dy)./ds;
+end
+
+% single layer
+if strcmpi(type,'s')
+  submat = chnk.helm1d.green(zkE,src,targ);
+end
+
+% normal derivative of double layer
+if strcmpi(type,'dprime')
+  targnorm = targinfo.n;
+  srcnorm = srcinfo.n;
+  [~,~,hess] = chnk.helm1d.green(zkE,src,targ);
+  nxsrc = repmat(srcnorm(1,:),nt,1);
+  nysrc = repmat(srcnorm(2,:),nt,1);
+  nxtarg = repmat((targnorm(1,:)).',1,ns);
+  nytarg = repmat((targnorm(2,:)).',1,ns);
+  submat = -(hess(:,:,1).*nxsrc.*nxtarg + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)...
+      + hess(:,:,3).*nysrc.*nytarg);
+end
+
+% Combined field 
+if strcmpi(type,'c')
+  srcnorm = srcinfo.n;
+  coef = ones(2,1);
+  if(nargin == 5); coef = varargin{1}; end
+  [submats,grad] = chnk.helm1d.green(zkE,src,targ);
+  nx = repmat(srcnorm(1,:),nt,1);
+  ny = repmat(srcnorm(2,:),nt,1);
+  submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny);
+  submat = coef(1)*submatd + coef(2)*submats;
+end
+
+% normal derivative of combined field
+if strcmpi(type,'cprime')
+  coef = ones(2,1);
+  if(nargin == 5); coef = varargin{1}; end
+  targnorm = targinfo.n;
+  srcnorm = srcinfo.n;
+  
+
+
+  % Get gradient and hessian info
+  [~,grad,hess] = chnk.helm1d.green(zkE,src,targ);
+
+  nxsrc = repmat(srcnorm(1,:),nt,1);
+  nysrc = repmat(srcnorm(2,:),nt,1);
+  nxtarg = repmat((targnorm(1,:)).',1,ns);
+  nytarg = repmat((targnorm(2,:)).',1,ns);
+
+  % D'
+  submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ...
+      + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)...
+      + hess(:,:,3).*nysrc.*nytarg);
+  % S'
+  submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg);
+
+  submat = coef(1)*submatdp + coef(2)*submatsp;
+end
+
+
+% Dirichlet and neumann data corresponding to combined field
+if strcmpi(type,'c2trans') 
+  coef = ones(2,1);
+  if(nargin == 5); coef = varargin{1}; end
+  targnorm = targinfo.n;
+  srcnorm = srcinfo.n;
+
+  % Get gradient and hessian info
+  [submats,grad,hess] = chnk.helm1d.green(zkE,src,targ);
+
+  nxsrc = repmat(srcnorm(1,:),nt,1);
+  nysrc = repmat(srcnorm(2,:),nt,1);
+  nxtarg = repmat((targnorm(1,:)).',1,ns);
+  nytarg = repmat((targnorm(2,:)).',1,ns);
+
+  % D'
+  submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ...
+      + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)...
+      + hess(:,:,3).*nysrc.*nytarg);
+  % S'
+  submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg);
+  % D
+  submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc);
+
+  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;
+
+end
+
+
+% all kernels, [c11 D, c12 S; c21 D', c22 S'] 
+if strcmpi(type,'all')
+ 
+  targnorm = targinfo.n;
+  srcnorm = srcinfo.n;
+  cc = varargin{1};
+  
+  submat = zeros(2*nt,2*ns);
+
+  % Get gradient and hessian info
+  [submats,grad,hess] = chnk.helm1d.green(zkE,src,targ);
+
+  nxsrc = repmat(srcnorm(1,:),nt,1);
+  nysrc = repmat(srcnorm(2,:),nt,1);
+  nxtarg = repmat((targnorm(1,:)).',1,ns);
+  nytarg = repmat((targnorm(2,:)).',1,ns);
+
+  submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ...
+      + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)...
+      + hess(:,:,3).*nysrc.*nytarg);
+  % S'
+  submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg);
+  % D
+  submatd  = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc);
+  
+  
+  submat(1:2:2*nt,1:2:2*ns) = submatd*cc(1,1);
+  submat(1:2:2*nt,2:2:2*ns) = submats*cc(1,2);
+  submat(2:2:2*nt,1:2:2*ns) = submatdp*cc(2,1);
+  submat(2:2:2*nt,2:2:2*ns) = submatsp*cc(2,2);
+end
+
+% Dirichlet data/potential correpsonding to transmission rep
+if strcmpi(type,'trans_rep') 
+
+  coef = ones(2,1);
+  if(nargin == 5); coef = varargin{1}; end;
+  srcnorm = srcinfo.n;
+  [submats,grad] = chnk.helm1d.green(zkE,src,targ);
+  nx = repmat(srcnorm(1,:),nt,1);
+  ny = repmat(srcnorm(2,:),nt,1);
+  submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny);
+
+  submat = zeros(1*nt,2*ns);
+  submat(1:1:1*nt,1:2:2*ns) = coef(1)*submatd;
+  submat(1:1:1*nt,2:2:2*ns) = coef(2)*submats;
+end
+
+% Neumann data corresponding to transmission rep
+if strcmpi(type,'trans_rep_prime')
+  coef = ones(2,1);
+  if(nargin == 5); coef = varargin{1}; end;
+  targnorm = targinfo.n;
+  srcnorm = srcinfo.n;
+  
+  submat = zeros(nt,ns);
+
+  % Get gradient and hessian info
+  [submats,grad,hess] = chnk.helm1d.green(zkE,src,targ);
+
+  nxsrc = repmat(srcnorm(1,:),nt,1);
+  nysrc = repmat(srcnorm(2,:),nt,1);
+  nxtarg = repmat((targnorm(1,:)).',1,ns);
+  nytarg = repmat((targnorm(2,:)).',1,ns);
+
+  % D'
+  submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ...
+      + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)...
+      + hess(:,:,3).*nysrc.*nytarg);
+  % S'
+  submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg);
+  submat = zeros(nt,2*ns)
+  submat(1:1:1*nt,1:2:2*ns) = coef(1)*submatdp;
+  submat(1:1:1*nt,2:2:2*ns) = coef(2)*submatsp;
+end
+
+
+% Gradient correpsonding to transmission rep
+if strcmpi(type,'trans_rep_grad')
+  coef = ones(2,1);
+  if(nargin == 5); coef = varargin{1}; end;
+  
+  srcnorm = srcinfo.n;
+  
+  submat = zeros(nt,ns,6);
+  % S
+  [submats,grad,hess] = chnk.helm1d.green(zkE,src,targ);
+
+  nxsrc = repmat(srcnorm(1,:),nt,1);
+  nysrc = repmat(srcnorm(2,:),nt,1);
+  % D
+  submatd  = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc);
+
+  submat = zeros(2*nt,2*ns);
+  
+  submat(1:2:2*nt,1:2:2*ns) = -coef(1)*(hess(:,:,1).*nxsrc + hess(:,:,2).*nysrc);
+  submat(1:2:2*nt,2:2:2*ns) = coef(2)*grad(:,:,1);
+    
+  submat(2:2:2*nt,1:2:2*ns) = -coef(1)*(hess(:,:,2).*nxsrc + hess(:,:,3).*nysrc);
+  submat(2:2:2*nt,2:2:2*ns) = coef(2)*grad(:,:,2);
+end
+
+
diff --git a/chunkie/+chnk/+helm1d/sweep.m b/chunkie/+chnk/+helm1d/sweep.m
new file mode 100644
index 0000000..d826ab9
--- /dev/null
+++ b/chunkie/+chnk/+helm1d/sweep.m
@@ -0,0 +1,28 @@
+function pot = sweep(uin,inds,ts,wts,zkE)
+%sweeping algorithm for preconditioner
+
+    vexpsp = exp(1i*diff(ts)*zkE);
+
+    u = zeros([numel(wts),1]);
+    u(inds) = uin;
+
+    nt = numel(u);
+    chrg = u.*(wts);
+    
+    voutp = zeros([nt,1]);
+    voutp(1) = chrg(1);
+    
+    for i=2:nt
+        voutp(i) = vexpsp(i-1)*voutp(i-1)+chrg(i);
+    end
+    
+    voutm = zeros([nt,1]);
+    chrg  = flipud(chrg);
+    vexpsp= flipud(vexpsp);
+    
+    for i=2:nt
+        voutm(i) = vexpsp(i-1)*(voutm(i-1)+chrg(i-1));
+    end
+    
+    pot = voutp + flipud(voutm);
+end
\ No newline at end of file
diff --git a/chunkie/demo/demo_KG.m b/chunkie/demo/demo_KG.m
new file mode 100644
index 0000000..64b62ee
--- /dev/null
+++ b/chunkie/demo/demo_KG.m
@@ -0,0 +1,296 @@
+% Defines a smooth, 1D interface and solves the corresponding 2D
+% time-harmonic Klein Gordon equation on either side of the interface.
+% We impose continuity and a jump in the normal derivative across the
+% interface. We also impose outgoing radiation conditions at infinity 
+% along the interface.
+
+
+clear all
+close('all')
+ifsvdplot = false;
+
+m=2;
+E=1;
+
+xmin=-60;
+xmax=60;
+
+cparams = [];
+cparams.ta = xmin;
+cparams.tb = xmax;
+cparams.ifclosed = 0;
+pref = [];
+pref.nchmax=100000;
+
+cparams.nchmin=24;
+
+a=1.0;
+b=1;
+c=0.0;
+d=1;
+[chnkr,ab] = chunkerfunc(@(t) nonflatinterface(t,a,b,c,d),cparams,pref);
+
+chnkr = sort(chnkr);
+
+n = chnkr.k*chnkr.nch;
+
+xs = chnkr.r(1,:);
+xs = xs(:);
+xx = xs(:);
+ys = chnkr.r(2,:);
+ys = ys(:);
+yy = ys(:);
+kh = 1j*sqrt(m^2-E^2);
+
+x0 = 0;
+y0 = 1.5;
+source = [x0 y0];
+rr = sqrt((xs-source(1)).^2+(ys-source(2)).^2);
+rr = rr(:);
+u_test = (1i/4)*besselh(0,1,kh*rr);
+
+decay_per_c = m*(xmax-xmin)/chnkr.nch;
+
+nchpad = 2*ceil(-log(10^(-16))/abs(decay_per_c));
+ich1 = nchpad;
+ich2 = chnkr.nch-nchpad;
+istart = (ich1-1)*chnkr.k +1;
+iend = (ich2-1)*chnkr.k;
+
+uu = -2*m*u_test(istart:iend);
+[sol_gmres,vdens] = fast_solve_wrap(uu,m,E,chnkr,istart,iend,ab);
+
+num_xts = 240;
+xts_min=-6;
+xts_max=6;
+num_yts = 120;
+yts_min=-2;
+yts_max=4;
+xts = linspace(xts_min,xts_max,num_xts);
+yts = linspace(yts_min,yts_max,num_yts);
+[XX,YY] = ndgrid(xts,yts);
+XX = XX(:);
+YY = YY(:);
+
+opts = [];
+opts.fac = 0.5;
+flags = flagnear(chnkr,[XX.';YY.'],opts);
+[rowloc,colloc] = find(flags);
+rowloc = reshape(rowloc,[],1);
+colloc = reshape(colloc,[],1);
+zk = 1j*sqrt(m^2-E^2);
+fkern = @(s,t) chnk.helm2d.kern(zk,s,t,'s',1);
+targs = [XX(rowloc).';YY(rowloc).'];
+if length(targs)>0
+    fints = chunkerkerneval(chnkr,fkern,vdens,targs,opts);
+end
+
+wts= weights(chnkr);
+wts= wts(:);
+
+rs = chnkr.r(1:2,:);
+srcinfo = [];
+srcinfo.sources = rs;
+srcinfo.charges = (vdens.*(wts(:))).';
+ns = numel(wts);
+
+rtargs = [XX.';YY.'];
+eps = 1e-12; %1e-12
+tic;
+pg = 0;
+pgt = 1;
+zk = 1j*sqrt(m^2-E^2);
+[v,~] = hfmm2d(eps,zk,srcinfo,pg,rtargs,pgt);
+vpot = v.pottarg;
+if length(targs)>0
+    vpot(rowloc) = 0;
+    vpot(rowloc) = fints;
+end
+
+XX   = reshape(XX,[numel(xts),numel(yts)]);
+YY   = reshape(YY,[numel(xts),numel(yts)]);
+vpot = reshape(vpot,[numel(xts),numel(yts)]);
+
+source = [x0 y0];
+rr = sqrt((XX-source(1)).^2+(YY-source(2)).^2);
+u_inc = (1i/4)*besselh(0,1,kh*rr);
+utot = -vpot + u_inc;
+
+maxsolution=max(abs(vpot(:)));
+figure
+h=pcolor(XX,YY,real(utot)); set(h,'EdgeColor','none'); daspect([1 1 1]);
+colormap(redblue)
+colorbar
+caxis([-maxsolution,maxsolution])
+
+
+
+
+
+function [sol_gmres,vdens] = fast_solve_wrap(uu,m,E,chnkr,istart,iend,ab)
+figs(1) = 0;
+
+zk = 1j*sqrt(m^2-E^2);
+
+%%%%%%%%% initialize fast apply parameters %%%%%%%%%%%%%%%
+
+tic;
+
+%%%% first: the single layer operator
+
+fkern = @(s,t) chnk.helm2d.kern(zk,s,t,'s',1);
+opdims(1) = 1; opdims(2) = 1;
+opts = [];
+opts.nonsmoothonly = true;
+tic; sys_nn = chunkermat(chnkr,fkern,opts); toc
+
+tic; [iinds,jinds,sinds] = find(sys_nn); toc
+
+tic;
+rs = chnkr.r(:,:);
+wts= weights(chnkr);
+wts= wts(:);
+
+rsrc = rs(:,jinds);
+rtar = rs(:,iinds);
+
+rdif = vecnorm(rtar-rsrc);
+h0 = besselh(0,1,zk*rdif);
+h1 = besselh(1,1,zk*rdif);
+vcors = 0.25*1i*(h0.').*(wts(jinds));
+vcors(iinds==jinds) = 0;
+
+smat = sparse(iinds,jinds,sinds-vcors);
+
+%%%% second: the preconditioner
+%obtain parametrization of the interface.
+t_coarse_1 = ab(1,:);
+t_coarse_2 = ab(2,:);
+chnklen=t_coarse_2 - t_coarse_1;
+[xs,~] = lege.exps(chnkr.k);
+tvals=zeros(chnkr.k,chnkr.nch);
+for i=1:chnkr.nch
+    tvals(:,i)=t_coarse_1(i)+chnklen(i)*(xs+1)/2;
+end
+ts=tvals(:);
+
+%now, discretize flat interface using the above nodes.
+cparamsflat = [];
+cparamsflat.ta = ab(1,1);
+cparamsflat.tb = ab(end,end);
+cparamsflat.tsplits = [ab(1,:),ab(end,end)];
+cparamsflat.ifclosed = 0;
+cparamsflat.nchmin = chnkr.nch;
+prefflat = [];
+prefflat.nchmax = chnkr.nch;
+chnkrflat = chunkerfunc(@(t) nonflatinterface(t,0,0,0,0),cparamsflat,prefflat);
+chnkrflat = sort(chnkrflat);
+
+zkE = E;
+fkern2 = @(s,t) chnk.helm1d.kern(zkE,s,t,'s');
+opdims(1) = 1; opdims(2) = 1;
+opts = [];
+opts.nonsmoothonly = true;
+opts.sing = 'hs'; 
+start = tic; sysmat_nn = m^2*chunkermat(chnkrflat,fkern2,opts)/zkE;
+t2 = toc(start)
+
+tic;
+
+wts= weights(chnkr);
+wts= wts(:);
+
+[iEnds,jEnds,sEnds] = find(sysmat_nn);
+tEis = ts(iEnds);
+tEjs = ts(jEnds);
+wEts = (wts(jEnds));
+vEnds = m^2*exp(1i*zkE*abs(tEis-tEjs)).*wEts/E;
+
+sEmat = sparse(iEnds,jEnds,sEnds-vEnds);
+
+%%%% wrap
+
+inds = istart:iend;
+
+fapplymat = @(u) apply_all(u,rs,wts,zk,smat,zkE,sEmat,...
+    ts,m,inds);
+
+toc
+
+%%%%%%%%% with a fast apply
+
+
+[sol_gmres,flag,relres,iter,resvec]  = gmres(fapplymat,uu,500,10^(-12),4);
+if (figs(1) == 1)
+    figure; plot(ts(istart:iend),log10(abs(sol_gmres)));
+end
+
+flag
+relres
+
+sol_gmres_large = zeros([numel(wts),1]);
+sol_gmres_large(inds) = sol_gmres;
+[vdens] = sol_gmres_large + ...
+    m^2/zkE*chnk.helm1d.sweep(sol_gmres,inds,ts,wts,zkE) + ...
+    sEmat*sol_gmres_large;
+
+end
+
+
+
+
+
+function [v] = apply_all(uin,rs,wts,zk,smat,zkE,sEmat,...
+    ts,dm,inds)
+
+vout = dm^2/zkE*chnk.helm1d.sweep(uin,inds,ts,wts,zkE);
+
+u = zeros([numel(wts),1]);
+u(inds) = uin;
+
+charges = u + vout + sEmat*u;
+
+%%%%%%%%%% apply single layer
+srcinfo = [];
+srcinfo.sources = rs;
+srcinfo.charges = (charges.*(wts(:))).';
+ns = numel(wts);
+
+targ = rs;
+eps = 1e-14;
+tic;
+pg = 1;
+[v,~] = hfmm2d(eps,zk,srcinfo,pg);
+
+v = (v.pot).'+smat*(charges);
+
+v = charges-2*dm*v;
+
+v = v(inds);
+
+end
+
+
+
+function [val] = get_gxi(xi,dm,de,dl,y0,x0,x,y)
+domeg = sqrt(dm*dm-de*de);
+zk = sqrt(xi^2+domeg^2);
+vmat = zeros(4,4);
+vmat(1,:)= [1,-exp(-zk*y0),-1,0];
+vmat(2,:)= [-zk,zk*exp(-zk*y0),-zk,0];
+vmat(3,:) = [0,1,exp(-zk*y0),-1];
+vmat(4,:) = [0,-zk,zk*exp(-zk*y0),-zk+2*dl];
+
+vrhs = zeros([4,1]);
+vrhs(1) = 0;
+vrhs(2) = 1/sqrt(2*pi);
+vrhs(3) = 0;
+vrhs(4) = 0;
+
+vcoef = vmat\vrhs;
+val=(double(y>y0).*vcoef(1).*exp(zk*(y0-y))...
+    +double(and(y>0,y<=y0)).*(vcoef(2).*exp(-zk*y)...
+    +vcoef(3).*exp(zk*(y-y0)))...
+    +double(y<=0).*vcoef(4).*exp(zk*y))...
+    .*exp(1i*xi*(x-x0))/sqrt(2*pi);
+end
\ No newline at end of file
diff --git a/chunkie/nonflatinterface.m b/chunkie/nonflatinterface.m
new file mode 100644
index 0000000..57b1daf
--- /dev/null
+++ b/chunkie/nonflatinterface.m
@@ -0,0 +1,32 @@
+function [r,d,d2] = nonflatinterface(t,a,b,c,d)
+%NONFLATINTERFACE curve parameterization of a perturbed interface given as 
+% the graph (x(t),y(t)) = (t,d e^(-at^2/2) sin(bt+c)) 
+% 
+% Syntax: 
+%    [r,d,d2] = nonflatinterface(t,a,b,c,d)
+%
+% Input: 
+%    t - array of t parameter values where curve coordinates 
+%    are desired
+%    a,b,c,d - curve parameters as in formula above 
+%
+% Output:
+%   r,d,d2 - 2 x numel(t) arrays of position, first and second 
+%   derivatives at desired t parameter values
+%
+
+% author: Solomon Quinn 
+
+   xs = t;
+   xp = ones(size(t));
+   xpp = zeros(size(t));
+   
+   
+    ys = d*exp(-a*t.^2/2).*sin(b*t+c);
+    yp = d*exp(-a*t.^2/2).*(b*cos(b*t+c) - a*t.*sin(b*t+c));
+    ypp = d*exp(-a*t.^2/2).*(-2*a*b*t.*cos(b*t+c) + (a*a*t.^2-a-b^2).*sin(b*t+c));
+    
+    r = [(xs(:)).'; (ys(:)).'];
+    d = [(xp(:)).'; (yp(:)).'];
+    d2 = [(xpp(:)).'; (ypp(:)).'];
+end
diff --git a/devtools/test/helm1d_greenTest.m b/devtools/test/helm1d_greenTest.m
new file mode 100644
index 0000000..64b11ee
--- /dev/null
+++ b/devtools/test/helm1d_greenTest.m
@@ -0,0 +1,308 @@
+clear all
+close('all')
+ifsvdplot = false;
+
+m=2;
+E=1;
+
+xmin=-60;
+xmax=60;
+
+cparams = [];
+cparams.ta = xmin;
+cparams.tb = xmax;
+cparams.ifclosed = 0;
+pref = [];
+pref.nchmax=100000;
+
+cparams.nchmin=24;
+
+a=0.0;
+b=0;
+c=0.0;
+d=0;
+[chnkr,ab] = chunkerfunc(@(t) nonflatinterface(t,a,b,c,d),cparams,pref);
+
+chnkr = sort(chnkr);
+
+n = chnkr.k*chnkr.nch;
+
+xs = chnkr.r(1,:);
+xs = xs(:);
+xx = xs(:);
+ys = chnkr.r(2,:);
+ys = ys(:);
+yy = ys(:);
+kh = 1j*sqrt(m^2-E^2);
+
+x0 = 0;
+y0 = 2.5;
+source = [x0 y0];
+rr = sqrt((xs-source(1)).^2+(ys-source(2)).^2);
+rr = rr(:);
+u_test = (1i/4)*besselh(0,1,kh*rr);
+
+decay_per_c = m*(xmax-xmin)/chnkr.nch;
+
+nchpad = 2*ceil(-log(10^(-16))/abs(decay_per_c));
+ich1 = nchpad;
+ich2 = chnkr.nch-nchpad;
+istart = (ich1-1)*chnkr.k +1;
+iend = (ich2-1)*chnkr.k;
+
+uu = -2*m*u_test(istart:iend);
+[sol_gmres,vdens] = fast_solve_wrap(uu,m,E,chnkr,istart,iend,ab);
+
+xts = [-2.0, 0.1];
+yts = [2.4, 2.8];
+[XX,YY] = ndgrid(xts,yts);
+XX = XX(:);
+YY = YY(:);
+
+opts = [];
+opts.fac = 0.5;
+flags = flagnear(chnkr,[XX.';YY.'],opts);
+[rowloc,colloc] = find(flags);
+rowloc = reshape(rowloc,[],1);
+colloc = reshape(colloc,[],1);
+zk = 1j*sqrt(m^2-E^2);
+fkern = @(s,t) chnk.helm2d.kern(zk,s,t,'s',1);
+targs = [XX(rowloc).';YY(rowloc).'];
+if length(targs)>0
+    fints = chunkerkerneval(chnkr,fkern,vdens,targs,opts);
+end
+
+wts= weights(chnkr);
+wts= wts(:);
+
+rs = chnkr.r(1:2,:);
+srcinfo = [];
+srcinfo.sources = rs;
+srcinfo.charges = (vdens.*(wts(:))).';
+ns = numel(wts);
+
+rtargs = [XX.';YY.'];
+eps = 1e-12; %1e-12
+tic;
+pg = 0;
+pgt = 1;
+zk = 1j*sqrt(m^2-E^2);
+[v,~] = hfmm2d(eps,zk,srcinfo,pg,rtargs,pgt);
+vpot = v.pottarg;
+if length(targs)>0
+    vpot(rowloc) = 0;
+    vpot(rowloc) = fints;
+end
+
+XX   = reshape(XX,[numel(xts),numel(yts)]);
+YY   = reshape(YY,[numel(xts),numel(yts)]);
+vpot = reshape(vpot,[numel(xts),numel(yts)]);
+
+source = [x0 y0];
+rr = sqrt((XX-source(1)).^2+(YY-source(2)).^2);
+u_inc = (1i/4)*besselh(0,1,kh*rr);
+utot = -vpot + u_inc;
+
+dm = m;
+de = E;
+
+domeg = sqrt(dm*dm-de*de);
+
+rr=0.25;
+dl = dm;
+
+xs=xts;
+ys=yts;
+[x,y] = ndgrid(xs,ys);
+        
+f1 = @(xi) get_gxi(xi,dm,de,dl,y0,x0,x,y);
+i1 = integral(f1,-200,-de-rr,'ArrayValued',true);
+
+f2 = @(xi) (rr*1i*exp(-1i*xi))*get_gxi(-de-rr*exp(-1i*xi),dm,de,dl,y0,x0,x,y);
+i2 = integral(f2,0,pi,'ArrayValued',true);
+
+i3 = integral(f1,-de+rr,de-rr,'ArrayValued',true);
+
+f2 = @(xi) (rr*1i*exp(-1i*xi))*get_gxi(de-rr*exp(-1i*xi),dm,de,dl,y0,x0,x,y);
+i4 = integral(f2,0,-pi,'ArrayValued',true);
+
+i5 = integral(f1,de+rr,200,'ArrayValued',true);
+
+vs=i1+i2+i3+i4+i5;
+
+relerr = abs(vs+utot)./abs(vs);
+
+assert(max(relerr(:)) < 1e-5)
+
+
+
+
+
+
+function [sol_gmres,vdens] = fast_solve_wrap(uu,m,E,chnkr,istart,iend,ab)
+figs(1) = 0;
+
+zk = 1j*sqrt(m^2-E^2);
+
+%%%%%%%%% initialize fast apply parameters %%%%%%%%%%%%%%%
+
+tic;
+
+%%%% first: the single layer operator
+
+fkern = @(s,t) chnk.helm2d.kern(zk,s,t,'s',1);
+opdims(1) = 1; opdims(2) = 1;
+opts = [];
+opts.nonsmoothonly = true;
+tic; sys_nn = chunkermat(chnkr,fkern,opts); toc
+
+tic; [iinds,jinds,sinds] = find(sys_nn); toc
+
+tic;
+rs = chnkr.r(:,:);
+wts= weights(chnkr);
+wts= wts(:);
+
+rsrc = rs(:,jinds);
+rtar = rs(:,iinds);
+
+rdif = vecnorm(rtar-rsrc);
+h0 = besselh(0,1,zk*rdif);
+h1 = besselh(1,1,zk*rdif);
+vcors = 0.25*1i*(h0.').*(wts(jinds));
+vcors(iinds==jinds) = 0;
+
+smat = sparse(iinds,jinds,sinds-vcors);
+
+%%%% second: the preconditioner
+%obtain parametrization of the interface.
+t_coarse_1 = ab(1,:);
+t_coarse_2 = ab(2,:);
+chnklen=t_coarse_2 - t_coarse_1;
+[xs,~] = lege.exps(chnkr.k);
+tvals=zeros(chnkr.k,chnkr.nch);
+for i=1:chnkr.nch
+    tvals(:,i)=t_coarse_1(i)+chnklen(i)*(xs+1)/2;
+end
+ts=tvals(:);
+
+%now, discretize flat interface using the above nodes.
+cparamsflat = [];
+cparamsflat.ta = ab(1,1);
+cparamsflat.tb = ab(end,end);
+cparamsflat.tsplits = [ab(1,:),ab(end,end)];
+cparamsflat.ifclosed = 0;
+cparamsflat.nchmin = chnkr.nch;
+prefflat = [];
+prefflat.nchmax = chnkr.nch;
+chnkrflat = chunkerfunc(@(t) nonflatinterface(t,0,0,0,0),cparamsflat,prefflat);
+chnkrflat = sort(chnkrflat);
+
+zkE = E;
+fkern2 = @(s,t) chnk.helm1d.kern(zkE,s,t,'s');
+opdims(1) = 1; opdims(2) = 1;
+opts = [];
+opts.nonsmoothonly = true;
+opts.sing = 'hs'; 
+start = tic; sysmat_nn = m^2*chunkermat(chnkrflat,fkern2,opts)/zkE;
+t2 = toc(start)
+
+tic;
+
+wts= weights(chnkr);
+wts= wts(:);
+
+[iEnds,jEnds,sEnds] = find(sysmat_nn);
+tEis = ts(iEnds);
+tEjs = ts(jEnds);
+wEts = (wts(jEnds));
+vEnds = m^2*exp(1i*zkE*abs(tEis-tEjs)).*wEts/E;
+
+sEmat = sparse(iEnds,jEnds,sEnds-vEnds);
+
+%%%% wrap
+
+inds = istart:iend;
+
+fapplymat = @(u) apply_all(u,rs,wts,zk,smat,zkE,sEmat,...
+    ts,m,inds);
+
+toc
+
+%%%%%%%%% with a fast apply
+
+
+[sol_gmres,flag,relres,iter,resvec]  = gmres(fapplymat,uu,500,10^(-12),4);
+if (figs(1) == 1)
+    figure; plot(ts(istart:iend),log10(abs(sol_gmres)));
+end
+
+flag
+relres
+
+sol_gmres_large = zeros([numel(wts),1]);
+sol_gmres_large(inds) = sol_gmres;
+[vdens] = sol_gmres_large + ...
+    m^2/zkE*chnk.helm1d.sweep(sol_gmres,inds,ts,wts,zkE) + ...
+    sEmat*sol_gmres_large;
+
+end
+
+
+
+
+
+function [v] = apply_all(uin,rs,wts,zk,smat,zkE,sEmat,...
+    ts,dm,inds)
+
+vout = dm^2/zkE*chnk.helm1d.sweep(uin,inds,ts,wts,zkE);
+
+u = zeros([numel(wts),1]);
+u(inds) = uin;
+
+charges = u + vout + sEmat*u;
+
+%%%%%%%%%% apply single layer
+srcinfo = [];
+srcinfo.sources = rs;
+srcinfo.charges = (charges.*(wts(:))).';
+ns = numel(wts);
+
+targ = rs;
+eps = 1e-14;
+tic;
+pg = 1;
+[v,~] = hfmm2d(eps,zk,srcinfo,pg);
+
+v = (v.pot).'+smat*(charges);
+
+v = charges-2*dm*v;
+
+v = v(inds);
+
+end
+
+
+
+function [val] = get_gxi(xi,dm,de,dl,y0,x0,x,y)
+domeg = sqrt(dm*dm-de*de);
+zk = sqrt(xi^2+domeg^2);
+vmat = zeros(4,4);
+vmat(1,:)= [1,-exp(-zk*y0),-1,0];
+vmat(2,:)= [-zk,zk*exp(-zk*y0),-zk,0];
+vmat(3,:) = [0,1,exp(-zk*y0),-1];
+vmat(4,:) = [0,-zk,zk*exp(-zk*y0),-zk+2*dl];
+
+vrhs = zeros([4,1]);
+vrhs(1) = 0;
+vrhs(2) = 1/sqrt(2*pi);
+vrhs(3) = 0;
+vrhs(4) = 0;
+
+vcoef = vmat\vrhs;
+val=(double(y>y0).*vcoef(1).*exp(zk*(y0-y))...
+    +double(and(y>0,y<=y0)).*(vcoef(2).*exp(-zk*y)...
+    +vcoef(3).*exp(zk*(y-y0)))...
+    +double(y<=0).*vcoef(4).*exp(zk*y))...
+    .*exp(1i*xi*(x-x0))/sqrt(2*pi);
+end
\ No newline at end of file