import mmf.math

def f(t, k, eta, theta, b, c):
   beta_m = 1.0/theta
   beta_mu = eta + beta_m
   if 1 == b:
      if 0 == c:
         num = np.sinh(beta_mu)
      else:
         num = np.cosh(beta_mu) + np.exp(-t-beta_m)
   elif 2 == b:
      if 0 == c:
         num = np.sinh(t + beta_m)*np.sinh(beta_mu)
      else:
         num = 1 + np.cosh(beta_mu)*np.cosh(t+beta_m)
   res = (t**k*np.sqrt(1 + t/beta_m/2)*
          num/(np.cosh(t+beta_m) + np.cosh(beta_mu))**b)
   return res

t_0 = np.linspace(0,5,100)
for eta in [0.01,0.5,1.0,2.0,5,100]:
  for theta in [0.01, 1.0, 100.0]:
    for b in [1,2]:
      for c in [0,1]:
        for k in [0.5,1.5,2.5]:
          a = k + 0.5
          t0 = a + mmf.math.LambertW(a*np.exp(eta - a))
          t1 = a + eta
          f0 = a*t0**(k-1)*np.sqrt(1 + t0*theta/2)*np.exp(eta - t0)
          f0 = f(t0, k, eta, theta, b, c)
          plt.plot(t_0, f(t_0*t0, k, eta, theta, b, c)/f0, '-b')
          plt.plot(t_0, f(t_0*t1, k, eta, theta, b, c)/f0, ':y')
plt.xlabel("t/t_0")
plt.ylabel("integrand/integrand(t_0)")