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