/Proj4/PJ_nzmg.c

  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
 41PJ_CVSID("$Id: PJ_nzmg.c,v 1.3 2002/12/14 19:37:29 warmerda Exp $");
 42
 43PROJ_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
 49bf[] = {
 50	{.7557853228,	0.0},
 51	{.249204646,	.003371507},
 52	{-.001541739,	.041058560},
 53	{-.10162907,	.01727609},
 54	{-.26623489,	-.36249218},
 55	{-.6870983,	-1.1651967} };
 56	static double
 57tphi[] = { 1.5627014243, .5185406398, -.03333098, -.1052906, -.0368594,
 58	.007317, .01220, .00394, -.0013 },
 59tpsi[] = { .6399175073, -.1358797613, .063294409, -.02526853, .0117879,
 60	-.0055161, .0026906, -.001333, .00067, -.00034 };
 61#define Nbf 5
 62#define Ntpsi 9
 63#define Ntphi 8
 64FORWARD(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}
 79INVERSE(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}
105FREEUP; if (P) pj_dalloc(P); }
106ENTRY0(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;
114ENDENTRY(P)