propagate() takes an array of photon offsets and a range of

offsets to load.  Now events with more photons than RNG states
can be propagated through multiple kernel calls.

Also lays the groundwork for consolidating photons between steps
to reduce the amount of propagation work required.

2 files changed, 24 insertions, 19 deletions

diff --git a/src/daq.cu b/src/daq.cu
index 2b5f9b4..7c5e6a5 100644
--- a/src/daq.cu
+++ b/src/daq.cu
@@ -30,6 +30,7 @@ __global__ void reset_earliest_time_int(float maxtime, int ntime_ints, unsigned
 
 __global__ void run_daq(curandState *s, unsigned int detection_state,
 			float time_rms,
+			int first_photon,
 			int nphotons, float *photon_times,
 			unsigned int *photon_histories,
 			int *last_hit_triangles, int *solid_map,
@@ -41,17 +42,17 @@ __global__ void run_daq(curandState *s, unsigned int detection_state,
 	if (id < nphotons)
 	{
 		curandState rng = s[id];
-
-		int triangle_id = last_hit_triangles[id];
+		int photon_id = id + first_photon;
+		int triangle_id = last_hit_triangles[photon_id];
 		
 		if (triangle_id > -1)
 		{
 			int solid_id = solid_map[triangle_id];
-			int history = photon_histories[id];
+			int history = photon_histories[photon_id];
 
 			if (solid_id < nsolids && (history & detection_state)) 
 			{
-			  float time = photon_times[id] + curand_normal(&rng) * time_rms;
+			  float time = photon_times[photon_id] + curand_normal(&rng) * time_rms;
 			  unsigned int time_int = float_to_sortable_int(time);
 			  
 			  atomicMin(earliest_time_int + solid_id, time_int);
diff --git a/src/kernel.cu b/src/kernel.cu
index 82efe88..17b829c 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -260,7 +260,9 @@ __global__ void ray_trace(int nthreads, float3 *positions, float3 *directions, i
 
 } // ray_trace
 
-__global__ void propagate(int nthreads, curandState *rng_states, float3 *positions, float3 *directions, float *wavelengths, float3 *polarizations, float *times, unsigned int *histories, int *last_hit_triangles, int max_steps)
+__global__ void propagate(int first_photon, int nthreads, int *photon_offsets, curandState *rng_states, 
+			  float3 *positions, float3 *directions, float *wavelengths, float3 *polarizations, float *times, 
+			  unsigned int *histories, int *last_hit_triangles, int max_steps)
 {
 	int id = blockIdx.x*blockDim.x + threadIdx.x;
 
@@ -269,16 +271,18 @@ __global__ void propagate(int nthreads, curandState *rng_states, float3 *positio
 
 	curandState rng = rng_states[id];
 
+	int photon_id = photon_offsets[first_photon + id];
+
 	Photon p;
-	p.position = positions[id];
-	p.direction = directions[id];
+	p.position = positions[photon_id];
+	p.direction = directions[photon_id];
 	p.direction /= norm(p.direction);
-	p.polarization = polarizations[id];
+	p.polarization = polarizations[photon_id];
 	p.polarization /= norm(p.polarization);
-	p.wavelength = wavelengths[id];
-	p.time = times[id];
-	p.last_hit_triangle = last_hit_triangles[id];
-	p.history = histories[id];
+	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;
@@ -321,13 +325,13 @@ __global__ void propagate(int nthreads, curandState *rng_states, float3 *positio
 	} // while (steps < max_steps)
 
 	rng_states[id] = rng;
-	positions[id] = p.position;
-	directions[id] = p.direction;
-	polarizations[id] = p.polarization;
-	wavelengths[id] = p.wavelength;
-	times[id] = p.time;
-	histories[id] = p.history;
-	last_hit_triangles[id] = p.last_hit_triangle;
+	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;
 
 } // propagate