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

uniform sampler2D tex1;
#define iChannel1 tex1

uniform vec4 param[5];

#define  sbulge  param[0].xyz // bulge size
#define  hbulge  param[0].w   // bulge height

#define  hgal    param[1].x   // galaxy thickness
#define  Wb      param[1].y   // twist parameter 1
#define  Wn      param[1].z   // twist parameter 2
#define  iAw     param[1].w   // arms thickness

#define  D0      param[2].xyz // Density between arms
#define  ihI     param[2].w   // ionized hydrogen intensity

#define  Ds      param[3].xyz // Density scale
#define  Ibulge  param[3].w   // Bulge Intensity

#define  bulgeFade param[4].x   // Bulge fade
#define  edgeFade  param[4].y   // Edge fade


const vec3 tilt  = vec3(1.,2.5,2.);
const vec3 scale = vec3(8.,16. ,1.6);
const float Istar= .5;
const float CloudsDensity=1.;

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

float noise_3( 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).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,float D) // x=interstellar gas,y=clouds,z=dust,w=density
{
  vec3 pg = p;
  pg.y/=hgal;
 float r = length(pg.xz);
  vec3  f=vec3(0.0);
  float dh=bulgeFade*(smoothstep(sbulge.x*0.4,sbulge.x*.8,r));
  float fade=smoothstep(0.0,edgeFade,1.-r)*
             max(0.0,1.-dh-abs(pg.y));
  fade *= smoothstep(0.0,0.5,length(p.xz/sbulge.xz));

  vec3  Dg=fade*fade*hbulge*(Ds*D+D0*.5);
  for (float k=1.;k<7.;k++)
  {
     vec3 k1 = k*scale;
     float n1 = noise_3( k1.x*p)/k;
     f.x += n1;
     f.y += noise_3( k1.y*p)/k;
     f.z += 1.0/n1;
     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.);
}

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  ashadow = vec2(1.0,0.0);
   float abulge  = 0.;
   float astar   = 0.;
   float t,t1;

   ivec2 ts=textureSize(tex1,0);
   float ratio=float(ts.x)/float(ts.y);
   vec2 tc   = ro.xz-(ro.y/rd.y)*rd.xz;
   vec4 tcol;
   tcol     = texture( iChannel1, vec2(tc.x,tc.y*ratio)*.5+.5);
   tcol.a   = 0.0;

   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);
   vec3 pos = ro+t*rd;
   float y0 = pos.y;
   for(;;)
   {
     if (t>=tt.x&&t<=tt.y)
     {
      float b = texture( iChannel1, vec2(pos.x,pos.z*ratio)*.5+.5).a;
      vec3  f = fbm(pos,b); // x=gas light,y=clouds+dust, z=density
      float   sh = max(0.0,ashadow.x-f.y*dt*CloudsDensity);
      if (pos.y*y0<=0&&pos.y!=y0)
      {
       tcol.a=mix(ashadow.x,sh,(0-y0)/(pos.y-y0))*16.0;
      }
      ashadow.x  = sh;
      float   c1 = f.x*ashadow.x,c2=c1*c1;
      astar     += f.x*f.z*stars(pos*7.)*ashadow.x;
      alight    += vec4(c2*c1,c2,c1,c2*c2)*vec4( 8.0, 4.0,1.8,2.0);
      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;
     y0   = pos.y;
     pos += dt*rd;
   }
   vec4 col=alight;  // interstellar gas light
   col.rgb+=tcol.rgb*tcol.a;
   col+=(astar*Istar)*vec4( .8, .8,.4,0.);    // star light
   col+=abulge*vec4(1.1, .95,.5,.01); // bulge light
   col.rgb+=ashadow.y*.16;           // dust+clouds reflected light
   return vec4(col.rgb*(dt*2.0),col.a);
 }
 return vec4(0.);
}

void main( )
{
  vec3 ro  =  camera.xyz*camera.w;
  vec3 rd  =  normalize(dir.xyz-camera.xyz);
  vec4 col =  castray(ro,rd);
  fragColor = col*dir.w;
}
