%Función principal del programa
function varargout = CalculoArea(varargin)
	%Comprobando como se realizó la llamada de la función
	switch nargin
		case 0
			[pRec] = EntradaPuntos();
		case 1
			pRec = varargin{1};
		case 2
			pRec = [varargin{1},varargin{2}];
	end
	%Comprobando que se pueda trazar la figura
	if length(pRec(:,1))<2
		fprintf('Debe introducir al menos 2 puntos.\n');
		return
	end
	%Calculando el centroide de los puntos
	centroide = CalculoCentroide(pRec);
	%Transformando los puntos a coordenadas polares
	[pPol] = RectangularAPolar(pRec,centroide);
	%Calculando el area de la figura
	I = AreaDPolar(pPol);
	%Formando la salida
	if nargout==0
		fprintf('El área es: %f\n',I);
	else
		varargout{1}=I;
	end
end

%Función que pide los datos al usuario
function [pRectangular] = EntradaPuntos(void)
	k=1;
	fin = 'fin';
	fprintf('Introduce los puntos conocidos:\n');
	while (1)
		fprintf('Punto %d.\n',k);
			xE = input('	x: ');
			if(xE == fin)
				break;
			end
			yE = input('	y: ');
			if(yE == fin)
				break;
			end
		if(k==1)
			x = xE;
			y = yE;
		else
			x = [x(:,:);xE];
			y = [y(:,:);yE];
		end
		k=k+1;
	end
	pRectangular = [x,y];
end

%Calcula el centroide de un conjunto de puntos
function [centroide] = CalculoCentroide(puntos)
	numPuntos = length(puntos);
	centroide = (ones(1,numPuntos)*puntos(:,1))/(numPuntos);
	centroide = [centroide,(ones(1,numPuntos)*puntos(:,2))/(numPuntos)];
end

%Transforma una serie de puntos de rectangular a polar, formando una figura cerrada
function [pPolar] = RectangularAPolar(pRectangular,centroide)
	radio = sqrt((pRectangular(:,1)-centroide(1)).^2 + (pRectangular(:,2)-centroide(2)).^2);
	angulo = atan2(pRectangular(:,2)-centroide(2),pRectangular(:,1)-centroide(1));
	for k=1:length(angulo)
		if(angulo(k)<0)
			angulo(k) = angulo(k)+(2*pi);
		end
	end
	pPolar = [angulo,radio];
	pPolar = sortrows(pPolar,1);
	pPolar = [pPolar;[pPolar(1,1)+(2*pi),(pPolar(1,2))]];
end

%Regla de Simpson 3/8 para calcular el area
function I = Simpson38(x, y, a, b, N) 
	h = abs(b-a)/N;
	suma =0;
	for k=0:N-1
    	x0 = a + h*k;
    	x1 = x0 + h/3;
    	x2 = x0 + 2*h/3;
    	x3 = x0 + h;
		suma = suma + f(x0,x,y) + 3*f(x1,x,y) + 3*f(x2,x,y) + f(x3,x,y);
	end;
	I = h*suma/8;
end

%Esta función es la interpolación lineal
function y = f(x0, x, y)
	N =length(x);
	for k=1:N-1
    	if x0 >= x(k) && x0 <= x(k+1) break;
    	end
	end
	y = y(k) + (y(k+1) - y(k))/(x(k+1)-x(k))*(x0-x(k));
end

%Regresa el area bajo una curva en el dominio polar
function [area] = AreaDPolar(p)
	area = Simpson38(p(:,1), p(:,2).^2, p(1,1), p(length(p),1) ,length(p))/2;
end