#include <stdio.h>
#include <math.h>

#include "tmp/c.h"
#include "tmp/postgres.h"

#include "utils/log.h"

#include "f2c.h"
#include "gctp.h"
#include "geoloc.h"

GEOLOC *
lazea_tostd(v)
	LAZEA *v;
{
	int i;
	GEOLOC *gv;
	gctpcvt cvt;
	gctpparam inspec;
	gctpparam outspec;

	/* get a return value */
	gv = (GEOLOC *) palloc(sizeof(GEOLOC));
	gv->geo_varlen = sizeof(GEOLOC);
	sprintf(&(gv->geo_name[0]), "latlon");

	cvt.cvt_insys = SYS_LAZEA;
	cvt.cvt_inunit = UNIT_METERS;
	cvt.cvt_inzone = 61;			/* XXX frew magic */

	cvt.cvt_outsys = SYS_LATLON;
	cvt.cvt_outunit = UNIT_DEGARC;
	cvt.cvt_outzone = 0;			/* XXX frew magic */

	inspec.parm_x = v->la_g.geo_x;
	inspec.parm_y = v->la_g.geo_y;

	for (i = 0; i < 15; i++) {
		outspec.parm_params[i] = 0.0;
		inspec.parm_params[i] = 0.0;
	}

	inspec.parm_params[0] = (double) v->la_sphcode;
	inspec.parm_params[4] = v->la_lon0;
	inspec.parm_params[5] = v->la_lat0;
	inspec.parm_params[6] = v->la_falseE; 
	inspec.parm_params[7] = v->la_falseN; 

	gctp(&cvt, &inspec, &outspec);

	gv->geo_prec = v->la_g.geo_prec;	/* XXX this is WRONG!!! */
	gv->geo_x = outspec.parm_x;
	gv->geo_y = outspec.parm_y;

	return (gv);
}

LAZEA *
lazea_fromstd(gv)
	GEOLOC *gv;
{
	int i;
	LAZEA *v;
	gctpcvt cvt;
	gctpparam inspec;
	gctpparam outspec;

	/* get a return value */
	v = (LAZEA *) palloc(sizeof(LAZEA));
	v->la_g.geo_varlen = sizeof(LAZEA);
	sprintf(&(v->la_g.geo_name[0]), "lazea");

	cvt.cvt_insys = SYS_LATLON;
	cvt.cvt_inunit = UNIT_DEGARC;
	cvt.cvt_inzone = 0;				/* XXX frew magic */

	cvt.cvt_outsys = SYS_LAZEA;
	cvt.cvt_outunit = UNIT_METERS;
	cvt.cvt_outzone = 61;				/* XXX frew magic */

	inspec.parm_x = gv->geo_x;
	inspec.parm_y = gv->geo_y;

	for (i = 0; i < 15; i++) {
		inspec.parm_params[i] = 0.0;
		outspec.parm_params[i] = 0.0;
	}

	/*
	 *  The output spec parameter vector shouldn't be hardwired,
	 *  but so far we haven't thought of a good way to allow users
	 *  to specify state vectors dynamically for the various projection
	 *  transformations that gctp does.  For the demo, we hardwire
	 *  spheroid code to 0 (Clarke 1866), origin to the center of the
	 *  conterminous US, and false eastings and northings to zero.
	 */

	outspec.parm_params[4] = -100000000.0;	/* origin lon */
	outspec.parm_params[5] =   45000000.0;	/* origin lat */

	gctp(&cvt, &inspec, &outspec);

	v->la_g.geo_prec = gv->geo_prec;	/* XXX this is WRONG!!! */
	v->la_g.geo_x = outspec.parm_x;
	v->la_g.geo_y = outspec.parm_y;

	v->la_sphcode = outspec.parm_params[0];
	v->la_lon0 = outspec.parm_params[4];
	v->la_lat0 = outspec.parm_params[5];
	v->la_falseE = outspec.parm_params[6];
	v->la_falseN = outspec.parm_params[7];

	return (v);
}

/*
 *  gctp() -- Wrapper routine to conceal FORTRAN parameter-passing
 *	      conventions from polite company.
 */

gctp(cvt, inspec, outspec)
	gctpcvt *cvt;
	gctpparam *inspec;
	gctpparam *outspec;
{
	integer iflg, ipr, jpr, lemsg, lparm;
	integer ln27, ln83;
	char fn27, fn83;
	integer spheroid;
	integer length;
	doublereal crdin[2];
	doublereal crdio[2];
	int status;

	crdin[0] = inspec->parm_x;
	crdin[1] = inspec->parm_y;

	spheroid = 0;				/* Clarke 1866 */
	ipr = 1;
	jpr = 1;
	lemsg = 0;
	lparm = 0;
	ln27 = ln83 = 0;
	fn27 = '\0'; fn83 = '\0';
	length = 0;

	/* just close your eyes and jump */
	status = gtpz0_(crdin, &cvt->cvt_insys, &cvt->cvt_inzone,
	       inspec->parm_params,
	       &cvt->cvt_inunit, &spheroid,
	       &ipr, &jpr, &lemsg, &lparm,
	       crdio, &cvt->cvt_outsys, &cvt->cvt_outzone,
	       outspec->parm_params, &cvt->cvt_outunit,
	       &ln27, &ln83, &fn27, &fn83, &length,
	       &iflg,
	       0, 0);

	outspec->parm_x = crdio[0];
	outspec->parm_y = crdio[1];
}

/************************************************************************
 *	From here down, the GCTP code was automatically translated	*
 *	from FORTRAN to C and then massaged by Frew to eliminate all	*
 *	dependencies on the FORTRAN libraries.  This code is opaque	*
 *	and you don't want to read it.					*
 ************************************************************************/

/* Common Block Declarations */

struct {
    integer ierror;
} errmz0_;

#define errmz0_1 errmz0_

struct {
    integer ipemsg, ipelun, ipparm, ipplun;
} prinz0_;

#define prinz0_1 prinz0_

struct {
    doublereal az, ez, esz, e0z, e1z, e2z, e3z, e4z;
} ellpz0_;

#define ellpz0_1 ellpz0_

struct {
    integer iproj;
} projz0_;

#define projz0_1 projz0_

struct {
    integer ispher, lu27, lu83, len, msys;
    char file27[32], file83[32];
} spcs_;

#define spcs_1 spcs_

struct {
    integer switch__[23];
} toggle_;

#define toggle_1 toggle_

struct {
    integer itype;
} pj02_;

#define pj02_1 pj02_

union {
    struct {
	doublereal a, lon0, x0, y0, c, e, es, ns, rh0;
    } _1;
    struct {
	doublereal a03, lon003, x003, y003, c, e03, es03, ns03, rh003;
    } _2;
} pj03_;

#define pj03_1 (pj03_._1)
#define pj03_2 (pj03_._2)

union {
    struct {
	doublereal a, lon0, x0, y0, e, f, ns, rh0;
    } _1;
    struct {
	doublereal a04, lon004, x004, y004, e04, f04, ns04, rh004;
    } _2;
} pj04_;

#define pj04_1 (pj04_._1)
#define pj04_2 (pj04_._2)

union {
    struct {
	doublereal a, lon0, x0, y0, e, m1;
    } _1;
    struct {
	doublereal a05, lon005, x005, y005, e05, m1;
    } _2;
} pj05_;

#define pj05_1 (pj05_._1)
#define pj05_2 (pj05_._2)

union {
    struct {
	doublereal a, lon0, x0, y0, e, e4, fac, mcs, tcs;
	integer ind;
    } _1;
    struct {
	doublereal a06, lon006, x006, y006, e06, e4, fac, mcs, tcs;
	integer ind06;
    } _2;
} pj06_;

#define pj06_1 (pj06_._1)
#define pj06_2 (pj06_._2)

union {
    struct {
	doublereal a, lon0, x0, y0, e, e0, e1, e2, e3, es, ml0;
    } _1;
    struct {
	doublereal a07, lon007, x007, y007, e07, e007, e107, e207, e307, es07,
		 ml007;
    } _2;
} pj07_;

#define pj07_1 (pj07_._1)
#define pj07_2 (pj07_._2)

union {
    struct {
	doublereal a, lon0, x0, y0, e0, e1, e2, e3, gl, ns, rh0;
    } _1;
    struct {
	doublereal a08, lon008, x008, y008, e008, e108, e208, e308, gl, ns08, 
		rh008;
    } _2;
} pj08_;

#define pj08_1 (pj08_._1)
#define pj08_2 (pj08_._2)

union {
    struct {
	doublereal a, lon0, x0, y0, es, esp, e0, e1, e2, e3, ks0, lat0, ml0;
	integer ind;
    } _1;
    struct {
	doublereal a09, lon009, x009, y009, es09, esp, e009, e109, e209, e309,
		 ks009, lat009, ml009;
	integer ind09;
    } _2;
} pj09_;

#define pj09_1 (pj09_._1)
#define pj09_2 (pj09_._2)

union {
    struct {
	doublereal a, lon0, x0, y0, cosph0, lat0, sinph0;
    } _1;
    struct {
	doublereal a10, lon010, x010, y010, cosp10, lat010, sinp10;
    } _2;
} pj10_;

#define pj10_1 (pj10_._1)
#define pj10_2 (pj10_._2)

union {
    struct {
	doublereal a, lon0, x0, y0, cosph0, lat0, sinph0;
    } _1;
    struct {
	doublereal a11, lon011, x011, y011, cosp11, lat011, sinp11;
    } _2;
} pj11_;

#define pj11_1 (pj11_._1)
#define pj11_2 (pj11_._2)

union {
    struct {
	doublereal a, lon0, x0, y0, cosph0, lat0, sinph0;
    } _1;
    struct {
	doublereal a12, lon012, x012, y012, cosp12, lat012, sinp12;
    } _2;
} pj12_;

#define pj12_1 (pj12_._1)
#define pj12_2 (pj12_._2)

union {
    struct {
	doublereal a, lon0, x0, y0, cosph0, lat0, sinph0;
    } _1;
    struct {
	doublereal a13, lon013, x013, y013, cosp13, lat013, sinp13;
    } _2;
} pj13_;

#define pj13_1 (pj13_._1)
#define pj13_2 (pj13_._2)

union {
    struct {
	doublereal a, lon0, x0, y0, cosph0, lat0, sinph0;
    } _1;
    struct {
	doublereal a14, lon014, x014, y014, cosp14, lat014, sinp14;
    } _2;
} pj14_;

#define pj14_1 (pj14_._1)
#define pj14_2 (pj14_._2)

union {
    struct {
	doublereal a, lon0, x0, y0, cosph0, lat0, p, sinph0;
    } _1;
    struct {
	doublereal a15, lon015, x015, y015, cosp15, lat015, p, sinp15;
    } _2;
} pj15_;

#define pj15_1 (pj15_._1)
#define pj15_2 (pj15_._2)

union {
    struct {
	doublereal a, lon0, x0, y0;
    } _1;
    struct {
	doublereal a16, lon016, x016, y016;
    } _2;
} pj16_;

#define pj16_1 (pj16_._1)
#define pj16_2 (pj16_._2)

union {
    struct {
	doublereal a, lon0, x0, y0, lat1;
    } _1;
    struct {
	doublereal a17, lon017, x017, y017, lat1;
    } _2;
} pj17_;

#define pj17_1 (pj17_._1)
#define pj17_2 (pj17_._2)

union {
    struct {
	doublereal a, lon0, x0, y0;
    } _1;
    struct {
	doublereal a18, lon018, x018, y018;
    } _2;
} pj18_;

#define pj18_1 (pj18_._1)
#define pj18_2 (pj18_._2)

union {
    struct {
	doublereal a, lon0, x0, y0;
    } _1;
    struct {
	doublereal a19, lon019, x019, y019;
    } _2;
} pj19_;

#define pj19_1 (pj19_._1)
#define pj19_2 (pj19_._2)

union {
    struct {
	doublereal lon0, x0, y0, al, bl, cosalf, cosgam, e, el, sinalf, 
		singam, u0;
    } _1;
    struct {
	doublereal lon020, x020, y020, al, bl, cosalf, cosgam, e20, el, 
		sinalf, singam, u0;
    } _2;
} pj20_;

#define pj20_1 (pj20_._1)
#define pj20_2 (pj20_._2)

union {
    struct {
	doublereal a, lon0, x0, y0, pr[20], xlr[20];
    } _1;
    struct {
	doublereal a21, lon021, x021, y021, pr[20], xlr[20];
    } _2;
} pj21_;

#define pj21_1 (pj21_._1)
#define pj21_2 (pj21_._2)

union {
    struct {
	doublereal q, t, u, w, es, p22, sa, ca, xj;
    } _1;
    struct {
	doublereal q, t, u, w, es22, p22, sa, ca, xj;
    } _2;
} norm_;

#define norm_1 (norm_._1)
#define norm_2 (norm_._2)

union {
    struct {
	doublereal a, x0, y0, a2, a4, b, c1, c3;
	integer land, path;
    } _1;
    struct {
	doublereal a22, x022, y022, a2, a4, b, c1, c3;
	integer land, path;
    } _2;
} pj22_;

#define pj22_1 (pj22_._1)
#define pj22_2 (pj22_._2)

union {
    struct {
	doublereal a, lon0, x0, y0, acoef[6], bcoef[6], ec, lat0, cchio, 
		schio;
	integer n;
    } _1;
    struct {
	doublereal a23, lon023, x023, y023, acoef[6], bcoef[6], ec, lat023, 
		cchio, schio;
	integer n;
    } _2;
} pj23_;

#define pj23_1 (pj23_._1)
#define pj23_2 (pj23_._2)

struct {
    doublereal azz;
} sphrz0_;

#define sphrz0_1 sphrz0_

/* Table of constant values */

static integer c__0 = 0;
static doublereal c_b34 = 4e7;
static shortint cs__1 = 1;
static shortint cs__0 = 0;
static doublereal c_b239 = 1.;

/*     GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE - VERSION 2.0.0 */
/*     FORTRAN 77 LANGUAGE FOR IBM, AMDAHL, GOULD, VAX, AND CONCURRENT */
/*     (P-E) COMPUTERS */
/*                   ADJLZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

doublereal adjlz0_(lon)
doublereal *lon;
{
    /* Initialized data */

    static doublereal two = 2.;
    static doublereal pi = 3.14159265358979323846;

    /* System generated locals */
    doublereal ret_val;

    /* Builtin functions */
    double d_sign();

    /* Local variables */
    static doublereal twopi;


/* FUNCTION TO ADJUST LONGITUDE ANGLE TO MODULO 180 DEGREES. */


L20:
    ret_val = *lon;
    if (abs(*lon) <= pi) {
	return ret_val;
    }
    twopi = two * pi;
    *lon -= d_sign(&twopi, lon);
    goto L20;

} /* adjlz0_ */

/*                   ASINZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

doublereal asinz0_(con)
doublereal *con;
{
    /* Initialized data */

    static doublereal one = 1.;

    /* System generated locals */
    doublereal ret_val;

    /* Builtin functions */
    double d_sign(), asin();


/* THIS FUNCTION ADJUSTS FOR ROUND-OFF ERRORS IN COMPUTING ARCSINE */


    if (abs(*con) > one) {
	*con = d_sign(&one, con);
    }
    ret_val = asin(*con);
    return ret_val;

} /* asinz0_ */

/*                   DMSPZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

doublereal dmspz0_(sgna, degs, mins, secs, sgna_len)
char *sgna;
integer *degs, *mins;
real *secs;
ftnlen sgna_len;
{
    /* Initialized data */

    static doublereal con1 = 1e6;
    static doublereal con2 = 1e3;
    static char neg[1+1] = "-";

    /* System generated locals */
    doublereal ret_val;

    /* Local variables */
    static doublereal con;


/* SUBROUTINE TO CONVERT UNPACKED DMS TO PACKED DMS ANGLE */
/* SGNA : SIGN OF ANGLE */
/* DEGS : DEGREES PORTION OF ANGLE */
/* MINS : MINUTES PORTION OF ANGLE */
/* SECS : SECONDS PORTION OF ANGLE */


    con = (doublereal) (*degs) * con1 + (doublereal) (*mins) * con2 + (
	    doublereal) (*secs);
    if (*sgna == neg[0]) {
	con = -con;
    }
    ret_val = con;
    return ret_val;

} /* dmspz0_ */

/*                   E0FNZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

doublereal e0fnz0_(eccnts)
doublereal *eccnts;
{
    /* Initialized data */

    static doublereal one = 1.;
    static doublereal oneq = 1.25;
    static doublereal three = 3.;
    static doublereal sixt = 16.;
    static doublereal quart = .25;

    /* System generated locals */
    doublereal ret_val;


/* FUNCTION TO COMPUTE CONSTANT (E0). */


    ret_val = one - quart * *eccnts * (one + *eccnts / sixt * (three + oneq * 
	    *eccnts));

    return ret_val;
} /* e0fnz0_ */

/*                   E1FNZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

doublereal e1fnz0_(eccnts)
doublereal *eccnts;
{
    /* Initialized data */

    static doublereal con1 = .375;
    static doublereal con2 = .25;
    static doublereal con3 = .46875;
    static doublereal one = 1.;

    /* System generated locals */
    doublereal ret_val;


/* FUNCTION TO COMPUTE CONSTANT (E1). */


    ret_val = con1 * *eccnts * (one + con2 * *eccnts * (one + con3 * *eccnts))
	    ;

    return ret_val;
} /* e1fnz0_ */

/*                   E2FNZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

doublereal e2fnz0_(eccnts)
doublereal *eccnts;
{
    /* Initialized data */

    static doublereal con1 = .05859375;
    static doublereal con2 = .75;
    static doublereal one = 1.;

    /* System generated locals */
    doublereal ret_val;


/* FUNCTION TO COMPUTE CONSTANT (E2). */


    ret_val = con1 * *eccnts * *eccnts * (one + con2 * *eccnts);

    return ret_val;
} /* e2fnz0_ */

/*                   E3FNZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

doublereal e3fnz0_(eccnts)
doublereal *eccnts;
{
    /* System generated locals */
    doublereal ret_val;


/* FUNCTION TO COMPUTE CONSTANT (E3). */


    ret_val = *eccnts * *eccnts * *eccnts * .011393229166666666;

    return ret_val;
} /* e3fnz0_ */

/*                   E4FNZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

doublereal e4fnz0_(eccent)
doublereal *eccent;
{
    /* Initialized data */

    static doublereal one = 1.;

    /* System generated locals */
    doublereal ret_val;

    /* Builtin functions */
    double pow_dd(), sqrt();

    /* Local variables */
    static doublereal com, con;


/* FUNCTION TO COMPUTE CONSTANT (E4). */


    con = one + *eccent;
    com = one - *eccent;
    ret_val = sqrt(pow_dd(&con, &con) * pow_dd(&com, &com));

    return ret_val;
} /* e4fnz0_ */

/*                   GTPZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

/* Subroutine */ int gtpz0_(crdin, insys, inzone, tparin, inunit, insph, ipr, 
	jpr, lemsg, lparm, crdio, iosys, iozone, tpario, iounit, ln27, ln83, 
	fn27, fn83, length, iflg, fn27_len, fn83_len)
doublereal *crdin;
integer *insys, *inzone;
doublereal *tparin;
integer *inunit, *insph, *ipr, *jpr, *lemsg, *lparm;
doublereal *crdio;
integer *iosys, *iozone;
doublereal *tpario;
integer *iounit, *ln27, *ln83;
char *fn27, *fn83;
integer *length, *iflg;
ftnlen fn27_len;
ftnlen fn83_len;
{
    /* Initialized data */

    static integer sysunt[24] = { 0,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
	    2,2 };
    static doublereal pdin[15] = { 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	    0. };
    static doublereal pdio[15] = { 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	    0. };
    static integer insp = 999;
    static integer inpj = 999;
    static integer inzn = 99999;
    static integer iosp = 999;
    static integer iopj = 999;
    static integer iozn = 99999;
    static integer iter = 0;
    static integer jflag = 0;
    static integer nad27[134] = { 101,102,5010,5300,201,202,203,301,302,401,
	    402,403,404,405,406,407,501,502,503,600,700,901,902,903,1001,1002,
	    5101,5102,5103,5104,5105,1101,1102,1103,1201,1202,1301,1302,1401,
	    1402,1501,1502,1601,1602,1701,1702,1703,1801,1802,1900,2001,2002,
	    2101,2102,2103,2111,2112,2113,2201,2202,2203,2301,2302,2401,2402,
	    2403,2501,2502,2503,2601,2602,2701,2702,2703,2800,2900,3001,3002,
	    3003,3101,3102,3103,3104,3200,3301,3302,3401,3402,3501,3502,3601,
	    3602,3701,3702,3800,3901,3902,4001,4002,4100,4201,4202,4203,4204,
	    4205,4301,4302,4303,4400,4501,4502,4601,4602,4701,4702,4801,4802,
	    4803,4901,4902,4903,4904,5001,5002,5003,5004,5005,5006,5007,5008,
	    5009,5201,5202,5400 };
    static integer nad83[134] = { 101,102,5010,5300,201,202,203,301,302,401,
	    402,403,404,405,406,0,501,502,503,600,700,901,902,903,1001,1002,
	    5101,5102,5103,5104,5105,1101,1102,1103,1201,1202,1301,1302,1401,
	    1402,1501,1502,1601,1602,1701,1702,1703,1801,1802,1900,2001,2002,
	    2101,2102,2103,2111,2112,2113,2201,2202,2203,2301,2302,2401,2402,
	    2403,2500,0,0,2600,0,2701,2702,2703,2800,2900,3001,3002,3003,3101,
	    3102,3103,3104,3200,3301,3302,3401,3402,3501,3502,3601,3602,3701,
	    3702,3800,3900,0,4001,4002,4100,4201,4202,4203,4204,4205,4301,
	    4302,4303,4400,4501,4502,4601,4602,4701,4702,4801,4802,4803,4901,
	    4902,4903,4904,5001,5002,5003,5004,5005,5006,5007,5008,5009,5200,
	    0,5400 };
    static integer nadut[134] = { 2,2,2,2,5,5,5,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
	    2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
	    2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
	    2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
	    2,2,2,2,2,2,2,2,2,2,2,2,2 };
    static integer sptype[134] = { 9,9,4,4,9,9,9,4,4,4,4,4,4,4,4,4,4,4,4,4,9,
	    9,9,4,9,9,9,9,9,9,9,9,9,9,9,9,9,9,4,4,4,4,4,4,4,4,4,9,9,4,4,4,9,9,
	    9,4,4,4,4,4,4,9,9,9,9,9,4,4,4,4,4,9,9,9,9,9,9,9,9,9,9,9,4,4,4,4,4,
	    4,4,4,4,4,4,4,9,4,4,4,4,4,4,4,4,4,4,4,4,4,9,4,4,4,4,4,4,4,4,4,9,9,
	    9,9,20,9,9,9,9,9,9,9,9,4,4,7 };

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    /* Subroutine */ int s_copy();
    double d_sign();

    /* Local variables */
    extern /* Subroutine */ int pj01z0_(), pj02z0_(), pj03z0_(), pj04z0_(), 
	    pj05z0_(), pj06z0_(), pj07z0_(), pj08z0_(), pj09z0_(), pj10z0_(), 
	    pj11z0_(), pj12z0_(), pj13z0_(), pj14z0_(), pj15z0_(), pj16z0_(), 
	    pj17z0_(), pj18z0_(), pj19z0_(), pj20z0_();
    static integer i;
    extern /* Subroutine */ int pj21z0_(), pj22z0_(), pj23z0_();
    static shortint inmod, iomod;
    static doublereal coord[2];
    static integer iosph, iunit;
    static doublereal dummy[15];
    extern doublereal adjlz0_();
    extern /* Subroutine */ int sphdz0_();
    static integer id;
    extern /* Subroutine */ int untfz0_();
    static doublereal factor;
    static integer inspcs, iospcs;
    extern /* Subroutine */ int pjinit_();
    static integer ind;


/* ********************************************************************** 
*/

/* INPUT **************************************************************** 
*/
/* CRDIN  : COORDINATES IN INPUT SYSTEM (2 DP WORDS ARRAY). */
/* INSYS  : CODE NUMBER OF INPUT COORDINATE SYSTEM (INTEGER). */
/*            =  0 , GEOGRAPHIC */
/*            =  1 , U T M */
/*            =  2 , STATE PLANE */
/*            =  3 , ALBERS CONICAL EQUAL-AREA */
/*            =  4 , LAMBERT CONFORMAL CONIC */
/*            =  5 , MERCATOR */
/*            =  6 , POLAR STEREOGRAPHIC */
/*            =  7 , POLYCONIC */
/*            =  8 , EQUIDISTANT CONIC */
/*            =  9 , TRANSVERSE MERCATOR */
/*            = 10 , STEREOGRAPHIC */
/*            = 11 , LAMBERT AZIMUTHAL EQUAL-AREA */
/*            = 12 , AZIMUTHAL EQUIDISTANT */
/*            = 13 , GNOMONIC */
/*            = 14 , ORTHOGRAPHIC */
/*            = 15 , GENERAL VERTICAL NEAR-SIDE PERSPECTIVE */
/*            = 16 , SINUSOIDAL */
/*            = 17 , EQUIRECTANGULAR (PLATE CARREE) */
/*            = 18 , MILLER CYLINDRICAL */
/*            = 19 , VAN DER GRINTEN I */
/*            = 20 , OBLIQUE MERCATOR (HOTINE) */
/*            = 21 , ROBINSON */
/*            = 22 , SPACE OBLIQUE MERCATOR */
/*            = 23 , MODIFIED-STEREOGRAPHIC CONFORMAL (ALASKA) */
/* INZONE : CODE NUMBER OF INPUT COORDINATE ZONE (INTEGER). */
/* TPARIN : PARAMETERS OF INPUT REFERENCE SYSTEM (15 DP WORDS ARRAY). */
/* INUNIT : CODE NUMBER OF UNITS OF MEASURE FOR INPUT COORDINATES (I*4) */

/*            = 0 , RADIANS. */
/*            = 1 , U.S. FEET. */
/*            = 2 , METERS. */
/*            = 3 , SECONDS OF ARC. */
/*            = 4 , DEGREES OF ARC. */
/*            = 5 , INTERNATIONAL FEET. */
/*            = 6 , USE LEGISLATED DISTANCE UNITS FROM NADUT TABLE */
/* INSPH  : INPUT SPHEROID CODE.  SEE SPHDZ0 FOR PROPER CODES. */
/* IPR    : PRINTOUT FLAG FOR ERROR MESSAGES. 0=YES, 1=NO */
/* JPR    : PRINTOUT FLAG FOR PROJECTION PARAMETERS 0=YES, 1=NO */
/* LEMSG  : LOGICAL UNIT FOR LISTING ERROR MESSAGES IF IPR = 0 */
/* LPARM  : LOGICAL UNIT FOR LISTING PROJECTION PARAMETERS IF JPR = 0 */
/* LN27   : LOGICAL UNIT FOR NAD 1927 SPCS PARAMETER FILE */
/* FN27   : FILE NAME OF NAD 1927 SPCS PARAMETERS */
/* LN83   : LOGICAL UNIT FOR NAD 1983 SPCS PARAMETER FILE */
/* FN83   : FILE NAME OF NAD 1983 SPCS PARAMETERS */
/* LENGTH : RECORD LENGTH OF NAD1927 AND NAD1983 PARAMETER FILES */
/* OUTPUT ***                                                       ***** 
*/
/* IOSYS  : CODE NUMBER OF OUTPUT COORDINATE SYSTEM (INTEGER). */
/* IOZONE : CODE NUMBER OF OUTPUT COORDINATE ZONE (INTEGER). */
/* TPARIO : PARAMETERS OF OUTPUT REFERENCE SYSTEM (15 DP WORDS ARRAY). */
/* IOUNIT : CODE NUMBER OF UNITS OF MEASURE FOR OUTPUT COORDINATES (I*4) 
*/
/* CRDIO  : COORDINATES IN OUTPUT REFERENCE SYSTEM (2 DP WORDS ARRAY). */
/* IFLG   : RETURN FLAG (INTEGER). */
/*            = 0 , SUCCESSFUL TRANSFORMATION. */
/*            = 1 , ILLEGAL INPUT SYSTEM CODE. */
/*            = 2 , ILLEGAL OUTPUT SYSTEM CODE. */
/*            = 3 , ILLEGAL INPUT UNIT CODE. */
/*            = 4 , ILLEGAL OUTPUT UNIT CODE. */
/*            = 5 , INCONSISTENT UNIT AND SYSTEM CODES FOR INPUT. */
/*            = 6 , INCONSISTENT UNIT AND SYSTEM CODES FOR OUTPUT. */
/*            = 7 , ILLEGAL INPUT ZONE CODE. */
/*            = 8 , ILLEGAL OUTPUT ZONE CODE. */
/*      OTHERWISE , ERROR CODE FROM PROJECTION COMPUTATIONAL MODULE. */


    /* Parameter adjustments */
    --tpario;
    --crdio;
    --tparin;
    --crdin;

    /* Function Body */



/*    TABLE OF UNIT CODES AS SPECIFIED BY STATE LAWS AS OF 7/18/89 */
/*    FOR NAD 1983 SPCS - 1 = U.S. FOOT, 2 = METER, 5 = INTERNATIONAL FT. 
*/


/*     TABLE OF STATE PLANE ZONE TYPES:  4 = LAMBERT, 7 = POLYCONIC, */
/*     9 = TRANSVERSE MERCATOR, AND 20 = OBLIQUE MERCATOR */


/*     SETUP */

    iosph = *insph;
    prinz0_1.ipemsg = *ipr;
    prinz0_1.ipparm = *jpr;
    prinz0_1.ipelun = *lemsg;
    prinz0_1.ipplun = *lparm;
    projz0_1.iproj = *insys;
    spcs_1.lu27 = *ln27;
    s_copy(spcs_1.file27, fn27, 32L, 32L);
    spcs_1.lu83 = *ln83;
    s_copy(spcs_1.file83, fn83, 32L, 32L);
    spcs_1.len = *length;

/*     INITIALIZE SWITCH FOR EACH PROJECTION TO ZERO */

    ++iter;
    if (iter <= 1) {
	for (i = 1; i <= 23; ++i) {
	    toggle_1.switch__[i - 1] = 0;
/* L10: */
	}
	spcs_1.msys = 2;
    }
    inspcs = 2;
    iospcs = 2;
    if (jflag != 0) {
	goto L11;
    }
    ellpz0_1.ez = 0.;
    ellpz0_1.esz = 0.;
    sphdz0_(&c__0, dummy);
    jflag = 1;

/*     CHECK FOR REPEAT OF INPUT SYSTEM */

L11:
    inmod = 1;
    if (*insph != insp && *insys > 0) {
	toggle_1.switch__[*insys - 1] = 0;
    }
    if (*insph != iosp && *insys > 0) {
	toggle_1.switch__[*insys - 1] = 0;
    }
    if (*insys == 2) {
	if (*inzone > 0) {
	    id = 0;
	    if (*insph == 0) {
		for (i = 1; i <= 134; ++i) {
		    if (*inzone == nad27[i - 1]) {
			id = i;
		    }
/* L12: */
		}
	    }
	    if (*insph == 8) {
		for (i = 1; i <= 134; ++i) {
		    if (*inzone == nad83[i - 1]) {
			id = i;
		    }
/* L13: */
		}
	    }
	    if (id != 0) {
		inspcs = sptype[id - 1];
	    }
	    if (*inzone != toggle_1.switch__[inspcs - 1]) {
		goto L15;
	    }
	}
    }
    if (insp != *insph) {
	goto L15;
    }
    if (inpj != *insys) {
	goto L15;
    }
    if (inzn != *inzone) {
	goto L15;
    }
    if (*insys >= 3) {
	for (i = 1; i <= 15; ++i) {
	    if (tparin[i] != pdin[i - 1]) {
		goto L15;
	    }
/* L14: */
	}
    }
    inmod = 0;
    goto L30;

/*     SAVE INPUT SYSTEM PARAMETERS */

L15:
    insp = *insph;
    inpj = *insys;
    inzn = *inzone;
    for (i = 1; i <= 15; ++i) {
/* L16: */
	pdin[i - 1] = tparin[i];
    }

/* CHECK VALIDITY OF CODES FOR UNITS OF MEASURE AND REFERENCE SYSTEMS. */

    if (*insys < 0 || *insys > 24) {
	*iflg = 1;
	return 0;
    }
    if (*inunit < 0 || *inunit > 6) {
	*iflg = 3;
	return 0;
    }

/*     CHECK FOR REPEAT OF OUTPUT SYSTEM */

L30:
    iomod = 1;
    if (*insph != iosp && *iosys > 0) {
	toggle_1.switch__[*iosys - 1] = 0;
    }
    if (*insph != insp && *iosys > 0) {
	toggle_1.switch__[*iosys - 1] = 0;
    }
    if (*iosys == 2) {
	if (*iozone > 0) {
	    id = 0;
	    if (iosph == 0) {
		for (i = 1; i <= 134; ++i) {
		    if (*iozone == nad27[i - 1]) {
			id = i;
		    }
/* L32: */
		}
	    }
	    if (iosph == 8) {
		for (i = 1; i <= 134; ++i) {
		    if (*iozone == nad83[i - 1]) {
			id = i;
		    }
/* L33: */
		}
	    }
	    if (id != 0) {
		iospcs = sptype[id - 1];
	    }
	    if (*iozone != toggle_1.switch__[inspcs - 1]) {
		goto L35;
	    }
	}
    }
    if (iosp != *insph) {
	goto L35;
    }
    if (iosp != iosph) {
	goto L35;
    }
    if (iopj != *iosys) {
	goto L35;
    }
    if (iozn != *iozone) {
	goto L35;
    }
    if (*iosys >= 3) {
	for (i = 1; i <= 15; ++i) {
	    if (tpario[i] != pdio[i - 1]) {
		goto L35;
	    }
/* L34: */
	}
    }
    iomod = 0;
    goto L80;

/*     SAVE OUTPUT SYSTEM PARAMETERS */

L35:
    iosp = *insph;
    iopj = *iosys;
    iozn = *iozone;
    for (i = 1; i <= 15; ++i) {
/* L36: */
	pdio[i - 1] = tpario[i];
    }

/* CHECK VALIDITY OF CODES FOR UNITS OF MEASURE AND REFERENCE SYSTEMS. */

    if (*iosys < 0 || *iosys > 24) {
	*iflg = 2;
	return 0;
    }
    if (*iounit < 0 || *iounit > 6) {
	*iflg = 4;
	return 0;
    }

/* CHECK CONSISTENCY BETWEEN UNITS OF MEASURE AND REFERENCE SYSTEM. */

L80:
    iunit = sysunt[*insys];

/*     CHANGE UNITS TO LEGISLATED UNITS USING TABLE */

    if (*insph == 0 && *insys == 2 && *inunit == 6) {
	*inunit = 1;
    }
    if (*insph == 8 && *insys == 2 && *inunit == 6) {
	ind = 0;
	for (i = 1; i <= 134; ++i) {
	    if (*inzone == nad83[i - 1]) {
		ind = i;
	    }
/* L90: */
	}
	if (ind != 0) {
	    *inunit = nadut[ind - 1];
	}
    }
    untfz0_(inunit, &iunit, &factor, iflg);
    if (*iflg == 0) {
	goto L100;
    }
    *iflg = 5;
    return 0;
L100:
    coord[0] = factor * crdin[1];
    coord[1] = factor * crdin[2];
    iunit = sysunt[*iosys];

/*     CHANGE UNITS TO LEGISLATED UNITS USING TABLE */

    if (*insph == 0 && *iosys == 2 && *iounit == 6) {
	*iounit = 1;
    }
    if (*insph == 8 && *iosys == 2 && *iounit == 6) {
	ind = 0;
	for (i = 1; i <= 134; ++i) {
	    if (*iozone == nad83[i - 1]) {
		ind = i;
	    }
/* L110: */
	}
	if (ind != 0) {
	    *iounit = nadut[ind - 1];
	}
    }
    untfz0_(&iunit, iounit, &factor, iflg);
    if (*iflg == 0) {
	goto L120;
    }
    *iflg = 6;
    return 0;
L120:
    if (*insys != *iosys || *inzone != *iozone || *inzone <= 0) {
	goto L140;
    }
    crdio[1] = factor * coord[0];
    crdio[2] = factor * coord[1];
    return 0;

/* COMPUTE TRANSFORMED COORDINATES AND ADJUST THEIR UNITS. */

L140:
    if (*insys == 0) {
	goto L520;
    }
    if (*inzone > 60 || *insys == 1) {
	goto L200;
    }
    *iflg = 7;
    return 0;

/* INVERSE TRANSFORMATION. */

L200:
    projz0_1.iproj = *insys;
    spcs_1.ispher = *insph;
    if (*insys >= 3) {
	sphdz0_(insph, &tparin[1]);
    }

/*     CHECK FOR CHANGE IN ZONE FROM LAST USE OF THE INPUT PROJECTION */

    if (*insys == 1 && *inzone != toggle_1.switch__[8]) {
	toggle_1.switch__[0] = 0;
	inmod = 1;
    }
    if (*insys == 2 && *inzone != toggle_1.switch__[inspcs - 1]) {
	toggle_1.switch__[1] = 0;
	inmod = 1;
    }
    if (*inzone != toggle_1.switch__[*insys - 1]) {
	toggle_1.switch__[*insys - 1] = 0;
	inmod = 1;
    }

    if (*insys == 1) {
	if (*inzone == 0 && tparin[1] != 0.) {
	    goto L211;
	}
	tparin[1] = (doublereal) (*inzone * 6 - 183) * 1e6;
	d__1 = (doublereal) (*inzone);
	tparin[2] = d_sign(&c_b34, &d__1);
L211:
	sphdz0_(insph, dummy);
	tparin[14] = dummy[0];
	tparin[15] = dummy[1];
	if (inmod != 0) {
	    pjinit_(insys, inzone, &tparin[1]);
	    if (errmz0_1.ierror != 0) {
		inzn = 99999;
	    }
	    if (errmz0_1.ierror != 0) {
		goto L500;
	    }
	}
	pj01z0_(coord, &crdio[1], &cs__1);
    }

    if (*insys > 1) {
	if (inmod != 0) {
	    spcs_1.msys = inspcs;
	    pjinit_(insys, inzone, &tparin[1]);
	    if (errmz0_1.ierror != 0) {
		inzn = 99999;
	    }
	    if (errmz0_1.ierror != 0) {
		goto L500;
	    }
	}
	if (*insys == 2) {
	    pj02z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 3) {
	    pj03z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 4) {
	    pj04z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 5) {
	    pj05z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 6) {
	    pj06z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 7) {
	    pj07z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 8) {
	    pj08z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 9) {
	    pj09z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 10) {
	    pj10z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 11) {
	    pj11z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 12) {
	    pj12z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 13) {
	    pj13z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 14) {
	    pj14z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 15) {
	    pj15z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 16) {
	    pj16z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 17) {
	    pj17z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 18) {
	    pj18z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 19) {
	    pj19z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 20) {
	    pj20z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 21) {
	    pj21z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 22) {
	    pj22z0_(coord, &crdio[1], &cs__1);
	}
	if (*insys == 23) {
	    pj23z0_(coord, &crdio[1], &cs__1);
	}
    }

L500:
    *iflg = errmz0_1.ierror;
    for (i = 1; i <= 15; ++i) {
/* L510: */
	tparin[i] = pdin[i - 1];
    }
    if (*iflg != 0) {
	return 0;
    }
    crdio[1] = adjlz0_(&crdio[1]);
    if (*iosys == 0) {
	goto L920;
    }
    coord[0] = crdio[1];
    coord[1] = crdio[2];
L520:
    if (*insys == 0 && *iosys == 0) {
	crdio[1] = coord[0];
	crdio[2] = coord[1];
	goto L920;
    }
    if (*iozone > 60 || *iosys == 1) {
	goto L540;
    }
    *iflg = 8;
    return 0;

/* FORWARD TRANSFORMATION. */

L540:
    projz0_1.iproj = *iosys;
    spcs_1.ispher = *insph;
    if (*iosys >= 3) {
	sphdz0_(insph, &tpario[1]);
    }

/*     CHECK FOR CHANGE IN ZONE FROM LAST USE OF THE OUTPUT PROJECTION */

    if (*iosys == 1 && *iozone != toggle_1.switch__[8]) {
	toggle_1.switch__[0] = 0;
	iomod = 1;
    }
    if (*iosys == 2 && *iozone != toggle_1.switch__[iospcs - 1]) {
	toggle_1.switch__[1] = 0;
	iomod = 1;
    }
    if (*iozone != toggle_1.switch__[*iosys - 1]) {
	toggle_1.switch__[*iosys - 1] = 0;
	iomod = 1;
    }

    if (*iosys == 1) {
	tpario[1] = coord[0];
	tpario[2] = coord[1];
	sphdz0_(insph, dummy);
	tpario[14] = dummy[0];
	tpario[15] = dummy[1];
	if (iomod != 0) {
	    pjinit_(iosys, iozone, &tpario[1]);
	    if (errmz0_1.ierror != 0) {
		iozn = 99999;
	    }
	    if (errmz0_1.ierror != 0) {
		goto L900;
	    }
	    inzn = *iozone;
	}
	pj01z0_(coord, &crdio[1], &cs__0);
    }

    if (*iosys > 1) {
	if (iomod != 0) {
	    spcs_1.msys = iospcs;
	    pjinit_(iosys, iozone, &tpario[1]);
	    if (errmz0_1.ierror != 0) {
		iozn = 99999;
	    }
	    if (errmz0_1.ierror != 0) {
		goto L900;
	    }
	    inzn = *iozone;
	}
	if (*iosys == 2) {
	    pj02z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 3) {
	    pj03z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 4) {
	    pj04z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 5) {
	    pj05z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 6) {
	    pj06z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 7) {
	    pj07z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 8) {
	    pj08z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 9) {
	    pj09z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 10) {
	    pj10z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 11) {
	    pj11z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 12) {
	    pj12z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 13) {
	    pj13z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 14) {
	    pj14z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 15) {
	    pj15z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 16) {
	    pj16z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 17) {
	    pj17z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 18) {
	    pj18z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 19) {
	    pj19z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 20) {
	    pj20z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 21) {
	    pj21z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 22) {
	    pj22z0_(coord, &crdio[1], &cs__0);
	}
	if (*iosys == 23) {
	    pj23z0_(coord, &crdio[1], &cs__0);
	}
    }

L900:
    *iflg = errmz0_1.ierror;
    for (i = 1; i <= 15; ++i) {
/* L910: */
	tpario[i] = pdio[i - 1];
    }
L920:
    crdio[1] = factor * crdio[1];
    crdio[2] = factor * crdio[2];
    return 0;

} /* gtpz0_ */

/*                   MLFNZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

doublereal mlfnz0_(e0, e1, e2, e3, phi)
doublereal *e0, *e1, *e2, *e3, *phi;
{
    /* Initialized data */

    static doublereal two = 2.;
    static doublereal four = 4.;
    static doublereal six = 6.;

    /* System generated locals */
    doublereal ret_val;

    /* Builtin functions */
    double sin();


/* FUNCTION TO COMPUTE CONSTANT (M). */


    ret_val = *e0 * *phi - *e1 * sin(two * *phi) + *e2 * sin(four * *phi) - *
	    e3 * sin(six * *phi);

    return ret_val;
} /* mlfnz0_ */

/*                   MSFNZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

doublereal msfnz0_(eccent, sinphi, cosphi)
doublereal *eccent, *sinphi, *cosphi;
{
    /* Initialized data */

    static doublereal one = 1.;

    /* System generated locals */
    doublereal ret_val;

    /* Builtin functions */
    double sqrt();

    /* Local variables */
    static doublereal con;


/* FUNCTION TO COMPUTE CONSTANT (SMALL M). */


    con = *eccent * *sinphi;
    ret_val = *cosphi / sqrt(one - con * con);

    return ret_val;
} /* msfnz0_ */

/*                   PAKCZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

doublereal pakcz0_(pak)
doublereal *pak;
{
    /* Initialized data */

    static doublereal zero = 0.;
    static doublereal con1 = 1e4;
    static doublereal con2 = 100.;
    static doublereal con3 = 1e6;
    static doublereal con4 = 1e3;
    static char iblank[1+1] = " ";
    static char neg[1+1] = "-";

    /* System generated locals */
    doublereal ret_val;

    /* Builtin functions */
    double d_mod();

    /* Local variables */
    static integer degs;
    static char sgna[1];
    static doublereal secs;
    static integer mins;
    static doublereal con;


/* SUBROUTINE TO CONVERT 2 DIGIT PACKED DMS TO 3 DIGIT PACKED DMS ANGLE. 
*/

/* SGNA : SIGN OF ANGLE */
/* DEGS : DEGREES PORTION OF ANGLE */
/* MINS : MINUTES PORTION OF ANGLE */
/* SECS : SECONDS PORTION OF ANGLE */


    *sgna = iblank[0];
    if (*pak < zero) {
	*sgna = neg[0];
    }
    con = abs(*pak);
    degs = (integer) (con / con1);
    con = d_mod(&con, &con1);
    mins = (integer) (con / con2);
    secs = d_mod(&con, &con2);

    con = (doublereal) degs * con3 + (doublereal) mins * con4 + secs;
    if (*sgna == neg[0]) {
	con = -con;
    }
    ret_val = con;
    return ret_val;

} /* pakcz0_ */

/*                   PAKDZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

/* Subroutine */ int pakdz0_(pak, sgna, degs, mins, secs, sgna_len)
doublereal *pak;
char *sgna;
integer *degs, *mins;
real *secs;
ftnlen sgna_len;
{
    /* Initialized data */

    static doublereal zero = 0.;
    static doublereal con1 = 1e6;
    static doublereal con2 = 1e3;
    static char iblank[1+1] = " ";
    static char neg[1+1] = "-";

    /* Builtin functions */
    double d_mod();

    /* Local variables */
    static doublereal con;


/* SUBROUTINE TO CONVERT PACKED DMS TO UNPACKED DMS ANGLE. */

/* SGNA : SIGN OF ANGLE */
/* DEGS : DEGREES PORTION OF ANGLE */
/* MINS : MINUTES PORTION OF ANGLE */
/* SECS : SECONDS PORTION OF ANGLE */


    *sgna = iblank[0];
    if (*pak < zero) {
	*sgna = neg[0];
    }
    con = abs(*pak);
    *degs = (integer) (con / con1);
    con = d_mod(&con, &con1);
    *mins = (integer) (con / con2);
    *secs = d_mod(&con, &con2);
    return 0;

} /* pakdz0_ */

/*                   PAKRZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

doublereal pakrz0_(ang)
doublereal *ang;
{
    /* Initialized data */

    static doublereal secrad = 4.848136811095359e-6;

    /* System generated locals */
    doublereal ret_val;

    /* Local variables */
    extern doublereal paksz0_();
    static doublereal sec;


/* FUNCTION TO CONVERT DMS PACKED ANGLE INTO RADIANS. */


/* CONVERT ANGLE TO SECONDS OF ARC */

    sec = paksz0_(ang);

/* CONVERT ANGLE TO RADIANS. */

    ret_val = sec * secrad;

    return ret_val;
} /* pakrz0_ */

/*                   PAKSZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

doublereal paksz0_(ang)
doublereal *ang;
{
    /* Initialized data */

    static doublereal code[2] = { 1e6,1e3 };
    static doublereal zero = 0.;
    static doublereal one = 1.;
    static doublereal c1 = 3600.;
    static doublereal c2 = 60.;

    /* System generated locals */
    doublereal ret_val;

    /* Local variables */
    static integer i;
    static doublereal factor, deg, sec, min__, tmp;


/* FUNCTION TO CONVERT DMS PACKED ANGLE INTO SECONDS OF ARC. */


/* SEPARATE DEGREE FIELD. */

    factor = one;
    if (*ang < zero) {
	factor = -one;
    }
    sec = abs(*ang);
    tmp = code[0];
    i = (integer) (sec / tmp);
    if (i > 360) {
	goto L20;
    }
    deg = (doublereal) i;

/* SEPARATE MINUTES FIELD. */

    sec -= deg * tmp;
    tmp = code[1];
    i = (integer) (sec / tmp);
    if (i > 60) {
	goto L20;
    }
    min__ = (doublereal) i;

/* SEPARATE SECONDS FIELD. */

    sec -= min__ * tmp;
    if (sec > c2) {
	goto L20;
    }
    sec = factor * (deg * c1 + min__ * c2 + sec);
    goto L40;

/* ERROR DETECTED IN DMS FORM. */

L20:

L40:
    ret_val = sec;

    return ret_val;
} /* paksz0_ */

/*                   PHI1Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

doublereal phi1z0_(eccent, qs)
doublereal *eccent, *qs;
{
    /* Initialized data */

    static doublereal half = .5;
    static doublereal one = 1.;
    static doublereal epsln = 1e-7;
    static doublereal tol = 1e-10;
    static integer nit = 15;

    /* System generated locals */
    integer i__1;
    doublereal ret_val, d__1;

    /* Builtin functions */
    double sin(), cos(), log();

    /* Local variables */
    static doublereal dphi, cospi, sinpi;
    extern doublereal asinz0_();
    static integer ii;
    static doublereal eccnts, com, con, phi;


/* FUNCTION TO COMPUTE LATITUDE ANGLE (PHI-1). */


    d__1 = half * *qs;
    ret_val = asinz0_(&d__1);
    if (*eccent < epsln) {
	return ret_val;
    }

    eccnts = *eccent * *eccent;
    phi = ret_val;
    i__1 = nit;
    for (ii = 1; ii <= i__1; ++ii) {
	sinpi = sin(phi);
	cospi = cos(phi);
	con = *eccent * sinpi;
	com = one - con * con;
	dphi = half * com * com / cospi * (*qs / (one - eccnts) - sinpi / com 
		+ half / *eccent * log((one - con) / (one + con)));
	phi += dphi;
	if (abs(dphi) > tol) {
	    goto L20;
	}
	ret_val = phi;
	return ret_val;
L20:
	;
    }

    errmz0_1.ierror = 1;
    return ret_val;

} /* phi1z0_ */

/*                   PHI2Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

doublereal phi2z0_(eccent, ts)
doublereal *eccent, *ts;
{
    /* Initialized data */

    static doublereal half = .5;
    static doublereal one = 1.;
    static doublereal two = 2.;
    static doublereal tol = 1e-10;
    static integer nit = 15;
    static doublereal halfpi = 1.5707963267948966;

    /* System generated locals */
    integer i__1;
    doublereal ret_val, d__1;

    /* Builtin functions */
    double atan(), sin(), pow_dd();

    /* Local variables */
    static doublereal dphi, sinpi;
    static integer ii;
    static doublereal eccnth, con, phi;


/* FUNCTION TO COMPUTE LATITUDE ANGLE (PHI-2). */


    eccnth = half * *eccent;
    phi = halfpi - two * atan(*ts);
    i__1 = nit;
    for (ii = 1; ii <= i__1; ++ii) {
	sinpi = sin(phi);
	con = *eccent * sinpi;
	d__1 = (one - con) / (one + con);
	dphi = halfpi - two * atan(*ts * pow_dd(&d__1, &eccnth)) - phi;
	phi += dphi;
	if (abs(dphi) > tol) {
	    goto L20;
	}
	ret_val = phi;
	return ret_val;
L20:
	;
    }

    errmz0_1.ierror = 2;
    ret_val = 0.;
    return ret_val;

} /* phi2z0_ */

/*                   PHI3Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

doublereal phi3z0_(ml, e0, e1, e2, e3)
doublereal *ml, *e0, *e1, *e2, *e3;
{
    /* Initialized data */

    static doublereal two = 2.;
    static doublereal four = 4.;
    static doublereal six = 6.;
    static doublereal tol = 1e-10;
    static integer nit = 15;

    /* System generated locals */
    integer i__1;
    doublereal ret_val;

    /* Builtin functions */
    double sin();

    /* Local variables */
    static doublereal dphi;
    static integer ii;
    static doublereal phi;


/* FUNCTION TO COMPUTE LATITUDE ANGLE (PHI-3). */


    phi = *ml;
    i__1 = nit;
    for (ii = 1; ii <= i__1; ++ii) {
	dphi = (*ml + *e1 * sin(two * phi) - *e2 * sin(four * phi) + *e3 * 
		sin(six * phi)) / *e0 - phi;
	phi += dphi;
	if (abs(dphi) > tol) {
	    goto L20;
	}
	ret_val = phi;
	return ret_val;
L20:
	;
    }

    errmz0_1.ierror = 3;
    ret_val = 0.;
    return ret_val;

} /* phi3z0_ */

/*                   PHI4Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

/* Subroutine */ int phi4z0_(eccnts, e0, e1, e2, e3, a, b, c, phi)
doublereal *eccnts, *e0, *e1, *e2, *e3, *a, *b, *c, *phi;
{
    /* Initialized data */

    static doublereal one = 1.;
    static doublereal two = 2.;
    static doublereal four = 4.;
    static doublereal six = 6.;
    static doublereal tol = 1e-10;
    static integer nit = 15;

    /* System generated locals */
    integer i__1;

    /* Builtin functions */
    double sin(), tan(), sqrt(), cos();

    /* Local variables */
    static doublereal dphi, sin2ph;
    static integer ii;
    static doublereal ml, tanphi, sinphi, mlp, con1, con2, con3;


/* FUNCTION TO COMPUTE LATITUDE ANGLE (PHI-4). */


    *phi = *a;
    i__1 = nit;
    for (ii = 1; ii <= i__1; ++ii) {
	sinphi = sin(*phi);
	tanphi = tan(*phi);
	*c = tanphi * sqrt(one - *eccnts * sinphi * sinphi);
	sin2ph = sin(two * *phi);
	ml = *e0 * *phi - *e1 * sin2ph + *e2 * sin(four * *phi) - *e3 * sin(
		six * *phi);
	mlp = *e0 - two * *e1 * cos(two * *phi) + four * *e2 * cos(four * *
		phi) - six * *e3 * cos(six * *phi);
	con1 = two * ml + *c * (ml * ml + *b) - two * *a * (*c * ml + one);
	con2 = *eccnts * sin2ph * (ml * ml + *b - two * *a * ml) / (two * *c);

	con3 = two * (*a - ml) * (*c * mlp - two / sin2ph) - two * mlp;
	dphi = con1 / (con2 + con3);
	*phi += dphi;
	if (abs(dphi) > tol) {
	    goto L20;
	}
	return 0;
L20:
	;
    }

    errmz0_1.ierror = 4;
    return 0;

} /* phi4z0_ */

/*                   PJ01Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */
/*                              *  U T M  * */
/* ********************************************************************** */

/* Subroutine */ int pj01z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    static doublereal geog[2], proj[2];
    extern /* Subroutine */ int pj09z0_();




/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */
    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[0] != 0) {
	    goto L140;
	}
	errmz0_1.ierror = 13;
	return 0;
L140:
	pj09z0_(geog, proj, &cs__0);
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[0] != 0) {
	    goto L160;
	}
	errmz0_1.ierror = 14;
	return 0;
L160:
	pj09z0_(proj, geog, &cs__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj01z0_ */

/*                   PJ02Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */
/*                           *  STATE PLANE  * */
/* ********************************************************************** */

/* Subroutine */ int pj02z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    static doublereal geog[2], proj[2];
    extern /* Subroutine */ int pj20z0_(), pj04z0_(), pj07z0_(), pj09z0_();




/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */
    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[1] == 0) {
	    errmz0_1.ierror = 23;
	    return 0;
	}

/*     TRANSVERSE MERCATOR PROJECTION */

	if (pj02_1.itype == 1) {
	    pj09z0_(geog, proj, &cs__0);
	}

/*     LAMBERT CONFORMAL PROJECTION */

	if (pj02_1.itype == 2) {
	    pj04z0_(geog, proj, &cs__0);
	}

/*     POLYCONIC PROJECTION */

	if (pj02_1.itype == 3) {
	    pj07z0_(geog, proj, &cs__0);
	}

/*     OBLIQUE MERCATOR PROJECTION */

	if (pj02_1.itype == 4) {
	    pj20z0_(geog, proj, &cs__0);
	}

	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[1] == 0) {
	    errmz0_1.ierror = 25;
	    return 0;
	}

/*     TRANSVERSE MERCATOR PROJECTION */

	if (pj02_1.itype == 1) {
	    pj09z0_(proj, geog, &cs__1);
	}

/*     LAMBERT CONFORMAL PROJECTION */

	if (pj02_1.itype == 2) {
	    pj04z0_(proj, geog, &cs__1);
	}

/*     POLYCONIC PROJECTION */

	if (pj02_1.itype == 3) {
	    pj07z0_(proj, geog, &cs__1);
	}

/*     OBLIQUE MERCATOR PROJECTION */

	if (pj02_1.itype == 4) {
	    pj20z0_(proj, geog, &cs__1);
	}

	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj02z0_ */

/*                   PJ03Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*                    *  ALBERS CONICAL EQUAL AREA  * */
/* ********************************************************************** */

/* Subroutine */ int pj03z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal tol = 1e-7;
    static doublereal halfpi = 1.5707963267948966;
    static doublereal zero = 0.;
    static doublereal half = .5;
    static doublereal one = 1.;

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double sin(), cos(), sqrt(), d_sign(), atan2(), log();

    /* Local variables */
    static doublereal geog[2], proj[2], x, y, theta;
    extern doublereal phi1z0_(), adjlz0_(), qsfnz0_();
    static doublereal rh, qs, cosphi, sinphi, con;


/* **** PARAMETERS **** A,E,ES,LAT1,LAT2,LON0,LAT0,X0,Y0,NS,C,RH0 ******* 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[2] != 0) {
	    goto L220;
	}
	errmz0_1.ierror = 33;
	return 0;
L220:
	sinphi = sin(geog[1]);
	cosphi = cos(geog[1]);
	qs = qsfnz0_(&pj03_1.e, &sinphi, &cosphi);
	rh = pj03_1.a * sqrt(pj03_1.c - pj03_1.ns * qs) / pj03_1.ns;
	d__1 = geog[0] - pj03_1.lon0;
	theta = pj03_1.ns * adjlz0_(&d__1);
	proj[0] = pj03_1.x0 + rh * sin(theta);
	proj[1] = pj03_1.y0 + pj03_1.rh0 - rh * cos(theta);
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[2] != 0) {
	    goto L240;
	}
	errmz0_1.ierror = 34;
	return 0;
L240:
	x = proj[0] - pj03_1.x0;
	y = pj03_1.rh0 - proj[1] + pj03_1.y0;
	d__1 = sqrt(x * x + y * y);
	rh = d_sign(&d__1, &pj03_1.ns);
	theta = zero;
	con = d_sign(&one, &pj03_1.ns);
	if (rh != zero) {
	    theta = atan2(con * x, con * y);
	}
	con = rh * pj03_1.ns / pj03_1.a;
	qs = (pj03_1.c - con * con) / pj03_1.ns;
	if (pj03_1.e < tol) {
	    goto L260;
	}
	con = one - half * (one - pj03_1.es) * log((one - pj03_1.e) / (one + 
		pj03_1.e)) / pj03_1.e;
	if (abs(con) - abs(qs) > tol) {
	    goto L260;
	}
	geog[1] = d_sign(&halfpi, &qs);
	goto L280;
L260:
	geog[1] = phi1z0_(&pj03_1.e, &qs);
	if (errmz0_1.ierror == 0) {
	    goto L280;
	}
	errmz0_1.ierror = 35;
	return 0;
L280:
	d__1 = theta / pj03_1.ns + pj03_1.lon0;
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj03z0_ */

/*                   PJ04Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*                     *  LAMBERT CONFORMAL CONIC  * */
/* ********************************************************************** */

/* Subroutine */ int pj04z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal halfpi = 1.5707963267948966;
    static doublereal epsln = 1e-10;
    static doublereal zero = 0.;
    static doublereal one = 1.;

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double sin(), pow_dd(), cos(), sqrt(), d_sign(), atan2();

    /* Local variables */
    static doublereal geog[2], proj[2], x, y, theta;
    extern doublereal phi2z0_(), adjlz0_(), tsfnz0_();
    static doublereal rh, ts, sinphi, con;


/* **** PARAMETERS **** A,E,ES,LAT1,LAT2,LON0,LAT0,X0,Y0,NS,F,RH0 ******* 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[3] != 0) {
	    goto L200;
	}
	errmz0_1.ierror = 43;
	return 0;
L200:
	con = (d__1 = abs(geog[1]) - halfpi, abs(d__1));
	if (con > epsln) {
	    goto L220;
	}
	con = geog[1] * pj04_1.ns;
	if (con > zero) {
	    goto L210;
	}
	errmz0_1.ierror = 44;
	return 0;
L210:
	rh = zero;
	goto L230;
L220:
	sinphi = sin(geog[1]);
	ts = tsfnz0_(&pj04_1.e, &geog[1], &sinphi);
	rh = pj04_1.a * pj04_1.f * pow_dd(&ts, &pj04_1.ns);
L230:
	d__1 = geog[0] - pj04_1.lon0;
	theta = pj04_1.ns * adjlz0_(&d__1);
	proj[0] = pj04_1.x0 + rh * sin(theta);
	proj[1] = pj04_1.y0 + pj04_1.rh0 - rh * cos(theta);
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[3] != 0) {
	    goto L240;
	}
	errmz0_1.ierror = 45;
	return 0;
L240:
	x = proj[0] - pj04_1.x0;
	y = pj04_1.rh0 - proj[1] + pj04_1.y0;
	d__1 = sqrt(x * x + y * y);
	rh = d_sign(&d__1, &pj04_1.ns);
	theta = zero;
	con = d_sign(&one, &pj04_1.ns);
	if (rh != zero) {
	    theta = atan2(con * x, con * y);
	}
	if (rh != zero || pj04_1.ns > zero) {
	    goto L250;
	}
	geog[1] = -halfpi;
	goto L260;
L250:
	con = one / pj04_1.ns;
	d__1 = rh / (pj04_1.a * pj04_1.f);
	ts = pow_dd(&d__1, &con);
	geog[1] = phi2z0_(&pj04_1.e, &ts);
	if (errmz0_1.ierror == 0) {
	    goto L260;
	}
	errmz0_1.ierror = 46;
	return 0;
L260:
	d__1 = theta / pj04_1.ns + pj04_1.lon0;
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj04z0_ */

/*                   PJ05Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*                            *  MERCATOR  * */
/* ********************************************************************** */

/* Subroutine */ int pj05z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal halfpi = 1.5707963267948966;
    static doublereal epsln = 1e-10;

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double sin(), log(), exp();

    /* Local variables */
    static doublereal geog[2], proj[2], x, y;
    extern doublereal phi2z0_(), adjlz0_(), tsfnz0_();
    static doublereal ts, sinphi;


/* **** PARAMETERS **** A,E,ES,LON0,X0,Y0,NS,F,RH0,LAT1,M1 ************** 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[4] != 0) {
	    goto L220;
	}
	errmz0_1.ierror = 52;
	return 0;
L220:
	if ((d__1 = abs(geog[1]) - halfpi, abs(d__1)) > epsln) {
	    goto L240;
	}
	errmz0_1.ierror = 53;
	return 0;
L240:
	sinphi = sin(geog[1]);
	ts = tsfnz0_(&pj05_1.e, &geog[1], &sinphi);
	d__1 = geog[0] - pj05_1.lon0;
	proj[0] = pj05_1.x0 + pj05_1.a * pj05_1.m1 * adjlz0_(&d__1);
	proj[1] = pj05_1.y0 - pj05_1.a * pj05_1.m1 * log(ts);
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[4] != 0) {
	    goto L260;
	}
	errmz0_1.ierror = 54;
	return 0;
L260:
	x = proj[0] - pj05_1.x0;
	y = proj[1] - pj05_1.y0;
	ts = exp(-y / (pj05_1.a * pj05_1.m1));
	geog[1] = phi2z0_(&pj05_1.e, &ts);
	if (errmz0_1.ierror == 0) {
	    goto L280;
	}
	errmz0_1.ierror = 55;
	return 0;
L280:
	d__1 = pj05_1.lon0 + x / (pj05_1.a * pj05_1.m1);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj05z0_ */

/*                   PJ06Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*                       *  POLAR STEREOGRAPHIC  * */
/* ********************************************************************** */

/* Subroutine */ int pj06z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal zero = 0.;
    static doublereal two = 2.;

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double sin(), cos(), sqrt(), atan2();

    /* Local variables */
    static doublereal geog[2], proj[2], x, y;
    extern doublereal phi2z0_(), adjlz0_(), tsfnz0_();
    static doublereal rh, ts, sinphi, con1, con2;


/* **** PARAMETERS **** A,E,ES,LON0,LATC,X0,Y0,E4,MCS,TCS,FAC,IND ******* 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[5] != 0) {
	    goto L220;
	}
	errmz0_1.ierror = 62;
	return 0;
L220:
	d__1 = geog[0] - pj06_1.lon0;
	con1 = pj06_1.fac * adjlz0_(&d__1);
	con2 = pj06_1.fac * geog[1];
	sinphi = sin(con2);
	ts = tsfnz0_(&pj06_1.e, &con2, &sinphi);
	if (pj06_1.ind == 0) {
	    goto L240;
	}
	rh = pj06_1.a * pj06_1.mcs * ts / pj06_1.tcs;
	goto L260;
L240:
	rh = two * pj06_1.a * ts / pj06_1.e4;
L260:
	proj[0] = pj06_1.x0 + pj06_1.fac * rh * sin(con1);
	proj[1] = pj06_1.y0 - pj06_1.fac * rh * cos(con1);
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[5] != 0) {
	    goto L320;
	}
	errmz0_1.ierror = 63;
	return 0;
L320:
	x = pj06_1.fac * (proj[0] - pj06_1.x0);
	y = pj06_1.fac * (proj[1] - pj06_1.y0);
	rh = sqrt(x * x + y * y);
	if (pj06_1.ind == 0) {
	    goto L340;
	}
	ts = rh * pj06_1.tcs / (pj06_1.a * pj06_1.mcs);
	goto L360;
L340:
	ts = rh * pj06_1.e4 / (two * pj06_1.a);
L360:
	geog[1] = pj06_1.fac * phi2z0_(&pj06_1.e, &ts);
	if (errmz0_1.ierror == 0) {
	    goto L380;
	}
	errmz0_1.ierror = 64;
	return 0;
L380:
	if (rh != zero) {
	    goto L400;
	}
	geog[0] = pj06_1.fac * pj06_1.lon0;
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L400:
	d__1 = pj06_1.fac * atan2(x, -y) + pj06_1.lon0;
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj06z0_ */

/*                   PJ07Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*                            *  POLYCONIC  * */
/* ********************************************************************** */

/* Subroutine */ int pj07z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal tol = 1e-7;
    static doublereal zero = 0.;
    static doublereal one = 1.;

    /* System generated locals */
    doublereal d__1, d__2;

    /* Builtin functions */
    double sin(), cos();

    /* Local variables */
    static doublereal geog[2], proj[2], b, c, x, y;
    extern /* Subroutine */ int phi4z0_();
    extern doublereal adjlz0_(), asinz0_(), mlfnz0_(), msfnz0_();
    static doublereal al, ml, ms, cosphi, sinphi, con;


/* **** PARAMETERS **** A,E,ES,LON0,LAT0,X0,Y0,E0,E1,E2,ML0 ************* 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[6] != 0) {
	    goto L220;
	}
	errmz0_1.ierror = 72;
	return 0;
L220:
	d__1 = geog[0] - pj07_1.lon0;
	con = adjlz0_(&d__1);
	if (abs(geog[1]) > tol) {
	    goto L240;
	}
	proj[0] = pj07_1.x0 + pj07_1.a * con;
	proj[1] = pj07_1.y0 - pj07_1.a * pj07_1.ml0;
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
L240:
	sinphi = sin(geog[1]);
	cosphi = cos(geog[1]);
	ml = mlfnz0_(&pj07_1.e0, &pj07_1.e1, &pj07_1.e2, &pj07_1.e3, &geog[1])
		;
	ms = msfnz0_(&pj07_1.e, &sinphi, &cosphi);
	con *= sinphi;
	proj[0] = pj07_1.x0 + pj07_1.a * ms * sin(con) / sinphi;
	proj[1] = pj07_1.y0 + pj07_1.a * (ml - pj07_1.ml0 + ms * (one - cos(
		con)) / sinphi);
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[6] != 0) {
	    goto L320;
	}
	errmz0_1.ierror = 73;
	return 0;
L320:
	x = proj[0] - pj07_1.x0;
	y = proj[1] - pj07_1.y0;
	al = pj07_1.ml0 + y / pj07_1.a;
	if (abs(al) > tol) {
	    goto L340;
	}
	geog[0] = x / pj07_1.a + pj07_1.lon0;
	geog[1] = zero;
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L340:
/* Computing 2nd power */
	d__1 = x / pj07_1.a;
	b = al * al + d__1 * d__1;
	phi4z0_(&pj07_1.es, &pj07_1.e0, &pj07_1.e1, &pj07_1.e2, &pj07_1.e3, &
		al, &b, &c, &geog[1]);
	if (errmz0_1.ierror == 0) {
	    goto L360;
	}
	errmz0_1.ierror = 74;
	return 0;
L360:
	d__2 = x * c / pj07_1.a;
	d__1 = asinz0_(&d__2) / sin(geog[1]) + pj07_1.lon0;
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj07z0_ */

/*                   PJ08Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*                        *  EQUIDISTANT CONIC  * */
/* ********************************************************************** */

/* Subroutine */ int pj08z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal zero = 0.;
    static doublereal one = 1.;

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double sin(), cos(), sqrt(), d_sign(), atan2();

    /* Local variables */
    static doublereal geog[2], proj[2], x, y, theta;
    extern doublereal phi3z0_(), adjlz0_(), mlfnz0_();
    static doublereal ml, rh, con;


/* ** PARAMETERS * A,E,ES,LAT1,LAT2,LON0,LAT0,X0,Y0,E0,E1,E2,E3,NS,GL,RH0 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[7] != 0) {
	    goto L300;
	}
	errmz0_1.ierror = 83;
	return 0;
L300:
	ml = mlfnz0_(&pj08_1.e0, &pj08_1.e1, &pj08_1.e2, &pj08_1.e3, &geog[1])
		;
	rh = pj08_1.a * (pj08_1.gl - ml);
	d__1 = geog[0] - pj08_1.lon0;
	theta = pj08_1.ns * adjlz0_(&d__1);
	proj[0] = pj08_1.x0 + rh * sin(theta);
	proj[1] = pj08_1.y0 + pj08_1.rh0 - rh * cos(theta);
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[7] != 0) {
	    goto L320;
	}
	errmz0_1.ierror = 84;
	return 0;
L320:
	x = proj[0] - pj08_1.x0;
	y = pj08_1.rh0 - proj[1] + pj08_1.y0;
	d__1 = sqrt(x * x + y * y);
	rh = d_sign(&d__1, &pj08_1.ns);
	theta = zero;
	con = d_sign(&one, &pj08_1.ns);
	if (rh != zero) {
	    theta = atan2(con * x, con * y);
	}
	ml = pj08_1.gl - rh / pj08_1.a;
	geog[1] = phi3z0_(&ml, &pj08_1.e0, &pj08_1.e1, &pj08_1.e2, &pj08_1.e3)
		;
	if (errmz0_1.ierror == 0) {
	    goto L340;
	}
	errmz0_1.ierror = 85;
	return 0;
L340:
	d__1 = pj08_1.lon0 + theta / pj08_1.ns;
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj08z0_ */

/*                   PJ09Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*                       *  TRANSVERSE MERCATOR  * */
/* ********************************************************************** */

/* Subroutine */ int pj09z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal zero = 0.;
    static doublereal half = .5;
    static doublereal one = 1.;
    static doublereal two = 2.;
    static doublereal three = 3.;
    static doublereal four = 4.;
    static doublereal five = 5.;
    static doublereal six = 6.;
    static doublereal eight = 8.;
    static doublereal nine = 9.;
    static doublereal halfpi = 1.5707963267948966;
    static doublereal ten = 10.;
    static doublereal epsln = 1e-10;
    static integer nit = 6;

    /* System generated locals */
    integer i__1;
    doublereal d__1;

    /* Builtin functions */
    double cos(), sin(), log(), sqrt(), acos(), tan(), exp(), atan2(), d_sign(
	    );

    /* Local variables */
    static doublereal geog[2], dphi, dlon, temp, proj[2], b, c, d, f, g, h;
    static integer i;
    static doublereal n, r, t, x, y;
    extern doublereal adjlz0_(), asinz0_(), mlfnz0_();
    static doublereal al, cs, ds, ml, tq, ts, tanphi, cosphi, sinphi, con, 
	    lat, als, phi;


/* **** PARAMETERS ** A,E,ES,KS0,LON0,LAT0,X0,Y0,E0,E1,E2,E3,ESP,ML0,IND 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[8] != 0) {
	    goto L220;
	}
	errmz0_1.ierror = 92;
	return 0;
L220:
	d__1 = geog[0] - pj09_1.lon0;
	dlon = adjlz0_(&d__1);
	lat = geog[1];
	if (pj09_1.ind == 0) {
	    goto L240;
	}
	cosphi = cos(lat);
	b = cosphi * sin(dlon);
	if ((d__1 = abs(b) - one, abs(d__1)) > epsln) {
	    goto L230;
	}
	errmz0_1.ierror = 93;
	return 0;
L230:
	proj[0] = half * pj09_1.a * pj09_1.ks0 * log((one + b) / (one - b));
	con = acos(cosphi * cos(dlon) / sqrt(one - b * b));
	if (lat < zero) {
	    con = -con;
	}
	proj[1] = pj09_1.a * pj09_1.ks0 * (con - pj09_1.lat0);
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;

L240:
	sinphi = sin(lat);
	cosphi = cos(lat);
	al = cosphi * dlon;
	als = al * al;
	c = pj09_1.esp * cosphi * cosphi;
	tq = tan(lat);
	t = tq * tq;
	n = pj09_1.a / sqrt(one - pj09_1.es * sinphi * sinphi);
	ml = pj09_1.a * mlfnz0_(&pj09_1.e0, &pj09_1.e1, &pj09_1.e2, &
		pj09_1.e3, &lat);
	proj[0] = pj09_1.ks0 * n * al * (one + als / six * (one - t + c + als 
		/ 20. * (five - t * 18. + t * t + c * 72. - pj09_1.esp * 58.))
		) + pj09_1.x0;
	proj[1] = pj09_1.ks0 * (ml - pj09_1.ml0 + n * tq * (als * (half + als 
		/ 24. * (five - t + nine * c + four * c * c + als / 30. * (
		61. - t * 58. + t * t + c * 600. - pj09_1.esp * 330.))))) + 
		pj09_1.y0;
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[8] != 0) {
	    goto L320;
	}
	errmz0_1.ierror = 94;
	return 0;
L320:
	x = proj[0] - pj09_1.x0;
	y = proj[1] - pj09_1.y0;
	if (pj09_1.ind == 0) {
	    goto L340;
	}
	f = exp(x / (pj09_1.a * pj09_1.ks0));
	g = half * (f - one / f);
	temp = pj09_1.lat0 + y / (pj09_1.a * pj09_1.ks0);
	h = cos(temp);
	con = sqrt((one - h * h) / (one + g * g));
	geog[1] = asinz0_(&con);
	if (temp < zero) {
	    geog[1] = -geog[1];
	}
	if (g != zero || h != zero) {
	    goto L330;
	}
	geog[0] = pj09_1.lon0;
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L330:
	d__1 = atan2(g, h) + pj09_1.lon0;
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;

L340:
	con = (pj09_1.ml0 + y / pj09_1.ks0) / pj09_1.a;
	phi = con;
	i__1 = nit;
	for (i = 1; i <= i__1; ++i) {
	    dphi = (con + pj09_1.e1 * sin(two * phi) - pj09_1.e2 * sin(four * 
		    phi) + pj09_1.e3 * sin(six * phi)) / pj09_1.e0 - phi;
	    phi += dphi;
	    if (abs(dphi) <= epsln) {
		goto L380;
	    }
/* L360: */
	}
	errmz0_1.ierror = 95;
	return 0;
L380:
	if (abs(phi) < halfpi) {
	    goto L400;
	}
	geog[1] = d_sign(&halfpi, &y);
	geog[0] = pj09_1.lon0;
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L400:
	sinphi = sin(phi);
	cosphi = cos(phi);
	tanphi = tan(phi);
	c = pj09_1.esp * cosphi * cosphi;
	cs = c * c;
	t = tanphi * tanphi;
	ts = t * t;
	con = one - pj09_1.es * sinphi * sinphi;
	n = pj09_1.a / sqrt(con);
	r = n * (one - pj09_1.es) / con;
	d = x / (n * pj09_1.ks0);
	ds = d * d;
	geog[1] = phi - n * tanphi * ds / r * (half - ds / 24. * (five + 
		three * t + ten * c - four * cs - nine * pj09_1.esp - ds / 
		30. * (t * 90. + 61. + c * 298. + ts * 45. - pj09_1.esp * 
		252. - three * cs)));
	d__1 = pj09_1.lon0 + d * (one - ds / six * (one + two * t + c - ds / 
		20. * (five - two * c + t * 28. - three * cs + eight * 
		pj09_1.esp + ts * 24.))) / cosphi;
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj09z0_ */

/*                   PJ10Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*                          *  STEREOGRAPHIC  * */
/* ********************************************************************** */

/* Subroutine */ int pj10z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal halfpi = 1.5707963267948966;
    static doublereal epsln = 1e-10;
    static doublereal zero = 0.;
    static doublereal one = 1.;
    static doublereal two = 2.;

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double sin(), cos(), sqrt(), atan(), atan2();

    /* Local variables */
    static doublereal geog[2], proj[2], cosz, sinz, g, x, y, z;
    extern doublereal adjlz0_(), asinz0_();
    static doublereal rh, cosphi, sinphi, coslon, con, lon, ksp;


/* **** PARAMETERS **** A,LON0,LAT0,X0,Y0,SINPH0,COSPH0 ***************** 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[9] != 0) {
	    goto L120;
	}
	errmz0_1.ierror = 102;
	return 0;
L120:
	d__1 = geog[0] - pj10_1.lon0;
	lon = adjlz0_(&d__1);
	sinphi = sin(geog[1]);
	cosphi = cos(geog[1]);
	coslon = cos(lon);
	g = pj10_1.sinph0 * sinphi + pj10_1.cosph0 * cosphi * coslon;
	if ((d__1 = g + one, abs(d__1)) > epsln) {
	    goto L140;
	}
	errmz0_1.ierror = 103;
	return 0;
L140:
	ksp = two / (one + g);
	proj[0] = pj10_1.x0 + pj10_1.a * ksp * cosphi * sin(lon);
	proj[1] = pj10_1.y0 + pj10_1.a * ksp * (pj10_1.cosph0 * sinphi - 
		pj10_1.sinph0 * cosphi * coslon);
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[9] != 0) {
	    goto L220;
	}
	errmz0_1.ierror = 104;
	return 0;
L220:
	x = proj[0] - pj10_1.x0;
	y = proj[1] - pj10_1.y0;
	rh = sqrt(x * x + y * y);
	z = two * atan(rh / (two * pj10_1.a));
	sinz = sin(z);
	cosz = cos(z);
	geog[0] = pj10_1.lon0;
	if (abs(rh) > epsln) {
	    goto L240;
	}
	geog[1] = pj10_1.lat0;
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L240:
	d__1 = cosz * pj10_1.sinph0 + y * sinz * pj10_1.cosph0 / rh;
	geog[1] = asinz0_(&d__1);
	con = abs(pj10_1.lat0) - halfpi;
	if (abs(con) > epsln) {
	    goto L260;
	}
	if (pj10_1.lat0 < zero) {
	    goto L250;
	}
	d__1 = pj10_1.lon0 + atan2(x, -y);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L250:
	d__1 = pj10_1.lon0 - atan2(-x, y);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L260:
	con = cosz - pj10_1.sinph0 * sin(geog[1]);
	if (abs(con) < epsln && abs(x) < epsln) {
	    return 0;
	}
	d__1 = pj10_1.lon0 + atan2(x * sinz * pj10_1.cosph0, con * rh);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj10z0_ */

/*                   PJ11Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*                   *  LAMBERT AZIMUTHAL EQUAL-AREA  * */
/* ********************************************************************** */

/* Subroutine */ int pj11z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal halfpi = 1.5707963267948966;
    static doublereal epsln = 1e-10;
    static doublereal zero = 0.;
    static doublereal one = 1.;
    static doublereal two = 2.;

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double sin(), cos(), sqrt(), atan2();

    /* Local variables */
    static doublereal geog[2], proj[2], cosz, sinz, g, x, y, z;
    extern doublereal adjlz0_(), asinz0_();
    static doublereal rh, cosphi, sinphi, coslon, con, lon, ksp;


/* **** PARAMETERS **** A,LON0,LAT0,X0,Y0,SINPH0,COSPH0 ***************** 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[10] != 0) {
	    goto L120;
	}
	errmz0_1.ierror = 112;
	return 0;
L120:
	d__1 = geog[0] - pj11_1.lon0;
	lon = adjlz0_(&d__1);
	sinphi = sin(geog[1]);
	cosphi = cos(geog[1]);
	coslon = cos(lon);
	g = pj11_1.sinph0 * sinphi + pj11_1.cosph0 * cosphi * coslon;
	if (g != -one) {
	    goto L140;
	}
	con = two * pj11_1.a;
	errmz0_1.ierror = 113;
	return 0;
L140:
	ksp = sqrt(two / (one + g));
	proj[0] = pj11_1.x0 + pj11_1.a * ksp * cosphi * sin(lon);
	proj[1] = pj11_1.y0 + pj11_1.a * ksp * (pj11_1.cosph0 * sinphi - 
		pj11_1.sinph0 * cosphi * coslon);
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[10] != 0) {
	    goto L220;
	}
	errmz0_1.ierror = 114;
	return 0;
L220:
	x = proj[0] - pj11_1.x0;
	y = proj[1] - pj11_1.y0;
	rh = sqrt(x * x + y * y);
	con = rh / (two * pj11_1.a);
	if (con <= one) {
	    goto L230;
	}
	errmz0_1.ierror = 115;
	return 0;
L230:
	z = two * asinz0_(&con);
	sinz = sin(z);
	cosz = cos(z);
	geog[0] = pj11_1.lon0;
	if (abs(rh) > epsln) {
	    goto L240;
	}
	geog[1] = pj11_1.lat0;
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L240:
	d__1 = cosz * pj11_1.sinph0 + y * sinz * pj11_1.cosph0 / rh;
	geog[1] = asinz0_(&d__1);
	con = abs(pj11_1.lat0) - halfpi;
	if (abs(con) > epsln) {
	    goto L260;
	}
	if (pj11_1.lat0 < zero) {
	    goto L250;
	}
	d__1 = pj11_1.lon0 + atan2(x, -y);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L250:
	d__1 = pj11_1.lon0 - atan2(-x, y);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L260:
	con = cosz - pj11_1.sinph0 * sin(geog[1]);
	if (con == zero) {
	    return 0;
	}
	d__1 = pj11_1.lon0 + atan2(x * sinz * pj11_1.cosph0, con * rh);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj11z0_ */

/*                   PJ12Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*                      *  AZIMUTHAL EQUIDISTANT  * */
/* ********************************************************************** */

/* Subroutine */ int pj12z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal halfpi = 1.5707963267948966;
    static doublereal epsln = 1e-10;
    static doublereal zero = 0.;
    static doublereal one = 1.;
    static doublereal two = 2.;

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double sin(), cos(), acos(), sqrt(), atan2();

    /* Local variables */
    static doublereal geog[2], proj[2], cosz, sinz, g, x, z, y;
    extern doublereal adjlz0_(), asinz0_();
    static doublereal rh, cosphi, sinphi, coslon, con, lon, ksp;


/* **** PARAMETERS **** A,LON0,LAT0,X0,Y0,SINPH0,COSPH0 ***************** 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[11] != 0) {
	    goto L120;
	}
	errmz0_1.ierror = 122;
	return 0;
L120:
	d__1 = geog[0] - pj12_1.lon0;
	lon = adjlz0_(&d__1);
	sinphi = sin(geog[1]);
	cosphi = cos(geog[1]);
	coslon = cos(lon);
	g = pj12_1.sinph0 * sinphi + pj12_1.cosph0 * cosphi * coslon;
	if ((d__1 = abs(g) - one, abs(d__1)) >= epsln) {
	    goto L140;
	}
	ksp = one;
	if (g >= zero) {
	    goto L160;
	}
	con = two * halfpi * pj12_1.a;
	errmz0_1.ierror = 123;
	return 0;
L140:
	z = acos(g);
	ksp = z / sin(z);
L160:
	proj[0] = pj12_1.x0 + pj12_1.a * ksp * cosphi * sin(lon);
	proj[1] = pj12_1.y0 + pj12_1.a * ksp * (pj12_1.cosph0 * sinphi - 
		pj12_1.sinph0 * cosphi * coslon);
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[11] != 0) {
	    goto L220;
	}
	errmz0_1.ierror = 124;
	return 0;
L220:
	x = proj[0] - pj12_1.x0;
	y = proj[1] - pj12_1.y0;
	rh = sqrt(x * x + y * y);
	if (rh <= two * halfpi * pj12_1.a) {
	    goto L230;
	}
	errmz0_1.ierror = 125;
	return 0;
L230:
	z = rh / pj12_1.a;
	sinz = sin(z);
	cosz = cos(z);
	geog[0] = pj12_1.lon0;
	if (abs(rh) > epsln) {
	    goto L240;
	}
	geog[1] = pj12_1.lat0;
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L240:
	d__1 = cosz * pj12_1.sinph0 + y * sinz * pj12_1.cosph0 / rh;
	geog[1] = asinz0_(&d__1);
	con = abs(pj12_1.lat0) - halfpi;
	if (abs(con) > epsln) {
	    goto L260;
	}
	if (pj12_1.lat0 < zero) {
	    goto L250;
	}
	d__1 = pj12_1.lon0 + atan2(x, -y);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L250:
	d__1 = pj12_1.lon0 - atan2(-x, y);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L260:
	con = cosz - pj12_1.sinph0 * sin(geog[1]);
	if (abs(con) < epsln && abs(x) < epsln) {
	    return 0;
	}
	d__1 = pj12_1.lon0 + atan2(x * sinz * pj12_1.cosph0, con * rh);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj12z0_ */

/*                   PJ13Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*                            *  GNOMONIC  * */
/* ********************************************************************** */

/* Subroutine */ int pj13z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal halfpi = 1.5707963267948966;
    static doublereal epsln = 1e-10;
    static doublereal zero = 0.;
    static doublereal one = 1.;

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double sin(), cos(), sqrt(), atan(), atan2();

    /* Local variables */
    static doublereal geog[2], proj[2], cosz, sinz, g, x, y, z;
    extern doublereal adjlz0_(), asinz0_();
    static doublereal rh, cosphi, sinphi, coslon, con, lon, ksp;


/* **** PARAMETERS **** A,LON0,LAT0,X0,Y0,SINPH0,COSPH0 ***************** 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[12] != 0) {
	    goto L120;
	}
	errmz0_1.ierror = 132;
	return 0;
L120:
	d__1 = geog[0] - pj13_1.lon0;
	lon = adjlz0_(&d__1);
	sinphi = sin(geog[1]);
	cosphi = cos(geog[1]);
	coslon = cos(lon);
	g = pj13_1.sinph0 * sinphi + pj13_1.cosph0 * cosphi * coslon;
	if (g > zero) {
	    goto L140;
	}
	errmz0_1.ierror = 133;
	return 0;
L140:
	ksp = one / g;
	proj[0] = pj13_1.x0 + pj13_1.a * ksp * cosphi * sin(lon);
	proj[1] = pj13_1.y0 + pj13_1.a * ksp * (pj13_1.cosph0 * sinphi - 
		pj13_1.sinph0 * cosphi * coslon);
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[12] != 0) {
	    goto L220;
	}
	errmz0_1.ierror = 134;
	return 0;
L220:
	x = proj[0] - pj13_1.x0;
	y = proj[1] - pj13_1.y0;
	rh = sqrt(x * x + y * y);
	z = atan(rh / pj13_1.a);
	sinz = sin(z);
	cosz = cos(z);
	geog[0] = pj13_1.lon0;
	if (abs(rh) > epsln) {
	    goto L240;
	}
	geog[1] = pj13_1.lat0;
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L240:
	d__1 = cosz * pj13_1.sinph0 + y * sinz * pj13_1.cosph0 / rh;
	geog[1] = asinz0_(&d__1);
	con = abs(pj13_1.lat0) - halfpi;
	if (abs(con) > epsln) {
	    goto L260;
	}
	if (pj13_1.lat0 < zero) {
	    goto L250;
	}
	d__1 = pj13_1.lon0 + atan2(x, -y);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L250:
	d__1 = pj13_1.lon0 - atan2(-x, y);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L260:
	con = cosz - pj13_1.sinph0 * sin(geog[1]);
	if (abs(con) < epsln && abs(x) < epsln) {
	    return 0;
	}
	d__1 = pj13_1.lon0 + atan2(x * sinz * pj13_1.cosph0, con * rh);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj13z0_ */

/*                   PJ14Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*                          *  ORTHOGRAPHIC  * */
/* ********************************************************************** */

/* Subroutine */ int pj14z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal halfpi = 1.5707963267948966;
    static doublereal epsln = 1e-10;
    static doublereal zero = 0.;
    static doublereal one = 1.;

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double sin(), cos(), sqrt(), atan2();

    /* Local variables */
    static doublereal geog[2], proj[2], cosz, sinz, g, x, y, z;
    extern doublereal adjlz0_(), asinz0_();
    static doublereal rh, cosphi, sinphi, coslon, con, lon, ksp;


/* **** PARAMETERS **** A,LON0,LAT0,X0,Y0,SINPH0,COSPH0 ***************** 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[13] != 0) {
	    goto L120;
	}
	errmz0_1.ierror = 142;
	return 0;
L120:
	d__1 = geog[0] - pj14_1.lon0;
	lon = adjlz0_(&d__1);
	sinphi = sin(geog[1]);
	cosphi = cos(geog[1]);
	coslon = cos(lon);
	g = pj14_1.sinph0 * sinphi + pj14_1.cosph0 * cosphi * coslon;
	ksp = one;
	if (g > zero || abs(g) <= epsln) {
	    goto L140;
	}
	errmz0_1.ierror = 143;
	return 0;
L140:
	proj[0] = pj14_1.x0 + pj14_1.a * ksp * cosphi * sin(lon);
	proj[1] = pj14_1.y0 + pj14_1.a * ksp * (pj14_1.cosph0 * sinphi - 
		pj14_1.sinph0 * cosphi * coslon);
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[13] != 0) {
	    goto L220;
	}
	errmz0_1.ierror = 144;
	return 0;
L220:
	x = proj[0] - pj14_1.x0;
	y = proj[1] - pj14_1.y0;
	rh = sqrt(x * x + y * y);
	if (rh <= pj14_1.a) {
	    goto L230;
	}
	errmz0_1.ierror = 145;
	return 0;
L230:
	d__1 = rh / pj14_1.a;
	z = asinz0_(&d__1);
	sinz = sin(z);
	cosz = cos(z);
	geog[0] = pj14_1.lon0;
	if (abs(rh) > epsln) {
	    goto L240;
	}
	geog[1] = pj14_1.lat0;
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L240:
	d__1 = cosz * pj14_1.sinph0 + y * sinz * pj14_1.cosph0 / rh;
	geog[1] = asinz0_(&d__1);
	con = abs(pj14_1.lat0) - halfpi;
	if (abs(con) > epsln) {
	    goto L260;
	}
	if (pj14_1.lat0 < zero) {
	    goto L250;
	}
	d__1 = pj14_1.lon0 + atan2(x, -y);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L250:
	d__1 = pj14_1.lon0 - atan2(-x, y);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L260:
	con = cosz - pj14_1.sinph0 * sin(geog[1]);
	if (abs(con) < epsln && abs(x) < epsln) {
	    return 0;
	}
	d__1 = pj14_1.lon0 + atan2(x * sinz * pj14_1.cosph0, con * rh);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj14z0_ */

/*                   PJ15Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*              *  GENERAL VERTICAL NEAR-SIDE PERSPECTIVE  * */
/* ********************************************************************** */

/* Subroutine */ int pj15z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal halfpi = 1.5707963267948966;
    static doublereal epsln = 1e-10;
    static doublereal zero = 0.;
    static doublereal one = 1.;

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double sin(), cos(), sqrt(), atan2();

    /* Local variables */
    static doublereal geog[2], proj[2], cosz, sinz, g, r, x, y, z;
    extern doublereal adjlz0_(), asinz0_();
    static doublereal rh, cosphi, sinphi, coslon, com, con, lon, ksp;


/* **** PARAMETERS **** A,P,LON0,LAT0,X0,Y0,SINPH0,COSPH0 *************** 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[14] != 0) {
	    goto L120;
	}
	errmz0_1.ierror = 152;
	return 0;
L120:
	d__1 = geog[0] - pj15_1.lon0;
	lon = adjlz0_(&d__1);
	sinphi = sin(geog[1]);
	cosphi = cos(geog[1]);
	coslon = cos(lon);
	g = pj15_1.sinph0 * sinphi + pj15_1.cosph0 * cosphi * coslon;
	if (g >= one / pj15_1.p) {
	    goto L140;
	}
	errmz0_1.ierror = 153;
	return 0;
L140:
	ksp = (pj15_1.p - one) / (pj15_1.p - g);
	proj[0] = pj15_1.x0 + pj15_1.a * ksp * cosphi * sin(lon);
	proj[1] = pj15_1.y0 + pj15_1.a * ksp * (pj15_1.cosph0 * sinphi - 
		pj15_1.sinph0 * cosphi * coslon);
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[14] != 0) {
	    goto L220;
	}
	errmz0_1.ierror = 154;
	return 0;
L220:
	x = proj[0] - pj15_1.x0;
	y = proj[1] - pj15_1.y0;
	rh = sqrt(x * x + y * y);
	r = rh / pj15_1.a;
	con = pj15_1.p - one;
	com = pj15_1.p + one;
	if (r <= sqrt(con / com)) {
	    goto L230;
	}
	errmz0_1.ierror = 155;
	return 0;
L230:
	sinz = (pj15_1.p - sqrt(one - r * r * com / con)) / (con / r + r / 
		con);
	z = asinz0_(&sinz);
	sinz = sin(z);
	cosz = cos(z);
	geog[0] = pj15_1.lon0;
	if (abs(rh) > epsln) {
	    goto L240;
	}
	geog[1] = pj15_1.lat0;
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L240:
	d__1 = cosz * pj15_1.sinph0 + y * sinz * pj15_1.cosph0 / rh;
	geog[1] = asinz0_(&d__1);
	con = abs(pj15_1.lat0) - halfpi;
	if (abs(con) > epsln) {
	    goto L260;
	}
	if (pj15_1.lat0 < zero) {
	    goto L250;
	}
	d__1 = pj15_1.lon0 + atan2(x, -y);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L250:
	d__1 = pj15_1.lon0 - atan2(-x, y);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L260:
	con = cosz - pj15_1.sinph0 * sin(geog[1]);
	if (abs(con) < epsln && abs(x) < epsln) {
	    return 0;
	}
	d__1 = pj15_1.lon0 + atan2(x * sinz * pj15_1.cosph0, con * rh);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj15z0_ */

/*                   PJ16Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*                           *  SINUSOIDAL  * */
/* ********************************************************************** */

/* Subroutine */ int pj16z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal halfpi = 1.5707963267948966;
    static doublereal epsln = 1e-10;

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double cos();

    /* Local variables */
    static doublereal geog[2], proj[2], x, y;
    extern doublereal adjlz0_();
    static doublereal con, lon;


/* **** PARAMETERS **** A,LON0,X0,Y0 ************************************ 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[15] != 0) {
	    goto L120;
	}
	errmz0_1.ierror = 162;
	return 0;
L120:
	d__1 = geog[0] - pj16_1.lon0;
	lon = adjlz0_(&d__1);
	proj[0] = pj16_1.x0 + pj16_1.a * lon * cos(geog[1]);
	proj[1] = pj16_1.y0 + pj16_1.a * geog[1];
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[15] != 0) {
	    goto L220;
	}
	errmz0_1.ierror = 163;
	return 0;
L220:
	x = proj[0] - pj16_1.x0;
	y = proj[1] - pj16_1.y0;
	geog[1] = y / pj16_1.a;
	if (abs(geog[1]) <= halfpi) {
	    goto L230;
	}
	errmz0_1.ierror = 164;
	return 0;
L230:
	con = abs(geog[1]) - halfpi;
	if (abs(con) > epsln) {
	    goto L240;
	}
	geog[0] = pj16_1.lon0;
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L240:
	d__1 = pj16_1.lon0 + x / (pj16_1.a * cos(geog[1]));
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj16z0_ */

/*                   PJ17Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*                  *  EQUIRECTANGULAR   * */
/* ********************************************************************** */

/* Subroutine */ int pj17z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal halfpi = 1.5707963267948966;

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double cos();

    /* Local variables */
    static doublereal geog[2], proj[2], x, y;
    extern doublereal adjlz0_();
    static doublereal lon;


/* **** PARAMETERS **** A,LON0,X0,Y0,LAT1 ******************************* 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[16] != 0) {
	    goto L120;
	}
	errmz0_1.ierror = 172;
	return 0;
L120:
	d__1 = geog[0] - pj17_1.lon0;
	lon = adjlz0_(&d__1);
	proj[0] = pj17_1.x0 + pj17_1.a * lon * cos(pj17_1.lat1);
	proj[1] = pj17_1.y0 + pj17_1.a * geog[1];
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[16] != 0) {
	    goto L220;
	}
	errmz0_1.ierror = 173;
	return 0;
L220:
	x = proj[0] - pj17_1.x0;
	y = proj[1] - pj17_1.y0;
	geog[1] = y / pj17_1.a;
	if (abs(geog[1]) <= halfpi) {
	    goto L240;
	}
	errmz0_1.ierror = 174;
	return 0;
L240:
	d__1 = pj17_1.lon0 + x / (pj17_1.a * cos(pj17_1.lat1));
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj17z0_ */

/*                   PJ18Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*                       *  MILLER CYLINDRICAL  * */
/* ********************************************************************** */

/* Subroutine */ int pj18z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal fortpi = .78539816339744833;
    static doublereal oneq = 1.25;
    static doublereal twoh = 2.5;

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double tan(), log(), exp(), atan();

    /* Local variables */
    static doublereal geog[2], proj[2], x, y;
    extern doublereal adjlz0_();
    static doublereal lon;


/* **** PARAMETERS **** A,LON0,X0,Y0 ************************************ 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[17] != 0) {
	    goto L120;
	}
	errmz0_1.ierror = 182;
	return 0;
L120:
	d__1 = geog[0] - pj18_1.lon0;
	lon = adjlz0_(&d__1);
	proj[0] = pj18_1.x0 + pj18_1.a * lon;
	proj[1] = pj18_1.y0 + pj18_1.a * log(tan(fortpi + geog[1] / twoh)) * 
		oneq;
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[17] != 0) {
	    goto L220;
	}
	errmz0_1.ierror = 183;
	return 0;
L220:
	x = proj[0] - pj18_1.x0;
	y = proj[1] - pj18_1.y0;
	d__1 = pj18_1.lon0 + x / pj18_1.a;
	geog[0] = adjlz0_(&d__1);
	geog[1] = twoh * atan(exp(y / pj18_1.a / oneq)) - fortpi * twoh;
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj18z0_ */

/*                   PJ19Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*                        *  VAN DER GRINTEN I  * */
/* ********************************************************************** */

/* Subroutine */ int pj19z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal pi = 3.14159265358979323846;
    static doublereal halfpi = 1.5707963267948966;
    static doublereal epsln = 1e-10;
    static doublereal half = .5;
    static doublereal one = 1.;
    static doublereal two = 2.;
    static doublereal three = 3.;

    /* System generated locals */
    doublereal d__1, d__2, d__3;

    /* Builtin functions */
    double tan(), d_sign(), sin(), cos(), sqrt(), acos();

    /* Local variables */
    static doublereal geog[2], proj[2], d, g, m, x, y, theta, a1, c1, c2, c3, 
	    m1;
    extern doublereal adjlz0_(), asinz0_();
    static doublereal al, xx, yy, costht, sintht, th1, con, lat, asq, lon, 
	    gsq, msq, xys;


/* **** PARAMETERS **** A,LON0,X0,Y0 ************************************ 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[18] != 0) {
	    goto L120;
	}
	errmz0_1.ierror = 192;
	return 0;
L120:
	d__1 = geog[0] - pj19_1.lon0;
	lon = adjlz0_(&d__1);
	lat = geog[1];
	if (abs(lat) > epsln) {
	    goto L140;
	}
	proj[0] = pj19_1.x0 + pj19_1.a * lon;
	proj[1] = pj19_1.y0;
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
L140:
/* Computing MIN */
	d__3 = (d__1 = lat / halfpi, abs(d__1));
	d__2 = min(d__3,one);
	theta = asinz0_(&d__2);
	if (abs(lon) > epsln && (d__1 = abs(lat) - halfpi, abs(d__1)) > epsln)
		 {
	    goto L160;
	}
	proj[0] = pj19_1.x0;
	d__1 = tan(half * theta);
	proj[1] = pj19_1.y0 + pi * pj19_1.a * d_sign(&d__1, &lat);
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
L160:
	al = half * (d__1 = pi / lon - lon / pi, abs(d__1));
	asq = al * al;
	sintht = sin(theta);
	costht = cos(theta);
	g = costht / (sintht + costht - one);
	gsq = g * g;
	m = g * (two / sintht - one);
	msq = m * m;
/* Computing 2nd power */
	d__1 = g - msq;
	con = pi * pj19_1.a * (al * (g - msq) + sqrt(asq * (d__1 * d__1) - (
		msq + asq) * (gsq - msq))) / (msq + asq);
	con = d_sign(&con, &lon);
	proj[0] = pj19_1.x0 + con;
	con = (d__1 = con / (pi * pj19_1.a), abs(d__1));
	d__1 = pi * pj19_1.a * sqrt(one - con * con - two * al * con);
	proj[1] = pj19_1.y0 + d_sign(&d__1, &lat);
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/*     ALGORITHM DEVELOPED BY D.P. RUBINCAM, THE AMERICAN CARTOGRAPHER, */

/*                1981, V. 8, NO. 2, P. 177-180. */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[18] != 0) {
	    goto L220;
	}
	errmz0_1.ierror = 193;
	return 0;
L220:
	x = proj[0] - pj19_1.x0;
	y = proj[1] - pj19_1.y0;
	con = pi * pj19_1.a;
	xx = x / con;
	yy = y / con;
	xys = xx * xx + yy * yy;
	c1 = -abs(yy) * (one + xys);
	c2 = c1 - two * yy * yy + xx * xx;
	c3 = -two * c1 + one + two * yy * yy + xys * xys;
	d = yy * yy / c3 + (two * c2 * c2 * c2 / c3 / c3 / c3 - c1 * 9. * c2 /
		 c3 / c3) / 27.;
	a1 = (c1 - c2 * c2 / three / c3) / c3;
	m1 = two * sqrt(-a1 / three);
	con = three * d / a1 / m1;
	if (abs(con) > one) {
	    con = d_sign(&one, &con);
	}
	th1 = acos(con) / three;
	geog[1] = (-m1 * cos(th1 + pi / three) - c2 / three / c3) * d_sign(&
		pi, &y);
	if (abs(xx) >= epsln) {
	    goto L230;
	}
	geog[0] = pj19_1.lon0;
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L230:
	geog[0] = pj19_1.lon0 + pi * (xys - one + sqrt(one + two * (xx * xx - 
		yy * yy) + xys * xys)) / two / xx;
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj19z0_ */

/*                   PJ20Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*                    *  OBLIQUE MERCATOR (HOTINE)  * */
/* ********************************************************************** */

/* Subroutine */ int pj20z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal pi = 3.14159265358979323846;
    static doublereal halfpi = 1.5707963267948966;
    static doublereal tol = 1e-7;
    static doublereal epsln = 1e-10;
    static doublereal zero = 0.;
    static doublereal half = .5;
    static doublereal one = 1.;

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double sin(), d_sign(), pow_dd(), cos(), atan(), log(), exp(), sqrt(), 
	    atan2();

    /* Local variables */
    static doublereal geog[2], dlon, proj[2], q, s, t, x, y;
    extern doublereal phi2z0_(), adjlz0_(), tsfnz0_();
    static doublereal ul, vl, us, ts, vs, sinphi, con, lon;


/* **** PARAMETERS **** A,E,ES,KS0,ALPHA,LONC,LON1,LAT1,LON2,LAT2,LAT0 ** 
*/
/* ********************** X0,Y0,GAMMA,LON0,AL,BL,EL ********************* 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[19] != 0) {
	    goto L220;
	}
	errmz0_1.ierror = 204;
	return 0;
L220:
	sinphi = sin(geog[1]);
	d__1 = geog[0] - pj20_1.lon0;
	dlon = adjlz0_(&d__1);
	vl = sin(pj20_1.bl * dlon);
	if ((d__1 = abs(geog[1]) - halfpi, abs(d__1)) > epsln) {
	    goto L230;
	}
	ul = pj20_1.singam * d_sign(&one, &geog[1]);
	us = pj20_1.al * geog[1] / pj20_1.bl;
	goto L250;
L230:
	ts = tsfnz0_(&pj20_1.e, &geog[1], &sinphi);
	q = pj20_1.el / pow_dd(&ts, &pj20_1.bl);
	s = half * (q - one / q);
	t = half * (q + one / q);
	ul = (s * pj20_1.singam - vl * pj20_1.cosgam) / t;
	con = cos(pj20_1.bl * dlon);
	if (abs(con) < tol) {
	    goto L240;
	}
	us = pj20_1.al * atan((s * pj20_1.cosgam + vl * pj20_1.singam) / con) 
		/ pj20_1.bl;
	if (con < zero) {
	    us += pi * pj20_1.al / pj20_1.bl;
	}
	goto L250;
L240:
	us = pj20_1.al * pj20_1.bl * dlon;
L250:
	if ((d__1 = abs(ul) - one, abs(d__1)) > epsln) {
	    goto L255;
	}
	errmz0_1.ierror = 205;
	return 0;
L255:
	vs = half * pj20_1.al * log((one - ul) / (one + ul)) / pj20_1.bl;
	us -= pj20_1.u0;
	proj[0] = pj20_1.x0 + vs * pj20_1.cosalf + us * pj20_1.sinalf;
	proj[1] = pj20_1.y0 + us * pj20_1.cosalf - vs * pj20_1.sinalf;
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[19] != 0) {
	    goto L280;
	}
	errmz0_1.ierror = 206;
	return 0;
L280:
	x = proj[0] - pj20_1.x0;
	y = proj[1] - pj20_1.y0;
	vs = x * pj20_1.cosalf - y * pj20_1.sinalf;
	us = y * pj20_1.cosalf + x * pj20_1.sinalf;
	us += pj20_1.u0;
	q = exp(-pj20_1.bl * vs / pj20_1.al);
	s = half * (q - one / q);
	t = half * (q + one / q);
	vl = sin(pj20_1.bl * us / pj20_1.al);
	ul = (vl * pj20_1.cosgam + s * pj20_1.singam) / t;
	if ((d__1 = abs(ul) - one, abs(d__1)) >= epsln) {
	    goto L300;
	}
	geog[0] = pj20_1.lon0;
	geog[1] = d_sign(&halfpi, &ul);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L300:
	con = one / pj20_1.bl;
	d__1 = pj20_1.el / sqrt((one + ul) / (one - ul));
	ts = pow_dd(&d__1, &con);
	geog[1] = phi2z0_(&pj20_1.e, &ts);
	con = cos(pj20_1.bl * us / pj20_1.al);
	lon = pj20_1.lon0 - atan2(s * pj20_1.cosgam - vl * pj20_1.singam, con)
		 / pj20_1.bl;
	geog[0] = adjlz0_(&lon);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj20z0_ */

/*                   PJ21Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*                       *       ROBINSON       * */
/* ********************************************************************** */

/* Subroutine */ int pj21z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal dg1 = .01745329252;
    static doublereal pi = 3.14159265358979323846;
    static doublereal epsln = 1e-10;

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double d_sign();

    /* Local variables */
    static doublereal geog[2], phid, proj[2], c, t, u, v, x, y, p2, _y1;
    extern doublereal adjlz0_();
    static integer nn;
    static doublereal yy;
    static integer ip1;
    static doublereal lon;


/* **** PARAMETERS **** A,LON0,X0,Y0 ************************************ 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[20] != 0) {
	    goto L120;
	}
	errmz0_1.ierror = 212;
	return 0;
L120:
	d__1 = geog[0] - pj21_1.lon0;
	lon = adjlz0_(&d__1);
	p2 = (d__1 = geog[1] / 5. / dg1, abs(d__1));
	ip1 = (integer) (p2 - epsln);

/*        STIRLING'S INTERPOLATION FORMULA (USING 2ND DIFF.) */
/*            USED WITH LOOKUP TABLE TO COMPUTE RECTANGULAR COORDINATE
S */
/*            FROM LAT/LONG. */

	p2 -= (doublereal) ip1;
	x = pj21_1.a * (pj21_1.xlr[ip1 + 1] + p2 * (pj21_1.xlr[ip1 + 2] - 
		pj21_1.xlr[ip1]) / 2. + p2 * p2 * (pj21_1.xlr[ip1 + 2] - 
		pj21_1.xlr[ip1 + 1] * 2. + pj21_1.xlr[ip1]) / 2.) * lon;
	y = pj21_1.a * (pj21_1.pr[ip1 + 1] + p2 * (pj21_1.pr[ip1 + 2] - 
		pj21_1.pr[ip1]) / 2. + p2 * p2 * (pj21_1.pr[ip1 + 2] - 
		pj21_1.pr[ip1 + 1] * 2. + pj21_1.pr[ip1]) / 2.) * pi / 2. * 
		d_sign(&c_b239, &geog[1]);
	proj[0] = pj21_1.x0 + x;
	proj[1] = pj21_1.y0 + y;
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[20] != 0) {
	    goto L220;
	}
	errmz0_1.ierror = 213;
	return 0;
L220:
	x = proj[0] - pj21_1.x0;
	y = proj[1] - pj21_1.y0;
	yy = y * 2. / pi / pj21_1.a;
	phid = yy * 90.;
	p2 = (d__1 = phid / 5., abs(d__1));
	ip1 = (integer) (p2 - epsln);
	if (ip1 == 0) {
	    ip1 = 1;
	}
	nn = 0;

/*     STIRLING'S INTERPOLATION FORMULA AS USED IN FORWARD TRANSFORMAT
ION */
/*     IS REVERSED FOR FIRST ESTIMATION OF LAT. FROM RECTANGULAR */
/*     COORDINATES.  LAT. IS THEN ADJUSTED BY ITERATION UNTIL USE OF 
*/
/*     FORWARD SERIES PROVIDES CORRECT VALUE OF Y WITHIN TOLERANCE. */


L230:
	u = pj21_1.pr[ip1 + 2] - pj21_1.pr[ip1];
	v = pj21_1.pr[ip1 + 2] - pj21_1.pr[ip1 + 1] * 2. + pj21_1.pr[ip1];
	t = (abs(yy) - pj21_1.pr[ip1 + 1]) * 2. / u;
	c = v / u;
	p2 = t * (1. - c * t * (1. - c * 2. * t));
	if (p2 < 0. && ip1 != 1) {
	    goto L240;
	}
	d__1 = (p2 + (doublereal) ip1) * 5.;
	phid = d_sign(&d__1, &y);
L235:
	p2 = (d__1 = phid / 5., abs(d__1));
	ip1 = (integer) (p2 - epsln);
	p2 -= (doublereal) ip1;
	_y1 = pj21_1.a * (pj21_1.pr[ip1 + 1] + p2 * (pj21_1.pr[ip1 + 2] - 
		pj21_1.pr[ip1]) / 2. + p2 * p2 * (pj21_1.pr[ip1 + 2] - 
		pj21_1.pr[ip1 + 1] * 2. + pj21_1.pr[ip1]) / 2.) * pi / 2. * 
		d_sign(&c_b239, &y);
	phid -= (_y1 - y) * 180. / pi / pj21_1.a;
	++nn;
	if (nn <= 20) {
	    goto L237;
	}
	errmz0_1.ierror = 214;
	return 0;
L237:
	if ((d__1 = _y1 - y, abs(d__1)) > 1e-5) {
	    goto L235;
	}
	goto L250;
L240:
	--ip1;
	goto L230;
/* L245: */
L250:
	geog[1] = phid * dg1;

/*        CALCULATE LONG. USING FINAL LAT. WITH TRANSPOSED FORWARD */
/*        STIRLING'S INTERPOLATION FORMULA. */

	geog[0] = pj21_1.lon0 + x / pj21_1.a / (pj21_1.xlr[ip1 + 1] + p2 * (
		pj21_1.xlr[ip1 + 2] - pj21_1.xlr[ip1]) / 2. + p2 * p2 * (
		pj21_1.xlr[ip1 + 2] - pj21_1.xlr[ip1 + 1] * 2. + pj21_1.xlr[
		ip1]) / 2.);
	geog[0] = adjlz0_(geog);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj21z0_ */

/*                   PJ22Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*                      *  SPACE OBLIQUE MERCATOR  * */
/* ********************************************************************** */

/* Subroutine */ int pj22z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal tol = 1e-7;
    static doublereal dg1 = .01745329252;
    static doublereal pi = 3.14159265358979323846;
    static doublereal one = 1.;
    static doublereal two = 2.;

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double cos(), d_sign(), sin(), tan(), atan(), sqrt(), log(), exp();

    /* Local variables */
    static doublereal geog[2], lamt, xlam, proj[2], sdsq, c, d;
    static integer l;
    static doublereal s, x, y, actan, lamdp, phidp, lampp, tanph, lamtp, 
	    sppsq;
    extern doublereal asinz0_();
    static doublereal dd, cl, sd;
    static integer nn;
    static doublereal sl, sp, fac, dif, lat, scl, lon, sav, rlm, spp, rlm2;


/* **** PARAMETERS **** A,E,ES,LON0,LATC,X0,Y0,MCS,TCS,FAC,IND ********** 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[21] != 0) {
	    goto L220;
	}
	errmz0_1.ierror = 222;
	return 0;
L220:
	if (pj22_1.land >= 4) {
	    goto L225;
	}
	lon = geog[0] - dg1 * 128.87 + pi * two / 251. * (doublereal) 
		pj22_1.path;
	goto L230;
L225:
	lon = geog[0] - dg1 * 129.3 + pi * two / 233. * (doublereal) 
		pj22_1.path;
L230:
	lat = geog[1];

/*        TEST FOR LAT. AND LONG. APPROACHING 90 DEGREES. */

	if (lat > 1.570796) {
	    lat = 1.570796;
	}
	if (lat < -1.570796) {
	    lat = -1.570796;
	}
	if (lat >= 0.) {
	    lampp = pi / two;
	}
	if (lat < 0.) {
	    lampp = pi * 1.5;
	}
	nn = 0;
L231:
	sav = lampp;
	l = 0;
	lamtp = lon + norm_1.p22 * lampp;
	cl = cos(lamtp);
	if (abs(cl) < tol) {
	    lamtp -= tol;
	}
	fac = lampp - d_sign(&one, &cl) * sin(lampp) * pi / two;
L232:
	lamt = lon + norm_1.p22 * sav;
	c = cos(lamt);
	if (abs(c) < tol) {
	    lamdp = sav;
	    goto L233;
	}
	xlam = ((one - norm_1.es) * tan(lat) * norm_1.sa + sin(lamt) * 
		norm_1.ca) / c;
	lamdp = atan(xlam);
	lamdp += fac;
	dif = abs(sav) - abs(lamdp);
	if (abs(dif) < tol) {
	    goto L233;
	}
	sav = lamdp;
	++l;
	if (l > 50) {
	    goto L234;
	}
	goto L232;

/*        ADJUST FOR LANDSAT ORIGIN. */

L233:
	rlm = pi * (one / 248. + .5161290322580645);
	rlm2 = rlm + two * pi;
	++nn;
	if (nn >= 3) {
	    goto L236;
	}
	if (lamdp > rlm && lamdp < rlm2) {
	    goto L236;
	}
	if (lamdp <= rlm) {
	    lampp = pi * 2.5;
	}
	if (lamdp >= rlm2) {
	    lampp = pi / two;
	}
	goto L231;
L234:
	errmz0_1.ierror = 223;
L236:

/*        LAMDP COMPUTED.  NOW COMPUTE PHIDP. */

	sp = sin(lat);
	d__1 = ((one - norm_1.es) * norm_1.ca * sp - norm_1.sa * cos(lat) * 
		sin(lamt)) / sqrt(one - norm_1.es * sp * sp);
	phidp = asinz0_(&d__1);

/*        COMPUTE X AND Y */

	tanph = log(tan(pi / 4. + phidp / two));
	sd = sin(lamdp);
	sdsq = sd * sd;
	s = norm_1.p22 * norm_1.sa * cos(lamdp) * sqrt((one + norm_1.t * sdsq)
		 / ((one + norm_1.w * sdsq) * (one + norm_1.q * sdsq)));
	d = sqrt(norm_1.xj * norm_1.xj + s * s);
	x = pj22_1.b * lamdp + pj22_1.a2 * sin(two * lamdp) + pj22_1.a4 * sin(
		lamdp * 4.) - tanph * s / d;
	x = pj22_1.a * x;
	y = pj22_1.c1 * sd + pj22_1.c3 * sin(lamdp * 3.) + tanph * norm_1.xj /
		 d;
	y = pj22_1.a * y;
	proj[0] = x + pj22_1.x0;
	proj[1] = y + pj22_1.y0;
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[21] != 0) {
	    goto L320;
	}
	errmz0_1.ierror = 224;
	return 0;
L320:
	x = proj[0] - pj22_1.x0;
	y = proj[1] - pj22_1.y0;

/*        COMPUTE TRANSFORMED LAT/LONG AND GEODETIC LAT/LONG, GIVEN X,
Y. */

/*        BEGIN INVERSE COMPUTATION WITH APPROXIMATION FOR LAMDP.  SOL
VE */
/*        FOR TRANSFORMED LONG. */

	lamdp = x / pj22_1.a / pj22_1.b;
	nn = 0;
L325:
	sav = lamdp;
	sd = sin(lamdp);
	sdsq = sd * sd;
	s = norm_1.p22 * norm_1.sa * cos(lamdp) * sqrt((one + norm_1.t * sdsq)
		 / ((one + norm_1.w * sdsq) * (one + norm_1.q * sdsq)));
	lamdp = x / pj22_1.a + y / pj22_1.a * s / norm_1.xj - pj22_1.a2 * sin(
		two * lamdp) - pj22_1.a4 * sin(lamdp * 4.) - s / norm_1.xj * (
		pj22_1.c1 * sin(lamdp) + pj22_1.c3 * sin(lamdp * 3.));
	lamdp /= pj22_1.b;
	dif = lamdp - sav;
	if (abs(dif) < tol) {
	    goto L330;
	}
	++nn;
	if (nn == 50) {
	    goto L330;
	}
	goto L325;

/*        COMPUTE TRANSFORMED LAT. */

L330:
	sl = sin(lamdp);
	fac = exp(sqrt(one + s * s / norm_1.xj / norm_1.xj) * (y / pj22_1.a - 
		pj22_1.c1 * sl - pj22_1.c3 * sin(lamdp * 3.)));
	actan = atan(fac);
	phidp = two * (actan - pi / 4.);

/*        COMPUTE GEODETIC LATITUDE. */

	dd = sl * sl;
	if ((d__1 = cos(lamdp), abs(d__1)) < tol) {
	    lamdp -= tol;
	}
	spp = sin(phidp);
	sppsq = spp * spp;
	lamt = atan(((one - sppsq / (one - norm_1.es)) * tan(lamdp) * 
		norm_1.ca - spp * norm_1.sa * sqrt((one + norm_1.q * dd) * (
		one - sppsq) - sppsq * norm_1.u) / cos(lamdp)) / (one - sppsq 
		* (one + norm_1.u)));

/*        CORRECT INVERSE QUADRANT. */

	if (lamt >= 0.) {
	    sl = one;
	}
	if (lamt < 0.) {
	    sl = -one;
	}
	if (cos(lamdp) >= 0.) {
	    scl = one;
	}
	if (cos(lamdp) < 0.) {
	    scl = -one;
	}
	lamt -= pi / two * (one - scl) * sl;
	lon = lamt - norm_1.p22 * lamdp;

/*        COMPUTE GEODETIC LATITUDE. */

	if (abs(norm_1.sa) < tol) {
	    d__1 = spp / sqrt((one - norm_1.es) * (one - norm_1.es) + 
		    norm_1.es * sppsq);
	    lat = asinz0_(&d__1);
	}
	if (abs(norm_1.sa) < tol) {
	    goto L335;
	}
	lat = atan((tan(lamdp) * cos(lamt) - norm_1.ca * sin(lamt)) / ((one - 
		norm_1.es) * norm_1.sa));
L335:
	if (pj22_1.land >= 4) {
	    goto L370;
	}
	geog[0] = lon + dg1 * 128.87 - pi * two / 251. * (doublereal) 
		pj22_1.path;
	goto L380;
L370:
	geog[0] = lon + dg1 * 129.3 - pi * two / 233. * (doublereal) 
		pj22_1.path;
L380:
	geog[1] = lat;
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

/*###192 [f77] Error on line 192 of pj22z0.f missing statement label 234%%
%*/
} /* pj22z0_ */

/*                   PJ23Z0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* **              MATHEMATICAL ANALYSIS BY JOHN SNYDER                ** */
/* ********************************************************************** */
/*            * MODIFIED-STEREOGRAPHIC CONFORMAL (FOR ALASKA) * */
/* ********************************************************************** */

/* Subroutine */ int pj23z0_(coord, crdio, indic)
doublereal *coord, *crdio;
shortint *indic;
{
    /* Initialized data */

    static doublereal halfpi = 1.5707963267948966;
    static doublereal epsln = 1e-10;
    static doublereal one = 1.;
    static doublereal two = 2.;

    /* System generated locals */
    integer i__1;
    doublereal d__1, d__2;

    /* Builtin functions */
    double sin(), cos(), tan(), pow_dd(), atan(), sqrt(), atan2();

    /* Local variables */
    static doublereal cchi, geog[2], dphi, schi, proj[2], cosz, fxyi, sinz, 
	    fxyr, g;
    static integer j;
    static doublereal r, s, x, y, z, esphi, fpxyi, fpxyr;
    extern doublereal adjlz0_(), asinz0_();
    static doublereal ai, bi, ci, di, ar, br, cr, dr, ds, rh;
    static integer nn;
    static doublereal xp, yp, coslon, sinlon, chi, den, ain, cin, arn, crn, 
	    phi, lon, dxp, dyp;


/* **** PARAMETERS **** A,E,ES,LON0,LAT0,X0,Y0,SINPH0,COSPH0 ************ 
*/
    /* Parameter adjustments */
    --crdio;
    --coord;

    /* Function Body */

/* ...................................................................... 
*/
/*                      .  FORWARD TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 0) {

	geog[0] = coord[1];
	geog[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[22] != 0) {
	    goto L120;
	}
	errmz0_1.ierror = 232;
	return 0;
L120:
	d__1 = geog[0] - pj23_1.lon0;
	lon = adjlz0_(&d__1);

/*     CALCULATE X-PRIME AND Y-PRIME FOR OBLIQUE STEREOGRAPHIC PROJ. 
*/
/*          FROM LAT/LONG. */

	sinlon = sin(lon);
	coslon = cos(lon);
	esphi = pj23_1.ec * sin(geog[1]);
	d__1 = (one - esphi) / (one + esphi);
	d__2 = pj23_1.ec / two;
	chi = two * atan(tan((halfpi + geog[1]) / two) * pow_dd(&d__1, &d__2))
		 - halfpi;
	schi = sin(chi);
	cchi = cos(chi);
	g = pj23_1.schio * schi + pj23_1.cchio * cchi * coslon;
	s = two / (one + g);
	xp = s * cchi * sinlon;
	yp = s * (pj23_1.cchio * schi - pj23_1.schio * cchi * coslon);

/*     USE KNUTH ALGORITHM FOR SUMMING COMPLEX TERMS, TO CONVERT */
/*     OBLIQUE STEREOGRAPHIC TO MODIFIED-STEREOGRAPHIC COORD. */

	r = xp + xp;
	s = xp * xp + yp * yp;
	ar = pj23_1.acoef[pj23_1.n - 1];
	ai = pj23_1.bcoef[pj23_1.n - 1];
	br = pj23_1.acoef[pj23_1.n - 2];
	bi = pj23_1.bcoef[pj23_1.n - 2];
	i__1 = pj23_1.n;
	for (j = 2; j <= i__1; ++j) {
	    arn = br + r * ar;
	    ain = bi + r * ai;
	    if (j == pj23_1.n) {
		goto L140;
	    }
	    br = pj23_1.acoef[pj23_1.n - j - 1] - s * ar;
	    bi = pj23_1.bcoef[pj23_1.n - j - 1] - s * ai;
	    ar = arn;
	    ai = ain;
L140:
	    ;
	}
	br = -s * ar;
	bi = -s * ai;
	ar = arn;
	ai = ain;
	x = xp * ar - yp * ai + br;
	y = yp * ar + xp * ai + bi;
	proj[0] = x * pj23_1.a + pj23_1.x0;
	proj[1] = y * pj23_1.a + pj23_1.y0;
	crdio[1] = proj[0];
	crdio[2] = proj[1];
	return 0;
    }

/* ...................................................................... 
*/
/*                      .  INVERSE TRANSFORMATION  . */
/* ...................................................................... 
*/

    if (*indic == 1) {

	proj[0] = coord[1];
	proj[1] = coord[2];
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[22] != 0) {
	    goto L220;
	}
	errmz0_1.ierror = 234;
	return 0;
L220:
	x = (proj[0] - pj23_1.x0) / pj23_1.a;
	y = (proj[1] - pj23_1.y0) / pj23_1.a;
	xp = x;
	yp = y;
	nn = 0;

/*     USE KNUTH ALGORITHM FOR SUMMING COMPLEX TERMS, TO CONVERT */
/*     MODIFIED-STEREOGRAPHIC CONFORMAL TO OBLIQUE STEREOGRAPHIC */
/*     COORDINATES (XP,YP). */

L225:
	r = xp + xp;
	s = xp * xp + yp * yp;
	ar = pj23_1.acoef[pj23_1.n - 1];
	ai = pj23_1.bcoef[pj23_1.n - 1];
	br = pj23_1.acoef[pj23_1.n - 2];
	bi = pj23_1.bcoef[pj23_1.n - 2];
	cr = (doublereal) pj23_1.n * ar;
	ci = (doublereal) pj23_1.n * ai;
	dr = (doublereal) (pj23_1.n - 1) * br;
	di = (doublereal) (pj23_1.n - 1) * bi;
	i__1 = pj23_1.n;
	for (j = 2; j <= i__1; ++j) {
	    arn = br + r * ar;
	    ain = bi + r * ai;
	    if (j == pj23_1.n) {
		goto L230;
	    }
	    br = pj23_1.acoef[pj23_1.n - j - 1] - s * ar;
	    bi = pj23_1.bcoef[pj23_1.n - j - 1] - s * ai;
	    ar = arn;
	    ai = ain;
	    crn = dr + r * cr;
	    cin = di + r * ci;
	    dr = (doublereal) (pj23_1.n - j) * pj23_1.acoef[pj23_1.n - j - 1] 
		    - s * cr;
	    di = (doublereal) (pj23_1.n - j) * pj23_1.bcoef[pj23_1.n - j - 1] 
		    - s * ci;
	    cr = crn;
	    ci = cin;
L230:
	    ;
	}
	br = -s * ar;
	bi = -s * ai;
	ar = arn;
	ai = ain;
	fxyr = xp * ar - yp * ai + br - x;
	fxyi = yp * ar + xp * ai + bi - y;
	fpxyr = xp * cr - yp * ci + dr;
	fpxyi = yp * cr + xp * ci + di;
	den = fpxyr * fpxyr + fpxyi * fpxyi;
	dxp = -(fxyr * fpxyr + fxyi * fpxyi) / den;
	dyp = -(fxyi * fpxyr - fxyr * fpxyi) / den;
	xp += dxp;
	yp += dyp;
	ds = abs(dxp) + abs(dyp);
	++nn;
	if (nn <= 20) {
	    goto L237;
	}
	errmz0_1.ierror = 235;
	goto L238;
L237:
	if (ds > epsln) {
	    goto L225;
	}

/*     CONVERT OBLIQUE STEREOGRAPHIC COORDINATES TO LAT/LONG. */

L238:
	rh = sqrt(xp * xp + yp * yp);
	z = two * atan(rh / two);
	sinz = sin(z);
	cosz = cos(z);
	geog[0] = pj23_1.lon0;
	if (abs(rh) > epsln) {
	    goto L240;
	}
	geog[1] = pj23_1.lat0;
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
L240:
	d__1 = cosz * pj23_1.schio + yp * sinz * pj23_1.cchio / rh;
	chi = asinz0_(&d__1);
	nn = 0;
	phi = chi;
L250:
	esphi = pj23_1.ec * sin(phi);
	d__1 = (one + esphi) / (one - esphi);
	d__2 = pj23_1.ec / two;
	dphi = two * atan(tan((halfpi + chi) / two) * pow_dd(&d__1, &d__2)) - 
		halfpi - phi;
	phi += dphi;
	++nn;
	if (nn <= 20) {
	    goto L257;
	}
	errmz0_1.ierror = 236;
	goto L260;
L257:
	if (abs(dphi) > epsln) {
	    goto L250;
	}
L260:
	geog[1] = phi;
	d__1 = pj23_1.lon0 + atan2(xp * sinz, rh * pj23_1.cchio * cosz - yp * 
		pj23_1.schio * sinz);
	geog[0] = adjlz0_(&d__1);
	crdio[1] = geog[0];
	crdio[2] = geog[1];
	return 0;
    }

    return 0;

} /* pj23z0_ */

/*                   PJINIT */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

/* Subroutine */ int pjinit_(isys, zone, data)
integer *isys, *zone;
doublereal *data;
{
    /* Initialized data */

    static doublereal pi = 3.14159265358979323846;
    static doublereal halfpi = 1.5707963267948966;
    static doublereal zero = 0.;
    static doublereal half = .5;
    static doublereal one = 1.;
    static doublereal two = 2.;
    static doublereal epsln = 1e-10;
    static doublereal tol = 1e-7;
    static doublereal tol09 = 1e-5;
    static doublereal nintyd = 9e7;
    static doublereal dg1 = .01745329252;
    static integer nad27[134] = { 101,102,5010,5300,201,202,203,301,302,401,
	    402,403,404,405,406,407,501,502,503,600,700,901,902,903,1001,1002,
	    5101,5102,5103,5104,5105,1101,1102,1103,1201,1202,1301,1302,1401,
	    1402,1501,1502,1601,1602,1701,1702,1703,1801,1802,1900,2001,2002,
	    2101,2102,2103,2111,2112,2113,2201,2202,2203,2301,2302,2401,2402,
	    2403,2501,2502,2503,2601,2602,2701,2702,2703,2800,2900,3001,3002,
	    3003,3101,3102,3103,3104,3200,3301,3302,3401,3402,3501,3502,3601,
	    3602,3701,3702,3800,3901,3902,4001,4002,4100,4201,4202,4203,4204,
	    4205,4301,4302,4303,4400,4501,4502,4601,4602,4701,4702,4801,4802,
	    4803,4901,4902,4903,4904,5001,5002,5003,5004,5005,5006,5007,5008,
	    5009,5201,5202,5400 };
    static integer nad83[134] = { 101,102,5010,5300,201,202,203,301,302,401,
	    402,403,404,405,406,0,501,502,503,600,700,901,902,903,1001,1002,
	    5101,5102,5103,5104,5105,1101,1102,1103,1201,1202,1301,1302,1401,
	    1402,1501,1502,1601,1602,1701,1702,1703,1801,1802,1900,2001,2002,
	    2101,2102,2103,2111,2112,2113,2201,2202,2203,2301,2302,2401,2402,
	    2403,2500,0,0,2600,0,2701,2702,2703,2800,2900,3001,3002,3003,3101,
	    3102,3103,3104,3200,3301,3302,3401,3402,3501,3502,3601,3602,3701,
	    3702,3800,3900,0,4001,4002,4100,4201,4202,4203,4204,4205,4301,
	    4302,4303,4400,4501,4502,4601,4602,4701,4702,4801,4802,4803,4901,
	    4902,4903,4904,5001,5002,5003,5004,5005,5006,5007,5008,5009,5200,
	    0,5400 };

    /* System generated locals */
    doublereal d__1, d__2;

    /* Builtin functions */
    /* Subroutine */ int s_copy();
    double sqrt(), sin(), cos(), log(), pow_dd(), d_sign(), tan(), atan();

    /* Local variables */
    static integer ind02, degs[5], mode;
    static doublereal latc;
    static char sgna[1*5];
    static doublereal chio, lonc;
    static real secs[5];
    static doublereal save, dlon;
    static integer mins[5];
    static doublereal sumb, lat003, save9, lat004, lat007, a, d, e, f, g;
    static integer i;
    static doublereal h, j, l, _gamma, suma2, table[9], alpha, suma4, buffl[15]
	    , sumc1, sumc3;
    extern doublereal e0fnz0_(), e1fnz0_(), e2fnz0_();
    static char datum[32];
    static integer itemp, limit, lunit;
    extern doublereal e3fnz0_(), e4fnz0_();
    static doublereal sinp03, cosp03, sinp04, cosp04, p1, p2, esphi;
    extern /* Subroutine */ int raddz0_();
    extern doublereal adjlz0_(), pakcz0_();
    static doublereal cosph0, sinph0;
    static integer id;
    static doublereal es;
    extern doublereal pakrz0_(), msfnz0_(), qsfnz0_(), tsfnz0_();
    static doublereal cosphi;
    extern doublereal mlfnz0_();
    static doublereal e09, sinphi;
    static integer keepzn, ind;
    static doublereal con, ms1, qs1, ms2, qs2, qs0, ts1, ts2, ts0, ml1, ml2, 
	    ml0, ks0, com;
    extern doublereal asinz0_();
    static doublereal alf, esc, ess, lam;
    extern /* Subroutine */ int seraz0_();
    static doublereal fb, fa2, fa4, fc1, fc3, con1, lat2, lat0, ec2, lon1, 
	    lon2;




    /* Parameter adjustments */
    --data;

    /* Function Body */


/* .................................................................... */

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                              .  U T M  . */
/* ...................................................................... 
*/

    if (*isys == 1) {

	errmz0_1.ierror = 0;
	if (toggle_1.switch__[0] != 0 && toggle_1.switch__[0] == *zone) {
	    return 0;
	}
	toggle_1.switch__[0] = *zone;
	if (toggle_1.switch__[8] != 0 && toggle_1.switch__[8] == *zone && 
		data[14] == save) {
	    return 0;
	}
	keepzn = *zone;
	*zone = abs(*zone);
	save = data[1];
	if (*zone == 0) {
	    *zone = (integer) ((data[1] * 180. / pi + tol09 / 3600.) / 6.);
	    ind = 1;
	    if (data[1] < zero) {
		ind = 0;
	    }
	    *zone = (*zone + 30) % 60 + ind;
	    keepzn = *zone;
	    if (data[2] < zero) {
		keepzn = -(*zone);
	    }
	}
	if (*zone < 1 || *zone > 60) {
	    errmz0_1.ierror = 11;
	    return 0;
	}
	buffl[0] = data[14];
	buffl[1] = data[15];
	buffl[2] = .9996;
	buffl[3] = zero;
	buffl[4] = (doublereal) (*zone * 6 - 183) * 1e6;
	buffl[5] = zero;
	buffl[6] = 5e5;
	buffl[7] = zero;
	if (data[2] < zero) {
	    buffl[7] = 1e7;
	}
	if (keepzn < 0) {
	    buffl[7] = 1e7;
	}
	if (buffl[0] != 0. && buffl[0] != save9) {
	    toggle_1.switch__[8] = 0;
	}
	save9 = buffl[0];
	itemp = prinz0_1.ipparm;
	prinz0_1.ipparm = 1;
	for (i = 1; i <= 8; ++i) {
	    data[i] = buffl[i - 1];
/* L145: */
	}
	ellpz0_1.az = data[14];
	ellpz0_1.ez = data[15];
	toggle_1.switch__[8] = 0;
	goto L900;
L170:
	prinz0_1.ipparm = itemp;
	buffl[0] = pj09_2.a09;
	buffl[1] = pj09_2.es09;
	*zone = keepzn;
	toggle_1.switch__[8] = *zone;
	if (errmz0_1.ierror != 0) {
	    return 0;
	}

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	toggle_1.switch__[0] = *zone;
	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                           .  STATE PLANE  . */
/* ...................................................................... 
*/

    if (*isys == 2) {

	errmz0_1.ierror = 0;
	if (toggle_1.switch__[1] != 0 && toggle_1.switch__[1] == *zone) {
	    return 0;
	}
	if (spcs_1.ispher != 0 && spcs_1.ispher != 8) {
	    errmz0_1.ierror = 20;
	    return 0;
	}
	if (*zone > 0) {
	    ind02 = 0;
	    if (spcs_1.ispher == 0) {
		for (i = 1; i <= 134; ++i) {
		    if (*zone == nad27[i - 1]) {
			ind02 = i;
		    }
/* L210: */
		}
	    }
	    if (spcs_1.ispher == 8) {
		for (i = 1; i <= 134; ++i) {
		    if (*zone == nad83[i - 1]) {
			ind02 = i;
		    }
/* L220: */
		}
	    }
	    if (ind02 == 0) {
		errmz0_1.ierror = 21;
		return 0;
	    }
	} else {
	    errmz0_1.ierror = 21;
	    return 0;
	}
	if (spcs_1.ispher == 0) {
	    lunit = spcs_1.lu27;
	    s_copy(datum, spcs_1.file27, 32L, 32L);
	}
	if (spcs_1.ispher == 8) {
	    lunit = spcs_1.lu83;
	    s_copy(datum, spcs_1.file83, 32L, 32L);
	}
	if (id <= 0) {
	    errmz0_1.ierror = 21;
	    return 0;
	}
	pj02_1.itype = id;
	ellpz0_1.az = table[0];
	es = table[1];
	ellpz0_1.esz = es;
	ellpz0_1.ez = sqrt(es);
	ellpz0_1.e0z = e0fnz0_(&es);
	ellpz0_1.e1z = e1fnz0_(&es);
	ellpz0_1.e2z = e2fnz0_(&es);
	ellpz0_1.e3z = e3fnz0_(&es);
	ellpz0_1.e4z = e4fnz0_(&ellpz0_1.ez);
	itemp = prinz0_1.ipparm;
	prinz0_1.ipparm = 1;

/*     TRANSVERSE MERCATOR PROJECTION */

	if (pj02_1.itype == 1) {
	    data[3] = table[3];
	    data[5] = pakcz0_(&table[2]);
	    data[6] = pakcz0_(&table[6]);
	    data[7] = table[7];
	    data[8] = table[8];
	    spcs_1.msys = 9;
	    toggle_1.switch__[spcs_1.msys - 1] = 0;
	    goto L900;
	}

/*     LAMBERT CONFORMAL PROJECTION */

	if (pj02_1.itype == 2) {
	    data[3] = pakcz0_(&table[5]);
	    data[4] = pakcz0_(&table[4]);
	    data[5] = pakcz0_(&table[2]);
	    data[6] = pakcz0_(&table[6]);
	    data[7] = table[7];
	    data[8] = table[8];
	    spcs_1.msys = 4;
	    toggle_1.switch__[spcs_1.msys - 1] = 0;
	    goto L400;
	}

/*     POLYCONIC PROJECTION */

	if (pj02_1.itype == 3) {
	    data[5] = pakcz0_(&table[2]);
	    data[6] = pakcz0_(&table[3]);
	    data[7] = table[4];
	    data[8] = table[5];
	    spcs_1.msys = 7;
	    toggle_1.switch__[spcs_1.msys - 1] = 0;
	    goto L700;
	}

/*     OBLIQUE MERCATOR PROJECTION */

	if (pj02_1.itype == 4) {
	    data[3] = table[3];
	    data[4] = pakcz0_(&table[5]);
	    data[5] = pakcz0_(&table[2]);
	    data[6] = pakcz0_(&table[6]);
	    data[7] = table[7];
	    data[8] = table[8];
	    data[13] = one;
	    spcs_1.msys = 20;
	    toggle_1.switch__[spcs_1.msys - 1] = 0;
	    goto L2000;
	}

L270:
	prinz0_1.ipparm = itemp;
	if (errmz0_1.ierror != 0) {
	    return 0;
	}

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	toggle_1.switch__[1] = *zone;
	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                    .  ALBERS CONICAL EQUAL AREA  . */
/* ...................................................................... 
*/

    if (*isys == 3) {

	errmz0_1.ierror = 0;
	if (toggle_1.switch__[2] != 0 && toggle_1.switch__[2] == *zone) {
	    return 0;
	}
	toggle_1.switch__[2] = 0;
	pj03_2.a03 = ellpz0_1.az;
	pj03_2.e03 = ellpz0_1.ez;
	pj03_2.es03 = ellpz0_1.esz;
	pj17_2.lat1 = pakrz0_(&data[3]);
	lat2 = pakrz0_(&data[4]);
	if ((d__1 = pj17_2.lat1 + lat2, abs(d__1)) < epsln) {
	    errmz0_1.ierror = 31;
	    return 0;
	}
	pj03_2.lon003 = pakrz0_(&data[5]);
	lat003 = pakrz0_(&data[6]);
	pj03_2.x003 = data[7];
	pj03_2.y003 = data[8];
	sinp03 = sin(pj17_2.lat1);
	con = sinp03;
	cosp03 = cos(pj17_2.lat1);
	ms1 = msfnz0_(&pj03_2.e03, &sinp03, &cosp03);
	qs1 = qsfnz0_(&pj03_2.e03, &sinp03, &cosp03);
	sinp03 = sin(lat2);
	cosp03 = cos(lat2);
	ms2 = msfnz0_(&pj03_2.e03, &sinp03, &cosp03);
	qs2 = qsfnz0_(&pj03_2.e03, &sinp03, &cosp03);
	sinp03 = sin(lat003);
	cosp03 = cos(lat003);
	qs0 = qsfnz0_(&pj03_2.e03, &sinp03, &cosp03);
	if ((d__1 = pj17_2.lat1 - lat2, abs(d__1)) >= epsln) {
	    pj03_2.ns03 = (ms1 * ms1 - ms2 * ms2) / (qs2 - qs1);
	} else {
	    pj03_2.ns03 = con;
	}
	pj03_2.c = ms1 * ms1 + pj03_2.ns03 * qs1;
	pj03_2.rh003 = pj03_2.a03 * sqrt(pj03_2.c - pj03_2.ns03 * qs0) / 
		pj03_2.ns03;

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	raddz0_(&pj17_2.lat1, sgna, degs, mins, secs, 1L);
	raddz0_(&lat2, sgna + 1, &degs[1], &mins[1], &secs[1], 1L);
	raddz0_(&pj03_2.lon003, sgna + 2, &degs[2], &mins[2], &secs[2], 1L);
	raddz0_(&lat003, sgna + 3, &degs[3], &mins[3], &secs[3], 1L);
	data[1] = pj03_2.a03;
	data[2] = pj03_2.es03;
	toggle_1.switch__[2] = *zone;
	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                     .  LAMBERT CONFORMAL CONIC  . */
/* ...................................................................... 
*/

    if (*isys == 4) {

L400:
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[3] != 0 && toggle_1.switch__[3] == *zone) {
	    return 0;
	}
	toggle_1.switch__[3] = 0;
	pj04_2.a04 = ellpz0_1.az;
	pj04_2.e04 = ellpz0_1.ez;
	es = ellpz0_1.esz;
	pj17_2.lat1 = pakrz0_(&data[3]);
	lat2 = pakrz0_(&data[4]);
	if ((d__1 = pj17_2.lat1 + lat2, abs(d__1)) < epsln) {
	    errmz0_1.ierror = 41;
	    return 0;
	}
	pj04_2.lon004 = pakrz0_(&data[5]);
	lat004 = pakrz0_(&data[6]);
	pj04_2.x004 = data[7];
	pj04_2.y004 = data[8];
	sinp04 = sin(pj17_2.lat1);
	con = sinp04;
	cosp04 = cos(pj17_2.lat1);
	ms1 = msfnz0_(&pj04_2.e04, &sinp04, &cosp04);
	ts1 = tsfnz0_(&pj04_2.e04, &pj17_2.lat1, &sinp04);
	sinp04 = sin(lat2);
	cosp04 = cos(lat2);
	ms2 = msfnz0_(&pj04_2.e04, &sinp04, &cosp04);
	ts2 = tsfnz0_(&pj04_2.e04, &lat2, &sinp04);
	sinp04 = sin(lat004);
	ts0 = tsfnz0_(&pj04_2.e04, &lat004, &sinp04);
	if ((d__1 = pj17_2.lat1 - lat2, abs(d__1)) >= epsln) {
	    pj04_2.ns04 = log(ms1 / ms2) / log(ts1 / ts2);
	} else {
	    pj04_2.ns04 = con;
	}
	pj04_2.f04 = ms1 / (pj04_2.ns04 * pow_dd(&ts1, &pj04_2.ns04));
	pj04_2.rh004 = pj04_2.a04 * pj04_2.f04 * pow_dd(&ts0, &pj04_2.ns04);

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	raddz0_(&pj17_2.lat1, sgna, degs, mins, secs, 1L);
	raddz0_(&lat2, sgna + 1, &degs[1], &mins[1], &secs[1], 1L);
	raddz0_(&pj04_2.lon004, sgna + 2, &degs[2], &mins[2], &secs[2], 1L);
	raddz0_(&lat004, sgna + 3, &degs[3], &mins[3], &secs[3], 1L);
	data[1] = pj04_2.a04;
	data[2] = es;
	toggle_1.switch__[3] = *zone;

/*     RETURN TO STATE PLANE INITIALIZATION IF NECESSARY */

	if (*isys == 2) {
	    goto L270;
	}

	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                            .  MERCATOR  . */
/* ...................................................................... 
*/

    if (*isys == 5) {

	errmz0_1.ierror = 0;
	if (toggle_1.switch__[4] != 0 && toggle_1.switch__[4] == *zone) {
	    return 0;
	}
	toggle_1.switch__[4] = 0;
	pj05_2.a05 = ellpz0_1.az;
	pj05_2.e05 = ellpz0_1.ez;
	es = ellpz0_1.esz;
	pj05_2.lon005 = pakrz0_(&data[5]);
	pj17_2.lat1 = pakrz0_(&data[6]);
/* Computing 2nd power */
	d__1 = sin(pj17_2.lat1);
	pj05_2.m1 = cos(pj17_2.lat1) / sqrt(one - es * (d__1 * d__1));
	pj05_2.x005 = data[7];
	pj05_2.y005 = data[8];

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	raddz0_(&pj17_2.lat1, sgna, degs, mins, secs, 1L);
	raddz0_(&pj05_2.lon005, sgna + 1, &degs[1], &mins[1], &secs[1], 1L);
	data[1] = pj05_2.a05;
	data[2] = es;
	toggle_1.switch__[4] = *zone;
	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                       .  POLAR STEREOGRAPHIC  . */
/* ...................................................................... 
*/

    if (*isys == 6) {

	errmz0_1.ierror = 0;
	if (toggle_1.switch__[5] != 0 && toggle_1.switch__[5] == *zone) {
	    return 0;
	}
	toggle_1.switch__[5] = 0;
	pj06_2.a06 = ellpz0_1.az;
	pj06_2.e06 = ellpz0_1.ez;
	es = ellpz0_1.esz;
	pj06_2.e4 = ellpz0_1.e4z;
	pj06_2.lon006 = pakrz0_(&data[5]);
	save = data[6];
	latc = pakrz0_(&save);
	pj06_2.x006 = data[7];
	pj06_2.y006 = data[8];
	pj06_2.fac = one;
	if (save < zero) {
	    pj06_2.fac = -one;
	}
	pj06_2.ind06 = 0;
	if (abs(save) != nintyd) {
	    pj06_2.ind06 = 1;
	    con1 = pj06_2.fac * latc;
	    sinphi = sin(con1);
	    cosphi = cos(con1);
	    pj06_2.mcs = msfnz0_(&pj06_2.e06, &sinphi, &cosphi);
	    pj06_2.tcs = tsfnz0_(&pj06_2.e06, &con1, &sinphi);
	}

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	raddz0_(&pj06_2.lon006, sgna, degs, mins, secs, 1L);
	raddz0_(&latc, sgna + 1, &degs[1], &mins[1], &secs[1], 1L);
	data[1] = pj06_2.a06;
	data[2] = es;
	toggle_1.switch__[5] = *zone;
	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                            .  POLYCONIC  . */
/* ...................................................................... 
*/

    if (*isys == 7) {

L700:
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[6] != 0 && toggle_1.switch__[6] == *zone) {
	    return 0;
	}
	toggle_1.switch__[6] = 0;
	pj07_2.a07 = ellpz0_1.az;
	pj07_2.e07 = ellpz0_1.ez;
	pj07_2.es07 = ellpz0_1.esz;
	pj07_2.e007 = ellpz0_1.e0z;
	pj07_2.e107 = ellpz0_1.e1z;
	pj07_2.e207 = ellpz0_1.e2z;
	pj07_2.e307 = ellpz0_1.e3z;
	pj07_2.lon007 = pakrz0_(&data[5]);
	lat007 = pakrz0_(&data[6]);
	pj07_2.x007 = data[7];
	pj07_2.y007 = data[8];
	pj07_2.ml007 = mlfnz0_(&pj07_2.e007, &pj07_2.e107, &pj07_2.e207, &
		pj07_2.e307, &lat007);

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	raddz0_(&pj07_2.lon007, sgna, degs, mins, secs, 1L);
	raddz0_(&lat007, sgna + 1, &degs[1], &mins[1], &secs[1], 1L);
	data[1] = pj07_2.a07;
	data[2] = pj07_2.es07;
	toggle_1.switch__[6] = *zone;

/*     RETURN TO STATE PLANE INITIALIZATION IF NECESSARY */

	if (*isys == 2) {
	    goto L270;
	}

	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                        .  EQUIDISTANT CONIC  . */
/* ...................................................................... 
*/

    if (*isys == 8) {

	errmz0_1.ierror = 0;
	if (toggle_1.switch__[7] != 0 && toggle_1.switch__[7] == *zone) {
	    return 0;
	}
	toggle_1.switch__[7] = 0;
	pj08_2.a08 = ellpz0_1.az;
	e = ellpz0_1.ez;
	es = ellpz0_1.esz;
	pj08_2.e008 = ellpz0_1.e0z;
	pj08_2.e108 = ellpz0_1.e1z;
	pj08_2.e208 = ellpz0_1.e2z;
	pj08_2.e308 = ellpz0_1.e3z;
	pj17_2.lat1 = pakrz0_(&data[3]);
	lat2 = pakrz0_(&data[4]);
	if ((d__1 = pj17_2.lat1 + lat2, abs(d__1)) < epsln) {
	    errmz0_1.ierror = 81;
	    return 0;
	}
	pj08_2.lon008 = pakrz0_(&data[5]);
	lat0 = pakrz0_(&data[6]);
	pj08_2.x008 = data[7];
	pj08_2.y008 = data[8];
	sinphi = sin(pj17_2.lat1);
	cosphi = cos(pj17_2.lat1);
	ms1 = msfnz0_(&e, &sinphi, &cosphi);
	ml1 = mlfnz0_(&pj08_2.e008, &pj08_2.e108, &pj08_2.e208, &pj08_2.e308, 
		&pj17_2.lat1);
	ind = 0;
	if (data[9] != zero) {
	    ind = 1;
	    sinphi = sin(lat2);
	    cosphi = cos(lat2);
	    ms2 = msfnz0_(&e, &sinphi, &cosphi);
	    ml2 = mlfnz0_(&pj08_2.e008, &pj08_2.e108, &pj08_2.e208, &
		    pj08_2.e308, &lat2);
	    if ((d__1 = pj17_2.lat1 - lat2, abs(d__1)) >= epsln) {
		pj08_2.ns08 = (ms1 - ms2) / (ml2 - ml1);
	    } else {
		pj08_2.ns08 = sinphi;
	    }
	} else {
	    pj08_2.ns08 = sinphi;
	}
	pj08_2.gl = ml1 + ms1 / pj08_2.ns08;
	ml0 = mlfnz0_(&pj08_2.e008, &pj08_2.e108, &pj08_2.e208, &pj08_2.e308, 
		&lat0);
	pj08_2.rh008 = pj08_2.a08 * (pj08_2.gl - ml0);

/*    LIST RESULTS OF PARAMETER INITIALIZATION. */

	raddz0_(&pj17_2.lat1, sgna, degs, mins, secs, 1L);
	raddz0_(&lat2, sgna + 1, &degs[1], &mins[1], &secs[1], 1L);
	raddz0_(&pj08_2.lon008, sgna + 2, &degs[2], &mins[2], &secs[2], 1L);
	raddz0_(&lat0, sgna + 3, &degs[3], &mins[3], &secs[3], 1L);
	data[1] = pj08_2.a08;
	data[2] = es;
	toggle_1.switch__[7] = *zone;
	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                       .  TRANSVERSE MERCATOR  . */
/* ...................................................................... 
*/

    if (*isys == 9) {

L900:
	errmz0_1.ierror = 0;
	if (data[1] != 0. && data[1] != save) {
	    toggle_1.switch__[8] = 0;
	}
	if (toggle_1.switch__[8] != 0 && toggle_1.switch__[8] == *zone) {
	    return 0;
	}
	toggle_1.switch__[8] = 0;
	save = data[1];
	pj09_2.a09 = ellpz0_1.az;
	e09 = ellpz0_1.ez;
	pj09_2.es09 = ellpz0_1.esz;
	pj09_2.e009 = ellpz0_1.e0z;
	pj09_2.e109 = ellpz0_1.e1z;
	pj09_2.e209 = ellpz0_1.e2z;
	pj09_2.e309 = ellpz0_1.e3z;
	pj09_2.ks009 = data[3];
	pj09_2.lon009 = pakrz0_(&data[5]);
	pj09_2.lat009 = pakrz0_(&data[6]);
	pj09_2.x009 = data[7];
	pj09_2.y009 = data[8];
	pj09_2.ml009 = pj09_2.a09 * mlfnz0_(&pj09_2.e009, &pj09_2.e109, &
		pj09_2.e209, &pj09_2.e309, &pj09_2.lat009);
	pj09_2.ind09 = 1;
	pj09_2.esp = pj09_2.es09;
	if (e09 >= tol09) {
	    pj09_2.ind09 = 0;
	    pj09_2.esp = pj09_2.es09 / (one - pj09_2.es09);
	}

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	raddz0_(&pj09_2.lon009, sgna, degs, mins, secs, 1L);
	raddz0_(&pj09_2.lat009, sgna + 1, &degs[1], &mins[1], &secs[1], 1L);
	data[1] = pj09_2.a09;
	data[2] = pj09_2.es09;
	toggle_1.switch__[8] = *zone;

/*     RETURN TO UTM PROJECTION INITIALIZATION IF NECESSARY */

	if (*isys == 1) {
	    goto L170;
	}

/*     RETURN TO STATE PLANE INITIALIZATION IF NECESSARY */

	if (*isys == 2) {
	    goto L270;
	}

	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                          .  STEREOGRAPHIC  . */
/* ...................................................................... 
*/

    if (*isys == 10) {

	errmz0_1.ierror = 0;
	if (toggle_1.switch__[9] != 0 && toggle_1.switch__[9] == *zone) {
	    return 0;
	}
	toggle_1.switch__[9] = 0;
	pj10_2.a10 = sphrz0_1.azz;
	pj10_2.lon010 = pakrz0_(&data[5]);
	pj10_2.lat010 = pakrz0_(&data[6]);
	pj10_2.x010 = data[7];
	pj10_2.y010 = data[8];
	pj10_2.sinp10 = sin(pj10_2.lat010);
	pj10_2.cosp10 = cos(pj10_2.lat010);

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	raddz0_(&pj10_2.lon010, sgna, degs, mins, secs, 1L);
	raddz0_(&pj10_2.lat010, sgna + 1, &degs[1], &mins[1], &secs[1], 1L);
	data[1] = pj10_2.a10;
	toggle_1.switch__[9] = *zone;
	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                   .  LAMBERT AZIMUTHAL EQUAL-AREA  . */
/* ...................................................................... 
*/

    if (*isys == 11) {

	errmz0_1.ierror = 0;
	if (toggle_1.switch__[10] != 0 && toggle_1.switch__[10] == *zone) {
	    return 0;
	}
	toggle_1.switch__[10] = 0;
	pj11_2.a11 = sphrz0_1.azz;
	pj11_2.lon011 = pakrz0_(&data[5]);
	pj11_2.lat011 = pakrz0_(&data[6]);
	pj11_2.x011 = data[7];
	pj11_2.y011 = data[8];
	pj11_2.sinp11 = sin(pj11_2.lat011);
	pj11_2.cosp11 = cos(pj11_2.lat011);

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	raddz0_(&pj11_2.lon011, sgna, degs, mins, secs, 1L);
	raddz0_(&pj11_2.lat011, sgna + 1, &degs[1], &mins[1], &secs[1], 1L);
	data[1] = pj11_2.a11;
	toggle_1.switch__[10] = *zone;
	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                      .  AZIMUTHAL EQUIDISTANT  . */
/* ...................................................................... 
*/

    if (*isys == 12) {

	errmz0_1.ierror = 0;
	if (toggle_1.switch__[11] != 0 && toggle_1.switch__[11] == *zone) {
	    return 0;
	}
	toggle_1.switch__[11] = 0;
	pj12_2.a12 = sphrz0_1.azz;
	pj12_2.lon012 = pakrz0_(&data[5]);
	pj12_2.lat012 = pakrz0_(&data[6]);
	pj12_2.x012 = data[7];
	pj12_2.y012 = data[8];
	pj12_2.sinp12 = sin(pj12_2.lat012);
	pj12_2.cosp12 = cos(pj12_2.lat012);

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	raddz0_(&pj12_2.lon012, sgna, degs, mins, secs, 1L);
	raddz0_(&pj12_2.lat012, sgna + 1, &degs[1], &mins[1], &secs[1], 1L);
	data[1] = pj12_2.a12;
	toggle_1.switch__[11] = *zone;
	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                            .  GNOMONIC  . */
/* ...................................................................... 
*/

    if (*isys == 13) {

	errmz0_1.ierror = 0;
	if (toggle_1.switch__[12] != 0 && toggle_1.switch__[12] == *zone) {
	    return 0;
	}
	toggle_1.switch__[12] = 0;
	pj13_2.a13 = sphrz0_1.azz;
	pj13_2.lon013 = pakrz0_(&data[5]);
	pj13_2.lat013 = pakrz0_(&data[6]);
	pj13_2.x013 = data[7];
	pj13_2.y013 = data[8];
	pj13_2.sinp13 = sin(pj13_2.lat013);
	pj13_2.cosp13 = cos(pj13_2.lat013);

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	raddz0_(&pj13_2.lon013, sgna, degs, mins, secs, 1L);
	raddz0_(&pj13_2.lat013, sgna + 1, &degs[1], &mins[1], &secs[1], 1L);
	data[1] = pj13_2.a13;
	toggle_1.switch__[12] = *zone;
	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                          .  ORTHOGRAPHIC  . */
/* ...................................................................... 
*/

    if (*isys == 14) {

	errmz0_1.ierror = 0;
	if (toggle_1.switch__[13] != 0 && toggle_1.switch__[13] == *zone) {
	    return 0;
	}
	toggle_1.switch__[13] = 0;
	pj14_2.a14 = sphrz0_1.azz;
	pj14_2.lon014 = pakrz0_(&data[5]);
	pj14_2.lat014 = pakrz0_(&data[6]);
	pj14_2.x014 = data[7];
	pj14_2.y014 = data[8];
	pj14_2.sinp14 = sin(pj14_2.lat014);
	pj14_2.cosp14 = cos(pj14_2.lat014);

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	raddz0_(&pj14_2.lon014, sgna, degs, mins, secs, 1L);
	raddz0_(&pj14_2.lat014, sgna + 1, &degs[1], &mins[1], &secs[1], 1L);
	data[1] = pj14_2.a14;
	toggle_1.switch__[13] = *zone;
	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*             .  GENERAL VERTICAL NEAR-SIDE PERSPECTIVE  . */
/* ...................................................................... 
*/

    if (*isys == 15) {

	errmz0_1.ierror = 0;
	if (toggle_1.switch__[14] != 0 && toggle_1.switch__[14] == *zone) {
	    return 0;
	}
	toggle_1.switch__[14] = 0;
	pj15_2.a15 = sphrz0_1.azz;
	pj15_2.p = one + data[3] / pj15_2.a15;
	pj15_2.lon015 = pakrz0_(&data[5]);
	pj15_2.lat015 = pakrz0_(&data[6]);
	pj15_2.x015 = data[7];
	pj15_2.y015 = data[8];
	pj15_2.sinp15 = sin(pj15_2.lat015);
	pj15_2.cosp15 = cos(pj15_2.lat015);

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	raddz0_(&pj15_2.lon015, sgna, degs, mins, secs, 1L);
	raddz0_(&pj15_2.lat015, sgna + 1, &degs[1], &mins[1], &secs[1], 1L);
	data[1] = pj15_2.a15;
	toggle_1.switch__[14] = *zone;
	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                           .  SINUSOIDAL  . */
/* ...................................................................... 
*/

    if (*isys == 16) {

	errmz0_1.ierror = 0;
	if (toggle_1.switch__[15] != 0 && toggle_1.switch__[15] == *zone) {
	    return 0;
	}
	toggle_1.switch__[15] = 0;
	pj16_2.a16 = sphrz0_1.azz;
	pj16_2.lon016 = pakrz0_(&data[5]);
	pj16_2.x016 = data[7];
	pj16_2.y016 = data[8];

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	raddz0_(&pj16_2.lon016, sgna, degs, mins, secs, 1L);
	data[1] = pj16_2.a16;
	toggle_1.switch__[15] = *zone;
	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                         .  EQUIRECTANGULAR  . */
/* ...................................................................... 
*/

    if (*isys == 17) {

	errmz0_1.ierror = 0;
	if (toggle_1.switch__[16] != 0 && toggle_1.switch__[16] == *zone) {
	    return 0;
	}
	toggle_1.switch__[16] = 0;
	pj17_2.a17 = sphrz0_1.azz;
	pj17_2.lat1 = pakrz0_(&data[6]);
	pj17_2.lon017 = pakrz0_(&data[5]);
	pj17_2.x017 = data[7];
	pj17_2.y017 = data[8];

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	raddz0_(&pj17_2.lat1, sgna, degs, mins, secs, 1L);
	raddz0_(&pj17_2.lon017, sgna + 1, &degs[1], &mins[1], &secs[1], 1L);
	data[1] = pj17_2.a17;
	toggle_1.switch__[16] = *zone;
	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                       .  MILLER CYLINDRICAL  . */
/* ...................................................................... 
*/

    if (*isys == 18) {

	errmz0_1.ierror = 0;
	if (toggle_1.switch__[17] != 0 && toggle_1.switch__[17] == *zone) {
	    return 0;
	}
	toggle_1.switch__[17] = 0;
	pj18_2.a18 = sphrz0_1.azz;
	pj18_2.lon018 = pakrz0_(&data[5]);
	pj18_2.x018 = data[7];
	pj18_2.y018 = data[8];

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	raddz0_(&pj18_2.lon018, sgna, degs, mins, secs, 1L);
	data[1] = pj18_2.a18;
	toggle_1.switch__[17] = *zone;
	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                        .  VAN DER GRINTEN I  . */
/* ...................................................................... 
*/

    if (*isys == 19) {

	errmz0_1.ierror = 0;
	if (toggle_1.switch__[18] != 0 && toggle_1.switch__[18] == *zone) {
	    return 0;
	}
	toggle_1.switch__[18] = 0;
	pj19_2.a19 = sphrz0_1.azz;
	pj19_2.lon019 = pakrz0_(&data[5]);
	pj19_2.x019 = data[7];
	pj19_2.y019 = data[8];

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	raddz0_(&pj19_2.lon019, sgna, degs, mins, secs, 1L);
	data[1] = pj19_2.a19;
	toggle_1.switch__[18] = *zone;
	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                    .  OBLIQUE MERCATOR (HOTINE)  . */
/* ...................................................................... 
*/

    if (*isys == 20) {

L2000:
	errmz0_1.ierror = 0;
	if (toggle_1.switch__[19] != 0 && toggle_1.switch__[19] == *zone) {
	    return 0;
	}
	toggle_1.switch__[19] = 0;
	mode = 0;
	if (data[13] != zero) {
	    mode = 1;
	}
	a = ellpz0_1.az;
	pj20_2.e20 = ellpz0_1.ez;
	es = ellpz0_1.esz;
	ks0 = data[3];
	lat0 = pakrz0_(&data[6]);
	pj20_2.x020 = data[7];
	pj20_2.y020 = data[8];
	sinph0 = sin(lat0);
	cosph0 = cos(lat0);
	con = one - es * sinph0 * sinph0;
	com = sqrt(one - es);
/* Computing 4th power */
	d__1 = cosph0, d__1 *= d__1;
	pj20_2.bl = sqrt(one + es * (d__1 * d__1) / (one - es));
	pj20_2.al = a * pj20_2.bl * ks0 * com / con;
	if (abs(lat0) < epsln) {
	    ts0 = 1.;
	}
	if (abs(lat0) < epsln) {
	    d = 1.;
	}
	if (abs(lat0) < epsln) {
	    pj20_2.el = 1.;
	}
	if (abs(lat0) >= epsln) {
	    ts0 = tsfnz0_(&pj20_2.e20, &lat0, &sinph0);
	    con = sqrt(con);
	    d = pj20_2.bl * com / (cosph0 * con);
/* Computing MAX */
	    d__2 = d * d - one;
	    d__1 = sqrt((max(d__2,0.)));
	    f = d + d_sign(&d__1, &lat0);
	    pj20_2.el = f * pow_dd(&ts0, &pj20_2.bl);
	}
	if (mode != 0) {
	    alpha = pakrz0_(&data[4]);
	    lonc = pakrz0_(&data[5]);
	    g = half * (f - one / f);
	    d__1 = sin(alpha) / d;
	    _gamma = asinz0_(&d__1);
	    d__1 = g * tan(_gamma);
	    pj20_2.lon020 = lonc - asinz0_(&d__1) / pj20_2.bl;

/*     LIST INITIALIZATION PARAMETERS (CASE B). */

	    raddz0_(&alpha, sgna, degs, mins, secs, 1L);
	    raddz0_(&lonc, sgna + 1, &degs[1], &mins[1], &secs[1], 1L);
	    raddz0_(&lat0, sgna + 2, &degs[2], &mins[2], &secs[2], 1L);
	    con = abs(lat0);
	    if (con > epsln && (d__1 = con - halfpi, abs(d__1)) > epsln) {
		pj20_2.singam = sin(_gamma);
		pj20_2.cosgam = cos(_gamma);
		pj20_2.sinalf = sin(alpha);
		pj20_2.cosalf = cos(alpha);
		d__1 = pj20_2.al / pj20_2.bl * atan(sqrt(d * d - one) / 
			pj20_2.cosalf);
		pj20_2.u0 = d_sign(&d__1, &lat0);
		data[1] = a;
		data[2] = es;
		toggle_1.switch__[19] = *zone;

/*     RETURN TO STATE PLANE INITIALIZATION IF NECESSARY */

		if (*isys == 2) {
		    goto L270;
		}

		return 0;
	    } else {
		errmz0_1.ierror = 201;
		return 0;
	    }
	}
	lon1 = pakrz0_(&data[9]);
	pj17_2.lat1 = pakrz0_(&data[10]);
	lon2 = pakrz0_(&data[11]);
	lat2 = pakrz0_(&data[12]);
	sinphi = sin(pj17_2.lat1);
	ts1 = tsfnz0_(&pj20_2.e20, &pj17_2.lat1, &sinphi);
	sinphi = sin(lat2);
	ts2 = tsfnz0_(&pj20_2.e20, &lat2, &sinphi);
	h = pow_dd(&ts1, &pj20_2.bl);
	l = pow_dd(&ts2, &pj20_2.bl);
	f = pj20_2.el / h;
	g = half * (f - one / f);
	j = (pj20_2.el * pj20_2.el - l * h) / (pj20_2.el * pj20_2.el + l * h);

	pj15_2.p = (l - h) / (l + h);
	raddz0_(&lon2, sgna + 2, &degs[2], &mins[2], &secs[2], 1L);
	dlon = lon1 - lon2;
	if (dlon < -pi) {
	    lon2 -= pi * 2.;
	}
	if (dlon > pi) {
	    lon2 += pi * 2.;
	}
	dlon = lon1 - lon2;
	pj20_2.lon020 = half * (lon1 + lon2) - atan(j * tan(half * pj20_2.bl *
		 dlon) / pj15_2.p) / pj20_2.bl;
	d__1 = lon1 - pj20_2.lon020;
	dlon = adjlz0_(&d__1);
	_gamma = atan(sin(pj20_2.bl * dlon) / g);
	d__1 = d * sin(_gamma);
	alpha = asinz0_(&d__1);
	raddz0_(&lon1, sgna, degs, mins, secs, 1L);
	raddz0_(&pj17_2.lat1, sgna + 1, &degs[1], &mins[1], &secs[1], 1L);
	raddz0_(&lat2, sgna + 3, &degs[3], &mins[3], &secs[3], 1L);
	raddz0_(&lat0, sgna + 4, &degs[4], &mins[4], &secs[4], 1L);
	if ((d__1 = pj17_2.lat1 - lat2, abs(d__1)) <= epsln) {
	    errmz0_1.ierror = 202;
	    return 0;
	} else {
	    con = abs(pj17_2.lat1);
	}
	if (con <= epsln || (d__1 = con - halfpi, abs(d__1)) <= epsln) {
	    errmz0_1.ierror = 202;
	    return 0;
	} else {
	    if ((d__1 = abs(lat0) - halfpi, abs(d__1)) <= epsln) {
		errmz0_1.ierror = 202;
		return 0;
	    }
	}
	pj20_2.singam = sin(_gamma);
	pj20_2.cosgam = cos(_gamma);
	pj20_2.sinalf = sin(alpha);
	pj20_2.cosalf = cos(alpha);
	d__1 = pj20_2.al / pj20_2.bl * atan(sqrt(d * d - one) / pj20_2.cosalf)
		;
	pj20_2.u0 = d_sign(&d__1, &lat0);
	data[1] = a;
	data[2] = es;
	toggle_1.switch__[19] = *zone;
	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                       .       ROBINSON       . */
/* ...................................................................... 
*/

    if (*isys == 21) {

	errmz0_1.ierror = 0;
	if (toggle_1.switch__[20] != 0 && toggle_1.switch__[20] == *zone) {
	    return 0;
	}
	toggle_1.switch__[20] = 0;
	pj21_2.a21 = sphrz0_1.azz;
	pj21_2.lon021 = pakrz0_(&data[5]);
	pj21_2.x021 = data[7];
	pj21_2.y021 = data[8];
	pj21_2.pr[0] = -.062;
	pj21_2.xlr[0] = .9986;
	pj21_2.pr[1] = 0.;
	pj21_2.xlr[1] = 1.;
	pj21_2.pr[2] = .062;
	pj21_2.xlr[2] = .9986;
	pj21_2.pr[3] = .124;
	pj21_2.xlr[3] = .9954;
	pj21_2.pr[4] = .186;
	pj21_2.xlr[4] = .99;
	pj21_2.pr[5] = .248;
	pj21_2.xlr[5] = .9822;
	pj21_2.pr[6] = .31;
	pj21_2.xlr[6] = .973;
	pj21_2.pr[7] = .372;
	pj21_2.xlr[7] = .96;
	pj21_2.pr[8] = .434;
	pj21_2.xlr[8] = .9427;
	pj21_2.pr[9] = .4958;
	pj21_2.xlr[9] = .9216;
	pj21_2.pr[10] = .5571;
	pj21_2.xlr[10] = .8962;
	pj21_2.pr[11] = .6176;
	pj21_2.xlr[11] = .8679;
	pj21_2.pr[12] = .6769;
	pj21_2.xlr[12] = .835;
	pj21_2.pr[13] = .7346;
	pj21_2.xlr[13] = .7986;
	pj21_2.pr[14] = .7903;
	pj21_2.xlr[14] = .7597;
	pj21_2.pr[15] = .8435;
	pj21_2.xlr[15] = .7186;
	pj21_2.pr[16] = .8936;
	pj21_2.xlr[16] = .6732;
	pj21_2.pr[17] = .9394;
	pj21_2.xlr[17] = .6213;
	pj21_2.pr[18] = .9761;
	pj21_2.xlr[18] = .5722;
	pj21_2.pr[19] = 1.;
	pj21_2.xlr[19] = .5322;
	for (i = 1; i <= 20; ++i) {
/* L2110: */
	    pj21_2.xlr[i - 1] *= .9858;
	}

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	raddz0_(&pj21_2.lon021, sgna, degs, mins, secs, 1L);
	data[1] = pj21_2.a21;
	toggle_1.switch__[20] = *zone;
	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*                      .  SPACE OBLIQUE MERCATOR  . */
/* ...................................................................... 
*/

    if (*isys == 22) {

	errmz0_1.ierror = 0;
	if (toggle_1.switch__[21] != 0 && toggle_1.switch__[21] == *zone) {
	    return 0;
	}
	toggle_1.switch__[21] = 0;
	pj22_2.a22 = ellpz0_1.az;
	e = ellpz0_1.ez;
	norm_2.es22 = ellpz0_1.esz;
	pj22_2.x022 = data[7];
	pj22_2.y022 = data[8];
	pj22_2.land = (integer) (data[3] + tol);
	pj22_2.path = (integer) (data[4] + tol);

/*        CHECK IF LANDSAT NUMBER IS WITHIN RANGE 1 - 5 */

	if (pj22_2.land > 0 && pj22_2.land <= 5) {
	    if (pj22_2.land <= 3) {
		limit = 251;
	    }
	    if (pj22_2.land >= 4) {
		limit = 233;
	    }
	} else {
	    errmz0_1.ierror = 221;
	    return 0;
	}

/*        CHECK IF PATH NUMBER IS WITHIN RANGE 1 - 251 FOR LANDSATS 1 
- 3 */
/*        OR RANGE 1 - 233 FOR LANDSATS 4 - 5 */

	if (pj22_2.path <= 0 || pj22_2.path > limit) {
	    errmz0_1.ierror = 221;
	    return 0;
	}
	p1 = 1440.;
	if (pj22_2.land <= 3) {
	    p2 = 103.2669323;
	    alf = dg1 * 99.092;
	} else {
	    p2 = 98.8841202;
	    alf = dg1 * 98.2;
	}
	norm_2.sa = sin(alf);
	norm_2.ca = cos(alf);
	if (abs(norm_2.ca) < 1e-9) {
	    norm_2.ca = 1e-9;
	}
	esc = norm_2.es22 * norm_2.ca * norm_2.ca;
	ess = norm_2.es22 * norm_2.sa * norm_2.sa;
	d__1 = (one - esc) / (one - norm_2.es22);
	norm_2.w = pow_dd(&d__1, &two) - one;
	norm_2.q = ess / (one - norm_2.es22);
	d__1 = one - norm_2.es22;
	norm_2.t = ess * (two - norm_2.es22) / pow_dd(&d__1, &two);
	norm_2.u = esc / (one - norm_2.es22);
/* Computing 3rd power */
	d__1 = one - norm_2.es22, d__2 = d__1;
	norm_2.xj = d__2 * (d__1 * d__1);
	norm_2.p22 = p2 / p1;

/*        COMPUTE FOURIER COEFFICIENTS.  LAM IS CURRENT VALUE OF */
/*        LAMBDA DOUBLE-PRIME. */

	lam = 0.;
	seraz0_(&fb, &fa2, &fa4, &fc1, &fc3, &lam);
	suma2 = fa2;
	suma4 = fa4;
	sumb = fb;
	sumc1 = fc1;
	sumc3 = fc3;
	for (i = 9; i <= 81; i += 18) {
	    lam = (doublereal) i;
	    seraz0_(&fb, &fa2, &fa4, &fc1, &fc3, &lam);
	    suma2 += fa2 * 4.;
	    suma4 += fa4 * 4.;
	    sumb += fb * 4.;
	    sumc1 += fc1 * 4.;
	    sumc3 += fc3 * 4.;
/* L2210: */
	}
	for (i = 18; i <= 72; i += 18) {
	    lam = (doublereal) i;
	    seraz0_(&fb, &fa2, &fa4, &fc1, &fc3, &lam);
	    suma2 += two * fa2;
	    suma4 += two * fa4;
	    sumb += two * fb;
	    sumc1 += two * fc1;
	    sumc3 += two * fc3;
/* L2220: */
	}
	lam = 90.;
	seraz0_(&fb, &fa2, &fa4, &fc1, &fc3, &lam);
	suma2 += fa2;
	suma4 += fa4;
	sumb += fb;
	sumc1 += fc1;
	sumc3 += fc3;

/*        THESE ARE THE VALUES OF FOURIER CONSTANTS. */

	pj22_2.a2 = suma2 / 30.;
	pj22_2.a4 = suma4 / 60.;
	pj22_2.b = sumb / 30.;
	pj22_2.c1 = sumc1 / 15.;
	pj22_2.c3 = sumc3 / 45.;

/*        LIST RESULTS OF PARAMETER INITIALIZATION. */

	data[1] = pj22_2.a22;
	data[2] = norm_2.es22;
	toggle_1.switch__[21] = *zone;
	return 0;
    }

/* ...................................................................... 
*/
/*             .  INITIALIZATION OF PROJECTION PARAMETERS  . */

/*          .  MODIFIED-STEREOGRAPHIC CONFORMAL (FOR ALASKA)  . */
/* ...................................................................... 
*/

    if (*isys == 23) {

	errmz0_1.ierror = 0;
	if (toggle_1.switch__[22] != 0 && toggle_1.switch__[22] == *zone) {
	    return 0;
	}
	toggle_1.switch__[22] = 0;
	pj23_2.a23 = ellpz0_1.az;
	ec2 = .006768657997291094;
	pj23_2.ec = sqrt(ec2);
	pj23_2.n = 6;
	pj23_2.lon023 = dg1 * -152.;
	pj23_2.lat023 = dg1 * 64.;
	pj23_2.x023 = data[7];
	pj23_2.y023 = data[8];
	pj23_2.acoef[0] = .9945303;
	pj23_2.acoef[1] = .0052083;
	pj23_2.acoef[2] = .0072721;
	pj23_2.acoef[3] = -.0151089;
	pj23_2.acoef[4] = .0642675;
	pj23_2.acoef[5] = .3582802;
	pj23_2.bcoef[0] = 0.;
	pj23_2.bcoef[1] = -.0027404;
	pj23_2.bcoef[2] = .0048181;
	pj23_2.bcoef[3] = -.1932526;
	pj23_2.bcoef[4] = -.1381226;
	pj23_2.bcoef[5] = -.2884586;
	esphi = pj23_2.ec * sin(pj23_2.lat023);
	d__1 = (one - esphi) / (one + esphi);
	d__2 = pj23_2.ec / two;
	chio = two * atan(tan((halfpi + pj23_2.lat023) / two) * pow_dd(&d__1, 
		&d__2)) - halfpi;
	pj23_2.schio = sin(chio);
	pj23_2.cchio = cos(chio);

/*     LIST RESULTS OF PARAMETER INITIALIZATION. */

	raddz0_(&pj23_2.lon023, sgna, degs, mins, secs, 1L);
	raddz0_(&pj23_2.lat023, sgna + 1, &degs[1], &mins[1], &secs[1], 1L);
	data[1] = pj23_2.a23;
	toggle_1.switch__[22] = *zone;
	return 0;
    }

    return 0;

/*     INITIALIZATION OF PROJECTION COMPLETED */

/*###1234 [f77] Warning on line 1234 of pjinit.f local variable pname neve
r used%%%*/
} /* pjinit_ */

/*                   QSFNZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

doublereal qsfnz0_(eccent, sinphi, cosphi)
doublereal *eccent, *sinphi, *cosphi;
{
    /* Initialized data */

    static doublereal half = .5;
    static doublereal one = 1.;
    static doublereal two = 2.;
    static doublereal epsln = 1e-7;

    /* System generated locals */
    doublereal ret_val;

    /* Builtin functions */
    double log();

    /* Local variables */
    static doublereal con;


/* FUNCTION TO COMPUTE CONSTANT (SMALL Q). */


    if (*eccent < epsln) {
	goto L20;
    }
    con = *eccent * *sinphi;
    ret_val = (one - *eccent * *eccent) * (*sinphi / (one - con * con) - half 
	    / *eccent * log((one - con) / (one + con)));
    return ret_val;

L20:
    ret_val = two * *sinphi;
    return ret_val;
} /* qsfnz0_ */

/*                   RADDZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

/* Subroutine */ int raddz0_(rad, sgna, degs, mins, secs, sgna_len)
doublereal *rad;
char *sgna;
integer *degs, *mins;
real *secs;
ftnlen sgna_len;
{
    /* Initialized data */

    static doublereal radsec = 206264.806247;
    static real zero = (float)0.;
    static char blank[1+1] = " ";
    static char neg[1+1] = "-";

    static integer isec, intg;
    static doublereal con;


/* SUBROUTINE TO CONVERT ANGLE FROM RADIANS TO SIGNED DMS */
/* SGNA : SIGN OF ANGLE */
/* DEGS : DEGREES PORTION OF ANGLE */
/* MINS : MINUTES PORTION OF ANGLE */
/* SECS : SECONDS PORTION OF ANGLE */


/* DETERMINE THE SIGN OF THE ANGLE. */

    *sgna = blank[0];
    if (*rad < zero) {
	*sgna = neg[0];
    }

/* CONVERT THE ANGLE TO THOUSANDTH OF SECONDS. */

    con = abs(*rad) * radsec;
    isec = (integer) con;

/* COMPUTE DEGREES PART OF THE ANGLE. */

    intg = isec / 3600;
    *degs = intg;
    isec = intg * 3600;
    con -= (doublereal) isec;
    isec = (integer) con;

/* COMPUTE MINUTES PART OF THE ANGLE. */

    *mins = isec / 60;
    isec = *mins * 60;
    con -= (doublereal) isec;

/* COMPUTE SECONDS PART OF THE ANGLE. */

    *secs = con;

/*     INCREASE MINS IF SECS CLOSE TO 60.000 */

    if (*secs < 59.9995) {
	return 0;
    }
    ++(*mins);
    *secs = (float)0.;

/*     INCREASE DEGS IF MINS EQUAL 60 */

    if (*mins <= 59) {
	return 0;
    }
    *mins = 0;
    ++(*degs);

    return 0;
} /* raddz0_ */

/*                   SERAZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

/* Subroutine */ int seraz0_(fb, fa2, fa4, fc1, fc3, lam)
doublereal *fb, *fa2, *fa4, *fc1, *fc3, *lam;
{
    /* Initialized data */

    static doublereal dg1 = .01745329252;
    static doublereal one = 1.;
    static doublereal two = 2.;

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double sin(), cos(), sqrt(), pow_dd();

    /* Local variables */
    static doublereal sdsq, h, s, fc, sd, sq;


/*     COMPUTES INTEGRAL FUNCTION OF TRANSFORMED LONG. FOR FOURIER */
/*     CONSTANTS A2, A4, B, C1, AND C3. */
/*     LAM IS INTEGRAL VALUE OF TRANSFORMED LONG. */

    *lam *= dg1;
    sd = sin(*lam);
    sdsq = sd * sd;
    s = norm_1.p22 * norm_1.sa * cos(*lam) * sqrt((one + norm_1.t * sdsq) / ((
	    one + norm_1.w * sdsq) * (one + norm_1.q * sdsq)));
    d__1 = one + norm_1.q * sdsq;
    h = sqrt((one + norm_1.q * sdsq) / (one + norm_1.w * sdsq)) * ((one + 
	    norm_1.w * sdsq) / pow_dd(&d__1, &two) - norm_1.p22 * norm_1.ca);
    sq = sqrt(norm_1.xj * norm_1.xj + s * s);
    *fb = (h * norm_1.xj - s * s) / sq;
    *fa2 = *fb * cos(two * *lam);
    *fa4 = *fb * cos(*lam * 4.);
    fc = s * (h + norm_1.xj) / sq;
    *fc1 = fc * cos(*lam);
    *fc3 = fc * cos(*lam * 3.);
    return 0;
} /* seraz0_ */

/*                   SPHDZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

/* Subroutine */ int sphdz0_(isph, parm)
integer *isph;
doublereal *parm;
{
    /* Initialized data */

    static doublereal zero = 0.;
    static doublereal one = 1.;
    static doublereal axis[20] = { 6378206.4,6378249.145,6377397.155,
	    6378157.5,6378388.,6378135.,6377276.3452,6378145.,6378137.,
	    6377563.396,6377304.063,6377340.189,6378137.,6378155.,6378160.,
	    6378245.,6378270.,6378166.,6378150.,6370997. };
    static doublereal bxis[20] = { 6356583.8,6356514.86955,6356078.96284,
	    6356772.2,6356911.94613,6356750.519915,6356075.4133,
	    6356759.769356,6356752.31414,6356256.91,6356103.039,6356034.448,
	    6356752.314245,6356773.3205,6356774.719,6356863.0188,
	    6356794.343479,6356784.283666,6356768.337303,6370997. };

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double sqrt();

    /* Local variables */
    static integer jsph;
    static doublereal a, b;
    extern doublereal e0fnz0_(), e1fnz0_(), e2fnz0_(), e3fnz0_(), e4fnz0_();
    static doublereal es;


/*     SUBROUTINE TO COMPUTE SPHEROID PARAMETERS */

/*     ISPH IS THE SPHEROID CODE FROM THE FOLLOWING LIST: */
/*     0 = CLARKE 1866           1 = CLARKE 1880 */
/*     2 = BESSEL                3 = NEW INTERNATIONAL 1967 */
/*     4 = INTERNATIONAL 1909    5 = WGS 72 */
/*     6 = EVEREST               7 = WGS 66 */
/*     8 = GRS 1980              9 = AIRY */
/*    10 = MODIFIED EVEREST     11 = MODIFIED AIRY */
/*    12 = WGS 84               13 = SOUTHEAST ASIA */
/*    14 = AUSTRALIAN NATIONAL  15 = KRASSOVSKY */
/*    16 = HOUGH                17 = MERCURY 1960 */
/*    18 = MODIFIED MERC 1968   19 = SPHERE OF RADIUS 6370997 M */

/*    PARM IS ARRAY OF PROJECTION PARAMETERS: */
/*       PARM(1) IS THE SEMI-MAJOR AXIS */
/*       PARM(2) IS THE ECCENTRICITY SQUARED */

/*     IF ISPH IS NEGATIVE, USER SPECIFIED PROJECTION PARAMETERS ARE TO */

/*     DEFINE THE RADIUS OF SPHERE OR ELLIPSOID CONSTANTS AS APPROPRIATE 
*/

/*     IF ISPH = 0 , THE DEFAULT IS RESET TO CLARKE 1866 */

/* ****                                                             ***** 
*/



    /* Parameter adjustments */
    --parm;

    /* Function Body */



    if (*isph >= 0) {
	goto L5;
    }

/*     INITIALIZE USER SPECIFIED SPHERE AND ELLIPSOID PARAMETERS */

    sphrz0_1.azz = zero;
    ellpz0_1.az = zero;
    ellpz0_1.ez = zero;
    ellpz0_1.esz = zero;
    ellpz0_1.e0z = zero;
    ellpz0_1.e1z = zero;
    ellpz0_1.e2z = zero;
    ellpz0_1.e3z = zero;
    ellpz0_1.e4z = zero;

/*     FETCH FIRST TWO USER SPECIFIED PROJECTION PARAMETERS */

    a = abs(parm[1]);
    b = abs(parm[2]);
    if (a > zero && b > zero) {
	goto L13;
    }
    if (a > zero && b <= zero) {
	goto L12;
    }
    if (a <= zero && b > zero) {
	goto L11;
    }

/*     DEFAULT NORMAL SPHERE AND CLARKE 1866 ELLIPSOID */

    jsph = 1;
    goto L10;

/*     DEFAULT CLARKE 1866 ELLIPSOID */

L11:
    a = axis[0];
    b = bxis[0];
    goto L14;

/*     USER SPECIFIED RADIUS OF SPHERE */

L12:
    sphrz0_1.azz = a;
    goto L15;

/*     USER SPECIFIED SEMI-MAJOR AND SEMI-MINOR AXES OF ELLIPSOID */

L13:
    if (b <= one) {
	goto L15;
    }
L14:
/* Computing 2nd power */
    d__1 = b / a;
    es = one - d__1 * d__1;
    goto L16;

/*     USER SPECIFIED SEMI-MAJOR AXIS AND ECCENTRICITY SQUARED */

L15:
    es = b;
L16:
    ellpz0_1.az = a;
    ellpz0_1.esz = es;
    ellpz0_1.ez = sqrt(es);
    ellpz0_1.e0z = e0fnz0_(&es);
    ellpz0_1.e1z = e1fnz0_(&es);
    ellpz0_1.e2z = e2fnz0_(&es);
    ellpz0_1.e3z = e3fnz0_(&es);
    ellpz0_1.e4z = e4fnz0_(&ellpz0_1.ez);
    parm[1] = a;
    parm[2] = es;
    return 0;

/*     CHECK FOR VALID SPHEROID SELECTION */

L5:
    if (parm[1] != zero && projz0_1.iproj != 1) {
	return 0;
    }
    jsph = abs(*isph) + 1;
    if (jsph <= 20) {
	goto L10;
    }
    errmz0_1.ierror = 999;
    *isph = 0;
    jsph = 1;

/*     RETRIEVE A AND B AXES FOR SELECTED SPHEROID */

L10:
    a = axis[jsph - 1];
    b = bxis[jsph - 1];
/* Computing 2nd power */
    d__1 = b / a;
    es = one - d__1 * d__1;

/*     SET COMMON BLOCK PARAMETERS FOR SELECTED SPHEROID */

    sphrz0_1.azz = 6370997.;
    ellpz0_1.ez = sqrt(es);
    ellpz0_1.e0z = e0fnz0_(&es);
    ellpz0_1.e1z = e1fnz0_(&es);
    ellpz0_1.e2z = e2fnz0_(&es);
    ellpz0_1.e3z = e3fnz0_(&es);
    ellpz0_1.e4z = e4fnz0_(&ellpz0_1.ez);
    ellpz0_1.az = a;
    ellpz0_1.esz = es;
    if (es == zero) {
	sphrz0_1.azz = a;
    }

    parm[1] = a;
    parm[2] = es;
    return 0;
} /* sphdz0_ */

/*                   TSFNZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

doublereal tsfnz0_(eccent, phi, sinphi)
doublereal *eccent, *phi, *sinphi;
{
    /* Initialized data */

    static doublereal half = .5;
    static doublereal one = 1.;
    static doublereal halfpi = 1.5707963267948966;

    /* System generated locals */
    doublereal ret_val, d__1;

    /* Builtin functions */
    double pow_dd(), tan();

    /* Local variables */
    static doublereal com, con;


/* FUNCTION TO COMPUTE CONSTANT (SMALL T). */


    con = *eccent * *sinphi;
    com = half * *eccent;
    d__1 = (one - con) / (one + con);
    con = pow_dd(&d__1, &com);
    ret_val = tan(half * (halfpi - *phi)) / con;

    return ret_val;
} /* tsfnz0_ */

/*                   UNTFZ0 */
/* ********************************************************************** */
/* ** GENERAL CARTOGRAPHIC TRANSFORMATION PACKAGE (GCTP) VERSION 2.0.0 ** */
/* ** U. S. GEOLOGICAL SURVEY - SNYDER, ELASSAL, AND LINCK    06/18/92 ** */
/* ********************************************************************** */

/* Subroutine */ int untfz0_(inunit, iounit, factor, iflg)
integer *inunit, *iounit;
doublereal *factor;
integer *iflg;
{
    /* Initialized data */

    static doublereal factrs[36]	/* was [6][6] */ = { 1.,0.,0.,
	    206264.8062470963,57.29577951308231,0.,0.,1.,.3048006096012192,0.,
	    0.,1.000002000004,0.,3.280833333333333,1.,0.,0.,3.280839895013124,
	    4.84813681109536e-6,0.,0.,1.,2.777777777777778e-4,0.,
	    .0174532925199433,0.,0.,3600.,1.,0.,0.,.999998,.3048,0.,0.,1. };


/* SUBROUTINE TO DETERMINE CONVERGENCE FACTOR BETWEEN TWO LINEAL UNITS */

/* * INPUT ........ */
/* * INUNIT * UNIT CODE OF SOURCE. */
/* * IOUNIT * UNIT CODE OF TARGET. */

/* * OUTPUT ....... */
/* * FACTOR * CONVERGENCE FACTOR FROM SOURCE TO TARGET. */
/* * IFLG   * RETURN FLAG .EQ. 0 , NORMAL RETURN. */
/*            RETURN FLAG .NE. 0 , ABNORMAL RETURN. */


    if (*inunit >= 0 && *inunit < 6 && *iounit >= 0 && *iounit < 6) {
	*factor = factrs[*iounit + 1 + (*inunit + 1) * 6 - 7];
	if (*factor != 0.) {
	    *iflg = 0;
	    return 0;
	} else {
	    *iflg = 12;
	    return 0;
	}
    } else {
	*iflg = 11;
	return 0;
    }

} /* untfz0_ */

#ifdef KR_headers
double d_mod(x,y) doublereal *x, *y;
#else
double d_mod(doublereal *x, doublereal *y)
#endif
{
#ifdef IEEE_drem
	double xa, ya, z;
	if ((ya = *y) < 0.)
		ya = -ya;
	z = drem(xa = *x, ya);
	if (xa > 0) {
		if (z < 0)
			z += ya;
		}
	else if (z > 0)
		z -= ya;
	return z;
#else
	double quotient;
	if( (quotient = *x / *y) >= 0)
		quotient = floor(quotient);
	else
		quotient = -floor(-quotient);
	return(*x - (*y) * quotient );
#endif
}

#ifdef KR_headers
double d_sign(a,b) doublereal *a, *b;
#else
double d_sign(doublereal *a, doublereal *b)
#endif
{
double x;
x = (*a >= 0 ? *a : - *a);
return( *b >= 0 ? x : -x);
}

#ifdef KR_headers
double pow_dd(ap, bp) doublereal *ap, *bp;
#else
double pow_dd(doublereal *ap, doublereal *bp)
#endif
{
return(pow(*ap, *bp) );
}

/* assign strings:  a = b */

#ifdef KR_headers
s_copy(a, b, la, lb) char *a, *b; ftnlen la, lb;
#else
s_copy(char *a, char *b, ftnlen la, ftnlen lb)
#endif
{
char *aend, *bend;

aend = a + la;

if(la <= lb)
	while(a < aend)
		*a++ = *b++;

else
	{
	bend = b + lb;
	while(b < bend)
		*a++ = *b++;
	while(a < aend)
		*a++ = ' ';
	}
}