root / org.gvsig.jcrs / libjni-proj4 / src / PJ_nzmg.c @ 36
History | View | Annotate | Download (3.56 KB)
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 |
/* */
|
38 |
#define PJ_LIB__
|
39 |
#include <projects.h> |
40 |
|
41 |
PJ_CVSID("$Id: PJ_nzmg.c,v 1.3 2002/12/14 19:37:29 warmerda Exp $");
|
42 |
|
43 |
PROJ_HEAD(nzmg, "New Zealand Map Grid") "\n\tfixed Earth"; |
44 |
|
45 |
#define EPSLN 1e-10 |
46 |
#define SEC5_TO_RAD 0.4848136811095359935899141023 |
47 |
#define RAD_TO_SEC5 2.062648062470963551564733573 |
48 |
static COMPLEX
|
49 |
bf[] = { |
50 |
{.7557853228, 0.0}, |
51 |
{.249204646, .003371507}, |
52 |
{-.001541739, .041058560}, |
53 |
{-.10162907, .01727609}, |
54 |
{-.26623489, -.36249218}, |
55 |
{-.6870983, -1.1651967} }; |
56 |
static double |
57 |
tphi[] = { 1.5627014243, .5185406398, -.03333098, -.1052906, -.0368594, |
58 |
.007317, .01220, .00394, -.0013 }, |
59 |
tpsi[] = { .6399175073, -.1358797613, .063294409, -.02526853, .0117879, |
60 |
-.0055161, .0026906, -.001333, .00067, -.00034 }; |
61 |
#define Nbf 5 |
62 |
#define Ntpsi 9 |
63 |
#define Ntphi 8 |
64 |
FORWARD(e_forward); /* ellipsoid */
|
65 |
COMPLEX p; |
66 |
double *C;
|
67 |
int i;
|
68 |
|
69 |
lp.phi = (lp.phi - P->phi0) * RAD_TO_SEC5; |
70 |
for (p.r = *(C = tpsi + (i = Ntpsi)); i ; --i)
|
71 |
p.r = *--C + lp.phi * p.r; |
72 |
p.r *= lp.phi; |
73 |
p.i = lp.lam; |
74 |
p = pj_zpoly1(p, bf, Nbf); |
75 |
xy.x = p.i; |
76 |
xy.y = p.r; |
77 |
return xy;
|
78 |
} |
79 |
INVERSE(e_inverse); /* ellipsoid */
|
80 |
int nn, i;
|
81 |
COMPLEX p, f, fp, dp; |
82 |
double den, *C;
|
83 |
|
84 |
p.r = xy.y; |
85 |
p.i = xy.x; |
86 |
for (nn = 20; nn ;--nn) { |
87 |
f = pj_zpolyd1(p, bf, Nbf, &fp); |
88 |
f.r -= xy.y; |
89 |
f.i -= xy.x; |
90 |
den = fp.r * fp.r + fp.i * fp.i; |
91 |
p.r += dp.r = -(f.r * fp.r + f.i * fp.i) / den; |
92 |
p.i += dp.i = -(f.i * fp.r - f.r * fp.i) / den; |
93 |
if ((fabs(dp.r) + fabs(dp.i)) <= EPSLN)
|
94 |
break;
|
95 |
} |
96 |
if (nn) {
|
97 |
lp.lam = p.i; |
98 |
for (lp.phi = *(C = tphi + (i = Ntphi)); i ; --i)
|
99 |
lp.phi = *--C + p.r * lp.phi; |
100 |
lp.phi = P->phi0 + p.r * lp.phi * SEC5_TO_RAD; |
101 |
} else
|
102 |
lp.lam = lp.phi = HUGE_VAL; |
103 |
return lp;
|
104 |
} |
105 |
FREEUP; if (P) pj_dalloc(P); }
|
106 |
ENTRY0(nzmg) |
107 |
/* force to International major axis */
|
108 |
P->ra = 1. / (P->a = 6378388.0); |
109 |
P->lam0 = DEG_TO_RAD * 173.; |
110 |
P->phi0 = DEG_TO_RAD * -41.; |
111 |
P->x0 = 2510000.; |
112 |
P->y0 = 6023150.; |
113 |
P->inv = e_inverse; P->fwd = e_forward; |
114 |
ENDENTRY(P) |