summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--chroma/cuda/geometry_types.h20
-rw-r--r--chroma/cuda/photon.h183
-rw-r--r--chroma/geometry.py4
3 files changed, 155 insertions, 52 deletions
diff --git a/chroma/cuda/geometry_types.h b/chroma/cuda/geometry_types.h
index 8836d91..5985be0 100644
--- a/chroma/cuda/geometry_types.h
+++ b/chroma/cuda/geometry_types.h
@@ -11,13 +11,33 @@ struct Material
float wavelength0;
};
+// surface models
+enum {
+ SURFACE_DEFAULT,
+ SURFACE_SPECULAR, // perfect specular reflector
+ SURFACE_DIFFUSE, // perfect diffuse reflector
+ SURFACE_COMBO, // both specular and diffuse components
+ SURFACE_MIRROR, // mirror including complex index
+ SURFACE_PHOTOCATHODE, // transmissive photocathode with complex index
+ SURFACE_TPB
+};
+
+// not all parameters are used by all surface models!
struct Surface
{
+ // process probabilities
float *detect;
float *absorb;
+ float *reemit;
+ float *reflect_tpb;
float *reflect_diffuse;
float *reflect_specular;
+
+ float *reemission_wavelength;
+ float *reemission_cdf;
+ unsigned int model;
unsigned int n;
+ unsigned int reemission_n;
float step;
float wavelength0;
};
diff --git a/chroma/cuda/photon.h b/chroma/cuda/photon.h
index 5677fc1..2e53621 100644
--- a/chroma/cuda/photon.h
+++ b/chroma/cuda/photon.h
@@ -50,6 +50,7 @@ enum
RAYLEIGH_SCATTER = 0x1 << 4,
REFLECT_DIFFUSE = 0x1 << 5,
REFLECT_SPECULAR = 0x1 << 6,
+ SURFACE_REEMIT = 0x1 << 7,
NAN_ABORT = 0x1 << 31
}; // processes
@@ -179,8 +180,9 @@ rayleigh_scatter(Photon &p, curandState &rng)
p.polarization /= norm(p.polarization);
} // scatter
-__device__ int propagate_to_boundary(Photon &p, State &s, curandState &rng,
- bool use_weights=false, int scatter_first=0)
+__device__
+int propagate_to_boundary(Photon &p, State &s, curandState &rng,
+ bool use_weights=false, int scatter_first=0)
{
float absorption_distance = -s.absorption_length*logf(curand_uniform(&rng));
float scattering_distance = -s.scattering_length*logf(curand_uniform(&rng));
@@ -317,77 +319,154 @@ propagate_at_boundary(Photon &p, State &s, curandState &rng)
} // propagate_at_boundary
__device__ int
-propagate_at_surface(Photon &p, State &s, curandState &rng, Geometry *geometry,
- bool use_weights=false)
+propagate_at_specular_reflector(Photon &p, State &s)
{
- Surface *surface = geometry->surfaces[s.surface_index];
+ float incident_angle = get_theta(s.surface_normal, -p.direction);
+ float3 incident_plane_normal = cross(p.direction, s.surface_normal);
+ incident_plane_normal /= norm(incident_plane_normal);
+
+ p.direction = rotate(s.surface_normal, incident_angle, incident_plane_normal);
+
+ p.history |= REFLECT_SPECULAR;
+
+ return CONTINUE;
+} // reflect_specular
+
+__device__ int
+propagate_at_diffuse_reflector(Photon &p, State &s, curandState &rng)
+{
+ p.direction = uniform_sphere(&rng);
+
+ if (dot(p.direction, s.surface_normal) < 0.0f)
+ p.direction = -p.direction;
+
+ p.polarization = cross(uniform_sphere(&rng), p.direction);
+ p.polarization /= norm(p.polarization);
+
+ p.history |= REFLECT_DIFFUSE;
+
+ return CONTINUE;
+} // reflect_diffuse
+
+__device__ int
+propagate_at_mirror(Photon &p, State &s, bool use_weights=false)
+{
+ // just a perfect specular reflector for now...
+ return propagate_at_specular_reflector(p, s);
+} // propagate_at_mirror
- /* since the surface properties are interpolated linearly, we are
- guaranteed that they still sum to 1.0. */
+__device__ int
+propagate_at_photocathode(Photon &p, State &s, curandState &rng, bool use_weights=false)
+{
+ return CONTINUE;
+} // propagate_at_photocathode
- float detect = interp_property(surface, p.wavelength, surface->detect);
+__device__ int
+propagate_at_tpb(Photon &p, State &s, curandState &rng, Surface *surface, bool use_weights=false)
+{
float absorb = interp_property(surface, p.wavelength, surface->absorb);
- float reflect_diffuse = interp_property(surface, p.wavelength, surface->reflect_diffuse);
- float reflect_specular = interp_property(surface, p.wavelength, surface->reflect_specular);
+ float reflect = interp_property(surface, p.wavelength, surface->reflect_tpb);
+ float reemit = interp_property(surface, p.wavelength, surface->reemit);
float uniform_sample = curand_uniform(&rng);
- if (use_weights && p.weight > WEIGHT_LOWER_THRESHOLD
- && absorb < (1.0f - WEIGHT_LOWER_THRESHOLD)) {
- // Prevent absorption and reweight accordingly
- float survive = 1.0f - absorb;
- absorb = 0.0f;
- p.weight *= survive;
-
- // Renormalize remaining probabilities
- detect /= survive;
- reflect_diffuse /= survive;
- reflect_specular /= survive;
- }
+ if (use_weights && p.weight > WEIGHT_LOWER_THRESHOLD && absorb < (1.0f - WEIGHT_LOWER_THRESHOLD)) {
+ // Prevent absorption and reweight accordingly
+ float survive = 1.0f - absorb;
+ absorb = 0.0f;
+ p.weight *= survive;
- if (use_weights && detect > 0.0f) {
- p.history |= SURFACE_DETECT;
- p.weight *= detect;
- return BREAK;
+ reflect /= survive;
+ reemit /= survive;
}
if (uniform_sample < absorb) {
- p.history |= SURFACE_ABSORB;
- return BREAK;
+ p.history |= SURFACE_ABSORB;
+ return BREAK;
}
- else if (uniform_sample < absorb + detect) {
- p.history |= SURFACE_DETECT;
- return BREAK;
- }
- else if (uniform_sample < absorb + detect + reflect_diffuse) {
- // diffusely reflect
- p.direction = uniform_sphere(&rng);
-
- if (dot(p.direction, s.surface_normal) < 0.0f)
- p.direction = -p.direction;
-
- // randomize polarization?
- p.polarization = cross(uniform_sphere(&rng), p.direction);
- p.polarization /= norm(p.polarization);
+ else if (uniform_sample < absorb + reemit) {
+ p.history |= SURFACE_ABSORB;
+ p.history |= SURFACE_REEMIT;
- p.history |= REFLECT_DIFFUSE;
+ // pick a new wavelength according to the reemission cdf
+ p.wavelength = sample_cdf(&rng, surface->reemission_n, surface->reemission_wavelength, surface->reemission_cdf);
- return CONTINUE;
+ // reemit isotropically (eh?)
+ return propagate_at_diffuse_reflector(p, s, rng);
}
else {
- // specularly reflect
- float incident_angle = get_theta(s.surface_normal,-p.direction);
- float3 incident_plane_normal = cross(p.direction, s.surface_normal);
- incident_plane_normal /= norm(incident_plane_normal);
+ return propagate_at_diffuse_reflector(p, s, rng);
+ }
- p.direction = rotate(s.surface_normal, incident_angle,
- incident_plane_normal);
+} // propagate_at_tpb
- p.history |= REFLECT_SPECULAR;
+__device__ int
+propagate_at_surface(Photon &p, State &s, curandState &rng, Geometry *geometry,
+ bool use_weights=false)
+{
+ Surface *surface = geometry->surfaces[s.surface_index];
- return CONTINUE;
+ if (surface->model == SURFACE_SPECULAR)
+ return propagate_at_specular_reflector(p, s);
+ else if (surface->model == SURFACE_DIFFUSE)
+ return propagate_at_diffuse_reflector(p, s, rng);
+ else if (surface->model == SURFACE_MIRROR)
+ return propagate_at_mirror(p, s);
+ else if (surface->model = SURFACE_PHOTOCATHODE) {
+ //propagate_at_photocathode();
+ p.history |= SURFACE_DETECT;
+ return BREAK;
+ }
+ else if (surface->model == SURFACE_TPB)
+ return propagate_at_tpb(p, s, rng, surface, use_weights);
+ else {
+ // if no surface model is specified, do a combination of specular and
+ // diffuse reflection, detection, and absorption based on relative
+ // probabilties (i.e. the "old" behavior)
+
+ // since the surface properties are interpolated linearly, we are
+ // guaranteed that they still sum to 1.0.
+ float detect = interp_property(surface, p.wavelength, surface->detect);
+ float absorb = interp_property(surface, p.wavelength, surface->absorb);
+ float reflect_diffuse = interp_property(surface, p.wavelength, surface->reflect_diffuse);
+ float reflect_specular = interp_property(surface, p.wavelength, surface->reflect_specular);
+
+ float uniform_sample = curand_uniform(&rng);
+
+ if (use_weights && p.weight > WEIGHT_LOWER_THRESHOLD
+ && absorb < (1.0f - WEIGHT_LOWER_THRESHOLD)) {
+ // Prevent absorption and reweight accordingly
+ float survive = 1.0f - absorb;
+ absorb = 0.0f;
+ p.weight *= survive;
+
+ // Renormalize remaining probabilities
+ detect /= survive;
+ reflect_diffuse /= survive;
+ reflect_specular /= survive;
+ }
+
+ if (use_weights && detect > 0.0f) {
+ p.history |= SURFACE_DETECT;
+ p.weight *= detect;
+ return BREAK;
+ }
+
+ if (uniform_sample < absorb) {
+ p.history |= SURFACE_ABSORB;
+ return BREAK;
+ }
+ else if (uniform_sample < absorb + detect) {
+ p.history |= SURFACE_DETECT;
+ return BREAK;
+ }
+ else if (uniform_sample < absorb + detect + reflect_diffuse)
+ return propagate_at_diffuse_reflector(p, s, rng);
+ else
+ return propagate_at_specular_reflector(p, s);
}
} // propagate_at_surface
#endif
+
diff --git a/chroma/geometry.py b/chroma/geometry.py
index 41b24a9..b40cda4 100644
--- a/chroma/geometry.py
+++ b/chroma/geometry.py
@@ -173,8 +173,12 @@ class Surface(object):
self.set('detect', 0)
self.set('absorb', 0)
+ self.set('reemit', 0)
+ self.set('reflect_tpb', 0)
self.set('reflect_diffuse', 0)
self.set('reflect_specular', 0)
+ self.set('reemission_wavelength', 0)
+ self.set('reemission_cdf', 0)
def set(self, name, value, wavelengths=standard_wavelengths):
if np.iterable(value):
generated by cgit (git 2.34.1) at 2025-05-07 22:52:28 +0000