layers

layers.glsl

#version 150

uniform float time;
uniform vec2 resolution;
uniform vec2 mouse;
uniform vec3 spectrum;

uniform sampler2D texture0;
uniform sampler2D texture1;
uniform sampler2D texture2;
uniform sampler2D texture3;
uniform sampler2D prevFrame;
uniform sampler2D prevPass;

in VertexData
{
    vec4 v_position;
    vec3 v_normal;
    vec2 v_texcoord;
} inData;

out vec4 fragColor;

// -----------------------------==>> COMMON <<==--------------------------------

#define PI              3.141592654
#define TAU             (2.0*PI)
#define TIME            time
#define RESOLUTION      resolution
#define SCA(a)          vec2(sin(a), cos(a))

#define LESS(a,b,c)     mix(a,b,step(0.,c))
#define SABS(x,k)       LESS((.5/(k))*(x)*(x)+(k)*.5,abs(x),abs(x)-(k))

#define MROT(a)         mat2(cos(a), sin(a), -sin(a), cos(a))

#define ORT(p)              vec2((p).y, -(p).x)

#define BPM             240.0

const vec2 hexcell_sz       = vec2(1.0, sqrt(3.0));
const vec2 hexcell_hsz      = 0.5*hexcell_sz;

const mat2 rot120           = MROT(TAU/3.0);
const mat2 rotSome          = MROT(1.0);

const vec2 hexcell_n1       = normalize(hexcell_sz.yx);
const vec2 hexcell_n2       = hexcell_n1*vec2(-1.0, 1.0);

const float hexcell_radius  = 0.5/sqrt(3.0);
const vec2 hexcell_off1     = 2.0*hexcell_radius*hexcell_n1;
const vec2 hexcell_off2     = rot120*hexcell_off1;
const vec2 hexcell_off3     = rot120*hexcell_off2;

const vec3 std_gamma        = vec3(2.2, 2.2, 2.2);

const float is3             = 1.0/sqrt(3.0);
const vec2 cdir             = normalize(vec2(1.0, is3));
const vec2 flipy            = vec2(1.0, -1.0);

const vec2 coords[6] = vec2[6](
  is3*cdir*1.0/3.0,
  is3*cdir*2.0/3.0,
  vec2(0.5, is3/6.0),
  vec2(0.5, -is3/6.0),
  is3*cdir*2.0/3.0*flipy,
  is3*cdir*1.0/3.0*flipy
  );

const vec2 dcoords[6] = vec2[6](
  ORT(cdir),
  ORT(cdir),
  vec2(-1.0, 0.0),
  vec2(-1.0, 0.0),
  ORT(-cdir*flipy),
  ORT(-cdir*flipy)
  );

const int corners[] = int[](
  0, 1, 2, 3, 4, 5, 
  0, 1, 2, 4, 3, 5, 
  0, 1, 2, 5, 3, 4, 
  0, 2, 1, 3, 4, 5, 
  0, 2, 1, 4, 3, 5, 
  0, 2, 1, 5, 3, 4, 
  0, 3, 1, 2, 4, 5, 
  0, 3, 1, 4, 2, 5, 
  0, 3, 1, 5, 2, 4, 
  0, 4, 1, 2, 3, 5, 
  0, 4, 1, 3, 2, 5, 
  0, 4, 1, 5, 2, 3, 
  0, 5, 1, 2, 3, 4, 
  0, 5, 1, 3, 2, 4, 
  0, 5, 1, 4, 2, 3
  );
const int noCorners = corners.length()/6;

vec3 hsv2rgb(vec3 c) {
  const vec4 K = vec4(1.0, 2.0 / 3.0, 1.0 / 3.0, 3.0);
  vec3 p = abs(fract(c.xxx + K.xyz) * 6.0 - K.www);
  return c.z * mix(K.xxx, clamp(p - K.xxx, 0.0, 1.0), c.y);
}

vec3 rgb2hsv(vec3 c) {
  const vec4 K = vec4(0.0, -1.0 / 3.0, 2.0 / 3.0, -1.0);
  vec4 p = mix(vec4(c.bg, K.wz), vec4(c.gb, K.xy), step(c.b, c.g));
  vec4 q = mix(vec4(p.xyw, c.r), vec4(c.r, p.yzx), step(p.x, c.r));

  float d = q.x - min(q.w, q.y);
  float e = 1.0e-10;
  return vec3(abs(q.z + (q.w - q.y) / (6.0 * d + e)), d / (q.x + e), q.x);
}

float hash(in float co) {
  return fract(sin(co*12.9898) * 13758.5453);
}

float hash(in vec2 co) {
  return fract(sin(dot(co, vec2(12.9898,58.233))) * 13758.5453);
}

float hash(in vec3 co) {
  return fract(sin(dot(co, vec3(12.9898,58.233, 12.9898+58.233))) * 13758.5453);
}

float psin(float a) {
  return 0.5 + 0.5*sin(a);
}

float pcos(float a) {
  return 0.5 + 0.5*cos(a);
}

vec2 toPolar(vec2 p) {
  return vec2(length(p), atan(p.y, p.x));
}

vec2 toRect(vec2 p) {
  return vec2(p.x*cos(p.y), p.x*sin(p.y));
}

float mod1(inout float p, float size) {
  float halfsize = size*0.5;
  float c = floor((p + halfsize)/size);
  p = mod(p + halfsize, size) - halfsize;
  return c;
}

vec2 mod2_1(inout vec2 p) {
  vec2 c = floor(p + 0.5);
  p = fract(p + 0.5) - 0.5;
  return c;
}


vec2 mod2(inout vec2 p, vec2 size) {
  vec2 c = floor((p + size*0.5)/size);
  p = mod(p + size*0.5,size) - size*0.5;
  return c;
}

vec2 modMirror2(inout vec2 p, vec2 size) {
  vec2 halfsize = size*0.5;
  vec2 c = floor((p + halfsize)/size);
  p = mod(p + halfsize, size) - halfsize;
  p *= mod(c,vec2(2))*2 - vec2(1);
  return c;
}

float modPolar(inout vec2 p, float repetitions) {
  float angle = 2.0*PI/repetitions;
  float a = atan(p.y, p.x) + angle/2.;
  float r = length(p);
  float c = floor(a/angle);
  a = mod(a,angle) - angle/2.;
  p = vec2(cos(a), sin(a))*r;
  // For an odd number of repetitions, fix cell index of the cell in -x direction
  // (cell index would be e.g. -5 and 5 in the two halves of the cell):
  if (abs(c) >= (repetitions/2.0)) c = abs(c);
  return c;
}

float hex(vec2 p, float r) {
  const vec3 k = vec3(-sqrt(3.0)/2.0,1.0/2.0,sqrt(3.0)/3.0);
  p = p.yx;
  p = abs(p);
  p -= 2.0*min(dot(k.xy,p),0.0)*k.xy;
  p -= vec2(clamp(p.x, -k.z*r, k.z*r), r);
  return length(p)*sign(p.y);
}

vec2 hextile(inout vec2 p) {
  // See Art of Code: Hexagonal Tiling Explained!
  // https://www.youtube.com/watch?v=VmrIDyYiJBA

  vec2 p1 = mod(p, hexcell_sz)-hexcell_hsz;
  vec2 p2 = mod(p - hexcell_hsz, hexcell_sz)-hexcell_hsz;
  vec2 p3 = mix(p2, p1, vec2(dot(p1, p1) < dot(p2, p2)));
  vec2 n = round((p3 - p + hexcell_hsz)/hexcell_hsz);
  
  p = p3;

  // Rounding to make hextile 0,0 well behaved
  return round(n*2.0)/2.0;
}

vec3 triTile(inout vec2 p, out vec3 bc, out vec3 de) {
  const mat2 m1 = mat2(1, -1.0/sqrt(3.0), 0,2./sqrt(3.0))*1;
  const mat2 m2 = mat2(1, 0.5, 0,sqrt(3.0)/2.0);

  vec2 u = p;
  vec2 v = p;
  
  u *= m1;
  
  vec3 g = vec3(u,1.-u.x-u.y);
  vec3 id = floor(g);
  
  g = fract(g);
  if (length(g)>1.) {
    g = 1.-g;
  }
  bc = g;
  de = 1.0 - abs(2.*fract(g)-1.);
  u = (g.xy-ceil(1.-g.z)/3.) * m2;

  p = u;
  
  return id;
}

float smoother(float d, float s) {
  return s*tanh(d/s);
}


float pmin(float a, float b, float k) {
  float h = clamp( 0.5+0.5*(b-a)/k, 0.0, 1.0 );
  return mix( b, a, h ) - k*h*(1.0-h);
}

float dot2(vec2 v) { return dot(v,v); }

// IQ Bezier: https://www.shadertoy.com/view/MlKcDD
float bezier(vec2 pos, vec2 A, vec2 B, vec2 C) {    
  const float sqrt3 = sqrt(3.0);
  vec2 a = B - A;
  vec2 b = A - 2.0*B + C;
  vec2 c = a * 2.0;
  vec2 d = A - pos;

  float kk = 1.0/dot(b,b);
  float kx = kk * dot(a,b);
  float ky = kk * (2.0*dot(a,a)+dot(d,b))/3.0;
  float kz = kk * dot(d,a);      

  float res = 0.0;

  float p = ky - kx*kx;
  float p3 = p*p*p;
  float q = kx*(2.0*kx*kx - 3.0*ky) + kz;
  float h = q*q + 4.0*p3;

  if(h>=0.0) {   // 1 root
      h = sqrt(h);
      vec2 x = (vec2(h,-h)-q)/2.0;
      vec2 uv = sign(x)*pow(abs(x), vec2(1.0/3.0));
      float t = clamp(uv.x+uv.y-kx, 0.0, 1.0);
      res = dot2(d+(c+b*t)*t);
  } else {   // 3 roots
      float z = sqrt(-p);
      float v = acos(q/(p*z*2.0))/3.0;
      float m = cos(v);
      float n = sin(v)*sqrt3;
      vec3  t = clamp(vec3(m+m,-n-m,n-m)*z-kx, 0.0, 1.0);
      res = min(dot2(d+(c+b*t.x)*t.x), dot2(d+(c+b*t.y)*t.y));
      // the third root cannot be the closest. See https://www.shadertoy.com/view/4dsfRS
      // res = min(res,dot2(d+(c+b*t.z)*t.z));
  }
  
  return sqrt(res);
}

float bezier2(vec2 p, float f, vec2 p0, vec2 dp0, vec2 p1, vec2 dp1) {
  float dist = length(p0 - p1);
  float hdist = 0.5*f*dist;
  vec2 mp0 = p0 + hdist*dp0;
  vec2 mp1 = p1 + hdist*dp1;
  vec2 jp = (mp0 + mp1)*0.5;
  float d0 = bezier(p, p0, mp0, jp);
  float d1 = bezier(p, p1, mp1, jp);
  
  float d = d0;
  d = min(d, d1);
  return d;
}

float onoise(vec2 x) {
  float a = sin(x.x);
  float b = sin(x.y);
  float c = mix(a, b, psin(TAU*tanh(a*b+a+b)));

  return c;
}

float snoise(vec2 x) {
  float a = sin(x.x);
  float b = sin(x.y);

  return a*b;
}

float vnoise(vec2 x) {
  vec2 i = floor(x);
  vec2 w = fract(x);
    
#if 1
  // quintic interpolation
  vec2 u = w*w*w*(w*(w*6.0-15.0)+10.0);
#else
  // cubic interpolation
  vec2 u = w*w*(3.0-2.0*w);
#endif    

  float a = hash(i+vec2(0.0,0.0));
  float b = hash(i+vec2(1.0,0.0));
  float c = hash(i+vec2(0.0,1.0));
  float d = hash(i+vec2(1.0,1.0));
    
  float k0 =   a;
  float k1 =   b - a;
  float k2 =   c - a;
  float k3 =   d - c + a - b;

  float aa = mix(a, b, u.x);
  float bb = mix(c, d, u.x);
  float cc = mix(aa, bb, u.y);
  
  return k0 + k1*u.x + k2*u.y + k3*u.x*u.y;
}


float circle(vec2 p, float r) {
  return length(p) - r;
}

float box(vec2 p, vec2 b) {
  vec2 d = abs(p)-b;
  return length(max(d,0.0)) + min(max(d.x,d.y),0.0);
}

float starn(vec2 p, float r, int n, float m) {
  // next 4 lines can be precomputed for a given shape
  float an = PI/float(n);
  float en = PI/m;  // m is between 2 and n
  vec2  acs = vec2(cos(an),sin(an));
  vec2  ecs = vec2(cos(en),sin(en)); // ecs=vec2(0,1) for regular polygon,

  // Uses "wrongly" SABS to smooth inner corners
  float bn = mod(atan(p.x,p.y),2.0*an) - an;
  p = length(p)*vec2(cos(bn),SABS(sin(bn), 0.15));
//  p = length(p)*vec2(cos(bn),abs(sin(bn)));
  p -= r*acs;
  p += ecs*clamp( -dot(p,ecs), 0.0, r*acs.y/ecs.y);
  return length(p)*sign(p.x);
}

vec4 alphaBlend(vec4 back, vec4 front) {
  vec3 colb = back.xyz;
  vec3 colf = front.xyz;
  vec3 xyz = mix(colb*back.w, colf.xyz, front.w);
  float w = mix(back.w, 1.0, front.w);
  return vec4(xyz, w);
}

vec3 alphaBlend(vec3 back, vec4 front) {
  vec3 colb = back.xyz;
  vec3 colf = front.xyz;
  vec3 xyz = mix(colb, colf.xyz, front.w);
  return xyz;
}


vec4 alphaBlendGamma(vec4 back, vec4 front, vec3 gamma) {
  vec3 colb = max(back.xyz, 0.0);
  vec3 colf = max(front.xyz, 0.0);
  colb = pow(colb, gamma);
  colf = pow(colf, gamma);
  vec3 xyz = mix(colb*back.w, colf.xyz, front.w);
  float w = mix(back.w, 1.0, front.w);
  return max(vec4(pow(xyz, 1.0/gamma), w), 0.0);
}

vec3 alphaBlendGamma(vec3 back, vec4 front, vec3 gamma) {
  vec3 colb = max(back.xyz, 0.0);
  vec3 colf = max(front.xyz, 0.0);;
  
  colb = pow(colb, gamma);
  colf = pow(colf, gamma);
  vec3 xyz = mix(colb, colf.xyz, front.w);
  return pow(xyz, 1.0/gamma);
}


void rot(inout vec2 p, float a) {
  float c = cos(a);
  float s = sin(a);
  p = vec2(c*p.x + s*p.y, -s*p.x + c*p.y);
}

vec3 offset_0(float z) {
  float a = z;
  vec2 p = vec2(0.0);
  return vec3(p, z);
}

vec3 offset_1(float z) {
  float a = z;
  vec2 p = -0.075*(vec2(cos(a), sin(a*sqrt(2.0))) + vec2(cos(a*sqrt(0.75)), sin(a*sqrt(0.5))));
  return vec3(p, z);
}

vec3 offset(float z) {
  return offset_1(z);
}

vec3 doffset(float z) {
  float eps = 0.1;
  return 0.5*(offset(z + eps) - offset(z - eps))/eps;
}

vec3 ddoffset(float z) {
  float eps = 0.1;
  return 0.125*(doffset(z + eps) - doffset(z - eps))/eps;
}

const float truchetCell_linew = 0.05;

vec2 truchetCell4_0(vec2 p) {
  return vec2(1E6, 1E6);  
}

vec2 truchetCell4_1(vec2 p) {
  float d0 = circle(p + vec2(0.5), 0.5);
  float d1 = circle(p - vec2(0.5), 0.5);
  float d2 = circle(p, 0.25);
  
  float d = d0;
  d = min(d, d1);
  
  return vec2(d, 1E6);  
}

vec2 truchetCell4_2(vec2 p) {
  float d0 = abs(p.x);
  float d1 = abs(p.y);
  float d2 = circle(p, 0.25);
  
  float d = d0;
  d = min(d, d1);
  d = min(d, d2);

  return vec2(d, d2+truchetCell_linew);  
}

vec2 truchet4(vec2 p, float n) {
  vec2 pt = p;

  vec2 nt = mod2_1(pt);
  float r = hash(nt+10.0*n);
  float r1 = fract(13.0*r);
  float r2 = fract(23.0*r);
  pt *= vec2(-1) + 2*vec2(step(0.5, r1), step(0.5, r2));
  vec2 dt = vec2(1E6);
  
  
  if (r < 0.5) {
      dt = truchetCell4_1(pt);
  } else {
      dt = truchetCell4_2(pt);
  }

  vec2 apt = abs(pt);
  vec2 d0 = vec2(length(apt - vec2(0.5, 0.0)), 1E6);
  vec2 d1 = vec2(length(apt - vec2(0.0, 0.5)), 1E6);

  vec2 d = d0;
  d = min(d, d1);
  d = min(d, dt);

  d.x = abs(d.x) - truchetCell_linew;

  return d;
}

vec2 tnoise4(vec2 p, float n) {
  vec2 d = truchet4(p, n);
  float s = 0.1;
  vec2 nn = smoothstep(-s, s, -d);
  return nn;
}

vec2 truchetCell6_0(vec2 p) {
  float d0 = abs(p.y);
  float d1 = abs(dot(hexcell_n1, p));
  float d2 = abs(dot(hexcell_n2, p));
  float d3 = circle(p, hexcell_radius);
 
  float d = 1E6;
  d = min(d, d0);
  d = min(d, d1);
  d = min(d, d2);
  d = min(d, d3);
   
  return vec2(d, -abs(d)+truchetCell_linew*2);
}

vec2 truchetCell6_1(vec2 p) {
  p = abs(p);
  float d0 = circle(p-hexcell_off1, hexcell_radius); 
  float d = d0;
  return vec2(d, -abs(d)+truchetCell_linew*2);
}

vec2 truchetCell6_2(vec2 p) {
  float d0 = circle(p-hexcell_off1, hexcell_radius);
  float d1 = circle(p-hexcell_off2, hexcell_radius);
  float d2 = circle(p-hexcell_off3, hexcell_radius);
   
  float d = 1E6;
  d = min(d, d0);
  d = min(d, d1);
  d = min(d, d2);
 
  return vec2(d, -abs(d)+truchetCell_linew*2);
}

vec2 truchet6(vec2 p, float n) {
   vec2 nn = hextile(p);
   float r1 = hash(nn+n*10);
   float r2 = fract(23.0*r1);
   rot(p, floor(r2*6.0)* TAU/6.0);
   vec2 d = vec2(1E6);
   if (r1 < 0.5) {
     d = truchetCell6_2(p);
   } else if (r1 < 0.8) {
     d = truchetCell6_1(p);
   } else {
     d = truchetCell6_0(p);
   }
    
   d.x = abs(d.x) - truchetCell_linew;

   return d;
}

vec2 tnoise6(vec2 p, float n) {
  vec2 d = truchet6(p, n);
  float s = 0.1;
  vec2 nn = smoothstep(-s, s, -d);
  return nn;
}


vec2 dtruchet6(vec2 p, float n) {
  vec2 hp = p;
  vec2 hn = hextile(hp);
  
  vec2 pp = toPolar(hp);
  float pn = mod1(pp.y, TAU/6.0);
  vec2 tp = toRect(pp);

  pn = mod(pn+3.0, 6.0);
  vec3 nn = vec3(hn, pn);

  float r = hash(nn+n);
  int sel = int(float(noCorners)*r);
  int off = sel*6;
  

  const float sw = 0.05;

  float d = 1E6;

//  float linew = truchetCell_linew*0.125*mix(0.25, 1.0, 1.0-psin(1*TIME+(0.5*p.x+p.y)*10));
  const float linew = truchetCell_linew*0.125;


  for (int i = 0; i < 3; ++i) {
    int c0 = corners[off + i*2 + 0];
    int c1 = corners[off + i*2 + 1];
    
    int c = max(c0, c1) - min(c0, c1);
    
    vec2 p0 = coords[c0];
    vec2 p1 = coords[c1];

    vec2 dp0 = dcoords[c0];
    vec2 dp1 = dcoords[c1];
    
    float mi = 0.5;    
    float mx = 0.5;    
    
    float rr = fract(r*27.0*float(i+1));
    switch(c) {
      case 1:
        mx = 1.75;
        break;
      case 2:
        mx = .95;
        break;
      case 3:
        mx = 1.5;
        break;
      case 4:
        mx = 0.75;
        break;
      case 5:
        mx = 1.95;
        break;
      default:
        break;
    }
    
    float f = mix(mi, mx, rr);
    
    float dd = (bezier2(tp, f, p0, dp0, p1, dp1)-linew);
    d = pmin(d, dd, linew*10);
    
  }
  
  float dd = abs(d) - linew;
  return vec2(d, -dd);
}

// -----------------------------==>> COMMON <<==--------------------------------

vec4 plane_0(vec3 ro, vec3 rd, vec3 pp, float aa, float n) {
  vec2 p = pp.xy;
  float z = pp.z;

  p -= (1+5*(pp.z - ro.z))*offset(z).xy;
  rot(p, n*sqrt(0.5));
  vec3 nor = normalize(vec3(0.0, 0.0, -1.0) + 0.10*vec3(psin(p.x*TAU*20), psin(p.y*TAU*sqrt(13)+TIME), 0.0));
  float s = 0.5;
  p *= s;
  vec2 nn = mod2_1(p);
  float r = hash(vec3(nn, n));
  rot(p, ((TAU*r+n)*0.25));
  float d = starn(p, 0.20, 3 + 2*int(3.0*r), 3);
  d -= 0.06;
  d/=s;
  
  
  vec3 lp1 = offset(ro.z-3)+vec3(-1.5, 2.0, 1.0);
  vec3 ld1 = normalize(pp - lp1);
  float diff1 = mix(0.5, 1.0, max(dot(-ld1, nor), 0.0));
  float hspe1 = pow(max(dot(-ld1, reflect(rd, nor)), 0.0), 500.);
  float lspe1 = pow(max(dot(-ld1, reflect(rd, nor)), 0.0), 50.);

  vec3 lp2 = offset(ro.z-5)+vec3(1.0, -1.5, 3.0);
  vec3 ld2 = normalize(pp - lp2);
  float diff2 = mix(0.5, 1.0, max(dot(-ld2, nor), 0.0));
  float hspe2 = pow(max(dot(-ld2, reflect(rd, nor)), 0.0), 500.);
  float lspe2 = pow(max(dot(-ld2, reflect(rd, nor)), 0.0), 50.);

  float db = abs(d) - (0.06);
  db = abs(db) - 0.03;
  db = abs(db) - 0.00;
  db = max(db, -d+0.03);
  const vec3 colb = 1.5*vec3(1.0, 0.7, 0.5);
  float tb = exp(-(db)*30*(1.0 - 10*aa));
  vec4 cb = vec4(colb, tb);

  float ds = -d+0.03;
  vec3 cols = hsv2rgb(vec3(337.0/360.0+0.1*sin(n*0.3), 0.8, 0.54+0.2*sin(n*0.3)));
  float ts = 1.0 - smoothstep(-aa, 0.0, ds);
  vec4 cs =  vec4(cols, ts*0.94);

  cs.xyz *= tanh(diff1*diff1 + diff2*diff2);

  float td = length(pp-ro);
  float hs = 1.0 - smoothstep(1.0*3.0, 3.0*3.0, td);
  vec4 ct = alphaBlendGamma(cs, cb, vec3(3.2, 1.5, 3.2));

  ct.xyz  += colb*hspe1;
  ct.xyz  += 0.25*colb.xzy*lspe1;

  ct.xyz  += colb*hspe2;
  ct.xyz  += 0.25*colb.xzy*lspe2;

  return ct;
}

vec2 df_0(vec2 p, float offz, float n, float nz) {
  return truchet4(p, n);
}

vec2 df_1(vec2 p, float offz, float n, float nz) {
  const float s = 0.25;
  float r = hash(n*10);
  rot(p, r*TAU);
  p -= nz;
  return truchet6(p*s, n)/s;
}

vec2 df_2(vec2 p, float offz, float n, float nz) {

  float tn = TAU*(n*0.05);
//  tn = 0;
//  rot(p, tn+TIME/30.0);
//  float d0 = box(p, vec2(p));
//  float d0 = circle(p, 1.0);
//  float d0 = box(p, vec2(5*psin(tn+ TAU*TIME/15.0)));
//  float d0 = starn(p, 4.5, 10,1+5*psin(tn));
  float d0 = starn(p, 4.5, 15,1.1);
//  float d0 = starn(p, 4.5, 15,1.1+20*psin(tn+ 0.*TAU*TIME/5.0));
//  float d0 = starn(p, 4.5, 15,1.1+5*psin(tn+ 0.*TAU*TIME/5.0));
  float per = 0.31;
  d0 = mod(d0, per) - per*0.5;
  d0 = abs(d0) - 0.125;

//  return vec2(d0, d0+0.025);
  return vec2(d0, -d0);
}

vec2 df_3(vec2 p, float offz, float n, float nz) {
  const float s = 0.25;
  float r = hash(n*10);
  rot(p, r*TAU);
  p -= nz;
  return dtruchet6(p*s, n)/s;
}

vec2 df_4(vec2 p, float offz, float n, float nz) {
  const float s = 0.25;
  p *= s;
  vec2 hp0 = p;
  vec2 hp1 = p+vec2(0., 1.0/sqrt(3.0));
  vec2 hn0 = hextile(hp0);
  vec2 hn1 = hextile(hp1);

  float sz = mix(0.25, 0.5, 1.0-pcos(TAU*nz+TAU*length(p)));
//  sz = 0.5;
  float d0 = hex(hp0, sz);
  d0 = abs(d0) - 0.025;
  
  float d1 = hex(hp1, sz);
  d1 = abs(d1) - 0.025;

  float d = d0;
  d = pmin(d, d1, 0.25);
  
  d /= s;

  const float linew = truchetCell_linew*0.5;
  float dd = abs(d) - linew;
  return vec2(d, -dd);
}


vec2 df(vec2 p, float offz, float n, float nz) {
  return df_4(p, offz, n, nz);
}

vec4 plane_1(vec3 ro, vec3 rd, vec3 pp, float aa, float n) {
  vec2 p = pp.xy;
  float z = pp.z;
  float nz = pp.z-ro.z;

  p -= (1+0*(pp.z - ro.z))*offset(z).xy;
//  rot(p, TAU*n/20.0+1.0*TAU*TIME/60.0);
  const float st = 5;
  vec2 dt = df(p*st, nz, n, nz);
  dt = -dt;
  dt /= st;

  vec2 pc = p;  
  float dc = -circle(pc, 0.5);
//  float dc = -starn(pc, 0.5, 15, mix(2.0, 10.0, tanh(0.5*nz)));

  vec3 pcol = vec3(1.0);
  vec3 bct = mix(vec3(0.0), pcol, smoothstep(-aa, aa, -dt.x));
  float tt = smoothstep(-aa, aa, -dt.y);

  float ss = aa;
  float ci = smoothstep(-ss, ss, sin(1.25*n+30.0/60.0*TAU*TIME+dc));
  vec3 bcc = vec3(1.0)*smoothstep(-aa, aa, -dc-0.02);
  bcc *= bct;
//  bcc = mix(bcc, 1.0-bcc, ci);
  float tc = max(smoothstep(-aa, aa, -dc)*tt, 1.0*smoothstep(-aa, aa, -abs(dc)+0.02));
  vec4 cc = vec4(bcc, tc);
//  cc = vec4(bct, tt)  ;
  return cc;
}

float globalCloudDensity(vec2 p, float n) {
  p*=0.1;
  float gcd = vnoise(p+10*hash(n)+100);
  return gcd;
}

float localCloudDensity(vec2 p, float n) {
  p*=1.0;
  const float aa = -0.45;
  const mat2 pp = 2.03*rotSome;
  float a = 0.5;
  float s = 0.0;
  p += 10*hash(n)+100;

  s += a*vnoise(p); a *= aa; p *= pp;
  s += a*vnoise(p); a *= aa; p *= pp;
  s += a*vnoise(p); a *= aa; p *= pp;
  s += a*vnoise(p); a *= aa; p *= pp;
  s += a*vnoise(p); a *= aa; p *= pp;
  return s*2.1-0.0;
}

vec4 plane_2(vec3 ro, vec3 rd, vec3 pp, float aa, float n) {
  vec2 p = pp.xy;
  float z = pp.z;
  float nz = pp.z-ro.z;

  float ds = 1E6;
  float r = 0.0;
  const float s = 1.5;
  p+=1.35;
  vec2 ps = p;
  ps *= s;
  for (int i = 0; i < 5; ++i) {
    ps *= rotSome;
    vec2 ips = ps;
    vec2 ipn = mod2_1(ips);
    float ir = hash(ipn+n*100+i*10);
    ips = ips -0.3*vec2(ir, fract(-23.0*ir));
    float ids = length(ips)-0.0025;
    if (ids < ds) {
      r = ir;
      ds = min(ds, ids);
    }
  }

  float gcd = globalCloudDensity(p, n);
  /*
  const float poff = -5.0;
  vec2 ppp = toPolar(p-poff);
  ppp.y += length(p-poff)*0.2;
  p = toRect(ppp)+poff;
*/
  float hues = mix(0.6, 0.8, r*r);
  float sats = mix(0.5, 0.0, sqrt(r));
  float bris = mix(0.5, 1.0, r);
  float ts = pow(max(1.0-ds, 0.0), mix(200.0, 100.0, gcd)/sqrt(s));
  vec3 cols = 3.0*hsv2rgb(vec3(hues, sats, bris));
  vec4 cs = vec4(cols, ts*sqrt(gcd));


  float cd = gcd*localCloudDensity(p, n);
  float cdo = gcd*localCloudDensity(p+vec2(0.125, 0.25), n);
  const float level0 = 0.0;
  const float level1 = 0.05;
  float cli = mix(-0.1, 1.0, 0.5 + 0.5*tanh(10*(cd-cdo)));

  float huec = (mix(-0.2, 0.05, (cd))+0.05)-0.15*pcos(2*pp.z);
  float tc = clamp(cd, 0.0, 1.0);
  float satc = 0.5;
  float bric = 1.0;
  vec3 colc = hsv2rgb(vec3(huec, satc, bric))+cli*0.65;
  vec4 cc = vec4(colc, tc);
  
  vec4 ct = alphaBlendGamma(cs, cc, std_gamma);
  return ct;
}

float twirlyHeight(vec2 p, float n, out vec2 diff) {
  const float aan = 0.45;
  const float ppn = 2.0+0.2;
  
  float s = 0.0;
  float d = 0.0;
  float an = 1.0;
  vec2 pn = 4*p+n*10;
  vec2 opn = pn;

  const int md = 1;
  const int mx = 4;
  
  for (int i = 0; i < md; ++i) {
    s += an*(vnoise(1*pn)-0.0); 
    d += abs(an);
    pn *= ppn*rotSome;
    an *= aan; 
//    pn += 10*s-time*0.5; 
  }

  for (int i = md; i < mx; ++i) {
    s += an*1*(vnoise(1*pn)-0.0); 
    d += abs(an);
    pn *= ppn*rotSome;
    an *= aan; 
    pn += (3*(i+1))*s-time*5.5; 
  }

  s /= d;

  diff = (pn - opn);

  return s;
}

vec4 plane_3(vec3 ro, vec3 rd, vec3 pp, float aa, float n) {
  vec2 p = pp.xy;
  float z = pp.z;
  float nz = pp.z-ro.z;
  
  vec2 diff;
  vec2 hp = p;
  rot(hp, -0.3);
  hp -=0.125*vec2(nz, 0.0);
  hp -= 2.3*n;
  hp *= 0.6;
  float h = twirlyHeight(hp, n, diff);
  float gh = vnoise(hp);
  
  h = abs(h);
  float s = h;
  vec3 col = vec3(0.0);
  col = vec3(0.0);
  vec2 pdiff = toPolar(diff);
  float huen = fract(pdiff.x/200.0+gh);
  float satn = 1.0;
  float brin = s;
  col = hsv2rgb(vec3(huen, satn, brin));
  
  float t = sqrt(s)*(smoothstep(0.0, 0.5, length(pp - ro)))*smoothstep(0.0, mix(0.4, 0.75, pow(psin(TIME*TAU*BPM/60.0), 4.0)), length(p));
  return vec4(col, t);
}


float ofbm(vec2 p) {
  const float aa = 0.6;
  const float ff = 2.08;
  float a = 1.0;
  float f = 1.0;
  float s = 0.0;
  float d = 0.0;

  const float ss = 0.5;
  float toff = +TIME*0.1;
  p *= 1.0;
  for (int i = 0; i < 4; ++i) {
    float n = onoise(f*p);
    s += a*(n);
    d += abs(a);
    a *= aa;
    f *= ff;
    p *= rotSome;
    p += ss*s;
    p += toff;
  }
  
  return (0.5 + 0.5*(s/d));
}

float oheight(vec2 p, float hm) {
  p *= 5.0;
  return ofbm(p)*hm;
}

vec3 onormal(vec2 p, float hm, float aa) {
  vec2 eps = vec2(aa*3.0, 0.0);
  
  vec3 n;
   
  n.x = oheight(p - eps.xy, hm) - oheight(p + eps.xy, hm);
  n.y = oheight(p - eps.yx, hm) - oheight(p + eps.yx, hm);
  n.z = 2.0*eps.x;
  
  return normalize(n);
}

vec4 plane_4(vec3 ro, vec3 rd, vec3 pp, float aa, float n) {
  float r = hash(n);
  vec2 p = pp.xy;
  float z = pp.z;
  float nz = pp.z-ro.z;
  
  float l = length(p);

  p -= (1+0*(pp.z - ro.z))*offset(z).xy;
  vec2 pps = toPolar(p);
  pps.y += PI*pps.x-0.25*TIME;
  p = toRect(pps);

  vec2 ps = p;
//  vec2 pps = toPolar(ps);
//  ps += 0.0*offset(z);
  float ds = starn(ps, 0.25, 9, mix(3.0, 9.0, tanh(0.25*nz)))-0.0;
  float hms = 1.0-exp(-10*ds);
  float ts = smoothstep(-aa*2, aa*2, ds);

  vec2 pn = p;
  pn += r*10;
  rot(pn, r*TAU);
  float hm = exp(-6*l*l);
  float h = oheight(pn, hm);
  vec3 nor = onormal(pn, hm, aa);

  vec3 lp = ro + 10*vec3(1.0, 1.0, 0.0);
  vec3 p3 = pp + vec3(0.0, 0.0, -h);
//  lp = p3 + 2*vec3(-2.0, -2.0, 1);
  vec3 ld = normalize(lp - p3);

  float diff = max(dot(nor, ld), 0.0);

  vec3 col = vec3(0.0);
  float edge = abs(1.0 - nor.z);
  float flatness = max(dot(nor, vec3(0.0, 0.0, 1.0)), 0.0);
  col += pow(vec3(diff), 1.0*vec3(15.0))+pow(h, 1.5)*vec3(1.0, 0.5, 0.5);
  col += pow(vec3(edge), 8.0*vec3(0.25, 0.75, 1.25));
  col += pow(vec3(edge), 10.0*vec3(1.0));
  col = tanh(col);
  vec3 hsv = rgb2hsv(col);
  hsv.x+=0.25*nz+0.5*l+r;
  hsv.y *= hms;
  hsv.z *= pow(1.0-flatness, 1.0);
  hsv.z = mix(1.0, hsv.z, hms);
  hsv.z *= ts;
  hsv = fract(hsv);
  col = hsv2rgb(hsv);
  
  float t = ts*mix(1.0, 1.0-0.95*exp(-2*(1.0-flatness)), hms*hms);

  vec4 fcol = vec4(col, t);
//  fcol.xyz += vec3(1.0)*step(ds, 0.0);
  fcol.w += 1.5*(1.0-smoothstep(0.0, 0.25, sqrt(-ds)))*step(ds, 0.0);

  return clamp(fcol, 0.0, 1.0);
}

vec4 plane_5(vec3 ro, vec3 rd, vec3 pp, float aa, float n) {
  float r = hash(n);
  vec2 p = pp.xy;
  float z = pp.z;
  vec2 op = p+sin(n)*1;
  float nz = pp.z-ro.z;
  p -= (1+2*(pp.z - ro.z))*offset(z).xy;

  vec2 ppp = toPolar(p);
//  ppp.y += ppp.x*0.5+n+nz*sqrt(r);
  ppp.y += 0.25*pow(nz, 1.5);
  ppp.x *= (-1.0+0.5*length(p));
  p = toRect(ppp);

  const float st = 5.0;
  vec2 dd = df(p*st, nz, n, nz);
  dd /= st;
  float d = dd.x;
  
  vec4 plainCol = vec4(1.0);
  plainCol.xyz = mix(plainCol.xyz, vec3(0.0), smoothstep(-aa, aa, -d));
  plainCol.w = smoothstep(-aa, aa, dd.y);
  
  d = abs(d)-0.005;
  d = abs(d)-0.0025;
  
  float df = 0.8*pow(1.0/(nz+1), 0.5);
  
  vec3 glowHsv = mix(vec3(0.0, 0.5*df, 1.0), vec3(1.0, 1.0*df, 1.0), 0.5 + 0.5*sin(-1.0*TIME+0.25*TAU*length(op)));
  vec3 glowCol = hsv2rgb(glowHsv);

  vec4 col = vec4(0.0);
  col = vec4(glowCol, exp(-d*30*df*df));  
  col.xyz += 0.5*(glowCol+vec3(0.2))*exp(-d*800*df*df);  
//  col.xyz = 1.0-col.xyz;
  float cc = sin(nz+-TIME*3+1*p.x*p.y);
  float c = 0.5 + 0.5*(sign(cc)*pow(abs(cc), 0.25));
  col = mix(col, plainCol, c);
  return col;
}

float noise1(vec2 p) {
  float s = 1.0;

  p *= tanh(0.1*length(p));

  float a = cos(p.x);
  float b = cos(p.y);

  float c = cos(p.x*sqrt(3.5));
  float d = cos(p.y*sqrt(1.5));

  return a*b*c*d;
}

float galaxy(vec2 p, float a, float z) {
  vec2 pp = toPolar(p);
  pp.y += pp.x*3.0 + a;
  p = toRect(pp);
  
  p *= z;
  
  return noise1(p);
}

float height(vec2 p) {
  if (length(p) > 1.1) return 0.0;
  
  p *= 0.75;
  float s = 0.0;
  float a = 1.0;
  float f = 15.0;
  float d = 0.0;
  for (int i = 0; i < 3; ++i) {
    float g = a*galaxy(p, -0.5*TIME + TIME*(0.2*pow(float(i), 1.0)), f);
    s += g;
    a *= sqrt(0.5);
    f *= mix(sqrt(2.0), s, psin(TIME));
    d += a;
  }
  
  float l = length(p);
  
  return SABS((-0.25+ s/d), 0.5)*exp(-5.5*l*l);
}

vec3 normal(vec2 p, float aa) {
  vec2 eps = vec2(aa*5, 0.0);
  
  vec3 n;
  
  n.x = height(p - eps.xy) - height(p + eps.xy);
  n.y = 2.0*eps.x;
  n.z = height(p - eps.yx) - height(p + eps.yx);
  
  return normalize(n);
}

vec4 plane_6(vec3 ro, vec3 rd, vec3 pp, float aa, float n) {
  float ss = 3.0;
  float rr = hash(n);
  vec2 p = pp.xy/ss;
  vec2 op = pp.xy;
  vec2 ppp = toPolar(p);
//  ppp.x *= mix(1.0, (-1.0+ppp.x*2), psin(TIME+length(op)+TAU*rr));
  p = toRect(ppp);
  p += fract(vec2(rr, 23.0*rr))*1232;
  vec2 hn = hextile(p);
  p *= ss;
  float r = hash(hn+n);
  p += (-1.0 + 2.0*fract(vec2(r, 23.0*r))*ss/4.0)*0;
  vec3 lp1 = vec3(-2.0, 0.5, 2.0);
  vec3 ld1 = normalize(lp1 - vec3(p.x, 0.0, p.y));
  vec3 lp2 = vec3(2.0, 0.5, 2.0);
  vec3 ld2 = normalize(lp2 - vec3(p.x, 0.0, p.y));

  float l = length(p);
  float tl = tanh(l);

  float h = height(p);
  float th = tanh(h);
  vec3 nn = normal(p, aa);
  
  float diff1 = max(dot(ld1, nn), 0.0);
  float diff2 = max(dot(ld2, nn), 0.0);

  const vec3 col1 = vec3(0.5, 0.5, 0.75).zyx;
  const vec3 col2 = vec3(0.5, 1.0, 1.5);
  const vec3 col3 = vec3(0.5, 1.0, 1.5).zyx;
  vec3 col = vec3(0.0);
  col += col1*h;
  col += 0.5*pow(diff1, 10.0);
  col += 0.25*pow(diff2, 2.0);
  col += pow(vec3(0.5)*h, nn.y*1.75*(mix(col2, col3, 1.25*tl)));
  col += 0.9*vec3(1.0, 0.9, 0.75)*exp(-10.0*l*l);

  float dd = length(ro - pp)*0.9;

  return vec4(1*col/(1.0+0.2*dd*dd), tanh(10*h));
}

vec4 plane(vec3 ro, vec3 rd, vec3 pp, float aa, float n) {
  vec4 col = vec4(0.0);
/*
  int mn = int(mod(n/6.0, 3.0));
  switch(mn) {
    case 0: col = plane_0(ro, rd, pp, aa, n); break;
    case 1: col = plane_1(ro, rd, pp, aa, n); break;
    case 2: col = plane_2(ro, rd, pp, aa, n); break;
    default: col = vec4(1.0, 0.0, 0.0, 0.5); break;
  }
  */
  col = plane_6(ro, rd, pp, aa, n);
  
  return col;
}

vec3 skyColor(vec3 ro, vec3 rd) {
  return mix(1.5*vec3(0.75, 0.75, 1.0), vec3(0.0), tanh(length(3.5*rd.xy)));
//  return mix(1.5*vec3(0.75, 0.75, 1.0), vec3(0.0), length(2.0*rd.xy));
  return vec3(0);
}

vec3 color(vec3 ww, vec3 uu, vec3 vv, vec3 ro, vec2 p) {
  float lp = length(p);
//  vec3 rd = normalize(p.x*uu + p.y*vv + (2.0)*ww);
  vec3 rd = normalize(p.x*uu + p.y*vv + (2.0+0.5*tanh(lp))*ww);
//  vec3 rd = normalize(p.x*uu + p.y*vv + (2.0-1.0*tanh(lp))*ww);
// Warp off!
//  float warpf = smoothstep(5.0, 20.0, TIME);
//  float planeDist = mix(6.0, 0.1, warpf);
//  vec3 rd = normalize(p.x*uu + p.y*vv + (2.0+mix(0.0, 15.0, warpf)*tanh(lp))*ww);

//  vec3 rd = normalize(p.x*uu + p.y*vv + (3.0 + 2*tanh(lp)- 0.75*tanh(lp)*sin(10*p.x+time)*sin(19*p.y+time))*ww);
//  vec3 rd = normalize(p.x*uu + p.y*vv + (3.0+1.0*tanh(lp)*sin(lp*3+2*TIME))*ww);
  
  const vec3 errorCol = vec3(1.0, 0.0, 0.0);

  float planeDist = 1.0;
  const int furthest = 6;
  const int fadeFrom = max(furthest-4, 0);

  float nz = floor(ro.z / planeDist);

  vec3 skyCol = skyColor(ro, rd);  
  
  vec3 col = skyCol;

  for (int i = furthest; i >= 1 ; --i) {
    float pz = planeDist*nz + planeDist*i;
    
    float pd = (pz - ro.z)/rd.z;
    
    if (pd > 0.0) {
      vec3 pp = ro + rd*pd;
   
      float aa = length(dFdy(pp));

      vec4 pcol = plane(ro, rd, pp, aa, nz+i);
      float nz = pp.z-ro.z;
      float fadeIn = (1.0-smoothstep(planeDist*fadeFrom, planeDist*furthest, nz));
      float fadeOut = smoothstep(0.0, planeDist*0.1, nz);
      pcol.xyz = mix(skyCol, pcol.xyz, (fadeIn));
//      pcol.w *= sqrt(fadeIn);
      pcol.w *= fadeOut;

#if 1
      col = alphaBlendGamma(col, pcol, std_gamma);
#else      
      col = alphaBlend(col, pcol);
#endif      
    } else {
      break;
    }
    
  }

//  col =  max(pow(col, 1.0/std_gamma), 0.0) ;
  return col;
}

vec3 postProcess(vec3 col, vec2 q)  {
  col=pow(clamp(col,0.0,1.0),vec3(0.75)); 
  col=col*0.6+0.4*col*col*(3.0-2.0*col);
  col=mix(col, vec3(dot(col, vec3(0.33))), -0.4);
  col*=0.5+0.5*pow(19.0*q.x*q.y*(1.0-q.x)*(1.0-q.y),0.7);
  return col;
}

vec3 effect(vec2 p, vec2 q) {
  float tm = time*0.5;
  vec3 ro   = offset(tm);
  vec3 dro  = doffset(tm);
  vec3 ddro = ddoffset(tm);

  vec3 ww = normalize(dro);
  vec3 uu = normalize(cross(normalize(vec3(0.0,1.0,0.0)+ddro), ww));
  vec3 vv = normalize(cross(ww, uu));
  
  vec3 col = color(ww, uu, vv, ro, p);
  col = postProcess(col, q);
  return col;
}


void main(void) {
  vec2 q = inData.v_texcoord;
  vec2 p = -1. + 2. * q;
  p.x *= RESOLUTION.x/RESOLUTION.y;
  
  vec3 col = effect(p, q);
  
  fragColor = vec4(col, 1.0);
}