svn-gvsig-desktop / tags / v1_9_Build_1241 / libraries / libjni-proj4 / src / PJ_gnom.c @ 34066
History | View | Annotate | Download (2.19 KB)
| 1 |
#ifndef lint
|
|---|---|
| 2 |
static const char SCCSID[]="@(#)PJ_gnom.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> |
| 10 |
PROJ_HEAD(gnom, "Gnomonic") "\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 |
| 16 |
FORWARD(s_forward); /* spheroid */
|
| 17 |
double coslam, cosphi, sinphi;
|
| 18 |
|
| 19 |
sinphi = sin(lp.phi); |
| 20 |
cosphi = cos(lp.phi); |
| 21 |
coslam = cos(lp.lam); |
| 22 |
switch (P->mode) {
|
| 23 |
case EQUIT:
|
| 24 |
xy.y = cosphi * coslam; |
| 25 |
break;
|
| 26 |
case OBLIQ:
|
| 27 |
xy.y = P->sinph0 * sinphi + P->cosph0 * cosphi * coslam; |
| 28 |
break;
|
| 29 |
case S_POLE:
|
| 30 |
xy.y = - sinphi; |
| 31 |
break;
|
| 32 |
case N_POLE:
|
| 33 |
xy.y = sinphi; |
| 34 |
break;
|
| 35 |
} |
| 36 |
if (xy.y <= EPS10) F_ERROR;
|
| 37 |
xy.x = (xy.y = 1. / xy.y) * cosphi * sin(lp.lam);
|
| 38 |
switch (P->mode) {
|
| 39 |
case EQUIT:
|
| 40 |
xy.y *= sinphi; |
| 41 |
break;
|
| 42 |
case OBLIQ:
|
| 43 |
xy.y *= P->cosph0 * sinphi - P->sinph0 * cosphi * coslam; |
| 44 |
break;
|
| 45 |
case N_POLE:
|
| 46 |
coslam = - coslam; |
| 47 |
case S_POLE:
|
| 48 |
xy.y *= cosphi * coslam; |
| 49 |
break;
|
| 50 |
} |
| 51 |
return (xy);
|
| 52 |
} |
| 53 |
INVERSE(s_inverse); /* spheroid */
|
| 54 |
double rh, cosz, sinz;
|
| 55 |
|
| 56 |
rh = hypot(xy.x, xy.y); |
| 57 |
sinz = sin(lp.phi = atan(rh)); |
| 58 |
cosz = sqrt(1. - sinz * sinz);
|
| 59 |
if (fabs(rh) <= EPS10) {
|
| 60 |
lp.phi = P->phi0; |
| 61 |
lp.lam = 0.;
|
| 62 |
} else {
|
| 63 |
switch (P->mode) {
|
| 64 |
case OBLIQ:
|
| 65 |
lp.phi = cosz * P->sinph0 + xy.y * sinz * P->cosph0 / rh; |
| 66 |
if (fabs(lp.phi) >= 1.) |
| 67 |
lp.phi = lp.phi > 0. ? HALFPI : - HALFPI;
|
| 68 |
else
|
| 69 |
lp.phi = asin(lp.phi); |
| 70 |
xy.y = (cosz - P->sinph0 * sin(lp.phi)) * rh; |
| 71 |
xy.x *= sinz * P->cosph0; |
| 72 |
break;
|
| 73 |
case EQUIT:
|
| 74 |
lp.phi = xy.y * sinz / rh; |
| 75 |
if (fabs(lp.phi) >= 1.) |
| 76 |
lp.phi = lp.phi > 0. ? HALFPI : - HALFPI;
|
| 77 |
else
|
| 78 |
lp.phi = asin(lp.phi); |
| 79 |
xy.y = cosz * rh; |
| 80 |
xy.x *= sinz; |
| 81 |
break;
|
| 82 |
case S_POLE:
|
| 83 |
lp.phi -= HALFPI; |
| 84 |
break;
|
| 85 |
case N_POLE:
|
| 86 |
lp.phi = HALFPI - lp.phi; |
| 87 |
xy.y = -xy.y; |
| 88 |
break;
|
| 89 |
} |
| 90 |
lp.lam = atan2(xy.x, xy.y); |
| 91 |
} |
| 92 |
return (lp);
|
| 93 |
} |
| 94 |
FREEUP; if (P) pj_dalloc(P); }
|
| 95 |
ENTRY0(gnom) |
| 96 |
if (fabs(fabs(P->phi0) - HALFPI) < EPS10)
|
| 97 |
P->mode = P->phi0 < 0. ? S_POLE : N_POLE;
|
| 98 |
else if (fabs(P->phi0) < EPS10) |
| 99 |
P->mode = EQUIT; |
| 100 |
else {
|
| 101 |
P->mode = OBLIQ; |
| 102 |
P->sinph0 = sin(P->phi0); |
| 103 |
P->cosph0 = cos(P->phi0); |
| 104 |
} |
| 105 |
P->inv = s_inverse; |
| 106 |
P->fwd = s_forward; |
| 107 |
P->es = 0.;
|
| 108 |
ENDENTRY(P) |