import numpy as np

def reza_band_pass_filter(
    data,
    fs,
    fc1,
    fc2,
    c=0.9,
    offset=1.0,
    d=None,
    initial_d=10.0,
    d_increment=5.0,
    threshold=1e-4,
    max_iter=200,
    axis=-1,
    return_params=False,
):
    """
    Reza Exponential BAND-PASS (frequency-domain magnitude template).
      hp_gain(f) = 1                                  for f >= fc1
                = exp( -c * ((fc1 - f) + offset)^d )  for f <  fc1
      lp_gain(f) = 1                                  for f <= fc2
                = exp( -c * ((f - fc2) + offset)^d )  for f >  fc2
      gain = hp_gain * lp_gain

    If d is None, d is auto-tuned using two-edge sharpness stabilization.
    """
    x = np.asarray(data, dtype=float)

    # ---- local helpers ----
    def _validate():
        if fs is None or fs <= 0:
            raise ValueError(f"fs must be positive. Got fs={fs}.")
        if not isinstance(axis, int):
            raise ValueError("axis must be an int.")
        if axis < -x.ndim or axis >= x.ndim:
            raise ValueError(f"axis={axis} out of bounds for ndim={x.ndim}.")
        nyq = 0.5 * fs
        if not (0.0 < fc1 < fc2 < nyq):
            raise ValueError(f"Require 0 < fc1 < fc2 < fs/2. Got fc1={fc1}, fc2={fc2}, fs/2={nyq}.")

    def _gain_bp(freqs, fc1_, fc2_, c_, off_, d_):
        hp = np.where(
            freqs >= fc1_,
            1.0,
            np.exp(-c_ * ((fc1_ - freqs) + off_) ** d_),
        )
        lp = np.where(
            freqs <= fc2_,
            1.0,
            np.exp(-c_ * ((freqs - fc2_) + off_) ** d_),
        )
        return hp * lp

    def _auto_tune_d(freqs):
        d_val = float(initial_d)
        last_sharp = 0.0
        sharp_impr = np.inf

        for _ in range(int(max_iter)):
            g = _gain_bp(freqs, fc1, fc2, c, offset, d_val)

            idx1 = int(np.searchsorted(freqs, fc1))
            idx2 = int(np.searchsorted(freqs, fc2))
            idx1 = min(max(idx1, 0), len(g) - 1)
            idx2 = min(max(idx2, 0), len(g) - 1)

            idx1_lo = max(idx1 - 1, 0)
            idx2_hi = min(idx2 + 1, len(g) - 1)

            edge1 = abs(g[idx1] - g[idx1_lo])   # near fc1
            edge2 = abs(g[idx2] - g[idx2_hi])   # near fc2
            sharp = 0.5 * (edge1 + edge2)

            sharp_impr = abs(sharp - last_sharp)
            if sharp_impr > threshold:
                last_sharp = sharp
                d_val += d_increment
            else:
                break

        return d_val

    # ---- main ----
    _validate()

    x_m = np.moveaxis(x, axis, -1)
    n = x_m.shape[-1]
    if n < 2:
        raise ValueError("data length along filtering axis must be >= 2.")

    freqs = np.fft.rfftfreq(n, d=1.0 / fs)
    d_used = float(d) if d is not None else _auto_tune_d(freqs)

    gain = _gain_bp(freqs, fc1, fc2, c, offset, d_used)
    X = np.fft.rfft(x_m, axis=-1)
    Y = X * gain
    y_m = np.fft.irfft(Y, n=n, axis=-1)
    y = np.moveaxis(y_m, -1, axis)

    if return_params:
        return y, {"d": d_used, "fc1": fc1, "fc2": fc2, "c": c, "offset": offset}
    return y
