Loop over the first axis applying Fowler processing.

    *fowlerdata* is assumed to be a 3D numpy.ndarray containing the
    result of a nIR observation in Fowler mode (Fowler and Gatley 1991).
    The shape of the array must be of the form 2N_p x M x N, with N_p being
    the number of pairs in Fowler mode.

    The output signal is just the mean value of the differences between the
    last N_p values (S_i) and the first N_p values (R-i).

    .. math::

        S_F = \\frac{1}{N_p}\\sum\\limits_{i=0}^{N_p-1} S_i - R_i


    If the source has a radiance F, then the measured signal is equivalent
    to:

    .. math::

        S_F = F T_I - F T_S (N_p -1) = F T_E

    being T_I the integration time (*ti*), the time since the first
    productive read to the last productive read for a given pixel and T_S the
    time between samples (*ts*). T_E is the time between correlated reads
    :math:`T_E = T_I - T_S (N_p - 1)`.

    The variance of the signnal is the sum of two terms, one for the readout
    noise:

    .. math::

        \\mathrm{var}(S_{F1}) =\\frac{2\sigma_R^2}{N_p}

    and other for the photon noise:

    .. math::

        \\mathrm{var}(S_{F2}) = F T_E - F T_S \\frac{1}{3}(N_p-\\frac{1}{N_p})
        = F T_I - F T_S (\\frac{4}{3} N_p -1 -  \\frac{1}{3N_p})


    :param fowlerdata: Convertible to a 3D numpy.ndarray with first axis even
    :param ti: Integration time.
    :param ts: Time between samples.
    :param gain: Detector gain.
    :param ron: Detector readout noise in counts.
    :param badpixels: An optional MxN mask of dtype 'uint8'.
    :param dtype: The dtype of the float outputs.
    :param saturation: The saturation level of the detector.
    :param blank: Invalid values in output are substituted by *blank*.
    :returns: A tuple of (signal, variance of the signal, numper of pixels used
        and badpixel mask.
    :raises: ValueError

def fowler_array(fowlerdata, ti=0.0, ts=0.0, gain=1.0, ron=1.0,
                 badpixels=None, dtype='float64',
                 saturation=65631, blank=0, normalize=False):
    """Loop over the first axis applying Fowler processing.

    *fowlerdata* is assumed to be a 3D numpy.ndarray containing the
    result of a nIR observation in Fowler mode (Fowler and Gatley 1991).
    The shape of the array must be of the form 2N_p x M x N, with N_p being
    the number of pairs in Fowler mode.

    The output signal is just the mean value of the differences between the
    last N_p values (S_i) and the first N_p values (R-i).

    .. math::

        S_F = \\frac{1}{N_p}\\sum\\limits_{i=0}^{N_p-1} S_i - R_i


    If the source has a radiance F, then the measured signal is equivalent
    to:

    .. math::

        S_F = F T_I - F T_S (N_p -1) = F T_E

    being T_I the integration time (*ti*), the time since the first
    productive read to the last productive read for a given pixel and T_S the
    time between samples (*ts*). T_E is the time between correlated reads
    :math:`T_E = T_I - T_S (N_p - 1)`.

    The variance of the signnal is the sum of two terms, one for the readout
    noise:

    .. math::

        \\mathrm{var}(S_{F1}) =\\frac{2\sigma_R^2}{N_p}

    and other for the photon noise:

    .. math::

        \\mathrm{var}(S_{F2}) = F T_E - F T_S \\frac{1}{3}(N_p-\\frac{1}{N_p})
        = F T_I - F T_S (\\frac{4}{3} N_p -1 -  \\frac{1}{3N_p})


    :param fowlerdata: Convertible to a 3D numpy.ndarray with first axis even
    :param ti: Integration time.
    :param ts: Time between samples.
    :param gain: Detector gain.
    :param ron: Detector readout noise in counts.
    :param badpixels: An optional MxN mask of dtype 'uint8'.
    :param dtype: The dtype of the float outputs.
    :param saturation: The saturation level of the detector.
    :param blank: Invalid values in output are substituted by *blank*.
    :returns: A tuple of (signal, variance of the signal, numper of pixels used
        and badpixel mask.
    :raises: ValueError

    """
    import numina.array._nirproc as _nirproc

    if gain <= 0:
        raise ValueError("invalid parameter, gain <= 0.0")

    if ron <= 0:
        raise ValueError("invalid parameter, ron < 0.0")

    if ti < 0:
        raise ValueError("invalid parameter, ti < 0.0")

    if ts < 0:
        raise ValueError("invalid parameter, ts < 0.0")

    if saturation <= 0:
        raise ValueError("invalid parameter, saturation <= 0")

    fowlerdata = numpy.asarray(fowlerdata)

    if fowlerdata.ndim != 3:
        raise ValueError('fowlerdata must be 3D')

    npairs = fowlerdata.shape[0] // 2
    if 2 * npairs != fowlerdata.shape[0]:
        raise ValueError('axis-0 in fowlerdata must be even')

    # change byteorder
    ndtype = fowlerdata.dtype.newbyteorder('=')
    fowlerdata = numpy.asarray(fowlerdata, dtype=ndtype)
    # type of the output
    fdtype = numpy.result_type(fowlerdata.dtype, dtype)
    # Type of the mask
    mdtype = numpy.dtype('uint8')

    fshape = (fowlerdata.shape[1], fowlerdata.shape[2])

    if badpixels is None:
        badpixels = numpy.zeros(fshape, dtype=mdtype)
    else:
        if badpixels.shape != fshape:
            raise ValueError('shape of badpixels is not '
                             'compatible with shape of fowlerdata')
        if badpixels.dtype != mdtype:
            raise ValueError('dtype of badpixels must be uint8')

    result = numpy.empty(fshape, dtype=fdtype)
    var = numpy.empty_like(result)
    npix = numpy.empty(fshape, dtype=mdtype)
    mask = badpixels.copy()

    _nirproc._process_fowler_intl(
        fowlerdata, ti, ts,  gain, ron,
        badpixels, saturation, blank,
        result, var, npix, mask
    )
    return result, var, npix, mask