in  vec4 dir;
out vec4 fragColor;
uniform vec4  sparam;
uniform vec4  param[2];
uniform vec4 camera; // xyz- camera origin(LY)  , w - 1/rad(LY)
uniform sampler2D tex0;
#define  iChannel0 tex0
float _time;

const float radius=2.8;
const vec3 star_pos=vec3(0.0);
const float star_power=1.0;

//vec4 param=vec4(1000.,2500.,2.0e7,0.1); // protoplanet nebula
// vec4 param=vec4(1500.,2000.,3.0e7,0.065); // supernova


#define   T0         param[0].x
#define   Ts         param[0].y
#define   Exp        param[0].z
#define  dust_level  param[0].w

#define  age         param[1].x
#define  seed        param[1].y

vec3 firePalette(float i)
{
    float T = T0 + Ts*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(-Exp/L)); // Exposure level. Set to "50." For "70," change the "5" to a "7," etc.
}

//--------------------------------------------------------------

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

float fbm( vec3 p )
{
   return noise(p*.06125)*.5 + noise(p*.125)*.25 + noise(p*.25)*.125 + noise(p*.4)*.2;
}

float Sphere( vec3 p, float r )
{
    return length(p)-r;
}

const float ilife=1./10.;

//==============================================================
// otaviogood's noise from https://www.shadertoy.com/view/ld2SzK
//--------------------------------------------------------------
// This spiral noise works by successively adding and rotating sin waves while increasing frequency.
// It should work the same on all computers since it's not based on a hash function like some other noises.
// It can be much faster than other noise functions if you're ok with some repetition.
const float nudge = 4.; // size of perpendicular vector
float normalizer = 1.0 / sqrt(1.0 + nudge*nudge); // pythagorean theorem on that perpendicular to maintain scale
float SpiralNoiseC(vec3 p)
{
    float n = 1.-_time; // noise amount
    float iter = 2.0;
    for (int i = 0; i < 8; i++)
    {
        // add sin and cos scaled inverse with the frequency
        n += -abs(sin(p.y*iter) + cos(p.x*iter)) / iter;  // abs for a ridged look
        // rotate by adding perpendicular and scaling down
        p.xy += vec2(p.y, -p.x) * nudge;
        p.xy *= normalizer;
        // rotate on other axis
        p.xz += vec2(p.z, -p.x) * nudge;
        p.xz *= normalizer;
        // increase the frequency
        iter *= 1.733733;
    }
    return n;
}

float VolumetricExplosion(vec3 p)
{
    float final = Sphere(p,4.);
    final += noise(p*20.)*.4;
    final += SpiralNoiseC(p.zxy*fbm(p*10.))*2.5; //1.25;

    return final;
}

//--------------------------------------------------------------


bool RaySphereIntersect(vec3 org, vec3 dir, out float near, out float far)
{
  float b = dot(dir, org);
  float c = dot(org, org) - radius*radius;
  float delta = b*b - c;
  if( delta < 0.0)
    return false;
  float deltasqrt = sqrt(delta);
  near = -b - deltasqrt;
  far = -b + deltasqrt;
  return far > 0.0;
}


void main()
{
 vec3 ro     = camera.xyz*camera.w*radius;
 vec3 rd     = normalize(dir.xyz-camera.xyz);
 float len_ro= length(ro);
 if (sparam.y>=0.0)
 {
  _time = min(10.0,sqrt(max(0.01,sparam.y-2.))*2.);
  ro*=(len_ro+100.-_time*10.)/len_ro;
 }
 else                 _time = age;

 float ld=0., td=0., w=0.;
 float d=1., t=0.;
 const float h = 0.1;
 vec4 sum = vec4(0.0);
 float min_dist, max_dist;
 if(!RaySphereIntersect(ro, rd, min_dist, max_dist))  discard; // { fragColor = vec4(0.1); return; }
 _time=mod(_time * ilife,-1.);

 t = max(0.0,min_dist);
 float df = 0.8+0.08*texture(iChannel0,rd.xy*vec2(280.,270.)).r;
 for (int i=0; i<86; i++)
 {
   vec3 pos = ro + t*rd;
   if (td>0.9 || sum.a > 0.99 || t>max_dist)  break;
   if ( t>2.*radius+min_dist) break;
   float d = VolumetricExplosion(pos/(1.+_time))*(1.+_time); // scale
   d=abs(d)+dust_level;
   vec3 ldst = star_pos-pos;
   float lDist = clamp(length(ldst)/star_power, 0.0001,2.0);
   vec3 lightColor=firePalette(1.-lDist*.25);
   sum.rgb+=(vec3(0.67,0.75,1.00)/(lDist*lDist*1.e4)); // star itself
   sum.rgb+=(lightColor/exp(lDist*lDist*lDist*.08)/30.); // bloom
   if (d<h)
   {
     ld = h - d;
     w = (1. - td) * ld;
     td += w + 1./200.;
     sum += sum.a * vec4(sum.rgb, 0.0) * 0.2 / lDist;
     sum.a +=td*0.2*(1.0 - sum.a);
   }
   td += 1./70.;
   d=abs(d)*df;
   t += max(d * 0.1 * max(min(length(ldst),len_ro),0.1), 0.01);
 }
 sum *= 1. / exp( ld * 0.2 ) * 0.8;
 sum = clamp( sum, 0.0, 1.0 );
 sum.xyz = sum.xyz*sum.xyz*(3.0-2.0*sum.xyz);
 fragColor = vec4(sum.xyz*sum.xyz*dir.w,0.0);
}
