##// 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 261 def getFmax(self):
262 262 PRF = 1. / (self.ippSeconds * self.nCohInt)
263
263 #print("ippsec",self.ippSeconds)
264 264 fmax = PRF
265 265 return fmax
266 266
Show file before | Show file after | Hide comments
@@ -256,6 +256,7 class Plot(Operation):
256 256 self.__throttle_plot = apply_throttle(self.throttle)
257 257 code = self.attr_data if self.attr_data else self.CODE
258 258 self.data = PlotterData(self.CODE, self.exp_code, self.localtime)
259 #self.EEJtype = kwargs.get('EEJtype', 2)
259 260
260 261 if self.server:
261 262 if not self.server.startswith('tcp://'):
Show file before | Show file after | Hide comments
@@ -75,7 +75,7 class SnrPlot(RTIPlot):
75 75 def update(self, dataOut):
76 76
77 77 data = {
78 'snr': 10*numpy.log10(dataOut.data_snr)
78 'snr': 10*numpy.log10(dataOut.data_snr)
79 79 }
80 80
81 81 return data, {}
@@ -91,7 +91,120 class DopplerPlot(RTIPlot):
91 91 def update(self, dataOut):
92 92
93 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 210 return data, {}
@@ -107,7 +220,7 class PowerPlot(RTIPlot):
107 220 def update(self, dataOut):
108 221
109 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 226 return data, {}
@@ -208,7 +321,7 class GenericRTIPlot(Plot):
208 321 meta = {}
209 322
210 323 return data, meta
211
324
212 325 def plot(self):
213 326 # self.data.normalize_heights()
214 327 self.x = self.data.times
Show file before | Show file after | Hide comments
@@ -39,9 +39,14 class SpectraPlot(Plot):
39 39
40 40 data = {}
41 41 meta = {}
42
42 43 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
44 #print("Spc: ",spc[0])
45 #exit(1)
43 46 data['spc'] = spc
44 47 data['rti'] = dataOut.getPower()
48 #print(data['rti'][0])
49 #exit(1)
45 50 #print("NormFactor: ",dataOut.normFactor)
46 51 #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
47 52 if hasattr(dataOut, 'LagPlot'): #Double Pulse
@@ -163,11 +168,22 class SpectraObliquePlot(Plot):
163 168 data['rti'] = dataOut.getPower()
164 169 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
165 170 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
166
167 data['shift1'] = dataOut.Oblique_params[0][1]
168 data['shift2'] = dataOut.Oblique_params[0][4]
169 data['shift1_error'] = dataOut.Oblique_param_errors[0][1]
170 data['shift2_error'] = dataOut.Oblique_param_errors[0][4]
171 '''
172 data['shift1'] = dataOut.Oblique_params[0,-2,:]
173 data['shift2'] = dataOut.Oblique_params[0,-1,:]
174 data['shift1_error'] = dataOut.Oblique_param_errors[0,-2,:]
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 188 return data, meta
173 189
@@ -187,15 +203,19 class SpectraObliquePlot(Plot):
187 203
188 204 y = self.data.yrange
189 205 self.y = y
190 z = self.data['spc']
206
207 data = self.data[-1]
208 z = data['spc']
191 209
192 210 for n, ax in enumerate(self.axes):
193 211 noise = self.data['noise'][n][-1]
194 shift1 = self.data['shift1']
195 shift2 = self.data['shift2']
196 err1 = self.data['shift1_error']
197 err2 = self.data['shift2_error']
212 shift1 = data['shift1']
213 #print(shift1)
214 shift2 = data['shift2']
215 err1 = data['shift1_error']
216 err2 = data['shift2_error']
198 217 if ax.firsttime:
218
199 219 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
200 220 self.xmin = self.xmin if self.xmin else -self.xmax
201 221 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
@@ -212,17 +232,20 class SpectraObliquePlot(Plot):
212 232 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
213 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)
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)
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)
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 239 else:
240 #print("else plotter1: ", self.ploterr1,shift1)
218 241 self.ploterr1.remove()
219 242 self.ploterr2.remove()
220 243 ax.plt.set_array(z[n].T.ravel())
221 244 if self.showprofile:
222 245 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
223 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)
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)
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)
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 250 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
228 251
@@ -671,6 +694,7 class RTIPlot(Plot):
671 694 data = {}
672 695 meta = {}
673 696 data['rti'] = dataOut.getPower()
697 #print(numpy.shape(data['rti']))
674 698
675 699 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
676 700
@@ -692,6 +716,7 class RTIPlot(Plot):
692 716 for n, ax in enumerate(self.axes):
693 717 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
694 718 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
719
695 720 if ax.firsttime:
696 721 ax.plt = ax.pcolormesh(x, y, z[n].T,
697 722 vmin=self.zmin,
Show file before | Show file after | Hide comments
@@ -213,7 +213,7 class DenRTIPlot(RTIPlot):
213 213 data = {}
214 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 218 return data, meta
219 219
Show file before | Show file after | Hide comments
@@ -70,7 +70,7 class ProcessingUnit(object):
70 70 def call(self, **kwargs):
71 71 '''
72 72 '''
73 #print("call")
73
74 74 try:
75 75 if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
76 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 125 runNextUnit = self.dataOut.isReady()
126 126 except:
127 127 runNextUnit = self.dataOut.isReady()
128 #exit(1)
128 129 #if not self.dataOut.isReady():
129 130 #return 'Error' if self.dataOut.error else input()
130 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 2 import math
3 3 from scipy import optimize, interpolate, signal, stats, ndimage
4 4 import scipy
5 from scipy.optimize import least_squares
5 6 import re
6 7 import datetime
7 8 import copy
@@ -91,6 +92,7 class ParametersProc(ProcessingUnit):
91 92 self.dataOut.timeInterval1 = self.dataIn.timeInterval
92 93 self.dataOut.heightList = self.dataIn.heightList
93 94 self.dataOut.frequency = self.dataIn.frequency
95 #self.dataOut.runNextUnit = self.dataIn.runNextUnit
94 96 #self.dataOut.noise = self.dataIn.noise
95 97
96 98 def run(self):
@@ -156,6 +158,7 class ParametersProc(ProcessingUnit):
156 158 if hasattr(self.dataIn, 'COFA'): #COFA
157 159 self.dataOut.COFA = self.dataIn.COFA
158 160
161 #self.dataOut.runNextUnit = self.dataIn.runNextUnit
159 162
160 163
161 164 #---------------------- Correlation Data ---------------------------
@@ -668,7 +671,9 class GaussianFit(Operation):
668 671 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.)
669 672
670 673 class Oblique_Gauss_Fit(Operation):
671
674 '''
675 Written by R. Flores
676 '''
672 677 def __init__(self):
673 678 Operation.__init__(self)
674 679
@@ -1060,6 +1065,17 class Oblique_Gauss_Fit(Operation):
1060 1065 val = a1 * numpy.exp(-z1**2/2) + a2 * numpy.exp(-y2**2/2)/(1-k2*z2) + d
1061 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 1079 def double_gaussian_double_skew(self,x, a1, b1, c1, k1, a2, b2, c2, k2, d):
1064 1080
1065 1081 z1 = (x-b1)/c1
@@ -1136,7 +1152,7 class Oblique_Gauss_Fit(Operation):
1136 1152 #return A1f, B1f, C1f, A2f, B2f, C2f, K2f, Df, doppler
1137 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 1157 from scipy.optimize import least_squares
1142 1158
@@ -1146,6 +1162,7 class Oblique_Gauss_Fit(Operation):
1146 1162 from scipy.signal import medfilt
1147 1163 Nincoh = 20
1148 1164 Nincoh = 80
1165 Nincoh = Nincoh
1149 1166 spcm = medfilt(spc,11)/numpy.sqrt(Nincoh)
1150 1167
1151 1168 # define a least squares function to optimize
@@ -1156,11 +1173,23 class Oblique_Gauss_Fit(Operation):
1156 1173 # bounds=([0,-460,0,0,-400,120,0],[numpy.inf,-340,50,numpy.inf,0,250,numpy.inf])
1157 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 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 1178 #print(bounds)
1161 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 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 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 1194 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1166 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 1222 #exit(1)
1194 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 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 1286 from scipy.optimize import least_squares
@@ -1264,7 +1352,124 class Oblique_Gauss_Fit(Operation):
1264 1352
1265 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 1474 pwcode = 1
1270 1475
@@ -1293,12 +1498,54 class Oblique_Gauss_Fit(Operation):
1293 1498 elif mode == 9:
1294 1499 dataOut.Oblique_params = numpy.ones((1,11,dataOut.nHeights))*numpy.NAN
1295 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 1509 dataOut.VelRange = x
1298 1510
1299 l1=range(22,36)
1511
1512
1513 #l1=range(22,36) #+62
1300 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 1550 if mode == 4:
1304 1551 '''
@@ -1321,9 +1568,10 class Oblique_Gauss_Fit(Operation):
1321 1568
1322 1569 for hei in itertools.chain(l1, l2):
1323 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 1574 try:
1326
1327 1575 spc = dataOut.data_spc[0,:,hei]
1328 1576
1329 1577 if mode == 6: #Skew Weighted Bounded
@@ -1340,14 +1588,39 class Oblique_Gauss_Fit(Operation):
1340 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 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)
1344 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,10,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei]))
1345 #print(hei)
1346 #print(dataOut.Oblique_params[0,10,hei])
1347 #print(dataOut.dplr_2_u[0,0,hei])
1348 #print("outside",dataOut.Oblique_param_errors[0,:,hei])
1349 #print("SUCCESSSSSSS")
1350 #exit(1)
1591 #if numpy.max(spc) <= 0:
1592 if x[numpy.argmax(spc)] <= 0:
1593 #print("EEJ")
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)
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:
1596 # dataOut.Oblique_params[0,:,hei] *= numpy.NAN
1597 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,10,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei]))
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 1625 else:
1353 1626 spc_fit, A1, B1, C1, D1 = self.Gauss_fit_2(spc,x,'first')
@@ -1393,7 +1666,34 class Oblique_Gauss_Fit(Operation):
1393 1666 ###dataOut.Oblique_params[0,:,hei] = dataOut.Oblique_params[0,:,hei]*numpy.NAN
1394 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 1698 return dataOut
1399 1699
@@ -3104,7 +3404,7 class SpectralFitting(Operation):
3104 3404 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
3105 3405 return chisq
3106 3406
3107 class WindProfiler(Operation):
3407 class WindProfiler_V0(Operation):
3108 3408
3109 3409 __isConfig = False
3110 3410
@@ -3656,7 +3956,8 class WindProfiler(Operation):
3656 3956 def run(self, dataOut, technique, nHours=1, hmin=70, hmax=110, **kwargs):
3657 3957
3658 3958 param = dataOut.data_param
3659 if dataOut.abscissaList != None:
3959 #if dataOut.abscissaList != None:
3960 if numpy.any(dataOut.abscissaList):
3660 3961 absc = dataOut.abscissaList[:-1]
3661 3962 # noise = dataOut.noise
3662 3963 heightList = dataOut.heightList
@@ -3821,11 +4122,60 class WindProfiler(Operation):
3821 4122
3822 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 4141 def __init__(self):
3827 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 4179 def __correctValues(self, heiRang, phi, velRadial, SNR):
3830 4180 listPhi = phi.tolist()
3831 4181 maxid = listPhi.index(max(listPhi))
@@ -3845,23 +4195,13 class EWDriftsEstimation(Operation):
3845 4195 for i in rango:
3846 4196 x = heiRang*math.cos(phi[i])
3847 4197 y1 = velRadial[i,:]
3848 vali= (numpy.isfinite(y1)==True).nonzero()
3849 y1=y1[vali]
3850 x = x[vali]
3851 f1 = interpolate.interp1d(x,y1,kind = 'cubic',bounds_error=False)
4198 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
3852 4199
3853 #heiRang1 = x*math.cos(phi[maxid])
3854 4200 x1 = heiRang1
3855 4201 y11 = f1(x1)
3856 4202
3857 4203 y2 = SNR[i,:]
3858 #print 'snr ', y2
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)
4204 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
3865 4205 y21 = f2(x1)
3866 4206
3867 4207 velRadial1[i,:] = y11
@@ -3869,36 +4209,716 class EWDriftsEstimation(Operation):
3869 4209
3870 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
3877 velRadial = dataOut.data_param[:,3,:]
3878 velRadialm = dataOut.data_param[:,2:4,:]*-1
4223 return velUVW
3879 4224
3880 rbufc=dataOut.data_paramC[:,:,0]
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
4225 # def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
3893 4226
3894 zenith = numpy.array(zenith)
3895 zenith -= zenithCorrection
3896 zenith *= numpy.pi/180
3897 if zenithCorrection != 0 :
3898 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
3899 else :
3900 heiRang1 = heiRang
3901 velRadial1 = velRadial
4227 def techniqueDBS(self, kwargs):
4228 """
4229 Function that implements Doppler Beam Swinging (DBS) technique.
4230
4231 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
4232 Direction correction (if necessary), Ranges and SNR
4233
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 4922 SNR1 = SNR
3903 4923
3904 4924 alp = zenith[0]
@@ -5542,6 +6562,9 class SMOperations():
5542 6562
5543 6563
5544 6564 class IGRFModel(Operation):
6565 '''
6566 Written by R. Flores
6567 '''
5545 6568 """Operation to calculate Geomagnetic parameters.
5546 6569
5547 6570 Parameters:
@@ -5597,10 +6620,13 class MergeProc(ProcessingUnit):
5597 6620 ProcessingUnit.__init__(self)
5598 6621
5599 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 6625 self.dataOut = getattr(self, self.inputs[0])
5602 6626 data_inputs = [getattr(self, attr) for attr in self.inputs]
5603 6627 #print(data_inputs)
6628 #print("Run: ",self.dataOut.runNextUnit)
6629 #exit(1)
5604 6630 #print(numpy.shape([getattr(data, attr_data) for data in data_inputs][1]))
5605 6631 #exit(1)
5606 6632 if mode==0:
@@ -5668,39 +6694,7 class MergeProc(ProcessingUnit):
5668 6694 #setattr(self.dataOut, 'data_acf', getattr(data_inputs[1], attr_data_5)) #LP
5669 6695 setattr(self.dataOut, 'data_acf', getattr(data_inputs[1], attr_data_5)) #LP
5670 6696 #print("Merge data_acf: ",self.dataOut.data_acf.shape)
5671 #exit(1)
5672 #print(self.dataOut.data_spc_LP.shape)
5673 #print("Exit")
5674 #exit(1)
5675 #setattr(self.dataOut, 'dataLag_cspc', [getattr(data, attr_data_2) for data in data_inputs][0])
5676 #setattr(self.dataOut, 'dataLag_cspc_LP', [getattr(data, attr_data_2) for data in data_inputs][1])
5677 #setattr(self.dataOut, 'nIncohInt', [getattr(data, attr_data_3) for data in data_inputs][0])
5678 #setattr(self.dataOut, 'nIncohInt_LP', [getattr(data, attr_data_3) for data in data_inputs][1])
5679 '''
5680 print(self.dataOut.dataLag_spc_LP.shape)
5681 print(self.dataOut.dataLag_cspc_LP.shape)
5682 exit(1)
5683 '''
5684 '''
5685 print(self.dataOut.dataLag_spc_LP[0,:,100])
5686 print(self.dataOut.dataLag_spc_LP[1,:,100])
5687 exit(1)
5688 '''
5689 #self.dataOut.dataLag_spc_LP = numpy.transpose(self.dataOut.dataLag_spc_LP[0],(2,0,1))
5690 #self.dataOut.dataLag_cspc_LP = numpy.transpose(self.dataOut.dataLag_cspc_LP,(3,1,2,0))
5691 '''
5692 print("Merge")
5693 print(numpy.shape(self.dataOut.dataLag_spc))
5694 print(numpy.shape(self.dataOut.dataLag_spc_LP))
5695 print(numpy.shape(self.dataOut.dataLag_cspc))
5696 print(numpy.shape(self.dataOut.dataLag_cspc_LP))
5697 exit(1)
5698 '''
5699 #print(numpy.sum(self.dataOut.dataLag_spc_LP[2,:,164])/128)
5700 #print(numpy.sum(self.dataOut.dataLag_cspc_LP[0,:,30,1])/128)
5701 #exit(1)
5702 #print(self.dataOut.NDP)
5703 #print(self.dataOut.nNoiseProfiles)
6697
5704 6698
5705 6699 #self.dataOut.nIncohInt_LP = 128
5706 6700 #self.dataOut.nProfiles_LP = 128#self.dataOut.nIncohInt_LP
@@ -5711,8 +6705,45 class MergeProc(ProcessingUnit):
5711 6705 #print("sahpi",self.dataOut.nIncohInt_LP)
5712 6706 #exit(1)
5713 6707 self.dataOut.NLAG = 16
6708 self.dataOut.NLAG = self.dataOut.data_acf.shape[1]
5714 6709 self.dataOut.NRANGE = self.dataOut.data_acf.shape[-1]
5715 6710
5716 6711 #print(numpy.shape(self.dataOut.data_spc))
5717 6712
5718 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 123 self.dataOut.flagShiftFFT = False
124 124
125 125 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False, runNextUnit = 0):
126
126
127 127 self.dataIn.runNextUnit = runNextUnit
128 128 if self.dataIn.type == "Spectra":
129
129 130 self.dataOut.copy(self.dataIn)
130 131 if shift_fft:
131 132 #desplaza a la derecha en el eje 2 determinadas posiciones
@@ -147,9 +148,15 class SpectraProc(ProcessingUnit):
147 148
148 149 if nProfiles == None:
149 150 nProfiles = nFFTPoints
150
151 #print(self.dataOut.ipp)
152 #exit(1)
151 153 if ippFactor == None:
152 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 161 self.dataOut.nFFTPoints = nFFTPoints
155 162
@@ -425,6 +432,46 class SpectraProc(ProcessingUnit):
425 432
426 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 475 class removeDC(Operation):
429 476
430 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 19 from schainpy.utils import log
20 20 from schainpy.model.data import _HS_algorithm
21 21
22 from schainpy.model.proc.jroproc_voltage import CleanCohEchoes
23
22 24 from time import time, mktime, strptime, gmtime, ctime
23 25
24 26
25 27 class SpectraLagProc(ProcessingUnit):
28 '''
29 Written by R. Flores
30 '''
26 31 def __init__(self):
27 32
28 33 ProcessingUnit.__init__(self)
@@ -308,7 +313,9 class SpectraLagProc(ProcessingUnit):
308 313 #print("after",self.dataOut.data_spc[0,:,20])
309 314
310 315 class removeDCLag(Operation):
311
316 '''
317 Written by R. Flores
318 '''
312 319 def remover(self,mode):
313 320 jspectra = self.dataOut.data_spc
314 321 jcspectra = self.dataOut.data_cspc
@@ -414,7 +421,9 class removeDCLag(Operation):
414 421
415 422
416 423 class removeDCLagFlip(Operation):
417
424 '''
425 Written by R. Flores
426 '''
418 427 #CHANGES MADE ONLY FOR MODE 2 AND NOT CONSIDERING CSPC
419 428
420 429 def remover(self,mode):
@@ -501,7 +510,7 class removeDCLagFlip(Operation):
501 510
502 511
503 512 def run(self, dataOut, mode=2):
504 #print("***********************************Remove DC***********************************")
513 print("***********************************Remove DC***********************************")
505 514 ##print(dataOut.FlipChannels)
506 515 #exit(1)
507 516 self.dataOut = dataOut
@@ -1098,7 +1107,9 class removeInterference(Operation):
1098 1107
1099 1108
1100 1109 class IntegrationFaradaySpectra(Operation):
1101
1110 '''
1111 Written by R. Flores
1112 '''
1102 1113 __profIndex = 0
1103 1114 __withOverapping = False
1104 1115
@@ -1441,7 +1452,9 class IntegrationFaradaySpectra(Operation):
1441 1452
1442 1453
1443 1454 class IntegrationFaradaySpectra2(Operation):
1444
1455 '''
1456 Written by R. Flores
1457 '''
1445 1458 __profIndex = 0
1446 1459 __withOverapping = False
1447 1460
@@ -2061,7 +2074,9 class IntegrationFaradaySpectra2(Operation):
2061 2074 return dataOut
2062 2075
2063 2076 class IntegrationFaradaySpectra3(Operation): #This class should manage data with no lags as well
2064
2077 '''
2078 Written by R. Flores
2079 '''
2065 2080 __profIndex = 0
2066 2081 __withOverapping = False
2067 2082
@@ -2815,7 +2830,9 class IntegrationFaradaySpectra3(Operation): #This class should manage data with
2815 2830 return dataOut
2816 2831
2817 2832 class IntegrationFaradaySpectraNoLags(Operation):
2818
2833 '''
2834 Written by R. Flores
2835 '''
2819 2836 __profIndex = 0
2820 2837 __withOverapping = False
2821 2838
@@ -3149,6 +3166,9 class IntegrationFaradaySpectraNoLags(Operation):
3149 3166 return dataOut
3150 3167
3151 3168 class HybridSelectSpectra(Operation):
3169 '''
3170 Written by R. Flores
3171 '''
3152 3172 """Operation to rearange and use selected channels of spectra data and pairs of cross-spectra data for Hybrid Experiment.
3153 3173
3154 3174 Parameters:
@@ -3213,7 +3233,9 class HybridSelectSpectra(Operation):
3213 3233
3214 3234
3215 3235 class IncohIntLag(Operation):
3216
3236 '''
3237 Written by R. Flores
3238 '''
3217 3239 __profIndex = 0
3218 3240 __withOverapping = False
3219 3241
@@ -3598,6 +3620,7 class IncohInt(Operation):
3598 3620 #exit(1)
3599 3621
3600 3622 if n == 1:
3623 dataOut.VelRange = dataOut.getVelRange(0)
3601 3624 return dataOut
3602 3625
3603 3626 dataOut.flagNoData = True
@@ -3619,6 +3642,9 class IncohInt(Operation):
3619 3642 dataOut.nIncohInt *= self.n
3620 3643 dataOut.utctime = avgdatatime
3621 3644 dataOut.flagNoData = False
3645
3646 dataOut.VelRange = dataOut.getVelRange(0)
3647
3622 3648 #print("Power",numpy.sum(dataOut.data_spc[0,:,20:30],axis=0))
3623 3649 #print("Power",numpy.sum(dataOut.data_spc[0,100:110,:],axis=1))
3624 3650 #exit(1)
@@ -3678,6 +3704,9 class IncohInt(Operation):
3678 3704 return dataOut
3679 3705
3680 3706 class SnrFaraday(Operation):
3707 '''
3708 Written by R. Flores
3709 '''
3681 3710 """Operation to use get SNR in Faraday processing.
3682 3711
3683 3712 Parameters:
@@ -3723,6 +3752,9 class SnrFaraday(Operation):
3723 3752 return dataOut
3724 3753
3725 3754 class SpectraDataToFaraday(Operation):
3755 '''
3756 Written by R. Flores
3757 '''
3726 3758 """Operation to use spectra data in Faraday processing.
3727 3759
3728 3760 Parameters:
@@ -3754,12 +3786,19 class SpectraDataToFaraday(Operation):
3754 3786 dataOut.year=dataOut.bd_time.tm_year+(dataOut.bd_time.tm_yday-1)/364.0
3755 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 3790 tmpx=numpy.zeros((dataOut.nHeights,dataOut.DPL,2),'float32')
3759 3791 tmpx_a2=numpy.zeros((dataOut.nHeights,dataOut.DPL,2),'float32')
3760 3792 tmpx_b2=numpy.zeros((dataOut.nHeights,dataOut.DPL,2),'float32')
3761 3793 tmpx_abr=numpy.zeros((dataOut.nHeights,dataOut.DPL,2),'float32')
3762 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 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 3804 dataOut.rnint2=numpy.zeros(dataOut.DPL,'float32')
@@ -3806,25 +3845,7 class SpectraDataToFaraday(Operation):
3806 3845 for l in range(dataOut.DPL):
3807 3846 dataOut.rnint2[l]=1.0/(dataOut.nIncohInt*dataOut.nProfiles)
3808 3847
3809
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
3848 def noise(self,dataOut):
3828 3849
3829 3850 dataOut.noise_lag = numpy.zeros((dataOut.nChannels,dataOut.DPL),'float32')
3830 3851 #print("Lags")
@@ -3878,6 +3899,114 class SpectraDataToFaraday(Operation):
3878 3899 dataOut.pan = dataOut.tnoise[0]
3879 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 4010 #dataOut.pan = numpy.mean(dataOut.pan)
3882 4011 #dataOut.pbn = numpy.mean(dataOut.pbn)
3883 4012 #print(dataOut.pan)
@@ -3888,12 +4017,18 class SpectraDataToFaraday(Operation):
3888 4017 #print("Noise dB: ",10*numpy.log10(dataOut.tnoise))
3889 4018 #exit(1)
3890 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 4022 print("done")
4023 #exit(1)
3892 4024 return dataOut
3893 4025
3894 4026
3895 4027
3896 4028 class SpectraDataToHybrid(SpectraDataToFaraday):
4029 '''
4030 Written by R. Flores
4031 '''
3897 4032 """Operation to use spectra data in Faraday processing.
3898 4033
3899 4034 Parameters:
@@ -4067,6 +4202,9 class SpectraDataToHybrid(SpectraDataToFaraday):
4067 4202 return dataOut
4068 4203
4069 4204 class SpectraDataToHybrid_V2(SpectraDataToFaraday):
4205 '''
4206 Written by R. Flores
4207 '''
4070 4208 """Operation to use spectra data in Faraday processing.
4071 4209
4072 4210 Parameters:
@@ -4092,7 +4230,7 class SpectraDataToHybrid_V2(SpectraDataToFaraday):
4092 4230 self.dataLag_cspc_LP=None
4093 4231 self.dataLag_dc_LP=None
4094 4232
4095 def noise(self,dataOut):
4233 def noise_v0(self,dataOut):
4096 4234
4097 4235 dataOut.data_spc = dataOut.dataLag_spc_LP.real
4098 4236 #print(dataOut.dataLag_spc.shape)
@@ -4115,9 +4253,104 class SpectraDataToHybrid_V2(SpectraDataToFaraday):
4115 4253 #print("pbn: ",dataOut.pbn)
4116 4254 #print(numpy.shape(dataOut.pnoise))
4117 4255 #exit(1)
4118 #print("pan: ",numpy.sum(dataOut.pan))
4256 #print("pan: ",dataOut.pan)
4257 #print("pbn: ",dataOut.pbn)
4119 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 4354 def ConvertDataLP(self,dataOut):
4122 4355
4123 4356 #print(dataOut.dataLag_spc[:,:,:,1]/dataOut.data_spc)
@@ -4161,6 +4394,7 class SpectraDataToHybrid_V2(SpectraDataToFaraday):
4161 4394 #dataOut.output_LP_integrated[:,:,3] *= float(dataOut.NSCAN/22)#(dataOut.nNoiseProfiles) #Corrects the zero padding
4162 4395
4163 4396 dataOut.nis=dataOut.NSCAN*dataOut.nIncohInt_LP*10
4397 dataOut.nis=dataOut.NSCAN*dataOut.nIncohInt_LP*dataOut.nProfiles_LP*10
4164 4398
4165 4399 self.ConvertData(dataOut)
4166 4400
@@ -4170,10 +4404,98 class SpectraDataToHybrid_V2(SpectraDataToFaraday):
4170 4404 dataOut.kabxys_integrated[10][:,(1,2,7,8,9,10),0] *= 2 #Corrects the zero padding
4171 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 4417 dataOut.NAVG=1#dataOut.rnint2[0] #CHECK THIS!
4176 4418 dataOut.nint=dataOut.nIncohInt
4177 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 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 29 self.flip = 1
30 30 self.setupReq = False
31 31
32 def run(self):
32 def run(self, runNextUnit = 0):
33 33
34 34 if self.dataIn.type == 'AMISR':
35 35 self.__updateObjFromAmisrInput()
36 36
37 37 if self.dataIn.type == 'Voltage':
38 38 self.dataOut.copy(self.dataIn)
39 self.dataOut.runNextUnit = runNextUnit
40 #print(self.dataOut.data.shape)
39 41
40 42 def __updateObjFromAmisrInput(self):
41 43
@@ -397,6 +399,9 class deFlip(Operation):
397 399 return dataOut
398 400
399 401 class deFlipHP(Operation):
402 '''
403 Written by R. Flores
404 '''
400 405 def __init__(self):
401 406
402 407 self.flip = 1
@@ -530,6 +535,9 class interpolateHeights(Operation):
530 535
531 536
532 537 class LagsReshape(Operation):
538 '''
539 Written by R. Flores
540 '''
533 541 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
534 542
535 543 Parameters:
@@ -621,6 +629,9 class LagsReshape(Operation):
621 629 return dataOut
622 630
623 631 class LagsReshape150(Operation):
632 '''
633 Written by R. Flores
634 '''
624 635 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
625 636
626 637 Parameters:
@@ -714,6 +725,9 class LagsReshape150(Operation):
714 725 return dataOut
715 726
716 727 class LagsReshapeHP(Operation):
728 '''
729 Written by R. Flores
730 '''
717 731 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
718 732
719 733 Parameters:
@@ -805,6 +819,9 class LagsReshapeHP(Operation):
805 819 return dataOut
806 820
807 821 class LagsReshapeDP(Operation):
822 '''
823 Written by R. Flores
824 '''
808 825 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
809 826
810 827 Parameters:
@@ -957,6 +974,9 class LagsReshapeDP(Operation):
957 974 return dataOut
958 975
959 976 class LagsReshapeDP_V2(Operation):
977 '''
978 Written by R. Flores
979 '''
960 980 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
961 981
962 982 Parameters:
@@ -1156,6 +1176,9 class LagsReshapeDP_V2(Operation):
1156 1176 return dataOut
1157 1177
1158 1178 class LagsReshapeLP(Operation):
1179 '''
1180 Written by R. Flores
1181 '''
1159 1182 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
1160 1183
1161 1184 Parameters:
@@ -1322,6 +1345,9 class LagsReshapeLP(Operation):
1322 1345 return dataOut
1323 1346
1324 1347 class LagsReshapeHP2(Operation):
1348 '''
1349 Written by R. Flores
1350 '''
1325 1351 """Operation to reshape input data into (Channels,Profiles(with same lag),Heights,Lags) and heights reconstruction.
1326 1352
1327 1353 Parameters:
@@ -1499,6 +1525,9 class LagsReshapeHP2(Operation):
1499 1525 return dataOut
1500 1526
1501 1527 class CrossProdDP(Operation):
1528 '''
1529 Written by R. Flores
1530 '''
1502 1531 """Operation to calculate cross products of the Double Pulse Experiment.
1503 1532
1504 1533 Parameters:
@@ -1983,6 +2012,9 class CrossProdDP(Operation):
1983 2012
1984 2013
1985 2014 class IntegrationDP(Operation):
2015 '''
2016 Written by R. Flores
2017 '''
1986 2018 """Operation to integrate the Double Pulse data.
1987 2019
1988 2020 Parameters:
@@ -2068,6 +2100,9 class IntegrationDP(Operation):
2068 2100
2069 2101
2070 2102 class SumFlips(Operation):
2103 '''
2104 Written by R. Flores
2105 '''
2071 2106 """Operation to sum the flip and unflip part of certain cross products of the Double Pulse.
2072 2107
2073 2108 Parameters:
@@ -2128,6 +2163,9 class SumFlips(Operation):
2128 2163
2129 2164
2130 2165 class FlagBadHeights(Operation):
2166 '''
2167 Written by R. Flores
2168 '''
2131 2169 """Operation to flag bad heights (bad data) of the Double Pulse.
2132 2170
2133 2171 Parameters:
@@ -2161,6 +2199,9 class FlagBadHeights(Operation):
2161 2199 return dataOut
2162 2200
2163 2201 class FlagBadHeightsSpectra(Operation):
2202 '''
2203 Written by R. Flores
2204 '''
2164 2205 """Operation to flag bad heights (bad data) of the Double Pulse.
2165 2206
2166 2207 Parameters:
@@ -2194,6 +2235,9 class FlagBadHeightsSpectra(Operation):
2194 2235 return dataOut
2195 2236
2196 2237 class CleanCohEchoes(Operation):
2238 '''
2239 Written by R. Flores
2240 '''
2197 2241 """Operation to clean coherent echoes.
2198 2242
2199 2243 Parameters:
@@ -2561,6 +2605,9 class CleanCohEchoes(Operation):
2561 2605 return dataOut
2562 2606
2563 2607 class CleanCohEchoesHP(Operation):
2608 '''
2609 Written by R. Flores
2610 '''
2564 2611 """Operation to clean coherent echoes.
2565 2612
2566 2613 Parameters:
@@ -2928,6 +2975,9 class CleanCohEchoesHP(Operation):
2928 2975 return dataOut
2929 2976
2930 2977 class NoisePower(Operation):
2978 '''
2979 Written by R. Flores
2980 '''
2931 2981 """Operation to get noise power from the integrated data of the Double Pulse.
2932 2982
2933 2983 Parameters:
@@ -3012,6 +3062,9 class NoisePower(Operation):
3012 3062
3013 3063
3014 3064 class DoublePulseACFs(Operation):
3065 '''
3066 Written by R. Flores
3067 '''
3015 3068 """Operation to get the ACFs of the Double Pulse.
3016 3069
3017 3070 Parameters:
@@ -3066,7 +3119,7 class DoublePulseACFs(Operation):
3066 3119 '''
3067 3120 #print("init 2.6",pa,dataOut.pan)
3068 3121 dataOut.p[i,j]=pa+pb-(dataOut.pan+dataOut.pbn)
3069
3122 #print(i,j,dataOut.p[i,j])
3070 3123 dataOut.sdp[i,j]=2*dataOut.rnint2[j]*((pa+pb)*(pa+pb))
3071 3124 ## ACF
3072 3125
@@ -3116,11 +3169,11 class DoublePulseACFs(Operation):
3116 3169 '''
3117 3170 #'''
3118 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 3173 dataOut.igcej[i,j]=1
3121 3174
3122 3175 elif ((pa/dataOut.pan-1.0)>2.25*(pb/dataOut.pbn-1.0)):
3123 #print("EJJ")
3176 #print(dataOut.heightList[i],"EJJ")
3124 3177 dataOut.igcej[i,j]=1
3125 3178 #'''
3126 3179 '''
@@ -3154,10 +3207,14 class DoublePulseACFs(Operation):
3154 3207 #print("p: ",dataOut.p[33,:])
3155 3208 #exit(1)
3156 3209 '''
3157
3210 #print(numpy.sum(dataOut.rhor))
3211 #exit(1)
3158 3212 return dataOut
3159 3213
3160 3214 class DoublePulseACFs_PerLag(Operation):
3215 '''
3216 Written by R. Flores
3217 '''
3161 3218 """Operation to get the ACFs of the Double Pulse.
3162 3219
3163 3220 Parameters:
@@ -3258,18 +3315,18 class DoublePulseACFs_PerLag(Operation):
3258 3315 '''
3259 3316 #'''
3260 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 3319 dataOut.igcej[i,j]=1
3263 3320
3264 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 3323 dataOut.igcej[i,j]=1
3267 3324 #'''
3268 3325 '''
3269 3326 if ((pa/dataOut.pan-1.0)>2.25*(pb/dataOut.pbn-1.0)):
3270 3327 #print("EJJ")
3271 3328 dataOut.igcej[i,j]=1
3272 '''
3329 #'''
3273 3330 '''
3274 3331 if i == 4:
3275 3332 exit(1)
@@ -3299,6 +3356,9 class DoublePulseACFs_PerLag(Operation):
3299 3356 return dataOut
3300 3357
3301 3358 class FaradayAngleAndDPPower(Operation):
3359 '''
3360 Written by R. Flores
3361 '''
3302 3362 """Operation to calculate Faraday angle and Double Pulse power.
3303 3363
3304 3364 Parameters:
@@ -3345,6 +3405,7 class FaradayAngleAndDPPower(Operation):
3345 3405 st=0.# // total signal
3346 3406 ibt=0# // bad lags
3347 3407 ns=0# // no. good lags
3408 #print(dataOut.heightList[j])
3348 3409 for l in range(dataOut.DPL):
3349 3410 #add in other lags if outside of e-jet contamination
3350 3411 if( (dataOut.igcej[j][l] == 0) and (dataOut.ibad[j][l] == 0) ):
@@ -3352,8 +3413,8 class FaradayAngleAndDPPower(Operation):
3352 3413 dataOut.ph2[j]+=dataOut.p[j][l]/dataOut.sdp[j][l]
3353 3414 dataOut.sdp2[j]=dataOut.sdp2[j]+1./dataOut.sdp[j][l]
3354 3415 ns+=1
3355
3356
3416 #if dataOut.igcej[j][l] != 0:
3417 #print(l)
3357 3418 pt+=dataOut.p[j][l]/dataOut.sdp[j][l]
3358 3419 st+=1./dataOut.sdp[j][l]
3359 3420 ibt|=dataOut.ibad[j][l];
@@ -3387,6 +3448,9 class FaradayAngleAndDPPower(Operation):
3387 3448
3388 3449
3389 3450 class ElectronDensityFaraday(Operation):
3451 '''
3452 Written by R. Flores
3453 '''
3390 3454 """Operation to calculate electron density from Faraday angle.
3391 3455
3392 3456 Parameters:
@@ -3433,6 +3497,7 class ElectronDensityFaraday(Operation):
3433 3497 ndphi=dataOut.NSHTS-4
3434 3498 #print(dataOut.phi)
3435 3499 #exit(1)
3500 #'''
3436 3501 if dataOut.flagSpreadF:
3437 3502 nanindex = numpy.argwhere(numpy.isnan(dataOut.phi))
3438 3503 i1 = nanindex[-1][0]
@@ -3443,6 +3508,7 class ElectronDensityFaraday(Operation):
3443 3508 else:
3444 3509 #dataOut.phi_uwrp = dataOut.phi.copy()
3445 3510 dataOut.phi[:]=numpy.unwrap(dataOut.phi[:]) #Better results
3511 #'''
3446 3512 #print(dataOut.phi)
3447 3513 #print(dataOut.ph2)
3448 3514 #exit(1)
@@ -3455,9 +3521,11 class ElectronDensityFaraday(Operation):
3455 3521 for i in range(2,dataOut.NSHTS-2):
3456 3522 fact=(-0.5/(dataOut.RATE*dataOut.DH))*dataOut.bki[i]
3457 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 3526 ##dataOut.dphi[i]=((((theta[i-2]-theta[i+2])+(8.0*(theta[i+1]-theta[i-1])))/thetai[i])).real/12.0
3460 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 3529 #dataOut.dphi_uc[i] = abs(dataOut.phi[i]*dataOut.bki[i]*(-0.5)/dataOut.DH)
3462 3530 dataOut.dphi[i]=abs(dataOut.dphi[i]*fact)
3463 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 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 3544 Parameters:
3474 3545 -----------
@@ -3602,6 +3673,9 class NormalizeDPPower(Operation):
3602 3673 return dataOut
3603 3674
3604 3675 class NormalizeDPPowerRoberto(Operation):
3676 '''
3677 Written by R. Flores
3678 '''
3605 3679 """Operation to normalize relative electron density from power with total electron density from Farday angle.
3606 3680
3607 3681 Parameters:
@@ -3738,6 +3812,9 class NormalizeDPPowerRoberto(Operation):
3738 3812 return dataOut
3739 3813
3740 3814 class NormalizeDPPowerRoberto_V2(Operation):
3815 '''
3816 Written by R. Flores
3817 '''
3741 3818 """Operation to normalize relative electron density from power with total electron density from Farday angle.
3742 3819
3743 3820 Parameters:
@@ -3799,7 +3876,7 class NormalizeDPPowerRoberto_V2(Operation):
3799 3876 self.aux=0
3800 3877
3801 3878 night_first=250.0
3802 night_first1= 400.0
3879 night_first1= 350.0
3803 3880 night_end= 450.0
3804 3881 day_first=220.0
3805 3882 day_end=400.0
@@ -3812,7 +3889,7 class NormalizeDPPowerRoberto_V2(Operation):
3812 3889 #print("EARLY")
3813 3890 i2=(night_end-dataOut.range1[0])/dataOut.DH
3814 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 3893 #print("NIGHT")
3817 3894 i2=(night_end-dataOut.range1[0])/dataOut.DH
3818 3895 i1=(night_first1 -dataOut.range1[0])/dataOut.DH
@@ -3847,15 +3924,17 class NormalizeDPPowerRoberto_V2(Operation):
3847 3924 #print("Flag: ",dataOut.flagTeTiCorrection)
3848 3925 #print(dataOut.dphi[i1::])
3849 3926 #print(dataOut.ph2[:])
3927
3850 3928 if dataOut.flagTeTiCorrection:
3851 3929 for i in range(dataOut.NSHTS):
3852 3930 dataOut.ph2[i]/=dataOut.cf
3853 3931 dataOut.sdp2[i]/=dataOut.cf
3932
3854 3933 #'''
3855 3934 if dataOut.flagSpreadF:
3856 3935 i2=int((620-dataOut.range1[0])/dataOut.DH)
3857 3936 nanindex = numpy.argwhere(numpy.isnan(dataOut.ph2))
3858 print(nanindex)
3937 print("nanindex",nanindex)
3859 3938 i1 = nanindex[-1][0] #VER CUANDO i1>i2
3860 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 3940 #print("i1, i2",i1,i2)
@@ -3887,11 +3966,24 class NormalizeDPPowerRoberto_V2(Operation):
3887 3966 except:
3888 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 3972 #if (time_text.hour == 5 and time_text.minute == 32): #Year: 2022, DOY:104
3892 3973 #if (time_text.hour == 0 and time_text.minute == 12): #Year: 2022, DOY:93
3893 #dataOut.cf = dataOut.cflast[0]
3894
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
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 3987 dataOut.cflast[0]=dataOut.cf
3896 3988 #print(dataOut.cf)
3897 3989
@@ -3953,6 +4045,9 class suppress_stdout_stderr(object):
3953 4045
3954 4046
3955 4047 class DPTemperaturesEstimation(Operation):
4048 '''
4049 Written by R. Flores
4050 '''
3956 4051 """Operation to estimate temperatures for Double Pulse data.
3957 4052
3958 4053 Parameters:
@@ -4172,7 +4267,9 class DPTemperaturesEstimation(Operation):
4172 4267
4173 4268
4174 4269 class DenCorrection(NormalizeDPPowerRoberto_V2):
4175
4270 '''
4271 Written by R. Flores
4272 '''
4176 4273 def __init__(self, **kwargs):
4177 4274
4178 4275 Operation.__init__(self, **kwargs)
@@ -4199,42 +4296,109 class DenCorrection(NormalizeDPPowerRoberto_V2):
4199 4296 bline=0.0
4200 4297 #bline=numpy.zeros(1,order='F',dtype='float32')
4201 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 4359 #print("**** ACF2 WRAPPER ***** ",fitacf_acf2.acf2.__doc__ )
4204 4360
4205 4361 for i in range(dataOut.NSHTS):
4206 4362 if dataOut.info2[i]==1:
4207 4363 angle=dataOut.thb[i]*0.01745
4208 4364 nue=nui[0]=nui[1]=nui[2]=0.0#nui[3]=0.0
4209 wion[0]=16
4210 wion[1]=1
4365 wion[0]=16 #O
4366 wion[1]=1 #H
4211 4367 wion[2]=4
4212 4368 tion[0]=tion[1]=tion[2]=dataOut.ti2[i]
4213 fion[0]=1.0-dataOut.phy2[i]
4214 fion[1]=dataOut.phy2[i]
4215 fion[2]=0.0
4369 fion[0]=1.0-dataOut.phy2[i] #1
4370 fion[1]=dataOut.phy2[i] #0
4371 fion[2]=0.0 #0
4216 4372 for j in range(dataOut.DPL):
4217 4373 tau=dataOut.alag[j]*1.0e-3
4218 4374
4219 4375 with suppress_stdout_stderr():
4220 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 4379 #if dataOut.ut_Faraday>11.0 and dataOut.range1[i]>150.0 and dataOut.range1[i]<400.0:
4223 4380
4224 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 4383 tau=0.0
4226 4384 with suppress_stdout_stderr():
4227 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 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 4393 #dataOut.sdp2[i]=cf*dataOut.sdp2[i] #in order to get smoother values of density
4232 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 4397 y[0]=dataOut.range1[i]+dataOut.DH
4235 4398
4236 4399
4237 4400 ratio = my_aux-1
4401 #ratio = dataOut.te2[:dataOut.NSHTS]/dataOut.ti2[:dataOut.NSHTS]
4238 4402 def lsq_func(params):
4239 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 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 4414 import matplotlib.pyplot as plt
4251 4415 plt.clf()
4252 plt.plot(aux,dataOut.heightList[:dataOut.NSHTS],label='Fitting')
4253 plt.plot(my_aux,dataOut.heightList[:dataOut.NSHTS],label='Ratio')
4416 plt.plot(aux,dataOut.heightList[:dataOut.NSHTS],'*:',label='Fitting')
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 4421 #plt.ylim(180)
4255 plt.title("{}".format(datetime.fromtimestamp(dataOut.utctime)))
4422 plt.title("{}".format(datetime.datetime.fromtimestamp(dataOut.utctime)))
4256 4423 plt.legend()
4257 4424 plt.grid()
4425 #plt.xlim(.99,1.25)
4258 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 4429 #plt.savefig("/home/roberto/Pictures/Faraday_TeTi_Test/Ratio/{}.png".format(dataOut.utctime))
4261 4430 '''
4262 4431 #print("inside correction",dataOut.ph2)
4432
4433 #print("heeere",aux)
4434 #exit(1)
4263 4435 dataOut.ph2[:dataOut.NSHTS]*=aux
4264 4436 dataOut.sdp2[:dataOut.NSHTS]*=aux
4437 #dataOut.ph2[:26]*=aux[:26]
4438 #dataOut.sdp2[:26]*=aux[:26]
4265 4439 #print(aux)
4266 4440 #print("inside correction",dataOut.ph2)
4267 4441
4268 4442 def run(self,dataOut):
4269 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 4445 #print("inside")
4272 4446 self.TeTiEstimation(dataOut)
4273 4447 dataOut.flagTeTiCorrection = True
@@ -4277,6 +4451,9 class DenCorrection(NormalizeDPPowerRoberto_V2):
4277 4451 return dataOut
4278 4452
4279 4453 class DataPlotCleaner(Operation):
4454 '''
4455 Written by R. Flores
4456 '''
4280 4457 def __init__(self, **kwargs):
4281 4458
4282 4459 Operation.__init__(self, **kwargs)
@@ -4341,6 +4518,9 class DataPlotCleaner(Operation):
4341 4518
4342 4519
4343 4520 class DataSaveCleaner(Operation):
4521 '''
4522 Written by R. Flores
4523 '''
4344 4524 def __init__(self, **kwargs):
4345 4525
4346 4526 Operation.__init__(self, **kwargs)
@@ -4350,6 +4530,7 class DataSaveCleaner(Operation):
4350 4530 #print(dataOut.heightList)
4351 4531 #exit(1)
4352 4532 dataOut.DensityFinal=numpy.zeros((1,dataOut.NDP))
4533 dataOut.dphiFinal=numpy.zeros((1,dataOut.NDP))
4353 4534 dataOut.EDensityFinal=numpy.zeros((1,dataOut.NDP))
4354 4535 dataOut.ElecTempFinal=numpy.zeros((1,dataOut.NDP))
4355 4536 dataOut.EElecTempFinal=numpy.zeros((1,dataOut.NDP))
@@ -4359,6 +4540,7 class DataSaveCleaner(Operation):
4359 4540 dataOut.EPhyFinal=numpy.zeros((1,dataOut.NDP))
4360 4541
4361 4542 dataOut.DensityFinal[0]=numpy.copy(dataOut.ph2)
4543 dataOut.dphiFinal[0]=numpy.copy(dataOut.dphi)
4362 4544 dataOut.EDensityFinal[0]=numpy.copy(dataOut.sdp2)
4363 4545 dataOut.ElecTempFinal[0,:dataOut.NSHTS]=numpy.copy(dataOut.te2)
4364 4546 dataOut.EElecTempFinal[0,:dataOut.NSHTS]=numpy.copy(dataOut.ete2)
@@ -4368,14 +4550,14 class DataSaveCleaner(Operation):
4368 4550 dataOut.EPhyFinal[0,:dataOut.NSHTS]=numpy.copy(dataOut.ephy2)
4369 4551
4370 4552 missing=numpy.nan
4371
4553 #print("den1: ",dataOut.DensityFinal)
4372 4554 temp_min=100.0
4373 4555 temp_max=3000.0#6000.0e
4374 4556 #print("Density: ",dataOut.DensityFinal[0])
4375 4557 #print("Error: ",dataOut.EDensityFinal[0])
4376 4558 #print(100*dataOut.EDensityFinal[0]/dataOut.DensityFinal[0])
4377 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 4561 for i in range(dataOut.NSHTS):
4380 4562
4381 4563 if den_err_percent[i] >= max_den_err_per:
@@ -4478,25 +4660,28 class DataSaveCleaner(Operation):
4478 4660 dataOut.EDensityFinal[0,i1:] = dataOut.DensityFinal[0,i1:] = missing
4479 4661 '''
4480 4662
4481 #'''
4663 '''
4664 #print("den2: ",dataOut.DensityFinal)
4482 4665 if gmtime(dataOut.utctime).tm_hour >= 23. or gmtime(dataOut.utctime).tm_hour < 5.: #18-00 LT
4483 4666 nanindex = numpy.argwhere(numpy.isnan(dataOut.DensityFinal[0,:33]))
4484 4667 #print(nanindex)
4485 4668 i1 = nanindex[-1][0]
4486 4669 #print("i1",i1)
4487 4670 dataOut.EDensityFinal[0,:i1] = dataOut.DensityFinal[0,:i1] = missing
4671 #print("den3: ",dataOut.DensityFinal)
4488 4672 elif gmtime(dataOut.utctime).tm_hour >= 6. or gmtime(dataOut.utctime).tm_hour < 11.: #18-00 LT
4489 4673 nanindex = numpy.argwhere(numpy.isnan(dataOut.DensityFinal[0,:20]))
4490 4674 #print(nanindex)
4491 4675 i1 = nanindex[-1][0]
4492 4676 #print("i1",i1)
4493 4677 dataOut.EDensityFinal[0,:i1] = dataOut.DensityFinal[0,:i1] = missing
4494 #'''
4678 '''
4495 4679 #print("den_nans: ",dataOut.DensityFinal[0,12:50])
4496 4680 if numpy.count_nonzero(~numpy.isnan(dataOut.DensityFinal[0,12:50]))<=5:
4497 4681 dataOut.DensityFinal[0,:]=dataOut.EDensityFinal[0,:]=missing
4498 4682 #for i in range(dataOut.NSHTS,dataOut.NDP):
4499 4683 #for i in range(40,dataOut.NDP):
4684 #print("den2: ",dataOut.DensityFinal)
4500 4685 dataOut.DensityFinal[0,dataOut.NSHTS:]=missing
4501 4686 dataOut.EDensityFinal[0,dataOut.NSHTS:]=missing
4502 4687 dataOut.ElecTempFinal[0,dataOut.NSHTS:]=missing
@@ -4510,6 +4695,7 class DataSaveCleaner(Operation):
4510 4695 '''
4511 4696 nanindex = numpy.argwhere(numpy.isnan(dataOut.DensityFinal))
4512 4697 i1 = nanindex[-1][0]
4698 print("i1",i1)
4513 4699 dataOut.DensityFinal[0,i1+1:]=missing
4514 4700 dataOut.EDensityFinal[0,i1+1:]=missing
4515 4701 dataOut.ElecTempFinal[0,i1+1:]=missing
@@ -4519,6 +4705,27 class DataSaveCleaner(Operation):
4519 4705 dataOut.PhyFinal[0,i1+1:]=missing
4520 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 4729 if gmtime(dataOut.utctime).tm_hour >= 11. or gmtime(dataOut.utctime).tm_hour < 23.: #06-18 LT
4523 4730 dataOut.DensityFinal[0,:13]=missing
4524 4731 dataOut.EDensityFinal[0,:13]=missing
@@ -4538,6 +4745,7 class DataSaveCleaner(Operation):
4538 4745 dataOut.EIonTempFinal[0,:12]=missing
4539 4746 dataOut.PhyFinal[0,:12]=missing
4540 4747 dataOut.EPhyFinal[0,:12]=missing
4748 '''
4541 4749 #print(dataOut.EDensityFinal)
4542 4750 #exit(1)
4543 4751 '''
@@ -4545,9 +4753,8 class DataSaveCleaner(Operation):
4545 4753 print(dataOut.heightList)
4546 4754 exit(1)
4547 4755 '''
4548
4549 '''
4550 4756 time_text = datetime.datetime.utcfromtimestamp(dataOut.utctime)
4757 #'''
4551 4758 #Agregar año y mes para no estar comentando esta parte!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4552 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 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 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 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 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 4777 dataOut.DensityFinal[0,:]=missing
4564 4778 dataOut.EDensityFinal[0,:]=missing
@@ -4569,30 +4783,139 class DataSaveCleaner(Operation):
4569 4783 dataOut.PhyFinal[0,:]=missing
4570 4784 dataOut.EPhyFinal[0,:]=missing
4571 4785
4572 dataOut.flagNoData = True
4573 '''
4786 dataOut.flagNoData = True #Remueve todo el perfil
4787 #'''
4574 4788 '''
4575 if (time_text.hour >= 7 and time_text.hour < 9): #Year: 2022, DOY:102
4576 dataOut.DensityFinal[0,33:]=missing
4577 dataOut.EDensityFinal[0,33:]=missing
4578 dataOut.ElecTempFinal[0,33:]=missing
4579 dataOut.EElecTempFinal[0,33:]=missing
4580 dataOut.IonTempFinal[0,33:]=missing
4581 dataOut.EIonTempFinal[0,33:]=missing
4582 dataOut.PhyFinal[0,33:]=missing
4583 dataOut.EPhyFinal[0,33:]=missing
4584
4585 #dataOut.flagNoData = True
4789 #if (time_text.hour >= 7 and time_text.hour < 9): #Year: 2022, DOY:102
4790 if (time_text.hour == 20 and time_text.minute == 8): #Year: 2022, DOY:243
4791 id_aux = 35
4792 dataOut.DensityFinal[0,id_aux:]=missing
4793 dataOut.EDensityFinal[0,id_aux:]=missing
4794 dataOut.ElecTempFinal[0,id_aux:]=missing
4795 dataOut.EElecTempFinal[0,id_aux:]=missing
4796 dataOut.IonTempFinal[0,id_aux:]=missing
4797 dataOut.EIonTempFinal[0,id_aux:]=missing
4798 dataOut.PhyFinal[0,id_aux:]=missing
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 4903 #print("den_final",dataOut.DensityFinal)
4589 4904
4590 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 4911 #print(dataOut.flagNoData)
4592 4912 return dataOut
4593 4913
4594 4914
4595 4915 class DataSaveCleanerHP(Operation):
4916 '''
4917 Written by R. Flores
4918 '''
4596 4919 def __init__(self, **kwargs):
4597 4920
4598 4921 Operation.__init__(self, **kwargs)
@@ -4770,6 +5093,9 class DataSaveCleanerHP(Operation):
4770 5093
4771 5094
4772 5095 class ACFs(Operation):
5096 '''
5097 Written by R. Flores
5098 '''
4773 5099 def __init__(self, **kwargs):
4774 5100
4775 5101 Operation.__init__(self, **kwargs)
@@ -5110,7 +5436,7 class CohInt(Operation):
5110 5436 if not self.isConfig:
5111 5437 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
5112 5438 self.isConfig = True
5113 print("inside")
5439 #print("inside")
5114 5440 if dataOut.flagDataAsBlock:
5115 5441 """
5116 5442 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
@@ -5141,6 +5467,9 class CohInt(Operation):
5141 5467 return dataOut
5142 5468
5143 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 5479 Operation.__init__(self, **kwargs)
5151 5480
5152
5153
5154 5481 def run(self,dataOut,code):
5155 5482
5156 5483 #code = numpy.repeat(code, repeats=osamp, axis=1)
@@ -5214,44 +5541,11 class TimesCode(Operation):
5214 5541
5215 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 5545 class RemoveDcHae(Operation):
5254
5546 '''
5547 Written by R. Flores
5548 '''
5255 5549 def __init__(self, **kwargs):
5256 5550
5257 5551 Operation.__init__(self, **kwargs)
@@ -5523,6 +5817,7 class Decoder(Operation):
5523 5817 profilesList = range(self.__nProfiles)
5524 5818 #print(numpy.shape(self.datadecTime))
5525 5819 #print(numpy.shape(data))
5820 #print(profilesList)
5526 5821 for i in range(self.__nChannels):
5527 5822 for j in profilesList:
5528 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 6988 self.lag_products_LP_median_estimates_aux=0
6694 6989
6695 6990
6696 #for i in range(dataOut.NLAG):
6697 my_list = ([0,1,2,3,4,5,6,7]) #hasta 7 funciona, en 6 ya no
6698 for i in my_list:
6991 for i in range(dataOut.NLAG):
6992 #my_list = ([0,1,2,3,4,5,6,7]) #hasta 7 funciona, en 6 ya no
6993 #for i in my_list:
6699 6994 for j in range(dataOut.NRANGE):
6700 6995 for l in range(4): #four outputs
6701 6996 for k in range(dataOut.NAVG):
@@ -6709,6 +7004,7 class CrossProdHybrid(CrossProdDP):
6709 7004 self.lagp2[i,j,:]=sorted(self.lagp2[i,j,:], key=lambda x: x.real) #sorted(self.lagp2[i,j,:].real)
6710 7005 if l==3:
6711 7006 self.lagp3[i,j,:]=sorted(self.lagp3[i,j,:], key=lambda x: x.real) #sorted(self.lagp3[i,j,:].real)
7007 '''
6712 7008 x = 2
6713 7009 if k>=x and k<dataOut.NAVG-dataOut.nkill/2:
6714 7010 if l==0:
@@ -6719,7 +7015,8 class CrossProdHybrid(CrossProdDP):
6719 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 7016 if l==3:
6721 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 7020 if k>=dataOut.nkill/2 and k<dataOut.NAVG-dataOut.nkill/2:
6724 7021 if l==0:
6725 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 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 7027 if l==3:
6731 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 7032 dataOut.output_LP=self.output
@@ -6796,10 +7093,10 class CrossProdHybrid(CrossProdDP):
6796 7093 #x02[i]=10.0*numpy.log10(x02[i])
6797 7094 return x00,x01,x02,x03
6798 7095
6799 def run(self, dataOut, NLAG=None, NRANGE=None, NCAL=None, DPL=None,
6800 NDN=None, NDT=None, NDP=None, NSCAN=None,
6801 lagind=None, lagfirst=None,
6802 NAVG=None, nkill=None):
7096 def run(self, dataOut, NLAG=16, NRANGE=200, NCAL=0, DPL=11,
7097 NDN=0, NDT=67, NDP=67, NSCAN=128,
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),
7099 NAVG=16, nkill=6):
6803 7100 #print(dataOut.data[1,:12,200:200+15])
6804 7101 #exit(1)
6805 7102 dataOut.NLAG=NLAG
@@ -6958,6 +7255,7 class IntegrationHP(IntegrationDP):
6958 7255 #print(dataOut.kabxys_integrated[8][53,6,0]+dataOut.kabxys_integrated[11][53,6,0])
6959 7256 #print(dataOut.kabxys_integrated[8][53,9,0]+dataOut.kabxys_integrated[11][53,9,0])
6960 7257 #exit(1)
7258 print(dataOut.flagNoData)
6961 7259 return dataOut
6962 7260
6963 7261 class SumFlipsHP(SumFlips):
@@ -6981,7 +7279,7 class SumFlipsHP(SumFlips):
6981 7279 def rint2HP(self,dataOut):
6982 7280
6983 7281 dataOut.rnint2=numpy.zeros(dataOut.DPL,'float32')
6984
7282 #print(dataOut.nint,dataOut.NAVG)
6985 7283 for l in range(dataOut.DPL):
6986 7284 if(l==0 or (l>=3 and l <=6)):
6987 7285 dataOut.rnint2[l]=0.5/float(dataOut.nint*dataOut.NAVG*16.0)
@@ -7009,6 +7307,10 class SumFlipsHP(SumFlips):
7009 7307 print((dataOut.kabxys_integrated[10][hei,lag,0]-dataOut.kabxys_integrated[9][hei,lag,0]))
7010 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 7314 return dataOut
7013 7315
7014 7316
General Comments 0
You need to be logged in to leave comments. Login now