From 2c146cb976d3c35a65fbea6d64098964d4f414e8 2023-06-01 14:38:29 From: Roberto Flores Date: 2023-06-01 14:38:29 Subject: [PATCH] last update --- diff --git a/schainf/acf2/cacai.f b/schainf/acf2/cacai.f deleted file mode 100644 index b12b057..0000000 --- a/schainf/acf2/cacai.f +++ /dev/null @@ -1,101 +0,0 @@ -*DECK CACAI - SUBROUTINE CACAI (Z, FNU, KODE, MR, N, Y, NZ, RL, TOL, ELIM, ALIM) -C***BEGIN PROLOGUE CACAI -C***SUBSIDIARY -C***PURPOSE Subsidiary to CAIRY -C***LIBRARY SLATEC -C***TYPE ALL (CACAI-A, ZACAI-A) -C***AUTHOR Amos, D. E., (SNL) -C***DESCRIPTION -C -C CACAI APPLIES THE ANALYTIC CONTINUATION FORMULA -C -C K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN) -C MP=PI*MR*CMPLX(0.0,1.0) -C -C TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT -C HALF Z PLANE FOR USE WITH CAIRY WHERE FNU=1/3 OR 2/3 AND N=1. -C CACAI IS THE SAME AS CACON WITH THE PARTS FOR LARGER ORDERS AND -C RECURRENCE REMOVED. A RECURSIVE CALL TO CACON CAN RESULT IF CACON -C IS CALLED FROM CAIRY. -C -C***SEE ALSO CAIRY -C***ROUTINES CALLED CASYI, CBKNU, CMLRI, CS1S2, CSERI, R1MACH -C***REVISION HISTORY (YYMMDD) -C 830501 DATE WRITTEN -C 910415 Prologue converted to Version 4.0 format. (BAB) -C***END PROLOGUE CACAI - COMPLEX CSGN, CSPN, C1, C2, Y, Z, ZN, CY - REAL ALIM, ARG, ASCLE, AZ, CPN, DFNU, ELIM, FMR, FNU, PI, RL, - * SGN, SPN, TOL, YY, R1MACH - INTEGER INU, IUF, KODE, MR, N, NN, NW, NZ - DIMENSION Y(N), CY(2) - DATA PI / 3.14159265358979324E0 / -C***FIRST EXECUTABLE STATEMENT CACAI - NZ = 0 - ZN = -Z - AZ = ABS(Z) - NN = N - DFNU = FNU + (N-1) - IF (AZ.LE.2.0E0) GO TO 10 - IF (AZ*AZ*0.25E0.GT.DFNU+1.0E0) GO TO 20 - 10 CONTINUE -C----------------------------------------------------------------------- -C POWER SERIES FOR THE I FUNCTION -C----------------------------------------------------------------------- - CALL CSERI(ZN, FNU, KODE, NN, Y, NW, TOL, ELIM, ALIM) - GO TO 40 - 20 CONTINUE - IF (AZ.LT.RL) GO TO 30 -C----------------------------------------------------------------------- -C ASYMPTOTIC EXPANSION FOR LARGE Z FOR THE I FUNCTION -C----------------------------------------------------------------------- - CALL CASYI(ZN, FNU, KODE, NN, Y, NW, RL, TOL, ELIM, ALIM) - IF (NW.LT.0) GO TO 70 - GO TO 40 - 30 CONTINUE -C----------------------------------------------------------------------- -C MILLER ALGORITHM NORMALIZED BY THE SERIES FOR THE I FUNCTION -C----------------------------------------------------------------------- - CALL CMLRI(ZN, FNU, KODE, NN, Y, NW, TOL) - IF(NW.LT.0) GO TO 70 - 40 CONTINUE -C----------------------------------------------------------------------- -C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION -C----------------------------------------------------------------------- - CALL CBKNU(ZN, FNU, KODE, 1, CY, NW, TOL, ELIM, ALIM) - IF (NW.NE.0) GO TO 70 - FMR = MR - SGN = -SIGN(PI,FMR) - CSGN = CMPLX(0.0E0,SGN) - IF (KODE.EQ.1) GO TO 50 - YY = -AIMAG(ZN) - CPN = COS(YY) - SPN = SIN(YY) - CSGN = CSGN*CMPLX(CPN,SPN) - 50 CONTINUE -C----------------------------------------------------------------------- -C CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE -C WHEN FNU IS LARGE -C----------------------------------------------------------------------- - INU = FNU - ARG = (FNU-INU)*SGN - CPN = COS(ARG) - SPN = SIN(ARG) - CSPN = CMPLX(CPN,SPN) - IF (MOD(INU,2).EQ.1) CSPN = -CSPN - C1 = CY(1) - C2 = Y(1) - IF (KODE.EQ.1) GO TO 60 - IUF = 0 - ASCLE = 1.0E+3*R1MACH(1)/TOL - CALL CS1S2(ZN, C1, C2, NW, ASCLE, ALIM, IUF) - NZ = NZ + NW - 60 CONTINUE - Y(1) = CSPN*C1 + CSGN*C2 - RETURN - 70 CONTINUE - NZ = -1 - IF(NW.EQ.(-2)) NZ=-2 - RETURN - END diff --git a/schainf/acf2/cairy.f b/schainf/acf2/cairy.f deleted file mode 100644 index 4fe2f7e..0000000 --- a/schainf/acf2/cairy.f +++ /dev/null @@ -1,342 +0,0 @@ -*DECK CAIRY - SUBROUTINE CAIRY (Z, ID, KODE, AI, NZ, IERR) -C***BEGIN PROLOGUE CAIRY -C***PURPOSE Compute the Airy function Ai(z) or its derivative dAi/dz -C for complex argument z. A scaling option is available -C to help avoid underflow and overflow. -C***LIBRARY SLATEC -C***CATEGORY C10D -C***TYPE COMPLEX (CAIRY-C, ZAIRY-C) -C***KEYWORDS AIRY FUNCTION, BESSEL FUNCTION OF ORDER ONE THIRD, -C BESSEL FUNCTION OF ORDER TWO THIRDS -C***AUTHOR Amos, D. E., (SNL) -C***DESCRIPTION -C -C On KODE=1, CAIRY computes the complex Airy function Ai(z) -C or its derivative dAi/dz on ID=0 or ID=1 respectively. On -C KODE=2, a scaling option exp(zeta)*Ai(z) or exp(zeta)*dAi/dz -C is provided to remove the exponential decay in -pi/31 and from power series when abs(z)<=1. -C -C In most complex variable computation, one must evaluate ele- -C mentary functions. When the magnitude of Z is large, losses -C of significance by argument reduction occur. Consequently, if -C the magnitude of ZETA=(2/3)*Z**(3/2) exceeds U1=SQRT(0.5/UR), -C then losses exceeding half precision are likely and an error -C flag IERR=3 is triggered where UR=R1MACH(4)=UNIT ROUNDOFF. -C Also, if the magnitude of ZETA is larger than U2=0.5/UR, then -C all significance is lost and IERR=4. In order to use the INT -C function, ZETA must be further restricted not to exceed -C U3=I1MACH(9)=LARGEST INTEGER. Thus, the magnitude of ZETA -C must be restricted by MIN(U2,U3). In IEEE arithmetic, U1,U2, -C and U3 are approximately 2.0E+3, 4.2E+6, 2.1E+9 in single -C precision and 4.7E+7, 2.3E+15, 2.1E+9 in double precision. -C This makes U2 limiting is single precision and U3 limiting -C in double precision. This means that the magnitude of Z -C cannot exceed approximately 3.4E+4 in single precision and -C 2.1E+6 in double precision. This also means that one can -C expect to retain, in the worst cases on 32-bit machines, -C no digits in single precision and only 6 digits in double -C precision. -C -C The approximate relative error in the magnitude of a complex -C Bessel function can be expressed as P*10**S where P=MAX(UNIT -C ROUNDOFF,1.0E-18) is the nominal precision and 10**S repre- -C sents the increase in error due to argument reduction in the -C elementary functions. Here, S=MAX(1,ABS(LOG10(ABS(Z))), -C ABS(LOG10(FNU))) approximately (i.e., S=MAX(1,ABS(EXPONENT OF -C ABS(Z),ABS(EXPONENT OF FNU)) ). However, the phase angle may -C have only absolute accuracy. This is most likely to occur -C when one component (in magnitude) is larger than the other by -C several orders of magnitude. If one component is 10**K larger -C than the other, then one can expect only MAX(ABS(LOG10(P))-K, -C 0) significant digits; or, stated another way, when K exceeds -C the exponent of P, no significant digits remain in the smaller -C component. However, the phase angle retains absolute accuracy -C because, in complex arithmetic with precision P, the smaller -C component will not (as a rule) decrease below P times the -C magnitude of the larger component. In these extreme cases, -C the principal phase angle is on the order of +P, -P, PI/2-P, -C or -PI/2+P. -C -C***REFERENCES 1. M. Abramowitz and I. A. Stegun, Handbook of Mathe- -C matical Functions, National Bureau of Standards -C Applied Mathematics Series 55, U. S. Department -C of Commerce, Tenth Printing (1972) or later. -C 2. D. E. Amos, Computation of Bessel Functions of -C Complex Argument and Large Order, Report SAND83-0643, -C Sandia National Laboratories, Albuquerque, NM, May -C 1983. -C 3. D. E. Amos, A Subroutine Package for Bessel Functions -C of a Complex Argument and Nonnegative Order, Report -C SAND85-1018, Sandia National Laboratory, Albuquerque, -C NM, May 1985. -C 4. D. E. Amos, A portable package for Bessel functions -C of a complex argument and nonnegative order, ACM -C Transactions on Mathematical Software, 12 (September -C 1986), pp. 265-273. -C -C***ROUTINES CALLED CACAI, CBKNU, I1MACH, R1MACH -C***REVISION HISTORY (YYMMDD) -C 830501 DATE WRITTEN -C 890801 REVISION DATE from Version 3.2 -C 910415 Prologue converted to Version 4.0 format. (BAB) -C 920128 Category corrected. (WRB) -C 920811 Prologue revised. (DWL) -C***END PROLOGUE CAIRY - COMPLEX AI, CONE, CSQ, CY, S1, S2, TRM1, TRM2, Z, ZTA, Z3 - REAL AA, AD, AK, ALIM, ATRM, AZ, AZ3, BK, CK, COEF, C1, C2, DIG, - * DK, D1, D2, ELIM, FID, FNU, RL, R1M5, SFAC, TOL, TTH, ZI, ZR, - * Z3I, Z3R, R1MACH, BB, ALAZ - INTEGER ID, IERR, IFLAG, K, KODE, K1, K2, MR, NN, NZ, I1MACH - DIMENSION CY(1) - DATA TTH, C1, C2, COEF /6.66666666666666667E-01, - * 3.55028053887817240E-01,2.58819403792806799E-01, - * 1.83776298473930683E-01/ - DATA CONE / (1.0E0,0.0E0) / -C***FIRST EXECUTABLE STATEMENT CAIRY - IERR = 0 - NZ=0 - IF (ID.LT.0 .OR. ID.GT.1) IERR=1 - IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 - IF (IERR.NE.0) RETURN - AZ = ABS(Z) - TOL = MAX(R1MACH(4),1.0E-18) - FID = ID - IF (AZ.GT.1.0E0) GO TO 60 -C----------------------------------------------------------------------- -C POWER SERIES FOR ABS(Z).LE.1. -C----------------------------------------------------------------------- - S1 = CONE - S2 = CONE - IF (AZ.LT.TOL) GO TO 160 - AA = AZ*AZ - IF (AA.LT.TOL/AZ) GO TO 40 - TRM1 = CONE - TRM2 = CONE - ATRM = 1.0E0 - Z3 = Z*Z*Z - AZ3 = AZ*AA - AK = 2.0E0 + FID - BK = 3.0E0 - FID - FID - CK = 4.0E0 - FID - DK = 3.0E0 + FID + FID - D1 = AK*DK - D2 = BK*CK - AD = MIN(D1,D2) - AK = 24.0E0 + 9.0E0*FID - BK = 30.0E0 - 9.0E0*FID - Z3R = REAL(Z3) - Z3I = AIMAG(Z3) - DO 30 K=1,25 - TRM1 = TRM1*CMPLX(Z3R/D1,Z3I/D1) - S1 = S1 + TRM1 - TRM2 = TRM2*CMPLX(Z3R/D2,Z3I/D2) - S2 = S2 + TRM2 - ATRM = ATRM*AZ3/AD - D1 = D1 + AK - D2 = D2 + BK - AD = MIN(D1,D2) - IF (ATRM.LT.TOL*AD) GO TO 40 - AK = AK + 18.0E0 - BK = BK + 18.0E0 - 30 CONTINUE - 40 CONTINUE - IF (ID.EQ.1) GO TO 50 - AI = S1*CMPLX(C1,0.0E0) - Z*S2*CMPLX(C2,0.0E0) - IF (KODE.EQ.1) RETURN - ZTA = Z*CSQRT(Z)*CMPLX(TTH,0.0E0) - AI = AI*CEXP(ZTA) - RETURN - 50 CONTINUE - AI = -S2*CMPLX(C2,0.0E0) - IF (AZ.GT.TOL) AI = AI + Z*Z*S1*CMPLX(C1/(1.0E0+FID),0.0E0) - IF (KODE.EQ.1) RETURN - ZTA = Z*CSQRT(Z)*CMPLX(TTH,0.0E0) - AI = AI*CEXP(ZTA) - RETURN -C----------------------------------------------------------------------- -C CASE FOR ABS(Z).GT.1.0 -C----------------------------------------------------------------------- - 60 CONTINUE - FNU = (1.0E0+FID)/3.0E0 -C----------------------------------------------------------------------- -C SET PARAMETERS RELATED TO MACHINE CONSTANTS. -C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. -C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. -C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND -C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR -C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. -C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. -C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). -C----------------------------------------------------------------------- - K1 = I1MACH(12) - K2 = I1MACH(13) - R1M5 = R1MACH(5) - K = MIN(ABS(K1),ABS(K2)) - ELIM = 2.303E0*(K*R1M5-3.0E0) - K1 = I1MACH(11) - 1 - AA = R1M5*K1 - DIG = MIN(AA,18.0E0) - AA = AA*2.303E0 - ALIM = ELIM + MAX(-AA,-41.45E0) - RL = 1.2E0*DIG + 3.0E0 - ALAZ=ALOG(AZ) -C----------------------------------------------------------------------- -C TEST FOR RANGE -C----------------------------------------------------------------------- - AA=0.5E0/TOL - BB=I1MACH(9)*0.5E0 - AA=MIN(AA,BB) - AA=AA**TTH - IF (AZ.GT.AA) GO TO 260 - AA=SQRT(AA) - IF (AZ.GT.AA) IERR=3 - CSQ=CSQRT(Z) - ZTA=Z*CSQ*CMPLX(TTH,0.0E0) -C----------------------------------------------------------------------- -C RE(ZTA).LE.0 WHEN RE(Z).LT.0, ESPECIALLY WHEN IM(Z) IS SMALL -C----------------------------------------------------------------------- - IFLAG = 0 - SFAC = 1.0E0 - ZI = AIMAG(Z) - ZR = REAL(Z) - AK = AIMAG(ZTA) - IF (ZR.GE.0.0E0) GO TO 70 - BK = REAL(ZTA) - CK = -ABS(BK) - ZTA = CMPLX(CK,AK) - 70 CONTINUE - IF (ZI.NE.0.0E0) GO TO 80 - IF (ZR.GT.0.0E0) GO TO 80 - ZTA = CMPLX(0.0E0,AK) - 80 CONTINUE - AA = REAL(ZTA) - IF (AA.GE.0.0E0 .AND. ZR.GT.0.0E0) GO TO 100 - IF (KODE.EQ.2) GO TO 90 -C----------------------------------------------------------------------- -C OVERFLOW TEST -C----------------------------------------------------------------------- - IF (AA.GT.(-ALIM)) GO TO 90 - AA = -AA + 0.25E0*ALAZ - IFLAG = 1 - SFAC = TOL - IF (AA.GT.ELIM) GO TO 240 - 90 CONTINUE -C----------------------------------------------------------------------- -C CBKNU AND CACAI RETURN EXP(ZTA)*K(FNU,ZTA) ON KODE=2 -C----------------------------------------------------------------------- - MR = 1 - IF (ZI.LT.0.0E0) MR = -1 - CALL CACAI(ZTA, FNU, KODE, MR, 1, CY, NN, RL, TOL, ELIM, ALIM) - IF (NN.LT.0) GO TO 250 - NZ = NZ + NN - GO TO 120 - 100 CONTINUE - IF (KODE.EQ.2) GO TO 110 -C----------------------------------------------------------------------- -C UNDERFLOW TEST -C----------------------------------------------------------------------- - IF (AA.LT.ALIM) GO TO 110 - AA = -AA - 0.25E0*ALAZ - IFLAG = 2 - SFAC = 1.0E0/TOL - IF (AA.LT.(-ELIM)) GO TO 180 - 110 CONTINUE - CALL CBKNU(ZTA, FNU, KODE, 1, CY, NZ, TOL, ELIM, ALIM) - 120 CONTINUE - S1 = CY(1)*CMPLX(COEF,0.0E0) - IF (IFLAG.NE.0) GO TO 140 - IF (ID.EQ.1) GO TO 130 - AI = CSQ*S1 - RETURN - 130 AI = -Z*S1 - RETURN - 140 CONTINUE - S1 = S1*CMPLX(SFAC,0.0E0) - IF (ID.EQ.1) GO TO 150 - S1 = S1*CSQ - AI = S1*CMPLX(1.0E0/SFAC,0.0E0) - RETURN - 150 CONTINUE - S1 = -S1*Z - AI = S1*CMPLX(1.0E0/SFAC,0.0E0) - RETURN - 160 CONTINUE - AA = 1.0E+3*R1MACH(1) - S1 = CMPLX(0.0E0,0.0E0) - IF (ID.EQ.1) GO TO 170 - IF (AZ.GT.AA) S1 = CMPLX(C2,0.0E0)*Z - AI = CMPLX(C1,0.0E0) - S1 - RETURN - 170 CONTINUE - AI = -CMPLX(C2,0.0E0) - AA = SQRT(AA) - IF (AZ.GT.AA) S1 = Z*Z*CMPLX(0.5E0,0.0E0) - AI = AI + S1*CMPLX(C1,0.0E0) - RETURN - 180 CONTINUE - NZ = 1 - AI = CMPLX(0.0E0,0.0E0) - RETURN - 240 CONTINUE - NZ = 0 - IERR=2 - RETURN - 250 CONTINUE - IF(NN.EQ.(-1)) GO TO 240 - NZ=0 - IERR=5 - RETURN - 260 CONTINUE - IERR=4 - NZ=0 - RETURN - END diff --git a/schainf/acf2/casyi.f b/schainf/acf2/casyi.f deleted file mode 100644 index fdb2ee2..0000000 --- a/schainf/acf2/casyi.f +++ /dev/null @@ -1,136 +0,0 @@ -*DECK CASYI - SUBROUTINE CASYI (Z, FNU, KODE, N, Y, NZ, RL, TOL, ELIM, ALIM) -C***BEGIN PROLOGUE CASYI -C***SUBSIDIARY -C***PURPOSE Subsidiary to CBESI and CBESK -C***LIBRARY SLATEC -C***TYPE ALL (CASYI-A, ZASYI-A) -C***AUTHOR Amos, D. E., (SNL) -C***DESCRIPTION -C -C CASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY -C MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE ABS(Z) IN THE -C REGION ABS(Z).GT.MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN. -C NZ.LT.0 INDICATES AN OVERFLOW ON KODE=1. -C -C***SEE ALSO CBESI, CBESK -C***ROUTINES CALLED R1MACH -C***REVISION HISTORY (YYMMDD) -C 830501 DATE WRITTEN -C 910415 Prologue converted to Version 4.0 format. (BAB) -C***END PROLOGUE CASYI - COMPLEX AK1, CK, CONE, CS1, CS2, CZ, CZERO, DK, EZ, P1, RZ, S2, - * Y, Z - REAL AA, ACZ, AEZ, AK, ALIM, ARG, ARM, ATOL, AZ, BB, BK, DFNU, - * DNU2, ELIM, FDN, FNU, PI, RL, RTPI, RTR1, S, SGN, SQK, TOL, X, - * YY, R1MACH - INTEGER I, IB, IL, INU, J, JL, K, KODE, KODED, M, N, NN, NZ - DIMENSION Y(N) - DATA PI, RTPI /3.14159265358979324E0 , 0.159154943091895336E0 / - DATA CZERO, CONE / (0.0E0,0.0E0), (1.0E0,0.0E0) / -C***FIRST EXECUTABLE STATEMENT CASYI - NZ = 0 - AZ = ABS(Z) - X = REAL(Z) - ARM = 1.0E+3*R1MACH(1) - RTR1 = SQRT(ARM) - IL = MIN(2,N) - DFNU = FNU + (N-IL) -C----------------------------------------------------------------------- -C OVERFLOW TEST -C----------------------------------------------------------------------- - AK1 = CMPLX(RTPI,0.0E0)/Z - AK1 = CSQRT(AK1) - CZ = Z - IF (KODE.EQ.2) CZ = Z - CMPLX(X,0.0E0) - ACZ = REAL(CZ) - IF (ABS(ACZ).GT.ELIM) GO TO 80 - DNU2 = DFNU + DFNU - KODED = 1 - IF ((ABS(ACZ).GT.ALIM) .AND. (N.GT.2)) GO TO 10 - KODED = 0 - AK1 = AK1*CEXP(CZ) - 10 CONTINUE - FDN = 0.0E0 - IF (DNU2.GT.RTR1) FDN = DNU2*DNU2 - EZ = Z*CMPLX(8.0E0,0.0E0) -C----------------------------------------------------------------------- -C WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE -C FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE -C EXPANSION FOR THE IMAGINARY PART. -C----------------------------------------------------------------------- - AEZ = 8.0E0*AZ - S = TOL/AEZ - JL = RL+RL + 2 - YY = AIMAG(Z) - P1 = CZERO - IF (YY.EQ.0.0E0) GO TO 20 -C----------------------------------------------------------------------- -C CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF -C SIGNIFICANCE WHEN FNU OR N IS LARGE -C----------------------------------------------------------------------- - INU = FNU - ARG = (FNU-INU)*PI - INU = INU + N - IL - AK = -SIN(ARG) - BK = COS(ARG) - IF (YY.LT.0.0E0) BK = -BK - P1 = CMPLX(AK,BK) - IF (MOD(INU,2).EQ.1) P1 = -P1 - 20 CONTINUE - DO 50 K=1,IL - SQK = FDN - 1.0E0 - ATOL = S*ABS(SQK) - SGN = 1.0E0 - CS1 = CONE - CS2 = CONE - CK = CONE - AK = 0.0E0 - AA = 1.0E0 - BB = AEZ - DK = EZ - DO 30 J=1,JL - CK = CK*CMPLX(SQK,0.0E0)/DK - CS2 = CS2 + CK - SGN = -SGN - CS1 = CS1 + CK*CMPLX(SGN,0.0E0) - DK = DK + EZ - AA = AA*ABS(SQK)/BB - BB = BB + AEZ - AK = AK + 8.0E0 - SQK = SQK - AK - IF (AA.LE.ATOL) GO TO 40 - 30 CONTINUE - GO TO 90 - 40 CONTINUE - S2 = CS1 - IF (X+X.LT.ELIM) S2 = S2 + P1*CS2*CEXP(-Z-Z) - FDN = FDN + 8.0E0*DFNU + 4.0E0 - P1 = -P1 - M = N - IL + K - Y(M) = S2*AK1 - 50 CONTINUE - IF (N.LE.2) RETURN - NN = N - K = NN - 2 - AK = K - RZ = (CONE+CONE)/Z - IB = 3 - DO 60 I=IB,NN - Y(K) = CMPLX(AK+FNU,0.0E0)*RZ*Y(K+1) + Y(K+2) - AK = AK - 1.0E0 - K = K - 1 - 60 CONTINUE - IF (KODED.EQ.0) RETURN - CK = CEXP(CZ) - DO 70 I=1,NN - Y(I) = Y(I)*CK - 70 CONTINUE - RETURN - 80 CONTINUE - NZ = -1 - RETURN - 90 CONTINUE - NZ=-2 - RETURN - END diff --git a/schainf/acf2/cbesi.f b/schainf/acf2/cbesi.f deleted file mode 100644 index f99bd30..0000000 --- a/schainf/acf2/cbesi.f +++ /dev/null @@ -1,261 +0,0 @@ -*DECK CBESI - SUBROUTINE CBESI (Z, FNU, KODE, N, CY, NZ, IERR) -C***BEGIN PROLOGUE CBESI -C***PURPOSE Compute a sequence of the Bessel functions I(a,z) for -C complex argument z and real nonnegative orders a=b,b+1, -C b+2,... where b>0. A scaling option is available to -C help avoid overflow. -C***LIBRARY SLATEC -C***CATEGORY C10B4 -C***TYPE COMPLEX (CBESI-C, ZBESI-C) -C***KEYWORDS BESSEL FUNCTIONS OF COMPLEX ARGUMENT, I BESSEL FUNCTIONS, -C MODIFIED BESSEL FUNCTIONS -C***AUTHOR Amos, D. E., (SNL) -C***DESCRIPTION -C -C On KODE=1, CBESI computes an N-member sequence of complex -C Bessel functions CY(L)=I(FNU+L-1,Z) for real nonnegative -C orders FNU+L-1, L=1,...,N and complex Z in the cut plane -C -pi=0 -C KODE - A parameter to indicate the scaling option -C KODE=1 returns -C CY(L)=I(FNU+L-1,Z), L=1,...,N -C =2 returns -C CY(L)=exp(-abs(X))*I(FNU+L-1,Z), L=1,...,N -C where X=Re(Z) -C N - Number of terms in the sequence, N>=1 -C -C Output -C CY - Result vector of type COMPLEX -C NZ - Number of underflows set to zero -C NZ=0 Normal return -C NZ>0 CY(L)=0, L=N-NZ+1,...,N -C IERR - Error flag -C IERR=0 Normal return - COMPUTATION COMPLETED -C IERR=1 Input error - NO COMPUTATION -C IERR=2 Overflow - NO COMPUTATION -C (Re(Z) too large on KODE=1) -C IERR=3 Precision warning - COMPUTATION COMPLETED -C (Result has half precision or less -C because abs(Z) or FNU+N-1 is large) -C IERR=4 Precision error - NO COMPUTATION -C (Result has no precision because -C abs(Z) or FNU+N-1 is too large) -C IERR=5 Algorithmic error - NO COMPUTATION -C (Termination condition not met) -C -C *Long Description: -C -C The computation of I(a,z) is carried out by the power series -C for small abs(z), the asymptotic expansion for large abs(z), -C the Miller algorithm normalized by the Wronskian and a -C Neumann series for intermediate magnitudes of z, and the -C uniform asymptotic expansions for I(a,z) and J(a,z) for -C large orders a. Backward recurrence is used to generate -C sequences or reduce orders when necessary. -C -C The calculations above are done in the right half plane and -C continued into the left half plane by the formula -C -C I(a,z*exp(t)) = exp(t*a)*I(a,z), Re(z)>0 -C t = i*pi or -i*pi -C -C For negative orders, the formula -C -C I(-a,z) = I(a,z) + (2/pi)*sin(pi*a)*K(a,z) -C -C can be used. However, for large orders close to integers the -C the function changes radically. When a is a large positive -C integer, the magnitude of I(-a,z)=I(a,z) is a large -C negative power of ten. But when a is not an integer, -C K(a,z) dominates in magnitude with a large positive power of -C ten and the most that the second term can be reduced is by -C unit roundoff from the coefficient. Thus, wide changes can -C occur within unit roundoff of a large integer for a. Here, -C large means a>abs(z). -C -C In most complex variable computation, one must evaluate ele- -C mentary functions. When the magnitude of Z or FNU+N-1 is -C large, losses of significance by argument reduction occur. -C Consequently, if either one exceeds U1=SQRT(0.5/UR), then -C losses exceeding half precision are likely and an error flag -C IERR=3 is triggered where UR=R1MACH(4)=UNIT ROUNDOFF. Also, -C if either is larger than U2=0.5/UR, then all significance is -C lost and IERR=4. In order to use the INT function, arguments -C must be further restricted not to exceed the largest machine -C integer, U3=I1MACH(9). Thus, the magnitude of Z and FNU+N-1 -C is restricted by MIN(U2,U3). In IEEE arithmetic, U1,U2, and -C U3 approximate 2.0E+3, 4.2E+6, 2.1E+9 in single precision -C and 4.7E+7, 2.3E+15 and 2.1E+9 in double precision. This -C makes U2 limiting in single precision and U3 limiting in -C double precision. This means that one can expect to retain, -C in the worst cases on IEEE machines, no digits in single pre- -C cision and only 6 digits in double precision. Similar con- -C siderations hold for other machines. -C -C The approximate relative error in the magnitude of a complex -C Bessel function can be expressed as P*10**S where P=MAX(UNIT -C ROUNDOFF,1.0E-18) is the nominal precision and 10**S repre- -C sents the increase in error due to argument reduction in the -C elementary functions. Here, S=MAX(1,ABS(LOG10(ABS(Z))), -C ABS(LOG10(FNU))) approximately (i.e., S=MAX(1,ABS(EXPONENT OF -C ABS(Z),ABS(EXPONENT OF FNU)) ). However, the phase angle may -C have only absolute accuracy. This is most likely to occur -C when one component (in magnitude) is larger than the other by -C several orders of magnitude. If one component is 10**K larger -C than the other, then one can expect only MAX(ABS(LOG10(P))-K, -C 0) significant digits; or, stated another way, when K exceeds -C the exponent of P, no significant digits remain in the smaller -C component. However, the phase angle retains absolute accuracy -C because, in complex arithmetic with precision P, the smaller -C component will not (as a rule) decrease below P times the -C magnitude of the larger component. In these extreme cases, -C the principal phase angle is on the order of +P, -P, PI/2-P, -C or -PI/2+P. -C -C***REFERENCES 1. M. Abramowitz and I. A. Stegun, Handbook of Mathe- -C matical Functions, National Bureau of Standards -C Applied Mathematics Series 55, U. S. Department -C of Commerce, Tenth Printing (1972) or later. -C 2. D. E. Amos, Computation of Bessel Functions of -C Complex Argument, Report SAND83-0086, Sandia National -C Laboratories, Albuquerque, NM, May 1983. -C 3. D. E. Amos, Computation of Bessel Functions of -C Complex Argument and Large Order, Report SAND83-0643, -C Sandia National Laboratories, Albuquerque, NM, May -C 1983. -C 4. D. E. Amos, A Subroutine Package for Bessel Functions -C of a Complex Argument and Nonnegative Order, Report -C SAND85-1018, Sandia National Laboratory, Albuquerque, -C NM, May 1985. -C 5. D. E. Amos, A portable package for Bessel functions -C of a complex argument and nonnegative order, ACM -C Transactions on Mathematical Software, 12 (September -C 1986), pp. 265-273. -C -C***ROUTINES CALLED CBINU, I1MACH, R1MACH -C***REVISION HISTORY (YYMMDD) -C 830501 DATE WRITTEN -C 890801 REVISION DATE from Version 3.2 -C 910415 Prologue converted to Version 4.0 format. (BAB) -C 920128 Category corrected. (WRB) -C 920811 Prologue revised. (DWL) -C***END PROLOGUE CBESI - COMPLEX CONE, CSGN, CY, Z, ZN - REAL AA, ALIM, ARG, DIG, ELIM, FNU, FNUL, PI, RL, R1M5, S1, S2, - * TOL, XX, YY, R1MACH, AZ, FN, BB, ASCLE, RTOL, ATOL - INTEGER I, IERR, INU, K, KODE, K1, K2, N, NN, NZ, I1MACH - DIMENSION CY(N) - DATA PI /3.14159265358979324E0/ - DATA CONE / (1.0E0,0.0E0) / -C -C***FIRST EXECUTABLE STATEMENT CBESI - IERR = 0 - NZ=0 - IF (FNU.LT.0.0E0) IERR=1 - IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 - IF (N.LT.1) IERR=1 - IF (IERR.NE.0) RETURN - XX = REAL(Z) - YY = AIMAG(Z) -C----------------------------------------------------------------------- -C SET PARAMETERS RELATED TO MACHINE CONSTANTS. -C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. -C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. -C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND -C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR -C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. -C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. -C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). -C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU. -C----------------------------------------------------------------------- - TOL = MAX(R1MACH(4),1.0E-18) - K1 = I1MACH(12) - K2 = I1MACH(13) - R1M5 = R1MACH(5) - K = MIN(ABS(K1),ABS(K2)) - ELIM = 2.303E0*(K*R1M5-3.0E0) - K1 = I1MACH(11) - 1 - AA = R1M5*K1 - DIG = MIN(AA,18.0E0) - AA = AA*2.303E0 - ALIM = ELIM + MAX(-AA,-41.45E0) - RL = 1.2E0*DIG + 3.0E0 - FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0) - AZ = ABS(Z) -C----------------------------------------------------------------------- -C TEST FOR RANGE -C----------------------------------------------------------------------- - AA = 0.5E0/TOL - BB=I1MACH(9)*0.5E0 - AA=MIN(AA,BB) - IF(AZ.GT.AA) GO TO 140 - FN=FNU+(N-1) - IF(FN.GT.AA) GO TO 140 - AA=SQRT(AA) - IF(AZ.GT.AA) IERR=3 - IF(FN.GT.AA) IERR=3 - ZN = Z - CSGN = CONE - IF (XX.GE.0.0E0) GO TO 40 - ZN = -Z -C----------------------------------------------------------------------- -C CALCULATE CSGN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE -C WHEN FNU IS LARGE -C----------------------------------------------------------------------- - INU = FNU - ARG = (FNU-INU)*PI - IF (YY.LT.0.0E0) ARG = -ARG - S1 = COS(ARG) - S2 = SIN(ARG) - CSGN = CMPLX(S1,S2) - IF (MOD(INU,2).EQ.1) CSGN = -CSGN - 40 CONTINUE - CALL CBINU(ZN, FNU, KODE, N, CY, NZ, RL, FNUL, TOL, ELIM, ALIM) - IF (NZ.LT.0) GO TO 120 - IF (XX.GE.0.0E0) RETURN -C----------------------------------------------------------------------- -C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE -C----------------------------------------------------------------------- - NN = N - NZ - IF (NN.EQ.0) RETURN - RTOL = 1.0E0/TOL - ASCLE = R1MACH(1)*RTOL*1.0E+3 - DO 50 I=1,NN -C CY(I) = CY(I)*CSGN - ZN=CY(I) - AA=REAL(ZN) - BB=AIMAG(ZN) - ATOL=1.0E0 - IF (MAX(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 55 - ZN = ZN*CMPLX(RTOL,0.0E0) - ATOL = TOL - 55 CONTINUE - ZN = ZN*CSGN - CY(I) = ZN*CMPLX(ATOL,0.0E0) - CSGN = -CSGN - 50 CONTINUE - RETURN - 120 CONTINUE - IF(NZ.EQ.(-2)) GO TO 130 - NZ = 0 - IERR=2 - RETURN - 130 CONTINUE - NZ=0 - IERR=5 - RETURN - 140 CONTINUE - NZ=0 - IERR=4 - RETURN - END diff --git a/schainf/acf2/cbinu.f b/schainf/acf2/cbinu.f deleted file mode 100644 index 8f0e830..0000000 --- a/schainf/acf2/cbinu.f +++ /dev/null @@ -1,115 +0,0 @@ -*DECK CBINU - SUBROUTINE CBINU (Z, FNU, KODE, N, CY, NZ, RL, FNUL, TOL, ELIM, - + ALIM) -C***BEGIN PROLOGUE CBINU -C***SUBSIDIARY -C***PURPOSE Subsidiary to CAIRY, CBESH, CBESI, CBESJ, CBESK and CBIRY -C***LIBRARY SLATEC -C***TYPE ALL (CBINU-A, ZBINU-A) -C***AUTHOR Amos, D. E., (SNL) -C***DESCRIPTION -C -C CBINU COMPUTES THE I FUNCTION IN THE RIGHT HALF Z PLANE -C -C***SEE ALSO CAIRY, CBESH, CBESI, CBESJ, CBESK, CBIRY -C***ROUTINES CALLED CASYI, CBUNI, CMLRI, CSERI, CUOIK, CWRSK -C***REVISION HISTORY (YYMMDD) -C 830501 DATE WRITTEN -C 910415 Prologue converted to Version 4.0 format. (BAB) -C***END PROLOGUE CBINU - COMPLEX CW, CY, CZERO, Z - REAL ALIM, AZ, DFNU, ELIM, FNU, FNUL, RL, TOL - INTEGER I, INW, KODE, N, NLAST, NN, NUI, NW, NZ - DIMENSION CY(N), CW(2) - DATA CZERO / (0.0E0,0.0E0) / -C***FIRST EXECUTABLE STATEMENT CBINU - NZ = 0 - AZ = ABS(Z) - NN = N - DFNU = FNU + (N-1) - IF (AZ.LE.2.0E0) GO TO 10 - IF (AZ*AZ*0.25E0.GT.DFNU+1.0E0) GO TO 20 - 10 CONTINUE -C----------------------------------------------------------------------- -C POWER SERIES -C----------------------------------------------------------------------- - CALL CSERI(Z, FNU, KODE, NN, CY, NW, TOL, ELIM, ALIM) - INW = ABS(NW) - NZ = NZ + INW - NN = NN - INW - IF (NN.EQ.0) RETURN - IF (NW.GE.0) GO TO 120 - DFNU = FNU + (NN-1) - 20 CONTINUE - IF (AZ.LT.RL) GO TO 40 - IF (DFNU.LE.1.0E0) GO TO 30 - IF (AZ+AZ.LT.DFNU*DFNU) GO TO 50 -C----------------------------------------------------------------------- -C ASYMPTOTIC EXPANSION FOR LARGE Z -C----------------------------------------------------------------------- - 30 CONTINUE - CALL CASYI(Z, FNU, KODE, NN, CY, NW, RL, TOL, ELIM, ALIM) - IF (NW.LT.0) GO TO 130 - GO TO 120 - 40 CONTINUE - IF (DFNU.LE.1.0E0) GO TO 70 - 50 CONTINUE -C----------------------------------------------------------------------- -C OVERFLOW AND UNDERFLOW TEST ON I SEQUENCE FOR MILLER ALGORITHM -C----------------------------------------------------------------------- - CALL CUOIK(Z, FNU, KODE, 1, NN, CY, NW, TOL, ELIM, ALIM) - IF (NW.LT.0) GO TO 130 - NZ = NZ + NW - NN = NN - NW - IF (NN.EQ.0) RETURN - DFNU = FNU+(NN-1) - IF (DFNU.GT.FNUL) GO TO 110 - IF (AZ.GT.FNUL) GO TO 110 - 60 CONTINUE - IF (AZ.GT.RL) GO TO 80 - 70 CONTINUE -C----------------------------------------------------------------------- -C MILLER ALGORITHM NORMALIZED BY THE SERIES -C----------------------------------------------------------------------- - CALL CMLRI(Z, FNU, KODE, NN, CY, NW, TOL) - IF(NW.LT.0) GO TO 130 - GO TO 120 - 80 CONTINUE -C----------------------------------------------------------------------- -C MILLER ALGORITHM NORMALIZED BY THE WRONSKIAN -C----------------------------------------------------------------------- -C----------------------------------------------------------------------- -C OVERFLOW TEST ON K FUNCTIONS USED IN WRONSKIAN -C----------------------------------------------------------------------- - CALL CUOIK(Z, FNU, KODE, 2, 2, CW, NW, TOL, ELIM, ALIM) - IF (NW.GE.0) GO TO 100 - NZ = NN - DO 90 I=1,NN - CY(I) = CZERO - 90 CONTINUE - RETURN - 100 CONTINUE - IF (NW.GT.0) GO TO 130 - CALL CWRSK(Z, FNU, KODE, NN, CY, NW, CW, TOL, ELIM, ALIM) - IF (NW.LT.0) GO TO 130 - GO TO 120 - 110 CONTINUE -C----------------------------------------------------------------------- -C INCREMENT FNU+NN-1 UP TO FNUL, COMPUTE AND RECUR BACKWARD -C----------------------------------------------------------------------- - NUI = FNUL-DFNU + 1 - NUI = MAX(NUI,0) - CALL CBUNI(Z, FNU, KODE, NN, CY, NW, NUI, NLAST, FNUL, TOL, ELIM, - * ALIM) - IF (NW.LT.0) GO TO 130 - NZ = NZ + NW - IF (NLAST.EQ.0) GO TO 120 - NN = NLAST - GO TO 60 - 120 CONTINUE - RETURN - 130 CONTINUE - NZ = -1 - IF(NW.EQ.(-2)) NZ=-2 - RETURN - END diff --git a/schainf/acf2/cbknu.f b/schainf/acf2/cbknu.f deleted file mode 100644 index 03ff0a1..0000000 --- a/schainf/acf2/cbknu.f +++ /dev/null @@ -1,466 +0,0 @@ -*DECK CBKNU - SUBROUTINE CBKNU (Z, FNU, KODE, N, Y, NZ, TOL, ELIM, ALIM) -C***BEGIN PROLOGUE CBKNU -C***SUBSIDIARY -C***PURPOSE Subsidiary to CAIRY, CBESH, CBESI and CBESK -C***LIBRARY SLATEC -C***TYPE ALL (CBKNU-A, ZBKNU-A) -C***AUTHOR Amos, D. E., (SNL) -C***DESCRIPTION -C -C CBKNU COMPUTES THE K BESSEL FUNCTION IN THE RIGHT HALF Z PLANE -C -C***SEE ALSO CAIRY, CBESH, CBESI, CBESK -C***ROUTINES CALLED CKSCL, CSHCH, CUCHK, GAMLN, I1MACH, R1MACH -C***REVISION HISTORY (YYMMDD) -C 830501 DATE WRITTEN -C 910415 Prologue converted to Version 4.0 format. (BAB) -C***END PROLOGUE CBKNU -C - COMPLEX CCH, CK, COEF, CONE, CRSC, CS, CSCL, CSH, CSR, CSS, CTWO, - * CZ, CZERO, F, FMU, P, PT, P1, P2, Q, RZ, SMU, ST, S1, S2, Y, Z, - * ZD, CELM, CY - REAL AA, AK, ALIM, ASCLE, A1, A2, BB, BK, BRY, CAZ, CC, DNU, - * DNU2, ELIM, ETEST, FC, FHS, FK, FKS, FNU, FPI, G1, G2, HPI, PI, - * P2I, P2M, P2R, RK, RTHPI, R1, S, SPI, TM, TOL, TTH, T1, T2, XX, - * YY, GAMLN, R1MACH, HELIM, ELM, XD, YD, ALAS, AS - INTEGER I, IDUM, IFLAG, INU, K, KFLAG, KK, KMAX, KODE, KODED, N, - * NZ, I1MACH, NW, J, IC, INUB - DIMENSION BRY(3), CC(8), CSS(3), CSR(3), Y(N), CY(2) -C - DATA KMAX / 30 / - DATA R1 / 2.0E0 / - DATA CZERO,CONE,CTWO /(0.0E0,0.0E0),(1.0E0,0.0E0),(2.0E0,0.0E0)/ -C - DATA PI, RTHPI, SPI ,HPI, FPI, TTH / - 1 3.14159265358979324E0, 1.25331413731550025E0, - 2 1.90985931710274403E0, 1.57079632679489662E0, - 3 1.89769999331517738E0, 6.66666666666666666E-01/ -C - DATA CC(1), CC(2), CC(3), CC(4), CC(5), CC(6), CC(7), CC(8)/ - 1 5.77215664901532861E-01, -4.20026350340952355E-02, - 2 -4.21977345555443367E-02, 7.21894324666309954E-03, - 3 -2.15241674114950973E-04, -2.01348547807882387E-05, - 4 1.13302723198169588E-06, 6.11609510448141582E-09/ -C -C***FIRST EXECUTABLE STATEMENT CBKNU - XX = REAL(Z) - YY = AIMAG(Z) - CAZ = ABS(Z) - CSCL = CMPLX(1.0E0/TOL,0.0E0) - CRSC = CMPLX(TOL,0.0E0) - CSS(1) = CSCL - CSS(2) = CONE - CSS(3) = CRSC - CSR(1) = CRSC - CSR(2) = CONE - CSR(3) = CSCL - BRY(1) = 1.0E+3*R1MACH(1)/TOL - BRY(2) = 1.0E0/BRY(1) - BRY(3) = R1MACH(2) - NZ = 0 - IFLAG = 0 - KODED = KODE - RZ = CTWO/Z - INU = FNU+0.5E0 - DNU = FNU - INU - IF (ABS(DNU).EQ.0.5E0) GO TO 110 - DNU2 = 0.0E0 - IF (ABS(DNU).GT.TOL) DNU2 = DNU*DNU - IF (CAZ.GT.R1) GO TO 110 -C----------------------------------------------------------------------- -C SERIES FOR ABS(Z).LE.R1 -C----------------------------------------------------------------------- - FC = 1.0E0 - SMU = CLOG(RZ) - FMU = SMU*CMPLX(DNU,0.0E0) - CALL CSHCH(FMU, CSH, CCH) - IF (DNU.EQ.0.0E0) GO TO 10 - FC = DNU*PI - FC = FC/SIN(FC) - SMU = CSH*CMPLX(1.0E0/DNU,0.0E0) - 10 CONTINUE - A2 = 1.0E0 + DNU -C----------------------------------------------------------------------- -C GAM(1-Z)*GAM(1+Z)=PI*Z/SIN(PI*Z), T1=1/GAM(1-DNU), T2=1/GAM(1+DNU) -C----------------------------------------------------------------------- - T2 = EXP(-GAMLN(A2,IDUM)) - T1 = 1.0E0/(T2*FC) - IF (ABS(DNU).GT.0.1E0) GO TO 40 -C----------------------------------------------------------------------- -C SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU) -C----------------------------------------------------------------------- - AK = 1.0E0 - S = CC(1) - DO 20 K=2,8 - AK = AK*DNU2 - TM = CC(K)*AK - S = S + TM - IF (ABS(TM).LT.TOL) GO TO 30 - 20 CONTINUE - 30 G1 = -S - GO TO 50 - 40 CONTINUE - G1 = (T1-T2)/(DNU+DNU) - 50 CONTINUE - G2 = 0.5E0*(T1+T2)*FC - G1 = G1*FC - F = CMPLX(G1,0.0E0)*CCH + SMU*CMPLX(G2,0.0E0) - PT = CEXP(FMU) - P = CMPLX(0.5E0/T2,0.0E0)*PT - Q = CMPLX(0.5E0/T1,0.0E0)/PT - S1 = F - S2 = P - AK = 1.0E0 - A1 = 1.0E0 - CK = CONE - BK = 1.0E0 - DNU2 - IF (INU.GT.0 .OR. N.GT.1) GO TO 80 -C----------------------------------------------------------------------- -C GENERATE K(FNU,Z), 0.0D0 .LE. FNU .LT. 0.5D0 AND N=1 -C----------------------------------------------------------------------- - IF (CAZ.LT.TOL) GO TO 70 - CZ = Z*Z*CMPLX(0.25E0,0.0E0) - T1 = 0.25E0*CAZ*CAZ - 60 CONTINUE - F = (F*CMPLX(AK,0.0E0)+P+Q)*CMPLX(1.0E0/BK,0.0E0) - P = P*CMPLX(1.0E0/(AK-DNU),0.0E0) - Q = Q*CMPLX(1.0E0/(AK+DNU),0.0E0) - RK = 1.0E0/AK - CK = CK*CZ*CMPLX(RK,0.0) - S1 = S1 + CK*F - A1 = A1*T1*RK - BK = BK + AK + AK + 1.0E0 - AK = AK + 1.0E0 - IF (A1.GT.TOL) GO TO 60 - 70 CONTINUE - Y(1) = S1 - IF (KODED.EQ.1) RETURN - Y(1) = S1*CEXP(Z) - RETURN -C----------------------------------------------------------------------- -C GENERATE K(DNU,Z) AND K(DNU+1,Z) FOR FORWARD RECURRENCE -C----------------------------------------------------------------------- - 80 CONTINUE - IF (CAZ.LT.TOL) GO TO 100 - CZ = Z*Z*CMPLX(0.25E0,0.0E0) - T1 = 0.25E0*CAZ*CAZ - 90 CONTINUE - F = (F*CMPLX(AK,0.0E0)+P+Q)*CMPLX(1.0E0/BK,0.0E0) - P = P*CMPLX(1.0E0/(AK-DNU),0.0E0) - Q = Q*CMPLX(1.0E0/(AK+DNU),0.0E0) - RK = 1.0E0/AK - CK = CK*CZ*CMPLX(RK,0.0E0) - S1 = S1 + CK*F - S2 = S2 + CK*(P-F*CMPLX(AK,0.0E0)) - A1 = A1*T1*RK - BK = BK + AK + AK + 1.0E0 - AK = AK + 1.0E0 - IF (A1.GT.TOL) GO TO 90 - 100 CONTINUE - KFLAG = 2 - BK = REAL(SMU) - A1 = FNU + 1.0E0 - AK = A1*ABS(BK) - IF (AK.GT.ALIM) KFLAG = 3 - P2 = S2*CSS(KFLAG) - S2 = P2*RZ - S1 = S1*CSS(KFLAG) - IF (KODED.EQ.1) GO TO 210 - F = CEXP(Z) - S1 = S1*F - S2 = S2*F - GO TO 210 -C----------------------------------------------------------------------- -C IFLAG=0 MEANS NO UNDERFLOW OCCURRED -C IFLAG=1 MEANS AN UNDERFLOW OCCURRED- COMPUTATION PROCEEDS WITH -C KODED=2 AND A TEST FOR ON SCALE VALUES IS MADE DURING FORWARD -C RECURSION -C----------------------------------------------------------------------- - 110 CONTINUE - COEF = CMPLX(RTHPI,0.0E0)/CSQRT(Z) - KFLAG = 2 - IF (KODED.EQ.2) GO TO 120 - IF (XX.GT.ALIM) GO TO 290 -C BLANK LINE - A1 = EXP(-XX)*REAL(CSS(KFLAG)) - PT = CMPLX(A1,0.0E0)*CMPLX(COS(YY),-SIN(YY)) - COEF = COEF*PT - 120 CONTINUE - IF (ABS(DNU).EQ.0.5E0) GO TO 300 -C----------------------------------------------------------------------- -C MILLER ALGORITHM FOR ABS(Z).GT.R1 -C----------------------------------------------------------------------- - AK = COS(PI*DNU) - AK = ABS(AK) - IF (AK.EQ.0.0E0) GO TO 300 - FHS = ABS(0.25E0-DNU2) - IF (FHS.EQ.0.0E0) GO TO 300 -C----------------------------------------------------------------------- -C COMPUTE R2=F(E). IF ABS(Z).GE.R2, USE FORWARD RECURRENCE TO -C DETERMINE THE BACKWARD INDEX K. R2=F(E) IS A STRAIGHT LINE ON -C 12.LE.E.LE.60. E IS COMPUTED FROM 2**(-E)=B**(1-I1MACH(11))= -C TOL WHERE B IS THE BASE OF THE ARITHMETIC. -C----------------------------------------------------------------------- - T1 = (I1MACH(11)-1)*R1MACH(5)*3.321928094E0 - T1 = MAX(T1,12.0E0) - T1 = MIN(T1,60.0E0) - T2 = TTH*T1 - 6.0E0 - IF (XX.NE.0.0E0) GO TO 130 - T1 = HPI - GO TO 140 - 130 CONTINUE - T1 = ATAN(YY/XX) - T1 = ABS(T1) - 140 CONTINUE - IF (T2.GT.CAZ) GO TO 170 -C----------------------------------------------------------------------- -C FORWARD RECURRENCE LOOP WHEN ABS(Z).GE.R2 -C----------------------------------------------------------------------- - ETEST = AK/(PI*CAZ*TOL) - FK = 1.0E0 - IF (ETEST.LT.1.0E0) GO TO 180 - FKS = 2.0E0 - RK = CAZ + CAZ + 2.0E0 - A1 = 0.0E0 - A2 = 1.0E0 - DO 150 I=1,KMAX - AK = FHS/FKS - BK = RK/(FK+1.0E0) - TM = A2 - A2 = BK*A2 - AK*A1 - A1 = TM - RK = RK + 2.0E0 - FKS = FKS + FK + FK + 2.0E0 - FHS = FHS + FK + FK - FK = FK + 1.0E0 - TM = ABS(A2)*FK - IF (ETEST.LT.TM) GO TO 160 - 150 CONTINUE - GO TO 310 - 160 CONTINUE - FK = FK + SPI*T1*SQRT(T2/CAZ) - FHS = ABS(0.25E0-DNU2) - GO TO 180 - 170 CONTINUE -C----------------------------------------------------------------------- -C COMPUTE BACKWARD INDEX K FOR ABS(Z).LT.R2 -C----------------------------------------------------------------------- - A2 = SQRT(CAZ) - AK = FPI*AK/(TOL*SQRT(A2)) - AA = 3.0E0*T1/(1.0E0+CAZ) - BB = 14.7E0*T1/(28.0E0+CAZ) - AK = (ALOG(AK)+CAZ*COS(AA)/(1.0E0+0.008E0*CAZ))/COS(BB) - FK = 0.12125E0*AK*AK/CAZ + 1.5E0 - 180 CONTINUE - K = FK -C----------------------------------------------------------------------- -C BACKWARD RECURRENCE LOOP FOR MILLER ALGORITHM -C----------------------------------------------------------------------- - FK = K - FKS = FK*FK - P1 = CZERO - P2 = CMPLX(TOL,0.0E0) - CS = P2 - DO 190 I=1,K - A1 = FKS - FK - A2 = (FKS+FK)/(A1+FHS) - RK = 2.0E0/(FK+1.0E0) - T1 = (FK+XX)*RK - T2 = YY*RK - PT = P2 - P2 = (P2*CMPLX(T1,T2)-P1)*CMPLX(A2,0.0E0) - P1 = PT - CS = CS + P2 - FKS = A1 - FK + 1.0E0 - FK = FK - 1.0E0 - 190 CONTINUE -C----------------------------------------------------------------------- -C COMPUTE (P2/CS)=(P2/ABS(CS))*(CONJG(CS)/ABS(CS)) FOR BETTER -C SCALING -C----------------------------------------------------------------------- - TM = ABS(CS) - PT = CMPLX(1.0E0/TM,0.0E0) - S1 = PT*P2 - CS = CONJG(CS)*PT - S1 = COEF*S1*CS - IF (INU.GT.0 .OR. N.GT.1) GO TO 200 - ZD = Z - IF(IFLAG.EQ.1) GO TO 270 - GO TO 240 - 200 CONTINUE -C----------------------------------------------------------------------- -C COMPUTE P1/P2=(P1/ABS(P2)*CONJG(P2)/ABS(P2) FOR SCALING -C----------------------------------------------------------------------- - TM = ABS(P2) - PT = CMPLX(1.0E0/TM,0.0E0) - P1 = PT*P1 - P2 = CONJG(P2)*PT - PT = P1*P2 - S2 = S1*(CONE+(CMPLX(DNU+0.5E0,0.0E0)-PT)/Z) -C----------------------------------------------------------------------- -C FORWARD RECURSION ON THE THREE TERM RECURSION RELATION WITH -C SCALING NEAR EXPONENT EXTREMES ON KFLAG=1 OR KFLAG=3 -C----------------------------------------------------------------------- - 210 CONTINUE - CK = CMPLX(DNU+1.0E0,0.0E0)*RZ - IF (N.EQ.1) INU = INU - 1 - IF (INU.GT.0) GO TO 220 - IF (N.EQ.1) S1=S2 - ZD = Z - IF(IFLAG.EQ.1) GO TO 270 - GO TO 240 - 220 CONTINUE - INUB = 1 - IF (IFLAG.EQ.1) GO TO 261 - 225 CONTINUE - P1 = CSR(KFLAG) - ASCLE = BRY(KFLAG) - DO 230 I=INUB,INU - ST = S2 - S2 = CK*S2 + S1 - S1 = ST - CK = CK + RZ - IF (KFLAG.GE.3) GO TO 230 - P2 = S2*P1 - P2R = REAL(P2) - P2I = AIMAG(P2) - P2R = ABS(P2R) - P2I = ABS(P2I) - P2M = MAX(P2R,P2I) - IF (P2M.LE.ASCLE) GO TO 230 - KFLAG = KFLAG + 1 - ASCLE = BRY(KFLAG) - S1 = S1*P1 - S2 = P2 - S1 = S1*CSS(KFLAG) - S2 = S2*CSS(KFLAG) - P1 = CSR(KFLAG) - 230 CONTINUE - IF (N.EQ.1) S1 = S2 - 240 CONTINUE - Y(1) = S1*CSR(KFLAG) - IF (N.EQ.1) RETURN - Y(2) = S2*CSR(KFLAG) - IF (N.EQ.2) RETURN - KK = 2 - 250 CONTINUE - KK = KK + 1 - IF (KK.GT.N) RETURN - P1 = CSR(KFLAG) - ASCLE = BRY(KFLAG) - DO 260 I=KK,N - P2 = S2 - S2 = CK*S2 + S1 - S1 = P2 - CK = CK + RZ - P2 = S2*P1 - Y(I) = P2 - IF (KFLAG.GE.3) GO TO 260 - P2R = REAL(P2) - P2I = AIMAG(P2) - P2R = ABS(P2R) - P2I = ABS(P2I) - P2M = MAX(P2R,P2I) - IF (P2M.LE.ASCLE) GO TO 260 - KFLAG = KFLAG + 1 - ASCLE = BRY(KFLAG) - S1 = S1*P1 - S2 = P2 - S1 = S1*CSS(KFLAG) - S2 = S2*CSS(KFLAG) - P1 = CSR(KFLAG) - 260 CONTINUE - RETURN -C----------------------------------------------------------------------- -C IFLAG=1 CASES, FORWARD RECURRENCE ON SCALED VALUES ON UNDERFLOW -C----------------------------------------------------------------------- - 261 CONTINUE - HELIM = 0.5E0*ELIM - ELM = EXP(-ELIM) - CELM = CMPLX(ELM,0.0) - ASCLE = BRY(1) - ZD = Z - XD = XX - YD = YY - IC = -1 - J = 2 - DO 262 I=1,INU - ST = S2 - S2 = CK*S2+S1 - S1 = ST - CK = CK+RZ - AS = ABS(S2) - ALAS = ALOG(AS) - P2R = -XD+ALAS - IF(P2R.LT.(-ELIM)) GO TO 263 - P2 = -ZD+CLOG(S2) - P2R = REAL(P2) - P2I = AIMAG(P2) - P2M = EXP(P2R)/TOL - P1 = CMPLX(P2M,0.0E0)*CMPLX(COS(P2I),SIN(P2I)) - CALL CUCHK(P1,NW,ASCLE,TOL) - IF(NW.NE.0) GO TO 263 - J=3-J - CY(J) = P1 - IF(IC.EQ.(I-1)) GO TO 264 - IC = I - GO TO 262 - 263 CONTINUE - IF(ALAS.LT.HELIM) GO TO 262 - XD = XD-ELIM - S1 = S1*CELM - S2 = S2*CELM - ZD = CMPLX(XD,YD) - 262 CONTINUE - IF(N.EQ.1) S1 = S2 - GO TO 270 - 264 CONTINUE - KFLAG = 1 - INUB = I+1 - S2 = CY(J) - J = 3 - J - S1 = CY(J) - IF(INUB.LE.INU) GO TO 225 - IF(N.EQ.1) S1 = S2 - GO TO 240 - 270 CONTINUE - Y(1) = S1 - IF (N.EQ.1) GO TO 280 - Y(2) = S2 - 280 CONTINUE - ASCLE = BRY(1) - CALL CKSCL(ZD, FNU, N, Y, NZ, RZ, ASCLE, TOL, ELIM) - INU = N - NZ - IF (INU.LE.0) RETURN - KK = NZ + 1 - S1 = Y(KK) - Y(KK) = S1*CSR(1) - IF (INU.EQ.1) RETURN - KK = NZ + 2 - S2 = Y(KK) - Y(KK) = S2*CSR(1) - IF (INU.EQ.2) RETURN - T2 = FNU + (KK-1) - CK = CMPLX(T2,0.0E0)*RZ - KFLAG = 1 - GO TO 250 - 290 CONTINUE -C----------------------------------------------------------------------- -C SCALE BY EXP(Z), IFLAG = 1 CASES -C----------------------------------------------------------------------- - KODED = 2 - IFLAG = 1 - KFLAG = 2 - GO TO 120 -C----------------------------------------------------------------------- -C FNU=HALF ODD INTEGER CASE, DNU=-0.5 -C----------------------------------------------------------------------- - 300 CONTINUE - S1 = COEF - S2 = COEF - GO TO 210 - 310 CONTINUE - NZ=-2 - RETURN - END diff --git a/schainf/acf2/cbuni.f b/schainf/acf2/cbuni.f deleted file mode 100644 index 629851c..0000000 --- a/schainf/acf2/cbuni.f +++ /dev/null @@ -1,169 +0,0 @@ -*DECK CBUNI - SUBROUTINE CBUNI (Z, FNU, KODE, N, Y, NZ, NUI, NLAST, FNUL, TOL, - + ELIM, ALIM) -C***BEGIN PROLOGUE CBUNI -C***SUBSIDIARY -C***PURPOSE Subsidiary to CBESI and CBESK -C***LIBRARY SLATEC -C***TYPE ALL (CBUNI-A, ZBUNI-A) -C***AUTHOR Amos, D. E., (SNL) -C***DESCRIPTION -C -C CBUNI COMPUTES THE I BESSEL FUNCTION FOR LARGE ABS(Z).GT. -C FNUL AND FNU+N-1.LT.FNUL. THE ORDER IS INCREASED FROM -C FNU+N-1 GREATER THAN FNUL BY ADDING NUI AND COMPUTING -C ACCORDING TO THE UNIFORM ASYMPTOTIC EXPANSION FOR I(FNU,Z) -C ON IFORM=1 AND THE EXPANSION FOR J(FNU,Z) ON IFORM=2 -C -C***SEE ALSO CBESI, CBESK -C***ROUTINES CALLED CUNI1, CUNI2, R1MACH -C***REVISION HISTORY (YYMMDD) -C 830501 DATE WRITTEN -C 910415 Prologue converted to Version 4.0 format. (BAB) -C***END PROLOGUE CBUNI - COMPLEX CSCL, CSCR, CY, RZ, ST, S1, S2, Y, Z - REAL ALIM, AX, AY, DFNU, ELIM, FNU, FNUI, FNUL, GNU, TOL, XX, YY, - * ASCLE, BRY, STR, STI, STM, R1MACH - INTEGER I, IFLAG, IFORM, K, KODE, N, NL, NLAST, NUI, NW, NZ - DIMENSION Y(N), CY(2), BRY(3) -C***FIRST EXECUTABLE STATEMENT CBUNI - NZ = 0 - XX = REAL(Z) - YY = AIMAG(Z) - AX = ABS(XX)*1.7321E0 - AY = ABS(YY) - IFORM = 1 - IF (AY.GT.AX) IFORM = 2 - IF (NUI.EQ.0) GO TO 60 - FNUI = NUI - DFNU = FNU + (N-1) - GNU = DFNU + FNUI - IF (IFORM.EQ.2) GO TO 10 -C----------------------------------------------------------------------- -C ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN -C -PI/3.LE.ARG(Z).LE.PI/3 -C----------------------------------------------------------------------- - CALL CUNI1(Z, GNU, KODE, 2, CY, NW, NLAST, FNUL, TOL, ELIM, ALIM) - GO TO 20 - 10 CONTINUE -C----------------------------------------------------------------------- -C ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU -C APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I -C AND HPI=PI/2 -C----------------------------------------------------------------------- - CALL CUNI2(Z, GNU, KODE, 2, CY, NW, NLAST, FNUL, TOL, ELIM, ALIM) - 20 CONTINUE - IF (NW.LT.0) GO TO 50 - IF (NW.NE.0) GO TO 90 - AY = ABS(CY(1)) -C---------------------------------------------------------------------- -C SCALE BACKWARD RECURRENCE, BRY(3) IS DEFINED BUT NEVER USED -C---------------------------------------------------------------------- - BRY(1) = 1.0E+3*R1MACH(1)/TOL - BRY(2) = 1.0E0/BRY(1) - BRY(3) = BRY(2) - IFLAG = 2 - ASCLE = BRY(2) - AX = 1.0E0 - CSCL = CMPLX(AX,0.0E0) - IF (AY.GT.BRY(1)) GO TO 21 - IFLAG = 1 - ASCLE = BRY(1) - AX = 1.0E0/TOL - CSCL = CMPLX(AX,0.0E0) - GO TO 25 - 21 CONTINUE - IF (AY.LT.BRY(2)) GO TO 25 - IFLAG = 3 - ASCLE = BRY(3) - AX = TOL - CSCL = CMPLX(AX,0.0E0) - 25 CONTINUE - AY = 1.0E0/AX - CSCR = CMPLX(AY,0.0E0) - S1 = CY(2)*CSCL - S2 = CY(1)*CSCL - RZ = CMPLX(2.0E0,0.0E0)/Z - DO 30 I=1,NUI - ST = S2 - S2 = CMPLX(DFNU+FNUI,0.0E0)*RZ*S2 + S1 - S1 = ST - FNUI = FNUI - 1.0E0 - IF (IFLAG.GE.3) GO TO 30 - ST = S2*CSCR - STR = REAL(ST) - STI = AIMAG(ST) - STR = ABS(STR) - STI = ABS(STI) - STM = MAX(STR,STI) - IF (STM.LE.ASCLE) GO TO 30 - IFLAG = IFLAG+1 - ASCLE = BRY(IFLAG) - S1 = S1*CSCR - S2 = ST - AX = AX*TOL - AY = 1.0E0/AX - CSCL = CMPLX(AX,0.0E0) - CSCR = CMPLX(AY,0.0E0) - S1 = S1*CSCL - S2 = S2*CSCL - 30 CONTINUE - Y(N) = S2*CSCR - IF (N.EQ.1) RETURN - NL = N - 1 - FNUI = NL - K = NL - DO 40 I=1,NL - ST = S2 - S2 = CMPLX(FNU+FNUI,0.0E0)*RZ*S2 + S1 - S1 = ST - ST = S2*CSCR - Y(K) = ST - FNUI = FNUI - 1.0E0 - K = K - 1 - IF (IFLAG.GE.3) GO TO 40 - STR = REAL(ST) - STI = AIMAG(ST) - STR = ABS(STR) - STI = ABS(STI) - STM = MAX(STR,STI) - IF (STM.LE.ASCLE) GO TO 40 - IFLAG = IFLAG+1 - ASCLE = BRY(IFLAG) - S1 = S1*CSCR - S2 = ST - AX = AX*TOL - AY = 1.0E0/AX - CSCL = CMPLX(AX,0.0E0) - CSCR = CMPLX(AY,0.0E0) - S1 = S1*CSCL - S2 = S2*CSCL - 40 CONTINUE - RETURN - 50 CONTINUE - NZ = -1 - IF(NW.EQ.(-2)) NZ=-2 - RETURN - 60 CONTINUE - IF (IFORM.EQ.2) GO TO 70 -C----------------------------------------------------------------------- -C ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN -C -PI/3.LE.ARG(Z).LE.PI/3 -C----------------------------------------------------------------------- - CALL CUNI1(Z, FNU, KODE, N, Y, NW, NLAST, FNUL, TOL, ELIM, ALIM) - GO TO 80 - 70 CONTINUE -C----------------------------------------------------------------------- -C ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU -C APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I -C AND HPI=PI/2 -C----------------------------------------------------------------------- - CALL CUNI2(Z, FNU, KODE, N, Y, NW, NLAST, FNUL, TOL, ELIM, ALIM) - 80 CONTINUE - IF (NW.LT.0) GO TO 50 - NZ = NW - RETURN - 90 CONTINUE - NLAST = N - RETURN - END diff --git a/schainpy/model/data/jrodata.py b/schainpy/model/data/jrodata.py index 213dfcb..e9ec76c 100644 --- a/schainpy/model/data/jrodata.py +++ b/schainpy/model/data/jrodata.py @@ -548,6 +548,9 @@ class Spectra(JROData): if self.flagDecodeData: pwcode = numpy.sum(self.code[0]**2) + #pwcode = 64 + #print("pwcode: ", pwcode) + #exit(1) #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter diff --git a/schainpy/model/graphics/jroplot_spectra.py b/schainpy/model/graphics/jroplot_spectra.py index ffd9a0e..9b9ed90 100644 --- a/schainpy/model/graphics/jroplot_spectra.py +++ b/schainpy/model/graphics/jroplot_spectra.py @@ -42,6 +42,9 @@ class SpectraPlot(Plot): meta = {} spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) + #print("dataOut.normFactor: ", dataOut.normFactor) + #print("spc: ", dataOut.data_spc[0,0,0]) + #spc = 10*numpy.log10(dataOut.data_spc) #print("Spc: ",spc[0]) #exit(1) data['spc'] = spc @@ -716,14 +719,15 @@ class RTIPlot(Plot): self.x = self.data.times self.y = self.data.yrange self.z = self.data[self.CODE] - + #print("Inside RTI: ", self.z) self.z = numpy.ma.masked_invalid(self.z) if self.decimation is None: x, y, z = self.fill_gaps(self.x, self.y, self.z) else: x, y, z = self.fill_gaps(*self.decimate()) - + #print("self.z: ", self.z) + #exit(1) ''' if not isinstance(self.zmin, collections.abc.Sequence): if not self.zmin: diff --git a/schainpy/model/graphics/jroplot_voltage_lags.py b/schainpy/model/graphics/jroplot_voltage_lags.py index 4699136..7ebf8c9 100644 --- a/schainpy/model/graphics/jroplot_voltage_lags.py +++ b/schainpy/model/graphics/jroplot_voltage_lags.py @@ -673,9 +673,9 @@ class EDensityPlot(Plot): data['den_Faraday'] = dataOut.dphi[:dataOut.NSHTS] data['den_error'] = dataOut.sdp2[:dataOut.NSHTS] #data['err_Faraday'] = dataOut.sdn1[:dataOut.NSHTS] - print(numpy.shape(data['den_power'])) - print(numpy.shape(data['den_Faraday'])) - print(numpy.shape(data['den_error'])) + #print(numpy.shape(data['den_power'])) + #print(numpy.shape(data['den_Faraday'])) + #print(numpy.shape(data['den_error'])) data['NSHTS'] = dataOut.NSHTS diff --git a/schainpy/model/io/jroIO_madrigal.py b/schainpy/model/io/jroIO_madrigal.py index 5dd5bff..58d82ab 100644 --- a/schainpy/model/io/jroIO_madrigal.py +++ b/schainpy/model/io/jroIO_madrigal.py @@ -49,6 +49,7 @@ DEF_HEADER = { MNEMONICS = { 10: 'jro', 11: 'jbr', + 14: 'jmp', #Added by R. Flores 840: 'jul', 13: 'jas', 1000: 'pbr', diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py index 71475b5..665e61a 100644 --- a/schainpy/model/proc/jroproc_parameters.py +++ b/schainpy/model/proc/jroproc_parameters.py @@ -1152,7 +1152,7 @@ class Oblique_Gauss_Fit(Operation): #return A1f, B1f, C1f, A2f, B2f, C2f, K2f, Df, doppler return A1f, B1f, C1f, A2f, B2f, C2f, K2f, Df, doppler - def Double_Gauss_Double_Skew_fit_weight_bound_no_inputs(self,spc,freq,Nincoh): + def Double_Gauss_Double_Skew_fit_weight_bound_no_inputs(self,spc,freq,Nincoh,hei): from scipy.optimize import least_squares @@ -1176,7 +1176,7 @@ class Oblique_Gauss_Fit(Operation): #print(a1,b1,c1,a2,b2,c2,k2,d) #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]) #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]) - bounds=([0,-numpy.inf,0,-5,0,-400,0,0,0],[numpy.inf,-300,numpy.inf,5,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf]) + 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]) #print(bounds) #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): x0_value = numpy.array([spc_max,dop1_x0,30,-.1,spc_max/4, dop2_x0,150,1,1.0e7]) #x0_value = numpy.array([spc_max,-400.5,30,-.1,spc_max/4,-200.5,150,1,1.0e7]) #popt = least_squares(lsq_func,[a1,b1,c1,a2,b2,c2,k2,d],x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=1) + ''' + print("INSIDE 1") + print("x0_value: ", x0_value) + print("boundaries: ", bounds) + import matplotlib.pyplot as plt + plt.plot(freq,spc) + 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])) + plt.title(hei) + plt.show() + ''' popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0) # popt = least_squares(lsq_func,[a1,b1,c1,a2,b2,c2,k2,d],x_scale=params_scale,verbose=1) #print(popt) - + #########print("INSIDE 2") J = popt.jac try: @@ -1223,6 +1233,8 @@ class Oblique_Gauss_Fit(Operation): doppler2 = freq[numpy.argmax(aux2)] #print("error",error) #exit(1) + + return A1f, B1f, C1f, K1f, A2f, B2f, C2f, K2f, Df, doppler1, doppler2, error def Double_Gauss_fit_weight_bound_no_inputs(self,spc,freq,Nincoh): @@ -1468,6 +1480,87 @@ class Oblique_Gauss_Fit(Operation): return A1f, B1f, C1f, Df, error + def clean_outliers(self,param): + + threshold = 700 + + param = numpy.where(param < -threshold, numpy.nan, param) + param = numpy.where(param > +threshold, numpy.nan, param) + + return param + + def windowing_single(self,spc,x,A,B,C,D,nFFTPoints): + from scipy.optimize import curve_fit,fmin + + def R_gaussian(x, a, b, c): + N = int(numpy.shape(x)[0]) + val = a * numpy.exp(-((x)*c*2*2*numpy.pi)**2 / (2))* numpy.exp(1.j*b*x*4*numpy.pi) + return val + + def T(x,N): + T = 1-abs(x)/N + return T + + def R_T_spc_fun(x, a, b, c, d, nFFTPoints): + + N = int(numpy.shape(x)[0]) + + x_max = x[-1] + + x_pos = x[int(nFFTPoints/2):] + x_neg = x[:int(nFFTPoints/2)] + + R_T_neg_1 = R_gaussian(x, a, b, c)[:int(nFFTPoints/2)]*T(x_neg,-x[0]) + R_T_pos_1 = R_gaussian(x, a, b, c)[int(nFFTPoints/2):]*T(x_pos,x[-1]) + R_T_sum_1 = R_T_pos_1 + R_T_neg_1 + R_T_spc_1 = numpy.fft.fft(R_T_sum_1).real + R_T_spc_1 = numpy.fft.fftshift(R_T_spc_1) + max_val_1 = numpy.max(R_T_spc_1) + R_T_spc_1 = R_T_spc_1*a/max_val_1 + + R_T_d = d*numpy.fft.fftshift(signal.unit_impulse(N)) + R_T_d_neg = R_T_d[:int(nFFTPoints/2)]*T(x_neg,-x[0]) + R_T_d_pos = R_T_d[int(nFFTPoints/2):]*T(x_pos,x[-1]) + R_T_d_sum = R_T_d_pos + R_T_d_neg + R_T_spc_3 = numpy.fft.fft(R_T_d_sum).real + R_T_spc_3 = numpy.fft.fftshift(R_T_spc_3) + + R_T_final = R_T_spc_1 + R_T_spc_3 + + return R_T_final + + y = spc#gaussian(x, a, meanY, sigmaY) + a*0.1*numpy.random.normal(0, 1, size=len(x)) + + from scipy.stats import norm + mean,std=norm.fit(spc) + + # estimate starting values from the data + a = A + b = B + c = C#numpy.std(spc) + d = D + ''' + ippSeconds = 250*20*1.e-6/3 + + x_t = ippSeconds * (numpy.arange(1600) -1600 / 2.) + + x_t = numpy.linspace(x_t[0],x_t[-1],3200) + + x_freq = numpy.fft.fftfreq(1600,d=ippSeconds) + x_freq = numpy.fft.fftshift(x_freq) + ''' + # define a least squares function to optimize + def minfunc(params): + 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) + + # fit + + popt_full = fmin(minfunc,[a,b,c,d],full_output=True) + #print("nIter", popt_full[2]) + popt = popt_full[0] + + #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] + return popt[0], popt[1], popt[2], popt[3] def run(self, dataOut, mode = 0, Hmin1 = None, Hmax1 = None, Hmin2 = None, Hmax2 = None, Dop = 'Shift'): @@ -1565,10 +1658,17 @@ class Oblique_Gauss_Fit(Operation): pass else: - + #print("After: ", dataOut.data_snr[0]) + #######import matplotlib.pyplot as plt + #######plt.plot(dataOut.data_snr[0],dataOut.heightList,marker='*',linestyle='--') + #######plt.show() + #print("l1: ", dataOut.heightList[l1]) + #print("l2: ", dataOut.heightList[l2]) for hei in itertools.chain(l1, l2): #for hei in range(79,81): - if numpy.isnan(dataOut.data_snr[0,hei]): + #if numpy.isnan(dataOut.data_snr[0,hei]) or numpy.isnan(numpy.log10(dataOut.data_snr[0,hei])): + if numpy.isnan(dataOut.snl[0,hei]) or dataOut.snl[0,hei]<.0: + continue #Avoids the analysis when there is only noise try: @@ -1592,8 +1692,11 @@ class Oblique_Gauss_Fit(Operation): from scipy.signal import medfilt spcm = medfilt(spc,11) if x[numpy.argmax(spcm)] <= 0: - #print("EEJ") - 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) + #print("EEJ", dataOut.heightList[hei], hei) + #if hei != 70: + #continue + #else: + 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]) #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: # dataOut.Oblique_params[0,:,hei] *= numpy.NAN 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): dataOut.paramInterval = dataOut.nProfiles*dataOut.nCohInt*dataOut.ippSeconds dataOut.lat=-11.95 dataOut.lon=-76.87 - + ''' + dataOut.Oblique_params = numpy.where(dataOut.Oblique_params<-700, numpy.nan, dop_t1) + dataOut.Oblique_params = numpy.where(dataOut.Oblique_params<+700, numpy.nan, dop_t1) + Aquí debo exceptuar las amplitudes + ''' if mode == 9: #Double Skew Gaussian #dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,-2,:] #Pos[Max_value] #dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,1,:] #Shift @@ -1703,9 +1810,154 @@ class Oblique_Gauss_Fit(Operation): dataOut.Err_Dop_EEJ_T2 = dataOut.Oblique_param_errors[:,4,:] dataOut.Err_Spec_W_T2 = dataOut.Oblique_param_errors[:,5,:] + #print("Before: ", dataOut.Dop_EEJ_T2) + dataOut.Spec_W_T1 = self.clean_outliers(dataOut.Spec_W_T1) + dataOut.Spec_W_T2 = self.clean_outliers(dataOut.Spec_W_T2) + dataOut.Dop_EEJ_T1 = self.clean_outliers(dataOut.Dop_EEJ_T1) + dataOut.Dop_EEJ_T2 = self.clean_outliers(dataOut.Dop_EEJ_T2) + #print("After: ", dataOut.Dop_EEJ_T2) + dataOut.Err_Spec_W_T1 = self.clean_outliers(dataOut.Err_Spec_W_T1) + dataOut.Err_Spec_W_T2 = self.clean_outliers(dataOut.Err_Spec_W_T2) + dataOut.Err_Dop_EEJ_T1 = self.clean_outliers(dataOut.Err_Dop_EEJ_T1) + dataOut.Err_Dop_EEJ_T2 = self.clean_outliers(dataOut.Err_Dop_EEJ_T2) + #print("Before data_snr: ", dataOut.data_snr) + #dataOut.data_snr = numpy.where(numpy.isnan(dataOut.Dop_EEJ_T1), numpy.nan, dataOut.data_snr) + dataOut.snl = numpy.where(numpy.isnan(dataOut.Dop_EEJ_T1), numpy.nan, dataOut.snl) + + #print("After data_snr: ", dataOut.data_snr) dataOut.mode = mode dataOut.flagNoData = numpy.all(numpy.isnan(dataOut.Dop_EEJ_T1)) #Si todos los valores son NaN no se prosigue - #dataOut.flagNoData = False #Descomentar solo para ploteo sino mantener comentado + ###dataOut.flagNoData = False #Descomentar solo para ploteo sino mantener comentado (para guardado) + + return dataOut + +class Gaussian_Windowed(Operation): + ''' + Written by R. Flores + ''' + def __init__(self): + Operation.__init__(self) + + def windowing_single(self,spc,x,A,B,C,D,nFFTPoints): + from scipy.optimize import curve_fit,fmin + + def gaussian(x, a, b, c, d): + val = a * numpy.exp(-(x - b)**2 / (2*c**2)) + d + return val + + def R_gaussian(x, a, b, c): + N = int(numpy.shape(x)[0]) + val = a * numpy.exp(-((x)*c*2*2*numpy.pi)**2 / (2))* numpy.exp(1.j*b*x*4*numpy.pi) + return val + + def T(x,N): + T = 1-abs(x)/N + return T + + def R_T_spc_fun(x, a, b, c, d, nFFTPoints): + + N = int(numpy.shape(x)[0]) + + x_max = x[-1] + + x_pos = x[nFFTPoints:] + x_neg = x[:nFFTPoints] + #print([int(nFFTPoints/2)) + #print("x: ", x) + #print("x_neg: ", x_neg) + #print("x_pos: ", x_pos) + + + R_T_neg_1 = R_gaussian(x, a, b, c)[:nFFTPoints]*T(x_neg,-x[0]) + R_T_pos_1 = R_gaussian(x, a, b, c)[nFFTPoints:]*T(x_pos,x[-1]) + #print(T(x_pos,x[-1]),x_pos,x[-1]) + #print(R_T_neg_1.shape,R_T_pos_1.shape) + R_T_sum_1 = R_T_pos_1 + R_T_neg_1 + R_T_spc_1 = numpy.fft.fft(R_T_sum_1).real + R_T_spc_1 = numpy.fft.fftshift(R_T_spc_1) + max_val_1 = numpy.max(R_T_spc_1) + R_T_spc_1 = R_T_spc_1*a/max_val_1 + + R_T_d = d*numpy.fft.fftshift(signal.unit_impulse(N)) + R_T_d_neg = R_T_d[:nFFTPoints]*T(x_neg,-x[0]) + R_T_d_pos = R_T_d[nFFTPoints:]*T(x_pos,x[-1]) + R_T_d_sum = R_T_d_pos + R_T_d_neg + R_T_spc_3 = numpy.fft.fft(R_T_d_sum).real + R_T_spc_3 = numpy.fft.fftshift(R_T_spc_3) + + R_T_final = R_T_spc_1 + R_T_spc_3 + + return R_T_final + + y = spc#gaussian(x, a, meanY, sigmaY) + a*0.1*numpy.random.normal(0, 1, size=len(x)) + + from scipy.stats import norm + mean,std=norm.fit(spc) + + # estimate starting values from the data + a = A + b = B + c = C#numpy.std(spc) + d = D + #''' + #ippSeconds = 250*20*1.e-6/3 + + #x_t = ippSeconds * (numpy.arange(nFFTPoints) - nFFTPoints / 2.) + + #x_t = numpy.linspace(x_t[0],x_t[-1],3200) + #print("x_t: ", x_t) + #print("nFFTPoints: ", nFFTPoints) + x_vel = numpy.linspace(x[0],x[-1],int(2*nFFTPoints)) + #print("x_vel: ", x_vel) + #x_freq = numpy.fft.fftfreq(1600,d=ippSeconds) + #x_freq = numpy.fft.fftshift(x_freq) + #''' + # define a least squares function to optimize + def minfunc(params): + #print("y.shape: ", numpy.shape(y)) + return sum((y-R_T_spc_fun(x_vel,params[0],params[1],params[2],params[3],nFFTPoints))**2/1)#y**2) + + # fit + + popt_full = fmin(minfunc,[a,b,c,d], disp=False) + #print("nIter", popt_full[2]) + popt = popt_full#[0] + + fun = gaussian(x, popt[0], popt[1], popt[2], popt[3]) + + #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] + return fun, popt[0], popt[1], popt[2], popt[3] + + def run(self, dataOut): + + from scipy.signal import medfilt + import matplotlib.pyplot as plt + dataOut.moments = numpy.ones((dataOut.nChannels,4,dataOut.nHeights))*numpy.NAN + dataOut.VelRange = dataOut.getVelRange(0) + for nChannel in range(dataOut.nChannels): + for hei in range(dataOut.heightList.shape[0]): + #print("ipp: ", dataOut.ippSeconds) + spc = numpy.copy(dataOut.data_spc[nChannel,:,hei]) + + #print(VelRange) + #print(dataOut.getFreqRange(64)) + spcm = medfilt(spc,11) + spc_max = numpy.max(spcm) + dop1_x0 = dataOut.VelRange[numpy.argmax(spcm)] + D = numpy.min(spcm) + + fun, A, B, C, D = self.windowing_single(spc,dataOut.VelRange,spc_max,dop1_x0,abs(dop1_x0),D,dataOut.nFFTPoints) + dataOut.moments[nChannel,0,hei] = A + dataOut.moments[nChannel,1,hei] = B + dataOut.moments[nChannel,2,hei] = C + dataOut.moments[nChannel,3,hei] = D + ''' + plt.figure() + plt.plot(VelRange,spc,marker='*',linestyle='') + plt.plot(VelRange,fun) + plt.title(dataOut.heightList[hei]) + plt.show() + ''' return dataOut @@ -2547,7 +2799,7 @@ class SpectralFitting(Operation): #def __DiffCoherent(self,snrth, spectra, cspectra, nProf, heights,nChan, nHei, nPairs, channels, noise, crosspairs): def __DiffCoherent(self, spectra, cspectra, dataOut, noise, snrth, coh_th, hei_th): - import matplotlib.pyplot as plt + #import matplotlib.pyplot as plt nProf = dataOut.nProfiles heights = dataOut.heightList nHei = len(heights) @@ -2573,20 +2825,19 @@ class SpectralFitting(Operation): if coh_th == None : coh_th = numpy.array([0.75,0.65,0.15]) # 0.65 if hei_th == None : hei_th = numpy.array([60,300,650]) - for ic in range(2): + for ic in range(nPairs): pair = crosspairs[ic] #si el SNR es mayor que el SNR threshold los datos se toman coherentes s_n0 = power[pair[0],:]/noise[pair[0]] s_n1 = power[pair[1],:]/noise[pair[1]] - valid1 =(s_n0>=snr_th).nonzero() valid2 = (s_n1>=snr_th).nonzero() - #valid = valid2 + valid1 #numpy.concatenate((valid1,valid2), axis=None) + valid1 = numpy.array(valid1[0]) valid2 = numpy.array(valid2[0]) valid = valid1 for iv in range(len(valid2)): - #for ivv in range(len(valid1)) : + indv = numpy.array((valid1 == valid2[iv]).nonzero()) if len(indv[0]) == 0 : valid = numpy.concatenate((valid,valid2[iv]), axis=None) @@ -2594,17 +2845,16 @@ class SpectralFitting(Operation): my_coh_aver[pair[0],valid]=1 my_coh_aver[pair[1],valid]=1 # si la coherencia es mayor a la coherencia threshold los datos se toman - #print my_coh_aver[0,:] + 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))) - #print('coh',numpy.absolute(coh)) + for ih in range(len(hei_th)): hvalid = (heights>hei_th[ih]).nonzero() hvalid = hvalid[0] if len(hvalid)>0: valid = (numpy.absolute(coh[hvalid])>coh_th[ih]).nonzero() valid = valid[0] - #print('hvalid:',hvalid) - #print('valid', valid) + if len(valid)>0: my_coh_aver[pair[0],hvalid[valid]] =1 my_coh_aver[pair[1],hvalid[valid]] =1 @@ -2620,7 +2870,7 @@ class SpectralFitting(Operation): my_incoh_aver[pair[1],incoh_echoes] = 1 - for ic in range(2): + for ic in range(nPairs): pair = crosspairs[ic] valid1 =(my_coh_aver[pair[0],:]==1 ).nonzero() @@ -2628,29 +2878,25 @@ class SpectralFitting(Operation): valid1 = numpy.array(valid1[0]) valid2 = numpy.array(valid2[0]) valid = valid1 - #print valid1 , valid2 + for iv in range(len(valid2)): - #for ivv in range(len(valid1)) : + indv = numpy.array((valid1 == valid2[iv]).nonzero()) if len(indv[0]) == 0 : valid = numpy.concatenate((valid,valid2[iv]), axis=None) - #print valid - #valid = numpy.concatenate((valid1,valid2), axis=None) valid1 =(my_coh_aver[pair[0],:] !=1 ).nonzero() valid2 = (my_coh_aver[pair[1],:] !=1).nonzero() valid1 = numpy.array(valid1[0]) valid2 = numpy.array(valid2[0]) incoh_echoes = valid1 - #print valid1, valid2 - #incoh_echoes= numpy.concatenate((valid1,valid2), axis=None) + for iv in range(len(valid2)): - #for ivv in range(len(valid1)) : + indv = numpy.array((valid1 == valid2[iv]).nonzero()) if len(indv[0]) == 0 : incoh_echoes = numpy.concatenate(( incoh_echoes,valid2[iv]), axis=None) - #print incoh_echoes + if len(valid)>0: - #print pair coh_spectra[pair[0],:,valid] = spectra[pair[0],:,valid] coh_spectra[pair[1],:,valid] = spectra[pair[1],:,valid] coh_cspectra[ic,:,valid] = cspectra[ic,:,valid] @@ -2662,20 +2908,12 @@ class SpectralFitting(Operation): incoh_cspectra[ic,:,incoh_echoes] = cspectra[ic,:,incoh_echoes] incoh_aver[pair[0],incoh_echoes]=1 incoh_aver[pair[1],incoh_echoes]=1 - #plt.imshow(spectra[0,:,:],vmin=20000000) - #plt.show() - #my_incoh_aver = my_incoh_aver+1 - - #spec = my_incoh_spectra.copy() - #cspec = my_incoh_cspectra.copy() - #print('######################', spec) - #print(self.numpy) - #return spec, cspec,coh_aver + 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 def __CleanCoherent(self,snrth, spectra, cspectra, coh_aver,dataOut, noise,clean_coh_echoes,index): - import matplotlib.pyplot as plt + #import matplotlib.pyplot as plt nProf = dataOut.nProfiles heights = dataOut.heightList nHei = len(heights) @@ -2684,15 +2922,9 @@ class SpectralFitting(Operation): crosspairs = dataOut.groupList nPairs = len(crosspairs) - #data = dataOut.data_pre[0] absc = dataOut.abscissaList[:-1] - #noise = dataOut.noise - #nChannel = data.shape[0] data_param = numpy.zeros((nChan, 4, spectra.shape[2])) - - #plt.plot(absc) - #plt.show() clean_coh_spectra = spectra.copy() clean_coh_cspectra = cspectra.copy() clean_coh_aver = coh_aver.copy() @@ -2703,17 +2935,14 @@ class SpectralFitting(Operation): rtime0 = [6,18] # periodo sin ESF rtime1 = [10.5,13.5] # periodo con alta coherencia y alto ancho espectral (esperado): SOL. - time = index*5./60 + time = index*5./60 # en base a 5 min de proceso if clean_coh_echoes == 1 : for ind in range(nChan): data_param[ind,:,:] = self.__calculateMoments( spectra[ind,:,:] , absc , noise[ind] ) - #print data_param[:,3] + spwd = data_param[:,3] - #print spwd.shape + # 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 - #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] - #spwd=numpy.array([spwd1,spwd1,spwd1,spwd1]) - #print spwd.shape, heights.shape,coh_aver.shape # para obtener spwd for ic in range(nPairs): pair = crosspairs[ic] @@ -2748,29 +2977,20 @@ class SpectralFitting(Operation): self.i=0 self.isConfig = False - def setup(self,nChan,nProf,nHei,nBlocks): self.__dataReady = False self.bloques = numpy.zeros([2, nProf, nHei,nBlocks], dtype= complex) self.bloque0 = numpy.zeros([nChan, nProf, nHei, nBlocks]) - #def CleanRayleigh(self,dataOut,spectra,cspectra,out_spectra,out_cspectra,sat_spectra,sat_cspectra,crosspairs,heights, channels, nProf,nHei,nChan,nPairs,nIncohInt,nBlocks): def CleanRayleigh(self,dataOut,spectra,cspectra,save_drifts): - #import matplotlib.pyplot as plt - #for k in range(149): - # self.bloque0[:,:,:,k] = spectra[:,:,0:nHei] - # self.bloques[:,:,:,k] = cspectra[:,:,0:nHei] - #if self.i==nBlocks: - # self.i==0 - rfunc = cspectra.copy() #self.bloques + rfunc = cspectra.copy() n_funct = len(rfunc[0,:,0,0]) - val_spc = spectra*0.0 #self.bloque0*0.0 - val_cspc = cspectra*0.0 #self.bloques*0.0 - in_sat_spectra = spectra.copy() #self.bloque0 - in_sat_cspectra = cspectra.copy() #self.bloques + val_spc = spectra*0.0 + val_cspc = cspectra*0.0 + in_sat_spectra = spectra.copy() + in_sat_cspectra = cspectra.copy() - #print( rfunc.shape) min_hei = 200 nProf = dataOut.nProfiles heights = dataOut.heightList @@ -2781,20 +3001,19 @@ class SpectralFitting(Operation): nPairs = len(crosspairs) hval=(heights >= min_hei).nonzero() ih=hval[0] - #print numpy.absolute(rfunc[:,0,0,14]) for ih in range(hval[0][0],nHei): for ifreq in range(nProf): for ii in range(n_funct): func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) - #print numpy.amin(func2clean) + val = (numpy.isfinite(func2clean)==True).nonzero() if len(val)>0: min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40) if min_val <= -40 : min_val = -40 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200 if max_val >= 200 : max_val = 200 - #print min_val, max_val + step = 1 #Getting bins and the histogram x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step @@ -2809,18 +3028,6 @@ class SpectralFitting(Operation): except: mode = mean stdv = sigma -# if ih == 14 and ii == 0 and ifreq ==0 : -# print x_dist.shape, y_dist.shape -# print x_dist, y_dist -# print min_val, max_val, binstep -# print func2clean -# print mean,sigma -# mean1,std = norm.fit(y_dist) -# print mean1, std, gauss_fit -# print fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2]) - # 7.84616 53.9307 3.61863 - #stdv = 3.61863 # 2.99089 - #mode = 53.9307 #7.79008 #Removing echoes greater than mode + 3*stdv factor_stdv = 2.5 @@ -2831,29 +3038,17 @@ class SpectralFitting(Operation): cross_pairs = crosspairs[ii] #Getting coherent echoes which are removed. if len(novall[0]) > 0: - #val_spc[(0,1),novall[a],ih] = 1 - #val_spc[,(2,3),novall[a],ih] = 1 val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1 val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1 val_cspc[novall[0],ii,ifreq,ih] = 1 - #print("OUT NOVALL 1") #Removing coherent from ISR data -# if ih == 17 and ii == 0 and ifreq ==0 : -# print spectra[:,cross_pairs[0],ifreq,ih] spectra[noval,cross_pairs[0],ifreq,ih] = numpy.nan spectra[noval,cross_pairs[1],ifreq,ih] = numpy.nan cspectra[noval,ii,ifreq,ih] = numpy.nan -# if ih == 17 and ii == 0 and ifreq ==0 : -# print spectra[:,cross_pairs[0],ifreq,ih] -# print noval, len(noval[0]) -# print novall, len(novall[0]) -# print factor_stdv*stdv -# print func2clean-mode -# print val_spc[:,cross_pairs[0],ifreq,ih] -# print spectra[:,cross_pairs[0],ifreq,ih] +# #no sale es para savedrifts >2 - ''' channels = channels - cross_pairs = cross_pairs + ''' channels = dataOut.channelList + cross_pairs = dataOut.groupList #print("OUT NOVALL 2") vcross0 = (cross_pairs[0] == channels[ii]).nonzero() @@ -2874,6 +3069,7 @@ class SpectralFitting(Operation): self.bloques[vcross,ifreq,ih,noval] = numpy.nan ''' #Getting average of the spectra and cross-spectra from incoherent echoes. + out_spectra = numpy.zeros([nChan,nProf,nHei], dtype=float) #+numpy.nan out_cspectra = numpy.zeros([nPairs,nProf,nHei], dtype=complex) #+numpy.nan for ih in range(nHei): @@ -2881,22 +3077,15 @@ class SpectralFitting(Operation): for ich in range(nChan): tmp = spectra[:,ich,ifreq,ih] valid = (numpy.isfinite(tmp[:])==True).nonzero() -# if ich == 0 and ifreq == 0 and ih == 17 : -# print tmp -# print valid -# print len(valid[0]) - #print('TMP',tmp) if len(valid[0]) >0 : out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0]) - #for icr in range(nPairs): + for icr in range(nPairs): tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih]) valid = (numpy.isfinite(tmp)==True).nonzero() if len(valid[0]) > 0: out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0]) - # print('##########################################################') #Removing fake coherent echoes (at least 4 points around the point) - val_spectra = numpy.sum(val_spc,0) val_cspectra = numpy.sum(val_cspc,0) @@ -2932,13 +3121,6 @@ class SpectralFitting(Operation): tmp_sat_cspectra = cspectra.copy() tmp_sat_cspectra = tmp_sat_cspectra*numpy.nan -# fig = plt.figure(figsize=(6,5)) -# left, bottom, width, height = 0.1, 0.1, 0.8, 0.8 -# ax = fig.add_axes([left, bottom, width, height]) -# cp = ax.contour(10*numpy.log10(numpy.absolute(spectra[0,0,:,:]))) -# ax.clabel(cp, inline=True,fontsize=10) -# plt.show() - val = (val_spc > 0).nonzero() if len(val[0]) > 0: tmp_sat_spectra[val] = in_sat_spectra[val] @@ -2963,52 +3145,29 @@ class SpectralFitting(Operation): valid = (numpy.isfinite(tmp)).nonzero() if len(valid[0]) > 0: sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0]) - #self.__dataReady= True - #sat_spectra, sat_cspectra= sat_spectra, sat_cspectra - #if not self.__dataReady: - #return None, None + return out_spectra, out_cspectra,sat_spectra,sat_cspectra def REM_ISOLATED_POINTS(self,array,rth): -# import matplotlib.pyplot as plt if rth == None : rth = 4 - num_prof = len(array[0,:,0]) num_hei = len(array[0,0,:]) n2d = len(array[:,0,0]) for ii in range(n2d) : - #print ii,n2d tmp = array[ii,:,:] - #print tmp.shape, array[ii,101,:],array[ii,102,:] - -# fig = plt.figure(figsize=(6,5)) -# left, bottom, width, height = 0.1, 0.1, 0.8, 0.8 -# ax = fig.add_axes([left, bottom, width, height]) -# x = range(num_prof) -# y = range(num_hei) -# cp = ax.contour(y,x,tmp) -# ax.clabel(cp, inline=True,fontsize=10) -# plt.show() - - #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs) tmp = numpy.reshape(tmp,num_prof*num_hei) indxs1 = (numpy.isfinite(tmp)==True).nonzero() indxs2 = (tmp > 0).nonzero() - indxs1 = (indxs1[0]) indxs2 = indxs2[0] - #indxs1 = numpy.array(indxs1[0]) - #indxs2 = numpy.array(indxs2[0]) indxs = None - #print indxs1 , indxs2 + for iv in range(len(indxs2)): indv = numpy.array((indxs1 == indxs2[iv]).nonzero()) - #print len(indxs2), indv if len(indv[0]) > 0 : indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None) -# print indxs + indxs = indxs[1:] - #print indxs, len(indxs) if len(indxs) < 4 : array[ii,:,:] = 0. return @@ -3016,7 +3175,6 @@ class SpectralFitting(Operation): xpos = numpy.mod(indxs ,num_hei) ypos = (indxs / num_hei) sx = numpy.argsort(xpos) # Ordering respect to "x" (time) - #print sx xpos = xpos[sx] ypos = ypos[sx] @@ -3024,36 +3182,25 @@ class SpectralFitting(Operation): ic = 0 while True : r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2))) - #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh) - #plt.plot(r) - #plt.show() + no_coh1 = (numpy.isfinite(r)==True).nonzero() no_coh2 = (r <= rth).nonzero() - #print r, no_coh1, no_coh2 no_coh1 = numpy.array(no_coh1[0]) no_coh2 = numpy.array(no_coh2[0]) no_coh = None - #print valid1 , valid2 for iv in range(len(no_coh2)): indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero()) if len(indv[0]) > 0 : no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None) no_coh = no_coh[1:] - #print len(no_coh), no_coh if len(no_coh) < 4 : - #print xpos[ic], ypos[ic], ic -# plt.plot(r) -# plt.show() xpos[ic] = numpy.nan ypos[ic] = numpy.nan ic = ic + 1 if (ic == len(indxs)) : break - #print( xpos, ypos) - indxs = (numpy.isfinite(list(xpos))==True).nonzero() - #print indxs[0] if len(indxs[0]) < 4 : array[ii,:,:] = 0. return @@ -3068,27 +3215,11 @@ class SpectralFitting(Operation): tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))] array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei)) - #print array.shape - #tmp = numpy.reshape(tmp,(num_prof,num_hei)) - #print tmp.shape - -# fig = plt.figure(figsize=(6,5)) -# left, bottom, width, height = 0.1, 0.1, 0.8, 0.8 -# ax = fig.add_axes([left, bottom, width, height]) -# x = range(num_prof) -# y = range(num_hei) -# cp = ax.contour(y,x,array[ii,:,:]) -# ax.clabel(cp, inline=True,fontsize=10) -# plt.show() return array def moments(self,doppler,yarray,npoints): ytemp = yarray - #val = WHERE(ytemp GT 0,cval) - #if cval == 0 : val = range(npoints-1) val = (ytemp > 0).nonzero() val = val[0] - #print('hvalid:',hvalid) - #print('valid', valid) if len(val) == 0 : val = range(npoints-1) ynew = 0.5*(ytemp[val[0]]+ytemp[val[len(val)-1]]) @@ -3102,19 +3233,289 @@ class SpectralFitting(Operation): fmom = numpy.sum(doppler*ytemp)/numpy.sum(ytemp)+(index-(npoints/2-1))*numpy.abs(doppler[1]-doppler[0]) smom = numpy.sum(doppler*doppler*ytemp)/numpy.sum(ytemp) return [fmom,numpy.sqrt(smom)] + + def windowing_single_old(self,spc,x,A,B,C,D,nFFTPoints): + ''' + Written by R. Flores + ''' + from scipy.optimize import curve_fit,fmin + + def gaussian(x, a, b, c, d): + val = a * numpy.exp(-(x - b)**2 / (2*c**2)) + d + return val + + def R_gaussian(x, a, b, c): + N = int(numpy.shape(x)[0]) + val = a * numpy.exp(-((x)*c*2*2*numpy.pi)**2 / (2))* numpy.exp(1.j*b*x*4*numpy.pi) + return val + + def T(x,N): + T = 1-abs(x)/N + return T + + def R_T_spc_fun(x, a, b, c, d, nFFTPoints): + + N = int(numpy.shape(x)[0]) + + x_max = x[-1] + + x_pos = x[nFFTPoints:] + x_neg = x[:nFFTPoints] + + R_T_neg_1 = R_gaussian(x, a, b, c)[:nFFTPoints]*T(x_neg,-x[0]) + R_T_pos_1 = R_gaussian(x, a, b, c)[nFFTPoints:]*T(x_pos,x[-1]) + #print(T(x_pos,x[-1]),x_pos,x[-1]) + #print(R_T_neg_1.shape,R_T_pos_1.shape) + R_T_sum_1 = R_T_pos_1 + R_T_neg_1 + R_T_spc_1 = numpy.fft.fft(R_T_sum_1).real + R_T_spc_1 = numpy.fft.fftshift(R_T_spc_1) + max_val_1 = numpy.max(R_T_spc_1) + R_T_spc_1 = R_T_spc_1*a/max_val_1 + print("R_T_spc_1: ", R_T_spc_1) + + R_T_d = d*numpy.fft.fftshift(signal.unit_impulse(N)) + R_T_d_neg = R_T_d[:nFFTPoints]*T(x_neg,-x[0]) + R_T_d_pos = R_T_d[nFFTPoints:]*T(x_pos,x[-1]) + R_T_d_sum = R_T_d_pos + R_T_d_neg + R_T_spc_3 = numpy.fft.fft(R_T_d_sum).real + R_T_spc_3 = numpy.fft.fftshift(R_T_spc_3) + + R_T_final = R_T_spc_1# + R_T_spc_3 + + return R_T_final + + y = spc#gaussian(x, a, meanY, sigmaY) + a*0.1*numpy.random.normal(0, 1, size=len(x)) + + from scipy.stats import norm + mean,std=norm.fit(spc) + + # estimate starting values from the data + print("A: ", A) + a = A-D + b = B + c = C#numpy.std(spc) #C + d = D + #''' + #ippSeconds = 250*20*1.e-6/3 + + #x_t = ippSeconds * (numpy.arange(nFFTPoints) - nFFTPoints / 2.) + + #x_t = numpy.linspace(x_t[0],x_t[-1],3200) + #print("x_t: ", x_t) + #print("nFFTPoints: ", nFFTPoints) + x_vel = numpy.linspace(x[0],x[-1],int(2*nFFTPoints)) + #print("x_vel: ", x_vel) + #x_freq = numpy.fft.fftfreq(1600,d=ippSeconds) + #x_freq = numpy.fft.fftshift(x_freq) + #''' + # define a least squares function to optimize + import matplotlib.pyplot as plt + aui = R_T_spc_fun(x_vel,a,b,c,d,nFFTPoints) + print("aux_max: ", numpy.nanmax(aui)) + #print(dataOut.heightList[hei]) + plt.figure() + plt.plot(x,spc,marker='*',linestyle='--') + plt.plot(x,gaussian(x, a, b, c, d),color='b',marker='^',linestyle='') + plt.plot(x,aui,color='k') + #plt.title(dataOut.heightList[hei]) + plt.show() + + def minfunc(params): + #print("y.shape: ", numpy.shape(y)) + return sum((y-R_T_spc_fun(x_vel,params[0],params[1],params[2],params[3],nFFTPoints))**2/1)#y**2) + + # fit + + popt_full = fmin(minfunc,[a,b,c,d], disp=False) + #print("nIter", popt_full[2]) + popt = popt_full#[0] + + fun = gaussian(x, popt[0], popt[1], popt[2], popt[3]) + print("pop1[0]: ", popt[0]) + #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] + return fun, popt[0], popt[1], popt[2], popt[3] + + def windowing_single(self,spc,x,A,B,C,D,nFFTPoints): + ''' + Written by R. Flores + ''' + from scipy.optimize import curve_fit,fmin + + def gaussian(x, a, b, c, d): + val = a * numpy.exp(-(x - b)**2 / (2*c**2)) + d + return val + + def R_gaussian(x, a, b, c): + N = int(numpy.shape(x)[0]) + + val = (a*numpy.exp((-(1/2)*x*(x*c**2 + 2*1.j*b)))/numpy.sqrt(1/c**2)) + + return val + + def T(x,N): + T = 1-abs(x)/N + return T + + def R_T_spc_fun(x, a, id_dop, c, d, nFFTPoints): + + N = int(numpy.shape(x)[0]) + b = 0 + x_max = x[-1] + + x_pos = x[nFFTPoints:] + x_neg = x[:nFFTPoints] + R_T_neg_1 = R_gaussian(x, a, b, c)[:nFFTPoints]*T(x_neg,-x[0]) + R_T_pos_1 = R_gaussian(x, a, b, c)[nFFTPoints:]*T(x_pos,x[-1]) + + R_T_sum_1 = R_T_pos_1 + R_T_neg_1 + R_T_spc_1 = numpy.fft.fft(R_T_sum_1).real + R_T_spc_1 = numpy.fft.fftshift(R_T_spc_1) + max_val_1 = numpy.max(R_T_spc_1) + R_T_spc_1 = R_T_spc_1*a/max_val_1 + #raise NotImplementedError + R_T_d = d*numpy.fft.fftshift(signal.unit_impulse(N)) + R_T_d_neg = R_T_d[:nFFTPoints]*T(x_neg,-x[0]) + R_T_d_pos = R_T_d[nFFTPoints:]*T(x_pos,x[-1]) + R_T_d_sum = R_T_d_pos + R_T_d_neg + R_T_spc_3 = numpy.fft.fft(R_T_d_sum).real + R_T_spc_3 = numpy.fft.fftshift(R_T_spc_3) + + R_T_final = R_T_spc_1 + R_T_spc_3 + + id_dop = int(id_dop) + + R_T_final = numpy.roll(R_T_final,-id_dop) + + return R_T_final + + y = spc#gaussian(x, a, meanY, sigmaY) + a*0.1*numpy.random.normal(0, 1, size=len(x)) + + from scipy.stats import norm + mean,std=norm.fit(spc) + + # estimate starting values from the data + a = A-D + b = B + c = C#numpy.std(spc) #C + d = D + + id_dop = numpy.argmax(spc) + id_dop = int(spc.shape[0]/2 - id_dop) + + x_vel = numpy.linspace(x[0],x[-1],int(2*nFFTPoints)) + + # define a least squares function to optimize + + def minfunc(params): + #print("y.shape: ", numpy.shape(y)) + return sum((y-R_T_spc_fun(x_vel,params[0],params[1],params[2],params[3],nFFTPoints))**2/1)#y**2) + + # fit + popt_full = fmin(minfunc,[a,id_dop,c,d], disp=False) + popt = popt_full#[0] + + fun = gaussian(x, a, 0, popt[2], popt[3]) + fun = numpy.roll(fun,-int(popt[1])) + + return fun, popt[0], popt[1], popt[2], popt[3] + + def windowing_single_direct(self,spc_mod,x,A,B,C,D,nFFTPoints,timeInterval): + ''' + Written by R. Flores + ''' + from scipy.optimize import curve_fit,fmin + + def gaussian(x, a, b, c, d): + val = a * numpy.exp(-(x - b)**2 / (2*c**2)) + d + return val + + def R_gaussian(x, a, b, c, d): + N = int(numpy.shape(x)[0]) + 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 + + return 2*val/numpy.shape(val)[0] + + def T(x,N): + T = 1-abs(x)/N + return T + + def R_T_spc_fun(x, a, b, c, d, nFFTPoints, timeInterval): #"x" should be time + + #timeInterval = 2 + x_double = numpy.linspace(0,timeInterval,nFFTPoints) + x_double_m = numpy.flip(x_double) + x_double_aux = numpy.linspace(0,x_double[-2],nFFTPoints) + x_double_t = numpy.concatenate((x_double_m,x_double_aux)) + x_double_t /= max(x_double_t) + + + R_T_sum_1 = R_gaussian(x, a, b, c, d) + + R_T_sum_1_flip = numpy.copy(numpy.flip(R_T_sum_1)) + R_T_sum_1_flip[-1] = R_T_sum_1_flip[0] + R_T_sum_1_flip = numpy.roll(R_T_sum_1_flip,1) + + R_T_sum_1_flip.imag *= -1 + + R_T_sum_1_total = numpy.concatenate((R_T_sum_1,R_T_sum_1_flip)) + R_T_sum_1_total *= x_double_t #times trian_fun + + R_T_sum_1_total = R_T_sum_1_total[:nFFTPoints] + R_T_sum_1_total[nFFTPoints:] + + R_T_spc_1 = numpy.fft.fft(R_T_sum_1_total).real + R_T_spc_1 = numpy.fft.fftshift(R_T_spc_1) + + freq = numpy.fft.fftfreq(nFFTPoints, d=timeInterval/nFFTPoints) + + freq = numpy.fft.fftshift(freq) + + freq *= 6/2 #lambda/2 + + return R_T_spc_1 + + y = spc_mod + + #from scipy.stats import norm + + # estimate starting values from the data + + a = A-D + b = B + c = C + d = D + + # define a least squares function to optimize + import matplotlib.pyplot as plt + #ippSeconds = 2 + t_range = numpy.linspace(0,timeInterval,nFFTPoints) + #aui = R_T_spc_fun(t_range,a,b,c,d,nFFTPoints,timeInterval) + + def minfunc(params): + return sum((y-R_T_spc_fun(t_range,params[0],params[1],params[2],params[3],nFFTPoints,timeInterval))**2/1)#y**2) + + # fit + popt_full = fmin(minfunc,[a,b,c,d], disp=False) + popt = popt_full + + fun = R_T_spc_fun(t_range,popt[0],popt[1],popt[2],popt[3],nFFTPoints,timeInterval) + + return fun, popt[0], popt[1], popt[2], popt[3] + # ********************************************************************************************** index = 0 fint = 0 buffer = 0 buffer2 = 0 buffer3 = 0 - def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None): + 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): nChannels = dataOut.nChannels nHeights= dataOut.heightList.size nProf = dataOut.nProfiles + if numpy.any(taver): taver=int(taver) + else : taver = 15 tini=time.localtime(dataOut.utctime) - if (tini.tm_min % 5) == 0 and (tini.tm_sec < 5 and self.fint==0): -# print tini.tm_min + if (tini.tm_min % taver) == 0 and (tini.tm_sec < 5 and self.fint==0): + self.index = 0 jspc = self.buffer jcspc = self.buffer2 @@ -3124,14 +3525,14 @@ class SpectralFitting(Operation): self.buffer3 = dataOut.noise self.fint = 1 if numpy.any(jspc) : - jspc= numpy.reshape(jspc,(int(len(jspc)/4),nChannels,nProf,nHeights)) - jcspc= numpy.reshape(jcspc,(int(len(jcspc)/2),2,nProf,nHeights)) - jnoise= numpy.reshape(jnoise,(int(len(jnoise)/4),nChannels)) + jspc= numpy.reshape(jspc,(int(len(jspc)/nChannels),nChannels,nProf,nHeights)) + jcspc= numpy.reshape(jcspc,(int(len(jcspc)/int(nChannels/2)),int(nChannels/2),nProf,nHeights)) + jnoise= numpy.reshape(jnoise,(int(len(jnoise)/nChannels),nChannels)) else: dataOut.flagNoData = True return dataOut else : - if (tini.tm_min % 5) == 0 : self.fint = 1 + if (tini.tm_min % taver) == 0 : self.fint = 1 else : self.fint = 0 self.index += 1 if numpy.any(self.buffer): @@ -3146,7 +3547,13 @@ class SpectralFitting(Operation): return dataOut if path != None: sys.path.append(path) - self.library = importlib.import_module(file) + try: + self.library = importlib.import_module(file) + except: + pass + if filec != None: + self.weightf = importlib.import_module(filec) + #self.weightf = importlib.import_module('weightfit') #To be inserted as a parameter groupArray = numpy.array(groupList) @@ -3162,8 +3569,11 @@ class SpectralFitting(Operation): dataOut.data_paramC = None #Set constants - constants = self.library.setConstants(dataOut) - dataOut.constants = constants + try: + constants = self.library.setConstants(dataOut) + dataOut.constants = constants + except: + pass M = dataOut.normFactor N = dataOut.nFFTPoints ippSeconds = dataOut.ippSeconds @@ -3190,210 +3600,300 @@ class SpectralFitting(Operation): if not self.isConfig: self.isConfig = True - index = tini.tm_hour*12+tini.tm_min/5 + index = tini.tm_hour*12+tini.tm_min/taver jspc = jspc/N/N jcspc = jcspc/N/N tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.CleanRayleigh(dataOut,jspc,jcspc,2) jspectra = tmp_spectra*len(jspc[:,0,0,0]) jcspectra = tmp_cspectra*len(jspc[:,0,0,0]) - my_incoh_spectra ,my_incoh_cspectra,my_incoh_aver,my_coh_aver, incoh_spectra, coh_spectra, incoh_cspectra, coh_cspectra, incoh_aver, coh_aver = self.__DiffCoherent(jspectra, jcspectra, dataOut, noise, snrth, None, None) + 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) clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index) dataOut.data_spc = incoh_spectra dataOut.data_cspc = incoh_cspectra + #dataOut.data_spc = tmp_spectra + #dataOut.data_cspc = tmp_cspectra clean_num_aver = incoh_aver*len(jspc[:,0,0,0]) coh_num_aver = clean_coh_aver*len(jspc[:,0,0,0]) + #clean_num_aver = (numpy.zeros([nChan, nHei])+1)*len(jspc[:,0,0,0]) + #coh_num_aver = numpy.zeros([nChan, nHei])*0*len(jspc[:,0,0,0]) #List of possible combinations listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2) indCross = numpy.zeros(len(list(listComb)), dtype = 'int') - - if getSNR: - listChannels = groupArray.reshape((groupArray.size)) - listChannels.sort() - dataOut.data_SNR = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise[listChannels]) - if dataOut.data_paramC is None: - dataOut.data_paramC = numpy.zeros((nGroups*4, nHeights,2))*numpy.nan - for i in range(nGroups): - coord = groupArray[i,:] - #Input data array - data = dataOut.data_spc[coord,:,:]/(M*N) - data = data.reshape((data.shape[0]*data.shape[1],data.shape[2])) - - #Cross Spectra data array for Covariance Matrixes - ind = 0 - for pairs in listComb: - pairsSel = numpy.array([coord[x],coord[y]]) - indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0]) - ind += 1 - dataCross = dataOut.data_cspc[indCross,:,:]/(M*N) - dataCross = dataCross**2 - nhei = nHeights - poweri = numpy.sum(dataOut.data_spc[:,1:nProf-0,:],axis=1)/clean_num_aver[:,:] - if i == 0 : my_noises = numpy.zeros(4,dtype=float) #FLTARR(4) - n0i = numpy.nanmin(poweri[0+i*2,0:nhei-0])/(nProf-1) - n1i = numpy.nanmin(poweri[1+i*2,0:nhei-0])/(nProf-1) - n0 = n0i - n1= n1i - my_noises[2*i+0] = n0 - my_noises[2*i+1] = n1 - snrth = -16.0 - snrth = 10**(snrth/10.0) - - for h in range(nHeights): - d = data[:,h] - smooth = clean_num_aver[i+1,h] #dataOut.data_spc[:,1:nProf-0,:] - signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth - signalpn1 = (dataOut.data_spc[i*2+1,1:(nProf-0),h])/smooth - signal0 = signalpn0-n0 - signal1 = signalpn1-n1 - snr0 = numpy.sum(signal0/n0)/(nProf-1) - snr1 = numpy.sum(signal1/n1)/(nProf-1) - if snr0 > snrth and snr1 > snrth and clean_num_aver[i+1,h] > 0 : - #Covariance Matrix - D = numpy.diag(d**2) - ind = 0 - for pairs in listComb: - #Coordinates in Covariance Matrix - x = pairs[0] - y = pairs[1] - #Channel Index - S12 = dataCross[ind,:,h] - D12 = numpy.diag(S12) - #Completing Covariance Matrix with Cross Spectras - D[x*N:(x+1)*N,y*N:(y+1)*N] = D12 - D[y*N:(y+1)*N,x*N:(x+1)*N] = D12 - ind += 1 - diagD = numpy.zeros(256) - if h == 17 : - for ii in range(256): diagD[ii] = D[ii,ii] - #Dinv=numpy.linalg.inv(D) - #L=numpy.linalg.cholesky(Dinv) - try: - Dinv=numpy.linalg.inv(D) - L=numpy.linalg.cholesky(Dinv) - except: - Dinv = D*numpy.nan - L= D*numpy.nan - LT=L.T - - dp = numpy.dot(LT,d) - - #Initial values - data_spc = dataOut.data_spc[coord,:,h] - - if (h>0)and(error1[3]<5): - p0 = dataOut.data_param[i,:,h-1] - else: - p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants))# sin el i(data_spc, constants, i) - try: - #Least Squares - #print (dp,LT,constants) - #value =self.__residFunction(p0,dp,LT,constants) - #print ("valueREADY",value.shape, type(value)) - #optimize.leastsq(value) - minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True) - #minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants)) - #Chi square error - #print(minp,covp.infodict,mesg,ier) - #print("REALIZA OPTIMIZ") - error0 = numpy.sum(infodict['fvec']**2)/(2*N) - #Error with Jacobian - error1 = self.library.errorFunction(minp,constants,LT) -# print self.__residFunction(p0,dp,LT, constants) -# print infodict['fvec'] -# print self.__residFunction(minp,dp,LT,constants) - - except: + if Gaussian_Windowed == 1: + #dataOut.data_spc = jspectra + ''' + Written by R. Flores + ''' + print("normFactor: ", dataOut.normFactor) + data_spc_aux = numpy.copy(dataOut.data_spc)#[:,0,:] + data_spc_aux[:,0,:] = (data_spc_aux[:,1,:]+data_spc_aux[:,-1,:])/2 + #''' + from scipy.signal import medfilt + import matplotlib.pyplot as plt + dataOut.moments = numpy.ones((dataOut.nChannels,4,dataOut.nHeights))*numpy.NAN + dataOut.VelRange = dataOut.getVelRange(0) + for nChannel in range(dataOut.nChannels): + for hei in range(dataOut.heightList.shape[0]): + #print("ipp: ", dataOut.ippSeconds) + #spc = numpy.copy(dataOut.data_spc[nChannel,:,hei]) + spc = data_spc_aux[nChannel,:,hei] + if spc.all() == 0.: + print("CONTINUE") + continue + #print(VelRange) + #print(dataOut.getFreqRange(64)) + #print("Hei: ", dataOut.heightList[hei]) + + spc_mod = numpy.copy(spc) + spcm = medfilt(spc_mod,11) + spc_max = numpy.max(spcm) + dop1_x0 = dataOut.VelRange[numpy.argmax(spcm)] + #D = numpy.min(spcm) + D_in = (numpy.mean(spcm[:15])+numpy.mean(spcm[-15:]))/2. + #print("spc_max: ", spc_max) + #print("dataOut.ippSeconds: ", dataOut.ippSeconds, dataOut.timeInterval) + ##fun, A, B, C, D = self.windowing_single(spc,dataOut.VelRange,spc_max,dop1_x0,abs(dop1_x0),D,dataOut.nFFTPoints) + #fun, A, B, C, D = self.windowing_single(spc,dataOut.VelRange,spc_max,dop1_x0,abs(dop1_x0),D,dataOut.nFFTPoints) + 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) + + dataOut.moments[nChannel,0,hei] = A + dataOut.moments[nChannel,1,hei] = B + dataOut.moments[nChannel,2,hei] = C + dataOut.moments[nChannel,3,hei] = D + ''' + if nChannel == 0: + print(dataOut.heightList[hei]) + plt.figure() + plt.plot(dataOut.VelRange,spc,marker='*',linestyle='--') + plt.plot(dataOut.VelRange,fun) + plt.title(dataOut.heightList[hei]) + plt.show() + ''' + #plt.show() + #''' + dataOut.data_spc = jspectra + print("SUCCESS") + return dataOut + + elif Gaussian_Windowed == 2: #Only to clean spc + dataOut.VelRange = dataOut.getVelRange(0) + return dataOut + else: + if getSNR: + listChannels = groupArray.reshape((groupArray.size)) + listChannels.sort() + dataOut.data_SNR = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise[listChannels]) + if dataOut.data_paramC is None: + dataOut.data_paramC = numpy.zeros((nGroups*4, nHeights,2))*numpy.nan + for i in range(nGroups): + coord = groupArray[i,:] + #Input data array + data = dataOut.data_spc[coord,:,:]/(M*N) + data = data.reshape((data.shape[0]*data.shape[1],data.shape[2])) + + #Cross Spectra data array for Covariance Matrixes + ind = 0 + for pairs in listComb: + pairsSel = numpy.array([coord[x],coord[y]]) + indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0]) + ind += 1 + dataCross = dataOut.data_cspc[indCross,:,:]/(M*N) + dataCross = dataCross**2 + nhei = nHeights + poweri = numpy.sum(dataOut.data_spc[:,1:nProf-0,:],axis=1)/clean_num_aver[:,:] + if i == 0 : my_noises = numpy.zeros(4,dtype=float) + n0i = numpy.nanmin(poweri[0+i*2,0:nhei-0])/(nProf-1) + n1i = numpy.nanmin(poweri[1+i*2,0:nhei-0])/(nProf-1) + n0 = n0i + n1= n1i + my_noises[2*i+0] = n0 + my_noises[2*i+1] = n1 + snrth = -25.0 # -4 + snrth = 10**(snrth/10.0) + jvelr = numpy.zeros(nHeights, dtype = 'float') + #snr0 = numpy.zeros(nHeights, dtype = 'float') + #snr1 = numpy.zeros(nHeights, dtype = 'float') + hvalid = [0] + + 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,:]) + + for h in range(nHeights): + smooth = clean_num_aver[i+1,h] + signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth + signalpn1 = (dataOut.data_spc[i*2+1,1:(nProf-0),h])/smooth + signal0 = signalpn0-n0 + signal1 = signalpn1-n1 + snr0 = numpy.sum(signal0/n0)/(nProf-1) + snr1 = numpy.sum(signal1/n1)/(nProf-1) + #jmax0 = MAX(signal0,maxp0) + #jmax1 = MAX(signal1,maxp1) + gamma = coh2[:,h] + + indxs = (numpy.isfinite(list(gamma))==True).nonzero() + + if len(indxs) >0: + if numpy.nanmean(gamma) > 0.07: + maxp0 = numpy.argmax(signal0*gamma) + maxp1 = numpy.argmax(signal1*gamma) + #print('usa gamma',numpy.nanmean(gamma)) + else: + maxp0 = numpy.argmax(signal0) + maxp1 = numpy.argmax(signal1) + jvelr[h] = (absc[maxp0]+absc[maxp1])/2. + else: jvelr[h] = absc[0] + if snr0 > 0.1 and snr1 > 0.1: hvalid = numpy.concatenate((hvalid,h), axis=None) + #print(maxp0,absc[maxp0],snr0,jvelr[h]) + + if len(hvalid)> 1: fd0 = numpy.median(jvelr[hvalid[1:]])*-1 + else: fd0 = numpy.nan + #print(fd0,hvalid) + for h in range(nHeights): + d = data[:,h] + smooth = clean_num_aver[i+1,h] #dataOut.data_spc[:,1:nProf-0,:] + signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth + signalpn1 = (dataOut.data_spc[i*2+1,1:(nProf-0),h])/smooth + signal0 = signalpn0-n0 + signal1 = signalpn1-n1 + snr0 = numpy.sum(signal0/n0)/(nProf-1) + snr1 = numpy.sum(signal1/n1)/(nProf-1) + + if snr0 > snrth and snr1 > snrth and clean_num_aver[i+1,h] > 0 : + #Covariance Matrix + D = numpy.diag(d**2) + ind = 0 + for pairs in listComb: + #Coordinates in Covariance Matrix + x = pairs[0] + y = pairs[1] + #Channel Index + S12 = dataCross[ind,:,h] + D12 = numpy.diag(S12) + #Completing Covariance Matrix with Cross Spectras + D[x*N:(x+1)*N,y*N:(y+1)*N] = D12 + D[y*N:(y+1)*N,x*N:(x+1)*N] = D12 + ind += 1 + diagD = numpy.zeros(256) + + #Dinv=numpy.linalg.inv(D) + #L=numpy.linalg.cholesky(Dinv) + try: + Dinv=numpy.linalg.inv(D) + L=numpy.linalg.cholesky(Dinv) + except: + Dinv = D*numpy.nan + L= D*numpy.nan + LT=L.T + + dp = numpy.dot(LT,d) + + #Initial values + data_spc = dataOut.data_spc[coord,:,h] + w = data_spc/data_spc + if filec != None: + w = self.weightf.weightfit(w,tini.tm_year,tini.tm_yday,index,h,i) + + if (h>6) and (error1[3]<25): + p0 = dataOut.data_param[i,:,h-1] + #print('usa anterior') + else: + p0 = numpy.array(self.library.initialValuesFunction(data_spc*w, constants))# sin el i(data_spc, constants, i) + + if filec != None: + p0 = self.weightf.Vrfit(p0,tini.tm_year,tini.tm_yday,index,h,i) + p0[3] = fd0 + #if index == 175 and i==1 and h>=27 and h<=35: p0[3]=30 + #if h >= 6 and i==1 and h<= 10: print(p0) + try: + #Least Squares + minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True) + #minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants)) + #Chi square error + error0 = numpy.sum(infodict['fvec']**2)/(2*N) + #Error with Jacobian + error1 = self.library.errorFunction(minp,constants,LT) + #if h >= 0 and h<= 10 and i ==0: print(p0,minp,error1) + #if i>=0 and h>=0: print(index,h,minp[3]) + # print self.__residFunction(p0,dp,LT, constants) + # print infodict['fvec'] + # print self.__residFunction(minp,dp,LT,constants) + + except: + minp = p0*numpy.nan + error0 = numpy.nan + error1 = p0*numpy.nan + # s_sq = (self.__residFunction(minp,dp,LT,constants)).sum()/(len(dp)-len(p0)) + # covp = covp*s_sq + # error = [] + # for ip in range(len(minp)): + # try: + # error.append(numpy.absolute(covp[ip][ip])**0.5) + # except: + # error.append( 0.00 ) + #if i==1 and h==11 and index == 139: print(p0, minp,data_spc) + else : + data_spc = dataOut.data_spc[coord,:,h] + p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants)) minp = p0*numpy.nan error0 = numpy.nan error1 = p0*numpy.nan - #print ("EXCEPT 0000000000") -# s_sq = (self.__residFunction(minp,dp,LT,constants)).sum()/(len(dp)-len(p0)) -# covp = covp*s_sq -# #print("TRY___________________________________________1") -# error = [] -# for ip in range(len(minp)): -# try: -# error.append(numpy.absolute(covp[ip][ip])**0.5) -# except: -# error.append( 0.00 ) - else : - data_spc = dataOut.data_spc[coord,:,h] - p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants)) - minp = p0*numpy.nan - error0 = numpy.nan - error1 = p0*numpy.nan - #Save - if dataOut.data_param is None: - dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan - dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan - - dataOut.data_error[i,:,h] = numpy.hstack((error0,error1)) - dataOut.data_param[i,:,h] = minp - - for ht in range(nHeights-1) : - smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam] - dataOut.data_paramC[4*i,ht,1] = smooth - signalpn0 = (coh_spectra[i*2 ,1:(nProf-0),ht])/smooth #coh_spectra - signalpn1 = (coh_spectra[i*2+1,1:(nProf-0),ht])/smooth - - #val0 = WHERE(signalpn0 > 0,cval0) - val0 = (signalpn0 > 0).nonzero() - val0 = val0[0] - #print('hvalid:',hvalid) - #print('valid', valid) - if len(val0) == 0 : val0_npoints = nProf - else : val0_npoints = len(val0) - - #val1 = WHERE(signalpn1 > 0,cval1) - val1 = (signalpn1 > 0).nonzero() - val1 = val1[0] - if len(val1) == 0 : val1_npoints = nProf - else : val1_npoints = len(val1) - - dataOut.data_paramC[0+4*i,ht,0] = numpy.sum((signalpn0/val0_npoints))/n0 - dataOut.data_paramC[1+4*i,ht,0] = numpy.sum((signalpn1/val1_npoints))/n1 - - signal0 = (signalpn0-n0) # > 0 - vali = (signal0 < 0).nonzero() - vali = vali[0] - if len(vali) > 0 : signal0[vali] = 0 - signal1 = (signalpn1-n1) #> 0 - vali = (signal1 < 0).nonzero() - vali = vali[0] - if len(vali) > 0 : signal1[vali] = 0 - snr0 = numpy.sum(signal0/n0)/(nProf-1) - snr1 = numpy.sum(signal1/n1)/(nProf-1) - doppler = absc[1:] - if snr0 >= snrth and snr1 >= snrth and smooth : - signalpn0_n0 = signalpn0 - signalpn0_n0[val0] = signalpn0[val0] - n0 - mom0 = self.moments(doppler,signalpn0-n0,nProf) -# sigtmp= numpy.transpose(numpy.tile(signalpn0, [4,1])) -# momt= self.__calculateMoments( sigtmp, doppler , n0 ) - signalpn1_n1 = signalpn1 - signalpn1_n1[val1] = signalpn1[val1] - n1 - mom1 = self.moments(doppler,signalpn1_n1,nProf) - dataOut.data_paramC[2+4*i,ht,0] = (mom0[0]+mom1[0])/2. - dataOut.data_paramC[3+4*i,ht,0] = (mom0[1]+mom1[1])/2. -# if graph == 1 : -# window, 13 -# plot,doppler,signalpn0 -# oplot,doppler,signalpn1,linest=1 -# oplot,mom0(0)*doppler/doppler,signalpn0 -# oplot,mom1(0)*doppler/doppler,signalpn1 -# print,interval/12.,beam,45+ht*15,snr0,snr1,mom0(0),mom1(0),mom0(1),mom1(1) - #ENDIF - #ENDIF - #ENDFOR End height - - dataOut.data_spc = jspectra - if getSNR: - listChannels = groupArray.reshape((groupArray.size)) - listChannels.sort() - - dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:], my_noises[listChannels]) - return dataOut + + if dataOut.data_param is None: + dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan + dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan + + dataOut.data_error[i,:,h] = numpy.hstack((error0,error1)) + dataOut.data_param[i,:,h] = minp + + for ht in range(nHeights-1) : + smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam] + dataOut.data_paramC[4*i,ht,1] = smooth + signalpn0 = (coh_spectra[i*2 ,1:(nProf-0),ht])/smooth #coh_spectra + signalpn1 = (coh_spectra[i*2+1,1:(nProf-0),ht])/smooth + + val0 = (signalpn0 > 0).nonzero() + val0 = val0[0] + + if len(val0) == 0 : val0_npoints = nProf + else : val0_npoints = len(val0) + + val1 = (signalpn1 > 0).nonzero() + val1 = val1[0] + if len(val1) == 0 : val1_npoints = nProf + else : val1_npoints = len(val1) + + dataOut.data_paramC[0+4*i,ht,0] = numpy.sum((signalpn0/val0_npoints))/n0 + dataOut.data_paramC[1+4*i,ht,0] = numpy.sum((signalpn1/val1_npoints))/n1 + + signal0 = (signalpn0-n0) + vali = (signal0 < 0).nonzero() + vali = vali[0] + if len(vali) > 0 : signal0[vali] = 0 + signal1 = (signalpn1-n1) + vali = (signal1 < 0).nonzero() + vali = vali[0] + if len(vali) > 0 : signal1[vali] = 0 + snr0 = numpy.sum(signal0/n0)/(nProf-1) + snr1 = numpy.sum(signal1/n1)/(nProf-1) + doppler = absc[1:] + if snr0 >= snrth and snr1 >= snrth and smooth : + signalpn0_n0 = signalpn0 + signalpn0_n0[val0] = signalpn0[val0] - n0 + mom0 = self.moments(doppler,signalpn0-n0,nProf) + + signalpn1_n1 = signalpn1 + signalpn1_n1[val1] = signalpn1[val1] - n1 + mom1 = self.moments(doppler,signalpn1_n1,nProf) + dataOut.data_paramC[2+4*i,ht,0] = (mom0[0]+mom1[0])/2. + dataOut.data_paramC[3+4*i,ht,0] = (mom0[1]+mom1[1])/2. + + dataOut.data_spc = jspectra + if getSNR: + listChannels = groupArray.reshape((groupArray.size)) + listChannels.sort() + + dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:], my_noises[listChannels]) + return dataOut def __residFunction(self, p, dp, LT, constants): @@ -6628,9 +7128,9 @@ class IGRFModel(Operation): mkfact_short_2020_2.mkfact(dataOut.year,dataOut.h,dataOut.bfm,dataOut.thb,dataOut.bki,dataOut.MAXNRANGENDT) #mkfact_short_2020.mkfact(dataOut.year,dataOut.h,dataOut.bfm,dataOut.thb,dataOut.bki,dataOut.MAXNRANGENDT) - print("bki: ", dataOut.bki[:10]) - print("thb: ", dataOut.thb[:10]) - print("bfm: ", dataOut.bfm[:10]) + #print("bki: ", dataOut.bki[:10]) + #print("thb: ", dataOut.thb[:10]) + #print("bfm: ", dataOut.bfm[:10]) #print("IDs: ", id(dataOut.bki)) return dataOut @@ -6794,8 +7294,10 @@ class MergeProc(ProcessingUnit): self.dataOut.NSHTS = int(numpy.shape(self.dataOut.ph2)[0]) dH = self.dataOut.heightList[1]-self.dataOut.heightList[0] dH /= self.dataOut.windowOfFilter - self.dataOut.heightList = numpy.arange(0,self.dataOut.NSHTS)*dH + self.dataOut.heightList = numpy.arange(0,self.dataOut.NSHTS)*dH + dH + #print("heightList: ", self.dataOut.heightList) self.dataOut.NDP = self.dataOut.NSHTS + #exit(1) #print(self.dataOut.heightList) class MST_Den_Conv(Operation): diff --git a/schainpy/model/proc/jroproc_spectra.py b/schainpy/model/proc/jroproc_spectra.py index 1b79952..b06cbc2 100644 --- a/schainpy/model/proc/jroproc_spectra.py +++ b/schainpy/model/proc/jroproc_spectra.py @@ -455,19 +455,43 @@ class GetSNR(Operation): def run(self,dataOut): - noise = dataOut.getNoise() - maxdB = 16 - - normFactor = 24 - + #noise = dataOut.getNoise() + noise = dataOut.getNoise(ymin_index=-10) #Región superior donde solo debería de haber ruido + #print("Noise: ", noise) + #print("Noise_dB: ", 10*numpy.log10(noise/dataOut.normFactor)) + #print("Heights: ", dataOut.heightList) #dataOut.data_snr = (dataOut.data_spc.sum(axis=1))/(noise[:,None]*dataOut.normFactor) - dataOut.data_snr = (dataOut.data_spc.sum(axis=1))/(noise[:,None]*dataOut.nFFTPoints) - - dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.5, numpy.nan, dataOut.data_snr) + ################dataOut.data_snr = (dataOut.data_spc.sum(axis=1))/(noise[:,None]*dataOut.nFFTPoints) #Before 12Jan2023 + #dataOut.data_snr = (dataOut.data_spc.sum(axis=1)-noise[:,None])/(noise[:,None]) + dataOut.data_snr = (dataOut.data_spc.sum(axis=1)-noise[:,None]*dataOut.nFFTPoints)/(noise[:,None]*dataOut.nFFTPoints) #It works apparently + dataOut.snl = numpy.log10(dataOut.data_snr) + #print("snl: ", dataOut.snl) + #exit(1) + #print(dataOut.heightList[-11]) + #print(numpy.shape(dataOut.heightList)) + #print(dataOut.data_snr) + #print(dataOut.data_snr[0,-11]) + #exit(1) + #dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.5, numpy.nan, dataOut.data_snr) + #dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.1, numpy.nan, dataOut.data_snr) + #dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.0, numpy.nan, dataOut.data_snr) + #dataOut.data_snr = numpy.where(dataOut.data_snr<.05, numpy.nan, dataOut.data_snr) + dataOut.snl = numpy.where(dataOut.data_snr<.01, numpy.nan, dataOut.snl) + ''' + import matplotlib.pyplot as plt + #plt.plot(10*numpy.log10(dataOut.data_snr[0]),dataOut.heightList) + plt.plot(dataOut.data_snr[0],dataOut.heightList)#,marker='*') + plt.xlim(-1,10) + plt.axvline(1,color='k') + plt.axvline(.1,color='k',linestyle='--') + plt.grid() + plt.show() + ''' #dataOut.data_snr = 10*numpy.log10(dataOut.data_snr) #dataOut.data_snr = numpy.expand_dims(dataOut.data_snr,axis=0) #print(dataOut.data_snr.shape) #exit(1) + #print("Before: ", dataOut.data_snr[0]) return dataOut diff --git a/schainpy/model/proc/jroproc_spectra_lags_faraday.py b/schainpy/model/proc/jroproc_spectra_lags_faraday.py index fd91c32..c546b37 100644 --- a/schainpy/model/proc/jroproc_spectra_lags_faraday.py +++ b/schainpy/model/proc/jroproc_spectra_lags_faraday.py @@ -3815,6 +3815,9 @@ class IncohInt(Operation): dataOut.flagNoData = False dataOut.VelRange = dataOut.getVelRange(0) + dataOut.FreqRange = dataOut.getFreqRange(0)/1000. + #print("VelRange: ", dataOut.VelRange) + #exit(1) #print("Power",numpy.sum(dataOut.data_spc[0,:,20:30],axis=0)) #print("Power",numpy.sum(dataOut.data_spc[0,100:110,:],axis=1)) @@ -3921,7 +3924,7 @@ class SnrFaraday(Operation): return dataOut -class SpectraDataToFaraday_07_11_2022(Operation): +class SpectraDataToFaraday(Operation): #ISR MODE ''' Written by R. Flores ''' @@ -4144,10 +4147,14 @@ class SpectraDataToFaraday_07_11_2022(Operation): #print(data_eej) index_eej = CleanCohEchoes.mad_based_outlier(self,data_eej[:17]) aux_eej = numpy.array(index_eej.nonzero()).ravel() + print("aux_eej: ", aux_eej) + if aux_eej != []: + dataOut.min_id_eej = aux_eej[-1] + else: + dataOut.min_id_eej = 12 - dataOut.min_id_eej = aux_eej[-1] - print(dataOut.min_id_eej) + #print("min_id_eej: ", dataOut.min_id_eej) #exit(1) def run(self,dataOut): @@ -4204,7 +4211,7 @@ class SpectraDataToFaraday_07_11_2022(Operation): #exit(1) return dataOut -class SpectraDataToFaraday(Operation): +class SpectraDataToFaraday_MST(Operation): #MST MODE """Operation to use spectra data in Faraday processing. Parameters: diff --git a/schainpy/model/proc/jroproc_voltage.py b/schainpy/model/proc/jroproc_voltage.py index b875ce5..4b7b497 100644 --- a/schainpy/model/proc/jroproc_voltage.py +++ b/schainpy/model/proc/jroproc_voltage.py @@ -38,7 +38,8 @@ class VoltageProc(ProcessingUnit): if self.dataIn.type == 'Voltage': self.dataOut.copy(self.dataIn) self.dataOut.runNextUnit = runNextUnit - #print(self.dataOut.data.shape) + #print("data shape: ", self.dataOut.data.shape) + #exit(1) def __updateObjFromAmisrInput(self): @@ -226,6 +227,7 @@ class selectHeights(Operation): #print(dataOut.nHeights) #exit(1) ''' + #print("select Heights: ", self.dataOut.heightList) return self.dataOut def selectHeightsByIndex(self, minIndex, maxIndex): @@ -1250,7 +1252,7 @@ class LagsReshapeLP(Operation): self.buffer[3,:,:,:] = numpy.transpose(numpy.conj(self.lagp3),(1,2,0)) #print(self.buffer[3,:,100,15]) - print(sum(self.buffer[3,:,199,2])) + print("Sum: ", sum(self.buffer[3,:,199,2])) #print(self.cnorm) exit(1) @@ -3348,7 +3350,7 @@ class DoublePulseACFs_PerLag(Operation): exit(1) ''' #print(pa) - #''' + ''' import matplotlib.pyplot as plt #plt.plot(dataOut.p[:,-1],dataOut.heightList) plt.plot(pa/dataOut.pan-1.,dataOut.heightList) @@ -3358,7 +3360,7 @@ class DoublePulseACFs_PerLag(Operation): plt.show() #print("p: ",dataOut.p[33,:]) #exit(1) - #''' + ''' return dataOut class FaradayAngleAndDPPower(Operation): @@ -3399,7 +3401,7 @@ class FaradayAngleAndDPPower(Operation): for i in range(dataOut.MAXNRANGENDT): dataOut.range1[i]=dataOut.H0 + i*dataOut.DH dataOut.h2[i]=dataOut.range1[i]**2 - print("shape ph2",numpy.shape(dataOut.ph2)) + #print("shape ph2",numpy.shape(dataOut.ph2)) for j in range(dataOut.NDP): dataOut.ph2[j]=0. dataOut.sdp2[j]=0. @@ -3505,7 +3507,7 @@ class ElectronDensityFaraday(Operation): ndphi=dataOut.NSHTS-4 #print(dataOut.phi) #exit(1) - ''' + #''' if hasattr(dataOut, 'flagSpreadF') and dataOut.flagSpreadF: #if dataOut.flagSpreadF: nanindex = numpy.argwhere(numpy.isnan(dataOut.phi)) @@ -3517,7 +3519,7 @@ class ElectronDensityFaraday(Operation): else: #dataOut.phi_uwrp = dataOut.phi.copy() dataOut.phi[:]=numpy.unwrap(dataOut.phi[:]) #Better results - ''' + #''' #print(dataOut.phi) #print(dataOut.ph2) #exit(1) @@ -3968,12 +3970,12 @@ class NormalizeDPPowerRoberto_V2(Operation): if hasattr(dataOut, 'flagSpreadF') and dataOut.flagSpreadF: i2=int((620-dataOut.range1[0])/dataOut.DH) nanindex = numpy.argwhere(numpy.isnan(dataOut.ph2)) - print("nanindex",nanindex) + #print("nanindex",nanindex) i1 = nanindex[-1][0] #VER CUANDO i1>i2 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 #print("i1, i2",i1,i2) - print(dataOut.heightList) - print("Bounds: ", dataOut.heightList[i1],dataOut.heightList[i2]) + #print(dataOut.heightList) + #print("Bounds: ", dataOut.heightList[i1],dataOut.heightList[i2]) #print(dataOut.dphi[33]) #print(dataOut.ph2[33]) #print(dataOut.dphi[i1::]) @@ -4008,9 +4010,9 @@ class NormalizeDPPowerRoberto_V2(Operation): #''' #if (time_text.hour == 5 and time_text.minute == 32): #Year: 2022, DOY:104 #if (time_text.hour == 0 and time_text.minute == 12): #Year: 2022, DOY:93 - 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 + #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 #if (time_text.hour == 1 and time_text.minute == 23) or (time_text.hour == 1 and time_text.minute == 44): #Year: 2022, DOY:243 - dataOut.cf = dataOut.cflast[0] + #dataOut.cf = dataOut.cflast[0] #if (time_text.hour == 0 and time_text.minute == 4): #Year: 2022, DOY:244 #dataOut.cf = 0.08 #print("here") @@ -4020,12 +4022,30 @@ class NormalizeDPPowerRoberto_V2(Operation): #dataOut.cf = 0.09 #if (time_text.hour == 3 and time_text.minute == 59) or (time_text.hour == 4 and time_text.minute == 20): #Year: 2022, DOY:244 #dataOut.cf = 0.09 + + #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 + #dataOut.cf = 0.029844088 + #if (time_text.hour == 0 and time_text.minute == 12) or (time_text.hour == 0 and time_text.minute == 22): #Year: 2023, DOY:028 + #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 + #dataOut.cf = dataOut.cflast[0] + #if (time_text.hour == 7 and time_text.minute == 29): #Year: 2023, DOY:029 + #dataOut.cf = 0.025670702 + #if (time_text.hour == 8 and time_text.minute == 12): #Year: 2023, DOY:029 + #dataOut.cf = 0.036969288 + #if (time_text.hour == 8 and time_text.minute == 22): #Year: 2023, DOY:029 + #dataOut.cf = 0.0418645 + 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 + dataOut.cf = dataOut.cflast[0] + #if (time_text.hour == 0 and time_text.minute == 22) or (time_text.hour == 0 and time_text.minute == 33): #Year: 2023, DOY:031 + #if (time_text.hour == 7 and time_text.minute == 29): #Year: 2023, DOY:032 + #dataOut.cf = dataOut.cflast[0] + #''' #dataOut.cf = 0.000057#0.0008136899 dataOut.cflast[0]=dataOut.cf - print("cf: ", dataOut.cf) + #print("cf: ", dataOut.cf) #print(dataOut.ph2) #print(dataOut.sdp2) @@ -4049,7 +4069,7 @@ class NormalizeDPPowerRoberto_V2(Operation): #print(dataOut.ph2) #print(dataOut.sdp2) #input() - print("shape before" ,numpy.shape(dataOut.ph2)) + #print("shape before" ,numpy.shape(dataOut.ph2)) return dataOut @@ -4227,15 +4247,18 @@ class DPTemperaturesEstimation(Operation): eb=numpy.resize(eb,10) dataOut.ifit=numpy.resize(dataOut.ifit,10) + #print("*********************FITACF_FIT*********************",dataOut.covinv,e,dataOut.params,eb,dataOut.m) + #exit(1) with suppress_stdout_stderr(): 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) # - + #exit(1) if dataOut.params[2]>dataOut.params[1]*1.05: dataOut.ifit[2]=0 dataOut.params[1]=dataOut.params[2]=t1 with suppress_stdout_stderr(): 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) # - + #print("*********************FIT SUCCESS*********************",dataOut.covinv,e,dataOut.params,eb,dataOut.m) + #exit(1) if (dataOut.ifit[2]==0): dataOut.params[2]=dataOut.params[1] if (dataOut.ifit[3]==0 and iflag==0): @@ -4306,7 +4329,7 @@ class DPTemperaturesEstimation(Operation): return dataOut -class DenCorrection(NormalizeDPPowerRoberto_V2): +class DenCorrection_bckp_Dec_2022(NormalizeDPPowerRoberto_V2): ''' Written by R. Flores ''' @@ -4322,6 +4345,7 @@ class DenCorrection(NormalizeDPPowerRoberto_V2): def TeTiEstimation(self,dataOut): y=numpy.zeros(dataOut.DPL,order='F',dtype='float32') + #y_aux = numpy.zeros(1,,dtype='float32') for i in range(dataOut.NSHTS): y[0]=y[1]=dataOut.range1[i] @@ -4423,7 +4447,7 @@ class DenCorrection(NormalizeDPPowerRoberto_V2): fion[2]=0.0 #0 for j in range(dataOut.DPL): tau=dataOut.alag[j]*1.0e-3 - + #print("***********ACF2***********") with suppress_stdout_stderr():#The smoothness in range of "y" depends on the smoothness of the input parameters 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) #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): return dataOut +class DenCorrection(NormalizeDPPowerRoberto_V2): + ''' + Written by R. Flores + ''' + def __init__(self, **kwargs): + + Operation.__init__(self, **kwargs) + self.aux = 0 + + def gaussian(self, x, a, b, c): + val = a * numpy.exp(-(x - b)**2 / (2*c**2)) + return val + + def TeTiEstimation(self,dataOut): + + dataOut.DPL = 2 #for MST + y=numpy.zeros(dataOut.DPL,order='F',dtype='float32') + + #y_aux = numpy.zeros(1,,dtype='float32') + for i in range(dataOut.NSHTS): + y[0]=y[1]=dataOut.range1[i] + + y = y.astype(dtype='float64',order='F') + three=int(3) + wl = 3.0 + tion=numpy.zeros(three,order='F',dtype='float32') + fion=numpy.zeros(three,order='F',dtype='float32') + nui=numpy.zeros(three,order='F',dtype='float32') + wion=numpy.zeros(three,order='F',dtype='int32') + bline=0.0 + #bline=numpy.zeros(1,order='F',dtype='float32') + my_aux = numpy.ones(dataOut.NSHTS,order='F',dtype='float32') + acf_Temps = numpy.ones(dataOut.NSHTS,order='F',dtype='float32')*numpy.nan + acf_no_Temps = numpy.ones(dataOut.NSHTS,order='F',dtype='float32')*numpy.nan + + from scipy import signal + + def func(params): + return (ratio2-self.gaussian(dataOut.heightList[:dataOut.NSHTS],params[0],params[1],params[2])) + + #print("Before loop") + dataOut.info2[0] = 1 + for i in range(dataOut.NSHTS): + if dataOut.info2[i]==1: + angle=dataOut.thb[i]*0.01745 + nue=nui[0]=nui[1]=nui[2]=0.0#nui[3]=0.0 + wion[0]=16 #O + wion[1]=1 #H + wion[2]=4 #He + tion[0]=tion[1]=tion[2]=dataOut.ti2[i] + #tion[0]=tion[1]=tion[2]=ti2_smooth[i] + fion[0]=1.0-dataOut.phy2[i] #1 + fion[1]=dataOut.phy2[i] #0 + fion[2]=0.0 #0 + for j in range(dataOut.DPL): + tau=dataOut.alag[j]*1.0e-3 + #print("***********ACF2***********") + with suppress_stdout_stderr():#The smoothness in range of "y" depends on the smoothness of the input parameters + 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) + #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) + #print("i: ", i, "j: ", j) + #if i == 0 and j == 1: + #print("Params: ",wl,tau,dataOut.te2[i],tion,fion,nue,nui,wion,angle,dataOut.ph2[i],dataOut.bfm[i],y[j],three) + #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) + #exit(1) + #if dataOut.ut_Faraday>11.0 and dataOut.range1[i]>150.0 and dataOut.range1[i]<400.0: + + if dataOut.ut_Faraday>11.0 and dataOut.range1[i]>150.0 and dataOut.range1[i]<300.0: + #if dataOut.ut_Faraday>11.0 and dataOut.range1[i]>150.0 and dataOut.range1[i]<400.0: + tau=0.0 + with suppress_stdout_stderr(): + bline=fitacf_acf2.acf2(wl,tau,tion,tion,fion,nue,nui,wion,angle,dataOut.ph2[i],dataOut.bfm[i],bline,three) + + #if i == 0 and j == 1: + #print("bline: ",bline) + #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) + #exit(1) + cf=min(1.2,max(1.0,bline/y[0])) #FACTOR DE EFICIENCIA + #cf = bline/y[0] + #cf=min(2.,max(1.0,bline/y[0])) + my_aux[i] = cf + acf_Temps[i] = y[0] + acf_no_Temps[i] = bline + #dataOut.ph2[i]=cf*dataOut.ph2[i] #Instead we adjust the curve "cf" into a Gaussian, + #dataOut.sdp2[i]=cf*dataOut.sdp2[i] #in order to get smoother values of density + for j in range(1,dataOut.DPL): + #y[j]=(y[j]/y[0])*dataOut.DH+dataOut.range1[i] + y[j]=min(max((y[j]/y[0]),-1.0),1.0)*dataOut.DH+dataOut.range1[i] + y[0]=dataOut.range1[i]+dataOut.DH + + + ratio = my_aux-1 + #ratio = dataOut.te2[:dataOut.NSHTS]/dataOut.ti2[:dataOut.NSHTS] + def lsq_func(params): + return (ratio-self.gaussian(dataOut.heightList[:dataOut.NSHTS],params[0],params[1],params[2])) + + x0_value = numpy.array([max(ratio),250,20]) + + popt = least_squares(lsq_func,x0=x0_value,verbose=0) + + A = popt.x[0]; B = popt.x[1]; C = popt.x[2] + + aux = self.gaussian(dataOut.heightList[:dataOut.NSHTS], A, B, C) + 1 #ratio + 1 + + dataOut.ph2[:dataOut.NSHTS]*=aux + dataOut.sdp2[:dataOut.NSHTS]*=aux + #dataOut.ph2[:26]*=aux[:26] + #dataOut.sdp2[:26]*=aux[:26] + #print(aux) + #print("inside correction",dataOut.ph2) + + def run(self,dataOut): + #print("hour",gmtime(dataOut.utctime).tm_hour) + if gmtime(dataOut.utctime).tm_hour < 24. and gmtime(dataOut.utctime).tm_hour >= 11.: + #print("inside") + self.TeTiEstimation(dataOut) + dataOut.flagTeTiCorrection = True + + self.normalize(dataOut) + + return dataOut + class DataPlotCleaner(Operation): ''' Written by R. Flores @@ -4825,10 +4971,16 @@ class DataSaveCleaner(Operation): #print(time_text.hour,time_text.minute) #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 - 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 + #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 #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 #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 #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 + #if (time_text.hour == 0 and time_text.minute==8) or (time_text.hour == 0 and time_text.minute==18): #Year: 2023, DOY:027 + #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 + #if (time_text.hour >=2 and time_text.hour<5) or (time_text.hour == 10): #Year: 2023, DOY:029 + 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 + #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 + #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 dataOut.DensityFinal[0,:]=missing dataOut.EDensityFinal[0,:]=missing @@ -4842,9 +4994,19 @@ class DataSaveCleaner(Operation): dataOut.flagNoData = True #Remueve todo el perfil #''' ''' - #if (time_text.hour >= 7 and time_text.hour < 9): #Year: 2022, DOY:102 - if (time_text.hour == 20 and time_text.minute == 8): #Year: 2022, DOY:243 - id_aux = 35 + if (time_text.hour == 5) or (time_text.hour == 6) or (time_text.hour == 7) or (time_text.hour == 9): #Year: 2023, DOY:032 + id_aux = 11 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + + if (time_text.hour == 6 and time_text.minute == 57) or (time_text.hour == 7 and time_text.minute == 29): #Year: 2023, DOY:032 + id_aux = 36 dataOut.DensityFinal[0,id_aux:]=missing dataOut.EDensityFinal[0,id_aux:]=missing dataOut.ElecTempFinal[0,id_aux:]=missing @@ -4853,8 +5015,9 @@ class DataSaveCleaner(Operation): dataOut.EIonTempFinal[0,id_aux:]=missing dataOut.PhyFinal[0,id_aux:]=missing dataOut.EPhyFinal[0,id_aux:]=missing - if (time_text.hour == 20 and time_text.minute == 19): #Year: 2022, DOY:243 - id_aux = 33 + + if (time_text.hour == 9 and time_text.minute == 24): #Year: 2023, DOY:032 + id_aux = 30 dataOut.DensityFinal[0,id_aux:]=missing dataOut.EDensityFinal[0,id_aux:]=missing dataOut.ElecTempFinal[0,id_aux:]=missing @@ -4863,8 +5026,9 @@ class DataSaveCleaner(Operation): dataOut.EIonTempFinal[0,id_aux:]=missing dataOut.PhyFinal[0,id_aux:]=missing dataOut.EPhyFinal[0,id_aux:]=missing - if (time_text.hour == 20 and time_text.minute == 29): #Year: 2022, DOY:243 - id_aux = 30 + + if (time_text.hour == 5 and time_text.minute == 32): #Year: 2023, DOY:032 + id_aux = 41 dataOut.DensityFinal[0,id_aux:]=missing dataOut.EDensityFinal[0,id_aux:]=missing dataOut.ElecTempFinal[0,id_aux:]=missing @@ -4873,8 +5037,9 @@ class DataSaveCleaner(Operation): dataOut.EIonTempFinal[0,id_aux:]=missing dataOut.PhyFinal[0,id_aux:]=missing dataOut.EPhyFinal[0,id_aux:]=missing - if (time_text.hour == 20 and time_text.minute == 44): #Year: 2022, DOY:243 - id_aux = 31 + + if (time_text.hour == 11 and time_text.minute == 14): #Year: 2023, DOY:032 + id_aux = 35 dataOut.DensityFinal[0,id_aux:]=missing dataOut.EDensityFinal[0,id_aux:]=missing dataOut.ElecTempFinal[0,id_aux:]=missing @@ -4883,8 +5048,10 @@ class DataSaveCleaner(Operation): dataOut.EIonTempFinal[0,id_aux:]=missing dataOut.PhyFinal[0,id_aux:]=missing dataOut.EPhyFinal[0,id_aux:]=missing - if (time_text.hour <= 8): #Year: 2022, DOY:243 - id_aux = 11 + ''' + ''' + if (time_text.hour == 23 and time_text.minute == 50) or (time_text.hour == 0): #Year: 2023, DOY:031 + id_aux = 18 dataOut.DensityFinal[0,:id_aux]=missing dataOut.EDensityFinal[0,:id_aux]=missing dataOut.ElecTempFinal[0,:id_aux]=missing @@ -4893,7 +5060,8 @@ class DataSaveCleaner(Operation): dataOut.EIonTempFinal[0,:id_aux]=missing dataOut.PhyFinal[0,:id_aux]=missing dataOut.EPhyFinal[0,:id_aux]=missing - if (time_text.hour == 23): #Year: 2022, DOY:243 + + 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 id_aux = 12 dataOut.DensityFinal[0,:id_aux]=missing dataOut.EDensityFinal[0,:id_aux]=missing @@ -4903,28 +5071,66 @@ class DataSaveCleaner(Operation): dataOut.EIonTempFinal[0,:id_aux]=missing dataOut.PhyFinal[0,:id_aux]=missing dataOut.EPhyFinal[0,:id_aux]=missing - if (time_text.hour == 5 and time_text.minute == 21): #Year: 2022, DOY:243 - id_aux = (36,37) - dataOut.DensityFinal[0,id_aux]=missing - dataOut.EDensityFinal[0,id_aux]=missing - dataOut.ElecTempFinal[0,id_aux]=missing - dataOut.EElecTempFinal[0,id_aux]=missing - dataOut.IonTempFinal[0,id_aux]=missing - dataOut.EIonTempFinal[0,id_aux]=missing - dataOut.PhyFinal[0,id_aux]=missing - dataOut.EPhyFinal[0,id_aux]=missing - if (time_text.hour == 5 and time_text.minute == 53): #Year: 2022, DOY:243 - id_aux = (37,38) - dataOut.DensityFinal[0,id_aux]=missing - dataOut.EDensityFinal[0,id_aux]=missing - dataOut.ElecTempFinal[0,id_aux]=missing - dataOut.EElecTempFinal[0,id_aux]=missing - dataOut.IonTempFinal[0,id_aux]=missing - dataOut.EIonTempFinal[0,id_aux]=missing - dataOut.PhyFinal[0,id_aux]=missing - dataOut.EPhyFinal[0,id_aux]=missing - if (time_text.hour == 6 and time_text.minute == 4): #Year: 2022, DOY:243 - id_aux = (38,39) + + if (time_text.hour == 11 and time_text.minute == 45): #Year: 2023, DOY:031 + id_aux = 8 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + + if (time_text.hour == 13): #Year: 2023, DOY:031 + id_aux = 37 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + + if (time_text.hour == 15 and time_text.minute == 18): #Year: 2023, DOY:031 + id_aux = 34 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + + if (time_text.hour == 10 and time_text.minute == 52): #Year: 2023, DOY:031 + id_aux = 29 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + + if (time_text.hour == 10 and time_text.minute == 9): #Year: 2023, DOY:031 + id_aux = 33 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + ''' + + #''' + if (time_text.hour == 0 and time_text.minute == 33): #Year: 2023, DOY:030 + id_aux = 38 dataOut.DensityFinal[0,id_aux]=missing dataOut.EDensityFinal[0,id_aux]=missing dataOut.ElecTempFinal[0,id_aux]=missing @@ -4933,8 +5139,9 @@ class DataSaveCleaner(Operation): dataOut.EIonTempFinal[0,id_aux]=missing dataOut.PhyFinal[0,id_aux]=missing dataOut.EPhyFinal[0,id_aux]=missing - if (time_text.hour == 12 and time_text.minute == 6): #Year: 2022, DOY:243 - id_aux = (29,30) + + if (time_text.hour == 18 and time_text.minute == 41): #Year: 2023, DOY:030 + id_aux = 49 dataOut.DensityFinal[0,id_aux]=missing dataOut.EDensityFinal[0,id_aux]=missing dataOut.ElecTempFinal[0,id_aux]=missing @@ -4943,8 +5150,9 @@ class DataSaveCleaner(Operation): dataOut.EIonTempFinal[0,id_aux]=missing dataOut.PhyFinal[0,id_aux]=missing dataOut.EPhyFinal[0,id_aux]=missing - if (time_text.hour == 14 and time_text.minute == 14): #Year: 2022, DOY:243 - id_aux = (35,36) + + if (time_text.hour == 14 and time_text.minute == 14): #Year: 2023, DOY:030 + id_aux = 34 dataOut.DensityFinal[0,id_aux]=missing dataOut.EDensityFinal[0,id_aux]=missing dataOut.ElecTempFinal[0,id_aux]=missing @@ -4953,8 +5161,9 @@ class DataSaveCleaner(Operation): dataOut.EIonTempFinal[0,id_aux]=missing dataOut.PhyFinal[0,id_aux]=missing dataOut.EPhyFinal[0,id_aux]=missing - if (time_text.hour == 23 and time_text.minute == 2): #Year: 2022, DOY:243 - id_aux = (41,42) + + if (time_text.hour == 16 and time_text.minute == 1): #Year: 2023, DOY:030 + id_aux = (40,41) dataOut.DensityFinal[0,id_aux]=missing dataOut.EDensityFinal[0,id_aux]=missing dataOut.ElecTempFinal[0,id_aux]=missing @@ -4963,8 +5172,9 @@ class DataSaveCleaner(Operation): dataOut.EIonTempFinal[0,id_aux]=missing dataOut.PhyFinal[0,id_aux]=missing dataOut.EPhyFinal[0,id_aux]=missing - if (time_text.hour == 0 and time_text.minute == 8): #Year: 2022, DOY:243 - id_aux = 33 + + if (time_text.hour == 5 and time_text.minute == 53) or (time_text.hour == 15 and time_text.minute == 40): #Year: 2023, DOY:030 + id_aux = 40 dataOut.DensityFinal[0,id_aux:]=missing dataOut.EDensityFinal[0,id_aux:]=missing dataOut.ElecTempFinal[0,id_aux:]=missing @@ -4973,8 +5183,434 @@ class DataSaveCleaner(Operation): dataOut.EIonTempFinal[0,id_aux:]=missing dataOut.PhyFinal[0,id_aux:]=missing dataOut.EPhyFinal[0,id_aux:]=missing - id_aux = 18 - dataOut.DensityFinal[0,id_aux]=missing + + 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 + id_aux = 37 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + + if (time_text.hour == 7 and time_text.minute == 40): #Year: 2023, DOY:030 + id_aux = 31 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + + if (time_text.hour == 11 and time_text.minute == 2): #Year: 2023, DOY:030 + id_aux = 24 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + + if (time_text.hour == 5) or (time_text.hour == 6) or (time_text.hour == 7) or (time_text.hour == 0): #Year: 2023, DOY:030 + id_aux = 12 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + + if (time_text.hour == 23 and time_text.minute == 40) or (time_text.hour == 23 and time_text.minute == 50): #Year: 2023, DOY:030 + id_aux = 13 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + + if (time_text.hour == 11 and time_text.minute == 56): #Year: 2023, DOY:030 + id_aux = 8 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + #''' + + ''' + if (time_text.hour == 14 and time_text.minute == 42): #Year: 2023, DOY:027 + id_aux = 31 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + + if (time_text.hour == 15 and time_text.minute == 25) or (time_text.hour == 15 and time_text.minute == 36): #Year: 2023, DOY:027 + id_aux = 35 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + + 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 + id_aux = 13 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + ''' + + ''' + 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 + id_aux = 14 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + + 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 + id_aux = 12 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + + if (time_text.hour == 5 and time_text.minute == 32): #Year: 2023, DOY:029 + id_aux = 16 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + + if (time_text.hour == 0 and time_text.minute > 22 and time_text.minute != 54): #Year: 2023, DOY:029 + id_aux = 22 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + + if (time_text.hour == 12 and time_text.minute == 49): #Year: 2023, DOY:029 + id_aux = 32 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + + 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 + id_aux = 37 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + + 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 + id_aux = 35 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + + if (time_text.hour == 7 and time_text.minute == 18): #Year: 2023, DOY:029 + id_aux = 36 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + + if (time_text.hour == 9 and time_text.minute == 26) or (time_text.hour == 9 and time_text.minute == 37): #Year: 2023, DOY:029 + id_aux = 28 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + + if (time_text.hour == 22 and time_text.minute == 14) or (time_text.hour == 22 and time_text.minute == 25): #Year: 2023, DOY:029 + id_aux = (37,38) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + ''' + + ''' + 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 + id_aux = 13 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + + if (time_text.hour == 9 and time_text.minute == 5): #Year: 2023, DOY:028 + id_aux = 11 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + + if (time_text.hour == 11 and time_text.minute == 13) or (time_text.hour == 12 and time_text.minute == 17): #Year: 2023, DOY:028 + id_aux = 34 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + + if (time_text.hour == 7 and time_text.minute == 29): #Year: 2023, DOY:028 + id_aux = 35 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + + 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 + id_aux = 37 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + + if (time_text.hour == 0 and time_text.minute == 22) or (time_text.hour == 0 and time_text.minute == 33): #Year: 2023, DOY:028 + id_aux = 22 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + + 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 + id_aux = 12 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + ''' + + ''' + #if (time_text.hour >= 7 and time_text.hour < 9): #Year: 2022, DOY:102 + if (time_text.hour == 20 and time_text.minute == 8): #Year: 2022, DOY:243 + id_aux = 35 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + if (time_text.hour == 20 and time_text.minute == 19): #Year: 2022, DOY:243 + id_aux = 33 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + if (time_text.hour == 20 and time_text.minute == 29): #Year: 2022, DOY:243 + id_aux = 30 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + if (time_text.hour == 20 and time_text.minute == 44): #Year: 2022, DOY:243 + id_aux = 31 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + if (time_text.hour <= 8): #Year: 2022, DOY:243 + id_aux = 11 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + if (time_text.hour == 23): #Year: 2022, DOY:243 + id_aux = 12 + dataOut.DensityFinal[0,:id_aux]=missing + dataOut.EDensityFinal[0,:id_aux]=missing + dataOut.ElecTempFinal[0,:id_aux]=missing + dataOut.EElecTempFinal[0,:id_aux]=missing + dataOut.IonTempFinal[0,:id_aux]=missing + dataOut.EIonTempFinal[0,:id_aux]=missing + dataOut.PhyFinal[0,:id_aux]=missing + dataOut.EPhyFinal[0,:id_aux]=missing + if (time_text.hour == 5 and time_text.minute == 21): #Year: 2022, DOY:243 + id_aux = (36,37) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 5 and time_text.minute == 53): #Year: 2022, DOY:243 + id_aux = (37,38) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 6 and time_text.minute == 4): #Year: 2022, DOY:243 + id_aux = (38,39) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 12 and time_text.minute == 6): #Year: 2022, DOY:243 + id_aux = (29,30) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 14 and time_text.minute == 14): #Year: 2022, DOY:243 + id_aux = (35,36) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 23 and time_text.minute == 2): #Year: 2022, DOY:243 + id_aux = (41,42) + dataOut.DensityFinal[0,id_aux]=missing + dataOut.EDensityFinal[0,id_aux]=missing + dataOut.ElecTempFinal[0,id_aux]=missing + dataOut.EElecTempFinal[0,id_aux]=missing + dataOut.IonTempFinal[0,id_aux]=missing + dataOut.EIonTempFinal[0,id_aux]=missing + dataOut.PhyFinal[0,id_aux]=missing + dataOut.EPhyFinal[0,id_aux]=missing + if (time_text.hour == 0 and time_text.minute == 8): #Year: 2022, DOY:243 + id_aux = 33 + dataOut.DensityFinal[0,id_aux:]=missing + dataOut.EDensityFinal[0,id_aux:]=missing + dataOut.ElecTempFinal[0,id_aux:]=missing + dataOut.EElecTempFinal[0,id_aux:]=missing + dataOut.IonTempFinal[0,id_aux:]=missing + dataOut.EIonTempFinal[0,id_aux:]=missing + dataOut.PhyFinal[0,id_aux:]=missing + dataOut.EPhyFinal[0,id_aux:]=missing + id_aux = 18 + dataOut.DensityFinal[0,id_aux]=missing dataOut.EDensityFinal[0,id_aux]=missing dataOut.ElecTempFinal[0,id_aux]=missing dataOut.EElecTempFinal[0,id_aux]=missing @@ -5371,7 +6007,7 @@ class DataSaveCleaner(Operation): dataOut.PhyFinal[0,id_aux:]=missing dataOut.EPhyFinal[0,id_aux:]=missing ''' - #''' + ''' #print("den: ", dataOut.DensityFinal[0,27]) if (time_text.hour == 5 and time_text.minute == 42): #Year: 2022, DOY:242 id_aux = 16 @@ -5493,13 +6129,13 @@ class DataSaveCleaner(Operation): dataOut.EIonTempFinal[0,id_aux]=missing dataOut.PhyFinal[0,id_aux]=missing dataOut.EPhyFinal[0,id_aux]=missing - #''' + ''' #print("den_final",dataOut.DensityFinal) dataOut.flagNoData = numpy.all(numpy.isnan(dataOut.DensityFinal)) #Si todos los valores son NaN no se prosigue - ####dataOut.flagNoData = False #Solo para ploteo + #dataOut.flagNoData = False #Descomentar solo para ploteo #Comentar para MADWriter dataOut.DensityFinal *= 1.e6 #Convert units to m^⁻3 dataOut.EDensityFinal *= 1.e6 #Convert units to m^⁻3 @@ -6319,7 +6955,7 @@ class removeDCHAE(Operation): return dataOut -class Decoder(Operation): +class Decoder_Original(Operation): isConfig = False __profIndex = 0 @@ -6355,10 +6991,13 @@ class Decoder(Operation): self.__nChannels = dataOut.nChannels self.__nProfiles = dataOut.nProfiles self.__nHeis = dataOut.nHeights - + #print("self.__nHeis: ", self.__nHeis) + #print("self.nBaud: ", self.nBaud) + #exit(1) if self.__nHeis < self.nBaud: raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud)) - + #print("JJE") + #exit(1) #Frequency __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex) @@ -6421,6 +7060,10 @@ class Decoder(Operation): #print(profilesList) for i in range(self.__nChannels): for j in profilesList: + #print("data.shape: ", data.shape) + #print("code_block: ", code_block.shape) + #print("corr: ", numpy.shape(numpy.correlate(data[i,j,:], code_block[j,:], mode='full'))) + #exit(1) self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:] return self.datadecTime @@ -6444,7 +7087,8 @@ class Decoder(Operation): if dataOut.flagDecodeData: print("This data is already decoded, recoding again ...") - + #print("code: ", numpy.shape(code)) + #exit(1) if not self.isConfig: if code is None: @@ -6522,6 +7166,234 @@ class Decoder(Operation): return dataOut # dataOut.flagDeflipData = True #asumo q la data no esta sin flip +class Decoder(Operation): + + isConfig = False + __profIndex = 0 + + code = None + + nCode = None + nBaud = None + + def __init__(self, **kwargs): + + Operation.__init__(self, **kwargs) + + self.times = None + self.osamp = None + # self.__setValues = False + self.isConfig = False + self.setupReq = False + def setup(self, code, osamp, dataOut): + + self.__profIndex = 0 + + self.code = code + + self.nCode = len(code) + #self.nBaud = len(code[0]) + self.nBaud = int(numpy.shape(code)[-1]) + + if (osamp != None) and (osamp >1): + self.osamp = osamp + self.code = numpy.repeat(code, repeats=self.osamp, axis=1) + self.nBaud = self.nBaud*self.osamp + + self.__nChannels = dataOut.nChannels + self.__nProfiles = dataOut.nProfiles + self.__nHeis = dataOut.nHeights + #print("self.__nHeis: ", self.__nHeis) + #print("self.nBaud: ", self.nBaud) + #exit(1) + if self.__nHeis < self.nBaud: + raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud)) + #print("JJE") + #exit(1) + ''' + #Frequency + __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex) + + __codeBuffer[:,0:self.nBaud] = self.code + + self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1)) + ''' + if dataOut.flagDataAsBlock: + + self.ndatadec = self.__nHeis #- self.nBaud + 1 + + self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex) + + else: + + #Time + self.ndatadec = self.__nHeis #- self.nBaud + 1 + + + self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex) + + def __convolutionInFreq(self, data): + + fft_code = self.fft_code[self.__profIndex].reshape(1,-1) + + fft_data = numpy.fft.fft(data, axis=1) + + conv = fft_data*fft_code + + data = numpy.fft.ifft(conv,axis=1) + + return data + + def __convolutionInFreqOpt(self, data): + + raise NotImplementedError + + def __convolutionInTime(self, data): + + code = self.code[self.__profIndex] + for i in range(self.__nChannels): + #aux=numpy.correlate(data[i,:], code, mode='full') + #print(numpy.shape(aux)) + #print(numpy.shape(data[i,:])) + #print(numpy.shape(code)) + #exit(1) + self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:] + + return self.datadecTime + + def __convolutionByBlockInTime(self, data, AutoDecod): + + if not AutoDecod: + repetitions = int(self.__nProfiles / self.nCode) + junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize)) + junk = junk.flatten() + code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud)) + else: + code_block = self.code + profilesList = range(self.__nProfiles) + #print(numpy.shape(self.datadecTime)) + #print(numpy.shape(data)) + #print(profilesList) + for i in range(self.__nChannels): + for j in profilesList: + #print("data.shape: ", data.shape) + #print("code_block: ", code_block.shape) + #print("corr: ", numpy.shape(numpy.correlate(data[i,j,:], code_block[j,:], mode='full'))) + #exit(1) + if not AutoDecod: + self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:] + else: + if i%2: + self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[0,j,:], mode='full')[self.nBaud-1:] + else: + self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[1,j,:], mode='full')[self.nBaud-1:] + + return self.datadecTime + + def __convolutionByBlockInFreq(self, data): + + raise NotImplementedError("Decoder by frequency fro Blocks not implemented") + + + fft_code = self.fft_code[self.__profIndex].reshape(1,-1) + + fft_data = numpy.fft.fft(data, axis=2) + + conv = fft_data*fft_code + + data = numpy.fft.ifft(conv,axis=2) + + return data + + + def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None, AutoDecod = 0): + + if dataOut.flagDecodeData: + print("This data is already decoded, recoding again ...") + #print("code: ", numpy.shape(code)) + #exit(1) + if not self.isConfig: + + if code is None and not AutoDecod: + if dataOut.code is None: + raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type) + + code = dataOut.code + else: + if not AutoDecod: + code = numpy.array(code).reshape(nCode,nBaud) + else: + po = 0 + print("***********AutoDecod***********") + code = dataOut.data[:2,:,po:po+64] + #print("AutoDecod Shape: ", numpy.shape(code)) + #exit(1) + self.setup(code, osamp, dataOut) + + self.isConfig = True + + if mode == 3: + sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode) + + if times != None: + sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n") + + if self.code is None: + print("Fail decoding: Code is not defined.") + return + + self.__nProfiles = dataOut.nProfiles + datadec = None + + if mode == 3: + mode = 0 + + if dataOut.flagDataAsBlock: + """ + Decoding when data have been read as block, + """ + + if mode == 0: + datadec = self.__convolutionByBlockInTime(dataOut.data, AutoDecod) + if mode == 1: + datadec = self.__convolutionByBlockInFreq(dataOut.data) + else: + """ + Decoding when data have been read profile by profile + """ + if mode == 0: + datadec = self.__convolutionInTime(dataOut.data) + + if mode == 1: + datadec = self.__convolutionInFreq(dataOut.data) + + if mode == 2: + datadec = self.__convolutionInFreqOpt(dataOut.data) + + if datadec is None: + raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode) + + dataOut.code = self.code + dataOut.nCode = self.nCode + dataOut.nBaud = self.nBaud + + dataOut.data = datadec#/self.nBaud + + #print("before",dataOut.heightList) + dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]] + #print("after",dataOut.heightList) + + dataOut.flagDecodeData = True #asumo q la data esta decodificada + + if self.__profIndex == self.nCode-1: + self.__profIndex = 0 + return dataOut + + self.__profIndex += 1 + + return dataOut + + class DecoderRoll(Operation): isConfig = False @@ -6883,7 +7755,19 @@ class ProfileSelector(Operation): dataOut.nProfiles = len(profileList) dataOut.profileIndex = dataOut.nProfiles - 1 dataOut.flagNoData = False - + #print("Shape after prof select: ", numpy.shape(dataOut.data)) + #print(dataOut.heightList) + #exit(1) + ''' + po = 5 + import matplotlib.pyplot as plt + this_data = dataOut.data[0,0,:] + plt.plot(this_data*numpy.conjugate(this_data),marker='*',linestyle='-') + plt.axvline(po) + plt.axvline(po+64-1) + plt.grid() + plt.show() + ''' return dataOut """ @@ -7520,7 +8404,7 @@ class CrossProdHybrid(CrossProdDP): #exit(1) if i==0: self.lagp0[n][j][self.bcounter-1]=numpy.sum(c[:dataOut.NSCAN]) - ww[n,j,:,self.bcounter-1]=c[:dataOut.NSCAN] + #ww[n,j,:,self.bcounter-1]=c[:dataOut.NSCAN] self.lagp3[n][j][self.bcounter-1]=numpy.sum(c[dataOut.NSCAN:]/self.cnorm) elif i==1: self.lagp1[n][j][self.bcounter-1]=numpy.sum(c[:dataOut.NSCAN]) @@ -7540,8 +8424,8 @@ class CrossProdHybrid(CrossProdDP): #print(self.cnorm) #exit(1) #print("self,lagp0: ", self.lagp0[0,0,self.bcounter-1]) - print(ww[:,0,0,self.bcounter-1]) - exit(1) + #print(ww[:,0,0,self.bcounter-1]) + #exit(1) def LP_median_estimates_original(self,dataOut): @@ -7998,9 +8882,9 @@ class LongPulseAnalysis(Operation): #print(anoise0) #print(anoise1) #exit(1) - print("nis: ", dataOut.nis) - print("pan: ", dataOut.pan) - print("pbn: ", dataOut.pbn) + #print("nis: ", dataOut.nis) + #print("pan: ", dataOut.pan) + #print("pbn: ", dataOut.pbn) #print(numpy.sum(dataOut.output_LP_integrated[0,:,0])) ''' import matplotlib.pyplot as plt @@ -8008,8 +8892,8 @@ class LongPulseAnalysis(Operation): plt.show() ''' #print(dataOut.output_LP_integrated[0,40,0]) - print(numpy.sum(dataOut.output_LP_integrated[:,0,0])) - exit(1) + #print(numpy.sum(dataOut.output_LP_integrated[:,0,0])) + #exit(1) #################### PROBAR MÁS INTEGRACIÓN, SINO MODIFICAR VALOR DE "NIS" #################### # VER dataOut.nProfiles_LP # @@ -8237,7 +9121,7 @@ class LongPulseAnalysis(Operation): print(numpy.sum(dataOut.te2)) exit(1) ''' - print("Success") + #print("Success 1") ###################Correlation pulse and itself #print(dataOut.NRANGE)