Added semi-analytical approximations to the dynamical Lorentz LF class.

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

src/external/libLFRelaxation/TIntegrator.cpp

@@ -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;
-}
-
-
-

src/external/libLFRelaxation/TIntegrator.h

@@ -15,7 +15,7 @@
 #include "Math/GSLIntegrator.h"
 #include "TMath.h"
 
-#include<vector>
+#include <vector>
 
 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);

src/external/libLFRelaxation/TLFRelaxation.cpp

@@ -10,17 +10,24 @@
 ***************************************************************************/
 
 #include <cassert>
+#include <sys/time.h>
+#include <iostream>
 
 #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<double> &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;
+  fDw = PI/fNSteps/fDt;
-  fC = 2.0*TMath::Log(double(fNSteps))/(double(fNSteps-1)*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);
+  fFFTplanFORW = fftw_plan_dft_r2c_1d(fNSteps, fFFTtime, fFFTfreq, FFTW_ESTIMATE);
-  fFFTplanBACK = fftw_plan_dft_c2r_1d(fNSteps, fFFTfreq, fFFTtime, FFTW_EXHAUSTIVE);
+  fFFTplanBACK = fftw_plan_dft_c2r_1d(fNSteps, fFFTfreq, fFFTtime, FFTW_ESTIMATE);
 }
 
 TLFDynGssKT::~TLFDynGssKT() {
@@ -167,39 +174,53 @@ double TLFDynGssKT::operator()(double t, const vector<double> &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(fabs(par[0])<0.00135538817){ // Zero Field
-    if(par[2]/par[1] > 5. || omegaL > 20.0*par[1]){
+    if(!nusq){
-      if(TMath::Abs(par[0])<0.00135538817){
+      double sigsqtsq(sigsq*t*t);
-        return TMath::Exp(-2.0*sigsq/nusq*(TMath::Exp(-par[2]*t)-1.0+par[2]*t)); // ZF Abragam Delta^2->2*Delta^2
+      return 0.33333333333333333333+0.66666666666666666667*(1.0-sigsqtsq)*exp(-0.5*sigsqtsq); // ZF static Gauss KT
-      }
-      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
     }
+    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<fNSteps; i++) {
         tt=double(i)*fDt;
-        fFFTtime[i]=(0.33333333333333333333+0.66666666666666666667*(1.0-sigsq*tt*tt)*TMath::Exp(-0.5*sigsq*tt*tt))*TMath::Exp(-(fC+par[2])*tt)*fDt;
+        sigsqtsq=sigsq*tt*tt;
+        fFFTtime[i]=(0.33333333333333333333+0.66666666666666666667*(1.0-sigsqtsq)*exp(-0.5*sigsqtsq))*exp(mcplusnu*tt)*fDt;
       }
     } else {
+      double mcplusnu(-(fC+par[2]));
+      double coeff1(2.0*sigsq/omegaLsq);
+      double coeff2(coeff1*sigsq/omegaL);
 
-      double coeff1(2.0*sigsq/(omegaL*omegaL));
+      fFFTtime[0] = fDt;
-      double coeff2(2.0*sigsq*sigsq/(omegaL*omegaL*omegaL));
-
-      fFFTtime[0] = 1.0*fDt;
 
       for(unsigned int i(1); i<fNSteps; i++) {
         tt=(double(i)-0.5)*fDt;
-        fFFTtime[i]=(TMath::Sin(omegaL*tt) * TMath::Exp(-0.5*sigsq*tt*tt))*fDt;
+        fFFTtime[i]=(gsl_sf_sin(omegaL*tt) * exp(-0.5*sigsq*tt*tt))*fDt;
       }
 
       double totoIntegrale(0.);
@@ -207,22 +228,27 @@ double TLFDynGssKT::operator()(double t, const vector<double> &par) const {
       for(unsigned int i(1); i<fNSteps; i++) {
         tt=double(i)*fDt;
         totoIntegrale+=fFFTtime[i];
-        fFFTtime[i]=(1.0-(coeff1*(1.0-TMath::Exp(-0.5*sigsq*tt*tt)*TMath::Cos(omegaL*tt)))+(coeff2*totoIntegrale))*TMath::Exp(-(fC+par[2])*tt)*fDt;
+        fFFTtime[i]=(1.0-(coeff1*(1.0-exp(-0.5*sigsq*tt*tt)*gsl_sf_cos(omegaL*tt)))+(coeff2*totoIntegrale))*exp(mcplusnu*tt)*fDt;
       }
     }
 
+    gettimeofday(&tv_stop, 0);
+    t1 = (tv_stop.tv_sec - tv_start.tv_sec)*1000.0 + (tv_stop.tv_usec - tv_start.tv_usec)/1000.0;
+
     // Transform to frequency domain
 
     fftw_execute(fFFTplanFORW);
 
     // calculate F(s)
 
-    double denom(1.0);
+    double denom(1.0), imagsq(0.0), oneMINrealnu(1.0);
 
     for (unsigned int i(0); i<fNSteps/2+1; i++) {
-      denom=(1.0-(par[2]*fFFTfreq[i][0]))*(1.0-(par[2]*fFFTfreq[i][0])) + (nusq*fFFTfreq[i][1]*fFFTfreq[i][1]);
+      imagsq=fFFTfreq[i][1]*fFFTfreq[i][1];
-      fFFTfreq[i][0] = (fFFTfreq[i][0]-(par[2]*fFFTfreq[i][0]*fFFTfreq[i][0])-(par[2]*fFFTfreq[i][1]*fFFTfreq[i][1]))/denom;
+      oneMINrealnu=1.0-(par[2]*fFFTfreq[i][0]);
-      fFFTfreq[i][1] = fFFTfreq[i][1]/denom;
+      denom=oneMINrealnu*oneMINrealnu + (nusq*imagsq);
+      fFFTfreq[i][0] = (oneMINrealnu*fFFTfreq[i][0]-(par[2]*imagsq))/denom;
+      fFFTfreq[i][1] /= denom;
     }
 
     // Transform back to time domain
@@ -234,14 +260,19 @@ double TLFDynGssKT::operator()(double t, const vector<double> &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<double> &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<double> &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<double> &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<double> &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<fNSteps/2+1; i++) {
-      denom=(1.0-(par[2]*fFFTfreq[i][0]))*(1.0-(par[2]*fFFTfreq[i][0])) + (nusq*fFFTfreq[i][1]*fFFTfreq[i][1]);
+      imagsq=fFFTfreq[i][1]*fFFTfreq[i][1];
-      fFFTfreq[i][0] = (fFFTfreq[i][0]-(par[2]*fFFTfreq[i][0]*fFFTfreq[i][0])-(par[2]*fFFTfreq[i][1]*fFFTfreq[i][1]))/denom;
+      oneMINrealnu=1.0-(par[2]*fFFTfreq[i][0]);
-      fFFTfreq[i][1] = fFFTfreq[i][1]/denom;
+      denom=oneMINrealnu*oneMINrealnu + (nusq*imagsq);
+      fFFTfreq[i][0] = (oneMINrealnu*fFFTfreq[i][0]-(par[2]*imagsq))/denom;
+      fFFTfreq[i][1] /= denom;
     }
 
     // Transform back to time domain

src/external/libLFRelaxation/TLFRelaxation.h

@@ -14,10 +14,17 @@
 
 #include<vector>
 #include<cstdio>
+#include<cmath>
 
 using namespace std;
 
-#include "TMath.h"
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_sf_exp.h>
+#include <gsl/gsl_sf_log.h>
+#include <gsl/gsl_sf_trig.h>
+#include <gsl/gsl_sf_bessel.h>
+
+//#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_