/* \file LliboutryDuval.cpp
 * \brief figure out B of ice for a certain temperature and water fraction or enthalpy
 */

#include <math.h>
#include "../Numerics/types.h"
#include "../Exceptions/exceptions.h"


/* get ice stiffness B from enthalpy, pressure and flow law exponent*/
IssmDouble LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure,IssmDouble n){
  /*Use parameterization for the rheology: Aschwanden 2012
   *
   *  A(H,p) = A0 exp(-Q/RT(H,p)), if H < H_s(p)
   *         = A0 exp(-Q/RTpmp) (1+181.25w(H,p)), if H_s(p) \le H < H_l(p)
   *  
   *  T(H,p) = Tref + H/c_i, if H < H_s(p)
   *         = Tpmp , if H_s(p) \le H \le H_l(p)
   *
   *  w(H,p) = 0, if H < H_s(p)
   *         = (H - H_s(p))/L
   *
   *  H_s(p) = c_i (Tpmp - Tref)
   *
   *  Tpmp   = T - beta p;
   *
   *  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
   *  
   *  Convert A to B :  B = A^(-1/n) */
  
  /*Some physical constants (Aschwanden 2012)*/
  /*TODO: get those constants from model*/
  IssmDouble beta=7.9*pow(10.,-8.);
  IssmDouble R=8.314;
  IssmDouble heatcapacity=2009; // J/kg/K
  IssmDouble Tref=253.15;
  IssmDouble latentheat=3.34*pow(10,5.); // from Aschwanden 2012

  /*Intermediaries*/
  IssmDouble A,B,Tstar,Tpmp,H_sp,waterfraction;
	
  _assert_(pressure>0);
  _assert_(enthalpy>0);
  Tpmp=273.15-beta*pressure; 
  H_sp=heatcapacity*(Tpmp - Tref);
  if (enthalpy < H_sp){
    Tstar = Tref + enthalpy/heatcapacity - beta*pressure;	
    waterfraction = 0;
  }
  else{
    Tstar=Tpmp;
    waterfraction=(enthalpy - H_sp)/latentheat;
  }

  /*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));
  }
  A*=(1 + 181.25*waterfraction);

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

  return B;
}

/*Get stiffness from temperature, waterfraction and depth*/
IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n){
	/*Use parameterization for the rheology: Greve and Blatter 2009
	 * get enthalpy from temperature and water fraction,
	 * and use LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure,IssmDouble n) */
  
  IssmDouble rho_ice=910; // kg/m^3
  IssmDouble g=9.81; //kg*m/s^2
  IssmDouble heatcapacity=2009; // J/kg/K
  IssmDouble Tref=253.15;
  IssmDouble beta=7.9*pow(10.,-8.);
  IssmDouble latentheat=3.34*pow(10,5.); // from Aschwanden 2012
  IssmDouble Tstar, enthalpy, pressure, B;
  _assert_(temperature>0);
  _assert_(waterfraction>0);
  _assert_(depth>0);

  /*get pressure*/
  pressure= rho_ice*g*depth;
  Tstar=temperature-beta*pressure; // TODO: check whether plus or minus
  /*get enthalpy*/
  if (Tstar < 273.15){
    enthalpy=heatcapacity*(Tstar - Tref);
  }
  else{
    enthalpy=heatcapacity*(273.15 - Tref) + waterfraction*latentheat;
  }
	
  B=LliboutryDuval(enthalpy, pressure, n);

  return B;
}

