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