From 76b2364169ca72ace57ccf231c74cbefbcb84654 Mon Sep 17 00:00:00 2001
From: nemu <nemu>
Date: Tue, 17 Aug 2010 11:54:59 +0000
Subject: [PATCH] proper implementation of the variance B

---
 src/tests/nonlocal/PPippard.cpp |  6 +++---
 src/tests/nonlocal/PPippard.h   |  4 ++--
 src/tests/nonlocal/nonlocal.cpp | 18 +++++++++++-------
 3 files changed, 16 insertions(+), 12 deletions(-)

diff --git a/src/tests/nonlocal/PPippard.cpp b/src/tests/nonlocal/PPippard.cpp
index ed639435..79edc591 100644
--- a/src/tests/nonlocal/PPippard.cpp
+++ b/src/tests/nonlocal/PPippard.cpp
@@ -395,9 +395,9 @@ void PPippard::SaveField()
   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\n", fParams.meanX);
-    fprintf(fp, "%%   Mean Field/Bext         = %lf\n", fParams.meanB);
-    fprintf(fp, "%%   2nd Moment Field/Bext^2 = %lf\n", fParams.secondMomentB);
+    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);
   }
   fprintf(fp, "%%\n");
 
diff --git a/src/tests/nonlocal/PPippard.h b/src/tests/nonlocal/PPippard.h
index 386a16c1..6903acec 100644
--- a/src/tests/nonlocal/PPippard.h
+++ b/src/tests/nonlocal/PPippard.h
@@ -32,7 +32,7 @@ typedef struct {
   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 secondMomentB; // int_x=0^infty B^2(x) n(x) dx / int_x=0^infty n(x) dx / Bext^2
+  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
 } PippardParams;
 
@@ -50,7 +50,7 @@ class PPippard
     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 SetSecondMomentB(Double_t BB) { fParams.secondMomentB = BB; }
+    virtual void SetVarB(Double_t BB) { fParams.varB = BB; }
 
     virtual void CalculateField();
 
diff --git a/src/tests/nonlocal/nonlocal.cpp b/src/tests/nonlocal/nonlocal.cpp
index d26499a7..0496fe2c 100644
--- a/src/tests/nonlocal/nonlocal.cpp
+++ b/src/tests/nonlocal/nonlocal.cpp
@@ -351,7 +351,7 @@ int main(int argc, char *argv[])
 
   pippard->CalculateField();
 
-  // check if it is necessary to calculate the <B(x)>
+  // check if it is necessary to calculate the <B(x)> and <[B(x)-<B(x)>]^2>
   if (x.size() > 0) {
     if ((params.b_ext == -1.0) || (params.deadLayer == -1.0)) {
       cout << endl << "**ERROR** Bext or deadLayer missing :-(" << endl;
@@ -365,25 +365,29 @@ int main(int argc, char *argv[])
     meanX *= (x[1]-x[0]);
 
     Double_t meanB = 0.0;
-    Double_t secondMomentB = 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];
-        secondMomentB += 1.0 * n[i];
       } else {
         BB = pippard->GetMagneticField(x[i]-params.deadLayer);
         meanB += BB * n[i];
-        secondMomentB += BB * BB * n[i];
       }
     }
     meanB *= (x[1]-x[0]);
-    secondMomentB *= (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]);
 
     cout << endl << ">> mean x = " << meanX << " (nm), mean field = " << params.b_ext * meanB << " (G)";
-    cout << " 2nd Moment B = " << secondMomentB << " (G^2)";
+    cout << " <(B-<B>)^2> = " << varB * params.b_ext * params.b_ext << " (G^2)";
     pippard->SetMeanX(meanX);
-    pippard->SetSecondMomentB(secondMomentB);
+    pippard->SetVarB(varB);
     pippard->SetMeanB(meanB);
   }