Tero's suggestion for direction

// Parameters:

//   ionPos        : original ionPosition or cordinates
//   bar           : neigboring high charge atom
//   normal        : unit normal vector
//   firstNeighDist: nearest neigbour
//   thetaDeg      : the polar angle in degrees.
//   phiDeg        : the azimuth angle in degrees.

Vect rollUpVelocityDirection(const Vect &ionPos, const Vect &bar, const Vect &normal,
                     double firstNeighDist)
{
  // Rolled-up ion position
  Vect ionPos_r = bar + normal * firstNeighDist;

  // Direction Vector from original ionPos to rolled-up ionPos
  Vect dirVec = ionPos_r - ionPos;

  // Component of dirVec perpendicular to the normal (i.e., parallel to surface)
  Vect dirVec_tangent = dirVec - normal * dot(dirVec, normal);

  return normalize(dirVec_tangent);
}