New

git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@1229 4e78687f-474d-0410-85f9-8d5e500ac6b2

3 files changed, 355 insertions, 0 deletions

diff --git a/src/PJ_geos.c b/src/PJ_geos.c
new file mode 100644
index 00000000..b0e1b707
--- /dev/null
+++ b/src/PJ_geos.c
@@ -0,0 +1,151 @@
+/*
+** libproj -- library of cartographic projections
+**
+** Copyright (c) 2004   Gerald I. Evenden
+*/
+static const char
+LIBPROJ_ID[] = "$Id$";
+/*
+** 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	h; \
+	double  radius_p; \
+	double  radius_p2; \
+	double  radius_p_inv2; \
+	double  radius_g; \
+	double  radius_g_1; \
+	double  C;
+#define PJ_LIB__
+#include	<projects.h>
+
+PROJ_HEAD(geos, "Geostationary Satellite View") "\n\tAzi, Sph&Ell\n\th=";
+
+FORWARD(s_forward); /* spheroid */
+	double Vx, Vy, Vz, tmp;
+
+/* Calculation of the three components of the vector from satellite to
+** position on earth surface (lon,lat).*/
+	tmp = cos(lp.phi);
+	Vx = cos (lp.lam) * tmp;
+	Vy = sin (lp.lam) * tmp;
+	Vz = sin (lp.phi);
+/* Check visibility.*/
+	if (((P->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz) < 0.) F_ERROR;
+/* Calculation based on view angles from satellite.*/
+	tmp = P->radius_g - Vx;
+	xy.x = P->radius_g_1 * atan(Vy / tmp);
+	xy.y = P->radius_g_1 * atan(Vz / hypot(Vy, tmp));
+	return (xy);
+}
+FORWARD(e_forward); /* ellipsoid */
+	double r, Vx, Vy, Vz, tmp;
+
+/* Calculation of geocentric latitude. */
+	lp.phi = atan (P->radius_p2 * tan (lp.phi));
+/* Calculation of the three components of the vector from satellite to
+** position on earth surface (lon,lat).*/
+	r = (P->radius_p) / hypot(P->radius_p * cos (lp.phi), sin (lp.phi));
+	Vx = r * cos (lp.lam) * cos (lp.phi);
+	Vy = r * sin (lp.lam) * cos (lp.phi);
+	Vz = r * sin (lp.phi);
+/* Check visibility. */
+	if (((P->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz * P->radius_p_inv2) < 0.)
+		F_ERROR;
+/* Calculation based on view angles from satellite. */
+	tmp = P->radius_g - Vx;
+	xy.x = P->radius_g_1 * atan (Vy / tmp);
+	xy.y = P->radius_g_1 * atan (Vz / hypot (Vy, tmp));
+	return (xy);
+}
+INVERSE(s_inverse); /* spheroid */
+	double Vx, Vy, Vz, a, b, c, det, k;
+
+/* Setting three components of vector from satellite to position.*/
+	Vx = -1.0;
+	Vy = tan (xy.x / (P->radius_g - 1.0));
+	Vz = tan (xy.y / (P->radius_g - 1.0)) * sqrt (1.0 + Vy * Vy);
+/* Calculation of terms in cubic equation and determinant.*/
+	a   = Vy * Vy + Vz * Vz + Vx * Vx;
+	b   = 2 * P->radius_g * Vx;
+	if ((det = (b * b) - 4 * a * P->C) < 0.) I_ERROR;
+/* Calculation of three components of vector from satellite to position.*/
+	k  = (-b - sqrt(det)) / (2 * a);
+	Vx = P->radius_g + k * Vx;
+	Vy *= k;
+	Vz *= k;
+/* Calculation of longitude and latitude.*/
+	lp.lam = atan2 (Vy, Vx);
+	lp.phi = atan (Vz * cos (lp.lam) / Vx);
+	return (lp);
+}
+INVERSE(e_inverse); /* ellipsoid */
+	double Vx, Vy, Vz, a, b, c, det, k;
+
+/* Setting three components of vector from satellite to position.*/
+	Vx = -1.0;
+	Vy = tan (xy.x / P->radius_g_1);
+	Vz = tan (xy.y / P->radius_g_1) * hypot(1.0, Vy);
+/* Calculation of terms in cubic equation and determinant.*/
+	a = Vz / P->radius_p;
+	a   = Vy * Vy + a * a + Vx * Vx;
+	b   = 2 * P->radius_g * Vx;
+	if ((det = (b * b) - 4 * a * P->C) < 0.) I_ERROR;
+/* Calculation of three components of vector from satellite to position.*/
+	k  = (-b - sqrt(det)) / (2. * a);
+	Vx = P->radius_g + k * Vx;
+	Vy *= k;
+	Vz *= k;
+/* Calculation of longitude and latitude.*/
+	lp.lam  = atan2 (Vy, Vx);
+	lp.phi = atan (Vz * cos (lp.lam) / Vx);
+	lp.phi = atan (P->radius_p_inv2 * tan (lp.phi));
+	return (lp);
+}
+FREEUP; if (P) free(P); }
+ENTRY0(geos)
+	if ((P->h = pj_param(P->params, "dh").f) <= 0.) E_ERROR(-30);
+	if (P->phi0) E_ERROR(-46);
+	P->radius_g = 1. + (P->radius_g_1 = P->h / P->a);
+	P->C  = P->radius_g * P->radius_g - 1.0;
+	if (P->es) {
+		P->radius_p      = sqrt (P->one_es);
+		P->radius_p2     = P->one_es;
+		P->radius_p_inv2 = P->rone_es;
+		P->inv = e_inverse;
+		P->fwd = e_forward;
+	} else {
+		P->radius_p = P->radius_p2 = P->radius_p_inv2 = 1.0;
+		P->inv = s_inverse;
+		P->fwd = s_forward;
+	}
+ENDENTRY(P)
+/*
+** $Log$
+** Revision 1.1  2004/10/20 17:04:00  fwarmerdam
+** New
+**
+** Revision 1.2  2004/07/14 18:08:57  gie
+** corrected P->phi_0 to P->phi0
+**
+** Revision 1.1  2004/07/12 17:58:25  gie
+** Initial revision
+**
+*/
diff --git a/src/PJ_sterea.c b/src/PJ_sterea.c
new file mode 100644
index 00000000..f1ab98fe
--- /dev/null
+++ b/src/PJ_sterea.c
@@ -0,0 +1,97 @@
+/*
+** libproj -- library of cartographic projections
+**
+** Copyright (c) 2003   Gerald I. Evenden
+*/
+static const char
+LIBPROJ_ID[] = "$Id$";
+/*
+** 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$
+** 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
+**
+*/
diff --git a/src/pj_gauss.c b/src/pj_gauss.c
new file mode 100644
index 00000000..f96d08fa
--- /dev/null
+++ b/src/pj_gauss.c
@@ -0,0 +1,107 @@
+/*
+** libproj -- library of cartographic projections
+**
+** Copyright (c) 2003   Gerald I. Evenden
+*/
+static const char
+LIBPROJ_ID[] = "$Id$";
+/*
+** 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 PJ_LIB__
+#include <projects.h>
+
+#define MAX_ITER 20
+
+struct GAUSS {
+	double C;
+	double K;
+	double e;
+	double ratexp;
+};
+#define EN ((struct GAUSS *)en)
+#define DEL_TOL 1e-14
+	static double
+srat(double esinp, double exp) {
+	return(pow((1.-esinp)/(1.+esinp), exp));
+}
+
+	void *
+pj_gauss_ini(double e, double phi0, double *chi, double *rc) {
+	double sphi, cphi, es;
+	struct GAUSS *en;
+
+	if ((en = (struct GAUSS *)malloc(sizeof(struct GAUSS))) == NULL)
+		return (NULL);
+	es = e * e;
+	EN->e = e;
+	sphi = sin(phi0);
+	cphi = cos(phi0);  cphi *= cphi;
+	*rc = sqrt(1. - es) / (1. - es * sphi * sphi);
+	EN->C = sqrt(1. + es * cphi * cphi / (1. - es));
+	*chi = asin(sphi / EN->C);
+	EN->ratexp = 0.5 * EN->C * e;
+	EN->K = tan(.5 * *chi + FORTPI) / (
+		pow(tan(.5 * phi0 + FORTPI), EN->C) *
+		srat(EN->e * sphi, EN->ratexp)  );
+	return ((void *)en);
+}
+	LP
+pj_gauss(LP elp, const void *en) {
+	LP slp;
+
+	slp.phi = 2. * atan( EN->K *
+		pow(tan(.5 * elp.phi + FORTPI), EN->C) *
+		srat(EN->e * sin(elp.phi), EN->ratexp) ) - HALFPI;
+	slp.lam = EN->C * (elp.lam);
+	return(slp);
+}
+	LP
+pj_inv_gauss(LP slp, const void *en) {
+	LP elp;
+	double num;
+	int i;
+
+	elp.lam = slp.lam / EN->C;
+	num = pow(tan(.5 * slp.phi + FORTPI)/EN->K, 1./EN->C);
+	for (i = MAX_ITER; i; --i) {
+		elp.phi = 2. * atan(num * srat(EN->e * sin(slp.phi), -.5 * EN->e))
+			- HALFPI;
+		if (fabs(elp.phi - slp.phi) < DEL_TOL) break;
+			slp.phi = elp.phi;
+	}	
+	/* convergence failed */
+	if (!i)
+		pj_errno = -17;
+	return (elp);
+}
+/* Revision Log:
+** $Log$
+** Revision 1.1  2004/10/20 17:04:00  fwarmerdam
+** New
+**
+** Revision 2.2  2004/03/15 16:07:42  gie
+** removed es from init structure
+**
+** Revision 2.1  2003/03/28 01:44:30  gie
+** Initial
+**
+*/
+