precision highp float;

#define lofi(i,j) (floor((i)/(j))*(j))
#define saturate(i) clamp(i,0.,1.)

// 200 is recommended, use 400 if your machine is a beef
const int ITER=400;

const float FOG_START=30.0;
const float FAR=40.0;
const vec3 FOG_COLOR=vec3(1,2,3);

const float PI=acos(-1.);
const float TAU=2.*PI;
const float SQRT3=sqrt(3.);

uniform vec2 resolution;
uniform float time;
uniform int frame;
uniform sampler2D backbuffer;

out vec4 outColor;

uvec3 seed;

uvec3 hash3u(uvec3 s){
  s=s*1145141919u+1919810u;
  s.x+=s.y*s.z;
  s.y+=s.z*s.x;
  s.z+=s.x*s.y;
  s^=s>>16;
  s.x+=s.y*s.z;
  s.y+=s.z*s.x;
  s.z+=s.x*s.y;
  return s;
}

vec3 hash3f(vec3 s){
  uvec3 r=hash3u(floatBitsToUint(s));
  return vec3(r)/float(0xffffffffu);
}

vec3 random3(){
  seed=hash3u(seed);
  return vec3(seed)/float(0xffffffffu);
}

float random(){
  return random3().x;
}

float easein(float x,float k){
  x=saturate(x);
  return (k+1.)*pow(x,k)-k*pow(x,k+1.);
}

float easeout(float x,float k){
  return 1.-easein(1.-x,k);
}

vec2 cis(float t){
  return vec2(cos(t),sin(t));
}

mat2 rotate2d(float t){
  float c=cos(t),s=sin(t);
  return mat2(c,s,-s,c);
}

mat3 orthbas(vec3 z){
  z=normalize(z);
  float isveryy=step(.999,abs(z.y));
  vec3 up=vec3(0,1.-isveryy,isveryy);
  vec3 x=normalize(cross(up,z));
  return mat3(x,cross(z,x),z);
}

vec2 dohex(vec2 p){
  mat2 MAT_SKEW=mat2(1,1./SQRT3,0,2./SQRT3);
  mat2 INV_MAT_SKEW=mat2(1,-.5,0,SQRT3/2.);

  vec2 pt=p*MAT_SKEW;

  vec2 cell=floor(pt);
  cell.y-=mod(cell.y+cell.x+2.,3.)-1.;

  vec2 coord=pt-cell;
  float istop=step(coord.x,coord.y);
  cell+=vec2(1.-istop,istop);

  return cell*INV_MAT_SKEW;
}

vec4 gridt(vec3 ro,vec3 rd){
  vec3 cell=vec3(dohex(ro.xy+rd.xy*1E-2),0);
  vec3 dice=hash3f(cell.xyy);

  float zoffset=dice.x;
  cell.z=lofi(ro.z+rd.z*1E-2-zoffset,1.)+.5+zoffset;

  mat4 dir=mat4(
    vec4(SQRT3*.5,.5,0,0),
    vec4(0,1,0,0),
    vec4(-SQRT3*.5,.5,0,0),
    vec4(0,0,1,0)
  );

  vec4 rdd=vec4(rd,0)*dir;
  vec4 src=-(vec4(ro-cell,0)*dir)/rdd;
  vec4 dst=abs(vec2(SQRT3/2.,.5).xxxy/rdd);

  vec4 bv=src+dst;
  float b=min(min(bv.x,bv.y),min(bv.z,bv.w));

  return vec4(cell,b);
}

vec4 isectmin(vec4 a,vec4 b){
  return mix(a,b,step(b.w,a.w));
}

vec4 isectbox(vec3 ro,vec3 rd,vec3 s){
  vec3 src=-ro/rd;
  vec3 dst=abs(s/rd);
  vec3 fv=src-dst;
  vec3 bv=src+dst;
  float f=max(max(fv.x,fv.y),fv.z);
  float b=min(min(bv.x,bv.y),bv.z);

  vec3 n=-sign(rd)*step(f,fv);

  float isvalid=step(f,b)*step(0.,f);
  return vec4(
    n,
    mix(FAR,f,isvalid)
  );
}

vec4 isecthex(vec3 ro,vec3 rd,float r,float l){
  mat4 dir=mat4(
    vec4(SQRT3*.5,.5,0,0),
    vec4(0,1,0,0),
    vec4(-SQRT3*.5,.5,0,0),
    vec4(0,0,1,0)
  );

  vec4 rdd=vec4(rd,0)*dir;
  vec4 src=-(vec4(ro,0)*dir)/rdd;
  vec4 dst=abs(vec2((SQRT3/2.)*r,l).xxxy/rdd);

  vec4 bv=src+dst;
  float b=min(min(bv.x,bv.y),min(bv.z,bv.w));

  vec4 fv=src-dst;
  float f=max(max(fv.x,fv.y),max(fv.z,fv.w));

  vec4 mask=step(f,fv);
  vec3 n=-(dir*mask).xyz*sign(dot(rdd,mask));

  float isvalid=step(f,b)*step(0.,f);
  return vec4(
    n,
    mix(FAR,f,isvalid)
  );
}

vec3 sampleggx(vec3 n,float aa){
  vec3 xi=random3();

  float costheta=sqrt((1.-xi.y)/(1.+(aa-1.)*xi.y));
  float sintheta=sqrt(1.-costheta*costheta);

  return orthbas(n)*vec3(
    sintheta*cis(TAU*xi.x),
    costheta
  );
}

void main() {
  seed=uvec3(gl_FragCoord.xy,frame);

  vec2 uv=gl_FragCoord.xy/resolution.xy;
  vec2 p=uv*2.-1.;
  p.x*=resolution.x/resolution.y;

  vec3 co;
  vec3 ro;
  vec3 rd;
  float rl;
  vec3 col=vec3(0);
  vec3 colrem=vec3(0);
  float samples=0.;

  for(int i=0;i<ITER;i++){
    const float gridzmod=256.;

    float mytime;

    if(dot(colrem,vec3(1))<.4){
      // reset
      float warptime=time/10.-0.8; // 3.8 for the sake of direction
      float warp1=easeout(6.*fract(warptime),6.);
      float warp0=easeout(6.*fract(warptime)-.1,6.);
      float warp=(warp1-warp0)/.1;
      warp=warp/(warp+1.);

      float prog=random();
      colrem=mix(
        vec3(1),
        4.*(.5-.5*cos(TAU*saturate(2.*prog-vec3(0,.5,1.)))),
        sqrt(warp)
      );
      mytime=time-.0166*prog;
      float warpt=mix(warp0,warp1,prog)+floor(warptime);
      mytime+=3.*warpt;

      samples++;

      // init camera
      vec3 cx=vec3(1,0,0);
      vec3 cy=vec3(0,1,0);
      cx.yz*=rotate2d(1.2*warpt);
      cy.yz*=rotate2d(1.2*warpt);
      cx.zx*=rotate2d(-.11*mytime);
      cy.zx*=rotate2d(-.11*mytime);
      mat3 cb=mat3(cx,cy,cross(cx,cy));

      co=cb*vec3(0,0,.6);
      co.z-=mod(4.*mytime,gridzmod);

      // init ray
      ro=co;
      vec2 pt=(p+2.*random3().xy/resolution.y)*rotate2d(.1*mytime+warpt);
      rd=cb*normalize(vec3(pt,-1.8+warp*(-.4-.2*length(p)*length(p))));
      rl=0.;

      // dof
      {
        vec3 fp=ro+rd*5.;
        vec3 dice=random3();
        vec2 disc=sqrt(dice.x)*cis(TAU*dice.y);
        ro+=0.01*vec3(disc,0)*cb;
        rd=normalize(fp-ro);
      }
    }

    // vars
    vec4 isect=vec4(FAR);
    float alpha=0.;
    float f0=.1;
    vec3 ecol=vec3(0);

    // scroll
    float scrdice=hash3f(dohex(ro.xy+rd.xy*1E-2).xyy).x;
    float scr=24.*smoothstep(-.2,.2,sin(.1*mytime+TAU*scrdice));
    ro.z-=scr;

    // traverse
    vec4 grid=gridt(ro,rd);

    // roll dice
    vec3 dice=hash3f(grid.xyy);
    dice.z=floor(mod(grid.z,gridzmod)*dice.z*dice.z);
    dice=hash3f(dice);

    // isect emissive hexes
    if(length(grid.xy)>1.&&dice.z<.1){
      float showz=grid.z-co.z+mix(0.,40.,dice.x);
      float size=easeout(showz/5.,5.)-easeout((showz-20.)/5.,3.);
      vec4 isect2=isecthex(ro-grid.xyz,rd,size*.7,size*.2);

      if(isect2.w<isect.w){
        isect=isect2;
        alpha=.001;
        f0=.1;
        ecol=4.*(.1+saturate(cos(PI*(fract(dice.z*100.)-vec3(.16,.5,.83)))));
      }
    }

    // lines
    if(length(grid.xy)>1.){
      float mt=PI/3.*floor(mod(dice.x*18.,6.));
      mat2 m=rotate2d(mt);

      vec3 rpt=ro-grid.xyz;
      rpt.xy*=m;

      vec3 rdt=rd;
      rdt.xy*=m;

      vec4 isect2=isectmin(
        isecthex(rpt-vec3(.72,.0,.0),rdt,.01,.49),
        isecthex(rpt-vec3(-.72,.0,.0),rdt,.01,.49)
      );
      isect2.xy*=rotate2d(-mt);

      if(isect2.w<isect.w){
        isect=isect2;
        alpha=1.;
        f0=.0;
        ecol=vec3(1.8,1.9,2);
      }
    }

    // reflectors
    if(length(grid.xy)>1.&&dice.z<.54){
      float mt=PI/3.*floor(dice.x*3.);
      mat2 m=rotate2d(mt);

      vec3 rpt=ro-grid.xyz;
      rpt.xy*=m;

      vec3 rdt=rd;
      rdt.xy*=m;

      vec4 isect2=isectmin(
        isectbox(rpt-vec3(0,0.8,0.),rdt,vec3(.5,0.02,.45)),
        isectbox(rpt-vec3(0,-0.8,0.),rdt,vec3(.5,0.02,.45))
      );
      bool isface=abs(isect2.y)>.5;
      isect2.xy*=rotate2d(-mt);

      if(isect2.w<isect.w){
        isect=isect2;
        alpha=isface?exp(-3.-3.*dice.y):1.;
        f0=isface?.5:.04;
        ecol=vec3(0);
      }
    }

    // revert the scroll
    ro.z+=scr;

    // if the ray misses, continue
    if(isect.w==FAR){
      ro+=grid.w*rd;
      rl+=grid.w;

      // if the ray goes way too far, grant fog instead
      if(rl>=FAR){
        col+=colrem*FOG_COLOR;
        colrem*=0.;
      }

      continue;
    }

    // update ray
    ro+=isect.w*rd;
    rl+=isect.w;

    // do the shading
    vec3 n=isect.xyz;

    // fog
    float fog=smoothstep(FOG_START,FAR,rl);
    col+=colrem*fog*FOG_COLOR;
    colrem*=1.-fog;

    // ggx
    vec3 h=sampleggx(n,alpha*alpha);
    float dotvh=dot(-rd,h);
    rd=reflect(rd,h);

    // fresnel
    float f=mix(f0,1.,pow(saturate(1.-dotvh),5.));

    // emissive
    col+=colrem*(1.-f)*ecol;

    // albedo
    colrem*=f*step(0.,dot(rd,n)); // fresnel, visibility
  }

  // post processing
  col/=samples; // sum -> average
  col=pow(col/(col+1.),vec3(.4545)); // reinhard, oetf
  // col*=1.-.2*length(p)*length(p); // vignette
  col=smoothstep(vec3(0.0,-0.1,-0.3),vec3(1.0,1.0,1.1),col); // color grading

  // temporal accumulation
  vec4 prev=texture(backbuffer,uv);
  outColor=mix(vec4(col,1),prev,.5*prev.w);

  // debug
  // outColor.xyz=cos(PI*(saturate(samples/10.)-vec3(.5,1.5,2.5)/3.)); // blue == enough samples, red == too few samples
}