#define STEPS 128.0
#define MDIST 225.0
#define pi 3.1415965
#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))
#define pmod(p,x) (mod(p,x) - (x)*0.5)

float box(vec3 p, vec3 b){
    vec3 d = abs(p)-b;
    return max(d.x,max(d.y,d.z));
}
float smin(float a,float b, float k){ 
    float h=max(0.,k-abs(a-b));
    return min(a,b)-h*h*.25/k;
}

//iq palette
vec3 pal( in float t, in vec3 a, in vec3 b, in vec3 c, in vec3 d ){
    return a + b*cos( 6.28318*(c*t+d) );
}
float sat(float a) {return clamp(a,0.0,1.0);}

float anim(){
    //return 0.0;
    return sat(sin(iTime*0.3)*2.5-1.0);
}
float rand(vec2 a){
    return fract(sin(dot(a,vec2(43.234,21.4343)))*94544.3434343)-0.5;
}
float ssRemap (float t,float s1, float s2, float c){
    return 0.5*(s2-s1)*(t-asin(cos(t*pi)/sqrt(c*c+1.0))/pi)+s1*t;  
}
float ball(vec3 p, float t){
    //Ball
    float scl = 100.0+anim()*50.0; 
    vec3 po = p;
    p.y+=(fract(t)*scl-0.5*scl);
    float a = length(p)-1.0-anim();
    return a;

}
float wave(vec3 p, float t){
    float dist = length(p.xz)-mod(t+0.5,1.0)*55.0+0.5;
    dist = min(dist,0.0);
    float wave = 1.*sin(dist)*exp(-abs(length(dist*0.2)));
    wave*=max(0.0,1.0-fract(t+0.5)*2.0);
    return wave;
}
vec2 map(vec3 p){
    vec2 a = vec2(1);
    vec2 b = vec2(2);
    
    float t = iTime*0.25;
    
    vec3 po2 = p;
    p.y= ssRemap(p.y*0.7,-0.2,0.6,.3);

    p.y=mix(po2.y,p.y,sat(anim()));
    vec3 po = p;
    
    float balls = 10.0;
    float wav = 0.0;
    for(float i = 0.0; i<balls; i++){
        p = po2;
        
        t+=1.0/balls;
        
        float mag = 15.0;
        p.x+=rand(vec2(floor(t),i))*mag;
        p.z+=rand(vec2(floor(t),i*1.5))*mag; 

        //a.x = smin(a.x,ball(p,t),1.0); 
        
        b = vec2(ball(p,t),i+3.0);
        
        a=(a.x<b.x)?a:b;
        
        wav+=wave(p,t)*1.5;
    }
        
    t=iTime*4.5;
    p=po;
    wav+=sin(p.z*0.8+t)*0.1+sin(p.x*0.8+t)*0.1;
    wav+=sin(p.z*0.3-t)*0.2+sin(p.x*0.2-t*2.0)*0.2;
    
    b.x = p.y+wav;
    b.y = 2.0;
    a.x=smin(a.x,p.y+wav,0.4);
    
    a=(a.x<b.x)?a:b;
    
    return a;
}
vec3 norm(vec3 p){
    vec2 e = vec2(0.01,0);
    return normalize(map(p).x-vec3(
    map(p-e.xyy).x,
    map(p-e.yxy).x,
    map(p-e.yyx).x));
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = (fragCoord-0.5*iResolution.xy)/iResolution.y;
    vec3 col = vec3(0);
    
    vec3 ro = vec3(0,15,-30);
    ro.y+=sin(iTime)*3.0;

    ro.xz*=rot(iTime*0.2);

    vec3 lk = vec3(0,0,0);
    vec3 f = normalize(lk-ro);
    vec3 r = normalize(cross(vec3(0,1,0),f));
    vec3 rd = normalize(f*1.0+uv.x*r+uv.y*cross(f,r));
    
    vec3 p = ro;
    float dO, shad;
    vec2 d;
    
    for(float i =0.0; i<STEPS; i++){
        p = ro+rd*dO;
        d = map(p);
        dO+=d.x*0.8;
        
        if((d.x)<0.01){
            shad = i/STEPS;
            break;
        }
        if(dO>MDIST){
            dO=MDIST;
            p = ro+rd*dO;
            d.y=0.0;
            break;
        }
        shad = 1.0;
    }
    vec3 n = norm(p);
    //vec3 ld = vec3(1,1,1);
    vec3 ld=normalize(vec3(1,1.3,2));
    vec3 h = normalize(ld - rd);
    float spec = pow(max(dot(n, h),0.), 20.);
    shad = 1.0-shad;
    
    col = vec3(shad);

    float sss=0.;

    //sss stolen from nusan
    for(float i=1.; i<20.; ++i){
        float dist = i*2.;
        sss += smoothstep(0.,1.,map(p+ld*dist).x/dist);
    }
 
    vec3 al;
    if(d.y==2.0)al = vec3(0.220,0.741,1.000);
    if(d.y>2.0)al = 1.4*pal(d.y/10.0, vec3(0.5),vec3(0.5),vec3(1),vec3(0,0.33,0.66) );
   
    
    
    //I'm sorry for butchering these lighting equations, I have no idea what I am doing
    col = shad*vec3(sss)*al;
    col+=spec*0.3;

    
    //Sky
    p.y-=5.0;
        
    vec3 sky = mix(vec3(0.012,0.471,0.647),vec3(0.145,0.118,0.522),clamp(p.y*0.05,0.0,1.0));
    
    col = mix(col,clamp(sky*3.0,0.0,2.0),dO/MDIST);
        
    fragColor = vec4(col,1.0);
}