Apply simple cubic bias correction in the side case

examples/plot_sorted.py

@@ -36,7 +36,7 @@ import os
 from scipy import stats
 from scipy import optimize
 
-dist_cutoff_ratio=1.5
+dist_cutoff_ratio=1.2
 ortho_max_fit = 0.2
 ortho_min_fit = 0.
 ortho_max = 2.5
@@ -79,14 +79,14 @@ limit_exact = ( dist_cutoff_ratio - 1. )
 print limit_exact    
 #names = ["side", "edge", "corner"]
 
-for orientation in range( 26 ):
-# for jjj in range(0,3):
-#     if jjj == 0:
-#         orientation = 0
-#     if jjj == 1:
-#         orientation = 1
-#     if jjj == 2:
-#         orientation = 4
+#for orientation in range( 26 ):
+for jjj in range(0,3):
+    if jjj == 0:
+        orientation = 0
+    if jjj == 1:
+        orientation = 1
+    if jjj == 2:
+        orientation = 4
          
     # Read Quickshed accelerations
     data=loadtxt( "interaction_dump_%d.dat"%orientation )
@@ -168,8 +168,8 @@ for orientation in range( 26 ):
     delta = sim / exact    
     # parabola_param1,pcov = optimize.curve_fit(parabola, xx, delta, p0=[1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
     # #print parabola_param1
-    # cubic_param1,pcov = optimize.curve_fit(cubic, xx, delta, p0=[1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
-    # print cubic_param1
+    cubic_param1,pcov = optimize.curve_fit(cubic, xx, delta, p0=[1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
+    print cubic_param1
     # quartic_param1,pcov = optimize.curve_fit(quartic, xx, delta, p0=[1,1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
     # #print quartic_param1
     # quintic_param1,pcov = optimize.curve_fit(quintic, xx, delta, p0=[1,1,1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
@@ -200,8 +200,8 @@ for orientation in range( 26 ):
     delta = sim / exact    
     # parabola_param2,pcov = optimize.curve_fit(parabola, xx, delta, p0=[1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)   
     # #print parabola_param2
-    # cubic_param2,pcov = optimize.curve_fit(cubic, xx, delta, p0=[1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)   
-    # print cubic_param2
+    cubic_param2,pcov = optimize.curve_fit(cubic, xx, delta, p0=[1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)   
+    print cubic_param2
     # quartic_param2,pcov = optimize.curve_fit(quartic, xx, delta, p0=[1,1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)   
     # #print quartic_param2
     # quintic_param2,pcov = optimize.curve_fit(quintic, xx, delta, p0=[1,1,1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)

examples/test_bh_sorted.c

@@ -42,7 +42,7 @@
 #define task_limit 1e8
 #define const_G 1    // 6.6738e-8
 #define dist_min 0.5 /* Used for legacy walk only */
-#define dist_cutoff_ratio 1.5
+#define dist_cutoff_ratio 1.2
 #define iact_pair_direct iact_pair_direct_sorted_multipole
 #define ICHECK -1
 #define NO_SANITY_CHECKS
@@ -1170,11 +1170,17 @@ float correction_coefs[6*4] =
     -0.04397042,  0.10367151, -0.06130978,  1.00754152,
     0.00536257, -0.06370355,  0.19954405,  0.81828777,
     -0.00677179, -0.01237365,  0.02703505,  0.99483005,
-    0.01082535, -0.0262587,   0.00342499,  1.0151183,
-    -0.00352327,  0.02770999, -0.02312047,  1.00282468
+    0.00701166, -0.0113404,  -0.0159732,   1.02330933,
+    0.00297062,  0.01928414, -0.01929454,  1.00190972
   };
 
 
+static inline float cubic(float x, float a, float b, float c, float d)
+{
+  return a*x*x*x + b*x*x + c*x + d;
+}
+
+
 /**
  * @brief Compute the best-fitting bias correction factor for the acceleration of a particle 
  *
@@ -1187,7 +1193,6 @@ static inline void gravity_bias_correction( float d, float axis[3], float corr[3
 
   int k, orientation;
   int is_zero[3];
-  float limit_exact = ( dist_cutoff_ratio - 1. );
 
   /* Start with no applied correction */
   corr[0] = corr[1] = corr[2] = 1.f;
@@ -1200,6 +1205,10 @@ static inline void gravity_bias_correction( float d, float axis[3], float corr[3
   }
 
   //message("%f %f", d, limit_exact);
+#if 1
+
+  float limit_exact = ( dist_cutoff_ratio - 1. );
+
 
   /* Apply the correction corresponding to this axis */
   switch(orientation)
@@ -1207,10 +1216,17 @@ static inline void gravity_bias_correction( float d, float axis[3], float corr[3
 
     case 1: /* side */
       {
-	if (d > limit_exact || d < -limit_exact)
+	if ( d > limit_exact)
 	  for( k = 0; k < 3; k++ )
-	    corr[k] = 0.;
+	    corr[k] = 1./cubic( d, correction_coefs[4*4+0], correction_coefs[4*4+1], correction_coefs[4*4+2], correction_coefs[4*4+3]);	
 	
+	else if ( d < -limit_exact)
+	  for( k = 0; k < 3; k++ )
+	    corr[k] = 1./cubic( d, correction_coefs[5*4+0], correction_coefs[5*4+1], correction_coefs[5*4+2], correction_coefs[5*4+3]);	
+
+	message("corr= %f %f %f", corr[0], corr[1], corr[2]);
+	
+
       }
       break;
 
@@ -1243,7 +1259,7 @@ static inline void gravity_bias_correction( float d, float axis[3], float corr[3
   for ( k = 0; k < 3; ++k )
     corr[k] = (is_zero[k] ? corr[k] : 1 );
 
-  
+#endif 
 }