1 /** 2 * \file AlbersEqualArea.cpp 3 * \brief Implementation for GeographicLib::AlbersEqualArea class 4 * 5 * Copyright (c) Charles Karney (2010-2021) <charles@karney.com> and licensed 6 * under the MIT/X11 License. For more information, see 7 * https://geographiclib.sourceforge.io/ 8 **********************************************************************/ 9 10 #include <GeographicLib/AlbersEqualArea.hpp> 11 12 #if defined(_MSC_VER) 13 // Squelch warnings about constant conditional expressions 14 # pragma warning (disable: 4127) 15 #endif 16 17 namespace GeographicLib { 18 19 using namespace std; 20 AlbersEqualArea(real a,real f,real stdlat,real k0)21 AlbersEqualArea::AlbersEqualArea(real a, real f, real stdlat, real k0) 22 : eps_(numeric_limits<real>::epsilon()) 23 , epsx_(Math::sq(eps_)) 24 , epsx2_(Math::sq(epsx_)) 25 , tol_(sqrt(eps_)) 26 , tol0_(tol_ * sqrt(sqrt(eps_))) 27 , _a(a) 28 , _f(f) 29 , _fm(1 - _f) 30 , _e2(_f * (2 - _f)) 31 , _e(sqrt(abs(_e2))) 32 , _e2m(1 - _e2) 33 , _qZ(1 + _e2m * atanhee(real(1))) 34 , _qx(_qZ / ( 2 * _e2m )) 35 { 36 if (!(isfinite(_a) && _a > 0)) 37 throw GeographicErr("Equatorial radius is not positive"); 38 if (!(isfinite(_f) && _f < 1)) 39 throw GeographicErr("Polar semi-axis is not positive"); 40 if (!(isfinite(k0) && k0 > 0)) 41 throw GeographicErr("Scale is not positive"); 42 if (!(abs(stdlat) <= 90)) 43 throw GeographicErr("Standard latitude not in [-90d, 90d]"); 44 real sphi, cphi; 45 Math::sincosd(stdlat, sphi, cphi); 46 Init(sphi, cphi, sphi, cphi, k0); 47 } 48 AlbersEqualArea(real a,real f,real stdlat1,real stdlat2,real k1)49 AlbersEqualArea::AlbersEqualArea(real a, real f, real stdlat1, real stdlat2, 50 real k1) 51 : eps_(numeric_limits<real>::epsilon()) 52 , epsx_(Math::sq(eps_)) 53 , epsx2_(Math::sq(epsx_)) 54 , tol_(sqrt(eps_)) 55 , tol0_(tol_ * sqrt(sqrt(eps_))) 56 , _a(a) 57 , _f(f) 58 , _fm(1 - _f) 59 , _e2(_f * (2 - _f)) 60 , _e(sqrt(abs(_e2))) 61 , _e2m(1 - _e2) 62 , _qZ(1 + _e2m * atanhee(real(1))) 63 , _qx(_qZ / ( 2 * _e2m )) 64 { 65 if (!(isfinite(_a) && _a > 0)) 66 throw GeographicErr("Equatorial radius is not positive"); 67 if (!(isfinite(_f) && _f < 1)) 68 throw GeographicErr("Polar semi-axis is not positive"); 69 if (!(isfinite(k1) && k1 > 0)) 70 throw GeographicErr("Scale is not positive"); 71 if (!(abs(stdlat1) <= 90)) 72 throw GeographicErr("Standard latitude 1 not in [-90d, 90d]"); 73 if (!(abs(stdlat2) <= 90)) 74 throw GeographicErr("Standard latitude 2 not in [-90d, 90d]"); 75 real sphi1, cphi1, sphi2, cphi2; 76 Math::sincosd(stdlat1, sphi1, cphi1); 77 Math::sincosd(stdlat2, sphi2, cphi2); 78 Init(sphi1, cphi1, sphi2, cphi2, k1); 79 } 80 AlbersEqualArea(real a,real f,real sinlat1,real coslat1,real sinlat2,real coslat2,real k1)81 AlbersEqualArea::AlbersEqualArea(real a, real f, 82 real sinlat1, real coslat1, 83 real sinlat2, real coslat2, 84 real k1) 85 : eps_(numeric_limits<real>::epsilon()) 86 , epsx_(Math::sq(eps_)) 87 , epsx2_(Math::sq(epsx_)) 88 , tol_(sqrt(eps_)) 89 , tol0_(tol_ * sqrt(sqrt(eps_))) 90 , _a(a) 91 , _f(f) 92 , _fm(1 - _f) 93 , _e2(_f * (2 - _f)) 94 , _e(sqrt(abs(_e2))) 95 , _e2m(1 - _e2) 96 , _qZ(1 + _e2m * atanhee(real(1))) 97 , _qx(_qZ / ( 2 * _e2m )) 98 { 99 if (!(isfinite(_a) && _a > 0)) 100 throw GeographicErr("Equatorial radius is not positive"); 101 if (!(isfinite(_f) && _f < 1)) 102 throw GeographicErr("Polar semi-axis is not positive"); 103 if (!(isfinite(k1) && k1 > 0)) 104 throw GeographicErr("Scale is not positive"); 105 if (!(coslat1 >= 0)) 106 throw GeographicErr("Standard latitude 1 not in [-90d, 90d]"); 107 if (!(coslat2 >= 0)) 108 throw GeographicErr("Standard latitude 2 not in [-90d, 90d]"); 109 if (!(abs(sinlat1) <= 1 && coslat1 <= 1) || (coslat1 == 0 && sinlat1 == 0)) 110 throw GeographicErr("Bad sine/cosine of standard latitude 1"); 111 if (!(abs(sinlat2) <= 1 && coslat2 <= 1) || (coslat2 == 0 && sinlat2 == 0)) 112 throw GeographicErr("Bad sine/cosine of standard latitude 2"); 113 if (coslat1 == 0 && coslat2 == 0 && sinlat1 * sinlat2 <= 0) 114 throw GeographicErr 115 ("Standard latitudes cannot be opposite poles"); 116 Init(sinlat1, coslat1, sinlat2, coslat2, k1); 117 } 118 Init(real sphi1,real cphi1,real sphi2,real cphi2,real k1)119 void AlbersEqualArea::Init(real sphi1, real cphi1, 120 real sphi2, real cphi2, real k1) { 121 { 122 real r; 123 r = hypot(sphi1, cphi1); 124 sphi1 /= r; cphi1 /= r; 125 r = hypot(sphi2, cphi2); 126 sphi2 /= r; cphi2 /= r; 127 } 128 bool polar = (cphi1 == 0); 129 cphi1 = max(epsx_, cphi1); // Avoid singularities at poles 130 cphi2 = max(epsx_, cphi2); 131 // Determine hemisphere of tangent latitude 132 _sign = sphi1 + sphi2 >= 0 ? 1 : -1; 133 // Internally work with tangent latitude positive 134 sphi1 *= _sign; sphi2 *= _sign; 135 if (sphi1 > sphi2) { 136 swap(sphi1, sphi2); swap(cphi1, cphi2); // Make phi1 < phi2 137 } 138 real 139 tphi1 = sphi1/cphi1, tphi2 = sphi2/cphi2; 140 141 // q = (1-e^2)*(sphi/(1-e^2*sphi^2) - atanhee(sphi)) 142 // qZ = q(pi/2) = (1 + (1-e^2)*atanhee(1)) 143 // atanhee(x) = atanh(e*x)/e 144 // q = sxi * qZ 145 // dq/dphi = 2*(1-e^2)*cphi/(1-e^2*sphi^2)^2 146 // 147 // n = (m1^2-m2^2)/(q2-q1) -> sin(phi0) for phi1, phi2 -> phi0 148 // C = m1^2 + n*q1 = (m1^2*q2-m2^2*q1)/(q2-q1) 149 // let 150 // rho(pi/2)/rho(-pi/2) = (1-s)/(1+s) 151 // s = n*qZ/C 152 // = qZ * (m1^2-m2^2)/(m1^2*q2-m2^2*q1) 153 // = qZ * (scbet2^2 - scbet1^2)/(scbet2^2*q2 - scbet1^2*q1) 154 // = (scbet2^2 - scbet1^2)/(scbet2^2*sxi2 - scbet1^2*sxi1) 155 // = (tbet2^2 - tbet1^2)/(scbet2^2*sxi2 - scbet1^2*sxi1) 156 // 1-s = -((1-sxi2)*scbet2^2 - (1-sxi1)*scbet1^2)/ 157 // (scbet2^2*sxi2 - scbet1^2*sxi1) 158 // 159 // Define phi0 to give same value of s, i.e., 160 // s = sphi0 * qZ / (m0^2 + sphi0*q0) 161 // = sphi0 * scbet0^2 / (1/qZ + sphi0 * scbet0^2 * sxi0) 162 163 real tphi0, C; 164 if (polar || tphi1 == tphi2) { 165 tphi0 = tphi2; 166 C = 1; // ignored 167 } else { 168 real 169 tbet1 = _fm * tphi1, scbet12 = 1 + Math::sq(tbet1), 170 tbet2 = _fm * tphi2, scbet22 = 1 + Math::sq(tbet2), 171 txi1 = txif(tphi1), cxi1 = 1/hyp(txi1), sxi1 = txi1 * cxi1, 172 txi2 = txif(tphi2), cxi2 = 1/hyp(txi2), sxi2 = txi2 * cxi2, 173 dtbet2 = _fm * (tbet1 + tbet2), 174 es1 = 1 - _e2 * Math::sq(sphi1), es2 = 1 - _e2 * Math::sq(sphi2), 175 /* 176 dsxi = ( (_e2 * sq(sphi2 + sphi1) + es2 + es1) / (2 * es2 * es1) + 177 Datanhee(sphi2, sphi1) ) * Dsn(tphi2, tphi1, sphi2, sphi1) / 178 ( 2 * _qx ), 179 */ 180 dsxi = ( (1 + _e2 * sphi1 * sphi2) / (es2 * es1) + 181 Datanhee(sphi2, sphi1) ) * Dsn(tphi2, tphi1, sphi2, sphi1) / 182 ( 2 * _qx ), 183 den = (sxi2 + sxi1) * dtbet2 + (scbet22 + scbet12) * dsxi, 184 // s = (sq(tbet2) - sq(tbet1)) / (scbet22*sxi2 - scbet12*sxi1) 185 s = 2 * dtbet2 / den, 186 // 1-s = -(sq(scbet2)*(1-sxi2) - sq(scbet1)*(1-sxi1)) / 187 // (scbet22*sxi2 - scbet12*sxi1) 188 // Write 189 // sq(scbet)*(1-sxi) = sq(scbet)*(1-sphi) * (1-sxi)/(1-sphi) 190 sm1 = -Dsn(tphi2, tphi1, sphi2, sphi1) * 191 ( -( ((sphi2 <= 0 ? (1 - sxi2) / (1 - sphi2) : 192 Math::sq(cxi2/cphi2) * (1 + sphi2) / (1 + sxi2)) + 193 (sphi1 <= 0 ? (1 - sxi1) / (1 - sphi1) : 194 Math::sq(cxi1/cphi1) * (1 + sphi1) / (1 + sxi1))) ) * 195 (1 + _e2 * (sphi1 + sphi2 + sphi1 * sphi2)) / 196 (1 + (sphi1 + sphi2 + sphi1 * sphi2)) + 197 (scbet22 * (sphi2 <= 0 ? 1 - sphi2 : 198 Math::sq(cphi2) / ( 1 + sphi2)) + 199 scbet12 * (sphi1 <= 0 ? 1 - sphi1 : Math::sq(cphi1) / ( 1 + sphi1))) 200 * (_e2 * (1 + sphi1 + sphi2 + _e2 * sphi1 * sphi2)/(es1 * es2) 201 +_e2m * DDatanhee(sphi1, sphi2) ) / _qZ ) / den; 202 // C = (scbet22*sxi2 - scbet12*sxi1) / (scbet22 * scbet12 * (sx2 - sx1)) 203 C = den / (2 * scbet12 * scbet22 * dsxi); 204 tphi0 = (tphi2 + tphi1)/2; 205 real stol = tol0_ * max(real(1), abs(tphi0)); 206 for (int i = 0; i < 2*numit0_ || GEOGRAPHICLIB_PANIC; ++i) { 207 // Solve (scbet0^2 * sphi0) / (1/qZ + scbet0^2 * sphi0 * sxi0) = s 208 // for tphi0 by Newton's method on 209 // v(tphi0) = (scbet0^2 * sphi0) - s * (1/qZ + scbet0^2 * sphi0 * sxi0) 210 // = 0 211 // Alt: 212 // (scbet0^2 * sphi0) / (1/qZ - scbet0^2 * sphi0 * (1-sxi0)) 213 // = s / (1-s) 214 // w(tphi0) = (1-s) * (scbet0^2 * sphi0) 215 // - s * (1/qZ - scbet0^2 * sphi0 * (1-sxi0)) 216 // = (1-s) * (scbet0^2 * sphi0) 217 // - S/qZ * (1 - scbet0^2 * sphi0 * (qZ-q0)) 218 // Now 219 // qZ-q0 = (1+e2*sphi0)*(1-sphi0)/(1-e2*sphi0^2) + 220 // (1-e2)*atanhee((1-sphi0)/(1-e2*sphi0)) 221 // In limit sphi0 -> 1, qZ-q0 -> 2*(1-sphi0)/(1-e2), so wrte 222 // qZ-q0 = 2*(1-sphi0)/(1-e2) + A + B 223 // A = (1-sphi0)*( (1+e2*sphi0)/(1-e2*sphi0^2) - (1+e2)/(1-e2) ) 224 // = -e2 *(1-sphi0)^2 * (2+(1+e2)*sphi0) / ((1-e2)*(1-e2*sphi0^2)) 225 // B = (1-e2)*atanhee((1-sphi0)/(1-e2*sphi0)) - (1-sphi0) 226 // = (1-sphi0)*(1-e2)/(1-e2*sphi0)* 227 // ((atanhee(x)/x-1) - e2*(1-sphi0)/(1-e2)) 228 // x = (1-sphi0)/(1-e2*sphi0), atanhee(x)/x = atanh(e*x)/(e*x) 229 // 230 // 1 - scbet0^2 * sphi0 * (qZ-q0) 231 // = 1 - scbet0^2 * sphi0 * (2*(1-sphi0)/(1-e2) + A + B) 232 // = D - scbet0^2 * sphi0 * (A + B) 233 // D = 1 - scbet0^2 * sphi0 * 2*(1-sphi0)/(1-e2) 234 // = (1-sphi0)*(1-e2*(1+2*sphi0*(1+sphi0)))/((1-e2)*(1+sphi0)) 235 // dD/dsphi0 = -2*(1-e2*sphi0^2*(2*sphi0+3))/((1-e2)*(1+sphi0)^2) 236 // d(A+B)/dsphi0 = 2*(1-sphi0^2)*e2*(2-e2*(1+sphi0^2))/ 237 // ((1-e2)*(1-e2*sphi0^2)^2) 238 239 real 240 scphi02 = 1 + Math::sq(tphi0), scphi0 = sqrt(scphi02), 241 // sphi0m = 1-sin(phi0) = 1/( sec(phi0) * (tan(phi0) + sec(phi0)) ) 242 sphi0 = tphi0 / scphi0, sphi0m = 1/(scphi0 * (tphi0 + scphi0)), 243 // scbet0^2 * sphi0 244 g = (1 + Math::sq( _fm * tphi0 )) * sphi0, 245 // dg/dsphi0 = dg/dtphi0 * scphi0^3 246 dg = _e2m * scphi02 * (1 + 2 * Math::sq(tphi0)) + _e2, 247 D = sphi0m * (1 - _e2*(1 + 2*sphi0*(1+sphi0))) / (_e2m * (1+sphi0)), 248 // dD/dsphi0 249 dD = -2 * (1 - _e2*Math::sq(sphi0) * (2*sphi0+3)) / 250 (_e2m * Math::sq(1+sphi0)), 251 A = -_e2 * Math::sq(sphi0m) * (2+(1+_e2)*sphi0) / 252 (_e2m*(1-_e2*Math::sq(sphi0))), 253 B = (sphi0m * _e2m / (1 - _e2*sphi0) * 254 (atanhxm1(_e2 * 255 Math::sq(sphi0m / (1-_e2*sphi0))) - _e2*sphi0m/_e2m)), 256 // d(A+B)/dsphi0 257 dAB = (2 * _e2 * (2 - _e2 * (1 + Math::sq(sphi0))) / 258 (_e2m * Math::sq(1 - _e2*Math::sq(sphi0)) * scphi02)), 259 u = sm1 * g - s/_qZ * ( D - g * (A + B) ), 260 // du/dsphi0 261 du = sm1 * dg - s/_qZ * (dD - dg * (A + B) - g * dAB), 262 dtu = -u/du * (scphi0 * scphi02); 263 tphi0 += dtu; 264 if (!(abs(dtu) >= stol)) 265 break; 266 } 267 } 268 _txi0 = txif(tphi0); _scxi0 = hyp(_txi0); _sxi0 = _txi0 / _scxi0; 269 _n0 = tphi0/hyp(tphi0); 270 _m02 = 1 / (1 + Math::sq(_fm * tphi0)); 271 _nrho0 = polar ? 0 : _a * sqrt(_m02); 272 _k0 = sqrt(tphi1 == tphi2 ? 1 : C / (_m02 + _n0 * _qZ * _sxi0)) * k1; 273 _k2 = Math::sq(_k0); 274 _lat0 = _sign * atan(tphi0)/Math::degree(); 275 } 276 CylindricalEqualArea()277 const AlbersEqualArea& AlbersEqualArea::CylindricalEqualArea() { 278 static const AlbersEqualArea 279 cylindricalequalarea(Constants::WGS84_a(), Constants::WGS84_f(), 280 real(0), real(1), real(0), real(1), real(1)); 281 return cylindricalequalarea; 282 } 283 AzimuthalEqualAreaNorth()284 const AlbersEqualArea& AlbersEqualArea::AzimuthalEqualAreaNorth() { 285 static const AlbersEqualArea 286 azimuthalequalareanorth(Constants::WGS84_a(), Constants::WGS84_f(), 287 real(1), real(0), real(1), real(0), real(1)); 288 return azimuthalequalareanorth; 289 } 290 AzimuthalEqualAreaSouth()291 const AlbersEqualArea& AlbersEqualArea::AzimuthalEqualAreaSouth() { 292 static const AlbersEqualArea 293 azimuthalequalareasouth(Constants::WGS84_a(), Constants::WGS84_f(), 294 real(-1), real(0), real(-1), real(0), real(1)); 295 return azimuthalequalareasouth; 296 } 297 txif(real tphi) const298 Math::real AlbersEqualArea::txif(real tphi) const { 299 // sxi = ( sphi/(1-e2*sphi^2) + atanhee(sphi) ) / 300 // ( 1/(1-e2) + atanhee(1) ) 301 // 302 // txi = ( sphi/(1-e2*sphi^2) + atanhee(sphi) ) / 303 // sqrt( ( (1+e2*sphi)*(1-sphi)/( (1-e2*sphi^2) * (1-e2) ) + 304 // atanhee((1-sphi)/(1-e2*sphi)) ) * 305 // ( (1-e2*sphi)*(1+sphi)/( (1-e2*sphi^2) * (1-e2) ) + 306 // atanhee((1+sphi)/(1+e2*sphi)) ) ) 307 // = ( tphi/(1-e2*sphi^2) + atanhee(sphi, e2)/cphi ) / 308 // sqrt( 309 // ( (1+e2*sphi)/( (1-e2*sphi^2) * (1-e2) ) + Datanhee(1, sphi) ) * 310 // ( (1-e2*sphi)/( (1-e2*sphi^2) * (1-e2) ) + Datanhee(1, -sphi) ) ) 311 // 312 // This function maintains odd parity 313 real 314 cphi = 1 / sqrt(1 + Math::sq(tphi)), 315 sphi = tphi * cphi, 316 es1 = _e2 * sphi, 317 es2m1 = 1 - es1 * sphi, // 1 - e2 * sphi^2 318 es2m1a = _e2m * es2m1; // (1 - e2 * sphi^2) * (1 - e2) 319 return ( tphi / es2m1 + atanhee(sphi) / cphi ) / 320 sqrt( ( (1 + es1) / es2m1a + Datanhee(1, sphi) ) * 321 ( (1 - es1) / es2m1a + Datanhee(1, -sphi) ) ); 322 } 323 tphif(real txi) const324 Math::real AlbersEqualArea::tphif(real txi) const { 325 real 326 tphi = txi, 327 stol = tol_ * max(real(1), abs(txi)); 328 // CHECK: min iterations = 1, max iterations = 2; mean = 1.99 329 for (int i = 0; i < numit_ || GEOGRAPHICLIB_PANIC; ++i) { 330 // dtxi/dtphi = (scxi/scphi)^3 * 2*(1-e^2)/(qZ*(1-e^2*sphi^2)^2) 331 real 332 txia = txif(tphi), 333 tphi2 = Math::sq(tphi), 334 scphi2 = 1 + tphi2, 335 scterm = scphi2/(1 + Math::sq(txia)), 336 dtphi = (txi - txia) * scterm * sqrt(scterm) * 337 _qx * Math::sq(1 - _e2 * tphi2 / scphi2); 338 tphi += dtphi; 339 if (!(abs(dtphi) >= stol)) 340 break; 341 } 342 return tphi; 343 } 344 345 // return atanh(sqrt(x))/sqrt(x) - 1 = x/3 + x^2/5 + x^3/7 + ... 346 // typical x < e^2 = 2*f atanhxm1(real x)347 Math::real AlbersEqualArea::atanhxm1(real x) { 348 real s = 0; 349 if (abs(x) < real(0.5)) { 350 static const real lg2eps_ = -log2(numeric_limits<real>::epsilon() / 2); 351 int e; 352 frexp(x, &e); 353 e = -e; 354 // x = [0.5,1) * 2^(-e) 355 // estimate n s.t. x^n/(2*n+1) < x/3 * epsilon/2 356 // a stronger condition is x^(n-1) < epsilon/2 357 // taking log2 of both sides, a stronger condition is 358 // (n-1)*(-e) < -lg2eps or (n-1)*e > lg2eps or n > ceiling(lg2eps/e)+1 359 int n = x == 0 ? 1 : int(ceil(lg2eps_ / e)) + 1; 360 while (n--) // iterating from n-1 down to 0 361 s = x * s + (n ? 1 : 0)/Math::real(2*n + 1); 362 } else { 363 real xs = sqrt(abs(x)); 364 s = (x > 0 ? atanh(xs) : atan(xs)) / xs - 1; 365 } 366 return s; 367 } 368 369 // return (Datanhee(1,y) - Datanhee(1,x))/(y-x) DDatanhee(real x,real y) const370 Math::real AlbersEqualArea::DDatanhee(real x, real y) const { 371 // This function is called with x = sphi1, y = sphi2, phi1 <= phi2, sphi2 372 // >= 0, abs(sphi1) <= phi2. However for safety's sake we enforce x <= y. 373 if (y < x) swap(x, y); // ensure that x <= y 374 real q1 = abs(_e2), 375 q2 = abs(2 * _e / _e2m * (1 - x)); 376 return 377 x <= 0 || !(min(q1, q2) < real(0.75)) ? DDatanhee0(x, y) : 378 (q1 < q2 ? DDatanhee1(x, y) : DDatanhee2(x, y)); 379 } 380 381 // Rearrange difference so that 1 - x is in the denominator, then do a 382 // straight divided difference. DDatanhee0(real x,real y) const383 Math::real AlbersEqualArea::DDatanhee0(real x, real y) const { 384 return (Datanhee(1, y) - Datanhee(x, y))/(1 - x); 385 } 386 387 // The expansion for e2 small DDatanhee1(real x,real y) const388 Math::real AlbersEqualArea::DDatanhee1(real x, real y) const { 389 // The series in e2 is 390 // sum( c[l] * e2^l, l, 1, N) 391 // where 392 // c[l] = sum( x^i * y^j; i >= 0, j >= 0, i+j < 2*l) / (2*l + 1) 393 // = ( (x-y) - (1-y) * x^(2*l+1) + (1-x) * y^(2*l+1) ) / 394 // ( (2*l+1) * (x-y) * (1-y) * (1-x) ) 395 // For x = y = 1, c[l] = l 396 // 397 // In the limit x,y -> 1, 398 // 399 // DDatanhee -> e2/(1-e2)^2 = sum(l * e2^l, l, 1, inf) 400 // 401 // Use if e2 is sufficiently small. 402 real s = 0; 403 real z = 1, k = 1, t = 0, c = 0, en = 1; 404 while (true) { 405 t = y * t + z; c += t; z *= x; 406 t = y * t + z; c += t; z *= x; 407 k += 2; en *= _e2; 408 // Here en[l] = e2^l, k[l] = 2*l + 1, 409 // c[l] = sum( x^i * y^j; i >= 0, j >= 0, i+j < 2*l) / (2*l + 1) 410 // Taylor expansion is 411 // s = sum( c[l] * e2^l, l, 1, N) 412 real ds = en * c / k; 413 s += ds; 414 if (!(abs(ds) > abs(s) * eps_/2)) 415 break; // Iterate until the added term is sufficiently small 416 } 417 return s; 418 } 419 420 // The expansion for x (and y) close to 1 DDatanhee2(real x,real y) const421 Math::real AlbersEqualArea::DDatanhee2(real x, real y) const { 422 // If x and y are both close to 1, expand in Taylor series in dx = 1-x and 423 // dy = 1-y: 424 // 425 // DDatanhee = sum(C_m * (dx^(m+1) - dy^(m+1)) / (dx - dy), m, 0, inf) 426 // 427 // where 428 // 429 // C_m = sum( (m+2)!! / (m+2-2*k)!! * 430 // ((m+1)/2)! / ((m+1)/2-k)! / 431 // (k! * (2*k-1)!!) * 432 // e2^((m+1)/2+k), 433 // k, 0, (m+1)/2) * (-1)^m / ((m+2) * (1-e2)^(m+2)) 434 // for m odd, and 435 // 436 // C_m = sum( 2 * (m+1)!! / (m+1-2*k)!! * 437 // (m/2+1)! / (m/2-k)! / 438 // (k! * (2*k+1)!!) * 439 // e2^(m/2+1+k), 440 // k, 0, m/2)) * (-1)^m / ((m+2) * (1-e2)^(m+2)) 441 // for m even. 442 // 443 // Here i!! is the double factorial extended to negative i with 444 // i!! = (i+2)!!/(i+2). 445 // 446 // Note that 447 // (dx^(m+1) - dy^(m+1)) / (dx - dy) = 448 // dx^m + dx^(m-1)*dy ... + dx*dy^(m-1) + dy^m 449 // 450 // Leading (m = 0) term is e2 / (1 - e2)^2 451 // 452 // Magnitude of mth term relative to the leading term scales as 453 // 454 // 2*(2*e/(1-e2)*dx)^m 455 // 456 // So use series if (2*e/(1-e2)*dx) is sufficiently small 457 real s, dx = 1 - x, dy = 1 - y, xy = 1, yy = 1, ee = _e2 / Math::sq(_e2m); 458 s = ee; 459 for (int m = 1; ; ++m) { 460 real c = m + 2, t = c; 461 yy *= dy; // yy = dy^m 462 xy = dx * xy + yy; 463 // Now xy = dx^m + dx^(m-1)*dy ... + dx*dy^(m-1) + dy^m 464 // = (dx^(m+1) - dy^(m+1)) / (dx - dy) 465 // max value = (m+1) * max(dx,dy)^m 466 ee /= -_e2m; 467 if (m % 2 == 0) ee *= _e2; 468 // Now ee = (-1)^m * e2^(floor(m/2)+1) / (1-e2)^(m+2) 469 int kmax = (m+1)/2; 470 for (int k = kmax - 1; k >= 0; --k) { 471 // max coeff is less than 2^(m+1) 472 c *= (k + 1) * (2 * (k + m - 2*kmax) + 3); 473 c /= (kmax - k) * (2 * (kmax - k) + 1); 474 // Horner sum for inner _e2 series 475 t = _e2 * t + c; 476 } 477 // Straight sum for outer m series 478 real ds = t * ee * xy / (m + 2); 479 s = s + ds; 480 if (!(abs(ds) > abs(s) * eps_/2)) 481 break; // Iterate until the added term is sufficiently small 482 } 483 return s; 484 } 485 Forward(real lon0,real lat,real lon,real & x,real & y,real & gamma,real & k) const486 void AlbersEqualArea::Forward(real lon0, real lat, real lon, 487 real& x, real& y, real& gamma, real& k) const { 488 lon = Math::AngDiff(lon0, lon); 489 lat *= _sign; 490 real sphi, cphi; 491 Math::sincosd(Math::LatFix(lat) * _sign, sphi, cphi); 492 cphi = max(epsx_, cphi); 493 real 494 lam = lon * Math::degree(), 495 tphi = sphi/cphi, txi = txif(tphi), sxi = txi/hyp(txi), 496 dq = _qZ * Dsn(txi, _txi0, sxi, _sxi0) * (txi - _txi0), 497 drho = - _a * dq / (sqrt(_m02 - _n0 * dq) + _nrho0 / _a), 498 theta = _k2 * _n0 * lam, stheta = sin(theta), ctheta = cos(theta), 499 t = _nrho0 + _n0 * drho; 500 x = t * (_n0 != 0 ? stheta / _n0 : _k2 * lam) / _k0; 501 y = (_nrho0 * 502 (_n0 != 0 ? 503 (ctheta < 0 ? 1 - ctheta : Math::sq(stheta)/(1 + ctheta)) / _n0 : 504 0) 505 - drho * ctheta) / _k0; 506 k = _k0 * (t != 0 ? t * hyp(_fm * tphi) / _a : 1); 507 y *= _sign; 508 gamma = _sign * theta / Math::degree(); 509 } 510 Reverse(real lon0,real x,real y,real & lat,real & lon,real & gamma,real & k) const511 void AlbersEqualArea::Reverse(real lon0, real x, real y, 512 real& lat, real& lon, 513 real& gamma, real& k) const { 514 y *= _sign; 515 real 516 nx = _k0 * _n0 * x, ny = _k0 * _n0 * y, y1 = _nrho0 - ny, 517 den = hypot(nx, y1) + _nrho0, // 0 implies origin with polar aspect 518 drho = den != 0 ? (_k0*x*nx - 2*_k0*y*_nrho0 + _k0*y*ny) / den : 0, 519 // dsxia = scxi0 * dsxi 520 dsxia = - _scxi0 * (2 * _nrho0 + _n0 * drho) * drho / 521 (Math::sq(_a) * _qZ), 522 txi = (_txi0 + dsxia) / sqrt(max(1 - dsxia * (2*_txi0 + dsxia), epsx2_)), 523 tphi = tphif(txi), 524 theta = atan2(nx, y1), 525 lam = _n0 != 0 ? theta / (_k2 * _n0) : x / (y1 * _k0); 526 gamma = _sign * theta / Math::degree(); 527 lat = Math::atand(_sign * tphi); 528 lon = lam / Math::degree(); 529 lon = Math::AngNormalize(lon + Math::AngNormalize(lon0)); 530 k = _k0 * (den != 0 ? (_nrho0 + _n0 * drho) * hyp(_fm * tphi) / _a : 1); 531 } 532 SetScale(real lat,real k)533 void AlbersEqualArea::SetScale(real lat, real k) { 534 if (!(isfinite(k) && k > 0)) 535 throw GeographicErr("Scale is not positive"); 536 if (!(abs(lat) < 90)) 537 throw GeographicErr("Latitude for SetScale not in (-90d, 90d)"); 538 real x, y, gamma, kold; 539 Forward(0, lat, 0, x, y, gamma, kold); 540 k /= kold; 541 _k0 *= k; 542 _k2 = Math::sq(_k0); 543 } 544 545 } // namespace GeographicLib 546