Kamil Sedlak
fcd5eea567
This is the first version of the muSR simulation code (musrSim) based on the merged codes of Kamil Sedlak and Toni Shiroka. It should be a running version of the simulation code, however it has not been very well tested, therefore it will probably need some further development.
161 lines
5.6 KiB
C++
161 lines
5.6 KiB
C++
// BLEngeFunction.hh
|
|
/*
|
|
This source file is part of G4beamline, http://g4beamline.muonsinc.com
|
|
Copyright (C) 2003,2004,2005,2006 by Tom Roberts, all rights reserved.
|
|
|
|
This program is free software; you can redistribute it and/or
|
|
modify it under the terms of the GNU General Public License
|
|
as published by the Free Software Foundation; either version 2
|
|
of the License, or (at your option) any later version.
|
|
|
|
This program is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
GNU General Public License for more details.
|
|
|
|
http://www.gnu.org/copyleft/gpl.html
|
|
*/
|
|
|
|
#ifndef BLENGEFUNCTION_HH
|
|
#define BLENGEFUNCTION_HH
|
|
|
|
#include <math.h>
|
|
|
|
/** enum BLEngeType specifies the default parameters of BLEngeFunction.
|
|
**/
|
|
enum BLEngeType { ENGE_BLOCK, ENGE_BEND, ENGE_QUAD, ENGE_OTHER };
|
|
|
|
/** class EngeFunction implements the Enge function for fringe fields.
|
|
*
|
|
* z is the distance from the nominal edge, with z=0 being the edge.
|
|
* it should be divided by the aperture diameter or full width/height.
|
|
* z<0 is inside, z>0 is outside.
|
|
*
|
|
* See the COSY reference manual (pp 32-35) for suggested values of
|
|
* a1-a6, or use the ENGE_BEND or ENGE_QUAD types (which come from there).
|
|
* http://cosy.pa.msu.edu/cosymanu/index.html
|
|
*
|
|
* Mathematica was used to compute the derivatives.
|
|
**/
|
|
class BLEngeFunction {
|
|
BLEngeType type;
|
|
double a1,a2,a3,a4,a5,a6;
|
|
public:
|
|
/// default constructor.
|
|
BLEngeFunction() { type=ENGE_BLOCK; set(0,0,0,0,0,0); }
|
|
/// constructor for common magnet types.
|
|
BLEngeFunction(BLEngeType t) {
|
|
switch(t) {
|
|
case ENGE_BLOCK:
|
|
case ENGE_OTHER:
|
|
set(0,0,0,0,0,0);
|
|
break;
|
|
case ENGE_BEND:
|
|
set(0.478959,1.911289,-1.185953,1.630554,-1.082657,0.318111);
|
|
break;
|
|
case ENGE_QUAD:
|
|
set(0.296471,4.533219,-2.270982,1.068627,-0.036391,0.022261);
|
|
break;
|
|
}
|
|
type = t;
|
|
}
|
|
/// general constructor.
|
|
BLEngeFunction(double _a1, double _a2, double _a3, double _a4,
|
|
double _a5, double _a6) {
|
|
set(_a1,_a2,_a3,_a4,_a5,_a6);
|
|
}
|
|
/// set the parameters.
|
|
void set(double _a1, double _a2, double _a3, double _a4,
|
|
double _a5, double _a6) {
|
|
a1=_a1; a2=_a2; a3=_a3; a4=_a4; a5=_a5; a6=_a6;
|
|
if(a1==0.0 && a2==0.0 && a3==0.0 && a4==0.0 && a5==0.0 &&
|
|
a6==0.0)
|
|
type = ENGE_BLOCK;
|
|
else
|
|
type = ENGE_OTHER;
|
|
}
|
|
/// evaluate the Enge function at z.
|
|
double operator()(double z) const {
|
|
if(type == ENGE_BLOCK) return (z<=0.0 ? 1.0 : 0.0);
|
|
if(z < -4.0) return 1.0;
|
|
if(z > 4.0) return 0.0;
|
|
return 1.0/(1.0+exp(a1+z*(a2+z*(a3+z*(a4+z*(a5+z*a6))))));
|
|
}
|
|
/// evaluate the derivative of the Enge function at z.
|
|
double prime(double z) const {
|
|
if(type == ENGE_BLOCK) return 0.0;
|
|
if(fabs(z) > 4.0) return 0.0;
|
|
double exp1 = exp(a1+z*(a2+z*(a3+z*(a4+z*(a5+z*a6)))));
|
|
return -exp1/(1.0+exp1)/(1.0+exp1)*
|
|
(a2+z*(2.0*a3+z*(3.0*a4+z*(4.0*a5+z*5.0*a6))));
|
|
}
|
|
double first(double z) { return prime(z); }
|
|
/// evaluate the second derivative of the Enge function at z.
|
|
double second(double z) const {
|
|
if(type == ENGE_BLOCK) return 0.0;
|
|
if(fabs(z) > 4.0) return 0.0;
|
|
double f1 = a1+z*(a2+z*(a3+z*(a4+z*(a5+z*a6))));
|
|
double f2 = (a2+2*a3*z+3*a4*z*z+4*a5*z*z*z+5*a6*z*z*z*z);
|
|
double f3 = (2*a3+6*a4*z+12*a5*z*z+20*a6*z*z*z);
|
|
double exp1 = exp(f1);
|
|
return exp1*((exp1-1.0)*f2*f2-(1.0+exp1)*f3)/
|
|
(1.0+exp1)/(1.0+exp1)/(1.0+exp1);
|
|
}
|
|
/// evaluate the third derivative of the Enge function at z.
|
|
double third(double z) const {
|
|
if(type == ENGE_BLOCK) return 0.0;
|
|
if(fabs(z) > 4.0) return 0.0;
|
|
double f1 = a1+z*(a2+z*(a3+z*(a4+z*(a5+z*a6))));
|
|
double f2 = a2+z*(2*a3+z*(3*a4+4*a5*z+5*a6*z*z));
|
|
double f3 = 2*(a3+z*(3*a4+2*z*(3*a5+5*a6*z)));
|
|
double f4 = a4+2.0*z*(2.0*a5+5.0*a6*z);
|
|
double exp1 = exp(f1);
|
|
double onepexp1 = 1.0 + exp1;
|
|
return -exp1*(6*exp1*exp1*f2*f2*f2-6*exp1*f2*(f2*f2+f3)*onepexp1
|
|
+(f2*f2*f2+3*f2*f3+6*f4)*onepexp1*onepexp1)
|
|
/(onepexp1*onepexp1*onepexp1*onepexp1);
|
|
}
|
|
/// evaluate the fourth derivative of the Enge function at z.
|
|
double fourth(double z) const {
|
|
if(type == ENGE_BLOCK) return 0.0;
|
|
if(fabs(z) > 4.0) return 0.0;
|
|
double f1 = a1+z*(a2+z*(a3+z*(a4+z*(a5+z*a6))));
|
|
double f2 = a2+z*(2*a3+z*(3*a4+4*a5*z+5*a6*z*z));
|
|
double f3 = 2*(a3+z*(3*a4+2*z*(3*a5+5*a6*z)));
|
|
double f4 = a4+2.0*z*(2.0*a5+5.0*a6*z);
|
|
double f5 = a5 + 5*a6*z;
|
|
double exp1 = exp(f1);
|
|
double onepexp1 = 1.0 + exp1;
|
|
return -exp1*(-24*exp1*exp1*exp1*f2*f2*f2*f2+onepexp1*
|
|
(36*exp1*exp1*f2*f2*(f2*f2+f3)-2*exp1*(7*f2*f2*f2*f2
|
|
+18*f2*f2*f3+3*f3*f3+24*f2*f4)*onepexp1
|
|
+(f2*f2*f2*f2+6*f2*f2*f3+3*f3*f3+24*f2*f4+24*f5)
|
|
*onepexp1*onepexp1))
|
|
/(onepexp1*onepexp1*onepexp1*onepexp1*onepexp1);
|
|
}
|
|
/// evaluate the fifth derivative of the Enge function at z.
|
|
double fifth(double z) const {
|
|
if(type == ENGE_BLOCK) return 0.0;
|
|
if(fabs(z) > 4.0) return 0.0;
|
|
double f1 = a1+z*(a2+z*(a3+z*(a4+z*(a5+z*a6))));
|
|
double f2 = a2+z*(2*a3+z*(3*a4+4*a5*z+5*a6*z*z));
|
|
double f3 = 2*(a3+z*(3*a4+2*z*(3*a5+5*a6*z)));
|
|
double f4 = a4+2.0*z*(2.0*a5+5.0*a6*z);
|
|
double f5 = a5 + 5*a6*z;
|
|
double exp1 = exp(f1);
|
|
double onepexp1 = 1.0 + exp1;
|
|
return -exp1/(onepexp1*onepexp1*onepexp1*onepexp1*onepexp1*onepexp1)
|
|
*(120*exp1*exp1*exp1*exp1*f2*f2*f2*f2*f2
|
|
-240*exp1*exp1*exp1*f2*f2*f2*(f2*f2+f3)*onepexp1
|
|
+onepexp1*onepexp1*(30*exp1*exp1*f2*(5*f2*f2*f2*f2
|
|
+12*f2*f2*f3+3*f3*f3+12*f2*f4)-10*exp1*(3*f2*f2*f2*f2*f2
|
|
+14*f2*f2*f2*f3+9*f2*f3*f3+36*f2*f2*f4+12*f3*f4+24*f2*f5)
|
|
*onepexp1+(120*a6+f2*f2*f2*f2*f2+10*f2*f2*f2*f3+60*f2*f2*f4
|
|
+60*f3*f4+15*f2*(f3*f3+8*f5))*onepexp1*onepexp1));
|
|
}
|
|
/// return the type of Enge function
|
|
BLEngeType getType() const { return type; }
|
|
};
|
|
|
|
#endif // BLENGEFUNCTION_HH
|