diff --git a/bfps/cpp/field.cpp b/bfps/cpp/field.cpp
index 0e05983875d3288987c8af153b961e163cac711f..bbf910abde5377b1fe31b028af8a3fe00835fd4a 100644
--- a/bfps/cpp/field.cpp
+++ b/bfps/cpp/field.cpp
@@ -1223,19 +1223,19 @@ int joint_rspace_PDF(
             ptrdiff_t *local_histm = local_histm_threaded.getMine();
 
             double mag1, mag2;
-            mag1 = 0;
-            mag2 = 0;
-
+            mag1 = 0.0;
+            mag2 = 0.0;
             for (unsigned int i=0; i<3; i++)
             {
                 double val1 = f1->rval(rindex, i);
                 mag1 += val1*val1;
                 int bin1 = int(floor((val1 + max_f1_estimate[i])/bin1size[i]));
+                mag2 = 0.0;
                 for (unsigned int j=0; j<3; j++)
                 {
-                    double val2 = f2->rval(rindex, i);
+                    double val2 = f2->rval(rindex, j);
                     mag2 += val2*val2;
-                    int bin2 = int(floor((val2 + max_f2_estimate[i])/bin2size[i]));
+                    int bin2 = int(floor((val2 + max_f2_estimate[j])/bin2size[j]));
                     if ((bin1 >= 0 && bin1 < nbins) &&
                         (bin2 >= 0 && bin2 < nbins))
                         local_histc[(bin1*nbins + bin2)*9 + i*3 + j]++;
diff --git a/tests/PP/test_joint_acc_vel.py b/tests/PP/test_joint_acc_vel.py
index d623af28efea50800b9a3f27d4632f73e0af1a34..e61ec32ed1d5ad5309a122f3fca33223e4cebde3 100644
--- a/tests/PP/test_joint_acc_vel.py
+++ b/tests/PP/test_joint_acc_vel.py
@@ -11,10 +11,25 @@ def main():
     df.close()
 
     f = plt.figure()
-    a = f.add_subplot(111)
+    a = f.add_subplot(211)
     a.plot(acc_hist[0, :, :3])
     a.plot(np.sum(acc_vel_histc[0, :, :, :, 0], axis = 1), dashes = (1, 1))
-    f.savefig('sanity_test.pdf')
+    a = f.add_subplot(212)
+    a.plot(acc_hist[0, :, 3])
+    a.plot(np.sum(acc_vel_histm[0, :, :], axis = 1), dashes = (1, 1))
+    f.tight_layout()
+    f.savefig('sanity_test_acceleration.pdf')
+    plt.close(f)
+
+    f = plt.figure()
+    a = f.add_subplot(211)
+    a.plot(vel_hist[0, :, :3])
+    a.plot(np.sum(acc_vel_histc[0, :, :, 0, :], axis = 0), dashes = (1, 1))
+    a = f.add_subplot(212)
+    a.plot(vel_hist[0, :, 3])
+    a.plot(np.sum(acc_vel_histm[0, :, :], axis = 0), dashes = (1, 1))
+    f.tight_layout()
+    f.savefig('sanity_test_velocity.pdf')
     plt.close(f)
     return None