/Proj4/PJ_robin.c

  1#ifndef lint
  2static const char SCCSID[]="@(#)PJ_robin.c	4.1 94/02/15     GIE     REL";
  3#endif
  4#define PJ_LIB__
  5#include	"projects.h"
  6PROJ_HEAD(robin, "Robinson") "\n\tPCyl., Sph.";
  7#define V(C,z) (C.c0 + z * (C.c1 + z * (C.c2 + z * C.c3)))
  8#define DV(C,z) (C.c1 + z * (C.c2 + C.c2 + z * 3. * C.c3))
  9/* note: following terms based upon 5 deg. intervals in degrees. */
 10static struct COEFS {
 11	float c0, c1, c2, c3;
 12} X[] = {
 131,	-5.67239e-12,	-7.15511e-05,	3.11028e-06,
 140.9986,	-0.000482241,	-2.4897e-05,	-1.33094e-06,
 150.9954,	-0.000831031,	-4.4861e-05,	-9.86588e-07,
 160.99,	-0.00135363,	-5.96598e-05,	3.67749e-06,
 170.9822,	-0.00167442,	-4.4975e-06,	-5.72394e-06,
 180.973,	-0.00214869,	-9.03565e-05,	1.88767e-08,
 190.96,	-0.00305084,	-9.00732e-05,	1.64869e-06,
 200.9427,	-0.00382792,	-6.53428e-05,	-2.61493e-06,
 210.9216,	-0.00467747,	-0.000104566,	4.8122e-06,
 220.8962,	-0.00536222,	-3.23834e-05,	-5.43445e-06,
 230.8679,	-0.00609364,	-0.0001139,	3.32521e-06,
 240.835,	-0.00698325,	-6.40219e-05,	9.34582e-07,
 250.7986,	-0.00755337,	-5.00038e-05,	9.35532e-07,
 260.7597,	-0.00798325,	-3.59716e-05,	-2.27604e-06,
 270.7186,	-0.00851366,	-7.0112e-05,	-8.63072e-06,
 280.6732,	-0.00986209,	-0.000199572,	1.91978e-05,
 290.6213,	-0.010418,	8.83948e-05,	6.24031e-06,
 300.5722,	-0.00906601,	0.000181999,	6.24033e-06,
 310.5322, 0.,0.,0.  },
 32Y[] = {
 330,	0.0124,	3.72529e-10,	1.15484e-09,
 340.062,	0.0124001,	1.76951e-08,	-5.92321e-09,
 350.124,	0.0123998,	-7.09668e-08,	2.25753e-08,
 360.186,	0.0124008,	2.66917e-07,	-8.44523e-08,
 370.248,	0.0123971,	-9.99682e-07,	3.15569e-07,
 380.31,	0.0124108,	3.73349e-06,	-1.1779e-06,
 390.372,	0.0123598,	-1.3935e-05,	4.39588e-06,
 400.434,	0.0125501,	5.20034e-05,	-1.00051e-05,
 410.4968,	0.0123198,	-9.80735e-05,	9.22397e-06,
 420.5571,	0.0120308,	4.02857e-05,	-5.2901e-06,
 430.6176,	0.0120369,	-3.90662e-05,	7.36117e-07,
 440.6769,	0.0117015,	-2.80246e-05,	-8.54283e-07,
 450.7346,	0.0113572,	-4.08389e-05,	-5.18524e-07,
 460.7903,	0.0109099,	-4.86169e-05,	-1.0718e-06,
 470.8435,	0.0103433,	-6.46934e-05,	5.36384e-09,
 480.8936,	0.00969679,	-6.46129e-05,	-8.54894e-06,
 490.9394,	0.00840949,	-0.000192847,	-4.21023e-06,
 500.9761,	0.00616525,	-0.000256001,	-4.21021e-06,
 511., 0.,0.,0 };
 52#define FXC	0.8487
 53#define FYC	1.3523
 54#define C1	11.45915590261646417544
 55#define RC1	0.08726646259971647884
 56#define NODES	18
 57#define ONEEPS	1.000001
 58#define EPS	1e-8
 59FORWARD(s_forward); /* spheroid */
 60	int i;
 61	double dphi;
 62
 63	i = floor((dphi = fabs(lp.phi)) * C1);
 64	if (i >= NODES) i = NODES - 1;
 65	dphi = RAD_TO_DEG * (dphi - RC1 * i);
 66	xy.x = V(X[i], dphi) * FXC * lp.lam;
 67	xy.y = V(Y[i], dphi) * FYC;
 68	if (lp.phi < 0.) xy.y = -xy.y;
 69	return (xy);
 70}
 71INVERSE(s_inverse); /* spheroid */
 72	int i;
 73	double t, t1;
 74	struct COEFS T;
 75
 76	lp.lam = xy.x / FXC;
 77	lp.phi = fabs(xy.y / FYC);
 78	if (lp.phi >= 1.) { /* simple pathologic cases */
 79		if (lp.phi > ONEEPS) I_ERROR
 80		else {
 81			lp.phi = xy.y < 0. ? -HALFPI : HALFPI;
 82			lp.lam /= X[NODES].c0;
 83		}
 84	} else { /* general problem */
 85		/* in Y space, reduce to table interval */
 86		for (i = floor(lp.phi * NODES);;) {
 87			if (Y[i].c0 > lp.phi) --i;
 88			else if (Y[i+1].c0 <= lp.phi) ++i;
 89			else break;
 90		}
 91		T = Y[i];
 92		/* first guess, linear interp */
 93		t = 5. * (lp.phi - T.c0)/(Y[i+1].c0 - T.c0);
 94		/* make into root */
 95		T.c0 -= lp.phi;
 96		for (;;) { /* Newton-Raphson reduction */
 97			t -= t1 = V(T,t) / DV(T,t);
 98			if (fabs(t1) < EPS)
 99				break;
100		}
101		lp.phi = (5 * i + t) * DEG_TO_RAD;
102		if (xy.y < 0.) lp.phi = -lp.phi;
103		lp.lam /= V(X[i], t);
104	}
105	return (lp);
106}
107FREEUP; if (P) pj_dalloc(P); }
108ENTRY0(robin) P->es = 0.; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)