schain Commit - r922:d680543828ae · Jicamarca Repository

Fix all PlotData, add SpectraMean, CrossSpectra plots, now Parameters extends Spectra fix bugs in ParametersProc

jespinoza -

r922:d680543828ae

parent child

Context file:

r922:d680543828ae

Collapse all files

schainpy/model/data/jrodata.py +6 -4

             '''
             $Author: murco $
             $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
             '''
             import copy
             import numpy
             import datetime
             from jroheaderIO import SystemHeader, RadarControllerHeader
             from schainpy import cSchain
             def getNumpyDtype(dataTypeCode):
                 if dataTypeCode == 0:
                     numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')])
                 elif dataTypeCode == 1:
                     numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')])
                 elif dataTypeCode == 2:
                     numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')])
                 elif dataTypeCode == 3:
                     numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')])
                 elif dataTypeCode == 4:
                     numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')])
                 elif dataTypeCode == 5:
                     numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')])
                 else:
                     raise ValueError, 'dataTypeCode was not defined'
                 return numpyDtype
             def getDataTypeCode(numpyDtype):
                 if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]):
                     datatype = 0
                 elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]):
                     datatype = 1
                 elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]):
                     datatype = 2
                 elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]):
                     datatype = 3
                 elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]):
                     datatype = 4
                 elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]):
                     datatype = 5
                 else:
                     datatype = None
                 return datatype
             def hildebrand_sekhon(data, navg):
                 """
                 This method is for the objective determination of the noise level in Doppler spectra. This
                 implementation technique is based on the fact that the standard deviation of the spectral
                 densities is equal to the mean spectral density for white Gaussian noise
                 Inputs:
                     Data    :    heights
                     navg    :    numbers of averages
                 Return:
                     -1        :    any error
                     anoise    :    noise's level
                 """
                 sortdata = numpy.sort(data,axis=None)
             #     lenOfData = len(sortdata)
             #     nums_min = lenOfData*0.2
             #
             #     if nums_min <= 5:
             #         nums_min = 5
             #
             #     sump = 0.
             #
             #     sumq = 0.
             #
             #     j = 0
             #
             #     cont = 1
             #
             #     while((cont==1)and(j<lenOfData)):
             #
             #         sump += sortdata[j]
             #
             #         sumq += sortdata[j]**2
             #
             #         if j > nums_min:
             #             rtest = float(j)/(j-1) + 1.0/navg
             #             if ((sumq*j) > (rtest*sump**2)):
             #                 j = j - 1
             #                 sump  = sump - sortdata[j]
             #                 sumq =  sumq - sortdata[j]**2
             #                 cont = 0
             #
             #         j += 1
             #
             #     lnoise = sump /j
             #
             #     return lnoise
                 return cSchain.hildebrand_sekhon(sortdata, navg)
             class Beam:
                 def __init__(self):
                     self.codeList = []
                     self.azimuthList = []
                     self.zenithList = []
             class GenericData(object):
                 flagNoData = True
                 def __init__(self):
                     raise NotImplementedError
                 def copy(self, inputObj=None):
                     if inputObj == None:
                         return copy.deepcopy(self)
                     for key in inputObj.__dict__.keys():
                         attribute = inputObj.__dict__[key]
                         #If this attribute is a tuple or list
                         if type(inputObj.__dict__[key]) in (tuple, list):
                             self.__dict__[key] = attribute[:]
                             continue
                         #If this attribute is another object or instance
                         if hasattr(attribute, '__dict__'):
                             self.__dict__[key] = attribute.copy()
                             continue
                         self.__dict__[key] = inputObj.__dict__[key]
                 def deepcopy(self):
                     return copy.deepcopy(self)
                 def isEmpty(self):
                     return self.flagNoData
             class JROData(GenericData):
             #    m_BasicHeader = BasicHeader()
             #    m_ProcessingHeader = ProcessingHeader()
                 systemHeaderObj = SystemHeader()
                 radarControllerHeaderObj = RadarControllerHeader()
             #    data = None
                 type = None
                 datatype = None #dtype but in string
             #     dtype = None
             #    nChannels = None
             #    nHeights = None
                 nProfiles = None
                 heightList = None
                 channelList = None
                 flagDiscontinuousBlock = False
                 useLocalTime = False
                 utctime = None
                 timeZone = None
                 dstFlag = None
                 errorCount = None
                 blocksize = None
             #     nCode = None
             #
             #     nBaud = None
             #
             #     code = None
                 flagDecodeData = False #asumo q la data no esta decodificada
                 flagDeflipData = False #asumo q la data no esta sin flip
                 flagShiftFFT = False
             #     ippSeconds = None
             #     timeInterval = None
                 nCohInt = None
             #     noise = None
                 windowOfFilter = 1
                 #Speed of ligth
                 C = 3e8
                 frequency = 49.92e6
                 realtime = False
                 beacon_heiIndexList = None
                 last_block = None
                 blocknow = None
                 azimuth = None
                 zenith = None
                 beam = Beam()
                 profileIndex = None
                 def __init__(self):
                     raise NotImplementedError
                 def getNoise(self):
                     raise NotImplementedError
                 def getNChannels(self):
                     return len(self.channelList)
                 def getChannelIndexList(self):
                     return range(self.nChannels)
                 def getNHeights(self):
                     return len(self.heightList)
                 def getHeiRange(self, extrapoints=0):
                     heis = self.heightList
             #        deltah = self.heightList[1] - self.heightList[0]
             #
             #        heis.append(self.heightList[-1])
                     return heis
                 def getDeltaH(self):
                     delta = self.heightList[1] - self.heightList[0]
                     return delta
                 def getltctime(self):
                     if self.useLocalTime:
                         return self.utctime - self.timeZone*60
                     return self.utctime
                 def getDatatime(self):
                     datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
                     return datatimeValue
                 def getTimeRange(self):
                     datatime = []
                     datatime.append(self.ltctime)
                     datatime.append(self.ltctime + self.timeInterval+1)
                     datatime = numpy.array(datatime)
                     return datatime
                 def getFmaxTimeResponse(self):
                     period = (10**-6)*self.getDeltaH()/(0.15)
                     PRF = 1./(period * self.nCohInt)
                     fmax = PRF
                     return fmax
                 def getFmax(self):
                     PRF = 1./(self.ippSeconds * self.nCohInt)
                     fmax = PRF
                     return fmax
                 def getVmax(self):
                     _lambda = self.C/self.frequency
                     vmax = self.getFmax() * _lambda/2
                     return vmax
                 def get_ippSeconds(self):
                     '''
                     '''
                     return self.radarControllerHeaderObj.ippSeconds
                 def set_ippSeconds(self, ippSeconds):
                     '''
                     '''
                     self.radarControllerHeaderObj.ippSeconds = ippSeconds
                     return
                 def get_dtype(self):
                     '''
                     '''
                     return getNumpyDtype(self.datatype)
                 def set_dtype(self, numpyDtype):
                     '''
                     '''
                     self.datatype = getDataTypeCode(numpyDtype)
                 def get_code(self):
                     '''
                     '''
                     return self.radarControllerHeaderObj.code
                 def set_code(self, code):
                     '''
                     '''
                     self.radarControllerHeaderObj.code = code
                     return
                 def get_ncode(self):
                     '''
                     '''
                     return self.radarControllerHeaderObj.nCode
                 def set_ncode(self, nCode):
                     '''
                     '''
                     self.radarControllerHeaderObj.nCode = nCode
                     return
                 def get_nbaud(self):
                     '''
                     '''
                     return self.radarControllerHeaderObj.nBaud
                 def set_nbaud(self, nBaud):
                     '''
                     '''
                     self.radarControllerHeaderObj.nBaud = nBaud
                     return
                 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
                 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
                 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
                 #noise = property(getNoise, "I'm the 'nHeights' property.")
                 datatime = property(getDatatime, "I'm the 'datatime' property")
                 ltctime = property(getltctime, "I'm the 'ltctime' property")
                 ippSeconds = property(get_ippSeconds, set_ippSeconds)
                 dtype = property(get_dtype, set_dtype)
             #     timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
                 code = property(get_code, set_code)
                 nCode = property(get_ncode, set_ncode)
                 nBaud = property(get_nbaud, set_nbaud)
             class Voltage(JROData):
                 #data es un numpy array de 2 dmensiones (canales, alturas)
                 data = None
                 def __init__(self):
                     '''
                     Constructor
                     '''
                     self.useLocalTime = True
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.systemHeaderObj = SystemHeader()
                     self.type = "Voltage"
                     self.data = None
             #         self.dtype = None
             #        self.nChannels = 0
             #        self.nHeights = 0
                     self.nProfiles = None
                     self.heightList = None
                     self.channelList = None
             #        self.channelIndexList = None
                     self.flagNoData = True
                     self.flagDiscontinuousBlock = False
                     self.utctime = None
                     self.timeZone = None
                     self.dstFlag = None
                     self.errorCount = None
                     self.nCohInt = None
                     self.blocksize = None
                     self.flagDecodeData = False #asumo q la data no esta decodificada
                     self.flagDeflipData = False #asumo q la data no esta sin flip
                     self.flagShiftFFT = False
                     self.flagDataAsBlock = False    #Asumo que la data es leida perfil a perfil
                     self.profileIndex = 0
                 def getNoisebyHildebrand(self, channel = None):
                     """
                     Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
                     Return:
                         noiselevel
                     """
                     if channel != None:
                         data = self.data[channel]
                         nChannels = 1
                     else:
                         data = self.data
                         nChannels = self.nChannels
                     noise = numpy.zeros(nChannels)
                     power = data * numpy.conjugate(data)
                     for thisChannel in range(nChannels):
                         if nChannels == 1:
                             daux = power[:].real
                         else:
                             daux = power[thisChannel,:].real
                         noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
                     return noise
                 def getNoise(self, type = 1, channel = None):
                     if type == 1:
                         noise = self.getNoisebyHildebrand(channel)
                     return noise
                 def getPower(self, channel = None):
                     if channel != None:
                         data = self.data[channel]
                     else:
                         data = self.data
                     power = data * numpy.conjugate(data)
                     powerdB = 10*numpy.log10(power.real)
                     powerdB = numpy.squeeze(powerdB)
                     return powerdB
                 def getTimeInterval(self):
                     timeInterval = self.ippSeconds * self.nCohInt
                     return timeInterval
                 noise = property(getNoise, "I'm the 'nHeights' property.")
                 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
             class Spectra(JROData):
                 #data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
                 data_spc = None
                 #data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
                 data_cspc = None
                 #data dc es un numpy array de 2 dmensiones (canales, alturas)
                 data_dc = None
                 #data power
                 data_pwr = None
                 nFFTPoints = None
             #     nPairs = None
                 pairsList = None
                 nIncohInt = None
                 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
                 nCohInt = None #se requiere para determinar el valor de timeInterval
                 ippFactor = None
                 profileIndex = 0
                 plotting = "spectra"
                 def __init__(self):
                     '''
                     Constructor
                     '''
                     self.useLocalTime = True
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.systemHeaderObj = SystemHeader()
                     self.type = "Spectra"
             #        self.data = None
             #         self.dtype = None
             #        self.nChannels = 0
             #        self.nHeights = 0
                     self.nProfiles = None
                     self.heightList = None
                     self.channelList = None
             #        self.channelIndexList = None
                     self.pairsList = None
                     self.flagNoData = True
                     self.flagDiscontinuousBlock = False
                     self.utctime = None
                     self.nCohInt = None
                     self.nIncohInt = None
                     self.blocksize = None
                     self.nFFTPoints = None
                     self.wavelength = None
                     self.flagDecodeData = False #asumo q la data no esta decodificada
                     self.flagDeflipData = False #asumo q la data no esta sin flip
                     self.flagShiftFFT = False
                     self.ippFactor = 1
                     #self.noise = None
                     self.beacon_heiIndexList = []
                     self.noise_estimation = None
                 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
                     """
                     Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
                     Return:
                         noiselevel
                     """
                     noise = numpy.zeros(self.nChannels)
                     for channel in range(self.nChannels):
                         daux = self.data_spc[channel,xmin_index:xmax_index,ymin_index:ymax_index]
                         noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
                     return noise
                 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
                     if self.noise_estimation is not None:
                         return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
                     else:
                         noise = self.getNoisebyHildebrand(xmin_index, xmax_index, ymin_index, ymax_index)
                         return noise
                 def getFreqRangeTimeResponse(self, extrapoints=0):
                     deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints*self.ippFactor)
                     freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
                     return freqrange
                 def getAcfRange(self, extrapoints=0):
                     deltafreq = 10./(self.getFmax() / (self.nFFTPoints*self.ippFactor))
                     freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
                     return freqrange
                 def getFreqRange(self, extrapoints=0):
                     deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
                     freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
                     return freqrange
                 def getVelRange(self, extrapoints=0):
                     deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
                     velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) #- deltav/2
                     return velrange
                 def getNPairs(self):
                     return len(self.pairsList)
                 def getPairsIndexList(self):
                     return range(self.nPairs)
                 def getNormFactor(self):
                     pwcode = 1
                     if self.flagDecodeData:
                         pwcode = numpy.sum(self.code[0]**2)
                     #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
                     normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
                     return normFactor
                 def getFlagCspc(self):
                     if self.data_cspc is None:
                         return True
                     return False
                 def getFlagDc(self):
                     if self.data_dc is None:
                         return True
                     return False
                 def getTimeInterval(self):
                     timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles
                     return timeInterval
                 def getPower(self):
                     factor = self.normFactor
                     z = self.data_spc/factor
                     z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
                     avg = numpy.average(z, axis=1)
                     return 10*numpy.log10(avg)
                 def getCoherence(self, pairsList=None, phase=False):
                     z = []
                     if pairsList is None:
                         pairsIndexList = self.pairsIndexList
                     else:
                         pairsIndexList = []
                         for pair in pairsList:
                             if pair not in self.pairsList:
                                 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
                             pairsIndexList.append(self.pairsList.index(pair))
                     for i in range(len(pairsIndexList)):
                         pair = self.pairsList[pairsIndexList[i]]
                         ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
                         powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
                         powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
                         avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
                         if phase:
                             data = numpy.arctan2(avgcoherenceComplex.imag,
                                                  avgcoherenceComplex.real)*180/numpy.pi
                         else:
                             data = numpy.abs(avgcoherenceComplex)
                         z.append(data)
                     return numpy.array(z)
                 def setValue(self, value):
                     print "This property should not be initialized"
                     return
                 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
                 pairsIndexList = property(getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
                 normFactor = property(getNormFactor, setValue, "I'm the 'getNormFactor' property.")
                 flag_cspc = property(getFlagCspc, setValue)
                 flag_dc = property(getFlagDc, setValue)
                 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
                 timeInterval = property(getTimeInterval, setValue, "I'm the 'timeInterval' property")
             class SpectraHeis(Spectra):
                 data_spc = None
                 data_cspc = None
                 data_dc = None
                 nFFTPoints = None
             #     nPairs = None
                 pairsList = None
                 nCohInt = None
                 nIncohInt = None
                 def __init__(self):
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.systemHeaderObj = SystemHeader()
                     self.type = "SpectraHeis"
             #         self.dtype = None
             #        self.nChannels = 0
             #        self.nHeights = 0
                     self.nProfiles = None
                     self.heightList = None
                     self.channelList = None
             #        self.channelIndexList = None
                     self.flagNoData = True
                     self.flagDiscontinuousBlock = False
             #         self.nPairs = 0
                     self.utctime = None
                     self.blocksize = None
                     self.profileIndex = 0
                     self.nCohInt = 1
                     self.nIncohInt = 1
                 def getNormFactor(self):
                     pwcode = 1
                     if self.flagDecodeData:
                         pwcode = numpy.sum(self.code[0]**2)
                     normFactor = self.nIncohInt*self.nCohInt*pwcode
                     return normFactor
                 def getTimeInterval(self):
                     timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
                     return timeInterval
                 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
                 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
             class Fits(JROData):
                 heightList = None
                 channelList = None
                 flagNoData = True
                 flagDiscontinuousBlock = False
                 useLocalTime = False
                 utctime = None
                 timeZone = None
             #     ippSeconds = None
             #     timeInterval = None
                 nCohInt = None
                 nIncohInt = None
                 noise = None
                 windowOfFilter = 1
                 #Speed of ligth
                 C = 3e8
                 frequency = 49.92e6
                 realtime = False
                 def __init__(self):
                     self.type = "Fits"
                     self.nProfiles = None
                     self.heightList = None
                     self.channelList = None
             #         self.channelIndexList = None
                     self.flagNoData = True
                     self.utctime = None
                     self.nCohInt = 1
                     self.nIncohInt = 1
                     self.useLocalTime = True
                     self.profileIndex = 0
             #         self.utctime = None
             #         self.timeZone = None
             #         self.ltctime = None
             #         self.timeInterval = None
             #         self.header = None
             #         self.data_header = None
             #         self.data = None
             #         self.datatime = None
             #         self.flagNoData = False
             #         self.expName = ''
             #         self.nChannels = None
             #         self.nSamples = None
             #         self.dataBlocksPerFile = None
             #         self.comments = ''
             #
                 def getltctime(self):
                     if self.useLocalTime:
                         return self.utctime - self.timeZone*60
                     return self.utctime
                 def getDatatime(self):
                     datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
                     return datatime
                 def getTimeRange(self):
                     datatime = []
                     datatime.append(self.ltctime)
                     datatime.append(self.ltctime + self.timeInterval)
                     datatime = numpy.array(datatime)
                     return datatime
                 def getHeiRange(self):
                     heis = self.heightList
                     return heis
                 def getNHeights(self):
                     return len(self.heightList)
                 def getNChannels(self):
                     return len(self.channelList)
                 def getChannelIndexList(self):
                     return range(self.nChannels)
                 def getNoise(self, type = 1):
                     #noise = numpy.zeros(self.nChannels)
                     if type == 1:
                         noise = self.getNoisebyHildebrand()
                     if type == 2:
                         noise = self.getNoisebySort()
                     if type == 3:
                         noise = self.getNoisebyWindow()
                     return noise
                 def getTimeInterval(self):
                     timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
                     return timeInterval
                 datatime = property(getDatatime, "I'm the 'datatime' property")
                 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
                 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
                 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
                 noise = property(getNoise, "I'm the 'nHeights' property.")
                 ltctime = property(getltctime, "I'm the 'ltctime' property")
                 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
             class Correlation(JROData):
                 noise = None
                 SNR = None
                 #--------------------------------------------------
                 mode = None
                 split = False
                 data_cf = None
                 lags = None
                 lagRange = None
                 pairsList = None
                 normFactor = None
                 #--------------------------------------------------
             #     calculateVelocity = None
                 nLags = None
                 nPairs = None
                 nAvg = None
                 def __init__(self):
                     '''
                     Constructor
                     '''
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.systemHeaderObj = SystemHeader()
                     self.type = "Correlation"
                     self.data = None
                     self.dtype = None
                     self.nProfiles = None
                     self.heightList = None
                     self.channelList = None
                     self.flagNoData = True
                     self.flagDiscontinuousBlock = False
                     self.utctime = None
                     self.timeZone = None
                     self.dstFlag = None
                     self.errorCount = None
                     self.blocksize = None
                     self.flagDecodeData = False #asumo q la data no esta decodificada
                     self.flagDeflipData = False #asumo q la data no esta sin flip
                     self.pairsList = None
                     self.nPoints = None
                 def getPairsList(self):
                     return self.pairsList
                 def getNoise(self, mode = 2):
                     indR = numpy.where(self.lagR == 0)[0][0]
                     indT = numpy.where(self.lagT == 0)[0][0]
                     jspectra0 = self.data_corr[:,:,indR,:]
                     jspectra = copy.copy(jspectra0)
                     num_chan = jspectra.shape[0]
                     num_hei = jspectra.shape[2]
                     freq_dc = jspectra.shape[1]/2
                     ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
                     if ind_vel[0]<0:
                         ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
                     if mode == 1:
                         jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
                     if mode == 2:
                         vel = numpy.array([-2,-1,1,2])
                         xx = numpy.zeros([4,4])
                         for fil in range(4):
                             xx[fil,:] = vel[fil]**numpy.asarray(range(4))
                         xx_inv = numpy.linalg.inv(xx)
                         xx_aux = xx_inv[0,:]
                         for ich in range(num_chan):
                             yy = jspectra[ich,ind_vel,:]
                             jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
                             junkid = jspectra[ich,freq_dc,:]<=0
                             cjunkid = sum(junkid)
                             if cjunkid.any():
                                 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
                     noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
                     return noise
                 def getTimeInterval(self):
                     timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
                     return timeInterval
                 def splitFunctions(self):
                     pairsList = self.pairsList
                     ccf_pairs = []
                     acf_pairs = []
                     ccf_ind = []
                     acf_ind = []
                     for l in range(len(pairsList)):
                         chan0 = pairsList[l][0]
                         chan1 = pairsList[l][1]
                         #Obteniendo pares de Autocorrelacion
                         if chan0 == chan1:
                             acf_pairs.append(chan0)
                             acf_ind.append(l)
                         else:
                             ccf_pairs.append(pairsList[l])
                             ccf_ind.append(l)
                     data_acf = self.data_cf[acf_ind]
                     data_ccf = self.data_cf[ccf_ind]
                     return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
                 def getNormFactor(self):
                     acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
                     acf_pairs = numpy.array(acf_pairs)
                     normFactor = numpy.zeros((self.nPairs,self.nHeights))
                     for p in range(self.nPairs):
                         pair = self.pairsList[p]
                         ch0 = pair[0]
                         ch1 = pair[1]
                         ch0_max = numpy.max(data_acf[acf_pairs==ch0,:,:], axis=1)
                         ch1_max = numpy.max(data_acf[acf_pairs==ch1,:,:], axis=1)
                         normFactor[p,:] = numpy.sqrt(ch0_max*ch1_max)
                     return normFactor
                 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
                 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
-            class Parameters(JROData):
+            class Parameters(Spectra):
                 experimentInfo = None   #Information about the experiment
                 #Information from previous data
                 inputUnit = None        #Type of data to be processed
                 operation = None        #Type of operation to parametrize
-                normFactor = None       #Normalization Factor
+                #normFactor = None       #Normalization Factor
                 groupList = None        #List of Pairs, Groups, etc
                 #Parameters
                 data_param = None       #Parameters obtained
                 data_pre = None         #Data Pre Parametrization
                 data_SNR = None         #Signal to Noise Ratio
             #     heightRange = None      #Heights
                 abscissaList = None    #Abscissa, can be velocities, lags or time
-                noise = None            #Noise Potency
+                #noise = None            #Noise Potency
                 utctimeInit = None      #Initial UTC time
                 paramInterval = None    #Time interval to calculate Parameters in seconds
                 useLocalTime = True
                 #Fitting
                 data_error = None       #Error of the estimation
                 constants = None
                 library = None
                 #Output signal
                 outputInterval = None     #Time interval to calculate output signal in seconds
                 data_output = None       #Out signal
                 nAvg = None
+                noise_estimation = None
                 def __init__(self):
                     '''
                     Constructor
                     '''
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.systemHeaderObj = SystemHeader()
                     self.type = "Parameters"
                 def getTimeRange1(self, interval):
                     datatime = []
                     if self.useLocalTime:
                         time1 = self.utctimeInit - self.timeZone*60
                     else:
                         time1 = self.utctimeInit
             #         datatime.append(self.utctimeInit)
             #         datatime.append(self.utctimeInit + self.outputInterval)
                     datatime.append(time1)
                     datatime.append(time1 + interval)
                     datatime = numpy.array(datatime)
                     return datatime

schainpy/model/graphics/jroplot_data.py +205 -27

             import os
             import zmq
             import time
             import numpy
             import datetime
             import numpy as np
             import matplotlib.pyplot as plt
             from mpl_toolkits.axes_grid1 import make_axes_locatable
             from matplotlib.ticker import FuncFormatter, LinearLocator
             from multiprocessing import Process
             from schainpy.model.proc.jroproc_base import Operation
             #plt.ion()
-            func = lambda x, pos: ('%s') %(datetime.datetime.utcfromtimestamp(x).strftime('%H:%M'))
+            func = lambda x, pos: ('%s') %(datetime.datetime.fromtimestamp(x).strftime('%H:%M'))
             d1970 = datetime.datetime(1970,1,1)
             class PlotData(Operation, Process):
                 CODE = 'Figure'
-                colormap = 'jet'
+                colormap = 'jro'
                 CONFLATE = True
                 __MAXNUMX = 80
                 __MAXNUMY = 80
                 __missing = 1E30
                 def __init__(self, **kwargs):
                     Operation.__init__(self, plot=True, **kwargs)
                     Process.__init__(self)
                     self.kwargs['code'] = self.CODE
                     self.mp = False
                     self.dataOut = None
                     self.isConfig = False
                     self.figure = None
                     self.axes = []
                     self.localtime = kwargs.pop('localtime', True)
                     self.show = kwargs.get('show', True)
                     self.save = kwargs.get('save', False)
                     self.colormap = kwargs.get('colormap', self.colormap)
+                    self.colormap_coh = kwargs.get('colormap_coh', 'jet')
+                    self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
                     self.showprofile = kwargs.get('showprofile', True)
                     self.title = kwargs.get('wintitle', '')
-                    self.xaxis = kwargs.get('xaxis', 'time')
+                    self.xaxis = kwargs.get('xaxis', 'frequency')
                     self.zmin = kwargs.get('zmin', None)
                     self.zmax = kwargs.get('zmax', None)
                     self.xmin = kwargs.get('xmin', None)
                     self.xmax = kwargs.get('xmax', None)
                     self.xrange = kwargs.get('xrange', 24)
                     self.ymin = kwargs.get('ymin', None)
                     self.ymax = kwargs.get('ymax', None)
                     self.throttle_value = 5
                 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
                     if x_buffer.shape[0] < 2:
                         return x_buffer, y_buffer, z_buffer
                     deltas = x_buffer[1:] - x_buffer[0:-1]
                     x_median = np.median(deltas)
                     index = np.where(deltas > 5*x_median)
                     if len(index[0]) != 0:
                         z_buffer[::, index[0], ::] = self.__missing
                         z_buffer = np.ma.masked_inside(z_buffer,
 .99*self.__missing,
 .01*self.__missing)
                     return x_buffer, y_buffer, z_buffer
                 def decimate(self):
                     # dx = int(len(self.x)/self.__MAXNUMX) + 1
                     dy = int(len(self.y)/self.__MAXNUMY) + 1
                     # x = self.x[::dx]
                     x = self.x
                     y = self.y[::dy]
                     z = self.z[::, ::, ::dy]
                     return x, y, z
                 def __plot(self):
                     print 'plotting...{}'.format(self.CODE)
                     if self.show:
                         self.figure.show()
                     self.plot()
-                    self.figure.suptitle('{} {} - Date:{}'.format(self.title, self.CODE.upper(),
+                    plt.tight_layout()
-                                                                  datetime.datetime.utcfromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')))
+                    self.figure.canvas.manager.set_window_title('{} {} - Date:{}'.format(self.title, self.CODE.upper(),
+                                                                                         datetime.datetime.fromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')))
                     if self.save:
                         figname = os.path.join(self.save, '{}_{}.png'.format(self.CODE,
-                                                                             datetime.datetime.utcfromtimestamp(self.times[0]).strftime('%y%m%d_%H%M%S')))
+                                                                             datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S')))
                         print 'Saving figure: {}'.format(figname)
                         self.figure.savefig(figname)
                     self.figure.canvas.draw()
                 def plot(self):
                     print 'plotting...{}'.format(self.CODE.upper())
                     return
                 def run(self):
                     print '[Starting] {}'.format(self.name)
                     context = zmq.Context()
                     receiver = context.socket(zmq.SUB)
                     receiver.setsockopt(zmq.SUBSCRIBE, '')
                     receiver.setsockopt(zmq.CONFLATE, self.CONFLATE)
                     receiver.connect("ipc:///tmp/zmq.plots")
                     while True:
                         try:
-                        #if True:
                             self.data = receiver.recv_pyobj(flags=zmq.NOBLOCK)
                             self.dataOut = self.data['dataOut']
                             self.times = self.data['times']
                             self.times.sort()
                             self.throttle_value = self.data['throttle']
                             self.min_time = self.times[0]
                             self.max_time = self.times[-1]
                             if self.isConfig is False:
                                 self.setup()
                                 self.isConfig = True
                             self.__plot()
                             if self.data['ENDED'] is True:
-                                # self.__plot()
                                 self.isConfig = False
                         except zmq.Again as e:
                             print 'Waiting for data...'
                             plt.pause(self.throttle_value)
-                            # time.sleep(3)
                 def close(self):
                     if self.dataOut:
-                        self._plot()
+                        self.__plot()
             class PlotSpectraData(PlotData):
                 CODE = 'spc'
                 colormap = 'jro'
                 CONFLATE = False
                 def setup(self):
                     ncolspan = 1
                     colspan = 1
                     self.ncols = int(numpy.sqrt(self.dataOut.nChannels)+0.9)
                     self.nrows = int(self.dataOut.nChannels*1./self.ncols + 0.9)
                     self.width = 3.6*self.ncols
                     self.height = 3.2*self.nrows
                     if self.showprofile:
                         ncolspan = 3
                         colspan = 2
                         self.width += 1.2*self.ncols
                     self.ylabel = 'Range [Km]'
                     self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList]
                     if self.figure is None:
                         self.figure = plt.figure(figsize=(self.width, self.height),
                                                  edgecolor='k',
                                                  facecolor='w')
                     else:
                         self.figure.clf()
                     n = 0
                     for y in range(self.nrows):
                         for x in range(self.ncols):
                             if n >= self.dataOut.nChannels:
                                 break
                             ax = plt.subplot2grid((self.nrows, self.ncols*ncolspan), (y, x*ncolspan), 1, colspan)
                             if self.showprofile:
                                 ax.ax_profile = plt.subplot2grid((self.nrows, self.ncols*ncolspan), (y, x*ncolspan+colspan), 1, 1)
                             ax.firsttime = True
                             self.axes.append(ax)
                             n += 1
-                    self.figure.subplots_adjust(left=0.1, right=0.95, bottom=0.15, top=0.85, wspace=0.9, hspace=0.5)
                 def plot(self):
                     if self.xaxis == "frequency":
                         x = self.dataOut.getFreqRange(1)/1000.
                         xlabel = "Frequency (kHz)"
                     elif self.xaxis == "time":
                         x = self.dataOut.getAcfRange(1)
                         xlabel = "Time (ms)"
                     else:
                         x = self.dataOut.getVelRange(1)
                         xlabel = "Velocity (m/s)"
                     y = self.dataOut.getHeiRange()
                     z = self.data[self.CODE]
                     for n, ax in enumerate(self.axes):
                         if ax.firsttime:
                             self.xmax = self.xmax if self.xmax else np.nanmax(x)
                             self.xmin = self.xmin if self.xmin else -self.xmax
                             self.ymin = self.ymin if self.ymin else np.nanmin(y)
                             self.ymax = self.ymax if self.ymax else np.nanmax(y)
                             self.zmin = self.zmin if self.zmin else np.nanmin(z)
                             self.zmax = self.zmax if self.zmax else np.nanmax(z)
                             ax.plot = ax.pcolormesh(x, y, z[n].T,
                                                     vmin=self.zmin,
                                                     vmax=self.zmax,
                                                     cmap=plt.get_cmap(self.colormap)
                                                    )
                             divider = make_axes_locatable(ax)
                             cax = divider.new_horizontal(size='3%', pad=0.05)
                             self.figure.add_axes(cax)
                             plt.colorbar(ax.plot, cax)
                             ax.set_xlim(self.xmin, self.xmax)
                             ax.set_ylim(self.ymin, self.ymax)
-                            ax.xaxis.set_major_locator(LinearLocator(5))
-                            #ax.yaxis.set_major_locator(LinearLocator(4))
                             ax.set_ylabel(self.ylabel)
                             ax.set_xlabel(xlabel)
                             ax.firsttime = False
                             if self.showprofile:
                                 ax.plot_profile= ax.ax_profile.plot(self.data['rti'][self.max_time][n], y)[0]
                                 ax.ax_profile.set_xlim(self.zmin, self.zmax)
                                 ax.ax_profile.set_ylim(self.ymin, self.ymax)
                                 ax.ax_profile.set_xlabel('dB')
                                 ax.ax_profile.grid(b=True, axis='x')
                                 ax.plot_noise = ax.ax_profile.plot(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y,
                                                                    color="k", linestyle="dashed", lw=2)[0]
                                 [tick.set_visible(False) for tick in ax.ax_profile.get_yticklabels()]
                         else:
                             ax.plot.set_array(z[n].T.ravel())
                             if self.showprofile:
                                 ax.plot_profile.set_data(self.data['rti'][self.max_time][n], y)
                                 ax.plot_noise.set_data(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y)
                         ax.set_title('{} - Noise: {:.2f} dB'.format(self.titles[n], self.data['noise'][self.max_time][n]),
                                      size=8)
+                        self.saveTime = self.max_time
+            class PlotCrossSpectraData(PlotData):
+                CODE = 'cspc'
+                zmin_coh = None
+                zmax_coh = None
+                zmin_phase = None
+                zmax_phase = None
+                CONFLATE = False
+                def setup(self):
+                    ncolspan = 1
+                    colspan = 1
+                    self.ncols = 2
+                    self.nrows = self.dataOut.nPairs
+                    self.width = 3.6*self.ncols
+                    self.height = 3.2*self.nrows
+                    self.ylabel = 'Range [Km]'
+                    self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList]
+                    if self.figure is None:
+                        self.figure = plt.figure(figsize=(self.width, self.height),
+                                                 edgecolor='k',
+                                                 facecolor='w')
+                    else:
+                        self.figure.clf()
+                    for y in range(self.nrows):
+                        for x in range(self.ncols):
+                            ax = plt.subplot2grid((self.nrows, self.ncols), (y, x), 1, 1)
+                            ax.firsttime = True
+                            self.axes.append(ax)
+                def plot(self):
+                    if self.xaxis == "frequency":
+                        x = self.dataOut.getFreqRange(1)/1000.
+                        xlabel = "Frequency (kHz)"
+                    elif self.xaxis == "time":
+                        x = self.dataOut.getAcfRange(1)
+                        xlabel = "Time (ms)"
+                    else:
+                        x = self.dataOut.getVelRange(1)
+                        xlabel = "Velocity (m/s)"
+                    y = self.dataOut.getHeiRange()
+                    z_coh = self.data['cspc_coh']
+                    z_phase = self.data['cspc_phase']
+                    for n in range(self.nrows):
+                        ax = self.axes[2*n]
+                        ax1 = self.axes[2*n+1]
+                        if ax.firsttime:
+                            self.xmax = self.xmax if self.xmax else np.nanmax(x)
+                            self.xmin = self.xmin if self.xmin else -self.xmax
+                            self.ymin = self.ymin if self.ymin else np.nanmin(y)
+                            self.ymax = self.ymax if self.ymax else np.nanmax(y)
+                            self.zmin_coh = self.zmin_coh if self.zmin_coh else 0.0
+                            self.zmax_coh = self.zmax_coh if self.zmax_coh else 1.0
+                            self.zmin_phase = self.zmin_phase if self.zmin_phase else -180
+                            self.zmax_phase = self.zmax_phase if self.zmax_phase else 180
+                            ax.plot = ax.pcolormesh(x, y, z_coh[n].T,
+                                                    vmin=self.zmin_coh,
+                                                    vmax=self.zmax_coh,
+                                                    cmap=plt.get_cmap(self.colormap_coh)
+                                                   )
+                            divider = make_axes_locatable(ax)
+                            cax = divider.new_horizontal(size='3%', pad=0.05)
+                            self.figure.add_axes(cax)
+                            plt.colorbar(ax.plot, cax)
+                            ax.set_xlim(self.xmin, self.xmax)
+                            ax.set_ylim(self.ymin, self.ymax)
+                            ax.set_ylabel(self.ylabel)
+                            ax.set_xlabel(xlabel)
+                            ax.firsttime = False
+                            ax1.plot = ax1.pcolormesh(x, y, z_phase[n].T,
+                                                    vmin=self.zmin_phase,
+                                                    vmax=self.zmax_phase,
+                                                    cmap=plt.get_cmap(self.colormap_phase)
+                                                   )
+                            divider = make_axes_locatable(ax1)
+                            cax = divider.new_horizontal(size='3%', pad=0.05)
+                            self.figure.add_axes(cax)
+                            plt.colorbar(ax1.plot, cax)
+                            ax1.set_xlim(self.xmin, self.xmax)
+                            ax1.set_ylim(self.ymin, self.ymax)
+                            ax1.set_ylabel(self.ylabel)
+                            ax1.set_xlabel(xlabel)
+                            ax1.firsttime = False
+                        else:
+                            ax.plot.set_array(z_coh[n].T.ravel())
+                            ax1.plot.set_array(z_phase[n].T.ravel())
+                        ax.set_title('Coherence Ch{} * Ch{}'.format(self.dataOut.pairsList[n][0], self.dataOut.pairsList[n][1]), size=8)
+                        ax1.set_title('Phase Ch{} * Ch{}'.format(self.dataOut.pairsList[n][0], self.dataOut.pairsList[n][1]), size=8)
+                        self.saveTime = self.max_time
+            class PlotSpectraMeanData(PlotSpectraData):
+                CODE = 'spc_mean'
+                colormap = 'jet'
+                def plot(self):
+                    if self.xaxis == "frequency":
+                        x = self.dataOut.getFreqRange(1)/1000.
+                        xlabel = "Frequency (kHz)"
+                    elif self.xaxis == "time":
+                        x = self.dataOut.getAcfRange(1)
+                        xlabel = "Time (ms)"
+                    else:
+                        x = self.dataOut.getVelRange(1)
+                        xlabel = "Velocity (m/s)"
+                    y = self.dataOut.getHeiRange()
+                    z = self.data['spc']
+                    mean = self.data['mean'][self.max_time]
+                    for n, ax in enumerate(self.axes):
+                        if ax.firsttime:
+                            self.xmax = self.xmax if self.xmax else np.nanmax(x)
+                            self.xmin = self.xmin if self.xmin else -self.xmax
+                            self.ymin = self.ymin if self.ymin else np.nanmin(y)
+                            self.ymax = self.ymax if self.ymax else np.nanmax(y)
+                            self.zmin = self.zmin if self.zmin else np.nanmin(z)
+                            self.zmax = self.zmax if self.zmax else np.nanmax(z)
+                            ax.plt = ax.pcolormesh(x, y, z[n].T,
+                                                   vmin=self.zmin,
+                                                   vmax=self.zmax,
+                                                   cmap=plt.get_cmap(self.colormap)
+                                                   )
+                            ax.plt_dop = ax.plot(mean[n], y,
+                                                 color='k')[0]
+                            divider = make_axes_locatable(ax)
+                            cax = divider.new_horizontal(size='3%', pad=0.05)
+                            self.figure.add_axes(cax)
+                            plt.colorbar(ax.plt, cax)
+                            ax.set_xlim(self.xmin, self.xmax)
+                            ax.set_ylim(self.ymin, self.ymax)
+                            ax.set_ylabel(self.ylabel)
+                            ax.set_xlabel(xlabel)
+                            ax.firsttime = False
+                            if self.showprofile:
+                                ax.plt_profile= ax.ax_profile.plot(self.data['rti'][self.max_time][n], y)[0]
+                                ax.ax_profile.set_xlim(self.zmin, self.zmax)
+                                ax.ax_profile.set_ylim(self.ymin, self.ymax)
+                                ax.ax_profile.set_xlabel('dB')
+                                ax.ax_profile.grid(b=True, axis='x')
+                                ax.plt_noise = ax.ax_profile.plot(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y,
+                                                                   color="k", linestyle="dashed", lw=2)[0]
+                                [tick.set_visible(False) for tick in ax.ax_profile.get_yticklabels()]
+                        else:
+                            ax.plt.set_array(z[n].T.ravel())
+                            ax.plt_dop.set_data(mean[n], y)
+                            if self.showprofile:
+                                ax.plt_profile.set_data(self.data['rti'][self.max_time][n], y)
+                                ax.plt_noise.set_data(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y)
+                        ax.set_title('{} - Noise: {:.2f} dB'.format(self.titles[n], self.data['noise'][self.max_time][n]),
+                                     size=8)
+                        self.saveTime = self.max_time
             class PlotRTIData(PlotData):
                 CODE = 'rti'
                 colormap = 'jro'
                 def setup(self):
                     self.ncols = 1
                     self.nrows = self.dataOut.nChannels
                     self.width = 10
-                    self.height = 2.2*self.nrows
+                    self.height = 2.2*self.nrows if self.nrows<6 else 12
                     if self.nrows==1:
                         self.height += 1
                     self.ylabel = 'Range [Km]'
                     self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList]
                     if self.figure is None:
                         self.figure = plt.figure(figsize=(self.width, self.height),
                                                  edgecolor='k',
                                                  facecolor='w')
                     else:
                         self.figure.clf()
                         self.axes = []
                     for n in range(self.nrows):
                         ax = self.figure.add_subplot(self.nrows, self.ncols, n+1)
                         ax.firsttime = True
                         self.axes.append(ax)
-                    self.figure.subplots_adjust(hspace=0.5)
                 def plot(self):
                     self.x = np.array(self.times)
                     self.y = self.dataOut.getHeiRange()
                     self.z = []
                     for ch in range(self.nrows):
                         self.z.append([self.data[self.CODE][t][ch] for t in self.times])
                     self.z = np.array(self.z)
                     for n, ax in enumerate(self.axes):
                         x, y, z = self.fill_gaps(*self.decimate())
                         xmin = self.min_time
                         xmax = xmin+self.xrange*60*60
                         if ax.firsttime:
                             self.ymin = self.ymin if self.ymin else np.nanmin(self.y)
                             self.ymax = self.ymax if self.ymax else np.nanmax(self.y)
                             self.zmin = self.zmin if self.zmin else np.nanmin(self.z)
                             self.zmax = self.zmax if self.zmax else np.nanmax(self.z)
                             plot = ax.pcolormesh(x, y, z[n].T,
                                                  vmin=self.zmin,
                                                  vmax=self.zmax,
                                                  cmap=plt.get_cmap(self.colormap)
                                                 )
                             divider = make_axes_locatable(ax)
                             cax = divider.new_horizontal(size='2%', pad=0.05)
                             self.figure.add_axes(cax)
                             plt.colorbar(plot, cax)
                             ax.set_ylim(self.ymin, self.ymax)
-                            if self.xaxis == 'time':
-                                ax.xaxis.set_major_formatter(FuncFormatter(func))
-                                ax.xaxis.set_major_locator(LinearLocator(6))
-                            # ax.yaxis.set_major_locator(LinearLocator(4))
+                            ax.xaxis.set_major_formatter(FuncFormatter(func))
+                            ax.xaxis.set_major_locator(LinearLocator(6))
                             ax.set_ylabel(self.ylabel)
                             # if self.xmin is None:
                             #     xmin = self.min_time
                             # else:
                             #     xmin = (datetime.datetime.combine(self.dataOut.datatime.date(),
                             #                                      datetime.time(self.xmin, 0, 0))-d1970).total_seconds()
                             ax.set_xlim(xmin, xmax)
                             ax.firsttime = False
                         else:
                             ax.collections.remove(ax.collections[0])
                             ax.set_xlim(xmin, xmax)
                             plot = ax.pcolormesh(x, y, z[n].T,
                                                  vmin=self.zmin,
                                                  vmax=self.zmax,
                                                  cmap=plt.get_cmap(self.colormap)
                                                 )
-                            ax.set_title('{} {}'.format(self.titles[n],
+                        ax.set_title('{} {}'.format(self.titles[n],
-                                                        datetime.datetime.utcfromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')),
+                                                    datetime.datetime.fromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')),
-                                         size=8)
+                                     size=8)
+                        self.saveTime = self.min_time
             class PlotCOHData(PlotRTIData):
                 CODE = 'coh'
                 def setup(self):
                     self.ncols = 1
                     self.nrows = self.dataOut.nPairs
                     self.width = 10
-                    self.height = 2.2*self.nrows
+                    self.height = 2.2*self.nrows if self.nrows<6 else 12
                     if self.nrows==1:
                         self.height += 1
                     self.ylabel = 'Range [Km]'
-                    self.titles = ['Channels {}'.format(x) for x in self.dataOut.pairsList]
+                    self.titles = ['{} Ch{} * Ch{}'.format(self.CODE.upper(), x[0], x[1]) for x in self.dataOut.pairsList]
                     if self.figure is None:
                         self.figure = plt.figure(figsize=(self.width, self.height),
                                                  edgecolor='k',
                                                  facecolor='w')
                     else:
                         self.figure.clf()
                         self.axes = []
                     for n in range(self.nrows):
                         ax = self.figure.add_subplot(self.nrows, self.ncols, n+1)
                         ax.firsttime = True
                         self.axes.append(ax)
-                    self.figure.subplots_adjust(hspace=0.5)
             class PlotNoiseData(PlotData):
                 CODE = 'noise'
                 def setup(self):
                     self.ncols = 1
                     self.nrows = 1
                     self.width = 10
                     self.height = 3.2
                     self.ylabel = 'Intensity [dB]'
                     self.titles = ['Noise']
                     if self.figure is None:
                         self.figure = plt.figure(figsize=(self.width, self.height),
                                                  edgecolor='k',
                                                  facecolor='w')
                     else:
                         self.figure.clf()
                         self.axes = []
                     self.ax = self.figure.add_subplot(self.nrows, self.ncols, 1)
                     self.ax.firsttime = True
                 def plot(self):
                     x = self.times
                     xmin = self.min_time
                     xmax = xmin+self.xrange*60*60
                     if self.ax.firsttime:
                         for ch in self.dataOut.channelList:
                             y = [self.data[self.CODE][t][ch] for t in self.times]
                             self.ax.plot(x, y, lw=1, label='Ch{}'.format(ch))
                         self.ax.firsttime = False
                         self.ax.xaxis.set_major_formatter(FuncFormatter(func))
                         self.ax.xaxis.set_major_locator(LinearLocator(6))
                         self.ax.set_ylabel(self.ylabel)
                         plt.legend()
                     else:
                         for ch in self.dataOut.channelList:
                             y = [self.data[self.CODE][t][ch] for t in self.times]
                             self.ax.lines[ch].set_data(x, y)
                     self.ax.set_xlim(xmin, xmax)
                     self.ax.set_ylim(min(y)-5, max(y)+5)
+                    self.saveTime = self.min_time
             class PlotSNRData(PlotRTIData):
                 CODE = 'snr'
+                colormap = 'jet'
             class PlotDOPData(PlotRTIData):
                 CODE = 'dop'
                 colormap = 'jet'
             class PlotPHASEData(PlotCOHData):
                 CODE = 'phase'
                 colormap = 'seismic'

schainpy/model/proc/jroproc_parameters.py +18 -10

             import numpy
             import math
             from scipy import optimize, interpolate, signal, stats, ndimage
             import re
             import datetime
             import copy
             import sys
             import importlib
             import itertools
             from jroproc_base import ProcessingUnit, Operation
             from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
             class ParametersProc(ProcessingUnit):
                 nSeconds = None
                 def __init__(self):
                     ProcessingUnit.__init__(self)
             #         self.objectDict = {}
                     self.buffer = None
                     self.firstdatatime = None
                     self.profIndex = 0
                     self.dataOut = Parameters()
                 def __updateObjFromInput(self):
                     self.dataOut.inputUnit = self.dataIn.type
                     self.dataOut.timeZone = self.dataIn.timeZone
                     self.dataOut.dstFlag = self.dataIn.dstFlag
                     self.dataOut.errorCount = self.dataIn.errorCount
                     self.dataOut.useLocalTime = self.dataIn.useLocalTime
                     self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
                     self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
                     self.dataOut.channelList = self.dataIn.channelList
                     self.dataOut.heightList = self.dataIn.heightList
                     self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
             #         self.dataOut.nHeights = self.dataIn.nHeights
             #         self.dataOut.nChannels = self.dataIn.nChannels
                     self.dataOut.nBaud = self.dataIn.nBaud
                     self.dataOut.nCode = self.dataIn.nCode
                     self.dataOut.code = self.dataIn.code
             #        self.dataOut.nProfiles = self.dataOut.nFFTPoints
                     self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
             #         self.dataOut.utctime = self.firstdatatime
                     self.dataOut.utctime = self.dataIn.utctime
                     self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
                     self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
                     self.dataOut.nCohInt = self.dataIn.nCohInt
             #        self.dataOut.nIncohInt = 1
                     self.dataOut.ippSeconds = self.dataIn.ippSeconds
             #        self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
-                    self.dataOut.timeInterval = self.dataIn.timeInterval
+            #        self.dataOut.timeInterval = self.dataIn.timeInterval
                     self.dataOut.heightList = self.dataIn.getHeiRange()
                     self.dataOut.frequency = self.dataIn.frequency
-                    self.dataOut.noise = self.dataIn.noise
+                    #self.dataOut.noise = self.dataIn.noise
                 def run(self):
                     #----------------------    Voltage Data    ---------------------------
                     if self.dataIn.type == "Voltage":
                         self.__updateObjFromInput()
                         self.dataOut.data_pre = self.dataIn.data.copy()
                         self.dataOut.flagNoData = False
                         self.dataOut.utctimeInit = self.dataIn.utctime
                         self.dataOut.paramInterval = self.dataIn.nProfiles*self.dataIn.nCohInt*self.dataIn.ippSeconds
                         return
                     #----------------------    Spectra Data    ---------------------------
                     if self.dataIn.type == "Spectra":
-                        self.dataOut.data_pre = (self.dataIn.data_spc,self.dataIn.data_cspc)
+                        self.dataOut.data_pre = (self.dataIn.data_spc, self.dataIn.data_cspc)
-                        self.dataOut.abscissaList = self.dataIn.getVelRange(1)
+                        self.dataOut.data_spc = self.dataIn.data_spc
-                        # self.dataOut.noise = self.dataIn.getNoise()
+                        self.dataOut.data_cspc = self.dataIn.data_cspc
-                        self.dataOut.normFactor = self.dataIn.normFactor
+                        self.dataOut.nProfiles = self.dataIn.nProfiles
+                        self.dataOut.nIncohInt = self.dataIn.nIncohInt
+                        self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
+                        self.dataOut.ippFactor = self.dataIn.ippFactor
+                        #self.dataOut.normFactor = self.dataIn.getNormFactor()
+                        self.dataOut.pairsList = self.dataIn.pairsList
                         self.dataOut.groupList = self.dataIn.pairsList
+                        self.dataOut.abscissaList = self.dataIn.getVelRange(1)
                         self.dataOut.flagNoData = False
                     #----------------------    Correlation Data    ---------------------------
                     if self.dataIn.type == "Correlation":
                         acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.dataIn.splitFunctions()
                         self.dataOut.data_pre = (self.dataIn.data_cf[acf_ind,:], self.dataIn.data_cf[ccf_ind,:,:])
                         self.dataOut.normFactor = (self.dataIn.normFactor[acf_ind,:], self.dataIn.normFactor[ccf_ind,:])
                         self.dataOut.groupList = (acf_pairs, ccf_pairs)
                         self.dataOut.abscissaList = self.dataIn.lagRange
                         self.dataOut.noise = self.dataIn.noise
                         self.dataOut.data_SNR = self.dataIn.SNR
                         self.dataOut.flagNoData = False
                         self.dataOut.nAvg = self.dataIn.nAvg
                     #----------------------    Parameters Data    ---------------------------
                     if self.dataIn.type == "Parameters":
                         self.dataOut.copy(self.dataIn)
                         self.dataOut.utctimeInit = self.dataIn.utctime
                         self.dataOut.flagNoData = False
                         return True
                     self.__updateObjFromInput()
                     self.dataOut.utctimeInit = self.dataIn.utctime
                     self.dataOut.paramInterval = self.dataIn.timeInterval
                     return
             class SpectralMoments(Operation):
                 '''
                     Function SpectralMoments()
                     Calculates moments (power, mean, standard deviation) and SNR of the signal
                     Type of dataIn:    Spectra
                     Configuration Parameters:
                         dirCosx    :     Cosine director in X axis
                         dirCosy    :     Cosine director in Y axis
                         elevation  :
                         azimuth    :
                     Input:
                         channelList    :    simple channel list to select e.g. [2,3,7]
                         self.dataOut.data_pre        :    Spectral data
                         self.dataOut.abscissaList    :    List of frequencies
                         self.dataOut.noise           :    Noise level per channel
                     Affected:
                         self.dataOut.data_param      :    Parameters per channel
                         self.dataOut.data_SNR        :    SNR per channel
                 '''
                 def run(self, dataOut):
-                    dataOut.data_pre = dataOut.data_pre[0]
+                    #dataOut.data_pre = dataOut.data_pre[0]
-                    data = dataOut.data_pre
+                    data = dataOut.data_pre[0]
                     absc = dataOut.abscissaList[:-1]
                     noise = dataOut.noise
                     nChannel = data.shape[0]
                     data_param = numpy.zeros((nChannel, 4, data.shape[2]))
                     for ind in range(nChannel):
                         data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])
                     dataOut.data_param = data_param[:,1:,:]
                     dataOut.data_SNR = data_param[:,0]
                     dataOut.data_DOP = data_param[:,1]
+                    dataOut.data_MEAN = data_param[:,2]
+                    dataOut.data_STD = data_param[:,3]
                     return
                 def __calculateMoments(self, oldspec, oldfreq, n0, nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
                     if (nicoh is None): nicoh = 1
                     if (graph is None): graph = 0
                     if (smooth is None): smooth = 0
                     elif (self.smooth < 3): smooth = 0
                     if (type1 is None): type1 = 0
                     if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
                     if (snrth is None): snrth = -3
                     if (dc is None): dc = 0
                     if (aliasing is None): aliasing = 0
                     if (oldfd is None): oldfd = 0
                     if (wwauto is None): wwauto = 0
                     if (n0 < 1.e-20):   n0 = 1.e-20
                     freq = oldfreq
                     vec_power = numpy.zeros(oldspec.shape[1])
                     vec_fd = numpy.zeros(oldspec.shape[1])
                     vec_w = numpy.zeros(oldspec.shape[1])
                     vec_snr = numpy.zeros(oldspec.shape[1])
                     for ind in range(oldspec.shape[1]):
                         spec = oldspec[:,ind]
                         aux = spec*fwindow
                         max_spec = aux.max()
                         m = list(aux).index(max_spec)
                         #Smooth
                         if (smooth == 0):   spec2 = spec
                         else:   spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
                         #    Calculo de Momentos
                         bb = spec2[range(m,spec2.size)]
                         bb = (bb<n0).nonzero()
                         bb = bb[0]
                         ss = spec2[range(0,m + 1)]
                         ss = (ss<n0).nonzero()
                         ss = ss[0]
                         if (bb.size == 0):
                             bb0 = spec.size - 1 - m
                         else:
                             bb0 = bb[0] - 1
                             if (bb0 < 0):
                                 bb0 = 0
                         if (ss.size == 0):   ss1 = 1
                         else: ss1 = max(ss) + 1
                         if (ss1 > m):   ss1 = m
                         valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
                         power = ((spec2[valid] - n0)*fwindow[valid]).sum()
                         fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
                         w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
                         snr = (spec2.mean()-n0)/n0
                         if (snr < 1.e-20) :
                             snr = 1.e-20
                         vec_power[ind] = power
                         vec_fd[ind] = fd
                         vec_w[ind] = w
                         vec_snr[ind] = snr
                     moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
                     return moments
                 #------------------    Get SA Parameters    --------------------------
                 def GetSAParameters(self):
                     #SA en frecuencia
                     pairslist = self.dataOut.groupList
                     num_pairs = len(pairslist)
                     vel = self.dataOut.abscissaList
-                    spectra = self.dataOut.data_pre
+                    spectra = self.dataOut.data_pre[0]
-                    cspectra = self.dataIn.data_cspc
+                    cspectra = self.dataOut.data_pre[1]
                     delta_v = vel[1] - vel[0]
                     #Calculating the power spectrum
                     spc_pow = numpy.sum(spectra, 3)*delta_v
                     #Normalizing Spectra
                     norm_spectra = spectra/spc_pow
                     #Calculating the norm_spectra at peak
                     max_spectra = numpy.max(norm_spectra, 3)
                     #Normalizing Cross Spectra
                     norm_cspectra = numpy.zeros(cspectra.shape)
                     for i in range(num_chan):
                         norm_cspectra[i,:,:] = cspectra[i,:,:]/numpy.sqrt(spc_pow[pairslist[i][0],:]*spc_pow[pairslist[i][1],:])
                     max_cspectra = numpy.max(norm_cspectra,2)
                     max_cspectra_index = numpy.argmax(norm_cspectra, 2)
                     for i in range(num_pairs):
                         cspc_par[i,:,:] = __calculateMoments(norm_cspectra)
                 #-------------------    Get Lags    ----------------------------------
             class SALags(Operation):
                 '''
                 Function GetMoments()
                 Input:
                     self.dataOut.data_pre
                     self.dataOut.abscissaList
                     self.dataOut.noise
                     self.dataOut.normFactor
                     self.dataOut.data_SNR
                     self.dataOut.groupList
                     self.dataOut.nChannels
                 Affected:
                     self.dataOut.data_param
                 '''
                 def run(self, dataOut):
                     data_acf = dataOut.data_pre[0]
                     data_ccf = dataOut.data_pre[1]
                     normFactor_acf = dataOut.normFactor[0]
                     normFactor_ccf = dataOut.normFactor[1]
                     pairs_acf = dataOut.groupList[0]
                     pairs_ccf = dataOut.groupList[1]
                     nHeights = dataOut.nHeights
                     absc = dataOut.abscissaList
                     noise = dataOut.noise
                     SNR = dataOut.data_SNR
                     nChannels = dataOut.nChannels
             #         pairsList = dataOut.groupList
             #         pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
                     for l in range(len(pairs_acf)):
                         data_acf[l,:,:] = data_acf[l,:,:]/normFactor_acf[l,:]
                     for l in range(len(pairs_ccf)):
                         data_ccf[l,:,:] = data_ccf[l,:,:]/normFactor_ccf[l,:]
                     dataOut.data_param = numpy.zeros((len(pairs_ccf)*2 + 1, nHeights))
                     dataOut.data_param[:-1,:] = self.__calculateTaus(data_acf, data_ccf, absc)
                     dataOut.data_param[-1,:] = self.__calculateLag1Phase(data_acf, absc)
                     return
             #     def __getPairsAutoCorr(self, pairsList, nChannels):
             #
             #         pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
             #
             #         for l in range(len(pairsList)):
             #             firstChannel = pairsList[l][0]
             #             secondChannel = pairsList[l][1]
             #
             #             #Obteniendo pares de Autocorrelacion
             #             if firstChannel == secondChannel:
             #                 pairsAutoCorr[firstChannel] = int(l)
             #
             #         pairsAutoCorr = pairsAutoCorr.astype(int)
             #
             #         pairsCrossCorr = range(len(pairsList))
             #         pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
             #
             #         return pairsAutoCorr, pairsCrossCorr
                 def __calculateTaus(self, data_acf, data_ccf, lagRange):
                     lag0 = data_acf.shape[1]/2
                     #Funcion de Autocorrelacion
                     mean_acf = stats.nanmean(data_acf, axis = 0)
                     #Obtencion Indice de TauCross
                     ind_ccf = data_ccf.argmax(axis = 1)
                     #Obtencion Indice de TauAuto
                     ind_acf = numpy.zeros(ind_ccf.shape,dtype = 'int')
                     ccf_lag0 = data_ccf[:,lag0,:]
                     for i in range(ccf_lag0.shape[0]):
                         ind_acf[i,:] = numpy.abs(mean_acf - ccf_lag0[i,:]).argmin(axis = 0)
                     #Obtencion de TauCross y TauAuto
                     tau_ccf = lagRange[ind_ccf]
                     tau_acf  = lagRange[ind_acf]
                     Nan1, Nan2 = numpy.where(tau_ccf == lagRange[0])
                     tau_ccf[Nan1,Nan2] = numpy.nan
                     tau_acf[Nan1,Nan2] = numpy.nan
                     tau = numpy.vstack((tau_ccf,tau_acf))
                     return tau
                 def __calculateLag1Phase(self, data, lagTRange):
                     data1 = stats.nanmean(data, axis = 0)
                     lag1 = numpy.where(lagTRange == 0)[0][0] + 1
                     phase = numpy.angle(data1[lag1,:])
                     return phase
             class SpectralFitting(Operation):
                 '''
                     Function GetMoments()
                     Input:
                     Output:
                     Variables modified:
                 '''
                 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None):
                     if path != None:
                         sys.path.append(path)
                     self.dataOut.library = importlib.import_module(file)
                     #To be inserted as a parameter
                     groupArray = numpy.array(groupList)
             #         groupArray = numpy.array([[0,1],[2,3]])
                     self.dataOut.groupList = groupArray
                     nGroups = groupArray.shape[0]
                     nChannels = self.dataIn.nChannels
                     nHeights=self.dataIn.heightList.size
                     #Parameters Array
                     self.dataOut.data_param = None
                     #Set constants
                     constants = self.dataOut.library.setConstants(self.dataIn)
                     self.dataOut.constants = constants
                     M = self.dataIn.normFactor
                     N = self.dataIn.nFFTPoints
                     ippSeconds = self.dataIn.ippSeconds
                     K = self.dataIn.nIncohInt
                     pairsArray = numpy.array(self.dataIn.pairsList)
                     #List of possible combinations
                     listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
                     indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
                     if getSNR:
                         listChannels = groupArray.reshape((groupArray.size))
                         listChannels.sort()
                         noise = self.dataIn.getNoise()
                         self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
                     for i in range(nGroups):
                         coord = groupArray[i,:]
                         #Input data array
                         data = self.dataIn.data_spc[coord,:,:]/(M*N)
                         data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
                         #Cross Spectra data array for Covariance Matrixes
                         ind = 0
                         for pairs in listComb:
                             pairsSel = numpy.array([coord[x],coord[y]])
                             indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
                             ind += 1
                         dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
                         dataCross = dataCross**2/K
                         for h in range(nHeights):
             #                 print self.dataOut.heightList[h]
                             #Input
                             d = data[:,h]
                             #Covariance Matrix
                             D = numpy.diag(d**2/K)
                             ind = 0
                             for pairs in listComb:
                                 #Coordinates in Covariance Matrix
                                 x = pairs[0]
                                 y = pairs[1]
                                 #Channel Index
                                 S12 = dataCross[ind,:,h]
                                 D12 = numpy.diag(S12)
                                 #Completing Covariance Matrix with Cross Spectras
                                 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
                                 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
                                 ind += 1
                             Dinv=numpy.linalg.inv(D)
                             L=numpy.linalg.cholesky(Dinv)
                             LT=L.T
                             dp = numpy.dot(LT,d)
                             #Initial values
                             data_spc = self.dataIn.data_spc[coord,:,h]
                             if (h>0)and(error1[3]<5):
                                 p0 = self.dataOut.data_param[i,:,h-1]
                             else:
                                 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
                             try:
                                 #Least Squares
                                 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
             #                   minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
                                 #Chi square error
                                 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
                                 #Error with Jacobian
                                 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
                             except:
                                 minp = p0*numpy.nan
                                 error0 = numpy.nan
                                 error1 = p0*numpy.nan
                             #Save
                             if self.dataOut.data_param is None:
                                 self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
                                 self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
                             self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
                             self.dataOut.data_param[i,:,h] = minp
                     return
                 def __residFunction(self, p, dp, LT, constants):
                     fm = self.dataOut.library.modelFunction(p, constants)
                     fmp=numpy.dot(LT,fm)
                     return  dp-fmp
                 def __getSNR(self, z, noise):
                     avg = numpy.average(z, axis=1)
                     SNR = (avg.T-noise)/noise
                     SNR = SNR.T
                     return SNR
                 def __chisq(p,chindex,hindex):
                     #similar to Resid but calculates CHI**2
                     [LT,d,fm]=setupLTdfm(p,chindex,hindex)
                     dp=numpy.dot(LT,d)
                     fmp=numpy.dot(LT,fm)
                     chisq=numpy.dot((dp-fmp).T,(dp-fmp))
                     return chisq
             class WindProfiler(Operation):
                 __isConfig = False
                 __initime = None
                 __lastdatatime = None
                 __integrationtime = None
                 __buffer = None
                 __dataReady = False
                 __firstdata = None
                 n = None
                 def __calculateCosDir(self, elev, azim):
                     zen = (90 - elev)*numpy.pi/180
                     azim = azim*numpy.pi/180
                     cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
                     cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
                     signX = numpy.sign(numpy.cos(azim))
                     signY = numpy.sign(numpy.sin(azim))
                     cosDirX = numpy.copysign(cosDirX, signX)
                     cosDirY = numpy.copysign(cosDirY, signY)
                     return cosDirX, cosDirY
                 def __calculateAngles(self, theta_x, theta_y, azimuth):
                     dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
                     zenith_arr = numpy.arccos(dir_cosw)
                     azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
                     dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
                     dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
                     return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
                 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
             #
                     if horOnly:
                         A = numpy.c_[dir_cosu,dir_cosv]
                     else:
                         A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
                     A = numpy.asmatrix(A)
                     A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
                     return A1
                 def __correctValues(self, heiRang, phi, velRadial, SNR):
                     listPhi = phi.tolist()
                     maxid = listPhi.index(max(listPhi))
                     minid = listPhi.index(min(listPhi))
                     rango = range(len(phi))
                #     rango = numpy.delete(rango,maxid)
                     heiRang1 = heiRang*math.cos(phi[maxid])
                     heiRangAux = heiRang*math.cos(phi[minid])
                     indOut = (heiRang1 < heiRangAux[0]).nonzero()
                     heiRang1 = numpy.delete(heiRang1,indOut)
                     velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
                     SNR1 = numpy.zeros([len(phi),len(heiRang1)])
                     for i in rango:
                         x = heiRang*math.cos(phi[i])
                         y1 = velRadial[i,:]
                         f1 = interpolate.interp1d(x,y1,kind = 'cubic')
                         x1 = heiRang1
                         y11 = f1(x1)
                         y2 = SNR[i,:]
                         f2 = interpolate.interp1d(x,y2,kind = 'cubic')
                         y21 = f2(x1)
                         velRadial1[i,:] = y11
                         SNR1[i,:] = y21
                     return heiRang1, velRadial1, SNR1
                 def __calculateVelUVW(self, A, velRadial):
                     #Operacion Matricial
             #         velUVW = numpy.zeros((velRadial.shape[1],3))
             #         for ind in range(velRadial.shape[1]):
             #             velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
             #         velUVW = velUVW.transpose()
                     velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
                     velUVW[:,:] = numpy.dot(A,velRadial)
                     return velUVW
             #     def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
                 def techniqueDBS(self, kwargs):
                     """
                     Function that implements Doppler Beam Swinging (DBS) technique.
                     Input:    Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
                                 Direction correction (if necessary), Ranges and SNR
                     Output:    Winds estimation (Zonal, Meridional and Vertical)
                     Parameters affected:    Winds, height range, SNR
                     """
                     velRadial0 = kwargs['velRadial']
                     heiRang = kwargs['heightList']
                     SNR0 = kwargs['SNR']
                     if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):
                         theta_x = numpy.array(kwargs['dirCosx'])
                         theta_y = numpy.array(kwargs['dirCosy'])
                     else:
                         elev = numpy.array(kwargs['elevation'])
                         azim = numpy.array(kwargs['azimuth'])
                         theta_x, theta_y = self.__calculateCosDir(elev, azim)
                     azimuth = kwargs['correctAzimuth']
                     if kwargs.has_key('horizontalOnly'):
                         horizontalOnly = kwargs['horizontalOnly']
                     else:   horizontalOnly = False
                     if kwargs.has_key('correctFactor'):
                         correctFactor = kwargs['correctFactor']
                     else:   correctFactor = 1
                     if kwargs.has_key('channelList'):
                         channelList = kwargs['channelList']
                         if len(channelList) == 2:
                             horizontalOnly = True
                         arrayChannel = numpy.array(channelList)
                         param = param[arrayChannel,:,:]
                         theta_x = theta_x[arrayChannel]
                         theta_y = theta_y[arrayChannel]
                     azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
                     heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0)
                     A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
                     #Calculo de Componentes de la velocidad con DBS
                     winds = self.__calculateVelUVW(A,velRadial1)
                     return winds, heiRang1, SNR1
                 def __calculateDistance(self, posx, posy, pairs_ccf, azimuth = None):
                     nPairs = len(pairs_ccf)
                     posx = numpy.asarray(posx)
                     posy = numpy.asarray(posy)
                     #Rotacion Inversa para alinear con el azimuth
                     if azimuth!= None:
                         azimuth = azimuth*math.pi/180
                         posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
                         posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
                     else:
                         posx1 = posx
                         posy1 = posy
                     #Calculo de Distancias
                     distx = numpy.zeros(nPairs)
                     disty = numpy.zeros(nPairs)
                     dist = numpy.zeros(nPairs)
                     ang = numpy.zeros(nPairs)
                     for i in range(nPairs):
                         distx[i] = posx1[pairs_ccf[i][1]] - posx1[pairs_ccf[i][0]]
                         disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]]
                         dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
                         ang[i] = numpy.arctan2(disty[i],distx[i])
                     return distx, disty, dist, ang
                     #Calculo de Matrices
             #         nPairs = len(pairs)
             #         ang1 = numpy.zeros((nPairs, 2, 1))
             #         dist1 = numpy.zeros((nPairs, 2, 1))
             #
             #         for j in range(nPairs):
             #             dist1[j,0,0] = dist[pairs[j][0]]
             #             dist1[j,1,0] = dist[pairs[j][1]]
             #             ang1[j,0,0] = ang[pairs[j][0]]
             #             ang1[j,1,0] = ang[pairs[j][1]]
             #
             #         return distx,disty, dist1,ang1
                 def __calculateVelVer(self, phase, lagTRange, _lambda):
                     Ts = lagTRange[1] - lagTRange[0]
                     velW = -_lambda*phase/(4*math.pi*Ts)
                     return velW
                 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
                     nPairs = tau1.shape[0]
                     nHeights = tau1.shape[1]
                     vel = numpy.zeros((nPairs,3,nHeights))
                     dist1 = numpy.reshape(dist, (dist.size,1))
                     angCos = numpy.cos(ang)
                     angSin = numpy.sin(ang)
                     vel0 = dist1*tau1/(2*tau2**2)
                     vel[:,0,:] = (vel0*angCos).sum(axis = 1)
                     vel[:,1,:] = (vel0*angSin).sum(axis = 1)
                     ind = numpy.where(numpy.isinf(vel))
                     vel[ind] = numpy.nan
                     return vel
             #     def __getPairsAutoCorr(self, pairsList, nChannels):
             #
             #         pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
             #
             #         for l in range(len(pairsList)):
             #             firstChannel = pairsList[l][0]
             #             secondChannel = pairsList[l][1]
             #
             #             #Obteniendo pares de Autocorrelacion
             #             if firstChannel == secondChannel:
             #                 pairsAutoCorr[firstChannel] = int(l)
             #
             #         pairsAutoCorr = pairsAutoCorr.astype(int)
             #
             #         pairsCrossCorr = range(len(pairsList))
             #         pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
             #
             #         return pairsAutoCorr, pairsCrossCorr
             #     def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
                 def techniqueSA(self, kwargs):
                     """
                     Function that implements Spaced Antenna (SA) technique.
                     Input:    Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
                                 Direction correction (if necessary), Ranges and SNR
                     Output:    Winds estimation (Zonal, Meridional and Vertical)
                     Parameters affected:    Winds
                     """
                     position_x = kwargs['positionX']
                     position_y = kwargs['positionY']
                     azimuth = kwargs['azimuth']
                     if kwargs.has_key('correctFactor'):
                         correctFactor = kwargs['correctFactor']
                     else:
                         correctFactor = 1
                     groupList = kwargs['groupList']
                     pairs_ccf = groupList[1]
                     tau = kwargs['tau']
                     _lambda = kwargs['_lambda']
                     #Cross Correlation pairs obtained
             #         pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairssList, nChannels)
             #         pairsArray = numpy.array(pairsList)[pairsCrossCorr]
             #         pairsSelArray = numpy.array(pairsSelected)
             #         pairs = []
             #
             #         #Wind estimation pairs obtained
             #         for i in range(pairsSelArray.shape[0]/2):
             #             ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
             #             ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
             #             pairs.append((ind1,ind2))
                     indtau = tau.shape[0]/2
                     tau1 = tau[:indtau,:]
                     tau2 = tau[indtau:-1,:]
             #         tau1 = tau1[pairs,:]
             #         tau2 = tau2[pairs,:]
                     phase1 = tau[-1,:]
                     #---------------------------------------------------------------------
                     #Metodo Directo
                     distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairs_ccf,azimuth)
                     winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
                     winds = stats.nanmean(winds, axis=0)
                     #---------------------------------------------------------------------
                     #Metodo General
             #         distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
             #         #Calculo Coeficientes de Funcion de Correlacion
             #         F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
             #         #Calculo de Velocidades
             #         winds = self.calculateVelUV(F,G,A,B,H)
                     #---------------------------------------------------------------------
                     winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
                     winds = correctFactor*winds
                     return winds
                 def __checkTime(self, currentTime, paramInterval, outputInterval):
                     dataTime = currentTime + paramInterval
                     deltaTime = dataTime - self.__initime
                     if deltaTime >= outputInterval or deltaTime < 0:
                         self.__dataReady = True
                     return
                 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax, binkm=2):
                     '''
                     Function that implements winds estimation technique with detected meteors.
                     Input:    Detected meteors, Minimum meteor quantity to wind estimation
                     Output:    Winds estimation (Zonal and Meridional)
                     Parameters affected:    Winds
                     '''
             #         print arrayMeteor.shape
                     #Settings
                     nInt = (heightMax - heightMin)/binkm
             #         print nInt
                     nInt = int(nInt)
             #         print nInt
                     winds = numpy.zeros((2,nInt))*numpy.nan
                     #Filter errors
                     error = numpy.where(arrayMeteor[:,-1] == 0)[0]
                     finalMeteor = arrayMeteor[error,:]
                     #Meteor Histogram
                     finalHeights = finalMeteor[:,2]
                     hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
                     nMeteorsPerI = hist[0]
                     heightPerI = hist[1]
                     #Sort of meteors
                     indSort = finalHeights.argsort()
                     finalMeteor2 = finalMeteor[indSort,:]
                     #    Calculating winds
                     ind1 = 0
                     ind2 = 0
                     for i in range(nInt):
                         nMet = nMeteorsPerI[i]
                         ind1 = ind2
                         ind2 = ind1 + nMet
                         meteorAux = finalMeteor2[ind1:ind2,:]
                         if meteorAux.shape[0] >= meteorThresh:
                             vel = meteorAux[:, 6]
                             zen = meteorAux[:, 4]*numpy.pi/180
                             azim = meteorAux[:, 3]*numpy.pi/180
                             n = numpy.cos(zen)
                     #         m = (1 - n**2)/(1 - numpy.tan(azim)**2)
                     #         l = m*numpy.tan(azim)
                             l = numpy.sin(zen)*numpy.sin(azim)
                             m = numpy.sin(zen)*numpy.cos(azim)
                             A = numpy.vstack((l, m)).transpose()
                             A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
                             windsAux = numpy.dot(A1, vel)
                             winds[0,i] = windsAux[0]
                             winds[1,i] = windsAux[1]
                     return winds, heightPerI[:-1]
                 def techniqueNSM_SA(self, **kwargs):
                     metArray = kwargs['metArray']
                     heightList = kwargs['heightList']
                     timeList = kwargs['timeList']
                     rx_location = kwargs['rx_location']
                     groupList = kwargs['groupList']
                     azimuth = kwargs['azimuth']
                     dfactor = kwargs['dfactor']
                     k = kwargs['k']
                     azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth)
                     d = dist*dfactor
                     #Phase calculation
                     metArray1 = self.__getPhaseSlope(metArray, heightList, timeList)
                     metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities
                     velEst = numpy.zeros((heightList.size,2))*numpy.nan
                     azimuth1 = azimuth1*numpy.pi/180
                     for i in range(heightList.size):
                         h = heightList[i]
                         indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0]
                         metHeight = metArray1[indH,:]
                         if metHeight.shape[0] >= 2:
                             velAux = numpy.asmatrix(metHeight[:,-2]).T    #Radial Velocities
                             iazim = metHeight[:,1].astype(int)
                             azimAux = numpy.asmatrix(azimuth1[iazim]).T    #Azimuths
                             A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux)))
                             A = numpy.asmatrix(A)
                             A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose()
                             velHor = numpy.dot(A1,velAux)
                             velEst[i,:] = numpy.squeeze(velHor)
                     return velEst
                 def __getPhaseSlope(self, metArray, heightList, timeList):
                     meteorList = []
                     #utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2
                     #Putting back together the meteor matrix
                     utctime = metArray[:,0]
                     uniqueTime = numpy.unique(utctime)
                     phaseDerThresh = 0.5
                     ippSeconds = timeList[1] - timeList[0]
                     sec = numpy.where(timeList>1)[0][0]
                     nPairs = metArray.shape[1] - 6
                     nHeights = len(heightList)
                     for t in uniqueTime:
                         metArray1 = metArray[utctime==t,:]
             #         phaseDerThresh = numpy.pi/4 #reducir Phase thresh
                         tmet = metArray1[:,1].astype(int)
                         hmet = metArray1[:,2].astype(int)
                         metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1))
                         metPhase[:,:] = numpy.nan
                         metPhase[:,hmet,tmet] = metArray1[:,6:].T
                         #Delete short trails
                         metBool = ~numpy.isnan(metPhase[0,:,:])
                         heightVect = numpy.sum(metBool, axis = 1)
                         metBool[heightVect<sec,:] = False
                         metPhase[:,heightVect<sec,:] = numpy.nan
                         #Derivative
                         metDer = numpy.abs(metPhase[:,:,1:] - metPhase[:,:,:-1])
                         phDerAux = numpy.dstack((numpy.full((nPairs,nHeights,1), False, dtype=bool),metDer > phaseDerThresh))
                         metPhase[phDerAux] = numpy.nan
                         #--------------------------METEOR DETECTION    -----------------------------------------
                         indMet = numpy.where(numpy.any(metBool,axis=1))[0]
                         for p in numpy.arange(nPairs):
                             phase = metPhase[p,:,:]
                             phDer = metDer[p,:,:]
                             for h in indMet:
                                 height = heightList[h]
                                 phase1 = phase[h,:] #82
                                 phDer1 = phDer[h,:]
                                 phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)])   #Unwrap
                                 indValid = numpy.where(~numpy.isnan(phase1))[0]
                                 initMet = indValid[0]
                                 endMet = 0
                                 for i in range(len(indValid)-1):
                                     #Time difference
                                     inow = indValid[i]
                                     inext = indValid[i+1]
                                     idiff = inext - inow
                                     #Phase difference
                                     phDiff = numpy.abs(phase1[inext] - phase1[inow])
                                     if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]:   #End of Meteor
                                         sizeTrail = inow - initMet + 1
                                         if sizeTrail>3*sec:  #Too short meteors
                                             x = numpy.arange(initMet,inow+1)*ippSeconds
                                             y = phase1[initMet:inow+1]
                                             ynnan = ~numpy.isnan(y)
                                             x = x[ynnan]
                                             y = y[ynnan]
                                             slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
                                             ylin = x*slope + intercept
                                             rsq = r_value**2
                                             if rsq > 0.5:
                                                 vel = slope#*height*1000/(k*d)
                                                 estAux = numpy.array([utctime,p,height, vel, rsq])
                                                 meteorList.append(estAux)
                                         initMet = inext
                     metArray2 = numpy.array(meteorList)
                     return metArray2
                 def __calculateAzimuth1(self, rx_location, pairslist, azimuth0):
                     azimuth1 = numpy.zeros(len(pairslist))
                     dist = numpy.zeros(len(pairslist))
                     for i in range(len(rx_location)):
                         ch0 = pairslist[i][0]
                         ch1 = pairslist[i][1]
                         diffX = rx_location[ch0][0] - rx_location[ch1][0]
                         diffY = rx_location[ch0][1] - rx_location[ch1][1]
                         azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi
                         dist[i] = numpy.sqrt(diffX**2 + diffY**2)
                     azimuth1 -= azimuth0
                     return azimuth1, dist
                 def techniqueNSM_DBS(self, **kwargs):
                     metArray = kwargs['metArray']
                     heightList = kwargs['heightList']
                     timeList = kwargs['timeList']
                     zenithList = kwargs['zenithList']
                     nChan = numpy.max(cmet) + 1
                     nHeights = len(heightList)
                     utctime = metArray[:,0]
                     cmet = metArray[:,1]
                     hmet = metArray1[:,3].astype(int)
                     h1met = heightList[hmet]*zenithList[cmet]
                     vmet = metArray1[:,5]
                     for i in range(nHeights - 1):
                         hmin = heightList[i]
                         hmax = heightList[i + 1]
                         vthisH = vmet[(h1met>=hmin) & (h1met<hmax)]
                     return data_output
                 def run(self, dataOut, technique, **kwargs):
                     param = dataOut.data_param
                     if dataOut.abscissaList != None:
                         absc = dataOut.abscissaList[:-1]
                     noise = dataOut.noise
                     heightList = dataOut.heightList
                     SNR = dataOut.data_SNR
                     if technique == 'DBS':
                         kwargs['velRadial'] = param[:,1,:] #Radial velocity
                         kwargs['heightList'] = heightList
                         kwargs['SNR'] = SNR
                         dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(kwargs) #DBS Function
                         dataOut.utctimeInit = dataOut.utctime
                         dataOut.outputInterval = dataOut.paramInterval
                     elif technique == 'SA':
                         #Parameters
             #             position_x = kwargs['positionX']
             #             position_y = kwargs['positionY']
             #             azimuth = kwargs['azimuth']
             #
             #             if kwargs.has_key('crosspairsList'):
             #                 pairs = kwargs['crosspairsList']
             #             else:
             #                 pairs = None
             #
             #             if kwargs.has_key('correctFactor'):
             #                 correctFactor = kwargs['correctFactor']
             #             else:
             #                 correctFactor = 1
             #             tau = dataOut.data_param
             #             _lambda = dataOut.C/dataOut.frequency
             #             pairsList = dataOut.groupList
             #             nChannels = dataOut.nChannels
                         kwargs['groupList'] = dataOut.groupList
                         kwargs['tau'] = dataOut.data_param
                         kwargs['_lambda'] = dataOut.C/dataOut.frequency
             #             dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
                         dataOut.data_output = self.techniqueSA(kwargs)
                         dataOut.utctimeInit = dataOut.utctime
                         dataOut.outputInterval = dataOut.timeInterval
                     elif technique == 'Meteors':
                         dataOut.flagNoData = True
                         self.__dataReady = False
                         if kwargs.has_key('nHours'):
                             nHours = kwargs['nHours']
                         else:
                             nHours = 1
                         if kwargs.has_key('meteorsPerBin'):
                             meteorThresh = kwargs['meteorsPerBin']
                         else:
                             meteorThresh = 6
                         if kwargs.has_key('hmin'):
                             hmin = kwargs['hmin']
                         else:   hmin = 70
                         if kwargs.has_key('hmax'):
                             hmax = kwargs['hmax']
                         else:   hmax = 110
                         if kwargs.has_key('BinKm'):
                             binkm = kwargs['BinKm']
                         else:
                             binkm = 2
                         dataOut.outputInterval = nHours*3600
                         if self.__isConfig == False:
             #                 self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
                             #Get Initial LTC time
                             self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
                             self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
                             self.__isConfig = True
                         if self.__buffer is None:
                             self.__buffer = dataOut.data_param
                             self.__firstdata = copy.copy(dataOut)
                         else:
                             self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
                         self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
                         if self.__dataReady:
                             dataOut.utctimeInit = self.__initime
                             self.__initime += dataOut.outputInterval #to erase time offset
                             dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax, binkm)
                             dataOut.flagNoData = False
                             self.__buffer = None
                     elif technique == 'Meteors1':
                         dataOut.flagNoData = True
                         self.__dataReady = False
                         if kwargs.has_key('nMins'):
                             nMins = kwargs['nMins']
                         else: nMins = 20
                         if kwargs.has_key('rx_location'):
                             rx_location = kwargs['rx_location']
                         else: rx_location = [(0,1),(1,1),(1,0)]
                         if kwargs.has_key('azimuth'):
                             azimuth = kwargs['azimuth']
                         else: azimuth = 51
                         if kwargs.has_key('dfactor'):
                             dfactor = kwargs['dfactor']
                         if kwargs.has_key('mode'):
                             mode = kwargs['mode']
                         else: mode = 'SA'
                         #Borrar luego esto
                         if dataOut.groupList is None:
                             dataOut.groupList = [(0,1),(0,2),(1,2)]
                         groupList = dataOut.groupList
                         C = 3e8
                         freq = 50e6
                         lamb = C/freq
                         k = 2*numpy.pi/lamb
                         timeList = dataOut.abscissaList
                         heightList = dataOut.heightList
                         if self.__isConfig == False:
                             dataOut.outputInterval = nMins*60
             #                 self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
                             #Get Initial LTC time
                             initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
                             minuteAux = initime.minute
                             minuteNew = int(numpy.floor(minuteAux/nMins)*nMins)
                             self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
                             self.__isConfig = True
                         if self.__buffer is None:
                             self.__buffer = dataOut.data_param
                             self.__firstdata = copy.copy(dataOut)
                         else:
                             self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
                         self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
                         if self.__dataReady:
                             dataOut.utctimeInit = self.__initime
                             self.__initime += dataOut.outputInterval #to erase time offset
                             metArray = self.__buffer
                             if mode == 'SA':
                                 dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)
                             elif mode == 'DBS':
                                 dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList)
                             dataOut.data_output = dataOut.data_output.T
                             dataOut.flagNoData = False
                             self.__buffer = None
                     return
             class EWDriftsEstimation(Operation):
                 def __correctValues(self, heiRang, phi, velRadial, SNR):
                     listPhi = phi.tolist()
                     maxid = listPhi.index(max(listPhi))
                     minid = listPhi.index(min(listPhi))
                     rango = range(len(phi))
                #     rango = numpy.delete(rango,maxid)
                     heiRang1 = heiRang*math.cos(phi[maxid])
                     heiRangAux = heiRang*math.cos(phi[minid])
                     indOut = (heiRang1 < heiRangAux[0]).nonzero()
                     heiRang1 = numpy.delete(heiRang1,indOut)
                     velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
                     SNR1 = numpy.zeros([len(phi),len(heiRang1)])
                     for i in rango:
                         x = heiRang*math.cos(phi[i])
                         y1 = velRadial[i,:]
                         f1 = interpolate.interp1d(x,y1,kind = 'cubic')
                         x1 = heiRang1
                         y11 = f1(x1)
                         y2 = SNR[i,:]
                         f2 = interpolate.interp1d(x,y2,kind = 'cubic')
                         y21 = f2(x1)
                         velRadial1[i,:] = y11
                         SNR1[i,:] = y21
                     return heiRang1, velRadial1, SNR1
                 def run(self, dataOut, zenith, zenithCorrection):
                     heiRang = dataOut.heightList
                     velRadial = dataOut.data_param[:,3,:]
                     SNR = dataOut.data_SNR
                     zenith = numpy.array(zenith)
                     zenith -= zenithCorrection
                     zenith *= numpy.pi/180
                     heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
                     alp = zenith[0]
                     bet = zenith[1]
                     w_w = velRadial1[0,:]
                     w_e = velRadial1[1,:]
                     w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
                     u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
                     winds = numpy.vstack((u,w))
                     dataOut.heightList = heiRang1
                     dataOut.data_output = winds
                     dataOut.data_SNR = SNR1
                     dataOut.utctimeInit = dataOut.utctime
                     dataOut.outputInterval = dataOut.timeInterval
                     return
             #---------------    Non Specular Meteor    ----------------
             class NonSpecularMeteorDetection(Operation):
                 def run(self, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False):
                     data_acf = self.dataOut.data_pre[0]
                     data_ccf = self.dataOut.data_pre[1]
                     lamb = self.dataOut.C/self.dataOut.frequency
                     tSamp = self.dataOut.ippSeconds*self.dataOut.nCohInt
                     paramInterval = self.dataOut.paramInterval
                     nChannels = data_acf.shape[0]
                     nLags = data_acf.shape[1]
                     nProfiles = data_acf.shape[2]
                     nHeights = self.dataOut.nHeights
                     nCohInt = self.dataOut.nCohInt
                     sec = numpy.round(nProfiles/self.dataOut.paramInterval)
                     heightList = self.dataOut.heightList
                     ippSeconds = self.dataOut.ippSeconds*self.dataOut.nCohInt*self.dataOut.nAvg
                     utctime = self.dataOut.utctime
                     self.dataOut.abscissaList = numpy.arange(0,paramInterval+ippSeconds,ippSeconds)
                     #------------------------    SNR    --------------------------------------
                     power = data_acf[:,0,:,:].real
                     noise = numpy.zeros(nChannels)
                     SNR = numpy.zeros(power.shape)
                     for i in range(nChannels):
                         noise[i] = hildebrand_sekhon(power[i,:], nCohInt)
                         SNR[i] = (power[i]-noise[i])/noise[i]
                     SNRm = numpy.nanmean(SNR, axis = 0)
                     SNRdB = 10*numpy.log10(SNR)
                     if mode == 'SA':
                         nPairs = data_ccf.shape[0]
                         #----------------------    Coherence and Phase   --------------------------
                         phase = numpy.zeros(data_ccf[:,0,:,:].shape)
             #             phase1 = numpy.copy(phase)
                         coh1 = numpy.zeros(data_ccf[:,0,:,:].shape)
                         for p in range(nPairs):
                             ch0 = self.dataOut.groupList[p][0]
                             ch1 = self.dataOut.groupList[p][1]
                             ccf = data_ccf[p,0,:,:]/numpy.sqrt(data_acf[ch0,0,:,:]*data_acf[ch1,0,:,:])
                             phase[p,:,:] = ndimage.median_filter(numpy.angle(ccf), size = (5,1)) #median filter
             #                 phase1[p,:,:] = numpy.angle(ccf) #median filter
                             coh1[p,:,:] = ndimage.median_filter(numpy.abs(ccf), 5) #median filter
             #                 coh1[p,:,:] = numpy.abs(ccf) #median filter
                         coh = numpy.nanmax(coh1, axis = 0)
             #             struc = numpy.ones((5,1))
             #             coh = ndimage.morphology.grey_dilation(coh, size=(10,1))
                         #----------------------    Radial Velocity    ----------------------------
                         phaseAux = numpy.mean(numpy.angle(data_acf[:,1,:,:]), axis = 0)
                         velRad = phaseAux*lamb/(4*numpy.pi*tSamp)
                         if allData:
                             boolMetFin = ~numpy.isnan(SNRm)
             #                 coh[:-1,:] = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
                         else:
                             #------------------------    Meteor mask    ---------------------------------
             #                 #SNR mask
             #                 boolMet = (SNRdB>SNRthresh)#|(~numpy.isnan(SNRdB))
             #
             #                 #Erase small objects
             #                 boolMet1 = self.__erase_small(boolMet, 2*sec, 5)
             #
             #                 auxEEJ = numpy.sum(boolMet1,axis=0)
             #                 indOver = auxEEJ>nProfiles*0.8  #Use this later
             #                 indEEJ = numpy.where(indOver)[0]
             #                 indNEEJ = numpy.where(~indOver)[0]
             #
             #                 boolMetFin = boolMet1
             #
             #                 if indEEJ.size > 0:
             #                     boolMet1[:,indEEJ] = False  #Erase heights with EEJ
             #
             #                     boolMet2 = coh > cohThresh
             #                     boolMet2 = self.__erase_small(boolMet2, 2*sec,5)
             #
             #                     #Final Meteor mask
             #                     boolMetFin = boolMet1|boolMet2
                             #Coherence mask
                             boolMet1 = coh > 0.75
                             struc = numpy.ones((30,1))
                             boolMet1 = ndimage.morphology.binary_dilation(boolMet1, structure=struc)
                             #Derivative mask
                             derPhase = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
                             boolMet2 = derPhase < 0.2
             #                 boolMet2 = ndimage.morphology.binary_opening(boolMet2)
             #                 boolMet2 = ndimage.morphology.binary_closing(boolMet2, structure = numpy.ones((10,1)))
                             boolMet2 = ndimage.median_filter(boolMet2,size=5)
                             boolMet2 = numpy.vstack((boolMet2,numpy.full((1,nHeights), True, dtype=bool)))
             #                 #Final mask
             #                 boolMetFin = boolMet2
                             boolMetFin = boolMet1&boolMet2
             #                 boolMetFin = ndimage.morphology.binary_dilation(boolMetFin)
                         #Creating data_param
                         coordMet = numpy.where(boolMetFin)
                         tmet = coordMet[0]
                         hmet = coordMet[1]
                         data_param = numpy.zeros((tmet.size, 6 + nPairs))
                         data_param[:,0] = utctime
                         data_param[:,1] = tmet
                         data_param[:,2] = hmet
                         data_param[:,3] = SNRm[tmet,hmet]
                         data_param[:,4] = velRad[tmet,hmet]
                         data_param[:,5] = coh[tmet,hmet]
                         data_param[:,6:] = phase[:,tmet,hmet].T
                     elif mode == 'DBS':
                         self.dataOut.groupList = numpy.arange(nChannels)
                         #Radial Velocities
             #             phase = numpy.angle(data_acf[:,1,:,:])
                         phase = ndimage.median_filter(numpy.angle(data_acf[:,1,:,:]), size = (1,5,1))
                         velRad = phase*lamb/(4*numpy.pi*tSamp)
                         #Spectral width
                         acf1 = ndimage.median_filter(numpy.abs(data_acf[:,1,:,:]), size = (1,5,1))
                         acf2 = ndimage.median_filter(numpy.abs(data_acf[:,2,:,:]), size = (1,5,1))
                         spcWidth = (lamb/(2*numpy.sqrt(6)*numpy.pi*tSamp))*numpy.sqrt(numpy.log(acf1/acf2))
             #             velRad = ndimage.median_filter(velRad, size = (1,5,1))
                         if allData:
                             boolMetFin = ~numpy.isnan(SNRdB)
                         else:
                             #SNR
                             boolMet1 = (SNRdB>SNRthresh) #SNR mask
                             boolMet1 = ndimage.median_filter(boolMet1, size=(1,5,5))
                             #Radial velocity
                             boolMet2 = numpy.abs(velRad) < 30
                             boolMet2 = ndimage.median_filter(boolMet2, (1,5,5))
                             #Spectral Width
                             boolMet3 = spcWidth < 30
                             boolMet3 = ndimage.median_filter(boolMet3, (1,5,5))
             #                 boolMetFin = self.__erase_small(boolMet1, 10,5)
                             boolMetFin = boolMet1&boolMet2&boolMet3
                         #Creating data_param
                         coordMet = numpy.where(boolMetFin)
                         cmet = coordMet[0]
                         tmet = coordMet[1]
                         hmet = coordMet[2]
                         data_param = numpy.zeros((tmet.size, 7))
                         data_param[:,0] = utctime
                         data_param[:,1] = cmet
                         data_param[:,2] = tmet
                         data_param[:,3] = hmet
                         data_param[:,4] = SNR[cmet,tmet,hmet].T
                         data_param[:,5] = velRad[cmet,tmet,hmet].T
                         data_param[:,6] = spcWidth[cmet,tmet,hmet].T
             #         self.dataOut.data_param = data_int
                     if len(data_param) == 0:
                         self.dataOut.flagNoData = True
                     else:
                         self.dataOut.data_param = data_param
                 def __erase_small(self, binArray, threshX, threshY):
                     labarray, numfeat = ndimage.measurements.label(binArray)
                     binArray1 = numpy.copy(binArray)
                     for i in range(1,numfeat + 1):
                         auxBin = (labarray==i)
                         auxSize = auxBin.sum()
                         x,y = numpy.where(auxBin)
                         widthX = x.max() - x.min()
                         widthY = y.max() - y.min()
                         #width X: 3 seg -> 12.5*3
                         #width Y:
                         if (auxSize < 50) or (widthX < threshX) or (widthY < threshY):
                             binArray1[auxBin] = False
                     return binArray1
             #---------------    Specular Meteor    ----------------
             class SMDetection(Operation):
                 '''
                     Function DetectMeteors()
                         Project developed with paper:
                         HOLDSWORTH ET AL. 2004
                     Input:
                         self.dataOut.data_pre
                         centerReceiverIndex:      From the channels, which is the center receiver
                         hei_ref:                  Height reference for the Beacon signal extraction
                         tauindex:
                         predefinedPhaseShifts:    Predefined phase offset for the voltge signals
                         cohDetection:             Whether to user Coherent detection or not
                         cohDet_timeStep:          Coherent Detection calculation time step
                         cohDet_thresh:            Coherent Detection phase threshold to correct phases
                         noise_timeStep:           Noise calculation time step
                         noise_multiple:           Noise multiple to define signal threshold
                         multDet_timeLimit:        Multiple Detection Removal time limit in seconds
                         multDet_rangeLimit:       Multiple Detection Removal range limit in km
                         phaseThresh:              Maximum phase difference between receiver to be consider a meteor
                         SNRThresh:                Minimum SNR threshold of the meteor signal to be consider a meteor
                         hmin:                     Minimum Height of the meteor to use it in the further wind estimations
                         hmax:                     Maximum Height of the meteor to use it in the further wind estimations
                         azimuth:                  Azimuth angle correction
                     Affected:
                         self.dataOut.data_param
                     Rejection Criteria (Errors):
 : No error; analysis OK
 : SNR < SNR threshold
 : angle of arrival (AOA) ambiguously determined
 : AOA estimate not feasible
 : Large difference in AOAs obtained from different antenna baselines
 : echo at start or end of time series
 : echo less than 5 examples long; too short for analysis
 : echo rise exceeds 0.3s
 : echo decay time less than twice rise time
 : large power level before echo
 : large power level after echo
 : poor fit to amplitude for estimation of decay time
 : poor fit to CCF phase variation for estimation of radial drift velocity
 : height unresolvable echo: not valid height within 70 to 110 km
 : height ambiguous echo: more then one possible height within 70 to 110 km
 : radial drift velocity or projected horizontal velocity exceeds 200 m/s
 : oscilatory echo, indicating event most likely not an underdense echo
 : phase difference in meteor Reestimation
                     Data Storage:
                         Meteors for Wind Estimation   (8):
                         Utc Time   |    Range    Height
                         Azimuth    Zenith    errorCosDir
                         VelRad    errorVelRad
                         Phase0 Phase1 Phase2 Phase3
                         TypeError
                      '''
                 def run(self, dataOut, hei_ref = None, tauindex = 0,
                                   phaseOffsets = None,
                                   cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
                                   noise_timeStep = 4, noise_multiple = 4,
                                   multDet_timeLimit = 1, multDet_rangeLimit = 3,
                                   phaseThresh = 20, SNRThresh = 5,
                                   hmin = 50, hmax=150, azimuth = 0,
                                   channelPositions = None) :
                     #Getting Pairslist
                     if channelPositions is None:
             #             channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)]   #T
                         channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)]   #Estrella
                     meteorOps = SMOperations()
                     pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
                     heiRang = dataOut.getHeiRange()
                     #Get Beacon signal - No Beacon signal anymore
             #         newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
             #
             #         if hei_ref != None:
             #             newheis = numpy.where(self.dataOut.heightList>hei_ref)
             #
                     #****************REMOVING HARDWARE PHASE DIFFERENCES***************
                     # see if the user put in pre defined phase shifts
                     voltsPShift = dataOut.data_pre.copy()
             #         if predefinedPhaseShifts != None:
             #             hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
             #
             # #         elif beaconPhaseShifts:
             # #             #get hardware phase shifts using beacon signal
             # #             hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
             # #             hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
             #
             #         else:
             #             hardwarePhaseShifts = numpy.zeros(5)
             #
             #         voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
             #         for i in range(self.dataOut.data_pre.shape[0]):
             #             voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
                     #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
                     #Remove DC
                     voltsDC = numpy.mean(voltsPShift,1)
                     voltsDC = numpy.mean(voltsDC,1)
                     for i in range(voltsDC.shape[0]):
                         voltsPShift[i] = voltsPShift[i] - voltsDC[i]
                     #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
             #         voltsPShift = voltsPShift[:,:,:newheis[0][0]]
                     #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
                     #Coherent Detection
                     if cohDetection:
                         #use coherent detection to get the net power
                         cohDet_thresh = cohDet_thresh*numpy.pi/180
                         voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, dataOut.timeInterval, pairslist0, cohDet_thresh)
                     #Non-coherent detection!
                     powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
                     #********** END OF COH/NON-COH POWER CALCULATION**********************
                     #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
                     #Get noise
                     noise, noise1 = self.__getNoise(powerNet, noise_timeStep, dataOut.timeInterval)
             #         noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
                     #Get signal threshold
                     signalThresh = noise_multiple*noise
                     #Meteor echoes detection
                     listMeteors = self.__findMeteors(powerNet, signalThresh)
                     #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
                     #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
                     #Parameters
                     heiRange = dataOut.getHeiRange()
                     rangeInterval = heiRange[1] - heiRange[0]
                     rangeLimit = multDet_rangeLimit/rangeInterval
                     timeLimit = multDet_timeLimit/dataOut.timeInterval
                     #Multiple detection removals
                     listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
                     #************ END OF REMOVE MULTIPLE DETECTIONS **********************
                     #*********************     METEOR REESTIMATION  (3.7, 3.8, 3.9, 3.10)   ********************
                     #Parameters
                     phaseThresh = phaseThresh*numpy.pi/180
                     thresh = [phaseThresh, noise_multiple, SNRThresh]
                     #Meteor reestimation  (Errors N 1, 6, 12, 17)
                     listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist0, thresh, noise, dataOut.timeInterval, dataOut.frequency)
             #         listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
                     #Estimation of decay times (Errors N 7, 8, 11)
                     listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, dataOut.timeInterval, dataOut.frequency)
                     #*******************     END OF METEOR REESTIMATION    *******************
                     #*********************    METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13)    **************************
                     #Calculating Radial Velocity (Error N 15)
                     radialStdThresh = 10
                     listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, dataOut.timeInterval)
                     if len(listMeteors4) > 0:
                         #Setting New Array
                         date = dataOut.utctime
                         arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
                         #Correcting phase offset
                         if phaseOffsets != None:
                             phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
                             arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
                         #Second Pairslist
                         pairsList = []
                         pairx = (0,1)
                         pairy = (2,3)
                         pairsList.append(pairx)
                         pairsList.append(pairy)
                         jph = numpy.array([0,0,0,0])
                         h = (hmin,hmax)
                         arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
             #             #Calculate AOA (Error N 3, 4)
             #             #JONES ET AL. 1998
             #             error = arrayParameters[:,-1]
             #             AOAthresh = numpy.pi/8
             #             phases = -arrayParameters[:,9:13]
             #             arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)
             #
             #             #Calculate Heights (Error N 13 and 14)
             #             error = arrayParameters[:,-1]
             #             Ranges = arrayParameters[:,2]
             #             zenith = arrayParameters[:,5]
             #             arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax)
             #             error = arrayParameters[:,-1]
                     #*********************    END OF PARAMETERS CALCULATION    **************************
                     #***************************+     PASS DATA TO NEXT STEP    **********************
             #             arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
                         dataOut.data_param = arrayParameters
                         if arrayParameters is None:
                             dataOut.flagNoData = True
                     else:
                         dataOut.flagNoData = True
                     return
                 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
                     minIndex = min(newheis[0])
                     maxIndex = max(newheis[0])
                     voltage = voltage0[:,:,minIndex:maxIndex+1]
                     nLength = voltage.shape[1]/n
                     nMin = 0
                     nMax = 0
                     phaseOffset = numpy.zeros((len(pairslist),n))
                     for i in range(n):
                         nMax += nLength
                         phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
                         phaseCCF = numpy.mean(phaseCCF, axis = 2)
                         phaseOffset[:,i] = phaseCCF.transpose()
                         nMin = nMax
             #         phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
                     #Remove Outliers
                     factor = 2
                     wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
                     dw = numpy.std(wt,axis = 1)
                     dw = dw.reshape((dw.size,1))
                     ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
                     phaseOffset[ind] = numpy.nan
                     phaseOffset = stats.nanmean(phaseOffset, axis=1)
                     return phaseOffset
                 def __shiftPhase(self, data, phaseShift):
                     #this will shift the phase of a complex number
                     dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
                     return dataShifted
                 def __estimatePhaseDifference(self, array, pairslist):
                     nChannel = array.shape[0]
                     nHeights = array.shape[2]
                     numPairs = len(pairslist)
             #         phaseCCF = numpy.zeros((nChannel, 5, nHeights))
                     phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
                     #Correct phases
                     derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
                     indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
                     if indDer[0].shape[0] > 0:
                         for i in range(indDer[0].shape[0]):
                             signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
                             phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
             #         for j in range(numSides):
             #             phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
             #             phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
             #
                     #Linear
                     phaseInt = numpy.zeros((numPairs,1))
                     angAllCCF = phaseCCF[:,[0,1,3,4],0]
                     for j in range(numPairs):
                         fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
                         phaseInt[j] = fit[1]
                     #Phase Differences
                     phaseDiff = phaseInt - phaseCCF[:,2,:]
                     phaseArrival = phaseInt.reshape(phaseInt.size)
                     #Dealias
                     phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival))
             #         indAlias = numpy.where(phaseArrival > numpy.pi)
             #         phaseArrival[indAlias] -= 2*numpy.pi
             #         indAlias = numpy.where(phaseArrival < -numpy.pi)
             #         phaseArrival[indAlias] += 2*numpy.pi
                     return phaseDiff, phaseArrival
                 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
                     #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
                     #find the phase shifts of each channel over 1 second intervals
                     #only look at ranges below the beacon signal
                     numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
                     numBlocks = int(volts.shape[1]/numProfPerBlock)
                     numHeights = volts.shape[2]
                     nChannel = volts.shape[0]
                     voltsCohDet = volts.copy()
                     pairsarray = numpy.array(pairslist)
                     indSides = pairsarray[:,1]
             #         indSides = numpy.array(range(nChannel))
             #         indSides = numpy.delete(indSides, indCenter)
             #
             #         listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
                     listBlocks = numpy.array_split(volts, numBlocks, 1)
                     startInd = 0
                     endInd = 0
                     for i in range(numBlocks):
                         startInd = endInd
                         endInd = endInd + listBlocks[i].shape[1]
                         arrayBlock = listBlocks[i]
             #             arrayBlockCenter = listCenter[i]
                         #Estimate the Phase Difference
                         phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
                         #Phase Difference RMS
                         arrayPhaseRMS = numpy.abs(phaseDiff)
                         phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
                         indPhase = numpy.where(phaseRMSaux==4)
                         #Shifting
                         if indPhase[0].shape[0] > 0:
                             for j in range(indSides.size):
                                 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
                             voltsCohDet[:,startInd:endInd,:] = arrayBlock
                     return voltsCohDet
                 def __calculateCCF(self, volts, pairslist ,laglist):
                     nHeights = volts.shape[2]
                     nPoints = volts.shape[1]
                     voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
                     for i in range(len(pairslist)):
                         volts1 = volts[pairslist[i][0]]
                         volts2 = volts[pairslist[i][1]]
                         for t in range(len(laglist)):
                             idxT = laglist[t]
                             if idxT >= 0:
                                 vStacked = numpy.vstack((volts2[idxT:,:],
                                                        numpy.zeros((idxT, nHeights),dtype='complex')))
                             else:
                                 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
                                                        volts2[:(nPoints + idxT),:]))
                             voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
                             vStacked = None
                     return voltsCCF
                 def __getNoise(self, power, timeSegment, timeInterval):
                     numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
                     numBlocks = int(power.shape[0]/numProfPerBlock)
                     numHeights = power.shape[1]
                     listPower = numpy.array_split(power, numBlocks, 0)
                     noise = numpy.zeros((power.shape[0], power.shape[1]))
                     noise1 = numpy.zeros((power.shape[0], power.shape[1]))
                     startInd = 0
                     endInd = 0
                     for i in range(numBlocks):             #split por canal
                         startInd = endInd
                         endInd = endInd + listPower[i].shape[0]
                         arrayBlock = listPower[i]
                         noiseAux = numpy.mean(arrayBlock, 0)
             #             noiseAux = numpy.median(noiseAux)
             #             noiseAux = numpy.mean(arrayBlock)
                         noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
                         noiseAux1 = numpy.mean(arrayBlock)
                         noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
                     return noise, noise1
                 def __findMeteors(self, power, thresh):
                     nProf = power.shape[0]
                     nHeights = power.shape[1]
                     listMeteors = []
                     for i in range(nHeights):
                         powerAux = power[:,i]
                         threshAux = thresh[:,i]
                         indUPthresh = numpy.where(powerAux > threshAux)[0]
                         indDNthresh = numpy.where(powerAux <= threshAux)[0]
                         j = 0
                         while (j < indUPthresh.size - 2):
                             if (indUPthresh[j + 2] == indUPthresh[j] + 2):
                                 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
                                 indDNthresh = indDNthresh[indDNAux]
                                 if (indDNthresh.size > 0):
                                     indEnd = indDNthresh[0] - 1
                                     indInit = indUPthresh[j] if isinstance(indUPthresh[j], (int, float)) else indUPthresh[j][0] ##CHECK!!!!
                                     meteor = powerAux[indInit:indEnd + 1]
                                     indPeak = meteor.argmax() + indInit
                                     FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
                                     listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
                                     j = numpy.where(indUPthresh == indEnd)[0] + 1
                                 else: j+=1
                             else: j+=1
                     return listMeteors
                 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
                     arrayMeteors = numpy.asarray(listMeteors)
                     listMeteors1 = []
                     while arrayMeteors.shape[0] > 0:
                         FLAs = arrayMeteors[:,4]
                         maxFLA = FLAs.argmax()
                         listMeteors1.append(arrayMeteors[maxFLA,:])
                         MeteorInitTime = arrayMeteors[maxFLA,1]
                         MeteorEndTime = arrayMeteors[maxFLA,3]
                         MeteorHeight = arrayMeteors[maxFLA,0]
                         #Check neighborhood
                         maxHeightIndex = MeteorHeight + rangeLimit
                         minHeightIndex = MeteorHeight - rangeLimit
                         minTimeIndex = MeteorInitTime - timeLimit
                         maxTimeIndex = MeteorEndTime + timeLimit
                         #Check Heights
                         indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
                         indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
                         indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
                         arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
                     return listMeteors1
                 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
                     numHeights = volts.shape[2]
                     nChannel = volts.shape[0]
                     thresholdPhase = thresh[0]
                     thresholdNoise = thresh[1]
                     thresholdDB = float(thresh[2])
                     thresholdDB1 = 10**(thresholdDB/10)
                     pairsarray = numpy.array(pairslist)
                     indSides = pairsarray[:,1]
                     pairslist1 = list(pairslist)
                     pairslist1.append((0,4))
                     pairslist1.append((1,3))
                     listMeteors1 = []
                     listPowerSeries = []
                     listVoltageSeries = []
                     #volts has the war data
                     if frequency == 30.175e6:
                         timeLag = 45*10**-3
                     else:
                         timeLag = 15*10**-3
                     lag = int(numpy.ceil(timeLag/timeInterval))
                     for i in range(len(listMeteors)):
                         ######################   3.6 - 3.7 PARAMETERS REESTIMATION    #########################
                         meteorAux = numpy.zeros(16)
                         #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
                         mHeight = int(listMeteors[i][0])
                         mStart = int(listMeteors[i][1])
                         mPeak = int(listMeteors[i][2])
                         mEnd = int(listMeteors[i][3])
                         #get the volt data between the start and end times of the meteor
                         meteorVolts = volts[:,mStart:mEnd+1,mHeight]
                         meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
                         #3.6. Phase Difference estimation
                         phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
                         #3.7. Phase difference removal & meteor start, peak and end times reestimated
                         #meteorVolts0.- all Channels, all Profiles
                         meteorVolts0 = volts[:,:,mHeight]
                         meteorThresh = noise[:,mHeight]*thresholdNoise
                         meteorNoise = noise[:,mHeight]
                         meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
                         powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0)  #Power
                         #Times reestimation
                         mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
                         if mStart1.size > 0:
                             mStart1 = mStart1[-1] + 1
                         else:
                             mStart1 = mPeak
                         mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
                         mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
                         if mEndDecayTime1.size == 0:
                             mEndDecayTime1 = powerNet0.size
                         else:
                             mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
             #             mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
                         #meteorVolts1.- all Channels, from start to end
                         meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
                         meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
                         if meteorVolts2.shape[1] == 0:
                             meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
                         meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
                         meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
                         #####################    END PARAMETERS REESTIMATION    #########################
                         #####################   3.8 PHASE DIFFERENCE REESTIMATION  ########################
             #             if mEnd1 - mStart1 > 4:       #Error Number 6: echo less than 5 samples long; too short for analysis
                         if meteorVolts2.shape[1] > 0:
                             #Phase Difference re-estimation
                             phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1)   #Phase Difference Estimation
             #                 phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
                             meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
                             phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
                             meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4])     #Phase Shifting
                             #Phase Difference RMS
                             phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
                             powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
                             #Data from Meteor
                             mPeak1 = powerNet1.argmax() + mStart1
                             mPeakPower1 = powerNet1.max()
                             noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
                             mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
                             Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
                             Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
                             PowerSeries  = powerNet0[mStart1:mEndDecayTime1 + 1]
                             #Vectorize
                             meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
                             meteorAux[7:11] = phaseDiffint[0:4]
                             #Rejection Criterions
                             if phaseRMS1 > thresholdPhase:  #Error Number 17: Phase variation
                                 meteorAux[-1] = 17
                             elif mSNR1 < thresholdDB1:      #Error Number 1: SNR < threshold dB
                                 meteorAux[-1] = 1
                         else:
                             meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
                             meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
                             PowerSeries = 0
                         listMeteors1.append(meteorAux)
                         listPowerSeries.append(PowerSeries)
                         listVoltageSeries.append(meteorVolts1)
                     return listMeteors1, listPowerSeries, listVoltageSeries
                 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
                     threshError = 10
                     #Depending if it is 30 or 50 MHz
                     if frequency == 30.175e6:
                         timeLag = 45*10**-3
                     else:
                         timeLag = 15*10**-3
                     lag = int(numpy.ceil(timeLag/timeInterval))
                     listMeteors1 = []
                     for i in range(len(listMeteors)):
                         meteorPower = listPower[i]
                         meteorAux = listMeteors[i]
                         if meteorAux[-1] == 0:
                             try:
                                 indmax = meteorPower.argmax()
                                 indlag = indmax + lag
                                 y = meteorPower[indlag:]
                                 x = numpy.arange(0, y.size)*timeLag
                                 #first guess
                                 a = y[0]
                                 tau = timeLag
                                 #exponential fit
                                 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
                                 y1 = self.__exponential_function(x, *popt)
                                 #error estimation
                                 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
                                 decayTime = popt[1]
                                 riseTime = indmax*timeInterval
                                 meteorAux[11:13] = [decayTime, error]
                                 #Table items 7, 8 and 11
                                 if (riseTime > 0.3):            #Number 7: Echo rise exceeds 0.3s
                                     meteorAux[-1] = 7
                                 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
                                     meteorAux[-1] = 8
                                 if (error > threshError):       #Number 11: Poor fit to amplitude for estimation of decay time
                                     meteorAux[-1] = 11
                             except:
                                 meteorAux[-1] = 11
                         listMeteors1.append(meteorAux)
                     return listMeteors1
                 #Exponential Function
                 def __exponential_function(self, x, a, tau):
                     y = a*numpy.exp(-x/tau)
                     return y
                 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist,  timeInterval):
                     pairslist1 = list(pairslist)
                     pairslist1.append((0,4))
                     pairslist1.append((1,3))
                     numPairs = len(pairslist1)
                     #Time Lag
                     timeLag = 45*10**-3
                     c = 3e8
                     lag = numpy.ceil(timeLag/timeInterval)
                     freq = 30.175e6
                     listMeteors1 = []
                     for i in range(len(listMeteors)):
                         meteorAux = listMeteors[i]
                         if meteorAux[-1] == 0:
                             mStart = listMeteors[i][1]
                             mPeak = listMeteors[i][2]
                             mLag = mPeak - mStart + lag
                             #get the volt data between the start and end times of the meteor
                             meteorVolts = listVolts[i]
                             meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
                             #Get CCF
                             allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
                             #Method 2
                             slopes = numpy.zeros(numPairs)
                             time = numpy.array([-2,-1,1,2])*timeInterval
                             angAllCCF = numpy.angle(allCCFs[:,[0,4,2,3],0])
                             #Correct phases
                             derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
                             indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
                             if indDer[0].shape[0] > 0:
                                 for i in range(indDer[0].shape[0]):
                                     signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
                                     angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
             #                     fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
                             for j in range(numPairs):
                                 fit = stats.linregress(time, angAllCCF[j,:])
                                 slopes[j] = fit[0]
                             #Remove Outlier
             #                 indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
             #                 slopes = numpy.delete(slopes,indOut)
             #                 indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
             #                 slopes = numpy.delete(slopes,indOut)
                             radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
                             radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
                             meteorAux[-2] = radialError
                             meteorAux[-3] = radialVelocity
                             #Setting Error
                             #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
                             if numpy.abs(radialVelocity) > 200:
                                 meteorAux[-1] = 15
                             #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
                             elif radialError > radialStdThresh:
                                 meteorAux[-1] = 12
                         listMeteors1.append(meteorAux)
                     return listMeteors1
                 def __setNewArrays(self, listMeteors, date, heiRang):
                     #New arrays
                     arrayMeteors = numpy.array(listMeteors)
                     arrayParameters = numpy.zeros((len(listMeteors), 13))
                     #Date inclusion
             #         date = re.findall(r'\((.*?)\)', date)
             #         date = date[0].split(',')
             #         date = map(int, date)
             #
             #         if len(date)<6:
             #             date.append(0)
             #
             #         date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
             #         arrayDate = numpy.tile(date, (len(listMeteors), 1))
                     arrayDate = numpy.tile(date, (len(listMeteors)))
                     #Meteor array
             #         arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
             #         arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
                     #Parameters Array
                     arrayParameters[:,0] = arrayDate #Date
                     arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range
                     arrayParameters[:,6:8] = arrayMeteors[:,-3:-1] #Radial velocity and its error
                     arrayParameters[:,8:12] = arrayMeteors[:,7:11]  #Phases
                     arrayParameters[:,-1] = arrayMeteors[:,-1]  #Error
                     return arrayParameters
             class CorrectSMPhases(Operation):
                 def run(self, dataOut, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None):
                     arrayParameters = dataOut.data_param
                     pairsList = []
                     pairx = (0,1)
                     pairy = (2,3)
                     pairsList.append(pairx)
                     pairsList.append(pairy)
                     jph = numpy.zeros(4)
                     phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
                 #         arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
                     arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))
                     meteorOps = SMOperations()
                     if channelPositions is None:
                 #             channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)]   #T
                         channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)]   #Estrella
                     pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
                     h = (hmin,hmax)
                     arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
                     dataOut.data_param = arrayParameters
                     return
             class SMPhaseCalibration(Operation):
                 __buffer = None
                 __initime = None
                 __dataReady = False
                 __isConfig = False
                 def __checkTime(self, currentTime, initTime, paramInterval, outputInterval):
                     dataTime = currentTime + paramInterval
                     deltaTime = dataTime - initTime
                     if deltaTime >= outputInterval or deltaTime < 0:
                         return True
                     return False
                 def __getGammas(self, pairs, d, phases):
                     gammas = numpy.zeros(2)
                     for i in range(len(pairs)):
                         pairi = pairs[i]
                         phip3 = phases[:,pairi[1]]
                         d3 = d[pairi[1]]
                         phip2 = phases[:,pairi[0]]
                         d2 = d[pairi[0]]
                         #Calculating gamma
             #             jdcos = alp1/(k*d1)
             #             jgamma = numpy.angle(numpy.exp(1j*(d0*alp1/d1 - alp0)))
                         jgamma = -phip2*d3/d2 - phip3
                         jgamma = numpy.angle(numpy.exp(1j*jgamma))
             #             jgamma[jgamma>numpy.pi] -= 2*numpy.pi
             #             jgamma[jgamma<-numpy.pi] += 2*numpy.pi
                         #Revised distribution
                         jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi))
                         #Histogram
                         nBins = 64.0
                         rmin = -0.5*numpy.pi
                         rmax = 0.5*numpy.pi
                         phaseHisto = numpy.histogram(jgammaArray, bins=nBins, range=(rmin,rmax))
                         meteorsY = phaseHisto[0]
                         phasesX = phaseHisto[1][:-1]
                         width = phasesX[1] - phasesX[0]
                         phasesX += width/2
                         #Gaussian aproximation
                         bpeak = meteorsY.argmax()
                         peak = meteorsY.max()
                         jmin = bpeak - 5
                         jmax = bpeak + 5 + 1
                         if jmin<0:
                             jmin = 0
                             jmax = 6
                         elif jmax > meteorsY.size:
                             jmin = meteorsY.size - 6
                             jmax = meteorsY.size
                         x0 = numpy.array([peak,bpeak,50])
                         coeff = optimize.leastsq(self.__residualFunction, x0, args=(meteorsY[jmin:jmax], phasesX[jmin:jmax]))
                         #Gammas
                         gammas[i] = coeff[0][1]
                     return gammas
                 def __residualFunction(self, coeffs, y, t):
                     return y - self.__gauss_function(t, coeffs)
                 def __gauss_function(self, t, coeffs):
                     return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2)
                 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
                     meteorOps = SMOperations()
                     nchan = 4
                     pairx = pairsList[0]
                     pairy = pairsList[1]
                     center_xangle = 0
                     center_yangle = 0
                     range_angle = numpy.array([10*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4])
                     ntimes = len(range_angle)
                     nstepsx = 20.0
                     nstepsy = 20.0
                     for iz in range(ntimes):
                         min_xangle = -range_angle[iz]/2 + center_xangle
                         max_xangle = range_angle[iz]/2 + center_xangle
                         min_yangle = -range_angle[iz]/2 + center_yangle
                         max_yangle = range_angle[iz]/2 + center_yangle
                         inc_x = (max_xangle-min_xangle)/nstepsx
                         inc_y = (max_yangle-min_yangle)/nstepsy
                         alpha_y = numpy.arange(nstepsy)*inc_y + min_yangle
                         alpha_x = numpy.arange(nstepsx)*inc_x + min_xangle
                         penalty = numpy.zeros((nstepsx,nstepsy))
                         jph_array = numpy.zeros((nchan,nstepsx,nstepsy))
                         jph = numpy.zeros(nchan)
                         # Iterations looking for the offset
                         for iy in range(int(nstepsy)):
                             for ix in range(int(nstepsx)):
                                 jph[pairy[1]] = alpha_y[iy]
                                 jph[pairy[0]] = -gammas[1] - alpha_y[iy]*d[pairy[1]]/d[pairy[0]]
                                 jph[pairx[1]] = alpha_x[ix]
                                 jph[pairx[0]] = -gammas[0] - alpha_x[ix]*d[pairx[1]]/d[pairx[0]]
                                 jph_array[:,ix,iy] = jph
                                 meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, d, jph)
                                 error = meteorsArray1[:,-1]
                                 ind1 = numpy.where(error==0)[0]
                                 penalty[ix,iy] = ind1.size
                         i,j = numpy.unravel_index(penalty.argmax(), penalty.shape)
                         phOffset = jph_array[:,i,j]
                         center_xangle = phOffset[pairx[1]]
                         center_yangle = phOffset[pairy[1]]
                     phOffset = numpy.angle(numpy.exp(1j*jph_array[:,i,j]))
                     phOffset = phOffset*180/numpy.pi
                     return phOffset
                 def run(self, dataOut, hmin, hmax, channelPositions=None, nHours = 1):
                     dataOut.flagNoData = True
                     self.__dataReady = False
                     dataOut.outputInterval = nHours*3600
                     if self.__isConfig == False:
             #                 self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
                         #Get Initial LTC time
                         self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
                         self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
                         self.__isConfig = True
                     if self.__buffer is None:
                         self.__buffer = dataOut.data_param.copy()
                     else:
                         self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
                     self.__dataReady = self.__checkTime(dataOut.utctime, self.__initime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
                     if self.__dataReady:
                         dataOut.utctimeInit = self.__initime
                         self.__initime += dataOut.outputInterval #to erase time offset
                         freq = dataOut.frequency
                         c = dataOut.C #m/s
                         lamb = c/freq
                         k = 2*numpy.pi/lamb
                         azimuth = 0
                         h = (hmin, hmax)
                         pairs = ((0,1),(2,3))
                         if channelPositions is None:
             #             channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)]   #T
                             channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)]   #Estrella
                         meteorOps = SMOperations()
                         pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
             #             distances1 = [-distances[0]*lamb, distances[1]*lamb, -distances[2]*lamb, distances[3]*lamb]
                         meteorsArray = self.__buffer
                         error = meteorsArray[:,-1]
                         boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14)
                         ind1 = numpy.where(boolError)[0]
                         meteorsArray = meteorsArray[ind1,:]
                         meteorsArray[:,-1] = 0
                         phases = meteorsArray[:,8:12]
                         #Calculate Gammas
                         gammas = self.__getGammas(pairs, distances, phases)
             #             gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180
                         #Calculate Phases
                         phasesOff = self.__getPhases(azimuth, h, pairs, distances, gammas, meteorsArray)
                         phasesOff = phasesOff.reshape((1,phasesOff.size))
                         dataOut.data_output = -phasesOff
                         dataOut.flagNoData = False
                         dataOut.channelList = pairslist0
                         self.__buffer = None
                     return
             class SMOperations():
                 def __init__(self):
                     return
                 def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, distances, jph):
                     arrayParameters = arrayParameters0.copy()
                     hmin = h[0]
                     hmax = h[1]
                     #Calculate AOA (Error N 3, 4)
                     #JONES ET AL. 1998
                     AOAthresh = numpy.pi/8
                     error = arrayParameters[:,-1]
                     phases = -arrayParameters[:,8:12] + jph
             #         phases = numpy.unwrap(phases)
                     arrayParameters[:,3:6], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, distances, error, AOAthresh, azimuth)
                     #Calculate Heights (Error N 13 and 14)
                     error = arrayParameters[:,-1]
                     Ranges = arrayParameters[:,1]
                     zenith = arrayParameters[:,4]
                     arrayParameters[:,2], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
                     #----------------------- Get Final data    ------------------------------------
             #         error = arrayParameters[:,-1]
             #         ind1 = numpy.where(error==0)[0]
             #         arrayParameters = arrayParameters[ind1,:]
                     return arrayParameters
                 def __getAOA(self, phases, pairsList, directions, error, AOAthresh, azimuth):
                     arrayAOA = numpy.zeros((phases.shape[0],3))
                     cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList,directions)
                     arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
                     cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
                     arrayAOA[:,2] = cosDirError
                     azimuthAngle = arrayAOA[:,0]
                     zenithAngle = arrayAOA[:,1]
                     #Setting Error
                     indError = numpy.where(numpy.logical_or(error == 3, error == 4))[0]
                     error[indError] = 0
                     #Number 3: AOA not fesible
                     indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
                     error[indInvalid] = 3
                     #Number 4: Large difference in AOAs obtained from different antenna baselines
                     indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
                     error[indInvalid] = 4
                     return arrayAOA, error
                 def __getDirectionCosines(self, arrayPhase, pairsList, distances):
                     #Initializing some variables
                     ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
                     ang_aux = ang_aux.reshape(1,ang_aux.size)
                     cosdir = numpy.zeros((arrayPhase.shape[0],2))
                     cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
                     for i in range(2):
                         ph0 = arrayPhase[:,pairsList[i][0]]
                         ph1 = arrayPhase[:,pairsList[i][1]]
                         d0 = distances[pairsList[i][0]]
                         d1 = distances[pairsList[i][1]]
                         ph0_aux = ph0 + ph1
                         ph0_aux = numpy.angle(numpy.exp(1j*ph0_aux))
             #             ph0_aux[ph0_aux > numpy.pi] -= 2*numpy.pi
             #             ph0_aux[ph0_aux < -numpy.pi] += 2*numpy.pi
                         #First Estimation
                         cosdir0[:,i] = (ph0_aux)/(2*numpy.pi*(d0 - d1))
                         #Most-Accurate Second Estimation
                         phi1_aux =  ph0 - ph1
                         phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
                         #Direction Cosine 1
                         cosdir1 = (phi1_aux + ang_aux)/(2*numpy.pi*(d0 + d1))
                         #Searching the correct Direction Cosine
                         cosdir0_aux = cosdir0[:,i]
                         cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
                         #Minimum Distance
                         cosDiff = (cosdir1 - cosdir0_aux)**2
                         indcos = cosDiff.argmin(axis = 1)
                         #Saving Value obtained
                         cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
                     return cosdir0, cosdir
                 def __calculateAOA(self, cosdir, azimuth):
                     cosdirX = cosdir[:,0]
                     cosdirY = cosdir[:,1]
                     zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
                     azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth#0 deg north, 90 deg east
                     angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
                     return angles
                 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
                     Ramb = 375  #Ramb = c/(2*PRF)
                     Re = 6371   #Earth Radius
                     heights = numpy.zeros(Ranges.shape)
                     R_aux = numpy.array([0,1,2])*Ramb
                     R_aux = R_aux.reshape(1,R_aux.size)
                     Ranges = Ranges.reshape(Ranges.size,1)
                     Ri = Ranges + R_aux
                     hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
                     #Check if there is a height between 70 and 110 km
                     h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
                     ind_h = numpy.where(h_bool == 1)[0]
                     hCorr = hi[ind_h, :]
                     ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
                     hCorr = hi[ind_hCorr]
                     heights[ind_h] = hCorr
                     #Setting Error
                     #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
                     #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
                     indError = numpy.where(numpy.logical_or(error == 13, error == 14))[0]
                     error[indError] = 0
                     indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
                     error[indInvalid2] = 14
                     indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
                     error[indInvalid1] = 13
                     return heights, error
                 def getPhasePairs(self, channelPositions):
                     chanPos = numpy.array(channelPositions)
                     listOper = list(itertools.combinations(range(5),2))
                     distances = numpy.zeros(4)
                     axisX = []
                     axisY = []
                     distX = numpy.zeros(3)
                     distY = numpy.zeros(3)
                     ix = 0
                     iy = 0
                     pairX = numpy.zeros((2,2))
                     pairY = numpy.zeros((2,2))
                     for i in range(len(listOper)):
                         pairi = listOper[i]
                         posDif = numpy.abs(chanPos[pairi[0],:] - chanPos[pairi[1],:])
                         if posDif[0] == 0:
                             axisY.append(pairi)
                             distY[iy] = posDif[1]
                             iy += 1
                         elif posDif[1] == 0:
                             axisX.append(pairi)
                             distX[ix] = posDif[0]
                             ix += 1
                     for i in range(2):
                         if i==0:
                             dist0 = distX
                             axis0 = axisX
                         else:
                             dist0 = distY
                             axis0 = axisY
                         side = numpy.argsort(dist0)[:-1]
                         axis0 = numpy.array(axis0)[side,:]
                         chanC = int(numpy.intersect1d(axis0[0,:], axis0[1,:])[0])
                         axis1 = numpy.unique(numpy.reshape(axis0,4))
                         side = axis1[axis1 != chanC]
                         diff1 = chanPos[chanC,i] - chanPos[side[0],i]
                         diff2 = chanPos[chanC,i] - chanPos[side[1],i]
                         if diff1<0:
                             chan2 = side[0]
                             d2 = numpy.abs(diff1)
                             chan1 = side[1]
                             d1 = numpy.abs(diff2)
                         else:
                             chan2 = side[1]
                             d2 = numpy.abs(diff2)
                             chan1 = side[0]
                             d1 = numpy.abs(diff1)
                         if i==0:
                             chanCX = chanC
                             chan1X = chan1
                             chan2X = chan2
                             distances[0:2] = numpy.array([d1,d2])
                         else:
                             chanCY = chanC
                             chan1Y = chan1
                             chan2Y = chan2
                             distances[2:4] = numpy.array([d1,d2])
             #         axisXsides = numpy.reshape(axisX[ix,:],4)
             #
             #         channelCentX = int(numpy.intersect1d(pairX[0,:], pairX[1,:])[0])
             #         channelCentY = int(numpy.intersect1d(pairY[0,:], pairY[1,:])[0])
             #
             #         ind25X = numpy.where(pairX[0,:] != channelCentX)[0][0]
             #         ind20X = numpy.where(pairX[1,:] != channelCentX)[0][0]
             #         channel25X = int(pairX[0,ind25X])
             #         channel20X = int(pairX[1,ind20X])
             #         ind25Y = numpy.where(pairY[0,:] != channelCentY)[0][0]
             #         ind20Y = numpy.where(pairY[1,:] != channelCentY)[0][0]
             #         channel25Y = int(pairY[0,ind25Y])
             #         channel20Y = int(pairY[1,ind20Y])
             #         pairslist = [(channelCentX, channel25X),(channelCentX, channel20X),(channelCentY,channel25Y),(channelCentY, channel20Y)]
                     pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]
                     return pairslist, distances
             #     def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
             #
             #         arrayAOA = numpy.zeros((phases.shape[0],3))
             #         cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
             #
             #         arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
             #         cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
             #         arrayAOA[:,2] = cosDirError
             #
             #         azimuthAngle = arrayAOA[:,0]
             #         zenithAngle = arrayAOA[:,1]
             #
             #         #Setting Error
             #         #Number 3: AOA not fesible
             #         indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
             #         error[indInvalid] = 3
             #         #Number 4: Large difference in AOAs obtained from different antenna baselines
             #         indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
             #         error[indInvalid] = 4
             #         return arrayAOA, error
             #
             #     def __getDirectionCosines(self, arrayPhase, pairsList):
             #
             #         #Initializing some variables
             #         ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
             #         ang_aux = ang_aux.reshape(1,ang_aux.size)
             #
             #         cosdir = numpy.zeros((arrayPhase.shape[0],2))
             #         cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
             #
             #
             #         for i in range(2):
             #             #First Estimation
             #             phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
             #             #Dealias
             #             indcsi = numpy.where(phi0_aux > numpy.pi)
             #             phi0_aux[indcsi] -= 2*numpy.pi
             #             indcsi = numpy.where(phi0_aux < -numpy.pi)
             #             phi0_aux[indcsi] += 2*numpy.pi
             #             #Direction Cosine 0
             #             cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
             #
             #             #Most-Accurate Second Estimation
             #             phi1_aux =  arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
             #             phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
             #             #Direction Cosine 1
             #             cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
             #
             #             #Searching the correct Direction Cosine
             #             cosdir0_aux = cosdir0[:,i]
             #             cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
             #             #Minimum Distance
             #             cosDiff = (cosdir1 - cosdir0_aux)**2
             #             indcos = cosDiff.argmin(axis = 1)
             #             #Saving Value obtained
             #             cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
             #
             #         return cosdir0, cosdir
             #
             #     def __calculateAOA(self, cosdir, azimuth):
             #         cosdirX = cosdir[:,0]
             #         cosdirY = cosdir[:,1]
             #
             #         zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
             #         azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
             #         angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
             #
             #         return angles
             #
             #     def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
             #
             #         Ramb = 375  #Ramb = c/(2*PRF)
             #         Re = 6371   #Earth Radius
             #         heights = numpy.zeros(Ranges.shape)
             #
             #         R_aux = numpy.array([0,1,2])*Ramb
             #         R_aux = R_aux.reshape(1,R_aux.size)
             #
             #         Ranges = Ranges.reshape(Ranges.size,1)
             #
             #         Ri = Ranges + R_aux
             #         hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
             #
             #         #Check if there is a height between 70 and 110 km
             #         h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
             #         ind_h = numpy.where(h_bool == 1)[0]
             #
             #         hCorr = hi[ind_h, :]
             #         ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
             #
             #         hCorr = hi[ind_hCorr]
             #         heights[ind_h] = hCorr
             #
             #         #Setting Error
             #         #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
             #         #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
             #
             #         indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
             #         error[indInvalid2] = 14
             #         indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
             #         error[indInvalid1] = 13
             #
             #         return heights, error

schainpy/model/utils/jroutils_publish.py +15 -3

             '''
             @author: Juan C. Espinoza
             '''
             import time
             import json
             import numpy
             import paho.mqtt.client as mqtt
             import zmq
             import cPickle as pickle
             import datetime
             from zmq.utils.monitor import recv_monitor_message
             from functools import wraps
             from threading import Thread
             from multiprocessing import Process
             from schainpy.model.proc.jroproc_base import Operation, ProcessingUnit
             MAXNUMX = 100
             MAXNUMY = 100
             class PrettyFloat(float):
                 def __repr__(self):
                     return '%.2f' % self
             def roundFloats(obj):
                 if isinstance(obj, list):
                     return map(roundFloats, obj)
                 elif isinstance(obj, float):
                     return round(obj, 2)
             def decimate(z):
                 # dx = int(len(self.x)/self.__MAXNUMX) + 1
                 dy = int(len(z[0])/MAXNUMY) + 1
                 return z[::, ::dy]
             class throttle(object):
                 """Decorator that prevents a function from being called more than once every
                 time period.
                 To create a function that cannot be called more than once a minute, but
                 will sleep until it can be called:
                 @throttle(minutes=1)
                 def foo():
                   pass
                 for i in range(10):
                   foo()
                   print "This function has run %s times." % i
                 """
                 def __init__(self, seconds=0, minutes=0, hours=0):
                     self.throttle_period = datetime.timedelta(
                         seconds=seconds, minutes=minutes, hours=hours
                     )
                     self.time_of_last_call = datetime.datetime.min
                 def __call__(self, fn):
                     @wraps(fn)
                     def wrapper(*args, **kwargs):
                         now = datetime.datetime.now()
                         time_since_last_call = now - self.time_of_last_call
                         time_left = self.throttle_period - time_since_last_call
                         if time_left > datetime.timedelta(seconds=0):
                             return
                         self.time_of_last_call = datetime.datetime.now()
                         return fn(*args, **kwargs)
                     return wrapper
             class PublishData(Operation):
                 """Clase publish."""
                 def __init__(self, **kwargs):
                     """Inicio."""
                     Operation.__init__(self, **kwargs)
                     self.isConfig = False
                     self.client = None
                     self.zeromq = None
                     self.mqtt = None
                 def on_disconnect(self, client, userdata, rc):
                     if rc != 0:
                         print("Unexpected disconnection.")
                     self.connect()
                 def connect(self):
                     print 'trying to connect'
                     try:
                         self.client.connect(
                             host=self.host,
                             port=self.port,
                             keepalive=60*10,
                             bind_address='')
                         self.client.loop_start()
                         # self.client.publish(
                         #     self.topic + 'SETUP',
                         #     json.dumps(setup),
                         #     retain=True
                         #     )
                     except:
                         print "MQTT Conection error."
                         self.client = False
                 def setup(self, port=1883, username=None, password=None, clientId="user", zeromq=1, **kwargs):
                     self.counter = 0
                     self.topic = kwargs.get('topic', 'schain')
                     self.delay = kwargs.get('delay', 0)
                     self.plottype = kwargs.get('plottype', 'spectra')
                     self.host = kwargs.get('host', "10.10.10.82")
                     self.port = kwargs.get('port', 3000)
                     self.clientId = clientId
                     self.cnt = 0
                     self.zeromq = zeromq
                     self.mqtt = kwargs.get('plottype', 0)
                     self.client = None
                     setup = []
                     if mqtt is 1:
                         self.client = mqtt.Client(
                             client_id=self.clientId + self.topic + 'SCHAIN',
                             clean_session=True)
                         self.client.on_disconnect = self.on_disconnect
                         self.connect()
                         for plot in self.plottype:
                             setup.append({
                                 'plot': plot,
                                 'topic': self.topic + plot,
                                 'title': getattr(self, plot + '_' + 'title', False),
                                 'xlabel': getattr(self, plot + '_' + 'xlabel', False),
                                 'ylabel': getattr(self, plot + '_' + 'ylabel', False),
                                 'xrange': getattr(self, plot + '_' + 'xrange', False),
                                 'yrange': getattr(self, plot + '_' + 'yrange', False),
                                 'zrange': getattr(self, plot + '_' + 'zrange', False),
                             })
                     if zeromq is 1:
                         context = zmq.Context()
                         self.zmq_socket = context.socket(zmq.PUSH)
                         server = kwargs.get('server', 'zmq.pipe')
                         if 'tcp://' in server:
                             address = server
                         else:
                             address = 'ipc:///tmp/%s' % server
                         self.zmq_socket.connect(address)
                         time.sleep(1)
                 def publish_data(self):
                     self.dataOut.finished = False
                     if self.mqtt is 1:
                         yData = self.dataOut.heightList[:2].tolist()
                         if self.plottype == 'spectra':
                             data = getattr(self.dataOut, 'data_spc')
                             z = data/self.dataOut.normFactor
                             zdB = 10*numpy.log10(z)
                             xlen, ylen = zdB[0].shape
                             dx = int(xlen/MAXNUMX) + 1
                             dy = int(ylen/MAXNUMY) + 1
                             Z = [0 for i in self.dataOut.channelList]
                             for i in self.dataOut.channelList:
                                 Z[i] = zdB[i][::dx, ::dy].tolist()
                             payload = {
                                 'timestamp': self.dataOut.utctime,
                                 'data': roundFloats(Z),
                                 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
                                 'interval': self.dataOut.getTimeInterval(),
                                 'type': self.plottype,
                                 'yData': yData
                             }
                             # print payload
                         elif self.plottype in ('rti', 'power'):
                             data = getattr(self.dataOut, 'data_spc')
                             z = data/self.dataOut.normFactor
                             avg = numpy.average(z, axis=1)
                             avgdB = 10*numpy.log10(avg)
                             xlen, ylen = z[0].shape
                             dy = numpy.floor(ylen/self.__MAXNUMY) + 1
                             AVG = [0 for i in self.dataOut.channelList]
                             for i in self.dataOut.channelList:
                                 AVG[i] = avgdB[i][::dy].tolist()
                             payload = {
                                 'timestamp': self.dataOut.utctime,
                                 'data': roundFloats(AVG),
                                 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
                                 'interval': self.dataOut.getTimeInterval(),
                                 'type': self.plottype,
                                 'yData': yData
                             }
                         elif self.plottype == 'noise':
                             noise = self.dataOut.getNoise()/self.dataOut.normFactor
                             noisedB = 10*numpy.log10(noise)
                             payload = {
                                 'timestamp': self.dataOut.utctime,
                                 'data': roundFloats(noisedB.reshape(-1, 1).tolist()),
                                 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
                                 'interval': self.dataOut.getTimeInterval(),
                                 'type': self.plottype,
                                 'yData': yData
                             }
                         elif self.plottype == 'snr':
                             data = getattr(self.dataOut, 'data_SNR')
                             avgdB = 10*numpy.log10(data)
                             ylen = data[0].size
                             dy = numpy.floor(ylen/self.__MAXNUMY) + 1
                             AVG = [0 for i in self.dataOut.channelList]
                             for i in self.dataOut.channelList:
                                 AVG[i] = avgdB[i][::dy].tolist()
                             payload = {
                                 'timestamp': self.dataOut.utctime,
                                 'data': roundFloats(AVG),
                                 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
                                 'type': self.plottype,
                                 'yData': yData
                             }
                         else:
                             print "Tipo de grafico invalido"
                             payload = {
                                 'data': 'None',
                                 'timestamp': 'None',
                                 'type': None
                             }
                             # print 'Publishing data to {}'.format(self.host)
                         self.client.publish(self.topic + self.plottype, json.dumps(payload), qos=0)
                     if self.zeromq is 1:
                         print '[Sending] {} - {}'.format(self.dataOut.type, self.dataOut.datatime)
                         self.zmq_socket.send_pyobj(self.dataOut)
                 def run(self, dataOut, **kwargs):
                     self.dataOut = dataOut
                     if not self.isConfig:
                         self.setup(**kwargs)
                         self.isConfig = True
                     self.publish_data()
                     time.sleep(self.delay)
                 def close(self):
                     if self.zeromq is 1:
                         self.dataOut.finished = True
                         self.zmq_socket.send_pyobj(self.dataOut)
                     if self.client:
                         self.client.loop_stop()
                         self.client.disconnect()
             class ReceiverData(ProcessingUnit, Process):
                 throttle_value = 5
                 def __init__(self, **kwargs):
                     ProcessingUnit.__init__(self, **kwargs)
                     Process.__init__(self)
                     self.mp = False
                     self.isConfig = False
                     self.isWebConfig = False
                     self.plottypes =[]
                     self.connections = 0
                     server = kwargs.get('server', 'zmq.pipe')
                     plot_server = kwargs.get('plot_server', 'zmq.web')
                     if 'tcp://' in server:
                         address = server
                     else:
                         address = 'ipc:///tmp/%s' % server
                     if 'tcp://' in plot_server:
                         plot_address = plot_server
                     else:
                         plot_address = 'ipc:///tmp/%s' % plot_server
                     self.address = address
                     self.plot_address = plot_address
                     self.plottypes = [s.strip() for s in kwargs.get('plottypes', 'rti').split(',')]
                     self.realtime = kwargs.get('realtime', False)
-                    self.throttle_value = kwargs.get('throttle', 10)
+                    self.throttle_value = kwargs.get('throttle', 5)
                     self.sendData = self.initThrottle(self.throttle_value)
                     self.setup()
                 def setup(self):
                     self.data = {}
                     self.data['times'] = []
                     for plottype in self.plottypes:
                         self.data[plottype] = {}
                     self.data['noise'] = {}
                     self.data['throttle'] = self.throttle_value
                     self.data['ENDED'] = False
                     self.isConfig = True
                     self.data_web = {}
                 def event_monitor(self, monitor):
                     events = {}
                     for name in dir(zmq):
                         if name.startswith('EVENT_'):
                             value = getattr(zmq, name)
                             events[value] = name
                     while monitor.poll():
                         evt = recv_monitor_message(monitor)
                         if evt['event'] == 32:
                             self.connections += 1
                         if evt['event'] == 512:
                             pass
                         if self.connections == 0 and self.started is True:
                             self.ended = True
                             # send('ENDED')
                         evt.update({'description': events[evt['event']]})
                         if evt['event'] == zmq.EVENT_MONITOR_STOPPED:
                             break
                     monitor.close()
                     print("event monitor thread done!")
                 def initThrottle(self, throttle_value):
                     @throttle(seconds=throttle_value)
                     def sendDataThrottled(fn_sender, data):
                         fn_sender(data)
                     return sendDataThrottled
                 def send(self, data):
                     # print '[sending] data=%s size=%s' % (data.keys(), len(data['times']))
                     self.sender.send_pyobj(data)
                 def update(self):
                     t = self.dataOut.utctime
                     self.data['times'].append(t)
                     self.data['dataOut'] = self.dataOut
                     for plottype in self.plottypes:
                         if plottype == 'spc':
                             z = self.dataOut.data_spc/self.dataOut.normFactor
                             self.data[plottype] = 10*numpy.log10(z)
                             self.data['noise'][t] = 10*numpy.log10(self.dataOut.getNoise()/self.dataOut.normFactor)
+                        if plottype == 'cspc':
+                            jcoherence = self.dataOut.data_cspc/numpy.sqrt(self.dataOut.data_spc*self.dataOut.data_spc)
+                            self.data['cspc_coh'] = numpy.abs(jcoherence)
+                            self.data['cspc_phase'] = numpy.arctan2(jcoherence.imag, jcoherence.real)*180/numpy.pi
                         if plottype == 'rti':
                             self.data[plottype][t] = self.dataOut.getPower()
                         if plottype == 'snr':
                             self.data[plottype][t] = 10*numpy.log10(self.dataOut.data_SNR)
                         if plottype == 'dop':
                             self.data[plottype][t] = 10*numpy.log10(self.dataOut.data_DOP)
+                        if plottype == 'mean':
+                            self.data[plottype][t] = self.dataOut.data_MEAN
+                        if plottype == 'std':
+                            self.data[plottype][t] = self.dataOut.data_STD
                         if plottype == 'coh':
                             self.data[plottype][t] = self.dataOut.getCoherence()
                         if plottype == 'phase':
                             self.data[plottype][t] = self.dataOut.getCoherence(phase=True)
                         if self.realtime:
-                            self.data_web[plottype] = roundFloats(decimate(self.data[plottype][t]).tolist())
+                            self.data_web['timestamp'] = t
-                            self.data_web['timestamp'] = t
                             if plottype == 'spc':
                                 self.data_web[plottype] = roundFloats(decimate(self.data[plottype]).tolist())
+                            elif plottype == 'cspc':
+                                self.data_web['cspc_coh'] = roundFloats(decimate(self.data['cspc_coh']).tolist())
+                                self.data_web['cspc_phase'] = roundFloats(decimate(self.data['cspc_phase']).tolist())
+                            elif plottype == 'noise':
+                                self.data_web['noise'] = roundFloats(self.data['noise'][t].tolist())
                             else:
                                 self.data_web[plottype] = roundFloats(decimate(self.data[plottype][t]).tolist())
                             self.data_web['interval'] = self.dataOut.getTimeInterval()
                             self.data_web['type'] = plottype
                 def run(self):
                     print '[Starting] {} from {}'.format(self.name, self.address)
                     self.context = zmq.Context()
                     self.receiver = self.context.socket(zmq.PULL)
                     self.receiver.bind(self.address)
                     monitor = self.receiver.get_monitor_socket()
                     self.sender = self.context.socket(zmq.PUB)
                     if self.realtime:
                         self.sender_web = self.context.socket(zmq.PUB)
                         self.sender_web.connect(self.plot_address)
                         time.sleep(1)
                     self.sender.bind("ipc:///tmp/zmq.plots")
                     t = Thread(target=self.event_monitor, args=(monitor,))
                     t.start()
                     while True:
                         self.dataOut = self.receiver.recv_pyobj()
                         # print '[Receiving] {} - {}'.format(self.dataOut.type,
                         #                                    self.dataOut.datatime.ctime())
                         self.update()
                         if self.dataOut.finished is True:
                             self.send(self.data)
                             self.connections -= 1
                             if self.connections == 0 and self.started:
                                 self.ended = True
                                 self.data['ENDED'] = True
                                 self.send(self.data)
                                 self.setup()
                         else:
                             if self.realtime:
                                 self.send(self.data)
                                 self.sender_web.send_string(json.dumps(self.data_web))
                             else:
                                 self.sendData(self.send, self.data)
                             self.started = True
                     return
                 def sendToWeb(self):
                     if not self.isWebConfig:
                         context = zmq.Context()
                         sender_web_config = context.socket(zmq.PUB)
                         if 'tcp://' in self.plot_address:
                             dum, address, port = self.plot_address.split(':')
                             conf_address = '{}:{}:{}'.format(dum, address, int(port)+1)
                         else:
                             conf_address = self.plot_address + '.config'
                         sender_web_config.bind(conf_address)
                         time.sleep(1)
                         for kwargs in self.operationKwargs.values():
                             if 'plot' in kwargs:
                                 print '[Sending] Config data to web for {}'.format(kwargs['code'].upper())
                                 sender_web_config.send_string(json.dumps(kwargs))
                         self.isWebConfig = True

General Comments 0

Write
Preview

You need to be logged in to leave comments. Login now

No TODOs yet

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages