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

- Extended the usage of xml-startup files in some plugin libraries and renamed the file itself.

Browse Source 
- Added some documentation to libLFRelaxation which in the present stage is mainly used for testing purposes.
This commit is contained in:
Bastian M. Wojek
2011-03-16 16:43:13 +00:00
parent 233ca10dff
commit 032ce4a1e4
18 changed files with 964 additions and 436 deletions
Download Patch File Download Diff File Expand all files Collapse all files

23
src/external/libLFRelaxation/Makefile.am vendored
View File

@ -1,31 +1,40 @@
## Process this file with automake to create Makefile.in
startup_path = ../BMWStartupHandler
integrator_path = ../BMWIntegrator
h_sources = \
TLFRelaxation.h \
../BMWIntegrator/BMWIntegrator.h
$(startup_path)/BMWStartupHandler.h \
$(integrator_path)/BMWIntegrator.h
h_linkdef = \
TLFRelaxationLinkDef.h
TLFRelaxationLinkDef.h \
$(startup_path)/BMWStartupHandlerLinkDef.h
dict_h_sources = \
TLFRelaxationDict.h
TLFRelaxationDict.h \
$(startup_path)/BMWStartupHandlerDict.h
cpp_sources = \
TLFRelaxation.cpp \
../BMWIntegrator/BMWIntegrator.cpp
$(startup_path)/BMWStartupHandler.cpp \
$(integrator_path)/BMWIntegrator.cpp
dict_cpp_sources = \
TLFRelaxationDict.cpp
TLFRelaxationDict.cpp \
$(startup_path)/BMWStartupHandlerDict.cpp
include_HEADERS = $(h_sources)
noinst_HEADERS = $(h_linkdef) $(dict_h_sources)
INCLUDES = -I$(top_srcdir)/src/include $(PMUSR_CFLAGS) $(FFTW3_CFLAGS) $(GSL_CFLAGS) $(ROOT_CFLAGS) $(CUBA_CFLAGS)
INCLUDES = -I$(top_srcdir)/src/include -I$(startup_path) $(PMUSR_CFLAGS) $(FFTW3_CFLAGS) $(GSL_CFLAGS) $(ROOT_CFLAGS) $(CUBA_CFLAGS)
AM_CXXFLAGS = $(LOCAL_LIB_CXXFLAGS)
BUILT_SOURCES = $(dict_cpp_sources) $(dict_h_sources)
AM_LDFLAGS = $(LOCAL_LIB_LDFLAGS) -L@ROOTLIBDIR@
CLEANFILES = *Dict.cpp *Dict.h *~ core
CLEANFILES = $(startup_path)/*Dict.cpp $(startup_path)/*Dict.h $(startup_path)/*~ \
$(integrator_path)/*~ *Dict.cpp *Dict.h *~ core
%Dict.cpp %Dict.h: %.h %LinkDef.h
@ROOTCINT@ -v -f $*Dict.cpp -c -p $(INCLUDES) $^

512
src/external/libLFRelaxation/TLFRelaxation.cpp vendored
View File

@ -5,7 +5,7 @@
Author: Bastian M. Wojek
e-mail: bastian.wojek@psi.ch
2008/12/04
$Id$
***************************************************************************/
@ -29,12 +29,22 @@
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <cassert>
#include <sys/time.h>
#include <iostream>
#include <cmath>
#ifdef HAVE_GOMP
#include <omp.h>
#endif
#include <TSAXParser.h>
#include "../BMWIntegrator/BMWIntegrator.h"
#include "BMWStartupHandler.h"
#include "TLFRelaxation.h"
using namespace std;
@ -44,23 +54,46 @@ using namespace std;
#define Nsteps 100
ClassImp(TLFStatGssKT)
ClassImp(TLFStatLorKT)
ClassImp(TLFStatExpKT)
ClassImp(TLFDynGssKT)
ClassImp(TLFDynLorKT)
ClassImp(TLFDynExpKT)
ClassImp(TLFDynSG)
ClassImp(TLFSGInterpolation)
// LF Static Gaussian KT
TLFStatGssKT::TLFStatGssKT() {
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
* <p>Constructor: Initialize variables and allocate memory for the integrator.
*/
TLFStatGssKT::TLFStatGssKT() : fCalcNeeded(true), fLastTime(0.) {
fIntSinGss = new TIntSinGss();
fPar.resize(2, 0.);
}
//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
/**
* <p>Destructor: Free memory of the integrator, clean up some vectors.
*/
TLFStatGssKT::~TLFStatGssKT() {
delete fIntSinGss;
fIntSinGss = 0;
delete fIntSinGss;
fIntSinGss = 0;
fIntValues.clear();
fPar.clear();
}
//--------------------------------------------------------------------------
// operator()
//--------------------------------------------------------------------------
/**
* <p>Method returning the function value at a given time for a given set of parameters.
*
* \param t time \htmlonly (&#956;s) \endhtmlonly \latexonly ($\mu\mathrm{s}$) \endlatexonly
* \param par parameters [\htmlonly &#957;<sub>L</sub>=<i>B</i>&#947;<sub>&#956;</sub>/2&#960; (MHz), &#963; (&#956;s<sup>-1</sup>)\endhtmlonly \latexonly $\nu_{\mathrm{L}}=B\gamma_{\mu}/2\pi~(\mathrm{MHz})$, $\sigma~(\mu\mathrm{s}^{-1})$ \endlatexonly]
*/
double TLFStatGssKT::operator()(double t, const vector<double> &par) const {
assert(par.size() == 2); // two parameters nu=gbar*B,sigma
@ -69,32 +102,77 @@ double TLFStatGssKT::operator()(double t, const vector<double> &par) const {
return 1.0;
double sigsq(par[1]*par[1]);
double sigsqtt(sigsq*t*t);
if(TMath::Abs(par[0])<0.00135538817){
return (0.33333333333333333333+0.66666666666666666667*(1.0-sigsq*t*t)*TMath::Exp(-0.5*sigsq*t*t));
return (0.33333333333333333333+0.66666666666666666667*(1.0-sigsqtt)*TMath::Exp(-0.5*sigsqtt));
}
fIntSinGss->SetParameters(par);
for (unsigned int i(0); i<2; ++i) {
if( fPar[i]-par[i] ) {
fPar[i] = par[i];
fCalcNeeded = true;
}
}
if (fCalcNeeded) {
fIntValues.clear();
fIntValues[0.] = 0.;
fLastTime = 0.;
} else {
if (t > fIntValues.rbegin()->first) {
fCalcNeeded = true;
fLastTime = fIntValues.rbegin()->first;
}
}
double omegaL(TMath::TwoPi()*par[0]);
double coeff1(2.0*sigsq/(omegaL*omegaL));
double coeff2(2.0*sigsq*sigsq/(omegaL*omegaL*omegaL));
double coeff2(coeff1*sigsq/omegaL);
return (1.0-(coeff1*(1.0-TMath::Exp(-0.5*sigsq*t*t)*TMath::Cos(omegaL*t)))+(coeff2*fIntSinGss->IntegrateFunc(0.,t)));
if (fCalcNeeded) {
fIntSinGss->SetParameters(par);
fIntValues[t] = fIntValues[fLastTime] + fIntSinGss->IntegrateFunc(fLastTime,t);
fCalcNeeded = false;
}
return (1.0-(coeff1*(1.0-TMath::Exp(-0.5*sigsqtt)*TMath::Cos(omegaL*t)))+(coeff2*fIntValues[t]));
}
// LF Static Lorentzian KT
TLFStatLorKT::TLFStatLorKT() {
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
* <p>Constructor: Initialize variables and allocate memory for the integrator.
*/
TLFStatExpKT::TLFStatExpKT() : fCalcNeeded(true), fLastTime(0.) {
fIntBesselJ0Exp = new TIntBesselJ0Exp();
fPar.resize(2, 0.);
}
TLFStatLorKT::~TLFStatLorKT() {
delete fIntBesselJ0Exp;
fIntBesselJ0Exp = 0;
//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
/**
* <p>Destructor: Free memory of the integrator, clean up some vectors.
*/
TLFStatExpKT::~TLFStatExpKT() {
delete fIntBesselJ0Exp;
fIntBesselJ0Exp = 0;
fIntValues.clear();
fPar.clear();
}
double TLFStatLorKT::operator()(double t, const vector<double> &par) const {
//--------------------------------------------------------------------------
// operator()
//--------------------------------------------------------------------------
/**
* <p>Method returning the function value at a given time for a given set of parameters.
*
* \param t time \htmlonly (&#956;s) \endhtmlonly \latexonly ($\mu\mathrm{s}$) \endlatexonly
* \param par parameters [\htmlonly &#957;<sub>L</sub>=<i>B</i>&#947;<sub>&#956;</sub>/2&#960; (MHz), <i>a</i> (&#956;s<sup>-1</sup>)\endhtmlonly \latexonly $\nu_{\mathrm{L}}=B\gamma_{\mu}/2\pi~(\mathrm{MHz})$, $a~(\mu\mathrm{s}^{-1})$ \endlatexonly]
*/
double TLFStatExpKT::operator()(double t, const vector<double> &par) const {
assert(par.size() == 2); // two parameters nu=gbar*B,rate
@ -105,7 +183,23 @@ double TLFStatLorKT::operator()(double t, const vector<double> &par) const {
return (0.33333333333333333333+0.66666666666666666667*(1.0-par[1]*t)*TMath::Exp(-par[1]*t));
}
fIntBesselJ0Exp->SetParameters(par);
for (unsigned int i(0); i<2; ++i) {
if( fPar[i]-par[i] ) {
fPar[i] = par[i];
fCalcNeeded = true;
}
}
if (fCalcNeeded) {
fIntValues.clear();
fIntValues[0.] = 0.;
fLastTime = 0.;
} else {
if (t > fIntValues.rbegin()->first) {
fCalcNeeded = true;
fLastTime = fIntValues.rbegin()->first;
}
}
double omegaL(TMath::TwoPi()*par[0]);
double coeff1(par[1]/omegaL);
@ -114,7 +208,7 @@ double TLFStatLorKT::operator()(double t, const vector<double> &par) const {
double w0t(omegaL*t);
double j1, j0;
if (fabs(w0t) < 0.001) { // check zero time limits of the spherical bessel functions j0(x) and j1(x)
if (fabs(w0t) < 1.e-5) { // check zero time limits of the spherical bessel functions j0(x) and j1(x)
j0 = 1.0;
j1 = 0.0;
} else {
@ -122,50 +216,113 @@ double TLFStatLorKT::operator()(double t, const vector<double> &par) const {
j1 = (TMath::Sin(w0t)-w0t*TMath::Cos(w0t))/(w0t*w0t);
}
return (1.0-(coeff1*TMath::Exp(-par[1]*t)*j1)-(coeff2*(j0*TMath::Exp(-par[1]*t)-1.0))-coeff3*fIntBesselJ0Exp->IntegrateFunc(0.,t));
if (fCalcNeeded) {
fIntBesselJ0Exp->SetParameters(par);
fIntValues[t] = fIntValues[fLastTime] + fIntBesselJ0Exp->IntegrateFunc(fLastTime,t);
fCalcNeeded = false;
}
return (1.0-(coeff1*TMath::Exp(-par[1]*t)*j1)-(coeff2*(j0*TMath::Exp(-par[1]*t)-1.0))-coeff3*fIntValues[t]);
}
// LF Dynamic Gaussian KT
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
* <p>Constructor
* - read XML startup file
* - initialize variables
* - read FFTW3 wisdom if desired
* - allocate memory for the Laplace transform
* - create FFTW3 plans for the Laplace transform
*/
TLFDynGssKT::TLFDynGssKT() : fCalcNeeded(true), fFirstCall(true), fCounter(0) {
#ifdef HAVE_LIBFFTW3F_THREADS
int init_threads(fftwf_init_threads());
if (!init_threads)
cout << "TLFDynGssKT::TLFDynGssKT: Couldn't initialize multiple FFTW-threads ..." << endl;
else {
#ifdef HAVE_GOMP
fftwf_plan_with_nthreads(omp_get_num_procs());
#else
fftwf_plan_with_nthreads(2);
#endif /* HAVE_GOMP */
}
#endif /* HAVE_LIBFFTW3F_THREADS */
// read startup file
string startup_path_name("BMW_startup.xml");
TSAXParser *saxParser = new TSAXParser();
BMWStartupHandler *startupHandler = new BMWStartupHandler();
saxParser->ConnectToHandler("BMWStartupHandler", startupHandler);
int status (saxParser->ParseFile(startup_path_name.c_str()));
// check for parse errors
if (status) { // error
cerr << endl << "**ERROR** reading/parsing " << startup_path_name << " failed." \
<< endl << "**ERROR** Please make sure that the file exists in the local directory and it is set up correctly!" \
<< endl;
assert(false);
}
fNSteps = startupHandler->GetNStepsLF();
fDt = startupHandler->GetDeltat();
fWisdom = startupHandler->GetWisdomFileFloat();
TLFDynGssKT::TLFDynGssKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("WordsOfWisdom.dat"), fNSteps(524288), fDt(0.00004), fCounter(0) {
// Calculate d_omega and C for given NFFT and dt
fDw = PI/fNSteps/fDt;
fDw = PI/(fNSteps*fDt);
fC = 2.0*TMath::Log(double(fNSteps))/(double(fNSteps-1)*fDt);
// Load FFTW Wisdom
fUseWisdom = true;
int wisdomLoaded(0);
FILE *wordsOfWisdomR;
wordsOfWisdomR = fopen(fWisdom.c_str(), "r");
if (wordsOfWisdomR == NULL) {
cout << "TLFDynGssKT::TLFDynGssKT: Couldn't open wisdom file ..." << endl;
fUseWisdom = false;
} else {
wisdomLoaded = fftwf_import_wisdom_from_file(wordsOfWisdomR);
fclose(wordsOfWisdomR);
}
if (!wisdomLoaded) {
cout << "TLFDynGssKT::TLFDynGssKT: No wisdom is imported..." << endl;
fUseWisdom = false;
}
// END of WisdomLoading
// allocating memory for the FFtransform pairs and create the FFT plans
fFFTtime = (float *)malloc(sizeof(float) * fNSteps);
fFFTfreq = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * (fNSteps/2+1));
fFFTplanFORW = fftwf_plan_dft_r2c_1d(fNSteps, fFFTtime, fFFTfreq, FFTW_ESTIMATE);
fFFTplanBACK = fftwf_plan_dft_c2r_1d(fNSteps, fFFTfreq, fFFTtime, FFTW_ESTIMATE);
if (fUseWisdom) {
fFFTplanFORW = fftwf_plan_dft_r2c_1d(fNSteps, fFFTtime, fFFTfreq, FFTW_EXHAUSTIVE);
fFFTplanBACK = fftwf_plan_dft_c2r_1d(fNSteps, fFFTfreq, fFFTtime, FFTW_EXHAUSTIVE);
} else {
fFFTplanFORW = fftwf_plan_dft_r2c_1d(fNSteps, fFFTtime, fFFTfreq, FFTW_ESTIMATE);
fFFTplanBACK = fftwf_plan_dft_c2r_1d(fNSteps, fFFTfreq, fFFTtime, FFTW_ESTIMATE);
}
}
//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
/**
* <p>Destructor
* - write FFTW3 wisdom if desired
* - free memory used for the Laplace transform
*/
TLFDynGssKT::~TLFDynGssKT() {
// export FFTW Wisdom so it has not to be checked for the FFT-plan next time
FILE *wordsOfWisdomW;
wordsOfWisdomW = fopen(fWisdom.c_str(), "w");
if (wordsOfWisdomW == NULL) {
cout << "TLFDynGssKT::~TLFDynGssKT: Could not open file ... No wisdom is exported..." << endl;
} else {
fftwf_export_wisdom_to_file(wordsOfWisdomW);
fclose(wordsOfWisdomW);
if (fUseWisdom) {
FILE *wordsOfWisdomW;
wordsOfWisdomW = fopen(fWisdom.c_str(), "w");
if (wordsOfWisdomW == NULL) {
cout << "TLFDynGssKT::~TLFDynGssKT(): Could not open file ... No wisdom is exported..." << endl;
} else {
fftwf_export_wisdom_to_file(wordsOfWisdomW);
fclose(wordsOfWisdomW);
}
}
// END of Wisdom Export
@ -177,9 +334,18 @@ TLFDynGssKT::~TLFDynGssKT() {
cout << "TLFDynGssKT::~TLFDynGssKT(): " << fCounter << " full FFT cycles needed..." << endl;
}
//--------------------------------------------------------------------------
// operator()
//--------------------------------------------------------------------------
/**
* <p>Method returning the function value at a given time for a given set of parameters.
*
* \param t time \htmlonly (&#956;s) \endhtmlonly \latexonly ($\mu\mathrm{s}$) \endlatexonly
* \param par parameters [\htmlonly &#957;<sub>L</sub>=<i>B</i>&#947;<sub>&#956;</sub>/2&#960; (MHz), &#963; (&#956;s<sup>-1</sup>), &#957; (MHz)\endhtmlonly \latexonly $\nu_{\mathrm{L}}=B\gamma_{\mu}/2\pi~(\mathrm{MHz})$, $\sigma~(\mu\mathrm{s}^{-1})$, $\nu~(\mathrm{MHz})$ \endlatexonly]
*/
double TLFDynGssKT::operator()(double t, const vector<double> &par) const {
assert(par.size() == 3 || par.size() == 4); // three parameters nuL=gbar*B,sigma,fluct.rate nu (fourth parameter: t[us] for SG-function)
assert(par.size() == 3); // three parameters nuL=gbar*B,sigma,fluct.rate nu (fourth parameter: t[us] for SG-function)
if(t<0.0)
return 1.0;
@ -192,7 +358,7 @@ double TLFDynGssKT::operator()(double t, const vector<double> &par) const {
fFirstCall=false;
}
for (unsigned int i(0); i<3; i++) {
for (unsigned int i(0); i<3; ++i) {
if( fPar[i]-par[i] ) {
fPar[i] = par[i];
fCalcNeeded=true;
@ -228,12 +394,16 @@ double TLFDynGssKT::operator()(double t, const vector<double> &par) const {
struct timeval tv_start, tv_stop;
gettimeofday(&tv_start, 0);
*/
double tt(0.),sigsqtsq(0.);
double tt(0.), sigsqtsq(0.);
int i;
if(fabs(par[0])<0.00135538817) {
double mcplusnu(-(fC+par[2]));
for(unsigned int i(0); i<fNSteps; i++) {
tt=double(i)*fDt;
#ifdef HAVE_GOMP
#pragma omp parallel for default(shared) private(i,tt) schedule(dynamic)
#endif
for(i = 0; i < static_cast<int>(fNSteps); ++i) {
tt=static_cast<double>(i)*fDt;
sigsqtsq=sigsq*tt*tt;
fFFTtime[i]=(0.33333333333333333333+0.66666666666666666667*(1.0-sigsqtsq)*exp(-0.5*sigsqtsq))*exp(mcplusnu*tt)*fDt;
}
@ -244,15 +414,18 @@ double TLFDynGssKT::operator()(double t, const vector<double> &par) const {
fFFTtime[0] = fDt;
for(unsigned int i(1); i<fNSteps; i++) {
tt=(double(i)-0.5)*fDt;
#ifdef HAVE_GOMP
#pragma omp parallel for default(shared) private(i,tt) schedule(dynamic)
#endif
for(i = 1; i < static_cast<int>(fNSteps); ++i) {
tt=(static_cast<double>(i)-0.5)*fDt;
fFFTtime[i]=(sin(omegaL*tt) * exp(-0.5*sigsq*tt*tt))*fDt;
}
double totoIntegrale(0.);
for(unsigned int i(1); i<fNSteps; i++) {
tt=double(i)*fDt;
for(i = 1; i < static_cast<int>(fNSteps); ++i) {
tt=static_cast<double>(i)*fDt;
totoIntegrale+=fFFTtime[i];
fFFTtime[i]=(1.0-(coeff1*(1.0-exp(-0.5*sigsq*tt*tt)*cos(omegaL*tt)))+(coeff2*totoIntegrale))*exp(mcplusnu*tt)*fDt;
}
@ -269,7 +442,10 @@ double TLFDynGssKT::operator()(double t, const vector<double> &par) const {
double denom(1.0), imagsq(0.0), oneMINrealnu(1.0);
for (unsigned int i(0); i<fNSteps/2+1; i++) {
#ifdef HAVE_GOMP
#pragma omp parallel for default(shared) private(i, imagsq, oneMINrealnu, denom) schedule(dynamic)
#endif
for (i = 0; i < static_cast<int>(fNSteps)/2+1; ++i) {
imagsq=fFFTfreq[i][1]*fFFTfreq[i][1];
oneMINrealnu=1.0-(par[2]*fFFTfreq[i][0]);
denom=oneMINrealnu*oneMINrealnu + (nusq*imagsq);
@ -285,7 +461,7 @@ double TLFDynGssKT::operator()(double t, const vector<double> &par) const {
// fFFTtime[i]=(fDw*TMath::Exp(fC*i*fDt)/TMath::Pi()*fFFTtime[i]);
// }
fCalcNeeded=false;
fCounter++;
++fCounter;
/*
gettimeofday(&tv_stop, 0);
t2 = (tv_stop.tv_sec - tv_start.tv_sec)*1000.0 + (tv_stop.tv_usec - tv_start.tv_usec)/1000.0;
@ -293,17 +469,27 @@ double TLFDynGssKT::operator()(double t, const vector<double> &par) const {
*/
}
// return fFFTtime[int(t/fDt)];
return fDw*exp(fC*t)/PI*fFFTtime[int(t/fDt)];
return fDw*exp(fC*t)/PI*fFFTtime[static_cast<int>(t/fDt)];
}
// LF Dynamic Spin Glass Relaxation
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
* <p>Constructor: create dynamic-Gaussian-depolarization objects for a subset of static widths and initialize the parameter vector.
*/
TLFDynSG::TLFDynSG() {
for(unsigned int i(0); i<Nsteps; ++i)
fLFDynGssIntegral.push_back(new TLFDynGssKT());
fPar.resize(3, 0.0);
}
//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
/**
* <p>Destructor: free the memory of the dynamic-Gaussian-depolarization objects and cleanup the parameter vector.
*/
TLFDynSG::~TLFDynSG() {
for(unsigned int i(0); i<Nsteps; ++i)
delete fLFDynGssIntegral[i];
@ -311,6 +497,15 @@ TLFDynSG::~TLFDynSG() {
fPar.clear();
}
//--------------------------------------------------------------------------
// operator()
//--------------------------------------------------------------------------
/**
* <p>Method returning the function value at a given time for a given set of parameters.
*
* \param t time \htmlonly (&#956;s) \endhtmlonly \latexonly ($\mu\mathrm{s}$) \endlatexonly
* \param par parameters [\htmlonly &#957;<sub>L</sub>=<i>B</i>&#947;<sub>&#956;</sub>/2&#960; (MHz), <i>a</i> (&#956;s<sup>-1</sup>), &#957; (MHz)\endhtmlonly \latexonly $\nu_{\mathrm{L}}=B\gamma_{\mu}/2\pi~(\mathrm{MHz})$, $a~(\mu\mathrm{s}^{-1})$, $\nu~(\mathrm{MHz})$ \endlatexonly]
*/
double TLFDynSG::operator()(double t, const vector<double> &par) const {
assert(par.size() == 3); // three parameters nuL=gbar*B, a ,fluct.rate nu
@ -321,19 +516,21 @@ double TLFDynSG::operator()(double t, const vector<double> &par) const {
if(t>20.0)
return 0.0;
const double upperCutOff(50.0); // times a
// Parameters for the integration
for (unsigned int i(0); i<3; ++i)
fPar[i] = par[i];
double sigma_step((40.0-0.1)*par[1]/static_cast<double>(Nsteps));
// at the moment, integrate between 0.1a and 40a (a "real upper limit" would be 10000a)
double sigma_step(upperCutOff*par[1]/static_cast<double>(Nsteps));
// at the moment, integrate between 0 and 50a (a "real upper limit" would be 10000a)
double integral(0.0);
double a_sigmaSq(0.0);
double rho(0.0);
double rhoNormalizer(0.0); // is [1/]sqrt(2./PI) when the integration is done between 0 and infinity
for (unsigned int i(0); i<Nsteps; ++i) {
fPar[1] = 0.1*par[1] + (static_cast<double>(i) + 0.5)*sigma_step;
fPar[1] = (static_cast<double>(i) + 0.5)*sigma_step;
a_sigmaSq = (par[1]/(fPar[1]*fPar[1]));
rho = a_sigmaSq*exp(-0.5*par[1]*a_sigmaSq);
rhoNormalizer += rho;
@ -347,47 +544,104 @@ double TLFDynSG::operator()(double t, const vector<double> &par) const {
}
// LF Dynamic Lorentz KT
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
* <p>Constructor
* - read XML startup file
* - initialize variables
* - read FFTW3 wisdom if desired
* - allocate memory for the Laplace transform
* - create FFTW3 plans for the Laplace transform
*/
TLFDynExpKT::TLFDynExpKT() : fCalcNeeded(true), fFirstCall(true), fCounter(0), fL1(0.0), fL2(0.0) {
#ifdef HAVE_LIBFFTW3F_THREADS
int init_threads(fftwf_init_threads());
if (!init_threads)
cout << "TLFDynExpKT::TLFDynExpKT: Couldn't initialize multiple FFTW-threads ..." << endl;
else {
#ifdef HAVE_GOMP
fftwf_plan_with_nthreads(omp_get_num_procs());
#else
fftwf_plan_with_nthreads(2);
#endif /* HAVE_GOMP */
}
#endif /* HAVE_LIBFFTW3F_THREADS */
// read startup file
string startup_path_name("BMW_startup.xml");
TSAXParser *saxParser = new TSAXParser();
BMWStartupHandler *startupHandler = new BMWStartupHandler();
saxParser->ConnectToHandler("BMWStartupHandler", startupHandler);
int status (saxParser->ParseFile(startup_path_name.c_str()));
// check for parse errors
if (status) { // error
cerr << endl << "**ERROR** reading/parsing " << startup_path_name << " failed." \
<< endl << "**ERROR** Please make sure that the file exists in the local directory and it is set up correctly!" \
<< endl;
assert(false);
}
fNSteps = startupHandler->GetNStepsLF();
fDt = startupHandler->GetDeltat();
fWisdom = startupHandler->GetWisdomFileFloat();
TLFDynLorKT::TLFDynLorKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("WordsOfWisdom.dat"), fNSteps(524288), fDt(0.000040), fCounter(0), fL1(0.0), fL2(0.0) {
// Calculate d_omega and C for given NFFT and dt
fDw = TMath::Pi()/fNSteps/fDt;
fDw = PI/(fNSteps*fDt);
fC = 2.0*TMath::Log(double(fNSteps))/(double(fNSteps-1)*fDt);
// Load FFTW Wisdom
fUseWisdom = true;
int wisdomLoaded(0);
FILE *wordsOfWisdomR;
wordsOfWisdomR = fopen(fWisdom.c_str(), "r");
if (wordsOfWisdomR == NULL) {
cout << "TLFDynLorKT::TLFDynLorKT: Couldn't open wisdom file ..." << endl;
fUseWisdom = false;
} else {
wisdomLoaded = fftwf_import_wisdom_from_file(wordsOfWisdomR);
fclose(wordsOfWisdomR);
}
if (!wisdomLoaded) {
cout << "TLFDynLorKT::TLFDynLorKT: No wisdom is imported..." << endl;
fUseWisdom = false;
}
// END of WisdomLoading
// allocating memory for the FFtransform pairs and create the FFT plans
fFFTtime = (float *)malloc(sizeof(float) * fNSteps);
fFFTfreq = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * (fNSteps/2+1));
fFFTplanFORW = fftwf_plan_dft_r2c_1d(fNSteps, fFFTtime, fFFTfreq, FFTW_ESTIMATE);
fFFTplanBACK = fftwf_plan_dft_c2r_1d(fNSteps, fFFTfreq, fFFTtime, FFTW_ESTIMATE);
if (fUseWisdom) {
fFFTplanFORW = fftwf_plan_dft_r2c_1d(fNSteps, fFFTtime, fFFTfreq, FFTW_EXHAUSTIVE);
fFFTplanBACK = fftwf_plan_dft_c2r_1d(fNSteps, fFFTfreq, fFFTtime, FFTW_EXHAUSTIVE);
} else {
fFFTplanFORW = fftwf_plan_dft_r2c_1d(fNSteps, fFFTtime, fFFTfreq, FFTW_ESTIMATE);
fFFTplanBACK = fftwf_plan_dft_c2r_1d(fNSteps, fFFTfreq, fFFTtime, FFTW_ESTIMATE);
}
}
TLFDynLorKT::~TLFDynLorKT() {
//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
/**
* <p>Destructor
* - write FFTW3 wisdom if desired
* - free memory used for the Laplace transform
*/
TLFDynExpKT::~TLFDynExpKT() {
// export FFTW Wisdom so it has not to be checked for the FFT-plan next time
FILE *wordsOfWisdomW;
wordsOfWisdomW = fopen(fWisdom.c_str(), "w");
if (wordsOfWisdomW == NULL) {
cout << "TLFDynLorKT::~TLFDynLorKT: Could not open file ... No wisdom is exported..." << endl;
} else {
fftwf_export_wisdom_to_file(wordsOfWisdomW);
fclose(wordsOfWisdomW);
if (fUseWisdom) {
FILE *wordsOfWisdomW;
wordsOfWisdomW = fopen(fWisdom.c_str(), "w");
if (wordsOfWisdomW == NULL) {
cout << "TLFDynExpKT::~TLFDynExpKT(): Could not open file ... No wisdom is exported..." << endl;
} else {
fftwf_export_wisdom_to_file(wordsOfWisdomW);
fclose(wordsOfWisdomW);
}
}
// END of Wisdom Export
@ -396,10 +650,19 @@ TLFDynLorKT::~TLFDynLorKT() {
fftwf_destroy_plan(fFFTplanBACK);
free(fFFTtime);
fftwf_free(fFFTfreq);
cout << "TLFDynLorKT::~TLFDynLorKT(): " << fCounter << " full FFT cyles needed..." << endl;
cout << "TLFDynExpKT::~TLFDynExpKT(): " << fCounter << " full FFT cyles needed..." << endl;
}
double TLFDynLorKT::operator()(double t, const vector<double> &par) const {
//--------------------------------------------------------------------------
// operator()
//--------------------------------------------------------------------------
/**
* <p>Method returning the function value at a given time for a given set of parameters.
*
* \param t time \htmlonly (&#956;s) \endhtmlonly \latexonly ($\mu\mathrm{s}$) \endlatexonly
* \param par parameters [\htmlonly &#957;<sub>L</sub>=<i>B</i>&#947;<sub>&#956;</sub>/2&#960; (MHz), <i>a</i> (&#956;s<sup>-1</sup>), &#957; (MHz)\endhtmlonly \latexonly $\nu_{\mathrm{L}}=B\gamma_{\mu}/2\pi~(\mathrm{MHz})$, $a~(\mu\mathrm{s}^{-1})$, $\nu~(\mathrm{MHz})$ \endlatexonly]
*/
double TLFDynExpKT::operator()(double t, const vector<double> &par) const {
assert(par.size() == 3); // three parameters nuL=gbar*B,sigma,fluct.rate nu
@ -445,7 +708,7 @@ double TLFDynLorKT::operator()(double t, const vector<double> &par) const {
double nuN[4];
w0N[0] = omegaL;
nuN[0] = nu;
for (unsigned int i=1; i<4; i++) {
for (unsigned int i=1; i<4; ++i) {
w0N[i] = omegaL * w0N[i-1];
nuN[i] = nu * nuN[i-1];
}
@ -458,7 +721,7 @@ double TLFDynLorKT::operator()(double t, const vector<double> &par) const {
double w0t(omegaL*t);
double j1, j0;
if (fabs(w0t) < 0.001) { // check zero time limits of the spherical bessel functions j0(x) and j1(x)
if (fabs(w0t) < 1.e-5) { // check zero time limits of the spherical bessel functions j0(x) and j1(x)
j0 = 1.0;
j1 = 0.0;
} else {
@ -476,10 +739,14 @@ double TLFDynLorKT::operator()(double t, const vector<double> &par) const {
if(fCalcNeeded){
double tt(0.);
int i;
if(TMath::Abs(par[0])<0.00135538817){
for(unsigned int i(0); i<fNSteps; i++) {
tt=double(i)*fDt;
#ifdef HAVE_GOMP
#pragma omp parallel for default(shared) private(i, tt) schedule(dynamic)
#endif
for(i = 0; i < static_cast<int>(fNSteps); ++i) {
tt=static_cast<double>(i)*fDt;
fFFTtime[i]=(0.33333333333333333333+0.66666666666666666667*(1.0-par[1]*tt)*TMath::Exp(-par[1]*tt))*TMath::Exp(-(fC+par[2])*tt)*fDt;
}
} else {
@ -490,7 +757,10 @@ double TLFDynLorKT::operator()(double t, const vector<double> &par) const {
fFFTtime[0] = 1.0*fDt;
for(unsigned int i(1); i<fNSteps; i++) {
#ifdef HAVE_GOMP
#pragma omp parallel for default(shared) private(i, tt) schedule(dynamic)
#endif
for(i = 1; i < static_cast<int>(fNSteps); ++i) {
tt=(double(i)-0.5)*fDt;
fFFTtime[i]=TMath::Sin(omegaL*tt)/(omegaL*tt)*TMath::Exp(-par[1]*tt)*fDt;
}
@ -498,7 +768,7 @@ double TLFDynLorKT::operator()(double t, const vector<double> &par) const {
double totoIntegrale(0.);
double w0t(1.0);
for(unsigned int i(1); i<fNSteps; i++) {
for(i = 1; i < static_cast<int>(fNSteps); ++i) {
tt=double(i)*fDt;
totoIntegrale+=fFFTtime[i];
w0t = omegaL*tt;
@ -515,7 +785,10 @@ double TLFDynLorKT::operator()(double t, const vector<double> &par) const {
double nusq(par[2]*par[2]); // nu^2
double denom(1.0), imagsq(0.0), oneMINrealnu(1.0);
for (unsigned int i(0); i<fNSteps/2+1; i++) {
#ifdef HAVE_GOMP
#pragma omp parallel for default(shared) private(i, imagsq, oneMINrealnu, denom) schedule(dynamic)
#endif
for (i = 0; i < static_cast<int>(fNSteps)/2+1; ++i) {
imagsq=fFFTfreq[i][1]*fFFTfreq[i][1];
oneMINrealnu=1.0-(par[2]*fFFTfreq[i][0]);
denom=oneMINrealnu*oneMINrealnu + (nusq*imagsq);
@ -531,47 +804,88 @@ double TLFDynLorKT::operator()(double t, const vector<double> &par) const {
// fFFTtime[i]=(fDw*TMath::Exp(fC*i*fDt)/TMath::Pi()*fFFTtime[i]);
// }
fCalcNeeded=false;
fCounter++;
++fCounter;
}
// return fFFTtime[int(t/fDt)];
return fDw*TMath::Exp(fC*t)/TMath::Pi()*fFFTtime[int(t/fDt)];
return fDw*TMath::Exp(fC*t)/TMath::Pi()*fFFTtime[static_cast<int>(t/fDt)];
}
// De Renzi Spin Glass Interpolation
TLFSGInterpolation::TLFSGInterpolation() {
TIntSGInterpolation *integral = new TIntSGInterpolation();
fIntegral = integral;
integral = 0;
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
* <p>Constructor: Initialize variables and allocate memory for the integrator.
*/
TLFSGInterpolation::TLFSGInterpolation() : fCalcNeeded(true), fLastTime(0.) {
fIntegral = new TIntSGInterpolation();
fPar.resize(4, 0.);
}
//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
/**
* <p>Destructor: Free memory of the integrator, clean up some vectors.
*/
TLFSGInterpolation::~TLFSGInterpolation() {
delete fIntegral;
fIntegral = 0;
delete fIntegral;
fIntegral = 0;
fPar.clear();
fIntValues.clear();
}
//--------------------------------------------------------------------------
// operator()
//--------------------------------------------------------------------------
/**
* <p>Method returning the function value at a given time for a given set of parameters.
*
* \param t time \htmlonly (&#956;s) \endhtmlonly \latexonly ($\mu\mathrm{s}$) \endlatexonly
* \param par parameters [\htmlonly &#957;<sub>L</sub>=<i>B</i>&#947;<sub>&#956;</sub>/2&#960; (MHz), <i>a</i> (&#956;s<sup>-1</sup>), &#955; (&#956;s<sup>-1</sup>), &#946; (1) \endhtmlonly \latexonly $\nu_{\mathrm{L}}=B\gamma_{\mu}/2\pi~(\mathrm{MHz})$, $a~(\mu\mathrm{s}^{-1})$, $\lambda~(\mu\mathrm{s}^{-1})$, $\beta~(1)$ \endlatexonly]
*/
double TLFSGInterpolation::operator()(double t, const vector<double> &par) const {
assert(par.size() == 4); // two parameters nu=gbar*B [MHz], a [1/us], lambda [1/us], beta [1]
assert(par.size() == 4); // four parameters nu=gbar*B [MHz], a [1/us], lambda [1/us], beta [1]
if((t <= 0.0) || (par[3] <= 0.0))
return 1.0;
double expo(0.5*par[1]*par[1]*t*t/par[3]+par[2]*t);
if(TMath::Abs(par[0])<0.00135538817){
if (TMath::Abs(par[0]) < 0.00135538817) {
return (0.33333333333333333333*TMath::Exp(-TMath::Power(par[2]*t,par[3])) + \
0.66666666666666666667*(1.0-par[1]*par[1]*t*t/TMath::Power(expo,(1.0-par[3])))*TMath::Exp(-TMath::Power(expo,par[3])));
}
vector<double> intpar(par);
intpar.push_back(t);
fIntegral->SetParameters(intpar);
for (unsigned int i(0); i<4; ++i) {
if( fPar[i]-par[i] ) {
fPar[i] = par[i];
fCalcNeeded = true;
}
}
if (fCalcNeeded) {
fIntValues.clear();
fIntValues[0.] = 0.;
fLastTime = 0.;
} else {
if (t > fIntValues.rbegin()->first) {
fCalcNeeded = true;
fLastTime = fIntValues.rbegin()->first;
}
}
double omegaL(TMath::TwoPi()*par[0]);
if (fCalcNeeded) {
vector<double> intpar(par);
intpar.push_back(t);
fIntegral->SetParameters(intpar);
fIntValues[t] = fIntValues[fLastTime] + fIntegral->IntegrateFunc(fLastTime,t);
intpar.clear();
}
return (TMath::Exp(-TMath::Power(par[2]*t,par[3])) + 2.0*par[1]*par[1]/(omegaL*omegaL) * \
((TMath::Cos(omegaL*t)-TMath::Sin(omegaL*t)/(omegaL*t))*TMath::Exp(-TMath::Power(expo,par[3]))/TMath::Power(expo,(1.0-par[3])) + \
fIntegral->IntegrateFunc(0.,t)));
fIntValues[t]));
}

194
src/external/libLFRelaxation/TLFRelaxation.h vendored
View File

@ -5,7 +5,7 @@
Author: Bastian M. Wojek
e-mail: bastian.wojek@psi.ch
2008/12/04
$Id$
***************************************************************************/
@ -35,6 +35,7 @@
#include<vector>
#include<cstdio>
#include<cmath>
#include<map>
using namespace std;
@ -50,6 +51,15 @@ using namespace std;
#include "fftw3.h"
#include "../BMWIntegrator/BMWIntegrator.h"
//-----------------------------------------------------------------------------------------------------------------
/**
* <p>User function for a static Gaussian depolarization in a longitudinal field using direct integration through GSL
* See also: R. S. Hayano, Y. J. Uemura, J. Imazato, N. Nishida, T. Yamazaki, and R. Kubo,
* \htmlonly Phys. Rev. B <b>20</b>, 850&#150;859 (1979), doi:<a href="http://dx.doi.org/10.1103/PhysRevB.20.850">10.1103/PhysRevB.20.850</a>
* \endhtmlonly
* \latexonly Phys. Rev. B \textbf{20}, 850--859 (1979), \texttt{http://dx.doi.org/10.1103/PhysRevB.20.850}
* \endlatexonly
*/
class TLFStatGssKT : public PUserFcnBase {
public:
@ -63,16 +73,29 @@ public:
double operator()(double, const vector<double>&) const;
private:
TIntSinGss *fIntSinGss;
TIntSinGss *fIntSinGss; ///< integrator
mutable map<double, double> fIntValues; ///< previously calculated integral values
mutable vector<double> fPar; ///< parameters of the function [\htmlonly &#957;<sub>L</sub>=<i>B</i>&#947;<sub>&#956;</sub>/2&#960; (MHz), &#963; (&#956;s<sup>-1</sup>)\endhtmlonly \latexonly $\nu_{\mathrm{L}}=B\gamma_{\mu}/2\pi~(\mathrm{MHz})$, $\sigma~(\mu\mathrm{s}^{-1})$ \endlatexonly]
mutable bool fCalcNeeded; ///< flag specifying if an integration has to be done
mutable double fLastTime; ///< time of the last calculated function value
ClassDef(TLFStatGssKT,1)
ClassDef(TLFStatGssKT,2)
};
class TLFStatLorKT : public PUserFcnBase {
//-----------------------------------------------------------------------------------------------------------------
/**
* <p>User function for a static exponential depolarization in a longitudinal field using direct integration through GSL
* See also: Y. J. Uemura, T. Yamazaki, D. R. Harshman, M. Senba, and E. J. Ansaldo,
* \htmlonly Phys. Rev. B <b>31</b>, 546&#150;563 (1985), doi:<a href="http://dx.doi.org/10.1103/PhysRevB.31.546">10.1103/PhysRevB.31.546</a>
* \endhtmlonly
* \latexonly Phys. Rev. B \textbf{31}, 546--563 (1985), \texttt{http://dx.doi.org/10.1103/PhysRevB.31.546}
* \endlatexonly
*/
class TLFStatExpKT : public PUserFcnBase {
public:
TLFStatLorKT();
~TLFStatLorKT();
TLFStatExpKT();
~TLFStatExpKT();
Bool_t NeedGlobalPart() const { return false; }
void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@ -81,11 +104,31 @@ public:
double operator()(double, const vector<double>&) const;
private:
TIntBesselJ0Exp *fIntBesselJ0Exp;
TIntBesselJ0Exp *fIntBesselJ0Exp; ///< integrator
mutable map<double, double> fIntValues; ///< previously calculated integral values
mutable vector<double> fPar; ///< parameters of the function [\htmlonly &#957;<sub>L</sub>=<i>B</i>&#947;<sub>&#956;</sub>/2&#960; (MHz), <i>a</i> (&#956;s<sup>-1</sup>)\endhtmlonly \latexonly $\nu_{\mathrm{L}}=B\gamma_{\mu}/2\pi~(\mathrm{MHz})$, $a~(\mu\mathrm{s}^{-1})$ \endlatexonly]
mutable bool fCalcNeeded; ///< flag specifying if an integration has to be done
mutable double fLastTime; ///< time of the last calculated function value
ClassDef(TLFStatLorKT,1)
ClassDef(TLFStatExpKT,2)
};
//-----------------------------------------------------------------------------------------------------------------
/**
* <p>User function for a dynamic Gaussian depolarization in a longitudinal field calculated using Laplace transforms
* See also: R. S. Hayano, Y. J. Uemura, J. Imazato, N. Nishida, T. Yamazaki, and R. Kubo,
* \htmlonly Phys. Rev. B <b>20</b>, 850&#150;859 (1979), doi:<a href="http://dx.doi.org/10.1103/PhysRevB.20.850">10.1103/PhysRevB.20.850</a>
* \endhtmlonly
* \latexonly Phys. Rev. B \textbf{20}, 850--859 (1979), \texttt{http://dx.doi.org/10.1103/PhysRevB.20.850}
* \endlatexonly
*
* For details regarding the numerical Laplace transform, see e.g.
* P. Moreno and A. Ramirez,
* \htmlonly IEEE Trans. Power Delivery <b>23</b>, 2599&#150;2609 (2008), doi:<a href="http://dx.doi.org/10.1109/TPWRD.2008.923404">10.1109/TPWRD.2008.923404</a>
* \endhtmlonly
* \latexonly IEEE Trans. Power Delivery \textbf{23}, 2599--2609 (2008), \texttt{http://dx.doi.org/10.1109/TPWRD.2008.923404}
* \endlatexonly
*/
class TLFDynGssKT : public PUserFcnBase {
public:
@ -99,28 +142,51 @@ public:
double operator()(double, const vector<double>&) const;
private:
mutable vector<double> fPar;
mutable bool fCalcNeeded;
mutable bool fFirstCall;
string fWisdom;
unsigned int fNSteps;
double fDt;
double fDw;
double fC;
fftwf_plan fFFTplanFORW;
fftwf_plan fFFTplanBACK;
float *fFFTtime;
fftwf_complex *fFFTfreq;
mutable unsigned int fCounter;
mutable vector<double> fPar; ///< parameters of the function [\htmlonly &#957;<sub>L</sub>=<i>B</i>&#947;<sub>&#956;</sub>/2&#960; (MHz), &#963; (&#956;s<sup>-1</sup>), &#957; (MHz)\endhtmlonly \latexonly $\nu_{\mathrm{L}}=B\gamma_{\mu}/2\pi~(\mathrm{MHz})$, $\sigma~(\mu\mathrm{s}^{-1})$, $\nu~(\mathrm{MHz})$ \endlatexonly]
mutable bool fCalcNeeded; ///< flag indicating if the expensive Laplace transform has to be done (e.g. after parameters have changed)
mutable bool fFirstCall; ///< flag indicating if the function is evaluated for the first time
bool fUseWisdom; ///< flag showing if a FFTW3-wisdom file is used
string fWisdom; ///< path to the wisdom file
unsigned int fNSteps; ///< length of the Laplace transform
double fDt; ///< time resolution
double fDw; ///< frequency resolution
double fC; ///< coefficient depending on the length of the Laplace transform and time resolution
fftwf_plan fFFTplanFORW; ///< FFTW3 plan: time domain \htmlonly &#8594; \endhtmlonly \latexonly $\rightarrow$ \endlatexonly frequency domain
fftwf_plan fFFTplanBACK; ///< FFTW3 plan: time domain \htmlonly &#8592; \endhtmlonly \latexonly $\leftarrow$ \endlatexonly frequency domain
float *fFFTtime; ///< time-domain array
fftwf_complex *fFFTfreq; ///< frequency-domain array
mutable unsigned int fCounter; ///< counter determining how many Laplace transforms are done (mainly for debugging purposes)
ClassDef(TLFDynGssKT,1)
ClassDef(TLFDynGssKT,2)
};
class TLFDynLorKT : public PUserFcnBase {
//-----------------------------------------------------------------------------------------------------------------
/**
* <p>User function for a dynamic exponential depolarization in a longitudinal field calculated using Laplace transforms
* See also: Y. J. Uemura, T. Yamazaki, D. R. Harshman, M. Senba, and E. J. Ansaldo,
* \htmlonly Phys. Rev. B <b>31</b>, 546&#150;563 (1985), doi:<a href="http://dx.doi.org/10.1103/PhysRevB.31.546">10.1103/PhysRevB.31.546</a>
* \endhtmlonly
* \latexonly Phys. Rev. B \textbf{31}, 546--563 (1985), \texttt{http://dx.doi.org/10.1103/PhysRevB.31.546}
* \endlatexonly
*
* and: R. S. Hayano, Y. J. Uemura, J. Imazato, N. Nishida, T. Yamazaki, and R. Kubo,
* \htmlonly Phys. Rev. B <b>20</b>, 850&#150;859 (1979), doi:<a href="http://dx.doi.org/10.1103/PhysRevB.20.850">10.1103/PhysRevB.20.850</a>
* \endhtmlonly
* \latexonly Phys. Rev. B \textbf{20}, 850--859 (1979), \texttt{http://dx.doi.org/10.1103/PhysRevB.20.850}
* \endlatexonly
*
* For details regarding the numerical Laplace transform, see e.g.
* P. Moreno and A. Ramirez,
* \htmlonly IEEE Trans. Power Delivery <b>23</b>, 2599&#150;2609 (2008), doi:<a href="http://dx.doi.org/10.1109/TPWRD.2008.923404">10.1109/TPWRD.2008.923404</a>
* \endhtmlonly
* \latexonly IEEE Trans. Power Delivery \textbf{23}, 2599--2609 (2008), \texttt{http://dx.doi.org/10.1109/TPWRD.2008.923404}
* \endlatexonly
*/
class TLFDynExpKT : public PUserFcnBase {
public:
TLFDynLorKT();
~TLFDynLorKT();
TLFDynExpKT();
~TLFDynExpKT();
Bool_t NeedGlobalPart() const { return false; }
void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@ -129,25 +195,48 @@ public:
double operator()(double, const vector<double>&) const;
private:
mutable vector<double> fPar;
mutable bool fCalcNeeded;
mutable bool fFirstCall;
string fWisdom;
unsigned int fNSteps;
double fDt;
double fDw;
double fC;
fftwf_plan fFFTplanFORW;
fftwf_plan fFFTplanBACK;
float *fFFTtime;
fftwf_complex *fFFTfreq;
mutable unsigned int fCounter;
mutable double fL1;
mutable double fL2;
mutable vector<double> fPar; ///< parameters of the function [\htmlonly &#957;<sub>L</sub>=<i>B</i>&#947;<sub>&#956;</sub>/2&#960; (MHz), <i>a</i> (&#956;s<sup>-1</sup>), &#957; (MHz)\endhtmlonly \latexonly $\nu_{\mathrm{L}}=B\gamma_{\mu}/2\pi~(\mathrm{MHz})$, $a~(\mu\mathrm{s}^{-1})$, $\nu~(\mathrm{MHz})$ \endlatexonly]
mutable bool fCalcNeeded; ///< flag indicating if the expensive Laplace transform has to be done (e.g. after parameters have changed)
mutable bool fFirstCall; ///< flag indicating if the function is evaluated for the first time
bool fUseWisdom; ///< flag showing if a FFTW3-wisdom file is used
string fWisdom; ///< path to the wisdom file
unsigned int fNSteps; ///< length of the Laplace transform
double fDt; ///< time resolution
double fDw; ///< frequency resolution
double fC; ///< coefficient depending on the length of the Laplace transform and time resolution
fftwf_plan fFFTplanFORW; ///< FFTW3 plan: time domain \htmlonly &#8594; \endhtmlonly \latexonly $\rightarrow$ \endlatexonly frequency domain
fftwf_plan fFFTplanBACK; ///< FFTW3 plan: time domain \htmlonly &#8592; \endhtmlonly \latexonly $\leftarrow$ \endlatexonly frequency domain
float *fFFTtime; ///< time-domain array
fftwf_complex *fFFTfreq; ///< frequency-domain array
mutable unsigned int fCounter; ///< counter determining how many Laplace transforms are done (mainly for debugging purposes)
mutable double fL1; ///< coefficient used for the high-field and high-hopping-rate approximation
mutable double fL2; ///< coefficient used for the high-field and high-hopping-rate approximation
ClassDef(TLFDynLorKT,1)
ClassDef(TLFDynExpKT,2)
};
//-----------------------------------------------------------------------------------------------------------------
/**
* <p>User function for a dynamic depolarization in a spin glass in a longitudinal field calculated using Laplace transforms
* See also: Y. J. Uemura, T. Yamazaki, D. R. Harshman, M. Senba, and E. J. Ansaldo,
* \htmlonly Phys. Rev. B <b>31</b>, 546&#150;563 (1985), doi:<a href="http://dx.doi.org/10.1103/PhysRevB.31.546">10.1103/PhysRevB.31.546</a>
* \endhtmlonly
* \latexonly Phys. Rev. B \textbf{31}, 546--563 (1985), \texttt{http://dx.doi.org/10.1103/PhysRevB.31.546}
* \endlatexonly
*
* and: R. S. Hayano, Y. J. Uemura, J. Imazato, N. Nishida, T. Yamazaki, and R. Kubo,
* \htmlonly Phys. Rev. B <b>20</b>, 850&#150;859 (1979), doi:<a href="http://dx.doi.org/10.1103/PhysRevB.20.850">10.1103/PhysRevB.20.850</a>
* \endhtmlonly
* \latexonly Phys. Rev. B \textbf{20}, 850--859 (1979), \texttt{http://dx.doi.org/10.1103/PhysRevB.20.850}
* \endlatexonly
*
* For details regarding the numerical Laplace transform, see e.g.
* P. Moreno and A. Ramirez,
* \htmlonly IEEE Trans. Power Delivery <b>23</b>, 2599&#150;2609 (2008), doi:<a href="http://dx.doi.org/10.1109/TPWRD.2008.923404">10.1109/TPWRD.2008.923404</a>
* \endhtmlonly
* \latexonly IEEE Trans. Power Delivery \textbf{23}, 2599--2609 (2008), \texttt{http://dx.doi.org/10.1109/TPWRD.2008.923404}
* \endlatexonly
*/
class TLFDynSG : public PUserFcnBase {
public:
@ -161,12 +250,21 @@ public:
double operator()(double, const vector<double>&) const;
private:
mutable vector<double> fPar;
vector<TLFDynGssKT*> fLFDynGssIntegral;
mutable vector<double> fPar; ///< parameters of the dynamic Gaussian depolarization function [\htmlonly &#957;<sub>L</sub>=<i>B</i>&#947;<sub>&#956;</sub>/2&#960; (MHz), &#963; (&#956;s<sup>-1</sup>), &#957; (MHz)\endhtmlonly \latexonly $\nu_{\mathrm{L}}=B\gamma_{\mu}/2\pi~(\mathrm{MHz})$, $\sigma~(\mu\mathrm{s}^{-1})$, $\nu~(\mathrm{MHz})$ \endlatexonly]
vector<TLFDynGssKT*> fLFDynGssIntegral; ///< vector of dynamic Gaussian depolarization functions for a subset of static widths
ClassDef(TLFDynSG,1)
};
//-----------------------------------------------------------------------------------------------------------------
/**
* <p>User function for a dynamic depolarization in a spin glass in a longitudinal field
* See also: R. De Renzi and S. Fanesi
* \htmlonly Physica B <b>289&#150;290</b>, 209&#150;212 (2000), doi:<a href="http://dx.doi.org/10.1016/S0921-4526(00)00368-9">10.1016/S0921-4526(00)00368-9</a>
* \endhtmlonly
* \latexonly Physica B \textbf{289--290}, 209--212 (2000), \texttt{http://dx.doi.org/10.1016/S0921-4526(00)00368-9}
* \endlatexonly
*/
class TLFSGInterpolation : public PUserFcnBase {
public:
@ -180,10 +278,14 @@ public:
double operator()(double, const vector<double>&) const;
private:
TIntSGInterpolation *fIntegral;
TIntSGInterpolation *fIntegral; ///< integrator
mutable map<double, double> fIntValues; ///< previously calculated integral values
mutable vector<double> fPar; ///< parameters of the function [\htmlonly &#957;<sub>L</sub>=<i>B</i>&#947;<sub>&#956;</sub>/2&#960; (MHz), <i>a</i> (&#956;s<sup>-1</sup>), &#955; (&#956;s<sup>-1</sup>), &#946; (1) \endhtmlonly \latexonly $\nu_{\mathrm{L}}=B\gamma_{\mu}/2\pi~(\mathrm{MHz})$, $a~(\mu\mathrm{s}^{-1})$, $\lambda~(\mu\mathrm{s}^{-1})$, $\beta~(1)$ \endlatexonly]
mutable bool fCalcNeeded; ///< flag specifying if an integration has to be done
mutable double fLastTime; ///< time of the last calculated function value
ClassDef(TLFSGInterpolation,1)
ClassDef(TLFSGInterpolation,2)
};
#endif //_LFRelaxation_H_
#endif //_TLFRelaxation_H_

6
src/external/libLFRelaxation/TLFRelaxationLinkDef.h vendored
View File

@ -5,7 +5,7 @@
Author: Bastian M. Wojek
e-mail: bastian.wojek@psi.ch
2008/12/04
$Id$
***************************************************************************/
@ -37,9 +37,9 @@
#pragma link off all functions;
#pragma link C++ class TLFStatGssKT+;
#pragma link C++ class TLFStatLorKT+;
#pragma link C++ class TLFStatExpKT+;
#pragma link C++ class TLFDynGssKT+;
#pragma link C++ class TLFDynLorKT+;
#pragma link C++ class TLFDynExpKT+;
#pragma link C++ class TLFDynSG+;
#pragma link C++ class TLFSGInterpolation+;