From b06c771637b7add0e186a167837dfd55571c7864 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 18 Feb 2010 17:09:16 +0000
Subject: [PATCH] stable implementation of sinc
[[Imported from SVN: r5588]]
---
src/unitvector.hh | 7 ++++++-
1 file changed, 6 insertions(+), 1 deletion(-)
diff --git a/src/unitvector.hh b/src/unitvector.hh
index aa6db5d7..6ca22338 100644
--- a/src/unitvector.hh
+++ b/src/unitvector.hh
@@ -6,6 +6,11 @@
template <int dim>
class UnitVector
{
+ /** \brief Computes sin(x/2) / x without getting unstable for small x */
+ static double sinc(const double& x) {
+ return (x < 1e-4) ? 1 + (x*x/6) : std::sin(x)/x;
+ }
+
public:
/** \brief The type used for coordinates */
@@ -29,7 +34,7 @@ public:
const double norm = v.two_norm();
UnitVector result = p;
result.data_ *= std::cos(norm);
- result.data_.axpy(std::sin(norm)/norm, v);
+ result.data_.axpy(sinc(norm), v);
return result;
}
--
GitLab