[27] | 1 | #ifndef MATHEMATICALTOOLS_SEEN |
---|
| 2 | #define MATHEMATICALTOOLS_SEEN |
---|
| 3 | #include <iostream> |
---|
| 4 | #include <fstream> |
---|
| 5 | #include <sstream> |
---|
| 6 | #include <vector> |
---|
| 7 | #include <cmath> |
---|
| 8 | |
---|
| 9 | |
---|
| 10 | |
---|
| 11 | using namespace std; |
---|
| 12 | |
---|
| 13 | |
---|
| 14 | |
---|
| 15 | |
---|
| 16 | class mathTools |
---|
| 17 | { |
---|
| 18 | |
---|
| 19 | |
---|
| 20 | public : |
---|
| 21 | // return the index K such that : |
---|
| 22 | // vec[K] <= xx < vec[K+1] |
---|
| 23 | // vec is assumed to be increasingly ordered |
---|
| 24 | // if xx is outside , return -1 |
---|
| 25 | static int locateInVector(const vector<double>& vec, double xx) |
---|
| 26 | { |
---|
| 27 | int n= vec.size(); |
---|
| 28 | int jl = 0; |
---|
| 29 | int ju = n+1; |
---|
| 30 | int jm; |
---|
| 31 | if ( vec[0] == xx ) return 0; |
---|
| 32 | if ( vec.back() <= xx || vec[0] > xx) return -1; |
---|
| 33 | while ( ju - jl > 1 ) |
---|
| 34 | { |
---|
| 35 | jm = (ju + jl)/2; |
---|
| 36 | if ( vec[jm-1] <= xx ) jl = jm; |
---|
| 37 | else ju = jm; |
---|
| 38 | } |
---|
| 39 | return jl - 1; |
---|
| 40 | } |
---|
| 41 | |
---|
| 42 | |
---|
| 43 | |
---|
| 44 | |
---|
| 45 | // algo. Boris (Birdsall p. 356) in 2D |
---|
| 46 | // 'tgArcMoitie' is equal to tg(theta/2), if theta is |
---|
| 47 | // the desired rotation angle. |
---|
| 48 | // if tg(theta/2) is given exactly, the rotation is 'exact' |
---|
| 49 | // |
---|
| 50 | static void borisBunemanRotation(double tgArcMoitie, double& vz, double& vx) |
---|
| 51 | { |
---|
| 52 | double t2 = tgArcMoitie * tgArcMoitie; |
---|
| 53 | double sinAngle = 2.0 * tgArcMoitie / ( 1 + t2 ); |
---|
| 54 | double vauxz = vz + vx*tgArcMoitie; |
---|
| 55 | vx = vx - vauxz*sinAngle; |
---|
| 56 | vz = vauxz + vx*tgArcMoitie; |
---|
| 57 | } |
---|
| 58 | |
---|
| 59 | }; |
---|
| 60 | |
---|
| 61 | |
---|
| 62 | |
---|
| 63 | |
---|
| 64 | class TRIDVECTOR |
---|
| 65 | { |
---|
| 66 | double vec_[3]; |
---|
| 67 | |
---|
| 68 | public : |
---|
| 69 | |
---|
| 70 | TRIDVECTOR() |
---|
| 71 | { |
---|
| 72 | int k; |
---|
| 73 | for (k=0; k<3; k++) vec_[k] = 0.0; |
---|
| 74 | } |
---|
| 75 | |
---|
| 76 | TRIDVECTOR( const TRIDVECTOR& triv) |
---|
| 77 | { |
---|
| 78 | setComponents( triv.vec_[0],triv.vec_[1], triv.vec_[2]); |
---|
| 79 | } |
---|
| 80 | |
---|
| 81 | |
---|
| 82 | |
---|
| 83 | TRIDVECTOR(double x,double y,double z) |
---|
| 84 | { |
---|
| 85 | vec_[0] = x; |
---|
| 86 | vec_[1] = y; |
---|
| 87 | vec_[2] = z; |
---|
| 88 | } |
---|
| 89 | |
---|
| 90 | inline double& operator() (int index) { return vec_[index]; } |
---|
| 91 | inline const double& operator() (int index) const { return vec_[index]; } |
---|
| 92 | |
---|
[166] | 93 | inline TRIDVECTOR& operator= (const TRIDVECTOR& triv) |
---|
[27] | 94 | { |
---|
[166] | 95 | setComponents( triv.vec_[0],triv.vec_[1], triv.vec_[2]); |
---|
| 96 | return *this; |
---|
[27] | 97 | } |
---|
| 98 | |
---|
| 99 | |
---|
| 100 | inline TRIDVECTOR& operator = (double value) |
---|
| 101 | { |
---|
| 102 | setComponents( value,value, value); |
---|
[48] | 103 | return *this; |
---|
[27] | 104 | } |
---|
| 105 | |
---|
| 106 | inline void operator *= (const double& factor) |
---|
| 107 | { |
---|
| 108 | int k; |
---|
| 109 | for (k=0; k < 3; k++) vec_[k] *= factor; |
---|
| 110 | // return *this; |
---|
| 111 | } |
---|
| 112 | |
---|
| 113 | |
---|
| 114 | |
---|
| 115 | inline TRIDVECTOR operator + (const TRIDVECTOR& v2) |
---|
| 116 | { |
---|
| 117 | // int k; |
---|
| 118 | return TRIDVECTOR(vec_[0] + v2.vec_[0], vec_[1] + v2.vec_[1], vec_[2] + v2.vec_[2]); |
---|
| 119 | // for (k=0; k < 3; k++) vec_[k] += v2.vec_[k]; |
---|
| 120 | // return *this; |
---|
| 121 | } |
---|
| 122 | |
---|
| 123 | inline void operator += (const TRIDVECTOR& v2) |
---|
| 124 | { |
---|
| 125 | int k; |
---|
| 126 | for (k=0; k < 3; k++) vec_[k] += v2.vec_[k]; |
---|
| 127 | // return *this; |
---|
| 128 | } |
---|
| 129 | |
---|
| 130 | |
---|
| 131 | inline TRIDVECTOR operator * (const double& facteur) const |
---|
| 132 | { |
---|
| 133 | return TRIDVECTOR( vec_[0] * facteur, vec_[1] * facteur, vec_[2] * facteur); |
---|
| 134 | } |
---|
| 135 | |
---|
| 136 | |
---|
| 137 | inline TRIDVECTOR operator * (const TRIDVECTOR& v2) |
---|
| 138 | { |
---|
| 139 | double auxx, auxy, auxz; |
---|
| 140 | // for ( k=0; k<3; k++) aux[k] = vec_[k]; |
---|
| 141 | auxx = vec_[1] * v2.vec_[2] - vec_[2] * v2.vec_[1]; |
---|
| 142 | auxy = vec_[2] * v2.vec_[0] - vec_[0] * v2.vec_[2]; |
---|
| 143 | auxz = vec_[0] * v2.vec_[1] - vec_[1] * v2.vec_[0]; |
---|
| 144 | // return *this; |
---|
| 145 | return TRIDVECTOR(auxx, auxy, auxz); |
---|
| 146 | } |
---|
| 147 | |
---|
| 148 | |
---|
| 149 | |
---|
| 150 | inline double getComponent(int index ) const {return vec_[index];} |
---|
| 151 | |
---|
| 152 | inline void setComponent(int index, double value ) {vec_[index] = value;} |
---|
| 153 | |
---|
| 154 | inline void incrementComponent(int index, double value ) {vec_[index] += value;} |
---|
| 155 | |
---|
| 156 | |
---|
| 157 | inline void getComponents(double& x,double& y,double& z) const |
---|
| 158 | { |
---|
| 159 | x = vec_[0]; |
---|
| 160 | y = vec_[1]; |
---|
| 161 | z = vec_[2]; |
---|
| 162 | } |
---|
| 163 | |
---|
| 164 | inline void setComponents(double x, double y, double z) |
---|
| 165 | { |
---|
| 166 | vec_[0] = x; |
---|
| 167 | vec_[1] = y; |
---|
| 168 | vec_[2] = z; |
---|
| 169 | } |
---|
| 170 | |
---|
| 171 | |
---|
| 172 | |
---|
| 173 | // algo. Boris (Birdsall p. 356) |
---|
| 174 | // the input vector 'tgArcMoitie' is in the direction of the axis of |
---|
| 175 | // rotation. Its module is equal to tg(theta/2), if theta is |
---|
| 176 | // the desired rotation angle. |
---|
| 177 | // if tg(theta/2) is given exactly, the rotation is 'exact' |
---|
| 178 | // |
---|
| 179 | void borisBunemanRotation(const TRIDVECTOR& tgArcMoitie) |
---|
| 180 | { |
---|
| 181 | double t2 = tgArcMoitie.norm2(); |
---|
| 182 | double fac = 2.0 / ( 1.0 + t2); |
---|
| 183 | TRIDVECTOR sinArc = tgArcMoitie * fac; |
---|
| 184 | TRIDVECTOR vecPrim(*this); |
---|
| 185 | vecPrim += (*this) * tgArcMoitie; |
---|
| 186 | (*this) += vecPrim * sinArc; |
---|
| 187 | } |
---|
| 188 | |
---|
| 189 | |
---|
| 190 | |
---|
| 191 | |
---|
| 192 | // rotation in ZX plane |
---|
| 193 | inline void rotateZXComponentsBuneman(double tgAngleMoitie ) |
---|
| 194 | { |
---|
| 195 | mathTools::borisBunemanRotation(tgAngleMoitie, vec_[2], vec_[0]); |
---|
| 196 | } |
---|
| 197 | |
---|
| 198 | // transform the vector Er, Etheta, Ez to a vector Ex,Ey,Ez, assuming the end |
---|
| 199 | // of the vector at (x,y) in the plane (X,Y) |
---|
| 200 | inline void fromCylindricalToCartesian(double costet, double sintet) |
---|
| 201 | { |
---|
| 202 | double auxr = vec_[0]; |
---|
| 203 | double auxtet = vec_[1]; |
---|
| 204 | vec_[0] = auxr * costet - auxtet * sintet; |
---|
| 205 | vec_[1] = auxr * sintet + auxtet * costet; |
---|
| 206 | } |
---|
| 207 | |
---|
| 208 | |
---|
| 209 | inline void clear() |
---|
| 210 | { |
---|
| 211 | vec_[0] = 0.0; |
---|
| 212 | vec_[1] = 0.0; |
---|
| 213 | vec_[2] = 0.0; |
---|
| 214 | } |
---|
| 215 | |
---|
| 216 | inline double norm2() const |
---|
| 217 | { |
---|
| 218 | return vec_[0]*vec_[0] + vec_[1]*vec_[1] + vec_[2]*vec_[2]; |
---|
| 219 | } |
---|
| 220 | inline double norm() const |
---|
| 221 | { |
---|
| 222 | return sqrt(abs(norm2())); |
---|
| 223 | } |
---|
| 224 | |
---|
| 225 | inline void renormalize() |
---|
| 226 | { |
---|
| 227 | int k; |
---|
| 228 | double normeInv = 1.0/norm(); |
---|
| 229 | for (k=0; k< 3 ; k++) vec_[k] *= normeInv; |
---|
| 230 | } |
---|
| 231 | |
---|
| 232 | inline void opposite() |
---|
| 233 | { |
---|
| 234 | int k; |
---|
| 235 | for (k=0; k< 3 ; k++) vec_[k] = -vec_[k]; |
---|
| 236 | } |
---|
| 237 | |
---|
| 238 | inline void print() const |
---|
| 239 | { |
---|
| 240 | cout << " x comp. = " << vec_[0] << " y comp. = " << vec_[1] << " z comp. = " << vec_[2] << endl; |
---|
| 241 | cout << " norme = " << norm() << endl; |
---|
| 242 | } |
---|
| 243 | |
---|
| 244 | string output_flow() const |
---|
| 245 | { |
---|
| 246 | ostringstream sortie; |
---|
| 247 | sortie << " " << vec_[0] << " " << vec_[1] << " " << vec_[2]; |
---|
| 248 | return sortie.str(); |
---|
| 249 | } |
---|
| 250 | |
---|
[52] | 251 | bool input_flow( ifstream& ifs) |
---|
[27] | 252 | { |
---|
| 253 | bool test; |
---|
[52] | 254 | if ( ifs >> vec_[0] >> vec_[1] >> vec_[2]) test= true; |
---|
[27] | 255 | else test = false; |
---|
| 256 | return test; |
---|
| 257 | } |
---|
| 258 | |
---|
| 259 | |
---|
| 260 | }; |
---|
| 261 | |
---|
| 262 | |
---|
| 263 | |
---|
| 264 | |
---|
| 265 | |
---|
| 266 | |
---|
| 267 | |
---|
| 268 | #endif |
---|