1 files changed, 100 insertions, 0 deletions
diff --git a/src/propagate.cu b/src/propagate.cu
new file mode 100644
index 0000000..2d5183e
--- /dev/null
+++ b/src/propagate.cu
@@ -0,0 +1,100 @@
+//-*-c-*-
+
+#include "linalg.h"
+#include "mesh.h"
+#include "geometry.h"
+#include "photon.h"
+
+extern "C"
+{
+
+__global__ void
+propagate(int first_photon, int nthreads, unsigned int *input_queue,
+	  unsigned int *output_queue, curandState *rng_states,
+	  float3 *positions, float3 *directions,
+	  float *wavelengths, float3 *polarizations,
+	  float *times, unsigned int *histories,
+	  int *last_hit_triangles, int max_steps,
+	  Geometry *geometry)
+{
+    int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+    if (id >= nthreads)
+	return;
+
+    curandState rng = rng_states[id];
+
+    int photon_id = input_queue[first_photon + id];
+
+    Photon p;
+    p.position = positions[photon_id];
+    p.direction = directions[photon_id];
+    p.direction /= norm(p.direction);
+    p.polarization = polarizations[photon_id];
+    p.polarization /= norm(p.polarization);
+    p.wavelength = wavelengths[photon_id];
+    p.time = times[photon_id];
+    p.last_hit_triangle = last_hit_triangles[photon_id];
+    p.history = histories[photon_id];
+
+    if (p.history & (NO_HIT | BULK_ABSORB | SURFACE_DETECT | SURFACE_ABSORB))
+	return;
+
+    State s;
+
+    int steps = 0;
+    while (steps < max_steps) {
+	steps++;
+
+	int command;
+
+	// check for NaN and fail
+	if (isnan(p.direction.x*p.direction.y*p.direction.z*p.position.x*p.position.y*p.position.z)) {
+	    p.history |= NO_HIT | NAN_ABORT;
+	    break;
+	}
+
+	fill_state(s, p, geometry);
+
+	if (p.last_hit_triangle == -1)
+	    break;
+
+	command = propagate_to_boundary(p, s, rng);
+
+	if (command == BREAK)
+	    break;
+
+	if (command == CONTINUE)
+	    continue;
+
+	if (s.surface_index != -1) {
+	    command = propagate_at_surface(p, s, rng, geometry);
+
+	    if (command == BREAK)
+		break;
+
+	    if (command == CONTINUE)
+		continue;
+	}
+
+	propagate_at_boundary(p, s, rng);
+
+    } // while (steps < max_steps)
+
+    rng_states[id] = rng;
+    positions[photon_id] = p.position;
+    directions[photon_id] = p.direction;
+    polarizations[photon_id] = p.polarization;
+    wavelengths[photon_id] = p.wavelength;
+    times[photon_id] = p.time;
+    histories[photon_id] = p.history;
+    last_hit_triangles[photon_id] = p.last_hit_triangle;
+	
+    // Not done, put photon in output queue
+    if ((p.history & (NO_HIT | BULK_ABSORB | SURFACE_DETECT | SURFACE_ABSORB)) == 0) {
+	int out_idx = atomicAdd(output_queue, 1);
+	output_queue[out_idx] = photon_id;
+    }
+} // propagate
+
+} // extern "C"