PageRenderTime 58ms CodeModel.GetById 24ms app.highlight 29ms RepoModel.GetById 1ms app.codeStats 0ms

/Proj4/pj_ell_set.c

http://github.com/route-me/route-me
C | 107 lines | 97 code | 4 blank | 6 comment | 37 complexity | 168563f4bd08b883dfc0c6ea4f079261 MD5 | raw file
  1/* set ellipsoid parameters a and es */
  2#ifndef lint
  3static const char SCCSID[]="@(#)pj_ell_set.c	4.5	93/06/12	GIE	REL";
  4#endif
  5#include "projects.h"
  6#include <string.h>
  7#define SIXTH .1666666666666666667 /* 1/6 */
  8#define RA4 .04722222222222222222 /* 17/360 */
  9#define RA6 .02215608465608465608 /* 67/3024 */
 10#define RV4 .06944444444444444444 /* 5/72 */
 11#define RV6 .04243827160493827160 /* 55/1296 */
 12	int /* initialize geographic shape parameters */
 13pj_ell_set(paralist *pl, double *a, double *es) {
 14	int i;
 15	double b=0.0, e;
 16	char *name;
 17	paralist *start = 0, *curr;
 18
 19		/* check for varying forms of ellipsoid input */
 20	*a = *es = 0.;
 21	/* R takes precedence */
 22	if (pj_param(pl, "tR").i)
 23		*a = pj_param(pl, "dR").f;
 24	else { /* probable elliptical figure */
 25
 26		/* check if ellps present and temporarily append its values to pl */
 27		name = pj_param(pl, "sellps").s;
 28		if (name && pl) {
 29			char *s;
 30
 31			for (start = pl; start && start->next ; start = start->next) ;
 32			curr = start;
 33			for (i = 0; (s = pj_ellps[i].id) && strcmp(name, s) ; ++i) ;
 34			if (!s) { pj_errno = -9; return 1; }
 35			curr->next = pj_mkparam(pj_ellps[i].major);
 36			curr = curr->next;
 37			curr->next = pj_mkparam(pj_ellps[i].ell);
 38		}
 39		*a = pj_param(pl, "da").f;
 40		if (pj_param(pl, "tes").i) /* eccentricity squared */
 41			*es = pj_param(pl, "des").f;
 42		else if (pj_param(pl, "te").i) { /* eccentricity */
 43			e = pj_param(pl, "de").f;
 44			*es = e * e;
 45		} else if (pj_param(pl, "trf").i) { /* recip flattening */
 46			*es = pj_param(pl, "drf").f;
 47			if (!*es) {
 48				pj_errno = -10;
 49				goto bomb;
 50			}
 51			*es = 1./ *es;
 52			*es = *es * (2. - *es);
 53		} else if (pj_param(pl, "tf").i) { /* flattening */
 54			*es = pj_param(pl, "df").f;
 55			*es = *es * (2. - *es);
 56		} else if (pj_param(pl, "tb").i) { /* minor axis */
 57			b = pj_param(pl, "db").f;
 58			*es = 1. - (b * b) / (*a * *a);
 59		}     /* else *es == 0. and sphere of radius *a */
 60		if (!b)
 61			b = *a * sqrt(1. - *es);
 62		/* following options turn ellipsoid into equivalent sphere */
 63		if (pj_param(pl, "bR_A").i) { /* sphere--area of ellipsoid */
 64			*a *= 1. - *es * (SIXTH + *es * (RA4 + *es * RA6));
 65			*es = 0.;
 66		} else if (pj_param(pl, "bR_V").i) { /* sphere--vol. of ellipsoid */
 67			*a *= 1. - *es * (SIXTH + *es * (RV4 + *es * RV6));
 68			*es = 0.;
 69		} else if (pj_param(pl, "bR_a").i) { /* sphere--arithmetic mean */
 70			*a = .5 * (*a + b);
 71			*es = 0.;
 72		} else if (pj_param(pl, "bR_g").i) { /* sphere--geometric mean */
 73			*a = sqrt(*a * b);
 74			*es = 0.;
 75		} else if (pj_param(pl, "bR_h").i) { /* sphere--harmonic mean */
 76			*a = 2. * *a * b / (*a + b);
 77			*es = 0.;
 78		} else if ((i = pj_param(pl, "tR_lat_a").i) || /* sphere--arith. */
 79			pj_param(pl, "tR_lat_g").i) { /* or geom. mean at latitude */
 80			double tmp;
 81
 82			tmp = sin(pj_param(pl, i ? "rR_lat_a" : "rR_lat_g").f);
 83			if (fabs(tmp) > HALFPI) {
 84				pj_errno = -11;
 85				goto bomb;
 86			}
 87			tmp = 1. - *es * tmp * tmp;
 88			*a *= i ? .5 * (1. - *es + tmp) / ( tmp * sqrt(tmp)) :
 89				sqrt(1. - *es) / tmp;
 90			*es = 0.;
 91		}
 92bomb:
 93		if (start) { /* clean up temporary extension of list */
 94			pj_dalloc(start->next->next);
 95			pj_dalloc(start->next);
 96			start->next = 0;
 97		}
 98		if (pj_errno)
 99			return 1;
100	}
101	/* some remaining checks */
102	if (*es < 0.)
103		{ pj_errno = -12; return 1; }
104	if (*a <= 0.)
105		{ pj_errno = -13; return 1; }
106	return 0;
107}