/*! \file lic.diderot
 *
 * \author John Reppy
 *
 * This Diderot program implements a lic on an iso surface that was discovered by
 * the rcast.diderot program.
 */

/*
 * This code is part of the Diderot Project (http://diderot-language.cs.uchicago.edu)
 *
 * COPYRIGHT (c) 2015 The University of Chicago
 * All rights reserved.
 */

input image(2)[] nDepth;  // depth output from RCast phase
field D = tent ⊛ nDepth;

field#1(2)[2] V = image("data/vfrhand-nohip-smooth.nrrd") ⊛ bspln3;
field#0(2)[] R = image("data/turb2d.nrrd") ⊛ tent;	// noise field

strand LIC (int xi, int yi) {

    vec3 pos0 = [xx, yy, -D([xx,yy])];  // viewing plane is X-Y plane looking toward -Z
    vec3 forw = pos0;
    vec3 back = pos0;
    int step = 0;
    output real sum = R(pos0);

/* what I'd like to be able to write:
    initially {
	if (isNan(sum)) {
	    stabilize;
	}
    }
*/

    update {
	if (isNan(sum)) {
	    sum = 0;
	    stabilize;
	}
        // Euler integration step
        // forw = forw + h*V(forw);
        // back = back - h*V(back);
        // Midpoint method step
        forw += h*V(forw + 0.5*h*V(forw));
        back -= h*V(back - 0.5*h*V(back));
        sum += R(forw) + R(back);
        step += 1;
        if (step == stepNum) {
          // modulate output by velocity at initial position
          sum = |V(pos0)|*sum/real(1 + 2*stepNum);
          stabilize;
        }
    }

} /* strand LIC */


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