 // Based on "Binary system" by Duke, www.shadertoy.com/view/Msd3DH
in  vec4 dir;
out vec4 fragColor;
uniform vec4  sparam;
uniform vec4 camera; // xyz- camera origin(LY)  , w - 1/galaxy_rad(LY)
uniform sampler2D tex0;
#define  iChannel0 tex0
#define  iGlobalTime sparam.x

const float jet_len=1.0,jet_rad=0.05,jscale=0.2,disk_rad=1.0,disk_thickness=0.2;
#define pi 3.14159265

vec2  R(vec2 p, float a) { return cos(a)*p+sin(a)*vec2(p.y, -p.x); }
//-------------------------Color palette------------------
// See "Combustible Voronoi"
// https://www.shadertoy.com/view/4tlSzl
vec3 firePalette(float i)
{
    float T = 1550. + 1400.*i; // Temperature range (in Kelvin).
    vec3 L = vec3(7.4, 5.6, 4.4); // Red, green, blue wavelengths (in hundreds of nanometers).
    L = pow(L,vec3(5.0)) * (exp(1.43876719683e5/(T*L))-1.0);
    return 1.0-exp(-5e8/L); // Exposure level. Set to "50." For "70," change the "5" to a "7," etc.
}

//------------------------- Ray-Cylinder Intersection --------------------------
bool RCI(in vec3 ro,const vec3 rd,const vec2 size/*w,h*/,out vec2 res)
{
  ro/=size.x;
  float a = dot(rd.xy,rd.xy);
  float b = dot(rd.xy,ro.xy)*2.;
  float c = dot(ro.xy,ro.xy)-1.;
  float d = b*b - 4.*a*c;
  if(d <= 0.0) return false;
  float h=size.y/size.x;
  d     = sqrt(d);
  res.x = (-b - d)/(2.*a);
  res.y = (-b + d)/(2.*a);
  vec2 z=ro.z+rd.z*res;
  if (z.x>h)  { if (z.y>h) return false;  res.x+=(res.y-res.x)*(h-z.x)/(z.y-z.x); z.x=h; } else
  if (z.y>h)  {                           res.y+=(res.x-res.y)*(h-z.y)/(z.x-z.y); z.y=h; }
  if (z.x<-h) { if (z.y<-h) return false; res.x+=(res.y-res.x)*(-h-z.x)/(z.y-z.x); } else
  if (z.y<-h) {                           res.y+=(res.x-res.y)*(-h-z.y)/(z.x-z.y); }
  res*=size.x;
  return true;
}

mat2 Spin(float angle){  return mat2(cos(angle),-sin(angle),sin(angle),cos(angle)); }

float smin( 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);
}

//-------------------------Noise--------------------------
// IQ's noise
float pn( in vec3 p )
{
    vec3 ip = floor(p);
    p = fract(p);
    p *= p*(3.0-2.0*p);
    vec2 uv = (ip.xy+vec2(37.0,17.0)*ip.z) + p.xy;
    uv = texture( iChannel0, (uv+ 0.5)/256.0, -100.0 ).yx;
    return mix( uv.x, uv.y, p.z );
}

// FBM
float fpn(vec3 p) {
    return pn(p*.06125)*.57 + pn(p*.125)*.28 + pn(p*.25)*.15;
}

float rand(vec2 co){// implementation found at: lumina.sourceforge.net/Tutorials/Noise.html
  return fract(sin(dot(co*0.123,vec2(12.9898,78.233))) * 43758.5453);
}

//-------------------------------------------
float jet(vec3 p)
{
 float z=abs(p.z)/jet_len;
 return length(vec3(p.xy*(1.0/jet_rad),z*z)) + pn(vec3(p.xy*80.,(p.z*40.-iGlobalTime*sign(p.z)*10.)))*0.1;
}

//-------------------------------------------
float disk(vec3 p)
{
 float t=iGlobalTime*.05;
 p.xy=R(p.xy, pi*iGlobalTime*0.02);
 //---------------------------
 float DiskD = fpn(vec3(Spin(t*0.25+p.z*.10)*p.xy*40.,p.z*20.-t)*5.0);
 float r = length(p.xy)*.25;
 float a = atan(p.y, p.x);
 float dOffset = cos( (a+r*r*120.0)*2.-iGlobalTime*2.) * .5;
 dOffset *= 1.0-r;
 return DiskD + dOffset*.45;
}

//-------------------------Main---------------------------
void main()
{
 vec3 ro  =  camera.xzy*camera.w*jet_len;
 vec3 rd  =  normalize(dir.xzy-camera.xzy);
 vec2 td,tj;
 bool fd=RCI(ro,rd,vec2(disk_rad,disk_thickness),td);
 bool fj=RCI(ro,rd,vec2(jet_rad,jet_len),tj);
 float cj=0.0,cd=0.;
 vec3 col=vec3(0.0);
 td+=0.02*rand(rd.xy + fract(iGlobalTime)); // dithering
 if (fj)
 {
  float te = fd ? min(td.x,tj.y) : tj.y;
  tj.x=max(tj.x,0.0);
  for (int n=0;n<100;n++)
  {
   if (cj>1./jscale||tj.x>=te) break;
   cj+=clamp(1.0-jet(ro+tj.x*rd),0.,1.);
   tj.x+=0.02;
  }
 }
 float shaded=0.;
 if (fd)
 {
  td.x=max(td.x,0.0);
  float ln = (td.y-td.x+1e-16);
  float dtd=.01*ln;
  for (int n=0;n<100;n++)
  {
   if (shaded>=1.0) break;
   vec3   p   = ro+td.x*rd;
   float dist =  disk(p);
   float dxy  = length(p.xy)/disk_rad;
   float dz   = pow(abs(p.z)/disk_thickness,.1);
   float fade = clamp((1.-dxy)*(1.-dz),0.,1.)*2.;
   float  Ld  = clamp(1. - dist,0.0,1.)*.2;
   cd += (1.-shaded) * (Ld*fade*fade);
   shaded += Ld*fade;
   td.x+=dtd;
  }
  cd/=ln;
  if (fj)
  {
   for (int n=0;n<100;n++)
   {
    if (cj>1./jscale||tj.x>=tj.y) break;
    cj+=(1.-shaded) * clamp(1.0-jet(ro+tj.x*rd),0.,1.);
    tj.x+=0.02;
   }
  }
  col+=mix(firePalette(cd*3.)*smoothstep(0.0,.2,cd),vec3(.5),shaded*.4);
}
 col +=vec3(cj*cj*(jscale*jscale));
 fragColor = vec4(col*col*dir.w, 0.0);
}
