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

first LGKT_LF. Might need to slightly improve Simpson.

Browse Source
This commit is contained in:
suter_a 2025-05-12 16:36:37 +02:00
parent 819d209863
commit 69ab03d5ae
6 changed files with 180 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
src/external/LF_GL/CMakeLists.txt vendored
View File

@ -58,6 +58,7 @@ message("")
add_executable(lf_gl
main.cpp
PGKT_LF.cpp
PLGKT_LF.cpp
)
set_property(TARGET lf_gl PROPERTY CXX_STANDARD 20)
set(CMAKE_CXX_STANDARD_REQUIRED ON)

5
src/external/LF_GL/PGKT_LF.cpp vendored
View File

@ -1,6 +1,6 @@
/***************************************************************************
PGKT_LF.h
PGKT_LF.cpp
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
@ -41,8 +41,7 @@
* @param tt time vector in (ns)
* @param pol polarization vector
*/
PGKT_LF::PGKT_LF(std::vector<double> &param, std::vector<double> &tt, std::vector<double> &pol) :
fParam(param), fTime(tt), fPol(pol)
PGKT_LF::PGKT_LF(std::vector<double> &param) : fParam(param)
{
fParam[0] *= fGammaMu;

2
src/external/LF_GL/PGKT_LF.h vendored
View File

@ -35,7 +35,7 @@
class PGKT_LF
{
public:
PGKT_LF(std::vector<double> &param, std::vector<double> &tt, std::vector<double> &pol);
PGKT_LF(std::vector<double> &param);
int DynamicGaussKTLF();
bool IsValid() { return fValid; }

109
src/external/LF_GL/PLGKT_LF.cpp vendored Normal file
View File

@ -0,0 +1,109 @@
/***************************************************************************
PLGKT_LF.cpp
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2025 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program; if not, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#include <iostream>
#include <cmath>
#include <algorithm>
#include "PLGKT_LF.h"
#include "PGKT_LF.h"
//-----------------------------------------------------------------------------
PLGKT_LF::PLGKT_LF(std::vector<double> &param) : fParam(param)
{
if (DynamicLGKTLF() == 0)
fValid = true;
}
//-----------------------------------------------------------------------------
/**
* <p>The weight density \f$\rho_{\Delta_{\rm L}}(\Delta_{\rm G}) dG_{\rm G} =
* \sqrt{\frac{2}{\pi}} \frac{1}{r^2} \exp\left(-\frac{1}{2 r^2}\right) dr\f$, with
* \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>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(r) = P_z(t, \Delta_{\rm G} = r\cdot\Delta_{\rm L}, B, \nu)
*
* @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> pp={fParam[0], 0.0, fParam[2]};
std::vector<double> pol;
double scale, rr2;
fTime.clear();
fPol.clear();
for (unsigned int i=0; i<rr.size(); i++) {
pp[1] = rr[i] * fParam[1];
PGKT_LF gkt_lf(pp);
if (!gkt_lf.IsValid()) {
std::cout << "**ERROR** in Gaussian LF calculation" << std::endl;
return 2;
}
if (fTime.size()==0) {
fTime = gkt_lf.GetTime();
fPol.resize(fTime.size());
}
rr2 = pow(rr[i],2.0);
scale = ww[i]*exp(-0.5/rr2)/rr2;
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;});
return 0;
}

52
src/external/LF_GL/PLGKT_LF.h vendored Normal file
View File

@ -0,0 +1,52 @@
/***************************************************************************
PLGKT_LF.h
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2025 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program; if not, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifndef _PLGKT_LF_
#define _PLGKT_LF_
#include <vector>
class PLGKT_LF {
public:
PLGKT_LF(std::vector<double> &param);
int DynamicLGKTLF();
bool IsValid() { return fValid; }
std::vector<double> GetTime() { return fTime; }
std::vector<double> GetPol() { return fPol; }
private:
bool fValid{false};
std::vector<double> fParam;
std::vector<double> fTime;
std::vector<double> fPol;
};
#endif // _PLGKT_LF_

19
src/external/LF_GL/main.cpp vendored
View File

@ -7,6 +7,7 @@
#include <numbers>
#include "PGKT_LF.h"
#include "PLGKT_LF.h"
// ----------------------------------------------------------------------------
void lf_gl_syntax()
@ -23,11 +24,15 @@ void lf_gl_syntax()
}
// ----------------------------------------------------------------------------
int lf_gl_write(const std::string fln, const std::vector<double> &param,
int lf_gl_write(const std::string fln, const bool gaussian_only, const std::vector<double> &param,
const std::vector<double> &tt, const std::vector<double> &pol)
{
std::ofstream fout(fln, std::ofstream::out);
if (gaussian_only)
fout << "# Gaussian only" << std::endl;
else
fout << "# Gaussian/Lorentzian" << std::endl;
fout << "# field=" << param[0] << " (G), width=" << param[1] << " (1/us), hopp=" << param[2] << " (1/us)" << std::endl;
fout << "# t (us), pol" << std::endl;
for (unsigned int i=0; i<tt.size(); i++) {
@ -115,7 +120,7 @@ int main(int argc, char *argv[])
std::vector<double> tt;
std::vector<double> pol;
if (gaussian_only) {
PGKT_LF gkt_lf(param, tt, pol);
PGKT_LF gkt_lf(param);
if (!gkt_lf.IsValid()) {
std::cout << "**ERROR** in Gaussian LF calculation" << std::endl;
return 2;
@ -123,10 +128,16 @@ int main(int argc, char *argv[])
tt = gkt_lf.GetTime();
pol = gkt_lf.GetPol();
} else {
std::cout << "not yet implemented." << std::endl;
PLGKT_LF lgkt_lf(param);
if (!lgkt_lf.IsValid()) {
std::cout << "**ERROR** in Gaussian/Lorentzian LF calculation" << std::endl;
return 2;
}
tt = lgkt_lf.GetTime();
pol = lgkt_lf.GetPol();
}
lf_gl_write(flnOut, param, tt, pol);
lf_gl_write(flnOut, gaussian_only, param, tt, pol);
return 0;
}