adding implementation for 3- and 4-point functions

still need to be tested
FlexibleSUSY · Mar 1, 2018 · efc64b5 · efc64b5
1 parent 58eac5f
commit efc64b5
Show file tree

Hide file tree

Showing 2 changed files with 112 additions and 0 deletions.
diff --git a/src/pv2.cpp b/src/pv2.cpp
@@ -252,4 +252,109 @@ double h0(double p2, double m12, double m22, double q2) noexcept
    return 4. * b22(p2,m12,m22,q2) + g0(p2,m12,m22,q2);
 }
 
+double c0(double m12, double m22, double m32) noexcept
+{
+  using std::log;
+
+  double ans = 0.;
+
+  if (is_close(m12,0.,EPSTOL) && is_close(m22,0.,EPSTOL) && is_close(m32,0.,EPSTOL)) {
+     // c0 is undefined for m1 == m2 == m3 == 0
+     ans = 0.;
+  } else if (is_close(m22,0.,EPSTOL) && is_close(m32,0.,EPSTOL)) {
+     // c0 is undefined for m2 == m3 == 0
+     ans = 0.;
+  } else if (is_close(m12,0.,EPSTOL) && is_close(m32,0.,EPSTOL)) {
+     // c0 is undefined for m1 == m3 == 0
+     ans = 0.;
+  } else if (is_close(m12,0.,EPSTOL) && is_close(m22,0.,EPSTOL)) {
+     // c0 is undefined for m1 == m2 == 0
+     ans= 0.;
+  } else if (is_close(m12,0.,EPSTOL)) {
+     if (is_close(m22,m32,EPSTOL)) {
+        ans = -1./m22;
+     } else {
+        ans = log(m32/m22)/(m22 - m32);
+     }
+  } else if (is_close(m22,0.,EPSTOL)) {
+     if (is_close(m12,m32,EPSTOL)) {
+        ans = -1./m12;
+     } else {
+        ans = log(m32/m12)/(m12 - m32);
+     }
+  } else if (is_close(m32,0.,EPSTOL)) {
+     if (is_close(m12,m22,EPSTOL)) {
+        ans = -1./m12;
+     } else {
+        ans = log(m22/m12)/(m12 - m22);
+     }
+  } else if (is_close(m22, m32, EPSTOL)) {
+    if (is_close(m12, m22, EPSTOL)) {
+      ans = ( - 0.5 / m22 );
+    } else {
+      ans = ( m12 / sqr(m12-m22) * log(m22/m12)
+               + 1.0 / (m12 - m22) );
+    }
+  } else if (is_close(m12, m22, EPSTOL)) {
+     ans = ( - ( 1.0 + m32 / (m22-m32) * log(m32/m22) )
+             / (m22-m32) );
+  } else if (is_close(m12, m32, EPSTOL)) {
+     ans = ( - (1.0 + m22 / (m32-m22) * log(m22/m32))
+             / (m32-m22) );
+  } else {
+     ans = (1.0 / (m22 - m32) *
+            (m22 / (m12 - m22) *
+             log(m22 / m12) -
+             m32 / (m12 - m32) *
+             log(m32 / m12)) );
+  }
+
+  return ans;
+}
+
+double d0(double m12, double m22, double m32, double m42) noexcept
+{
+  using std::log;
+
+  if (is_close(m12, m22, EPSTOL)) {
+    if (is_close(m22,0.,EPSTOL)) {
+       // d0 is undefined for m1 == m2 == 0
+       return 0.;
+    } else if (is_close(m32,0.,EPSTOL)) {
+       return (-m22 + m42 - m22 * log(m42/m22))/
+          (m22 * sqr(m22 - m42));
+    } else if (is_close(m42,0.,EPSTOL)) {
+       return (-m22 + m32 - m22 * log(m32/m22))/
+          (m22 * sqr(m22 - m32));
+    } else if (is_close(m22, m32, EPSTOL) && is_close(m22, m42, EPSTOL)) {
+      return 1.0 / (6.0 * sqr(m22));
+    } else if (is_close(m22, m32, EPSTOL)) {
+      return (sqr(m22) - sqr(m42) + 2.0 * m42 * m22 * log(m42 / m22)) /
+	(2.0 * m22 * sqr(m22 - m42) * (m22 - m42));
+    } else if (is_close(m22, m42, EPSTOL)) {
+      return (sqr(m22) - sqr(m32) + 2.0 * m32 * m22 * log(m32 / m22)) /
+	(2.0 * m22 * sqr(m22 - m32) * (m22 - m32));
+    } else if (is_close(m32, m42, EPSTOL)) {
+      return -1.0 / sqr(m22 - m32) *
+	((m22 + m32) / (m22 - m32) * log(m32 / m22) + 2.0);
+    }
+
+    return
+      (m42 / sqr(m22 - m42) * log(m42 / m22) +
+       m42 / (m22 * (m22 - m42)) -
+       m32 / sqr(m22 - m32) * log(m32 / m22) -
+       m32 / (m22 * (m22 - m32))) / (m32 - m42);
+  }
+  return (c0(m12, m32, m42) - c0(m22, m32, m42)) / (m12 - m22);
+}
+
+double d27(double m12, double m22, double m32, double m42) noexcept
+{
+  if (is_close(m12, m22, EPSTOL))
+    m12 += TOL * 0.01;
+
+  return (m12 * c0(m12, m32, m42) - m22 * c0(m22, m32, m42))
+    / (4.0 * (m12 - m22));
+}
+
 } // namespace flexiblesusy
diff --git a/src/pv2.hpp b/src/pv2.hpp
@@ -37,6 +37,13 @@ double f0(double p2, double m12, double m22, double q2) noexcept;
 double g0(double p2, double m12, double m22, double q2) noexcept;
 double h0(double p2, double m12, double m22, double q2) noexcept;
 
+// 3-point function for p^2 = 0
+double c0(double m12, double m22, double m32) noexcept;
+
+// 4-point function for p^2 = 0
+double d27(double m12, double m22, double m32, double m42) noexcept;
+double d0(double m12, double m22, double m32, double m42) noexcept;
+
 } // namespace flexiblesusy
 
 #endif