/Proj4/PJ_hatano.c

http://github.com/route-me/route-me · C · 53 lines · 51 code · 2 blank · 0 comment · 7 complexity · 81c232be867fa58ca05ed206d20a905d MD5 · raw file 

    

  1. #ifndef lint
    2. 

  2. static const char SCCSID[]="@(#)PJ_hatano.c	4.1 94/02/15     GIE     REL";
    3. 

  3. #endif
    4. 

  4. #define PJ_LIB__
    5. 

  5. #include	"projects.h"
    6. 

  6. PROJ_HEAD(hatano, "Hatano Asymmetrical Equal Area") "\n\tPCyl, Sph.";
    7. 

  7. #define NITER	20
    8. 

  8. #define EPS	1e-7
    9. 

  9. #define ONETOL 1.000001
    10. 

  10. #define CN	2.67595
    11. 

  11. #define CS	2.43763
    12. 

  12. #define RCN	0.37369906014686373063
    13. 

  13. #define RCS	0.41023453108141924738
    14. 

  14. #define FYCN	1.75859
    15. 

  15. #define FYCS	1.93052
    16. 

  16. #define RYCN	0.56863737426006061674
    17. 

  17. #define RYCS	0.51799515156538134803
    18. 

  18. #define FXC	0.85
    19. 

  19. #define RXC	1.17647058823529411764
    20. 

  20. FORWARD(s_forward); /* spheroid */
    21. 

  21. 	double th1, c;
    22. 

  22. 	int i;
    23. 

  23. 

  24. 	c = sin(lp.phi) * (lp.phi < 0. ? CS : CN);
    25. 

  25. 	for (i = NITER; i; --i) {
    26. 

  26. 		lp.phi -= th1 = (lp.phi + sin(lp.phi) - c) / (1. + cos(lp.phi));
    27. 

  27. 		if (fabs(th1) < EPS) break;
    28. 

  28. 	}
    29. 

  29. 	xy.x = FXC * lp.lam * cos(lp.phi *= .5);
    30. 

  30. 	xy.y = sin(lp.phi) * (lp.phi < 0. ? FYCS : FYCN);
    31. 

  31. 	return (xy);
    32. 

  32. }
    33. 

  33. INVERSE(s_inverse); /* spheroid */
    34. 

  34. 	double th;
    35. 

  35. 

  36. 	th = xy.y * ( xy.y < 0. ? RYCS : RYCN);
    37. 

  37. 	if (fabs(th) > 1.)
    38. 

  38. 		if (fabs(th) > ONETOL)	I_ERROR
    39. 

  39. 		else			th = th > 0. ? HALFPI : - HALFPI;
    40. 

  40. 	else
    41. 

  41. 		th = asin(th);
    42. 

  42. 	lp.lam = RXC * xy.x / cos(th);
    43. 

  43. 	th += th;
    44. 

  44. 	lp.phi = (th + sin(th)) * (xy.y < 0. ? RCS : RCN);
    45. 

  45. 	if (fabs(lp.phi) > 1.)
    46. 

  46. 		if (fabs(lp.phi) > ONETOL)	I_ERROR
    47. 

  47. 		else			lp.phi = lp.phi > 0. ? HALFPI : - HALFPI;
    48. 

  48. 	else
    49. 

  49. 		lp.phi = asin(lp.phi);
    50. 

  50. 	return (lp);
    51. 

  51. }
    52. 

  52. FREEUP; if (P) pj_dalloc(P); }
    53. 

  53. ENTRY0(hatano) P->es = 0.; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)
    54.