diff options
| author | Even Rouault <even.rouault@spatialys.com> | 2018-12-18 20:24:11 +0100 |
|---|---|---|
| committer | Even Rouault <even.rouault@spatialys.com> | 2018-12-26 10:08:53 +0100 |
| commit | 610957f7035242f15743c399ffd429b92bc36206 (patch) | |
| tree | 73f0d51147e2f4860c4bfc875f7a4bf9359386d4 /src/proj_mdist.cpp | |
| parent | 355d681ed88019e97742344bd642c2fd97e700a1 (diff) | |
| download | PROJ-610957f7035242f15743c399ffd429b92bc36206.tar.gz PROJ-610957f7035242f15743c399ffd429b92bc36206.zip | |
cpp conversion: minimal steps to fix compilation errors, not warnings
Diffstat (limited to 'src/proj_mdist.cpp')
| -rw-r--r-- | src/proj_mdist.cpp | 127 |
1 files changed, 127 insertions, 0 deletions
diff --git a/src/proj_mdist.cpp b/src/proj_mdist.cpp new file mode 100644 index 00000000..777f704d --- /dev/null +++ b/src/proj_mdist.cpp @@ -0,0 +1,127 @@ +/* +** libproj -- library of cartographic projections +** +** Copyright (c) 2003, 2006 Gerald I. Evenden +*/ +/* +** Permission is hereby granted, free of charge, to any person obtaining +** a copy of this software and associated documentation files (the +** "Software"), to deal in the Software without restriction, including +** without limitation the rights to use, copy, modify, merge, publish, +** distribute, sublicense, and/or sell copies of the Software, and to +** permit persons to whom the Software is furnished to do so, subject to +** the following conditions: +** +** The above copyright notice and this permission notice shall be +** included in all copies or substantial portions of the Software. +** +** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY +** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, +** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE +** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. +*/ +/* Computes distance from equator along the meridian to latitude phi +** and inverse on unit ellipsoid. +** Precision commensurate with double precision. +*/ +#define PJ_LIB__ + +#include <math.h> +#include <stdlib.h> + +#include "projects.h" + +#define MAX_ITER 20 +#define TOL 1e-14 + +struct MDIST { + int nb; + double es; + double E; + double b[1]; +}; + void * +proj_mdist_ini(double es) { + double numf, numfi, twon1, denf, denfi, ens, T, twon; + double den, El, Es; + double E[MAX_ITER]; + struct MDIST *b; + int i, j; + +/* generate E(e^2) and its terms E[] */ + ens = es; + numf = twon1 = denfi = 1.; + denf = 1.; + twon = 4.; + Es = El = E[0] = 1.; + for (i = 1; i < MAX_ITER ; ++i) { + numf *= (twon1 * twon1); + den = twon * denf * denf * twon1; + T = numf/den; + Es -= (E[i] = T * ens); + ens *= es; + twon *= 4.; + denf *= ++denfi; + twon1 += 2.; + if (Es == El) /* jump out if no change */ + break; + El = Es; + } + if ((b = (struct MDIST *)malloc(sizeof(struct MDIST)+ + (i*sizeof(double)))) == NULL) + return(NULL); + b->nb = i - 1; + b->es = es; + b->E = Es; + /* generate b_n coefficients--note: collapse with prefix ratios */ + b->b[0] = Es = 1. - Es; + numf = denf = 1.; + numfi = 2.; + denfi = 3.; + for (j = 1; j < i; ++j) { + Es -= E[j]; + numf *= numfi; + denf *= denfi; + b->b[j] = Es * numf / denf; + numfi += 2.; + denfi += 2.; + } + return (b); +} + double +proj_mdist(double phi, double sphi, double cphi, const void *data) { + const struct MDIST *b = (const struct MDIST *)data; + double sc, sum, sphi2, D; + int i; + + sc = sphi * cphi; + sphi2 = sphi * sphi; + D = phi * b->E - b->es * sc / sqrt(1. - b->es * sphi2); + sum = b->b[i = b->nb]; + while (i) sum = b->b[--i] + sphi2 * sum; + return(D + sc * sum); +} + double +proj_inv_mdist(projCtx ctx, double dist, const void *data) { + const struct MDIST *b = (const struct MDIST *)data; + double s, t, phi, k; + int i; + + k = 1./(1.- b->es); + i = MAX_ITER; + phi = dist; + while ( i-- ) { + s = sin(phi); + t = 1. - b->es * s * s; + phi -= t = (proj_mdist(phi, s, cos(phi), b) - dist) * + (t * sqrt(t)) * k; + if (fabs(t) < TOL) /* that is no change */ + return phi; + } + /* convergence failed */ + pj_ctx_set_errno(ctx, PJD_ERR_NON_CONV_INV_MERI_DIST); + return phi; +} |
