// Demonstrate slow TF1 constructor
//
// 20210525  Michael Kelsey

#include "TF1.h"
#include <stdlib.h>
#include <time.h>
#include <iostream>
using namespace std;


double T5diff(double a, double b) {     // a^5-b^5 avoid cancellations
  double b2=b*b, b3=b2*b, b4=b3*b;
  return (a-b) * (a*(a*(a*(a+b) + b2) + b3) + b4);
}

double Tdot(double* x, double* par) { // x = IoverIc
  const double& IoverIc = x[0];

  const double& A = par[0];
  const double& Ksubst = par[1];
  const double& Tsubst = par[2];
  const double& Tw = par[3];
  const double& Tc = par[4];
  const double& Ic = par[5];
  const double& Rn = par[6];
  const double& heatcap = par[7];
  const double& nsc = par[8];
  const double& PP = par[9];
  const int verboseLevel = int(par[10]);

  if (heatcap <= 0.) return 0.;

  double Ti = M_LN2*Tw*A+Tc*pow(std::max(0.,1.-fabs(IoverIc)),nsc);
  double KdeltaT5 = Ksubst * T5diff(Tsubst,Ti);
  double I2R = IoverIc*Ic*IoverIc*Ic*Rn/(1+exp(-2.*A));

  return (KdeltaT5 + I2R + PP)/(M_LN2*Tw*heatcap);
}

int main() {
  cout << "Demonstrating slow TF1 constructor" << endl;

  time_t my_time;

  my_time = time(NULL);
  cout << "Before ctor: " << ctime(&my_time) << endl;

  TF1 f1("Function", Tdot, 0., 6., 12);

  my_time = time(NULL);
  cout << "After  ctor: " << ctime(&my_time) << endl;
}