Untitled

 avatar
unknown
c_cpp
15 days ago
2.1 kB
9
Indexable
// Parameters:
//   bar           : neigboring high charge atom
//   ionPos        : original ionPosition or cordinates
//   normal        : unit normal vector
//   firstNeighDist: nearest neigbour
//   thetaDeg      : the polar angle in degrees.
//   phiDeg        : the azimuth angle in degrees.
Vect rollUpSpherical(const Vect &ionPos, const Vect &bar, const Vect &normal,
                     double firstNeighDist, double thetaDeg, double phiDeg)
{

  // 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_local-axis as normalized direction vector
  Vect zLocal = normalize(dirVec);

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

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

  // 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);

  // Transform the local vector components in local coord sys to the global coordinate system.
  Vect rotatedDirVec{
      compX * xLocal.x + compY * yLocal.x + compZ * zLocal.x,
      compX * xLocal.y + compY * yLocal.y + compZ * zLocal.y,
      compX * xLocal.z + compY * yLocal.z + compZ * zLocal.z};

  // Add rotatedDirVec to the original ionPos.
  Vect rotatedIonPos = ionPos + rotatedDirVec;

  return rotatedIonPos;
}
Editor is loading...
Leave a Comment