def multiInterp(x, xp, fp, extrap='bounded'):
    j = (np.searchsorted(xp, x) - 1)
    dd = np.zeros(len(x))
    bOK = np.logical_and((j >= 0), (j < (len(xp) - 1)))
    bLower = (j < 0)
    bUpper = (j >= (len(xp) - 1))
    jOK = j[bOK]
    dd[bOK] = ((x[bOK] - xp[jOK]) / (xp[(jOK + 1)] - xp[jOK]))
    jBef = j
    jAft = (j + 1)
    jAft[bUpper] = (len(xp) - 1)
    jBef[bUpper] = (len(xp) - 1)
    jAft[bLower] = 0
    jBef[bLower] = 0
    if (extrap == 'bounded'):
        pass
    elif (extrap == 'nan'):
        dd[(~ bOK)] = np.nan
    else:
        raise NotImplementedError()
    return (((1 - dd) * fp[:, jBef]) + (fp[:, jAft] * dd))
def multiInterp(x, xp, fp, extrap='bounded'):
    j = (np.searchsorted(xp, x) - 1)
    dd = np.zeros(len(x))
    bOK = np.logical_and((j >= 0), (j < (len(xp) - 1)))
    bLower = (j < 0)
    bUpper = (j >= (len(xp) - 1))
    jOK = j[bOK]
    dd[bOK] = ((x[bOK] - xp[jOK]) / (xp[(jOK + 1)] - xp[jOK]))
    jBef = j
    jAft = (j + 1)
    jAft[bUpper] = (len(xp) - 1)
    jBef[bUpper] = (len(xp) - 1)
    jAft[bLower] = 0
    jBef[bLower] = 0
    if (extrap == 'bounded'):
        pass
    elif (extrap == 'nan'):
        dd[(~ bOK)] = np.nan
    else:
        raise NotImplementedError()
    return (((1 - dd) * fp[:, jBef]) + (fp[:, jAft] * dd))