Horizontal grid shift: fix issue on iterative inverse computation when switching between (sub)grids (fixes #1663)

Given in.txt with
53.999759140 5.144478208 252.6995

Before the fix,
cct  -t 0 -d 4 +proj=pipeline +step +proj=axisswap +order=2,1,3,4 +step +proj=hgridshift +inv +grids=rdtrans2018.gsb +step +proj=vgridshift +grids=naptrans2018.gtx +step +proj=sterea +lat_0=52.156160556 +lon_0=5.387638889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel in.txt

returned:
  139079.8814    668306.0302      212.1724        0.0000

It now returns:
  139079.8850    668306.0458      212.1724        0.0000

which meets with the 1mm accuracy the expected result of test point
```
30010049    53.999759140    5.144478208     252.6995        139079.8850     668306.0460     212.1723
```

1 files changed, 32 insertions, 11 deletions

diff --git a/src/nad_cvt.cpp b/src/nad_cvt.cpp
index 79441d0a..e8b8e9b7 100644
--- a/src/nad_cvt.cpp
+++ b/src/nad_cvt.cpp
@@ -7,10 +7,12 @@
 #include "proj_internal.h"
 #include <math.h>
 
+#include <limits>
+
 #define MAX_ITERATIONS 10
 #define TOL 1e-12
 
-PJ_LP nad_cvt(PJ_LP in, int inverse, struct CTABLE *ct) {
+PJ_LP nad_cvt(projCtx ctx, PJ_LP in, int inverse, struct CTABLE *ct, int grid_count, PJ_GRIDINFO **tables) {
     PJ_LP t, tb,del, dif;
     int i = MAX_ITERATIONS;
     const double toltol = TOL*TOL;
@@ -40,18 +42,37 @@ PJ_LP nad_cvt(PJ_LP in, int inverse, struct CTABLE *ct) {
     do {
         del = nad_intr(t, ct);
 
-        /* This case used to return failure, but I have
-           changed it to return the first order approximation
-           of the inverse shift.  This avoids cases where the
-           grid shift *into* this grid came from another grid.
-           While we aren't returning optimally correct results
-           I feel a close result in this case is better than
-           no result.  NFW
-           To demonstrate use -112.5839956 49.4914451 against
-           the NTv2 grid shift file from Canada. */
-        if (del.lam == HUGE_VAL)
+        /* In case we are (or have switched) on the null grid, exit now */
+        if( del.lam == 0 && del.phi == 0 )
             break;
 
+        /* We can possibly go outside of the initial guessed grid, so try */
+        /* to fetch a new grid into which iterate... */
+        if (del.lam == HUGE_VAL)
+        {
+            if( grid_count == 0 )
+                break;
+            PJ_LP lp;
+            lp.lam = t.lam + ct->ll.lam;
+            lp.phi = t.phi + ct->ll.phi;
+            auto newCt = find_ctable(ctx, lp, grid_count, tables);
+            if( newCt == nullptr || newCt == ct )
+                break;
+            pj_log(ctx, PJ_LOG_DEBUG_MINOR, "Switching from grid %s to grid %s",
+                   ct->id,
+                   newCt->id);
+            ct = newCt;
+            t.lam = lp.lam - ct->ll.lam;
+            t.phi = lp.phi - ct->ll.phi;
+            tb = in;
+            tb.lam -= ct->ll.lam;
+            tb.phi -= ct->ll.phi;
+            tb.lam = adjlon (tb.lam - M_PI) + M_PI;
+            dif.lam = std::numeric_limits<double>::max();
+            dif.phi = std::numeric_limits<double>::max();
+            continue;
+        }
+
         dif.lam = t.lam - del.lam - tb.lam;
            dif.phi = t.phi + del.phi - tb.phi;
         t.lam -= dif.lam;