float waterline = 1.0;

vec4 hash4( vec4 n ) { return fract(abs(sin(n))*1399763.5453123); }
vec3 hash3( vec3 n ) { return fract(abs(sin(n))*1399763.5453123); }
vec2 NC0=vec2(1.0,136.0);
vec3 NC1=vec3(1.0,136.0,31.0);
float seed = 0.0;

// Simple 2d interpolated noise for mercator projection, using global seed.
float sphereNoise2D(vec2 bl,float scale)
{
    vec2 x = bl;
    vec2 xc = vec2(floor(x.y*scale)) + vec2(1.0,0.0);
    xc = floor(sin(3.1415 * xc/scale)* scale);
    vec4 p = vec4(xc.x*x.x,xc.y*x.x,x.y*scale,x.y*scale+1.0);
    vec4 f = fract(p);
    f=f*f*(3.0-2.0*f);
    p = floor(p);
    vec4 pp = vec4(p.x+1.0>=xc.x? -xc.x:p.x+1.0,p.y+1.0>=xc.y? -xc.y:p.y+1.0,p.z+scale,p.w+scale);
    p.zw +=  vec2(scale);

    vec4 n= hash4(vec4(dot(p.yz,NC0),dot(pp.yz,NC0),dot(p.xw,NC0),dot(pp.xw,NC0))+seed);
    return mix(mix(n.x,n.y,f.y),mix(n.z,n.w,f.x),f.z);
}

// Simple 2d interpolated noise for mercator projection, using local seed.
float sphereNoise2Ds(vec2 bl,float scale,float s)
{
    vec2 x = bl;
    vec2 xc = vec2(floor(x.y*scale)) + vec2(1.0,0.0);
    xc = floor(sin(3.1415 * xc/scale)* scale);
    vec4 p = vec4(xc.x*x.x,xc.y*x.x,x.y*scale,x.y*scale+1.0);
    vec4 f = fract(p);
    f=f*f*(3.0-2.0*f);
    p = floor(p);
    vec4 pp = vec4(p.x+1.0>=xc.x? -xc.x:p.x+1.0,p.y+1.0>=xc.y? -xc.y:p.y+1.0,p.z+scale,p.w+scale);
    p.zw +=  vec2(scale);

    vec4 n= hash4(vec4(dot(p.yz,NC0),dot(pp.yz,NC0),dot(p.xw,NC0),dot(pp.xw,NC0))+s);
    return mix(mix(n.x,n.y,f.y),mix(n.z,n.w,f.x),f.z);
}

// Simple 2d interpolated noise for mercator projection, using global seed animated.
float sphereNoise2D(vec3 blt,float scale)
{
    float z = floor(blt.z)*54.0;
    float f = fract(blt.z);
    //f=f*f*(3.0-2.0*f);
    return mix(sphereNoise2Ds(blt.xy,scale,seed+z),
    sphereNoise2Ds(blt.xy,scale,seed+z+54.0),f);
}

// Simple 2d interpolated distortion noise for mercator projection, using global seed.
float sphereNoise2Do(vec2 bl,float scale,float force)
{
    bl = bl+vec2(sphereNoise2Ds(bl,scale*2.0f,3.0),sphereNoise2Ds(bl,scale*2.0f,20.0))*force/scale;
    bl.x = bl.x>1.0 ? bl.x - 2.0 : bl.x;
    bl.y = clamp(bl.y,0.0,1.0);
    return sphereNoise2D(bl,scale);
}

float sphereNoise2Do2(vec2 bl,float scale,float force, float s)
{
    bl = bl+vec2(sphereNoise2Ds(bl,scale*s,3.0),sphereNoise2Ds(bl,scale*s,20.0))*force/scale;
    bl.x = bl.x>1.0 ? bl.x - 2.0 : bl.x;
    bl.y = clamp(bl.y,0.0,1.0);
    return sphereNoise2D(bl,scale);
}

vec2 sphereIntersect(vec3 ray,vec3 pos,float r)
{
    float r2=r*r;
    float cd = dot(pos,ray);
    vec3 p = ray * cd-pos;
    float t2 = dot(p,p);
    if (t2>=r2) return vec2(-1.0, -1.0);
    float d = sqrt(r2-t2);
    return cd+vec2(- d,d);
}

float sphereIntersectIn(vec3 ray,vec3 pos,float r)
{
    float r2=r*r;
    float cd = dot(pos,ray);
    vec3 p = ray * cd-pos;
    float t2 = dot(p,p);
    if (t2>=r2) return -1.0;
    return cd + sqrt(r2-t2);
}

vec3 sphereNormal(vec3 pos, vec3 surface)
{
    return normalize(surface-pos);
}

float spherePerlin(vec2 bl,float scale)
{
    return (sphereNoise2D(bl,512.0*scale)*0.8+sphereNoise2D(bl,256.0*scale)*0.8
     +sphereNoise2D(bl,128.0*scale)+sphereNoise2D(bl,64.0*scale)
     +sphereNoise2D(bl,32.0*scale)+sphereNoise2D(bl,16.0*scale)
     +sphereNoise2D(bl,8.0*scale)+sphereNoise2D(bl,4.0*scale)
     +sphereNoise2D(bl,2.0*scale)*1.2+sphereNoise2D(bl,scale)*1.2)*0.1;
}

float spherePerlinA(vec3 blt,float scale)
{
    return (sphereNoise2D(blt,16.0*scale)
     +sphereNoise2D(blt,8.0*scale)+sphereNoise2D(blt,4.0*scale)
     +sphereNoise2D(blt,2.0*scale)*1.2+sphereNoise2D(blt,scale)*1.2)*0.2;
}

float spherePerlinHalf(vec2 bl,float scale)
{
    return (sphereNoise2D(bl,16.0*scale)+sphereNoise2D(bl,8.0*scale)+sphereNoise2D(bl,4.0*scale)
     +sphereNoise2D(bl,2.0*scale)+sphereNoise2D(bl,scale))*0.2;
}

// Perlin distortion noise
float spherePerlinHalfd(vec2 bl,float scale,float force)
{
    return (sphereNoise2Do(bl,16.0*scale,force)+sphereNoise2Do(bl,8.0*scale,force)
     +sphereNoise2Do(bl,4.0*scale,force)
     +sphereNoise2Do(bl,2.0*scale,force)+sphereNoise2Do(bl,scale,force))*0.2;
}

// Perlin distortion noise with increasing coefficient
float spherePerlinHalfds(vec2 bl,float scale,float force)
{
    return sphereNoise2Do(bl,16.0*scale,force)*0.05
    +sphereNoise2Do(bl,8.0*scale,force)*0.15+sphereNoise2Do(bl,4.0*scale,force)*0.2
    +sphereNoise2Do(bl,2.0*scale,force)*0.25+sphereNoise2Do(bl,scale,force)*0.35;
}

float spherePerlinHalfds2(vec2 bl,float scale,float force,float s)
{
    return sphereNoise2Do2(bl,16.0*scale,force,16.*s)*0.05
    +sphereNoise2Do2(bl,8.0*scale,force,8.*s)*0.15+sphereNoise2Do2(bl,4.0*scale,force,4.*s)*0.2
    +sphereNoise2Do2(bl,2.0*scale,force,2.*s)*0.25+sphereNoise2Do2(bl,scale,force,s)*0.35;
}

vec2 twirls3D(vec3 pos,float rings,float skip)
{
    vec3 xv = normalize(cross(pos,vec3(0.,1.,0.)));
    vec3 yv = normalize(cross(pos,xv));
    vec3 f = fract(pos)-0.5;
    pos = floor(pos);
    float r = 0.0;
    float s = 0.0;

    for (int z = 0;z<2;z++)
    for (int y = 0;y<2;y++)
    for (int x = 0;x<2;x++) {
        vec3 off = vec3(float(x),float(y),float(z));
        vec3 xyz = pos + off;
        float p = dot(xyz,NC1);
        vec3 add = hash3(vec3(p,p+30.0,p+23.0));
        if (add.z<skip) continue;
        vec3 fxyz = f + add*0.8-off;
        float l = length(fxyz);
        float a = abs(0.5-fract(l*rings*(0.8+add.y*0.2)*2.0+atan(dot(fxyz,xv),dot(fxyz,yv))/3.1415*0.5));//+add.x*10.0

        float al = l*2.0;
        float sf = clamp(1.0-al*al*al,0.,1.);
        r += mix(1.0,a,clamp(l*6.0,0.,1.))*sf;
        s += sf;
    }
    return clamp(vec2(r,s),vec2(0.),vec2(1.));
}

float atmosphere(vec3 pos, vec3 atmosSurface, vec3 ray, vec3 lnorm, float planetRadius, float atmosRadius, float atmosDepth)
{
    float wall = sqrt(atmosRadius*atmosRadius - planetRadius*planetRadius);
    float acc = clamp(atmosDepth/wall,0.0,1.0);
    acc = acc*clamp(dot(lnorm,normalize(atmosSurface-pos)),0.0,1.0);
    acc = acc*acc;
    acc = acc*acc;

    return acc;
}

vec3 sphereTexture(vec3 realSurface,vec3 ray, vec3 snorm,vec3 lnorm,float time)
{
    seed =1.0;
    vec2 bl = vec2(atan(realSurface.x,realSurface.z),atan(length(realSurface.xz),realSurface.y))/3.1415;

    seed = 2.0;
    float d = spherePerlinHalfds(bl,8.0,0.5);
    d = clamp(1.0-abs(d-0.5)*10.0,0.0,1.0);
    seed = 10.0;
    float mc = clamp(spherePerlinHalfds(bl,6.0,0.5)*6.0-2.4, 0.0, 1.0);

    seed = 20.0;
    float c = spherePerlinHalfds(bl,2.0,0.5);
    c = c-d*0.05;

    seed = 5.0;
    float fa = (spherePerlinA(vec3(bl,time),200.0)+spherePerlinA(vec3(bl,time+2.5),200.0))*0.5;

    vec3 water = mix(vec3(0.2,0.0,0.7),vec3(0.0,0.7,1.0),c+mc*0.2);

    float light = clamp(dot(snorm,lnorm),0.0,1.0);
    float ref = clamp(dot(reflect(ray,normalize(snorm+0.6*cross(ray,snorm)*(0.5-fa))),lnorm),0.0,1.0);
    ref = ref*0.9;

    vec3 col = water;

    float l = pow(ref,8.0)*(fa+1.0);
    float wl = clamp((light+pow(ref,6.0)) * mix((abs(0.5-fa)*2.0),0.0,0.3),0.0,1.0);

    vec3 lcolor = l*mix(vec3(1.0),vec3(1.0,0.9,0.8),clamp(dot(ray,lnorm),0.,1.));
    return clamp(mix(col*light+lcolor,vec3(1.0),wl*1.0*light),vec3(0.,0.,0.),vec3(1.,1.,1.));
}

float cloudTexture(vec3 realSurface,float time)
{
   seed = 14.0;
   vec2 bl = vec2(atan(realSurface.x,realSurface.z),atan(length(realSurface.xz),realSurface.y))/3.1415;

   float d = spherePerlinHalfd(bl, 1.0, 1.5);

   float f = clamp((d*2.0-spherePerlinHalfd(bl,20.0,0.6)*0.4-0.53)*4.0,0.0,1.0);
   seed = 24.0;

   float c = spherePerlinA(vec3(bl,time),50.0);
   vec2 t = twirls3D(realSurface*2.,1.,0.);
   t += twirls3D(realSurface*1.,1.5,0.);
   t += twirls3D(realSurface*6.,1.5,0.5)*0.6;
   f = clamp(mix(mix(f,0.0,c),t.x*mix(f,0.0,c),t.y),0.,1.);
   return f;
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    float res = 1.0 / iResolution.y;
  vec2 p = (-iResolution.xy + 2.0*fragCoord.xy) *res;

    float time=iTime*1.0;
    vec3 ray = normalize(vec3(p,2.0));
    vec3 col = vec3(0.0);

  float mx = iMouse.z>0.0 ? iMouse.x/iResolution.x*10.0:0.0;
    float my = iMouse.z>0.0 ? iMouse.y/iResolution.y*4.0-2.0:-0.6;
    vec2 rotate = vec2(mx,my);
    float protate = time*0.025;
    vec2 sins=sin(rotate);
    vec2 coss=cos(rotate);
    vec2 pcs=vec2(cos(protate),sin(protate));
    mat3 mr=mat3(vec3(coss.x,0.0,sins.x),vec3(0.0,1.0,0.0),vec3(-sins.x,0.0,coss.x));
    mr=mat3(vec3(1.0,0.0,0.0),vec3(0.0,coss.y,sins.y),vec3(0.0,-sins.y,coss.y))*mr;
    mat3 pmr=mat3(vec3(pcs.x,0.0,pcs.y),vec3(0.0,1.0,0.0),vec3(-pcs.y,0.0,pcs.x));
    mat3 cmr=mat3(vec3(pcs.x,0.0,pcs.y),vec3(0.0,1.0,0.0),vec3(-pcs.y,0.0,pcs.x));

    float dist=3.0;
    float rplanet = 1.0;
    float ratmos = 1.15;
    float rcloud = 1.005;
    vec3 spos = vec3(0.0,0.0,3.0);
    vec3 lpos = vec3(1.5,1.5,0.0);
    lpos = mr*(lpos-spos) + spos;
    vec3 lnorm = normalize(lpos - spos);

    float dplanet = sphereIntersect(ray,spos,rplanet).x;
    vec2 dcloud = sphereIntersect(ray,spos,rcloud);
    float datmos = sphereIntersect(ray,spos,ratmos).x;
    if (dplanet>0.0) {
        vec3 ssurface = ray*dplanet;
        vec3 snorm = sphereNormal(spos,ssurface);
        vec3 lsnorm = normalize(lpos - ssurface);
        float dcs = sphereIntersectIn(lnorm,spos-ssurface,rcloud);
        vec3 scloud = ssurface+lnorm*dcs;
        float cs = clamp(1.0-cloudTexture((scloud-spos)*mr*pmr*cmr,time*0.2)*2.6,0.0,1.0);

        col = sphereTexture((ssurface-spos)*mr*pmr,ray,snorm,lsnorm,time);
        col = mix(col,col*cs,0.7);
    }
    if (datmos>0.0) {
        float adepth = dplanet>0.0 ? (dplanet-datmos) : abs(dot(ray,spos)-datmos);
        vec3 ssurface = ray*datmos;
        float ap = atmosphere(spos, ssurface, ray, lnorm, rplanet, ratmos, adepth);
        col = mix(col,vec3(0.3,0.76,1.0),ap);
    }
    if (dcloud.x>0.0) {
        vec3 sback = ray*dcloud.y;
        if (dplanet<0.0) {
            vec3 bnorm = sphereNormal(spos,sback);
            col = mix(col, vec3(0.0), cloudTexture((sback-spos)*mr*pmr*cmr,time*0.2) * clamp(dot(bnorm,lnorm),0.0,1.0));
        }
        vec3 ssurface = ray*dcloud.x;
        vec3 snorm = sphereNormal(spos,ssurface);
        vec3 lsnorm = normalize(lpos - ssurface);
        float a = cloudTexture((ssurface-spos)*mr*pmr*cmr,time*0.2);
        float light = clamp(dot(snorm,lsnorm),0.0,1.0);
        col = mix(col+a*light,mix(vec3(0.93,0.9,1.0),vec3(0.94,1.0,1.0),clamp(1.0*a,0.,1.)),clamp(1.5*a*light,0.,1.));
    }
    // Output to screen
    fragColor = vec4(col,1.0);
}