 // License Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
// Created by AVR 2016
// Galaxy morphology based on http://iopscience.iop.org/0004-637X/783/2/138/pdf/0004-637X_783_2_138.pdf

//----------------------------------------
const float hbulge   = 1.5; // [1;2]
const float  hgal  = 0.1; // galaxy thickness
const vec3  sbulge = vec3(.5,hgal*hbulge,.5); //  bulge size
//----------------------------------
const float arms = 2.0; //[1;4]
const float Wb   = 0.3;  // twist parameter 1
const float Wn   = 12.0; // twist parameter 2
const float iAw  = 4.;   // arms thickness
//------- x=interstellar gas,y=clouds,z=dust ------
const vec3 tilt  = vec3(1.,2.5,2.);
const vec3 scale = vec3(1.,2. ,4.)*16.;
const vec3 D0    = vec3(2.,8. ,1.2)*.5;    // Density between arms
const vec3 Ds    = vec3(2.,6. ,3.2);       // Density scale
//----------------------------------
const float ihI=.6;   // ionized hydrogen intensity
const float CloudsDensity=2.;
const float Ibulge = 1.5; // Bulge Intensity

const mat3 mx = mat3(0.,0.8,0.6,-0.8,0.36,-0.48,-0.6,-0.48,0.64);
const mat2 m2 = mat2(.8,.6,-.6,.8);

mat2 rot(float a)
{
  float sa = sin(a);
  float ca = cos(a);
  return mat2(ca,sa,-sa,ca);
}

float noise3( 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 = texture2D( iChannel0, (uv+0.5)/256.0).yx;
  float v=mix( rg.x, rg.y, f.z );
    return v;
}

float noise(in vec3 p)
{
  p *= 2.0;
  float res = 0.0;
  float f = 1.0;
  for( int i=0; i<3; i++ )
  {
    p.xy = m2*p.xy;
    p = p.zxy*f+.6;
    f *= 1.15;
    res += sin(p.y+1.3*sin(1.2*p.x)+1.7*sin(1.7*p.z));
  }
  return res / 3.0;
}

vec3 fbm( vec3 p) // x=interstellar gas,y=clouds,z=dust,w=density
{
  vec3 pg = p;
  pg.y*=1./hgal;
  float t = atan(pg.z,pg.x);
  float r = length(pg.xz);
  float angle = (atan(exp(1./r)/Wb)*2.0*Wn-t)*arms;
  vec3  f=vec3(0.0);
  float dh=smoothstep(0.2,.4,r)*(1.-1./hbulge);
  float fade=smoothstep(0.0,.3,1.-r)*max(0.0,1.-abs(pg.y)-dh)*smoothstep(0.0,0.1,r);
  float D = pow(abs(cos(angle)),iAw);
  vec3  Dg=fade*fade*hbulge*(Ds*D+D0);
  for (float k=1.;k<7.;k++)
  {
     vec3 k1 = k*scale;
     f.x +=      noise3( k1.x*p)/k;
     f.y +=      noise3( k1.y*p)/k;
     f.z += 1.0/(noise3( k1.z*p)/k);
     p = mx*p;
  }
  f.x=.2/f.x;
  f.xyz=pow(f.xyz,tilt);
  f.yz=exp(vec2(-f.y,-1.e-4*abs(f.z)))*Dg.yz;
  f.y+=f.z;
  f.x*=Dg.x;
  f.z=Dg.x;
  return f;
}

float stars( vec3 p )
{
  p = p;
  float f = 1.0;
  float r = 0.0;
  for(int i = 1;i<5;i++)
  {
    r += noise( p*(20.+3.*f) )/f;
    p.xz *= m2;
    f += 1.0;
  }
  return pow(abs(r),8.);
}

// http://adsabs.harvard.edu/abs/1987A&A...175....1M
float bulge(vec3 p)
{
  p/=sbulge;
  float r=length(p);
  float I=pow(r,-0.855)*exp(-pow(r,.25));
  return smoothstep(0.0,.5,.5-r)*(I*Ibulge);
}

bool RSI(in vec3 ro,in vec3 rd,in vec3 scale,out vec2 res)
{
  ro*=scale;
  rd*=scale;
  float k=1.0/length(rd);
  rd *=k;
  float b = dot(ro,rd);
  float c = dot(ro,ro) - 1.0;
  float d = b*b - c;

  if(d <= 0.0) return false;
  d = sqrt(d);
  res.x=(-b - d)*k;
  res.y=(-b + d)*k;
  return true;
}


vec4 castray( in vec3 ro, vec3 rd)
{
  vec2 tt,tt1;
  bool fg=RSI(ro,rd,vec3(1.,1./hgal,1.),tt);
  bool fb=RSI(ro,rd,2./sbulge,tt1);
  if (fg||fb)
  {
   const float dt = 0.01;
   vec4  alight  = vec4(0.);
   vec2  ahydro  = vec2(0.);
   vec2  ashadow = vec2(1.0,0.0);
   float abulge  = 0.;
   float astar   = 0.;
   float t,t1;
   if (fg)
       if (fb) { t = min(tt.x,tt1.x); t1= max(tt.y,tt1.y); }
       else    { t = tt.x; t1=tt.y; }
   else        { t = tt1.x; t1=tt1.y; }
   t=max(t,0.0);
   for(float i=0.;i<3.;i+=dt  )
   {
     vec3 pos = ro+t*rd;
     if (t>=tt.x&&t<=tt.y)
     {
      vec3  f=fbm(pos); // x=gas light,y=clouds+dust, z=density
      ashadow.x  = max(0.0,ashadow.x-f.y*dt*CloudsDensity);
      float   c1 = f.x*ashadow.x,c2=c1*c1;
      vec2    hd = vec2(smoothstep(ihI-.3,ihI,f.x),
                        smoothstep(ihI,ihI+.1,f.x))*ashadow.x;
      astar     += f.z*stars(pos*7.)*.02*ashadow.x;
      alight    += vec4(c2*c1,c2,c1,c2*c2);
      hd*=hd;
      ahydro    += hd*hd;
      ashadow.y += f.y*c1; // dust+clouds reflected light
     }
     if (t>=tt1.x&&t<=tt1.y) abulge += bulge(pos)*ashadow.x;
     t+=dt;
     if (t>t1) break;
   }
   vec3 col=astar*vec3( 8., 8.,4.);    // star light
   col+=alight.rgb*vec3( 8.0, 4.0,1.8);  // interstellar gas light
   col+=abulge*vec3(1.1, .95,.5); // bulge light
   col+=ashadow.y*.16;           // dust+clouds reflected light
   col  +=ahydro.x*vec3(5.1,2.78,1.88);
   col.r+=ahydro.y*10.; // glow of ionized hydrogen
   return vec4(col*dt*2.0,ashadow.x);
 }
 return vec4(0.,0.,0.,1.);
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
  vec2 pos= (fragCoord.xy - 0.5*iResolution.xy )/ iResolution.y;
  vec3 ro = vec3(0.,0.,-2.165);
  vec3 rd = normalize(vec3(pos,2.));
  mat2 mx=rot(-iGlobalTime*.1);
  mat2 my=rot(-.5);
  ro.yz*=my;
  ro.xz*=mx;
  rd.yz*=my;
  rd.xz*=mx;
  vec4 col = castray(ro,rd);
  col.rgb +=(.003*col.a*stars(rd*8.));
  fragColor = vec4( pow(col.rgb,vec3(.4545)), 1.0/*-col.a*/ );
}

