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;