#!/usr/bin/env python

"""
@package ion_functions.qc_functions
@file ion_functions/qc_functions.py
@author Christopher Mueller
@brief Module containing QC functions ported from matlab samples in DPS documents
"""
from ion_functions.qc.qc_extensions import stuckvalues, spikevalues, gradientvalues, ntp_to_month 

import time
import numpy as np
import numexpr as ne
from scipy.interpolate import LinearNDInterpolator
from ion_functions import utils
from ion_functions.utils import fill_value

# try to load the OOI logging module, using default Python logging module if
# unavailable
try:
    from ooi.logging import log
except ImportError:
    import logging
    log = logging.getLogger('ion-functions')

def is_fill(arr):
    return np.atleast_1d(arr)[-1] == -9999. # Not the normal fill value, it's hardcoded to the QC params
def is_none(arr):
    return arr is None or (np.atleast_1d(arr)[-1] == None)

def dataqc_globalrangetest_minmax(dat, dat_min, dat_max, strict_validation=False):
    '''
    Python wrapper for dataqc_globalrangetest
    Combines the min/max arguments into list for dataqc_globalrangetest
    '''
    if is_none(dat_min) or is_none(dat_max) or is_fill(dat_min) or is_fill(dat_max):
        out = np.empty(dat.shape, dtype=np.int8)
        out.fill(-99)
        return out
    return dataqc_globalrangetest(dat, [np.atleast_1d(dat_min)[-1], np.atleast_1d(dat_max)[-1]], strict_validation=strict_validation)


def dataqc_globalrangetest(dat, datlim, strict_validation=False):
    """
    Description:

        Data quality control algorithm testing if measurements fall into a
        user-defined valid range. Returns 1 for presumably good data and 0 for
        data presumed bad.

    Implemented by:

        2010-07-28: DPS authored by Mathias Lankhorst. Example code provided
        for Matlab.
        2013-04-06: Christopher Wingard. Initial python implementation.
        2013-05-30: Christopher Mueller. Performance improvements by adding
            strict_validation flag.

    Usage:

        qcflag = dataqc_globalrangetest(dat, datlim, strict_validation)

            where

        qcflag = Boolean, 0 if value is outside range, else = 1.

        dat = Input dataset, any scalar or vector. Must be numeric and real.
        datlim = Two-element vector with the minimum and maximum values
            considered to be valid.
        strict_validation = Flag (default is False) to assert testing of input
            types (e.g. isreal, isnumeric)

    References:

        OOI (2012). Data Product Specification for Global Range Test. Document
            Control Number 1341-10004. https://alfresco.oceanobservatories.org
            (See: Company Home >> OOI >> Controlled >> 1000 System Level >>
            1341-10004_Data_Product_SPEC_GLBLRNG_OOI.pdf)
    """
    dat = np.atleast_1d(dat)
    datlim = np.atleast_1d(datlim)

    if strict_validation:
        if not utils.isnumeric(dat).all():
            raise ValueError('\'dat\' must be numeric')

        if not utils.isreal(dat).all():
            raise ValueError('\'dat\' must be real')

        if not utils.isnumeric(datlim).all():
            raise ValueError('\'datlim\' must be numeric')

        if not utils.isreal(datlim).all():
            raise ValueError('\'datlim\' must be real')

        if len(datlim) < 2:  # Must have at least 2 elements
            raise ValueError('\'datlim\' must have at least 2 elements')

    return (datlim.min() <= dat) & (dat <= datlim.max()).astype('int8')

def dataqc_localrangetest_wrapper(dat, datlim, datlimz, dims, pval_callback):
    if is_none(datlim) or np.all(np.atleast_1d(datlim).flatten() == -9999):
        out = np.empty(dat.shape, dtype=np.int8)
        out.fill(-99)
        return out
    if is_none(datlimz) or np.all(np.atleast_1d(datlim).flatten() == -9999):
        out = np.empty(dat.shape, dtype=np.int8)
        out.fill(-99)
        return out
    if is_none(dims):
        out = np.empty(dat.shape, dtype=np.int8)
        out.fill(-99)
        return out
    if is_none(pval_callback):
        out = np.empty(dat.shape, dtype=np.int8)
        out.fill(-99)
        return out

    z = []
    for dim in dims:
        if dim == 'month':
            # Convert time vector to vector of months
            v = pval_callback('time')
            v = np.asanyarray(v, dtype=np.float)
            v = ntp_to_month(v)
            z.append(v)
        else:
            # Fetch the dimension from the callback method
            v = pval_callback(dim)
            z.append(v)
    if len(dims)>1:
        z = np.column_stack(z)
    else:
        z = z[0]
        datlimz = datlimz[:,0]
    return dataqc_localrangetest(dat, z, datlim, datlimz)


def dataqc_localrangetest(dat, z, datlim, datlimz, strict_validation=False):
    """
    Description:

        Data quality control algorithm testing if measurements fall into a
        user-defined valid range. This range is not constant but varies with
        measurement location. Returns 1 for presumably good data and 0 for data
        presumed bad.

    Implemented by:

        2012-07-17: DPS authored by Mathias Lankhorst. Example code provided
        for Matlab.
        2013-04-06: Christopher Wingard. Initial python implementation.

    Usage:

        qcflag = dataqc_localrangetest(dat, z, datlim, datlimz)

            where

        qcflag = Boolean, 0 if value is outside range, else = 1.

        dat = input data set, a numeric real scalar or column vector.
        z = location of measurement dat. must have same # of rows as dat and
            same # of columns as datlimz
        datlim = two column array with the minimum (column 1) and maximum
            (column 2) values considered valid.
        datlimz = array with the locations where datlim is given. must have
            same # of rows as datlim and same # of columns as z.

    References:

        OOI (2012). Data Product Specification for Local Range Test. Document
            Control Number 1341-10005. https://alfresco.oceanobservatories.org/
            (See: Company Home >> OOI >> Controlled >> 1000 System Level >>
            1341-10005_Data_Product_SPEC_LOCLRNG_OOI.pdf)
    """
    if strict_validation:
        # check if dat and datlim are matrices
        if not utils.isvector(dat):
            raise ValueError('\'dat\' must be a matrix')

        if not utils.ismatrix(datlim):
            raise ValueError('\'datlim\' must be a matrix')

        # check if all inputs are numeric and real
        for k, arg in {'dat': dat, 'z': z, 'datlim': datlim,
                       'datlimz': datlimz}.iteritems():
            if not utils.isnumeric(arg).all():
                raise ValueError('\'{0}\' must be numeric'.format(k))

            if not utils.isreal(arg).all():
                raise ValueError('\'{0}\' must be real'.format(k))

    if len(datlim.shape) == 3 and datlim.shape[0] == 1:
        datlim = datlim.reshape(datlim.shape[1:])

    if len(datlimz.shape) == 3 and datlimz.shape[0] == 1:
        datlimz = datlimz.reshape(datlimz.shape[1:])

    # test size and shape of the input arrays datlimz and datlim, setting test
    # variables.
    array_size = datlimz.shape
    if len(array_size) == 1:
        numlim = array_size[0]
        ndim = 1
    else:
        numlim = array_size[0]
        ndim = array_size[1]

    array_size = datlim.shape
    tmp1 = array_size[0]
    tmp2 = array_size[1]
    if tmp1 != numlim:
        raise ValueError('\'datlim\' and \'datlimz\' must '
                         'have the same number of rows.')

    if tmp2 != 2:
        raise ValueError('\'datlim\' must be structured as 2-D array '
                         'with exactly 2 columns and 1 through N rows.')

    # test the size and shape of the z input array
    array_size = z.shape
    if len(array_size) == 1:
        num = array_size[0]
        tmp2 = 1
    else:
        num = array_size[0]
        tmp2 = array_size[1]

    if tmp2 != ndim:
        raise ValueError('\'z\' must have the same number of columns '
                         'as \'datlimz\'.')

    if num != dat.size:
        raise ValueError('Len of \'dat\' must match number of '
                         'rows in \'z\'')

    # test datlim, values in column 2 must be greater than those in column 1
    if not all(datlim[:, 1] > datlim[:, 0]):
        raise ValueError('Second column values of \'datlim\' should be '
                         'greater than first column values.')

    # calculate the upper and lower limits for the data set
    if ndim == 1:
        # determine the lower limits using linear interpolation
        lim1 = np.interp(z, datlimz, datlim[:, 0], left=np.nan, right=np.nan)
        # determine the upper limits using linear interpolation
        lim2 = np.interp(z, datlimz, datlim[:, 1], left=np.nan, right=np.nan)
    else:
        # Compute Delaunay Triangulation and use linear interpolation to
        # determine the N-dimensional lower limits
        F = LinearNDInterpolator(datlimz, datlim[:, 0].reshape(numlim, 1))
        lim1 = F(z).reshape(dat.size)

        # Compute Delaunay Triangulation and use linear interpolation to
        # determine the N-dimensional upper limits
        F = LinearNDInterpolator(datlimz, datlim[:, 1].reshape(numlim, 1))
        lim2 = F(z).reshape(dat.size)

    # replace NaNs from above interpolations
    ff = (np.isnan(lim1)) | (np.isnan(lim2))
    lim1[ff] = np.max(datlim[:, 1])
    lim2[ff] = np.min(datlim[:, 0])

    # compute the qcflags
    qcflag = (dat >= lim1) & (dat <= lim2)
    return qcflag.astype('int8')


def dataqc_spiketest_wrapper(dat, acc, N, L, strict_validation=False):
    if is_none(acc) or is_fill(acc) or is_none(N) or is_fill(N) or is_none(L) or is_fill(L):
        out = np.empty(dat.shape, dtype=np.int8)
        out.fill(-99)
        return out
    return dataqc_spiketest(dat, np.atleast_1d(acc)[-1], np.atleast_1d(N)[-1], np.atleast_1d(L)[-1], strict_validation=strict_validation)

def dataqc_spiketest(dat, acc, N=5, L=5, strict_validation=False):
    """
    Description:

        Data quality control algorithm testing a time series for spikes.
        Returns 1 for presumably good data and 0 for data presumed bad.

        The time series is divided into windows of len L (an odd integer
        number). Then, window by window, each value is compared to its (L-1)
        neighboring values: a range R of these (L-1) values is computed (max.
        minus min.), and replaced with the measurement accuracy ACC if ACC>R. A
        value is presumed to be good, i.e. no spike, if it deviates from the
        mean of the (L-1) peers by less than a multiple of the range,
        N*max(R,ACC).

        Further than (L-1)/2 values from the start or end points, the peer
        values are symmetrically before and after the test value. Within that
        range of the start and end, the peers are the first/last L values
        (without the test value itself).

        The purpose of ACC is to restrict spike detection to deviations
        exceeding a minimum threshold value (N*ACC) even if the data have
        little variability. Use ACC=0 to disable this behavior.

    Implemented by:

        2012-07-28: DPS authored by Mathias Lankhorst. Example code provided
        for Matlab.
        2013-04-06: Christopher Wingard. Initial python implementation.
        2013-05-30: Christopher Mueller. Performance optimizations.

    Usage:

        qcflag = dataqc_spiketest(dat, acc, N, L)

            where

        qcflag = Boolean, 0 if value is outside range, else = 1.

        dat = input data set, a numeric, real vector.
        acc = Accuracy of any input measurement.
        N = (optional, defaults to 5) Range multiplier, cf. above
        L = (optional, defaults to 5) Window len, cf. above

    References:

        OOI (2012). Data Product Specification for Spike Test. Document
            Control Number 1341-10006. https://alfresco.oceanobservatories.org/
            (See: Company Home >> OOI >> Controlled >> 1000 System Level >>
            1341-10006_Data_Product_SPEC_SPKETST_OOI.pdf)
    """
    dat = np.atleast_1d(dat)

    if strict_validation:
        if not utils.isnumeric(dat).all():
            raise ValueError('\'dat\' must be numeric')

        if not utils.isreal(dat).all():
            raise ValueError('\'dat\' must be real')

        if not utils.isvector(dat):
            raise ValueError('\'dat\' must be a vector')

        for k, arg in {'acc': acc, 'N': N, 'L': L}.iteritems():
            if not utils.isnumeric(arg).all():
                raise ValueError('\'{0}\' must be numeric'.format(k))

            if not utils.isreal(arg).all():
                raise ValueError('\'{0}\' must be real'.format(k))
    dat = np.asanyarray(dat, dtype=np.float)
    
    out = spikevalues(dat, L, N, acc)
    return out


def dataqc_polytrendtest_wrapper(dat, t, ord_n, nstd, strict_validation=False):
    if is_none(ord_n) or is_fill(ord_n) or is_none(nstd) or is_fill(ord_n):
        out = np.empty(dat.shape, dtype=np.int8)
        out.fill(-99)
        return out
    return dataqc_polytrendtest(dat, t, np.atleast_1d(ord_n)[-1], np.atleast_1d(nstd)[-1], strict_validation=strict_validation)

def dataqc_polytrendtest(dat, t, ord_n=1, nstd=3, strict_validation=False):
    """
    Description:

        Data quality control algorithm testing if measurements contain a
        significant portion of a polynomial. Returns 1 if this is not the case,
        else 0.

        The purpose of this test is to check if a significant fraction of the
        variability in a time series can be explained by a drift, possibly
        interpreted as a sensor drift. This drift is assumed to be a polynomial
        of order ORD. Use ORD=1 to consider a linear drift

        The time series dat is passed to MatLab's POLYFIT routine to obtain a
        polynomial fit PP to dat, and the difference dat-PP is compared to the
        original dat. If the standard deviation of (dat-PP) is less than that
        of dat by a factor of NSTD, the time series is assumed to contain a
        significant trend (output will be 0), else not (output will be 1).

    Implemented by:

        2012-10-29: DPS authored by Mathias Lankhorst. Example code provided
        for Matlab.
        2013-04-06: Christopher Wingard. Initial python implementation.
        2013-05-30: Christopher Mueller. Performance optimizations.

    Usage:

        qcflag = dataqc_polytrendtest(dat, t, ord_n, nstd, strict_validation)

            where

        qcflag = Boolean, 0 a trend is detected, 1 elsewhere.
        dat = Input dataset, a numeric real vector.
        t = time record associated with dat
        ord_n (optional, defaults to 1) = Polynomial order.
        nstd (optional, defaults to 3) = Factor by how much the standard
            deviation must be reduced before qcflag switches from 1 to 0
        strict_validation (optional, defaults to False) = Flag asserting
            testing of inputs.

    References:

        OOI (2012). Data Product Specification for Trend Test. Document
            Control Number 1341-10007. https://alfresco.oceanobservatories.org/
            (See: Company Home >> OOI >> Controlled >> 1000 System Level >>
            1341-10007_Data_Product_SPEC_TRNDTST_OOI.pdf)
    """
    dat = np.atleast_1d(dat)
    t = np.atleast_1d(t)

    if strict_validation:
        for k, arg in {'dat': dat, 't': t, 'ord_n': ord_n, 'nstd': nstd}.iteritems():
            if not utils.isnumeric(arg).all():
                raise ValueError('\'{0}\' must be numeric'.format(k))

            if not utils.isreal(arg).all():
                raise ValueError('\'{0}\' must be real'.format(k))

        for k, arg in {'dat': dat, 't': t}.iteritems():
            if not utils.isvector(arg):
                raise ValueError('\'{0}\' must be a vector'.format(k))

        for k, arg in {'ord_n': ord_n, 'nstd': nstd}.iteritems():
            if not utils.isscalar(arg):
                raise ValueError('\'{0}\' must be a scalar'.format(k))

    ord_n = int(round(abs(ord_n)))
    nstd = int(abs(nstd))

    ll = len(dat)
    # Not needed because time is incorporated as 't'
    # t = range(ll)

    pp = np.polyfit(t, dat, ord_n)
    datpp = np.polyval(pp, t)

    # test for a trend
    if np.atleast_1d((np.std(dat - datpp) * nstd) < np.std(dat)).all():
        trndtst = 0
    else:
        trndtst = 1

    # insure output size equals input, even though test yields a single value.
    qcflag = np.ones(dat.shape).astype('int8') * trndtst
    return qcflag


def dataqc_stuckvaluetest_wrapper(x, reso, num, strict_validation=False):
    if is_none(reso) or is_fill(reso) or is_none(num) or is_fill(num):
        out = np.empty(x.shape, np.int8)
        out.fill(-99)
        return out
    return  dataqc_stuckvaluetest(x, np.atleast_1d(reso)[-1], np.atleast_1d(num)[-1], strict_validation=strict_validation)

def dataqc_stuckvaluetest(x, reso, num=10, strict_validation=False):
    """
    Description:

        Data quality control algorithm testing a time series for "stuck
        values", i.e. repeated occurences of one value. Returns 1 for
        presumably good data and 0 for data presumed bad.

    Implemented by:

        2012-10-29: DPS authored by Mathias Lankhorst. Example code provided
        for Matlab.
        2013-04-06: Christopher Wingard. Initial python implementation.

    Usage:

        qcflag = =dataqc_stuckvaluetest(x, RESO, NUM);

            where

        qcflag = Boolean output: 0 where stuck values are found, 1 elsewhere.
        x = Input time series (vector, numeric).
        reso = Resolution; repeat values less than reso apart will be
            considered "stuck values".
        num = Minimum number of successive values within reso of each other
            that will trigger the "stuck value". num is optional and defaults
            to 10 if omitted or empty.

    References:

        OOI (2012). Data Product Specification for Stuck Value Test. Document
            Control Number 1341-10008. https://alfresco.oceanobservatories.org/
            (See: Company Home >> OOI >> Controlled >> 1000 System Level >>
            1341-10008_Data_Product_SPEC_STUCKVL_OOI.pdf)
    """
    dat = np.atleast_1d(x)

    if strict_validation:
        if not utils.isnumeric(dat).all():
            raise ValueError('\'x\' must be numeric')

        if not utils.isvector(dat):
            raise ValueError('\'x\' must be a vector')

        if not utils.isreal(dat).all():
            raise ValueError('\'x\' must be real')

        for k, arg in {'reso': reso, 'num': num}.iteritems():
            if not utils.isnumeric(arg).all():
                raise ValueError('\'{0}\' must be numeric'.format(k))

            if not utils.isscalar(arg):
                raise ValueError('\'{0}\' must be a scalar'.format(k))

            if not utils.isreal(arg).all():
                raise ValueError('\'{0}\' must be real'.format(k))

    num = np.abs(num)
    dat = np.asanyarray(dat, dtype=np.float)
    ll = len(x)
    if ll < num:
        # Warn - 'num' is greater than len(x), returning zeros
        out = np.zeros(dat.size, dtype='int8')
    else:
        out = stuckvalues(dat, reso, num)

    return out

def dataqc_gradienttest_wrapper(dat, x, ddatdx, mindx, startdat, toldat, strict_validation=False):
    if is_none(ddatdx) or is_fill(ddatdx) or is_none(mindx) or is_fill(mindx) or is_none(startdat) or is_fill(startdat) or is_none(toldat) or is_fill(toldat):
        out = np.empty(dat.shape, dtype=np.int8)
        out.fill(-99)
        return out
    outqc = dataqc_gradienttest(dat, x, [-np.atleast_1d(ddatdx)[-1], np.atleast_1d(ddatdx)[-1]], np.atleast_1d(mindx)[-1], np.atleast_1d(startdat)[-1], np.atleast_1d(toldat)[-1], strict_validation=strict_validation)
    return outqc


def dataqc_gradienttest(dat, x, ddatdx, mindx, startdat, toldat, strict_validation=False):
    """
    Description

        Data quality control algorithm testing if changes between successive
        data points fall within a certain range.

        Input data dat are given as a function of coordinate x. The algorithm
        will flag dat values as bad if the change deltaDAT/deltaX between
        successive dat values exceeds thresholds given in ddatdx. Once the
        threshold is exceeded, following dat are considered bad until a dat
        value returns to within toldat of the last known good value.

        It is possible to remove data points that are too close together in x
        coordinates (use mindx).

        By default, the first value of dat is considered good. To change this,
        use startdat and toldat to set as the first good data point the first
        one that comes within toldat of startdat.

    Implemented by:

        2012-07-17: DPS authored by Mathias Lankhorst. Example code provided
        for Matlab.
        2013-04-06: Christopher Wingard. Initial python implementation.

    Usage:

        outdat, outx, outqc = dataqc_gradienttest(dat, x, ddatdx, mindx,
                                                  startdat, toldat);

            where

        outdat = same as dat except that NaNs and values not meeting mindx are
            removed.
        outx = same as x except that NaNs and values not meeting mindx are
            removed.
        outqc = output quality control flags for outdat. 0 means bad data, 1
            means good data.

        dat = input dataset, a numeric real vector.
        x = coordinate (e.g. time, distance) along which dat is given. Must be
            of the same size as dat and strictly increasing.
        ddatdx = two-element vector defining the valid range of ddat/dx
            from one point to the next.
        mindx = scalar. minimum dx for which this test will be applied (data
            that are less than mindx apart will be deleted). defaults to zero
            if NaN/empty.
        startdat = start value (scalar) of dat that is presumed good. defaults
            to first non-NaN value of dat if NaN/empty.
        toldat = tolerance value (scalar) for dat; threshold to within which
            dat must return to be counted as good, after exceeding a ddatdx
            threshold detected bad data.

    References:

        OOI (2012). Data Product Specification for Gradient Test. Document
            Control Number 1341-100010.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI
            >> Controlled >> 1000 System Level >>
            1341-10010_Data_Product_SPEC_GRDTEST_OOI.pdf)
    """

    if strict_validation:
        if not utils.isvector(dat) or not utils.isvector(x):
            raise ValueError('\'dat\' and \'x\' must be vectors')

        if len(dat) != len(x):
            raise ValueError('\'dat\' and \'x\' must be of equal len')

        if not all(np.diff(x) > 0):
            raise ValueError('\'x\' must be montonically increasing')

    dat = np.asanyarray(dat, dtype=np.float).flatten()
    x = np.asanyarray(x, dtype=np.float).flatten()

    if np.isnan(mindx):
        mindx = 0
    mindx = mindx or 0
    if np.isnan(startdat):
        startdat = 0
    startdat = startdat or 0
    
    # No strict validation here, they are scalards and they must be validated
    # before going into the C-layer
    if not utils.isscalar(mindx):
        raise ValueError("'mindx' must be scalar, NaN, or empty.")
    if not utils.isscalar(startdat):
        raise ValueError("'startdat' must be scalar, NaN, or empty.")


    # Confirm that there are still data points left, else abort:
    if np.abs(x[0] - x[-1]) < mindx:
        out = np.zeros(x.shape)
        out.fill(1)
        log.warn('Too few values to inspect')
        return out


    grad_min = ddatdx[0]
    grad_max = ddatdx[1]
    out = gradientvalues(dat, x, grad_min, grad_max, mindx, startdat, toldat)

    return out


def dataqc_solarelevation(lon, lat, dt):
    """
    Description

        Computes instantaneous no-sky solar radiation and altitude from date
        and time stamp and position data. It is put together from expressions
        taken from Appendix E in the 1978 edition of Almanac for Computers,
        Nautical Almanac Office, U.S. Naval Observatory. They are reduced
        accuracy expressions valid for the years 1800-2100. Solar declination
        computed from these expressions is accurate to at least 1'. The solar
        constant (1368.0 W/m^2) represents a mean of satellite measurements
        made over the last sunspot cycle (1979-1995) taken from Coffey et al
        (1995), Earth System Monitor, 6, 6-10.

        This code is a python implementation of soradna1.m available in Air-Sea
        Toolbox.

    Implemented by:

        1997-03-08: Version 1.0 (author unknown) of soradna1.m.
        1998-08-28: Version 1.1 (author unknown) of soradna1.m.
        1999-08-05: Version 2.0 (author unknown) of soradna1.m.

        2013-04-07: Christopher Wingard. Initial python implementation. Note,
        this function is derived from old, unmaintained code. More robust
        implementations exist (e.g. PyEphem and PySolar) that will probably
        calculate these values more accurately.

    Usage:

        z, sorad = dataqc_solarelevation(lon, lat, dt)

            where

        z = solar altitude [degrees]
        sorad = no atmosphere solar radiation [W m^-2]

        lon = longitude (east is positive) [decimal degress]
        lat = latitude [decimal degrees]
        dt = date and time stamp in UTC [seconds since 1970-01-01]

    Examples

        dt = 1329177600     # 2012-02-14 00:00:00
        z, sorad = dataqc_solarelevation(120, 30, dt)
        z = 15.1566, sorad = 366.8129

        OOI (2012). Data Product Specification for Solar Elevation. Document
            Control Number 1341-100011.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI
            >> Controlled >> 1000 System Level >>
            1341-10011_Data_Product_SPEC_SOLRELV_OOI.pdf)
    """
    # Test lengths and types of inputs. Latitude and longitude must be the same
    # size and can either be a scalar or a vecotr. The date and time stamp
    # can also be either a scalar or a vector. If all three inputs are vectors,
    # they must be of the same length.
    if len(lon) != len(lat):
        raise ValueError('\'lon\' and \'lat\' must be the same size')

    if utils.isvector(lon) and utils.isvector(lat) and utils.isvector(dt):
        # test their lengths
        if not len(lon) == len(lat) == len(dt):
            raise ValueError('If all inputs are vectors, these must all '
                             'be of the same length')

    # set constants (using values from as_consts.m)
    # ------ short-wave flux calculations
    #   the solar constant [W m^-2] represents a mean of satellite measurements
    #   made over the last sunspot cycle (1979-1995), taken from Coffey et al.
    #   (1995), Earth System Monitor, 6, 6-10.
    solar_const = 1368.0

    # Create a time tuple in UTC from the Epoch time input, and then create
    # scalars or numpy arrays of time elements for subsequent calculations.
    ldt = len(dt)
    yy = np.zeros(ldt, dtype=np.int)
    mn = np.zeros(ldt, dtype=np.int)
    dd = np.zeros(ldt, dtype=np.int)
    hh = np.zeros(ldt, dtype=np.int)
    mm = np.zeros(ldt, dtype=np.int)
    ss = np.zeros(ldt, dtype=np.int)
    for i in range(ldt):
        # create time tuple in UTC
        gtime = time.gmtime(dt[i])
        # create scalar elements
        yy[i] = gtime[0]
        mn[i] = gtime[1]
        dd[i] = gtime[2]
        hh[i] = gtime[3]
        mm[i] = gtime[4]
        ss[i] = gtime[5]

    #constants used in function
    deg2rad = np.pi / 180.0
    rad2deg = 1 / deg2rad

    # compute Universal Time in hours
    utime = hh + (mm + ss / 60.0) / 60.0

    # compute Julian ephemeris date in days (Day 1 is 1 Jan 4713 B.C. which
    # equals -4712 Jan 1)
    jed = (367.0 * yy - np.fix(7.0*(yy+np.fix((mn+9)/12.0))/4.0)
           + np.fix(275.0*mn/9.0) + dd + 1721013 + utime / 24.0)

    # compute interval in Julian centuries since 1900
    jc_int = (jed - 2415020.0) / 36525.0

    # compute mean anomaly of the sun
    ma_sun = 358.475833 + 35999.049750 * jc_int - 0.000150 * jc_int**2
    ma_sun = (ma_sun - np.fix(ma_sun/360.0) * 360.0) * deg2rad

    # compute mean longitude of sun
    ml_sun = 279.696678 + 36000.768920 * jc_int + 0.000303 * jc_int**2
    ml_sun = (ml_sun - np.fix(ml_sun/360.0) * 360.0) * deg2rad

    # compute mean anomaly of Jupiter
    ma_jup = 225.444651 + 2880.0 * jc_int + 154.906654 * jc_int
    ma_jup = (ma_jup - np.fix(ma_jup/360.0) * 360.0) * deg2rad

    # compute longitude of the ascending node of the moon's orbit
    an_moon = (259.183275 - 1800 * jc_int - 134.142008 * jc_int
               + 0.002078 * jc_int**2)
    an_moon = (an_moon - np.fix(an_moon/360.0) * 360.0 + 360.0) * deg2rad

    # compute mean anomaly of Venus
    ma_ven = (212.603219 + 58320 * jc_int + 197.803875 * jc_int
              + 0.001286 * jc_int**2)
    ma_ven = (ma_ven - np.fix(ma_ven/360.0) * 360.0) * deg2rad

    # compute sun theta
    theta = (0.397930 * np.sin(ml_sun) + 0.009999 * np.sin(ma_sun-ml_sun)
             + 0.003334 * np.sin(ma_sun+ml_sun) - 0.000208 * jc_int
             * np.sin(ml_sun) + 0.000042 * np.sin(2*ma_sun+ml_sun) - 0.000040
             * np.cos(ml_sun) - 0.000039 * np.sin(an_moon-ml_sun) - 0.000030
             * jc_int * np.sin(ma_sun-ml_sun) - 0.000014
             * np.sin(2*ma_sun-ml_sun) - 0.000010
             * np.cos(ma_sun-ml_sun-ma_jup) - 0.000010 * jc_int
             * np.sin(ma_sun+ml_sun))

    # compute sun rho
    rho = (1.000421 - 0.033503 * np.cos(ma_sun) - 0.000140 * np.cos(2*ma_sun)
           + 0.000084 * jc_int * np.cos(ma_sun) - 0.000033
           * np.sin(ma_sun-ma_jup) + 0.000027 * np.sin(2.*ma_sun-2.*ma_ven))

    # compute declination
    decln = np.arcsin(theta/np.sqrt(rho))

    # compute equation of time (in seconds of time)
    l = 276.697 + 0.98564734 * (jed-2415020.0)
    l = (l - 360.0 * np.fix(l/360.0)) * deg2rad
    eqt = (-97.8 * np.sin(l) - 431.3 * np.cos(l) + 596.6 * np.sin(2*l)
           - 1.9 * np.cos(2*l) + 4.0 * np.sin(3*l) + 19.3 * np.cos(3*l)
           - 12.7 * np.sin(4*l))
    eqt = eqt / 60.0

    # compute local hour angle from global hour angle
    gha = 15.0 * (utime-12) + 15.0 * eqt / 60.0
    lha = gha - lon

    # compute radius vector
    rv = np.sqrt(rho)

    # compute solar altitude
    sz = (np.sin(deg2rad*lat) * np.sin(decln) + np.cos(deg2rad*lat)
          * np.cos(decln) * np.cos(deg2rad*lha))
    z = rad2deg * np.arcsin(sz)

    # compute solar radiation outside atmosphere (defaults to 0 when solar
    # altitude is below the horizon)
    sorad = (solar_const / rv**2) * np.sin(deg2rad * z)
    sorad[z < 0] = 0

    return (z, sorad)

def dataqc_propagateflags_wrapper(strict_validation=False, *args):
    '''
    This is a function that wraps dataqc_propagateflags for use in ION
    It accepts a variable number of vector arguments (of the same shape) and calls dataqc_propagateflags
    '''

    if not strict_validation:

        shapes = np.array([i.shape[0] for i in args])
        if not (shapes == shapes[0]).all():
            raise ValueError('Input vectors are not the same shape')

    return dataqc_propagateflags(np.array(args), strict_validation=strict_validation)



def dataqc_propagateflags(inflags, strict_validation=False):
    """
    Description:

        Propagate "bad" qc flags (from an arbitrary number of source datasets)
        to another (derived) dataset.

        Consider data from an oceanographic CTD (conductivity, temperature, and
        pressure) instrument. From these three time series, you want to compute
        salinity. If any of the three source data (conductivity, temperature,
        pressure) is of bad quality, the salinity will be bad as well. You can
        feed your QC assessment of the former three into this routine, which
        will then give you the combined assessment for the derived (here:
        salinity) property.

    Implemented by:

        2012-07-17: DPS authored by Mathias Lankhorst. Example code provided
        for Matlab.

        2013-04-06: Christopher Wingard. Initial python implementation.

    Usage:

        outflag = dataqc_propagateflags(inflags)

            where

        outflag = a 1-by-N boolean vector that contains 1 where all of the
            inflags are 1, and 0 otherwise.

        inflags = an M-by-N boolean matrix, where each of the M rows contains
            flags of an independent data set such that "0" means bad data and
            "1" means good data.

    References:

        OOI (2012). Data Product Specification for Combined QC Flags. Document
            Control Number 1341-100012.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI
            >> Controlled >> 1000 System Level >>
            1341-10012_Data_Product_SPEC_CMBNFLG_OOI.pdf)
    """

    if strict_validation:
        if not utils.islogical(inflags):
            raise ValueError('\'inflags\' must be \'0\' or \'1\' '
                             'integer flag array')

    array_size = inflags.shape
    nrows = array_size[0]
    if nrows < 2:
        error('\'inflags\' must be at least a two-dimensional array')

    outflag = np.all(inflags, 0)
    return outflag.astype('int8')


def dataqc_condcompress(p_orig, p_new, c_orig, cpcor=-9.57e-8):
    """
    Description:

        Implementation of the Sea-Bird conductivity compressibility correction,
        scaling the input conductivity based on ratio of the original pressure
        and the updated pressure.

    Implemented by:

        2013-04-07: Christopher Wingard. Initial python implementation.

    Usage:

        c_new = dataqc_condcompress(p_orig, p_new, c_orig, cpcor)

            where

        c_new = updated conductivity record [S/m]

        p_orig = original pressure used to calculate original conductivity,
            this typically the L1a PRESWAT [dbar]
        p_new = updated pressure, typically L1b PRESWAT [dbar]
        c_orig = original conductivty record, typically L1a CONDWAT [S/m]
        cpcor = pressure correction coefficient used to calculate original
            conductivity, default is -9.57e-8

    References:

        OOI (2012). Data Product Specification for Conductivity Compressibility
            Correction. Document Control Number 1341-10030.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI
            >> Controlled >> 1000 System Level >>
            1341-10030_Data_Product_SPEC_CNDCMPR_OOI.pdf)
    """
    c_new = c_orig * (1 + cpcor * p_orig) / (1 + cpcor * p_new)
    return c_new

[((80, 10, 80, 28), 'numpy.atleast_1d', 'np.atleast_1d', ({(80, 24, 80, 27): 'dat'}, {}), '(dat)', True, 'import numpy as np\n'), ((81, 13, 81, 34), 'numpy.atleast_1d', 'np.atleast_1d', ({(81, 27, 81, 33): 'datlim'}, {}), '(datlim)', True, 'import numpy as np\n'), ((262, 15, 262, 35), 'numpy.max', 'np.max', ({(262, 22, 262, 34): 'datlim[:, (1)]'}, {}), '(datlim[:, (1)])', True, 'import numpy as np\n'), ((263, 15, 263, 35), 'numpy.min', 'np.min', ({(263, 22, 263, 34): 'datlim[:, (0)]'}, {}), '(datlim[:, (0)])', True, 'import numpy as np\n'), ((328, 10, 328, 28), 'numpy.atleast_1d', 'np.atleast_1d', ({(328, 24, 328, 27): 'dat'}, {}), '(dat)', True, 'import numpy as np\n'), ((346, 10, 346, 44), 'numpy.asanyarray', 'np.asanyarray', (), '', True, 'import numpy as np\n'), ((348, 10, 348, 37), 'ion_functions.qc.qc_extensions.spikevalues', 'spikevalues', ({(348, 22, 348, 25): 'dat', (348, 27, 348, 28): 'L', (348, 30, 348, 31): 'N', (348, 33, 348, 36): 'acc'}, {}), '(dat, L, N, acc)', False, 'from ion_functions.qc.qc_extensions import stuckvalues, spikevalues, gradientvalues, ntp_to_month\n'), ((407, 10, 407, 28), 'numpy.atleast_1d', 'np.atleast_1d', ({(407, 24, 407, 27): 'dat'}, {}), '(dat)', True, 'import numpy as np\n'), ((408, 8, 408, 24), 'numpy.atleast_1d', 'np.atleast_1d', ({(408, 22, 408, 23): 't'}, {}), '(t)', True, 'import numpy as np\n'), ((433, 9, 433, 34), 'numpy.polyfit', 'np.polyfit', ({(433, 20, 433, 21): 't', (433, 23, 433, 26): 'dat', (433, 28, 433, 33): 'ord_n'}, {}), '(t, dat, ord_n)', True, 'import numpy as np\n'), ((434, 12, 434, 29), 'numpy.polyval', 'np.polyval', ({(434, 23, 434, 25): 'pp', (434, 27, 434, 28): 't'}, {}), '(pp, t)', True, 'import numpy as np\n'), ((489, 10, 489, 26), 'numpy.atleast_1d', 'np.atleast_1d', ({(489, 24, 489, 25): 'x'}, {}), '(x)', True, 'import numpy as np\n'), ((511, 10, 511, 21), 'numpy.abs', 'np.abs', ({(511, 17, 511, 20): 'num'}, {}), '(num)', True, 'import numpy as np\n'), ((512, 10, 512, 44), 'numpy.asanyarray', 'np.asanyarray', (), '', True, 'import numpy as np\n'), ((607, 7, 607, 22), 'numpy.isnan', 'np.isnan', ({(607, 16, 607, 21): 'mindx'}, {}), '(mindx)', True, 'import numpy as np\n'), ((610, 7, 610, 25), 'numpy.isnan', 'np.isnan', ({(610, 16, 610, 24): 'startdat'}, {}), '(startdat)', True, 'import numpy as np\n'), ((632, 10, 632, 77), 'ion_functions.qc.qc_extensions.gradientvalues', 'gradientvalues', ({(632, 25, 632, 28): 'dat', (632, 30, 632, 31): 'x', (632, 33, 632, 41): 'grad_min', (632, 43, 632, 51): 'grad_max', (632, 53, 632, 58): 'mindx', (632, 60, 632, 68): 'startdat', (632, 70, 632, 76): 'toldat'}, {}), '(dat, x, grad_min, grad_max, mindx, startdat, toldat)', False, 'from ion_functions.qc.qc_extensions import stuckvalues, spikevalues, gradientvalues, ntp_to_month\n'), ((713, 9, 713, 36), 'numpy.zeros', 'np.zeros', (), '', True, 'import numpy as np\n'), ((714, 9, 714, 36), 'numpy.zeros', 'np.zeros', (), '', True, 'import numpy as np\n'), ((715, 9, 715, 36), 'numpy.zeros', 'np.zeros', (), '', True, 'import numpy as np\n'), ((716, 9, 716, 36), 'numpy.zeros', 'np.zeros', (), '', True, 'import numpy as np\n'), ((717, 9, 717, 36), 'numpy.zeros', 'np.zeros', (), '', True, 'import numpy as np\n'), ((718, 9, 718, 36), 'numpy.zeros', 'np.zeros', (), '', True, 'import numpy as np\n'), ((798, 9, 798, 21), 'numpy.sqrt', 'np.sqrt', ({(798, 17, 798, 20): 'rho'}, {}), '(rho)', True, 'import numpy as np\n'), ((882, 14, 882, 32), 'numpy.all', 'np.all', ({(882, 21, 882, 28): 'inflags', (882, 30, 882, 31): '0'}, {}), '(inflags, 0)', True, 'import numpy as np\n'), ((24, 10, 24, 44), 'logging.getLogger', 'logging.getLogger', ({(24, 28, 24, 43): '"""ion-functions"""'}, {}), "('ion-functions')", False, 'import logging\n'), ((37, 14, 37, 48), 'numpy.empty', 'np.empty', (), '', True, 'import numpy as np\n'), ((103, 14, 103, 48), 'numpy.empty', 'np.empty', (), '', True, 'import numpy as np\n'), ((107, 14, 107, 48), 'numpy.empty', 'np.empty', (), '', True, 'import numpy as np\n'), ((111, 14, 111, 48), 'numpy.empty', 'np.empty', (), '', True, 'import numpy as np\n'), ((115, 14, 115, 48), 'numpy.empty', 'np.empty', (), '', True, 'import numpy as np\n'), ((132, 12, 132, 30), 'numpy.column_stack', 'np.column_stack', ({(132, 28, 132, 29): 'z'}, {}), '(z)', True, 'import numpy as np\n'), ((246, 15, 246, 77), 'numpy.interp', 'np.interp', (), '', True, 'import numpy as np\n'), ((248, 15, 248, 77), 'numpy.interp', 'np.interp', (), '', True, 'import numpy as np\n'), ((261, 10, 261, 24), 'numpy.isnan', 'np.isnan', ({(261, 19, 261, 23): 'lim1'}, {}), '(lim1)', True, 'import numpy as np\n'), ((261, 29, 261, 43), 'numpy.isnan', 'np.isnan', ({(261, 38, 261, 42): 'lim2'}, {}), '(lim2)', True, 'import numpy as np\n'), ((272, 14, 272, 48), 'numpy.empty', 'np.empty', (), '', True, 'import numpy as np\n'), ((354, 14, 354, 48), 'numpy.empty', 'np.empty', (), '', True, 'import numpy as np\n'), ((449, 14, 449, 40), 'numpy.empty', 'np.empty', ({(449, 23, 449, 30): 'x.shape', (449, 32, 449, 39): 'np.int8'}, {}), '(x.shape, np.int8)', True, 'import numpy as np\n'), ((516, 14, 516, 46), 'numpy.zeros', 'np.zeros', (), '', True, 'import numpy as np\n'), ((518, 14, 518, 41), 'ion_functions.qc.qc_extensions.stuckvalues', 'stuckvalues', ({(518, 26, 518, 29): 'dat', (518, 31, 518, 35): 'reso', (518, 37, 518, 40): 'num'}, {}), '(dat, reso, num)', False, 'from ion_functions.qc.qc_extensions import stuckvalues, spikevalues, gradientvalues, ntp_to_month\n'), ((524, 14, 524, 48), 'numpy.empty', 'np.empty', (), '', True, 'import numpy as np\n'), ((616, 11, 616, 32), 'ion_functions.utils.isscalar', 'utils.isscalar', ({(616, 26, 616, 31): 'mindx'}, {}), '(mindx)', False, 'from ion_functions import utils\n'), ((618, 11, 618, 35), 'ion_functions.utils.isscalar', 'utils.isscalar', ({(618, 26, 618, 34): 'startdat'}, {}), '(startdat)', False, 'from ion_functions import utils\n'), ((623, 7, 623, 27), 'numpy.abs', 'np.abs', ({(623, 14, 623, 26): '(x[0] - x[-1])'}, {}), '(x[0] - x[-1])', True, 'import numpy as np\n'), ((624, 14, 624, 31), 'numpy.zeros', 'np.zeros', ({(624, 23, 624, 30): 'x.shape'}, {}), '(x.shape)', True, 'import numpy as np\n'), ((697, 7, 697, 26), 'ion_functions.utils.isvector', 'utils.isvector', ({(697, 22, 697, 25): 'lon'}, {}), '(lon)', False, 'from ion_functions import utils\n'), ((697, 31, 697, 50), 'ion_functions.utils.isvector', 'utils.isvector', ({(697, 46, 697, 49): 'lat'}, {}), '(lat)', False, 'from ion_functions import utils\n'), ((697, 55, 697, 73), 'ion_functions.utils.isvector', 'utils.isvector', ({(697, 70, 697, 72): 'dt'}, {}), '(dt)', False, 'from ion_functions import utils\n'), ((721, 16, 721, 34), 'time.gmtime', 'time.gmtime', ({(721, 28, 721, 33): 'dt[i]'}, {}), '(dt[i])', False, 'import time\n'), ((803, 18, 803, 31), 'numpy.arcsin', 'np.arcsin', ({(803, 28, 803, 30): 'sz'}, {}), '(sz)', True, 'import numpy as np\n'), ((807, 36, 807, 55), 'numpy.sin', 'np.sin', ({(807, 43, 807, 54): '(deg2rad * z)'}, {}), '(deg2rad * z)', True, 'import numpy as np\n'), ((820, 17, 820, 53), 'numpy.array', 'np.array', ({(820, 26, 820, 52): '[i.shape[0] for i in args]'}, {}), '([i.shape[0] for i in args])', True, 'import numpy as np\n'), ((824, 33, 824, 47), 'numpy.array', 'np.array', ({(824, 42, 824, 46): 'args'}, {}), '(args)', True, 'import numpy as np\n'), ((27, 11, 27, 29), 'numpy.atleast_1d', 'np.atleast_1d', ({(27, 25, 27, 28): 'arr'}, {}), '(arr)', True, 'import numpy as np\n'), ((124, 16, 124, 48), 'numpy.asanyarray', 'np.asanyarray', (), '', True, 'import numpy as np\n'), ((125, 16, 125, 31), 'ion_functions.qc.qc_extensions.ntp_to_month', 'ntp_to_month', ({(125, 29, 125, 30): 'v'}, {}), '(v)', False, 'from ion_functions.qc.qc_extensions import stuckvalues, spikevalues, gradientvalues, ntp_to_month\n'), ((179, 15, 179, 34), 'ion_functions.utils.isvector', 'utils.isvector', ({(179, 30, 179, 33): 'dat'}, {}), '(dat)', False, 'from ion_functions import utils\n'), ((182, 15, 182, 37), 'ion_functions.utils.ismatrix', 'utils.ismatrix', ({(182, 30, 182, 36): 'datlim'}, {}), '(datlim)', False, 'from ion_functions import utils\n'), ((275, 33, 275, 51), 'numpy.atleast_1d', 'np.atleast_1d', ({(275, 47, 275, 50): 'acc'}, {}), '(acc)', True, 'import numpy as np\n'), ((275, 57, 275, 73), 'numpy.atleast_1d', 'np.atleast_1d', ({(275, 71, 275, 72): 'N'}, {}), '(N)', True, 'import numpy as np\n'), ((275, 79, 275, 95), 'numpy.atleast_1d', 'np.atleast_1d', ({(275, 93, 275, 94): 'L'}, {}), '(L)', True, 'import numpy as np\n'), ((337, 15, 337, 34), 'ion_functions.utils.isvector', 'utils.isvector', ({(337, 30, 337, 33): 'dat'}, {}), '(dat)', False, 'from ion_functions import utils\n'), ((357, 40, 357, 60), 'numpy.atleast_1d', 'np.atleast_1d', ({(357, 54, 357, 59): 'ord_n'}, {}), '(ord_n)', True, 'import numpy as np\n'), ((357, 66, 357, 85), 'numpy.atleast_1d', 'np.atleast_1d', ({(357, 80, 357, 84): 'nstd'}, {}), '(nstd)', True, 'import numpy as np\n'), ((452, 37, 452, 56), 'numpy.atleast_1d', 'np.atleast_1d', ({(452, 51, 452, 55): 'reso'}, {}), '(reso)', True, 'import numpy as np\n'), ((452, 62, 452, 80), 'numpy.atleast_1d', 'np.atleast_1d', ({(452, 76, 452, 79): 'num'}, {}), '(num)', True, 'import numpy as np\n'), ((495, 15, 495, 34), 'ion_functions.utils.isvector', 'utils.isvector', ({(495, 30, 495, 33): 'dat'}, {}), '(dat)', False, 'from ion_functions import utils\n'), ((527, 97, 527, 117), 'numpy.atleast_1d', 'np.atleast_1d', ({(527, 111, 527, 116): 'mindx'}, {}), '(mindx)', True, 'import numpy as np\n'), ((527, 123, 527, 146), 'numpy.atleast_1d', 'np.atleast_1d', ({(527, 137, 527, 145): 'startdat'}, {}), '(startdat)', True, 'import numpy as np\n'), ((527, 152, 527, 173), 'numpy.atleast_1d', 'np.atleast_1d', ({(527, 166, 527, 172): 'toldat'}, {}), '(toldat)', True, 'import numpy as np\n'), ((604, 10, 604, 44), 'numpy.asanyarray', 'np.asanyarray', (), '', True, 'import numpy as np\n'), ((605, 8, 605, 40), 'numpy.asanyarray', 'np.asanyarray', (), '', True, 'import numpy as np\n'), ((775, 15, 775, 36), 'numpy.sin', 'np.sin', ({(775, 22, 775, 35): '(ma_sun + ml_sun)'}, {}), '(ma_sun + ml_sun)', True, 'import numpy as np\n'), ((780, 48, 780, 75), 'numpy.sin', 'np.sin', ({(780, 55, 780, 74): '(2.0 * ma_sun - 2.0 * ma_ven)'}, {}), '(2.0 * ma_sun - 2.0 * ma_ven)', True, 'import numpy as np\n'), ((783, 28, 783, 40), 'numpy.sqrt', 'np.sqrt', ({(783, 36, 783, 39): 'rho'}, {}), '(rho)', True, 'import numpy as np\n'), ((790, 20, 790, 31), 'numpy.sin', 'np.sin', ({(790, 27, 790, 30): '(4 * l)'}, {}), '(4 * l)', True, 'import numpy as np\n'), ((801, 10, 801, 29), 'numpy.sin', 'np.sin', ({(801, 17, 801, 28): '(deg2rad * lat)'}, {}), '(deg2rad * lat)', True, 'import numpy as np\n'), ((801, 32, 801, 45), 'numpy.sin', 'np.sin', ({(801, 39, 801, 44): 'decln'}, {}), '(decln)', True, 'import numpy as np\n'), ((802, 28, 802, 47), 'numpy.cos', 'np.cos', ({(802, 35, 802, 46): '(deg2rad * lha)'}, {}), '(deg2rad * lha)', True, 'import numpy as np\n'), ((873, 15, 873, 39), 'ion_functions.utils.islogical', 'utils.islogical', ({(873, 31, 873, 38): 'inflags'}, {}), '(inflags)', False, 'from ion_functions import utils\n'), ((29, 27, 29, 45), 'numpy.atleast_1d', 'np.atleast_1d', ({(29, 41, 29, 44): 'arr'}, {}), '(arr)', True, 'import numpy as np\n'), ((40, 40, 40, 62), 'numpy.atleast_1d', 'np.atleast_1d', ({(40, 54, 40, 61): 'dat_min'}, {}), '(dat_min)', True, 'import numpy as np\n'), ((40, 68, 40, 90), 'numpy.atleast_1d', 'np.atleast_1d', ({(40, 82, 40, 89): 'dat_max'}, {}), '(dat_max)', True, 'import numpy as np\n'), ((419, 19, 419, 38), 'ion_functions.utils.isvector', 'utils.isvector', ({(419, 34, 419, 37): 'arg'}, {}), '(arg)', False, 'from ion_functions import utils\n'), ((423, 19, 423, 38), 'ion_functions.utils.isscalar', 'utils.isscalar', ({(423, 34, 423, 37): 'arg'}, {}), '(arg)', False, 'from ion_functions import utils\n'), ((443, 13, 443, 31), 'numpy.ones', 'np.ones', ({(443, 21, 443, 30): 'dat.shape'}, {}), '(dat.shape)', True, 'import numpy as np\n'), ((505, 19, 505, 38), 'ion_functions.utils.isscalar', 'utils.isscalar', ({(505, 34, 505, 37): 'arg'}, {}), '(arg)', False, 'from ion_functions import utils\n'), ((527, 69, 527, 90), 'numpy.atleast_1d', 'np.atleast_1d', ({(527, 83, 527, 89): 'ddatdx'}, {}), '(ddatdx)', True, 'import numpy as np\n'), ((595, 15, 595, 34), 'ion_functions.utils.isvector', 'utils.isvector', ({(595, 30, 595, 33): 'dat'}, {}), '(dat)', False, 'from ion_functions import utils\n'), ((595, 42, 595, 59), 'ion_functions.utils.isvector', 'utils.isvector', ({(595, 57, 595, 58): 'x'}, {}), '(x)', False, 'from ion_functions import utils\n'), ((747, 23, 747, 43), 'numpy.fix', 'np.fix', ({(747, 30, 747, 42): '(ma_sun / 360.0)'}, {}), '(ma_sun / 360.0)', True, 'import numpy as np\n'), ((751, 23, 751, 43), 'numpy.fix', 'np.fix', ({(751, 30, 751, 42): '(ml_sun / 360.0)'}, {}), '(ml_sun / 360.0)', True, 'import numpy as np\n'), ((755, 23, 755, 43), 'numpy.fix', 'np.fix', ({(755, 30, 755, 42): '(ma_jup / 360.0)'}, {}), '(ma_jup / 360.0)', True, 'import numpy as np\n'), ((765, 23, 765, 43), 'numpy.fix', 'np.fix', ({(765, 30, 765, 42): '(ma_ven / 360.0)'}, {}), '(ma_ven / 360.0)', True, 'import numpy as np\n'), ((774, 15, 774, 43), 'numpy.cos', 'np.cos', ({(774, 22, 774, 42): '(ma_sun - ml_sun - ma_jup)'}, {}), '(ma_sun - ml_sun - ma_jup)', True, 'import numpy as np\n'), ((780, 13, 780, 34), 'numpy.sin', 'np.sin', ({(780, 20, 780, 33): '(ma_sun - ma_jup)'}, {}), '(ma_sun - ma_jup)', True, 'import numpy as np\n'), ((787, 21, 787, 36), 'numpy.fix', 'np.fix', ({(787, 28, 787, 35): '(l / 360.0)'}, {}), '(l / 360.0)', True, 'import numpy as np\n'), ((789, 60, 789, 71), 'numpy.cos', 'np.cos', ({(789, 67, 789, 70): '(3 * l)'}, {}), '(3 * l)', True, 'import numpy as np\n'), ((801, 48, 801, 67), 'numpy.cos', 'np.cos', ({(801, 55, 801, 66): '(deg2rad * lat)'}, {}), '(deg2rad * lat)', True, 'import numpy as np\n'), ((802, 12, 802, 25), 'numpy.cos', 'np.cos', ({(802, 19, 802, 24): 'decln'}, {}), '(decln)', True, 'import numpy as np\n'), ((84, 15, 84, 35), 'ion_functions.utils.isnumeric', 'utils.isnumeric', ({(84, 31, 84, 34): 'dat'}, {}), '(dat)', False, 'from ion_functions import utils\n'), ((87, 15, 87, 32), 'ion_functions.utils.isreal', 'utils.isreal', ({(87, 28, 87, 31): 'dat'}, {}), '(dat)', False, 'from ion_functions import utils\n'), ((90, 15, 90, 38), 'ion_functions.utils.isnumeric', 'utils.isnumeric', ({(90, 31, 90, 37): 'datlim'}, {}), '(datlim)', False, 'from ion_functions import utils\n'), ((93, 15, 93, 35), 'ion_functions.utils.isreal', 'utils.isreal', ({(93, 28, 93, 34): 'datlim'}, {}), '(datlim)', False, 'from ion_functions import utils\n'), ((331, 15, 331, 35), 'ion_functions.utils.isnumeric', 'utils.isnumeric', ({(331, 31, 331, 34): 'dat'}, {}), '(dat)', False, 'from ion_functions import utils\n'), ((334, 15, 334, 32), 'ion_functions.utils.isreal', 'utils.isreal', ({(334, 28, 334, 31): 'dat'}, {}), '(dat)', False, 'from ion_functions import utils\n'), ((437, 52, 437, 63), 'numpy.std', 'np.std', ({(437, 59, 437, 62): 'dat'}, {}), '(dat)', True, 'import numpy as np\n'), ((492, 15, 492, 35), 'ion_functions.utils.isnumeric', 'utils.isnumeric', ({(492, 31, 492, 34): 'dat'}, {}), '(dat)', False, 'from ion_functions import utils\n'), ((498, 15, 498, 32), 'ion_functions.utils.isreal', 'utils.isreal', ({(498, 28, 498, 31): 'dat'}, {}), '(dat)', False, 'from ion_functions import utils\n'), ((527, 42, 527, 63), 'numpy.atleast_1d', 'np.atleast_1d', ({(527, 56, 527, 62): 'ddatdx'}, {}), '(ddatdx)', True, 'import numpy as np\n'), ((601, 19, 601, 29), 'numpy.diff', 'np.diff', ({(601, 27, 601, 28): 'x'}, {}), '(x)', True, 'import numpy as np\n'), ((740, 13, 740, 33), 'numpy.fix', 'np.fix', ({(740, 20, 740, 32): '(275.0 * mn / 9.0)'}, {}), '(275.0 * mn / 9.0)', True, 'import numpy as np\n'), ((760, 25, 760, 46), 'numpy.fix', 'np.fix', ({(760, 32, 760, 45): '(an_moon / 360.0)'}, {}), '(an_moon / 360.0)', True, 'import numpy as np\n'), ((773, 15, 773, 38), 'numpy.sin', 'np.sin', ({(773, 22, 773, 37): '(2 * ma_sun - ml_sun)'}, {}), '(2 * ma_sun - ml_sun)', True, 'import numpy as np\n'), ((779, 33, 779, 47), 'numpy.cos', 'np.cos', ({(779, 40, 779, 46): 'ma_sun'}, {}), '(ma_sun)', True, 'import numpy as np\n'), ((789, 39, 789, 50), 'numpy.sin', 'np.sin', ({(789, 46, 789, 49): '(3 * l)'}, {}), '(3 * l)', True, 'import numpy as np\n'), ((102, 33, 102, 54), 'numpy.atleast_1d', 'np.atleast_1d', ({(102, 47, 102, 53): 'datlim'}, {}), '(datlim)', True, 'import numpy as np\n'), ((106, 34, 106, 55), 'numpy.atleast_1d', 'np.atleast_1d', ({(106, 48, 106, 54): 'datlim'}, {}), '(datlim)', True, 'import numpy as np\n'), ((188, 19, 188, 39), 'ion_functions.utils.isnumeric', 'utils.isnumeric', ({(188, 35, 188, 38): 'arg'}, {}), '(arg)', False, 'from ion_functions import utils\n'), ((191, 19, 191, 36), 'ion_functions.utils.isreal', 'utils.isreal', ({(191, 32, 191, 35): 'arg'}, {}), '(arg)', False, 'from ion_functions import utils\n'), ((341, 19, 341, 39), 'ion_functions.utils.isnumeric', 'utils.isnumeric', ({(341, 35, 341, 38): 'arg'}, {}), '(arg)', False, 'from ion_functions import utils\n'), ((344, 19, 344, 36), 'ion_functions.utils.isreal', 'utils.isreal', ({(344, 32, 344, 35): 'arg'}, {}), '(arg)', False, 'from ion_functions import utils\n'), ((412, 19, 412, 39), 'ion_functions.utils.isnumeric', 'utils.isnumeric', ({(412, 35, 412, 38): 'arg'}, {}), '(arg)', False, 'from ion_functions import utils\n'), ((415, 19, 415, 36), 'ion_functions.utils.isreal', 'utils.isreal', ({(415, 32, 415, 35): 'arg'}, {}), '(arg)', False, 'from ion_functions import utils\n'), ((437, 22, 437, 41), 'numpy.std', 'np.std', ({(437, 29, 437, 40): '(dat - datpp)'}, {}), '(dat - datpp)', True, 'import numpy as np\n'), ((502, 19, 502, 39), 'ion_functions.utils.isnumeric', 'utils.isnumeric', ({(502, 35, 502, 38): 'arg'}, {}), '(arg)', False, 'from ion_functions import utils\n'), ((508, 19, 508, 36), 'ion_functions.utils.isreal', 'utils.isreal', ({(508, 32, 508, 35): 'arg'}, {}), '(arg)', False, 'from ion_functions import utils\n'), ((772, 24, 772, 45), 'numpy.sin', 'np.sin', ({(772, 31, 772, 44): '(ma_sun - ml_sun)'}, {}), '(ma_sun - ml_sun)', True, 'import numpy as np\n'), ((778, 61, 778, 77), 'numpy.cos', 'np.cos', ({(778, 68, 778, 76): '(2 * ma_sun)'}, {}), '(2 * ma_sun)', True, 'import numpy as np\n'), ((789, 19, 789, 30), 'numpy.cos', 'np.cos', ({(789, 26, 789, 29): '(2 * l)'}, {}), '(2 * l)', True, 'import numpy as np\n'), ((771, 43, 771, 65), 'numpy.sin', 'np.sin', ({(771, 50, 771, 64): '(an_moon - ml_sun)'}, {}), '(an_moon - ml_sun)', True, 'import numpy as np\n'), ((778, 33, 778, 47), 'numpy.cos', 'np.cos', ({(778, 40, 778, 46): 'ma_sun'}, {}), '(ma_sun)', True, 'import numpy as np\n'), ((788, 59, 788, 70), 'numpy.sin', 'np.sin', ({(788, 66, 788, 69): '(2 * l)'}, {}), '(2 * l)', True, 'import numpy as np\n'), ((771, 15, 771, 29), 'numpy.cos', 'np.cos', ({(771, 22, 771, 28): 'ml_sun'}, {}), '(ml_sun)', True, 'import numpy as np\n'), ((788, 19, 788, 28), 'numpy.sin', 'np.sin', ({(788, 26, 788, 27): 'l'}, {}), '(l)', True, 'import numpy as np\n'), ((788, 39, 788, 48), 'numpy.cos', 'np.cos', ({(788, 46, 788, 47): 'l'}, {}), '(l)', True, 'import numpy as np\n'), ((770, 43, 770, 66), 'numpy.sin', 'np.sin', ({(770, 50, 770, 65): '(2 * ma_sun + ml_sun)'}, {}), '(2 * ma_sun + ml_sun)', True, 'import numpy as np\n'), ((739, 39, 739, 58), 'numpy.fix', 'np.fix', ({(739, 46, 739, 57): '((mn + 9) / 12.0)'}, {}), '((mn + 9) / 12.0)', True, 'import numpy as np\n'), ((770, 15, 770, 29), 'numpy.sin', 'np.sin', ({(770, 22, 770, 28): 'ml_sun'}, {}), '(ml_sun)', True, 'import numpy as np\n'), ((769, 26, 769, 47), 'numpy.sin', 'np.sin', ({(769, 33, 769, 46): '(ma_sun + ml_sun)'}, {}), '(ma_sun + ml_sun)', True, 'import numpy as np\n'), ((768, 24, 768, 38), 'numpy.sin', 'np.sin', ({(768, 31, 768, 37): 'ml_sun'}, {}), '(ml_sun)', True, 'import numpy as np\n'), ((768, 52, 768, 73), 'numpy.sin', 'np.sin', ({(768, 59, 768, 72): '(ma_sun - ml_sun)'}, {}), '(ma_sun - ml_sun)', True, 'import numpy as np\n')]