/*
 * Copyright 2011, Blender Foundation.
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software Foundation,
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

 
/* 
 * Lyapunov Shader - in memory of great mathematician Aleksandr Mikhailovich Lyapunov
 * Original code: Sylvio Sell - maitag.de - Rostock Germany 2013
 * recoded 2015 for better node integration and artists control
 * devel task: https://developer.blender.org/T32305
 * thanks to Thomas Dinges for initial OSL port
 */

 
/* ------- function that describes behavior of dynamic system ------- */

float func(float x, float p, float p1, float p2)
{
    return p1*sin(x+p)*sin(x+p)+p2;
}
float deriv(float x, float abc, float p1, float p2)
{
    return 2.0*p1*sin(x+abc)*cos(x+abc);
}


/* ------- sequence (todo: customizing by parameter;) ------- */

void pre_step(output float x, point p, float p1, float p2)
{
    for(int i = 0; i < 3; i++) {
        x = func(x,p[i],p1,p2);
    }
}

void main_step(output float index, output int iter, output float x, point p, float p1, float p2, float balance)
{
    for(int i = 0; i < 3; i++) {
        x = func(x,p[i],p1,p2);
        index += (  log(fabs(deriv(x,p[i],p1,p2)))*balance + deriv(x,p[i],p1,p2)*(1.0-balance)  ) / 2.0;
        iter++;
    }
}


/* ------- calculates Lyapunov index ------- */

float lyapunov(point p, float x_start, float iteration_pre, float iteration_main, float p1, float p2, float balance)
{
    float x = x_start;
	int iter_pre =  (int)floor(iteration_pre);
	int iter_main = (int)floor(iteration_main);
	float nabla_pre = iteration_pre - (float)iter_pre;
	float nabla_main = iteration_main - (float)iter_main;

	float index = 0.0;
	float ableitung = 0.0;
	int iter = 0;

	/* Pre-iteration */
	for(int i = 0; i < iter_pre; i++) {
		pre_step(x,p,p1,p2);
	}

	if (nabla_pre != 0.0) {
		float x_pre = x;
		pre_step(x,p,p1,p2);
		x = x*nabla_pre + x_pre*(1.0-nabla_pre);
	}

	/* Main-iteration */
	for(int i = 0; i < iter_main; i++) {
		main_step(index, iter, x, p, p1, p2, balance);
	}

	if (nabla_main == 0.0) {
		index = (iter != 0) ? index/(float)(iter) : 0.0;
	}
	else {
		float index_pre = (iter != 0) ? index/(float)(iter) : 0.0;

		main_step(index, iter, x, p, p1, p2, balance);

		index = (iter != 0) ? index/(float)(iter) : 0.0;
		index = index*nabla_main + index_pre*(1.0-nabla_main);
	}

	return index;
}


/* ------------------ SHADER MAIN ---------------------- */

shader node_lyapunov(
	float Shift = 0.0,
	float Balance = 0.0,
	float Pre_Iteration = 0.0,
	float Main_Iteration = 1.0,
	float Param1 = 2.0,
	float Param2 = 2.0,
	float Scale = 1.0,
	point Pos = P,
	output float Fac = 0.0,
	output float Positiv = 0.0,
	output float Negativ = 0.0)
{
	/* Calculate Lyapunov Index */
	float index = lyapunov(Pos*Scale, Shift, Pre_Iteration, Main_Iteration, Param1, Param2, Balance);

	/* separate output */
	if (index > 0.0) {
		Positiv = index;
	}
	else {
		Negativ = -index;
	}

    Fac = 0.5 + index; 

}