source: issm/branches/trunk-jpl-damage/src/c/shared/Elements/Paterson.cpp@ 12878

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

#include <math.h>

#include "../../include/include.h"

IssmDouble Paterson(IssmDouble temperature){
        
        /*output: */
        IssmDouble B;
        IssmDouble T;

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

//      %The routine below is equivalent to:
//      % n=3; T=temperature-273;
//      % %From Paterson,
//      % Temp=[0;-2;-5;-10;-15;-20;-25;-30;-35;-40;-45;-50];
//      % A=[6.8*10^-15;2.4*10^-15;1.6*10^-15;4.9*10^-16;2.9*10^-16;1.7*10^-16;9.4*
//      % 10^-17;5.1*10^-17;2.7*10^-17;1.4*10^-17;7.3*10^-18;3.6*10^-18];;%s-1(kPa-3)
//      % %Convert into B B
//      % B=A.^(-1/n)*10^3; %s^(1/3)Pa
//      % %Now, do a cubic fit between Temp and B: 
//      % fittedmodel=fit(Temp,B,'cubicspline');
//      % B=fittedmodel(temperature);


        if(T<=-45.0){
                B=pow((IssmPDouble)10,(IssmPDouble)8)*(-0.000292866376675*pow(T+50,3)+ 0.011672640664130*pow(T+50,2)  -0.325004442485481*(T+50)+  6.524779401948101);
        }
        else if((T>=-45.0) && (T<=-40.0)){
                B=pow((IssmPDouble)10,(IssmPDouble)8)*(-0.000292866376675*pow(T+45,3)+ 0.007279645014004*pow(T+45,2)  -0.230243014094813*(T+45)+  5.154964909039554);
        }
        else if((T>=-40.0) && (T<=-35.0)){
                B=pow((IssmPDouble)10,(IssmPDouble)8)*(0.000072737147457*pow(T+40,3)+  0.002886649363879*pow(T+40,2)  -0.179411542205399*(T+40)+  4.149132666831214);
        }
        else if((T>=-35.0) && (T<=-30.0)){
                B=pow((IssmPDouble)10,(IssmPDouble)8)*(-0.000086144770023*pow(T+35,3)+ 0.003977706575736*pow(T+35,2)  -0.145089762507325*(T+35)+  3.333333333333331);
        }
        else if((T>=-30.0) && (T<=-25.0)){
                B=pow((IssmPDouble)10,(IssmPDouble)8)*(-0.000043984685769*pow(T+30,3)+ 0.002685535025386*pow(T+30,2)  -0.111773554501713*(T+30)+  2.696559088937191);
        }
        else if((T>=-25.0) && (T<=-20.0)){
                B=pow((IssmPDouble)10,(IssmPDouble)8)*(-0.000029799523463*pow(T+25,3)+ 0.002025764738854*pow(T+25,2)  -0.088217055680511*(T+25)+  2.199331606342181);
        }
        else if((T>=-20.0) && (T<=-15.0)){
                B=pow((IssmPDouble)10,(IssmPDouble)8)*(0.000136920904777*pow(T+20,3)+  0.001578771886910*pow(T+20,2)  -0.070194372551690*(T+20)+  1.805165505978111);
        }
        else if((T>=-15.0) && (T<=-10.0)){
                B=pow((IssmPDouble)10,(IssmPDouble)8)*(-0.000899763781026*pow(T+15,3)+ 0.003632585458564*pow(T+15,2)  -0.044137585824322*(T+15)+  1.510778053489523);
        }
        else if((T>=-10.0) && (T<=-5.0)){
                B=pow((IssmPDouble)10,(IssmPDouble)8)*(0.001676964325070*pow(T+10,3)-  0.009863871256831*pow(T+10,2)  -0.075294014815659*(T+10)+  1.268434288203714);
        }
        else if((T>=-5.0) && (T<=-2.0)){
                B=pow((IssmPDouble)10,(IssmPDouble)8)*(-0.003748937622487*pow(T+5,3)+0.015290593619213*pow(T+5,2)  -0.048160403003748*(T+5)+  0.854987973338348);
        }
        else if(T>=-2.0){
                B=pow((IssmPDouble)10,(IssmPDouble)8)*(-0.003748937622488*pow(T+2,3)-0.018449844983174*pow(T+2,2)  -0.057638157095631*(T+2)+  0.746900791092860);
        }

        /*B cannot be negative!*/
        if(B<0) B=pow((IssmPDouble)10,(IssmPDouble)6);

        return B;
}