Merge branch 'feature/p2p_rod_typos' into develop

cpp/particles/p2p/p2p_ghost_collisions.hpp

@@ -218,7 +218,7 @@ public:
                              real_number /*rhs_part2*/[],
                              const real_number dist_pow2,
                              const real_number /*cutoff*/,
-                             const real_number xseparation,
+                             const real_number xseparation, /* This separation is x1-x2 */
                              const real_number yseparation,
                              const real_number zseparation){
         switch(this->current_particle_shape)
@@ -245,13 +245,13 @@ public:
 
                         /* p and q are the orientation vectors of the first and second particles. */
                         /* pq, xp, xq are scalar products of these vectors with x, relative position */
-                        pq = pos_part1[3] * pos_part2[3] + pos_part1[4] * pos_part2[4] + pos_part1[5] * pos_part2[5];
-                        xp = x * pos_part1[3] + y * pos_part1[4] + z * pos_part1[5];
-                        xq = x * pos_part2[3] + y * pos_part2[4] + z * pos_part2[5];
+                        pq = pos_part1[IDXC_X+3] * pos_part2[IDXC_X+3] + pos_part1[IDXC_Y+3] * pos_part2[IDXC_Y+3] + pos_part1[IDXC_Z+3] * pos_part2[IDXC_Z+3];
+                        xp = x * pos_part1[IDXC_X+3] + y * pos_part1[IDXC_Y+3] + z * pos_part1[IDXC_Z+3];
+                        xq = x * pos_part2[IDXC_X+3] + y * pos_part2[IDXC_Y+3] + z * pos_part2[IDXC_Z+3];
                         /* t and s parametrize the two rods. Find min distance: */
                         assert(this->cylinder_length > 0);
-                        t = 2.0/(this->cylinder_length*(pq*pq-1.0))*(-xp+pq*xq);
-                        s = 2.0/(this->cylinder_length*(pq*pq-1.0))*(-pq*xp+xq);
+                        t = 2.0/(this->cylinder_length*(pq*pq-1.0))*(xp-pq*xq);
+                        s = 2.0/(this->cylinder_length*(pq*pq-1.0))*(pq*xp-xq);
                         /* Test if -1<s<1 and -1<t<1 */
                         if( abs(t)<=1.0 and abs(s)<=1.0 )
                         {
@@ -266,7 +266,7 @@ public:
                             min_distance = this->cylinder_length;
                             /* t fixed at 1, find min along s */
                             t = 1.0;
-                            s = t*pq-2.0/this->cylinder_length*xq;
+                            s = t*pq+2.0/this->cylinder_length*xq;
                             if( abs(s)>1.0 ) { s = s / abs(s) ;}
                             x_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_X+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_X+3]+x;
                             y_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Y+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Y+3]+y;
@@ -274,7 +274,7 @@ public:
                             min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance );
                             /* t fixed at -1, find min along s */
                             t = -1.0;
-                            s = t*pq-2.0/this->cylinder_length*xq;
+                            s = t*pq+2.0/this->cylinder_length*xq;
                             if( abs(s)>1.0 ) { s = s / abs(s) ;}
                             x_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_X+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_X+3]+x;
                             y_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Y+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Y+3]+y;
@@ -282,7 +282,7 @@ public:
                             min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance );
                             /* s fixed at 1, find min along t */
                             s = 1.0;
-                            t = s*pq+2.0/this->cylinder_length*xp;
+                            t = s*pq-2.0/this->cylinder_length*xp;
                             if( abs(t)>1.0 ) { t = t / abs(t) ;}
                             x_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_X+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_X+3]+x;
                             y_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Y+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Y+3]+y;
@@ -290,7 +290,7 @@ public:
                             min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance );
                             /* s fixed at -1, find min along t */
                             s = -1.0;
-                            t = s*pq+2.0/this->cylinder_length*xp;
+                            t = s*pq-2.0/this->cylinder_length*xp;
                             if( abs(t)>1.0 ) { t = t / abs(t) ;}
                             x_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_X+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_X+3]+x;
                             y_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Y+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Y+3]+y;