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

adopted to newest Eigen version and at the same time added higher functionality

Browse Source
This commit is contained in:
suter_a
2011-10-18 08:34:46 +00:00
parent d9a83c6c19
commit 1f0ee0662d
4 changed files with 256 additions and 131 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
src/tests/nonlocal/Makefile
View File

@ -87,6 +87,9 @@ nonlocal$(SUFFIX): $(OBJS)
#
clean:; @rm -f $(OBJS)
@echo "---> removing $(OBJS)"
install:
cp nonlocal $(HOME)/bin
#
$(OBJS): %.o: %.cpp

185
src/tests/nonlocal/PPippard.cpp
View File

@ -3,6 +3,7 @@
// $Id$
// -----------------------------------------------------------------------
#include <cassert>
#include <cstdio>
#include <cmath>
@ -25,9 +26,6 @@ PPippard::PPippard(const PippardParams &params) : fParams(params)
fFieldq = 0;
fFieldB = 0;
fSecondDerivativeMatrix = 0;
fKernelMatrix = 0;
fBoundaryCondition = 0;
fFieldDiffuse = 0;
f_dx = 0.02;
@ -48,29 +46,48 @@ PPippard::~PPippard()
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;
}
}
//-----------------------------------------------------------------------------------------------------------
/**
*
*/
void PPippard::DumpParams()
{
cout << endl << ">> Parameters: ";
cout << endl << ">> reduced temperature: " << fParams.t;
cout << endl << ">> lambdaL : " << fParams.lambdaL << " (nm)";
cout << endl << ">> xi0 : " << fParams.xi0 << " (nm)";
cout << endl << ">> mean free path : " << fParams.meanFreePath << " (nm)";
cout << endl << ">> film thickness : " << fParams.filmThickness << " (nm)";
if (fParams.specular)
cout << endl << ">> specular scattering boundary conditions";
else
cout << endl << ">> diffuse scattering boundary conditions";
cout << endl << ">> Bext : " << fParams.b_ext << " (G)";
cout << endl << ">> dead layer : " << fParams.deadLayer << " (nm)";
for (UInt_t i=0; i<fParams.rgeFileName.size(); i++)
cout << endl << ">> rge-file-name : " << fParams.rgeFileName[i].Data();
if (fParams.outputFileName.Length() > 0)
cout << endl << ">> output file name : " << fParams.outputFileName.Data();
else
cout << endl << ">> output file name : n/a";
if (fParams.outputFileNameBmean.Length() > 0)
cout << endl << ">> output file name Bmean : " << fParams.outputFileNameBmean.Data();
else
cout << endl << ">> output file name Bmean : n/a";
cout << endl << ">> input file name : " << fParams.inputFileName.Data();
cout << endl << endl;
}
//-----------------------------------------------------------------------------------------------------------
/**
*
@ -248,8 +265,6 @@ void PPippard::CalculateFieldSpecular()
}
}
cout << endl << "debug> fShift = " << fShift;
for (Int_t i=0; i<PippardFourierPoints; i++) {
fFieldB[i][1] /= norm;
}
@ -282,69 +297,64 @@ cout << endl << "debug> fShift = " << fShift;
*/
void PPippard::CalculateFieldDiffuse()
{
MatrixXd fSecondDerivativeMatrix(PippardDiffusePoints+1, PippardDiffusePoints+1); // 2nd derivative matrix
MatrixXd fKernelMatrix(PippardDiffusePoints+1, PippardDiffusePoints+1); // kernel matrix
VectorXd fBoundaryCondition(PippardDiffusePoints+1); // boundary condition vector
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(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(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;
cout << endl << ">> f_dz = " << f_dz << ", invL = " << invL;
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;
// 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;
}
fSecondDerivativeMatrix.setZero();
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);
}
fKernelMatrix.setZero();
// 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;
}
fBoundaryCondition.setZero();
fBoundaryCondition(0) = 1.0;
//cout << endl << "fBoundaryCondition = " << *fBoundaryCondition << endl;
if (fFieldDiffuse == 0) {
fFieldDiffuse = new VectorXd(PippardDiffusePoints+1);
fFieldDiffuse->setZero(PippardDiffusePoints+1);
if (fFieldDiffuse != 0) {
delete fFieldDiffuse;
fFieldDiffuse = 0;
}
fFieldDiffuse = new VectorXd(PippardDiffusePoints+1);
fFieldDiffuse->setZero();
// solve equation
*fSecondDerivativeMatrix = (*fSecondDerivativeMatrix)-(*fKernelMatrix);
fSecondDerivativeMatrix->lu().solve(*fBoundaryCondition, fFieldDiffuse);
fSecondDerivativeMatrix = fSecondDerivativeMatrix-fKernelMatrix;
*fFieldDiffuse = fSecondDerivativeMatrix.colPivHouseholderQr().solve(fBoundaryCondition);
// normalize field
Double_t norm = 0.0;
@ -352,6 +362,8 @@ cout << endl << ">> ampl = " << ampl << endl;
if (norm < (*fFieldDiffuse)(i))
norm = (*fFieldDiffuse)(i);
assert(norm != 0);
for (Int_t i=0; i<PippardDiffusePoints+1; i++)
(*fFieldDiffuse)(i) /= norm;
@ -392,12 +404,15 @@ void PPippard::SaveField()
fprintf(fp, "%% Boundary Conditions: Diffuse\n");
fprintf(fp, "%% Bext = %lf (G)\n", fParams.b_ext);
fprintf(fp, "%% deadLayer = %lf (nm)\n", fParams.deadLayer);
if (fParams.rgeFileName.Length() > 0)
fprintf(fp, "%% rge file name : %s\n", fParams.rgeFileName.Data());
if (fParams.meanB != 0.0) {
fprintf(fp, "%% Mean Distance = %lf (nm)\n", fParams.meanX);
fprintf(fp, "%% Mean Field/Bext = %lf, Mean Field = %lf (G)\n", fParams.meanB, fParams.meanB * fParams.b_ext);
fprintf(fp, "%% Var Field/Bext^2 = %lf, Var Field = %lf (G^2)\n", fParams.varB, fParams.varB * fParams.b_ext * fParams.b_ext);
for (UInt_t i=0; i<fParams.rgeFileName.size(); i++) {
fprintf(fp, "%% rge file name : %s\n", fParams.rgeFileName[i].Data());
}
for (UInt_t i=0; i<fParams.energy.size(); i++) {
fprintf(fp, "%%------------------------------------------\n");
fprintf(fp, "%% Implantation Energy = %lf (eV)\n", fParams.energy[i]);
fprintf(fp, "%% Mean Distance = %lf (nm)\n", fParams.meanX[i]);
fprintf(fp, "%% Mean Field/Bext = %lf, Mean Field = %lf (G)\n", fParams.meanB[i], fParams.meanB[i] * fParams.b_ext);
fprintf(fp, "%% Var Field/Bext^2 = %lf, Var Field = %lf (G^2)\n", fParams.varB[i], fParams.varB[i] * fParams.b_ext * fParams.b_ext);
}
fprintf(fp, "%%\n");
@ -417,6 +432,38 @@ void PPippard::SaveField()
fclose(fp);
}
//-----------------------------------------------------------------------------------------------------------
/**
*
*/
void PPippard::SaveBmean()
{
FILE *fp;
fp = fopen(fParams.outputFileNameBmean.Data(), "w");
if (fp == NULL) {
cout << endl << "Coudln't open " << fParams.outputFileNameBmean.Data() << " for writting, sorry ...";
cout << endl << endl;
return;
}
// write header
fprintf(fp, "%% Header -----------------------\n");
fprintf(fp, "%% generated from input file: %s\n", fParams.inputFileName.Data());
fprintf(fp, "%%\n");
fprintf(fp, "%% Data -------------------------\n");
fprintf(fp, "%% energy (eV), meanZ (nm), meanB/Bext, meanB (G), varB/Bext^2, varB (G^2)\n");
// write data
for (UInt_t i=0; i<fParams.energy.size(); i++) {
fprintf(fp, "%lf, %lf, %lf, %lf, %lf, %lf\n", fParams.energy[i], fParams.meanX[i],
fParams.meanB[i], fParams.b_ext*fParams.meanB[i],
fParams.varB[i], (fParams.b_ext*fParams.b_ext)*fParams.varB[i]);
}
fclose(fp);
}
//-----------------------------------------------------------------------------------------------------------
/**
*

30
src/tests/nonlocal/PPippard.h
View File

@ -6,13 +6,15 @@
#ifndef _PPIPPARD_H_
#define _PPIPPARD_H_
#include <vector>
using namespace std;
#include <fftw3.h>
#include <Rtypes.h>
#include <TString.h>
#include <Eigen/Core>
#include <Eigen/LU>
#include <Eigen/Dense>
using namespace Eigen;
const Int_t PippardFourierPoints = 200000;
@ -29,11 +31,14 @@ typedef struct {
Double_t specularIntegral; // int_x=0^infty B(x) dx / Bext
Double_t b_ext; // external field in (G)
Double_t deadLayer; // dead layer in (nm)
TString rgeFileName;
TString inputFileName;
vector<TString> rgeFileName;
TString outputFileName;
Double_t meanB; // int_x=0^infty B(x) n(x) dx / int_x=0^infty n(x) dx / Bext
Double_t varB; // int_x=0^infty (B(x)-<B>)^2 n(x) dx / int_x=0^infty n(x) dx / Bext^2
Double_t meanX; // int_x=0^infty x n(x) dx / int_x=0^infty n(x) dx
TString outputFileNameBmean;
vector<Double_t> energy;
vector<Double_t> meanB; // int_x=0^infty B(x) n(x) dx / int_x=0^infty n(x) dx / Bext
vector<Double_t> varB; // int_x=0^infty (B(x)-<B>)^2 n(x) dx / int_x=0^infty n(x) dx / Bext^2
vector<Double_t> meanX; // int_x=0^infty x n(x) dx / int_x=0^infty n(x) dx
} PippardParams;
class PPippard
@ -42,15 +47,18 @@ class PPippard
PPippard(const PippardParams &params);
virtual ~PPippard();
virtual void DumpParams();
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 SetMeanX(Double_t meanX) { fParams.meanX = meanX; }
virtual void SetMeanB(Double_t meanB) { fParams.meanB = meanB; }
virtual void SetVarB(Double_t BB) { fParams.varB = BB; }
virtual void SetEnergy(Double_t energy) { fParams.energy.push_back(energy); }
virtual void SetMeanX(Double_t meanX) { fParams.meanX.push_back(meanX); }
virtual void SetMeanB(Double_t meanB) { fParams.meanB.push_back(meanB); }
virtual void SetVarB(Double_t BB) { fParams.varB.push_back(BB); }
virtual void CalculateField();
@ -58,6 +66,7 @@ class PPippard
virtual Double_t GetMagneticField(const Double_t x) const;
virtual void SaveField();
virtual void SaveBmean();
private:
PippardParams fParams;
@ -72,9 +81,6 @@ class PPippard
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;

169
src/tests/nonlocal/nonlocal.cpp
View File

@ -10,6 +10,12 @@ using namespace std;
#include "PPippard.h"
typedef struct {
Double_t energy;
vector<Double_t> x;
vector<Double_t> n;
} Pnofx;
//-----------------------------------------------------------------------------------------
/**
*
@ -29,8 +35,9 @@ void syntax()
cout << endl << "% specular : 1=specular, 0=diffuse";
cout << endl << "% [Bext] : external magnetic field (G), i.e. this is optional";
cout << endl << "% [deadLayer] : dead layer (nm), i.e. this is optional";
cout << endl << "% [rgeFileName] : rge input file name, i.e. this is optional";
cout << endl << "% [rgeFileName] : rge input file name, i.e. this is optional. Multiple rge-files possible";
cout << endl << "% [outputFileName] : output file name, i.e. this is optional";
cout << endl << "% [outputFileNameBmean] : output file name for <B> data";
cout << endl;
cout << endl << "reducedTemp = 0.8287";
cout << endl << "lambdaL = 30.0";
@ -40,8 +47,10 @@ void syntax()
cout << endl << "specular = 1";
cout << endl << "Bext = 47.14";
cout << endl << "deadLayer = 0.4";
cout << endl << "rgeFileName = /afs/psi.ch/project/nemu/analysis/2009/Nonlocal/trimsp/InSne060.rge";
cout << endl << "rgeFileName = /afs/psi.ch/project/nemu/analysis/2009/Nonlocal/trimsp/InSn_E6000.rge";
cout << endl << "rgeFileName = /afs/psi.ch/project/nemu/analysis/2009/Nonlocal/trimsp/InSn_E14100.rge";
cout << endl << "outputFileName = In_37_T2P83.dat";
cout << endl << "outputFileNameBmean = In_37_T2P83_Bmean.dat";
cout << endl << endl;
cout << endl << "lines starting with a '%' or a '#' are comment lines.";
cout << endl << "empty lines are ignored.";
@ -183,7 +192,7 @@ int readInputFile(char *fln, PippardParams &param)
param.valid = false;
break;
}
param.rgeFileName = TString(str);
param.rgeFileName.push_back(TString(str));
}
if (tstr.BeginsWith("outputFileName =")) {
@ -197,6 +206,18 @@ int readInputFile(char *fln, PippardParams &param)
}
param.outputFileName = TString(str);
}
if (tstr.BeginsWith("outputFileNameBmean =")) {
result = sscanf(line, "outputFileNameBmean = %s", str);
if (result != 1) {
cout << endl << "**ERROR** while trying to extracted the Bmean output file name";
cout << endl << "line = " << line;
cout << endl;
param.valid = false;
break;
}
param.outputFileNameBmean = TString(str);
}
}
fclose(fp);
@ -208,10 +229,28 @@ int readInputFile(char *fln, PippardParams &param)
/**
*
*/
int readRgeFile(const char *fln, vector<Double_t> &x, vector<Double_t> &n)
int readRgeFile(const char *fln, vector<Double_t> &x, vector<Double_t> &n, Double_t &energy)
{
FILE *fp;
x.clear();
n.clear();
TString str = TString(fln);
str.Remove(0, str.Index("_E")+2);
str.Remove(str.Index(".rge"));
if (str.Length() > 0) {
if (str.IsFloat()) {
energy = str.Atof();
} else {
cerr << endl << "**ERROR** couldn't extract energy from the file name '" << fln << "', will quit." << endl << endl;
return 0;
}
} else {
cerr << endl << "**ERROR** couldn't extract energy from the file name '" << fln << "', will quit." << endl << endl;
return 0;
}
fp = fopen(fln, "r");
if (fp == NULL) {
cout << endl << "Couldn't read rge file " << fln;
@ -253,15 +292,6 @@ int readRgeFile(const char *fln, vector<Double_t> &x, vector<Double_t> &n)
n[i] /= norm * (x[1]-x[0]);
}
/*
Double_t xx=0.0;
for (unsigned int i=0; i<x.size(); i++) {
xx += n[i];
}
xx *= (x[1]-x[0]);
cout << "debug> xx=" << xx << endl;
*/
fclose(fp);
return 1;
@ -291,10 +321,12 @@ int main(int argc, char *argv[])
params.specularIntegral = -1.0;
params.b_ext = -1.0;
params.deadLayer = -1.0;
params.rgeFileName = "";
params.rgeFileName.clear();
params.inputFileName = TString(argv[1]);
params.outputFileName = "";
params.meanX = 0.0;
params.meanB = 0.0;
params.outputFileNameBmean = "";
params.meanX.clear();
params.meanB.clear();
if (!readInputFile(argv[1], params)) {
cout << endl << "**ERROR** Couldn't open the input file " << argv[1] << ". Will quit.";
@ -335,12 +367,34 @@ int main(int argc, char *argv[])
}
int result;
Pnofx rgeSet;
vector<Pnofx> nOfx;
double energy;
vector<Double_t> x, n;
if (params.rgeFileName.Length() > 0) {
result = readRgeFile(params.rgeFileName.Data(), x, n);
if (result == 0)
params.rgeFileName = "";
for (UInt_t i=0; i<params.rgeFileName.size(); i++) {
x.clear();
n.clear();
energy = 0.0;
rgeSet.energy = 0.0;
rgeSet.x.clear();
rgeSet.n.clear();
result = readRgeFile(params.rgeFileName[i].Data(), x, n, energy);
if (result == 0) {
cerr << endl << "**ERROR** Couldn't read rge-file " << params.rgeFileName[i].Data() << ", will quit." << endl << endl;
return 0;
} else {
cout << endl << ">> read rge-file '" << params.rgeFileName[i].Data() << "'";
}
cout << endl << ">> energy = " << energy << " eV";
rgeSet.energy = energy;
rgeSet.x = x;
rgeSet.n = n;
nOfx.push_back(rgeSet);
}
cout << endl << ">> found " << nOfx.size() << " rge-files.";
cout << endl << ">> temp = " << params.t;
cout << endl << ">> lambdaL = " << params.lambdaL;
@ -348,60 +402,75 @@ int main(int argc, char *argv[])
cout << endl << ">> meanFreePath = " << params.meanFreePath;
cout << endl << ">> filmThickness = " << params.filmThickness;
if (params.specular)
cout << endl << ">> specular = true";
cout << endl << ">> specular = true";
else
cout << endl << ">> specular = false";
cout << endl << ">> specular = false";
cout << endl;
PPippard *pippard = new PPippard(params);
// pippard->DumpParams();
pippard->CalculateField();
// check if it is necessary to calculate the <B(x)> and <[B(x)-<B(x)>]^2>
if (x.size() > 0) {
if (nOfx.size() > 0) {
if ((params.b_ext == -1.0) || (params.deadLayer == -1.0)) {
cout << endl << "**ERROR** Bext or deadLayer missing :-(" << endl;
return 0;
}
Double_t meanX = 0.0;
for (unsigned int i=0; i<x.size(); i++) {
meanX += x[i]*n[i];
}
meanX *= (x[1]-x[0]);
Double_t meanB = 0.0;
Double_t varB = 0.0;
Double_t BB = 0.0;
for (unsigned int i=0; i<x.size()-1; i++) {
if (x[i] <= params.deadLayer) {
meanB += 1.0 * n[i];
} else if (x[i] > params.deadLayer+params.filmThickness) {
meanB += 1.0 * n[i];
} else {
BB = pippard->GetMagneticField(x[i]-params.deadLayer);
meanB += BB * n[i];
for (UInt_t i=0; i<nOfx.size(); i++) {
meanX = 0.0;
for (unsigned int j=0; j<nOfx[i].x.size(); j++) {
meanX += nOfx[i].x[j]*nOfx[i].n[j];
}
}
meanB *= (x[1]-x[0]);
for (unsigned int i=0; i<x.size()-1; i++) {
if (x[i] > params.deadLayer) {
BB = pippard->GetMagneticField(x[i]-params.deadLayer);
varB += (BB-meanB) * (BB-meanB) * n[i];
}
}
varB *= (x[1]-x[0]);
meanX *= (nOfx[i].x[1]-nOfx[i].x[0]);
cout << endl << ">> mean x = " << meanX << " (nm), mean field = " << params.b_ext * meanB << " (G)";
cout << " <(B-<B>)^2> = " << varB * params.b_ext * params.b_ext << " (G^2)";
pippard->SetMeanX(meanX);
pippard->SetVarB(varB);
pippard->SetMeanB(meanB);
meanB = 0.0;
varB = 0.0;
BB = 0.0;
for (unsigned int j=0; j<nOfx[i].x.size()-1; j++) {
if (nOfx[i].x[j] <= params.deadLayer) {
meanB += 1.0 * nOfx[i].n[j];
} else if (nOfx[i].x[j] > params.deadLayer+params.filmThickness) {
meanB += 1.0 * nOfx[i].n[j];
} else {
BB = pippard->GetMagneticField(nOfx[i].x[j]-params.deadLayer);
meanB += BB * nOfx[i].n[j];
}
}
meanB *= (nOfx[i].x[1]-nOfx[i].x[0]);
for (unsigned int j=0; j<nOfx[i].x.size()-1; j++) {
if (nOfx[i].x[j] > params.deadLayer) {
BB = pippard->GetMagneticField(nOfx[i].x[j]-params.deadLayer);
varB += (BB-meanB) * (BB-meanB) * nOfx[i].n[j];
}
}
varB *= (nOfx[i].x[1]-nOfx[i].x[0]);
cout << endl << "----------------------------------------";
cout << endl << ">> energy = " << nOfx[i].energy << "eV";
cout << endl << ">> mean x = " << meanX << " (nm), mean field = " << params.b_ext * meanB << " (G),";
cout << " <(B-<B>)^2> = " << varB * params.b_ext * params.b_ext << " (G^2)";
pippard->SetEnergy(nOfx[i].energy);
pippard->SetMeanX(meanX);
pippard->SetVarB(varB);
pippard->SetMeanB(meanB);
}
}
if (params.outputFileName.Length() > 0)
pippard->SaveField();
if (params.outputFileNameBmean.Length() > 0)
pippard->SaveBmean();
cout << endl << endl;
if (pippard) {