/* \file Arrhenius.cpp
 * \brief figure out B of ice for a certain temperature
 */

#include "./elements.h"
#include <math.h>

IssmDouble Arrhenius(IssmDouble temperature,IssmDouble depth,IssmDouble n){
	/*Use EISMINT Parameterization for the rheology: Payne2000
	 *
	 *  A(T*) = A0 exp(-Q/RT*)
	 *
	 *  A0 constant of proportionality
	 *     = 3.61 * 10^-13   if T*<263.15K
	 *     = 1.73 * 10^3     if T*>263.15K
	 *  Q  Activation energy for creep
	 *     = 6.0  * 10^4     if T*<263.15K
	 *     = 13.9 * 10^4     if T*>263.15K
	 *  R  Universal gas constant
	 *     = 8.314
	 *  T* Absolute temperature corrected for the dependence of Tpmp on P
	 *     = T - beta (s-z)
	 *
	 *  Convert A to B :  B = A^(-1/n) */

	/*Some physical constants (Payne2000)*/
	IssmDouble beta=8.66*pow(10.,-4.);
	IssmDouble R=8.314;

	/*Intermediaries*/
	IssmDouble A,B,Tstar;

	/*convert temperature to absolute temperature*/
	_assert_(depth);
	Tstar=temperature-beta*depth;
	_assert_(Tstar);

	/*Get A*/
	if(Tstar<263.15){
		A=3.61*pow(10.,-13.) * exp(  -6.*pow(10.,4.)/(R*Tstar));
	}
	else{
		A=1.73*pow(10.,  3.) * exp(-13.9*pow(10.,4.)/(R*Tstar));
	}

	/*Convert to B*/
	B=pow(A,-1./n);

	return B;
}
