// Jacquemet Matthieu

// Contants --------------------------------------------------------------------------------------------
const float PI   = 3.141592653589793238462643383279502884197169;
const float PI_2 = 1.570796326794896619231321691639751442098585;

#define TONEMAP_LINEAR 0
#define TONEMAP_FILMIC 1
#define TONEMAP_REINHARD 2
#define TONEMAP_FILMIC_REINHARD 3
#define TONEMAP_UNCHARTED2 4
#define TONEMAP_ACES 5

// Configuration ----------------------------------------------------------------------------------------

const int MAX_STEPS = 200; // Number of steps
const float EPSILON = 0.01; // Marching epsilon
const float FLOOR_LEVEL = -5.0; // floor height
const float TIME_SCALE = 0.05; // time speed

#define TONEMAP_MODE TONEMAP_ACES // Tone mapping mode
#define _DEBUG 0                  // set 1 to enable debuging



#if _DEBUG
    vec3 _debug_color;
    bool _is_debug = false;
    #define DEBUG(color) if (_is_debug) _debug_color = color; // use this to debug a color
    #define CATCH_DEBUG(expr) _is_debug = true; expr ; _is_debug = false;
#else
    #define DEBUG(color)
    #define CATCH_DEBUG(expr) expr // any call to DEBUG(color) in this macro will output de color
#endif

// TODO: moon shadowing

// Primitive functions -----------------------------------------------------------------------------------

// Hashing function
// Returns a random number in [-1,1]
float Hash(float seed)
{
    return fract(sin(seed)*43758.5453 );
}

// Map value from [Imin, Imax] to [Omin, Omax]
float map(float value, float Imin, float Imax, float Omin, float Omax)
{
  return Omin + (value - Imin) * (Omax - Omin) / (Imax - Imin);
}

// Translate point p
vec3 Translate(vec3 pos, vec3 p) {

    return p - pos;
}

// Scale point p
vec3 Scale(vec3 scale, vec3 p) {

    return vec3(p.x/scale.x, p.y/scale.y, p.z/scale.z);
}

// Rotate point p around X axis (radians)
vec3 RotateX(float theta, vec3 p)
{
    float _sin = sin(theta);
    float _cos = cos(theta);

    mat3 M = mat3(  1,     0,     0,
                    0,  _cos, -_sin,
                    0,  _sin,  _cos);

    return M*p;
}

// Rotate point p around Y axis (radians)
vec3 RotateY(float theta, vec3 p)
{
    float _sin = sin(theta);
    float _cos = cos(theta);

    mat3 M = mat3( _cos,    0,  -_sin,
                    0,      1,      0,
                   _sin,    0,   _cos);

    return M*p;
}

// Rotate point p around Z axis (radians)
vec3 RotateZ(float theta, vec3 p)
{
    float _sin = sin(theta);
    float _cos = cos(theta);

    mat3 M = mat3(  _cos,  -_sin,   0,
                    _sin,   _cos,   0,
                    0,      0,      1);

    return M*p;
}

// Rotate point p
vec3 Rotate(vec3 rot, vec3 p)
{
    p = RotateX(rot.x, p);
    p = RotateY(rot.y, p);
    p = RotateZ(rot.z, p);

    return p;
}


// Create scaling matrix
mat3 Scaling(vec3 scale) {

    return mat3(1.0/scale.x,0,          0,
                0,          1.0/scale.y,0,
                0,          0,          1.0/scale.z);
}

// Create rotation matrix for theta angle around X axis (radians)
mat3 RotationX(float theta)
{
    float _sin = sin(theta);
    float _cos = cos(theta);

    return mat3(    1,     0,     0,
                    0,  _cos, -_sin,
                    0,  _sin,  _cos);
}

// Create rotation matrix for theta angle around Y axis (radians)
mat3 RotationY(float theta)
{
    float _sin = sin(theta);
    float _cos = cos(theta);

    return mat3( _cos,     0, -_sin,
                    0,     1,     0,
                 _sin,     0,  _cos);
}

// Create rotation matrix for theta angle around Z axis (radians)
mat3 RotationZ(float theta)
{
    float _sin = sin(theta);
    float _cos = cos(theta);

    return mat3(_cos, -_sin,  0,
                _sin,  _cos,  0,
                0,     0,     1);
}

// Create rotation matrix for the 3 axes (radians)
mat3 Rotation(float x, float y, float z)
{
    return RotationZ(z) * RotationY(y) * RotationX(x);
}

// Create rotation matrix for the 3 axes (radians)
mat3 Rotation(vec3 rot)
{
    return Rotation(rot.x, rot.y, rot.y);
}

// Create translation matrix
mat4 Translation(vec3 trans) {

    mat4 M = mat4(1.0);
    M[3] = vec4(-trans, 1.0);

    return M;
}

// Create transform matrix
mat4 Transform(vec3 scale, vec3 rot, vec3 trans)
{
    return mat4(Scaling(scale) * Rotation(rot)) * Translation(trans);
}

// Transform of point p
vec3 Transform(vec3 scale, vec3 rot, vec3 trans, vec3 p)
{
    p = Scale(scale, p);
    p = Rotate(rot, p);
    p = Translate(trans, p);

    return p;
}


void Ray(in vec2 m, in vec2 p,out vec3 ro,out vec3 rd)
{
    // focal length
   	float le = 2.0;

    // position camera
    ro=vec3(-40.0,0.0,0.0);

    // reset camera position
    // shadertoy initialize mouse pos at (0,0)
    if (m == vec2(0))
        m = vec2(0.9, 0.51);

    m = (m*2.0 - vec2(1.0))*3.0;
    m.y = clamp(-m.y, -PI_2 + 0.1, PI_2 - 0.1); // clamp camera y rotation

    ro = RotateY(m.y, ro);
    ro = RotateZ(m.x, ro);

    vec3 ww = normalize(-ro);
    vec3 uu = normalize( cross(ww,vec3(0.0,0.0,1.0) ) );
    vec3 vv = normalize( cross(uu,ww));

	rd = normalize( p.x*uu + p.y*vv + le*ww );
}


vec3 Cosine( in float seed, in vec3 n)
{
    float u = Hash( 78.233 + seed);
    float v = Hash( 10.873 + seed);

    // Method by fizzer: http://www.amietia.com/lambertnotangent.html
    float a = 6.2831853 * v;
    u = 2.0 * u - 1.0;
    return normalize( n + vec3(sqrt(1.0-u*u) * vec2(cos(a), sin(a)), u) );
}


// SDF COMBINATION ---------------------------------------------------------------------------

// Union
// a : Field function of left sub-tree,
// b : Field function of right sub-tree
float Union(float a,float b)
{
    return min(a,b);
}

// difference
// a : Field function of left sub-tree,
// b : Field function of right sub-tree
float Diff(float a,float b)
{
    return max(a,-b);
}

// Intersection
// a : Field function of left sub-tree,
// b : Field function of right sub-tree
float Inter(float a, float b)
{
    return max(a,b);
}



// https://www.iquilezles.org/www/articles/distfunctions/distfunctions.html

// Union with smoothing
// a : Field function of left sub-tree,
// b : Field function of right sub-tree
float SmoothUnion( float a, float b, float k )
{
    float h = clamp( 0.5 + 0.5*(b-a)/k, 0.0, 1.0 );
    return mix( b, a, h ) - k*h*(1.0-h);
}

// Difference with smoothing
// a : Field function of left sub-tree,
// b : Field function of right sub-tree
float SmoothDiff( float a, float b, float k )
{
    float h = clamp( 0.5 - 0.5*(a+b)/k, 0.0, 1.0 );
    return mix( a, -b, h ) + k*h*(1.0-h);
}

// Intersection with smoothing
// a : field function of left sub-tree,
// b : field function of right sub-tree
float SmoothInter( float a, float b, float k )
{
    float h = clamp( 0.5 - 0.5*(a-b)/k, 0.0, 1.0 );
    return mix( a, b, h ) + k*h*(1.0-h);
}

// SDF Primitives --------------------------------------------------------------------------

float Plane(vec3 point, vec3 nor, vec3 pos) {

    return dot(point - pos, nor);
}

float Sphere(vec3 point, vec3 pos, float radius) {

    return length(point - pos) - radius;
}

float Circle(vec3 point, float radius) {

    vec2 p;
    p.x = length(point.xy) - radius;
    p.y = point.z;

    return length(p);
}

float Box(vec3 point, vec3 box) {

  vec3 q = abs(point) - box;
  return length(max(q,0.0)) + min(max(q.x,max(q.y,q.z)),0.0);
}


float Torus(vec3 point, float radius, float r) {

    return Circle(point, radius) - r;
}

float Segment(vec3 point, vec3 a, vec3 b) {

    vec3 ba = b-a;
    vec3 pa = point-a;
    float t = dot(pa, ba) / dot(ba, ba);
    vec3 c = ba*clamp(t, 0.0, 1.0);

    return length(pa - c);
}


float Capsule(vec3 point, vec3 a, vec3 b, float radius) {

    return Segment(point, a, b) - radius;
}

float Cylinder(vec3 point, vec3 a, vec3 b, float radius) {

    float v = Diff( Capsule(point, a, b, radius),
                    Plane(point, normalize(a-b), b));

    v = Diff(v, Plane(point, normalize(b-a), a));
    return v;
}


// Returns a smooth noise in [-1,1]
// x : Point
float Noise(in vec3 x) {

    vec3 p = floor(x);
    vec3 f = fract(x);

    p*=1.0;
    f*=1.0;

    f = f*f*(3.-2.*f);
    float n = p.x+p.y*157.+113.*p.z;

    float nn = mix(
        mix(
            mix(Hash(n+0.),   Hash(n+1.0),f.x),
            mix(Hash(n+157.0),Hash(n+158.0),f.x),f.y
        ),
        mix(
            mix(Hash(n+113.0),Hash(n+114.0),f.x),
            mix(Hash(n+270.0),Hash(n+271.0),f.x),f.y
        ),f.z);

    return -1.0 + 2.0*nn;
}


// Create a column
// len : Column's length
// radius : Column's radius
float Column(vec3 p, float len, float radius) {

    const float base_height = 0.2;
    const float top_height = 0.2;
    const float deco0_radius = 0.3;
    const float deco1_radius = 0.2;

    float half_base_height = base_height*0.5;
    float half_top_height = top_height*0.5;

    vec3 base_size = vec3(radius + 0.2, radius + 0.3, half_base_height);
    vec3 top_size = vec3(radius + 0.2, radius + 0.2, half_top_height);

    vec3 base_p = vec3(p.x, p.y, p.z - half_base_height - .15);
    vec3 top_p = vec3(p.x, p.y, p.z - (len-half_top_height));

    vec3 deco0_pos = vec3(0,0,base_height+deco0_radius+0.5);
    vec3 deco1_pos = deco0_pos + vec3(0,0,deco0_radius*0.9+deco1_radius);

    vec3 col_start = vec3(0,0,base_height);
    vec3 col_end = vec3(0,0,len-top_height);

    // column's base block
    float base = Box(base_p, base_size) - 0.1;
    float v = base;


    float offset = pow(abs(p.z-len*0.45)*0.2,2.0)*-0.15 + 0.15;
    float cut = pow(abs(fract(p.z*0.4 + 0.33)-0.5), 15.0)*2000.0;

    float column = Cylinder(p, col_start, col_end, radius+offset-cut+0.05);

    v = SmoothUnion(column, base, 0.7);


    // loower column decoration
    float base0_deco = Torus(p-deco0_pos, radius, deco0_radius);
    float base1_deco = Torus(p-deco1_pos, radius, deco1_radius);

    float base_deco = SmoothUnion(base0_deco, base1_deco,0.1);

    v = SmoothUnion(v, base_deco,0.05);


    // cannelure autour de la colonne
    const float theta = PI/8.0;

    float angle = atan(p.x,p.y);
    angle = fract(angle/theta) - 0.5;
    vec3 cp = vec3(sin(angle)*0.5, length(p.xy) - (radius+0.05)*1.3, p.z);
    cp.y -= offset;

    vec3 start = vec3(0, 0, col_start.z + 2.0*deco0_radius + 1.3);
    vec3 end = vec3(0, 0, col_end.z - 0.6);
    v = SmoothDiff(v, Capsule(cp, start, end, theta*radius*1.2), 0.02);


    float top = Box(top_p, top_size);

    // upper column decoration
    vec3 top_a = vec3( top_size.x,  top_size.x, len-top_height);
    vec3 top_b = vec3(-top_size.x,  top_size.x, len-top_height);

    vec3 tp = p;
    tp.y = abs(tp.y);
    top = Union(top, Cylinder(tp, top_a, top_b, top_height));

    v = Union(v, top);

    return v;
}

// p : Point
// size : Size of the stairs
// num_step : number of stair's steps
float TempleStairs(vec3 p, vec3 size, int num_step) {


    float step_height = size.z / float(num_step);
    float half_step_height = step_height*0.5;
    const float k = 0.1;

    vec3 step_size = size;
    step_size.z = half_step_height - k;

    vec3 pos = p;
    pos.z -= half_step_height;
    float v = Box(pos, step_size) - k;

    for (int i=1; i<num_step; ++i) {
        pos.z -= step_height;
        step_size.xy -= step_height*1.5;
        // one step
        v = Union(v, Box(pos, step_size) - k);
    }
    return v;
}

float TempleRoof(vec3 p, vec3 size) {

    size.z *= 0.5;
    const float k = 0.1;

    float top_size_x = size.x - 0.5;
    vec3 left_normal = normalize(vec3(0,1,size.z));
    vec3 right_normal = left_normal;
    right_normal.y *= -1.0;

    vec3 tp = p - vec3(0,0,1.9);

    // left
    float v = Plane(tp, left_normal, vec3(0,size.y,0));
    // right
    v = SmoothInter(v, Plane(tp, right_normal, vec3(0,-size.y,0)), k);
    // forward
    v = SmoothInter(v, Plane(tp, vec3(1,0,0), vec3(top_size_x,0,0)), k);
    // backward
    v = SmoothInter(v, Plane(tp, vec3(-1,0,0), vec3(-top_size_x,0,0)), k);
    // underneath
    v = SmoothInter(v, Plane(tp, vec3(0,0,-1), vec3(0,0,0)), k);

    vec3 lower_part_size = size;

    const float lp_h = 1.0;
    lower_part_size.z = lp_h;
    float lower_part = Box(p-vec3(0,0,1.12), lower_part_size) - 0.01;

    float cap_z = lower_part_size.z * -0.0;
    float cap_x = lower_part_size.x + lower_part_size.z*2.3;
    float cap_y = lower_part_size.y + lower_part_size.z*2.3;

    vec3 a = vec3( cap_x,  cap_y, cap_z);
    vec3 b = vec3(-cap_x,  cap_y, cap_z);
    vec3 c = vec3( cap_x, -cap_y, cap_z);

    float cap_size = lower_part_size.z*3.0;
    vec3 cap_pos = p;
    cap_pos.xy = abs(cap_pos.xy);
    lower_part = SmoothDiff(lower_part, Capsule(cap_pos, a, b, cap_size),0.1);
    lower_part = SmoothDiff(lower_part, Capsule(cap_pos, a, c, cap_size),0.1);


    return SmoothUnion(v, lower_part-0.1, 0.1);
}

// Potential field of the object -> p : point
float object(vec3 p)
{

    const float stars_height = 2.0;
    const float columns_level = FLOOR_LEVEL + stars_height;
    const float columns_height = 12.0;
    const float columns_radius = 0.7;
    const float roof_level = columns_level + columns_height;

    // destruction
    vec3 destroy_pos = (p + vec3(9,0,0))*0.10 + vec3(4.4,3.7,12.09);
    float destroy_noise = Noise(destroy_pos)*2.0;
    // float destroy_factor = max(-p.z*0.05 + 0.1, 0.0)*10.0;
    // destroy_noise -= destroy_factor;

    // fractal destruction
    for (int i=2; i<5; ++i) {
        float k = float(i);
        destroy_noise += Noise(p*k)/(15.0*k);
    }

    // columns part
    vec3 pc = p;
    pc.xy = fract(pc.xy*0.25) * 4.0 - vec2(2.0,2.0);
    pc.yz -= vec2(0,columns_level);

    float columns = Column(pc, columns_height, columns_radius);

    columns = Inter(columns, Box(p, vec3(11.0+columns_radius, 7.0+columns_radius, columns_height)));
    columns = Diff(columns, Box(p, vec3(9.0-columns_radius, 5.0-columns_radius, columns_height)));

    vec3 fallen_col_p = Rotate(vec3(-PI_2,1.0,-3.2), p);
    fallen_col_p = Translate(vec3(1.0,-2.0,-10.0),fallen_col_p);
    fallen_col_p.z = 13.0 - fallen_col_p.z;

    float fallen_col = Column(fallen_col_p, columns_height, columns_radius);
    columns = Union(columns, fallen_col);

    // roof part
    float roof = TempleRoof(p-vec3(0,0, roof_level), vec3(12.7,8.9,5));
    float v = SmoothInter(Union(roof, columns), destroy_noise, 0.2);


    // stairs part
    vec3 sp = Translate(vec3(0,0,FLOOR_LEVEL), p);
    float stairs = TempleStairs(sp,vec3(14,10,stars_height), 4);

    v = Union(v, stairs);

    return v;
}


// Analysis of the scalar field --------------------------------------------------------------------------

// Calculate object normal
// p : Point
vec3 ObjectNormal(vec3 p)
{
    float eps = 0.002;
    vec3 n;
    float v = object(p);
    n.x = object( vec3(p.x+eps, p.y, p.z) ) - v;
    n.y = object( vec3(p.x, p.y+eps, p.z) ) - v;
    n.z = object( vec3(p.x, p.y, p.z+eps) ) - v;

    return normalize(n);
}

// Trace ray using ray marching
// o : ray origin
// u : ray direction
// e : Maximum distance
// h : hit
// s : Number of steps
float SphereTrace(vec3 o, vec3 u, float e,out bool h,out int s)
{
    h = false;

    // Don't start at the origin, instead move a little bit forward
    float t=0.0;

    for(int i=0; i<MAX_STEPS; i++) {
        s=i;
        vec3 p = o+t*u;
        float v = object(p);
        // Hit object
        if (v < 0.0) {
            s=i;
            h = true;
            break;
        }
        // Move along ray
        t += max(EPSILON,v);
        // Escape marched too far away
        if (t>e)
            break;
    }
    return t;
}

// Lighting --------------------------------------------------------------------

struct PointLight {
    vec3 position;
    vec3 color;
    float energy;
    float radius;
    float shadow_dist;
};

struct DirectionalLight {
    vec3 direction;
    vec3 color;
    float energy;
    float shadow_dist;
};

const int max_light = 16;
int num_lights = 0;

PointLight lights[max_light];

// Add point light to the scene
int addPointLight(vec3 position, vec3 color, float energy, float radius, float shadow_dist)
{
    lights[num_lights] = PointLight(position, color, energy, radius, shadow_dist);
    ++num_lights;

    return num_lights;
}

// Shading ---------------------------------------------------------------------


// Get random vector in [0,1]
vec3 Hash3(vec3 vec, in float seed)
{
    float x = dot(vec, vec3(12.9898,78.233,81.564));
    float y = dot(vec, vec3(54.6515,25.445,13.848));
    float z = dot(vec, vec3(89.5405,95.540,99.244));

    const vec3 seed3 = vec3(43758.5453, 15605.8033, 84053.4468);

    return fract(sin(vec3(x,y,z)*seed3)*seed);
}



float Stars(vec3 viewDir, float grid, float seed)
{
    vec3 scaled = viewDir*grid;
    vec3 local = fract(scaled);
    vec3 cell = floor(scaled);

    // Random vector and scalar for each grid cell
    vec3 randv = Hash3(cell, seed);
    float f = Hash(dot(randv, cell))*0.5 + 0.5;

    float radius = dot(viewDir, normalize(cell+randv));
    radius = sqrt(1.0 - radius*radius) / f*700.0;
    radius = clamp(1.0-radius*radius, 0.0, 1.0);

    return radius * pow(f,2.0)* 0.5;
}


vec3 SkyStars(vec3 viewDir)
{

    const float grid = 50.0; // spacing between stars
    float light = 0.0;

    light += Stars(viewDir, grid, 123456.789);
    light += Stars(RotateY(321.0654,viewDir), grid, 345678.912);
    light += Stars(RotateZ(876.2523, viewDir), grid, 541214.541);

    return vec3(light);
}





// Atmospheric scattering based on preetham's analytical model
// https://github.com/vorg/pragmatic-pbr/blob/master/local_modules/glsl-sky/index.glsl


const float turbidity = 10.0;
const float reileighCoefficient = 2.0;
const float mieCoefficient = 0.005;
const float mieDirectionalG = 0.75;

// constants for atmospheric scattering

const float n = 1.0003; // refractive index of air
const float N = 2.545E25; // number of molecules per unit volume for air at
// 288.15K and 1013mb (sea level -45 celsius)
const float pn = 0.035; // depolatization factor for standard air

// wavelength of used primaries, according to preetham
const vec3 lambda = vec3(680E-9, 550E-9, 450E-9);

// mie stuff
// K coefficient for the primaries
const vec3 K = vec3(0.686, 0.678, 0.666);
const float V = 4.0;

// optical length at zenith for molecules
const float rayleighZenithLength = 8.4E3;
const float mieZenithLength = 1.25E3;
const vec3 up = vec3(0.0, 0.0, 1.0);

const float EE = 1000.0;
const float sunAngularDiameterCos = 0.999956676946448443553574619906976478926848692873900859324;
const float moonAngularDiameterCos = 0.99996192303;
// 66 arc seconds -> degrees, and the cosine of that

// earth shadow hack
const float cutoffAngle = PI/1.95;
const float steepness = 1.5;


vec3 totalRayleigh(vec3 lambda)
{
    return (8.0 * pow(PI, 3.0) * pow(pow(n, 2.0) - 1.0, 2.0) * (6.0 + 3.0 * pn)) / (3.0 * N * pow(lambda, vec3(4.0)) * (6.0 - 7.0 * pn));
}

float rayleighPhase(float cosTheta)
{
    return (3.0 / (16.0*PI)) * (1.0 + pow(cosTheta, 2.0));
    // return (1.0 / (3.0*PI)) * (1.0 + pow(cosTheta, 2.0));
    // return (3.0 / 4.0) * (1.0 + pow(cosTheta, 2.0));
}

vec3 totalMie(vec3 lambda, vec3 K, float T)
{
    float c = (0.2 * T ) * 10E-18;
    return 0.434 * c * PI * pow((2.0 * PI) / lambda, vec3(V - 2.0)) * K;
}

float hgPhase(float cosTheta, float g)
{
    return (1.0 / (4.0*PI)) * ((1.0 - pow(g, 2.0)) / pow(1.0 - 2.0*g*cosTheta + pow(g, 2.0), 1.5));
}

float sunIntensity(float zenithAngleCos)
{
    return max(0.0, 1.0 - exp(-((cutoffAngle - acos(zenithAngleCos))/steepness)));
}



void AtmosphericScattering(DirectionalLight light, vec3 worldNormal,
    out float cosTheta, out vec3 Lin, out vec3 Fex)
{

    vec3 lightDirection = light.direction;
    float lightEnergy = light.energy;

    float sunfade = 1.0-clamp(1.0- exp(light.direction.z / 450000.0) ,0.0,1.0);

    float reileigh = reileighCoefficient - (1.0-sunfade);

    // extinction (absorbtion + out scattering)
    // rayleigh coefficients
    vec3 betaR = totalRayleigh(lambda) * reileigh;

    // mie coefficients
    vec3 betaM = totalMie(lambda, K, turbidity) * mieCoefficient;

    // optical length
    // cutoff angle at 90 to avoid singularity in next formula.
    //float zenithAngle = acos(max(0.0, dot(up, normalize(vWorldPosition - cameraPos))));
    float zenithAngle = acos(max(0.0, dot(up, worldNormal)));
    float sR = rayleighZenithLength / (cos(zenithAngle) + 0.15 * pow(93.885 - ((zenithAngle * 180.0) / PI), -1.253));
    float sM = mieZenithLength / (cos(zenithAngle) + 0.15 * pow(93.885 - ((zenithAngle * 180.0) / PI), -1.253));


    // combined extinction factor
    Fex = exp(-(betaR * sR + betaM * sM));

    // in scattering
    cosTheta = dot(worldNormal, lightDirection);

    float rPhase = rayleighPhase(cosTheta*0.5+0.5);
    vec3 betaRTheta = betaR * rPhase;

    float mPhase = hgPhase(cosTheta, mieDirectionalG);
    vec3 betaMTheta = betaM * mPhase;


    Lin = pow(lightEnergy * ((betaRTheta + betaMTheta) / (betaR + betaM)) * (1.0 - Fex),vec3(1.5));
    Lin *= mix(vec3(1.0),pow(lightEnergy * ((betaRTheta + betaMTheta) / (betaR + betaM)) * Fex,vec3(1.0/2.0)),clamp(pow(1.0-dot(up, lightDirection),5.0),0.0,1.0));
}


vec3 AtmosphericScattering(DirectionalLight sun, DirectionalLight moon, vec3 viewDir) {

    float cosTheta;
    float moonCosTheta;
    vec3 sunLin;
    vec3 moonLin;
    vec3 sunFex;
    vec3 moonFex;

    AtmosphericScattering(sun, viewDir, cosTheta, sunLin, sunFex);
    AtmosphericScattering(moon, viewDir, cosTheta, moonLin, moonFex);


    vec3 texColor = moonLin + sunLin;
    texColor += (sunFex + moonFex)*0.1;
    texColor *= 0.04;
    texColor += vec3(0.0,0.001,0.0025)*0.3;

    return texColor;
}

vec3 GetSkyFex(DirectionalLight light)
{

    float sunfade = 1.0-clamp(1.0-exp(light.direction.z),0.0,1.0);

    float reileigh = reileighCoefficient - (1.0-sunfade);

    // rayleigh coefficients
    vec3 betaR = totalRayleigh(lambda) * reileigh;

    // mie coefficients
    vec3 betaM = totalMie(lambda, K, turbidity) * mieCoefficient;

    // sun optical length
    float zenithAngle = acos(max(0.0, dot(up, light.direction)));
    float sR = rayleighZenithLength / (cos(zenithAngle) + 0.15 * pow(93.885 - ((zenithAngle * 180.0) / PI), -1.253));
    float sM = mieZenithLength / (cos(zenithAngle) + 0.15 * pow(93.885 - ((zenithAngle * 180.0) / PI), -1.253));

    // combined extinction factor
    return exp(-(betaR * sR + betaM * sM));
}


// Get sky ambient color
// sunDirection : Sun direction
// worldNormal : Ray direction
vec3 SkyAmbient(DirectionalLight sun) {

    float cosTheta;
    vec3 Lin;
    vec3 Fex;

    vec3 normal = normalize(sun.direction*1.5 + up);

    AtmosphericScattering(sun, normal, cosTheta, Lin, Fex);

    vec3 L0 = Fex * 0.1;

    vec3 texColor = (Lin+L0) * 0.04;
    texColor += vec3(0.0,0.001,0.0025)*0.3;


    return texColor;
}


// Get sky color
// sunPosition : Sun direction
// worldNormal : Ray direction
vec3 Sky(DirectionalLight sun, DirectionalLight moon, vec3 viewDir, vec3 stallarViewDir) {

    float sunCosTheta;
    float moonCosTheta;
    vec3 sunLin;
    vec3 moonLin;
    vec3 sunFex;
    vec3 moonFex;

    // sun scattering
    AtmosphericScattering(sun, viewDir, sunCosTheta, sunLin, sunFex);

    // moon scattering
    AtmosphericScattering(moon, viewDir, moonCosTheta, moonLin, moonFex);


    //nightsky
    vec3 L0 = SkyStars(stallarViewDir) * moonFex;


    float sundisk = smoothstep(sunAngularDiameterCos,sunAngularDiameterCos+0.00002,sunCosTheta);
    L0 += sun.energy * 19000.0 * sundisk * sunFex;

    float moondisk = smoothstep(moonAngularDiameterCos,moonAngularDiameterCos+0.00002,moonCosTheta);
    L0 += moon.energy * 20.0 * moondisk * moonFex;

    vec3 texColor = (moonLin + sunLin + L0) * 0.04;
    texColor += vec3(0.0,0.001,0.0025)*0.3;


    return texColor;
}




// Cast soft shadow based on https://www.shadertoy.com/view/tlBcDK
// p : Point
// n : Normal at point
// l : Point to light vector
// d : Max tracing distance
float Shadow(vec3 p,vec3 n,vec3 l, float d)
{
    // float v;
    // int s;
    // bool h;
    // SphereTrace(p+0.1*n,l,d,h,s);

    // if (!h)
    //  return 1.0;

    // return 0.0;

    const float k = 20.0;

    float res = 1.0;
    float t = 0.1;

    for (int i = 0; i < MAX_STEPS; ++i) {

        if (res < 0.0 || t > d)
            break;

        float h = object(p+t*l);

        res = min(res, k * h / t);
        t += h;
    }

    return clamp(res, 0.0, 1.0);
}


// Compute ambient occlusion based on https://www.shadertoy.com/view/3lXyWs
// p : Point
// n : Normal at point
float AmbientOcclusion(vec3 p,vec3 n) {

    const int AO_STEPS = 4;
    const float AO_MAX_DIST = 3.0;

    const float SCALE = AO_MAX_DIST / pow(2.0, float(AO_STEPS))*2.0;
    float ocl = 0.0;

    for(int i = 1; i <= AO_STEPS; ++i) {
        float dist = pow(2.0, float(i)) * SCALE;
        ocl += 1.0 - (max(0.0, object(p + n * dist)) / dist);
    }

    return min(1.0-(ocl / float(AO_STEPS)),1.0);
    // return pow(abs(object(p + 2.0 * n)),0.9);
}


// Background color
// r : View ray direction
// cd : Celestial ray direction
vec3 background(vec3 r, vec3 cd, DirectionalLight sun,  DirectionalLight moon)
{
    return Sky(sun, moon, r, cd);

    // return mix(vec3(0.452,0.551,0.995),vec3(0.652,0.697,0.995), d.z*0.5+0.5);
}

// Compute Blinn-Phong specular
// l : Vector to light
// n : Normal at point
// r : View ray direction
// k : glossyness
float Specular(vec3 l, vec3 n, vec3 r, float k)
{
    vec3 half_dir = normalize(l + r);
    float spec_angle = max(dot(half_dir, n), 0.0);
    return pow(spec_angle, k);

    // Phong
//     vec3 ref = reflect(r, n);
//     float c = max(dot(ref, r), 0.0);
//     return pow(c, k/4.0);
}


// Compute point light
// p : Point
// n : Normal at point
// r : View ray direction
// light : light's data
// specular : Specular output
vec3 ComputePointLight(vec3 p,vec3 n, vec3 r, PointLight light, inout vec3 specular) {

    vec3 l = light.position - p;
    float distance = length(l);
    vec3 l_nor = l/distance;

    float lambertian = max(dot(n, l_nor),0.0);

    if (lambertian == 0.0)
        return vec3(0);

    float attenuation =  1.0/(distance*distance + 1.0);
    float s = 1.0;

    if (distance < light.shadow_dist)
        s = Shadow(p, n ,l_nor, distance);

    vec3 c = s * light.energy * light.color * attenuation;

    specular += Specular(l_nor, n, -r, 100.0) * c;

    return lambertian * c;

}

// Compute sun light
// p : Point
// n : Normal at point
// r : View ray direction
// color : color of light
// shadow_dist : Max shadow tracing distance
// specular : Specular output
vec3 ComputeDirectionalLight(vec3 p, vec3 n, vec3 r, DirectionalLight light, inout vec3 specular) {

    float s = Shadow(p, n, light.direction, light.shadow_dist);

    vec3 c = light.color*s;

    specular += Specular(light.direction, n, -r, 100.0) * c;

    return c*max(dot(n , light.direction), 0.0);
}

// Shading and lighting
// p : Point
// n : Normal at point
// r : View ray direction
// sun_dir : Sun direction
vec3 Shade(vec3 p, vec3 n, vec3 r, DirectionalLight sun, DirectionalLight moon)
{
     // Point light
    const vec3 lp = vec3(-20.0, -20.0, 10.0);

    // Light direction to point light
    vec3 l = normalize(lp - p);

    // Ambient color
    vec3 ambient = SkyAmbient(sun) * 0.3;

    // Ambient occlusion
    ambient *= AmbientOcclusion(p, n);

    vec3 surface = vec3(1);

    // Lambert diffuse
    vec3 diffuse = vec3(0);
    vec3 specular = vec3(0);

    // Compute all points lights
    for (int i=0; i<num_lights; ++i) {

        PointLight light = lights[i];

        diffuse += ComputePointLight(p, n, r, light, specular);
    }

    // Compute sun and moon lighting
    if (sun.direction.z > 0.0)
        diffuse += ComputeDirectionalLight(p, n, r, sun, specular);

    if (moon.direction.z > 0.0)
        diffuse += ComputeDirectionalLight(p, n, r, moon, specular);

    return (ambient + diffuse)*surface + specular;
}


// light scattering based on https://ijdykeman.github.io/graphics/simple_fog_shader
// ray_dir : View ray direction
// light_dir : View to light vector
// ray_dist : View ray length
// radius : radius of the light
float Scatter(vec3 ray_dir, vec3 light_dir,float ray_dist, float radius)
{
    const float anisotropy = 0.3;

    float a = dot(-light_dir, ray_dir);
    float c = dot(light_dir, light_dir);
    float h = 1.0/max(sqrt(c-a*a), radius);

    float b = a+ray_dist;

    // approximate mie scattering, but not physically accurate
    h = pow(h, anisotropy + 1.0);

    return (atan(b*h) - atan(a*h))*h;
}

// Gather in-scattering of all lights
// rd : View ray direction
// ro : View ray origin
// t  : View ray length
vec3 Scattering(vec3 rd, vec3 ro, float t) {

    const float density = 0.02;
    vec3 total_scatter = vec3(0);

    for (int i=0; i<num_lights; ++i) {

        PointLight light = lights[i];
        vec3 lp = light.position;
        float radius = light.radius;
        vec3 lv = lp-ro;

        float scatter = Scatter(rd, lv, t, radius);

        total_scatter += scatter*light.color*light.energy;
    }

    return total_scatter * density;
}


// Shading according to the number of steps in sphere tracing
// n : Number of steps
vec3 ShadeSteps(int n)
{
   float t=float(n)/(float(MAX_STEPS-1));
   return 0.5+mix(vec3(0.05,0.05,0.5),vec3(0.65,0.39,0.65),t);
}

// Picture in picture
// pixel : Pixel
// pip : Boolean, true if pixel was in sub-picture zone
vec2 Pip(in vec2 pixel, out bool pip)
{
    // Pixel coordinates
    vec2 p = (-iResolution.xy + 2.0*pixel)/iResolution.y;
   if (pip==true)
   {
    const float fraction=1.0/4.0;
    // Recompute pixel coordinates in sub-picture
    if ((pixel.x<iResolution.x*fraction) && (pixel.y<iResolution.y*fraction))
    {
        p=(-iResolution.xy*fraction + 2.0*pixel)/(iResolution.y*fraction);
        pip=true;
    }
       else
       {
           pip=false;
       }
   }
   return p;
}


// Tone mappin -------------------------------------------------------------------

// based on https://www.shadertoy.com/view/ldcSRN

const float W =11.2; // white scale

// filmic (John Hable)


const float A = 0.22; // shoulder strength
const float B = 0.3; // linear strength
const float C = 0.1; // linear angle
const float D = 0.20; // toe strength
const float E = 0.01; // toe numerator
const float F = 0.30; // toe denominator

vec3 LinearToSRGB(vec3 x)
{
    vec3 t = step(x,vec3(0.0031308));
    return mix(1.055*pow(x, vec3(1./2.4)) - 0.055, 12.92*x, t);
}

vec3 Gamma(vec3 color, float gamma)
{
    return pow(color, vec3(gamma));
}

vec3 Uncharted2Curve(vec3 x)
{
    float A = 0.15;
    float B = 0.50;
    float C = 0.10;
    float D = 0.20;
    float E = 0.02;
    float F = 0.30;

    return ((x*(A*x+C*B)+D*E)/(x*(A*x+B)+D*F))-E/F;
}

vec3 Uncharted2(vec3 color)
{
    vec3 white_scale = Uncharted2Curve(vec3(W));
    return Uncharted2Curve(color) / white_scale;
}


vec3 ReinhardCurve (vec3 x)
{
	return x / (1.0 + x);
}

vec3 Reinhard(vec3 color)
{
    vec3 white_scale = ReinhardCurve(vec3(W));
    return ReinhardCurve(color) / white_scale;
}


vec3 FilmicReinhardCurve (vec3 x)
{
    const float T = 0.01;
    vec3 q = (T + 1.0)*x*x;
	return q / (q + x + T);
}

vec3 FilmicReinhard(vec3 color)
{
    vec3 white_scale = FilmicReinhardCurve(vec3(W));
    return FilmicReinhardCurve(color) / white_scale;
}


vec3 FilmicCurve(vec3 x)
{
	return ((x*(0.22*x+0.1*0.3)+0.2*0.01)/(x*(0.22*x+0.3)+0.2*0.3))-0.01/0.3;
}

vec3 Filmic(vec3 color)
{
    vec3 white_scale = FilmicCurve(vec3(W));
    return FilmicCurve(color) / white_scale;
}


vec3 ACESFitted(vec3 color) {

    color = pow(color, vec3(0.833));
    color *= 1.07;

    const mat3 ACESInput = mat3(
        0.59719, 0.35458, 0.04823,
        0.07600, 0.90834, 0.01566,
        0.02840, 0.13383, 0.83777
    );

    const mat3 ACESOutput = mat3(
        1.60475, -0.53108, -0.07367,
        -0.10208,  1.10813, -0.00605,
        -0.00327, -0.07276,  1.07602
    );


    color = color * ACESInput;

    // Apply RRT and ODT
    vec3 a = color * (color + 0.0245786) - 0.000090537;
    vec3 b = color * (0.983729 * color + 0.4329510) + 0.38081;
    color = a/b;

    return color * ACESOutput;
}


vec3 ToneMapping(vec3 color) {

    color = Gamma(color,1.0);

    #if TONEMAP_MODE == TONEMAP_FILMIC
        color = Filmic(color);
    #elif TONEMAP_MODE == TONEMAP_REINHARD
        color = Reinhard(color);
    #elif TONEMAP_MODE == TONEMAP_FILMIC_REINHARD
        color = FilmicReinhard(color);
    #elif TONEMAP_MODE == TONEMAP_UNCHARTED2
        color = Uncharted2(color);
    #elif TONEMAP_MODE == TONEMAP_ACES
        color = ACESFitted(color);
    #endif

    color = clamp(LinearToSRGB(color), 0.0, 1.0);

    return color;
}


float smoothRand(float seed) {

    float amplitude = 1.;
    float frequency = 1.;
    float y = sin((iTime+seed) * frequency);
    float t = 0.01*(-(iTime+seed)*130.0);
    y += sin(frequency*2.1 + t)*4.5;
    y += sin(frequency*1.72 + t*1.121)*4.0;
    y += sin(frequency*2.221 + t*0.437)*5.0;
    y += sin(frequency*3.1122+ t*4.269)*2.5;
    y *= amplitude*0.06;

    return y;
}

// Image
void mainImage( out vec4 color, in vec2 pxy )
{
    // central blue light
    float x, y = 0.0;
    float z = 2.0;

    float r = 0.0;
    float g = 0.5;
    float b = 1.0;

    addPointLight(vec3(x,y,z), vec3(r,g,b), 1.0, 0.3, 0.0);

    // Picture in picture on
    bool pip=true;

    // Pixel
    vec2 pixel=Pip(pxy, pip);

    // Mouse
    vec2 m=iMouse.xy/iResolution.xy;

    // Camera
    vec3 ro,rd;
    Ray(m,pixel,ro,rd);


    // Hit and number of steps
    bool hit;
    int s;


    CATCH_DEBUG(float t = SphereTrace(ro, rd, 1000.0, hit, s));

    const float latitude = 0.0;
    const float orbital_axis = 0.41015;

    // Position
    vec3 pt = ro + t * rd;

    DirectionalLight sun;
    DirectionalLight moon;

    const float sunAzimuth = 1.5;
    const float moonAzimuth = 1.0;

    float curr_time = iTime*TIME_SCALE + PI;

    // curr_time = 6.7; // day
    // curr_time = 9.3; // night
    // curr_time = 1.6; // sunset

    float earthRotTheta = curr_time;

    mat3 celestial_rot = RotationX(earthRotTheta);
    celestial_rot = RotationZ(earthRotTheta) * celestial_rot;
    celestial_rot = RotationX(orbital_axis) * celestial_rot;

    const vec3 sun_pos = vec3(cos(sunAzimuth), sin(sunAzimuth),0);

    float moonRotTheta = curr_time * 0.0338983;
    mat3 moon_rot = RotationZ(moonRotTheta); // moon around earth
    moon_rot *= RotationX(0.08970992); // moon orbit
    vec3 moon_pos = moon_rot[0];

    sun.direction = celestial_rot * sun_pos;
    moon.direction = celestial_rot * moon_pos;

    sun.energy = sunIntensity(sun.direction.z) * EE;
    moon.energy = sunIntensity(moon.direction.z) * 5.0;

    float elevation = acos(sun.direction.z) / PI_2;

    vec3 rgb;

    moon.color = GetSkyFex(moon)*0.005;
    sun.color = GetSkyFex(sun)*19.0;

    // Some optimization for preventing useless shadow computation
    if (sun.direction.z < 0.0) {
        sun.shadow_dist = 0.0;
        moon.shadow_dist = 100.0;
    } else {
        sun.shadow_dist = 100.0;
        moon.shadow_dist = 0.0;
    }


    if (hit || rd.z < 0.0) {
        vec3 n;

        if (hit)
            n = ObjectNormal(pt);

        else { // infinite ground plane
            n = up;

            float co = rd.z;
            float si = sqrt(1.0-co*co);
            float ta = si/co;

            float a = ro.z - FLOOR_LEVEL;
            float b = ta*a;

            t = length(vec2(a,b));
            pt = ro + t * rd;
        }

        // Shade object with light
        rgb = Shade(pt, n, rd, sun, moon);

        vec3 atmosphere = AtmosphericScattering(sun, moon, rd);
        float fog = clamp(smoothstep(40.0,1000.0,t), 0.0, 1.0);
        rgb = mix(rgb, atmosphere, fog);
    }
    else {
        // Shade background
        CATCH_DEBUG(rgb = background(rd, rd*celestial_rot, sun, moon));
    }

    // Volumetric lighting
    rgb += Scattering(rd, ro, t);


    // Auto exposure
    float sun_ext = clamp(map(elevation, 0.98, 1.2, 0.0, 1.0), 0.0, 1.0);
    sun_ext = smoothstep(0.0, 1.0, sun_ext);
    float exposure_bias = mix(0.2, 20.0, sun_ext);
    rgb *= exposure_bias;

    rgb = ToneMapping(rgb);

    // Uncomment this line to shade image with false colors representing the number of steps
    if (pip==true)
        rgb = ShadeSteps(s);

#if _DEBUG
    rgb = _debug_color;
#endif

    color = vec4(rgb, 1.0);
}