Untitled

// Parameters:
//   bar           : neigboring high charge atom
//   ionPos        : original ionPosition or coordinates
//   normal        : unit normal vector
//   firstNeighDist: nearest neigbor
//   thetaDeg      : the polar (tilt) angle in degrees.
//   phiDeg        : the azimuth (roll) angle in degrees.
Vect rollUpSpherical(const Vect &ionPos, const Vect &bar, const Vect &normal,
                     double firstNeighDist, double thetaDeg, double phiDeg)
{
  // Vect normal = normalize(ionPos);

  // Rolled-up ion position
  Vect ionPos_r = bar + normal * firstNeighDist;

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

  // Its magnitude
  double r = magnitude(dirVec);

  // Defining the new z'-axis as normalized direction vector
  Vect zPrime = normalize(dirVec);

  // Construct a local orthonormal frame (x', y', z') where:
  //    - zPrime is as defined above.
  //    - xPrime is chosen by crossing an arbitrary reference vector with zPrime.
  Vect reference(1.0, 0.0, 0.0);
  if (std::fabs(dot(zPrime, reference)) > 0.99)
  {
    reference = Vect(0.0, 1.0, 0.0);
  }
  Vect xPrime = normalize(cross(reference, zPrime));
  Vect yPrime = normalize(cross(zPrime, xPrime));

  double theta = toRadians(thetaDeg);
  double phi = toRadians(phiDeg);

  // 6) Calculatecomponents of the rotated direction vector (in the local frame).
  //   - Component along zPrime: r * cos(theta)
  //   - Component in the perpendicular (x'y' plane): r * sin(theta), further split using phi.
  double compX = r * sin(theta) * cos(phi);
  double compY = r * sin(theta) * sin(phi);
  double compZ = r * cos(theta);

  // 7) Transform the local vector components to the global coordinate system.
  Vect rotatedDirVec{
      compX * xPrime.x + compY * yPrime.x + compZ * zPrime.x,
      compX * xPrime.y + compY * yPrime.y + compZ * zPrime.y,
      compX * xPrime.z + compY * yPrime.z + compZ * zPrime.z};

  // 8) Add rotatedDirVec to the ionPos.
  Vect perturbedIonPos = ionPos + rotatedDirVec;

  return perturbedIonPos;
}