wcs.c

            /* Refine the solution interval. */
            while (lng0 > span[0]) {
               lng0 -= step;
               if (lng0 < span[0]) lng0 = span[0];
               world[wcs->lng] = lng0;
               if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
                          prj, imgcrd, lin, pixcrd))) {
                  return err;
               }
               d0 = pixcrd[mixpix] - pixmix;

               /* Check for a solution. */
               dabs = fabs(d0);
               if (dabs < tol) return 0;

               /* Record the point of closest approach. */
               if (dabs < dmin) {
                  lmin = lng0;
                  dmin = dabs;
               }

               /* Check for a crossing point. */
               if (wcs_signbit(d0) != wcs_signbit(d1)) {
                  crossed = 2;
                  dx = d0;
                  break;
               }

               /* Advance to the next subinterval. */
               lng1 = lng0;
               d1 = d0;
            }

            if (crossed) {
               /* A crossing point was found. */
               for (iter = 0; iter < niter; iter++) {
                  /* Use regula falsi division of the interval. */
                  lambda = d0/(d0-d1);
                  if (lambda < 0.1) {
                     lambda = 0.1;
                  } else if (lambda > 0.9) {
                     lambda = 0.9;
                  }

                  lng = lng0 + lambda*(lng1 - lng0);
                  world[wcs->lng] = lng;
                  if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
                                   prj, imgcrd, lin, pixcrd))) {
                     return err;
                  }
                  d = pixcrd[mixpix] - pixmix;

                  /* Check for a solution. */
                  dabs = fabs(d);
                  if (dabs < tol) return 0;

                  /* Record the point of closest approach. */
                  if (dabs < dmin) {
                     lmin = lng;
                     dmin = dabs;
                  }

                  if (wcs_signbit(d0) == wcs_signbit(d)) {
                     lng0 = lng;
                     d0 = d;
                  } else {
                     lng1 = lng;
                     d1 = d;
                  }
               }

               /* No convergence, must have been a discontinuity. */
               if (crossed == 1) lng0 = span[1];
               lng1 = lng0;
               d1 = dx;
               crossed = 0;

            } else {
               /* No crossing point; look for a tangent point. */
               if (lmin == span[0]) break;
               if (lmin == span[1]) break;

               lng = lmin;
               lng0 = lng - step;
               if (lng0 < span[0]) lng0 = span[0];
               lng1 = lng + step;
               if (lng1 > span[1]) lng1 = span[1];

               world[wcs->lng] = lng0;
               if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
                                prj, imgcrd, lin, pixcrd))) {
                  return err;
               }
               d0 = fabs(pixcrd[mixpix] - pixmix);

               d  = dmin;

               world[wcs->lng] = lng1;
               if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
                                prj, imgcrd, lin, pixcrd))) {
                  return err;
               }
               d1 = fabs(pixcrd[mixpix] - pixmix);

               for (iter = 0; iter < niter; iter++) {
                  lng0m = (lng0 + lng)/2.0;
                  world[wcs->lng] = lng0m;
                  if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
                                   prj, imgcrd, lin, pixcrd))) {
                     return err;
                  }
                  d0m = fabs(pixcrd[mixpix] - pixmix);

                  if (d0m < tol) return 0;

                  lng1m = (lng1 + lng)/2.0;
                  world[wcs->lng] = lng1m;
                  if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
                                   prj, imgcrd, lin, pixcrd))) {
                     return err;
                  }
                  d1m = fabs(pixcrd[mixpix] - pixmix);

                  if (d1m < tol) return 0;

                  if (d0m < d && d0m <= d1m) {
                     lng1 = lng;
                     d1   = d;
                     lng  = lng0m;
                     d    = d0m;
                  } else if (d1m < d) {
                     lng0 = lng;
                     d0   = d;
                     lng  = lng1m;
                     d    = d1m;
                  } else {
                     lng0 = lng0m;
                     d0   = d0m;
                     lng1 = lng1m;
                     d1   = d1m;
                  }
               }
            }
         }
      }
   }


   /* Set cel0 to the unity transformation. */
   cel0.flag = CELSET;
   cel0.ref[0] = cel->ref[0];
   cel0.ref[1] = cel->ref[1];
   cel0.ref[2] = cel->ref[2];
   cel0.ref[3] = cel->ref[3];
   cel0.euler[0] = -90.0;
   cel0.euler[1] =   0.0;
   cel0.euler[2] =  90.0;
   cel0.euler[3] =   1.0;
   cel0.euler[4] =   0.0;
   cel0.prjfwd = cel->prjfwd;
   cel0.prjrev = cel->prjrev;

   /* No convergence, check for aberrant behaviour at a native pole. */
   *theta = -90.0;
   for (j = 1; j <= 2; j++) {
      /* Could the celestial coordinate element map to a native pole? */
      *theta = -*theta;
      err = sphrev(0.0, *theta, cel->euler, &lng, &lat);

      if (mixcel == 1) {
         if (fabs(fmod(world[wcs->lng]-lng,360.0)) > tol) continue;
         if (lat < span[0]) continue;
         if (lat > span[1]) continue;
         world[wcs->lat] = lat;
      } else {
         if (fabs(world[wcs->lat]-lat) > tol) continue;
         if (lng < span[0]) lng += 360.0;
         if (lng > span[1]) lng -= 360.0;
         if (lng < span[0]) continue;
         if (lng > span[1]) continue;
         world[wcs->lng] = lng;
      }

      /* Is there a solution for the given pixel coordinate element? */
      lng = world[wcs->lng];
      lat = world[wcs->lat];

      /* Feed native coordinates to wcsfwd() with cel0 set to unity. */
      world[wcs->lng] = -180.0;
      world[wcs->lat] = *theta;
      if ((err = wcsfwd(ctype, wcs, world, crval, &cel0, phi, theta, prj,
                       imgcrd, lin, pixcrd))) {
         return err;
      }
      d0 = pixcrd[mixpix] - pixmix;

      /* Check for a solution. */
      if (fabs(d0) < tol) {
         /* Recall saved world coordinates. */
         world[wcs->lng] = lng;
         world[wcs->lat] = lat;
         return 0;
      }

      /* Search for a crossing interval. */
      phi0 = -180.0;
      for (k = -179; k <= 180; k++) {
         phi1 = (float) k;
         world[wcs->lng] = phi1;
         if ((err = wcsfwd(ctype, wcs, world, crval, &cel0, phi, theta, prj,
                          imgcrd, lin, pixcrd))) {
            return err;
         }
         d1 = pixcrd[mixpix] - pixmix;

         /* Check for a solution. */
         dabs = fabs(d1);
         if (dabs < tol) {
            /* Recall saved world coordinates. */
            world[wcs->lng] = lng;
            world[wcs->lat] = lat;
            return 0;
         }

         /* Is it a crossing interval? */
         if (wcs_signbit(d0) != wcs_signbit(d1)) break;

         phi0 = phi1;
         d0 = d1;
      }

      for (iter = 1; iter <= niter; iter++) {
         /* Use regula falsi division of the interval. */
         lambda = d0/(d0-d1);
         if (lambda < 0.1) {
            lambda = 0.1;
         } else if (lambda > 0.9) {
            lambda = 0.9;
         }
 
         world[wcs->lng] = phi0 + lambda*(phi1 - phi0);
         if ((err = wcsfwd(ctype, wcs, world, crval, &cel0, phi, theta, prj,
                          imgcrd, lin, pixcrd))) {
            return err;
         }
         d = pixcrd[mixpix] - pixmix;
 
         /* Check for a solution. */
         dabs = fabs(d);
         if (dabs < tol) {
            /* Recall saved world coordinates. */
            world[wcs->lng] = lng;
            world[wcs->lat] = lat;
            return 0;
         }
 
         if (wcs_signbit(d0) == wcs_signbit(d)) {
            phi0 = world[wcs->lng];
            d0 = d;
         } else {
            phi1 = world[wcs->lng];
            d1 = d;
         }
      }
   }


   /* No solution. */
   return 5;

}