From e7c58593cea0f5b57b3dd5babfe6367abf9f1bc2 Mon Sep 17 00:00:00 2001
From: Manas Rachh <mrachh@flatironinstitute.org>
Date: Tue, 2 Apr 2024 23:06:01 -0400
Subject: [PATCH 1/7] fixing more chnkr.h issues and adding test for
 arclengthfun

---
 chunkie/+chnk/+quadba/buildmattd.m |  1 -
 chunkie/+chnk/funcuni.m            |  8 ++++---
 chunkie/@chunker/arclengthdens.m   |  5 ++---
 chunkie/@chunker/arclengthfun.m    |  5 ++---
 chunkie/@chunker/move.m            |  2 +-
 chunkie/@chunker/nearest.m         |  1 -
 chunkie/@chunker/whts.m            |  8 +++----
 chunkie/chunkerkernevalmat.m       |  3 +--
 chunkie/chunkermat.m               |  3 +--
 devtools/test/arclengthfunTest.m   | 36 ++++++++++++++++++++++++++++++
 10 files changed, 52 insertions(+), 20 deletions(-)
 create mode 100644 devtools/test/arclengthfunTest.m

diff --git a/chunkie/+chnk/+quadba/buildmattd.m b/chunkie/+chnk/+quadba/buildmattd.m
index dfe942e..07069d5 100644
--- a/chunkie/+chnk/+quadba/buildmattd.m
+++ b/chunkie/+chnk/+quadba/buildmattd.m
@@ -11,7 +11,6 @@
 adj = chnkr.adj;
 d = chnkr.d;
 d2 = chnkr.d2;
-h = chnkr.h;
 
 [~,~,u] = lege.exps(k);
 
diff --git a/chunkie/+chnk/funcuni.m b/chunkie/+chnk/funcuni.m
index ac93bd3..b11dffb 100644
--- a/chunkie/+chnk/funcuni.m
+++ b/chunkie/+chnk/funcuni.m
@@ -104,16 +104,18 @@
 for i = 1:nch
     a=ab(1,i);
     b=ab(2,i);
+    h = (b-a)/2;
     
     ts = a + (b-a)*(xs+1)/2;
     [out{:}] = fcurve(ts);
     chnkr.r(:,:,i) = reshape(out{1},dim,k);
-    chnkr.d(:,:,i) = reshape(out{2},dim,k);
-    chnkr.d2(:,:,i) = reshape(out{3},dim,k);
-    chnkr.h(i) = (b-a)/2;
+    chnkr.d(:,:,i) = reshape(out{2},dim,k)*h;
+    chnkr.d2(:,:,i) = reshape(out{3},dim,k)*h*h;
 end
 
 chnkr.adj = adjs(:,1:nch);
+chnkr.wtsstor(:,1:nch) = weights(chnkr);
+chnkr.nstor(:,:,1:nch) = normals(chnkr);
 
 end
 
diff --git a/chunkie/@chunker/arclengthdens.m b/chunkie/@chunker/arclengthdens.m
index cecc8d2..470e567 100644
--- a/chunkie/@chunker/arclengthdens.m
+++ b/chunkie/@chunker/arclengthdens.m
@@ -1,10 +1,9 @@
 function ds = arclengthdens(chnkr)
 %ARCLENGTHDENS arc length density on chunks
 %
-% warning: this takes the length of the ith chunk in parameter space to 
-%   be 2*chnkr.h(i) as opposed to 2. Thus the smooth integration rule
+%   The smooth integration rule
 %   on the ith chunk is 
-%       ds(:,i).*w*chnkr.h(i)
+%       ds(:,i).*w
 %   where w are the standard Legendre weights of appropriate order and 
 %   ds is the output of this routine
 %
diff --git a/chunkie/@chunker/arclengthfun.m b/chunkie/@chunker/arclengthfun.m
index 092da6b..63f7724 100644
--- a/chunkie/@chunker/arclengthfun.m
+++ b/chunkie/@chunker/arclengthfun.m
@@ -13,10 +13,9 @@
 istart = 1;
 
 A = lege.intmat(chnkr.k);
-[~,w] = lege.exps(chnkr.k);
 ds = arclengthdens(chnkr);
-chunklens = sum((w(:)*(chnkr.h(:).')).*ds,1);
-s = (A*ds).*(chnkr.h(:).');
+chunklens = sum(chnkr.wts,1);
+s = A*ds;
 
 for i = 1:ncomp
     nch = nchs(i);
diff --git a/chunkie/@chunker/move.m b/chunkie/@chunker/move.m
index 8cbe34a..9484aef 100644
--- a/chunkie/@chunker/move.m
+++ b/chunkie/@chunker/move.m
@@ -19,4 +19,4 @@
 
 obj.wts = weights(obj);
  
-end
\ No newline at end of file
+end
diff --git a/chunkie/@chunker/nearest.m b/chunkie/@chunker/nearest.m
index d585f3e..7f1bbdb 100644
--- a/chunkie/@chunker/nearest.m
+++ b/chunkie/@chunker/nearest.m
@@ -71,7 +71,6 @@
     r = chnkr.rstor(:,:,ii);
     d = chnkr.dstor(:,:,ii);
     d2 = chnkr.d2stor(:,:,ii);
-    h = chnkr.hstor(ii);
     rc = u*(r.');
     dc = u*(d.');
     d2c = u*(d2.');
diff --git a/chunkie/@chunker/whts.m b/chunkie/@chunker/whts.m
index 1da1df0..3a7fc54 100644
--- a/chunkie/@chunker/whts.m
+++ b/chunkie/@chunker/whts.m
@@ -1,12 +1,12 @@
-function wchnk = whts(chnkr)
+function wts = whts(chnkr)
 
 
 warning('whts is deprecated and will be removed, use weights instead');
 
   k = chnkr.k;
   nch = chnkr.nch;
-  [~,w] = lege.exps(k);
-  wchnk = reshape(sqrt(sum((chnkr.d).^2,1)),k,nch);
-  wchnk = wchnk.*bsxfun(@times,w(:),(chnkr.h(:)).');
+  w = chnkr.wstor;
+  wts = reshape(sqrt(sum((chnkr.d).^2,1)),k,nch);
+  wts = wts.*w(:);
 
 end
diff --git a/chunkie/chunkerkernevalmat.m b/chunkie/chunkerkernevalmat.m
index 4169d53..293e75b 100644
--- a/chunkie/chunkerkernevalmat.m
+++ b/chunkie/chunkerkernevalmat.m
@@ -273,13 +273,12 @@
     d = chnkr.d;
     n = chnkr.n;
     d2 = chnkr.d2;
-    h = chnkr.h;
     
     for i = 1:nch
         jmat = 1 + (i-1)*k*opdims(2);
         jmatend = i*k*opdims(2);
                         
-        mat(:,jmat:jmatend) =  chnk.adapgausswts(r,d,n,d2,h,ct,bw,i,targs, ...
+        mat(:,jmat:jmatend) =  chnk.adapgausswts(r,d,n,d2,ct,bw,i,targs, ...
                     targd,targd2,kern,opdims,t,w,opts);
                 
         js1 = jmat:jmatend;
diff --git a/chunkie/chunkermat.m b/chunkie/chunkermat.m
index f3fd7bd..e694208 100644
--- a/chunkie/chunkermat.m
+++ b/chunkie/chunkermat.m
@@ -573,14 +573,13 @@
 d = chnkr.d;
 n = chnkr.n;
 d2 = chnkr.d2;
-h = chnkr.h;
 
 for i = 1:nch
     jmat = 1 + (i-1)*k*opdims(2);
     jmatend = i*k*opdims(2);
                     
     [ji] = find(flag(:,i));
-    mat1 =  chnk.adapgausswts(r,d,n,d2,h,ct,bw,i,targs(:,ji), ...
+    mat1 =  chnk.adapgausswts(r,d,n,d2,ct,bw,i,targs(:,ji), ...
                 targd(:,ji),targn(:,ji),targd2(:,ji),kernev,opdims,t,w,opts);
             
     js1 = jmat:jmatend;
diff --git a/devtools/test/arclengthfunTest.m b/devtools/test/arclengthfunTest.m
new file mode 100644
index 0000000..1d11d70
--- /dev/null
+++ b/devtools/test/arclengthfunTest.m
@@ -0,0 +1,36 @@
+%ARCLENGTHFUNTEST tests the arclengthfun routine for computing distance
+% along length of curve
+% 
+
+clearvars; close all;
+seed = 8675309;
+rng(seed);
+addpaths_loc();
+
+% geometry parameters and construction
+
+
+cparams = [];
+cparams.eps = 1.0e-9;
+pref = []; 
+pref.k = 16;
+narms = 0;
+amp = 0.0;
+start = tic; chnkr = chunkerfunc(@(t) starfish(t,narms,amp),cparams,pref);
+
+[s, ~, chnkr] = arclengthfun(chnkr);
+% Compare to analytic arclength for circle
+ts = squeeze(atan2(chnkr.r(2,:,:), chnkr.r(1,:,:)));
+ts(ts<0) = ts(ts<0) + 2*pi;
+assert(norm(s-ts) < 1e-12);
+
+% Now test two circles
+chnkrs(2) = chunker();
+chnkrs(1) = chnkr;
+rfac = 1.1;
+chnkrs(2) = move(chnkr, [0;0], [3;0], 0, rfac);
+chnkrtotal = merge(chnkrs);
+[s, nchs, ~] = arclengthfun(chnkrtotal);
+
+assert(norm(s(:,1:nchs(1)) - ts) < 1e-12);
+assert(norm(s(:,nchs(1)+1:end) - rfac*ts) < 1e-12);
\ No newline at end of file

From 85111b60b6145f8920457f147d0e12b0c7b6bb9c Mon Sep 17 00:00:00 2001
From: Travis Askham <askhamwhat@gmail.com>
Date: Thu, 4 Apr 2024 14:50:58 -0400
Subject: [PATCH 2/7] add test of chunkerkerneval that triggers "all adaptive"
 option. clean some docs

---
 chunkie/+chnk/adapgausswts.m                    | 3 +--
 devtools/test/chunkerkernevalmat_greenlapTest.m | 9 +++++++--
 2 files changed, 8 insertions(+), 4 deletions(-)

diff --git a/chunkie/+chnk/adapgausswts.m b/chunkie/+chnk/adapgausswts.m
index 3fd69be..7afc58b 100644
--- a/chunkie/+chnk/adapgausswts.m
+++ b/chunkie/+chnk/adapgausswts.m
@@ -8,7 +8,7 @@
 % WARNING: this routine currently assumes that the kernel function is
 %   translation invariant to recenter (improves stability). 
 %
-% Syntax: [mat,maxrecs,numints,iers] = adapgausswts(r,d,d2,h,ct,bw,j, ...
+% Syntax: [mat,maxrecs,numints,iers] = adapgausswts(r,d,d2,ct,bw,j, ...
 %   rt,dt,d2t,kern,opdims,t,w,opts)
 %
 % Input:
@@ -16,7 +16,6 @@
 %   d - chnkr derivatives at nodes
 %   n - chnkr normals at nodes
 %   d2 - chnkr 2nd derivatives at nodes
-%   h - lengths of chunks in parameter space
 %   ct - Legendre nodes at order of chunker
 %   bw - barycentric interpolation weights for Legendre nodes at order of
 %   chunker
diff --git a/devtools/test/chunkerkernevalmat_greenlapTest.m b/devtools/test/chunkerkernevalmat_greenlapTest.m
index ed4913f..4baa848 100644
--- a/devtools/test/chunkerkernevalmat_greenlapTest.m
+++ b/devtools/test/chunkerkernevalmat_greenlapTest.m
@@ -6,7 +6,7 @@
 clearvars; close all;
 seed = 8675309;
 rng(seed);
-%addpaths_loc();
+addpaths_loc();
 
 doadap = false;
 
@@ -81,11 +81,16 @@
 start=tic; Dmat = chunkerkernevalmat(chnkr,kernd,targets,opts); 
 Du = Dmat*densu;
 toc(start)
+opts=[]; opts.forceadap=true;
+start=tic; Dmat = chunkerkernevalmat(chnkr,kernd,targets,opts); 
+Du2 = Dmat*densu;
+toc(start)
+opts=[];
 start=tic; Smat = chunkerkernevalmat(chnkr,kerns,targets,opts); 
 Sun = Smat*densun;
 toc(start)
 
-utarg2 = Sun-Du;
+utarg2 = Sun-Du2;
 
 %
 

From 31506ac05076d0706f0bd7c43c62efbaa5abcec9 Mon Sep 17 00:00:00 2001
From: Travis Askham <askhamwhat@gmail.com>
Date: Thu, 4 Apr 2024 14:54:33 -0400
Subject: [PATCH 3/7] fix bug in adaptive for everything

---
 chunkie/@chunker/arclengthfun.m                 |  2 +-
 chunkie/chunkerkernevalmat.m                    | 16 +---------------
 devtools/test/chunkerkernevalmat_greenlapTest.m |  6 +++++-
 3 files changed, 7 insertions(+), 17 deletions(-)

diff --git a/chunkie/@chunker/arclengthfun.m b/chunkie/@chunker/arclengthfun.m
index 63f7724..18dd033 100644
--- a/chunkie/@chunker/arclengthfun.m
+++ b/chunkie/@chunker/arclengthfun.m
@@ -14,7 +14,7 @@
 
 A = lege.intmat(chnkr.k);
 ds = arclengthdens(chnkr);
-chunklens = sum(chnkr.wts,1);
+chunklens = chunklen(chnkr);
 s = A*ds;
 
 for i = 1:ncomp
diff --git a/chunkie/chunkerkernevalmat.m b/chunkie/chunkerkernevalmat.m
index 293e75b..a6f5d1e 100644
--- a/chunkie/chunkerkernevalmat.m
+++ b/chunkie/chunkerkernevalmat.m
@@ -279,21 +279,7 @@
         jmatend = i*k*opdims(2);
                         
         mat(:,jmat:jmatend) =  chnk.adapgausswts(r,d,n,d2,ct,bw,i,targs, ...
-                    targd,targd2,kern,opdims,t,w,opts);
-                
-        js1 = jmat:jmatend;
-        js1 = repmat( (js1(:)).',1,opdims(1)*numel(ji));
-                
-        indji = (ji-1)*opdims(1);
-        indji = repmat( (indji(:)).', opdims(1),1) + ( (1:opdims(1)).');
-        indji = indji(:);
-        indji = repmat(indji,1,opdims(2)*k);
-        
-        iend = istart+numel(mat1)-1;
-        is(istart:iend) = indji(:);
-        js(istart:iend) = js1(:);
-        vs(istart:iend) = mat1(:);
-        istart = iend+1;
+                    targd,targn,targd2,kern,opdims,t,w,opts);
     end
     
 else
diff --git a/devtools/test/chunkerkernevalmat_greenlapTest.m b/devtools/test/chunkerkernevalmat_greenlapTest.m
index 4baa848..4ffa782 100644
--- a/devtools/test/chunkerkernevalmat_greenlapTest.m
+++ b/devtools/test/chunkerkernevalmat_greenlapTest.m
@@ -90,13 +90,17 @@
 Sun = Smat*densun;
 toc(start)
 
-utarg2 = Sun-Du2;
+utarg2 = Sun-Du;
+utarg3 = Sun-Du2;
 
 %
 
 relerr = norm(utarg-utarg2,'fro')/norm(utarg,'fro');
+assert(relerr < 1e-11);
 
 fprintf('relative frobenius error %5.2e\n',relerr);
 
+relerr = norm(utarg-utarg3,'fro')/norm(utarg,'fro');
 assert(relerr < 1e-11);
 
+

From 1531d87d83d52c55458f2c11020f736ca1555c6b Mon Sep 17 00:00:00 2001
From: Travis Askham <askhamwhat@gmail.com>
Date: Thu, 4 Apr 2024 14:55:37 -0400
Subject: [PATCH 4/7] pointing whts to weights

---
 chunkie/@chunker/whts.m | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/chunkie/@chunker/whts.m b/chunkie/@chunker/whts.m
index 3a7fc54..43bd30c 100644
--- a/chunkie/@chunker/whts.m
+++ b/chunkie/@chunker/whts.m
@@ -3,10 +3,6 @@
 
 warning('whts is deprecated and will be removed, use weights instead');
 
-  k = chnkr.k;
-  nch = chnkr.nch;
-  w = chnkr.wstor;
-  wts = reshape(sqrt(sum((chnkr.d).^2,1)),k,nch);
-  wts = wts.*w(:);
+    wts = weights(chnkr);
 
 end

From f4fe61a6578aa1eb3628eab0864f62ff2e417c1e Mon Sep 17 00:00:00 2001
From: Travis Askham <askhamwhat@gmail.com>
Date: Thu, 4 Apr 2024 15:09:36 -0400
Subject: [PATCH 5/7] fix .h bug in chunkerinterior for axissym, add test of
 chunkerinterior

---
 chunkie/chunkerinterior.m           |  1 -
 devtools/test/chunkerinteriorTest.m | 10 ++++++++--
 2 files changed, 8 insertions(+), 3 deletions(-)

diff --git a/chunkie/chunkerinterior.m b/chunkie/chunkerinterior.m
index 4a23291..f0c29fe 100644
--- a/chunkie/chunkerinterior.m
+++ b/chunkie/chunkerinterior.m
@@ -99,7 +99,6 @@
     chnkr.d2(:,:,istart:iend) = fliplr(chnkr.d2(:,:,1:nch));
     chnkr.d2(1,:,istart:iend) = chnkr.d2(1,:,istart:iend);
     
-    chnkr.h(istart:iend) = chnkr.h(1:nch);
     chnkr.wts = weights(chnkr);
     chnkr.n = normals(chnkr);
     chnkr.adj(1,:) = 0:chnkr.nch-1;
diff --git a/devtools/test/chunkerinteriorTest.m b/devtools/test/chunkerinteriorTest.m
index 2538be6..281d205 100644
--- a/devtools/test/chunkerinteriorTest.m
+++ b/devtools/test/chunkerinteriorTest.m
@@ -82,5 +82,11 @@
 fprintf('%5.2e s : time for chunkerinterior (chunker, with fmm)\n',t5);
 assert(all(in5(:) == 1));
 
-
-
+% test axissym option 
+chnkr = chunkerfunc(@(t) starfish(t),struct('ta',-pi/2,'tb',pi/2,'ifclosed',0));
+nt = 1000;
+ttarg = -pi/2+pi*rand(nt,1); scal = 2*rand(1,nt);
+targs = starfish(ttarg).*scal;
+
+in = chunkerinterior(chnkr,targs,struct('axissym',true));
+assert(all(in(:) == (scal(:) <= 1)));
\ No newline at end of file

From 2ddf268ff786029736ad689bd1a9d4e9f9b7bcd0 Mon Sep 17 00:00:00 2001
From: Travis Askham <askhamwhat@gmail.com>
Date: Fri, 5 Apr 2024 10:05:59 -0400
Subject: [PATCH 6/7] fix .h issue in nearest/nearparam, fixes other bugs, and
 removes duplicate code.

---
 chunkie/+chnk/chunk_nearparam.m |  33 +++++----
 chunkie/@chunker/nearest.m      | 119 ++------------------------------
 2 files changed, 24 insertions(+), 128 deletions(-)

diff --git a/chunkie/+chnk/chunk_nearparam.m b/chunkie/+chnk/chunk_nearparam.m
index efdf2fd..9d670cd 100644
--- a/chunkie/+chnk/chunk_nearparam.m
+++ b/chunkie/+chnk/chunk_nearparam.m
@@ -24,6 +24,7 @@
 
 maxnewt = 15;
 thresh0 = 1.0d-14;
+iextra = 3;
 
 if nargin < 4 || isempty(t) || nargin < 5 || isempty(u)
     [t,~,u] = lege.exps(chnkr.k);
@@ -40,12 +41,6 @@
     thresh0 = opts.thresh;
 end
 
-dist = Inf;
-
-if nargin < 3
-    ich = 1:nch;
-end
-
 dim = size(pts,1);
 npts = numel(pts)/dim;
 pts = reshape(pts,dim,npts);
@@ -92,9 +87,10 @@
 
     rdiff0 = r0(:) - ref;
     dprime = rdiff0.'*d0(:);
-    dprime2 = d0(:).'*d0(:) + d0(:).'*d20(:);
+    dprime2 = (d0(:).'*d0(:) + rdiff0(:).'*d20(:));
 
     newtsuccess = false;
+    ifend = 0;
 
     for iii = 1:maxnewt
 
@@ -102,14 +98,17 @@
         
         deltat = - dprime/dprime2;
 
-        t0 = t0+deltat;
+        t1 = t0 + deltat;
         
-        if (t0 > 1.0)
-            t0 = 1;
+        if (t1 > 1.0)
+            t1 = 1;
         end
-        if (t0 < -1.0)
-            t0 = -1;
+        if (t1 < -1.0)
+            t1 = -1;
         end
+    
+        deltat = t1-t0;
+        t0 = t1;
 
         all0 = (lege.exev(t0,cfs));
         r0 = all0(1:dim);
@@ -118,9 +117,13 @@
 
         rdiff0 = r0(:)-ref(:);
         dprime = rdiff0.'*d0(:);
-        dprime2 = d0(:).'*d0(:) + d0(:).'*d20(:);
+        dprime2 = d0(:).'*d0(:) + rdiff0(:).'*d20(:);
+
+        if min(abs(dprime),abs(deltat)) < thresh
+            ifend = ifend+1;
+        end
 
-        if abs(dprime) < thresh
+        if (ifend >= iextra)
             newtsuccess = true;
             break;
         end
@@ -196,7 +199,7 @@
                 d20 = d21;
                 rdiff0 = rdiff1;
                 dprime = rdiff0.'*d0(:);
-                dprime2 = d0(:).'*d0(:) + d0(:).'*d20(:);
+                dprime2 = d0(:).'*d0(:);
                 dist0 = dist1;
                 lam = lam/3;
                 if abs(dprime) < thresh
diff --git a/chunkie/@chunker/nearest.m b/chunkie/@chunker/nearest.m
index 7f1bbdb..f87c7d3 100644
--- a/chunkie/@chunker/nearest.m
+++ b/chunkie/@chunker/nearest.m
@@ -40,6 +40,7 @@
 maxnewt = 15;
 thresh0 = 1.0d-9;
 iextra = 3;
+nch = chnkr.nch;
 
 if nargin < 5 || isempty(u)
     [~,~,u] = lege.exps(chnkr.k);
@@ -69,123 +70,15 @@
     % grab chunk data
     ii = ich(i);
     r = chnkr.rstor(:,:,ii);
-    d = chnkr.dstor(:,:,ii);
-    d2 = chnkr.d2stor(:,:,ii);
-    rc = u*(r.');
-    dc = u*(d.');
-    d2c = u*(d2.');
 
-    % starting guess on over sampled grid
-    rover = (ainterpover*(r.')).';
-    rdist0 = sqrt(sum((rover-ref(:)).^2,1));
-    [disti,ind] = min(rdist0);
-    tt = xover(ind);
-    
-    % check the endpoints
-    ts = [-1;1];
-    rs = (lege.exev(ts,rc)).';
-    rsdist = sqrt(sum((rs-ref).^2,1));
-    [dists,inds] = min(rsdist);
-    if dists < disti
-        % continue if endpoints closer
-        ri = rs(:,inds);
-        if dists < dist
-            dist = dists;
-            ds = (lege.exev(ts,dc)).';
-            d2s = (lege.exev(ts,d2c)).';
-            di = ds(:,inds);
-            d2i = d2s(:,inds);
-            rn = ri; dn = di; d2n = d2i;
-            tn = ts(inds);
-            ichn = ii;
-        end
-        continue;
-    end
-        
-    % run levenberg
-
-    t0 = tt;
-
-    ifend = 0;
-
-    r0 = (lege.exev(t0,rc)).';
-
-    rdiff0 = ref(:) - r0(:);
-
-    drc = lege.derpol(rc);
-    %d2rc = lege.derpol(drc);
-    thresh = thresh0;
-    
-    for iii = 1:maxnewt
-        %d0 = (lege.exev(t0,dc)).'*h;
-        d0 = (lege.exev(t0,drc)).'; % actual derivative of chunk
-%        d20 = (lege.exev(t0,d2rc)).'; % actual 2nd derivative of chunk
-
-        % newton info
-%         
-%         dprime = rdiff0.'*d0(:);
-%         dprime2 = (d0(:).'*d0(:) + d0(:).'*d20(:));
-
-        % levenberg instead
-        
-        deltat = d0(:)\(rdiff0(:));
-        t1 = t0+deltat;
-        
-        if (t1 > 1.0)
-            t1 = (1.0+t0)/2;
-        end
-        if (t1 < -1.0)
-            t1 = (-1.0+t0)/2;
-        end
-
-        r1 = (lege.exev(t1,rc)).';
-        rdiff1 = ref(:)-r1(:);
-
-        err = min(abs(deltat),abs(sum(d0(:).*rdiff0(:))));
-        if err < thresh
-            ifend = ifend+1;
-        end
-        
-        %sum(rdiff0(:).*d0(:))
-        %abs(deltat)
-        
-        t0 = t1;
-        rdiff0 = rdiff1;
-
-%         if (err < thresh)
-%             ifend = ifend + 1;
-%         end
-% 
-        %t0 = t1;
-
-        if (ifend >= iextra)
-            break;
-        end
-    end
-    
-    % done levenberg
-    
-    ri = (lege.exev(t0,rc)).';
-    disti = sqrt(sum((ri(:)-ref(:)).^2,1));
-    if (ifend < iextra)
-        warning('newton failed in nearest');
-%          deltat
-%          disti
-%          jac
-%          t0
-%          rhs
-%          ifend
-%          sum(jac.*rhs)
-    else
-%        fprintf('success\n')
-    end
-    
+    [ti,ri,di,d2i,disti] = chnk.chunk_nearparam(r,ref,opts,chnkr.tstor,u);
+    disti = sqrt(disti);
     if disti < dist
         dist = disti;
         rn = ri;
-        dn = (lege.exev(t1,dc)).'*h;
-        d2n = (lege.exev(t1,d2c)).'*h*h;
-        tn = t1;
+        dn = di;
+        d2n = d2i;
+        tn = ti;
         ichn = ii;
     end
 end

From 6668915723533d114b8833fff0d5d324313d8dbb Mon Sep 17 00:00:00 2001
From: Travis Askham <askhamwhat@gmail.com>
Date: Fri, 5 Apr 2024 10:08:01 -0400
Subject: [PATCH 7/7] adds a test for nearest utilities

---
 devtools/test/chunker_nearestTest.m | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)
 create mode 100644 devtools/test/chunker_nearestTest.m

diff --git a/devtools/test/chunker_nearestTest.m b/devtools/test/chunker_nearestTest.m
new file mode 100644
index 0000000..54ab6c5
--- /dev/null
+++ b/devtools/test/chunker_nearestTest.m
@@ -0,0 +1,24 @@
+%chunker.nearest test
+%
+
+iseed = 1234;
+rng(iseed);
+
+% compare closest computed point on circle geometry to exact answer
+
+circ = chunkerfunc(@(t) chnk.curves.bymode(t,1));
+nt = 1000;
+thetas = -pi+2*pi*rand(1,nt);  
+
+% the 0.1 here keeps targets away from coordinate singularity of circle
+scal = 0.1+2*rand(1,nt);
+targs = [cos(thetas); sin(thetas)].*scal;
+
+for j = 1:nt
+    [rn,dn,dist,tn,ichn] = nearest(circ,targs(:,j));
+    th1 = atan2(targs(2,j),targs(1,j));
+    th2 = atan2(rn(2),rn(1));
+    err = min(abs(th1-th2),abs(th1-th2+2*pi));
+    err = min(err,abs(th1-th2-2*pi));
+    assert(err < 1e-12);
+end
\ No newline at end of file