/*
 * Copyright (C) 2007 askywhale
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License (version 2)
 *  as published by the Free Software Foundation.
 *
 *  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.
 *
 */
package com.askywhale.rings;

import net.sf.j3d.util.graphics.RGB;
import net.sf.j3d.util.Vector;
/**
 * Rings Raytracer related class
 * 
 * Propose some simple 3d fractals (change kind attribute)
 * 
 * @author askywhale
 *
 */
public class FractalIsoSurface extends IsoSurface {
    
	public static final int MANDLEBROT_Z_CYL = 0;
	public static final int JULIA_DENDRITE_Z_CYL = 2;
	public static final int JULIA_DOUADYS_RABBIT_Z_CYL = 3;
	public static final int JULIA_SAN_MARCO_Z_CYL = 4;
	public static final int JULIA_SIEGEL_DISK_Z_CYL = 5;
	public static final int JULIA_QUATER_D_0_1_0_0 = 6;
	public static final int JULIA_HYPER_D_0_1_0_0 = 7;
		
	public static final int DEFAULT_ITERATION_EXCESS = 1;
	public static final double DEFAULT_THRESOLD = .5; 

	private int kind;
	private int iterationsExcess, iterationsMax;
	private double thresold;
	
    public FractalIsoSurface(Vector location, double size, RGB color, int kind) {
        super(location, size, color);
        this.kind = kind;
        setIterationsExcess(DEFAULT_ITERATION_EXCESS);
        thresold = DEFAULT_THRESOLD;
        discreteShadeMode = true;
    }

    protected double getValue(Vector point) {
    	//double x = point.getX(), y = point.getY(), z = point.getZ();
    	double x = point.x, y = point.y, z = point.z;
    	switch (kind) {
    	case JULIA_DENDRITE_Z_CYL:
    		return thresold-getJulia(0,1,z,Math.hypot(x,y));
    	case JULIA_DOUADYS_RABBIT_Z_CYL:
    		return thresold-getJulia(-.123,.745,z,Math.hypot(x,y));
    	case JULIA_SAN_MARCO_Z_CYL:
    		return thresold-getJulia(0,-.75,z,Math.hypot(x,y));
    	case JULIA_SIEGEL_DISK_Z_CYL:
    		return thresold-getJulia(-.391,-.587,z,Math.hypot(x,y));
    	case JULIA_QUATER_D_0_1_0_0: {
    		double [] c = {0,1,0,0};    	    
    		double [] z0 = {x,y,z,0};
    		return thresold-getJuliaQuaternions(c,z0); }
    	case JULIA_HYPER_D_0_1_0_0: {
    		double [] c = {0,1,0,0};
    		double [] z0 = {x,y,z,0};
    		return thresold-getJuliaHypercomplex(c,z0); }
    	case MANDLEBROT_Z_CYL:
		default:
			return thresold-getMandelbrot(-z,Math.hypot(x,y));
		}    	  	
    }

	public int getKind() {
		return kind;
	}

	public void setKind(int kind) {
		this.kind = kind;
	}
	
	public int getIterationsExcess() {
		return iterationsExcess;
	}

	public void setIterationsExcess(int iterationsExcess) {
		this.iterationsExcess = iterationsExcess;
		recountIterationsMax(); 
	}
	
	public double getThresold() {
		return thresold;
	}

	public void setThresold(double thresold) {
		this.thresold = thresold;
		recountIterationsMax();
	}
	
	protected void recountIterationsMax() {
		iterationsMax = (int)(Math.floor(thresold)+iterationsExcess+1);
	}
	
	protected int getMandelbrot(double cx, double cy){
		int i;
		double zx = 0, zy = 0, nzx, nzy; 
		for(i=0;i<iterationsMax && 
				zx*zx+zy*zy<4;i++) {
			// nz = nē+c
			nzx = zx*zx-zy*zy+cx;
			nzy = 2.0*zx*zy+cy;
			zx = nzx;
			zy = nzy;
		}
		return i;
	}

	// http://framy.free.fr/travaux_mathematiques.htm
	// http://mathworld.wolfram.com/JuliaSet.html
	// http://www.povray.org/documentation/view/3.6.1/280/
	protected int getJulia(double cx, double cy, double z0x, double z0y){
		int i;
		double zx = z0x, zy = z0y, nzx, nzy;
		for(i=0;i<iterationsMax && 
				zx*zx+zy*zy<4;i++) {
			// nz = nē+c
			nzx = zx*zx-zy*zy+cx;
			nzy = 2.0*zx*zy+cy;
			zx = nzx;
			zy = nzy;
		}
		return i;
	}
	
	protected int getJuliaQuaternions(double [] c, double [] z0){
		int i;
		//double [] z = z0.clone();
		double [] z = new double[4];
		z[0] = z0[0];
		z[1] = z0[1];
		z[2] = z0[2];
		z[3] = z0[3];
		double [] nz = new double[4];
		for(i=0;i< iterationsMax &&	z[0]*z[0]+z[1]*z[1]+z[2]*z[2]+z[3]*z[3]<4;i++) {			
			// nz = nē+c 
			/* z1 * z2 = a*e - b*f - c*g- d*h + 
				i (a*f + b*e + c*h - d*g) + 
				j (a*g - b*h + c*e + d*f) + 
				k (a*h + b*g - c*f + d*e)
			*/
			nz[0] = z[0]*z[0] - z[1]*z[1] - z[2]*z[2] - z[3]*z[3] + c[0];
			nz[1] = z[0]*z[1] + z[1]*z[0] + z[2]*z[3] - z[3]*z[2] + c[1];
			nz[2] = z[0]*z[2] - z[1]*z[3] + z[2]*z[0] + z[3]*z[1] + c[2];
			nz[3] = z[0]*z[3] + z[1]*z[2] - z[2]*z[1] + z[3]*z[0] + c[3];
			//z = nz.clone();
			z[0] = nz[0];
			z[1] = nz[1];
			z[2] = nz[2];
			z[3] = nz[3];
		}
		return i;
	}
	
	//http://home.usit.net/~cmdaven/hyprcplx.htm
	protected int getJuliaHypercomplex(double [] c, double [] z0){
		int i;
		//double [] z = z0.clone();
		double [] z = new double[4];
		z[0] = z0[0];
		z[1] = z0[1];
		z[2] = z0[2];
		z[3] = z0[3];
		double [] nz = new double[4];
		for(i=0;i< iterationsMax && 
				z[0]*z[0]+z[1]*z[1]+z[2]*z[2]+z[3]*z[3]<4;i++) {			
			// nz = nē+c
			nz[0] = z[0]*z[0] - z[1]*z[1] - z[2]*z[2] + z[3]*z[3] + c[0];
			nz[1] = z[0]*z[1] + z[1]*z[0] - z[2]*z[3] - z[3]*z[2] + c[1];
			nz[2] = z[0]*z[2] - z[1]*z[3] + z[2]*z[0] - z[3]*z[1] + c[2];
			nz[3] = z[0]*z[3] + z[1]*z[2] + z[2]*z[1] + z[3]*z[0] + c[3];
			//z = nz.clone();
			z[0] = nz[0];
			z[1] = nz[1];
			z[2] = nz[2];
			z[3] = nz[3];
		}
		return i;
	}
}
