Index: ../trunk-jpl/src/m/materials/nye.py =================================================================== --- ../trunk-jpl/src/m/materials/nye.py (revision 24793) +++ ../trunk-jpl/src/m/materials/nye.py (revision 24794) @@ -1,51 +1,53 @@ import numpy as np +from warnings import warn + def nye(temperature, ice_type): - """ - NYE - figure out the rigidity of ice (either CO2 or H2O) for a given temperature - rigidity (in s^(1/n)Pa) is the flow law parameter in the flow law sigma=B*e(1/n) (Nye, p2000). - temperature is in Kelvin degrees + """ + NYE - figure out the rigidity of ice (either CO2 or H2O) for a given temperature + rigidity (in s^(1/n)Pa) is the flow law parameter in the flow law sigma=B*e(1/n) (Nye, p2000). + temperature is in Kelvin degrees Usage: - rigidity=nye(temperature,ice_type) % ice_type = 1: CO2 ice // ice_type = 2: H2O ice - """ + rigidity=nye(temperature,ice_type) % ice_type = 1: CO2 ice // ice_type = 2: H2O ice + """ - # Declaring temperature and rigidity arrays - if np.ndim(temperature)==2: - T=temperature.flatten() - elif isinstance(temperature,float) or isinstance(temperature,int): - T=np.array([temperature]) - else: - T=temperature - rigidity=np.zeros_like(T) + # Declaring temperature and rigidity arrays + if np.ndim(temperature) == 2: + T = temperature.flatten() + elif isinstance(temperature, float) or isinstance(temperature, int): + T = np.array([temperature]) + else: + T = temperature + rigidity = np.zeros_like(T) - # Beyond-melting-point cases - if (ice_type==1): - for i in range(len(T)): - if (200=220).any()): - warnings.warn('CO2 ICE - GUARANTEED MELTING. Some temperature values are beyond 220K.') - elif (ice_type==2) and ((T>273.15).any()): - warnings.warn('H2O ICE - GUARANTEED MELTING. Some temperature values are beyond 273.15K.') + # Beyond-melting-point cases + if (ice_type == 1): + for i in range(len(T)): + if (200 < T[i] < 220): + warn('CO2 ICE - POSSIBLE MELTING. Some temperature values are between 200K and 220K.') + break + if ((T >= 220).any()): + warn('CO2 ICE - GUARANTEED MELTING. Some temperature values are beyond 220K.') + elif (ice_type == 2) and ((T > 273.15).any()): + warn('H2O ICE - GUARANTEED MELTING. Some temperature values are beyond 273.15K.') - Rg = 8.3144598 # J mol^-1 K^-1 + Rg = 8.3144598 # J mol^-1 K^-1 - if ice_type == 1: # CO2 ice - A_const = 10**(13.0) # s^-1 MPa - Q = 66900. # J mol^-1 - n = 8. # Glen's exponent - elif ice_type == 2: # H2O ice - A_const = 9*10**4 # s^-1 MPa - Q = 60000. # J mol^-1 - n = 3. # Glen's exponent - else: - raise RuntimeError('Ice type not supported') + if ice_type == 1: # CO2 ice + A_const = 1.0e13 # s^-1 MPa + Q = 66900. # J mol^-1 + n = 8. # Glen's exponent + elif ice_type == 2: # H2O ice + A_const = 9 * 1.0e4 # s^-1 MPa + Q = 60000. # J mol^-1 + n = 3. # Glen's exponent + else: + raise RuntimeError('Ice type not supported') - # Arrhenius Law - A=A_const*np.exp(-1*Q/(T*Rg)) # s^-1 MPa - rigidity=A**(-1/n)*10**6 # s^(1/n) Pa + # Arrhenius Law + A = A_const * np.exp(-1 * Q / (T * Rg)) # s^-1 MPa + rigidity = A**(-1 / n) * 1.0e6 # s^(1/n) Pa - # Return output - return rigidity + # Return output + return rigidity