##// 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 653 elif not ok:
654 654 break
655 655 #print("****************************************************end")
656 #exit(1)
656 657 if n == 0:
657 658 err = True
658 659
Show file before | Show file after | Hide comments
@@ -42,6 +42,7 class SpectraPlot(Plot):
42 42 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
43 43 data['spc'] = spc
44 44 data['rti'] = dataOut.getPower()
45 #print("NormFactor: ",dataOut.normFactor)
45 46 #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
46 47 if hasattr(dataOut, 'LagPlot'): #Double Pulse
47 48 max_hei_id = dataOut.nHeights - 2*dataOut.LagPlot
@@ -101,10 +102,11 class SpectraPlot(Plot):
101 102 self.xmin = self.xmin if self.xmin else numpy.nanmin(x)#-self.xmax
102 103 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
103 104 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
105
104 106 ax.plt = ax.pcolormesh(x, y, z[n].T,
105 107 vmin=self.zmin,
106 108 vmax=self.zmax,
107 cmap=plt.get_cmap(self.colormap)
109 cmap=plt.get_cmap(self.colormap),
108 110 )
109 111
110 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 266 class ETempRTIPlot(RTIPlot):
269 267
270 268 '''
@@ -351,7 +349,6 class ETempRTIPlot(RTIPlot):
351 349 )
352 350
353 351
354
355 352 class ITempRTIPlot(ETempRTIPlot):
356 353
357 354 '''
@@ -372,7 +369,6 class ITempRTIPlot(ETempRTIPlot):
372 369 return data, meta
373 370
374 371
375
376 372 class HFracRTIPlot(ETempRTIPlot):
377 373
378 374 '''
@@ -545,6 +541,8 class TempsHPPlot(Plot):
545 541 ax.errorbar(Ti, self.y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
546 542 plt.legend(loc='lower right')
547 543 ax.yaxis.set_minor_locator(MultipleLocator(15))
544 ax.grid(which='minor')
545
548 546
549 547 class FracsHPPlot(Plot):
550 548 '''
Show file before | Show file after | Hide comments
@@ -70,14 +70,16 class ProcessingUnit(object):
70 70 def call(self, **kwargs):
71 71 '''
72 72 '''
73
73 #print("call")
74 74 try:
75 #print("try")
76 75 if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
77 #print("Try")
78 #print(self.dataIn)
79 #print(self.dataIn.flagNoData)
80 return self.dataIn.isReady()
76 #if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error and not self.dataIn.runNextUnit:
77 if self.dataIn.runNextUnit:
78 #print("SUCCESSSSSSS")
79 #exit(1)
80 return not self.dataIn.isReady()
81 else:
82 return self.dataIn.isReady()
81 83 elif self.dataIn is None or not self.dataIn.error:
82 84 #print([getattr(self, at) for at in self.inputs])
83 85 #print("Elif 1")
@@ -112,6 +114,7 class ProcessingUnit(object):
112 114 #elif optype == 'external' and self.dataOut.isReady():
113 115 #op.queue.put(copy.deepcopy(self.dataOut))
114 116 #print(not self.dataOut.isReady())
117
115 118 try:
116 119 if self.dataOut.runNextUnit:
117 120 runNextUnit = self.dataOut.runNextUnit
@@ -125,6 +128,7 class ProcessingUnit(object):
125 128 #if not self.dataOut.isReady():
126 129 #return 'Error' if self.dataOut.error else input()
127 130 #print("NexT",runNextUnit)
131 #print("error: ",self.dataOut.error)
128 132 return 'Error' if self.dataOut.error else runNextUnit# self.dataOut.isReady()
129 133
130 134 def setup(self):
Show file before | Show file after | Hide comments
@@ -1164,18 +1164,34 class Oblique_Gauss_Fit(Operation):
1164 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 1165 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1166 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 1179 A1f = popt.x[0]; B1f = popt.x[1]; C1f = popt.x[2]; K1f = popt.x[3]
1169 1180 A2f = popt.x[4]; B2f = popt.x[5]; C2f = popt.x[6]; K2f = popt.x[7]
1170 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 1187 aux1 = self.gaussian_skew(freq, A1f, B1f, C1f, K1f, Df)
1173 1188 doppler1 = freq[numpy.argmax(aux1)]
1174 1189
1175 1190 aux2 = self.gaussian_skew(freq, A2f, B2f, C2f, K2f, Df)
1176 1191 doppler2 = freq[numpy.argmax(aux2)]
1177
1178 return A1f, B1f, C1f, K1f, A2f, B2f, C2f, K2f, Df, doppler1, doppler2
1192 #print("error",error)
1193 #exit(1)
1194 return A1f, B1f, C1f, K1f, A2f, B2f, C2f, K2f, Df, doppler1, doppler2, error
1179 1195
1180 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 1292 dataOut.Oblique_params = numpy.ones((1,10,dataOut.nHeights))*numpy.NAN
1277 1293 elif mode == 9:
1278 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 1297 dataOut.VelRange = x
1281 1298
@@ -1303,6 +1320,7 class Oblique_Gauss_Fit(Operation):
1303 1320 else:
1304 1321
1305 1322 for hei in itertools.chain(l1, l2):
1323 #for hei in range(79,81):
1306 1324
1307 1325 try:
1308 1326
@@ -1322,11 +1340,14 class Oblique_Gauss_Fit(Operation):
1322 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 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 1344 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,10,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei]))
1327 1345 #print(hei)
1328 1346 #print(dataOut.Oblique_params[0,10,hei])
1329 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 1352 else:
1332 1353 spc_fit, A1, B1, C1, D1 = self.Gauss_fit_2(spc,x,'first')
@@ -5575,7 +5596,7 class MergeProc(ProcessingUnit):
5575 5596 def __init__(self):
5576 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 5601 self.dataOut = getattr(self, self.inputs[0])
5581 5602 data_inputs = [getattr(self, attr) for attr in self.inputs]
@@ -5636,3 +5657,62 class MergeProc(ProcessingUnit):
5636 5657 self.dataOut.freqRange = self.dataOut.getFreqRange(1)/1000.
5637 5658
5638 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 4 from schainpy.model.data.jrodata import Spectra
5 5 from schainpy.model.data.jrodata import hildebrand_sekhon
6 6
7 class SpectraAFCProc(ProcessingUnit):
7 class SpectraAFCProc_V0(ProcessingUnit):
8 8
9 9 def __init__(self, **kwargs):
10 10
@@ -140,12 +140,16 class SpectraAFCProc(ProcessingUnit):
140 140
141 141 if self.dataIn.type == "Spectra":
142 142 self.dataOut.copy(self.dataIn)
143 #print(self.dataOut.data.shape)
144 #exit(1)
143 145 spc = self.dataOut.data_spc
144 146 data = numpy.fft.fftshift( spc, axes=(1,))
145 147 data = numpy.fft.ifft(data, axis=1)
146 #data = numpy.fft.fftshift( data, axes=(1,))
147 #acf = numpy.abs(data) # Autocorrelacion LLAMAR A ESTE VALOR ACF
148 acf = data
148 data = numpy.fft.fftshift( data, axes=(1,))
149 acf = numpy.abs(data) # Autocorrelacion LLAMAR A ESTE VALOR ACF
150 acf = data #Comentarlo?
151 print("acf",acf[0,:,150])
152 exit(1)
149 153 #'''
150 154 if real:
151 155 acf = data.real
@@ -167,7 +171,7 class SpectraAFCProc(ProcessingUnit):
167 171 acf[i,:,j]= acf[i,:,j] / numpy.max(numpy.abs(acf[i,:,j]))
168 172 '''
169 173 self.dataOut.data_acf = acf
170 self.dataOut.data_spc = acf
174 self.dataOut.data_spc = acf.real
171 175 #print(self.dataOut.data_acf[0,:,0])
172 176 #exit(1)
173 177 '''
@@ -183,6 +187,7 class SpectraAFCProc(ProcessingUnit):
183 187 #print(self.dataOut.data_spc.shape)
184 188 #print(self.dataOut.data_acf.shape)
185 189 '''
190 '''
186 191 import matplotlib.pyplot as plt
187 192 hei = 10
188 193 #print(self.dataOut.heightList)
@@ -197,7 +202,851 class SpectraAFCProc(ProcessingUnit):
197 202 plt.ylim(0,1000)
198 203 plt.show()
199 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 1050 return True
202 1051
203 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 17 from schainpy.model.data.jrodata import Spectra
18 18 from schainpy.model.data.jrodata import hildebrand_sekhon
19 19 from schainpy.utils import log
20 from schainpy.model.data import _HS_algorithm
20 21
21 22 from time import time, mktime, strptime, gmtime, ctime
22 23
@@ -87,6 +88,7 class SpectraLagProc(ProcessingUnit):
87 88 """
88 89 #print(self.buffer[1,:,0])
89 90 #exit(1)
91 #print("buffer shape",self.buffer.shape)
90 92 fft_volt = numpy.fft.fft(
91 93 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
92 94 fft_volt = fft_volt.astype(numpy.dtype('complex'))
@@ -208,6 +210,17 class SpectraLagProc(ProcessingUnit):
208 210 self.dataOut.ByLags=ByLags
209 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 224 if self.dataIn.type == "Spectra":
212 225 self.dataOut.copy(self.dataIn)
213 226 if shift_fft:
@@ -224,6 +237,7 class SpectraLagProc(ProcessingUnit):
224 237 elif self.dataIn.type == "Voltage":
225 238
226 239 if not self.dataOut.ByLags:
240 #self.dataOut.data = self.dataIn.data
227 241 self.VoltageType(nFFTPoints,nProfiles,ippFactor,pairsList)
228 242 else:
229 243 self.dataOut.nLags = nLags
@@ -1232,7 +1246,9 class IntegrationFaradaySpectra(Operation):
1232 1246 if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
1233 1247 continue
1234 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 1253 if i==1 and l==0 and k==37:
1238 1254 print("j",j)
@@ -1632,7 +1648,9 class IntegrationFaradaySpectra2(Operation):
1632 1648 print(buffer)
1633 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 1655 indexes.append(index)
1638 1656 #sortIDs.append(sortID)
@@ -1771,7 +1789,7 class IntegrationFaradaySpectra2(Operation):
1771 1789 exit(1)
1772 1790 '''
1773 1791
1774 print(self.nLags)
1792 #print(self.nLags)
1775 1793 '''
1776 1794 if self.nLags == 16:
1777 1795 self.nLags = 3
@@ -1821,7 +1839,9 class IntegrationFaradaySpectra2(Operation):
1821 1839 print(buffer)
1822 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 1846 indexes.append(index)
1827 1847 #sortIDs.append(sortID)
@@ -2251,7 +2271,9 class IntegrationFaradaySpectra3(Operation): #This class should manage data with
2251 2271 print(buffer)
2252 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 2278 indexes.append(index)
2257 2279 #sortIDs.append(sortID)
@@ -2440,7 +2462,9 class IntegrationFaradaySpectra3(Operation): #This class should manage data with
2440 2462 print(buffer)
2441 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 2469 indexes.append(index)
2446 2470 #sortIDs.append(sortID)
@@ -2572,7 +2596,9 class IntegrationFaradaySpectra3(Operation): #This class should manage data with
2572 2596 #buffer=buffer1[:,j]
2573 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 2603 indexes.append(index)
2578 2604 #sortIDs.append(sortID)
@@ -2876,7 +2902,7 class IntegrationFaradaySpectraNoLags(Operation):
2876 2902 self.__profIndex += 1
2877 2903
2878 2904 return
2879
2905 #'''
2880 2906 def hildebrand_sekhon_Integration(self,data,navg):
2881 2907
2882 2908 sortdata = numpy.sort(data, axis=None)
@@ -2894,6 +2920,8 class IntegrationFaradaySpectraNoLags(Operation):
2894 2920 sumq += sortdata[j]**2
2895 2921 if j > nums_min:
2896 2922 rtest = float(j)/(j-1) + 1.0/navg
2923 #print(rtest)
2924 #print(sump)
2897 2925 if ((sumq*j) > (rtest*sump**2)):
2898 2926 j = j - 1
2899 2927 sump = sump - sortdata[j]
@@ -2903,7 +2931,7 class IntegrationFaradaySpectraNoLags(Operation):
2903 2931 #lnoise = sump / j
2904 2932
2905 2933 return j,sortID
2906
2934 #'''
2907 2935 def pushData(self):
2908 2936 """
2909 2937 Return the sum of the last profiles and the profiles used in the sum.
@@ -2938,12 +2966,20 class IntegrationFaradaySpectraNoLags(Operation):
2938 2966 if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
2939 2967 continue
2940 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 2978 indexes.append(index)
2944 2979 #sortIDs.append(sortID)
2945 2980 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
2946
2981 #if k == 100:
2982 # exit(1)
2947 2983 outliers_IDs=numpy.array(outliers_IDs)
2948 2984 outliers_IDs=outliers_IDs.ravel()
2949 2985 outliers_IDs=numpy.unique(outliers_IDs)
@@ -3332,7 +3368,9 class IncohIntLag(Operation):
3332 3368 return dataOut
3333 3369
3334 3370 dataOut.flagNoData = True
3335
3371 #print("incohint")
3372 #print("IncohInt",dataOut.data_spc.shape)
3373 #print("IncohInt",dataOut.data_cspc.shape)
3336 3374 if not self.isConfig:
3337 3375 self.setup(n, timeInterval, overlapping)
3338 3376 self.isConfig = True
@@ -3352,7 +3390,7 class IncohIntLag(Operation):
3352 3390 dataOut.dataLag_spc,
3353 3391 dataOut.dataLag_cspc,
3354 3392 dataOut.dataLag_dc)
3355
3393 #print("Incoh Int: ",self.__profIndex,n)
3356 3394 if self.__dataReady:
3357 3395
3358 3396 if not dataOut.ByLags:
@@ -3369,7 +3407,8 class IncohIntLag(Operation):
3369 3407 #print(numpy.sum(dataOut.dataLag_spc[1,:,100,2]))
3370 3408 #exit(1)
3371 3409
3372 #print("HERE")
3410 #print("INCOH INT DONE")
3411 #exit(1)
3373 3412 '''
3374 3413 print(numpy.sum(dataOut.dataLag_spc[0,:,20,10])/32)
3375 3414 print(numpy.sum(dataOut.dataLag_spc[1,:,20,10])/32)
@@ -3580,6 +3619,8 class IncohInt(Operation):
3580 3619 dataOut.nIncohInt *= self.n
3581 3620 dataOut.utctime = avgdatatime
3582 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 3624 #exit(1)
3584 3625 #print(numpy.sum(dataOut.data_spc[0,:,53]*numpy.conjugate(dataOut.data_spc[0,:,53])))
3585 3626 #print(numpy.sum(dataOut.data_spc[1,:,53]*numpy.conjugate(dataOut.data_spc[1,:,53])))
@@ -3773,7 +3814,6 class SpectraDataToFaraday(Operation):
3773 3814 dataOut.lat=-11.95
3774 3815 dataOut.lon=-76.87
3775 3816
3776
3777 3817 self.normFactor(dataOut)
3778 3818
3779 3819 dataOut.NDP=dataOut.nHeights
@@ -3848,7 +3888,7 class SpectraDataToFaraday(Operation):
3848 3888 #print("Noise dB: ",10*numpy.log10(dataOut.tnoise))
3849 3889 #exit(1)
3850 3890 #dataOut.pan=dataOut.tnoise[0]/float(dataOut.nProfiles_LP*dataOut.nIncohInt)
3851
3891 print("done")
3852 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 3948 #print(dataOut.dataLag_spc[:,:,:,1]/dataOut.data_spc)
3909 3949 #exit(1)
@@ -3922,6 +3962,20 class SpectraDataToHybrid(SpectraDataToFaraday):
3922 3962
3923 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 3979 def normFactor(self,dataOut):
3926 3980 dataOut.rnint2=numpy.zeros(dataOut.DPL,'float32')
3927 3981 for l in range(dataOut.DPL):
@@ -4011,3 +4065,115 class SpectraDataToHybrid(SpectraDataToFaraday):
4011 4065 #exit(1)
4012 4066
4013 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 442 if not channelList:
443 443 data[:,:] = data[:,:]*self.flip
444 444 else:
445 channelList=[1]
445 #channelList=[1]
446 446
447 447 for thisChannel in channelList:
448 448 if thisChannel not in dataOut.channelList:
449 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 462 self.flip *= -1.
454 463
@@ -1093,28 +1102,56 class LagsReshapeDP_V2(Operation):
1093 1102 #print(dataOut.data[1,:12,:15])
1094 1103 #exit(1)
1095 1104 #print(numpy.shape(dataOut.data))
1096 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")
1105 #print(dataOut.profileIndex)
1108 1106
1109 dataOut.data = dataOut.data[:,:,200:]
1110 self.data_buffer = []
1111 #dataOut.flagNoData = False
1107 if not dataOut.flagDataAsBlock:
1112 1108
1113 else:
1114 self.data_buffer.append(dataOut.data)
1115 #print(numpy.shape(dataOut.data))
1116 #exit(1)
1109 dataOut.flagNoData = True
1110 #print("nProfiles: ",dataOut.nProfiles)
1111 #if dataOut.profileIndex == 140:
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 1156 return dataOut
1120 1157
@@ -2523,8 +2560,8 class CleanCohEchoes(Operation):
2523 2560
2524 2561 return dataOut
2525 2562
2526 class NoisePower(Operation):
2527 """Operation to get noise power from the integrated data of the Double Pulse.
2563 class CleanCohEchoesHP(Operation):
2564 """Operation to clean coherent echoes.
2528 2565
2529 2566 Parameters:
2530 2567 -----------
@@ -2533,7 +2570,7 class NoisePower(Operation):
2533 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 2579 Operation.__init__(self, **kwargs)
2543 2580
2544 def hildebrand(self,dataOut,data):
2545
2546 divider=10 # divider was originally 10
2547 noise=0.0
2548 n1=0
2549 n2=int(dataOut.NDP/2)
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
2581 def remove_coh(self,pow):
2582 #print("pow inside: ",pow)
2583 #print(pow.shape)
2584 q75,q25 = numpy.percentile(pow,[75,25],axis=0)
2585 #print(q75,q25)
2586 intr_qr = q75-q25
2575 2587
2576 noise= sump/j
2577 stdv=numpy.sqrt((sumq- noise*noise)/(j-1))
2578 return noise
2588 max = q75+(1.5*intr_qr)
2589 min = q25-(1.5*intr_qr)
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')
2583 av=numpy.zeros(dataOut.NDP,'float32')
2584 dataOut.pnoise=numpy.zeros(dataOut.NR,'float32')
2593 #print("Max: ",max)
2594 #print("Min: ",min)
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
2587 p[1,:,:]=dataOut.kabxys_integrated[6][:,:,0]+dataOut.kabxys_integrated[7][:,:,0] #total power for channel 1
2596 return pow
2588 2597
2589 for i in range(dataOut.NR):
2590 dataOut.pnoise[i]=0.0
2591 for k in range(dataOut.DPL):
2592 dataOut.pnoise[i]+= self.hildebrand(dataOut,p[i,:,k])
2598 def mad_based_outlier_V0(self, points, thresh=3.5):
2599 #print("points: ",points)
2600 if len(points.shape) == 1:
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
2598 dataOut.pbn=1.0*dataOut.pnoise[1] # weights could change
2599 print("pan: ",dataOut.pan)
2600 print("pbn: ",dataOut.pbn)
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
2613 median = numpy.nanmedian(points)
2614 diff = (points - median)**2
2615 diff = numpy.sqrt(diff)
2616 med_abs_deviation = numpy.nanmedian(diff)
2606 2617
2618 modified_z_score = 0.6745 * diff / med_abs_deviation
2607 2619
2608 class DoublePulseACFs(Operation):
2609 """Operation to get the ACFs of the Double Pulse.
2620 return modified_z_score > thresh
2610 2621
2611 Parameters:
2612 -----------
2613 None
2622 def removeSpreadF_V0(self,dataOut):
2623 for i in range(11):
2624 print("BEFORE Chb: ",i,dataOut.kabxys_integrated[6][:,i,0])
2625 #exit(1)
2614 2626
2615 Example
2616 --------
2627 #Removing echoes greater than 35 dB
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)
2625 self.aux=1
2696 (uniq, freq) = (numpy.unique(outliers_IDs, return_counts=True))
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:
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)
2718 dataOut.flagSpreadF = True
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):
2649 for j in range(dataOut.DPL):
2650 ################# Total power
2651 pa=numpy.abs(dataOut.kabxys_integrated[4][i,j,0]+dataOut.kabxys_integrated[5][i,j,0])
2652 pb=numpy.abs(dataOut.kabxys_integrated[6][i,j,0]+dataOut.kabxys_integrated[7][i,j,0])
2653 st4=pa*pb
2654 '''
2655 if i > id[0]:
2656 dataOut.p[i,j] pa-dataOut.pan
2657 else:
2658 dataOut.p[i,j]=pa+pb-(dataOut.pan+dataOut.pbn)
2659 '''
2660 dataOut.p[i,j]=pa+pb-(dataOut.pan+dataOut.pbn)
2661 dataOut.sdp[i,j]=2*dataOut.rnint2[j]*((pa+pb)*(pa+pb))
2725 '''
2726 for lag in range(11):
2727 #print("Lag: ",lag)
2728 outliers = self.mad_based_outlier(dataOut.kabxys_integrated[6][minIndex:maxIndex,lag,0], thresh=2.)
2729 dataOut.kabxys_integrated[6][minIndex:maxIndex,lag,0][outliers == True] = numpy.nan
2730 '''
2731 #dataOut.kabxys_integrated[4][minIndex:,:,0] = self.remove_coh(dataOut.kabxys_integrated[4][minIndex:,:,0])
2732 '''
2733 import matplotlib.pyplot as plt
2734 for i in range(11):
2735 plt.plot(dataOut.kabxys_integrated[6][:,i,0],dataOut.heightList,label='{}'.format(i))
2736 plt.xlim(0,2000)
2737 plt.legend()
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 3071 ## ACF
3072
2663 3073 rhorp=dataOut.kabxys_integrated[8][i,j,0]+dataOut.kabxys_integrated[11][i,j,0]
2664 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 3154 #print("p: ",dataOut.p[33,:])
2745 3155 #exit(1)
2746 3156 '''
3157
2747 3158 return dataOut
2748 3159
2749 3160 class DoublePulseACFs_PerLag(Operation):
@@ -3040,6 +3451,7 class ElectronDensityFaraday(Operation):
3040 3451 plt.plot(dataOut.bki)
3041 3452 plt.show()
3042 3453 '''
3454
3043 3455 for i in range(2,dataOut.NSHTS-2):
3044 3456 fact=(-0.5/(dataOut.RATE*dataOut.DH))*dataOut.bki[i]
3045 3457 #four-point derivative, no phase unwrapping necessary
@@ -3755,6 +4167,7 class DPTemperaturesEstimation(Operation):
3755 4167
3756 4168 plt.show()
3757 4169 '''
4170
3758 4171 return dataOut
3759 4172
3760 4173
@@ -4981,8 +5394,13 class SSheightProfiles(Operation):
4981 5394 dataOut.flagDataAsBlock = True
4982 5395 dataOut.ippSeconds = ippSeconds
4983 5396 dataOut.step = self.step
5397
5398
5399 #print("SSH")
4984 5400 #print(numpy.shape(dataOut.data))
4985 5401 #exit(1)
5402 #print(dataOut.data[0,:,150])
5403 #exit(1)
4986 5404
4987 5405 return dataOut
4988 5406
@@ -6476,125 +6894,451 class IntegrationHP(IntegrationDP):
6476 6894 nint : int
6477 6895 Number of integrations.
6478 6896
6479 Example
6480 --------
6897 Example
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')
6483 op.addParameter(name='nint', value='30', format='int')
7166 #estimate sample variances for long-pulse power profile
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:
6497 dataOut.tnoise=numpy.zeros((dataOut.NR),dtype='float32')
7182 if dataOut.ut_Faraday>14.0:
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:
6504 dataOut.output_LP_integrated=numpy.zeros((dataOut.NLAG,dataOut.NRANGE,dataOut.NR),order='F',dtype='complex128')
7193 for i in range(range1_nnls):
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
6513 dataOut.paramInterval=0#int(dataOut.nint*dataOut.header[7][0]*2 )
6514 dataOut.lat=-11.95
6515 dataOut.lon=-76.87
7206 self.perror[:range2_nnls]=0.0
7207 self.perror[:range2_nnls]=numpy.matmul(1./(self.sigma[dataOut.IBITS+ioff:range1_nnls+dataOut.IBITS+ioff]),g_nnlswrap**2)
7208 self.perror[:range1_nnls]+=(alpha_nnlswrap**2)/(self.sigma[dataOut.IBITS+ioff:range1_nnls+dataOut.IBITS+ioff])
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:
6522 dataOut.nis=dataOut.NSCAN*dataOut.NAVG*dataOut.nint*10
6523 #print(dataOut.tnoise)
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])
7217 self.powerb[:range2_nnls]=x_nnlswrap
7218 #print(self.powerb[40])
7219 #print(self.powerb[66])
6542 7220 #exit(1)
6543 return dataOut
6544
6545 class SumFlipsHP(SumFlips):
6546 """Operation to sum the flip and unflip part of certain cross products of the Double Pulse.
7221 #############END nnlswrap#############
7222 #print(numpy.sum(numpy.sqrt(self.perror[0:dataOut.NACF])))
7223 #print(self.powerb[0:dataOut.NACF])
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:
6549 -----------
6550 None
7232 for i in range(dataOut.IBITS,dataOut.NACF):
7233 self.dpulse[i]=self.lpulse[i]=0.0
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
6553 --------
7242 #find scale factor that best merges profiles
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):
6568 if(l==0 or (l>=3 and l <=6)):
6569 dataOut.rnint2[l]=0.5/float(dataOut.nint*dataOut.NAVG*16.0)
6570 else:
6571 dataOut.rnint2[l]=0.5/float(dataOut.nint*dataOut.NAVG*8.0)
7272 for i in range(dataOut.NACF):
7273 for j in range(i,dataOut.NACF):
7274 if j-i>=dataOut.IBITS:
7275 self.u[i,j]=0.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)
6576 self.SumLags(dataOut)
7283 #now error analyis for lag product matrix (diag), place in acf_err
6577 7284
6578 hei = 2
6579 lag = 0
7285 for i in range(dataOut.NACF):
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):
6582 print("hei",hei)
6583 print(dataOut.kabxys_integrated[8][hei,:,0]+dataOut.kabxys_integrated[11][hei,:,0])
6584 print(dataOut.kabxys_integrated[10][hei,:,0]-dataOut.kabxys_integrated[9][hei,:,0])
7292 print(numpy.sum(dataOut.output_LP_integrated))
7293 print(numpy.sum(dataOut.errors))
7294 print(numpy.sum(self.powerb))
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 7314 exit(1)
6586 7315 '''
6587 7316 '''
6588 print("b",(dataOut.kabxys_integrated[4][hei,lag,0]+dataOut.kabxys_integrated[5][hei,lag,0]))
6589 print((dataOut.kabxys_integrated[6][hei,lag,0]+dataOut.kabxys_integrated[7][hei,lag,0]))
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]))
7317 print(dataOut.te2[13:16])
7318 print(numpy.sum(dataOut.te2))
6592 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 7342 """Operation to estimate ACFs, temperatures, total electron density and Hydrogen/Helium fractions from the Long Pulse data.
6599 7343
6600 7344 Parameters:
@@ -6853,36 +7597,7 class LongPulseAnalysis(Operation):
6853 7597 dataOut.errors[0,i]=numpy.sqrt(self.u[i,i])
6854 7598 else:
6855 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 '''
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 '''
7600
6886 7601 print("Success")
6887 7602 with suppress_stdout_stderr():
6888 7603 #pass
@@ -6900,7 +7615,6 class LongPulseAnalysis(Operation):
6900 7615
6901 7616 return dataOut
6902 7617
6903
6904 7618 class PulsePairVoltage(Operation):
6905 7619 '''
6906 7620 Function PulsePair(Signal Power, Velocity)
Show file before | Show file after | Hide comments
@@ -74,6 +74,7 setup(
74 74 cmdclass = {'build_ext': build_ext},
75 75 ext_modules=[
76 76 Extension("schainpy.model.data._noise", ["schainc/_noise.c"]),
77 Extension("schainpy.model.data._HS_algorithm", ["schainc/_HS_algorithm.c"]),
77 78 ],
78 79 setup_requires = ["numpy"],
79 80 install_requires = [
General Comments 0
You need to be logged in to leave comments. Login now