1 //
2 // SHCalc.cpp
3 // pb_solvers_code
4 //
5 // Created by David Brookes on 9/25/15.
6 // Copyright © 2015 David Brookes. All rights reserved.
7 //
8
9 #include "SHCalc.h"
10
SHCalcConstants(const int N)11 SHCalcConstants::SHCalcConstants(const int N)
12 :numVals_(N), legConsts1_(N, N), legConsts2_(N, N),
13 shConsts_(N, N), dubFac_(N)
14 {
15 vector<double> temp;
16 temp.reserve(2 * numVals_);
17 temp.push_back(1.0);
18 int i, n, m;
19 for (i = 1; i < 2 * numVals_; i++)
20 {
21 temp.push_back(temp[i-1] * sqrt(i));
22 }
23
24 for (n = 0; n < numVals_; n++)
25 {
26 for (m = 0; m <= n; m++)
27 {
28 if ((n-m) != 0 )
29 {
30 legConsts1_.set_val(n, m, (2*n-1) / (double) (n-m));
31 legConsts2_.set_val(n, m, (n+m-1) / (double) (n-m));
32 }
33 else
34 {
35 legConsts1_.set_val(n, m, 0.0);
36 legConsts2_.set_val(n, m, 0.0);
37 }
38 shConsts_.set_val(n, m, temp[n-m] / temp[n+m]);
39 }
40 }
41
42 dubFac_[0] = 1.0;
43 dubFac_[1] = 1.0;
44 for (i = 2; i < numVals_; i++)
45 {
46 dubFac_[i] = dubFac_[i-1] * (2*i - 1);
47 }
48
49 }
50
SHCalc(const int num_vals,shared_ptr<SHCalcConstants> _consts)51 SHCalc::SHCalc(const int num_vals, shared_ptr<SHCalcConstants> _consts)
52 :numVals_(num_vals), P_( num_vals, num_vals), _consts_(_consts),
53 Y_( num_vals, num_vals)
54 {
55 assert (_consts_->get_n() == numVals_);
56 }
57
58 /*
59
60 Calculate the Legendre polynomial for the input theta using the
61 recursion functions for the polynomials, which are as follows:
62 Pl,l (x) = (-1)^l * (2l-1)!! * (1-x^2)^(l/2) (1)
63 Pl,l+1 (x) = x * (2l+1) * Pl,l(x) (2)
64 Pl,m (x) = x * (2l-1)/(l-m) * Pl-1,m(x) - (l+m-1)/(l-m) * Pl-2,m(x) (3)
65
66 This sets the member P_ to the results
67 */
calc_legendre(const double theta)68 void SHCalc::calc_legendre(const double theta)
69 {
70 double x = cos(theta);
71 P_.set_val(0, 0, 1.0); // base value for recursion
72 P_.set_val(1, 0, x);
73
74 int l, m;
75 double val = 0.0;
76 for (l = 0; l < numVals_; l++)
77 {
78 for (m = 0; m <= l; m++)
79 {
80 if ((l == 0 && m == 0) || (l == 1 && m == 0)) continue; //skip base values
81 else if (l == m)
82 {
83 double dblL = (double) l;
84 val = pow(-1.0, dblL) * _consts_->get_dub_fac_val(l)
85 * pow(1.0-x*x, dblL/2.0); // (1) in doc string
86 }
87 else if (m == l + 1)
88 {
89 val = x * (2*l + 1) * P_(l, l); // (2)
90 }
91 else if (m < l)
92 {
93 val = _consts_->get_leg_consts1_val(l, m) * x * P_(l-1, m); // (3)
94 val -= _consts_->get_leg_consts2_val(l, m) * P_(l-2, m); // (3)
95 }
96 P_.set_val(l, m, val);
97 }
98 }
99 }
100
101
102 /*
103 Return the results of the legendre calculation for an n, m.
104 */
get_legendre_result(int n,int m)105 double SHCalc::get_legendre_result( int n, int m )
106 {
107 return P_( n, m);
108 }
109
110 /*
111 Calculate the spherical harmonics according to the equation:
112
113 Y_(n,m)(theta, phi) = (-1)^m * sqrt((n-m)! / (n + m)!)
114 * P_(n,m)(cos(theta)) * exp(i*m*phi)
115 where P_(n, m) are the associated Legendre polynomials.
116
117 */
calc_sh(const double theta,const double phi)118 void SHCalc::calc_sh(const double theta, const double phi)
119 {
120 calc_legendre(theta); // first calculate legendre
121 cmplx iu (0, 1.0); // complex unit
122 int n, m;
123 cmplx val, mcomp;
124 double shc; // constant value
125 for (n = 0; n < numVals_; n++)
126 {
127 for (m = 0; m < numVals_; m++)
128 {
129 shc = _consts_->get_sh_consts_val(n, m);
130 double dblM = (double) m;
131 mcomp = complex<double> (dblM, 0.0);
132 val = pow(-1.0, dblM) * shc * P_(n, m) * exp(iu * mcomp * phi);
133 Y_.set_val(n, m, val);
134 }
135 }
136
137 }
138
139 /*
140 Return the results of the spherical harmonic calculation for an n, m.
141 If m is negative, then we return the complex conjugate of the calculated
142 value for the positive value of m
143 */
get_result(const int n,const int m)144 cmplx SHCalc::get_result(const int n, const int m)
145 {
146 if (m < 0)
147 {
148 return conj(Y_(n, -m)); // complex conjugate
149 }
150 else
151 {
152 return Y_(n, m);
153 }
154 }
155
156
157
calc_and_add(Pt p,shared_ptr<SHCalc> sh)158 void PreCalcSH::calc_and_add(Pt p, shared_ptr<SHCalc> sh)
159 {
160 sh->calc_sh(p.theta(), p.phi());
161 map_[p] = sh->get_full_result();
162 }
163
get_sh(Pt p,int n,int m)164 cmplx PreCalcSH::get_sh(Pt p, int n, int m)
165 {
166 if (m < 0) return conj(map_[p](n, -m)); // complex conjugate
167 else return map_[p](n, m);
168 }
169
170
171