/*
 * 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
 * coded by Sylvio Sell - maitag.de - Rostock Germany 2013
 * devel task: https://developer.blender.org/T32305
 * for blender internal renderer equivalent look at: source/blender/render/intern/source/render_texture.c
 */

CCL_NAMESPACE_BEGIN

/* Lyapunov */

ccl_device float svm_lyapunov(float3 p, float lyap_iteration_pre, float lyap_iteration_main, float lyap_param1, float lyap_param2, float lyap_shift, float lyap_balance)
{
	float a = p.x;
	float b = p.y;
	float c = p.z;

	int iter_pre =  floor_to_int(lyap_iteration_pre);
	int iter_main = floor_to_int(lyap_iteration_main);
	float nabla_pre = lyap_iteration_pre - (float)iter_pre;
	float nabla_main = lyap_iteration_main - (float)iter_main;
	float p1 = lyap_param1;
	float p2 = lyap_param2;
	float x = lyap_shift;
	float bal = lyap_balance;
	float x_pre;
	float index = 0.0f;
	float index_pre;
	float ableitung;

	int i;
	int iter = 0;

	/* pre-iteration */
	for(i = 0; i < iter_pre; i++)
	{
			x = p1*sinf(x+a)*sinf(x+a)+p2;
			x = p1*sinf(x+b)*sinf(x+b)+p2;
			x = p1*sinf(x+c)*sinf(x+c)+p2;
	}

	if (nabla_pre != 0.0f)
	{
		x_pre = x;
		x = p1*sinf(x+a)*sinf(x+a)+p2;
		x = p1*sinf(x+b)*sinf(x+b)+p2;
		x = p1*sinf(x+c)*sinf(x+c)+p2;
		x = x*nabla_pre + x_pre*(1.0f-nabla_pre);
	}

	/* main-iteration */
	for(i = 0; i < iter_main; i++)
	{
		x = p1*sinf(x+a)*sinf(x+a)+p2;
		ableitung = 2.0f*p1*sinf(x+a)*cosf(x+a);
		if (ableitung != 0.0f) { index += (logf(fabsf(ableitung))*bal + ableitung*(1.0f-bal)) / 2.0f; iter++; }

		x = p1*sinf(x+b)*sinf(x+b)+p2;
		ableitung = 2.0f*p1*sinf(x+b)*cosf(x+b);
		if (ableitung != 0.0f) { index += (logf(fabsf(ableitung))*bal + ableitung*(1.0f-bal)) / 2.0f; iter++; }

		x = p1*sinf(x+c)*sinf(x+c)+p2;
		ableitung = 2.0f*p1*sinf(x+c)*cosf(x+c);
		if (ableitung != 0.0f) { index += (logf(fabsf(ableitung))*bal + ableitung*(1.0f-bal)) / 2.0f; iter++; }
	}


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

		x = p1*sinf(x+a)*sinf(x+a)+p2;
		ableitung = 2.0f*p1*sinf(x+a)*cosf(x+a);
		if (ableitung != 0.0f) { index += (logf(fabsf(ableitung))*bal + ableitung*(1.0f-bal)) / 2.0f; iter++; }

		x = p1*sinf(x+b)*sinf(x+b)+p2;
		ableitung = 2.0f*p1*sinf(x+b)*cosf(x+b);
		if (ableitung != 0.0f) { index += (logf(fabsf(ableitung))*bal + ableitung*(1.0f-bal)) / 2.0f; iter++; }

		x = p1*sinf(x+c)*sinf(x+c)+p2;
		ableitung = 2.0f*p1*sinf(x+c)*cosf(x+c);
		if (ableitung != 0.0f) { index += (logf(fabsf(ableitung))*bal + ableitung*(1.0f-bal)) / 2.0f; iter++; }

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

	return(index);
}

ccl_device void svm_node_tex_lyapunov(KernelGlobals *kg, ShaderData *sd, float *stack, uint4 node, int *offset)
{
	uint4 node2 = read_node(kg, offset);
	uint4 node3 = read_node(kg, offset);

	/* Input and Output Sockets */
	uint co_offset, fac_offset, positiv_offset, negativ_offset;
	uint shift_offset, balance_offset, iteration_pre_offset, iteration_main_offset;
	uint param1_offset, param2_offset, scale_offset;

	decode_node_uchar4(node.y, &co_offset, &fac_offset, &positiv_offset, &negativ_offset);
	decode_node_uchar4(node.z, &shift_offset, &balance_offset, &iteration_pre_offset, &iteration_main_offset);
	decode_node_uchar4(node.w, &param1_offset, &param2_offset, &scale_offset, NULL);

	float3 co = stack_load_float3(stack, co_offset);
	
	float shift = stack_load_float_default(stack, shift_offset, node2.x);
	float balance = stack_load_float_default(stack, balance_offset, node2.y);
	float iteration_pre = stack_load_float_default(stack, iteration_pre_offset, node2.z);
	float iteration_main = stack_load_float_default(stack, iteration_main_offset, node2.w);

	float param1 = stack_load_float_default(stack, param1_offset, node3.x);
	float param2 = stack_load_float_default(stack, param2_offset, node3.y);
	float scale = stack_load_float_default(stack, scale_offset, node3.z);

	float positiv = 0.0f;
	float negativ = 0.0f;
	float index = svm_lyapunov(co * scale * 4.0f, iteration_pre, iteration_main, param1, param2, shift, balance);
	if (index > 0.0f) positiv = index; else negativ = -index;

	if(stack_valid(fac_offset))
		stack_store_float(stack, fac_offset, 0.5f + index);
	if(stack_valid(positiv_offset))
		stack_store_float(stack, positiv_offset, positiv);
	if(stack_valid(negativ_offset))
		stack_store_float(stack, negativ_offset, negativ);
}

CCL_NAMESPACE_END

