diff options
Diffstat (limited to 'src/3rdparty/proj/PJ_mbtfpq.c')
| -rw-r--r-- | src/3rdparty/proj/PJ_mbtfpq.c | 47 |
1 files changed, 47 insertions, 0 deletions
diff --git a/src/3rdparty/proj/PJ_mbtfpq.c b/src/3rdparty/proj/PJ_mbtfpq.c new file mode 100644 index 00000000..3298db4f --- /dev/null +++ b/src/3rdparty/proj/PJ_mbtfpq.c @@ -0,0 +1,47 @@ +#define PJ_LIB__ +#include <projects.h> +PROJ_HEAD(mbtfpq, "McBryde-Thomas Flat-Polar Quartic") "\n\tCyl., Sph."; +#define NITER 20 +#define EPS 1e-7 +#define ONETOL 1.000001 +#define C 1.70710678118654752440 +#define RC 0.58578643762690495119 +#define FYC 1.87475828462269495505 +#define RYC 0.53340209679417701685 +#define FXC 0.31245971410378249250 +#define RXC 3.20041258076506210122 +FORWARD(s_forward); /* spheroid */ + double th1, c; + int i; + + c = C * sin(lp.phi); + for (i = NITER; i; --i) { + lp.phi -= th1 = (sin(.5*lp.phi) + sin(lp.phi) - c) / + (.5*cos(.5*lp.phi) + cos(lp.phi)); + if (fabs(th1) < EPS) break; + } + xy.x = FXC * lp.lam * (1.0 + 2. * cos(lp.phi)/cos(0.5 * lp.phi)); + xy.y = FYC * sin(0.5 * lp.phi); + return (xy); +} +INVERSE(s_inverse); /* spheroid */ + double t; + + lp.phi = RYC * xy.y; + if (fabs(lp.phi) > 1.) { + if (fabs(lp.phi) > ONETOL) I_ERROR + else if (lp.phi < 0.) { t = -1.; lp.phi = -PI; } + else { t = 1.; lp.phi = PI; } + } else + lp.phi = 2. * asin(t = lp.phi); + lp.lam = RXC * xy.x / (1. + 2. * cos(lp.phi)/cos(0.5 * lp.phi)); + lp.phi = RC * (t + sin(lp.phi)); + if (fabs(lp.phi) > 1.) + if (fabs(lp.phi) > ONETOL) I_ERROR + else lp.phi = lp.phi < 0. ? -HALFPI : HALFPI; + else + lp.phi = asin(lp.phi); + return (lp); +} +FREEUP; if (P) pj_dalloc(P); } +ENTRY0(mbtfpq) P->es = 0.; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P) |