prepare advance electric field according to path length

not actually advance, to ensure tests run at this point

Author: ichinii

Intern/rayx-core/src/Element/Behaviour.cpp

@@ -128,7 +128,7 @@ Behaviour makeRZPBehaviour(const DesignElement& dele) {
 
     // designEnergy = designEnergy; // if Auto == true, take energy of Source
     // (param sourceEnergy), else designEnergy = designEnergy
-    auto designWavelength = designEnergy == 0 ? 0 : hvlam(designEnergy);
+    auto designWavelength = designEnergy == 0 ? 0 : energyToWaveLength(designEnergy);
     auto additionalOrder = double(dele.getAdditionalOrder());
 
     auto imageType = dele.getImageType();
@@ -160,4 +160,4 @@ Behaviour makeFoil(const DesignElement& dele) {
     });
 }
 
-}  // namespace RAYX
+}  // namespace RAYX

Intern/rayx-core/src/Shader/Behave.cpp

@@ -28,7 +28,7 @@ Ray behaveCrystal(Ray r, const Behaviour behaviour, [[maybe_unused]] Collision c
     double polFactorS = 1.0;
     double polFactorP = std::fabs(cos(2 * bragg));
 
-    double wavelength = inm2eV / r.m_energy;
+    double wavelength = energyToWaveLength(r.m_energy);
     double gamma = getDiffractionPrefactor(wavelength, b.m_unitCellVolume);
 
     std::complex<double> F0(b.m_structureFactorReF0, b.m_structureFactorImF0);
@@ -80,7 +80,7 @@ Ray behaveSlit(Ray r, const Behaviour behaviour, Rand& __restrict rand) {
 
     double dPhi = 0;
     double dPsi = 0;
-    double wavelength = hvlam(r.m_energy);
+    double wavelength = energyToWaveLength(r.m_energy);
 
     // this was previously called "diffraction"
     if (wavelength > 0) {
@@ -109,7 +109,7 @@ RAYX_FN_ACC
 Ray behaveRZP(Ray r, const Behaviour behaviour, const Collision col, Rand& __restrict rand) {
     RZPBehaviour b = deserializeRZP(behaviour);
 
-    double WL = hvlam(r.m_energy);
+    double WL = energyToWaveLength(r.m_energy);
     double Ord = b.m_orderOfDiffraction;
     int additional_order = int(b.m_additionalOrder);
 
@@ -137,7 +137,7 @@ Ray behaveGrating(Ray r, const Behaviour behaviour, const Collision col) {
     GratingBehaviour b = deserializeGrating(behaviour);
 
     // vls parameters passed in q.elementParams
-    double WL = hvlam(r.m_energy);
+    double WL = energyToWaveLength(r.m_energy);
     double lineDensity = b.m_lineDensity;
     double orderOfDiffraction = b.m_orderOfDiffraction;
 
@@ -177,7 +177,7 @@ RAYX_FN_ACC
 Ray behaveFoil(Ray r, const Behaviour behaviour, const Collision col, const int material, const int* __restrict materialIndices,
                const double* __restrict materialTable) {
     FoilBehaviour f = deserializeFoil(behaviour);
-    const double wavelength = hvlam(r.m_energy);
+    const double wavelength = energyToWaveLength(r.m_energy);
 
     const auto indexVacuum = complex::Complex(1., 0.);
     const auto indexMaterial = getRefractiveIndex(r.m_energy, material, materialIndices, materialTable);

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

@@ -8,12 +8,9 @@
 
 namespace RAYX {
 
-// inverse nanometer-electron volt relationship / reziprocal factor
+/// inverse nanometer-electron volt relationship / reziprocal factor
 constexpr double INV_NM_TO_EVOLT = 1239.841984332002622;
 
-// inverse nanometer-electron volt relationship / reziprocal factor
-constexpr double inm2eV = 1.239852e3;  // TODO: Check, find better name maybe
-
 constexpr double PI = 3.14159265358979323846264338327950;
 // 141592653589793238462643383279502884197169399
 
@@ -39,7 +36,6 @@ constexpr double ELECTRON_MASS = 9.1093837015e-31;  // Checked 2019-7-25, PB, NI
 /// RAY-UI shortcut for getFactorElectronEnergy
 constexpr double FACTOR_ELECTRON_ENERGY_SC = 1957;
 
-//
 constexpr double FINE_STRUCTURE_CONSTANT = 7.2973525693e-3;  // Checked 2020-2-18, PB, NIST: fine-structure constant, no units
 
 /// \f$\epsilon_0\f$ [As/Vm]. vacuum electric permittivity.
@@ -53,6 +49,6 @@ constexpr double FACTOR_SCHWINGER_RAY = 1.2556937e15;
 constexpr double ELECTRIC_PERMITIVITY_MULTIPLIES_SPEED_OF_LIGHT = 2.6544187279929558624e-3;
 
 constexpr double ELECTRON_RADIUS = 2.8179403205e-6;  // Classical electron radius (nm)
-
 // constexpr double ELECTRON_RADIUS = (ELEMENTARY_CHARGE / (SPEED_OF_LIGHT * SPEED_OF_LIGHT)) * 1.e9;  // Classical electron radius in nm
+
 }  // namespace RAYX

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

@@ -59,7 +59,7 @@ inline double getTheta(Ray r, glm::dvec3 normal, double offsetAngle) {
 RAYX_FN_ACC
 inline double getBraggAngle(double energy, double dSpacing2) {
     int order = 1;
-    double wavelength = hvlam(energy);
+    double wavelength = energyToWaveLength(energy);
     double theta_factor = (order * wavelength) / dSpacing2;
 
     // Check for physical validity

Intern/rayx-core/src/Shader/DynamicElements.cpp

@@ -2,6 +2,7 @@
 
 #include "Behave.h"
 #include "Collision.h"
+#include "RecordEvent.h"
 #include "Utils.h"
 
 namespace RAYX {
@@ -64,7 +65,7 @@ void dynamicElements(const int gid, const InvState& inv, OutputEvents& outputEve
 
         // write ray in local element coordinates to global memory
         if (numRecorded < inv.maxEvents && (!inv.recordMask || inv.recordMask[col.elementIndex])) {
-            outputEvents.events[gid * inv.maxEvents + numRecorded] = ray;
+            recordEvent(outputEvents.events, ray, getRecordIndex(gid, numRecorded, inv.maxEvents));
             ++numRecorded;
         }
 
@@ -75,9 +76,9 @@ void dynamicElements(const int gid, const InvState& inv, OutputEvents& outputEve
     // check if the number of events exceeds capacity
     if (!colNotFound && inv.sequential == Sequential::No && isRayActive(ray.m_eventType)) {
         Collision col = findCollisionNonSequential(ray.m_position, ray.m_direction, inv.elements, inv.numElements, rand);
-        if (col.found) {
+        if (col.found && (!inv.recordMask || inv.recordMask[col.elementIndex])) {
             ray = terminateRay(ray, EventType::TooManyEvents);
-            outputEvents.events[gid * inv.maxEvents + inv.maxEvents - 1] = ray;
+            recordEvent(outputEvents.events, ray, getRecordIndex(gid, inv.maxEvents, inv.maxEvents - 1));
         }
     }
 

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

@@ -27,6 +27,35 @@ inline double intensity(const Stokes stokes) { return stokes.x; }
 RAYX_FN_ACC
 inline double degreeOfPolarization(const Stokes stokes) { return glm::length(glm::vec3(stokes.y, stokes.z, stokes.w)) / stokes.x; }
 
+/**
+ * Advances an electric field propagating in a certain media by a certain distance
+ * @param field electric field of incident photon
+ * @param waveLength wavelength of incident photon
+ * @param distance distance traveled in media
+ * @param ior_real real component of complex index of refraction of media
+ * @return advanced electric field
+ * @TODO: check if waveLength and distance use the same
+ */
+RAYX_FN_ACC
+inline ElectricField advanceElectricField(const ElectricField field, double waveLength, const double distance, const float ior_real) {
+    // bring wavelength from nanometers into millimeters
+    waveLength /= 1e6;
+
+    // compute wave number (2π * n / λ)
+    const double waveNumber = 2.0 * PI * ior_real / waveLength;
+
+    // reduce the distance modulo wavelength to avoid large angle errors
+    const double reducedDistance = std::fmod(distance, waveLength);
+
+    // compute the phase shift using the reduced distance
+    const double deltaPhi = waveNumber * reducedDistance;
+
+    // apply the complex exponential of the phase shift
+    const auto phaseShift = complex::exp(complex::Complex(0.0, deltaPhi));
+
+    return field * phaseShift;
+}
+
 /*
  *  rotation of electric field
  */

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

@@ -40,6 +40,8 @@ struct RAYX_API Ray {
     /// The complex electric field
     ElectricField m_field;
 
+    // TODO: until we support transmission in media with refractive index != 1.0, path length equivalent to optical path length. in future, we should
+    // consider renaming this attribute to optical path length
     /// The distance that this ray has already traveled (in mm).
     double m_pathLength;
 

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

@@ -0,0 +1,23 @@
+#pragma once
+
+#include "Ray.h"
+#include "Utils.h"
+
+namespace RAYX {
+
+RAYX_FN_ACC
+inline int getRecordIndex(const int pathId, const int bounceId, const int maxEvents) { return pathId * maxEvents + bounceId; }
+
+RAYX_FN_ACC
+inline void recordEvent(Ray* __restrict output, Ray& __restrict ray, const int recordIndex) {
+    // TODO: this code is currently disabled to ensure all tests run with the changes that were made. in the next commit this code will be enabled and
+    // seeded will be override, due to changes in the electric field of written rays
+
+    // we use a temporary here, to avoid copying the whole ray when writing the ray with the correct field value
+    // const auto tmpField = ray.m_field;
+    // ray.m_field = advanceElectricField(ray.m_field, energyToWaveLength(ray.m_energy), ray.m_pathLength, 1.0);
+    output[recordIndex] = ray;
+    // ray.m_field = tmpField;
+}
+
+}  // namespace RAYX

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

@@ -4,15 +4,14 @@
 
 namespace RAYX {
 
-// converts energy in eV to wavelength in nm
 RAYX_FN_ACC
-double RAYX_API hvlam(double x) {
-    if (x == 0) {
-        return 0.0;
-    }
+double RAYX_API energyToWaveLength(double x) {
     return INV_NM_TO_EVOLT / x;
 }
 
+RAYX_FN_ACC
+double waveLengthToEnergy(const double waveLength) { return INV_NM_TO_EVOLT / waveLength; }
+
 // originally we calculcated 1/0, which might be UB (https://stackoverflow.com/questions/5802351/what-happens-when-you-divide-by-0-in-a-shader).
 // `return 1e+308 * 2.0;` correctly returns positive infinity on my machine (Rudi), but we have no guarantee that it always does.
 // So instead we return a sufficiently large number to be used like positive infinity.

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

@@ -5,8 +5,19 @@
 
 namespace RAYX {
 
-// converts energy in eV to wavelength in nm
-RAYX_FN_ACC double RAYX_API hvlam(double x);
+/**
+ * converts energy (eV) to wavelength (nm)
+ * @param energy energy of a photon in eV. must not be 0
+ * @return wavelength of photon in nm
+ */
+RAYX_FN_ACC double RAYX_API energyToWaveLength(double energy);
+
+/**
+ * Convert photon wavelength (nm) to energy (eV)
+ * @param wavelength of photon in nm. must not be 0
+ * @return energy energy of a photon in eV
+ */
+RAYX_FN_ACC double RAYX_API waveLengthToEnergy(const double waveLength);
 
 RAYX_FN_ACC double infinity();