% td4_ul3- Calcul de l'aire délimitée par une courbe
%          et calcul d'un volume d'une ellipsoide
%          par la méthode de Monte Carlo.
%          Licence de Physique, 1997-1998
%          Université Paul Sabatier, JLA
%	   Version 3.0

help td4_ul3; clear all
hold off
clf
clg

% Initialisation
% --------------

	initiales='jluc'
	n_total=1000

% menu pour les differents cas
% ----------------------------

k=menu('Choix de calcul', 'Aire d"Ellipse', 'Aire de gaussienne', ...
	'Volume d"ellipsoide', 'fin');



if k==1
% a)
% calcul de l'aire d'un  cercle  (à compléter)
% ou d'une d'ellipse. A la recherche de pi.......
% ----------------------------------------------------

          n_interne=0;

          a=?;
	  b=?;
hold on
          for i=1:n_total

        	x(i)=rand*a;
	        y(i)=rand*b;
		plot(x(i),y(i));

	    if ????????????????
                n_interne=n_interne+1;
		plot(x(i),y(i), 'c');               
            end
    	     aire(i)=4*n_interne/i;
	  end
	    aire_totale=aire(n_total)              % Aire complete

% Graphes
% -------

   z=0:0.01:a;
   y=b*sqrt(1-z.^2/a^2)
   plot(z,y)
   title([initiales, ' Aire calculee=', num2str(aire_totale), ' Aire vraie=', ...
   num2str(pi*a*b),' a,b=',', ',num2str(a),', ',num2str(b)])
   xlabel('Nombre de tirages')
   ylabel('Valeur de l"aire')
   hold off

   pause


   plot(aire)
   title([initiales, ' Ellipse entiere: Calcul=', num2str(aire_totale),...
    ' a,b=',', ',num2str(a),', ',num2str(b)])
   xlabel('Nombre de tirages')
   xlabel('Nombre de tirages')
   ylabel('Valeur de l"aire')

        
% Calcul de la moyenne
% --------------------

	abs_nouv=n_total/10  
        moyenne=4*sum(aire(abs_nouv:n_total))/(n_total-abs_nouv)*a*b
end

% --------------------------------------------------------------------

if k==2
% b)
% Calcul de l'aire de la gaussienne normalisée
% entre -sigma et +sigma puis -2s et +2s
% --------------------------------------------
% initialisation (a completer)
% ----------------
	sigma=?;

       	valeur=2*sigma/n_total;         % Nombre de points pour la gaussienne

% Définition de la 1/2 gaussienne
% ---------------------------
	hold on
	A=1/(sqrt(2*pi)*sigma);               % Amplitude de la gaussienne

	ik=0;
	for i=0:valeur:sigma

	 ik=ik+1;
         gauss(ik)=A*exp(-i^2/(2*sigma^2));
         axe_x(ik)=i;

	end

          mini=A*exp(-0.5);                % mini pour x=sigma
          plot(axe_x, gauss)
          n_interne=0;

          for i=1:n_total

 		x=?;            % tirage selon x
                y=?;		% tirage selon y
                plot(x,y)

            if ??????                              % interieure de la gaussienne
            plot(x,y, 'c')            
              n_interne=n_interne+1;
            end

    	     aire(i)=n_interne/i;
	  end

         aire_exacte=?                                           ;% Aire de la gaussienne
         aire_totale=aire(n_total)                         % valeur de l'aire complete

% Graphes
% -------

        plot(axe_x, gauss)
	title([initiales, ' Aire Calculee=', num2str(aire_totale/2)])
	xlabel('x')
	ylabel('y')
        hold off
        pause

%----------- Gaussienne entiere

	extreme=2*sigma
	ik=0;
	for i=0:valeur:extreme
	 ik=ik+1;
         gauss(ik)=A*exp(-i^2/(2*sigma^2));
         axe_x(ik)=i;
	end
         ik2=ik
	for i=-extreme:valeur:0
	 ik=ik+1;
         gauss(ik)=A*exp(-i^2/(2*sigma^2));
         axe_x(ik)=i;
	end


	subplot(211)
	hold on
	plot(axe_x, gauss, [-sigma, -sigma], [gauss(ik2), gauss(1)],[sigma, sigma],[gauss(ik2), gauss(1)])
        hold off
	subplot(212)
	plot(aire)
	title([initiales, ' Calcul=', num2str(aire_totale)])
	xlabel('Nombre de tirages')
	ylabel('Aire de la gaussienne')

 % moyenne (absisse nouvelle pour moyenne)
 % ---------------------------------------

	abs_nouv=n_total/10
        moyenne=2*sum(aire(abs_nouv:n_total))/(n_total-abs_nouv)
	
end

if k==3
% d)
% calcul du volume d'une sphere de rayon 1
% calcul du volume d'une ellipsoide
% ---------------------------------
clg

% Parametre de l'ellipsoide (à completer)
% --------------------------
            a=?;
	    b=?;
            c=?;

            n_interne=0;
          for i=1:n_total
		
		x(i)=?;                % tirage selon x
		y(i)=?;		       % tirage selon y
		z(i)=?;		       % tirage selon z
                
            if ???????????
              n_interne=n_interne+1;
            end

    	     volume(i)=n_interne/i;
	  end

        volume=?;                  % volume de l'ellipsoide
	volume_total=volume(n_total)            % valeur finale du volume
	
% Graphes
% -------

	plot(volume)
	title([initiales, ' Calcul=', num2str(volume_total), '/ a,b,c=', num2str(a),...
             ', ', num2str(b),', ', num2str(c)])
	xlabel('Nombre de tirages')
	ylabel('Volume ellipsoide')
        pause


 % moyenne
 % -------
	abs_nouv=n_total/10
        moyenne=sum(volume(abs_nouv:n_total))/(n_total-abs_nouv)


% Etude de la convergence
% -----------------------

	for k=1:n_total-100
         ecart_type(k)=std(volume(k:k+100));
	end


	plot(ecart_type)
        title([initiales, 'Ecart type ellipsoide', '/ a,b,c=', num2str(a),...
             ', ', num2str(b),', ', num2str(c)])
	xlabel('Nombre de tirages')
	ylabel('Ecart_Type')
        hold off
	
end
end




























