
package attracsimple;

public class AttracSimple {
   public static void main(String[] args) {
      //***************************************//
      // This parameters can be changed:       //
      double e = 8.0, s = 3.0; // s > 0        //
      //***************************************//
      double radi, dr, x;
      //Initial conditions (x, y, z)
      double x1 = 0.01, y1 = -0.02, z1 = 0.04;
      double x2 = -0.04, y2 = -0.05, z2 = 0.01;
      double x3 = 0.05, y3 = -0.04, z3 = -0.01;
      double rq12, r12, scalar12, dx12 = 1.0, dy12 = -1.0, dz12 = 1.0;
      double rq13, r13, scalar13, dx13 = -1.0, dy13 = 1.0, dz13 = -1.0;
      double eist = 1.0;
      int N = 15000; //steps
      //Iteration loop----------------------------------------------------------
      for(int i = 0; i < N; i++) {
         rq12 = (x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1);
         r12 = Math.sqrt(rq12);
         scalar12 = dx12*(x2-x1)+dy12*(y2-y1)+dz12*(z2-z1);
         if(eist > e && r12 > s) {
            radi = (rq12/(s*s)-1.0)*(1.0-(scalar12*scalar12)/rq12);
            if(radi < 0.0) radi= 0.0;
            x = -scalar12/r12+Math.sqrt(radi);
            dx12 = dx12-x*(x1-x2)/r12;
            dy12 = dy12-x*(y1-y2)/r12;
            dz12 = dz12-x*(z1-z2)/r12;
            dr = Math.sqrt(dx12*dx12+dy12*dy12+dz12*dz12);
            dx12 = dx12/dr;
            dy12 = dy12/dr;
            dz12 = dz12/dr;
         }
         rq13 = (x3-x1)*(x3-x1)+(y3-y1)*(y3-y1)+(z3-z1)*(z3-z1);
         r13 = Math.sqrt(rq13);
         scalar13 = dx13*(x3-x1)+dy13*(y3-y1)+dz13*(z3-z1);
         if(eist > e && r13 > s) {
            radi = (rq13/(s*s)-1.0)*(1.0-(scalar13*scalar13)/rq13);
            if(radi < 0.0) radi = 0.0;
            x = -scalar13/r13+Math.sqrt(radi);
            dx13 = dx13-x*(x1-x3)/r13;
            dy13 = dy13-x*(y1-y3)/r13;
            dz13 = dz13-x*(z1-z3)/r13;
            dr = Math.sqrt(dx13*dx13+dy13*dy13+dz13*dz13);
            dx13 = dx13/dr;
            dy13 = dy13/dr;
            dz13 = dz13/dr;
         }
      //Boundary value
      eist = r12+r13;
      //Superposition
      x1 = x1+dx12+dx13;
      y1 = y1+dy12+dy13;
      z1 = z1+dz12+dz13;
      x2 = x2-dx12;
      y2 = y2-dy12;
      z2 = z2-dz12;
      x3 = x3-dx13;
      y3 = y3-dy13;
      z3 = z3-dz13;
      //Position of node1, node2 and node3
      System.out.format("%10f %10f %10f %10f %10f %10f %n",
      x1, y1, x2, y2, x3, y3);
      }
   }
}