##// END OF EJS Templates
schain
RSS
Tested window 20
Christianpl -
r1782:c73c874ee2f8
parent child
Show More
Side by Side Unified
Context file:
r1782:c73c874ee2f8 Loading diff...:
Show file before | Show file after | Hide comments
@@ -78,7 +78,7 def hildebrand_sekhon(data, navg):
78 78 sortdata = numpy.sort(data, axis=None)
79 79 #print(numpy.shape(data))
80 80 #exit()
81 '''
81
82 82 lenOfData = len(sortdata)
83 83 nums_min = lenOfData*0.2
84 84
@@ -110,8 +110,8 def hildebrand_sekhon(data, navg):
110 110 lnoise = sump / j
111 111
112 112 return lnoise
113 '''
114 return _noise.hildebrand_sekhon(sortdata, navg)
113
114 #return _noise.hildebrand_sekhon(sortdata, navg)
115 115
116 116
117 117 class Beam:
@@ -491,6 +491,7 class Spectra(JROData):
491 491 #exit(1)
492 492 daux = self.data_spc[channel,
493 493 xmin_index:xmax_index, ymin_index:ymax_index]
494 #print("daux",daux)
494 495 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
495 496
496 497 return noise
Show file before | Show file after | Hide comments
@@ -54,15 +54,14 class SpectraPlot(Plot):
54 54 #print("NormFactor: ",dataOut.normFactor)
55 55 #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
56 56 if hasattr(dataOut, 'LagPlot'): #Double Pulse
57 max_hei_id = dataOut.nHeights - 2*dataOut.LagPlot
58 #data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=46,ymax_index=max_hei_id)/dataOut.normFactor)
59 #data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=40,ymax_index=max_hei_id)/dataOut.normFactor)
60 data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=53,ymax_index=max_hei_id)/dataOut.normFactor)
61 data['noise'][0] = 10*numpy.log10(dataOut.getNoise(ymin_index=53)[0]/dataOut.normFactor)
57 ymin_index = numpy.abs(dataOut.heightList - 800).argmin()
58 max_hei_id = dataOut.nHeights - dataOut.TxLagRate*dataOut.LagPlot
59 data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=ymin_index,ymax_index=max_hei_id)/dataOut.normFactor)
60 data['noise'][0] = 10*numpy.log10(dataOut.getNoise(ymin_index=ymin_index)[0]/dataOut.normFactor)
62 61 #data['noise'][1] = 22.035507
63 62 else:
64 63 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
65 #data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=26,ymax_index=44)/dataOut.normFactor)
64
66 65 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
67 66
68 67 if self.CODE == 'spc_moments':
Show file before | Show file after | Hide comments
@@ -1008,7 +1008,7 class ACFsPlot(Plot):
1008 1008 BadHei2 = data['Height_contaminated_2']
1009 1009
1010 1010 self.xmin = 0.0
1011 self.xmax = 2.0
1011 #self.xmax = 2.0
1012 1012 self.y = ACFs
1013 1013
1014 1014 ax = self.axes[0]
Show file before | Show file after | Hide comments
@@ -171,8 +171,24
171 171 {"year": 2025, "doy": 24, "initial_time": [0,0], "final_time": [5,0], "aux_index": [ null, 13]},
172 172
173 173 {"year": 2025, "doy": 24, "initial_time": [5,0], "final_time": [23,55], "aux_index": [ null, 13]},
174 {"year": 2025, "doy": 25, "initial_time": [0,0], "final_time": [4,55], "aux_index": [ null, 27]},
174 {"year": 2025, "doy": 24, "initial_time": [5,0], "final_time": [5,5], "aux_index": [ null, 34]},
175 {"year": 2025, "doy": 24, "initial_time": [5,10], "final_time": [5,10], "aux_index": [ 24, 35]},
176 {"year": 2025, "doy": 24, "initial_time": [5,15], "final_time": [5,15], "aux_index": [ 25, 33]},
177 {"year": 2025, "doy": 24, "initial_time": [5,20], "final_time": [5,20], "aux_index": [ 26, 33]},
178 {"year": 2025, "doy": 24, "initial_time": [5,25], "final_time": [5,25], "aux_index": [ 27, 32]},
179 {"year": 2025, "doy": 24, "initial_time": [5,30], "final_time": [5,30], "aux_index": [ 23, 24]},
180 {"year": 2025, "doy": 24, "initial_time": [5,35], "final_time": [5,35], "aux_index": [ 22, 27]},
181 {"year": 2025, "doy": 24, "initial_time": [5,40], "final_time": [5,55], "aux_index": [ 22, 31]},
182 {"year": 2025, "doy": 24, "initial_time": [6,0], "final_time": [6,0], "aux_index": [ 22, 33]},
183 {"year": 2025, "doy": 24, "initial_time": [6,30], "final_time": [6,30], "aux_index": [ 22, 25]},
184 {"year": 2025, "doy": 24, "initial_time": [6,35], "final_time": [6,35], "aux_index": [ 22, 32]},
185 {"year": 2025, "doy": 24, "initial_time": [6,35], "final_time": [6,40], "aux_index": [ 22, 32]},
186 {"year": 2025, "doy": 24, "initial_time": [6,45], "final_time": [6,45], "aux_index": [ 23, 32]},
187 {"year": 2025, "doy": 24, "initial_time": [6,50], "final_time": [6,50], "aux_index": [ 24, 34]},
188 {"year": 2025, "doy": 24, "initial_time": [6,55], "final_time": [6,55], "aux_index": [ 25, 32]},
189 {"year": 2025, "doy": 25, "initial_time": [0,50], "final_time": [4,55], "aux_index": [ null, 27]},
175 190 {"year": 2025, "doy": 25, "initial_time": [2,50], "final_time": [4,55], "aux_index": [ null, 30]},
191 {"year": 2025, "doy": 25, "initial_time": [3,5], "final_time": [3,45], "aux_index": [ null, 33]},
176 192 {"year": 2025, "doy": 25, "initial_time": [0,0], "final_time": [5,0], "aux_index": [ null, 13]},
177 193
178 194 {"year": 2025, "doy": 25, "initial_time": [5,0], "final_time": [23,55], "aux_index": [ null, 13]},
@@ -183,9 +199,43
183 199 {"year": 2025, "doy": 26, "initial_time": [0,0], "final_time": [5,0], "aux_index": [ null, 13]},
184 200
185 201 {"year": 2025, "doy": 26, "initial_time": [5,0], "final_time": [23,55], "aux_index": [ null, 13]},
202 {"year": 2025, "doy": 26, "initial_time": [5,20], "final_time": [5,20], "aux_index": [ null, 28]},
203 {"year": 2025, "doy": 26, "initial_time": [5,25], "final_time": [5,25], "aux_index": [ 25, 34]},
186 204 {"year": 2025, "doy": 26, "initial_time": [0,55], "final_time": [0,55], "aux_index": [ null, null]},
187 205 {"year": 2025, "doy": 26, "initial_time": [0,55], "final_time": [4,55], "aux_index": [ null, 24]},
188 {"year": 2025, "doy": 27, "initial_time": [0,0], "final_time": [5,0], "aux_index": [ null, 13]}
206 {"year": 2025, "doy": 27, "initial_time": [0,0], "final_time": [5,0], "aux_index": [ null, 13]},
207
208 {"year": 2025, "doy": 41, "initial_time": [23,0], "final_time": [23,0], "aux_index": [ null, null]},
209 {"year": 2025, "doy": 42, "initial_time": [0,0], "final_time": [5,0], "aux_index": [ null, 14]},
210 {"year": 2025, "doy": 42, "initial_time": [4,24], "final_time": [4,24], "aux_index": [ 35, 38]},
211 {"year": 2025, "doy": 42, "initial_time": [1,30], "final_time": [5,0], "aux_index": [ null, 26]},
212 {"year": 2025, "doy": 42, "initial_time": [2,12], "final_time": [4,0], "aux_index": [ null, null]},
213 {"year": 2025, "doy": 42, "initial_time": [1,24], "final_time": [1,24], "aux_index": [ 39, 40]},
214 {"year": 2025, "doy": 42, "initial_time": [1,18], "final_time": [1,18], "aux_index": [ 48, 50]},
215
216 {"year": 2025, "doy": 42, "initial_time": [5,0], "final_time": [12,0], "aux_index": [ null, 7]},
217 {"year": 2025, "doy": 42, "initial_time": [5,0], "final_time": [7,30], "aux_index": [ null, 13]},
218 {"year": 2025, "doy": 42, "initial_time": [5,0], "final_time": [11,0], "aux_index": [ null, 10]},
219 {"year": 2025, "doy": 42, "initial_time": [5,0], "final_time": [5,0], "aux_index": [ 20, 22]},
220
221 {"year": 2025, "doy": 43, "initial_time": [0,0], "final_time": [5,0], "aux_index": [ null, 19]},
222 {"year": 2025, "doy": 43, "initial_time": [2,25], "final_time": [2,38], "aux_index": [ null, 24]},
223 {"year": 2025, "doy": 43, "initial_time": [2,2], "final_time": [2,2], "aux_index": [ null, null]},
224
225 {"year": 2025, "doy": 43, "initial_time": [5,0], "final_time": [10,50], "aux_index": [ null, 12]},
226 {"year": 2025, "doy": 43, "initial_time": [5,0], "final_time": [23,0], "aux_index": [ null, 8]},
227 {"year": 2025, "doy": 44, "initial_time": [0,50], "final_time": [0,50], "aux_index": [ 47, null]},
228 {"year": 2025, "doy": 44, "initial_time": [1,0], "final_time": [1,15], "aux_index": [ 39, null]},
229 {"year": 2025, "doy": 44, "initial_time": [3,14], "final_time": [3,40], "aux_index": [ null, 30]},
230 {"year": 2025, "doy": 44, "initial_time": [3,26], "final_time": [3,26], "aux_index": [ null, 32]},
231 {"year": 2025, "doy": 44, "initial_time": [0,14], "final_time": [3,40], "aux_index": [ null, 22]},
232 {"year": 2025, "doy": 44, "initial_time": [0,40], "final_time": [3,40], "aux_index": [ null, 25]},
233 {"year": 2025, "doy": 44, "initial_time": [4,0], "final_time": [5,0], "aux_index": [ null, null]},
234
235 {"year": 2025, "doy": 44, "initial_time": [5,50], "final_time": [5,50], "aux_index": [ null, null]},
236 {"year": 2025, "doy": 44, "initial_time": [7,50], "final_time": [8,38], "aux_index": [ null, null]},
237 {"year": 2025, "doy": 44, "initial_time": [5,0], "final_time": [9,14], "aux_index": [ null, 14]}
238
189 239
190 240 ]}
191 241
Show file before | Show file after | Hide comments
@@ -1827,7 +1827,7 class Oblique_Gauss_Fit(Operation):
1827 1827 #print("After data_snr: ", dataOut.data_snr)
1828 1828 dataOut.mode = mode
1829 1829 dataOut.flagNoData = numpy.all(numpy.isnan(dataOut.Dop_EEJ_T1)) #Si todos los valores son NaN no se prosigue
1830 dataOut.flagNoData = False #Descomentar solo para ploteo sino mantener comentado (para guardado)
1830 #dataOut.flagNoData = False #Descomentar solo para ploteo sino mantener comentado (para guardado)
1831 1831
1832 1832 return dataOut
1833 1833
Show file before | Show file after | Hide comments
@@ -84,6 +84,7 class SpectraLagProc(ProcessingUnit):
84 84 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
85 85 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
86 86 self.dataOut.runNextUnit = self.dataIn.runNextUnit
87 self.dataOut.TxLagRate = self.dataIn.TxLagRate
87 88 try:
88 89 self.dataOut.final_noise = self.dataIn.final_noise
89 90 except:
@@ -1403,8 +1404,6 class IntegrationFaradaySpectra(Operation):
1403 1404 for p in list(outliers_IDs):
1404 1405 buffer1[p,:] = avg
1405 1406
1406 else:
1407 print("No average!!")
1408 1407 '''Reasignate the buffer 1 edition to the original buffer'''
1409 1408 self.__buffer_spc[:,i,:,k,l]=numpy.copy(buffer1)
1410 1409 ###cspc IDs
@@ -2772,7 +2771,7 class SpectraDataToFaraday(Operation): #ISR MODE
2772 2771 def noise(self,dataOut):
2773 2772
2774 2773 dataOut.noise_lag = numpy.zeros((dataOut.nChannels,dataOut.DPL),'float32')
2775 #print("Lags")
2774
2776 2775 '''
2777 2776 for lag in range(dataOut.DPL):
2778 2777 #print(lag)
@@ -2780,43 +2779,57 class SpectraDataToFaraday(Operation): #ISR MODE
2780 2779 dataOut.noise_lag[:,lag] = dataOut.getNoise(ymin_index=46)
2781 2780 #dataOut.noise_lag[:,lag] = dataOut.getNoise(ymin_index=33,ymax_index=46)
2782 2781 '''
2783 #print(dataOut.NDP)
2784 #exit(1)
2782
2785 2783 #Channel B
2784
2785 ymin_index = numpy.abs(dataOut.heightList - 800).argmin() # Index more near to 800 km # ymin_index = 53
2786 2786 for lag in range(dataOut.DPL):
2787 #print(lag)
2787
2788 2788 dataOut.data_spc = dataOut.dataLag_spc[:,:,:,lag]
2789 max_hei_id = dataOut.NDP - 2*lag
2790 #if lag < 6:
2791 dataOut.noise_lag[1,lag] = dataOut.getNoise(ymin_index=53,ymax_index=max_hei_id)[1]
2792 #else:
2793 #dataOut.noise_lag[1,lag] = numpy.mean(dataOut.noise_lag[1,:6])
2794 #dataOut.noise_lag[:,lag] = dataOut.getNoise(ymin_index=33,ymax_index=46)
2789 max_hei_id = dataOut.NDP - dataOut.TxLagRate*lag # - 2*lag #max_hei_id = 32
2790 if max_hei_id < ymin_index: break_lag = lag; continue
2791
2792 #dataOut.noise_lag[1,lag] = dataOut.getNoise(ymin_index=ymin_index,ymax_index=max_hei_id)[1] #H_S algorithm choise
2793
2794 daux = dataOut.data_spc[1,:, ymin_index:max_hei_id] # Median choise
2795 sortdata = numpy.sort(daux, axis=None)
2796 dataOut.noise_lag[1,lag] = numpy.median(sortdata)
2797
2798 dataOut.noise_lag[1,break_lag:] = numpy.mean(dataOut.noise_lag[1,:break_lag]) #El ruido de lags que no pueden calcularse
2799 #su piso de ruido se determina a partir
2800 # del promedio de los lags limpios
2801
2795 2802 #Channel A
2803
2796 2804 for lag in range(dataOut.DPL):
2797 #print(lag)
2798 2805 dataOut.data_spc = dataOut.dataLag_spc[:,:,:,lag]
2799 dataOut.noise_lag[0,lag] = dataOut.getNoise(ymin_index=53)[0]
2806 #dataOut.noise_lag[0,lag] = dataOut.getNoise(ymin_index=ymin_index)[0]
2807 daux = dataOut.data_spc[0,:, ymin_index:]
2808 sortdata = numpy.sort(daux, axis=None)
2809 dataOut.noise_lag[0,lag] = numpy.median(sortdata)
2800 2810
2811 print("dataOut.noise_lag", dataOut.noise_lag[1,:])
2801 2812 nanindex = numpy.argwhere(numpy.isnan(numpy.log10(dataOut.noise_lag[1,:])))
2813 print(nanindex)
2814
2815 try:
2816 i1 = nanindex[0][0]
2817 dataOut.noise_lag[1,i1:] = numpy.mean(dataOut.noise_lag[1,:i1])
2818 except:
2819 pass
2820
2821 '''try:
2802 2822 i1 = nanindex[0][0]
2803 2823 dataOut.noise_lag[1,i1:] = numpy.mean(dataOut.noise_lag[1,:i1]) #El ruido de lags contaminados se
2804 2824 #determina a partir del promedio del
2805 2825 #ruido de los lags limpios
2806 '''
2826
2827 except:
2807 2828 dataOut.noise_lag[1,:] = dataOut.noise_lag[1,0] #El ruido de los lags diferentes de cero para
2808 2829 #el canal B es contaminado por el Tx y EEJ
2809 2830 #del siguiente perfil, por ello se asigna el ruido
2810 #del lag 0 a todos los lags
2811 '''
2812 #print("Noise lag: ", 10*numpy.log10(dataOut.noise_lag/dataOut.normFactor))
2813 #exit(1)
2814 '''
2815 dataOut.tnoise = dataOut.getNoise(ymin_index=46)
2816 dataOut.tnoise /= float(dataOut.nProfiles*dataOut.nIncohInt)
2817 dataOut.pan = dataOut.tnoise[0]
2818 dataOut.pbn = dataOut.tnoise[1]
2819 '''
2831 #del lag 0 a todos los lags'''
2832
2820 2833
2821 2834 dataOut.tnoise = dataOut.noise_lag/float(dataOut.nProfiles*dataOut.nIncohInt)
2822 2835 #dataOut.tnoise /= float(dataOut.nProfiles*dataOut.nIncohInt)
@@ -2892,11 +2905,12 class SpectraDataToFaraday(Operation): #ISR MODE
2892 2905
2893 2906 def get_eej_index(self,data_to_remov_eej,dataOut):
2894 2907
2895 dataOut.data_spc = data_to_remov_eej
2908 max_eej = numpy.abs(dataOut.heightList - 255).argmin()
2896 2909
2910 dataOut.data_spc = data_to_remov_eej
2897 2911 data_eej = dataOut.getPower()[0]
2898 2912 #print(data_eej)
2899 index_eej = CleanCohEchoes.mad_based_outlier(self,data_eej[:17])
2913 index_eej = CleanCohEchoes.mad_based_outlier(self,data_eej[:max_eej])
2900 2914 aux_eej = numpy.array(index_eej.nonzero()).ravel()
2901 2915 print("aux_eej: ", aux_eej)
2902 2916 if aux_eej != []:
Show file before | Show file after | Hide comments
@@ -647,25 +647,44 class LagsReshape(Operation):
647 647 self.buffer_HRonelag = numpy.zeros((int(dataOut.NSCAN/dataOut.DPL),
648 648 dataOut.nHeights),
649 649 dtype='complex')
650
651 for i in range(self.buffer_HRonelag.shape[0]):
652 for j in range(dataOut.nHeights):
653 if j+int(2*whichlag)<dataOut.nHeights:
654 self.buffer_HRonelag[i,j]=dataOut.datalags[1,i,j+2*whichlag,whichlag]
650 TxLagRate = dataOut.TxLagRate
651 for i in range(self.buffer_HRonelag.shape[0]): #perfil
652 for j in range(dataOut.nHeights): # height
653 if j+int(TxLagRate*whichlag)<dataOut.nHeights:
654 self.buffer_HRonelag[i,j]=dataOut.datalags[1,i,j+TxLagRate*whichlag,whichlag]
655 655 else:
656 656 if whichlag!=10:
657 self.buffer_HRonelag[i,j]=dataOut.datalags[1,i,(j+2*whichlag)%dataOut.nHeights,whichlag+1]
657 self.buffer_HRonelag[i,j]=dataOut.datalags[1,i,(j+TxLagRate*whichlag)%dataOut.nHeights,whichlag+1]
658 658 else:
659 659 if i+2<self.buffer_HRonelag.shape[0]:
660 self.buffer_HRonelag[i,j]=dataOut.datalags[1,i+2,(j+2*whichlag)%dataOut.nHeights,0]
660 self.buffer_HRonelag[i,j]=dataOut.datalags[1,i+2,(j+TxLagRate*whichlag)%dataOut.nHeights,0]
661 661 else: #i+1==self.buffer_HRonelag.shape[0]:
662 self.buffer_HRonelag[i,j]=dataOut.datalags[1,i,(j+2*whichlag)%dataOut.nHeights,whichlag]
662 self.buffer_HRonelag[i,j]=dataOut.datalags[1,i,(j+TxLagRate*whichlag)%dataOut.nHeights,whichlag] #1, 198,64 = 1,198, 0, 10
663 663
664 664 return self.buffer_HRonelag
665 665
666 666
667 667
668 def run(self,dataOut,DPL=11,NSCAN=132):
668 def run(self,dataOut,DPL=11,NSCAN=132, TxLagRate=2):
669 dataOut.TxLagRate = TxLagRate
670 '''PA = dataOut.data[0,:,:]
671 PB = dataOut.data[1,:,:]
672 import matplotlib.pyplot as plt
673 fig, axes = plt.subplots(2, 11, figsize=(18, 6), sharex=True, sharey=True)
674
675 for i in range(11):
676 axes[0,i].plot(PA[i, :], dataOut.heightList, label=f'PA {i+1}')
677 axes[0, i].set_title(f'Lag {i+1}')
678 #axes[0, i].set_xscale("log") # Log scale for y-axis
679 #axes[0, i].set_xlim([0,1e+7])
680 axes[1,i].plot(PB[i, :], dataOut.heightList, label=f'PB {i+1}')
681 #axes[1, i].set_xscale("log") # Log scale for y-axis
682 #axes[1, i].set_xlim([0,1e+7])
683
684
685 plt.tight_layout()
686 plt.show()'''
687
669 688
670 689 dataOut.DPL=DPL
671 690 dataOut.NSCAN=NSCAN
@@ -677,6 +696,25 class LagsReshape(Operation):
677 696 dataOut.datalags=numpy.copy(self.LagDistribution(dataOut))
678 697 dataOut.datalags[1,:,:,:]=self.HeightReconstruction(dataOut)
679 698
699 '''
700 PA = dataOut.datalags[0,5,:,:]
701 PB = dataOut.datalags[1,5,:,:]
702 import matplotlib.pyplot as plt
703 fig, axes = plt.subplots(2, 11, figsize=(18, 6), sharex=True, sharey=True)
704
705 for i in range(11):
706 axes[0,i].plot(PA[:, i], dataOut.heightList, label=f'PA {i+1}')
707 axes[0, i].set_title(f'Lag {i+1}')
708 #axes[0, i].set_xscale("log") # Log scale for y-axis
709 #axes[0, i].set_xlim([0,1e+7])
710 axes[1,i].plot(PB[:, i], dataOut.heightList, label=f'PB {i+1}')
711 #axes[1, i].set_xscale("log") # Log scale for y-axis
712 #axes[1, i].set_xlim([0,1e+7])
713
714
715 plt.tight_layout()
716 plt.show()#'''
717
680 718 return dataOut
681 719
682 720
@@ -1351,7 +1389,7 class FlagBadHeights(Operation):
1351 1389 dataOut.ibad[j][l]=1
1352 1390 else:
1353 1391 dataOut.ibad[j][l]=0
1354
1392 #print("dataOut.ibad",dataOut.ibad)
1355 1393 return dataOut
1356 1394
1357 1395 class FlagBadHeightsSpectra(Operation):
@@ -2074,7 +2112,7 class DoublePulseACFs_PerLag(Operation):
2074 2112 # Stores lags for which ACFs are calculated
2075 2113 dataOut.alag=numpy.zeros(dataOut.NDP,'float32')
2076 2114 for l in range(dataOut.DPL):
2077 dataOut.alag[l]=l*dataOut.DH*2.0/150.0
2115 dataOut.alag[l]=l*dataOut.DH*dataOut.TxLagRate/150.0
2078 2116 self.aux=0
2079 2117 # dataOut.pan.- Power noise level of channel A - definned in SpectraDataToFaraday
2080 2118 # Signal noise
@@ -2085,6 +2123,9 class DoublePulseACFs_PerLag(Operation):
2085 2123
2086 2124 id = numpy.where(dataOut.heightList>700)[0]
2087 2125
2126 PA = numpy.zeros((dataOut.NDP, dataOut.DPL), dtype=float)
2127 PB = numpy.zeros((dataOut.NDP, dataOut.DPL), dtype=float)
2128
2088 2129 for i in range(dataOut.NDP): #Heights
2089 2130 for j in range(dataOut.DPL): # Lags
2090 2131 ################# Total power
@@ -2104,7 +2145,7 class DoublePulseACFs_PerLag(Operation):
2104 2145 ## ACF
2105 2146 rhorp=dataOut.kabxys_integrated[8][i,j,0]+dataOut.kabxys_integrated[11][i,j,0]
2106 2147 rhoip=dataOut.kabxys_integrated[10][i,j,0]-dataOut.kabxys_integrated[9][i,j,0]
2107
2148 #PA[i,j] = pa; PB[i,j] = pb
2108 2149 if ((pa>dataOut.pan[j])&(pb>dataOut.pbn[j])):
2109 2150 # panrm is RMS of power, used to normalize ACFs
2110 2151 ss4=numpy.abs((pa-dataOut.pan[j])*(pb-dataOut.pbn[j]))
@@ -2160,11 +2201,25 class DoublePulseACFs_PerLag(Operation):
2160 2201 #print("EJJ")
2161 2202 dataOut.igcej[i,j]=1
2162 2203 #'''
2163 '''
2164 if i == 4:
2165 exit(1)
2166 '''
2167 print("dataOut.alag",dataOut.alag)
2204 '''import matplotlib.pyplot as plt
2205 fig, axes = plt.subplots(2, dataOut.DPL, figsize=(18, 6), sharex=True, sharey=True)
2206
2207 for i in range(dataOut.DPL):
2208 axes[0,i].plot(PA[:, i], dataOut.heightList, label=f'PA {i+1}')
2209 axes[0, i].axvline(dataOut.pan[i], color='gray', linestyle='--', linewidth=1)
2210 axes[0, i].set_title(f'Lag {i+1}')
2211 axes[0, i].set_xscale("log") # Log scale for y-axis
2212 axes[0, i].set_xlim([0,1e+7])
2213 axes[1,i].plot(PB[:, i], dataOut.heightList, label=f'PB {i+1}')
2214 axes[1, i].axvline(dataOut.pbn[i], color='gray', linestyle='--', linewidth=1)
2215 axes[1, i].set_xscale("log") # Log scale for y-axis
2216 axes[1, i].set_xlim([0,1e+7])
2217
2218
2219 plt.tight_layout()
2220 plt.show()'''
2221
2222
2168 2223 #print("dataOut.p",datetime.datetime.utcfromtimestamp(dataOut.utctime), dataOut.p)
2169 2224
2170 2225 #print(numpy.sum(dataOut.kabxys_integrated[8][:,:,0]+dataOut.kabxys_integrated[11][:,:,0]))
@@ -2285,19 +2340,20 class FaradayAngleAndDPPower(Operation):
2285 2340 dataOut.flagTeTiCorrection = False
2286 2341 #print("ph2: ", numpy.sum(dataOut.ph2[:16]))
2287 2342 #print("ph2: ", numpy.sum(dataOut.ph2[16:32]))
2288 '''
2289 import matplotlib.pyplot as plt
2343
2344 '''import matplotlib.pyplot as plt
2290 2345 #plt.plot(numpy.abs(dataOut.kabxys_integrated[4][:,j,0]+dataOut.kabxys_integrated[5][:,j,0])+numpy.abs(dataOut.kabxys_integrated[6][:,j,0]+dataOut.kabxys_integrated[7][:,j,0]),dataOut.heightList)
2291 2346 #plt.axvline((dataOut.pan+dataOut.pbn))
2292 2347 #print(numpy.shape(dataOut.p))
2293 2348 plt.plot(dataOut.ph2,dataOut.heightList)
2349 plt.plot(dataOut.phi,dataOut.heightList)
2294 2350
2295 2351 plt.xlim(1000,1000000000)
2296 2352 #plt.ylim(50,400)
2297 2353 plt.grid()
2298 2354 plt.show()
2299 #exit(1)
2300 '''
2355 #exit(1)'''
2356
2301 2357 return dataOut
2302 2358
2303 2359 class ElectronDensityFaraday(Operation):
@@ -2790,6 +2846,35 class DPTemperaturesEstimation(Operation):
2790 2846 op.addParameter(name='IBITS', value='16', format='int')
2791 2847
2792 2848 """
2849 '''
2850 NSTHS Number of sample heights (input), for temperature processing
2851 NDP Number of Data Points (nHeights dependent)
2852 NSTHS < NDP
2853
2854 Input/Output Data:
2855
2856 te2, ti2: Estimated electron and ion temperatures.
2857
2858 ete2, eti2: Errors in the estimated temperatures.
2859
2860 phy2, ephy2: Physical parameter and its error.
2861
2862 Fitting Process:
2863
2864 ifit: Flags for which parameters are being fitted.
2865
2866 params: Initial guesses and fitted parameters.
2867
2868 cov, covinv: Covariance matrix and its inverse for uncertainty estimation.
2869
2870 Metadata/Status:
2871
2872 m: Status or counter for the fitting process.
2873
2874 info2: Additional information about the fitting results.
2875
2876 '''
2877
2793 2878
2794 2879 def __init__(self, **kwargs):
2795 2880
@@ -2819,8 +2904,8 class DPTemperaturesEstimation(Operation):
2819 2904
2820 2905 #null_fd = os.open(os.devnull, os.O_RDWR)
2821 2906 #os.dup2(null_fd, 1)
2822
2823 for i in range(10,dataOut.NSHTS): #no point below 150 km
2907 ymin_index = numpy.abs(dataOut.heightList - 150).argmin() #no point below 150 km
2908 for i in range(ymin_index,dataOut.NSHTS):
2824 2909
2825 2910 #some definitions
2826 2911 iflag=0 # inicializado a cero?
@@ -2834,7 +2919,11 class DPTemperaturesEstimation(Operation):
2834 2919 t1=numpy.zeros(1,order='F',dtype='float32')
2835 2920 t2=numpy.zeros(1,order='F',dtype='float32')
2836 2921
2837 if i>10 and l1>=0:
2922 '''
2923 x lag time y correlation e their uncertanities
2924 t1 t2 initial guesses eb errors
2925 '''
2926 if i>ymin_index and l1>=0:
2838 2927 if l1==0:
2839 2928 l1=1
2840 2929
@@ -2864,7 +2953,7 class DPTemperaturesEstimation(Operation):
2864 2953 for l in range(0+1,dataOut.DPL):
2865 2954 if dataOut.igcej[i][l]==0 and dataOut.ibad[i][l]==0:
2866 2955 y[l1]=dataOut.rhor[i][l]*cc + dataOut.rhoi[i][l]*ss
2867 x[l1]=dataOut.alag[l]*1.0e-3
2956 x[l1]=dataOut.alag[l]*1.0e-3 # *1.0e-3
2868 2957 dataOut.sd[i][l]=dataOut.sd[i][l]/((acfm)**2)# important
2869 2958 e[l1]=dataOut.sd[i][l] #this is the variance, not the st. dev.
2870 2959 l1=l1+1
@@ -2890,8 +2979,8 class DPTemperaturesEstimation(Operation):
2890 2979
2891 2980 if True: #len(y)!=0:
2892 2981 with suppress_stdout_stderr():
2893 fitacf_guess.guess(y,x,zero,depth,t1,t2,len(y))
2894 t2=t1/t2
2982 fitacf_guess.guess(y,x,zero,depth,t1,t2,len(y)) #t1 = te , t2 = tr = te/ti
2983 t2=t1/t2 # ti
2895 2984
2896 2985 if (t1<5000.0 and t1> 600.0):
2897 2986 dataOut.params[1]=t1
@@ -2899,7 +2988,7 class DPTemperaturesEstimation(Operation):
2899 2988 dataOut.ifit[1]=dataOut.ifit[2]=1
2900 2989 dataOut.ifit[0]=dataOut.ifit[3]=dataOut.ifit[4]=0
2901 2990
2902 if dataOut.ut_Faraday<10.0 and dataOut.ut_Faraday>=0.5:
2991 if dataOut.ut_Faraday<10.0 and dataOut.ut_Faraday>=0.5: # 6 30 pm to 5 am LT
2903 2992 dataOut.ifit[2]=0
2904 2993
2905 2994 den=dataOut.ph2[i]
@@ -3514,7 +3603,7 class DataSaveCleaner(Operation):
3514 3603 dataOut.PhyFinal[0,25:]=missing
3515 3604 dataOut.EPhyFinal[0, 25:] = missing
3516 3605 '''
3517 # 8 Sep 24
3606 '''# 8 Sep 24
3518 3607 if True: #06-18 LT
3519 3608 #dataOut.DensityFinal[0,27:]=missing
3520 3609 #dataOut.EDensityFinal[0,27:]=missing
@@ -3523,7 +3612,17 class DataSaveCleaner(Operation):
3523 3612 dataOut.IonTempFinal[0,36:]=missing
3524 3613 dataOut.EIonTempFinal[0,36:]=missing
3525 3614 dataOut.PhyFinal[0,36:]=missing
3526 dataOut.EPhyFinal[0, 36:] = missing
3615 dataOut.EPhyFinal[0, 36:] = missing'''
3616 '''# 24 Jan 25
3617 if (time_text.hour >= 5 ) and (time_text.hour <= 7):
3618 #dataOut.DensityFinal[0,27:]=missing
3619 #dataOut.EDensityFinal[0,27:]=missing
3620 dataOut.ElecTempFinal[0,:]=missing
3621 dataOut.EElecTempFinal[0,:]=missing
3622 dataOut.IonTempFinal[0,:]=missing
3623 dataOut.EIonTempFinal[0,:]=missing
3624 dataOut.PhyFinal[0,:]=missing
3625 dataOut.EPhyFinal[0, :] = missing'''
3527 3626 start = time()
3528 3627 flagcleandata = True
3529 3628 if flagcleandata:
@@ -3636,7 +3735,7 class DataSaveCleaner(Operation):
3636 3735 f.writerow(cf)
3637 3736 file.close()
3638 3737 # for plot
3639 dataOut.flagNoData = False #Descomentar solo para ploteo #Comentar para MADWriter
3738 #dataOut.flagNoData = False #Descomentar solo para ploteo #Comentar para MADWriter
3640 3739
3641 3740 dataOut.DensityFinal *= 1.e6 #Convert units to m^⁻3
3642 3741 dataOut.EDensityFinal *= 1.e6 #Convert units to m^⁻3
@@ -3738,7 +3837,8 class ACFs(Operation):
3738 3837 dataOut.y_ibad_to_plot[i,:]=numpy.copy(self.y_ibad)
3739 3838
3740 3839 missing=numpy.nan#-32767
3741
3840 #print("dataOut.igcej",dataOut.igcej)
3841 #print("dataOut.ibad",dataOut.ibad)
3742 3842 for i in range(dataOut.NSHTS,dataOut.NDP):
3743 3843 for j in range(dataOut.DPL):
3744 3844 dataOut.acfs_to_save[i,j]=missing
Show file before | Show file after | Hide comments
@@ -1804,188 +1804,6 class Integration(Operation):
1804 1804
1805 1805
1806 1806
1807
1808
1809
1810 class SumLagProducts_Old(Operation):
1811 def __init__(self, **kwargs):
1812
1813 Operation.__init__(self, **kwargs)
1814 #dataOut.rnint2=numpy.zeros(dataOut.nlags_array,'float32')
1815
1816
1817 def run(self,dataOut):
1818
1819 if dataOut.AUX: #Solo cuando ya hizo la intregacion se ejecuta
1820
1821
1822 dataOut.rnint2=numpy.zeros(dataOut.header[17][0],'float32')
1823 #print(dataOut.experiment)
1824 if dataOut.experiment=="DP":
1825 for l in range(dataOut.header[17][0]):
1826 dataOut.rnint2[l]=1.0/(dataOut.nint*dataOut.header[7][0]*12.0)
1827
1828
1829 if dataOut.experiment=="HP":
1830 for l in range(dataOut.header[17][0]):
1831 if(l==0 or (l>=3 and l <=6)):
1832 dataOut.rnint2[l]=0.5/float(dataOut.nint*dataOut.header[7][0]*16.0)
1833 else:
1834 dataOut.rnint2[l]=0.5/float(dataOut.nint*dataOut.header[7][0]*8.0)
1835 #print(dataOut.rnint2)
1836 for l in range(dataOut.header[17][0]):
1837
1838 dataOut.kabxys_integrated[4][:,l,0]=(dataOut.kabxys_integrated[4][:,l,0]+dataOut.kabxys_integrated[4][:,l,1])*dataOut.rnint2[l]
1839 dataOut.kabxys_integrated[5][:,l,0]=(dataOut.kabxys_integrated[5][:,l,0]+dataOut.kabxys_integrated[5][:,l,1])*dataOut.rnint2[l]
1840 dataOut.kabxys_integrated[6][:,l,0]=(dataOut.kabxys_integrated[6][:,l,0]+dataOut.kabxys_integrated[6][:,l,1])*dataOut.rnint2[l]
1841 dataOut.kabxys_integrated[7][:,l,0]=(dataOut.kabxys_integrated[7][:,l,0]+dataOut.kabxys_integrated[7][:,l,1])*dataOut.rnint2[l]
1842
1843 dataOut.kabxys_integrated[8][:,l,0]=(dataOut.kabxys_integrated[8][:,l,0]-dataOut.kabxys_integrated[8][:,l,1])*dataOut.rnint2[l]
1844 dataOut.kabxys_integrated[9][:,l,0]=(dataOut.kabxys_integrated[9][:,l,0]-dataOut.kabxys_integrated[9][:,l,1])*dataOut.rnint2[l]
1845 dataOut.kabxys_integrated[10][:,l,0]=(dataOut.kabxys_integrated[10][:,l,0]-dataOut.kabxys_integrated[10][:,l,1])*dataOut.rnint2[l]
1846 dataOut.kabxys_integrated[11][:,l,0]=(dataOut.kabxys_integrated[11][:,l,0]-dataOut.kabxys_integrated[11][:,l,1])*dataOut.rnint2[l]
1847
1848
1849 #print("Final Integration: ",dataOut.kabxys_integrated[4][:,l,0])
1850
1851
1852
1853
1854
1855
1856 return dataOut
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866 class BadHeights_Old(Operation):
1867 def __init__(self, **kwargs):
1868
1869 Operation.__init__(self, **kwargs)
1870
1871
1872
1873 def run(self,dataOut):
1874
1875
1876 if dataOut.AUX==1:
1877 dataOut.ibad=numpy.zeros((dataOut.header[15][0],dataOut.header[17][0]),'int32')
1878
1879 for j in range(dataOut.header[15][0]):
1880 for l in range(dataOut.header[17][0]):
1881 ip1=j+dataOut.header[15][0]*(0+2*l)
1882 if( (dataOut.kabxys_integrated[5][j,l,0] <= 0.) or (dataOut.kabxys_integrated[4][j,l,0] <= 0.) or (dataOut.kabxys_integrated[7][j,l,0] <= 0.) or (dataOut.kabxys_integrated[6][j,l,0] <= 0.)):
1883 dataOut.ibad[j][l]=1
1884 else:
1885 dataOut.ibad[j][l]=1
1886 #print("ibad: ",dataOut.ibad)
1887
1888
1889
1890 return dataOut
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907 class NoisePower_old(Operation):
1908 def __init__(self, **kwargs):
1909
1910 Operation.__init__(self, **kwargs)
1911
1912 def hildebrand(self,dataOut,data):
1913 #print("data ",data )
1914 divider=10 # divider was originally 10
1915 noise=0.0
1916 n1=0
1917 n2=int(dataOut.header[15][0]/2)
1918 sorts= sorted(data)
1919
1920 nums_min= dataOut.header[15][0]/divider
1921 if((dataOut.header[15][0]/divider)> 2):
1922 nums_min= int(dataOut.header[15][0]/divider)
1923 else:
1924 nums_min=2
1925 sump=0.0
1926 sumq=0.0
1927 j=0
1928 cont=1
1929 while( (cont==1) and (j<n2)):
1930 sump+=sorts[j+n1]
1931 sumq+= sorts[j+n1]*sorts[j+n1]
1932 t3= sump/(j+1)
1933 j=j+1
1934 if(j> nums_min):
1935 rtest= float(j/(j-1)) +1.0/dataOut.header[7][0]
1936 t1= (sumq*j)
1937 t2=(rtest*sump*sump)
1938 if( (t1/t2) > 0.990):
1939 j=j-1
1940 sump-= sorts[j+n1]
1941 sumq-=sorts[j+n1]*sorts[j+n1]
1942 cont= 0
1943
1944 noise= sump/j
1945 stdv=numpy.sqrt((sumq- noise*noise)/(j-1))
1946 return noise
1947
1948 def run(self,dataOut):
1949
1950 if dataOut.AUX==1:
1951
1952 #print("ax2 shape ",ax2.shape)
1953 p=numpy.zeros((dataOut.header[2][0],dataOut.header[15][0],dataOut.header[17][0]),'float32')
1954 av=numpy.zeros(dataOut.header[15][0],'float32')
1955 dataOut.pnoise=numpy.zeros(dataOut.header[2][0],'float32')
1956
1957 p[0,:,:]=dataOut.kabxys_integrated[4][:,:,0]+dataOut.kabxys_integrated[5][:,:,0] #total power for channel 0, just pulse with non-flip
1958 p[1,:,:]=dataOut.kabxys_integrated[6][:,:,0]+dataOut.kabxys_integrated[7][:,:,0] #total power for channel 1
1959
1960 #print("p[0,:,:] ",p[0,:,:])
1961 #print("p[1,:,:] ",p[1,:,:])
1962
1963 for i in range(dataOut.header[2][0]):
1964 dataOut.pnoise[i]=0.0
1965 for k in range(dataOut.header[17][0]):
1966 dataOut.pnoise[i]+= self.hildebrand(dataOut,p[i,:,k])
1967 #print("dpl ",k, "pnoise[",i,"] ",pnoise[i] )
1968 dataOut.pnoise[i]=dataOut.pnoise[i]/dataOut.header[17][0]
1969
1970
1971 #print("POWERNOISE: ",dataOut.pnoise)
1972 dataOut.pan=1.0*dataOut.pnoise[0] # weights could change
1973 dataOut.pbn=1.0*dataOut.pnoise[1] # weights could change
1974 #print("dataOut.pan ",dataOut.pan, " dataOut.pbn ",dataOut.pbn)
1975 #print("AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAa")
1976
1977 #print("POWERNOISE: ",dataOut.pnoise)
1978
1979
1980 return dataOut
1981
1982
1983
1984
1985
1986
1987
1988
1989 1807 class double_pulse_ACFs(Operation):
1990 1808 def __init__(self, **kwargs):
1991 1809
@@ -2025,8 +1843,8 class double_pulse_ACFs(Operation):
2025 1843 dataOut.p[i,j]=pa+pb-(dataOut.pan+dataOut.pbn)
2026 1844 dataOut.sdp[i,j]=2*dataOut.rnint2[j]*((pa+pb)*(pa+pb))
2027 1845 ## ACF
2028 rhorp=dataOut.kabxys_integrated[8][i,j,0]+dataOut.kabxys_integrated[11][i,j,0]
2029 rhoip=dataOut.kabxys_integrated[10][i,j,0]-dataOut.kabxys_integrated[9][i,j,0]
1846 rhorp=dataOut.kabxys_integrated[8][i,j,0]+dataOut.kabxys_integrated[11][i,j,0] # cspc real +
1847 rhoip=dataOut.kabxys_integrated[10][i,j,0]-dataOut.kabxys_integrated[9][i,j,0] # cspc imag +
2030 1848 if ((pa>dataOut.pan)&(pb>dataOut.pbn)):
2031 1849 #print("dataOut.pnoise[0]: ",dataOut.pnoise[0])
2032 1850 #print("dataOut.pnoise[1]: ",dataOut.pnoise[1])
Show file before | Show file after | Hide comments
@@ -468,12 +468,48
468 468 "time": [[4,0]],
469 469 "cf": "dataOut.cflast[0]"},
470 470
471 {"year": 2025,"doy": 24,
472 "time": [[6,50]],
473 "cf": "dataOut.cflast[0]"},
471 474 {"year": 2025,"doy": 25,
472 475 "time": [[3,5],[3,30]],
473 476 "cf": "dataOut.cflast[0]"},
474 477
475 478 {"year": 2025,"doy": 26,
476 479 "time": [[4,15],[4,20]],
480 "cf": "dataOut.cflast[0]"},
481
482 {"year": 2025,"doy": 26,
483 "time": [[8,0],[10,20],[10,30]],
484 "cf": "dataOut.cflast[0]"},
485
486 {"year": 2025,"doy": 42,
487 "time": [[4,24]],
488 "cf": 0.0076},
489
490 {"year": 2025,"doy": 42,
491 "time": [[10,30]],
492 "cf": "dataOut.cflast[0]"},
493
494 {"year": 2025,"doy": 43,
495 "time": [[2,26]],
496 "cf": 0.008362},
497 {"year": 2025,"doy": 43,
498 "time": [[3,26]],
499 "cf": 0.007769},
500 {"year": 2025,"doy": 43,
501 "time": [[4,26]],
502 "cf": 0.004193},
503 {"year": 2025,"doy": 43,
504 "time": [[4,38]],
505 "cf": 0.005492},
506
507 {"year": 2025,"doy": 43,
508 "time": [[7,2],[7,14]],
509 "cf": "dataOut.cflast[0]"},
510
511 {"year": 2025,"doy": 44,
512 "time": [[6,26],[6,38],[6,50],[7,2]],
477 513 "cf": "dataOut.cflast[0]"}
478 514
479 515
General Comments 0
You need to be logged in to leave comments. Login now