Phi*ToR(): rounding-errors will be attended, because the domain of atanh is (-1,1)

AMDiS/src/MathFunctions.h

@@ -16,8 +16,12 @@ namespace AMDiS {
   }
 
   //convert Phi1 to r
-  inline double Phi1ToR(double p1, double eps) { return atanh(1-2*p1)*eps/3; }
+  inline double Phi1ToR(double p1, double eps) {
+	return eps/3.0 * atanh( std::max(-1+1.e-14, std::min(1-1.e-14, 1-2*p1)) );
+  }
   //convert Phi2 to r
-  inline double Phi2ToR(double p2, double eps) { return atanh(1+2*p2)*eps/3; }
+  inline double Phi2ToR(double p2, double eps) {
+	return eps/3.0 * atanh( std::max(-1+1.e-14, std::min(1-1.e-14, 1+2*p2) );
+  }
 }
 #endif