/Proj4/pj_gauss.c

  1/*
  2** libproj -- library of cartographic projections
  3**
  4** Copyright (c) 2003   Gerald I. Evenden
  5*/
  6static const char
  7LIBPROJ_ID[] = "$Id: pj_gauss.c,v 1.1 2004/10/20 17:04:00 fwarmerdam Exp $";
  8/*
  9** Permission is hereby granted, free of charge, to any person obtaining
 10** a copy of this software and associated documentation files (the
 11** "Software"), to deal in the Software without restriction, including
 12** without limitation the rights to use, copy, modify, merge, publish,
 13** distribute, sublicense, and/or sell copies of the Software, and to
 14** permit persons to whom the Software is furnished to do so, subject to
 15** the following conditions:
 16**
 17** The above copyright notice and this permission notice shall be
 18** included in all copies or substantial portions of the Software.
 19**
 20** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 21** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 22** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 23** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
 24** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
 25** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
 26** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 27*/
 28#define PJ_LIB__
 29#include "projects.h"
 30
 31#define MAX_ITER 20
 32
 33struct GAUSS {
 34	double C;
 35	double K;
 36	double e;
 37	double ratexp;
 38};
 39#define EN ((struct GAUSS *)en)
 40#define DEL_TOL 1e-14
 41	static double
 42srat(double esinp, double exp) {
 43	return(pow((1.-esinp)/(1.+esinp), exp));
 44}
 45
 46	void *
 47pj_gauss_ini(double e, double phi0, double *chi, double *rc) {
 48	double sphi, cphi, es;
 49	struct GAUSS *en;
 50
 51	if ((en = (struct GAUSS *)malloc(sizeof(struct GAUSS))) == NULL)
 52		return (NULL);
 53	es = e * e;
 54	EN->e = e;
 55	sphi = sin(phi0);
 56	cphi = cos(phi0);  cphi *= cphi;
 57	*rc = sqrt(1. - es) / (1. - es * sphi * sphi);
 58	EN->C = sqrt(1. + es * cphi * cphi / (1. - es));
 59	*chi = asin(sphi / EN->C);
 60	EN->ratexp = 0.5 * EN->C * e;
 61	EN->K = tan(.5 * *chi + FORTPI) / (
 62		pow(tan(.5 * phi0 + FORTPI), EN->C) *
 63		srat(EN->e * sphi, EN->ratexp)  );
 64	return ((void *)en);
 65}
 66	LP
 67pj_gauss(LP elp, const void *en) {
 68	LP slp;
 69
 70	slp.phi = 2. * atan( EN->K *
 71		pow(tan(.5 * elp.phi + FORTPI), EN->C) *
 72		srat(EN->e * sin(elp.phi), EN->ratexp) ) - HALFPI;
 73	slp.lam = EN->C * (elp.lam);
 74	return(slp);
 75}
 76	LP
 77pj_inv_gauss(LP slp, const void *en) {
 78	LP elp;
 79	double num;
 80	int i;
 81
 82	elp.lam = slp.lam / EN->C;
 83	num = pow(tan(.5 * slp.phi + FORTPI)/EN->K, 1./EN->C);
 84	for (i = MAX_ITER; i; --i) {
 85		elp.phi = 2. * atan(num * srat(EN->e * sin(slp.phi), -.5 * EN->e))
 86			- HALFPI;
 87		if (fabs(elp.phi - slp.phi) < DEL_TOL) break;
 88			slp.phi = elp.phi;
 89	}	
 90	/* convergence failed */
 91	if (!i)
 92		pj_errno = -17;
 93	return (elp);
 94}
 95/* Revision Log:
 96** $Log: pj_gauss.c,v $
 97** Revision 1.1  2004/10/20 17:04:00  fwarmerdam
 98** New
 99**
100** Revision 2.2  2004/03/15 16:07:42  gie
101** removed es from init structure
102**
103** Revision 2.1  2003/03/28 01:44:30  gie
104** Initial
105**
106*/
107