/*
 * 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.
 */

CCL_NAMESPACE_BEGIN

/* Lyapunov */

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

	int iter_pre =  (int)floor(lyap_iteration_pre);
	int iter_main = (int)floor(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 = 0.0f;
	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)); 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)); 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)); 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)); 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)); 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)); iter++; }

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

	return(index);
}

__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 scale_offset, lyap_iteration_pre_offset, lyap_iteration_main_offset, vector_offset, fac_offset, color_offset;
	uint lyap_pos_color_offset, lyap_mid_color_offset, lyap_neg_color_offset;
	uint lyap_pos_scale_offset, lyap_neg_scale_offset, lyap_param1_offset, lyap_param2_offset;

	/* RNA properties */
	uint render_type, fac_type;


	decode_node_uchar4(node.y, &fac_type, &render_type, &color_offset, &fac_offset);
	decode_node_uchar4(node.z, &vector_offset, &scale_offset, &lyap_iteration_pre_offset, &lyap_iteration_main_offset);
	decode_node_uchar4(node.w, &lyap_pos_color_offset, &lyap_mid_color_offset, &lyap_neg_color_offset, NULL);

	decode_node_uchar4(node2.x, &lyap_pos_scale_offset, &lyap_neg_scale_offset, &lyap_param1_offset, &lyap_param2_offset);
	
	float3 vo = stack_load_float3(stack, vector_offset);

	float3 lyap_pos_color = stack_load_float3(stack, lyap_pos_color_offset);
	float3 lyap_mid_color = stack_load_float3(stack, lyap_mid_color_offset);
	float3 lyap_neg_color = stack_load_float3(stack, lyap_neg_color_offset);

	float scale = stack_load_float_default(stack, scale_offset, node2.y);
	float lyap_iteration_pre = stack_load_float_default(stack, lyap_iteration_pre_offset, node2.z);
	float lyap_iteration_main = stack_load_float_default(stack, lyap_iteration_main_offset, node2.w);
	float lyap_pos_scale = stack_load_float_default(stack, lyap_pos_scale_offset, node3.x);
	float lyap_neg_scale = stack_load_float_default(stack, lyap_neg_scale_offset, node3.y);
	float lyap_param1 = stack_load_float_default(stack, lyap_param1_offset, node3.z);
	float lyap_param2 = stack_load_float_default(stack, lyap_param2_offset, node3.w);

	NodeLyapunovFacType lyap_fac_type = (NodeLyapunovFacType)fac_type;
	NodeLyapunovRenderType lyap_render_type = (NodeLyapunovRenderType)render_type;
	float lyap_pos_color_r = lyap_pos_color.x;
	float lyap_pos_color_g = lyap_pos_color.y;
	float lyap_pos_color_b = lyap_pos_color.z;
	float lyap_mid_color_r = lyap_mid_color.x;
	float lyap_mid_color_g = lyap_mid_color.y;
	float lyap_mid_color_b = lyap_mid_color.z;
	float lyap_neg_color_r = lyap_neg_color.x;
	float lyap_neg_color_g = lyap_neg_color.y;
	float lyap_neg_color_b = lyap_neg_color.z;

	float index = svm_lyapunov(vo*scale*4, lyap_iteration_pre, lyap_iteration_main, lyap_param1, lyap_param2);
	float r,g,b;

	// calculate RGBA Values
	if ( index > 0.0f && (lyap_render_type != NODE_LYAPUNOV_NEG)) {
		index *= lyap_pos_scale;
		if (index > 1.0f) { index = 1.0f; }
		r = (lyap_pos_color_r - lyap_mid_color_r)*index+lyap_mid_color_r;
		g = (lyap_pos_color_g - lyap_mid_color_g)*index+lyap_mid_color_g;
		b = (lyap_pos_color_b - lyap_mid_color_b)*index+lyap_mid_color_b;
		//a = (lyap_pos_color_a - lyap_mid_color_a)*index+lyap_mid_color_a;
	}
	else if ( index < 0.0f && (lyap_render_type != NODE_LYAPUNOV_POS)) {
		index *= lyap_neg_scale;
		if (index < -1.0f) { index = -1.0f; }
		r = (lyap_mid_color_r - lyap_neg_color_r)*index+lyap_mid_color_r;
		g = (lyap_mid_color_g - lyap_neg_color_g)*index+lyap_mid_color_g;
		b = (lyap_mid_color_b - lyap_neg_color_b)*index+lyap_mid_color_b;
		//a = (lyap_mid_color_a - lyap_neg_color_a)*index+lyap_mid_color_a;
	}
	else {
		r = lyap_mid_color_r;
		g = lyap_mid_color_g;
		b = lyap_mid_color_b;
		//a = lyap_mid_color_a;
	}
	
	if (lyap_fac_type == NODE_LYAPUNOV_SPREAD) {
		index = 0.5f + index/2.0f;
	}
	else if (lyap_fac_type == NODE_LYAPUNOV_ABS) {
		index = fabsf(index);
	}
	else if (lyap_fac_type == NODE_LYAPUNOV_COLOR) {
		index = (r+g+b)/3.0f;
	}

	if(stack_valid(fac_offset))
		stack_store_float(stack, fac_offset, index);
	if(stack_valid(color_offset))
		stack_store_float3(stack, color_offset, make_float3(r, g, b));
}

CCL_NAMESPACE_END

