diff --git a/qpms/ewald.c b/qpms/ewald.c index 07f3892..e35942c 100644 --- a/qpms/ewald.c +++ b/qpms/ewald.c @@ -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 diff --git a/tests/ewalds.c b/tests/ewalds.c index faf27ec..3bc0eda 100644 --- a/tests/ewalds.c +++ b/tests/ewalds.c @@ -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];