diff --git a/chunkie/+chnk/+axissymhelm2d/asym_helm_data2.mat b/chunkie/+chnk/+axissymhelm2d/asym_helm_data2.mat
new file mode 100644
index 00000000..e7398172
Binary files /dev/null and b/chunkie/+chnk/+axissymhelm2d/asym_helm_data2.mat differ
diff --git a/chunkie/+chnk/+axissymhelm2d/green.m b/chunkie/+chnk/+axissymhelm2d/green.m
index fe576152..9100f616 100644
--- a/chunkie/+chnk/+axissymhelm2d/green.m
+++ b/chunkie/+chnk/+axissymhelm2d/green.m
@@ -32,6 +32,13 @@
     asym_tables = chnk.axissymhelm2d.load_asym_tables();
 end
 
+persistent asym_tables_sub   
+if isempty(asym_tables_sub)
+    asym_tables_sub = chnk.axissymhelm2d.load_asym_tables_sub();
+end
+
+htables = asym_tables;
+
 [~, ns] = size(src);
 [~, nt] = size(targ);
 
@@ -45,12 +52,19 @@
     ifun = 2;
 end
 
-if ifdiff
+if (ifdiff == 1)
     if abs(zi) > eps
         error('AXISSYMHELM2D.green: Difference of greens function only supported for real wave numbers\n');
     end
     ifun = 3;
 end
+if (ifdiff == 2)
+    if abs(zi) > eps
+        error('AXISSYMHELM2D.green: Difference of greens function only supported for real wave numbers\n');
+    end
+    htables = asym_tables_sub;
+    ifun = 5;
+end
 
 over2pi = 1/2/pi;
 kabs = abs(k);
@@ -61,7 +75,7 @@
 r  = (rt + origin(1))*kabs;
 dz = dz*kabs;
 dr = (rs-rt)*kabs;
-dout = chnk.axissymhelm2d.helm_axi_all(r, dr, dz, asym_tables,ifun);
+dout = chnk.axissymhelm2d.helm_axi_all(r, dr, dz, htables,ifun);
 int = dout.int;
 intdz = dout.intdz;
 intdq = dout.intdq;
diff --git a/chunkie/+chnk/+axissymhelm2d/helm_axi_all.m b/chunkie/+chnk/+axissymhelm2d/helm_axi_all.m
index a5dd026f..cd58377c 100644
--- a/chunkie/+chnk/+axissymhelm2d/helm_axi_all.m
+++ b/chunkie/+chnk/+axissymhelm2d/helm_axi_all.m
@@ -5,10 +5,10 @@
 
     int   = zeros(size(alphs));
 
-    iflag_rk = (ifun == 1 || ifun == 4);
+    iflag_rk = (ifun == 1 || ifun == 4 || ifun == 5);
     iflag_ik = (ifun == 2 || ifun == 4);
-    iflag_dk = (ifun == 3 || ifun == 4);
-
+    iflag_dk = (ifun == 3 || ifun == 4 || ifun == 5);
+    
     if (iflag_rk)
     rk = cell(6);
     for ii=1:6
@@ -186,6 +186,10 @@
         [dsdiff] = proc_kerns(rs,drs,dzs,dk);
         varargout{1} = dsdiff;
     end
+    if(ifun==5)
+        [dsdiff] = proc_kerns(rs,drs,dzs,dk);
+        varargout{1} = dsdiff;
+    end
 
 end
 
diff --git a/chunkie/+chnk/+axissymhelm2d/helm_axi_smooth.m b/chunkie/+chnk/+axissymhelm2d/helm_axi_smooth.m
index 740631f7..99af3498 100644
--- a/chunkie/+chnk/+axissymhelm2d/helm_axi_smooth.m
+++ b/chunkie/+chnk/+axissymhelm2d/helm_axi_smooth.m
@@ -68,7 +68,7 @@
     s{6} = kernsdaa;
     s{4} = kernsdkk;
 
-    if (ifun == 1 || ifun==4)
+    if (ifun == 1 || ifun==4 || ifun==5)
         sout.rk = s;
     end
     if (ifun ==2)
@@ -106,6 +106,36 @@
         sout.dk = dk;
     end
 
+
+    if (ifun == 5)
+        r0t = 0;
+        efac = exp(1i*r0t.*rts)./rts.*wle;
+        r2   = rts.^2;
+        k2   = 1i*r0t.*rts;
+        kern2_v = sum(asymint_v(rts,r0t,efac),1);
+        kern2dk_v = sum(asymintdk_v(rts,r0t,efac),1);
+        kern2da_v = sum(asymintda_v(xle,rts,r0t,efac),1);
+        kern2daa_v = sum(asymintdaa_v(xle,rts,r0t,efac,k2,r2),1);
+        kern2dak_v = sum(asymintdak_v(xle,rts,r0t,efac),1);
+        kern2dkk_v = sum(asymintdkk_v(xle,rts,r0t,efac),1);
+
+        ik = {};
+        ik{1} = kern2_v;
+        ik{2} = kern2da_v;
+        ik{3} = 0*kern2dk_v;
+        ik{6} = kern2daa_v;
+        ik{5} = 0*kern2dak_v;
+        ik{4} = 0*kern2dkk_v;
+        sout.k0 = ik;
+
+        dk = {};
+        for ii=1:6
+            dk{ii} = sout.rk{ii}-sout.k0{ii};
+        end
+        sout.dk = dk;
+    end
+
+
 end
 
 function [val] = asymint_v(r,k,efac)
diff --git a/chunkie/+chnk/+axissymhelm2d/kern.m b/chunkie/+chnk/+axissymhelm2d/kern.m
index e4694444..d46446f9 100644
--- a/chunkie/+chnk/+axissymhelm2d/kern.m
+++ b/chunkie/+chnk/+axissymhelm2d/kern.m
@@ -1,4 +1,4 @@
-function submat = kern(zk, srcinfo, targinfo, origin, type, varargin)
+function submat = kern(zks, srcinfo, targinfo, origin, type, varargin)
 %CHNK.AXISSYMHELM2D.KERN axissymmetric Helmholtz layer potential kernels in 2D
 % 
 % Syntax: submat = chnk.axissymhelm2d.kern(zk,srcinfo,targingo,type,htables)
@@ -61,6 +61,13 @@
 [~, ns] = size(src);
 [~, nt] = size(targ);
 
+if (numel(zks) == 1)
+    zk = zks(1);
+else
+    zk1 = zks(1);
+    zk2 = zks(2);
+end
+
 if strcmpi(type, 'd')
     srcnorm = srcinfo.n;
     [~, grad] = chnk.axissymhelm2d.green(zk, src, targ, origin);
@@ -183,6 +190,31 @@
               - hess(:,:,6).*nxsrc.*nytarg + hess(:,:,3).*nysrc.*nytarg;
 end
 
+if strcmpi(type, 'dprime_re_diff')
+
+  targnorm = targinfo.n;
+  srcnorm = srcinfo.n;
+  ifdiff = 2;
+  [~,~,hess] = chnk.axissymhelm2d.green(zk1, src, targ, origin,ifdiff);
+  nxsrc = repmat(srcnorm(1,:),nt,1);
+  nysrc = repmat(srcnorm(2,:),nt,1);
+  nxtarg = repmat((targnorm(1,:)).',1,ns);
+  nytarg = repmat((targnorm(2,:)).',1,ns);
+  submat1 = hess(:,:,4).*nxsrc.*nxtarg - hess(:,:,5).*nysrc.*nxtarg ... 
+              - hess(:,:,6).*nxsrc.*nytarg + hess(:,:,3).*nysrc.*nytarg;
+  ifdiff = 2;
+  [~,~,hess] = chnk.axissymhelm2d.green(zk2, src, targ, origin,ifdiff);
+  nxsrc = repmat(srcnorm(1,:),nt,1);
+  nysrc = repmat(srcnorm(2,:),nt,1);
+  nxtarg = repmat((targnorm(1,:)).',1,ns);
+  nytarg = repmat((targnorm(2,:)).',1,ns);
+  submat2 = hess(:,:,4).*nxsrc.*nxtarg - hess(:,:,5).*nysrc.*nxtarg ... 
+              - hess(:,:,6).*nxsrc.*nytarg + hess(:,:,3).*nysrc.*nytarg;
+
+  submat = submat1-submat2;
+
+end
+
 
 
 if strcmpi(type, 'dprimediff')
diff --git a/chunkie/+chnk/+axissymhelm2d/load_asym_tables_sub.m b/chunkie/+chnk/+axissymhelm2d/load_asym_tables_sub.m
new file mode 100644
index 00000000..3b6c0b48
--- /dev/null
+++ b/chunkie/+chnk/+axissymhelm2d/load_asym_tables_sub.m
@@ -0,0 +1,95 @@
+function asym_tables = load_asym_tables_sub()
+%CHNK.AXISSYMHELM2D.load_asym_tables loads the precomputed tables
+% for axissymmetric Helmholtz kernels
+    p = mfilename('fullpath');
+    dirname = dir([p ,'.m']).folder;
+    fname = [dirname '/asym_helm_data2.mat'];
+    load(fname,'allvs');
+    msizes = size(allvs);
+    ncheb = msizes(1);
+    nks = msizes(5);
+    nas = msizes(4);
+    nkers = msizes(3);
+    xcheb = (0:(ncheb-1))/(ncheb-1)*pi;
+
+    [N,X]=meshgrid(0:(ncheb-1),xcheb);
+    Tc2v = cos(N.*X);
+    Tv2c = inv(Tc2v);
+    Tfull = kron(Tv2c,Tv2c);
+    Tc2vd = N.*sin(N.*X)./sqrt(1-cos(X).^2);
+    Tc2vd(1,:) = (0:(ncheb-1)).^2;
+    Tc2vd(end,:) = (-1).^(1:ncheb).*(0:(ncheb-1)).^2;
+    Tc2cd = Tv2c*Tc2vd;
+
+    xcheb = cos(xcheb)';
+
+
+    allda = allvs;
+    alldk = allvs;
+    alldaa = allvs;
+    alldak = allvs;
+    alldkk = allvs;
+
+    for kk=1:nks
+        for jj=1:nas
+            for ii=1:nkers
+%                vmat = squeeze(allvs(:,:,ii,jj,kk));
+                if (ii >=1)
+                    vmat = squeeze(allvs(:,:,ii,jj,kk));
+                else
+                    vmat = squeeze(allvs(:,:,ii,jj,kk));
+                    vmat = reshape(Tfull*vmat(:),[12,12]);
+                    allvs(:,:,ii,jj,kk) = vmat;
+                end
+                ia = jj;
+                vda = (2)^(ia-1)*20*Tc2cd*vmat;
+                vdaa= (2)^(ia-1)*20*Tc2cd*vda;
+                vdk = 2/(pi)*Tc2cd*vmat.';
+                vdkk= (2)/pi*Tc2cd*vdk;
+                vdak = 2/(pi)*Tc2cd*vda.';
+                allda(:,:,ii,jj,kk) = vda;
+                alldaa(:,:,ii,jj,kk)= vdaa;
+                alldk(:,:,ii,jj,kk) = vdk.';
+                alldak(:,:,ii,jj,kk)= vdak.';
+                alldkk(:,:,ii,jj,kk)= vdkk.';
+            end
+        end
+    end    
+
+    asym_tables =[];
+    asym_tables.allvs = allvs;
+    asym_tables.allda = allda;
+    asym_tables.alldaa= alldaa;
+    asym_tables.alldk = alldk;
+    asym_tables.alldak= alldak;
+    asym_tables.alldkk= alldkk;
+    asym_tables.ncheb = ncheb;
+
+    nlege = 500;
+    [xlege,wlege,~,~] = lege.exps(nlege);
+    xlege = (pi*(xlege+1)/2);
+    wlege = wlege*pi/2;
+    asym_tables.xlege = xlege;
+    asym_tables.wlege = wlege;
+    
+   	nlege = 100;
+    [xlege,wlege,~,~] = lege.exps(nlege);
+    xlege = (pi*(xlege+1)/2);
+    wlege = wlege*pi/2;
+    asym_tables.xlege_mid = xlege;
+    asym_tables.wlege_mid = wlege;
+    
+  	nlege = 50;
+    [xlege,wlege,~,~] = lege.exps(nlege);
+    xlege = (pi*(xlege+1)/2);
+    wlege = wlege*pi/2;
+    asym_tables.xlege_midnear = xlege;
+    asym_tables.wlege_midnear = wlege;
+    
+  	nlege = 20;
+    [xlege,wlege,~,~] = lege.exps(nlege);
+    xlege = (pi*(xlege+1)/2);
+    wlege = wlege*pi/2;
+    asym_tables.xlege_midnearnear = xlege;
+    asym_tables.wlege_midnearnear = wlege;
+end
diff --git a/chunkie/@kernel/axissymhelm2ddiff.m b/chunkie/@kernel/axissymhelm2ddiff.m
index 3c7b1a72..d4cec2b5 100644
--- a/chunkie/@kernel/axissymhelm2ddiff.m
+++ b/chunkie/@kernel/axissymhelm2ddiff.m
@@ -11,8 +11,8 @@
 %
 %  COEFS is optional. If not given then COEFS is set to [1 1].
 %
-%  NOTES: currently only supports coefs = [c c], and zks(1), real, and
-%   zks(2) = I*zks(1)
+%  NOTES: currently only supports coefs = [c c], and zks(1:2) both real or
+%   zks(2) = I*zks(1), with zks(1) real.
 %
 % See also CHNK.AXISSYMHELM2D.KERN.
 
@@ -27,8 +27,12 @@
 zr1 = real(zks(1)); zi1 = imag(zks(1));
 zr2 = real(zks(2)); zi2 = imag(zks(2));
 
+ifreal = 0;
 if zi1 ~=0 || zr2 ~=0 || zr1 ~= zi2
-    error('Wave numbers must be of the form zk, 1i*zk with zk real');
+    ifreal = 1;
+    if (zi1~= 0 || zi2 ~= 0)
+        error('Wave numbers must be of the form zk, 1i*zk with zk real');
+    end
 end
 
 obj = kernel();
@@ -89,11 +93,19 @@
             error('Coefs must be of form [c c]');   
         end
         obj.type = 'dp';
+        if(~ifreal)
         obj.eval = @(s,t) coefs(1)*chnk.axissymhelm2d.kern(real(zks(1)),  ... 
                                                 s, t, [0,0], 'dprimediff');
         obj.shifted_eval = @(s,t,o) coefs(1)*chnk.axissymhelm2d.kern(real(zks(1)),  ... 
                                                 s, t, o, 'dprimediff');
         obj.fmm = [];
+        else
+        obj.eval = @(s,t) coefs(1)*chnk.axissymhelm2d.kern(zks,  ... 
+                                                s, t, [0,0], 'dprime_re_diff');
+        obj.shifted_eval = @(s,t,o) coefs(1)*chnk.axissymhelm2d.kern(zks,  ... 
+                                                s, t, o, 'dprime_re_diff');
+        obj.fmm = [];
+        end
         if ( abs(coefs(1)-coefs(2)) < eps )
             obj.sing = 'log';
         else
diff --git a/chunkie/demo/demo_axissym_neumann.m b/chunkie/demo/demo_axissym_neumann.m
index db289f4e..380a4798 100644
--- a/chunkie/demo/demo_axissym_neumann.m
+++ b/chunkie/demo/demo_axissym_neumann.m
@@ -1,4 +1,4 @@
-addpaths_loc();
+%addpaths_loc();
 clear();
 
 rng(pi)
diff --git a/devtools/test/chunkermat_axissymhelm_transmissionTest.m b/devtools/test/chunkermat_axissymhelm_transmissionTest.m
index 82ef0deb..a699359b 100644
--- a/devtools/test/chunkermat_axissymhelm_transmissionTest.m
+++ b/devtools/test/chunkermat_axissymhelm_transmissionTest.m
@@ -2,12 +2,13 @@
 iseed = 8675309;
 rng(iseed);
 
-addpaths_loc();
+%addpaths_loc();
 
 zk = 5.1;
 
-% type = 'cgrph';
-type = 'chnkr-torus';
+type = 'chnkr-star';
+% type = 'chnkr-torus';
+%type = 'cgrph';
 
 pref = [];
 pref.k = 16;
@@ -19,6 +20,8 @@
 
 [chnkr, sources, targets] = get_geometry(type, pref, ns, nt, maxchunklen);
 wts = chnkr.wts; wts = wts(:);
+sources(1,:) = abs(sources(1,:));
+targets(1,:) = abs(targets(1,:));
 
 l2scale = false;
 fprintf('Done building geometry\n');
diff --git a/devtools/test/chunkermat_axissymhelm_transmissionTest_re.m b/devtools/test/chunkermat_axissymhelm_transmissionTest_re.m
new file mode 100644
index 00000000..586c7ad3
--- /dev/null
+++ b/devtools/test/chunkermat_axissymhelm_transmissionTest_re.m
@@ -0,0 +1,355 @@
+clearvars; close all;
+iseed = 8675309;
+rng(iseed);
+
+%addpaths_loc();
+
+zk = 10.1;
+
+type = 'chnkr-star';
+%type = 'chnkr-torus';
+type = 'cgrph';
+
+pref = [];
+pref.k = 16;
+ns = 10;
+nt = 10;
+ppw = 80;   % points per wavelength;
+maxchunklen = pref.k/ppw/real(zk)*2*pi;
+maxchunklen = 0.5;
+
+[chnkr, sources, targets] = get_geometry(type, pref, ns, nt, maxchunklen);
+wts = chnkr.wts; wts = wts(:);
+
+l2scale = false;
+fprintf('Done building geometry\n');
+
+% source strengths
+strengths_in = randn(ns, 1);
+strengths_out = randn(nt, 1);
+
+
+% targets
+sources(1,:) = abs(sources(1,:));
+targets(1,:) = abs(targets(1,:));
+sr = sources;
+sources = targets;
+targets = sr;
+
+% Plot everything
+
+figure(1)
+clf
+hold off
+plot(chnkr)
+hold on
+scatter(sources(1,:), sources(2,:), 'o')
+scatter(targets(1,:), targets(2,:), 'x')
+axis equal
+
+
+
+% For solving the transmission boundary value problem, we
+% use the repesentation 
+% u_{i} = D_{k_{i}} [\sigma] - S_{k_{i}}[\tau]
+% 
+% and impose jumps in u and du/dn as exterior-interior respectively
+%
+
+% Set up kernels
+Skdiff     = kernel('axissymhelmdiff', 's', [zk zk/2]);
+Dkdiff     = kernel('axissymhelmdiff', 'd', [zk zk/2]);
+Skpdiff    = kernel('axissymhelmdiff', 'sprime', [zk zk/2]);
+Dkpdiff    = kernel('axissymhelmdiff', 'dprime', [zk zk/2]);
+
+K = [Dkdiff  -1*Skdiff;
+     Dkpdiff    -1*Skpdiff];
+K = kernel(K);
+Dk    = kernel('axissymhelm', 'd', zk);
+Sk    = kernel('axissymhelm', 's', zk);
+Skp    = kernel('axissymhelm', 'sprime', zk);
+
+Dik    = kernel('axissymhelm', 'd', zk/2);
+Sik    = kernel('axissymhelm', 's', zk/2);
+Sikp    = kernel('axissymhelm', 'sprime', zk/2);
+
+Kouteval =  kernel([Dk, -1*Sk]);
+Kineval =  kernel([Dik, -1*Sik]);
+
+
+% Set up boundary data
+
+srcinfo  = []; srcinfo.r = sources; 
+targinfo = []; targinfo.r = chnkr.r(:,:); targinfo.n = chnkr.n(:,:);
+kernmats = Sk.eval(srcinfo, targinfo);
+kernmatsp = Skp.eval(srcinfo, targinfo);
+ubdry = kernmats*strengths_in;
+dudnbdry = kernmatsp*strengths_in;
+
+
+
+srcinfo  = []; srcinfo.r = targets; 
+targinfo = []; targinfo.r = chnkr.r(:,:); targinfo.n = chnkr.n(:,:);
+kernmats = Sik.eval(srcinfo, targinfo);
+kernmatsp = Sikp.eval(srcinfo, targinfo);
+ubdry = ubdry - kernmats*strengths_out;
+dudnbdry = dudnbdry - kernmatsp*strengths_out;
+
+
+npts = chnkr.npt;
+nsys = K.opdims(1)*npts;
+rhs = zeros(nsys, 1);
+
+
+if(l2scale)
+    rhs(1:2:end) = ubdry.*sqrt(wts);
+    rhs(2:2:end) = dudnbdry.*sqrt(wts);
+else
+    rhs(1:2:end) = ubdry;
+    rhs(2:2:end) = dudnbdry;
+end
+
+% Form matrix
+opts = [];
+opts.l2scale = l2scale;
+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_out = kernmatstarg*strengths_in;
+
+% Compute exact solution
+srcinfo  = []; srcinfo.r  = targets;
+targinfo = []; targinfo.r = sources;
+kernmatstarg = Sik.eval(srcinfo, targinfo);
+utarg_in = kernmatstarg*strengths_out;
+
+
+
+% Compute solution using chunkerkerneval
+% evaluate at targets and compare
+
+opts.forceadap = true;
+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_out = chunkerkerneval(chnkrtotal, Kouteval, sol, targets, opts);
+Dsol_in =  chunkerkerneval(chnkrtotal, Kineval, sol, sources, 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_out-Dsol_out) / (sqrt(chnkrtotal.nch)*norm(utarg_out));
+relerr2 = norm(utarg_out-Dsol_out, 'inf') / dot(abs(sol(:)), wchnkr(:));
+fprintf('relative frobenius error in exterior %5.2e\n', relerr);
+fprintf('relative l_inf/l_1 error in exterior %5.2e\n', relerr2);
+
+relerr  = norm(utarg_in-Dsol_in) / (sqrt(chnkrtotal.nch)*norm(utarg_in));
+relerr2 = norm(utarg_in-Dsol_in, 'inf') / dot(abs(sol(:)), wchnkr(:));
+fprintf('relative frobenius error in interior %5.2e\n', relerr);
+fprintf('relative l_inf/l_1 error in interior %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/(2*pi);
+    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;
+
+    chnkobj = chunkgraph(verts, edge2verts, fchnks, cparams, pref);
+    chnkobj = balance(chnkobj);
+    
+    ts = 0.0+2*pi*rand(ns,1);
+    sources = 3.0*[cos(ts)';sin(ts)'];
+    
+    ts = 0.0+2*pi*rand(nt,1);
+    targets = 0.2*[cos(ts)'; sin(ts)'];
+
+
+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 = sources .* (0.7 - 0.4.*repmat(rand(1,ns), 2, 1));
+
+    ts = -pi/2 + pi*rand(nt, 1);
+    targets = starfish(ts, narms, amp);
+    targets = targets .* (1.2 + 3.0*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 = sources .* (0.7 - 0.4.*repmat(rand(1,ns), 2, 1));
+    sources(1,:) = sources(1,:) + ctr(1);
+
+    ts = -pi/2 + pi*rand(nt, 1);
+    targets = starfish(ts, narms, amp);
+    targets = targets .* (1.2 + 3.0*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
+