/Proj4/PJ_lcca.c

http://github.com/route-me/route-me · C · 76 lines · 64 code · 6 blank · 6 comment · 9 complexity · 343d16310027e7fdf17c1c91c36adeb3 MD5 · raw file 

    

  1. static const char RCS_ID[] = "$Id: PJ_lcca.c,v 1.1 2003/03/04 02:59:41 warmerda Exp $";
    2. 

  2. /* PROJ.4 Cartographic Projection System -- Revision Log:
    3. 

  3. **$Log: PJ_lcca.c,v $
    4. 

  4. **Revision 1.1  2003/03/04 02:59:41  warmerda
    5. 

  5. **New
    6. 

  6. **
    7. 

  7. */
    8. 

  8. #define MAX_ITER 10
    9. 

  9. #define DEL_TOL 1e-12
    10. 

  10. #define PROJ_PARMS__ \
    11. 

  11. 	double	*en; \
    12. 

  12. 	double	r0, l, M0; \
    13. 

  13. 	double	C;
    14. 

  14. #define PJ_LIB__
    15. 

  15. #include	"projects.h"
    16. 

  16. 

  17. PROJ_HEAD(lcca, "Lambert Conformal Conic Alternative")
    18. 

  18. 	"\n\tConic, Sph&Ell\n\tlat_0=";
    19. 

  19. 

  20. 	static double /* func to compute dr */
    21. 

  21. fS(double S, double C) {
    22. 

  22. 		return(S * ( 1. + S * S * C));
    23. 

  23. }
    24. 

  24. 	static double /* deriv of fs */
    25. 

  25. fSp(double S, double C) {
    26. 

  26. 	return(1. + 3.* S * S * C);
    27. 

  27. }
    28. 

  28. FORWARD(e_forward); /* ellipsoid */
    29. 

  29. 	double S, S3, r, dr;
    30. 

  30. 	
    31. 

  31. 	S = pj_mlfn(lp.phi, sin(lp.phi), cos(lp.phi), P->en) - P->M0;
    32. 

  32. 	dr = fS(S, P->C);
    33. 

  33. 	r = P->r0 - dr;
    34. 

  34. 	xy.x = P->k0 * (r * sin( lp.lam *= P->l ) );
    35. 

  35. 	xy.y = P->k0 * (P->r0 - r * cos(lp.lam) );
    36. 

  36. 	return (xy);
    37. 

  37. }
    38. 

  38. INVERSE(e_inverse); /* ellipsoid & spheroid */
    39. 

  39. 	double theta, dr, S, dif;
    40. 

  40. 	int i;
    41. 

  41. 

  42. 	xy.x /= P->k0;
    43. 

  43. 	xy.y /= P->k0;
    44. 

  44. 	theta = atan2(xy.x , P->r0 - xy.y);
    45. 

  45. 	dr = xy.y - xy.x * tan(0.5 * theta);
    46. 

  46. 	lp.lam = theta / P->l;
    47. 

  47. 	S = dr;
    48. 

  48. 	for (i = MAX_ITER; i ; --i) {
    49. 

  49. 		S -= (dif = (fS(S, P->C) - dr) / fSp(S, P->C));
    50. 

  50. 		if (fabs(dif) < DEL_TOL) break;
    51. 

  51. 	}
    52. 

  52. 	if (!i) I_ERROR
    53. 

  53. 	lp.phi = pj_inv_mlfn(S + P->M0, P->es, P->en);
    54. 

  54. 	return (lp);
    55. 

  55. }
    56. 

  56. FREEUP; if (P) { if (P->en) pj_dalloc(P->en); pj_dalloc(P); } }
    57. 

  57. ENTRY0(lcca)
    58. 

  58. 	double s2p0, N0, R0, tan0, tan20;
    59. 

  59. 

  60. 	if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
    61. 

  61. 	if (!pj_param(P->params, "tlat_0").i) E_ERROR(50);
    62. 

  62. 	if (P->phi0 == 0.) E_ERROR(51);
    63. 

  63. 	P->l = sin(P->phi0);
    64. 

  64. 	P->M0 = pj_mlfn(P->phi0, P->l, cos(P->phi0), P->en);
    65. 

  65. 	s2p0 = P->l * P->l;
    66. 

  66. 	R0 = 1. / (1. - P->es * s2p0);
    67. 

  67. 	N0 = sqrt(R0);
    68. 

  68. 	R0 *= P->one_es * N0;
    69. 

  69. 	tan0 = tan(P->phi0);
    70. 

  70. 	tan20 = tan0 * tan0;
    71. 

  71. 	P->r0 = N0 / tan0;
    72. 

  72. 	P->C = 1. / (6. * R0 * N0);
    73. 

  73. 	P->inv = e_inverse;
    74. 

  74. 	P->fwd = e_forward;
    75. 

  75. ENDENTRY(P)
    76.