From 60aeab2239d47a5c0b68c5aa20a5596b90180bde Mon Sep 17 00:00:00 2001
From: Andreas Suter <andreas.suter@psi.ch>
Date: Wed, 28 May 2025 12:46:37 +0200
Subject: [PATCH] changed the integral approximation approach for
 Gaussian/Lorentzian.

---
 src/external/LF_GL/PLGKT_LF.cpp | 43 +++++++++------------------------
 1 file changed, 11 insertions(+), 32 deletions(-)

diff --git a/src/external/LF_GL/PLGKT_LF.cpp b/src/external/LF_GL/PLGKT_LF.cpp
index 66023aaf..63bb458f 100644
--- a/src/external/LF_GL/PLGKT_LF.cpp
+++ b/src/external/LF_GL/PLGKT_LF.cpp
@@ -31,6 +31,7 @@
 #include <cmath>
 #include <algorithm>
 #include <numbers>
+#include <chrono>
 
 #include "PLGKT_LF.h"
 #include "PGKT_LF.h"
@@ -49,43 +50,17 @@ PLGKT_LF::PLGKT_LF(std::vector<double> &param, const double tmax) : fParam(param
  * \f$r = \frac{\Delta_{\rm G}}{\Delta_{\rm L}}\f$
  * (see A. Yaouanc and P. Dalmas de Reotier, ``Muon Spin Rotation, Relaxation, and Resonance'', p.129)
  *
- * <p>Integration method: Simpson on the \f$r-\f$ intervals:
- * \f$[0.2, 0.6], [0.6, 1.0], [1.0, 3.0], [3.0, 5.0], [5.0, 7.5], [7.5, 10.0], [10.0, 55.0], [55.0, 100.0]\f$
+ * <p>Integration method is described in the docu directory.
  *
- * <p>This leads to
- * \f{eqnarray*}{
- *   P_z(t, \Delta_{\rm L}, B, \nu) &=& \frac{0.4}{6} \left[f(0.2) + 4 f(0.4) + f(0.6) \right] +
- *                                      \frac{0.4}{6} \left[f(0.6) + 4 f(0.8) + f(1.0) \right] \\
- *                                  & & \frac{2}{6} \left[f(1.0) + 4 f(2.0) + f(3.0) \right] +
- *                                      \frac{2}{6} \left[f(3.0) + 4 f(4.0) + f(5.0) \right] \\
- *                                  & & \frac{2.5}{6} \left[f(5.0) + 4 f(6.25) + f(7.5) \right] +
- *                                      \frac{2.5}{6} \left[f(7.5) + 4 f(8.75) + f(10.0) \right] \\
- *                                  & & \frac{45}{6} \left[f(10.0) + 4 f(32.5) + f(55.0) \right] +
- *                                      \frac{45}{6} \left[f(55.0) + 4 f(77.5) + f(100.0) \right] \\
- *                                  &=& \frac{1}{6} \left[
- *                                      0.4 f(0.2) + 1.6 f(0.4) + 0.8 f(0.6) + 1.6 f(0.8) + 2.4 f(1.0) +\\
- *                                  & & 8.0 f(2.0) + 4 f(3.0) + 8 f(4.0) + 4.5 f(5.0) + 10.0 f(6.25) +\\
- *                                  & & 5 f(7.5) + 10.0 f(8.75) + 47.5 f(10.0) + 180.0 f(32.5) + 90 f(55.0) +\\
- *                                  & & 180 f(77.5) + 45 f(100.0)
- *                                  \right].
- * \f}
- *
- * <p>where
- * \f[
- *   f(r) = P_z(t, \Delta_{\rm G} = r\cdot\Delta_{\rm L}, B, \nu) \cdot \sqrt{\frac{2}{\pi}} \frac{1}{r^2} \, \exp(-\frac{1}{2 r^2})
- * \f]
  * @return 0 on success, >0 otherwise
  */
 int PLGKT_LF::DynamicLGKTLF()
 {
-//  std::vector<double> rr={0.2, 0.4, 0.6, 0.8, 1.0, 2.0, 3.0, 4.0, 5.0, 6.25, 7.5, 8.75, 10.0, 32.5, 55.0, 77.5, 100.0};
-//  std::vector<double> ww={0.4, 1.6, 0.8, 1.6, 2.4, 8.0, 4.0, 8.0, 4.5, 10.0, 5.0, 10.0, 47.5, 180.0, 90.0, 180.0, 45.0};
   std::vector<double> rr={0.2, 0.4, 0.6, 0.8, 1.0, 1.25, 1.5, 1.75, 2.0, 2.5, 3.0, 4.0, 5.0, 7.5, 10.0, 12.8125, 15.625, 18.4375, 21.25, 26.875, 32.5, 43.75, 55.0, 77.5, 100.0};
-  std::vector<double> ww={0.4, 1.6, 0.8, 1.6, 0.9, 2.0, 1.0, 2.0, 1.5, 4.0, 3.0, 8.0, 7.0, 20.0, 10.625, 22.5, 11.25, 22.5, 16.875, 45.0, 33.75, 90.0, 67.5, 180.0, 45.0};
 
   std::vector<double> pp={fParam[0], 0.0, fParam[2]};
   std::vector<double> pol;
-  double scale, rr2;
+  double scale, up{0.0}, low{-1.0};
   fTime.clear();
   fPol.clear();
 
@@ -105,6 +80,8 @@ int PLGKT_LF::DynamicLGKTLF()
     return 0;
   }
 
+  auto t_start = std::chrono::high_resolution_clock::now();
+  double sqrtTwoInv = 1.0/sqrt(2.0);
   for (unsigned int i=0; i<rr.size(); i++) {
     pp[1] = rr[i] * fParam[1];
     PGKT_LF gkt_lf(pp, fTmax);
@@ -116,16 +93,18 @@ int PLGKT_LF::DynamicLGKTLF()
       fTime = gkt_lf.GetTime();
       fPol.resize(fTime.size());
     }
-    rr2 = pow(rr[i],2.0);
-    scale = ww[i]*exp(-0.5/rr2)/rr2;
+    up = -std::erf(sqrtTwoInv/rr[i]);
+    scale = up - low;
+    low = up;
+
     pol.clear();
     pol = gkt_lf.GetPol();
     for (unsigned int j=0; j<fPol.size(); j++)
       fPol[j] += scale*pol[j];
   }
 
-  scale = pow(2.0/std::numbers::pi_v<double>, 0.5)/6.0;
-  std::transform(fPol.begin(), fPol.end(), fPol.begin(), [&scale](double el) { return el *= scale;});
+  auto t_end = std::chrono::high_resolution_clock::now();
+  std::cout << ">> time used: " << std::chrono::duration<double, std::milli>(t_end-t_start).count() << " ms." << std::endl;
 
   return 0;
 }