spcl.py

# -*- coding: utf-8 -*-
"""
Spectral Class Library:
    The classes Spectrum, TimeSeries and

@author: Martin Rabe, 2016-2018
"""
import numpy as _np
import matplotlib.pyplot as _plt
import copy

class Spectrum(object):
    """ Basic spectral data (X: wavelength, frequency, energy, ... and 
    Y: intensity, absorbance, ...), time stamps and methods for interpretation.
    t: should be datetime.datetime object for full functionality."""
    # ToDo :      
        
    def __init__(self, x, intensity, t = 0, xlabel = '', ylabel = '', 
                 FileName = '', comment='' ):
        
        #: x and and intensity should be instances of numpy.array         
        if isinstance(x, _np.ndarray):
            self.x = x
        else:
            self.x = _np.array(x, dtype = _np.float64)
                
        if isinstance(intensity, _np.ndarray):
            self.intensity = intensity
        else:
            self.intensity = _np.array(intensity, dtype=_np.float64)
                
        self.__linehandles=[]
        self.t = t
        self.xlabel = xlabel
        self.ylabel = ylabel
        
        #:create dataset with increasing x
        SortingInd = _np.argsort(self.x)
        self.x = self.x[SortingInd]
        self.intensity = self.intensity[SortingInd]
        self.filename = FileName
        self.comment=comment
    
    
    def __add__(self, other):
        '''Addition of a scalar to Spectrum returns a spectrum offset in 
        intensity by the scalar value. Addition of another Spectrum works only 
        for identical x and returns a Spectrum whose intensity is the sum of 
        the intensities of the input Spectrum objects.'''
        if _np.isscalar(other):
            linespec = Spectrum(self.x, 
                                _np.ones(_np.size(self.intensity))*other)
            sumSpec = self + linespec
        else:
            sumSpec = self.__addspectrum(other)
        
        return sumSpec
    
    

    def __radd__(self, other):
        return self + other


    
    def __addspectrum(self, other):
        #easiest (quickest) case when all x are equal:
        sumSpec = copy.deepcopy(self)
        if len(self.x)==len(other.x) and (self.x == other.x).all():
            newInt = self.intensity + other.intensity
        else:
            raise ValueError
        sumSpec.intensity = newInt
        return sumSpec
    
    def __mul__(self, other):
        '''multiplication only with scalar'''
        if _np.isscalar(other):
            prodSpec = copy.deepcopy(self)
            newInt = other*self.intensity
        else:
            raise ValueError
        prodSpec.intensity = newInt
        return prodSpec
    
    def __rmul__(self,other):
        return self*other
            
    def __sub__(self, other):
        return self+(-1*other)
    
    def __rsub__(self,other):
        return (-1*self) + other
    
    
    def __truediv__(self, other):
        return self*(1/other)
    #probably a plot function should not be part of a basic library (?)
    def plot(self, ax , plotColor='blue', showComment=False, 
             plotBLcorrectedPeak=False, labl=[]):
        """Plot the spectrum in the matplotlib.axes ax. Probably one should not use this"""
        
        
        if plotBLcorrectedPeak: #this is not really clean because the bl is only defined at the peak!
            self.__linehandles.append(ax.plot(self.__Peak[:,0], self.__Peak[:,1]-self.__Peak[:,2], 
                                              color=plotColor, label=labl))
        else:
            self.__linehandles.append(ax.plot(self.x, self.intensity, 
                                              color=plotColor, label=labl))
        
        if showComment:
            _plt.text(0.8,0.8,self.comment,
                     horizontalalignment='center',
                     verticalalignment='center',
                     transform = ax.transAxes)
            
        ax.yaxis.set_label_text(self.ylabel)
        ax.xaxis.set_label_text(self.xlabel)
    
    def get_t(self):
        """Return time stamp"""
        return self.t

    def get_x(self):
        """Return x as numpy array."""
        return self.x

    def get_intensity(self):
        """Return intensity as numpy array."""
        return self.intensity
        
    def get_filename(self):
        """Return the file name."""
        return self.filename

    def get_peakbaseline(self):
        """Return a linear baseline in the x-range defined by set_peaklimit()"""
        return self.__PeakBL
        
    def get_peak(self):
        """Return array containig the x-range, the intensity and the linear 
        baseline defined by set_peaklimit()"""
        return self.__Peak
    
    def get_peakarea(self):
        """Return area under peak (baseline corrected) in the x-range 
        defined by set_peaklimit()"""
        return self.__PeakArea    
    
    def set_peaklimit(self, XLim):
        """Define x-limits for peak area calculation."""
        self.peakarea_calc(min(XLim), max(XLim))    
    
    def get_intat(self, xValue, PointsToAverage=1):
        """Return mean intensity at xValue"""
        XIndex = self.get_closestxind(xValue)
        if PointsToAverage % 2 == 0:
            PointsToAverage+=1
        return _np.mean(self.intensity[int(XIndex-((PointsToAverage-1)/2)):
                                         int(1+XIndex+((PointsToAverage-1)/2))])

    def get_intat_blc(self, xValue, PointsToAverage=1):
        """Return mean intensity at xValue corrected by baseline values. 
        Attention: baseline should be defined in the region, for instance
        by set_peaklimit method. Otherwise the returned values are not useful.""" 
        xBL = self.get_peak()[:,0]
        xIndex = abs(xBL-xValue).argmin()
        if PointsToAverage % 2 == 0:
            PointsToAverage+=1
        y = self.get_intat(xValue, PointsToAverage)
        if PointsToAverage == 1:
            indexrange=PointsToAverage
        else:
            indexrange=range(int(xIndex-((PointsToAverage-1)/2)), 
                             int(xIndex+((PointsToAverage-1)/2)))
        
        return y - _np.mean(self.__PeakBL[indexrange])
        
    def get_closestxind(self, XValue):
        """Return index of value in x closest to XValue"""
        index=abs(self.x-XValue).argmin()
        return index
    
    def get_numintegral(self):
        """Return numerical integral calculated by trapezoidal rule in the 
        in the x-range defined by set_peaklimit()""" 
        return self.__AUC
    
        
    def normalize(self, normalizationValue=1):
        """tbd"""
        pass
    
    def peakarea_calc(self, xmin, xmax):
        """Return area under baseline corrected intensity between xmin and xmax"""

        self._LimIndices = abs(self.x-xmin).argmin(), abs(self.x-xmax).argmin()
        LimIndices=self._LimIndices

        # calculate area under the curve        
        self.__AUC = _np.trapz(self.intensity[LimIndices[0]:LimIndices[1]], 
                              x=self.x[LimIndices[0]:LimIndices[1]])

        # calculate area under linear baseline                              
        self.__AUBL = _np.trapz((self.intensity[LimIndices[0]], 
                                self.intensity[LimIndices[1]]), 
                                x=(self.x[LimIndices[0]], 
                                   self.x[LimIndices[1]]) )

        # calculate linear Baseline
        m=((self.intensity[LimIndices[1]]-self.intensity[LimIndices[0]])/
           (self.x[LimIndices[1]]-self.x[LimIndices[0]]))
        n=self.intensity[LimIndices[1]]-(m*self.x[LimIndices[1]])
        self.m=m
        self.n=n
        self.__PeakBL = m*self.x[LimIndices[0]:(LimIndices[1]+1)]+n
        
        #return peak area
        self.__PeakArea = self.__AUC - self.__AUBL
      
        self.__Peak = _np.array([self.x[LimIndices[0]:(LimIndices[1]+1)],
                                self.intensity[LimIndices[0]:(LimIndices[1]+1)], 
                                self.__PeakBL]).T
        
        return self.__PeakArea
    
    def get_max(self, xrange):
        """Return x and intensity of maximum intensity in xrange."""
        XInds=[self.get_closestxind(min(xrange)),
               self.get_closestxind(max(xrange))]
        X = self.x[XInds[0]:XInds[1]]
        Y = self.intensity[XInds[0]:XInds[1]]
        Xmax = X[Y.argmax()]
        Ymax = Y.max()
        maxInd = Y.argmax() #this is probably not very useful and is thus not returned anymmore
        return Xmax, Ymax 
    
    def get_min(self, xrange):
        """Return x and intensity of minimum intensity in xrange."""
        XInds=[self.get_closestxind(min(xrange)),
               self.get_closestxind(max(xrange))]
        X = self.x[XInds[0]:XInds[1]]
        Y = self.intensity[XInds[0]:XInds[1]]
        Xmin = X[Y.argmin()]
        Ymin = Y.min()
        minInd = Y.argmin() #this is probably not very useful and is thus not returned anymmore
        return Xmin, Ymin 


    def set_peakfitresults(self, peakseries, fitresults):
        """set the results of peak fitting peakseries should be SpectralSeries 
        object"""
        if isinstance(peakseries, SpectralSeries):
            self.fittedpeaks = peakseries
            self.fitresults = fitresults
        else:
            raise(TypeError)
    


class SpectralSum(object):
    """Merged Spectra whose sum is another full spectrum. Can for instance be 
    used for peak fitting results, or concatenating spectra over different 
    Energy ranges. 
    TODO: needs to be programmed, should have the same methods and attributes
    as Spectrum Class (+other specific methods & attributes) so it can be used 
    in TimeSeries Class and others."""
    def __init__(self, spectra_list):
        self.__list = spectra_list

class SeriesIterator:
   ''' Iterator class for TimeSeries and Spectralseries '''
   def __init__(self, series):
       # Team object reference
       self._series = series
       # member variable to keep track of current index
       self._index = 0
   def __next__(self):
       ''''Returns the next value from series object's lists '''
       if self._index < len(self._series.get_spectra()):
           result = self._series.get_spectra()[self._index] 
           self._index +=1
           return result
       # End of Iteration
       raise StopIteration

class TimeSeries(object):
    ''' Merged Spectrum or SpectralSum objects, giving access to interpretation 
    of time dependent spectral variations. '''

# ToDo : check input time, if 0 for all spec create arbitrary time vector        
    
    def __init__(self, ListofSpectra, FileNames = [],t0Idx = 0, TimeUnit = ''):
        
        self.__List = ListofSpectra
        self.x = ListofSpectra[0].get_x()
        self.__NumDataSets = len(self.__List)
        self.__NumDataPts = ListofSpectra[0].get_intensity().shape[0]
        self.intarray = _np.zeros((self.__NumDataPts ,self.__NumDataSets))
        self.t_abs = _np.zeros(self.__NumDataSets)        
        self.filenames = FileNames
        
        for idx, Spec in enumerate(self.__List):
            
            try:
                self.t_abs[idx] = Spec.get_t()
            except TypeError:
                self.t_abs[idx] = Spec.get_t().timestamp()
                    
            self.intarray[:,idx] = Spec.get_intensity()
            
            if any(self.x != Spec.get_x()): # x consistency check
                 raise Exception('X-range must be equal in all spectra!')

        self.calc_startdatetime(self.__List[int(t0Idx)].t)
        self.t_unit = TimeUnit
    
    def __len__(self):
        return len(self.__List)
    
    def __iter__(self):
        ''' Returns the Iterator object '''
        return SeriesIterator(self)
    
    def __getitem__(self, key):
        return self.__List[key]
    
    
    def temporalmean(self, boxwidth):
        """Return new TimeSeries object with averaged spectra. The number of 
        spectra to be averaged is determined by boxwidth. Attention: Comments 
        of the original data set are lost in this operation."""
        newspeclist = []
        for firstoldindex in _np.arange(0, len(self), boxwidth):
            specinds = _np.arange(firstoldindex, firstoldindex+boxwidth)
            if _np.isin(specinds, _np.arange(0, len(self))).all():
                newint = self.intarray[:,specinds].mean(axis=1) 
                newspeclist.append(Spectrum(self.x, newint, 
                            t = self.get_spectra()[firstoldindex].t,
                            xlabel=self.get_spectra()[firstoldindex].xlabel,
                            ylabel=self.get_spectra()[firstoldindex].ylabel,
                            comment='Automatically generated by spcl.TimeSeries.temporalmean(...)'))
        return TimeSeries(newspeclist, FileNames=self.filenames)
#        return newspeclist
    def get_spectra(self):
        """Return numpy.array of all spectrum objects"""
        return _np.array(self.__List)
        

    def get_x(self):
        """Return list of all spectrum objects"""
        return self.x
        

    def get_size(self):
        """Return number of spectrum objects"""
        return self.__NumDataSets
        

    def get_xsize(self):
        """Return number of x values"""
        return self.__NumDataPts
        

    def get_intarr(self):
        """Return array of intensities."""
        return self.intarray
        

    def get_tabs(self):
        """Return numpy array of absolute time stamps"""
        return self.t_abs
        

    def get_trel(self):
        """Return numpy array of elapsed time stamps"""
        return self.t_rel
        
    def get_filenames(self):
        """Return filenames"""
        return self.filenames
    
    def calc_startdatetime(self, t0):
        self.t_rel = self.t_abs - t0.timestamp()
        self.startdate = t0.date()
        self.starttime = t0.time()
        self.startdatetime = t0
        return
    
    def set_t0idx(self, t0Idx = 0):
        """Set the index of the initial measurement (0 by default) and return 
        numpy array of elapsed time stamps"""        
        
        self.calc_startdatetime(self.get_spectra()[int(t0Idx)].t)
        return self.t_rel        
        

    def set_t0 (self, t0):
        """Set t0 manually and return numpy array of elapsed time stamps. t0
        should be datetime.datetime object."""   
        self.calc_startdatetime(t0)
        return self.t_rel
        
    def trel_rangebool(self, t_range):
        """Return boolean mask that is true for intervall t_range in t_rel and 
        false elsewhere."""
        bols = (self.t_rel >= min(t_range)) & (self.t_rel <= max(t_range))
        return bols

    def tabs_rangebool(self, t_range):
        """Return boolean mask that is true for intervall t_range in t_abs and 
        false elsewhere. t_range maybe list of datetime.datetime objects, 
        or timestamps. TODO: this has not been properly tested,yet"""
        try:
            bols = ((self.t_abs >= min(t_range).timestamp()) & 
                    (self.t_abs <= max(t_range).timestamp()))  
        except AttributeError:
            bols = ((self.t_abs >= min(t_range)) & 
                    (self.t_abs <= max(t_range)))

        return bols
    
        
    def trel_ind(self, t):
        """Return index of t_rel value closest to t"""
        return _np.argmin(abs(self.t_rel-t))
    
    def trel_rangeind(self, t_range):
        """Return indeces for intervall t_range in t_rel"""
        return _np.where(self.trel_rangebool(t_range))[0]
    
    def tabs_rangeind(self, t_range):
        """Return indeces for intervall t_range in t_abs. t_range maybe list of
        datetime.datetime objects, or timestamps. TODO: this has not been 
        properly tested,yet"""        
        return _np.where(self.tabs_rangebool(t_range))[0]
    
    def meanint_trel(self, t_range):
        """Return mean intensities of the time interval t_range."""  
        return _np.mean(self.intarray[:,self.trel_rangebool(t_range)], axis=1)
    
    def peak_areas(self, xrange):
        """Return areas of the baseline corrected peaks between xmin and xmax"""
        PAs = _np.zeros((self.__NumDataSets))
        for idx, Spec in enumerate(self.__List):
            Spec.set_peaklimit(xrange)
            PAs[idx] = Spec.get_peakarea()
        return PAs
        
    def integrals(self, xrange):
        """Return numerical integrals in xrange""" 
        PIs = _np.zeros((self.__NumDataSets))
        for idx, Spec in enumerate(self.__List):
            Spec.set_peaklimit(xrange)
            PIs[idx] = Spec.get_numintegral()
        return PIs
    
    def intensities_at(self, x, PointsToAverage = 1):
        """Return mean intensity at x"""
        IsA = _np.zeros((self.__NumDataSets))
        for idx, Spec in enumerate(self.__List):
            IsA[idx] = Spec.get_intat(x, PointsToAverage)
        return IsA
    
    def intensities_at_blc(self, x, PointsToAverage = 1):
        """Return mean intensities at xValue corrected by baseline values. 
        Attention: baselines should be defined in the region, for instance
        by set_peaklimit method. Otherwise an error is raised or the returned 
        values are not useful.""" 
        IsABlc = _np.zeros((self.__NumDataSets))
        for idx, Spec in enumerate(self.__List):
            IsABlc[idx] = Spec.get_intat_blc(x, PointsToAverage)
        return IsABlc    
    
    def max_values(self, xrange):
        """Return maximum intensities in xrange."""
        MVs = _np.zeros((self.__NumDataSets))
        for idx, Spec in enumerate(self.__List):
            MVs[idx] = Spec.get_max(xrange)[1]
        return MVs

    def max_positions(self, xrange):
        """Return x values of maximum intensities in xrange."""
        MPs = _np.zeros((self.__NumDataSets))
        for idx, Spec in enumerate(self.__List):
            MPs[idx] = Spec.get_max(xrange)[0]
        return MPs    

    def min_values(self, xrange):
        """Return minimum intensities in xrange."""
        MVs = _np.zeros((self.__NumDataSets))
        for idx, Spec in enumerate(self.__List):
            MVs[idx] = Spec.get_min(xrange)[1]
        return MVs

    def min_positions(self, xrange):
        """Return x values of minimum intensities in xrange."""
        MPs = _np.zeros((self.__NumDataSets))
        for idx, Spec in enumerate(self.__List):
            MPs[idx] = Spec.get_min(xrange)[0]
        return MPs    

    
    def get_tunit(self):
        """Return unit of time. Not necessary when timestamps are of class 
        datetime."""
        return self.t_unit
        

    def set_tunit(self, tUnit):
        """Set unit of time. Not necessary when timestamps are of class 
        datetime."""
        self.t_unit = tUnit

class SpectralSeries(object):
    ''' Merged Spectrum or SpectralSum objects, giving access to interpretation 
    of spectral variations depending on a control variable c. Similar to 
    TimeSeries, but more general '''

    # TODO : not yet finished. For instance one could implement __add__, 
    #__mul__ etc. functions for all Series classes
    
    def __init__(self, ListofSpectra, c, FileNames = [], cName = 'c'
                 ,cUnit = ''):
        
        self.__List = ListofSpectra
        self.x = ListofSpectra[0].get_x()
        self.__NumDataSets = len(self.__List)
        self.__NumDataPts = ListofSpectra[0].get_intensity().shape[0]
        self.intarray = _np.zeros((self.__NumDataPts ,self.__NumDataSets))
        self.t_abs = _np.zeros(self.__NumDataSets)        
        self.filenames = FileNames
        self.c = c
        
        for idx, Spec in enumerate(self.__List):
            
            try:
                self.t_abs[idx] = Spec.get_t()
            except TypeError:
                self.t_abs[idx] = Spec.get_t().timestamp()
                    
            self.intarray[:,idx] = Spec.get_intensity()
            
            if any(self.x != Spec.get_x()): # x consistency check
                 raise Exception('X-range must be equal in all spectra!')

        #self.calc_startdatetime(self.__List[int(c0Idx)].t)
        self.c_name = cName
        self.c_unit = cUnit
        
        #calculate the sum of all specs:
        
        
    def __len__(self):
        return len(self.__List)
    
    def __getitem__(self, key):
        return self.__List[key]
    
    def __iter__(self):
        ''' Returns the Iterator object '''
        return SeriesIterator(self)
    
    def cmean(self, boxwidth):
        """Return new SpectralSeries object with averaged spectra. The number of 
        spectra to be averaged is determined by boxwidth. Attention: Comments 
        of the original data set are lost in this operation."""
        #TODO: needs to be programmed, below is the code for the same TimeSeries function
        
        # newspeclist = []
        # for firstoldindex in _np.arange(0, len(self), boxwidth):
        #     specinds = _np.arange(firstoldindex, firstoldindex+boxwidth)
        #     if _np.isin(specinds, _np.arange(0, len(self))).all():
        #         newint = self.intarray[:,specinds].mean(axis=1) 
        #         newspeclist.append(Spectrum(self.x, newint, 
        #                     t = self.get_spectra()[firstoldindex].t,
        #                     xlabel=self.get_spectra()[firstoldindex].xlabel,
        #                     ylabel=self.get_spectra()[firstoldindex].ylabel,
        #                     comment='Automatically generated by spcl.TimeSeries.temporalmean(...)'))
        # return TimeSeries(newspeclist, FileNames=self.filenames)

    def get_spectra(self):
        """Return numpy.array of all spectrum objects"""
        return _np.array(self.__List)
        

    def get_x(self):
        """Return list of all spectrum objects"""
        return self.x
        

    def get_size(self):
        """Return number of spectrum objects"""
        return self.__NumDataSets
        

    def get_xsize(self):
        """Return number of x values"""
        return self.__NumDataPts
        

    def get_intarr(self):
        """Return array of intensities."""
        return self.intarray
        

    def get_tabs(self):
        """Return numpy array of absolute time stamps"""
        return self.t_abs
        

    def get_trel(self):
        """Return numpy array of elapsed time stamps"""
        return self.t_rel
        
    def get_filenames(self):
        """Return filenames"""
        return self.filenames
   
    def set_peaklimit(self, xrange):
        """set peaklimits xrange= [xmin, xmax]"""
        for idx, Spec in enumerate(self.__List):
            Spec.set_peaklimit(xrange)    
    
    def get_peakbaselines(self):
        """Return numpy array with wavenumbers and numpy array with baselines"""
        BL1 = self.__List[0].get_peakbaseline()
        BLs = _np.zeros((BL1.shape[0], self.__NumDataSets))
        for idx, Spec in enumerate(self.__List):
            BLs[:, idx] = Spec.get_peakbaseline()
        return self.__List[0].get_peak()[:,0], BLs
    
    def get_peaks(self):
        """Return numpy array with wavenumbers and numpy array with peaks"""
        peak1 = self.__List[0].get_peak()
        peaks = _np.zeros((peak1.shape[0], self.__NumDataSets))
        for idx, Spec in enumerate(self.__List):
            peaks[:, idx] = Spec.get_peak()[:,1]
        return peak1[:,0], peaks
    
    def peak_areas(self, xrange):
        """Return areas of the baseline corrected peaks between xmin and xmax"""
        PAs = _np.zeros((self.__NumDataSets))
        for idx, Spec in enumerate(self.__List):
            Spec.set_peaklimit(xrange)
            PAs[idx] = Spec.get_peakarea()
        return PAs
        
    def integrals(self, xrange):
        """Return numerical integrals in xrange""" 
        PIs = _np.zeros((self.__NumDataSets))
        for idx, Spec in enumerate(self.__List):
            Spec.set_peaklimit(xrange)
            PIs[idx] = Spec.get_numintegral()
        return PIs
    
    def intensities_at(self, x, PointsToAverage = 1):
        """Return mean intensity at x"""
        IsA = _np.zeros((self.__NumDataSets))
        for idx, Spec in enumerate(self.__List):
            IsA[idx] = Spec.get_intat(x, PointsToAverage)
        return IsA
    
    def intensities_at_blc(self, x, PointsToAverage = 1):
        """Return mean intensities at xValue corrected by baseline values. 
        Attention: baselines should be defined in the region, for instance
        by set_peaklimit method. Otherwise an error is raised or the returned 
        values are not useful.""" 
        IsABlc = _np.zeros((self.__NumDataSets))
        for idx, Spec in enumerate(self.__List):
            IsABlc[idx] = Spec.get_intat_blc(x, PointsToAverage)
        return IsABlc    
    
    def max_values(self, xrange):
        """Return maximum intensities in xrange."""
        MVs = _np.zeros((self.__NumDataSets))
        for idx, Spec in enumerate(self.__List):
            MVs[idx] = Spec.get_max(xrange)[1]
        return MVs

    def max_positions(self, xrange):
        """Return x values of maximum intensities in xrange."""
        MPs = _np.zeros((self.__NumDataSets))
        for idx, Spec in enumerate(self.__List):
            MPs[idx] = Spec.get_max(xrange)[0]
        return MPs    

    def min_values(self, xrange):
        """Return minimum intensities in xrange."""
        MVs = _np.zeros((self.__NumDataSets))
        for idx, Spec in enumerate(self.__List):
            MVs[idx] = Spec.get_min(xrange)[1]
        return MVs

    def min_positions(self, xrange):
        """Return x values of minimum intensities in xrange."""
        MPs = _np.zeros((self.__NumDataSets))
        for idx, Spec in enumerate(self.__List):
            MPs[idx] = Spec.get_min(xrange)[0]
        return MPs    

    
    def get_cunit(self):
        """Return unit of time. Not necessary when timestamps are of class 
        datetime."""
        return self.c_unit
        

    def set_cunit(self, cUnit):
        """Set unit of time. Not necessary when timestamps are of class 
        datetime."""
        self.c_unit = cUnit

    def get_cname(self):
        """Return unit of time. Not necessary when timestamps are of class 
        datetime."""
        return self.c_name
        

    def set_cname(self, cName):
        """Set unit of time. Not necessary when timestamps are of class 
        datetime."""
        self.c_name = cName

# All the commented functions below probably don't make much sense for this class
    # def calc_startdatetime(self, t0):
    #     self.t_rel = self.t_abs - t0.timestamp()
    #     self.startdate = t0.date()
    #     self.starttime = t0.time()
    #     self.startdatetime = t0
    #     return
    
    # def set_c0idx(self, c0Idx = 0):
    #     """Set the index of the initial measurement (0 by default) and return 
    #     numpy array of elapsed time stamps"""        
        
    #     self.calc_startdatetime(self.get_spectra()[int(c0Idx)].t)
    #     return self.t_rel        

    # def set_t0 (self, t0):
    #     """Set t0 manually and return numpy array of elapsed time stamps. t0
    #     should be datetime.datetime object."""   
        
    #     self.calc_startdatetime(t0)
    #     return self.t_rel
        
    # def trel_rangebool(self, t_range):
    #     """Return boolean mask that is true for intervall t_range in t_rel and 
    #     false elsewhere."""
    #     bols = (self.t_rel >= min(t_range)) & (self.t_rel <= max(t_range))
    #     return bols

    # def tabs_rangebool(self, t_range):
    #     """Return boolean mask that is true for intervall t_range in t_abs and 
    #     false elsewhere. t_range maybe list of datetime.datetime objects, 
    #     or timestamps. TODO: this has not been properly tested,yet"""
    #     try:
    #         bols = ((self.t_abs >= min(t_range).timestamp()) & 
    #                 (self.t_abs <= max(t_range).timestamp()))  
    #     except AttributeError:
    #         bols = ((self.t_abs >= min(t_range)) & 
    #                 (self.t_abs <= max(t_range)))

    #     return bols
    
        
    # def trel_ind(self, t):
    #     """Return index of t_rel value closest to t"""
    #     return _np.argmin(abs(self.t_rel-t))
    
    # def trel_rangeind(self, t_range):
    #     """Return indeces for intervall t_range in t_rel"""
    #     return _np.where(self.trel_rangebool(t_range))[0]
    
    # def tabs_rangeind(self, t_range):
    #     """Return indeces for intervall t_range in t_abs. t_range maybe list of
    #     datetime.datetime objects, or timestamps. TODO: this has not been 
    #     properly tested,yet"""        
    #     return _np.where(self.tabs_rangebool(t_range))[0]
    
    # def meanint_trel(self, t_range):
    #     """Return mean intensities of the time interval t_range."""  
    #     return _np.mean(self.intarray[:,self.trel_rangebool(t_range)], axis=1)