From f9971d4858836f7fc2447bd06d3c1126ad0e4acb Mon Sep 17 00:00:00 2001 From: "Bastian M. Wojek" Date: Wed, 4 Feb 2009 15:07:45 +0000 Subject: [PATCH] Added semi-analytical approximations to the dynamical Lorentz LF class. --- .../libLFRelaxation/Makefile.libLFRelaxation | 1 - src/external/libLFRelaxation/TIntegrator.cpp | 30 ---- src/external/libLFRelaxation/TIntegrator.h | 15 +- .../libLFRelaxation/TLFRelaxation.cpp | 154 +++++++++++++----- src/external/libLFRelaxation/TLFRelaxation.h | 15 +- 5 files changed, 145 insertions(+), 70 deletions(-) delete mode 100644 src/external/libLFRelaxation/TIntegrator.cpp diff --git a/src/external/libLFRelaxation/Makefile.libLFRelaxation b/src/external/libLFRelaxation/Makefile.libLFRelaxation index 35d53b1f..bd752a6a 100644 --- a/src/external/libLFRelaxation/Makefile.libLFRelaxation +++ b/src/external/libLFRelaxation/Makefile.libLFRelaxation @@ -24,7 +24,6 @@ LDFLAGS += # some definitions: headers (used to generate *Dict* stuff), sources, objects,... OBJS = -OBJS += TIntegrator.o OBJS += TLFRelaxation.o TLFRelaxationDict.o SHLIB = libLFRelaxation.so diff --git a/src/external/libLFRelaxation/TIntegrator.cpp b/src/external/libLFRelaxation/TIntegrator.cpp deleted file mode 100644 index 3ad5859a..00000000 --- a/src/external/libLFRelaxation/TIntegrator.cpp +++ /dev/null @@ -1,30 +0,0 @@ -/*************************************************************************** - - TIntegrator.cpp - - Author: Bastian M. Wojek - e-mail: bastian.wojek@psi.ch - - 2008/12/24 - -***************************************************************************/ - -#include "TIntegrator.h" - -using namespace std; - -TIntegrator::TIntegrator() : fFunc(0) { - ROOT::Math::GSLIntegrator *integrator = new ROOT::Math::GSLIntegrator(ROOT::Math::Integration::kADAPTIVE,ROOT::Math::Integration::kGAUSS51); - fIntegrator = integrator; - integrator = 0; - delete integrator; -} - -TIntegrator::~TIntegrator(){ - delete fIntegrator; - fIntegrator=0; - fFunc=0; -} - - - diff --git a/src/external/libLFRelaxation/TIntegrator.h b/src/external/libLFRelaxation/TIntegrator.h index 752fc1ee..ea96d2e1 100644 --- a/src/external/libLFRelaxation/TIntegrator.h +++ b/src/external/libLFRelaxation/TIntegrator.h @@ -15,7 +15,7 @@ #include "Math/GSLIntegrator.h" #include "TMath.h" -#include +#include using namespace std; @@ -38,6 +38,19 @@ class TIntegrator { mutable double (*fFunc)(double, void *); }; +inline TIntegrator::TIntegrator() : fFunc(0) { + ROOT::Math::GSLIntegrator *integrator = new ROOT::Math::GSLIntegrator(ROOT::Math::Integration::kADAPTIVE,ROOT::Math::Integration::kGAUSS31); + fIntegrator = integrator; + integrator = 0; + delete integrator; +} + +inline TIntegrator::~TIntegrator(){ + delete fIntegrator; + fIntegrator=0; + fFunc=0; +} + inline double TIntegrator::FuncAtXgsl(double x, void *obj) { return ((TIntegrator*)obj)->FuncAtX(x); diff --git a/src/external/libLFRelaxation/TLFRelaxation.cpp b/src/external/libLFRelaxation/TLFRelaxation.cpp index 2c02b49c..9fde97ce 100644 --- a/src/external/libLFRelaxation/TLFRelaxation.cpp +++ b/src/external/libLFRelaxation/TLFRelaxation.cpp @@ -10,17 +10,24 @@ ***************************************************************************/ #include +#include +#include #include "TIntegrator.h" #include "TLFRelaxation.h" using namespace std; +#define PI 3.14159265358979323846 +#define TWOPI 6.28318530717958647692 + + ClassImp(TLFStatGssKT) ClassImp(TLFStatLorKT) ClassImp(TLFDynGssKT) ClassImp(TLFDynLorKT) + // LF Static Gaussian KT TLFStatGssKT::TLFStatGssKT() { @@ -94,10 +101,10 @@ double TLFStatLorKT::operator()(double t, const vector &par) const { // LF Dynamic Gaussian KT -TLFDynGssKT::TLFDynGssKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("/home/l_wojek/analysis/WordsOfWisdom.dat"), fNSteps(524288), fDt(0.000040), fCounter(0) { +TLFDynGssKT::TLFDynGssKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("/home/l_wojek/analysis/WordsOfWisdom.dat"), fNSteps(524288), fDt(0.00004), fCounter(0) { // Calculate d_omega and C for given NFFT and dt - fDw = TMath::Pi()/fNSteps/fDt; - fC = 2.0*TMath::Log(double(fNSteps))/(double(fNSteps-1)*fDt); + fDw = PI/fNSteps/fDt; + fC = 2.0*gsl_sf_log(double(fNSteps))/(double(fNSteps-1)*fDt); // Load FFTW Wisdom int wisdomLoaded(0); @@ -120,8 +127,8 @@ TLFDynGssKT::TLFDynGssKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("/home fFFTtime = (double *)malloc(sizeof(double) * fNSteps); fFFTfreq = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (fNSteps/2+1)); - fFFTplanFORW = fftw_plan_dft_r2c_1d(fNSteps, fFFTtime, fFFTfreq, FFTW_EXHAUSTIVE); - fFFTplanBACK = fftw_plan_dft_c2r_1d(fNSteps, fFFTfreq, fFFTtime, FFTW_EXHAUSTIVE); + fFFTplanFORW = fftw_plan_dft_r2c_1d(fNSteps, fFFTtime, fFFTfreq, FFTW_ESTIMATE); + fFFTplanBACK = fftw_plan_dft_c2r_1d(fNSteps, fFFTfreq, fFFTtime, FFTW_ESTIMATE); } TLFDynGssKT::~TLFDynGssKT() { @@ -167,39 +174,53 @@ double TLFDynGssKT::operator()(double t, const vector &par) const { } double sigsq(par[1]*par[1]); // sigma^2 - double omegaL(TMath::TwoPi()*par[0]); // Larmor frequency + double omegaL(TWOPI*par[0]); // Larmor frequency double nusq(par[2]*par[2]); // nu^2 + double nut(par[2]*t); // nu*t + double omegaLsq(omegaL*omegaL); // omega^2 + double omegaLnusqp(omegaLsq+nusq); //omega^2+nu^2 + double omegaLnusqm(omegaLsq-nusq); //omega^2-nu^2 + double wt(omegaL*t); // omega*t + double enut(exp(-nut)); // exp(-nu*t) - if(par[1]){ - if(par[2]/par[1] > 5. || omegaL > 20.0*par[1]){ - if(TMath::Abs(par[0])<0.00135538817){ - return TMath::Exp(-2.0*sigsq/nusq*(TMath::Exp(-par[2]*t)-1.0+par[2]*t)); // ZF Abragam Delta^2->2*Delta^2 - } - double omegaLnusqp(omegaL*omegaL+nusq); - double omegaLnusqm(omegaL*omegaL-nusq); - return TMath::Exp(-2.0*sigsq/(omegaLnusqp*omegaLnusqp)*(omegaLnusqp*par[2]*t+omegaLnusqm*(1.0-TMath::Exp(-par[2]*t)*TMath::Cos(omegaL*t))-2.0*par[2]*omegaL*TMath::Exp(-par[2]*t)*TMath::Sin(omegaL*t))); // Keren + if(fabs(par[0])<0.00135538817){ // Zero Field + if(!nusq){ + double sigsqtsq(sigsq*t*t); + return 0.33333333333333333333+0.66666666666666666667*(1.0-sigsqtsq)*exp(-0.5*sigsqtsq); // ZF static Gauss KT } + return exp(-2.0*sigsq/nusq*(enut-1.0+nut)); // ZF Abragam Delta^2->2*Delta^2 } - if(fCalcNeeded){ + if((par[2] >= 5.0*par[1]) || (omegaL >= 10.0*par[1])) { + return exp(-2.0*sigsq/(omegaLnusqp*omegaLnusqp)*(omegaLnusqp*nut+omegaLnusqm*(1.0-enut*gsl_sf_cos(wt))-2.0*par[2]*omegaL*enut*gsl_sf_sin(wt))); // Keren + } - double tt(0.); + if(fCalcNeeded) { - if(TMath::Abs(par[0])<0.00135538817){ + double t1,t2; + // get start time + struct timeval tv_start, tv_stop; + gettimeofday(&tv_start, 0); + + double tt(0.),sigsqtsq(0.); + + if(fabs(par[0])<0.00135538817) { + double mcplusnu(-(fC+par[2])); for(unsigned int i(0); i &par) const { for(unsigned int i(1); i &par) const { // } fCalcNeeded=false; 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; + cout << "# Calculation times: " << t1 << " (ms), " << t2 << " (ms)" << endl; + } // return fFFTtime[int(t/fDt)]; - return fDw*TMath::Exp(fC*t)/TMath::Pi()*fFFTtime[int(t/fDt)]; + return fDw*exp(fC*t)/PI*fFFTtime[int(t/fDt)]; } // LF Dynamic Lorentz KT -TLFDynLorKT::TLFDynLorKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("/home/l_wojek/analysis/WordsOfWisdom.dat"), fNSteps(524288), fDt(0.000040), fCounter(0) { +TLFDynLorKT::TLFDynLorKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("/home/l_wojek/analysis/WordsOfWisdom.dat"), fNSteps(524288), fDt(0.000040), fCounter(0), fA(1.0), fL1(0.0), fL2(0.0) { // Calculate d_omega and C for given NFFT and dt fDw = TMath::Pi()/fNSteps/fDt; fC = 2.0*TMath::Log(double(fNSteps))/(double(fNSteps-1)*fDt); @@ -252,14 +283,14 @@ TLFDynLorKT::TLFDynLorKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("/home FILE *wordsOfWisdomR; wordsOfWisdomR = fopen(fWisdom.c_str(), "r"); if (wordsOfWisdomR == NULL) { - cout << "TLFDynGssKT::TLFDynGssKT: Couldn't open wisdom file ..." << endl; + cout << "TLFDynLorKT::TLFDynGssKT: Couldn't open wisdom file ..." << endl; } else { wisdomLoaded = fftw_import_wisdom_from_file(wordsOfWisdomR); fclose(wordsOfWisdomR); } if (!wisdomLoaded) { - cout << "TLFDynGssKT::TLFDynGssKT: No wisdom is imported..." << endl; + cout << "TLFDynLorKT::TLFDynGssKT: No wisdom is imported..." << endl; } // END of WisdomLoading @@ -291,6 +322,10 @@ TLFDynLorKT::~TLFDynLorKT() { cout << "TLFDynLorKT::~TLFDynLorKT(): " << fCounter << " full FFT cyles needed..." << endl; } +const double TLFDynLorKT::fX[16]={1.61849, 1.27059, 0.890315, 6.72608, 0.327258, 0.981975, 1.5964, 2.31023, 2.37689, 5.63945, 0.235734, 1.10617, 1.1664, 0.218402, 0.266562, 0.471331}; + +const double TLFDynLorKT::fY[20]={1.54699, 0.990321, 0.324923, 0.928399, 2.11155, 0.615106, 1.49907, 0.679423, 6.17621, 1.38907, 0.979608, 0.0593631, 0.806427, 3.33144, 1.05059, 1.29922, 0.254661, 0.284348, 0.991419, 0.375213}; + double TLFDynLorKT::operator()(double t, const vector &par) const { assert(par.size() == 3); // three parameters nuL=gbar*B,sigma,fluct.rate nu @@ -313,6 +348,50 @@ double TLFDynLorKT::operator()(double t, const vector &par) const { } } + double omegaL(TWOPI*par[0]); // Larmor frequency + + if(fabs(par[0])<0.00135538817){ // Zero Field + if(!par[2]){ // nu=0 + double at(par[1]*t); + return 0.33333333333333333333+0.66666666666666666667*(1.0-at)*exp(-at); // ZF static Lorentz KT + } + } + + // semi-analytical approximations for {nu >= 4a} and {omegaL >= 10a} + + if(par[2]){ // nu!=0 + if(par[2] >= 4.0*par[1]){ // nu >= 4a + if(fCalcNeeded){ + double wLnu(omegaL/par[2]); // omegaL/nu + double wLnusq(wLnu*wLnu); // (omegaL/nu)^2 + double anu(par[1]/par[2]); // a/nu + + fA=1.0+wLnu*anu*(-fX[0]/(wLnusq+fX[1])+fX[2]/(wLnusq*wLnu+fX[3])+fX[4]/(wLnusq*wLnusq+fX[5])); + fL1=par[1]*(fX[6]/(wLnu+fX[7])+fX[8]/(wLnusq+fX[9])+fX[10]/(wLnusq*wLnusq+fX[11])); + fL2=par[2]*(fX[12]*tanh(fX[13]*wLnu)+fX[14]*wLnu+fX[15]); + + fCalcNeeded=false; + } + return fA*exp(-fL1*t)+(1.0-fA)*exp(-fL2*t)*cos(omegaL*t); + } + if(omegaL >= 10.0*par[1]){ // omegaL >= 10a + if(fCalcNeeded){ + double wLnu(omegaL/par[2]); // omegaL/nu + double wLnusq(wLnu*wLnu); // (omegaL/nu)^2 + double anu(par[1]/par[2]); // a/nu + + fA=1.0-fY[0]*pow(anu,fY[1])/(wLnu+fY[2])+fY[3]*pow(anu,fY[4])/(wLnusq-fY[5])+fY[6]*pow(anu,fY[7])/(wLnusq*wLnu+fY[8]); + fL1=par[1]*(fY[9]/(pow(wLnu,fY[10])+fY[11])-fY[12]/(pow(wLnu,fY[13])+fY[14])); + fL2=par[2]*(fY[15]*tanh(fY[16]*wLnu)+fY[17]*pow(wLnu,fY[18])+fY[19]); + + fCalcNeeded=false; + } + return fA*exp(-fL1*t)+(1.0-fA)*exp(-fL2*t)*cos(omegaL*t); + } + } + + // if no approximation can be used -> Laplace transform + if(fCalcNeeded){ double tt(0.); @@ -324,7 +403,6 @@ double TLFDynLorKT::operator()(double t, const vector &par) const { } } else { - double omegaL(TMath::TwoPi()*par[0]); // Larmor frequency double coeff1(par[1]/omegaL); double coeff2(coeff1*coeff1); double coeff3((1.0+coeff2)*par[1]); @@ -351,13 +429,15 @@ double TLFDynLorKT::operator()(double t, const vector &par) const { // calculate F(s) - double denom(1.0); double nusq(par[2]*par[2]); // nu^2 + double denom(1.0), imagsq(0.0), oneMINrealnu(1.0); for (unsigned int i(0); i #include +#include using namespace std; -#include "TMath.h" +#include +#include +#include +#include +#include + +//#include "TMath.h" #include "PUserFcnBase.h" #include "fftw3.h" #include "TIntegrator.h" @@ -98,8 +105,14 @@ private: double *fFFTtime; fftw_complex *fFFTfreq; mutable unsigned int fCounter; + static const double fX[16]; + static const double fY[20]; + mutable double fA; + mutable double fL1; + mutable double fL2; ClassDef(TLFDynLorKT,1) }; + #endif //_LFRelaxation_H_