Untitled

Vect calculateDirVec(const Vect &ionPos, double thetaDeg, double phiDeg)
{

  double theta = thetaDeg * M_PI / 180.0;
  double phi = phiDeg * M_PI / 180.0;

  Vect zLocal = normalize(ionPos);

  // To align phi = 0 towards global z-axis (0,0,1) we project it onto the plane perpendiuclar to zLocal
  Vect globalZ{0.0, 0.0, 1.0};
  Vect referenceX = globalZ - zLocal * dot(zLocal, globalZ);
  // Handle the degenerate case: if ionPos nearly aligns with globalZ (if proj is negligible), use a fallback.
  if (mag(referenceX) < 1e-5)
  {
    referenceX = Vect{1.0, 0.0, 0.0};
  }

  Vect xLocal = normalize(referenceX);
  Vect yLocal = normalize(cross(zLocal, xLocal)); // local RH coord system completed

  // Convert spherical coordinates to cartesian in the local coordinate system.
  double localX = std::sin(theta) * std::cos(phi);
  double localY = std::sin(theta) * std::sin(phi);
  double localZ = std::cos(theta);

  // Create the local direction vector
  Vect dirLocal{localX, localY, localZ};

  // Transform the local dir vector into the global coordinate system.
  Vect dirGlobal = xLocal * dirLocal.x + yLocal * dirLocal.y + zLocal * dirLocal.z;

  // Normalize
  return normalize(dirGlobal);
}