Index: /issm/trunk-jpl/src/m/materials/nye.m
===================================================================
--- /issm/trunk-jpl/src/m/materials/nye.m	(revision 24071)
+++ /issm/trunk-jpl/src/m/materials/nye.m	(revision 24072)
@@ -1,8 +1,11 @@
 function rigidity = nye(temperature,ice_type)
-	%NYE - figure out the rigidity of ice (either CO2 or H2O) for a given temperature
-	%   rigidity (in s^(1/3)Pa) is the flow law parameter in the flow law sigma=B*e(1/3) (Nye, p2000). 
-	%   temperature is in Kelvin degrees
-	%   Usage:
-	%      rigidity=nye(temperature,ice_type) % ice_type = 1: CO2 ice // ice_type = 2: H2O ice
+%NYE - Nye viscosity coefficient
+%
+%   Compute rigidity of ice (either CO2 or H2O) for a given temperature
+%   rigidity (in s^(1/3)Pa) is the flow law parameter in the flow law
+%   sigma=B*e(1/3) (Nye, p2000).  temperature is in Kelvin degrees
+%
+%   Usage:
+%      rigidity=nye(temperature,ice_type) % ice_type = 1: CO2 ice // ice_type = 2: H2O ice
 
 	if (any(temperature < 130) || any(temperature > 273))
@@ -12,26 +15,21 @@
 	Rg = 8.3144598; % J mol^-1 K^-1
 
-	if ice_type == 1 % CO2 ice
+	if(ice_type==1) %CO2 ice
 		A_const = 10^(10.8); % s^-1 MPa
-		Q = 63000; % J mol^-1
-		n = 7; % Glen's exponent
-		T = 250:-5:100; % K
-		A = A_const*exp(-Q./(T*Rg)); % s^-1 MPa
-		B = A.^(-1/n)*1e6; % s^(1/n) Pa
+		Q       = 63000;     % J mol^-1
+		n       = 7;         % Glen's exponent
+
+	elseif(ice_type==2) % H2O ice
+		A_const = 9e4;   % s^-1 MPa
+		Q       = 60000; % J mol^-1
+		n       = 3;     % Glen's exponent
+
+	else
+		error('ice type not supported');
 	end
 
-	if ice_type == 2 % H2O ice
-		A_const = 9e4; % s^-1 MPa
-		Q = 60000; %  J mol^-1
-		n = 3; % Glen's exponent
-		T = -100:5:0; % Celsius 
-		T = T + 273; % K
-		A = A_const*exp(-Q./(T*Rg)); % s^-1 MPa
-		B = A.^(-1/n)*1e6; % s^(1/n) Pa
-	end
-
-	% Now, do a cubic fit between Temp and B: 
-	fittedmodel = fit(T',B','cubicspline');
-	rigidity = fittedmodel(temperature);
+	%Arhenius law
+	A = A_const*exp(-Q./(temperature*Rg)); % s^-1 MPa
+	B = A.^(-1/n)*1e6; % s^(1/n) Pa
 
 end