Skip to content

Instantly share code, notes, and snippets.

@a-voel
Created April 24, 2021 15:24
Show Gist options
  • Save a-voel/59a6b17966bebb5da9061c51a07a264d to your computer and use it in GitHub Desktop.
Save a-voel/59a6b17966bebb5da9061c51a07a264d to your computer and use it in GitHub Desktop.
#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);
}
fclose(outfile);
}
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){
/* https://en.wikipedia.org/wiki/Snell%27s_law#Vector_form */
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,
r*dir.y+nfac*normal.y,
r*dir.z+nfac*normal.z));
}
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,
origin.y+len*dir.y,
origin.z+len*dir.z);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment