diff --git a/chunkie/+chnk/+flam/kernbyindexr.m b/chunkie/+chnk/+flam/kernbyindexr.m
index 8f6b1eb..b280657 100644
--- a/chunkie/+chnk/+flam/kernbyindexr.m
+++ b/chunkie/+chnk/+flam/kernbyindexr.m
@@ -1,4 +1,4 @@
-function mat = kernbyindexr(i,j,targs,chnkr,kern,opdims,spmat)
+function mat = kernbyindexr(i,j,targobj,chnkr,kern,opdims,spmat)
 %% evaluate system matrix by entry index utility function for 
 % general kernels, with replacement for specific entries and the
 % ability to add a low-rank modification, rectangular version
@@ -17,6 +17,10 @@
 %
 % i - array of row indices to compute
 % j - array of col indices to compute
+% targobj - object describing the target points, can be specified as
+%       * array of points
+%       * chunker object
+%       * chunkgraph object
 % chnkr - chunker object describing boundary
 % whts - smooth integration weights on chnkr
 % kern - kernel function of the form kern(s,t,stau,ttau) where s and t 
@@ -42,15 +46,32 @@
 [juni,~,ijuni] = unique(jpts);
 
 % matrix-valued entries of kernel for unique points
-
-ri = targs(:,iuni); rj = chnkr.r(:,juni);
+rj = chnkr.r(:,juni);
 dj = chnkr.d(:,juni); d2j = chnkr.d2(:,juni);
 nj = chnkr.n(:,juni);
-%di = bsxfun(@rdivide,di,sqrt(sum(di.^2,1)));
-%dj = bsxfun(@rdivide,dj,sqrt(sum(dj.^2,1)));
 srcinfo = []; srcinfo.r = rj; srcinfo.d = dj; srcinfo.d2 = d2j;
 srcinfo.n = nj;
-targinfo = []; targinfo.r = ri;
+% Assign appropriate object to targinfo
+targinfo = [];
+if isa(targobj, "chunker") || isa(targobj, "chunkgraph")
+    targinfo.r = targobj.r(:,iuni);
+    targinfo.d = targobj.d(:,iuni);
+    targinfo.d2 = targobj.d2(:,iuni);
+    targinfo.n = targobj.n(:,iuni);
+elseif isstruct(targobj)
+    if isfield(targobj,'d')
+        targinfo.d = targobj.d(:,iuni);
+    end
+    if isfield(targobj,'d2')
+        targinfo.d = targobj.d(:,iuni);
+    end
+    if isfield(targobj,'n')
+        targinfo.n = targobj.n(:,iuni);
+    end
+    targinfo.r = targobj.r(:,iuni);
+else
+    targinfo.r = targobj(:,iuni);
+end
 
 matuni = kern(srcinfo,targinfo);
 
diff --git a/chunkie/+chnk/+flam/proxy_square_pts.m b/chunkie/+chnk/+flam/proxy_square_pts.m
index 1d7bb7f..3f1d362 100644
--- a/chunkie/+chnk/+flam/proxy_square_pts.m
+++ b/chunkie/+chnk/+flam/proxy_square_pts.m
@@ -1,4 +1,4 @@
-function [pr,ptau,pw,pin] = proxy_square_pts(porder)
+function [pr,ptau,pw,pin] = proxy_square_pts(porder, opts)
 %PROXY_SQUARE_PTS return the proxy points, unit tangents, weights, and 
 % function handle for determining if points are within proxy surface.
 % This function is for the proxy surface around a unit box centered at 
@@ -16,12 +16,31 @@
     porder = 64;
 end
 
+if nargin < 2
+    opts = [];
+end
+
+iflege = 0;
+if isfield(opts, 'iflege')
+    iflege = opts.iflege;
+end
+    
+
 assert(mod(porder,4) == 0, ...
     'number of proxy points on square should be multiple of 4')
 
 po4 = porder/4;
 
-pts = -1.5+3*(0:(po4-1))*1.0/po4;
+if ~iflege
+    pts = -1.5+3*(0:(po4-1))*1.0/po4;
+    pw = ones(porder,1)*(4.0*3.0/porder);
+else
+    [xleg, wleg] = lege.exps(po4);
+    pts = xleg.'*1.5;
+    wleg = wleg*1.5;   
+    pw = [wleg; wleg; wleg; wleg];
+end
+
 one = ones(size(pts));
 
 pr = [pts, one*1.5, -pts, -1.5*one;
@@ -30,8 +49,6 @@
 ptau = [one, one*0, -one, one*0;
     one*0, one, one*0, -one];
 
-pw = ones(porder,1)*(4.0*3.0/porder);
-
 pin = @(x) max(abs(x),[],1) < 1.5;
 
 end
\ No newline at end of file
diff --git a/chunkie/+chnk/+helm2d/fmm.m b/chunkie/+chnk/+helm2d/fmm.m
index 4d50789..d7eceb0 100644
--- a/chunkie/+chnk/+helm2d/fmm.m
+++ b/chunkie/+chnk/+helm2d/fmm.m
@@ -38,6 +38,8 @@
 %                   note that targinfo must contain targ.n
 %                type = 'dprime' normal derivative of double layer,
 %                   note that targinfo must contain targ.n
+%                type = 'cprime' normal derivative of combined layer,
+%                   note that targinfo must contain targ.n
 %   sigma - density
 %   varargin{1} - coef in the combined layer formula, otherwise
 %                does nothing
@@ -61,7 +63,7 @@
     case {'d', 'dprime'}
         srcuse.dipstr = sigma(:).';
         srcuse.dipvec = srcinfo.n(1:2,:);
-    case 'c'
+    case {'c', 'cprime'}
         coefs = varargin{1};
         srcuse.charges = coefs(2)*sigma(:).';
         srcuse.dipstr  = coefs(1)*sigma(:).';
@@ -83,7 +85,7 @@
     switch lower(type)
         case {'s', 'd', 'c'}
             varargout{1} = U.pottarg.';
-        case {'sprime', 'dprime'}
+        case {'sprime', 'dprime', 'cprime'}
             if ( ~isfield(targinfo, 'n') )
                 error('CHUNKIE:helm2d:fmm:normals', ...
                     'Targets require normal info when evaluating Helmholtz kernel ''%s''.', type);
diff --git a/chunkie/+chnk/+helm2d/kern.m b/chunkie/+chnk/+helm2d/kern.m
index a90e636..d697efd 100644
--- a/chunkie/+chnk/+helm2d/kern.m
+++ b/chunkie/+chnk/+helm2d/kern.m
@@ -133,14 +133,14 @@
 % normal derivative of combined field
 if strcmpi(type,'cprime')
   coef = ones(2,1);
-  if(nargin == 5); coef = varargin{1}; end;
+  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.helm2d.green(zk,src,targ);
+  [~,grad,hess] = chnk.helm2d.green(zk,src,targ);
 
   nxsrc = repmat(srcnorm(1,:),nt,1);
   nysrc = repmat(srcnorm(2,:),nt,1);
@@ -161,7 +161,7 @@
 % Dirichlet and neumann data corresponding to combined field
 if strcmpi(type,'c2trans') 
   coef = ones(2,1);
-  if(nargin == 5); coef = varargin{1}; end;
+  if(nargin == 5); coef = varargin{1}; end
   targnorm = targinfo.n;
   srcnorm = srcinfo.n;
 
diff --git a/chunkie/+chnk/+quadnative/buildmat.m b/chunkie/+chnk/+quadnative/buildmat.m
index c75873a..7dbbd01 100644
--- a/chunkie/+chnk/+quadnative/buildmat.m
+++ b/chunkie/+chnk/+quadnative/buildmat.m
@@ -36,12 +36,7 @@
 targinfo = []; targinfo.r = rt; targinfo.d = dt; targinfo.d2 = d2t; targinfo.n = nt;
 hs = h(j); ht = h(i);
 
-dsnrms = sqrt(sum((ds).^2,1));
-%taus = bsxfun(@rdivide,ds,dsnrms);
-
-%dtnrms = sqrt(sum(abs(dt).^2,1));
-%taut = bsxfun(@rdivide,dt,dtnrms);
-
+dsnrms = sqrt(sum(ds.^2,1));
 ws = kron(hs(:),wts(:));
 
 dsdt = dsnrms(:).*ws;
diff --git a/chunkie/+chnk/+rcip/IPinit.m b/chunkie/+chnk/+rcip/IPinit.m
index c261171..c2828d4 100644
--- a/chunkie/+chnk/+rcip/IPinit.m
+++ b/chunkie/+chnk/+rcip/IPinit.m
@@ -1,10 +1,16 @@
   function [IP,IPW]=IPinit(T,W)
+  %CHNK.RCIP.IPinit 
+  %
   % construct the prolongation matrix IP that maps function values
   % on n_{gl} Gauss-Legendre nodes on [-1,1] to function values at 
   % 2n_{gl} Gauss-Legendre, with shifted and scaled n_{gl}
   % Gauss-Legendre nodes on each subinterval [-1,0], [0,1], respectively.
   %
   % IPW is the weighted prolongation matrix acted on the left side. 
+  %
+  
+  % author: Johan Helsing (part of RCIP tutorial)
+  
   ngl = length(T);
   A=ones(ngl);
   AA=ones(2*ngl,ngl);
diff --git a/chunkie/+chnk/+rcip/Rcompchunk.m b/chunkie/+chnk/+rcip/Rcompchunk.m
index d4b1614..e106bb1 100644
--- a/chunkie/+chnk/+rcip/Rcompchunk.m
+++ b/chunkie/+chnk/+rcip/Rcompchunk.m
@@ -1,13 +1,11 @@
-function [R]=Rcompchunk(chnkr,iedgechunks,fkern,ndim, ...
-    Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,...
-    glxs,sbcmat,lvmat,u,opts)
+function [R,rcipsav]=Rcompchunk(chnkr,iedgechunks,fkern,ndim, ...
+    Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,starL1,circL1,...
+    sbclmat,sbcrmat,lvmat,rvmat,u,opts)
 %CHNK.RCIP.Rcompchunk carry out the forward recursion for computing
 % the preconditioner R where geometry is described as a chunker
 %
 % This routine is not intended to be user-callable 
 %
-% Adapted from Shidong Jiang's RCIP implementation
-%
 % Function is passed as a handle, number of equations is given by
 % ndim
 %
@@ -15,32 +13,61 @@
 % 
 % Note that matrix must be scaled to have identity on the diagonal,
 % will not work with scaled version of identity
-  
+%
+
+% author: Shidong Jiang
+% modified: Jeremy Hoskins, Manas Rachh
+
+
 k = chnkr.k;  
 dim = chnkr.dim;
+nedge = size(iedgechunks,2);
 
+glxs = chnkr.tstor;
 glws = chnkr.wstor;
 
-if nargin < 14
+% return what's needed to interpolate from coarse
+
+rcipsav = [];
+rcipsav.k = k;
+rcipsav.ndim = ndim;
+rcipsav.nedge = nedge;
+rcipsav.Pbc = Pbc;
+rcipsav.PWbc = PWbc;
+rcipsav.starL = starL;
+rcipsav.starL1 = starL1;
+rcipsav.starS = starS;
+rcipsav.circL = circL;
+rcipsav.circL1 = circL1;
+rcipsav.circS = circS;
+rcipsav.ilist = ilist;
+rcipsav.nsub = nsub;
+
+if (nargin < 15 || isempty(sbclmat) || isempty(sbcrmat) || ...
+        isempty(lvmat) || isempty(rvmat) || isempty(u))
     [sbclmat,sbcrmat,lvmat,rvmat,u] = chnk.rcip.shiftedlegbasismats(k); 
 end
 
-if nargin < 17
+if nargin < 20
     opts = [];
 end
 
-nedge = size(iedgechunks,2);
+savedeep = false;
+if isfield(opts,'savedeep')
+    savedeep = opts.savedeep;
+end
 
-ileftright = zeros(nedge,1);
-nextchunk = zeros(nedge,1);
+rcipsav.savedeep = savedeep;
+
+if savedeep
+    rcipsav.R = cell(nsub+1,1);
+    rcipsav.MAT = cell(nsub,1);
+    rcipsav.chnkrlocals = cell(nsub,1);
+end
+
+
+% grab only those kernels relevant to this vertex
 
-km1 = k-1;
-rcs = zeros(km1,dim,nedge);
-dcs = zeros(k,dim,nedge);
-d2cs = zeros(k,dim,nedge);
-dscal = zeros(nedge,1);
-d2scal = zeros(nedge,1);
-ctr = zeros(dim,nedge);
 if(size(fkern)==1)
     fkernlocal = fkern;
 else
@@ -55,6 +82,20 @@
 
 end
 
+rcipsav.fkernlocal = fkernlocal;
+
+% get coefficients of recentered edge chunks and figure out orientation
+
+km1 = k-1;
+rcs = zeros(km1,dim,nedge);
+dcs = zeros(k,dim,nedge);
+d2cs = zeros(k,dim,nedge);
+dscal = zeros(nedge,1);
+d2scal = zeros(nedge,1);
+ctr = zeros(dim,nedge);
+
+ileftright = zeros(nedge,1);
+nextchunk = zeros(nedge,1);
 
 for i = 1:nedge
     ic = iedgechunks(1,i);
@@ -92,10 +133,19 @@
     end
 end
 
+rcipsav.ctr = ctr;
+rcipsav.rcs = rcs;
+rcipsav.dcs = dcs;
+rcipsav.d2cs = d2cs;
+rcipsav.dscal = dscal;
+rcipsav.d2scal = d2scal;
+rcipsav.ileftright = ileftright;
+rcipsav.glxs = glxs;
+rcipsav.glws = glws;
 
 pref = []; 
 pref.k = k;
-pref.nchmax = 5;
+pref.nchstor = 5;
 
 R = [];
 
@@ -108,6 +158,9 @@
 ts = cell(nedge,1);
 chnkrlocal(1,nedge) = chunker();
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% begin recursion proper
+
 h0=ones(nedge,1);
 for level=1:nsub
     h = h0/2^(nsub-level);
@@ -176,8 +229,16 @@
     if level==1    %  Dumb and lazy initializer for R, for now
   %R=eye(nR); 
         R = inv(MAT(starL,starL));
+        if savedeep
+            rcipsav.R{1} = R;
+        end
     end
     R=chnk.rcip.SchurBana(Pbc,PWbc,MAT,R,starL,circL,starS,circS);   
+    if savedeep
+        rcipsav.R{level+1} = R;
+        rcipsav.MAT{level} = MAT(starL,circL);
+        rcipsav.chnkrlocals{level} = merge(chnkrlocal);
+    end
 end
 
 end
diff --git a/chunkie/+chnk/+rcip/SchurBana.m b/chunkie/+chnk/+rcip/SchurBana.m
index 37643b8..8d02b71 100644
--- a/chunkie/+chnk/+rcip/SchurBana.m
+++ b/chunkie/+chnk/+rcip/SchurBana.m
@@ -1,10 +1,36 @@
   function A=SchurBana(P,PW,K,A,starL,circL,starS,circS)
-  % use matrix block inversion formula to recursively compute the
+  %CHNK.RCIP.SCHURBANA block inverse used in RCIP
+  % 
+  % uses a matrix block inversion formula to recursively compute the
   % preconditioner R.
+  %
+  % R_{i} =   [PW^T 0] [A^(-1)   U]^(-1) [P  0]
+  %           [ 0   I] [V        D]      [0  I]  
+  %
+  %           2N x 3N      3N x 3N      3N x 2N
+  %
+  % where A = R_{i-1}, U = K(bad,good), V = K(good,bad), D = K(good,good) 
+  % and bad refers to the two close panels on a type b mesh and good the 
+  % larger far panel. 
   % 
+  % This routine uses the Schur-Banachiewicz speedup as suggested by
+  % Helsing (see e.g. the RCIP Tutorial eq 31):
+  %
+  % R_{i} =   [M1 M2]
+  %           [M3 M4]
+  %
+  % M1 = PW^T A P + PW^T A U (D - VAU)^(-1) VAP
+  % M2 = -PW^TAU(D-VAU)^(-1) 
+  % M3 = -(D-VAU)^(-1)VAP
+  % M4 = (D-VAU)^(-1)
+  %
+  % this reduces a 3Nx3N matrix inverse to a NxN matrix inverse and
+  % several mat-mats
+  %
   % inputs: 
   % P, PW - nontrivial part (i.e., other than the identity matrices) 
-  %         prolongation and weighted prolongation matrices
+  %         prolongation and weighted prolongation matrices from type c to
+  %         type b meshes
   % K - the system matrix on a type-b mesh along each edge.
   %     the type-b mesh contains three chunks, with two small chunks
   %     close to the corner and one big chunk (twice of the small chunk
@@ -34,10 +60,13 @@
   % output:
   % A - R_{i}, the ith iteration of R.
   %
+  
+  % author: Johan Helsing (part of the RCIP tutorial)
+    
   VA=K(circL,starL)*A;
   PTA=PW'*A;
   PTAU=PTA*K(starL,circL);
-  DVAUI=chnk.rcip.myinv(K(circL,circL)-VA*K(starL,circL));
+  DVAUI=inv(K(circL,circL)-VA*K(starL,circL));
   DVAUIVAP=DVAUI*(VA*P);
   A(starS,starS)=PTA*P+PTAU*DVAUIVAP;
   A(circS,circS)=DVAUI;
diff --git a/chunkie/+chnk/+rcip/myinv.m b/chunkie/+chnk/+rcip/myinv.m
index f3fe439..7ebd1d2 100644
--- a/chunkie/+chnk/+rcip/myinv.m
+++ b/chunkie/+chnk/+rcip/myinv.m
@@ -1,5 +1,6 @@
   function Mi=myinv(M)
 % *** equation (9) in Henderson & Searle ***
+
   np=length(M)/2;
   A=M(1:np,1:np);
   U=M(1:np,np+1:2*np);
diff --git a/chunkie/+chnk/+rcip/rhohatInterp.m b/chunkie/+chnk/+rcip/rhohatInterp.m
new file mode 100644
index 0000000..a573015
--- /dev/null
+++ b/chunkie/+chnk/+rcip/rhohatInterp.m
@@ -0,0 +1,140 @@
+function [rhohatinterp,srcinfo,wts] = rhohatInterp(rhohat,rcipsav,ndepth)
+%CHNK.RCIP.RHOHATINTERP interpolate the (weight-corrected) coarse level
+% density rhohat to the requested depth using the backward recursion
+%
+% When using an RCIP preconditioner (in the style of eq 34 of the RCIP 
+% tutorial), the resulting coarse level density is 
+% accurate for interactions separated from the vertex 
+% *at the coarse scale*. 
+% By running the recursion for the RCIP preconditioner backwards, an 
+% equivalent density can be reconstructed on the fine mesh which is 
+% appropriate for closer interactions (at a distance well-separated from
+% the vertex *at the finest level in the mesh*).
+%
+% - Let Gamma_i be the portion of the boundary in the vicinity of
+% the vertex at level i (with level nsub the coarse level portion of the
+% boundary and level 1 the lowest level). 
+% - Let rhohat_i be the weight corrected density on a type c mesh on 
+% Gamma_i 
+% - Let the matrix I+K_i be a discretization of the integral operator on a
+% type b mesh on Gamma_i
+% - Let starL denote the "bad indices" of K_i and circL the "good indices";
+% the bad correspond to the pairs of small close panels on a type b mesh
+% - Let starS denote the "bad indices" of R_i and circS the "good indices";
+% the bad correspond to the close panels on a type c mesh
+% - Let P be the non-trivial part of an interpolator from type c to type b
+% meshes.
+% 
+% Then, to high precision 
+%
+% rhohat_{i-1} = R_{i-1}[ [P 0] R_i^(-1) rhohat_i -
+%                            (I+K_i)(starL,circL)rhohat_i(circS)]
+%
+% INPUT
+%
+% rhohat - entries of density rhohat corresponding to RCIP block in matrix.
+%          These come from the two nearest chunks to the corner and 
+%          should be ordered so that rhohat(starS) and rhohat(circS)
+%          correspond to the density at the bad and good points (near and
+%          far chunk) in the correct order (accounting for chunk
+%          orientation)
+% rcipsav - the quantities obtained from a previous call to Rcompchunk 
+% ndepth - (default: rcip.nsub) compute up to rhohat_{nsub-ndepth+1} 
+%
+% OUTPUT
+%
+% rhohatinterp - array of interpolated values, 
+% [rhohat_nsub(circS); rhohat_{nsub-1}(circS); rhohat_{nsub-2}(circS) ... 
+%                      rhohat_1(circS); rhohat_1(starS)]
+%                note that these are mostly the "good" indices except at
+%                the lowest level (naturally)
+% srcinfo - struct, containing points (srcinfo.r), normals (srcinfo.n), etc
+%                 in the same order as rhohatinterp.
+% wts - a set of weights for integrating functions sampled at these
+%           points
+
+
+rhohatinterp = [];
+srcinfo = [];
+wts = [];
+
+k = rcipsav.k;
+ndim = rcipsav.ndim;
+nedge = rcipsav.nedge;
+Pbc = rcipsav.Pbc;
+PWbc = rcipsav.PWbc;
+starL = rcipsav.starL;
+starL1 = rcipsav.starL1;
+starS = rcipsav.starS;
+circL = rcipsav.circL;
+circL1 = rcipsav.circL1;
+circS = rcipsav.circS;
+ilist = rcipsav.ilist;
+nsub = rcipsav.nsub;
+
+if nargin < 3
+    ndepth = nsub;
+end
+
+if ndepth > nsub
+    msg = "depth requested deeper than RCIP recursion performed\n " + ...
+        "going to depth %d instead";
+    warning(msg,nsub);
+    ndepth = nsub;
+end
+
+savedeep = rcipsav.savedeep;
+
+if savedeep
+    
+    % all relevant quantities are stored, just run backward recursion
+    
+    rhohat0 = rhohat;
+    rhohatinterp = [rhohatinterp; rhohat0(circS)];
+    r = [];
+    d = [];
+    d2 = [];
+    n = [];
+    h = [];
+    cl = rcipsav.chnkrlocals{nsub};
+    wt = weights(cl);
+    r = [r, cl.rstor(:,circL1)];
+    d = [d, cl.dstor(:,circL1)];
+    d2 = [d2, cl.d2stor(:,circL1)];
+    n = [n, cl.nstor(:,circL1)];
+    wts = [wts; wt(circL1(:))];
+    
+    R0 = rcipsav.R{nsub+1};
+    for i = 1:ndepth
+        R1 = rcipsav.R{nsub-i+1};
+        MAT = rcipsav.MAT{nsub-i+1};
+        rhotemp = R0\rhohat0;
+        rhohat0 = R1*(Pbc*rhotemp(starS) - MAT*rhohat0(circS));
+        if i == ndepth
+            rhohatinterp = [rhohatinterp; rhohat0];
+            r = [r, cl.rstor(:,starL1)];
+            d = [d, cl.dstor(:,starL1)];
+            d2 = [d2, cl.d2stor(:,starL1)];
+            n = [n, cl.nstor(:,starL1)];
+            wt = weights(cl);
+
+            wts = [wts; wt(starL1(:))];
+        else
+            cl = rcipsav.chnkrlocals{nsub-i};
+            rhohatinterp = [rhohatinterp; rhohat0(circS)];
+            r = [r, cl.rstor(:,circL1)];
+            d = [d, cl.dstor(:,circL1)];
+            d2 = [d2, cl.d2stor(:,circL1)];
+            n = [n, cl.nstor(:,circL1)];
+            wt = weights(cl);
+            wts = [wts; wt(circL1(:))];
+        end
+        R0 = R1;
+    end
+    
+    srcinfo.r = r;
+    srcinfo.d = d;
+    srcinfo.d2 = d2;
+    srcinfo.n = n;
+
+end
diff --git a/chunkie/+chnk/+rcip/setup.m b/chunkie/+chnk/+rcip/setup.m
index 6939737..35d3323 100644
--- a/chunkie/+chnk/+rcip/setup.m
+++ b/chunkie/+chnk/+rcip/setup.m
@@ -1,5 +1,7 @@
-  
-  function [Pbc,PWbc,starL,circL,starS,circS,ilist] = setup(ngl,ndim,nedge,isstart) 
+function [Pbc,PWbc,starL,circL,starS,circS,ilist,starL1,circL1] = ...
+      setup(ngl,ndim,nedge,isstart) 
+  %CHNK.RCIP.SETUP
+  %
   % setup for the RCIP forward recursion
   % inputs:
   % ngl - number of Gauss-Legendre nodes
@@ -15,6 +17,7 @@
   % Pbc - prolongation matrix
   % PWbc - weighted prolongation matrix
   % starL, circL - bad and good indices for the local system matrix
+  % starL1, circL1 - bad and good indices for arrays of nodes, normals, etc
   % starS, circS - bad and good indices for the preconditioner R
   %
   [T,W] = lege.exps(ngl);  
@@ -27,20 +30,30 @@
   % circL - good indices for the system matrix M
   starL = [];
   circL = [];
+  starL1 = [];
+  circL1 = [];
   indg1 = 2*ngl*ndim + (1:ngl*ndim);
   indb1 = 1:2*ngl*ndim;
+  indg11 = 2*ngl+ (1:ngl);
+  indb11 = 1:2*ngl;
   
   indg0 = 1:ngl*ndim;
   indb0 = ngl*ndim + (1:2*ngl*ndim);
+  indg01 = 1:ngl;
+  indb01 = ngl + (1:2*ngl);
   
   for i=1:nedge
     if isstart(i) 
       starL = [starL indb1+3*(i-1)*ngl*ndim];
       circL = [circL indg1+3*(i-1)*ngl*ndim];
+      starL1 = [starL1 indb11+3*(i-1)*ngl];
+      circL1 = [circL1 indg11+3*(i-1)*ngl];
       ilist(:,i) = [1, 2];
     else
       starL = [starL indb0+3*(i-1)*ngl*ndim];
       circL = [circL indg0+3*(i-1)*ngl*ndim];
+      starL1 = [starL1 indb01+3*(i-1)*ngl];
+      circL1 = [circL1 indg01+3*(i-1)*ngl];
       ilist(:,i) = [2, 3];
     end
   end
diff --git a/chunkie/@chunker/merge.m b/chunkie/@chunker/merge.m
index 0dcff92..ddc9587 100644
--- a/chunkie/@chunker/merge.m
+++ b/chunkie/@chunker/merge.m
@@ -1,21 +1,35 @@
-function chnkrout = merge(chnkrs)
+function chnkrout = merge(chnkrs,pref)
+%MERGE combine array of chunker objects into one chunker
+% 
+% input:
+%   chnkrs - array of chunker objects, must all have same order chunks
+%   pref - optional, chunkerpref object 
+% output:
+%   chnkrout - chunker object containing all nodes in chunker array. the
+%       ordering of nodes in chnkrout has the nodes from chnkrs(1) first,
+%       then chnkrs(2), etc. Adjacency information is updated as
+%       appropriate to the new indices. chunker data rows are copied as
+%       well. if chnkrs have different numbers of data rows, then those
+%       with fewer data rows are padded with zeros on merge. 
+%
 
 if isempty(chnkrs)
   chnkrout = chunker();
   return
 end
 assert(isa(chnkrs,'chunker'), 'input must be of chunker type');
+if nargin < 2
+    pref = [];
+end
 
-
-pref = [];
+% mandatory setting
 pref.k = chnkrs(1).k;
 chnkrout = chunker(pref);
 
-
 for i = 1:length(chnkrs)
   chnkrtemp = chnkrs(i);
-  assert(chnkrtemp.dim == chnkrout.dim,...
-      'chunkers to merge must be in same dimension');
+  assert(chnkrtemp.dim == chnkrout.dim && chnkrtemp.k == chnkrout.k,...
+      'chunkers to merge must be in same dimension and same order');
   nch = chnkrtemp.nch;
   nchold = chnkrout.nch;
   istart = nchold+1;
diff --git a/chunkie/@chunkgraph/chunkgraph.m b/chunkie/@chunkgraph/chunkgraph.m
index c9eaab0..e3eeb23 100644
--- a/chunkie/@chunkgraph/chunkgraph.m
+++ b/chunkie/@chunkgraph/chunkgraph.m
@@ -1,9 +1,16 @@
 classdef chunkgraph
 %CHUNKGRAPH chunk graph class for storing complex domains
 %
-% We describe a complex domain by its edges (smooth, regular
-% curves) and vertices (free ends of edges or points where the ends of 
-% one or more edges meets)
+% We describe a complex domain by a planar graph. Smooth boundary components
+% (discretized by chunkers) specify the edges of the graph and the
+% coordinates of free ends of the edges or points where the ends of multiple
+% edges meet specify the vertices. The interiors of distinct edges are not 
+% allowed to intersect (i.e. intersecting curves should be split at the
+% point where they intersect. 
+%
+% The chunkgraph is specified by an array of chunkers called echnks (the
+% edges), an array of points called verts (the vertices), and an array of
+% indices specifying the vertices at the left and right ends of any 
 %
 % 
     properties(SetAccess=public)
diff --git a/chunkie/@kernel/helm2d.m b/chunkie/@kernel/helm2d.m
index 64bed00..73cc8ce 100644
--- a/chunkie/@kernel/helm2d.m
+++ b/chunkie/@kernel/helm2d.m
@@ -14,6 +14,10 @@
 %   parameter ETA, i.e., COEFS(1)*KERNEL.HELM2D('d', ZK) + 
 %   COEFS(2)*KERNEL.HELM2D('s', ZK).
 %
+%   KERNEL.HELM2D('cp', ZK, COEFS) or KERNEL.HELM2D('combined', ZK, COEFS)
+%   constructs the derivative of the combined-layer Helmholtz kernel with 
+%   wavenumber ZK and  parameter ETA, i.e., COEFS(1)*KERNEL.HELM2D('dp', ZK) + 
+%   COEFS(2)*KERNEL.HELM2D('sp', ZK).
 % See also CHNK.HELM2D.KERN.
 
 if ( nargin < 1 )
@@ -66,6 +70,17 @@
         obj.fmm  = @(eps,s,t,sigma) chnk.helm2d.fmm(eps, zk, s, t, 'c', sigma, coefs);
         obj.sing = 'log';
 
+    case {'cp', 'cprime'}
+        if ( nargin < 3 )
+            warning('Missing combined layer coefficients. Defaulting to [1,1i].');
+            coefs = [1,1i];
+        end
+        obj.type = 'cprime';
+        obj.params.coefs = coefs;
+        obj.eval = @(s,t) chnk.helm2d.kern(zk, s, t, 'cprime', coefs);
+        obj.fmm  = @(eps,s,t,sigma) chnk.helm2d.fmm(eps, zk, s, t, 'cprime', sigma, coefs);
+        obj.sing = 'hs';
+
     otherwise
         error('Unknown Helmholtz kernel type ''%s''.', type);
 
diff --git a/chunkie/@kernel/kernel.m b/chunkie/@kernel/kernel.m
index f824811..30803a6 100644
--- a/chunkie/@kernel/kernel.m
+++ b/chunkie/@kernel/kernel.m
@@ -9,6 +9,7 @@
 %      ----                                         ----
 %      'laplace'    ('lap', 'l')                    's', 'd', 'sp', 'c'
 %      'helmholtz'  ('helm', 'h')                   's', 'd', 'sp', 'dp', 'c'
+%                                                   'cp'
 %      'helmholtz difference' ('helmdiff', 'hdiff') 's', 'd', 'sp', 'dp'
 %      'elasticity' ('elast', 'e')                  's', 'strac', 'd', 'dalt'
 %      'stokes'     ('stok', 's')                   'svel', 'spres', 'strac',
diff --git a/chunkie/FLAM b/chunkie/FLAM
index 2c9361e..a83698e 160000
--- a/chunkie/FLAM
+++ b/chunkie/FLAM
@@ -1 +1 @@
-Subproject commit 2c9361e7bc8fb8c3255f33275c0f2f7b7bbb7572
+Subproject commit a83698efdd58293c0e4c0c52f59498ffb9da9a99
diff --git a/chunkie/chunkerinterior.m b/chunkie/chunkerinterior.m
index 9563985..cde9605 100644
--- a/chunkie/chunkerinterior.m
+++ b/chunkie/chunkerinterior.m
@@ -1,16 +1,18 @@
-function [in] = chunkerinterior(chnkr,pts,opts)
+function [in] = chunkerinterior(chnkobj,ptsobj,opts)
 %CHUNKERINTERIOR returns an array indicating whether each point specified 
 % % by pts is inside the domain. Assumes the domain is closed.
 %
-% Syntax: in = chunkerinterior(chnkr,pts,opts)
-%         in = chunkerinterior(chnkr,{x,y},opts) % meshgrid version
+% Syntax: in = chunkerinterior(chnkobj,pts,opts)
+%         in = chunkerinterior(chnkobj,{x,y},opts) % meshgrid version
 %
 % Input:
-%   chnkr - chunker object describing geometry
-%   pts - (chnkr.dim,:) array of points to test
-% 
-%   {x,y} - length 2 cell array. the points checked then have the
-%       coordinates of a mesh grid [xx,yy] = meshgrid(x,y)
+%   chnkobj - chunker object or chunkgraph object describing geometry
+%   ptsobj - object describing the target points, can be specified as
+%       * (chnkr.dim,:) array of points to test
+%       * {x,y} - length 2 cell array. the points checked then have the
+%           coordinates of a mesh grid [xx,yy] = meshgrid(x,y)
+%       * chunker object, in which case it uses chunker.r(:,:) 
+%       * chunkgraph object, in which case it uses chunkgraph.r(:,:) 
 %
 % Optional input:
 %   opts - options structure with entries:
@@ -34,23 +36,39 @@
 
 grid = false;
 
-assert(chnkr.dim == 2,'interior only well-defined for 2D');
 
 if nargin < 3
     opts = [];
 end
 
-if isa(pts,"cell")
-    assert(length(pts)==2,'second input should be either 2xnpts array or length 2 cell array');
-    x = pts{1};
-    y = pts{2};
-    grid = true;
+% Assign appropriate object to chnkr
+if class(chnkobj) == "chunker"
+   chnkr = chnkobj;
+elseif class(chnkobj) == "chunkgraph"
+   chnkr = merge(chnkobj.echnks);
+else
+    msg = "Unsupported object in chunkerinterior";
+    error(msg)
 end
 
-if grid
+assert(chnkr.dim == 2,'interior only well-defined for 2D');
+
+% Assign appropriate object to pts
+if isa(ptsobj, "cell")
+    assert(length(ptsobj)==2,'second input should be either 2xnpts array or length 2 cell array');
+    x = ptsobj{1};
+    y = ptsobj{2};
+    grid = true;
     [xx,yy] = meshgrid(x,y);
     pts = [xx(:).'; yy(:).'];
-end    
+elseif isa(ptsobj, "chunker")
+    pts = ptsobj.r(:,:);
+elseif isa(ptsobj, "chunkgraph")
+    pts = ptsobj.r(:,:);
+else
+    pts = ptsobj;
+end
+
 
 usefmm = true;
 if isfield(opts,'fmm')
diff --git a/chunkie/chunkerkerneval.m b/chunkie/chunkerkerneval.m
index 24d068c..fb10f27 100644
--- a/chunkie/chunkerkerneval.m
+++ b/chunkie/chunkerkerneval.m
@@ -1,17 +1,20 @@
-function fints = chunkerkerneval(chnkr,kern,dens,targs,opts)
+function fints = chunkerkerneval(chnkobj,kern,dens,targobj,opts)
 %CHUNKERKERNEVAL compute the convolution of the integral kernel with
 % the density defined on the chunk geometry. 
 %
 % Syntax: fints = chunkerkerneval(chnkr,kern,dens,targs,opts)
 %
 % Input:
-%   chnkr - chunker object description of curve
+%   chnkobj - chunker object or chunkgraph object description of curve
 %   kern - kernel class object or kernel function taking inputs 
 %                      kern(srcinfo,targinfo) 
 %   dens - density on boundary, should have size opdims(2) x k x nch
 %          where k = chnkr.k, nch = chnkr.nch, where opdims is the 
 %           size of kern for a single src,targ pair
-%   targs - targ(1:2,i) gives the coords of the ith target
+%   targobj - object describing the target points, can be specified as
+%       * array of points
+%       * chunker object
+%       * chunkgraph object
 %
 % Optional input:
 %   opts - structure for setting various parameters
@@ -54,6 +57,16 @@
     kerneval = kern;
 end
 
+% Assign appropriate object to chnkr
+if class(chnkobj) == "chunker"
+   chnkr = chnkobj;
+elseif class(chnkobj) == "chunkgraph"
+   chnkr = merge(chnkobj.echnks);
+else
+    msg = "Unsupported object in chunkerkerneval";
+    error(msg)
+end
+
 srcinfo = []; targinfo = [];
 srcinfo.r = chnkr.r(:,1); srcinfo.d = chnkr.d(:,1); 
 srcinfo.n = chnkr.n(:,1); srcinfo.d2 = chnkr.d2(:,1);
@@ -86,61 +99,77 @@
 if isfield(opts,'fac'); opts_use.fac = opts.fac; end
 if isfield(opts,'eps'); opts_use.eps = opts.eps; end
 
-[dim,~] = size(targs);
+% Assign appropriate object to targinfo
+targinfo = [];
+if isa(targobj, "chunker")
+    targinfo.r = targobj.r(:,:);
+    targinfo.d = targobj.d(:,:);
+    targinfo.d2 = targobj.d2(:,:);
+    targinfo.n = targobj.n(:,:);
+elseif isa(targobj, "chunkgraph")
+    targinfo.r = targobj.r(:,:);
+    targinfo.d = targobj.d(:,:);
+    targinfo.d2 = targobj.d2(:,:);
+    targinfo.n = targobj.n(:,:);
+else
+    targinfo.r = targobj;
+end
+
+
+
+[dim,~] = size(targinfo.r);
 
 if (dim ~= 2); warning('only dimension two tested'); end
 
 if opts_use.forcesmooth
-    fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens,targs, ...
+    fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens,targinfo, ...
         [],opts_use);
     return
 end
 
 if opts_use.forceadap
     fints = chunkerkerneval_adap(chnkr,kern,opdims,dens, ...
-        targs,[],opts_use);
+        targinfo,[],opts_use);
     return
 end
 
 
 if opts_use.forcepquad
     optsflag = []; optsflag.fac = opts_use.fac;
-    flag = flagnear(chnkr,targs,optsflag);
-    fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens,targs, ...
+    flag = flagnear(chnkr,targinfo.r,optsflag);
+    fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens,targinfo, ...
         flag,opts_use);
 
     fints = fints + chunkerkerneval_ho(chnkr,kern,opdims,dens, ...
-        targs,flag,opts_use);
+        targinfo,flag,opts_use);
 
     return
 end
 
 % smooth for sufficiently far, adaptive otherwise
 
-%optsflag = []; optsflag.fac = opts_use.fac;
 rho = 1.8;
 optsflag = [];  optsflag.rho = rho;
-flag = flagnear_rectangle(chnkr,targs,optsflag);
+flag = flagnear_rectangle(chnkr,targinfo.r,optsflag);
 
 npoly = chnkr.k*2;
 nlegnew = chnk.ellipse_oversample(rho,npoly,opts_use.eps);
 nlegnew = max(nlegnew,chnkr.k);
 
 [chnkr2,dens2] = upsample(chnkr,nlegnew,dens);
-fints = chunkerkerneval_smooth(chnkr2,kern,opdims,dens2,targs, ...
+fints = chunkerkerneval_smooth(chnkr2,kern,opdims,dens2,targinfo, ...
     flag,opts_use);
 
 fints = fints + chunkerkerneval_adap(chnkr,kern,opdims,dens, ...
-        targs,flag,opts_use);
-
+        targinfo,flag,opts_use);
 
 end
 
 function fints = chunkerkerneval_ho(chnkr,kern,opdims,dens, ...
-    targs,flag,opts)
+    targinfo,flag,opts)
 
 % target
-[~,nt] = size(targs);
+[~,nt] = size(targinfo.r);
 fints = zeros(opdims(1)*nt,1);
 k = size(chnkr.r,2);
 [t,w] = lege.exps(2*k);
@@ -155,8 +184,22 @@
 % interpolation matrix
 intp = lege.matrin(k,t);          % interpolation from k to 2*k
 intp_ab = lege.matrin(k,[-1;1]);  % interpolation from k to end points
+
+targs = targinfo.r;
+
 targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt);
 targn = zeros(chnkr.dim,nt);
+if isfield(targinfo, 'd')
+    targd = targinfo.d;
+end
+
+if isfield(targinfo, 'd2')
+    targd2 = targinfo.d2;
+end
+
+if isfield(targinfo, 'n')
+    targn = targinfo.n;
+end
 for j=1:size(chnkr.r,3)
     [ji] = find(flag(:,j));
     if(~isempty(ji))
@@ -174,7 +217,7 @@
 
 
 function fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens, ...
-    targs,flag,opts)
+    targinfo,flag,opts)
 
 if isa(kern,'kernel')
     kerneval = kern.eval;
@@ -203,12 +246,10 @@
 dens = reshape(dens,opdims(2),k,nch);
 
 [~,w] = lege.exps(k);
-[~,nt] = size(targs);
+[~,nt] = size(targinfo.r);
 
 fints = zeros(opdims(1)*nt,1);
 
-targinfo = []; targinfo.r = targs;
-
 % assume smooth weights are good enough
 
 % Sequence of checks, first see ifflam is set as it supercedes
@@ -272,24 +313,6 @@
     end
 else
 
-%     % hack to ignore interactions: create sparse array with very small
-%     % number instead of zero in ignored entries. kernbyindexr overwrites
-%     %   with that small number
-%     [i,j] = find(flag);
-%     
-%     % targets correspond to multiple outputs (if matrix valued kernel)
-%     inew = (opdims(1)*(i(:)-1)).' + (1:opdims(1)).'; inew = inew(:);
-%     jnew = (j(:)).' + 0*(1:opdims(1)).'; jnew = jnew(:);
-%     
-%     % source chunks have multiple points and can be multi-dimensional
-%     jnew = (opdims(2)*chnkr.k*(jnew(:)-1)).' + (1:(chnkr.k*opdims(2))).';
-%     jnew = jnew(:);
-%     inew = (inew(:)).' + 0*(1:(chnkr.k*opdims(2))).';
-%     inew = inew(:);
-%     
-%     mm = nt*opdims(1); nn = chnkr.npt*opdims(2);
-%     v = 1e-300*ones(length(inew),1); sp = sparse(inew,jnew,v,mm,nn);
-
     wts = chnkr.wts;
     wts = wts(:);
     
@@ -297,14 +320,29 @@
         xflam1 = chnkr.r(:,:);
         xflam1 = repmat(xflam1,opdims(2),1);
         xflam1 = reshape(xflam1,chnkr.dim,numel(xflam1)/chnkr.dim);
-        targsflam = repmat(targs(:,:),opdims(1),1);
-        targsflam = reshape(targsflam,chnkr.dim,numel(targsflam)/chnkr.dim);
-        matfun = @(i,j) chnk.flam.kernbyindexr(i,j,targs,chnkr,kerneval, ...
-            opdims);
+
+        targinfo_flam = [];
+        targinfo_flam.r = repelem(targinfo.r(:,:),1,opdims(1));
+        if isfield(targinfo, 'd')
+            targinfo_flam.d = repelem(targinfo.d(:,:),1,opdims(1));
+        end
+        
+        if isfield(targinfo, 'd2')
+            targinfo_flam.d2 = repelem(targinfo.d2(:,:),1,opdims(1));
+        end
+        
+        if isfield(targinfo, 'n')
+            targinfo_flam.n = repelem(targinfo.n(:,:),1,opdims(1));
+        end
+
+% TODO: Pull through data?
+
+        matfun = @(i,j) chnk.flam.kernbyindexr(i, j, targinfo_flam, ...,
+                           chnkr, kerneval, opdims);
     
 
         width = max(abs(max(chnkr)-min(chnkr)))/3;
-        tmax = max(targs(:,:),[],2); tmin = min(targs(:,:),[],2);
+        tmax = max(targinfo.r(:,:),[],2); tmin = min(targinfo.r(:,:),[],2);
         wmax = max(abs(tmax-tmin));
         width = max(width,wmax/3);  
         npxy = chnk.flam.nproxy_square(kerneval,width);
@@ -314,16 +352,16 @@
             ctr,chnkr,kerneval,opdims,pr,ptau,pw,pin);
 
         optsifmm=[]; optsifmm.Tmax=Inf;
-        F = ifmm(matfun,targsflam,xflam1,200,1e-14,pxyfun,optsifmm);
+        F = ifmm(matfun,targinfo_flam.r,xflam1,200,1e-14,pxyfun,optsifmm);
         fints = ifmm_mv(F,dens(:),matfun);
     else
         wts2 = repmat(wts(:).', opdims(2), 1);
         sigma = wts2(:).*dens(:);
-        fints = kern.fmm(1e-14, chnkr, targs(:,:), sigma);
+        fints = kern.fmm(1e-14, chnkr, targinfo.r(:,:), sigma);
     end
     % delete interactions in flag array (possibly unstable approach)
     
-    targinfo = [];
+    
     if ~isempty(flag)
         for i = 1:nch
             densvals = dens(:,:,i); densvals = densvals(:);
@@ -338,9 +376,24 @@
             delsmooth = find(flag(:,i)); 
             delsmoothrow = (opdims(1)*(delsmooth(:)-1)).' + (1:opdims(1)).';
             delsmoothrow = delsmoothrow(:);
-            targinfo.r = targs(:,delsmooth);
 
-            kernmat = kerneval(srcinfo,targinfo);
+            targinfo_use = [];
+            targinfo_use.r = targinfo.r(:,delsmooth);
+
+            if isfield(targinfo, 'd')
+                targinfo_use.d = targinfo.d(:,delsmooth);
+            end
+        
+            if isfield(targinfo, 'd2')
+                targinfo_use.d2 = targinfo.d2(:,delsmooth);
+            end
+        
+            if isfield(targinfo, 'n')
+                targinfo_use.n = targinfo.n(:,delsmooth);
+            end
+
+
+            kernmat = kerneval(srcinfo,targinfo_use);
             fints(delsmoothrow) = fints(delsmoothrow) - kernmat*densvals;
         end
     end    
@@ -349,7 +402,7 @@
 end
 
 function fints = chunkerkerneval_adap(chnkr,kern,opdims,dens, ...
-    targs,flag,opts)
+    targinfo,flag,opts)
 
 if isa(kern,'kernel')
     kerneval = kern.eval;
@@ -370,7 +423,23 @@
 assert(numel(dens) == opdims(2)*k*nch,'dens not of appropriate size')
 dens = reshape(dens,opdims(2),k,nch);
 
+% Extract target info
+targs = targinfo.r;
 [~,nt] = size(targs);
+targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt);
+targn = zeros(chnkr.dim,nt);
+if isfield(targinfo, 'd')
+    targd = targinfo.d;
+end
+
+if isfield(targinfo, 'd2')
+    targd2 = targinfo.d2;
+end
+
+if isfield(targinfo, 'n')
+    targn = targinfo.n;
+end
+
 
 fints = zeros(opdims(1)*nt,1);
 
@@ -384,8 +453,8 @@
 n = chnkr.n;
 d2 = chnkr.d2;
 h = chnkr.h;
-targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt);
-targn = zeros(chnkr.dim,nt);
+
+
 
 if isempty(flag) % do all to all adaptive
     for i = 1:nch
diff --git a/chunkie/chunkerkernevalmat.m b/chunkie/chunkerkernevalmat.m
index c0a96c9..a697ac2 100644
--- a/chunkie/chunkerkernevalmat.m
+++ b/chunkie/chunkerkernevalmat.m
@@ -1,4 +1,4 @@
-function mat = chunkerkernevalmat(chnkr,kern,targs,opts)
+function mat = chunkerkernevalmat(chnkr,kern,targobj,opts)
 %CHUNKERKERNEVALMAT compute the matrix which maps density values on 
 % the chunk geometry to the value of the convolution of the given
 % integral kernel with the density at the specified target points
@@ -6,9 +6,21 @@
 % Syntax: mat = chunkerkernevalmat(chnkr,kern,targs,opts)
 %
 % Input:
-%   chnkr - chunker object description of curve
-%   kern - integral kernel taking inputs kern(srcinfo,targinfo) 
-%   targs - targ(1:2,i) gives the coords of the ith target
+%   chnkr - chunker object describing boundary, currently
+%              only supports chunkers, and not chunkgraphs
+%   kern  - kernel function. By default, this should be a function handle
+%           accepting input of the form kern(srcinfo,targinfo), where srcinfo
+%           and targinfo are in the ptinfo struct format, i.e.
+%                ptinfo.r - positions (2,:) array
+%                ptinfo.d - first derivative in underlying
+%                     parameterization (2,:)
+%                ptinfo.n - unit normals (2,:)
+%                ptinfo.d2 - second derivative in underlying
+%                     parameterization (2,:)
+%   targobj - object describing the target points, can be specified as
+%       * array of points
+%       * chunker object
+%       * chunkgraph object
 %
 % Optional input:
 %   opts - structure for setting various parameters
@@ -38,6 +50,37 @@
 
 % author: Travis Askham (askhamwhat@gmail.com)
 
+% convert kernel to kernel object, put in singularity info 
+% opts.sing provides a default value for singularities if not 
+% defined for kernels
+
+if isa(kern,'function_handle')
+    kern2 = kernel(kern);
+    kern = kern2;
+elseif isa(kern,'cell')
+    sz = size(kern);
+    kern2(sz(1),sz(2)) = kernel();
+    for j = 1:sz(2)
+        for i = 1:sz(1)
+            if isa(kern{i,j},'function_handle')
+                kern2(i,j) = kernel(kern{i,j});
+            elseif isa(kern{i,j},'kernel')
+                kern2(i,j) = kern{i,j};
+            else
+                msg = "Second input is not a kernel object, function handle, " ...
+                    + "or cell array";
+                error(msg);
+            end
+        end
+    end
+    kern = kern2;
+    
+elseif ~isa(kern,'kernel')
+    msg = "Second input is not a kernel object, function handle, " ...
+                + "or cell array";
+    error(msg);
+end
+    
 % determine operator dimensions using first two points
 
 
@@ -48,7 +91,9 @@
 targinfo.r = chnkr.r(:,2); targinfo.d = chnkr.d(:,2); 
 targinfo.d2 = chnkr.d2(:,2); targinfo.n = chnkr.n(:,2);
 
-ftemp = kern(srcinfo,targinfo);
+ftmp = kern.eval;
+
+ftemp = ftmp(srcinfo,targinfo);
 opdims = size(ftemp);
 
 if nargin < 4
@@ -57,18 +102,34 @@
 
 forcesmooth = false;
 forceadap = false;
-forcepqud = false;
+forcepquad = false;
 nonsmoothonly = false;
 fac = 1.0;
 eps = 1e-12;
 if isfield(opts,'forcesmooth'); forcesmooth = opts.forcesmooth; end
 if isfield(opts,'forceadap'); forceadap = opts.forceadap; end
-if isfield(opts,'forcepquad'); forcepqud = opts.forcepquad; end
+if isfield(opts,'forcepquad'); forcepquad = opts.forcepquad; end
 if isfield(opts,'nonsmoothonly'); nonsmoothonly = opts.nonsmoothonly; end
 if isfield(opts,'fac'); fac = opts.fac; end
 if isfield(opts,'eps'); eps = opts.eps; end
 
-[dim,~] = size(targs);
+% Assign appropriate object to targinfo
+targinfo = [];
+if isa(targobj, "chunker")
+    targinfo.r = targobj.r(:,:);
+    targinfo.d = targobj.d(:,:);
+    targinfo.d2 = targobj.d2(:,:);
+    targinfo.n = targobj.n(:,:);
+elseif isa(targobj, "chunkgraph")
+    targinfo.r = targobj.r(:,:);
+    targinfo.d = targobj.d(:,:);
+    targinfo.d2 = targobj.d2(:,:);
+    targinfo.n = targobj.n(:,:);
+else
+    targinfo.r = targobj;
+end
+
+[dim,~] = size(targinfo.r);
 
 
 if (dim ~= 2); warning('only dimension two tested'); end
@@ -77,24 +138,26 @@
 optsadap = []; 
 optsadap.eps = eps;
 
+
+
 if forcesmooth
-    mat = chunkerkernevalmat_smooth(chnkr,kern,opdims,targs, ...
+    mat = chunkerkernevalmat_smooth(chnkr,ftmp,opdims,targinfo, ...
         [],optssmooth);
     return
 end
 
 if forceadap
-    mat = chunkerkernevalmat_adap(chnkr,kern,opdims, ...
-        targs,[],optsadap);
+    mat = chunkerkernevalmat_adap(chnkr,ftmp,opdims, ...
+        targinfo,[],optsadap);
     return
 end
 
-if forcepqud
+if forcepquad
     optsflag = []; optsflag.fac = fac;
-    flag = flagnear(chnkr,targs,optsflag);
-    spmat = chunkerkernevalmat_ho(chnkr,kern,opdims, ...
-        targs,flag,optsadap);
-    mat = chunkerkernevalmat_smooth(chnkr,kern,opdims,targs, ...
+    flag = flagnear(chnkr,targinfo.r,optsflag);
+    spmat = chunkerkernevalmat_ho(chnkr,ftmp,opdims, ...
+        targinfo,flag,optsadap);
+    mat = chunkerkernevalmat_smooth(chnkr,ftmp,opdims,targinfo, ...
         flag,opts);
     mat = mat + spmat;
     return
@@ -102,17 +165,20 @@
 
 % smooth for sufficiently far, adaptive otherwise
 
+% TODO: change to chunkerkerneval system, need routine to generate
+% upsampling matrix.
+
 optsflag = []; optsflag.fac = fac;
-flag = flagnear(chnkr,targs,optsflag);
-spmat = chunkerkernevalmat_adap(chnkr,kern,opdims, ...
-        targs,flag,optsadap);
+flag = flagnear(chnkr,targinfo.r,optsflag);
+spmat = chunkerkernevalmat_adap(chnkr,ftmp,opdims, ...
+        targinfo,flag,optsadap);
 
 if nonsmoothonly
     mat = spmat;
     return;
 end
 
-mat = chunkerkernevalmat_smooth(chnkr,kern,opdims,targs, ...
+mat = chunkerkernevalmat_smooth(chnkr,ftmp,opdims,targinfo, ...
     flag,opts);
 
 mat = mat + spmat;
@@ -123,7 +189,7 @@
 
 
 function mat = chunkerkernevalmat_smooth(chnkr,kern,opdims, ...
-    targs,flag,opts)
+    targinfo,flag,opts)
 
 if nargin < 6
     flag = [];
@@ -135,7 +201,6 @@
 k = chnkr.k;
 nch = chnkr.nch;
 
-targinfo = []; targinfo.r = targs;
 srcinfo = []; srcinfo.r = chnkr.r(:,:); srcinfo.n = chnkr.n(:,:);
 srcinfo.d = chnkr.d(:,:); srcinfo.d2 = chnkr.d2(:,:);
 
@@ -163,7 +228,7 @@
 end
 
 function mat = chunkerkernevalmat_adap(chnkr,kern,opdims, ...
-    targs,flag,opts)
+    targinfo,flag,opts)
 
 k = chnkr.k;
 nch = chnkr.nch;
@@ -175,12 +240,30 @@
     opts = [];
 end
 
+% Extract target info
+targs = targinfo.r;
 [~,nt] = size(targs);
+targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt);
+targn = zeros(chnkr.dim,nt);
+if isfield(targinfo, 'd')
+    targd = targinfo.d;
+end
+
+if isfield(targinfo, 'd2')
+    targd2 = targinfo.d2;
+end
+
+if isfield(targinfo, 'n')
+    targn = targinfo.n;
+end
+
 
 % using adaptive quadrature
 
 
 if isempty(flag)
+
+
     mat = zeros(opdims(1)*nt,opdims(2)*chnkr.npt);
 
     [t,w] = lege.exps(2*k+1);
@@ -191,7 +274,7 @@
     n = chnkr.n;
     d2 = chnkr.d2;
     h = chnkr.h;
-    targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt);    
+    
     for i = 1:nch
         jmat = 1 + (i-1)*k*opdims(2);
         jmatend = i*k*opdims(2);
@@ -228,8 +311,6 @@
     n = chnkr.n;
     d2 = chnkr.d2;
     h = chnkr.h;
-    targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt);
-    targn = zeros(chnkr.dim,nt);
     for i = 1:nch
         jmat = 1 + (i-1)*k*opdims(2);
         jmatend = i*k*opdims(2);
@@ -260,7 +341,7 @@
 end
 
 function mat = chunkerkernevalmat_ho(chnkr,kern,opdims, ...
-    targs,flag,opts)
+    targinfo,flag,opts)
 
 k = chnkr.k;
 nch = chnkr.nch;
@@ -272,7 +353,24 @@
     opts = [];
 end
 
+% Extract target info
+targs = targinfo.r;
 [~,nt] = size(targs);
+targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt);
+targn = zeros(chnkr.dim,nt);
+if isfield(targinfo, 'd')
+    targd = targinfo.d;
+end
+
+if isfield(targinfo, 'd2')
+    targd2 = targinfo.d2;
+end
+
+if isfield(targinfo, 'n')
+    targn = targinfo.n;
+end
+
+
 
 % using Helsing-Ojala quadrature
 if isempty(flag) % figure out what is this flag for in adaptive routine
@@ -295,8 +393,7 @@
     % interpolation matrix 
     intp = lege.matrin(k,t);          % interpolation from k to 2*k
     intp_ab = lege.matrin(k,[-1;1]);  % interpolation from k to end points
-    targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt);
-    targn = zeros(chnkr.dim,nt);
+    
     for i = 1:nch
         jmat = 1 + (i-1)*k*opdims(2);
         jmatend = i*k*opdims(2);
diff --git a/chunkie/chunkermat.m b/chunkie/chunkermat.m
index 5e1897f..25f8bf2 100644
--- a/chunkie/chunkermat.m
+++ b/chunkie/chunkermat.m
@@ -7,7 +7,7 @@
 % Syntax: sysmat = chunkermat(chnkr,kern,opts)
 %
 % Input:
-%   chnkr - chunker object describing boundary
+%   chnkobj - chunker object describing boundary
 %   kern  - kernel function. By default, this should be a function handle
 %           accepting input of the form kern(srcinfo,targinfo), where srcinfo
 %           and targinfo are in the ptinfo struct format, i.e.
@@ -438,6 +438,7 @@
 
 
 if(icgrph && isrcip)
+    [sbclmat,sbcrmat,lvmat,rvmat,u] = chnk.rcip.shiftedlegbasismats(k);
     nch_all = horzcat(chnkobj.echnks.nch);
     [~,nv] = size(chnkobj.verts);
     ngl = chnkrs(1).k;
@@ -467,7 +468,7 @@
         
         if(norm(opdims_test(:)-ndim)>=0.5) 
             fprintf('in chunkermat: rcip: opdims did not match up for for vertex =%d\n',ivert)
-            fprtinf('returning without doing any rcip correction\m');
+            fprintf('returning without doing any rcip correction\n');
             break
         end
         starind = zeros(1,2*ngl*ndim*nedge);
@@ -482,13 +483,13 @@
         end
         
 
-        [Pbc,PWbc,starL,circL,starS,circS,ilist] = chnk.rcip.setup(ngl,ndim, ...
-          nedge,isstart);
+        [Pbc,PWbc,starL,circL,starS,circS,ilist,starL1,circL1] = ...
+            chnk.rcip.setup(ngl,ndim,nedge,isstart);
         
         % this might need to be fixed in triple junction case
         R = chnk.rcip.Rcompchunk(chnkrs,iedgechunks,kern,ndim, ...
-            Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,... 
-            glxs);
+            Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,starL1,circL1,... 
+            sbclmat,sbcrmat,lvmat,rvmat,u,opts);
        
         sysmat_tmp = inv(R) - eye(2*ngl*nedge*ndim);
         if (~nonsmoothonly)
diff --git a/chunkie/trapperkerneval.m b/chunkie/trapperkerneval.m
index d425b26..aaf7939 100644
--- a/chunkie/trapperkerneval.m
+++ b/chunkie/trapperkerneval.m
@@ -1,10 +1,21 @@
-function fints = trapperkerneval(trap,kern,dens,targs,opts)
+function fints = trapperkerneval(trap,kern,dens,targobj,opts)
 %TRAPPERINTKERN compute the convolution of the integral kernel with
 
 wts = trap.wts;
 srcinfo = []; srcinfo.r = trap.r; srcinfo.d = trap.d; srcinfo.d2 = trap.d2;
 srcinfo.n = trap.n;
-targinfo = []; targinfo.r = targs;
+
+% Assign appropriate object to targinfo
+targinfo = [];
+if isa(targobj, "trapper")
+    targinfo.r = targobj.r(:,:);
+    targinfo.d = targobj.d(:,:);
+    targinfo.d2 = targobj.d2(:,:);
+    targinfo.n = targobj.n(:,:);
+else
+    targinfo.r = targobj;
+end
+
 mat = kern(srcinfo,targinfo);
 fints = mat*diag(wts)*dens;
 
diff --git a/devtools/test/chunkerinteriorTest.m b/devtools/test/chunkerinteriorTest.m
index 65a6a91..2538be6 100644
--- a/devtools/test/chunkerinteriorTest.m
+++ b/devtools/test/chunkerinteriorTest.m
@@ -67,5 +67,20 @@
 opts = [];
 opts.fmm = true;
 opts.flam = false;
-start = tic; in3 = chunkerinterior(chnkr,{x1,x1},opts); t4 = toc(start);
-fprintf('%5.2e s : time for chunkerinterior (meshgrid, with fmm)\n',t4);
\ No newline at end of file
+start = tic; in4 = chunkerinterior(chnkr,{x1,x1},opts); t4 = toc(start);
+fprintf('%5.2e s : time for chunkerinterior (meshgrid, with fmm)\n',t4);
+
+% Test targets specified as chunker
+narms = 3;
+amp = 0.1;
+chnkr2 = chunkerfunc(@(t) 0.3*starfish(t,narms,amp),cparams,pref); 
+
+opts = [];
+opts.fmm = true;
+opts.flam = false;
+start = tic; in5 = chunkerinterior(chnkr,chnkr2,opts); t5 = toc(start);
+fprintf('%5.2e s : time for chunkerinterior (chunker, with fmm)\n',t5);
+assert(all(in5(:) == 1));
+
+
+
diff --git a/devtools/test/chunkerkerneval_greenlapTest.m b/devtools/test/chunkerkerneval_greenlapTest.m
index 9cccf63..ecf5c92 100644
--- a/devtools/test/chunkerkerneval_greenlapTest.m
+++ b/devtools/test/chunkerkerneval_greenlapTest.m
@@ -23,7 +23,7 @@
 
 % sources
 
-ns = 10;
+ns = 100;
 ts = 0.0+2*pi*rand(ns,1);
 sources = starfish(ts,narms,amp);
 sources = 3.0*sources;
@@ -31,7 +31,7 @@
 
 % targets
 
-nt = 100;
+nt = 300;
 ts = 0.0+2*pi*rand(nt,1);
 targets = starfish(ts,narms,amp);
 targets = targets.*repmat(rand(1,nt),2,1);
@@ -52,9 +52,9 @@
 
 % kernel defs
 
-kernd = @(s,t) chnk.lap2d.kern(s,t,'d');
-kerns = @(s,t) chnk.lap2d.kern(s,t,'s');
-kernsprime = @(s,t) chnk.lap2d.kern(s,t,'sprime');
+kernd = kernel('lap','d');
+kerns = kernel('lap','s');
+kernsprime = kernel('lap','sprime');
 
 opdims = [1 1];
 
@@ -63,21 +63,22 @@
 srcinfo = []; srcinfo.r = sources; 
 targinfo = []; targinfo.r = chnkr.r(:,:); 
 targinfo.d = chnkr.d(:,:);
-kernmats = kerns(srcinfo,targinfo);
-kernmatsprime = kernsprime(srcinfo,targinfo);
+kernmats = kerns.eval(srcinfo,targinfo);
+kernmatsprime = kernsprime.eval(srcinfo,targinfo);
 densu = kernmats*strengths;
 densun = kernmatsprime*strengths;
 
 % eval u at targets
 
 targinfo = []; targinfo.r = targets;
-kernmatstarg = kerns(srcinfo,targinfo);
+kernmatstarg = kerns.eval(srcinfo,targinfo);
 utarg = kernmatstarg*strengths;
 
 
-% test green's id
+% test green's id. we use FLAM, direct and FMM for smooth work
 
-opts.quadkgparams = {'RelTol',1.0e-13,'AbsTol',1.0e-13};
+opts = [];
+opts.flam = true;
 start=tic; Du = chunkerkerneval(chnkr,kernd,densu,targets,opts); 
 toc(start)
 start=tic; Sun = chunkerkerneval(chnkr,kerns,densun,targets,opts); 
@@ -89,7 +90,41 @@
 
 relerr = norm(utarg-utarg2,'fro')/norm(utarg,'fro');
 
-fprintf('relative frobenius error %5.2e\n',relerr);
+fprintf('relative frobenius error, forcing flam %5.2e\n',relerr);
+
+assert(relerr < 1e-11);
+
+opts = [];
+opts.accel = false;
+start=tic; Du = chunkerkerneval(chnkr,kernd,densu,targets,opts); 
+toc(start)
+start=tic; Sun = chunkerkerneval(chnkr,kerns,densun,targets,opts); 
+toc(start)
+
+utarg2 = Sun-Du;
+
+%
+
+relerr = norm(utarg-utarg2,'fro')/norm(utarg,'fro');
+
+fprintf('relative frobenius error, forcing direct %5.2e\n',relerr);
+
+assert(relerr < 1e-11);
+
+opts = [];
+opts.forcefmm = true;
+start=tic; Du = chunkerkerneval(chnkr,kernd,densu,targets,opts); 
+toc(start)
+start=tic; Sun = chunkerkerneval(chnkr,kerns,densun,targets,opts); 
+toc(start)
+
+utarg2 = Sun-Du;
+
+%
+
+relerr = norm(utarg-utarg2,'fro')/norm(utarg,'fro');
+
+fprintf('relative frobenius error, forcing fmm %5.2e\n',relerr);
 
 assert(relerr < 1e-11);
 
diff --git a/devtools/test/chunkgrphconstructTest.m b/devtools/test/chunkgrphconstructTest.m
index d829e77..6a6e50f 100644
--- a/devtools/test/chunkgrphconstructTest.m
+++ b/devtools/test/chunkgrphconstructTest.m
@@ -110,7 +110,7 @@
 
 prefs      = [];
 % prefs.chsmall = 1d-4;
-[cgrph] = chunkgraphinit(verts,edge2verts,fchnks,prefs);
+[cgrph] = chunkgraph(verts,edge2verts,fchnks,prefs);
 
 vstruc = procverts(cgrph);
 rgns = findregions(cgrph);
diff --git a/devtools/test/chunkgrphrcipTest.m b/devtools/test/chunkgrphrcipTest.m
index 154b123..391316a 100644
--- a/devtools/test/chunkgrphrcipTest.m
+++ b/devtools/test/chunkgrphrcipTest.m
@@ -62,7 +62,7 @@
 
 prefs      = [];
 % prefs.chsmall = 1d-4;
-[cgrph] = chunkgraphinit(verts,edge2verts,fchnks,prefs);
+[cgrph] = chunkgraph(verts,edge2verts,fchnks,prefs);
 
 vstruc = procverts(cgrph);
 rgns = findregions(cgrph);
@@ -110,7 +110,7 @@
 
 prefs      = [];
 % prefs.chsmall = 1d-4;
-[cgrph] = chunkgraphinit(verts,edge2verts,fchnks,prefs);
+[cgrph] = chunkgraph(verts,edge2verts,fchnks,prefs);
 
 vstruc = procverts(cgrph);
 rgns = findregions(cgrph);
diff --git a/devtools/test/chunkgrphrcipTransmissionTest.m b/devtools/test/chunkgrphrcipTransmissionTest.m
index 5547c62..7f1c4cf 100644
--- a/devtools/test/chunkgrphrcipTransmissionTest.m
+++ b/devtools/test/chunkgrphrcipTransmissionTest.m
@@ -55,7 +55,7 @@
 
 prefs      = [];
 % prefs.chsmall = 1d-4;
-[cgrph] = chunkgraphinit(verts,edge2verts,fchnks,prefs);
+[cgrph] = chunkgraph(verts,edge2verts,fchnks,prefs);
 
 
 vstruc = procverts(cgrph);
diff --git a/devtools/test/chunkrgrphOpdimTest.m b/devtools/test/chunkrgrphOpdimTest.m
index 4fb2652..38911c5 100644
--- a/devtools/test/chunkrgrphOpdimTest.m
+++ b/devtools/test/chunkrgrphOpdimTest.m
@@ -19,7 +19,7 @@
 fchnks = {}; % this by default gives me straight lines
 
 prefs = struct('maxchunklen',0.5);
-[cgrph] = chunkgraphinit(verts, edge2verts, fchnks, prefs);
+[cgrph] = chunkgraph(verts, edge2verts, fchnks, prefs);
 
 % figure(1); clf; hold on; axis equal; axis off;
 % plot(cgrph);
diff --git a/devtools/test/flamutilitiesTest.m b/devtools/test/flamutilitiesTest.m
index 0698b35..cd8d7bb 100644
--- a/devtools/test/flamutilitiesTest.m
+++ b/devtools/test/flamutilitiesTest.m
@@ -31,7 +31,7 @@
 end
 cparams = [];
 cparams.nover = 2;
-[cgrph] = chunkgraphinit(verts,edge2verts,fchnks,cparams);
+[cgrph] = chunkgraph(verts,edge2verts,fchnks,cparams);
 
 vstruc = procverts(cgrph);
 rgns = findregions(cgrph);
diff --git a/devtools/test/rcipTest.m b/devtools/test/rcipTest.m
index 30db8e5..3de9a28 100644
--- a/devtools/test/rcipTest.m
+++ b/devtools/test/rcipTest.m
@@ -13,10 +13,7 @@
 % set wave number
 zk = 1.1;
 
-
-
-
-nch = 5*ones(1,ncurve);
+nch = 8*ones(1,ncurve);
 
 a = -1.0;
 b = 1.0;
@@ -70,7 +67,7 @@
 
 fkern = @(s,t) -2*chnk.helm2d.kern(zk,s,t,'D');
 
-%%
+%
 
 % sources
 
@@ -86,6 +83,10 @@
 ts = 0.0+2*pi*rand(1,nt);
 targets = [cos(ts);sin(ts)];
 targets = 0.2*targets;
+targets(:,1) = [0.95;0];
+targets(:,2) = [0,0.36];
+targets(:,3) = [-0.95;0];
+targets(:,2) = [0,0.36];
 
 scatter(sources(1,:),sources(2,:),'o');
 scatter(targets(1,:),targets(2,:),'x');
@@ -133,6 +134,8 @@
 ncorner = 2;
 corners= cell(1,ncorner);
 R = cell(1,ncorner);
+rcipsav = cell(1,ncorner);
+starinds = cell(1,ncorner);
 
 
 corners{1}.clist = [1,2];
@@ -147,10 +150,11 @@
 
 ndim = 1;
 
-nsub = 100;
+nsub = 40;
 
 opts = [];
 
+
 for icorner=1:ncorner
     clist = corners{icorner}.clist;
     isstart = corners{icorner}.isstart;
@@ -170,30 +174,63 @@
         end
     end
     
+    starinds{icorner} = starind;
     
-    [Pbc,PWbc,starL,circL,starS,circS,ilist] = chnk.rcip.setup(ngl,ndim, ...
+    [Pbc,PWbc,starL,circL,starS,circS,ilist,starL1,circL1] = chnk.rcip.setup(ngl,ndim, ...
       nedge,isstart);
     opts_use = [];
     
     iedgechunks = corners{icorner}.iedgechunks;
-    tic; R{icorner} = chnk.rcip.Rcompchunk(chnkr,iedgechunks,fkern,ndim, ...
-        Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,... 
-        glxs);
+    optsrcip = [];
+    optsrcip.savedeep = true;
+    tic; [R{icorner},rcipsav{icorner}] = chnk.rcip.Rcompchunk(chnkr,iedgechunks,fkern,ndim, ...
+        Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,starL1,circL1,...,
+        [],[],[],[],[],optsrcip);
     toc
     
     sysmat(starind,starind) = inv(R{icorner}) - eye(2*ngl*nedge*ndim);
 end
 sysmat = sysmat + eye(np);
 
-sol = gmres(sysmat,ubdry,np,eps*20,np);
+[sol,flag,relres,iter] = gmres(sysmat,ubdry,np,eps*20,np);
+
+%
+
+% interpolate to fine grid
+
+ndepth = 10;
+cor = cell(1,ncorner);
+
+for icorner = 1:ncorner
+    solhat = sol(starinds{icorner});
+    [solhatinterp,srcinfo,wts] = chnk.rcip.rhohatInterp(solhat,rcipsav{icorner},ndepth);
 
-opts.usesmooth=true;
+    targtemp = targets(:,:) - rcipsav{icorner}.ctr(:,1);
+    targinfo = [];
+    targinfo.r = targtemp;
+    
+    cmat = fkern(srcinfo,targinfo);
+    mu = solhatinterp(:).*wts(:);
+    cor{icorner} = cmat*mu;
+end
+
+%
+
+soltemp = sol;
+for icorner = 1:ncorner
+    soltemp(starinds{icorner}) = 0;
+end
+    
+
+opts = [];
 opts.verb=false;
-opts.quadkgparams = {'RelTol',1e-16,'AbsTol',1.0e-16};
-start=tic; Dsol = chunkerkerneval(chnkrtotal,fkern,sol,targets,opts); 
+start=tic; Dsol = chunkerkerneval(chnkrtotal,fkern,soltemp,targets,opts); 
 t1 = toc(start);
 fprintf('%5.2e s : time to eval at targs (slow, adaptive routine)\n',t1)
 
+for icorner = 1:ncorner
+    Dsol = Dsol + cor{icorner};
+end
 %
 
 wchnkr = chnkrtotal.wts;