def psi0_psi1(lhi, rh_rt, zft):
    zfh = rh_rt * zft
    if lhi in psi0.index.values:  # use calculated data.
        # psi0
        psi0i = psi0.loc[lhi] * rh_rt
        psi0i_use = interp_psi(psi0i, zfh, kind='quadratic')
        # psi1
        psi1i = psi1.loc[lhi] * rh_rt ** 3
        psi1i_use = interp_psi(psi1i, zfh, kind='quadratic')
    else:  # fit the trend of parameters.
        #         itp_kind = 'linear' if zfh < 0.05 else 'quadratic'
        # psi0
        int_a0 = interp1d(psi0.columns.values, fit_psi0[:, 0], kind='quadratic',
                          fill_value='extrapolate')
        int_a1 = interp1d(psi0.columns.values, fit_psi0[:, 1], kind='quadratic',
                          fill_value='extrapolate')
        psi0i_use = func_psi0(lhi, int_a0(zfh), int_a1(zfh)) * rh_rt
        # psi1
        int_a0 = interp1d(psi1.columns.values, fit_psi1[:, 0], kind='quadratic',
                          fill_value='extrapolate')
        int_a1 = interp1d(psi1.columns.values, fit_psi1[:, 1], kind='quadratic',
                          fill_value='extrapolate')
        psi1i_use = func_psi1(lhi, int_a0(zfh), int_a1(zfh)) * rh_rt ** 3
    #         if zft == 0:
    #             psi1i_use =  func_psi1(lhi, *fit_psi1[0, :]) * rh_rt ** 3
    #         else:
    #             t1 = 1 / zfh
    #             psi1i_use = 4 * np.pi * (t1**2) / (t1**2 - 1) * lhi * rh_rt ** 3
    return psi0i_use, psi1i_use


def intp_psi_tail(psi, phi, chi, zfi):
    tpsi = psi.loc[phi][zfi]
    tx = tpsi.index
    ty = tpsi.values
    idxi = np.isfinite(tx) & np.isfinite(ty)
    tint = interp1d(tx[idxi], ty[idxi], kind='quadratic', fill_value='extrapolate')
    return tint(chi)


def interp_psi(psi, zf, kind='quadratic'):
    def interp_once(psi, zf, kind='quadratic'):
        tx = psi.index
        ty = psi.values
        idxi = np.isfinite(tx) & np.isfinite(ty)
        int_psi = interp1d(tx[idxi], ty[idxi], kind=kind, fill_value='extrapolate')
        psi_zf = int_psi(zf)
        return psi_zf

    if isinstance(psi, pd.core.series.Series):
        psi_zf = interp_once(psi, zf, kind)
    elif isinstance(psi, pd.core.frame.DataFrame):
        psi_zf = psi.copy()
        for phi, chi in psi.index.values:
            t_psi = psi.loc[phi].loc[chi]
            t_psi_zf = interp_once(t_psi, zf, kind)
            psi_zf.loc[phi].loc[chi][:] = t_psi_zf
    return psi_zf


def fit_psi_tail(psi, phi, chi, zfi, chmin, chmax):
    x = psi.loc[phi][zfi].index
    y = psi.loc[phi][zfi].values
    idx = np.array(x > chmin) & np.array(x < chmax) & np.isfinite(x) & np.isfinite(y)
    popt_line, pcov = curve_fit(spf.func_line, x[idx], y[idx], maxfev=10000)
    return spf.func_line(chi, *popt_line)


def fit_psi20_psi60_old(psi, phi, chi, chmin=3, chmax=np.inf):
    def func_psi0(x, a0, a1):
        # y = a0 * x / np.log(x) + a1
        y = a0 * np.log(np.exp(x) - x) + a1
        return y

    # mathematically func_psi0 and func_psi0_2 are same, but numerically, func_psi0 is better for
    #    fitting process, while func_psi0_2 is better for calculate larger chi
    def func_psi0_2(x, a0, a1):
        y = a0 * x / np.log(x) + a1
        return y

    tpsi = psi.loc[phi][0]
    tx = tpsi.index
    ty = tpsi.values
    idxi = np.isfinite(tx) & np.isfinite(ty) & (tx > chmin) & (tx < chmax)
    fit_psi0_line = spf.fit_line(None, tx[idxi], ty[idxi], 3, np.inf, ifprint=False)
    fit_psi0, pcov = curve_fit(func_psi0, tx[idxi], ty[idxi], p0=fit_psi0_line,
                               bounds=((0, 0), (np.inf, np.inf)), maxfev=10000)
    return func_psi0_2(chi, *fit_psi0)


def fit_psi20_psi60(psi, phi, chi, chmin=3, chmax=np.inf):
    # try to fit x/ln(x) = b0*y + b1, so we have {a0 = 1/b0, a1 = -b1/b0}
    #   where y = a0 * x / np.log(x) + a1.
    tpsi = psi.loc[phi][0]
    tx = tpsi.index
    ty = tpsi.values
    idxi = np.isfinite(tx) & np.isfinite(ty) & (chmin <= tx) & (tx <= chmax)
    fit_x = ty[idxi]
    fit_y = tx[idxi] / np.log(tx[idxi])
    b0, b1 = spf.fit_line(None, fit_x, fit_y, 0, np.inf, ifprint=False)
    a0 = 1 / b0
    a1 = -b1 / b0
    return a0 * chi / np.log(chi) + a1


def psi20_psi30_psi60(phi, chi):
    def func_psi_tail_0(x, a0, a1):
        # fit function of psi2 and psi6 in bulk fluid
        # y = a0 + a1 * x / np.log(x)
        y = a0 + a1 * np.log(np.exp(x) - x)
        return y

    if psi2.index.levels[1].values.min() <= chi <= psi2.index.levels[1].values.max():
        # in the range, intepolate
        psi2i_0 = intp_psi_tail(psi2, phi, chi, 0)
        psi3i_0 = intp_psi_tail(psi3, phi, chi, 0)
        psi6i_0 = -intp_psi_tail(-psi6, phi, chi, 0)
    else:
        # out of the range, fit
        psi2i_0 = fit_psi20_psi60(psi2, phi, chi, 3, np.inf)
        psi3i_0 = fit_psi_tail(psi3, phi, chi, 0, 3, np.inf)
        psi6i_0 = -fit_psi20_psi60(-psi6, phi, chi, 3, np.inf)
    return psi2i_0, psi3i_0, psi6i_0


def psi2_psi3_psi6(phi, chi, zft):
    # intepolate or fit psi2, psi3, psi6
    if psi2.index.levels[1].values.min() <= chi <= psi2.index.levels[1].values.max():
        # in the range, intepolate
        psi2i_use = intp_psi_tail(psi2, phi, chi, zft)
        psi3i_use = intp_psi_tail(psi3, phi, chi, zft)
        psi6i_use = intp_psi_tail(psi6, phi, chi, zft)
    else:
        # out of the range, fit
        psi2i_use = fit_psi_tail(psi2, phi, chi, zft, 4, 8)
        psi3i_use = fit_psi_tail(psi3, phi, chi, zft, 4, 8)
        psi6i_use = fit_psi_tail(psi6, phi, chi, zft, 4, 8)
    #     else:
    #         # out of the range, fit
    #         psi2i_use = fit_psi_tail(psi2, phi, chi, zft, 1, 3)
    #         psi3i_use = fit_psi_tail(psi3, phi, chi, zft, 1, 3)
    #         psi6i_use = fit_psi_tail(psi6, phi, chi, zft, 1, 3)
    return psi2i_use, psi3i_use, psi6i_use

def psi0_psi1(lhi, rh_rt, zft):
    zfh = rh_rt * zft
    if lhi in psi0.index.values:  # use calculated data.
        # psi0
        psi0i = psi0.loc[lhi] * rh_rt
        psi0i_use = interp_psi(psi0i, zfh, kind='quadratic')
        # psi1
        psi1i = psi1.loc[lhi] * rh_rt ** 3
        psi1i_use = interp_psi(psi1i, zfh, kind='quadratic')
    else:  # fit the trend of parameters.
        #         itp_kind = 'linear' if zfh < 0.05 else 'quadratic'
        # psi0
        int_a0 = interp1d(psi0.columns.values, fit_psi0[:, 0], kind='quadratic',
                          fill_value='extrapolate')
        int_a1 = interp1d(psi0.columns.values, fit_psi0[:, 1], kind='quadratic',
                          fill_value='extrapolate')
        psi0i_use = func_psi0(lhi, int_a0(zfh), int_a1(zfh)) * rh_rt
        # psi1
        int_a0 = interp1d(psi1.columns.values, fit_psi1[:, 0], kind='quadratic',
                          fill_value='extrapolate')
        int_a1 = interp1d(psi1.columns.values, fit_psi1[:, 1], kind='quadratic',
                          fill_value='extrapolate')
        psi1i_use = func_psi1(lhi, int_a0(zfh), int_a1(zfh)) * rh_rt ** 3
    #         if zft == 0:
    #             psi1i_use =  func_psi1(lhi, *fit_psi1[0, :]) * rh_rt ** 3
    #         else:
    #             t1 = 1 / zfh
    #             psi1i_use = 4 * np.pi * (t1**2) / (t1**2 - 1) * lhi * rh_rt ** 3
    return psi0i_use, psi1i_use


def intp_psi_tail(psi, phi, chi, zfi):
    tpsi = psi.loc[phi][zfi]
    tx = tpsi.index
    ty = tpsi.values
    idxi = np.isfinite(tx) & np.isfinite(ty)
    tint = interp1d(tx[idxi], ty[idxi], kind='quadratic', fill_value='extrapolate')
    return tint(chi)


def interp_psi(psi, zf, kind='quadratic'):
    def interp_once(psi, zf, kind='quadratic'):
        tx = psi.index
        ty = psi.values
        idxi = np.isfinite(tx) & np.isfinite(ty)
        int_psi = interp1d(tx[idxi], ty[idxi], kind=kind, fill_value='extrapolate')
        psi_zf = int_psi(zf)
        return psi_zf

    if isinstance(psi, pd.core.series.Series):
        psi_zf = interp_once(psi, zf, kind)
    elif isinstance(psi, pd.core.frame.DataFrame):
        psi_zf = psi.copy()
        for phi, chi in psi.index.values:
            t_psi = psi.loc[phi].loc[chi]
            t_psi_zf = interp_once(t_psi, zf, kind)
            psi_zf.loc[phi].loc[chi][:] = t_psi_zf
    return psi_zf


def fit_psi_tail(psi, phi, chi, zfi, chmin, chmax):
    x = psi.loc[phi][zfi].index
    y = psi.loc[phi][zfi].values
    idx = np.array(x > chmin) & np.array(x < chmax) & np.isfinite(x) & np.isfinite(y)
    popt_line, pcov = curve_fit(spf.func_line, x[idx], y[idx], maxfev=10000)
    return spf.func_line(chi, *popt_line)


def fit_psi20_psi60_old(psi, phi, chi, chmin=3, chmax=np.inf):
    def func_psi0(x, a0, a1):
        # y = a0 * x / np.log(x) + a1
        y = a0 * np.log(np.exp(x) - x) + a1
        return y

    # mathematically func_psi0 and func_psi0_2 are same, but numerically, func_psi0 is better for
    #    fitting process, while func_psi0_2 is better for calculate larger chi
    def func_psi0_2(x, a0, a1):
        y = a0 * x / np.log(x) + a1
        return y

    tpsi = psi.loc[phi][0]
    tx = tpsi.index
    ty = tpsi.values
    idxi = np.isfinite(tx) & np.isfinite(ty) & (tx > chmin) & (tx < chmax)
    fit_psi0_line = spf.fit_line(None, tx[idxi], ty[idxi], 3, np.inf, ifprint=False)
    fit_psi0, pcov = curve_fit(func_psi0, tx[idxi], ty[idxi], p0=fit_psi0_line,
                               bounds=((0, 0), (np.inf, np.inf)), maxfev=10000)
    return func_psi0_2(chi, *fit_psi0)


def fit_psi20_psi60(psi, phi, chi, chmin=3, chmax=np.inf):
    # try to fit x/ln(x) = b0*y + b1, so we have {a0 = 1/b0, a1 = -b1/b0}
    #   where y = a0 * x / np.log(x) + a1.
    tpsi = psi.loc[phi][0]
    tx = tpsi.index
    ty = tpsi.values
    idxi = np.isfinite(tx) & np.isfinite(ty) & (chmin <= tx) & (tx <= chmax)
    fit_x = ty[idxi]
    fit_y = tx[idxi] / np.log(tx[idxi])
    b0, b1 = spf.fit_line(None, fit_x, fit_y, 0, np.inf, ifprint=False)
    a0 = 1 / b0
    a1 = -b1 / b0
    return a0 * chi / np.log(chi) + a1


def psi20_psi30_psi60(phi, chi):
    def func_psi_tail_0(x, a0, a1):
        # fit function of psi2 and psi6 in bulk fluid
        # y = a0 + a1 * x / np.log(x)
        y = a0 + a1 * np.log(np.exp(x) - x)
        return y

    if psi2.index.levels[1].values.min() <= chi <= psi2.index.levels[1].values.max():
        # in the range, intepolate
        psi2i_0 = intp_psi_tail(psi2, phi, chi, 0)
        psi3i_0 = intp_psi_tail(psi3, phi, chi, 0)
        psi6i_0 = -intp_psi_tail(-psi6, phi, chi, 0)
    else:
        # out of the range, fit
        psi2i_0 = fit_psi20_psi60(psi2, phi, chi, 3, np.inf)
        psi3i_0 = fit_psi_tail(psi3, phi, chi, 0, 3, np.inf)
        psi6i_0 = -fit_psi20_psi60(-psi6, phi, chi, 3, np.inf)
    return psi2i_0, psi3i_0, psi6i_0


def psi2_psi3_psi6(phi, chi, zft):
    # intepolate or fit psi2, psi3, psi6
    if psi2.index.levels[1].values.min() <= chi <= psi2.index.levels[1].values.max():
        # in the range, intepolate
        psi2i_use = intp_psi_tail(psi2, phi, chi, zft)
        psi3i_use = intp_psi_tail(psi3, phi, chi, zft)
        psi6i_use = intp_psi_tail(psi6, phi, chi, zft)
    else:
        # out of the range, fit
        psi2i_use = fit_psi_tail(psi2, phi, chi, zft, 4, 8)
        psi3i_use = fit_psi_tail(psi3, phi, chi, zft, 4, 8)
        psi6i_use = fit_psi_tail(psi6, phi, chi, zft, 4, 8)
    #     else:
    #         # out of the range, fit
    #         psi2i_use = fit_psi_tail(psi2, phi, chi, zft, 1, 3)
    #         psi3i_use = fit_psi_tail(psi3, phi, chi, zft, 1, 3)
    #         psi6i_use = fit_psi_tail(psi6, phi, chi, zft, 1, 3)
    return psi2i_use, psi3i_use, psi6i_use

# use calculated data. # psi0 # psi1 # fit the trend of parameters. #         itp_kind = 'linear' if zfh < 0.05 else 'quadratic' # psi0 # psi1 #         if zft == 0: #             psi1i_use =  func_psi1(lhi, *fit_psi1[0, :]) * rh_rt ** 3 #         else: #             t1 = 1 / zfh #             psi1i_use = 4 * np.pi * (t1**2) / (t1**2 - 1) * lhi * rh_rt ** 3 # y = a0 * x / np.log(x) + a1 # mathematically func_psi0 and func_psi0_2 are same, but numerically, func_psi0 is better for #    fitting process, while func_psi0_2 is better for calculate larger chi # try to fit x/ln(x) = b0*y + b1, so we have {a0 = 1/b0, a1 = -b1/b0} #   where y = a0 * x / np.log(x) + a1. # fit function of psi2 and psi6 in bulk fluid # y = a0 + a1 * x / np.log(x) # in the range, intepolate # out of the range, fit # intepolate or fit psi2, psi3, psi6 # in the range, intepolate # out of the range, fit #     else: #         # out of the range, fit #         psi2i_use = fit_psi_tail(psi2, phi, chi, zft, 1, 3) #         psi3i_use = fit_psi_tail(psi3, phi, chi, zft, 1, 3) #         psi6i_use = fit_psi_tail(psi6, phi, chi, zft, 1, 3)