changed to file input handling
This commit is contained in:
nemu
21
src/tests/nonlocal/In_37_T2P83_E14P1.inp
Normal file
21
src/tests/nonlocal/In_37_T2P83_E14P1.inp
Normal file
@ -0,0 +1,21 @@
|
||||
% input file for nonlocal
|
||||
% --------------------------------
|
||||
% input parameters:
|
||||
% reduced temperature
|
||||
% lambda London
|
||||
% xi0
|
||||
% mean free path
|
||||
% film thickness
|
||||
% boundary conditions
|
||||
% [rge input file name]
|
||||
% [output file name]
|
||||
|
||||
reducedTemp = 0.8287
|
||||
lambdaL = 25.0
|
||||
xi0 = 380.0
|
||||
meanFreePath = 12000.0
|
||||
filmThickness = 500000.0
|
||||
specular = 1
|
||||
rgeFileName = /afs/psi.ch/project/nemu/analysis/2009/Nonlocal/trimsp/InSne141.rge
|
||||
outputFileName = In_37_T2P83_E14P1.dat
|
||||
|
||||
@ -91,9 +91,3 @@ clean:; @rm -f $(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
|
||||
|
||||
@ -19,8 +19,7 @@ using namespace std;
|
||||
/**
|
||||
*
|
||||
*/
|
||||
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)
|
||||
PPippard::PPippard(const PippardParams ¶ms) : fParams(params)
|
||||
{
|
||||
fPlanPresent = false;
|
||||
fFieldq = 0;
|
||||
@ -32,7 +31,7 @@ PPippard::PPippard(Double_t temp, Double_t lambdaL, Double_t xi0, Double_t meanF
|
||||
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)
|
||||
f_dz = XiP_T(fParams.t)*TMath::TwoPi()/PippardFourierPoints/f_dx; // see lab-book p.137, used for specular reflection boundary conditions (default)
|
||||
fShift = 0;
|
||||
}
|
||||
|
||||
@ -80,7 +79,7 @@ Double_t PPippard::GetMagneticField(const Double_t z) const
|
||||
{
|
||||
Double_t result = -1.0;
|
||||
|
||||
if (fSpecular) {
|
||||
if (fParams.specular) {
|
||||
|
||||
if (fFieldB == 0)
|
||||
return -1.0;
|
||||
@ -103,10 +102,10 @@ Double_t PPippard::GetMagneticField(const Double_t z) const
|
||||
if (z < 0.0)
|
||||
return 1.0;
|
||||
|
||||
if (z > PippardDiffusePoints * f_dz * XiP_T(fTemp))
|
||||
if (z > PippardDiffusePoints * f_dz * XiP_T(fParams.xi0))
|
||||
return 0.0;
|
||||
|
||||
Int_t bin = (Int_t)(z/(f_dz*XiP_T(fTemp)));
|
||||
Int_t bin = (Int_t)(z/(f_dz*XiP_T(fParams.xi0)));
|
||||
|
||||
result = (*fFieldDiffuse)(bin);
|
||||
}
|
||||
@ -160,7 +159,7 @@ Double_t PPippard::DeltaBCS(const Double_t t) const
|
||||
*/
|
||||
Double_t PPippard::LambdaL_T(const Double_t t) const
|
||||
{
|
||||
return fLambdaL/sqrt(1.0-pow(t,4.0));
|
||||
return fParams.lambdaL/sqrt(1.0-pow(t,4.0));
|
||||
}
|
||||
|
||||
//-----------------------------------------------------------------------------------------------------------
|
||||
@ -175,7 +174,7 @@ Double_t PPippard::XiP_T(Double_t t) const
|
||||
|
||||
Double_t J0T = DeltaBCS(t)/(1.0-pow(t,2.0)) * tanh(0.881925 * DeltaBCS(t) / t);
|
||||
|
||||
return fXi0*fMeanFreePath/(fMeanFreePath*J0T+fXi0);
|
||||
return fParams.xi0*fParams.meanFreePath/(fParams.meanFreePath*J0T+fParams.xi0);
|
||||
}
|
||||
|
||||
//-----------------------------------------------------------------------------------------------------------
|
||||
@ -185,7 +184,7 @@ Double_t PPippard::XiP_T(Double_t t) const
|
||||
void PPippard::CalculateField()
|
||||
{
|
||||
// calculate the field
|
||||
if (fSpecular)
|
||||
if (fParams.specular)
|
||||
CalculateFieldSpecular();
|
||||
else
|
||||
CalculateFieldDiffuse();
|
||||
@ -216,8 +215,8 @@ void PPippard::CalculateFieldSpecular()
|
||||
}
|
||||
|
||||
// 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;
|
||||
Double_t xiP = XiP_T(fParams.t);
|
||||
Double_t preFactor = pow(xiP/(LambdaL_T(fParams.t)),2.0)*xiP/fParams.xi0;
|
||||
|
||||
|
||||
// calculate the fFieldq vector, which is x/(x^2 + alpha k(x)), with alpha = xiP(T)^3/(lambdaL(T)^2 xiP(0)), and
|
||||
@ -255,9 +254,9 @@ cout << endl << "fShift = " << fShift;
|
||||
fFieldB[i][1] /= norm;
|
||||
}
|
||||
|
||||
if (fFilmThickness < PippardFourierPoints/2.0*f_dz) {
|
||||
if (fParams.filmThickness < 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);
|
||||
Int_t idx = (Int_t)(fParams.filmThickness/f_dz);
|
||||
norm = 1.0 + fFieldB[idx+fShift][1];
|
||||
for (Int_t i=0; i<PippardFourierPoints; i++) {
|
||||
fFieldB[i][0] = 0.0;
|
||||
@ -273,7 +272,8 @@ cout << endl << "fShift = " << fShift;
|
||||
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;
|
||||
cout << endl << ">> specular Integral = " << integral*f_dz;
|
||||
fParams.specularIntegral = integral*f_dz;
|
||||
}
|
||||
|
||||
//-----------------------------------------------------------------------------------------------------------
|
||||
@ -282,12 +282,12 @@ cout << endl << "fShift = " << fShift;
|
||||
*/
|
||||
void PPippard::CalculateFieldDiffuse()
|
||||
{
|
||||
f_dz = 5.0/XiP_T(fTemp);
|
||||
f_dz = 5.0/XiP_T(fParams.t);
|
||||
|
||||
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)));
|
||||
Double_t ampl = 1.0/pow(f_dz,2.0)/(3.0/4.0*pow(XiP_T(fParams.t),3.0)/(fParams.xi0*pow(LambdaL_T(fParams.t),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/alpha = " << 1.0/(3.0/4.0*pow(XiP_T(fParams.t),3.0)/(fParams.xi0*pow(LambdaL_T(fParams.t),2.0)));
|
||||
cout << endl << ">> 1/l^2 = " << 1.0/pow(f_dz,2.0);
|
||||
cout << endl << ">> ampl = " << ampl << endl;
|
||||
|
||||
@ -358,20 +358,20 @@ cout << endl << ">> ampl = " << ampl << endl;
|
||||
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);
|
||||
cout << endl << "Diffuse Integral = " << integral*f_dz*XiP_T(fParams.t);
|
||||
}
|
||||
|
||||
//-----------------------------------------------------------------------------------------------------------
|
||||
/**
|
||||
*
|
||||
*/
|
||||
void PPippard::SaveField(const char *fileName)
|
||||
void PPippard::SaveField()
|
||||
{
|
||||
FILE *fp;
|
||||
|
||||
fp = fopen(fileName, "w");
|
||||
fp = fopen(fParams.outputFileName.Data(), "w");
|
||||
if (fp == NULL) {
|
||||
cout << endl << "Coudln't open " << fileName << " for writting, sorry ...";
|
||||
cout << endl << "Coudln't open " << fParams.outputFileName.Data() << " for writting, sorry ...";
|
||||
cout << endl << endl;
|
||||
return;
|
||||
}
|
||||
@ -379,26 +379,33 @@ void PPippard::SaveField(const char *fileName)
|
||||
// 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, "%% Reduced Temperature = %lf\n", fParams.t);
|
||||
fprintf(fp, "%% LambdaL(0) = %lf, LambdaL(t) = %lf\n", fParams.lambdaL, LambdaL_T(fParams.t));
|
||||
if (fParams.specularIntegral > 0.0)
|
||||
fprintf(fp, "%% int_x=0^infty B(x) dx / Bext = %lf\n", fParams.specularIntegral);
|
||||
fprintf(fp, "%% xiP(0) = %lf, xiP(t) = %lf\n", fParams.xi0, XiP_T(fParams.t));
|
||||
fprintf(fp, "%% Mean Free Path = %lf\n", fParams.meanFreePath);
|
||||
fprintf(fp, "%% Film Thickness = %lf\n", fParams.filmThickness);
|
||||
if (fParams.specular)
|
||||
fprintf(fp, "%% Boundary Conditions: Specular\n");
|
||||
else
|
||||
fprintf(fp, "%% Boundary Conditions: Diffuse\n");
|
||||
if (fParams.rgeFileName.Length() > 0)
|
||||
fprintf(fp, "%% rge file name : %s\n", fParams.rgeFileName.Data());
|
||||
if (fParams.meanB != 0.0)
|
||||
fprintf(fp, "%% Mean Field/Bext = %lf\n", fParams.meanB);
|
||||
fprintf(fp, "%%\n");
|
||||
|
||||
// write data
|
||||
fprintf(fp, "%% Data --------------------------------------\n");
|
||||
fprintf(fp, "%% z (nm), B/B_0 \n");
|
||||
if (fSpecular) {
|
||||
if (fParams.specular) {
|
||||
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));
|
||||
fprintf(fp, "%lf, %lf\n", f_dz * XiP_T(fParams.t) * (Double_t)i, (*fFieldDiffuse)(i));
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -9,6 +9,7 @@
|
||||
#include <fftw3.h>
|
||||
|
||||
#include <Rtypes.h>
|
||||
#include <TString.h>
|
||||
|
||||
#include <Eigen/Core>
|
||||
#include <Eigen/LU>
|
||||
@ -17,32 +18,43 @@ using namespace Eigen;
|
||||
const Int_t PippardFourierPoints = 200000;
|
||||
const Int_t PippardDiffusePoints = 500;
|
||||
|
||||
typedef struct {
|
||||
Bool_t valid;
|
||||
Double_t t; // reduced temperature
|
||||
Double_t lambdaL; // lambda London
|
||||
Double_t xi0; // Pippard coherence length
|
||||
Double_t meanFreePath; // mean free path
|
||||
Double_t filmThickness; // film thickness
|
||||
Bool_t specular; // true -> specular reflection boundary conditions, false -> diffuse scattering boundary conditions
|
||||
Double_t specularIntegral; // int_x=0^infty B(x) dx / Bext
|
||||
TString rgeFileName;
|
||||
TString outputFileName;
|
||||
Double_t meanB; // int_x=0^infty B(x) n(x) dx / int_x=0^infty n(x) dx / Bext
|
||||
} PippardParams;
|
||||
|
||||
class PPippard
|
||||
{
|
||||
public:
|
||||
PPippard(Double_t temp, Double_t lambdaL, Double_t xi0, Double_t meanFreePath, Double_t filmThickness, Bool_t specular = true);
|
||||
PPippard(const PippardParams ¶ms);
|
||||
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 SetTemp(Double_t temp) { fParams.t = temp; }
|
||||
virtual void SetLambdaL(Double_t lambdaL) { fParams.lambdaL = lambdaL; }
|
||||
virtual void SetXi0(Double_t xi0) { fParams.xi0 = xi0; }
|
||||
virtual void SetMeanFreePath(Double_t meanFreePath) { fParams.meanFreePath = meanFreePath; }
|
||||
virtual void SetFilmThickness(Double_t thickness) { fParams.filmThickness = thickness; }
|
||||
virtual void SetSpecular(Bool_t specular) { fParams.specular = specular; }
|
||||
virtual void SetMeanB(Double_t meanB) { fParams.meanB = meanB; }
|
||||
|
||||
virtual void CalculateField();
|
||||
|
||||
virtual Double_t GetStepDistance() { return f_dz; }
|
||||
virtual Double_t GetMagneticField(const Double_t x) const;
|
||||
|
||||
virtual void SaveField(const char *fileName);
|
||||
virtual void SaveField();
|
||||
|
||||
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
|
||||
PippardParams fParams;
|
||||
|
||||
Double_t f_dx; // dx = xiPT dq
|
||||
Double_t f_dz; // spatial step size
|
||||
|
||||
@ -10,65 +10,351 @@ 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 << "usage: nonlocal <inputFile>";
|
||||
cout << endl << "The input file has the following structure:";
|
||||
cout << endl << "% input file for nonlocal";
|
||||
cout << endl << "% --------------------------------";
|
||||
cout << endl << "% input parameters:";
|
||||
cout << endl << "% reduced temperature";
|
||||
cout << endl << "% lambda London";
|
||||
cout << endl << "% xi0";
|
||||
cout << endl << "% mean free path";
|
||||
cout << endl << "% film thickness";
|
||||
cout << endl << "% boundary conditions";
|
||||
cout << endl << "% [rge input file name], i.e. this is optional";
|
||||
cout << endl << "% [output file name], i.e. this is optional";
|
||||
cout << endl;
|
||||
cout << endl << "reducedTemp = 0.8287";
|
||||
cout << endl << "lambdaL = 30.0";
|
||||
cout << endl << "xi0 = 380.0";
|
||||
cout << endl << "meanFreePath = 12000.0";
|
||||
cout << endl << "filmThickness = 5000.0";
|
||||
cout << endl << "specular = 1";
|
||||
cout << endl << "rgeFileName = /afs/psi.ch/project/nemu/analysis/2009/Nonlocal/trimsp/InSne060.rge";
|
||||
cout << endl << "outputFileName = In_37_T2P83.dat";
|
||||
cout << endl << endl;
|
||||
cout << endl << "lines starting with a '%' or a '#' are comment lines.";
|
||||
cout << endl << "empty lines are ignored.";
|
||||
cout << endl << "all lengths given in (nm)";
|
||||
cout << endl << " if <fileName> is given, the field profile will be saved";
|
||||
cout << endl << endl;
|
||||
}
|
||||
|
||||
//-----------------------------------------------------------------------------------------
|
||||
/**
|
||||
*
|
||||
*/
|
||||
int readInputFile(char *fln, PippardParams ¶m)
|
||||
{
|
||||
FILE *fp;
|
||||
|
||||
fp = fopen(fln, "r");
|
||||
if (fp == NULL) {
|
||||
cout << endl << "Couldn't read input file " << fln;
|
||||
cout << endl;
|
||||
return 0;
|
||||
}
|
||||
|
||||
char line[128], str[128];
|
||||
TString tstr;
|
||||
Double_t dval;
|
||||
Int_t result, ival;
|
||||
while (!feof(fp)) {
|
||||
fgets(line, sizeof(line), fp);
|
||||
tstr = line;
|
||||
if (tstr.BeginsWith("#") || tstr.BeginsWith("%") || tstr.IsWhitespace()) {
|
||||
continue;
|
||||
}
|
||||
|
||||
if (tstr.BeginsWith("reducedTemp =")) {
|
||||
result = sscanf(line, "reducedTemp = %lf", &dval);
|
||||
if (result != 1) {
|
||||
cout << endl << "**ERROR** while trying to extracted the reduced temperature";
|
||||
cout << endl << "line = " << line;
|
||||
cout << endl;
|
||||
param.valid = false;
|
||||
break;
|
||||
}
|
||||
param.t = dval;
|
||||
}
|
||||
|
||||
if (tstr.BeginsWith("lambdaL =")) {
|
||||
result = sscanf(line, "lambdaL = %lf", &dval);
|
||||
if (result != 1) {
|
||||
cout << endl << "**ERROR** while trying to extracted the London penetration length";
|
||||
cout << endl << "line = " << line;
|
||||
cout << endl;
|
||||
param.valid = false;
|
||||
break;
|
||||
}
|
||||
param.lambdaL = dval;
|
||||
}
|
||||
|
||||
if (tstr.BeginsWith("xi0 =")) {
|
||||
result = sscanf(line, "xi0 = %lf", &dval);
|
||||
if (result != 1) {
|
||||
cout << endl << "**ERROR** while trying to extracted the Pippard coherence length";
|
||||
cout << endl << "line = " << line;
|
||||
cout << endl;
|
||||
param.valid = false;
|
||||
break;
|
||||
}
|
||||
param.xi0 = dval;
|
||||
}
|
||||
|
||||
if (tstr.BeginsWith("meanFreePath =")) {
|
||||
result = sscanf(line, "meanFreePath = %lf", &dval);
|
||||
if (result != 1) {
|
||||
cout << endl << "**ERROR** while trying to extracted the mean free path";
|
||||
cout << endl << "line = " << line;
|
||||
cout << endl;
|
||||
param.valid = false;
|
||||
break;
|
||||
}
|
||||
param.meanFreePath = dval;
|
||||
}
|
||||
|
||||
if (tstr.BeginsWith("filmThickness =")) {
|
||||
result = sscanf(line, "filmThickness = %lf", &dval);
|
||||
if (result != 1) {
|
||||
cout << endl << "**ERROR** while trying to extracted the film thickness";
|
||||
cout << endl << "line = " << line;
|
||||
cout << endl;
|
||||
param.valid = false;
|
||||
break;
|
||||
}
|
||||
param.filmThickness = dval;
|
||||
}
|
||||
|
||||
if (tstr.BeginsWith("specular =")) {
|
||||
result = sscanf(line, "specular = %d", &ival);
|
||||
if (result != 1) {
|
||||
cout << endl << "**ERROR** while trying to extracted the boundary conditions";
|
||||
cout << endl << "line = " << line;
|
||||
cout << endl;
|
||||
param.valid = false;
|
||||
break;
|
||||
}
|
||||
if (ival == 1)
|
||||
param.specular = true;
|
||||
else
|
||||
param.specular = false;
|
||||
}
|
||||
|
||||
if (tstr.BeginsWith("rgeFileName =")) {
|
||||
result = sscanf(line, "rgeFileName = %s", str);
|
||||
if (result != 1) {
|
||||
cout << endl << "**ERROR** while trying to extracted the rge file name";
|
||||
cout << endl << "line = " << line;
|
||||
cout << endl;
|
||||
param.valid = false;
|
||||
break;
|
||||
}
|
||||
param.rgeFileName = TString(str);
|
||||
}
|
||||
|
||||
if (tstr.BeginsWith("outputFileName =")) {
|
||||
result = sscanf(line, "outputFileName = %s", str);
|
||||
if (result != 1) {
|
||||
cout << endl << "**ERROR** while trying to extracted the output file name";
|
||||
cout << endl << "line = " << line;
|
||||
cout << endl;
|
||||
param.valid = false;
|
||||
break;
|
||||
}
|
||||
param.outputFileName = TString(str);
|
||||
}
|
||||
}
|
||||
|
||||
fclose(fp);
|
||||
|
||||
return 1;
|
||||
}
|
||||
|
||||
//-----------------------------------------------------------------------------------------
|
||||
/**
|
||||
*
|
||||
*/
|
||||
int readRgeFile(const char *fln, vector<Double_t> &x, vector<Double_t> &n)
|
||||
{
|
||||
FILE *fp;
|
||||
|
||||
fp = fopen(fln, "r");
|
||||
if (fp == NULL) {
|
||||
cout << endl << "Couldn't read rge file " << fln;
|
||||
cout << endl;
|
||||
return 0;
|
||||
}
|
||||
|
||||
char line[128];
|
||||
int lineNo = 1;
|
||||
int result;
|
||||
double valX, valN;
|
||||
// read first line and ignore it since it is the header line
|
||||
fgets(line, sizeof(line), fp);
|
||||
while (!feof(fp)) {
|
||||
fgets(line, sizeof(line), fp);
|
||||
result = sscanf(line, "%lf %lf", &valX, &valN);
|
||||
if (result != 2) {
|
||||
cout << endl << "**ERROR** in line no " << lineNo;
|
||||
cout << endl << "Couldn't extract x and n(x), will quit.";
|
||||
cout << endl << endl;
|
||||
fclose(fp);
|
||||
return 0;
|
||||
}
|
||||
x.push_back(valX/10.0); // x in nm
|
||||
n.push_back(valN);
|
||||
lineNo++;
|
||||
}
|
||||
|
||||
// normalize n(z)
|
||||
Double_t norm = 0.0;
|
||||
for (unsigned int i=0; i<n.size(); i++) {
|
||||
norm += n[i];
|
||||
}
|
||||
for (unsigned int i=0; i<n.size(); i++) {
|
||||
n[i] /= norm * (x[1]-x[0]);
|
||||
}
|
||||
|
||||
/*
|
||||
for (unsigned int i=0; i<x.size(); i++) {
|
||||
cout << endl << i << ": " << x[i] << ", " << n[i];
|
||||
}
|
||||
cout << endl;
|
||||
*/
|
||||
|
||||
fclose(fp);
|
||||
|
||||
return 1;
|
||||
}
|
||||
|
||||
//-----------------------------------------------------------------------------------------
|
||||
/**
|
||||
*
|
||||
*/
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
if ((argc < 7) || (argc>8)){
|
||||
if (argc != 2) {
|
||||
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);
|
||||
PippardParams params;
|
||||
|
||||
char fln[1024];
|
||||
if (argc == 8) {
|
||||
strncpy(fln, argv[7], sizeof(fln));
|
||||
} else {
|
||||
strncpy(fln, "", sizeof(fln));
|
||||
}
|
||||
// init params
|
||||
params.valid = true;
|
||||
params.t = -1.0;
|
||||
params.lambdaL = -1.0;
|
||||
params.xi0 = -1.0;
|
||||
params.meanFreePath = -1.0;
|
||||
params.filmThickness = -1.0;
|
||||
params.specular = true;
|
||||
params.specularIntegral = -1.0;
|
||||
params.rgeFileName = "";
|
||||
params.outputFileName = "";
|
||||
params.meanB = 0.0;
|
||||
|
||||
|
||||
Bool_t specular;
|
||||
if (!strcmp(argv[6], "1"))
|
||||
specular = true;
|
||||
else if (!strcmp(argv[6], "0"))
|
||||
specular = false;
|
||||
else {
|
||||
syntax();
|
||||
if (!readInputFile(argv[1], params)) {
|
||||
cout << endl << "**ERROR** Couldn't open the input file " << argv[1] << ". Will quit.";
|
||||
cout << endl << endl;
|
||||
return 0;
|
||||
}
|
||||
|
||||
cout << endl << ">> temp = " << temp;
|
||||
cout << endl << ">> lambdaL = " << lambdaL;
|
||||
cout << endl << ">> xi0 = " << xi0;
|
||||
cout << endl << ">> meanFreePath = " << meanFreePath;
|
||||
cout << endl << ">> filmThickness = " << filmThickness;
|
||||
if (specular)
|
||||
// check params
|
||||
if (!params.valid) {
|
||||
cout << endl << "**ERROR** parameters are not valid will quit.";
|
||||
cout << endl << endl;
|
||||
return 0;
|
||||
}
|
||||
if (params.t < 0.0) {
|
||||
cout << endl << "**ERROR** reduced temperature not defined, will quit.";
|
||||
cout << endl << endl;
|
||||
return 0;
|
||||
}
|
||||
if (params.lambdaL < 0.0) {
|
||||
cout << endl << "**ERROR** London pentration depth not defined, will quit.";
|
||||
cout << endl << endl;
|
||||
return 0;
|
||||
}
|
||||
if (params.xi0 < 0.0) {
|
||||
cout << endl << "**ERROR** Pippard coherence length not defined, will quit.";
|
||||
cout << endl << endl;
|
||||
return 0;
|
||||
}
|
||||
if (params.meanFreePath < 0.0) {
|
||||
cout << endl << "**ERROR** mean free path not defined, will set it to 1.0e6 nm.";
|
||||
cout << endl << endl;
|
||||
params.meanFreePath = 1.0e6;
|
||||
}
|
||||
if (params.filmThickness < 0.0) {
|
||||
cout << endl << "**ERROR** film thickness not defined, will set it to 1.0e6 nm.";
|
||||
cout << endl << endl;
|
||||
params.filmThickness = 1.0e6;
|
||||
}
|
||||
|
||||
int result;
|
||||
vector<Double_t> x, n;
|
||||
if (params.rgeFileName.Length() > 0) {
|
||||
result = readRgeFile(params.rgeFileName.Data(), x, n);
|
||||
if (result == 0)
|
||||
params.rgeFileName = "";
|
||||
}
|
||||
|
||||
cout << endl << ">> temp = " << params.t;
|
||||
cout << endl << ">> lambdaL = " << params.lambdaL;
|
||||
cout << endl << ">> xi0 = " << params.xi0;
|
||||
cout << endl << ">> meanFreePath = " << params.meanFreePath;
|
||||
cout << endl << ">> filmThickness = " << params.filmThickness;
|
||||
if (params.specular)
|
||||
cout << endl << ">> specular = true";
|
||||
else
|
||||
cout << endl << ">> specular = false";
|
||||
cout << endl;
|
||||
|
||||
PPippard *pippard = new PPippard(temp, lambdaL, xi0, meanFreePath, filmThickness, specular);
|
||||
PPippard *pippard = new PPippard(params);
|
||||
|
||||
pippard->CalculateField();
|
||||
|
||||
if (strlen(fln) > 0)
|
||||
pippard->SaveField(fln);
|
||||
// check if it is necessary to calculate the <B(x)>
|
||||
if (x.size() > 0) {
|
||||
Double_t dz = pippard->GetStepDistance();
|
||||
Int_t max = (Int_t)((Double_t)(x[x.size()-1])/dz);
|
||||
//cout << endl << "max = " << max;
|
||||
Double_t mean = 0.0;
|
||||
Double_t nn;
|
||||
for (Int_t i=0; i<max; i++) {
|
||||
// find propper z value index of n(z)
|
||||
Int_t idx = -1;
|
||||
for (unsigned int j=0; j<x.size(); j++) {
|
||||
if (x[j] >= i*dz) {
|
||||
idx = j;
|
||||
break;
|
||||
}
|
||||
}
|
||||
// get n(z) at the proper index
|
||||
if (idx == -1) {
|
||||
nn = 0.0;
|
||||
} else if (idx == 0) {
|
||||
nn = n[0];
|
||||
} else {
|
||||
nn = n[idx-1] + (n[idx]-n[idx-1]) * (i*dz-x[idx-1])/(x[idx]-x[idx-1]);
|
||||
}
|
||||
//cout << endl << i << ", idx = " << idx << ", i*dz = " << i*dz << ", x[idx] = " << x[idx] << ", x[idx-1] = " << x[idx-1] << ", nn = " << nn;
|
||||
mean += pippard->GetMagneticField(i*dz) * nn;
|
||||
}
|
||||
mean *= dz/(x[1]-x[0]);
|
||||
|
||||
cout << endl << ">> mean field = " << mean;
|
||||
pippard->SetMeanB(mean);
|
||||
}
|
||||
|
||||
if (params.outputFileName.Length() > 0)
|
||||
pippard->SaveField();
|
||||
|
||||
cout << endl << ">> magnetic field = " << pippard->GetMagneticField(0.0);
|
||||
cout << endl;
|
||||
|
||||
Reference in New Issue
Block a user