/Proj4/pj_mlfn.c

http://github.com/route-me/route-me · C · 61 lines · 54 code · 2 blank · 5 comment · 3 complexity · 0a33f005f884f2ca1a56a9b9ece4d2ca MD5 · raw file 

    

  1. #ifndef lint
    2. 

  2. static const char SCCSID[]="@(#)pj_mlfn.c	4.5	95/07/06	GIE	REL";
    3. 

  3. #endif
    4. 

  4. #include "projects.h"
    5. 

  5. /* meridinal distance for ellipsoid and inverse
    6. 

  6. **	8th degree - accurate to < 1e-5 meters when used in conjuction
    7. 

  7. **		with typical major axis values.
    8. 

  8. **	Inverse determines phi to EPS (1e-11) radians, about 1e-6 seconds.
    9. 

  9. */
    10. 

  10. #define C00 1.
    11. 

  11. #define C02 .25
    12. 

  12. #define C04 .046875
    13. 

  13. #define C06 .01953125
    14. 

  14. #define C08 .01068115234375
    15. 

  15. #define C22 .75
    16. 

  16. #define C44 .46875
    17. 

  17. #define C46 .01302083333333333333
    18. 

  18. #define C48 .00712076822916666666
    19. 

  19. #define C66 .36458333333333333333
    20. 

  20. #define C68 .00569661458333333333
    21. 

  21. #define C88 .3076171875
    22. 

  22. #define EPS 1e-11
    23. 

  23. #define MAX_ITER 10
    24. 

  24. #define EN_SIZE 5
    25. 

  25. 	double *
    26. 

  26. pj_enfn(double es) {
    27. 

  27. 	double t, *en;
    28. 

  28. 

  29. 	en = (double *)pj_malloc(EN_SIZE * sizeof(double));
    30. 

  30. 	if (en) {
    31. 

  31. 		en[0] = C00 - es * (C02 + es * (C04 + es * (C06 + es * C08)));
    32. 

  32. 		en[1] = es * (C22 - es * (C04 + es * (C06 + es * C08)));
    33. 

  33. 		en[2] = (t = es * es) * (C44 - es * (C46 + es * C48));
    34. 

  34. 		en[3] = (t *= es) * (C66 - es * C68);
    35. 

  35. 		en[4] = t * es * C88;
    36. 

  36. 	} /* else return NULL if unable to allocate memory */
    37. 

  37. 	return en;
    38. 

  38. }
    39. 

  39. 	double
    40. 

  40. pj_mlfn(double phi, double sphi, double cphi, double *en) {
    41. 

  41. 	cphi *= sphi;
    42. 

  42. 	sphi *= sphi;
    43. 

  43. 	return(en[0] * phi - cphi * (en[1] + sphi*(en[2]
    44. 

  44. 		+ sphi*(en[3] + sphi*en[4]))));
    45. 

  45. }
    46. 

  46. 	double
    47. 

  47. pj_inv_mlfn(double arg, double es, double *en) {
    48. 

  48. 	double s, t, phi, k = 1./(1.-es);
    49. 

  49. 	int i;
    50. 

  50. 

  51. 	phi = arg;
    52. 

  52. 	for (i = MAX_ITER; i ; --i) { /* rarely goes over 2 iterations */
    53. 

  53. 		s = sin(phi);
    54. 

  54. 		t = 1. - es * s * s;
    55. 

  55. 		phi -= t = (pj_mlfn(phi, s, cos(phi), en) - arg) * (t * sqrt(t)) * k;
    56. 

  56. 		if (fabs(t) < EPS)
    57. 

  57. 			return phi;
    58. 

  58. 	}
    59. 

  59. 	pj_errno = -17;
    60. 

  60. 	return phi;
    61. 

  61. }
    62.