// Problème 3

// Définition de fonctions (problème 1)

// Cône de hauteur 4 de rayon de base 1 (x entre 0 et 4)
// Commande pour l'intégrale : intMoy(pb1a,N,0,4)
// Volume théorique : (1/3) * pi * 1 * 4 = 4/3 * pi
function y = pb1a(x)
  y = 1-(x/4);
  y = %pi*y*y;
endfunction

// Sphère de rayon 2 (x entre -2 et 2)
// Commande pour l'intégrale : intMoy(pb1b,N,-2,2)
// Volume théorique : (4/3 * pi * 2^3) = (32/3) pi
function y = pb1b(x)
  y = sqrt(4 - x*x);
  y = %pi*y*y;
endfunction

// Paraboloïde (x entre 2 et 7)
// Commande pour l'intégrale : intMoy(pb1c,N,2,7)
function y = pb1c(x)
  y = sqrt(2*x - 4);
  y = %pi*y*y;
endfunction

// Définition de domaines (problème 2)
// Ces fonctions renvoient true si le point donné correspond
// au domaine décrit, false sinon.

// Tore de section carrée (x entre 1 et 2, y entre 2 et 3)
// Commande pour l'intégrale : intRA(pb2a,9,N,1,2)
function res = pb2a(x,y)
  res = (y >= 2*2) & (y <= 3*3);
endfunction

// Solide entre x² et x^1/3 (x entre 0 et 1)
// Commande pour l'intégrale : intRA(pb2b,1,N,0,1)
function res = pb2b(x,y)
  res = (y >= x*x*x*x) & (y <= x^(2/3));
endfunction

// Méthodes d'intégration

// Méthode de la moyenne
function res = intMoy(f,n,a,b)
  s = 0;
  rand('seed',getdate('s')); // Génération d'un nouveau seed
  for i = 1:n
    r = (b-a) * rand();
    s = s + f(a + r);
  end
  res = (b-a)*(s/n);
endfunction

// Méthode rejet-acceptation
function res = intRA(f,m,n,a,b)
  cpt = 0;
  rand('seed',getdate('s'));
  for i = 1:n
    y = m*rand();
    x = a + (b-a) * rand();
    if f(x,y)
      cpt = cpt+1;
    end
  end
  res = %pi*m*(b-a)*(cpt/n);
endfunction

// Problème 4

// Définition de fonctions (2 variables)
// x entre 0 et 1, y entre 0 et 1
// Commande pour l'intégrale : intMoyXY(pb4,N,0,1,0,1)
function res = pb4(x,y)
  res = (1 - y) * (2 + sin(5*%pi*x));
endfunction

// Méthode moyenne (2 variables)
function res = intMoyXY(f,n,a,b,c,d)
  s = 0;
  rand('seed',getdate('s')); // Génération d'un nouveau seed
  for i = 1:n
    s = s + f(a + (b-a) * rand(), c + (d-c) * rand());
  end
  res = (b-a)*(d-c)*(s/n);
endfunction

// Problème 5

// Définition de volumes décrits
// Renvoie true si un point est compris dans le volume
// false sinon

// Sphère centrée en 0
// Commande pour l'intégrale : intRAxyz(pb5a,N,-r,r,-r,r,-r,r)
// Pour r = 3 : intRAxyz(pb5a,N,-3,3,-3,3,-3,3)
// (r le rayon voulu)
function res = sphere(x,y,z,r)
  res = (x*x + y*y + z*z - r*r <= 0);
endfunction

function res = pb5a(x,y,z)
  res = sphere(x,y,z,3);
endfunction

// Intersection cylindre + sphère
// z entre -2 et 4, x entre 0 et 4, y entre -2 et 2
// Commande pour l'intégrale : intRAxyz(pb5b,N,0,4,-2,2,-2,4)
function res = pb5b(x,y,z)
  res = (x-2)*(x-2) + y*y - 4 <= 0; // Cercle de centre (x,y) = (2,0), rayon 2
  res = res & sphere(x,y,z,3); // Intersection (rayon 3)
endfunction
  

// Méthode rejet-acceptation (problème 5) 
function res = intRAxyz(f,n,x1,x2,y1,y2,z1,z2)
  cpt = 0;
  rand('seed',getdate('s'));
  for i = 1:n
    x = x1 + (x2-x1) * rand();
    y = y1 + (y2-y1) * rand();
    z = z1 + (z2-z1) * rand();
    if f(x,y,z)
      cpt = cpt+1;
    end
  end
  res = (x2-x1)*(y2-y1)*(z2-z1)*(cpt/n);
endfunction