| 45 | 
  | 
    return; | 
| 46 | 
  | 
  } | 
| 47 | 
  | 
   | 
| 48 | 
– | 
   | 
| 48 | 
  | 
  vect r_ab; /*the vector whose origin is a and end is b */ | 
| 49 | 
  | 
  double force; /* the force scaling factor. */ | 
| 50 | 
  | 
  double Fab_x; /*the x,y, and z components of the force */ | 
| 51 | 
  | 
  double Fab_y; | 
| 52 | 
  | 
  double Fab_z; | 
| 53 | 
  | 
 | 
| 54 | 
+ | 
  double aR[3], bR[3]; | 
| 55 | 
+ | 
  double aF[3], bF[3]; | 
| 56 | 
+ | 
 | 
| 57 | 
  | 
  /* initialize the vector */ | 
| 58 | 
+ | 
   | 
| 59 | 
+ | 
  c_p_a->getPos(aR); | 
| 60 | 
+ | 
  c_p_b->getPos(bR); | 
| 61 | 
  | 
 | 
| 62 | 
< | 
  r_ab.x = c_p_b->getX() - c_p_a->getX(); | 
| 63 | 
< | 
  r_ab.y = c_p_b->getY() - c_p_a->getY(); | 
| 64 | 
< | 
  r_ab.z = c_p_b->getZ() - c_p_a->getZ(); | 
| 62 | 
> | 
  r_ab.x = bR[0] - aR[0]; | 
| 63 | 
> | 
  r_ab.y = bR[1] - aR[1]; | 
| 64 | 
> | 
  r_ab.z = bR[2] - aR[2]; | 
| 65 | 
> | 
 | 
| 66 | 
  | 
  r_ab.length = sqrt((r_ab.x * r_ab.x + r_ab.y * r_ab.y + r_ab.z * r_ab.z)); | 
| 67 | 
  | 
   | 
| 68 | 
  | 
  /* calculate the force here */ | 
| 73 | 
  | 
  Fab_y = -force *  r_ab.y / r_ab.length; | 
| 74 | 
  | 
  Fab_z = -force *  r_ab.z / r_ab.length; | 
| 75 | 
  | 
 | 
| 76 | 
< | 
  c_p_a->addFx(Fab_x); | 
| 77 | 
< | 
  c_p_a->addFy(Fab_y); | 
| 78 | 
< | 
  c_p_a->addFz(Fab_z); | 
| 76 | 
> | 
  aF[0] = Fab_x; | 
| 77 | 
> | 
  aF[1] = Fab_y; | 
| 78 | 
> | 
  aF[2] = Fab_z; | 
| 79 | 
  | 
 | 
| 80 | 
< | 
  c_p_b->addFx(-Fab_x); | 
| 81 | 
< | 
  c_p_b->addFy(-Fab_y); | 
| 82 | 
< | 
  c_p_b->addFz(-Fab_z); | 
| 80 | 
> | 
  bF[0] = -Fab_x; | 
| 81 | 
> | 
  bF[1] = -Fab_y; | 
| 82 | 
> | 
  bF[2] = -Fab_z; | 
| 83 | 
  | 
 | 
| 84 | 
+ | 
  c_p_a->addFrc(aF); | 
| 85 | 
+ | 
  c_p_b->addFrc(bF); | 
| 86 | 
+ | 
 | 
| 87 | 
  | 
  return; | 
| 88 | 
  | 
} |