some more work towards nonlocal fitting. Not functional yet.
This commit is contained in:
nemu
parent
edd3481f1e
commit
5a2f1cb036
245
src/external/Nonlocal/PNL_PippardFitter.cpp
vendored
245
src/external/Nonlocal/PNL_PippardFitter.cpp
vendored
@ -30,11 +30,13 @@
|
||||
***************************************************************************/
|
||||
|
||||
#include <cassert>
|
||||
#include <cmath>
|
||||
|
||||
#include <iostream>
|
||||
using namespace std;
|
||||
|
||||
#include <TSAXParser.h>
|
||||
#include <TMath.h>
|
||||
|
||||
#include "PNL_PippardFitter.h"
|
||||
|
||||
@ -78,13 +80,28 @@ PNL_PippardFitter::PNL_PippardFitter()
|
||||
}
|
||||
|
||||
// check if everything went fine with the startup handler
|
||||
if (!fStartupHandler->IsValid())
|
||||
if (!fStartupHandler->IsValid()) {
|
||||
cout << endl << "PNL_PippardFitter::PNL_PippardFitter **PANIC ERROR**";
|
||||
cout << endl << " startup handler too unhappy. Will terminate unfriendly, sorry.";
|
||||
cout << endl;
|
||||
assert(false);
|
||||
}
|
||||
|
||||
// load all the TRIM.SP rge-files
|
||||
fRgeHandler = new PNL_RgeHandler(fStartupHandler->GetTrimSpDataPathList());
|
||||
if (!fRgeHandler->IsValid())
|
||||
if (!fRgeHandler->IsValid()) {
|
||||
cout << endl << "PNL_PippardFitter::PNL_PippardFitter **PANIC ERROR**";
|
||||
cout << endl << " rge data handler too unhappy. Will terminate unfriendly, sorry.";
|
||||
cout << endl;
|
||||
assert(false);
|
||||
}
|
||||
|
||||
fPlanPresent = false;
|
||||
fFieldq = 0;
|
||||
fFieldB = 0;
|
||||
fShift = 0;
|
||||
|
||||
f_dx = 0.02;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
@ -99,6 +116,25 @@ PNL_PippardFitter::~PNL_PippardFitter()
|
||||
delete fStartupHandler;
|
||||
fStartupHandler = 0;
|
||||
}
|
||||
if (fRgeHandler) {
|
||||
delete fRgeHandler;
|
||||
fRgeHandler = 0;
|
||||
}
|
||||
|
||||
fPreviousParam.clear();
|
||||
|
||||
if (fPlanPresent) {
|
||||
fftw_destroy_plan(fPlan);
|
||||
}
|
||||
if (fFieldq) {
|
||||
fftw_free(fFieldq);
|
||||
fFieldq = 0;
|
||||
}
|
||||
|
||||
if (fFieldB) {
|
||||
fftw_free(fFieldq);
|
||||
fFieldB = 0;
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
@ -109,8 +145,209 @@ PNL_PippardFitter::~PNL_PippardFitter()
|
||||
*/
|
||||
Double_t PNL_PippardFitter::operator()(Double_t t, const std::vector<Double_t> ¶m) const
|
||||
{
|
||||
// expected parameters: energy, temp, thickness, meanFreePath, xi0, lambdaL
|
||||
assert(param.size() == 6);
|
||||
// expected parameters: energy, temp, thickness, meanFreePath, xi0, lambdaL, phase
|
||||
assert(param.size() == 7);
|
||||
|
||||
if (NewParameters(param)) { // new parameters, hence B(z), P(t), ..., needs to be calculated
|
||||
// keep parameters
|
||||
for (UInt_t i=0; i<param.size(); i++)
|
||||
fPreviousParam[i] = param[i];
|
||||
CalculateField(param);
|
||||
CalculatePolarization(param);
|
||||
}
|
||||
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// NewParameters
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
*
|
||||
*/
|
||||
Bool_t PNL_PippardFitter::NewParameters(const std::vector<Double_t> ¶m) const
|
||||
{
|
||||
if (fPreviousParam.size() == 0) {
|
||||
for (UInt_t i=0; i<param.size(); i++)
|
||||
fPreviousParam.push_back(param[i]);
|
||||
}
|
||||
|
||||
assert(param.size() == fPreviousParam.size());
|
||||
|
||||
Bool_t result = false;
|
||||
|
||||
for (UInt_t i=0; i<param.size(); i++) {
|
||||
if (param[i] != fPreviousParam[i]) {
|
||||
result = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// CalculateField
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
*
|
||||
*/
|
||||
void PNL_PippardFitter::CalculateField(const std::vector<Double_t> ¶m) const
|
||||
{
|
||||
// param: [0] energy, [1] temp, [2] thickness, [3] meanFreePath, [4] xi0, [5] lambdaL, [6] phase
|
||||
|
||||
Int_t pippardFourierPoints = fStartupHandler->GetFourierPoints();
|
||||
|
||||
f_dz = XiP_T(param[4], param[3], param[1])*TMath::TwoPi()/pippardFourierPoints/f_dx; // see lab-book p.137, used for specular reflection boundary conditions (default)
|
||||
|
||||
// check if it is necessary to allocate memory
|
||||
if (fFieldq == 0) {
|
||||
fFieldq = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * pippardFourierPoints);
|
||||
if (fFieldq == 0) {
|
||||
cout << endl << "PPippard::CalculateField(): **ERROR** couldn't allocate memory for fFieldq";
|
||||
cout << endl;
|
||||
return;
|
||||
}
|
||||
}
|
||||
if (fFieldB == 0) {
|
||||
fFieldB = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * pippardFourierPoints);
|
||||
if (fFieldB == 0) {
|
||||
cout << endl << "PPippard::CalculateField(): **ERROR** couldn't allocate memory for fFieldB";
|
||||
cout << endl;
|
||||
return;
|
||||
}
|
||||
}
|
||||
|
||||
// calculate the prefactor of the reduced kernel
|
||||
Double_t xiP = XiP_T(param[4], param[3], param[1]);
|
||||
Double_t preFactor = pow(xiP/(LambdaL_T(param[5], param[1])),2.0)*xiP/param[4];
|
||||
|
||||
|
||||
// calculate the fFieldq vector, which is x/(x^2 + alpha k(x)), with alpha = xiP(T)^3/(lambdaL(T)^2 xiP(0)), and
|
||||
// k(x) = 3/2 [(1+x^2) arctan(x) - x]/x^3, see lab-book p.137
|
||||
Double_t x;
|
||||
fFieldq[0][0] = 0.0;
|
||||
fFieldq[0][1] = 0.0;
|
||||
for (Int_t i=1; i<pippardFourierPoints; i++) {
|
||||
x = i * f_dx;
|
||||
fFieldq[i][0] = x/(pow(x,2.0)+preFactor*(1.5*((1.0+pow(x,2.0))*atan(x)-x)/pow(x,3.0)));
|
||||
fFieldq[i][1] = 0.0;
|
||||
}
|
||||
|
||||
// Fourier transform
|
||||
if (!fPlanPresent) {
|
||||
fPlan = fftw_plan_dft_1d(pippardFourierPoints, fFieldq, fFieldB, FFTW_FORWARD, FFTW_EXHAUSTIVE);
|
||||
fPlanPresent = true;
|
||||
}
|
||||
|
||||
fftw_execute(fPlan);
|
||||
|
||||
// normalize fFieldB
|
||||
Double_t norm = 0.0;
|
||||
fShift=0;
|
||||
for (Int_t i=0; i<pippardFourierPoints/2; i++) {
|
||||
if (fabs(fFieldB[i][1]) > fabs(norm)) {
|
||||
norm = fFieldB[i][1];
|
||||
fShift = i;
|
||||
}
|
||||
}
|
||||
|
||||
for (Int_t i=0; i<pippardFourierPoints; i++) {
|
||||
fFieldB[i][1] /= norm;
|
||||
}
|
||||
|
||||
if (param[2] < pippardFourierPoints/2.0*f_dz) {
|
||||
// B(z) = b(z)+b(D-z)/(1+b(D)) is the B(z) result
|
||||
Int_t idx = (Int_t)(param[2]/f_dz);
|
||||
norm = 1.0 + fFieldB[idx+fShift][1];
|
||||
for (Int_t i=0; i<pippardFourierPoints; i++) {
|
||||
fFieldB[i][0] = 0.0;
|
||||
}
|
||||
for (Int_t i=fShift; i<idx+fShift; i++) {
|
||||
fFieldB[i][0] = (fFieldB[i][1] + fFieldB[idx-i+2*fShift][1])/norm;
|
||||
}
|
||||
for (Int_t i=0; i<pippardFourierPoints; i++) {
|
||||
fFieldB[i][1] = fFieldB[i][0];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// CalculatePolarization
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
*
|
||||
*/
|
||||
void PNL_PippardFitter::CalculatePolarization(const std::vector<Double_t> ¶m) const
|
||||
{
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// DeltaBCS
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
*
|
||||
*/
|
||||
Double_t PNL_PippardFitter::DeltaBCS(const Double_t t) const
|
||||
{
|
||||
Double_t result = 0.0;
|
||||
|
||||
// reduced temperature table
|
||||
Double_t tt[] = {1, 0.98, 0.96, 0.94, 0.92, 0.9, 0.88, 0.86, 0.84, 0.82,
|
||||
0.8, 0.78, 0.76, 0.74, 0.72, 0.7, 0.68, 0.66, 0.64, 0.62,
|
||||
0.6, 0.58, 0.56, 0.54, 0.52, 0.5, 0.48, 0.46, 0.44, 0.42,
|
||||
0.4, 0.38, 0.36, 0.34, 0.32, 0.3, 0.28, 0.26, 0.24, 0.22,
|
||||
0.2, 0.18, 0.16, 0.14};
|
||||
|
||||
// gap table from Muehlschlegel
|
||||
Double_t ee[] = {0, 0.2436, 0.3416, 0.4148, 0.4749, 0.5263, 0.5715, 0.6117,
|
||||
0.648, 0.681, 0.711, 0.7386, 0.764, 0.7874, 0.8089, 0.8288,
|
||||
0.8471, 0.864, 0.8796, 0.8939, 0.907, 0.919, 0.9299, 0.9399,
|
||||
0.9488, 0.9569, 0.9641, 0.9704, 0.976, 0.9809, 0.985, 0.9885,
|
||||
0.9915, 0.9938, 0.9957, 0.9971, 0.9982, 0.9989, 0.9994, 0.9997,
|
||||
0.9999, 1, 1, 1, 1};
|
||||
|
||||
if (t>1.0)
|
||||
result = 0.0;
|
||||
else if (t<0.14)
|
||||
result = 1.0;
|
||||
else {
|
||||
// find correct bin for t
|
||||
UInt_t i=0;
|
||||
for (i=0; i<sizeof(tt); i++) {
|
||||
if (tt[i]<=t) break;
|
||||
}
|
||||
// interpolate linear between 2 bins
|
||||
result = ee[i]-(tt[i]-t)*(ee[i]-ee[i-1])/(tt[i]-tt[i-1]);
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// LambdaL_T
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
*
|
||||
*/
|
||||
Double_t PNL_PippardFitter::LambdaL_T(const Double_t lambdaL, const Double_t t) const
|
||||
{
|
||||
return lambdaL/sqrt(1.0-pow(t,4.0));
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// XiP_T
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p> Approximated xi_P(T). The main approximation is that (lamdaL(T)/lambdaL(0))^2 = 1/(1-t^2). This way
|
||||
* xi_P(T) is close the the BCS xi_BCS(T).
|
||||
*/
|
||||
Double_t PNL_PippardFitter::XiP_T(const Double_t xi0, const Double_t meanFreePath, Double_t t) const
|
||||
{
|
||||
if (t>0.96)
|
||||
t=0.96;
|
||||
|
||||
Double_t J0T = DeltaBCS(t)/(1.0-pow(t,2.0)) * tanh(0.881925 * DeltaBCS(t) / t);
|
||||
|
||||
return xi0*meanFreePath/(meanFreePath*J0T+xi0);
|
||||
}
|
||||
|
||||
21
src/external/Nonlocal/PNL_PippardFitter.h
vendored
21
src/external/Nonlocal/PNL_PippardFitter.h
vendored
@ -32,6 +32,8 @@
|
||||
#ifndef _PNL_PIPPARDFITTER_H_
|
||||
#define _PNL_PIPPARDFITTER_H_
|
||||
|
||||
#include <fftw3.h>
|
||||
|
||||
#include "PUserFcnBase.h"
|
||||
#include "PNL_StartupHandler.h"
|
||||
#include "PNL_RgeHandler.h"
|
||||
@ -48,6 +50,25 @@ class PNL_PippardFitter : public PUserFcnBase
|
||||
PNL_StartupHandler *fStartupHandler;
|
||||
PNL_RgeHandler *fRgeHandler;
|
||||
|
||||
mutable std::vector<Double_t> fPreviousParam;
|
||||
|
||||
Double_t f_dx; // dx = xiPT dq
|
||||
mutable Double_t f_dz; // spatial step size
|
||||
|
||||
mutable Bool_t fPlanPresent;
|
||||
mutable fftw_plan fPlan;
|
||||
mutable fftw_complex *fFieldq; // (xiPT x)/(x^2 + xiPT^2 K(x,T)), x = q xiPT
|
||||
mutable fftw_complex *fFieldB; // field calculated for specular boundary conditions
|
||||
mutable Int_t fShift; // shift needed to pick up fFieldB at the maximum for B->0
|
||||
|
||||
virtual Bool_t NewParameters(const std::vector<Double_t> ¶m) const;
|
||||
virtual void CalculateField(const std::vector<Double_t> ¶m) const;
|
||||
virtual void CalculatePolarization(const std::vector<Double_t> ¶m) const;
|
||||
|
||||
virtual Double_t DeltaBCS(const Double_t t) const;
|
||||
virtual Double_t LambdaL_T(const Double_t lambdaL, const Double_t t) const;
|
||||
virtual Double_t XiP_T(const Double_t xi0, const Double_t meanFreePath, Double_t t) const;
|
||||
|
||||
ClassDef(PNL_PippardFitter, 1)
|
||||
};
|
||||
|
||||
|
||||
3
src/external/Nonlocal/PNL_RgeHandler.cpp
vendored
3
src/external/Nonlocal/PNL_RgeHandler.cpp
vendored
@ -147,6 +147,9 @@ Bool_t PNL_RgeHandler::LoadRgeData(const PStringVector &rgeDataPathList)
|
||||
// open rge-file for reading
|
||||
fin.open(rgeDataPathList[i].Data(), iostream::in);
|
||||
if (!fin.is_open()) {
|
||||
cout << endl << "PNL_RgeHandler::LoadRgeData **ERROR**";
|
||||
cout << endl << " Could not open file " << rgeDataPathList[i].Data();
|
||||
cout << endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
4
src/external/Nonlocal/PNL_StartupHandler.cpp
vendored
4
src/external/Nonlocal/PNL_StartupHandler.cpp
vendored
@ -29,6 +29,8 @@
|
||||
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
|
||||
***************************************************************************/
|
||||
|
||||
#include <cmath>
|
||||
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
using namespace std;
|
||||
@ -180,7 +182,7 @@ void PNL_StartupHandler::OnCharacters(const char *str)
|
||||
fIsValid = false;
|
||||
}
|
||||
tstr = fTrimSpDataPath;
|
||||
sprintf(sstr, "%03d", (int)(dval*10.0));
|
||||
sprintf(sstr, "%03d", (int)(round(dval*10.0)));
|
||||
tstr += sstr;
|
||||
tstr += ".rge";
|
||||
fTrimSpDataPathList.push_back(tstr);
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user