@@ -75,7 +75,8 @@ namespace simple_physics_solver
7575// cout << "SimplePhysicsSolver::evaluate_internal_e_field(t=" << t << ")" << endl;
7676 double r = double (0.0 ),
7777 r3 = double (0.0 ),
78- dist = double (0.0 );
78+ dist2 = double (0.0 );
79+ double * dist = new double [DIM];
7980
8081 // computing forces on particle i
8182 for (size_t i = 0 ; i < num_particles; ++i) {
@@ -92,18 +93,19 @@ namespace simple_physics_solver
9293// cout << " skipping" << endl;
9394 continue ;
9495 }
95- dist = double (0.0 );
96+ dist2 = double (0.0 );
9697 for (size_t d = 0 ; d < DIM; ++d) {
97- dist += (positions[i * DIM + d] - positions[j * DIM + d]) * (positions[i * DIM + d] - positions[j * DIM + d]);
98+ dist[d] = positions[i * DIM + d] - positions[j * DIM + d];
99+ dist2 += dist[d] * dist[d];
98100 }
99101// cout << " dist = " << dist << endl;
100- r = sqrt (dist * dist + config->sigma2 );
102+ r = sqrt (dist2 + config->sigma2 );
101103// cout << " r = " << r << " (= sqrt(dist^2+" << config->sigma2 << ")" << endl;
102- phis[i] += charges[j] / r;
104+ // phis[i] += charges[j] / r;
103105
104106 r3 = r * r * r;
105107 for (size_t d = 0 ; d < DIM; ++d) {
106- exyz[i * DIM + d] += positions[j * DIM + d] / r3 * charges[j];
108+ exyz[i * DIM + d] += dist[ d] / r3 * charges[j];
107109// cout << " exyz[" << i * DIM + d << "] += " << positions[j * DIM + d] << " / " << r3 << " * " << charges[j] << "(q) => " << exyz[i * DIM + d] << endl;
108110 }
109111 }
@@ -113,6 +115,7 @@ namespace simple_physics_solver
113115// cout << endl
114116// << " phi_i = " << phis[i] << endl;
115117 }
118+ delete[] dist;
116119// cout << "DONE SimplePhysicsSolver::evaluate_internal_e_field(t=" << t << ")" << endl;
117120 }
118121
0 commit comments