#include "physics.hpp"
//------------------------
// Brooks and Corey model
//------------------------

namespace Physics{
  double g=9.81;
  double rho=1000;
  double phi=0.3;
  double k0=3e-5;
  double nivrivsat=0.01;
  double (*s)(double)=&s_BC;
  double (*ds)(double)=&ds_BC;
  void (*s_ds)(double,double&,double&)=&s_ds_BC;
  double (*kr)(double)=&kr_BC;
  double (*dkr)(double)=&dkr_BC;
  void (*kr_dkr)(double,double&,double&)=&kr_dkr_BC;
  double model_data[6]={-2000,0,3,11,0,0};
  
#define Psat model_data[0]
#define sres model_data[1]
#define lambda model_data[2]
#define alpha model_data[3]

  void setModel(Model model){
    switch(model){
    case BrooksCorey:
      s=&s_BC;
      ds=&ds_BC;
      s_ds=&s_ds_BC;
      kr=&kr_BC;
      dkr=&dkr_BC;
      kr_dkr=&kr_dkr_BC;
      Psat=-2000;
      sres=0;
      lambda=3;
      alpha=11;
      break;
    default:
      assert(false);
    };
  }
  
  
  double
  s_BC(double P){
    if(P>=Psat) return 1;
    return sres+(1-sres)*pow(Psat/P,lambda);
  }
  
  double
  ds_BC(double P){
    if(P>=Psat) return 0;
    return ((sres-1)*lambda*pow(Psat/P,lambda))/P;
  }
  
  void
  s_ds_BC(double P,double& v,double& dv){
    if(P>=Psat){
      v=1;
      dv=0;
    }
    else{
      double t=(1-sres)*pow(Psat/P,lambda);
      v=sres+t;
      dv=-(lambda*t)/P;
    }
  }
  
  double
  kr_BC(double P){
    if(P>=Psat) return 1;
    return pow(Psat/P,alpha);
  }

  double
  dkr_BC(double P){
    if(P>=Psat) return 0;
    return -alpha*pow(Psat/P,alpha)/P;
  }

  void
  kr_dkr_BC(double P,double& v,double& dv){
    if(P>=Psat){
      v=1;
      dv=0;
    }
    else{
      double t=pow(Psat/P,alpha);
      v=t;
    dv=-(alpha*t)/P;
    }
  }
}