| @@ -548,6 +548,9 class Spectra(JROData): | |||
|
|
548 | 548 | |
|
|
549 | 549 | if self.flagDecodeData: |
|
|
550 | 550 | pwcode = numpy.sum(self.code[0]**2) |
|
|
551 | #pwcode = 64 | |
|
|
552 | #print("pwcode: ", pwcode) | |
|
|
553 | #exit(1) | |
|
|
551 | 554 | #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter |
|
|
552 | 555 | normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter |
|
|
553 | 556 | |
| @@ -42,6 +42,9 class SpectraPlot(Plot): | |||
|
|
42 | 42 | meta = {} |
|
|
43 | 43 | |
|
|
44 | 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 | 48 | #print("Spc: ",spc[0]) |
|
|
46 | 49 | #exit(1) |
|
|
47 | 50 | data['spc'] = spc |
| @@ -716,14 +719,15 class RTIPlot(Plot): | |||
|
|
716 | 719 | self.x = self.data.times |
|
|
717 | 720 | self.y = self.data.yrange |
|
|
718 | 721 | self.z = self.data[self.CODE] |
|
|
719 | ||
|
|
722 | #print("Inside RTI: ", self.z) | |
|
|
720 | 723 | self.z = numpy.ma.masked_invalid(self.z) |
|
|
721 | 724 | |
|
|
722 | 725 | if self.decimation is None: |
|
|
723 | 726 | x, y, z = self.fill_gaps(self.x, self.y, self.z) |
|
|
724 | 727 | else: |
|
|
725 | 728 | x, y, z = self.fill_gaps(*self.decimate()) |
|
|
726 | ||
|
|
729 | #print("self.z: ", self.z) | |
|
|
730 | #exit(1) | |
|
|
727 | 731 | ''' |
|
|
728 | 732 | if not isinstance(self.zmin, collections.abc.Sequence): |
|
|
729 | 733 | if not self.zmin: |
| @@ -673,9 +673,9 class EDensityPlot(Plot): | |||
|
|
673 | 673 | data['den_Faraday'] = dataOut.dphi[:dataOut.NSHTS] |
|
|
674 | 674 | data['den_error'] = dataOut.sdp2[:dataOut.NSHTS] |
|
|
675 | 675 | #data['err_Faraday'] = dataOut.sdn1[:dataOut.NSHTS] |
|
|
676 | print(numpy.shape(data['den_power'])) | |
|
|
677 | print(numpy.shape(data['den_Faraday'])) | |
|
|
678 | print(numpy.shape(data['den_error'])) | |
|
|
676 | #print(numpy.shape(data['den_power'])) | |
|
|
677 | #print(numpy.shape(data['den_Faraday'])) | |
|
|
678 | #print(numpy.shape(data['den_error'])) | |
|
|
679 | 679 | |
|
|
680 | 680 | data['NSHTS'] = dataOut.NSHTS |
|
|
681 | 681 | |
| @@ -49,6 +49,7 DEF_HEADER = { | |||
|
|
49 | 49 | MNEMONICS = { |
|
|
50 | 50 | 10: 'jro', |
|
|
51 | 51 | 11: 'jbr', |
|
|
52 | 14: 'jmp', #Added by R. Flores | |
|
|
52 | 53 | 840: 'jul', |
|
|
53 | 54 | 13: 'jas', |
|
|
54 | 55 | 1000: 'pbr', |
| 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 | 1152 | #return A1f, B1f, C1f, A2f, B2f, C2f, K2f, Df, doppler |
|
|
1153 | 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 | 1157 | from scipy.optimize import least_squares |
|
|
1158 | 1158 | |
| @@ -1176,7 +1176,7 class Oblique_Gauss_Fit(Operation): | |||
|
|
1176 | 1176 | #print(a1,b1,c1,a2,b2,c2,k2,d) |
|
|
1177 | 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 | 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,- |
|
|
|
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 | 1181 | #print(bounds) |
|
|
1182 | 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 | 1194 | x0_value = numpy.array([spc_max,dop1_x0,30,-.1,spc_max/4, dop2_x0,150,1,1.0e7]) |
|
|
1195 | 1195 | #x0_value = numpy.array([spc_max,-400.5,30,-.1,spc_max/4,-200.5,150,1,1.0e7]) |
|
|
1196 | 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 | 1207 | popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0) |
|
|
1198 | 1208 | # popt = least_squares(lsq_func,[a1,b1,c1,a2,b2,c2,k2,d],x_scale=params_scale,verbose=1) |
|
|
1199 | 1209 | #print(popt) |
|
|
1200 | ||
|
|
1210 | #########print("INSIDE 2") | |
|
|
1201 | 1211 | J = popt.jac |
|
|
1202 | 1212 | |
|
|
1203 | 1213 | try: |
| @@ -1223,6 +1233,8 class Oblique_Gauss_Fit(Operation): | |||
|
|
1223 | 1233 | doppler2 = freq[numpy.argmax(aux2)] |
|
|
1224 | 1234 | #print("error",error) |
|
|
1225 | 1235 | #exit(1) |
|
|
1236 | ||
|
|
1237 | ||
|
|
1226 | 1238 | return A1f, B1f, C1f, K1f, A2f, B2f, C2f, K2f, Df, doppler1, doppler2, error |
|
|
1227 | 1239 | |
|
|
1228 | 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 | 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 | 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 | 1658 | pass |
|
|
1566 | 1659 | |
|
|
1567 | 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 | 1667 | for hei in itertools.chain(l1, l2): |
|
|
1570 | 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 | 1672 | continue #Avoids the analysis when there is only noise |
|
|
1573 | 1673 | |
|
|
1574 | 1674 | try: |
| @@ -1592,8 +1692,11 class Oblique_Gauss_Fit(Operation): | |||
|
|
1592 | 1692 | from scipy.signal import medfilt |
|
|
1593 | 1693 | spcm = medfilt(spc,11) |
|
|
1594 | 1694 | if x[numpy.argmax(spcm)] <= 0: |
|
|
1595 | #print("EEJ") | |
|
|
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) | |
|
|
1695 | #print("EEJ", dataOut.heightList[hei], hei) | |
|
|
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 | 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 | 1701 | # dataOut.Oblique_params[0,:,hei] *= numpy.NAN |
|
|
1599 | 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 | 1775 | dataOut.paramInterval = dataOut.nProfiles*dataOut.nCohInt*dataOut.ippSeconds |
|
|
1673 | 1776 | dataOut.lat=-11.95 |
|
|
1674 | 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 | 1783 | if mode == 9: #Double Skew Gaussian |
|
|
1677 | 1784 | #dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,-2,:] #Pos[Max_value] |
|
|
1678 | 1785 | #dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,1,:] #Shift |
| @@ -1703,9 +1810,154 class Oblique_Gauss_Fit(Operation): | |||
|
|
1703 | 1810 | dataOut.Err_Dop_EEJ_T2 = dataOut.Oblique_param_errors[:,4,:] |
|
|
1704 | 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 | 1828 | dataOut.mode = mode |
|
|
1707 | 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 | 1962 | return dataOut |
|
|
1711 | 1963 | |
| @@ -2547,7 +2799,7 class SpectralFitting(Operation): | |||
|
|
2547 | 2799 | #def __DiffCoherent(self,snrth, spectra, cspectra, nProf, heights,nChan, nHei, nPairs, channels, noise, crosspairs): |
|
|
2548 | 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 | 2803 | nProf = dataOut.nProfiles |
|
|
2552 | 2804 | heights = dataOut.heightList |
|
|
2553 | 2805 | nHei = len(heights) |
| @@ -2573,20 +2825,19 class SpectralFitting(Operation): | |||
|
|
2573 | 2825 | |
|
|
2574 | 2826 | if coh_th == None : coh_th = numpy.array([0.75,0.65,0.15]) # 0.65 |
|
|
2575 | 2827 | if hei_th == None : hei_th = numpy.array([60,300,650]) |
|
|
2576 |
for ic in range( |
|
|
|
2828 | for ic in range(nPairs): | |
|
|
2577 | 2829 | pair = crosspairs[ic] |
|
|
2578 | 2830 | #si el SNR es mayor que el SNR threshold los datos se toman coherentes |
|
|
2579 | 2831 | s_n0 = power[pair[0],:]/noise[pair[0]] |
|
|
2580 | 2832 | s_n1 = power[pair[1],:]/noise[pair[1]] |
|
|
2581 | ||
|
|
2582 | 2833 | valid1 =(s_n0>=snr_th).nonzero() |
|
|
2583 | 2834 | valid2 = (s_n1>=snr_th).nonzero() |
|
|
2584 | #valid = valid2 + valid1 #numpy.concatenate((valid1,valid2), axis=None) | |
|
|
2835 | ||
|
|
2585 | 2836 | valid1 = numpy.array(valid1[0]) |
|
|
2586 | 2837 | valid2 = numpy.array(valid2[0]) |
|
|
2587 | 2838 | valid = valid1 |
|
|
2588 | 2839 | for iv in range(len(valid2)): |
|
|
2589 | #for ivv in range(len(valid1)) : | |
|
|
2840 | ||
|
|
2590 | 2841 | indv = numpy.array((valid1 == valid2[iv]).nonzero()) |
|
|
2591 | 2842 | if len(indv[0]) == 0 : |
|
|
2592 | 2843 | valid = numpy.concatenate((valid,valid2[iv]), axis=None) |
| @@ -2594,17 +2845,16 class SpectralFitting(Operation): | |||
|
|
2594 | 2845 | my_coh_aver[pair[0],valid]=1 |
|
|
2595 | 2846 | my_coh_aver[pair[1],valid]=1 |
|
|
2596 | 2847 | # si la coherencia es mayor a la coherencia threshold los datos se toman |
|
|
2597 | #print my_coh_aver[0,:] | |
|
|
2848 | ||
|
|
2598 | 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 | 2851 | for ih in range(len(hei_th)): |
|
|
2601 | 2852 | hvalid = (heights>hei_th[ih]).nonzero() |
|
|
2602 | 2853 | hvalid = hvalid[0] |
|
|
2603 | 2854 | if len(hvalid)>0: |
|
|
2604 | 2855 | valid = (numpy.absolute(coh[hvalid])>coh_th[ih]).nonzero() |
|
|
2605 | 2856 | valid = valid[0] |
|
|
2606 | #print('hvalid:',hvalid) | |
|
|
2607 | #print('valid', valid) | |
|
|
2857 | ||
|
|
2608 | 2858 | if len(valid)>0: |
|
|
2609 | 2859 | my_coh_aver[pair[0],hvalid[valid]] =1 |
|
|
2610 | 2860 | my_coh_aver[pair[1],hvalid[valid]] =1 |
| @@ -2620,7 +2870,7 class SpectralFitting(Operation): | |||
|
|
2620 | 2870 | my_incoh_aver[pair[1],incoh_echoes] = 1 |
|
|
2621 | 2871 | |
|
|
2622 | 2872 | |
|
|
2623 |
for ic in range( |
|
|
|
2873 | for ic in range(nPairs): | |
|
|
2624 | 2874 | pair = crosspairs[ic] |
|
|
2625 | 2875 | |
|
|
2626 | 2876 | valid1 =(my_coh_aver[pair[0],:]==1 ).nonzero() |
| @@ -2628,29 +2878,25 class SpectralFitting(Operation): | |||
|
|
2628 | 2878 | valid1 = numpy.array(valid1[0]) |
|
|
2629 | 2879 | valid2 = numpy.array(valid2[0]) |
|
|
2630 | 2880 | valid = valid1 |
|
|
2631 | #print valid1 , valid2 | |
|
|
2881 | ||
|
|
2632 | 2882 | for iv in range(len(valid2)): |
|
|
2633 | #for ivv in range(len(valid1)) : | |
|
|
2883 | ||
|
|
2634 | 2884 | indv = numpy.array((valid1 == valid2[iv]).nonzero()) |
|
|
2635 | 2885 | if len(indv[0]) == 0 : |
|
|
2636 | 2886 | valid = numpy.concatenate((valid,valid2[iv]), axis=None) |
|
|
2637 | #print valid | |
|
|
2638 | #valid = numpy.concatenate((valid1,valid2), axis=None) | |
|
|
2639 | 2887 | valid1 =(my_coh_aver[pair[0],:] !=1 ).nonzero() |
|
|
2640 | 2888 | valid2 = (my_coh_aver[pair[1],:] !=1).nonzero() |
|
|
2641 | 2889 | valid1 = numpy.array(valid1[0]) |
|
|
2642 | 2890 | valid2 = numpy.array(valid2[0]) |
|
|
2643 | 2891 | incoh_echoes = valid1 |
|
|
2644 | #print valid1, valid2 | |
|
|
2645 | #incoh_echoes= numpy.concatenate((valid1,valid2), axis=None) | |
|
|
2892 | ||
|
|
2646 | 2893 | for iv in range(len(valid2)): |
|
|
2647 | #for ivv in range(len(valid1)) : | |
|
|
2894 | ||
|
|
2648 | 2895 | indv = numpy.array((valid1 == valid2[iv]).nonzero()) |
|
|
2649 | 2896 | if len(indv[0]) == 0 : |
|
|
2650 | 2897 | incoh_echoes = numpy.concatenate(( incoh_echoes,valid2[iv]), axis=None) |
|
|
2651 | #print incoh_echoes | |
|
|
2898 | ||
|
|
2652 | 2899 | if len(valid)>0: |
|
|
2653 | #print pair | |
|
|
2654 | 2900 | coh_spectra[pair[0],:,valid] = spectra[pair[0],:,valid] |
|
|
2655 | 2901 | coh_spectra[pair[1],:,valid] = spectra[pair[1],:,valid] |
|
|
2656 | 2902 | coh_cspectra[ic,:,valid] = cspectra[ic,:,valid] |
| @@ -2662,20 +2908,12 class SpectralFitting(Operation): | |||
|
|
2662 | 2908 | incoh_cspectra[ic,:,incoh_echoes] = cspectra[ic,:,incoh_echoes] |
|
|
2663 | 2909 | incoh_aver[pair[0],incoh_echoes]=1 |
|
|
2664 | 2910 | incoh_aver[pair[1],incoh_echoes]=1 |
|
|
2665 | #plt.imshow(spectra[0,:,:],vmin=20000000) | |
|
|
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 | |
|
|
2911 | ||
|
|
2674 | 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 | 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 | 2917 | nProf = dataOut.nProfiles |
|
|
2680 | 2918 | heights = dataOut.heightList |
|
|
2681 | 2919 | nHei = len(heights) |
| @@ -2684,15 +2922,9 class SpectralFitting(Operation): | |||
|
|
2684 | 2922 | crosspairs = dataOut.groupList |
|
|
2685 | 2923 | nPairs = len(crosspairs) |
|
|
2686 | 2924 | |
|
|
2687 | #data = dataOut.data_pre[0] | |
|
|
2688 | 2925 | absc = dataOut.abscissaList[:-1] |
|
|
2689 | #noise = dataOut.noise | |
|
|
2690 | #nChannel = data.shape[0] | |
|
|
2691 | 2926 | data_param = numpy.zeros((nChan, 4, spectra.shape[2])) |
|
|
2692 | 2927 | |
|
|
2693 | ||
|
|
2694 | #plt.plot(absc) | |
|
|
2695 | #plt.show() | |
|
|
2696 | 2928 | clean_coh_spectra = spectra.copy() |
|
|
2697 | 2929 | clean_coh_cspectra = cspectra.copy() |
|
|
2698 | 2930 | clean_coh_aver = coh_aver.copy() |
| @@ -2703,17 +2935,14 class SpectralFitting(Operation): | |||
|
|
2703 | 2935 | rtime0 = [6,18] # periodo sin ESF |
|
|
2704 | 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 | 2939 | if clean_coh_echoes == 1 : |
|
|
2708 | 2940 | for ind in range(nChan): |
|
|
2709 | 2941 | data_param[ind,:,:] = self.__calculateMoments( spectra[ind,:,:] , absc , noise[ind] ) |
|
|
2710 | #print data_param[:,3] | |
|
|
2942 | ||
|
|
2711 | 2943 | spwd = data_param[:,3] |
|
|
2712 | #print spwd.shape | |
|
|
2944 | ||
|
|
2713 | 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 | 2946 | # para obtener spwd |
|
|
2718 | 2947 | for ic in range(nPairs): |
|
|
2719 | 2948 | pair = crosspairs[ic] |
| @@ -2748,29 +2977,20 class SpectralFitting(Operation): | |||
|
|
2748 | 2977 | self.i=0 |
|
|
2749 | 2978 | self.isConfig = False |
|
|
2750 | 2979 | |
|
|
2751 | ||
|
|
2752 | 2980 | def setup(self,nChan,nProf,nHei,nBlocks): |
|
|
2753 | 2981 | self.__dataReady = False |
|
|
2754 | 2982 | self.bloques = numpy.zeros([2, nProf, nHei,nBlocks], dtype= complex) |
|
|
2755 | 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 | 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] | |
|
|
2763 | # self.bloques[:,:,:,k] = cspectra[:,:,0:nHei] | |
|
|
2764 | #if self.i==nBlocks: | |
|
|
2765 | # self.i==0 | |
|
|
2766 | rfunc = cspectra.copy() #self.bloques | |
|
|
2987 | rfunc = cspectra.copy() | |
|
|
2767 | 2988 | n_funct = len(rfunc[0,:,0,0]) |
|
|
2768 |
val_spc = spectra*0.0 |
|
|
|
2769 |
val_cspc = cspectra*0.0 |
|
|
|
2770 |
in_sat_spectra = spectra.copy() |
|
|
|
2771 |
in_sat_cspectra = cspectra.copy() |
|
|
|
2989 | val_spc = spectra*0.0 | |
|
|
2990 | val_cspc = cspectra*0.0 | |
|
|
2991 | in_sat_spectra = spectra.copy() | |
|
|
2992 | in_sat_cspectra = cspectra.copy() | |
|
|
2772 | 2993 | |
|
|
2773 | #print( rfunc.shape) | |
|
|
2774 | 2994 | min_hei = 200 |
|
|
2775 | 2995 | nProf = dataOut.nProfiles |
|
|
2776 | 2996 | heights = dataOut.heightList |
| @@ -2781,20 +3001,19 class SpectralFitting(Operation): | |||
|
|
2781 | 3001 | nPairs = len(crosspairs) |
|
|
2782 | 3002 | hval=(heights >= min_hei).nonzero() |
|
|
2783 | 3003 | ih=hval[0] |
|
|
2784 | #print numpy.absolute(rfunc[:,0,0,14]) | |
|
|
2785 | 3004 | for ih in range(hval[0][0],nHei): |
|
|
2786 | 3005 | for ifreq in range(nProf): |
|
|
2787 | 3006 | for ii in range(n_funct): |
|
|
2788 | 3007 | |
|
|
2789 | 3008 | func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) |
|
|
2790 | #print numpy.amin(func2clean) | |
|
|
3009 | ||
|
|
2791 | 3010 | val = (numpy.isfinite(func2clean)==True).nonzero() |
|
|
2792 | 3011 | if len(val)>0: |
|
|
2793 | 3012 | min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40) |
|
|
2794 | 3013 | if min_val <= -40 : min_val = -40 |
|
|
2795 | 3014 | max_val = numpy.around(numpy.amax(func2clean)+2) #< 200 |
|
|
2796 | 3015 | if max_val >= 200 : max_val = 200 |
|
|
2797 | #print min_val, max_val | |
|
|
3016 | ||
|
|
2798 | 3017 | step = 1 |
|
|
2799 | 3018 | #Getting bins and the histogram |
|
|
2800 | 3019 | x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step |
| @@ -2809,18 +3028,6 class SpectralFitting(Operation): | |||
|
|
2809 | 3028 | except: |
|
|
2810 | 3029 | mode = mean |
|
|
2811 | 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 | 3032 | #Removing echoes greater than mode + 3*stdv |
|
|
2826 | 3033 | factor_stdv = 2.5 |
| @@ -2831,29 +3038,17 class SpectralFitting(Operation): | |||
|
|
2831 | 3038 | cross_pairs = crosspairs[ii] |
|
|
2832 | 3039 | #Getting coherent echoes which are removed. |
|
|
2833 | 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 | 3041 | val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1 |
|
|
2837 | 3042 | val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1 |
|
|
2838 | 3043 | val_cspc[novall[0],ii,ifreq,ih] = 1 |
|
|
2839 | #print("OUT NOVALL 1") | |
|
|
2840 | 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 | 3045 | spectra[noval,cross_pairs[0],ifreq,ih] = numpy.nan |
|
|
2844 | 3046 | spectra[noval,cross_pairs[1],ifreq,ih] = numpy.nan |
|
|
2845 | 3047 | cspectra[noval,ii,ifreq,ih] = numpy.nan |
|
|
2846 | # if ih == 17 and ii == 0 and ifreq ==0 : | |
|
|
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] | |
|
|
3048 | # | |
|
|
2854 | 3049 | #no sale es para savedrifts >2 |
|
|
2855 | ''' channels = channels | |
|
|
2856 |
cross_pairs = |
|
|
|
3050 | ''' channels = dataOut.channelList | |
|
|
3051 | cross_pairs = dataOut.groupList | |
|
|
2857 | 3052 | #print("OUT NOVALL 2") |
|
|
2858 | 3053 | |
|
|
2859 | 3054 | vcross0 = (cross_pairs[0] == channels[ii]).nonzero() |
| @@ -2874,6 +3069,7 class SpectralFitting(Operation): | |||
|
|
2874 | 3069 | self.bloques[vcross,ifreq,ih,noval] = numpy.nan |
|
|
2875 | 3070 | ''' |
|
|
2876 | 3071 | #Getting average of the spectra and cross-spectra from incoherent echoes. |
|
|
3072 | ||
|
|
2877 | 3073 | out_spectra = numpy.zeros([nChan,nProf,nHei], dtype=float) #+numpy.nan |
|
|
2878 | 3074 | out_cspectra = numpy.zeros([nPairs,nProf,nHei], dtype=complex) #+numpy.nan |
|
|
2879 | 3075 | for ih in range(nHei): |
| @@ -2881,22 +3077,15 class SpectralFitting(Operation): | |||
|
|
2881 | 3077 | for ich in range(nChan): |
|
|
2882 | 3078 | tmp = spectra[:,ich,ifreq,ih] |
|
|
2883 | 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 | 3080 | if len(valid[0]) >0 : |
|
|
2890 | 3081 | out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0]) |
|
|
2891 | #for icr in range(nPairs): | |
|
|
3082 | ||
|
|
2892 | 3083 | for icr in range(nPairs): |
|
|
2893 | 3084 | tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih]) |
|
|
2894 | 3085 | valid = (numpy.isfinite(tmp)==True).nonzero() |
|
|
2895 | 3086 | if len(valid[0]) > 0: |
|
|
2896 | 3087 | out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0]) |
|
|
2897 | # print('##########################################################') | |
|
|
2898 | 3088 | #Removing fake coherent echoes (at least 4 points around the point) |
|
|
2899 | ||
|
|
2900 | 3089 | val_spectra = numpy.sum(val_spc,0) |
|
|
2901 | 3090 | val_cspectra = numpy.sum(val_cspc,0) |
|
|
2902 | 3091 | |
| @@ -2932,13 +3121,6 class SpectralFitting(Operation): | |||
|
|
2932 | 3121 | tmp_sat_cspectra = cspectra.copy() |
|
|
2933 | 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 | 3124 | val = (val_spc > 0).nonzero() |
|
|
2943 | 3125 | if len(val[0]) > 0: |
|
|
2944 | 3126 | tmp_sat_spectra[val] = in_sat_spectra[val] |
| @@ -2963,52 +3145,29 class SpectralFitting(Operation): | |||
|
|
2963 | 3145 | valid = (numpy.isfinite(tmp)).nonzero() |
|
|
2964 | 3146 | if len(valid[0]) > 0: |
|
|
2965 | 3147 | sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0]) |
|
|
2966 | #self.__dataReady= True | |
|
|
2967 | #sat_spectra, sat_cspectra= sat_spectra, sat_cspectra | |
|
|
2968 | #if not self.__dataReady: | |
|
|
2969 | #return None, None | |
|
|
3148 | ||
|
|
2970 | 3149 | return out_spectra, out_cspectra,sat_spectra,sat_cspectra |
|
|
2971 | 3150 | def REM_ISOLATED_POINTS(self,array,rth): |
|
|
2972 | # import matplotlib.pyplot as plt | |
|
|
2973 | 3151 | if rth == None : rth = 4 |
|
|
2974 | ||
|
|
2975 | 3152 | num_prof = len(array[0,:,0]) |
|
|
2976 | 3153 | num_hei = len(array[0,0,:]) |
|
|
2977 | 3154 | n2d = len(array[:,0,0]) |
|
|
2978 | 3155 | |
|
|
2979 | 3156 | for ii in range(n2d) : |
|
|
2980 | #print ii,n2d | |
|
|
2981 | 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 | 3158 | tmp = numpy.reshape(tmp,num_prof*num_hei) |
|
|
2995 | 3159 | indxs1 = (numpy.isfinite(tmp)==True).nonzero() |
|
|
2996 | 3160 | indxs2 = (tmp > 0).nonzero() |
|
|
2997 | ||
|
|
2998 | 3161 | indxs1 = (indxs1[0]) |
|
|
2999 | 3162 | indxs2 = indxs2[0] |
|
|
3000 | #indxs1 = numpy.array(indxs1[0]) | |
|
|
3001 | #indxs2 = numpy.array(indxs2[0]) | |
|
|
3002 | 3163 | indxs = None |
|
|
3003 | #print indxs1 , indxs2 | |
|
|
3164 | ||
|
|
3004 | 3165 | for iv in range(len(indxs2)): |
|
|
3005 | 3166 | indv = numpy.array((indxs1 == indxs2[iv]).nonzero()) |
|
|
3006 | #print len(indxs2), indv | |
|
|
3007 | 3167 | if len(indv[0]) > 0 : |
|
|
3008 | 3168 | indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None) |
|
|
3009 | # print indxs | |
|
|
3169 | ||
|
|
3010 | 3170 | indxs = indxs[1:] |
|
|
3011 | #print indxs, len(indxs) | |
|
|
3012 | 3171 | if len(indxs) < 4 : |
|
|
3013 | 3172 | array[ii,:,:] = 0. |
|
|
3014 | 3173 | return |
| @@ -3016,7 +3175,6 class SpectralFitting(Operation): | |||
|
|
3016 | 3175 | xpos = numpy.mod(indxs ,num_hei) |
|
|
3017 | 3176 | ypos = (indxs / num_hei) |
|
|
3018 | 3177 | sx = numpy.argsort(xpos) # Ordering respect to "x" (time) |
|
|
3019 | #print sx | |
|
|
3020 | 3178 | xpos = xpos[sx] |
|
|
3021 | 3179 | ypos = ypos[sx] |
|
|
3022 | 3180 | |
| @@ -3024,36 +3182,25 class SpectralFitting(Operation): | |||
|
|
3024 | 3182 | ic = 0 |
|
|
3025 | 3183 | while True : |
|
|
3026 | 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) | |
|
|
3028 | #plt.plot(r) | |
|
|
3029 | #plt.show() | |
|
|
3185 | ||
|
|
3030 | 3186 | no_coh1 = (numpy.isfinite(r)==True).nonzero() |
|
|
3031 | 3187 | no_coh2 = (r <= rth).nonzero() |
|
|
3032 | #print r, no_coh1, no_coh2 | |
|
|
3033 | 3188 | no_coh1 = numpy.array(no_coh1[0]) |
|
|
3034 | 3189 | no_coh2 = numpy.array(no_coh2[0]) |
|
|
3035 | 3190 | no_coh = None |
|
|
3036 | #print valid1 , valid2 | |
|
|
3037 | 3191 | for iv in range(len(no_coh2)): |
|
|
3038 | 3192 | indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero()) |
|
|
3039 | 3193 | if len(indv[0]) > 0 : |
|
|
3040 | 3194 | no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None) |
|
|
3041 | 3195 | no_coh = no_coh[1:] |
|
|
3042 | #print len(no_coh), no_coh | |
|
|
3043 | 3196 | if len(no_coh) < 4 : |
|
|
3044 | #print xpos[ic], ypos[ic], ic | |
|
|
3045 | # plt.plot(r) | |
|
|
3046 | # plt.show() | |
|
|
3047 | 3197 | xpos[ic] = numpy.nan |
|
|
3048 | 3198 | ypos[ic] = numpy.nan |
|
|
3049 | 3199 | |
|
|
3050 | 3200 | ic = ic + 1 |
|
|
3051 | 3201 | if (ic == len(indxs)) : |
|
|
3052 | 3202 | break |
|
|
3053 | #print( xpos, ypos) | |
|
|
3054 | ||
|
|
3055 | 3203 | indxs = (numpy.isfinite(list(xpos))==True).nonzero() |
|
|
3056 | #print indxs[0] | |
|
|
3057 | 3204 | if len(indxs[0]) < 4 : |
|
|
3058 | 3205 | array[ii,:,:] = 0. |
|
|
3059 | 3206 | return |
| @@ -3068,27 +3215,11 class SpectralFitting(Operation): | |||
|
|
3068 | 3215 | tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))] |
|
|
3069 | 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 | 3218 | return array |
|
|
3084 | 3219 | def moments(self,doppler,yarray,npoints): |
|
|
3085 | 3220 | ytemp = yarray |
|
|
3086 | #val = WHERE(ytemp GT 0,cval) | |
|
|
3087 | #if cval == 0 : val = range(npoints-1) | |
|
|
3088 | 3221 | val = (ytemp > 0).nonzero() |
|
|
3089 | 3222 | val = val[0] |
|
|
3090 | #print('hvalid:',hvalid) | |
|
|
3091 | #print('valid', valid) | |
|
|
3092 | 3223 | if len(val) == 0 : val = range(npoints-1) |
|
|
3093 | 3224 | |
|
|
3094 | 3225 | ynew = 0.5*(ytemp[val[0]]+ytemp[val[len(val)-1]]) |
| @@ -3102,19 +3233,289 class SpectralFitting(Operation): | |||
|
|
3102 | 3233 | fmom = numpy.sum(doppler*ytemp)/numpy.sum(ytemp)+(index-(npoints/2-1))*numpy.abs(doppler[1]-doppler[0]) |
|
|
3103 | 3234 | smom = numpy.sum(doppler*doppler*ytemp)/numpy.sum(ytemp) |
|
|
3104 | 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 | 3505 | index = 0 |
|
|
3107 | 3506 | fint = 0 |
|
|
3108 | 3507 | buffer = 0 |
|
|
3109 | 3508 | buffer2 = 0 |
|
|
3110 | 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 | 3511 | nChannels = dataOut.nChannels |
|
|
3113 | 3512 | nHeights= dataOut.heightList.size |
|
|
3114 | 3513 | nProf = dataOut.nProfiles |
|
|
3514 | if numpy.any(taver): taver=int(taver) | |
|
|
3515 | else : taver = 15 | |
|
|
3115 | 3516 | tini=time.localtime(dataOut.utctime) |
|
|
3116 |
if (tini.tm_min % |
|
|
|
3117 | # print tini.tm_min | |
|
|
3517 | if (tini.tm_min % taver) == 0 and (tini.tm_sec < 5 and self.fint==0): | |
|
|
3518 | ||
|
|
3118 | 3519 | self.index = 0 |
|
|
3119 | 3520 | jspc = self.buffer |
|
|
3120 | 3521 | jcspc = self.buffer2 |
| @@ -3124,14 +3525,14 class SpectralFitting(Operation): | |||
|
|
3124 | 3525 | self.buffer3 = dataOut.noise |
|
|
3125 | 3526 | self.fint = 1 |
|
|
3126 | 3527 | if numpy.any(jspc) : |
|
|
3127 |
jspc= numpy.reshape(jspc,(int(len(jspc)/ |
|
|
|
3128 | jcspc= numpy.reshape(jcspc,(int(len(jcspc)/2),2,nProf,nHeights)) | |
|
|
3129 |
jnoise= numpy.reshape(jnoise,(int(len(jnoise)/ |
|
|
|
3528 | jspc= numpy.reshape(jspc,(int(len(jspc)/nChannels),nChannels,nProf,nHeights)) | |
|
|
3529 | jcspc= numpy.reshape(jcspc,(int(len(jcspc)/int(nChannels/2)),int(nChannels/2),nProf,nHeights)) | |
|
|
3530 | jnoise= numpy.reshape(jnoise,(int(len(jnoise)/nChannels),nChannels)) | |
|
|
3130 | 3531 | else: |
|
|
3131 | 3532 | dataOut.flagNoData = True |
|
|
3132 | 3533 | return dataOut |
|
|
3133 | 3534 | else : |
|
|
3134 |
if (tini.tm_min % |
|
|
|
3535 | if (tini.tm_min % taver) == 0 : self.fint = 1 | |
|
|
3135 | 3536 | else : self.fint = 0 |
|
|
3136 | 3537 | self.index += 1 |
|
|
3137 | 3538 | if numpy.any(self.buffer): |
| @@ -3146,7 +3547,13 class SpectralFitting(Operation): | |||
|
|
3146 | 3547 | return dataOut |
|
|
3147 | 3548 | if path != None: |
|
|
3148 | 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 | 3558 | #To be inserted as a parameter |
|
|
3152 | 3559 | groupArray = numpy.array(groupList) |
| @@ -3162,8 +3569,11 class SpectralFitting(Operation): | |||
|
|
3162 | 3569 | dataOut.data_paramC = None |
|
|
3163 | 3570 | |
|
|
3164 | 3571 | #Set constants |
|
|
3165 | constants = self.library.setConstants(dataOut) | |
|
|
3166 | dataOut.constants = constants | |
|
|
3572 | try: | |
|
|
3573 | constants = self.library.setConstants(dataOut) | |
|
|
3574 | dataOut.constants = constants | |
|
|
3575 | except: | |
|
|
3576 | pass | |
|
|
3167 | 3577 | M = dataOut.normFactor |
|
|
3168 | 3578 | N = dataOut.nFFTPoints |
|
|
3169 | 3579 | ippSeconds = dataOut.ippSeconds |
| @@ -3190,210 +3600,300 class SpectralFitting(Operation): | |||
|
|
3190 | 3600 | if not self.isConfig: |
|
|
3191 | 3601 | self.isConfig = True |
|
|
3192 | 3602 | |
|
|
3193 |
index = tini.tm_hour*12+tini.tm_min/ |
|
|
|
3603 | index = tini.tm_hour*12+tini.tm_min/taver | |
|
|
3194 | 3604 | jspc = jspc/N/N |
|
|
3195 | 3605 | jcspc = jcspc/N/N |
|
|
3196 | 3606 | tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.CleanRayleigh(dataOut,jspc,jcspc,2) |
|
|
3197 | 3607 | jspectra = tmp_spectra*len(jspc[:,0,0,0]) |
|
|
3198 | 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, |
|
|
|
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 | 3610 | clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index) |
|
|
3201 | 3611 | dataOut.data_spc = incoh_spectra |
|
|
3202 | 3612 | dataOut.data_cspc = incoh_cspectra |
|
|
3613 | #dataOut.data_spc = tmp_spectra | |
|
|
3614 | #dataOut.data_cspc = tmp_cspectra | |
|
|
3203 | 3615 | |
|
|
3204 | 3616 | clean_num_aver = incoh_aver*len(jspc[:,0,0,0]) |
|
|
3205 | 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 | 3620 | #List of possible combinations |
|
|
3207 | 3621 | listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2) |
|
|
3208 | 3622 | indCross = numpy.zeros(len(list(listComb)), dtype = 'int') |
|
|
3209 | ||
|
|
3210 | if getSNR: | |
|
|
3211 | listChannels = groupArray.reshape((groupArray.size)) | |
|
|
3212 | listChannels.sort() | |
|
|
3213 | dataOut.data_SNR = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise[listChannels]) | |
|
|
3214 | if dataOut.data_paramC is None: | |
|
|
3215 | dataOut.data_paramC = numpy.zeros((nGroups*4, nHeights,2))*numpy.nan | |
|
|
3216 | for i in range(nGroups): | |
|
|
3217 | coord = groupArray[i,:] | |
|
|
3218 | #Input data array | |
|
|
3219 | data = dataOut.data_spc[coord,:,:]/(M*N) | |
|
|
3220 | data = data.reshape((data.shape[0]*data.shape[1],data.shape[2])) | |
|
|
3221 | ||
|
|
3222 | #Cross Spectra data array for Covariance Matrixes | |
|
|
3223 | ind = 0 | |
|
|
3224 | for pairs in listComb: | |
|
|
3225 | pairsSel = numpy.array([coord[x],coord[y]]) | |
|
|
3226 | indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0]) | |
|
|
3227 |
|
|
|
|
3228 | dataCross = dataOut.data_cspc[indCross,:,:]/(M*N) | |
|
|
3229 | dataCross = dataCross**2 | |
|
|
3230 | nhei = nHeights | |
|
|
3231 | poweri = numpy.sum(dataOut.data_spc[:,1:nProf-0,:],axis=1)/clean_num_aver[:,:] | |
|
|
3232 | if i == 0 : my_noises = numpy.zeros(4,dtype=float) #FLTARR(4) | |
|
|
3233 | n0i = numpy.nanmin(poweri[0+i*2,0:nhei-0])/(nProf-1) | |
|
|
3234 | n1i = numpy.nanmin(poweri[1+i*2,0:nhei-0])/(nProf-1) | |
|
|
3235 | n0 = n0i | |
|
|
3236 | n1= n1i | |
|
|
3237 | my_noises[2*i+0] = n0 | |
|
|
3238 | my_noises[2*i+1] = n1 | |
|
|
3239 | snrth = -16.0 | |
|
|
3240 | snrth = 10**(snrth/10.0) | |
|
|
3241 | ||
|
|
3242 | for h in range(nHeights): | |
|
|
3243 | d = data[:,h] | |
|
|
3244 | smooth = clean_num_aver[i+1,h] #dataOut.data_spc[:,1:nProf-0,:] | |
|
|
3245 | signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth | |
|
|
3246 | signalpn1 = (dataOut.data_spc[i*2+1,1:(nProf-0),h])/smooth | |
|
|
3247 | signal0 = signalpn0-n0 | |
|
|
3248 | signal1 = signalpn1-n1 | |
|
|
3249 | snr0 = numpy.sum(signal0/n0)/(nProf-1) | |
|
|
3250 | snr1 = numpy.sum(signal1/n1)/(nProf-1) | |
|
|
3251 | if snr0 > snrth and snr1 > snrth and clean_num_aver[i+1,h] > 0 : | |
|
|
3252 | #Covariance Matrix | |
|
|
3253 |
|
|
|
|
3254 | ind = 0 | |
|
|
3255 | for pairs in listComb: | |
|
|
3256 | #Coordinates in Covariance Matrix | |
|
|
3257 |
|
|
|
|
3258 |
|
|
|
|
3259 |
# |
|
|
|
3260 | S12 = dataCross[ind,:,h] | |
|
|
3261 | D12 = numpy.diag(S12) | |
|
|
3262 | #Completing Covariance Matrix with Cross Spectras | |
|
|
3263 | D[x*N:(x+1)*N,y*N:(y+1)*N] = D12 | |
|
|
3264 | D[y*N:(y+1)*N,x*N:(x+1)*N] = D12 | |
|
|
3265 | ind += 1 | |
|
|
3266 | diagD = numpy.zeros(256) | |
|
|
3267 | if h == 17 : | |
|
|
3268 | for ii in range(256): diagD[ii] = D[ii,ii] | |
|
|
3269 | #Dinv=numpy.linalg.inv(D) | |
|
|
3270 | #L=numpy.linalg.cholesky(Dinv) | |
|
|
3271 |
|
|
|
|
3272 | Dinv=numpy.linalg.inv(D) | |
|
|
3273 | L=numpy.linalg.cholesky(Dinv) | |
|
|
3274 | except: | |
|
|
3275 | Dinv = D*numpy.nan | |
|
|
3276 | L= D*numpy.nan | |
|
|
3277 |
|
|
|
|
3278 | ||
|
|
3279 | dp = numpy.dot(LT,d) | |
|
|
3280 | ||
|
|
3281 | #Initial values | |
|
|
3282 | data_spc = dataOut.data_spc[coord,:,h] | |
|
|
3283 | ||
|
|
3284 | if (h>0)and(error1[3]<5): | |
|
|
3285 | p0 = dataOut.data_param[i,:,h-1] | |
|
|
3286 |
|
|
|
|
3287 | p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants))# sin el i(data_spc, constants, i) | |
|
|
3288 |
|
|
|
|
3289 |
|
|
|
|
3290 | #print (dp,LT,constants) | |
|
|
3291 | #value =self.__residFunction(p0,dp,LT,constants) | |
|
|
3292 | #print ("valueREADY",value.shape, type(value)) | |
|
|
3293 | #optimize.leastsq(value) | |
|
|
3294 | minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True) | |
|
|
3295 | #minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants)) | |
|
|
3296 | #Chi square error | |
|
|
3297 | #print(minp,covp.infodict,mesg,ier) | |
|
|
3298 | #print("REALIZA OPTIMIZ") | |
|
|
3299 | error0 = numpy.sum(infodict['fvec']**2)/(2*N) | |
|
|
3300 | #Error with Jacobian | |
|
|
3301 | error1 = self.library.errorFunction(minp,constants,LT) | |
|
|
3302 | # print self.__residFunction(p0,dp,LT, constants) | |
|
|
3303 | # print infodict['fvec'] | |
|
|
3304 | # print self.__residFunction(minp,dp,LT,constants) | |
|
|
3305 | ||
|
|
3306 | except: | |
|
|
3623 | if Gaussian_Windowed == 1: | |
|
|
3624 | #dataOut.data_spc = jspectra | |
|
|
3625 | ''' | |
|
|
3626 | Written by R. Flores | |
|
|
3627 | ''' | |
|
|
3628 | print("normFactor: ", dataOut.normFactor) | |
|
|
3629 | data_spc_aux = numpy.copy(dataOut.data_spc)#[:,0,:] | |
|
|
3630 | data_spc_aux[:,0,:] = (data_spc_aux[:,1,:]+data_spc_aux[:,-1,:])/2 | |
|
|
3631 | #''' | |
|
|
3632 | from scipy.signal import medfilt | |
|
|
3633 | import matplotlib.pyplot as plt | |
|
|
3634 | dataOut.moments = numpy.ones((dataOut.nChannels,4,dataOut.nHeights))*numpy.NAN | |
|
|
3635 | dataOut.VelRange = dataOut.getVelRange(0) | |
|
|
3636 | for nChannel in range(dataOut.nChannels): | |
|
|
3637 | for hei in range(dataOut.heightList.shape[0]): | |
|
|
3638 | #print("ipp: ", dataOut.ippSeconds) | |
|
|
3639 | #spc = numpy.copy(dataOut.data_spc[nChannel,:,hei]) | |
|
|
3640 | spc = data_spc_aux[nChannel,:,hei] | |
|
|
3641 | if spc.all() == 0.: | |
|
|
3642 | print("CONTINUE") | |
|
|
3643 | continue | |
|
|
3644 | #print(VelRange) | |
|
|
3645 | #print(dataOut.getFreqRange(64)) | |
|
|
3646 | #print("Hei: ", dataOut.heightList[hei]) | |
|
|
3647 | ||
|
|
3648 | spc_mod = numpy.copy(spc) | |
|
|
3649 | spcm = medfilt(spc_mod,11) | |
|
|
3650 | spc_max = numpy.max(spcm) | |
|
|
3651 | dop1_x0 = dataOut.VelRange[numpy.argmax(spcm)] | |
|
|
3652 | #D = numpy.min(spcm) | |
|
|
3653 | D_in = (numpy.mean(spcm[:15])+numpy.mean(spcm[-15:]))/2. | |
|
|
3654 | #print("spc_max: ", spc_max) | |
|
|
3655 | #print("dataOut.ippSeconds: ", dataOut.ippSeconds, dataOut.timeInterval) | |
|
|
3656 | ##fun, A, B, C, D = self.windowing_single(spc,dataOut.VelRange,spc_max,dop1_x0,abs(dop1_x0),D,dataOut.nFFTPoints) | |
|
|
3657 | #fun, A, B, C, D = self.windowing_single(spc,dataOut.VelRange,spc_max,dop1_x0,abs(dop1_x0),D,dataOut.nFFTPoints) | |
|
|
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) | |
|
|
3659 | ||
|
|
3660 | dataOut.moments[nChannel,0,hei] = A | |
|
|
3661 | dataOut.moments[nChannel,1,hei] = B | |
|
|
3662 | dataOut.moments[nChannel,2,hei] = C | |
|
|
3663 | dataOut.moments[nChannel,3,hei] = D | |
|
|
3664 | ''' | |
|
|
3665 | if nChannel == 0: | |
|
|
3666 | print(dataOut.heightList[hei]) | |
|
|
3667 | plt.figure() | |
|
|
3668 | plt.plot(dataOut.VelRange,spc,marker='*',linestyle='--') | |
|
|
3669 | plt.plot(dataOut.VelRange,fun) | |
|
|
3670 | plt.title(dataOut.heightList[hei]) | |
|
|
3671 | plt.show() | |
|
|
3672 | ''' | |
|
|
3673 | #plt.show() | |
|
|
3674 | #''' | |
|
|
3675 | dataOut.data_spc = jspectra | |
|
|
3676 | print("SUCCESS") | |
|
|
3677 | return dataOut | |
|
|
3678 | ||
|
|
3679 | elif Gaussian_Windowed == 2: #Only to clean spc | |
|
|
3680 | dataOut.VelRange = dataOut.getVelRange(0) | |
|
|
3681 | return dataOut | |
|
|
3682 | else: | |
|
|
3683 | if getSNR: | |
|
|
3684 | listChannels = groupArray.reshape((groupArray.size)) | |
|
|
3685 | listChannels.sort() | |
|
|
3686 | dataOut.data_SNR = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise[listChannels]) | |
|
|
3687 | if dataOut.data_paramC is None: | |
|
|
3688 | dataOut.data_paramC = numpy.zeros((nGroups*4, nHeights,2))*numpy.nan | |
|
|
3689 | for i in range(nGroups): | |
|
|
3690 | coord = groupArray[i,:] | |
|
|
3691 | #Input data array | |
|
|
3692 | data = dataOut.data_spc[coord,:,:]/(M*N) | |
|
|
3693 | data = data.reshape((data.shape[0]*data.shape[1],data.shape[2])) | |
|
|
3694 | ||
|
|
3695 | #Cross Spectra data array for Covariance Matrixes | |
|
|
3696 | ind = 0 | |
|
|
3697 | for pairs in listComb: | |
|
|
3698 | pairsSel = numpy.array([coord[x],coord[y]]) | |
|
|
3699 | indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0]) | |
|
|
3700 | ind += 1 | |
|
|
3701 | dataCross = dataOut.data_cspc[indCross,:,:]/(M*N) | |
|
|
3702 | dataCross = dataCross**2 | |
|
|
3703 | nhei = nHeights | |
|
|
3704 | poweri = numpy.sum(dataOut.data_spc[:,1:nProf-0,:],axis=1)/clean_num_aver[:,:] | |
|
|
3705 | if i == 0 : my_noises = numpy.zeros(4,dtype=float) | |
|
|
3706 | n0i = numpy.nanmin(poweri[0+i*2,0:nhei-0])/(nProf-1) | |
|
|
3707 | n1i = numpy.nanmin(poweri[1+i*2,0:nhei-0])/(nProf-1) | |
|
|
3708 | n0 = n0i | |
|
|
3709 | n1= n1i | |
|
|
3710 | my_noises[2*i+0] = n0 | |
|
|
3711 | my_noises[2*i+1] = n1 | |
|
|
3712 | snrth = -25.0 # -4 | |
|
|
3713 | snrth = 10**(snrth/10.0) | |
|
|
3714 | jvelr = numpy.zeros(nHeights, dtype = 'float') | |
|
|
3715 | #snr0 = numpy.zeros(nHeights, dtype = 'float') | |
|
|
3716 | #snr1 = numpy.zeros(nHeights, dtype = 'float') | |
|
|
3717 | hvalid = [0] | |
|
|
3718 | ||
|
|
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,:]) | |
|
|
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 | 3837 | minp = p0*numpy.nan |
|
|
3308 | 3838 | error0 = numpy.nan |
|
|
3309 | 3839 | error1 = p0*numpy.nan |
|
|
3310 | #print ("EXCEPT 0000000000") | |
|
|
3311 | # s_sq = (self.__residFunction(minp,dp,LT,constants)).sum()/(len(dp)-len(p0)) | |
|
|
3312 | # covp = covp*s_sq | |
|
|
3313 | # #print("TRY___________________________________________1") | |
|
|
3314 | # error = [] | |
|
|
3315 | # for ip in range(len(minp)): | |
|
|
3316 | # try: | |
|
|
3317 | # error.append(numpy.absolute(covp[ip][ip])**0.5) | |
|
|
3318 | # except: | |
|
|
3319 | # error.append( 0.00 ) | |
|
|
3320 | else : | |
|
|
3321 | data_spc = dataOut.data_spc[coord,:,h] | |
|
|
3322 | p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants)) | |
|
|
3323 | minp = p0*numpy.nan | |
|
|
3324 |
|
|
|
|
3325 |
|
|
|
|
3326 | #Save | |
|
|
3327 | if dataOut.data_param is None: | |
|
|
3328 | dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan | |
|
|
3329 | dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan | |
|
|
3330 | ||
|
|
3331 | dataOut.data_error[i,:,h] = numpy.hstack((error0,error1)) | |
|
|
3332 | dataOut.data_param[i,:,h] = minp | |
|
|
3333 | ||
|
|
3334 | for ht in range(nHeights-1) : | |
|
|
3335 | smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam] | |
|
|
3336 |
dataOut.data_paramC[4*i,ht, |
|
|
|
3337 | signalpn0 = (coh_spectra[i*2 ,1:(nProf-0),ht])/smooth #coh_spectra | |
|
|
3338 | signalpn1 = (coh_spectra[i*2+1,1:(nProf-0),ht])/smooth | |
|
|
3339 | ||
|
|
3340 |
|
|
|
|
3341 |
|
|
|
|
3342 |
|
|
|
|
3343 | #print('hvalid:',hvalid) | |
|
|
3344 | #print('valid', valid) | |
|
|
3345 |
if len(val |
|
|
|
3346 | else : val0_npoints = len(val0) | |
|
|
3347 | ||
|
|
3348 | #val1 = WHERE(signalpn1 > 0,cval1) | |
|
|
3349 | val1 = (signalpn1 > 0).nonzero() | |
|
|
3350 |
|
|
|
|
3351 | if len(val1) == 0 : val1_npoints = nProf | |
|
|
3352 | else : val1_npoints = len(val1) | |
|
|
3353 | ||
|
|
3354 | dataOut.data_paramC[0+4*i,ht,0] = numpy.sum((signalpn0/val0_npoints))/n0 | |
|
|
3355 | dataOut.data_paramC[1+4*i,ht,0] = numpy.sum((signalpn1/val1_npoints))/n1 | |
|
|
3356 | ||
|
|
3357 | signal0 = (signalpn0-n0) # > 0 | |
|
|
3358 | vali = (signal0 < 0).nonzero() | |
|
|
3359 | vali = vali[0] | |
|
|
3360 | if len(vali) > 0 : signal0[vali] = 0 | |
|
|
3361 | signal1 = (signalpn1-n1) #> 0 | |
|
|
3362 | vali = (signal1 < 0).nonzero() | |
|
|
3363 |
|
|
|
|
3364 | if len(vali) > 0 : signal1[vali] = 0 | |
|
|
3365 | snr0 = numpy.sum(signal0/n0)/(nProf-1) | |
|
|
3366 | snr1 = numpy.sum(signal1/n1)/(nProf-1) | |
|
|
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 | |
|
|
3840 | ||
|
|
3841 | if dataOut.data_param is None: | |
|
|
3842 | dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan | |
|
|
3843 | dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan | |
|
|
3844 | ||
|
|
3845 | dataOut.data_error[i,:,h] = numpy.hstack((error0,error1)) | |
|
|
3846 | dataOut.data_param[i,:,h] = minp | |
|
|
3847 | ||
|
|
3848 | for ht in range(nHeights-1) : | |
|
|
3849 | smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam] | |
|
|
3850 | dataOut.data_paramC[4*i,ht,1] = smooth | |
|
|
3851 | signalpn0 = (coh_spectra[i*2 ,1:(nProf-0),ht])/smooth #coh_spectra | |
|
|
3852 | signalpn1 = (coh_spectra[i*2+1,1:(nProf-0),ht])/smooth | |
|
|
3853 | ||
|
|
3854 | val0 = (signalpn0 > 0).nonzero() | |
|
|
3855 | val0 = val0[0] | |
|
|
3856 | ||
|
|
3857 | if len(val0) == 0 : val0_npoints = nProf | |
|
|
3858 | else : val0_npoints = len(val0) | |
|
|
3859 | ||
|
|
3860 | val1 = (signalpn1 > 0).nonzero() | |
|
|
3861 | val1 = val1[0] | |
|
|
3862 | if len(val1) == 0 : val1_npoints = nProf | |
|
|
3863 | else : val1_npoints = len(val1) | |
|
|
3864 | ||
|
|
3865 | dataOut.data_paramC[0+4*i,ht,0] = numpy.sum((signalpn0/val0_npoints))/n0 | |
|
|
3866 | dataOut.data_paramC[1+4*i,ht,0] = numpy.sum((signalpn1/val1_npoints))/n1 | |
|
|
3867 | ||
|
|
3868 | signal0 = (signalpn0-n0) | |
|
|
3869 | vali = (signal0 < 0).nonzero() | |
|
|
3870 | vali = vali[0] | |
|
|
3871 | if len(vali) > 0 : signal0[vali] = 0 | |
|
|
3872 | signal1 = (signalpn1-n1) | |
|
|
3873 | vali = (signal1 < 0).nonzero() | |
|
|
3874 | vali = vali[0] | |
|
|
3875 | if len(vali) > 0 : signal1[vali] = 0 | |
|
|
3876 | snr0 = numpy.sum(signal0/n0)/(nProf-1) | |
|
|
3877 | snr1 = numpy.sum(signal1/n1)/(nProf-1) | |
|
|
3878 | doppler = absc[1:] | |
|
|
3879 | if snr0 >= snrth and snr1 >= snrth and smooth : | |
|
|
3880 | signalpn0_n0 = signalpn0 | |
|
|
3881 | signalpn0_n0[val0] = signalpn0[val0] - n0 | |
|
|
3882 | mom0 = self.moments(doppler,signalpn0-n0,nProf) | |
|
|
3883 | ||
|
|
3884 | signalpn1_n1 = signalpn1 | |
|
|
3885 | signalpn1_n1[val1] = signalpn1[val1] - n1 | |
|
|
3886 | mom1 = self.moments(doppler,signalpn1_n1,nProf) | |
|
|
3887 | dataOut.data_paramC[2+4*i,ht,0] = (mom0[0]+mom1[0])/2. | |
|
|
3888 | dataOut.data_paramC[3+4*i,ht,0] = (mom0[1]+mom1[1])/2. | |
|
|
3889 | ||
|
|
3890 | dataOut.data_spc = jspectra | |
|
|
3891 | if getSNR: | |
|
|
3892 | listChannels = groupArray.reshape((groupArray.size)) | |
|
|
3893 | listChannels.sort() | |
|
|
3894 | ||
|
|
3895 | dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:], my_noises[listChannels]) | |
|
|
3896 | return dataOut | |
|
|
3397 | 3897 | |
|
|
3398 | 3898 | def __residFunction(self, p, dp, LT, constants): |
|
|
3399 | 3899 | |
| @@ -6628,9 +7128,9 class IGRFModel(Operation): | |||
|
|
6628 | 7128 | mkfact_short_2020_2.mkfact(dataOut.year,dataOut.h,dataOut.bfm,dataOut.thb,dataOut.bki,dataOut.MAXNRANGENDT) |
|
|
6629 | 7129 | |
|
|
6630 | 7130 | #mkfact_short_2020.mkfact(dataOut.year,dataOut.h,dataOut.bfm,dataOut.thb,dataOut.bki,dataOut.MAXNRANGENDT) |
|
|
6631 | print("bki: ", dataOut.bki[:10]) | |
|
|
6632 | print("thb: ", dataOut.thb[:10]) | |
|
|
6633 | print("bfm: ", dataOut.bfm[:10]) | |
|
|
7131 | #print("bki: ", dataOut.bki[:10]) | |
|
|
7132 | #print("thb: ", dataOut.thb[:10]) | |
|
|
7133 | #print("bfm: ", dataOut.bfm[:10]) | |
|
|
6634 | 7134 | #print("IDs: ", id(dataOut.bki)) |
|
|
6635 | 7135 | |
|
|
6636 | 7136 | return dataOut |
| @@ -6794,8 +7294,10 class MergeProc(ProcessingUnit): | |||
|
|
6794 | 7294 | self.dataOut.NSHTS = int(numpy.shape(self.dataOut.ph2)[0]) |
|
|
6795 | 7295 | dH = self.dataOut.heightList[1]-self.dataOut.heightList[0] |
|
|
6796 | 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 | 7299 | self.dataOut.NDP = self.dataOut.NSHTS |
|
|
7300 | #exit(1) | |
|
|
6799 | 7301 | #print(self.dataOut.heightList) |
|
|
6800 | 7302 | |
|
|
6801 | 7303 | class MST_Den_Conv(Operation): |
| @@ -455,19 +455,43 class GetSNR(Operation): | |||
|
|
455 | 455 | |
|
|
456 | 456 | def run(self,dataOut): |
|
|
457 | 457 | |
|
|
458 | noise = dataOut.getNoise() | |
|
|
459 | maxdB = 16 | |
|
|
460 | ||
|
|
461 | normFactor = 24 | |
|
|
462 | ||
|
|
458 | #noise = dataOut.getNoise() | |
|
|
459 | noise = dataOut.getNoise(ymin_index=-10) #Región superior donde solo debería de haber ruido | |
|
|
460 | #print("Noise: ", noise) | |
|
|
461 | #print("Noise_dB: ", 10*numpy.log10(noise/dataOut.normFactor)) | |
|
|
462 | #print("Heights: ", dataOut.heightList) | |
|
|
463 | 463 | #dataOut.data_snr = (dataOut.data_spc.sum(axis=1))/(noise[:,None]*dataOut.normFactor) |
|
|
464 | dataOut.data_snr = (dataOut.data_spc.sum(axis=1))/(noise[:,None]*dataOut.nFFTPoints) | |
|
|
465 | ||
|
|
466 | dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.5, numpy.nan, dataOut.data_snr) | |
|
|
464 | ################dataOut.data_snr = (dataOut.data_spc.sum(axis=1))/(noise[:,None]*dataOut.nFFTPoints) #Before 12Jan2023 | |
|
|
465 | #dataOut.data_snr = (dataOut.data_spc.sum(axis=1)-noise[:,None])/(noise[:,None]) | |
|
|
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 | 490 | #dataOut.data_snr = 10*numpy.log10(dataOut.data_snr) |
|
|
468 | 491 | #dataOut.data_snr = numpy.expand_dims(dataOut.data_snr,axis=0) |
|
|
469 | 492 | #print(dataOut.data_snr.shape) |
|
|
470 | 493 | #exit(1) |
|
|
494 | #print("Before: ", dataOut.data_snr[0]) | |
|
|
471 | 495 | |
|
|
472 | 496 | |
|
|
473 | 497 | return dataOut |
| @@ -3815,6 +3815,9 class IncohInt(Operation): | |||
|
|
3815 | 3815 | dataOut.flagNoData = False |
|
|
3816 | 3816 | |
|
|
3817 | 3817 | dataOut.VelRange = dataOut.getVelRange(0) |
|
|
3818 | dataOut.FreqRange = dataOut.getFreqRange(0)/1000. | |
|
|
3819 | #print("VelRange: ", dataOut.VelRange) | |
|
|
3820 | #exit(1) | |
|
|
3818 | 3821 | |
|
|
3819 | 3822 | #print("Power",numpy.sum(dataOut.data_spc[0,:,20:30],axis=0)) |
|
|
3820 | 3823 | #print("Power",numpy.sum(dataOut.data_spc[0,100:110,:],axis=1)) |
| @@ -3921,7 +3924,7 class SnrFaraday(Operation): | |||
|
|
3921 | 3924 | |
|
|
3922 | 3925 | return dataOut |
|
|
3923 | 3926 | |
|
|
3924 |
class SpectraDataToFaraday |
|
|
|
3927 | class SpectraDataToFaraday(Operation): #ISR MODE | |
|
|
3925 | 3928 | ''' |
|
|
3926 | 3929 | Written by R. Flores |
|
|
3927 | 3930 | ''' |
| @@ -4144,10 +4147,14 class SpectraDataToFaraday_07_11_2022(Operation): | |||
|
|
4144 | 4147 | #print(data_eej) |
|
|
4145 | 4148 | index_eej = CleanCohEchoes.mad_based_outlier(self,data_eej[:17]) |
|
|
4146 | 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 | 4158 | #exit(1) |
|
|
4152 | 4159 | |
|
|
4153 | 4160 | def run(self,dataOut): |
| @@ -4204,7 +4211,7 class SpectraDataToFaraday_07_11_2022(Operation): | |||
|
|
4204 | 4211 | #exit(1) |
|
|
4205 | 4212 | return dataOut |
|
|
4206 | 4213 | |
|
|
4207 | class SpectraDataToFaraday(Operation): | |
|
|
4214 | class SpectraDataToFaraday_MST(Operation): #MST MODE | |
|
|
4208 | 4215 | """Operation to use spectra data in Faraday processing. |
|
|
4209 | 4216 | |
|
|
4210 | 4217 | Parameters: |
| 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 | 38 | if self.dataIn.type == 'Voltage': |
|
|
39 | 39 | self.dataOut.copy(self.dataIn) |
|
|
40 | 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 | 44 | def __updateObjFromAmisrInput(self): |
|
|
44 | 45 | |
| @@ -226,6 +227,7 class selectHeights(Operation): | |||
|
|
226 | 227 | #print(dataOut.nHeights) |
|
|
227 | 228 | #exit(1) |
|
|
228 | 229 | ''' |
|
|
230 | #print("select Heights: ", self.dataOut.heightList) | |
|
|
229 | 231 | return self.dataOut |
|
|
230 | 232 | |
|
|
231 | 233 | def selectHeightsByIndex(self, minIndex, maxIndex): |
| @@ -1250,7 +1252,7 class LagsReshapeLP(Operation): | |||
|
|
1250 | 1252 | self.buffer[3,:,:,:] = numpy.transpose(numpy.conj(self.lagp3),(1,2,0)) |
|
|
1251 | 1253 | |
|
|
1252 | 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 | 1256 | #print(self.cnorm) |
|
|
1255 | 1257 | exit(1) |
|
|
1256 | 1258 | |
| @@ -3348,7 +3350,7 class DoublePulseACFs_PerLag(Operation): | |||
|
|
3348 | 3350 | exit(1) |
|
|
3349 | 3351 | ''' |
|
|
3350 | 3352 | #print(pa) |
|
|
3351 |
|
|
|
|
3353 | ''' | |
|
|
3352 | 3354 | import matplotlib.pyplot as plt |
|
|
3353 | 3355 | #plt.plot(dataOut.p[:,-1],dataOut.heightList) |
|
|
3354 | 3356 | plt.plot(pa/dataOut.pan-1.,dataOut.heightList) |
| @@ -3358,7 +3360,7 class DoublePulseACFs_PerLag(Operation): | |||
|
|
3358 | 3360 | plt.show() |
|
|
3359 | 3361 | #print("p: ",dataOut.p[33,:]) |
|
|
3360 | 3362 | #exit(1) |
|
|
3361 |
|
|
|
|
3363 | ''' | |
|
|
3362 | 3364 | return dataOut |
|
|
3363 | 3365 | |
|
|
3364 | 3366 | class FaradayAngleAndDPPower(Operation): |
| @@ -3399,7 +3401,7 class FaradayAngleAndDPPower(Operation): | |||
|
|
3399 | 3401 | for i in range(dataOut.MAXNRANGENDT): |
|
|
3400 | 3402 | dataOut.range1[i]=dataOut.H0 + i*dataOut.DH |
|
|
3401 | 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 | 3405 | for j in range(dataOut.NDP): |
|
|
3404 | 3406 | dataOut.ph2[j]=0. |
|
|
3405 | 3407 | dataOut.sdp2[j]=0. |
| @@ -3505,7 +3507,7 class ElectronDensityFaraday(Operation): | |||
|
|
3505 | 3507 | ndphi=dataOut.NSHTS-4 |
|
|
3506 | 3508 | #print(dataOut.phi) |
|
|
3507 | 3509 | #exit(1) |
|
|
3508 | ''' | |
|
|
3510 | #''' | |
|
|
3509 | 3511 | if hasattr(dataOut, 'flagSpreadF') and dataOut.flagSpreadF: |
|
|
3510 | 3512 | #if dataOut.flagSpreadF: |
|
|
3511 | 3513 | nanindex = numpy.argwhere(numpy.isnan(dataOut.phi)) |
| @@ -3517,7 +3519,7 class ElectronDensityFaraday(Operation): | |||
|
|
3517 | 3519 | else: |
|
|
3518 | 3520 | #dataOut.phi_uwrp = dataOut.phi.copy() |
|
|
3519 | 3521 | dataOut.phi[:]=numpy.unwrap(dataOut.phi[:]) #Better results |
|
|
3520 | ''' | |
|
|
3522 | #''' | |
|
|
3521 | 3523 | #print(dataOut.phi) |
|
|
3522 | 3524 | #print(dataOut.ph2) |
|
|
3523 | 3525 | #exit(1) |
| @@ -3968,12 +3970,12 class NormalizeDPPowerRoberto_V2(Operation): | |||
|
|
3968 | 3970 | if hasattr(dataOut, 'flagSpreadF') and dataOut.flagSpreadF: |
|
|
3969 | 3971 | i2=int((620-dataOut.range1[0])/dataOut.DH) |
|
|
3970 | 3972 | nanindex = numpy.argwhere(numpy.isnan(dataOut.ph2)) |
|
|
3971 | print("nanindex",nanindex) | |
|
|
3973 | #print("nanindex",nanindex) | |
|
|
3972 | 3974 | i1 = nanindex[-1][0] #VER CUANDO i1>i2 |
|
|
3973 | 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 | 3976 | #print("i1, i2",i1,i2) |
|
|
3975 | print(dataOut.heightList) | |
|
|
3976 | print("Bounds: ", dataOut.heightList[i1],dataOut.heightList[i2]) | |
|
|
3977 | #print(dataOut.heightList) | |
|
|
3978 | #print("Bounds: ", dataOut.heightList[i1],dataOut.heightList[i2]) | |
|
|
3977 | 3979 | #print(dataOut.dphi[33]) |
|
|
3978 | 3980 | #print(dataOut.ph2[33]) |
|
|
3979 | 3981 | #print(dataOut.dphi[i1::]) |
| @@ -4008,9 +4010,9 class NormalizeDPPowerRoberto_V2(Operation): | |||
|
|
4008 | 4010 | #''' |
|
|
4009 | 4011 | #if (time_text.hour == 5 and time_text.minute == 32): #Year: 2022, DOY:104 |
|
|
4010 | 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 | 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 | 4016 | #if (time_text.hour == 0 and time_text.minute == 4): #Year: 2022, DOY:244 |
|
|
4015 | 4017 | #dataOut.cf = 0.08 |
|
|
4016 | 4018 | #print("here") |
| @@ -4020,12 +4022,30 class NormalizeDPPowerRoberto_V2(Operation): | |||
|
|
4020 | 4022 | #dataOut.cf = 0.09 |
|
|
4021 | 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 | 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 | 4044 | #dataOut.cf = 0.000057#0.0008136899 |
|
|
4025 | 4045 | |
|
|
4026 | 4046 | |
|
|
4027 | 4047 | dataOut.cflast[0]=dataOut.cf |
|
|
4028 | print("cf: ", dataOut.cf) | |
|
|
4048 | #print("cf: ", dataOut.cf) | |
|
|
4029 | 4049 | |
|
|
4030 | 4050 | #print(dataOut.ph2) |
|
|
4031 | 4051 | #print(dataOut.sdp2) |
| @@ -4049,7 +4069,7 class NormalizeDPPowerRoberto_V2(Operation): | |||
|
|
4049 | 4069 | #print(dataOut.ph2) |
|
|
4050 | 4070 | #print(dataOut.sdp2) |
|
|
4051 | 4071 | #input() |
|
|
4052 | print("shape before" ,numpy.shape(dataOut.ph2)) | |
|
|
4072 | #print("shape before" ,numpy.shape(dataOut.ph2)) | |
|
|
4053 | 4073 | |
|
|
4054 | 4074 | return dataOut |
|
|
4055 | 4075 | |
| @@ -4227,15 +4247,18 class DPTemperaturesEstimation(Operation): | |||
|
|
4227 | 4247 | |
|
|
4228 | 4248 | eb=numpy.resize(eb,10) |
|
|
4229 | 4249 | dataOut.ifit=numpy.resize(dataOut.ifit,10) |
|
|
4250 | #print("*********************FITACF_FIT*********************",dataOut.covinv,e,dataOut.params,eb,dataOut.m) | |
|
|
4251 | #exit(1) | |
|
|
4230 | 4252 | with suppress_stdout_stderr(): |
|
|
4231 | 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 | 4255 | if dataOut.params[2]>dataOut.params[1]*1.05: |
|
|
4234 | 4256 | dataOut.ifit[2]=0 |
|
|
4235 | 4257 | dataOut.params[1]=dataOut.params[2]=t1 |
|
|
4236 | 4258 | with suppress_stdout_stderr(): |
|
|
4237 | 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 | 4262 | if (dataOut.ifit[2]==0): |
|
|
4240 | 4263 | dataOut.params[2]=dataOut.params[1] |
|
|
4241 | 4264 | if (dataOut.ifit[3]==0 and iflag==0): |
| @@ -4306,7 +4329,7 class DPTemperaturesEstimation(Operation): | |||
|
|
4306 | 4329 | return dataOut |
|
|
4307 | 4330 | |
|
|
4308 | 4331 | |
|
|
4309 | class DenCorrection(NormalizeDPPowerRoberto_V2): | |
|
|
4332 | class DenCorrection_bckp_Dec_2022(NormalizeDPPowerRoberto_V2): | |
|
|
4310 | 4333 | ''' |
|
|
4311 | 4334 | Written by R. Flores |
|
|
4312 | 4335 | ''' |
| @@ -4322,6 +4345,7 class DenCorrection(NormalizeDPPowerRoberto_V2): | |||
|
|
4322 | 4345 | def TeTiEstimation(self,dataOut): |
|
|
4323 | 4346 | |
|
|
4324 | 4347 | y=numpy.zeros(dataOut.DPL,order='F',dtype='float32') |
|
|
4348 | ||
|
|
4325 | 4349 | #y_aux = numpy.zeros(1,,dtype='float32') |
|
|
4326 | 4350 | for i in range(dataOut.NSHTS): |
|
|
4327 | 4351 | y[0]=y[1]=dataOut.range1[i] |
| @@ -4423,7 +4447,7 class DenCorrection(NormalizeDPPowerRoberto_V2): | |||
|
|
4423 | 4447 | fion[2]=0.0 #0 |
|
|
4424 | 4448 | for j in range(dataOut.DPL): |
|
|
4425 | 4449 | tau=dataOut.alag[j]*1.0e-3 |
|
|
4426 | ||
|
|
4450 | #print("***********ACF2***********") | |
|
|
4427 | 4451 | with suppress_stdout_stderr():#The smoothness in range of "y" depends on the smoothness of the input parameters |
|
|
4428 | 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 | 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 | 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 | 4655 | class DataPlotCleaner(Operation): |
|
|
4510 | 4656 | ''' |
|
|
4511 | 4657 | Written by R. Flores |
| @@ -4825,10 +4971,16 class DataSaveCleaner(Operation): | |||
|
|
4825 | 4971 | |
|
|
4826 | 4972 | #print(time_text.hour,time_text.minute) |
|
|
4827 | 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 | 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 | 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 | 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 | 4985 | dataOut.DensityFinal[0,:]=missing |
|
|
4834 | 4986 | dataOut.EDensityFinal[0,:]=missing |
| @@ -4842,9 +4994,19 class DataSaveCleaner(Operation): | |||
|
|
4842 | 4994 | dataOut.flagNoData = True #Remueve todo el perfil |
|
|
4843 | 4995 | #''' |
|
|
4844 | 4996 | ''' |
|
|
4845 |
|
|
|
|
4846 | if (time_text.hour == 20 and time_text.minute == 8): #Year: 2022, DOY:243 | |
|
|
4847 | id_aux = 35 | |
|
|
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 | |
|
|
4998 | id_aux = 11 | |
|
|
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 | 5010 | dataOut.DensityFinal[0,id_aux:]=missing |
|
|
4849 | 5011 | dataOut.EDensityFinal[0,id_aux:]=missing |
|
|
4850 | 5012 | dataOut.ElecTempFinal[0,id_aux:]=missing |
| @@ -4853,8 +5015,9 class DataSaveCleaner(Operation): | |||
|
|
4853 | 5015 | dataOut.EIonTempFinal[0,id_aux:]=missing |
|
|
4854 | 5016 | dataOut.PhyFinal[0,id_aux:]=missing |
|
|
4855 | 5017 | dataOut.EPhyFinal[0,id_aux:]=missing |
|
|
4856 | if (time_text.hour == 20 and time_text.minute == 19): #Year: 2022, DOY:243 | |
|
|
4857 | id_aux = 33 | |
|
|
5018 | ||
|
|
5019 | if (time_text.hour == 9 and time_text.minute == 24): #Year: 2023, DOY:032 | |
|
|
5020 | id_aux = 30 | |
|
|
4858 | 5021 | dataOut.DensityFinal[0,id_aux:]=missing |
|
|
4859 | 5022 | dataOut.EDensityFinal[0,id_aux:]=missing |
|
|
4860 | 5023 | dataOut.ElecTempFinal[0,id_aux:]=missing |
| @@ -4863,8 +5026,9 class DataSaveCleaner(Operation): | |||
|
|
4863 | 5026 | dataOut.EIonTempFinal[0,id_aux:]=missing |
|
|
4864 | 5027 | dataOut.PhyFinal[0,id_aux:]=missing |
|
|
4865 | 5028 | dataOut.EPhyFinal[0,id_aux:]=missing |
|
|
4866 | if (time_text.hour == 20 and time_text.minute == 29): #Year: 2022, DOY:243 | |
|
|
4867 | id_aux = 30 | |
|
|
5029 | ||
|
|
5030 | if (time_text.hour == 5 and time_text.minute == 32): #Year: 2023, DOY:032 | |
|
|
5031 | id_aux = 41 | |
|
|
4868 | 5032 | dataOut.DensityFinal[0,id_aux:]=missing |
|
|
4869 | 5033 | dataOut.EDensityFinal[0,id_aux:]=missing |
|
|
4870 | 5034 | dataOut.ElecTempFinal[0,id_aux:]=missing |
| @@ -4873,8 +5037,9 class DataSaveCleaner(Operation): | |||
|
|
4873 | 5037 | dataOut.EIonTempFinal[0,id_aux:]=missing |
|
|
4874 | 5038 | dataOut.PhyFinal[0,id_aux:]=missing |
|
|
4875 | 5039 | dataOut.EPhyFinal[0,id_aux:]=missing |
|
|
4876 | if (time_text.hour == 20 and time_text.minute == 44): #Year: 2022, DOY:243 | |
|
|
4877 | id_aux = 31 | |
|
|
5040 | ||
|
|
5041 | if (time_text.hour == 11 and time_text.minute == 14): #Year: 2023, DOY:032 | |
|
|
5042 | id_aux = 35 | |
|
|
4878 | 5043 | dataOut.DensityFinal[0,id_aux:]=missing |
|
|
4879 | 5044 | dataOut.EDensityFinal[0,id_aux:]=missing |
|
|
4880 | 5045 | dataOut.ElecTempFinal[0,id_aux:]=missing |
| @@ -4883,8 +5048,10 class DataSaveCleaner(Operation): | |||
|
|
4883 | 5048 | dataOut.EIonTempFinal[0,id_aux:]=missing |
|
|
4884 | 5049 | dataOut.PhyFinal[0,id_aux:]=missing |
|
|
4885 | 5050 | dataOut.EPhyFinal[0,id_aux:]=missing |
|
|
4886 | if (time_text.hour <= 8): #Year: 2022, DOY:243 | |
|
|
4887 | id_aux = 11 | |
|
|
5051 | ''' | |
|
|
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 | 5055 | dataOut.DensityFinal[0,:id_aux]=missing |
|
|
4889 | 5056 | dataOut.EDensityFinal[0,:id_aux]=missing |
|
|
4890 | 5057 | dataOut.ElecTempFinal[0,:id_aux]=missing |
| @@ -4893,7 +5060,8 class DataSaveCleaner(Operation): | |||
|
|
4893 | 5060 | dataOut.EIonTempFinal[0,:id_aux]=missing |
|
|
4894 | 5061 | dataOut.PhyFinal[0,:id_aux]=missing |
|
|
4895 | 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 | 5065 | id_aux = 12 |
|
|
4898 | 5066 | dataOut.DensityFinal[0,:id_aux]=missing |
|
|
4899 | 5067 | dataOut.EDensityFinal[0,:id_aux]=missing |
| @@ -4903,28 +5071,66 class DataSaveCleaner(Operation): | |||
|
|
4903 | 5071 | dataOut.EIonTempFinal[0,:id_aux]=missing |
|
|
4904 | 5072 | dataOut.PhyFinal[0,:id_aux]=missing |
|
|
4905 | 5073 | dataOut.EPhyFinal[0,:id_aux]=missing |
|
|
4906 | if (time_text.hour == 5 and time_text.minute == 21): #Year: 2022, DOY:243 | |
|
|
4907 | id_aux = (36,37) | |
|
|
4908 | dataOut.DensityFinal[0,id_aux]=missing | |
|
|
4909 |
dataOut. |
|
|
|
4910 |
dataOut.E |
|
|
|
4911 |
dataOut. |
|
|
|
4912 |
dataOut. |
|
|
|
4913 |
dataOut. |
|
|
|
4914 |
dataOut. |
|
|
|
4915 |
dataOut. |
|
|
|
4916 | if (time_text.hour == 5 and time_text.minute == 53): #Year: 2022, DOY:243 | |
|
|
4917 | id_aux = (37,38) | |
|
|
4918 | dataOut.DensityFinal[0,id_aux]=missing | |
|
|
4919 | dataOut.EDensityFinal[0,id_aux]=missing | |
|
|
4920 |
dataOut. |
|
|
|
4921 |
dataOut.E |
|
|
|
4922 |
dataOut. |
|
|
|
4923 |
dataOut.E |
|
|
|
4924 |
dataOut. |
|
|
|
4925 |
dataOut.E |
|
|
|
4926 | if (time_text.hour == 6 and time_text.minute == 4): #Year: 2022, DOY:243 | |
|
|
4927 | id_aux = (38,39) | |
|
|
5074 | ||
|
|
5075 | if (time_text.hour == 11 and time_text.minute == 45): #Year: 2023, DOY:031 | |
|
|
5076 | id_aux = 8 | |
|
|
5077 | dataOut.DensityFinal[0,:id_aux]=missing | |
|
|
5078 | dataOut.EDensityFinal[0,:id_aux]=missing | |
|
|
5079 | dataOut.ElecTempFinal[0,:id_aux]=missing | |
|
|
5080 | dataOut.EElecTempFinal[0,:id_aux]=missing | |
|
|
5081 | dataOut.IonTempFinal[0,:id_aux]=missing | |
|
|
5082 | dataOut.EIonTempFinal[0,:id_aux]=missing | |
|
|
5083 | dataOut.PhyFinal[0,:id_aux]=missing | |
|
|
5084 | dataOut.EPhyFinal[0,:id_aux]=missing | |
|
|
5085 | ||
|
|
5086 | if (time_text.hour == 13): #Year: 2023, DOY:031 | |
|
|
5087 | id_aux = 37 | |
|
|
5088 | dataOut.DensityFinal[0,id_aux:]=missing | |
|
|
5089 | dataOut.EDensityFinal[0,id_aux:]=missing | |
|
|
5090 | dataOut.ElecTempFinal[0,id_aux:]=missing | |
|
|
5091 | dataOut.EElecTempFinal[0,id_aux:]=missing | |
|
|
5092 | dataOut.IonTempFinal[0,id_aux:]=missing | |
|
|
5093 | dataOut.EIonTempFinal[0,id_aux:]=missing | |
|
|
5094 | dataOut.PhyFinal[0,id_aux:]=missing | |
|
|
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 | 5134 | dataOut.DensityFinal[0,id_aux]=missing |
|
|
4929 | 5135 | dataOut.EDensityFinal[0,id_aux]=missing |
|
|
4930 | 5136 | dataOut.ElecTempFinal[0,id_aux]=missing |
| @@ -4933,8 +5139,9 class DataSaveCleaner(Operation): | |||
|
|
4933 | 5139 | dataOut.EIonTempFinal[0,id_aux]=missing |
|
|
4934 | 5140 | dataOut.PhyFinal[0,id_aux]=missing |
|
|
4935 | 5141 | dataOut.EPhyFinal[0,id_aux]=missing |
|
|
4936 | if (time_text.hour == 12 and time_text.minute == 6): #Year: 2022, DOY:243 | |
|
|
4937 | id_aux = (29,30) | |
|
|
5142 | ||
|
|
5143 | if (time_text.hour == 18 and time_text.minute == 41): #Year: 2023, DOY:030 | |
|
|
5144 | id_aux = 49 | |
|
|
4938 | 5145 | dataOut.DensityFinal[0,id_aux]=missing |
|
|
4939 | 5146 | dataOut.EDensityFinal[0,id_aux]=missing |
|
|
4940 | 5147 | dataOut.ElecTempFinal[0,id_aux]=missing |
| @@ -4943,8 +5150,9 class DataSaveCleaner(Operation): | |||
|
|
4943 | 5150 | dataOut.EIonTempFinal[0,id_aux]=missing |
|
|
4944 | 5151 | dataOut.PhyFinal[0,id_aux]=missing |
|
|
4945 | 5152 | dataOut.EPhyFinal[0,id_aux]=missing |
|
|
4946 | if (time_text.hour == 14 and time_text.minute == 14): #Year: 2022, DOY:243 | |
|
|
4947 | id_aux = (35,36) | |
|
|
5153 | ||
|
|
5154 | if (time_text.hour == 14 and time_text.minute == 14): #Year: 2023, DOY:030 | |
|
|
5155 | id_aux = 34 | |
|
|
4948 | 5156 | dataOut.DensityFinal[0,id_aux]=missing |
|
|
4949 | 5157 | dataOut.EDensityFinal[0,id_aux]=missing |
|
|
4950 | 5158 | dataOut.ElecTempFinal[0,id_aux]=missing |
| @@ -4953,8 +5161,9 class DataSaveCleaner(Operation): | |||
|
|
4953 | 5161 | dataOut.EIonTempFinal[0,id_aux]=missing |
|
|
4954 | 5162 | dataOut.PhyFinal[0,id_aux]=missing |
|
|
4955 | 5163 | dataOut.EPhyFinal[0,id_aux]=missing |
|
|
4956 | if (time_text.hour == 23 and time_text.minute == 2): #Year: 2022, DOY:243 | |
|
|
4957 | id_aux = (41,42) | |
|
|
5164 | ||
|
|
5165 | if (time_text.hour == 16 and time_text.minute == 1): #Year: 2023, DOY:030 | |
|
|
5166 | id_aux = (40,41) | |
|
|
4958 | 5167 | dataOut.DensityFinal[0,id_aux]=missing |
|
|
4959 | 5168 | dataOut.EDensityFinal[0,id_aux]=missing |
|
|
4960 | 5169 | dataOut.ElecTempFinal[0,id_aux]=missing |
| @@ -4963,8 +5172,9 class DataSaveCleaner(Operation): | |||
|
|
4963 | 5172 | dataOut.EIonTempFinal[0,id_aux]=missing |
|
|
4964 | 5173 | dataOut.PhyFinal[0,id_aux]=missing |
|
|
4965 | 5174 | dataOut.EPhyFinal[0,id_aux]=missing |
|
|
4966 | if (time_text.hour == 0 and time_text.minute == 8): #Year: 2022, DOY:243 | |
|
|
4967 | id_aux = 33 | |
|
|
5175 | ||
|
|
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 | 5178 | dataOut.DensityFinal[0,id_aux:]=missing |
|
|
4969 | 5179 | dataOut.EDensityFinal[0,id_aux:]=missing |
|
|
4970 | 5180 | dataOut.ElecTempFinal[0,id_aux:]=missing |
| @@ -4973,8 +5183,434 class DataSaveCleaner(Operation): | |||
|
|
4973 | 5183 | dataOut.EIonTempFinal[0,id_aux:]=missing |
|
|
4974 | 5184 | dataOut.PhyFinal[0,id_aux:]=missing |
|
|
4975 | 5185 | dataOut.EPhyFinal[0,id_aux:]=missing |
|
|
4976 | id_aux = 18 | |
|
|
4977 | dataOut.DensityFinal[0,id_aux]=missing | |
|
|
5186 | ||
|
|
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 | 5614 | dataOut.EDensityFinal[0,id_aux]=missing |
|
|
4979 | 5615 | dataOut.ElecTempFinal[0,id_aux]=missing |
|
|
4980 | 5616 | dataOut.EElecTempFinal[0,id_aux]=missing |
| @@ -5371,7 +6007,7 class DataSaveCleaner(Operation): | |||
|
|
5371 | 6007 | dataOut.PhyFinal[0,id_aux:]=missing |
|
|
5372 | 6008 | dataOut.EPhyFinal[0,id_aux:]=missing |
|
|
5373 | 6009 | ''' |
|
|
5374 |
|
|
|
|
6010 | ''' | |
|
|
5375 | 6011 | #print("den: ", dataOut.DensityFinal[0,27]) |
|
|
5376 | 6012 | if (time_text.hour == 5 and time_text.minute == 42): #Year: 2022, DOY:242 |
|
|
5377 | 6013 | id_aux = 16 |
| @@ -5493,13 +6129,13 class DataSaveCleaner(Operation): | |||
|
|
5493 | 6129 | dataOut.EIonTempFinal[0,id_aux]=missing |
|
|
5494 | 6130 | dataOut.PhyFinal[0,id_aux]=missing |
|
|
5495 | 6131 | dataOut.EPhyFinal[0,id_aux]=missing |
|
|
5496 |
|
|
|
|
6132 | ''' | |
|
|
5497 | 6133 | #print("den_final",dataOut.DensityFinal) |
|
|
5498 | 6134 | |
|
|
5499 | 6135 | |
|
|
5500 | 6136 | dataOut.flagNoData = numpy.all(numpy.isnan(dataOut.DensityFinal)) #Si todos los valores son NaN no se prosigue |
|
|
5501 | 6137 | |
|
|
5502 |
|
|
|
|
6138 | #dataOut.flagNoData = False #Descomentar solo para ploteo #Comentar para MADWriter | |
|
|
5503 | 6139 | |
|
|
5504 | 6140 | dataOut.DensityFinal *= 1.e6 #Convert units to m^⁻3 |
|
|
5505 | 6141 | dataOut.EDensityFinal *= 1.e6 #Convert units to m^⁻3 |
| @@ -6319,7 +6955,7 class removeDCHAE(Operation): | |||
|
|
6319 | 6955 | |
|
|
6320 | 6956 | return dataOut |
|
|
6321 | 6957 | |
|
|
6322 | class Decoder(Operation): | |
|
|
6958 | class Decoder_Original(Operation): | |
|
|
6323 | 6959 | |
|
|
6324 | 6960 | isConfig = False |
|
|
6325 | 6961 | __profIndex = 0 |
| @@ -6355,10 +6991,13 class Decoder(Operation): | |||
|
|
6355 | 6991 | self.__nChannels = dataOut.nChannels |
|
|
6356 | 6992 | self.__nProfiles = dataOut.nProfiles |
|
|
6357 | 6993 | self.__nHeis = dataOut.nHeights |
|
|
6358 | ||
|
|
6994 | #print("self.__nHeis: ", self.__nHeis) | |
|
|
6995 | #print("self.nBaud: ", self.nBaud) | |
|
|
6996 | #exit(1) | |
|
|
6359 | 6997 | if self.__nHeis < self.nBaud: |
|
|
6360 | 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 | 7001 | #Frequency |
|
|
6363 | 7002 | __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex) |
|
|
6364 | 7003 | |
| @@ -6421,6 +7060,10 class Decoder(Operation): | |||
|
|
6421 | 7060 | #print(profilesList) |
|
|
6422 | 7061 | for i in range(self.__nChannels): |
|
|
6423 | 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 | 7067 | self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:] |
|
|
6425 | 7068 | return self.datadecTime |
|
|
6426 | 7069 | |
| @@ -6444,7 +7087,8 class Decoder(Operation): | |||
|
|
6444 | 7087 | |
|
|
6445 | 7088 | if dataOut.flagDecodeData: |
|
|
6446 | 7089 | print("This data is already decoded, recoding again ...") |
|
|
6447 | ||
|
|
7090 | #print("code: ", numpy.shape(code)) | |
|
|
7091 | #exit(1) | |
|
|
6448 | 7092 | if not self.isConfig: |
|
|
6449 | 7093 | |
|
|
6450 | 7094 | if code is None: |
| @@ -6522,6 +7166,234 class Decoder(Operation): | |||
|
|
6522 | 7166 | return dataOut |
|
|
6523 | 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 | 7397 | class DecoderRoll(Operation): |
|
|
6526 | 7398 | |
|
|
6527 | 7399 | isConfig = False |
| @@ -6883,7 +7755,19 class ProfileSelector(Operation): | |||
|
|
6883 | 7755 | dataOut.nProfiles = len(profileList) |
|
|
6884 | 7756 | dataOut.profileIndex = dataOut.nProfiles - 1 |
|
|
6885 | 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 | 7771 | return dataOut |
|
|
6888 | 7772 | |
|
|
6889 | 7773 | """ |
| @@ -7520,7 +8404,7 class CrossProdHybrid(CrossProdDP): | |||
|
|
7520 | 8404 | #exit(1) |
|
|
7521 | 8405 | if i==0: |
|
|
7522 | 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 | 8408 | self.lagp3[n][j][self.bcounter-1]=numpy.sum(c[dataOut.NSCAN:]/self.cnorm) |
|
|
7525 | 8409 | elif i==1: |
|
|
7526 | 8410 | self.lagp1[n][j][self.bcounter-1]=numpy.sum(c[:dataOut.NSCAN]) |
| @@ -7540,8 +8424,8 class CrossProdHybrid(CrossProdDP): | |||
|
|
7540 | 8424 | #print(self.cnorm) |
|
|
7541 | 8425 | #exit(1) |
|
|
7542 | 8426 | #print("self,lagp0: ", self.lagp0[0,0,self.bcounter-1]) |
|
|
7543 | print(ww[:,0,0,self.bcounter-1]) | |
|
|
7544 | exit(1) | |
|
|
8427 | #print(ww[:,0,0,self.bcounter-1]) | |
|
|
8428 | #exit(1) | |
|
|
7545 | 8429 | |
|
|
7546 | 8430 | |
|
|
7547 | 8431 | def LP_median_estimates_original(self,dataOut): |
| @@ -7998,9 +8882,9 class LongPulseAnalysis(Operation): | |||
|
|
7998 | 8882 | #print(anoise0) |
|
|
7999 | 8883 | #print(anoise1) |
|
|
8000 | 8884 | #exit(1) |
|
|
8001 | print("nis: ", dataOut.nis) | |
|
|
8002 | print("pan: ", dataOut.pan) | |
|
|
8003 | print("pbn: ", dataOut.pbn) | |
|
|
8885 | #print("nis: ", dataOut.nis) | |
|
|
8886 | #print("pan: ", dataOut.pan) | |
|
|
8887 | #print("pbn: ", dataOut.pbn) | |
|
|
8004 | 8888 | #print(numpy.sum(dataOut.output_LP_integrated[0,:,0])) |
|
|
8005 | 8889 | ''' |
|
|
8006 | 8890 | import matplotlib.pyplot as plt |
| @@ -8008,8 +8892,8 class LongPulseAnalysis(Operation): | |||
|
|
8008 | 8892 | plt.show() |
|
|
8009 | 8893 | ''' |
|
|
8010 | 8894 | #print(dataOut.output_LP_integrated[0,40,0]) |
|
|
8011 | print(numpy.sum(dataOut.output_LP_integrated[:,0,0])) | |
|
|
8012 | exit(1) | |
|
|
8895 | #print(numpy.sum(dataOut.output_LP_integrated[:,0,0])) | |
|
|
8896 | #exit(1) | |
|
|
8013 | 8897 | |
|
|
8014 | 8898 | #################### PROBAR MÁS INTEGRACIÓN, SINO MODIFICAR VALOR DE "NIS" #################### |
|
|
8015 | 8899 | # VER dataOut.nProfiles_LP # |
| @@ -8237,7 +9121,7 class LongPulseAnalysis(Operation): | |||
|
|
8237 | 9121 | print(numpy.sum(dataOut.te2)) |
|
|
8238 | 9122 | exit(1) |
|
|
8239 | 9123 | ''' |
|
|
8240 | print("Success") | |
|
|
9124 | #print("Success 1") | |
|
|
8241 | 9125 | ###################Correlation pulse and itself |
|
|
8242 | 9126 | |
|
|
8243 | 9127 | #print(dataOut.NRANGE) |
|
|
1 | NO CONTENT: file was removed |
|
|
1 | NO CONTENT: file was removed |
|
|
1 | NO CONTENT: file was removed |
|
|
1 | NO CONTENT: file was removed |
|
|
1 | NO CONTENT: file was removed |
|
|
1 | NO CONTENT: file was removed |
|
|
1 | NO CONTENT: file was removed |
General Comments 0
You need to be logged in to leave comments.
Login now