QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

konstantisace; dudopráce 

Former-commit-id: a146fe83f9dbf17044cc723c325ba70230969e44
This commit is contained in:
Marek Nečada 2018-09-10 05:10:13 +00:00
parent bad8105c9b
commit 0426524842
2 changed files with 13 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
qpms/ewald.c
View File

@ -286,7 +286,7 @@ int ewald32_sigma_short_shiftedpoints(
gsl_integration_workspace_alloc(INTEGRATION_WORKSPACE_LIMIT);
// Manual init of the ewald summation targets
complex double *target_c = calloc(nelem_sc, sizeof(complex double));
complex double * const target_c = calloc(nelem_sc, sizeof(complex double));
memset(target, 0, nelem_sc * sizeof(complex double));
double *err_c = NULL;
if (err) {
@ -297,17 +297,16 @@ int ewald32_sigma_short_shiftedpoints(
// CHOOSE POINT BEGIN
for (size_t i = 0; i < npoints; ++i) { // BEGIN POINT LOOP
point2d Rpq_shifted = Rpoints_plus_particle_shift[i];
double r_pq_shifted = cart2norm(Rpq_shifted);
const point2d Rpq_shifted = Rpoints_plus_particle_shift[i];
const double r_pq_shifted = cart2norm(Rpq_shifted);
// CHOOSE POINT END
// This should be imaginary, but gcc does not support:
double Rpq_shifted_arg = atan2(Rpq_shifted.x, Rpq_shifted.y); // POINT-DEPENDENT
complex double e_beta_Rpq = cexp(I*cart2_dot(beta, Rpq_shifted)); // POINT-DEPENDENT
const double Rpq_shifted_arg = atan2(Rpq_shifted.y, Rpq_shifted.x); // POINT-DEPENDENT
const complex double e_beta_Rpq = cexp(I*cart2_dot(beta, Rpq_shifted)); // POINT-DEPENDENT
for(qpms_l_t n = 0; n <= lMax; ++n) {
double complex prefacn = - I * pow(2./k, n+1) * M_2_SQRTPI / 2; // TODO put outside the R-loop and multiply in the end
double R_pq_pown = pow(r_pq_shifted, n);
const double complex prefacn = - I * pow(2./k, n+1) * M_2_SQRTPI / 2; // TODO put outside the R-loop and multiply in the end
const double R_pq_pown = pow(r_pq_shifted, n);
// TODO the integral here
double intres, interr;
int retval = ewald32_sr_integral(r_pq_shifted, k, n, eta,
@ -316,9 +315,9 @@ int ewald32_sigma_short_shiftedpoints(
for (qpms_m_t m = -n; m <= n; ++m){
if((m+n) % 2 != 0) // odd coefficients are zero.
continue; // nothing needed, already done by memset
complex double e_imf = cexp(I*m*Rpq_shifted_arg);
double leg = c->legendre0[gsl_sf_legendre_array_index(n, abs(m))];
qpms_y_t y = qpms_mn2y_sc(m,n);
const complex double e_imf = cexp(I*m*Rpq_shifted_arg);
const double leg = c->legendre0[gsl_sf_legendre_array_index(n, abs(m))];
const qpms_y_t y = qpms_mn2y_sc(m,n);
if(err)
kahanadd(err + y, err_c + y, cabs(leg * (prefacn / I) * R_pq_pown
* interr)); // TODO include also other errors

6
tests/ewalds.c
View File

@ -39,7 +39,7 @@ typedef struct ewaldtest_triang_results {
ewaldtest_triang_params paramslist[] = {
{ 3, {1.1, 0.23}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_HORIZONTAL},
{ 3, {1.1, 0.23}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL},
{ 3, {1.1, 0.23}, 2.3, 0.97, 0.5, 30, 30, 1., TRIANGULAR_VERTICAL},
{ 3, {1.1, 0.23}, 2.3, 0.97, 0.5, 30, 30, -1., TRIANGULAR_VERTICAL},
{ 3, {1.1, 0.23}, 2.3, 0.97, 0.9, 30, 30, 1., TRIANGULAR_VERTICAL},
{ 3, {1.1, 0.23}, 2.3, 0.97, 1.3, 30, 30, 1., TRIANGULAR_VERTICAL},
{ 6, {1.1, 0.23}, 2.3, 0.97, 1.9, 30, 30, 1., TRIANGULAR_VERTICAL},
@ -75,7 +75,7 @@ int main() {
qpms_y_t y = qpms_mn2y_sc(m,n);
qpms_y_t y_conj = qpms_mn2y_sc(-m,n);
// y n m sigma_total (err), regsigmas_416 regsigmas_415_recon
printf("%zd %d %d: T:%.16g%+.16gj(%.3g) L:%.16g%+.16gj(%.3g) S:%.16g%+.16gj(%.3g) \n| predict %.16g%+.16gj \n| actual %.16g%+.16gj\n",
printf("%zd %d %d: T:%.16g%+.16gj(%.3g) L:%.16g%+.16gj(%.3g) S:%.16g%+.16gj(%.3g) \n| predict %.16g%+.16gj \n| actual %.16g%+.16gj\n",
y, n, m, creal(r->sigmas_total[y]), cimag(r->sigmas_total[y]),
r->err_sigmas_total[y],
creal(r->sigmas_long[y]), cimag(r->sigmas_long[y]),
@ -109,7 +109,7 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) {
triangular_lattice_gen_t *Klg = triangular_lattice_gen_init(K_len, reverseTriangularLatticeOrientation(p.orientation), true, 0);
triangular_lattice_gen_extend_to_r(Klg, p.maxK + K_len);
point2d *Rpoints = Rlg->ps.base, *Kpoints = Klg->ps.base;
point2d *Rpoints = Rlg->ps.base; //point2d *Kpoints = Klg->ps.base;
size_t nR = Rlg->ps.r_offsets[Rlg->ps.nrs],
nK = Klg->ps.r_offsets[Klg->ps.nrs];