/Proj4/PJ_robin.c

http://github.com/route-me/route-me · C · 108 lines · 102 code · 2 blank · 4 comment · 15 complexity · c682333968ce57df501c64420025f135 MD5 · raw file 

    

  1. #ifndef lint
    2. 

  2. static const char SCCSID[]="@(#)PJ_robin.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(robin, "Robinson") "\n\tPCyl., Sph.";
    7. 

  7. #define V(C,z) (C.c0 + z * (C.c1 + z * (C.c2 + z * C.c3)))
    8. 

  8. #define DV(C,z) (C.c1 + z * (C.c2 + C.c2 + z * 3. * C.c3))
    9. 

  9. /* note: following terms based upon 5 deg. intervals in degrees. */
    10. 

  10. static struct COEFS {
    11. 

  11. 	float c0, c1, c2, c3;
    12. 

  12. } X[] = {
    13. 

  13. 1,	-5.67239e-12,	-7.15511e-05,	3.11028e-06,
    14. 

  14. 0.9986,	-0.000482241,	-2.4897e-05,	-1.33094e-06,
    15. 

  15. 0.9954,	-0.000831031,	-4.4861e-05,	-9.86588e-07,
    16. 

  16. 0.99,	-0.00135363,	-5.96598e-05,	3.67749e-06,
    17. 

  17. 0.9822,	-0.00167442,	-4.4975e-06,	-5.72394e-06,
    18. 

  18. 0.973,	-0.00214869,	-9.03565e-05,	1.88767e-08,
    19. 

  19. 0.96,	-0.00305084,	-9.00732e-05,	1.64869e-06,
    20. 

  20. 0.9427,	-0.00382792,	-6.53428e-05,	-2.61493e-06,
    21. 

  21. 0.9216,	-0.00467747,	-0.000104566,	4.8122e-06,
    22. 

  22. 0.8962,	-0.00536222,	-3.23834e-05,	-5.43445e-06,
    23. 

  23. 0.8679,	-0.00609364,	-0.0001139,	3.32521e-06,
    24. 

  24. 0.835,	-0.00698325,	-6.40219e-05,	9.34582e-07,
    25. 

  25. 0.7986,	-0.00755337,	-5.00038e-05,	9.35532e-07,
    26. 

  26. 0.7597,	-0.00798325,	-3.59716e-05,	-2.27604e-06,
    27. 

  27. 0.7186,	-0.00851366,	-7.0112e-05,	-8.63072e-06,
    28. 

  28. 0.6732,	-0.00986209,	-0.000199572,	1.91978e-05,
    29. 

  29. 0.6213,	-0.010418,	8.83948e-05,	6.24031e-06,
    30. 

  30. 0.5722,	-0.00906601,	0.000181999,	6.24033e-06,
    31. 

  31. 0.5322, 0.,0.,0.  },
    32. 

  32. Y[] = {
    33. 

  33. 0,	0.0124,	3.72529e-10,	1.15484e-09,
    34. 

  34. 0.062,	0.0124001,	1.76951e-08,	-5.92321e-09,
    35. 

  35. 0.124,	0.0123998,	-7.09668e-08,	2.25753e-08,
    36. 

  36. 0.186,	0.0124008,	2.66917e-07,	-8.44523e-08,
    37. 

  37. 0.248,	0.0123971,	-9.99682e-07,	3.15569e-07,
    38. 

  38. 0.31,	0.0124108,	3.73349e-06,	-1.1779e-06,
    39. 

  39. 0.372,	0.0123598,	-1.3935e-05,	4.39588e-06,
    40. 

  40. 0.434,	0.0125501,	5.20034e-05,	-1.00051e-05,
    41. 

  41. 0.4968,	0.0123198,	-9.80735e-05,	9.22397e-06,
    42. 

  42. 0.5571,	0.0120308,	4.02857e-05,	-5.2901e-06,
    43. 

  43. 0.6176,	0.0120369,	-3.90662e-05,	7.36117e-07,
    44. 

  44. 0.6769,	0.0117015,	-2.80246e-05,	-8.54283e-07,
    45. 

  45. 0.7346,	0.0113572,	-4.08389e-05,	-5.18524e-07,
    46. 

  46. 0.7903,	0.0109099,	-4.86169e-05,	-1.0718e-06,
    47. 

  47. 0.8435,	0.0103433,	-6.46934e-05,	5.36384e-09,
    48. 

  48. 0.8936,	0.00969679,	-6.46129e-05,	-8.54894e-06,
    49. 

  49. 0.9394,	0.00840949,	-0.000192847,	-4.21023e-06,
    50. 

  50. 0.9761,	0.00616525,	-0.000256001,	-4.21021e-06,
    51. 

  51. 1., 0.,0.,0 };
    52. 

  52. #define FXC	0.8487
    53. 

  53. #define FYC	1.3523
    54. 

  54. #define C1	11.45915590261646417544
    55. 

  55. #define RC1	0.08726646259971647884
    56. 

  56. #define NODES	18
    57. 

  57. #define ONEEPS	1.000001
    58. 

  58. #define EPS	1e-8
    59. 

  59. FORWARD(s_forward); /* spheroid */
    60. 

  60. 	int i;
    61. 

  61. 	double dphi;
    62. 

  62. 

  63. 	i = floor((dphi = fabs(lp.phi)) * C1);
    64. 

  64. 	if (i >= NODES) i = NODES - 1;
    65. 

  65. 	dphi = RAD_TO_DEG * (dphi - RC1 * i);
    66. 

  66. 	xy.x = V(X[i], dphi) * FXC * lp.lam;
    67. 

  67. 	xy.y = V(Y[i], dphi) * FYC;
    68. 

  68. 	if (lp.phi < 0.) xy.y = -xy.y;
    69. 

  69. 	return (xy);
    70. 

  70. }
    71. 

  71. INVERSE(s_inverse); /* spheroid */
    72. 

  72. 	int i;
    73. 

  73. 	double t, t1;
    74. 

  74. 	struct COEFS T;
    75. 

  75. 

  76. 	lp.lam = xy.x / FXC;
    77. 

  77. 	lp.phi = fabs(xy.y / FYC);
    78. 

  78. 	if (lp.phi >= 1.) { /* simple pathologic cases */
    79. 

  79. 		if (lp.phi > ONEEPS) I_ERROR
    80. 

  80. 		else {
    81. 

  81. 			lp.phi = xy.y < 0. ? -HALFPI : HALFPI;
    82. 

  82. 			lp.lam /= X[NODES].c0;
    83. 

  83. 		}
    84. 

  84. 	} else { /* general problem */
    85. 

  85. 		/* in Y space, reduce to table interval */
    86. 

  86. 		for (i = floor(lp.phi * NODES);;) {
    87. 

  87. 			if (Y[i].c0 > lp.phi) --i;
    88. 

  88. 			else if (Y[i+1].c0 <= lp.phi) ++i;
    89. 

  89. 			else break;
    90. 

  90. 		}
    91. 

  91. 		T = Y[i];
    92. 

  92. 		/* first guess, linear interp */
    93. 

  93. 		t = 5. * (lp.phi - T.c0)/(Y[i+1].c0 - T.c0);
    94. 

  94. 		/* make into root */
    95. 

  95. 		T.c0 -= lp.phi;
    96. 

  96. 		for (;;) { /* Newton-Raphson reduction */
    97. 

  97. 			t -= t1 = V(T,t) / DV(T,t);
    98. 

  98. 			if (fabs(t1) < EPS)
    99. 

  99. 				break;
    100. 

  100. 		}
    101. 

  101. 		lp.phi = (5 * i + t) * DEG_TO_RAD;
    102. 

  102. 		if (xy.y < 0.) lp.phi = -lp.phi;
    103. 

  103. 		lp.lam /= V(X[i], t);
    104. 

  104. 	}
    105. 

  105. 	return (lp);
    106. 

  106. }
    107. 

  107. FREEUP; if (P) pj_dalloc(P); }
    108. 

  108. ENTRY0(robin) P->es = 0.; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)
    109.