/Proj4/pj_gauss.c

http://github.com/route-me/route-me · C · 107 lines · 62 code · 7 blank · 38 comment · 5 complexity · 6664bdff4c492a935f3b4b40da7ef1b1 MD5 · raw file

  1. /*
  2. ** libproj -- library of cartographic projections
  3. **
  4. ** Copyright (c) 2003 Gerald I. Evenden
  5. */
  6. static const char
  7. LIBPROJ_ID[] = "$Id: pj_gauss.c,v 1.1 2004/10/20 17:04:00 fwarmerdam Exp $";
  8. /*
  9. ** Permission is hereby granted, free of charge, to any person obtaining
  10. ** a copy of this software and associated documentation files (the
  11. ** "Software"), to deal in the Software without restriction, including
  12. ** without limitation the rights to use, copy, modify, merge, publish,
  13. ** distribute, sublicense, and/or sell copies of the Software, and to
  14. ** permit persons to whom the Software is furnished to do so, subject to
  15. ** the following conditions:
  16. **
  17. ** The above copyright notice and this permission notice shall be
  18. ** included in all copies or substantial portions of the Software.
  19. **
  20. ** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  21. ** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  22. ** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
  23. ** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
  24. ** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
  25. ** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
  26. ** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  27. */
  28. #define PJ_LIB__
  29. #include "projects.h"
  30. #define MAX_ITER 20
  31. struct GAUSS {
  32. double C;
  33. double K;
  34. double e;
  35. double ratexp;
  36. };
  37. #define EN ((struct GAUSS *)en)
  38. #define DEL_TOL 1e-14
  39. static double
  40. srat(double esinp, double exp) {
  41. return(pow((1.-esinp)/(1.+esinp), exp));
  42. }
  43. void *
  44. pj_gauss_ini(double e, double phi0, double *chi, double *rc) {
  45. double sphi, cphi, es;
  46. struct GAUSS *en;
  47. if ((en = (struct GAUSS *)malloc(sizeof(struct GAUSS))) == NULL)
  48. return (NULL);
  49. es = e * e;
  50. EN->e = e;
  51. sphi = sin(phi0);
  52. cphi = cos(phi0); cphi *= cphi;
  53. *rc = sqrt(1. - es) / (1. - es * sphi * sphi);
  54. EN->C = sqrt(1. + es * cphi * cphi / (1. - es));
  55. *chi = asin(sphi / EN->C);
  56. EN->ratexp = 0.5 * EN->C * e;
  57. EN->K = tan(.5 * *chi + FORTPI) / (
  58. pow(tan(.5 * phi0 + FORTPI), EN->C) *
  59. srat(EN->e * sphi, EN->ratexp) );
  60. return ((void *)en);
  61. }
  62. LP
  63. pj_gauss(LP elp, const void *en) {
  64. LP slp;
  65. slp.phi = 2. * atan( EN->K *
  66. pow(tan(.5 * elp.phi + FORTPI), EN->C) *
  67. srat(EN->e * sin(elp.phi), EN->ratexp) ) - HALFPI;
  68. slp.lam = EN->C * (elp.lam);
  69. return(slp);
  70. }
  71. LP
  72. pj_inv_gauss(LP slp, const void *en) {
  73. LP elp;
  74. double num;
  75. int i;
  76. elp.lam = slp.lam / EN->C;
  77. num = pow(tan(.5 * slp.phi + FORTPI)/EN->K, 1./EN->C);
  78. for (i = MAX_ITER; i; --i) {
  79. elp.phi = 2. * atan(num * srat(EN->e * sin(slp.phi), -.5 * EN->e))
  80. - HALFPI;
  81. if (fabs(elp.phi - slp.phi) < DEL_TOL) break;
  82. slp.phi = elp.phi;
  83. }
  84. /* convergence failed */
  85. if (!i)
  86. pj_errno = -17;
  87. return (elp);
  88. }
  89. /* Revision Log:
  90. ** $Log: pj_gauss.c,v $
  91. ** Revision 1.1 2004/10/20 17:04:00 fwarmerdam
  92. ** New
  93. **
  94. ** Revision 2.2 2004/03/15 16:07:42 gie
  95. ** removed es from init structure
  96. **
  97. ** Revision 2.1 2003/03/28 01:44:30 gie
  98. ** Initial
  99. **
  100. */