Implement the higher-order gravity softened derivatives

src/kernel_gravity.h

@@ -168,8 +168,8 @@ __attribute__((always_inline)) INLINE static void kernel_grav_eval_force_double(
 /* Derivatives of softening kernel used for FMM */
 /************************************************/
 
-__attribute__((always_inline)) INLINE static float D_soft_1(float u,
-                                                            float u_inv) {
+__attribute__((const)) INLINE static float D_soft_1(const float u,
+                                                    const float u_inv) {
 
   /* phi(u) = 3u^7 - 15u^6 + 28u^5 - 21u^4 + 7u^2 - 3 */
   float phi = 3.f * u - 15.f;
@@ -183,8 +183,8 @@ __attribute__((always_inline)) INLINE static float D_soft_1(float u,
   return phi;
 }
 
-__attribute__((always_inline)) INLINE static float D_soft_3(float u,
-                                                            float u_inv) {
+__attribute__((const)) INLINE static float D_soft_3(const float u,
+                                                    const float u_inv) {
 
   /* phi'(u)/u = 21u^5 - 90u^4 + 140u^3 - 84u^2 + 14 */
   float phi = 21.f * u - 90.f;
@@ -196,8 +196,8 @@ __attribute__((always_inline)) INLINE static float D_soft_3(float u,
   return phi;
 }
 
-__attribute__((always_inline)) INLINE static float D_soft_5(float u,
-                                                            float u_inv) {
+__attribute__((const)) INLINE static float D_soft_5(const float u,
+                                                    const float u_inv) {
 
   /* (phi'(u)/u)'/u = 105u^3 - 360u^2 + 420u - 168 */
   float phi = 105.f * u - 360.f;
@@ -207,18 +207,31 @@ __attribute__((always_inline)) INLINE static float D_soft_5(float u,
   return phi;
 }
 
-__attribute__((always_inline)) INLINE static float D_soft_7(float u,
-                                                            float u_inv) {
-  return 0.f;
+__attribute__((const)) INLINE static float D_soft_7(const float u,
+                                                    const float u_inv) {
+  return 315.f * u - 720.f + 420.f * u_inv;
 }
 
-__attribute__((always_inline)) INLINE static float D_soft_9(float u,
-                                                            float u_inv) {
-  return 0.f;
+__attribute__((const)) INLINE static float D_soft_9(const float u,
+                                                    const float u_inv) {
+  /* -315 u^-1 + 420 u^-3 */
+  /* -315 u_inv + 420 u_inv^3 */
+  float phi = 420.f * u_inv;
+  phi = phi * u_inv - 315.f;
+  phi = phi * u_inv;
+  return phi;
 }
 
-__attribute__((always_inline)) INLINE static float D_soft_11(float u,
-                                                             float u_inv) {
+__attribute__((const)) INLINE static float D_soft_11(const float u,
+                                                     const float u_inv) {
+  /* 315 u^-3 - 1260 u^-5 */
+  /* 315 u_inv^3 - 1260 u_inv^5 */
+  float phi = -1260.f * u_inv;
+  phi = phi * u_inv + 315.f;
+  phi = phi * u_inv;
+  phi = phi * u_inv;
+  phi = phi * u_inv;
+
   return 0.f;
 }