OpenFPM  5.2.0
Project that contain the implementation of distributed structures
SphericalHarmonics.hpp
1 //
2 // Created by Abhinav Singh on 03.11.20.
3 //
4 
5 #ifndef SPHERICALHARMONICS_HPP_
6 #define SPHERICALHARMONICS_HPP_
7 //#include "util/util_debug.hpp"
8 #include <boost/math/special_functions/spherical_harmonic.hpp>
9 
10 //type used for dictionary arguments of a spherical harmonic coordinates
11 typedef std::tuple <int, int> lm;
12 
16 struct key_hash
17 {
18  using argument_type = lm;
19  using result_type = std::size_t;
20 
21  std::size_t operator()(const lm& k) const
22  {
23  return std::get<0>(k) ^ std::get<1>(k);
24  }
25 };
26 
30 struct key_equal
31 {
32  using first_argument_type = lm;
33  using second_argument_type = lm;
34  using result_type = bool;
35 
36  bool operator()(const lm& v0, const lm& v1) const
37  {
38  return (
39  std::get<0>(v0) == std::get<0>(v1) &&
40  std::get<1>(v0) == std::get<1>(v1)
41  );
42  }
43 };
44 
45 
46 
47 
48 namespace openfpm {
49  namespace math {
50 
51  namespace detail {
52  template<class T, class Policy>
53  inline T spherical_harmonic_prefix_raw(unsigned n, unsigned m, T theta, const Policy &pol) {
54  BOOST_MATH_STD_USING
55 
56  if (m > n)
57  return 0;
58  T prefix;
59  if (m==0){
60  prefix=1.0;
61  }
62  else{
63  prefix = boost::math::tgamma_delta_ratio(static_cast<T>(n - m + 1), static_cast<T>(2 * m), pol);
64  }
65  prefix *= (2 * n + 1) / (4 * boost::math::constants::pi<T>());
66  prefix = sqrt(prefix);
67  return prefix;
68  }
69 
70  template<class T, class Policy>
71  T Y(unsigned n, int m, T theta, T phi, const Policy &pol) {
72  BOOST_MATH_STD_USING // ADL of std functions
73 
74  bool sign = false;
75  if (m < 0) {
76  // Reflect and adjust sign if m < 0:
77  sign = m & 1;
78  m = abs(m);
79  }
80  /* if (m & 1) {
81  // Check phase if theta is outside [0, PI]:
82  T mod = boost::math::tools::fmod_workaround(theta, T(2 * constants::pi<T>()));
83  if (mod < 0)
84  mod += 2 * constants::pi<T>();
85  if (mod > constants::pi<T>())
86  sign = !sign;
87  }*/
88  // Get the value and adjust sign as required:
89  T prefix = spherical_harmonic_prefix_raw(n, m, theta, pol);
90  //T sin_theta = sin(theta);
91  T x = cos(theta);
92  T leg = boost::math::legendre_p(n, m, x, pol);
93  if (m != 0)
94  prefix *= sqrt(2);
95  prefix *= leg;
96  return sign ? prefix * sin(m * phi) : prefix * cos(m * phi);
97  }
98 
99  template<class T, class Policy>
100  T DYdTheta(unsigned n, int m, T theta, T phi, const Policy &pol) {
101  BOOST_MATH_STD_USING // ADL of std functions
102 
103  bool sign = false;
104  if (m < 0) {
105  // Reflect and adjust sign if m < 0:
106  sign = m & 1;
107  m = abs(m);
108  }
109 /* if (m & 1) {
110  // Check phase if theta is outside [0, PI]:
111  T mod = boost::math::tools::fmod_workaround(theta, T(2 * constants::pi<T>()));
112  if (mod < 0)
113  mod += 2 * constants::pi<T>();
114  if (mod > constants::pi<T>())
115  sign = !sign;
116  }*/
117  // Get the value and adjust sign as required:
118  T prefix = spherical_harmonic_prefix_raw(n, m, theta, pol);
119  //T sin_theta = sin(theta);
120  T x = cos(theta);
121  T leg2,leg1 = boost::math::legendre_p(n, m + 1, x, pol);
122  if(n+m==0||n-m==-1){
123  leg2=0.0;
124  }
125  else{
126  leg2 = boost::math::legendre_p(n, m - 1, x, pol);
127  }
128  if (m != 0)
129  prefix *= sqrt(2);
130  prefix *= (0.5 * leg1 - 0.5 * (n + m) * (n - m + 1) * leg2);
131  return sign ? prefix * sin(m * phi) : prefix * cos(m * phi);
132  }
133 
134 
135  template<class T, class Policy>
136  T DYdPhi(unsigned n, int m, T theta, T phi, const Policy &pol) {
137  BOOST_MATH_STD_USING // ADL of std functions
138 
139  bool sign = false;
140  if (m == 0)
141  return 0;
142  if (m < 0) {
143  // Reflect and adjust sign if m < 0:
144  sign = m & 1;
145  m = abs(m);
146  }
147  /* if (m & 1) {
148  // Check phase if theta is outside [0, PI]:
149  T mod = boost::math::tools::fmod_workaround(theta, T(2 * constants::pi<T>()));
150  if (mod < 0)
151  mod += 2 * constants::pi<T>();
152  if (mod > constants::pi<T>())
153  sign = !sign;
154  }*/
155  // Get the value and adjust sign as required:
156  T prefix = spherical_harmonic_prefix_raw(n, m, theta, pol);
157  T sin_theta = sin(theta);
158  T x = cos(theta);
159  T leg = boost::math::legendre_p(n, m, x, pol);
160  prefix *= sqrt(2) * leg;
161  return sign ? m * prefix * cos(m * phi) : -m * prefix * sin(m * phi);
162  }
163 
164  }
165 
166  template<class T1, class T2, class Policy>
167  inline typename boost::math::tools::promote_args<T1, T2>::type
168  Y(unsigned n, int m, T1 theta, T2 phi, const Policy &pol) {
169  typedef typename boost::math::tools::promote_args<T1, T2>::type result_type;
170  typedef typename boost::math::policies::evaluation<result_type, Policy>::type value_type;
171  return boost::math::policies::checked_narrowing_cast<result_type, Policy>(
172  detail::Y(n, m, static_cast<value_type>(theta), static_cast<value_type>(phi), pol),
173  "openfpm::math::Y<%1%>(unsigned, int, %1%, %1%)");
174  }
175 
176  template<class T1, class T2>
177  inline typename boost::math::tools::promote_args<T1, T2>::type
178  Y(unsigned n, int m, T1 theta, T2 phi) {
179  return openfpm::math::Y(n, m, theta, phi, boost::math::policies::policy<>());
180  }
181 
182 
183  template<class T1, class T2, class Policy>
184  inline typename boost::math::tools::promote_args<T1, T2>::type
185  DYdTheta(unsigned n, int m, T1 theta, T2 phi, const Policy &pol) {
186  typedef typename boost::math::tools::promote_args<T1, T2>::type result_type;
187  typedef typename boost::math::policies::evaluation<result_type, Policy>::type value_type;
188  return boost::math::policies::checked_narrowing_cast<result_type, Policy>(
189  detail::DYdTheta(n, m, static_cast<value_type>(theta), static_cast<value_type>(phi), pol),
190  "openfpm::math::DYdTheta<%1%>(unsigned, int, %1%, %1%)");
191  }
192 
193  template<class T1, class T2>
194  inline typename boost::math::tools::promote_args<T1, T2>::type
195  DYdTheta(unsigned n, int m, T1 theta, T2 phi) {
196  return openfpm::math::DYdTheta(n, m, theta, phi, boost::math::policies::policy<>());
197  }
198 
199  template<class T1, class T2, class Policy>
200  inline typename boost::math::tools::promote_args<T1, T2>::type
201  DYdPhi(unsigned n, int m, T1 theta, T2 phi, const Policy &pol) {
202  typedef typename boost::math::tools::promote_args<T1, T2>::type result_type;
203  typedef typename boost::math::policies::evaluation<result_type, Policy>::type value_type;
204  return boost::math::policies::checked_narrowing_cast<result_type, Policy>(
205  detail::DYdPhi(n, m, static_cast<value_type>(theta), static_cast<value_type>(phi), pol),
206  "openfpm::math::DYdPhi<%1%>(unsigned, int, %1%, %1%)");
207  }
208 
209  template<class T1, class T2>
210  inline typename boost::math::tools::promote_args<T1, T2>::type
211  DYdPhi(unsigned n, int m, T1 theta, T2 phi) {
212  return openfpm::math::DYdPhi(n, m, theta, phi, boost::math::policies::policy<>());
213  }
214 
215  inline double sph_A1(int l,int m,double v1, double vr) {
216  return 0.5 * (1 + l) * (l * v1 - vr);
217  }
218 
219  inline double sph_A2(int l,int m,double v1, double vr) {
220  return 0.5 * ((1 + l) * (-l) * v1 + (l + 3) * vr);
221  }
222 
223  inline double sph_B(int l, int m,double v2) {
224  return v2;
225  }
226 
227  inline double sph_A3(int l,int m,double v1, double vr) {
228  if (m == 1){
229  return 0.5 *l* ((1 + l)*v1 - vr)-1.5*sph_A2(l,m,v1,vr);
230  }
231  else{
232  return 0.5 *l* ((1 + l)*v1 - vr);
233  }
234  }
235 
236  inline double sph_A4(int l,int m,double v1, double vr) {
237  if (m == 1){
238  return 0.5* (-l*(1 + l)*v1 + (2-l)*vr)+0.5*sph_A2(l,m,v1,vr);
239  }
240  else{
241  return 0.5* (-l*(1 + l)*v1 + (2-l)*vr);
242  }
243  }
244 
257  inline std::vector<double> sph_anasol_u(double nu,int l,int m,double vr,double v1,double v2,double r) {
258  double ur,u1,u2,p;
259  if(l==0)
260  {
261  if(r==0){
262  u2=sph_B(l,m,v2);
263  return {0,0,u2,0.0};
264  }
265  ur=sph_A1(l,m,v1,vr)*r+sph_A2(l,m,v1,vr)/r;
266  u2=sph_B(l,m,v2);
267  return {ur,0,u2,0};
268  }
269  BOOST_MATH_STD_USING
270  //double ur,u1,u2,p;
271  ur=sph_A1(l,m,v1,vr)*pow(r,l+1)+sph_A2(l,m,v1,vr)*pow(r,l-1);//+sph_A3(l,m,v1,vr)*pow(r,-l)+sph_A4(l,m,v1,vr)*pow(r,-l-2);
272  //std::cout<<sph_A1(l,m,v1,vr)<<":"<<sph_A2(l,m,v1,vr)<<":"<<sph_A3(l,m,v1,vr)<<":"<<sph_A4(l,m,v1,vr)<<std::endl;
273  //std::cout<<pow(r,l+1)<<":"<<pow(r,l-1)<<":"<<pow(r,-l)<<":"<<pow(r,-l-2)<<std::endl;
274  u1=1.0/double(l*(l+1))*((l+3.0)*sph_A1(l,m,v1,vr)*pow(r,l+1)+(l+1)*sph_A2(l,m,v1,vr)*pow(r,l-1));//-(l-2)*sph_A3(l,m,v1,vr)*pow(r,-l)-l*sph_A4(l,m,v1,vr)*pow(r,-l-2);
275  u2=sph_B(l,m,v2)*(pow(r,l));//+pow(r,-l-1));
276  p= nu*((4.0*l+6.0)/double(l)*sph_A1(l,m,v1,vr)*pow(r,l));//(4.0*l-2.0)/double(l+1)*sph_A3(l,m,v1,vr)*pow(r,-l-1));
277  return {ur,u1,u2,p};
278  }
279 
291  template<unsigned int k>
292  std::vector<double> sumY(double r, double theta, double phi,const std::unordered_map<const lm,double,key_hash,key_equal> &Vr,const std::unordered_map<const lm,double,key_hash,key_equal> &V1,const std::unordered_map<const lm,double,key_hash,key_equal> &V2) {
293  double Sum1 = 0.0;
294  double Sum2 = 0.0;
295  double Sum3 = 0.0;
296  for (int l = 0; l <= k; l++) {
297  for (int m = -l; m <= l; m++) {
298  auto Er= Vr.find(std::make_tuple(l,m));
299  auto E1= V1.find(std::make_tuple(l,m));
300  auto E2= V2.find(std::make_tuple(l,m));
301  Sum1 += Er->second * openfpm::math::Y(l, m, theta, phi);
302  double DYdPhi=openfpm::math::DYdPhi(l, m, theta, phi);
303  double DYdTheta=openfpm::math::DYdTheta(l, m, theta, phi);
304  /*if (DYdPhi==0 ||E2->second==0){
305  Sum2 += E1->second * DYdTheta;
306  }
307  else{*/
308  //assert(theta!=0);
309  Sum2 += E1->second * DYdTheta -
310  E2->second / sin(theta) * DYdPhi;
311  /*}*/
312  Sum3 += E2->second * DYdTheta +
313  E1->second/sin(theta) * DYdPhi;
314 
315  /* Sum3 += E2->second *sin(theta)* DYdTheta +
316  E1->second * DYdPhi;*/
317 
318  /*if (DYdPhi==0 ||E1->second==0){
319  Sum3 += E2->second * DYdTheta;
320  }
321  else{
322  Sum3 += E2->second * DYdTheta +
323  E1->second/sin(theta) * DYdPhi;
324  }*/
325 
326  }
327  }
328  double x=Sum2*cos(theta)*cos(phi)-Sum3*sin(phi)+Sum1*cos(phi)*sin(theta);
329  double y=Sum3*cos(phi)+Sum2*cos(theta)*sin(phi)+Sum1*sin(phi)*sin(theta);
330  double z=Sum1*cos(theta)-Sum2*sin(theta);
331 
332  return {x,y,z};
333 
334  }
344  template<unsigned int k>
345  double sumY_Scalar(double r, double theta, double phi,const std::unordered_map<const lm,double,key_hash,key_equal> &Vr) {
346  double Sum1 = 0.0;
347  for (int l = 0; l <= k; l++) {
348  for (int m = -l; m <= l; m++) {
349  auto Er= Vr.find(std::make_tuple(l,m));
350  Sum1 += Er->second * openfpm::math::Y(l, m, theta, phi);
351  }
352  }
353  return Sum1;
354 
355  }
356 
357 /* double PsiTheta(unsigned l, int m, double theta, double phi) {
358  return DYdTheta(l, m, theta, phi);
359  }
360 
361  double PsiPhi(unsigned l, int m, double theta, double phi) {
362  return DYdPhi(l, m, theta, phi);
363  }
364  double PhiTheta(unsigned l, int m, double theta, double phi) {
365  return -1/sin(theta)*DYdPhi(l, m, theta, phi);
366  }
367  double PhiPhi(unsigned l, int m, double theta, double phi) {
368  return sin(theta)*DYdTheta(l, m, theta, phi);
369  }*/
370 
371 
372 }
373 }
374 
375 #endif /* SPHERICALHARMONICS_HPP_ */
key_hash
Structure required for the Sph Harmonic amplitude dictionary arguments.
Definition: SphericalHarmonics.hpp:16
openfpm
convert a type into constant type
Definition: aggregate.hpp:301
key_equal
Structure required for the Sph Harmonic amplitude dictionary arguments.
Definition: SphericalHarmonics.hpp:30