#include <Multipole.h>

Inheritance diagram for Multipole:

Public Member Functions
	Multipole (int p, NumericalVector<> *squaredFactorials, Vec3D location)

virtual	~Multipole ()

virtual void	computeMultipoleExpansion ()

NumericalVector< std::complex< Mdouble > >	TranslateMultipoleExpansionTo (Vec3D location)

NumericalVector< std::complex< Mdouble > >	convertMultipoleToLocal (Vec3D location)

void	addMultipoleCoefficients (NumericalVector< std::complex< Mdouble >> multipoleExpansionCoefficients)

NumericalVector< std::complex< Mdouble > >	getExpansionCoefficients ()

void	setExpansionCoefficients (NumericalVector< std::complex< Mdouble >> multipoleExpansionCoefficients)

NumericalVector *	getSquaredFactorials ()

int	getP ()

Protected Attributes
int	p_

NumericalVector *	squaredFactorials_

Vec3D	location_

NumericalVector< std::complex< Mdouble > >	multipoleExpansionCoefficients_

Constructor & Destructor Documentation

◆ Multipole()

Multipole::Multipole	(	int	p,
		NumericalVector<> *	squaredFactorials,
		Vec3D	location
	)

                                                                                 :
         p_(p),
         squaredFactorials_(squaredFactorials),
         location_(location)
 {
 }

◆ ~Multipole()

Multipole::~Multipole ( )

virtualdefault

Member Function Documentation

◆ addMultipoleCoefficients()

MERCURYDPM_DEPRECATED void Multipole::addMultipoleCoefficients ( NumericalVector< std::complex< Mdouble >> multipoleExpansionCoefficients )

Adds multipole coefficients to an existing multipole

Deprecated:: Remove this function; it should not be required anymore:

 {
     if (multipoleExpansionCoefficients.size() > multipoleExpansionCoefficients_.size())
     {
         logger(ERROR, "Multipole expansion coefficient sizes are not correct.");
     }
     
     multipoleExpansionCoefficients_ += multipoleExpansionCoefficients;
 }

References ERROR, logger, and multipoleExpansionCoefficients_.

Referenced by Panel::computeMultipoleExpansion(), and Panel::translateMultipoleExpansion().

◆ computeMultipoleExpansion()

void Multipole::computeMultipoleExpansion ( )

virtual

Reimplemented in Dipole.

 {
     int nTerms = 0.5 * (p_ + 1) * (2 * p_ + 2);
     multipoleExpansionCoefficients_(nTerms);
 }

References multipoleExpansionCoefficients_, and p_.

Referenced by Panel::Panel().

◆ convertMultipoleToLocal()

NumericalVector< std::complex< Mdouble > > Multipole::convertMultipoleToLocal ( Vec3D location )

 {
     int nTerms = 0.5 * (p_ + 1) * (2 * p_ + 2);
     NumericalVector<std::complex<Mdouble>> localExpansionCoefficients(nTerms);
     
     //Compute alpha and beta;
     //Todo: use quarternions to do this shite
     Mdouble rho = 1.0;
     Mdouble alpha = 1.0;
     Mdouble beta = 1.0;
     
     NumericalVector<std::complex<Mdouble>> sphericalHarmonics = sphericalHarmonics::sphericalHarmonics(2 * p_, alpha,
                                                                                                        beta);
     
     for (int j = 0; j <= p_; j++)
     {
         for (int k = -j; k <= j; k++)
         {
             std::complex<Mdouble> result = 0.0;
             int location = j * j + (k + j);
             for (int n = 0; n <= p_; n++)
             {
                 for (int m = -n; m <= n; m++)
                 {
                     int location_A1 = n * n + (m + n);
                     int location_A2 = j * j + (j + k);
                     int location_A3 = (j + n) * (j + n) + ((m - k) + (j + n));
                     int location_Y = location_A3;
                     int location_O = location_A1;
                     std::complex<Mdouble> J = std::pow(constants::i, std::abs(k - m) - std::abs(k) - std::abs(m));
                     //\todo TW note: a warning says += cannot be done here
                     result += multipoleExpansionCoefficients_[location_O] * J * (*squaredFactorials_)(location_A1) *
                               (*squaredFactorials_)(location_A2) * sphericalHarmonics[location_Y] /
                               ((*squaredFactorials_)(location_A3) * std::pow(rho, (j + n + 1)));
     
                 }
             }
             localExpansionCoefficients[location] = result;
         }
     }
     return localExpansionCoefficients;
 }

References abs(), alpha, beta, constants::i, J, j, k, m, multipoleExpansionCoefficients_, n, p_, Eigen::bfloat16_impl::pow(), sphericalHarmonics::sphericalHarmonics(), and squaredFactorials_.

◆ getExpansionCoefficients()

NumericalVector<std::complex<Mdouble> > Multipole::getExpansionCoefficients ( )

inline

     {
         return multipoleExpansionCoefficients_;
     }

References multipoleExpansionCoefficients_.

◆ getP()

int Multipole::getP ( )

inline

     {
         return p_;
     }

References p_.

◆ getSquaredFactorials()

NumericalVector* Multipole::getSquaredFactorials ( )

inline

     {
         return squaredFactorials_;
     }

References squaredFactorials_.

◆ setExpansionCoefficients()

void Multipole::setExpansionCoefficients ( NumericalVector< std::complex< Mdouble >> multipoleExpansionCoefficients )

inline

     {
         multipoleExpansionCoefficients_ = multipoleExpansionCoefficients;
     }

References multipoleExpansionCoefficients_.

◆ TranslateMultipoleExpansionTo()

NumericalVector< std::complex< Mdouble > > Multipole::TranslateMultipoleExpansionTo ( Vec3D location )

 {
     //todo: Find a better name for A/squaredFactorials
     int nTerms = 0.5 * (p_ + 1) * (2 * p_ + 2);
     NumericalVector<std::complex<Mdouble>> translatedMultipoleCoefficients(nTerms);
     
     //Check if a multipole expansion has taken place
     if (multipoleExpansionCoefficients_.size() == 0)
     {
         logger(ERROR, "Multipole is not yet expanded.");
     }
     
     //Determine rho, alpha and beta
     //Vec3D distance = location_ - location;
     
     //Todo: Add a quarternion step in here with distance as input, to avoid NaN values in the angles.
     Mdouble rho = 1.0;
     Mdouble alpha = 1.0;
     Mdouble beta = 1.0;
  
 /*  std::cout << "rho=" << rho << std::endl;
     std::cout << "alpha=" << alpha << std::endl;
     std::cout << "beta=" << beta <<std::endl;*/
     
     //Calculate spherical harmonics for alpha and beta
     NumericalVector<std::complex<Mdouble>> sphericalHarmonics = sphericalHarmonics::sphericalHarmonics(p_, alpha, beta);
     
     //Compute the transfered multipole coefficients
     for (int j = 0; j <= p_; j++)
     {
         for (int k = -j; k <= j; k++)
         {
             std::complex<Mdouble> result = 0.0;
             int location = j * j + (k + j);
             for (int n = 0; n <= j; n++)
             {
                 int a = std::max(k + n - j, -n);
                 int b = std::min(k + j - n, n);
                 for (int m = a; m <= b; m++)
                 {
                     int location_O = (j - n) * (j - n) + ((k - m) + (j - n));
                     int location_A1 = n * n + (m + n);
                     int location_Y = n * n + (-m + n);
                     int location_A2 = location_O;
                     int location_A3 = location;
                     result += multipoleExpansionCoefficients_[location_O] *
                               std::pow(constants::i, (std::abs(k) - std::abs(m) - std::abs(k - m)))
                               * (*squaredFactorials_)(location_A1) * (*squaredFactorials_)(location_A2) *
                               std::pow(rho, n) * sphericalHarmonics[location_Y] / (*squaredFactorials_)(location_A3);
                 }
             }
             translatedMultipoleCoefficients[location] = result;
         }
     }
     
     return translatedMultipoleCoefficients;
 }