From dfa384744de746ef7533edf83c7fde4e3150ba15 Mon Sep 17 00:00:00 2001
From: Mike O'Neil <oneil@cims.nyu.edu>
Date: Tue, 3 Dec 2024 15:26:18 -0500
Subject: [PATCH] testing file for modal green's functions

---
 chunkie/+chnk/+axissymlap2d/g0funcall.m | 356 +++++++++++++-----------
 chunkie/+chnk/+axissymlap2d/gfunc.m     |  29 +-
 chunkie/+chnk/+axissymlap2d/gfuncall.m  |   2 +-
 chunkie/+chnk/+axissymlap2d/gkm_brute.m |  31 +++
 chunkie/+chnk/+axissymlap2d/green.m     |  42 +--
 chunkie/+chnk/+axissymlap2d/qleg_half.m |   1 +
 devtools/test/g0funcallTest.m           |  75 +++++
 7 files changed, 345 insertions(+), 191 deletions(-)
 create mode 100644 chunkie/+chnk/+axissymlap2d/gkm_brute.m
 create mode 100644 devtools/test/g0funcallTest.m

diff --git a/chunkie/+chnk/+axissymlap2d/g0funcall.m b/chunkie/+chnk/+axissymlap2d/g0funcall.m
index 610316a..a95e69f 100644
--- a/chunkie/+chnk/+axissymlap2d/g0funcall.m
+++ b/chunkie/+chnk/+axissymlap2d/g0funcall.m
@@ -1,4 +1,5 @@
 function [gvals, gdzs, gdrs, gdrps] = g0funcall(r, rp, dr, z, zp, dz, maxm)
+
 %
 % chnk.axissymlap2d.g0funcall evaluates a collection of axisymmetric Laplace
 % Green's functions, defined by the expression:
@@ -12,12 +13,13 @@
 % The above scaling should be consistent with what is in
 % chnk.axissymlap2d.gfunc, which is for merely the zero-mode
 % 
-
+    
     twopi = 2*pi;
     fourpi = 4*pi;
     done = 1.0;
     ima = 1i;
 
+    
 
     %r = targ(1)
     %z = targ(2)
@@ -28,21 +30,22 @@
     x = 1/alpha;
     xminus = (dr*dr + dz*dz)/2/r/r0;
 
-    dxdr = (r^2 - r0^2 - (dz)^2)/2/r0/r^2
-    dxdz = 2*(dz)/2/r/r0
-    dxdr0 = (r0^2 - r^2 - (dz)^2)/2/r/r0^2
-    dxdz0 = -dxdz
+    dxdr = (r^2 - r0^2 - (dz)^2)/2/r0/r^2;
+    dxdz = 2*(dz)/2/r/r0;
+    dxdr0 = (r0^2 - r^2 - (dz)^2)/2/r/r0^2;
+    dxdz0 = -dxdz;
   
     %!
     %! if xminus is very small, use the forward recurrence
     %!
 
     %!!!!print *, 'inside g0mall, x = ', x
-  
-  
-    iffwd = 0
+
+
+    
+    iffwd = 0;
     if (x < 1.005)
-        iffwd = 1
+        iffwd = 1;
         if ((x >= 1.0005d0) && (maxm > 163))
             iffwd = 0;
         end
@@ -65,181 +68,222 @@
         
     end
 
-  
-    if (iffwd = 1)
 
-        %!!xminus = .000005d0
-        %!!x = 1 + xminus
+    %iffwd
     
-        call axi_q2lege01(x, xminus, q0, q1, dq0, dq1)
 
-        done = 1
-        half = done/2
-        %pi = 4*atan(done)
+    gvals = zeros(maxm+1,1);
+    gdzs = zeros(maxm+1,1);
+    gdrs = zeros(maxm+1,1);
+    gdrps = zeros(maxm+1,1);
     
-        fac = done/sqrt(r*r0)/4/pi^2
-        vals(0) = q0
-        vals(1) = q1
+    if (iffwd == 1)
 
-        derprev = dq0
-        der = dq1
+        [q0, q1, dq0] = chnk.axissymlap2d.qleg_half(xminus);
+        dq1 = (-q0 + x*q1)/2/(x+1)/xminus;
+        
+        %call axi_q2lege01(x, xminus, q0, q1, dq0, dq1)
+
+        half = done/2;
+    
+        %fac = done/sqrt(r*r0)/4/pi^2;
+        gvals(1) = 2*pi*sqrt(rp/r)*q0;
+        gvals(2) = 2*pi*sqrt(rp/r)*q1;
 
-        grads(1,0) = (dq0*dxdr - q0/2/r)
-        grads(1,1) = (dq1*dxdr - q1/2/r)
+        %derprev = dq0
+        %der = dq1
+
+        
+        
+        %grads(1,0) = (dq0*dxdr - q0/2/r)
+        %grads(1,1) = (dq1*dxdr - q1/2/r)
         
-        grads(2,0) = dq0*dxdz
-        grads(2,1) = dq1*dxdz
+        %grads(2,0) = dq0*dxdz
+        %grads(2,1) = dq1*dxdz
         
-        grad0s(1,0) = (dq0*dxdr0 - q0/2/r0)
-        grad0s(1,1) = (dq1*dxdr0 - q1/2/r0)
+        %grad0s(1,0) = (dq0*dxdr0 - q0/2/r0)
+        %grad0s(1,1) = (dq1*dxdr0 - q1/2/r0)
         
-        grad0s(2,0) = dq0*dxdz0
-        grad0s(2,1) = dq1*dxdz0
+        %grad0s(2,0) = dq0*dxdz0
+        %grad0s(2,1) = dq1*dxdz0
         
         for i = 1:(maxm-1)
-            vals(i+1) = (2*i*x*vals(i) - (i-half)*vals(i-1))/(i+half)
-            dernext = (2*i*(vals(i)+x*der) - (i-half)*derprev)/(i+half)
-            grads(1,i+1) = (dernext*dxdr - vals(i+1)/2/r)
-            grads(2,i+1) = dernext*dxdz
-            grad0s(1,i+1) = (dernext*dxdr0 - vals(i+1)/2/r0)
-            grad0s(2,i+1) = dernext*dxdz0
-            derprev = der
-            der = dernext
-            %! if (abs(vals(i+1)) > abs(vals(i))) then
-            %!   print *
-            %!   print *
-            %!   call prin2('x = *', x, 1)
-            %!   call prinf('i = *', i, 1)
-            %!   call prin2('vals(i+1) = *', vals(i+1), 1)
-            %!   call prin2('vals(i) = *', vals(i), 1)
-            %!   stop
-            %! endif
+            gvals(i+2) = (2*i*x*gvals(i+1) - (i-half)*gvals(i))/(i+half);
+            %dernext = (2*i*(vals(i)+x*der) - (i-half)*derprev)/(i+half)
+            %grads(1,i+1) = (dernext*dxdr - vals(i+1)/2/r)
+            %grads(2,i+1) = dernext*dxdz
+            %grad0s(1,i+1) = (dernext*dxdr0 - vals(i+1)/2/r0)
+            %grad0s(2,i+1) = dernext*dxdz0
+            %derprev = der
+            %der = dernext
         end
 
-        for i = 0:maxm
-            vals(i) = vals(i)*fac
-            grads(1,i) = grads(1,i)*fac
-            grads(2,i) = grads(2,i)*fac
-            grad0s(1,i) = grad0s(1,i)*fac
-            grad0s(2,i) = grad0s(2,i)*fac
-            vals(-i) = vals(i)
-            grads(1,-i) = grads(1,i)
-            grads(2,-i) = grads(2,i)
-            grad0s(1,-i) = grad0s(1,i)
-            grad0s(2,-i) = grad0s(2,i)
+        for i = 1:(maxm+1)
+            %gvals(i) = gvals(i)*fac;
+            %grads(1,i) = grads(1,i)*fac
+            %grads(2,i) = grads(2,i)*fac
+            %grad0s(1,i) = grad0s(1,i)*fac
+            %grad0s(2,i) = grad0s(2,i)*fac
+            %vals(-i) = vals(i)
+            %grads(1,-i) = grads(1,i)
+            %grads(2,-i) = grads(2,i)
+            %grad0s(1,-i) = grad0s(1,i)
+            %grad0s(2,-i) = grad0s(2,i)
         end
 
+        %gvals
+
+
+
+        % check the scaling here correctly
+
+        %mode = 5;
+        
+        % m = 4000;
+        % dint = 0;
+        % h = 2*pi/m;
+        % for i = 1:m
+        %     s = (i-1)*h;
+        %     x = r;
+        %     y = 0;
+        %     xp = rp*cos(s);
+        %     yp = rp*sin(s);
+        %     dist = sqrt( (x-xp)^2 + (y-yp)^2 + (z-zp)^2 );
+        %     dint = dint + exp(-ima*mode*s)*h/dist;
+        % end
+
+        %dint = dint*pi*rp;
+        
+        %zk = 0;
+        %nq = 2000;
+        %dint = chnk.axissymlap2d.gkm_brute(r, rp, z, zp, zk, mode, nq);
+        %dint = dint*4*pi*pi*rp;
+
+        %error = dint-dint2
+    
+        %disp(['gval = ' num2str(gvals(mode+1))]);
+        %disp(['dint = ' num2str(dint)]);
+        %disp(['error = ' num2str(abs((gvals(mode+1) - dint)/dint))]);
+        %disp(['ratio = ' num2str(gvals(mode+1)/dint)]);
+        
+
+
+        
         return
     end
 
-    %!
-    %! if here, xminus > .005, so run forward and backward recurrence
-    %!
 
-    %!
-    %! run the recurrence, starting from maxm, until it has exploded
-    %! for BOTH the values and derivatives
-    %!
-    done = 1
-    half = done/2
-    f = 1
-    fprev = 0
-    der = 1
-    derprev = 0
-    maxiter = 100000
-    upbound = 1.0d19
+
+
+    
+  %   %!
+  %   %! if here, xminus > .005, so run forward and backward recurrence
+  %   %!
+
+  %   %!
+  %   %! run the recurrence, starting from maxm, until it has exploded
+  %   %! for BOTH the values and derivatives
+  %   %!
+  %   done = 1
+  %   half = done/2
+  %   f = 1
+  %   fprev = 0
+  %   der = 1
+  %   derprev = 0
+  %   maxiter = 100000
+  %   upbound = 1.0d19
   
-    for i = maxm:maxiter
-        fnext = (2*i*x*f - (i-half)*fprev)/(i+half)
-        dernext = (2*i*(x*der+f) - (i-half)*derprev)/(i+half)
-        if (abs(fnext) .ge. upbound) 
-            if (abs(dernext) .ge. upbound)
-                nterms = i+1
-                exit
-            end
-        end
-        fprev = f
-        f = fnext
-        derprev = der
-        der = dernext
-    end
+  %   for i = maxm:maxiter
+  %       fnext = (2*i*x*f - (i-half)*fprev)/(i+half)
+  %       dernext = (2*i*(x*der+f) - (i-half)*derprev)/(i+half)
+  %       if (abs(fnext) .ge. upbound) 
+  %           if (abs(dernext) .ge. upbound)
+  %               nterms = i+1
+  %               exit
+  %           end
+  %       end
+  %       fprev = f
+  %       f = fnext
+  %       derprev = der
+  %       der = dernext
+  %   end
     
-    %!
-    %! now start at nterms and recurse down to maxm
-    %!
+  %   %!
+  %   %! now start at nterms and recurse down to maxm
+  %   %!
 
-    if (nterms .lt. 10) then
-        nterms = 10
-    end
+  %   if (nterms .lt. 10) then
+  %       nterms = 10
+  %   end
 
-    fnext = 0
-    f = 1
-    dernext = 0
-    der = 1
+  %   fnext = 0
+  %   f = 1
+  %   dernext = 0
+  %   der = 1
 
 
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%555
-    % Make correct downwrd recurrence
+  %   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
+  %   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%555
+  %   % Make correct downwrd recurrence
     
-  for i = nterm,maxm,-1
-    fprev = (2*i*x*f - (i+half)*fnext)/(i-half)
-    fnext = f
-    f = fprev
-    derprev = (2*i*(x*der+f) - (i+half)*dernext)/(i-half)
-    dernext = der
-    der = derprev
-  enddo
-
-  vals(maxm-1) = f
-  vals(maxm) = fnext
-
-  ders(maxm-1) = der
-  ders(maxm) = dernext
+  % for i = nterm,maxm,-1
+  %   fprev = (2*i*x*f - (i+half)*fnext)/(i-half)
+  %   fnext = f
+  %   f = fprev
+  %   derprev = (2*i*(x*der+f) - (i+half)*dernext)/(i-half)
+  %   dernext = der
+  %   der = derprev
+  % enddo
+
+  % vals(maxm-1) = f
+  % vals(maxm) = fnext
+
+  % ders(maxm-1) = der
+  % ders(maxm) = dernext
   
-  do i = maxm-1,1,-1
-    vals(i-1) = (2*i*x*vals(i) - (i+half)*vals(i+1))/(i-half)
-    ders(i-1) = (2*i*(x*ders(i)+vals(i)) &
-        - (i+half)*ders(i+1))/(i-half)
-  end do
-
-
-  !
-  ! normalize the values, and use a formula for the derivatives
-  !
-  call axi_q2lege01(x, xminus, q0, q1, dq0, dq1)
-
-  done = 1
-  pi = 4*atan(done)
-  ratio = q0/vals(0)
-
-  do i = 0,maxm
-    vals(i) = vals(i)*ratio
-  enddo
-
-  ders(0) = dq0
-  ders(1) = dq1
-  do i = 2,maxm
-    ders(i) = -(i-.5d0)*(vals(i-1) - x*vals(i))/(1+x)/xminus
-  end do
-
-  !
-  ! and scale everyone...
-  !
-  fac = 1/sqrt(r*r0)/4/pi^2
+  % do i = maxm-1,1,-1
+  %   vals(i-1) = (2*i*x*vals(i) - (i+half)*vals(i+1))/(i-half)
+  %   ders(i-1) = (2*i*(x*ders(i)+vals(i)) &
+  %       - (i+half)*ders(i+1))/(i-half)
+  % end do
+
+
+  % !
+  % ! normalize the values, and use a formula for the derivatives
+  % !
+  % call axi_q2lege01(x, xminus, q0, q1, dq0, dq1)
+
+  % done = 1
+  % pi = 4*atan(done)
+  % ratio = q0/vals(0)
+
+  % do i = 0,maxm
+  %   vals(i) = vals(i)*ratio
+  % enddo
+
+  % ders(0) = dq0
+  % ders(1) = dq1
+  % do i = 2,maxm
+  %   ders(i) = -(i-.5d0)*(vals(i-1) - x*vals(i))/(1+x)/xminus
+  % end do
+
+  % !
+  % ! and scale everyone...
+  % !
+  % fac = 1/sqrt(r*r0)/4/pi^2
   
-  do i = 0,maxm
-    grads(1,i) = (ders(i)*dxdr - vals(i)/2/r)*fac
-    grads(2,i) = ders(i)*dxdz*fac
-    grad0s(1,i) = (ders(i)*dxdr0 - vals(i)/2/r0)*fac
-    grad0s(2,i) = ders(i)*dxdz0*fac
-    vals(i) = vals(i)*fac
-    vals(-i) = vals(i)
-    grads(1,-i) = grads(1,i)
-    grads(2,-i) = grads(2,i)
-    grad0s(1,-i) = grad0s(1,i)
-    grad0s(2,-i) = grad0s(2,i)
-  end do
+  % do i = 0,maxm
+  %   grads(1,i) = (ders(i)*dxdr - vals(i)/2/r)*fac
+  %   grads(2,i) = ders(i)*dxdz*fac
+  %   grad0s(1,i) = (ders(i)*dxdr0 - vals(i)/2/r0)*fac
+  %   grad0s(2,i) = ders(i)*dxdz0*fac
+  %   vals(i) = vals(i)*fac
+  %   vals(-i) = vals(i)
+  %   grads(1,-i) = grads(1,i)
+  %   grads(2,-i) = grads(2,i)
+  %   grad0s(1,-i) = grad0s(1,i)
+  %   grad0s(2,-i) = grad0s(2,i)
+  % end do
     
 
 
diff --git a/chunkie/+chnk/+axissymlap2d/gfunc.m b/chunkie/+chnk/+axissymlap2d/gfunc.m
index cd58942..1381c88 100644
--- a/chunkie/+chnk/+axissymlap2d/gfunc.m
+++ b/chunkie/+chnk/+axissymlap2d/gfunc.m
@@ -1,9 +1,10 @@
-function [gval,gdz,gdr,gdrp] = gfunc(r, rp, dr, z, zp, dz)
+function [gval, gdz, gdr, gdrp] = gfunc(r, rp, dr, z, zp, dz)
+
 %
 % chnk.axissymlap2d.gfunc evaluates the zeroth order axisymmetric Laplace
 % Green's funcion, defined by the expression:
 %
-%     gfunc = pi * rp * \int_0^{2\pi} 1/|x - x'| d\theta
+%     gfunc = pi * rp * \int_0^{2\pi} 1/|x - x'| d\theta'
 %
 % The extra factor of rp out front makes subsequent interfacing with RCIP
 % slightly easier
@@ -11,23 +12,23 @@
 % 
     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;
-
+    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;
-
+    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;
-
-
+    gdr  = domega3*(rp./r).^((n-2)/2)./(rp.*r) ...
+           .*rfac;
+    
+    
     % % check the scaling here correctly
     % m = 1000;
     % dint = 0;
diff --git a/chunkie/+chnk/+axissymlap2d/gfuncall.m b/chunkie/+chnk/+axissymlap2d/gfuncall.m
index ac7ce93..657832a 100644
--- a/chunkie/+chnk/+axissymlap2d/gfuncall.m
+++ b/chunkie/+chnk/+axissymlap2d/gfuncall.m
@@ -1,6 +1,6 @@
 function [gvals, gdzs, gdrs, gdrps] = gfuncall(zk, r, rp, dr, z, zp, dz, maxm)
 %
-% chnk.axissymlap2d.gfuncall evaluates a collection of axisymmetric Laplace
+% chnk.axissymlap2d.gfuncall evaluates a collection of axisymmetric Helmholtz
 % Green's functions, defined by the expression:
 %
 %     gfunc(n) = pi*rp * \int_0^{2\pi} e^(i*k*|x-x'|)/|x - x'| e^(-i n t) dt 
diff --git a/chunkie/+chnk/+axissymlap2d/gkm_brute.m b/chunkie/+chnk/+axissymlap2d/gkm_brute.m
new file mode 100644
index 0000000..698f0d5
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/gkm_brute.m
@@ -0,0 +1,31 @@
+function gval = gkm_brute(r, rp, z, zp, zk, mode, n)
+
+%
+% This routine evaluates the modal Green's function of mode using the
+% trapezoidal rule with n points:
+%
+%     gval = \int_0^{2\pi} e^{i k |x-x'|} / (4 pi |x-x'|) e^(-i mode t) dt 
+%
+% It is assumed that
+%     (x,y,z) = (r,0,z) in cylindrical coordinates
+%     (x',y',z') = (rp*cos(t), rp*sin(t), z') in cylindrical coordinates
+%
+    ima = 1j;
+    
+    gval = 0;
+    h = 2*pi/n;
+    for i = 1:n
+        s = (i-1)*h;
+        x = r;
+        y = 0;
+        xp = rp*cos(s);
+        yp = rp*sin(s);
+        dist = sqrt( (x-xp)^2 + (y-yp)^2 + (z-zp)^2 );
+        gval = gval + exp(-ima*mode*s)*h/dist;
+    end
+
+    gval = gval/4/pi;
+    
+
+end
+
diff --git a/chunkie/+chnk/+axissymlap2d/green.m b/chunkie/+chnk/+axissymlap2d/green.m
index 564d6ab..89fa498 100644
--- a/chunkie/+chnk/+axissymlap2d/green.m
+++ b/chunkie/+chnk/+axissymlap2d/green.m
@@ -1,5 +1,6 @@
 function [val, grad] = green(src, targ, origin)
-%CHNK.AXISSYMHELM2D.GREEN evaluate the Laplace green's function
+%
+% 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.
@@ -12,25 +13,26 @@
 % 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);
+    [~, ns] = size(src);
+    [~, nt] = size(targ);
 
-val = g;
-gtmp(:,:,1) = gdr;
-gtmp(:,:,2) = gdrp;
-gtmp(:,:,3) = gdz;
-grad = gtmp;
+    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;
+    
 
+end
diff --git a/chunkie/+chnk/+axissymlap2d/qleg_half.m b/chunkie/+chnk/+axissymlap2d/qleg_half.m
index df1d1d5..e33acc7 100644
--- a/chunkie/+chnk/+axissymlap2d/qleg_half.m
+++ b/chunkie/+chnk/+axissymlap2d/qleg_half.m
@@ -1,4 +1,5 @@
 function [q0,q1,q0d] = qleg_half(t)
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
 %   caveat utilitor: this function evaluates Q_{-1/2} (and its
diff --git a/devtools/test/g0funcallTest.m b/devtools/test/g0funcallTest.m
new file mode 100644
index 0000000..c9b5028
--- /dev/null
+++ b/devtools/test/g0funcallTest.m
@@ -0,0 +1,75 @@
+
+
+%
+% verify that the modeal Green's functions for Laplace are working properly
+%
+
+h = 0.01;
+
+r = 1.0;
+rp = r + h;
+dr = -h;
+
+z = 0.5;
+zp = z-h;
+dz = h;
+
+maxm = 30;
+
+exact = zeros(maxm+1,1);
+exact_gdz = zeros(maxm+1,1);
+exact_gdr = zeros(maxm+1,1);
+exact_gdrp = zeros(maxm+1,1);
+
+zk = 0;
+nq = 2000;
+hhh = 0.00001
+
+for m = 1:(maxm+1)
+    dint = chnk.axissymlap2d.gkm_brute(r, rp, z, zp, zk, m-1, nq);
+    exact(m) = dint*4*pi*pi*rp;
+
+    % calculate z derivative using finite differences
+    dint2 = chnk.axissymlap2d.gkm_brute(r, rp, z+hhh, zp, zk, m-1, nq);
+    dint1 = chnk.axissymlap2d.gkm_brute(r, rp, z-hhh, zp, zk, m-1, nq);
+    exact_gdz(m) = (dint2-dint1)/(2*hhh)*4*pi*pi*rp;
+
+    % calculate r derivative using finite differences
+    dint2 = chnk.axissymlap2d.gkm_brute(r+hhh, rp, z, zp, zk, m-1, nq);
+    dint1 = chnk.axissymlap2d.gkm_brute(r-hhh, rp, z, zp, zk, m-1, nq);
+    exact_gdr(m) = (dint2-dint1)/(2*hhh)*4*pi*pi*rp;
+
+    % calculate rp derivative using finite differences
+    dint2 = chnk.axissymlap2d.gkm_brute(r, rp+hhh, z, zp, zk, m-1, nq);
+    dint1 = chnk.axissymlap2d.gkm_brute(r, rp-hhh, z, zp, zk, m-1, nq);
+    exact_gdrp(m) = (dint2-dint1)/(2*hhh)*4*pi*pi*rp;
+
+end
+
+exact
+
+exact_gdz
+
+exact_gdr
+
+exact_gdrp
+
+
+%
+% compare with routine for just zero mode
+%
+[gval, gdz, gdr, gdrp] = chnk.axissymlap2d.gfunc(r, rp, dr, z, zp, dz);
+
+disp(['from 0 mode gfunc, gval = ' num2str(gval) ' error = ' num2str(abs(gval-exact(1)))]);
+disp(['from 0 mode gfunc, gdz = ' num2str(gdz) ' error = ' num2str(abs(gdz-exact_gdz(1)))]);
+disp(['from 0 mode gfunc, gdr = ' num2str(gdr) ' error = ' num2str(abs(gdr-exact_gdr(1)))]);
+disp(['from 0 mode gfunc, gdrp = ' num2str(gdrp) ' error = ' num2str(abs(gdrp-exact_gdrp(1)))]);
+
+
+
+
+% evaluate all modes now
+[gvals, gdzs, gdrs, gdrps] = chnk.axissymlap2d.g0funcall(r, rp, dr, z, zp, dz, maxm);
+disp('errors in all modes')
+errors = abs(gvals-exact)
+