function mu2=GetViscosity2(matice,epsilon)
%GETVISCOSITY2 - returns viscosity for control method
%
%   Return viscosity in 2d, for triaelem, used in the inverse control methods.
%   Careful, this is not the visocsity used by the diagnostic equations!
%	                                  2*(1-n)/2n
%	 mu2= -------------------------------------------------------------------
%	      2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(3n-1)/2n]
%	 
%	 where mu2 is the viscotiy, (u,v) the velocity 
%	 vector, and n the flow law exponent.
%	 
%	 If epsilon is NULL, it means this is the first time Emg is being run, and we 
%	 return 10^14, initial viscosity.
%	
%   usage:
%      mu2=getviscosity2(matice,epsilon)

eps0=10^-27;
if matice.n==1,
	error('control method not implemented yet for viscous case');
else
	%Non linear creep
	if norm(epsilon,2)~=0, 

		n=matice.n;

		exx=epsilon(1,:);
		eyy=epsilon(2,:);
		exy=epsilon(3,:);

		%Build viscosity: mu2=2*B/(2*A^e):
		A=exx.^2+eyy.^2+exy.^2+exx.*eyy+eps0^2;
		e=(3*n-1)/2/n;
		mu2=(2*(1-n)/2/n*ones(1,size(epsilon,2)))./(2*A.^e);
	else
		mu2=10^14*size(epsilon,2);
	end
end

