1
- // CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
2
1
#include " CartesianR6PoissonCircular.h"
3
2
#include < stdlib.h>
4
3
#include < math.h>
4
+ #include < stdint.h>
5
5
6
6
7
7
/* ........................................*/
@@ -133,22 +133,22 @@ void CartesianR6PoissonCircular::J_tt(double r, std::vector<double> const& theta
133
133
/* ........................................*/
134
134
double CartesianR6PoissonCircular::J_xs (double r, double theta, double unused_1, double unused_2, double Rmax) const
135
135
{
136
- return (-pow (sin (theta), 2.0 )) / cos (theta) + pow (cos (theta), (double )((-1 )));
136
+ return (-pow (sin (theta), 2.0 )) / cos (theta) + pow (cos (theta), (double )((-INT64_C ( 1 ) )));
137
137
}
138
138
/* ........................................*/
139
139
void CartesianR6PoissonCircular::J_xs (std::vector<double > const & r, double theta, double unused_1, double unused_2, double Rmax, std::vector<double >& sol) const
140
140
{
141
141
for (std::size_t i=0 ; i < sol.size (); ++i)
142
142
{
143
- sol[i] = (-pow (sin (theta), 2.0 )) / cos (theta) + pow (cos (theta), (double )((-1 )));
143
+ sol[i] = (-pow (sin (theta), 2.0 )) / cos (theta) + pow (cos (theta), (double )((-INT64_C ( 1 ) )));
144
144
}
145
145
}
146
146
/* ........................................*/
147
147
void CartesianR6PoissonCircular::J_xs (double r, std::vector<double > const & theta, double unused_1, double unused_2, double Rmax, std::vector<double >& sol, std::vector<double >& sin_theta, std::vector<double >& cos_theta) const
148
148
{
149
149
for (std::size_t i=0 ; i < sol.size (); ++i)
150
150
{
151
- sol[i] = (-pow (sin_theta[i], 2.0 )) / cos_theta[i] + pow (cos_theta[i], (double )((-1 )));
151
+ sol[i] = (-pow (sin_theta[i], 2.0 )) / cos_theta[i] + pow (cos_theta[i], (double )((-INT64_C ( 1 ) )));
152
152
}
153
153
}
154
154
/* ........................................*/
@@ -217,22 +217,22 @@ void CartesianR6PoissonCircular::J_yt(double r, std::vector<double> const& theta
217
217
/* ........................................*/
218
218
double CartesianR6PoissonCircular::rho_glob (double r, double theta, double unused_1, double unused_2, double Rmax) const
219
219
{
220
- return 0.0 ;
220
+ return (-((-1.6384) * (M_PI * M_PI) * (r/Rmax) * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 6.0) * pow(sin(theta), 2.0) * sin(2.0 * M_PI * (r/Rmax) * sin(theta)) * cos(2.0 * M_PI * (r/Rmax) * cos(theta)) + 3.2768 * (M_PI * M_PI) * (r/Rmax) * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 6.0) * sin(theta) * sin(2.0 * M_PI * (r/Rmax) * cos(theta)) * cos(theta) * cos(2.0 * M_PI * (r/Rmax) * sin(theta)) - 1.6384 * (M_PI * M_PI) * (r/Rmax) * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 6.0) * sin(2.0 * M_PI * (r/Rmax) * sin(theta)) * pow(cos(theta), 2.0) * cos(2.0 * M_PI * (r/Rmax) * cos(theta)) + 1.0 * (r/Rmax) * ((-1.6384) * (M_PI * M_PI) * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 6.0) * pow(sin(theta), 2.0) * sin(2.0 * M_PI * (r/Rmax) * sin(theta)) * cos(2.0 * M_PI * (r/Rmax) * cos(theta)) - 3.2768 * (M_PI * M_PI) * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 6.0) * sin(theta) * sin(2.0 * M_PI * (r/Rmax) * cos(theta)) * cos(theta) * cos(2.0 * M_PI * (r/Rmax) * sin(theta)) - 1.6384 * (M_PI * M_PI) * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 6.0) * sin(2.0 * M_PI * (r/Rmax) * sin(theta)) * pow(cos(theta), 2.0) * cos(2.0 * M_PI * (r/Rmax) * cos(theta)) + 9.8304 * M_PI * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 5.0) * sin(theta) * cos(2.0 * M_PI * (r/Rmax) * sin(theta)) * cos(2.0 * M_PI * (r/Rmax) * cos(theta)) - 9.8304 * M_PI * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 5.0) * sin(2.0 * M_PI * (r/Rmax) * sin(theta)) * sin(2.0 * M_PI * (r/Rmax) * cos(theta)) * cos(theta) + 12.288 * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 4.0) * sin(2.0 * M_PI * (r/Rmax) * sin(theta)) * cos(2.0 * M_PI * (r/Rmax) * cos(theta)) + 9.8304 * M_PI * pow(((r/Rmax) - 1.0), 5.0) * pow(((r/Rmax) + 1.0), 6.0) * sin(theta) * cos(2.0 * M_PI * (r/Rmax) * sin(theta)) * cos(2.0 * M_PI * (r/Rmax) * cos(theta)) - 9.8304 * M_PI * pow(((r/Rmax) - 1.0), 5.0) * pow(((r/Rmax) + 1.0), 6.0) * sin(2.0 * M_PI * (r/Rmax) * sin(theta)) * sin(2.0 * M_PI * (r/Rmax) * cos(theta)) * cos(theta) + 29.4912 * pow(((r/Rmax) - 1.0), 5.0) * pow(((r/Rmax) + 1.0), 5.0) * sin(2.0 * M_PI * (r/Rmax) * sin(theta)) * cos(2.0 * M_PI * (r/Rmax) * cos(theta)) + 12.288 * pow(((r/Rmax) - 1.0), 4.0) * pow(((r/Rmax) + 1.0), 6.0) * sin(2.0 * M_PI * (r/Rmax) * sin(theta)) * cos(2.0 * M_PI * (r/Rmax) * cos(theta))) + 2.4576 * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 5.0) * sin(2.0 * M_PI * (r/Rmax) * sin(theta)) * cos(2.0 * M_PI * (r/Rmax) * cos(theta)) + 2.4576 * pow(((r/Rmax) - 1.0), 5.0) * pow(((r/Rmax) + 1.0), 6.0) * sin(2.0 * M_PI * (r/Rmax) * sin(theta)) * cos(2.0 * M_PI * (r/Rmax) * cos(theta)))) / (r/Rmax);
221
221
}
222
222
/* ........................................*/
223
223
void CartesianR6PoissonCircular::rho_glob (std::vector<double > const & r, double theta, double unused_1, double unused_2, double Rmax, std::vector<double >& sol) const
224
224
{
225
225
for (std::size_t i=0 ; i < sol.size (); ++i)
226
226
{
227
- sol[i] = 0.0 ;
227
+ sol[i] = (-((-1.6384) * (M_PI * M_PI) * (r[i]/Rmax) * pow(((r[i]/Rmax) - 1.0), 6.0) * pow(((r[i]/Rmax) + 1.0), 6.0) * pow(sin(theta), 2.0) * sin(2.0 * M_PI * (r[i]/Rmax) * sin(theta)) * cos(2.0 * M_PI * (r[i]/Rmax) * cos(theta)) + 3.2768 * (M_PI * M_PI) * (r[i]/Rmax) * pow(((r[i]/Rmax) - 1.0), 6.0) * pow(((r[i]/Rmax) + 1.0), 6.0) * sin(theta) * sin(2.0 * M_PI * (r[i]/Rmax) * cos(theta)) * cos(theta) * cos(2.0 * M_PI * (r[i]/Rmax) * sin(theta)) - 1.6384 * (M_PI * M_PI) * (r[i]/Rmax) * pow(((r[i]/Rmax) - 1.0), 6.0) * pow(((r[i]/Rmax) + 1.0), 6.0) * sin(2.0 * M_PI * (r[i]/Rmax) * sin(theta)) * pow(cos(theta), 2.0) * cos(2.0 * M_PI * (r[i]/Rmax) * cos(theta)) + 1.0 * (r[i]/Rmax) * ((-1.6384) * (M_PI * M_PI) * pow(((r[i]/Rmax) - 1.0), 6.0) * pow(((r[i]/Rmax) + 1.0), 6.0) * pow(sin(theta), 2.0) * sin(2.0 * M_PI * (r[i]/Rmax) * sin(theta)) * cos(2.0 * M_PI * (r[i]/Rmax) * cos(theta)) - 3.2768 * (M_PI * M_PI) * pow(((r[i]/Rmax) - 1.0), 6.0) * pow(((r[i]/Rmax) + 1.0), 6.0) * sin(theta) * sin(2.0 * M_PI * (r[i]/Rmax) * cos(theta)) * cos(theta) * cos(2.0 * M_PI * (r[i]/Rmax) * sin(theta)) - 1.6384 * (M_PI * M_PI) * pow(((r[i]/Rmax) - 1.0), 6.0) * pow(((r[i]/Rmax) + 1.0), 6.0) * sin(2.0 * M_PI * (r[i]/Rmax) * sin(theta)) * pow(cos(theta), 2.0) * cos(2.0 * M_PI * (r[i]/Rmax) * cos(theta)) + 9.8304 * M_PI * pow(((r[i]/Rmax) - 1.0), 6.0) * pow(((r[i]/Rmax) + 1.0), 5.0) * sin(theta) * cos(2.0 * M_PI * (r[i]/Rmax) * sin(theta)) * cos(2.0 * M_PI * (r[i]/Rmax) * cos(theta)) - 9.8304 * M_PI * pow(((r[i]/Rmax) - 1.0), 6.0) * pow(((r[i]/Rmax) + 1.0), 5.0) * sin(2.0 * M_PI * (r[i]/Rmax) * sin(theta)) * sin(2.0 * M_PI * (r[i]/Rmax) * cos(theta)) * cos(theta) + 12.288 * pow(((r[i]/Rmax) - 1.0), 6.0) * pow(((r[i]/Rmax) + 1.0), 4.0) * sin(2.0 * M_PI * (r[i]/Rmax) * sin(theta)) * cos(2.0 * M_PI * (r[i]/Rmax) * cos(theta)) + 9.8304 * M_PI * pow(((r[i]/Rmax) - 1.0), 5.0) * pow(((r[i]/Rmax) + 1.0), 6.0) * sin(theta) * cos(2.0 * M_PI * (r[i]/Rmax) * sin(theta)) * cos(2.0 * M_PI * (r[i]/Rmax) * cos(theta)) - 9.8304 * M_PI * pow(((r[i]/Rmax) - 1.0), 5.0) * pow(((r[i]/Rmax) + 1.0), 6.0) * sin(2.0 * M_PI * (r[i]/Rmax) * sin(theta)) * sin(2.0 * M_PI * (r[i]/Rmax) * cos(theta)) * cos(theta) + 29.4912 * pow(((r[i]/Rmax) - 1.0), 5.0) * pow(((r[i]/Rmax) + 1.0), 5.0) * sin(2.0 * M_PI * (r[i]/Rmax) * sin(theta)) * cos(2.0 * M_PI * (r[i]/Rmax) * cos(theta)) + 12.288 * pow(((r[i]/Rmax) - 1.0), 4.0) * pow(((r[i]/Rmax) + 1.0), 6.0) * sin(2.0 * M_PI * (r[i]/Rmax) * sin(theta)) * cos(2.0 * M_PI * (r[i]/Rmax) * cos(theta))) + 2.4576 * pow(((r[i]/Rmax) - 1.0), 6.0) * pow(((r[i]/Rmax) + 1.0), 5.0) * sin(2.0 * M_PI * (r[i]/Rmax) * sin(theta)) * cos(2.0 * M_PI * (r[i]/Rmax) * cos(theta)) + 2.4576 * pow(((r[i]/Rmax) - 1.0), 5.0) * pow(((r[i]/Rmax) + 1.0), 6.0) * sin(2.0 * M_PI * (r[i]/Rmax) * sin(theta)) * cos(2.0 * M_PI * (r[i]/Rmax) * cos(theta)))) / (r[i]/Rmax);
228
228
}
229
229
}
230
230
/* ........................................*/
231
231
void CartesianR6PoissonCircular::rho_glob (double r, std::vector<double > const & theta, double unused_1, double unused_2, double Rmax, std::vector<double >& sol, std::vector<double >& sin_theta, std::vector<double >& cos_theta) const
232
232
{
233
233
for (std::size_t i=0 ; i < sol.size (); ++i)
234
234
{
235
- sol[i] = 0.0 ;
235
+ sol[i] = (-((-1.6384) * (M_PI * M_PI) * (r/Rmax) * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 6.0) * pow(sin_theta[i], 2.0) * sin(2.0 * M_PI * (r/Rmax) * sin_theta[i]) * cos(2.0 * M_PI * (r/Rmax) * cos_theta[i]) + 3.2768 * (M_PI * M_PI) * (r/Rmax) * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 6.0) * sin_theta[i] * sin(2.0 * M_PI * (r/Rmax) * cos_theta[i]) * cos_theta[i] * cos(2.0 * M_PI * (r/Rmax) * sin_theta[i]) - 1.6384 * (M_PI * M_PI) * (r/Rmax) * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 6.0) * sin(2.0 * M_PI * (r/Rmax) * sin_theta[i]) * pow(cos_theta[i], 2.0) * cos(2.0 * M_PI * (r/Rmax) * cos_theta[i]) + 1.0 * (r/Rmax) * ((-1.6384) * (M_PI * M_PI) * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 6.0) * pow(sin_theta[i], 2.0) * sin(2.0 * M_PI * (r/Rmax) * sin_theta[i]) * cos(2.0 * M_PI * (r/Rmax) * cos_theta[i]) - 3.2768 * (M_PI * M_PI) * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 6.0) * sin_theta[i] * sin(2.0 * M_PI * (r/Rmax) * cos_theta[i]) * cos_theta[i] * cos(2.0 * M_PI * (r/Rmax) * sin_theta[i]) - 1.6384 * (M_PI * M_PI) * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 6.0) * sin(2.0 * M_PI * (r/Rmax) * sin_theta[i]) * pow(cos_theta[i], 2.0) * cos(2.0 * M_PI * (r/Rmax) * cos_theta[i]) + 9.8304 * M_PI * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 5.0) * sin_theta[i] * cos(2.0 * M_PI * (r/Rmax) * sin_theta[i]) * cos(2.0 * M_PI * (r/Rmax) * cos_theta[i]) - 9.8304 * M_PI * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 5.0) * sin(2.0 * M_PI * (r/Rmax) * sin_theta[i]) * sin(2.0 * M_PI * (r/Rmax) * cos_theta[i]) * cos_theta[i] + 12.288 * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 4.0) * sin(2.0 * M_PI * (r/Rmax) * sin_theta[i]) * cos(2.0 * M_PI * (r/Rmax) * cos_theta[i]) + 9.8304 * M_PI * pow(((r/Rmax) - 1.0), 5.0) * pow(((r/Rmax) + 1.0), 6.0) * sin_theta[i] * cos(2.0 * M_PI * (r/Rmax) * sin_theta[i]) * cos(2.0 * M_PI * (r/Rmax) * cos_theta[i]) - 9.8304 * M_PI * pow(((r/Rmax) - 1.0), 5.0) * pow(((r/Rmax) + 1.0), 6.0) * sin(2.0 * M_PI * (r/Rmax) * sin_theta[i]) * sin(2.0 * M_PI * (r/Rmax) * cos_theta[i]) * cos_theta[i] + 29.4912 * pow(((r/Rmax) - 1.0), 5.0) * pow(((r/Rmax) + 1.0), 5.0) * sin(2.0 * M_PI * (r/Rmax) * sin_theta[i]) * cos(2.0 * M_PI * (r/Rmax) * cos_theta[i]) + 12.288 * pow(((r/Rmax) - 1.0), 4.0) * pow(((r/Rmax) + 1.0), 6.0) * sin(2.0 * M_PI * (r/Rmax) * sin_theta[i]) * cos(2.0 * M_PI * (r/Rmax) * cos_theta[i])) + 2.4576 * pow(((r/Rmax) - 1.0), 6.0) * pow(((r/Rmax) + 1.0), 5.0) * sin(2.0 * M_PI * (r/Rmax) * sin_theta[i]) * cos(2.0 * M_PI * (r/Rmax) * cos_theta[i]) + 2.4576 * pow(((r/Rmax) - 1.0), 5.0) * pow(((r/Rmax) + 1.0), 6.0) * sin(2.0 * M_PI * (r/Rmax) * sin_theta[i]) * cos(2.0 * M_PI * (r/Rmax) * cos_theta[i]))) / (r/Rmax);
236
236
}
237
237
}
238
238
/* ........................................*/
@@ -259,14 +259,14 @@ void CartesianR6PoissonCircular::rho_pole(double r, std::vector<double> const& t
259
259
/* ........................................*/
260
260
double CartesianR6PoissonCircular::coeffs1 (double r, double Rmax) const
261
261
{
262
- return 0 .0 ;
262
+ return 1 .0 ;
263
263
}
264
264
/* ........................................*/
265
265
void CartesianR6PoissonCircular::coeffs1 (std::vector<double > const & r, double Rmax, std::vector<double >& sol) const
266
266
{
267
267
for (std::size_t i=0 ; i < sol.size (); ++i)
268
268
{
269
- sol[i] = 0 .0 ;
269
+ sol[i] = 1 .0 ;
270
270
}
271
271
}
272
272
/* ........................................*/
0 commit comments