/Proj4/PJ_aeqd.c
C | 277 lines | 232 code | 11 blank | 34 comment | 49 complexity | 734bf86c1287f4407c696e528eacc50b MD5 | raw file
1/****************************************************************************** 2 * $Id: PJ_aeqd.c,v 1.3 2002/12/14 19:27:06 warmerda Exp $ 3 * 4 * Project: PROJ.4 5 * Purpose: Implementation of the aeqd (Azimuthal Equidistant) projection. 6 * Author: Gerald Evenden 7 * 8 ****************************************************************************** 9 * Copyright (c) 1995, Gerald Evenden 10 * 11 * Permission is hereby granted, free of charge, to any person obtaining a 12 * copy of this software and associated documentation files (the "Software"), 13 * to deal in the Software without restriction, including without limitation 14 * the rights to use, copy, modify, merge, publish, distribute, sublicense, 15 * and/or sell copies of the Software, and to permit persons to whom the 16 * Software is furnished to do so, subject to the following conditions: 17 * 18 * The above copyright notice and this permission notice shall be included 19 * in all copies or substantial portions of the Software. 20 * 21 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 22 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 23 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 24 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 25 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 26 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 27 * DEALINGS IN THE SOFTWARE. 28 ****************************************************************************** 29 * 30 * $Log: PJ_aeqd.c,v $ 31 * Revision 1.3 2002/12/14 19:27:06 warmerda 32 * updated header 33 * 34 */ 35 36#define PROJ_PARMS__ \ 37 double sinph0; \ 38 double cosph0; \ 39 double *en; \ 40 double M1; \ 41 double N1; \ 42 double Mp; \ 43 double He; \ 44 double G; \ 45 int mode; 46#define PJ_LIB__ 47#include "projects.h" 48 49PJ_CVSID("$Id: PJ_aeqd.c,v 1.3 2002/12/14 19:27:06 warmerda Exp $"); 50 51PROJ_HEAD(aeqd, "Azimuthal Equidistant") "\n\tAzi, Sph&Ell\n\tlat_0 guam"; 52 53#define EPS10 1.e-10 54#define TOL 1.e-14 55 56#define N_POLE 0 57#define S_POLE 1 58#define EQUIT 2 59#define OBLIQ 3 60FORWARD(e_guam_fwd); /* Guam elliptical */ 61 double cosphi, sinphi, t; 62 63 cosphi = cos(lp.phi); 64 sinphi = sin(lp.phi); 65 t = 1. / sqrt(1. - P->es * sinphi * sinphi); 66 xy.x = lp.lam * cosphi * t; 67 xy.y = pj_mlfn(lp.phi, sinphi, cosphi, P->en) - P->M1 + 68 .5 * lp.lam * lp.lam * cosphi * sinphi * t; 69 return (xy); 70} 71FORWARD(e_forward); /* elliptical */ 72 double coslam, cosphi, sinphi, rho, s, H, H2, c, Az, t, ct, st, cA, sA; 73 74 coslam = cos(lp.lam); 75 cosphi = cos(lp.phi); 76 sinphi = sin(lp.phi); 77 switch (P->mode) { 78 case N_POLE: 79 coslam = - coslam; 80 case S_POLE: 81 xy.x = (rho = fabs(P->Mp - pj_mlfn(lp.phi, sinphi, cosphi, P->en))) * 82 sin(lp.lam); 83 xy.y = rho * coslam; 84 break; 85 case EQUIT: 86 case OBLIQ: 87 if (fabs(lp.lam) < EPS10 && fabs(lp.phi - P->phi0) < EPS10) { 88 xy.x = xy.y = 0.; 89 break; 90 } 91 t = atan2(P->one_es * sinphi + P->es * P->N1 * P->sinph0 * 92 sqrt(1. - P->es * sinphi * sinphi), cosphi); 93 ct = cos(t); st = sin(t); 94 Az = atan2(sin(lp.lam) * ct, P->cosph0 * st - P->sinph0 * coslam * ct); 95 cA = cos(Az); sA = sin(Az); 96 s = aasin( fabs(sA) < TOL ? 97 (P->cosph0 * st - P->sinph0 * coslam * ct) / cA : 98 sin(lp.lam) * ct / sA ); 99 H = P->He * cA; 100 H2 = H * H; 101 c = P->N1 * s * (1. + s * s * (- H2 * (1. - H2)/6. + 102 s * ( P->G * H * (1. - 2. * H2 * H2) / 8. + 103 s * ((H2 * (4. - 7. * H2) - 3. * P->G * P->G * (1. - 7. * H2)) / 104 120. - s * P->G * H / 48.)))); 105 xy.x = c * sA; 106 xy.y = c * cA; 107 break; 108 } 109 return (xy); 110} 111FORWARD(s_forward); /* spherical */ 112 double coslam, cosphi, sinphi; 113 114 sinphi = sin(lp.phi); 115 cosphi = cos(lp.phi); 116 coslam = cos(lp.lam); 117 switch (P->mode) { 118 case EQUIT: 119 xy.y = cosphi * coslam; 120 goto oblcon; 121 case OBLIQ: 122 xy.y = P->sinph0 * sinphi + P->cosph0 * cosphi * coslam; 123oblcon: 124 if (fabs(fabs(xy.y) - 1.) < TOL) 125 if (xy.y < 0.) 126 F_ERROR 127 else 128 xy.x = xy.y = 0.; 129 else { 130 xy.y = acos(xy.y); 131 xy.y /= sin(xy.y); 132 xy.x = xy.y * cosphi * sin(lp.lam); 133 xy.y *= (P->mode == EQUIT) ? sinphi : 134 P->cosph0 * sinphi - P->sinph0 * cosphi * coslam; 135 } 136 break; 137 case N_POLE: 138 lp.phi = -lp.phi; 139 coslam = -coslam; 140 case S_POLE: 141 if (fabs(lp.phi - HALFPI) < EPS10) F_ERROR; 142 xy.x = (xy.y = (HALFPI + lp.phi)) * sin(lp.lam); 143 xy.y *= coslam; 144 break; 145 } 146 return (xy); 147} 148INVERSE(e_guam_inv); /* Guam elliptical */ 149 double x2, t; 150 int i; 151 152 x2 = 0.5 * xy.x * xy.x; 153 lp.phi = P->phi0; 154 for (i = 0; i < 3; ++i) { 155 t = P->e * sin(lp.phi); 156 lp.phi = pj_inv_mlfn(P->M1 + xy.y - 157 x2 * tan(lp.phi) * (t = sqrt(1. - t * t)), P->es, P->en); 158 } 159 lp.lam = xy.x * t / cos(lp.phi); 160 return (lp); 161} 162INVERSE(e_inverse); /* elliptical */ 163 double c, Az, cosAz, A, B, D, E, F, psi, t; 164 165 if ((c = hypot(xy.x, xy.y)) < EPS10) { 166 lp.phi = P->phi0; 167 lp.lam = 0.; 168 return (lp); 169 } 170 if (P->mode == OBLIQ || P->mode == EQUIT) { 171 cosAz = cos(Az = atan2(xy.x, xy.y)); 172 t = P->cosph0 * cosAz; 173 B = P->es * t / P->one_es; 174 A = - B * t; 175 B *= 3. * (1. - A) * P->sinph0; 176 D = c / P->N1; 177 E = D * (1. - D * D * (A * (1. + A) / 6. + B * (1. + 3.*A) * D / 24.)); 178 F = 1. - E * E * (A / 2. + B * E / 6.); 179 psi = aasin(P->sinph0 * cos(E) + t * sin(E)); 180 lp.lam = aasin(sin(Az) * sin(E) / cos(psi)); 181 if ((t = fabs(psi)) < EPS10) 182 lp.phi = 0.; 183 else if (fabs(t - HALFPI) < 0.) 184 lp.phi = HALFPI; 185 else 186 lp.phi = atan((1. - P->es * F * P->sinph0 / sin(psi)) * tan(psi) / 187 P->one_es); 188 } else { /* Polar */ 189 lp.phi = pj_inv_mlfn(P->mode == N_POLE ? P->Mp - c : P->Mp + c, 190 P->es, P->en); 191 lp.lam = atan2(xy.x, P->mode == N_POLE ? -xy.y : xy.y); 192 } 193 return (lp); 194} 195INVERSE(s_inverse); /* spherical */ 196 double cosc, c_rh, sinc; 197 198 if ((c_rh = hypot(xy.x, xy.y)) > PI) { 199 if (c_rh - EPS10 > PI) I_ERROR; 200 c_rh = PI; 201 } else if (c_rh < EPS10) { 202 lp.phi = P->phi0; 203 lp.lam = 0.; 204 return (lp); 205 } 206 if (P->mode == OBLIQ || P->mode == EQUIT) { 207 sinc = sin(c_rh); 208 cosc = cos(c_rh); 209 if (P->mode == EQUIT) { 210 lp.phi = aasin(xy.y * sinc / c_rh); 211 xy.x *= sinc; 212 xy.y = cosc * c_rh; 213 } else { 214 lp.phi = aasin(cosc * P->sinph0 + xy.y * sinc * P->cosph0 / 215 c_rh); 216 xy.y = (cosc - P->sinph0 * sin(lp.phi)) * c_rh; 217 xy.x *= sinc * P->cosph0; 218 } 219 lp.lam = xy.y == 0. ? 0. : atan2(xy.x, xy.y); 220 } else if (P->mode == N_POLE) { 221 lp.phi = HALFPI - c_rh; 222 lp.lam = atan2(xy.x, -xy.y); 223 } else { 224 lp.phi = c_rh - HALFPI; 225 lp.lam = atan2(xy.x, xy.y); 226 } 227 return (lp); 228} 229FREEUP; 230 if (P) { 231 if (P->en) 232 pj_dalloc(P->en); 233 pj_dalloc(P); 234 } 235} 236ENTRY1(aeqd, en) 237 P->phi0 = pj_param(P->params, "rlat_0").f; 238 if (fabs(fabs(P->phi0) - HALFPI) < EPS10) { 239 P->mode = P->phi0 < 0. ? S_POLE : N_POLE; 240 P->sinph0 = P->phi0 < 0. ? -1. : 1.; 241 P->cosph0 = 0.; 242 } else if (fabs(P->phi0) < EPS10) { 243 P->mode = EQUIT; 244 P->sinph0 = 0.; 245 P->cosph0 = 1.; 246 } else { 247 P->mode = OBLIQ; 248 P->sinph0 = sin(P->phi0); 249 P->cosph0 = cos(P->phi0); 250 } 251 if (! P->es) { 252 P->inv = s_inverse; P->fwd = s_forward; 253 } else { 254 if (!(P->en = pj_enfn(P->es))) E_ERROR_0; 255 if (pj_param(P->params, "bguam").i) { 256 P->M1 = pj_mlfn(P->phi0, P->sinph0, P->cosph0, P->en); 257 P->inv = e_guam_inv; P->fwd = e_guam_fwd; 258 } else { 259 switch (P->mode) { 260 case N_POLE: 261 P->Mp = pj_mlfn(HALFPI, 1., 0., P->en); 262 break; 263 case S_POLE: 264 P->Mp = pj_mlfn(-HALFPI, -1., 0., P->en); 265 break; 266 case EQUIT: 267 case OBLIQ: 268 P->inv = e_inverse; P->fwd = e_forward; 269 P->N1 = 1. / sqrt(1. - P->es * P->sinph0 * P->sinph0); 270 P->G = P->sinph0 * (P->He = P->e / sqrt(P->one_es)); 271 P->He *= P->cosph0; 272 break; 273 } 274 P->inv = e_inverse; P->fwd = e_forward; 275 } 276 } 277ENDENTRY(P)