Created
April 24, 2021 15:24
-
-
Save a-voel/59a6b17966bebb5da9061c51a07a264d to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#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