PageRenderTime 341ms CodeModel.GetById 326ms RepoModel.GetById 0ms app.codeStats 0ms

/Proj4/PJ_nzmg.c

http://github.com/route-me/route-me
C | 114 lines | 71 code | 6 blank | 37 comment | 6 complexity | da6f497d7eb9244acf1b797c06f512af MD5 | raw file
  1. /******************************************************************************
  2. * $Id: PJ_nzmg.c,v 1.3 2002/12/14 19:37:29 warmerda Exp $
  3. *
  4. * Project: PROJ.4
  5. * Purpose: Implementation of the nzmg (New Zealand Map Grid) projection.
  6. * Very loosely based upon DMA code by Bradford W. Drew
  7. * Author: Gerald Evenden
  8. *
  9. ******************************************************************************
  10. * Copyright (c) 1995, Gerald Evenden
  11. *
  12. * Permission is hereby granted, free of charge, to any person obtaining a
  13. * copy of this software and associated documentation files (the "Software"),
  14. * to deal in the Software without restriction, including without limitation
  15. * the rights to use, copy, modify, merge, publish, distribute, sublicense,
  16. * and/or sell copies of the Software, and to permit persons to whom the
  17. * Software is furnished to do so, subject to the following conditions:
  18. *
  19. * The above copyright notice and this permission notice shall be included
  20. * in all copies or substantial portions of the Software.
  21. *
  22. * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  23. * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  24. * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  25. * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  26. * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  27. * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  28. * DEALINGS IN THE SOFTWARE.
  29. ******************************************************************************
  30. *
  31. * $Log: PJ_nzmg.c,v $
  32. * Revision 1.3 2002/12/14 19:37:29 warmerda
  33. * updated headers
  34. *
  35. */
  36. /* */
  37. #define PJ_LIB__
  38. #include "projects.h"
  39. PJ_CVSID("$Id: PJ_nzmg.c,v 1.3 2002/12/14 19:37:29 warmerda Exp $");
  40. PROJ_HEAD(nzmg, "New Zealand Map Grid") "\n\tfixed Earth";
  41. #define EPSLN 1e-10
  42. #define SEC5_TO_RAD 0.4848136811095359935899141023
  43. #define RAD_TO_SEC5 2.062648062470963551564733573
  44. static COMPLEX
  45. bf[] = {
  46. {.7557853228, 0.0},
  47. {.249204646, .003371507},
  48. {-.001541739, .041058560},
  49. {-.10162907, .01727609},
  50. {-.26623489, -.36249218},
  51. {-.6870983, -1.1651967} };
  52. static double
  53. tphi[] = { 1.5627014243, .5185406398, -.03333098, -.1052906, -.0368594,
  54. .007317, .01220, .00394, -.0013 },
  55. tpsi[] = { .6399175073, -.1358797613, .063294409, -.02526853, .0117879,
  56. -.0055161, .0026906, -.001333, .00067, -.00034 };
  57. #define Nbf 5
  58. #define Ntpsi 9
  59. #define Ntphi 8
  60. FORWARD(e_forward); /* ellipsoid */
  61. COMPLEX p;
  62. double *C;
  63. int i;
  64. lp.phi = (lp.phi - P->phi0) * RAD_TO_SEC5;
  65. for (p.r = *(C = tpsi + (i = Ntpsi)); i ; --i)
  66. p.r = *--C + lp.phi * p.r;
  67. p.r *= lp.phi;
  68. p.i = lp.lam;
  69. p = pj_zpoly1(p, bf, Nbf);
  70. xy.x = p.i;
  71. xy.y = p.r;
  72. return xy;
  73. }
  74. INVERSE(e_inverse); /* ellipsoid */
  75. int nn, i;
  76. COMPLEX p, f, fp, dp;
  77. double den, *C;
  78. p.r = xy.y;
  79. p.i = xy.x;
  80. for (nn = 20; nn ;--nn) {
  81. f = pj_zpolyd1(p, bf, Nbf, &fp);
  82. f.r -= xy.y;
  83. f.i -= xy.x;
  84. den = fp.r * fp.r + fp.i * fp.i;
  85. p.r += dp.r = -(f.r * fp.r + f.i * fp.i) / den;
  86. p.i += dp.i = -(f.i * fp.r - f.r * fp.i) / den;
  87. if ((fabs(dp.r) + fabs(dp.i)) <= EPSLN)
  88. break;
  89. }
  90. if (nn) {
  91. lp.lam = p.i;
  92. for (lp.phi = *(C = tphi + (i = Ntphi)); i ; --i)
  93. lp.phi = *--C + p.r * lp.phi;
  94. lp.phi = P->phi0 + p.r * lp.phi * SEC5_TO_RAD;
  95. } else
  96. lp.lam = lp.phi = HUGE_VAL;
  97. return lp;
  98. }
  99. FREEUP; if (P) pj_dalloc(P); }
  100. ENTRY0(nzmg)
  101. /* force to International major axis */
  102. P->ra = 1. / (P->a = 6378388.0);
  103. P->lam0 = DEG_TO_RAD * 173.;
  104. P->phi0 = DEG_TO_RAD * -41.;
  105. P->x0 = 2510000.;
  106. P->y0 = 6023150.;
  107. P->inv = e_inverse; P->fwd = e_forward;
  108. ENDENTRY(P)