/*
=======================================================

SIGGRAPH 2019 submission                     paper_0306

=======================================================

This code is for reviewing only.

              !!! Do not distribute !!!

=======================================================

HowTo:
------

Copy-paste the code into https://shadertoy.com/new/

=> Please be careful not to save the shader in shadertoy

Click on the play button at the bottom left of the 
code window. 

Comment/uncomment lines at the bottom to try different
noises and profiles.

Thank you!

=======================================================
*/

# define M_PI 3.14159265358979323846

//phasor noise parameters
float _f = 25.0;
float _b = 8.0;
float _o = 1.0;
float _kr;
int _impPerKernel = 32;
int _seed = 1;


vec2 phasor(vec2 x, float f, float b, float o, float phi)
{
    float a = exp(-M_PI * (b * b) * ((x.x * x.x) + (x.y * x.y)));
    float s = sin (2.0* M_PI * f  * (x.x*cos( o) + x.y * sin(o)) + phi);
    float c = cos (2.0* M_PI * f  * (x.x*cos( o) + x.y * sin(o)) + phi);
    return vec2(a*c,a*s);
}


///////////////////////////////////////////////
//prng
///////////////////////////////////////////////

int N = 15487469;
int x_;
void seed(int s){x_ = s;}
int next() { x_ *= 3039177861; x_ = x_ % N;return x_; }
float uni_0_1() {return  float(next()) / float(N);}
float uni(float min, float max){ return min + (uni_0_1() * (max - min));}
int poisson(float mean)
{
	float g_ = exp(-mean);
	int em = 0;
	float t = uni_0_1();
	while (t > g_) {
	  ++em;
	  t *= uni_0_1();
	}
	return em;
}


int morton(int x, int y)
{
  int z = 0;
  for (int i = 0 ; i < 32* 4 ; i++) {
    z |= ((x & (1 << i)) << i) | ((y & (1 << i)) << (i + 1));
  }
  return z;
}


void init_noise()
{
    _kr = sqrt(-log(0.05) / M_PI) / _b;
}


vec2 cell(ivec2 ij, vec2 uv, float f, float b, float o)
{
	int s= morton(ij.x,ij.y) + 333;
	s = s==0? 1: s +_seed;
	seed(s);
	int impulse  =0;
	int nImpulse = _impPerKernel;
	float  cellsz = 2.0 * _kr;
	vec2 noise = vec2(0.0);
	while (impulse <= nImpulse){
		vec2 impulse_centre = vec2(uni_0_1(),uni_0_1());
		vec2 d = (uv - impulse_centre) *cellsz;
		float rp = uni(0.0,2.0*M_PI);
		noise += phasor(d, f, b ,o,rp );
		impulse++;
	}
	return noise;
}

vec2 eval_noise(vec2 uv, float f, float b, float o)
{   
	float cellsz = 2.0 *_kr;
	vec2 _ij = uv / cellsz;
	ivec2  ij = ivec2(_ij);
	vec2  fij = _ij - vec2(ij);
	vec2 noise = vec2(0.0);
	for (int j = -2; j <= 2; j++) {
		for (int i = -2; i <= 2; i++) {
			ivec2 nij = ivec2(i, j);
			noise += cell(ij + nij , fij - vec2(nij),f,b,o );
		}
	}
    return noise;
}

float PWM(float x, float r)
{
	return mod(x,2.0*M_PI)> 2.0*M_PI *r ? 1.0 : 0.0; 
}

float square(float x)
{
  return PWM(x,0.5);   
}

float sawTooth(float x)
{
	return mod(x,2.0*M_PI)/(2.0*M_PI);
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
  vec2 uv = fragCoord/iResolution.y;
  uv.y=-uv.y;
  
  init_noise();
  
  float o = _o + uv.x*M_PI/2.0;
  vec2 phasorNoise = eval_noise(uv,_f,_b,o);
  
	float Phi = atan(phasorNoise.y,phasorNoise.x);
  
  fragColor = vec4(vec3(PWM(Phi,0.5+0.4*cos(iTime)*(0.1+0.9*uv.x))),1.0); // animated PWM
  //fragColor = vec4(vec3(sin(Phi)*0.4+0.5),1.0); // sinewave  
  //fragColor = vec4(vec3(sawTooth(Phi)),1.0); //sawtooth
  //fragColor = vec4(vec3(square(Phi)),1.0); // square wave
}