Fixed force

include/force_simple.hpp

@@ -30,16 +30,13 @@ class ForceSimple: public Force<value_t, params_t, vars_t, buffer_t>{
 
                         value_t part1 = - drum_params.omega_sq * drum_vars.x_temp; //Hook restore
                         value_t part2 = - drum_params.c * drum_vars.xdot_temp; //damping
-                        value_t part3 =   drum_params.coulomb_prefactor * (drum_params.gapinv_sq - drum_params.gapinv_cb*drum_vars.x_temp) * drum_vars.V*drum_vars.V; //Coulomb force
-                        //interaction TODO: This is probably (slightly) incorrect for now. But it (most likely) has the same complexity as the real deal.
+                        value_t part3 =   drum_params.coulomb_prefactor * drum_vars.V*drum_vars.V * (1. + drum_vars.x_temp / drum_params.gap); //Coulomb drive
+                        //interaction DONE: This is probably (slightly) incorrect for now. But it (most likely) has the same complexity as the real deal.
                         value_t part4 = 0.;
-                        //TODO: this branching can be ommitted if using ni_vars.ti instead of drum_vars.ti. Let's time it.
-                        if(n0_vars.t0 != value_t(0.)) [[likely]]
-                                part4 += drum_vars.t0*drum_vars.t0*(drum_params.gapinv_sq - drum_params.gapinv_cb*(drum_params.gap - drum_vars.x_temp - n0_vars.x_temp)); //neighbour 0
-                        if(n0_vars.t1 != value_t(0.)) [[likely]]
-                                part4 += drum_vars.t1*drum_vars.t1*(drum_params.gapinv_sq - drum_params.gapinv_cb*(drum_params.gap - drum_vars.x_temp - n1_vars.x_temp)); //neighbour 1
-                        if(n0_vars.t2 != value_t(0.)) [[likely]]
-                                part4 += drum_vars.t2*drum_vars.t2*(drum_params.gapinv_sq - drum_params.gapinv_cb*(drum_params.gap - drum_vars.x_temp - n2_vars.x_temp)); //neighbour 2
+                        //DONE: branching can be ommitted if using ni_vars.ti instead of drum_vars.ti. Let's time it.
+                        part4 += n0_vars.t0*drum_vars.t0 * (1. + (drum_vars.x_temp + n0_vars.x_temp) / drum_params.gap);
+                        part4 += n1_vars.t0*drum_vars.t0 * (1. + (drum_vars.x_temp + n1_vars.x_temp) / drum_params.gap);
+                        part4 += n2_vars.t0*drum_vars.t0 * (1. + (drum_vars.x_temp + n2_vars.x_temp) / drum_params.gap);
                         part4 *= drum_params.coulomb_prefactor;
 
                         return part1 + part2 + part3 + part4;