O P A R - Open Architecture Particle in Cell Simulation - Version 3.0
Plasma simulations with dust particles
 All Classes Files Functions Variables Friends Macros Groups Pages
dimensionality.h
Go to the documentation of this file.
1 
6 #ifndef DIMENSIONALITY_H
7 #define DIMENSIONALITY_H
8 #include "normalization.h"
9 #include <iostream>
10 #include <istream>
11 #include <ostream>
12 // The dimensionality of the simulation
13 //#define ONE_DIMENSIONAL
14 //#define TWO_DIMENSIONAL
15 #define THREE_DIMENSIONAL
16 #define OMIT_CODE
17 //---------------------------------------------------------------------------------------------------------------------
18 template<class T> class Position_;
19 // Position is an n-dimensional 'vector' of double
20 typedef Position_<double> Position;
21 // GridPosition is an n-dimensional 'vector' of int
22 typedef Position_<int> GridPosition;
23 //---------------------------------------------------------------------------------------------------------------------
24 #ifdef ONE_DIMENSIONAL
25 const int Dimension = 1;
26 const int TwoPowDim = 2;
27 template<class T>
28 class Position_ {
29  public:
30  T P[1];
31  T& x() {return P[0];};
32  Position_() {}
33  Position_(T ax) {P[0]=ax;}
34  operator T*() {return P;}
35  T volume() const {return P[0];}
36  T sum() const {return P[0];}
37  double length() const {return fabs(P[0]);}
38  void stretch(const Position_<T>& A) {P[0]*=A.P[0];}
39  void limitto(const Position_<T>& pos) {if (P[0] > pos.P[0]) P[0] = pos.P[0];}
40  T arraypos(const Position_<T>& pos) const {return P[0];}
41  template<class T2>
42  operator Position_<T2>() const {return Position_<T2>(T2(P[0]));}
43  T& operator[](const int i) {return P[i];}
44  T operator[](const int i) const {return P[i];}
45 };
46 template<class T>
47 Position_<T> operator-(const Position_<T>& A, T b) {return Position_<T>(A.P[0]-b);}
48 template<class T>
49 Position_<T> operator-(T a, const Position_<T>& B) {return Position_<T>(a-B.P[0]);}
50 template<class T>
51 Position_<T> operator+(const Position_<T>& A, T b) {return Position_<T>(A.P[0]+b);}
52 template<class T>
53 Position_<T> operator-(const Position_<T>& A, const Position_<T>& B) {return Position_<T>(A.P[0]-B.P[0]);}
54 template<class T>
55 Position_<T> operator+(const Position_<T>& A, const Position_<T>& B) {return Position_<T>(A.P[0]+B.P[0]);}
56 template<class T1, class T2>
57 double operator*(const Position_<T1>& A, const Position_<T2>& B) {return A.P[0]*B.P[0];}
58 template<class T1, class T2>
59 Position_<T1>& operator*=(Position_<T1>& A, T2 B) {
60  A.P[0] *= B;
61  return A;
62 }
63 template<class T>
64 Position_<T>& operator*=(Position_<T>& A, const Position_<T>& B) {
65  A.P[0] *= B.P[0];
66  return A;
67 }
68 template<class T1, class T2>
69  Position_<T1>& operator+=(Position_<T1>& A, const Position_<T2>& B) {
70  A.P[0] += B.P[0];
71  return A;
72 }
73 template<class T>
74 Position_<T> operator*(const Position_<T>& A, T B) {return Position_<T>(A.P[0]*B);}
75 template<class T>
76 Position_<T> operator/(const Position_<T>& A, T B) {return Position_<T>(A.P[0]/B);}
77 template<class T>
78 Position_<T> operator/(const Position_<T>& A, const Position_<T>& B) {return Position_<T>(A.P[0]/B.P[0]);}
79 template<class T1, class T2>
80 bool operator<(const Position_<T1>& A, const Position_<T2>& B) {return A.P[0]<B.P[0];}
81 template<class T1, class T2>
82 bool operator>(const Position_<T1>& A, const Position_<T2>& B) {return A.P[0]>B.P[0];}
83 template<class T1, class T2>
84 bool operator==(const Position_<T1>& A, const Position_<T2>& B) {return A.P[0]==B.P[0];}
85 template<class T>
86 std::ostream& operator<<(std::ostream &O, const Position_<T>& A) {return O << A.P[0];}
87 template<class T>
88 std::istream& operator>>(std::istream &I, Position_<T>& A) {return I >> A.P[0];}
89 template<class T1, class T2>
90 int indexgreater(const Position_<T1>& A, const Position_<T2>& B) {if (A.P[0]>=B.P[0]) return 0;}
91 template<class G>
92 void WeightOnGrid(G& g, const Position& p, double v) {
93  GridPosition gp(p);
94  Position Diff = p - Position(gp);
95  Position diff = Position(1) - Diff;
96  g[gp] += diff.P[0]*v;
97  ++gp.P[0];
98  g[gp] += Diff.P[0]*v;
99 }
100 const GridPosition ZeroGrid(0);
101 const GridPosition UnitGrid(1);
102 #endif
103 //---------------------------------------------------------------------------------------------------------------------
104 #ifdef TWO_DIMENSIONAL
105 const int Dimension = 2;
106 const int TwoPowDim = 4;
107 template<class T>
108 class Position_ {
109  public:
111  T P[2];
113  T& x() {return P[0];};
115  T& y() {return P[1];};
117  const T& x() const {return P[0];};
119  const T& y() const {return P[1];};
121  Position_() {}
123  Position_(T ax, T ay) {P[0]=ax; P[1]=ay;}
125  Position_(T x) {P[0]=x; P[1]=x;}
127  operator T*() {return P;}
129  T volume() const {return P[0]*P[1];}
131  T sum() const {return P[0]+P[1];}
133  double length() const {return sqrt(P[0]*P[0]+P[1]*P[1]);}
135  void stretch(const Position_<T>& A) {P[0] *= A.P[0]; P[1] *= A.P[1];}
136 
140  void limitto(const Position_<T>& pos) {
141  if (P[0]>pos.P[0]) P[0]=pos.P[0];
142  if (P[1]>pos.P[1]) P[1]=pos.P[1];
143  }
144 
149  T arraypos(const Position_<T>& pos) const {return P[0] + pos.P[0]*P[1];}
150  // Typecast into a Position_<T2> where T2 is different from T
151  template<class T2>
152  operator Position_<T2>() const {return Position_<T2>(T2(P[0]), T2(P[1]));}
154  T& operator[](const int i) {return P[i];}
155  T operator[](const int i) const {return P[i];}
157 
158 };
159 
177 // @see PosOperators
178 template<class T>
179 Position_<T> operator-(const Position_<T>& A, T b) {return Position_<T>(A.P[0]-b, A.P[1]-b);}
181 template<class T>
182 Position_<T> operator-(T a, const Position_<T>& B) {return Position_<T>(a-B.P[0], a-B.P[1]);}
184 template<class T>
185 Position_<T> operator+(const Position_<T>& A, T b) {return Position_<T>(A.P[0]+b, A.P[1]+b);}
187 template<class T>
188 Position_<T> operator-(const Position_<T>& A, const Position_<T>& B) {
189  return Position_<T>(A.P[0]-B.P[0], A.P[1]-B.P[1]);
190 }
192 template<class T>
193 Position_<T> operator+(const Position_<T>& A, const Position_<T>& B) {
194  return Position_<T>(A.P[0]+B.P[0], A.P[1]+B.P[1]);
195 }
197 template<class T1, class T2>
198 double operator*(const Position_<T1>& A, const Position_<T2>& B) {return A.P[0]*B.P[0] + A.P[1]*B.P[1];}
200 template<class T1, class T2>
201  Position_<T1>& operator*=(Position_<T1>& A, T2 B) {
202  A.P[0] *= B;
203  A.P[1] *= B;
204  return A;
205 }
207 template<class T>
208 Position_<T>& operator*=(Position_<T>& A, const Position_<T>& B) {
209  A.P[0] *= B.P[0];
210  A.P[1] *= B.P[1];
211  return A;
212 }
214 template<class T1, class T2>
215  Position_<T1>& operator*=(Position_<T1>& A, const Position_<T2>& B) {
216  A.P[0] *= B.P[0];
217  A.P[1] *= B.P[1];
218  return A;
219 }
221 template<class T1, class T2>
222  Position_<T1>& operator+=(Position_<T1>& A, const Position_<T2>& B) {
223  A.P[0] += B.P[0];
224  A.P[1] += B.P[1];
225  return A;
226 }
228 template<class T>
229 Position_<T> operator*(const Position_<T>& A, T B) {return Position_<T>(A.P[0]*B, A.P[1]*B);}
231 template<class T>
232 Position_<T> operator/(const Position_<T>& A, T B) {return Position_<T>(A.P[0]/B, A.P[1]/B);}
234 template<class T>
235 Position_<T> operator/(const Position_<T>& A, const Position_<T>& B) {
236  return Position_<T>(A.P[0]/B.P[0], A.P[1]/B.P[1]);
237 }
239 template<class T1, class T2>
240 bool operator<(const Position_<T1>& A, const Position_<T2>& B) {return (A.P[0]<B.P[0]) || (A.P[1]<B.P[1]);}
242 template<class T1, class T2>
243 bool operator>(const Position_<T1>& A, const Position_<T2>& B) {return (A.P[0]>B.P[0]) || (A.P[1]>B.P[1]);}
245 template<class T1, class T2>
246 bool operator==(const Position_<T1>& A, const Position_<T2>& B) {return (A.P[0]==B.P[0]) && (A.P[1]==B.P[1]);}
248 
249 template<class T>
250 std::ostream& operator<<(std::ostream &O, const Position_<T>& A) {return O << A.P[0] << ' ' << A.P[1];}
252 template<class T>
253 std::istream& operator>>(std::istream &I, Position_<T>& A) {return I >> A.P[0] >> A.P[1];}
254 
259 template<class T1, class T2>
260 int indexgreater(const Position_<T1>& A, const Position_<T2>& B) {
261  if (A.P[0]>=B.P[0]) return 0;
262  if (A.P[1]>=B.P[1]) return 1;
263 }
264 
269 template<class G>
270 void WeightOnGrid(G& g, const Position& p, double v) {
271  Position j1 = Position(1,1);
272  GridPosition gp(p);
273  Position Diff = p - Position(gp);
274  Position diff = Position(1,1) - Diff;
275  g[gp] += diff[0]*diff[1]*v;
276  ++gp.P[0];
277  g[gp] += Diff[0]*diff[1]*v;
278  ++gp.P[1];
279  g[gp] += Diff[0]*Diff[1]*v;
280  --gp.P[0];
281  g[gp] += diff[0]*Diff[1]*v;
282 }
284 const GridPosition ZeroGrid(0,0);
285 const GridPosition UnitGrid(1,1);
286 const GridPosition UnitX(1,0);
287 const GridPosition UnitY(0,1);
288 const GridPosition UnitnX(-1,0);
289 const GridPosition UnitnY(0,-1);
291 //@]
292 #endif
293 //---------------------------------------------------------------------------------------------------------------------
294 #ifdef THREE_DIMENSIONAL
295 const int Dimension = 3;
296 const int TwoPowDim = 8;
297 // Position_ holds the template of an n-dimensional 'vector'
298 template<class T>
299 class Position_ {
300  public:
302  T P[3];
304  T& x() {return P[0];};
306  T& y() {return P[1];};
308  T& z() {return P[2];};
310  const T& x() const {return P[0];};
312  const T& y() const {return P[1];};
314  const T& z() const {return P[2];};
316  Position_() {}
318  Position_(T ax, T ay, T az) {P[0]=ax; P[1]=ay; P[2]=az;}
320  Position_(T x) {P[0]=x; P[1]=x; P[2]=x;}
322  operator T*() {return P;}
324  T volume() const {return P[0]*P[1]*P[2];}
326  T sum() const {return P[0]+P[1]+P[2];}
328  double length() const {return sqrt(P[0]*P[0]+P[1]*P[1]+P[2]*P[2]);}
330  void stretch(const Position_<T>& A) {P[0] *= A.P[0]; P[1] *= A.P[1]; P[2] *= A.P[2];}
331 
335  void limitto(const Position_<T>& pos) {
336  if (P[0]>pos.P[0]) P[0]=pos.P[0];
337  if (P[1]>pos.P[1]) P[1]=pos.P[1];
338  if (P[2]>pos.P[2]) P[2]=pos.P[2];
339  }
340 
345  T arraypos(const Position_<T>& pos) const {return P[0] + (pos.P[0])*(P[1] + (pos.P[1])*P[2]);}
346  // Typecast into a Position_<T2> where T3 is different from T
347  template<class T3>
348  operator Position_<T3>() const {return Position_<T3>(T3(P[0]), T3(P[1]), T3(P[2]));}
350  T& operator[](const int i) {return P[i];}
351  T operator[](const int i) const {return P[i];}
353 
354 };
355 
373 // @see PosOperators
374 template<class T>
375 Position_<T> operator-(const Position_<T>& A, T b) {return Position_<T>(A.P[0]-b, A.P[1]-b, A.P[2]-b);}
377 template<class T>
378 Position_<T> operator-(T a, const Position_<T>& B) {return Position_<T>(a-B.P[0], a-B.P[1], a-B.P[2]);}
380 template<class T>
381 Position_<T> operator+(const Position_<T>& A, T b) {return Position_<T>(A.P[0]+b, A.P[1]+b, A.P[2]+b);}
383 template<class T>
384 Position_<T> operator-(const Position_<T>& A, const Position_<T>& B) {
385  return Position_<T>(A.P[0]-B.P[0], A.P[1]-B.P[1], A.P[2]-B.P[2]);
386 }
388 template<class T>
389 Position_<T> operator+(const Position_<T>& A, const Position_<T>& B) {
390  return Position_<T>(A.P[0]+B.P[0], A.P[1]+B.P[1], A.P[2]+B.P[2]);
391 }
393 template<class T1, class T2>
394 double operator*(const Position_<T1>& A, const Position_<T2>& B) {return A.P[0]*B.P[0] + A.P[1]*B.P[1] + A.P[2]*B.P[2];}
396 template<class T1, class T2>
397  Position_<T1>& operator*=(Position_<T1>& A, T2 B) {
398  A.P[0] *= B;
399  A.P[1] *= B;
400  A.P[2] *= B;
401  return A;
402 }
404 template<class T>
405 Position_<T>& operator*=(Position_<T>& A, const Position_<T>& B) {
406  A.P[0] *= B.P[0];
407  A.P[1] *= B.P[1];
408  A.P[2] *= B.P[2];
409  return A;
410 }
412 template<class T1, class T2>
413  Position_<T1>& operator*=(Position_<T1>& A, const Position_<T2>& B) {
414  A.P[0] *= B.P[0];
415  A.P[1] *= B.P[1];
416  A.P[2] *= B.P[2];
417  return A;
418 }
420 template<class T1, class T2>
421  Position_<T1>& operator+=(Position_<T1>& A, const Position_<T2>& B) {
422  A.P[0] += B.P[0];
423  A.P[1] += B.P[1];
424  A.P[2] += B.P[2];
425  return A;
426 }
428 template<class T>
429 Position_<T> operator*(const Position_<T>& A, T B) {return Position_<T>(A.P[0]*B, A.P[1]*B, A.P[2]*B);}
431 template<class T>
432 Position_<T> operator/(const Position_<T>& A, T B) {return Position_<T>(A.P[0]/B, A.P[1]/B, A.P[2]/B);}
434 template<class T>
435 Position_<T> operator/(const Position_<T>& A, const Position_<T>& B) {
436  return Position_<T>(A.P[0]/B.P[0], A.P[1]/B.P[1], A.P[2]/B.P[2]);
437 }
439 template<class T1, class T2>
440  bool operator<(const Position_<T1>& A, const Position_<T2>& B) {
441  return (A.P[0]<B.P[0]) || (A.P[1]<B.P[1]) || (A.P[2]<B.P[2]);
442 }
444 template<class T1, class T2>
445  bool operator>(const Position_<T1>& A, const Position_<T2>& B) {
446  return (A.P[0]>B.P[0]) || (A.P[1]>B.P[1]) || (A.P[2]>B.P[2]);
447 }
449 template<class T1, class T2>
450 bool operator==(const Position_<T1>& A, const Position_<T2>& B) {
451  return (A.P[0]==B.P[0]) && (A.P[1]==B.P[1]) && (A.P[2]==B.P[2]);
452 }
454 
455 template<class T>
456 std::ostream& operator<<(std::ostream &O, const Position_<T>& A) {return O << A.P[0] << ' ' << A.P[1] << ' ' << A.P[2];}
458 template<class T>
459 std::istream& operator>>(std::istream &I, Position_<T>& A) {return I >> A.P[0] >> A.P[1] >> A.P[2];}
460 
465 template<class T1, class T2>
466 int indexgreater(const Position_<T1>& A, const Position_<T2>& B) {
467  if (A.P[0]>=B.P[0]) return 0;
468  if (A.P[1]>=B.P[1]) return 1;
469  if (A.P[2]>=B.P[2]) return 2;
470 }
471 template<class G>
472 void WeightOnGrid(G& g, const Position& p, double v) {
473  GridPosition gp(p);
474  Position Diff = p - Position(gp);
475  Position diff = Position(1,1,1) - Diff;
476  g[gp] += diff[0]*diff[1]*diff[2]*v;
477  gp[0] += 1; // ++gp.P[0];
478  g[gp] += Diff[0]*diff[1]*diff[2]*v;
479  gp[1] += 1; // ++gp.P[1];
480  g[gp] += Diff[0]*Diff[1]*diff[2]*v;
481  gp[0] -= 1; // --gp.P[0];
482  g[gp] += diff[0]*Diff[1]*diff[2]*v;
483  gp[2] += 1; // ++gp.P[2];
484  g[gp] += diff[0]*Diff[1]*Diff[2]*v;
485  gp[0] += 1; // ++gp.P[0];
486  g[gp] += Diff[0]*Diff[1]*Diff[2]*v;
487  gp[1] -= 1; // --gp.P[1];
488  g[gp] += Diff[0]*diff[1]*Diff[2]*v;
489  gp[0] -= 1; // --gp.P[0];
490  g[gp] += diff[0]*diff[1]*Diff[2]*v;
491 }
492 const GridPosition ZeroGrid(0,0,0);
493 const GridPosition UnitGrid(1,1,1);
494 const GridPosition UnitX(1,0,0);
495 const GridPosition UnitY(0,1,0);
496 const GridPosition UnitZ(0,0,1);
497 const GridPosition UnitnX(-1,0,0);
498 const GridPosition UnitnY(0,-1,0);
499 const GridPosition UnitnZ(0,0,-1);
500 //typedef Velocity_ Position_;
501 #endif
502 //---------------------------------------------------------------------------------------------------------------------
504 template<class T>
505 class Velocity_ {
506  public:
508  T P[3];
510  T& x() {return P[0];};
512  T& y() {return P[1];};
514  T& z() {return P[2];};
516  Velocity_() { }
518  Velocity_(T ax, T ay, T az) {P[0]=ax; P[1]=ay; P[2]=az;}
520  Velocity_(T x) {P[0]=x; P[1]=x; P[2]=x;}
522  operator T*() {return P;}
524  T volume() const {return P[0]*P[1]*P[2];}
526  T sum() const {return P[0]+P[1]+P[2];}
528  double length() const {return sqrt(P[0]*P[0]+P[1]*P[1]+P[2]*P[2]);}
530  void stretch(const Position_<T>& A) {P[0] *= A.P[0]; P[1] *= A.P[1]; P[2] *= A.P[2];}
531 
535  void limitto(const Velocity_<T>& pos) {
536  if (P[0]>pos.P[0]) P[0]=pos.P[0];
537  if (P[1]>pos.P[1]) P[1]=pos.P[1];
538  if (P[2]>pos.P[2]) P[2]=pos.P[2];
539  }
540 
545  T arraypos(const Velocity_<T>& pos) const{return P[0] + (pos.P[0]+1)*(P[1] + (pos.P[1]+1)*P[2]);}
547  template<class T2>
548  operator Velocity_<T2>() const {return Velocity_<T2>(T2(P[0]), T2(P[1]), T2(P[2]));}
550  Velocity_<T> scale(const Velocity_<T> &V) const {return Velocity_<T>(P[0]*V.P[0], P[1]*V.P[1], P[2]*V.P[2]);}
552  double sqr() {return P[0]*P[0]+P[1]*P[1]+P[2]*P[2];}
554 
555  T& operator[](const int i) {return P[i];}
556  T operator[](const int i) const {return P[i];}
558 };
559 
577 
578 template<class T>
579 Velocity_<T> operator-(const Velocity_<T>& A, T b) {return Velocity_<T>(A.P[0]-b, A.P[1]-b, A.P[2]-b);}
581 template<class T>
582 Velocity_<T> operator+(const Velocity_<T>& A, T b) {return Velocity_<T>( A.P[0]+b, A.P[1]+b, A.P[2]+b);}
584 template<class T>
585 Velocity_<T> operator*(const Velocity_<T>& A, T b) {return Velocity_<T>( A.P[0]*b, A.P[1]*b, A.P[2]*b);}
587 template<class T>
588 Velocity_<T> operator/(const Velocity_<T>& A, T b) {return Velocity_<T>( A.P[0]/b, A.P[1]/b, A.P[2]/b);}
590 template<class T>
591 Velocity_<T> operator/(const Velocity_<T>& A, const Velocity_<T>& B) {
592  return Velocity_<T>( A.P[0]/B.P[0], A.P[1]/B.P[1], A.P[2]/B.P[2]);
593 }
595 template<class T>
596 Velocity_<T> operator-(const Velocity_<T>& A, const Velocity_<T>& B) {
597  return Velocity_<T>(A.P[0]-B.P[0], A.P[1]-B.P[1], A.P[2]-B.P[2]);
598 }
600 template<class T>
601 Velocity_<T> operator+(const Velocity_<T>& A, const Velocity_<T>& B) {
602  return Velocity_<T>(A.P[0]+B.P[0], A.P[1]+B.P[1], A.P[2]+B.P[2]);
603 }
605 template<class T1, class T2>
606 double operator*(const Velocity_<T1>& A, const Velocity_<T2>& B) {
607  return A.P[0]*B.P[0] + A.P[1]*B.P[1] + A.P[2]*B.P[2];
608 }
610 template<class T>
611 Velocity_<T>& operator*=(Velocity_<T>& A, const Velocity_<T>& B) {
612  A.P[0] *= B.P[0];
613  A.P[1] *= B.P[1];
614  A.P[2] *= B.P[2];
615  return A;
616 }
618 template<class T1, class T2>
619 Velocity_<T1>& operator*=(Velocity_<T1>& A, const T2& B) {
620  A.P[0] *= B;
621  A.P[1] *= B;
622  A.P[2] *= B;
623  return A;
624 }
626 template<class T>
627 Velocity_<T>& operator/=(Velocity_<T>& A, const T& B) {
628  A.P[0] /= B;
629  A.P[1] /= B;
630  A.P[2] /= B;
631  return A;
632 }
634 template<class T1, class T2>
635 Velocity_<T1>& operator+=(Velocity_<T1>& A, const Velocity_<T2>& B) {
636  A.P[0] += B.P[0];
637  A.P[1] += B.P[1];
638  A.P[2] += B.P[2];
639  return A;
640 }
642 template<class T>
643 Velocity_<T>& operator+=(Velocity_<T>& A, const T& B) {
644  A.P[0] += B;
645  A.P[1] += B;
646  A.P[2] += B;
647  return A;
648 }
650 template<class T>
651 Velocity_<T> operator/(double A, const Velocity_<T>& B) {return Velocity_<T>(A/B.P[0], A/B.P[1], A/B.P[2]);}
653 template<class T1, class T2>
654 bool operator<=(const Velocity_<T1>& A, const Velocity_<T2>& B) {
655  return (A.P[0]<=B.P[0]) || (A.P[1]<=B.P[1]) || (A.P[2]<=B.P[2]);
656 }
658 template<class T1, class T2>
659 bool operator>=(const Velocity_<T1>& A, const Velocity_<T2>& B) {
660  return (A.P[0]>=B.P[0]) || (A.P[1]>=B.P[1]) || (A.P[2]>=B.P[2]);
661 }
663 
668 template<class T1, class T2>
669 int indexgreater(const Velocity_<T1>& A, const Velocity_<T2>& B) {
670  if (A.P[0]>=B.P[0]) return 0;
671  if (A.P[1]>=B.P[1]) return 1;
672  if (A.P[2]>=B.P[2]) return 2;
673 }
675 template<class T>
676 std::ostream& operator<<(std::ostream &O, const Velocity_<T>& A) {return O << A.P[0] << ' ' << A.P[1] << ' ' << A.P[2];}
678 template<class T>
679 std::istream& operator>>(std::istream &I, Velocity_<T>& A) {return I >> A.P[0] >> A.P[1] >> A.P[2];}
680 
685 template<class G, class T>
686 T WeightedSum(const G& g, Velocity_<T> p1, Velocity_<T> p2, const Velocity_<T>& wa, const Velocity_<T>& wb) {
687  T sum = g[p1]*wa.P[0]*wa.P[1]*wa.P[2] + g[p2]*wb.P[0]*wb.P[1]*wb.P[2];
688  ++p1.P[0];
689  --p2.P[0];
690  sum += g[p1]*wb.P[0]*wa.P[1]*wa.P[2] + g[p2]*wa.P[0]*wb.P[1]*wb.P[2];
691  ++p1.P[1];
692  --p2.P[1];
693  sum += g[p1]*wb.P[0]*wb.P[1]*wa.P[2] + g[p2]*wa.P[0]*wa.P[1]*wb.P[2];
694  --p1.P[0];
695  ++p2.P[0];
696  return sum + g[p1]*wa.P[0]*wb.P[1]*wa.P[2] + g[p2]*wb.P[0]*wa.P[1]*wb.P[2];
697 }
698 typedef Velocity_<double> Velocity;
699 const Velocity ZeroVelocity(0,0,0);
700 
704 Position frandPosition();
705 #endif
Generated by   doxygen 1.8.1.1