diff --git a/examples/plot_sorted.py b/examples/plot_sorted.py
index 2caf184b597627575784e4bbfdf9a9101a9681a9..2ba84411af2bb0ff06c7d7030fda0165bc69b5aa 100644
--- a/examples/plot_sorted.py
+++ b/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)
diff --git a/examples/test_bh_sorted.c b/examples/test_bh_sorted.c
index 2345aa11f22bfea1a807f62f02777aeac67bc6de..ab5319a00a2d4f994759cd88b96a981b639b88fb 100644
--- a/examples/test_bh_sorted.c
+++ b/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
}