def initialize_bps(
    V, pixelsperkm, timestep, p_par=None, p_perp=None, randstate=None, seed=None
):
    """
    Initialize the motion field perturbator described in :cite:`BPS2006`.
    For simplicity, the bias adjustment procedure described there has not been
    implemented. The perturbator generates a field whose magnitude increases
    with respect to lead time.

    Parameters
    ----------
    V: array_like
      Array of shape (2,m,n) containing the x- and y-components of the m*n
      motion field to perturb.
    p_par: tuple
      Tuple containing the parameters a,b and c for the standard deviation of
      the perturbations in the direction parallel to the motion vectors. The
      standard deviations are modeled by the function f_par(t) = a*t**b+c,
      where t is lead time. The default values are taken from :cite:`BPS2006`.
    p_perp: tuple
      Tuple containing the parameters a,b and c for the standard deviation of
      the perturbations in the direction perpendicular to the motion vectors.
      The standard deviations are modeled by the function f_par(t) = a*t**b+c,
      where t is lead time. The default values are taken from :cite:`BPS2006`.
    pixelsperkm: float
      Spatial resolution of the motion field (pixels/kilometer).
    timestep: float
      Time step for the motion vectors (minutes).
    randstate: mtrand.RandomState
      Optional random generator to use. If set to None, use numpy.random.
    seed: int
      Optional seed number for the random generator.

    Returns
    -------
    out: dict
      A dictionary containing the perturbator that can be supplied to
      generate_motion_perturbations_bps.

    See also
    --------
    pysteps.noise.motion.generate_bps

    """
    if len(V.shape) != 3:
        raise ValueError("V is not a three-dimensional array")
    if V.shape[0] != 2:
        raise ValueError("the first dimension of V is not 2")

    if p_par is None:
        p_par = get_default_params_bps_par()
    if p_perp is None:
        p_perp = get_default_params_bps_perp()

    if len(p_par) != 3:
        raise ValueError("the length of p_par is not 3")
    if len(p_perp) != 3:
        raise ValueError("the length of p_perp is not 3")

    perturbator = {}
    if randstate is None:
        randstate = np.random

    if seed is not None:
        randstate.seed(seed)

    eps_par = randstate.laplace(scale=1.0 / np.sqrt(2))
    eps_perp = randstate.laplace(scale=1.0 / np.sqrt(2))

    # scale factor for converting the unit of the advection velocities
    # into km/h
    vsf = 60.0 / (timestep * pixelsperkm)

    N = linalg.norm(V, axis=0)
    mask = N > 1e-12
    V_n = np.empty(V.shape)
    V_n[:, mask] = V[:, mask] / np.stack([N[mask], N[mask]])
    V_n[:, ~mask] = 0.0

    perturbator["randstate"] = randstate
    perturbator["vsf"] = vsf
    perturbator["p_par"] = p_par
    perturbator["p_perp"] = p_perp
    perturbator["eps_par"] = eps_par
    perturbator["eps_perp"] = eps_perp
    perturbator["V_par"] = V_n
    perturbator["V_perp"] = np.stack([-V_n[1, :, :], V_n[0, :, :]])

    return perturbator
Initialize the motion field perturbator described in :cite:`BPS2006`.
For simplicity, the bias adjustment procedure described there has not been
implemented. The perturbator generates a field whose magnitude increases
with respect to lead time.

Parameters
----------
V: array_like
Array of shape (2,m,n) containing the x- and y-components of the m*n
motion field to perturb.
p_par: tuple
Tuple containing the parameters a,b and c for the standard deviation of
the perturbations in the direction parallel to the motion vectors. The
standard deviations are modeled by the function f_par(t) = a*t**b+c,
where t is lead time. The default values are taken from :cite:`BPS2006`.
p_perp: tuple
Tuple containing the parameters a,b and c for the standard deviation of
the perturbations in the direction perpendicular to the motion vectors.
The standard deviations are modeled by the function f_par(t) = a*t**b+c,
where t is lead time. The default values are taken from :cite:`BPS2006`.
pixelsperkm: float
Spatial resolution of the motion field (pixels/kilometer).
timestep: float
Time step for the motion vectors (minutes).
randstate: mtrand.RandomState
Optional random generator to use. If set to None, use numpy.random.
seed: int
Optional seed number for the random generator.

Returns
-------
out: dict
A dictionary containing the perturbator that can be supplied to
generate_motion_perturbations_bps.

See also
--------
pysteps.noise.motion.generate_bps
pysteps/noise/motion.py