1 files changed, 273 insertions, 0 deletions
diff --git a/src/photon.h b/src/photon.h
new file mode 100644
index 0000000..da866bc
--- /dev/null
+++ b/src/photon.h
@@ -0,0 +1,273 @@
+#ifndef __PHOTON_H__
+#define __PHOTON_H__
+
+#include "linalg.h"
+#include "materials.h"
+#include "rotate.h"
+#include "random.h"
+#include "physical_constants.h"
+#include "mesh.h"
+
+struct Photon
+{
+	float3 position;
+	float3 direction;
+	float3 polarization;
+	float wavelength;
+	float time;
+
+	//bool alive;
+	//unsigned int last_process;
+	unsigned int history;
+
+	int last_hit_triangle;
+
+	//int photon_state;
+
+	curandState rng;
+};
+
+struct State
+{
+	float3 surface_normal;
+	//Material material1, material2;
+	//Surface surface;
+
+	float refractive_index1, refractive_index2;
+	float absorption_length;
+	float scattering_length;
+	//float absorption_distance;
+	//float scattering_distance;
+
+	int surface_index;
+
+	//int last_hit_triangle;
+	float distance_to_boundary;
+};
+
+enum
+{
+	//DEBUG = -2,
+	//INIT = -1,
+	NO_HIT           = 0x1 << 0,
+	BULK_ABSORB      = 0x1 << 1,
+	SURFACE_DETECT   = 0x1 << 2,
+	SURFACE_ABSORB   = 0x1 << 3,
+	RAYLEIGH_SCATTER = 0x1 << 4,
+	REFLECT_DIFFUSE  = 0x1 << 5,
+	REFLECT_SPECULAR = 0x1 << 6
+}; // processes
+
+enum {BREAK, CONTINUE, PASS}; // return value from propagate_to_boundary
+
+__device__ int fill_state(State &s, Photon &p)
+{
+	p.last_hit_triangle = intersect_mesh(p.position, p.direction, s.distance_to_boundary, p.last_hit_triangle);
+
+	if (p.last_hit_triangle == -1)
+	{
+		p.history |= NO_HIT;
+		return BREAK;
+	}
+
+	uint4 triangle_data = g_triangles[p.last_hit_triangle];
+
+	float3 v0 = g_vertices[triangle_data.x];
+	float3 v1 = g_vertices[triangle_data.y];
+	float3 v2 = g_vertices[triangle_data.z];
+
+	int inner_material_index = convert(0xFF & (triangle_data.w >> 24));
+	int outer_material_index = convert(0xFF & (triangle_data.w >> 16));
+	s.surface_index = convert(0xFF & (triangle_data.w >> 8));
+
+	s.surface_normal = cross(v1-v0, v2-v1);
+	s.surface_normal /= norm(s.surface_normal);
+
+	Material material1, material2;
+	if (dot(s.surface_normal,-p.direction) > 0.0f)
+	{
+		// outside to inside
+		material1 = materials[outer_material_index];
+		material2 = materials[inner_material_index];
+	}
+	else
+	{
+		// inside to outside
+		material1 = materials[inner_material_index];
+		material2 = materials[outer_material_index];
+		s.surface_normal = -s.surface_normal;
+	}
+
+	s.refractive_index1 = interp_property(p.wavelength, material1.refractive_index);
+	s.refractive_index2 = interp_property(p.wavelength, material2.refractive_index);
+	s.absorption_length = interp_property(p.wavelength, material1.absorption_length);
+	s.scattering_length = interp_property(p.wavelength, material1.scattering_length);
+
+	return PASS;
+
+} // fill_state
+
+__device__ void rayleigh_scatter(Photon &p)
+{
+	float theta, y;
+
+	while (true)
+	{
+		y = curand_uniform(&p.rng);
+		theta = uniform(&p.rng, 0, 2*PI);
+
+		if (y < powf(cosf(theta),2))
+			break;
+	}
+
+	float phi = uniform(&p.rng, 0, 2*PI);
+
+	float3 b = cross(p.polarization, p.direction);
+	float3 c = p.polarization;
+
+	p.direction = rotate(p.direction, theta, b);
+	p.direction = rotate(p.direction, phi, c);
+
+	p.polarization = rotate(p.polarization, theta, b);
+	p.polarization = rotate(p.polarization, phi, c);
+
+} // scatter
+
+__device__ int propagate_to_boundary(Photon &p, State &s)
+{
+	float absorption_distance = -s.absorption_length*logf(curand_uniform(&p.rng));
+	float scattering_distance = -s.scattering_length*logf(curand_uniform(&p.rng));
+
+	if (absorption_distance <= scattering_distance)
+	{
+		if (absorption_distance <= s.distance_to_boundary)
+		{
+			p.time += absorption_distance/(SPEED_OF_LIGHT/s.refractive_index1);
+			p.position += absorption_distance*p.direction;
+			p.history |= BULK_ABSORB;
+
+			p.last_hit_triangle = -1;
+
+			return BREAK;
+		} // photon is absorbed in material1
+	}
+	else
+	{
+		if (scattering_distance <= s.distance_to_boundary)
+		{
+			p.time += scattering_distance/(SPEED_OF_LIGHT/s.refractive_index1);
+			p.position += scattering_distance*p.direction;
+
+			rayleigh_scatter(p);
+
+			p.history |= RAYLEIGH_SCATTER;
+
+			p.last_hit_triangle = -1;
+
+			return CONTINUE;
+		} // photon is scattered in material1
+	} // if scattering_distance < absorption_distance
+
+	p.position += s.distance_to_boundary*p.direction;
+	p.time += s.distance_to_boundary/(SPEED_OF_LIGHT/s.refractive_index1);
+
+	return PASS;
+
+} // propagate_to_boundary
+
+__device__ void propagate_at_boundary(Photon &p, State &s)
+{
+	float incident_angle = acosf(dot(s.surface_normal, -p.direction));
+	float refracted_angle = asinf(sinf(incident_angle)*s.refractive_index1/s.refractive_index2);
+
+	float3 incident_plane_normal = cross(p.direction, s.surface_normal);
+	incident_plane_normal /= norm(incident_plane_normal);
+
+	float normal_coefficient = dot(p.polarization, incident_plane_normal);
+	float normal_probability = normal_coefficient*normal_coefficient;
+
+	float reflection_coefficient;
+	if (curand_uniform(&p.rng) < normal_probability)
+	{
+		reflection_coefficient = -sinf(incident_angle-refracted_angle)/sinf(incident_angle+refracted_angle);
+
+		p.polarization = incident_plane_normal;
+	} // photon polarization normal to plane of incidence
+	else
+	{
+		reflection_coefficient = tanf(incident_angle-refracted_angle)/tanf(incident_angle+refracted_angle);
+
+		p.polarization = cross(incident_plane_normal, p.direction);
+		p.polarization /= norm(p.polarization);
+	} // photon polarization parallel to plane of incidence
+
+	if ((curand_uniform(&p.rng) < reflection_coefficient*reflection_coefficient) || isnan(refracted_angle))
+	{
+		p.direction = rotate(s.surface_normal, incident_angle, incident_plane_normal);
+
+		p.history |= REFLECT_SPECULAR;
+	}
+	else
+	{
+		p.direction = rotate(s.surface_normal, PI-refracted_angle, incident_plane_normal);
+	}
+
+} // propagate_at_boundary
+
+__device__ int propagate_at_surface(Photon &p, State &s)
+{
+	Surface surface = surfaces[s.surface_index];
+
+	float detect = interp_property(p.wavelength, surface.detect);
+	float absorb = interp_property(p.wavelength, surface.absorb);
+	float reflect_diffuse = interp_property(p.wavelength, surface.reflect_diffuse);
+	float reflect_specular = interp_property(p.wavelength, surface.reflect_specular);
+
+	// since the surface properties are interpolated linearly, we are
+	// guaranteed that they still sum to 1.0.
+
+	float uniform_sample = curand_uniform(&p.rng);
+
+	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)
+	{
+		// diffusely reflect
+		p.direction = uniform_sphere(&p.rng);
+
+		if (dot(p.direction, s.surface_normal) < 0.0f)
+			p.direction = -p.direction;
+
+		// randomize polarization?
+		p.polarization = cross(uniform_sphere(&p.rng), p.direction);
+		p.polarization /= norm(p.polarization);
+
+		p.history |= REFLECT_DIFFUSE;
+
+		return CONTINUE;
+	}
+	else
+	{
+		// specularly reflect
+		float incident_angle = acosf(dot(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;
+	}
+
+} // propagate_at_surface
+
+#endif