/Proj4/PJ_nsper.c

  1#ifndef lint
  2static const char SCCSID[]="@(#)PJ_nsper.c	4.1	94/02/15	GIE	REL";
  3#endif
  4#define PROJ_PARMS__ \
  5	double	height; \
  6	double	sinph0; \
  7	double	cosph0; \
  8	double	p; \
  9	double	rp; \
 10	double	pn1; \
 11	double	pfact; \
 12	double	h; \
 13	double	cg; \
 14	double	sg; \
 15	double	sw; \
 16	double	cw; \
 17	int		mode; \
 18	int		tilt;
 19#define PJ_LIB__
 20#include	"projects.h"
 21PROJ_HEAD(nsper, "Near-sided perspective") "\n\tAzi, Sph\n\th=";
 22PROJ_HEAD(tpers, "Tilted perspective") "\n\tAzi, Sph\n\ttilt= azi= h=";
 23# define EPS10 1.e-10
 24# define N_POLE	0
 25# define S_POLE 1
 26# define EQUIT	2
 27# define OBLIQ	3
 28FORWARD(s_forward); /* spheroid */
 29	double  coslam, cosphi, sinphi;
 30
 31	sinphi = sin(lp.phi);
 32	cosphi = cos(lp.phi);
 33	coslam = cos(lp.lam);
 34	switch (P->mode) {
 35	case OBLIQ:
 36		xy.y = P->sinph0 * sinphi + P->cosph0 * cosphi * coslam;
 37		break;
 38	case EQUIT:
 39		xy.y = cosphi * coslam;
 40		break;
 41	case S_POLE:
 42		xy.y = - sinphi;
 43		break;
 44	case N_POLE:
 45		xy.y = sinphi;
 46		break;
 47	}
 48	if (xy.y < P->rp) F_ERROR;
 49	xy.y = P->pn1 / (P->p - xy.y);
 50	xy.x = xy.y * cosphi * sin(lp.lam);
 51	switch (P->mode) {
 52	case OBLIQ:
 53		xy.y *= (P->cosph0 * sinphi -
 54		   P->sinph0 * cosphi * coslam);
 55		break;
 56	case EQUIT:
 57		xy.y *= sinphi;
 58		break;
 59	case N_POLE:
 60		coslam = - coslam;
 61	case S_POLE:
 62		xy.y *= cosphi * coslam;
 63		break;
 64	}
 65	if (P->tilt) {
 66		double yt, ba;
 67
 68		yt = xy.y * P->cg + xy.x * P->sg;
 69		ba = 1. / (yt * P->sw * P->h + P->cw);
 70		xy.x = (xy.x * P->cg - xy.y * P->sg) * P->cw * ba;
 71		xy.y = yt * ba;
 72	}
 73	return (xy);
 74}
 75INVERSE(s_inverse); /* spheroid */
 76	double  rh, cosz, sinz;
 77
 78	if (P->tilt) {
 79		double bm, bq, yt;
 80
 81		yt = 1./(P->pn1 - xy.y * P->sw);
 82		bm = P->pn1 * xy.x * yt;
 83		bq = P->pn1 * xy.y * P->cw * yt;
 84		xy.x = bm * P->cg + bq * P->sg;
 85		xy.y = bq * P->cg - bm * P->sg;
 86	}
 87	rh = hypot(xy.x, xy.y);
 88	if ((sinz = 1. - rh * rh * P->pfact) < 0.) I_ERROR;
 89	sinz = (P->p - sqrt(sinz)) / (P->pn1 / rh + rh / P->pn1);
 90	cosz = sqrt(1. - sinz * sinz);
 91	if (fabs(rh) <= EPS10) {
 92		lp.lam = 0.;
 93		lp.phi = P->phi0;
 94	} else {
 95		switch (P->mode) {
 96		case OBLIQ:
 97			lp.phi = asin(cosz * P->sinph0 + xy.y * sinz * P->cosph0 / rh);
 98			xy.y = (cosz - P->sinph0 * sin(lp.phi)) * rh;
 99			xy.x *= sinz * P->cosph0;
100			break;
101		case EQUIT:
102			lp.phi = asin(xy.y * sinz / rh);
103			xy.y = cosz * rh;
104			xy.x *= sinz;
105			break;
106		case N_POLE:
107			lp.phi = asin(cosz);
108			xy.y = -xy.y;
109			break;
110		case S_POLE:
111			lp.phi = - asin(cosz);
112			break;
113		}
114		lp.lam = atan2(xy.x, xy.y);
115	}
116	return (lp);
117}
118FREEUP; if (P) pj_dalloc(P); }
119	static PJ *
120setup(PJ *P) {
121	if ((P->height = pj_param(P->params, "dh").f) <= 0.) E_ERROR(-30);
122	if (fabs(fabs(P->phi0) - HALFPI) < EPS10)
123		P->mode = P->phi0 < 0. ? S_POLE : N_POLE;
124	else if (fabs(P->phi0) < EPS10)
125		P->mode = EQUIT;
126	else {
127		P->mode = OBLIQ;
128		P->sinph0 = sin(P->phi0);
129		P->cosph0 = cos(P->phi0);
130	}
131	P->pn1 = P->height / P->a; /* normalize by radius */
132	P->p = 1. + P->pn1;
133	P->rp = 1. / P->p;
134	P->h = 1. / P->pn1;
135	P->pfact = (P->p + 1.) * P->h;
136	P->inv = s_inverse;
137	P->fwd = s_forward;
138	P->es = 0.;
139	return P;
140}
141ENTRY0(nsper)
142	P->tilt = 0;
143ENDENTRY(setup(P))
144ENTRY0(tpers)
145	double omega, gamma;
146
147	omega = pj_param(P->params, "dtilt").f * DEG_TO_RAD;
148	gamma = pj_param(P->params, "dazi").f * DEG_TO_RAD;
149	P->tilt = 1;
150	P->cg = cos(gamma); P->sg = sin(gamma);
151	P->cw = cos(omega); P->sw = sin(omega);
152ENDENTRY(setup(P))