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