1212
1313namespace gpl {
1414
15- FFT::FFT (int binCntX , int binCntY , float binSizeX , float binSizeY )
16- : binCntX_(binCntX ),
17- binCntY_ (binCntY ),
18- binSizeX_(binSizeX ),
19- binSizeY_(binSizeY )
15+ FFT::FFT (int bin_cnt_x , int bin_cnt_y , float bin_size_x , float bin_size_y )
16+ : bin_cnt_X_(bin_cnt_x ),
17+ bin_cnt_y_ (bin_cnt_y ),
18+ bin_size_x_(bin_size_x ),
19+ bin_size_y_(bin_size_y )
2020{
21- binDensity_ = new float *[binCntX_ ];
22- electroPhi_ = new float *[binCntX_ ];
23- electroForceX_ = new float *[binCntX_ ];
24- electroForceY_ = new float *[binCntX_ ];
25-
26- for (int i = 0 ; i < binCntX_ ; i++) {
27- binDensity_ [i] = new float [binCntY_ ];
28- electroPhi_ [i] = new float [binCntY_ ];
29- electroForceX_ [i] = new float [binCntY_ ];
30- electroForceY_ [i] = new float [binCntY_ ];
31-
32- for (int j = 0 ; j < binCntY_ ; j++) {
33- binDensity_ [i][j] = electroPhi_ [i][j] = electroForceX_ [i][j]
34- = electroForceY_ [i][j] = 0 .0f ;
21+ bin_density_ = new float *[bin_cnt_X_ ];
22+ electro_phi_ = new float *[bin_cnt_X_ ];
23+ electro_field_x_ = new float *[bin_cnt_X_ ];
24+ electro_field_y_ = new float *[bin_cnt_X_ ];
25+
26+ for (int i = 0 ; i < bin_cnt_X_ ; i++) {
27+ bin_density_ [i] = new float [bin_cnt_y_ ];
28+ electro_phi_ [i] = new float [bin_cnt_y_ ];
29+ electro_field_x_ [i] = new float [bin_cnt_y_ ];
30+ electro_field_y_ [i] = new float [bin_cnt_y_ ];
31+
32+ for (int j = 0 ; j < bin_cnt_y_ ; j++) {
33+ bin_density_ [i][j] = electro_phi_ [i][j] = electro_field_x_ [i][j]
34+ = electro_field_y_ [i][j] = 0 .0f ;
3535 }
3636 }
3737
38- csTable_ .resize (std::max (binCntX_, binCntY_ ) * 3 / 2 , 0 );
38+ cs_table_ .resize (std::max (bin_cnt_X_, bin_cnt_y_ ) * 3 / 2 , 0 );
3939
40- wx_.resize (binCntX_ , 0 );
41- wxSquare_ .resize (binCntX_ , 0 );
42- wy_.resize (binCntY_ , 0 );
43- wySquare_ .resize (binCntY_ , 0 );
40+ wx_.resize (bin_cnt_X_ , 0 );
41+ wx_square_ .resize (bin_cnt_X_ , 0 );
42+ wy_.resize (bin_cnt_y_ , 0 );
43+ wy_square_ .resize (bin_cnt_y_ , 0 );
4444
45- workArea_ .resize (round (sqrt (std::max (binCntX_, binCntY_ ))) + 2 , 0 );
45+ work_area_ .resize (round (sqrt (std::max (bin_cnt_X_, bin_cnt_y_ ))) + 2 , 0 );
4646
4747 constexpr auto kPi = std::numbers::pi_v<long double >;
4848
49- for (int i = 0 ; i < binCntX_ ; i++) {
50- wx_[i] = kPi * static_cast <float >(i) / static_cast <float >(binCntX_ );
51- wxSquare_ [i] = wx_[i] * wx_[i];
49+ for (int i = 0 ; i < bin_cnt_X_ ; i++) {
50+ wx_[i] = kPi * static_cast <float >(i) / static_cast <float >(bin_cnt_X_ );
51+ wx_square_ [i] = wx_[i] * wx_[i];
5252 }
5353
54- for (int i = 0 ; i < binCntY_ ; i++) {
55- wy_[i] = kPi * static_cast <float >(i) / static_cast <float >(binCntY_ )
56- * binSizeY_ / binSizeX_ ;
57- wySquare_ [i] = wy_[i] * wy_[i];
54+ for (int i = 0 ; i < bin_cnt_y_ ; i++) {
55+ wy_[i] = kPi * static_cast <float >(i) / static_cast <float >(bin_cnt_y_ )
56+ * bin_size_y_ / bin_size_x_ ;
57+ wy_square_ [i] = wy_[i] * wy_[i];
5858 }
5959}
6060
6161FFT ::~FFT ()
6262{
6363 using std::vector;
64- for (int i = 0 ; i < binCntX_ ; i++) {
65- delete[] binDensity_ [i];
66- delete[] electroPhi_ [i];
67- delete[] electroForceX_ [i];
68- delete[] electroForceY_ [i];
64+ for (int i = 0 ; i < bin_cnt_X_ ; i++) {
65+ delete[] bin_density_ [i];
66+ delete[] electro_phi_ [i];
67+ delete[] electro_field_x_ [i];
68+ delete[] electro_field_y_ [i];
6969 }
70- delete[] binDensity_ ;
71- delete[] electroPhi_ ;
72- delete[] electroForceX_ ;
73- delete[] electroForceY_ ;
70+ delete[] bin_density_ ;
71+ delete[] electro_phi_ ;
72+ delete[] electro_field_x_ ;
73+ delete[] electro_field_y_ ;
7474
75- csTable_ .clear ();
75+ cs_table_ .clear ();
7676 wx_.clear ();
77- wxSquare_ .clear ();
77+ wx_square_ .clear ();
7878 wy_.clear ();
79- wySquare_ .clear ();
79+ wy_square_ .clear ();
8080
81- workArea_ .clear ();
81+ work_area_ .clear ();
8282}
8383
8484void FFT::updateDensity (int x, int y, float density)
8585{
86- binDensity_ [x][y] = density;
86+ bin_density_ [x][y] = density;
8787}
8888
89- std::pair<float , float > FFT::getElectroForce (int x, int y) const
89+ std::pair<float , float > FFT::getElectroField (int x, int y) const
9090{
91- return std::make_pair (electroForceX_ [x][y], electroForceY_ [x][y]);
91+ return std::make_pair (electro_field_x_ [x][y], electro_field_y_ [x][y]);
9292}
9393
9494float FFT::getElectroPhi (int x, int y) const
9595{
96- return electroPhi_ [x][y];
96+ return electro_phi_ [x][y];
9797}
9898
9999void FFT::doFFT ()
100100{
101- ddct2d (binCntX_ ,
102- binCntY_ ,
101+ ddct2d (bin_cnt_X_ ,
102+ bin_cnt_y_ ,
103103 -1 ,
104- binDensity_ ,
104+ bin_density_ ,
105105 nullptr ,
106- workArea_ .data (),
107- csTable_ .data ());
106+ work_area_ .data (),
107+ cs_table_ .data ());
108108
109- for (int i = 0 ; i < binCntX_; i++) {
110- binDensity_[i][0 ] *= 0.5 ;
109+ // Normalizations required to perform the inverse operation
110+ for (int i = 1 ; i < bin_cnt_X_; i++) {
111+ bin_density_[i][0 ] *= 0.5 ;
111112 }
112-
113- for (int i = 0 ; i < binCntY_; i++) {
114- binDensity_[0 ][i] *= 0.5 ;
113+ for (int i = 1 ; i < bin_cnt_y_; i++) {
114+ bin_density_[0 ][i] *= 0.5 ;
115115 }
116-
117- for (int i = 0 ; i < binCntX_; i++) {
118- for (int j = 0 ; j < binCntY_; j++) {
119- binDensity_[i][j] *= 4.0 / binCntX_ / binCntY_;
116+ for (int i = 0 ; i < bin_cnt_X_; i++) {
117+ for (int j = 0 ; j < bin_cnt_y_; j++) {
118+ bin_density_[i][j] *= 4.0 / bin_cnt_X_ / bin_cnt_y_;
120119 }
121120 }
122121
123- for (int i = 0 ; i < binCntX_; i++) {
122+ // Solve the PDE in the new basis
123+ for (int i = 0 ; i < bin_cnt_X_; i++) {
124124 float wx = wx_[i];
125- float wx2 = wxSquare_ [i];
125+ float wx2 = wx_square_ [i];
126126
127- for (int j = 0 ; j < binCntY_ ; j++) {
127+ for (int j = 0 ; j < bin_cnt_y_ ; j++) {
128128 float wy = wy_[j];
129- float wy2 = wySquare_ [j];
129+ float wy2 = wy_square_ [j];
130130
131- float density = binDensity_ [i][j];
131+ float density = bin_density_ [i][j];
132132 float phi = 0 ;
133- float electroX = 0 , electroY = 0 ;
133+ float electro_x = 0 , electro_y = 0 ;
134134
135135 if (i == 0 && j == 0 ) {
136- phi = electroX = electroY = 0 .0f ;
136+ // Removes the DC component
137+ phi = electro_x = electro_y = 0 .0f ;
137138 } else {
138139 // ////////// lutong
139140 // denom =
@@ -147,36 +148,38 @@ void FFT::doFFT()
147148 // a_ey = a_phi * wy / 2.0 ;
148149 // /////////
149150 phi = density / (wx2 + wy2);
150- electroX = phi * wx;
151- electroY = phi * wy;
151+ electro_x = phi * wx;
152+ electro_y = phi * wy;
152153 }
153- electroPhi_[i][j] = phi;
154- electroForceX_[i][j] = electroX;
155- electroForceY_[i][j] = electroY;
154+
155+ electro_phi_[i][j] = phi;
156+ electro_field_x_[i][j] = electro_x;
157+ electro_field_y_[i][j] = electro_y;
156158 }
157159 }
160+
158161 // Inverse DCT
159- ddct2d (binCntX_ ,
160- binCntY_ ,
162+ ddct2d (bin_cnt_X_ ,
163+ bin_cnt_y_ ,
161164 1 ,
162- electroPhi_ ,
165+ electro_phi_ ,
163166 nullptr ,
164- workArea_ .data (),
165- csTable_ .data ());
166- ddsct2d (binCntX_ ,
167- binCntY_ ,
167+ work_area_ .data (),
168+ cs_table_ .data ());
169+ ddsct2d (bin_cnt_X_ ,
170+ bin_cnt_y_ ,
168171 1 ,
169- electroForceX_ ,
172+ electro_field_x_ ,
170173 nullptr ,
171- workArea_ .data (),
172- csTable_ .data ());
173- ddcst2d (binCntX_ ,
174- binCntY_ ,
174+ work_area_ .data (),
175+ cs_table_ .data ());
176+ ddcst2d (bin_cnt_X_ ,
177+ bin_cnt_y_ ,
175178 1 ,
176- electroForceY_ ,
179+ electro_field_y_ ,
177180 nullptr ,
178- workArea_ .data (),
179- csTable_ .data ());
181+ work_area_ .data (),
182+ cs_table_ .data ());
180183}
181184
182185} // namespace gpl
0 commit comments