/Proj4/geod_inv.c

 1#ifndef lint
 2static const char SCCSID[]="@(#)geod_inv.c	4.5	95/09/23	GIE	REL";
 3#endif
 4# include "projects.h"
 5# include "geodesic.h"
 6# define DTOL	1e-12
 7	void
 8geod_inv(void) {
 9	double	th1,th2,thm,dthm,dlamm,dlam,sindlamm,costhm,sinthm,cosdthm,
10		sindthm,L,E,cosd,d,X,Y,T,sind,tandlammp,u,v,D,A,B;
11
12	if (ellipse) {
13		th1 = atan(onef * tan(phi1));
14		th2 = atan(onef * tan(phi2));
15	} else {
16		th1 = phi1;
17		th2 = phi2;
18	}
19	thm = .5 * (th1 + th2);
20	dthm = .5 * (th2 - th1);
21	dlamm = .5 * ( dlam = adjlon(lam2 - lam1) );
22	if (fabs(dlam) < DTOL && fabs(dthm) < DTOL) {
23		al12 =  al21 = geod_S = 0.;
24		return;
25	}
26	sindlamm = sin(dlamm);
27	costhm = cos(thm);	sinthm = sin(thm);
28	cosdthm = cos(dthm);	sindthm = sin(dthm);
29	L = sindthm * sindthm + (cosdthm * cosdthm - sinthm * sinthm)
30		* sindlamm * sindlamm;
31	d = acos(cosd = 1 - L - L);
32	if (ellipse) {
33		E = cosd + cosd;
34		sind = sin( d );
35		Y = sinthm * cosdthm;
36		Y *= (Y + Y) / (1. - L);
37		T = sindthm * costhm;
38		T *= (T + T) / L;
39		X = Y + T;
40		Y -= T;
41		T = d / sind;
42		D = 4. * T * T;
43		A = D * E;
44		B = D + D;
45		geod_S = geod_a * sind * (T - f4 * (T * X - Y) +
46			f64 * (X * (A + (T - .5 * (A - E)) * X) -
47			Y * (B + E * Y) + D * X * Y));
48		tandlammp = tan(.5 * (dlam - .25 * (Y + Y - E * (4. - X)) *
49			(f2 * T + f64 * (32. * T - (20. * T - A)
50			* X - (B + 4.) * Y)) * tan(dlam)));
51	} else {
52		geod_S = geod_a * d;
53		tandlammp = tan(dlamm);
54	}
55	u = atan2(sindthm , (tandlammp * costhm));
56	v = atan2(cosdthm , (tandlammp * sinthm));
57	al12 = adjlon(TWOPI + v - u);
58	al21 = adjlon(TWOPI - v - u);
59}