/Proj4/geod_inv.c

http://github.com/route-me/route-me · C · 59 lines · 58 code · 1 blank · 0 comment · 6 complexity · 61e5e07b2affee7b11b7124bf8808cd1 MD5 · raw file 

    

  1. #ifndef lint
    2. 

  2. static const char SCCSID[]="@(#)geod_inv.c	4.5	95/09/23	GIE	REL";
    3. 

  3. #endif
    4. 

  4. # include "projects.h"
    5. 

  5. # include "geodesic.h"
    6. 

  6. # define DTOL	1e-12
    7. 

  7. 	void
    8. 

  8. geod_inv(void) {
    9. 

  9. 	double	th1,th2,thm,dthm,dlamm,dlam,sindlamm,costhm,sinthm,cosdthm,
    10. 

  10. 		sindthm,L,E,cosd,d,X,Y,T,sind,tandlammp,u,v,D,A,B;
    11. 

  11. 

  12. 	if (ellipse) {
    13. 

  13. 		th1 = atan(onef * tan(phi1));
    14. 

  14. 		th2 = atan(onef * tan(phi2));
    15. 

  15. 	} else {
    16. 

  16. 		th1 = phi1;
    17. 

  17. 		th2 = phi2;
    18. 

  18. 	}
    19. 

  19. 	thm = .5 * (th1 + th2);
    20. 

  20. 	dthm = .5 * (th2 - th1);
    21. 

  21. 	dlamm = .5 * ( dlam = adjlon(lam2 - lam1) );
    22. 

  22. 	if (fabs(dlam) < DTOL && fabs(dthm) < DTOL) {
    23. 

  23. 		al12 =  al21 = geod_S = 0.;
    24. 

  24. 		return;
    25. 

  25. 	}
    26. 

  26. 	sindlamm = sin(dlamm);
    27. 

  27. 	costhm = cos(thm);	sinthm = sin(thm);
    28. 

  28. 	cosdthm = cos(dthm);	sindthm = sin(dthm);
    29. 

  29. 	L = sindthm * sindthm + (cosdthm * cosdthm - sinthm * sinthm)
    30. 

  30. 		* sindlamm * sindlamm;
    31. 

  31. 	d = acos(cosd = 1 - L - L);
    32. 

  32. 	if (ellipse) {
    33. 

  33. 		E = cosd + cosd;
    34. 

  34. 		sind = sin( d );
    35. 

  35. 		Y = sinthm * cosdthm;
    36. 

  36. 		Y *= (Y + Y) / (1. - L);
    37. 

  37. 		T = sindthm * costhm;
    38. 

  38. 		T *= (T + T) / L;
    39. 

  39. 		X = Y + T;
    40. 

  40. 		Y -= T;
    41. 

  41. 		T = d / sind;
    42. 

  42. 		D = 4. * T * T;
    43. 

  43. 		A = D * E;
    44. 

  44. 		B = D + D;
    45. 

  45. 		geod_S = geod_a * sind * (T - f4 * (T * X - Y) +
    46. 

  46. 			f64 * (X * (A + (T - .5 * (A - E)) * X) -
    47. 

  47. 			Y * (B + E * Y) + D * X * Y));
    48. 

  48. 		tandlammp = tan(.5 * (dlam - .25 * (Y + Y - E * (4. - X)) *
    49. 

  49. 			(f2 * T + f64 * (32. * T - (20. * T - A)
    50. 

  50. 			* X - (B + 4.) * Y)) * tan(dlam)));
    51. 

  51. 	} else {
    52. 

  52. 		geod_S = geod_a * d;
    53. 

  53. 		tandlammp = tan(dlamm);
    54. 

  54. 	}
    55. 

  55. 	u = atan2(sindthm , (tandlammp * costhm));
    56. 

  56. 	v = atan2(cosdthm , (tandlammp * sinthm));
    57. 

  57. 	al12 = adjlon(TWOPI + v - u);
    58. 

  58. 	al21 = adjlon(TWOPI - v - u);
    59. 

  59. }
    60.