
layout(location = 0, index = 0) out vec4 outcol;

uniform sampler2D tex0;
uniform vec4 time_s;
uniform vec4 parms;
uniform vec4 parms2;
uniform vec4 parms3;

uniform vec4 camctr;
uniform vec4 campos;

in vec2 texcoord0;

const float PI = 3.14159265358979f;


const float NEAR_CLIP = 2.8;
const float FAR_CLIP = 30.0;

const int NUM_ITERATIONS = 128;
const float NUM_ITERATIONS_F = float(NUM_ITERATIONS);
const float TERM_DIST = 0.1;
const float STEP_MULT = 0.25;
//const vec3 GRAD_EPS = vec3( 0.025, 0, 0 );

//const vec2 ASPECT = vec2( 16.0/9.0, 1.0f );//TODO: input

//TODO: fade to this
const vec2 ASPECT = vec2( 2.35, 1.0f ); //TODO: input
const vec2 ASPECT0 = vec2( 16.0f/9.0f, 1.0f );

// ============================================================
// most primitive distance-functions honestly stolen from http://iquilezles.org/www/articles/distfunctions/distfunctions.htm

float sat( float t )
{
	return clamp( t, 0, 1 );
}
vec3 sat( vec3 t )
{
	return clamp( t, 0, 1 );
}

//remaps interval [a;b] to [0;1]
float remap( float t, float a, float b ) {
	return clamp( (t - a) / (b - a), 0, 1 );
}

// ====
vec4 select( vec4 a, vec4 b ) {
	return (a.x<b.x) ? a : b;
}

// ====
// subtract
vec4 select_max( vec4 a, vec4 b ) {
	return (a.x>b.x) ? a : b;
}

// ====
// union
vec4 select_min( vec4 a, vec4 b ) {
	return (a.x<b.x) ? a : b;
}

// ====
vec2 rot2d( vec2 p, float a )
{
	vec2 sc = vec2(sin(a),cos(a));
	vec2 ret;
	ret.x = dot( p, sc.yx*vec2(1,-1) );
	ret.y = dot( p, sc.xy );
	return ret;
}

// ====
vec3 mirror( vec3 p, vec3 pn_ctr, vec3 pn_norm )
{
	vec3 v0 = p-pn_ctr; 
	float dp0 = dot( v0, pn_norm );
	if ( dp0 > 0 ) {
		return p;
	}
	return p - 2 * dp0 * pn_norm;
}

// ====
//note: local sphere, radius r
float sdSphere( vec3 p, float r ) {
	return length( p ) - r;
}

// ====
vec3 spherical_uv( vec3 v ) {
	float rad = length( v );
	float theta = atan( v.y, v.x );
	float phi = acos( v.z / rad );
	return vec3( rad, theta, phi );
}

// ====
//note: local box, bounds b, radius r
float udRoundBox( vec3 p, vec3 b, float r )
{
	return length(max(abs(p)-b,0.0))-r;
	//vec2 uv = vec2(0.0f);
	//const float matid = 1;
	//return vec4( dist, uv, matid );
}

// ====
float sdBox( vec3 p, vec3 b ) {
	vec3  di = abs(p) - b;
	float mc = max( di.x, max(di.y, di.z) );
	float dist = min(mc, length(max(di,0.0)));
	return dist;
}

float sdPlane( vec3 p, vec3 pointonplane, vec3 norm )
{
	return dot( norm, p-pointonplane);
}

float sdTube( vec3 p, float r0, float r1 )
{
	float ld0 = length(p.xz) - r0;
	float ld1 = length(p.xz) - r1;
	float linedist = max( -ld0, ld1 );

	float dp0 = sdPlane( p, vec3(0,-1,0), vec3( 0, -1, 0 ) );
	float dp1 = sdPlane( p, vec3(0,1,0), vec3( 0, 1, 0 ) );
	return max( linedist, max( dp0, dp1 ) );
}

// ====
vec4 plane_uv( vec3 p, vec3 pointonplane, vec3 norm ) {
	float dist = sdPlane( p, pointonplane, norm );
	vec2 uv = vec2( 0, 0 ); //TODO
	float matid = 3;
	return vec4( dist, uv, matid );
}

// ====
// \/ [0;1]
float saw( float t ) {
	return abs(mod(abs(t), 2.0f)-1.0f);
}

float svinus( float t ) {
	return 0.5+0.5*sin((t+0.5)*PI);
}

//====
// normalized, y=[-1;1], r
float sdCylinder( vec3 p, float r )
{
	float l = length( p.xz ) - r;
	float dp0 = sdPlane( p, vec3(0,-1,0), vec3( 0, -1, 0 ) ); //TODO: -p.z+1
	float dp1 = sdPlane( p, vec3(0, 1,0), vec3( 0,  1, 0 ) ); //TODO:  p.z-1
	return max( l, max( dp0, dp1 ) );
}

//note: from http://www.iquilezles.org/www/articles/distfunctions/distfunctions.htm
float sdHexPrism( vec3 p, vec2 h )
{
    vec3 q = abs(p);
    return max(q.z-h.y,max(q.x+q.y*0.57735,q.y*1.1547)-h.x);
}

//note: [-1;1]
//note: honestly stolen from iq: https://www.shadertoy.com/view/Xsl3Dl
vec3 shash3( vec3 p )
{
	p = vec3( dot(p,vec3(127.1,311.7, 74.7)),
			  dot(p,vec3(269.5,183.3,246.1)),
			  dot(p,vec3(113.5,271.9,124.6)));

	return -1.0 + 2.0*fract(sin(p)*43758.5453123);
}
//note: [-1;1]
//float snoise( in vec3 p )
//{
//    vec3 i = floor( p );
//    vec3 f = fract( p );
//
//	vec3 u = f*f*(3.0-2.0*f);
//
//    return mix( mix( mix( dot( shash3( i + vec3(0.0,0.0,0.0) ), f - vec3(0.0,0.0,0.0) ), 
//                          dot( shash3( i + vec3(1.0,0.0,0.0) ), f - vec3(1.0,0.0,0.0) ), u.x),
//                     mix( dot( shash3( i + vec3(0.0,1.0,0.0) ), f - vec3(0.0,1.0,0.0) ), 
//                          dot( shash3( i + vec3(1.0,1.0,0.0) ), f - vec3(1.0,1.0,0.0) ), u.x), u.y),
//                mix( mix( dot( shash3( i + vec3(0.0,0.0,1.0) ), f - vec3(0.0,0.0,1.0) ), 
//                          dot( shash3( i + vec3(1.0,0.0,1.0) ), f - vec3(1.0,0.0,1.0) ), u.x),
//                     mix( dot( shash3( i + vec3(0.0,1.0,1.0) ), f - vec3(0.0,1.0,1.0) ), 
//                          dot( shash3( i + vec3(1.0,1.0,1.0) ), f - vec3(1.0,1.0,1.0) ), u.x), u.y), u.z );
//}
//float nnoise( in vec3 p )
//{
//	return 0.5 + 0.5 * snoise( p );
//}

//const mat3 m = mat3( 0.00,  0.80,  0.60,
//                    -0.80,  0.36, -0.48,
//                    -0.60, -0.48,  0.64 );
//
//float nfbm( vec3 q )
//{
//	float f = 0;
//	f  = 0.5000*nnoise( q ); q = m*q*2.01;
//	f += 0.2500*nnoise( q ); q = m*q*2.02;
//	f += 0.1250*nnoise( q ); q = m*q*2.03;
//	f += 0.0625*nnoise( q ); q = m*q*2.01;
//	
//	return f;
//}
//float sfbm( vec3 q )
//{
//	float f = 0;
//	f  = 0.5000*snoise( q ); q = m*q*2.01;
//	f += 0.2500*snoise( q ); q = m*q*2.02;
//	f += 0.1250*snoise( q ); q = m*q*2.03;
//	f += 0.0625*snoise( q ); q = m*q*2.01;
//	
//	return f;
//}

float trunc( float x, float num_levels )
{
	return floor(x*num_levels) / num_levels;
}
vec3 trunc( vec3 x, vec3 num_levels )
{
	return floor(x*num_levels) / num_levels;
}

float sdHalfCircle( vec3 p )
{
	vec3 ofs = shash3( trunc( p, vec3(13)));
	p += 0.02  * shash3( trunc( p, vec3(32) + ofs ) );
	p += 0.01  * shash3( trunc( p, vec3(16) ) );

	const float spherer = 1.0;
	float ds = sdSphere( p, spherer );
	ds = max( ds, -sdSphere( p, 0.9 * spherer ));
	ds = max( ds,  p.x-0.1);  //cut sides
	ds = max( ds, -p.x-0.1); // -
	
	return ds + parms2.z;
}

vec4 cylinder_radialrepeat( vec3 p, float r )
{
	float d = FAR_CLIP;
	
	//lpt logo
	if ( parms3.x > -0.5 )
	{
		const vec3 logoctr = vec3( -0.05, 0.98, 0 );
		const float scl = 0.25;
		const vec3 box0_pos = logoctr + scl * vec3(  0.5, 0, -0.4 );
		const vec3 box1_pos = logoctr + scl * vec3(  0.5, 0, -0.1 );
		const vec3 box2_pos = logoctr + scl * vec3(  0.5, 0,  0.2 );
		const vec3 box4_pos = logoctr + scl * vec3(  0.2, 0,  0.2 );
		const vec3 box3_pos = logoctr + scl * vec3(  0.2, 0, -0.1 );
		const vec3 box5_pos = logoctr + scl * vec3( -0.1, 0,  0.2 );
		vec3 box_siz = scl*vec3( parms3.x, parms3.y, parms3.x );
    
		float dlpt = sdBox( p-box0_pos, box_siz );
		dlpt = min( dlpt, sdBox( p-box1_pos, box_siz ) );
		dlpt = min( dlpt, sdBox( p-box2_pos, box_siz ) );
		dlpt = min( dlpt, sdBox( p-box3_pos, box_siz ) );
		dlpt = min( dlpt, sdBox( p-box4_pos, box_siz ) );
		dlpt = min( dlpt, sdBox( p-box5_pos, box_siz ) );
    
		d = min( d, dlpt );
		
		//d = max(-d, p.y-1); //cut top
	}

	p *= parms2.w;
	{
		//vec2 dp2d = rot2d( p.xz, time_s.x );
		vec3 dp = p; //vec3( dp2d.x, p.y, dp2d.y );

		dp.zx = rot2d( dp.xz, 0.3*time_s.x );
		float ds = sdHalfCircle( dp );

		dp.xy = rot2d( dp.xy, -0.7 * time_s.x );
		ds = min ( ds, sdHalfCircle( dp.xyz*1.3 ) );

		dp.xz = rot2d( dp.xz, 1.1*time_s.x );
		ds = min ( ds, sdHalfCircle( dp.xyz*1.8 ) );

		dp.xz = rot2d( dp.xz, -1.3*time_s.x );
		ds = min ( ds, sdHalfCircle( dp.xyz*2.8 ) );

		d = min( d, ds );
	}

	return vec4( d, vec2(0,0), 0 );
}

// ============================================================

// ====
//note: ray_dir should be normalized
void isect_RayPlane( vec3 ray_org, vec3 ray_dir, vec4 plane, out float t )
{
	vec3 n = plane.xyz;
	float d = plane.w;
	t = -( dot(ray_org, n) + d ) / dot( ray_dir, n );
}

vec3 isect_RayPlane( vec3 ray_org, vec3 ray_dir, vec4 plane )
{
	float t;
	isect_RayPlane( ray_org, ray_dir, plane, t );
	return ray_org + t * ray_dir;
}

void init1( out vec3 cam_eye, out vec3 cam_lookat, out vec3 cam_up ) {
	//cam_eye    = vec3( 5*cos(0.2*time.x), 5, 5*sin(0.2*time.x) );
	cam_eye    = vec3( 10, 10, 10 );
	cam_lookat = vec3( 0, 0, 0 );
	cam_up     = normalize( vec3( 0.2, 1, 0 ) );
}

vec4 scene1( vec3 p )
{
	vec4 d = vec4( FAR_CLIP, 0.0f, 0.0f, -1.0f ); //note: background

	vec3 cyl_pos = vec3( 0, -10, 0 );
	float cyl_scl = 10.0;
	float cyl_r = 0.1;
	vec4 crr = cylinder_radialrepeat( (p-cyl_pos)/cyl_scl, cyl_r*0.5 );
	crr.x *= cyl_scl;
	d = select( d, crr );

	return d;
}

// ===
void init( out vec3 cam_eye, out vec3 cam_lookat, out vec3 cam_up )
{
	init1( cam_eye, cam_lookat, cam_up );
}

// ====
vec4 scene( vec3 p )
{
	return scene1( p );
}

// ============================================================

// ====
//vec3 grad( vec3 p, float d )
//{
//	return normalize( vec3( scene( p + GRAD_EPS.xyz ).x - d,
//							scene( p + GRAD_EPS.zxy ).x - d,
//							scene( p + GRAD_EPS.yzx ).x - d ) );
//}

//note: more expensive version
//vec3 grad( vec3 p )
//{
//	return normalize( vec3( scene( p + GRAD_EPS.xyz ).x - scene( p - GRAD_EPS.xyz ).x,
//							scene( p + GRAD_EPS.zxy ).x - scene( p - GRAD_EPS.zxy ).x,
//							scene( p + GRAD_EPS.yzx ).x - scene( p - GRAD_EPS.yzx ).x ) );
//}


// ====
float fakeao( vec3 p ) {
	const vec3 sphere_ctr = vec3( 0, 3, 0 );
	const float sphere_rad = 4;
	float l = length( sphere_ctr - p ) / 5;
	return smoothstep( 0, 1, pow( l, 2) );
}

// ====
float ao( vec3 p, vec3 n, float max_dist ) {
	const int num_samples = 5;
	const float num_samples_f = float(num_samples);
	float stepsiz = max_dist / num_samples_f;
	const float falloff = 3.0f;

	float d = 1.0f;
	for ( int i=1; i<num_samples; i++ )
	{
		float t = float(i) / num_samples_f;
		float s = t * max_dist;
		float val = scene( p + s * n ).x - s;
		d *= 1.0f - val / pow( falloff, float(i) );
	}

	return clamp( d, 0, 1 );
}

// ====
float topReflect( vec3 norm, vec3 dir )
{
	float v = step( norm.x, norm.y ) * step( norm.z, norm.y );
	const float falloff = 4;
	const float intensity = 3;
	float fresnel = pow( 1 - dot( norm, -dir ), falloff ); //TODO: schlick
	return (0.15+v) * fresnel * intensity;
}

//vec3 get_material_col( vec4 d ) {
//
//	vec4 matid = d.wwww;
//
//	vec3 mat0_col = texture( tex0, d.yz*vec2(2,1)-vec2(0.5f, 0) ).rgb;
//	vec3 mat1_col = mix( vec3( 0.9f, 0.75f, 0.1f ), vec3(0.0f), clamp(parms.z,0,1) );
//	//vec3 mat2_col = vec3(0.8f); //1.25f*vec3( 0.5f, 0.6f, 0.7f );
//	vec3 mat2_col = vec3(0.05f);
//
//	vec4 limits = vec4(-0.5f, 0.5f, 1.5f, 2.5f);
//	vec4 matid_less = step( limits, matid );
//	vec4 matid_great = step( matid, limits+1);
//	vec4 matid_mask = matid_less * matid_great;
//
//	//return mat2_col;
//
//	return matid_mask.x * mat0_col
//		 + matid_mask.y * mat1_col
//		 + matid_mask.z * mat2_col;
//}

// ====
//note: mostly used thingy...
//[0;1]
float nrand( vec2 n ) {
	return fract(sin(dot(n.xy, vec2(12.9898, 78.233)))* 43758.5453);
}

// ============================================================

vec4 raymarch( inout vec3 p, inout vec3 dir, out int out_steps )
{
	int i=0;
	vec4 d;
	float rdt = 0;
	for ( i=0; i<NUM_ITERATIONS; i++ )
	{
		d = scene( p );

		if ( (abs(d.x) < TERM_DIST ) || (rdt > FAR_CLIP) ) {
			break;
		}
		else
		{
			float dt = STEP_MULT * d.x; //note: constant-multiply to compensate for distorted space, actual dist < dist - could use gradient-approximation instead? (see iq)
			p += dir * dt;
			rdt += dt;
		}
	}

	out_steps = i;
	return d;

}

//float domainrepeat( float x, float period )
//{
//	return mod( x + 0.5 * period, period ) - 0.5 * period;
//}

// ====
//TODO: animated twist-modifier
//TODO: blob intensities
//void densityField( vec3 p, out float out_dens, out vec3 out_col )
//{
//	const vec3 itsspherectr = vec3(0, 0, 0);
//	float itssphererad = parms3.w;
//	
//	//float sphereits = (sdSphere( p-itsspherectr, itssphererad ) < 0) ? 1.0 : 0.0;
//	float sphereits = sat(-sdSphere( p-itsspherectr, itssphererad ));
//	//float sphereits = sin( sdSphere( p-itsspherectr, itssphererad ) );
//
//	vec2 rotctr = vec2(1,2);
//	vec3 mp1 = p;
//	mp1.xz = rotctr+rot2d( mp1.xz-rotctr, parms3.z * length(rotctr-p.xz) );
//	vec3 modp = mp1 / sqrt( parms3.x );
//	modp.x = domainrepeat ( modp.x, 0.5 );
//	modp.y = domainrepeat ( modp.y, 0.5 );
//	modp.z = domainrepeat ( modp.z, 0.5 );
//	
//	const float sphere_r = 0.125;
//	float spheredist = sdSphere( modp, sphere_r );
//	float planedist = sdPlane( p, vec3(0,-1.0,0), vec3(0,1,0) );
//
//	float vits = sat( 2*sin( p.y ) );
//	
//	float d = max( spheredist, planedist );
//	float spheredensity = 0;
//	//spheredensity += spheredist < 0 ? -0.5*spheredist : 0;
//	spheredensity += (d<0) ? vits * 0.4 * sphereits : 0.0;
//	
//	out_dens = spheredensity;
//	out_col = vec3(0.7,0.5, 1) * spheredensity;
//}


// ====
//vec4 volumemarch( vec3 in_porg, vec3 in_p )
//{
//	const int NUM_VOL_ITER = 128;
//	const float NUM_VOL_ITER_F = float(NUM_VOL_ITER);
//	
//	vec3 p = in_porg;
//	vec3 stepvec = (in_p - in_porg) / NUM_VOL_ITER_F;
//	float stepveclen = length( stepvec );
//
//	float accum_density = 0;
//	vec3 accum_color = vec3(0);
//	for ( int i=0; i<NUM_VOL_ITER; ++i )
//	{
//		float point_dens;
//		vec3 point_col;
//		
//		densityField( p, point_dens, point_col );
//		
//		//TODO: density should obscure bg?
//		accum_color += stepveclen * point_col * (1-accum_density );
//		accum_density += stepveclen * point_dens; //note: assume uniform density within a step
//		
//		p += stepvec;
//	}
//
//	return vec4( vec3(accum_density), 1 );
//	//return vec4( vec3(accum_color), 1 );
//}


// ====
void main()
{
	vec3 cam_eye;
	vec3 cam_lookat;
	vec3 cam_up;

	init( cam_eye, cam_lookat, cam_up );

	cam_eye = campos.xyz;
	cam_lookat = camctr.xyz;


	vec3 negz = normalize( cam_lookat - cam_eye );
	vec3 u = cross( negz, cam_up );
	vec3 v = cross( u, negz );

	u *= texcoord0.x;
	v *= texcoord0.y;

	vec2 aspect = mix( ASPECT, ASPECT0, parms2.x );
	
	vec2 dim = 0.25f * aspect; //controls fov
	u *= dim.xxx;
	v *= dim.yyy;

	float dist = 0.25f; //...also controls fov :p
	vec3 dir = dist * negz + u + v;
	dir = normalize( dir );

	vec3 p_org = cam_eye + NEAR_CLIP * dir;
	vec3 p = p_org;

	//note: scene-specific optimization
	//note: project to top plane to begin with...
	//vec3 p_plane = isect_RayPlane( p, dir, vec4(0,1,0,-2) ); //note: -2 should be 0... offset is for iteration-glow :D
	//p = p_plane;
	
	float rnd = nrand( texcoord0.xy + 0.01f * fract( time_s.x ) );

	
	//note: jitter startpos
	p -= 4*TERM_DIST * dir * rnd;

	int i;
	vec4 d = raymarch( p, dir, i );

	
	bool valid = d.x < 10 * TERM_DIST ;


	vec3 matcol = vec3( 0.5, 0.6, 0.7 );

	outcol = vec4( 0,0,0, valid ? 1 : 0 );


	
	float globelo = valid ? (-parms.x * (p.y+1.0)) : 0;
	globelo *= globelo;

	outcol.rgb *= outcol.a;
	outcol.rgb = max( vec3(0), outcol.rgb );
	
	//note: iteration-glow
	float it_f = float(i) + 0.5f - rnd;
	float iterations = it_f / NUM_ITERATIONS_F;
	float glowits = parms2.y * pow(iterations,2.0);
	const vec3 gc0 = vec3(104,79,255)/255.0f;
	vec3 glowcol = mix( 0.5*gc0, 2*gc0, parms.x );
	//outcol.rgb += max(vec3(0),glowits * glowcol);
	outcol.rgb = 3*glowits * glowcol;
	outcol.rgb += 0.3*globelo*gc0;

	//note: vignette
	float vign = 1.5-length(texcoord0.xy*vec2(2.35,1));
	vign = pow(vign, 2);
	outcol.rgb *= vec3( sat(vign) );

	//float m = mod( gl_FragCoord.y, 2.0f );
	//outcol.rgb *= 0.5 + 0.5*m;
	

	//note: angel overlay
	if ( parms.y > 0 )
	{
		float s0 = 1 + 0.25 * sin( 89 * texcoord0.x + 1.1*time_s.y - 1.89*time_s.x );
		float s1 = 1 + 0.25 * sin( 97 * texcoord0.x - 1.3*time_s.y + 1.13*time_s.x );
		float s2 = 0.5 + 0.5 * sin ( sin( 17 * texcoord0.x - 2 * time_s.y + 1.07*time_s.x ) * cos ( 3 * texcoord0.x + time_s.y - 1.71*time_s.x ) );
		float s0_h = pow ( remap( texcoord0.y, -0.2, 0.5), 3 );
		float s1_h = pow ( remap( texcoord0.y, -0.5, 0.5), 1.75 );
		float s2_h = pow ( remap( texcoord0.y, -0.45, 0.5), 2 );
		vec3 aits = vec3( 0.3 * (s0_h*s0 + s1_h*s1) + s2_h*s2 );
		//aits = sat( aits );
		float lits = smoothstep( 0, 1, 1-abs(texcoord0.x-parms.z) );
		aits *= lits;

		//screen
		outcol.rgb = mix( outcol.rgb, 1-(1-parms.y*gc0*vec3(aits))*(1-outcol.rgb), 0.25 );

		//color dodge
		outcol.rgb = mix ( outcol.rgb, outcol.rgb / (1.0f-parms.y*sat(aits)), 0.35 );
	}

	//add
	//outcol.rgb += vec3( parms.y * aits );

	//dbg
	//outcol.rgb = vec3( parms.y * aits );
}
