proper implementation of the variance B
This commit is contained in:
@ -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");
|
||||
|
||||
|
@ -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();
|
||||
|
||||
|
@ -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);
|
||||
}
|
||||
|
||||
|
Reference in New Issue
Block a user