/Proj4/pj_factors.c

 1/* projection scale factors */
 2#ifndef lint
 3static const char SCCSID[]="@(#)pj_factors.c	4.9	94/03/17	GIE	REL";
 4#endif
 5#define PJ_LIB__
 6#include "projects.h"
 7#include <errno.h>
 8#ifndef DEFAULT_H
 9#define DEFAULT_H   1e-5    /* radian default for numeric h */
10#endif
11#define EPS 1.0e-12
12	int
13pj_factors(LP lp, PJ *P, double h, struct FACTORS *fac) {
14	struct DERIVS der;
15	double cosphi, t, n, r;
16
17	/* check for forward and latitude or longitude overange */
18	t = fabs(lp.phi)-HALFPI;
19	if (t > EPS || fabs(lp.lam) > 10.) {
20		pj_errno = -14;
21		return 1;
22	} else { /* proceed */
23		errno = pj_errno = 0;
24		if (h < EPS)
25			h = DEFAULT_H;
26		if (fabs(lp.phi) > (HALFPI - h)) 
27                /* adjust to value around pi/2 where derived still exists*/
28		        lp.phi = lp.phi < 0. ? (-HALFPI+h) : (HALFPI-h);
29		else if (P->geoc)
30			lp.phi = atan(P->rone_es * tan(lp.phi));
31		lp.lam -= P->lam0;	/* compute del lp.lam */
32		if (!P->over)
33			lp.lam = adjlon(lp.lam); /* adjust del longitude */
34		if (P->spc)	/* get what projection analytic values */
35			P->spc(lp, P, fac);
36		if (((fac->code & (IS_ANAL_XL_YL+IS_ANAL_XP_YP)) !=
37			  (IS_ANAL_XL_YL+IS_ANAL_XP_YP)) &&
38			  pj_deriv(lp, h, P, &der))
39			return 1;
40		if (!(fac->code & IS_ANAL_XL_YL)) {
41			fac->der.x_l = der.x_l;
42			fac->der.y_l = der.y_l;
43		}
44		if (!(fac->code & IS_ANAL_XP_YP)) {
45			fac->der.x_p = der.x_p;
46			fac->der.y_p = der.y_p;
47		}
48		cosphi = cos(lp.phi);
49		if (!(fac->code & IS_ANAL_HK)) {
50			fac->h = hypot(fac->der.x_p, fac->der.y_p);
51			fac->k = hypot(fac->der.x_l, fac->der.y_l) / cosphi;
52			if (P->es) {
53				t = sin(lp.phi);
54				t = 1. - P->es * t * t;
55				n = sqrt(t);
56				fac->h *= t * n / P->one_es;
57				fac->k *= n;
58				r = t * t / P->one_es;
59			} else
60				r = 1.;
61		} else if (P->es) {
62			r = sin(lp.phi);
63			r = 1. - P->es * r * r;
64			r = r * r / P->one_es;
65		} else
66			r = 1.;
67		/* convergence */
68		if (!(fac->code & IS_ANAL_CONV)) {
69			fac->conv = - atan2(fac->der.y_l, fac->der.x_l);
70			if (fac->code & IS_ANAL_XL_YL)
71				fac->code |= IS_ANAL_CONV;
72		}
73		/* areal scale factor */
74		fac->s = (fac->der.y_p * fac->der.x_l - fac->der.x_p * fac->der.y_l) *
75			r / cosphi;
76		/* meridian-parallel angle theta prime */
77		fac->thetap = aasin(fac->s / (fac->h * fac->k));
78		/* Tissot ellips axis */
79		t = fac->k * fac->k + fac->h * fac->h;
80		fac->a = sqrt(t + 2. * fac->s);
81		t = (t = t - 2. * fac->s) <= 0. ? 0. : sqrt(t);
82		fac->b = 0.5 * (fac->a - t);
83		fac->a = 0.5 * (fac->a + t);
84		/* omega */
85		fac->omega = 2. * aasin((fac->a - fac->b)/(fac->a + fac->b));
86	}
87	return 0;
88}