some improvement callculating <B>

src/tests/nonlocal/In_37_T2P83_E14P1.inp

@@ -7,15 +7,19 @@
 %   mean free path
 %   film thickness
 %   boundary conditions
+%   [Bext in (G)]
+%   [deadLayer in (nm)]
 %   [rge input file name]
 %   [output file name]
 
-reducedTemp = 0.8287
+reducedTemp = 0.8492
-lambdaL = 25.0
+lambdaL = 33.6
 xi0 = 380.0
 meanFreePath = 12000.0
 filmThickness = 500000.0
 specular = 1
+Bext = 41.17
+deadLayer = 0.42
 rgeFileName = /afs/psi.ch/project/nemu/analysis/2009/Nonlocal/trimsp/InSne141.rge
 outputFileName =  In_37_T2P83_E14P1.dat
 

src/tests/nonlocal/PPippard.cpp

@@ -392,8 +392,10 @@ void PPippard::SaveField()
     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)
+  if (fParams.meanB != 0.0) {
+    fprintf(fp, "%%   Mean Distance       = %lf\n", fParams.meanX);
     fprintf(fp, "%%   Mean Field/Bext     = %lf\n", fParams.meanB);
+  }
   fprintf(fp, "%%\n");
 
   // write data

src/tests/nonlocal/PPippard.h

@@ -27,9 +27,12 @@ typedef struct {
   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
+  Double_t b_ext;         // external field in (G)
+  Double_t deadLayer;     // dead layer in (nm)
   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 meanX;         // int_x=0^infty x n(x) dx / int_x=0^infty n(x) dx
 } PippardParams;
 
 class PPippard
@@ -44,6 +47,7 @@ class PPippard
     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 CalculateField();

src/tests/nonlocal/nonlocal.cpp

@@ -27,6 +27,8 @@ void syntax()
   cout << endl << "%   mean free path";
   cout << endl << "%   film thickness";
   cout << endl << "%   boundary conditions";
+  cout << endl << "%   [Bext in (G)], i.e. this is optional";
+  cout << endl << "%   [dead layer in (nm)], i.e. this is optional";
   cout << endl << "%   [rge input file name], i.e. this is optional";
   cout << endl << "%   [output file name], i.e. this is optional";
   cout << endl;
@@ -36,6 +38,8 @@ void syntax()
   cout << endl << "meanFreePath = 12000.0";
   cout << endl << "filmThickness = 5000.0";
   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 << "outputFileName = In_37_T2P83.dat";
   cout << endl << endl;
@@ -146,6 +150,30 @@ int readInputFile(char *fln, PippardParams &param)
         param.specular = false;
     }
 
+    if (tstr.BeginsWith("Bext = ")) {
+      result = sscanf(line, "Bext = %lf", &dval);
+      if (result != 1) {
+        cout << endl << "**ERROR** while trying to extracted the external field value.";
+        cout << endl << "line = " << line;
+        cout << endl;
+        param.valid = false;
+        break;
+      }
+      param.b_ext = dval;
+    }
+
+    if (tstr.BeginsWith("deadLayer = ")) {
+      result = sscanf(line, "deadLayer = %lf", &dval);
+      if (result != 1) {
+        cout << endl << "**ERROR** while trying to extracted the dead layer value.";
+        cout << endl << "line = " << line;
+        cout << endl;
+        param.valid = false;
+        break;
+      }
+      param.deadLayer = dval;
+    }
+
     if (tstr.BeginsWith("rgeFileName =")) {
       result = sscanf(line, "rgeFileName = %s", str);
       if (result != 1) {
@@ -255,8 +283,11 @@ int main(int argc, char *argv[])
   params.filmThickness = -1.0;
   params.specular = true;
   params.specularIntegral = -1.0;
+  params.b_ext = -1.0;
+  params.deadLayer = -1.0;
   params.rgeFileName = "";
   params.outputFileName = "";
+  params.meanX = 0.0;
   params.meanB = 0.0;
 
   if (!readInputFile(argv[1], params)) {
@@ -322,35 +353,29 @@ int main(int argc, char *argv[])
 
   // check if it is necessary to calculate the <B(x)>
   if (x.size() > 0) {
-    Double_t dz = pippard->GetStepDistance();
+    if ((params.b_ext == -1.0) || (params.deadLayer == -1.0)) {
-    Int_t max = (Int_t)((Double_t)(x[x.size()-1])/dz);
+      cout << endl << "**ERROR** Bext or deadLayer missing :-(" << endl;
-//cout << endl << "max = " << max;
+      return 0;
-    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;
+    Double_t meanX = 0.0;
-    pippard->SetMeanB(mean);
+    for (unsigned int i=0; i<x.size(); i++) {
+      meanX += x[i]*n[i];
+    }
+    meanX *= (x[1]-x[0]);
+
+    Double_t meanB = 0.0;
+    for (unsigned int i=0; i<x.size()-1; i++) {
+      if (x[i] <= params.deadLayer)
+        meanB += 1.0 * n[i];
+      else
+        meanB += pippard->GetMagneticField(x[i]-params.deadLayer) * n[i];
+    }
+    meanB *= (x[1]-x[0]);
+
+    cout << endl << ">> mean x = " << meanX << ", mean field = " << params.b_ext * meanB;
+    pippard->SetMeanX(meanX);
+    pippard->SetMeanB(meanB);
   }
 
   if (params.outputFileName.Length() > 0)