##// END OF EJS Templates
schain
RSS
last update
rflores -
r1600:2c146cb976d3
parent child
Show More
Side by Side Unified
Context file:
r1600:2c146cb976d3 Loading diff...:
Show file before | Show file after | Hide comments
@@ -548,6 +548,9 class Spectra(JROData):
548
548
549 if self.flagDecodeData:
549 if self.flagDecodeData:
550 pwcode = numpy.sum(self.code[0]**2)
550 pwcode = numpy.sum(self.code[0]**2)
551 #pwcode = 64
552 #print("pwcode: ", pwcode)
553 #exit(1)
551 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
554 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
552 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
555 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
553
556
Show file before | Show file after | Hide comments
@@ -42,6 +42,9 class SpectraPlot(Plot):
42 meta = {}
42 meta = {}
43
43
44 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
44 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
45 #print("dataOut.normFactor: ", dataOut.normFactor)
46 #print("spc: ", dataOut.data_spc[0,0,0])
47 #spc = 10*numpy.log10(dataOut.data_spc)
45 #print("Spc: ",spc[0])
48 #print("Spc: ",spc[0])
46 #exit(1)
49 #exit(1)
47 data['spc'] = spc
50 data['spc'] = spc
@@ -716,14 +719,15 class RTIPlot(Plot):
716 self.x = self.data.times
719 self.x = self.data.times
717 self.y = self.data.yrange
720 self.y = self.data.yrange
718 self.z = self.data[self.CODE]
721 self.z = self.data[self.CODE]
719
722 #print("Inside RTI: ", self.z)
720 self.z = numpy.ma.masked_invalid(self.z)
723 self.z = numpy.ma.masked_invalid(self.z)
721
724
722 if self.decimation is None:
725 if self.decimation is None:
723 x, y, z = self.fill_gaps(self.x, self.y, self.z)
726 x, y, z = self.fill_gaps(self.x, self.y, self.z)
724 else:
727 else:
725 x, y, z = self.fill_gaps(*self.decimate())
728 x, y, z = self.fill_gaps(*self.decimate())
726
729 #print("self.z: ", self.z)
730 #exit(1)
727 '''
731 '''
728 if not isinstance(self.zmin, collections.abc.Sequence):
732 if not isinstance(self.zmin, collections.abc.Sequence):
729 if not self.zmin:
733 if not self.zmin:
Show file before | Show file after | Hide comments
@@ -673,9 +673,9 class EDensityPlot(Plot):
673 data['den_Faraday'] = dataOut.dphi[:dataOut.NSHTS]
673 data['den_Faraday'] = dataOut.dphi[:dataOut.NSHTS]
674 data['den_error'] = dataOut.sdp2[:dataOut.NSHTS]
674 data['den_error'] = dataOut.sdp2[:dataOut.NSHTS]
675 #data['err_Faraday'] = dataOut.sdn1[:dataOut.NSHTS]
675 #data['err_Faraday'] = dataOut.sdn1[:dataOut.NSHTS]
676 print(numpy.shape(data['den_power']))
676 #print(numpy.shape(data['den_power']))
677 print(numpy.shape(data['den_Faraday']))
677 #print(numpy.shape(data['den_Faraday']))
678 print(numpy.shape(data['den_error']))
678 #print(numpy.shape(data['den_error']))
679
679
680 data['NSHTS'] = dataOut.NSHTS
680 data['NSHTS'] = dataOut.NSHTS
681
681
Show file before | Show file after | Hide comments
@@ -49,6 +49,7 DEF_HEADER = {
49 MNEMONICS = {
49 MNEMONICS = {
50 10: 'jro',
50 10: 'jro',
51 11: 'jbr',
51 11: 'jbr',
52 14: 'jmp', #Added by R. Flores
52 840: 'jul',
53 840: 'jul',
53 13: 'jas',
54 13: 'jas',
54 1000: 'pbr',
55 1000: 'pbr',
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (1232 lines changed) Show them Hide them
@@ -1152,7 +1152,7 class Oblique_Gauss_Fit(Operation):
1152 #return A1f, B1f, C1f, A2f, B2f, C2f, K2f, Df, doppler
1152 #return A1f, B1f, C1f, A2f, B2f, C2f, K2f, Df, doppler
1153 return A1f, B1f, C1f, A2f, B2f, C2f, K2f, Df, doppler
1153 return A1f, B1f, C1f, A2f, B2f, C2f, K2f, Df, doppler
1154
1154
1155 def Double_Gauss_Double_Skew_fit_weight_bound_no_inputs(self,spc,freq,Nincoh):
1155 def Double_Gauss_Double_Skew_fit_weight_bound_no_inputs(self,spc,freq,Nincoh,hei):
1156
1156
1157 from scipy.optimize import least_squares
1157 from scipy.optimize import least_squares
1158
1158
@@ -1176,7 +1176,7 class Oblique_Gauss_Fit(Operation):
1176 #print(a1,b1,c1,a2,b2,c2,k2,d)
1176 #print(a1,b1,c1,a2,b2,c2,k2,d)
1177 #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,-340,numpy.inf,0,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf])
1178 #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])
1178 #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])
1179 bounds=([0,-numpy.inf,0,-5,0,-400,0,0,0],[numpy.inf,-300,numpy.inf,5,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf])
1179 bounds=([0,-numpy.inf,0,-5,0,-400,0,0,0],[numpy.inf,-200,numpy.inf,5,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf])
1180
1180
1181 #print(bounds)
1181 #print(bounds)
1182 #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])
1182 #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])
@@ -1194,10 +1194,20 class Oblique_Gauss_Fit(Operation):
1194 x0_value = numpy.array([spc_max,dop1_x0,30,-.1,spc_max/4, dop2_x0,150,1,1.0e7])
1194 x0_value = numpy.array([spc_max,dop1_x0,30,-.1,spc_max/4, dop2_x0,150,1,1.0e7])
1195 #x0_value = numpy.array([spc_max,-400.5,30,-.1,spc_max/4,-200.5,150,1,1.0e7])
1195 #x0_value = numpy.array([spc_max,-400.5,30,-.1,spc_max/4,-200.5,150,1,1.0e7])
1196 #popt = least_squares(lsq_func,[a1,b1,c1,a2,b2,c2,k2,d],x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=1)
1196 #popt = least_squares(lsq_func,[a1,b1,c1,a2,b2,c2,k2,d],x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=1)
1197 '''
1198 print("INSIDE 1")
1199 print("x0_value: ", x0_value)
1200 print("boundaries: ", bounds)
1201 import matplotlib.pyplot as plt
1202 plt.plot(freq,spc)
1203 plt.plot(freq,self.double_gaussian_double_skew(freq,x0_value[0],x0_value[1],x0_value[2],x0_value[3],x0_value[4],x0_value[5],x0_value[6],x0_value[7],x0_value[8]))
1204 plt.title(hei)
1205 plt.show()
1206 '''
1197 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1207 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1198 # popt = least_squares(lsq_func,[a1,b1,c1,a2,b2,c2,k2,d],x_scale=params_scale,verbose=1)
1208 # popt = least_squares(lsq_func,[a1,b1,c1,a2,b2,c2,k2,d],x_scale=params_scale,verbose=1)
1199 #print(popt)
1209 #print(popt)
1200
1210 #########print("INSIDE 2")
1201 J = popt.jac
1211 J = popt.jac
1202
1212
1203 try:
1213 try:
@@ -1223,6 +1233,8 class Oblique_Gauss_Fit(Operation):
1223 doppler2 = freq[numpy.argmax(aux2)]
1233 doppler2 = freq[numpy.argmax(aux2)]
1224 #print("error",error)
1234 #print("error",error)
1225 #exit(1)
1235 #exit(1)
1236
1237
1226 return A1f, B1f, C1f, K1f, A2f, B2f, C2f, K2f, Df, doppler1, doppler2, error
1238 return A1f, B1f, C1f, K1f, A2f, B2f, C2f, K2f, Df, doppler1, doppler2, error
1227
1239
1228 def Double_Gauss_fit_weight_bound_no_inputs(self,spc,freq,Nincoh):
1240 def Double_Gauss_fit_weight_bound_no_inputs(self,spc,freq,Nincoh):
@@ -1468,6 +1480,87 class Oblique_Gauss_Fit(Operation):
1468
1480
1469 return A1f, B1f, C1f, Df, error
1481 return A1f, B1f, C1f, Df, error
1470
1482
1483 def clean_outliers(self,param):
1484
1485 threshold = 700
1486
1487 param = numpy.where(param < -threshold, numpy.nan, param)
1488 param = numpy.where(param > +threshold, numpy.nan, param)
1489
1490 return param
1491
1492 def windowing_single(self,spc,x,A,B,C,D,nFFTPoints):
1493 from scipy.optimize import curve_fit,fmin
1494
1495 def R_gaussian(x, a, b, c):
1496 N = int(numpy.shape(x)[0])
1497 val = a * numpy.exp(-((x)*c*2*2*numpy.pi)**2 / (2))* numpy.exp(1.j*b*x*4*numpy.pi)
1498 return val
1499
1500 def T(x,N):
1501 T = 1-abs(x)/N
1502 return T
1503
1504 def R_T_spc_fun(x, a, b, c, d, nFFTPoints):
1505
1506 N = int(numpy.shape(x)[0])
1507
1508 x_max = x[-1]
1509
1510 x_pos = x[int(nFFTPoints/2):]
1511 x_neg = x[:int(nFFTPoints/2)]
1512
1513 R_T_neg_1 = R_gaussian(x, a, b, c)[:int(nFFTPoints/2)]*T(x_neg,-x[0])
1514 R_T_pos_1 = R_gaussian(x, a, b, c)[int(nFFTPoints/2):]*T(x_pos,x[-1])
1515 R_T_sum_1 = R_T_pos_1 + R_T_neg_1
1516 R_T_spc_1 = numpy.fft.fft(R_T_sum_1).real
1517 R_T_spc_1 = numpy.fft.fftshift(R_T_spc_1)
1518 max_val_1 = numpy.max(R_T_spc_1)
1519 R_T_spc_1 = R_T_spc_1*a/max_val_1
1520
1521 R_T_d = d*numpy.fft.fftshift(signal.unit_impulse(N))
1522 R_T_d_neg = R_T_d[:int(nFFTPoints/2)]*T(x_neg,-x[0])
1523 R_T_d_pos = R_T_d[int(nFFTPoints/2):]*T(x_pos,x[-1])
1524 R_T_d_sum = R_T_d_pos + R_T_d_neg
1525 R_T_spc_3 = numpy.fft.fft(R_T_d_sum).real
1526 R_T_spc_3 = numpy.fft.fftshift(R_T_spc_3)
1527
1528 R_T_final = R_T_spc_1 + R_T_spc_3
1529
1530 return R_T_final
1531
1532 y = spc#gaussian(x, a, meanY, sigmaY) + a*0.1*numpy.random.normal(0, 1, size=len(x))
1533
1534 from scipy.stats import norm
1535 mean,std=norm.fit(spc)
1536
1537 # estimate starting values from the data
1538 a = A
1539 b = B
1540 c = C#numpy.std(spc)
1541 d = D
1542 '''
1543 ippSeconds = 250*20*1.e-6/3
1544
1545 x_t = ippSeconds * (numpy.arange(1600) -1600 / 2.)
1546
1547 x_t = numpy.linspace(x_t[0],x_t[-1],3200)
1548
1549 x_freq = numpy.fft.fftfreq(1600,d=ippSeconds)
1550 x_freq = numpy.fft.fftshift(x_freq)
1551 '''
1552 # define a least squares function to optimize
1553 def minfunc(params):
1554 return sum((y-R_T_spc_fun(x,params[0],params[1],params[2],params[3],params[4],params[5],params[6]))**2/1)#y**2)
1555
1556 # fit
1557
1558 popt_full = fmin(minfunc,[a,b,c,d],full_output=True)
1559 #print("nIter", popt_full[2])
1560 popt = popt_full[0]
1561
1562 #return R_T_spc_fun(x_t,popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6]), popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6]
1563 return popt[0], popt[1], popt[2], popt[3]
1471
1564
1472 def run(self, dataOut, mode = 0, Hmin1 = None, Hmax1 = None, Hmin2 = None, Hmax2 = None, Dop = 'Shift'):
1565 def run(self, dataOut, mode = 0, Hmin1 = None, Hmax1 = None, Hmin2 = None, Hmax2 = None, Dop = 'Shift'):
1473
1566
@@ -1565,10 +1658,17 class Oblique_Gauss_Fit(Operation):
1565 pass
1658 pass
1566
1659
1567 else:
1660 else:
1568
1661 #print("After: ", dataOut.data_snr[0])
1662 #######import matplotlib.pyplot as plt
1663 #######plt.plot(dataOut.data_snr[0],dataOut.heightList,marker='*',linestyle='--')
1664 #######plt.show()
1665 #print("l1: ", dataOut.heightList[l1])
1666 #print("l2: ", dataOut.heightList[l2])
1569 for hei in itertools.chain(l1, l2):
1667 for hei in itertools.chain(l1, l2):
1570 #for hei in range(79,81):
1668 #for hei in range(79,81):
1571 if numpy.isnan(dataOut.data_snr[0,hei]):
1669 #if numpy.isnan(dataOut.data_snr[0,hei]) or numpy.isnan(numpy.log10(dataOut.data_snr[0,hei])):
1670 if numpy.isnan(dataOut.snl[0,hei]) or dataOut.snl[0,hei]<.0:
1671
1572 continue #Avoids the analysis when there is only noise
1672 continue #Avoids the analysis when there is only noise
1573
1673
1574 try:
1674 try:
@@ -1592,8 +1692,11 class Oblique_Gauss_Fit(Operation):
1592 from scipy.signal import medfilt
1692 from scipy.signal import medfilt
1593 spcm = medfilt(spc,11)
1693 spcm = medfilt(spc,11)
1594 if x[numpy.argmax(spcm)] <= 0:
1694 if x[numpy.argmax(spcm)] <= 0:
1595 #print("EEJ")
1695 #print("EEJ", dataOut.heightList[hei], hei)
1596 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(spcm,x,dataOut.nIncohInt)
1696 #if hei != 70:
1697 #continue
1698 #else:
1699 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(spcm,x,dataOut.nIncohInt,dataOut.heightList[hei])
1597 #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:
1700 #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:
1598 # dataOut.Oblique_params[0,:,hei] *= numpy.NAN
1701 # dataOut.Oblique_params[0,:,hei] *= numpy.NAN
1599 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,10,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei]))
1702 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,10,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei]))
@@ -1672,7 +1775,11 class Oblique_Gauss_Fit(Operation):
1672 dataOut.paramInterval = dataOut.nProfiles*dataOut.nCohInt*dataOut.ippSeconds
1775 dataOut.paramInterval = dataOut.nProfiles*dataOut.nCohInt*dataOut.ippSeconds
1673 dataOut.lat=-11.95
1776 dataOut.lat=-11.95
1674 dataOut.lon=-76.87
1777 dataOut.lon=-76.87
1675
1778 '''
1779 dataOut.Oblique_params = numpy.where(dataOut.Oblique_params<-700, numpy.nan, dop_t1)
1780 dataOut.Oblique_params = numpy.where(dataOut.Oblique_params<+700, numpy.nan, dop_t1)
1781 Aquí debo exceptuar las amplitudes
1782 '''
1676 if mode == 9: #Double Skew Gaussian
1783 if mode == 9: #Double Skew Gaussian
1677 #dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,-2,:] #Pos[Max_value]
1784 #dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,-2,:] #Pos[Max_value]
1678 #dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,1,:] #Shift
1785 #dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,1,:] #Shift
@@ -1703,9 +1810,154 class Oblique_Gauss_Fit(Operation):
1703 dataOut.Err_Dop_EEJ_T2 = dataOut.Oblique_param_errors[:,4,:]
1810 dataOut.Err_Dop_EEJ_T2 = dataOut.Oblique_param_errors[:,4,:]
1704 dataOut.Err_Spec_W_T2 = dataOut.Oblique_param_errors[:,5,:]
1811 dataOut.Err_Spec_W_T2 = dataOut.Oblique_param_errors[:,5,:]
1705
1812
1813 #print("Before: ", dataOut.Dop_EEJ_T2)
1814 dataOut.Spec_W_T1 = self.clean_outliers(dataOut.Spec_W_T1)
1815 dataOut.Spec_W_T2 = self.clean_outliers(dataOut.Spec_W_T2)
1816 dataOut.Dop_EEJ_T1 = self.clean_outliers(dataOut.Dop_EEJ_T1)
1817 dataOut.Dop_EEJ_T2 = self.clean_outliers(dataOut.Dop_EEJ_T2)
1818 #print("After: ", dataOut.Dop_EEJ_T2)
1819 dataOut.Err_Spec_W_T1 = self.clean_outliers(dataOut.Err_Spec_W_T1)
1820 dataOut.Err_Spec_W_T2 = self.clean_outliers(dataOut.Err_Spec_W_T2)
1821 dataOut.Err_Dop_EEJ_T1 = self.clean_outliers(dataOut.Err_Dop_EEJ_T1)
1822 dataOut.Err_Dop_EEJ_T2 = self.clean_outliers(dataOut.Err_Dop_EEJ_T2)
1823 #print("Before data_snr: ", dataOut.data_snr)
1824 #dataOut.data_snr = numpy.where(numpy.isnan(dataOut.Dop_EEJ_T1), numpy.nan, dataOut.data_snr)
1825 dataOut.snl = numpy.where(numpy.isnan(dataOut.Dop_EEJ_T1), numpy.nan, dataOut.snl)
1826
1827 #print("After data_snr: ", dataOut.data_snr)
1706 dataOut.mode = mode
1828 dataOut.mode = mode
1707 dataOut.flagNoData = numpy.all(numpy.isnan(dataOut.Dop_EEJ_T1)) #Si todos los valores son NaN no se prosigue
1829 dataOut.flagNoData = numpy.all(numpy.isnan(dataOut.Dop_EEJ_T1)) #Si todos los valores son NaN no se prosigue
1708 #dataOut.flagNoData = False #Descomentar solo para ploteo sino mantener comentado
1830 ###dataOut.flagNoData = False #Descomentar solo para ploteo sino mantener comentado (para guardado)
1831
1832 return dataOut
1833
1834 class Gaussian_Windowed(Operation):
1835 '''
1836 Written by R. Flores
1837 '''
1838 def __init__(self):
1839 Operation.__init__(self)
1840
1841 def windowing_single(self,spc,x,A,B,C,D,nFFTPoints):
1842 from scipy.optimize import curve_fit,fmin
1843
1844 def gaussian(x, a, b, c, d):
1845 val = a * numpy.exp(-(x - b)**2 / (2*c**2)) + d
1846 return val
1847
1848 def R_gaussian(x, a, b, c):
1849 N = int(numpy.shape(x)[0])
1850 val = a * numpy.exp(-((x)*c*2*2*numpy.pi)**2 / (2))* numpy.exp(1.j*b*x*4*numpy.pi)
1851 return val
1852
1853 def T(x,N):
1854 T = 1-abs(x)/N
1855 return T
1856
1857 def R_T_spc_fun(x, a, b, c, d, nFFTPoints):
1858
1859 N = int(numpy.shape(x)[0])
1860
1861 x_max = x[-1]
1862
1863 x_pos = x[nFFTPoints:]
1864 x_neg = x[:nFFTPoints]
1865 #print([int(nFFTPoints/2))
1866 #print("x: ", x)
1867 #print("x_neg: ", x_neg)
1868 #print("x_pos: ", x_pos)
1869
1870
1871 R_T_neg_1 = R_gaussian(x, a, b, c)[:nFFTPoints]*T(x_neg,-x[0])
1872 R_T_pos_1 = R_gaussian(x, a, b, c)[nFFTPoints:]*T(x_pos,x[-1])
1873 #print(T(x_pos,x[-1]),x_pos,x[-1])
1874 #print(R_T_neg_1.shape,R_T_pos_1.shape)
1875 R_T_sum_1 = R_T_pos_1 + R_T_neg_1
1876 R_T_spc_1 = numpy.fft.fft(R_T_sum_1).real
1877 R_T_spc_1 = numpy.fft.fftshift(R_T_spc_1)
1878 max_val_1 = numpy.max(R_T_spc_1)
1879 R_T_spc_1 = R_T_spc_1*a/max_val_1
1880
1881 R_T_d = d*numpy.fft.fftshift(signal.unit_impulse(N))
1882 R_T_d_neg = R_T_d[:nFFTPoints]*T(x_neg,-x[0])
1883 R_T_d_pos = R_T_d[nFFTPoints:]*T(x_pos,x[-1])
1884 R_T_d_sum = R_T_d_pos + R_T_d_neg
1885 R_T_spc_3 = numpy.fft.fft(R_T_d_sum).real
1886 R_T_spc_3 = numpy.fft.fftshift(R_T_spc_3)
1887
1888 R_T_final = R_T_spc_1 + R_T_spc_3
1889
1890 return R_T_final
1891
1892 y = spc#gaussian(x, a, meanY, sigmaY) + a*0.1*numpy.random.normal(0, 1, size=len(x))
1893
1894 from scipy.stats import norm
1895 mean,std=norm.fit(spc)
1896
1897 # estimate starting values from the data
1898 a = A
1899 b = B
1900 c = C#numpy.std(spc)
1901 d = D
1902 #'''
1903 #ippSeconds = 250*20*1.e-6/3
1904
1905 #x_t = ippSeconds * (numpy.arange(nFFTPoints) - nFFTPoints / 2.)
1906
1907 #x_t = numpy.linspace(x_t[0],x_t[-1],3200)
1908 #print("x_t: ", x_t)
1909 #print("nFFTPoints: ", nFFTPoints)
1910 x_vel = numpy.linspace(x[0],x[-1],int(2*nFFTPoints))
1911 #print("x_vel: ", x_vel)
1912 #x_freq = numpy.fft.fftfreq(1600,d=ippSeconds)
1913 #x_freq = numpy.fft.fftshift(x_freq)
1914 #'''
1915 # define a least squares function to optimize
1916 def minfunc(params):
1917 #print("y.shape: ", numpy.shape(y))
1918 return sum((y-R_T_spc_fun(x_vel,params[0],params[1],params[2],params[3],nFFTPoints))**2/1)#y**2)
1919
1920 # fit
1921
1922 popt_full = fmin(minfunc,[a,b,c,d], disp=False)
1923 #print("nIter", popt_full[2])
1924 popt = popt_full#[0]
1925
1926 fun = gaussian(x, popt[0], popt[1], popt[2], popt[3])
1927
1928 #return R_T_spc_fun(x_t,popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6]), popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6]
1929 return fun, popt[0], popt[1], popt[2], popt[3]
1930
1931 def run(self, dataOut):
1932
1933 from scipy.signal import medfilt
1934 import matplotlib.pyplot as plt
1935 dataOut.moments = numpy.ones((dataOut.nChannels,4,dataOut.nHeights))*numpy.NAN
1936 dataOut.VelRange = dataOut.getVelRange(0)
1937 for nChannel in range(dataOut.nChannels):
1938 for hei in range(dataOut.heightList.shape[0]):
1939 #print("ipp: ", dataOut.ippSeconds)
1940 spc = numpy.copy(dataOut.data_spc[nChannel,:,hei])
1941
1942 #print(VelRange)
1943 #print(dataOut.getFreqRange(64))
1944 spcm = medfilt(spc,11)
1945 spc_max = numpy.max(spcm)
1946 dop1_x0 = dataOut.VelRange[numpy.argmax(spcm)]
1947 D = numpy.min(spcm)
1948
1949 fun, A, B, C, D = self.windowing_single(spc,dataOut.VelRange,spc_max,dop1_x0,abs(dop1_x0),D,dataOut.nFFTPoints)
1950 dataOut.moments[nChannel,0,hei] = A
1951 dataOut.moments[nChannel,1,hei] = B
1952 dataOut.moments[nChannel,2,hei] = C
1953 dataOut.moments[nChannel,3,hei] = D
1954 '''
1955 plt.figure()
1956 plt.plot(VelRange,spc,marker='*',linestyle='')
1957 plt.plot(VelRange,fun)
1958 plt.title(dataOut.heightList[hei])
1959 plt.show()
1960 '''
1709
1961
1710 return dataOut
1962 return dataOut
1711
1963
@@ -2547,7 +2799,7 class SpectralFitting(Operation):
2547 #def __DiffCoherent(self,snrth, spectra, cspectra, nProf, heights,nChan, nHei, nPairs, channels, noise, crosspairs):
2799 #def __DiffCoherent(self,snrth, spectra, cspectra, nProf, heights,nChan, nHei, nPairs, channels, noise, crosspairs):
2548 def __DiffCoherent(self, spectra, cspectra, dataOut, noise, snrth, coh_th, hei_th):
2800 def __DiffCoherent(self, spectra, cspectra, dataOut, noise, snrth, coh_th, hei_th):
2549
2801
2550 import matplotlib.pyplot as plt
2802 #import matplotlib.pyplot as plt
2551 nProf = dataOut.nProfiles
2803 nProf = dataOut.nProfiles
2552 heights = dataOut.heightList
2804 heights = dataOut.heightList
2553 nHei = len(heights)
2805 nHei = len(heights)
@@ -2573,20 +2825,19 class SpectralFitting(Operation):
2573
2825
2574 if coh_th == None : coh_th = numpy.array([0.75,0.65,0.15]) # 0.65
2826 if coh_th == None : coh_th = numpy.array([0.75,0.65,0.15]) # 0.65
2575 if hei_th == None : hei_th = numpy.array([60,300,650])
2827 if hei_th == None : hei_th = numpy.array([60,300,650])
2576 for ic in range(2):
2828 for ic in range(nPairs):
2577 pair = crosspairs[ic]
2829 pair = crosspairs[ic]
2578 #si el SNR es mayor que el SNR threshold los datos se toman coherentes
2830 #si el SNR es mayor que el SNR threshold los datos se toman coherentes
2579 s_n0 = power[pair[0],:]/noise[pair[0]]
2831 s_n0 = power[pair[0],:]/noise[pair[0]]
2580 s_n1 = power[pair[1],:]/noise[pair[1]]
2832 s_n1 = power[pair[1],:]/noise[pair[1]]
2581
2582 valid1 =(s_n0>=snr_th).nonzero()
2833 valid1 =(s_n0>=snr_th).nonzero()
2583 valid2 = (s_n1>=snr_th).nonzero()
2834 valid2 = (s_n1>=snr_th).nonzero()
2584 #valid = valid2 + valid1 #numpy.concatenate((valid1,valid2), axis=None)
2835
2585 valid1 = numpy.array(valid1[0])
2836 valid1 = numpy.array(valid1[0])
2586 valid2 = numpy.array(valid2[0])
2837 valid2 = numpy.array(valid2[0])
2587 valid = valid1
2838 valid = valid1
2588 for iv in range(len(valid2)):
2839 for iv in range(len(valid2)):
2589 #for ivv in range(len(valid1)) :
2840
2590 indv = numpy.array((valid1 == valid2[iv]).nonzero())
2841 indv = numpy.array((valid1 == valid2[iv]).nonzero())
2591 if len(indv[0]) == 0 :
2842 if len(indv[0]) == 0 :
2592 valid = numpy.concatenate((valid,valid2[iv]), axis=None)
2843 valid = numpy.concatenate((valid,valid2[iv]), axis=None)
@@ -2594,17 +2845,16 class SpectralFitting(Operation):
2594 my_coh_aver[pair[0],valid]=1
2845 my_coh_aver[pair[0],valid]=1
2595 my_coh_aver[pair[1],valid]=1
2846 my_coh_aver[pair[1],valid]=1
2596 # si la coherencia es mayor a la coherencia threshold los datos se toman
2847 # si la coherencia es mayor a la coherencia threshold los datos se toman
2597 #print my_coh_aver[0,:]
2848
2598 coh = numpy.squeeze(numpy.nansum(cspectra[ic,:,:], axis=0)/numpy.sqrt(numpy.nansum(spectra[pair[0],:,:], axis=0)*numpy.nansum(spectra[pair[1],:,:], axis=0)))
2849 coh = numpy.squeeze(numpy.nansum(cspectra[ic,:,:], axis=0)/numpy.sqrt(numpy.nansum(spectra[pair[0],:,:], axis=0)*numpy.nansum(spectra[pair[1],:,:], axis=0)))
2599 #print('coh',numpy.absolute(coh))
2850
2600 for ih in range(len(hei_th)):
2851 for ih in range(len(hei_th)):
2601 hvalid = (heights>hei_th[ih]).nonzero()
2852 hvalid = (heights>hei_th[ih]).nonzero()
2602 hvalid = hvalid[0]
2853 hvalid = hvalid[0]
2603 if len(hvalid)>0:
2854 if len(hvalid)>0:
2604 valid = (numpy.absolute(coh[hvalid])>coh_th[ih]).nonzero()
2855 valid = (numpy.absolute(coh[hvalid])>coh_th[ih]).nonzero()
2605 valid = valid[0]
2856 valid = valid[0]
2606 #print('hvalid:',hvalid)
2857
2607 #print('valid', valid)
2608 if len(valid)>0:
2858 if len(valid)>0:
2609 my_coh_aver[pair[0],hvalid[valid]] =1
2859 my_coh_aver[pair[0],hvalid[valid]] =1
2610 my_coh_aver[pair[1],hvalid[valid]] =1
2860 my_coh_aver[pair[1],hvalid[valid]] =1
@@ -2620,7 +2870,7 class SpectralFitting(Operation):
2620 my_incoh_aver[pair[1],incoh_echoes] = 1
2870 my_incoh_aver[pair[1],incoh_echoes] = 1
2621
2871
2622
2872
2623 for ic in range(2):
2873 for ic in range(nPairs):
2624 pair = crosspairs[ic]
2874 pair = crosspairs[ic]
2625
2875
2626 valid1 =(my_coh_aver[pair[0],:]==1 ).nonzero()
2876 valid1 =(my_coh_aver[pair[0],:]==1 ).nonzero()
@@ -2628,29 +2878,25 class SpectralFitting(Operation):
2628 valid1 = numpy.array(valid1[0])
2878 valid1 = numpy.array(valid1[0])
2629 valid2 = numpy.array(valid2[0])
2879 valid2 = numpy.array(valid2[0])
2630 valid = valid1
2880 valid = valid1
2631 #print valid1 , valid2
2881
2632 for iv in range(len(valid2)):
2882 for iv in range(len(valid2)):
2633 #for ivv in range(len(valid1)) :
2883
2634 indv = numpy.array((valid1 == valid2[iv]).nonzero())
2884 indv = numpy.array((valid1 == valid2[iv]).nonzero())
2635 if len(indv[0]) == 0 :
2885 if len(indv[0]) == 0 :
2636 valid = numpy.concatenate((valid,valid2[iv]), axis=None)
2886 valid = numpy.concatenate((valid,valid2[iv]), axis=None)
2637 #print valid
2638 #valid = numpy.concatenate((valid1,valid2), axis=None)
2639 valid1 =(my_coh_aver[pair[0],:] !=1 ).nonzero()
2887 valid1 =(my_coh_aver[pair[0],:] !=1 ).nonzero()
2640 valid2 = (my_coh_aver[pair[1],:] !=1).nonzero()
2888 valid2 = (my_coh_aver[pair[1],:] !=1).nonzero()
2641 valid1 = numpy.array(valid1[0])
2889 valid1 = numpy.array(valid1[0])
2642 valid2 = numpy.array(valid2[0])
2890 valid2 = numpy.array(valid2[0])
2643 incoh_echoes = valid1
2891 incoh_echoes = valid1
2644 #print valid1, valid2
2892
2645 #incoh_echoes= numpy.concatenate((valid1,valid2), axis=None)
2646 for iv in range(len(valid2)):
2893 for iv in range(len(valid2)):
2647 #for ivv in range(len(valid1)) :
2894
2648 indv = numpy.array((valid1 == valid2[iv]).nonzero())
2895 indv = numpy.array((valid1 == valid2[iv]).nonzero())
2649 if len(indv[0]) == 0 :
2896 if len(indv[0]) == 0 :
2650 incoh_echoes = numpy.concatenate(( incoh_echoes,valid2[iv]), axis=None)
2897 incoh_echoes = numpy.concatenate(( incoh_echoes,valid2[iv]), axis=None)
2651 #print incoh_echoes
2898
2652 if len(valid)>0:
2899 if len(valid)>0:
2653 #print pair
2654 coh_spectra[pair[0],:,valid] = spectra[pair[0],:,valid]
2900 coh_spectra[pair[0],:,valid] = spectra[pair[0],:,valid]
2655 coh_spectra[pair[1],:,valid] = spectra[pair[1],:,valid]
2901 coh_spectra[pair[1],:,valid] = spectra[pair[1],:,valid]
2656 coh_cspectra[ic,:,valid] = cspectra[ic,:,valid]
2902 coh_cspectra[ic,:,valid] = cspectra[ic,:,valid]
@@ -2662,20 +2908,12 class SpectralFitting(Operation):
2662 incoh_cspectra[ic,:,incoh_echoes] = cspectra[ic,:,incoh_echoes]
2908 incoh_cspectra[ic,:,incoh_echoes] = cspectra[ic,:,incoh_echoes]
2663 incoh_aver[pair[0],incoh_echoes]=1
2909 incoh_aver[pair[0],incoh_echoes]=1
2664 incoh_aver[pair[1],incoh_echoes]=1
2910 incoh_aver[pair[1],incoh_echoes]=1
2665 #plt.imshow(spectra[0,:,:],vmin=20000000)
2911
2666 #plt.show()
2667 #my_incoh_aver = my_incoh_aver+1
2668
2669 #spec = my_incoh_spectra.copy()
2670 #cspec = my_incoh_cspectra.copy()
2671 #print('######################', spec)
2672 #print(self.numpy)
2673 #return spec, cspec,coh_aver
2674 return my_incoh_spectra ,my_incoh_cspectra,my_incoh_aver,my_coh_aver, incoh_spectra, coh_spectra, incoh_cspectra, coh_cspectra, incoh_aver, coh_aver
2912 return my_incoh_spectra ,my_incoh_cspectra,my_incoh_aver,my_coh_aver, incoh_spectra, coh_spectra, incoh_cspectra, coh_cspectra, incoh_aver, coh_aver
2675
2913
2676 def __CleanCoherent(self,snrth, spectra, cspectra, coh_aver,dataOut, noise,clean_coh_echoes,index):
2914 def __CleanCoherent(self,snrth, spectra, cspectra, coh_aver,dataOut, noise,clean_coh_echoes,index):
2677
2915
2678 import matplotlib.pyplot as plt
2916 #import matplotlib.pyplot as plt
2679 nProf = dataOut.nProfiles
2917 nProf = dataOut.nProfiles
2680 heights = dataOut.heightList
2918 heights = dataOut.heightList
2681 nHei = len(heights)
2919 nHei = len(heights)
@@ -2684,15 +2922,9 class SpectralFitting(Operation):
2684 crosspairs = dataOut.groupList
2922 crosspairs = dataOut.groupList
2685 nPairs = len(crosspairs)
2923 nPairs = len(crosspairs)
2686
2924
2687 #data = dataOut.data_pre[0]
2688 absc = dataOut.abscissaList[:-1]
2925 absc = dataOut.abscissaList[:-1]
2689 #noise = dataOut.noise
2690 #nChannel = data.shape[0]
2691 data_param = numpy.zeros((nChan, 4, spectra.shape[2]))
2926 data_param = numpy.zeros((nChan, 4, spectra.shape[2]))
2692
2927
2693
2694 #plt.plot(absc)
2695 #plt.show()
2696 clean_coh_spectra = spectra.copy()
2928 clean_coh_spectra = spectra.copy()
2697 clean_coh_cspectra = cspectra.copy()
2929 clean_coh_cspectra = cspectra.copy()
2698 clean_coh_aver = coh_aver.copy()
2930 clean_coh_aver = coh_aver.copy()
@@ -2703,17 +2935,14 class SpectralFitting(Operation):
2703 rtime0 = [6,18] # periodo sin ESF
2935 rtime0 = [6,18] # periodo sin ESF
2704 rtime1 = [10.5,13.5] # periodo con alta coherencia y alto ancho espectral (esperado): SOL.
2936 rtime1 = [10.5,13.5] # periodo con alta coherencia y alto ancho espectral (esperado): SOL.
2705
2937
2706 time = index*5./60
2938 time = index*5./60 # en base a 5 min de proceso
2707 if clean_coh_echoes == 1 :
2939 if clean_coh_echoes == 1 :
2708 for ind in range(nChan):
2940 for ind in range(nChan):
2709 data_param[ind,:,:] = self.__calculateMoments( spectra[ind,:,:] , absc , noise[ind] )
2941 data_param[ind,:,:] = self.__calculateMoments( spectra[ind,:,:] , absc , noise[ind] )
2710 #print data_param[:,3]
2942
2711 spwd = data_param[:,3]
2943 spwd = data_param[:,3]
2712 #print spwd.shape
2944
2713 # SPECB_JULIA,header=anal_header,jspectra=spectra,vel=velocities,hei=heights, num_aver=1, mode_fit=0,smoothing=smoothing,jvelr=velr,jspwd=spwd,jsnr=snr,jnoise=noise,jstdvnoise=stdvnoise
2945 # SPECB_JULIA,header=anal_header,jspectra=spectra,vel=velocities,hei=heights, num_aver=1, mode_fit=0,smoothing=smoothing,jvelr=velr,jspwd=spwd,jsnr=snr,jnoise=noise,jstdvnoise=stdvnoise
2714 #spwd1=[ 1.65607, 1.43416, 0.500373, 0.208361, 0.000000, 26.7767, 22.5936, 26.7530, 20.6962, 29.1098, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 28.0300, 27.0511, 27.8810, 26.3126, 27.8445, 24.6181, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000]
2715 #spwd=numpy.array([spwd1,spwd1,spwd1,spwd1])
2716 #print spwd.shape, heights.shape,coh_aver.shape
2717 # para obtener spwd
2946 # para obtener spwd
2718 for ic in range(nPairs):
2947 for ic in range(nPairs):
2719 pair = crosspairs[ic]
2948 pair = crosspairs[ic]
@@ -2748,29 +2977,20 class SpectralFitting(Operation):
2748 self.i=0
2977 self.i=0
2749 self.isConfig = False
2978 self.isConfig = False
2750
2979
2751
2752 def setup(self,nChan,nProf,nHei,nBlocks):
2980 def setup(self,nChan,nProf,nHei,nBlocks):
2753 self.__dataReady = False
2981 self.__dataReady = False
2754 self.bloques = numpy.zeros([2, nProf, nHei,nBlocks], dtype= complex)
2982 self.bloques = numpy.zeros([2, nProf, nHei,nBlocks], dtype= complex)
2755 self.bloque0 = numpy.zeros([nChan, nProf, nHei, nBlocks])
2983 self.bloque0 = numpy.zeros([nChan, nProf, nHei, nBlocks])
2756
2984
2757 #def CleanRayleigh(self,dataOut,spectra,cspectra,out_spectra,out_cspectra,sat_spectra,sat_cspectra,crosspairs,heights, channels, nProf,nHei,nChan,nPairs,nIncohInt,nBlocks):
2758 def CleanRayleigh(self,dataOut,spectra,cspectra,save_drifts):
2985 def CleanRayleigh(self,dataOut,spectra,cspectra,save_drifts):
2759 #import matplotlib.pyplot as plt
2760 #for k in range(149):
2761
2986
2762 # self.bloque0[:,:,:,k] = spectra[:,:,0:nHei]
2987 rfunc = cspectra.copy()
2763 # self.bloques[:,:,:,k] = cspectra[:,:,0:nHei]
2764 #if self.i==nBlocks:
2765 # self.i==0
2766 rfunc = cspectra.copy() #self.bloques
2767 n_funct = len(rfunc[0,:,0,0])
2988 n_funct = len(rfunc[0,:,0,0])
2768 val_spc = spectra*0.0 #self.bloque0*0.0
2989 val_spc = spectra*0.0
2769 val_cspc = cspectra*0.0 #self.bloques*0.0
2990 val_cspc = cspectra*0.0
2770 in_sat_spectra = spectra.copy() #self.bloque0
2991 in_sat_spectra = spectra.copy()
2771 in_sat_cspectra = cspectra.copy() #self.bloques
2992 in_sat_cspectra = cspectra.copy()
2772
2993
2773 #print( rfunc.shape)
2774 min_hei = 200
2994 min_hei = 200
2775 nProf = dataOut.nProfiles
2995 nProf = dataOut.nProfiles
2776 heights = dataOut.heightList
2996 heights = dataOut.heightList
@@ -2781,20 +3001,19 class SpectralFitting(Operation):
2781 nPairs = len(crosspairs)
3001 nPairs = len(crosspairs)
2782 hval=(heights >= min_hei).nonzero()
3002 hval=(heights >= min_hei).nonzero()
2783 ih=hval[0]
3003 ih=hval[0]
2784 #print numpy.absolute(rfunc[:,0,0,14])
2785 for ih in range(hval[0][0],nHei):
3004 for ih in range(hval[0][0],nHei):
2786 for ifreq in range(nProf):
3005 for ifreq in range(nProf):
2787 for ii in range(n_funct):
3006 for ii in range(n_funct):
2788
3007
2789 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih]))
3008 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih]))
2790 #print numpy.amin(func2clean)
3009
2791 val = (numpy.isfinite(func2clean)==True).nonzero()
3010 val = (numpy.isfinite(func2clean)==True).nonzero()
2792 if len(val)>0:
3011 if len(val)>0:
2793 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
3012 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
2794 if min_val <= -40 : min_val = -40
3013 if min_val <= -40 : min_val = -40
2795 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
3014 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
2796 if max_val >= 200 : max_val = 200
3015 if max_val >= 200 : max_val = 200
2797 #print min_val, max_val
3016
2798 step = 1
3017 step = 1
2799 #Getting bins and the histogram
3018 #Getting bins and the histogram
2800 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
3019 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
@@ -2809,18 +3028,6 class SpectralFitting(Operation):
2809 except:
3028 except:
2810 mode = mean
3029 mode = mean
2811 stdv = sigma
3030 stdv = sigma
2812 # if ih == 14 and ii == 0 and ifreq ==0 :
2813 # print x_dist.shape, y_dist.shape
2814 # print x_dist, y_dist
2815 # print min_val, max_val, binstep
2816 # print func2clean
2817 # print mean,sigma
2818 # mean1,std = norm.fit(y_dist)
2819 # print mean1, std, gauss_fit
2820 # print fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
2821 # 7.84616 53.9307 3.61863
2822 #stdv = 3.61863 # 2.99089
2823 #mode = 53.9307 #7.79008
2824
3031
2825 #Removing echoes greater than mode + 3*stdv
3032 #Removing echoes greater than mode + 3*stdv
2826 factor_stdv = 2.5
3033 factor_stdv = 2.5
@@ -2831,29 +3038,17 class SpectralFitting(Operation):
2831 cross_pairs = crosspairs[ii]
3038 cross_pairs = crosspairs[ii]
2832 #Getting coherent echoes which are removed.
3039 #Getting coherent echoes which are removed.
2833 if len(novall[0]) > 0:
3040 if len(novall[0]) > 0:
2834 #val_spc[(0,1),novall[a],ih] = 1
2835 #val_spc[,(2,3),novall[a],ih] = 1
2836 val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
3041 val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
2837 val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
3042 val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
2838 val_cspc[novall[0],ii,ifreq,ih] = 1
3043 val_cspc[novall[0],ii,ifreq,ih] = 1
2839 #print("OUT NOVALL 1")
2840 #Removing coherent from ISR data
3044 #Removing coherent from ISR data
2841 # if ih == 17 and ii == 0 and ifreq ==0 :
2842 # print spectra[:,cross_pairs[0],ifreq,ih]
2843 spectra[noval,cross_pairs[0],ifreq,ih] = numpy.nan
3045 spectra[noval,cross_pairs[0],ifreq,ih] = numpy.nan
2844 spectra[noval,cross_pairs[1],ifreq,ih] = numpy.nan
3046 spectra[noval,cross_pairs[1],ifreq,ih] = numpy.nan
2845 cspectra[noval,ii,ifreq,ih] = numpy.nan
3047 cspectra[noval,ii,ifreq,ih] = numpy.nan
2846 # if ih == 17 and ii == 0 and ifreq ==0 :
3048 #
2847 # print spectra[:,cross_pairs[0],ifreq,ih]
2848 # print noval, len(noval[0])
2849 # print novall, len(novall[0])
2850 # print factor_stdv*stdv
2851 # print func2clean-mode
2852 # print val_spc[:,cross_pairs[0],ifreq,ih]
2853 # print spectra[:,cross_pairs[0],ifreq,ih]
2854 #no sale es para savedrifts >2
3049 #no sale es para savedrifts >2
2855 ''' channels = channels
3050 ''' channels = dataOut.channelList
2856 cross_pairs = cross_pairs
3051 cross_pairs = dataOut.groupList
2857 #print("OUT NOVALL 2")
3052 #print("OUT NOVALL 2")
2858
3053
2859 vcross0 = (cross_pairs[0] == channels[ii]).nonzero()
3054 vcross0 = (cross_pairs[0] == channels[ii]).nonzero()
@@ -2874,6 +3069,7 class SpectralFitting(Operation):
2874 self.bloques[vcross,ifreq,ih,noval] = numpy.nan
3069 self.bloques[vcross,ifreq,ih,noval] = numpy.nan
2875 '''
3070 '''
2876 #Getting average of the spectra and cross-spectra from incoherent echoes.
3071 #Getting average of the spectra and cross-spectra from incoherent echoes.
3072
2877 out_spectra = numpy.zeros([nChan,nProf,nHei], dtype=float) #+numpy.nan
3073 out_spectra = numpy.zeros([nChan,nProf,nHei], dtype=float) #+numpy.nan
2878 out_cspectra = numpy.zeros([nPairs,nProf,nHei], dtype=complex) #+numpy.nan
3074 out_cspectra = numpy.zeros([nPairs,nProf,nHei], dtype=complex) #+numpy.nan
2879 for ih in range(nHei):
3075 for ih in range(nHei):
@@ -2881,22 +3077,15 class SpectralFitting(Operation):
2881 for ich in range(nChan):
3077 for ich in range(nChan):
2882 tmp = spectra[:,ich,ifreq,ih]
3078 tmp = spectra[:,ich,ifreq,ih]
2883 valid = (numpy.isfinite(tmp[:])==True).nonzero()
3079 valid = (numpy.isfinite(tmp[:])==True).nonzero()
2884 # if ich == 0 and ifreq == 0 and ih == 17 :
2885 # print tmp
2886 # print valid
2887 # print len(valid[0])
2888 #print('TMP',tmp)
2889 if len(valid[0]) >0 :
3080 if len(valid[0]) >0 :
2890 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
3081 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
2891 #for icr in range(nPairs):
3082
2892 for icr in range(nPairs):
3083 for icr in range(nPairs):
2893 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
3084 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
2894 valid = (numpy.isfinite(tmp)==True).nonzero()
3085 valid = (numpy.isfinite(tmp)==True).nonzero()
2895 if len(valid[0]) > 0:
3086 if len(valid[0]) > 0:
2896 out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
3087 out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
2897 # print('##########################################################')
2898 #Removing fake coherent echoes (at least 4 points around the point)
3088 #Removing fake coherent echoes (at least 4 points around the point)
2899
2900 val_spectra = numpy.sum(val_spc,0)
3089 val_spectra = numpy.sum(val_spc,0)
2901 val_cspectra = numpy.sum(val_cspc,0)
3090 val_cspectra = numpy.sum(val_cspc,0)
2902
3091
@@ -2932,13 +3121,6 class SpectralFitting(Operation):
2932 tmp_sat_cspectra = cspectra.copy()
3121 tmp_sat_cspectra = cspectra.copy()
2933 tmp_sat_cspectra = tmp_sat_cspectra*numpy.nan
3122 tmp_sat_cspectra = tmp_sat_cspectra*numpy.nan
2934
3123
2935 # fig = plt.figure(figsize=(6,5))
2936 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
2937 # ax = fig.add_axes([left, bottom, width, height])
2938 # cp = ax.contour(10*numpy.log10(numpy.absolute(spectra[0,0,:,:])))
2939 # ax.clabel(cp, inline=True,fontsize=10)
2940 # plt.show()
2941
2942 val = (val_spc > 0).nonzero()
3124 val = (val_spc > 0).nonzero()
2943 if len(val[0]) > 0:
3125 if len(val[0]) > 0:
2944 tmp_sat_spectra[val] = in_sat_spectra[val]
3126 tmp_sat_spectra[val] = in_sat_spectra[val]
@@ -2963,52 +3145,29 class SpectralFitting(Operation):
2963 valid = (numpy.isfinite(tmp)).nonzero()
3145 valid = (numpy.isfinite(tmp)).nonzero()
2964 if len(valid[0]) > 0:
3146 if len(valid[0]) > 0:
2965 sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
3147 sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
2966 #self.__dataReady= True
3148
2967 #sat_spectra, sat_cspectra= sat_spectra, sat_cspectra
2968 #if not self.__dataReady:
2969 #return None, None
2970 return out_spectra, out_cspectra,sat_spectra,sat_cspectra
3149 return out_spectra, out_cspectra,sat_spectra,sat_cspectra
2971 def REM_ISOLATED_POINTS(self,array,rth):
3150 def REM_ISOLATED_POINTS(self,array,rth):
2972 # import matplotlib.pyplot as plt
2973 if rth == None : rth = 4
3151 if rth == None : rth = 4
2974
2975 num_prof = len(array[0,:,0])
3152 num_prof = len(array[0,:,0])
2976 num_hei = len(array[0,0,:])
3153 num_hei = len(array[0,0,:])
2977 n2d = len(array[:,0,0])
3154 n2d = len(array[:,0,0])
2978
3155
2979 for ii in range(n2d) :
3156 for ii in range(n2d) :
2980 #print ii,n2d
2981 tmp = array[ii,:,:]
3157 tmp = array[ii,:,:]
2982 #print tmp.shape, array[ii,101,:],array[ii,102,:]
2983
2984 # fig = plt.figure(figsize=(6,5))
2985 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
2986 # ax = fig.add_axes([left, bottom, width, height])
2987 # x = range(num_prof)
2988 # y = range(num_hei)
2989 # cp = ax.contour(y,x,tmp)
2990 # ax.clabel(cp, inline=True,fontsize=10)
2991 # plt.show()
2992
2993 #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
2994 tmp = numpy.reshape(tmp,num_prof*num_hei)
3158 tmp = numpy.reshape(tmp,num_prof*num_hei)
2995 indxs1 = (numpy.isfinite(tmp)==True).nonzero()
3159 indxs1 = (numpy.isfinite(tmp)==True).nonzero()
2996 indxs2 = (tmp > 0).nonzero()
3160 indxs2 = (tmp > 0).nonzero()
2997
2998 indxs1 = (indxs1[0])
3161 indxs1 = (indxs1[0])
2999 indxs2 = indxs2[0]
3162 indxs2 = indxs2[0]
3000 #indxs1 = numpy.array(indxs1[0])
3001 #indxs2 = numpy.array(indxs2[0])
3002 indxs = None
3163 indxs = None
3003 #print indxs1 , indxs2
3164
3004 for iv in range(len(indxs2)):
3165 for iv in range(len(indxs2)):
3005 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
3166 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
3006 #print len(indxs2), indv
3007 if len(indv[0]) > 0 :
3167 if len(indv[0]) > 0 :
3008 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
3168 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
3009 # print indxs
3169
3010 indxs = indxs[1:]
3170 indxs = indxs[1:]
3011 #print indxs, len(indxs)
3012 if len(indxs) < 4 :
3171 if len(indxs) < 4 :
3013 array[ii,:,:] = 0.
3172 array[ii,:,:] = 0.
3014 return
3173 return
@@ -3016,7 +3175,6 class SpectralFitting(Operation):
3016 xpos = numpy.mod(indxs ,num_hei)
3175 xpos = numpy.mod(indxs ,num_hei)
3017 ypos = (indxs / num_hei)
3176 ypos = (indxs / num_hei)
3018 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
3177 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
3019 #print sx
3020 xpos = xpos[sx]
3178 xpos = xpos[sx]
3021 ypos = ypos[sx]
3179 ypos = ypos[sx]
3022
3180
@@ -3024,36 +3182,25 class SpectralFitting(Operation):
3024 ic = 0
3182 ic = 0
3025 while True :
3183 while True :
3026 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
3184 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
3027 #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
3185
3028 #plt.plot(r)
3029 #plt.show()
3030 no_coh1 = (numpy.isfinite(r)==True).nonzero()
3186 no_coh1 = (numpy.isfinite(r)==True).nonzero()
3031 no_coh2 = (r <= rth).nonzero()
3187 no_coh2 = (r <= rth).nonzero()
3032 #print r, no_coh1, no_coh2
3033 no_coh1 = numpy.array(no_coh1[0])
3188 no_coh1 = numpy.array(no_coh1[0])
3034 no_coh2 = numpy.array(no_coh2[0])
3189 no_coh2 = numpy.array(no_coh2[0])
3035 no_coh = None
3190 no_coh = None
3036 #print valid1 , valid2
3037 for iv in range(len(no_coh2)):
3191 for iv in range(len(no_coh2)):
3038 indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
3192 indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
3039 if len(indv[0]) > 0 :
3193 if len(indv[0]) > 0 :
3040 no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
3194 no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
3041 no_coh = no_coh[1:]
3195 no_coh = no_coh[1:]
3042 #print len(no_coh), no_coh
3043 if len(no_coh) < 4 :
3196 if len(no_coh) < 4 :
3044 #print xpos[ic], ypos[ic], ic
3045 # plt.plot(r)
3046 # plt.show()
3047 xpos[ic] = numpy.nan
3197 xpos[ic] = numpy.nan
3048 ypos[ic] = numpy.nan
3198 ypos[ic] = numpy.nan
3049
3199
3050 ic = ic + 1
3200 ic = ic + 1
3051 if (ic == len(indxs)) :
3201 if (ic == len(indxs)) :
3052 break
3202 break
3053 #print( xpos, ypos)
3054
3055 indxs = (numpy.isfinite(list(xpos))==True).nonzero()
3203 indxs = (numpy.isfinite(list(xpos))==True).nonzero()
3056 #print indxs[0]
3057 if len(indxs[0]) < 4 :
3204 if len(indxs[0]) < 4 :
3058 array[ii,:,:] = 0.
3205 array[ii,:,:] = 0.
3059 return
3206 return
@@ -3068,27 +3215,11 class SpectralFitting(Operation):
3068 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
3215 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
3069 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
3216 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
3070
3217
3071 #print array.shape
3072 #tmp = numpy.reshape(tmp,(num_prof,num_hei))
3073 #print tmp.shape
3074
3075 # fig = plt.figure(figsize=(6,5))
3076 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
3077 # ax = fig.add_axes([left, bottom, width, height])
3078 # x = range(num_prof)
3079 # y = range(num_hei)
3080 # cp = ax.contour(y,x,array[ii,:,:])
3081 # ax.clabel(cp, inline=True,fontsize=10)
3082 # plt.show()
3083 return array
3218 return array
3084 def moments(self,doppler,yarray,npoints):
3219 def moments(self,doppler,yarray,npoints):
3085 ytemp = yarray
3220 ytemp = yarray
3086 #val = WHERE(ytemp GT 0,cval)
3087 #if cval == 0 : val = range(npoints-1)
3088 val = (ytemp > 0).nonzero()
3221 val = (ytemp > 0).nonzero()
3089 val = val[0]
3222 val = val[0]
3090 #print('hvalid:',hvalid)
3091 #print('valid', valid)
3092 if len(val) == 0 : val = range(npoints-1)
3223 if len(val) == 0 : val = range(npoints-1)
3093
3224
3094 ynew = 0.5*(ytemp[val[0]]+ytemp[val[len(val)-1]])
3225 ynew = 0.5*(ytemp[val[0]]+ytemp[val[len(val)-1]])
@@ -3102,19 +3233,289 class SpectralFitting(Operation):
3102 fmom = numpy.sum(doppler*ytemp)/numpy.sum(ytemp)+(index-(npoints/2-1))*numpy.abs(doppler[1]-doppler[0])
3233 fmom = numpy.sum(doppler*ytemp)/numpy.sum(ytemp)+(index-(npoints/2-1))*numpy.abs(doppler[1]-doppler[0])
3103 smom = numpy.sum(doppler*doppler*ytemp)/numpy.sum(ytemp)
3234 smom = numpy.sum(doppler*doppler*ytemp)/numpy.sum(ytemp)
3104 return [fmom,numpy.sqrt(smom)]
3235 return [fmom,numpy.sqrt(smom)]
3236
3237 def windowing_single_old(self,spc,x,A,B,C,D,nFFTPoints):
3238 '''
3239 Written by R. Flores
3240 '''
3241 from scipy.optimize import curve_fit,fmin
3242
3243 def gaussian(x, a, b, c, d):
3244 val = a * numpy.exp(-(x - b)**2 / (2*c**2)) + d
3245 return val
3246
3247 def R_gaussian(x, a, b, c):
3248 N = int(numpy.shape(x)[0])
3249 val = a * numpy.exp(-((x)*c*2*2*numpy.pi)**2 / (2))* numpy.exp(1.j*b*x*4*numpy.pi)
3250 return val
3251
3252 def T(x,N):
3253 T = 1-abs(x)/N
3254 return T
3255
3256 def R_T_spc_fun(x, a, b, c, d, nFFTPoints):
3257
3258 N = int(numpy.shape(x)[0])
3259
3260 x_max = x[-1]
3261
3262 x_pos = x[nFFTPoints:]
3263 x_neg = x[:nFFTPoints]
3264
3265 R_T_neg_1 = R_gaussian(x, a, b, c)[:nFFTPoints]*T(x_neg,-x[0])
3266 R_T_pos_1 = R_gaussian(x, a, b, c)[nFFTPoints:]*T(x_pos,x[-1])
3267 #print(T(x_pos,x[-1]),x_pos,x[-1])
3268 #print(R_T_neg_1.shape,R_T_pos_1.shape)
3269 R_T_sum_1 = R_T_pos_1 + R_T_neg_1
3270 R_T_spc_1 = numpy.fft.fft(R_T_sum_1).real
3271 R_T_spc_1 = numpy.fft.fftshift(R_T_spc_1)
3272 max_val_1 = numpy.max(R_T_spc_1)
3273 R_T_spc_1 = R_T_spc_1*a/max_val_1
3274 print("R_T_spc_1: ", R_T_spc_1)
3275
3276 R_T_d = d*numpy.fft.fftshift(signal.unit_impulse(N))
3277 R_T_d_neg = R_T_d[:nFFTPoints]*T(x_neg,-x[0])
3278 R_T_d_pos = R_T_d[nFFTPoints:]*T(x_pos,x[-1])
3279 R_T_d_sum = R_T_d_pos + R_T_d_neg
3280 R_T_spc_3 = numpy.fft.fft(R_T_d_sum).real
3281 R_T_spc_3 = numpy.fft.fftshift(R_T_spc_3)
3282
3283 R_T_final = R_T_spc_1# + R_T_spc_3
3284
3285 return R_T_final
3286
3287 y = spc#gaussian(x, a, meanY, sigmaY) + a*0.1*numpy.random.normal(0, 1, size=len(x))
3288
3289 from scipy.stats import norm
3290 mean,std=norm.fit(spc)
3291
3292 # estimate starting values from the data
3293 print("A: ", A)
3294 a = A-D
3295 b = B
3296 c = C#numpy.std(spc) #C
3297 d = D
3298 #'''
3299 #ippSeconds = 250*20*1.e-6/3
3300
3301 #x_t = ippSeconds * (numpy.arange(nFFTPoints) - nFFTPoints / 2.)
3302
3303 #x_t = numpy.linspace(x_t[0],x_t[-1],3200)
3304 #print("x_t: ", x_t)
3305 #print("nFFTPoints: ", nFFTPoints)
3306 x_vel = numpy.linspace(x[0],x[-1],int(2*nFFTPoints))
3307 #print("x_vel: ", x_vel)
3308 #x_freq = numpy.fft.fftfreq(1600,d=ippSeconds)
3309 #x_freq = numpy.fft.fftshift(x_freq)
3310 #'''
3311 # define a least squares function to optimize
3312 import matplotlib.pyplot as plt
3313 aui = R_T_spc_fun(x_vel,a,b,c,d,nFFTPoints)
3314 print("aux_max: ", numpy.nanmax(aui))
3315 #print(dataOut.heightList[hei])
3316 plt.figure()
3317 plt.plot(x,spc,marker='*',linestyle='--')
3318 plt.plot(x,gaussian(x, a, b, c, d),color='b',marker='^',linestyle='')
3319 plt.plot(x,aui,color='k')
3320 #plt.title(dataOut.heightList[hei])
3321 plt.show()
3322
3323 def minfunc(params):
3324 #print("y.shape: ", numpy.shape(y))
3325 return sum((y-R_T_spc_fun(x_vel,params[0],params[1],params[2],params[3],nFFTPoints))**2/1)#y**2)
3326
3327 # fit
3328
3329 popt_full = fmin(minfunc,[a,b,c,d], disp=False)
3330 #print("nIter", popt_full[2])
3331 popt = popt_full#[0]
3332
3333 fun = gaussian(x, popt[0], popt[1], popt[2], popt[3])
3334 print("pop1[0]: ", popt[0])
3335 #return R_T_spc_fun(x_t,popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6]), popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6]
3336 return fun, popt[0], popt[1], popt[2], popt[3]
3337
3338 def windowing_single(self,spc,x,A,B,C,D,nFFTPoints):
3339 '''
3340 Written by R. Flores
3341 '''
3342 from scipy.optimize import curve_fit,fmin
3343
3344 def gaussian(x, a, b, c, d):
3345 val = a * numpy.exp(-(x - b)**2 / (2*c**2)) + d
3346 return val
3347
3348 def R_gaussian(x, a, b, c):
3349 N = int(numpy.shape(x)[0])
3350
3351 val = (a*numpy.exp((-(1/2)*x*(x*c**2 + 2*1.j*b)))/numpy.sqrt(1/c**2))
3352
3353 return val
3354
3355 def T(x,N):
3356 T = 1-abs(x)/N
3357 return T
3358
3359 def R_T_spc_fun(x, a, id_dop, c, d, nFFTPoints):
3360
3361 N = int(numpy.shape(x)[0])
3362 b = 0
3363 x_max = x[-1]
3364
3365 x_pos = x[nFFTPoints:]
3366 x_neg = x[:nFFTPoints]
3367 R_T_neg_1 = R_gaussian(x, a, b, c)[:nFFTPoints]*T(x_neg,-x[0])
3368 R_T_pos_1 = R_gaussian(x, a, b, c)[nFFTPoints:]*T(x_pos,x[-1])
3369
3370 R_T_sum_1 = R_T_pos_1 + R_T_neg_1
3371 R_T_spc_1 = numpy.fft.fft(R_T_sum_1).real
3372 R_T_spc_1 = numpy.fft.fftshift(R_T_spc_1)
3373 max_val_1 = numpy.max(R_T_spc_1)
3374 R_T_spc_1 = R_T_spc_1*a/max_val_1
3375 #raise NotImplementedError
3376 R_T_d = d*numpy.fft.fftshift(signal.unit_impulse(N))
3377 R_T_d_neg = R_T_d[:nFFTPoints]*T(x_neg,-x[0])
3378 R_T_d_pos = R_T_d[nFFTPoints:]*T(x_pos,x[-1])
3379 R_T_d_sum = R_T_d_pos + R_T_d_neg
3380 R_T_spc_3 = numpy.fft.fft(R_T_d_sum).real
3381 R_T_spc_3 = numpy.fft.fftshift(R_T_spc_3)
3382
3383 R_T_final = R_T_spc_1 + R_T_spc_3
3384
3385 id_dop = int(id_dop)
3386
3387 R_T_final = numpy.roll(R_T_final,-id_dop)
3388
3389 return R_T_final
3390
3391 y = spc#gaussian(x, a, meanY, sigmaY) + a*0.1*numpy.random.normal(0, 1, size=len(x))
3392
3393 from scipy.stats import norm
3394 mean,std=norm.fit(spc)
3395
3396 # estimate starting values from the data
3397 a = A-D
3398 b = B
3399 c = C#numpy.std(spc) #C
3400 d = D
3401
3402 id_dop = numpy.argmax(spc)
3403 id_dop = int(spc.shape[0]/2 - id_dop)
3404
3405 x_vel = numpy.linspace(x[0],x[-1],int(2*nFFTPoints))
3406
3407 # define a least squares function to optimize
3408
3409 def minfunc(params):
3410 #print("y.shape: ", numpy.shape(y))
3411 return sum((y-R_T_spc_fun(x_vel,params[0],params[1],params[2],params[3],nFFTPoints))**2/1)#y**2)
3412
3413 # fit
3414 popt_full = fmin(minfunc,[a,id_dop,c,d], disp=False)
3415 popt = popt_full#[0]
3416
3417 fun = gaussian(x, a, 0, popt[2], popt[3])
3418 fun = numpy.roll(fun,-int(popt[1]))
3419
3420 return fun, popt[0], popt[1], popt[2], popt[3]
3421
3422 def windowing_single_direct(self,spc_mod,x,A,B,C,D,nFFTPoints,timeInterval):
3423 '''
3424 Written by R. Flores
3425 '''
3426 from scipy.optimize import curve_fit,fmin
3427
3428 def gaussian(x, a, b, c, d):
3429 val = a * numpy.exp(-(x - b)**2 / (2*c**2)) + d
3430 return val
3431
3432 def R_gaussian(x, a, b, c, d):
3433 N = int(numpy.shape(x)[0])
3434 val = (a*numpy.exp(-2*c**2*x**2 + 2*x*1.j*b))*(numpy.sqrt(2*numpy.pi)*c)/((numpy.pi)) + d*signal.unit_impulse(N)*numpy.shape(x)[0]/2
3435
3436 return 2*val/numpy.shape(val)[0]
3437
3438 def T(x,N):
3439 T = 1-abs(x)/N
3440 return T
3441
3442 def R_T_spc_fun(x, a, b, c, d, nFFTPoints, timeInterval): #"x" should be time
3443
3444 #timeInterval = 2
3445 x_double = numpy.linspace(0,timeInterval,nFFTPoints)
3446 x_double_m = numpy.flip(x_double)
3447 x_double_aux = numpy.linspace(0,x_double[-2],nFFTPoints)
3448 x_double_t = numpy.concatenate((x_double_m,x_double_aux))
3449 x_double_t /= max(x_double_t)
3450
3451
3452 R_T_sum_1 = R_gaussian(x, a, b, c, d)
3453
3454 R_T_sum_1_flip = numpy.copy(numpy.flip(R_T_sum_1))
3455 R_T_sum_1_flip[-1] = R_T_sum_1_flip[0]
3456 R_T_sum_1_flip = numpy.roll(R_T_sum_1_flip,1)
3457
3458 R_T_sum_1_flip.imag *= -1
3459
3460 R_T_sum_1_total = numpy.concatenate((R_T_sum_1,R_T_sum_1_flip))
3461 R_T_sum_1_total *= x_double_t #times trian_fun
3462
3463 R_T_sum_1_total = R_T_sum_1_total[:nFFTPoints] + R_T_sum_1_total[nFFTPoints:]
3464
3465 R_T_spc_1 = numpy.fft.fft(R_T_sum_1_total).real
3466 R_T_spc_1 = numpy.fft.fftshift(R_T_spc_1)
3467
3468 freq = numpy.fft.fftfreq(nFFTPoints, d=timeInterval/nFFTPoints)
3469
3470 freq = numpy.fft.fftshift(freq)
3471
3472 freq *= 6/2 #lambda/2
3473
3474 return R_T_spc_1
3475
3476 y = spc_mod
3477
3478 #from scipy.stats import norm
3479
3480 # estimate starting values from the data
3481
3482 a = A-D
3483 b = B
3484 c = C
3485 d = D
3486
3487 # define a least squares function to optimize
3488 import matplotlib.pyplot as plt
3489 #ippSeconds = 2
3490 t_range = numpy.linspace(0,timeInterval,nFFTPoints)
3491 #aui = R_T_spc_fun(t_range,a,b,c,d,nFFTPoints,timeInterval)
3492
3493 def minfunc(params):
3494 return sum((y-R_T_spc_fun(t_range,params[0],params[1],params[2],params[3],nFFTPoints,timeInterval))**2/1)#y**2)
3495
3496 # fit
3497 popt_full = fmin(minfunc,[a,b,c,d], disp=False)
3498 popt = popt_full
3499
3500 fun = R_T_spc_fun(t_range,popt[0],popt[1],popt[2],popt[3],nFFTPoints,timeInterval)
3501
3502 return fun, popt[0], popt[1], popt[2], popt[3]
3503
3105 # **********************************************************************************************
3504 # **********************************************************************************************
3106 index = 0
3505 index = 0
3107 fint = 0
3506 fint = 0
3108 buffer = 0
3507 buffer = 0
3109 buffer2 = 0
3508 buffer2 = 0
3110 buffer3 = 0
3509 buffer3 = 0
3111 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None):
3510 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None, filec=None,coh_th=None, hei_th=None,taver=None,Gaussian_Windowed=0):
3112 nChannels = dataOut.nChannels
3511 nChannels = dataOut.nChannels
3113 nHeights= dataOut.heightList.size
3512 nHeights= dataOut.heightList.size
3114 nProf = dataOut.nProfiles
3513 nProf = dataOut.nProfiles
3514 if numpy.any(taver): taver=int(taver)
3515 else : taver = 15
3115 tini=time.localtime(dataOut.utctime)
3516 tini=time.localtime(dataOut.utctime)
3116 if (tini.tm_min % 5) == 0 and (tini.tm_sec < 5 and self.fint==0):
3517 if (tini.tm_min % taver) == 0 and (tini.tm_sec < 5 and self.fint==0):
3117 # print tini.tm_min
3518
3118 self.index = 0
3519 self.index = 0
3119 jspc = self.buffer
3520 jspc = self.buffer
3120 jcspc = self.buffer2
3521 jcspc = self.buffer2
@@ -3124,14 +3525,14 class SpectralFitting(Operation):
3124 self.buffer3 = dataOut.noise
3525 self.buffer3 = dataOut.noise
3125 self.fint = 1
3526 self.fint = 1
3126 if numpy.any(jspc) :
3527 if numpy.any(jspc) :
3127 jspc= numpy.reshape(jspc,(int(len(jspc)/4),nChannels,nProf,nHeights))
3528 jspc= numpy.reshape(jspc,(int(len(jspc)/nChannels),nChannels,nProf,nHeights))
3128 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/2),2,nProf,nHeights))
3529 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/int(nChannels/2)),int(nChannels/2),nProf,nHeights))
3129 jnoise= numpy.reshape(jnoise,(int(len(jnoise)/4),nChannels))
3530 jnoise= numpy.reshape(jnoise,(int(len(jnoise)/nChannels),nChannels))
3130 else:
3531 else:
3131 dataOut.flagNoData = True
3532 dataOut.flagNoData = True
3132 return dataOut
3533 return dataOut
3133 else :
3534 else :
3134 if (tini.tm_min % 5) == 0 : self.fint = 1
3535 if (tini.tm_min % taver) == 0 : self.fint = 1
3135 else : self.fint = 0
3536 else : self.fint = 0
3136 self.index += 1
3537 self.index += 1
3137 if numpy.any(self.buffer):
3538 if numpy.any(self.buffer):
@@ -3146,7 +3547,13 class SpectralFitting(Operation):
3146 return dataOut
3547 return dataOut
3147 if path != None:
3548 if path != None:
3148 sys.path.append(path)
3549 sys.path.append(path)
3149 self.library = importlib.import_module(file)
3550 try:
3551 self.library = importlib.import_module(file)
3552 except:
3553 pass
3554 if filec != None:
3555 self.weightf = importlib.import_module(filec)
3556 #self.weightf = importlib.import_module('weightfit')
3150
3557
3151 #To be inserted as a parameter
3558 #To be inserted as a parameter
3152 groupArray = numpy.array(groupList)
3559 groupArray = numpy.array(groupList)
@@ -3162,8 +3569,11 class SpectralFitting(Operation):
3162 dataOut.data_paramC = None
3569 dataOut.data_paramC = None
3163
3570
3164 #Set constants
3571 #Set constants
3165 constants = self.library.setConstants(dataOut)
3572 try:
3166 dataOut.constants = constants
3573 constants = self.library.setConstants(dataOut)
3574 dataOut.constants = constants
3575 except:
3576 pass
3167 M = dataOut.normFactor
3577 M = dataOut.normFactor
3168 N = dataOut.nFFTPoints
3578 N = dataOut.nFFTPoints
3169 ippSeconds = dataOut.ippSeconds
3579 ippSeconds = dataOut.ippSeconds
@@ -3190,210 +3600,300 class SpectralFitting(Operation):
3190 if not self.isConfig:
3600 if not self.isConfig:
3191 self.isConfig = True
3601 self.isConfig = True
3192
3602
3193 index = tini.tm_hour*12+tini.tm_min/5
3603 index = tini.tm_hour*12+tini.tm_min/taver
3194 jspc = jspc/N/N
3604 jspc = jspc/N/N
3195 jcspc = jcspc/N/N
3605 jcspc = jcspc/N/N
3196 tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.CleanRayleigh(dataOut,jspc,jcspc,2)
3606 tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.CleanRayleigh(dataOut,jspc,jcspc,2)
3197 jspectra = tmp_spectra*len(jspc[:,0,0,0])
3607 jspectra = tmp_spectra*len(jspc[:,0,0,0])
3198 jcspectra = tmp_cspectra*len(jspc[:,0,0,0])
3608 jcspectra = tmp_cspectra*len(jspc[:,0,0,0])
3199 my_incoh_spectra ,my_incoh_cspectra,my_incoh_aver,my_coh_aver, incoh_spectra, coh_spectra, incoh_cspectra, coh_cspectra, incoh_aver, coh_aver = self.__DiffCoherent(jspectra, jcspectra, dataOut, noise, snrth, None, None)
3609 my_incoh_spectra ,my_incoh_cspectra,my_incoh_aver,my_coh_aver, incoh_spectra, coh_spectra, incoh_cspectra, coh_cspectra, incoh_aver, coh_aver = self.__DiffCoherent(jspectra, jcspectra, dataOut, noise, snrth,coh_th, hei_th)
3200 clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index)
3610 clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index)
3201 dataOut.data_spc = incoh_spectra
3611 dataOut.data_spc = incoh_spectra
3202 dataOut.data_cspc = incoh_cspectra
3612 dataOut.data_cspc = incoh_cspectra
3613 #dataOut.data_spc = tmp_spectra
3614 #dataOut.data_cspc = tmp_cspectra
3203
3615
3204 clean_num_aver = incoh_aver*len(jspc[:,0,0,0])
3616 clean_num_aver = incoh_aver*len(jspc[:,0,0,0])
3205 coh_num_aver = clean_coh_aver*len(jspc[:,0,0,0])
3617 coh_num_aver = clean_coh_aver*len(jspc[:,0,0,0])
3618 #clean_num_aver = (numpy.zeros([nChan, nHei])+1)*len(jspc[:,0,0,0])
3619 #coh_num_aver = numpy.zeros([nChan, nHei])*0*len(jspc[:,0,0,0])
3206 #List of possible combinations
3620 #List of possible combinations
3207 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
3621 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
3208 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
3622 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
3209
3623 if Gaussian_Windowed == 1:
3210 if getSNR:
3624 #dataOut.data_spc = jspectra
3211 listChannels = groupArray.reshape((groupArray.size))
3625 '''
3212 listChannels.sort()
3626 Written by R. Flores
3213 dataOut.data_SNR = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise[listChannels])
3627 '''
3214 if dataOut.data_paramC is None:
3628 print("normFactor: ", dataOut.normFactor)
3215 dataOut.data_paramC = numpy.zeros((nGroups*4, nHeights,2))*numpy.nan
3629 data_spc_aux = numpy.copy(dataOut.data_spc)#[:,0,:]
3216 for i in range(nGroups):
3630 data_spc_aux[:,0,:] = (data_spc_aux[:,1,:]+data_spc_aux[:,-1,:])/2
3217 coord = groupArray[i,:]
3631 #'''
3218 #Input data array
3632 from scipy.signal import medfilt
3219 data = dataOut.data_spc[coord,:,:]/(M*N)
3633 import matplotlib.pyplot as plt
3220 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
3634 dataOut.moments = numpy.ones((dataOut.nChannels,4,dataOut.nHeights))*numpy.NAN
3221
3635 dataOut.VelRange = dataOut.getVelRange(0)
3222 #Cross Spectra data array for Covariance Matrixes
3636 for nChannel in range(dataOut.nChannels):
3223 ind = 0
3637 for hei in range(dataOut.heightList.shape[0]):
3224 for pairs in listComb:
3638 #print("ipp: ", dataOut.ippSeconds)
3225 pairsSel = numpy.array([coord[x],coord[y]])
3639 #spc = numpy.copy(dataOut.data_spc[nChannel,:,hei])
3226 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
3640 spc = data_spc_aux[nChannel,:,hei]
3227 ind += 1
3641 if spc.all() == 0.:
3228 dataCross = dataOut.data_cspc[indCross,:,:]/(M*N)
3642 print("CONTINUE")
3229 dataCross = dataCross**2
3643 continue
3230 nhei = nHeights
3644 #print(VelRange)
3231 poweri = numpy.sum(dataOut.data_spc[:,1:nProf-0,:],axis=1)/clean_num_aver[:,:]
3645 #print(dataOut.getFreqRange(64))
3232 if i == 0 : my_noises = numpy.zeros(4,dtype=float) #FLTARR(4)
3646 #print("Hei: ", dataOut.heightList[hei])
3233 n0i = numpy.nanmin(poweri[0+i*2,0:nhei-0])/(nProf-1)
3647
3234 n1i = numpy.nanmin(poweri[1+i*2,0:nhei-0])/(nProf-1)
3648 spc_mod = numpy.copy(spc)
3235 n0 = n0i
3649 spcm = medfilt(spc_mod,11)
3236 n1= n1i
3650 spc_max = numpy.max(spcm)
3237 my_noises[2*i+0] = n0
3651 dop1_x0 = dataOut.VelRange[numpy.argmax(spcm)]
3238 my_noises[2*i+1] = n1
3652 #D = numpy.min(spcm)
3239 snrth = -16.0
3653 D_in = (numpy.mean(spcm[:15])+numpy.mean(spcm[-15:]))/2.
3240 snrth = 10**(snrth/10.0)
3654 #print("spc_max: ", spc_max)
3241
3655 #print("dataOut.ippSeconds: ", dataOut.ippSeconds, dataOut.timeInterval)
3242 for h in range(nHeights):
3656 ##fun, A, B, C, D = self.windowing_single(spc,dataOut.VelRange,spc_max,dop1_x0,abs(dop1_x0),D,dataOut.nFFTPoints)
3243 d = data[:,h]
3657 #fun, A, B, C, D = self.windowing_single(spc,dataOut.VelRange,spc_max,dop1_x0,abs(dop1_x0),D,dataOut.nFFTPoints)
3244 smooth = clean_num_aver[i+1,h] #dataOut.data_spc[:,1:nProf-0,:]
3658 fun, A, B, C, D = self.windowing_single_direct(spc_mod,dataOut.VelRange,spc_max,dop1_x0,abs(dop1_x0/5),D_in,dataOut.nFFTPoints,dataOut.timeInterval)
3245 signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth
3659
3246 signalpn1 = (dataOut.data_spc[i*2+1,1:(nProf-0),h])/smooth
3660 dataOut.moments[nChannel,0,hei] = A
3247 signal0 = signalpn0-n0
3661 dataOut.moments[nChannel,1,hei] = B
3248 signal1 = signalpn1-n1
3662 dataOut.moments[nChannel,2,hei] = C
3249 snr0 = numpy.sum(signal0/n0)/(nProf-1)
3663 dataOut.moments[nChannel,3,hei] = D
3250 snr1 = numpy.sum(signal1/n1)/(nProf-1)
3664 '''
3251 if snr0 > snrth and snr1 > snrth and clean_num_aver[i+1,h] > 0 :
3665 if nChannel == 0:
3252 #Covariance Matrix
3666 print(dataOut.heightList[hei])
3253 D = numpy.diag(d**2)
3667 plt.figure()
3254 ind = 0
3668 plt.plot(dataOut.VelRange,spc,marker='*',linestyle='--')
3255 for pairs in listComb:
3669 plt.plot(dataOut.VelRange,fun)
3256 #Coordinates in Covariance Matrix
3670 plt.title(dataOut.heightList[hei])
3257 x = pairs[0]
3671 plt.show()
3258 y = pairs[1]
3672 '''
3259 #Channel Index
3673 #plt.show()
3260 S12 = dataCross[ind,:,h]
3674 #'''
3261 D12 = numpy.diag(S12)
3675 dataOut.data_spc = jspectra
3262 #Completing Covariance Matrix with Cross Spectras
3676 print("SUCCESS")
3263 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
3677 return dataOut
3264 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
3678
3265 ind += 1
3679 elif Gaussian_Windowed == 2: #Only to clean spc
3266 diagD = numpy.zeros(256)
3680 dataOut.VelRange = dataOut.getVelRange(0)
3267 if h == 17 :
3681 return dataOut
3268 for ii in range(256): diagD[ii] = D[ii,ii]
3682 else:
3269 #Dinv=numpy.linalg.inv(D)
3683 if getSNR:
3270 #L=numpy.linalg.cholesky(Dinv)
3684 listChannels = groupArray.reshape((groupArray.size))
3271 try:
3685 listChannels.sort()
3272 Dinv=numpy.linalg.inv(D)
3686 dataOut.data_SNR = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise[listChannels])
3273 L=numpy.linalg.cholesky(Dinv)
3687 if dataOut.data_paramC is None:
3274 except:
3688 dataOut.data_paramC = numpy.zeros((nGroups*4, nHeights,2))*numpy.nan
3275 Dinv = D*numpy.nan
3689 for i in range(nGroups):
3276 L= D*numpy.nan
3690 coord = groupArray[i,:]
3277 LT=L.T
3691 #Input data array
3278
3692 data = dataOut.data_spc[coord,:,:]/(M*N)
3279 dp = numpy.dot(LT,d)
3693 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
3280
3694
3281 #Initial values
3695 #Cross Spectra data array for Covariance Matrixes
3282 data_spc = dataOut.data_spc[coord,:,h]
3696 ind = 0
3283
3697 for pairs in listComb:
3284 if (h>0)and(error1[3]<5):
3698 pairsSel = numpy.array([coord[x],coord[y]])
3285 p0 = dataOut.data_param[i,:,h-1]
3699 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
3286 else:
3700 ind += 1
3287 p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants))# sin el i(data_spc, constants, i)
3701 dataCross = dataOut.data_cspc[indCross,:,:]/(M*N)
3288 try:
3702 dataCross = dataCross**2
3289 #Least Squares
3703 nhei = nHeights
3290 #print (dp,LT,constants)
3704 poweri = numpy.sum(dataOut.data_spc[:,1:nProf-0,:],axis=1)/clean_num_aver[:,:]
3291 #value =self.__residFunction(p0,dp,LT,constants)
3705 if i == 0 : my_noises = numpy.zeros(4,dtype=float)
3292 #print ("valueREADY",value.shape, type(value))
3706 n0i = numpy.nanmin(poweri[0+i*2,0:nhei-0])/(nProf-1)
3293 #optimize.leastsq(value)
3707 n1i = numpy.nanmin(poweri[1+i*2,0:nhei-0])/(nProf-1)
3294 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
3708 n0 = n0i
3295 #minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
3709 n1= n1i
3296 #Chi square error
3710 my_noises[2*i+0] = n0
3297 #print(minp,covp.infodict,mesg,ier)
3711 my_noises[2*i+1] = n1
3298 #print("REALIZA OPTIMIZ")
3712 snrth = -25.0 # -4
3299 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
3713 snrth = 10**(snrth/10.0)
3300 #Error with Jacobian
3714 jvelr = numpy.zeros(nHeights, dtype = 'float')
3301 error1 = self.library.errorFunction(minp,constants,LT)
3715 #snr0 = numpy.zeros(nHeights, dtype = 'float')
3302 # print self.__residFunction(p0,dp,LT, constants)
3716 #snr1 = numpy.zeros(nHeights, dtype = 'float')
3303 # print infodict['fvec']
3717 hvalid = [0]
3304 # print self.__residFunction(minp,dp,LT,constants)
3718
3305
3719 coh2 = abs(dataOut.data_cspc[i,1:nProf,:])**2/(dataOut.data_spc[0+i*2,1:nProf-0,:]*dataOut.data_spc[1+i*2,1:nProf-0,:])
3306 except:
3720
3721 for h in range(nHeights):
3722 smooth = clean_num_aver[i+1,h]
3723 signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth
3724 signalpn1 = (dataOut.data_spc[i*2+1,1:(nProf-0),h])/smooth
3725 signal0 = signalpn0-n0
3726 signal1 = signalpn1-n1
3727 snr0 = numpy.sum(signal0/n0)/(nProf-1)
3728 snr1 = numpy.sum(signal1/n1)/(nProf-1)
3729 #jmax0 = MAX(signal0,maxp0)
3730 #jmax1 = MAX(signal1,maxp1)
3731 gamma = coh2[:,h]
3732
3733 indxs = (numpy.isfinite(list(gamma))==True).nonzero()
3734
3735 if len(indxs) >0:
3736 if numpy.nanmean(gamma) > 0.07:
3737 maxp0 = numpy.argmax(signal0*gamma)
3738 maxp1 = numpy.argmax(signal1*gamma)
3739 #print('usa gamma',numpy.nanmean(gamma))
3740 else:
3741 maxp0 = numpy.argmax(signal0)
3742 maxp1 = numpy.argmax(signal1)
3743 jvelr[h] = (absc[maxp0]+absc[maxp1])/2.
3744 else: jvelr[h] = absc[0]
3745 if snr0 > 0.1 and snr1 > 0.1: hvalid = numpy.concatenate((hvalid,h), axis=None)
3746 #print(maxp0,absc[maxp0],snr0,jvelr[h])
3747
3748 if len(hvalid)> 1: fd0 = numpy.median(jvelr[hvalid[1:]])*-1
3749 else: fd0 = numpy.nan
3750 #print(fd0,hvalid)
3751 for h in range(nHeights):
3752 d = data[:,h]
3753 smooth = clean_num_aver[i+1,h] #dataOut.data_spc[:,1:nProf-0,:]
3754 signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth
3755 signalpn1 = (dataOut.data_spc[i*2+1,1:(nProf-0),h])/smooth
3756 signal0 = signalpn0-n0
3757 signal1 = signalpn1-n1
3758 snr0 = numpy.sum(signal0/n0)/(nProf-1)
3759 snr1 = numpy.sum(signal1/n1)/(nProf-1)
3760
3761 if snr0 > snrth and snr1 > snrth and clean_num_aver[i+1,h] > 0 :
3762 #Covariance Matrix
3763 D = numpy.diag(d**2)
3764 ind = 0
3765 for pairs in listComb:
3766 #Coordinates in Covariance Matrix
3767 x = pairs[0]
3768 y = pairs[1]
3769 #Channel Index
3770 S12 = dataCross[ind,:,h]
3771 D12 = numpy.diag(S12)
3772 #Completing Covariance Matrix with Cross Spectras
3773 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
3774 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
3775 ind += 1
3776 diagD = numpy.zeros(256)
3777
3778 #Dinv=numpy.linalg.inv(D)
3779 #L=numpy.linalg.cholesky(Dinv)
3780 try:
3781 Dinv=numpy.linalg.inv(D)
3782 L=numpy.linalg.cholesky(Dinv)
3783 except:
3784 Dinv = D*numpy.nan
3785 L= D*numpy.nan
3786 LT=L.T
3787
3788 dp = numpy.dot(LT,d)
3789
3790 #Initial values
3791 data_spc = dataOut.data_spc[coord,:,h]
3792 w = data_spc/data_spc
3793 if filec != None:
3794 w = self.weightf.weightfit(w,tini.tm_year,tini.tm_yday,index,h,i)
3795
3796 if (h>6) and (error1[3]<25):
3797 p0 = dataOut.data_param[i,:,h-1]
3798 #print('usa anterior')
3799 else:
3800 p0 = numpy.array(self.library.initialValuesFunction(data_spc*w, constants))# sin el i(data_spc, constants, i)
3801
3802 if filec != None:
3803 p0 = self.weightf.Vrfit(p0,tini.tm_year,tini.tm_yday,index,h,i)
3804 p0[3] = fd0
3805 #if index == 175 and i==1 and h>=27 and h<=35: p0[3]=30
3806 #if h >= 6 and i==1 and h<= 10: print(p0)
3807 try:
3808 #Least Squares
3809 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
3810 #minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
3811 #Chi square error
3812 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
3813 #Error with Jacobian
3814 error1 = self.library.errorFunction(minp,constants,LT)
3815 #if h >= 0 and h<= 10 and i ==0: print(p0,minp,error1)
3816 #if i>=0 and h>=0: print(index,h,minp[3])
3817 # print self.__residFunction(p0,dp,LT, constants)
3818 # print infodict['fvec']
3819 # print self.__residFunction(minp,dp,LT,constants)
3820
3821 except:
3822 minp = p0*numpy.nan
3823 error0 = numpy.nan
3824 error1 = p0*numpy.nan
3825 # s_sq = (self.__residFunction(minp,dp,LT,constants)).sum()/(len(dp)-len(p0))
3826 # covp = covp*s_sq
3827 # error = []
3828 # for ip in range(len(minp)):
3829 # try:
3830 # error.append(numpy.absolute(covp[ip][ip])**0.5)
3831 # except:
3832 # error.append( 0.00 )
3833 #if i==1 and h==11 and index == 139: print(p0, minp,data_spc)
3834 else :
3835 data_spc = dataOut.data_spc[coord,:,h]
3836 p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants))
3307 minp = p0*numpy.nan
3837 minp = p0*numpy.nan
3308 error0 = numpy.nan
3838 error0 = numpy.nan
3309 error1 = p0*numpy.nan
3839 error1 = p0*numpy.nan
3310 #print ("EXCEPT 0000000000")
3840
3311 # s_sq = (self.__residFunction(minp,dp,LT,constants)).sum()/(len(dp)-len(p0))
3841 if dataOut.data_param is None:
3312 # covp = covp*s_sq
3842 dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
3313 # #print("TRY___________________________________________1")
3843 dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
3314 # error = []
3844
3315 # for ip in range(len(minp)):
3845 dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
3316 # try:
3846 dataOut.data_param[i,:,h] = minp
3317 # error.append(numpy.absolute(covp[ip][ip])**0.5)
3847
3318 # except:
3848 for ht in range(nHeights-1) :
3319 # error.append( 0.00 )
3849 smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam]
3320 else :
3850 dataOut.data_paramC[4*i,ht,1] = smooth
3321 data_spc = dataOut.data_spc[coord,:,h]
3851 signalpn0 = (coh_spectra[i*2 ,1:(nProf-0),ht])/smooth #coh_spectra
3322 p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants))
3852 signalpn1 = (coh_spectra[i*2+1,1:(nProf-0),ht])/smooth
3323 minp = p0*numpy.nan
3853
3324 error0 = numpy.nan
3854 val0 = (signalpn0 > 0).nonzero()
3325 error1 = p0*numpy.nan
3855 val0 = val0[0]
3326 #Save
3856
3327 if dataOut.data_param is None:
3857 if len(val0) == 0 : val0_npoints = nProf
3328 dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
3858 else : val0_npoints = len(val0)
3329 dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
3859
3330
3860 val1 = (signalpn1 > 0).nonzero()
3331 dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
3861 val1 = val1[0]
3332 dataOut.data_param[i,:,h] = minp
3862 if len(val1) == 0 : val1_npoints = nProf
3333
3863 else : val1_npoints = len(val1)
3334 for ht in range(nHeights-1) :
3864
3335 smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam]
3865 dataOut.data_paramC[0+4*i,ht,0] = numpy.sum((signalpn0/val0_npoints))/n0
3336 dataOut.data_paramC[4*i,ht,1] = smooth
3866 dataOut.data_paramC[1+4*i,ht,0] = numpy.sum((signalpn1/val1_npoints))/n1
3337 signalpn0 = (coh_spectra[i*2 ,1:(nProf-0),ht])/smooth #coh_spectra
3867
3338 signalpn1 = (coh_spectra[i*2+1,1:(nProf-0),ht])/smooth
3868 signal0 = (signalpn0-n0)
3339
3869 vali = (signal0 < 0).nonzero()
3340 #val0 = WHERE(signalpn0 > 0,cval0)
3870 vali = vali[0]
3341 val0 = (signalpn0 > 0).nonzero()
3871 if len(vali) > 0 : signal0[vali] = 0
3342 val0 = val0[0]
3872 signal1 = (signalpn1-n1)
3343 #print('hvalid:',hvalid)
3873 vali = (signal1 < 0).nonzero()
3344 #print('valid', valid)
3874 vali = vali[0]
3345 if len(val0) == 0 : val0_npoints = nProf
3875 if len(vali) > 0 : signal1[vali] = 0
3346 else : val0_npoints = len(val0)
3876 snr0 = numpy.sum(signal0/n0)/(nProf-1)
3347
3877 snr1 = numpy.sum(signal1/n1)/(nProf-1)
3348 #val1 = WHERE(signalpn1 > 0,cval1)
3878 doppler = absc[1:]
3349 val1 = (signalpn1 > 0).nonzero()
3879 if snr0 >= snrth and snr1 >= snrth and smooth :
3350 val1 = val1[0]
3880 signalpn0_n0 = signalpn0
3351 if len(val1) == 0 : val1_npoints = nProf
3881 signalpn0_n0[val0] = signalpn0[val0] - n0
3352 else : val1_npoints = len(val1)
3882 mom0 = self.moments(doppler,signalpn0-n0,nProf)
3353
3883
3354 dataOut.data_paramC[0+4*i,ht,0] = numpy.sum((signalpn0/val0_npoints))/n0
3884 signalpn1_n1 = signalpn1
3355 dataOut.data_paramC[1+4*i,ht,0] = numpy.sum((signalpn1/val1_npoints))/n1
3885 signalpn1_n1[val1] = signalpn1[val1] - n1
3356
3886 mom1 = self.moments(doppler,signalpn1_n1,nProf)
3357 signal0 = (signalpn0-n0) # > 0
3887 dataOut.data_paramC[2+4*i,ht,0] = (mom0[0]+mom1[0])/2.
3358 vali = (signal0 < 0).nonzero()
3888 dataOut.data_paramC[3+4*i,ht,0] = (mom0[1]+mom1[1])/2.
3359 vali = vali[0]
3889
3360 if len(vali) > 0 : signal0[vali] = 0
3890 dataOut.data_spc = jspectra
3361 signal1 = (signalpn1-n1) #> 0
3891 if getSNR:
3362 vali = (signal1 < 0).nonzero()
3892 listChannels = groupArray.reshape((groupArray.size))
3363 vali = vali[0]
3893 listChannels.sort()
3364 if len(vali) > 0 : signal1[vali] = 0
3894
3365 snr0 = numpy.sum(signal0/n0)/(nProf-1)
3895 dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:], my_noises[listChannels])
3366 snr1 = numpy.sum(signal1/n1)/(nProf-1)
3896 return dataOut
3367 doppler = absc[1:]
3368 if snr0 >= snrth and snr1 >= snrth and smooth :
3369 signalpn0_n0 = signalpn0
3370 signalpn0_n0[val0] = signalpn0[val0] - n0
3371 mom0 = self.moments(doppler,signalpn0-n0,nProf)
3372 # sigtmp= numpy.transpose(numpy.tile(signalpn0, [4,1]))
3373 # momt= self.__calculateMoments( sigtmp, doppler , n0 )
3374 signalpn1_n1 = signalpn1
3375 signalpn1_n1[val1] = signalpn1[val1] - n1
3376 mom1 = self.moments(doppler,signalpn1_n1,nProf)
3377 dataOut.data_paramC[2+4*i,ht,0] = (mom0[0]+mom1[0])/2.
3378 dataOut.data_paramC[3+4*i,ht,0] = (mom0[1]+mom1[1])/2.
3379 # if graph == 1 :
3380 # window, 13
3381 # plot,doppler,signalpn0
3382 # oplot,doppler,signalpn1,linest=1
3383 # oplot,mom0(0)*doppler/doppler,signalpn0
3384 # oplot,mom1(0)*doppler/doppler,signalpn1
3385 # print,interval/12.,beam,45+ht*15,snr0,snr1,mom0(0),mom1(0),mom0(1),mom1(1)
3386 #ENDIF
3387 #ENDIF
3388 #ENDFOR End height
3389
3390 dataOut.data_spc = jspectra
3391 if getSNR:
3392 listChannels = groupArray.reshape((groupArray.size))
3393 listChannels.sort()
3394
3395 dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:], my_noises[listChannels])
3396 return dataOut
3397
3897
3398 def __residFunction(self, p, dp, LT, constants):
3898 def __residFunction(self, p, dp, LT, constants):
3399
3899
@@ -6628,9 +7128,9 class IGRFModel(Operation):
6628 mkfact_short_2020_2.mkfact(dataOut.year,dataOut.h,dataOut.bfm,dataOut.thb,dataOut.bki,dataOut.MAXNRANGENDT)
7128 mkfact_short_2020_2.mkfact(dataOut.year,dataOut.h,dataOut.bfm,dataOut.thb,dataOut.bki,dataOut.MAXNRANGENDT)
6629
7129
6630 #mkfact_short_2020.mkfact(dataOut.year,dataOut.h,dataOut.bfm,dataOut.thb,dataOut.bki,dataOut.MAXNRANGENDT)
7130 #mkfact_short_2020.mkfact(dataOut.year,dataOut.h,dataOut.bfm,dataOut.thb,dataOut.bki,dataOut.MAXNRANGENDT)
6631 print("bki: ", dataOut.bki[:10])
7131 #print("bki: ", dataOut.bki[:10])
6632 print("thb: ", dataOut.thb[:10])
7132 #print("thb: ", dataOut.thb[:10])
6633 print("bfm: ", dataOut.bfm[:10])
7133 #print("bfm: ", dataOut.bfm[:10])
6634 #print("IDs: ", id(dataOut.bki))
7134 #print("IDs: ", id(dataOut.bki))
6635
7135
6636 return dataOut
7136 return dataOut
@@ -6794,8 +7294,10 class MergeProc(ProcessingUnit):
6794 self.dataOut.NSHTS = int(numpy.shape(self.dataOut.ph2)[0])
7294 self.dataOut.NSHTS = int(numpy.shape(self.dataOut.ph2)[0])
6795 dH = self.dataOut.heightList[1]-self.dataOut.heightList[0]
7295 dH = self.dataOut.heightList[1]-self.dataOut.heightList[0]
6796 dH /= self.dataOut.windowOfFilter
7296 dH /= self.dataOut.windowOfFilter
6797 self.dataOut.heightList = numpy.arange(0,self.dataOut.NSHTS)*dH
7297 self.dataOut.heightList = numpy.arange(0,self.dataOut.NSHTS)*dH + dH
7298 #print("heightList: ", self.dataOut.heightList)
6798 self.dataOut.NDP = self.dataOut.NSHTS
7299 self.dataOut.NDP = self.dataOut.NSHTS
7300 #exit(1)
6799 #print(self.dataOut.heightList)
7301 #print(self.dataOut.heightList)
6800
7302
6801 class MST_Den_Conv(Operation):
7303 class MST_Den_Conv(Operation):
Show file before | Show file after | Hide comments
@@ -455,19 +455,43 class GetSNR(Operation):
455
455
456 def run(self,dataOut):
456 def run(self,dataOut):
457
457
458 noise = dataOut.getNoise()
458 #noise = dataOut.getNoise()
459 maxdB = 16
459 noise = dataOut.getNoise(ymin_index=-10) #Región superior donde solo debería de haber ruido
460
460 #print("Noise: ", noise)
461 normFactor = 24
461 #print("Noise_dB: ", 10*numpy.log10(noise/dataOut.normFactor))
462
462 #print("Heights: ", dataOut.heightList)
463 #dataOut.data_snr = (dataOut.data_spc.sum(axis=1))/(noise[:,None]*dataOut.normFactor)
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)
464 ################dataOut.data_snr = (dataOut.data_spc.sum(axis=1))/(noise[:,None]*dataOut.nFFTPoints) #Before 12Jan2023
465
465 #dataOut.data_snr = (dataOut.data_spc.sum(axis=1)-noise[:,None])/(noise[:,None])
466 dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.5, numpy.nan, dataOut.data_snr)
466 dataOut.data_snr = (dataOut.data_spc.sum(axis=1)-noise[:,None]*dataOut.nFFTPoints)/(noise[:,None]*dataOut.nFFTPoints) #It works apparently
467 dataOut.snl = numpy.log10(dataOut.data_snr)
468 #print("snl: ", dataOut.snl)
469 #exit(1)
470 #print(dataOut.heightList[-11])
471 #print(numpy.shape(dataOut.heightList))
472 #print(dataOut.data_snr)
473 #print(dataOut.data_snr[0,-11])
474 #exit(1)
475 #dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.5, numpy.nan, dataOut.data_snr)
476 #dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.1, numpy.nan, dataOut.data_snr)
477 #dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.0, numpy.nan, dataOut.data_snr)
478 #dataOut.data_snr = numpy.where(dataOut.data_snr<.05, numpy.nan, dataOut.data_snr)
479 dataOut.snl = numpy.where(dataOut.data_snr<.01, numpy.nan, dataOut.snl)
480 '''
481 import matplotlib.pyplot as plt
482 #plt.plot(10*numpy.log10(dataOut.data_snr[0]),dataOut.heightList)
483 plt.plot(dataOut.data_snr[0],dataOut.heightList)#,marker='*')
484 plt.xlim(-1,10)
485 plt.axvline(1,color='k')
486 plt.axvline(.1,color='k',linestyle='--')
487 plt.grid()
488 plt.show()
489 '''
467 #dataOut.data_snr = 10*numpy.log10(dataOut.data_snr)
490 #dataOut.data_snr = 10*numpy.log10(dataOut.data_snr)
468 #dataOut.data_snr = numpy.expand_dims(dataOut.data_snr,axis=0)
491 #dataOut.data_snr = numpy.expand_dims(dataOut.data_snr,axis=0)
469 #print(dataOut.data_snr.shape)
492 #print(dataOut.data_snr.shape)
470 #exit(1)
493 #exit(1)
494 #print("Before: ", dataOut.data_snr[0])
471
495
472
496
473 return dataOut
497 return dataOut
Show file before | Show file after | Hide comments
@@ -3815,6 +3815,9 class IncohInt(Operation):
3815 dataOut.flagNoData = False
3815 dataOut.flagNoData = False
3816
3816
3817 dataOut.VelRange = dataOut.getVelRange(0)
3817 dataOut.VelRange = dataOut.getVelRange(0)
3818 dataOut.FreqRange = dataOut.getFreqRange(0)/1000.
3819 #print("VelRange: ", dataOut.VelRange)
3820 #exit(1)
3818
3821
3819 #print("Power",numpy.sum(dataOut.data_spc[0,:,20:30],axis=0))
3822 #print("Power",numpy.sum(dataOut.data_spc[0,:,20:30],axis=0))
3820 #print("Power",numpy.sum(dataOut.data_spc[0,100:110,:],axis=1))
3823 #print("Power",numpy.sum(dataOut.data_spc[0,100:110,:],axis=1))
@@ -3921,7 +3924,7 class SnrFaraday(Operation):
3921
3924
3922 return dataOut
3925 return dataOut
3923
3926
3924 class SpectraDataToFaraday_07_11_2022(Operation):
3927 class SpectraDataToFaraday(Operation): #ISR MODE
3925 '''
3928 '''
3926 Written by R. Flores
3929 Written by R. Flores
3927 '''
3930 '''
@@ -4144,10 +4147,14 class SpectraDataToFaraday_07_11_2022(Operation):
4144 #print(data_eej)
4147 #print(data_eej)
4145 index_eej = CleanCohEchoes.mad_based_outlier(self,data_eej[:17])
4148 index_eej = CleanCohEchoes.mad_based_outlier(self,data_eej[:17])
4146 aux_eej = numpy.array(index_eej.nonzero()).ravel()
4149 aux_eej = numpy.array(index_eej.nonzero()).ravel()
4150 print("aux_eej: ", aux_eej)
4151 if aux_eej != []:
4152 dataOut.min_id_eej = aux_eej[-1]
4153 else:
4154 dataOut.min_id_eej = 12
4147
4155
4148 dataOut.min_id_eej = aux_eej[-1]
4149
4156
4150 print(dataOut.min_id_eej)
4157 #print("min_id_eej: ", dataOut.min_id_eej)
4151 #exit(1)
4158 #exit(1)
4152
4159
4153 def run(self,dataOut):
4160 def run(self,dataOut):
@@ -4204,7 +4211,7 class SpectraDataToFaraday_07_11_2022(Operation):
4204 #exit(1)
4211 #exit(1)
4205 return dataOut
4212 return dataOut
4206
4213
4207 class SpectraDataToFaraday(Operation):
4214 class SpectraDataToFaraday_MST(Operation): #MST MODE
4208 """Operation to use spectra data in Faraday processing.
4215 """Operation to use spectra data in Faraday processing.
4209
4216
4210 Parameters:
4217 Parameters:
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (1044 lines changed) Show them Hide them
@@ -38,7 +38,8 class VoltageProc(ProcessingUnit):
38 if self.dataIn.type == 'Voltage':
38 if self.dataIn.type == 'Voltage':
39 self.dataOut.copy(self.dataIn)
39 self.dataOut.copy(self.dataIn)
40 self.dataOut.runNextUnit = runNextUnit
40 self.dataOut.runNextUnit = runNextUnit
41 #print(self.dataOut.data.shape)
41 #print("data shape: ", self.dataOut.data.shape)
42 #exit(1)
42
43
43 def __updateObjFromAmisrInput(self):
44 def __updateObjFromAmisrInput(self):
44
45
@@ -226,6 +227,7 class selectHeights(Operation):
226 #print(dataOut.nHeights)
227 #print(dataOut.nHeights)
227 #exit(1)
228 #exit(1)
228 '''
229 '''
230 #print("select Heights: ", self.dataOut.heightList)
229 return self.dataOut
231 return self.dataOut
230
232
231 def selectHeightsByIndex(self, minIndex, maxIndex):
233 def selectHeightsByIndex(self, minIndex, maxIndex):
@@ -1250,7 +1252,7 class LagsReshapeLP(Operation):
1250 self.buffer[3,:,:,:] = numpy.transpose(numpy.conj(self.lagp3),(1,2,0))
1252 self.buffer[3,:,:,:] = numpy.transpose(numpy.conj(self.lagp3),(1,2,0))
1251
1253
1252 #print(self.buffer[3,:,100,15])
1254 #print(self.buffer[3,:,100,15])
1253 print(sum(self.buffer[3,:,199,2]))
1255 print("Sum: ", sum(self.buffer[3,:,199,2]))
1254 #print(self.cnorm)
1256 #print(self.cnorm)
1255 exit(1)
1257 exit(1)
1256
1258
@@ -3348,7 +3350,7 class DoublePulseACFs_PerLag(Operation):
3348 exit(1)
3350 exit(1)
3349 '''
3351 '''
3350 #print(pa)
3352 #print(pa)
3351 #'''
3353 '''
3352 import matplotlib.pyplot as plt
3354 import matplotlib.pyplot as plt
3353 #plt.plot(dataOut.p[:,-1],dataOut.heightList)
3355 #plt.plot(dataOut.p[:,-1],dataOut.heightList)
3354 plt.plot(pa/dataOut.pan-1.,dataOut.heightList)
3356 plt.plot(pa/dataOut.pan-1.,dataOut.heightList)
@@ -3358,7 +3360,7 class DoublePulseACFs_PerLag(Operation):
3358 plt.show()
3360 plt.show()
3359 #print("p: ",dataOut.p[33,:])
3361 #print("p: ",dataOut.p[33,:])
3360 #exit(1)
3362 #exit(1)
3361 #'''
3363 '''
3362 return dataOut
3364 return dataOut
3363
3365
3364 class FaradayAngleAndDPPower(Operation):
3366 class FaradayAngleAndDPPower(Operation):
@@ -3399,7 +3401,7 class FaradayAngleAndDPPower(Operation):
3399 for i in range(dataOut.MAXNRANGENDT):
3401 for i in range(dataOut.MAXNRANGENDT):
3400 dataOut.range1[i]=dataOut.H0 + i*dataOut.DH
3402 dataOut.range1[i]=dataOut.H0 + i*dataOut.DH
3401 dataOut.h2[i]=dataOut.range1[i]**2
3403 dataOut.h2[i]=dataOut.range1[i]**2
3402 print("shape ph2",numpy.shape(dataOut.ph2))
3404 #print("shape ph2",numpy.shape(dataOut.ph2))
3403 for j in range(dataOut.NDP):
3405 for j in range(dataOut.NDP):
3404 dataOut.ph2[j]=0.
3406 dataOut.ph2[j]=0.
3405 dataOut.sdp2[j]=0.
3407 dataOut.sdp2[j]=0.
@@ -3505,7 +3507,7 class ElectronDensityFaraday(Operation):
3505 ndphi=dataOut.NSHTS-4
3507 ndphi=dataOut.NSHTS-4
3506 #print(dataOut.phi)
3508 #print(dataOut.phi)
3507 #exit(1)
3509 #exit(1)
3508 '''
3510 #'''
3509 if hasattr(dataOut, 'flagSpreadF') and dataOut.flagSpreadF:
3511 if hasattr(dataOut, 'flagSpreadF') and dataOut.flagSpreadF:
3510 #if dataOut.flagSpreadF:
3512 #if dataOut.flagSpreadF:
3511 nanindex = numpy.argwhere(numpy.isnan(dataOut.phi))
3513 nanindex = numpy.argwhere(numpy.isnan(dataOut.phi))
@@ -3517,7 +3519,7 class ElectronDensityFaraday(Operation):
3517 else:
3519 else:
3518 #dataOut.phi_uwrp = dataOut.phi.copy()
3520 #dataOut.phi_uwrp = dataOut.phi.copy()
3519 dataOut.phi[:]=numpy.unwrap(dataOut.phi[:]) #Better results
3521 dataOut.phi[:]=numpy.unwrap(dataOut.phi[:]) #Better results
3520 '''
3522 #'''
3521 #print(dataOut.phi)
3523 #print(dataOut.phi)
3522 #print(dataOut.ph2)
3524 #print(dataOut.ph2)
3523 #exit(1)
3525 #exit(1)
@@ -3968,12 +3970,12 class NormalizeDPPowerRoberto_V2(Operation):
3968 if hasattr(dataOut, 'flagSpreadF') and dataOut.flagSpreadF:
3970 if hasattr(dataOut, 'flagSpreadF') and dataOut.flagSpreadF:
3969 i2=int((620-dataOut.range1[0])/dataOut.DH)
3971 i2=int((620-dataOut.range1[0])/dataOut.DH)
3970 nanindex = numpy.argwhere(numpy.isnan(dataOut.ph2))
3972 nanindex = numpy.argwhere(numpy.isnan(dataOut.ph2))
3971 print("nanindex",nanindex)
3973 #print("nanindex",nanindex)
3972 i1 = nanindex[-1][0] #VER CUANDO i1>i2
3974 i1 = nanindex[-1][0] #VER CUANDO i1>i2
3973 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
3975 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
3974 #print("i1, i2",i1,i2)
3976 #print("i1, i2",i1,i2)
3975 print(dataOut.heightList)
3977 #print(dataOut.heightList)
3976 print("Bounds: ", dataOut.heightList[i1],dataOut.heightList[i2])
3978 #print("Bounds: ", dataOut.heightList[i1],dataOut.heightList[i2])
3977 #print(dataOut.dphi[33])
3979 #print(dataOut.dphi[33])
3978 #print(dataOut.ph2[33])
3980 #print(dataOut.ph2[33])
3979 #print(dataOut.dphi[i1::])
3981 #print(dataOut.dphi[i1::])
@@ -4008,9 +4010,9 class NormalizeDPPowerRoberto_V2(Operation):
4008 #'''
4010 #'''
4009 #if (time_text.hour == 5 and time_text.minute == 32): #Year: 2022, DOY:104
4011 #if (time_text.hour == 5 and time_text.minute == 32): #Year: 2022, DOY:104
4010 #if (time_text.hour == 0 and time_text.minute == 12): #Year: 2022, DOY:93
4012 #if (time_text.hour == 0 and time_text.minute == 12): #Year: 2022, DOY:93
4011 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
4013 #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
4012 #if (time_text.hour == 1 and time_text.minute == 23) or (time_text.hour == 1 and time_text.minute == 44): #Year: 2022, DOY:243
4014 #if (time_text.hour == 1 and time_text.minute == 23) or (time_text.hour == 1 and time_text.minute == 44): #Year: 2022, DOY:243
4013 dataOut.cf = dataOut.cflast[0]
4015 #dataOut.cf = dataOut.cflast[0]
4014 #if (time_text.hour == 0 and time_text.minute == 4): #Year: 2022, DOY:244
4016 #if (time_text.hour == 0 and time_text.minute == 4): #Year: 2022, DOY:244
4015 #dataOut.cf = 0.08
4017 #dataOut.cf = 0.08
4016 #print("here")
4018 #print("here")
@@ -4020,12 +4022,30 class NormalizeDPPowerRoberto_V2(Operation):
4020 #dataOut.cf = 0.09
4022 #dataOut.cf = 0.09
4021 #if (time_text.hour == 3 and time_text.minute == 59) or (time_text.hour == 4 and time_text.minute == 20): #Year: 2022, DOY:244
4023 #if (time_text.hour == 3 and time_text.minute == 59) or (time_text.hour == 4 and time_text.minute == 20): #Year: 2022, DOY:244
4022 #dataOut.cf = 0.09
4024 #dataOut.cf = 0.09
4025
4026 #if (time_text.hour == 7 and time_text.minute == 18) or (time_text.hour == 7 and time_text.minute == 40) or (time_text.hour == 7 and time_text.minute == 50): #Year: 2023, DOY:028
4027 #dataOut.cf = 0.029844088
4028 #if (time_text.hour == 0 and time_text.minute == 12) or (time_text.hour == 0 and time_text.minute == 22): #Year: 2023, DOY:028
4029 #if (time_text.hour == 7 and time_text.minute == 8) 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.minute == 1) or (time_text.hour == 9 and time_text.minute == 48) or (time_text.hour == 9 and time_text.minute == 58): #Year: 2023, DOY:029
4030 #dataOut.cf = dataOut.cflast[0]
4031 #if (time_text.hour == 7 and time_text.minute == 29): #Year: 2023, DOY:029
4032 #dataOut.cf = 0.025670702
4033 #if (time_text.hour == 8 and time_text.minute == 12): #Year: 2023, DOY:029
4034 #dataOut.cf = 0.036969288
4035 #if (time_text.hour == 8 and time_text.minute == 22): #Year: 2023, DOY:029
4036 #dataOut.cf = 0.0418645
4037 if (time_text.hour == 5 and time_text.minute == 10) or (time_text.hour == 5 and time_text.minute == 42) or (time_text.hour == 0 and time_text.minute == 22) or (time_text.hour == 0 and time_text.minute == 33): #Year: 2023, DOY:030
4038 dataOut.cf = dataOut.cflast[0]
4039 #if (time_text.hour == 0 and time_text.minute == 22) or (time_text.hour == 0 and time_text.minute == 33): #Year: 2023, DOY:031
4040 #if (time_text.hour == 7 and time_text.minute == 29): #Year: 2023, DOY:032
4041 #dataOut.cf = dataOut.cflast[0]
4042
4023 #'''
4043 #'''
4024 #dataOut.cf = 0.000057#0.0008136899
4044 #dataOut.cf = 0.000057#0.0008136899
4025
4045
4026
4046
4027 dataOut.cflast[0]=dataOut.cf
4047 dataOut.cflast[0]=dataOut.cf
4028 print("cf: ", dataOut.cf)
4048 #print("cf: ", dataOut.cf)
4029
4049
4030 #print(dataOut.ph2)
4050 #print(dataOut.ph2)
4031 #print(dataOut.sdp2)
4051 #print(dataOut.sdp2)
@@ -4049,7 +4069,7 class NormalizeDPPowerRoberto_V2(Operation):
4049 #print(dataOut.ph2)
4069 #print(dataOut.ph2)
4050 #print(dataOut.sdp2)
4070 #print(dataOut.sdp2)
4051 #input()
4071 #input()
4052 print("shape before" ,numpy.shape(dataOut.ph2))
4072 #print("shape before" ,numpy.shape(dataOut.ph2))
4053
4073
4054 return dataOut
4074 return dataOut
4055
4075
@@ -4227,15 +4247,18 class DPTemperaturesEstimation(Operation):
4227
4247
4228 eb=numpy.resize(eb,10)
4248 eb=numpy.resize(eb,10)
4229 dataOut.ifit=numpy.resize(dataOut.ifit,10)
4249 dataOut.ifit=numpy.resize(dataOut.ifit,10)
4250 #print("*********************FITACF_FIT*********************",dataOut.covinv,e,dataOut.params,eb,dataOut.m)
4251 #exit(1)
4230 with suppress_stdout_stderr():
4252 with suppress_stdout_stderr():
4231 dataOut.covinv,e,dataOut.params,eb,dataOut.m=fitacf_fit_short.fit(wl,x,y,dataOut.cov,dataOut.covinv,e,dataOut.params,bm,angle,den,dataOut.range1[i],dataOut.year,dataOut.ifit,dataOut.m,l1) #
4253 dataOut.covinv,e,dataOut.params,eb,dataOut.m=fitacf_fit_short.fit(wl,x,y,dataOut.cov,dataOut.covinv,e,dataOut.params,bm,angle,den,dataOut.range1[i],dataOut.year,dataOut.ifit,dataOut.m,l1) #
4232
4254 #exit(1)
4233 if dataOut.params[2]>dataOut.params[1]*1.05:
4255 if dataOut.params[2]>dataOut.params[1]*1.05:
4234 dataOut.ifit[2]=0
4256 dataOut.ifit[2]=0
4235 dataOut.params[1]=dataOut.params[2]=t1
4257 dataOut.params[1]=dataOut.params[2]=t1
4236 with suppress_stdout_stderr():
4258 with suppress_stdout_stderr():
4237 dataOut.covinv,e,dataOut.params,eb,dataOut.m=fitacf_fit_short.fit(wl,x,y,dataOut.cov,dataOut.covinv,e,dataOut.params,bm,angle,den,dataOut.range1[i],dataOut.year,dataOut.ifit,dataOut.m,l1) #
4259 dataOut.covinv,e,dataOut.params,eb,dataOut.m=fitacf_fit_short.fit(wl,x,y,dataOut.cov,dataOut.covinv,e,dataOut.params,bm,angle,den,dataOut.range1[i],dataOut.year,dataOut.ifit,dataOut.m,l1) #
4238
4260 #print("*********************FIT SUCCESS*********************",dataOut.covinv,e,dataOut.params,eb,dataOut.m)
4261 #exit(1)
4239 if (dataOut.ifit[2]==0):
4262 if (dataOut.ifit[2]==0):
4240 dataOut.params[2]=dataOut.params[1]
4263 dataOut.params[2]=dataOut.params[1]
4241 if (dataOut.ifit[3]==0 and iflag==0):
4264 if (dataOut.ifit[3]==0 and iflag==0):
@@ -4306,7 +4329,7 class DPTemperaturesEstimation(Operation):
4306 return dataOut
4329 return dataOut
4307
4330
4308
4331
4309 class DenCorrection(NormalizeDPPowerRoberto_V2):
4332 class DenCorrection_bckp_Dec_2022(NormalizeDPPowerRoberto_V2):
4310 '''
4333 '''
4311 Written by R. Flores
4334 Written by R. Flores
4312 '''
4335 '''
@@ -4322,6 +4345,7 class DenCorrection(NormalizeDPPowerRoberto_V2):
4322 def TeTiEstimation(self,dataOut):
4345 def TeTiEstimation(self,dataOut):
4323
4346
4324 y=numpy.zeros(dataOut.DPL,order='F',dtype='float32')
4347 y=numpy.zeros(dataOut.DPL,order='F',dtype='float32')
4348
4325 #y_aux = numpy.zeros(1,,dtype='float32')
4349 #y_aux = numpy.zeros(1,,dtype='float32')
4326 for i in range(dataOut.NSHTS):
4350 for i in range(dataOut.NSHTS):
4327 y[0]=y[1]=dataOut.range1[i]
4351 y[0]=y[1]=dataOut.range1[i]
@@ -4423,7 +4447,7 class DenCorrection(NormalizeDPPowerRoberto_V2):
4423 fion[2]=0.0 #0
4447 fion[2]=0.0 #0
4424 for j in range(dataOut.DPL):
4448 for j in range(dataOut.DPL):
4425 tau=dataOut.alag[j]*1.0e-3
4449 tau=dataOut.alag[j]*1.0e-3
4426
4450 #print("***********ACF2***********")
4427 with suppress_stdout_stderr():#The smoothness in range of "y" depends on the smoothness of the input parameters
4451 with suppress_stdout_stderr():#The smoothness in range of "y" depends on the smoothness of the input parameters
4428 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)
4452 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)
4429 #y[j]=fitacf_acf2.acf2(wl,tau,te2_smooth[i],tion,fion,nue,nui,wion,angle,dataOut.ph2[i],dataOut.bfm[i],y[j],three)
4453 #y[j]=fitacf_acf2.acf2(wl,tau,te2_smooth[i],tion,fion,nue,nui,wion,angle,dataOut.ph2[i],dataOut.bfm[i],y[j],three)
@@ -4506,6 +4530,128 class DenCorrection(NormalizeDPPowerRoberto_V2):
4506
4530
4507 return dataOut
4531 return dataOut
4508
4532
4533 class DenCorrection(NormalizeDPPowerRoberto_V2):
4534 '''
4535 Written by R. Flores
4536 '''
4537 def __init__(self, **kwargs):
4538
4539 Operation.__init__(self, **kwargs)
4540 self.aux = 0
4541
4542 def gaussian(self, x, a, b, c):
4543 val = a * numpy.exp(-(x - b)**2 / (2*c**2))
4544 return val
4545
4546 def TeTiEstimation(self,dataOut):
4547
4548 dataOut.DPL = 2 #for MST
4549 y=numpy.zeros(dataOut.DPL,order='F',dtype='float32')
4550
4551 #y_aux = numpy.zeros(1,,dtype='float32')
4552 for i in range(dataOut.NSHTS):
4553 y[0]=y[1]=dataOut.range1[i]
4554
4555 y = y.astype(dtype='float64',order='F')
4556 three=int(3)
4557 wl = 3.0
4558 tion=numpy.zeros(three,order='F',dtype='float32')
4559 fion=numpy.zeros(three,order='F',dtype='float32')
4560 nui=numpy.zeros(three,order='F',dtype='float32')
4561 wion=numpy.zeros(three,order='F',dtype='int32')
4562 bline=0.0
4563 #bline=numpy.zeros(1,order='F',dtype='float32')
4564 my_aux = numpy.ones(dataOut.NSHTS,order='F',dtype='float32')
4565 acf_Temps = numpy.ones(dataOut.NSHTS,order='F',dtype='float32')*numpy.nan
4566 acf_no_Temps = numpy.ones(dataOut.NSHTS,order='F',dtype='float32')*numpy.nan
4567
4568 from scipy import signal
4569
4570 def func(params):
4571 return (ratio2-self.gaussian(dataOut.heightList[:dataOut.NSHTS],params[0],params[1],params[2]))
4572
4573 #print("Before loop")
4574 dataOut.info2[0] = 1
4575 for i in range(dataOut.NSHTS):
4576 if dataOut.info2[i]==1:
4577 angle=dataOut.thb[i]*0.01745
4578 nue=nui[0]=nui[1]=nui[2]=0.0#nui[3]=0.0
4579 wion[0]=16 #O
4580 wion[1]=1 #H
4581 wion[2]=4 #He
4582 tion[0]=tion[1]=tion[2]=dataOut.ti2[i]
4583 #tion[0]=tion[1]=tion[2]=ti2_smooth[i]
4584 fion[0]=1.0-dataOut.phy2[i] #1
4585 fion[1]=dataOut.phy2[i] #0
4586 fion[2]=0.0 #0
4587 for j in range(dataOut.DPL):
4588 tau=dataOut.alag[j]*1.0e-3
4589 #print("***********ACF2***********")
4590 with suppress_stdout_stderr():#The smoothness in range of "y" depends on the smoothness of the input parameters
4591 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)
4592 #y[j]=fitacf_acf2.acf2(wl,tau,te2_smooth[i],tion,fion,nue,nui,wion,angle,dataOut.ph2[i],dataOut.bfm[i],y[j],three)
4593 #print("i: ", i, "j: ", j)
4594 #if i == 0 and j == 1:
4595 #print("Params: ",wl,tau,dataOut.te2[i],tion,fion,nue,nui,wion,angle,dataOut.ph2[i],dataOut.bfm[i],y[j],three)
4596 #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)
4597 #exit(1)
4598 #if dataOut.ut_Faraday>11.0 and dataOut.range1[i]>150.0 and dataOut.range1[i]<400.0:
4599
4600 if dataOut.ut_Faraday>11.0 and dataOut.range1[i]>150.0 and dataOut.range1[i]<300.0:
4601 #if dataOut.ut_Faraday>11.0 and dataOut.range1[i]>150.0 and dataOut.range1[i]<400.0:
4602 tau=0.0
4603 with suppress_stdout_stderr():
4604 bline=fitacf_acf2.acf2(wl,tau,tion,tion,fion,nue,nui,wion,angle,dataOut.ph2[i],dataOut.bfm[i],bline,three)
4605
4606 #if i == 0 and j == 1:
4607 #print("bline: ",bline)
4608 #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)
4609 #exit(1)
4610 cf=min(1.2,max(1.0,bline/y[0])) #FACTOR DE EFICIENCIA
4611 #cf = bline/y[0]
4612 #cf=min(2.,max(1.0,bline/y[0]))
4613 my_aux[i] = cf
4614 acf_Temps[i] = y[0]
4615 acf_no_Temps[i] = bline
4616 #dataOut.ph2[i]=cf*dataOut.ph2[i] #Instead we adjust the curve "cf" into a Gaussian,
4617 #dataOut.sdp2[i]=cf*dataOut.sdp2[i] #in order to get smoother values of density
4618 for j in range(1,dataOut.DPL):
4619 #y[j]=(y[j]/y[0])*dataOut.DH+dataOut.range1[i]
4620 y[j]=min(max((y[j]/y[0]),-1.0),1.0)*dataOut.DH+dataOut.range1[i]
4621 y[0]=dataOut.range1[i]+dataOut.DH
4622
4623
4624 ratio = my_aux-1
4625 #ratio = dataOut.te2[:dataOut.NSHTS]/dataOut.ti2[:dataOut.NSHTS]
4626 def lsq_func(params):
4627 return (ratio-self.gaussian(dataOut.heightList[:dataOut.NSHTS],params[0],params[1],params[2]))
4628
4629 x0_value = numpy.array([max(ratio),250,20])
4630
4631 popt = least_squares(lsq_func,x0=x0_value,verbose=0)
4632
4633 A = popt.x[0]; B = popt.x[1]; C = popt.x[2]
4634
4635 aux = self.gaussian(dataOut.heightList[:dataOut.NSHTS], A, B, C) + 1 #ratio + 1
4636
4637 dataOut.ph2[:dataOut.NSHTS]*=aux
4638 dataOut.sdp2[:dataOut.NSHTS]*=aux
4639 #dataOut.ph2[:26]*=aux[:26]
4640 #dataOut.sdp2[:26]*=aux[:26]
4641 #print(aux)
4642 #print("inside correction",dataOut.ph2)
4643
4644 def run(self,dataOut):
4645 #print("hour",gmtime(dataOut.utctime).tm_hour)
4646 if gmtime(dataOut.utctime).tm_hour < 24. and gmtime(dataOut.utctime).tm_hour >= 11.:
4647 #print("inside")
4648 self.TeTiEstimation(dataOut)
4649 dataOut.flagTeTiCorrection = True
4650
4651 self.normalize(dataOut)
4652
4653 return dataOut
4654
4509 class DataPlotCleaner(Operation):
4655 class DataPlotCleaner(Operation):
4510 '''
4656 '''
4511 Written by R. Flores
4657 Written by R. Flores
@@ -4825,10 +4971,16 class DataSaveCleaner(Operation):
4825
4971
4826 #print(time_text.hour,time_text.minute)
4972 #print(time_text.hour,time_text.minute)
4827 #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
4973 #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
4828 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
4974 #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
4829 #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
4975 #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
4830 #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
4976 #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
4831 #if (time_text.hour >= 8 and time_text.hour < 11) or (time_text.hour == 1) or (time_text.hour == 0 and time_text.minute==15) 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) or (time_text.hour == 3 and time_text.minute==16) or (time_text.hour == 3 and time_text.minute==27): #Year: 2022, DOY:244
4977 #if (time_text.hour >= 8 and time_text.hour < 11) or (time_text.hour == 1) or (time_text.hour == 0 and time_text.minute==15) 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) or (time_text.hour == 3 and time_text.minute==16) or (time_text.hour == 3 and time_text.minute==27): #Year: 2022, DOY:244
4978 #if (time_text.hour == 0 and time_text.minute==8) or (time_text.hour == 0 and time_text.minute==18): #Year: 2023, DOY:027
4979 #if (time_text.hour >= 3 and time_text.hour<5) or (time_text.hour == 6) or (time_text.hour == 7 and time_text.minute==8) or (time_text.hour == 8 and time_text.minute==1) 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 == 9 and time_text.minute>=16) or (time_text.hour == 10) or (time_text.hour == 11 and time_text.minute==2): #Year: 2023, DOY:028
4980 #if (time_text.hour >=2 and time_text.hour<5) or (time_text.hour == 10): #Year: 2023, DOY:029
4981 if (time_text.hour == 6 and time_text.minute==36) or (time_text.hour == 6 and time_text.minute==46) or (time_text.hour == 6 and time_text.minute==57) or (time_text.hour == 7 and time_text.minute==8) or (time_text.hour == 8) or (time_text.hour == 9) or (time_text.hour == 10) or (time_text.hour == 14 and time_text.minute==25) or (time_text.hour == 14 and time_text.minute==57): #Year: 2023, DOY:030
4982 #if (time_text.hour == 0 and time_text.minute==54) or (time_text.hour >= 1 and time_text.hour<5) or (time_text.hour == 5) or (time_text.hour == 6) or (time_text.hour == 7) or (time_text.hour == 8) or (time_text.hour == 9 and time_text.minute != 58): #Year: 2023, DOY:031
4983 #if (time_text.hour == 5 and time_text.minute<=21) or (time_text.hour == 8 and time_text.minute==1) or (time_text.hour == 9 and time_text.minute==2) or (time_text.hour == 10) or (time_text.hour == 12 and time_text.minute==31): #Year: 2023, DOY:032
4832
4984
4833 dataOut.DensityFinal[0,:]=missing
4985 dataOut.DensityFinal[0,:]=missing
4834 dataOut.EDensityFinal[0,:]=missing
4986 dataOut.EDensityFinal[0,:]=missing
@@ -4842,9 +4994,19 class DataSaveCleaner(Operation):
4842 dataOut.flagNoData = True #Remueve todo el perfil
4994 dataOut.flagNoData = True #Remueve todo el perfil
4843 #'''
4995 #'''
4844 '''
4996 '''
4845 #if (time_text.hour >= 7 and time_text.hour < 9): #Year: 2022, DOY:102
4997 if (time_text.hour == 5) or (time_text.hour == 6) or (time_text.hour == 7) or (time_text.hour == 9): #Year: 2023, DOY:032
4846 if (time_text.hour == 20 and time_text.minute == 8): #Year: 2022, DOY:243
4998 id_aux = 11
4847 id_aux = 35
4999 dataOut.DensityFinal[0,:id_aux]=missing
5000 dataOut.EDensityFinal[0,:id_aux]=missing
5001 dataOut.ElecTempFinal[0,:id_aux]=missing
5002 dataOut.EElecTempFinal[0,:id_aux]=missing
5003 dataOut.IonTempFinal[0,:id_aux]=missing
5004 dataOut.EIonTempFinal[0,:id_aux]=missing
5005 dataOut.PhyFinal[0,:id_aux]=missing
5006 dataOut.EPhyFinal[0,:id_aux]=missing
5007
5008 if (time_text.hour == 6 and time_text.minute == 57) or (time_text.hour == 7 and time_text.minute == 29): #Year: 2023, DOY:032
5009 id_aux = 36
4848 dataOut.DensityFinal[0,id_aux:]=missing
5010 dataOut.DensityFinal[0,id_aux:]=missing
4849 dataOut.EDensityFinal[0,id_aux:]=missing
5011 dataOut.EDensityFinal[0,id_aux:]=missing
4850 dataOut.ElecTempFinal[0,id_aux:]=missing
5012 dataOut.ElecTempFinal[0,id_aux:]=missing
@@ -4853,8 +5015,9 class DataSaveCleaner(Operation):
4853 dataOut.EIonTempFinal[0,id_aux:]=missing
5015 dataOut.EIonTempFinal[0,id_aux:]=missing
4854 dataOut.PhyFinal[0,id_aux:]=missing
5016 dataOut.PhyFinal[0,id_aux:]=missing
4855 dataOut.EPhyFinal[0,id_aux:]=missing
5017 dataOut.EPhyFinal[0,id_aux:]=missing
4856 if (time_text.hour == 20 and time_text.minute == 19): #Year: 2022, DOY:243
5018
4857 id_aux = 33
5019 if (time_text.hour == 9 and time_text.minute == 24): #Year: 2023, DOY:032
5020 id_aux = 30
4858 dataOut.DensityFinal[0,id_aux:]=missing
5021 dataOut.DensityFinal[0,id_aux:]=missing
4859 dataOut.EDensityFinal[0,id_aux:]=missing
5022 dataOut.EDensityFinal[0,id_aux:]=missing
4860 dataOut.ElecTempFinal[0,id_aux:]=missing
5023 dataOut.ElecTempFinal[0,id_aux:]=missing
@@ -4863,8 +5026,9 class DataSaveCleaner(Operation):
4863 dataOut.EIonTempFinal[0,id_aux:]=missing
5026 dataOut.EIonTempFinal[0,id_aux:]=missing
4864 dataOut.PhyFinal[0,id_aux:]=missing
5027 dataOut.PhyFinal[0,id_aux:]=missing
4865 dataOut.EPhyFinal[0,id_aux:]=missing
5028 dataOut.EPhyFinal[0,id_aux:]=missing
4866 if (time_text.hour == 20 and time_text.minute == 29): #Year: 2022, DOY:243
5029
4867 id_aux = 30
5030 if (time_text.hour == 5 and time_text.minute == 32): #Year: 2023, DOY:032
5031 id_aux = 41
4868 dataOut.DensityFinal[0,id_aux:]=missing
5032 dataOut.DensityFinal[0,id_aux:]=missing
4869 dataOut.EDensityFinal[0,id_aux:]=missing
5033 dataOut.EDensityFinal[0,id_aux:]=missing
4870 dataOut.ElecTempFinal[0,id_aux:]=missing
5034 dataOut.ElecTempFinal[0,id_aux:]=missing
@@ -4873,8 +5037,9 class DataSaveCleaner(Operation):
4873 dataOut.EIonTempFinal[0,id_aux:]=missing
5037 dataOut.EIonTempFinal[0,id_aux:]=missing
4874 dataOut.PhyFinal[0,id_aux:]=missing
5038 dataOut.PhyFinal[0,id_aux:]=missing
4875 dataOut.EPhyFinal[0,id_aux:]=missing
5039 dataOut.EPhyFinal[0,id_aux:]=missing
4876 if (time_text.hour == 20 and time_text.minute == 44): #Year: 2022, DOY:243
5040
4877 id_aux = 31
5041 if (time_text.hour == 11 and time_text.minute == 14): #Year: 2023, DOY:032
5042 id_aux = 35
4878 dataOut.DensityFinal[0,id_aux:]=missing
5043 dataOut.DensityFinal[0,id_aux:]=missing
4879 dataOut.EDensityFinal[0,id_aux:]=missing
5044 dataOut.EDensityFinal[0,id_aux:]=missing
4880 dataOut.ElecTempFinal[0,id_aux:]=missing
5045 dataOut.ElecTempFinal[0,id_aux:]=missing
@@ -4883,8 +5048,10 class DataSaveCleaner(Operation):
4883 dataOut.EIonTempFinal[0,id_aux:]=missing
5048 dataOut.EIonTempFinal[0,id_aux:]=missing
4884 dataOut.PhyFinal[0,id_aux:]=missing
5049 dataOut.PhyFinal[0,id_aux:]=missing
4885 dataOut.EPhyFinal[0,id_aux:]=missing
5050 dataOut.EPhyFinal[0,id_aux:]=missing
4886 if (time_text.hour <= 8): #Year: 2022, DOY:243
5051 '''
4887 id_aux = 11
5052 '''
5053 if (time_text.hour == 23 and time_text.minute == 50) or (time_text.hour == 0): #Year: 2023, DOY:031
5054 id_aux = 18
4888 dataOut.DensityFinal[0,:id_aux]=missing
5055 dataOut.DensityFinal[0,:id_aux]=missing
4889 dataOut.EDensityFinal[0,:id_aux]=missing
5056 dataOut.EDensityFinal[0,:id_aux]=missing
4890 dataOut.ElecTempFinal[0,:id_aux]=missing
5057 dataOut.ElecTempFinal[0,:id_aux]=missing
@@ -4893,7 +5060,8 class DataSaveCleaner(Operation):
4893 dataOut.EIonTempFinal[0,:id_aux]=missing
5060 dataOut.EIonTempFinal[0,:id_aux]=missing
4894 dataOut.PhyFinal[0,:id_aux]=missing
5061 dataOut.PhyFinal[0,:id_aux]=missing
4895 dataOut.EPhyFinal[0,:id_aux]=missing
5062 dataOut.EPhyFinal[0,:id_aux]=missing
4896 if (time_text.hour == 23): #Year: 2022, DOY:243
5063
5064 if (time_text.hour == 23 and time_text.minute == 40) or (time_text.hour == 9 and time_text.minute == 58) or (time_text.hour == 10 and time_text.minute == 20): #Year: 2023, DOY:031
4897 id_aux = 12
5065 id_aux = 12
4898 dataOut.DensityFinal[0,:id_aux]=missing
5066 dataOut.DensityFinal[0,:id_aux]=missing
4899 dataOut.EDensityFinal[0,:id_aux]=missing
5067 dataOut.EDensityFinal[0,:id_aux]=missing
@@ -4903,28 +5071,66 class DataSaveCleaner(Operation):
4903 dataOut.EIonTempFinal[0,:id_aux]=missing
5071 dataOut.EIonTempFinal[0,:id_aux]=missing
4904 dataOut.PhyFinal[0,:id_aux]=missing
5072 dataOut.PhyFinal[0,:id_aux]=missing
4905 dataOut.EPhyFinal[0,:id_aux]=missing
5073 dataOut.EPhyFinal[0,:id_aux]=missing
4906 if (time_text.hour == 5 and time_text.minute == 21): #Year: 2022, DOY:243
5074
4907 id_aux = (36,37)
5075 if (time_text.hour == 11 and time_text.minute == 45): #Year: 2023, DOY:031
4908 dataOut.DensityFinal[0,id_aux]=missing
5076 id_aux = 8
4909 dataOut.EDensityFinal[0,id_aux]=missing
5077 dataOut.DensityFinal[0,:id_aux]=missing
4910 dataOut.ElecTempFinal[0,id_aux]=missing
5078 dataOut.EDensityFinal[0,:id_aux]=missing
4911 dataOut.EElecTempFinal[0,id_aux]=missing
5079 dataOut.ElecTempFinal[0,:id_aux]=missing
4912 dataOut.IonTempFinal[0,id_aux]=missing
5080 dataOut.EElecTempFinal[0,:id_aux]=missing
4913 dataOut.EIonTempFinal[0,id_aux]=missing
5081 dataOut.IonTempFinal[0,:id_aux]=missing
4914 dataOut.PhyFinal[0,id_aux]=missing
5082 dataOut.EIonTempFinal[0,:id_aux]=missing
4915 dataOut.EPhyFinal[0,id_aux]=missing
5083 dataOut.PhyFinal[0,:id_aux]=missing
4916 if (time_text.hour == 5 and time_text.minute == 53): #Year: 2022, DOY:243
5084 dataOut.EPhyFinal[0,:id_aux]=missing
4917 id_aux = (37,38)
5085
4918 dataOut.DensityFinal[0,id_aux]=missing
5086 if (time_text.hour == 13): #Year: 2023, DOY:031
4919 dataOut.EDensityFinal[0,id_aux]=missing
5087 id_aux = 37
4920 dataOut.ElecTempFinal[0,id_aux]=missing
5088 dataOut.DensityFinal[0,id_aux:]=missing
4921 dataOut.EElecTempFinal[0,id_aux]=missing
5089 dataOut.EDensityFinal[0,id_aux:]=missing
4922 dataOut.IonTempFinal[0,id_aux]=missing
5090 dataOut.ElecTempFinal[0,id_aux:]=missing
4923 dataOut.EIonTempFinal[0,id_aux]=missing
5091 dataOut.EElecTempFinal[0,id_aux:]=missing
4924 dataOut.PhyFinal[0,id_aux]=missing
5092 dataOut.IonTempFinal[0,id_aux:]=missing
4925 dataOut.EPhyFinal[0,id_aux]=missing
5093 dataOut.EIonTempFinal[0,id_aux:]=missing
4926 if (time_text.hour == 6 and time_text.minute == 4): #Year: 2022, DOY:243
5094 dataOut.PhyFinal[0,id_aux:]=missing
4927 id_aux = (38,39)
5095 dataOut.EPhyFinal[0,id_aux:]=missing
5096
5097 if (time_text.hour == 15 and time_text.minute == 18): #Year: 2023, DOY:031
5098 id_aux = 34
5099 dataOut.DensityFinal[0,id_aux:]=missing
5100 dataOut.EDensityFinal[0,id_aux:]=missing
5101 dataOut.ElecTempFinal[0,id_aux:]=missing
5102 dataOut.EElecTempFinal[0,id_aux:]=missing
5103 dataOut.IonTempFinal[0,id_aux:]=missing
5104 dataOut.EIonTempFinal[0,id_aux:]=missing
5105 dataOut.PhyFinal[0,id_aux:]=missing
5106 dataOut.EPhyFinal[0,id_aux:]=missing
5107
5108 if (time_text.hour == 10 and time_text.minute == 52): #Year: 2023, DOY:031
5109 id_aux = 29
5110 dataOut.DensityFinal[0,id_aux:]=missing
5111 dataOut.EDensityFinal[0,id_aux:]=missing
5112 dataOut.ElecTempFinal[0,id_aux:]=missing
5113 dataOut.EElecTempFinal[0,id_aux:]=missing
5114 dataOut.IonTempFinal[0,id_aux:]=missing
5115 dataOut.EIonTempFinal[0,id_aux:]=missing
5116 dataOut.PhyFinal[0,id_aux:]=missing
5117 dataOut.EPhyFinal[0,id_aux:]=missing
5118
5119 if (time_text.hour == 10 and time_text.minute == 9): #Year: 2023, DOY:031
5120 id_aux = 33
5121 dataOut.DensityFinal[0,id_aux:]=missing
5122 dataOut.EDensityFinal[0,id_aux:]=missing
5123 dataOut.ElecTempFinal[0,id_aux:]=missing
5124 dataOut.EElecTempFinal[0,id_aux:]=missing
5125 dataOut.IonTempFinal[0,id_aux:]=missing
5126 dataOut.EIonTempFinal[0,id_aux:]=missing
5127 dataOut.PhyFinal[0,id_aux:]=missing
5128 dataOut.EPhyFinal[0,id_aux:]=missing
5129 '''
5130
5131 #'''
5132 if (time_text.hour == 0 and time_text.minute == 33): #Year: 2023, DOY:030
5133 id_aux = 38
4928 dataOut.DensityFinal[0,id_aux]=missing
5134 dataOut.DensityFinal[0,id_aux]=missing
4929 dataOut.EDensityFinal[0,id_aux]=missing
5135 dataOut.EDensityFinal[0,id_aux]=missing
4930 dataOut.ElecTempFinal[0,id_aux]=missing
5136 dataOut.ElecTempFinal[0,id_aux]=missing
@@ -4933,8 +5139,9 class DataSaveCleaner(Operation):
4933 dataOut.EIonTempFinal[0,id_aux]=missing
5139 dataOut.EIonTempFinal[0,id_aux]=missing
4934 dataOut.PhyFinal[0,id_aux]=missing
5140 dataOut.PhyFinal[0,id_aux]=missing
4935 dataOut.EPhyFinal[0,id_aux]=missing
5141 dataOut.EPhyFinal[0,id_aux]=missing
4936 if (time_text.hour == 12 and time_text.minute == 6): #Year: 2022, DOY:243
5142
4937 id_aux = (29,30)
5143 if (time_text.hour == 18 and time_text.minute == 41): #Year: 2023, DOY:030
5144 id_aux = 49
4938 dataOut.DensityFinal[0,id_aux]=missing
5145 dataOut.DensityFinal[0,id_aux]=missing
4939 dataOut.EDensityFinal[0,id_aux]=missing
5146 dataOut.EDensityFinal[0,id_aux]=missing
4940 dataOut.ElecTempFinal[0,id_aux]=missing
5147 dataOut.ElecTempFinal[0,id_aux]=missing
@@ -4943,8 +5150,9 class DataSaveCleaner(Operation):
4943 dataOut.EIonTempFinal[0,id_aux]=missing
5150 dataOut.EIonTempFinal[0,id_aux]=missing
4944 dataOut.PhyFinal[0,id_aux]=missing
5151 dataOut.PhyFinal[0,id_aux]=missing
4945 dataOut.EPhyFinal[0,id_aux]=missing
5152 dataOut.EPhyFinal[0,id_aux]=missing
4946 if (time_text.hour == 14 and time_text.minute == 14): #Year: 2022, DOY:243
5153
4947 id_aux = (35,36)
5154 if (time_text.hour == 14 and time_text.minute == 14): #Year: 2023, DOY:030
5155 id_aux = 34
4948 dataOut.DensityFinal[0,id_aux]=missing
5156 dataOut.DensityFinal[0,id_aux]=missing
4949 dataOut.EDensityFinal[0,id_aux]=missing
5157 dataOut.EDensityFinal[0,id_aux]=missing
4950 dataOut.ElecTempFinal[0,id_aux]=missing
5158 dataOut.ElecTempFinal[0,id_aux]=missing
@@ -4953,8 +5161,9 class DataSaveCleaner(Operation):
4953 dataOut.EIonTempFinal[0,id_aux]=missing
5161 dataOut.EIonTempFinal[0,id_aux]=missing
4954 dataOut.PhyFinal[0,id_aux]=missing
5162 dataOut.PhyFinal[0,id_aux]=missing
4955 dataOut.EPhyFinal[0,id_aux]=missing
5163 dataOut.EPhyFinal[0,id_aux]=missing
4956 if (time_text.hour == 23 and time_text.minute == 2): #Year: 2022, DOY:243
5164
4957 id_aux = (41,42)
5165 if (time_text.hour == 16 and time_text.minute == 1): #Year: 2023, DOY:030
5166 id_aux = (40,41)
4958 dataOut.DensityFinal[0,id_aux]=missing
5167 dataOut.DensityFinal[0,id_aux]=missing
4959 dataOut.EDensityFinal[0,id_aux]=missing
5168 dataOut.EDensityFinal[0,id_aux]=missing
4960 dataOut.ElecTempFinal[0,id_aux]=missing
5169 dataOut.ElecTempFinal[0,id_aux]=missing
@@ -4963,8 +5172,9 class DataSaveCleaner(Operation):
4963 dataOut.EIonTempFinal[0,id_aux]=missing
5172 dataOut.EIonTempFinal[0,id_aux]=missing
4964 dataOut.PhyFinal[0,id_aux]=missing
5173 dataOut.PhyFinal[0,id_aux]=missing
4965 dataOut.EPhyFinal[0,id_aux]=missing
5174 dataOut.EPhyFinal[0,id_aux]=missing
4966 if (time_text.hour == 0 and time_text.minute == 8): #Year: 2022, DOY:243
5175
4967 id_aux = 33
5176 if (time_text.hour == 5 and time_text.minute == 53) or (time_text.hour == 15 and time_text.minute == 40): #Year: 2023, DOY:030
5177 id_aux = 40
4968 dataOut.DensityFinal[0,id_aux:]=missing
5178 dataOut.DensityFinal[0,id_aux:]=missing
4969 dataOut.EDensityFinal[0,id_aux:]=missing
5179 dataOut.EDensityFinal[0,id_aux:]=missing
4970 dataOut.ElecTempFinal[0,id_aux:]=missing
5180 dataOut.ElecTempFinal[0,id_aux:]=missing
@@ -4973,8 +5183,434 class DataSaveCleaner(Operation):
4973 dataOut.EIonTempFinal[0,id_aux:]=missing
5183 dataOut.EIonTempFinal[0,id_aux:]=missing
4974 dataOut.PhyFinal[0,id_aux:]=missing
5184 dataOut.PhyFinal[0,id_aux:]=missing
4975 dataOut.EPhyFinal[0,id_aux:]=missing
5185 dataOut.EPhyFinal[0,id_aux:]=missing
4976 id_aux = 18
5186
4977 dataOut.DensityFinal[0,id_aux]=missing
5187 if (time_text.hour == 6 and time_text.minute == 25) or (time_text.hour == 15 and time_text.minute == 8) or (time_text.hour == 13 and time_text.minute == 42) or (time_text.hour == 13 and time_text.minute == 0): #Year: 2023, DOY:030
5188 id_aux = 37
5189 dataOut.DensityFinal[0,id_aux:]=missing
5190 dataOut.EDensityFinal[0,id_aux:]=missing
5191 dataOut.ElecTempFinal[0,id_aux:]=missing
5192 dataOut.EElecTempFinal[0,id_aux:]=missing
5193 dataOut.IonTempFinal[0,id_aux:]=missing
5194 dataOut.EIonTempFinal[0,id_aux:]=missing
5195 dataOut.PhyFinal[0,id_aux:]=missing
5196 dataOut.EPhyFinal[0,id_aux:]=missing
5197
5198 if (time_text.hour == 7 and time_text.minute == 40): #Year: 2023, DOY:030
5199 id_aux = 31
5200 dataOut.DensityFinal[0,id_aux:]=missing
5201 dataOut.EDensityFinal[0,id_aux:]=missing
5202 dataOut.ElecTempFinal[0,id_aux:]=missing
5203 dataOut.EElecTempFinal[0,id_aux:]=missing
5204 dataOut.IonTempFinal[0,id_aux:]=missing
5205 dataOut.EIonTempFinal[0,id_aux:]=missing
5206 dataOut.PhyFinal[0,id_aux:]=missing
5207 dataOut.EPhyFinal[0,id_aux:]=missing
5208
5209 if (time_text.hour == 11 and time_text.minute == 2): #Year: 2023, DOY:030
5210 id_aux = 24
5211 dataOut.DensityFinal[0,id_aux:]=missing
5212 dataOut.EDensityFinal[0,id_aux:]=missing
5213 dataOut.ElecTempFinal[0,id_aux:]=missing
5214 dataOut.EElecTempFinal[0,id_aux:]=missing
5215 dataOut.IonTempFinal[0,id_aux:]=missing
5216 dataOut.EIonTempFinal[0,id_aux:]=missing
5217 dataOut.PhyFinal[0,id_aux:]=missing
5218 dataOut.EPhyFinal[0,id_aux:]=missing
5219
5220 if (time_text.hour == 5) or (time_text.hour == 6) or (time_text.hour == 7) or (time_text.hour == 0): #Year: 2023, DOY:030
5221 id_aux = 12
5222 dataOut.DensityFinal[0,:id_aux]=missing
5223 dataOut.EDensityFinal[0,:id_aux]=missing
5224 dataOut.ElecTempFinal[0,:id_aux]=missing
5225 dataOut.EElecTempFinal[0,:id_aux]=missing
5226 dataOut.IonTempFinal[0,:id_aux]=missing
5227 dataOut.EIonTempFinal[0,:id_aux]=missing
5228 dataOut.PhyFinal[0,:id_aux]=missing
5229 dataOut.EPhyFinal[0,:id_aux]=missing
5230
5231 if (time_text.hour == 23 and time_text.minute == 40) or (time_text.hour == 23 and time_text.minute == 50): #Year: 2023, DOY:030
5232 id_aux = 13
5233 dataOut.DensityFinal[0,:id_aux]=missing
5234 dataOut.EDensityFinal[0,:id_aux]=missing
5235 dataOut.ElecTempFinal[0,:id_aux]=missing
5236 dataOut.EElecTempFinal[0,:id_aux]=missing
5237 dataOut.IonTempFinal[0,:id_aux]=missing
5238 dataOut.EIonTempFinal[0,:id_aux]=missing
5239 dataOut.PhyFinal[0,:id_aux]=missing
5240 dataOut.EPhyFinal[0,:id_aux]=missing
5241
5242 if (time_text.hour == 11 and time_text.minute == 56): #Year: 2023, DOY:030
5243 id_aux = 8
5244 dataOut.DensityFinal[0,:id_aux]=missing
5245 dataOut.EDensityFinal[0,:id_aux]=missing
5246 dataOut.ElecTempFinal[0,:id_aux]=missing
5247 dataOut.EElecTempFinal[0,:id_aux]=missing
5248 dataOut.IonTempFinal[0,:id_aux]=missing
5249 dataOut.EIonTempFinal[0,:id_aux]=missing
5250 dataOut.PhyFinal[0,:id_aux]=missing
5251 dataOut.EPhyFinal[0,:id_aux]=missing
5252 #'''
5253
5254 '''
5255 if (time_text.hour == 14 and time_text.minute == 42): #Year: 2023, DOY:027
5256 id_aux = 31
5257 dataOut.DensityFinal[0,id_aux:]=missing
5258 dataOut.EDensityFinal[0,id_aux:]=missing
5259 dataOut.ElecTempFinal[0,id_aux:]=missing
5260 dataOut.EElecTempFinal[0,id_aux:]=missing
5261 dataOut.IonTempFinal[0,id_aux:]=missing
5262 dataOut.EIonTempFinal[0,id_aux:]=missing
5263 dataOut.PhyFinal[0,id_aux:]=missing
5264 dataOut.EPhyFinal[0,id_aux:]=missing
5265
5266 if (time_text.hour == 15 and time_text.minute == 25) or (time_text.hour == 15 and time_text.minute == 36): #Year: 2023, DOY:027
5267 id_aux = 35
5268 dataOut.DensityFinal[0,id_aux:]=missing
5269 dataOut.EDensityFinal[0,id_aux:]=missing
5270 dataOut.ElecTempFinal[0,id_aux:]=missing
5271 dataOut.EElecTempFinal[0,id_aux:]=missing
5272 dataOut.IonTempFinal[0,id_aux:]=missing
5273 dataOut.EIonTempFinal[0,id_aux:]=missing
5274 dataOut.PhyFinal[0,id_aux:]=missing
5275 dataOut.EPhyFinal[0,id_aux:]=missing
5276
5277 if (time_text.hour == 23 and time_text.minute == 36) or (time_text.hour == 23 and time_text.minute == 46) or (time_text.hour == 23 and time_text.minute == 57): #Year: 2023, DOY:027
5278 id_aux = 13
5279 dataOut.DensityFinal[0,:id_aux]=missing
5280 dataOut.EDensityFinal[0,:id_aux]=missing
5281 dataOut.ElecTempFinal[0,:id_aux]=missing
5282 dataOut.EElecTempFinal[0,:id_aux]=missing
5283 dataOut.IonTempFinal[0,:id_aux]=missing
5284 dataOut.EIonTempFinal[0,:id_aux]=missing
5285 dataOut.PhyFinal[0,:id_aux]=missing
5286 dataOut.EPhyFinal[0,:id_aux]=missing
5287 '''
5288
5289 '''
5290 if (time_text.hour == 23 and time_text.minute == 40) or (time_text.hour == 23 and time_text.minute == 50) or (time_text.hour == 0 and time_text.minute <= 22) or (time_text.hour == 5 and time_text.minute == 10) or (time_text.hour == 6) or (time_text.hour == 7): #Year: 2023, DOY:029
5291 id_aux = 14
5292 dataOut.DensityFinal[0,:id_aux]=missing
5293 dataOut.EDensityFinal[0,:id_aux]=missing
5294 dataOut.ElecTempFinal[0,:id_aux]=missing
5295 dataOut.EElecTempFinal[0,:id_aux]=missing
5296 dataOut.IonTempFinal[0,:id_aux]=missing
5297 dataOut.EIonTempFinal[0,:id_aux]=missing
5298 dataOut.PhyFinal[0,:id_aux]=missing
5299 dataOut.EPhyFinal[0,:id_aux]=missing
5300
5301 if (time_text.hour == 8) or (time_text.hour == 9 and time_text.minute != 58) or (time_text.hour == 11 and time_text.minute == 2) or (time_text.hour == 5 and time_text.minute == 21) or (time_text.hour == 5 and time_text.minute == 42) or (time_text.hour == 5 and time_text.minute == 53): #Year: 2023, DOY:029
5302 id_aux = 12
5303 dataOut.DensityFinal[0,:id_aux]=missing
5304 dataOut.EDensityFinal[0,:id_aux]=missing
5305 dataOut.ElecTempFinal[0,:id_aux]=missing
5306 dataOut.EElecTempFinal[0,:id_aux]=missing
5307 dataOut.IonTempFinal[0,:id_aux]=missing
5308 dataOut.EIonTempFinal[0,:id_aux]=missing
5309 dataOut.PhyFinal[0,:id_aux]=missing
5310 dataOut.EPhyFinal[0,:id_aux]=missing
5311
5312 if (time_text.hour == 5 and time_text.minute == 32): #Year: 2023, DOY:029
5313 id_aux = 16
5314 dataOut.DensityFinal[0,:id_aux]=missing
5315 dataOut.EDensityFinal[0,:id_aux]=missing
5316 dataOut.ElecTempFinal[0,:id_aux]=missing
5317 dataOut.EElecTempFinal[0,:id_aux]=missing
5318 dataOut.IonTempFinal[0,:id_aux]=missing
5319 dataOut.EIonTempFinal[0,:id_aux]=missing
5320 dataOut.PhyFinal[0,:id_aux]=missing
5321 dataOut.EPhyFinal[0,:id_aux]=missing
5322
5323 if (time_text.hour == 0 and time_text.minute > 22 and time_text.minute != 54): #Year: 2023, DOY:029
5324 id_aux = 22
5325 dataOut.DensityFinal[0,:id_aux]=missing
5326 dataOut.EDensityFinal[0,:id_aux]=missing
5327 dataOut.ElecTempFinal[0,:id_aux]=missing
5328 dataOut.EElecTempFinal[0,:id_aux]=missing
5329 dataOut.IonTempFinal[0,:id_aux]=missing
5330 dataOut.EIonTempFinal[0,:id_aux]=missing
5331 dataOut.PhyFinal[0,:id_aux]=missing
5332 dataOut.EPhyFinal[0,:id_aux]=missing
5333
5334 if (time_text.hour == 12 and time_text.minute == 49): #Year: 2023, DOY:029
5335 id_aux = 32
5336 dataOut.DensityFinal[0,id_aux:]=missing
5337 dataOut.EDensityFinal[0,id_aux:]=missing
5338 dataOut.ElecTempFinal[0,id_aux:]=missing
5339 dataOut.EElecTempFinal[0,id_aux:]=missing
5340 dataOut.IonTempFinal[0,id_aux:]=missing
5341 dataOut.EIonTempFinal[0,id_aux:]=missing
5342 dataOut.PhyFinal[0,id_aux:]=missing
5343 dataOut.EPhyFinal[0,id_aux:]=missing
5344
5345 if (time_text.hour == 13 and time_text.minute == 21) or (time_text.hour == 14 and time_text.minute == 25) or (time_text.hour == 14 and time_text.minute == 36) or (time_text.hour == 14 and time_text.minute == 57): #Year: 2023, DOY:029
5346 id_aux = 37
5347 dataOut.DensityFinal[0,id_aux:]=missing
5348 dataOut.EDensityFinal[0,id_aux:]=missing
5349 dataOut.ElecTempFinal[0,id_aux:]=missing
5350 dataOut.EElecTempFinal[0,id_aux:]=missing
5351 dataOut.IonTempFinal[0,id_aux:]=missing
5352 dataOut.EIonTempFinal[0,id_aux:]=missing
5353 dataOut.PhyFinal[0,id_aux:]=missing
5354 dataOut.EPhyFinal[0,id_aux:]=missing
5355
5356 if (time_text.hour == 6 and time_text.minute == 46) or (time_text.hour == 6 and time_text.minute == 57) or (time_text.hour == 7 and time_text.minute == 50) or (time_text.hour == 6 and time_text.minute == 14): #Year: 2023, DOY:029
5357 id_aux = 35
5358 dataOut.DensityFinal[0,id_aux:]=missing
5359 dataOut.EDensityFinal[0,id_aux:]=missing
5360 dataOut.ElecTempFinal[0,id_aux:]=missing
5361 dataOut.EElecTempFinal[0,id_aux:]=missing
5362 dataOut.IonTempFinal[0,id_aux:]=missing
5363 dataOut.EIonTempFinal[0,id_aux:]=missing
5364 dataOut.PhyFinal[0,id_aux:]=missing
5365 dataOut.EPhyFinal[0,id_aux:]=missing
5366
5367 if (time_text.hour == 7 and time_text.minute == 18): #Year: 2023, DOY:029
5368 id_aux = 36
5369 dataOut.DensityFinal[0,id_aux:]=missing
5370 dataOut.EDensityFinal[0,id_aux:]=missing
5371 dataOut.ElecTempFinal[0,id_aux:]=missing
5372 dataOut.EElecTempFinal[0,id_aux:]=missing
5373 dataOut.IonTempFinal[0,id_aux:]=missing
5374 dataOut.EIonTempFinal[0,id_aux:]=missing
5375 dataOut.PhyFinal[0,id_aux:]=missing
5376 dataOut.EPhyFinal[0,id_aux:]=missing
5377
5378 if (time_text.hour == 9 and time_text.minute == 26) or (time_text.hour == 9 and time_text.minute == 37): #Year: 2023, DOY:029
5379 id_aux = 28
5380 dataOut.DensityFinal[0,id_aux:]=missing
5381 dataOut.EDensityFinal[0,id_aux:]=missing
5382 dataOut.ElecTempFinal[0,id_aux:]=missing
5383 dataOut.EElecTempFinal[0,id_aux:]=missing
5384 dataOut.IonTempFinal[0,id_aux:]=missing
5385 dataOut.EIonTempFinal[0,id_aux:]=missing
5386 dataOut.PhyFinal[0,id_aux:]=missing
5387 dataOut.EPhyFinal[0,id_aux:]=missing
5388
5389 if (time_text.hour == 22 and time_text.minute == 14) or (time_text.hour == 22 and time_text.minute == 25): #Year: 2023, DOY:029
5390 id_aux = (37,38)
5391 dataOut.DensityFinal[0,id_aux]=missing
5392 dataOut.EDensityFinal[0,id_aux]=missing
5393 dataOut.ElecTempFinal[0,id_aux]=missing
5394 dataOut.EElecTempFinal[0,id_aux]=missing
5395 dataOut.IonTempFinal[0,id_aux]=missing
5396 dataOut.EIonTempFinal[0,id_aux]=missing
5397 dataOut.PhyFinal[0,id_aux]=missing
5398 dataOut.EPhyFinal[0,id_aux]=missing
5399 '''
5400
5401 '''
5402 if (time_text.hour == 5 and time_text.minute == 32) or (time_text.hour == 5 and time_text.minute == 42) or (time_text.hour == 5 and time_text.minute == 53) or (time_text.hour == 7 and time_text.minute == 18) or (time_text.hour == 7 and time_text.minute == 29) or (time_text.hour == 7 and time_text.minute == 40) or (time_text.hour == 7 and time_text.minute == 50): #Year: 2023, DOY:028
5403 id_aux = 13
5404 dataOut.DensityFinal[0,:id_aux]=missing
5405 dataOut.EDensityFinal[0,:id_aux]=missing
5406 dataOut.ElecTempFinal[0,:id_aux]=missing
5407 dataOut.EElecTempFinal[0,:id_aux]=missing
5408 dataOut.IonTempFinal[0,:id_aux]=missing
5409 dataOut.EIonTempFinal[0,:id_aux]=missing
5410 dataOut.PhyFinal[0,:id_aux]=missing
5411 dataOut.EPhyFinal[0,:id_aux]=missing
5412
5413 if (time_text.hour == 9 and time_text.minute == 5): #Year: 2023, DOY:028
5414 id_aux = 11
5415 dataOut.DensityFinal[0,:id_aux]=missing
5416 dataOut.EDensityFinal[0,:id_aux]=missing
5417 dataOut.ElecTempFinal[0,:id_aux]=missing
5418 dataOut.EElecTempFinal[0,:id_aux]=missing
5419 dataOut.IonTempFinal[0,:id_aux]=missing
5420 dataOut.EIonTempFinal[0,:id_aux]=missing
5421 dataOut.PhyFinal[0,:id_aux]=missing
5422 dataOut.EPhyFinal[0,:id_aux]=missing
5423
5424 if (time_text.hour == 11 and time_text.minute == 13) or (time_text.hour == 12 and time_text.minute == 17): #Year: 2023, DOY:028
5425 id_aux = 34
5426 dataOut.DensityFinal[0,id_aux:]=missing
5427 dataOut.EDensityFinal[0,id_aux:]=missing
5428 dataOut.ElecTempFinal[0,id_aux:]=missing
5429 dataOut.EElecTempFinal[0,id_aux:]=missing
5430 dataOut.IonTempFinal[0,id_aux:]=missing
5431 dataOut.EIonTempFinal[0,id_aux:]=missing
5432 dataOut.PhyFinal[0,id_aux:]=missing
5433 dataOut.EPhyFinal[0,id_aux:]=missing
5434
5435 if (time_text.hour == 7 and time_text.minute == 29): #Year: 2023, DOY:028
5436 id_aux = 35
5437 dataOut.DensityFinal[0,id_aux:]=missing
5438 dataOut.EDensityFinal[0,id_aux:]=missing
5439 dataOut.ElecTempFinal[0,id_aux:]=missing
5440 dataOut.EElecTempFinal[0,id_aux:]=missing
5441 dataOut.IonTempFinal[0,id_aux:]=missing
5442 dataOut.EIonTempFinal[0,id_aux:]=missing
5443 dataOut.PhyFinal[0,id_aux:]=missing
5444 dataOut.EPhyFinal[0,id_aux:]=missing
5445
5446 if (time_text.hour == 13 and time_text.minute == 53) or (time_text.hour == 14 and time_text.minute == 4) or (time_text.hour == 14 and time_text.minute == 36): #Year: 2023, DOY:028
5447 id_aux = 37
5448 dataOut.DensityFinal[0,id_aux:]=missing
5449 dataOut.EDensityFinal[0,id_aux:]=missing
5450 dataOut.ElecTempFinal[0,id_aux:]=missing
5451 dataOut.EElecTempFinal[0,id_aux:]=missing
5452 dataOut.IonTempFinal[0,id_aux:]=missing
5453 dataOut.EIonTempFinal[0,id_aux:]=missing
5454 dataOut.PhyFinal[0,id_aux:]=missing
5455 dataOut.EPhyFinal[0,id_aux:]=missing
5456
5457 if (time_text.hour == 0 and time_text.minute == 22) or (time_text.hour == 0 and time_text.minute == 33): #Year: 2023, DOY:028
5458 id_aux = 22
5459 dataOut.DensityFinal[0,:id_aux]=missing
5460 dataOut.EDensityFinal[0,:id_aux]=missing
5461 dataOut.ElecTempFinal[0,:id_aux]=missing
5462 dataOut.EElecTempFinal[0,:id_aux]=missing
5463 dataOut.IonTempFinal[0,:id_aux]=missing
5464 dataOut.EIonTempFinal[0,:id_aux]=missing
5465 dataOut.PhyFinal[0,:id_aux]=missing
5466 dataOut.EPhyFinal[0,:id_aux]=missing
5467
5468 if (time_text.hour == 23 and time_text.minute == 40) or (time_text.hour == 23 and time_text.minute == 50) or (time_text.hour == 0 and time_text.minute == 12): #Year: 2023, DOY:028
5469 id_aux = 12
5470 dataOut.DensityFinal[0,:id_aux]=missing
5471 dataOut.EDensityFinal[0,:id_aux]=missing
5472 dataOut.ElecTempFinal[0,:id_aux]=missing
5473 dataOut.EElecTempFinal[0,:id_aux]=missing
5474 dataOut.IonTempFinal[0,:id_aux]=missing
5475 dataOut.EIonTempFinal[0,:id_aux]=missing
5476 dataOut.PhyFinal[0,:id_aux]=missing
5477 dataOut.EPhyFinal[0,:id_aux]=missing
5478 '''
5479
5480 '''
5481 #if (time_text.hour >= 7 and time_text.hour < 9): #Year: 2022, DOY:102
5482 if (time_text.hour == 20 and time_text.minute == 8): #Year: 2022, DOY:243
5483 id_aux = 35
5484 dataOut.DensityFinal[0,id_aux:]=missing
5485 dataOut.EDensityFinal[0,id_aux:]=missing
5486 dataOut.ElecTempFinal[0,id_aux:]=missing
5487 dataOut.EElecTempFinal[0,id_aux:]=missing
5488 dataOut.IonTempFinal[0,id_aux:]=missing
5489 dataOut.EIonTempFinal[0,id_aux:]=missing
5490 dataOut.PhyFinal[0,id_aux:]=missing
5491 dataOut.EPhyFinal[0,id_aux:]=missing
5492 if (time_text.hour == 20 and time_text.minute == 19): #Year: 2022, DOY:243
5493 id_aux = 33
5494 dataOut.DensityFinal[0,id_aux:]=missing
5495 dataOut.EDensityFinal[0,id_aux:]=missing
5496 dataOut.ElecTempFinal[0,id_aux:]=missing
5497 dataOut.EElecTempFinal[0,id_aux:]=missing
5498 dataOut.IonTempFinal[0,id_aux:]=missing
5499 dataOut.EIonTempFinal[0,id_aux:]=missing
5500 dataOut.PhyFinal[0,id_aux:]=missing
5501 dataOut.EPhyFinal[0,id_aux:]=missing
5502 if (time_text.hour == 20 and time_text.minute == 29): #Year: 2022, DOY:243
5503 id_aux = 30
5504 dataOut.DensityFinal[0,id_aux:]=missing
5505 dataOut.EDensityFinal[0,id_aux:]=missing
5506 dataOut.ElecTempFinal[0,id_aux:]=missing
5507 dataOut.EElecTempFinal[0,id_aux:]=missing
5508 dataOut.IonTempFinal[0,id_aux:]=missing
5509 dataOut.EIonTempFinal[0,id_aux:]=missing
5510 dataOut.PhyFinal[0,id_aux:]=missing
5511 dataOut.EPhyFinal[0,id_aux:]=missing
5512 if (time_text.hour == 20 and time_text.minute == 44): #Year: 2022, DOY:243
5513 id_aux = 31
5514 dataOut.DensityFinal[0,id_aux:]=missing
5515 dataOut.EDensityFinal[0,id_aux:]=missing
5516 dataOut.ElecTempFinal[0,id_aux:]=missing
5517 dataOut.EElecTempFinal[0,id_aux:]=missing
5518 dataOut.IonTempFinal[0,id_aux:]=missing
5519 dataOut.EIonTempFinal[0,id_aux:]=missing
5520 dataOut.PhyFinal[0,id_aux:]=missing
5521 dataOut.EPhyFinal[0,id_aux:]=missing
5522 if (time_text.hour <= 8): #Year: 2022, DOY:243
5523 id_aux = 11
5524 dataOut.DensityFinal[0,:id_aux]=missing
5525 dataOut.EDensityFinal[0,:id_aux]=missing
5526 dataOut.ElecTempFinal[0,:id_aux]=missing
5527 dataOut.EElecTempFinal[0,:id_aux]=missing
5528 dataOut.IonTempFinal[0,:id_aux]=missing
5529 dataOut.EIonTempFinal[0,:id_aux]=missing
5530 dataOut.PhyFinal[0,:id_aux]=missing
5531 dataOut.EPhyFinal[0,:id_aux]=missing
5532 if (time_text.hour == 23): #Year: 2022, DOY:243
5533 id_aux = 12
5534 dataOut.DensityFinal[0,:id_aux]=missing
5535 dataOut.EDensityFinal[0,:id_aux]=missing
5536 dataOut.ElecTempFinal[0,:id_aux]=missing
5537 dataOut.EElecTempFinal[0,:id_aux]=missing
5538 dataOut.IonTempFinal[0,:id_aux]=missing
5539 dataOut.EIonTempFinal[0,:id_aux]=missing
5540 dataOut.PhyFinal[0,:id_aux]=missing
5541 dataOut.EPhyFinal[0,:id_aux]=missing
5542 if (time_text.hour == 5 and time_text.minute == 21): #Year: 2022, DOY:243
5543 id_aux = (36,37)
5544 dataOut.DensityFinal[0,id_aux]=missing
5545 dataOut.EDensityFinal[0,id_aux]=missing
5546 dataOut.ElecTempFinal[0,id_aux]=missing
5547 dataOut.EElecTempFinal[0,id_aux]=missing
5548 dataOut.IonTempFinal[0,id_aux]=missing
5549 dataOut.EIonTempFinal[0,id_aux]=missing
5550 dataOut.PhyFinal[0,id_aux]=missing
5551 dataOut.EPhyFinal[0,id_aux]=missing
5552 if (time_text.hour == 5 and time_text.minute == 53): #Year: 2022, DOY:243
5553 id_aux = (37,38)
5554 dataOut.DensityFinal[0,id_aux]=missing
5555 dataOut.EDensityFinal[0,id_aux]=missing
5556 dataOut.ElecTempFinal[0,id_aux]=missing
5557 dataOut.EElecTempFinal[0,id_aux]=missing
5558 dataOut.IonTempFinal[0,id_aux]=missing
5559 dataOut.EIonTempFinal[0,id_aux]=missing
5560 dataOut.PhyFinal[0,id_aux]=missing
5561 dataOut.EPhyFinal[0,id_aux]=missing
5562 if (time_text.hour == 6 and time_text.minute == 4): #Year: 2022, DOY:243
5563 id_aux = (38,39)
5564 dataOut.DensityFinal[0,id_aux]=missing
5565 dataOut.EDensityFinal[0,id_aux]=missing
5566 dataOut.ElecTempFinal[0,id_aux]=missing
5567 dataOut.EElecTempFinal[0,id_aux]=missing
5568 dataOut.IonTempFinal[0,id_aux]=missing
5569 dataOut.EIonTempFinal[0,id_aux]=missing
5570 dataOut.PhyFinal[0,id_aux]=missing
5571 dataOut.EPhyFinal[0,id_aux]=missing
5572 if (time_text.hour == 12 and time_text.minute == 6): #Year: 2022, DOY:243
5573 id_aux = (29,30)
5574 dataOut.DensityFinal[0,id_aux]=missing
5575 dataOut.EDensityFinal[0,id_aux]=missing
5576 dataOut.ElecTempFinal[0,id_aux]=missing
5577 dataOut.EElecTempFinal[0,id_aux]=missing
5578 dataOut.IonTempFinal[0,id_aux]=missing
5579 dataOut.EIonTempFinal[0,id_aux]=missing
5580 dataOut.PhyFinal[0,id_aux]=missing
5581 dataOut.EPhyFinal[0,id_aux]=missing
5582 if (time_text.hour == 14 and time_text.minute == 14): #Year: 2022, DOY:243
5583 id_aux = (35,36)
5584 dataOut.DensityFinal[0,id_aux]=missing
5585 dataOut.EDensityFinal[0,id_aux]=missing
5586 dataOut.ElecTempFinal[0,id_aux]=missing
5587 dataOut.EElecTempFinal[0,id_aux]=missing
5588 dataOut.IonTempFinal[0,id_aux]=missing
5589 dataOut.EIonTempFinal[0,id_aux]=missing
5590 dataOut.PhyFinal[0,id_aux]=missing
5591 dataOut.EPhyFinal[0,id_aux]=missing
5592 if (time_text.hour == 23 and time_text.minute == 2): #Year: 2022, DOY:243
5593 id_aux = (41,42)
5594 dataOut.DensityFinal[0,id_aux]=missing
5595 dataOut.EDensityFinal[0,id_aux]=missing
5596 dataOut.ElecTempFinal[0,id_aux]=missing
5597 dataOut.EElecTempFinal[0,id_aux]=missing
5598 dataOut.IonTempFinal[0,id_aux]=missing
5599 dataOut.EIonTempFinal[0,id_aux]=missing
5600 dataOut.PhyFinal[0,id_aux]=missing
5601 dataOut.EPhyFinal[0,id_aux]=missing
5602 if (time_text.hour == 0 and time_text.minute == 8): #Year: 2022, DOY:243
5603 id_aux = 33
5604 dataOut.DensityFinal[0,id_aux:]=missing
5605 dataOut.EDensityFinal[0,id_aux:]=missing
5606 dataOut.ElecTempFinal[0,id_aux:]=missing
5607 dataOut.EElecTempFinal[0,id_aux:]=missing
5608 dataOut.IonTempFinal[0,id_aux:]=missing
5609 dataOut.EIonTempFinal[0,id_aux:]=missing
5610 dataOut.PhyFinal[0,id_aux:]=missing
5611 dataOut.EPhyFinal[0,id_aux:]=missing
5612 id_aux = 18
5613 dataOut.DensityFinal[0,id_aux]=missing
4978 dataOut.EDensityFinal[0,id_aux]=missing
5614 dataOut.EDensityFinal[0,id_aux]=missing
4979 dataOut.ElecTempFinal[0,id_aux]=missing
5615 dataOut.ElecTempFinal[0,id_aux]=missing
4980 dataOut.EElecTempFinal[0,id_aux]=missing
5616 dataOut.EElecTempFinal[0,id_aux]=missing
@@ -5371,7 +6007,7 class DataSaveCleaner(Operation):
5371 dataOut.PhyFinal[0,id_aux:]=missing
6007 dataOut.PhyFinal[0,id_aux:]=missing
5372 dataOut.EPhyFinal[0,id_aux:]=missing
6008 dataOut.EPhyFinal[0,id_aux:]=missing
5373 '''
6009 '''
5374 #'''
6010 '''
5375 #print("den: ", dataOut.DensityFinal[0,27])
6011 #print("den: ", dataOut.DensityFinal[0,27])
5376 if (time_text.hour == 5 and time_text.minute == 42): #Year: 2022, DOY:242
6012 if (time_text.hour == 5 and time_text.minute == 42): #Year: 2022, DOY:242
5377 id_aux = 16
6013 id_aux = 16
@@ -5493,13 +6129,13 class DataSaveCleaner(Operation):
5493 dataOut.EIonTempFinal[0,id_aux]=missing
6129 dataOut.EIonTempFinal[0,id_aux]=missing
5494 dataOut.PhyFinal[0,id_aux]=missing
6130 dataOut.PhyFinal[0,id_aux]=missing
5495 dataOut.EPhyFinal[0,id_aux]=missing
6131 dataOut.EPhyFinal[0,id_aux]=missing
5496 #'''
6132 '''
5497 #print("den_final",dataOut.DensityFinal)
6133 #print("den_final",dataOut.DensityFinal)
5498
6134
5499
6135
5500 dataOut.flagNoData = numpy.all(numpy.isnan(dataOut.DensityFinal)) #Si todos los valores son NaN no se prosigue
6136 dataOut.flagNoData = numpy.all(numpy.isnan(dataOut.DensityFinal)) #Si todos los valores son NaN no se prosigue
5501
6137
5502 ####dataOut.flagNoData = False #Solo para ploteo
6138 #dataOut.flagNoData = False #Descomentar solo para ploteo #Comentar para MADWriter
5503
6139
5504 dataOut.DensityFinal *= 1.e6 #Convert units to m^⁻3
6140 dataOut.DensityFinal *= 1.e6 #Convert units to m^⁻3
5505 dataOut.EDensityFinal *= 1.e6 #Convert units to m^⁻3
6141 dataOut.EDensityFinal *= 1.e6 #Convert units to m^⁻3
@@ -6319,7 +6955,7 class removeDCHAE(Operation):
6319
6955
6320 return dataOut
6956 return dataOut
6321
6957
6322 class Decoder(Operation):
6958 class Decoder_Original(Operation):
6323
6959
6324 isConfig = False
6960 isConfig = False
6325 __profIndex = 0
6961 __profIndex = 0
@@ -6355,10 +6991,13 class Decoder(Operation):
6355 self.__nChannels = dataOut.nChannels
6991 self.__nChannels = dataOut.nChannels
6356 self.__nProfiles = dataOut.nProfiles
6992 self.__nProfiles = dataOut.nProfiles
6357 self.__nHeis = dataOut.nHeights
6993 self.__nHeis = dataOut.nHeights
6358
6994 #print("self.__nHeis: ", self.__nHeis)
6995 #print("self.nBaud: ", self.nBaud)
6996 #exit(1)
6359 if self.__nHeis < self.nBaud:
6997 if self.__nHeis < self.nBaud:
6360 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
6998 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
6361
6999 #print("JJE")
7000 #exit(1)
6362 #Frequency
7001 #Frequency
6363 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
7002 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
6364
7003
@@ -6421,6 +7060,10 class Decoder(Operation):
6421 #print(profilesList)
7060 #print(profilesList)
6422 for i in range(self.__nChannels):
7061 for i in range(self.__nChannels):
6423 for j in profilesList:
7062 for j in profilesList:
7063 #print("data.shape: ", data.shape)
7064 #print("code_block: ", code_block.shape)
7065 #print("corr: ", numpy.shape(numpy.correlate(data[i,j,:], code_block[j,:], mode='full')))
7066 #exit(1)
6424 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
7067 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
6425 return self.datadecTime
7068 return self.datadecTime
6426
7069
@@ -6444,7 +7087,8 class Decoder(Operation):
6444
7087
6445 if dataOut.flagDecodeData:
7088 if dataOut.flagDecodeData:
6446 print("This data is already decoded, recoding again ...")
7089 print("This data is already decoded, recoding again ...")
6447
7090 #print("code: ", numpy.shape(code))
7091 #exit(1)
6448 if not self.isConfig:
7092 if not self.isConfig:
6449
7093
6450 if code is None:
7094 if code is None:
@@ -6522,6 +7166,234 class Decoder(Operation):
6522 return dataOut
7166 return dataOut
6523 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
7167 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
6524
7168
7169 class Decoder(Operation):
7170
7171 isConfig = False
7172 __profIndex = 0
7173
7174 code = None
7175
7176 nCode = None
7177 nBaud = None
7178
7179 def __init__(self, **kwargs):
7180
7181 Operation.__init__(self, **kwargs)
7182
7183 self.times = None
7184 self.osamp = None
7185 # self.__setValues = False
7186 self.isConfig = False
7187 self.setupReq = False
7188 def setup(self, code, osamp, dataOut):
7189
7190 self.__profIndex = 0
7191
7192 self.code = code
7193
7194 self.nCode = len(code)
7195 #self.nBaud = len(code[0])
7196 self.nBaud = int(numpy.shape(code)[-1])
7197
7198 if (osamp != None) and (osamp >1):
7199 self.osamp = osamp
7200 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
7201 self.nBaud = self.nBaud*self.osamp
7202
7203 self.__nChannels = dataOut.nChannels
7204 self.__nProfiles = dataOut.nProfiles
7205 self.__nHeis = dataOut.nHeights
7206 #print("self.__nHeis: ", self.__nHeis)
7207 #print("self.nBaud: ", self.nBaud)
7208 #exit(1)
7209 if self.__nHeis < self.nBaud:
7210 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
7211 #print("JJE")
7212 #exit(1)
7213 '''
7214 #Frequency
7215 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
7216
7217 __codeBuffer[:,0:self.nBaud] = self.code
7218
7219 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
7220 '''
7221 if dataOut.flagDataAsBlock:
7222
7223 self.ndatadec = self.__nHeis #- self.nBaud + 1
7224
7225 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
7226
7227 else:
7228
7229 #Time
7230 self.ndatadec = self.__nHeis #- self.nBaud + 1
7231
7232
7233 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
7234
7235 def __convolutionInFreq(self, data):
7236
7237 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
7238
7239 fft_data = numpy.fft.fft(data, axis=1)
7240
7241 conv = fft_data*fft_code
7242
7243 data = numpy.fft.ifft(conv,axis=1)
7244
7245 return data
7246
7247 def __convolutionInFreqOpt(self, data):
7248
7249 raise NotImplementedError
7250
7251 def __convolutionInTime(self, data):
7252
7253 code = self.code[self.__profIndex]
7254 for i in range(self.__nChannels):
7255 #aux=numpy.correlate(data[i,:], code, mode='full')
7256 #print(numpy.shape(aux))
7257 #print(numpy.shape(data[i,:]))
7258 #print(numpy.shape(code))
7259 #exit(1)
7260 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
7261
7262 return self.datadecTime
7263
7264 def __convolutionByBlockInTime(self, data, AutoDecod):
7265
7266 if not AutoDecod:
7267 repetitions = int(self.__nProfiles / self.nCode)
7268 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
7269 junk = junk.flatten()
7270 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
7271 else:
7272 code_block = self.code
7273 profilesList = range(self.__nProfiles)
7274 #print(numpy.shape(self.datadecTime))
7275 #print(numpy.shape(data))
7276 #print(profilesList)
7277 for i in range(self.__nChannels):
7278 for j in profilesList:
7279 #print("data.shape: ", data.shape)
7280 #print("code_block: ", code_block.shape)
7281 #print("corr: ", numpy.shape(numpy.correlate(data[i,j,:], code_block[j,:], mode='full')))
7282 #exit(1)
7283 if not AutoDecod:
7284 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
7285 else:
7286 if i%2:
7287 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[0,j,:], mode='full')[self.nBaud-1:]
7288 else:
7289 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[1,j,:], mode='full')[self.nBaud-1:]
7290
7291 return self.datadecTime
7292
7293 def __convolutionByBlockInFreq(self, data):
7294
7295 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
7296
7297
7298 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
7299
7300 fft_data = numpy.fft.fft(data, axis=2)
7301
7302 conv = fft_data*fft_code
7303
7304 data = numpy.fft.ifft(conv,axis=2)
7305
7306 return data
7307
7308
7309 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None, AutoDecod = 0):
7310
7311 if dataOut.flagDecodeData:
7312 print("This data is already decoded, recoding again ...")
7313 #print("code: ", numpy.shape(code))
7314 #exit(1)
7315 if not self.isConfig:
7316
7317 if code is None and not AutoDecod:
7318 if dataOut.code is None:
7319 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
7320
7321 code = dataOut.code
7322 else:
7323 if not AutoDecod:
7324 code = numpy.array(code).reshape(nCode,nBaud)
7325 else:
7326 po = 0
7327 print("***********AutoDecod***********")
7328 code = dataOut.data[:2,:,po:po+64]
7329 #print("AutoDecod Shape: ", numpy.shape(code))
7330 #exit(1)
7331 self.setup(code, osamp, dataOut)
7332
7333 self.isConfig = True
7334
7335 if mode == 3:
7336 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
7337
7338 if times != None:
7339 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
7340
7341 if self.code is None:
7342 print("Fail decoding: Code is not defined.")
7343 return
7344
7345 self.__nProfiles = dataOut.nProfiles
7346 datadec = None
7347
7348 if mode == 3:
7349 mode = 0
7350
7351 if dataOut.flagDataAsBlock:
7352 """
7353 Decoding when data have been read as block,
7354 """
7355
7356 if mode == 0:
7357 datadec = self.__convolutionByBlockInTime(dataOut.data, AutoDecod)
7358 if mode == 1:
7359 datadec = self.__convolutionByBlockInFreq(dataOut.data)
7360 else:
7361 """
7362 Decoding when data have been read profile by profile
7363 """
7364 if mode == 0:
7365 datadec = self.__convolutionInTime(dataOut.data)
7366
7367 if mode == 1:
7368 datadec = self.__convolutionInFreq(dataOut.data)
7369
7370 if mode == 2:
7371 datadec = self.__convolutionInFreqOpt(dataOut.data)
7372
7373 if datadec is None:
7374 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
7375
7376 dataOut.code = self.code
7377 dataOut.nCode = self.nCode
7378 dataOut.nBaud = self.nBaud
7379
7380 dataOut.data = datadec#/self.nBaud
7381
7382 #print("before",dataOut.heightList)
7383 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
7384 #print("after",dataOut.heightList)
7385
7386 dataOut.flagDecodeData = True #asumo q la data esta decodificada
7387
7388 if self.__profIndex == self.nCode-1:
7389 self.__profIndex = 0
7390 return dataOut
7391
7392 self.__profIndex += 1
7393
7394 return dataOut
7395
7396
6525 class DecoderRoll(Operation):
7397 class DecoderRoll(Operation):
6526
7398
6527 isConfig = False
7399 isConfig = False
@@ -6883,7 +7755,19 class ProfileSelector(Operation):
6883 dataOut.nProfiles = len(profileList)
7755 dataOut.nProfiles = len(profileList)
6884 dataOut.profileIndex = dataOut.nProfiles - 1
7756 dataOut.profileIndex = dataOut.nProfiles - 1
6885 dataOut.flagNoData = False
7757 dataOut.flagNoData = False
6886
7758 #print("Shape after prof select: ", numpy.shape(dataOut.data))
7759 #print(dataOut.heightList)
7760 #exit(1)
7761 '''
7762 po = 5
7763 import matplotlib.pyplot as plt
7764 this_data = dataOut.data[0,0,:]
7765 plt.plot(this_data*numpy.conjugate(this_data),marker='*',linestyle='-')
7766 plt.axvline(po)
7767 plt.axvline(po+64-1)
7768 plt.grid()
7769 plt.show()
7770 '''
6887 return dataOut
7771 return dataOut
6888
7772
6889 """
7773 """
@@ -7520,7 +8404,7 class CrossProdHybrid(CrossProdDP):
7520 #exit(1)
8404 #exit(1)
7521 if i==0:
8405 if i==0:
7522 self.lagp0[n][j][self.bcounter-1]=numpy.sum(c[:dataOut.NSCAN])
8406 self.lagp0[n][j][self.bcounter-1]=numpy.sum(c[:dataOut.NSCAN])
7523 ww[n,j,:,self.bcounter-1]=c[:dataOut.NSCAN]
8407 #ww[n,j,:,self.bcounter-1]=c[:dataOut.NSCAN]
7524 self.lagp3[n][j][self.bcounter-1]=numpy.sum(c[dataOut.NSCAN:]/self.cnorm)
8408 self.lagp3[n][j][self.bcounter-1]=numpy.sum(c[dataOut.NSCAN:]/self.cnorm)
7525 elif i==1:
8409 elif i==1:
7526 self.lagp1[n][j][self.bcounter-1]=numpy.sum(c[:dataOut.NSCAN])
8410 self.lagp1[n][j][self.bcounter-1]=numpy.sum(c[:dataOut.NSCAN])
@@ -7540,8 +8424,8 class CrossProdHybrid(CrossProdDP):
7540 #print(self.cnorm)
8424 #print(self.cnorm)
7541 #exit(1)
8425 #exit(1)
7542 #print("self,lagp0: ", self.lagp0[0,0,self.bcounter-1])
8426 #print("self,lagp0: ", self.lagp0[0,0,self.bcounter-1])
7543 print(ww[:,0,0,self.bcounter-1])
8427 #print(ww[:,0,0,self.bcounter-1])
7544 exit(1)
8428 #exit(1)
7545
8429
7546
8430
7547 def LP_median_estimates_original(self,dataOut):
8431 def LP_median_estimates_original(self,dataOut):
@@ -7998,9 +8882,9 class LongPulseAnalysis(Operation):
7998 #print(anoise0)
8882 #print(anoise0)
7999 #print(anoise1)
8883 #print(anoise1)
8000 #exit(1)
8884 #exit(1)
8001 print("nis: ", dataOut.nis)
8885 #print("nis: ", dataOut.nis)
8002 print("pan: ", dataOut.pan)
8886 #print("pan: ", dataOut.pan)
8003 print("pbn: ", dataOut.pbn)
8887 #print("pbn: ", dataOut.pbn)
8004 #print(numpy.sum(dataOut.output_LP_integrated[0,:,0]))
8888 #print(numpy.sum(dataOut.output_LP_integrated[0,:,0]))
8005 '''
8889 '''
8006 import matplotlib.pyplot as plt
8890 import matplotlib.pyplot as plt
@@ -8008,8 +8892,8 class LongPulseAnalysis(Operation):
8008 plt.show()
8892 plt.show()
8009 '''
8893 '''
8010 #print(dataOut.output_LP_integrated[0,40,0])
8894 #print(dataOut.output_LP_integrated[0,40,0])
8011 print(numpy.sum(dataOut.output_LP_integrated[:,0,0]))
8895 #print(numpy.sum(dataOut.output_LP_integrated[:,0,0]))
8012 exit(1)
8896 #exit(1)
8013
8897
8014 #################### PROBAR MÁS INTEGRACIÓN, SINO MODIFICAR VALOR DE "NIS" ####################
8898 #################### PROBAR MÁS INTEGRACIÓN, SINO MODIFICAR VALOR DE "NIS" ####################
8015 # VER dataOut.nProfiles_LP #
8899 # VER dataOut.nProfiles_LP #
@@ -8237,7 +9121,7 class LongPulseAnalysis(Operation):
8237 print(numpy.sum(dataOut.te2))
9121 print(numpy.sum(dataOut.te2))
8238 exit(1)
9122 exit(1)
8239 '''
9123 '''
8240 print("Success")
9124 #print("Success 1")
8241 ###################Correlation pulse and itself
9125 ###################Correlation pulse and itself
8242
9126
8243 #print(dataOut.NRANGE)
9127 #print(dataOut.NRANGE)
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
NO CONTENT: file was removed
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
NO CONTENT: file was removed
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
NO CONTENT: file was removed
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
NO CONTENT: file was removed
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
NO CONTENT: file was removed
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
NO CONTENT: file was removed
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
NO CONTENT: file was removed
General Comments 0
You need to be logged in to leave comments. Login now