19 , epsx_(
Math::sq(eps_))
20 , ahypover_(
Math::digits() *
log(
real(numeric_limits<
real>::radix)) + 2)
25 , _es((_f < 0 ? -1 : 1) *
sqrt(
abs(_e2)))
33 if (!(
abs(stdlat) <= 90))
37 Init(sphi, cphi, sphi, cphi,
k0);
44 , epsx_(
Math::sq(eps_))
45 , ahypover_(
Math::digits() *
log(
real(numeric_limits<
real>::radix)) + 2)
50 , _es((_f < 0 ? -1 : 1) *
sqrt(
abs(_e2)))
58 if (!(
abs(stdlat1) <= 90))
59 throw GeographicErr(
"Standard latitude 1 not in [-90d, 90d]");
60 if (!(
abs(stdlat2) <= 90))
61 throw GeographicErr(
"Standard latitude 2 not in [-90d, 90d]");
62 real sphi1, cphi1, sphi2, cphi2;
65 Init(sphi1, cphi1, sphi2, cphi2,
k1);
73 , epsx_(
Math::sq(eps_))
74 , ahypover_(
Math::digits() *
log(
real(numeric_limits<
real>::radix)) + 2)
79 , _es((_f < 0 ? -1 : 1) *
sqrt(
abs(_e2)))
88 throw GeographicErr(
"Standard latitude 1 not in [-90d, 90d]");
90 throw GeographicErr(
"Standard latitude 2 not in [-90d, 90d]");
91 if (!(
abs(sinlat1) <= 1 && coslat1 <= 1) || (coslat1 == 0 && sinlat1 == 0))
92 throw GeographicErr(
"Bad sine/cosine of standard latitude 1");
93 if (!(
abs(sinlat2) <= 1 && coslat2 <= 1) || (coslat2 == 0 && sinlat2 == 0))
94 throw GeographicErr(
"Bad sine/cosine of standard latitude 2");
95 if (coslat1 == 0 || coslat2 == 0)
96 if (!(coslat1 == coslat2 && sinlat1 == sinlat2))
98 (
"Standard latitudes must be equal is either is a pole");
99 Init(sinlat1, coslat1, sinlat2, coslat2,
k1);
107 sphi1 /= r; cphi1 /= r;
109 sphi2 /= r; cphi2 /= r;
111 bool polar = (cphi1 == 0);
115 _sign = sphi1 + sphi2 >= 0 ? 1 : -1;
119 swap(sphi1, sphi2);
swap(cphi1, cphi2);
122 tphi1 = sphi1/cphi1, tphi2 = sphi2/cphi2, tphi0;
142 tbet1 =
_fm * tphi1, scbet1 =
hyp(tbet1),
143 tbet2 =
_fm * tphi2, scbet2 =
hyp(tbet2);
147 tchi1 = chxi1 * tphi1 - shxi1 * scphi1, scchi1 =
hyp(tchi1),
150 tchi2 = chxi2 * tphi2 - shxi2 * scphi2, scchi2 =
hyp(tchi2),
152 if (tphi2 - tphi1 != 0) {
156 *
Dhyp(tbet2, tbet1, scbet2, scbet1) *
_fm;
158 real den =
Dasinh(tphi2, tphi1, scphi2, scphi1)
159 -
Deatanhe(sphi2, sphi1) *
Dsn(tphi2, tphi1, sphi2, sphi1);
185 s1 = (tphi1 * (2 * shxi1 * chxi1 * scphi1 -
_e2 * tphi1) -
187 s2 = (tphi2 * (2 * shxi2 * chxi2 * scphi2 -
_e2 * tphi2) -
190 t1 = tchi1 < 0 ? scbet1 - tchi1 : (s1 + 1)/(scbet1 + tchi1),
191 t2 = tchi2 < 0 ? scbet2 - tchi2 : (s2 + 1)/(scbet2 + tchi2),
192 a2 = -(s2 / (scbet2 + scchi2) + t2) / (2 * scbet2),
193 a1 = -(s1 / (scbet1 + scchi1) + t1) / (2 * scbet1);
197 t *= ( ( (tchi2 >= 0 ? scchi2 + tchi2 : 1/(scchi2 - tchi2)) +
198 (tchi1 >= 0 ? scchi1 + tchi1 : 1/(scchi1 - tchi1)) ) /
199 (4 * scbet1 * scbet2) ) *
_fm;
206 real tbm = ( ((tbet1 > 0 ? 1/(scbet1+tbet1) : scbet1 - tbet1) +
207 (tbet2 > 0 ? 1/(scbet2+tbet2) : scbet2 - tbet2)) /
219 dtchi = den /
Dasinh(tchi2, tchi1, scchi2, scchi1),
221 dbet = (
_e2/
_fm) * ( 1 / (scbet2 +
_fm * scphi2) +
222 1 / (scbet1 +
_fm * scphi1) );
238 shxiZ =
sinh(xiZ), chxiZ =
hyp(shxiZ),
241 dxiZ1 =
Deatanhe(
real(1), sphi1)/(scphi1*(tphi1+scphi1)),
242 dxiZ2 =
Deatanhe(
real(1), sphi2)/(scphi2*(tphi2+scphi2)),
243 dshxiZ1 =
Dsinh(xiZ, xi1, shxiZ, shxi1, chxiZ, chxi1) * dxiZ1,
244 dshxiZ2 =
Dsinh(xiZ, xi2, shxiZ, shxi2, chxiZ, chxi2) * dxiZ2,
245 dchxiZ1 =
Dhyp(shxiZ, shxi1, chxiZ, chxi1) * dshxiZ1,
246 dchxiZ2 =
Dhyp(shxiZ, shxi2, chxiZ, chxi2) * dshxiZ2,
248 amu12 = (- scphi1 * dchxiZ1 + tphi1 * dshxiZ1
249 - scphi2 * dchxiZ2 + tphi2 * dshxiZ2),
251 dxi =
Deatanhe(sphi1, sphi2) *
Dsn(tphi2, tphi1, sphi2, sphi1),
254 ( (
_f * 4 * scphi2 * dshxiZ2 >
_f * scphi1 * dshxiZ1 ?
256 (dshxiZ1 + dshxiZ2)/2 *
Dhyp(tphi1, tphi2, scphi1, scphi2)
257 - ( (scphi1 + scphi2)/2
258 *
Dsinh(xi1, xi2, shxi1, shxi2, chxi1, chxi2) * dxi ) :
260 (scphi2 * dshxiZ2 - scphi1 * dshxiZ1)/(tphi2 - tphi1))
261 + ( (tphi1 + tphi2)/2 *
Dhyp(shxi1, shxi2, chxi1, chxi2)
262 *
Dsinh(xi1, xi2, shxi1, shxi2, chxi1, chxi2) * dxi )
263 - (dchxiZ1 + dchxiZ2)/2 ),
265 dchia = (amu12 - dnu12 * (scphi2 + scphi1)),
266 tam = (dchia - dtchi * dbet) / (scchi1 + scchi2);
297 * (tchi1 >= 0 ? scchi1 + tchi1 : 1 / (scchi1 - tchi1));
304 * (tchi1 >= 0 ? scchi1 + tchi1 : 1 / (scchi1 - tchi1)) /
310 sphi = -1, cphi =
epsx_,
313 tchi =
hyp(shxi) * tphi - shxi * scphi, scchi =
hyp(tchi),
318 (tchi > 0 ? 1/(scchi + tchi) : (scchi - tchi))
350 tphi = sphi/cphi, scbet =
hyp(
_fm * tphi),
352 tchi =
hyp(shxi) * tphi - shxi * scphi, scchi =
hyp(tchi),
354 theta =
_n * lam, stheta =
sin(theta), ctheta =
cos(theta),
356 drho = -
_scale * (2 *
_nc < 1 && dpsi != 0 ?
358 (tchi > 0 ? 1/(scchi + tchi) : (scchi - tchi))
364 (ctheta < 0 ? 1 - ctheta :
Math::sq(stheta)/(1 + ctheta)) /
_n : 0)
368 * (tchi >= 0 ? scchi + tchi : 1 / (scchi - tchi)) / (
_scchi0 +
_tchi0));
395 ? (
x*nx +
y * (ny - 2*
_nrho0)) / den
402 dpsi = (den == 0 ? 0 :
420 tn = tnm1 + 1 == 0 ?
epsx_ : tnm1 + 1,
423 tchi = sh * (tn + 1/tn)/2 -
hyp(sh) * (tnm1 * (tn + 1)/tn)/2;
430 scbet =
hyp(
_fm * tphi), scchi =
hyp(tchi),
437 * (tchi >= 0 ? scchi + tchi : 1 / (scchi - tchi)) / (
_scchi0 +
_tchi0));
445 throw GeographicErr(
"Latitude for SetScale not in [-90d, 90d]");
447 throw GeographicErr(
"Incompatible polar latitude in SetScale");