/* \file Cuffey.cpp
 * \brief figure out B of ice for a certain temperature
 *	  INPUT function B=Cuffey(temperature)
 *    where rigidigty (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Cuffey, p75). 
 */

#include <math.h>

#include "../Numerics/types.h"

IssmDouble Cuffey(IssmDouble temperature){

	/*output: */
	IssmDouble B,T;

	/*Switch to celsius from Kelvin: */
	T=temperature-273.15;

	if(T<=-45.0){
		B=1.e+8*(-0.000396645116301*pow(T+50.,3)+ 0.013345579471334*pow(T+50.,2) -0.356868703259105*(T+50.)+7.272363035371383);
	}
	else if((T>=-45.0) && (T<=-40.0)){
		B=1.e+8*(-0.000396645116301*pow(T+45.,3)+ 0.007395902726819*pow(T+45.,2) -0.253161292268336*(T+45.)+5.772078366321591);
	}
	else if((T>=-40.0) && (T<=-35.0)){
		B=1.e+8*(0.000408322072669*pow(T+40.,3)+  0.001446225982305*pow(T+40.,2)  -0.208950648722716*(T+40.)+4.641588833612773);
	}
	else if((T>=-35.0) && (T<=-30.0)){
		B=1.e+8*(-0.000423888728124*pow(T+35.,3)+ 0.007571057072334*pow(T+35.,2)  -0.163864233449525*(T+35.)+3.684031498640382);
	}
	else if((T>=-30.0) && (T<=-25.0)){
		B=1.e+8*(0.000147154327025*pow(T+30.,3)+ 0.001212726150476*pow(T+30.,2) -0.119945317335478*(T+30.)+3.001000667185614);
	}
	else if((T>=-25.0) && (T<=-20.0)){
		B=1.e+8*(-0.000193435838672*pow(T+25.,3)+ 0.003420041055847*pow(T+25.,2)  -0.096781481303861*(T+25.)+2.449986525148220);
	}
	else if((T>=-20.0) && (T<=-15.0)){
		B=1.e+8*(0.000219771255067*pow(T+20.,3)+  0.000518503475772*pow(T+20.,2)  -0.077088758645767*(T+20.)+2.027400665191131);
	}
	else if((T>=-15.0) && (T<=-10.0)){
		B=1.e+8*(-0.000653438900191*pow(T+15.,3)+ 0.003815072301777*pow(T+15.,2)  -0.055420879758021*(T+15.)+1.682390865739973);
	}
	else if((T>=-10.0) && (T<=-5.0)){
		B=1.e+8*(0.000692439419762*pow(T+10.,3) -0.005986511201093 *pow(T+10.,2) -0.066278074254598*(T+10.)+1.418983411970382);
	}
	else if((T>=-5.0) && (T<=-2.0)){
		B=1.e+8*(-0.000132282004110*pow(T+5.,3) +0.004400080095332*pow(T+5.,2)    -0.074210229783403*(T+5.)+ 1.024485188140279);
	}
	else{
		B=1.e+8*(-0.000132282004110*pow(T+2.,3) +0.003209542058346*pow(T+2.,2)    -0.051381363322371*(T+2.)+ 0.837883605537096);
	}

	/*B cannot be negative!*/
	if(B<0) B=1.e+6;

	return B;
}