// The maximum number of steps a ray can take during marching before giving up
// and colliding with nothing. This prevents scenes from taking infinite time to render.
#define MAX_STEPS 300
// The maximum distance a ray can travel before we give up and just say it collides
// with nothing. This helps prevent the background from appearing warped by the foreground
// due to rays which march close to a foreground object run out of steps before
// reaching their destination when slightly farther rays do reach their target.
#define MAX_DIST 50.
// The minimum distance between two points before they are considered the same point.
// Setting a minimum distance prevents graphical glitches when ray marching parallel
// to a surface, where the ray does not intersect an object, but comes close enough
// that the march becomes so slow that it fails to reach its actual destination.
//
// This is arbitrarily chosen to be 2^(-9) meters, or ~2 mm.
// (decreased by experimentation)
#define MIN_DIST     (0.001953125/4.)
// The distance between samples when estimating a surface's normal.
// This is arbitrarily chosen to be 2x the MIN_DIST.
#define NORMAL_DELTA (MIN_DIST*2.)
//#define NORMAL_DELTA (0.00390625)
// The maximum number of recursive steps taken while calculating light.
#define LIGHT_STEPS 4
#define LIGHT_MIN_DIST MIN_DIST*2.
#define SAMPLES 1

// borrowed from IQ: https://www.shadertoy.com/view/4sfGzS
//------------------------------------------------------------------
// oldschool rand() from Visual Studio
//------------------------------------------------------------------
int  seed = 1;
void srand(int s ) { seed = s; }
int  rand(void) { seed = seed*0x343fd+0x269ec3; return (seed>>16)&32767; }
float frand(void) { return float(rand())/32767.0; }
//------------------------------------------------------------------
// hash to initialize the random sequence (copied from Hugo Elias)
//------------------------------------------------------------------
int hash( int n )
{
	n = (n << 13) ^ n;
    return n * (n * n * 15731 + 789221) + 1376312589;
}
//------------------------------------------------------------------
// end borrowed from IQ

// by fizzer via IQ: http://www.amietia.com/lambertnotangent.html
vec3 cosine_direction(vec3 norm) {
    float u = frand();
    float v = frand();

    // method 3 
    float a = 6.2831853 * v;
    u = 2.0*u - 1.0;
    return normalize(norm + vec3(sqrt(1.0 - u*u)*vec2(cos(a), sin(a)), u));
}
// end by fizzer

// The distance from a point to the nearest object in the scene.
float dist(vec3 pos) {
    float sphere = distance(pos, vec3(0.0, -0.2, 8.)) - 1.;
    float plane = pos.y + 1.0;
    return min(sphere, plane);
}

// Estimate the angle from a point to the nearest surface.
vec3 normal(vec3 pos) {
    vec2 delta = vec2(NORMAL_DELTA, 0.);
    vec3 dq = (dist(pos) - vec3(
        dist(pos - delta.xyy),
        dist(pos - delta.yxy),
        dist(pos - delta.yyx)
    )) / delta.x;
    
    // This function is concerningly incorrect without normalize.
    // Maybe I can find a better approximation?
    return normalize(dq);
}

// Approximate the distance to the nearest object along a ray.
#if 1
float march(vec3 origin, vec3 direction) {
    float total_dist = 0.;
    float delta = 1./0.;
    int steps = 0;
    for (; steps < MAX_STEPS && total_dist < MAX_DIST && delta >= MIN_DIST; steps++) {
        vec3 pos = origin + direction * total_dist;
        delta = dist(pos);
        total_dist += delta;
    }
    return delta < MIN_DIST ? total_dist : 1./0.;
}
#else
float march(vec3 origin, vec3 direction) {
    float total_dist = 0.;
    float delta = 1./0.;
    int steps = 0;
    vec3 pos = origin;
    for (; steps < MAX_STEPS && distance(pos, origin) < MAX_DIST && delta >= MIN_DIST; steps++) {
        delta = dist(pos);
        pos += direction * delta;
    }
    return distance(pos, origin) < MAX_DIST && steps < MAX_STEPS ? distance(pos, origin) : 1./0.;
}
#endif

float light(vec3 pos, vec3 eye) {
    vec3 source = vec3(1.5, 1.5, 6.);
    
    float light = 0.;
    int steps = 1;
    for (; steps < LIGHT_STEPS + 1; steps++) {
        vec3 dir = normalize(source - pos);
        
        vec3 norm = normal(pos);

        // Diffuse lighting
        float diff = dot(norm, dir);
        diff *= max(0., dot(norm, -eye));
        diff = max(0., diff);

        // Specular lighting
        vec3 refl = normalize(2.*dot(dir, norm)*norm - dir);
        float spec = pow(max(0., dot(refl, -eye)), 5.);
        
        float direct = diff + spec;
        float shadowed = distance(source, pos) - march(source, -dir);
        // magic number 16 determined by experimentation
        if (shadowed < MIN_DIST*16.) light += direct;
        
        // Indirect lighting
        eye = cosine_direction(norm);
        pos += eye*LIGHT_MIN_DIST;
        float next = march(pos, eye);
        if (isinf(next)) break;
        pos += eye * next;
    }
    
    return min(1., light);
}

vec3 project(vec2 uv) {
    return vec3(uv, sqrt(1. - uv.x*uv.x - uv.y*uv.y));
}

vec3 linear2srgb(vec3 linear_rgb) {
    // copied from somewhere on the internet
    bvec3 cutoff = lessThan(linear_rgb, vec3(0.0031308));
    vec3 higher = vec3(1.055)*pow(linear_rgb, vec3(1.0/2.4)) - vec3(0.055);
    vec3 lower = linear_rgb * vec3(12.92);

    return mix(higher, lower, cutoff);
    // end copied from somewhere on the internet
}

void mainImage(out vec4 fragColor, in vec2 fragCoord) {
    // borrowed from IQ
    // init randoms
    ivec2 q = ivec2(fragCoord);
    srand( hash(q.x+hash(q.y+hash(iFrame))));
    // end borrowed from IQ

    vec3 col = vec3(0.);
    for (int i = 0; i < SAMPLES; i++) {
        vec2 coord = fragCoord + vec2(frand(), frand() - 0.5);
        vec2 uv = ((coord/iResolution.xy) - 0.5) * 2.;
        // Preserve aspect ratio when x > y.
        uv = vec2(uv.x, uv.y * iResolution.y / iResolution.x);
        // FOV
        uv /= 2.;
        
        vec3 ray = project(uv);
        float dist = march(vec3(0.), ray);
        if (isinf(dist)) continue;
        
        vec3 pos = dist * ray;
        vec3 norm = normal(pos);
        
        col += light(pos, ray);
    }
    
    col /= float(SAMPLES);
    
    
    // Transform color space from linear RGB to sRGB.
    col = linear2srgb(col);
    fragColor = vec4(col, 1.0);
}