#include <stdlib.h>
#include <stdio.h>

#define _USE_MATH_DEFINES // so we can get M_PI
#include <math.h>


// in this raytracer the following basis convection is used:
// we imagine the x axis and y axis as usual: the x axis to the right, the y axis upward
// we place the z axis "towards" ourselves

// random float in [0,1]
float randf()
{
    return (float)rand()/(float)(RAND_MAX);
}

static inline int maxi(const int a, const int b) 
{
    return a > b ? a : b;
}

static inline int mini(const int a, const int b) 
{
    return a < b ? a : b;
}

static inline int clampi(const int a, const int min, const int max) 
{
    return mini(maxi(a, min), max);
}

static inline float maxf(const float a, const float b) 
{
    return a > b ? a : b;
}

static inline float minf(const float a, const float b) 
{
    return a < b ? a : b;
}

struct vec3
{
    float x, y, z;
};
typedef struct vec3 vec3; // allows us to write just 'vec3' instead of 'struct vec3'

static inline vec3 vec3_add(const vec3 a, const vec3 b)
{
    return (vec3){a.x + b.x, a.y + b.y, a.z + b.z};
}

static inline vec3 vec3_sub(const vec3 a, const vec3 b)
{
    return (vec3){a.x - b.x, a.y - b.y, a.z - b.z};
}

static inline vec3 vec3_mul(const vec3 a, const float b)
{
    return (vec3){a.x * b, a.y * b, a.z * b};
}

static inline vec3 vec3_mulv(const vec3 a, const vec3 b)
{
    return (vec3){a.x * b.x, a.y * b.y, a.z * b.z};
}

static inline vec3 vec3_div(const vec3 a, const float b)
{
    return (vec3){a.x / b, a.y / b, a.z / b};
}

static inline float vec3_dot(const vec3 a, const vec3 b)
{
    return a.x * b.x + a.y * b.y + a.z * b.z; 
}

static inline float vec3_length(const vec3 a)
{
    return sqrtf(a.x * a.x + a.y * a.y + a.z * a.z);
}

static inline vec3 vec3_normalize(const vec3 a)
{
    return vec3_div(a, vec3_length(a));
}

static inline vec3 vec3_cross(const vec3 a, const vec3 b)
{
    return (vec3){
        a.y * b.z - a.z * b.y, 
        a.z * b.x - a.x * b.z, 
        a.x * b.y - a.y * b.x
    };
}

static inline vec3 vec3_reflect(const vec3 incidence, const vec3 normal)
{
    return vec3_sub(incidence, vec3_mul(normal, 2.0 * vec3_dot(incidence, normal)));
}

static inline float vec3_distance(const vec3 a, const vec3 b)
{
    return vec3_length(vec3_sub(a,b));
}

// sphere is defined as all points x for which
// dis(x, center) = rad
//
// the sphere_intersect function solves
// dis(ori + dir * t, cen) = rad for the smallest positive t possible
//
// define del = ori - cen
// we can rewrite the above as
// ||del + dir * t||^2 = rad^2
// dot(del,del) + 2 * dot(del, dir) * t + dot(dir, dir) * t^2 = rad^2
//
// this is a quadratic equation of the form
// C + B * t + A * t^2 = 0
// with C = dot(del,del) - rad^2
// B = 2 * dot(del, dir)
// A = dot(dir, dir)
//
// so then use the 'abc-formula'
float sphere_intersect(
    const vec3 center, 
    const float radius, 
    const vec3 origin, 
    const vec3 direction
)
{
    const vec3 delta = vec3_sub(origin, center);

    const float A = vec3_dot(direction, direction);
    const float B = 2.0 * vec3_dot(delta, direction);
    const float C = vec3_dot(delta, delta) - radius * radius;

    const float dis = B * B - 4.0 * A * C;

    if(dis < 0.0) NAN; // no solutions

    const float sqrt_disc = sqrtf(dis);

    const float t1 = (- B - sqrt_disc) / (2.0 * A);
    const float t2 = (- B + sqrt_disc) / (2.0 * A);

    if( t2 < 0.0 ) return NAN; // no positive solutions

    if ( t1 > 0.0 ) return t1;

    return t2;
}

// plane is defined as all point x for which
// dot(x - center, normal) = 0
//
// the plane_intersect function solves
// dot(ori + dir * t - center, normal) = 0 for the smallest positive t possible
//
// define del = ori - cen
// we can rewrite the above as
// dot(del, nor) + t * dot(dir, nor) = 0
// t = - dot(del, nor) / dot(dir, nor)
float plane_intersect(
    const vec3 center,
    const vec3 normal,
    const vec3 origin, 
    const vec3 direction
)
{
    const vec3 delta = vec3_sub(origin, center);
    const float t = - vec3_dot(delta, normal) / vec3_dot(direction, normal);

    if( t < 0.0 ) return NAN; // no positive solution

    return t;
}



struct hit_info
{
    float t;
    vec3 color, normal, position, reflectance;
};
typedef struct hit_info hit_info; // allows us to write just 'hit_info' instead of 'struct hit_info'


hit_info scene_intersect( 
    const vec3 origin, 
    const vec3 direction
)
{
    const vec3 sphere1_center = {0.0, 0.0, -2.0};
    const float sphere1_radius = 1.0;
    const vec3 sphere1_color = {0.75, 0.25, 0.25};
    const vec3 sphere1_reflectance = {0.05, 0.05, 0.05};

    const vec3 sphere2_center = {0.0, 0.0, 2.0};
    const float sphere2_radius = 1.0;
    const vec3 sphere2_color = {0.75, 0.25, 0.25};
    const vec3 sphere2_reflectance = {1.0, 1.0, 1.0};

    const vec3 plane_center = {0.0, -1.0, 0.0};
    const vec3 plane_normal = {0.0, 1.0, 0.0};
    const vec3 plane_color1 = {0.75, 0.75, 0.75};
    const vec3 plane_color2 = {0.25, 0.25, 0.25};
    const vec3 plane_reflectance = {0.0, 0.0, 0.0};

    hit_info hit = {
        INFINITY, 
        (vec3){0.0,0.0,0.0}, 
        (vec3){0.0,0.0,0.0}, 
        (vec3){0.0,0.0,0.0},
        (vec3){0.0,0.0,0.0}
    };

    float t;
    t = sphere_intersect(sphere1_center, sphere1_radius, origin, direction);

    if(!isnan(t) && t < hit.t)
    {
        hit.t = t;
        hit.position = vec3_add(origin, vec3_mul(direction, t));
        hit.normal = vec3_div(vec3_sub(hit.position, sphere1_center), sphere1_radius);
        hit.color = sphere1_color;
        hit.reflectance = sphere1_reflectance;
    }

    t = sphere_intersect(sphere2_center, sphere2_radius, origin, direction);

    if(!isnan(t) && t < hit.t)
    {
        hit.t = t;
        hit.position = vec3_add(origin, vec3_mul(direction, t));
        hit.normal = vec3_div(vec3_sub(hit.position, sphere2_center), sphere2_radius);
        hit.color = sphere2_color;
        hit.reflectance = sphere2_reflectance;
    }

    t = plane_intersect(plane_center, plane_normal, origin, direction);
    if(!isnan(t) && t < hit.t)
    {
        hit.t = t;
        hit.position = vec3_add(origin, vec3_mul(direction, t));
        hit.normal = plane_normal;
        hit.reflectance = plane_reflectance;

        if( cos( 0.25 * hit.position.x * M_PI) + cos( 0.25 * hit.position.z * M_PI) < 0.0 )
        {
            hit.color = plane_color1;
        }else
        {
            hit.color = plane_color2;
        }
    }

    return hit;
}

vec3 schlick(
    const float costheta, 
    const vec3 reflectance
)
{
    const float f = powf(1.0 - costheta, 5.0);
    return vec3_add((vec3){f,f,f}, vec3_mul(reflectance, 1.0 - f) );
}

const int MAX_DEPTH = 10;

vec3 raytrace(
    const vec3 origin,
    const vec3 direction,
    int depth
)
{    
    const vec3 light_position = {10.0, 10.0, 10.0};
    const vec3 light_color = {200.0, 200.0, 200.0};

    const vec3 sky_color = {0.25, 0.6, 0.8};

    vec3 light = {0.0, 0.0, 0.0};

    hit_info hit = scene_intersect(origin, direction);

    if(depth > MAX_DEPTH) return light;

    if(isfinite(hit.t))
    {
        vec3 light_delta = vec3_sub(light_position, hit.position);
        float light_distance = vec3_length(light_delta);
        vec3 light_direction = vec3_div(light_delta, light_distance);

        float incidence = vec3_dot(vec3_mul(direction, -1.0), hit.normal);
        vec3 reflectance = schlick(incidence, hit.reflectance);

        float dot = vec3_dot(hit.normal, light_direction);
        float diffuse = maxf(dot, 0.0);

        float shadow = 1.0;
        if(diffuse > 0.0) // if we even has a chance if hitting the light source
        {
            float epsilon = 1e-5;
            vec3 shadow_ray_origin = vec3_add(hit.position, vec3_mul(hit.normal, epsilon)); // start the shadow ray a little bit away from the surface
            hit_info shadow_check = scene_intersect(shadow_ray_origin, light_direction);
            if(shadow_check.t < light_distance)
            {
                shadow = 0.0;
            }
        }

        vec3 light_from_diffuse = vec3_mulv(hit.color, vec3_add(sky_color, vec3_mul(light_color, shadow * diffuse / (light_distance * light_distance))));
        
        float epsilon = 1e-5;
        vec3 reflect_ray_origin = vec3_add(hit.position, vec3_mul(hit.normal, epsilon));
        vec3 reflect_ray_direction = vec3_reflect(direction, hit.normal);

        vec3 light_from_reflect = raytrace(reflect_ray_origin, reflect_ray_direction, depth+1);

        light = vec3_add( 
            vec3_mulv(light_from_diffuse, vec3_sub((vec3){1.0,1.0,1.0}, reflectance)), 
            vec3_mulv(light_from_reflect, reflectance)
        );
    }
    else
    {
        light = sky_color;
    }

    return light;
}

int main()
{
    const int width = 2048;
    const int height = 1024;

    const int min_res = maxi(width, height);
    const int max_res = maxi(width, height);

    const vec3 world_up = {0.0, 1.0, 0.0};

    const vec3 camera_pos = {5.0, 1.0, 0.0};
    const vec3 camera_look = {0.0, 0.0, 0.0};

    const vec3 camera_dir = vec3_normalize(vec3_sub(camera_look, camera_pos));
    const vec3 camera_right = vec3_normalize(vec3_cross(camera_dir, world_up));
    const vec3 camera_upward = vec3_normalize(vec3_cross(camera_right, camera_dir));

    const float focal_length = 5.0;
    const float aperture = 0.1;

    printf("P3 %i %i 255 ", width, height);

    const int rays_per_pixel = 64;

    for (int y = height - 1; y > 0; y--)
    for (int x = 0; x < width; x++)
    {

        vec3 accumulated = {0.0, 0.0, 0.0};

        for(int i = 0; i < rays_per_pixel; i++)
        {
            // random offset for subpixel antialiasing
            const float rx = randf() - 0.5;
            const float ry = randf() - 0.5;

            const float sx = (float)(2 * (x + rx) - width) / (float)(min_res);
            const float sy = (float)(2 * (y + ry) - height) / (float)(min_res);

            // pos + (dir + right * sx + upward * sy) * focal_length
            const vec3 focus = vec3_add(camera_pos, vec3_mul(vec3_add(vec3_add(camera_dir, vec3_mul(camera_right, sx)), vec3_mul(camera_upward, sy)), focal_length));
            
            const float random_radius = sqrtf(randf()) * aperture;
            const float random_angle = randf() * 2.0 *  M_PI;
            const float x = cos(random_angle) * random_radius;
            const float y = sin(random_angle) * random_radius;

            // right * x + right * y
            vec3 random_aperture_sample = vec3_add(vec3_mul(camera_right, x), vec3_mul(camera_upward, y));

            vec3 origin = vec3_add(camera_pos, random_aperture_sample);

            vec3 direction = vec3_normalize(vec3_sub(focus, origin));

            vec3 color = raytrace(origin, direction, 0);

            accumulated = vec3_add(accumulated, color);
        }

        vec3 average = vec3_div(accumulated, (float)rays_per_pixel);

        int r = clampi((int)(average.x*255.0), 0, 255);
        int g = clampi((int)(average.y*255.0), 0, 255);
        int b = clampi((int)(average.z*255.0), 0, 255);

        printf("%i %i %i ", r, g, b);
    }
}