← Back to branch summary

~max-mcginley/third-order-lie-map/generator

« back to all changes in this revision

Viewing changes to solenoid_models.cc

  • Committer: Max McGinley
  • Date: 2015-09-10 17:13:14 UTC
  • Revision ID: max.mcginley@stfc.ac.uk-20150910171314-dsy6iqs2u3qhm3u4
More tests added

Show diffs side-by-side

added added

removed removed

Lines of Context:
solenoid_models.cc
 
1
#ifndef SOLENOID_MODELS_CC
 
2
#define SOLENOID_MODELS_CC
 
3
 
 
4
#include <iostream>
 
5
#include "solenoid_models.hh"
 
6
#include <cmath>
 
7
#include "FieldTools/DerivativesSolenoid.hh"
 
8
 
 
9
//*************BUMP FUNCTION based on tanh**********************
 
10
//3 params to pass: lambda, peak bField, and length of solenoid
 
11
 
 
12
double bumpFunction(double z, void* params) {
 
13
        void** args = static_cast<void**>(params);
 
14
        double lambda = *static_cast<double*>(args[0]);
 
15
        double b0 = *static_cast<double*>(args[1]);
 
16
        double L = *static_cast<double*>(args[2]);
 
17
        double origin = *static_cast<double*>(args[3]);
 
18
        if (lambda == 0.0) return (std::abs(z - origin) <= L*0.5) ? b0 : 0.0;
 
19
        return b0*(tanh((z+L*0.5 - origin)/lambda) + tanh((-z+L*0.5+origin)/lambda))/2.0;
 
20
}
 
21
 
 
22
double bumpSecondDer(double z, void* params) {
 
23
        void** args = static_cast<void**>(params);
 
24
        double lambda = *static_cast<double*>(args[0]);
 
25
        double b0 = *static_cast<double*>(args[1]);
 
26
        double L = *static_cast<double*>(args[2]);
 
27
        double origin = *static_cast<double*>(args[3]);
 
28
        z -= origin;
 
29
        if (lambda == 0.0) {
 
30
                std::cerr << "Cannot calculate second derivative of bump function with lambda = 0" <<std::endl;
 
31
                return nan("");
 
32
        }
 
33
        return -b0*(tanh((z+L*0.5)/lambda)/(cosh((z+L*0.5)/lambda)*cosh((z+L*0.5)/lambda)) + tanh((-z+L*0.5)/lambda)/(cosh((-z+L*0.5)/lambda)*cosh((-z+L*0.5)/lambda)))/(lambda*lambda);
 
34
}
 
35
 
 
36
//***************Glaser Model Lorentzian************************
 
37
//2 params to pass: half-width and peak bField
 
38
 
 
39
double glaserFunction(double z, void* params) {
 
40
        void** args = static_cast<void**>(params);
 
41
        double half = *static_cast<double*>(args[0]);
 
42
        double b0 = *static_cast<double*>(args[1]);
 
43
        return b0/(1. + (z*z/(half*half)));
 
44
}
 
45
 
 
46
double glaserSecondDer(double z, void* params) {
 
47
        void** args = static_cast<void**>(params);
 
48
        double half = *static_cast<double*>(args[0]);
 
49
        double b0 = *static_cast<double*>(args[1]);
 
50
        double halfSq = half*half;
 
51
        return 2.0*b0*halfSq*(3.*z*z-halfSq)/pow(halfSq + z*z,3);
 
52
}
 
53
 
 
54
//*****************Unshielded Amperian Model**********************
 
55
//3 params to pass: radius of solenoid, peak bField, and length of solenoid
 
56
 
 
57
double ampFunction(double z, void* params) {
 
58
        void** args = static_cast<void**>(params);
 
59
        double R = *static_cast<double*>(args[0]);
 
60
        double b0 = *static_cast<double*>(args[1]);
 
61
        double L = *static_cast<double*>(args[2]);
 
62
        return b0*(((z+L*0.5)/sqrt(pow(z+L*0.5,2)+R*R)) - ((z-L*0.5)/sqrt(pow(z-L*0.5,2)+R*R)))/2.0;
 
63
}
 
64
 
 
65
double ampSecondDer(double z, void* params) {
 
66
        void** args = static_cast<void**>(params);
 
67
        double R = *static_cast<double*>(args[0]);
 
68
        double b0 = *static_cast<double*>(args[1]);
 
69
        double L = *static_cast<double*>(args[2]);
 
70
        double denom1 = sqrt(pow(z - L*0.5,2) + R*R);
 
71
        double denom2 = sqrt(pow(z + L*0.5,2) + R*R);
 
72
        return 1.5*b0*(z-L*0.5)*(pow(denom1,-3) - pow(z-L*0.5,2)*pow(denom1,-5)) - 1.5*b0*(z+L*0.5)*(pow(denom2,-3) - pow(z+L*0.5,2)*pow(denom2,-5)); 
 
73
}
 
74
 
 
75
//*****************BUMP LATTICE based on tanh*********************
 
76
 
 
77
static const int mausFieldCutoffCellNumber = 3;
 
78
 
 
79
BumpLattice::BumpLattice(int repeats, double cellLength, double lambda, double b0, double solLength) {
 
80
        _numberOfFunctions = repeats*2;
 
81
        _data = new double*[_numberOfFunctions];
 
82
        _fieldGroup = new BTFieldGroup();
 
83
        _params = new void**[_numberOfFunctions];
 
84
        for (int i = 0; i < repeats; i++) {
 
85
                _data[i*2] = getParams(i,0,repeats,cellLength,lambda,b0,solLength);
 
86
                _data[i*2 + 1] = getParams(i,1,repeats,cellLength,lambda,b0,solLength);
 
87
                
 
88
                _params[i*2] = new void*[4];
 
89
                for (int j = 0; j < 4; j++) _params[i*2][j] = static_cast<void*>(_data[i*2] + j);
 
90
                _params[i*2 + 1] = new void*[4];
 
91
                for (int j = 0; j < 4; j++) _params[i*2 + 1][j] = static_cast<void*>(_data[i*2 + 1] + j);
 
92
                
 
93
                BTMultipole::TanhEndField* endFieldPos = new BTMultipole::TanhEndField(solLength/2.,lambda,3);
 
94
                MAUS::DerivativesSolenoid* solPos = new MAUS::DerivativesSolenoid(b0,500.,cellLength*mausFieldCutoffCellNumber,2,endFieldPos);
 
95
                BTMultipole::TanhEndField* endFieldNeg = new BTMultipole::TanhEndField(solLength/2.,lambda,3);
 
96
                MAUS::DerivativesSolenoid* solNeg = new MAUS::DerivativesSolenoid(-b0,500.,cellLength*mausFieldCutoffCellNumber,2,endFieldNeg);
 
97
                _fieldGroup->AddField(solPos,Hep3Vector(0.0,0.0,(i - repeats/2 + 1)*cellLength - cellLength/4.));
 
98
                _fieldGroup->AddField(solNeg,Hep3Vector(0.0,0.0,(i - repeats/2 + 1)*cellLength + cellLength/4.));
 
99
        }
 
100
}
 
101
 
 
102
double* BumpLattice::getParams(int cellNumber, int tanhNumber, int totalRepeats, double cellLength, double lambda, double b0, double solLength) {
 
103
        double* ret = new double[4];
 
104
        ret[0] = lambda;
 
105
        ret[1] = (tanhNumber == 1) ? -b0 : b0;
 
106
        ret[2] = solLength;
 
107
        ret[3] = (cellNumber - totalRepeats/2 + 1)*cellLength + ((tanhNumber == 1) ? cellLength/4. : -cellLength/4.);
 
108
        return ret;
 
109
}
 
110
 
 
111
double BumpLattice::getBField(double z) {
 
112
        double ret = 0.0;
 
113
        for (int i = 0; i < _numberOfFunctions; i++) {
 
114
                ret += bumpFunction(z,_params[i]);
 
115
        }
 
116
        return ret;
 
117
}
 
118
 
 
119
double BumpLattice::getSecondDer(double z) {
 
120
        double ret = 0.0;
 
121
        for (int i = 0; i < _numberOfFunctions; i++) {
 
122
                ret += bumpSecondDer(z,_params[i]);
 
123
        }
 
124
        return ret;
 
125
}
 
126
 
 
127
BumpLattice::~BumpLattice() {
 
128
        for (int i = 0; i < _numberOfFunctions; i++) {
 
129
                delete[] _params[i];
 
130
        }
 
131
        delete[] _params;
 
132
        for (int i = 0; i < _numberOfFunctions; i++) {
 
133
                delete[] _data[i];
 
134
        }
 
135
        delete[] _data;
 
136
}
 
137
 
 
138
#endif