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

/* 
 * Mandelbulb - in memory of great mathematician Benoit B. Mandelbrot
 * coded by Sylvio Sell - maitag.de - Rostock Germany 2017
 */

CCL_NAMESPACE_BEGIN

ccl_device float svm_mandelbulb(float3 p, float max_iter, float max_pow, float n)
{
    int iter = 0;
	
	float r, theta, phi;	
	float tmp_pow_r_n, tmp_sin_theta_n;

    float x = p.x;
    float y = p.y;
    float z = p.z;
    
    float tmp_xyz = x*x + y*y + z*z;

	while ( tmp_xyz < max_pow && iter < max_iter)
	{
	
		r = sqrtf(tmp_xyz);
		theta = atan2f(sqrtf(x*x + y*y) , z);
		phi = atan2f(y,x);
		
		tmp_pow_r_n = powf(r,n);
		tmp_sin_theta_n = sinf(theta*n);
		
		x = tmp_pow_r_n * tmp_sin_theta_n * cosf(phi*n) + p.x;
		y = tmp_pow_r_n * tmp_sin_theta_n * sinf(phi*n) + p.y;
		z = tmp_pow_r_n * cosf(theta*n) + p.z;

		tmp_xyz = x*x + y*y + z*z;

		iter++;
	}

	return( iter / max_iter );
}

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

	/* Input and Output Sockets */
	uint co_offset, fac_offset, iteration_offset, symmetry_offset;
	uint power_offset, scale_offset;

	decode_node_uchar4(node.y, &co_offset, &fac_offset, &iteration_offset, &symmetry_offset);
	decode_node_uchar4(node.z, &power_offset, &scale_offset, NULL, NULL);

	float3 co = stack_load_float3(stack, co_offset);
	
	float iteration = stack_load_float_default(stack, iteration_offset, node2.x);
	float symmetry  = stack_load_float_default(stack, symmetry_offset,  node2.y);
	float power = stack_load_float_default(stack, power_offset, node2.z);
	float scale = stack_load_float_default(stack, scale_offset, node2.w);
	
	float intensity = svm_mandelbulb(co * scale, iteration, power, symmetry);
	
	if(stack_valid(fac_offset))
		stack_store_float(stack, fac_offset, intensity);
}

CCL_NAMESPACE_END

