Application: gvSIG desktop » gvSIG 3D

<?xml version="1.0"?>

<!DOCTYPE CodeBlocks_project_file>

<CodeBlocks_project_file>

	<FileVersion major="1" minor="1"/>

	<Project>

		<Option title="crsjniproj"/>

		<Option makefile="Makefile"/>

		<Option makefile_is_custom="0"/>

		<Option default_target="-1"/>

		<Option compiler="0"/>

		<Build>

			<Target title="default">

				<Option output="dist\crsjniproj.dll"/>

				<Option working_dir="."/>

				<Option object_output="dist"/>

				<Option deps_output=".deps"/>

				<Option type="3"/>

				<Option compiler="0"/>

				<Option createDefFile="1"/>

			</Target>

		</Build>

		<Compiler>

			<Add option="-DBUILDING_DLL=1"/>

			<Add directory="include"/>

			<Add directory="C:\j2sdk1.4.2_11\include\win32"/>

			<Add directory="C:\j2sdk1.4.2_11\include"/>

		</Compiler>

		<Unit filename="src\crsJniProj.c">

			<Option compilerVar="CC"/>

			<Option target="default"/>

		</Unit>

	</Project>

</CodeBlocks_project_file>

/* arc sin, cosine, tan2 and sqrt that will NOT fail */

#ifndef lint

static const char SCCSID[]="@(#)aasincos.c	4.6	93/12/12	GIE	REL";

#endif

#include <projects.h>

#define ONE_TOL	 1.00000000000001

#define TOL	0.000000001

#define ATOL 1e-50

	double

aasin(double v) {

	double av;

	if ((av = fabs(v)) >= 1.) {

		if (av > ONE_TOL)

			pj_errno = -19;

		return (v < 0. ? -HALFPI : HALFPI);

	}

	return asin(v);

}

	double

aacos(double v) {

	double av;

	if ((av = fabs(v)) >= 1.) {

		if (av > ONE_TOL)

			pj_errno = -19;

		return (v < 0. ? PI : 0.);

	}

	return acos(v);

}

	double

asqrt(double v) { return ((v <= 0) ? 0. : sqrt(v)); }

	double

aatan2(double n, double d) {

	return ((fabs(n) < ATOL && fabs(d) < ATOL) ? 0. : atan2(n,d));

}

/* dervative of (*P->fwd) projection */

#ifndef lint

static const char SCCSID[]="@(#)pj_deriv.c	4.4	93/06/12	GIE	REL";

#endif

#define PJ_LIB__

#include "projects.h"

	int

pj_deriv(LP lp, double h, PJ *P, struct DERIVS *der) {

	XY t;

	lp.lam += h;

	lp.phi += h;

	if (fabs(lp.phi) > HALFPI) return 1;

	h += h;

	t = (*P->fwd)(lp, P);

	if (t.x == HUGE_VAL) return 1;

	der->x_l = t.x; der->y_p = t.y; der->x_p = -t.x; der->y_l = -t.y;

	lp.phi -= h;

	if (fabs(lp.phi) > HALFPI) return 1;

	t = (*P->fwd)(lp, P);

	if (t.x == HUGE_VAL) return 1;

	der->x_l += t.x; der->y_p -= t.y; der->x_p += t.x; der->y_l -= t.y;

	lp.lam -= h;

	t = (*P->fwd)(lp, P);

	if (t.x == HUGE_VAL) return 1;

	der->x_l -= t.x; der->y_p -= t.y; der->x_p += t.x; der->y_l += t.y;

	lp.phi += h;

	t = (*P->fwd)(lp, P);

	if (t.x == HUGE_VAL) return 1;

	der->x_l -= t.x; der->y_p += t.y; der->x_p -= t.x; der->y_l += t.y;

	der->x_l /= (h += h);

	der->y_p /= h;

	der->x_p /= h;

	der->y_l /= h;

	return 0;

}

#ifndef lint

static const char SCCSID[]="@(#)PJ_urmfps.c	4.1	94/02/15	GIE	REL";

#endif

#define PROJ_PARMS__ \

	double	n, C_y;

#define PJ_LIB__

#include	<projects.h>

PROJ_HEAD(urmfps, "Urmaev Flat-Polar Sinusoidal") "\n\tPCyl, Sph.\n\tn=";

PROJ_HEAD(wag1, "Wagner I (Kavraisky VI)") "\n\tPCyl, Sph.";

#define C_x 0.8773826753

#define Cy 1.139753528477

FORWARD(s_forward); /* sphere */

	lp.phi = aasin(P->n * sin(lp.phi));

	xy.x = C_x * lp.lam * cos(lp.phi);

	xy.y = P->C_y * lp.phi;

	return (xy);

}

INVERSE(s_inverse); /* sphere */

	xy.y /= P->C_y;

	lp.phi = aasin(sin(xy.y) / P->n);

	lp.lam = xy.x / (C_x * cos(xy.y));

	return (lp);

}

FREEUP; if (P) pj_dalloc(P); }

	static PJ *

setup(PJ *P) {

	P->C_y = Cy / P->n;

	P->es = 0.;

	P->inv = s_inverse;

	P->fwd = s_forward;

	return P;

}

ENTRY0(urmfps)

	if (pj_param(P->params, "tn").i) {

		P->n = pj_param(P->params, "dn").f;

		if (P->n <= 0. || P->n > 1.)

			E_ERROR(-40)

	} else

		E_ERROR(-40)

ENDENTRY(setup(P))

ENTRY0(wag1)

	P->n = 0.8660254037844386467637231707;

ENDENTRY(setup(P))

/* print projection's list of parameters */

#ifndef lint

static const char SCCSID[]="@(#)pj_pr_list.c	4.6   94/03/19 GIE REL";

#endif

#include <projects.h>

#include <stdio.h>

#include <string.h>

#define LINE_LEN 72

	static int

pr_list(PJ *P, int not_used) {

	paralist *t;

	int l, n = 1, flag = 0;

	(void)putchar('#');

	for (t = P->params; t; t = t->next)

		if ((!not_used && t->used) || (not_used && !t->used)) {

			l = strlen(t->param) + 1;

			if (n + l > LINE_LEN) {

				(void)fputs("\n#", stdout);

				n = 2;

			}

			(void)putchar(' ');

			if (*(t->param) != '+')

				(void)putchar('+');

			(void)fputs(t->param, stdout);

			n += l;

		} else

			flag = 1;

	if (n > 1)

		(void)putchar('\n');

	return flag;

}

	void /* print link list of projection parameters */

pj_pr_list(PJ *P) {

	char const *s;

	(void)putchar('#');

	for (s = P->descr; *s ; ++s) {

		(void)putchar(*s);

		if (*s == '\n')

			(void)putchar('#');

	}

	(void)putchar('\n');

	if (pr_list(P, 0)) {

		(void)fputs("#--- following specified but NOT used\n", stdout);

		(void)pr_list(P, 1);

	}

}

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

/*                             pj_get_def()                             */

/*                                                                      */

/*      Returns the PROJ.4 command string that would produce this       */

/*      definition expanded as much as possible.  For instance,         */

/*      +init= calls and +datum= defintions would be expanded.          */

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

char *pj_get_def( PJ *P, int options )

{

    paralist *t;

    int l;

    char *definition;

    int  def_max = 10;

    definition = (char *) pj_malloc(def_max);

    definition[0] = '\0';

    for (t = P->params; t; t = t->next)

    {

        /* skip unused parameters ... mostly appended defaults and stuff */

        if (!t->used)

            continue;

        /* grow the resulting string if needed */

        l = strlen(t->param) + 1;

        if( strlen(definition) + l + 5 > def_max )

        {

            char *def2;

            def_max = def_max * 2 + l + 5;

            def2 = (char *) pj_malloc(def_max);

            strcpy( def2, definition );

            pj_dalloc( definition );

            definition = def2;

        }

        /* append this parameter */

        strcat( definition, " +" );

        strcat( definition, t->param );

    }

    return definition;

}

#ifndef lint

static const char SCCSID[]="@(#)PJ_putp5.c	4.1	94/02/15	GIE	REL";

#endif

#define PROJ_PARMS__ \

	double	A, B;

#define PJ_LIB__

# include	<projects.h>

PROJ_HEAD(putp5, "Putnins P5") "\n\tPCyl., Sph.";

PROJ_HEAD(putp5p, "Putnins P5'") "\n\tPCyl., Sph.";

#define C	1.01346

#define D	1.2158542

FORWARD(s_forward); /* spheroid */

	xy.x = C * lp.lam * (P->A - P->B * sqrt(1. + D * lp.phi * lp.phi));

	xy.y = C * lp.phi;

	return (xy);

}

INVERSE(s_inverse); /* spheroid */

	lp.phi = xy.y / C;

	lp.lam = xy.x / (C * (P->A - P->B * sqrt(1. + D * lp.phi * lp.phi)));

	return (lp);

}

FREEUP; if (P) pj_dalloc(P); }

	static PJ *

setup(PJ *P) {

	P->es = 0.; P->inv = s_inverse; P->fwd = s_forward;

	return P;

}

ENTRY0(putp5) P->A = 2.; P->B = 1.; ENDENTRY(setup(P))

ENTRY0(putp5p) P->A = 1.5; P->B = 0.5; ENDENTRY(setup(P))

/* definition of standard geoids */

#ifndef lint

static const char SCCSID[]="@(#)pj_ellps.c	4.6	95/08/25	GIE	REL";

#endif

#define PJ_ELLPS__

#include "projects.h"

C_NAMESPACE struct PJ_ELLPS

pj_ellps[] = {

"MERIT",	"a=6378137.0", "rf=298.257", "MERIT 1983",

"SGS85",	"a=6378136.0", "rf=298.257",  "Soviet Geodetic System 85",

"GRS80",	"a=6378137.0", "rf=298.257222101", "GRS 1980(IUGG, 1980)",

"IAU76",	"a=6378140.0", "rf=298.257", "IAU 1976",

"airy",		"a=6377563.396", "b=6356256.910", "Airy 1830",

"APL4.9",	"a=6378137.0.",  "rf=298.25", "Appl. Physics. 1965",

"NWL9D",	"a=6378145.0.",  "rf=298.25", "Naval Weapons Lab., 1965",

"mod_airy",	"a=6377340.189", "b=6356034.446", "Modified Airy",

"andrae",	"a=6377104.43",  "rf=300.0", 	"Andrae 1876 (Den., Iclnd.)",

"aust_SA",	"a=6378160.0", "rf=298.25", "Australian Natl & S. Amer. 1969",

"GRS67",	"a=6378160.0", "rf=298.2471674270", "GRS 67(IUGG 1967)",

"bessel",	"a=6377397.155", "rf=299.1528128", "Bessel 1841",

"bess_nam",	"a=6377483.865", "rf=299.1528128", "Bessel 1841 (Namibia)",

"clrk66",	"a=6378206.4", "b=6356583.8", "Clarke 1866",

"clrk80",	"a=6378249.145", "rf=293.4663", "Clarke 1880 mod.",

"CPM",  	"a=6375738.7", "rf=334.29", "Comm. des Poids et Mesures 1799",

"delmbr",	"a=6376428.",  "rf=311.5", "Delambre 1810 (Belgium)",

"engelis",	"a=6378136.05", "rf=298.2566", "Engelis 1985",

"evrst30",  "a=6377276.345", "rf=300.8017",  "Everest 1830",

"evrst48",  "a=6377304.063", "rf=300.8017",  "Everest 1948",

"evrst56",  "a=6377301.243", "rf=300.8017",  "Everest 1956",

"evrst69",  "a=6377295.664", "rf=300.8017",  "Everest 1969",

"evrstSS",  "a=6377298.556", "rf=300.8017",  "Everest (Sabah & Sarawak)",

"fschr60",  "a=6378166.",   "rf=298.3", "Fischer (Mercury Datum) 1960",

"fschr60m", "a=6378155.",   "rf=298.3", "Modified Fischer 1960",

"fschr68",  "a=6378150.",   "rf=298.3", "Fischer 1968",

"helmert",  "a=6378200.",   "rf=298.3", "Helmert 1906",

"hough",	"a=6378270.0", "rf=297.", "Hough",

"intl",		"a=6378388.0", "rf=297.", "International 1909 (Hayford)",

"krass",	"a=6378245.0", "rf=298.3", "Krassovsky, 1942",

"kaula",	"a=6378163.",  "rf=298.24", "Kaula 1961",

"lerch",	"a=6378139.",  "rf=298.257", "Lerch 1979",

"mprts",	"a=6397300.",  "rf=191.", "Maupertius 1738",

"new_intl",	"a=6378157.5", "b=6356772.2", "New International 1967",

"plessis",	"a=6376523.",  "b=6355863.", "Plessis 1817 (France)",

"SEasia",	"a=6378155.0", "b=6356773.3205", "Southeast Asia",

"walbeck",	"a=6376896.0", "b=6355834.8467", "Walbeck",

"WGS60",    "a=6378165.0",  "rf=298.3", "WGS 60",

"WGS66",	"a=6378145.0", "rf=298.25", "WGS 66",

"WGS72",	"a=6378135.0", "rf=298.26", "WGS 72",

"WGS84",    "a=6378137.0",  "rf=298.257223563", "WGS 84",

"sphere",   "a=6370997.0",  "b=6370997.0", "Normal Sphere (r=6370997)",

0,0,0,0

};

struct PJ_ELLPS *pj_get_ellps_ref()

{

    return pj_ellps;

}

#ifndef lint

static const char SCCSID[]="@(#)PJ_stere.c	4.1	94/02/15	GIE	REL";

#endif

#define PROJ_PARMS__ \

	double phits; \

	double sinX1; \

	double cosX1; \

	double akm1; \

	int	mode;

#define PJ_LIB__

#include	<projects.h>

PROJ_HEAD(stere, "Stereographic") "\n\tAzi, Sph&Ell\n\tlat_ts=";

PROJ_HEAD(ups, "Universal Polar Stereographic") "\n\tAzi, Sph&Ell\n\tsouth";

#define sinph0	P->sinX1

#define cosph0	P->cosX1

#define EPS10	1.e-10

#define TOL	1.e-8

#define NITER	8

#define CONV	1.e-10

#define S_POLE	0

#define N_POLE	1

#define OBLIQ	2

#define EQUIT	3

	static double

ssfn_(double phit, double sinphi, double eccen) {

	sinphi *= eccen;

	return (tan (.5 * (HALFPI + phit)) *

	   pow((1. - sinphi) / (1. + sinphi), .5 * eccen));

}

FORWARD(e_forward); /* ellipsoid */

	double coslam, sinlam, sinX=0.0, cosX=0.0, X, A, sinphi;

	coslam = cos(lp.lam);

	sinlam = sin(lp.lam);

	sinphi = sin(lp.phi);

	if (P->mode == OBLIQ || P->mode == EQUIT) {

		sinX = sin(X = 2. * atan(ssfn_(lp.phi, sinphi, P->e)) - HALFPI);

		cosX = cos(X);

	}

	switch (P->mode) {

	case OBLIQ:

		A = P->akm1 / (P->cosX1 * (1. + P->sinX1 * sinX +

		   P->cosX1 * cosX * coslam));

		xy.y = A * (P->cosX1 * sinX - P->sinX1 * cosX * coslam);

		goto xmul;

	case EQUIT:

		A = 2. * P->akm1 / (1. + cosX * coslam);

		xy.y = A * sinX;

xmul:

		xy.x = A * cosX;

		break;

	case S_POLE:

		lp.phi = -lp.phi;

		coslam = - coslam;

		sinphi = -sinphi;

	case N_POLE:

		xy.x = P->akm1 * pj_tsfn(lp.phi, sinphi, P->e);

		xy.y = - xy.x * coslam;

		break;

	}

	xy.x = xy.x * sinlam;

	return (xy);

}

FORWARD(s_forward); /* spheroid */

	double  sinphi, cosphi, coslam, sinlam;

	sinphi = sin(lp.phi);

	cosphi = cos(lp.phi);

	coslam = cos(lp.lam);

	sinlam = sin(lp.lam);

	switch (P->mode) {

	case EQUIT:

		xy.y = 1. + cosphi * coslam;

		goto oblcon;

	case OBLIQ:

		xy.y = 1. + sinph0 * sinphi + cosph0 * cosphi * coslam;

oblcon:

		if (xy.y <= EPS10) F_ERROR;

		xy.x = (xy.y = P->akm1 / xy.y) * cosphi * sinlam;

		xy.y *= (P->mode == EQUIT) ? sinphi :

		   cosph0 * sinphi - sinph0 * cosphi * coslam;

		break;

	case N_POLE:

		coslam = - coslam;

		lp.phi = - lp.phi;

	case S_POLE:

		if (fabs(lp.phi - HALFPI) < TOL) F_ERROR;

		xy.x = sinlam * ( xy.y = P->akm1 * tan(FORTPI + .5 * lp.phi) );

		xy.y *= coslam;

		break;

	}

	return (xy);

}

INVERSE(e_inverse); /* ellipsoid */

	double cosphi, sinphi, tp=0.0, phi_l=0.0, rho, halfe=0.0, halfpi=0.0;

	int i;

	rho = hypot(xy.x, xy.y);

	switch (P->mode) {

	case OBLIQ:

	case EQUIT:

		cosphi = cos( tp = 2. * atan2(rho * P->cosX1 , P->akm1) );

		sinphi = sin(tp);

                if( rho == 0.0 )

		    phi_l = asin(cosphi * P->sinX1);

                else

		    phi_l = asin(cosphi * P->sinX1 + (xy.y * sinphi * P->cosX1 / rho));

		tp = tan(.5 * (HALFPI + phi_l));

		xy.x *= sinphi;

		xy.y = rho * P->cosX1 * cosphi - xy.y * P->sinX1* sinphi;

		halfpi = HALFPI;

		halfe = .5 * P->e;

		break;

	case N_POLE:

		xy.y = -xy.y;

	case S_POLE:

		phi_l = HALFPI - 2. * atan(tp = - rho / P->akm1);

		halfpi = -HALFPI;

		halfe = -.5 * P->e;

		break;

	}

	for (i = NITER; i--; phi_l = lp.phi) {

		sinphi = P->e * sin(phi_l);

		lp.phi = 2. * atan(tp * pow((1.+sinphi)/(1.-sinphi),

		   halfe)) - halfpi;

		if (fabs(phi_l - lp.phi) < CONV) {

			if (P->mode == S_POLE)

				lp.phi = -lp.phi;

			lp.lam = (xy.x == 0. && xy.y == 0.) ? 0. : atan2(xy.x, xy.y);

			return (lp);

		}

	}

	I_ERROR;

}

INVERSE(s_inverse); /* spheroid */

	double  c, rh, sinc, cosc;

	sinc = sin(c = 2. * atan((rh = hypot(xy.x, xy.y)) / P->akm1));

	cosc = cos(c);

	lp.lam = 0.;

	switch (P->mode) {

	case EQUIT:

		if (fabs(rh) <= EPS10)

			lp.phi = 0.;

		else

			lp.phi = asin(xy.y * sinc / rh);

		if (cosc != 0. || xy.x != 0.)

			lp.lam = atan2(xy.x * sinc, cosc * rh);

		break;

	case OBLIQ:

		if (fabs(rh) <= EPS10)

			lp.phi = P->phi0;

		else

			lp.phi = asin(cosc * sinph0 + xy.y * sinc * cosph0 / rh);

		if ((c = cosc - sinph0 * sin(lp.phi)) != 0. || xy.x != 0.)

			lp.lam = atan2(xy.x * sinc * cosph0, c * rh);

		break;

	case N_POLE:

		xy.y = -xy.y;

	case S_POLE:

		if (fabs(rh) <= EPS10)

			lp.phi = P->phi0;

		else

			lp.phi = asin(P->mode == S_POLE ? - cosc : cosc);

		lp.lam = (xy.x == 0. && xy.y == 0.) ? 0. : atan2(xy.x, xy.y);

		break;

	}

	return (lp);

}

FREEUP; if (P) pj_dalloc(P); }

	static PJ *

setup(PJ *P) { /* general initialization */

	double t;

	if (fabs((t = fabs(P->phi0)) - HALFPI) < EPS10)

		P->mode = P->phi0 < 0. ? S_POLE : N_POLE;

	else

		P->mode = t > EPS10 ? OBLIQ : EQUIT;

	P->phits = fabs(P->phits);

	if (P->es) {

		double X;

		switch (P->mode) {

		case N_POLE:

		case S_POLE:

			if (fabs(P->phits - HALFPI) < EPS10)

				P->akm1 = 2. * P->k0 /

				   sqrt(pow(1+P->e,1+P->e)*pow(1-P->e,1-P->e));

			else {

				P->akm1 = cos(P->phits) /

				   pj_tsfn(P->phits, t = sin(P->phits), P->e);

				t *= P->e;

				P->akm1 /= sqrt(1. - t * t);

			}

			break;

		case EQUIT:

			P->akm1 = 2. * P->k0;

			break;

		case OBLIQ:

			t = sin(P->phi0);

			X = 2. * atan(ssfn_(P->phi0, t, P->e)) - HALFPI;

			t *= P->e;

			P->akm1 = 2. * P->k0 * cos(P->phi0) / sqrt(1. - t * t);

			P->sinX1 = sin(X);

			P->cosX1 = cos(X);

			break;

		}

		P->inv = e_inverse;

		P->fwd = e_forward;

	} else {

		switch (P->mode) {

		case OBLIQ:

			sinph0 = sin(P->phi0);

			cosph0 = cos(P->phi0);

		case EQUIT:

			P->akm1 = 2. * P->k0;

			break;

		case S_POLE:

		case N_POLE:

			P->akm1 = fabs(P->phits - HALFPI) >= EPS10 ?

			   cos(P->phits) / tan(FORTPI - .5 * P->phits) :

			   2. * P->k0 ;

			break;

		}

		P->inv = s_inverse;

		P->fwd = s_forward;

	}

	return P;

}

ENTRY0(stere)

	P->phits = pj_param(P->params, "tlat_ts").i ?

		P->phits = pj_param(P->params, "rlat_ts").f : HALFPI;

ENDENTRY(setup(P))

ENTRY0(ups)

	/* International Ellipsoid */

	P->phi0 = pj_param(P->params, "bsouth").i ? - HALFPI: HALFPI;

	if (!P->es) E_ERROR(-34);

	P->k0 = .994;

	P->x0 = 2000000.;

	P->y0 = 2000000.;

	P->phits = HALFPI;

	P->lam0 = 0.;

ENDENTRY(setup(P))

/*

** libproj -- library of cartographic projections

**

** Copyright (c) 2003   Gerald I. Evenden

*/

static const char

LIBPROJ_ID[] = "$Id: PJ_sterea.c,v 1.1 2004/10/20 17:04:00 fwarmerdam Exp $";

/*

** Permission is hereby granted, free of charge, to any person obtaining

** a copy of this software and associated documentation files (the

** "Software"), to deal in the Software without restriction, including

** without limitation the rights to use, copy, modify, merge, publish,

** distribute, sublicense, and/or sell copies of the Software, and to

** permit persons to whom the Software is furnished to do so, subject to

** the following conditions:

**

** The above copyright notice and this permission notice shall be

** included in all copies or substantial portions of the Software.

**

** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,

** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF

** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.

** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY

** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,

** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE

** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

*/

#define PROJ_PARMS__ \

	double phic0; \

	double cosc0, sinc0; \

	double R2; \

	void *en;

#define PJ_LIB__

#include	<projects.h>

PROJ_HEAD(sterea, "Oblique Stereographic Alternative")

	"\n\tAzimuthal, Sph&Ell";

# define DEL_TOL	1.e-14

# define MAX_ITER	10

FORWARD(e_forward); /* ellipsoid */

	double cosc, sinc, cosl, k;

	lp = pj_gauss(lp, P->en);

	sinc = sin(lp.phi);

	cosc = cos(lp.phi);

	cosl = cos(lp.lam);

	k = P->k0 * P->R2 / (1. + P->sinc0 * sinc + P->cosc0 * cosc * cosl);

	xy.x = k * cosc * sin(lp.lam);

	xy.y = k * (P->cosc0 * sinc - P->sinc0 * cosc * cosl);

	return (xy);

}

INVERSE(e_inverse); /* ellipsoid */

	double rho, c, sinc, cosc;

	xy.x /= P->k0;

	xy.y /= P->k0;

	if((rho = hypot(xy.x, xy.y))) {

		c = 2. * atan2(rho, P->R2);

		sinc = sin(c);

		cosc = cos(c);

		lp.phi = asin(cosc * P->sinc0 + xy.y * sinc * P->cosc0 / rho);

		lp.lam = atan2(xy.x * sinc, rho * P->cosc0 * cosc -

			xy.y * P->sinc0 * sinc);

	} else {

		lp.phi = P->phic0;

		lp.lam = 0.;

	}

	return(pj_inv_gauss(lp, P->en));

}

FREEUP; if (P) { if (P->en) free(P->en); free(P); } }

ENTRY0(sterea)

	double R;

	if (!(P->en = pj_gauss_ini(P->e, P->phi0, &(P->phic0), &R))) E_ERROR_0;

	P->sinc0 = sin(P->phic0);

	P->cosc0 = cos(P->phic0);

	P->R2 = 2. * R;

	P->inv = e_inverse;

	P->fwd = e_forward;

ENDENTRY(P)

/*

** $Log: PJ_sterea.c,v $

** Revision 1.1  2004/10/20 17:04:00  fwarmerdam

** New

**

** Revision 2.3  2004/04/07 17:18:32  gie

** corrected comment stamp

**

** Revision 2.2  2003/08/05 00:15:09  gie

** corrected 0 rho on inverse.

**

** Revision 2.1  2003/03/28 01:46:02  gie

** Initial

**

*/

#ifndef lint

static const char SCCSID[]="@(#)PJ_eck1.c	4.1 94/02/15     GIE     REL";

#endif

#define PJ_LIB__

#include	<projects.h>

PROJ_HEAD(eck1, "Eckert I") "\n\tPCyl., Sph.";

#define FC	.92131773192356127802

#define RP	.31830988618379067154

FORWARD(s_forward); /* spheroid */

	xy.x = FC * lp.lam * (1. - RP * fabs(lp.phi));

	xy.y = FC * lp.phi;

	return (xy);

}

INVERSE(s_inverse); /* spheroid */

	lp.phi = xy.y / FC;

	lp.lam = xy.x / (FC * (1. - RP * fabs(lp.phi)));

	return (lp);

}

FREEUP; if (P) pj_dalloc(P); }

ENTRY0(eck1)

	P->es = 0.; P->inv = s_inverse; P->fwd = s_forward;

ENDENTRY(P)

#ifndef lint

static const char SCCSID[]="@(#)PJ_labrd.c	4.1	94/02/15	GIE	REL";

#endif

#define PROJ_PARMS__ \

	double	Az, kRg, p0s, A, C, Ca, Cb, Cc, Cd; \

	int		rot;

#define PJ_LIB__

#include	<projects.h>

PROJ_HEAD(labrd, "Laborde") "\n\tCyl, Sph\n\tSpecial for Madagascar";

#define EPS	1.e-10

FORWARD(e_forward);

	double V1, V2, ps, sinps, cosps, sinps2, cosps2, I1, I2, I3, I4, I5, I6,

		x2, y2, t;

	V1 = P->A * log( tan(FORTPI + .5 * lp.phi) );

	t = P->e * sin(lp.phi);

	V2 = .5 * P->e * P->A * log ((1. + t)/(1. - t));

	ps = 2. * (atan(exp(V1 - V2 + P->C)) - FORTPI);

	I1 = ps - P->p0s;

	cosps = cos(ps);	cosps2 = cosps * cosps;

	sinps = sin(ps);	sinps2 = sinps * sinps;

	I4 = P->A * cosps;

	I2 = .5 * P->A * I4 * sinps;

	I3 = I2 * P->A * P->A * (5. * cosps2 - sinps2) / 12.;

	I6 = I4 * P->A * P->A;

	I5 = I6 * (cosps2 - sinps2) / 6.;

	I6 *= P->A * P->A *

		(5. * cosps2 * cosps2 + sinps2 * (sinps2 - 18. * cosps2)) / 120.;

	t = lp.lam * lp.lam;

	xy.x = P->kRg * lp.lam * (I4 + t * (I5 + t * I6));

	xy.y = P->kRg * (I1 + t * (I2 + t * I3));

	x2 = xy.x * xy.x;

	y2 = xy.y * xy.y;

	V1 = 3. * xy.x * y2 - xy.x * x2;

	V2 = xy.y * y2 - 3. * x2 * xy.y;

	xy.x += P->Ca * V1 + P->Cb * V2;

	xy.y += P->Ca * V2 - P->Cb * V1;

	return (xy);

}

INVERSE(e_inverse); /* ellipsoid & spheroid */

	double x2, y2, V1, V2, V3, V4, t, t2, ps, pe, tpe, s,

		I7, I8, I9, I10, I11, d, Re;

	int i;

	x2 = xy.x * xy.x;

	y2 = xy.y * xy.y;

	V1 = 3. * xy.x * y2 - xy.x * x2;

	V2 = xy.y * y2 - 3. * x2 * xy.y;

	V3 = xy.x * (5. * y2 * y2 + x2 * (-10. * y2 + x2 ));

	V4 = xy.y * (5. * x2 * x2 + y2 * (-10. * x2 + y2 ));

	xy.x += - P->Ca * V1 - P->Cb * V2 + P->Cc * V3 + P->Cd * V4;

	xy.y +=   P->Cb * V1 - P->Ca * V2 - P->Cd * V3 + P->Cc * V4;

	ps = P->p0s + xy.y / P->kRg;

	pe = ps + P->phi0 - P->p0s;

	for ( i = 20; i; --i) {

		V1 = P->A * log(tan(FORTPI + .5 * pe));

		tpe = P->e * sin(pe);

		V2 = .5 * P->e * P->A * log((1. + tpe)/(1. - tpe));

		t = ps - 2. * (atan(exp(V1 - V2 + P->C)) - FORTPI);

		pe += t;

		if (fabs(t) < EPS)

			break;

	}

/*

	if (!i) {

	} else {

	}

*/

	t = P->e * sin(pe);

	t = 1. - t * t;

	Re = P->one_es / ( t * sqrt(t) );

	t = tan(ps);

	t2 = t * t;

	s = P->kRg * P->kRg;

	d = Re * P->k0 * P->kRg;

	I7 = t / (2. * d);

	I8 = t * (5. + 3. * t2) / (24. * d * s);

	d = cos(ps) * P->kRg * P->A;

	I9 = 1. / d;

	d *= s;

	I10 = (1. + 2. * t2) / (6. * d);

	I11 = (5. + t2 * (28. + 24. * t2)) / (120. * d * s);

	x2 = xy.x * xy.x;

	lp.phi = pe + x2 * (-I7 + I8 * x2);

	lp.lam = xy.x * (I9 + x2 * (-I10 + x2 * I11));

	return (lp);

}

FREEUP; if (P) pj_dalloc(P); }

ENTRY0(labrd)

	double Az, sinp, R, N, t;

	P->rot	= pj_param(P->params, "bno_rot").i == 0;

	Az = pj_param(P->params, "razi").f;

	sinp = sin(P->phi0);

	t = 1. - P->es * sinp * sinp;

	N = 1. / sqrt(t);

	R = P->one_es * N / t;

	P->kRg = P->k0 * sqrt( N * R );

	P->p0s = atan( sqrt(R / N) * tan(P->phi0) );

	P->A = sinp / sin(P->p0s);

	t = P->e * sinp;

	P->C = .5 * P->e * P->A * log((1. + t)/(1. - t)) +

		- P->A * log( tan(FORTPI + .5 * P->phi0))

		+ log( tan(FORTPI + .5 * P->p0s));

	t = Az + Az;

	P->Ca = (1. - cos(t)) * ( P->Cb = 1. / (12. * P->kRg * P->kRg) );

	P->Cb *= sin(t);

	P->Cc = 3. * (P->Ca * P->Ca - P->Cb * P->Cb);

	P->Cd = 6. * P->Ca * P->Cb;

	P->inv = e_inverse;

	P->fwd = e_forward;

ENDENTRY(P)

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

 * $Id: pj_geocent.c,v 1.1 2002/12/14 20:13:21 warmerda Exp $

 *

 * Project:  PROJ.4

 * Purpose:  Stub projection for geocentric.  The transformation isn't

 *           really done here since this code is 2D.  The real transformation

 *           is handled by pj_transform.c.

 * Author:   Frank Warmerdam, warmerdam@pobox.com

 *

 ******************************************************************************

 * Copyright (c) 2002, Frank Warmerdam

 *

 * Permission is hereby granted, free of charge, to any person obtaining a

 * copy of this software and associated documentation files (the "Software"),

 * to deal in the Software without restriction, including without limitation

 * the rights to use, copy, modify, merge, publish, distribute, sublicense,

 * and/or sell copies of the Software, and to permit persons to whom the

 * Software is furnished to do so, subject to the following conditions:

 *

 * The above copyright notice and this permission notice shall be included

 * in all copies or substantial portions of the Software.

 *

 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER

 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING

 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER

 * DEALINGS IN THE SOFTWARE.

 ******************************************************************************

 *

 * $Log: pj_geocent.c,v $

 * Revision 1.1  2002/12/14 20:13:21  warmerda

 * New

 *

 */

#define PJ_LIB__

#include <projects.h>

PJ_CVSID("$Id: pj_geocent.c,v 1.1 2002/12/14 20:13:21 warmerda Exp $");

PROJ_HEAD(geocent, "Geocentric")  "\n\t";

FORWARD(forward);

        xy.x = lp.lam;

        xy.y = lp.phi;

        return xy;

}

INVERSE(inverse);

        lp.phi = xy.y;

        lp.lam = xy.x;

        return lp;

}

FREEUP; if (P) pj_dalloc(P); }

ENTRY0(geocent)

    P->is_geocent = 1; 

    P->x0 = 0.0;

    P->y0 = 0.0;

    P->inv = inverse; P->fwd = forward;

ENDENTRY(P)

/* For full ANSI compliance of global variable */

#ifndef lint

static const char SCCSID[]="@(#)pj_errno.c	4.3	95/06/03	GIE	REL";

#endif

#include <projects.h>

C_NAMESPACE int pj_errno = 0;

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

/*                          pj_get_errno_ref()                          */

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

int *pj_get_errno_ref()

{

    return &pj_errno;

}

/* end */

#ifndef lint

static const char SCCSID[]="@(#)PJ_eck5.c	4.1 94/02/15     GIE     REL";

#endif

#define PJ_LIB__

# include	<projects.h>

PROJ_HEAD(eck5, "Eckert V") "\n\tPCyl, Sph.";

#define XF	0.44101277172455148219

#define RXF	2.26750802723822639137

#define YF	0.88202554344910296438

#define RYF	1.13375401361911319568

FORWARD(s_forward); /* spheroid */

	xy.x = XF * (1. + cos(lp.phi)) * lp.lam;

	xy.y = YF * lp.phi;

	return (xy);

}

INVERSE(s_inverse); /* spheroid */

	lp.lam = RXF * xy.x / (1. + cos( lp.phi = RYF * xy.y));

	return (lp);

}

FREEUP; if (P) pj_dalloc(P); }

ENTRY0(eck5); P->es = 0.; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)

/* set ellipsoid parameters a and es */

#ifndef lint

static const char SCCSID[]="@(#)pj_ell_set.c	4.5	93/06/12	GIE	REL";

#endif

#include <projects.h>

#include <string.h>

#define SIXTH .1666666666666666667 /* 1/6 */

#define RA4 .04722222222222222222 /* 17/360 */

#define RA6 .02215608465608465608 /* 67/3024 */

#define RV4 .06944444444444444444 /* 5/72 */

#define RV6 .04243827160493827160 /* 55/1296 */

	int /* initialize geographic shape parameters */

pj_ell_set(paralist *pl, double *a, double *es) {

	int i;

	double b=0.0, e;

	char *name;

	paralist *start = 0, *curr;

		/* check for varying forms of ellipsoid input */

	*a = *es = 0.;

	/* R takes precedence */

	if (pj_param(pl, "tR").i)

		*a = pj_param(pl, "dR").f;

	else { /* probable elliptical figure */

		/* check if ellps present and temporarily append its values to pl */

		if (name = pj_param(pl, "sellps").s) {

			char *s;

			for (start = pl; start && start->next ; start = start->next) ;

			curr = start;

			for (i = 0; (s = pj_ellps[i].id) && strcmp(name, s) ; ++i) ;

			if (!s) { pj_errno = -9; return 1; }

			curr = curr->next = pj_mkparam(pj_ellps[i].major);

			curr = curr->next = pj_mkparam(pj_ellps[i].ell);

		}

		*a = pj_param(pl, "da").f;

		if (pj_param(pl, "tes").i) /* eccentricity squared */

			*es = pj_param(pl, "des").f;

		else if (pj_param(pl, "te").i) { /* eccentricity */

			e = pj_param(pl, "de").f;

			*es = e * e;

		} else if (pj_param(pl, "trf").i) { /* recip flattening */

			*es = pj_param(pl, "drf").f;

			if (!*es) {

				pj_errno = -10;

				goto bomb;

			}

			*es = 1./ *es;

			*es = *es * (2. - *es);

		} else if (pj_param(pl, "tf").i) { /* flattening */

			*es = pj_param(pl, "df").f;

			*es = *es * (2. - *es);

		} else if (pj_param(pl, "tb").i) { /* minor axis */

			b = pj_param(pl, "db").f;

			*es = 1. - (b * b) / (*a * *a);

		}     /* else *es == 0. and sphere of radius *a */

		if (!b)

			b = *a * sqrt(1. - *es);

		/* following options turn ellipsoid into equivalent sphere */

		if (pj_param(pl, "bR_A").i) { /* sphere--area of ellipsoid */

			*a *= 1. - *es * (SIXTH + *es * (RA4 + *es * RA6));

			*es = 0.;

		} else if (pj_param(pl, "bR_V").i) { /* sphere--vol. of ellipsoid */

			*a *= 1. - *es * (SIXTH + *es * (RV4 + *es * RV6));

			*es = 0.;

		} else if (pj_param(pl, "bR_a").i) { /* sphere--arithmetic mean */

			*a = .5 * (*a + b);

			*es = 0.;

		} else if (pj_param(pl, "bR_g").i) { /* sphere--geometric mean */

			*a = sqrt(*a * b);

			*es = 0.;

		} else if (pj_param(pl, "bR_h").i) { /* sphere--harmonic mean */

			*a = 2. * *a * b / (*a + b);

			*es = 0.;

		} else if ((i = pj_param(pl, "tR_lat_a").i) || /* sphere--arith. */

			pj_param(pl, "tR_lat_g").i) { /* or geom. mean at latitude */

			double tmp;

			tmp = sin(pj_param(pl, i ? "rR_lat_a" : "rR_lat_g").f);

			if (fabs(tmp) > HALFPI) {

				pj_errno = -11;

				goto bomb;

			}

			tmp = 1. - *es * tmp * tmp;

			*a *= i ? .5 * (1. - *es + tmp) / ( tmp * sqrt(tmp)) :

				sqrt(1. - *es) / tmp;

			*es = 0.;

		}

bomb:

		if (start) { /* clean up temporary extension of list */

			pj_dalloc(start->next->next);

			pj_dalloc(start->next);

			start->next = 0;

		}

		if (pj_errno)

			return 1;

	}

	/* some remaining checks */

	if (*es < 0.)

		{ pj_errno = -12; return 1; }

	if (*a <= 0.)

		{ pj_errno = -13; return 1; }

	return 0;

}

#ifndef lint

static const char SCCSID[]="@(#)PJ_wink2.c	4.1 94/02/15     GIE     REL";

#endif

#define PROJ_PARMS__ \

	double	cosphi1;

#define PJ_LIB__

# include	<projects.h>

PROJ_HEAD(wink2, "Winkel II") "\n\tPCyl., Sph., no inv.\n\tlat_1=";

#define MAX_ITER    10

#define LOOP_TOL    1e-7

#define TWO_D_PI 0.636619772367581343

FORWARD(s_forward); /* spheroid */

	double k, V;

	int i;

	xy.y = lp.phi * TWO_D_PI;

	k = PI * sin(lp.phi);

	lp.phi *= 1.8;

	for (i = MAX_ITER; i ; --i) {

		lp.phi -= V = (lp.phi + sin(lp.phi) - k) /

			(1. + cos(lp.phi));

		if (fabs(V) < LOOP_TOL)

			break;

	}

	if (!i)

		lp.phi = (lp.phi < 0.) ? -HALFPI : HALFPI;

	else

		lp.phi *= 0.5;

	xy.x = 0.5 * lp.lam * (cos(lp.phi) + P->cosphi1);

	xy.y = FORTPI * (sin(lp.phi) + xy.y);

	return (xy);

}

FREEUP; if (P) pj_dalloc(P); }

ENTRY0(wink2)

	P->cosphi1 = cos(pj_param(P->params, "rlat_1").f);

	P->es = 0.; P->inv = 0; P->fwd = s_forward;

ENDENTRY(P)

/* make storage for one and two dimensional matricies */

#ifndef lint

static const char SCCSID[]="@(#)vector1.c	4.4	94/03/22	GIE	REL";

#endif

#include <stdlib.h>

#include <projects.h>

	void * /* one dimension array */

vector1(int nvals, int size) { return((void *)pj_malloc(size * nvals)); }

	void /* free 2D array */

freev2(void **v, int nrows) {

	if (v) {

		for (v += nrows; nrows > 0; --nrows)

			pj_dalloc(*--v);

		pj_dalloc(v);

	}

}

	void ** /* two dimension array */

vector2(int nrows, int ncols, int size) {

	void **s;

	if (s = (void **)pj_malloc(sizeof(void *) * nrows)) {

		int rsize, i;

		rsize = size * ncols;

		for (i = 0; i < nrows; ++i)

			if (!(s[i] = pj_malloc(rsize))) {

				freev2(s, i);

				return (void **)0;

			}

	}

	return s;

}

#ifndef lint

static const char SCCSID[]="@(#)PJ_nell.c	4.1	94/02/15	GIE	REL";

#endif

#define PJ_LIB__

#include	<projects.h>

PROJ_HEAD(nell, "Nell") "\n\tPCyl., Sph.";

#define MAX_ITER	10

#define LOOP_TOL	1e-7

FORWARD(s_forward); /* spheroid */

	double k, V;

	int i;

	k = 2. * sin(lp.phi);

	V = lp.phi * lp.phi;

	lp.phi *= 1.00371 + V * (-0.0935382 + V * -0.011412);

	for (i = MAX_ITER; i ; --i) {

		lp.phi -= V = (lp.phi + sin(lp.phi) - k) /

			(1. + cos(lp.phi));

		if (fabs(V) < LOOP_TOL)

			break;

	}

	xy.x = 0.5 * lp.lam * (1. + cos(lp.phi));

	xy.y = lp.phi;

	return (xy);

}

INVERSE(s_inverse); /* spheroid */

	double th, s;

	lp.lam = 2. * xy.x / (1. + cos(xy.y));

	lp.phi = aasin(0.5 * (xy.y + sin(xy.y)));

	return (lp);

}

FREEUP; if (P) pj_dalloc(P); }

ENTRY0(nell) P->es = 0; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)

#ifndef lint

static const char SCCSID[]="@(#)PJ_ortho.c	4.1	94/02/15	GIE	REL";

#endif

#define PROJ_PARMS__ \

	double	sinph0; \

	double	cosph0; \

	int		mode;

#define PJ_LIB__

#include	<projects.h>

PROJ_HEAD(ortho, "Orthographic") "\n\tAzi, Sph.";

#define EPS10 1.e-10

#define N_POLE	0

#define S_POLE 1

#define EQUIT	2

#define OBLIQ	3

FORWARD(s_forward); /* spheroid */

	double  coslam, cosphi, sinphi;

	cosphi = cos(lp.phi);

	coslam = cos(lp.lam);

	switch (P->mode) {

	case EQUIT:

		if (cosphi * coslam < - EPS10) F_ERROR;

		xy.y = sin(lp.phi);

		break;

	case OBLIQ:

		if (P->sinph0 * (sinphi = sin(lp.phi)) +

		   P->cosph0 * cosphi * coslam < - EPS10) F_ERROR;

		xy.y = P->cosph0 * sinphi - P->sinph0 * cosphi * coslam;

		break;

	case N_POLE:

		coslam = - coslam;

	case S_POLE:

		if (fabs(lp.phi - P->phi0) - EPS10 > HALFPI) F_ERROR;

		xy.y = cosphi * coslam;

		break;

	}

	xy.x = cosphi * sin(lp.lam);

	return (xy);

}

INVERSE(s_inverse); /* spheroid */

    double  rh, cosc, sinc;

    if ((sinc = (rh = hypot(xy.x, xy.y))) > 1.) {

        if ((sinc - 1.) > EPS10) I_ERROR;

        sinc = 1.;

    }

    cosc = sqrt(1. - sinc * sinc); /* in this range OK */

    if (fabs(rh) <= EPS10) {

        lp.phi = P->phi0;

        lp.lam = 0.0;

    } else {

        switch (P->mode) {

        case N_POLE:

            xy.y = -xy.y;

            lp.phi = acos(sinc);

            break;

        case S_POLE:

            lp.phi = - acos(sinc);

            break;

        case EQUIT:

            lp.phi = xy.y * sinc / rh;

            xy.x *= sinc;

            xy.y = cosc * rh;

            goto sinchk;

        case OBLIQ:

            lp.phi = cosc * P->sinph0 + xy.y * sinc * P->cosph0 /rh;

            xy.y = (cosc - P->sinph0 * lp.phi) * rh;

            xy.x *= sinc * P->cosph0;

        sinchk:

            if (fabs(lp.phi) >= 1.)

                lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;

            else

                lp.phi = asin(lp.phi);

            break;

        }

        lp.lam = (xy.y == 0. && (P->mode == OBLIQ || P->mode == EQUIT))

             ? (xy.x == 0. ? 0. : xy.x < 0. ? -HALFPI : HALFPI)

                           : atan2(xy.x, xy.y);

    }

    return (lp);

}

FREEUP; if (P) pj_dalloc(P); }

ENTRY0(ortho)

	if (fabs(fabs(P->phi0) - HALFPI) <= EPS10)

		P->mode = P->phi0 < 0. ? S_POLE : N_POLE;

	else if (fabs(P->phi0) > EPS10) {

		P->mode = OBLIQ;

		P->sinph0 = sin(P->phi0);

		P->cosph0 = cos(P->phi0);

	} else

		P->mode = EQUIT;

	P->inv = s_inverse;

	P->fwd = s_forward;

	P->es = 0.;

ENDENTRY(P)

/* generate double bivariate Chebychev polynomial */

#ifndef lint

static const char SCCSID[]="@(#)bchgen.c	4.5	94/03/22	GIE	REL";

#endif

#include <projects.h>

	int

bchgen(projUV a, projUV b, int nu, int nv, projUV **f, projUV(*func)(projUV)) {

	int i, j, k;

	projUV arg, *t, bma, bpa, *c;

	double d, fac;

	bma.u = 0.5 * (b.u - a.u); bma.v = 0.5 * (b.v - a.v);

	bpa.u = 0.5 * (b.u + a.u); bpa.v = 0.5 * (b.v + a.v);

	for ( i = 0; i < nu; ++i) {

		arg.u = cos(PI * (i + 0.5) / nu) * bma.u + bpa.u;

		for ( j = 0; j < nv; ++j) {

			arg.v = cos(PI * (j + 0.5) / nv) * bma.v + bpa.v;

			f[i][j] = (*func)(arg);

			if ((f[i][j]).u == HUGE_VAL)

				return(1);

		}

	}

	if (!(c = (projUV *) vector1(nu, sizeof(projUV)))) return 1;

	fac = 2. / nu;

	for ( j = 0; j < nv ; ++j) {

		for ( i = 0; i < nu; ++i) {

			arg.u = arg.v = 0.;

			for (k = 0; k < nu; ++k) {

				d = cos(PI * i * (k + .5) / nu);

				arg.u += f[k][j].u * d;

				arg.v += f[k][j].v * d;

			}

			arg.u *= fac;

			arg.v *= fac;

			c[i] = arg;

		}

		for (i = 0; i < nu; ++i)

			f[i][j] = c[i];

	}

	pj_dalloc(c);

	if (!(c = (projUV*) vector1(nv, sizeof(projUV)))) return 1;

	fac = 2. / nv;

	for ( i = 0; i < nu; ++i) {

		t = f[i];

		for (j = 0; j < nv; ++j) {

			arg.u = arg.v = 0.;

			for (k = 0; k < nv; ++k) {

				d = cos(PI * j * (k + .5) / nv);

				arg.u += t[k].u * d;

				arg.v += t[k].v * d;

			}

			arg.u *= fac;

			arg.v *= fac;

			c[j] = arg;

		}

		f[i] = c;

		c = t;

	}

	pj_dalloc(c);

	return(0);

}

/* evaluate complex polynomial */

#ifndef lint

static const char SCCSID[]="@(#)pj_zpoly1.c	4.3	93/06/12	GIE	REL";

#endif

#include <projects.h>

/* note: coefficients are always from C_1 to C_n

**	i.e. C_0 == (0., 0)

**	n should always be >= 1 though no checks are made

*/

	COMPLEX

pj_zpoly1(COMPLEX z, COMPLEX *C, int n) {

	COMPLEX a;

	double t;

	a = *(C += n);

	while (n-- > 0) {

		a.r = (--C)->r + z.r * (t = a.r) - z.i * a.i;

		a.i = C->i + z.r * a.i + z.i * t;

	}

	a.r = z.r * (t = a.r) - z.i * a.i;

	a.i = z.r * a.i + z.i * t;

	return a;

}

/* evaluate complex polynomial and derivative */

	COMPLEX

pj_zpolyd1(COMPLEX z, COMPLEX *C, int n, COMPLEX *der) {

	COMPLEX a, b;

	double t;

	int first = 1;

	a = *(C += n);

	while (n-- > 0) {

		if (first) {

			first = 0;

			b = a;

		} else {

			b.r = a.r + z.r * (t = b.r) - z.i * b.i;

			b.i = a.i + z.r * b.i + z.i * t;

		}

		a.r = (--C)->r + z.r * (t = a.r) - z.i * a.i;

		a.i = C->i + z.r * a.i + z.i * t;

	}

	b.r = a.r + z.r * (t = b.r) - z.i * b.i;

	b.i = a.i + z.r * b.i + z.i * t;

	a.r = z.r * (t = a.r) - z.i * a.i;

	a.i = z.r * a.i + z.i * t;

	*der = b;

	return a;

}

#ifndef lint

static const char SCCSID[]="@(#)PJ_vandg.c	4.1	94/02/15	GIE	REL";

#endif

#define PJ_LIB__

# include	<projects.h>

PROJ_HEAD(vandg, "van der Grinten (I)") "\n\tMisc Sph";

# define TOL		1.e-10

# define THIRD		.33333333333333333333

# define TWO_THRD	.66666666666666666666

# define C2_27		.07407407407407407407

# define PI4_3		4.18879020478639098458

# define PISQ		9.86960440108935861869

# define TPISQ		19.73920880217871723738

# define HPISQ		4.93480220054467930934

FORWARD(s_forward); /* spheroid */

	double  al, al2, g, g2, p2;

	p2 = fabs(lp.phi / HALFPI);

	if ((p2 - TOL) > 1.) F_ERROR;

	if (p2 > 1.)

		p2 = 1.;

	if (fabs(lp.phi) <= TOL) {

		xy.x = lp.lam;

		xy.y = 0.;

	} else if (fabs(lp.lam) <= TOL || fabs(p2 - 1.) < TOL) {

		xy.x = 0.;

		xy.y = PI * tan(.5 * asin(p2));

		if (lp.phi < 0.) xy.y = -xy.y;

	} else {

		al = .5 * fabs(PI / lp.lam - lp.lam / PI);

		al2 = al * al;

		g = sqrt(1. - p2 * p2);

		g = g / (p2 + g - 1.);

		g2 = g * g;

		p2 = g * (2. / p2 - 1.);

		p2 = p2 * p2;

		xy.x = g - p2; g = p2 + al2;

		xy.x = PI * (al * xy.x + sqrt(al2 * xy.x * xy.x - g * (g2 - p2))) / g;

		if (lp.lam < 0.) xy.x = -xy.x;

		xy.y = fabs(xy.x / PI);

		xy.y = 1. - xy.y * (xy.y + 2. * al);

		if (xy.y < -TOL) F_ERROR;

		if (xy.y < 0.)	xy.y = 0.;

		else		xy.y = sqrt(xy.y) * (lp.phi < 0. ? -PI : PI);

	}

	return (xy);

}

INVERSE(s_inverse); /* spheroid */

	double t, c0, c1, c2, c3, al, r2, r, m, d, ay, x2, y2;

	x2 = xy.x * xy.x;

	if ((ay = fabs(xy.y)) < TOL) {

		lp.phi = 0.;

		t = x2 * x2 + TPISQ * (x2 + HPISQ);

		lp.lam = fabs(xy.x) <= TOL ? 0. :

		   .5 * (x2 - PISQ + sqrt(t)) / xy.x;

		return (lp);

	}

	y2 = xy.y * xy.y;

	r = x2 + y2;	r2 = r * r;

	c1 = - PI * ay * (r + PISQ);

	c3 = r2 + TWOPI * (ay * r + PI * (y2 + PI * (ay + HALFPI)));

	c2 = c1 + PISQ * (r - 3. *  y2);

	c0 = PI * ay;

	c2 /= c3;

	al = c1 / c3 - THIRD * c2 * c2;