MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
LocalExpansion.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 "LocalExpansion.h"
26 #include "Math/ExtendedMath.h"
27 #include "Math/NumericalVector.h"
28 #include "Math/Vector.h"
29 #include "Multipole.h"
30 #include <complex>
31 #include <vector>
32 
33 LocalExpansion::LocalExpansion(int p, NumericalVector<>* squaredFactorials, Vec3D location) :
34  p_(p),
35  squaredFactorials_(squaredFactorials),
36  location_(location)
37 {
38 }
39 
41 = default;
42 
44 {
45  size_t nTerms = 0.5 * (p_ + 1) * (2 * p_ + 2);
47 }
48 
50 {
51  int nTerms = 0.5 * (p_ + 1) * (2 * p_ + 2);
52  NumericalVector<std::complex<Mdouble>> translatedLocalExpansionCoefficients(nTerms);
53  //std::cout << "size: " << nTerms << std::endl;
54 
55  //compute angles and distance in new framework
56  //todo: fix this rubble with quaternions
57  Mdouble rho = 1.0;
58  Mdouble alpha = 1.0;
59  Mdouble beta = 1.0;
60 
61 
62  //Compute translated local expansion coefficients
64 
65  for (int j = 0; j <= p_; j++)
66  {
67  for (int k = -j; k <= j; k++)
68  {
69  std::complex<Mdouble> result = {0.0, 0.0};
70  int location = j * j + (k + j);
71  for (int n = j; n <= p_; n++)
72  {
73  for (int m = (k - j + n); m <= (k - n + j); m++)
74  {
75  int location_O = n * n + (m + n);
76  int location_A1 = (n - j) * (n - j) + ((m - k) + (n - j));
77  int location_A2 = location;
78  int location_Y = location_A1;
79  int location_A3 = location_O;
80  std::complex<Mdouble> J = std::pow(constants::i, std::abs(m) - std::abs(m - k) - std::abs(k));
81  /* std::cout << "location_A1: " << location_O << std::endl;
82  std::cout << "locationExpansion: " << localExpansionCoefficients_[location_O] << std::endl;
83  std::cout << "A1: " << (*squaredFactorials_)(location_A1) << std::endl;
84  std::cout << "A2: " << (*squaredFactorials_)(location_A2) << std::endl;
85  std::cout << "A3: " << (*squaredFactorials_)(location_A3) << std::endl;
86  std::cout << "Spherical: " << sphericalHarmonics[location_Y] << std::endl;*/
87  //\todo TW note: a warning says += cannot be done here
88  result += localExpansionCoefficients_[location_O] * J * (*squaredFactorials_)(location_A1) *
89  (*squaredFactorials_)(location_A2) * sphericalHarmonics[location_Y] *
90  std::pow(rho, n - j) / (std::pow(-1, n + j) * (*squaredFactorials_)(location_A3));
91  }
92  }
93  //std::cout << "location: " << location << std::endl;
94  translatedLocalExpansionCoefficients[location] = result;
95  }
96  }
97  return translatedLocalExpansionCoefficients;
98 }
99 
100 void LocalExpansion::addLocalExpansionCoefficients(NumericalVector<std::complex<Mdouble>> localExpansionCoefficients)
101 {
102  if (localExpansionCoefficients.size() > localExpansionCoefficients_.size())
103  {
104  std::cout << "Multipole expansion coefficient sizes are not correct." << std::endl;
105  std::exit(-1);
106  }
107 
108  localExpansionCoefficients_ += localExpansionCoefficients;
109 
110 }
double Mdouble
Definition: GeneralDefine.h:34
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
void addLocalExpansionCoefficients(NumericalVector< std::complex< Mdouble >> localExpansionCoefficients)
LocalExpansion(int p, NumericalVector<> *squaredFactorials, Vec3D location)
NumericalVector * squaredFactorials_
NumericalVector< std::complex< Mdouble > > localExpansionCoefficients_
NumericalVector< std::complex< Mdouble > > sphericalHarmonics(int p, Mdouble theta, Mdouble phi)
virtual ~LocalExpansion()
void initialiseLocalExpansion()
NumericalVector< std::complex< Mdouble > > translateLocalExpansion(Vec3D location)
Definition: Vector.h:49
std::size_t size() const