MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Multipole.cc
Go to the documentation of this file.
1 //Copyright (c) 2013-2020, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 #include "Math/ExtendedMath.h"
26 #include "Math/NumericalVector.h"
27 #include "Math/Vector.h"
28 #include "Multipole.h"
29 #include <cmath>
30 #include <iostream>
31 #include <vector>
32 
33 Multipole::Multipole(int p, NumericalVector<>* squaredFactorials, Vec3D location) :
34  p_(p),
35  squaredFactorials_(squaredFactorials),
36  location_(location)
37 {
38 }
39 
41 = default;
42 
44 {
45  int nTerms = 0.5 * (p_ + 1) * (2 * p_ + 2);
47 }
48 
50 {
51  //todo: Find a better name for A/squaredFactorials
52  int nTerms = 0.5 * (p_ + 1) * (2 * p_ + 2);
53  NumericalVector<std::complex<Mdouble>> translatedMultipoleCoefficients(nTerms);
54 
55  //Check if a multipole expansion has taken place
57  {
58  std::cout << "Multipole is not yet expanded." << std::endl;
59  std::exit(-1);
60  }
61 
62  //Determine rho, alpha and beta
63  Vec3D distance = location_ - location;
64 
65  //Todo: Add a quarternion step in here with distance as input, to avoid NaN values in the angles.
66  Mdouble rho = 1.0;
67  Mdouble alpha = 1.0;
68  Mdouble beta = 1.0;
69 
70 /* std::cout << "rho=" << rho << std::endl;
71  std::cout << "alpha=" << alpha << std::endl;
72  std::cout << "beta=" << beta <<std::endl;*/
73 
74  //Calculate spherical harmonics for alpha and beta
76 
77  //Compute the transfered multipole coefficients
78  for (int j = 0; j <= p_; j++)
79  {
80  for (int k = -j; k <= j; k++)
81  {
82  std::complex<Mdouble> result = 0.0;
83  int location = j * j + (k + j);
84  for (int n = 0; n <= j; n++)
85  {
86  int a = std::max(k + n - j, -n);
87  int b = std::min(k + j - n, n);
88  for (int m = a; m <= b; m++)
89  {
90  int location_O = (j - n) * (j - n) + ((k - m) + (j - n));
91  int location_A1 = n * n + (m + n);
92  int location_Y = n * n + (-m + n);
93  int location_A2 = location_O;
94  int location_A3 = location;
95  result += multipoleExpansionCoefficients_[location_O] *
96  std::pow(constants::i, (std::abs(k) - std::abs(m) - std::abs(k - m)))
97  * (*squaredFactorials_)(location_A1) * (*squaredFactorials_)(location_A2) *
98  std::pow(rho, n) * sphericalHarmonics[location_Y] / (*squaredFactorials_)(location_A3);
99  }
100  }
101  translatedMultipoleCoefficients[location] = result;
102  }
103  }
104 
105  return translatedMultipoleCoefficients;
106 }
107 
109 {
110  int nTerms = 0.5 * (p_ + 1) * (2 * p_ + 2);
111  NumericalVector<std::complex<Mdouble>> localExpansionCoefficients(nTerms);
112 
113  //Compute alpha and beta;
114  //Todo: use quarternions to do this shite
115  Mdouble rho = 1.0;
116  Mdouble alpha = 1.0;
117  Mdouble beta = 1.0;
118 
120  beta);
121 
122  for (int j = 0; j <= p_; j++)
123  {
124  for (int k = -j; k <= j; k++)
125  {
126  std::complex<Mdouble> result = 0.0;
127  int location = j * j + (k + j);
128  for (int n = 0; n <= p_; n++)
129  {
130  for (int m = -n; m <= n; m++)
131  {
132  int location_A1 = n * n + (m + n);
133  int location_A2 = j * j + (j + k);
134  int location_A3 = (j + n) * (j + n) + ((m - k) + (j + n));
135  int location_Y = location_A3;
136  int location_O = location_A1;
137  std::complex<Mdouble> J = std::pow(constants::i, std::abs(k - m) - std::abs(k) - std::abs(m));
138  //\todo TW note: a warning says += cannot be done here
139  result += multipoleExpansionCoefficients_[location_O] * J * (*squaredFactorials_)(location_A1) *
140  (*squaredFactorials_)(location_A2) * sphericalHarmonics[location_Y] /
141  ((*squaredFactorials_)(location_A3) * std::pow(rho, (j + n + 1)));
142 
143  }
144  }
145  localExpansionCoefficients[location] = result;
146  }
147  }
148  return localExpansionCoefficients;
149 }
150 
153 void Multipole::addMultipoleCoefficients(NumericalVector<std::complex<Mdouble>> multipoleExpansionCoefficients)
154 {
155  if (multipoleExpansionCoefficients.size() > multipoleExpansionCoefficients_.size())
156  {
157  std::cout << "Multipole expansion coefficient sizes are not correct." << std::endl;
158  std::exit(-1);
159  }
160 
161  multipoleExpansionCoefficients_ += multipoleExpansionCoefficients;
162 }
virtual void computeMultipoleExpansion()
Definition: Multipole.cc:43
NumericalVector * squaredFactorials_
Definition: Multipole.h:72
Vec3D location_
Definition: Multipole.h:73
double Mdouble
Definition: GeneralDefine.h:34
void addMultipoleCoefficients(NumericalVector< std::complex< Mdouble >> multipoleExpansionCoefficients)
Adds multipole coefficients to an existing multipole todo: remove this function; it should not be req...
Definition: Multipole.cc:153
Mdouble beta(Mdouble z, Mdouble w)
This is the beta function, returns the approximation based on cmath's implementation of ln(gamma) ...
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
NumericalVector< std::complex< Mdouble > > convertMultipoleToLocal(Vec3D location)
Definition: Multipole.cc:108
int p_
Definition: Multipole.h:71
NumericalVector< std::complex< Mdouble > > TranslateMultipoleExpansionTo(Vec3D location)
Definition: Multipole.cc:49
virtual ~Multipole()
NumericalVector< std::complex< Mdouble > > multipoleExpansionCoefficients_
Definition: Multipole.h:74
NumericalVector< std::complex< Mdouble > > sphericalHarmonics(int p, Mdouble theta, Mdouble phi)
Multipole(int p, NumericalVector<> *squaredFactorials, Vec3D location)
Definition: Multipole.cc:33
Definition: Vector.h:49
std::size_t size() const