##// END OF EJS Templates
schain
RSS
DP+LP Join Spc+Voltage
rflores -
r1549:c70288c93a07
parent child
Show More
Side by Side Unified
Context file:
r1549:c70288c93a07 Loading diff...:
Show file before | Show file after | Hide comments
@@ -260,7 +260,7 class JROData(GenericData):
260
260
261 def getFmax(self):
261 def getFmax(self):
262 PRF = 1. / (self.ippSeconds * self.nCohInt)
262 PRF = 1. / (self.ippSeconds * self.nCohInt)
263
263 #print("ippsec",self.ippSeconds)
264 fmax = PRF
264 fmax = PRF
265 return fmax
265 return fmax
266
266
Show file before | Show file after | Hide comments
@@ -256,6 +256,7 class Plot(Operation):
256 self.__throttle_plot = apply_throttle(self.throttle)
256 self.__throttle_plot = apply_throttle(self.throttle)
257 code = self.attr_data if self.attr_data else self.CODE
257 code = self.attr_data if self.attr_data else self.CODE
258 self.data = PlotterData(self.CODE, self.exp_code, self.localtime)
258 self.data = PlotterData(self.CODE, self.exp_code, self.localtime)
259 #self.EEJtype = kwargs.get('EEJtype', 2)
259
260
260 if self.server:
261 if self.server:
261 if not self.server.startswith('tcp://'):
262 if not self.server.startswith('tcp://'):
Show file before | Show file after | Hide comments
@@ -75,7 +75,7 class SnrPlot(RTIPlot):
75 def update(self, dataOut):
75 def update(self, dataOut):
76
76
77 data = {
77 data = {
78 'snr': 10*numpy.log10(dataOut.data_snr)
78 'snr': 10*numpy.log10(dataOut.data_snr)
79 }
79 }
80
80
81 return data, {}
81 return data, {}
@@ -91,7 +91,120 class DopplerPlot(RTIPlot):
91 def update(self, dataOut):
91 def update(self, dataOut):
92
92
93 data = {
93 data = {
94 'dop': 10*numpy.log10(dataOut.data_dop)
94 'dop': 10*numpy.log10(dataOut.data_dop)
95 }
96
97 return data, {}
98
99 class DopplerEEJPlot_V0(RTIPlot):
100 '''
101 Written by R. Flores
102 '''
103 '''
104 Plot for EEJ
105 '''
106
107 CODE = 'dop'
108 colormap = 'RdBu_r'
109 colormap = 'jet'
110
111 def setup(self):
112
113 self.xaxis = 'time'
114 self.ncols = 1
115 self.nrows = len(self.data.channels)
116 self.nplots = len(self.data.channels)
117 self.ylabel = 'Range [km]'
118 self.xlabel = 'Time'
119 self.cb_label = '(m/s)'
120 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
121 self.titles = ['{} Channel {}'.format(
122 self.CODE.upper(), x) for x in range(self.nrows)]
123
124 def update(self, dataOut):
125 #print(self.EEJtype)
126
127 if self.EEJtype == 1:
128 data = {
129 'dop': dataOut.Oblique_params[:,-2,:]
130 }
131 elif self.EEJtype == 2:
132 data = {
133 'dop': dataOut.Oblique_params[:,-1,:]
134 }
135
136 return data, {}
137
138 class DopplerEEJPlot(RTIPlot):
139 '''
140 Written by R. Flores
141 '''
142 '''
143 Plot for Doppler Shift EEJ
144 '''
145
146 CODE = 'dop'
147 colormap = 'RdBu_r'
148 #colormap = 'jet'
149
150 def setup(self):
151
152 self.xaxis = 'time'
153 self.ncols = 1
154 self.nrows = 2
155 self.nplots = 2
156 self.ylabel = 'Range [km]'
157 self.xlabel = 'Time'
158 self.cb_label = '(m/s)'
159 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
160 self.titles = ['{} EJJ Type {} /'.format(
161 self.CODE.upper(), x) for x in range(1,1+self.nrows)]
162
163 def update(self, dataOut):
164
165 if dataOut.mode == 11: #Double Gaussian
166 doppler = numpy.append(dataOut.Oblique_params[:,1,:],dataOut.Oblique_params[:,4,:],axis=0)
167 elif dataOut.mode == 9: #Double Skew Gaussian
168 doppler = numpy.append(dataOut.Oblique_params[:,-2,:],dataOut.Oblique_params[:,-1,:],axis=0)
169 data = {
170 'dop': doppler
171 }
172
173 return data, {}
174
175 class SpcWidthEEJPlot(RTIPlot):
176 '''
177 Written by R. Flores
178 '''
179 '''
180 Plot for EEJ Spectral Width
181 '''
182
183 CODE = 'width'
184 colormap = 'RdBu_r'
185 colormap = 'jet'
186
187 def setup(self):
188
189 self.xaxis = 'time'
190 self.ncols = 1
191 self.nrows = 2
192 self.nplots = 2
193 self.ylabel = 'Range [km]'
194 self.xlabel = 'Time'
195 self.cb_label = '(m/s)'
196 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
197 self.titles = ['{} EJJ Type {} /'.format(
198 self.CODE.upper(), x) for x in range(1,1+self.nrows)]
199
200 def update(self, dataOut):
201
202 if dataOut.mode == 11: #Double Gaussian
203 width = numpy.append(dataOut.Oblique_params[:,2,:],dataOut.Oblique_params[:,5,:],axis=0)
204 elif dataOut.mode == 9: #Double Skew Gaussian
205 width = numpy.append(dataOut.Oblique_params[:,2,:],dataOut.Oblique_params[:,6,:],axis=0)
206 data = {
207 'width': width
95 }
208 }
96
209
97 return data, {}
210 return data, {}
@@ -107,7 +220,7 class PowerPlot(RTIPlot):
107 def update(self, dataOut):
220 def update(self, dataOut):
108
221
109 data = {
222 data = {
110 'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
223 'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
111 }
224 }
112
225
113 return data, {}
226 return data, {}
@@ -208,7 +321,7 class GenericRTIPlot(Plot):
208 meta = {}
321 meta = {}
209
322
210 return data, meta
323 return data, meta
211
324
212 def plot(self):
325 def plot(self):
213 # self.data.normalize_heights()
326 # self.data.normalize_heights()
214 self.x = self.data.times
327 self.x = self.data.times
Show file before | Show file after | Hide comments
@@ -39,9 +39,14 class SpectraPlot(Plot):
39
39
40 data = {}
40 data = {}
41 meta = {}
41 meta = {}
42
42 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
43 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
44 #print("Spc: ",spc[0])
45 #exit(1)
43 data['spc'] = spc
46 data['spc'] = spc
44 data['rti'] = dataOut.getPower()
47 data['rti'] = dataOut.getPower()
48 #print(data['rti'][0])
49 #exit(1)
45 #print("NormFactor: ",dataOut.normFactor)
50 #print("NormFactor: ",dataOut.normFactor)
46 #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
51 #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
47 if hasattr(dataOut, 'LagPlot'): #Double Pulse
52 if hasattr(dataOut, 'LagPlot'): #Double Pulse
@@ -163,11 +168,22 class SpectraObliquePlot(Plot):
163 data['rti'] = dataOut.getPower()
168 data['rti'] = dataOut.getPower()
164 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
169 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
165 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
170 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
166
171 '''
167 data['shift1'] = dataOut.Oblique_params[0][1]
172 data['shift1'] = dataOut.Oblique_params[0,-2,:]
168 data['shift2'] = dataOut.Oblique_params[0][4]
173 data['shift2'] = dataOut.Oblique_params[0,-1,:]
169 data['shift1_error'] = dataOut.Oblique_param_errors[0][1]
174 data['shift1_error'] = dataOut.Oblique_param_errors[0,-2,:]
170 data['shift2_error'] = dataOut.Oblique_param_errors[0][4]
175 data['shift2_error'] = dataOut.Oblique_param_errors[0,-1,:]
176 '''
177 '''
178 data['shift1'] = dataOut.Oblique_params[0,1,:]
179 data['shift2'] = dataOut.Oblique_params[0,4,:]
180 data['shift1_error'] = dataOut.Oblique_param_errors[0,1,:]
181 data['shift2_error'] = dataOut.Oblique_param_errors[0,4,:]
182 '''
183 data['shift1'] = dataOut.Dop_EEJ_T1[0]
184 data['shift2'] = dataOut.Dop_EEJ_T2[0]
185 data['shift1_error'] = dataOut.Err_Dop_EEJ_T1[0]
186 data['shift2_error'] = dataOut.Err_Dop_EEJ_T2[0]
171
187
172 return data, meta
188 return data, meta
173
189
@@ -187,15 +203,19 class SpectraObliquePlot(Plot):
187
203
188 y = self.data.yrange
204 y = self.data.yrange
189 self.y = y
205 self.y = y
190 z = self.data['spc']
206
207 data = self.data[-1]
208 z = data['spc']
191
209
192 for n, ax in enumerate(self.axes):
210 for n, ax in enumerate(self.axes):
193 noise = self.data['noise'][n][-1]
211 noise = self.data['noise'][n][-1]
194 shift1 = self.data['shift1']
212 shift1 = data['shift1']
195 shift2 = self.data['shift2']
213 #print(shift1)
196 err1 = self.data['shift1_error']
214 shift2 = data['shift2']
197 err2 = self.data['shift2_error']
215 err1 = data['shift1_error']
216 err2 = data['shift2_error']
198 if ax.firsttime:
217 if ax.firsttime:
218
199 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
219 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
200 self.xmin = self.xmin if self.xmin else -self.xmax
220 self.xmin = self.xmin if self.xmin else -self.xmax
201 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
221 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
@@ -212,17 +232,20 class SpectraObliquePlot(Plot):
212 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
232 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
213 color="k", linestyle="dashed", lw=1)[0]
233 color="k", linestyle="dashed", lw=1)[0]
214
234
215 self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^', elinewidth=0.2, marker='x', linestyle='None',markersize=0.5,capsize=0.3,markeredgewidth=0.2)
235 self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^', elinewidth=2.2, marker='o', linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
216 self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=0.2,marker='x',linestyle='None',markersize=0.5,capsize=0.3,markeredgewidth=0.2)
236 self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
237 #print("plotter1: ", self.ploterr1,shift1)
238
217 else:
239 else:
240 #print("else plotter1: ", self.ploterr1,shift1)
218 self.ploterr1.remove()
241 self.ploterr1.remove()
219 self.ploterr2.remove()
242 self.ploterr2.remove()
220 ax.plt.set_array(z[n].T.ravel())
243 ax.plt.set_array(z[n].T.ravel())
221 if self.showprofile:
244 if self.showprofile:
222 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
245 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
223 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
246 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
224 self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^',elinewidth=0.2,marker='x',linestyle='None',markersize=0.5,capsize=0.3,markeredgewidth=0.2)
247 self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^', elinewidth=2.2, marker='o', linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
225 self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=0.2,marker='x',linestyle='None',markersize=0.5,capsize=0.3,markeredgewidth=0.2)
248 self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
226
249
227 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
250 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
228
251
@@ -671,6 +694,7 class RTIPlot(Plot):
671 data = {}
694 data = {}
672 meta = {}
695 meta = {}
673 data['rti'] = dataOut.getPower()
696 data['rti'] = dataOut.getPower()
697 #print(numpy.shape(data['rti']))
674
698
675 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
699 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
676
700
@@ -692,6 +716,7 class RTIPlot(Plot):
692 for n, ax in enumerate(self.axes):
716 for n, ax in enumerate(self.axes):
693 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
717 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
694 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
718 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
719
695 if ax.firsttime:
720 if ax.firsttime:
696 ax.plt = ax.pcolormesh(x, y, z[n].T,
721 ax.plt = ax.pcolormesh(x, y, z[n].T,
697 vmin=self.zmin,
722 vmin=self.zmin,
Show file before | Show file after | Hide comments
@@ -213,7 +213,7 class DenRTIPlot(RTIPlot):
213 data = {}
213 data = {}
214 meta = {}
214 meta = {}
215
215
216 data['denrti'] = dataOut.DensityFinal
216 data['denrti'] = dataOut.DensityFinal*1.e-6 #To Plot in cm^-3
217
217
218 return data, meta
218 return data, meta
219
219
Show file before | Show file after | Hide comments
@@ -70,7 +70,7 class ProcessingUnit(object):
70 def call(self, **kwargs):
70 def call(self, **kwargs):
71 '''
71 '''
72 '''
72 '''
73 #print("call")
73
74 try:
74 try:
75 if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
75 if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
76 #if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error and not self.dataIn.runNextUnit:
76 #if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error and not self.dataIn.runNextUnit:
@@ -125,6 +125,7 class ProcessingUnit(object):
125 runNextUnit = self.dataOut.isReady()
125 runNextUnit = self.dataOut.isReady()
126 except:
126 except:
127 runNextUnit = self.dataOut.isReady()
127 runNextUnit = self.dataOut.isReady()
128 #exit(1)
128 #if not self.dataOut.isReady():
129 #if not self.dataOut.isReady():
129 #return 'Error' if self.dataOut.error else input()
130 #return 'Error' if self.dataOut.error else input()
130 #print("NexT",runNextUnit)
131 #print("NexT",runNextUnit)
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (1211 lines changed) Show them Hide them
@@ -2,6 +2,7 import numpy
2 import math
2 import math
3 from scipy import optimize, interpolate, signal, stats, ndimage
3 from scipy import optimize, interpolate, signal, stats, ndimage
4 import scipy
4 import scipy
5 from scipy.optimize import least_squares
5 import re
6 import re
6 import datetime
7 import datetime
7 import copy
8 import copy
@@ -91,6 +92,7 class ParametersProc(ProcessingUnit):
91 self.dataOut.timeInterval1 = self.dataIn.timeInterval
92 self.dataOut.timeInterval1 = self.dataIn.timeInterval
92 self.dataOut.heightList = self.dataIn.heightList
93 self.dataOut.heightList = self.dataIn.heightList
93 self.dataOut.frequency = self.dataIn.frequency
94 self.dataOut.frequency = self.dataIn.frequency
95 #self.dataOut.runNextUnit = self.dataIn.runNextUnit
94 #self.dataOut.noise = self.dataIn.noise
96 #self.dataOut.noise = self.dataIn.noise
95
97
96 def run(self):
98 def run(self):
@@ -156,6 +158,7 class ParametersProc(ProcessingUnit):
156 if hasattr(self.dataIn, 'COFA'): #COFA
158 if hasattr(self.dataIn, 'COFA'): #COFA
157 self.dataOut.COFA = self.dataIn.COFA
159 self.dataOut.COFA = self.dataIn.COFA
158
160
161 #self.dataOut.runNextUnit = self.dataIn.runNextUnit
159
162
160
163
161 #---------------------- Correlation Data ---------------------------
164 #---------------------- Correlation Data ---------------------------
@@ -668,7 +671,9 class GaussianFit(Operation):
668 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.)
671 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.)
669
672
670 class Oblique_Gauss_Fit(Operation):
673 class Oblique_Gauss_Fit(Operation):
671
674 '''
675 Written by R. Flores
676 '''
672 def __init__(self):
677 def __init__(self):
673 Operation.__init__(self)
678 Operation.__init__(self)
674
679
@@ -1060,6 +1065,17 class Oblique_Gauss_Fit(Operation):
1060 val = a1 * numpy.exp(-z1**2/2) + a2 * numpy.exp(-y2**2/2)/(1-k2*z2) + d
1065 val = a1 * numpy.exp(-z1**2/2) + a2 * numpy.exp(-y2**2/2)/(1-k2*z2) + d
1061 return val
1066 return val
1062
1067
1068 def gaussian(self, x, a, b, c, d):
1069 z = (x-b)/c
1070 val = a * numpy.exp(-z**2/2) + d
1071 return val
1072
1073 def double_gaussian(self, x, a1, b1, c1, a2, b2, c2, d):
1074 z1 = (x-b1)/c1
1075 z2 = (x-b2)/c2
1076 val = a1 * numpy.exp(-z1**2/2) + a2 * numpy.exp(-z2**2/2) + d
1077 return val
1078
1063 def double_gaussian_double_skew(self,x, a1, b1, c1, k1, a2, b2, c2, k2, d):
1079 def double_gaussian_double_skew(self,x, a1, b1, c1, k1, a2, b2, c2, k2, d):
1064
1080
1065 z1 = (x-b1)/c1
1081 z1 = (x-b1)/c1
@@ -1136,7 +1152,7 class Oblique_Gauss_Fit(Operation):
1136 #return A1f, B1f, C1f, A2f, B2f, C2f, K2f, Df, doppler
1152 #return A1f, B1f, C1f, A2f, B2f, C2f, K2f, Df, doppler
1137 return A1f, B1f, C1f, A2f, B2f, C2f, K2f, Df, doppler
1153 return A1f, B1f, C1f, A2f, B2f, C2f, K2f, Df, doppler
1138
1154
1139 def Double_Gauss_Double_Skew_fit_weight_bound_no_inputs(self,spc,freq):
1155 def Double_Gauss_Double_Skew_fit_weight_bound_no_inputs(self,spc,freq,Nincoh):
1140
1156
1141 from scipy.optimize import least_squares
1157 from scipy.optimize import least_squares
1142
1158
@@ -1146,6 +1162,7 class Oblique_Gauss_Fit(Operation):
1146 from scipy.signal import medfilt
1162 from scipy.signal import medfilt
1147 Nincoh = 20
1163 Nincoh = 20
1148 Nincoh = 80
1164 Nincoh = 80
1165 Nincoh = Nincoh
1149 spcm = medfilt(spc,11)/numpy.sqrt(Nincoh)
1166 spcm = medfilt(spc,11)/numpy.sqrt(Nincoh)
1150
1167
1151 # define a least squares function to optimize
1168 # define a least squares function to optimize
@@ -1156,11 +1173,23 class Oblique_Gauss_Fit(Operation):
1156 # bounds=([0,-460,0,0,-400,120,0],[numpy.inf,-340,50,numpy.inf,0,250,numpy.inf])
1173 # bounds=([0,-460,0,0,-400,120,0],[numpy.inf,-340,50,numpy.inf,0,250,numpy.inf])
1157 # bounds=([0,-numpy.inf,0,0,-numpy.inf,0,-numpy.inf,0],[numpy.inf,-200,numpy.inf,numpy.inf,0,numpy.inf,0,numpy.inf])
1174 # bounds=([0,-numpy.inf,0,0,-numpy.inf,0,-numpy.inf,0],[numpy.inf,-200,numpy.inf,numpy.inf,0,numpy.inf,0,numpy.inf])
1158 #print(a1,b1,c1,a2,b2,c2,k2,d)
1175 #print(a1,b1,c1,a2,b2,c2,k2,d)
1159 bounds=([0,-numpy.inf,0,-numpy.inf,0,-400,0,0,0],[numpy.inf,-340,numpy.inf,0,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf])
1176 #bounds=([0,-numpy.inf,0,-numpy.inf,0,-400,0,0,0],[numpy.inf,-340,numpy.inf,0,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf])
1177 bounds=([0,-numpy.inf,0,-numpy.inf,0,-400,0,0,0],[numpy.inf,-140,numpy.inf,0,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf])
1160 #print(bounds)
1178 #print(bounds)
1161 #bounds=([0,-numpy.inf,0,0,-numpy.inf,0,0,0],[numpy.inf,-200,numpy.inf,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf])
1179 #bounds=([0,-numpy.inf,0,0,-numpy.inf,0,0,0],[numpy.inf,-200,numpy.inf,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf])
1162 params_scale = [spc_max,freq_max,freq_max,1,spc_max,freq_max,freq_max,1,spc_max]
1180 params_scale = [spc_max,freq_max,freq_max,1,spc_max,freq_max,freq_max,1,spc_max]
1163 x0_value = numpy.array([spc_max,-400,30,-.1,spc_max/4,-200,150,1,1.0e7])
1181 ####################x0_value = numpy.array([spc_max,-400,30,-.1,spc_max/4,-200,150,1,1.0e7])
1182
1183 dop1_x0 = freq[numpy.argmax(spc)]
1184 ####dop1_x0 = freq[numpy.argmax(spcm)]
1185 if dop1_x0 < 0:
1186 dop2_x0 = dop1_x0 + 100
1187 if dop1_x0 > 0:
1188 dop2_x0 = dop1_x0 - 100
1189
1190 ###########x0_value = numpy.array([spc_max,-200.5,30,-.1,spc_max/4,-100.5,150,1,1.0e7])
1191 x0_value = numpy.array([spc_max,dop1_x0,30,-.1,spc_max/4, dop2_x0,150,1,1.0e7])
1192 #x0_value = numpy.array([spc_max,-400.5,30,-.1,spc_max/4,-200.5,150,1,1.0e7])
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)
1193 #popt = least_squares(lsq_func,[a1,b1,c1,a2,b2,c2,k2,d],x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=1)
1165 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1194 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1166 # popt = least_squares(lsq_func,[a1,b1,c1,a2,b2,c2,k2,d],x_scale=params_scale,verbose=1)
1195 # popt = least_squares(lsq_func,[a1,b1,c1,a2,b2,c2,k2,d],x_scale=params_scale,verbose=1)
@@ -1193,6 +1222,65 class Oblique_Gauss_Fit(Operation):
1193 #exit(1)
1222 #exit(1)
1194 return A1f, B1f, C1f, K1f, A2f, B2f, C2f, K2f, Df, doppler1, doppler2, error
1223 return A1f, B1f, C1f, K1f, A2f, B2f, C2f, K2f, Df, doppler1, doppler2, error
1195
1224
1225 def Double_Gauss_fit_weight_bound_no_inputs(self,spc,freq,Nincoh):
1226
1227 from scipy.optimize import least_squares
1228
1229 freq_max = numpy.max(numpy.abs(freq))
1230 spc_max = numpy.max(spc)
1231
1232 from scipy.signal import medfilt
1233 Nincoh = 20
1234 Nincoh = 80
1235 Nincoh = Nincoh
1236 spcm = medfilt(spc,11)/numpy.sqrt(Nincoh)
1237
1238 # define a least squares function to optimize
1239 def lsq_func(params):
1240 return (spc-self.double_gaussian(freq,params[0],params[1],params[2],params[3],params[4],params[5],params[6]))/spcm
1241
1242 # fit
1243 # bounds=([0,-460,0,0,-400,120,0],[numpy.inf,-340,50,numpy.inf,0,250,numpy.inf])
1244 # bounds=([0,-numpy.inf,0,0,-numpy.inf,0,-numpy.inf,0],[numpy.inf,-200,numpy.inf,numpy.inf,0,numpy.inf,0,numpy.inf])
1245 #print(a1,b1,c1,a2,b2,c2,k2,d)
1246
1247 dop1_x0 = freq[numpy.argmax(spcm)]
1248
1249 #####bounds=([0,-numpy.inf,0,0,-400,0,0],[numpy.inf,-340,numpy.inf,numpy.inf,0,numpy.inf,numpy.inf])
1250 #####bounds=([0,-numpy.inf,0,0,dop1_x0-50,0,0],[numpy.inf,-340,numpy.inf,numpy.inf,0,numpy.inf,numpy.inf])
1251 bounds=([0,-numpy.inf,0,0,dop1_x0-50,0,0],[numpy.inf,-300,numpy.inf,numpy.inf,0,numpy.inf,numpy.inf])
1252 #####bounds=([0,-numpy.inf,0,0,-500,0,0],[numpy.inf,-340,numpy.inf,numpy.inf,0,numpy.inf,numpy.inf])
1253 #bounds=([0,-numpy.inf,0,-numpy.inf,0,-500,0,0,0],[numpy.inf,-240,numpy.inf,0,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf])
1254 #print(bounds)
1255 #bounds=([0,-numpy.inf,0,0,-numpy.inf,0,0,0],[numpy.inf,-200,numpy.inf,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf])
1256 params_scale = [spc_max,freq_max,freq_max,spc_max,freq_max,freq_max,spc_max]
1257 #x0_value = numpy.array([spc_max,-400.5,30,spc_max/4,-200.5,150,1.0e7])
1258 x0_value = numpy.array([spc_max,-400.5,30,spc_max/4,dop1_x0,150,1.0e7])
1259 #x0_value = numpy.array([spc_max,-420.5,30,-.1,spc_max/4,-50,150,.1,numpy.mean(spc[-50:])])
1260 #print("before popt")
1261 #print(x0_value)
1262 #print("freq: ",freq)
1263 #popt = least_squares(lsq_func,[a1,b1,c1,a2,b2,c2,k2,d],x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=1)
1264 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1265 # popt = least_squares(lsq_func,[a1,b1,c1,a2,b2,c2,k2,d],x_scale=params_scale,verbose=1)
1266 #print("after popt")
1267 J = popt.jac
1268
1269 try:
1270 cov = numpy.linalg.inv(J.T.dot(J))
1271 error = numpy.sqrt(numpy.diagonal(cov))
1272 except:
1273 error = numpy.ones((7))*numpy.NAN
1274
1275 A1f = popt.x[0]; B1f = popt.x[1]; C1f = popt.x[2]
1276 A2f = popt.x[3]; B2f = popt.x[4]; C2f = popt.x[5]
1277 Df = popt.x[6]
1278 #print("before return")
1279 return A1f, B1f, C1f, A2f, B2f, C2f, Df, error
1280
1281
1282
1283
1196 def Double_Gauss_Double_Skew_fit_weight_bound_with_inputs(self, spc, freq, a1, b1, c1, a2, b2, c2, k2, d):
1284 def Double_Gauss_Double_Skew_fit_weight_bound_with_inputs(self, spc, freq, a1, b1, c1, a2, b2, c2, k2, d):
1197
1285
1198 from scipy.optimize import least_squares
1286 from scipy.optimize import least_squares
@@ -1264,7 +1352,124 class Oblique_Gauss_Fit(Operation):
1264
1352
1265 return A1f, B1f, C1f, A2f, B2f, C2f, K2f, A3f, B3f, C3f, K3f, Df, doppler
1353 return A1f, B1f, C1f, A2f, B2f, C2f, K2f, A3f, B3f, C3f, K3f, Df, doppler
1266
1354
1267 def run(self, dataOut, mode = 0):
1355 def CEEJ_Skew_fit_weight_bound_no_inputs(self,spc,freq,Nincoh):
1356
1357 from scipy.optimize import least_squares
1358
1359 freq_max = numpy.max(numpy.abs(freq))
1360 spc_max = numpy.max(spc)
1361
1362 from scipy.signal import medfilt
1363 Nincoh = 20
1364 Nincoh = 80
1365 Nincoh = Nincoh
1366 spcm = medfilt(spc,11)/numpy.sqrt(Nincoh)
1367
1368 # define a least squares function to optimize
1369 def lsq_func(params):
1370 return (spc-self.gaussian_skew(freq,params[0],params[1],params[2],params[3],params[4]))#/spcm
1371
1372
1373 bounds=([0,0,0,-numpy.inf,0],[numpy.inf,numpy.inf,numpy.inf,0,numpy.inf])
1374
1375 params_scale = [spc_max,freq_max,freq_max,1,spc_max]
1376
1377 x0_value = numpy.array([spc_max,freq[numpy.argmax(spc)],30,-.1,numpy.mean(spc[:50])])
1378
1379 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1380
1381 J = popt.jac
1382
1383 try:
1384 error = numpy.ones((9))*numpy.NAN
1385 cov = numpy.linalg.inv(J.T.dot(J))
1386 error[:4] = numpy.sqrt(numpy.diagonal(cov))[:4]
1387 error[-1] = numpy.sqrt(numpy.diagonal(cov))[-1]
1388 except:
1389 error = numpy.ones((9))*numpy.NAN
1390
1391 A1f = popt.x[0]; B1f = popt.x[1]; C1f = popt.x[2]; K1f = popt.x[3]
1392 Df = popt.x[4]
1393
1394 aux1 = self.gaussian_skew(freq, A1f, B1f, C1f, K1f, Df)
1395 doppler1 = freq[numpy.argmax(aux1)]
1396 #print("CEEJ ERROR:",error)
1397
1398 return A1f, B1f, C1f, K1f, numpy.NAN, numpy.NAN, numpy.NAN, numpy.NAN, Df, doppler1, numpy.NAN, error
1399
1400 def CEEJ_fit_weight_bound_no_inputs(self,spc,freq,Nincoh):
1401
1402 from scipy.optimize import least_squares
1403
1404 freq_max = numpy.max(numpy.abs(freq))
1405 spc_max = numpy.max(spc)
1406
1407 from scipy.signal import medfilt
1408 Nincoh = 20
1409 Nincoh = 80
1410 Nincoh = Nincoh
1411 spcm = medfilt(spc,11)/numpy.sqrt(Nincoh)
1412
1413 # define a least squares function to optimize
1414 def lsq_func(params):
1415 return (spc-self.gaussian(freq,params[0],params[1],params[2],params[3]))#/spcm
1416
1417
1418 bounds=([0,0,0,0],[numpy.inf,numpy.inf,numpy.inf,numpy.inf])
1419
1420 params_scale = [spc_max,freq_max,freq_max,spc_max]
1421
1422 x0_value = numpy.array([spc_max,freq[numpy.argmax(spcm)],30,numpy.mean(spc[:50])])
1423
1424 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1425
1426 J = popt.jac
1427
1428 try:
1429 error = numpy.ones((4))*numpy.NAN
1430 cov = numpy.linalg.inv(J.T.dot(J))
1431 error = numpy.sqrt(numpy.diagonal(cov))
1432 except:
1433 error = numpy.ones((4))*numpy.NAN
1434
1435 A1f = popt.x[0]; B1f = popt.x[1]; C1f = popt.x[2]
1436 Df = popt.x[3]
1437
1438 return A1f, B1f, C1f, Df, error
1439
1440 def Simple_fit_bound(self,spc,freq,Nincoh):
1441
1442 freq_max = numpy.max(numpy.abs(freq))
1443 spc_max = numpy.max(spc)
1444
1445 Nincoh = Nincoh
1446
1447 def lsq_func(params):
1448 return (spc-self.gaussian(freq,params[0],params[1],params[2],params[3]))
1449
1450 bounds=([0,-50,0,0],[numpy.inf,+50,numpy.inf,numpy.inf])
1451
1452 params_scale = [spc_max,freq_max,freq_max,spc_max]
1453
1454 x0_value = numpy.array([spc_max,-20.5,5,1.0e7])
1455
1456 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1457
1458 J = popt.jac
1459
1460 try:
1461 cov = numpy.linalg.inv(J.T.dot(J))
1462 error = numpy.sqrt(numpy.diagonal(cov))
1463 except:
1464 error = numpy.ones((4))*numpy.NAN
1465
1466 A1f = popt.x[0]; B1f = popt.x[1]; C1f = popt.x[2]
1467 Df = popt.x[3]
1468
1469 return A1f, B1f, C1f, Df, error
1470
1471
1472 def run(self, dataOut, mode = 0, Hmin1 = None, Hmax1 = None, Hmin2 = None, Hmax2 = None):
1268
1473
1269 pwcode = 1
1474 pwcode = 1
1270
1475
@@ -1293,12 +1498,54 class Oblique_Gauss_Fit(Operation):
1293 elif mode == 9:
1498 elif mode == 9:
1294 dataOut.Oblique_params = numpy.ones((1,11,dataOut.nHeights))*numpy.NAN
1499 dataOut.Oblique_params = numpy.ones((1,11,dataOut.nHeights))*numpy.NAN
1295 dataOut.Oblique_param_errors = numpy.ones((1,9,dataOut.nHeights))*numpy.NAN
1500 dataOut.Oblique_param_errors = numpy.ones((1,9,dataOut.nHeights))*numpy.NAN
1501 elif mode == 11:
1502 dataOut.Oblique_params = numpy.ones((1,7,dataOut.nHeights))*numpy.NAN
1503 dataOut.Oblique_param_errors = numpy.ones((1,7,dataOut.nHeights))*numpy.NAN
1504 elif mode == 10: #150 km
1505 dataOut.Oblique_params = numpy.ones((1,4,dataOut.nHeights))*numpy.NAN
1506 dataOut.Oblique_param_errors = numpy.ones((1,4,dataOut.nHeights))*numpy.NAN
1507 dataOut.snr_log10 = numpy.ones((1,dataOut.nHeights))*numpy.NAN
1296
1508
1297 dataOut.VelRange = x
1509 dataOut.VelRange = x
1298
1510
1299 l1=range(22,36)
1511
1512
1513 #l1=range(22,36) #+62
1300 #l1=range(32,36)
1514 #l1=range(32,36)
1301 l2=range(58,99)
1515 #l2=range(58,99) #+62
1516
1517 #if Hmin1 == None or Hmax1 == None or Hmin2 == None or Hmax2 == None:
1518
1519 minHei1 = 105.
1520 maxHei1 = 122.5
1521 maxHei1 = 130.5
1522
1523 if mode == 10: #150 km
1524 minHei1 = 100
1525 maxHei1 = 100
1526
1527 inda1 = numpy.where(dataOut.heightList >= minHei1)
1528 indb1 = numpy.where(dataOut.heightList <= maxHei1)
1529
1530 minIndex1 = inda1[0][0]
1531 maxIndex1 = indb1[0][-1]
1532
1533 minHei2 = 150.
1534 maxHei2 = 201.25
1535 maxHei2 = 225.3
1536
1537 if mode == 10: #150 km
1538 minHei2 = 110
1539 maxHei2 = 165
1540
1541 inda2 = numpy.where(dataOut.heightList >= minHei2)
1542 indb2 = numpy.where(dataOut.heightList <= maxHei2)
1543
1544 minIndex2 = inda2[0][0]
1545 maxIndex2 = indb2[0][-1]
1546
1547 l1=range(minIndex1,maxIndex1)
1548 l2=range(minIndex2,maxIndex2)
1302
1549
1303 if mode == 4:
1550 if mode == 4:
1304 '''
1551 '''
@@ -1321,9 +1568,10 class Oblique_Gauss_Fit(Operation):
1321
1568
1322 for hei in itertools.chain(l1, l2):
1569 for hei in itertools.chain(l1, l2):
1323 #for hei in range(79,81):
1570 #for hei in range(79,81):
1571 if numpy.isnan(dataOut.data_snr[0,hei]):
1572 continue #Avoids the analysis when there is only noise
1324
1573
1325 try:
1574 try:
1326
1327 spc = dataOut.data_spc[0,:,hei]
1575 spc = dataOut.data_spc[0,:,hei]
1328
1576
1329 if mode == 6: #Skew Weighted Bounded
1577 if mode == 6: #Skew Weighted Bounded
@@ -1340,14 +1588,39 class Oblique_Gauss_Fit(Operation):
1340 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,9,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei]))
1588 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,9,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei]))
1341
1589
1342 elif mode == 9: #Double Skewed Weighted Bounded no inputs
1590 elif mode == 9: #Double Skewed Weighted Bounded no inputs
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)
1591 #if numpy.max(spc) <= 0:
1344 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,10,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei]))
1592 if x[numpy.argmax(spc)] <= 0:
1345 #print(hei)
1593 #print("EEJ")
1346 #print(dataOut.Oblique_params[0,10,hei])
1594 dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_params[0,4,hei],dataOut.Oblique_params[0,5,hei],dataOut.Oblique_params[0,6,hei],dataOut.Oblique_params[0,7,hei],dataOut.Oblique_params[0,8,hei],dataOut.Oblique_params[0,9,hei],dataOut.Oblique_params[0,10,hei],dataOut.Oblique_param_errors[0,:,hei] = self.Double_Gauss_Double_Skew_fit_weight_bound_no_inputs(spc,x,dataOut.nIncohInt)
1347 #print(dataOut.dplr_2_u[0,0,hei])
1595 #if dataOut.Oblique_params[0,-2,hei] < -500 or dataOut.Oblique_params[0,-2,hei] > 500 or dataOut.Oblique_params[0,-1,hei] < -500 or dataOut.Oblique_params[0,-1,hei] > 500:
1348 #print("outside",dataOut.Oblique_param_errors[0,:,hei])
1596 # dataOut.Oblique_params[0,:,hei] *= numpy.NAN
1349 #print("SUCCESSSSSSS")
1597 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,10,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei]))
1350 #exit(1)
1598
1599 else:
1600 #print("CEEJ")
1601 dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_params[0,4,hei],dataOut.Oblique_params[0,5,hei],dataOut.Oblique_params[0,6,hei],dataOut.Oblique_params[0,7,hei],dataOut.Oblique_params[0,8,hei],dataOut.Oblique_params[0,9,hei],dataOut.Oblique_params[0,10,hei],dataOut.Oblique_param_errors[0,:,hei] = self.CEEJ_Skew_fit_weight_bound_no_inputs(spc,x,dataOut.nIncohInt)
1602 #if dataOut.Oblique_params[0,-2,hei] < -500 or dataOut.Oblique_params[0,-2,hei] > 500 or dataOut.Oblique_params[0,-1,hei] < -500 or dataOut.Oblique_params[0,-1,hei] > 500:
1603 # dataOut.Oblique_params[0,:,hei] *= numpy.NAN
1604 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,10,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei]))
1605 elif mode == 11: #Double Weighted Bounded no inputs
1606 #if numpy.max(spc) <= 0:
1607 from scipy.signal import medfilt
1608 spcm = medfilt(spc,11)
1609
1610 if x[numpy.argmax(spcm)] <= 0:
1611 #print("EEJ")
1612 #print("EEJ",dataOut.heightList[hei])
1613 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_param_errors[0,:,hei] = self.Double_Gauss_fit_weight_bound_no_inputs(spc,x,dataOut.nIncohInt)
1614 #if dataOut.Oblique_params[0,-2,hei] < -500 or dataOut.Oblique_params[0,-2,hei] > 500 or dataOut.Oblique_params[0,-1,hei] < -500 or dataOut.Oblique_params[0,-1,hei] > 500:
1615 # dataOut.Oblique_params[0,:,hei] *= numpy.NAN
1616 else:
1617 #print("CEEJ",dataOut.heightList[hei])
1618 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_param_errors[0,:,hei] = self.CEEJ_fit_weight_bound_no_inputs(spc,x,dataOut.nIncohInt)
1619
1620 elif mode == 10: #150km
1621 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_param_errors[0,:,hei] = self.Simple_fit_bound(spc,x,dataOut.nIncohInt)
1622 snr = (dataOut.power[0,hei]*factor - dataOut.Oblique_params[0,3,hei])/dataOut.Oblique_params[0,3,hei]
1623 dataOut.snr_log10[0,hei] = numpy.log10(snr)
1351
1624
1352 else:
1625 else:
1353 spc_fit, A1, B1, C1, D1 = self.Gauss_fit_2(spc,x,'first')
1626 spc_fit, A1, B1, C1, D1 = self.Gauss_fit_2(spc,x,'first')
@@ -1393,7 +1666,34 class Oblique_Gauss_Fit(Operation):
1393 ###dataOut.Oblique_params[0,:,hei] = dataOut.Oblique_params[0,:,hei]*numpy.NAN
1666 ###dataOut.Oblique_params[0,:,hei] = dataOut.Oblique_params[0,:,hei]*numpy.NAN
1394 pass
1667 pass
1395
1668
1396 #exit(1)
1669 #exit(1)
1670 dataOut.paramInterval = dataOut.nProfiles*dataOut.nCohInt*dataOut.ippSeconds
1671 dataOut.lat=-11.95
1672 dataOut.lon=-76.87
1673
1674 if mode == 9: #Double Skew Gaussian
1675 dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,-2,:]
1676 dataOut.Spec_W_T1 = dataOut.Oblique_params[:,2,:]
1677 dataOut.Dop_EEJ_T2 = dataOut.Oblique_params[:,-1,:]
1678 dataOut.Spec_W_T2 = dataOut.Oblique_params[:,6,:]
1679
1680 dataOut.Err_Dop_EEJ_T1 = dataOut.Oblique_param_errors[:,1,:] #En realidad este es el error?
1681 dataOut.Err_Spec_W_T1 = dataOut.Oblique_param_errors[:,2,:]
1682 dataOut.Err_Dop_EEJ_T2 = dataOut.Oblique_param_errors[:,5,:] #En realidad este es el error?
1683 dataOut.Err_Spec_W_T2 = dataOut.Oblique_param_errors[:,6,:]
1684
1685 elif mode == 11: #Double Gaussian
1686 dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,1,:]
1687 dataOut.Spec_W_T1 = dataOut.Oblique_params[:,2,:]
1688 dataOut.Dop_EEJ_T2 = dataOut.Oblique_params[:,4,:]
1689 dataOut.Spec_W_T2 = dataOut.Oblique_params[:,5,:]
1690
1691 dataOut.Err_Dop_EEJ_T1 = dataOut.Oblique_param_errors[:,1,:]
1692 dataOut.Err_Spec_W_T1 = dataOut.Oblique_param_errors[:,2,:]
1693 dataOut.Err_Dop_EEJ_T2 = dataOut.Oblique_param_errors[:,4,:]
1694 dataOut.Err_Spec_W_T2 = dataOut.Oblique_param_errors[:,5,:]
1695
1696 dataOut.mode = mode
1397
1697
1398 return dataOut
1698 return dataOut
1399
1699
@@ -3104,7 +3404,7 class SpectralFitting(Operation):
3104 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
3404 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
3105 return chisq
3405 return chisq
3106
3406
3107 class WindProfiler(Operation):
3407 class WindProfiler_V0(Operation):
3108
3408
3109 __isConfig = False
3409 __isConfig = False
3110
3410
@@ -3656,7 +3956,8 class WindProfiler(Operation):
3656 def run(self, dataOut, technique, nHours=1, hmin=70, hmax=110, **kwargs):
3956 def run(self, dataOut, technique, nHours=1, hmin=70, hmax=110, **kwargs):
3657
3957
3658 param = dataOut.data_param
3958 param = dataOut.data_param
3659 if dataOut.abscissaList != None:
3959 #if dataOut.abscissaList != None:
3960 if numpy.any(dataOut.abscissaList):
3660 absc = dataOut.abscissaList[:-1]
3961 absc = dataOut.abscissaList[:-1]
3661 # noise = dataOut.noise
3962 # noise = dataOut.noise
3662 heightList = dataOut.heightList
3963 heightList = dataOut.heightList
@@ -3821,11 +4122,60 class WindProfiler(Operation):
3821
4122
3822 return
4123 return
3823
4124
3824 class EWDriftsEstimation(Operation):
4125 class WindProfiler(Operation):
4126
4127 __isConfig = False
4128
4129 __initime = None
4130 __lastdatatime = None
4131 __integrationtime = None
4132
4133 __buffer = None
4134
4135 __dataReady = False
4136
4137 __firstdata = None
4138
4139 n = None
3825
4140
3826 def __init__(self):
4141 def __init__(self):
3827 Operation.__init__(self)
4142 Operation.__init__(self)
3828
4143
4144 def __calculateCosDir(self, elev, azim):
4145 zen = (90 - elev)*numpy.pi/180
4146 azim = azim*numpy.pi/180
4147 cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
4148 cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
4149
4150 signX = numpy.sign(numpy.cos(azim))
4151 signY = numpy.sign(numpy.sin(azim))
4152
4153 cosDirX = numpy.copysign(cosDirX, signX)
4154 cosDirY = numpy.copysign(cosDirY, signY)
4155 return cosDirX, cosDirY
4156
4157 def __calculateAngles(self, theta_x, theta_y, azimuth):
4158
4159 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
4160 zenith_arr = numpy.arccos(dir_cosw)
4161 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
4162
4163 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
4164 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
4165
4166 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
4167
4168 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
4169
4170 if horOnly:
4171 A = numpy.c_[dir_cosu,dir_cosv]
4172 else:
4173 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
4174 A = numpy.asmatrix(A)
4175 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
4176
4177 return A1
4178
3829 def __correctValues(self, heiRang, phi, velRadial, SNR):
4179 def __correctValues(self, heiRang, phi, velRadial, SNR):
3830 listPhi = phi.tolist()
4180 listPhi = phi.tolist()
3831 maxid = listPhi.index(max(listPhi))
4181 maxid = listPhi.index(max(listPhi))
@@ -3845,23 +4195,13 class EWDriftsEstimation(Operation):
3845 for i in rango:
4195 for i in rango:
3846 x = heiRang*math.cos(phi[i])
4196 x = heiRang*math.cos(phi[i])
3847 y1 = velRadial[i,:]
4197 y1 = velRadial[i,:]
3848 vali= (numpy.isfinite(y1)==True).nonzero()
4198 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
3849 y1=y1[vali]
3850 x = x[vali]
3851 f1 = interpolate.interp1d(x,y1,kind = 'cubic',bounds_error=False)
3852
4199
3853 #heiRang1 = x*math.cos(phi[maxid])
3854 x1 = heiRang1
4200 x1 = heiRang1
3855 y11 = f1(x1)
4201 y11 = f1(x1)
3856
4202
3857 y2 = SNR[i,:]
4203 y2 = SNR[i,:]
3858 #print 'snr ', y2
4204 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
3859 x = heiRang*math.cos(phi[i])
3860 vali= (y2 != -1).nonzero()
3861 y2 = y2[vali]
3862 x = x[vali]
3863 #print 'snr ',y2
3864 f2 = interpolate.interp1d(x,y2,kind = 'cubic',bounds_error=False)
3865 y21 = f2(x1)
4205 y21 = f2(x1)
3866
4206
3867 velRadial1[i,:] = y11
4207 velRadial1[i,:] = y11
@@ -3869,36 +4209,716 class EWDriftsEstimation(Operation):
3869
4209
3870 return heiRang1, velRadial1, SNR1
4210 return heiRang1, velRadial1, SNR1
3871
4211
4212 def __calculateVelUVW(self, A, velRadial):
3872
4213
4214 #Operacion Matricial
4215 # velUVW = numpy.zeros((velRadial.shape[1],3))
4216 # for ind in range(velRadial.shape[1]):
4217 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
4218 # velUVW = velUVW.transpose()
4219 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
4220 velUVW[:,:] = numpy.dot(A,velRadial)
3873
4221
3874 def run(self, dataOut, zenith, zenithCorrection):
3875
4222
3876 heiRang = dataOut.heightList
4223 return velUVW
3877 velRadial = dataOut.data_param[:,3,:]
3878 velRadialm = dataOut.data_param[:,2:4,:]*-1
3879
4224
3880 rbufc=dataOut.data_paramC[:,:,0]
4225 # def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
3881 ebufc=dataOut.data_paramC[:,:,1]
3882 SNR = dataOut.data_snr
3883 velRerr = dataOut.data_error[:,4,:]
3884 moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]],[velRadialm[1,:]],[velRadialm[1,:]]))
3885 dataOut.moments=moments
3886 # Coherent
3887 smooth_wC = ebufc[0,:]
3888 p_w0C = rbufc[0,:]
3889 p_w1C = rbufc[1,:]
3890 w_wC = rbufc[2,:]*-1 #*radial_sign(radial EQ 1)
3891 t_wC = rbufc[3,:]
3892 my_nbeams = 2
3893
4226
3894 zenith = numpy.array(zenith)
4227 def techniqueDBS(self, kwargs):
3895 zenith -= zenithCorrection
4228 """
3896 zenith *= numpy.pi/180
4229 Function that implements Doppler Beam Swinging (DBS) technique.
3897 if zenithCorrection != 0 :
4230
3898 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
4231 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
3899 else :
4232 Direction correction (if necessary), Ranges and SNR
3900 heiRang1 = heiRang
4233
3901 velRadial1 = velRadial
4234 Output: Winds estimation (Zonal, Meridional and Vertical)
4235
4236 Parameters affected: Winds, height range, SNR
4237 """
4238 velRadial0 = kwargs['velRadial']
4239 heiRang = kwargs['heightList']
4240 SNR0 = kwargs['SNR']
4241
4242 if 'dirCosx' in kwargs and 'dirCosy' in kwargs:
4243 theta_x = numpy.array(kwargs['dirCosx'])
4244 theta_y = numpy.array(kwargs['dirCosy'])
4245 else:
4246 elev = numpy.array(kwargs['elevation'])
4247 azim = numpy.array(kwargs['azimuth'])
4248 theta_x, theta_y = self.__calculateCosDir(elev, azim)
4249 azimuth = kwargs['correctAzimuth']
4250 if 'horizontalOnly' in kwargs:
4251 horizontalOnly = kwargs['horizontalOnly']
4252 else: horizontalOnly = False
4253 if 'correctFactor' in kwargs:
4254 correctFactor = kwargs['correctFactor']
4255 else: correctFactor = 1
4256 if 'channelList' in kwargs:
4257 channelList = kwargs['channelList']
4258 if len(channelList) == 2:
4259 horizontalOnly = True
4260 arrayChannel = numpy.array(channelList)
4261 param = param[arrayChannel,:,:]
4262 theta_x = theta_x[arrayChannel]
4263 theta_y = theta_y[arrayChannel]
4264
4265 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
4266 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0)
4267 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
4268
4269 #Calculo de Componentes de la velocidad con DBS
4270 winds = self.__calculateVelUVW(A,velRadial1)
4271
4272 return winds, heiRang1, SNR1
4273
4274 def __calculateDistance(self, posx, posy, pairs_ccf, azimuth = None):
4275
4276 nPairs = len(pairs_ccf)
4277 posx = numpy.asarray(posx)
4278 posy = numpy.asarray(posy)
4279
4280 #Rotacion Inversa para alinear con el azimuth
4281 if azimuth!= None:
4282 azimuth = azimuth*math.pi/180
4283 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
4284 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
4285 else:
4286 posx1 = posx
4287 posy1 = posy
4288
4289 #Calculo de Distancias
4290 distx = numpy.zeros(nPairs)
4291 disty = numpy.zeros(nPairs)
4292 dist = numpy.zeros(nPairs)
4293 ang = numpy.zeros(nPairs)
4294
4295 for i in range(nPairs):
4296 distx[i] = posx1[pairs_ccf[i][1]] - posx1[pairs_ccf[i][0]]
4297 disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]]
4298 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
4299 ang[i] = numpy.arctan2(disty[i],distx[i])
4300
4301 return distx, disty, dist, ang
4302 #Calculo de Matrices
4303 # nPairs = len(pairs)
4304 # ang1 = numpy.zeros((nPairs, 2, 1))
4305 # dist1 = numpy.zeros((nPairs, 2, 1))
4306 #
4307 # for j in range(nPairs):
4308 # dist1[j,0,0] = dist[pairs[j][0]]
4309 # dist1[j,1,0] = dist[pairs[j][1]]
4310 # ang1[j,0,0] = ang[pairs[j][0]]
4311 # ang1[j,1,0] = ang[pairs[j][1]]
4312 #
4313 # return distx,disty, dist1,ang1
4314
4315
4316 def __calculateVelVer(self, phase, lagTRange, _lambda):
4317
4318 Ts = lagTRange[1] - lagTRange[0]
4319 velW = -_lambda*phase/(4*math.pi*Ts)
4320
4321 return velW
4322
4323 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
4324 nPairs = tau1.shape[0]
4325 nHeights = tau1.shape[1]
4326 vel = numpy.zeros((nPairs,3,nHeights))
4327 dist1 = numpy.reshape(dist, (dist.size,1))
4328
4329 angCos = numpy.cos(ang)
4330 angSin = numpy.sin(ang)
4331
4332 vel0 = dist1*tau1/(2*tau2**2)
4333 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
4334 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
4335
4336 ind = numpy.where(numpy.isinf(vel))
4337 vel[ind] = numpy.nan
4338
4339 return vel
4340
4341 # def __getPairsAutoCorr(self, pairsList, nChannels):
4342 #
4343 # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
4344 #
4345 # for l in range(len(pairsList)):
4346 # firstChannel = pairsList[l][0]
4347 # secondChannel = pairsList[l][1]
4348 #
4349 # #Obteniendo pares de Autocorrelacion
4350 # if firstChannel == secondChannel:
4351 # pairsAutoCorr[firstChannel] = int(l)
4352 #
4353 # pairsAutoCorr = pairsAutoCorr.astype(int)
4354 #
4355 # pairsCrossCorr = range(len(pairsList))
4356 # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
4357 #
4358 # return pairsAutoCorr, pairsCrossCorr
4359
4360 # def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
4361 def techniqueSA(self, kwargs):
4362
4363 """
4364 Function that implements Spaced Antenna (SA) technique.
4365
4366 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
4367 Direction correction (if necessary), Ranges and SNR
4368
4369 Output: Winds estimation (Zonal, Meridional and Vertical)
4370
4371 Parameters affected: Winds
4372 """
4373 position_x = kwargs['positionX']
4374 position_y = kwargs['positionY']
4375 azimuth = kwargs['azimuth']
4376
4377 if 'correctFactor' in kwargs:
4378 correctFactor = kwargs['correctFactor']
4379 else:
4380 correctFactor = 1
4381
4382 groupList = kwargs['groupList']
4383 pairs_ccf = groupList[1]
4384 tau = kwargs['tau']
4385 _lambda = kwargs['_lambda']
4386
4387 #Cross Correlation pairs obtained
4388 # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairssList, nChannels)
4389 # pairsArray = numpy.array(pairsList)[pairsCrossCorr]
4390 # pairsSelArray = numpy.array(pairsSelected)
4391 # pairs = []
4392 #
4393 # #Wind estimation pairs obtained
4394 # for i in range(pairsSelArray.shape[0]/2):
4395 # ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
4396 # ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
4397 # pairs.append((ind1,ind2))
4398
4399 indtau = tau.shape[0]/2
4400 tau1 = tau[:indtau,:]
4401 tau2 = tau[indtau:-1,:]
4402 # tau1 = tau1[pairs,:]
4403 # tau2 = tau2[pairs,:]
4404 phase1 = tau[-1,:]
4405
4406 #---------------------------------------------------------------------
4407 #Metodo Directo
4408 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairs_ccf,azimuth)
4409 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
4410 winds = stats.nanmean(winds, axis=0)
4411 #---------------------------------------------------------------------
4412 #Metodo General
4413 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
4414 # #Calculo Coeficientes de Funcion de Correlacion
4415 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
4416 # #Calculo de Velocidades
4417 # winds = self.calculateVelUV(F,G,A,B,H)
4418
4419 #---------------------------------------------------------------------
4420 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
4421 winds = correctFactor*winds
4422 return winds
4423
4424 def __checkTime(self, currentTime, paramInterval, outputInterval):
4425
4426 dataTime = currentTime + paramInterval
4427 deltaTime = dataTime - self.__initime
4428
4429 if deltaTime >= outputInterval or deltaTime < 0:
4430 self.__dataReady = True
4431 return
4432
4433 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
4434 '''
4435 Function that implements winds estimation technique with detected meteors.
4436
4437 Input: Detected meteors, Minimum meteor quantity to wind estimation
4438
4439 Output: Winds estimation (Zonal and Meridional)
4440
4441 Parameters affected: Winds
4442 '''
4443 #Settings
4444 nInt = (heightMax - heightMin)/2
4445 nInt = int(nInt)
4446 winds = numpy.zeros((2,nInt))*numpy.nan
4447
4448 #Filter errors
4449 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
4450 finalMeteor = arrayMeteor[error,:]
4451
4452 #Meteor Histogram
4453 finalHeights = finalMeteor[:,2]
4454 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
4455 nMeteorsPerI = hist[0]
4456 heightPerI = hist[1]
4457
4458 #Sort of meteors
4459 indSort = finalHeights.argsort()
4460 finalMeteor2 = finalMeteor[indSort,:]
4461
4462 # Calculating winds
4463 ind1 = 0
4464 ind2 = 0
4465
4466 for i in range(nInt):
4467 nMet = nMeteorsPerI[i]
4468 ind1 = ind2
4469 ind2 = ind1 + nMet
4470
4471 meteorAux = finalMeteor2[ind1:ind2,:]
4472
4473 if meteorAux.shape[0] >= meteorThresh:
4474 vel = meteorAux[:, 6]
4475 zen = meteorAux[:, 4]*numpy.pi/180
4476 azim = meteorAux[:, 3]*numpy.pi/180
4477
4478 n = numpy.cos(zen)
4479 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
4480 # l = m*numpy.tan(azim)
4481 l = numpy.sin(zen)*numpy.sin(azim)
4482 m = numpy.sin(zen)*numpy.cos(azim)
4483
4484 A = numpy.vstack((l, m)).transpose()
4485 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
4486 windsAux = numpy.dot(A1, vel)
4487
4488 winds[0,i] = windsAux[0]
4489 winds[1,i] = windsAux[1]
4490
4491 return winds, heightPerI[:-1]
4492
4493 def techniqueNSM_SA(self, **kwargs):
4494 metArray = kwargs['metArray']
4495 heightList = kwargs['heightList']
4496 timeList = kwargs['timeList']
4497
4498 rx_location = kwargs['rx_location']
4499 groupList = kwargs['groupList']
4500 azimuth = kwargs['azimuth']
4501 dfactor = kwargs['dfactor']
4502 k = kwargs['k']
4503
4504 azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth)
4505 d = dist*dfactor
4506 #Phase calculation
4507 metArray1 = self.__getPhaseSlope(metArray, heightList, timeList)
4508
4509 metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities
4510
4511 velEst = numpy.zeros((heightList.size,2))*numpy.nan
4512 azimuth1 = azimuth1*numpy.pi/180
4513
4514 for i in range(heightList.size):
4515 h = heightList[i]
4516 indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0]
4517 metHeight = metArray1[indH,:]
4518 if metHeight.shape[0] >= 2:
4519 velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities
4520 iazim = metHeight[:,1].astype(int)
4521 azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths
4522 A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux)))
4523 A = numpy.asmatrix(A)
4524 A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose()
4525 velHor = numpy.dot(A1,velAux)
4526
4527 velEst[i,:] = numpy.squeeze(velHor)
4528 return velEst
4529
4530 def __getPhaseSlope(self, metArray, heightList, timeList):
4531 meteorList = []
4532 #utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2
4533 #Putting back together the meteor matrix
4534 utctime = metArray[:,0]
4535 uniqueTime = numpy.unique(utctime)
4536
4537 phaseDerThresh = 0.5
4538 ippSeconds = timeList[1] - timeList[0]
4539 sec = numpy.where(timeList>1)[0][0]
4540 nPairs = metArray.shape[1] - 6
4541 nHeights = len(heightList)
4542
4543 for t in uniqueTime:
4544 metArray1 = metArray[utctime==t,:]
4545 # phaseDerThresh = numpy.pi/4 #reducir Phase thresh
4546 tmet = metArray1[:,1].astype(int)
4547 hmet = metArray1[:,2].astype(int)
4548
4549 metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1))
4550 metPhase[:,:] = numpy.nan
4551 metPhase[:,hmet,tmet] = metArray1[:,6:].T
4552
4553 #Delete short trails
4554 metBool = ~numpy.isnan(metPhase[0,:,:])
4555 heightVect = numpy.sum(metBool, axis = 1)
4556 metBool[heightVect<sec,:] = False
4557 metPhase[:,heightVect<sec,:] = numpy.nan
4558
4559 #Derivative
4560 metDer = numpy.abs(metPhase[:,:,1:] - metPhase[:,:,:-1])
4561 phDerAux = numpy.dstack((numpy.full((nPairs,nHeights,1), False, dtype=bool),metDer > phaseDerThresh))
4562 metPhase[phDerAux] = numpy.nan
4563
4564 #--------------------------METEOR DETECTION -----------------------------------------
4565 indMet = numpy.where(numpy.any(metBool,axis=1))[0]
4566
4567 for p in numpy.arange(nPairs):
4568 phase = metPhase[p,:,:]
4569 phDer = metDer[p,:,:]
4570
4571 for h in indMet:
4572 height = heightList[h]
4573 phase1 = phase[h,:] #82
4574 phDer1 = phDer[h,:]
4575
4576 phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap
4577
4578 indValid = numpy.where(~numpy.isnan(phase1))[0]
4579 initMet = indValid[0]
4580 endMet = 0
4581
4582 for i in range(len(indValid)-1):
4583
4584 #Time difference
4585 inow = indValid[i]
4586 inext = indValid[i+1]
4587 idiff = inext - inow
4588 #Phase difference
4589 phDiff = numpy.abs(phase1[inext] - phase1[inow])
4590
4591 if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor
4592 sizeTrail = inow - initMet + 1
4593 if sizeTrail>3*sec: #Too short meteors
4594 x = numpy.arange(initMet,inow+1)*ippSeconds
4595 y = phase1[initMet:inow+1]
4596 ynnan = ~numpy.isnan(y)
4597 x = x[ynnan]
4598 y = y[ynnan]
4599 slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
4600 ylin = x*slope + intercept
4601 rsq = r_value**2
4602 if rsq > 0.5:
4603 vel = slope#*height*1000/(k*d)
4604 estAux = numpy.array([utctime,p,height, vel, rsq])
4605 meteorList.append(estAux)
4606 initMet = inext
4607 metArray2 = numpy.array(meteorList)
4608
4609 return metArray2
4610
4611 def __calculateAzimuth1(self, rx_location, pairslist, azimuth0):
4612
4613 azimuth1 = numpy.zeros(len(pairslist))
4614 dist = numpy.zeros(len(pairslist))
4615
4616 for i in range(len(rx_location)):
4617 ch0 = pairslist[i][0]
4618 ch1 = pairslist[i][1]
4619
4620 diffX = rx_location[ch0][0] - rx_location[ch1][0]
4621 diffY = rx_location[ch0][1] - rx_location[ch1][1]
4622 azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi
4623 dist[i] = numpy.sqrt(diffX**2 + diffY**2)
4624
4625 azimuth1 -= azimuth0
4626 return azimuth1, dist
4627
4628 def techniqueNSM_DBS(self, **kwargs):
4629 metArray = kwargs['metArray']
4630 heightList = kwargs['heightList']
4631 timeList = kwargs['timeList']
4632 azimuth = kwargs['azimuth']
4633 theta_x = numpy.array(kwargs['theta_x'])
4634 theta_y = numpy.array(kwargs['theta_y'])
4635
4636 utctime = metArray[:,0]
4637 cmet = metArray[:,1].astype(int)
4638 hmet = metArray[:,3].astype(int)
4639 SNRmet = metArray[:,4]
4640 vmet = metArray[:,5]
4641 spcmet = metArray[:,6]
4642
4643 nChan = numpy.max(cmet) + 1
4644 nHeights = len(heightList)
4645
4646 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
4647 hmet = heightList[hmet]
4648 h1met = hmet*numpy.cos(zenith_arr[cmet]) #Corrected heights
4649
4650 velEst = numpy.zeros((heightList.size,2))*numpy.nan
4651
4652 for i in range(nHeights - 1):
4653 hmin = heightList[i]
4654 hmax = heightList[i + 1]
4655
4656 thisH = (h1met>=hmin) & (h1met<hmax) & (cmet!=2) & (SNRmet>8) & (vmet<50) & (spcmet<10)
4657 indthisH = numpy.where(thisH)
4658
4659 if numpy.size(indthisH) > 3:
4660
4661 vel_aux = vmet[thisH]
4662 chan_aux = cmet[thisH]
4663 cosu_aux = dir_cosu[chan_aux]
4664 cosv_aux = dir_cosv[chan_aux]
4665 cosw_aux = dir_cosw[chan_aux]
4666
4667 nch = numpy.size(numpy.unique(chan_aux))
4668 if nch > 1:
4669 A = self.__calculateMatA(cosu_aux, cosv_aux, cosw_aux, True)
4670 velEst[i,:] = numpy.dot(A,vel_aux)
4671
4672 return velEst
4673
4674 def run(self, dataOut, technique, nHours=1, hmin=70, hmax=110, **kwargs):
4675
4676 param = dataOut.moments
4677 #param = dataOut.data_param
4678 #if dataOut.abscissaList != None:
4679 if numpy.any(dataOut.abscissaList) :
4680 absc = dataOut.abscissaList[:-1]
4681 # noise = dataOut.noise
4682 heightList = dataOut.heightList
4683 SNR = dataOut.data_snr
4684
4685 if technique == 'DBS':
4686
4687 kwargs['velRadial'] = param[:,1,:] #Radial velocity
4688 kwargs['heightList'] = heightList
4689 kwargs['SNR'] = SNR
4690
4691 dataOut.data_output, dataOut.heightList, dataOut.data_snr = self.techniqueDBS(kwargs) #DBS Function
4692 dataOut.utctimeInit = dataOut.utctime
4693 dataOut.outputInterval = dataOut.paramInterval
4694
4695 elif technique == 'SA':
4696
4697 #Parameters
4698 # position_x = kwargs['positionX']
4699 # position_y = kwargs['positionY']
4700 # azimuth = kwargs['azimuth']
4701 #
4702 # if kwargs.has_key('crosspairsList'):
4703 # pairs = kwargs['crosspairsList']
4704 # else:
4705 # pairs = None
4706 #
4707 # if kwargs.has_key('correctFactor'):
4708 # correctFactor = kwargs['correctFactor']
4709 # else:
4710 # correctFactor = 1
4711
4712 # tau = dataOut.data_param
4713 # _lambda = dataOut.C/dataOut.frequency
4714 # pairsList = dataOut.groupList
4715 # nChannels = dataOut.nChannels
4716
4717 kwargs['groupList'] = dataOut.groupList
4718 kwargs['tau'] = dataOut.data_param
4719 kwargs['_lambda'] = dataOut.C/dataOut.frequency
4720 # dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
4721 dataOut.data_output = self.techniqueSA(kwargs)
4722 dataOut.utctimeInit = dataOut.utctime
4723 dataOut.outputInterval = dataOut.timeInterval
4724
4725 elif technique == 'Meteors':
4726 dataOut.flagNoData = True
4727 self.__dataReady = False
4728
4729 if 'nHours' in kwargs:
4730 nHours = kwargs['nHours']
4731 else:
4732 nHours = 1
4733
4734 if 'meteorsPerBin' in kwargs:
4735 meteorThresh = kwargs['meteorsPerBin']
4736 else:
4737 meteorThresh = 6
4738
4739 if 'hmin' in kwargs:
4740 hmin = kwargs['hmin']
4741 else: hmin = 70
4742 if 'hmax' in kwargs:
4743 hmax = kwargs['hmax']
4744 else: hmax = 110
4745
4746 dataOut.outputInterval = nHours*3600
4747
4748 if self.__isConfig == False:
4749 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
4750 #Get Initial LTC time
4751 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
4752 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
4753
4754 self.__isConfig = True
4755
4756 if self.__buffer is None:
4757 self.__buffer = dataOut.data_param
4758 self.__firstdata = copy.copy(dataOut)
4759
4760 else:
4761 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
4762
4763 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
4764
4765 if self.__dataReady:
4766 dataOut.utctimeInit = self.__initime
4767
4768 self.__initime += dataOut.outputInterval #to erase time offset
4769
4770 dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
4771 dataOut.flagNoData = False
4772 self.__buffer = None
4773
4774 elif technique == 'Meteors1':
4775 dataOut.flagNoData = True
4776 self.__dataReady = False
4777
4778 if 'nMins' in kwargs:
4779 nMins = kwargs['nMins']
4780 else: nMins = 20
4781 if 'rx_location' in kwargs:
4782 rx_location = kwargs['rx_location']
4783 else: rx_location = [(0,1),(1,1),(1,0)]
4784 if 'azimuth' in kwargs:
4785 azimuth = kwargs['azimuth']
4786 else: azimuth = 51.06
4787 if 'dfactor' in kwargs:
4788 dfactor = kwargs['dfactor']
4789 if 'mode' in kwargs:
4790 mode = kwargs['mode']
4791 if 'theta_x' in kwargs:
4792 theta_x = kwargs['theta_x']
4793 if 'theta_y' in kwargs:
4794 theta_y = kwargs['theta_y']
4795 else: mode = 'SA'
4796
4797 #Borrar luego esto
4798 if dataOut.groupList is None:
4799 dataOut.groupList = [(0,1),(0,2),(1,2)]
4800 groupList = dataOut.groupList
4801 C = 3e8
4802 freq = 50e6
4803 lamb = C/freq
4804 k = 2*numpy.pi/lamb
4805
4806 timeList = dataOut.abscissaList
4807 heightList = dataOut.heightList
4808
4809 if self.__isConfig == False:
4810 dataOut.outputInterval = nMins*60
4811 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
4812 #Get Initial LTC time
4813 initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
4814 minuteAux = initime.minute
4815 minuteNew = int(numpy.floor(minuteAux/nMins)*nMins)
4816 self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
4817
4818 self.__isConfig = True
4819
4820 if self.__buffer is None:
4821 self.__buffer = dataOut.data_param
4822 self.__firstdata = copy.copy(dataOut)
4823
4824 else:
4825 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
4826
4827 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
4828
4829 if self.__dataReady:
4830 dataOut.utctimeInit = self.__initime
4831 self.__initime += dataOut.outputInterval #to erase time offset
4832
4833 metArray = self.__buffer
4834 if mode == 'SA':
4835 dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)
4836 elif mode == 'DBS':
4837 dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList, azimuth=azimuth, theta_x=theta_x, theta_y=theta_y)
4838 dataOut.data_output = dataOut.data_output.T
4839 dataOut.flagNoData = False
4840 self.__buffer = None
4841 #print("ENDDD")
4842 return dataOut
4843
4844 class EWDriftsEstimation(Operation):
4845
4846 def __init__(self):
4847 Operation.__init__(self)
4848
4849 def __correctValues(self, heiRang, phi, velRadial, SNR):
4850 listPhi = phi.tolist()
4851 maxid = listPhi.index(max(listPhi))
4852 minid = listPhi.index(min(listPhi))
4853
4854 rango = list(range(len(phi)))
4855 # rango = numpy.delete(rango,maxid)
4856
4857 heiRang1 = heiRang*math.cos(phi[maxid])
4858 heiRangAux = heiRang*math.cos(phi[minid])
4859 indOut = (heiRang1 < heiRangAux[0]).nonzero()
4860 heiRang1 = numpy.delete(heiRang1,indOut)
4861
4862 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
4863 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
4864
4865 for i in rango:
4866 x = heiRang*math.cos(phi[i])
4867 y1 = velRadial[i,:]
4868 vali= (numpy.isfinite(y1)==True).nonzero()
4869 y1=y1[vali]
4870 x = x[vali]
4871 f1 = interpolate.interp1d(x,y1,kind = 'cubic',bounds_error=False)
4872
4873 #heiRang1 = x*math.cos(phi[maxid])
4874 x1 = heiRang1
4875 y11 = f1(x1)
4876
4877 y2 = SNR[i,:]
4878 #print 'snr ', y2
4879 x = heiRang*math.cos(phi[i])
4880 vali= (y2 != -1).nonzero()
4881 y2 = y2[vali]
4882 x = x[vali]
4883 #print 'snr ',y2
4884 f2 = interpolate.interp1d(x,y2,kind = 'cubic',bounds_error=False)
4885 y21 = f2(x1)
4886
4887 velRadial1[i,:] = y11
4888 SNR1[i,:] = y21
4889
4890 return heiRang1, velRadial1, SNR1
4891
4892
4893
4894 def run(self, dataOut, zenith, zenithCorrection):
4895
4896 heiRang = dataOut.heightList
4897 velRadial = dataOut.data_param[:,3,:]
4898 velRadialm = dataOut.data_param[:,2:4,:]*-1
4899
4900 rbufc=dataOut.data_paramC[:,:,0]
4901 ebufc=dataOut.data_paramC[:,:,1]
4902 SNR = dataOut.data_snr
4903 velRerr = dataOut.data_error[:,4,:]
4904 moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]],[velRadialm[1,:]],[velRadialm[1,:]]))
4905 dataOut.moments=moments
4906 # Coherent
4907 smooth_wC = ebufc[0,:]
4908 p_w0C = rbufc[0,:]
4909 p_w1C = rbufc[1,:]
4910 w_wC = rbufc[2,:]*-1 #*radial_sign(radial EQ 1)
4911 t_wC = rbufc[3,:]
4912 my_nbeams = 2
4913
4914 zenith = numpy.array(zenith)
4915 zenith -= zenithCorrection
4916 zenith *= numpy.pi/180
4917 if zenithCorrection != 0 :
4918 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
4919 else :
4920 heiRang1 = heiRang
4921 velRadial1 = velRadial
3902 SNR1 = SNR
4922 SNR1 = SNR
3903
4923
3904 alp = zenith[0]
4924 alp = zenith[0]
@@ -5542,6 +6562,9 class SMOperations():
5542
6562
5543
6563
5544 class IGRFModel(Operation):
6564 class IGRFModel(Operation):
6565 '''
6566 Written by R. Flores
6567 '''
5545 """Operation to calculate Geomagnetic parameters.
6568 """Operation to calculate Geomagnetic parameters.
5546
6569
5547 Parameters:
6570 Parameters:
@@ -5597,10 +6620,13 class MergeProc(ProcessingUnit):
5597 ProcessingUnit.__init__(self)
6620 ProcessingUnit.__init__(self)
5598
6621
5599 def run(self, attr_data, attr_data_2 = None, attr_data_3 = None, attr_data_4 = None, attr_data_5 = None, mode=0):
6622 def run(self, attr_data, attr_data_2 = None, attr_data_3 = None, attr_data_4 = None, attr_data_5 = None, mode=0):
6623 #print("*****************************Merge***************")
5600
6624
5601 self.dataOut = getattr(self, self.inputs[0])
6625 self.dataOut = getattr(self, self.inputs[0])
5602 data_inputs = [getattr(self, attr) for attr in self.inputs]
6626 data_inputs = [getattr(self, attr) for attr in self.inputs]
5603 #print(data_inputs)
6627 #print(data_inputs)
6628 #print("Run: ",self.dataOut.runNextUnit)
6629 #exit(1)
5604 #print(numpy.shape([getattr(data, attr_data) for data in data_inputs][1]))
6630 #print(numpy.shape([getattr(data, attr_data) for data in data_inputs][1]))
5605 #exit(1)
6631 #exit(1)
5606 if mode==0:
6632 if mode==0:
@@ -5668,39 +6694,7 class MergeProc(ProcessingUnit):
5668 #setattr(self.dataOut, 'data_acf', getattr(data_inputs[1], attr_data_5)) #LP
6694 #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
6695 setattr(self.dataOut, 'data_acf', getattr(data_inputs[1], attr_data_5)) #LP
5670 #print("Merge data_acf: ",self.dataOut.data_acf.shape)
6696 #print("Merge data_acf: ",self.dataOut.data_acf.shape)
5671 #exit(1)
6697
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
6698
5705 #self.dataOut.nIncohInt_LP = 128
6699 #self.dataOut.nIncohInt_LP = 128
5706 #self.dataOut.nProfiles_LP = 128#self.dataOut.nIncohInt_LP
6700 #self.dataOut.nProfiles_LP = 128#self.dataOut.nIncohInt_LP
@@ -5711,8 +6705,45 class MergeProc(ProcessingUnit):
5711 #print("sahpi",self.dataOut.nIncohInt_LP)
6705 #print("sahpi",self.dataOut.nIncohInt_LP)
5712 #exit(1)
6706 #exit(1)
5713 self.dataOut.NLAG = 16
6707 self.dataOut.NLAG = 16
6708 self.dataOut.NLAG = self.dataOut.data_acf.shape[1]
5714 self.dataOut.NRANGE = self.dataOut.data_acf.shape[-1]
6709 self.dataOut.NRANGE = self.dataOut.data_acf.shape[-1]
5715
6710
5716 #print(numpy.shape(self.dataOut.data_spc))
6711 #print(numpy.shape(self.dataOut.data_spc))
5717
6712
5718 #exit(1)
6713 #exit(1)
6714 if mode==5:
6715 data = numpy.concatenate([getattr(data, attr_data) for data in data_inputs])
6716 setattr(self.dataOut, attr_data, data)
6717 data = numpy.concatenate([getattr(data, attr_data_2) for data in data_inputs])
6718 setattr(self.dataOut, attr_data_2, data)
6719
6720 if mode==6: #Hybrid Spectra-Voltage
6721 #data = numpy.concatenate([getattr(data, attr_data) for data in data_inputs],axis=1)
6722 #setattr(self.dataOut, attr_data, data)
6723 setattr(self.dataOut, 'dataLag_spc', getattr(data_inputs[1], attr_data)) #DP
6724 setattr(self.dataOut, 'dataLag_cspc', getattr(data_inputs[1], attr_data_2)) #DP
6725 setattr(self.dataOut, 'output_LP_integrated', getattr(data_inputs[0], attr_data_3)) #LP
6726 #setattr(self.dataOut, 'dataLag_cspc_LP', getattr(data_inputs[1], attr_data_4)) #LP
6727 #setattr(self.dataOut, 'data_acf', getattr(data_inputs[1], attr_data_5)) #LP
6728 #setattr(self.dataOut, 'data_acf', getattr(data_inputs[1], attr_data_5)) #LP
6729 #print("Merge data_acf: ",self.dataOut.data_acf.shape)
6730 #print(self.dataOut.NSCAN)
6731 self.dataOut.nIncohInt = int(self.dataOut.NAVG * self.dataOut.nint)
6732 #print(self.dataOut.dataLag_spc.shape)
6733 self.dataOut.nProfiles = self.dataOut.nProfiles_DP = self.dataOut.dataLag_spc.shape[1]
6734 '''
6735 #self.dataOut.nIncohInt_LP = 128
6736 #self.dataOut.nProfiles_LP = 128#self.dataOut.nIncohInt_LP
6737 self.dataOut.nProfiles_LP = 16#28#self.dataOut.nIncohInt_LP
6738 self.dataOut.nProfiles_LP = self.dataOut.data_acf.shape[1]#28#self.dataOut.nIncohInt_LP
6739 self.dataOut.NSCAN = 128
6740 self.dataOut.nIncohInt_LP = self.dataOut.nIncohInt*self.dataOut.NSCAN
6741 #print("sahpi",self.dataOut.nIncohInt_LP)
6742 #exit(1)
6743 self.dataOut.NLAG = 16
6744 self.dataOut.NLAG = self.dataOut.data_acf.shape[1]
6745 self.dataOut.NRANGE = self.dataOut.data_acf.shape[-1]
6746 '''
6747 #print(numpy.shape(self.dataOut.data_spc))
6748 print("*************************GOOD*************************")
6749 #exit(1)
Show file before | Show file after | Hide comments
@@ -123,9 +123,10 class SpectraProc(ProcessingUnit):
123 self.dataOut.flagShiftFFT = False
123 self.dataOut.flagShiftFFT = False
124
124
125 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False, runNextUnit = 0):
125 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False, runNextUnit = 0):
126
126
127 self.dataIn.runNextUnit = runNextUnit
127 self.dataIn.runNextUnit = runNextUnit
128 if self.dataIn.type == "Spectra":
128 if self.dataIn.type == "Spectra":
129
129 self.dataOut.copy(self.dataIn)
130 self.dataOut.copy(self.dataIn)
130 if shift_fft:
131 if shift_fft:
131 #desplaza a la derecha en el eje 2 determinadas posiciones
132 #desplaza a la derecha en el eje 2 determinadas posiciones
@@ -147,9 +148,15 class SpectraProc(ProcessingUnit):
147
148
148 if nProfiles == None:
149 if nProfiles == None:
149 nProfiles = nFFTPoints
150 nProfiles = nFFTPoints
150
151 #print(self.dataOut.ipp)
152 #exit(1)
151 if ippFactor == None:
153 if ippFactor == None:
152 self.dataOut.ippFactor = 1
154 self.dataOut.ippFactor = 1
155 #if ippFactor is not None:
156 #self.dataOut.ippFactor = ippFactor
157 #print(ippFactor)
158 #print(self.dataOut.ippFactor)
159 #exit(1)
153
160
154 self.dataOut.nFFTPoints = nFFTPoints
161 self.dataOut.nFFTPoints = nFFTPoints
155
162
@@ -425,6 +432,46 class SpectraProc(ProcessingUnit):
425
432
426 return 1
433 return 1
427
434
435 class GetSNR(Operation):
436 '''
437 Written by R. Flores
438 '''
439 """Operation to get SNR.
440
441 Parameters:
442 -----------
443
444 Example
445 --------
446
447 op = proc_unit.addOperation(name='GetSNR', optype='other')
448
449 """
450
451 def __init__(self, **kwargs):
452
453 Operation.__init__(self, **kwargs)
454
455
456 def run(self,dataOut):
457
458 noise = dataOut.getNoise()
459 maxdB = 16
460
461 normFactor = 24
462
463 #dataOut.data_snr = (dataOut.data_spc.sum(axis=1))/(noise[:,None]*dataOut.normFactor)
464 dataOut.data_snr = (dataOut.data_spc.sum(axis=1))/(noise[:,None]*dataOut.nFFTPoints)
465
466 dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.5, numpy.nan, dataOut.data_snr)
467 #dataOut.data_snr = 10*numpy.log10(dataOut.data_snr)
468 #dataOut.data_snr = numpy.expand_dims(dataOut.data_snr,axis=0)
469 #print(dataOut.data_snr.shape)
470 #exit(1)
471
472
473 return dataOut
474
428 class removeDC(Operation):
475 class removeDC(Operation):
429
476
430 def run(self, dataOut, mode=2):
477 def run(self, dataOut, mode=2):
Show file before | Show file after | Hide comments
@@ -19,10 +19,15 from schainpy.model.data.jrodata import hildebrand_sekhon
19 from schainpy.utils import log
19 from schainpy.utils import log
20 from schainpy.model.data import _HS_algorithm
20 from schainpy.model.data import _HS_algorithm
21
21
22 from schainpy.model.proc.jroproc_voltage import CleanCohEchoes
23
22 from time import time, mktime, strptime, gmtime, ctime
24 from time import time, mktime, strptime, gmtime, ctime
23
25
24
26
25 class SpectraLagProc(ProcessingUnit):
27 class SpectraLagProc(ProcessingUnit):
28 '''
29 Written by R. Flores
30 '''
26 def __init__(self):
31 def __init__(self):
27
32
28 ProcessingUnit.__init__(self)
33 ProcessingUnit.__init__(self)
@@ -308,7 +313,9 class SpectraLagProc(ProcessingUnit):
308 #print("after",self.dataOut.data_spc[0,:,20])
313 #print("after",self.dataOut.data_spc[0,:,20])
309
314
310 class removeDCLag(Operation):
315 class removeDCLag(Operation):
311
316 '''
317 Written by R. Flores
318 '''
312 def remover(self,mode):
319 def remover(self,mode):
313 jspectra = self.dataOut.data_spc
320 jspectra = self.dataOut.data_spc
314 jcspectra = self.dataOut.data_cspc
321 jcspectra = self.dataOut.data_cspc
@@ -414,7 +421,9 class removeDCLag(Operation):
414
421
415
422
416 class removeDCLagFlip(Operation):
423 class removeDCLagFlip(Operation):
417
424 '''
425 Written by R. Flores
426 '''
418 #CHANGES MADE ONLY FOR MODE 2 AND NOT CONSIDERING CSPC
427 #CHANGES MADE ONLY FOR MODE 2 AND NOT CONSIDERING CSPC
419
428
420 def remover(self,mode):
429 def remover(self,mode):
@@ -501,7 +510,7 class removeDCLagFlip(Operation):
501
510
502
511
503 def run(self, dataOut, mode=2):
512 def run(self, dataOut, mode=2):
504 #print("***********************************Remove DC***********************************")
513 print("***********************************Remove DC***********************************")
505 ##print(dataOut.FlipChannels)
514 ##print(dataOut.FlipChannels)
506 #exit(1)
515 #exit(1)
507 self.dataOut = dataOut
516 self.dataOut = dataOut
@@ -1098,7 +1107,9 class removeInterference(Operation):
1098
1107
1099
1108
1100 class IntegrationFaradaySpectra(Operation):
1109 class IntegrationFaradaySpectra(Operation):
1101
1110 '''
1111 Written by R. Flores
1112 '''
1102 __profIndex = 0
1113 __profIndex = 0
1103 __withOverapping = False
1114 __withOverapping = False
1104
1115
@@ -1441,7 +1452,9 class IntegrationFaradaySpectra(Operation):
1441
1452
1442
1453
1443 class IntegrationFaradaySpectra2(Operation):
1454 class IntegrationFaradaySpectra2(Operation):
1444
1455 '''
1456 Written by R. Flores
1457 '''
1445 __profIndex = 0
1458 __profIndex = 0
1446 __withOverapping = False
1459 __withOverapping = False
1447
1460
@@ -2061,7 +2074,9 class IntegrationFaradaySpectra2(Operation):
2061 return dataOut
2074 return dataOut
2062
2075
2063 class IntegrationFaradaySpectra3(Operation): #This class should manage data with no lags as well
2076 class IntegrationFaradaySpectra3(Operation): #This class should manage data with no lags as well
2064
2077 '''
2078 Written by R. Flores
2079 '''
2065 __profIndex = 0
2080 __profIndex = 0
2066 __withOverapping = False
2081 __withOverapping = False
2067
2082
@@ -2815,7 +2830,9 class IntegrationFaradaySpectra3(Operation): #This class should manage data with
2815 return dataOut
2830 return dataOut
2816
2831
2817 class IntegrationFaradaySpectraNoLags(Operation):
2832 class IntegrationFaradaySpectraNoLags(Operation):
2818
2833 '''
2834 Written by R. Flores
2835 '''
2819 __profIndex = 0
2836 __profIndex = 0
2820 __withOverapping = False
2837 __withOverapping = False
2821
2838
@@ -3149,6 +3166,9 class IntegrationFaradaySpectraNoLags(Operation):
3149 return dataOut
3166 return dataOut
3150
3167
3151 class HybridSelectSpectra(Operation):
3168 class HybridSelectSpectra(Operation):
3169 '''
3170 Written by R. Flores
3171 '''
3152 """Operation to rearange and use selected channels of spectra data and pairs of cross-spectra data for Hybrid Experiment.
3172 """Operation to rearange and use selected channels of spectra data and pairs of cross-spectra data for Hybrid Experiment.
3153
3173
3154 Parameters:
3174 Parameters:
@@ -3213,7 +3233,9 class HybridSelectSpectra(Operation):
3213
3233
3214
3234
3215 class IncohIntLag(Operation):
3235 class IncohIntLag(Operation):
3216
3236 '''
3237 Written by R. Flores
3238 '''
3217 __profIndex = 0
3239 __profIndex = 0
3218 __withOverapping = False
3240 __withOverapping = False
3219
3241
@@ -3598,6 +3620,7 class IncohInt(Operation):
3598 #exit(1)
3620 #exit(1)
3599
3621
3600 if n == 1:
3622 if n == 1:
3623 dataOut.VelRange = dataOut.getVelRange(0)
3601 return dataOut
3624 return dataOut
3602
3625
3603 dataOut.flagNoData = True
3626 dataOut.flagNoData = True
@@ -3619,6 +3642,9 class IncohInt(Operation):
3619 dataOut.nIncohInt *= self.n
3642 dataOut.nIncohInt *= self.n
3620 dataOut.utctime = avgdatatime
3643 dataOut.utctime = avgdatatime
3621 dataOut.flagNoData = False
3644 dataOut.flagNoData = False
3645
3646 dataOut.VelRange = dataOut.getVelRange(0)
3647
3622 #print("Power",numpy.sum(dataOut.data_spc[0,:,20:30],axis=0))
3648 #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))
3649 #print("Power",numpy.sum(dataOut.data_spc[0,100:110,:],axis=1))
3624 #exit(1)
3650 #exit(1)
@@ -3678,6 +3704,9 class IncohInt(Operation):
3678 return dataOut
3704 return dataOut
3679
3705
3680 class SnrFaraday(Operation):
3706 class SnrFaraday(Operation):
3707 '''
3708 Written by R. Flores
3709 '''
3681 """Operation to use get SNR in Faraday processing.
3710 """Operation to use get SNR in Faraday processing.
3682
3711
3683 Parameters:
3712 Parameters:
@@ -3723,6 +3752,9 class SnrFaraday(Operation):
3723 return dataOut
3752 return dataOut
3724
3753
3725 class SpectraDataToFaraday(Operation):
3754 class SpectraDataToFaraday(Operation):
3755 '''
3756 Written by R. Flores
3757 '''
3726 """Operation to use spectra data in Faraday processing.
3758 """Operation to use spectra data in Faraday processing.
3727
3759
3728 Parameters:
3760 Parameters:
@@ -3754,12 +3786,19 class SpectraDataToFaraday(Operation):
3754 dataOut.year=dataOut.bd_time.tm_year+(dataOut.bd_time.tm_yday-1)/364.0
3786 dataOut.year=dataOut.bd_time.tm_year+(dataOut.bd_time.tm_yday-1)/364.0
3755 dataOut.ut_Faraday=dataOut.bd_time.tm_hour+dataOut.bd_time.tm_min/60.0+dataOut.bd_time.tm_sec/3600.0
3787 dataOut.ut_Faraday=dataOut.bd_time.tm_hour+dataOut.bd_time.tm_min/60.0+dataOut.bd_time.tm_sec/3600.0
3756
3788
3757
3789 '''
3758 tmpx=numpy.zeros((dataOut.nHeights,dataOut.DPL,2),'float32')
3790 tmpx=numpy.zeros((dataOut.nHeights,dataOut.DPL,2),'float32')
3759 tmpx_a2=numpy.zeros((dataOut.nHeights,dataOut.DPL,2),'float32')
3791 tmpx_a2=numpy.zeros((dataOut.nHeights,dataOut.DPL,2),'float32')
3760 tmpx_b2=numpy.zeros((dataOut.nHeights,dataOut.DPL,2),'float32')
3792 tmpx_b2=numpy.zeros((dataOut.nHeights,dataOut.DPL,2),'float32')
3761 tmpx_abr=numpy.zeros((dataOut.nHeights,dataOut.DPL,2),'float32')
3793 tmpx_abr=numpy.zeros((dataOut.nHeights,dataOut.DPL,2),'float32')
3762 tmpx_abi=numpy.zeros((dataOut.nHeights,dataOut.DPL,2),'float32')
3794 tmpx_abi=numpy.zeros((dataOut.nHeights,dataOut.DPL,2),'float32')
3795 '''
3796 #print("NDP",dataOut.NDP)
3797 tmpx=numpy.zeros((dataOut.NDP,dataOut.DPL,2),'float32')
3798 tmpx_a2=numpy.zeros((dataOut.NDP,dataOut.DPL,2),'float32')
3799 tmpx_b2=numpy.zeros((dataOut.NDP,dataOut.DPL,2),'float32')
3800 tmpx_abr=numpy.zeros((dataOut.NDP,dataOut.DPL,2),'float32')
3801 tmpx_abi=numpy.zeros((dataOut.NDP,dataOut.DPL,2),'float32')
3763 dataOut.kabxys_integrated=[tmpx,tmpx,tmpx,tmpx,tmpx_a2,tmpx,tmpx_b2,tmpx,tmpx_abr,tmpx,tmpx_abi,tmpx,tmpx,tmpx]
3802 dataOut.kabxys_integrated=[tmpx,tmpx,tmpx,tmpx,tmpx_a2,tmpx,tmpx_b2,tmpx,tmpx_abr,tmpx,tmpx_abi,tmpx,tmpx,tmpx]
3764 '''
3803 '''
3765 dataOut.rnint2=numpy.zeros(dataOut.DPL,'float32')
3804 dataOut.rnint2=numpy.zeros(dataOut.DPL,'float32')
@@ -3806,25 +3845,7 class SpectraDataToFaraday(Operation):
3806 for l in range(dataOut.DPL):
3845 for l in range(dataOut.DPL):
3807 dataOut.rnint2[l]=1.0/(dataOut.nIncohInt*dataOut.nProfiles)
3846 dataOut.rnint2[l]=1.0/(dataOut.nIncohInt*dataOut.nProfiles)
3808
3847
3809
3848 def noise(self,dataOut):
3810 def run(self,dataOut):
3811 #print(dataOut.nIncohInt)
3812 #exit(1)
3813 dataOut.paramInterval=dataOut.nIncohInt*2*2#nIncohInt*numero de fft/nprofiles*segundos de cada muestra
3814 dataOut.lat=-11.95
3815 dataOut.lon=-76.87
3816
3817 self.normFactor(dataOut)
3818
3819 dataOut.NDP=dataOut.nHeights
3820 dataOut.NR=len(dataOut.channelList)
3821 dataOut.DH=dataOut.heightList[1]-dataOut.heightList[0]
3822 dataOut.H0=int(dataOut.heightList[0])
3823
3824 self.ConvertData(dataOut)
3825
3826 dataOut.NAVG=16#dataOut.rnint2[0] #CHECK THIS!
3827 dataOut.MAXNRANGENDT=dataOut.NDP
3828
3849
3829 dataOut.noise_lag = numpy.zeros((dataOut.nChannels,dataOut.DPL),'float32')
3850 dataOut.noise_lag = numpy.zeros((dataOut.nChannels,dataOut.DPL),'float32')
3830 #print("Lags")
3851 #print("Lags")
@@ -3878,6 +3899,114 class SpectraDataToFaraday(Operation):
3878 dataOut.pan = dataOut.tnoise[0]
3899 dataOut.pan = dataOut.tnoise[0]
3879 dataOut.pbn = dataOut.tnoise[1]
3900 dataOut.pbn = dataOut.tnoise[1]
3880
3901
3902 def get_eej_index_V0(self,data_to_remov_eej,dataOut):
3903
3904 dataOut.data_spc = data_to_remov_eej
3905 #print(dataOut.data_spc)
3906 data_eej = dataOut.getPower()[1]
3907 print("data_eej: ", data_eej)
3908 #exit(1)
3909 index_eej = CleanCohEchoes.mad_based_outlier(self,data_eej[:20])
3910 aux_eej = numpy.array(index_eej.nonzero()).ravel()
3911
3912 index2 = CleanCohEchoes.mad_based_outlier(self,data_eej[aux_eej[-1]+1:aux_eej[-1]+1+20])
3913 aux2 = numpy.array(index2.nonzero()).ravel()
3914 if aux2.size > 0:
3915 #print(aux2)
3916 #print(aux2[-1])
3917 #print(arr[aux[-1]+aux2[-1]+1])
3918 dataOut.min_id_eej = aux_eej[-1]+aux2[-1]+1
3919 else:
3920 dataOut.min_id_eej = aux_eej[-1]
3921
3922
3923 print(dataOut.min_id_eej)
3924 exit(1)
3925
3926 def get_eej_index_V1(self,data_to_remov_eej,dataOut):
3927
3928 dataOut.data_spc = data_to_remov_eej
3929 outliers_IDs = []
3930 #print(dataOut.data_spc)
3931 for ich in range(dataOut.nChannels):
3932
3933 data_eej = dataOut.getPower()[ich]
3934 #print("data_eej: ", data_eej)
3935 #exit(1)
3936 index_eej = CleanCohEchoes.mad_based_outlier(self,data_eej[:20])
3937 aux_eej = numpy.array(index_eej.nonzero()).ravel()
3938
3939 #index2 = CleanCohEchoes.mad_based_outlier(self,data_eej[aux_eej[-1]+1:aux_eej[-1]+1+20])
3940 index2 = CleanCohEchoes.mad_based_outlier(self,data_eej[aux_eej[-1]+1:aux_eej[-1]+1+10],thresh=1.)
3941 aux2 = numpy.array(index2.nonzero()).ravel()
3942 if aux2.size > 0:
3943 #min_id_eej = aux_eej[-1]+aux2[-1]+1
3944 ids = numpy.concatenate((aux_eej,aux2+aux_eej[-1]+1))
3945 else:
3946 ids = aux_eej
3947
3948 outliers_IDs=numpy.append(outliers_IDs,ids)
3949
3950 outliers_IDs=numpy.array(outliers_IDs)
3951 outliers_IDs=outliers_IDs.astype(numpy.dtype('int64'))
3952
3953 (uniq, freq) = (numpy.unique(outliers_IDs, return_counts=True))
3954 aux_arr = numpy.column_stack((uniq,freq))
3955
3956 final_index = []
3957 for i in range(aux_arr.shape[0]):
3958 if aux_arr[i,1] == 2:
3959 final_index.append(aux_arr[i,0])
3960
3961 if final_index != []:
3962 dataOut.min_id_eej = final_index[-1]
3963 else:
3964 print("CHECKKKKK!!!!!!!!!!!!!!!")
3965
3966 print(dataOut.min_id_eej)
3967 exit(1)
3968
3969 def get_eej_index(self,data_to_remov_eej,dataOut):
3970
3971 dataOut.data_spc = data_to_remov_eej
3972
3973 data_eej = dataOut.getPower()[0]
3974 #print(data_eej)
3975 index_eej = CleanCohEchoes.mad_based_outlier(self,data_eej[:17])
3976 aux_eej = numpy.array(index_eej.nonzero()).ravel()
3977
3978 dataOut.min_id_eej = aux_eej[-1]
3979
3980 print(dataOut.min_id_eej)
3981 #exit(1)
3982
3983 def run(self,dataOut):
3984 #print(dataOut.nIncohInt)
3985 #exit(1)
3986 dataOut.paramInterval=dataOut.nIncohInt*2*2#nIncohInt*numero de fft/nprofiles*segundos de cada muestra
3987 dataOut.lat=-11.95
3988 dataOut.lon=-76.87
3989
3990 data_to_remov_eej = dataOut.dataLag_spc[:,:,:,0]
3991
3992 self.normFactor(dataOut)
3993
3994 dataOut.NDP=dataOut.nHeights
3995 dataOut.NR=len(dataOut.channelList)
3996 dataOut.DH=dataOut.heightList[1]-dataOut.heightList[0]
3997 dataOut.H0=int(dataOut.heightList[0])
3998
3999 self.ConvertData(dataOut)
4000
4001 dataOut.NAVG=16#dataOut.rnint2[0] #CHECK THIS!
4002 if hasattr(dataOut, 'NRANGE'):
4003 dataOut.MAXNRANGENDT = max(dataOut.NRANGE,dataOut.NDT)
4004 else:
4005 dataOut.MAXNRANGENDT = dataOut.NDP
4006
4007 if not hasattr(dataOut, 'tnoise'):
4008 self.noise(dataOut)
4009
3881 #dataOut.pan = numpy.mean(dataOut.pan)
4010 #dataOut.pan = numpy.mean(dataOut.pan)
3882 #dataOut.pbn = numpy.mean(dataOut.pbn)
4011 #dataOut.pbn = numpy.mean(dataOut.pbn)
3883 #print(dataOut.pan)
4012 #print(dataOut.pan)
@@ -3888,12 +4017,18 class SpectraDataToFaraday(Operation):
3888 #print("Noise dB: ",10*numpy.log10(dataOut.tnoise))
4017 #print("Noise dB: ",10*numpy.log10(dataOut.tnoise))
3889 #exit(1)
4018 #exit(1)
3890 #dataOut.pan=dataOut.tnoise[0]/float(dataOut.nProfiles_LP*dataOut.nIncohInt)
4019 #dataOut.pan=dataOut.tnoise[0]/float(dataOut.nProfiles_LP*dataOut.nIncohInt)
4020 if gmtime(dataOut.utctime).tm_hour >= 22. or gmtime(dataOut.utctime).tm_hour < 12.:
4021 self.get_eej_index(data_to_remov_eej,dataOut)
3891 print("done")
4022 print("done")
4023 #exit(1)
3892 return dataOut
4024 return dataOut
3893
4025
3894
4026
3895
4027
3896 class SpectraDataToHybrid(SpectraDataToFaraday):
4028 class SpectraDataToHybrid(SpectraDataToFaraday):
4029 '''
4030 Written by R. Flores
4031 '''
3897 """Operation to use spectra data in Faraday processing.
4032 """Operation to use spectra data in Faraday processing.
3898
4033
3899 Parameters:
4034 Parameters:
@@ -4067,6 +4202,9 class SpectraDataToHybrid(SpectraDataToFaraday):
4067 return dataOut
4202 return dataOut
4068
4203
4069 class SpectraDataToHybrid_V2(SpectraDataToFaraday):
4204 class SpectraDataToHybrid_V2(SpectraDataToFaraday):
4205 '''
4206 Written by R. Flores
4207 '''
4070 """Operation to use spectra data in Faraday processing.
4208 """Operation to use spectra data in Faraday processing.
4071
4209
4072 Parameters:
4210 Parameters:
@@ -4092,7 +4230,7 class SpectraDataToHybrid_V2(SpectraDataToFaraday):
4092 self.dataLag_cspc_LP=None
4230 self.dataLag_cspc_LP=None
4093 self.dataLag_dc_LP=None
4231 self.dataLag_dc_LP=None
4094
4232
4095 def noise(self,dataOut):
4233 def noise_v0(self,dataOut):
4096
4234
4097 dataOut.data_spc = dataOut.dataLag_spc_LP.real
4235 dataOut.data_spc = dataOut.dataLag_spc_LP.real
4098 #print(dataOut.dataLag_spc.shape)
4236 #print(dataOut.dataLag_spc.shape)
@@ -4115,9 +4253,104 class SpectraDataToHybrid_V2(SpectraDataToFaraday):
4115 #print("pbn: ",dataOut.pbn)
4253 #print("pbn: ",dataOut.pbn)
4116 #print(numpy.shape(dataOut.pnoise))
4254 #print(numpy.shape(dataOut.pnoise))
4117 #exit(1)
4255 #exit(1)
4118 #print("pan: ",numpy.sum(dataOut.pan))
4256 #print("pan: ",dataOut.pan)
4257 #print("pbn: ",dataOut.pbn)
4119 #exit(1)
4258 #exit(1)
4120
4259
4260 def noise_v0_aux(self,dataOut):
4261
4262 dataOut.data_spc = dataOut.dataLag_spc
4263 #print(dataOut.dataLag_spc.shape)
4264 #exit(1)
4265 #dataOut.data_spc = dataOut.dataLag_spc[:,:,:,0].real
4266 #print("spc noise shape: ",dataOut.data_spc.shape)
4267 tnoise = dataOut.getNoise(ymin_index=100,ymax_index=166)
4268 #print("Noise LP: ",10*numpy.log10(dataOut.tnoise))
4269 #exit(1)
4270 #dataOut.tnoise[0]*=0.995#0.976
4271 #dataOut.tnoise[1]*=0.995
4272 #print(dataOut.nProfiles)
4273 #dataOut.pan=dataOut.tnoise[0]/float(dataOut.nProfiles_LP*dataOut.nIncohInt)
4274 #dataOut.pbn=dataOut.tnoise[1]/float(dataOut.nProfiles_LP*dataOut.nIncohInt)
4275 dataOut.pan=tnoise[0]/float(dataOut.nProfiles*dataOut.nIncohInt)
4276 dataOut.pbn=tnoise[1]/float(dataOut.nProfiles*dataOut.nIncohInt)
4277
4278 def noise(self,dataOut):
4279
4280 dataOut.noise_lag = numpy.zeros((dataOut.nChannels,dataOut.DPL),'float32')
4281 #print("Lags")
4282 '''
4283 for lag in range(dataOut.DPL):
4284 #print(lag)
4285 dataOut.data_spc = dataOut.dataLag_spc[:,:,:,lag]
4286 dataOut.noise_lag[:,lag] = dataOut.getNoise(ymin_index=46)
4287 #dataOut.noise_lag[:,lag] = dataOut.getNoise(ymin_index=33,ymax_index=46)
4288 '''
4289 #print(dataOut.NDP)
4290 #exit(1)
4291 #Channel B
4292 for lag in range(dataOut.DPL):
4293 #print(lag)
4294 dataOut.data_spc = dataOut.dataLag_spc[:,:,:,lag]
4295 max_hei_id = dataOut.NDP - 2*lag
4296 #if lag < 6:
4297 dataOut.noise_lag[1,lag] = dataOut.getNoise(ymin_index=53,ymax_index=max_hei_id)[1]
4298 #else:
4299 #dataOut.noise_lag[1,lag] = numpy.mean(dataOut.noise_lag[1,:6])
4300 #dataOut.noise_lag[:,lag] = dataOut.getNoise(ymin_index=33,ymax_index=46)
4301 #Channel A
4302 for lag in range(dataOut.DPL):
4303 #print(lag)
4304 dataOut.data_spc = dataOut.dataLag_spc[:,:,:,lag]
4305 dataOut.noise_lag[0,lag] = dataOut.getNoise(ymin_index=53)[0]
4306
4307 nanindex = numpy.argwhere(numpy.isnan(numpy.log10(dataOut.noise_lag[1,:])))
4308 i1 = nanindex[0][0]
4309 dataOut.noise_lag[1,(1,2,7,8,9,10)] *= 2 #Correction LP
4310 dataOut.noise_lag[1,i1:] = numpy.mean(dataOut.noise_lag[1,:i1]) #El ruido de lags contaminados se
4311 #determina a partir del promedio del
4312 #ruido de los lags limpios
4313 '''
4314 dataOut.noise_lag[1,:] = dataOut.noise_lag[1,0] #El ruido de los lags diferentes de cero para
4315 #el canal B es contaminado por el Tx y EEJ
4316 #del siguiente perfil, por ello se asigna el ruido
4317 #del lag 0 a todos los lags
4318 '''
4319 #print("Noise lag: ", 10*numpy.log10(dataOut.noise_lag/dataOut.normFactor))
4320 #exit(1)
4321 '''
4322 dataOut.tnoise = dataOut.getNoise(ymin_index=46)
4323 dataOut.tnoise /= float(dataOut.nProfiles*dataOut.nIncohInt)
4324 dataOut.pan = dataOut.tnoise[0]
4325 dataOut.pbn = dataOut.tnoise[1]
4326 '''
4327 #print("i1: ", i1)
4328 #exit(1)
4329 tnoise = dataOut.noise_lag/float(dataOut.nProfiles*dataOut.nIncohInt)
4330 #dataOut.tnoise /= float(dataOut.nProfiles*dataOut.nIncohInt)
4331 dataOut.pan = tnoise[0]
4332 dataOut.pbn = tnoise[1]
4333
4334 def noise_LP(self,dataOut):
4335
4336 dataOut.data_spc = dataOut.dataLag_spc_LP.real
4337 #print(dataOut.dataLag_spc.shape)
4338 #exit(1)
4339 #dataOut.data_spc = dataOut.dataLag_spc[:,:,:,0].real
4340 #print("spc noise shape: ",dataOut.data_spc.shape)
4341 dataOut.tnoise = dataOut.getNoise(ymin_index=100,ymax_index=166)
4342 #print("Noise LP: ",10*numpy.log10(dataOut.tnoise))
4343 #exit(1)
4344 #dataOut.tnoise[0]*=0.995#0.976
4345 #dataOut.tnoise[1]*=0.995
4346 #print(dataOut.nProfiles)
4347 #dataOut.pan=dataOut.tnoise[0]/float(dataOut.nProfiles_LP*dataOut.nIncohInt)
4348 #dataOut.pbn=dataOut.tnoise[1]/float(dataOut.nProfiles_LP*dataOut.nIncohInt)
4349 ######dataOut.pan=dataOut.tnoise[0]/float(dataOut.nProfiles_LP*dataOut.nIncohInt_LP)
4350 ######dataOut.pbn=dataOut.tnoise[1]/float(dataOut.nProfiles_LP*dataOut.nIncohInt_LP)
4351 dataOut.pan_LP=dataOut.tnoise[0]/float(dataOut.nProfiles_LP*dataOut.nIncohInt_LP)
4352 dataOut.pbn_LP=dataOut.tnoise[1]/float(dataOut.nProfiles_LP*dataOut.nIncohInt_LP)
4353
4121 def ConvertDataLP(self,dataOut):
4354 def ConvertDataLP(self,dataOut):
4122
4355
4123 #print(dataOut.dataLag_spc[:,:,:,1]/dataOut.data_spc)
4356 #print(dataOut.dataLag_spc[:,:,:,1]/dataOut.data_spc)
@@ -4161,6 +4394,7 class SpectraDataToHybrid_V2(SpectraDataToFaraday):
4161 #dataOut.output_LP_integrated[:,:,3] *= float(dataOut.NSCAN/22)#(dataOut.nNoiseProfiles) #Corrects the zero padding
4394 #dataOut.output_LP_integrated[:,:,3] *= float(dataOut.NSCAN/22)#(dataOut.nNoiseProfiles) #Corrects the zero padding
4162
4395
4163 dataOut.nis=dataOut.NSCAN*dataOut.nIncohInt_LP*10
4396 dataOut.nis=dataOut.NSCAN*dataOut.nIncohInt_LP*10
4397 dataOut.nis=dataOut.NSCAN*dataOut.nIncohInt_LP*dataOut.nProfiles_LP*10
4164
4398
4165 self.ConvertData(dataOut)
4399 self.ConvertData(dataOut)
4166
4400
@@ -4170,10 +4404,98 class SpectraDataToHybrid_V2(SpectraDataToFaraday):
4170 dataOut.kabxys_integrated[10][:,(1,2,7,8,9,10),0] *= 2 #Corrects the zero padding
4404 dataOut.kabxys_integrated[10][:,(1,2,7,8,9,10),0] *= 2 #Corrects the zero padding
4171 hei = 2
4405 hei = 2
4172
4406
4173 self.noise(dataOut)
4407 self.noise(dataOut) #Noise for DP Profiles
4408 dataOut.pan[[1,2,7,8,9,10]] *= 2 #Corrects the zero padding
4409 #dataOut.pbn[[1,2,7,8,9,10]] *= 2 #Corrects the zero padding #Chequear debido a que se están mezclando lags en self.noise()
4410 self.noise_LP(dataOut) #Noise for LP Profiles
4411
4412 print("pan: , pan_LP: ",dataOut.pan,dataOut.pan_LP)
4413 print("pbn: , pbn_LP: ",dataOut.pbn,dataOut.pbn_LP)
4414
4415
4174
4416
4175 dataOut.NAVG=1#dataOut.rnint2[0] #CHECK THIS!
4417 dataOut.NAVG=1#dataOut.rnint2[0] #CHECK THIS!
4176 dataOut.nint=dataOut.nIncohInt
4418 dataOut.nint=dataOut.nIncohInt
4177 dataOut.MAXNRANGENDT=dataOut.output_LP_integrated.shape[1]
4419 dataOut.MAXNRANGENDT=dataOut.output_LP_integrated.shape[1]
4178
4420
4421 '''
4422 range_aux=numpy.zeros(dataOut.MAXNRANGENDT,order='F',dtype='float32')
4423 range_aux_dp=numpy.zeros(dataOut.NDT,order='F',dtype='float32')
4424 for i in range(dataOut.MAXNRANGENDT):
4425 range_aux[i]=dataOut.H0 + i*dataOut.DH
4426 for i in range(dataOut.NDT):
4427 range_aux_dp[i]=dataOut.H0 + i*dataOut.DH
4428 import matplotlib.pyplot as plt
4429 #plt.plot(10*numpy.log10(dataOut.output_LP_integrated.real[0,:,0]),range_aux)
4430 plt.plot(10*numpy.log10(dataOut.output_LP_integrated.real[0,:,0]),range_aux_dp)
4431 #plt.plot(10*numpy.log10(dataOut.output_LP_integrated.real[0,:,0]/dataOut.nProfiles_LP),dataOut.range1)
4432 plt.axvline(10*numpy.log10(dataOut.tnoise[0]),color='k',linestyle='dashed')
4433 plt.grid()
4434 plt.xlim(20,100)
4435 plt.show()
4436 exit(1)
4437 '''
4438 return dataOut
4439
4440 class SpcVoltageDataToHybrid(SpectraDataToFaraday):
4441 '''
4442 Written by R. Flores
4443 '''
4444 """Operation to use spectra data in Faraday processing.
4445
4446 Parameters:
4447 -----------
4448 nint : int
4449 Number of integrations.
4450
4451 Example
4452 --------
4453
4454 op = proc_unit.addOperation(name='SpcVoltageDataToHybrid', optype='other')
4455
4456 """
4457
4458 def __init__(self, **kwargs):
4459
4460 Operation.__init__(self, **kwargs)
4461
4462 self.dataLag_spc=None
4463 self.dataLag_cspc=None
4464 self.dataLag_dc=None
4465
4466 def normFactor(self,dataOut):
4467 dataOut.rnint2=numpy.zeros(dataOut.DPL,'float32')
4468 #print(dataOut.nIncohInt,dataOut.nProfiles)
4469 for l in range(dataOut.DPL):
4470 if(l==0 or (l>=3 and l <=6)):
4471 dataOut.rnint2[l]=1.0/(dataOut.nIncohInt*dataOut.nProfiles_DP)
4472 else:
4473 dataOut.rnint2[l]=2*(1.0/(dataOut.nIncohInt*dataOut.nProfiles_DP))
4474
4475 def run(self,dataOut):
4476
4477 dataOut.paramInterval=0#int(dataOut.nint*dataOut.header[7][0]*2 )
4478 dataOut.lat=-11.95
4479 dataOut.lon=-76.87
4480
4481 #dataOut.NDP=dataOut.nHeights
4482 #dataOut.NR=len(dataOut.channelList)
4483 #dataOut.DH=dataOut.heightList[1]-dataOut.heightList[0]
4484 #dataOut.H0=int(dataOut.heightList[0])
4485
4486 self.normFactor(dataOut)
4487
4488 dataOut.nis=dataOut.NSCAN*dataOut.NAVG*dataOut.nint*10
4489
4490 self.ConvertData(dataOut)
4491
4492 dataOut.kabxys_integrated[4][:,(1,2,7,8,9,10),0] *= 2 #Corrects the zero padding
4493 dataOut.kabxys_integrated[6][:,(1,2,7,8,9,10),0] *= 2 #Corrects the zero padding
4494 dataOut.kabxys_integrated[8][:,(1,2,7,8,9,10),0] *= 2 #Corrects the zero padding
4495 dataOut.kabxys_integrated[10][:,(1,2,7,8,9,10),0] *= 2 #Corrects the zero padding
4496 #print(numpy.sum(dataOut.kabxys_integrated[4][:,1,0]))
4497 dataOut.MAXNRANGENDT = max(dataOut.NRANGE,dataOut.NDP)
4498 #print(dataOut.rnint2)
4499 #print(dataOut.nis)
4500 #exit(1)
4179 return dataOut
4501 return dataOut
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (508 lines changed) Show them Hide them
@@ -29,13 +29,15 class VoltageProc(ProcessingUnit):
29 self.flip = 1
29 self.flip = 1
30 self.setupReq = False
30 self.setupReq = False
31
31
32 def run(self):
32 def run(self, runNextUnit = 0):
33
33
34 if self.dataIn.type == 'AMISR':
34 if self.dataIn.type == 'AMISR':
35 self.__updateObjFromAmisrInput()
35 self.__updateObjFromAmisrInput()
36
36
37 if self.dataIn.type == 'Voltage':
37 if self.dataIn.type == 'Voltage':
38 self.dataOut.copy(self.dataIn)
38 self.dataOut.copy(self.dataIn)
39 self.dataOut.runNextUnit = runNextUnit
40 #print(self.dataOut.data.shape)
39
41
40 def __updateObjFromAmisrInput(self):
42 def __updateObjFromAmisrInput(self):
41
43
@@ -397,6 +399,9 class deFlip(Operation):
397 return dataOut
399 return dataOut
398
400
399 class deFlipHP(Operation):
401 class deFlipHP(Operation):
402 '''
403 Written by R. Flores
404 '''
400 def __init__(self):
405 def __init__(self):
401
406
402 self.flip = 1
407 self.flip = 1
@@ -530,6 +535,9 class interpolateHeights(Operation):
530
535
531
536
532 class LagsReshape(Operation):
537 class LagsReshape(Operation):
538 '''
539 Written by R. Flores
540 '''
533 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
541 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
534
542
535 Parameters:
543 Parameters:
@@ -621,6 +629,9 class LagsReshape(Operation):
621 return dataOut
629 return dataOut
622
630
623 class LagsReshape150(Operation):
631 class LagsReshape150(Operation):
632 '''
633 Written by R. Flores
634 '''
624 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
635 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
625
636
626 Parameters:
637 Parameters:
@@ -714,6 +725,9 class LagsReshape150(Operation):
714 return dataOut
725 return dataOut
715
726
716 class LagsReshapeHP(Operation):
727 class LagsReshapeHP(Operation):
728 '''
729 Written by R. Flores
730 '''
717 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
731 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
718
732
719 Parameters:
733 Parameters:
@@ -805,6 +819,9 class LagsReshapeHP(Operation):
805 return dataOut
819 return dataOut
806
820
807 class LagsReshapeDP(Operation):
821 class LagsReshapeDP(Operation):
822 '''
823 Written by R. Flores
824 '''
808 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
825 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
809
826
810 Parameters:
827 Parameters:
@@ -957,6 +974,9 class LagsReshapeDP(Operation):
957 return dataOut
974 return dataOut
958
975
959 class LagsReshapeDP_V2(Operation):
976 class LagsReshapeDP_V2(Operation):
977 '''
978 Written by R. Flores
979 '''
960 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
980 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
961
981
962 Parameters:
982 Parameters:
@@ -1156,6 +1176,9 class LagsReshapeDP_V2(Operation):
1156 return dataOut
1176 return dataOut
1157
1177
1158 class LagsReshapeLP(Operation):
1178 class LagsReshapeLP(Operation):
1179 '''
1180 Written by R. Flores
1181 '''
1159 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
1182 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
1160
1183
1161 Parameters:
1184 Parameters:
@@ -1322,6 +1345,9 class LagsReshapeLP(Operation):
1322 return dataOut
1345 return dataOut
1323
1346
1324 class LagsReshapeHP2(Operation):
1347 class LagsReshapeHP2(Operation):
1348 '''
1349 Written by R. Flores
1350 '''
1325 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
1351 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
1326
1352
1327 Parameters:
1353 Parameters:
@@ -1499,6 +1525,9 class LagsReshapeHP2(Operation):
1499 return dataOut
1525 return dataOut
1500
1526
1501 class CrossProdDP(Operation):
1527 class CrossProdDP(Operation):
1528 '''
1529 Written by R. Flores
1530 '''
1502 """Operation to calculate cross products of the Double Pulse Experiment.
1531 """Operation to calculate cross products of the Double Pulse Experiment.
1503
1532
1504 Parameters:
1533 Parameters:
@@ -1983,6 +2012,9 class CrossProdDP(Operation):
1983
2012
1984
2013
1985 class IntegrationDP(Operation):
2014 class IntegrationDP(Operation):
2015 '''
2016 Written by R. Flores
2017 '''
1986 """Operation to integrate the Double Pulse data.
2018 """Operation to integrate the Double Pulse data.
1987
2019
1988 Parameters:
2020 Parameters:
@@ -2068,6 +2100,9 class IntegrationDP(Operation):
2068
2100
2069
2101
2070 class SumFlips(Operation):
2102 class SumFlips(Operation):
2103 '''
2104 Written by R. Flores
2105 '''
2071 """Operation to sum the flip and unflip part of certain cross products of the Double Pulse.
2106 """Operation to sum the flip and unflip part of certain cross products of the Double Pulse.
2072
2107
2073 Parameters:
2108 Parameters:
@@ -2128,6 +2163,9 class SumFlips(Operation):
2128
2163
2129
2164
2130 class FlagBadHeights(Operation):
2165 class FlagBadHeights(Operation):
2166 '''
2167 Written by R. Flores
2168 '''
2131 """Operation to flag bad heights (bad data) of the Double Pulse.
2169 """Operation to flag bad heights (bad data) of the Double Pulse.
2132
2170
2133 Parameters:
2171 Parameters:
@@ -2161,6 +2199,9 class FlagBadHeights(Operation):
2161 return dataOut
2199 return dataOut
2162
2200
2163 class FlagBadHeightsSpectra(Operation):
2201 class FlagBadHeightsSpectra(Operation):
2202 '''
2203 Written by R. Flores
2204 '''
2164 """Operation to flag bad heights (bad data) of the Double Pulse.
2205 """Operation to flag bad heights (bad data) of the Double Pulse.
2165
2206
2166 Parameters:
2207 Parameters:
@@ -2194,6 +2235,9 class FlagBadHeightsSpectra(Operation):
2194 return dataOut
2235 return dataOut
2195
2236
2196 class CleanCohEchoes(Operation):
2237 class CleanCohEchoes(Operation):
2238 '''
2239 Written by R. Flores
2240 '''
2197 """Operation to clean coherent echoes.
2241 """Operation to clean coherent echoes.
2198
2242
2199 Parameters:
2243 Parameters:
@@ -2561,6 +2605,9 class CleanCohEchoes(Operation):
2561 return dataOut
2605 return dataOut
2562
2606
2563 class CleanCohEchoesHP(Operation):
2607 class CleanCohEchoesHP(Operation):
2608 '''
2609 Written by R. Flores
2610 '''
2564 """Operation to clean coherent echoes.
2611 """Operation to clean coherent echoes.
2565
2612
2566 Parameters:
2613 Parameters:
@@ -2928,6 +2975,9 class CleanCohEchoesHP(Operation):
2928 return dataOut
2975 return dataOut
2929
2976
2930 class NoisePower(Operation):
2977 class NoisePower(Operation):
2978 '''
2979 Written by R. Flores
2980 '''
2931 """Operation to get noise power from the integrated data of the Double Pulse.
2981 """Operation to get noise power from the integrated data of the Double Pulse.
2932
2982
2933 Parameters:
2983 Parameters:
@@ -3012,6 +3062,9 class NoisePower(Operation):
3012
3062
3013
3063
3014 class DoublePulseACFs(Operation):
3064 class DoublePulseACFs(Operation):
3065 '''
3066 Written by R. Flores
3067 '''
3015 """Operation to get the ACFs of the Double Pulse.
3068 """Operation to get the ACFs of the Double Pulse.
3016
3069
3017 Parameters:
3070 Parameters:
@@ -3066,7 +3119,7 class DoublePulseACFs(Operation):
3066 '''
3119 '''
3067 #print("init 2.6",pa,dataOut.pan)
3120 #print("init 2.6",pa,dataOut.pan)
3068 dataOut.p[i,j]=pa+pb-(dataOut.pan+dataOut.pbn)
3121 dataOut.p[i,j]=pa+pb-(dataOut.pan+dataOut.pbn)
3069
3122 #print(i,j,dataOut.p[i,j])
3070 dataOut.sdp[i,j]=2*dataOut.rnint2[j]*((pa+pb)*(pa+pb))
3123 dataOut.sdp[i,j]=2*dataOut.rnint2[j]*((pa+pb)*(pa+pb))
3071 ## ACF
3124 ## ACF
3072
3125
@@ -3116,11 +3169,11 class DoublePulseACFs(Operation):
3116 '''
3169 '''
3117 #'''
3170 #'''
3118 if ((pb/dataOut.pbn-1.0)>2.25*(pa/dataOut.pan-1.0)): #To flag bad points from the pulse and EEJ for lags != 0 for Channel B
3171 if ((pb/dataOut.pbn-1.0)>2.25*(pa/dataOut.pan-1.0)): #To flag bad points from the pulse and EEJ for lags != 0 for Channel B
3119 #print("EJJ")
3172 #print(dataOut.heightList[i],"EJJ")
3120 dataOut.igcej[i,j]=1
3173 dataOut.igcej[i,j]=1
3121
3174
3122 elif ((pa/dataOut.pan-1.0)>2.25*(pb/dataOut.pbn-1.0)):
3175 elif ((pa/dataOut.pan-1.0)>2.25*(pb/dataOut.pbn-1.0)):
3123 #print("EJJ")
3176 #print(dataOut.heightList[i],"EJJ")
3124 dataOut.igcej[i,j]=1
3177 dataOut.igcej[i,j]=1
3125 #'''
3178 #'''
3126 '''
3179 '''
@@ -3154,10 +3207,14 class DoublePulseACFs(Operation):
3154 #print("p: ",dataOut.p[33,:])
3207 #print("p: ",dataOut.p[33,:])
3155 #exit(1)
3208 #exit(1)
3156 '''
3209 '''
3157
3210 #print(numpy.sum(dataOut.rhor))
3211 #exit(1)
3158 return dataOut
3212 return dataOut
3159
3213
3160 class DoublePulseACFs_PerLag(Operation):
3214 class DoublePulseACFs_PerLag(Operation):
3215 '''
3216 Written by R. Flores
3217 '''
3161 """Operation to get the ACFs of the Double Pulse.
3218 """Operation to get the ACFs of the Double Pulse.
3162
3219
3163 Parameters:
3220 Parameters:
@@ -3258,18 +3315,18 class DoublePulseACFs_PerLag(Operation):
3258 '''
3315 '''
3259 #'''
3316 #'''
3260 if ((pb/dataOut.pbn[j]-1.0)>2.25*(pa/dataOut.pan[j]-1.0)): #To flag bad points from the pulse and EEJ for lags != 0 for Channel B
3317 if ((pb/dataOut.pbn[j]-1.0)>2.25*(pa/dataOut.pan[j]-1.0)): #To flag bad points from the pulse and EEJ for lags != 0 for Channel B
3261 #print("EJJ")
3318 #print(dataOut.heightList[i],j,"EJJ")
3262 dataOut.igcej[i,j]=1
3319 dataOut.igcej[i,j]=1
3263
3320
3264 elif ((pa/dataOut.pan[j]-1.0)>2.25*(pb/dataOut.pbn[j]-1.0)):
3321 elif ((pa/dataOut.pan[j]-1.0)>2.25*(pb/dataOut.pbn[j]-1.0)):
3265 #print("EJJ")
3322 #print(dataOut.heightList[i],j,"EJJ")
3266 dataOut.igcej[i,j]=1
3323 dataOut.igcej[i,j]=1
3267 #'''
3324 #'''
3268 '''
3325 '''
3269 if ((pa/dataOut.pan-1.0)>2.25*(pb/dataOut.pbn-1.0)):
3326 if ((pa/dataOut.pan-1.0)>2.25*(pb/dataOut.pbn-1.0)):
3270 #print("EJJ")
3327 #print("EJJ")
3271 dataOut.igcej[i,j]=1
3328 dataOut.igcej[i,j]=1
3272 '''
3329 #'''
3273 '''
3330 '''
3274 if i == 4:
3331 if i == 4:
3275 exit(1)
3332 exit(1)
@@ -3299,6 +3356,9 class DoublePulseACFs_PerLag(Operation):
3299 return dataOut
3356 return dataOut
3300
3357
3301 class FaradayAngleAndDPPower(Operation):
3358 class FaradayAngleAndDPPower(Operation):
3359 '''
3360 Written by R. Flores
3361 '''
3302 """Operation to calculate Faraday angle and Double Pulse power.
3362 """Operation to calculate Faraday angle and Double Pulse power.
3303
3363
3304 Parameters:
3364 Parameters:
@@ -3345,6 +3405,7 class FaradayAngleAndDPPower(Operation):
3345 st=0.# // total signal
3405 st=0.# // total signal
3346 ibt=0# // bad lags
3406 ibt=0# // bad lags
3347 ns=0# // no. good lags
3407 ns=0# // no. good lags
3408 #print(dataOut.heightList[j])
3348 for l in range(dataOut.DPL):
3409 for l in range(dataOut.DPL):
3349 #add in other lags if outside of e-jet contamination
3410 #add in other lags if outside of e-jet contamination
3350 if( (dataOut.igcej[j][l] == 0) and (dataOut.ibad[j][l] == 0) ):
3411 if( (dataOut.igcej[j][l] == 0) and (dataOut.ibad[j][l] == 0) ):
@@ -3352,8 +3413,8 class FaradayAngleAndDPPower(Operation):
3352 dataOut.ph2[j]+=dataOut.p[j][l]/dataOut.sdp[j][l]
3413 dataOut.ph2[j]+=dataOut.p[j][l]/dataOut.sdp[j][l]
3353 dataOut.sdp2[j]=dataOut.sdp2[j]+1./dataOut.sdp[j][l]
3414 dataOut.sdp2[j]=dataOut.sdp2[j]+1./dataOut.sdp[j][l]
3354 ns+=1
3415 ns+=1
3355
3416 #if dataOut.igcej[j][l] != 0:
3356
3417 #print(l)
3357 pt+=dataOut.p[j][l]/dataOut.sdp[j][l]
3418 pt+=dataOut.p[j][l]/dataOut.sdp[j][l]
3358 st+=1./dataOut.sdp[j][l]
3419 st+=1./dataOut.sdp[j][l]
3359 ibt|=dataOut.ibad[j][l];
3420 ibt|=dataOut.ibad[j][l];
@@ -3387,6 +3448,9 class FaradayAngleAndDPPower(Operation):
3387
3448
3388
3449
3389 class ElectronDensityFaraday(Operation):
3450 class ElectronDensityFaraday(Operation):
3451 '''
3452 Written by R. Flores
3453 '''
3390 """Operation to calculate electron density from Faraday angle.
3454 """Operation to calculate electron density from Faraday angle.
3391
3455
3392 Parameters:
3456 Parameters:
@@ -3433,6 +3497,7 class ElectronDensityFaraday(Operation):
3433 ndphi=dataOut.NSHTS-4
3497 ndphi=dataOut.NSHTS-4
3434 #print(dataOut.phi)
3498 #print(dataOut.phi)
3435 #exit(1)
3499 #exit(1)
3500 #'''
3436 if dataOut.flagSpreadF:
3501 if dataOut.flagSpreadF:
3437 nanindex = numpy.argwhere(numpy.isnan(dataOut.phi))
3502 nanindex = numpy.argwhere(numpy.isnan(dataOut.phi))
3438 i1 = nanindex[-1][0]
3503 i1 = nanindex[-1][0]
@@ -3443,6 +3508,7 class ElectronDensityFaraday(Operation):
3443 else:
3508 else:
3444 #dataOut.phi_uwrp = dataOut.phi.copy()
3509 #dataOut.phi_uwrp = dataOut.phi.copy()
3445 dataOut.phi[:]=numpy.unwrap(dataOut.phi[:]) #Better results
3510 dataOut.phi[:]=numpy.unwrap(dataOut.phi[:]) #Better results
3511 #'''
3446 #print(dataOut.phi)
3512 #print(dataOut.phi)
3447 #print(dataOut.ph2)
3513 #print(dataOut.ph2)
3448 #exit(1)
3514 #exit(1)
@@ -3455,9 +3521,11 class ElectronDensityFaraday(Operation):
3455 for i in range(2,dataOut.NSHTS-2):
3521 for i in range(2,dataOut.NSHTS-2):
3456 fact=(-0.5/(dataOut.RATE*dataOut.DH))*dataOut.bki[i]
3522 fact=(-0.5/(dataOut.RATE*dataOut.DH))*dataOut.bki[i]
3457 #four-point derivative, no phase unwrapping necessary
3523 #four-point derivative, no phase unwrapping necessary
3458 ####dataOut.dphi[i]=((((theta[i+1]-theta[i-1])+(2.0*(theta[i+2]-theta[i-2])))/thetai[i])).real/10.0
3524 #####dataOut.dphi[i]=((((theta[i+1]-theta[i-1])+(2.0*(theta[i+2]-theta[i-2])))/thetai[i])).real/10.0 #Original from C program
3525
3459 ##dataOut.dphi[i]=((((theta[i-2]-theta[i+2])+(8.0*(theta[i+1]-theta[i-1])))/thetai[i])).real/12.0
3526 ##dataOut.dphi[i]=((((theta[i-2]-theta[i+2])+(8.0*(theta[i+1]-theta[i-1])))/thetai[i])).real/12.0
3460 dataOut.dphi[i]=((dataOut.phi[i+1]-dataOut.phi[i-1])+(2.0*(dataOut.phi[i+2]-dataOut.phi[i-2])))/10.0 #Better results
3527 dataOut.dphi[i]=((dataOut.phi[i+1]-dataOut.phi[i-1])+(2.0*(dataOut.phi[i+2]-dataOut.phi[i-2])))/10.0 #Better results
3528
3461 #dataOut.dphi_uc[i] = abs(dataOut.phi[i]*dataOut.bki[i]*(-0.5)/dataOut.DH)
3529 #dataOut.dphi_uc[i] = abs(dataOut.phi[i]*dataOut.bki[i]*(-0.5)/dataOut.DH)
3462 dataOut.dphi[i]=abs(dataOut.dphi[i]*fact)
3530 dataOut.dphi[i]=abs(dataOut.dphi[i]*fact)
3463 dataOut.sdn1[i]=(4.*(dataOut.sdn2[i-2]+dataOut.sdn2[i+2])+dataOut.sdn2[i-1]+dataOut.sdn2[i+1])
3531 dataOut.sdn1[i]=(4.*(dataOut.sdn2[i-2]+dataOut.sdn2[i+2])+dataOut.sdn2[i-1]+dataOut.sdn2[i+1])
@@ -3468,7 +3536,10 class ElectronDensityFaraday(Operation):
3468
3536
3469
3537
3470 class NormalizeDPPower(Operation):
3538 class NormalizeDPPower(Operation):
3471 """Operation to normalize relative electron density from power with total electron density from Farday angle.
3539 '''
3540 Written by R. Flores
3541 '''
3542 """Operation to normalize relative electron density from power with total electron density from Faraday angle.
3472
3543
3473 Parameters:
3544 Parameters:
3474 -----------
3545 -----------
@@ -3602,6 +3673,9 class NormalizeDPPower(Operation):
3602 return dataOut
3673 return dataOut
3603
3674
3604 class NormalizeDPPowerRoberto(Operation):
3675 class NormalizeDPPowerRoberto(Operation):
3676 '''
3677 Written by R. Flores
3678 '''
3605 """Operation to normalize relative electron density from power with total electron density from Farday angle.
3679 """Operation to normalize relative electron density from power with total electron density from Farday angle.
3606
3680
3607 Parameters:
3681 Parameters:
@@ -3738,6 +3812,9 class NormalizeDPPowerRoberto(Operation):
3738 return dataOut
3812 return dataOut
3739
3813
3740 class NormalizeDPPowerRoberto_V2(Operation):
3814 class NormalizeDPPowerRoberto_V2(Operation):
3815 '''
3816 Written by R. Flores
3817 '''
3741 """Operation to normalize relative electron density from power with total electron density from Farday angle.
3818 """Operation to normalize relative electron density from power with total electron density from Farday angle.
3742
3819
3743 Parameters:
3820 Parameters:
@@ -3799,7 +3876,7 class NormalizeDPPowerRoberto_V2(Operation):
3799 self.aux=0
3876 self.aux=0
3800
3877
3801 night_first=250.0
3878 night_first=250.0
3802 night_first1= 400.0
3879 night_first1= 350.0
3803 night_end= 450.0
3880 night_end= 450.0
3804 day_first=220.0
3881 day_first=220.0
3805 day_end=400.0
3882 day_end=400.0
@@ -3812,7 +3889,7 class NormalizeDPPowerRoberto_V2(Operation):
3812 #print("EARLY")
3889 #print("EARLY")
3813 i2=(night_end-dataOut.range1[0])/dataOut.DH
3890 i2=(night_end-dataOut.range1[0])/dataOut.DH
3814 i1=(night_first -dataOut.range1[0])/dataOut.DH
3891 i1=(night_first -dataOut.range1[0])/dataOut.DH
3815 elif (dataOut.ut_Faraday>=23.0 or dataOut.ut_Faraday<2.0):
3892 elif (dataOut.ut_Faraday>=23.0 or dataOut.ut_Faraday<=2.0):
3816 #print("NIGHT")
3893 #print("NIGHT")
3817 i2=(night_end-dataOut.range1[0])/dataOut.DH
3894 i2=(night_end-dataOut.range1[0])/dataOut.DH
3818 i1=(night_first1 -dataOut.range1[0])/dataOut.DH
3895 i1=(night_first1 -dataOut.range1[0])/dataOut.DH
@@ -3847,15 +3924,17 class NormalizeDPPowerRoberto_V2(Operation):
3847 #print("Flag: ",dataOut.flagTeTiCorrection)
3924 #print("Flag: ",dataOut.flagTeTiCorrection)
3848 #print(dataOut.dphi[i1::])
3925 #print(dataOut.dphi[i1::])
3849 #print(dataOut.ph2[:])
3926 #print(dataOut.ph2[:])
3927
3850 if dataOut.flagTeTiCorrection:
3928 if dataOut.flagTeTiCorrection:
3851 for i in range(dataOut.NSHTS):
3929 for i in range(dataOut.NSHTS):
3852 dataOut.ph2[i]/=dataOut.cf
3930 dataOut.ph2[i]/=dataOut.cf
3853 dataOut.sdp2[i]/=dataOut.cf
3931 dataOut.sdp2[i]/=dataOut.cf
3932
3854 #'''
3933 #'''
3855 if dataOut.flagSpreadF:
3934 if dataOut.flagSpreadF:
3856 i2=int((620-dataOut.range1[0])/dataOut.DH)
3935 i2=int((620-dataOut.range1[0])/dataOut.DH)
3857 nanindex = numpy.argwhere(numpy.isnan(dataOut.ph2))
3936 nanindex = numpy.argwhere(numpy.isnan(dataOut.ph2))
3858 print(nanindex)
3937 print("nanindex",nanindex)
3859 i1 = nanindex[-1][0] #VER CUANDO i1>i2
3938 i1 = nanindex[-1][0] #VER CUANDO i1>i2
3860 i1 += 1+2 #Se suma uno para no tomar el nan, se suma 2 para no tomar datos nan de "phi" debido al calculo de la derivada
3939 i1 += 1+2 #Se suma uno para no tomar el nan, se suma 2 para no tomar datos nan de "phi" debido al calculo de la derivada
3861 #print("i1, i2",i1,i2)
3940 #print("i1, i2",i1,i2)
@@ -3887,11 +3966,24 class NormalizeDPPowerRoberto_V2(Operation):
3887 except:
3966 except:
3888 pass
3967 pass
3889
3968
3890 #time_text = datetime.datetime.utcfromtimestamp(dataOut.utctime)
3969 #print(dataOut.cf,dataOut.cflast[0])
3970 time_text = datetime.datetime.utcfromtimestamp(dataOut.utctime)
3971 #'''
3891 #if (time_text.hour == 5 and time_text.minute == 32): #Year: 2022, DOY:104
3972 #if (time_text.hour == 5 and time_text.minute == 32): #Year: 2022, DOY:104
3892 #if (time_text.hour == 0 and time_text.minute == 12): #Year: 2022, DOY:93
3973 #if (time_text.hour == 0 and time_text.minute == 12): #Year: 2022, DOY:93
3893 #dataOut.cf = dataOut.cflast[0]
3974 #if (time_text.hour == 0 and time_text.minute == 22) or (time_text.hour == 0 and time_text.minute == 54) or (time_text.hour == 1 and time_text.minute == 48): #Year: 2022, DOY:242
3894
3975 #if (time_text.hour == 1 and time_text.minute == 23) or (time_text.hour == 1 and time_text.minute == 44): #Year: 2022, DOY:243
3976 if (time_text.hour == 0 and time_text.minute == 4): #Year: 2022, DOY:244
3977 dataOut.cf = dataOut.cflast[0]
3978 #dataOut.cf = 0.08
3979 #print("here")
3980 if (time_text.hour == 2 and time_text.minute == 23): #Year: 2022, DOY:244
3981 dataOut.cf = 0.08
3982 if (time_text.hour == 2 and time_text.minute == 33): #Year: 2022, DOY:244
3983 dataOut.cf = 0.09
3984 if (time_text.hour == 3 and time_text.minute == 59) or (time_text.hour == 4 and time_text.minute == 20): #Year: 2022, DOY:244
3985 dataOut.cf = 0.09
3986 #'''
3895 dataOut.cflast[0]=dataOut.cf
3987 dataOut.cflast[0]=dataOut.cf
3896 #print(dataOut.cf)
3988 #print(dataOut.cf)
3897
3989
@@ -3953,6 +4045,9 class suppress_stdout_stderr(object):
3953
4045
3954
4046
3955 class DPTemperaturesEstimation(Operation):
4047 class DPTemperaturesEstimation(Operation):
4048 '''
4049 Written by R. Flores
4050 '''
3956 """Operation to estimate temperatures for Double Pulse data.
4051 """Operation to estimate temperatures for Double Pulse data.
3957
4052
3958 Parameters:
4053 Parameters:
@@ -4172,7 +4267,9 class DPTemperaturesEstimation(Operation):
4172
4267
4173
4268
4174 class DenCorrection(NormalizeDPPowerRoberto_V2):
4269 class DenCorrection(NormalizeDPPowerRoberto_V2):
4175
4270 '''
4271 Written by R. Flores
4272 '''
4176 def __init__(self, **kwargs):
4273 def __init__(self, **kwargs):
4177
4274
4178 Operation.__init__(self, **kwargs)
4275 Operation.__init__(self, **kwargs)
@@ -4199,42 +4296,109 class DenCorrection(NormalizeDPPowerRoberto_V2):
4199 bline=0.0
4296 bline=0.0
4200 #bline=numpy.zeros(1,order='F',dtype='float32')
4297 #bline=numpy.zeros(1,order='F',dtype='float32')
4201 my_aux = numpy.ones(dataOut.NSHTS,order='F',dtype='float32')
4298 my_aux = numpy.ones(dataOut.NSHTS,order='F',dtype='float32')
4299 acf_Temps = numpy.ones(dataOut.NSHTS,order='F',dtype='float32')*numpy.nan
4300 acf_no_Temps = numpy.ones(dataOut.NSHTS,order='F',dtype='float32')*numpy.nan
4301
4302 from scipy import signal
4303
4304
4305 teti = numpy.ones(dataOut.NSHTS,order='F',dtype='float32')
4306
4307 for i in range(10,26):
4308 teti[i] = dataOut.te2[i]/dataOut.ti2[i]
4309
4310 ratio2 = teti-1
4311
4312 ratio2 = signal.medfilt(ratio2)
4313
4314 def func(params):
4315 return (ratio2-self.gaussian(dataOut.heightList[:dataOut.NSHTS],params[0],params[1],params[2]))
4316 #print(ratio2)
4317
4318 #plt.show()
4319 x0_value = numpy.array([.5,250,20])
4320
4321 popt = least_squares(func,x0=x0_value,verbose=0)
4322
4323 A = popt.x[0]; B = popt.x[1]; C = popt.x[2]
4324
4325 ratio2 = self.gaussian(dataOut.heightList[:dataOut.NSHTS], A, B, C) + 1 #Te/Ti + 1
4326 '''
4327 import matplotlib.pyplot as plt
4328 plt.clf()
4329 plt.plot(teti,dataOut.heightList[:dataOut.NSHTS])
4330 plt.plot(ratio2,dataOut.heightList[:dataOut.NSHTS])
4331 plt.title("{}".format(datetime.datetime.fromtimestamp(dataOut.utctime)))
4332 plt.xlim(.99,3)
4333 plt.grid()
4334 plt.savefig("/home/roberto/Pictures/Density_Comparison/TeTi_from_temps/{}.png".format(dataOut.utctime))
4335 '''
4336
4337 my_te2 = dataOut.ti2*ratio2
4338 #'''
4339 def func(params):
4340 return (dataOut.te2-self.gaussian(dataOut.heightList[:dataOut.NSHTS],params[0],params[1],params[2]))
4341 x0_value = numpy.array([2000,250,20])
4342 popt = least_squares(func,x0=x0_value,verbose=0)
4343 A = popt.x[0]; B = popt.x[1]; C = popt.x[2]
4344 te2_smooth = self.gaussian(dataOut.heightList[:dataOut.NSHTS], A, B, C)
4345 #'''
4202
4346
4347 '''
4348 import matplotlib.pyplot as plt
4349 plt.clf()
4350 plt.plot(te2_smooth,dataOut.heightList[:dataOut.NSHTS],'*-',label = 'My Te')
4351 plt.plot(dataOut.te2,dataOut.heightList[:dataOut.NSHTS],'*-',label = 'Te')
4352 #plt.plot(signal.medfilt(dataOut.te2),dataOut.heightList[:dataOut.NSHTS],'*-',label = 'Te')
4353 plt.title("{}".format(datetime.datetime.fromtimestamp(dataOut.utctime)))
4354 plt.xlim(-50,3000)
4355 plt.grid()
4356 plt.legend()
4357 plt.savefig("/home/roberto/Pictures/Density_Comparison/Te/{}.png".format(dataOut.utctime))
4358 '''
4203 #print("**** ACF2 WRAPPER ***** ",fitacf_acf2.acf2.__doc__ )
4359 #print("**** ACF2 WRAPPER ***** ",fitacf_acf2.acf2.__doc__ )
4204
4360
4205 for i in range(dataOut.NSHTS):
4361 for i in range(dataOut.NSHTS):
4206 if dataOut.info2[i]==1:
4362 if dataOut.info2[i]==1:
4207 angle=dataOut.thb[i]*0.01745
4363 angle=dataOut.thb[i]*0.01745
4208 nue=nui[0]=nui[1]=nui[2]=0.0#nui[3]=0.0
4364 nue=nui[0]=nui[1]=nui[2]=0.0#nui[3]=0.0
4209 wion[0]=16
4365 wion[0]=16 #O
4210 wion[1]=1
4366 wion[1]=1 #H
4211 wion[2]=4
4367 wion[2]=4
4212 tion[0]=tion[1]=tion[2]=dataOut.ti2[i]
4368 tion[0]=tion[1]=tion[2]=dataOut.ti2[i]
4213 fion[0]=1.0-dataOut.phy2[i]
4369 fion[0]=1.0-dataOut.phy2[i] #1
4214 fion[1]=dataOut.phy2[i]
4370 fion[1]=dataOut.phy2[i] #0
4215 fion[2]=0.0
4371 fion[2]=0.0 #0
4216 for j in range(dataOut.DPL):
4372 for j in range(dataOut.DPL):
4217 tau=dataOut.alag[j]*1.0e-3
4373 tau=dataOut.alag[j]*1.0e-3
4218
4374
4219 with suppress_stdout_stderr():
4375 with suppress_stdout_stderr():
4220 y[j]=fitacf_acf2.acf2(wl,tau,dataOut.te2[i],tion,fion,nue,nui,wion,angle,dataOut.ph2[i],dataOut.bfm[i],y[j],three)
4376 y[j]=fitacf_acf2.acf2(wl,tau,dataOut.te2[i],tion,fion,nue,nui,wion,angle,dataOut.ph2[i],dataOut.bfm[i],y[j],three)
4377 #y[j]=fitacf_acf2.acf2(wl,tau,my_te2[i],tion,fion,nue,nui,wion,angle,dataOut.ph2[i],dataOut.bfm[i],y[j],three)
4221
4378
4222 #if dataOut.ut_Faraday>11.0 and dataOut.range1[i]>150.0 and dataOut.range1[i]<400.0:
4379 #if dataOut.ut_Faraday>11.0 and dataOut.range1[i]>150.0 and dataOut.range1[i]<400.0:
4223
4380
4224 if dataOut.ut_Faraday>11.0 and dataOut.range1[i]>150.0 and dataOut.range1[i]<300.0:
4381 if dataOut.ut_Faraday>11.0 and dataOut.range1[i]>150.0 and dataOut.range1[i]<300.0:
4382 #if dataOut.ut_Faraday>11.0 and dataOut.range1[i]>150.0 and dataOut.range1[i]<400.0:
4225 tau=0.0
4383 tau=0.0
4226 with suppress_stdout_stderr():
4384 with suppress_stdout_stderr():
4227 bline=fitacf_acf2.acf2(wl,tau,tion,tion,fion,nue,nui,wion,angle,dataOut.ph2[i],dataOut.bfm[i],bline,three)
4385 bline=fitacf_acf2.acf2(wl,tau,tion,tion,fion,nue,nui,wion,angle,dataOut.ph2[i],dataOut.bfm[i],bline,three)
4228 cf=min(1.2,max(1.0,bline/y[0]))
4386 cf=min(1.2,max(1.0,bline/y[0])) #FACTOR DE EFICIENCIA
4387 #cf = bline/y[0]
4388 #cf=min(2.,max(1.0,bline/y[0]))
4229 my_aux[i] = cf
4389 my_aux[i] = cf
4230 #dataOut.ph2[i]=cf*dataOut.ph2[i] #Now we adjust the curve "cf" into a Gaussian,
4390 acf_Temps[i] = y[0]
4391 acf_no_Temps[i] = bline
4392 #dataOut.ph2[i]=cf*dataOut.ph2[i] #Instead we adjust the curve "cf" into a Gaussian,
4231 #dataOut.sdp2[i]=cf*dataOut.sdp2[i] #in order to get smoother values of density
4393 #dataOut.sdp2[i]=cf*dataOut.sdp2[i] #in order to get smoother values of density
4232 for j in range(1,dataOut.DPL):
4394 for j in range(1,dataOut.DPL):
4233 y[j]=(y[j]/y[0])*dataOut.DH+dataOut.range1[i]
4395 #y[j]=(y[j]/y[0])*dataOut.DH+dataOut.range1[i]
4396 y[j]=min(max((y[j]/y[0]),-1.0),1.0)*dataOut.DH+dataOut.range1[i]
4234 y[0]=dataOut.range1[i]+dataOut.DH
4397 y[0]=dataOut.range1[i]+dataOut.DH
4235
4398
4236
4399
4237 ratio = my_aux-1
4400 ratio = my_aux-1
4401 #ratio = dataOut.te2[:dataOut.NSHTS]/dataOut.ti2[:dataOut.NSHTS]
4238 def lsq_func(params):
4402 def lsq_func(params):
4239 return (ratio-self.gaussian(dataOut.heightList[:dataOut.NSHTS],params[0],params[1],params[2]))
4403 return (ratio-self.gaussian(dataOut.heightList[:dataOut.NSHTS],params[0],params[1],params[2]))
4240
4404
@@ -4244,30 +4408,40 class DenCorrection(NormalizeDPPowerRoberto_V2):
4244
4408
4245 A = popt.x[0]; B = popt.x[1]; C = popt.x[2]
4409 A = popt.x[0]; B = popt.x[1]; C = popt.x[2]
4246
4410
4247 aux = self.gaussian(dataOut.heightList[:dataOut.NSHTS], A, B, C) + 1
4411 aux = self.gaussian(dataOut.heightList[:dataOut.NSHTS], A, B, C) + 1 #ratio + 1
4248
4412
4249 '''
4413 '''
4250 import matplotlib.pyplot as plt
4414 import matplotlib.pyplot as plt
4251 plt.clf()
4415 plt.clf()
4252 plt.plot(aux,dataOut.heightList[:dataOut.NSHTS],label='Fitting')
4416 plt.plot(aux,dataOut.heightList[:dataOut.NSHTS],'*:',label='Fitting')
4253 plt.plot(my_aux,dataOut.heightList[:dataOut.NSHTS],label='Ratio')
4417 plt.plot(my_aux,dataOut.heightList[:dataOut.NSHTS],'*:',label='Ratio')
4418 #plt.plot(acf_Temps,dataOut.heightList[:dataOut.NSHTS],'b*:',label='Temps')
4419 #plt.plot(acf_no_Temps,dataOut.heightList[:dataOut.NSHTS],'k*:',label='No Temps')
4420 #plt.plot(dataOut.te2[:dataOut.NSHTS]/dataOut.ti2[:dataOut.NSHTS],dataOut.heightList[:dataOut.NSHTS],label='Ratio')
4254 #plt.ylim(180)
4421 #plt.ylim(180)
4255 plt.title("{}".format(datetime.fromtimestamp(dataOut.utctime)))
4422 plt.title("{}".format(datetime.datetime.fromtimestamp(dataOut.utctime)))
4256 plt.legend()
4423 plt.legend()
4257 plt.grid()
4424 plt.grid()
4425 #plt.xlim(.99,1.25)
4258 #plt.show()
4426 #plt.show()
4259 plt.savefig("/home/roberto/Pictures/Faraday_TeTi_Test/Ratio_at_400km/{}.png".format(dataOut.utctime))
4427 #plt.savefig("/home/roberto/Pictures/Density_Comparison/FactorEf_NoLimits/{}.png".format(dataOut.utctime))
4428 plt.savefig("/home/roberto/Pictures/Faraday/2022/08/Density_Comparison/FactorEf/{}.png".format(dataOut.utctime))
4260 #plt.savefig("/home/roberto/Pictures/Faraday_TeTi_Test/Ratio/{}.png".format(dataOut.utctime))
4429 #plt.savefig("/home/roberto/Pictures/Faraday_TeTi_Test/Ratio/{}.png".format(dataOut.utctime))
4261 '''
4430 '''
4262 #print("inside correction",dataOut.ph2)
4431 #print("inside correction",dataOut.ph2)
4432
4433 #print("heeere",aux)
4434 #exit(1)
4263 dataOut.ph2[:dataOut.NSHTS]*=aux
4435 dataOut.ph2[:dataOut.NSHTS]*=aux
4264 dataOut.sdp2[:dataOut.NSHTS]*=aux
4436 dataOut.sdp2[:dataOut.NSHTS]*=aux
4437 #dataOut.ph2[:26]*=aux[:26]
4438 #dataOut.sdp2[:26]*=aux[:26]
4265 #print(aux)
4439 #print(aux)
4266 #print("inside correction",dataOut.ph2)
4440 #print("inside correction",dataOut.ph2)
4267
4441
4268 def run(self,dataOut):
4442 def run(self,dataOut):
4269 #print("hour",gmtime(dataOut.utctime).tm_hour)
4443 #print("hour",gmtime(dataOut.utctime).tm_hour)
4270 if gmtime(dataOut.utctime).tm_hour < 23. and gmtime(dataOut.utctime).tm_hour >= 11.:
4444 if gmtime(dataOut.utctime).tm_hour < 24. and gmtime(dataOut.utctime).tm_hour >= 11.:
4271 #print("inside")
4445 #print("inside")
4272 self.TeTiEstimation(dataOut)
4446 self.TeTiEstimation(dataOut)
4273 dataOut.flagTeTiCorrection = True
4447 dataOut.flagTeTiCorrection = True
@@ -4277,6 +4451,9 class DenCorrection(NormalizeDPPowerRoberto_V2):
4277 return dataOut
4451 return dataOut
4278
4452
4279 class DataPlotCleaner(Operation):
4453 class DataPlotCleaner(Operation):
4454 '''
4455 Written by R. Flores
4456 '''
4280 def __init__(self, **kwargs):
4457 def __init__(self, **kwargs):
4281
4458
4282 Operation.__init__(self, **kwargs)
4459 Operation.__init__(self, **kwargs)
@@ -4341,6 +4518,9 class DataPlotCleaner(Operation):
4341
4518
4342
4519
4343 class DataSaveCleaner(Operation):
4520 class DataSaveCleaner(Operation):
4521 '''
4522 Written by R. Flores
4523 '''
4344 def __init__(self, **kwargs):
4524 def __init__(self, **kwargs):
4345
4525
4346 Operation.__init__(self, **kwargs)
4526 Operation.__init__(self, **kwargs)
@@ -4350,6 +4530,7 class DataSaveCleaner(Operation):
4350 #print(dataOut.heightList)
4530 #print(dataOut.heightList)
4351 #exit(1)
4531 #exit(1)
4352 dataOut.DensityFinal=numpy.zeros((1,dataOut.NDP))
4532 dataOut.DensityFinal=numpy.zeros((1,dataOut.NDP))
4533 dataOut.dphiFinal=numpy.zeros((1,dataOut.NDP))
4353 dataOut.EDensityFinal=numpy.zeros((1,dataOut.NDP))
4534 dataOut.EDensityFinal=numpy.zeros((1,dataOut.NDP))
4354 dataOut.ElecTempFinal=numpy.zeros((1,dataOut.NDP))
4535 dataOut.ElecTempFinal=numpy.zeros((1,dataOut.NDP))
4355 dataOut.EElecTempFinal=numpy.zeros((1,dataOut.NDP))
4536 dataOut.EElecTempFinal=numpy.zeros((1,dataOut.NDP))
@@ -4359,6 +4540,7 class DataSaveCleaner(Operation):
4359 dataOut.EPhyFinal=numpy.zeros((1,dataOut.NDP))
4540 dataOut.EPhyFinal=numpy.zeros((1,dataOut.NDP))
4360
4541
4361 dataOut.DensityFinal[0]=numpy.copy(dataOut.ph2)
4542 dataOut.DensityFinal[0]=numpy.copy(dataOut.ph2)
4543 dataOut.dphiFinal[0]=numpy.copy(dataOut.dphi)
4362 dataOut.EDensityFinal[0]=numpy.copy(dataOut.sdp2)
4544 dataOut.EDensityFinal[0]=numpy.copy(dataOut.sdp2)
4363 dataOut.ElecTempFinal[0,:dataOut.NSHTS]=numpy.copy(dataOut.te2)
4545 dataOut.ElecTempFinal[0,:dataOut.NSHTS]=numpy.copy(dataOut.te2)
4364 dataOut.EElecTempFinal[0,:dataOut.NSHTS]=numpy.copy(dataOut.ete2)
4546 dataOut.EElecTempFinal[0,:dataOut.NSHTS]=numpy.copy(dataOut.ete2)
@@ -4368,14 +4550,14 class DataSaveCleaner(Operation):
4368 dataOut.EPhyFinal[0,:dataOut.NSHTS]=numpy.copy(dataOut.ephy2)
4550 dataOut.EPhyFinal[0,:dataOut.NSHTS]=numpy.copy(dataOut.ephy2)
4369
4551
4370 missing=numpy.nan
4552 missing=numpy.nan
4371
4553 #print("den1: ",dataOut.DensityFinal)
4372 temp_min=100.0
4554 temp_min=100.0
4373 temp_max=3000.0#6000.0e
4555 temp_max=3000.0#6000.0e
4374 #print("Density: ",dataOut.DensityFinal[0])
4556 #print("Density: ",dataOut.DensityFinal[0])
4375 #print("Error: ",dataOut.EDensityFinal[0])
4557 #print("Error: ",dataOut.EDensityFinal[0])
4376 #print(100*dataOut.EDensityFinal[0]/dataOut.DensityFinal[0])
4558 #print(100*dataOut.EDensityFinal[0]/dataOut.DensityFinal[0])
4377 den_err_percent = 100*dataOut.EDensityFinal[0]/dataOut.DensityFinal[0]
4559 den_err_percent = 100*dataOut.EDensityFinal[0]/dataOut.DensityFinal[0]
4378 max_den_err_per = 30#30 #Densidades con error mayor al 30% se setean en NaN
4560 max_den_err_per = 35#30 #Densidades con error mayor al 30% se setean en NaN
4379 for i in range(dataOut.NSHTS):
4561 for i in range(dataOut.NSHTS):
4380
4562
4381 if den_err_percent[i] >= max_den_err_per:
4563 if den_err_percent[i] >= max_den_err_per:
@@ -4478,25 +4660,28 class DataSaveCleaner(Operation):
4478 dataOut.EDensityFinal[0,i1:] = dataOut.DensityFinal[0,i1:] = missing
4660 dataOut.EDensityFinal[0,i1:] = dataOut.DensityFinal[0,i1:] = missing
4479 '''
4661 '''
4480
4662
4481 #'''
4663 '''
4664 #print("den2: ",dataOut.DensityFinal)
4482 if gmtime(dataOut.utctime).tm_hour >= 23. or gmtime(dataOut.utctime).tm_hour < 5.: #18-00 LT
4665 if gmtime(dataOut.utctime).tm_hour >= 23. or gmtime(dataOut.utctime).tm_hour < 5.: #18-00 LT
4483 nanindex = numpy.argwhere(numpy.isnan(dataOut.DensityFinal[0,:33]))
4666 nanindex = numpy.argwhere(numpy.isnan(dataOut.DensityFinal[0,:33]))
4484 #print(nanindex)
4667 #print(nanindex)
4485 i1 = nanindex[-1][0]
4668 i1 = nanindex[-1][0]
4486 #print("i1",i1)
4669 #print("i1",i1)
4487 dataOut.EDensityFinal[0,:i1] = dataOut.DensityFinal[0,:i1] = missing
4670 dataOut.EDensityFinal[0,:i1] = dataOut.DensityFinal[0,:i1] = missing
4671 #print("den3: ",dataOut.DensityFinal)
4488 elif gmtime(dataOut.utctime).tm_hour >= 6. or gmtime(dataOut.utctime).tm_hour < 11.: #18-00 LT
4672 elif gmtime(dataOut.utctime).tm_hour >= 6. or gmtime(dataOut.utctime).tm_hour < 11.: #18-00 LT
4489 nanindex = numpy.argwhere(numpy.isnan(dataOut.DensityFinal[0,:20]))
4673 nanindex = numpy.argwhere(numpy.isnan(dataOut.DensityFinal[0,:20]))
4490 #print(nanindex)
4674 #print(nanindex)
4491 i1 = nanindex[-1][0]
4675 i1 = nanindex[-1][0]
4492 #print("i1",i1)
4676 #print("i1",i1)
4493 dataOut.EDensityFinal[0,:i1] = dataOut.DensityFinal[0,:i1] = missing
4677 dataOut.EDensityFinal[0,:i1] = dataOut.DensityFinal[0,:i1] = missing
4494 #'''
4678 '''
4495 #print("den_nans: ",dataOut.DensityFinal[0,12:50])
4679 #print("den_nans: ",dataOut.DensityFinal[0,12:50])
4496 if numpy.count_nonzero(~numpy.isnan(dataOut.DensityFinal[0,12:50]))<=5:
4680 if numpy.count_nonzero(~numpy.isnan(dataOut.DensityFinal[0,12:50]))<=5:
4497 dataOut.DensityFinal[0,:]=dataOut.EDensityFinal[0,:]=missing
4681 dataOut.DensityFinal[0,:]=dataOut.EDensityFinal[0,:]=missing
4498 #for i in range(dataOut.NSHTS,dataOut.NDP):
4682 #for i in range(dataOut.NSHTS,dataOut.NDP):
4499 #for i in range(40,dataOut.NDP):
4683 #for i in range(40,dataOut.NDP):
4684 #print("den2: ",dataOut.DensityFinal)
4500 dataOut.DensityFinal[0,dataOut.NSHTS:]=missing
4685 dataOut.DensityFinal[0,dataOut.NSHTS:]=missing
4501 dataOut.EDensityFinal[0,dataOut.NSHTS:]=missing
4686 dataOut.EDensityFinal[0,dataOut.NSHTS:]=missing
4502 dataOut.ElecTempFinal[0,dataOut.NSHTS:]=missing
4687 dataOut.ElecTempFinal[0,dataOut.NSHTS:]=missing
@@ -4510,6 +4695,7 class DataSaveCleaner(Operation):
4510 '''
4695 '''
4511 nanindex = numpy.argwhere(numpy.isnan(dataOut.DensityFinal))
4696 nanindex = numpy.argwhere(numpy.isnan(dataOut.DensityFinal))
4512 i1 = nanindex[-1][0]
4697 i1 = nanindex[-1][0]
4698 print("i1",i1)
4513 dataOut.DensityFinal[0,i1+1:]=missing
4699 dataOut.DensityFinal[0,i1+1:]=missing
4514 dataOut.EDensityFinal[0,i1+1:]=missing
4700 dataOut.EDensityFinal[0,i1+1:]=missing
4515 dataOut.ElecTempFinal[0,i1+1:]=missing
4701 dataOut.ElecTempFinal[0,i1+1:]=missing
@@ -4519,6 +4705,27 class DataSaveCleaner(Operation):
4519 dataOut.PhyFinal[0,i1+1:]=missing
4705 dataOut.PhyFinal[0,i1+1:]=missing
4520 dataOut.EPhyFinal[0,i1+1:]=missing
4706 dataOut.EPhyFinal[0,i1+1:]=missing
4521 '''
4707 '''
4708 #'''
4709 if gmtime(dataOut.utctime).tm_hour >= 12. and gmtime(dataOut.utctime).tm_hour < 22.: #07-17 LT
4710 dataOut.DensityFinal[0,:13]=missing
4711 dataOut.EDensityFinal[0,:13]=missing
4712 dataOut.ElecTempFinal[0,:13]=missing
4713 dataOut.EElecTempFinal[0,:13]=missing
4714 dataOut.IonTempFinal[0,:13]=missing
4715 dataOut.EIonTempFinal[0,:13]=missing
4716 dataOut.PhyFinal[0,:13]=missing
4717 dataOut.EPhyFinal[0,:13]=missing
4718 #'''
4719 else:
4720 dataOut.DensityFinal[0,:dataOut.min_id_eej+1]=missing
4721 dataOut.EDensityFinal[0,:dataOut.min_id_eej+1]=missing
4722 dataOut.ElecTempFinal[0,:dataOut.min_id_eej+1]=missing
4723 dataOut.EElecTempFinal[0,:dataOut.min_id_eej+1]=missing
4724 dataOut.IonTempFinal[0,:dataOut.min_id_eej+1]=missing
4725 dataOut.EIonTempFinal[0,:dataOut.min_id_eej+1]=missing
4726 dataOut.PhyFinal[0,:dataOut.min_id_eej+1]=missing
4727 dataOut.EPhyFinal[0,:dataOut.min_id_eej+1]=missing
4728 '''
4522 if gmtime(dataOut.utctime).tm_hour >= 11. or gmtime(dataOut.utctime).tm_hour < 23.: #06-18 LT
4729 if gmtime(dataOut.utctime).tm_hour >= 11. or gmtime(dataOut.utctime).tm_hour < 23.: #06-18 LT
4523 dataOut.DensityFinal[0,:13]=missing
4730 dataOut.DensityFinal[0,:13]=missing
4524 dataOut.EDensityFinal[0,:13]=missing
4731 dataOut.EDensityFinal[0,:13]=missing
@@ -4538,6 +4745,7 class DataSaveCleaner(Operation):
4538 dataOut.EIonTempFinal[0,:12]=missing
4745 dataOut.EIonTempFinal[0,:12]=missing
4539 dataOut.PhyFinal[0,:12]=missing
4746 dataOut.PhyFinal[0,:12]=missing
4540 dataOut.EPhyFinal[0,:12]=missing
4747 dataOut.EPhyFinal[0,:12]=missing
4748 '''
4541 #print(dataOut.EDensityFinal)
4749 #print(dataOut.EDensityFinal)
4542 #exit(1)
4750 #exit(1)
4543 '''
4751 '''
@@ -4545,9 +4753,8 class DataSaveCleaner(Operation):
4545 print(dataOut.heightList)
4753 print(dataOut.heightList)
4546 exit(1)
4754 exit(1)
4547 '''
4755 '''
4548
4549 '''
4550 time_text = datetime.datetime.utcfromtimestamp(dataOut.utctime)
4756 time_text = datetime.datetime.utcfromtimestamp(dataOut.utctime)
4757 #'''
4551 #Agregar año y mes para no estar comentando esta parte!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4758 #Agregar año y mes para no estar comentando esta parte!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4552 #if (time_text.hour == 13 and time_text.minute == 20) or (time_text.hour == 0 and time_text.minute > 40) or (time_text.hour == 2 and time_text.minute == 8) or (time_text.hour == 2 and time_text.minute == 19): #Year: 2022, DOY:101
4759 #if (time_text.hour == 13 and time_text.minute == 20) or (time_text.hour == 0 and time_text.minute > 40) or (time_text.hour == 2 and time_text.minute == 8) or (time_text.hour == 2 and time_text.minute == 19): #Year: 2022, DOY:101
4553 #if (time_text.hour == 17 and time_text.minute == 5) or (time_text.hour == 8 and time_text.minute >= 22) or (time_text.hour == 9) or (time_text.hour == 11 and time_text.minute == 2): #Year: 2022, DOY:102
4760 #if (time_text.hour == 17 and time_text.minute == 5) or (time_text.hour == 8 and time_text.minute >= 22) or (time_text.hour == 9) or (time_text.hour == 11 and time_text.minute == 2): #Year: 2022, DOY:102
@@ -4558,7 +4765,14 class DataSaveCleaner(Operation):
4558 #if (time_text.hour == 6 and time_text.minute>=34) or (time_text.hour >= 7 and time_text.hour<11) or (time_text.hour == 11 and time_text.minute<47): #Year: 2022, DOY:94
4765 #if (time_text.hour == 6 and time_text.minute>=34) or (time_text.hour >= 7 and time_text.hour<11) or (time_text.hour == 11 and time_text.minute<47): #Year: 2022, DOY:94
4559 #if (time_text.hour == 7 and time_text.minute>=18) or (time_text.hour >= 8 and time_text.hour<11) or (time_text.hour == 11 and time_text.minute==2) or (time_text.hour == 1) or (time_text.hour == 2 and time_text.minute<=20) or (time_text.hour >= 3 and time_text.hour<=4): #Year: 2022, DOY:92
4766 #if (time_text.hour == 7 and time_text.minute>=18) or (time_text.hour >= 8 and time_text.hour<11) or (time_text.hour == 11 and time_text.minute==2) or (time_text.hour == 1) or (time_text.hour == 2 and time_text.minute<=20) or (time_text.hour >= 3 and time_text.hour<=4): #Year: 2022, DOY:92
4560 #if (time_text.hour < 11 and time_text.hour >=5 ) or (time_text.hour == 11 and time_text.minute<=34): #Year: 2022, DOY:93
4767 #if (time_text.hour < 11 and time_text.hour >=5 ) or (time_text.hour == 11 and time_text.minute<=34): #Year: 2022, DOY:93
4561 if (time_text.hour == 7 and time_text.minute>=18) or (time_text.hour >= 8 and time_text.hour <=11 ): #Year: 2022, DOY:91
4768 #if (time_text.hour == 7 and time_text.minute>=18) or (time_text.hour >= 8 and time_text.hour <=11 ): #Year: 2022, DOY:91
4769
4770 #print(time_text.hour,time_text.minute)
4771 #if (time_text.hour == 16 and time_text.minute==48) or (time_text.hour == 19 and time_text.minute ==49 ) or (time_text.hour >= 0 and time_text.hour < 5): #Year: 2022, DOY:241
4772 #if (time_text.hour == 5 and time_text.minute==21) or (time_text.hour == 19 and time_text.minute ==49 ) or (time_text.hour == 7 and time_text.minute==40) or (time_text.hour == 7 and time_text.minute==50) or (time_text.hour >= 8 and time_text.hour < 11) or (time_text.hour == 11 and time_text.minute==2) or (time_text.hour == 11 and time_text.minute==13): #Year: 2022, DOY:242
4773 #if (time_text.hour >= 8 and time_text.hour < 11) or (time_text.hour == 11 and time_text.minute==2) or (time_text.hour == 11 and time_text.minute==13) or (time_text.hour == 11 and time_text.minute==24): #Year: 2022, DOY:243
4774 #if (time_text.hour >= 9 and time_text.hour < 11) or (time_text.hour == 8 and time_text.minute==12) or (time_text.hour == 8 and time_text.minute==22) or (time_text.hour == 8 and time_text.minute==33) or (time_text.hour == 8 and time_text.minute==44) or (time_text.hour == 8 and time_text.minute==54) or (time_text.hour == 11 and time_text.minute==2) or (time_text.hour == 11 and time_text.minute==13): #Year: 2022, DOY:245
4775 if (time_text.hour >= 8 and time_text.hour < 11) or (time_text.hour == 1) or (time_text.hour == 0 and time_text.minute==25) or (time_text.hour == 0 and time_text.minute==36) or (time_text.hour == 0 and time_text.minute==47) or (time_text.hour == 0 and time_text.minute==57) or (time_text.hour == 2 and time_text.minute==1) or (time_text.hour == 11 and time_text.minute==2) or (time_text.hour == 11 and time_text.minute==13) or (time_text.hour == 11 and time_text.minute==24) or (time_text.hour == 7 and time_text.minute==40) or (time_text.hour == 7 and time_text.minute==50) or (time_text.hour == 3 and time_text.minute==5): #Year: 2022, DOY:244
4562
4776
4563 dataOut.DensityFinal[0,:]=missing
4777 dataOut.DensityFinal[0,:]=missing
4564 dataOut.EDensityFinal[0,:]=missing
4778 dataOut.EDensityFinal[0,:]=missing
@@ -4569,30 +4783,139 class DataSaveCleaner(Operation):
4569 dataOut.PhyFinal[0,:]=missing
4783 dataOut.PhyFinal[0,:]=missing
4570 dataOut.EPhyFinal[0,:]=missing
4784 dataOut.EPhyFinal[0,:]=missing
4571
4785
4572 dataOut.flagNoData = True
4786 dataOut.flagNoData = True #Remueve todo el perfil
4573 '''
4787 #'''
4574 '''
4788 '''
4575 if (time_text.hour >= 7 and time_text.hour < 9): #Year: 2022, DOY:102
4789 #if (time_text.hour >= 7 and time_text.hour < 9): #Year: 2022, DOY:102
4576 dataOut.DensityFinal[0,33:]=missing
4790 if (time_text.hour == 20 and time_text.minute == 8): #Year: 2022, DOY:243
4577 dataOut.EDensityFinal[0,33:]=missing
4791 id_aux = 35
4578 dataOut.ElecTempFinal[0,33:]=missing
4792 dataOut.DensityFinal[0,id_aux:]=missing
4579 dataOut.EElecTempFinal[0,33:]=missing
4793 dataOut.EDensityFinal[0,id_aux:]=missing
4580 dataOut.IonTempFinal[0,33:]=missing
4794 dataOut.ElecTempFinal[0,id_aux:]=missing
4581 dataOut.EIonTempFinal[0,33:]=missing
4795 dataOut.EElecTempFinal[0,id_aux:]=missing
4582 dataOut.PhyFinal[0,33:]=missing
4796 dataOut.IonTempFinal[0,id_aux:]=missing
4583 dataOut.EPhyFinal[0,33:]=missing
4797 dataOut.EIonTempFinal[0,id_aux:]=missing
4584
4798 dataOut.PhyFinal[0,id_aux:]=missing
4585 #dataOut.flagNoData = True
4799 dataOut.EPhyFinal[0,id_aux:]=missing
4800 if (time_text.hour == 20 and time_text.minute == 19): #Year: 2022, DOY:243
4801 id_aux = 33
4802 dataOut.DensityFinal[0,id_aux:]=missing
4803 dataOut.EDensityFinal[0,id_aux:]=missing
4804 dataOut.ElecTempFinal[0,id_aux:]=missing
4805 dataOut.EElecTempFinal[0,id_aux:]=missing
4806 dataOut.IonTempFinal[0,id_aux:]=missing
4807 dataOut.EIonTempFinal[0,id_aux:]=missing
4808 dataOut.PhyFinal[0,id_aux:]=missing
4809 dataOut.EPhyFinal[0,id_aux:]=missing
4810 if (time_text.hour == 20 and time_text.minute == 29) or (time_text.hour == 20 and time_text.minute == 44): #Year: 2022, DOY:243
4811 id_aux = 31
4812 dataOut.DensityFinal[0,id_aux:]=missing
4813 dataOut.EDensityFinal[0,id_aux:]=missing
4814 dataOut.ElecTempFinal[0,id_aux:]=missing
4815 dataOut.EElecTempFinal[0,id_aux:]=missing
4816 dataOut.IonTempFinal[0,id_aux:]=missing
4817 dataOut.EIonTempFinal[0,id_aux:]=missing
4818 dataOut.PhyFinal[0,id_aux:]=missing
4819 dataOut.EPhyFinal[0,id_aux:]=missing
4586 '''
4820 '''
4587
4821 #'''
4822 if (time_text.hour == 2 and time_text.minute == 23): #Year: 2022, DOY:244
4823 id_aux = 12
4824 dataOut.DensityFinal[0,:id_aux]=missing
4825 dataOut.EDensityFinal[0,:id_aux]=missing
4826 dataOut.ElecTempFinal[0,:id_aux]=missing
4827 dataOut.EElecTempFinal[0,:id_aux]=missing
4828 dataOut.IonTempFinal[0,:id_aux]=missing
4829 dataOut.EIonTempFinal[0,:id_aux]=missing
4830 dataOut.PhyFinal[0,:id_aux]=missing
4831 dataOut.EPhyFinal[0,:id_aux]=missing
4832 if (time_text.hour == 20 and time_text.minute == 52): #Year: 2022, DOY:244
4833 id_aux = 25
4834 dataOut.DensityFinal[0,id_aux:]=missing
4835 dataOut.EDensityFinal[0,id_aux:]=missing
4836 dataOut.ElecTempFinal[0,id_aux:]=missing
4837 dataOut.EElecTempFinal[0,id_aux:]=missing
4838 dataOut.IonTempFinal[0,id_aux:]=missing
4839 dataOut.EIonTempFinal[0,id_aux:]=missing
4840 dataOut.PhyFinal[0,id_aux:]=missing
4841 dataOut.EPhyFinal[0,id_aux:]=missing
4842 if (time_text.hour == 21 and time_text.minute == 3): #Year: 2022, DOY:244
4843 id_aux = 38
4844 dataOut.DensityFinal[0,id_aux:]=missing
4845 dataOut.EDensityFinal[0,id_aux:]=missing
4846 dataOut.ElecTempFinal[0,id_aux:]=missing
4847 dataOut.EElecTempFinal[0,id_aux:]=missing
4848 dataOut.IonTempFinal[0,id_aux:]=missing
4849 dataOut.EIonTempFinal[0,id_aux:]=missing
4850 dataOut.PhyFinal[0,id_aux:]=missing
4851 dataOut.EPhyFinal[0,id_aux:]=missing
4852 if (time_text.hour == 21 and time_text.minute == 13): #Year: 2022, DOY:244
4853 id_aux = 26
4854 dataOut.DensityFinal[0,id_aux:]=missing
4855 dataOut.EDensityFinal[0,id_aux:]=missing
4856 dataOut.ElecTempFinal[0,id_aux:]=missing
4857 dataOut.EElecTempFinal[0,id_aux:]=missing
4858 dataOut.IonTempFinal[0,id_aux:]=missing
4859 dataOut.EIonTempFinal[0,id_aux:]=missing
4860 dataOut.PhyFinal[0,id_aux:]=missing
4861 dataOut.EPhyFinal[0,id_aux:]=missing
4862 if (time_text.hour == 21 and time_text.minute == 24): #Year: 2022, DOY:244
4863 id_aux = 37
4864 dataOut.DensityFinal[0,id_aux:]=missing
4865 dataOut.EDensityFinal[0,id_aux:]=missing
4866 dataOut.ElecTempFinal[0,id_aux:]=missing
4867 dataOut.EElecTempFinal[0,id_aux:]=missing
4868 dataOut.IonTempFinal[0,id_aux:]=missing
4869 dataOut.EIonTempFinal[0,id_aux:]=missing
4870 dataOut.PhyFinal[0,id_aux:]=missing
4871 dataOut.EPhyFinal[0,id_aux:]=missing
4872 if (time_text.hour == 21 and time_text.minute == 35): #Year: 2022, DOY:244
4873 id_aux = 37
4874 dataOut.DensityFinal[0,id_aux:]=missing
4875 dataOut.EDensityFinal[0,id_aux:]=missing
4876 dataOut.ElecTempFinal[0,id_aux:]=missing
4877 dataOut.EElecTempFinal[0,id_aux:]=missing
4878 dataOut.IonTempFinal[0,id_aux:]=missing
4879 dataOut.EIonTempFinal[0,id_aux:]=missing
4880 dataOut.PhyFinal[0,id_aux:]=missing
4881 dataOut.EPhyFinal[0,id_aux:]=missing
4882 if (time_text.hour == 21 and time_text.minute == 45): #Year: 2022, DOY:244
4883 id_aux = 36
4884 dataOut.DensityFinal[0,id_aux:]=missing
4885 dataOut.EDensityFinal[0,id_aux:]=missing
4886 dataOut.ElecTempFinal[0,id_aux:]=missing
4887 dataOut.EElecTempFinal[0,id_aux:]=missing
4888 dataOut.IonTempFinal[0,id_aux:]=missing
4889 dataOut.EIonTempFinal[0,id_aux:]=missing
4890 dataOut.PhyFinal[0,id_aux:]=missing
4891 dataOut.EPhyFinal[0,id_aux:]=missing
4892 if (time_text.hour == 3 and time_text.minute == 59): #Year: 2022, DOY:244
4893 id_aux = 36
4894 dataOut.DensityFinal[0,id_aux:]=missing
4895 dataOut.EDensityFinal[0,id_aux:]=missing
4896 dataOut.ElecTempFinal[0,id_aux:]=missing
4897 dataOut.EElecTempFinal[0,id_aux:]=missing
4898 dataOut.IonTempFinal[0,id_aux:]=missing
4899 dataOut.EIonTempFinal[0,id_aux:]=missing
4900 dataOut.PhyFinal[0,id_aux:]=missing
4901 dataOut.EPhyFinal[0,id_aux:]=missing
4902 #'''
4588 #print("den_final",dataOut.DensityFinal)
4903 #print("den_final",dataOut.DensityFinal)
4589
4904
4590 dataOut.flagNoData = numpy.all(numpy.isnan(dataOut.DensityFinal)) #Si todos los valores son NaN no se prosigue
4905 dataOut.flagNoData = numpy.all(numpy.isnan(dataOut.DensityFinal)) #Si todos los valores son NaN no se prosigue
4906
4907 ####dataOut.flagNoData = False #Solo para ploteo
4908
4909 dataOut.DensityFinal *= 1.e6 #Convert units to m^⁻3
4910 dataOut.EDensityFinal *= 1.e6 #Convert units to m^⁻3
4591 #print(dataOut.flagNoData)
4911 #print(dataOut.flagNoData)
4592 return dataOut
4912 return dataOut
4593
4913
4594
4914
4595 class DataSaveCleanerHP(Operation):
4915 class DataSaveCleanerHP(Operation):
4916 '''
4917 Written by R. Flores
4918 '''
4596 def __init__(self, **kwargs):
4919 def __init__(self, **kwargs):
4597
4920
4598 Operation.__init__(self, **kwargs)
4921 Operation.__init__(self, **kwargs)
@@ -4770,6 +5093,9 class DataSaveCleanerHP(Operation):
4770
5093
4771
5094
4772 class ACFs(Operation):
5095 class ACFs(Operation):
5096 '''
5097 Written by R. Flores
5098 '''
4773 def __init__(self, **kwargs):
5099 def __init__(self, **kwargs):
4774
5100
4775 Operation.__init__(self, **kwargs)
5101 Operation.__init__(self, **kwargs)
@@ -5110,7 +5436,7 class CohInt(Operation):
5110 if not self.isConfig:
5436 if not self.isConfig:
5111 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
5437 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
5112 self.isConfig = True
5438 self.isConfig = True
5113 print("inside")
5439 #print("inside")
5114 if dataOut.flagDataAsBlock:
5440 if dataOut.flagDataAsBlock:
5115 """
5441 """
5116 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
5442 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
@@ -5141,6 +5467,9 class CohInt(Operation):
5141 return dataOut
5467 return dataOut
5142
5468
5143 class TimesCode(Operation):
5469 class TimesCode(Operation):
5470 '''
5471 Written by R. Flores
5472 '''
5144 """
5473 """
5145
5474
5146 """
5475 """
@@ -5149,8 +5478,6 class TimesCode(Operation):
5149
5478
5150 Operation.__init__(self, **kwargs)
5479 Operation.__init__(self, **kwargs)
5151
5480
5152
5153
5154 def run(self,dataOut,code):
5481 def run(self,dataOut,code):
5155
5482
5156 #code = numpy.repeat(code, repeats=osamp, axis=1)
5483 #code = numpy.repeat(code, repeats=osamp, axis=1)
@@ -5214,44 +5541,11 class TimesCode(Operation):
5214
5541
5215 return dataOut
5542 return dataOut
5216
5543
5217 '''
5218 class Spectrogram(Operation):
5219 """
5220
5221 """
5222
5223 def __init__(self, **kwargs):
5224
5225 Operation.__init__(self, **kwargs)
5226
5227
5228
5229 def run(self,dataOut):
5230
5231 import scipy
5232
5233
5234
5235 fs = 3200*1e-6
5236 fs = fs/64
5237 fs = 1/fs
5238
5239 nperseg=64
5240 noverlap=48
5241
5242 f, t, Sxx = signal.spectrogram(x, fs, return_onesided=False, nperseg=nperseg, noverlap=noverlap, mode='complex')
5243
5244
5245 for ich in range(dataOut.nChannels):
5246 for ihe in range(nheicode):
5247
5248
5249 return dataOut
5250 '''
5251
5252
5544
5253 class RemoveDcHae(Operation):
5545 class RemoveDcHae(Operation):
5254
5546 '''
5547 Written by R. Flores
5548 '''
5255 def __init__(self, **kwargs):
5549 def __init__(self, **kwargs):
5256
5550
5257 Operation.__init__(self, **kwargs)
5551 Operation.__init__(self, **kwargs)
@@ -5523,6 +5817,7 class Decoder(Operation):
5523 profilesList = range(self.__nProfiles)
5817 profilesList = range(self.__nProfiles)
5524 #print(numpy.shape(self.datadecTime))
5818 #print(numpy.shape(self.datadecTime))
5525 #print(numpy.shape(data))
5819 #print(numpy.shape(data))
5820 #print(profilesList)
5526 for i in range(self.__nChannels):
5821 for i in range(self.__nChannels):
5527 for j in profilesList:
5822 for j in profilesList:
5528 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
5823 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
@@ -6693,9 +6988,9 class CrossProdHybrid(CrossProdDP):
6693 self.lag_products_LP_median_estimates_aux=0
6988 self.lag_products_LP_median_estimates_aux=0
6694
6989
6695
6990
6696 #for i in range(dataOut.NLAG):
6991 for i in range(dataOut.NLAG):
6697 my_list = ([0,1,2,3,4,5,6,7]) #hasta 7 funciona, en 6 ya no
6992 #my_list = ([0,1,2,3,4,5,6,7]) #hasta 7 funciona, en 6 ya no
6698 for i in my_list:
6993 #for i in my_list:
6699 for j in range(dataOut.NRANGE):
6994 for j in range(dataOut.NRANGE):
6700 for l in range(4): #four outputs
6995 for l in range(4): #four outputs
6701 for k in range(dataOut.NAVG):
6996 for k in range(dataOut.NAVG):
@@ -6709,6 +7004,7 class CrossProdHybrid(CrossProdDP):
6709 self.lagp2[i,j,:]=sorted(self.lagp2[i,j,:], key=lambda x: x.real) #sorted(self.lagp2[i,j,:].real)
7004 self.lagp2[i,j,:]=sorted(self.lagp2[i,j,:], key=lambda x: x.real) #sorted(self.lagp2[i,j,:].real)
6710 if l==3:
7005 if l==3:
6711 self.lagp3[i,j,:]=sorted(self.lagp3[i,j,:], key=lambda x: x.real) #sorted(self.lagp3[i,j,:].real)
7006 self.lagp3[i,j,:]=sorted(self.lagp3[i,j,:], key=lambda x: x.real) #sorted(self.lagp3[i,j,:].real)
7007 '''
6712 x = 2
7008 x = 2
6713 if k>=x and k<dataOut.NAVG-dataOut.nkill/2:
7009 if k>=x and k<dataOut.NAVG-dataOut.nkill/2:
6714 if l==0:
7010 if l==0:
@@ -6719,7 +7015,8 class CrossProdHybrid(CrossProdDP):
6719 self.output[i,j,l]=self.output[i,j,l]+((float(dataOut.NAVG)/(float)(dataOut.NAVG-x-dataOut.nkill/2))*self.lagp2[i,j,k])
7015 self.output[i,j,l]=self.output[i,j,l]+((float(dataOut.NAVG)/(float)(dataOut.NAVG-x-dataOut.nkill/2))*self.lagp2[i,j,k])
6720 if l==3:
7016 if l==3:
6721 self.output[i,j,l]=self.output[i,j,l]+((float(dataOut.NAVG)/(float)(dataOut.NAVG-x-dataOut.nkill/2))*self.lagp3[i,j,k])
7017 self.output[i,j,l]=self.output[i,j,l]+((float(dataOut.NAVG)/(float)(dataOut.NAVG-x-dataOut.nkill/2))*self.lagp3[i,j,k])
6722 '''
7018 '''
7019 #'''
6723 if k>=dataOut.nkill/2 and k<dataOut.NAVG-dataOut.nkill/2:
7020 if k>=dataOut.nkill/2 and k<dataOut.NAVG-dataOut.nkill/2:
6724 if l==0:
7021 if l==0:
6725 self.output[i,j,l]=self.output[i,j,l]+((float(dataOut.NAVG)/(float)(dataOut.NAVG-dataOut.nkill))*self.lagp0[i,j,k])
7022 self.output[i,j,l]=self.output[i,j,l]+((float(dataOut.NAVG)/(float)(dataOut.NAVG-dataOut.nkill))*self.lagp0[i,j,k])
@@ -6729,7 +7026,7 class CrossProdHybrid(CrossProdDP):
6729 self.output[i,j,l]=self.output[i,j,l]+((float(dataOut.NAVG)/(float)(dataOut.NAVG-dataOut.nkill))*self.lagp2[i,j,k])
7026 self.output[i,j,l]=self.output[i,j,l]+((float(dataOut.NAVG)/(float)(dataOut.NAVG-dataOut.nkill))*self.lagp2[i,j,k])
6730 if l==3:
7027 if l==3:
6731 self.output[i,j,l]=self.output[i,j,l]+((float(dataOut.NAVG)/(float)(dataOut.NAVG-dataOut.nkill))*self.lagp3[i,j,k])
7028 self.output[i,j,l]=self.output[i,j,l]+((float(dataOut.NAVG)/(float)(dataOut.NAVG-dataOut.nkill))*self.lagp3[i,j,k])
6732 '''
7029 #'''
6733
7030
6734
7031
6735 dataOut.output_LP=self.output
7032 dataOut.output_LP=self.output
@@ -6796,10 +7093,10 class CrossProdHybrid(CrossProdDP):
6796 #x02[i]=10.0*numpy.log10(x02[i])
7093 #x02[i]=10.0*numpy.log10(x02[i])
6797 return x00,x01,x02,x03
7094 return x00,x01,x02,x03
6798
7095
6799 def run(self, dataOut, NLAG=None, NRANGE=None, NCAL=None, DPL=None,
7096 def run(self, dataOut, NLAG=16, NRANGE=200, NCAL=0, DPL=11,
6800 NDN=None, NDT=None, NDP=None, NSCAN=None,
7097 NDN=0, NDT=67, NDP=67, NSCAN=128,
6801 lagind=None, lagfirst=None,
7098 lagind=(0,1,2,3,4,5,6,7,0,3,4,5,6,8,9,10), lagfirst=(1,1,1,1,1,1,1,1,0,0,0,0,0,1,1,1),
6802 NAVG=None, nkill=None):
7099 NAVG=16, nkill=6):
6803 #print(dataOut.data[1,:12,200:200+15])
7100 #print(dataOut.data[1,:12,200:200+15])
6804 #exit(1)
7101 #exit(1)
6805 dataOut.NLAG=NLAG
7102 dataOut.NLAG=NLAG
@@ -6958,6 +7255,7 class IntegrationHP(IntegrationDP):
6958 #print(dataOut.kabxys_integrated[8][53,6,0]+dataOut.kabxys_integrated[11][53,6,0])
7255 #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])
7256 #print(dataOut.kabxys_integrated[8][53,9,0]+dataOut.kabxys_integrated[11][53,9,0])
6960 #exit(1)
7257 #exit(1)
7258 print(dataOut.flagNoData)
6961 return dataOut
7259 return dataOut
6962
7260
6963 class SumFlipsHP(SumFlips):
7261 class SumFlipsHP(SumFlips):
@@ -6981,7 +7279,7 class SumFlipsHP(SumFlips):
6981 def rint2HP(self,dataOut):
7279 def rint2HP(self,dataOut):
6982
7280
6983 dataOut.rnint2=numpy.zeros(dataOut.DPL,'float32')
7281 dataOut.rnint2=numpy.zeros(dataOut.DPL,'float32')
6984
7282 #print(dataOut.nint,dataOut.NAVG)
6985 for l in range(dataOut.DPL):
7283 for l in range(dataOut.DPL):
6986 if(l==0 or (l>=3 and l <=6)):
7284 if(l==0 or (l>=3 and l <=6)):
6987 dataOut.rnint2[l]=0.5/float(dataOut.nint*dataOut.NAVG*16.0)
7285 dataOut.rnint2[l]=0.5/float(dataOut.nint*dataOut.NAVG*16.0)
@@ -7009,6 +7307,10 class SumFlipsHP(SumFlips):
7009 print((dataOut.kabxys_integrated[10][hei,lag,0]-dataOut.kabxys_integrated[9][hei,lag,0]))
7307 print((dataOut.kabxys_integrated[10][hei,lag,0]-dataOut.kabxys_integrated[9][hei,lag,0]))
7010 exit(1)
7308 exit(1)
7011 '''
7309 '''
7310 #print(dataOut.rnint2)
7311 #print(numpy.sum(dataOut.kabxys_integrated[4][:,1,0]+dataOut.kabxys_integrated[5][:,1,0]))
7312 #print(dataOut.nis)
7313 #exit(1)
7012 return dataOut
7314 return dataOut
7013
7315
7014
7316
General Comments 0
You need to be logged in to leave comments. Login now