eic-opticks/qudarap/qwls.h

0001 #pragma once
0002 /**
0003 qwls.h : GPU-side Wavelength Shifting
0004 =========================================
0005
0006 Device-side struct for WLS wavelength sampling via ICDF texture lookup.
0007 Supports multiple WLS materials indexed by material ID.
0008
0009 The ICDF texture layout:
0010 - Each WLS material occupies 3 rows (standard, LHS HD, RHS HD)
0011 - material_map[mat_idx] gives the base row for that material (-1 = no WLS)
0012 - time_constants[wls_idx] gives the re-emission time constant in ns
0013
0014 Wavelength sampling uses the same HD (high-definition) technique as qscint.h:
0015 - hd_factor=20: 20x resolution at extremes (u < 0.05 or u > 0.95)
0016 - Normalized texture coordinates with linear interpolation
0017
0018 **/
0019
0020 #if defined(__CUDACC__) || defined(__CUDABE__)
0021 #define QWLS_METHOD __device__
0022 #else
0023 #define QWLS_METHOD
0024 #endif
0025
0026 struct qwls
0027 {
0028     cudaTextureObject_t wls_tex; // ICDF texture: (num_wls*3, 4096, 1)
0029     unsigned hd_factor;          // 0, 10, or 20
0030     int *material_map;           // device ptr: mat_idx -> base ICDF row (-1 = no WLS)
0031     float *time_constants;       // device ptr: per-WLS-material time constant (ns)
0032     unsigned num_wls;            // number of WLS materials
0033     unsigned tex_height;         // total rows in texture = num_wls * 3
0034
0035 #if defined(__CUDACC__) || defined(__CUDABE__) || defined(MOCK_CURAND) || defined(MOCK_CUDA)
0036
0037     QWLS_METHOD bool has_wls(unsigned mat_idx) const;
0038     QWLS_METHOD float wavelength(unsigned mat_idx, const float &u0) const;
0039     QWLS_METHOD float wavelength_at_row(unsigned base_row, const float &u0) const;
0040     QWLS_METHOD float time_constant(unsigned mat_idx) const;
0041
0042 #endif
0043 };
0044
0045 #if defined(__CUDACC__) || defined(__CUDABE__) || defined(MOCK_CURAND) || defined(MOCK_CUDA)
0046
0047 /**
0048 qwls::has_wls
0049 ---------------
0050
0051 Returns true if material at mat_idx has WLS properties.
0052 The material_map holds -1 for non-WLS materials.
0053
0054 **/
0055
0056 inline QWLS_METHOD bool qwls::has_wls(unsigned mat_idx) const
0057 {
0058     return material_map[mat_idx] >= 0;
0059 }
0060
0061 /**
0062 qwls::time_constant
0063 ---------------------
0064
0065 Returns the WLS re-emission time constant in ns for the given material.
0066 Returns 0.f if material has no WLS (instant re-emission / no delay).
0067
0068 **/
0069
0070 inline QWLS_METHOD float qwls::time_constant(unsigned mat_idx) const
0071 {
0072     int base_row = material_map[mat_idx];
0073     if (base_row < 0)
0074         return 0.f;
0075     unsigned wls_idx = base_row / 3;
0076     return time_constants[wls_idx];
0077 }
0078
0079 /**
0080 qwls::wavelength
0081 -------------------
0082
0083 Sample a re-emitted wavelength from the WLS emission spectrum ICDF
0084 for the material at mat_idx, using uniform random u0 in [0,1).
0085
0086 Returns 0.f if material has no WLS (should not happen in practice
0087 as callers check has_wls first).
0088
0089 **/
0090
0091 inline QWLS_METHOD float qwls::wavelength(unsigned mat_idx, const float &u0) const
0092 {
0093     int base_row = material_map[mat_idx];
0094     if (base_row < 0)
0095         return 0.f;
0096     return wavelength_at_row(unsigned(base_row), u0);
0097 }
0098
0099 /**
0100 qwls::wavelength_at_row
0101 --------------------------
0102
0103 ICDF texture lookup with HD (high-definition) support.
0104 base_row is the first of 3 rows for this WLS material:
0105   row 0: standard resolution  (full CDF range)
0106   row 1: LHS high-res         (0.00 -> 0.05 for hd_factor=20)
0107   row 2: RHS high-res         (0.95 -> 1.00 for hd_factor=20)
0108
0109 Uses normalized texture coordinates with linear interpolation,
0110 matching the qscint.h implementation.
0111
0112 **/
0113
0114 inline QWLS_METHOD float qwls::wavelength_at_row(unsigned base_row, const float &u0) const
0115 {
0116     float y0 = (float(base_row) + 0.5f) / float(tex_height);
0117     float y1 = (float(base_row + 1) + 0.5f) / float(tex_height);
0118     float y2 = (float(base_row + 2) + 0.5f) / float(tex_height);
0119
0120     float wl;
0121
0122     if (hd_factor == 0)
0123     {
0124         wl = tex2D<float>(wls_tex, u0, y0);
0125     }
0126     else if (hd_factor == 10)
0127     {
0128         if (u0 < 0.1f)
0129             wl = tex2D<float>(wls_tex, u0 * 10.f, y1);
0130         else if (u0 > 0.9f)
0131             wl = tex2D<float>(wls_tex, (u0 - 0.9f) * 10.f, y2);
0132         else
0133             wl = tex2D<float>(wls_tex, u0, y0);
0134     }
0135     else // hd_factor == 20
0136     {
0137         if (u0 < 0.05f)
0138             wl = tex2D<float>(wls_tex, u0 * 20.f, y1);
0139         else if (u0 > 0.95f)
0140             wl = tex2D<float>(wls_tex, (u0 - 0.95f) * 20.f, y2);
0141         else
0142             wl = tex2D<float>(wls_tex, u0, y0);
0143     }
0144
0145     return wl;
0146 }
0147
0148 #endif