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

newly added

Browse Source
This commit is contained in:
nemu 2009-06-26 13:32:31 +00:00
parent bad93a0db1
commit d3c5c7841a
6 changed files with 953 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

99
src/tests/nonlocal/Makefile Normal file
View File

@ -0,0 +1,99 @@
#---------------------------------------------------
# Makefile
#
# Author: Andreas Suter
# e-mail: andreas.suter@psi.ch
#
#---------------------------------------------------
#---------------------------------------------------
# get compilation and library flags from root-config
ROOTCFLAGS = $(shell $(ROOTSYS)/bin/root-config --cflags)
ROOTLIBS = $(shell $(ROOTSYS)/bin/root-config --libs)
ROOTGLIBS = $(shell $(ROOTSYS)/bin/root-config --glibs)
#---------------------------------------------------
# depending on the architecture, choose the compiler,
# linker, and the flags to use
#
ARCH = $(shell $(ROOTSYS)/bin/root-config --arch)
ifeq ($(ARCH),linux)
OS = LINUX
endif
ifeq ($(ARCH),linuxx8664gcc)
OS = LINUX
endif
ifeq ($(ARCH),win32gcc)
OS = WIN32GCC
endif
ifeq ($(ARCH),macosx)
OS = DARWIN
endif
# -- Linux
ifeq ($(OS),LINUX)
CXX = g++
CXXFLAGS = -O3 -Wall -fPIC
PMUSRPATH = ./include
MNPATH = $(ROOTSYS)/include
GSLPATH = /usr/include/gsl
EIGEN2PATH = /usr/local/include/eigen2
INCLUDES = -I$(PMUSRPATH) -I$(MNPATH) -I$(GSLPATH) -I$(EIGEN2PATH)
LD = g++
LDFLAGS = -O
INSTALLPATH = $(ROOTSYS)/bin
EXEC = nonlocal
SUFFIX =
endif
# the output from the root-config script:
CXXFLAGS += $(ROOTCFLAGS)
LDFLAGS +=
# the ROOT libraries (G = graphic)
LIBS = $(ROOTLIBS) -lXMLParser
GLIBS = $(ROOTGLIBS) -lXMLParser
# PSI libs
PSILIBS = -L$(ROOTSYS)/lib -lTLemRunHeader -lPMusr
# Minuit2 lib
MNLIB = -L$(ROOTSYS)/lib -lMinuit2
# MathMore lib
MMLIB = -L$(ROOTSYS)/lib -lMathMore
# GSL lib
GSLLIB = -L/usr/lib -lgsl
# some definitions: headers, sources, objects,...
OBJS =
OBJS += nonlocal.o PPippard.o
# make the executable:
#
all: $(EXEC)
nonlocal$(SUFFIX): $(OBJS)
@echo "---> Building $@ ..."
/bin/rm -f $@
$(LD) $< -o $@ PPippard.o $(LDFLAGS) $(GLIBS) $(PSILIBS) $(MNLIB) $(MMLIB) $(GSLLIB)
@echo "done"
# clean up: remove all object file (and core files)
# semicolon needed to tell make there is no source
# for this target!
#
clean:; @rm -f $(OBJS)
@echo "---> removing $(OBJS)"
#
$(OBJS): %.o: %.cpp
$(CXX) $(INCLUDES) $(CXXFLAGS) -c $<
install: all
cp -fvp $(EXEC) $(INSTALLPATH)
cp -fvp musrfit_startup.xml $(INSTALLPATH)
cp -fvp external/scripts/msr2data $(INSTALLPATH)
chmod 755 $(INSTALLPATH)/msr2data

188
src/tests/nonlocal/NonlocalDiffuseFunc_v_w.nb Normal file
View File

@ -0,0 +1,188 @@
(************** Content-type: application/mathematica **************
CreatedBy='Mathematica 5.2'
Mathematica-Compatible Notebook
This notebook can be used with any Mathematica-compatible
application, such as Mathematica, MathReader or Publicon. The data
for the notebook starts with the line containing stars above.
To get the notebook into a Mathematica-compatible application, do
one of the following:
* Save the data starting with the line of stars above into a file
with a name ending in .nb, then open the file inside the
application;
* Copy the data starting with the line of stars above to the
clipboard, then use the Paste menu command inside the application.
Data for notebooks contains only printable 7-bit ASCII and can be
sent directly in email or through ftp in text mode. Newlines can be
CR, LF or CRLF (Unix, Macintosh or MS-DOS style).
NOTE: If you modify the data for this notebook not in a Mathematica-
compatible application, you must delete the line below containing
the word CacheID, otherwise Mathematica-compatible applications may
try to use invalid cache data.
For more information on notebooks and Mathematica-compatible
applications, contact Wolfram Research:
web: http://www.wolfram.com
email: info@wolfram.com
phone: +1-217-398-0700 (U.S.)
Notebook reader applications are available free of charge from
Wolfram Research.
*******************************************************************)
(*CacheID: 232*)
(*NotebookFileLineBreakTest
NotebookFileLineBreakTest*)
(*NotebookOptionsPosition[ 3636, 129]*)
(*NotebookOutlinePosition[ 4267, 151]*)
(* CellTagsIndexPosition[ 4223, 147]*)
(*WindowFrame->Normal*)
Notebook[{
Cell[BoxData[
\(v\ = \
2/3 + 1/6 Exp[\(-s\)] \((s \((s - 1)\) - 4)\) -
1/6 s \((s\^2 - 6)\) Gamma[0, s]\)], "Input"],
Cell[CellGroupData[{
Cell[BoxData[
\(w\ = \ \((8 s - 3)\)/12 +
1/24 Exp[\(-s\)] \((s\^3 - s\^2 - 10 s + 6)\) -
1/24 \( s\^2\) \((s\^2 - 12)\) Gamma[0, s]\)], "Input"],
Cell[BoxData[
\(1\/12\ \((\(-3\) + 8\ s)\) +
1\/24\ \[ExponentialE]\^\(-s\)\ \((6 - 10\ s - s\^2 + s\^3)\) -
1\/24\ s\^2\ \((\(-12\) + s\^2)\)\ Gamma[0, s]\)], "Output"]
}, Open ]],
Cell[CellGroupData[{
Cell[BoxData[
\(Limit[w, s \[Rule] 0]\)], "Input"],
Cell[BoxData[
\(0\)], "Output"]
}, Open ]],
Cell[CellGroupData[{
Cell[BoxData[
\(Series[v, {s, 0, 2}]\)], "Input"],
Cell[BoxData[
InterpretationBox[
RowBox[{\(\((1\/2 - EulerGamma - Log[s])\)\ s\), "+", \(s\^2\), "+",
InterpretationBox[\(O[s]\^3\),
SeriesData[ s, 0, {}, 1, 3, 1],
Editable->False]}],
SeriesData[ s, 0, {
Plus[
Rational[ 1, 2],
Times[ -1, EulerGamma],
Times[ -1,
Log[ s]]], 1}, 1, 3, 1],
Editable->False]], "Output"]
}, Open ]],
Cell[CellGroupData[{
Cell[BoxData[
\(N[EulerGamma]\)], "Input"],
Cell[BoxData[
\(0.5772156649015329`\)], "Output"]
}, Open ]],
Cell[CellGroupData[{
Cell[BoxData[
\(Series[w, {s, 0, 2}]\)], "Input"],
Cell[BoxData[
InterpretationBox[
RowBox[{\(\((1\/2 - EulerGamma\/2 - Log[s]\/2)\)\ s\^2\), "+",
InterpretationBox[\(O[s]\^3\),
SeriesData[ s, 0, {}, 2, 3, 1],
Editable->False]}],
SeriesData[ s, 0, {
Plus[
Rational[ 1, 2],
Times[
Rational[ -1, 2], EulerGamma],
Times[
Rational[ -1, 2],
Log[ s]]]}, 2, 3, 1],
Editable->False]], "Output"]
}, Open ]]
},
FrontEndVersion->"5.2 for X",
ScreenRectangle->{{0, 1280}, {0, 1024}},
WindowSize->{520, 600},
WindowMargins->{{150, Automatic}, {Automatic, 52}}
]
(*******************************************************************
Cached data follows. If you edit this Notebook file directly, not
using Mathematica, you must remove the line containing CacheID at
the top of the file. The cache data will then be recreated when
you save this file from within Mathematica.
*******************************************************************)
(*CellTagsOutline
CellTagsIndex->{}
*)
(*CellTagsIndex
CellTagsIndex->{}
*)
(*NotebookFileOutline
Notebook[{
Cell[1754, 51, 139, 3, 47, "Input"],
Cell[CellGroupData[{
Cell[1918, 58, 171, 3, 51, "Input"],
Cell[2092, 63, 187, 3, 80, "Output"]
}, Open ]],
Cell[CellGroupData[{
Cell[2316, 71, 54, 1, 27, "Input"],
Cell[2373, 74, 35, 1, 27, "Output"]
}, Open ]],
Cell[CellGroupData[{
Cell[2445, 80, 53, 1, 27, "Input"],
Cell[2501, 83, 421, 12, 44, "Output"]
}, Open ]],
Cell[CellGroupData[{
Cell[2959, 100, 46, 1, 27, "Input"],
Cell[3008, 103, 53, 1, 27, "Output"]
}, Open ]],
Cell[CellGroupData[{
Cell[3098, 109, 53, 1, 27, "Input"],
Cell[3154, 112, 466, 14, 44, "Output"]
}, Open ]]
}
]
*)
(*******************************************************************
End of Mathematica Notebook file.
*******************************************************************)

480
src/tests/nonlocal/PPippard.cpp Normal file
View File

@ -0,0 +1,480 @@
// -----------------------------------------------------------------------
// Author: Andreas Suter
// $Id$
// -----------------------------------------------------------------------
#include <cstdio>
#include <cmath>
#include <iostream>
using namespace std;
#include <gsl_sf_gamma.h>
#include <TMath.h>
#include "PPippard.h"
//-----------------------------------------------------------------------------------------------------------
/**
*
*/
PPippard::PPippard(Double_t temp, Double_t lambdaL, Double_t xi0, Double_t meanFreePath, Double_t filmThickness, Bool_t specular) :
fTemp(temp), fLambdaL(lambdaL), fXi0(xi0), fMeanFreePath(meanFreePath), fFilmThickness(filmThickness), fSpecular(specular)
{
fPlanPresent = false;
fFieldq = 0;
fFieldB = 0;
fSecondDerivativeMatrix = 0;
fKernelMatrix = 0;
fBoundaryCondition = 0;
fFieldDiffuse = 0;
f_dx = 0.02;
f_dz = XiP_T(fTemp)*TMath::TwoPi()/PippardFourierPoints/f_dx; // see lab-book p.137, used for specular reflection boundary conditions (default)
fShift = 0;
}
//-----------------------------------------------------------------------------------------------------------
/**
*
*/
PPippard::~PPippard()
{
if (fPlanPresent) {
fftw_destroy_plan(fPlan);
}
if (fFieldq) {
fftw_free(fFieldq);
fFieldq = 0;
}
if (fFieldB) {
fftw_free(fFieldq);
fFieldB = 0;
}
if (fSecondDerivativeMatrix) {
delete fSecondDerivativeMatrix;
fSecondDerivativeMatrix = 0;
}
if (fKernelMatrix) {
delete fKernelMatrix;
fKernelMatrix = 0;
}
if (fBoundaryCondition) {
delete fBoundaryCondition;
fBoundaryCondition = 0;
}
if (fFieldDiffuse) {
delete fFieldDiffuse;
fFieldDiffuse = 0;
}
}
//-----------------------------------------------------------------------------------------------------------
/**
*
*/
Double_t PPippard::GetMagneticField(const Double_t z) const
{
Double_t result = -1.0;
if (fSpecular) {
if (fFieldB == 0)
return -1.0;
if (z < 0.0)
return 1.0;
if (z > f_dz*PippardFourierPoints/2.0)
return 0.0;
Int_t bin = (Int_t)(z/f_dz);
result = fFieldB[bin+fShift][1];
} else { // diffuse
if (fFieldDiffuse == 0)
return -1.0;
if (z < 0.0)
return 1.0;
if (z > PippardDiffusePoints * f_dz * XiP_T(fTemp))
return 0.0;
Int_t bin = (Int_t)(z/(f_dz*XiP_T(fTemp)));
result = (*fFieldDiffuse)(bin);
}
return result;
}
//-----------------------------------------------------------------------------------------------------------
/**
*
*/
Double_t PPippard::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;
}
//-----------------------------------------------------------------------------------------------------------
/**
*
*/
Double_t PPippard::LambdaL_T(const Double_t t) const
{
return fLambdaL/sqrt(1.0-pow(t,4.0));
}
//-----------------------------------------------------------------------------------------------------------
/**
* <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 PPippard::XiP_T(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 fXi0*fMeanFreePath/(fMeanFreePath*J0T+fXi0);
}
//-----------------------------------------------------------------------------------------------------------
/**
*
*/
void PPippard::CalculateField()
{
// calculate the field
if (fSpecular)
CalculateFieldSpecular();
else
CalculateFieldDiffuse();
}
//-----------------------------------------------------------------------------------------------------------
/**
*
*/
void PPippard::CalculateFieldSpecular()
{
// 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(fTemp);
Double_t preFactor = pow(xiP/(LambdaL_T(fTemp)),2.0)*xiP/fXi0;
// 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_ESTIMATE);
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;
}
}
cout << endl << "fShift = " << fShift;
for (Int_t i=0; i<PippardFourierPoints; i++) {
fFieldB[i][1] /= norm;
}
if (fFilmThickness < 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)(fFilmThickness/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];
}
}
Double_t integral = 0.0;
for (Int_t i=fShift; i<PippardFourierPoints/2; i++)
integral += fFieldB[i][1];
cout << endl << "specular Integral = " << integral*f_dz;
}
//-----------------------------------------------------------------------------------------------------------
/**
*
*/
void PPippard::CalculateFieldDiffuse()
{
f_dz = 5.0/XiP_T(fTemp);
Double_t invL = 1/f_dz;
Double_t ampl = 1.0/pow(f_dz,2.0)/(3.0/4.0*pow(XiP_T(fTemp),3.0)/(fXi0*pow(LambdaL_T(fTemp),2.0)));
cout << endl << ">> 1/alpha = " << 1.0/(3.0/4.0*pow(XiP_T(fTemp),3.0)/(fXi0*pow(LambdaL_T(fTemp),2.0)));
cout << endl << ">> 1/l^2 = " << 1.0/pow(f_dz,2.0);
cout << endl << ">> ampl = " << ampl << endl;
// 2nd derivative matrix
if (fSecondDerivativeMatrix == 0) { // first time call, hence generate the 2nd derivative matrix
fSecondDerivativeMatrix = new MatrixXd(PippardDiffusePoints+1, PippardDiffusePoints+1);
fSecondDerivativeMatrix->setZero(PippardDiffusePoints+1, PippardDiffusePoints+1);
for (Int_t i=1; i<PippardDiffusePoints; i++) {
(*fSecondDerivativeMatrix)(i,i-1) = ampl;
(*fSecondDerivativeMatrix)(i,i) = -2.0*ampl;
(*fSecondDerivativeMatrix)(i,i+1) = ampl;
}
}
//cout << endl << "fSecondDerivativeMatrix = \n" << *fSecondDerivativeMatrix << endl;
// kernel matrix
if (fKernelMatrix == 0) { // first time call, hence generate the kernel matrix
fKernelMatrix = new MatrixXd(PippardDiffusePoints+1, PippardDiffusePoints+1);
fKernelMatrix->setZero(PippardDiffusePoints+1, PippardDiffusePoints+1);
// 1st line (dealing with boundary conditions)
(*fKernelMatrix)(0,0) = -1.5*invL;
(*fKernelMatrix)(0,1) = 2.0*invL;
(*fKernelMatrix)(0,2) = -0.5*invL;
// Nth line (dealing with boundary conditions)
(*fKernelMatrix)(PippardDiffusePoints,PippardDiffusePoints-2) = 0.5*invL;
(*fKernelMatrix)(PippardDiffusePoints,PippardDiffusePoints-1) = -2.0*invL;
(*fKernelMatrix)(PippardDiffusePoints,PippardDiffusePoints) = 1.5*invL;
// the real kernel
for (Int_t i=1; i<PippardDiffusePoints; i++) {
(*fKernelMatrix)(i,0)=Calc_a(i);
(*fKernelMatrix)(i,PippardDiffusePoints)=Calc_b(i);
for (Int_t j=1; j<PippardDiffusePoints; j++) {
(*fKernelMatrix)(i,j) = Calc_c(i,j);
}
}
}
//cout << endl << "fKernelMatrix = \n" << *fKernelMatrix << endl;
// boundary condition vector
if (fBoundaryCondition == 0) {
fBoundaryCondition = new VectorXd(PippardDiffusePoints+1);
fBoundaryCondition->setZero(PippardDiffusePoints+1);
(*fBoundaryCondition)(0) = 1.0;
}
//cout << endl << "fBoundaryCondition = " << *fBoundaryCondition << endl;
if (fFieldDiffuse == 0) {
fFieldDiffuse = new VectorXd(PippardDiffusePoints+1);
fFieldDiffuse->setZero(PippardDiffusePoints+1);
}
// solve equation
*fSecondDerivativeMatrix = (*fSecondDerivativeMatrix)-(*fKernelMatrix);
fSecondDerivativeMatrix->lu().solve(*fBoundaryCondition, fFieldDiffuse);
// normalize field
Double_t norm = 0.0;
for (Int_t i=0; i<PippardDiffusePoints+1; i++)
if (norm < (*fFieldDiffuse)(i))
norm = (*fFieldDiffuse)(i);
for (Int_t i=0; i<PippardDiffusePoints+1; i++)
(*fFieldDiffuse)(i) /= norm;
Double_t integral = 0.0;
for (Int_t i=0; i<PippardDiffusePoints+1; i++)
integral += (*fFieldDiffuse)(i);
cout << endl << "Diffuse Integral = " << integral*f_dz*XiP_T(fTemp);
}
//-----------------------------------------------------------------------------------------------------------
/**
*
*/
void PPippard::SaveField(const char *fileName)
{
FILE *fp;
fp = fopen(fileName, "w");
if (fp == NULL) {
cout << endl << "Coudln't open " << fileName << " for writting, sorry ...";
cout << endl << endl;
return;
}
// write header
fprintf(fp, "%% Header ------------------------------------\n");
fprintf(fp, "%% Parameters:\n");
fprintf(fp, "%% Reduced Temperature = %lf\n", fTemp);
fprintf(fp, "%% LambdaL(0) = %lf, LambdaL(t) = %lf\n", fLambdaL, LambdaL_T(fTemp));
fprintf(fp, "%% xiP(0) = %lf, xiP(t) = %lf\n", fXi0, XiP_T(fTemp));
fprintf(fp, "%% Mean Free Path = %lf\n", fMeanFreePath);
if (fSpecular)
fprintf(fp, "%% Boundary Conditions: Specular\n");
else
fprintf(fp, "%% Boundary Conditions: Diffuse\n");
fprintf(fp, "%%\n");
// write data
fprintf(fp, "%% Data --------------------------------------\n");
fprintf(fp, "%% z (nm), B/B_0 \n");
if (fSpecular) {
for (Int_t i=0; i<PippardFourierPoints/2; i++) {
fprintf(fp, "%lf, %lf\n", f_dz*(Double_t)i, fFieldB[i+fShift][1]);
}
} else {
for (Int_t i=0; i<PippardDiffusePoints; i++) {
fprintf(fp, "%lf, %lf\n", f_dz * XiP_T(fTemp) * (Double_t)i, (*fFieldDiffuse)(i));
}
}
fclose(fp);
}
//-----------------------------------------------------------------------------------------------------------
/**
*
*/
Double_t PPippard::Calc_v(const Double_t s) const
{
if (s == 0.0)
return 0.0;
Double_t s2 = pow(s,2.0);
Double_t ss = fabs(s);
Double_t v;
if (ss < 0.001) {
v = (-0.0772157-log(ss))*ss+s2; // series expansion in s up to 2nd order
} else {
v = 1/6.0*(exp(-ss)*(s2-ss-4.0) + 4.0 - ss*(s2-6.0) * gsl_sf_gamma_inc(0.0, ss));
}
if (s < 0)
v = -v;
return v;
}
//-----------------------------------------------------------------------------------------------------------
/**
*
*/
Double_t PPippard::Calc_w(const Double_t s) const
{
if (s == 0.0)
return 0.0;
Double_t s2 = pow(s,2.0);
Double_t ss = fabs(s);
Double_t w;
if (ss < 0.001) {
w = (0.211392 - 0.5*log(ss))*s2; // series expansion in s up to 2nd order
} else {
w = 1/24.0*(exp(-ss)*(6.0-10.0*ss-s2+ss*s2)+16.0*ss-6.0-s2*(s2-12.0)*gsl_sf_gamma_inc(0.0, ss));
}
return w;
}
//-----------------------------------------------------------------------------------------------------------
/**
*
*/
Double_t PPippard::Calc_a(const Int_t i) const
{
return -Calc_v(-i*f_dz) + (Calc_w((1-i)*f_dz)-Calc_w(-i*f_dz))/f_dz;
}
//-----------------------------------------------------------------------------------------------------------
/**
*
*/
Double_t PPippard::Calc_b(const Int_t i) const
{
return Calc_v((PippardDiffusePoints-i)*f_dz) - (Calc_w((PippardDiffusePoints-i)*f_dz)-Calc_w((PippardDiffusePoints-1-i)*f_dz))/f_dz;
}
//-----------------------------------------------------------------------------------------------------------
/**
*
*/
Double_t PPippard::Calc_c(const Int_t i, const Int_t j) const
{
return (Calc_w((i+1-j)*f_dz)-2.0*Calc_w((i-j)*f_dz)+Calc_w((i-1-j)*f_dz))/f_dz;
}

77
src/tests/nonlocal/PPippard.h Normal file
View File

@ -0,0 +1,77 @@
// -----------------------------------------------------------------------
// Author: Andreas Suter
// $Id$
// -----------------------------------------------------------------------
#ifndef _PPIPPARD_H_
#define _PPIPPARD_H_
#include <fftw3.h>
#include <Rtypes.h>
#include <Eigen/Core>
#include <Eigen/LU>
using namespace Eigen;
const Int_t PippardFourierPoints = 200000;
const Int_t PippardDiffusePoints = 500;
class PPippard
{
public:
PPippard(Double_t temp, Double_t lambdaL, Double_t xi0, Double_t meanFreePath, Double_t filmThickness, Bool_t specular = true);
virtual ~PPippard();
virtual void SetTemp(Double_t temp) { fTemp = temp; }
virtual void SetLambdaL(Double_t lambdaL) { fLambdaL = lambdaL; }
virtual void SetXi0(Double_t xi0) { fXi0 = xi0; }
virtual void SetMeanFreePath(Double_t meanFreePath) { fMeanFreePath = meanFreePath; }
virtual void SetFilmThickness(Double_t thickness) { fFilmThickness = thickness; }
virtual void SetSpecular(Bool_t specular) { fSpecular = specular; }
virtual void CalculateField();
virtual Double_t GetMagneticField(const Double_t x) const;
virtual void SaveField(const char *fileName);
private:
Double_t fTemp; // reduced temperature, i.e. t = T/T_c
Double_t fLambdaL; // lambdaL in (nm)
Double_t fXi0; // xi0 in (nm)
Double_t fMeanFreePath; // mean free path in (nm)
Double_t fFilmThickness; // film thickness in (nm)
Bool_t fSpecular; // = 1 -> specular, 0 -> diffuse
Double_t f_dx; // dx = xiPT dq
Double_t f_dz; // spatial step size
bool fPlanPresent;
fftw_plan fPlan;
fftw_complex *fFieldq; // (xiPT x)/(x^2 + xiPT^2 K(x,T)), x = q xiPT
fftw_complex *fFieldB; // field calculated for specular boundary conditions
Int_t fShift; // shift needed to pick up fFieldB at the maximum for B->0
MatrixXd *fSecondDerivativeMatrix; // 2nd derivative matrix
MatrixXd *fKernelMatrix; // kernel matrix
VectorXd *fBoundaryCondition; // boundary condition vector
VectorXd *fFieldDiffuse; // resulting B(z)/B(0) field
virtual Double_t DeltaBCS(const Double_t t) const;
virtual Double_t LambdaL_T(const Double_t t) const;
virtual Double_t XiP_T(Double_t t) const;
virtual void CalculateFieldSpecular();
virtual void CalculateFieldDiffuse();
virtual Double_t Calc_v(const Double_t s) const; // see 'A.Suter, Memorandum June 17, 2004' Eq.(13)
virtual Double_t Calc_w(const Double_t s) const; // see 'A.Suter, Memorandum June 17, 2004' Eq.(14)
virtual Double_t Calc_a(const Int_t i) const; // see 'A.Suter, Memorandum June 17, 2004' Eq.(17)
virtual Double_t Calc_b(const Int_t i) const; // see 'A.Suter, Memorandum June 17, 2004' Eq.(17)
virtual Double_t Calc_c(const Int_t i, const Int_t j) const; // see 'A.Suter, Memorandum June 17, 2004' Eq.(17)
};
#endif // _PPIPPARD_H_

82
src/tests/nonlocal/nonlocal.cpp Normal file
View File

@ -0,0 +1,82 @@
// -----------------------------------------------------------------------
// Author: Andreas Suter
// $Id$
// -----------------------------------------------------------------------
#include <cstdlib>
#include <iostream>
using namespace std;
#include "PPippard.h"
void syntax()
{
cout << endl << "usage: nonlocal temp lambdaL xi0 meanFreePath filmThickness specular [<fileName>]";
cout << endl << " temp : reduced temperature, i.e. t = T/T_c";
cout << endl << " specular: 1 = specular, 0 = diffuse";
cout << endl << " all lengths given in (nm)";
cout << endl << " if <fileName> is given, the field profile will be saved";
cout << endl << endl;
}
int main(int argc, char *argv[])
{
if ((argc < 7) || (argc>8)){
syntax();
return 0;
}
Double_t temp = strtod(argv[1], (char **)NULL);
Double_t lambdaL = strtod(argv[2], (char **)NULL);
Double_t xi0 = strtod(argv[3], (char **)NULL);
Double_t meanFreePath = strtod(argv[4], (char **)NULL);
Double_t filmThickness = strtod(argv[5], (char **)NULL);
char fln[1024];
if (argc == 8) {
strncpy(fln, argv[7], sizeof(fln));
} else {
strncpy(fln, "", sizeof(fln));
}
Bool_t specular;
if (!strcmp(argv[6], "1"))
specular = true;
else if (!strcmp(argv[6], "0"))
specular = false;
else {
syntax();
return 0;
}
cout << endl << ">> temp = " << temp;
cout << endl << ">> lambdaL = " << lambdaL;
cout << endl << ">> xi0 = " << xi0;
cout << endl << ">> meanFreePath = " << meanFreePath;
cout << endl << ">> filmThickness = " << filmThickness;
if (specular)
cout << endl << ">> specular = true";
else
cout << endl << ">> specular = false";
cout << endl;
PPippard *pippard = new PPippard(temp, lambdaL, xi0, meanFreePath, filmThickness, specular);
pippard->CalculateField();
if (strlen(fln) > 0)
pippard->SaveField(fln);
cout << endl << ">> magnetic field = " << pippard->GetMagneticField(0.0);
cout << endl;
if (pippard) {
delete pippard;
pippard = 0;
}
return 1;
}

27
src/tests/nonlocal/test.C Normal file
View File

@ -0,0 +1,27 @@
{
// 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};
TF1 *userFunc = new TF1("User Func", "[0]*pow(1.0-pow(x,[1]),[2])", 0.0, 1.0);
userFunc->SetParName(0, "A");
userFunc->SetParName(1, "p1");
userFunc->SetParName(2, "p2");
userFunc->SetParameters(1.0, 2.0, 0.5);
TGraph *gr = new TGraph(45, tt, ee);
gr->Draw("ap");
}