diff --git a/chunkie/+chnk/+axissymlap2d/gaus_agm.m b/chunkie/+chnk/+axissymlap2d/gaus_agm.m
new file mode 100644
index 0000000..dca93ef
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/gaus_agm.m
@@ -0,0 +1,38 @@
+function [rk0,re0] = gaus_agm(x)
+
+eps  = 1E-15;
+a    = sqrt(2./(x+1));
+delt = 1./sqrt(1-a.*a);
+aa0  = delt + sqrt(delt.*delt-1);
+bb0   = 1./(delt+sqrt(delt.*delt-1));
+a0  = ones(size(delt));
+b0  = 1./delt;
+
+
+fact = ((a0+b0)/2).^2;
+
+for i=1:1000
+    a1 = (a0+b0)/2;
+    b1 = sqrt(a0.*b0);
+
+    aa1 = (aa0+bb0)/2;
+    bb1 = sqrt(aa0.*bb0);
+    a0 = a1;
+    b0 = b1;
+    aa0 = aa1;
+    bb0 = bb1;
+
+    c0 = (a1-b1)/2;
+    fact = fact-(c0.*c0)*2^(i);
+    drel = abs(a0-b0)./abs(a0);
+    drel2= abs(aa0-bb0)./abs(aa0);
+    if (max(drel+drel2) <2*eps)
+        break
+    end
+
+end
+
+rk0 = pi./(2*aa0.*sqrt(1-a.*a));
+re0 = pi*fact./(2*a0);
+
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/gfunc.m b/chunkie/+chnk/+axissymlap2d/gfunc.m
new file mode 100644
index 0000000..951bfb8
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/gfunc.m
@@ -0,0 +1,18 @@
+function [gval,gdz,gdr,gdrp] = gfunc(r,rp,dr,z,zp,dz)
+    domega3 = 2*pi;
+    n       = 3;
+
+    t = (dz.^2+dr.^2)./(2.*r.*rp);
+    [q0,q1,q0d] = chnk.axissymlap2d.qleg_half(t);
+
+    gval = domega3*(rp./r).^((n-2)/2).*q0;
+    gdz  =-domega3*(rp./r).^((n-2)/2).*q0d ...
+           ./(rp.*r).*dz;
+
+    rfac = -r*(n-2)/2.*q0+(-(1+t).*r+rp).*q0d;
+    gdrp  = -domega3*(rp./r).^((n-2)/2)./(rp.*r) ...
+          .*rfac;
+    rfac = -rp*(n-2)/2.*q0+(-(1+t).*rp+r).*q0d;
+    gdr  = -domega3*(rp./r).^((n-2)/2)./(rp.*r) ...
+          .*rfac;
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/gfunc_on_axis.m b/chunkie/+chnk/+axissymlap2d/gfunc_on_axis.m
new file mode 100644
index 0000000..c7f1519
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/gfunc_on_axis.m
@@ -0,0 +1,18 @@
+function [gval,gdz,gdr,gdrp] = gfunc_on_axis(r,rp,dr,z,zp,dz)
+    domega3 = 2*pi;
+    n       = 3;
+
+    t = (dz.^2+dr.^2)./(2.*r.*rp);
+    [q0,q1,q0d] = chnk.axissymlap2d.qleg_half(t);
+
+    gval = domega3*(rp./r).^((n-2)/2).*q0;
+    gdz  =-domega3*(rp./r).^((n-2)/2).*q0d ...
+           ./(rp.*r).*dz;
+
+    rfac = -r*(n-2)/2.*q0+(-(1+t).*r+rp).*q0d;
+    gdrp  = -domega3*(rp./r).^((n-2)/2)./(rp.*r) ...
+          .*rfac;
+    rfac = -rp*(n-2)/2.*q0+(-(1+t).*rp+r).*q0d;
+    gdr  = -domega3*(rp./r).^((n-2)/2)./(rp.*r) ...
+          .*rfac;
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/green.m b/chunkie/+chnk/+axissymlap2d/green.m
new file mode 100644
index 0000000..564d6ab
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/green.m
@@ -0,0 +1,36 @@
+function [val, grad] = green(src, targ, origin)
+%CHNK.AXISSYMHELM2D.GREEN evaluate the Laplace green's function
+% for the given sources and targets. 
+%
+% Note: that the first coordinate is r, and the second z.
+% The code relies on precomputed tables and hence loops are required for 
+% computing various pairwise interactions.
+% Finally, the code is not efficient in the sense that val, grad, hess 
+% are always internally computed independent of nargout
+%
+% Since the kernels are not translationally invariant in r, the size
+% of the gradient is 3, for d/dr, d/dr', and d/dz
+%
+
+[~, ns] = size(src);
+[~, nt] = size(targ);
+
+gtmp = zeros(nt, ns, 3);
+
+rt = repmat(targ(1,:).',1,ns);
+rs = repmat(src(1,:),nt,1);
+dz = repmat(src(2,:),nt,1)-repmat(targ(2,:).',1,ns);
+r  = (rt + origin(1));
+rp = (rs + origin(1));
+dr = (rs-rt);
+z  = zeros(size(rt));
+zp = zeros(size(rt));
+[g,gdz,gdr,gdrp] = chnk.axissymlap2d.gfunc(r,rp,dr,z,zp,dz);
+
+val = g;
+gtmp(:,:,1) = gdr;
+gtmp(:,:,2) = gdrp;
+gtmp(:,:,3) = gdz;
+grad = gtmp;
+
+
diff --git a/chunkie/+chnk/+axissymlap2d/kern.m b/chunkie/+chnk/+axissymlap2d/kern.m
new file mode 100644
index 0000000..4e5f095
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/kern.m
@@ -0,0 +1,81 @@
+function submat = kern(srcinfo, targinfo, origin, type)
+%CHNK.AXISSYMLAP2D.KERN axissymmetric Laplace layer potential kernels in 2D
+% 
+% Syntax: submat = chnk.axissymlap2d.kern(srcinfo,targingo,type)
+%
+% 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  
+%
+% Here the first and second components correspond to the r and z
+% coordinates respectively. 
+%
+% Kernels based on G(x,y) = \int_{0}^{\pi} 1/(d(t)) \, dt \, 
+% where d(t) = \sqrt(r^2 + r'^2 - 2rr' \cos(t) + (z-z')^2) with
+% x = (r,z), and y = (r',z')
+%
+% 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:
+%   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'
+%
+%
+% 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  
+%
+%
+
+src = srcinfo.r; 
+targ = targinfo.r;
+
+[~, ns] = size(src);
+[~, nt] = size(targ);
+
+if strcmpi(type, 'd')
+    srcnorm = srcinfo.n;
+    [~, grad] = chnk.axissymlap2d.green(src, targ, origin);
+    nx = repmat(srcnorm(1,:), nt, 1);
+    ny = repmat(srcnorm(2,:), nt, 1);
+    % Due to lack of translation invariance in r, no sign flip needed, 
+    % as gradient is computed with repsect to r'
+    submat = (grad(:,:,2).*nx + grad(:,:,3).*ny);
+end
+
+if strcmpi(type, 'sprime')
+    targnorm = targinfo.n;
+    [~, grad] = chnk.axissymlap2d.green(src, targ, origin);
+
+    nx = repmat((targnorm(1,:)).',1,ns);
+    ny = repmat((targnorm(2,:)).',1,ns);
+    submat = (grad(:,:,1).*nx - grad(:,:,3).*ny);
+
+end
+
+if strcmpi(type, 's')
+    submat = chnk.axissymlap2d.green(src, targ, origin);
+    
+end
+
+
diff --git a/chunkie/+chnk/+axissymlap2d/qeval0_far.m b/chunkie/+chnk/+axissymlap2d/qeval0_far.m
new file mode 100644
index 0000000..3c3727c
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/qeval0_far.m
@@ -0,0 +1,24 @@
+function [q0] = qeval0_far(x)
+
+        coefs = ...
+       [2.2214414690791831235E0,
+        0,                    
+        0.41652027545234683566E0,                    
+        0,                    
+        0.22778452563800217575E0,                    
+        0,                    
+        0.15660186137612649583E0,                    
+        0,                    
+        0.11928657409509635424E0]; 
+
+
+        t = 1./x;
+        t0 = 1./sqrt(x);
+        q0 = zeros(size(t));
+
+        for i=1:9
+                q0 = q0 + t0*coefs(i);
+                t0 = t0.*t;        
+        end
+
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/qeval0_near.m b/chunkie/+chnk/+axissymlap2d/qeval0_near.m
new file mode 100644
index 0000000..185d47e
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/qeval0_near.m
@@ -0,0 +1,34 @@
+function [q0] = qeval0_near(t)
+
+        coefs = ...
+             [1.7328679513998632735E0,
+            -0.5000000000000000000E0,
+           -0.09160849392498290919E0,
+           0.062500000000000000000E0,
+           0.019905513916401443210E0,
+          -0.017578125000000000000E0,
+          -0.006097834693194945559E0,
+           0.006103515625000000000E0,
+          0.0021674343381175963469E0,
+         -0.0023365020751953125000E0,
+         -0.0008357538695841108955E0,
+          0.0009462833404541015625E0,
+          0.0003390851133264318725E0,
+        -0.00039757043123245239258E0,
+        -0.00014242013770157621830E0,
+         0.00017140153795480728149E0,
+         0.00006133159221782055041E0,
+        -0.00007532294148404616863E0];
+
+
+        b = log(t);
+        v = 1;
+        
+        q0 = zeros(size(t));
+
+        for i=1:9
+                q0 = q0 + v*coefs(2*i-1)+v.*b*coefs(2*i);
+                v  = v.*t;        
+        end
+
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/qeval0der_near.m b/chunkie/+chnk/+axissymlap2d/qeval0der_near.m
new file mode 100644
index 0000000..0ab6254
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/qeval0der_near.m
@@ -0,0 +1,36 @@
+function [q0d] = qeval0der_near(t)
+
+        coefs = ...
+             [1.7328679513998632735d0,
+            -0.5000000000000000000d0,
+           -0.09160849392498290919d0,
+           0.062500000000000000000d0,
+           0.019905513916401443210d0,
+          -0.017578125000000000000d0,
+          -0.006097834693194945559d0,
+           0.006103515625000000000d0,
+          0.0021674343381175963469d0,
+         -0.0023365020751953125000d0,
+         -0.0008357538695841108955d0,
+          0.0009462833404541015625d0,
+          0.0003390851133264318725d0,
+        -0.00039757043123245239258d0,
+        -0.00014242013770157621830d0,
+         0.00017140153795480728149d0,
+         0.00006133159221782055041d0,
+        -0.00007532294148404616863d0];
+
+        b = log(t);
+        v = 1;
+
+        q0d = coefs(2)./t;
+
+        for i=2:9
+                q0d = q0d + (i-1)* ...
+                   (v*coefs(2*i-1)+v.*b*coefs(2*i)) ...
+                       +v*coefs(2*i);
+                v  = v.*t;        
+        end
+     
+
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/qeval1_far.m b/chunkie/+chnk/+axissymlap2d/qeval1_far.m
new file mode 100644
index 0000000..27286bc
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/qeval1_far.m
@@ -0,0 +1,24 @@
+function [q1] = qeval1_far(x)
+
+        coefs = ...
+       [0.55536036726979578088E0,                    
+        0,                    
+        0.26032517215771677229E0,                    
+        0,                    
+        0.17083839422850163181E0,                    
+        0,                    
+        0.12723901236810277786E0,                    
+        0,                    
+        0.10139358798083190111E0];  
+
+
+        t = 1./x;
+        t0 = 1./sqrt(x).^3;
+        q1 = zeros(size(t));
+
+        for i=1:9
+                q1 = q1 + t0*coefs(i);
+                t0 = t0.*t;        
+        end
+
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/qeval1_near.m b/chunkie/+chnk/+axissymlap2d/qeval1_near.m
new file mode 100644
index 0000000..cd50647
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/qeval1_near.m
@@ -0,0 +1,34 @@
+function [q1] = qeval1_near(t)
+
+        coefs = ...
+            [-0.2671320486001367265d0,
+            -0.5000000000000000000d0,
+             0.5248254817749487276d0,
+           -0.18750000000000000000d0,
+           -0.04098835652733573868d0,
+           0.029296875000000000000d0, 
+           0.009513531070472923783d0,
+          -0.008544921875000000000d0,
+         -0.0029774361551467310174d0, 
+          0.0030040740966796875000d0, 
+          0.0010682069932178195667d0,
+         -0.0011565685272216796875d0,
+         -0.0004138797762860294822d0, 
+         0.00046985596418380737305d0, 
+          0.00016838776925222835322d0,
+        -0.00019777100533246994019d0,
+        -0.00007084821236213522235d0, 
+         0.00008536600034858565778d0];
+
+
+        b = log(t);
+        v = 1;
+        
+        q1 = zeros(size(t));
+
+        for i=1:9
+                q1 = q1 + v*coefs(2*i-1)+v.*b*coefs(2*i);
+                v  = v.*t;        
+        end
+
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/qget_zero.m b/chunkie/+chnk/+axissymlap2d/qget_zero.m
new file mode 100644
index 0000000..0e56a58
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/qget_zero.m
@@ -0,0 +1,16 @@
+function [qm,qmp] = qget_zero(x)
+
+        a = sqrt(2./(x+1));
+        [fF,fE] = chnk.axissymlap2d.gaus_agm(x);
+
+        q0 = a.*fF;
+        q1 = x.*q0-2*fE./(a);
+
+        qa = q0;
+        qb = q1;
+
+        qmm = 0;
+        qm  = q0;
+        qmp = q1;
+
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/qleg_half.m b/chunkie/+chnk/+axissymlap2d/qleg_half.m
new file mode 100644
index 0000000..6cb085e
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/qleg_half.m
@@ -0,0 +1,34 @@
+function [q0,q1,q0d] = qleg_half(t)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%   caveat utilitor: this function evaluates Q_{-1/2} and 
+%                    Q_{1/2} at t+1;
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    inear = find(t<0.01);
+    ifar  = find(t>100);
+    imid  = find((t>=0.01).*(t<=100));
+    [qmm,qmpm] = chnk.axissymlap2d.qget_zero(1+t(imid));
+    [qmn] = chnk.axissymlap2d.qeval0_near(t(inear));
+    [qmpn] = chnk.axissymlap2d.qeval1_near(t(inear));
+    [qmdn] = chnk.axissymlap2d.qeval0der_near(t(inear));
+    [qmf] = chnk.axissymlap2d.qeval0_far(1+t(ifar));
+    [qmpf] = chnk.axissymlap2d.qeval1_far(1+t(ifar));
+
+
+    
+    q0 = zeros(size(t));
+    q0(imid)  = qmm;
+    q0(inear) = qmn;
+    q0(ifar)  = qmf;
+
+    q1 = zeros(size(t));
+    q1(imid)  = qmpm;
+    q1(inear) = qmpn;
+    q1(ifar)  = qmpf;
+
+    q0d = -(1/2)*q1 + (1/2)*(1+t).*q0;
+    q0d = q0d./(-t.*(t+2));
+    q0d(inear) = qmdn;
+end
\ No newline at end of file
diff --git a/chunkie/@kernel/axissymlap2d.m b/chunkie/@kernel/axissymlap2d.m
new file mode 100644
index 0000000..c6187ce
--- /dev/null
+++ b/chunkie/@kernel/axissymlap2d.m
@@ -0,0 +1,63 @@
+function obj = axissymlap2d(type)
+%KERNEL.AXISSYMHELM2D   Construct the axissymmetric Helmholtz kernel.
+%   KERNEL.AXISSYMHELM2D('s', ZK) or KERNEL.AXISSYMHELM2D('single', ZK) 
+%   constructs the axissymmetric single-layer Helmholtz kernel with 
+%   wavenumber ZK.
+%
+%   KERNEL.AXISSYMHELM2D('d', ZK) or KERNEL.AXISSYMHELM2D('double', ZK) 
+%   constructs the axissymmetric double-layer Helmholtz kernel with 
+%   wavenumber ZK.
+%
+%   KERNEL.AXISSYMHELM2D('sp', ZK) or KERNEL.AXISSYMHELM2D('sprime', ZK) 
+%   constructs the normal derivative of the axissymmetric single-layer 
+%   Helmholtz kernel with wavenumber ZK.
+%
+%   KERNEL.AXISSYMHELM2D('c', ZK, COEFS) or 
+%   KERNEL.AXISSYMHELM2D('combined', ZK, COEFS)
+%   constructs the combined-layer axissymmetric Helmholtz kernel with 
+%   wavenumber ZK and parameter COEFS, 
+%   i.e., COEFS(1)*KERNEL.AXISSYMHELM2D('d', ZK) + 
+%         COEFS(2)*KERNEL.AXISSYMHELM2D('s', ZK).
+%
+% NOTES: The axissymetric kernels are currently supported only for purely
+% real or purely imaginary wave numbers
+%
+% See also CHNK.AXISSYMHELM2D.KERN.
+
+if ( nargin < 1 )
+    error('Missing Laplace kernel type.');
+end
+
+obj = kernel();
+obj.name = 'axissymlaplace';
+obj.opdims = [1 1];
+
+switch lower(type)
+
+    case {'s', 'single'}
+        obj.type = 's';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 's');
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 's');
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    case {'d', 'double'}
+        obj.type = 'd';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'd');
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'd');
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    case {'sp', 'sprime'}
+        obj.type = 'sp';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'sprime');
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'sprime');
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    otherwise
+        error('Unknown axissym Laplace kernel type ''%s''.', type);
+
+end
+
+end
diff --git a/chunkie/@kernel/kernel.m b/chunkie/@kernel/kernel.m
index d5ec312..5a06b7c 100644
--- a/chunkie/@kernel/kernel.m
+++ b/chunkie/@kernel/kernel.m
@@ -100,6 +100,8 @@
                       obj = kernel.nans(varargin{:});
                   case {'axis sym helmholtz', 'axissymh', 'axissymhelm'}
                       obj = kernel.axissymhelm2d(varargin{:});
+                  case {'axissymlaplace'}
+                      obj = kernel.axissymlap2d(varargin{:});
                   case {'axis sym helmholtz difference', 'axissymhdiff' ...
                            'axissymhelmdiff'}
                       obj = kernel.axissymhelm2ddiff(varargin{:});    
@@ -147,6 +149,7 @@
         obj = stok2d(varargin);
         obj = elast2d(varargin);
         obj = axissymhelm2d(varargin);
+        obj = axissymlap2d(varargin);
         obj = axissymhelm2ddiff(varargin);
         obj = zeros(varargin);
         obj = nans(varargin);
diff --git a/devtools/test/chunkermat_axissymlap_DirichletTest.m b/devtools/test/chunkermat_axissymlap_DirichletTest.m
new file mode 100644
index 0000000..5e8740b
--- /dev/null
+++ b/devtools/test/chunkermat_axissymlap_DirichletTest.m
@@ -0,0 +1,472 @@
+%clearvars; close all;
+iseed = 8675309;
+rng(iseed);
+
+%addpaths_loc();
+
+type = 'chnkr-star';
+type = 'chnkr-torus';
+ type = 'cgrph';
+%type = 'cgrph-sphere';
+
+% irep = 'rpcomb';
+irep = 'sk';
+
+pref = [];
+pref.k = 16;
+ns = 100;
+nt = 200;
+%ns = 1;
+%nt = 1;
+ppw = 10;   % points per wavelength;
+maxchunklen = 0.5;
+
+[chnkr, sources,targets] = get_geometry(type, pref, ns, nt, maxchunklen);
+[chnkr, targets, sources] = get_geometry(type, pref, nt, ns, maxchunklen);
+wts = chnkr.wts; wts = wts(:);
+
+
+%targets = [0.4,0.3;-0.2,0.2];
+%sources = [0.01;1.4];
+
+l2scale = false;
+fprintf('Done building geometry\n');
+
+% source strengths
+strengths = randn(ns, 1);
+
+% targets
+
+% Plot everything
+
+if (true)
+figure(1)
+clf
+hold off
+plot(chnkr, 'k.');
+hold on
+quiver(chnkr);
+scatter(sources(1,:), sources(2,:), 'o','green','filled')
+scatter(targets(1,:), targets(2,:), 'x','blue')
+axis equal
+
+end
+
+
+% For solving exterior the Neumann boundary value problem, we use the
+% following integral equation
+%
+% u = \beta(S_{k} + i \alpha D_{k} S_{ik}) \sigma
+%
+% with \beta = -1.0/(0.5 + 1i*0.25*alpha)
+% 
+% On imposing the boundary conditions, we get the following integral 
+% equation
+%
+%  du/dn = I + \beta S_{k}'[\sigma] + 
+%          i\beta \alpha(D'_{k} - D'_{ik})(S_{ik}[\sigma]) + 
+%          i\beta \alpha(S_{ik}')^2 [\sigma];
+% 
+% Setting -S_{ik}[\sigma] + \tau = 0, and - S_{ik}'[\sigma] + \mu=0, 
+% we have the following system of integral equations
+
+% Set up kernels
+D      = kernel('axissymlaplace','d');
+S      = kernel('axissymlaplace','s');
+
+K = 1/(2*pi^2)*D;
+Keval = K;
+
+opts = [];
+%opts.l2scale = l2scale;
+opts.rcip = false;
+opts.nsub_or_tol = 40;
+npts = chnkr.npt;
+nsys = K.opdims(1)*npts;
+tic, A = chunkermat(chnkr, K, opts) + eye(nsys); toc
+
+srcinfo  = []; srcinfo.r = sources; 
+targinfo = []; targinfo.r = chnkr.r(:,:); targinfo.n = chnkr.n(:,:);
+kernmats = S.eval(srcinfo, targinfo);
+ubdry = kernmats*strengths;
+rhs = ubdry;
+
+start = tic;
+sol = A\rhs;
+%sol = gmres(A, rhs, [], 1e-14, 200);
+t1 = toc(start);
+
+% Compute exact solution
+srcinfo  = []; srcinfo.r  = sources;
+targinfo = []; targinfo.r = targets;
+kernmatstarg = S.eval(srcinfo, targinfo);
+utarg = kernmatstarg*strengths;
+
+
+% Compute solution using chunkerkerneval
+% evaluate at targets and compare
+
+opts.forcesmooth = false;
+opts.verb = false;
+opts.quadkgparams = {'RelTol', 1e-15, 'AbsTol', 1.0e-15};
+
+if isa(chnkr, 'chunkgraph')
+    % Collapse cgrph into chnkrtotal
+    chnkrs = chnkr.echnks;
+    chnkrtotal = merge(chnkrs);
+else
+    chnkrtotal = chnkr;
+end
+
+
+
+start = tic;
+Dsol = chunkerkerneval(chnkrtotal, Keval, sol, targets, opts);
+t2 = toc(start);
+fprintf('%5.2e s : time to eval at targs (slow, adaptive routine)\n', t2)
+
+
+wchnkr = chnkrtotal.wts;
+wchnkr = repmat(wchnkr(:).', K.opdims(1), 1);
+relerr  = norm(utarg-Dsol) / (sqrt(chnkrtotal.nch)*norm(utarg));
+relerr2 = norm(utarg-Dsol, 'inf') / dot(abs(sol(:)), wchnkr(:));
+fprintf('relative frobenius error %5.2e\n', relerr);
+fprintf('relative l_inf/l_1 error %5.2e\n', relerr2);
+
+
+
+return
+
+
+
+
+
+if strcmpi(irep,'rpcomb')
+    Sik    = kernel('axissymhelm', 's', 1i*zk);
+    Sikp   = kernel('axissymhelm', 'sprime', 1i*zk);
+    Dkdiff = kernel('axissymhelmdiff', 'dprime', [zk 1i*zk]);
+    alpha = 1;
+    c1 = -1/(0.5 + 1i*alpha*0.25);
+    c2 = -1i*alpha/(0.5 + 1i*alpha*0.25);
+    c3 = -1;
+    K = [ c1*Skp  c2*Dkdiff c2*Sikp ;
+       c3*Sik  Z        Z        ;
+       c3*Sikp Z        Z        ];
+    K = kernel(K);
+    Keval = c1*kernel([Sk 1i*alpha*Dk Z]);
+else
+    K = -2*Skp;
+    Keval = -2*Sk;
+end
+
+% Set up boundary data
+
+srcinfo  = []; srcinfo.r = sources; 
+targinfo = []; targinfo.r = chnkr.r(:,:); targinfo.n = chnkr.n(:,:);
+kernmats = Skp.eval(srcinfo, targinfo);
+ubdry = kernmats*strengths;
+
+% rvals = chnkr.r(1,:);
+% ubdry = ubdry.*rvals(:);
+
+npts = chnkr.npt;
+nsys = K.opdims(1)*npts;
+rhs = zeros(nsys, 1);
+
+
+if(l2scale)
+    rhs(1:K.opdims(1):end) = ubdry.*sqrt(wts);
+else
+    rhs(1:K.opdims(1):end) = ubdry;
+end
+
+% Form matrix
+opts = [];
+opts.l2scale = l2scale;
+opts.rcip = true;
+opts.nsub_or_tol = 20;
+tic, A = chunkermat(chnkr, K, opts) + eye(nsys); toc
+start = tic;
+sol = gmres(A, rhs, [], 1e-14, 200);
+t1 = toc(start);
+
+% Compute exact solution
+srcinfo  = []; srcinfo.r  = sources;
+targinfo = []; targinfo.r = targets;
+kernmatstarg = Sk.eval(srcinfo, targinfo);
+utarg = kernmatstarg*strengths;
+
+% Compute solution using chunkerkerneval
+% evaluate at targets and compare
+
+opts.forcesmooth = false;
+opts.verb = false;
+opts.quadkgparams = {'RelTol', 1e-8, 'AbsTol', 1.0e-8};
+
+if(l2scale)
+    wts_rep = repmat(wts(:).', K.opdims(1),1);
+    wts_rep = wts_rep(:);
+    sol = sol./sqrt(wts_rep);
+end
+
+if isa(chnkr, 'chunkgraph')
+    % Collapse cgrph into chnkrtotal
+    chnkrs = chnkr.echnks;
+    chnkrtotal = merge(chnkrs);
+else
+    chnkrtotal = chnkr;
+end
+
+
+
+start = tic;
+Dsol = chunkerkerneval(chnkrtotal, Keval, sol, targets, opts);
+t2 = toc(start);
+fprintf('%5.2e s : time to eval at targs (slow, adaptive routine)\n', t2)
+
+
+wchnkr = chnkrtotal.wts;
+wchnkr = repmat(wchnkr(:).', K.opdims(1), 1);
+relerr  = norm(utarg-Dsol) / (sqrt(chnkrtotal.nch)*norm(utarg));
+relerr2 = norm(utarg-Dsol, 'inf') / dot(abs(sol(:)), wchnkr(:));
+fprintf('relative frobenius error %5.2e\n', relerr);
+fprintf('relative l_inf/l_1 error %5.2e\n', relerr2);
+
+return
+
+
+
+
+% Test fast direct solver interfaces
+
+% build sparse tridiag part
+
+opts.nonsmoothonly = true;
+opts.rcip = true;
+start = tic; spmat = chunkermat(chnkr, K, opts); t1 = toc(start);
+fprintf('%5.2e s : time to build tridiag\n',t1)
+
+
+spmat = spmat + speye(nsys);
+
+% test matrix entry evaluator
+start = tic; 
+opdims = K.opdims;
+sys2 = chnk.flam.kernbyindex(1:nsys, 1:nsys, chnkr, K, opdims, ...
+    spmat, l2scale);
+
+
+t1 = toc(start);
+
+fprintf('%5.2e s : time for mat entry eval on whole mat\n',t1)
+
+err2 = norm(sys2-A,'fro')/norm(A,'fro');
+fprintf('%5.2e   : fro error of build \n',err2);
+
+% test fast direct solver
+opts.ifproxy = false;
+F = chunkerflam(chnkr,K,1.0,opts);
+
+start = tic; sol2 = rskelf_sv(F,rhs); t1 = toc(start);
+
+if(l2scale)
+    wts_rep = repmat(wts(:).', K.opdims(1),1);
+    wts_rep = wts_rep(:);
+    sol2 = sol2./sqrt(wts_rep);
+end
+
+
+fprintf('%5.2e s : time for rskelf_sv \n',t1)
+
+err = norm(sol-sol2,'fro')/norm(sol,'fro');
+
+fprintf('difference between fast-direct and iterative %5.2e\n',err)
+
+
+function [chnkobj, sources, targets] = get_geometry(type, pref, ns, nt, maxchunklen)
+
+if nargin == 0
+    type = 'chnkr';
+end
+
+if nargin <= 1
+    pref = [];
+    pref.k = 16;
+end
+
+if nargin <= 2
+    ns = 10;
+end
+
+if nargin <= 3
+    nt = 10;
+end
+
+
+if nargin <= 4
+    maxchunklen = 1.0;
+end
+    
+
+if strcmpi(type, 'cgrph')
+    
+    
+    nverts = 3; 
+    verts = exp(-1i*pi/2 + 1i*pi*(0:(nverts-1))/(nverts-1));
+    verts = [real(verts);imag(verts)];
+
+
+    iind = 1:(nverts-1);
+    jind = 1:(nverts-1);
+
+    iind = [iind iind];
+    jind = [jind jind + 1];
+    jind(jind>nverts) = 1;
+    svals = [-ones(1,nverts-1) ones(1,nverts-1)];
+    edge2verts = sparse(iind, jind, svals, nverts-1, nverts);
+
+    amp = 0.1;
+    frq = 2;
+    fchnks    = cell(1,size(edge2verts,1));
+    for icurve = 1:size(edge2verts,1)
+        fchnks{icurve} = @(t) sinearc(t, amp, frq);
+    end
+    cparams = [];
+    cparams.nover = 2;
+    cparams.maxchunklen = maxchunklen;
+    cparams.ta = 0;
+    cparams.tb = 1;
+
+    chnkobj = chunkgraph(verts, edge2verts, fchnks, cparams, pref);
+    chnkobj = balance(chnkobj);
+    
+    ts = -pi/2 + pi*rand(ns,1);
+    sources = 0.2*[cos(ts)';sin(ts)'];
+    
+    ts = -pi/2 + pi*rand(nt,1);
+    targets = 3.0*[cos(ts)'; sin(ts)'];
+
+
+
+elseif strcmpi(type,'cgrph-sphere')
+    
+    
+    nverts = 2; 
+    verts = exp(-1i*pi/2 + 1i*pi*(0:(nverts-1))/(nverts-1));
+    verts = [real(verts);imag(verts)];
+
+
+    iind = 1:(nverts-1);
+    jind = 1:(nverts-1);
+
+    iind = [iind iind];
+    jind = [jind jind + 1];
+    jind(jind>nverts) = 1;
+    svals = [-ones(1,nverts-1) ones(1,nverts-1)];
+    edge2verts = sparse(iind, jind, svals, nverts-1, nverts);
+
+    cparams = [];
+    cparams.eps = 1.0e-10;
+    cparams.nover = 1;
+    cparams.ifclosed = false;
+    cparams.ta = 0;
+    cparams.tb = 1;
+    cparams.maxchunklen = maxchunklen;
+    narms = 0;
+    amp = 0.0;
+    
+    fchnks{1} = @(t) circle(t);
+
+    
+    chnkobj = chunkgraph(verts, edge2verts, fchnks, cparams, pref);
+    chnkobj = balance(chnkobj);
+    
+    ts = -pi/2 + pi*rand(ns, 1);
+    sources = starfish(ts, narms, amp);
+    sources = 1.2*sources;
+
+    ts = -pi/2 + pi*rand(nt, 1);
+    targets = starfish(ts, narms, amp);
+    targets = targets .* (0.8*repmat(rand(1, nt), 2, 1));   
+
+    
+    
+elseif strcmpi(type,'chnkr-star')
+    cparams = [];
+    cparams.eps = 1.0e-10;
+    cparams.nover = 1;
+    cparams.ifclosed = false;
+    cparams.ta = -pi/2;
+    cparams.tb = pi/2;
+    cparams.maxchunklen = maxchunklen;
+    narms = 0;
+    amp = 0.25;
+    chnkobj = chunkerfunc(@(t) starfish(t, narms, amp), cparams, pref); 
+    chnkobj = sort(chnkobj);
+    
+    ts = -pi/2 + pi*rand(ns, 1);
+    sources = starfish(ts, narms, amp);
+    sources = 0.5*sources;
+
+    ts = -pi/2 + pi*rand(nt, 1);
+    targets = starfish(ts, narms, amp);
+    targets = targets .* (1 + 0.5*repmat(rand(1, nt), 2, 1));   
+
+elseif strcmpi(type,'chnkr-torus')
+    cparams = [];
+    cparams.eps = 1.0e-10;
+    cparams.nover = 1;
+    cparams.ifclosed = true;
+    cparams.ta = 0;
+    cparams.tb = 2*pi;
+    cparams.maxchunklen = maxchunklen;
+    narms = 0;
+    amp = 0.25;
+    ctr = [3 0];
+    chnkobj = chunkerfunc(@(t) starfish(t, narms, amp, ctr), cparams, pref); 
+    chnkobj = sort(chnkobj);
+    
+    ts = -pi/2 + pi*rand(ns, 1);
+    sources = starfish(ts, narms, amp);
+    sources = 0.5*sources;
+    sources(1,:) = sources(1,:) + ctr(1);
+
+    ts = -pi/2 + pi*rand(nt, 1);
+    targets = starfish(ts, narms, amp);
+    targets = targets .* (1 + 0.5*repmat(rand(1, nt), 2, 1));    
+    targets(1,:) = targets(1,:) + ctr(1);
+end
+
+
+
+end
+
+
+function [r,d,d2] = sinearc(t,amp,frq)
+xs = t;
+ys = amp*sin(frq*t);
+xp = ones(size(t));
+yp = amp*frq*cos(frq*t);
+xpp = zeros(size(t));
+ypp = -frq*frq*amp*sin(t);
+
+r = [(xs(:)).'; (ys(:)).'];
+d = [(xp(:)).'; (yp(:)).'];
+d2 = [(xpp(:)).'; (ypp(:)).'];
+end
+
+
+function [r,d,d2] = circle(t)
+xs = cos(-pi/2 + pi*t);
+ys = sin(-pi/2 + pi*t);
+xp = -pi*ys;
+yp = pi*xs;
+xpp = -pi*pi*xs;
+ypp = -pi*pi*ys;
+r = [(xs(:)).'; (ys(:)).'];
+d = [(xp(:)).'; (yp(:)).'];
+d2 = [(xpp(:)).'; (ypp(:)).'];
+end
+
+
+