/Proj4/PJ_ortho.c

 1#ifndef lint
 2static const char SCCSID[]="@(#)PJ_ortho.c	4.1	94/02/15	GIE	REL";
 3#endif
 4#define PROJ_PARMS__ \
 5	double	sinph0; \
 6	double	cosph0; \
 7	int		mode;
 8#define PJ_LIB__
 9#include	"projects.h"
10PROJ_HEAD(ortho, "Orthographic") "\n\tAzi, Sph.";
11#define EPS10 1.e-10
12#define N_POLE	0
13#define S_POLE 1
14#define EQUIT	2
15#define OBLIQ	3
16FORWARD(s_forward); /* spheroid */
17	double  coslam, cosphi, sinphi;
18
19	cosphi = cos(lp.phi);
20	coslam = cos(lp.lam);
21	switch (P->mode) {
22	case EQUIT:
23		if (cosphi * coslam < - EPS10) F_ERROR;
24		xy.y = sin(lp.phi);
25		break;
26	case OBLIQ:
27		if (P->sinph0 * (sinphi = sin(lp.phi)) +
28		   P->cosph0 * cosphi * coslam < - EPS10) F_ERROR;
29		xy.y = P->cosph0 * sinphi - P->sinph0 * cosphi * coslam;
30		break;
31	case N_POLE:
32		coslam = - coslam;
33	case S_POLE:
34		if (fabs(lp.phi - P->phi0) - EPS10 > HALFPI) F_ERROR;
35		xy.y = cosphi * coslam;
36		break;
37	}
38	xy.x = cosphi * sin(lp.lam);
39	return (xy);
40}
41
42INVERSE(s_inverse); /* spheroid */
43    double  rh, cosc, sinc;
44
45    if ((sinc = (rh = hypot(xy.x, xy.y))) > 1.) {
46        if ((sinc - 1.) > EPS10) I_ERROR;
47        sinc = 1.;
48    }
49    cosc = sqrt(1. - sinc * sinc); /* in this range OK */
50    if (fabs(rh) <= EPS10) {
51        lp.phi = P->phi0;
52        lp.lam = 0.0;
53    } else {
54        switch (P->mode) {
55        case N_POLE:
56            xy.y = -xy.y;
57            lp.phi = acos(sinc);
58            break;
59        case S_POLE:
60            lp.phi = - acos(sinc);
61            break;
62        case EQUIT:
63            lp.phi = xy.y * sinc / rh;
64            xy.x *= sinc;
65            xy.y = cosc * rh;
66            goto sinchk;
67        case OBLIQ:
68            lp.phi = cosc * P->sinph0 + xy.y * sinc * P->cosph0 /rh;
69            xy.y = (cosc - P->sinph0 * lp.phi) * rh;
70            xy.x *= sinc * P->cosph0;
71        sinchk:
72            if (fabs(lp.phi) >= 1.)
73                lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
74            else
75                lp.phi = asin(lp.phi);
76            break;
77        }
78        lp.lam = (xy.y == 0. && (P->mode == OBLIQ || P->mode == EQUIT))
79             ? (xy.x == 0. ? 0. : xy.x < 0. ? -HALFPI : HALFPI)
80                           : atan2(xy.x, xy.y);
81    }
82    return (lp);
83}
84
85FREEUP; if (P) pj_dalloc(P); }
86ENTRY0(ortho)
87	if (fabs(fabs(P->phi0) - HALFPI) <= EPS10)
88		P->mode = P->phi0 < 0. ? S_POLE : N_POLE;
89	else if (fabs(P->phi0) > EPS10) {
90		P->mode = OBLIQ;
91		P->sinph0 = sin(P->phi0);
92		P->cosph0 = cos(P->phi0);
93	} else
94		P->mode = EQUIT;
95	P->inv = s_inverse;
96	P->fwd = s_forward;
97	P->es = 0.;
98ENDENTRY(P)