← Back to branch summary

~t-carlisle1/maus/test

« back to all changes in this revision

Viewing changes to tests/integration/optics/src/Tensor3.cc

  • Committer: Timothy Carlisle (MICE 2008)
  • Date: 2011-06-22 16:59:28 UTC
  • mfrom: (486.3.56 maus)
  • Revision ID: carlisle@pplxint6.physics.ox.ac.uk-20110622165928-4mr32fa3anlw2odu
bleu

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/integration/optics/src/Tensor3.cc
1
 
#include "Tensor3.hh"
2
 
 
3
 
Tensor3::Tensor3(int a, int b, int c, double val) : Tensor(0, 9), _a(a), _b(b), _c(c)
4
 
{
5
 
        for(int i=0; i<c; i++)
6
 
        {
7
 
                _data.push_back(CLHEP::HepMatrix(a,b,0));
8
 
                if(i<a && i<b && i<c)
9
 
                        _data[i][i][i] = val;
10
 
        }
11
 
}
12
 
 
13
 
Tensor3&  Tensor3::operator*=(double t)
14
 
{
15
 
        for(int ia=0; ia<_a; ia++)
16
 
                for(int ib=0; ib<_b; ib++)
17
 
                        for(int ic=0; ic<_c; ic++)
18
 
                                _data[ia][ib][ic]*=t;
19
 
        return *this;
20
 
}
21
 
 
22
 
Tensor3&  Tensor3::operator/=(double t)
23
 
{
24
 
        for(int ia=0; ia<_a; ia++)
25
 
                for(int ib=0; ib<_b; ib++)
26
 
                        for(int ic=0; ic<_c; ic++)
27
 
                                _data[ia][ib][ic]/=t;
28
 
        return *this;
29
 
}
30
 
 
31
 
Tensor3&  Tensor3::operator+=(double t)
32
 
{
33
 
        for(int ia=0; ia<_a; ia++)
34
 
                for(int ib=0; ib<_b; ib++)
35
 
                        for(int ic=0; ic<_c; ic++)
36
 
                                _data[ia][ib][ic]+=t;
37
 
        return *this;
38
 
}
39
 
 
40
 
Tensor3&  Tensor3::operator-=(double t)
41
 
{
42
 
        for(int ia=0; ia<_a; ia++)
43
 
                for(int ib=0; ib<_b; ib++)
44
 
                        for(int ic=0; ic<_c; ic++)
45
 
                                _data[ia][ib][ic]-=t;
46
 
        return *this;
47
 
}
48
 
 
49
 
Tensor3  Tensor3::operator-() const
50
 
{
51
 
        Tensor3 out(*this);
52
 
        out *= -1;
53
 
        return out;
54
 
}
55
 
 
56
 
Tensor3 operator/(const Tensor3& t1,    const double d1)
57
 
{
58
 
        Tensor3          out  = t1;
59
 
        std::vector<int> size = out.Size();
60
 
        for(int ia=0; ia<size[0]; ia++)
61
 
                for(int ib=0; ib<size[1]; ib++)
62
 
                        for(int ic=0; ic<size[2]; ic++)
63
 
                                out[ia][ib][ic]/=d1;
64
 
        return out;
65
 
}
66
 
 
67
 
CLHEP::HepVector Tensor3::Poisson(CLHEP::HepVector v) const
68
 
{
69
 
        //find :T:v
70
 
        //loop over upper diagonal of \mathbf{T}, add the partial derivative wrt \vec{v}.
71
 
        int    dimension = v.num_row();
72
 
        CLHEP::HepVector vOut(dimension);
73
 
        for(int i=0; i<dimension; i++)
74
 
                for(int j=0; j<dimension; j++)
75
 
                        vOut[i] += dGdUi(j, v)*Delta(i, GetConjugate(j))*pow(-1, j+1);
76
 
        return vOut;
77
 
}
78
 
 
79
 
std::vector<Tensor*> Tensor3::Exponential(int maxOrder) const
80
 
{
81
 
//:T_i:^n is of order n*(i-2)+1
82
 
//so :T_3:^n is of order n+1
83
 
        int     maxPower  = maxOrder-1;
84
 
        int     nFact     = 1; //n!
85
 
        std::vector<Tensor*> out; // sum(:t:^n)
86
 
        out.push_back(this->Clone());
87
 
        for(int n=2; n<=maxPower; n++)
88
 
        {
89
 
                nFact          *= n;
90
 
                Tensor* current = out.back(); 
91
 
                out.push_back(new Tensor(current->Poisson(*this)/nFact));
92
 
        }
93
 
        return out;
94
 
}
95
 
 
96
 
CLHEP::HepVector Tensor3::ExponentialPoisson(CLHEP::HepVector v, int maxOrder) const
97
 
{
98
 
        std::vector<Tensor*> ExpT = Exponential(maxOrder);
99
 
        CLHEP::HepVector out;
100
 
        for(unsigned int i=0; i<ExpT.size(); i++)
101
 
                out += ExpT[i]->Poisson(v);
102
 
        return out;
103
 
}
104
 
 
105
 
double Tensor3::dGdUi(int i, CLHEP::HepVector v) const
106
 
{
107
 
        double dGdUi     = 0;
108
 
        int    dimension = v.num_row();
109
 
        for(int p=0; p<dimension; p++)
110
 
                for(int q=p; q<dimension; q++)
111
 
                        for(int r=q; r<dimension; r++)
112
 
                                dGdUi += _data[p][q][r]*(v[q]*v[r]*Delta(i,p) + v[p]*v[r]*Delta(i,q) + v[p]*v[q]*Delta(i,r));
113
 
        return dGdUi;
114
 
}
115
 
 
116
 
void Tensor3::Print(std::ostream& out) const
117
 
{
118
 
        for(int i=0; i<_c; i++)
119
 
        out << _data[i] << "\n";
120
 
}
121
 
 
122
 
bool Tensor3::IsZero() const
123
 
{
124
 
        for(int i=0; i< _a; i++)
125
 
                for(int j=0; i<_b; j++)
126
 
                        for(int k=0; k<_c; k++)
127
 
                                if(_data[i][j][k] > 1e-10 || _data[i][j][k] < -1e-10) return false;
128
 
        return true;
129
 
}
130
 
 
131
 
Tensor3 operator+(const Tensor3& t1, const Tensor3& t2)
132
 
{
133
 
        Tensor3 t3(t1);
134
 
        std::vector<int> size = t3.Size();
135
 
        for(int ia=0; ia<size[0]; ia++)
136
 
                for(int ib=0; ib<size[1]; ib++)
137
 
                        for(int ic=0; ic<size[2]; ic++)
138
 
                                t3[ia][ib][ic] = t1[ia][ib][ic]+t2[ia][ib][ic];
139
 
        return t3;
140
 
}
141
 
 
142
 
Tensor3 operator-(const Tensor3& t1, const Tensor3& t2)
143
 
{
144
 
        Tensor3 t3(t1);
145
 
        std::vector<int> size = t3.Size();
146
 
        for(int ia=0; ia<size[0]; ia++)
147
 
                for(int ib=0; ib<size[1]; ib++)
148
 
                        for(int ic=0; ic<size[2]; ic++)
149
 
                                t3[ia][ib][ic] = t1[ia][ib][ic]-t2[ia][ib][ic];
150
 
        return t3;
151
 
}
152
 
 
153
 
std::ostream& operator<<(std::ostream& out, const Tensor3& t1)
154
 
{
155
 
        t1.Print(out);
156
 
        return out;
157
 
}
158
 
 
159
 
void Tensor3::SetAllFields(double value)
160
 
{
161
 
        std::vector<int> size = Size();
162
 
        for(int ia=0; ia<size[0]; ia++)
163
 
                for(int ib=0; ib<size[1]; ib++)
164
 
                        for(int ic=0; ic<size[2]; ic++)
165
 
                                (*this)[ia][ib][ic] = value;
166
 
}
167
 
 
168
 
std::vector<int> Tensor3::Size() const
169
 
{
170
 
        std::vector<int> size(3);
171
 
        size[0] = _a;
172
 
        size[1] = _b;
173
 
        size[2] = _c;
174
 
        return size;
175
 
}