diff --git a/schainpy/model/graphics/jroplot_base.py b/schainpy/model/graphics/jroplot_base.py index 1302d87..a491605 100644 --- a/schainpy/model/graphics/jroplot_base.py +++ b/schainpy/model/graphics/jroplot_base.py @@ -557,7 +557,7 @@ class Plot(Operation): self.sender_time = last_time - attrs = ['titles', 'zmin', 'zmax', 'tag', 'ymin', 'ymax'] + attrs = ['titles', 'zmin', 'zmax', 'tag', 'ymin', 'ymax', 'zlimits'] for attr in attrs: value = getattr(self, attr) if value: @@ -660,11 +660,12 @@ class Plot(Operation): self.poll.register(self.socket, zmq.POLLIN) tm = getattr(dataOut, self.attr_time) - if self.data and 'time' in self.xaxis and (tm - self.tmin) >= self.xrange*60*60: self.save_time = tm self.__plot() - self.tmin += self.xrange*60*60 + #self.tmin += self.xrange*60*60 #Modified by R. Flores + self.tmin += 24*60*60 #Modified by R. Flores + self.data.setup() self.clear_figures() @@ -675,6 +676,7 @@ class Plot(Operation): self.isPlotConfig = True if self.xaxis == 'time': dt = self.getDateTime(tm) + if self.xmin is None: self.tmin = tm self.xmin = dt.hour diff --git a/schainpy/model/graphics/jroplot_spectra.py b/schainpy/model/graphics/jroplot_spectra.py index 38740d2..b995962 100644 --- a/schainpy/model/graphics/jroplot_spectra.py +++ b/schainpy/model/graphics/jroplot_spectra.py @@ -8,7 +8,7 @@ import os import numpy -import collections.abc +#import collections.abc from schainpy.model.graphics.jroplot_base import Plot, plt, log @@ -186,6 +186,7 @@ class SpectraObliquePlot(Plot): ''' data['shift1'] = dataOut.Dop_EEJ_T1[0] data['shift2'] = dataOut.Dop_EEJ_T2[0] + data['max_val_2'] = dataOut.Oblique_params[0,-1,:] data['shift1_error'] = dataOut.Err_Dop_EEJ_T1[0] data['shift2_error'] = dataOut.Err_Dop_EEJ_T2[0] @@ -216,6 +217,7 @@ class SpectraObliquePlot(Plot): shift1 = data['shift1'] #print(shift1) shift2 = data['shift2'] + max_val_2 = data['max_val_2'] err1 = data['shift1_error'] err2 = data['shift2_error'] if ax.firsttime: @@ -238,18 +240,22 @@ class SpectraObliquePlot(Plot): self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^', elinewidth=2.2, marker='o', linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2) self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2) + self.ploterr3 = ax.errorbar(max_val_2, y, xerr=0, fmt='g^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2) + #print("plotter1: ", self.ploterr1,shift1) else: #print("else plotter1: ", self.ploterr1,shift1) self.ploterr1.remove() self.ploterr2.remove() + self.ploterr3.remove() ax.plt.set_array(z[n].T.ravel()) if self.showprofile: ax.plt_profile.set_data(self.data['rti'][n][-1], y) ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y) self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^', elinewidth=2.2, marker='o', linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2) self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2) + self.ploterr3 = ax.errorbar(max_val_2, y, xerr=0, fmt='g^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2) self.titles.append('CH {}: {:3.2f}dB'.format(n, noise)) diff --git a/schainpy/model/graphics/jroplot_voltage_lags.py b/schainpy/model/graphics/jroplot_voltage_lags.py index 1b24ebd..4c9dc99 100644 --- a/schainpy/model/graphics/jroplot_voltage_lags.py +++ b/schainpy/model/graphics/jroplot_voltage_lags.py @@ -4,7 +4,7 @@ import time import math import datetime import numpy -import collections.abc + from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator #YONG from .jroplot_spectra import RTIPlot, NoisePlot @@ -749,9 +749,9 @@ class FaradayAnglePlot(Plot): data['angle'] = numpy.degrees(dataOut.phi) #''' - print(dataOut.phi_uwrp) - print(data['angle']) - exit(1) + #print(dataOut.phi_uwrp) + #print(data['angle']) + #exit(1) #''' data['dphi'] = dataOut.dphi_uc*10 #print(dataOut.dphi) diff --git a/schainpy/model/io/jroIO_base.py b/schainpy/model/io/jroIO_base.py index d0777f4..5ba2b48 100644 --- a/schainpy/model/io/jroIO_base.py +++ b/schainpy/model/io/jroIO_base.py @@ -522,9 +522,8 @@ class Reader(object): def find_files(self, folders, ext, filefmt, startDate=None, endDate=None, expLabel='', last=False): - for path in folders: - files = glob.glob1(path, '*{}'.format(ext)) + files = glob.glob1(path+'/'+expLabel, '*{}'.format(ext)) files.sort() if last: if files: @@ -567,6 +566,7 @@ class Reader(object): if walk: folders = self.find_folders( path, startDate, endDate, folderfmt) + #print("folders: ", folders) else: folders = path.split(',') @@ -928,7 +928,6 @@ class JRODataReader(Reader): self.lastUTTime = self.basicHeaderObj.utc self.flagDiscontinuousBlock = 0 - if deltaTime > self.maxTimeStep: self.flagDiscontinuousBlock = 1 diff --git a/schainpy/model/io/jroIO_madrigal.py b/schainpy/model/io/jroIO_madrigal.py index f8bdbf7..5dd5bff 100644 --- a/schainpy/model/io/jroIO_madrigal.py +++ b/schainpy/model/io/jroIO_madrigal.py @@ -384,6 +384,7 @@ Inputs: __attrs__ = ['path', 'oneDDict', 'ind2DList', 'twoDDict','metadata', 'format', 'blocks'] missing = -32767 + currentDay = None def __init__(self): @@ -608,12 +609,29 @@ Inputs: header.createHeader(**self.header) header.write() + def timeFlag(self): + currentTime = self.dataOut.utctime + timeTuple = time.localtime(currentTime) + dataDay = timeTuple.tm_yday + + if self.currentDay is None: + self.currentDay = dataDay + return False + + #Si el dia es diferente + if dataDay != self.currentDay: + self.currentDay = dataDay + return True + + else: + return False + def putData(self): if self.dataOut.flagNoData: return 0 - if self.dataOut.flagDiscontinuousBlock or self.counter == self.blocks: + if self.dataOut.flagDiscontinuousBlock or self.counter == self.blocks or self.timeFlag(): if self.counter > 0: self.setHeader() self.counter = 0 diff --git a/schainpy/model/io/jroIO_spectra.py b/schainpy/model/io/jroIO_spectra.py index 589c9b6..145b8d4 100644 --- a/schainpy/model/io/jroIO_spectra.py +++ b/schainpy/model/io/jroIO_spectra.py @@ -75,7 +75,7 @@ class SpectraReader(JRODataReader, ProcessingUnit): self.pts2read_SelfSpectra = 0 self.pts2read_CrossSpectra = 0 - self.pts2read_DCchannels = 0 + self.pts2read_DCchannels = 0 self.ext = ".pdata" self.optchar = "P" self.basicHeaderObj = BasicHeader(LOCALTIME) @@ -162,7 +162,7 @@ class SpectraReader(JRODataReader, ProcessingUnit): Exceptions: Si un bloque leido no es un bloque valido """ - + fpointer = self.fp.tell() spc = numpy.fromfile( self.fp, self.dtype[0], self.pts2read_SelfSpectra ) @@ -364,7 +364,7 @@ class SpectraWriter(JRODataWriter, Operation): data.tofile(self.fp) if self.data_cspc is not None: - + cspc = numpy.transpose( self.data_cspc, (0,2,1) ) data = numpy.zeros( numpy.shape(cspc), self.dtype ) #print 'data.shape', self.shape_cspc_Buffer @@ -376,7 +376,7 @@ class SpectraWriter(JRODataWriter, Operation): data.tofile(self.fp) if self.data_dc is not None: - + dc = self.data_dc data = numpy.zeros( numpy.shape(dc), self.dtype ) data['real'] = dc.real @@ -524,4 +524,4 @@ class SpectraWriter(JRODataWriter, Operation): self.processingHeaderObj.processFlags = self.getProcessFlags() - self.setBasicHeader() \ No newline at end of file + self.setBasicHeader() diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py index c6fe251..69a515a 100644 --- a/schainpy/model/proc/jroproc_parameters.py +++ b/schainpy/model/proc/jroproc_parameters.py @@ -1159,11 +1159,12 @@ class Oblique_Gauss_Fit(Operation): freq_max = numpy.max(numpy.abs(freq)) spc_max = numpy.max(spc) - from scipy.signal import medfilt - Nincoh = 20 - Nincoh = 80 + #from scipy.signal import medfilt + #Nincoh = 20 + #Nincoh = 80 Nincoh = Nincoh - spcm = medfilt(spc,11)/numpy.sqrt(Nincoh) + #spcm = medfilt(spc,11)/numpy.sqrt(Nincoh) + spcm = spc/numpy.sqrt(Nincoh) # define a least squares function to optimize def lsq_func(params): @@ -1174,7 +1175,9 @@ class Oblique_Gauss_Fit(Operation): # bounds=([0,-numpy.inf,0,0,-numpy.inf,0,-numpy.inf,0],[numpy.inf,-200,numpy.inf,numpy.inf,0,numpy.inf,0,numpy.inf]) #print(a1,b1,c1,a2,b2,c2,k2,d) #bounds=([0,-numpy.inf,0,-numpy.inf,0,-400,0,0,0],[numpy.inf,-340,numpy.inf,0,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf]) - bounds=([0,-numpy.inf,0,-numpy.inf,0,-400,0,0,0],[numpy.inf,-140,numpy.inf,0,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf]) + #bounds=([0,-numpy.inf,0,-numpy.inf,0,-400,0,0,0],[numpy.inf,-140,numpy.inf,0,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf]) + bounds=([0,-numpy.inf,0,-5,0,-400,0,0,0],[numpy.inf,-300,numpy.inf,5,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf]) + #print(bounds) #bounds=([0,-numpy.inf,0,0,-numpy.inf,0,0,0],[numpy.inf,-200,numpy.inf,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf]) params_scale = [spc_max,freq_max,freq_max,1,spc_max,freq_max,freq_max,1,spc_max] @@ -1278,9 +1281,6 @@ class Oblique_Gauss_Fit(Operation): #print("before return") return A1f, B1f, C1f, A2f, B2f, C2f, Df, error - - - def Double_Gauss_Double_Skew_fit_weight_bound_with_inputs(self, spc, freq, a1, b1, c1, a2, b2, c2, k2, d): from scipy.optimize import least_squares @@ -1469,7 +1469,7 @@ class Oblique_Gauss_Fit(Operation): return A1f, B1f, C1f, Df, error - def run(self, dataOut, mode = 0, Hmin1 = None, Hmax1 = None, Hmin2 = None, Hmax2 = None): + def run(self, dataOut, mode = 0, Hmin1 = None, Hmax1 = None, Hmin2 = None, Hmax2 = None, Dop = 'Shift'): pwcode = 1 @@ -1589,16 +1589,18 @@ class Oblique_Gauss_Fit(Operation): elif mode == 9: #Double Skewed Weighted Bounded no inputs #if numpy.max(spc) <= 0: - if x[numpy.argmax(spc)] <= 0: + from scipy.signal import medfilt + spcm = medfilt(spc,11) + if x[numpy.argmax(spcm)] <= 0: #print("EEJ") - dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_params[0,4,hei],dataOut.Oblique_params[0,5,hei],dataOut.Oblique_params[0,6,hei],dataOut.Oblique_params[0,7,hei],dataOut.Oblique_params[0,8,hei],dataOut.Oblique_params[0,9,hei],dataOut.Oblique_params[0,10,hei],dataOut.Oblique_param_errors[0,:,hei] = self.Double_Gauss_Double_Skew_fit_weight_bound_no_inputs(spc,x,dataOut.nIncohInt) + dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_params[0,4,hei],dataOut.Oblique_params[0,5,hei],dataOut.Oblique_params[0,6,hei],dataOut.Oblique_params[0,7,hei],dataOut.Oblique_params[0,8,hei],dataOut.Oblique_params[0,9,hei],dataOut.Oblique_params[0,10,hei],dataOut.Oblique_param_errors[0,:,hei] = self.Double_Gauss_Double_Skew_fit_weight_bound_no_inputs(spcm,x,dataOut.nIncohInt) #if dataOut.Oblique_params[0,-2,hei] < -500 or dataOut.Oblique_params[0,-2,hei] > 500 or dataOut.Oblique_params[0,-1,hei] < -500 or dataOut.Oblique_params[0,-1,hei] > 500: # dataOut.Oblique_params[0,:,hei] *= numpy.NAN dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,10,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei])) else: #print("CEEJ") - dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_params[0,4,hei],dataOut.Oblique_params[0,5,hei],dataOut.Oblique_params[0,6,hei],dataOut.Oblique_params[0,7,hei],dataOut.Oblique_params[0,8,hei],dataOut.Oblique_params[0,9,hei],dataOut.Oblique_params[0,10,hei],dataOut.Oblique_param_errors[0,:,hei] = self.CEEJ_Skew_fit_weight_bound_no_inputs(spc,x,dataOut.nIncohInt) + dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_params[0,4,hei],dataOut.Oblique_params[0,5,hei],dataOut.Oblique_params[0,6,hei],dataOut.Oblique_params[0,7,hei],dataOut.Oblique_params[0,8,hei],dataOut.Oblique_params[0,9,hei],dataOut.Oblique_params[0,10,hei],dataOut.Oblique_param_errors[0,:,hei] = self.CEEJ_Skew_fit_weight_bound_no_inputs(spcm,x,dataOut.nIncohInt) #if dataOut.Oblique_params[0,-2,hei] < -500 or dataOut.Oblique_params[0,-2,hei] > 500 or dataOut.Oblique_params[0,-1,hei] < -500 or dataOut.Oblique_params[0,-1,hei] > 500: # dataOut.Oblique_params[0,:,hei] *= numpy.NAN dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,10,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei])) @@ -1672,10 +1674,18 @@ class Oblique_Gauss_Fit(Operation): dataOut.lon=-76.87 if mode == 9: #Double Skew Gaussian - dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,-2,:] + #dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,-2,:] #Pos[Max_value] + #dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,1,:] #Shift dataOut.Spec_W_T1 = dataOut.Oblique_params[:,2,:] - dataOut.Dop_EEJ_T2 = dataOut.Oblique_params[:,-1,:] + #dataOut.Dop_EEJ_T2 = dataOut.Oblique_params[:,-1,:] #Pos[Max_value] + #dataOut.Dop_EEJ_T2 = dataOut.Oblique_params[:,5,:] #Shift dataOut.Spec_W_T2 = dataOut.Oblique_params[:,6,:] + if Dop == 'Shift': + dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,1,:] #Shift + dataOut.Dop_EEJ_T2 = dataOut.Oblique_params[:,5,:] #Shift + elif Dop == 'Max': + dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,-2,:] #Pos[Max_value] + dataOut.Dop_EEJ_T2 = dataOut.Oblique_params[:,-1,:] #Pos[Max_value] dataOut.Err_Dop_EEJ_T1 = dataOut.Oblique_param_errors[:,1,:] #En realidad este es el error? dataOut.Err_Spec_W_T1 = dataOut.Oblique_param_errors[:,2,:] @@ -1694,6 +1704,8 @@ class Oblique_Gauss_Fit(Operation): dataOut.Err_Spec_W_T2 = dataOut.Oblique_param_errors[:,5,:] dataOut.mode = mode + dataOut.flagNoData = numpy.all(numpy.isnan(dataOut.Dop_EEJ_T1)) #Si todos los valores son NaN no se prosigue + #dataOut.flagNoData = False #Descomentar solo para ploteo sino mantener comentado return dataOut @@ -6601,8 +6613,9 @@ class IGRFModel(Operation): dataOut.ut=dataOut.bd_time.tm_hour+dataOut.bd_time.tm_min/60.0+dataOut.bd_time.tm_sec/3600.0 self.aux=0 - - dataOut.h=numpy.arange(0.0,15.0*dataOut.MAXNRANGENDT,15.0,dtype='float32') + dh = dataOut.heightList[1]-dataOut.heightList[0] + #dataOut.h=numpy.arange(0.0,15.0*dataOut.MAXNRANGENDT,15.0,dtype='float32') + dataOut.h=numpy.arange(0.0,dh*dataOut.MAXNRANGENDT,dh,dtype='float32') dataOut.bfm=numpy.zeros(dataOut.MAXNRANGENDT,dtype='float32') dataOut.bfm=numpy.array(dataOut.bfm,order='F') dataOut.thb=numpy.zeros(dataOut.MAXNRANGENDT,dtype='float32') @@ -6627,7 +6640,7 @@ class MergeProc(ProcessingUnit): #print(data_inputs) #print("Run: ",self.dataOut.runNextUnit) #exit(1) - #print(numpy.shape([getattr(data, attr_data) for data in data_inputs][1])) + #print("a:", [getattr(data, attr_data) for data in data_inputs][1]) #exit(1) if mode==0: data = numpy.concatenate([getattr(data, attr_data) for data in data_inputs]) @@ -6747,3 +6760,13 @@ class MergeProc(ProcessingUnit): #print(numpy.shape(self.dataOut.data_spc)) #print("*************************GOOD*************************") #exit(1) + + if mode==11: #MST ISR + #data = numpy.concatenate([getattr(data, attr_data) for data in data_inputs],axis=1) + #setattr(self.dataOut, attr_data, data) + setattr(self.dataOut, 'ph2', [getattr(data, attr_data) for data in data_inputs][1]) + print("MST Density", numpy.shape(self.dataOut.ph2)) + print("cf MST: ", self.dataOut.cf) + exit(1) + self.dataOut.ph2 *= self.dataOut.cf + self.dataOut.sdp2 *= self.dataOut.cf diff --git a/schainpy/model/proc/jroproc_spectra_acf.py b/schainpy/model/proc/jroproc_spectra_acf.py index 047890a..ff1edec 100644 --- a/schainpy/model/proc/jroproc_spectra_acf.py +++ b/schainpy/model/proc/jroproc_spectra_acf.py @@ -945,6 +945,7 @@ class SpectraAFCProc(ProcessingUnit): #data = numpy.fft.ifft(data, axis=1, n = 32) #data = numpy.fft.fftshift( data, axes=(1,)) #acf = numpy.abs(data) + #print("data", data[0,0,0]) acf = data[:,:16,:] #acf = data[:,16:,:] #print("SUM: ",numpy.sum(acf)) diff --git a/schainpy/model/proc/jroproc_spectra_lags_faraday.py b/schainpy/model/proc/jroproc_spectra_lags_faraday.py index 52753a2..271b79f 100644 --- a/schainpy/model/proc/jroproc_spectra_lags_faraday.py +++ b/schainpy/model/proc/jroproc_spectra_lags_faraday.py @@ -216,6 +216,7 @@ class SpectraLagProc(ProcessingUnit): self.dataOut.LagPlot=LagPlot #print(self.dataIn.data.shape) + #exit(1) ''' try: print(self.dataIn.data.shape) @@ -243,6 +244,10 @@ class SpectraLagProc(ProcessingUnit): if not self.dataOut.ByLags: #self.dataOut.data = self.dataIn.data + try: + self.dataOut.FlipChannels=self.dataIn.FlipChannels + except: pass + self.dataOut.TimeBlockSeconds=self.dataIn.TimeBlockSeconds self.VoltageType(nFFTPoints,nProfiles,ippFactor,pairsList) else: self.dataOut.nLags = nLags @@ -1810,7 +1815,7 @@ class IntegrationFaradaySpectra2(Operation): ''' #print(self.nHeights) #exit(1) - for l in range(self.nLags):#dataOut.DPL): + for l in range(self.nLags):#dataOut.DPL): #if DP --> nLags=11, elif HP --> nLags=16 #breakFlag=False for k in range(7,self.nHeights): if self.__buffer_cspc is not None: @@ -1818,7 +1823,7 @@ class IntegrationFaradaySpectra2(Operation): outliers_IDs_cspc=[] cspc_outliers_exist=False #indexmin_cspc=0 - for i in range(2): + for i in range(2): #Solo nos interesa los 2 primeros canales que son los canales con señal #for i in range(self.nChannels):#dataOut.nChannels): #if self.TrueLags: #print("HERE") @@ -2871,7 +2876,8 @@ class IntegrationFaradaySpectraNoLags(Operation): self.__lastdatatime = 0 self.__buffer_spc = [] - self.__buffer_cspc = [] + #self.__buffer_cspc = [] + self.__buffer_cspc = None self.__buffer_dc = 0 self.__profIndex = 0 @@ -2949,7 +2955,7 @@ class IntegrationFaradaySpectraNoLags(Operation): return j,sortID #''' - def pushData(self): + def pushData_original_09_11_22(self): """ Return the sum of the last profiles and the profiles used in the sum. @@ -2963,16 +2969,20 @@ class IntegrationFaradaySpectraNoLags(Operation): buffer1=None buffer_cspc=None self.__buffer_spc=numpy.array(self.__buffer_spc) - self.__buffer_cspc=numpy.array(self.__buffer_cspc) + #self.__buffer_cspc=numpy.array(self.__buffer_cspc) + if self.__buffer_cspc is not None: + self.__buffer_cspc=numpy.array(self.__buffer_cspc) freq_dc = int(self.__buffer_spc.shape[2] / 2) #print("FREQ_DC",freq_dc,self.__buffer_spc.shape,self.nHeights) for k in range(7,self.nHeights): - buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k]) - outliers_IDs_cspc=[] - cspc_outliers_exist=False - for i in range(self.nChannels):#dataOut.nChannels): - + if self.__buffer_cspc is not None: + buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k]) + outliers_IDs_cspc=[] + cspc_outliers_exist=False + #for i in range(self.nChannels):#dataOut.nChannels): + for i in range(2):#dataOut.nChannels): buffer1=numpy.copy(self.__buffer_spc[:,i,:,k]) + indexes=[] #sortIDs=[] outliers_IDs=[] @@ -2984,8 +2994,9 @@ class IntegrationFaradaySpectraNoLags(Operation): continue buffer=buffer1[:,j] #if k != 100: - index=int(_HS_algorithm.HS_algorithm(numpy.sort(buffer, axis=None),1)) - sortID = buffer.argsort() + #index=int(_HS_algorithm.HS_algorithm(numpy.sort(buffer, axis=None),1)) + index,sortID=self.hildebrand_sekhon_Integration(buffer,1) + #sortID = buffer.argsort() #else: #index,sortID=self.hildebrand_sekhon_Integration(buffer,1) #if k == 100: @@ -3008,29 +3019,35 @@ class IntegrationFaradaySpectraNoLags(Operation): cspc_outliers_exist=True ###sortdata=numpy.sort(buffer1,axis=0) ###avg2=numpy.mean(sortdata[:indexmin,:],axis=0) - lt=outliers_IDs - avg=numpy.mean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0) + lt=outliers_IDs + #print("buffer1: ", numpy.sum(buffer1)) + #print("outliers: ", buffer1[lt]) + #print("outliers_IDs: ", outliers_IDs) + avg=numpy.nanmean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0) + #print("avg: ", avg) for p in list(outliers_IDs): buffer1[p,:]=avg self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1) ###cspc IDs #indexmin_cspc+=indexmin_cspc - outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs) + if self.__buffer_cspc is not None: + outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs) #if not breakFlag: - outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64')) - if cspc_outliers_exist: - #sortdata=numpy.sort(buffer_cspc,axis=0) - #avg=numpy.mean(sortdata[:indexmin_cpsc,:],axis=0) - lt=outliers_IDs_cspc + if self.__buffer_cspc is not None: + outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64')) + if cspc_outliers_exist: + #sortdata=numpy.sort(buffer_cspc,axis=0) + #avg=numpy.mean(sortdata[:indexmin_cpsc,:],axis=0) + lt=outliers_IDs_cspc - avg=numpy.mean(buffer_cspc[[t for t in range(buffer_cspc.shape[0]) if t not in lt],:],axis=0) - for p in list(outliers_IDs_cspc): - buffer_cspc[p,:]=avg + avg=numpy.mean(buffer_cspc[[t for t in range(buffer_cspc.shape[0]) if t not in lt],:],axis=0) + for p in list(outliers_IDs_cspc): + buffer_cspc[p,:]=avg - self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc) + self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc) #else: #break @@ -3050,7 +3067,15 @@ class IntegrationFaradaySpectraNoLags(Operation): #print(self.__buffer_spc[:,1,3,20,0]) #print(self.__buffer_spc[:,1,5,37,0]) data_spc = numpy.sum(self.__buffer_spc,axis=0) - data_cspc = numpy.sum(self.__buffer_cspc,axis=0) + print("data_spc: ", data_spc[0,:,0]) + print("data_spc: ", data_spc[0,:,7]) + print("shape: ", numpy.shape(data_spc)) + #exit(1) + #data_cspc = numpy.sum(self.__buffer_cspc,axis=0) + if self.__buffer_cspc is not None: + data_cspc = numpy.sum(self.__buffer_cspc,axis=0) + else: + data_cspc = None #print(numpy.shape(data_spc)) #data_spc[1,4,20,0]=numpy.nan @@ -3060,7 +3085,150 @@ class IntegrationFaradaySpectraNoLags(Operation): n = self.__profIndex self.__buffer_spc = [] - self.__buffer_cspc = [] + #self.__buffer_cspc = [] + self.__buffer_cspc = None + self.__buffer_dc = 0 + self.__profIndex = 0 + + return data_spc, data_cspc, data_dc, n + + def pushData(self): + """ + Return the sum of the last profiles and the profiles used in the sum. + + Affected: + + self.__profileIndex + + """ + bufferH=None + buffer=None + buffer1=None + buffer_cspc=None + self.__buffer_spc=numpy.array(self.__buffer_spc) + #self.__buffer_cspc=numpy.array(self.__buffer_cspc) + if self.__buffer_cspc is not None: + self.__buffer_cspc=numpy.array(self.__buffer_cspc) + freq_dc = int(self.__buffer_spc.shape[2] / 2) + #print("FREQ_DC",freq_dc,self.__buffer_spc.shape,self.nHeights) + for k in range(7,self.nHeights): + if self.__buffer_cspc is not None: + buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k]) + outliers_IDs_cspc=[] + cspc_outliers_exist=False + #for i in range(self.nChannels):#dataOut.nChannels): + for i in range(2):#dataOut.nChannels): + buffer1=numpy.copy(self.__buffer_spc[:,i,:,k]) + + indexes=[] + #sortIDs=[] + outliers_IDs=[] + + for j in range(self.nProfiles): + if i==0 and j==freq_dc: #NOT CONSIDERING DC PROFILE AT CHANNEL 0 + continue + if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1 + continue + buffer=buffer1[:,j] + #if k != 100: + index=int(_HS_algorithm.HS_algorithm(numpy.sort(buffer, axis=None),1)) + #index,sortID=self.hildebrand_sekhon_Integration(buffer,1) + sortID = buffer.argsort() + #else: + #index,sortID=self.hildebrand_sekhon_Integration(buffer,1) + #if k == 100: + # print(k,index,sortID) + # exit(1) + #print("index: ", index) + indexes.append(index) + out_IDs = sortID[index:] + avg=numpy.nanmean(buffer1[[t for t in range(buffer1.shape[0]) if t not in out_IDs],:],axis=0) + #print("avg: ", avg) + for p in list(out_IDs): + buffer1[p]=avg + #sortIDs.append(sortID) + outliers_IDs=numpy.append(outliers_IDs,sortID[index:]) + #if k == 100: + # exit(1) + outliers_IDs=numpy.array(outliers_IDs) + outliers_IDs=outliers_IDs.ravel() + outliers_IDs=numpy.unique(outliers_IDs) + outliers_IDs=outliers_IDs.astype(numpy.dtype('int64')) + indexes=numpy.array(indexes) + indexmin=numpy.min(indexes) + ''' + if indexmin != buffer1.shape[0]: + cspc_outliers_exist=True + ###sortdata=numpy.sort(buffer1,axis=0) + ###avg2=numpy.mean(sortdata[:indexmin,:],axis=0) + + lt=outliers_IDs + #print("buffer1: ", numpy.sum(buffer1)) + #print("outliers: ", buffer1[lt]) + #print("outliers_IDs: ", outliers_IDs) + avg=numpy.nanmean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0) + #print("avg: ", avg) + for p in list(outliers_IDs): + buffer1[p,:]=avg + ''' + self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1) + ###cspc IDs + #indexmin_cspc+=indexmin_cspc + if self.__buffer_cspc is not None: + outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs) + + #if not breakFlag: + if self.__buffer_cspc is not None: + outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64')) + if cspc_outliers_exist: + #sortdata=numpy.sort(buffer_cspc,axis=0) + #avg=numpy.mean(sortdata[:indexmin_cpsc,:],axis=0) + lt=outliers_IDs_cspc + + avg=numpy.mean(buffer_cspc[[t for t in range(buffer_cspc.shape[0]) if t not in lt],:],axis=0) + for p in list(outliers_IDs_cspc): + buffer_cspc[p,:]=avg + + self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc) + #else: + #break + + + + + buffer=None + bufferH=None + buffer1=None + buffer_cspc=None + + #print("cpsc",self.__buffer_cspc[:,0,0,0,0]) + #print(self.__profIndex) + #exit() + + buffer=None + #print(self.__buffer_spc[:,1,3,20,0]) + #print(self.__buffer_spc[:,1,5,37,0]) + data_spc = numpy.sum(self.__buffer_spc,axis=0) + print("data_spc: ", data_spc[0,:,0]) + print("data_spc: ", data_spc[0,:,7]) + print("shape: ", numpy.shape(data_spc)) + #exit(1) + #data_cspc = numpy.sum(self.__buffer_cspc,axis=0) + if self.__buffer_cspc is not None: + data_cspc = numpy.sum(self.__buffer_cspc,axis=0) + else: + data_cspc = None + + #print(numpy.shape(data_spc)) + #data_spc[1,4,20,0]=numpy.nan + + #data_cspc = self.__buffer_cspc + data_dc = self.__buffer_dc + n = self.__profIndex + + self.__buffer_spc = [] + #self.__buffer_cspc = [] + self.__buffer_cspc = None self.__buffer_dc = 0 self.__profIndex = 0 @@ -3120,7 +3288,10 @@ class IntegrationFaradaySpectraNoLags(Operation): return dataOut dataOut.flagNoData = True - + #print(numpy.shape(dataOut.data_spc)) + #print(numpy.shape(dataOut.data_cspc)) + #exit(1) + dataOut.data_cspc = None if not self.isConfig: self.setup(dataOut, n, timeInterval, overlapping,DPL ) self.isConfig = True @@ -3751,7 +3922,7 @@ class SnrFaraday(Operation): return dataOut -class SpectraDataToFaraday(Operation): +class SpectraDataToFaraday_07_11_2022(Operation): ''' Written by R. Flores ''' @@ -3990,23 +4161,34 @@ class SpectraDataToFaraday(Operation): data_to_remov_eej = dataOut.dataLag_spc[:,:,:,0] self.normFactor(dataOut) - + #print(dataOut.NDP) dataOut.NDP=dataOut.nHeights dataOut.NR=len(dataOut.channelList) dataOut.DH=dataOut.heightList[1]-dataOut.heightList[0] dataOut.H0=int(dataOut.heightList[0]) self.ConvertData(dataOut) - + #print(dataOut.NDP) + #exit(1) dataOut.NAVG=16#dataOut.rnint2[0] #CHECK THIS! if hasattr(dataOut, 'NRANGE'): dataOut.MAXNRANGENDT = max(dataOut.NRANGE,dataOut.NDT) else: dataOut.MAXNRANGENDT = dataOut.NDP + + #if hasattr(dataOut, 'HP'): + #pass + #else: + self.noise(dataOut) + + ''' if not hasattr(dataOut, 'tnoise'): + print("noise") self.noise(dataOut) - + else: + delattr(dataOut, 'tnoise') + ''' #dataOut.pan = numpy.mean(dataOut.pan) #dataOut.pbn = numpy.mean(dataOut.pbn) #print(dataOut.pan) @@ -4023,7 +4205,197 @@ class SpectraDataToFaraday(Operation): #exit(1) return dataOut +class SpectraDataToFaraday(Operation): + """Operation to use spectra data in Faraday processing. + + Parameters: + ----------- + nint : int + Number of integrations. + + Example + -------- + + op = proc_unit.addOperation(name='SpectraDataToFaraday', optype='other') + + """ + + def __init__(self, **kwargs): + + Operation.__init__(self, **kwargs) + + self.dataLag_spc=None + self.dataLag_cspc=None + self.dataLag_dc=None + + def noise_original(self,dataOut): + + dataOut.noise_lag = numpy.zeros((dataOut.nChannels,dataOut.DPL),'float32') + #print("Lags") + ''' + for lag in range(dataOut.DPL): + #print(lag) + dataOut.data_spc = dataOut.dataLag_spc[:,:,:,lag] + dataOut.noise_lag[:,lag] = dataOut.getNoise(ymin_index=46) + #dataOut.noise_lag[:,lag] = dataOut.getNoise(ymin_index=33,ymax_index=46) + ''' + #print(dataOut.NDP) + #exit(1) + #Channel B + for lag in range(dataOut.DPL): + #print(lag) + dataOut.data_spc = dataOut.dataLag_spc[:,:,:,lag] + max_hei_id = dataOut.NDP - 2*lag + #if lag < 6: + dataOut.noise_lag[1,lag] = dataOut.getNoise(ymin_index=53,ymax_index=max_hei_id)[1] + #else: + #dataOut.noise_lag[1,lag] = numpy.mean(dataOut.noise_lag[1,:6]) + #dataOut.noise_lag[:,lag] = dataOut.getNoise(ymin_index=33,ymax_index=46) + #Channel A + for lag in range(dataOut.DPL): + #print(lag) + dataOut.data_spc = dataOut.dataLag_spc[:,:,:,lag] + dataOut.noise_lag[0,lag] = dataOut.getNoise(ymin_index=53)[0] + nanindex = numpy.argwhere(numpy.isnan(numpy.log10(dataOut.noise_lag[1,:]))) + i1 = nanindex[0][0] + dataOut.noise_lag[1,i1:] = numpy.mean(dataOut.noise_lag[1,:i1]) #El ruido de lags contaminados se + #determina a partir del promedio del + #ruido de los lags limpios + ''' + dataOut.noise_lag[1,:] = dataOut.noise_lag[1,0] #El ruido de los lags diferentes de cero para + #el canal B es contaminado por el Tx y EEJ + #del siguiente perfil, por ello se asigna el ruido + #del lag 0 a todos los lags + ''' + #print("Noise lag: ", 10*numpy.log10(dataOut.noise_lag/dataOut.normFactor)) + #exit(1) + ''' + dataOut.tnoise = dataOut.getNoise(ymin_index=46) + dataOut.tnoise /= float(dataOut.nProfiles*dataOut.nIncohInt) + dataOut.pan = dataOut.tnoise[0] + dataOut.pbn = dataOut.tnoise[1] + ''' + + dataOut.tnoise = dataOut.noise_lag/float(dataOut.nProfiles*dataOut.nIncohInt) + #dataOut.tnoise /= float(dataOut.nProfiles*dataOut.nIncohInt) + dataOut.pan = dataOut.tnoise[0] + dataOut.pbn = dataOut.tnoise[1] + + def noise(self,dataOut): + + dataOut.noise_lag = numpy.zeros((dataOut.nChannels),'float32') + #print("Lags") + ''' + for lag in range(dataOut.DPL): + #print(lag) + dataOut.data_spc = dataOut.dataLag_spc[:,:,:,lag] + dataOut.noise_lag[:,lag] = dataOut.getNoise(ymin_index=46) + #dataOut.noise_lag[:,lag] = dataOut.getNoise(ymin_index=33,ymax_index=46) + ''' + #print(dataOut.NDP) + #exit(1) + #Channel B + + #print(lag) + dataOut.data_spc = dataOut.dataLag_spc[:,:,:,0] + max_hei_id = dataOut.NDP - 2*0 + #if lag < 6: + #dataOut.noise_lag[1] = dataOut.getNoise(ymin_index=80,ymax_index=106)[1] + dataOut.noise_lag[1] = dataOut.getNoise()[1] + #else: + #dataOut.noise_lag[1,lag] = numpy.mean(dataOut.noise_lag[1,:6]) + #dataOut.noise_lag[:,lag] = dataOut.getNoise(ymin_index=33,ymax_index=46) + #Channel A + + #print(lag) + dataOut.data_spc = dataOut.dataLag_spc[:,:,:,0] + dataOut.noise_lag[0] = dataOut.getNoise()[0] + + dataOut.tnoise = dataOut.noise_lag/float(dataOut.nProfiles*dataOut.nIncohInt) + #dataOut.tnoise /= float(dataOut.nProfiles*dataOut.nIncohInt) + dataOut.pan = dataOut.tnoise[0]#*.95 + dataOut.pbn = dataOut.tnoise[1]#*.95 + + def ConvertData(self,dataOut): + + dataOut.TimeBlockSeconds_for_dp_power=dataOut.utctime + dataOut.bd_time=gmtime(dataOut.TimeBlockSeconds_for_dp_power) + 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 + + + tmpx=numpy.zeros((dataOut.nHeights,dataOut.DPL,2),'float32') + tmpx_a2=numpy.zeros((dataOut.nHeights,dataOut.DPL,2),'float32') + tmpx_b2=numpy.zeros((dataOut.nHeights,dataOut.DPL,2),'float32') + tmpx_abr=numpy.zeros((dataOut.nHeights,dataOut.DPL,2),'float32') + tmpx_abi=numpy.zeros((dataOut.nHeights,dataOut.DPL,2),'float32') + dataOut.kabxys_integrated=[tmpx,tmpx,tmpx,tmpx,tmpx_a2,tmpx,tmpx_b2,tmpx,tmpx_abr,tmpx,tmpx_abi,tmpx,tmpx,tmpx] + + dataOut.rnint2=numpy.zeros(dataOut.DPL,'float32') + for l in range(dataOut.DPL): + dataOut.rnint2[l]=1.0/(dataOut.nIncohInt*dataOut.nProfiles)#*dataOut.nProfiles + + + + self.dataLag_spc=(dataOut.dataLag_spc.sum(axis=1))*(dataOut.rnint2[0]/dataOut.nProfiles) + self.dataLag_cspc=(dataOut.dataLag_cspc.sum(axis=1))*(dataOut.rnint2[0]/dataOut.nProfiles) + #self.dataLag_dc=dataOut.dataLag_dc.sum(axis=1)/dataOut.rnint2[0] + + + dataOut.kabxys_integrated[4][:,:,0]=self.dataLag_spc[0,:,:].real + #dataOut.kabxys_integrated[5][:,:,0]+=self.dataLag_spc[0,:,:].imag + dataOut.kabxys_integrated[6][:,:,0]=self.dataLag_spc[1,:,:].real + #dataOut.kabxys_integrated[7][:,:,0]+=self.dataLag_spc[1,:,:].imag + + dataOut.kabxys_integrated[8][:,:,0]=self.dataLag_cspc[0,:,:].real + dataOut.kabxys_integrated[10][:,:,0]=self.dataLag_cspc[0,:,:].imag + + ''' + print(dataOut.kabxys_integrated[4][:,0,0]) + print(dataOut.kabxys_integrated[6][:,0,0]) + print("times 12") + print(dataOut.kabxys_integrated[4][:,0,0]*dataOut.nProfiles) + print(dataOut.kabxys_integrated[6][:,0,0]*dataOut.nProfiles) + print(dataOut.rnint2[0]) + input() + ''' + + + + + + + def run(self,dataOut): + + dataOut.paramInterval=0#int(dataOut.nint*dataOut.header[7][0]*2 ) + dataOut.lat=-11.95 + dataOut.lon=-76.87 + + dataOut.NDP=dataOut.nHeights + dataOut.NR=len(dataOut.channelList) + dataOut.DH=dataOut.heightList[1]-dataOut.heightList[0] + dataOut.H0=int(dataOut.heightList[0]) + #print(dataOut.data_spc.shape) + dataOut.dataLag_spc = numpy.stack((dataOut.data_spc, dataOut.data_spc), axis=-1) + dataOut.dataLag_cspc = numpy.stack((dataOut.data_cspc, dataOut.data_cspc), axis=-1) + #print(dataOut.dataLag_spc.shape) + dataOut.DPL = numpy.shape(dataOut.dataLag_spc)[-1] + #exit(1) + self.ConvertData(dataOut) + self.noise(dataOut) + dataOut.NAVG=16#dataOut.rnint2[0] #CHECK THIS! + dataOut.MAXNRANGENDT=dataOut.NDP + ''' + print(dataOut.kabxys_integrated[4][:,0,0]) + import matplotlib.pyplot as plt + plt.plot(dataOut.kabxys_integrated[4][:,0,0],dataOut.heightList) + plt.axvline(dataOut.pan) + plt.xlim(0,1.e3) + plt.show() + ''' + dataOut.DPL = 1 + return dataOut class SpectraDataToHybrid(SpectraDataToFaraday): ''' @@ -4352,7 +4724,7 @@ class SpectraDataToHybrid_V2(SpectraDataToFaraday): dataOut.pbn_LP=dataOut.tnoise[1]/float(dataOut.nProfiles_LP*dataOut.nIncohInt_LP) def ConvertDataLP(self,dataOut): - + #print(numpy.shape(dataOut.data_acf)) #print(dataOut.dataLag_spc[:,:,:,1]/dataOut.data_spc) #exit(1) self.normfactor_LP=1.0/(dataOut.nIncohInt_LP*dataOut.nProfiles_LP)#*dataOut.nProfiles @@ -4360,10 +4732,10 @@ class SpectraDataToHybrid_V2(SpectraDataToFaraday): #print("Power: ",numpy.mean(dataOut.dataLag_spc_LP[0,:,100])) #buffer = dataOut.data_acf*(1./(normfactor*dataOut.nProfiles_LP)) #buffer = dataOut.data_acf*(1./(normfactor)) - buffer = dataOut.data_acf#*(self.normfactor_LP) + buffer = dataOut.data_acf#*(self.normfactor_LP) #nChannels x nProfiles (nLags) x nHeights #print("acf: ",numpy.sum(buffer)) - dataOut.output_LP_integrated = numpy.transpose(buffer,(1,2,0)) + dataOut.output_LP_integrated = numpy.transpose(buffer,(1,2,0)) #nProfiles (nLags) x nHeights x nChannels def normFactor(self,dataOut): dataOut.rnint2=numpy.zeros(dataOut.DPL,'float32') @@ -4394,7 +4766,9 @@ class SpectraDataToHybrid_V2(SpectraDataToFaraday): #dataOut.output_LP_integrated[:,:,3] *= float(dataOut.NSCAN/22)#(dataOut.nNoiseProfiles) #Corrects the zero padding dataOut.nis=dataOut.NSCAN*dataOut.nIncohInt_LP*10 + #print("nis/10: ", dataOut.NSCAN,dataOut.nIncohInt_LP,dataOut.nProfiles_LP) dataOut.nis=dataOut.NSCAN*dataOut.nIncohInt_LP*dataOut.nProfiles_LP*10 + dataOut.nis=dataOut.nIncohInt_LP*dataOut.nProfiles_LP*10 #Removemos NSCAN debido a que está incluido en nProfiles_LP self.ConvertData(dataOut) diff --git a/schainpy/model/proc/jroproc_voltage.py b/schainpy/model/proc/jroproc_voltage.py index d0b8777..f6531bb 100644 --- a/schainpy/model/proc/jroproc_voltage.py +++ b/schainpy/model/proc/jroproc_voltage.py @@ -2543,12 +2543,11 @@ class CleanCohEchoes(Operation): dataOut.flagSpreadF = True #Removing echoes greater than 35 dB - if isinstance(dataOut.pbn, collections.abc.Sequence): - maxdB = 10*numpy.log10(dataOut.pbn[0]) + 10 #Lag 0 NOise + if hasattr(dataOut.pbn, "__len__"): + maxdB = 10*numpy.log10(dataOut.pbn[0]) + 10 #Lag 0 Noise else: - maxdB = 10*numpy.log10(dataOut.pbn) + 10 #Lag 0 NOise - #maxdB = 35 #DEBERÍA SER NOISE+ALGO!!!!!!!!!!!!!!!!!!!!!! - #print("noise: ",maxdB - 10) + maxdB = 10*numpy.log10(dataOut.pbn) + 10 + #print(dataOut.kabxys_integrated[6][:,0,0]) data = numpy.copy(10*numpy.log10(dataOut.kabxys_integrated[6][:,0,0])) #Lag0 ChB #print("data: ",data) @@ -3052,8 +3051,8 @@ class NoisePower(Operation): 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 + dataOut.pan=.8*dataOut.pnoise[0] # weights could change + dataOut.pbn=.8*dataOut.pnoise[1] # weights could change ''' print("pan: ",dataOut.pan) print("pbn: ",dataOut.pbn) @@ -3107,11 +3106,12 @@ class DoublePulseACFs(Operation): panrm=numpy.zeros((dataOut.NDP,dataOut.DPL), dtype=float) id = numpy.where(dataOut.heightList>700)[0] - + #print("kabxys: ", numpy.shape(dataOut.kabxys_integrated)) 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]) + #print("pa::",pa) pb=numpy.abs(dataOut.kabxys_integrated[6][i,j,0]+dataOut.kabxys_integrated[7][i,j,0]) st4=pa*pb @@ -3200,17 +3200,19 @@ class DoublePulseACFs(Operation): exit(1) ''' #print(pa) + #print("pa: ", numpy.shape(pa)) + #print(numpy.shape(dataOut.heightList)) ''' import matplotlib.pyplot as plt - #plt.plot(dataOut.p[:,-1],dataOut.heightList) - plt.plot(pa/dataOut.pan-1.,dataOut.heightList) - plt.plot(pb/dataOut.pbn-1.,dataOut.heightList) + plt.plot(dataOut.p[:,-1],dataOut.heightList) + #plt.plot(pa/dataOut.pan-1.,dataOut.heightList) + #plt.plot(pb/dataOut.pbn-1.,dataOut.heightList) plt.grid() - #plt.xlim(0,1e5) + plt.xlim(0,1e5) plt.show() #print("p: ",dataOut.p[33,:]) #exit(1) - ''' + #''' #print(numpy.sum(dataOut.rhor)) #exit(1) return dataOut @@ -3346,7 +3348,7 @@ class DoublePulseACFs_PerLag(Operation): exit(1) ''' #print(pa) - ''' + #''' import matplotlib.pyplot as plt #plt.plot(dataOut.p[:,-1],dataOut.heightList) plt.plot(pa/dataOut.pan-1.,dataOut.heightList) @@ -3356,7 +3358,7 @@ class DoublePulseACFs_PerLag(Operation): plt.show() #print("p: ",dataOut.p[33,:]) #exit(1) - ''' + #''' return dataOut class FaradayAngleAndDPPower(Operation): @@ -3446,6 +3448,7 @@ class FaradayAngleAndDPPower(Operation): dataOut.flagTeTiCorrection = False #print(dataOut.ph2) + #exit(1) return dataOut @@ -3519,10 +3522,11 @@ class ElectronDensityFaraday(Operation): #exit(1) ''' import matplotlib.pyplot as plt - plt.plot(dataOut.bki) + plt.plot(dataOut.phi,dataOut.heightList) plt.show() ''' - + #print(dataOut.bki) + print(dataOut.NDP) 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 @@ -3535,8 +3539,15 @@ class ElectronDensityFaraday(Operation): 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 + ''' #print(dataOut.dphi) #exit(1) + import matplotlib.pyplot as plt + plt.plot(dataOut.dphi,dataOut.heightList) + plt.grid() + plt.xlim(0,1e7) + plt.show() + ''' return dataOut @@ -3962,6 +3973,7 @@ class NormalizeDPPowerRoberto_V2(Operation): #print(dataOut.ph2) #input() # in case of spread F, normalize much higher + #print("dens: ", dataOut.dphi,dataOut.ph2) if(dataOut.cf11.0 and dataOut.range1[i]>150.0 and dataOut.range1[i]<400.0: if dataOut.ut_Faraday>11.0 and dataOut.range1[i]>150.0 and dataOut.range1[i]<300.0: @@ -4419,7 +4446,7 @@ class DenCorrection(NormalizeDPPowerRoberto_V2): ''' import matplotlib.pyplot as plt plt.clf() - plt.plot(aux,dataOut.heightList[:dataOut.NSHTS],'*:',label='Fitting') + #plt.plot(aux,dataOut.heightList[:dataOut.NSHTS],'*:',label='Fitting') plt.plot(my_aux,dataOut.heightList[:dataOut.NSHTS],'*:',label='Ratio') #plt.plot(acf_Temps,dataOut.heightList[:dataOut.NSHTS],'b*:',label='Temps') #plt.plot(acf_no_Temps,dataOut.heightList[:dataOut.NSHTS],'k*:',label='No Temps') @@ -4431,7 +4458,9 @@ class DenCorrection(NormalizeDPPowerRoberto_V2): #plt.xlim(.99,1.25) #plt.show() #plt.savefig("/home/roberto/Pictures/Density_Comparison/FactorEf_NoLimits/{}.png".format(dataOut.utctime)) - plt.savefig("/home/roberto/Pictures/Faraday/2022/08/Density_Comparison/FactorEf/{}.png".format(dataOut.utctime)) + #plt.savefig("/home/roberto/Pictures/Density_Comparison/V2/TeTi_cte/Te2/{}.png".format(dataOut.utctime)) + #plt.savefig("/home/roberto/Pictures/Density_Comparison/V2/bline/Te2/{}.png".format(dataOut.utctime)) + plt.savefig("/home/roberto/Pictures/Density_Comparison/V3/{}.png".format(dataOut.utctime)) #plt.savefig("/home/roberto/Pictures/Faraday_TeTi_Test/Ratio/{}.png".format(dataOut.utctime)) ''' #print("inside correction",dataOut.ph2) @@ -4775,10 +4804,10 @@ class DataSaveCleaner(Operation): #print(time_text.hour,time_text.minute) #if (time_text.hour == 16 and time_text.minute==48) or (time_text.hour == 19 and time_text.minute ==49 ) or (time_text.hour >= 0 and time_text.hour < 5): #Year: 2022, DOY:241 - #if (time_text.hour == 5 and time_text.minute==21) or (time_text.hour == 19 and time_text.minute ==49 ) or (time_text.hour == 7 and time_text.minute==40) or (time_text.hour == 7 and time_text.minute==50) or (time_text.hour >= 8 and time_text.hour < 11) or (time_text.hour == 11 and time_text.minute==2) or (time_text.hour == 11 and time_text.minute==13): #Year: 2022, DOY:242 + if (time_text.hour == 5 and time_text.minute==21) or (time_text.hour == 19 and time_text.minute ==49 ) or (time_text.hour == 7 and time_text.minute==40) or (time_text.hour == 7 and time_text.minute==50) or (time_text.hour >= 8 and time_text.hour < 11) or (time_text.hour == 11 and time_text.minute==2) or (time_text.hour == 11 and time_text.minute==13): #Year: 2022, DOY:242 #if (time_text.hour >= 8 and time_text.hour < 11) or (time_text.hour == 11 and time_text.minute==2) or (time_text.hour == 11 and time_text.minute==13) or (time_text.hour == 11 and time_text.minute==24): #Year: 2022, DOY:243 #if (time_text.hour >= 9 and time_text.hour < 11) or (time_text.hour == 8 and time_text.minute==12) or (time_text.hour == 8 and time_text.minute==22) or (time_text.hour == 8 and time_text.minute==33) or (time_text.hour == 8 and time_text.minute==44) or (time_text.hour == 8 and time_text.minute==54) or (time_text.hour == 11 and time_text.minute==2) or (time_text.hour == 11 and time_text.minute==13): #Year: 2022, DOY:245 - if (time_text.hour >= 8 and time_text.hour < 11) or (time_text.hour == 1) or (time_text.hour == 0 and time_text.minute==25) or (time_text.hour == 0 and time_text.minute==36) or (time_text.hour == 0 and time_text.minute==47) or (time_text.hour == 0 and time_text.minute==57) or (time_text.hour == 2 and time_text.minute==1) or (time_text.hour == 11 and time_text.minute==2) or (time_text.hour == 11 and time_text.minute==13) or (time_text.hour == 11 and time_text.minute==24) or (time_text.hour == 7 and time_text.minute==40) or (time_text.hour == 7 and time_text.minute==50) or (time_text.hour == 3 and time_text.minute==5): #Year: 2022, DOY:244 + #if (time_text.hour >= 8 and time_text.hour < 11) or (time_text.hour == 1) or (time_text.hour == 0 and time_text.minute==15) or (time_text.hour == 0 and time_text.minute==25) or (time_text.hour == 0 and time_text.minute==36) or (time_text.hour == 0 and time_text.minute==47) or (time_text.hour == 0 and time_text.minute==57) or (time_text.hour == 2 and time_text.minute==1) or (time_text.hour == 11 and time_text.minute==2) or (time_text.hour == 11 and time_text.minute==13) or (time_text.hour == 11 and time_text.minute==24) or (time_text.hour == 7 and time_text.minute==40) or (time_text.hour == 7 and time_text.minute==50) or (time_text.hour == 3 and time_text.minute==5) or (time_text.hour == 3 and time_text.minute==16) or (time_text.hour == 3 and time_text.minute==27): #Year: 2022, DOY:244 dataOut.DensityFinal[0,:]=missing dataOut.EDensityFinal[0,:]=missing @@ -4813,7 +4842,17 @@ class DataSaveCleaner(Operation): dataOut.EIonTempFinal[0,id_aux:]=missing dataOut.PhyFinal[0,id_aux:]=missing dataOut.EPhyFinal[0,id_aux:]=missing - if (time_text.hour == 20 and time_text.minute == 29) or (time_text.hour == 20 and time_text.minute == 44): #Year: 2022, DOY:243 + if (time_text.hour == 20 and time_text.minute == 29): #Year: 2022, DOY:243 + id_aux = 30 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + if (time_text.hour == 20 and time_text.minute == 44): #Year: 2022, DOY:243 id_aux = 31 dataOut.DensityFinal[0,id_aux:]=missing dataOut.EDensityFinal[0,id_aux:]=missing @@ -4823,8 +4862,187 @@ class DataSaveCleaner(Operation): dataOut.EIonTempFinal[0,id_aux:]=missing dataOut.PhyFinal[0,id_aux:]=missing dataOut.EPhyFinal[0,id_aux:]=missing + if (time_text.hour <= 8): #Year: 2022, DOY:243 + id_aux = 11 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + if (time_text.hour == 23): #Year: 2022, DOY:243 + id_aux = 12 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + if (time_text.hour == 5 and time_text.minute == 21): #Year: 2022, DOY:243 + id_aux = (36,37) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 5 and time_text.minute == 53): #Year: 2022, DOY:243 + id_aux = (37,38) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 6 and time_text.minute == 4): #Year: 2022, DOY:243 + id_aux = (38,39) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 12 and time_text.minute == 6): #Year: 2022, DOY:243 + id_aux = (29,30) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 14 and time_text.minute == 14): #Year: 2022, DOY:243 + id_aux = (35,36) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 23 and time_text.minute == 2): #Year: 2022, DOY:243 + id_aux = (41,42) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 0 and time_text.minute == 8): #Year: 2022, DOY:243 + id_aux = 33 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + id_aux = 18 + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 23 and time_text.minute == 26): #Year: 2022, DOY:243 + id_aux = (12,13,14) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 23 and time_text.minute == 36): #Year: 2022, DOY:243 + id_aux = (14,15,16) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 2 and time_text.minute == 6): #Year: 2022, DOY:243 + id_aux = (36,37,38) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 2 and time_text.minute == 16): #Year: 2022, DOY:243 + id_aux = (34,35) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 2 and time_text.minute == 38): #Year: 2022, DOY:243 + id_aux = (35,36) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 3 and time_text.minute == 20): #Year: 2022, DOY:243 + id_aux = (33,34) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 3 and time_text.minute == 42): #Year: 2022, DOY:243 + id_aux = 34 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + if (time_text.hour == 4 and time_text.minute == 35): #Year: 2022, DOY:243 + id_aux = (36,37) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing ''' - #''' + ''' if (time_text.hour == 2 and time_text.minute == 23): #Year: 2022, DOY:244 id_aux = 12 dataOut.DensityFinal[0,:id_aux]=missing @@ -4905,16 +5123,367 @@ class DataSaveCleaner(Operation): dataOut.EIonTempFinal[0,id_aux:]=missing dataOut.PhyFinal[0,id_aux:]=missing dataOut.EPhyFinal[0,id_aux:]=missing + if (time_text.hour <= 8): #Year: 2022, DOY:244 + id_aux = 12 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + if (time_text.hour == 23): #Year: 2022, DOY:244 + id_aux = 12 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + if (time_text.hour == 5 and time_text.minute == 42): #Year: 2022, DOY:244 + id_aux = (32,33) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 11 and time_text.minute == 56): #Year: 2022, DOY:244 + id_aux = (39,40) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 12 and time_text.minute == 52): #Year: 2022, DOY:244 + id_aux = (36,37) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 13 and time_text.minute == 3): #Year: 2022, DOY:244 + id_aux = (37,38) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 23 and time_text.minute == 11): #Year: 2022, DOY:244 + id_aux = (40,41) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 23 and time_text.minute == 21): #Year: 2022, DOY:244 + id_aux = (12,13,39,40,41) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 23 and time_text.minute == 53): #Year: 2022, DOY:244 + id_aux = (15,16,17,18) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 2 and time_text.minute == 44): #Year: 2022, DOY:244 + id_aux = (40,41,42) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 3 and time_text.minute == 37): #Year: 2022, DOY:244 + id_aux = (36,37) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 4 and time_text.minute == 9): #Year: 2022, DOY:244 + id_aux = (32,33,34) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 4 and time_text.minute == 20): #Year: 2022, DOY:244 + id_aux = 37 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + if (time_text.hour == 4 and time_text.minute == 31): #Year: 2022, DOY:244 + id_aux = 33 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + ''' + ''' + if (time_text.hour <= 10): #Year: 2022, DOY:245 + id_aux = 11 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + if (time_text.hour == 5 and time_text.minute == 10): #Year: 2022, DOY:245 + id_aux = 35 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + if (time_text.hour == 5 and time_text.minute == 21): #Year: 2022, DOY:245 + id_aux = 36 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + if (time_text.hour == 11 and time_text.minute == 45): #Year: 2022, DOY:245 + id_aux = 7 + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + ''' + ''' + if (time_text.hour == 23 and time_text.minute > 30): #Year: 2022, DOY:241 + id_aux = 17 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + + if (time_text.hour == 13 and time_text.minute == 36): #Year: 2022, DOY:241 + id_aux = 33 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + if (time_text.hour == 13 and time_text.minute == 47): #Year: 2022, DOY:241 + id_aux = 36 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + + if (time_text.hour == 13 and time_text.minute == 57): #Year: 2022, DOY:241 + id_aux = 36 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + ''' + #''' + #print("den: ", dataOut.DensityFinal[0,27]) + if (time_text.hour == 5 and time_text.minute == 42): #Year: 2022, DOY:242 + id_aux = 16 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + if (time_text.hour == 5 and time_text.minute == 53): #Year: 2022, DOY:242 + id_aux = 9 + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 6): #Year: 2022, DOY:242 + id_aux = 9 + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 6 and time_text.minute == 36): #Year: 2022, DOY:242 + id_aux = (10,36,37) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 7): #Year: 2022, DOY:242 + id_aux = 9 + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 13 and time_text.minute == 32): #Year: 2022, DOY:242 + id_aux = (36,37) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 23 or time_text.hour <= 4): #Year: 2022, DOY:242 + id_aux = 15 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + if (time_text.hour == 3 and time_text.minute == 13): #Year: 2022, DOY:242 + id_aux = (37,38) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 3 and time_text.minute == 34): #Year: 2022, DOY:242 + id_aux = (35,36) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 4 and time_text.minute == 17): #Year: 2022, DOY:242 + id_aux = (34,35) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 18 and time_text.minute == 30): #Year: 2022, DOY:242 + id_aux = (26,27) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 14 and time_text.minute == 14): #Year: 2022, DOY:242 + id_aux = (35,36) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing #''' #print("den_final",dataOut.DensityFinal) + dataOut.flagNoData = numpy.all(numpy.isnan(dataOut.DensityFinal)) #Si todos los valores son NaN no se prosigue ####dataOut.flagNoData = False #Solo para ploteo dataOut.DensityFinal *= 1.e6 #Convert units to m^⁻3 dataOut.EDensityFinal *= 1.e6 #Convert units to m^⁻3 - #print(dataOut.flagNoData) + print("Save Cleaner: ", dataOut.flagNoData) + #print("den: ", dataOut.DensityFinal[0,27]) return dataOut @@ -5704,6 +6273,8 @@ class SSheightProfiles(Operation): #exit(1) #print(dataOut.data[0,:,150]) #exit(1) + #print(dataOut.data[0,:,0]*numpy.conjugate(dataOut.data[0,0,0])) + #exit(1) return dataOut @@ -6905,6 +7476,7 @@ class CrossProdHybrid(CrossProdDP): #self.dataOut.nptsfft2=150 self.cnorm=float((dataOut.nProfiles-dataOut.NSCAN)/dataOut.NSCAN) self.lagp0=numpy.zeros((dataOut.NLAG,dataOut.NRANGE,dataOut.NAVG),'complex128') + ww=numpy.zeros((dataOut.NLAG,dataOut.NRANGE,dataOut.NSCAN,dataOut.NAVG),'complex128') self.lagp1=numpy.zeros((dataOut.NLAG,dataOut.NRANGE,dataOut.NAVG),'complex128') self.lagp2=numpy.zeros((dataOut.NLAG,dataOut.NRANGE,dataOut.NAVG),'complex128') self.lagp3=numpy.zeros((dataOut.NLAG,dataOut.NRANGE,dataOut.NAVG),'complex128') @@ -6927,6 +7499,7 @@ class CrossProdHybrid(CrossProdDP): #exit(1) if i==0: self.lagp0[n][j][self.bcounter-1]=numpy.sum(c[:dataOut.NSCAN]) + ww[n,j,:,self.bcounter-1]=c[:dataOut.NSCAN] self.lagp3[n][j][self.bcounter-1]=numpy.sum(c[dataOut.NSCAN:]/self.cnorm) elif i==1: self.lagp1[n][j][self.bcounter-1]=numpy.sum(c[:dataOut.NSCAN]) @@ -6945,6 +7518,9 @@ class CrossProdHybrid(CrossProdDP): #print(sum(self.buffer[3,:,199,2])) #print(self.cnorm) #exit(1) + #print("self,lagp0: ", self.lagp0[0,0,self.bcounter-1]) + print(ww[:,0,0,self.bcounter-1]) + exit(1) def LP_median_estimates_original(self,dataOut): @@ -6955,7 +7531,7 @@ class CrossProdHybrid(CrossProdDP): self.output=numpy.zeros((dataOut.NLAG,dataOut.NRANGE,dataOut.NR),'complex128') self.lag_products_LP_median_estimates_aux=0 - + #print("self,lagp0: ", numpy.sum(self.lagp0[0,0,:])) for i in range(dataOut.NLAG): for j in range(dataOut.NRANGE): for l in range(4): #four outputs @@ -6995,7 +7571,7 @@ class CrossProdHybrid(CrossProdDP): if self.lag_products_LP_median_estimates_aux==1: self.output=numpy.zeros((dataOut.NLAG,dataOut.NRANGE,dataOut.NR),'complex128') self.lag_products_LP_median_estimates_aux=0 - + #print("self,lagp0: ", numpy.sum(self.lagp0[0,0,:])) for i in range(dataOut.NLAG): #my_list = ([0,1,2,3,4,5,6,7]) #hasta 7 funciona, en 6 ya no @@ -7388,12 +7964,12 @@ class LongPulseAnalysis(Operation): self.aux=0 dataOut.cut=30 - for i in range(30,15,-1): + for i in range(30,15,-1): #Aquí se calcula en donde se unirá DP y LP en la parte final if numpy.nanmax(dataOut.acfs_error_to_plot[i,:])>=10 or dataOut.info2[i]==0: dataOut.cut=i-1 for i in range(dataOut.NLAG): - self.cal[i]=sum(dataOut.output_LP_integrated[i,:,3].real) + self.cal[i]=sum(dataOut.output_LP_integrated[i,:,3].real) #Lag x Height x Channel #print(numpy.sum(self.cal)) #Coinciden #exit(1) @@ -7401,7 +7977,18 @@ class LongPulseAnalysis(Operation): #print(anoise0) #print(anoise1) #exit(1) - + print("nis: ", dataOut.nis) + print("pan: ", dataOut.pan) + print("pbn: ", dataOut.pbn) + #print(numpy.sum(dataOut.output_LP_integrated[0,:,0])) + ''' + import matplotlib.pyplot as plt + plt.plot(dataOut.output_LP_integrated[:,40,0]) + plt.show() + ''' + #print(dataOut.output_LP_integrated[0,40,0]) + print(numpy.sum(dataOut.output_LP_integrated[:,0,0])) + exit(1) #################### PROBAR MÁS INTEGRACIÓN, SINO MODIFICAR VALOR DE "NIS" #################### # VER dataOut.nProfiles_LP # @@ -7427,7 +8014,7 @@ class LongPulseAnalysis(Operation): 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] + self.powera[j]=dataOut.output_LP_integrated.real[0,k,0] #Lag0 and Channel 0 ## examine drifts here - based on 60 'indep.' estimates #print(numpy.sum(self.powera)) @@ -7520,7 +8107,7 @@ class LongPulseAnalysis(Operation): self.perror[:range2_nnls]=1.00/self.perror[:range2_nnls] b_nnlswrap=numpy.zeros(range2_nnls,'float64') - b_nnlswrap[:]=numpy.matmul(self.powera[dataOut.IBITS+ioff:range1_nnls+dataOut.IBITS+ioff],g_nnlswrap) + b_nnlswrap[:]=numpy.matmul(self.powera[dataOut.IBITS+ioff:range1_nnls+dataOut.IBITS+ioff],g_nnlswrap) #match filter alturas x_nnlswrap=numpy.zeros(range2_nnls,'float64') x_nnlswrap[:]=nnls(a_nnlswrap,b_nnlswrap)[0] @@ -7630,6 +8217,8 @@ class LongPulseAnalysis(Operation): exit(1) ''' print("Success") + ###################Correlation pulse and itself + #print(dataOut.NRANGE) with suppress_stdout_stderr(): #pass @@ -7649,6 +8238,303 @@ class LongPulseAnalysis(Operation): return dataOut +class LongPulseAnalysisSpectra(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 + #print(anoise0) + #exit(1) + 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): #Aquí se calcula en donde se unirá DP y LP en la parte final + if numpy.nanmax(dataOut.acfs_error_to_plot[i,:])>=10 or dataOut.info2[i]==0: + dataOut.cut=i-1 + + for i in range(dataOut.NLAG): + self.cal[i]=sum(dataOut.output_LP_integrated[i,:,3].real) #Lag x Height x Channel + + #print(numpy.sum(self.cal)) #Coinciden + #exit(1) + self.cal/=float(dataOut.NRANGE) + + + #################### PROBAR MÁS INTEGRACIÓN, SINO MODIFICAR VALOR DE "NIS" #################### + # VER dataOut.nProfiles_LP # + + ''' + #PLOTEAR POTENCIA VS RUIDO, QUIZA SE ESTA REMOVIENDO MUCHA SEÑAL + #print(dataOut.heightList) + import matplotlib.pyplot as plt + plt.plot(10*numpy.log10(dataOut.output_LP_integrated.real[0,:,0]),dataOut.range1) + #plt.plot(10*numpy.log10(dataOut.output_LP_integrated.real[0,:,0]/dataOut.nProfiles_LP),dataOut.range1) + plt.axvline(10*numpy.log10(anoise0),color='k',linestyle='dashed') + plt.grid() + plt.xlim(20,100) + plt.show() + ''' + + + 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] #Lag0 and Channel 0 + + ## examine drifts here - based on 60 'indep.' estimates + #print(numpy.sum(self.powera)) + #exit(1) + #nis=dataOut.NSCAN*dataOut.NAVG*dataOut.nint*10 + nis = dataOut.nis + #print("nis",nis) + 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 + #print(gamma,beta) + #exit(1) + for i in range(1,3): + gamma=3.0/(2.0*numpy.pi*dataOut.lags_LP[i]*1.0e-3) + #print("gamma",gamma) + 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]) + #print("1",dataOut.output_LP_integrated.imag[i,j,0]) + #print("2",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]) + + #print(dataOut.dphi2) + #exit(1) + + 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 + nis = dataOut.nis/10 + #print("nis",nis) + + self.sigma[:dataOut.NACF+2*dataOut.IBITS+2]=((anoise0+self.powera[:dataOut.NACF+2*dataOut.IBITS+2])**2)/float(nis) + #print(self.sigma) + #exit(1) + 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 + #print(dataOut.h2) + #print(numpy.sum(alpha)) + #print(numpy.sum(dataOut.ph2)) + 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 + #print(numpy.sum(self.powerb)) + #print(numpy.sum(dataOut.ene)) + #print(numpy.sum(dataOut.ne)) + #exit(1) + #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)) + + print("Success") + #print(dataOut.NRANGE) + with suppress_stdout_stderr(): + pass + #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) + + print("status: ",dataOut.status) + + 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 LongPulseAnalysis_V2(Operation): """Operation to estimate ACFs, temperatures, total electron density and Hydrogen/Helium fractions from the Long Pulse data.