[algorithms] Improve the robustness of thomas_inverse formula.

For lat equal to pi/2 or -pi/2 avoid calculating "infinite" tangent result.
This commit is contained in:
Adam Wulkiewicz
2015-01-25 19:49:44 +01:00
parent 8c04a59076
commit 4520fbc373

View File

@@ -59,10 +59,18 @@ public:
CT const one_minus_f = CT(1) - f;
CT const tan_theta1 = one_minus_f * tan(lat1);
CT const tan_theta2 = one_minus_f * tan(lat2);
CT const theta1 = atan(tan_theta1);
CT const theta2 = atan(tan_theta2);
// CT const tan_theta1 = one_minus_f * tan(lat1);
// CT const tan_theta2 = one_minus_f * tan(lat2);
// CT const theta1 = atan(tan_theta1);
// CT const theta2 = atan(tan_theta2);
CT const pi_half = math::pi<CT>() / CT(2);
CT const theta1 = math::equals(lat1, pi_half) ? lat1 :
math::equals(lat1, -pi_half) ? lat1 :
atan(one_minus_f * tan(lat1));
CT const theta2 = math::equals(lat2, pi_half) ? lat2 :
math::equals(lat2, -pi_half) ? lat2 :
atan(one_minus_f * tan(lat2));
CT const theta_m = (theta1 + theta2) / CT(2);
CT const d_theta_m = (theta2 - theta1) / CT(2);