source: PSPA/Interface_Web/trunk/pspaWT/sources/controler/include/mathematicalTools.h @ 442

Last change on this file since 442 was 442, checked in by lemeur, 11 years ago

ajout traitement utilisateur

File size: 5.4 KB
Line 
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
11using namespace std;
12
13 
14
15
16class 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
64class 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
76TRIDVECTOR( const TRIDVECTOR& triv) 
77{
78  setComponents( triv.vec_[0],triv.vec_[1], triv.vec_[2]); 
79}
80
81
82
83TRIDVECTOR(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
93inline TRIDVECTOR& operator= (const TRIDVECTOR& triv) 
94{
95  setComponents( triv.vec_[0],triv.vec_[1], triv.vec_[2]);
96  return *this;
97}
98
99
100inline TRIDVECTOR& operator = (double value) 
101{
102  setComponents( value,value, value); 
103  return *this;
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;
145return TRIDVECTOR(auxx, auxy, auxz);
146 }
147
148
149
150inline double getComponent(int index ) const {return vec_[index];}
151
152 inline double& component(int index ) {return vec_[index];}
153
154 inline void setComponent(int index, double value )  {vec_[index] = value;}
155
156 inline void incrementComponent(int index, double value )  {vec_[index] += value;}
157
158
159 inline void getComponents(double& x,double& y,double& z) const
160   {
161     x = vec_[0];
162     y = vec_[1];
163     z = vec_[2];
164   }
165
166 inline void setComponents(double x, double y, double z) 
167   {
168     vec_[0] = x;
169     vec_[1] = y;
170     vec_[2] = z;
171   }
172
173
174
175  // algo. Boris (Birdsall p. 356)
176  // the input vector 'tgArcMoitie' is in the direction of the axis of
177  // rotation. Its module is equal to tg(theta/2), if theta is
178  // the desired rotation angle.
179  // if tg(theta/2) is given exactly, the rotation is 'exact'
180  //
181  void borisBunemanRotation(const TRIDVECTOR& tgArcMoitie)
182  {
183    double t2 = tgArcMoitie.norm2();
184    double fac = 2.0 / ( 1.0 + t2);
185    TRIDVECTOR sinArc = tgArcMoitie * fac;
186    TRIDVECTOR vecPrim(*this);
187    vecPrim += (*this) * tgArcMoitie;
188    (*this) += vecPrim * sinArc;
189  }
190
191
192
193
194// rotation in ZX plane
195inline void rotateZXComponentsBuneman(double tgAngleMoitie )
196{
197  mathTools::borisBunemanRotation(tgAngleMoitie, vec_[2], vec_[0]);
198}
199
200// transform the vector Er, Etheta, Ez to a vector Ex,Ey,Ez, assuming the end
201// of the vector at (x,y) in the plane (X,Y)
202  inline void fromCylindricalToCartesian(double costet, double sintet)
203  {
204     double auxr = vec_[0];
205     double auxtet = vec_[1];
206     vec_[0] = auxr * costet - auxtet * sintet;
207     vec_[1] = auxr * sintet + auxtet * costet;
208  }
209
210
211inline void clear()
212{
213     vec_[0] = 0.0;
214     vec_[1] = 0.0;
215     vec_[2] = 0.0;
216}
217
218inline double norm2() const
219{
220  return vec_[0]*vec_[0] + vec_[1]*vec_[1] + vec_[2]*vec_[2];
221}
222inline double norm() const
223{
224  return sqrt(abs(norm2()));
225}
226
227inline void renormalize()
228{
229  int k;
230  double normeInv = 1.0/norm();
231  for (k=0; k< 3 ; k++) vec_[k] *= normeInv;
232}
233
234inline void opposite()
235{
236  int k;
237  for (k=0; k< 3 ; k++) vec_[k] = -vec_[k];
238}
239
240inline void print() const
241{
242  cout << " x comp. = " << vec_[0] << " y comp. = " << vec_[1] << " z comp. = " << vec_[2] << endl;
243  cout << " norme = " << norm() << endl;
244}
245
246  string output_flow() const
247  {
248    ostringstream sortie;
249    sortie << " "  << vec_[0] << " " <<  vec_[1] << "  " << vec_[2];
250    return sortie.str();
251  }
252
253bool input_flow( ifstream& ifs)
254 {
255    bool test;
256    if ( ifs >> vec_[0] >> vec_[1] >> vec_[2]) test= true;
257    else test =  false;
258    return test;
259 }
260
261
262};
263
264
265
266
267
268
269
270#endif
Note: See TracBrowser for help on using the repository browser.