← Back to branch summary

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

« back to all changes in this revision

Viewing changes to EmittanceCalculator.cc

  • Committer: Max McGinley
  • Date: 2015-09-07 10:20:10 UTC
  • Revision ID: max.mcginley@stfc.ac.uk-20150907102010-0bot9845lf546xwc
Added a AnalyticLattice class in solenoid_models which uses an equation which has an equivalent MAUS::BTField. Also added function for tracking progress of particles round unit cell

Show diffs side-by-side

added added

removed removed

Lines of Context:
EmittanceCalculator.cc
2
2
#define EMITTANCE_CALCULATOR_CC
3
3
 
4
4
#include "EmittanceCalculator.hh"
 
5
#include "max_globals.hh"
5
6
 
6
7
namespace MAUS {
7
8
 
8
 
Matrix<double> EmittanceCalculator::getEmittanceMatrixFromMap(PolynomialMap map) const {
 
9
std::vector<double> EmittanceCalculator::getEmittancesFromMaps(std::vector<PolynomialMap> maps) {
 
10
        std::vector<double> ret(maps.size());
 
11
        for (int i = 0; i < maps.size(); i++) {
 
12
                ret[i] = getEmittanceFromMap(maps[i]);
 
13
        }
 
14
        return ret;
 
15
}
 
16
 
 
17
EmittanceCalculatorDistribution::EmittanceCalculatorDistribution(double emit, DistributionFunction dist) : EmittanceCalculator(emit) {
 
18
        switch(dist) {
 
19
                case distributionUniform:
 
20
                        _beam = BeamDistribution::generateUniformBeam(emit,1.0);
 
21
                        break;
 
22
                default:
 
23
                        throw std::runtime_error("EmittanceCalculatorDistribution::EmittanceCalculatorDistribution() - beam distribution not recognized");
 
24
        }
 
25
} 
 
26
 
 
27
EmittanceCalculatorDistribution::~EmittanceCalculatorDistribution() {
 
28
        delete _beam;
 
29
}
 
30
 
 
31
Matrix<double> EmittanceCalculatorDistribution::getEmittanceMatrixFromMap(PolynomialMap map) const {
9
32
        std::vector<PolynomialMap::PolynomialCoefficient> coeffs = map.GetCoefficientsAsVector();
10
33
        std::vector<PolynomialMap::PolynomialCoefficient> thirdCoeffs;
11
34
        Matrix<double> firstOrder(4,4,0.0);
21
44
        return _beam->getSigmaMatrixAfterTransform(firstOrder,thirdOrder);
22
45
}
23
46
 
 
47
Matrix<double> EmittanceCalculator::getEllipseFromMap(const Matrix<double>& mat) {
 
48
        std::pair<Vector<complex>, Matrix<complex> > eig = eigensystem(mat);
 
49
        Matrix<double> reOrtho = real(eig.second);
 
50
        Matrix<double> imOrtho = imag(eig.second);
 
51
        Vector<double> reEig = real(eig.first);
 
52
        Vector<double> imEig = imag(eig.first);
 
53
        for (int i = 1; i <= 4; i++)
 
54
                for (int j = 1; j <= 4; j++)
 
55
                        if (imOrtho(i,j) != 0.0) {
 
56
                                std::cerr << "Eigenvalues: " << eig.first << std::endl;
 
57
                                std::cerr << "Eigenvectors: " << eig.second << std::endl;
 
58
                                throw std::runtime_error("EmittanceCalculator::getEllipseFromMap - eigenvectors are not pure real");
 
59
                        }
 
60
        for (int i = 1; i <= 4; i++)
 
61
                if (imEig(i) != 0.0) throw std::runtime_error("EmittanceCalculator::getEllipseFromMap - eigenvalues are not pure real");
 
62
        Matrix<double> diag = transpose(reOrtho) * mat * reOrtho;
 
63
        std::cerr << "Diagonal: " << diag << std::endl;
 
64
        return diag;
 
65
}
 
66
 
 
67
//**************************Discrete******************************8
 
68
 
 
69
//TODO: Make this actually produce a beam with desired emittance
 
70
EmittanceCalculatorDiscrete::EmittanceCalculatorDiscrete(double emit, double posMomRatio, int noParticles) : EmittanceCalculator(emit), _particleNum(noParticles), _particles(new double*[noParticles]) {
 
71
        for (int i = 0; i < _particleNum; i++) {
 
72
                _particles[i] = new double[4];
 
73
                double scale = std::sqrt(6. * emit);
 
74
                double x1 = 1.0, x2 = 1.0, x3 = 1.0, x4 = 1.0;
 
75
                do {
 
76
                        x1 = static_cast <double> ((rand()) / (static_cast <float> (RAND_MAX)) - 0.5) * 2.0;
 
77
                        x2 = static_cast <double> ((rand()) / (static_cast <float> (RAND_MAX)) - 0.5) * 2.0;
 
78
                        x3 = static_cast <double> ((rand()) / (static_cast <float> (RAND_MAX)) - 0.5) * 2.0;
 
79
                        x4 = static_cast <double> ((rand()) / (static_cast <float> (RAND_MAX)) - 0.5) * 2.0;
 
80
                } while(x1*x1 + x2*x2 +x3*x3 +x4*x4 > 1.0);
 
81
                _particles[i][0] = x1 * scale * posMomRatio / sigma;
 
82
                _particles[i][1] = x2 * scale * posMomRatio / sigma;
 
83
                _particles[i][2] = x3 * scale / (posMomRatio*refMomentum);
 
84
                _particles[i][3] = x4 * scale / (posMomRatio*refMomentum);
 
85
        }
 
86
}
 
87
 
 
88
EmittanceCalculatorDiscrete::EmittanceCalculatorDiscrete(double maxPos, double maxMom, int noParticles, bool ellipse) : EmittanceCalculator(std::sqrt(maxPos*maxMom/6.)), _particleNum(noParticles), _particles(new double*[noParticles])  {
 
89
        for (int i = 0; i < _particleNum; i++) {
 
90
                _particles[i] = new double[4];
 
91
                double x1 = 1.0, x2 = 1.0, x3 = 1.0, x4 = 1.0;
 
92
                do {
 
93
                        x1 = static_cast <double> ((rand()) / (static_cast <float> (RAND_MAX)) - 0.5) * 2.0;
 
94
                        x2 = static_cast <double> ((rand()) / (static_cast <float> (RAND_MAX)) - 0.5) * 2.0;
 
95
                        x3 = static_cast <double> ((rand()) / (static_cast <float> (RAND_MAX)) - 0.5) * 2.0;
 
96
                        x4 = static_cast <double> ((rand()) / (static_cast <float> (RAND_MAX)) - 0.5) * 2.0;
 
97
                } while(x1*x1 + x2*x2 +x3*x3 +x4*x4 > 1.0 && ellipse);
 
98
                _particles[i][0] = x1 * maxPos / sigma;
 
99
                _particles[i][1] = x2 * maxPos / sigma;
 
100
                _particles[i][2] = x3 * maxMom / (refMomentum);
 
101
                _particles[i][3] = x4 * maxMom / (refMomentum);
 
102
        }
 
103
}
 
104
 
 
105
EmittanceCalculatorDiscrete::~EmittanceCalculatorDiscrete() {
 
106
        for (int i = 0; i < _particleNum; i++) {
 
107
                delete[] _particles[i];
 
108
        }
 
109
        delete _particles;
 
110
}
 
111
 
 
112
Matrix<double> EmittanceCalculatorDiscrete::getEmittanceFromParticles(double const * const * const particles, int number) {
 
113
        double* firstMoments = new double[4];
 
114
        Matrix<double> sigmas(4,4,0.0);
 
115
        for (int p = 0; p < number; p++) {
 
116
                for (int i = 0; i < 4; i++) {
 
117
                        firstMoments[i] += particles[p][i]/number;
 
118
                }
 
119
        }
 
120
        for (int p = 0; p < number; p++) {
 
121
                for (int i = 0; i < 4; i++) {
 
122
                        for (int j = 0; j < 4; j++) {
 
123
                                sigmas(i+1,j+1) += (particles[p][i] - firstMoments[i])*(particles[p][j] - firstMoments[j])/number;
 
124
                        }
 
125
                }
 
126
        }
 
127
        delete[] firstMoments;
 
128
        return sigmas;
 
129
}
 
130
 
 
131
Matrix<double> EmittanceCalculatorDiscrete::getEmittanceMatrixFromMap(PolynomialMap map) const {
 
132
        double ** out = new double*[_particleNum];
 
133
        for (int i = 0; i < _particleNum; i++) {
 
134
                out[i] = new double[4];
 
135
                map.F(_particles[i],out[i]);
 
136
        }
 
137
        Matrix<double> sig = getEmittanceFromParticles(out,_particleNum);
 
138
        for (int i = 0; i < _particleNum; i++)
 
139
                delete[] out[i];
 
140
        delete[] out;
 
141
        return sig;
 
142
}
24
143
 
25
144
}// namespace MAUS
26
145