LMU/musrfit
5
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases 4 Activity

changed the integral approximation approach for Gaussian/Lorentzian.

Browse Source
This commit is contained in:
suter_a
2025-05-28 12:46:37 +02:00
parent 8135035899
commit 60aeab2239
1 changed files with 11 additions and 32 deletions
Download Patch File Download Diff File Expand all files Collapse all files

43
src/external/LF_GL/PLGKT_LF.cpp vendored
View File

@ -31,6 +31,7 @@
#include <cmath> #include <cmath>
#include <algorithm> #include <algorithm>
#include <numbers> #include <numbers>
#include <chrono>
#include "PLGKT_LF.h" #include "PLGKT_LF.h"
#include "PGKT_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$ * \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) * (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: * <p>Integration method is described in the docu directory.
* \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>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 * @return 0 on success, >0 otherwise
*/ */
int PLGKT_LF::DynamicLGKTLF() 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> 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> pp={fParam[0], 0.0, fParam[2]};
std::vector<double> pol; std::vector<double> pol;
double scale, rr2; double scale, up{0.0}, low{-1.0};
fTime.clear(); fTime.clear();
fPol.clear(); fPol.clear();
@ -105,6 +80,8 @@ int PLGKT_LF::DynamicLGKTLF()
return 0; 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++) { for (unsigned int i=0; i<rr.size(); i++) {
pp[1] = rr[i] * fParam[1]; pp[1] = rr[i] * fParam[1];
PGKT_LF gkt_lf(pp, fTmax); PGKT_LF gkt_lf(pp, fTmax);
@ -116,16 +93,18 @@ int PLGKT_LF::DynamicLGKTLF()
fTime = gkt_lf.GetTime(); fTime = gkt_lf.GetTime();
fPol.resize(fTime.size()); fPol.resize(fTime.size());
} }
rr2 = pow(rr[i],2.0); up = -std::erf(sqrtTwoInv/rr[i]);
scale = ww[i]*exp(-0.5/rr2)/rr2; scale = up - low;
low = up;
pol.clear(); pol.clear();
pol = gkt_lf.GetPol(); pol = gkt_lf.GetPol();
for (unsigned int j=0; j<fPol.size(); j++) for (unsigned int j=0; j<fPol.size(); j++)
fPol[j] += scale*pol[j]; fPol[j] += scale*pol[j];
} }
scale = pow(2.0/std::numbers::pi_v<double>, 0.5)/6.0; auto t_end = std::chrono::high_resolution_clock::now();
std::transform(fPol.begin(), fPol.end(), fPol.begin(), [&scale](double el) { return el *= scale;}); std::cout << ">> time used: " << std::chrono::duration<double, std::milli>(t_end-t_start).count() << " ms." << std::endl;
return 0; return 0;
} }