#version 330 core
layout (location = 0) out vec4 prtPosition;
layout (location = 1) out vec4 prtVelocity;

in vec2 TexCoords;

uniform sampler2D prtPosTex;
uniform sampler2D prtVelTex;

uniform float delta_time;

uniform int prtAddSt;
uniform int prtAddEnd_e;
uniform int prtAddEnd_s;
uniform int prtMax;
uniform int prtTexWH;

uniform vec3 prtPosRng;
uniform vec3 prtPosOffset;
uniform vec3 prtVelRng;
uniform vec3 prtVelOffset;
uniform vec3 prtPosRngSt;
uniform vec3 prtPosOffsetSt;
uniform vec3 prtVelRngSt;
uniform vec3 prtVelOffsetSt;
uniform vec3 invPrtPosRng;
uniform vec3 invPrtVelRng;
uniform vec3 prtAcc;

//---------------------------------------------------------------------------------
// random function taken from Spatial answer to "Random / noise functions for GLSL"
// https://stackoverflow.com/a/17479300/10981777
//
// A single iteration of Bob Jenkins' One-At-A-Time hashing algorithm.
uint hash( uint x ) {
    x += ( x << 10u );
    x ^= ( x >>  6u );
    x += ( x <<  3u );
    x ^= ( x >> 11u );
    x += ( x << 15u );
    return x;
}

// Compound versions of the hashing algorithm I whipped together.
uint hash( uvec2 v ) { return hash( v.x ^ hash(v.y)                         ); }
uint hash( uvec3 v ) { return hash( v.x ^ hash(v.y) ^ hash(v.z)             ); }
uint hash( uvec4 v ) { return hash( v.x ^ hash(v.y) ^ hash(v.z) ^ hash(v.w) ); }

// Construct a float with half-open range [0:1] using low 23 bits.
// All zeroes yields 0.0, all ones yields the next smallest representable value below 1.0.
float floatConstruct( uint m ) {
    const uint ieeeMantissa = 0x007FFFFFu; // binary32 mantissa bitmask
    const uint ieeeOne      = 0x3F800000u; // 1.0 in IEEE binary32

    m &= ieeeMantissa;                     // Keep only mantissa bits (fractional part)
    m |= ieeeOne;                          // Add fractional part to 1.0

    float  f = uintBitsToFloat( m );       // Range [1:2]
    return f - 1.0;                        // Range [0:1]
}

// Pseudo-random value in half-open range [0:1].
float random( float x ) { return floatConstruct(hash(floatBitsToUint(x))); }
float random( vec2  v ) { return floatConstruct(hash(floatBitsToUint(v))); }
float random( vec3  v ) { return floatConstruct(hash(floatBitsToUint(v))); }
float random( vec4  v ) { return floatConstruct(hash(floatBitsToUint(v))); }

// From: Shadertoy user scholarius in 2019-05-06
// https://www.shadertoy.com/view/wts3RX
float cbrt( float x ){
	float y = sign(x) * uintBitsToFloat( floatBitsToUint( abs(x) ) / 3u + 0x2a514067u );
	y = ( 2. * y + x / ( y * y ) ) * .333333333;
    float y3 = y * y * y;
    y *= ( y3 + 2. * x ) / ( 2. * y3 + x );
    return y;
}

vec3 box_to_sphere( float r1, float r2, float r3 ){
    float phi = r1 * 6.283185307;
    float theta = acos( r2 * 2.0 - 1.0 );
    float r = cbrt( r3 );
    float r_sin_th = r * sin( theta );
    return vec3( r_sin_th * cos( phi ), r_sin_th * sin( phi ), r * cos( theta ) );
}

void main(){
    ivec2 TexCoord_px = ivec2(TexCoords * prtTexWH);
    int prtNmb = TexCoord_px.x + TexCoord_px.y * prtTexWH;
    if( prtNmb < prtMax ){
        vec3 prtPos;
        vec3 prtVel;
        vec3 prtPos_wrld;
        vec3 prtVel_wrld;
        if( ( prtNmb >= prtAddSt && prtNmb < prtAddEnd_e ) || ( prtNmb <= prtAddEnd_s ) ){
            vec3 rand_inputs = vec3( gl_FragCoord.xy, delta_time );
            prtPos.x = random( rand_inputs ); // [0:1)
            rand_inputs = vec3( gl_FragCoord.x, prtPos.x, delta_time );
            prtPos.y = random( rand_inputs ); // [0:1)
            rand_inputs = vec3( gl_FragCoord.y, prtPos.y, delta_time );
            prtPos.z = random( rand_inputs ); // [0:1)
            rand_inputs = vec3( gl_FragCoord.x, prtPos.z, delta_time );
            prtVel.x = random( rand_inputs ); // [0:1)
            rand_inputs = vec3( gl_FragCoord.x, prtVel.x, delta_time );
            prtVel.y = random( rand_inputs ); // [0:1)
            rand_inputs = vec3( gl_FragCoord.y, prtVel.y, delta_time );
            prtVel.z = random( rand_inputs ); // [0:1)
            // turn the box dists into spheres
            prtPos = box_to_sphere( prtPos.x, prtPos.y, prtPos.z );
            prtVel = box_to_sphere( prtVel.x, prtVel.y, prtVel.z );

            prtPos_wrld = prtPosRngSt * prtPos - prtPosOffsetSt;
            prtVel_wrld = prtVelRngSt * prtVel - prtVelOffsetSt;
        } else {
            prtPos      = texture( prtPosTex, TexCoords ).xyz;
            prtVel      = texture( prtVelTex, TexCoords ).xyz;
            prtPos_wrld = prtPosRng * prtPos - prtPosOffset;
            prtVel_wrld = prtVelRng * prtVel - prtVelOffset;
        }
        vec3 v0      = prtVel_wrld;
        prtVel_wrld += prtAcc * delta_time;
        prtPos_wrld +=  0.5f * ( v0 + prtVel_wrld ) * delta_time;
        prtPosition  = vec4( ( prtPos_wrld + prtPosOffset ) * invPrtPosRng, 1.0 );
        prtVelocity  = vec4( ( prtVel_wrld + prtVelOffset ) * invPrtVelRng, 1.0 );
    } else {
        discard;
    }
}