Fix/efficiency

Intern/rayx-core/src/Beamline/Objects/CircleSource.cpp

@@ -57,8 +57,7 @@ std::vector<Ray> CircleSource::getRays([[maybe_unused]] int thread_count) const
         // main ray (main ray: xDir=0,yDir=0,zDir=1 for phi=psi=0)
         glm::dvec3 direction = getDirection();
 
-        const auto rotation = glm::dmat3(m_orientation);
-        const auto field = rotation * stokesToElectricField(m_stokes);
+        const auto field = stokesToElectricField(m_stokes, direction);
 
         Ray r = {position, ETYPE_UNINIT, direction, en, field, 0.0, 0.0, -1.0, -1.0};
 

Intern/rayx-core/src/Beamline/Objects/DipoleSource.cpp

@@ -131,8 +131,7 @@ std::vector<Ray> DipoleSource::getRays(int thread_count) const {
         glm::dvec4 tempDir = m_orientation * glm::dvec4(direction, 0.0);
         direction = glm::dvec3(tempDir.x, tempDir.y, tempDir.z);
 
-        const auto rotation = glm::dmat3(m_orientation);
-        const auto field = rotation * stokesToElectricField(psiandstokes.stokes);
+        const auto field = stokesToElectricField(psiandstokes.stokes, direction);
 
         Ray r = {position, ETYPE_UNINIT, direction, en, field, 0.0, 0.0, -1.0, -1.0};
 #if defined(DIPOLE_OMP)

Intern/rayx-core/src/Beamline/Objects/MatrixSource.cpp

@@ -59,8 +59,7 @@ std::vector<Ray> MatrixSource::getRays([[maybe_unused]] int thread_count) const
             glm::dvec4 tempDir = m_orientation * glm::dvec4(direction, 0.0);
             direction = glm::dvec3(tempDir.x, tempDir.y, tempDir.z);
 
-            const auto rotation = glm::dmat3(m_orientation);
-            const auto field = rotation * stokesToElectricField(m_pol);
+            const auto field = stokesToElectricField(m_pol, direction);
 
             Ray r = {position, ETYPE_UNINIT, direction, en, field, 0.0, 0.0, -1.0, -1.0};
 

Intern/rayx-core/src/Beamline/Objects/PixelSource.cpp

@@ -74,7 +74,7 @@ std::vector<Ray> PixelSource::getRays([[maybe_unused]] int thread_count) const {
         direction = glm::dvec3(tempDir.x, tempDir.y, tempDir.z);
 
         const auto rotation = glm::dmat3(m_orientation);
-        const auto field = rotation * stokesToElectricField(m_pol);
+        const auto field = stokesToElectricField(m_pol, rotation);
 
         Ray r = {position, ETYPE_UNINIT, direction, en, field, 0.0, 0.0, -1.0, -1.0};
 

Intern/rayx-core/src/Beamline/Objects/PointSource.cpp

@@ -92,8 +92,7 @@ std::vector<Ray> PointSource::getRays(int thread_count) const {
         glm::dvec4 tempDir = m_orientation * glm::dvec4(direction, 0.0);
         direction = glm::dvec3(tempDir.x, tempDir.y, tempDir.z);
 
-        // const auto rotation = rotationMatrix(direction);
-        const auto field = /* rotation *  */ stokesToElectricField(m_pol);
+        const auto field = stokesToElectricField(m_pol, direction);
 
         Ray r = {position, ETYPE_UNINIT, direction, en, field, 0.0, 0.0, -1.0, -1.0};
 #if defined(DIPOLE_OMP)

Intern/rayx-core/src/Beamline/Objects/SimpleUndulatorSource.cpp

@@ -72,8 +72,7 @@ std::vector<Ray> SimpleUndulatorSource::getRays([[maybe_unused]] int thread_coun
         glm::dvec4 tempDir = m_orientation * glm::dvec4(direction, 0.0);
         direction = glm::dvec3(tempDir.x, tempDir.y, tempDir.z);
 
-        const auto rotation = glm::dmat3(m_orientation);
-        const auto field = rotation * stokesToElectricField(m_pol);
+        const auto field = stokesToElectricField(m_pol, direction);
 
         Ray r = {position, ETYPE_UNINIT, direction, en, field, 0.0, 0.0, -1.0, -1.0};
 

Intern/rayx-core/src/Shader/Efficiency.h

@@ -1,3 +1,9 @@
+/*
+ * This code implements several formulas from the following book:
+ * Russel A. Chipman, Wai-Sze Tiffany Lam, Garam Young "Polarized Light and Optical Systems" (2019).
+ * For specific references, see the inline comments in the respective functions.
+ */
+
 #pragma once
 
 #include <glm.h>
@@ -8,34 +14,43 @@
 
 namespace RAYX {
 
-using Stokes = glm::dvec4;
-using ElectricField = cvec3;
-using LocalElectricField = cvec2;
+using Stokes = glm::dvec4;         // Stokes vector represents the polarization state of light in the ray local plane
+using LocalElectricField = cvec2;  // Electric field in the ray local plane, as a complex 2D vector
+using ElectricField = cvec3;       // Electric field as a complex 3D vector
 
 struct FresnelCoeffs {
-    double s;
-    double p;
+    double s;  // Real component of Fresnel reflection/refraction coefficient for s-polarized light
+    double p;  // Real component of Fresnel reflection/refraction coefficient for p-polarized light
 };
 
 struct ComplexFresnelCoeffs {
-    complex::Complex s;
-    complex::Complex p;
+    complex::Complex s;  // Complex Fresnel coefficient for s-polarized light
+    complex::Complex p;  // Complex Fresnel coefficient for p-polarized light
 };
 
+// Computes the angle between two unit vectors in 3D space
 RAYX_FN_ACC
 inline double angleBetweenUnitVectors(glm::dvec3 a, glm::dvec3 b) { return glm::acos(glm::dot(a, b)); }
 
+// Calculates the refracted angle using Snell's Law, based on complex index of refraction
+// Source: "Polarized Light and Optical Systems", p. 297, Formula 8.3
 RAYX_FN_ACC
 inline complex::Complex calcRefractAngle(const complex::Complex incidentAngle, const complex::Complex iorI, const complex::Complex iorT) {
     return complex::asin((iorI / iorT) * complex::sin(incidentAngle));
 }
 
+// Calculates Brewster's angle, the angle at which no reflection occurs for p-polarized light
+// Source: "Polarized Light and Optical Systems", p. 305, Formula 8.26
 RAYX_FN_ACC
 inline complex::Complex calcBrewstersAngle(const complex::Complex iorI, const complex::Complex iorT) { return complex::atan(iorT / iorI); }
 
+// Calculates the critical angle for total internal reflection
+// Source: "Polarized Light and Optical Systems", p. 306, Formula 8.27
 RAYX_FN_ACC
 inline complex::Complex calcCriticalAngle(const complex::Complex iorI, const complex::Complex iorT) { return complex::asin(iorT / iorI); }
 
+// Computes the Fresnel reflection amplitude for s and p polarization states
+// Source: "Polarized Light and Optical Systems", p. 300, Formula 8.8 and 8.10
 RAYX_FN_ACC
 inline ComplexFresnelCoeffs calcReflectAmplitude(const complex::Complex incidentAngle, const complex::Complex refractAngle,
                                                  const complex::Complex iorI, const complex::Complex iorT) {
@@ -45,12 +60,11 @@ inline ComplexFresnelCoeffs calcReflectAmplitude(const complex::Complex incident
     const auto s = (iorI * cos_i - iorT * cos_t) / (iorI * cos_i + iorT * cos_t);
     const auto p = (iorT * cos_i - iorI * cos_t) / (iorT * cos_i + iorI * cos_t);
 
-    return {
-        .s = s,
-        .p = p,
-    };
+    return {.s = s, .p = p};
 }
 
+// Computes the Fresnel transmission (refracted) amplitude for s and p polarization states
+// Source: "Polarized Light and Optical Systems", p. 300, Formula 8.9 and 8.11
 RAYX_FN_ACC
 inline ComplexFresnelCoeffs calcRefractAmplitude(const complex::Complex incidentAngle, const complex::Complex refractAngle,
                                                  const complex::Complex iorI, const complex::Complex iorT) {
@@ -60,44 +74,39 @@ inline ComplexFresnelCoeffs calcRefractAmplitude(const complex::Complex incident
     const auto s = (2.0 * iorI * cos_i) / (iorI * cos_i + iorT * cos_t);
     const auto p = (2.0 * iorI * cos_i) / (iorT * cos_i + iorI * cos_t);
 
-    return {
-        .s = s,
-        .p = p,
-    };
+    return {.s = s, .p = p};
 }
 
+// Computes the reflectance (intensity) based on the Fresnel reflection amplitudes
+// Source: "Polarized Light and Optical Systems", p. 301, Formula 8.16 and 8.17
 RAYX_FN_ACC
 inline FresnelCoeffs calcReflectIntensity(const ComplexFresnelCoeffs reflectAmplitude) {
-    const auto s = (reflectAmplitude.s * complex::conj(reflectAmplitude.s)).real();
-    const auto p = (reflectAmplitude.p * complex::conj(reflectAmplitude.p)).real();
+    const auto s = (reflectAmplitude.s * complex::conj(reflectAmplitude.s)).real();  // s-polarized intensity
+    const auto p = (reflectAmplitude.p * complex::conj(reflectAmplitude.p)).real();  // p-polarized intensity
 
-    return {
-        .s = s,
-        .p = p,
-    };
+    return {.s = s, .p = p};
 }
 
+// Computes the transmittance (intensity) based on the Fresnel transmission amplitudes
+// Source: "Polarized Light and Optical Systems", p. 302, Formula 8.19 and 8.20
 RAYX_FN_ACC
 inline FresnelCoeffs calcRefractIntensity(const ComplexFresnelCoeffs refract_amplitude, const complex::Complex incidentAngle,
                                           const complex::Complex refractAngle, const complex::Complex iorI, const complex::Complex iorT) {
     const auto r = ((iorT * complex::cos(refractAngle)) / (iorI * complex::cos(incidentAngle))).real();
 
-    const auto s = r * (refract_amplitude.s * complex::conj(refract_amplitude.s)).real();
-    const auto p = r * (refract_amplitude.p * complex::conj(refract_amplitude.p)).real();
+    const auto s = r * (refract_amplitude.s * complex::conj(refract_amplitude.s)).real();  // s-polarized intensity
+    const auto p = r * (refract_amplitude.p * complex::conj(refract_amplitude.p)).real();  // p-polarized intensity
 
-    return {
-        .s = s,
-        .p = p,
-    };
+    return {.s = s, .p = p};
 }
 
+// Computes the Jones matrix, representing the effect on the polarization state, based on Fresnel coefficients
+// Source: "Polarized Light and Optical Systems", p. 385, Formula 10.26
 RAYX_FN_ACC
-inline cmat3 calcJonesMatrix(const ComplexFresnelCoeffs amplitude) {
-    return {
-        amplitude.s, 0, 0, 0, amplitude.p, 0, 0, 0, 1,
-    };
-}
+inline cmat3 calcJonesMatrix(const ComplexFresnelCoeffs amplitude) { return {amplitude.s, 0, 0, 0, amplitude.p, 0, 0, 0, 1}; }
 
+// Computes the polarization transformation matrix for reflection/refraction
+// Source: "Polarized Light and Optical Systems", p. 386, Formula 10.22 and 10.31
 RAYX_FN_ACC
 inline cmat3 calcPolaririzationMatrix(const glm::dvec3 incidentVec, const glm::dvec3 reflectOrRefractVec, const glm::dvec3 normalVec,
                                       const ComplexFresnelCoeffs amplitude) {
@@ -107,24 +116,22 @@ inline cmat3 calcPolaririzationMatrix(const glm::dvec3 incidentVec, const glm::d
     const auto p1 = glm::cross(reflectOrRefractVec, s0);
 
     const auto out = glm::dmat3(s1, p1, reflectOrRefractVec);
-
     const auto in = glm::dmat3(s0.x, p0.x, incidentVec.x, s0.y, p0.y, incidentVec.y, s0.z, p0.z, incidentVec.z);
 
     const auto jonesMatrix = calcJonesMatrix(amplitude);
-
     return out * jonesMatrix * in;
 }
 
+// Computes the polarization matrix for reflection at normal incidence (simplified case)
+// Source: "Polarized Light and Optical Systems", p. 387, Formula 10.35
 RAYX_FN_ACC
 inline cmat3 calcReflectPolarizationMatrixAtNormalIncidence(const ComplexFresnelCoeffs amplitude) {
-    // since no plane of incidence is defined at normal incidence,
-    // s and p components are equal and only contain the base reflectivity and a phase shift of 180 degrees
-    // here we apply the base reflectivity and phase shift independent of the ray direction to all components
     return {
-        amplitude.s, 0, 0, 0, amplitude.s, 0, 0, 0, amplitude.s,
+        amplitude.s, 0, 0, 0, amplitude.s, 0, 0, 0, amplitude.s,  // All components equal since it's normal incidence
     };
 }
 
+// Calculates the reflected electric field after intercepting a surface
 RAYX_FN_ACC
 inline ElectricField interceptReflect(const ElectricField incidentElectricField, const glm::dvec3 incidentVec, const glm::dvec3 reflectVec,
                                       const glm::dvec3 normalVec, const complex::Complex iorI, const complex::Complex iorT) {
@@ -133,45 +140,55 @@ inline ElectricField interceptReflect(const ElectricField incidentElectricField,
 
     const auto reflectAmplitude = calcReflectAmplitude(incidentAngle, refractAngle, iorI, iorT);
 
-    // TODO: make this more robust
-    const auto isNormalIncidence = incidentVec == -normalVec;
+    const auto isNormalIncidence = incidentVec == -normalVec;  // Check for normal incidence
     const auto reflectPolarizationMatrix = isNormalIncidence ? calcReflectPolarizationMatrixAtNormalIncidence(reflectAmplitude)
                                                              : calcPolaririzationMatrix(incidentVec, reflectVec, normalVec, reflectAmplitude);
 
-    const auto reflectElectricField = reflectPolarizationMatrix * incidentElectricField;
+    const auto reflectElectricField = reflectPolarizationMatrix * incidentElectricField;  // Reflect electric field
     return reflectElectricField;
 }
 
+// Computes the intensity from a 2D electric field (local)
+// Source: "Polarized Light and Optical Systems", p. 76, Formula 3.25
 RAYX_FN_ACC
 inline double intensity(const LocalElectricField field) {
     const auto mag = complex::abs(field);
     return glm::dot(mag, mag);
 }
 
+// Computes the intensity from a 3D electric field (global)
+// Source: "Polarized Light and Optical Systems", p. 76, Formula 3.25
 RAYX_FN_ACC
 inline double intensity(const ElectricField field) {
     const auto mag = complex::abs(field);
     return glm::dot(mag, mag);
 }
 
+// Returns the intensity from a Stokes vector
+// Source: https://en.wikipedia.org/wiki/Stokes_parameters
 RAYX_FN_ACC
 inline double intensity(const Stokes stokes) { return stokes.x; }
 
+// Computes the degree of polarization from a Stokes vector
+// Source: https://en.wikipedia.org/wiki/Stokes_parameters
 RAYX_FN_ACC
 inline double degreeOfPolarization(const Stokes stokes) { return glm::length(glm::vec3(stokes.y, stokes.z, stokes.w)) / stokes.x; }
 
+// Converts a local electric field to a Stokes vector
+// Source: "Polarized Light and Optical Systems", p. 76, Formula 3.25
 RAYX_FN_ACC
-inline Stokes fieldToStokes(const LocalElectricField field) {
+inline Stokes localElectricFieldToStokes(const LocalElectricField field) {
+    // note that we omit the magnitude scale with the impedance of free space (https://en.wikipedia.org/wiki/Impedance_of_free_space)
+
     const auto mag = complex::abs(field);
     const auto theta = complex::arg(field);
 
     return Stokes(mag.x * mag.x + mag.y * mag.y, mag.x * mag.x - mag.y * mag.y, 2.0 * mag.x * mag.y * glm::cos(theta.x - theta.y),
                   2.0 * mag.x * mag.y * glm::sin(theta.x - theta.y));
 }
 
-RAYX_FN_ACC
-inline Stokes fieldToStokes(const ElectricField field) { return fieldToStokes(LocalElectricField(field)); }
-
+// Converts a Stokes vector into a local electric field
+// Source: "Polarized Light and Optical Systems", p. 77, Formula 3.28
 RAYX_FN_ACC
 inline LocalElectricField stokesToLocalElectricField(const Stokes stokes) {
     const auto x_real = glm::sqrt((stokes.x + stokes.y) / 2.0);
@@ -183,30 +200,94 @@ inline LocalElectricField stokesToLocalElectricField(const Stokes stokes) {
     return LocalElectricField({x_real, 0}, y);
 }
 
-RAYX_FN_ACC
-inline ElectricField stokesToElectricField(const Stokes stokes) { return ElectricField(stokesToLocalElectricField(stokes), complex::Complex(0, 0)); }
+struct RotationBase {
+    glm::dvec3 right;
+    glm::dvec3 up;
+    glm::dvec3 forward;
+};
 
+// Computes a base given a forward vector
+// A convention for the up and right vectors is implemented, making this function well defined and for all forward directions
+// TODO(Sven): this convention should be exchanged with one that does not branch
 RAYX_FN_ACC
-inline glm::dmat3 rotationMatrix(const glm::dvec3 forward) {
+inline RotationBase forwardVectorToBaseConvention(glm::dvec3 forward) {
     auto up = glm::dvec3(0, 1, 0);
     glm::dvec3 right;
 
+    // If the forward vector is not nearly vertical, we use the up vector (0, 1, 0), othwewise use the right vector (1, 0, 0)
     if (glm::abs(glm::dot(forward, up)) < .5) {
-        right = glm::normalize(glm::cross(forward, up));
-        up = glm::normalize(glm::cross(right, forward));
+        right = glm::normalize(glm::cross(up, forward));
+        up = glm::cross(forward, right);
     } else {
         right = glm::dvec3(1, 0, 0);
         up = glm::normalize(glm::cross(forward, right));
-        right = glm::normalize(glm::cross(forward, up));
+        right = glm::cross(up, forward);
     }
 
-    return glm::dmat3(right, up, forward);
+    return RotationBase{
+        .right = right,
+        .up = up,
+        .forward = forward,
+    };
+}
+
+// Computes a rotation matrix given a forward vector. This matrix can be used to align an object with a direction
+RAYX_FN_ACC
+inline glm::dmat3 rotationMatrix(const glm::dvec3 forward) {
+    const auto base = forwardVectorToBaseConvention(forward);
+    return glm::dmat3(base.right, base.up, base.forward);
 }
 
+// Computes a rotation matrix given both forward and up vectors, used to align objects in space
 RAYX_FN_ACC
 inline glm::dmat3 rotationMatrix(const glm::dvec3 forward, const glm::dvec3 up) {
     const auto right = glm::cross(forward, up);
     return glm::dmat3(right, up, forward);
 }
 
+RAYX_FN_ACC
+inline ElectricField localToGlobalElectricField(LocalElectricField localField, glm::dvec3 forward) {
+    return rotationMatrix(forward) * ElectricField(localField, complex::Complex{0, 0});
+}
+
+RAYX_FN_ACC
+inline ElectricField localToGlobalElectricField(LocalElectricField localField, glm::dvec3 forward, glm::dvec3 up) {
+    return rotationMatrix(forward, up) * ElectricField(localField, complex::Complex{0, 0});
+}
+
+RAYX_FN_ACC
+inline LocalElectricField globalToLocalElectricField(ElectricField field, glm::dvec3 forward) {
+    return glm::transpose(rotationMatrix(forward)) * field;
+}
+
+RAYX_FN_ACC
+inline LocalElectricField globalToLocalElectricField(ElectricField field, glm::dvec3 forward, glm::vec3 up) {
+    return glm::transpose(rotationMatrix(forward, up)) * field;
+}
+
+RAYX_FN_ACC
+inline ElectricField stokesToElectricField(const Stokes stokes, const glm::dvec3 forward) {
+    return localToGlobalElectricField(stokesToLocalElectricField(stokes), forward);
+}
+
+RAYX_FN_ACC
+inline ElectricField stokesToElectricField(const Stokes stokes, const glm::dvec3 forward, const glm::dvec3 up) {
+    return localToGlobalElectricField(stokesToLocalElectricField(stokes), forward, up);
+}
+
+RAYX_FN_ACC
+inline ElectricField stokesToElectricField(const Stokes stokes, const glm::dmat3 rotation) {
+    return rotation * ElectricField(stokesToLocalElectricField(stokes), complex::Complex{0, 0});
+}
+
+RAYX_FN_ACC
+inline ElectricField electricFieldToStokes(const ElectricField field, const glm::dvec3 forward) {
+    return localElectricFieldToStokes(globalToLocalElectricField(field, forward));
+}
+
+RAYX_FN_ACC
+inline ElectricField electricFieldToStokes(const ElectricField field, const glm::dvec3 forward, const glm::dvec3 up) {
+    return localElectricFieldToStokes(globalToLocalElectricField(field, forward, up));
+}
+
 }  // namespace RAYX

Intern/rayx-core/src/Shader/Utils.h

@@ -16,6 +16,7 @@ RAYX_FN_ACC
 inline Ray RAYX_API rayMatrixMult(Ray r, const glm::dmat4 m) {
     r.m_position = glm::dvec3(m * glm::dvec4(r.m_position, 1));
     r.m_direction = glm::dvec3(m * glm::dvec4(r.m_direction, 0));
+    r.m_field = glm::dmat3(m) * r.m_field;
     return r;
 }