const float ArmCount = 5.;
const float ArmCurvature = 1.4f;
const float StarDensity = 30.;
const float StarSize = .2;

vec3 mod289(vec3 x) {
  return x - floor(x * (1.0 / 289.0)) * 289.0;
}

vec4 mod289(vec4 x) {
  return x - floor(x * (1.0 / 289.0)) * 289.0;
}

vec4 permute(vec4 x) {
     return mod289(((x*34.0)+1.0)*x);
}

vec4 taylorInvSqrt(vec4 r)
{
  return 1.79284291400159 - 0.85373472095314 * r;
}

float permutes(float vx)
{
      return mod((34.0 * vx + 1.0) * vx, 289.0);
}

float snoise(vec3 v) {
      const vec2  C = vec2(1.0/6.0, 1.0/3.0) ;
      const vec4  D = vec4(0.0, 0.5, 1.0, 2.0);

      // First corner
      vec3 i  = floor(v + dot(v, C.yyy) );
      vec3 x0 =   v - i + dot(i, C.xxx) ;

      // Other corners
      vec3 g = step(x0.yzx, x0.xyz);
      vec3 l = 1.0 - g;
      vec3 i1 = min( g.xyz, l.zxy );
      vec3 i2 = max( g.xyz, l.zxy );

      //   x0 = x0 - 0.0 + 0.0 * C.xxx;
      //   x1 = x0 - i1  + 1.0 * C.xxx;
      //   x2 = x0 - i2  + 2.0 * C.xxx;
      //   x3 = x0 - 1.0 + 3.0 * C.xxx;
      vec3 x1 = x0 - i1 + C.xxx;
      vec3 x2 = x0 - i2 + C.yyy; // 2.0*C.x = 1/3 = C.y
      vec3 x3 = x0 - D.yyy;      // -1.0+3.0*C.x = -0.5 = -D.y

      // Permutations
      i = mod289(i);
      vec4 p = permute( permute( permute(
                i.z + vec4(0.0, i1.z, i2.z, 1.0 ))
              + i.y + vec4(0.0, i1.y, i2.y, 1.0 ))
              + i.x + vec4(0.0, i1.x, i2.x, 1.0 ));

      // Gradients: 7x7 points over a square, mapped onto an octahedron.
      // The ring size 17*17 = 289 is close to a multiple of 49 (49*6 = 294)
      float n_ = 0.142857142857; // 1.0/7.0
      vec3  ns = n_ * D.wyz - D.xzx;

      vec4 j = p - 49.0 * floor(p * ns.z * ns.z);  //  mod(p,7*7)

      vec4 x_ = floor(j * ns.z);
      vec4 y_ = floor(j - 7.0 * x_ );    // mod(j,N)

      vec4 x = x_ *ns.x + ns.yyyy;
      vec4 y = y_ *ns.x + ns.yyyy;
      vec4 h = 1.0 - abs(x) - abs(y);

      vec4 b0 = vec4( x.xy, y.xy );
      vec4 b1 = vec4( x.zw, y.zw );

      //vec4 s0 = vec4(lessThan(b0,0.0))*2.0 - 1.0;
      //vec4 s1 = vec4(lessThan(b1,0.0))*2.0 - 1.0;
      vec4 s0 = floor(b0)*2.0 + 1.0;
      vec4 s1 = floor(b1)*2.0 + 1.0;
      vec4 sh = -step(h, vec4(0.0));

      vec4 a0 = b0.xzyw + s0.xzyw*sh.xxyy ;
      vec4 a1 = b1.xzyw + s1.xzyw*sh.zzww ;

      vec3 p0 = vec3(a0.xy,h.x);
      vec3 p1 = vec3(a0.zw,h.y);
      vec3 p2 = vec3(a1.xy,h.z);
      vec3 p3 = vec3(a1.zw,h.w);

      //Normalise gradients
      vec4 norm = taylorInvSqrt(vec4(dot(p0,p0), dot(p1,p1), dot(p2, p2), dot(p3,p3)));
      p0 *= norm.x;
      p1 *= norm.y;
      p2 *= norm.z;
      p3 *= norm.w;

      // Mix final noise value
      vec4 m = max(0.6 - vec4(dot(x0,x0), dot(x1,x1), dot(x2,x2), dot(x3,x3)), 0.0);
      m = m * m;
      return 42.0 * dot( m*m, vec4( dot(p0,x0), dot(p1,x1), dot(p2,x2), dot(p3,x3) ) );
}

float fBm (vec3 P, float octaves, float lacunarity, float gain)
{
  float amp = 1.0;
  float sum = 0.0;
  float i;

  for (i = 0.0;  i < octaves;  i += 1.0) {
    sum += amp * snoise (P);
    amp *= gain;
    P *= lacunarity;
  }
  return sum;
}

float Turbulence (vec3 P, float octaves, float lacunarity, float gain)
{
  float amp = 1.0;
  float sum = 0.0;
  float i;

  for (i = 0.0;  i < octaves;  i += 1.0) {
    sum += abs(amp * snoise (P));
    amp *= gain;
    P *= lacunarity;
  }
  return sum;
}

uniform vec4 mainColor;

vec3 colorAtTime(float armMask, vec2 polar, float t) {
    float time0 = polar.y + t;
    vec3 noiseCoord = 0.4 * vec3(
        cos(polar.x) * time0,
        sin(polar.x) * time0,
        iTime * .05
    );

    vec3 color1 = vec3(47., 110., 168.)/255.;
    float layer1 = abs(fBm(noiseCoord * 3., 6., 2.1, 0.5));
    vec3 color2 = vec3(179., 97., 199.)/255.;
    float layer2 = pow(Turbulence(noiseCoord * 4., 6., 2.1, 0.45)*.6, 1.5);
    vec3 col = layer1 * color1 + layer2 * color2;
    return col * armMask;
}

// A single iteration of Bob Jenkins' One-At-A-Time hashing algorithm.
uint hash( uint x ) {
    x += ( x << 10u );
    x ^= ( x >>  6u );
    x += ( x <<  3u );
    x ^= ( x >> 11u );
    x += ( x << 15u );
    return x;
}



// Compound versions of the hashing algorithm I whipped together.
uint hash( uvec2 v ) { return hash( v.x ^ hash(v.y)                         ); }
uint hash( uvec3 v ) { return hash( v.x ^ hash(v.y) ^ hash(v.z)             ); }
uint hash( uvec4 v ) { return hash( v.x ^ hash(v.y) ^ hash(v.z) ^ hash(v.w) ); }



// Construct a float with half-open range [0:1] using low 23 bits.
// All zeroes yields 0.0, all ones yields the next smallest representable value below 1.0.
float floatConstruct( uint m ) {
    const uint ieeeMantissa = 0x007FFFFFu; // binary32 mantissa bitmask
    const uint ieeeOne      = 0x3F800000u; // 1.0 in IEEE binary32

    m &= ieeeMantissa;                     // Keep only mantissa bits (fractional part)
    m |= ieeeOne;                          // Add fractional part to 1.0

    float  f = uintBitsToFloat( m );       // Range [1:2]
    return f - 1.0;                        // Range [0:1]
}



// Pseudo-random value in half-open range [0:1].
float random( float x ) { return floatConstruct(hash(floatBitsToUint(x))); }
float random( vec2  v ) { return floatConstruct(hash(floatBitsToUint(v))); }
float random( vec3  v ) { return floatConstruct(hash(floatBitsToUint(v))); }
float random( vec4  v ) { return floatConstruct(hash(floatBitsToUint(v))); }

vec2 GetCellCenter(vec2 uv) {
    return vec2(
        mod(mod(uv, 1.) * StarDensity, 1.) * 2. - 1.
    );
}

vec2 RandomCellOffset(vec2 uv) {
    vec3 grid = vec3(floor(uv*StarDensity)/StarDensity, 0.);
    vec3 gridNoise = vec3(
        random(grid+vec3(20., 50., -10.)),
        random(grid+vec3(10., -30., -26.)),
        random(grid+vec3(-78., -29., 50.))
    );

    return normalize(gridNoise.xy * 2. - 1.) *
        clamp(0., 1., gridNoise.z);
}

vec2 RandomActualCellOffset(vec2 uv) {
    float d = 1./StarDensity;
    vec2 center = GetCellCenter(uv);
    vec2 c0 = RandomCellOffset(uv);
    for (int i = 0; i < 9; i++) {
        int dx = i % 3 - 1; // -1, 0 or 1
        int dy = i / 3 - 1; // -1, 0 or 1
        vec2 offset = vec2(dx, dy);
        vec2 c1 = RandomCellOffset(uv + d*offset) + 2.*offset;
        if (distance(center, c0) > distance(center, c1)) {
            c0 = c1;
        }
    }
    return c0;
}

void mainImage( out vec4 fragColor, in vec2 fragCoord ) {
    // -1 to 1
    //vec2 uv = 2.0 * ((fragCoord + vec2(0.0, iResolution.x*.20))/iResolution.x - .5);

    vec2 uv = fragCoord/iResolution.x + vec2(0., .15);


    vec2 backUV = uv;
    {
        // We center texture coordinates around 0.0
        // Perspective plane equation
        float pcenter=0.1;
        float pmult=1.6;
        float vc=(uv.y-pcenter)*pmult;
        float uu=((uv.x-0.5)*1.0)/(1.0-vc);
        float vv=((vc-0.6)/(1.0-vc));
        uv=vec2(uu,vv);
    }

    vec2 middle = (uv + .5) * 2. - 1.;
    //backUV += normalize(middle) *
    //    (.2*pow(length(middle), .3)) * .5;




    vec2 blackHoleUV = uv * .125;

    float yShift = max(0.0, 1.8 - pow(length(blackHoleUV),0.2));
    blackHoleUV.y += yShift*.4;

    vec2 blackholePolar = vec2(
        atan(blackHoleUV.y,blackHoleUV.x),
        length(blackHoleUV)
    );

    //uv.y += yShift;

    vec2 polar = vec2(
        atan(uv.y, uv.x), // The angle (radians)
        length(uv)        // The distance (0 - sqrt(2))
    );

    // Bend the arms
    float curveFactor = pow(polar.y, 0.7) * ArmCurvature;
    curveFactor -= pow(max(0.0, 1.3 - length(uv)*2.),10.);
    polar.x += curveFactor;

    blackholePolar.x += pow(blackholePolar.y, 0.7) * ArmCurvature;
    blackholePolar.x += curveFactor;

    // Masks out the arms (0 - 1)
    float armU = mod(polar.x*ArmCount, 6.28)/6.28; // 0 to 1
    armU = armU * 2. - 1.;                   // -1 to 1
    armU = abs(armU);                        // 0 to 1
    armU = smoothstep(0.0, 1.0, armU);       // 0 to 1

    float barmU = mod(blackholePolar.x*ArmCount, 6.28)/6.28; // 0 to 1
    barmU = barmU * 2. - 1.;                   // -1 to 1
    barmU = abs(barmU);                        // 0 to 1
    barmU = smoothstep(0.0, 1.0, barmU);       // 0 to 1

    // Distance along the arm (0 - inf)
    float armV = curveFactor * armU;

    // Which arm? (0.0, 0.25, 0.5, 0.75, 1.0)
    float armW = floor(
        mod(polar.x+3.14, 6.28)*ArmCount / 6.28
    ) / (ArmCount - 1.0);

    // Noise Coordinate
    float time0 = polar.y/10.0 + mod(iTime, 5.0) / 5.0;
    float time1 = polar.y/10.0 + mod(iTime + 2.5, 5.0) / 5.0;
    float t     = abs(mod(iTime, 5.0) * 0.4 - 1.0);
    vec3 col0 = colorAtTime(armU, polar, time0);
    vec3 col1 = colorAtTime(armU, polar, time1);
    vec3 col = mix(col0, col1, t);

    // Noise Coordinate
    float btime0 = blackholePolar.y/10.0 + mod(iTime, 5.0) / 5.0;
    float btime1 = blackholePolar.y/10.0 + mod(iTime + 2.5, 5.0) / 5.0;
    vec3 bcol0 = colorAtTime(barmU, blackholePolar, btime0);
    vec3 bcol1 = colorAtTime(barmU, blackholePolar, btime1);
    vec3 bcol = mix(bcol0, bcol1, t);

    float hole = clamp(pow(length(uv)*5.,3.), 0.0, 1.0);
    float light = clamp(pow(1.0-length(uv),3.), 0.0, 1.0);

    float galaxyMask = clamp(1.0 - length(uv)*0.3, 0.0, 1.0);
    vec3 galaxy = galaxyMask * clamp(
        vec3(pow(hole,1.3)*(pow(armU,1.5)*light+col)),
        vec3(0.0), vec3(1.0));



    float backNoise0 = pow(abs(fBm(vec3(backUV*3., -1.),
        6., 2.4, 0.4)), 1.5);
    float backNoise1 = clamp(abs(Turbulence(vec3(backUV, 5.),
        7., 2.2, 0.4)-0.3), 0., 1.);
    float backNoise2 = clamp(1.0-pow(abs(Turbulence(vec3(backUV*3., -1.),
        6., 2.4, 0.4)), 0.5), 0.0, 1.0);

    vec3 backCol =
        //vec3(39., 78., 145.)/255. * backNoise0 * .7 +
        vec3(121., 34., 168.)/255. * backNoise1 * 0.4 +
        vec3(49., 187., 212.)/255. * backNoise2 * .7
    ;

    fragColor = vec4(backCol + galaxy, 1.0);



    vec2 gridNoise = RandomActualCellOffset(backUV);

    vec2 cellCenter = GetCellCenter(backUV);

    float starSizeT = abs(2.*mod(iTime*.1+random(RandomCellOffset(backUV)), 1.0)-1.);
    float starSize = starSizeT * StarSize + 0.001;
    float stars = (1.0/starSize)*(clamp(1.-distance(gridNoise, cellCenter),1.-starSize,1.)-(1.-starSize));

    stars = pow(stars, 5.);

    fragColor += vec4(vec3(stars), 1.);

    //fragColor = vec4(vec3(vec2(stars), 0.), 1.);
    //fragColor = vec4(vec3(gridNoise.xy, 0.), 1.);
}