////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////
//// ENSEMBLE DE MANDELBROT ET DE JULIA
////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
	
	xdel(winsid())
	tohome
	clear
	clc
	
	global Fractale
	
////////////////////////////////////////////////////////////////////////////////
//
//   FONCTION D'ACTUALISATION
//
////////////////////////////////////////////////////////////////////////////////

function actualisation()
	global Fractale, bandeau_information
	
	select Fractale.Fonction
	case "Carré"               then Fractale.TexteFonction = [' de la fontion z -> z^2'];
	case "Cube"                then Fractale.TexteFonction = [' de la fontion z -> z^3'];
	case "Quadrique"           then Fractale.TexteFonction = [' de la fontion z -> z^4'];
	case "Quinte"              then Fractale.TexteFonction = [' de la fontion z -> z^5'];
	
	case "Cardan"              then Fractale.TexteFonction = [' pour la suite de Cardan'];
	
	case "Sinus"               then Fractale.TexteFonction = [' de la fonction sinus'];
	case "Cosinus"             then Fractale.TexteFonction = [' de la fonction cosinus'];
	case "SinusHyperbolique"   then Fractale.TexteFonction = [' de la fonction sinus hyperbolique'];
	case "CosinusHyperbolique" then Fractale.TexteFonction = [' de la fonction cosinus hyperbolique'];
	
	case "Exponentielle"       then Fractale.TexteFonction = [' de la fonction exponentielle'];
	case "Logarithme"          then Fractale.TexteFonction = [' de la fonction logarithme'];
	case "Autopuissance"       then Fractale.TexteFonction = [' de la fonction z -> z^z'];	
	end
	
	//select Fractale.Fonction
	//case "Carré"               then Fractale.NomFichier = '2';
	//case "Cube"                then Fractale.NomFichier = '3';
	//case "Quadrique"           then Fractale.NomFichier = '4';
	//case "Quinte"              then Fractale.NomFichier = '5';
	
	//case "Cardan"              then Fractale.NomFichier = 'Cardan';
	
	//case "Sinus"               then Fractale.NomFichier = 'Sin';
	//case "Cosinus"             then Fractale.NomFichier = 'Cos';
	//case "SinusHyperbolique"   then Fractale.NomFichier = 'Sh';
	//case "CosinusHyperbolique" then Fractale.NomFichier = 'Ch';
	
	//case "Exponentielle"       then Fractale.NomFichier = 'Exp';
	//case "Logarithme"          then Fractale.NomFichier = 'Ln';
	//case "Autopuissance"       then Fractale.NomFichier = 'zZ';
	//end
	
	//bandeau_information.String = "Ensemble de Mandelbrot" + Fractale.TexteFonction;
	//bandeau_information.String = bandeau_information.String + string(Fractale.Dimension) + " x " + string(Fractale.Dimension) + " points | ";
	//bandeau_information.String = bandeau_information.String + string(Fractale.Iterations) + " itérations";
	bandeau_information.String = string(Fractale.Dimension) + " x " + string(Fractale.Dimension) + " points | " + string(Fractale.Iterations) + " itérations";
endfunction

////////////////////////////////////////////////////////////////////////////////
//
//	INITIALISATION DES PARAMETRES ANALYTIQUES
//
////////////////////////////////////////////////////////////////////////////////
	
	// Fonction d'itération par défaut
	Fractale.Ensemble = 'Mandelbrot';
	
	// Fonction d'itération par défaut
	Fractale.Fonction = 'Carré';
	
	// Nombre d'itérations par défaut
	Fractale.Iterations = 1000;
	
	// Dimension de la matrice des affixes par défaut
	Fractale.Dimension = 400;
	
	// Valeurs de la constante (C)
if (Fractale.Ensemble = 'Julia')
	c = [-1.3 + %i*0.1; 0.27334 - %i*0.00742; 0.1103103 - %i*0.67037; - 0.488679 - %i*0.56790; - 0.561321 + %i*0.64100; - 0.745429 + %i*0.11308];
	
	select Fractale.Fonction 
	case "Carré"             then (c = [-1.3 + %i*0.1; 0.27334 - %i*0.00742; 0.1103103 - %i*0.67037; - 0.488679 - %i*0.56790; - 0.561321 + %i*0.64100; - 0.745429 + %i*0.11308])
	case "Cube"              then (c = - 0.5 + %i*0.13)
	
	case "Sinus"             then (c = - 0.01; e = [0; %i*0.1; %i*0.4])
	case "SinusHyperbolique" then (c = - 0.01; e = [0; %i*0.1; %i*0.4])
	
	case "Exponentielle"     then (c = - 0.01; e = [0; %i*0.1; %i*0.4])
	case "Autopuissance"     then (c = - 0.01)
	end
end
	
	// Bornes du tracé
if (Fractale.Ensemble = 'Julia')
	select Fractale.Fonction
	case "Carré"				then x_min = - 2;		x_max = 1;		y_min = - 1.5;		y_max = 1.5;
	case "Cube"					then x_min = - 1.2;		x_max = 1.2;	y_min = - 1.5;		y_max = 1.5;
	case "Quinte"				then x_min = - 1.2;		x_max = 1.2;	y_min = - 1.5;		y_max = 1.5;
	case "Monôme"				then x_min = - 1.2;		x_max = 1.2;	y_min = - 1.5;		y_max = 1.5;
	
	case "Héron"				then x_min = - 2;		x_max = 1;		y_min = - 1.5;		y_max = 1.5;
	case "Cardan"				then x_min = - 1.2;		x_max = 1.2;	y_min = - 1.2;		y_max = 1.2;
	case "Rationnel"			then x_min = - 1.2;		x_max = 1.2;	y_min = - 1.5;		y_max = 1.5;
	
	case "Sinus"				then x_min = - 2;		x_max = 1;		y_min = - 1.5;		y_max = 1.5;
	case "Cosinus"				then x_min = - 2;		x_max = 1;		y_min = - 1.5;		y_max = 1.5;
	case "SinusHyperbolique"	then x_min = - 2;		x_max = 1;		y_min = - 1.5;		y_max = 1.5;
	case "CosinusHyperbolique"	then x_min = - 2;		x_max = 1;		y_min = - 1.5;		y_max = 1.5;
	
	case "Exponentielle"		then x_min = - 2;		x_max = 1;		y_min = - 1.5;		y_max = 1.5;
	case "Logarithme"			then x_min = - 2;		x_max = 1;		y_min = - 1.5;		y_max = 1.5;
	case "Autopuissance"		then x_min = - 2;		x_max = 1;		y_min = - 1.5;		y_max = 1.5;
	
	else							x_min = - 2;		x_max = 1;		y_min = - 1.5;		y_max = 1.5;
	//else							x_min = 0.2675;		x_max = 0.2685;	y_min = 0.591;		y_max = 0.592;
	end
end
	
	// Valeur de la matrice des constantes (C)
	x_min = - 2;		x_max = 1;		y_min = - 1.5;		y_max = 1.5;
	
	cx = linspace(x_min, x_max, Fractale.Dimension);
	cy = linspace(y_min, y_max, Fractale.Dimension);
	[CX, CY] = meshgrid(cx, cy);
	C = CX + %i*CY;
	
	// Valeur initiale de la matice des affixes (Z0)
	x0_min = -1.3;		x0_max = 1.3;	y0_min = - 1.3;		y0_max = 1.3;
	
	x0 = linspace(x0_min, x0_max, Fractale.Dimension);
	y0 = linspace(y0_min, y0_max, Fractale.Dimension);
	[X0, Y0] = meshgrid(x0, y0);
	Z0 = X0 + %i*Y0;
	
	stacksize("max");

////////////////////////////////////////////////////////////////////////////////
//
//	INITIALISATION DES PARAMETRES GRAPHIQUES
//
////////////////////////////////////////////////////////////////////////////////
	
	Fractale.NomFigure = "Ensembles de Mandelbrot et de Julia";
	Fractale.CouleurFond = [0.9 0.9 0.9];	
	
	Fractale.Algorithme = 'Affixes';
	//Fractale.Algorithme = 'Indices';
	
	//Fractale.Graphe = 'Spectre';
	Fractale.Graphe = 'Contraste';
	//Fractale.Graphe = 'Gradation';
	Fractale.Couleurs = 1000;	
	
	Fractale.NomFichier = '';
	
//	BOGUE : bandeau_information non encore déclaré-affecté !
//	actualisation();
//	mprintf(Fractale.NomFigure + Fractale.TexteFonction + '\n\n');

////////////////////////////////////////////////////////////////////////////////
//
//   FONCTION D'ITERATION
//
////////////////////////////////////////////////////////////////////////////////

function Z1 = f(Z, C)
	global Fractale
	
	select Fractale.Fonction
	case "Carré"				then Z1 = Z.^2 + C;
	case "Cube"					then Z1 = Z.^3 + C;
	case "Quadrique"			then Z1 = Z.^4 + C;
	case "Quinte"				then Z1 = Z.^5 + C;
	
	case "Cardan"				then Z1 = Z.^3 - Z + C;
	
	case "Sinus"				then Z1 = sin(Z) + C;
	case "Cosinus"				then Z1 = cos(Z) + C;
	case "SinusHyperbolique"	then Z1 = sinh(Z) + C;
	case "CosinusHyperbolique"	then Z1 = cosh(Z) + C;
	
	case "Exponentielle"		then Z1 = exp(Z) + C;
	case "Logarithme"			then Z1 = log(Z) + C;
	case "Autopuissance"		then Z1 = Z.^Z + C;
	end
endfunction

////////////////////////////////////////////////////////////////////////////////
//
//	CALCUL DE LA GRILLE DES INDICES
//
//	La matrice-ligne J stocke les indices des points dans l'ensemble de Mandelbrot.
//	Au démarrage J contient les indices de 1 à Fractale.Dimension*Fractale.Dimension, 
//	ce qui signifie que tous les points sont présumés dans l'ensemble.
//
//	A chaque itération les points dont la suite aura dépassé le module critique  
//	seront otés en affectant l'entrée correspondant de J à zéro.
//	Les indices des points dont la suite n'a pas dépassé le module critique 
//	sont calculés avec l'assertion L = J(J>0)
//
//	Au cours du déroulement de l'algorithme le nombre de points diminuera 
//	comme la taille de la matrice L.
//	Les calculs n'étant effectués que sur les entrées définies dans L, 
//	chaque itération sera donc plus rapide que la précédente.
//
////////////////////////////////////////////////////////////////////////////////	

function R = indices(C, N_Iterations)
	// Matrices des affixes Z dans l'expression Z = Z²+C
	Z = zeros(Fractale.Dimension, Fractale.Dimension);
	
	// Matrice des modules
	W = zeros(Fractale.Dimension, Fractale.Dimension);
	
	// Matrice des indices k pour lesquels le module excède 2
	R = - ones(Fractale.Dimension, Fractale.Dimension);
	
	// Indices des suites n'ayant pas divergé
	J = 1:Fractale.Dimension*Fractale.Dimension;
	
	// Itération
	for k = 0:N_Iterations
		// Actualise la matrice des affixes pour les suites n'ayant pas divergé
		L = J(J > 0);
		Z(L) = f(Z(L), C(L));
		
		// Détermine les suites ayant divergé (module de Z(k)supérieur à 2)
		W(L) = abs(Z(L));
		M = find(W(L)>2);
		R(L(M)) = k;
		J(L(M)) = 0;
		
		// Indicateur de progression
		mprintf(['\tNombre d''iterations effectuees: %d/%d\n'], k, Fractale.Iterations)
	end
	// R = R';
endfunction

////////////////////////////////////////////////////////////////////////////////
//
//	CALCUL DE LA GRILLE DES AFFIXES
//
////////////////////////////////////////////////////////////////////////////////

function Z = affixes(C)
	J = 1:Fractale.Dimension*Fractale.Dimension;
	Z = zeros(Fractale.Dimension, Fractale.Dimension);
	
	for k = 1:Fractale.Iterations
		// Ecrètement des affixes ( accélération du calcul )
		Liste = J(J > 0);
		
		// Traceur
		//disp(size(Z)); disp(size(C));
		
		// Calcul des nouvelles affixes
		Z(Liste) = f(Z(Liste), C(Liste));
		
		// Ecrètement des affixes ( accélération du calcul )
		M = find(abs(Z(Liste)) > 2);
		J(Liste(M)) = 0;
		
		// Indicateur de progression
		//mprintf(['\tNombre d''iterations effectuees: %d/%d\n'], k, Fractale.Iterations)
		waitbar(k/Fractale.Iterations, winH);
	end
endfunction

////////////////////////////////////////////////////////////////////////////////
//
//   TRACE DU BASSIN DE CONVERGENCE
//
////////////////////////////////////////////////////////////////////////////////

function traceFractale(R, c_max)
	f.color_map = jetcolormap(c_max);
	f.children.isoview = "on";
	f.children.axes_visible = ["off" "off" "off"];
	
	a = get("traceMatrice");
	
	D = R;
	
	k = find(R <> -1);
	D(k) = floor(R(k)/max(R(k))*c_max);
	
	k = find(R == -1);
	D(k) = c_max;
	
	Matplot(D);
endfunction
function traceFractale1(X, Y, Z)
	a = get("traceMatrice");
	
	a.children(1).data = abs(Z);
endfunction
function traceFractale2(X, Y, Z)
	a = get("traceMatrice");
	
	a.children(1).data = exp(- abs(Z))
	// shading flat
endfunction

///////////////////////////////////////////////////////////////////////////////
//
//	Fonction de tracé de l'ensemble de Mandelbrot
//
///////////////////////////////////////////////////////////////////////////////

function mandelbrot_trace()
	global Fractale
	
	mprintf(['\tTrace de l''ensemble de Mandelbrot\n'])
	
	// Matrice d'affixes de départ
	cx = linspace(x_min, x_max, Fractale.Dimension);
	cy = linspace(y_min, y_max, Fractale.Dimension);
	[CX, CY] = meshgrid(cx, cy);
	C = CX + %i*CY;
	
	// Calcul de la matrice d'affixes
	winH = waitbar("Calcul des affixes");
	
	tic();	
	if (Fractale.Algorithme == "Affixes"), Z = affixes(C); end
	if (Fractale.Algorithme == "Indices"), R = indices(C); end
	
	Fractale.DureeCalcul = toc();
	
	// mprintf("\t\n");
	// mprintf("\tTemps écoulé = %f (s)\n", Fractale.DureeCalcul);
	
	// Fractale.PPS = floor(Fractale.Dimension*Fractale.Dimension/Fractale.DureeCalcul);
	// mprintf("\t%d pixels à la seconde\n", Fractale.PPS);
	
	close(winH);
	
	// Sauvegarde de la matrice calculée
	if (Fractale.Algorithme == "Affixes")
		save("Mandelbrot" + Fractale.NomFichier + ".sod", "X", "Y", "Z")
	end
	
	// Tracé de l'ensemble de Mandelbrot
	a = get("traceMatrice");
	
	//a.data_bounds = [x_min, y_min; x_max, y_max];
	//a.data_bounds = [1, Fractale.Dimension; 1, Fractale.Dimension];
	a.children(1).data = abs(Z);
	
	//select Fractale.Graphe
	//case "Spectre"			then traceFractale(R, Fractale.Couleurs)
	//case "Contraste"			then traceFractale1(X, Y, Z)
	//case "Gradation"			then traceFractale2(X, Y, Z)
	//end
endfunction

///////////////////////////////////////////////////////////////////////////////
//
//	Fonction de tracé de l'ensemble de Julia
//
///////////////////////////////////////////////////////////////////////////////
function julia_trace()
	global Fractale
	
	mprintf(['\tSélection d''une affixe pour tracer son ensemble de Julia\n'])
	
	a = get("traceCourbe");
	
	[ibutton, xcoord, ycoord, iwin, cbmenu] = xclick()
	
	plot(xcoord, ycoord, 'rx')
	
	mprintf(['\tTrace de l''ensemble de Mandelbrot\n'])
endfunction

///////////////////////////////////////////////////////////////////////////////
//
//	Fonction d'export de la figure
//
///////////////////////////////////////////////////////////////////////////////

function figure_impression()
	global Fractale
	
	xs2png(interface, "Mandelbrot" + Fractale.NomFichier)
	// xs2png(interface, "Mandelbrot" + Fractale.NomFichier + "_Gradation")
endfunction	

///////////////////////////////////////////////////////////////////////////////
//
//	Initialisation de l'interface
//	figure > axes > tracé (courbe(s), matrice ...)
//
///////////////////////////////////////////////////////////////////////////////

//	Création de la figure
	interface = figure('figure_position', [0, 0])
	interface.Figure_size = [700, 700];
	interface.Auto_resize = 'on';
	interface.Figure_name = Fractale.NomFigure;
	interface.Backgroundcolor = Fractale.CouleurFond;
	
	delmenu(interface.figure_id, gettext('File'))
	delmenu(interface.figure_id, gettext('?'))
	delmenu(interface.figure_id, gettext('Tools'))
	toolbar(interface.figure_id, 'off')
	
//	Menu de sélection de la fonction d'itération
	mF = uimenu(interface, 'label', 'Fonction');
	
	mFZ2 = uimenu(mF, 'label', 'z -> z^2 + c',     'callback', "Fractale.Fonction = ''Carré'';               actualisation()")
	mFZ3 = uimenu(mF, 'label', 'z -> z^3 + c',     'callback', "Fractale.Fonction = ''Cube'';                actualisation()")
	mFZ4 = uimenu(mF, 'label', 'z -> z^4 + c',     'callback', "Fractale.Fonction = ''Quadrique'';           actualisation()")
	mFZ5 = uimenu(mF, 'label', 'z -> z^5 + c',     'callback', "Fractale.Fonction = ''Quinte'';              actualisation()")
	
	mF01 = uimenu(mF, 'label', ' ')
	
	mFCa = uimenu(mF, 'label', 'z -> z^3 - z + c', 'callback', "Fractale.Fonction = ''Cardan'';              actualisation()")
	
	mF02 = uimenu(mF, 'label', ' ')
	
	mFS =  uimenu(mF, 'label', 'z -> sin(z) + c',  'callback', "Fractale.Fonction = ''Sinus'';               actualisation()")
	mFC =  uimenu(mF, 'label', 'z -> cos(z) + c',  'callback', "Fractale.Fonction = ''Cosinus'';             actualisation()")
	mFSh = uimenu(mF, 'label', 'z -> sh(z) + c',   'callback', "Fractale.Fonction = ''SinusHyperbolique'';   actualisation()")
	mFCh = uimenu(mF, 'label', 'z -> ch(z) + c',   'callback', "Fractale.Fonction = ''CosinusHyperbolique''; actualisation()")
	
	mF03 = uimenu(mF, 'label', ' ')
	
	mFEx = uimenu(mF, 'label', 'z -> exp(z) + c',  'callback', "Fractale.Fonction = ''Exponentielle'';       actualisation()")
	mFLn = uimenu(mF, 'label', 'z -> ln(z) + c',   'callback', "Fractale.Fonction = ''Logarithme'';          actualisation()")
	mFZz = uimenu(mF, 'label', 'z -> z^z + c',     'callback', "Fractale.Fonction = ''Autopuissance'';       actualisation()")
	
//	Menu de sélection du nombre de points
	mD = uimenu(interface, 'label', 'Dimensions');
	mD1 = uimenu(mD, 'label', ' 200 x  200', 'callback', "Fractale.Dimension =  200; actualisation()")
	mD2 = uimenu(mD, 'label', ' 400 x  400', 'callback', "Fractale.Dimension =  400; actualisation()")
	mD3 = uimenu(mD, 'label', ' 600 x  600', 'callback', "Fractale.Dimension =  600; actualisation()")
	mD4 = uimenu(mD, 'label', ' 800 x  800', 'callback', "Fractale.Dimension =  800; actualisation()")
	mD5 = uimenu(mD, 'label', '1000 x 1000', 'callback', "Fractale.Dimension = 1000; actualisation()")

//	Menu de sélection du nombre d'itérations
	mI = uimenu(interface, 'label', 'Itérations');
	mI1 = uimenu(mI, 'label', '  20', 'callback', "Fractale.Iterations =   20; actualisation()")
	mI2 = uimenu(mI, 'label', '  50', 'callback', "Fractale.Iterations =   50; actualisation()")
	mI3 = uimenu(mI, 'label', ' 100', 'callback', "Fractale.Iterations =  100; actualisation()")
	mI4 = uimenu(mI, 'label', '1000', 'callback', "Fractale.Iterations = 1000; actualisation()")

//	Initialisation du tracé
	a = newaxes();
	
	a.axes_bounds = [0, 0, 0.75, 0.75];
	a.margins = [0, 0, 0, 0];
	
	a.data_bounds = [x_min, y_min; x_max, y_max];
	a.axes_visible = "off";
	a.x_location = "top";
	a.y_location = "left";
	a.box = "off";
	
	a.tag = "traceMatrice";
	
	Matplot(zeros(Fractale.Dimension, Fractale.Dimension), "021")
	
	
//	Bandeau d'information
	bandeau_information = uicontrol(interface, "style", "text");
	bandeau_information.Position = [50 50 385 50];
	bandeau_information.FontSize = 15;
	bandeau_information.FontWeight = "bold";
	bandeau_information.HorizontalAlignment = "center";
	bandeau_information.String = string(Fractale.Dimension) + "x" + string(Fractale.Dimension) + " points | ";
	bandeau_information.String = bandeau_information.String + string(Fractale.Iterations) + " itérations";
	bandeau_information.BackgroundColor = [0.7 0.7 0.7];
	
	global bandeau_information
	
//	Commande du tracé de l'ensemble de Mandelbrot
	mandelbrot_bouton = uicontrol(interface, "style", "pushbutton");
	mandelbrot_bouton.Units = "normalized";
	mandelbrot_bouton.Position = [0.8 0.9 0.2 0.1];
	mandelbrot_bouton.FontSize = 15;
	mandelbrot_bouton.String = "Mandelbrot";
	mandelbrot_bouton.BackgroundColor = Fractale.CouleurFond;
	mandelbrot_bouton.Callback = "mandelbrot_trace";
	mandelbrot_bouton.Relief = "raised";

//	Agrandissement de l'ensemble de Mandelbrot
	agrandissement_bouton = uicontrol(interface, "style", "pushbutton");
	agrandissement_bouton.Units = "normalized";
	agrandissement_bouton.Position = [0.8 0.7 0.2 0.1];
	agrandissement_bouton.FontSize = 15;
	agrandissement_bouton.String = "Agrandir";
	agrandissement_bouton.BackgroundColor = Fractale.CouleurFond;
	agrandissement_bouton.Callback = "mandelbrot_trace";
	agrandissement_bouton.Relief = "raised";
	agrandissement_bouton.Enable = "off";
	
//	Commande du tracé de l'ensemble de Julia
	julia_bouton = uicontrol(interface, "style", "pushbutton");
	julia_bouton.Units = "normalized";
	julia_bouton.Position = [0.8 0.5 0.2 0.1];
	julia_bouton.FontSize = 15;
	julia_bouton.String = "Julia";
	julia_bouton.BackgroundColor = Fractale.CouleurFond;
	julia_bouton.Callback = "julia_trace";
	julia_bouton.Relief = "raised";
	julia_bouton.Enable = "off";
	
//	Commande de sauvegarde de la figure
	impression_bouton = uicontrol(interface, "style", "pushbutton");
	impression_bouton.Units = "normalized";
	impression_bouton.Position = [0.8 0.3 0.2 0.1];
	impression_bouton.FontSize = 15;
	impression_bouton.String = "Enregister";
	impression_bouton.BackgroundColor = Fractale.CouleurFond;
	impression_bouton.Callback = "figure_impression";
	impression_bouton.Relief = "raised";

//	Commande de sortie du programme
	sortie_bouton = uicontrol(interface, "style", "pushbutton");
	sortie_bouton.Units = "normalized";
	sortie_bouton.Position = [0.8 0.1 0.2 0.1];
	sortie_bouton.FontSize = 15;
	sortie_bouton.String = "Sortir";
	sortie_bouton.BackgroundColor = Fractale.CouleurFond;
	sortie_bouton.Callback = "xdel(winsid())";
	sortie_bouton.Relief = "raised";
