struct CTABLE *nad_ctable_init( projCtx ctx, FILE * fid )
{
    struct CTABLE *ct;
    int		id_end;

    /* read the table header */
    ct = (struct CTABLE *) pj_malloc(sizeof(struct CTABLE));
    if( ct == NULL 
        || fread( ct, sizeof(struct CTABLE), 1, fid ) != 1 )
    {
        pj_ctx_set_errno( ctx, -38 );
        return NULL;
    }

    /* do some minimal validation to ensure the structure isn't corrupt */
    if( ct->lim.lam < 1 || ct->lim.lam > 100000 
        || ct->lim.phi < 1 || ct->lim.phi > 100000 )
    {
        pj_ctx_set_errno( ctx, -38 );
        return NULL;
    }
    
    /* trim white space and newlines off id */
    for( id_end = strlen(ct->id)-1; id_end > 0; id_end-- )
    {
        if( ct->id[id_end] == '\n' || ct->id[id_end] == ' ' )
            ct->id[id_end] = '\0';
        else
            break;
    }

    ct->cvs = NULL;

    return ct;
}

int nad_ctable2_load( projCtx ctx, struct CTABLE *ct, FILE *fid )

{
    int  a_size;

    fseek( fid, 160, SEEK_SET );

    /* read all the actual shift values */
    a_size = ct->lim.lam * ct->lim.phi;
    ct->cvs = (FLP *) pj_malloc(sizeof(FLP) * a_size);
    if( ct->cvs == NULL 
        || fread(ct->cvs, sizeof(FLP), a_size, fid) != a_size )
    {
        pj_dalloc( ct->cvs );
        ct->cvs = NULL;

        if( getenv("PROJ_DEBUG") != NULL )
        {
            fprintf( stderr,
            "ctable2 loading failed on fread() - binary incompatible?\n" );
        }

        pj_ctx_set_errno( ctx, -38 );
        return 0;
    }

    if( !IS_LSB )
    {
        swap_words( ct->cvs, 4, a_size * 2 );
    }

    return 1;
} 

	XY /* forward projection entry */
pj_fwd(LP lp, PJ *P) {
	XY xy;
	double t;

	/* check for forward and latitude or longitude overange */
	if ((t = fabs(lp.phi)-HALFPI) > EPS || fabs(lp.lam) > 10.) {
		xy.x = xy.y = HUGE_VAL;
		pj_ctx_set_errno( P->ctx, -14);
	} else { /* proceed with projection */
                P->ctx->last_errno = 0;
                pj_errno = 0;
                errno = 0;

		if (fabs(t) <= EPS)
			lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
		else if (P->geoc)
			lp.phi = atan(P->rone_es * tan(lp.phi));
		lp.lam -= P->lam0;	/* compute del lp.lam */
		if (!P->over)
			lp.lam = adjlon(lp.lam); /* adjust del longitude */
		xy = (*P->fwd)(lp, P); /* project */
		if ( P->ctx->last_errno )
			xy.x = xy.y = HUGE_VAL;
		/* adjust for major axis and easting/northings */
		else {
			xy.x = P->fr_meter * (P->a * xy.x + P->x0);
			xy.y = P->fr_meter * (P->a * xy.y + P->y0);
		}
	}
	return xy;
}

static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
    XY xy = {0.0,0.0};
    double b, cosphi;

    /*
     * Fail if our longitude is more than 90 degrees from the
     * central meridian since the results are essentially garbage.
     * Is error -20 really an appropriate return value?
     *
     *  http://trac.osgeo.org/proj/ticket/5
     */
    if( lp.lam < -HALFPI || lp.lam > HALFPI ) {
        xy.x = HUGE_VAL;
        xy.y = HUGE_VAL;
        pj_ctx_set_errno( P->ctx, -14 );
        return xy;
    }

    cosphi = cos(lp.phi);
    b = cosphi * sin (lp.lam);
    if (fabs (fabs (b) - 1.) <= EPS10)
        F_ERROR;

    xy.x = P->opaque->ml0 * log ((1. + b) / (1. - b));
    xy.y = cosphi * cos (lp.lam) / sqrt (1. - b * b);

    b = fabs ( xy.y );
    if (b >= 1.) {
        if ((b - 1.) > EPS10)
            F_ERROR
        else xy.y = 0.;
    } else

double aacos(projCtx ctx, double v) {
	double av;

	if ((av = fabs(v)) >= 1.) {
		if (av > ONE_TOL)
			pj_ctx_set_errno(ctx, -19);
		return (v < 0. ? PI : 0.);
	}
	return acos(v);
}

double aasin(projCtx ctx, double v) {
	double av;

	if ((av = fabs(v)) >= 1.) {
		if (av > ONE_TOL)
			pj_ctx_set_errno(ctx, -19);
		return (v < 0. ? -HALFPI : HALFPI);
	}
	return asin(v);
}

static PJ_LP e_inverse (PJ_XY xy, PJ *P) {           /* Ellipsoidal/spheroidal, inverse */
    PJ_LP lp = {0.0,0.0};
    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
    double yc, y2, y6;
    int i;

    /* Adjusting x and y for authalic radius */
    xy.x /= Q->rqda;
    xy.y /= Q->rqda;

    /* Make sure y is inside valid range */
    if (xy.y > MAX_Y)
        xy.y = MAX_Y;
    else if (xy.y < -MAX_Y)
        xy.y = -MAX_Y;

    yc = xy.y;

    /* Newton-Raphson */
    for (i = MAX_ITER; i ; --i) {
        double f, fder, tol;

        y2 = yc * yc;
        y6 = y2 * y2 * y2;

        f    = yc * (A1 + A2 * y2 + y6 * (A3 + A4 * y2)) - xy.y;
        fder = A1 + 3 * A2 * y2 + y6 * (7 * A3 + 9 * A4 * y2);

        tol  = f / fder;
        yc  -= tol;

        if (fabs(tol) < EPS)
            break;
    }

    if( i == 0 ) {
        pj_ctx_set_errno( P->ctx, PJD_ERR_NON_CONVERGENT );
        return lp;
    }

    /* Longitude */
    y2 = yc * yc;
    y6 = y2 * y2 * y2;

    lp.lam = M * xy.x * (A1 + 3 * A2 * y2 + y6 * (7 * A3 + 9 * A4 * y2)) / cos(yc);

    /* Latitude (for spheroidal case, this is latitude */
    lp.phi = asin(sin(yc) / M);

    /* Ellipsoidal case, converting auth. latitude */
    if (P->es != 0.0)
       lp.phi = pj_authlat(lp.phi, Q->apa);

    return lp;
}

static LP s_healpix_inverse(XY xy, PJ *P) { /* sphere */
    LP lp = {0.0,0.0};

    /* Check whether (x, y) lies in the HEALPix image */
    if (in_image(xy.x, xy.y, 0, 0, 0) == 0) {
        lp.lam = HUGE_VAL;
        lp.phi = HUGE_VAL;
        pj_ctx_set_errno(P->ctx, -15);
        return lp;
    }
    return healpix_sphere_inverse(xy);
}

static LP e_healpix_inverse(XY xy, PJ *P) { /* ellipsoid */
    LP lp = {0.0,0.0};

    /* Check whether (x, y) lies in the HEALPix image. */
    if (in_image(xy.x, xy.y, 0, 0, 0) == 0) {
        lp.lam = HUGE_VAL;
        lp.phi = HUGE_VAL;
        pj_ctx_set_errno(P->ctx, -15);
        return lp;
    }
    lp = healpix_sphere_inverse(xy);
    lp.phi = auth_lat(P, lp.phi, 1);
    return lp;
}

static LP s_rhealpix_inverse(XY xy, PJ *P) { /* sphere */
    struct pj_opaque *Q = P->opaque;
    LP lp = {0.0,0.0};

    /* Check whether (x, y) lies in the rHEALPix image. */
    if (in_image(xy.x, xy.y, 1, Q->north_square, Q->south_square) == 0) {
        lp.lam = HUGE_VAL;
        lp.phi = HUGE_VAL;
        pj_ctx_set_errno(P->ctx, -15);
        return lp;
    }
    xy = combine_caps(xy.x, xy.y, Q->north_square, Q->south_square, 1);
    return healpix_sphere_inverse(xy);
}

static PJ_LP krovak_e_inverse (PJ_XY xy, PJ *P) {                /* Ellipsoidal, inverse */
    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
    PJ_LP lp = {0.0,0.0};

    double u, deltav, s, d, eps, rho, fi1, xy0;
    int i;

    xy0 = xy.x;
    xy.x = xy.y;
    xy.y = xy0;

    xy.x *= Q->czech;
    xy.y *= Q->czech;

    rho = sqrt(xy.x * xy.x + xy.y * xy.y);
    eps = atan2(xy.y, xy.x);

    d = eps / sin(S0);
    if( rho == 0.0 ) {
        s = M_PI_2;
    }
    else {
        s = 2. * (atan(  pow(Q->rho0 / rho, 1. / Q->n) * tan(S0 / 2. + M_PI_4)) - M_PI_4);
    }

    u = asin(cos(Q->ad) * sin(s) - sin(Q->ad) * cos(s) * cos(d));
    deltav = asin(cos(s) * sin(d) / cos(u));

    lp.lam = P->lam0 - deltav / Q->alpha;

    /* ITERATION FOR lp.phi */
    fi1 = u;

    for (i = MAX_ITER; i ; --i) {
        lp.phi = 2. * ( atan( pow( Q->k, -1. / Q->alpha)  *
                              pow( tan(u / 2. + M_PI_4) , 1. / Q->alpha)  *
                              pow( (1. + P->e * sin(fi1)) / (1. - P->e * sin(fi1)) , P->e / 2.)
                            )  - M_PI_4);

        if (fabs(fi1 - lp.phi) < EPS)
            break;
        fi1 = lp.phi;
    }
    if( i == 0 )
        pj_ctx_set_errno( P->ctx, PJD_ERR_NON_CONVERGENT );

   lp.lam -= P->lam0;

   return lp;
}

static PJ_LP s_inverse (PJ_XY xy, PJ *P) {           /* Spheroidal, inverse */
    PJ_LP lp = {0.0,0.0};
    long i;
    double t, t1;
    struct COEFS T;
    int iters;

    lp.lam = xy.x / FXC;
    lp.phi = fabs(xy.y / FYC);
    if (lp.phi >= 1.) { /* simple pathologic cases */
        if (lp.phi > ONEEPS) {
            proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
            return lp;
        }
        else {
            lp.phi = xy.y < 0. ? -M_HALFPI : M_HALFPI;
            lp.lam /= X[NODES].c0;
        }
    } else { /* general problem */
        /* in Y space, reduce to table interval */
        i = isnan(lp.phi) ? -1 : lround(floor(lp.phi * NODES));
        if( i < 0 || i >= NODES ) {
            proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
            return lp;
        }
        for (;;) {
            if (Y[i].c0 > lp.phi) --i;
            else if (Y[i+1].c0 <= lp.phi) ++i;
            else break;
        }
        T = Y[i];
        /* first guess, linear interp */
        t = 5. * (lp.phi - T.c0)/(Y[i+1].c0 - T.c0);
        /* make into root */
        T.c0 = (float)(T.c0 - lp.phi);
        for (iters = MAX_ITER; iters ; --iters) { /* Newton-Raphson */
            t -= t1 = V(T,t) / DV(T,t);
            if (fabs(t1) < EPS)
                break;
        }
        if( iters == 0 )
            pj_ctx_set_errno( P->ctx, PJD_ERR_NON_CONVERGENT );
        lp.phi = (5 * i + t) * DEG_TO_RAD;
        if (xy.y < 0.) lp.phi = -lp.phi;
        lp.lam /= V(X[i], t);
    }
    return lp;
}

static paralist *
get_init(projCtx ctx, paralist **start, paralist *next, char *name) {
    char fname[MAX_PATH_FILENAME+ID_TAG_MAX+3], *opt;
    FILE *fid;
    paralist *init_items = NULL;
    const paralist *orig_next = next;

    (void)strncpy(fname, name, MAX_PATH_FILENAME + ID_TAG_MAX + 1);
	
    /* 
    ** Search for file/key pair in cache 
    */
	
    init_items = pj_search_initcache( name );
    if( init_items != NULL )
    {
        next->next = init_items;
        while( next->next != NULL )
            next = next->next;
        return next;
    }

    /*
    ** Otherwise we try to open the file and search for it.
    */
    if ((opt = strrchr(fname, ':')) != NULL)
        *opt++ = '\0';
    else { pj_ctx_set_errno(ctx,-3); return NULL; }

    if ( (fid = pj_open_lib(ctx,fname, "rt")) != NULL)
        next = get_opt(ctx, start, fid, opt, next);
    else
        return NULL;
    (void)fclose(fid);
    if (errno == 25)
        errno = 0; /* unknown problem with some sys errno<-25 */

    /* 
    ** If we seem to have gotten a result, insert it into the 
    ** init file cache.
    */
    if( next != NULL && next != orig_next )
        pj_insert_initcache( name, orig_next->next );

    return next;
}

	double
pj_phi2(projCtx ctx, double ts, double e) {
	double eccnth, Phi, con, dphi;
	int i;

	eccnth = .5 * e;
	Phi = HALFPI - 2. * atan (ts);
	i = N_ITER;
	do {
		con = e * sin (Phi);
		dphi = HALFPI - 2. * atan (ts * pow((1. - con) /
		   (1. + con), eccnth)) - Phi;
		Phi += dphi;
	} while ( fabs(dphi) > TOL && --i);
	if (i <= 0)
		pj_ctx_set_errno( ctx, -18 );
	return Phi;
}

double proj_inv_mdist(projCtx ctx, double dist, const void *b) {
	double s, t, phi, k;
	int i;

	k = 1. / (1. - B->es);
	i = MAX_ITER;
	phi = dist;
	while (i--) {
		s = sin(phi);
		t = 1. - B->es * s * s;
		phi -= t = (proj_mdist(phi, s, cos(phi), b) - dist) * (t * sqrt(t)) * k;
		if (fabs(t) < TOL) /* that is no change */
			return phi;
	}
	/* convergence failed */
	pj_ctx_set_errno(ctx, -17);
	return phi;
}

static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
    XY xy = {0.0, 0.0};
    struct pj_opaque *Q = P->opaque;
    double al, als, n, cosphi, sinphi, t;

    /*
     * Fail if our longitude is more than 90 degrees from the
     * central meridian since the results are essentially garbage.
     * Is error -20 really an appropriate return value?
     *
     *  http://trac.osgeo.org/proj/ticket/5
     */
    if( lp.lam < -HALFPI || lp.lam > HALFPI ) {
        xy.x = HUGE_VAL;
        xy.y = HUGE_VAL;
        pj_ctx_set_errno( P->ctx, -14 );
        return xy;
    }

    sinphi = sin (lp.phi);
    cosphi = cos (lp.phi);
    t = fabs (cosphi) > 1e-10 ? sinphi/cosphi : 0.;
    t *= t;
    al = cosphi * lp.lam;
    als = al * al;
    al /= sqrt (1. - P->es * sinphi * sinphi);
    n = Q->esp * cosphi * cosphi;
    xy.x = P->k0 * al * (FC1 +
        FC3 * als * (1. - t + n +
        FC5 * als * (5. + t * (t - 18.) + n * (14. - 58. * t)
        + FC7 * als * (61. + t * ( t * (179. - t) - 479. ) )
        )));
    xy.y = P->k0 * (pj_mlfn(lp.phi, sinphi, cosphi, Q->en) - Q->ml0 +
        sinphi * al * lp.lam * FC2 * ( 1. +
        FC4 * als * (5. - t + n * (9. + 4. * n) +
        FC6 * als * (61. + t * (t - 58.) + n * (270. - 330 * t)
        + FC8 * als * (1385. + t * ( t * (543. - t) - 3111.) )
        ))));
    return (xy);
}

	LP /* inverse projection entry */
pj_inv(XY xy, PJ *P) {
	LP lp;

	/* can't do as much preliminary checking as with forward */
	if (xy.x == HUGE_VAL || xy.y == HUGE_VAL) {
		lp.lam = lp.phi = HUGE_VAL;
		pj_ctx_set_errno( P->ctx, -15);
                return lp;
	}

	errno = pj_errno = 0;
        P->ctx->last_errno = 0;

	xy.x = (xy.x * P->to_meter - P->x0) * P->ra; /* descale and de-offset */
	xy.y = (xy.y * P->to_meter - P->y0) * P->ra;

        //Check for NULL pointer
        if (P->inv != NULL)
        {
	    lp = (*P->inv)(xy, P); /* inverse project */
	    if (P->ctx->last_errno )
		lp.lam = lp.phi = HUGE_VAL;
	    else {
		lp.lam += P->lam0; /* reduce from del lp.lam */
		if (!P->over)
			lp.lam = adjlon(lp.lam); /* adjust longitude to CM */
		if (P->geoc && fabs(fabs(lp.phi)-HALFPI) > EPS)
			lp.phi = atan(P->one_es * tan(lp.phi));
	    }
        }
        else
        {
           lp.lam = lp.phi = HUGE_VAL;
        }
	return lp;
}

int proj_errno_reset (const PJ *P) {
/******************************************************************************
    Clears errno in the context and thread local levels
    through the low level pj_ctx interface.

    Returns the previous value of the errno, for convenient reset/restore
    operations:

    int foo (PJ *P) {
        // errno may be set on entry, but we need to reset it to be able to
        // check for errors from "do_something_with_P(P)"
        int last_errno = proj_errno_reset (P);

        // local failure
        if (0==P)
            return proj_errno_set (P, 42);

        // call to function that may fail
        do_something_with_P (P);

        // failure in do_something_with_P? - keep latest error status
        if (proj_errno(P))
            return proj_errno (P);

        // success - restore previous error status, return 0
        return proj_errno_restore (P, last_errno);
    }
******************************************************************************/
    int last_errno;
    last_errno = proj_errno (P);

    pj_ctx_set_errno (pj_get_ctx ((PJ *) P), 0);
    errno = 0;
    pj_errno = 0;
    return last_errno;
}

int nad_ctable_load( projCtx ctx, struct CTABLE *ct, FILE *fid )

{
    int  a_size;

    fseek( fid, sizeof(struct CTABLE), SEEK_SET );

    /* read all the actual shift values */
    a_size = ct->lim.lam * ct->lim.phi;
    ct->cvs = (FLP *) pj_malloc(sizeof(FLP) * a_size);
    if( ct->cvs == NULL 
        || fread(ct->cvs, sizeof(FLP), a_size, fid) != a_size )
    {
        pj_dalloc( ct->cvs );
        ct->cvs = NULL;

        pj_log( ctx, PJ_LOG_ERROR, 
                "ctable loading failed on fread() - binary incompatible?\n" );
        pj_ctx_set_errno( ctx, -38 );
        return 0;
    }

    return 1;
} 

static int pj_gridinfo_init_ntv1( projCtx ctx, PAFile fid, PJ_GRIDINFO *gi )

{
    unsigned char header[176];
    struct CTABLE *ct;
    LP		ur;

    assert( sizeof(int) == 4 );
    assert( sizeof(double) == 8 );
    if( sizeof(int) != 4 || sizeof(double) != 8 )
    {
        pj_log( ctx, PJ_LOG_ERROR,
                 "basic types of inappropraiate size in nad_load_ntv1()" );
        pj_ctx_set_errno( ctx, -38 );
        return 0;
    }

/* -------------------------------------------------------------------- */
/*      Read the header.                                                */
/* -------------------------------------------------------------------- */
    if( pj_ctx_fread( ctx, header, sizeof(header), 1, fid ) != 1 )
    {
        pj_ctx_set_errno( ctx, -38 );
        return 0;
    }

/* -------------------------------------------------------------------- */
/*      Regularize fields of interest.                                  */
/* -------------------------------------------------------------------- */
    if( IS_LSB )
    {
        swap_words( header+8, 4, 1 );
        swap_words( header+24, 8, 1 );
        swap_words( header+40, 8, 1 );
        swap_words( header+56, 8, 1 );
        swap_words( header+72, 8, 1 );
        swap_words( header+88, 8, 1 );
        swap_words( header+104, 8, 1 );
    }

    if( *((int *) (header+8)) != 12 )
    {
        pj_log( ctx, PJ_LOG_ERROR,
                "NTv1 grid shift file has wrong record count, corrupt?" );
        pj_ctx_set_errno( ctx, -38 );
        return 0;
    }

/* -------------------------------------------------------------------- */
/*      Fill in CTABLE structure.                                       */
/* -------------------------------------------------------------------- */
    ct = (struct CTABLE *) pj_malloc(sizeof(struct CTABLE));
    strcpy( ct->id, "NTv1 Grid Shift File" );

    ct->ll.lam = - *((double *) (header+72));
    ct->ll.phi = *((double *) (header+24));
    ur.lam = - *((double *) (header+56));
    ur.phi = *((double *) (header+40));
    ct->del.lam = *((double *) (header+104));
    ct->del.phi = *((double *) (header+88));
    ct->lim.lam = (int) (fabs(ur.lam-ct->ll.lam)/ct->del.lam + 0.5) + 1;
    ct->lim.phi = (int) (fabs(ur.phi-ct->ll.phi)/ct->del.phi + 0.5) + 1;

    pj_log( ctx, PJ_LOG_DEBUG_MINOR,
            "NTv1 %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)",
            ct->lim.lam, ct->lim.phi,
            ct->ll.lam, ct->ll.phi, ur.lam, ur.phi );

    ct->ll.lam *= DEG_TO_RAD;
    ct->ll.phi *= DEG_TO_RAD;
    ct->del.lam *= DEG_TO_RAD;
    ct->del.phi *= DEG_TO_RAD;
    ct->cvs = NULL;

    gi->ct = ct;
    gi->grid_offset = pj_ctx_ftell( ctx, fid );
    gi->format = "ntv1";

    return 1;
}