import sys import numpy,math from scipy import interpolate from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon from schainpy.utils import log from time import time, mktime, strptime, gmtime, ctime import os class VoltageProc(ProcessingUnit): def __init__(self): ProcessingUnit.__init__(self) self.dataOut = Voltage() self.flip = 1 self.setupReq = False #self.dataOut.test=1 def run(self): #import time #time.sleep(3) if self.dataIn.type == 'AMISR': self.__updateObjFromAmisrInput() if self.dataIn.type == 'Voltage': self.dataOut.copy(self.dataIn) #self.dataOut.flagNoData=True #print(self.dataOut.data[-1,:]) #print(ctime(self.dataOut.utctime)) #print(self.dataOut.heightList) #print(self.dataOut.nHeights) #exit(1) #print(self.dataOut.data[6,:32]) #print(self.dataOut.data[0,320-5:320+5-5]) ##print(self.dataOut.heightList[-20:]) #print(numpy.shape(self.dataOut.data)) #print(self.dataOut.code) #print(numpy.shape(self.dataOut.code)) #exit(1) #print(self.dataOut.CurrentBlock) #print(self.dataOut.data[0,:,0]) #print(numpy.shape(self.dataOut.data)) #print(self.dataOut.data[0,:,1666:1666+320]) #exit(1) #print(self.dataOut.utctime) #self.dataOut.test+=1 def __updateObjFromAmisrInput(self): 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.flagNoData = self.dataIn.flagNoData self.dataOut.data = self.dataIn.data self.dataOut.utctime = self.dataIn.utctime self.dataOut.channelList = self.dataIn.channelList #self.dataOut.timeInterval = self.dataIn.timeInterval self.dataOut.heightList = self.dataIn.heightList self.dataOut.nProfiles = self.dataIn.nProfiles self.dataOut.nCohInt = self.dataIn.nCohInt self.dataOut.ippSeconds = self.dataIn.ippSeconds self.dataOut.frequency = self.dataIn.frequency self.dataOut.azimuth = self.dataIn.azimuth self.dataOut.zenith = self.dataIn.zenith self.dataOut.beam.codeList = self.dataIn.beam.codeList self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList self.dataOut.beam.zenithList = self.dataIn.beam.zenithList class selectChannels(Operation): def run(self, dataOut, channelList): channelIndexList = [] self.dataOut = dataOut for channel in channelList: if channel not in self.dataOut.channelList: raise ValueError("Channel %d is not in %s" %(channel, str(self.dataOut.channelList))) index = self.dataOut.channelList.index(channel) channelIndexList.append(index) self.selectChannelsByIndex(channelIndexList) return self.dataOut def selectChannelsByIndex(self, channelIndexList): """ Selecciona un bloque de datos en base a canales segun el channelIndexList Input: channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7] Affected: self.dataOut.data self.dataOut.channelIndexList self.dataOut.nChannels self.dataOut.m_ProcessingHeader.totalSpectra self.dataOut.systemHeaderObj.numChannels self.dataOut.m_ProcessingHeader.blockSize Return: None """ for channelIndex in channelIndexList: if channelIndex not in self.dataOut.channelIndexList: raise ValueError("The value %d in channelIndexList is not valid" %channelIndex) if self.dataOut.type == 'Voltage': if self.dataOut.flagDataAsBlock: """ Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis] """ data = self.dataOut.data[channelIndexList,:,:] else: data = self.dataOut.data[channelIndexList,:] self.dataOut.data = data # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList] self.dataOut.channelList = range(len(channelIndexList)) elif self.dataOut.type == 'Spectra': data_spc = self.dataOut.data_spc[channelIndexList, :] data_dc = self.dataOut.data_dc[channelIndexList, :] self.dataOut.data_spc = data_spc self.dataOut.data_dc = data_dc # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList] self.dataOut.channelList = range(len(channelIndexList)) self.__selectPairsByChannel(channelIndexList) return 1 def __selectPairsByChannel(self, channelList=None): if channelList == None: return pairsIndexListSelected = [] for pairIndex in self.dataOut.pairsIndexList: # First pair if self.dataOut.pairsList[pairIndex][0] not in channelList: continue # Second pair if self.dataOut.pairsList[pairIndex][1] not in channelList: continue pairsIndexListSelected.append(pairIndex) if not pairsIndexListSelected: self.dataOut.data_cspc = None self.dataOut.pairsList = [] return self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected] self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected] return class selectHeights(Operation): def run(self, dataOut, minHei=None, maxHei=None): """ Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango minHei <= height <= maxHei Input: minHei : valor minimo de altura a considerar maxHei : valor maximo de altura a considerar Affected: Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex Return: 1 si el metodo se ejecuto con exito caso contrario devuelve 0 """ self.dataOut = dataOut if minHei == None: minHei = self.dataOut.heightList[0] if maxHei == None: maxHei = self.dataOut.heightList[-1] if (minHei < self.dataOut.heightList[0]): minHei = self.dataOut.heightList[0] if (maxHei > self.dataOut.heightList[-1]): maxHei = self.dataOut.heightList[-1] minIndex = 0 maxIndex = 0 heights = self.dataOut.heightList inda = numpy.where(heights >= minHei) indb = numpy.where(heights <= maxHei) try: minIndex = inda[0][0] except: minIndex = 0 try: maxIndex = indb[0][-1] except: maxIndex = len(heights) self.selectHeightsByIndex(minIndex, maxIndex) #print(self.dataOut.nHeights) return self.dataOut def selectHeightsByIndex(self, minIndex, maxIndex): """ Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango minIndex <= index <= maxIndex Input: minIndex : valor de indice minimo de altura a considerar maxIndex : valor de indice maximo de altura a considerar Affected: self.dataOut.data self.dataOut.heightList Return: 1 si el metodo se ejecuto con exito caso contrario devuelve 0 """ if self.dataOut.type == 'Voltage': if (minIndex < 0) or (minIndex > maxIndex): raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex)) if (maxIndex >= self.dataOut.nHeights): maxIndex = self.dataOut.nHeights #voltage if self.dataOut.flagDataAsBlock: """ Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis] """ data = self.dataOut.data[:,:, minIndex:maxIndex] else: data = self.dataOut.data[:, minIndex:maxIndex] # firstHeight = self.dataOut.heightList[minIndex] self.dataOut.data = data self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex] if self.dataOut.nHeights <= 1: raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights)) elif self.dataOut.type == 'Spectra': if (minIndex < 0) or (minIndex > maxIndex): raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % ( minIndex, maxIndex)) if (maxIndex >= self.dataOut.nHeights): maxIndex = self.dataOut.nHeights - 1 # Spectra data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1] data_cspc = None if self.dataOut.data_cspc is not None: data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1] data_dc = None if self.dataOut.data_dc is not None: data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1] self.dataOut.data_spc = data_spc self.dataOut.data_cspc = data_cspc self.dataOut.data_dc = data_dc self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1] return 1 class filterByHeights(Operation): def run(self, dataOut, window): deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] if window == None: window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight newdelta = deltaHeight * window r = dataOut.nHeights % window newheights = (dataOut.nHeights-r)/window if newheights <= 1: raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window)) if dataOut.flagDataAsBlock: """ Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis] """ buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)] buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window) buffer = numpy.sum(buffer,3) else: buffer = dataOut.data[:,0:int(dataOut.nHeights-r)] buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window)) buffer = numpy.sum(buffer,2) dataOut.data = buffer dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta dataOut.windowOfFilter = window return dataOut class setH0(Operation): def run(self, dataOut, h0, deltaHeight = None): if not deltaHeight: deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] nHeights = dataOut.nHeights newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight dataOut.heightList = newHeiRange return dataOut class deFlip(Operation): def __init__(self): self.flip = 1 def run(self, dataOut, channelList = []): data = dataOut.data.copy() #print(dataOut.channelList) #exit() if channelList==1: #PARCHE channelList=[1] dataOut.FlipChannels=channelList if dataOut.flagDataAsBlock: flip = self.flip profileList = list(range(dataOut.nProfiles)) if not channelList: for thisProfile in profileList: data[:,thisProfile,:] = data[:,thisProfile,:]*flip flip *= -1.0 else: for thisChannel in channelList: if thisChannel not in dataOut.channelList: continue for thisProfile in profileList: data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip flip *= -1.0 self.flip = flip else: if not channelList: data[:,:] = data[:,:]*self.flip else: channelList=[1] #print(self.flip) for thisChannel in channelList: if thisChannel not in dataOut.channelList: continue data[thisChannel,:] = data[thisChannel,:]*self.flip self.flip *= -1. dataOut.data = data return dataOut class setAttribute(Operation): ''' Set an arbitrary attribute(s) to dataOut ''' def __init__(self): Operation.__init__(self) self._ready = False def run(self, dataOut, **kwargs): for key, value in kwargs.items(): setattr(dataOut, key, value) return dataOut @MPDecorator class printAttribute(Operation): ''' Print an arbitrary attribute of dataOut ''' def __init__(self): Operation.__init__(self) def run(self, dataOut, attributes): if isinstance(attributes, str): attributes = [attributes] for attr in attributes: if hasattr(dataOut, attr): log.log(getattr(dataOut, attr), attr) class interpolateHeights(Operation): def run(self, dataOut, topLim, botLim): #69 al 72 para julia #82-84 para meteoros if len(numpy.shape(dataOut.data))==2: sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1))) #dataOut.data[:,botLim:limSup+1] = sampInterp dataOut.data[:,botLim:topLim+1] = sampInterp else: nHeights = dataOut.data.shape[2] x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights))) y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))] f = interpolate.interp1d(x, y, axis = 2) xnew = numpy.arange(botLim,topLim+1) ynew = f(xnew) dataOut.data[:,:,botLim:topLim+1] = ynew return dataOut class LagsReshape(Operation): """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction. Parameters: ----------- Example -------- op = proc_unit.addOperation(name='LagsReshape') """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.buffer=None self.buffer_HR=None self.buffer_HRonelag=None def LagDistribution(self,dataOut): dataOut.datapure=numpy.copy(dataOut.data[:,0:dataOut.NSCAN,:]) self.buffer = numpy.zeros((dataOut.nChannels, int(dataOut.NSCAN/dataOut.DPL), dataOut.nHeights,dataOut.DPL), dtype='complex') for j in range(int(self.buffer.shape[1]/2)): for i in range(dataOut.DPL): if j+1==int(self.buffer.shape[1]/2) and i+1==dataOut.DPL: self.buffer[:,2*j:,:,i]=dataOut.datapure[:,2*i+int(2*j*dataOut.DPL):,:] else: self.buffer[:,2*j:2*(j+1),:,i]=dataOut.datapure[:,2*i+int(2*j*dataOut.DPL):2*(i+1)+int(2*j*dataOut.DPL),:] return self.buffer def HeightReconstruction(self,dataOut): self.buffer_HR = numpy.zeros((int(dataOut.NSCAN/dataOut.DPL), dataOut.nHeights,dataOut.DPL), dtype='complex') #self.buffer_HR[0,:,:,:]=dataOut.datalags[0,:,:,:] #No Lags for i in range(int(dataOut.DPL)): #Only channel B if i==0: self.buffer_HR[:,:,i]=dataOut.datalags[1,:,:,i] else: self.buffer_HR[:,:,i]=self.HRonelag(dataOut,i) return self.buffer_HR def HRonelag(self,dataOut,whichlag): self.buffer_HRonelag = numpy.zeros((int(dataOut.NSCAN/dataOut.DPL), dataOut.nHeights), dtype='complex') for i in range(self.buffer_HRonelag.shape[0]): for j in range(dataOut.nHeights): if j+int(2*whichlag)=nkill/2 and k 2): nums_min= int(ndata/divider) else: nums_min=2 sump=0.0 sumq=0.0 j=0 cont=1 while ( (cont==1) and (j nums_min): rtest= float(j/(j-1)) +1.0/ndata if( (sumq*j) > (rtest*sump*sump ) ): j=j-1 sump-= data[j] sumq-=data[j]*data[j] cont= 0 noise= (sump/j) return noise def run(self, dataOut, NLAG=16, NRANGE=0, NCAL=0, DPL=11, NDN=0, NDT=66, NDP=66, NSCAN=132, flags_array=(0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300), NAVG=16, nkill=6, **kwargs): dataOut.NLAG=NLAG dataOut.NR=len(dataOut.channelList) dataOut.NRANGE=NRANGE dataOut.NCAL=NCAL dataOut.DPL=DPL dataOut.NDN=NDN dataOut.NDT=NDT dataOut.NDP=NDP dataOut.NSCAN=NSCAN dataOut.DH=dataOut.heightList[1]-dataOut.heightList[0] dataOut.H0=int(dataOut.heightList[0]) dataOut.flags_array=flags_array dataOut.NAVG=NAVG dataOut.nkill=nkill dataOut.flagNoData = True self.get_dc(dataOut) self.get_products_cabxys(dataOut) self.cabxys_navg(dataOut) self.noise_estimation4x_DP(dataOut) self.kabxys(dataOut) return dataOut class IntegrationDP(Operation): """Operation to integrate the Double Pulse data. Parameters: ----------- nint : int Number of integrations. Example -------- op = proc_unit.addOperation(name='IntegrationDP', optype='other') op.addParameter(name='nint', value='30', format='int') """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.counter=0 self.aux=0 self.init_time=None def integration_for_double_pulse(self,dataOut): #print("inside") #print(self.aux) if self.aux==1: #print("CurrentBlockBBBBB: ",dataOut.CurrentBlock) #print(dataOut.datatime) #dataOut.TimeBlockDate_for_dp_power=dataOut.TimeBlockDate ########dataOut.TimeBlockSeconds_for_dp_power=dataOut.LastAVGDate #print("Date: ",dataOut.TimeBlockDate_for_dp_power) #dataOut.TimeBlockSeconds_for_dp_power=mktime(strptime(dataOut.TimeBlockDate_for_dp_power)) dataOut.TimeBlockSeconds_for_dp_power=dataOut.utctime#dataOut.TimeBlockSeconds-18000 #dataOut.TimeBlockSeconds_for_dp_power=dataOut.LastAVGDate #print("Seconds: ",dataOut.TimeBlockSeconds_for_dp_power) dataOut.bd_time=gmtime(dataOut.TimeBlockSeconds_for_dp_power) #print(dataOut.bd_time) #exit() dataOut.year=dataOut.bd_time.tm_year+(dataOut.bd_time.tm_yday-1)/364.0 dataOut.ut_Faraday=dataOut.bd_time.tm_hour+dataOut.bd_time.tm_min/60.0+dataOut.bd_time.tm_sec/3600.0 #print("date: ", dataOut.TimeBlockDate) self.aux=0 #print("after") if self.counter==0: tmpx=numpy.zeros((dataOut.NDP,dataOut.DPL,2),'float32') dataOut.kabxys_integrated=[tmpx,tmpx,tmpx,tmpx,tmpx,tmpx,tmpx,tmpx,tmpx,tmpx,tmpx,tmpx,tmpx,tmpx] self.init_time=dataOut.utctime if self.counter < dataOut.nint: #print("HERE") dataOut.final_cross_products=[dataOut.kax,dataOut.kay,dataOut.kbx,dataOut.kby,dataOut.kax2,dataOut.kay2,dataOut.kbx2,dataOut.kby2,dataOut.kaxbx,dataOut.kaxby,dataOut.kaybx,dataOut.kayby,dataOut.kaxay,dataOut.kbxby] for ind in range(len(dataOut.kabxys_integrated)): #final cross products dataOut.kabxys_integrated[ind]=dataOut.kabxys_integrated[ind]+dataOut.final_cross_products[ind] self.counter+=1 if self.counter==dataOut.nint-1: self.aux=1 #dataOut.TimeBlockDate_for_dp_power=dataOut.TimeBlockDate if self.counter==dataOut.nint: dataOut.flagNoData=False dataOut.utctime=self.init_time self.counter=0 def run(self,dataOut,nint=20): dataOut.flagNoData=True dataOut.nint=nint dataOut.paramInterval=0#int(dataOut.nint*dataOut.header[7][0]*2 ) dataOut.lat=-11.95 dataOut.lon=-76.87 self.integration_for_double_pulse(dataOut) return dataOut class SumFlips(Operation): """Operation to sum the flip and unflip part of certain cross products of the Double Pulse. Parameters: ----------- None Example -------- op = proc_unit.addOperation(name='SumFlips', optype='other') """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) def rint2DP(self,dataOut): dataOut.rnint2=numpy.zeros(dataOut.DPL,'float32') for l in range(dataOut.DPL): dataOut.rnint2[l]=1.0/(dataOut.nint*dataOut.NAVG*12.0) def SumLags(self,dataOut): for l in range(dataOut.DPL): dataOut.kabxys_integrated[4][:,l,0]=(dataOut.kabxys_integrated[4][:,l,0]+dataOut.kabxys_integrated[4][:,l,1])*dataOut.rnint2[l] dataOut.kabxys_integrated[5][:,l,0]=(dataOut.kabxys_integrated[5][:,l,0]+dataOut.kabxys_integrated[5][:,l,1])*dataOut.rnint2[l] dataOut.kabxys_integrated[6][:,l,0]=(dataOut.kabxys_integrated[6][:,l,0]+dataOut.kabxys_integrated[6][:,l,1])*dataOut.rnint2[l] dataOut.kabxys_integrated[7][:,l,0]=(dataOut.kabxys_integrated[7][:,l,0]+dataOut.kabxys_integrated[7][:,l,1])*dataOut.rnint2[l] dataOut.kabxys_integrated[8][:,l,0]=(dataOut.kabxys_integrated[8][:,l,0]-dataOut.kabxys_integrated[8][:,l,1])*dataOut.rnint2[l] dataOut.kabxys_integrated[9][:,l,0]=(dataOut.kabxys_integrated[9][:,l,0]-dataOut.kabxys_integrated[9][:,l,1])*dataOut.rnint2[l] dataOut.kabxys_integrated[10][:,l,0]=(dataOut.kabxys_integrated[10][:,l,0]-dataOut.kabxys_integrated[10][:,l,1])*dataOut.rnint2[l] dataOut.kabxys_integrated[11][:,l,0]=(dataOut.kabxys_integrated[11][:,l,0]-dataOut.kabxys_integrated[11][:,l,1])*dataOut.rnint2[l] def run(self,dataOut): self.rint2DP(dataOut) self.SumLags(dataOut) return dataOut class FlagBadHeights(Operation): """Operation to flag bad heights (bad data) of the Double Pulse. Parameters: ----------- None Example -------- op = proc_unit.addOperation(name='FlagBadHeights', optype='other') """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) def run(self,dataOut): dataOut.ibad=numpy.zeros((dataOut.NDP,dataOut.DPL),'int32') for j in range(dataOut.NDP): for l in range(dataOut.DPL): ip1=j+dataOut.NDP*(0+2*l) if( (dataOut.kabxys_integrated[5][j,l,0] <= 0.) or (dataOut.kabxys_integrated[4][j,l,0] <= 0.) or (dataOut.kabxys_integrated[7][j,l,0] <= 0.) or (dataOut.kabxys_integrated[6][j,l,0] <= 0.)): dataOut.ibad[j][l]=1 else: dataOut.ibad[j][l]=0 return dataOut class FlagBadHeightsSpectra(Operation): """Operation to flag bad heights (bad data) of the Double Pulse. Parameters: ----------- None Example -------- op = proc_unit.addOperation(name='FlagBadHeightsSpectra', optype='other') """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) def run(self,dataOut): dataOut.ibad=numpy.zeros((dataOut.NDP,dataOut.DPL),'int32') for j in range(dataOut.NDP): for l in range(dataOut.DPL): ip1=j+dataOut.NDP*(0+2*l) if( (dataOut.kabxys_integrated[4][j,l,0] <= 0.) or (dataOut.kabxys_integrated[6][j,l,0] <= 0.)): dataOut.ibad[j][l]=1 else: dataOut.ibad[j][l]=0 return dataOut class NoisePower(Operation): """Operation to get noise power from the integrated data of the Double Pulse. Parameters: ----------- None Example -------- op = proc_unit.addOperation(name='NoisePower', optype='other') """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) def hildebrand(self,dataOut,data): divider=10 # divider was originally 10 noise=0.0 n1=0 n2=int(dataOut.NDP/2) sorts= sorted(data) nums_min= dataOut.NDP/divider if((dataOut.NDP/divider)> 2): nums_min= int(dataOut.NDP/divider) else: nums_min=2 sump=0.0 sumq=0.0 j=0 cont=1 while( (cont==1) and (j nums_min): rtest= float(j/(j-1)) +1.0/dataOut.NAVG t1= (sumq*j) t2=(rtest*sump*sump) if( (t1/t2) > 0.990): j=j-1 sump-= sorts[j+n1] sumq-=sorts[j+n1]*sorts[j+n1] cont= 0 noise= sump/j stdv=numpy.sqrt((sumq- noise*noise)/(j-1)) return noise def run(self,dataOut): p=numpy.zeros((dataOut.NR,dataOut.NDP,dataOut.DPL),'float32') av=numpy.zeros(dataOut.NDP,'float32') dataOut.pnoise=numpy.zeros(dataOut.NR,'float32') p[0,:,:]=dataOut.kabxys_integrated[4][:,:,0]+dataOut.kabxys_integrated[5][:,:,0] #total power for channel 0, just pulse with non-flip p[1,:,:]=dataOut.kabxys_integrated[6][:,:,0]+dataOut.kabxys_integrated[7][:,:,0] #total power for channel 1 for i in range(dataOut.NR): dataOut.pnoise[i]=0.0 for k in range(dataOut.DPL): dataOut.pnoise[i]+= self.hildebrand(dataOut,p[i,:,k]) dataOut.pnoise[i]=dataOut.pnoise[i]/dataOut.DPL dataOut.pan=1.0*dataOut.pnoise[0] # weights could change dataOut.pbn=1.0*dataOut.pnoise[1] # weights could change return dataOut class DoublePulseACFs(Operation): """Operation to get the ACFs of the Double Pulse. Parameters: ----------- None Example -------- op = proc_unit.addOperation(name='DoublePulseACFs', optype='other') """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.aux=1 def run(self,dataOut): dataOut.igcej=numpy.zeros((dataOut.NDP,dataOut.DPL),'int32') if self.aux==1: dataOut.rhor=numpy.zeros((dataOut.NDP,dataOut.DPL), dtype=float) dataOut.rhoi=numpy.zeros((dataOut.NDP,dataOut.DPL), dtype=float) dataOut.sdp=numpy.zeros((dataOut.NDP,dataOut.DPL), dtype=float) dataOut.sd=numpy.zeros((dataOut.NDP,dataOut.DPL), dtype=float) dataOut.p=numpy.zeros((dataOut.NDP,dataOut.DPL), dtype=float) dataOut.alag=numpy.zeros(dataOut.NDP,'float32') for l in range(dataOut.DPL): dataOut.alag[l]=l*dataOut.DH*2.0/150.0 self.aux=0 sn4=dataOut.pan*dataOut.pbn rhorn=0 rhoin=0 panrm=numpy.zeros((dataOut.NDP,dataOut.DPL), dtype=float) for i in range(dataOut.NDP): for j in range(dataOut.DPL): ################# Total power pa=numpy.abs(dataOut.kabxys_integrated[4][i,j,0]+dataOut.kabxys_integrated[5][i,j,0]) pb=numpy.abs(dataOut.kabxys_integrated[6][i,j,0]+dataOut.kabxys_integrated[7][i,j,0]) st4=pa*pb dataOut.p[i,j]=pa+pb-(dataOut.pan+dataOut.pbn) dataOut.sdp[i,j]=2*dataOut.rnint2[j]*((pa+pb)*(pa+pb)) ## ACF rhorp=dataOut.kabxys_integrated[8][i,j,0]+dataOut.kabxys_integrated[11][i,j,0] rhoip=dataOut.kabxys_integrated[10][i,j,0]-dataOut.kabxys_integrated[9][i,j,0] if ((pa>dataOut.pan)&(pb>dataOut.pbn)): ss4=numpy.abs((pa-dataOut.pan)*(pb-dataOut.pbn)) panrm[i,j]=math.sqrt(ss4) rnorm=1/panrm[i,j] ## ACF dataOut.rhor[i,j]=rhorp*rnorm dataOut.rhoi[i,j]=rhoip*rnorm ############# Compute standard error for ACF stoss4=st4/ss4 snoss4=sn4/ss4 rp2=((rhorp*rhorp)+(rhoip*rhoip))/st4 rn2=((rhorn*rhorn)+(rhoin*rhoin))/sn4 rs2=(dataOut.rhor[i,j]*dataOut.rhor[i,j])+(dataOut.rhoi[i,j]*dataOut.rhoi[i,j]) st=1.0+rs2*(stoss4-(2*math.sqrt(stoss4*snoss4))) stn=1.0+rs2*(snoss4-(2*math.sqrt(stoss4*snoss4))) dataOut.sd[i,j]=((stoss4*((1.0+rp2)*st+(2.0*rp2*rs2*snoss4)-4.0*math.sqrt(rs2*rp2)))+(0.25*snoss4*((1.0+rn2)*stn+(2.0*rn2*rs2*stoss4)-4.0*math.sqrt(rs2*rn2))))*dataOut.rnint2[j] dataOut.sd[i,j]=numpy.abs(dataOut.sd[i,j]) else: #default values for bad points rnorm=1/math.sqrt(st4) dataOut.sd[i,j]=1.e30 dataOut.ibad[i,j]=4 dataOut.rhor[i,j]=rhorp*rnorm dataOut.rhoi[i,j]=rhoip*rnorm if ((pa/dataOut.pan-1.0)>2.25*(pb/dataOut.pbn-1.0)): dataOut.igcej[i,j]=1 return dataOut class FaradayAngleAndDPPower(Operation): """Operation to calculate Faraday angle and Double Pulse power. Parameters: ----------- None Example -------- op = proc_unit.addOperation(name='FaradayAngleAndDPPower', optype='other') """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.aux=1 def run(self,dataOut): if self.aux==1: dataOut.h2=numpy.zeros(dataOut.MAXNRANGENDT,'float32') dataOut.range1=numpy.zeros(dataOut.MAXNRANGENDT,order='F',dtype='float32') dataOut.sdn2=numpy.zeros(dataOut.NDP,'float32') dataOut.ph2=numpy.zeros(dataOut.NDP,'float32') dataOut.sdp2=numpy.zeros(dataOut.NDP,'float32') dataOut.ibd=numpy.zeros(dataOut.NDP,'float32') dataOut.phi=numpy.zeros(dataOut.NDP,'float32') self.aux=0 for i in range(dataOut.MAXNRANGENDT): dataOut.range1[i]=dataOut.H0 + i*dataOut.DH dataOut.h2[i]=dataOut.range1[i]**2 for j in range(dataOut.NDP): dataOut.ph2[j]=0. dataOut.sdp2[j]=0. ri=dataOut.rhoi[j][0]/dataOut.sd[j][0] rr=dataOut.rhor[j][0]/dataOut.sd[j][0] dataOut.sdn2[j]=1./dataOut.sd[j][0] pt=0.# // total power st=0.# // total signal ibt=0# // bad lags ns=0# // no. good lags for l in range(dataOut.DPL): #add in other lags if outside of e-jet contamination if( (dataOut.igcej[j][l] == 0) and (dataOut.ibad[j][l] == 0) ): dataOut.ph2[j]+=dataOut.p[j][l]/dataOut.sdp[j][l] dataOut.sdp2[j]=dataOut.sdp2[j]+1./dataOut.sdp[j][l] ns+=1 pt+=dataOut.p[j][l]/dataOut.sdp[j][l] st+=1./dataOut.sdp[j][l] ibt|=dataOut.ibad[j][l]; if(ns!= 0): dataOut.ibd[j]=0 dataOut.ph2[j]=dataOut.ph2[j]/dataOut.sdp2[j] dataOut.sdp2[j]=1./dataOut.sdp2[j] else: dataOut.ibd[j]=ibt dataOut.ph2[j]=pt/st dataOut.sdp2[j]=1./st dataOut.ph2[j]=dataOut.ph2[j]*dataOut.h2[j] dataOut.sdp2[j]=numpy.sqrt(dataOut.sdp2[j])*dataOut.h2[j] rr=rr/dataOut.sdn2[j] ri=ri/dataOut.sdn2[j] #rm[j]=np.sqrt(rr*rr + ri*ri) it is not used in c program dataOut.sdn2[j]=1./(dataOut.sdn2[j]*(rr*rr + ri*ri)) if( (ri == 0.) and (rr == 0.) ): dataOut.phi[j]=0. else: dataOut.phi[j]=math.atan2( ri , rr ) return dataOut class ElectronDensityFaraday(Operation): """Operation to calculate electron density from Faraday angle. Parameters: ----------- NSHTS : int .* RATE : float .* Example -------- op = proc_unit.addOperation(name='ElectronDensityFaraday', optype='other') op.addParameter(name='NSHTS', value='50', format='int') op.addParameter(name='RATE', value='1.8978873e-6', format='float') """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.aux=1 def run(self,dataOut,NSHTS=50,RATE=1.8978873e-6): #print(ctime(dataOut.utctime)) #3print("Faraday Angle",dataOut.phi) dataOut.NSHTS=NSHTS dataOut.RATE=RATE if self.aux==1: dataOut.dphi=numpy.zeros(dataOut.NDP,'float32') dataOut.sdn1=numpy.zeros(dataOut.NDP,'float32') self.aux=0 theta=numpy.zeros(dataOut.NDP,dtype=numpy.complex_) thetai=numpy.zeros(dataOut.NDP,dtype=numpy.complex_) # use complex numbers for phase for i in range(dataOut.NSHTS): theta[i]=math.cos(dataOut.phi[i])+math.sin(dataOut.phi[i])*1j thetai[i]=-math.sin(dataOut.phi[i])+math.cos(dataOut.phi[i])*1j # differentiate and convert to number density ndphi=dataOut.NSHTS-4 for i in range(2,dataOut.NSHTS-2): fact=(-0.5/(dataOut.RATE*dataOut.DH))*dataOut.bki[i] #four-point derivative, no phase unwrapping necessary ####dataOut.dphi[i]=((((theta[i+1]-theta[i-1])+(2.0*(theta[i+2]-theta[i-2])))/thetai[i])).real/10.0 dataOut.dphi[i]=((((theta[i-2]-theta[i+2])+(8.0*(theta[i+1]-theta[i-1])))/thetai[i])).real/12.0 dataOut.dphi[i]=abs(dataOut.dphi[i]*fact) dataOut.sdn1[i]=(4.*(dataOut.sdn2[i-2]+dataOut.sdn2[i+2])+dataOut.sdn2[i-1]+dataOut.sdn2[i+1]) dataOut.sdn1[i]=numpy.sqrt(dataOut.sdn1[i])*fact return dataOut class ElectronDensityRobertoTestFaraday(Operation): """Operation to calculate electron density from Faraday angle. Parameters: ----------- NSHTS : int .* RATE : float .* Example -------- op = proc_unit.addOperation(name='ElectronDensityFaraday', optype='other') op.addParameter(name='NSHTS', value='50', format='int') op.addParameter(name='RATE', value='1.8978873e-6', format='float') """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.aux=1 def run(self,dataOut,NSHTS=50,RATE=1.8978873e-6): #print(ctime(dataOut.utctime)) #print("Faraday Angle",dataOut.phi) dataOut.NSHTS=NSHTS dataOut.RATE=RATE if self.aux==1: dataOut.dphi=numpy.zeros(dataOut.NDP,'float32') dataOut.sdn1=numpy.zeros(dataOut.NDP,'float32') self.aux=0 theta=numpy.zeros(dataOut.NDP,dtype=numpy.complex_) thetai=numpy.zeros(dataOut.NDP,dtype=numpy.complex_) # use complex numbers for phase ''' for i in range(dataOut.NSHTS): theta[i]=math.cos(dataOut.phi[i])+math.sin(dataOut.phi[i])*1j thetai[i]=-math.sin(dataOut.phi[i])+math.cos(dataOut.phi[i])*1j ''' # differentiate and convert to number density ndphi=dataOut.NSHTS-4 dataOut.phi=numpy.unwrap(dataOut.phi) for i in range(2,dataOut.NSHTS-2): fact=(-0.5/(dataOut.RATE*dataOut.DH))*dataOut.bki[i] #four-point derivative, no phase unwrapping necessary ####dataOut.dphi[i]=((((theta[i+1]-theta[i-1])+(2.0*(theta[i+2]-theta[i-2])))/thetai[i])).real/10.0 ##dataOut.dphi[i]=((((theta[i-2]-theta[i+2])+(8.0*(theta[i+1]-theta[i-1])))/thetai[i])).real/12.0 dataOut.dphi[i]=((dataOut.phi[i+1]-dataOut.phi[i-1])+(2.0*(dataOut.phi[i+2]-dataOut.phi[i-2])))/10.0 dataOut.dphi[i]=abs(dataOut.dphi[i]*fact) dataOut.sdn1[i]=(4.*(dataOut.sdn2[i-2]+dataOut.sdn2[i+2])+dataOut.sdn2[i-1]+dataOut.sdn2[i+1]) dataOut.sdn1[i]=numpy.sqrt(dataOut.sdn1[i])*fact return dataOut class ElectronDensityRobertoTest2Faraday(Operation): """Operation to calculate electron density from Faraday angle. Parameters: ----------- NSHTS : int .* RATE : float .* Example -------- op = proc_unit.addOperation(name='ElectronDensityFaraday', optype='other') op.addParameter(name='NSHTS', value='50', format='int') op.addParameter(name='RATE', value='1.8978873e-6', format='float') """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.aux=1 def run(self,dataOut,NSHTS=50,RATE=1.8978873e-6): #print(ctime(dataOut.utctime)) #print("Faraday Angle",dataOut.phi) dataOut.NSHTS=NSHTS dataOut.RATE=RATE if self.aux==1: dataOut.dphi=numpy.zeros(dataOut.NDP,'float32') dataOut.sdn1=numpy.zeros(dataOut.NDP,'float32') self.aux=0 theta=numpy.zeros(dataOut.NDP,dtype=numpy.complex_) thetai=numpy.zeros(dataOut.NDP,dtype=numpy.complex_) # use complex numbers for phase ''' for i in range(dataOut.NSHTS): theta[i]=math.cos(dataOut.phi[i])+math.sin(dataOut.phi[i])*1j thetai[i]=-math.sin(dataOut.phi[i])+math.cos(dataOut.phi[i])*1j ''' # differentiate and convert to number density ndphi=dataOut.NSHTS-4 #dataOut.phi=numpy.unwrap(dataOut.phi) f1=numpy.exp((dataOut.phi*1.j)/10) f2=numpy.exp((dataOut.phi*2.j)/10) for i in range(2,dataOut.NSHTS-2): fact=(-0.5/(dataOut.RATE*dataOut.DH))*dataOut.bki[i] #four-point derivative, no phase unwrapping necessary ####dataOut.dphi[i]=((((theta[i+1]-theta[i-1])+(2.0*(theta[i+2]-theta[i-2])))/thetai[i])).real/10.0 ##dataOut.dphi[i]=((((theta[i-2]-theta[i+2])+(8.0*(theta[i+1]-theta[i-1])))/thetai[i])).real/12.0 ##dataOut.dphi[i]=((dataOut.phi[i+1]-dataOut.phi[i-1])+(2.0*(dataOut.phi[i+2]-dataOut.phi[i-2])))/10.0 dataOut.dphi[i]=numpy.angle(f1[i+1]*numpy.conjugate(f1[i-1])*f2[i+2]*numpy.conjugate(f2[i-2])) dataOut.dphi[i]=abs(dataOut.dphi[i]*fact) dataOut.sdn1[i]=(4.*(dataOut.sdn2[i-2]+dataOut.sdn2[i+2])+dataOut.sdn2[i-1]+dataOut.sdn2[i+1]) dataOut.sdn1[i]=numpy.sqrt(dataOut.sdn1[i])*fact return dataOut class NormalizeDPPower(Operation): """Operation to normalize relative electron density from power with total electron density from Farday angle. Parameters: ----------- None Example -------- op = proc_unit.addOperation(name='NormalizeDPPower', optype='other') """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.aux=1 def normal(self,a,b,n,m): chmin=1.0e30 chisq=numpy.zeros(150,'float32') temp=numpy.zeros(150,'float32') for i in range(2*m-1): an=al=be=chisq[i]=0.0 for j in range(int(n/m)): k=int(j+i*n/(2*m)) if(a[k]>0.0 and b[k]>0.0): al+=a[k]*b[k] be+=b[k]*b[k] if(be>0.0): temp[i]=al/be else: temp[i]=1.0 for j in range(int(n/m)): k=int(j+i*n/(2*m)) if(a[k]>0.0 and b[k]>0.0): chisq[i]+=(numpy.log10(b[k]*temp[i]/a[k]))**2 an=an+1 if(chisq[i]>0.0): chisq[i]/=an for i in range(int(2*m-1)): if(chisq[i]1.0e-6): chmin=chisq[i] cf=temp[i] return cf def normalize(self,dataOut): if self.aux==1: dataOut.cf=numpy.zeros(1,'float32') dataOut.cflast=numpy.zeros(1,'float32') self.aux=0 night_first=300.0 night_first1= 310.0 night_end= 450.0 day_first=250.0 day_end=400.0 day_first_sunrise=190.0 day_end_sunrise=280.0 print(dataOut.ut_Faraday) if(dataOut.ut_Faraday>4.0 and dataOut.ut_Faraday<11.0): #early print("EARLY") i2=(night_end-dataOut.range1[0])/dataOut.DH i1=(night_first -dataOut.range1[0])/dataOut.DH elif (dataOut.ut_Faraday>0.0 and dataOut.ut_Faraday<4.0): #night print("NIGHT") i2=(night_end-dataOut.range1[0])/dataOut.DH i1=(night_first1 -dataOut.range1[0])/dataOut.DH elif (dataOut.ut_Faraday>=11.0 and dataOut.ut_Faraday<13.5): #sunrise print("SUNRISE") i2=( day_end_sunrise-dataOut.range1[0])/dataOut.DH i1=(day_first_sunrise - dataOut.range1[0])/dataOut.DH else: print("ELSE") i2=(day_end-dataOut.range1[0])/dataOut.DH i1=(day_first -dataOut.range1[0])/dataOut.DH #print(i1*dataOut.DH) #print(i2*dataOut.DH) i1=int(i1) i2=int(i2) try: dataOut.cf=self.normal(dataOut.dphi[i1::], dataOut.ph2[i1::], i2-i1, 1) except: pass #print(dataOut.ph2) #input() # in case of spread F, normalize much higher if(dataOut.cf0.0 and b[k]>0.0): al+=a[k]*b[k] be+=b[k]*b[k] if(be>0.0): temp[i]=al/be else: temp[i]=1.0 for j in range(int(n/m)): k=int(j+i*n/(2*m)) if(a[k]>0.0 and b[k]>0.0): chisq[i]+=(numpy.log10(b[k]*temp[i]/a[k]))**2 an=an+1 if(chisq[i]>0.0): chisq[i]/=an for i in range(int(2*m-1)): if(chisq[i]1.0e-6): chmin=chisq[i] cf=temp[i] return cf def normalize(self,dataOut): if self.aux==1: dataOut.cf=numpy.zeros(1,'float32') dataOut.cflast=numpy.zeros(1,'float32') self.aux=0 night_first=300.0 night_first1= 310.0 night_end= 450.0 day_first=250.0 day_end=400.0 day_first_sunrise=190.0 day_end_sunrise=350.0 print(dataOut.ut_Faraday) ''' if(dataOut.ut_Faraday>4.0 and dataOut.ut_Faraday<11.0): #early print("EARLY") i2=(night_end-dataOut.range1[0])/dataOut.DH i1=(night_first -dataOut.range1[0])/dataOut.DH elif (dataOut.ut_Faraday>0.0 and dataOut.ut_Faraday<4.0): #night print("NIGHT") i2=(night_end-dataOut.range1[0])/dataOut.DH i1=(night_first1 -dataOut.range1[0])/dataOut.DH elif (dataOut.ut_Faraday>=11.0 and dataOut.ut_Faraday<13.5): #sunrise print("SUNRISE") i2=( day_end_sunrise-dataOut.range1[0])/dataOut.DH i1=(day_first_sunrise - dataOut.range1[0])/dataOut.DH else: print("ELSE") i2=(day_end-dataOut.range1[0])/dataOut.DH i1=(day_first -dataOut.range1[0])/dataOut.DH ''' i2=(420-dataOut.range1[0])/dataOut.DH i1=(200 -dataOut.range1[0])/dataOut.DH print(i1*dataOut.DH) print(i2*dataOut.DH) i1=int(i1) i2=int(i2) try: dataOut.cf=self.normal(dataOut.dphi[i1::], dataOut.ph2[i1::], i2-i1, 1) except: pass #print(dataOut.ph2) #input() # in case of spread F, normalize much higher if(dataOut.cf10 and l1>=0: if l1==0: l1=1 dataOut.cov=numpy.reshape(dataOut.cov,l1*l1) dataOut.cov=numpy.resize(dataOut.cov,dataOut.DPL*dataOut.DPL) dataOut.covinv=numpy.reshape(dataOut.covinv,l1*l1) dataOut.covinv=numpy.resize(dataOut.covinv,dataOut.DPL*dataOut.DPL) for l in range(dataOut.DPL*dataOut.DPL): dataOut.cov[l]=0.0 acfm= (dataOut.rhor[i][0])**2 + (dataOut.rhoi[i][0])**2 if acfm> 0.0: cc=dataOut.rhor[i][0]/acfm ss=dataOut.rhoi[i][0]/acfm else: cc=1. ss=0. # keep only uncontaminated data, don't pass zero lag to fitter l1=0 for l in range(0+1,dataOut.DPL): if dataOut.igcej[i][l]==0 and dataOut.ibad[i][l]==0: y[l1]=dataOut.rhor[i][l]*cc + dataOut.rhoi[i][l]*ss x[l1]=dataOut.alag[l]*1.0e-3 dataOut.sd[i][l]=dataOut.sd[i][l]/((acfm)**2)# important e[l1]=dataOut.sd[i][l] #this is the variance, not the st. dev. l1=l1+1 for l in range(l1*(l1+1)): dataOut.cov[l]=0.0 for l in range(l1): dataOut.cov[l*(1+l1)]=e[l] angle=dataOut.thb[i]*0.01745 bm=dataOut.bfm[i] dataOut.params[0]=1.0 #norm dataOut.params[1]=1000.0 #te dataOut.params[2]=800.0 #ti dataOut.params[3]=0.00 #ph dataOut.params[4]=0.00 #phe if l1!=0: x=numpy.resize(x,l1) y=numpy.resize(y,l1) else: x=numpy.resize(x,1) y=numpy.resize(y,1) if True: #len(y)!=0: fitacf_guess.guess(y,x,zero,depth,t1,t2,len(y)) t2=t1/t2 if (t1<5000.0 and t1> 600.0): dataOut.params[1]=t1 dataOut.params[2]=min(t2,t1) dataOut.ifit[1]=dataOut.ifit[2]=1 dataOut.ifit[0]=dataOut.ifit[3]=dataOut.ifit[4]=0 #print(dataOut.ut_Faraday) if dataOut.ut_Faraday<10.0 and dataOut.ut_Faraday>=0.5: dataOut.ifit[2]=0 den=dataOut.ph2[i] if l1!=0: dataOut.covinv=dataOut.covinv[0:l1*l1].reshape((l1,l1)) dataOut.cov=dataOut.cov[0:l1*l1].reshape((l1,l1)) e=numpy.resize(e,l1) else: dataOut.covinv=numpy.resize(dataOut.covinv,1) dataOut.cov=numpy.resize(dataOut.cov,1) e=numpy.resize(e,1) eb=numpy.resize(eb,10) dataOut.ifit=numpy.resize(dataOut.ifit,10) dataOut.covinv,e,dataOut.params,eb,dataOut.m=fitacf_fit_short.fit(wl,x,y,dataOut.cov,dataOut.covinv,e,dataOut.params,bm,angle,den,dataOut.range1[i],dataOut.year,dataOut.ifit,dataOut.m,l1) # #exit() if dataOut.params[2]>dataOut.params[1]*1.05: dataOut.ifit[2]=0 dataOut.params[1]=dataOut.params[2]=t1 dataOut.covinv,e,dataOut.params,eb,dataOut.m=fitacf_fit_short.fit(wl,x,y,dataOut.cov,dataOut.covinv,e,dataOut.params,bm,angle,den,dataOut.range1[i],dataOut.year,dataOut.ifit,dataOut.m,l1) # if (dataOut.ifit[2]==0): dataOut.params[2]=dataOut.params[1] if (dataOut.ifit[3]==0 and iflag==0): dataOut.params[3]=0.0 if (dataOut.ifit[4]==0): dataOut.params[4]=0.0 dataOut.te2[i]=dataOut.params[1] dataOut.ti2[i]=dataOut.params[2] dataOut.ete2[i]=eb[1] dataOut.eti2[i]=eb[2] if dataOut.eti2[i]==0: dataOut.eti2[i]=dataOut.ete2[i] dataOut.phy2[i]=dataOut.params[3] dataOut.ephy2[i]=eb[3] if(iflag==1): dataOut.ephy2[i]=0.0 if (dataOut.m<=3 and dataOut.m!= 0 and dataOut.te2[i]>400.0): dataOut.info2[i]=1 else: dataOut.info2[i]=0 def run(self,dataOut,IBITS=16): dataOut.IBITS = IBITS self.Estimation(dataOut) return dataOut class NeTeTiRecal(NormalizeDPPower,DPTemperaturesEstimation): def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.aux=0 def run(self,dataOut): for i in range(dataOut.NSHTS): print("H: ",i*15) print(1+(dataOut.te2[i]/dataOut.ti2[i])) dataOut.ph2[i]*=1+(dataOut.te2[i]/dataOut.ti2[i]) self.normalize(dataOut) self.Estimation(dataOut) return dataOut from schainpy.model.proc import fitacf_acf2 class DenCorrection(Operation): def __init__(self, **kwargs): Operation.__init__(self, **kwargs) def run(self,dataOut): y=numpy.zeros(dataOut.DPL,order='F',dtype='float32') #y_aux = numpy.zeros(1,,dtype='float32') for i in range(dataOut.NSHTS): y[0]=y[1]=dataOut.range1[i] y = y.astype(dtype='float64',order='F') three=int(3) wl = 3.0 tion=numpy.zeros(three,order='F',dtype='float32') fion=numpy.zeros(three,order='F',dtype='float32') nui=numpy.zeros(three,order='F',dtype='float32') wion=numpy.zeros(three,order='F',dtype='int32') bline=0.0 #bline=numpy.zeros(1,order='F',dtype='float32') #print("**** ACF2 WRAPPER ***** ",fitacf_acf2.acf2.__doc__ ) print("BEFORE",dataOut.ph2[10:35]) for i in range(dataOut.NSHTS): if dataOut.info2[i]==1: angle=dataOut.thb[i]*0.01745 nue=nui[0]=nui[1]=nui[2]=0.0#nui[3]=0.0 wion[0]=16 wion[1]=1 wion[2]=4 tion[0]=tion[1]=tion[2]=dataOut.ti2[i] fion[0]=1.0-dataOut.phy2[i] fion[1]=dataOut.phy2[i] fion[2]=0.0 for j in range(dataOut.DPL): tau=dataOut.alag[j]*1.0e-3 ''' print("**** input from acf2 ***** ") print("wl ",wl) print("tau ",tau) print("te2[i] ",dataOut.te2[i]) print("tion ",tion) print("fion ",fion) print("nue ",nue) print("nui ",nui) print("wion ",wion) print("angle ",angle) print("ph2[i] ",dataOut.ph2[i]) print("bfm[i] ",dataOut.bfm[i]) print("y[j] ",y[j]) ''' print("Before y[j] ",y[j]) #with suppress_stdout_stderr(): y[j]=fitacf_acf2.acf2(wl,tau,dataOut.te2[i],tion,fion,nue,nui,wion,angle,dataOut.ph2[i],dataOut.bfm[i],y[j],three) #print("l",l) print("After y[j] ",y[j]) ''' print("**** output from acf2 ***** ") print("wl ",wl) print("tau ",tau) print("te2[i] ",dataOut.te2[i]) print("tion ",tion) print("fion ",fion) print("nue ",nue) print("nui ",nui) print("wion ",wion) print("angle ",angle) print("ph2[i] ",dataOut.ph2[i]) print("bfm[i] ",dataOut.bfm[i]) print("y[j] ",y[j]) print("i ",i , " j ",j , "y[j] ",y[j]) ''' #exit(1) if dataOut.ut_Faraday>11.0 and dataOut.range1[i]>150.0 and dataOut.range1[i]<400.0: tau=0.0 #with suppress_stdout_stderr(): bline=fitacf_acf2.acf2(wl,tau,tion,tion,fion,nue,nui,wion,angle,dataOut.ph2[i],dataOut.bfm[i],bline,three) cf=min(1.2,max(1.0,bline/y[0])) print("bline: ",bline) if cf != 1.0: print("bline: ",bline) print("cf: ",cf) #exit(1) #print("cf: ",cf) dataOut.ph2[i]=cf*dataOut.ph2[i] dataOut.sdp2[i]=cf*dataOut.sdp2[i] for j in range(1,dataOut.DPL): y[j]=(y[j]/y[0])*dataOut.DH+dataOut.range1[i] y[0]=dataOut.range1[i]+dataOut.DH #exit(1) #exit(1) print("AFTER",dataOut.ph2[10:35]) #exit(1) return dataOut class DataPlotCleaner(Operation): def __init__(self, **kwargs): Operation.__init__(self, **kwargs) def run(self,dataOut): THRESH_MIN_POW=10000 THRESH_MAX_POW=10000000 THRESH_MIN_TEMP=500 THRESH_MAX_TEMP=4000 dataOut.DensityClean=numpy.zeros((1,dataOut.NDP)) dataOut.EDensityClean=numpy.zeros((1,dataOut.NDP)) dataOut.ElecTempClean=numpy.zeros((1,dataOut.NDP)) dataOut.EElecTempClean=numpy.zeros((1,dataOut.NDP)) dataOut.IonTempClean=numpy.zeros((1,dataOut.NDP)) dataOut.EIonTempClean=numpy.zeros((1,dataOut.NDP)) dataOut.DensityClean[0]=numpy.copy(dataOut.ph2) dataOut.EDensityClean[0]=numpy.copy(dataOut.sdp2) dataOut.ElecTempClean[0,:dataOut.NSHTS]=numpy.copy(dataOut.te2) dataOut.EElecTempClean[0,:dataOut.NSHTS]=numpy.copy(dataOut.ete2) dataOut.IonTempClean[0,:dataOut.NSHTS]=numpy.copy(dataOut.ti2) dataOut.EIonTempClean[0,:dataOut.NSHTS]=numpy.copy(dataOut.eti2) for i in range(dataOut.NDP): if dataOut.DensityClean[0,i]THRESH_MAX_POW: dataOut.DensityClean[0,i]=THRESH_MAX_POW for i in range(dataOut.NSHTS): dataOut.ElecTempClean[0,i]=(max(1.0, dataOut.ElecTempClean[0,i])) dataOut.IonTempClean[0,i]=(max(1.0, dataOut.IonTempClean[0,i])) for i in range(dataOut.NSHTS): if dataOut.ElecTempClean[0,i]THRESH_MAX_TEMP: dataOut.ElecTempClean[0,i]=THRESH_MAX_TEMP if dataOut.IonTempClean[0,i]>THRESH_MAX_TEMP: dataOut.IonTempClean[0,i]=THRESH_MAX_TEMP for i in range(dataOut.NSHTS): if dataOut.EElecTempClean[0,i]>500:# dataOut.ElecTempClean[0,i]=500 if dataOut.EIonTempClean[0,i]>500:# dataOut.IonTempClean[0,i]=500 missing=numpy.nan for i in range(dataOut.NSHTS,dataOut.NDP): dataOut.ElecTempClean[0,i]=missing dataOut.EElecTempClean[0,i]=missing dataOut.IonTempClean[0,i]=missing dataOut.EIonTempClean[0,i]=missing return dataOut class DataSaveCleaner(Operation): def __init__(self, **kwargs): Operation.__init__(self, **kwargs) def run(self,dataOut): dataOut.DensityFinal=numpy.zeros((1,dataOut.NDP)) dataOut.EDensityFinal=numpy.zeros((1,dataOut.NDP)) dataOut.ElecTempFinal=numpy.zeros((1,dataOut.NDP)) dataOut.EElecTempFinal=numpy.zeros((1,dataOut.NDP)) dataOut.IonTempFinal=numpy.zeros((1,dataOut.NDP)) dataOut.EIonTempFinal=numpy.zeros((1,dataOut.NDP)) dataOut.PhyFinal=numpy.zeros((1,dataOut.NDP)) dataOut.EPhyFinal=numpy.zeros((1,dataOut.NDP)) dataOut.DensityFinal[0]=numpy.copy(dataOut.ph2) dataOut.EDensityFinal[0]=numpy.copy(dataOut.sdp2) dataOut.ElecTempFinal[0,:dataOut.NSHTS]=numpy.copy(dataOut.te2) dataOut.EElecTempFinal[0,:dataOut.NSHTS]=numpy.copy(dataOut.ete2) dataOut.IonTempFinal[0,:dataOut.NSHTS]=numpy.copy(dataOut.ti2) dataOut.EIonTempFinal[0,:dataOut.NSHTS]=numpy.copy(dataOut.eti2) dataOut.PhyFinal[0,:dataOut.NSHTS]=numpy.copy(dataOut.phy2) dataOut.EPhyFinal[0,:dataOut.NSHTS]=numpy.copy(dataOut.ephy2) missing=numpy.nan temp_min=100.0 temp_max=3000.0#6000.0e for i in range(dataOut.NSHTS): if dataOut.info2[i]!=1: dataOut.ElecTempFinal[0,i]=dataOut.EElecTempFinal[0,i]=dataOut.IonTempFinal[0,i]=dataOut.EIonTempFinal[0,i]=missing if dataOut.ElecTempFinal[0,i]<=temp_min or dataOut.ElecTempFinal[0,i]>temp_max or dataOut.EElecTempFinal[0,i]>temp_max: dataOut.ElecTempFinal[0,i]=dataOut.EElecTempFinal[0,i]=missing if dataOut.IonTempFinal[0,i]<=temp_min or dataOut.IonTempFinal[0,i]>temp_max or dataOut.EIonTempFinal[0,i]>temp_max: dataOut.IonTempFinal[0,i]=dataOut.EIonTempFinal[0,i]=missing if dataOut.lags_to_plot[i,:][~numpy.isnan(dataOut.lags_to_plot[i,:])].shape[0]<6: dataOut.ElecTempFinal[0,i]=dataOut.EElecTempFinal[0,i]=dataOut.IonTempFinal[0,i]=dataOut.EIonTempFinal[0,i]=missing if dataOut.ut_Faraday>4 and dataOut.ut_Faraday<11: if numpy.nanmax(dataOut.acfs_error_to_plot[i,:])>=10: dataOut.ElecTempFinal[0,i]=dataOut.EElecTempFinal[0,i]=dataOut.IonTempFinal[0,i]=dataOut.EIonTempFinal[0,i]=missing if dataOut.EPhyFinal[0,i]<0.0 or dataOut.EPhyFinal[0,i]>1.0: dataOut.PhyFinal[0,i]=dataOut.EPhyFinal[0,i]=missing if dataOut.EDensityFinal[0,i]>0.0 and dataOut.DensityFinal[0,i]>0.0 and dataOut.DensityFinal[0,i]<9.9e6: dataOut.EDensityFinal[0,i]=max(dataOut.EDensityFinal[0,i],1000.0) else: dataOut.DensityFinal[0,i]=dataOut.EDensityFinal[0,i]=missing if dataOut.PhyFinal[0,i]==0 or dataOut.PhyFinal[0,i]>0.4: dataOut.PhyFinal[0,i]=dataOut.EPhyFinal[0,i]=missing if dataOut.ElecTempFinal[0,i]==dataOut.IonTempFinal[0,i]: dataOut.EElecTempFinal[0,i]=dataOut.EIonTempFinal[0,i] if numpy.isnan(dataOut.ElecTempFinal[0,i]): dataOut.EElecTempFinal[0,i]=missing if numpy.isnan(dataOut.IonTempFinal[0,i]): dataOut.EIonTempFinal[0,i]=missing if numpy.isnan(dataOut.ElecTempFinal[0,i]) or numpy.isnan(dataOut.EElecTempFinal[0,i]): dataOut.ElecTempFinal[0,i]=dataOut.EElecTempFinal[0,i]=dataOut.IonTempFinal[0,i]=dataOut.EIonTempFinal[0,i]=missing for i in range(12,dataOut.NSHTS-1): if numpy.isnan(dataOut.ElecTempFinal[0,i-1]) and numpy.isnan(dataOut.ElecTempFinal[0,i+1]): dataOut.ElecTempFinal[0,i]=dataOut.EElecTempFinal[0,i]=missing if numpy.isnan(dataOut.IonTempFinal[0,i-1]) and numpy.isnan(dataOut.IonTempFinal[0,i+1]): dataOut.IonTempFinal[0,i]=dataOut.EIonTempFinal[0,i]=missing if dataOut.ut_Faraday>4 and dataOut.ut_Faraday<11: if numpy.isnan(dataOut.ElecTempFinal[0,i-1]) and numpy.isnan(dataOut.ElecTempFinal[0,i-2]) and numpy.isnan(dataOut.ElecTempFinal[0,i+2]) and numpy.isnan(dataOut.ElecTempFinal[0,i+3]): #and numpy.isnan(dataOut.ElecTempFinal[0,i-5]): dataOut.ElecTempFinal[0,i]=dataOut.EElecTempFinal[0,i]=missing if numpy.isnan(dataOut.IonTempFinal[0,i-1]) and numpy.isnan(dataOut.IonTempFinal[0,i-2]) and numpy.isnan(dataOut.IonTempFinal[0,i+2]) and numpy.isnan(dataOut.IonTempFinal[0,i+3]): #and numpy.isnan(dataOut.IonTempFinal[0,i-5]): dataOut.IonTempFinal[0,i]=dataOut.EIonTempFinal[0,i]=missing if i>25: if numpy.isnan(dataOut.ElecTempFinal[0,i-1]) and numpy.isnan(dataOut.ElecTempFinal[0,i-2]) and numpy.isnan(dataOut.ElecTempFinal[0,i-3]) and numpy.isnan(dataOut.ElecTempFinal[0,i-4]): #and numpy.isnan(dataOut.ElecTempFinal[0,i-5]): dataOut.ElecTempFinal[0,i]=dataOut.EElecTempFinal[0,i]=missing if numpy.isnan(dataOut.IonTempFinal[0,i-1]) and numpy.isnan(dataOut.IonTempFinal[0,i-2]) and numpy.isnan(dataOut.IonTempFinal[0,i-3]) and numpy.isnan(dataOut.IonTempFinal[0,i-4]): #and numpy.isnan(dataOut.IonTempFinal[0,i-5]): dataOut.IonTempFinal[0,i]=dataOut.EIonTempFinal[0,i]=missing if numpy.isnan(dataOut.ElecTempFinal[0,i]) or numpy.isnan(dataOut.EElecTempFinal[0,i]): dataOut.ElecTempFinal[0,i]=dataOut.EElecTempFinal[0,i]=dataOut.IonTempFinal[0,i]=dataOut.EIonTempFinal[0,i]=missing for i in range(12,dataOut.NSHTS-1): if numpy.isnan(dataOut.ElecTempFinal[0,i-1]) and numpy.isnan(dataOut.ElecTempFinal[0,i+1]): dataOut.ElecTempFinal[0,i]=dataOut.EElecTempFinal[0,i]=missing if numpy.isnan(dataOut.IonTempFinal[0,i-1]) and numpy.isnan(dataOut.IonTempFinal[0,i+1]): dataOut.IonTempFinal[0,i]=dataOut.EIonTempFinal[0,i]=missing if numpy.isnan(dataOut.ElecTempFinal[0,i]) or numpy.isnan(dataOut.EElecTempFinal[0,i]): dataOut.ElecTempFinal[0,i]=dataOut.EElecTempFinal[0,i]=dataOut.IonTempFinal[0,i]=dataOut.EIonTempFinal[0,i]=missing if numpy.count_nonzero(~numpy.isnan(dataOut.ElecTempFinal[0,12:50]))<5: dataOut.ElecTempFinal[0,:]=dataOut.EElecTempFinal[0,:]=missing if numpy.count_nonzero(~numpy.isnan(dataOut.IonTempFinal[0,12:50]))<5: dataOut.IonTempFinal[0,:]=dataOut.EIonTempFinal[0,:]=missing for i in range(dataOut.NSHTS,dataOut.NDP): dataOut.ElecTempFinal[0,i]=missing dataOut.EElecTempFinal[0,i]=missing dataOut.IonTempFinal[0,i]=missing dataOut.EIonTempFinal[0,i]=missing dataOut.PhyFinal[0,i]=missing dataOut.EPhyFinal[0,i]=missing return dataOut class DataSaveCleanerHP(Operation): def __init__(self, **kwargs): Operation.__init__(self, **kwargs) def run(self,dataOut): dataOut.Density_DP=numpy.zeros(dataOut.cut) dataOut.EDensity_DP=numpy.zeros(dataOut.cut) dataOut.ElecTemp_DP=numpy.zeros(dataOut.cut) dataOut.EElecTemp_DP=numpy.zeros(dataOut.cut) dataOut.IonTemp_DP=numpy.zeros(dataOut.cut) dataOut.EIonTemp_DP=numpy.zeros(dataOut.cut) dataOut.Phy_DP=numpy.zeros(dataOut.cut) dataOut.EPhy_DP=numpy.zeros(dataOut.cut) dataOut.Phe_DP=numpy.empty(dataOut.cut) dataOut.EPhe_DP=numpy.empty(dataOut.cut) dataOut.Density_DP[:]=numpy.copy(dataOut.ph2[:dataOut.cut]) dataOut.EDensity_DP[:]=numpy.copy(dataOut.sdp2[:dataOut.cut]) dataOut.ElecTemp_DP[:]=numpy.copy(dataOut.te2[:dataOut.cut]) dataOut.EElecTemp_DP[:]=numpy.copy(dataOut.ete2[:dataOut.cut]) dataOut.IonTemp_DP[:]=numpy.copy(dataOut.ti2[:dataOut.cut]) dataOut.EIonTemp_DP[:]=numpy.copy(dataOut.eti2[:dataOut.cut]) dataOut.Phy_DP[:]=numpy.copy(dataOut.phy2[:dataOut.cut]) dataOut.EPhy_DP[:]=numpy.copy(dataOut.ephy2[:dataOut.cut]) dataOut.Phe_DP[:]=numpy.nan dataOut.EPhe_DP[:]=numpy.nan missing=numpy.nan temp_min=100.0 temp_max_dp=3000.0 for i in range(dataOut.cut): if dataOut.info2[i]!=1: dataOut.ElecTemp_DP[i]=dataOut.EElecTemp_DP[i]=dataOut.IonTemp_DP[i]=dataOut.EIonTemp_DP[i]=missing if dataOut.ElecTemp_DP[i]<=temp_min or dataOut.ElecTemp_DP[i]>temp_max_dp or dataOut.EElecTemp_DP[i]>temp_max_dp: dataOut.ElecTemp_DP[i]=dataOut.EElecTemp_DP[i]=missing if dataOut.IonTemp_DP[i]<=temp_min or dataOut.IonTemp_DP[i]>temp_max_dp or dataOut.EIonTemp_DP[i]>temp_max_dp: dataOut.IonTemp_DP[i]=dataOut.EIonTemp_DP[i]=missing ####################################################################################### CHECK THIS if dataOut.lags_to_plot[i,:][~numpy.isnan(dataOut.lags_to_plot[i,:])].shape[0]<6: dataOut.ElecTemp_DP[i]=dataOut.EElecTemp_DP[i]=dataOut.IonTemp_DP[i]=dataOut.EIonTemp_DP[i]=missing if dataOut.ut_Faraday>4 and dataOut.ut_Faraday<11: if numpy.nanmax(dataOut.acfs_error_to_plot[i,:])>=10: dataOut.ElecTemp_DP[i]=dataOut.EElecTemp_DP[i]=dataOut.IonTemp_DP[i]=dataOut.EIonTemp_DP[i]=missing ####################################################################################### if dataOut.EPhy_DP[i]<0.0 or dataOut.EPhy_DP[i]>1.0: dataOut.Phy_DP[i]=dataOut.EPhy_DP[i]=missing if dataOut.EDensity_DP[i]>0.0 and dataOut.Density_DP[i]>0.0 and dataOut.Density_DP[i]<9.9e6: dataOut.EDensity_DP[i]=max(dataOut.EDensity_DP[i],1000.0) else: dataOut.Density_DP[i]=dataOut.EDensity_DP[i]=missing if dataOut.Phy_DP[i]==0 or dataOut.Phy_DP[i]>0.4: dataOut.Phy_DP[i]=dataOut.EPhy_DP[i]=missing if dataOut.ElecTemp_DP[i]==dataOut.IonTemp_DP[i]: dataOut.EElecTemp_DP[i]=dataOut.EIonTemp_DP[i] if numpy.isnan(dataOut.ElecTemp_DP[i]): dataOut.EElecTemp_DP[i]=missing if numpy.isnan(dataOut.IonTemp_DP[i]): dataOut.EIonTemp_DP[i]=missing if numpy.isnan(dataOut.ElecTemp_DP[i]) or numpy.isnan(dataOut.EElecTemp_DP[i]): dataOut.ElecTemp_DP[i]=dataOut.EElecTemp_DP[i]=dataOut.IonTemp_DP[i]=dataOut.EIonTemp_DP[i]=missing dataOut.Density_LP=numpy.zeros(dataOut.NACF-dataOut.cut) dataOut.EDensity_LP=numpy.zeros(dataOut.NACF-dataOut.cut) dataOut.ElecTemp_LP=numpy.zeros(dataOut.NACF-dataOut.cut) dataOut.EElecTemp_LP=numpy.zeros(dataOut.NACF-dataOut.cut) dataOut.IonTemp_LP=numpy.zeros(dataOut.NACF-dataOut.cut) dataOut.EIonTemp_LP=numpy.zeros(dataOut.NACF-dataOut.cut) dataOut.Phy_LP=numpy.zeros(dataOut.NACF-dataOut.cut) dataOut.EPhy_LP=numpy.zeros(dataOut.NACF-dataOut.cut) dataOut.Phe_LP=numpy.zeros(dataOut.NACF-dataOut.cut) dataOut.EPhe_LP=numpy.zeros(dataOut.NACF-dataOut.cut) dataOut.Density_LP[:]=numpy.copy(dataOut.ne[dataOut.cut:dataOut.NACF]) dataOut.EDensity_LP[:]=numpy.copy(dataOut.ene[dataOut.cut:dataOut.NACF]) dataOut.ElecTemp_LP[:]=numpy.copy(dataOut.te[dataOut.cut:dataOut.NACF]) dataOut.EElecTemp_LP[:]=numpy.copy(dataOut.ete[dataOut.cut:dataOut.NACF]) dataOut.IonTemp_LP[:]=numpy.copy(dataOut.ti[dataOut.cut:dataOut.NACF]) dataOut.EIonTemp_LP[:]=numpy.copy(dataOut.eti[dataOut.cut:dataOut.NACF]) dataOut.Phy_LP[:]=numpy.copy(dataOut.ph[dataOut.cut:dataOut.NACF]) dataOut.EPhy_LP[:]=numpy.copy(dataOut.eph[dataOut.cut:dataOut.NACF]) dataOut.Phe_LP[:]=numpy.copy(dataOut.phe[dataOut.cut:dataOut.NACF]) dataOut.EPhe_LP[:]=numpy.copy(dataOut.ephe[dataOut.cut:dataOut.NACF]) temp_max_lp=6000.0 for i in range(dataOut.NACF-dataOut.cut): if dataOut.ElecTemp_LP[i]<=temp_min or dataOut.ElecTemp_LP[i]>temp_max_lp or dataOut.EElecTemp_LP[i]>temp_max_lp: dataOut.ElecTemp_LP[i]=dataOut.EElecTemp_LP[i]=missing if dataOut.IonTemp_LP[i]<=temp_min or dataOut.IonTemp_LP[i]>temp_max_lp or dataOut.EIonTemp_LP[i]>temp_max_lp: dataOut.IonTemp_LP[i]=dataOut.EIonTemp_LP[i]=missing if dataOut.EPhy_LP[i]<0.0 or dataOut.EPhy_LP[i]>1.0: dataOut.Phy_LP[i]=dataOut.EPhy_LP[i]=missing if dataOut.EPhe_LP[i]<0.0 or dataOut.EPhe_LP[i]>1.0: dataOut.Phe_LP[i]=dataOut.EPhe_LP[i]=missing if dataOut.EDensity_LP[i]>0.0 and dataOut.Density_LP[i]>0.0 and dataOut.Density_LP[i]<9.9e6 and dataOut.EDensity_LP[i]*dataOut.Density_LP[i]<9.9e6: dataOut.EDensity_LP[i]=max(dataOut.EDensity_LP[i],1000.0/dataOut.Density_LP[i]) else: dataOut.Density_LP[i]=missing dataOut.EDensity_LP[i]=1.0 if numpy.isnan(dataOut.Phy_LP[i]): dataOut.EPhy_LP[i]=missing if numpy.isnan(dataOut.Phe_LP[i]): dataOut.EPhe_LP[i]=missing if dataOut.ElecTemp_LP[i]==dataOut.IonTemp_LP[i]: dataOut.EElecTemp_LP[i]=dataOut.EIonTemp_LP[i] if numpy.isnan(dataOut.ElecTemp_LP[i]): dataOut.EElecTemp_LP[i]=missing if numpy.isnan(dataOut.IonTemp_LP[i]): dataOut.EIonTemp_LP[i]=missing if numpy.isnan(dataOut.ElecTemp_LP[i]) or numpy.isnan(dataOut.EElecTemp_LP[i]): dataOut.ElecTemp_LP[i]=dataOut.EElecTemp_LP[i]=dataOut.IonTemp_LP[i]=dataOut.EIonTemp_LP[i]=missing dataOut.DensityFinal=numpy.reshape(numpy.concatenate((dataOut.Density_DP,dataOut.Density_LP)),(1,-1)) dataOut.EDensityFinal=numpy.reshape(numpy.concatenate((dataOut.EDensity_DP,dataOut.EDensity_LP)),(1,-1)) dataOut.ElecTempFinal=numpy.reshape(numpy.concatenate((dataOut.ElecTemp_DP,dataOut.ElecTemp_LP)),(1,-1)) dataOut.EElecTempFinal=numpy.reshape(numpy.concatenate((dataOut.EElecTemp_DP,dataOut.EElecTemp_LP)),(1,-1)) dataOut.IonTempFinal=numpy.reshape(numpy.concatenate((dataOut.IonTemp_DP,dataOut.IonTemp_LP)),(1,-1)) dataOut.EIonTempFinal=numpy.reshape(numpy.concatenate((dataOut.EIonTemp_DP,dataOut.EIonTemp_LP)),(1,-1)) dataOut.PhyFinal=numpy.reshape(numpy.concatenate((dataOut.Phy_DP,dataOut.Phy_LP)),(1,-1)) dataOut.EPhyFinal=numpy.reshape(numpy.concatenate((dataOut.EPhy_DP,dataOut.EPhy_LP)),(1,-1)) dataOut.PheFinal=numpy.reshape(numpy.concatenate((dataOut.Phe_DP,dataOut.Phe_LP)),(1,-1)) dataOut.EPheFinal=numpy.reshape(numpy.concatenate((dataOut.EPhe_DP,dataOut.EPhe_LP)),(1,-1)) nan_array_2=numpy.empty(dataOut.NACF-dataOut.NDP) nan_array_2[:]=numpy.nan dataOut.acfs_DP=numpy.zeros((dataOut.NACF,dataOut.DPL),'float32') dataOut.acfs_error_DP=numpy.zeros((dataOut.NACF,dataOut.DPL),'float32') acfs_dp_aux=dataOut.acfs_to_save.transpose() acfs_error_dp_aux=dataOut.acfs_error_to_save.transpose() for i in range(dataOut.DPL): dataOut.acfs_DP[:,i]=numpy.concatenate((acfs_dp_aux[:,i],nan_array_2)) dataOut.acfs_error_DP[:,i]=numpy.concatenate((acfs_error_dp_aux[:,i],nan_array_2)) dataOut.acfs_DP=dataOut.acfs_DP.transpose() dataOut.acfs_error_DP=dataOut.acfs_error_DP.transpose() dataOut.acfs_LP=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32') dataOut.acfs_error_LP=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32') for i in range(dataOut.NACF): for j in range(dataOut.IBITS): if numpy.abs(dataOut.errors[j,i]/dataOut.output_LP_integrated.real[0,i,0])<1.0: dataOut.acfs_LP[i,j]=dataOut.output_LP_integrated.real[j,i,0]/dataOut.output_LP_integrated.real[0,i,0] dataOut.acfs_LP[i,j]=max(min(dataOut.acfs_LP[i,j],1.0),-1.0) dataOut.acfs_error_LP[i,j]=dataOut.errors[j,i]/dataOut.output_LP_integrated.real[0,i,0] else: dataOut.acfs_LP[i,j]=numpy.nan dataOut.acfs_error_LP[i,j]=numpy.nan dataOut.acfs_LP=dataOut.acfs_LP.transpose() dataOut.acfs_error_LP=dataOut.acfs_error_LP.transpose() return dataOut class ACFs(Operation): def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.aux=1 def run(self,dataOut): if self.aux: self.taup=numpy.zeros(dataOut.DPL,'float32') self.pacf=numpy.zeros(dataOut.DPL,'float32') self.sacf=numpy.zeros(dataOut.DPL,'float32') self.taup_full=numpy.zeros(dataOut.DPL,'float32') self.pacf_full=numpy.zeros(dataOut.DPL,'float32') self.sacf_full=numpy.zeros(dataOut.DPL,'float32') self.x_igcej=numpy.zeros(dataOut.DPL,'float32') self.y_igcej=numpy.zeros(dataOut.DPL,'float32') self.x_ibad=numpy.zeros(dataOut.DPL,'float32') self.y_ibad=numpy.zeros(dataOut.DPL,'float32') self.aux=0 dataOut.acfs_to_plot=numpy.zeros((dataOut.NDP,dataOut.DPL),'float32') dataOut.acfs_to_save=numpy.zeros((dataOut.NDP,dataOut.DPL),'float32') dataOut.acfs_error_to_plot=numpy.zeros((dataOut.NDP,dataOut.DPL),'float32') dataOut.acfs_error_to_save=numpy.zeros((dataOut.NDP,dataOut.DPL),'float32') dataOut.lags_to_plot=numpy.zeros((dataOut.NDP,dataOut.DPL),'float32') dataOut.x_igcej_to_plot=numpy.zeros((dataOut.NDP,dataOut.DPL),'float32') dataOut.x_ibad_to_plot=numpy.zeros((dataOut.NDP,dataOut.DPL),'float32') dataOut.y_igcej_to_plot=numpy.zeros((dataOut.NDP,dataOut.DPL),'float32') dataOut.y_ibad_to_plot=numpy.zeros((dataOut.NDP,dataOut.DPL),'float32') for i in range(dataOut.NSHTS): acfm=dataOut.rhor[i][0]**2+dataOut.rhoi[i][0]**2 if acfm>0: cc=dataOut.rhor[i][0]/acfm ss=dataOut.rhoi[i][0]/acfm else: cc=1. ss=0. # keep only uncontaminated data for l in range(dataOut.DPL): fact=dataOut.DH if (dataOut.igcej[i][l]==0 and dataOut.ibad[i][l]==0): self.pacf_full[l]=min(1.0,max(-1.0,(dataOut.rhor[i][l]*cc + dataOut.rhoi[i][l]*ss)))*fact+dataOut.range1[i] self.sacf_full[l]=min(1.0,numpy.sqrt(dataOut.sd[i][l]))*fact self.taup_full[l]=dataOut.alag[l] self.x_igcej[l]=numpy.nan self.y_igcej[l]=numpy.nan self.x_ibad[l]=numpy.nan self.y_ibad[l]=numpy.nan else: self.pacf_full[l]=numpy.nan self.sacf_full[l]=numpy.nan self.taup_full[l]=numpy.nan if dataOut.igcej[i][l]: self.x_igcej[l]=dataOut.alag[l] self.y_igcej[l]=dataOut.range1[i] self.x_ibad[l]=numpy.nan self.y_ibad[l]=numpy.nan if dataOut.ibad[i][l]: self.x_igcej[l]=numpy.nan self.y_igcej[l]=numpy.nan self.x_ibad[l]=dataOut.alag[l] self.y_ibad[l]=dataOut.range1[i] pacf_new=numpy.copy((self.pacf_full-dataOut.range1[i])/dataOut.DH) sacf_new=numpy.copy(self.sacf_full/dataOut.DH) dataOut.acfs_to_save[i,:]=numpy.copy(pacf_new) dataOut.acfs_error_to_save[i,:]=numpy.copy(sacf_new) dataOut.acfs_to_plot[i,:]=numpy.copy(self.pacf_full) dataOut.acfs_error_to_plot[i,:]=numpy.copy(self.sacf_full) dataOut.lags_to_plot[i,:]=numpy.copy(self.taup_full) dataOut.x_igcej_to_plot[i,:]=numpy.copy(self.x_igcej) dataOut.x_ibad_to_plot[i,:]=numpy.copy(self.x_ibad) dataOut.y_igcej_to_plot[i,:]=numpy.copy(self.y_igcej) dataOut.y_ibad_to_plot[i,:]=numpy.copy(self.y_ibad) missing=numpy.nan#-32767 for i in range(dataOut.NSHTS,dataOut.NDP): for j in range(dataOut.DPL): dataOut.acfs_to_save[i,j]=missing dataOut.acfs_error_to_save[i,j]=missing dataOut.acfs_to_plot[i,j]=missing dataOut.acfs_error_to_plot[i,j]=missing dataOut.lags_to_plot[i,j]=missing dataOut.x_igcej_to_plot[i,j]=missing dataOut.x_ibad_to_plot[i,j]=missing dataOut.y_igcej_to_plot[i,j]=missing dataOut.y_ibad_to_plot[i,j]=missing dataOut.acfs_to_save=dataOut.acfs_to_save.transpose() dataOut.acfs_error_to_save=dataOut.acfs_error_to_save.transpose() return dataOut class CohInt(Operation): isConfig = False __profIndex = 0 __byTime = False __initime = None __lastdatatime = None __integrationtime = None __buffer = None __bufferStride = [] __dataReady = False __profIndexStride = 0 __dataToPutStride = False n = None def __init__(self, **kwargs): Operation.__init__(self, **kwargs) # self.isConfig = False def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False): """ Set the parameters of the integration class. Inputs: n : Number of coherent integrations timeInterval : Time of integration. If the parameter "n" is selected this one does not work overlapping : """ self.__initime = None self.__lastdatatime = 0 self.__buffer = None self.__dataReady = False self.byblock = byblock self.stride = stride if n == None and timeInterval == None: raise ValueError("n or timeInterval should be specified ...") if n != None: self.n = n self.__byTime = False else: self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line self.n = 9999 self.__byTime = True if overlapping: self.__withOverlapping = True self.__buffer = None else: self.__withOverlapping = False self.__buffer = 0 self.__profIndex = 0 def putData(self, data): """ Add a profile to the __buffer and increase in one the __profileIndex """ if not self.__withOverlapping: self.__buffer += data.copy() self.__profIndex += 1 return #Overlapping data nChannels, nHeis = data.shape data = numpy.reshape(data, (1, nChannels, nHeis)) #If the buffer is empty then it takes the data value if self.__buffer is None: self.__buffer = data self.__profIndex += 1 return #If the buffer length is lower than n then stakcing the data value if self.__profIndex < self.n: self.__buffer = numpy.vstack((self.__buffer, data)) self.__profIndex += 1 return #If the buffer length is equal to n then replacing the last buffer value with the data value self.__buffer = numpy.roll(self.__buffer, -1, axis=0) self.__buffer[self.n-1] = data self.__profIndex = self.n return def pushData(self): """ Return the sum of the last profiles and the profiles used in the sum. Affected: self.__profileIndex """ if not self.__withOverlapping: data = self.__buffer n = self.__profIndex self.__buffer = 0 self.__profIndex = 0 return data, n #Integration with Overlapping data = numpy.sum(self.__buffer, axis=0) # print data # raise n = self.__profIndex return data, n def byProfiles(self, data): self.__dataReady = False avgdata = None # n = None # print data # raise self.putData(data) if self.__profIndex == self.n: avgdata, n = self.pushData() self.__dataReady = True return avgdata def byTime(self, data, datatime): self.__dataReady = False avgdata = None n = None self.putData(data) if (datatime - self.__initime) >= self.__integrationtime: avgdata, n = self.pushData() self.n = n self.__dataReady = True return avgdata def integrateByStride(self, data, datatime): # print data if self.__profIndex == 0: self.__buffer = [[data.copy(), datatime]] else: self.__buffer.append([data.copy(),datatime]) self.__profIndex += 1 self.__dataReady = False if self.__profIndex == self.n * self.stride : self.__dataToPutStride = True self.__profIndexStride = 0 self.__profIndex = 0 self.__bufferStride = [] for i in range(self.stride): current = self.__buffer[i::self.stride] data = numpy.sum([t[0] for t in current], axis=0) avgdatatime = numpy.average([t[1] for t in current]) # print data self.__bufferStride.append((data, avgdatatime)) if self.__dataToPutStride: self.__dataReady = True self.__profIndexStride += 1 if self.__profIndexStride == self.stride: self.__dataToPutStride = False # print self.__bufferStride[self.__profIndexStride - 1] # raise return self.__bufferStride[self.__profIndexStride - 1] return None, None def integrate(self, data, datatime=None): if self.__initime == None: self.__initime = datatime if self.__byTime: avgdata = self.byTime(data, datatime) else: avgdata = self.byProfiles(data) self.__lastdatatime = datatime if avgdata is None: return None, None avgdatatime = self.__initime deltatime = datatime - self.__lastdatatime if not self.__withOverlapping: self.__initime = datatime else: self.__initime += deltatime return avgdata, avgdatatime def integrateByBlock(self, dataOut): times = int(dataOut.data.shape[1]/self.n) avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex) id_min = 0 id_max = self.n for i in range(times): junk = dataOut.data[:,id_min:id_max,:] avgdata[:,i,:] = junk.sum(axis=1) id_min += self.n id_max += self.n timeInterval = dataOut.ippSeconds*self.n avgdatatime = (times - 1) * timeInterval + dataOut.utctime self.__dataReady = True return avgdata, avgdatatime def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs): if not self.isConfig: self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs) self.isConfig = True print("inside") if dataOut.flagDataAsBlock: """ Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis] """ avgdata, avgdatatime = self.integrateByBlock(dataOut) dataOut.nProfiles /= self.n else: if stride is None: avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime) else: avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime) # dataOut.timeInterval *= n dataOut.flagNoData = True if self.__dataReady: dataOut.data = avgdata if not dataOut.flagCohInt: dataOut.nCohInt *= self.n dataOut.flagCohInt = True dataOut.utctime = avgdatatime # print avgdata, avgdatatime # raise # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt dataOut.flagNoData = False return dataOut class TimesCode(Operation): """ """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) def run(self,dataOut,code): #code = numpy.repeat(code, repeats=osamp, axis=1) nCodes = numpy.shape(code)[1] #nprofcode = dataOut.nProfiles//nCodes code = numpy.array(code) #print("nHeights",dataOut.nHeights) #print("nheicode",nheicode) #print("Code.Shape",numpy.shape(code)) #print("Code",code[0,:]) nheicode = dataOut.nHeights//nCodes res = dataOut.nHeights%nCodes ''' buffer = numpy.zeros((dataOut.nChannels, nprofcode, nCodes, ndataOut.nHeights), dtype='complex') ''' #exit(1) #for ipr in range(dataOut.nProfiles): #print(dataOut.nHeights) #print(dataOut.data[0,384-2:]) #print(dataOut.profileIndex) #print(dataOut.data[0,:2]) #print(dataOut.data[0,0:64]) #print(dataOut.data[0,64:64+64]) #exit(1) for ich in range(dataOut.nChannels): for ihe in range(nheicode): #print(ihe*nCodes) #print((ihe+1)*nCodes) #dataOut.data[ich,ipr,ihe*nCodes:nCodes*(ihe+1)] #code[ipr,:] #print("before",dataOut.data[ich,ipr,ihe*nCodes:nCodes*(ihe+1)]) #dataOut.data[ich,ipr,ihe*nCodes:nCodes*(ihe+1)] = numpy.prod([dataOut.data[ich,ipr,ihe*nCodes:nCodes*(ihe+1)],code[ipr,:]],axis=0) dataOut.data[ich,ihe*nCodes:nCodes*(ihe+1)] = numpy.prod([dataOut.data[ich,ihe*nCodes:nCodes*(ihe+1)],code[dataOut.profileIndex,:]],axis=0) #print("after",dataOut.data[ich,ipr,ihe*nCodes:nCodes*(ihe+1)]) #exit(1) #print(dataOut.data[0,:2]) #exit(1) #print(nheicode) #print((nheicode)*nCodes) #print(((nheicode)*nCodes)+res) if res != 0: for ich in range(dataOut.nChannels): dataOut.data[ich,nheicode*nCodes:] = numpy.prod([dataOut.data[ich,nheicode*nCodes:],code[dataOut.profileIndex,:res]],axis=0) #pass #print(dataOut.data[0,384-2:]) #exit(1) #dataOut.data = numpy.mean(buffer,axis=1) #print(numpy.shape(dataOut.data)) #print(dataOut.nHeights) #dataOut.heightList = dataOut.heightList[0:nheicode] #print(dataOut.nHeights) #dataOut.nHeights = numpy.shape(dataOut.data)[2] #print(numpy.shape(dataOut.data)) #exit(1) return dataOut ''' class Spectrogram(Operation): """ """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) def run(self,dataOut): import scipy fs = 3200*1e-6 fs = fs/64 fs = 1/fs nperseg=64 noverlap=48 f, t, Sxx = signal.spectrogram(x, fs, return_onesided=False, nperseg=nperseg, noverlap=noverlap, mode='complex') for ich in range(dataOut.nChannels): for ihe in range(nheicode): return dataOut ''' class RemoveDcHae(Operation): def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.DcCounter = 0 def run(self, dataOut): if self.DcCounter == 0: dataOut.DcHae = numpy.zeros((dataOut.data.shape[0],320),dtype='complex') #dataOut.DcHae = [] self.DcCounter = 1 dataOut.dataaux = numpy.copy(dataOut.data) #dataOut.DcHae += dataOut.dataaux[:,1666:1666+320] dataOut.DcHae += dataOut.dataaux[:,0:0+320] hei = 1666 hei = 2000 hei = 1000 hei = 0 #dataOut.DcHae = numpy.concatenate([dataOut.DcHae,dataOut.dataaux[0,hei]],axis = None) return dataOut class SSheightProfiles(Operation): step = None nsamples = None bufferShape = None profileShape = None sshProfiles = None profileIndex = None def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.isConfig = False def setup(self,dataOut ,step = None , nsamples = None): if step == None and nsamples == None: #pass raise ValueError("step or nheights should be specified ...") self.step = step self.nsamples = nsamples self.__nChannels = dataOut.nChannels self.__nProfiles = dataOut.nProfiles self.__nHeis = dataOut.nHeights shape = dataOut.data.shape #nchannels, nprofiles, nsamples ''' print "input nChannels",self.__nChannels print "input nProfiles",self.__nProfiles print "input nHeis",self.__nHeis print "input Shape",shape ''' residue = (shape[1] - self.nsamples) % self.step if residue != 0: print("The residue is %d, step=%d should be multiple of %d to avoid loss of %d samples"%(residue,step,shape[1] - self.nsamples,residue)) deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] numberProfile = self.nsamples numberSamples = (shape[1] - self.nsamples)/self.step ''' print "new numberProfile",numberProfile print "new numberSamples",numberSamples print "New number of profile: %d, number of height: %d, Resolution %f Km"%(numberProfile,numberSamples,deltaHeight*self.step) ''' self.bufferShape = int(shape[0]), int(numberSamples), int(numberProfile) # nchannels, nsamples , nprofiles self.profileShape = int(shape[0]), int(numberProfile), int(numberSamples) # nchannels, nprofiles, nsamples self.buffer = numpy.zeros(self.bufferShape , dtype=numpy.complex) self.sshProfiles = numpy.zeros(self.profileShape, dtype=numpy.complex) def run(self, dataOut, step, nsamples, code = None, repeat = None): #print(dataOut.profileIndex) dataOut.flagNoData = True dataOut.flagDataAsBlock = False profileIndex = None #code = numpy.array(code) #print(dataOut.data[0,:]) #exit(1) if not self.isConfig: #print("STEP",step) self.setup(dataOut, step=step , nsamples=nsamples) self.isConfig = True #print(code[dataOut.profileIndex,:]) #DC_Hae = numpy.array([0.398+0.588j, -0.926+0.306j, -0.536-0.682j, -0.072+0.53j, 0.368-0.356j, 0.996+0.362j]) DC_Hae = numpy.array([ 0.001025 +0.0516375j, 0.03485 +0.20923125j, -0.168 -0.02720625j, -0.1105375 +0.0707125j, -0.20309375-0.09670625j, 0.189775 +0.02716875j])*(-3.5) DC_Hae = numpy.array([ -32.26 +8.66j, -32.26 +8.66j]) DC_Hae = numpy.array([-2.78500000e-01 -1.39175j, -6.63237294e+02+210.4268625j]) #print(dataOut.data[0,13:15]) dataOut.data = dataOut.data - DC_Hae[:,None] #print(dataOut.data[0,13:15]) #exit(1) code = numpy.array(code) roll = 0 code = numpy.roll(code,roll,axis=0) code = numpy.reshape(code,(5,100,64)) block = dataOut.CurrentBlock%5 #print(block) #code_block = code[block-1-2,:,:] day_dif = 1 #day_12 code_block = code[block-1-0,:,:] if repeat is not None: code_block = numpy.repeat(code_block, repeats=repeat, axis=1) #print(dataOut.data[0:2,13]) for i in range(self.buffer.shape[1]): #self.buffer[:,i] = numpy.flip(dataOut.data[:,i*self.step:i*self.step + self.nsamples]) ''' print(dataOut.profileIndex) print(code[dataOut.profileIndex,:]) print("before",dataOut.data[:,i*self.step:i*self.step + self.nsamples]) print("after",dataOut.data[:,i*self.step:i*self.step + self.nsamples]*code[dataOut.profileIndex,:]) exit(1) ''' #dif = numpy.copy(code) if code is not None: ''' code = numpy.array(code) #print(code[0,:]) #print("There is Code") #exit(1) #code = dataOut.code #print(code[0,:]) #exit(1) roll = 0 code = numpy.roll(code,roll,axis=0) code = numpy.reshape(code,(5,100,64)) block = dataOut.CurrentBlock%5 #print(block) #code_block = code[block-1-2,:,:] day_dif = 1 #day_12 code_block = code[block-1-0,:,:] if repeat is not None: code_block = numpy.repeat(code_block, repeats=repeat, axis=1) ''' #code_block = code[0,:,:] #print(code_block[2,:]) #for l in range(dataOut.data.shape[1]): #dataOut.data[:,l] = dataOut.data[:,l] - numpy.array([0.398+0.588j, -0.926+0.306j, -0.536-0.682j, -0.072+0.53j, 0.368-0.356j, 0.996+0.362j]) ##DC_Hae = numpy.array([0.398+0.588j, -0.926+0.306j, -0.536-0.682j, -0.072+0.53j, 0.368-0.356j, 0.996+0.362j]) #print(dataOut.data[0:2,13]) ##dataOut.data = dataOut.data - DC_Hae[:,None] #print(dataOut.data[0:2,13]) #exit(1) #print(dataOut.data[0,i*self.step:i*self.step + self.nsamples]) #print(dataOut.data[1,i*self.step:i*self.step + self.nsamples]) #print(code_block[dataOut.profileIndex,:]) #print(numpy.shape(code_block[dataOut.profileIndex,:])) #exit(1) ###aux = numpy.mean(dataOut.data[:,i*self.step:i*self.step + self.nsamples],axis=1) ###self.buffer[:,i] = (dataOut.data[:,i*self.step:i*self.step + self.nsamples]-aux[:,None])*code_block[dataOut.profileIndex,:] ''' if i == 18: buffer = dataOut.data[0,i*self.step:i*self.step + self.nsamples] import matplotlib.pyplot as plt fig, axes = plt.subplots(figsize=(14, 10)) x = numpy.linspace(0,20,numpy.shape(buffer)[0]) x = numpy.fft.fftfreq(numpy.shape(buffer)[0],0.00005) x = numpy.fft.fftshift(x) plt.plot(x,buffer) plt.show() import time time.sleep(50) ''' #for k in range(dataOut.nChannels): self.buffer[:,i] = dataOut.data[:,i*self.step:i*self.step + self.nsamples]*code_block[dataOut.profileIndex,:] #print(dataOut.data[0,:]) #print(code_block[0,:]) #print(self.buffer[1,i]) #exit(1) else: #print("There is no Code") #exit(1) self.buffer[:,i] = dataOut.data[:,i*self.step:i*self.step + self.nsamples]#*code[dataOut.profileIndex,:] #self.buffer[:,j,self.__nHeis-j*self.step - self.nheights:self.__nHeis-j*self.step] = numpy.flip(dataOut.data[:,j*self.step:j*self.step + self.nheights]) for j in range(self.buffer.shape[0]): self.sshProfiles[j] = numpy.transpose(self.buffer[j]) profileIndex = self.nsamples deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] ippSeconds = (deltaHeight*1.0e-6)/(0.15) #print "ippSeconds",ippSeconds try: if dataOut.concat_m is not None: ippSeconds= ippSeconds/float(dataOut.concat_m) #print "Profile concat %d"%dataOut.concat_m except: pass dataOut.data = self.sshProfiles dataOut.flagNoData = False dataOut.heightList = numpy.arange(self.buffer.shape[1]) *self.step*deltaHeight + dataOut.heightList[0] dataOut.nProfiles = int(dataOut.nProfiles*self.nsamples) ''' print(dataOut.profileIndex) if dataOut.profileIndex == 0: dataOut.data = dataOut.data*1.e5 buffer_prom = ''' #dataOut.utctime = dataOut.utctime - dataOut.profileIndex #print(dataOut.profileIndex) #print(dataOut.data[0,0,0]) ''' if dataOut.profileIndex == 0: self.buffer_prom = numpy.copy(dataOut.data) else: self.buffer_prom = dataOut.data+self.buffer_prom if dataOut.profileIndex == 99: dataOut.data = self.buffer_prom/100 ''' #print(dataOut.data[0,0,0]) #print(dataOut.profileIndex) dataOut.profileIndex = profileIndex dataOut.flagDataAsBlock = True dataOut.ippSeconds = ippSeconds dataOut.step = self.step #print(dataOut.profileIndex) #print(dataOut.heightList) #exit(1) #print(dataOut.times) return dataOut class Decoder(Operation): isConfig = False __profIndex = 0 code = None nCode = None nBaud = None def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.times = None self.osamp = None # self.__setValues = False self.isConfig = False self.setupReq = False def setup(self, code, osamp, dataOut): self.__profIndex = 0 self.code = code self.nCode = len(code) self.nBaud = len(code[0]) if (osamp != None) and (osamp >1): self.osamp = osamp self.code = numpy.repeat(code, repeats=self.osamp, axis=1) self.nBaud = self.nBaud*self.osamp self.__nChannels = dataOut.nChannels self.__nProfiles = dataOut.nProfiles self.__nHeis = dataOut.nHeights if self.__nHeis < self.nBaud: raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud)) #Frequency __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex) __codeBuffer[:,0:self.nBaud] = self.code self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1)) if dataOut.flagDataAsBlock: self.ndatadec = self.__nHeis #- self.nBaud + 1 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex) else: #Time self.ndatadec = self.__nHeis #- self.nBaud + 1 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex) def __convolutionInFreq(self, data): fft_code = self.fft_code[self.__profIndex].reshape(1,-1) fft_data = numpy.fft.fft(data, axis=1) conv = fft_data*fft_code data = numpy.fft.ifft(conv,axis=1) return data def __convolutionInFreqOpt(self, data): raise NotImplementedError def __convolutionInTime(self, data): code = self.code[self.__profIndex] for i in range(self.__nChannels): #aux=numpy.correlate(data[i,:], code, mode='full') #print(numpy.shape(aux)) #print(numpy.shape(data[i,:])) #print(numpy.shape(code)) #exit(1) self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:] return self.datadecTime def __convolutionByBlockInTime(self, data): repetitions = int(self.__nProfiles / self.nCode) junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize)) junk = junk.flatten() code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud)) profilesList = range(self.__nProfiles) #print(numpy.shape(self.datadecTime)) #print(numpy.shape(data)) for i in range(self.__nChannels): for j in profilesList: self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:] return self.datadecTime def __convolutionByBlockInFreq(self, data): raise NotImplementedError("Decoder by frequency fro Blocks not implemented") fft_code = self.fft_code[self.__profIndex].reshape(1,-1) fft_data = numpy.fft.fft(data, axis=2) conv = fft_data*fft_code data = numpy.fft.ifft(conv,axis=2) return data def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None): if dataOut.flagDecodeData: print("This data is already decoded, recoding again ...") if not self.isConfig: if code is None: if dataOut.code is None: raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type) code = dataOut.code else: code = numpy.array(code).reshape(nCode,nBaud) self.setup(code, osamp, dataOut) self.isConfig = True if mode == 3: sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode) if times != None: sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n") if self.code is None: print("Fail decoding: Code is not defined.") return self.__nProfiles = dataOut.nProfiles datadec = None if mode == 3: mode = 0 if dataOut.flagDataAsBlock: """ Decoding when data have been read as block, """ if mode == 0: datadec = self.__convolutionByBlockInTime(dataOut.data) if mode == 1: datadec = self.__convolutionByBlockInFreq(dataOut.data) else: """ Decoding when data have been read profile by profile """ if mode == 0: datadec = self.__convolutionInTime(dataOut.data) if mode == 1: datadec = self.__convolutionInFreq(dataOut.data) if mode == 2: datadec = self.__convolutionInFreqOpt(dataOut.data) if datadec is None: raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode) dataOut.code = self.code dataOut.nCode = self.nCode dataOut.nBaud = self.nBaud dataOut.data = datadec #print("before",dataOut.heightList) dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]] #print("after",dataOut.heightList) dataOut.flagDecodeData = True #asumo q la data esta decodificada if self.__profIndex == self.nCode-1: self.__profIndex = 0 return dataOut self.__profIndex += 1 #print("SHAPE",numpy.shape(dataOut.data)) return dataOut # dataOut.flagDeflipData = True #asumo q la data no esta sin flip class DecoderRoll(Operation): isConfig = False __profIndex = 0 code = None nCode = None nBaud = None def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.times = None self.osamp = None # self.__setValues = False self.isConfig = False self.setupReq = False def setup(self, code, osamp, dataOut): self.__profIndex = 0 self.code = code self.nCode = len(code) self.nBaud = len(code[0]) if (osamp != None) and (osamp >1): self.osamp = osamp self.code = numpy.repeat(code, repeats=self.osamp, axis=1) self.nBaud = self.nBaud*self.osamp self.__nChannels = dataOut.nChannels self.__nProfiles = dataOut.nProfiles self.__nHeis = dataOut.nHeights if self.__nHeis < self.nBaud: raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud)) #Frequency __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex) __codeBuffer[:,0:self.nBaud] = self.code self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1)) if dataOut.flagDataAsBlock: self.ndatadec = self.__nHeis #- self.nBaud + 1 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex) else: #Time self.ndatadec = self.__nHeis #- self.nBaud + 1 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex) def __convolutionInFreq(self, data): fft_code = self.fft_code[self.__profIndex].reshape(1,-1) fft_data = numpy.fft.fft(data, axis=1) conv = fft_data*fft_code data = numpy.fft.ifft(conv,axis=1) return data def __convolutionInFreqOpt(self, data): raise NotImplementedError def __convolutionInTime(self, data): code = self.code[self.__profIndex] #print("code",code[0,0]) for i in range(self.__nChannels): #aux=numpy.correlate(data[i,:], code, mode='full') #print(numpy.shape(aux)) #print(numpy.shape(data[i,:])) #print(numpy.shape(code)) #exit(1) self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:] return self.datadecTime def __convolutionByBlockInTime(self, data): repetitions = int(self.__nProfiles / self.nCode) junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize)) junk = junk.flatten() code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud)) profilesList = range(self.__nProfiles) #print(numpy.shape(self.datadecTime)) #print(numpy.shape(data)) for i in range(self.__nChannels): for j in profilesList: self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:] return self.datadecTime def __convolutionByBlockInFreq(self, data): raise NotImplementedError("Decoder by frequency fro Blocks not implemented") fft_code = self.fft_code[self.__profIndex].reshape(1,-1) fft_data = numpy.fft.fft(data, axis=2) conv = fft_data*fft_code data = numpy.fft.ifft(conv,axis=2) return data def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None): if dataOut.flagDecodeData: print("This data is already decoded, recoding again ...") #print(dataOut.ippSeconds) #exit(1) roll = 0 if self.isConfig: code = numpy.array(code) #roll = 29 code = numpy.roll(code,roll,axis=0) code = numpy.reshape(code,(5,100,64)) block = dataOut.CurrentBlock%5 #code = code[block-1,:,:] #NormalizeDPPower code = code[block-1-1,:,:] #Next Day self.code = numpy.repeat(code, repeats=self.osamp, axis=1) if not self.isConfig: if code is None: if dataOut.code is None: raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type) code = dataOut.code else: code = numpy.array(code) #roll = 29 code = numpy.roll(code,roll,axis=0) code = numpy.reshape(code,(5,100,64)) block = dataOut.CurrentBlock%5 code = code[block-1-1,:,:] #print(code.shape()) #exit(1) code = numpy.array(code).reshape(nCode,nBaud) self.setup(code, osamp, dataOut) self.isConfig = True if mode == 3: sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode) if times != None: sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n") if self.code is None: print("Fail decoding: Code is not defined.") return self.__nProfiles = dataOut.nProfiles datadec = None if mode == 3: mode = 0 if dataOut.flagDataAsBlock: """ Decoding when data have been read as block, """ if mode == 0: datadec = self.__convolutionByBlockInTime(dataOut.data) if mode == 1: datadec = self.__convolutionByBlockInFreq(dataOut.data) else: """ Decoding when data have been read profile by profile """ if mode == 0: datadec = self.__convolutionInTime(dataOut.data) if mode == 1: datadec = self.__convolutionInFreq(dataOut.data) if mode == 2: datadec = self.__convolutionInFreqOpt(dataOut.data) if datadec is None: raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode) dataOut.code = self.code dataOut.nCode = self.nCode dataOut.nBaud = self.nBaud dataOut.data = datadec #print("before",dataOut.heightList) dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]] #print("after",dataOut.heightList) dataOut.flagDecodeData = True #asumo q la data esta decodificada if self.__profIndex == self.nCode-1: self.__profIndex = 0 return dataOut self.__profIndex += 1 #print("SHAPE",numpy.shape(dataOut.data)) return dataOut class ProfileConcat(Operation): isConfig = False buffer = None def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.profileIndex = 0 def reset(self): self.buffer = numpy.zeros_like(self.buffer) self.start_index = 0 self.times = 1 def setup(self, data, m, n=1): self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0])) self.nHeights = data.shape[1]#.nHeights self.start_index = 0 self.times = 1 def concat(self, data): self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy() self.start_index = self.start_index + self.nHeights def run(self, dataOut, m): dataOut.flagNoData = True if not self.isConfig: self.setup(dataOut.data, m, 1) self.isConfig = True if dataOut.flagDataAsBlock: raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False") else: self.concat(dataOut.data) self.times += 1 if self.times > m: dataOut.data = self.buffer self.reset() dataOut.flagNoData = False # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight) dataOut.ippSeconds *= m return dataOut class ProfileSelector(Operation): profileIndex = None # Tamanho total de los perfiles nProfiles = None def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.profileIndex = 0 def incProfileIndex(self): self.profileIndex += 1 if self.profileIndex >= self.nProfiles: self.profileIndex = 0 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex): if profileIndex < minIndex: return False if profileIndex > maxIndex: return False return True def isThisProfileInList(self, profileIndex, profileList): if profileIndex not in profileList: return False return True def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None): """ ProfileSelector: Inputs: profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8) profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30) rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256)) """ if rangeList is not None: if type(rangeList[0]) not in (tuple, list): rangeList = [rangeList] dataOut.flagNoData = True if dataOut.flagDataAsBlock: """ data dimension = [nChannels, nProfiles, nHeis] """ if profileList != None: dataOut.data = dataOut.data[:,profileList,:] if profileRangeList != None: minIndex = profileRangeList[0] maxIndex = profileRangeList[1] profileList = list(range(minIndex, maxIndex+1)) dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:] if rangeList != None: profileList = [] for thisRange in rangeList: minIndex = thisRange[0] maxIndex = thisRange[1] profileList.extend(list(range(minIndex, maxIndex+1))) dataOut.data = dataOut.data[:,profileList,:] dataOut.nProfiles = len(profileList) dataOut.profileIndex = dataOut.nProfiles - 1 dataOut.flagNoData = False return dataOut """ data dimension = [nChannels, nHeis] """ if profileList != None: if self.isThisProfileInList(dataOut.profileIndex, profileList): self.nProfiles = len(profileList) dataOut.nProfiles = self.nProfiles dataOut.profileIndex = self.profileIndex dataOut.flagNoData = False self.incProfileIndex() return dataOut if profileRangeList != None: minIndex = profileRangeList[0] maxIndex = profileRangeList[1] if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex): self.nProfiles = maxIndex - minIndex + 1 dataOut.nProfiles = self.nProfiles dataOut.profileIndex = self.profileIndex dataOut.flagNoData = False self.incProfileIndex() return dataOut if rangeList != None: nProfiles = 0 for thisRange in rangeList: minIndex = thisRange[0] maxIndex = thisRange[1] nProfiles += maxIndex - minIndex + 1 for thisRange in rangeList: minIndex = thisRange[0] maxIndex = thisRange[1] if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex): self.nProfiles = nProfiles dataOut.nProfiles = self.nProfiles dataOut.profileIndex = self.profileIndex dataOut.flagNoData = False self.incProfileIndex() break return dataOut if beam != None: #beam is only for AMISR data if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]): dataOut.flagNoData = False dataOut.profileIndex = self.profileIndex self.incProfileIndex() return dataOut raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter") #return False return dataOut class Reshaper(Operation): def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.__buffer = None self.__nitems = 0 def __appendProfile(self, dataOut, nTxs): if self.__buffer is None: shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) ) self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype) ini = dataOut.nHeights * self.__nitems end = ini + dataOut.nHeights self.__buffer[:, ini:end] = dataOut.data self.__nitems += 1 return int(self.__nitems*nTxs) def __getBuffer(self): if self.__nitems == int(1./self.__nTxs): self.__nitems = 0 return self.__buffer.copy() return None def __checkInputs(self, dataOut, shape, nTxs): if shape is None and nTxs is None: raise ValueError("Reshaper: shape of factor should be defined") if nTxs: if nTxs < 0: raise ValueError("nTxs should be greater than 0") if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0: raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs))) shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs] return shape, nTxs if len(shape) != 2 and len(shape) != 3: raise ValueError("shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights)) if len(shape) == 2: shape_tuple = [dataOut.nChannels] shape_tuple.extend(shape) else: shape_tuple = list(shape) nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles return shape_tuple, nTxs def run(self, dataOut, shape=None, nTxs=None): shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs) dataOut.flagNoData = True profileIndex = None if dataOut.flagDataAsBlock: dataOut.data = numpy.reshape(dataOut.data, shape_tuple) dataOut.flagNoData = False profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1 else: if self.__nTxs < 1: self.__appendProfile(dataOut, self.__nTxs) new_data = self.__getBuffer() if new_data is not None: dataOut.data = new_data dataOut.flagNoData = False profileIndex = dataOut.profileIndex*nTxs else: raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)") deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0] dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs) dataOut.profileIndex = profileIndex dataOut.ippSeconds /= self.__nTxs return dataOut class SplitProfiles(Operation): def __init__(self, **kwargs): Operation.__init__(self, **kwargs) def run(self, dataOut, n): dataOut.flagNoData = True profileIndex = None if dataOut.flagDataAsBlock: #nchannels, nprofiles, nsamples shape = dataOut.data.shape if shape[2] % n != 0: raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2])) new_shape = shape[0], shape[1]*n, int(shape[2]/n) dataOut.data = numpy.reshape(dataOut.data, new_shape) dataOut.flagNoData = False profileIndex = int(dataOut.nProfiles/n) - 1 else: raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)") deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0] dataOut.nProfiles = int(dataOut.nProfiles*n) dataOut.profileIndex = profileIndex dataOut.ippSeconds /= n return dataOut class CombineProfiles(Operation): def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.__remData = None self.__profileIndex = 0 def run(self, dataOut, n): dataOut.flagNoData = True profileIndex = None if dataOut.flagDataAsBlock: #nchannels, nprofiles, nsamples shape = dataOut.data.shape new_shape = shape[0], shape[1]/n, shape[2]*n if shape[1] % n != 0: raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1])) dataOut.data = numpy.reshape(dataOut.data, new_shape) dataOut.flagNoData = False profileIndex = int(dataOut.nProfiles*n) - 1 else: #nchannels, nsamples if self.__remData is None: newData = dataOut.data else: newData = numpy.concatenate((self.__remData, dataOut.data), axis=1) self.__profileIndex += 1 if self.__profileIndex < n: self.__remData = newData #continue return self.__profileIndex = 0 self.__remData = None dataOut.data = newData dataOut.flagNoData = False profileIndex = dataOut.profileIndex/n deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0] dataOut.nProfiles = int(dataOut.nProfiles/n) dataOut.profileIndex = profileIndex dataOut.ippSeconds *= n return dataOut # import collections # from scipy.stats import mode # # class Synchronize(Operation): # # isConfig = False # __profIndex = 0 # # def __init__(self, **kwargs): # # Operation.__init__(self, **kwargs) # # self.isConfig = False # self.__powBuffer = None # self.__startIndex = 0 # self.__pulseFound = False # # def __findTxPulse(self, dataOut, channel=0, pulse_with = None): # # #Read data # # powerdB = dataOut.getPower(channel = channel) # noisedB = dataOut.getNoise(channel = channel)[0] # # self.__powBuffer.extend(powerdB.flatten()) # # dataArray = numpy.array(self.__powBuffer) # # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same") # # maxValue = numpy.nanmax(filteredPower) # # if maxValue < noisedB + 10: # #No se encuentra ningun pulso de transmision # return None # # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0] # # if len(maxValuesIndex) < 2: # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX # return None # # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples # # #Seleccionar solo valores con un espaciamiento de nSamples # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex) # # if len(pulseIndex) < 2: # #Solo se encontro un pulso de transmision con ancho mayor a 1 # return None # # spacing = pulseIndex[1:] - pulseIndex[:-1] # # #remover senales que se distancien menos de 10 unidades o muestras # #(No deberian existir IPP menor a 10 unidades) # # realIndex = numpy.where(spacing > 10 )[0] # # if len(realIndex) < 2: # #Solo se encontro un pulso de transmision con ancho mayor a 1 # return None # # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs) # realPulseIndex = pulseIndex[realIndex] # # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0] # # print "IPP = %d samples" %period # # self.__newNSamples = dataOut.nHeights #int(period) # self.__startIndex = int(realPulseIndex[0]) # # return 1 # # # def setup(self, nSamples, nChannels, buffer_size = 4): # # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float), # maxlen = buffer_size*nSamples) # # bufferList = [] # # for i in range(nChannels): # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN, # maxlen = buffer_size*nSamples) # # bufferList.append(bufferByChannel) # # self.__nSamples = nSamples # self.__nChannels = nChannels # self.__bufferList = bufferList # # def run(self, dataOut, channel = 0): # # if not self.isConfig: # nSamples = dataOut.nHeights # nChannels = dataOut.nChannels # self.setup(nSamples, nChannels) # self.isConfig = True # # #Append new data to internal buffer # for thisChannel in range(self.__nChannels): # bufferByChannel = self.__bufferList[thisChannel] # bufferByChannel.extend(dataOut.data[thisChannel]) # # if self.__pulseFound: # self.__startIndex -= self.__nSamples # # #Finding Tx Pulse # if not self.__pulseFound: # indexFound = self.__findTxPulse(dataOut, channel) # # if indexFound == None: # dataOut.flagNoData = True # return # # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex) # self.__pulseFound = True # self.__startIndex = indexFound # # #If pulse was found ... # for thisChannel in range(self.__nChannels): # bufferByChannel = self.__bufferList[thisChannel] # #print self.__startIndex # x = numpy.array(bufferByChannel) # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples] # # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6 # # dataOut.data = self.__arrayBuffer # # self.__startIndex += self.__newNSamples # # return ##############################LONG PULSE############################## class CrossProdHybrid(CrossProdDP): """Operation to calculate cross products of the Hybrid Experiment. Parameters: ----------- NLAG : int Number of lags for Long Pulse. NRANGE : int Number of samples (heights) for Long Pulse. NCAL : int .* DPL : int Number of lags for Double Pulse. NDN : int .* NDT : int Number of heights for Double Pulse.* NDP : int Number of heights for Double Pulse.* NSCAN : int Number of profiles when the transmitter is on. lagind : intlist .* lagfirst : intlist .* NAVG : int Number of blocks to be "averaged". nkill : int Number of blocks not to be considered when averaging. Example -------- op = proc_unit.addOperation(name='CrossProdHybrid', optype='other') op.addParameter(name='NLAG', value='16', format='int') op.addParameter(name='NRANGE', value='200', format='int') op.addParameter(name='NCAL', value='0', format='int') op.addParameter(name='DPL', value='11', format='int') op.addParameter(name='NDN', value='0', format='int') op.addParameter(name='NDT', value='67', format='int') op.addParameter(name='NDP', value='67', format='int') op.addParameter(name='NSCAN', value='128', format='int') op.addParameter(name='lagind', value='(0,1,2,3,4,5,6,7,0,3,4,5,6,8,9,10)', format='intlist') op.addParameter(name='lagfirst', value='(1,1,1,1,1,1,1,1,0,0,0,0,0,1,1,1)', format='intlist') op.addParameter(name='NAVG', value='16', format='int') op.addParameter(name='nkill', value='6', format='int') """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.bcounter=0 self.aux=1 self.aux_cross_lp=1 self.lag_products_LP_median_estimates_aux=1 def get_products_cabxys_HP(self,dataOut): if self.aux==1: self.set_header_output(dataOut) self.aux=0 self.cax=numpy.zeros((dataOut.NDP,dataOut.DPL,2))# hp:67x11x2 dp: 66x11x2 self.cay=numpy.zeros((dataOut.NDP,dataOut.DPL,2)) self.cbx=numpy.zeros((dataOut.NDP,dataOut.DPL,2)) self.cby=numpy.zeros((dataOut.NDP,dataOut.DPL,2)) self.cax2=numpy.zeros((dataOut.NDP,dataOut.DPL,2)) self.cay2=numpy.zeros((dataOut.NDP,dataOut.DPL,2)) self.cbx2=numpy.zeros((dataOut.NDP,dataOut.DPL,2)) self.cby2=numpy.zeros((dataOut.NDP,dataOut.DPL,2)) self.caxbx=numpy.zeros((dataOut.NDP,dataOut.DPL,2)) self.caxby=numpy.zeros((dataOut.NDP,dataOut.DPL,2)) self.caybx=numpy.zeros((dataOut.NDP,dataOut.DPL,2)) self.cayby=numpy.zeros((dataOut.NDP,dataOut.DPL,2)) self.caxay=numpy.zeros((dataOut.NDP,dataOut.DPL,2)) self.cbxby=numpy.zeros((dataOut.NDP,dataOut.DPL,2)) for i in range(2): # flipped and unflipped for j in range(dataOut.NDP): # loop over true ranges # 67 for k in range(int(dataOut.NSCAN)): # 128 n=dataOut.lagind[k%dataOut.NLAG] # 128=16x8 ax=dataOut.data[0,k,dataOut.NRANGE+dataOut.NCAL+j+i*dataOut.NDT].real-dataOut.dc.real[0] ay=dataOut.data[0,k,dataOut.NRANGE+dataOut.NCAL+j+i*dataOut.NDT].imag-dataOut.dc.imag[0] if dataOut.NRANGE+dataOut.NCAL+j+i*dataOut.NDT+2*n=dataOut.nkill/2 and k=dataOut.nkill/2 and k=dataOut.nkill/2 and k=dataOut.nkill/2 and k ((debris[i-12]+debris[i-11]+debris[i-10]+debris[i-9]+ debris[i+12]+debris[i+11]+debris[i+10]+debris[i+9])/2.0+ thresh)): dataOut.flagNoData=True print("LP Debris detected at",i*15,"km") debris=numpy.zeros(dataOut.NDP,dtype='float32') Range=numpy.arange(0,3000,15) for k in range(2): #flip for i in range(dataOut.NDP): # debris[i]+=numpy.sqrt((dataOut.kaxbx[i,0,k]+dataOut.kayby[i,0,k])**2+(dataOut.kaybx[i,0,k]-dataOut.kaxby[i,0,k])**2) if gmtime(dataOut.utctime).tm_hour > 11: for i in range(2,dataOut.NDP-2): if (debris[i]>3.0*debris[i-2] and debris[i]>3.0*debris[i+2] and Range[i]>200.0 and Range[i]<=540.0): dataOut.flagNoData=True print("DP Debris detected at",i*15,"km") print("inside debris",dataOut.flagNoData) return dataOut class IntegrationHP(IntegrationDP): """Operation to integrate Double Pulse and Long Pulse data. Parameters: ----------- nint : int Number of integrations. Example -------- op = proc_unit.addOperation(name='IntegrationHP', optype='other') op.addParameter(name='nint', value='30', format='int') """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.counter = 0 self.aux = 0 def integration_noise(self,dataOut): if self.counter == 0: dataOut.tnoise=numpy.zeros((dataOut.NR),dtype='float32') dataOut.tnoise+=dataOut.noise_final def integration_for_long_pulse(self,dataOut): if self.counter == 0: dataOut.output_LP_integrated=numpy.zeros((dataOut.NLAG,dataOut.NRANGE,dataOut.NR),order='F',dtype='complex64') dataOut.output_LP_integrated+=dataOut.output_LP def run(self,dataOut,nint=None): dataOut.flagNoData=True #print("flag_inside",dataOut.flagNoData) dataOut.nint=nint dataOut.paramInterval=0#int(dataOut.nint*dataOut.header[7][0]*2 ) dataOut.lat=-11.95 dataOut.lon=-76.87 self.integration_for_long_pulse(dataOut) self.integration_noise(dataOut) if self.counter==dataOut.nint-1: dataOut.tnoise[0]*=0.995 dataOut.tnoise[1]*=0.995 dataOut.pan=dataOut.tnoise[0]/float(dataOut.NSCAN*dataOut.nint*dataOut.NAVG) dataOut.pbn=dataOut.tnoise[1]/float(dataOut.NSCAN*dataOut.nint*dataOut.NAVG) self.integration_for_double_pulse(dataOut) return dataOut class IntegrationLP(Operation): """Operation to integrate Double Pulse and Long Pulse data. Parameters: ----------- nint : int Number of integrations. Example -------- op = proc_unit.addOperation(name='IntegrationHP', optype='other') op.addParameter(name='nint', value='30', format='int') """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.counter = 0 self.aux = 0 def integration_noise(self,dataOut): if self.counter == 0: dataOut.tnoise=numpy.zeros((dataOut.NR),dtype='float32') dataOut.tnoise+=dataOut.noise_final ''' def integration_for_long_pulse(self,dataOut): if self.counter == 0: dataOut.output_LP_integrated=numpy.zeros((dataOut.NLAG,dataOut.NRANGE,dataOut.NR),order='F',dtype='complex64') dataOut.output_LP_integrated+=dataOut.output_LP ''' def integration_for_long_pulse(self,dataOut): #print("inside") #print(self.aux) if self.counter == 0: dataOut.output_LP_integrated=numpy.zeros((dataOut.NLAG,dataOut.NRANGE,dataOut.NR),order='F',dtype='complex64') dataOut.output_LP_integrated+=dataOut.output_LP if self.aux==1: #print("CurrentBlockBBBBB: ",dataOut.CurrentBlock) #print(dataOut.datatime) #dataOut.TimeBlockDate_for_dp_power=dataOut.TimeBlockDate ########dataOut.TimeBlockSeconds_for_dp_power=dataOut.LastAVGDate #print("Date: ",dataOut.TimeBlockDate_for_dp_power) #dataOut.TimeBlockSeconds_for_dp_power=mktime(strptime(dataOut.TimeBlockDate_for_dp_power)) dataOut.TimeBlockSeconds_for_dp_power=dataOut.utctime#dataOut.TimeBlockSeconds-18000 #dataOut.TimeBlockSeconds_for_dp_power=dataOut.LastAVGDate #print("Seconds: ",dataOut.TimeBlockSeconds_for_dp_power) dataOut.bd_time=gmtime(dataOut.TimeBlockSeconds_for_dp_power) #print(dataOut.bd_time) #exit() dataOut.year=dataOut.bd_time.tm_year+(dataOut.bd_time.tm_yday-1)/364.0 dataOut.ut_Faraday=dataOut.bd_time.tm_hour+dataOut.bd_time.tm_min/60.0+dataOut.bd_time.tm_sec/3600.0 #print("date: ", dataOut.TimeBlockDate) self.aux=0 #print("after") self.integration_noise(dataOut) if self.counter==0: self.init_time=dataOut.utctime if self.counter < dataOut.nint: #print("HERE") self.counter+=1 if self.counter==dataOut.nint-1: self.aux=1 #dataOut.TimeBlockDate_for_dp_power=dataOut.TimeBlockDate if self.counter==dataOut.nint: dataOut.flagNoData=False dataOut.utctime=self.init_time self.counter=0 def run(self,dataOut,nint=None): dataOut.flagNoData=True #print("flag_inside",dataOut.flagNoData) dataOut.nint=nint dataOut.paramInterval=0#int(dataOut.nint*dataOut.header[7][0]*2 ) dataOut.lat=-11.95 dataOut.lon=-76.87 self.integration_for_long_pulse(dataOut) if self.counter==dataOut.nint: dataOut.tnoise[0]*=0.995 dataOut.tnoise[1]*=0.995 dataOut.pan=dataOut.tnoise[0]/float(dataOut.NSCAN*dataOut.nint*dataOut.NAVG) dataOut.pbn=dataOut.tnoise[1]/float(dataOut.NSCAN*dataOut.nint*dataOut.NAVG) #self.integration_for_double_pulse(dataOut) print("HERE2") return dataOut class SumFlipsHP(SumFlips): """Operation to sum the flip and unflip part of certain cross products of the Double Pulse. Parameters: ----------- None Example -------- op = proc_unit.addOperation(name='SumFlipsHP', optype='other') """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) def rint2HP(self,dataOut): dataOut.rnint2=numpy.zeros(dataOut.DPL,'float32') for l in range(dataOut.DPL): if(l==0 or (l>=3 and l <=6)): dataOut.rnint2[l]=0.5/float(dataOut.nint*dataOut.NAVG*16.0) else: dataOut.rnint2[l]=0.5/float(dataOut.nint*dataOut.NAVG*8.0) def run(self,dataOut): self.rint2HP(dataOut) self.SumLags(dataOut) return dataOut from schainpy.model.proc import full_profile_profile from scipy.optimize import nnls class LongPulseAnalysis(Operation): """Operation to estimate ACFs, temperatures, total electron density and Hydrogen/Helium fractions from the Long Pulse data. Parameters: ----------- NACF : int .* Example -------- op = proc_unit.addOperation(name='LongPulseAnalysis', optype='other') op.addParameter(name='NACF', value='16', format='int') """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.aux=1 def run(self,dataOut,NACF): dataOut.NACF=NACF dataOut.heightList=dataOut.DH*(numpy.arange(dataOut.NACF)) anoise0=dataOut.tnoise[0] anoise1=anoise0*0.0 #seems to be noise in 1st lag 0.015 before '14 if self.aux: #dataOut.cut=31#26#height=31*15=465 self.cal=numpy.zeros((dataOut.NLAG),'float32') self.drift=numpy.zeros((200),'float32') self.rdrift=numpy.zeros((200),'float32') self.ddrift=numpy.zeros((200),'float32') self.sigma=numpy.zeros((dataOut.NRANGE),order='F',dtype='float32') self.powera=numpy.zeros((dataOut.NRANGE),order='F',dtype='float32') self.powerb=numpy.zeros((dataOut.NRANGE),order='F',dtype='float32') self.perror=numpy.zeros((dataOut.NRANGE),order='F',dtype='float32') dataOut.ene=numpy.zeros((dataOut.NRANGE),'float32') self.dpulse=numpy.zeros((dataOut.NACF),'float32') self.lpulse=numpy.zeros((dataOut.NACF),'float32') dataOut.lags_LP=numpy.zeros((dataOut.IBITS),order='F',dtype='float32') self.lagp=numpy.zeros((dataOut.NACF),'float32') self.u=numpy.zeros((2*dataOut.NACF,2*dataOut.NACF),'float32') dataOut.ne=numpy.zeros((dataOut.NRANGE),order='F',dtype='float32') dataOut.te=numpy.zeros((dataOut.NACF),order='F',dtype='float32') dataOut.ete=numpy.zeros((dataOut.NACF),order='F',dtype='float32') dataOut.ti=numpy.zeros((dataOut.NACF),order='F',dtype='float32') dataOut.eti=numpy.zeros((dataOut.NACF),order='F',dtype='float32') dataOut.ph=numpy.zeros((dataOut.NACF),order='F',dtype='float32') dataOut.eph=numpy.zeros((dataOut.NACF),order='F',dtype='float32') dataOut.phe=numpy.zeros((dataOut.NACF),order='F',dtype='float32') dataOut.ephe=numpy.zeros((dataOut.NACF),order='F',dtype='float32') dataOut.errors=numpy.zeros((dataOut.IBITS,max(dataOut.NRANGE,dataOut.NSHTS)),order='F',dtype='float32') dataOut.fit_array_real=numpy.zeros((max(dataOut.NRANGE,dataOut.NSHTS),dataOut.NLAG),order='F',dtype='float32') dataOut.status=numpy.zeros(1,'float32') dataOut.tx=240.0 #debería provenir del header #hybrid for i in range(dataOut.IBITS): dataOut.lags_LP[i]=float(i)*(dataOut.tx/150.0)/float(dataOut.IBITS) # (float)i*(header.tx/150.0)/(float)IBITS; self.aux=0 dataOut.cut=30 for i in range(30,15,-1): if numpy.nanmax(dataOut.acfs_error_to_plot[i,:])>=10 or dataOut.info2[i]==0: dataOut.cut=i-1 #print(dataOut.cut) #print(dataOut.info2[:]) #print(dataOut.te2[:]) #print(dataOut.ti2[:]) for i in range(dataOut.NLAG): self.cal[i]=sum(dataOut.output_LP_integrated[i,:,3].real) self.cal/=float(dataOut.NRANGE) for j in range(dataOut.NACF+2*dataOut.IBITS+2): dataOut.output_LP_integrated.real[0,j,0]-=anoise0 #lag0 ch0 dataOut.output_LP_integrated.real[1,j,0]-=anoise1 #lag1 ch0 for i in range(1,dataOut.NLAG): #remove cal data from certain lags dataOut.output_LP_integrated.real[i,j,0]-=self.cal[i] k=max(j,26) #constant power below range 26 self.powera[j]=dataOut.output_LP_integrated.real[0,k,0] ## examine drifts here - based on 60 'indep.' estimates nis=dataOut.NSCAN*dataOut.NAVG*dataOut.nint*10 alpha=beta=delta=0.0 nest=0 gamma=3.0/(2.0*numpy.pi*dataOut.lags_LP[1]*1.0e-3) beta=gamma*(math.atan2(dataOut.output_LP_integrated.imag[14,0,2],dataOut.output_LP_integrated.real[14,0,2])-math.atan2(dataOut.output_LP_integrated.imag[1,0,2],dataOut.output_LP_integrated.real[1,0,2]))/13.0 for i in range(1,3): gamma=3.0/(2.0*numpy.pi*dataOut.lags_LP[i]*1.0e-3) for j in range(34,44): rho2=numpy.abs(dataOut.output_LP_integrated[i,j,0])/numpy.abs(dataOut.output_LP_integrated[0,j,0]) dataOut.dphi2=(1.0/rho2-1.0)/(float(2*nis)) dataOut.dphi2*=gamma**2 pest=gamma*math.atan(dataOut.output_LP_integrated.imag[i,j,0]/dataOut.output_LP_integrated.real[i,j,0]) self.drift[nest]=pest self.ddrift[nest]=dataOut.dphi2 self.rdrift[nest]=float(nest) nest+=1 sorted(self.drift[:nest]) for j in range(int(nest/4),int(3*nest/4)): #i=int(self.rdrift[j]) alpha+=self.drift[j]/self.ddrift[j] delta+=1.0/self.ddrift[j] alpha/=delta delta=1./numpy.sqrt(delta) vdrift=alpha-beta dvdrift=delta #need to develop estimate of complete density profile using all #available data #estimate sample variances for long-pulse power profile nis=dataOut.NSCAN*dataOut.NAVG*dataOut.nint self.sigma[:dataOut.NACF+2*dataOut.IBITS+2]=((anoise0+self.powera[:dataOut.NACF+2*dataOut.IBITS+2])**2)/float(nis) ioff=1 #deconvolve rectangular pulse shape from profile ==> powerb, perror ############# START nnlswrap############# if dataOut.ut_Faraday>14.0: alpha_nnlswrap=20.0 else: alpha_nnlswrap=30.0 range1_nnls=dataOut.NACF range2_nnls=dataOut.NACF+dataOut.IBITS-1 g_nnlswrap=numpy.zeros((range1_nnls,range2_nnls),'float32') a_nnlswrap=numpy.zeros((range2_nnls,range2_nnls),'float64') for i in range(range1_nnls): for j in range(range2_nnls): if j>=i and j16: self.dpulse[i]+=dataOut.ph2[k]/dataOut.h2[k] elif k>=36-aux: self.lpulse[i]+=self.powerb[k] self.lagp[i]=self.powera[i] #find scale factor that best merges profiles qi=sum(self.dpulse[32:dataOut.NACF]**2/(self.lagp[32:dataOut.NACF]+anoise0)**2) ri=sum((self.dpulse[32:dataOut.NACF]*self.lpulse[32:dataOut.NACF])/(self.lagp[32:dataOut.NACF]+anoise0)**2) si=sum((self.dpulse[32:dataOut.NACF]*self.lagp[32:dataOut.NACF])/(self.lagp[32:dataOut.NACF]+anoise0)**2) ui=sum(self.lpulse[32:dataOut.NACF]**2/(self.lagp[32:dataOut.NACF]+anoise0)**2) vi=sum((self.lpulse[32:dataOut.NACF]*self.lagp[32:dataOut.NACF])/(self.lagp[32:dataOut.NACF]+anoise0)**2) alpha=(si*ui-vi*ri)/(qi*ui-ri*ri) beta=(qi*vi-ri*si)/(qi*ui-ri*ri) #form density profile estimate, merging rescaled power profiles self.powerb[16:36-aux]=alpha*dataOut.ph2[16:36-aux]/dataOut.h2[16:36-aux] self.powerb[36-aux:dataOut.NACF]*=beta #form Ne estimate, fill in error estimate at low altitudes dataOut.ene[0:36-aux]=dataOut.sdp2[0:36-aux]/dataOut.ph2[0:36-aux] dataOut.ne[:dataOut.NACF]=self.powerb[:dataOut.NACF]*dataOut.h2[:dataOut.NACF]/alpha #now do error propagation: store zero lag error covariance in u nis=dataOut.NSCAN*dataOut.NAVG*dataOut.nint/1 # DLH serious debris removal for i in range(dataOut.NACF): for j in range(i,dataOut.NACF): if j-i>=dataOut.IBITS: self.u[i,j]=0.0 else: self.u[i,j]=dataOut.output_LP_integrated.real[j-i,i,0]**2/float(nis) self.u[i,j]*=(anoise0+dataOut.output_LP_integrated.real[0,i,0])/dataOut.output_LP_integrated.real[0,i,0] self.u[i,j]*=(anoise0+dataOut.output_LP_integrated.real[0,j,0])/dataOut.output_LP_integrated.real[0,j,0] self.u[j,i]=self.u[i,j] #now error analyis for lag product matrix (diag), place in acf_err for i in range(dataOut.NACF): for j in range(dataOut.IBITS): if j==0: dataOut.errors[0,i]=numpy.sqrt(self.u[i,i]) else: dataOut.errors[j,i]=numpy.sqrt(((dataOut.output_LP_integrated.real[0,i,0]+anoise0)*(dataOut.output_LP_integrated.real[0,i+j,0]+anoise0)+dataOut.output_LP_integrated.real[j,i,0]**2)/float(2*nis)) #with suppress_stdout_stderr(): #full_profile_profile.profile(numpy.transpose(dataOut.output_LP_integrated,(2,1,0)),numpy.transpose(dataOut.errors),self.powerb,dataOut.ne,dataOut.lags_LP,dataOut.thb,dataOut.bfm,dataOut.te,dataOut.ete,dataOut.ti,dataOut.eti,dataOut.ph,dataOut.eph,dataOut.phe,dataOut.ephe,dataOut.range1,dataOut.ut,dataOut.NACF,dataOut.fit_array_real,dataOut.status,dataOut.NRANGE,dataOut.IBITS) if dataOut.status>=3.5: dataOut.te[:]=numpy.nan dataOut.ete[:]=numpy.nan dataOut.ti[:]=numpy.nan dataOut.eti[:]=numpy.nan dataOut.ph[:]=numpy.nan dataOut.eph[:]=numpy.nan dataOut.phe[:]=numpy.nan dataOut.ephe[:]=numpy.nan return dataOut class LongPulseAnalysisLP(Operation): """Operation to estimate ACFs, temperatures, total electron density and Hydrogen/Helium fractions from the Long Pulse data. Parameters: ----------- NACF : int .* Example -------- op = proc_unit.addOperation(name='LongPulseAnalysis', optype='other') op.addParameter(name='NACF', value='16', format='int') """ def __init__(self, **kwargs): Operation.__init__(self, **kwargs) self.aux=1 def run(self,dataOut,NACF=None): dataOut.IBITS = 64 dataOut.NACF = dataOut.nHeights# - (2*dataOut.IBITS+5) #print(dataOut.heightList[int(dataOut.NACF)]) #exit(1) #dataOut.heightList=dataOut.DH*(numpy.arange(dataOut.NACF)) anoise0=dataOut.tnoise[0] anoise1=anoise0*0.0 #seems to be noise in 1st lag 0.015 before '14 if self.aux: #dataOut.cut=31#26#height=31*15=465 self.cal=numpy.zeros((dataOut.NLAG),'float32') self.drift=numpy.zeros((200),'float32') self.rdrift=numpy.zeros((200),'float32') self.ddrift=numpy.zeros((200),'float32') self.sigma=numpy.zeros((dataOut.NRANGE),order='F',dtype='float32') self.powera=numpy.zeros((dataOut.NRANGE),order='F',dtype='float32') self.powerb=numpy.zeros((dataOut.NRANGE),order='F',dtype='float32') self.perror=numpy.zeros((dataOut.NRANGE),order='F',dtype='float32') dataOut.ene=numpy.zeros((dataOut.NRANGE),'float32') self.dpulse=numpy.zeros((dataOut.NACF),'float32') self.lpulse=numpy.zeros((dataOut.NACF),'float32') dataOut.lags_LP=numpy.zeros((dataOut.IBITS),order='F',dtype='float32') self.lagp=numpy.zeros((dataOut.NACF),'float32') self.u=numpy.zeros((2*dataOut.NACF,2*dataOut.NACF),'float32') dataOut.ne=numpy.zeros((dataOut.NRANGE),order='F',dtype='float32') dataOut.te=numpy.zeros((dataOut.NACF),order='F',dtype='float32') dataOut.ete=numpy.zeros((dataOut.NACF),order='F',dtype='float32') dataOut.ti=numpy.zeros((dataOut.NACF),order='F',dtype='float32') dataOut.eti=numpy.zeros((dataOut.NACF),order='F',dtype='float32') dataOut.ph=numpy.zeros((dataOut.NACF),order='F',dtype='float32') dataOut.eph=numpy.zeros((dataOut.NACF),order='F',dtype='float32') dataOut.phe=numpy.zeros((dataOut.NACF),order='F',dtype='float32') dataOut.ephe=numpy.zeros((dataOut.NACF),order='F',dtype='float32') dataOut.errors=numpy.zeros((dataOut.IBITS,max(dataOut.NRANGE,1)),order='F',dtype='float32') dataOut.fit_array_real=numpy.zeros((max(dataOut.NRANGE,1),dataOut.NLAG),order='F',dtype='float32') dataOut.status=numpy.zeros(1,'float32') dataOut.tx=480.0 #debería provenir del header #HAE dataOut.h2=numpy.zeros(dataOut.MAXNRANGENDT,'float32') dataOut.range1=numpy.zeros(dataOut.MAXNRANGENDT,order='F',dtype='float32') for i in range(dataOut.IBITS): dataOut.lags_LP[i]=float(i)*(dataOut.tx/150.0)/float(dataOut.IBITS) # (float)i*(header.tx/150.0)/(float)IBITS; self.aux=0 for i in range(dataOut.MAXNRANGENDT): dataOut.range1[i]=dataOut.H0 + i*dataOut.DH dataOut.h2[i]=dataOut.range1[i]**2 dataOut.cut=30 #for i in range(30,15,-1): # if numpy.nanmax(dataOut.acfs_error_to_plot[i,:])>=10 or dataOut.info2[i]==0: # dataOut.cut=i-1 #print(dataOut.cut) #print(dataOut.info2[:]) #print(dataOut.te2[:]) #print(dataOut.ti2[:]) #for i in range(dataOut.NLAG): # self.cal[i]=sum(dataOut.output_LP_integrated[i,:,3].real) #self.cal/=float(dataOut.NRANGE) for j in range(dataOut.NACF):#+2*dataOut.IBITS+2): self.powera[j]=dataOut.output_LP_integrated.real[0,j,0] print(dataOut.heightList[:dataOut.NACF]) import matplotlib.pyplot as plt fig, axes = plt.subplots(figsize=(14, 10)) axes.plot(self.powera[:dataOut.NACF]*dataOut.h2[:dataOut.NACF],dataOut.heightList[:dataOut.NACF]) axes.set_xscale("log", nonposx='clip') #axes.set_xlim(1e18,2e19) axes.set_ylim(180,470) import time plt.title(time.ctime(dataOut.utctime)) plt.show() time.sleep(50) exit(1) ''' for j in range(dataOut.NACF+2*dataOut.IBITS+2): dataOut.output_LP_integrated.real[0,j,0]-=anoise0 #lag0 ch0 dataOut.output_LP_integrated.real[1,j,0]-=anoise1 #lag1 ch0 #for i in range(1,dataOut.NLAG): #remove cal data from certain lags # dataOut.output_LP_integrated.real[i,j,0]-=self.cal[i] k=max(j,26) #constant power below range 26 self.powera[j]=dataOut.output_LP_integrated.real[0,k,0] ## examine drifts here - based on 60 'indep.' estimates nis=dataOut.NSCAN*dataOut.NAVG*dataOut.nint*10 alpha=beta=delta=0.0 nest=0 gamma=3.0/(2.0*numpy.pi*dataOut.lags_LP[1]*1.0e-3) beta=gamma*(math.atan2(dataOut.output_LP_integrated.imag[14,0,2],dataOut.output_LP_integrated.real[14,0,2])-math.atan2(dataOut.output_LP_integrated.imag[1,0,2],dataOut.output_LP_integrated.real[1,0,2]))/13.0 for i in range(1,3): gamma=3.0/(2.0*numpy.pi*dataOut.lags_LP[i]*1.0e-3) for j in range(34,44): rho2=numpy.abs(dataOut.output_LP_integrated[i,j,0])/numpy.abs(dataOut.output_LP_integrated[0,j,0]) dataOut.dphi2=(1.0/rho2-1.0)/(float(2*nis)) dataOut.dphi2*=gamma**2 pest=gamma*math.atan(dataOut.output_LP_integrated.imag[i,j,0]/dataOut.output_LP_integrated.real[i,j,0]) self.drift[nest]=pest self.ddrift[nest]=dataOut.dphi2 self.rdrift[nest]=float(nest) nest+=1 sorted(self.drift[:nest]) for j in range(int(nest/4),int(3*nest/4)): #i=int(self.rdrift[j]) alpha+=self.drift[j]/self.ddrift[j] delta+=1.0/self.ddrift[j] alpha/=delta delta=1./numpy.sqrt(delta) vdrift=alpha-beta dvdrift=delta #need to develop estimate of complete density profile using all #available data #estimate sample variances for long-pulse power profile nis=dataOut.NSCAN*dataOut.NAVG*dataOut.nint self.sigma[:dataOut.NACF+2*dataOut.IBITS+2]=((anoise0+self.powera[:dataOut.NACF+2*dataOut.IBITS+2])**2)/float(nis) ''' ioff=1 #deconvolve rectangular pulse shape from profile ==> powerb, perror ############# START nnlswrap############# if dataOut.ut_Faraday>14.0: alpha_nnlswrap=20.0 else: alpha_nnlswrap=30.0 range1_nnls=dataOut.NACF range2_nnls=dataOut.NACF+dataOut.IBITS-1 g_nnlswrap=numpy.zeros((range1_nnls,range2_nnls),'float32') a_nnlswrap=numpy.zeros((range2_nnls,range2_nnls),'float64') for i in range(range1_nnls): for j in range(range2_nnls): if j>=i and j16: self.dpulse[i]+=dataOut.ph2[k]/dataOut.h2[k] elif k>=36-aux: self.lpulse[i]+=self.powerb[k] self.lagp[i]=self.powera[i] #find scale factor that best merges profiles qi=sum(self.dpulse[32:dataOut.NACF]**2/(self.lagp[32:dataOut.NACF]+anoise0)**2) ri=sum((self.dpulse[32:dataOut.NACF]*self.lpulse[32:dataOut.NACF])/(self.lagp[32:dataOut.NACF]+anoise0)**2) si=sum((self.dpulse[32:dataOut.NACF]*self.lagp[32:dataOut.NACF])/(self.lagp[32:dataOut.NACF]+anoise0)**2) ui=sum(self.lpulse[32:dataOut.NACF]**2/(self.lagp[32:dataOut.NACF]+anoise0)**2) vi=sum((self.lpulse[32:dataOut.NACF]*self.lagp[32:dataOut.NACF])/(self.lagp[32:dataOut.NACF]+anoise0)**2) alpha=(si*ui-vi*ri)/(qi*ui-ri*ri) beta=(qi*vi-ri*si)/(qi*ui-ri*ri) #form density profile estimate, merging rescaled power profiles self.powerb[16:36-aux]=alpha*dataOut.ph2[16:36-aux]/dataOut.h2[16:36-aux] self.powerb[36-aux:dataOut.NACF]*=beta #form Ne estimate, fill in error estimate at low altitudes dataOut.ene[0:36-aux]=dataOut.sdp2[0:36-aux]/dataOut.ph2[0:36-aux] dataOut.ne[:dataOut.NACF]=self.powerb[:dataOut.NACF]*dataOut.h2[:dataOut.NACF]/alpha #now do error propagation: store zero lag error covariance in u nis=dataOut.NSCAN*dataOut.NAVG*dataOut.nint/1 # DLH serious debris removal for i in range(dataOut.NACF): for j in range(i,dataOut.NACF): if j-i>=dataOut.IBITS: self.u[i,j]=0.0 else: self.u[i,j]=dataOut.output_LP_integrated.real[j-i,i,0]**2/float(nis) self.u[i,j]*=(anoise0+dataOut.output_LP_integrated.real[0,i,0])/dataOut.output_LP_integrated.real[0,i,0] self.u[i,j]*=(anoise0+dataOut.output_LP_integrated.real[0,j,0])/dataOut.output_LP_integrated.real[0,j,0] self.u[j,i]=self.u[i,j] #now error analyis for lag product matrix (diag), place in acf_err for i in range(dataOut.NACF): for j in range(dataOut.IBITS): if j==0: dataOut.errors[0,i]=numpy.sqrt(self.u[i,i]) else: dataOut.errors[j,i]=numpy.sqrt(((dataOut.output_LP_integrated.real[0,i,0]+anoise0)*(dataOut.output_LP_integrated.real[0,i+j,0]+anoise0)+dataOut.output_LP_integrated.real[j,i,0]**2)/float(2*nis)) #with suppress_stdout_stderr(): #full_profile_profile.profile(numpy.transpose(dataOut.output_LP_integrated,(2,1,0)),numpy.transpose(dataOut.errors),self.powerb,dataOut.ne,dataOut.lags_LP,dataOut.thb,dataOut.bfm,dataOut.te,dataOut.ete,dataOut.ti,dataOut.eti,dataOut.ph,dataOut.eph,dataOut.phe,dataOut.ephe,dataOut.range1,dataOut.ut,dataOut.NACF,dataOut.fit_array_real,dataOut.status,dataOut.NRANGE,dataOut.IBITS) ''' if dataOut.status>=3.5: dataOut.te[:]=numpy.nan dataOut.ete[:]=numpy.nan dataOut.ti[:]=numpy.nan dataOut.eti[:]=numpy.nan dataOut.ph[:]=numpy.nan dataOut.eph[:]=numpy.nan dataOut.phe[:]=numpy.nan dataOut.ephe[:]=numpy.nan ''' return dataOut class PulsePairVoltage(Operation): ''' Function PulsePair(Signal Power, Velocity) The real component of Lag[0] provides Intensity Information The imag component of Lag[1] Phase provides Velocity Information Configuration Parameters: nPRF = Number of Several PRF theta = Degree Azimuth angel Boundaries Input: self.dataOut lag[N] Affected: self.dataOut.spc ''' isConfig = False __profIndex = 0 __initime = None __lastdatatime = None __buffer = None noise = None __dataReady = False n = None __nch = 0 __nHeis = 0 removeDC = False ipp = None lambda_ = 0 def __init__(self,**kwargs): Operation.__init__(self,**kwargs) def setup(self, dataOut, n = None, removeDC=False): ''' n= Numero de PRF's de entrada ''' self.__initime = None self.__lastdatatime = 0 self.__dataReady = False self.__buffer = 0 self.__profIndex = 0 self.noise = None self.__nch = dataOut.nChannels self.__nHeis = dataOut.nHeights self.removeDC = removeDC self.lambda_ = 3.0e8/(9345.0e6) self.ippSec = dataOut.ippSeconds self.nCohInt = dataOut.nCohInt print("IPPseconds",dataOut.ippSeconds) print("ELVALOR DE n es:", n) if n == None: raise ValueError("n should be specified.") if n != None: if n<2: raise ValueError("n should be greater than 2") self.n = n self.__nProf = n self.__buffer = numpy.zeros((dataOut.nChannels, n, dataOut.nHeights), dtype='complex') def putData(self,data): ''' Add a profile to he __buffer and increase in one the __profiel Index ''' self.__buffer[:,self.__profIndex,:]= data self.__profIndex += 1 return def pushData(self,dataOut): ''' Return the PULSEPAIR and the profiles used in the operation Affected : self.__profileIndex ''' #----------------- Remove DC----------------------------------- if self.removeDC==True: mean = numpy.mean(self.__buffer,1) tmp = mean.reshape(self.__nch,1,self.__nHeis) dc= numpy.tile(tmp,[1,self.__nProf,1]) self.__buffer = self.__buffer - dc #------------------Calculo de Potencia ------------------------ pair0 = self.__buffer*numpy.conj(self.__buffer) pair0 = pair0.real lag_0 = numpy.sum(pair0,1) #------------------Calculo de Ruido x canal-------------------- self.noise = numpy.zeros(self.__nch) for i in range(self.__nch): daux = numpy.sort(pair0[i,:,:],axis= None) self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt) self.noise = self.noise.reshape(self.__nch,1) self.noise = numpy.tile(self.noise,[1,self.__nHeis]) noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis) noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1]) #------------------ Potencia recibida= P , Potencia senal = S , Ruido= N-- #------------------ P= S+N ,P=lag_0/N --------------------------------- #-------------------- Power -------------------------------------------------- data_power = lag_0/(self.n*self.nCohInt) #------------------ Senal --------------------------------------------------- data_intensity = pair0 - noise_buffer data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt) #data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt) for i in range(self.__nch): for j in range(self.__nHeis): if data_intensity[i][j] < 0: data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j])) #----------------- Calculo de Frecuencia y Velocidad doppler-------- pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:]) lag_1 = numpy.sum(pair1,1) data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1) data_velocity = (self.lambda_/2.0)*data_freq #---------------- Potencia promedio estimada de la Senal----------- lag_0 = lag_0/self.n S = lag_0-self.noise #---------------- Frecuencia Doppler promedio --------------------- lag_1 = lag_1/(self.n-1) R1 = numpy.abs(lag_1) #---------------- Calculo del SNR---------------------------------- data_snrPP = S/self.noise for i in range(self.__nch): for j in range(self.__nHeis): if data_snrPP[i][j] < 1.e-20: data_snrPP[i][j] = 1.e-20 #----------------- Calculo del ancho espectral ---------------------- L = S/R1 L = numpy.where(L<0,1,L) L = numpy.log(L) tmp = numpy.sqrt(numpy.absolute(L)) data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L) n = self.__profIndex self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex') self.__profIndex = 0 return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,n def pulsePairbyProfiles(self,dataOut): self.__dataReady = False data_power = None data_intensity = None data_velocity = None data_specwidth = None data_snrPP = None self.putData(data=dataOut.data) if self.__profIndex == self.n: data_power,data_intensity, data_velocity,data_snrPP,data_specwidth, n = self.pushData(dataOut=dataOut) self.__dataReady = True return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth def pulsePairOp(self, dataOut, datatime= None): if self.__initime == None: self.__initime = datatime data_power, data_intensity, data_velocity, data_snrPP, data_specwidth = self.pulsePairbyProfiles(dataOut) self.__lastdatatime = datatime if data_power is None: return None, None, None,None,None,None avgdatatime = self.__initime deltatime = datatime - self.__lastdatatime self.__initime = datatime return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth, avgdatatime def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs): if not self.isConfig: self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs) self.isConfig = True data_power, data_intensity, data_velocity,data_snrPP,data_specwidth, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime) dataOut.flagNoData = True if self.__dataReady: dataOut.nCohInt *= self.n dataOut.dataPP_POW = data_intensity # S dataOut.dataPP_POWER = data_power # P dataOut.dataPP_DOP = data_velocity dataOut.dataPP_SNR = data_snrPP dataOut.dataPP_WIDTH = data_specwidth dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo. dataOut.utctime = avgdatatime dataOut.flagNoData = False return dataOut # import collections # from scipy.stats import mode # # class Synchronize(Operation): # # isConfig = False # __profIndex = 0 # # def __init__(self, **kwargs): # # Operation.__init__(self, **kwargs) # # self.isConfig = False # self.__powBuffer = None # self.__startIndex = 0 # self.__pulseFound = False # # def __findTxPulse(self, dataOut, channel=0, pulse_with = None): # # #Read data # # powerdB = dataOut.getPower(channel = channel) # noisedB = dataOut.getNoise(channel = channel)[0] # # self.__powBuffer.extend(powerdB.flatten()) # # dataArray = numpy.array(self.__powBuffer) # # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same") # # maxValue = numpy.nanmax(filteredPower) # # if maxValue < noisedB + 10: # #No se encuentra ningun pulso de transmision # return None # # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0] # # if len(maxValuesIndex) < 2: # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX # return None # # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples # # #Seleccionar solo valores con un espaciamiento de nSamples # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex) # # if len(pulseIndex) < 2: # #Solo se encontro un pulso de transmision con ancho mayor a 1 # return None # # spacing = pulseIndex[1:] - pulseIndex[:-1] # # #remover senales que se distancien menos de 10 unidades o muestras # #(No deberian existir IPP menor a 10 unidades) # # realIndex = numpy.where(spacing > 10 )[0] # # if len(realIndex) < 2: # #Solo se encontro un pulso de transmision con ancho mayor a 1 # return None # # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs) # realPulseIndex = pulseIndex[realIndex] # # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0] # # print "IPP = %d samples" %period # # self.__newNSamples = dataOut.nHeights #int(period) # self.__startIndex = int(realPulseIndex[0]) # # return 1 # # # def setup(self, nSamples, nChannels, buffer_size = 4): # # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float), # maxlen = buffer_size*nSamples) # # bufferList = [] # # for i in range(nChannels): # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN, # maxlen = buffer_size*nSamples) # # bufferList.append(bufferByChannel) # # self.__nSamples = nSamples # self.__nChannels = nChannels # self.__bufferList = bufferList # # def run(self, dataOut, channel = 0): # # if not self.isConfig: # nSamples = dataOut.nHeights # nChannels = dataOut.nChannels # self.setup(nSamples, nChannels) # self.isConfig = True # # #Append new data to internal buffer # for thisChannel in range(self.__nChannels): # bufferByChannel = self.__bufferList[thisChannel] # bufferByChannel.extend(dataOut.data[thisChannel]) # # if self.__pulseFound: # self.__startIndex -= self.__nSamples # # #Finding Tx Pulse # if not self.__pulseFound: # indexFound = self.__findTxPulse(dataOut, channel) # # if indexFound == None: # dataOut.flagNoData = True # return # # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex) # self.__pulseFound = True # self.__startIndex = indexFound # # #If pulse was found ... # for thisChannel in range(self.__nChannels): # bufferByChannel = self.__bufferList[thisChannel] # #print self.__startIndex # x = numpy.array(bufferByChannel) # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples] # # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6 # # dataOut.data = self.__arrayBuffer # # self.__startIndex += self.__newNSamples # # return