/*
XLiFE++ is an extended library of finite elements written in C++
Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see .
*/
/*!
\file maxwell3d.cpp
\author D. Martin, E. Lunéville
\since 09 oct 2010
\date 8 sept 2017
\brief Implementation of exact solutions of Maxwell 3d problems
*/
#include "exactSolutions.hpp"
#include "../specialFunctions/specialFunctions.hpp"
namespace xlifepp
{
/*!
Exact solution to the Maxwell 3D problem outside a sphere :
Field scattered by solid sphere of radius Rad of an incoming vector field
Phi_w= (1,0,0)*exp^{(i*K*z}).
with electric field boundary condition Exn=0 on the sphere
Some computed values for K=1
X_1 X_2 X_3 Phi_1 Phi_2 Phi_3
1.0 0.0 0.0 | (1.75532926384 , 1.62107824182 ) (0 , 0 ) (-6.02796306309e-18,-1.58122117268e-18)
0.0 1.0 0.0 | (-1 , 4.81409564708e-17) (0 , 0 ) (0 , 0 )
0.0 0.0 1.0 | (-5.40302305868e-01,-8.41470984808e-01) (0 , 0 ) (0 , 0 )
1.0 1.0 1.0 | (-1.69184878965e-01, 1.61131192686e-01) (2.06589901137e-01, 2.52914553648e-01) (1.19848010392e-01, 3.28542217311e-01)
1.0 1.0 1.2 | (-1.86022424304e-01, 9.18508880491e-02) (1.49911272089e-01, 2.01342400004e-01) (9.37804387189e-02, 3.03355869671e-01)
1.0 1.0 1.4 | (-1.91775499206e-01, 3.34280295554e-02) (1.07142066815e-01, 1.59632591445e-01) (6.64556013889e-02, 2.71335010720e-01)
1.0 1.0 1.6 | (-1.88260942630e-01,-1.51180741005e-02) (7.54265673653e-02, 1.27220194457e-01) (4.11641109484e-02, 2.38176631400e-01)
1.0 1.0 1.8 | (-1.77261848423e-01,-5.51624601093e-02) (5.20175373132e-02, 1.02356587880e-01) (1.91974583103e-02, 2.06737287161e-01)
1.0 1.0 2.0 | (-1.60341338903e-01,-8.78244555191e-02) (3.47085585869e-02, 8.32071602821e-02) (8.24999561215e-04, 1.78106610047e-01)
1.2 1.4 1.0 | (-1.93790114095e-01, 7.12577504506e-02) (1.44884948487e-01, 2.06529528004e-01) (1.39532076671e-02, 1.97067304047e-01)
1.2 1.4 1.2 | (-1.90729686988e-01, 2.88712762800e-02) (1.12819966807e-01, 1.81439122245e-01) (9.10642428927e-03, 1.94073277948e-01)
1.2 1.4 1.4 | (-1.83137187727e-01,-1.04437244804e-02) (8.50160197427e-02, 1.57350273961e-01) (7.19355647585e-04, 1.84905785693e-01)
1.2 1.4 1.6 | (-1.70990951385e-01,-4.55518271208e-02) (6.17400731503e-02, 1.35416184077e-01) (-9.27523672141e-03, 1.71746716026e-01)
1.2 1.4 1.8 | (-1.54657449587e-01,-7.58575264387e-02) (4.27052811485e-02, 1.16028225273e-01) (-1.95184324637e-02, 1.56344795623e-01)
1.2 1.4 2.0 | (-1.34728108413e-01,-1.01074590219e-01) (2.73910156223e-02, 9.91350581457e-02) (-2.91430511034e-02, 1.39939464554e-01)
*/
Vector scatteredFieldMaxwellExn(const Point& p, Parameters& param)
{
// We first recover the current value of the wave number K and radius of sphere
real_t k = 1., rad = 1.;
if (param.contains("k"))
{
k = real((param)("k"));
}
if (param.contains("radius"))
{
rad = real(param("radius"));
}
const real_t kr(k * rad);
const int nmax = 50;
// spherical Bessel functions up to order nmax for argument K*Rad
std::vector js(nmax + 1), ys(nmax + 1);
js = sphericalbesselJ0N(kr, nmax);
ys = sphericalbesselY0N(kr, nmax);
// computation of coefficients an(K*Rad) and bn(K*Rad) until a minimum n is found such that
// | an(K*Rad) | a(nmax), b(nmax);
// H_{0}(K*Rad)
complex_t hn(js[0], ys[0]);
// Use recurrence relation to compute derivative : d/dz[H_{0}](z)=- H_{1}(z)
// d/dz[H_{0}](z)+H_{0}(z)=-H_{1}(z)+H_{0}(z)
complex_t hnpphn = -complex_t(js[1], ys[1]) + hn;
a[0] = -hnpphn.real() / hnpphn;
b[0] = -js[0] / hn;
while (std::abs(a[n]) > theZeroThreshold && n < nmax)
{
n++;
nx2p1 += 2.;
nm1++;
np1++;
hn = complex_t(js[n], ys[n]);
// Use recurrence relation to compute derivative d/dz[H_{n}](z) :
// (2n+1) d/dz[H_{n}](z)=n*H_{n-1}(z)-(n+1)*H_{n+1}(z)
// http://www.math.sfu.ca/~cbm/aands/page_439.htm (formula 10.1.20)
hnpphn = n * complex_t(js[nm1], ys[nm1]) - np1 * complex_t(js[np1], ys[np1]) + nx2p1 * hn;
a[n] = -in * nx2p1 * (hnpphn.real() / hnpphn);
b[n] = -in * nx2p1 * (js[n] / hn);
in *= i1;
}
int nf = n;
if (nf >= nmax)
{
where("scatteredFieldMaxwell_Exn");
error("conv_failed",nmax);
abort();
}
// spherical coordinates rho, theta and phi
real_t r(p[0]*p[0] + p[1]*p[1]), rho(std::sqrt(r + p[2]*p[2]));
r = std::sqrt(r);
real_t cost(p[2] / rho), sint(r / rho), cosp(1.), sinp(0.);
if (r > theEpsilon)
{
cosp = p[0] / r;
sinp = p[1] / r;
}
const dimen_t spaceDim(p.size());
Vector erho(spaceDim), etheta(spaceDim), ephi(spaceDim, 0.);
erho[0] = sint * cosp;
erho[1] = sint * sinp;
erho[2] = cost;
etheta[0] = cost * cosp;
etheta[1] = cost * sinp;
etheta[2] = -sint;
ephi[0] = -sinp;
ephi[1] = cosp;
// Compute partial sum of the series (including Legendre polynomials)
rho *= k;
// spherical Bessel functions up to order nf for argument K*rho
js = sphericalbesselJ0N(rho, nf);
ys = sphericalbesselY0N(rho, nf);
// initialize computation of Legendre polynomials and (2n+1)*cosp
real_t pcnm1 = 1., pcn = cost, pcnp1, n2p1cost = cost, twocost = 2 * cost, psnm1 = 0., psn = sint, psnp1, p0nm1 = 0., p0n = 1., p0np1, d0n = cost;
complex_t crho(0.), ctheta(0.), cphi(0.);
// initialize n-1, n+1, 2n+1, \sum_{1..n} (2i)
nm1 = 0;
np1 = 2;
nx2p1 = 3.;
real_t sumnx2 = 2.;
for (int n = 1; n < nf; nm1++, np1++, nx2p1 += 2, n++)
{
hn = complex_t(js[n], ys[n]);
// Use recurrence relation to compute derivative d/dz[H_{n}](z) :
// (2n+1) d/dz[H_{n}](z)=n*H_{n-1}(z)-(n+1)*H_{n+1}(z)
// to compute hnpphn=d/dz[H_{n}](K*rho)+H_{n}(K*rho)/(K*rho)
hnpphn = (n * complex_t(js[nm1], ys[nm1]) - np1 * complex_t(js[np1], ys[np1])) / nx2p1 + hn / rho;
crho += (a[n] * hn / rho) * psn;
hn *= i1 * b[n];
hnpphn *= a[n];
ctheta += (hnpphn * d0n + hn * p0n) / sumnx2;
cphi -= (hnpphn * p0n + hn * d0n) / sumnx2;
// Recurrence relation for Legendre polynomials
// (n+1) P_{n+1}(cosp)=(2n+1)*cosp*P_{n}(cosp)-n*P_{n-1}(cosp)
n2p1cost += twocost;
sumnx2 += 2 * np1;
pcnp1 = (n2p1cost * pcn - n * pcnm1) / np1;
pcnm1 = pcn;
pcn = pcnp1;
psnp1 = (n2p1cost * psn - np1 * psnm1) / n;
psnm1 = psn;
psn = psnp1;
p0np1 = (n2p1cost * p0n - np1 * p0nm1) / n;
p0nm1 = p0n;
p0n = p0np1;
d0n = sumnx2 * pcn - cost * p0n;
}
return (cosp * (crho * erho + ctheta * etheta) + sinp * cphi * ephi);
}
/*!
This function computes spherical harmonics Y_{l}^{m} up to order n on the unit sphere;
order n is computed by n=ylm.size()-1 and spherical harmonics ylm are computed for l=0,..,n and m=0,..,l.
Y_{l}^{m}=\f$(-1)^{(|m|-m) /2} \sqrt{ (2l+1)/4\pi*(l-|m|)!/(l+|m|)! } P_{l}^{|m|}(\cos(\theta)) e^{im \phi}\f$
which yields for m positive
Y_{l}^{m}=\f$\sqrt{ (2l+1)/4\pi*(l-m)!/(l+m)! } P_{l}^{m}(\cos(\theta)) e^{im \phi}\f$
Argument ylm is a "triangular 2d array" defined as a std::vector of vectors
the 1st std::vector carries the values : Y_{0}^{0}
the 2nd one carries the values : Y_{1}^{0} Y_{1}^{1}
and for l=2,..,n
the l-th std::vector carries the values: Y_{l}^{0} Y_{l}^{1} .. Y_{l}^{l}
if the size of std::vector ylm[l] is l+1.
For a given int n, container ylm can be declared (and defined) prior
to call to this function by
std::vector > ylm(n+1);
for (int l=0; l<=n; l++) ylm[l]=std::vector(l+1);
*/
void sphericalHarmonics(const Point& x, std::vector >& ylm)
{
// spherical coordinates rho, theta and phi
real_t rho(std::sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]));
real_t cost = x[2] / rho; //=cos(theta)
real_t phi = 0.;
if (std::abs(x[0]) > theEpsilon)
{
phi = std::atan(x[1] / x[0]);
}
//
// Compute Legendre functions: P_l^m (m=0,..,l; l=0,..,n)
//
int n = ylm.size() - 1;
// define 2d triangular array container for Legendre functions
std::vector > plm(n + 1);
for (int m = 0; m <= n; m++)
{
plm[m] = std::vector(n + 1 - m);
}
// compute associated Legendre functions
legendreFunctions(cost, plm);
// Compute normalized spherical harmonics: Y_l^m (m=0,..,l; l=0,..,n)
std::vector >::iterator it_yl(ylm.begin());
std::vector::iterator it_l, it_ylm;
int l = 0;
for (it_yl = ylm.begin(); it_yl != ylm.end(); it_yl++, l++)
{
real_t factRatio(1.), over4pi2lp1((2 * l + 1)*over4pi_);
// compute Y_{l}^{0}
*((*it_yl).begin()) = std::sqrt(over4pi2lp1) * plm[0][l];
int m(1);
if (l > 0)
{
factRatio = 1.;
int lpm(l + 1), lpm1(l);
for (it_ylm = (*it_yl).begin() + 1; it_ylm != (*it_yl).end() - 1; m++, lpm++, lpm1--, it_ylm++)
{
// compute ratio of factorials: (l-m)!/(l+m)! for m=1,..., l-1
factRatio /= lpm * lpm1;
// compute Y_{l}^{m} for m=1,..., l-1
*it_ylm = std::sqrt(over4pi2lp1 * factRatio) * plm[m][l - m] * std::exp(complex_t(0, m) * phi);
}
// compute Y_{l}^{l} for l > 0
factRatio /= 2 * l;
*((*it_yl).end() - 1) = std::sqrt(over4pi2lp1 * factRatio) * plm[l][0] * std::exp(complex_t(0, l) * phi);
}
}
}
/*!
This function computes spherical harmonics surface gradients gradY_{l}^{m} up to order n on the unit sphere;
order n is computed by n=ylm.size()-1 and spherical harmonics ylm are
computed for l=0,..,n and m=0,..,l.
Argument gradylm is a "triangular 2d array of 3d-vectors" defined as a vector of vectors of complex_t-Vectors.
the 1st std::vector carries the 3d-vectorial values : gradY_{0}^{0}
the 2nd one carries the 3d-vectorial values : gradY_{1}^{0} gradY_{1}^{1}
and for l=2,..,n
the l-th std::vector carries the 3d-vectorial values: gradY_{l}^{0} gradY_{l}^{1} .. gradY_{l}^{l}
where the size of std::vector gradylm[l] is l+1.
For a given int n, container gradylm can be declared (and defined) prior
to call to this function by
std::vector > > gradylm(n+1);
for (int l=0; l<=n; l++) {
gradylm[l]=std::vector >(l+1);
for (int m=0; m<=l; m++) gradylm[l][m]=Vector(spaceDim);
}
where "spaceDim" is previously defined as the dimension space (necessarily 3 here).
TELL DANIEL : We may put a test that stops the program if spaceDim is not 3.
*/
void sphericalHarmonicsSurfaceGrad(const Point& x, std::vector > >& gradylm)
{
// spherical coordinates rho, theta and phi
real_t rho = std::sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
real_t cost = x[2] / rho; //=cos(theta)
real_t sint = std::sqrt(1. - cost * cost);
real_t phi(0.), cosp(1.), sinp(0.);
if (std::abs(x[0]) > theEpsilon)
{
phi = std::atan(x[1] / x[0]);
cosp = std::cos(phi);
sinp = std::sin(phi);
}
const dimen_t spaceDim(x.size());
Vector etheta(spaceDim), ephi(spaceDim, 0.);
etheta[0] = cost * cosp;
etheta[1] = cost * sinp;
etheta[2] = -sint;
ephi[0] = -sinp;
ephi[1] = cosp;
// Compute Legendre functions: P_l^m (m=0,..,l; l=0,..,n)
int n = gradylm.size() - 1;
// define 2d triangular array container for Legendre functions
std::vector > plm(n + 1);
for (int m = 0; m <= n; m++)
{
plm[m] = std::vector(n + 1 - m);
}
// compute associated Legendre functions
legendreFunctions(cost, plm);
// Compute derivative of Legendre functions: P_l^m (m=0,..,l; l=0,..,n)
real_t factRatio(1.), over4pi2lp1;
complex_t coeffPhi;// coeff_theta, coeff
// define 2d triangular array container for Legendre functions derivatives
std::vector > dplm(n + 1);
for (int l = 0; l <= n; l++)
{
dplm[l] = std::vector(l + 1);
}
// compute derivative of associated Legendre functions
legendreFunctionsDerivative(cost, plm, dplm);
std::vector >::iterator it_dpl(dplm.begin());
std::vector::iterator it_dplm;
// Compute normalized spherical harmonics derivatives: gradY_l^m (m=0,..,l; l=0,..,n)
std::vector > >::iterator it_gradyl(gradylm.begin());
std::vector >::iterator it_gradylm;
if (std::abs(sint) < theEpsilon)
{
*((*it_gradyl).begin()) = complex_t(0.); // Case l=0 where gradylm is equal to 0.
int l = 1;
real_t costl1 = cost;
it_dpl = dplm.begin() + 1;
for (it_gradyl = gradylm.begin() + 1; it_gradyl != gradylm.end(); it_gradyl++, l++, it_dpl++)
{
int llp1 = l * (l + 1);
factRatio = 1. / llp1;
over4pi2lp1 = (2 * l + 1) * over4pi_;
*((*it_gradyl).begin()) = complex_t(0.);
costl1 *= cost;
// compute gradY_{l}^{1} -- the other are 0
// gradY_{l}^{1}=\sqrt{ (2l+1)/4*pi*(l-1)!/(l+1)! }*e^{i phi} *
// (P'_{l}^{1}(\cos(theta)) e_theta-i (cos(theta))^{l+1} cos(phi) l(l+1)/2*e_phi)
// where e_theta and e_phi are the unit tangent std::vector on the unit sphere.
//
it_gradylm = (*it_gradyl).begin() + 1;
it_dplm = (*it_dpl).begin() + 1;
coeffPhi = - complex_t(0, 0.5 * costl1 * cosp * llp1);
*it_gradylm = std::sqrt(over4pi2lp1 * factRatio) * std::exp(complex_t(0, 1) * phi) * (*it_dplm * etheta + coeffPhi * ephi);
for (it_gradylm = (*it_gradyl).begin() + 2; it_gradylm != (*it_gradyl).end(); it_gradylm++)
{
*it_gradylm = complex_t(0.); // Case l>1 where gradylm is 0.
}
}
}
else
{
*((*it_gradyl).begin()) = complex_t(0.); // Case l=0 where gradylm is 0.
int l = 1;
it_dpl = dplm.begin() + 1;
for (it_gradyl = gradylm.begin() + 1; it_gradyl != gradylm.end(); it_gradyl++, l++, it_dpl++)
{
int lpm = l + 1, lpm1 = l, m = 0;
factRatio = 1.;
over4pi2lp1 = (2 * l + 1) * over4pi_;
it_dplm = (*it_dpl).begin();
for (it_gradylm = (*it_gradyl).begin(); it_gradylm != (*it_gradyl).end() - 1; m++, lpm++, lpm1--, it_gradylm++, it_dplm++)
{
// compute gradY_{l}^{m} for m=1,..., l
// gradY_{l}^{m}=sqrt{ (2l+1)/4*pi*(l-m)!/(l+m)! }*e^{i m phi} *
// (P'_{l}^{|m|}(cos(theta)) e_theta+im*P_{l}^{|m|}(cos(theta)) e_phi)
// where e_theta and e_phi are the unit tangent std::vector on the unit sphere.
//
// coeff=sqrt(over4pi2lp1*factRatio)*exp(complex_t(0,m)*phi);
// coeff_theta=*it_dplm;
coeffPhi = complex_t(0, m) * plm[m][l - m] / sint;
*it_gradylm = std::sqrt(over4pi2lp1 * factRatio) * std::exp(complex_t(0, m) * phi) * (*it_dplm * etheta + coeffPhi * ephi);
// compute ratio of factorials: (l-(m+1))!/(l+(m+1))! for m=O,..., l-1
factRatio /= lpm * lpm1;
}
it_gradylm = (*it_gradyl).end() - 1;
it_dplm = (*it_dpl).end() - 1;
coeffPhi = complex_t(0, l * plm[l][0] / sint);
*it_gradylm = std::sqrt(over4pi2lp1 * factRatio) * std::exp(complex_t(0, l) * phi) * (*it_dplm * etheta + coeffPhi * ephi);
}
}
}
// FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST //
void sphericalHarmonicsTest(const Point& x, const number_t n, std::ostream& out)
{
// Output function for Spherical Harmonics
std::vector > yml(n + 1);
for (number_t l = 0; l <= yml.size() - 1; l++)
{
yml[l] = std::vector(l + 1);
}
sphericalHarmonics(x, yml);
out << " Spherical harmonics Y_l^m up to order " << yml.size() - 1 << " (m=0,..,l; l=0,..," << yml.size() - 1 << ")";
out << " for point(" << x[0] << ", " << x[1] << ", " << x[2] << ")" << std::endl;
out.setf(std::ios::scientific);
int l = 0;
for (std::vector >::iterator it_yl = yml.begin(); it_yl != yml.end(); it_yl++, l++)
{
int m = 0;
out << " l=" << l << std::endl;
for (std::vector::iterator it_yml = (*it_yl).begin(); it_yml != (*it_yl).end(); it_yml++, m++)
{
out << " Y_" << l << "^" << m << "=" << std::setw(19) << std::setprecision(12) << std::endl;
}
}
out.unsetf(std::ios::scientific);
}
void sphericalHarmonicsSurfaceGradTest(const Point& x, const number_t n, std::ostream& out)
{
// Output function for Spherical Harmonics surface gradients
const dimen_t spaceDim(x.size());
std::vector > > gradyml(n + 1);
for (number_t l = 0; l <= gradyml.size() - 1; l++)
{
gradyml[l] = std::vector >(l + 1);
for (number_t m = 0; m <= gradyml[l].size() - 1; m++)
{
gradyml[l][m] = Vector(spaceDim);
}
}
sphericalHarmonicsSurfaceGrad(x, gradyml);
out << " Spherical harmonics derivatives gradY_l^m up to order " << gradyml.size() - 1 << " (m=0,..,l; l=0,..," << gradyml.size() - 1 << ")";
out << " for point(" << x[0] << ", " << x[1] << ", " << x[2] << ")" << std::endl;
out.setf(std::ios::scientific);
int l = 0;
for (std::vector > >::iterator it_gradyl = gradyml.begin(); it_gradyl != gradyml.end(); it_gradyl++, l++)
{
int m = 0;
out << " l=" << l << std::endl;
for (std::vector >::iterator grad_it_yml = (*it_gradyl).begin(); grad_it_yml != (*it_gradyl).end(); grad_it_yml++, m++)
{
out << " gradY_" << l << "^" << m << "=("
<< std::setw(19) << std::setprecision(12) << (*grad_it_yml)[0] << " , " << (*grad_it_yml)[1] << " , " << (*grad_it_yml)[2] << "). " << std::endl;
}
}
out.unsetf(std::ios::scientific);
}
} // end of namespace xlifepp