SlowDown Algorithmic Regularization (SDAR)
Algorithmic Regularization with slowdown method for integrating few-body motions
matrix.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include<cmath>
4 #include<iomanip>
5 
6 namespace COMM{
7  template<class T>
8  class Vector3{
9  public:
10  T x,y,z;
11 
12  Vector3() {}
13  Vector3(const T _x, const T _y, const T _z): x(_x), y(_y), z(_z) {}
14 
15  T operator * (const Vector3<T> & _v) const {
16  return x*_v.x + y*_v.y + z*_v.z;
17  }
18 
19  Vector3 operator ^ (const Vector3 & rhs) const{
20  return Vector3( (y * rhs.z - z * rhs.y),
21  (z * rhs.x - x * rhs.z),
22  (x * rhs.y - y * rhs.x) );
23  }
24 
25  };
26 
27 
28  template<class T>
29  class Matrix3{
30  public:
31  T xx, xy, xz, yx, yy, yz, zx, zy, zz;
32  Matrix3() : xx(T(0)), xy(T(0)), xz(T(0)), yx(T(0)), yy(T(0)), yz(T(0)), zx(T(0)), zy(T(0)), zz(T(0)) {}
33 
34  Matrix3(const T c) : xx(c), xy(c), xz(c), yx(c), yy(c), yz(c), zx(c), zy(c), zz(c) {}
35 
36  Matrix3(const Matrix3 & c) : xx(c.xx), xy(c.xy), xz(c.xz), yx(c.yx), yy(c.yy), yz(c.yz), zx(c.zx), zy(c.zy), zz(c.zz) {}
37 
38 
39  Matrix3(const T _xx, const T _xy, const T _xz,
40  const T _yx, const T _yy, const T _yz,
41  const T _zx, const T _zy, const T _zz)
42  : xx(_xx), xy(_xy), xz(_xz), yx(_yx), yy(_yy), yz(_yz), zx(_zx), zy(_zy), zz(_zz) {}
43 
44  template <typename U>
45  operator Matrix3<U> () const {
46  return Matrix3<U>( static_cast<U>(xx), static_cast<U>(xy), static_cast<U>(xz),
47  static_cast<U>(yx), static_cast<U>(yy), static_cast<U>(yz),
48  static_cast<U>(zx), static_cast<U>(zy), static_cast<U>(zz) );
49  }
50 
51  void rotation(const T I, const T OMEGA, const T omega){
52  const T cosomega = cos(omega);
53  const T sinomega = sin(omega);
54  const T cosOMEGA = cos(OMEGA);
55  const T sinOMEGA = sin(OMEGA);
56  const T cosinc = cos(I);
57  const T sininc = sin(I);
58  xx = cosomega*cosOMEGA - sinomega*sinOMEGA*cosinc;
59  xy = -sinomega*cosOMEGA - cosomega*sinOMEGA*cosinc;
60  xz = sinOMEGA*sininc;
61 
62  yx = cosomega*sinOMEGA + sinomega*cosOMEGA*cosinc;
63  yy = -sinomega*sinOMEGA + cosomega*cosOMEGA*cosinc;
64  yz = -cosOMEGA*sininc;
65 
66  zx = sinomega*sininc;
67  zy = cosomega*sininc;
68  zz = cosinc;
69  }
70 
71 
72  template<class Tvec>
73  Tvec operator * (const Tvec & vec) const {
74  Tvec ret;
75  ret.x = xx*vec.x + xy*vec.y + xz*vec.z;
76  ret.y = yx*vec.x + yy*vec.y + yz*vec.z;
77  ret.z = zx*vec.x + zy*vec.y + zz*vec.z;
78  return ret;
79  }
80  /*
81  template<class Tvec>
82  friend Tvec operator * (const Tvec & vec, const Matrix3 & mat){
83  Tvec ret;
84  ret.x = mat.xx*vec.x + mat.xy*vec.y + mat.xz*vec.z;
85  ret.y = mat.yx*vec.x + mat.yy*vec.y + mat.yz*vec.z;
86  ret.z = mat.zx*vec.x + mat.zy*vec.y + mat.zz*vec.z;
87  return ret;
88  }
89  */
90 
91  friend std::ostream & operator << (std::ostream & c, Matrix3 & mtmp){
92  c<<std::setprecision(15)<<mtmp.xx<<" "<<mtmp.xy<<" "<<mtmp.xz<<std::endl;
93  c<<std::setprecision(15)<<mtmp.yx<<" "<<mtmp.yy<<" "<<mtmp.yz<<std::endl;
94  c<<std::setprecision(15)<<mtmp.zx<<" "<<mtmp.zy<<" "<<mtmp.zz<<std::endl;
95  return c;
96  }
97 
98  /*
99  const T * operator[](const int & i)const{return m[i];}
100 
101  T * operator[](const int & i){return m[i];}
102 
103 
104 
105  Matrix3& operator=(const Matrix3 &c)const{
106  m[0][0] = c.m[0][0]; m[0][1] = c.m[0][1]; m[0][2] = c.m[0][2];
107  m[1][0] = c.m[1][0]; m[1][1] = c.m[1][1]; m[1][2] = c.m[1][2];
108  m[2][0] = c.m[2][0]; m[2][1] = c.m[2][1]; m[2][2] = c.m[2][2];
109  return *this;
110  }
111 
112  Matrix3 & operator+(const Matrix3 &a)const{
113  Matrix3 x;
114  x.m[0][0]=m[0][0]+a.m[0][0]; x.m[0][1]=m[0][1]+a.m[0][1]; x.m[0][2]=m[0][2]+a.m[0][2];
115  x.m[1][0]=m[1][0]+a.m[1][0]; x.m[1][1]=m[1][1]+a.m[1][1]; x.m[1][2]=m[1][2]+a.m[1][2];
116  x.m[2][0]=m[2][0]+a.m[2][0]; x.m[2][1]=m[2][1]+a.m[2][1]; x.m[2][2]=m[2][2]+a.m[2][2];
117  return x;
118  }
119 
120  matrix3& operator-(const matrix3 &a)const{
121  matrix3 x;
122  x.m[0][0]=m[0][0]-a.m[0][0]; x.m[0][1]=m[0][1]-a.m[0][1]; x.m[0][2]=m[0][2]-a.m[0][2];
123  x.m[1][0]=m[1][0]-a.m[1][0]; x.m[1][1]=m[1][1]-a.m[1][1]; x.m[1][2]=m[1][2]-a.m[1][2];
124  x.m[2][0]=m[2][0]-a.m[2][0]; x.m[2][1]=m[2][1]-a.m[2][1]; x.m[2][2]=m[2][2]-a.m[2][2];
125  return x;
126  }
127 
128 
129  matrix3& operator*(const matrix3 &a)const{
130  matrix3 x;
131  x.m[0][0] = m[0][0]*a.m[0][0] + m[0][1]*a.m[1][0] + m[0][2]*a.m[2][0];
132  x.m[0][1] = m[0][0]*a.m[0][1] + m[0][1]*a.m[1][1] + m[0][2]*a.m[2][1];
133  x.m[0][2] = m[0][0]*a.m[0][2] + m[0][1]*a.m[1][2] + m[0][2]*a.m[2][2];
134 
135  x.m[1][0] = m[1][0]*a.m[0][0] + m[1][1]*a.m[1][0] + m[1][2]*a.m[2][0];
136  x.m[1][1] = m[1][0]*a.m[0][1] + m[1][1]*a.m[1][1] + m[1][2]*a.m[2][1];
137  x.m[1][2] = m[1][0]*a.m[0][2] + m[1][1]*a.m[1][2] + m[1][2]*a.m[2][2];
138 
139  x.m[2][0] = m[2][0]*a.m[0][0] + m[2][1]*a.m[1][0] + m[2][2]*a.m[2][0];
140  x.m[2][1] = m[2][0]*a.m[0][1] + m[2][1]*a.m[1][1] + m[2][2]*a.m[2][1];
141  x.m[2][2] = m[2][0]*a.m[0][2] + m[2][1]*a.m[1][2] + m[2][2]*a.m[2][2];
142  return x;
143  }
144 
145 
146  vector3<T>& operator * (const vector3<T> &a)const{
147  vector3<T> x;
148  x[0] = m[0][0]*a[0] + m[0][1]*a[1] + m[0][2]*a[2];
149  x[1] = m[1][0]*a[0] + m[1][1]*a[1] + m[1][2]*a[2];
150  x[2] = m[2][0]*a[0] + m[2][1]*a[1] + m[2][2]*a[2];
151  return x;
152  }
153 
154  matrix3& operator * (const T &s) const{
155  matrix3 x;
156  x.m[0][0]=m[0][0]*s; x.m[0][1]=m[0][1]*s; x.m[0][2]=m[0][2]*s;
157  x.m[1][0]=m[1][0]*s; x.m[1][1]=m[1][1]*s; x.m[1][2]=m[1][2]*s;
158  x.m[2][0]=m[2][0]*s; x.m[2][1]=m[2][1]*s; x.m[2][2]=m[2][2]*s;
159  return x;
160  }
161 
162  matrix3& transposed()const{
163  matrix3 mtmp;
164  for(int i=0; i<3; i++){
165  for(int j=0; j<3; j++){
166  mtmp.m[j][i] = m[i][j];
167  }
168  }
169  return mtmp;
170  }
171 
172  void unit(){
173  m[0][1] = m[0][2] = m[1][0] = m[1][2] = m[2][0] = m[2][1] = 0.0;
174  m[0][0] = m[1][1] = m[2][2] = 1.0;
175  }
176 
177  matrix3& operator += (const matrix3 &a){
178  m[0][0] += a.m[0][0]; m[0][1] += a.m[0][1]; m[0][2] += a.m[0][2];
179  m[1][0] += a.m[1][0]; m[1][1] += a.m[1][1]; m[1][2] += a.m[1][2];
180  m[2][0] += a.m[2][0]; m[2][1] += a.m[2][1]; m[2][2] +=a.m[2][2];
181  return *this;
182  }
183 
184  matrix3& operator -= (const matrix3 &a){
185  m[0][0] -= a.m[0][0]; m[0][1] -= a.m[0][1]; m[0][2] -= a.m[0][2];
186  m[1][0] -= a.m[1][0]; m[1][1] -= a.m[1][1]; m[1][2] -= a.m[1][2];
187  m[2][0] -= a.m[2][0]; m[2][1] -= a.m[2][1]; m[2][2] -=a.m[2][2];
188  return *this;
189  }
190 
191  matrix3& operator *= (const T &a){
192  m[0][0] *= a; m[0][1] *= a; m[0][2] *= a;
193  m[1][0] *= a; m[1][1] *= a; m[1][2] *= a;
194  m[2][0] *= a; m[2][1] *= a; m[2][2] *= a;
195  return *this;
196  }
197 
198  T& determinant()const{
199  T det = 0.0;
200  det = m[0][0] * m[1][1] * m[2][2]
201  + m[0][1] * m[1][2] * m[2][0]
202  + m[0][2] * m[1][0] * m[2][1]
203  - m[0][2] * m[1][1] * m[2][0]
204  - m[0][1] * m[1][0] * m[2][2]
205  - m[0][0] * m[1][2] * m[2][1];
206  return det;
207  }
208 
209  friend matrix3 outer_product(const vector3<T> &a, const vector3<T> &b){
210  matrix3 x;
211  x[0][0]=a[0]*b[0]; x[0][1]=a[0]*b[1]; x[0][2]=a[0]*b[2];
212  x[1][0]=a[1]*b[0]; x[1][1]=a[1]*b[1]; x[1][2]=a[1]*b[2];
213  x[2][0]=a[2]*b[0]; x[2][1]=a[2]*b[1]; x[2][2]=a[2]*b[2];
214  return x;
215  }
216 
217  friend matrix3 quadrupole(const vector3<T> &a){
218  matrix3 x=3.0*outer_product(a,a);
219  T a_squared=a*a;
220  matrix3 delta;
221  delta.unit();
222  return x-a_squared*delta;
223  }
224 
225 friend matrix3 inertiamoment(const T mass, const vector3<T> &pos){
226 matrix3 x = mass * outer_product(pos, pos);
227 return x;
228 }
229 
230 friend vector3<T> operator*(const vector3<T> &a, const matrix3 &c){
231 vector3<T> x;
232 x[0]=a[0]*c.m[0][0] + a[1]*c.m[1][0] + a[2]*c.m[2][0];
233 x[1]=a[0]*c.m[0][1] + a[1]*c.m[1][1] + a[2]*c.m[2][1];
234 x[2]=a[0]*c.m[0][2] + a[1]*c.m[1][2] + a[2]*c.m[2][2];
235 return x;
236 }
237 
238 friend matrix3 operator * (const T &s, const matrix3 &c){
239 matrix3 x;
240 x[0][0]=c.m[0][0]*s; x[0][1]=c.m[0][1]*s; x[0][2]=c.m[0][2]*s;
241 x[1][0]=c.m[1][0]*s; x[1][1]=c.m[1][1]*s; x[1][2]=c.m[1][2]*s;
242 x[2][0]=c.m[2][0]*s; x[2][1]=c.m[2][1]*s; x[2][2]=c.m[2][2]*s;
243 return x;
244 }
245  */
246 
247 
248  };
249 
250 }
COMM
Definition: binary_tree.h:8
COMM::Matrix3::zx
T zx
Definition: matrix.h:31
COMM::Vector3
Definition: matrix.h:8
COMM::Matrix3
Definition: matrix.h:29
COMM::Vector3::Vector3
Vector3()
Definition: matrix.h:12
COMM::Matrix3::xz
T xz
Definition: matrix.h:31
COMM::Matrix3::zy
T zy
Definition: matrix.h:31
COMM::Matrix3::rotation
void rotation(const T I, const T OMEGA, const T omega)
Definition: matrix.h:51
COMM::Matrix3::xy
T xy
Definition: matrix.h:31
COMM::Matrix3::yz
T yz
Definition: matrix.h:31
COMM::Matrix3::Matrix3
Matrix3(const Matrix3 &c)
Definition: matrix.h:36
COMM::Matrix3::Matrix3
Matrix3(const T _xx, const T _xy, const T _xz, const T _yx, const T _yy, const T _yz, const T _zx, const T _zy, const T _zz)
Definition: matrix.h:39
COMM::Matrix3::operator*
Tvec operator*(const Tvec &vec) const
Definition: matrix.h:73
COMM::Matrix3::yx
T yx
Definition: matrix.h:31
COMM::Vector3::y
T y
Definition: matrix.h:10
COMM::Matrix3::zz
T zz
Definition: matrix.h:31
COMM::Matrix3::xx
T xx
Definition: matrix.h:31
COMM::Vector3::operator^
Vector3 operator^(const Vector3 &rhs) const
Definition: matrix.h:19
COMM::Vector3::z
T z
Definition: matrix.h:10
COMM::Vector3::Vector3
Vector3(const T _x, const T _y, const T _z)
Definition: matrix.h:13
COMM::Matrix3::yy
T yy
Definition: matrix.h:31
COMM::Matrix3::Matrix3
Matrix3(const T c)
Definition: matrix.h:34
COMM::Vector3::x
T x
Definition: matrix.h:10
COMM::Matrix3::Matrix3
Matrix3()
Definition: matrix.h:32
COMM::Vector3::operator*
T operator*(const Vector3< T > &_v) const
Definition: matrix.h:15
COMM::Matrix3::operator<<
friend std::ostream & operator<<(std::ostream &c, Matrix3 &mtmp)
Definition: matrix.h:91