From 6e596915c82c4e7889e9defac9f3c35557422f42 Mon Sep 17 00:00:00 2001 From: nemu Date: Mon, 6 Jul 2009 12:29:36 +0000 Subject: [PATCH] changed to file input handling --- src/tests/nonlocal/In_37_T2P83_E14P1.inp | 21 ++ src/tests/nonlocal/Makefile | 6 - src/tests/nonlocal/PPippard.cpp | 63 ++-- src/tests/nonlocal/PPippard.h | 40 ++- src/tests/nonlocal/nonlocal.cpp | 356 ++++++++++++++++++++--- 5 files changed, 403 insertions(+), 83 deletions(-) create mode 100644 src/tests/nonlocal/In_37_T2P83_E14P1.inp diff --git a/src/tests/nonlocal/In_37_T2P83_E14P1.inp b/src/tests/nonlocal/In_37_T2P83_E14P1.inp new file mode 100644 index 00000000..7a4afe6d --- /dev/null +++ b/src/tests/nonlocal/In_37_T2P83_E14P1.inp @@ -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 + diff --git a/src/tests/nonlocal/Makefile b/src/tests/nonlocal/Makefile index 661bbb60..bae852ff 100644 --- a/src/tests/nonlocal/Makefile +++ b/src/tests/nonlocal/Makefile @@ -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 diff --git a/src/tests/nonlocal/PPippard.cpp b/src/tests/nonlocal/PPippard.cpp index 5390a755..83357782 100644 --- a/src/tests/nonlocal/PPippard.cpp +++ b/src/tests/nonlocal/PPippard.cpp @@ -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> 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 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 #include +#include #include #include @@ -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 diff --git a/src/tests/nonlocal/nonlocal.cpp b/src/tests/nonlocal/nonlocal.cpp index 1a62ce4f..af56fd02 100644 --- a/src/tests/nonlocal/nonlocal.cpp +++ b/src/tests/nonlocal/nonlocal.cpp @@ -10,65 +10,351 @@ using namespace std; #include "PPippard.h" - +//----------------------------------------------------------------------------------------- +/** + * + */ void syntax() { - cout << endl << "usage: nonlocal temp lambdaL xi0 meanFreePath filmThickness specular []"; - 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 is given, the field profile will be saved"; + cout << endl << "usage: nonlocal "; + 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 << 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 &x, vector &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; i8)){ + 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 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 + 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= 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;