clear;
stacksize(150000000);
driver("X11");
clf(1);
scf(1);
//xset("wpos",600,0)
xset("wpos",300,0)

f=gcf();
f.pixmap="on";
//xset("wdim",1000,1000);
xset("wdim",620,630);
rand("uniform");

//clf(2);
//scf(2);
////xset("wpos",600,0)
//xset("wpos",300,0)
//xselect();
//f2=gcf();
//f2.pixmap="on";
////xset("wdim",1000,1000);
//xset("wdim",640,630);

%symbol=1; 
%times=2; 
%times_italic=3; 
%times_bold=4; 
%helvetica= 6; 
%helvetica_italic=7; 
%helvetica_bold=8; 
%helvetica_bold_italique=9; 
%8pts=0;%10pts=1;%12pts=2; 
%14pts=3;%18pts=4;%24pts=5; 


a=1.7;
b=0.01;
deff('y=fct1(x)','y=x.*(3-x).*(3-x).*((x-a).*(x-a)+b).*( ( ones(x)./( (x-2.6).*(x-2.6) +1 ) ) +ones(x) ).*( ( ones(x)./( (x-0.9).*(x-0.9) +1 ) ) +ones(x) )./(3.1-x)^2');//
C=1./intg(0,3,fct1);
deff('y=fct_densite(x)','y=C*fct1(x)');

Nabs=500;
x_abs=linspace(0,3,Nabs+1);
y_dens=fct_densite(x_abs);
y_rep=[0];
Netapes=5;
ncolhisto=20;
nbre_histo_fluct=10;

for i=1:Nabs
	y_rep=[y_rep,intg(0,3*i/Nabs,fct_densite)];
end


function x=inv_rep(y)
	z=abs(y_rep-y);
	[m,k]=min(z);
	x=3*k/Nabs;
endfunction

ymax=1.05*max(y_dens);
rect_bas=[0,0,3,ymax];
rect_haut=[0,0,3,1];
xmin=.1;


function affichage_haut() 
	plot2d(x_abs,y_rep,style=5,rect=rect_haut,frameflag=1);
	a = gca(); // a est le handle sur ces axes par defaut 
	a.title.font_size = %24pts; //(=%18pts;)
	a.title.font_style = 6;  
	a.title.text="Fonction de repartition";
	//a.thickness=2; 
	xselect();
	show_pixmap();
endfunction 

function affichage_bas() 
	plot2d(x_abs,y_dens,style=5,rect=rect_bas,frameflag=1);
	a = gca(); // a est le handle sur ces axes par defaut 
	a.title.font_size = %24pts; //(=%18pts;)
	a.title.font_style = 6;  
	a.title.text="Densite";
	//a.thickness=2;
	xselect();
	show_pixmap();
endfunction 



subplot(2,1,1);

affichage_haut();

subplot(2,1,2);

affichage_bas();

subplot(2,1,1);

vect_choix=[1,100,1000];
vect_choix_string=[];
for i=1:length(vect_choix)
	vect_choix_string=[vect_choix_string;"+"+string(vect_choix(i))];
end

vect_choix_string=[vect_choix_string;"Histogramme"];

Nbre_sim=0;
ech=[];
last_to_erase=%F;
last_to_erasebas=%F;

factice=0;

while %T
	Choix=x_choose(vect_choix_string,string(Nbre_sim)+" simulations","Quitter");
	xselect();
	
	select Choix
	
	case 0 then break;
	
	case 1 then
		if last_to_erase then delete(p); end
		if factice==0 then u=0.7+0.05*rand(1,1); factice=1; //petite tricherie pour les besoins de la clarté de l'expose en amphi.
		elseif factice==1 then u=0.25+0.05*rand(1,1); factice=2;
		else u=rand(1,1); end
		xpoly(0.001,u);
		qa=get("hdl");
		qa.mark_foreground=5;
		qa.mark_background=5;
		qa.mark_style=8;
		xselect();
		show_pixmap();
		xclick();
		sim=inv_rep(u);
		abs_p=[linspace(0,sim,Netapes)';sim*ones(Netapes,1)];
		ord_p=[u*ones(Netapes,1);linspace(u,0,Netapes)'];
		xpoly(0,u);
		p=get("hdl");
		p.foreground=2;
		xselect();
		show_pixmap;
		inv_rep(u);
		for i=2:2*Netapes
			p.data=[abs_p(1:i),ord_p(1:i)];
			xselect();
			show_pixmap();
		end //i=2:2*Netapes
		xpoly(sim,0.001);
		q=get("hdl");
		q.mark_foreground=13;
		q.mark_background=13;
		q.mark_style=8;
		xselect();
		show_pixmap();
		Nbre_sim=Nbre_sim+1;
		ech=[ech,sim];
		last_to_erase=%T;
	
	case length(vect_choix)+1 then
		subplot(2,1,2);
		a=gca();
		a.title.text="Densite et histogramme des valeurs obtenues";
		if last_to_erasebas then delete(a.children(1)); end
		
		for j=1:nbre_histo_fluct
			
			if j>1 then delete(a.children(1)); end
			
			n=max(round(Nbre_sim*j/nbre_histo_fluct),2);
			histplot(ncolhisto,ech(1:n),rect=rect_bas, frameflag=1);
			xselect();
			show_pixmap();
			
		end // pour for j=1:nbre_histo_fluct
		last_to_erasebas=%T;
		subplot(2,1,1);
		
	case 2
		U=rand(1,vect_choix(Choix));
		for k=1:vect_choix(Choix)
			if last_to_erase then delete(p); end
			u=U(k);
			xpoly(0.001,u);
			qa=get("hdl");
			qa.mark_foreground=5;
			qa.mark_background=5;
			qa.mark_style=8;
			sim=inv_rep(u);
			abs_p=[0;sim;sim];
			ord_p=[u;u;0];
			xpoly(abs_p,ord_p);
			p=get("hdl");
			p.foreground=2;
			xpoly(sim,0.001);
			q=get("hdl");
			q.mark_foreground=13;
			q.mark_background=13;
			q.mark_style=8;
			xselect();
			show_pixmap();
			ech=[ech,sim];
			last_to_erase=%T;
		end // pour for k=1:vect_choix(Choix)
		Nbre_sim=Nbre_sim+vect_choix(Choix);
	
	else
		U=rand(1,vect_choix(Choix));
		for k=1:vect_choix(Choix)/50
			if last_to_erase then delete(p); end
			u=U((k-1)*50+1);
			//plot(0.01,u,"r.");
			sim=inv_rep(u);
			abs_p=[0;sim;sim];
			ord_p=[u;u;0];
			xpoly(abs_p,ord_p);
			p=get("hdl");
			p.foreground=2;
			//xpoly(sim,0.001);
			//q=get("hdl");
			//q.mark_foreground=13;
			//q.mark_style=8;
			xselect();
			show_pixmap();
			ech=[ech,sim];
			last_to_erase=%T;
			
			for j=2:50
				u=U((k-1)*50+j);
				sim=inv_rep(u);
				ech=[ech,sim];
			end // pour for j=2:50
			xpoly(0.001*ones(50,1),U((k-1)*50+1:k*50)');
			qa=get("hdl");
			qa.mark_foreground=5;
			qa.mark_background=5;
			qa.mark_style=8;
			qa.line_mode="off";
			xpoly(ech($-49:$)',0.001*ones(50,1));
			qa=get("hdl");
			qa.mark_foreground=13;
			qa.mark_background=13;
			qa.mark_style=8;
			qa.line_mode="off";
		end // pour for k=1:vect_choix(Choix)
		Nbre_sim=Nbre_sim+vect_choix(Choix);

	
	end //pour select Choix
end // pour le while %T
f.pixmap="off";
driver("Rec");