* Last processed by NICE on 23-SEP-1991 14:18:37 * Customized for : VAX VMS SUBROUTINE POLYNO(LPLOT,ERROR,LAST) C---------------------------------------------------------------------- C SAS Internal routine for command C BASE [N] [/PLOT] C Computes and remove a polynomial baseline C This version uses NAG routine E02ADF instead of VADEME C Arguments : C LPLOT L Plot the baseline or not Input C ERROR L Logical error flag C LAST L Use last baseline determination. C (2-Jan-1985) C 17-Apr-1988 Use single precision NAG routines. C---------------------------------------------------------------------- LOGICAL LPLOT,ERROR,LAST * INCLUDE 'inc:pi.inc' INCLUDE 'inc:constant.inc' INCLUDE 'rdata.inc' INCLUDE 'structure.par' INCLUDE 'setup.inc' EXTERNAL IPLOT * REAL*4 Z(MDATA),SMIN,VV,S0,S1,S2,VIT,DEL,ALIM,S12,FILLIN INTEGER IW,I,NNN,MMM, IFAIL, NP, I1, I2, LCHAIN, I3, IMIN INTEGER W1(MDATA),W2 INTEGER MMAX PARAMETER (MMAX=30) REAL*4 A(MMAX,MMAX),WORK2(2,MMAX),S(MMAX),A1(MMAX) REAL*4 W(MDATA),X(MDATA),Y(MDATA),XMIN,XMAX,YMIN,YMAX INTEGER WORK1,IPW,POINTER,MEMORY(1) CHARACTER*80 CHAIN * SAVE * * First set weights NP = 0 IF (RKIND.NE.KIND_ONOFF) THEN IF (RDATAX(CNCHAN).GT.RDATAX(1)) THEN I1 = 1 I2 = CNCHAN I3 = 1 ELSE I1 = CNCHAN I2 = 1 I3 = -1 ENDIF DO 30 I=I1, I2, I3 W1(I) = 0 IF (I.LT.CIMIN .OR. I.GT.CIMAX) GOTO 30 IF (RNWIND.GT.0) THEN VV = RDATAX(I) DO IW = 1, RNWIND IF ((VV-RW1(IW))*(VV-RW2(IW)) .LT. 0.) THEN IF (I.NE.CIMIN .AND. I.NE.CIMAX) THEN GOTO 30 ELSEIF (.NOT.LAST .AND. RDEG.GT.1) THEN CALL MESSAGE(6,2,'POLYNO', $ 'Baseline extrapolation is hazardous.') GOTO 30 ELSE GOTO 30 ENDIF ENDIF ENDDO ENDIF * * Avoid bad channels IF (RDATAY(I).NE.CBAD) THEN NP = NP + 1 W(NP)= +1.0 X(NP) = RDATAX(I) Y(NP) = RDATAY(I) ENDIF W1(I) = 1 30 CONTINUE * Handle the special case of continuum on/off: group the data points in pairs, * and ignore windows. ELSE NP = 1 X(1) = RDATAX(1) Y(1) = (RDATAY(1)+RDATAY(2))/2 W(1) = EPSR4 DO I=1,CNCHAN/2 IF (RDATAY(2*I-1).NE.CBAD.AND.RDATAY(2*I).NE.CBAD) THEN NP = NP + 1 W(NP)= +1.0 X(NP) = (RDATAX(2*I-1)+RDATAX(2*I))/2 Y(NP) = (RDATAY(2*I-1)+RDATAY(2*I))/2 ENDIF W1(2*I) = 1 W1(2*I-1) = 1 ENDDO NP = NP+1 X(NP) = RDATAX(CNCHAN) Y(NP) = (RDATAY(CNCHAN)+RDATAY(CNCHAN-1))/2 W(NP) = EPSR4 ENDIF * * Determine the baseline IF (.NOT.LAST) THEN MMM = RDEG XMIN = X(1) XMAX = X(NP) IFAIL = 1 CALL GET_VM(3*MDATA,WORK1,*99) IPW = POINTER(WORK1,MEMORY) CALL E02ADE(NP,MMM+1,MMAX,X,Y,W,MEMORY(IPW),WORK2,A,S,IFAIL) CALL FREE_VM(3*MDATA,WORK1) IF (IFAIL.NE.0) THEN WRITE(CHAIN,'(A,I)')'NAG Error in E02ADF, ifail = ',IFAIL CALL SIC_BLANC(CHAIN,LCHAIN) CALL MESSAGE(8,4,'POLYNO',CHAIN(1:LCHAIN)) ERROR = .TRUE. RETURN ENDIF IMIN = 1 SMIN = S(1) DO I=1, MMM+1 IF (S(I).LE.SMIN) THEN SMIN = S(I) IMIN = I ENDIF A1(I) = A(MMM+1,I) ENDDO IF (IMIN.LT.MMM) THEN WRITE(CHAIN,101) IMIN LCHAIN = 31 CALL MESSAGE(6,2,'POLYNO',CHAIN(1:LCHAIN)) ENDIF ENDIF * IF (MMM.GT.1) THEN * Degree > 1 * Compute polynomial in fitting range, constant value outside that range. IFAIL = 1 X(1) = -1.0 CALL E02AEE(MMM+1, A1, X(1), YMIN, IFAIL) IFAIL = 1 X(1) = 1.0 CALL E02AEE(MMM+1, A1, X(1), YMAX, IFAIL) DO I=1, CNCHAN X(1) = ((RDATAX(I)-XMIN)-(XMAX-RDATAX(I)))/(XMAX-XMIN) IF (X(1).LE.-1.0) THEN Y(1) = YMIN ELSEIF (X(1).GE.1.0) THEN Y(1) = YMAX ELSE IFAIL = 1 CALL E02AEE(MMM+1, A1, X(1), Y(1), IFAIL) IF (IFAIL.NE.0) THEN WRITE(CHAIN,'(A,I)') 'NAG Error in E02AEF, ifail = ' $ ,IFAIL CALL SIC_BLANC(CHAIN,LCHAIN) CALL MESSAGE(8,4,'POLYNO',CHAIN(1:LCHAIN)) ERROR = .TRUE. RETURN ENDIF ENDIF Z(I) = Y(1) ENDDO ELSE * Degree 0 or 1 : Compute polynomial everywhere IFAIL = 1 X(1) = -1.0 CALL E02AEE(MMM+1, A1, X(1), YMIN, IFAIL) IFAIL = 1 X(1) = 1.0 CALL E02AEE(MMM+1, A1, X(1), YMAX, IFAIL) DO I=1, CNCHAN X(1) = ((RDATAX(I)-XMIN)-(XMAX-RDATAX(I)))/(XMAX-XMIN) Y(1) = YMIN + (X(1)+1.0)*(YMAX-YMIN)/2.0 Z(I) = Y(1) ENDDO ENDIF * * Plot the baseline if needed IF (LPLOT) THEN NNN = CNCHAN IF (RDEG.LE.1) THEN CALL IPLOT(1.,Z(1),3) CALL IPLOT(FLOAT(NNN),Z(NNN),2) ELSE CALL CONNE2(1.,1.,1.,Z,NNN,IPLOT) ENDIF ENDIF * * Subtract the baseline and find the sigma S0 = 0. S1 = 0. S2 = 0. DO I = 1, CNCHAN * * Avoid bad channels IF (RDATAY(I).NE.CBAD) THEN RDATAY(I) = RDATAY(I)-Z(I) IF (W1(I).NE.0) THEN S0 = S0 + W1(I) S1 = S1 + W1(I)*RDATAY(I) S2 = S2 + W1(I)*RDATAY(I)**2 ENDIF ENDIF ENDDO IF (S0.NE.0) THEN S1 = S1 / S0 S2 = S2 / S0 RSIGFI = SQRT(S2-S1**2) WRITE(6,*) 'Polynomial baseline of order ',RDEG, $ ' r.m.s. is ',RSIGFI,' K' ELSE WRITE(6,*) 'Cannot measure r.m.s noise' ENDIF * * Aire, Vitesse, Largeur * RAIRE = 0. VIT = 0. DEL = 0. ALIM = RSIGFI * 3.0 S0 = 0. S1 = 0. S2 = 0. DO I = CIMIN, CIMAX * * Avoid bad channels W2 = 1 - W1(I) IF (W2.NE.0) THEN IF (RDATAY(I).NE.CBAD) THEN RAIRE = RAIRE + RDATAY(I) * W2 IF (RDATAY(I).GE.ALIM) THEN S0 = S0 + RDATAY(I) * W2 S1 = S1 + RDATAY(I) * I * W2 S2 = S2 + RDATAY(I) * FLOAT(I)**2 * W2 ENDIF ELSE RAIRE = RAIRE + FILLIN(RDATAY,I,CIMIN,CIMAX,CBAD) * W2 ENDIF ENDIF ENDDO * * Type out result according to Kind IF (RKIND.EQ.KIND_SPEC) THEN RAIRE = RAIRE * ABS (RVRES) IF (S0.NE.0) THEN S1 = S1 / S0 VIT = (S1-RRCHAN) * RVRES + RVOFF S2 = S2 / S0 S12 = S1**2 IF (S2.LE.S12) S2 = S12 DEL = ABS (RVRES) * SQRT((S2-S12)*8.*ALOG(2.)) ENDIF WRITE(6,*) 'Area ',RAIRE,' Vel.',VIT,' Width ',DEL ELSE RAIRE = RAIRE * ABS (RARES) IF (S0.NE.0) THEN S1 = S1 / S0 VIT = (S1-RRPOIN) * RARES + RAREF S2 = S2 / S0 S12 = S1**2 IF (S2.LE.S12) S2 = S12 DEL = ABS (RARES) * SQRT((S2-S12)*8.*ALOG(2.)) ENDIF WRITE(6,*) 'Area ',RAIRE,' Pos.',VIT,' Width ',DEL ENDIF RETURN * 99 CALL MESSAGE(8,3,'BASE','Insufficient memory for work space') ERROR = .TRUE. 101 FORMAT('Degree ',I2,' would be even better.') END * FUNCTION FILLIN(R,IVAL,IMIN,IMAX,BAD) C---------------------------------------------------------------------- C SAS Internal routine C Interpolate bad channel (if possible) C Arguments C R R*4(*) Array to be interpolated C IVAL I Pixel to be interpolated C IMIN I First pixel in array C IMAX I Last pixel in array C BAD R*4 Blanking value C Created 21-Mar-1986 S.Guilloteau C---------------------------------------------------------------------- INTEGER IVAL,IMIN,IMAX REAL R(IMAX),FILLIN,BAD INTEGER I,IF,IL * DO I=IVAL-1,IMIN,-1 IF (R(I).NE.BAD) GOTO 10 ENDDO * * Rien en dessous, tout au dessus DO I=IVAL+1,IMAX-1 IF (R(I).NE.BAD) GOTO 10 ENDDO * * Tout mauvais, sauf peut-etre FILLIN = R(IMAX) RETURN * * Un bon au dessus ou au dessous 10 IF = I DO I=MAX(IF+1,IVAL+1),IMAX IF (R(I).NE.BAD) GOTO 20 ENDDO * Un seul bon IF (IF.GT.IVAL .OR. IF.EQ.IMIN) THEN FILLIN = R(IF) RETURN ENDIF DO I=IF-1,IMIN,-1 IF (R(I).NE.BAD) GOTO 20 ENDDO FILLIN = R(IF) RETURN * * Deux bons 20 IL = I FILLIN = ( (IVAL-IF)*R(IL) + (IL-IVAL)*R(IF) ) / FLOAT(IL-IF) RETURN END