##// END OF EJS Templates
schain
RSS
LP Faraday update
rflores -
r1542:9f6529c5475c
parent child
Show More
Side by Side Unified
Context file:
r1542:9f6529c5475c Loading diff...:
Show file before | Show file after | Hide comments
@@ -653,6 +653,7 class Project(Process):
653 elif not ok:
653 elif not ok:
654 break
654 break
655 #print("****************************************************end")
655 #print("****************************************************end")
656 #exit(1)
656 if n == 0:
657 if n == 0:
657 err = True
658 err = True
658
659
Show file before | Show file after | Hide comments
@@ -42,6 +42,7 class SpectraPlot(Plot):
42 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
42 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
43 data['spc'] = spc
43 data['spc'] = spc
44 data['rti'] = dataOut.getPower()
44 data['rti'] = dataOut.getPower()
45 #print("NormFactor: ",dataOut.normFactor)
45 #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
46 #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
46 if hasattr(dataOut, 'LagPlot'): #Double Pulse
47 if hasattr(dataOut, 'LagPlot'): #Double Pulse
47 max_hei_id = dataOut.nHeights - 2*dataOut.LagPlot
48 max_hei_id = dataOut.nHeights - 2*dataOut.LagPlot
@@ -101,10 +102,11 class SpectraPlot(Plot):
101 self.xmin = self.xmin if self.xmin else numpy.nanmin(x)#-self.xmax
102 self.xmin = self.xmin if self.xmin else numpy.nanmin(x)#-self.xmax
102 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
103 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
103 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
104 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
105
104 ax.plt = ax.pcolormesh(x, y, z[n].T,
106 ax.plt = ax.pcolormesh(x, y, z[n].T,
105 vmin=self.zmin,
107 vmin=self.zmin,
106 vmax=self.zmax,
108 vmax=self.zmax,
107 cmap=plt.get_cmap(self.colormap)
109 cmap=plt.get_cmap(self.colormap),
108 )
110 )
109
111
110 if self.showprofile:
112 if self.showprofile:
Show file before | Show file after | Hide comments
@@ -263,8 +263,6 class DenRTIPlot(RTIPlot):
263 )
263 )
264
264
265
265
266
267
268 class ETempRTIPlot(RTIPlot):
266 class ETempRTIPlot(RTIPlot):
269
267
270 '''
268 '''
@@ -351,7 +349,6 class ETempRTIPlot(RTIPlot):
351 )
349 )
352
350
353
351
354
355 class ITempRTIPlot(ETempRTIPlot):
352 class ITempRTIPlot(ETempRTIPlot):
356
353
357 '''
354 '''
@@ -372,7 +369,6 class ITempRTIPlot(ETempRTIPlot):
372 return data, meta
369 return data, meta
373
370
374
371
375
376 class HFracRTIPlot(ETempRTIPlot):
372 class HFracRTIPlot(ETempRTIPlot):
377
373
378 '''
374 '''
@@ -545,6 +541,8 class TempsHPPlot(Plot):
545 ax.errorbar(Ti, self.y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
541 ax.errorbar(Ti, self.y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
546 plt.legend(loc='lower right')
542 plt.legend(loc='lower right')
547 ax.yaxis.set_minor_locator(MultipleLocator(15))
543 ax.yaxis.set_minor_locator(MultipleLocator(15))
544 ax.grid(which='minor')
545
548
546
549 class FracsHPPlot(Plot):
547 class FracsHPPlot(Plot):
550 '''
548 '''
Show file before | Show file after | Hide comments
@@ -70,14 +70,16 class ProcessingUnit(object):
70 def call(self, **kwargs):
70 def call(self, **kwargs):
71 '''
71 '''
72 '''
72 '''
73
73 #print("call")
74 try:
74 try:
75 #print("try")
76 if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
75 if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
77 #print("Try")
76 #if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error and not self.dataIn.runNextUnit:
78 #print(self.dataIn)
77 if self.dataIn.runNextUnit:
79 #print(self.dataIn.flagNoData)
78 #print("SUCCESSSSSSS")
80 return self.dataIn.isReady()
79 #exit(1)
80 return not self.dataIn.isReady()
81 else:
82 return self.dataIn.isReady()
81 elif self.dataIn is None or not self.dataIn.error:
83 elif self.dataIn is None or not self.dataIn.error:
82 #print([getattr(self, at) for at in self.inputs])
84 #print([getattr(self, at) for at in self.inputs])
83 #print("Elif 1")
85 #print("Elif 1")
@@ -112,6 +114,7 class ProcessingUnit(object):
112 #elif optype == 'external' and self.dataOut.isReady():
114 #elif optype == 'external' and self.dataOut.isReady():
113 #op.queue.put(copy.deepcopy(self.dataOut))
115 #op.queue.put(copy.deepcopy(self.dataOut))
114 #print(not self.dataOut.isReady())
116 #print(not self.dataOut.isReady())
117
115 try:
118 try:
116 if self.dataOut.runNextUnit:
119 if self.dataOut.runNextUnit:
117 runNextUnit = self.dataOut.runNextUnit
120 runNextUnit = self.dataOut.runNextUnit
@@ -125,6 +128,7 class ProcessingUnit(object):
125 #if not self.dataOut.isReady():
128 #if not self.dataOut.isReady():
126 #return 'Error' if self.dataOut.error else input()
129 #return 'Error' if self.dataOut.error else input()
127 #print("NexT",runNextUnit)
130 #print("NexT",runNextUnit)
131 #print("error: ",self.dataOut.error)
128 return 'Error' if self.dataOut.error else runNextUnit# self.dataOut.isReady()
132 return 'Error' if self.dataOut.error else runNextUnit# self.dataOut.isReady()
129
133
130 def setup(self):
134 def setup(self):
Show file before | Show file after | Hide comments
@@ -1164,18 +1164,34 class Oblique_Gauss_Fit(Operation):
1164 #popt = least_squares(lsq_func,[a1,b1,c1,a2,b2,c2,k2,d],x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=1)
1164 #popt = least_squares(lsq_func,[a1,b1,c1,a2,b2,c2,k2,d],x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=1)
1165 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1165 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1166 # popt = least_squares(lsq_func,[a1,b1,c1,a2,b2,c2,k2,d],x_scale=params_scale,verbose=1)
1166 # popt = least_squares(lsq_func,[a1,b1,c1,a2,b2,c2,k2,d],x_scale=params_scale,verbose=1)
1167 #print(popt)
1168
1169 J = popt.jac
1170
1171 try:
1172 cov = numpy.linalg.inv(J.T.dot(J))
1173 error = numpy.sqrt(numpy.diagonal(cov))
1174 except:
1175 error = numpy.ones((9))*numpy.NAN
1176 #print("error_inside",error)
1177 #exit(1)
1167
1178
1168 A1f = popt.x[0]; B1f = popt.x[1]; C1f = popt.x[2]; K1f = popt.x[3]
1179 A1f = popt.x[0]; B1f = popt.x[1]; C1f = popt.x[2]; K1f = popt.x[3]
1169 A2f = popt.x[4]; B2f = popt.x[5]; C2f = popt.x[6]; K2f = popt.x[7]
1180 A2f = popt.x[4]; B2f = popt.x[5]; C2f = popt.x[6]; K2f = popt.x[7]
1170 Df = popt.x[8]
1181 Df = popt.x[8]
1171
1182 '''
1183 A1f_err = error.x[0]; B1f_err= error.x[1]; C1f_err = error.x[2]; K1f_err = error.x[3]
1184 A2f_err = error.x[4]; B2f_err = error.x[5]; C2f_err = error.x[6]; K2f_err = error.x[7]
1185 Df_err = error.x[8]
1186 '''
1172 aux1 = self.gaussian_skew(freq, A1f, B1f, C1f, K1f, Df)
1187 aux1 = self.gaussian_skew(freq, A1f, B1f, C1f, K1f, Df)
1173 doppler1 = freq[numpy.argmax(aux1)]
1188 doppler1 = freq[numpy.argmax(aux1)]
1174
1189
1175 aux2 = self.gaussian_skew(freq, A2f, B2f, C2f, K2f, Df)
1190 aux2 = self.gaussian_skew(freq, A2f, B2f, C2f, K2f, Df)
1176 doppler2 = freq[numpy.argmax(aux2)]
1191 doppler2 = freq[numpy.argmax(aux2)]
1177
1192 #print("error",error)
1178 return A1f, B1f, C1f, K1f, A2f, B2f, C2f, K2f, Df, doppler1, doppler2
1193 #exit(1)
1194 return A1f, B1f, C1f, K1f, A2f, B2f, C2f, K2f, Df, doppler1, doppler2, error
1179
1195
1180 def Double_Gauss_Double_Skew_fit_weight_bound_with_inputs(self, spc, freq, a1, b1, c1, a2, b2, c2, k2, d):
1196 def Double_Gauss_Double_Skew_fit_weight_bound_with_inputs(self, spc, freq, a1, b1, c1, a2, b2, c2, k2, d):
1181
1197
@@ -1276,6 +1292,7 class Oblique_Gauss_Fit(Operation):
1276 dataOut.Oblique_params = numpy.ones((1,10,dataOut.nHeights))*numpy.NAN
1292 dataOut.Oblique_params = numpy.ones((1,10,dataOut.nHeights))*numpy.NAN
1277 elif mode == 9:
1293 elif mode == 9:
1278 dataOut.Oblique_params = numpy.ones((1,11,dataOut.nHeights))*numpy.NAN
1294 dataOut.Oblique_params = numpy.ones((1,11,dataOut.nHeights))*numpy.NAN
1295 dataOut.Oblique_param_errors = numpy.ones((1,9,dataOut.nHeights))*numpy.NAN
1279
1296
1280 dataOut.VelRange = x
1297 dataOut.VelRange = x
1281
1298
@@ -1303,6 +1320,7 class Oblique_Gauss_Fit(Operation):
1303 else:
1320 else:
1304
1321
1305 for hei in itertools.chain(l1, l2):
1322 for hei in itertools.chain(l1, l2):
1323 #for hei in range(79,81):
1306
1324
1307 try:
1325 try:
1308
1326
@@ -1322,11 +1340,14 class Oblique_Gauss_Fit(Operation):
1322 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,9,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei]))
1340 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,9,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei]))
1323
1341
1324 elif mode == 9: #Double Skewed Weighted Bounded no inputs
1342 elif mode == 9: #Double Skewed Weighted Bounded no inputs
1325 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] = self.Double_Gauss_Double_Skew_fit_weight_bound_no_inputs(spc,x)
1343 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)
1326 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,10,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei]))
1344 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,10,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei]))
1327 #print(hei)
1345 #print(hei)
1328 #print(dataOut.Oblique_params[0,10,hei])
1346 #print(dataOut.Oblique_params[0,10,hei])
1329 #print(dataOut.dplr_2_u[0,0,hei])
1347 #print(dataOut.dplr_2_u[0,0,hei])
1348 #print("outside",dataOut.Oblique_param_errors[0,:,hei])
1349 #print("SUCCESSSSSSS")
1350 #exit(1)
1330
1351
1331 else:
1352 else:
1332 spc_fit, A1, B1, C1, D1 = self.Gauss_fit_2(spc,x,'first')
1353 spc_fit, A1, B1, C1, D1 = self.Gauss_fit_2(spc,x,'first')
@@ -5575,7 +5596,7 class MergeProc(ProcessingUnit):
5575 def __init__(self):
5596 def __init__(self):
5576 ProcessingUnit.__init__(self)
5597 ProcessingUnit.__init__(self)
5577
5598
5578 def run(self, attr_data, attr_data_2 = None, mode=0):
5599 def run(self, attr_data, attr_data_2 = None, attr_data_3 = None, attr_data_4 = None, attr_data_5 = None, mode=0):
5579
5600
5580 self.dataOut = getattr(self, self.inputs[0])
5601 self.dataOut = getattr(self, self.inputs[0])
5581 data_inputs = [getattr(self, attr) for attr in self.inputs]
5602 data_inputs = [getattr(self, attr) for attr in self.inputs]
@@ -5636,3 +5657,62 class MergeProc(ProcessingUnit):
5636 self.dataOut.freqRange = self.dataOut.getFreqRange(1)/1000.
5657 self.dataOut.freqRange = self.dataOut.getFreqRange(1)/1000.
5637
5658
5638 #exit(1)
5659 #exit(1)
5660
5661 if mode==4: #Hybrid LP-SSheightProfiles
5662 #data = numpy.concatenate([getattr(data, attr_data) for data in data_inputs],axis=1)
5663 #setattr(self.dataOut, attr_data, data)
5664 setattr(self.dataOut, 'dataLag_spc', getattr(data_inputs[0], attr_data)) #DP
5665 setattr(self.dataOut, 'dataLag_cspc', getattr(data_inputs[0], attr_data_2)) #DP
5666 setattr(self.dataOut, 'dataLag_spc_LP', getattr(data_inputs[1], attr_data_3)) #LP
5667 #setattr(self.dataOut, 'dataLag_cspc_LP', getattr(data_inputs[1], attr_data_4)) #LP
5668 #setattr(self.dataOut, 'data_acf', getattr(data_inputs[1], attr_data_5)) #LP
5669 setattr(self.dataOut, 'data_acf', getattr(data_inputs[1], attr_data_5)) #LP
5670 #print("Merge data_acf: ",self.dataOut.data_acf.shape)
5671 #exit(1)
5672 #print(self.dataOut.data_spc_LP.shape)
5673 #print("Exit")
5674 #exit(1)
5675 #setattr(self.dataOut, 'dataLag_cspc', [getattr(data, attr_data_2) for data in data_inputs][0])
5676 #setattr(self.dataOut, 'dataLag_cspc_LP', [getattr(data, attr_data_2) for data in data_inputs][1])
5677 #setattr(self.dataOut, 'nIncohInt', [getattr(data, attr_data_3) for data in data_inputs][0])
5678 #setattr(self.dataOut, 'nIncohInt_LP', [getattr(data, attr_data_3) for data in data_inputs][1])
5679 '''
5680 print(self.dataOut.dataLag_spc_LP.shape)
5681 print(self.dataOut.dataLag_cspc_LP.shape)
5682 exit(1)
5683 '''
5684 '''
5685 print(self.dataOut.dataLag_spc_LP[0,:,100])
5686 print(self.dataOut.dataLag_spc_LP[1,:,100])
5687 exit(1)
5688 '''
5689 #self.dataOut.dataLag_spc_LP = numpy.transpose(self.dataOut.dataLag_spc_LP[0],(2,0,1))
5690 #self.dataOut.dataLag_cspc_LP = numpy.transpose(self.dataOut.dataLag_cspc_LP,(3,1,2,0))
5691 '''
5692 print("Merge")
5693 print(numpy.shape(self.dataOut.dataLag_spc))
5694 print(numpy.shape(self.dataOut.dataLag_spc_LP))
5695 print(numpy.shape(self.dataOut.dataLag_cspc))
5696 print(numpy.shape(self.dataOut.dataLag_cspc_LP))
5697 exit(1)
5698 '''
5699 #print(numpy.sum(self.dataOut.dataLag_spc_LP[2,:,164])/128)
5700 #print(numpy.sum(self.dataOut.dataLag_cspc_LP[0,:,30,1])/128)
5701 #exit(1)
5702 #print(self.dataOut.NDP)
5703 #print(self.dataOut.nNoiseProfiles)
5704
5705 #self.dataOut.nIncohInt_LP = 128
5706 #self.dataOut.nProfiles_LP = 128#self.dataOut.nIncohInt_LP
5707 self.dataOut.nProfiles_LP = 16#28#self.dataOut.nIncohInt_LP
5708 self.dataOut.nProfiles_LP = self.dataOut.data_acf.shape[1]#28#self.dataOut.nIncohInt_LP
5709 self.dataOut.NSCAN = 128
5710 self.dataOut.nIncohInt_LP = self.dataOut.nIncohInt*self.dataOut.NSCAN
5711 #print("sahpi",self.dataOut.nIncohInt_LP)
5712 #exit(1)
5713 self.dataOut.NLAG = 16
5714 self.dataOut.NRANGE = self.dataOut.data_acf.shape[-1]
5715
5716 #print(numpy.shape(self.dataOut.data_spc))
5717
5718 #exit(1)
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (859 lines changed) Show them Hide them
@@ -4,7 +4,7 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
4 from schainpy.model.data.jrodata import Spectra
4 from schainpy.model.data.jrodata import Spectra
5 from schainpy.model.data.jrodata import hildebrand_sekhon
5 from schainpy.model.data.jrodata import hildebrand_sekhon
6
6
7 class SpectraAFCProc(ProcessingUnit):
7 class SpectraAFCProc_V0(ProcessingUnit):
8
8
9 def __init__(self, **kwargs):
9 def __init__(self, **kwargs):
10
10
@@ -140,12 +140,16 class SpectraAFCProc(ProcessingUnit):
140
140
141 if self.dataIn.type == "Spectra":
141 if self.dataIn.type == "Spectra":
142 self.dataOut.copy(self.dataIn)
142 self.dataOut.copy(self.dataIn)
143 #print(self.dataOut.data.shape)
144 #exit(1)
143 spc = self.dataOut.data_spc
145 spc = self.dataOut.data_spc
144 data = numpy.fft.fftshift( spc, axes=(1,))
146 data = numpy.fft.fftshift( spc, axes=(1,))
145 data = numpy.fft.ifft(data, axis=1)
147 data = numpy.fft.ifft(data, axis=1)
146 #data = numpy.fft.fftshift( data, axes=(1,))
148 data = numpy.fft.fftshift( data, axes=(1,))
147 #acf = numpy.abs(data) # Autocorrelacion LLAMAR A ESTE VALOR ACF
149 acf = numpy.abs(data) # Autocorrelacion LLAMAR A ESTE VALOR ACF
148 acf = data
150 acf = data #Comentarlo?
151 print("acf",acf[0,:,150])
152 exit(1)
149 #'''
153 #'''
150 if real:
154 if real:
151 acf = data.real
155 acf = data.real
@@ -167,7 +171,7 class SpectraAFCProc(ProcessingUnit):
167 acf[i,:,j]= acf[i,:,j] / numpy.max(numpy.abs(acf[i,:,j]))
171 acf[i,:,j]= acf[i,:,j] / numpy.max(numpy.abs(acf[i,:,j]))
168 '''
172 '''
169 self.dataOut.data_acf = acf
173 self.dataOut.data_acf = acf
170 self.dataOut.data_spc = acf
174 self.dataOut.data_spc = acf.real
171 #print(self.dataOut.data_acf[0,:,0])
175 #print(self.dataOut.data_acf[0,:,0])
172 #exit(1)
176 #exit(1)
173 '''
177 '''
@@ -183,6 +187,7 class SpectraAFCProc(ProcessingUnit):
183 #print(self.dataOut.data_spc.shape)
187 #print(self.dataOut.data_spc.shape)
184 #print(self.dataOut.data_acf.shape)
188 #print(self.dataOut.data_acf.shape)
185 '''
189 '''
190 '''
186 import matplotlib.pyplot as plt
191 import matplotlib.pyplot as plt
187 hei = 10
192 hei = 10
188 #print(self.dataOut.heightList)
193 #print(self.dataOut.heightList)
@@ -197,7 +202,851 class SpectraAFCProc(ProcessingUnit):
197 plt.ylim(0,1000)
202 plt.ylim(0,1000)
198 plt.show()
203 plt.show()
199 exit(1)
204 exit(1)
205 '''
206 return True
207
208 if code is not None:
209 self.code = numpy.array(code).reshape(nCode,nBaud)
210 else:
211 self.code = None
212
213 if self.dataIn.type == "Voltage":
214
215 if nFFTPoints == None:
216 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
217
218 if nProfiles == None:
219 nProfiles = nFFTPoints
220
221 self.dataOut.ippFactor = 1
222
223 self.dataOut.nFFTPoints = nFFTPoints
224 self.dataOut.nProfiles = nProfiles
225 self.dataOut.pairsList = pairsList
226
227 # if self.buffer is None:
228 # self.buffer = numpy.zeros( (self.dataIn.nChannels, nProfiles, self.dataIn.nHeights),
229 # dtype='complex')
230
231 if not self.dataIn.flagDataAsBlock:
232 self.buffer = self.dataIn.data.copy()
233
234 # for i in range(self.dataIn.nHeights):
235 # self.buffer[:, self.profIndex, self.profIndex:] = voltage_data[:,:self.dataIn.nHeights - self.profIndex]
236 #
237 # self.profIndex += 1
238
239 else:
240 raise ValueError("")
241
242 self.firstdatatime = self.dataIn.utctime
243
244 self.profIndex == nProfiles
245
246 self.__updateSpecFromVoltage()
247
248 self.__getFft()
249
250 self.dataOut.flagNoData = False
251
252 return True
253
254 raise ValueError("The type of input object '%s' is not valid"%(self.dataIn.type))
255
256 def __selectPairs(self, pairsList):
257
258 if channelList == None:
259 return
260
261 pairsIndexListSelected = []
262
263 for thisPair in pairsList:
264
265 if thisPair not in self.dataOut.pairsList:
266 continue
267
268 pairIndex = self.dataOut.pairsList.index(thisPair)
269
270 pairsIndexListSelected.append(pairIndex)
271
272 if not pairsIndexListSelected:
273 self.dataOut.data_cspc = None
274 self.dataOut.pairsList = []
275 return
276
277 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
278 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
279
280 return
281
282 def __selectPairsByChannel(self, channelList=None):
283
284 if channelList == None:
285 return
286
287 pairsIndexListSelected = []
288 for pairIndex in self.dataOut.pairsIndexList:
289 #First pair
290 if self.dataOut.pairsList[pairIndex][0] not in channelList:
291 continue
292 #Second pair
293 if self.dataOut.pairsList[pairIndex][1] not in channelList:
294 continue
295
296 pairsIndexListSelected.append(pairIndex)
297
298 if not pairsIndexListSelected:
299 self.dataOut.data_cspc = None
300 self.dataOut.pairsList = []
301 return
302
303 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
304 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
305
306 return
307
308 def selectChannels(self, channelList):
309
310 channelIndexList = []
311
312 for channel in channelList:
313 if channel not in self.dataOut.channelList:
314 raise ValueError("Error selecting channels, Channel %d is not valid.\nAvailable channels = %s" %(channel, str(self.dataOut.channelList)))
315
316 index = self.dataOut.channelList.index(channel)
317 channelIndexList.append(index)
318
319 self.selectChannelsByIndex(channelIndexList)
320
321 def selectChannelsByIndex(self, channelIndexList):
322 """
323 Selecciona un bloque de datos en base a canales segun el channelIndexList
324
325 Input:
326 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
327
328 Affected:
329 self.dataOut.data_spc
330 self.dataOut.channelIndexList
331 self.dataOut.nChannels
332
333 Return:
334 None
335 """
336
337 for channelIndex in channelIndexList:
338 if channelIndex not in self.dataOut.channelIndexList:
339 raise ValueError("Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " %(channelIndex, self.dataOut.channelIndexList))
340
341 # nChannels = len(channelIndexList)
342
343 data_spc = self.dataOut.data_spc[channelIndexList,:]
344 data_dc = self.dataOut.data_dc[channelIndexList,:]
345
346 self.dataOut.data_spc = data_spc
347 self.dataOut.data_dc = data_dc
348
349 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
350 # self.dataOut.nChannels = nChannels
351
352 self.__selectPairsByChannel(self.dataOut.channelList)
353
354 return 1
355
356 def selectHeights(self, minHei, maxHei):
357 """
358 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
359 minHei <= height <= maxHei
360
361 Input:
362 minHei : valor minimo de altura a considerar
363 maxHei : valor maximo de altura a considerar
364
365 Affected:
366 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
367
368 Return:
369 1 si el metodo se ejecuto con exito caso contrario devuelve 0
370 """
371
372 if (minHei > maxHei):
373 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei))
374
375 if (minHei < self.dataOut.heightList[0]):
376 minHei = self.dataOut.heightList[0]
377
378 if (maxHei > self.dataOut.heightList[-1]):
379 maxHei = self.dataOut.heightList[-1]
380
381 minIndex = 0
382 maxIndex = 0
383 heights = self.dataOut.heightList
384
385 inda = numpy.where(heights >= minHei)
386 indb = numpy.where(heights <= maxHei)
387
388 try:
389 minIndex = inda[0][0]
390 except:
391 minIndex = 0
392
393 try:
394 maxIndex = indb[0][-1]
395 except:
396 maxIndex = len(heights)
397
398 self.selectHeightsByIndex(minIndex, maxIndex)
399
400 return 1
401
402 def getBeaconSignal(self, tauindex = 0, channelindex = 0, hei_ref=None):
403 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
404
405 if hei_ref != None:
406 newheis = numpy.where(self.dataOut.heightList>hei_ref)
407
408 minIndex = min(newheis[0])
409 maxIndex = max(newheis[0])
410 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
411 heightList = self.dataOut.heightList[minIndex:maxIndex+1]
412
413 # determina indices
414 nheis = int(self.dataOut.radarControllerHeaderObj.txB/(self.dataOut.heightList[1]-self.dataOut.heightList[0]))
415 avg_dB = 10*numpy.log10(numpy.sum(data_spc[channelindex,:,:],axis=0))
416 beacon_dB = numpy.sort(avg_dB)[-nheis:]
417 beacon_heiIndexList = []
418 for val in avg_dB.tolist():
419 if val >= beacon_dB[0]:
420 beacon_heiIndexList.append(avg_dB.tolist().index(val))
421
422 #data_spc = data_spc[:,:,beacon_heiIndexList]
423 data_cspc = None
424 if self.dataOut.data_cspc is not None:
425 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
426 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
427
428 data_dc = None
429 if self.dataOut.data_dc is not None:
430 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
431 #data_dc = data_dc[:,beacon_heiIndexList]
432
433 self.dataOut.data_spc = data_spc
434 self.dataOut.data_cspc = data_cspc
435 self.dataOut.data_dc = data_dc
436 self.dataOut.heightList = heightList
437 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
438
439 return 1
440
441
442 def selectHeightsByIndex(self, minIndex, maxIndex):
443 """
444 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
445 minIndex <= index <= maxIndex
446
447 Input:
448 minIndex : valor de indice minimo de altura a considerar
449 maxIndex : valor de indice maximo de altura a considerar
450
451 Affected:
452 self.dataOut.data_spc
453 self.dataOut.data_cspc
454 self.dataOut.data_dc
455 self.dataOut.heightList
200
456
457 Return:
458 1 si el metodo se ejecuto con exito caso contrario devuelve 0
459 """
460
461 if (minIndex < 0) or (minIndex > maxIndex):
462 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
463
464 if (maxIndex >= self.dataOut.nHeights):
465 maxIndex = self.dataOut.nHeights-1
466
467 #Spectra
468 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
469
470 data_cspc = None
471 if self.dataOut.data_cspc is not None:
472 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
473
474 data_dc = None
475 if self.dataOut.data_dc is not None:
476 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
477
478 self.dataOut.data_spc = data_spc
479 self.dataOut.data_cspc = data_cspc
480 self.dataOut.data_dc = data_dc
481
482 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
483
484 return 1
485
486 def removeDC(self, mode = 2):
487 jspectra = self.dataOut.data_spc
488 jcspectra = self.dataOut.data_cspc
489
490
491 num_chan = jspectra.shape[0]
492 num_hei = jspectra.shape[2]
493
494 if jcspectra is not None:
495 jcspectraExist = True
496 num_pairs = jcspectra.shape[0]
497 else: jcspectraExist = False
498
499 freq_dc = jspectra.shape[1]/2
500 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
501
502 if ind_vel[0]<0:
503 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
504
505 if mode == 1:
506 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
507
508 if jcspectraExist:
509 jcspectra[:,freq_dc,:] = (jcspectra[:,ind_vel[1],:] + jcspectra[:,ind_vel[2],:])/2
510
511 if mode == 2:
512
513 vel = numpy.array([-2,-1,1,2])
514 xx = numpy.zeros([4,4])
515
516 for fil in range(4):
517 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
518
519 xx_inv = numpy.linalg.inv(xx)
520 xx_aux = xx_inv[0,:]
521
522 for ich in range(num_chan):
523 yy = jspectra[ich,ind_vel,:]
524 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
525
526 junkid = jspectra[ich,freq_dc,:]<=0
527 cjunkid = sum(junkid)
528
529 if cjunkid.any():
530 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
531
532 if jcspectraExist:
533 for ip in range(num_pairs):
534 yy = jcspectra[ip,ind_vel,:]
535 jcspectra[ip,freq_dc,:] = numpy.dot(xx_aux,yy)
536
537
538 self.dataOut.data_spc = jspectra
539 self.dataOut.data_cspc = jcspectra
540
541 return 1
542
543 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
544
545 jspectra = self.dataOut.data_spc
546 jcspectra = self.dataOut.data_cspc
547 jnoise = self.dataOut.getNoise()
548 num_incoh = self.dataOut.nIncohInt
549
550 num_channel = jspectra.shape[0]
551 num_prof = jspectra.shape[1]
552 num_hei = jspectra.shape[2]
553
554 #hei_interf
555 if hei_interf is None:
556 count_hei = num_hei/2 #Como es entero no importa
557 hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
558 hei_interf = numpy.asarray(hei_interf)[0]
559 #nhei_interf
560 if (nhei_interf == None):
561 nhei_interf = 5
562 if (nhei_interf < 1):
563 nhei_interf = 1
564 if (nhei_interf > count_hei):
565 nhei_interf = count_hei
566 if (offhei_interf == None):
567 offhei_interf = 0
568
569 ind_hei = range(num_hei)
570 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
571 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
572 mask_prof = numpy.asarray(range(num_prof))
573 num_mask_prof = mask_prof.size
574 comp_mask_prof = [0, num_prof/2]
575
576
577 #noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
578 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
579 jnoise = numpy.nan
580 noise_exist = jnoise[0] < numpy.Inf
581
582 #Subrutina de Remocion de la Interferencia
583 for ich in range(num_channel):
584 #Se ordena los espectros segun su potencia (menor a mayor)
585 power = jspectra[ich,mask_prof,:]
586 power = power[:,hei_interf]
587 power = power.sum(axis = 0)
588 psort = power.ravel().argsort()
589
590 #Se estima la interferencia promedio en los Espectros de Potencia empleando
591 junkspc_interf = jspectra[ich,:,hei_interf[psort[range(offhei_interf, nhei_interf + offhei_interf)]]]
592
593 if noise_exist:
594 # tmp_noise = jnoise[ich] / num_prof
595 tmp_noise = jnoise[ich]
596 junkspc_interf = junkspc_interf - tmp_noise
597 #junkspc_interf[:,comp_mask_prof] = 0
598
599 jspc_interf = junkspc_interf.sum(axis = 0) / nhei_interf
600 jspc_interf = jspc_interf.transpose()
601 #Calculando el espectro de interferencia promedio
602 noiseid = numpy.where(jspc_interf <= tmp_noise/ numpy.sqrt(num_incoh))
603 noiseid = noiseid[0]
604 cnoiseid = noiseid.size
605 interfid = numpy.where(jspc_interf > tmp_noise/ numpy.sqrt(num_incoh))
606 interfid = interfid[0]
607 cinterfid = interfid.size
608
609 if (cnoiseid > 0): jspc_interf[noiseid] = 0
610
611 #Expandiendo los perfiles a limpiar
612 if (cinterfid > 0):
613 new_interfid = (numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof)%num_prof
614 new_interfid = numpy.asarray(new_interfid)
615 new_interfid = {x for x in new_interfid}
616 new_interfid = numpy.array(list(new_interfid))
617 new_cinterfid = new_interfid.size
618 else: new_cinterfid = 0
619
620 for ip in range(new_cinterfid):
621 ind = junkspc_interf[:,new_interfid[ip]].ravel().argsort()
622 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf/2],new_interfid[ip]]
623
624
625 jspectra[ich,:,ind_hei] = jspectra[ich,:,ind_hei] - jspc_interf #Corregir indices
626
627 #Removiendo la interferencia del punto de mayor interferencia
628 ListAux = jspc_interf[mask_prof].tolist()
629 maxid = ListAux.index(max(ListAux))
630
631
632 if cinterfid > 0:
633 for ip in range(cinterfid*(interf == 2) - 1):
634 ind = (jspectra[ich,interfid[ip],:] < tmp_noise*(1 + 1/numpy.sqrt(num_incoh))).nonzero()
635 cind = len(ind)
636
637 if (cind > 0):
638 jspectra[ich,interfid[ip],ind] = tmp_noise*(1 + (numpy.random.uniform(cind) - 0.5)/numpy.sqrt(num_incoh))
639
640 ind = numpy.array([-2,-1,1,2])
641 xx = numpy.zeros([4,4])
642
643 for id1 in range(4):
644 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
645
646 xx_inv = numpy.linalg.inv(xx)
647 xx = xx_inv[:,0]
648 ind = (ind + maxid + num_mask_prof)%num_mask_prof
649 yy = jspectra[ich,mask_prof[ind],:]
650 jspectra[ich,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
651
652
653 indAux = (jspectra[ich,:,:] < tmp_noise*(1-1/numpy.sqrt(num_incoh))).nonzero()
654 jspectra[ich,indAux[0],indAux[1]] = tmp_noise * (1 - 1/numpy.sqrt(num_incoh))
655
656 #Remocion de Interferencia en el Cross Spectra
657 if jcspectra is None: return jspectra, jcspectra
658 num_pairs = jcspectra.size/(num_prof*num_hei)
659 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
660
661 for ip in range(num_pairs):
662
663 #-------------------------------------------
664
665 cspower = numpy.abs(jcspectra[ip,mask_prof,:])
666 cspower = cspower[:,hei_interf]
667 cspower = cspower.sum(axis = 0)
668
669 cspsort = cspower.ravel().argsort()
670 junkcspc_interf = jcspectra[ip,:,hei_interf[cspsort[range(offhei_interf, nhei_interf + offhei_interf)]]]
671 junkcspc_interf = junkcspc_interf.transpose()
672 jcspc_interf = junkcspc_interf.sum(axis = 1)/nhei_interf
673
674 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
675
676 median_real = numpy.median(numpy.real(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
677 median_imag = numpy.median(numpy.imag(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
678 junkcspc_interf[comp_mask_prof,:] = numpy.complex(median_real, median_imag)
679
680 for iprof in range(num_prof):
681 ind = numpy.abs(junkcspc_interf[iprof,:]).ravel().argsort()
682 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf/2]]
683
684 #Removiendo la Interferencia
685 jcspectra[ip,:,ind_hei] = jcspectra[ip,:,ind_hei] - jcspc_interf
686
687 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
688 maxid = ListAux.index(max(ListAux))
689
690 ind = numpy.array([-2,-1,1,2])
691 xx = numpy.zeros([4,4])
692
693 for id1 in range(4):
694 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
695
696 xx_inv = numpy.linalg.inv(xx)
697 xx = xx_inv[:,0]
698
699 ind = (ind + maxid + num_mask_prof)%num_mask_prof
700 yy = jcspectra[ip,mask_prof[ind],:]
701 jcspectra[ip,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
702
703 #Guardar Resultados
704 self.dataOut.data_spc = jspectra
705 self.dataOut.data_cspc = jcspectra
706
707 return 1
708
709 def setRadarFrequency(self, frequency=None):
710
711 if frequency != None:
712 self.dataOut.frequency = frequency
713
714 return 1
715
716 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
717 #validacion de rango
718 if minHei == None:
719 minHei = self.dataOut.heightList[0]
720
721 if maxHei == None:
722 maxHei = self.dataOut.heightList[-1]
723
724 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
725 print('minHei: %.2f is out of the heights range'%(minHei))
726 print('minHei is setting to %.2f'%(self.dataOut.heightList[0]))
727 minHei = self.dataOut.heightList[0]
728
729 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
730 print('maxHei: %.2f is out of the heights range'%(maxHei))
731 print('maxHei is setting to %.2f'%(self.dataOut.heightList[-1]))
732 maxHei = self.dataOut.heightList[-1]
733
734 # validacion de velocidades
735 velrange = self.dataOut.getVelRange(1)
736
737 if minVel == None:
738 minVel = velrange[0]
739
740 if maxVel == None:
741 maxVel = velrange[-1]
742
743 if (minVel < velrange[0]) or (minVel > maxVel):
744 print('minVel: %.2f is out of the velocity range'%(minVel))
745 print('minVel is setting to %.2f'%(velrange[0]))
746 minVel = velrange[0]
747
748 if (maxVel > velrange[-1]) or (maxVel < minVel):
749 print('maxVel: %.2f is out of the velocity range'%(maxVel))
750 print('maxVel is setting to %.2f'%(velrange[-1]))
751 maxVel = velrange[-1]
752
753 # seleccion de indices para rango
754 minIndex = 0
755 maxIndex = 0
756 heights = self.dataOut.heightList
757
758 inda = numpy.where(heights >= minHei)
759 indb = numpy.where(heights <= maxHei)
760
761 try:
762 minIndex = inda[0][0]
763 except:
764 minIndex = 0
765
766 try:
767 maxIndex = indb[0][-1]
768 except:
769 maxIndex = len(heights)
770
771 if (minIndex < 0) or (minIndex > maxIndex):
772 raise ValueError("some value in (%d,%d) is not valid" % (minIndex, maxIndex))
773
774 if (maxIndex >= self.dataOut.nHeights):
775 maxIndex = self.dataOut.nHeights-1
776
777 # seleccion de indices para velocidades
778 indminvel = numpy.where(velrange >= minVel)
779 indmaxvel = numpy.where(velrange <= maxVel)
780 try:
781 minIndexVel = indminvel[0][0]
782 except:
783 minIndexVel = 0
784
785 try:
786 maxIndexVel = indmaxvel[0][-1]
787 except:
788 maxIndexVel = len(velrange)
789
790 #seleccion del espectro
791 data_spc = self.dataOut.data_spc[:,minIndexVel:maxIndexVel+1,minIndex:maxIndex+1]
792 #estimacion de ruido
793 noise = numpy.zeros(self.dataOut.nChannels)
794
795 for channel in range(self.dataOut.nChannels):
796 daux = data_spc[channel,:,:]
797 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
798
799 self.dataOut.noise_estimation = noise.copy()
800
801 return 1
802
803 class SpectraAFCProc(ProcessingUnit):
804
805 def __init__(self, **kwargs):
806
807 ProcessingUnit.__init__(self, **kwargs)
808
809 self.buffer = None
810 self.firstdatatime = None
811 self.profIndex = 0
812 self.dataOut = Spectra()
813 self.id_min = None
814 self.id_max = None
815
816 def __updateSpecFromVoltage(self):
817
818 self.dataOut.plotting = "spectra_acf"
819
820 self.dataOut.timeZone = self.dataIn.timeZone
821 self.dataOut.dstFlag = self.dataIn.dstFlag
822 self.dataOut.errorCount = self.dataIn.errorCount
823 self.dataOut.useLocalTime = self.dataIn.useLocalTime
824
825 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
826 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
827 self.dataOut.ippSeconds = self.dataIn.getDeltaH()*(10**-6)/0.15
828
829 self.dataOut.channelList = self.dataIn.channelList
830 self.dataOut.heightList = self.dataIn.heightList
831 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
832
833 self.dataOut.nBaud = self.dataIn.nBaud
834 self.dataOut.nCode = self.dataIn.nCode
835 self.dataOut.code = self.dataIn.code
836 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
837
838 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
839 self.dataOut.utctime = self.firstdatatime
840 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
841 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
842 self.dataOut.flagShiftFFT = False
843
844 self.dataOut.nCohInt = self.dataIn.nCohInt
845 self.dataOut.nIncohInt = 1
846
847 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
848
849 self.dataOut.frequency = self.dataIn.frequency
850 self.dataOut.realtime = self.dataIn.realtime
851
852 self.dataOut.azimuth = self.dataIn.azimuth
853 self.dataOut.zenith = self.dataIn.zenith
854
855 self.dataOut.beam.codeList = self.dataIn.beam.codeList
856 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
857 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
858
859 def __decodeData(self, nProfiles, code):
860
861 if code is None:
862 return
863
864 for i in range(nProfiles):
865 self.buffer[:,i,:] = self.buffer[:,i,:]*code[0][i]
866
867 def __getFft(self):
868 """
869 Convierte valores de Voltaje a Spectra
870
871 Affected:
872 self.dataOut.data_spc
873 self.dataOut.data_cspc
874 self.dataOut.data_dc
875 self.dataOut.heightList
876 self.profIndex
877 self.buffer
878 self.dataOut.flagNoData
879 """
880 nsegments = self.dataOut.nHeights
881
882 _fft_buffer = numpy.zeros((self.dataOut.nChannels, self.dataOut.nProfiles, nsegments), dtype='complex')
883
884 for i in range(nsegments):
885 try:
886 _fft_buffer[:,:,i] = self.buffer[:,i:i+self.dataOut.nProfiles]
887
888 if self.code is not None:
889 _fft_buffer[:,:,i] = _fft_buffer[:,:,i]*self.code[0]
890 except:
891 pass
892
893 fft_volt = numpy.fft.fft(_fft_buffer, n=self.dataOut.nFFTPoints, axis=1)
894 fft_volt = fft_volt.astype(numpy.dtype('complex'))
895 dc = fft_volt[:,0,:]
896
897 #calculo de self-spectra
898 spc = fft_volt * numpy.conjugate(fft_volt)
899 data = numpy.fft.ifft(spc, axis=1)
900 data = numpy.fft.fftshift(data, axes=(1,))
901
902 spc = data
903
904 blocksize = 0
905 blocksize += dc.size
906 blocksize += spc.size
907
908 cspc = None
909 pairIndex = 0
910
911 if self.dataOut.pairsList != None:
912 #calculo de cross-spectra
913 cspc = numpy.zeros((self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
914 for pair in self.dataOut.pairsList:
915 if pair[0] not in self.dataOut.channelList:
916 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList)))
917 if pair[1] not in self.dataOut.channelList:
918 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList)))
919
920 chan_index0 = self.dataOut.channelList.index(pair[0])
921 chan_index1 = self.dataOut.channelList.index(pair[1])
922
923 cspc[pairIndex,:,:] = fft_volt[chan_index0,:,:] * numpy.conjugate(fft_volt[chan_index1,:,:])
924 pairIndex += 1
925 blocksize += cspc.size
926
927 self.dataOut.data_spc = spc
928 self.dataOut.data_cspc = cspc
929 self.dataOut.data_dc = dc
930 self.dataOut.blockSize = blocksize
931 self.dataOut.flagShiftFFT = True
932
933 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], code=None, nCode=1, nBaud=1):
934
935 #self.dataOut.flagNoData = True
936
937 #self.dataIn.runNextUnit = runNextUnit
938 #print("here")
939 if self.dataIn.type == "Spectra":
940 self.dataOut.copy(self.dataIn)
941
942 spc = self.dataOut.data_spc
943 data = numpy.fft.fftshift( spc, axes=(1,))
944 data = numpy.fft.ifft(data, axis=1)
945 #data = numpy.fft.ifft(data, axis=1, n = 32)
946 #data = numpy.fft.fftshift( data, axes=(1,))
947 #acf = numpy.abs(data)
948 acf = data[:,:16,:]
949 #acf = data[:,16:,:]
950 #print("SUM: ",numpy.sum(acf))
951 #print(acf.shape)
952 '''
953 hei_id = 35
954
955 aux2 = numpy.fft.fft(self.dataOut.data[0,:,hei_id],n = 16)
956 aux2 = numpy.fft.fftshift(aux2)
957 aux2 = aux2*numpy.conjugate(aux2)
958 aux2 = aux2.real #Este valor es el que da SCh
959 print("spc 2: ",numpy.sum(aux2))
960 aux2 = numpy.fft.fftshift(aux2)
961 aux2 = numpy.fft.ifft(aux2)
962 print("Rate_Right?: ",aux2[0]/corr[0,0,hei_id])
963
964 print("AFC sum: ",numpy.sum(acf[0,:,hei_id]))
965 print("aux2 sum: ",numpy.sum(aux2))
966 print("Rate: ",acf[0,0,hei_id]/corr[0,0,hei_id])
967 print("Rate aux: ",acf[0,0,hei_id]/corr_aux[0,-1,hei_id])
968 '''
969 #print(acf[0,:,150])
970 #exit(1)
971 '''
972 if real:
973 acf = data.real
974 if imag:
975 acf = data.imag
976 '''
977 #shape = acf.shape # nchannels, nprofiles, nsamples //nchannles, lags , alturas
978 '''
979 for j in range(shape[0]):
980 for i in range(shape[2]):
981 tmp = int(shape[1]/2)
982 #print(i,j,tmp)
983 value = (acf[j,:,i][tmp-1]+acf[j,:,i][tmp+1])/2.0
984 acf[j,:,i][tmp] = value
985 '''
986 #acf = numpy.fft.fftshift( acf, axes=(1,))
987 '''
988 # Normalizando
989 for i in range(shape[0]):
990 for j in range(shape[2]):
991 acf[i,:,j]= acf[i,:,j] / numpy.max(numpy.abs(acf[i,:,j]))
992 '''
993
994 #self.dataOut.data_acf = acf[:,:16,:]#*2
995 #self.dataOut.data_spc = acf[:,:16,:].real#*2
996
997 self.dataOut.data_acf = acf
998 '''
999 self.dataOut.data_spc = data.imag
1000
1001 print("Real: ",data[0,:,26].real)
1002 print("Real dB:",10*numpy.log10(data[0,:,26].real))
1003 print("Imag: ",data[0,:,26].imag)
1004 print("Imag dB:",10*numpy.log10(data[0,:,26].imag))
1005 exit(1)
1006 '''
1007 #print("AFC",self.dataOut.flagNoData)
1008 #'''
1009 #print("acf 0: ", self.dataOut.data_acf[0,0,100])
1010 #print("spc: ",numpy.mean(self.dataOut.data_spc[0,:,100]))
1011 #print("spc 0: ",numpy.fft.fftshift(self.dataOut.data_spc[0,:,100])[0])
1012 #exit(1)
1013 #self.dataOut.data_spc = acf.real
1014 #print(self.dataOut.data_acf[0,:,0])
1015 #exit(1)
1016 #'''
1017 '''
1018 shape = self.dataOut.data_acf.shape
1019 resFactor = 5
1020 z = self.dataOut.data_acf.copy()
1021 min = numpy.min(z[0,:,0])
1022 max =numpy.max(z[0,:,0])
1023 deltaHeight = self.dataOut.heightList[1]-self.dataOut.heightList[0]
1024 for i in range(shape[0]):
1025 for j in range(shape[2]):
1026 z[i,:,j]= (((z[i,:,j]-min)/(max-min))*deltaHeight*resFactor + j*deltaHeight)
1027 #print(self.dataOut.data_spc.shape)
1028 #print(self.dataOut.data_acf.shape)
1029 '''
1030 '''
1031 import matplotlib.pyplot as plt
1032 #hei = 10
1033 #print(self.dataOut.heightList)
1034 #print(self.dataOut.heightList[hei])
1035 #plt.plot(z[0,0,:],self.dataOut.heightList)
1036 aux = self.dataOut.data_acf[0,:,100]
1037 power = aux*numpy.conjugate(aux)
1038 #print(power)
1039 powerdb = numpy.log10(power)
1040 #plt.plot(powerdb,self.dataOut.heightList)
1041 #plt.plot(self.dataOut.data_acf[0,:,33])
1042 #plt.plot(corr)
1043 #plt.imshow(corr.real[0])
1044 plt.imshow(self.dataOut.data_acf.real[0])
1045 #plt.ylim(0,1000)
1046 plt.grid()
1047 plt.show()
1048 exit(1)
1049 '''
201 return True
1050 return True
202
1051
203 if code is not None:
1052 if code is not None:
Show file before | Show file after | Hide comments
@@ -17,6 +17,7 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operat
17 from schainpy.model.data.jrodata import Spectra
17 from schainpy.model.data.jrodata import Spectra
18 from schainpy.model.data.jrodata import hildebrand_sekhon
18 from schainpy.model.data.jrodata import hildebrand_sekhon
19 from schainpy.utils import log
19 from schainpy.utils import log
20 from schainpy.model.data import _HS_algorithm
20
21
21 from time import time, mktime, strptime, gmtime, ctime
22 from time import time, mktime, strptime, gmtime, ctime
22
23
@@ -87,6 +88,7 class SpectraLagProc(ProcessingUnit):
87 """
88 """
88 #print(self.buffer[1,:,0])
89 #print(self.buffer[1,:,0])
89 #exit(1)
90 #exit(1)
91 #print("buffer shape",self.buffer.shape)
90 fft_volt = numpy.fft.fft(
92 fft_volt = numpy.fft.fft(
91 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
93 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
92 fft_volt = fft_volt.astype(numpy.dtype('complex'))
94 fft_volt = fft_volt.astype(numpy.dtype('complex'))
@@ -208,6 +210,17 class SpectraLagProc(ProcessingUnit):
208 self.dataOut.ByLags=ByLags
210 self.dataOut.ByLags=ByLags
209 self.dataOut.LagPlot=LagPlot
211 self.dataOut.LagPlot=LagPlot
210
212
213 #print(self.dataIn.data.shape)
214 '''
215 try:
216 print(self.dataIn.data.shape)
217 except:
218 print("datalags",self.dataIn.datalags.shape)
219 try:
220 print("datalags",self.dataIn.datalags.shape)
221 except:
222 pass
223 '''
211 if self.dataIn.type == "Spectra":
224 if self.dataIn.type == "Spectra":
212 self.dataOut.copy(self.dataIn)
225 self.dataOut.copy(self.dataIn)
213 if shift_fft:
226 if shift_fft:
@@ -224,6 +237,7 class SpectraLagProc(ProcessingUnit):
224 elif self.dataIn.type == "Voltage":
237 elif self.dataIn.type == "Voltage":
225
238
226 if not self.dataOut.ByLags:
239 if not self.dataOut.ByLags:
240 #self.dataOut.data = self.dataIn.data
227 self.VoltageType(nFFTPoints,nProfiles,ippFactor,pairsList)
241 self.VoltageType(nFFTPoints,nProfiles,ippFactor,pairsList)
228 else:
242 else:
229 self.dataOut.nLags = nLags
243 self.dataOut.nLags = nLags
@@ -1232,7 +1246,9 class IntegrationFaradaySpectra(Operation):
1232 if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
1246 if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
1233 continue
1247 continue
1234 buffer=buffer1[:,j]
1248 buffer=buffer1[:,j]
1235 index,sortID=self.hildebrand_sekhon_Integration(buffer,1)
1249 #index,sortID=self.hildebrand_sekhon_Integration(buffer,1)
1250 index=int(_HS_algorithm.HS_algorithm(numpy.sort(buffer, axis=None),1))
1251 sortID = buffer.argsort()
1236 '''
1252 '''
1237 if i==1 and l==0 and k==37:
1253 if i==1 and l==0 and k==37:
1238 print("j",j)
1254 print("j",j)
@@ -1632,7 +1648,9 class IntegrationFaradaySpectra2(Operation):
1632 print(buffer)
1648 print(buffer)
1633 exit(1)
1649 exit(1)
1634 '''
1650 '''
1635 index,sortID=self.hildebrand_sekhon_Integration(buffer,1)
1651 #index,sortID=self.hildebrand_sekhon_Integration(buffer,1)
1652 index=int(_HS_algorithm.HS_algorithm(numpy.sort(buffer, axis=None),1))
1653 sortID = buffer.argsort()
1636
1654
1637 indexes.append(index)
1655 indexes.append(index)
1638 #sortIDs.append(sortID)
1656 #sortIDs.append(sortID)
@@ -1771,7 +1789,7 class IntegrationFaradaySpectra2(Operation):
1771 exit(1)
1789 exit(1)
1772 '''
1790 '''
1773
1791
1774 print(self.nLags)
1792 #print(self.nLags)
1775 '''
1793 '''
1776 if self.nLags == 16:
1794 if self.nLags == 16:
1777 self.nLags = 3
1795 self.nLags = 3
@@ -1821,7 +1839,9 class IntegrationFaradaySpectra2(Operation):
1821 print(buffer)
1839 print(buffer)
1822 exit(1)
1840 exit(1)
1823 '''
1841 '''
1824 index,sortID=self.hildebrand_sekhon_Integration(numpy.abs(buffer),1)
1842 #index,sortID=self.hildebrand_sekhon_Integration(numpy.abs(buffer),1)
1843 index=int(_HS_algorithm.HS_algorithm(numpy.sort(buffer, axis=None),1))
1844 sortID = buffer.argsort()
1825
1845
1826 indexes.append(index)
1846 indexes.append(index)
1827 #sortIDs.append(sortID)
1847 #sortIDs.append(sortID)
@@ -2251,7 +2271,9 class IntegrationFaradaySpectra3(Operation): #This class should manage data with
2251 print(buffer)
2271 print(buffer)
2252 exit(1)
2272 exit(1)
2253 '''
2273 '''
2254 index,sortID=self.hildebrand_sekhon_Integration(buffer,1)
2274 #index,sortID=self.hildebrand_sekhon_Integration(buffer,1)
2275 index=int(_HS_algorithm.HS_algorithm(numpy.sort(buffer, axis=None),1))
2276 sortID = buffer.argsort()
2255
2277
2256 indexes.append(index)
2278 indexes.append(index)
2257 #sortIDs.append(sortID)
2279 #sortIDs.append(sortID)
@@ -2440,7 +2462,9 class IntegrationFaradaySpectra3(Operation): #This class should manage data with
2440 print(buffer)
2462 print(buffer)
2441 exit(1)
2463 exit(1)
2442 '''
2464 '''
2443 index,sortID=self.hildebrand_sekhon_Integration(numpy.abs(buffer),1)
2465 #index,sortID=self.hildebrand_sekhon_Integration(numpy.abs(buffer),1)
2466 index=int(_HS_algorithm.HS_algorithm(numpy.sort(buffer, axis=None),1))
2467 sortID = buffer.argsort()
2444
2468
2445 indexes.append(index)
2469 indexes.append(index)
2446 #sortIDs.append(sortID)
2470 #sortIDs.append(sortID)
@@ -2572,7 +2596,9 class IntegrationFaradaySpectra3(Operation): #This class should manage data with
2572 #buffer=buffer1[:,j]
2596 #buffer=buffer1[:,j]
2573 buffer=(buffer1[:,j])
2597 buffer=(buffer1[:,j])
2574
2598
2575 index,sortID=self.hildebrand_sekhon_Integration(numpy.abs(buffer),1)
2599 #index,sortID=self.hildebrand_sekhon_Integration(numpy.abs(buffer),1)
2600 index=int(_HS_algorithm.HS_algorithm(numpy.sort(buffer, axis=None),1))
2601 sortID = buffer.argsort()
2576
2602
2577 indexes.append(index)
2603 indexes.append(index)
2578 #sortIDs.append(sortID)
2604 #sortIDs.append(sortID)
@@ -2876,7 +2902,7 class IntegrationFaradaySpectraNoLags(Operation):
2876 self.__profIndex += 1
2902 self.__profIndex += 1
2877
2903
2878 return
2904 return
2879
2905 #'''
2880 def hildebrand_sekhon_Integration(self,data,navg):
2906 def hildebrand_sekhon_Integration(self,data,navg):
2881
2907
2882 sortdata = numpy.sort(data, axis=None)
2908 sortdata = numpy.sort(data, axis=None)
@@ -2894,6 +2920,8 class IntegrationFaradaySpectraNoLags(Operation):
2894 sumq += sortdata[j]**2
2920 sumq += sortdata[j]**2
2895 if j > nums_min:
2921 if j > nums_min:
2896 rtest = float(j)/(j-1) + 1.0/navg
2922 rtest = float(j)/(j-1) + 1.0/navg
2923 #print(rtest)
2924 #print(sump)
2897 if ((sumq*j) > (rtest*sump**2)):
2925 if ((sumq*j) > (rtest*sump**2)):
2898 j = j - 1
2926 j = j - 1
2899 sump = sump - sortdata[j]
2927 sump = sump - sortdata[j]
@@ -2903,7 +2931,7 class IntegrationFaradaySpectraNoLags(Operation):
2903 #lnoise = sump / j
2931 #lnoise = sump / j
2904
2932
2905 return j,sortID
2933 return j,sortID
2906
2934 #'''
2907 def pushData(self):
2935 def pushData(self):
2908 """
2936 """
2909 Return the sum of the last profiles and the profiles used in the sum.
2937 Return the sum of the last profiles and the profiles used in the sum.
@@ -2938,12 +2966,20 class IntegrationFaradaySpectraNoLags(Operation):
2938 if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
2966 if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
2939 continue
2967 continue
2940 buffer=buffer1[:,j]
2968 buffer=buffer1[:,j]
2941 index,sortID=self.hildebrand_sekhon_Integration(buffer,1)
2969 #if k != 100:
2970 index=int(_HS_algorithm.HS_algorithm(numpy.sort(buffer, axis=None),1))
2971 sortID = buffer.argsort()
2972 #else:
2973 #index,sortID=self.hildebrand_sekhon_Integration(buffer,1)
2974 #if k == 100:
2975 # print(k,index,sortID)
2976 # exit(1)
2942
2977
2943 indexes.append(index)
2978 indexes.append(index)
2944 #sortIDs.append(sortID)
2979 #sortIDs.append(sortID)
2945 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
2980 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
2946
2981 #if k == 100:
2982 # exit(1)
2947 outliers_IDs=numpy.array(outliers_IDs)
2983 outliers_IDs=numpy.array(outliers_IDs)
2948 outliers_IDs=outliers_IDs.ravel()
2984 outliers_IDs=outliers_IDs.ravel()
2949 outliers_IDs=numpy.unique(outliers_IDs)
2985 outliers_IDs=numpy.unique(outliers_IDs)
@@ -3332,7 +3368,9 class IncohIntLag(Operation):
3332 return dataOut
3368 return dataOut
3333
3369
3334 dataOut.flagNoData = True
3370 dataOut.flagNoData = True
3335
3371 #print("incohint")
3372 #print("IncohInt",dataOut.data_spc.shape)
3373 #print("IncohInt",dataOut.data_cspc.shape)
3336 if not self.isConfig:
3374 if not self.isConfig:
3337 self.setup(n, timeInterval, overlapping)
3375 self.setup(n, timeInterval, overlapping)
3338 self.isConfig = True
3376 self.isConfig = True
@@ -3352,7 +3390,7 class IncohIntLag(Operation):
3352 dataOut.dataLag_spc,
3390 dataOut.dataLag_spc,
3353 dataOut.dataLag_cspc,
3391 dataOut.dataLag_cspc,
3354 dataOut.dataLag_dc)
3392 dataOut.dataLag_dc)
3355
3393 #print("Incoh Int: ",self.__profIndex,n)
3356 if self.__dataReady:
3394 if self.__dataReady:
3357
3395
3358 if not dataOut.ByLags:
3396 if not dataOut.ByLags:
@@ -3369,7 +3407,8 class IncohIntLag(Operation):
3369 #print(numpy.sum(dataOut.dataLag_spc[1,:,100,2]))
3407 #print(numpy.sum(dataOut.dataLag_spc[1,:,100,2]))
3370 #exit(1)
3408 #exit(1)
3371
3409
3372 #print("HERE")
3410 #print("INCOH INT DONE")
3411 #exit(1)
3373 '''
3412 '''
3374 print(numpy.sum(dataOut.dataLag_spc[0,:,20,10])/32)
3413 print(numpy.sum(dataOut.dataLag_spc[0,:,20,10])/32)
3375 print(numpy.sum(dataOut.dataLag_spc[1,:,20,10])/32)
3414 print(numpy.sum(dataOut.dataLag_spc[1,:,20,10])/32)
@@ -3580,6 +3619,8 class IncohInt(Operation):
3580 dataOut.nIncohInt *= self.n
3619 dataOut.nIncohInt *= self.n
3581 dataOut.utctime = avgdatatime
3620 dataOut.utctime = avgdatatime
3582 dataOut.flagNoData = False
3621 dataOut.flagNoData = False
3622 #print("Power",numpy.sum(dataOut.data_spc[0,:,20:30],axis=0))
3623 #print("Power",numpy.sum(dataOut.data_spc[0,100:110,:],axis=1))
3583 #exit(1)
3624 #exit(1)
3584 #print(numpy.sum(dataOut.data_spc[0,:,53]*numpy.conjugate(dataOut.data_spc[0,:,53])))
3625 #print(numpy.sum(dataOut.data_spc[0,:,53]*numpy.conjugate(dataOut.data_spc[0,:,53])))
3585 #print(numpy.sum(dataOut.data_spc[1,:,53]*numpy.conjugate(dataOut.data_spc[1,:,53])))
3626 #print(numpy.sum(dataOut.data_spc[1,:,53]*numpy.conjugate(dataOut.data_spc[1,:,53])))
@@ -3773,7 +3814,6 class SpectraDataToFaraday(Operation):
3773 dataOut.lat=-11.95
3814 dataOut.lat=-11.95
3774 dataOut.lon=-76.87
3815 dataOut.lon=-76.87
3775
3816
3776
3777 self.normFactor(dataOut)
3817 self.normFactor(dataOut)
3778
3818
3779 dataOut.NDP=dataOut.nHeights
3819 dataOut.NDP=dataOut.nHeights
@@ -3848,7 +3888,7 class SpectraDataToFaraday(Operation):
3848 #print("Noise dB: ",10*numpy.log10(dataOut.tnoise))
3888 #print("Noise dB: ",10*numpy.log10(dataOut.tnoise))
3849 #exit(1)
3889 #exit(1)
3850 #dataOut.pan=dataOut.tnoise[0]/float(dataOut.nProfiles_LP*dataOut.nIncohInt)
3890 #dataOut.pan=dataOut.tnoise[0]/float(dataOut.nProfiles_LP*dataOut.nIncohInt)
3851
3891 print("done")
3852 return dataOut
3892 return dataOut
3853
3893
3854
3894
@@ -3903,7 +3943,7 class SpectraDataToHybrid(SpectraDataToFaraday):
3903
3943
3904
3944
3905
3945
3906 def ConvertDataLP(self,dataOut):
3946 def ConvertDataLP_V0(self,dataOut):
3907
3947
3908 #print(dataOut.dataLag_spc[:,:,:,1]/dataOut.data_spc)
3948 #print(dataOut.dataLag_spc[:,:,:,1]/dataOut.data_spc)
3909 #exit(1)
3949 #exit(1)
@@ -3922,6 +3962,20 class SpectraDataToHybrid(SpectraDataToFaraday):
3922
3962
3923 #exit(1)
3963 #exit(1)
3924
3964
3965 def ConvertDataLP(self,dataOut):
3966
3967 #print(dataOut.dataLag_spc[:,:,:,1]/dataOut.data_spc)
3968 #exit(1)
3969 normfactor=1.0/(dataOut.nIncohInt_LP*dataOut.nProfiles_LP)#*dataOut.nProfiles
3970
3971 buffer = self.dataLag_spc_LP=dataOut.dataLag_spc_LP
3972 ##self.dataLag_cspc_LP=(dataOut.dataLag_cspc_LP.sum(axis=1))*(1./dataOut.nProfiles_LP)
3973 #self.dataLag_dc=dataOut.dataLag_dc.sum(axis=1)/dataOut.rnint2[0]
3974 #aux=numpy.expand_dims(self.dataLag_spc_LP, axis=2)
3975 #print(aux.shape)
3976 ##buffer = numpy.concatenate((numpy.expand_dims(self.dataLag_spc_LP, axis=2),self.dataLag_cspc_LP),axis=2)
3977 dataOut.output_LP_integrated = numpy.transpose(buffer,(1,2,0))
3978
3925 def normFactor(self,dataOut):
3979 def normFactor(self,dataOut):
3926 dataOut.rnint2=numpy.zeros(dataOut.DPL,'float32')
3980 dataOut.rnint2=numpy.zeros(dataOut.DPL,'float32')
3927 for l in range(dataOut.DPL):
3981 for l in range(dataOut.DPL):
@@ -4011,3 +4065,115 class SpectraDataToHybrid(SpectraDataToFaraday):
4011 #exit(1)
4065 #exit(1)
4012
4066
4013 return dataOut
4067 return dataOut
4068
4069 class SpectraDataToHybrid_V2(SpectraDataToFaraday):
4070 """Operation to use spectra data in Faraday processing.
4071
4072 Parameters:
4073 -----------
4074 nint : int
4075 Number of integrations.
4076
4077 Example
4078 --------
4079
4080 op = proc_unit.addOperation(name='SpectraDataToFaraday', optype='other')
4081
4082 """
4083
4084 def __init__(self, **kwargs):
4085
4086 Operation.__init__(self, **kwargs)
4087
4088 self.dataLag_spc=None
4089 self.dataLag_cspc=None
4090 self.dataLag_dc=None
4091 self.dataLag_spc_LP=None
4092 self.dataLag_cspc_LP=None
4093 self.dataLag_dc_LP=None
4094
4095 def noise(self,dataOut):
4096
4097 dataOut.data_spc = dataOut.dataLag_spc_LP.real
4098 #print(dataOut.dataLag_spc.shape)
4099 #exit(1)
4100 #dataOut.data_spc = dataOut.dataLag_spc[:,:,:,0].real
4101 #print("spc noise shape: ",dataOut.data_spc.shape)
4102 dataOut.tnoise = dataOut.getNoise(ymin_index=100,ymax_index=166)
4103 #print("Noise LP: ",10*numpy.log10(dataOut.tnoise))
4104 #exit(1)
4105 #dataOut.tnoise[0]*=0.995#0.976
4106 #dataOut.tnoise[1]*=0.995
4107 #print(dataOut.nProfiles)
4108 #dataOut.pan=dataOut.tnoise[0]/float(dataOut.nProfiles_LP*dataOut.nIncohInt)
4109 #dataOut.pbn=dataOut.tnoise[1]/float(dataOut.nProfiles_LP*dataOut.nIncohInt)
4110 dataOut.pan=dataOut.tnoise[0]/float(dataOut.nProfiles_LP*dataOut.nIncohInt_LP)
4111 dataOut.pbn=dataOut.tnoise[1]/float(dataOut.nProfiles_LP*dataOut.nIncohInt_LP)
4112 ##dataOut.pan=dataOut.tnoise[0]*float(self.normfactor_LP)
4113 ##dataOut.pbn=dataOut.tnoise[1]*float(self.normfactor_LP)
4114 #print("pan: ",10*numpy.log10(dataOut.pan))
4115 #print("pbn: ",dataOut.pbn)
4116 #print(numpy.shape(dataOut.pnoise))
4117 #exit(1)
4118 #print("pan: ",numpy.sum(dataOut.pan))
4119 #exit(1)
4120
4121 def ConvertDataLP(self,dataOut):
4122
4123 #print(dataOut.dataLag_spc[:,:,:,1]/dataOut.data_spc)
4124 #exit(1)
4125 self.normfactor_LP=1.0/(dataOut.nIncohInt_LP*dataOut.nProfiles_LP)#*dataOut.nProfiles
4126 #print("acf: ",dataOut.data_acf[0,0,100])
4127 #print("Power: ",numpy.mean(dataOut.dataLag_spc_LP[0,:,100]))
4128 #buffer = dataOut.data_acf*(1./(normfactor*dataOut.nProfiles_LP))
4129 #buffer = dataOut.data_acf*(1./(normfactor))
4130 buffer = dataOut.data_acf#*(self.normfactor_LP)
4131 #print("acf: ",numpy.sum(buffer))
4132
4133 dataOut.output_LP_integrated = numpy.transpose(buffer,(1,2,0))
4134
4135 def normFactor(self,dataOut):
4136 dataOut.rnint2=numpy.zeros(dataOut.DPL,'float32')
4137 for l in range(dataOut.DPL):
4138 if(l==0 or (l>=3 and l <=6)):
4139 dataOut.rnint2[l]=1.0/(dataOut.nIncohInt*dataOut.nProfiles)
4140 else:
4141 dataOut.rnint2[l]=2*(1.0/(dataOut.nIncohInt*dataOut.nProfiles))
4142
4143 def run(self,dataOut):
4144
4145 dataOut.paramInterval=0#int(dataOut.nint*dataOut.header[7][0]*2 )
4146 dataOut.lat=-11.95
4147 dataOut.lon=-76.87
4148
4149 dataOut.NDP=dataOut.nHeights
4150 dataOut.NR=len(dataOut.channelList)
4151 dataOut.DH=dataOut.heightList[1]-dataOut.heightList[0]
4152 dataOut.H0=int(dataOut.heightList[0])
4153
4154 self.normFactor(dataOut)
4155
4156 #Probar sin comentar lo siguiente y comentando
4157 #dataOut.data_acf *= 16 #Corrects the zero padding
4158 #dataOut.dataLag_spc_LP *= 16 #Corrects the zero padding
4159 self.ConvertDataLP(dataOut)
4160 #dataOut.dataLag_spc_LP *= 2
4161 #dataOut.output_LP_integrated[:,:,3] *= float(dataOut.NSCAN/22)#(dataOut.nNoiseProfiles) #Corrects the zero padding
4162
4163 dataOut.nis=dataOut.NSCAN*dataOut.nIncohInt_LP*10
4164
4165 self.ConvertData(dataOut)
4166
4167 dataOut.kabxys_integrated[4][:,(1,2,7,8,9,10),0] *= 2 #Corrects the zero padding
4168 dataOut.kabxys_integrated[6][:,(1,2,7,8,9,10),0] *= 2 #Corrects the zero padding
4169 dataOut.kabxys_integrated[8][:,(1,2,7,8,9,10),0] *= 2 #Corrects the zero padding
4170 dataOut.kabxys_integrated[10][:,(1,2,7,8,9,10),0] *= 2 #Corrects the zero padding
4171 hei = 2
4172
4173 self.noise(dataOut)
4174
4175 dataOut.NAVG=1#dataOut.rnint2[0] #CHECK THIS!
4176 dataOut.nint=dataOut.nIncohInt
4177 dataOut.MAXNRANGENDT=dataOut.output_LP_integrated.shape[1]
4178
4179 return dataOut
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (1180 lines changed) Show them Hide them
@@ -442,13 +442,22 class deFlipHP(Operation):
442 if not channelList:
442 if not channelList:
443 data[:,:] = data[:,:]*self.flip
443 data[:,:] = data[:,:]*self.flip
444 else:
444 else:
445 channelList=[1]
445 #channelList=[1]
446
446
447 for thisChannel in channelList:
447 for thisChannel in channelList:
448 if thisChannel not in dataOut.channelList:
448 if thisChannel not in dataOut.channelList:
449 continue
449 continue
450
450
451 data[thisChannel,:] = data[thisChannel,:]*self.flip
451 if not byHeights:
452 data[thisChannel,:] = data[thisChannel,:]*flip
453
454 else:
455 firstHeight = HeiRangeList[0]
456 lastHeight = HeiRangeList[1]+1
457 flip = -1.0
458 data[thisChannel,firstHeight:lastHeight] = data[thisChannel,firstHeight:lastHeight]*flip
459
460 #data[thisChannel,:] = data[thisChannel,:]*self.flip
452
461
453 self.flip *= -1.
462 self.flip *= -1.
454
463
@@ -1093,28 +1102,56 class LagsReshapeDP_V2(Operation):
1093 #print(dataOut.data[1,:12,:15])
1102 #print(dataOut.data[1,:12,:15])
1094 #exit(1)
1103 #exit(1)
1095 #print(numpy.shape(dataOut.data))
1104 #print(numpy.shape(dataOut.data))
1096 print(dataOut.profileIndex)
1105 #print(dataOut.profileIndex)
1097 #dataOut.flagNoData = True
1098 if dataOut.profileIndex == 140:
1099 print("here")
1100 print(dataOut.data.shape)
1101
1102 dataOut.data = numpy.transpose(numpy.array(self.data_buffer),(1,0,2))
1103 print(dataOut.data.shape)
1104 dataOut.datalags = numpy.copy(self.LagDistribution(dataOut))
1105 #print(numpy.shape(dataOut.datalags))
1106 #exit(1)
1107 #print("AFTER RESHAPE DP")
1108
1106
1109 dataOut.data = dataOut.data[:,:,200:]
1107 if not dataOut.flagDataAsBlock:
1110 self.data_buffer = []
1111 #dataOut.flagNoData = False
1112
1108
1113 else:
1109 dataOut.flagNoData = True
1114 self.data_buffer.append(dataOut.data)
1110 #print("nProfiles: ",dataOut.nProfiles)
1115 #print(numpy.shape(dataOut.data))
1111 #if dataOut.profileIndex == 140:
1116 #exit(1)
1112 #print("id: ",dataOut.profileIndex)
1113 if dataOut.profileIndex == dataOut.nProfiles-1:
1114 #print("here")
1115 #print(dataOut.data.shape)
1116 self.data_buffer.append(dataOut.data)
1117 dataOut.data = numpy.transpose(numpy.array(self.data_buffer),(1,0,2))
1118 #print(dataOut.data.shape)
1119 #print(numpy.sum(dataOut.data))
1120 #print(dataOut.data[1,100,:])
1121 #exit(1)
1122 dataOut.datalags = numpy.copy(self.LagDistribution(dataOut))
1123 #print(numpy.shape(dataOut.datalags))
1124 #exit(1)
1125 #print("AFTER RESHAPE DP")
1126
1127 dataOut.data = dataOut.data[:,:,200:]
1128 self.data_buffer = []
1129 dataOut.flagDataAsBlock = True
1130 dataOut.flagNoData = False
1131
1132 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1133 dataOut.heightList = numpy.arange(dataOut.NDP) *deltaHeight# + dataOut.heightList[0]
1134 #exit(1)
1135 #print(numpy.sum(dataOut.datalags))
1136 #exit(1)
1117
1137
1138 else:
1139 self.data_buffer.append(dataOut.data)
1140 #print(numpy.shape(dataOut.data))
1141 #exit(1)
1142 else:
1143 #print(dataOut.data.shape)
1144 #print(numpy.sum(dataOut.data))
1145 #print(dataOut.data[1,100,:])
1146 #exit(1)
1147 dataOut.datalags = numpy.copy(self.LagDistribution(dataOut))
1148 #print(dataOut.datalags.shape)
1149 dataOut.data = dataOut.data[:,:,200:]
1150 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1151 dataOut.heightList = numpy.arange(dataOut.NDP) * deltaHeight# + dataOut.heightList[0]
1152 #print(dataOut.nHeights)
1153 #print(numpy.sum(dataOut.datalags))
1154 #exit(1)
1118
1155
1119 return dataOut
1156 return dataOut
1120
1157
@@ -2523,8 +2560,8 class CleanCohEchoes(Operation):
2523
2560
2524 return dataOut
2561 return dataOut
2525
2562
2526 class NoisePower(Operation):
2563 class CleanCohEchoesHP(Operation):
2527 """Operation to get noise power from the integrated data of the Double Pulse.
2564 """Operation to clean coherent echoes.
2528
2565
2529 Parameters:
2566 Parameters:
2530 -----------
2567 -----------
@@ -2533,7 +2570,7 class NoisePower(Operation):
2533 Example
2570 Example
2534 --------
2571 --------
2535
2572
2536 op = proc_unit.addOperation(name='NoisePower', optype='other')
2573 op = proc_unit.addOperation(name='CleanCohEchoes')
2537
2574
2538 """
2575 """
2539
2576
@@ -2541,125 +2578,498 class NoisePower(Operation):
2541
2578
2542 Operation.__init__(self, **kwargs)
2579 Operation.__init__(self, **kwargs)
2543
2580
2544 def hildebrand(self,dataOut,data):
2581 def remove_coh(self,pow):
2545
2582 #print("pow inside: ",pow)
2546 divider=10 # divider was originally 10
2583 #print(pow.shape)
2547 noise=0.0
2584 q75,q25 = numpy.percentile(pow,[75,25],axis=0)
2548 n1=0
2585 #print(q75,q25)
2549 n2=int(dataOut.NDP/2)
2586 intr_qr = q75-q25
2550 sorts= sorted(data)
2551 nums_min= dataOut.NDP/divider
2552 if((dataOut.NDP/divider)> 2):
2553 nums_min= int(dataOut.NDP/divider)
2554
2555 else:
2556 nums_min=2
2557 sump=0.0
2558 sumq=0.0
2559 j=0
2560 cont=1
2561 while( (cont==1) and (j<n2)):
2562 sump+=sorts[j+n1]
2563 sumq+= sorts[j+n1]*sorts[j+n1]
2564 t3= sump/(j+1)
2565 j=j+1
2566 if(j> nums_min):
2567 rtest= float(j/(j-1)) +1.0/dataOut.NAVG
2568 t1= (sumq*j)
2569 t2=(rtest*sump*sump)
2570 if( (t1/t2) > 0.990):
2571 j=j-1
2572 sump-= sorts[j+n1]
2573 sumq-=sorts[j+n1]*sorts[j+n1]
2574 cont= 0
2575
2587
2576 noise= sump/j
2588 max = q75+(1.5*intr_qr)
2577 stdv=numpy.sqrt((sumq- noise*noise)/(j-1))
2589 min = q25-(1.5*intr_qr)
2578 return noise
2579
2590
2580 def run(self,dataOut):
2591 pow[pow > max] = numpy.nan
2581
2592
2582 p=numpy.zeros((dataOut.NR,dataOut.NDP,dataOut.DPL),'float32')
2593 #print("Max: ",max)
2583 av=numpy.zeros(dataOut.NDP,'float32')
2594 #print("Min: ",min)
2584 dataOut.pnoise=numpy.zeros(dataOut.NR,'float32')
2585
2595
2586 p[0,:,:]=dataOut.kabxys_integrated[4][:,:,0]+dataOut.kabxys_integrated[5][:,:,0] #total power for channel 0, just pulse with non-flip
2596 return pow
2587 p[1,:,:]=dataOut.kabxys_integrated[6][:,:,0]+dataOut.kabxys_integrated[7][:,:,0] #total power for channel 1
2588
2597
2589 for i in range(dataOut.NR):
2598 def mad_based_outlier_V0(self, points, thresh=3.5):
2590 dataOut.pnoise[i]=0.0
2599 #print("points: ",points)
2591 for k in range(dataOut.DPL):
2600 if len(points.shape) == 1:
2592 dataOut.pnoise[i]+= self.hildebrand(dataOut,p[i,:,k])
2601 points = points[:,None]
2602 median = numpy.nanmedian(points, axis=0)
2603 diff = numpy.nansum((points - median)**2, axis=-1)
2604 diff = numpy.sqrt(diff)
2605 med_abs_deviation = numpy.nanmedian(diff)
2593
2606
2594 dataOut.pnoise[i]=dataOut.pnoise[i]/dataOut.DPL
2607 modified_z_score = 0.6745 * diff / med_abs_deviation
2608 #print(modified_z_score)
2609 return modified_z_score > thresh
2595
2610
2611 def mad_based_outlier(self, points, thresh=3.5):
2596
2612
2597 dataOut.pan=1.0*dataOut.pnoise[0] # weights could change
2613 median = numpy.nanmedian(points)
2598 dataOut.pbn=1.0*dataOut.pnoise[1] # weights could change
2614 diff = (points - median)**2
2599 print("pan: ",dataOut.pan)
2615 diff = numpy.sqrt(diff)
2600 print("pbn: ",dataOut.pbn)
2616 med_abs_deviation = numpy.nanmedian(diff)
2601 print("pan dB: ",10*numpy.log10(dataOut.pan))
2602 print("pbn dB: ",10*numpy.log10(dataOut.pbn))
2603 exit(1)
2604 dataOut.power = dataOut.getPower()
2605 return dataOut
2606
2617
2618 modified_z_score = 0.6745 * diff / med_abs_deviation
2607
2619
2608 class DoublePulseACFs(Operation):
2620 return modified_z_score > thresh
2609 """Operation to get the ACFs of the Double Pulse.
2610
2621
2611 Parameters:
2622 def removeSpreadF_V0(self,dataOut):
2612 -----------
2623 for i in range(11):
2613 None
2624 print("BEFORE Chb: ",i,dataOut.kabxys_integrated[6][:,i,0])
2625 #exit(1)
2614
2626
2615 Example
2627 #Removing echoes greater than 35 dB
2616 --------
2628 maxdB = 35 #DEBERÍA SER NOISE+ALGO!!!!!!!!!!!!!!!!!!!!!!
2629 #print(dataOut.kabxys_integrated[6][:,0,0])
2630 data = numpy.copy(10*numpy.log10(dataOut.kabxys_integrated[6][:,0,0])) #Lag0 ChB
2631 #print(data)
2632 for i in range(12,data.shape[0]):
2633 #for j in range(data.shape[1]):
2634 if data[i]>maxdB:
2635 dataOut.kabxys_integrated[4][i-2:i+3,:,0] = numpy.nan #Debido a que estos ecos son intensos, se
2636 dataOut.kabxys_integrated[6][i-2:i+3,:,0] = numpy.nan #remueve además dos muestras antes y después
2637 #dataOut.kabxys_integrated[4][i-1,:,0] = numpy.nan
2638 #dataOut.kabxys_integrated[6][i-1,:,0] = numpy.nan
2639 #dataOut.kabxys_integrated[4][i+1,:,0] = numpy.nan
2640 #dataOut.kabxys_integrated[6][i+1,:,0] = numpy.nan
2641 dataOut.flagSpreadF = True
2642 print("Removing Threshold",i)
2643 #print("i: ",i)
2617
2644
2618 op = proc_unit.addOperation(name='DoublePulseACFs', optype='other')
2645 #print("BEFORE Chb: ",dataOut.kabxys_integrated[6][:,0,0])
2646 #exit(1)
2619
2647
2620 """
2648 #Removing outliers from the profile
2649 nlag = 9
2650 minHei = 180
2651 #maxHei = 600
2652 maxHei = 525
2653 inda = numpy.where(dataOut.heightList >= minHei)
2654 indb = numpy.where(dataOut.heightList <= maxHei)
2655 minIndex = inda[0][0]
2656 maxIndex = indb[0][-1]
2657 #l0 = 0
2658 #print("BEFORE Cha: ",dataOut.kabxys_integrated[4][:,l0,0])
2659 #print("BEFORE Chb: ",dataOut.kabxys_integrated[6][:,l0,0])
2660 #exit(1)
2661 #'''
2662 l0 = 0
2663 #print("BEFORE Cha: ",dataOut.kabxys_integrated[4][:,l0,0])
2664 #print("BEFORE Chb: ",dataOut.kabxys_integrated[6][:,l0,0])
2621
2665
2622 def __init__(self, **kwargs):
2666 import matplotlib.pyplot as plt
2667 for i in range(l0,l0+11):
2668 plt.plot(dataOut.kabxys_integrated[6][:,i,0],dataOut.heightList,label='{}'.format(i))
2669 #plt.xlim(1.e5,1.e8)
2670 plt.legend()
2671 plt.xlim(0,2000)
2672 plt.show()
2673 #'''
2674 #dataOut.kabxys_integrated[4][minIndex:,:,0] = self.remove_coh(dataOut.kabxys_integrated[4][minIndex:,:,0
2675 outliers_IDs = []
2676 '''
2677 for lag in range(11):
2678 outliers = self.mad_based_outlier(dataOut.kabxys_integrated[4][minIndex:,lag,0], thresh=3.)
2679 #print("Outliers: ",outliers)
2680 #indexes.append(outliers.nonzero())
2681 #numpy.concatenate((outliers))
2682 #dataOut.kabxys_integrated[4][minIndex:,lag,0][outliers == True] = numpy.nan
2683 outliers_IDs=numpy.append(outliers_IDs,outliers.nonzero())
2684 '''
2685 for lag in range(11):
2686 #outliers = self.mad_based_outlier(dataOut.kabxys_integrated[6][minIndex:maxIndex,lag,0], thresh=2.)
2687 outliers = self.mad_based_outlier(dataOut.kabxys_integrated[6][minIndex:maxIndex,lag,0])
2688 outliers_IDs=numpy.append(outliers_IDs,outliers.nonzero())
2689 #print(outliers_IDs)
2690 #exit(1)
2691 if outliers_IDs != []:
2692 outliers_IDs=numpy.array(outliers_IDs)
2693 outliers_IDs=outliers_IDs.ravel()
2694 outliers_IDs=outliers_IDs.astype(numpy.dtype('int64'))
2623
2695
2624 Operation.__init__(self, **kwargs)
2696 (uniq, freq) = (numpy.unique(outliers_IDs, return_counts=True))
2625 self.aux=1
2697 aux_arr = numpy.column_stack((uniq,freq))
2698 #print("repetitions: ",aux_arr)
2626
2699
2627 def run(self,dataOut):
2700 #if aux_arr != []:
2701 final_index = []
2702 for i in range(aux_arr.shape[0]):
2703 if aux_arr[i,1] >= 10:
2704 final_index.append(aux_arr[i,0])
2628
2705
2629 dataOut.igcej=numpy.zeros((dataOut.NDP,dataOut.DPL),'int32')
2706 if final_index != [] and len(final_index) > 1:
2707 final_index += minIndex
2708 #print("final_index: ",final_index)
2709 following_index = final_index[-1]+1 #Remove following index to ensure we remove remaining SpreadF
2710 previous_index = final_index[0]-1 #Remove previous index to ensure we remove remaning SpreadF
2711 final_index = numpy.concatenate(([previous_index],final_index,[following_index]))
2712 final_index = numpy.unique(final_index) #If there was only one outlier
2713 #print("final_index: ",final_index)
2714 #exit(1)
2715 dataOut.kabxys_integrated[4][final_index,:,0] = numpy.nan
2716 dataOut.kabxys_integrated[6][final_index,:,0] = numpy.nan
2630
2717
2631 if self.aux==1:
2718 dataOut.flagSpreadF = True
2632 dataOut.rhor=numpy.zeros((dataOut.NDP,dataOut.DPL), dtype=float)
2633 dataOut.rhoi=numpy.zeros((dataOut.NDP,dataOut.DPL), dtype=float)
2634 dataOut.sdp=numpy.zeros((dataOut.NDP,dataOut.DPL), dtype=float)
2635 dataOut.sd=numpy.zeros((dataOut.NDP,dataOut.DPL), dtype=float)
2636 dataOut.p=numpy.zeros((dataOut.NDP,dataOut.DPL), dtype=float)
2637 dataOut.alag=numpy.zeros(dataOut.NDP,'float32')
2638 for l in range(dataOut.DPL):
2639 dataOut.alag[l]=l*dataOut.DH*2.0/150.0
2640 self.aux=0
2641 sn4=dataOut.pan*dataOut.pbn
2642 rhorn=0
2643 rhoin=0
2644 panrm=numpy.zeros((dataOut.NDP,dataOut.DPL), dtype=float)
2645
2719
2646 id = numpy.where(dataOut.heightList>700)[0]
2720 #print(final_index+minIndex)
2721 #print(outliers_IDs)
2722 #exit(1)
2723 #print("flagSpreadF",dataOut.flagSpreadF)
2647
2724
2648 for i in range(dataOut.NDP):
2725 '''
2649 for j in range(dataOut.DPL):
2726 for lag in range(11):
2650 ################# Total power
2727 #print("Lag: ",lag)
2651 pa=numpy.abs(dataOut.kabxys_integrated[4][i,j,0]+dataOut.kabxys_integrated[5][i,j,0])
2728 outliers = self.mad_based_outlier(dataOut.kabxys_integrated[6][minIndex:maxIndex,lag,0], thresh=2.)
2652 pb=numpy.abs(dataOut.kabxys_integrated[6][i,j,0]+dataOut.kabxys_integrated[7][i,j,0])
2729 dataOut.kabxys_integrated[6][minIndex:maxIndex,lag,0][outliers == True] = numpy.nan
2653 st4=pa*pb
2730 '''
2654 '''
2731 #dataOut.kabxys_integrated[4][minIndex:,:,0] = self.remove_coh(dataOut.kabxys_integrated[4][minIndex:,:,0])
2655 if i > id[0]:
2732 '''
2656 dataOut.p[i,j] pa-dataOut.pan
2733 import matplotlib.pyplot as plt
2657 else:
2734 for i in range(11):
2658 dataOut.p[i,j]=pa+pb-(dataOut.pan+dataOut.pbn)
2735 plt.plot(dataOut.kabxys_integrated[6][:,i,0],dataOut.heightList,label='{}'.format(i))
2659 '''
2736 plt.xlim(0,2000)
2660 dataOut.p[i,j]=pa+pb-(dataOut.pan+dataOut.pbn)
2737 plt.legend()
2661 dataOut.sdp[i,j]=2*dataOut.rnint2[j]*((pa+pb)*(pa+pb))
2738 plt.grid()
2739 plt.show()
2740 '''
2741 '''
2742 for nlag in range(11):
2743 print("BEFORE",dataOut.kabxys_integrated[6][:,nlag,0])
2744 #exit(1)
2745 '''
2746 #dataOut.kabxys_integrated[6][minIndex:,:,0] = self.remove_coh(dataOut.kabxys_integrated[6][minIndex:,:,0])
2747
2748
2749 '''
2750 for nlag in range(11):
2751 print("AFTER",dataOut.kabxys_integrated[6][:,nlag,0])
2752 exit(1)
2753 '''
2754 #print("AFTER",dataOut.kabxys_integrated[4][33,:,0])
2755 #print("AFTER",dataOut.kabxys_integrated[6][33,:,0])
2756 #exit(1)
2757
2758 def removeSpreadF(self,dataOut):
2759 #for i in range(11):
2760 #print("BEFORE Chb: ",i,dataOut.kabxys_integrated[6][:,i,0])
2761 #exit(1)
2762
2763 #for i in range(12,data.shape[0]):
2764 #for j in range(data.shape[1]):
2765 #if data[i]>maxdB:
2766 #dataOut.kabxys_integrated[4][i-2:i+3,:,0] = numpy.nan #Debido a que estos ecos son intensos, se
2767 #dataOut.kabxys_integrated[6][i-2:i+3,:,0] = numpy.nan #remueven además dos muestras antes y después
2768 #dataOut.flagSpreadF = True
2769 #print("Removing Threshold",i)
2770 #print("i: ",i)
2771
2772 #print("BEFORE Chb: ",dataOut.kabxys_integrated[6][:,0,0])
2773 #exit(1)
2774
2775 #Removing outliers from the profile
2776 nlag = 9
2777 minHei = 180
2778 #maxHei = 600
2779 maxHei = 525
2780 inda = numpy.where(dataOut.heightList >= minHei)
2781 indb = numpy.where(dataOut.heightList <= maxHei)
2782 minIndex = inda[0][0]
2783 maxIndex = indb[0][-1]
2784 #l0 = 0
2785 #print("BEFORE Cha: ",dataOut.kabxys_integrated[4][:,l0,0])
2786 #print("BEFORE Chb: ",dataOut.kabxys_integrated[6][:,l0,0])
2787 #exit(1)
2788 '''
2789 l0 = 0
2790 #print("BEFORE Cha: ",dataOut.kabxys_integrated[4][:,l0,0])
2791 #print("BEFORE Chb: ",dataOut.kabxys_integrated[6][:,l0,0])
2792
2793 import matplotlib.pyplot as plt
2794 for i in range(l0,l0+11):
2795 plt.plot(dataOut.kabxys_integrated[6][:,i,0],dataOut.heightList,label='{}'.format(i))
2796 #plt.xlim(1.e5,1.e8)
2797 plt.legend()
2798 plt.xlim(0,2000)
2799 plt.show()
2800 '''
2801 #dataOut.kabxys_integrated[4][minIndex:,:,0] = self.remove_coh(dataOut.kabxys_integrated[4][minIndex:,:,0
2802 outliers_IDs = []
2803 '''
2804 for lag in range(11):
2805 outliers = self.mad_based_outlier(dataOut.kabxys_integrated[4][minIndex:,lag,0], thresh=3.)
2806 #print("Outliers: ",outliers)
2807 #indexes.append(outliers.nonzero())
2808 #numpy.concatenate((outliers))
2809 #dataOut.kabxys_integrated[4][minIndex:,lag,0][outliers == True] = numpy.nan
2810 outliers_IDs=numpy.append(outliers_IDs,outliers.nonzero())
2811 '''
2812 '''
2813 for lag in range(11):
2814 #outliers = self.mad_based_outlier(dataOut.kabxys_integrated[6][minIndex:maxIndex,lag,0], thresh=2.)
2815 outliers = self.mad_based_outlier(dataOut.kabxys_integrated[6][minIndex:maxIndex,lag,0])
2816 outliers_IDs=numpy.append(outliers_IDs,outliers.nonzero())
2817 '''
2818
2819 for i in range(15):
2820 minIndex = 12+i#12
2821 #maxIndex = 22+i#35
2822 if gmtime(dataOut.utctime).tm_hour >= 23. or gmtime(dataOut.utctime).tm_hour < 3.:
2823 maxIndex = 31+i#35
2824 else:
2825 maxIndex = 22+i#35
2826 for lag in range(11):
2827 #outliers = mad_based_outlier(pow_clean3[12:27], thresh=2.)
2828 #print("Cuts: ",first_cut*15, last_cut*15)
2829 outliers = self.mad_based_outlier(dataOut.kabxys_integrated[6][minIndex:maxIndex,lag,0])
2830 aux = minIndex+numpy.array(outliers.nonzero()).ravel()
2831 outliers_IDs=numpy.append(outliers_IDs,aux)
2832 #print(minIndex+numpy.array(outliers.nonzero()).ravel())
2833 #print(outliers_IDs)
2834 #exit(1)
2835 if outliers_IDs != []:
2836 outliers_IDs=numpy.array(outliers_IDs)
2837 #outliers_IDs=outliers_IDs.ravel()
2838 outliers_IDs=outliers_IDs.astype(numpy.dtype('int64'))
2839 #print(outliers_IDs)
2840 #exit(1)
2841
2842 (uniq, freq) = (numpy.unique(outliers_IDs, return_counts=True))
2843 aux_arr = numpy.column_stack((uniq,freq))
2844 #print("repetitions: ",aux_arr)
2845 #exit(1)
2846
2847 #if aux_arr != []:
2848 final_index = []
2849 for i in range(aux_arr.shape[0]):
2850 if aux_arr[i,1] >= 3*11:
2851 final_index.append(aux_arr[i,0])
2852
2853 if final_index != []:# and len(final_index) > 1:
2854 #final_index += minIndex
2855 #print("final_index: ",final_index)
2856 following_index = final_index[-1]+1 #Remove following index to ensure we remove remaining SpreadF
2857 previous_index = final_index[0]-1 #Remove previous index to ensure we remove remaning SpreadF
2858 final_index = numpy.concatenate(([previous_index],final_index,[following_index]))
2859 final_index = numpy.unique(final_index) #If there was only one outlier
2860 #print("final_index: ",final_index)
2861 #exit(1)
2862 dataOut.kabxys_integrated[4][final_index,:,0] = numpy.nan
2863 dataOut.kabxys_integrated[6][final_index,:,0] = numpy.nan
2864
2865 dataOut.flagSpreadF = True
2866
2867 #Removing echoes greater than 35 dB
2868 maxdB = 10*numpy.log10(dataOut.pbn[0]) + 10 #Lag 0 NOise
2869 #maxdB = 35 #DEBERÍA SER NOISE+ALGO!!!!!!!!!!!!!!!!!!!!!!
2870 #print("noise: ",maxdB - 10)
2871 #print(dataOut.kabxys_integrated[6][:,0,0])
2872 data = numpy.copy(10*numpy.log10(dataOut.kabxys_integrated[6][:,0,0])) #Lag0 ChB
2873 #print("data: ",data)
2874
2875 for i in range(12,data.shape[0]):
2876 #for j in range(data.shape[1]):
2877 if data[i]>maxdB:
2878 dataOut.kabxys_integrated[4][i-2:i+3,:,0] = numpy.nan #Debido a que estos ecos son intensos, se
2879 dataOut.kabxys_integrated[6][i-2:i+3,:,0] = numpy.nan #remueven además dos muestras antes y después
2880 dataOut.flagSpreadF = True
2881 #print("Removing Threshold",i)
2882
2883 #print(final_index+minIndex)
2884 #print(outliers_IDs)
2885 #exit(1)
2886 #print("flagSpreadF",dataOut.flagSpreadF)
2887
2888 '''
2889 for lag in range(11):
2890 #print("Lag: ",lag)
2891 outliers = self.mad_based_outlier(dataOut.kabxys_integrated[6][minIndex:maxIndex,lag,0], thresh=2.)
2892 dataOut.kabxys_integrated[6][minIndex:maxIndex,lag,0][outliers == True] = numpy.nan
2893 '''
2894 #dataOut.kabxys_integrated[4][minIndex:,:,0] = self.remove_coh(dataOut.kabxys_integrated[4][minIndex:,:,0])
2895 '''
2896 import matplotlib.pyplot as plt
2897 for i in range(11):
2898 plt.plot(dataOut.kabxys_integrated[6][:,i,0],dataOut.heightList,label='{}'.format(i))
2899 plt.xlim(0,2000)
2900 plt.legend()
2901 plt.grid()
2902 plt.show()
2903 '''
2904 '''
2905 for nlag in range(11):
2906 print("BEFORE",dataOut.kabxys_integrated[6][:,nlag,0])
2907 #exit(1)
2908 '''
2909 #dataOut.kabxys_integrated[6][minIndex:,:,0] = self.remove_coh(dataOut.kabxys_integrated[6][minIndex:,:,0])
2910
2911
2912 '''
2913 for nlag in range(11):
2914 print("AFTER",dataOut.kabxys_integrated[6][:,nlag,0])
2915 exit(1)
2916 '''
2917
2918 def run(self,dataOut):
2919 dataOut.flagSpreadF = False
2920 #print(gmtime(dataOut.utctime).tm_hour)
2921 #print(dataOut.ut_Faraday)
2922 #exit(1)
2923 if gmtime(dataOut.utctime).tm_hour >= 23. or gmtime(dataOut.utctime).tm_hour < 11.: #18-06 LT
2924 #print("Inside if we are in SpreadF Time: ",gmtime(dataOut.utctime).tm_hour)
2925 self.removeSpreadF(dataOut)
2926 #exit(1)
2927
2928 return dataOut
2929
2930 class NoisePower(Operation):
2931 """Operation to get noise power from the integrated data of the Double Pulse.
2932
2933 Parameters:
2934 -----------
2935 None
2936
2937 Example
2938 --------
2939
2940 op = proc_unit.addOperation(name='NoisePower', optype='other')
2941
2942 """
2943
2944 def __init__(self, **kwargs):
2945
2946 Operation.__init__(self, **kwargs)
2947
2948 def hildebrand(self,dataOut,data):
2949
2950 divider=10 # divider was originally 10
2951 noise=0.0
2952 n1=0
2953 n2=int(dataOut.NDP/2)
2954 sorts= sorted(data)
2955 nums_min= dataOut.NDP/divider
2956 if((dataOut.NDP/divider)> 2):
2957 nums_min= int(dataOut.NDP/divider)
2958
2959 else:
2960 nums_min=2
2961 sump=0.0
2962 sumq=0.0
2963 j=0
2964 cont=1
2965 while( (cont==1) and (j<n2)):
2966 sump+=sorts[j+n1]
2967 sumq+= sorts[j+n1]*sorts[j+n1]
2968 t3= sump/(j+1)
2969 j=j+1
2970 if(j> nums_min):
2971 rtest= float(j/(j-1)) +1.0/dataOut.NAVG
2972 t1= (sumq*j)
2973 t2=(rtest*sump*sump)
2974 if( (t1/t2) > 0.990):
2975 j=j-1
2976 sump-= sorts[j+n1]
2977 sumq-=sorts[j+n1]*sorts[j+n1]
2978 cont= 0
2979
2980 noise= sump/j
2981 stdv=numpy.sqrt((sumq- noise*noise)/(j-1))
2982 return noise
2983
2984 def run(self,dataOut):
2985
2986 p=numpy.zeros((dataOut.NR,dataOut.NDP,dataOut.DPL),'float32')
2987 av=numpy.zeros(dataOut.NDP,'float32')
2988 dataOut.pnoise=numpy.zeros(dataOut.NR,'float32')
2989
2990 p[0,:,:]=dataOut.kabxys_integrated[4][:,:,0]+dataOut.kabxys_integrated[5][:,:,0] #total power for channel 0, just pulse with non-flip
2991 p[1,:,:]=dataOut.kabxys_integrated[6][:,:,0]+dataOut.kabxys_integrated[7][:,:,0] #total power for channel 1
2992
2993 for i in range(dataOut.NR):
2994 dataOut.pnoise[i]=0.0
2995 for k in range(dataOut.DPL):
2996 dataOut.pnoise[i]+= self.hildebrand(dataOut,p[i,:,k])
2997
2998 dataOut.pnoise[i]=dataOut.pnoise[i]/dataOut.DPL
2999
3000
3001 dataOut.pan=1.0*dataOut.pnoise[0] # weights could change
3002 dataOut.pbn=1.0*dataOut.pnoise[1] # weights could change
3003 '''
3004 print("pan: ",dataOut.pan)
3005 print("pbn: ",dataOut.pbn)
3006 print("pan dB: ",10*numpy.log10(dataOut.pan))
3007 print("pbn dB: ",10*numpy.log10(dataOut.pbn))
3008 exit(1)
3009 '''
3010 dataOut.power = dataOut.getPower()
3011 return dataOut
3012
3013
3014 class DoublePulseACFs(Operation):
3015 """Operation to get the ACFs of the Double Pulse.
3016
3017 Parameters:
3018 -----------
3019 None
3020
3021 Example
3022 --------
3023
3024 op = proc_unit.addOperation(name='DoublePulseACFs', optype='other')
3025
3026 """
3027
3028 def __init__(self, **kwargs):
3029
3030 Operation.__init__(self, **kwargs)
3031 self.aux=1
3032
3033 def run(self,dataOut):
3034
3035 dataOut.igcej=numpy.zeros((dataOut.NDP,dataOut.DPL),'int32')
3036 #print("init")
3037 if self.aux==1:
3038 dataOut.rhor=numpy.zeros((dataOut.NDP,dataOut.DPL), dtype=float)
3039 dataOut.rhoi=numpy.zeros((dataOut.NDP,dataOut.DPL), dtype=float)
3040 dataOut.sdp=numpy.zeros((dataOut.NDP,dataOut.DPL), dtype=float)
3041 dataOut.sd=numpy.zeros((dataOut.NDP,dataOut.DPL), dtype=float)
3042 dataOut.p=numpy.zeros((dataOut.NDP,dataOut.DPL), dtype=float)
3043 dataOut.alag=numpy.zeros(dataOut.NDP,'float32')
3044 for l in range(dataOut.DPL):
3045 dataOut.alag[l]=l*dataOut.DH*2.0/150.0
3046 self.aux=0
3047 sn4=dataOut.pan*dataOut.pbn
3048 rhorn=0
3049 rhoin=0
3050 panrm=numpy.zeros((dataOut.NDP,dataOut.DPL), dtype=float)
3051
3052 id = numpy.where(dataOut.heightList>700)[0]
3053
3054 for i in range(dataOut.NDP):
3055 for j in range(dataOut.DPL):
3056 ################# Total power
3057 pa=numpy.abs(dataOut.kabxys_integrated[4][i,j,0]+dataOut.kabxys_integrated[5][i,j,0])
3058 pb=numpy.abs(dataOut.kabxys_integrated[6][i,j,0]+dataOut.kabxys_integrated[7][i,j,0])
3059 st4=pa*pb
3060
3061 '''
3062 if i > id[0]:
3063 dataOut.p[i,j] pa-dataOut.pan
3064 else:
3065 dataOut.p[i,j]=pa+pb-(dataOut.pan+dataOut.pbn)
3066 '''
3067 #print("init 2.6",pa,dataOut.pan)
3068 dataOut.p[i,j]=pa+pb-(dataOut.pan+dataOut.pbn)
3069
3070 dataOut.sdp[i,j]=2*dataOut.rnint2[j]*((pa+pb)*(pa+pb))
2662 ## ACF
3071 ## ACF
3072
2663 rhorp=dataOut.kabxys_integrated[8][i,j,0]+dataOut.kabxys_integrated[11][i,j,0]
3073 rhorp=dataOut.kabxys_integrated[8][i,j,0]+dataOut.kabxys_integrated[11][i,j,0]
2664 rhoip=dataOut.kabxys_integrated[10][i,j,0]-dataOut.kabxys_integrated[9][i,j,0]
3074 rhoip=dataOut.kabxys_integrated[10][i,j,0]-dataOut.kabxys_integrated[9][i,j,0]
2665
3075
@@ -2744,6 +3154,7 class DoublePulseACFs(Operation):
2744 #print("p: ",dataOut.p[33,:])
3154 #print("p: ",dataOut.p[33,:])
2745 #exit(1)
3155 #exit(1)
2746 '''
3156 '''
3157
2747 return dataOut
3158 return dataOut
2748
3159
2749 class DoublePulseACFs_PerLag(Operation):
3160 class DoublePulseACFs_PerLag(Operation):
@@ -3040,6 +3451,7 class ElectronDensityFaraday(Operation):
3040 plt.plot(dataOut.bki)
3451 plt.plot(dataOut.bki)
3041 plt.show()
3452 plt.show()
3042 '''
3453 '''
3454
3043 for i in range(2,dataOut.NSHTS-2):
3455 for i in range(2,dataOut.NSHTS-2):
3044 fact=(-0.5/(dataOut.RATE*dataOut.DH))*dataOut.bki[i]
3456 fact=(-0.5/(dataOut.RATE*dataOut.DH))*dataOut.bki[i]
3045 #four-point derivative, no phase unwrapping necessary
3457 #four-point derivative, no phase unwrapping necessary
@@ -3755,6 +4167,7 class DPTemperaturesEstimation(Operation):
3755
4167
3756 plt.show()
4168 plt.show()
3757 '''
4169 '''
4170
3758 return dataOut
4171 return dataOut
3759
4172
3760
4173
@@ -4981,8 +5394,13 class SSheightProfiles(Operation):
4981 dataOut.flagDataAsBlock = True
5394 dataOut.flagDataAsBlock = True
4982 dataOut.ippSeconds = ippSeconds
5395 dataOut.ippSeconds = ippSeconds
4983 dataOut.step = self.step
5396 dataOut.step = self.step
5397
5398
5399 #print("SSH")
4984 #print(numpy.shape(dataOut.data))
5400 #print(numpy.shape(dataOut.data))
4985 #exit(1)
5401 #exit(1)
5402 #print(dataOut.data[0,:,150])
5403 #exit(1)
4986
5404
4987 return dataOut
5405 return dataOut
4988
5406
@@ -6476,125 +6894,451 class IntegrationHP(IntegrationDP):
6476 nint : int
6894 nint : int
6477 Number of integrations.
6895 Number of integrations.
6478
6896
6479 Example
6897 Example
6480 --------
6898 --------
6899
6900 op = proc_unit.addOperation(name='IntegrationHP', optype='other')
6901 op.addParameter(name='nint', value='30', format='int')
6902
6903 """
6904
6905 def __init__(self, **kwargs):
6906
6907 Operation.__init__(self, **kwargs)
6908
6909 self.counter = 0
6910 self.aux = 0
6911
6912 def integration_noise(self,dataOut):
6913
6914 if self.counter == 0:
6915 dataOut.tnoise=numpy.zeros((dataOut.NR),dtype='float32')
6916
6917 dataOut.tnoise+=dataOut.noise_final
6918
6919 def integration_for_long_pulse(self,dataOut):
6920
6921 if self.counter == 0:
6922 dataOut.output_LP_integrated=numpy.zeros((dataOut.NLAG,dataOut.NRANGE,dataOut.NR),order='F',dtype='complex128')
6923
6924 dataOut.output_LP_integrated+=dataOut.output_LP
6925
6926 def run(self,dataOut,nint=None):
6927
6928 dataOut.flagNoData=True
6929
6930 dataOut.nint=nint
6931 dataOut.paramInterval=0#int(dataOut.nint*dataOut.header[7][0]*2 )
6932 dataOut.lat=-11.95
6933 dataOut.lon=-76.87
6934
6935 self.integration_for_long_pulse(dataOut)
6936
6937 self.integration_noise(dataOut)
6938
6939 if self.counter==dataOut.nint-1:
6940 dataOut.nis=dataOut.NSCAN*dataOut.NAVG*dataOut.nint*10
6941 #print(dataOut.tnoise)
6942 #exit(1)
6943 dataOut.tnoise[0]*=0.995
6944 dataOut.tnoise[1]*=0.995
6945 dataOut.pan=dataOut.tnoise[0]/float(dataOut.NSCAN*dataOut.nint*dataOut.NAVG)
6946 dataOut.pbn=dataOut.tnoise[1]/float(dataOut.NSCAN*dataOut.nint*dataOut.NAVG)
6947 #print(self.counter) ToDo: Fix when nint = 1
6948 '''
6949 print("pan: ",dataOut.pan)
6950 print("pbn: ",dataOut.pbn)
6951 #print("tnoise: ",dataOut.tnoise)
6952 exit(1)
6953 '''
6954 #print(dataOut.output_LP_integrated[0,30,2])
6955 #exit(1)
6956 self.integration_for_double_pulse(dataOut)
6957 #if self.counter==dataOut.nint:
6958 #print(dataOut.kabxys_integrated[8][53,6,0]+dataOut.kabxys_integrated[11][53,6,0])
6959 #print(dataOut.kabxys_integrated[8][53,9,0]+dataOut.kabxys_integrated[11][53,9,0])
6960 #exit(1)
6961 return dataOut
6962
6963 class SumFlipsHP(SumFlips):
6964 """Operation to sum the flip and unflip part of certain cross products of the Double Pulse.
6965
6966 Parameters:
6967 -----------
6968 None
6969
6970 Example
6971 --------
6972
6973 op = proc_unit.addOperation(name='SumFlipsHP', optype='other')
6974
6975 """
6976
6977 def __init__(self, **kwargs):
6978
6979 Operation.__init__(self, **kwargs)
6980
6981 def rint2HP(self,dataOut):
6982
6983 dataOut.rnint2=numpy.zeros(dataOut.DPL,'float32')
6984
6985 for l in range(dataOut.DPL):
6986 if(l==0 or (l>=3 and l <=6)):
6987 dataOut.rnint2[l]=0.5/float(dataOut.nint*dataOut.NAVG*16.0)
6988 else:
6989 dataOut.rnint2[l]=0.5/float(dataOut.nint*dataOut.NAVG*8.0)
6990
6991 def run(self,dataOut):
6992
6993 self.rint2HP(dataOut)
6994 self.SumLags(dataOut)
6995
6996 hei = 2
6997 lag = 0
6998 '''
6999 for hei in range(67):
7000 print("hei",hei)
7001 print(dataOut.kabxys_integrated[8][hei,:,0]+dataOut.kabxys_integrated[11][hei,:,0])
7002 print(dataOut.kabxys_integrated[10][hei,:,0]-dataOut.kabxys_integrated[9][hei,:,0])
7003 exit(1)
7004 '''
7005 '''
7006 print("b",(dataOut.kabxys_integrated[4][hei,lag,0]+dataOut.kabxys_integrated[5][hei,lag,0]))
7007 print((dataOut.kabxys_integrated[6][hei,lag,0]+dataOut.kabxys_integrated[7][hei,lag,0]))
7008 print("c",(dataOut.kabxys_integrated[8][hei,lag,0]+dataOut.kabxys_integrated[11][hei,lag,0]))
7009 print((dataOut.kabxys_integrated[10][hei,lag,0]-dataOut.kabxys_integrated[9][hei,lag,0]))
7010 exit(1)
7011 '''
7012 return dataOut
7013
7014
7015 class LongPulseAnalysis(Operation):
7016 """Operation to estimate ACFs, temperatures, total electron density and Hydrogen/Helium fractions from the Long Pulse data.
7017
7018 Parameters:
7019 -----------
7020 NACF : int
7021 .*
7022
7023 Example
7024 --------
7025
7026 op = proc_unit.addOperation(name='LongPulseAnalysis', optype='other')
7027 op.addParameter(name='NACF', value='16', format='int')
7028
7029 """
7030
7031 def __init__(self, **kwargs):
7032
7033 Operation.__init__(self, **kwargs)
7034 self.aux=1
7035
7036 def run(self,dataOut,NACF):
7037
7038 dataOut.NACF=NACF
7039 dataOut.heightList=dataOut.DH*(numpy.arange(dataOut.NACF))
7040 anoise0=dataOut.tnoise[0]
7041 anoise1=anoise0*0.0 #seems to be noise in 1st lag 0.015 before '14
7042 #print(anoise0)
7043 #exit(1)
7044 if self.aux:
7045 #dataOut.cut=31#26#height=31*15=465
7046 self.cal=numpy.zeros((dataOut.NLAG),'float32')
7047 self.drift=numpy.zeros((200),'float32')
7048 self.rdrift=numpy.zeros((200),'float32')
7049 self.ddrift=numpy.zeros((200),'float32')
7050 self.sigma=numpy.zeros((dataOut.NRANGE),order='F',dtype='float32')
7051 self.powera=numpy.zeros((dataOut.NRANGE),order='F',dtype='float32')
7052 self.powerb=numpy.zeros((dataOut.NRANGE),order='F',dtype='float32')
7053 self.perror=numpy.zeros((dataOut.NRANGE),order='F',dtype='float32')
7054 dataOut.ene=numpy.zeros((dataOut.NRANGE),'float32')
7055 self.dpulse=numpy.zeros((dataOut.NACF),'float32')
7056 self.lpulse=numpy.zeros((dataOut.NACF),'float32')
7057 dataOut.lags_LP=numpy.zeros((dataOut.IBITS),order='F',dtype='float32')
7058 self.lagp=numpy.zeros((dataOut.NACF),'float32')
7059 self.u=numpy.zeros((2*dataOut.NACF,2*dataOut.NACF),'float32')
7060 dataOut.ne=numpy.zeros((dataOut.NRANGE),order='F',dtype='float32')
7061 dataOut.te=numpy.zeros((dataOut.NACF),order='F',dtype='float32')
7062 dataOut.ete=numpy.zeros((dataOut.NACF),order='F',dtype='float32')
7063 dataOut.ti=numpy.zeros((dataOut.NACF),order='F',dtype='float32')
7064 dataOut.eti=numpy.zeros((dataOut.NACF),order='F',dtype='float32')
7065 dataOut.ph=numpy.zeros((dataOut.NACF),order='F',dtype='float32')
7066 dataOut.eph=numpy.zeros((dataOut.NACF),order='F',dtype='float32')
7067 dataOut.phe=numpy.zeros((dataOut.NACF),order='F',dtype='float32')
7068 dataOut.ephe=numpy.zeros((dataOut.NACF),order='F',dtype='float32')
7069 dataOut.errors=numpy.zeros((dataOut.IBITS,max(dataOut.NRANGE,dataOut.NSHTS)),order='F',dtype='float32')
7070 dataOut.fit_array_real=numpy.zeros((max(dataOut.NRANGE,dataOut.NSHTS),dataOut.NLAG),order='F',dtype='float32')
7071 dataOut.status=numpy.zeros(1,'float32')
7072 dataOut.tx=240.0 #debería provenir del header #hybrid
7073
7074 for i in range(dataOut.IBITS):
7075 dataOut.lags_LP[i]=float(i)*(dataOut.tx/150.0)/float(dataOut.IBITS) # (float)i*(header.tx/150.0)/(float)IBITS;
7076
7077 self.aux=0
7078
7079 dataOut.cut=30
7080 for i in range(30,15,-1):
7081 if numpy.nanmax(dataOut.acfs_error_to_plot[i,:])>=10 or dataOut.info2[i]==0:
7082 dataOut.cut=i-1
7083
7084 for i in range(dataOut.NLAG):
7085 self.cal[i]=sum(dataOut.output_LP_integrated[i,:,3].real)
7086
7087 #print(numpy.sum(self.cal)) #Coinciden
7088 #exit(1)
7089 self.cal/=float(dataOut.NRANGE)
7090 #print(anoise0)
7091 #print(anoise1)
7092 #exit(1)
7093
7094
7095 #################### PROBAR MÁS INTEGRACIÓN, SINO MODIFICAR VALOR DE "NIS" ####################
7096 # VER dataOut.nProfiles_LP #
7097
7098 '''
7099 #PLOTEAR POTENCIA VS RUIDO, QUIZA SE ESTA REMOVIENDO MUCHA SEÑAL
7100 #print(dataOut.heightList)
7101 import matplotlib.pyplot as plt
7102 plt.plot(10*numpy.log10(dataOut.output_LP_integrated.real[0,:,0]),dataOut.range1)
7103 #plt.plot(10*numpy.log10(dataOut.output_LP_integrated.real[0,:,0]/dataOut.nProfiles_LP),dataOut.range1)
7104 plt.axvline(10*numpy.log10(anoise0),color='k',linestyle='dashed')
7105 plt.grid()
7106 plt.xlim(20,100)
7107 plt.show()
7108 '''
7109
7110
7111 for j in range(dataOut.NACF+2*dataOut.IBITS+2):
7112
7113 dataOut.output_LP_integrated.real[0,j,0]-=anoise0 #lag0 ch0
7114 dataOut.output_LP_integrated.real[1,j,0]-=anoise1 #lag1 ch0
7115
7116 for i in range(1,dataOut.NLAG): #remove cal data from certain lags
7117 dataOut.output_LP_integrated.real[i,j,0]-=self.cal[i]
7118 k=max(j,26) #constant power below range 26
7119 self.powera[j]=dataOut.output_LP_integrated.real[0,k,0]
7120
7121 ## examine drifts here - based on 60 'indep.' estimates
7122 #print(numpy.sum(self.powera))
7123 #exit(1)
7124 #nis=dataOut.NSCAN*dataOut.NAVG*dataOut.nint*10
7125 nis = dataOut.nis
7126 #print("nis",nis)
7127 alpha=beta=delta=0.0
7128 nest=0
7129 gamma=3.0/(2.0*numpy.pi*dataOut.lags_LP[1]*1.0e-3)
7130 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
7131 #print(gamma,beta)
7132 #exit(1)
7133 for i in range(1,3):
7134 gamma=3.0/(2.0*numpy.pi*dataOut.lags_LP[i]*1.0e-3)
7135 #print("gamma",gamma)
7136 for j in range(34,44):
7137 rho2=numpy.abs(dataOut.output_LP_integrated[i,j,0])/numpy.abs(dataOut.output_LP_integrated[0,j,0])
7138 dataOut.dphi2=(1.0/rho2-1.0)/(float(2*nis))
7139 dataOut.dphi2*=gamma**2
7140 pest=gamma*math.atan(dataOut.output_LP_integrated.imag[i,j,0]/dataOut.output_LP_integrated.real[i,j,0])
7141 #print("1",dataOut.output_LP_integrated.imag[i,j,0])
7142 #print("2",dataOut.output_LP_integrated.real[i,j,0])
7143 self.drift[nest]=pest
7144 self.ddrift[nest]=dataOut.dphi2
7145 self.rdrift[nest]=float(nest)
7146 nest+=1
7147
7148 sorted(self.drift[:nest])
7149
7150 #print(dataOut.dphi2)
7151 #exit(1)
7152
7153 for j in range(int(nest/4),int(3*nest/4)):
7154 #i=int(self.rdrift[j])
7155 alpha+=self.drift[j]/self.ddrift[j]
7156 delta+=1.0/self.ddrift[j]
7157
7158 alpha/=delta
7159 delta=1./numpy.sqrt(delta)
7160 vdrift=alpha-beta
7161 dvdrift=delta
7162
7163 #need to develop estimate of complete density profile using all
7164 #available data
6481
7165
6482 op = proc_unit.addOperation(name='IntegrationHP', optype='other')
7166 #estimate sample variances for long-pulse power profile
6483 op.addParameter(name='nint', value='30', format='int')
6484
7167
6485 """
7168 #nis=dataOut.NSCAN*dataOut.NAVG*dataOut.nint
7169 nis = dataOut.nis/10
7170 #print("nis",nis)
6486
7171
6487 def __init__(self, **kwargs):
7172 self.sigma[:dataOut.NACF+2*dataOut.IBITS+2]=((anoise0+self.powera[:dataOut.NACF+2*dataOut.IBITS+2])**2)/float(nis)
7173 #print(self.sigma)
7174 #exit(1)
7175 ioff=1
6488
7176
6489 Operation.__init__(self, **kwargs)
7177 #deconvolve rectangular pulse shape from profile ==> powerb, perror
6490
7178
6491 self.counter = 0
6492 self.aux = 0
6493
7179
6494 def integration_noise(self,dataOut):
7180 ############# START nnlswrap#############
6495
7181
6496 if self.counter == 0:
7182 if dataOut.ut_Faraday>14.0:
6497 dataOut.tnoise=numpy.zeros((dataOut.NR),dtype='float32')
7183 alpha_nnlswrap=20.0
7184 else:
7185 alpha_nnlswrap=30.0
6498
7186
6499 dataOut.tnoise+=dataOut.noise_final
7187 range1_nnls=dataOut.NACF
7188 range2_nnls=dataOut.NACF+dataOut.IBITS-1
6500
7189
6501 def integration_for_long_pulse(self,dataOut):
7190 g_nnlswrap=numpy.zeros((range1_nnls,range2_nnls),'float32')
7191 a_nnlswrap=numpy.zeros((range2_nnls,range2_nnls),'float64')
6502
7192
6503 if self.counter == 0:
7193 for i in range(range1_nnls):
6504 dataOut.output_LP_integrated=numpy.zeros((dataOut.NLAG,dataOut.NRANGE,dataOut.NR),order='F',dtype='complex128')
7194 for j in range(range2_nnls):
7195 if j>=i and j<i+dataOut.IBITS:
7196 g_nnlswrap[i,j]=1.0
7197 else:
7198 g_nnlswrap[i,j]=0.0
6505
7199
6506 dataOut.output_LP_integrated+=dataOut.output_LP
7200 a_nnlswrap[:]=numpy.matmul(numpy.transpose(g_nnlswrap),g_nnlswrap)
6507
7201
6508 def run(self,dataOut,nint=None):
7202 numpy.fill_diagonal(a_nnlswrap,a_nnlswrap.diagonal()+alpha_nnlswrap**2)
6509
7203
6510 dataOut.flagNoData=True
7204 #ERROR ANALYSIS#
6511
7205
6512 dataOut.nint=nint
7206 self.perror[:range2_nnls]=0.0
6513 dataOut.paramInterval=0#int(dataOut.nint*dataOut.header[7][0]*2 )
7207 self.perror[:range2_nnls]=numpy.matmul(1./(self.sigma[dataOut.IBITS+ioff:range1_nnls+dataOut.IBITS+ioff]),g_nnlswrap**2)
6514 dataOut.lat=-11.95
7208 self.perror[:range1_nnls]+=(alpha_nnlswrap**2)/(self.sigma[dataOut.IBITS+ioff:range1_nnls+dataOut.IBITS+ioff])
6515 dataOut.lon=-76.87
7209 self.perror[:range2_nnls]=1.00/self.perror[:range2_nnls]
6516
7210
6517 self.integration_for_long_pulse(dataOut)
7211 b_nnlswrap=numpy.zeros(range2_nnls,'float64')
7212 b_nnlswrap[:]=numpy.matmul(self.powera[dataOut.IBITS+ioff:range1_nnls+dataOut.IBITS+ioff],g_nnlswrap)
6518
7213
6519 self.integration_noise(dataOut)
7214 x_nnlswrap=numpy.zeros(range2_nnls,'float64')
7215 x_nnlswrap[:]=nnls(a_nnlswrap,b_nnlswrap)[0]
6520
7216
6521 if self.counter==dataOut.nint-1:
7217 self.powerb[:range2_nnls]=x_nnlswrap
6522 dataOut.nis=dataOut.NSCAN*dataOut.NAVG*dataOut.nint*10
7218 #print(self.powerb[40])
6523 #print(dataOut.tnoise)
7219 #print(self.powerb[66])
6524 #exit(1)
6525 dataOut.tnoise[0]*=0.995
6526 dataOut.tnoise[1]*=0.995
6527 dataOut.pan=dataOut.tnoise[0]/float(dataOut.NSCAN*dataOut.nint*dataOut.NAVG)
6528 dataOut.pbn=dataOut.tnoise[1]/float(dataOut.NSCAN*dataOut.nint*dataOut.NAVG)
6529 #print(self.counter) ToDo: Fix when nint = 1
6530 '''
6531 print("pan: ",dataOut.pan)
6532 print("pbn: ",dataOut.pbn)
6533 #print("tnoise: ",dataOut.tnoise)
6534 exit(1)
6535 '''
6536 #print(dataOut.output_LP_integrated[0,30,2])
6537 #exit(1)
6538 self.integration_for_double_pulse(dataOut)
6539 #if self.counter==dataOut.nint:
6540 #print(dataOut.kabxys_integrated[8][53,6,0]+dataOut.kabxys_integrated[11][53,6,0])
6541 #print(dataOut.kabxys_integrated[8][53,9,0]+dataOut.kabxys_integrated[11][53,9,0])
6542 #exit(1)
7220 #exit(1)
6543 return dataOut
7221 #############END nnlswrap#############
6544
7222 #print(numpy.sum(numpy.sqrt(self.perror[0:dataOut.NACF])))
6545 class SumFlipsHP(SumFlips):
7223 #print(self.powerb[0:dataOut.NACF])
6546 """Operation to sum the flip and unflip part of certain cross products of the Double Pulse.
7224 #exit(1)
7225 #estimate relative error for deconvolved profile (scaling irrelevant)
7226 #print(dataOut.NACF)
7227 dataOut.ene[0:dataOut.NACF]=numpy.sqrt(self.perror[0:dataOut.NACF])/self.powerb[0:dataOut.NACF]
7228 #print(numpy.sum(dataOut.ene))
7229 #exit(1)
7230 aux=0
6547
7231
6548 Parameters:
7232 for i in range(dataOut.IBITS,dataOut.NACF):
6549 -----------
7233 self.dpulse[i]=self.lpulse[i]=0.0
6550 None
7234 for j in range(dataOut.IBITS):
7235 k=int(i-j)
7236 if k<36-aux and k>16:
7237 self.dpulse[i]+=dataOut.ph2[k]/dataOut.h2[k]
7238 elif k>=36-aux:
7239 self.lpulse[i]+=self.powerb[k]
7240 self.lagp[i]=self.powera[i]
6551
7241
6552 Example
7242 #find scale factor that best merges profiles
6553 --------
6554
7243
6555 op = proc_unit.addOperation(name='SumFlipsHP', optype='other')
7244 qi=sum(self.dpulse[32:dataOut.NACF]**2/(self.lagp[32:dataOut.NACF]+anoise0)**2)
7245 ri=sum((self.dpulse[32:dataOut.NACF]*self.lpulse[32:dataOut.NACF])/(self.lagp[32:dataOut.NACF]+anoise0)**2)
7246 si=sum((self.dpulse[32:dataOut.NACF]*self.lagp[32:dataOut.NACF])/(self.lagp[32:dataOut.NACF]+anoise0)**2)
7247 ui=sum(self.lpulse[32:dataOut.NACF]**2/(self.lagp[32:dataOut.NACF]+anoise0)**2)
7248 vi=sum((self.lpulse[32:dataOut.NACF]*self.lagp[32:dataOut.NACF])/(self.lagp[32:dataOut.NACF]+anoise0)**2)
6556
7249
6557 """
7250 alpha=(si*ui-vi*ri)/(qi*ui-ri*ri)
7251 beta=(qi*vi-ri*si)/(qi*ui-ri*ri)
6558
7252
6559 def __init__(self, **kwargs):
7253 #form density profile estimate, merging rescaled power profiles
7254 #print(dataOut.h2)
7255 #print(numpy.sum(alpha))
7256 #print(numpy.sum(dataOut.ph2))
7257 self.powerb[16:36-aux]=alpha*dataOut.ph2[16:36-aux]/dataOut.h2[16:36-aux]
7258 self.powerb[36-aux:dataOut.NACF]*=beta
6560
7259
6561 Operation.__init__(self, **kwargs)
7260 #form Ne estimate, fill in error estimate at low altitudes
6562
7261
6563 def rint2HP(self,dataOut):
7262 dataOut.ene[0:36-aux]=dataOut.sdp2[0:36-aux]/dataOut.ph2[0:36-aux]
7263 dataOut.ne[:dataOut.NACF]=self.powerb[:dataOut.NACF]*dataOut.h2[:dataOut.NACF]/alpha
7264 #print(numpy.sum(self.powerb))
7265 #print(numpy.sum(dataOut.ene))
7266 #print(numpy.sum(dataOut.ne))
7267 #exit(1)
7268 #now do error propagation: store zero lag error covariance in u
6564
7269
6565 dataOut.rnint2=numpy.zeros(dataOut.DPL,'float32')
7270 nis=dataOut.NSCAN*dataOut.NAVG*dataOut.nint/1 # DLH serious debris removal
6566
7271
6567 for l in range(dataOut.DPL):
7272 for i in range(dataOut.NACF):
6568 if(l==0 or (l>=3 and l <=6)):
7273 for j in range(i,dataOut.NACF):
6569 dataOut.rnint2[l]=0.5/float(dataOut.nint*dataOut.NAVG*16.0)
7274 if j-i>=dataOut.IBITS:
6570 else:
7275 self.u[i,j]=0.0
6571 dataOut.rnint2[l]=0.5/float(dataOut.nint*dataOut.NAVG*8.0)
7276 else:
7277 self.u[i,j]=dataOut.output_LP_integrated.real[j-i,i,0]**2/float(nis)
7278 self.u[i,j]*=(anoise0+dataOut.output_LP_integrated.real[0,i,0])/dataOut.output_LP_integrated.real[0,i,0]
7279 self.u[i,j]*=(anoise0+dataOut.output_LP_integrated.real[0,j,0])/dataOut.output_LP_integrated.real[0,j,0]
6572
7280
6573 def run(self,dataOut):
7281 self.u[j,i]=self.u[i,j]
6574
7282
6575 self.rint2HP(dataOut)
7283 #now error analyis for lag product matrix (diag), place in acf_err
6576 self.SumLags(dataOut)
6577
7284
6578 hei = 2
7285 for i in range(dataOut.NACF):
6579 lag = 0
7286 for j in range(dataOut.IBITS):
7287 if j==0:
7288 dataOut.errors[0,i]=numpy.sqrt(self.u[i,i])
7289 else:
7290 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))
6580 '''
7291 '''
6581 for hei in range(67):
7292 print(numpy.sum(dataOut.output_LP_integrated))
6582 print("hei",hei)
7293 print(numpy.sum(dataOut.errors))
6583 print(dataOut.kabxys_integrated[8][hei,:,0]+dataOut.kabxys_integrated[11][hei,:,0])
7294 print(numpy.sum(self.powerb))
6584 print(dataOut.kabxys_integrated[10][hei,:,0]-dataOut.kabxys_integrated[9][hei,:,0])
7295 print(numpy.sum(dataOut.ne))
7296 print(numpy.sum(dataOut.lags_LP))
7297 print(numpy.sum(dataOut.thb))
7298 print(numpy.sum(dataOut.bfm))
7299 print(numpy.sum(dataOut.te))
7300 print(numpy.sum(dataOut.ete))
7301 print(numpy.sum(dataOut.ti))
7302 print(numpy.sum(dataOut.eti))
7303 print(numpy.sum(dataOut.ph))
7304 print(numpy.sum(dataOut.eph))
7305 print(numpy.sum(dataOut.phe))
7306 print(numpy.sum(dataOut.ephe))
7307 print(numpy.sum(dataOut.range1))
7308 print(numpy.sum(dataOut.ut))
7309 print(numpy.sum(dataOut.NACF))
7310 print(numpy.sum(dataOut.fit_array_real))
7311 print(numpy.sum(dataOut.status))
7312 print(numpy.sum(dataOut.NRANGE))
7313 print(numpy.sum(dataOut.IBITS))
6585 exit(1)
7314 exit(1)
6586 '''
7315 '''
6587 '''
7316 '''
6588 print("b",(dataOut.kabxys_integrated[4][hei,lag,0]+dataOut.kabxys_integrated[5][hei,lag,0]))
7317 print(dataOut.te2[13:16])
6589 print((dataOut.kabxys_integrated[6][hei,lag,0]+dataOut.kabxys_integrated[7][hei,lag,0]))
7318 print(numpy.sum(dataOut.te2))
6590 print("c",(dataOut.kabxys_integrated[8][hei,lag,0]+dataOut.kabxys_integrated[11][hei,lag,0]))
6591 print((dataOut.kabxys_integrated[10][hei,lag,0]-dataOut.kabxys_integrated[9][hei,lag,0]))
6592 exit(1)
7319 exit(1)
6593 '''
7320 '''
6594 return dataOut
7321 print("Success")
7322 #print(dataOut.NRANGE)
7323 with suppress_stdout_stderr():
7324 #pass
7325 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)
6595
7326
7327 print("status: ",dataOut.status)
6596
7328
6597 class LongPulseAnalysis(Operation):
7329 if dataOut.status>=3.5:
7330 dataOut.te[:]=numpy.nan
7331 dataOut.ete[:]=numpy.nan
7332 dataOut.ti[:]=numpy.nan
7333 dataOut.eti[:]=numpy.nan
7334 dataOut.ph[:]=numpy.nan
7335 dataOut.eph[:]=numpy.nan
7336 dataOut.phe[:]=numpy.nan
7337 dataOut.ephe[:]=numpy.nan
7338
7339 return dataOut
7340
7341 class LongPulseAnalysis_V2(Operation):
6598 """Operation to estimate ACFs, temperatures, total electron density and Hydrogen/Helium fractions from the Long Pulse data.
7342 """Operation to estimate ACFs, temperatures, total electron density and Hydrogen/Helium fractions from the Long Pulse data.
6599
7343
6600 Parameters:
7344 Parameters:
@@ -6853,36 +7597,7 class LongPulseAnalysis(Operation):
6853 dataOut.errors[0,i]=numpy.sqrt(self.u[i,i])
7597 dataOut.errors[0,i]=numpy.sqrt(self.u[i,i])
6854 else:
7598 else:
6855 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))
7599 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))
6856 '''
7600
6857 print(numpy.sum(dataOut.output_LP_integrated))
6858 print(numpy.sum(dataOut.errors))
6859 print(numpy.sum(self.powerb))
6860 print(numpy.sum(dataOut.ne))
6861 print(numpy.sum(dataOut.lags_LP))
6862 print(numpy.sum(dataOut.thb))
6863 print(numpy.sum(dataOut.bfm))
6864 print(numpy.sum(dataOut.te))
6865 print(numpy.sum(dataOut.ete))
6866 print(numpy.sum(dataOut.ti))
6867 print(numpy.sum(dataOut.eti))
6868 print(numpy.sum(dataOut.ph))
6869 print(numpy.sum(dataOut.eph))
6870 print(numpy.sum(dataOut.phe))
6871 print(numpy.sum(dataOut.ephe))
6872 print(numpy.sum(dataOut.range1))
6873 print(numpy.sum(dataOut.ut))
6874 print(numpy.sum(dataOut.NACF))
6875 print(numpy.sum(dataOut.fit_array_real))
6876 print(numpy.sum(dataOut.status))
6877 print(numpy.sum(dataOut.NRANGE))
6878 print(numpy.sum(dataOut.IBITS))
6879 exit(1)
6880 '''
6881 '''
6882 print(dataOut.te2[13:16])
6883 print(numpy.sum(dataOut.te2))
6884 exit(1)
6885 '''
6886 print("Success")
7601 print("Success")
6887 with suppress_stdout_stderr():
7602 with suppress_stdout_stderr():
6888 #pass
7603 #pass
@@ -6900,7 +7615,6 class LongPulseAnalysis(Operation):
6900
7615
6901 return dataOut
7616 return dataOut
6902
7617
6903
6904 class PulsePairVoltage(Operation):
7618 class PulsePairVoltage(Operation):
6905 '''
7619 '''
6906 Function PulsePair(Signal Power, Velocity)
7620 Function PulsePair(Signal Power, Velocity)
Show file before | Show file after | Hide comments
@@ -74,6 +74,7 setup(
74 cmdclass = {'build_ext': build_ext},
74 cmdclass = {'build_ext': build_ext},
75 ext_modules=[
75 ext_modules=[
76 Extension("schainpy.model.data._noise", ["schainc/_noise.c"]),
76 Extension("schainpy.model.data._noise", ["schainc/_noise.c"]),
77 Extension("schainpy.model.data._HS_algorithm", ["schainc/_HS_algorithm.c"]),
77 ],
78 ],
78 setup_requires = ["numpy"],
79 setup_requires = ["numpy"],
79 install_requires = [
80 install_requires = [
General Comments 0
You need to be logged in to leave comments. Login now