/Proj4/PJ_cea.c

http://github.com/route-me/route-me · C · 64 lines · 62 code · 2 blank · 0 comment · 9 complexity · abf061afab73939b2af8fdde43ada069 MD5 · raw file 

    

  1. #ifndef lint
    2. 

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

  3. #endif
    4. 

  4. #define PROJ_PARMS__ \
    5. 

  5. 	double qp; \
    6. 

  6. 	double *apa;
    7. 

  7. #define PJ_LIB__
    8. 

  8. # include	"projects.h"
    9. 

  9. PROJ_HEAD(cea, "Equal Area Cylindrical") "\n\tCyl, Sph&Ell\n\tlat_ts=";
    10. 

  10. # define EPS	1e-10
    11. 

  11. FORWARD(e_forward); /* spheroid */
    12. 

  12. 	xy.x = P->k0 * lp.lam;
    13. 

  13. 	xy.y = .5 * pj_qsfn(sin(lp.phi), P->e, P->one_es) / P->k0;
    14. 

  14. 	return (xy);
    15. 

  15. }
    16. 

  16. FORWARD(s_forward); /* spheroid */
    17. 

  17. 	xy.x = P->k0 * lp.lam;
    18. 

  18. 	xy.y = sin(lp.phi) / P->k0;
    19. 

  19. 	return (xy);
    20. 

  20. }
    21. 

  21. INVERSE(e_inverse); /* spheroid */
    22. 

  22. 	lp.phi = pj_authlat(asin( 2. * xy.y * P->k0 / P->qp), P->apa);
    23. 

  23. 	lp.lam = xy.x / P->k0;
    24. 

  24. 	return (lp);
    25. 

  25. }
    26. 

  26. INVERSE(s_inverse); /* spheroid */
    27. 

  27. 	double t;
    28. 

  28. 

  29. 	if ((t = fabs(xy.y *= P->k0)) - EPS <= 1.) {
    30. 

  30. 		if (t >= 1.)
    31. 

  31. 			lp.phi = xy.y < 0. ? -HALFPI : HALFPI;
    32. 

  32. 		else
    33. 

  33. 			lp.phi = asin(xy.y);
    34. 

  34. 		lp.lam = xy.x / P->k0;
    35. 

  35. 	} else I_ERROR;
    36. 

  36. 	return (lp);
    37. 

  37. }
    38. 

  38. FREEUP;
    39. 

  39. 	if (P) {
    40. 

  40. 		if (P->apa)
    41. 

  41. 			pj_dalloc(P->apa);
    42. 

  42. 		pj_dalloc(P);
    43. 

  43. 	}
    44. 

  44. }
    45. 

  45. ENTRY1(cea, apa)
    46. 

  46. 	double t;
    47. 

  47. 

  48. 	if (pj_param(P->params, "tlat_ts").i &&
    49. 

  49. 		(P->k0 = cos(t = pj_param(P->params, "rlat_ts").f)) < 0.) E_ERROR(-24)
    50. 

  50. 	else
    51. 

  51. 		t = 0.;
    52. 

  52. 	if (P->es) {
    53. 

  53. 		t = sin(t);
    54. 

  54. 		P->k0 /= sqrt(1. - P->es * t * t);
    55. 

  55. 		P->e = sqrt(P->es);
    56. 

  56. 		if (!(P->apa = pj_authset(P->es))) E_ERROR_0;
    57. 

  57. 		P->qp = pj_qsfn(1., P->e, P->one_es);
    58. 

  58. 		P->inv = e_inverse;
    59. 

  59. 		P->fwd = e_forward;
    60. 

  60. 	} else {
    61. 

  61. 		P->inv = s_inverse;
    62. 

  62. 		P->fwd = s_forward;
    63. 

  63. 	}
    64. 

  64. ENDENTRY(P)
    65.