/*
 * Copyright 2015, 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 OSL Shader - by Sylvio Sell - maitag.de - Rostock Germany 2015
 * More information: https://github.com/maitag/blender-osl-shader
 */

float mandelbulb_intensity (
    point p,
    int max_iter,  // max number of iterations
    float max_pow, // max power
    float n        // power
   )
{
    // -----------------------------
    int iter = 0;
    
	float r, theta, phi;	
	float tmp_pow_r_n, tmp_sin_theta_n;

    float x = p[0];
    float y = p[1];
    float z = p[2];
    
    float tmp_xyz = x*x + y*y + z*z;

	while ( tmp_xyz < max_pow && iter < max_iter)
	{
	
		r = sqrt(tmp_xyz);
		theta = atan2(sqrt(x*x + y*y) , z);
		phi = atan2(y,x);
		
		tmp_pow_r_n = pow(r,n);
		tmp_sin_theta_n = sin(theta*n);
		
		x = tmp_pow_r_n * tmp_sin_theta_n * cos(phi*n) + p[0];
		y = tmp_pow_r_n * tmp_sin_theta_n * sin(phi*n) + p[1];
		z = tmp_pow_r_n * cos(theta*n) + p[2];

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

		iter++;
	}
    
    return(iter);
}

// ----------------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------------

shader node_mandelbulb (
	int iterations = 12, // number of iterations
	float power = 2.0,   // max power
	float symmetry = 3.0,    // power
	
    float scale = 1.0,

    point position = P,
    output float intensity = 0.0 )
{
    // -----------------------------

    intensity = mandelbulb_intensity( position*scale, iterations, power, symmetry ) / iterations;
}