def stetson_jindex(ftimes, fmags, ferrs, weightbytimediff=False):
ndet = len(fmags)<EOL>if ndet > <NUM_LIT:9>:<EOL><INDENT>medmag = npmedian(fmags)<EOL>delta_prefactor = (ndet/(ndet - <NUM_LIT:1>))<EOL>sigma_i = delta_prefactor*(fmags - medmag)/ferrs<EOL>sigma_j = nproll(sigma_i,<NUM_LIT:1>)<EOL>if weightbytimediff:<EOL><INDENT>difft = npdiff(ftimes)<EOL>deltat = npmedian(difft)<EOL>weights_i = npexp(- difft/deltat )<EOL>products = (weights_i*sigma_i[<NUM_LIT:1>:]*sigma_j[<NUM_LIT:1>:])<EOL><DEDENT>else:<EOL><INDENT>products = (sigma_i*sigma_j)[<NUM_LIT:1>:]<EOL><DEDENT>stetsonj = (<EOL>npsum(npsign(products) * npsqrt(npabs(products)))<EOL>) / ndet<EOL>return stetsonj<EOL><DEDENT>else:<EOL><INDENT>LOGERROR('<STR_LIT>'<EOL>'<STR_LIT>')<EOL>return npnan<EOL><DEDENT>
This calculates the Stetson index for the magseries, based on consecutive
    pairs of observations.

    Based on Nicole Loncke's work for her Planets and Life certificate at
    Princeton in 2014.

    Parameters
    ----------

    ftimes,fmags,ferrs : np.array
        The input mag/flux time-series with all non-finite elements removed.

    weightbytimediff : bool
        If this is True, the Stetson index for any pair of mags will be
        reweighted by the difference in times between them using the scheme in
        Fruth+ 2012 and Zhange+ 2003 (as seen in Sokolovsky+ 2017)::

            w_i = exp(- (t_i+1 - t_i)/ delta_t )

    Returns
    -------

    float
        The calculated Stetson J variability index.