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

#define pi 3.14159265
#define R(p, a) p=cos(a)*p+sin(a)*vec2(p.y, -p.x)

const float spin_speed=0.2;
mat2 Spin(float angle){
  return mat2(cos(angle),-sin(angle),sin(angle),cos(angle));
}

//vec4 param=vec4(1500.,1400.,9.0e8,30.);
//vec4 param=vec4(2500.,3500.,6.0e6,55.);

// #define   T0         param.x
// #define   Ts         param.y
// #define   Exp        param.z
// #define  sun_density param.w

#define   T0          1500.
#define   Ts          1400.
#define   Exp         9.0e8
#define  sun_density  30.






// -- Common Tools ----------------------------------------------------

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.
}



//-------------------------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 = texture2D( 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);
}

//-------------------------Disk---------------------------

float length2( vec2 p )
{
  return sqrt( p.x*p.x + p.y*p.y );
}

float length8( vec2 p )
{
  p = p*p; p = p*p; p = p*p;
  return pow( p.x + p.y, 1.0/8.0 );
}

float Disk( vec3 p, vec3 t )
{
    vec2 q = vec2(length2(p.xy)-t.x,p.z*0.5);
    return max(length8(q)-t.y, abs(p.z) - t.z);
}

//-----------------------------------------------------
float hash(in float n) { return fract(sin(n)*43758.5453123); }
float noise(in vec3 x) {
    vec3 p = floor(x), f = fract(x);
  f = f*f*(3.-2.*f);
  vec2 rg = texture2D(iChannel0, (((p.xy+vec2(37.,17.)*p.z) + f.xy)+.5)/256., -100.).yx;
  return mix(rg.x, rg.y, f.z);
}
const mat3 msun = mat3(0., .8, .6, -.8, .36, -.48, -.6, -.48, .64);

float smoothNoise(in vec3 q)
{
  float   f  = .5000*noise(q); q=msun*q*2.01;
          f += .2500*noise(q); q=msun*q*2.02;
          f += .1250*noise(q);
  return f;
}

//-------------------------Star---------------------------
// adaptation of "Pyroclastic sun" https://www.shadertoy.com/view/XtlXR2
// Noise Settings
#define NoiseSteps 4
#define NoiseAmplitude 0.032

float Turbulence(vec3 position, float minFreq, float maxFreq, float qWidth)
{
    float value = 0.0;
    float cutoff = clamp(0.5/qWidth, 0.0, maxFreq);
    float fade;
    float fOut = minFreq;
    for(int i=NoiseSteps ; i>=0 ; i--)
    {
      if(fOut >= 0.5 * cutoff) break;
      fOut *= 2.0;
      value += abs(pn(position * fOut))/fOut;
    }
    fade = clamp(2.0 * (cutoff-fOut)/cutoff, 0.0, 1.0);
    value += fade * abs(fpn(position * fOut))/fOut;
    return 1.0-value;
}

float SphereDist(vec3 position, vec3 rposition, float radius)
{
    return length(position - rposition) - radius;
}

float Star(vec3 position, vec3 rotdir, vec3 rposition, float radius)
{
    float distance;
    float t=spin_speed*iGlobalTime/radius;
    float noise = Turbulence(vec3(Spin(t*0.25+rotdir.y*.25)*rotdir.zx*50.,rotdir.y*50.+t), 0.1, 1.5, 0.03) * NoiseAmplitude*radius;
    noise = clamp((abs(noise)),0.0,1.0);
    distance = SphereDist(position, rposition, radius) - noise;
    return distance;//noise;
}

//-------------------------Vortex-------------------------
// this part from "Tornado 1" https://www.shadertoy.com/view/4tSSDK
float ridged(float f){
  return 1.0-2.0*abs(f);
}

float Vortex(vec3 q) // the isosurface shape function, the surface is at o(q)=0
{
    q.y += -100.;
    q.x += 45.;
    float h = 80.0;
  float t=spin_speed*iGlobalTime;
    //q.z=abs(q.z);
  if(q.y<0.0) return length(q)-5.75;
    vec3 spin_pos=vec3(Spin(t-sqrt(q.y))*q.xz,q.y-t*5.0);
  float zcurve=pow(q.y,1.5)*0.03;
  float v=abs(length(q.xz)-zcurve)-5.5-clamp(zcurve*0.2,0.1,1.0)*pn(spin_pos*vec3(0.1,0.1,0.1))*5.0;
  //v=v-ridged(pn(vec3(Spin(t*1.5+0.1*q.y)*q.xz,q.y-t*4.0)*0.3))*1.2;
  return max(v, q.y - h);
}

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

//-------------------------Map----------------------------
float map(vec3 p)
{
    float t=spin_speed*iGlobalTime/.25;
    float d1 = Star(p,p,vec3(-1.0,6.5,0.0),2.5);
    float d2 = Vortex(p/0.02)*0.02;
    float d3 = Disk(p,vec3(1.0,0.98,0.322)) + fpn(vec3(Spin(t*0.125+p.z*.80)*p.xy*80.,p.z*40.-t)*6.0) * 0.545;
    float d4 = Star(p,p.yzx,vec3(0.0,0.0,0.0),0.3);
    float d=smin(d3,d4,0.85);
    d=smin(d2,d,1.0);
    d=smin(d,d1,1.5);
    return d;
}

vec4 raymarch(vec3 ro,vec3 rd)
{
    float ld=0., td=0., w;
    float d=1., t=0.;
    const float h = .1;
    float tc = 0.;
    vec2 seed = rd.xy + fract(iGlobalTime);

    // rm loop
    for (int i=0; i<48; i++)
    {
      if(td>(1.-1./200.) || d<0.001*t || t>22.*8.) break;
      d = map(ro+t*rd);
      ld = (h - d) * step(d, h);
      w = (1. - td) * ld;
      tc += w*w + 1./50.;  // Different weight distribution.
      td += w + 1./50.;
      d=abs(d)*(.8+0.28*rand(seed*vec2(i)));
      d = max(d, 0.04);
      t += d*0.5;
    }
    if (tc<=0.3) return vec4(0.0,0.,.3,0.);
    vec3 col = firePalette(clamp(tc,0.,1.));
    return vec4(col, 1.0)*smoothstep(0.3,0.4,tc);
}







//-------------------------Main---------------------------
void main()
{
 vec3 ro  =  camera.xyz*camera.w*8.;
 vec3 rd  =  normalize(dir.xyz-camera.xyz);
 fragColor=raymarch(ro,rd);
}
