%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of NS2DDV.                                                      %
%                                                                                   %
% Copyright(C) 2011-2018    C. Calgaro  (caterina.calgaro@math.univ-lille1.fr)      %
%                           E. Creusé   (emmanuel.creuse@math.univ-lille1.fr)       %
%                           T. Goudon   (thierry.goudon@inria.fr)                   %
%                           A. Mouton   (alexandre.mouton@math.univ-lille1.fr)      %
%                                                                                   %
% NS2DDV is free software: you can redistribute it and/or modify it under the terms %
% of the GNU General Public License as published by the Free Software Foundation,   %
% either version 3 of the License, or (at your option) any later version.           %
%                                                                                   %
% NS2DDV is distributed in the hope that it will be useful, but WITHOUT ANY         %
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A   %
% PARTICULAR PURPOSE. See the GNU General Public License for more details.          %
%                                                                                   %
% You should have received a copy of the GNU General Public License along with      %
% NS2DDV. If not, see <http://www.gnu.org/licenses/>.                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[res] = gamma_ratio(MESH, PARAMETERS)

	% INPUT
	% MESH			Mesh metadata
	% PARAMETERS	List of simulation parameters
	%
	% OUTPUT
	% res			The inverse of the maximal value of tau for defining tau-limiters

	gammaloc = zeros(MESH.nt1,1);

	for ke=1:MESH.nt1
		A = MESH.p1(:,MESH.t1(1,ke));
		B = MESH.p1(:,MESH.t1(2,ke));
		C = MESH.p1(:,MESH.t1(3,ke));

		if (strcmp(PARAMETERS.DOMAIN.GEOMETRY, 'RECTANGLE') && ...
				strcmp(PARAMETERS.MESH.GENERATION, 'NS2DDV') && strcmp(PARAMETERS.FV.CELLS_DESIGN, 'SQUARES'))
			% The cell centers are the circumcenters
			D = 2.*(A(1)*(B(2)-C(2)) + B(1)*(C(2)-A(2)) + C(1)*(A(2)-B(2)));
			G = zeros(2,1);
			G(1) = ( (A(1)^2+A(2)^2)*(B(2)-C(2)) + (B(1)^2+B(2)^2)*(C(2)-A(2)) + ...
							(C(1)^2+C(2)^2)*(A(2)-B(2)) ) / D;
			G(2) = ( (A(1)^2+A(2)^2)*(C(1)-B(1)) + (B(1)^2+B(2)^2)*(A(1)-C(1)) + ...
							(C(1)^2+C(2)^2)*(B(1)-A(1)) ) / D;
		else
			% The cell centers are the centroids
			G = (A+B+C)/3.;
		end

		Ap = 0.5*(B+C);
		Bp = 0.5*(C+A);
		Cp = 0.5*(A+B);

		App = 0.5*(Ap+G);
		Bpp = 0.5*(Bp+G);
		Cpp = 0.5*(Cp+G);

		A1ppp = intersection(A, Bpp, B, C);
		A2ppp = intersection(A, Cpp, B, C);
		B1ppp = intersection(B, Cpp, C, A);
		B2ppp = intersection(B, App, C, A);
		C1ppp = intersection(C, App, A, B);
		C2ppp = intersection(C, Bpp, A, B);

		ABpp = sqrt((Bpp(1)-A(1))^2+(Bpp(2)-A(2))^2);
		ACpp = sqrt((Cpp(1)-A(1))^2+(Cpp(2)-A(2))^2);
		BCpp = sqrt((Cpp(1)-B(1))^2+(Cpp(2)-B(2))^2);
		BApp = sqrt((App(1)-B(1))^2+(App(2)-B(2))^2);
		CApp = sqrt((App(1)-C(1))^2+(App(2)-C(2))^2);
		CBpp = sqrt((Bpp(1)-C(1))^2+(Bpp(2)-C(2))^2);
		AA1ppp = sqrt((A1ppp(1)-A(1))^2+(A1ppp(2)-A(2))^2);
		AA2ppp = sqrt((A2ppp(1)-A(1))^2+(A2ppp(2)-A(2))^2);
		BB1ppp = sqrt((B1ppp(1)-B(1))^2+(B1ppp(2)-B(2))^2);
		BB2ppp = sqrt((B2ppp(1)-B(1))^2+(B2ppp(2)-B(2))^2);
		CC1ppp = sqrt((C1ppp(1)-C(1))^2+(C1ppp(2)-C(2))^2);
		CC2ppp = sqrt((C2ppp(1)-C(1))^2+(C2ppp(2)-C(2))^2);

		gammaloc(ke) = max([ABpp/AA1ppp, ACpp/AA2ppp, BCpp/BB1ppp, BApp/BB2ppp, CApp/CC1ppp, CBpp/CC2ppp]);
	end

	res = max(gammaloc);

end