/Proj4/geod_for.c

  1#ifndef lint
  2static const char SCCSID[]="@(#)geod_for.c	4.6	95/09/23	GIE	REL";
  3#endif
  4# include "projects.h"
  5# include "geodesic.h"
  6# define MERI_TOL 1e-9
  7	static double
  8th1,costh1,sinth1,sina12,cosa12,M,N,c1,c2,D,P,s1;
  9	static int
 10merid, signS;
 11	void
 12geod_pre(void) {
 13	al12 = adjlon(al12); /* reduce to  +- 0-PI */
 14	signS = fabs(al12) > HALFPI ? 1 : 0;
 15	th1 = ellipse ? atan(onef * tan(phi1)) : phi1;
 16	costh1 = cos(th1);
 17	sinth1 = sin(th1);
 18	if ((merid = fabs(sina12 = sin(al12)) < MERI_TOL)) {
 19		sina12 = 0.;
 20		cosa12 = fabs(al12) < HALFPI ? 1. : -1.;
 21		M = 0.;
 22	} else {
 23		cosa12 = cos(al12);
 24		M = costh1 * sina12;
 25	}
 26	N = costh1 * cosa12;
 27	if (ellipse) {
 28		if (merid) {
 29			c1 = 0.;
 30			c2 = f4;
 31			D = 1. - c2;
 32			D *= D;
 33			P = c2 / D;
 34		} else {
 35			c1 = geod_f * M;
 36			c2 = f4 * (1. - M * M);
 37			D = (1. - c2)*(1. - c2 - c1 * M);
 38			P = (1. + .5 * c1 * M) * c2 / D;
 39		}
 40	}
 41	if (merid) s1 = HALFPI - th1;
 42	else {
 43		s1 = (fabs(M) >= 1.) ? 0. : acos(M);
 44		s1 =  sinth1 / sin(s1);
 45		s1 = (fabs(s1) >= 1.) ? 0. : acos(s1);
 46	}
 47}
 48	void
 49geod_for(void) {
 50	double d,sind,u,V,X,ds,cosds,sinds,ss,de;
 51	ss = 0.;
 52	
 53	if (ellipse) {
 54		d = geod_S / (D * geod_a);
 55		if (signS) d = -d;
 56		u = 2. * (s1 - d);
 57		V = cos(u + d);
 58		X = c2 * c2 * (sind = sin(d)) * cos(d) * (2. * V * V - 1.);
 59		ds = d + X - 2. * P * V * (1. - 2. * P * cos(u)) * sind;
 60		ss = s1 + s1 - ds;
 61	} else {
 62		ds = geod_S / geod_a;
 63		if (signS) ds = - ds;
 64	}
 65	cosds = cos(ds);
 66	sinds = sin(ds);
 67	if (signS) sinds = - sinds;
 68	al21 = N * cosds - sinth1 * sinds;
 69	if (merid) {
 70		phi2 = atan( tan(HALFPI + s1 - ds) / onef);
 71		if (al21 > 0.) {
 72			al21 = PI;
 73			if (signS)
 74				de = PI;
 75			else {
 76				phi2 = - phi2;
 77				de = 0.;
 78			}
 79		} else {
 80			al21 = 0.;
 81			if (signS) {
 82				phi2 = - phi2;
 83				de = 0;
 84			} else
 85				de = PI;
 86		}
 87	} else {
 88		al21 = atan(M / al21);
 89		if (al21 > 0)
 90			al21 += PI;
 91		if (al12 < 0.)
 92			al21 -= PI;
 93		al21 = adjlon(al21);
 94		phi2 = atan(-(sinth1 * cosds + N * sinds) * sin(al21) /
 95			(ellipse ? onef * M : M));
 96		de = atan2(sinds * sina12 ,
 97			(costh1 * cosds - sinth1 * sinds * cosa12));
 98		if (ellipse)
 99			if (signS)
100				de += c1 * ((1. - c2) * ds +
101					c2 * sinds * cos(ss));
102			else
103				de -= c1 * ((1. - c2) * ds -
104					c2 * sinds * cos(ss));
105	}
106	lam2 = adjlon( lam1 + de );
107}