Skip to content

Instantly share code, notes, and snippets.

Created April 24, 2021 15:24
What would you like to do?
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define NUM_SAMPLES 100000
typedef struct{
float x, y, z;
} vec3;
/*** Vector functions ***/
vec3 V3 (float x, float y, float z);
vec3 normalize(vec3 v);
float dot (vec3 a, vec3 b);
/*** math helper ***/
float sq (float x);
/* max x with to a*x*x+b*x+c==0 */
float solvequad(float a, float b, float c);
/*** ray handling ***/
vec3 refract (vec3 dir, vec3 normal, float n_in, float n_out);
vec3 ray_hit (vec3 origin, vec3 dir);
int main(int argc, char** argv){
if (argc < 3){
fprintf(stderr, "Usage: granular-transmission <refraction index> <output file>\n");
return 1;
float n = atof (argv[1]);
FILE* outfile = fopen(argv[2], "wb");
for(int i = 0; i != NUM_SAMPLES; i++){
/* Don't take samples at x = g_r. This would need more careful handling of edgecases */
float x = 1.0f*i/NUM_SAMPLES;
vec3 ray_incoming = V3(0, 0, 1);
vec3 ray_incoming_hit = V3(x, 0, -sqrtf(1-sq(x)));
vec3 normal_incoming = normalize(ray_incoming_hit);
vec3 ray_internal = refract(ray_incoming, normal_incoming, 1.0f, n);
vec3 ray_outgoing_hit = ray_hit(ray_incoming_hit, ray_internal);
vec3 normal_outgoing = normalize(ray_outgoing_hit);
vec3 inormal_outgoing = V3(-normal_outgoing.x, -normal_outgoing.y, -normal_outgoing.z);
vec3 ray_outgoing = refract(ray_internal, inormal_outgoing, n, 1.0f);
float angle_outgoing = 3.141592-acos(-ray_outgoing.z);
fprintf(outfile, "%f\n", angle_outgoing);
vec3 V3(float x, float y, float z){
vec3 v;
v.x = x;
v.y = y;
v.z = z;
return v;
vec3 normalize(vec3 v){
float fac = 1.0f/sqrtf(dot(v, v));
return V3(v.x*fac, v.y*fac, v.z*fac);
float dot(vec3 a, vec3 b){
return a.x*b.x+a.y*b.y+a.z*b.z;
float sq(float x){
return x*x;
float solvequad(float a, float b, float c){
/* This could be done faster! */
float rt = sqrtf(b*b-4*a*c);
float s1 = (-b+rt)/(2*a);
float s2 = (-b-rt)/(2*a);
return s1>s2 ? s1 : s2;
vec3 refract(vec3 dir, vec3 normal, float n_in, float n_out){
/* */
float r = n_in/n_out;
float c = -dot(normal, dir);
float nfac = r*c-sqrt(1-sq(r)*(1-sq(c)));
return normalize(V3(r*dir.x+nfac*normal.x,
vec3 ray_hit(vec3 origin, vec3 dir){
/* solve (origin+dir)^2 == 1^2 */
float len = solvequad(dot(dir, dir), 2*dot(dir, origin), dot(origin, origin)-1);
return V3(origin.x+len*dir.x,
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment