//Spring FauxBoxBulbs by eiffie

//#define EUCLIDEAN
//#define TAXICAB
#define CHEBYSHEV

#ifdef EUCLIDEAN
 #define Pi 3.14159
 #define Cos cos
 #define Atan atan
 #define Length length
#endif
#ifdef TAXICAB
 #define Pi (2.0*sqrt(2.0))
 #define Cos t_cos
 #define Atan t_atan
 #define Length t_length
#endif
#ifdef CHEBYSHEV
 #define Pi 4.0
 #define Cos c_cos
 #define Atan c_atan
 #define Length c_length
#endif
float t_cos(float a){return 2.0*abs(mod(a,2.0*Pi)-Pi)/Pi-1.0;}
float c_cos(float a){return clamp(abs(mod(a,2.0*Pi)-Pi)-Pi/2.0,-1.0,1.0);}
float t_atan(float y, float x){//atan is always complicated by the quadrant (probably a simpler way to write these)
 float a=x-y,b=x+y,res;
 if(b==0.0)res=(a>0.0?7.0:3.0);
 float d=a/b;
 if(abs(d)<1.0){
  if(b>0.0)res=1.0-d;
  else res=5.0-d;
 }else {
  d=b/a;
  if(a>0.0)res=7.0+d;
  else res=3.0+d;
 }
 return res*0.25*Pi;
}
float c_atan(float y, float x){
 if(y==0.0)return (x>0.0?0.0:4.0);
 float a=x/y;
 if(abs(a)<1.0){
  if(y>0.0)return 2.0-a;
  else return 6.0-a;
 }else {
  a=y/x;
  if(x>0.0)return mod(a,8.0);
  else return 4.0+a;
 }
}
float t_length(vec2 p){return abs(p.x)+abs(p.y);}//==(x^1+y^1)^(1/1) -abs(x & y) is assumed for each
float t_length(vec3 p){return abs(p.x)+abs(p.y)+abs(p.z);}//==(x^1+y^1)^(1/1) -abs(x & y) is assumed for each
float c_length(vec2 p){return max(abs(p.x),abs(p.y));}//==(x^inf+y^inf)^(1/inf)
float c_length(vec3 p){return max(abs(p.x),max(abs(p.y),abs(p.z)));}//==(x^inf+y^inf)^(1/inf)


uniform float u[32];
float Sin(float a){return Cos(a-Pi*.5);}
vec3 cmap(float a){return fract(vec3(a*13.23,a*1.89,a*3.66));}
vec3 mcol=vec3(0.0);
float DE(vec3 p){
  float d=p.z,pr=mod(floor(p.x/3.)*floor(p.y/3.),6.),m=3.+2.*pr,s=m,a,b;
  p.z+=Cos(pr+iTime+length(p.xy)*.1);p.xy=mod(p.xy,3.0)-1.5;
  for(int i=0;i<3;i++){
    a=Atan(p.y,p.x)*s;b=Atan(Length(p.xy),p.z)*s+iTime*1.7+pr;
    p+=vec3(Cos(b)*vec2(Cos(a),Sin(a)),Sin(b))/s;
    s*=m;
  }
  float r=(Length(p)-(pr==0.?.8:1.))/3.;
  if(mcol.x>0.){if(d<r)mcol+=cmap(pr);else mcol+=abs(cmap(pr+1.)+.5*Cos((p.y+p.x+p.z)*50.))*.75+.25;}
  return min(d,r);
}
vec3 normal(vec3 p, float d){//from dr2
  vec2 e=vec2(d,-d);vec4 v=vec4(DE(p+e.xxx),DE(p+e.xyy),DE(p+e.yxy),DE(p+e.yyx));
  return normalize(2.*v.yzw+vec3(v.x-v.y-v.z-v.w));
}
vec3 sky(vec3 rd, vec3 L){
  float d=max(0.1,rd.z)+max(0.1,dot(rd,L));
  return vec3(d);
}
float rnd;
void randomize(in vec2 p){rnd=fract(float(iTime)+sin(dot(p,vec2(13.3145,117.7391)))*42317.7654321);}

float ShadAO(in vec3 ro, in vec3 rd){
 float t=0.01*rnd,s=1.0,d,mn=0.01;
 for(int i=0;i<12;i++){
  d=max(DE(ro+rd*t)*1.5,mn);
  s=min(s,d/t+t*0.5);
  t+=d;
 }
 return s;
}
vec3 scene(vec3 ro, vec3 rd){
  float t=DE(ro)*rnd,d,px=1.0/iResolution.x;
  for(int i=0;i<64;i++){
    t+=d=DE(ro+rd*t);
    if(t>100.0 || d<px*t)break;
  }
  vec3 L=normalize(vec3(0.4,0.25,0.5));
  vec3 col=sky(rd,L);
  if(d<px*t*5.0){
    mcol=vec3(0.001);
    vec3 so=ro+rd*t;
    vec3 N=normal(so,d);if(N!=N)N=-rd;if(dot(rd,N)>0.)N=-N;
    vec3 scol=mcol*0.25;
    float dif=0.5+0.5*dot(N,L);
    float vis=clamp(dot(N,-rd),0.05,1.0);
    float fr=pow(1.-vis,5.0);
    float shad=ShadAO(so,L);
    col=mix((scol*dif+shad*fr*sky(reflect(rd,N),L))*shad,col,t*t/10000.);
  }
  return col;
}
mat3 lookat(vec3 fw){fw=normalize(fw);vec3 rt=normalize(cross(fw,vec3(0,0,1)));
  return mat3(rt,cross(rt,fw),fw);
}
void mainImage(out vec4 O, in vec2 U){
  vec2 uv=vec2(U-0.5*iResolution.xy)/iResolution.x;
  randomize(U);
  float t=iTime*.05+5.;
  vec3 ro=vec3(t_cos(t)*100.+30.,t_cos(t*.3-2.)*150.+74.,10.+7.*t_cos(t));
  vec3 rd=lookat(vec3(30.,70.,-200.+10.*ro.z)-ro)*normalize(vec3(uv.xy,1.0));
  O=vec4(scene(ro,rd),1.0);
}