// lic-turb.diderot
// 
// demo of line integral convolution (LIC) on 2D turbulent flow
//
// process output with:
// unu reshape -i lic-turb2d.txt -s 3 1020 561 | unu quantize -b 8  -o lic-turb2d.ppm

int imgSizeX = 1020;
int imgSizeY = 561;
real h = 0.005;    // step size of integration
int stepNum = 10; // take this many steps both upstream and downstream

field#2(2)[2] boundary = load("../../data/turb2d.nrrd") ⊛ bspln3;
field#2(2)[2] V = load("../../data/v.nrrd") ⊛ bspln3;
field#0(2)[] R = load("../../data/R.nrrd") ⊛ tent;

strand LIC (int xi, int yi) {
	real xx = lerp(0.04, 6.76, 0.0, real(xi), real(imgSizeX)-1.0);
	real yy = lerp(0.04, 3.7, 0.0, real(yi), real(imgSizeY)-1.0);
	vec2 pos0 = [xx,yy];
	vec2 forw = pos0;
	vec2 back = pos0;
	real sum = R(pos0);
	output vec3 RGB = [0.0, 0.0, 0.0];
	real l = 0.0;
	int step = 0;
	int num = 1;
	real crl = 0.0;

	update {
		if (inside(forw, boundary)) {
			forw += h*normalize(V(forw + 0.5*h*normalize(V(forw))));
			sum += R(forw);
			num += 1;
		}
		if (inside(back, boundary)) {
			back -= h*normalize(V(back - 0.5*h*normalize(V(back))));
			sum += R(back);
			num += 1;
		}
		step += 1;
		if ((!inside(forw, boundary) && !inside(back, boundary)) || step == stepNum) {
			sum = pow(|V(pos0)|,0.3)*sum/real(num);
			tensor[2,2] M = ∇⊗V(pos0);
			crl = M[0,1] - M[1,0];
			real s = lerp(0.0, 1.0, 0.0, |crl|, 132.0);
			l = lerp(0.0, 1.0, -1.6, sum, 2.6);
			real c = (1.0 - |2.0 * l - 1.0|) * s;
			if(crl < 0.0) {
				real x = c;
				real m = l - c / 2.0;
				RGB = [c + m, 0.0 + m, x + m];
			}
			else {
				real x = 0.0;
				real m = l - c / 2.0;
				RGB = [0.0 + m, c + m, x + m];
			}
			stabilize;
		}
	}
}

initially [ LIC(xi, yi) | yi in 0..(imgSizeY-1), xi in 0..(imgSizeX-1) ];