#version 330 core
out vec4 FragColor; in vec2 vUV; uniform float iTime; uniform vec2 iMouse; uniform vec2 iResolution; uniform float uAlpha;

// Aurora / Boreal shader adapted from Michele Morrone (BorealStarryNight)

#define COMET

float contrast(float valImg, float contrastVal) { return clamp(contrastVal*(valImg-.5)+.5, 0., 1.); }
vec3  contrast(vec3 valImg, float contrastVal)  { return clamp(contrastVal*(valImg-.5)+.5, 0., 1.); }
float gammaCorrection(float imgVal, float gVal)  { return pow(imgVal, 1./gVal); }
vec3  gammaCorrection(vec3 imgVal, float gVal)   { return pow(imgVal, vec3(1./gVal)); }
float cIntensity(vec3 col) { return dot(col, vec3(.299, .587, .114)); }

float hash( float n ) { return fract(sin(n)*758.5453); }
float rand(vec2 p) { return fract(sin(dot(p ,vec2(1552.9898,78.233))) * 43758.5453); }

float noise( in vec3 x )
{
    vec3 p = floor(x);
    vec3 f = fract(x);
    float n = p.x + p.y*57.0 + p.z*800.0;
    float res = mix(mix(mix( hash(n+  0.0), hash(n+  1.0),f.x), mix( hash(n+ 57.0), hash(n+ 58.0),f.x),f.y),
            mix(mix( hash(n+800.0), hash(n+801.0),f.x), mix( hash(n+857.0), hash(n+858.0),f.x),f.y),f.z);
    return res;
}

// simple pseudo-random 3-channel generator used for starfield
vec3 nrand3( vec2 co )
{
   vec3 a = fract( cos( co.x*8.3e-3 + co.y )*vec3(1.3e5, 4.7e5, 2.9e5) );
   vec3 b = fract( sin( co.x*0.3e-3 + co.y )*vec3(8.1e5, 1.0e5, 0.1e5) );
   return mix(a, b, 0.5);
}

#ifdef PULSED_STARS
float starField(vec2 p)
{
   vec3 rnd = nrand3(p * iResolution.x);
   float intensity = pow((1.+sin((iTime+27.)*rnd.x))*.5, 7.) ;
   return max(rnd.x * pow(rnd.y,7.) * intensity, 0.);
}
#else
float starField(vec2 p)
{
   vec3 rnd = nrand3(p * iResolution.x);
   return pow(abs(rnd.x + rnd.y + rnd.z)/2.93,9.7);
}
#endif

float fbm(vec3 p)
{
    float f  = 0.50000*noise( p ); p *= 2.02;
          f += 0.25000*noise( p ); p *= 2.03;
          f += 0.12500*noise( p ); p *= 2.01;
          f += 0.06250*noise( p ); p *= 2.04;
          f += 0.03125*noise( p );
    return f*1.032258;
}

float fbm2(vec3 p)
{
    float f  = 0.50000*noise( p ); p *= 2.021;
          f += 0.25000*noise( p ); p *= 2.027;
          f += 0.12500*noise( p ); p *= 2.01;
          f += 0.06250*noise( p ); p *= 2.03;
          f += 0.03125*noise( p ); p *= 4.01;
          f += 0.015625*noise( p );p *= 8.04;
          f += 0.00753125*noise( p );
    return f*1.05;
}

#define iterations 16
#define formuparam 0.53 //77

#define volsteps 4
#define stepsize 0.00733

#define zoom   1.2700
#define tile   .850
#define speed  0.000

#define brightness 0.0007
#define darkmatter .1700
#define distfading 1.75
#define saturation .250


vec3 starr(vec2 UV)
{
    float ratio = iResolution.y/iResolution.x;
    //get coords and direction
    vec2 uv=UV*.3723+vec2(0.,-.085);
    uv.y*= ratio;
    vec3 dir=vec3(uv*zoom/ratio,1.);

    dir.xz*=mat2(.803, .565, .565, .803);
    dir.xy*=mat2(.9935,.0998,.0998,.9935);

    vec3 from=vec3(-0.4299,-0.7817,-0.3568);

    //volumetric rendering
    float s=0.0902,fade=.7;
    float v=0.;
    for (int r=0; r<volsteps; r++) {
        vec3 p=from+s*dir*-9.9;
        p = abs(vec3(tile)-mod(p,vec3(tile*2.)));
        float pa,a=pa=0.;

        for (int i=0; i<iterations; i++) {
            p=abs(p)/dot(p,p)-formuparam;
            a+=abs(length(p)-pa);
            pa=length(p);
        }
        a*=a*a;
        v+=fade;
        v+=s*a*brightness*fade;
        fade*=distfading;
        s+=stepsize;
    }
    v = contrast(v *.009, .95)-.05;
    vec3 col = vec3(.43,.57,.97) * 1.7 * v;

    col = gammaCorrection(col, .7);

    return col + vec3(.67,.83,.97) * vec3(starField(UV)) * .9;
}

float borealCloud(vec3 p)
{
    p+=fbm(vec3(p.x,p.y,0.0)*0.5)*2.25;
    float a = smoothstep(.0, .9, fbm(p*2.)*2.2-1.1);
    return a<0.0 ? 0.0 : a;
}

vec3 smoothCloud(vec3 c, vec2 pos)
{
    c*=0.75-length(pos-0.5)*0.75;
    float w=length(c);
    c=mix(c*vec3(1.0,1.2,1.6),vec3(w)*vec3(1.,1.2,1.),w*1.25-.25);
    return clamp(c,0.,1.);
}

float shoothingStarLine(vec2 p, vec2 o, vec2 b, float sw, float ew){
    float d = distance(o, b);
    vec2  n = normalize(b - o);
    vec2 l = vec2( max(abs(dot(p - o, n.yx * vec2(-1.0, 1.0))), 0.0),
                   max(abs(dot(p - o, n) - d * 0.5) - d * 0.5, 0.0));
    return smoothstep( mix(sw, ew, 1.-distance(b,p)/d) , 0., l.x+l.y);
}

vec3 comet(vec2 p)
{
    const float modu = 4.;
    const float endPointY = -.1;
    vec2 cmtVel = mod(iTime/modu+modu*.5, 2.) > 1. ? vec2(2., 1.4)*.5 : vec2(-2., 1.4)*.5;
    vec2 cmtLen = vec2(.25)*cmtVel;
    vec2 cmtPt = 1. - mod(iTime*cmtVel, modu);
    cmtPt.x +=1.;
    vec2 cmtStartPt, cmtEndPt;
    if(cmtPt.y < endPointY) {
        cmtEndPt   = cmtPt + cmtLen;
        if(cmtEndPt.y > endPointY) cmtStartPt = vec2(cmtPt.x + cmtLen.x*((endPointY - cmtPt.y)/cmtLen.y), endPointY);
        else                       return vec3(.0);
    }
    else {
        cmtStartPt = cmtPt;
        cmtEndPt = cmtStartPt+cmtLen;
    }

    float bright = clamp(smoothstep(-.2,.65,distance(cmtStartPt, cmtEndPt)),0.,1.);

    vec2 dlt = vec2(.003) * cmtVel;

    float q = clamp( (p.y+.2)*2., 0., 1.);

    return  ( bright  * .75 *  (smoothstep(0.993, 0.999, 1. - length(p - cmtStartPt)) + shoothingStarLine(p, cmtStartPt, cmtStartPt+vec2(.06)*cmtVel,  0.009, 0.003)) +
             vec3(1., .7, .2) * .33 * shoothingStarLine(p, cmtStartPt,         cmtEndPt,        0.003, .0003) +
             vec3(1., .5, .1) * .33 * shoothingStarLine(p, cmtStartPt+dlt,     cmtEndPt+dlt*2., 0.002 ,.0002) +
             vec3(1., .3, .0) * .33 * shoothingStarLine(p, cmtStartPt+dlt+dlt, cmtEndPt+dlt*4., 0.001, .0001)
            ) * (bright) * q;
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 mouseLeft = iMouse / iResolution;
    vec2 mouse = mouseLeft;
    vec2 position = ( fragCoord / iResolution );
    vec2 uv = 2. * position - 1.;
    position.y+=0.2;
    vec2 coord= vec2((position.x-0.5)/position.y,1./(position.y+.2));
    float tm = iTime * .75 * .5;
    coord+=tm*0.0275 + vec2(1. - mouse.x, mouse.y);
    vec2 coord1=coord - tm*0.0275 + vec2(1. - mouse.x, mouse.y);
    vec2 ratio = iResolution.xy / max(iResolution.x, iResolution.y);
    vec2 uvs = uv * ratio;
    vec3 boreal = vec3 (vec3(.1,1.,.5 )  * borealCloud(vec3(coord*vec2(1.2,1.), tm*0.22)) * .9  +
                        vec3(.1,.9,.7) * borealCloud(vec3(coord*vec2(1.,.7)  , tm*0.27)) *  .9 +
                        vec3(.75,.3,.99) * borealCloud(vec3(coord1*vec2(.8,.6)  , tm*0.29)) *  .5 +
                        vec3(.0,.99,.99)  * borealCloud(vec3(coord1*vec2(.9,.5)  , tm*0.20)) *  .57);
    boreal = smoothCloud(boreal, position);
    boreal = gammaCorrection(boreal,1.3);
    // reduce ambient/ground illumination so background remains dark
    boreal *= 0.25;
    float ground = smoothstep(0.07,0.0722,fbm2(vec3((position.x-mouse.x)*5.1, position.y,.0))*(position.y-0.158)*.9+0.01);
    vec3 sunrise = vec3(.6,.3,.99) * (1.0 - position.y > .45 ? (1.0 - position.y - .45) * .65  : 0.) +  vec3(.0,.7,.99) * (1.0 - position.y >= .585 ? (1.0 - position.y - .585  ) * .6  : 0.) + vec3(.5,.99,.99)* (1.0 - position.y >= .67  ? (1.0 - position.y - .67) * .99 : 0.);
    vec3 starColor = starr(uv);
    float ib = clamp(1.0-cIntensity(boreal)*3.,0.,1.);
    ib*=ib;
    float im = 1.0-cIntensity(vec3(.99, .7, .3));
    float ih = 1.0-cIntensity(vec3(.0,  .2,  .7))*8.;
    float is = 1.0-cIntensity(sunrise)*6.;
    vec3 cometA = comet(uvs);
    vec3 sky = vec3(0.);
    sky += ((vec3(.99, .7, .3) + vec3(.0,  .2,  .7))*ib*ib) * is * is + sunrise + cometA + starColor  * ib *ib  * im * im * ih * ih * is * is;
    fragColor =  vec4(boreal + sky, 1.) * ground;
}

void main(){ vec2 fragCoord = gl_FragCoord.xy; mainImage(FragColor, fragCoord); FragColor.a *= uAlpha; }
