PageRenderTime 14ms CodeModel.GetById 2ms app.highlight 9ms RepoModel.GetById 2ms app.codeStats 0ms

/Proj4/proj_mdist.c

http://github.com/route-me/route-me
C | 132 lines | 87 code | 4 blank | 41 comment | 9 complexity | 57a36df679a7f7e08818a60f5216408c MD5 | raw file
  1/*
  2** libproj -- library of cartographic projections
  3**
  4** Copyright (c) 2003, 2006   Gerald I. Evenden
  5*/
  6static const char
  7LIBPROJ_ID[] = "$Id: proj_mdist.c,v 1.1 2006/10/18 05:21:31 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/* Computes distance from equator along the meridian to latitude phi
 29** and inverse on unit ellipsoid.
 30** Precision commensurate with double precision.
 31*/
 32#define PROJ_LIB__
 33#include "projects.h"
 34#define MAX_ITER 20
 35#define TOL 1e-14
 36
 37struct MDIST {
 38	int nb;
 39	double es;
 40	double E;
 41	double b[1];
 42};
 43#define B ((struct MDIST *)b)
 44	void *
 45proj_mdist_ini(double es) {
 46	double numf, numfi, twon1, denf, denfi, ens, T, twon;
 47	double den, El, Es;
 48	double E[MAX_ITER];
 49	struct MDIST *b;
 50	int i, j;
 51
 52/* generate E(e^2) and its terms E[] */
 53	ens = es;
 54	numf = twon1 = denfi = 1.;
 55	denf = 1.;
 56	twon = 4.;
 57	Es = El = E[0] = 1.;
 58	for (i = 1; i < MAX_ITER ; ++i) {
 59		numf *= (twon1 * twon1);
 60		den = twon * denf * denf * twon1;
 61		T = numf/den;
 62		Es -= (E[i] = T * ens);
 63		ens *= es;
 64		twon *= 4.;
 65		denf *= ++denfi;
 66		twon1 += 2.;
 67		if (Es == El) /* jump out if no change */
 68			break;
 69		El = Es;
 70	}
 71	if ((b = (struct MDIST *)malloc(sizeof(struct MDIST)+
 72		(i*sizeof(double)))) == NULL)
 73		return(NULL);
 74	b->nb = i - 1;
 75	b->es = es;
 76	b->E = Es;
 77	/* generate b_n coefficients--note: collapse with prefix ratios */
 78	b->b[0] = Es = 1. - Es;
 79	numf = denf = 1.;
 80	numfi = 2.;
 81	denfi = 3.;
 82	for (j = 1; j < i; ++j) {
 83		Es -= E[j];
 84		numf *= numfi;
 85		denf *= denfi;
 86		b->b[j] = Es * numf / denf;
 87		numfi += 2.;
 88		denfi += 2.;
 89	}
 90	return (b);
 91}
 92	double
 93proj_mdist(double phi, double sphi, double cphi, const void *b) {
 94	double sc, sum, sphi2, D;
 95	int i;
 96
 97	sc = sphi * cphi;
 98	sphi2 = sphi * sphi;
 99	D = phi * B->E - B->es * sc / sqrt(1. - B->es * sphi2);
100	sum = B->b[i = B->nb];
101	while (i) sum = B->b[--i] + sphi2 * sum;
102	return(D + sc * sum);
103}
104	double
105proj_inv_mdist(double dist, const void *b) {
106	double s, t, phi, k;
107	int i;
108
109	k = 1./(1.- B->es);
110	i = MAX_ITER;
111	phi = dist;
112	while ( i-- ) {
113		s = sin(phi);
114		t = 1. - B->es * s * s;
115		phi -= t = (proj_mdist(phi, s, cos(phi), b) - dist) *
116			(t * sqrt(t)) * k;
117		if (fabs(t) < TOL) /* that is no change */
118			return phi;
119	}
120		/* convergence failed */
121	pj_errno = -17;
122	return phi;
123}
124/* Revision Log:
125** $Log: proj_mdist.c,v $
126** Revision 1.1  2006/10/18 05:21:31  fwarmerdam
127** New
128**
129** Revision 3.1  2006/01/11 01:38:18  gie
130** Initial
131**
132*/