float func(float i, float n ,float uParama,float uParamb) {
	return (uParama*pow(sin(i+n),2.0))+uParamb;
}
float deriv(float i, float n ,float uParama,float uParamb) {
	return (((uParama*2.0)*pow(sin(n+i),2.0))*cos(n+i))/sin(n+i);
}

shader node_lyapunov(
	float Start_Index = 1.05691056910569,
	float Pre_Iteration = 1.70731707317073,
	float Main_Iteration = 38.2113821138211,

	float Param_a = 2.5,
	float Param_b = 3.29268292682927,
	
	
	float Balance = 1.4390243902439,
	
	color PosColor = color (0.223529411764706, 0, 0.603921568627451),
	color MidColor = color (0.294117647058824, 0.533333333333333, 0),
	color NegColor = color (0.905882352941176, 0.568627450980392, 0),
	
	point Pos = P,			
	float Scale = 1.0,
	
	output color Color = color (0.0, 0.0, 0.0),		
	output float Index = 0.0,
	output float PosIndex = 0.0,
	output float NegIndex = 0.0			
	)
{
	/* Calculate Lyapunov Index */

	point p = Pos * Scale;
	float i = Start_Index;
	
	int iter_pre =  (int)floor(Pre_Iteration);
	int iter_main = (int)floor(Main_Iteration);
	float iter_main_full = floor(Main_Iteration) * 3.0;
	
	float nabla_pre = Pre_Iteration - (float)iter_pre;
	float nabla_main = Main_Iteration - (float)iter_main;

	/* Pre-iteration */
	
	for(int j = 0; j < iter_pre; j++) {
		i =  func(i, p[0],Param_a,Param_b);
	i =  func(i, -p[1],Param_a,Param_b);
	i =  func(i, p[2],Param_a,Param_b);
	
	}

	if (nabla_pre != 0.0) {
		float i_pre = i;

		i =  func(i, p[0],Param_a,Param_b);
	i =  func(i, -p[1],Param_a,Param_b);
	i =  func(i, p[2],Param_a,Param_b);
	

		i = i*nabla_pre + i_pre*(1.0-nabla_pre);
	}

	/* Main-iteration */
	
	for(int j = 0; j < iter_main; j++) {
		i =  func(i, p[0],Param_a,Param_b);
	Index += (log(abs(deriv(i, p[0],Param_a,Param_b)))*Balance + (deriv(i, p[0],Param_a,Param_b))*(1.0-Balance)) / 2.0;
i =  func(i, -p[1],Param_a,Param_b);
	Index += (log(abs(deriv(i, -p[1],Param_a,Param_b)))*Balance + (deriv(i, -p[1],Param_a,Param_b))*(1.0-Balance)) / 2.0;
i =  func(i, p[2],Param_a,Param_b);
	Index += (log(abs(deriv(i, p[2],Param_a,Param_b)))*Balance + (deriv(i, p[2],Param_a,Param_b))*(1.0-Balance)) / 2.0;
			
	}
	

	if (nabla_main == 0.0) {
		Index = (iter_main_full != 0) ? Index/iter_main_full : 0.0;
	}
	else {
		float index_pre = (iter_main_full != 0) ? Index/iter_main_full : 0.0;
		i =  func(i, p[0],Param_a,Param_b);
	Index += (log(abs(deriv(i, p[0],Param_a,Param_b)))*Balance + (deriv(i, p[0],Param_a,Param_b))*(1.0-Balance)) / 2.0;
i =  func(i, -p[1],Param_a,Param_b);
	Index += (log(abs(deriv(i, -p[1],Param_a,Param_b)))*Balance + (deriv(i, -p[1],Param_a,Param_b))*(1.0-Balance)) / 2.0;
i =  func(i, p[2],Param_a,Param_b);
	Index += (log(abs(deriv(i, p[2],Param_a,Param_b)))*Balance + (deriv(i, p[2],Param_a,Param_b))*(1.0-Balance)) / 2.0;


		iter_main_full = iter_main_full + 3.0;

		Index = (iter_main_full != 0) ? Index/iter_main_full : 0.0;
		Index = Index*nabla_main + index_pre*(1.0-nabla_main);
	}

	/* separate output */
	if (Index > 0.0) {
		PosIndex = Index;
		Color = Index * (PosColor - MidColor) + MidColor;
	}
	else {
		NegIndex = -Index;
		Color = Index * (MidColor - NegColor) + MidColor;
	}		
}