in  vec4 dir;
out vec4 fragColor;
uniform vec4 sparam;
uniform vec4 camera; // xyz- camera origin(LY)  , w - 1/object_rad(LY)
uniform sampler2D tex0;
#define  iChannel0 tex0
#define  iGlobalTime sparam.x

const float jet_len=3.0,jet_rad=0.2,jscale=0.04,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); }


#if 1 // HIGH_QUALITY
#define STEPS 100
#define ALPHA_WEIGHT 0.025
#define BASE_STEP 0.05
#else
#define STEPS 50
#define ALPHA_WEIGHT 0.05
#define BASE_STEP 0.1
#endif

const float freq=1.3;
float K;
float hash21(in vec2 n){ return fract(sin(dot(n, vec2(12.9898, 4.1414))) * 43758.5453); }
float noise(in vec3 p) //iq's ubiquitous 3d noise
{
	vec3 ip = floor(p), f = fract(p);
    #ifdef HIGH_QUALITY
    f = f*f*f*(f*(f*6. - 15.) + 10.); //Quintic smoothing
    #else
    f = f*f*(3.0-2.0*f); //Cubic smoothing
    #endif
	vec2 uv = (ip.xy+vec2(37.0,17.0)*ip.z) + f.xy;
	vec2 rg = texture2D( iChannel0, (uv+ 0.5)/256.0, -100.0 ).yx;
	return mix(rg.x, rg.y, f.z);
}

float fbm(in vec3 p, in float sr)
{
    p *= 3.5;
    float rz = 0., z = 1.;
    for(int i=0;i<4;i++)
    {
        float n = noise(p-iGlobalTime*.5);
        rz += (sin(n*4.4)-.45)*z;
        z *= .5;
        p *= 3.;
    }
    return rz;
}

vec4 map(in vec3 p)
{
    float dtp = dot(p,p);
    p = .5*p/(dtp + .2);
    p.xy = R(p.xy, p.x*3.);
    p.zx = R(p.zx, p.x*4.);

    float dtp2 = dot(p, p);
    p = (freq + .6)*3.*p/(dtp2 - 5.);
    float r = clamp(fbm(p, dtp*0.1)*1.5-dtp*K, 0. ,1.);
    float r2=r*r;
    return vec4(r2*r2*.6,r2*r*.8,r,r*.5);
}

//------------------------- Ray-Sphere Intersection --------------------------
bool RSI(in vec3 ro,in vec3 rd,float rad,out vec2 res)
{
  float b = dot(ro,rd);
  float c = dot(ro,ro) - rad*rad;
  float d = b*b - c;
  if(d <= 0.0) return false;
  d = sqrt(d);
  res.x=(-b - d);
  res.y=(-b + d);
  return true;
}



vec4 vmarch(in vec3 ro, in vec3 rd)
{
 vec2  ts;
 vec4 rz = vec4(0);
 if (RSI(ro,rd,2.5,ts))
 {
  float t = ts.x + 0.03*hash21(gl_FragCoord.xy);
  K=.35-sin(iGlobalTime)*0.15;
	for(int i=0; i<STEPS; i++)
	{
    if(rz.a > 0.99 || t>=ts.y )break;
		vec3 pos = ro + t*rd;
    vec4 col = map(pos);
    float den = col.a;
    col.a *= ALPHA_WEIGHT;
		col.rgb *= col.a*1.7;
		rz += col*(1. - rz.a);
    t += BASE_STEP - den*(BASE_STEP-BASE_STEP*0.015);
	}
 }
 return rz;
}



//------------------------- 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;
}



//-------------------------------------------
float jet(vec3 p)
{
 return length(vec3(p.xy*(1.0/jet_rad),abs(p.z)*(.5/jet_len)));
}

 // http://adsabs.harvard.edu/abs/1987A&A...175....1M
float bulge(vec3 p)
{
 float  r=length(p);
 return clamp(pow(r,-0.855)*exp(-pow(r,.25))*.1-.05,0.,1.)*.7;
}

//-------------------------Main---------------------------
void main( )
{
 vec3 ro  =  camera.xyz*camera.w*jet_len;
 vec3 rd  =  normalize(dir.xyz-camera.xyz);
 vec2  ts,tj;
 float cj=0.0,cs=0.;
 float dist=length(ro)*.02;
 float T=dist>1.? 5.0:0.1;
 if (RSI(ro,rd,1.,ts))
 {
   ts.x=max(ts.x,0.0);
   for (int n=0;n<100;n++)
   {
    if (cs>1.||ts.x>=ts.y) break;
    float b=bulge(ro+ts.x*rd);
    cs += b*b;
    ts.x+=.1;
   }
 }
 float a  = fract(iGlobalTime*T)*pi*2.;
 float aa =  0.3*pi;
 rd.xy = R(rd.xy, a );
 ro.xy = R(ro.xy, a );
 rd.yz = R(rd.yz, aa);
 ro.yz = R(ro.yz, aa);
 if (RCI(ro,rd,vec2(jet_rad,jet_len),tj))
 {
  tj.x=max(tj.x,0.0);
  for (int n=0;n<100;n++)
  {
   if (cj>1./jscale||tj.x>=tj.y) break;
   cj+=clamp(1.0-jet(ro+tj.x*rd),0.,1.);
   tj.x+=0.02;
  }
 }
 vec3 col = clamp(vec3(9., 9.,10.)*cs,0.,1.) +
            vec3(.8,.8,1.) * (cj*cj*jscale*jscale);
 if (dist<1.0) col += vmarch(ro,rd).rgb*(1.0-dist);
 fragColor = vec4(col*col*dir.w, 0.0);
}
