/*original https://www.shadertoy.com/view/tllfRX original https://www.shadertoy.com/view/dl3XDS  https://www.shadertoy.com/view/ldySDh*/

#define NUM_LAYERS 8.
#define TAU 6.28318
#define PI 3.141592
#define Velocity .025 //modified value to increse or decrease speed, negative value travel backwards
#define StarGlow 0.525
#define StarSize 02.
#define CanvasView 20.


#define MULTISAMPLES 3 // Max 4
#define R(p,a,r)mix(a*dot(p,a),p,cos(r))+sin(r)*cross(p,a)
#define H(h)(cos((h)*6.3+vec3(0,23,21))*.5+.5)
const float N = 1.33;
const float zoom = 2.0;
const int max_intersections = 12;

const float eyedistance = 7.5; // Note: These depend on each other
const float min_distance = 3.0;
const float max_distance = 10.5;
const float min_stepsize = 0.25;
const int maxsteps = 30;

const float pi = 3.1415926536;

vec4 sphere1;
vec4 sphere2;
vec4 sphere3;

float sq(float x) { return x * x; }
float sq(vec3 x) { return dot(x, x); }

float fresnel(float n1, float n2, float cos_theta)
{
    float r = sq((n1 - n2) / (n1 + n2));
    return r + (1.0 - r) * pow(1.0 - clamp(cos_theta, 0.0, 1.0), 5.0);
}



float f(vec3 p)
{
    return 1.0 - (
        sphere1.w / sq(sphere1.xyz - p) +
        sphere2.w / sq(sphere2.xyz - p) +
        sphere3.w / sq(sphere3.xyz - p));
}

vec3 fd(vec3 p)
{
    vec3 d1 = sphere1.xyz - p;
    vec3 d2 = sphere2.xyz - p;
    vec3 d3 = sphere3.xyz - p;
    return 2.0 * (
        sphere1.w * d1 / sq(sq(d1)) +
        sphere2.w * d2 / sq(sq(d2)) +
        sphere3.w * d3 / sq(sq(d3)));
}

float stepsize(vec3 p)
{
    float md = sqrt(min(min(
        sq(p - sphere1.xyz),
        sq(p - sphere2.xyz)),
        sq(p - sphere3.xyz)));
    return max(min_stepsize, abs(md - 1.0) * 0.667);
}

vec4 ray(vec3 p, vec3 d)
{
    float k = min_distance;
    float nf = 1.0;
    vec4 c = vec4(0.0);
    float cr = 1.0;
    for (int j = 0; j < max_intersections; ++j)
    {
        for (int i = 0; i < maxsteps; ++i)
        {
            if (k > max_distance)
                return c  * cr;
            float ss = stepsize(p + d * k);
            if (f(p + d * (k + ss)) * nf < 0.0)
            {
                k += ss - min_stepsize * 0.5;
                k += f(p + d * k) / dot(d, fd(p + d * k));
                k += f(p + d * k) / dot(d, fd(p + d * k));
                p += d * k;

                vec3 n = -normalize(fd(p)) * nf;
                vec3 r = refract(d, n, nf > 0.0 ? 1.0 / N : N);

                if (nf < 0.0)
                {
                    float fa = k * 0.025;
                    c += vec4(0.5, 0.75, 1.0, 1.0) * fa * cr;
                    cr *= 1.0 - fa;
                }

                if (r == vec3(0.0))
                {
	                d = reflect(d, n);
                }
                else
                {
                    float f = nf > 0.0 ?
                        fresnel(1.0, N, dot(-d, n)) :
                    	fresnel(N, 1.0, dot(-d, n));
                    if (f > 0.5)
                    {
                        c +=   (1.0 - f) * cr;
                        cr *= f;
                        d = reflect(d, n);
                    }
                    else
                    {
                        c +=  f * cr;
                        cr *= 1.0 - f;
                        d = r;
                        nf *= -1.0;
                    }
                }
                k = 0.0;
                break;
            }
            k += ss;
        }
    }
    return c  * cr;
}

float Star(vec2 uv, float flare){
    float d = length(uv);
  	float m = sin(StarGlow*1.2)/d;
    float rays = max(0., .5-abs(uv.x*uv.y*1000.));
    m += (rays*flare)*2.;
    m *= smoothstep(1., .1, d);
    return m;
}

float Hash21(vec2 p){
    p = fract(p*vec2(123.34, 456.21));
    p += dot(p, p+45.32);
    return fract(p.x*p.y);
}
#define Rot(a) mat2(cos(a),-sin(a),sin(a),cos(a))

vec3 StarLayer(vec2 uv){
    vec3 col = vec3(0);
    vec2 gv = fract(uv);
    vec2 id = floor(uv);
    for(int y=-1;y<=1;y++){
        for(int x=-1; x<=1; x++){
            vec2 offs = vec2(x,y);
            float n = Hash21(id+offs);
            float size = fract(n);
            float star = Star(gv-offs-vec2(n, fract(n*34.))+.5, smoothstep(.1,.9,size)*.46);
            vec3 color = sin(vec3(.2,.3,.9)*fract(n*2345.2)*TAU)*.25+.75;
            color = color*vec3(.9,.59,.9+size);
            star *= sin(iTime*.6+n*TAU)*.5+.5;
            col += star*size*color;
        }
    }
    return col;
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = (fragCoord-.5*iResolution.xy)/iResolution.y;
	vec2 M = vec2(0);
        float t = iTime;

    vec4 vs1 = cos(t * vec4(0.87, 1.13, 1.2, 1.0) + vec4(0.0, 3.32, 0.97, 2.85)) * vec4(-1.7, 2.1, 2.37, -1.9);
    vec4 vs2 = cos(t * vec4(1.07, 0.93, 1.1, 0.81) + vec4(0.3, 3.02, 1.15, 2.97)) * vec4(1.77, -1.81, 1.47, 1.9);
    vec4 O = fragColor;
    vec2 C = fragCoord;
 O=vec4(0);
    vec3 p,q,r=iResolution,
    d=normalize(vec3((C*2.-r.xy)/r.y,1));

    for(float i=0.,a,s,e,g=0.;
        ++i<110.;
        O.xyz+=mix(vec3(1),H(g*.1),sin(.8))*1./e/8e3
    )
    {



        p=g*d;
        float t = mod(iTime*0.5,4.);

    if(t<1.){
       p.z+=tan(iTime*0.5);
       p.yx*=Rot(-iTime);

    } else if (t>=1. && t<2.){
       p.x+=tan(iTime*0.5);
       p.xz*=Rot(-iTime);
    } else if (t>=2. && t<3.){
         p.y+=tan(iTime*0.5);
          p.xz*=Rot(iTime);
    } else if (t>=3. && t<4.){
       p.y+=tan(iTime*0.5);
   p.yz*=Rot(-iTime);

        }


         float t2 = mod(iTime*0.7,4.);
          if(t2<1.){
          a=17.;
        p=mod(p-a,a*2.)-a;
        s=2.;

    } else if (t2>=1. && t2<2.){

        s=3.;
         a=15.;
       }

        else if (t2>=2. && t2<3.){

        s=4.;
        a=14.;

        } else if (t2>=3. && t2<4.){

        s=5.;
        a=11.;
       }



        for(int i=0;i++<8;){
            p=.3-abs(p);

            p.x<p.z?p=p.zyx:p;
            p.z<p.y?p=p.xzy:p;

            s*=e=1.4+sin(iTime*.234)*.1;
            p=abs(p)*e-
                vec3(
                    5.+cos(iTime*.3+.5*(iTime*.3))*3.,
                    120,
                    8.+cos(iTime*.5)*5.
                 );
         }
         g+=e=length(p.yz)/s;
    }
    sphere1 = vec4(vs1.x, 0.0, vs1.y, 1.0);
	sphere2 = vec4(vs1.z, vs1.w, vs2.z, 0.9);
	sphere3 = vec4(vs2.x, vs2.y, vs2.w, 0.8);

    vec2 r2 = -iMouse.yx / iResolution.yx * pi * 2.0;

    vec4 cs = cos(vec4(r2.y, r2.x, r2.y - pi * 0.5, r2.x - pi * 0.5));
    vec3 forward = -vec3(cs.x * cs.y, cs.w, cs.z * cs.y);
	vec3 up = vec3(cs.x * cs.w, -cs.y, cs.z * cs.w);
	vec3 left = cross(up, forward);
    vec3 eye = -forward * eyedistance;


    vec3 dir = normalize(vec3(forward + uv.y * up + uv.x * left));
    vec4 color = ray(eye, dir);
#if MULTISAMPLES > 1
    vec2 uvh = zoom * vec2(0.5) / iResolution.x;
    color += ray(eye, normalize(forward + (uv.y + uvh.y) * up + (uv.x + uvh.x) * left));
#if MULTISAMPLES > 2
    color += ray(eye, normalize(forward + (uv.y + uvh.y) * up  + uv.x * left));
#if MULTISAMPLES > 3
    color += ray(eye, normalize(forward + uv.y * up + (uv.x + uvh.x) * left));
#endif
#endif
    color /= float(MULTISAMPLES);
#endif
    M -= vec2(M.x+sin(iTime*0.22), M.y-cos(iTime*0.22));
    M +=(iMouse.xy-iResolution.xy*.5)/iResolution.y;
    float t2 = iTime*Velocity;
    vec3 col = vec3(0);
    for(float i=0.; i<1.; i+=1./NUM_LAYERS){
        float depth = fract(i+t2);
        float scale = mix(CanvasView, .5, depth);
        float fade = depth*smoothstep(1.,.9,depth);
        col += StarLayer(uv*scale+i*453.2-iTime*.05+M)*color.xyz+O.xyz*fade;}
    fragColor = vec4(col,1.0);
}