#define SHADOWS
//#define AA

#define STEPS 300
#define FAR 60.0

// Number of spheres
const int n = 61;

const float nf = float(n);
const float radius = 3./nf;
vec4 sphere[n];

const float firstImage = 8.;
const float speed = 2.;

const vec3 fog = vec3(.1,.1,.15);
vec3 moonDirection = normalize(vec3(1,1,-2));
vec3 moonCol = vec3(.1,.1,.15);
const float fov = 35.0;
const float fovFactor = tan(0.5*fov*0.01745);

float time;
const float PI = 3.14159;

// Hash functions from Dave_Hoskins
// https://www.shadertoy.com/view/4djSRW
float hash11(float p)
{
    p = fract(p * .1031);
    p *= p + 33.33;
    p *= p + p;
    return fract(p);
}
float hash12(vec2 p) {
	vec3 p3  = fract(vec3(p.xyx) * .1031);
    p3 += dot(p3, p3.yzx + 33.33);
    return fract((p3.x + p3.y) * p3.z);
}
float noise(vec2 p) {
    vec2 fl = floor(p);
 	vec2 fr = p - fl;
    fr = fr*fr*(3.-2.*fr);
    
    const vec2 b = vec2(0,1);
    return mix(mix(hash12(fl+b.xx), hash12(fl+b.yx), fr.x),
               mix(hash12(fl+b.xy), hash12(fl+b.yy), fr.x), fr.y);
}
float fbm(vec2 p) {
    float h = 0., f=1., a=.5;
    h += a*noise(f*p); f *= 2.; a *=.5;
    h += a*noise(f*p); f *= 2.; a *=.5;
    h += a*noise(f*p);
    return h;
}
// Smoothmin taken from vgs
// https://www.shadertoy.com/view/Ml3Gz8
// Polynomial smooth min (for copying and pasting into your shaders)
float smin(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);
}

mat3 viewMatrix(vec3 forward, vec3 up) {
 	vec3 w = -normalize(forward);
    vec3 u = normalize(cross(up, w));
    vec3 v = cross(w, u);
    
    return mat3(u,v,w);
}

// Camera path
vec3 camPosition(float t) {
    t *= speed;
    return vec3(1.*cos(t),.2*sin(2.*t),-3.*t);
}

vec3 camDirection(float t) {
    return normalize(camPosition(t+.1)-camPosition(t)+vec3(0,.12,0));
}
mat3 camMatrix(float t) {
    vec3 f = camDirection(t);
    vec3 u = vec3(0,1,0);
    return viewMatrix(f,u);
}

// We define here the images.
// 2D position of sphere number *i* on image number *im*
vec2 spherePositionOnImage(float i, float im) {
    return .8*(.0+1.*i/nf)*vec2(cos(2.*PI*i*im/nf),sin(2.*PI*i*im/nf));
}
vec3 sphereColor(float i) {
    return .5+.5*sin(.02*vec3(4,5,6)*i+vec3(3,4,5));
}

// Time at which image number *im* is seen on the screen  
float meetingTime(float im) {
    return 3.*PI*(im-firstImage)/speed;
}

// Calculates 3D sphere positions so that they project to the 
// 2D positions defined in spherePositionOnImage.
// Returns 3D position in .xyz channels, radius in .a channel.
vec4 spherePosition(float i, float im) {
    vec2 q = spherePositionOnImage(i, im);
    // Random depth
    float d = 5.+.2*i+10.*hash11(i+im);
    
    float r = fovFactor*radius*d;
    
    float t = meetingTime(im);
    vec3 eye = camPosition(t);
    mat3 mEye = camMatrix(t);
    
    return vec4(eye + d*mEye*vec3(fovFactor*q,-1.), r);
}

vec4 interpolate(vec4 a, vec4 b, float t) {
    return (1.-t)*(1.-t)*(1.-t)*a +
        3.*(1.-t)*(1.-t)*t*vec4(-a.x,-a.y,a.z,0)+ //.1*a.xy,a.z-10.,a.w) +
        3.*(1.-t)*t*t*vec4(5.*b.xy,b.z+10.,b.w) +
        t*t*t*b;
}

// Calculates 3D sphere positions by interpolating between
// the positions computed by spherePosition(i, im)
vec4 spherePosition(float i) {
    float fr = speed*time/(3.*PI)+firstImage;
    float image = floor(fr);
    fr -= image;
    
    fr = smoothstep(.0,.5,fr);
    
    vec4 a = spherePosition(i, image);
    vec4 b = spherePosition(i, image+1.);
    
    return interpolate(a,b,fr);
}

// Calculate light coming from spheres along a ray.
// Discard spheres farther than given depth.

// Returns color in .rgb channels, light quantity in .a channel
vec4 spheresLight(vec3 start, vec3 dir, float depth) {
    vec4 col = vec4(0);
    float minSphereDepth = 1e6;
        for(int i=0; i<n; i++) {
            // Compute sphere depth and distance to the ray
            vec4 q = sphere[i];
            q.xyz -= start;
            float sphereDepth = dot(q.xyz,dir);
            float d = sphereDepth * length(normalize(q.xyz)-dir);
            d -= q.a;
            
            if(sphereDepth<depth) {// Discard if sphere behind obstacle
                if(d<0. && sphereDepth<minSphereDepth) {
                    // Closest sphere 
                    col.xyz = 4.*sphereColor(float(i));
                    minSphereDepth = sphereDepth;
                    col.a = 1.;
                } else {
                    // Glow
                    float k = pow(max(0.,1.-2.*d),3.);
                    col.xyz += k*sphereColor(float(i));
                    col.a += k;

                }
            }
        }
    return col;
}

// Snow bumps
float bump(vec2 p) {
    return fbm(3.*p);
}
// Sd to floor
float sd(vec3 p) {
    // A tunnel arond camera path...
    float t = -p.z/3.;
    float dx = abs(p.x-1.*cos(t));
    float h = .2*sin(2.*t)+.1*dx*dx-.1;
    // ...with bumps
    h += .2*bump(p.xz)-.1;
   
    h = smin(h,5.,1.);
    
    return .7*(p.y-h);
}

vec3 normal(vec3 p) {
    const vec2 delta = vec2(0,.001);
    float d = sd(p);
    return normalize(vec3(
        sd(p+delta.yxx),
        sd(p+delta.xyx),
        sd(p+delta.xxy))-d);
}

float march(vec3 start, vec3 dir) {
	float total = 0., d = 1.;
    float epsilon = fovFactor*.5/iResolution.y;
    int i=0;
    for(; i<STEPS; i++) {
        if(d<epsilon*total) return total;
        if(total>FAR) break;
        d = sd(start + total*dir);
        total += d;
    }
    return FAR+1.;
}

float lightMarch(vec3 start, vec3 dir, float lightDist) {
	float total = 0., d = 1.;
    float epsilon = fovFactor*.5/iResolution.y;
    float minD = 1000.0;
    for(int i=0; i<STEPS; i++) {
        if(d<epsilon*total || total>lightDist) break;
        d = sd(start + total*dir);
        total += d;
        minD = min(minD,d);
    }
    return minD;
}
float shadow(vec3 p, vec3 toLight, float lightDist) {
    float minD = lightMarch(p+0.05*toLight, toLight, lightDist-.05);
    return smoothstep(0.0,0.04,minD);
}
vec3 rayColor(vec3 start, vec3 dir) {
    vec3 color = vec3(0);
    
    float minD;
    float d = march(start, dir);

    if(d > FAR) {// SKY
        // Sky
        vec3 sky = vec3(0,0,.05*pow(1.-dir.y,2.));
        sky = mix(fog,sky,dir.y/(.01+dir.y));
        
        // Stars
        const vec3 star = vec3(1,1,.7);
        vec2 q = 40.*dir.xy;
        vec2 fl = floor(q);
        fl += .1+.8*hash12(fl);
        // No star with probability .9
        float bit = step(.9,hash12(fl+315.812));
        float r = min(.1*hash12(fl+153.81),bit);
        color = mix(star,sky,smoothstep(.5*r,r,length(q-fl)));
        
        // Moon
        float l = length(dir-moonDirection);
        float s = 1./(1.+1e5*pow(l,4.));
        // Moon glow
        color = mix(color,2.*moonCol,s);
        // Moon
        color = mix(30.*moonCol*fbm(50.*dir.xy),color,smoothstep(.048,.05,l));
    } else {
        vec3 p = start + d * dir;

        vec3 normal = normal(p);
                
        // Fake subsurface: more light at the top of bumps
        float b = .5-bump(p.xz);
        const vec3 ice = vec3(.3,.8,1);
        vec3 subsurface = ice*exp(-3.*b);
        {
            // Diffuse
            float diff = max(dot(moonDirection, normal), 0.0);

            // Specular
            vec3 h = normalize(moonDirection-dir);
            float spec = 20. * pow(max(dot(h,normal),0.0), 300.0);

            #ifdef SHADOWS
            // Shadow
            float sh = shadow(p, moonDirection, 10.);
            #else
            float sh = 1.;
            #endif
            color += sh*(subsurface + diff + spec)*moonCol;
        }
        for(int i=0; i<n; i++) {
            vec3 lightPosition = sphere[i].xyz;
            vec3 lightCol = 30.*sphereColor(float(i))*sphere[i].a;
            
            vec3 lightDirection = lightPosition-p;
            float lightDist = length(lightDirection);
            lightDirection /= lightDist;


            // Diffuse
            float diff = max(dot(lightDirection, normal), 0.0);
            
            // Specular
            vec3 h = normalize(lightDirection-dir);
            float spec = 20. * pow(max(dot(h,normal),0.0), 300.0);

            color += (subsurface + diff + spec)*lightCol/(1.+20.*lightDist*lightDist);
        }
       // color = subsurface;
        color = mix(fog,color,smoothstep(100.,90.,minD));//exp(-.0003*d*d));
    }

    vec4 l = spheresLight(start, dir, d);
    return mix(color, l.xyz, l.a);
}

void initSpheres() {
    for(int i=0; i<n; i++) {
        sphere[i] = spherePosition(float(i));
    }
}
void mainImage( out vec4 fragColor, in vec2 fragCoord ) {
    time = iTime;//+3.*PI;
    vec3 cam = camPosition(time);
        mat3 m = camMatrix(time);

    if(iMouse.z>0.0) {
        float a = -1.5*PI*(2.*iMouse.x/iResolution.x-1.);
        float b = .5*PI-.53*PI*iMouse.y/iResolution.y;
        
        vec3 center = vec3(0,3,cam.z-23.);
        vec3 relativeCam = 40.*vec3(sin(a)*cos(b),sin(b),cos(a)*cos(b));
        cam = center + relativeCam;
        
        m = viewMatrix(-relativeCam,vec3(0,1,0));
       // time = 100.*iMouse.x/iResolution.x;
    }
    //moonDirection = normalize(vec3(cos(.1*iTime),sin(.1*iTime),-1.5));
    initSpheres();
    
    vec3 color = vec3(0.0);

    vec2 uv;
    #ifdef AA
    for(float i=-0.25; i<0.5; i+=0.5) {
        for(float j=-0.25; j<0.5; j+=0.5) {
            uv = (2.0*(fragCoord + vec2(i,j)) - iResolution.xy)/iResolution.y;
            vec3 pix = vec3(tan(0.5*fov*0.01745)*uv,-1.0);
    
            vec3 dir = normalize(m*pix);
    
            color += rayColor(cam, dir);
        }
    }
    color /= 4.0;
    #else
    uv = (2.*fragCoord - iResolution.xy)/iResolution.y;
    vec3 pix = vec3(fovFactor*uv,-1.0);
    vec3 dir = normalize(m*pix);

    color = rayColor(cam, dir);
    #endif
        
  //  color = color/(1.+color);
    color = 1.-exp(-color);
    // Vignette
    uv = fragCoord.xy / iResolution.xy;
    uv *=  1. - uv.yx;
    color *= pow(uv.x*uv.y * 15.0, 0.25);
    
    // Gamma
    color = pow(color, vec3(1./2.2));
        
    fragColor = vec4(color,1.0);
}