Add dielectric material

src/material.c

@@ -1,5 +1,24 @@
 #include "material.h"
 
+#include <math.h>
+#include <stdio.h>
+
+static vec3_t reflect(vec3_t v, vec3_t n)
+{
+	const vec3_t dn = vec3_scale(n, -2 * vec3_dot(n, v));
+	return vec3_unit(vec3_add(v, dn));
+}
+
+static vec3_t refract(vec3_t v, vec3_t n, double eta_ratio)
+{
+	const double cos_theta = fmin(-vec3_dot(v, n), 1.0);
+	const vec3_t perp
+	    = vec3_scale(vec3_add(v, vec3_scale(n, cos_theta)), eta_ratio);
+	const vec3_t para
+	    = vec3_scale(n, -sqrt(fabs(1.0 - vec3_len_squared(perp))));
+	return vec3_unit(vec3_add(perp, para));
+}
+
 bool scatter_lambertian(
     material_params_t params, hit_t hit, rng_t *rng, ray_t *ray,
     vec3_t *atten_out)
@@ -17,10 +36,30 @@ bool scatter_reflective(
 {
 	(void)rng;
 
-	const vec3_t dn
-	    = vec3_scale(hit.normal, -2 * vec3_dot(hit.normal, ray->dir));
 	ray->orig = hit.point;
-	ray->dir = vec3_unit(vec3_add(ray->dir, dn));
+	ray->dir = reflect(ray->dir, hit.normal);
 	*atten_out = params.reflective.tint;
 	return true;
 }
+
+bool scatter_dielectric(
+    material_params_t params, hit_t hit, rng_t *rng, ray_t *ray,
+    vec3_t *atten_out)
+{
+	const double eta_ratio
+	    = hit.front ? 1.0 / params.dielectric.eta : params.dielectric.eta;
+	const vec3_t n = hit.front ? hit.normal : vec3_scale(hit.normal, -1);
+	const double cos_theta = fmin(-vec3_dot(ray->dir, n), 1.0);
+	const double sin_theta = sqrt(1.0 - cos_theta * cos_theta);
+	const double r = pow((1 - eta_ratio) / (1 + eta_ratio), 2);
+	const double reflectance = r + (1 - r) * pow(1 - cos_theta, 5);
+
+	ray->orig = hit.point;
+	if (sin_theta * eta_ratio > 1.0 || rng_canon(rng) < reflectance)
+		ray->dir = reflect(ray->dir, n);
+	else
+		ray->dir = refract(ray->dir, n, eta_ratio);
+	*atten_out = (vec3_t) { 1.0, 1.0, 1.0 };
+
+	return true;
+}

src/rng.c

@@ -18,14 +18,9 @@ static uint32_t next(rng_t *rng)
 	return rng->state = x;
 }
 
-static double canonical(rng_t *rng)
-{
-	return (double)next(rng) / (double)UINT32_MAX;
-}
-
 static double disc(rng_t *rng)
 {
-	return 2.0 * canonical(rng) - 1.0;
+	return 2.0 * rng_canon(rng) - 1.0;
 }
 
 rng_t rng_init(unsigned seed)
@@ -35,6 +30,11 @@ rng_t rng_init(unsigned seed)
 	return (rng_t) { .state = tv.tv_usec + seed };
 }
 
+double rng_canon(rng_t *rng)
+{
+	return (double)next(rng) / (double)UINT32_MAX;
+}
+
 vec3_t rng_vec3(rng_t *rng)
 {
 	const vec3_t v = { disc(rng), disc(rng), disc(rng) };
@@ -43,10 +43,10 @@ vec3_t rng_vec3(rng_t *rng)
 
 vec3_t rng_gaussian_xy(rng_t *rng, double stddev)
 {
-	const double r1 = canonical(rng);
+	const double r1 = rng_canon(rng);
 	double r2;
 	do
-		r2 = canonical(rng);
+		r2 = rng_canon(rng);
 	while (r2 == 0);
 
 	const double theta = 2.0 * M_PI * r1;

src/vec3.c

@@ -49,9 +49,14 @@ vec3_t vec3_hadamard(vec3_t v, vec3_t u)
 	};
 }
 
+double vec3_len_squared(vec3_t v)
+{
+	return v.x * v.x + v.y * v.y + v.z * v.z;
+}
+
 double vec3_len(vec3_t v)
 {
-	return sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
+	return sqrt(vec3_len_squared(v));
 }
 
 vec3_t vec3_unit(vec3_t v)