// OSL version Mads Drøschler
// WebGL version by Witek https://www.shadertoy.com/view/ltKyRR)
// License Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License

#define vec3 vector

float fract(float x) { return x - floor(x); }
//------------------------------------------
struct vec2{ float x, y; };
vec2 __operator__add__(vec2 a, vec2 b) { return vec2(a.x+b.x, a.y+b.y); }
vec2 __operator__div__(vec2 a, float b) { return vec2(a.x/b, a.y/b); }
// ------------------------------------------
struct vec4{ float x, y, z, w; };
vec4 __operator__mul__(vec4 a, vec4 b) { return vec4(a.x * b.x, a.y * b.y, a.z * b.z, a.w * b.w); }
vec4 __operator__sub__(float a, vec4 b) { return vec4(a - b.x, a - b.y, a - b.z, a - b.w); }
vec4 max(vec4 x, float v) { return vec4( max(x.x, v), max(x.y, v), max(x.z, v), max(x.w, v) ); }
float dot(vec4 a, vec4 b) { return a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w; }
// ------------------------------------------
vec3 random3(vec3 c)
{
    float j = 4096.0*sin(dot(c, vec3(17.0, 59.4, 15.0)));
    vec3 r;
    r[2] = fract(512.0 * j); 
    j *= 0.125;
    r[0] = fract(512.0 * j); 
    j *= 0.125;
    r[1] = fract(512.0 * j);
   
    return r - 0.5;
}

float simplex3d(vec3 p)
{
     float F3 = 0.3333333;
	  float G3 = 0.1666667;
     vec3 s = floor(p + dot(p, vec3(F3)));
     vec3 x = p - s + dot(s, vec3(G3));
     vec3 e = step(vec3(0.0), x - vec3(x[1],x[2],x[0]));
     vec3 i1 = e*(1.0 - vec3(e[2],e[0],e[1]));
     vec3 i2 = 1.0 - vec3(e[2],e[0],e[1])*(1.0 - e);
     vec3 x1 = x - i1 + G3;
     vec3 x2 = x - i2 + 2.0*G3;
     vec3 x3 = x - 1.0 + 3.0*G3;
     vec4 w, d;
     w.x = dot(x, x);
     w.y = dot(x1, x1);
     w.z = dot(x2, x2);
     w.w = dot(x3, x3);
     w = max(0.6 - w, 0.0);
     d.x = dot(random3(s), x);
     d.y = dot(random3(s + i1), x1);
     d.z = dot(random3(s + i2), x2);
     d.w = dot(random3(s + 1.0), x3);
     w *= w;
     w *= w;
     d *= w;
     return dot(d, vec4(52.0,52.0,52.0,52.0));
}

float simplex3d_fractal(vec3 m, int steps)
{
    float sum = 0.0;
    for (int i = 0; i < steps; ++i)
    {
        float scale = pow(2.0, float(i));
        sum += simplex3d(scale * m) / scale;
    }
    return sum;
}

vec3 flow_texture( vec3 p, float Time, int steps)
{
    // animate initial coordinates
    vec3 p1 = 0.1 * p + vec3(1.0 + Time * 0.0023, 2.0 - Time * 0.0017, 4.0 + Time * 0.0005);
    // distort noise sampling coordinates using the same noise function
    vec3 p2 = p + 8.1 * simplex3d_fractal(p1,steps) + 0.5;
    vec3 p3 = p2 + 4.13 * simplex3d_fractal(0.5 * p2 + vec3(5.0, 4.0, 8.0 + Time * 0.07),steps) + 0.5;

    vec3 ret;
    ret[0] = simplex3d_fractal(p3 + vec3(0.0, 0.0, 0.0 + Time * 0.3),steps);
    ret[1] = simplex3d_fractal(p3 + vec3(0.0, 0.0, 0.2 + Time * 0.3),steps);
    ret[2] = simplex3d_fractal(p3 + vec3(0.0, 0.0, 0.3 + Time * 0.3),steps);

    // scale output & map
    ret = 0.5 + 0.5 * ret;
    ret = smoothstep(vec3(0.15), vec3(0.85), ret);
    return ret;
}

shader SubstanceFlow
(
	point Po = P,
	float Time = 0,
	int steps = 12,
	float scale = 1,
	output vector Out = 0,

)
{
	float iResolution = scale; 
	
	int numSamples = 1;
    vec3 result = vec3(0.0);
	vec2 fragCoord = vec2(Po[0],Po[1]);
     // cheap AA
    for (int x = 0; x < numSamples; ++x)
    {
        for (int y = 0; y < numSamples; ++y)
        {
            vec2 offset = vec2(float(x), float(y)) / float(numSamples);
            vec2 uv = (fragCoord + offset) / iResolution;
            vec3 p = vec3(uv.x,uv.y, Time*0.001);
            result += flow_texture(p * 6.0,Time,steps);
            
        }
    }

    result /= float(numSamples * numSamples);
    Out = vec3(sqrt(result));
    Out = pow(Out,2.2);
}