
//p[0]:  A        A
//p[1]:  tb1      t_{b,1}
//p[2]:  omega1   \omega_1
//p[3]:  alpha1   \alpha_1
//p[4]:  alpha2   \alpha_2
//p[5]:  alpha0   \alpha_0
//p[6]:  alpha3   \alpha_3
//p[7]:  tb2      t_{b,2}
//p[8]:  omega2   \omega_2
//p[9]:  tb0      t_{b,0}
//p[10]: tstar    t^\aster


long double lfunc_f12(double *x, double *p) {
  long double tstar = p[10];
  long double t = x[0] - tstar;
  if (t<0.01) return 0.0L;

  long double A = p[0];
  long double tb1 = p[1];
  long double omega1 = p[2];
  long double alpha1 = p[3];
  long double alpha2 = p[4];

  return A*exp(
     log(pow(t/tb1,-omega1*alpha1) + pow(t/tb1,-omega1*alpha2))*
    (-1.0L/omega1));
}


long double lfunc_f0(double *x, double *p) {
  long double tstar = p[10];
  long double t = x[0] - tstar;
  if (t<0.01) return 0.0L;

  long double alpha0 = p[5];
  long double tb0 = p[9];

  double xp[1];
  xp[0] = tb0 + tstar;
  long double f12 = lfunc_f12(xp,p);

  return f12*pow(t/tb0,alpha0);
}


long double lfunc_f3(double *x, double *p) {
  long double tstar = p[10];
  long double t = x[0] - tstar;
  if (t<0.01) return 0.0L;

  long double alpha3 = p[6];
  long double tb2 = p[7];

  double xp[1];
  xp[0] = tb2 + tstar;
  long double f12 = lfunc_f12(xp,p);

  return f12*exp(log(t/tb2)*alpha3);
}


long double lfunc_f123(double *x, double *p) {
  long double tstar = p[10];
  long double t = x[0] - tstar;
  if (t<0.01) return 0.0L;

  long double omega2 = p[8];

  long double f12 = lfunc_f12(x,p);
  long double f3 = lfunc_f3(x,p);

  if (f12<=f3) return f12*pow(1.0L+pow(f3/f12,-omega2),-1.0L/omega2);
  else return f3*pow(pow(f12/f3,-omega2)+1.0L,-1.0L/omega2);
}


long double lfunc_f012(double *x, double *p) {
  long double tstar = p[10];
  long double t = x[0] - tstar;
  if (t<0.01) return 0.0L;

  long double tb0 = p[9];

  if (t<tb0) return lfunc_f0(x,p);
  else return lfunc_f12(x,p);
}


long double lfunc_f0123(double *x, double *p) {
  long double tstar = p[10];
  long double t = x[0] - tstar;
  if (t<0.01) return 0.0L;

  long double omega2 = p[8];

  long double f012 = lfunc_f012(x,p);
  long double f3 = lfunc_f3(x,p);

  if (f012<=f3) return f012*pow(1.0L+pow(f3/f012,-omega2),-1.0L/omega2);
  else return f3*pow(pow(f012/f3,-omega2)+1.0L,-1.0L/omega2);
}


double func_f0(double *x, double *p) { return lfunc_f0(x,p); }
double func_f3(double *x, double *p) { return lfunc_f3(x,p); }
double func_f12(double *x, double *p) { return lfunc_f12(x,p); }
double func_f012(double *x, double *p) { return lfunc_f012(x,p); }
double func_f123(double *x, double *p) { return lfunc_f123(x,p); }
double func_f0123(double *x, double *p) { return lfunc_f0123(x,p); }


