SUBROUTINE EVLBND(SYSFIL,BNDFIL,OUTFIL,ERRCD,ERRMSG) C C FUNCTION: CF CF C USAGE: CU CU C INPUTS: CI CI C OUTPUTS: CO CO C ALGORITHM: CA CA C MACHINE DEPENDENCIES: CM CM C HISTORY: CH CH written by: CH date: CH current version: CH modifications: CH added dpcom: 7/16/88 jdb CH C ROUTINES CALLED: CC CC C COMMON MEMORY USED: CM CM DPCOM -- see dpcommon.f and dpcom.f CM C---------------------------------------------------------------------- C written for: The CASCADE Project C Oak Ridge National Laboratory C U.S. Department of Energy C contract number DE-AC05-840R21400 C subcontract number 37B-07685C S13 C organization: The University of Tennessee C---------------------------------------------------------------------- C THIS SOFTWARE IS IN THE PUBLIC DOMAIN C NO RESTRICTIONS ON ITS USE ARE IMPLIED C---------------------------------------------------------------------- C C INCLUDE 'Parameter.f' C CHARACTER*(*) SYSFIL CHARACTER*(*) BNDFIL CHARACTER*(*) OUTFIL CHARACTER*(*) ERRMSG C CHARACTER*3 SVX CHARACTER*80 TITLE C DOUBLE PRECISION A(SIZE,SIZE) DOUBLE PRECISION B(SIZE,SIZE) DOUBLE PRECISION C(SIZE,SIZE) DOUBLE PRECISION D(SIZE,SIZE) C DOUBLE PRECISION ASAV(SIZE,SIZE) DOUBLE PRECISION BSAV(SIZE,SIZE) DOUBLE PRECISION CSAV(SIZE,SIZE) DOUBLE PRECISION DSAV(SIZE,SIZE) C DOUBLE PRECISION MINPLT DOUBLE PRECISION MAXPLT DOUBLE PRECISION WCMAX DOUBLE PRECISION TEMP DOUBLE PRECISION LBFPR(2,SIZE2) DOUBLE PRECISION HBFPR(2,SIZE2) DOUBLE PRECISION WL, WH C DOUBLE PRECISION FREQ(NPTPLT) DOUBLE PRECISION SIGMIN(NPTPLT) DOUBLE PRECISION SIGMAX(NPTPLT) C DOUBLE PRECISION WRCGML DOUBLE PRECISION WWRSCL DOUBLE PRECISION WRCGMH DOUBLE PRECISION WWRSCH DOUBLE PRECISION AWCMAX DOUBLE PRECISION AWCMIN DOUBLE PRECISION SMAXSG DOUBLE PRECISION SMINSG DOUBLE PRECISION GHMINF DOUBLE PRECISION GHMAXF C DOUBLE COMPLEX G(SIZE,SIZE) DOUBLE COMPLEX AINVB(SIZE,SIZE) DOUBLE COMPLEX H(SIZE,SIZE) DOUBLE COMPLEX ZC(SIZE) DOUBLE PRECISION ZR(SIZE) DOUBLE PRECISION EVRE(SIZE) DOUBLE PRECISION EVIM(SIZE) INTEGER IPVT(SIZE) C INTEGER NSTATS INTEGER NINPS INTEGER NOUTS INTEGER ERRCD C INTEGER NW1 INTEGER NB INTEGER MINI C INTEGER NBNDL INTEGER NBNDH INTEGER IBND INTEGER JBND C INCLUDE 'dpcom.f' C C--get system description C CALL INSYS(SYSFIL,NINPS,NOUTS,NSTATS, 2 SIZE,SIZE,SIZE,ASAV,BSAV,CSAV,DSAV,ERRCD) IF (ERRCD .NE. 0) THEN ERRMSG = 'EVLBND: Fatal error from INSYS '// 1 'accessing '//SYSFIL CLOSE (UNIT=UNIT1) RETURN END IF CLOSE (UNIT=UNIT1) C CALL SAVE (SIZE,SIZE,NSTATS,NSTATS,ASAV,A) CALL SAVE (SIZE,SIZE,NSTATS,NINPS, BSAV,B) CALL SAVE (SIZE,SIZE,NOUTS, NSTATS,CSAV,C) CALL SAVE (SIZE,SIZE,NOUTS, NINPS, DSAV,D) C C--read and sort bounds C OPEN (UNIT=UNIT1,FILE=BNDFIL,STATUS='OLD',ERR=9991) C READ (UNIT1,*,END=9991,ERR=9991) MINPLT,MAXPLT READ (UNIT1,*,END=9991,ERR=9991) NBNDL C IF (NBNDL .GT. SIZE2) THEN ERRCD = 100 ERRMSG = 'EVLBND: Fatal error; Too many'// 1 ' low frequency bounds.' CLOSE (UNIT=UNIT1) RETURN END IF C IF (NBNDL .LT. 1) NBNDL = 1 DO 10, IBND = 1, NBNDL READ (UNIT1,*,END=9991,ERR=9991) 2 (LBFPR(INDEX,IBND),INDEX=1,2) 10 CONTINUE READ (UNIT1,*,END=9991,ERR=9991) NBNDH C IF (NBNDH .GT. SIZE2) THEN ERRCD = 200 ERRMSG = 'EVLBND: Fatal error; Too many high'// 1 ' frequency bounds.' CLOSE (UNIT=UNIT1) RETURN END IF C IF (NBNDH .LT. 1) NBNDH = 1 DO 20, IBND = 1, NBNDH READ (UNIT1,*,END=9991,ERR=9991) 2 (HBFPR(INDEX,IBND),INDEX=1,2) 20 CONTINUE READ (1,*,END=9991,ERR=9991) WCMAX C CLOSE (UNIT=UNIT1) C C--SORT BY FREQUENCY C IF (NBNDL .GT. 1) THEN DO 30, IBND = 1, NBNDL - 1 DO 30, JBND = IBND + 1, NBNDL IF (LBFPR(2,IBND) .GT. 2 LBFPR(2,JBND)) THEN C TEMP = LBFPR(2,IBND) LBFPR(2,IBND) = 2 LBFPR(2,JBND) LBFPR(2,JBND) = TEMP C TEMP = LBFPR(1,IBND) LBFPR(1,IBND) = 2 LBFPR(1,JBND) LBFPR(1,JBND) = TEMP END IF 30 CONTINUE END IF C IF (NBNDH .GT. 1) THEN DO 40, IBND = 1, NBNDH - 1 DO 40, JBND = IBND + 1, NBNDH IF (HBFPR(2,IBND) .GT. 2 HBFPR(2,JBND)) THEN C TEMP = HBFPR(2,IBND) HBFPR(2,IBND) = 2 HBFPR(2,JBND) HBFPR(2,JBND) = TEMP C TEMP = HBFPR(1,IBND) HBFPR(1,IBND) = 2 HBFPR(1,JBND) HBFPR(1,JBND) = TEMP END IF 40 CONTINUE END IF C C--modify min and max plot frequencies if necessary C TEMP = MIN(LBFPR(2,1), 1 HBFPR(2,1)) IF (MINPLT .GT. TEMP) MINPLT = TEMP C TEMP = MAX (LBFPR(2,NBNDL), 1 HBFPR(2,NBNDH)) IF (MAXPLT .LT. TEMP) MAXPLT = TEMP C C--calculate min and max singular values of return difference C--over range minplt ... lbfpr(2,n_bounds_l) C OPEN (UNIT=UNIT8,STATUS='SCRATCH',ERR=9992) C TEMP = LBFPR(2,NBNDL) IF (TEMP .LE. MINPLT) TEMP = MINPLT * 10.D0 C CALL HFRQPT(SIZE,NSTATS,NINPS,NOUTS, . A,B,C,D,G,ZR,EVRE,EVIM, . AINVB,H,ZC,IPVT,NPTPLT,MINPLT,TEMP,ERRCD,ERRMSG) IF (ERRCD .GT. 0) THEN ERRMSG = 'EVLBND Fatal Error: '//ERRMSG CLOSE (UNIT=UNIT8) RETURN END IF C SVX = 'SVD' C CALL SRSPG (SVX,ERRCD,ERRMSG) IF (ERRCD .GT. 0) THEN ERRMSG = 'EVLBND Fatal Error: '//ERRMSG CLOSE (UNIT=UNIT8) RETURN END IF C C--read in the min and max singular values, and associated frequencies, C--from unit 8, which was left open and rewound by srspg C READ (UNIT8,1003,END=9992,ERR=9992) TITLE 1003 FORMAT (A50) C READ (UNIT8,1004,END=9992,ERR=9992) WL, WH, NW1 1004 FORMAT (2F20.6,I5) C READ (UNIT8,1006,END=9992,ERR=9992) (FREQ(I),SIGMIN(I), 2 SIGMAX(I),I=1, NPTPLT) 1006 FORMAT (3F14.0) C REWIND UNIT8 C C--find worst case low frequency gain margin and frequency at C--which it occurs. C WRCGML = 1.D30 WWRSCL = 0.D0 C NB = 1 C DO 500, I = 1, NPTPLT 50 IF (NB .GT. NBNDL) GO TO 501 IF (FREQ(I) .GT. LBFPR(2,NB)) THEN NB = NB + 1 GO TO 50 END IF C C--convert singular value to dB C IF (SIGMIN(I).GT.0.D0) THEN SIGMIN(I) = 20.D0*DLOG10(SIGMIN(I)) ELSE SIGMIN(I) = -1000.D0 END IF C TEMP = SIGMIN(I) - LBFPR(1,NB) IF (TEMP .LT. WRCGML) THEN WRCGML = TEMP WWRSCL = FREQ(I) END IF C 500 CONTINUE 501 CONTINUE C C--find largest frequency at which minimum singular value of return C--difference is greater than 2, guaranteeing that the minimum C--singular value of the return ratio is greater than one. C I = NPTPLT C 60 IF (SIGMIN(I) .LT. 2.D0 .AND. I .GT. 1) THEN I = I - 1 GO TO 60 END IF C GHMINF = FREQ(I) C C--calculate min and max singular values of inverse return difference C--over range hbfpr(2,1) ... maxplt C TEMP = HBFPR(2,1) IF (TEMP .GE. MAXPLT) TEMP = MAXPLT / 10.D0 C CALL SAVE (SIZE,SIZE,NSTATS,NSTATS,ASAV,A) CALL SAVE (SIZE,SIZE,NSTATS,NINPS, BSAV,B) CALL SAVE (SIZE,SIZE,NOUTS, NSTATS,CSAV,C) CALL SAVE (SIZE,SIZE,NOUTS, NINPS, DSAV,D) C CALL HFRQPT(SIZE,NSTATS,NINPS,NOUTS, . A,B,C,D,G,ZR,EVRE,EVIM, . AINVB,H,ZC,IPVT,NPTPLT,TEMP,MAXPLT,ERRCD,ERRMSG) IF (ERRCD .GT. 0) THEN ERRMSG = 'EVLBND Fatal Error: '//ERRMSG CLOSE (UNIT=UNIT8) RETURN END IF C SVX = 'SVI' C CALL SRSPG (SVX,ERRCD,ERRMSG) IF (ERRCD .GT. 0) THEN ERRMSG = 'EVLBND Fatal Error: '//ERRMSG CLOSE (UNIT=UNIT8) RETURN END IF C C--read in the min and max singular values, and associated frequencies, C--from unit 8, which was left open and rewound by srspg C READ (UNIT8,1003,END=9992,ERR=9992) TITLE C READ (UNIT8,1004,END=9992,ERR=9992) WL, WH, NW1 C READ (UNIT8,1006,END=9992,ERR=9992) (FREQ(I),SIGMIN(I), 2 SIGMAX(I),I=1, NPTPLT) C REWIND UNIT8 C C--find worst case high frequency gain margin and frequency at C--WHICH IT OCCURS. C WRCGMH = 1.D30 WWRSCH = 0.D0 C NB = NBNDH C DO 505, I = NPTPLT, 1, -1 70 IF (NB .LT. 1) GO TO 506 IF (FREQ(I) .LT. HBFPR(2,NB)) THEN NB = NB - 1 GO TO 70 END IF C C--convert singular value to dB C IF (SIGMIN(I).GT.0.D0) THEN SIGMIN(I) = 20.D0*DLOG10(SIGMIN(I)) ELSE SIGMIN(I) = -1000.D0 END IF C TEMP = SIGMIN(I) - HBFPR(1,NB) IF (TEMP .LT. WRCGMH) THEN WRCGMH = TEMP WWRSCH = FREQ(I) END IF C 505 CONTINUE 506 CONTINUE C C--find smallest frequency at which minimum singular value of inverse C--return difference is greater than 2, guaranteeing that the maximum C--singular value of the return ratio is less than one. C I = 1 C 80 IF (SIGMIN(I) .LT. 2.D0 .AND. I .LT. NPTPLT) THEN I = I + 1 GO TO 80 END IF C GHMAXF = FREQ(I) C C--calculate frequencies at which minimum and maximum singular values C--cross 0 dB; these occur between gh_min_freq and gh_max_freq. C 600 IF (GHMINF .GE. GHMAXF) THEN IF (MINPLT .LT. MAXPLT) THEN GHMINF = MINPLT GHMAXF = MAXPLT ELSE IF (GHMINF .NE. 0.D0 .AND. GHMAXF .NE. 0.D0) 1 THEN GHMINF = GHMINF / 3.5D0 GHMAXF = GHMAXF * 3.5D0 GO TO 600 ELSE GHMINF = .001D0 GHMAXF = 1000.D0 END IF END IF END IF C CALL SAVE (SIZE,SIZE,NSTATS,NSTATS,ASAV,A) CALL SAVE (SIZE,SIZE,NSTATS,NINPS, BSAV,B) CALL SAVE (SIZE,SIZE,NOUTS, NSTATS,CSAV,C) CALL SAVE (SIZE,SIZE,NOUTS, NINPS, DSAV,D) C CALL HFRQPT(SIZE,NSTATS,NINPS,NOUTS, . A,B,C,D,G,ZR,EVRE,EVIM, . AINVB,H,ZC,IPVT,NPTPLT,GHMINF,GHMAXF,ERRCD,ERRMSG) IF (ERRCD .GT. 0) THEN ERRMSG = 'EVLBND Fatal Error: '//ERRMSG CLOSE (UNIT=UNIT8) RETURN END IF C SVX = 'SVG' C CALL SRSPG (SVX,ERRCD,ERRMSG) IF (ERRCD .GT. 0) THEN ERRMSG = 'EVLBND Fatal Error: '//ERRMSG CLOSE (UNIT=UNIT8) RETURN END IF C C--read in the min and max singular values, and associated frequencies, C--from unit 8, which was left open and rewound by srspg C READ (UNIT8,1003,END=9992,ERR=9992) TITLE C READ (UNIT8,1004,END=9992,ERR=9992) WL, WH, NW1 C READ (UNIT8,1006,END=9992,ERR=9992) (FREQ(I),SIGMIN(I), 2 SIGMAX(I),I=1, NPTPLT) C C--this close deletes the temporary file C CLOSE (UNIT=UNIT8) C C--find frequency at which minimum singular value of return ratio C--IS NEAREST ONE C TEMP = 1.D30 C DO 90, I = 1, NPTPLT IF (ABS(SIGMIN(I)-1.D0) .LT. TEMP) THEN TEMP = ABS(SIGMIN(I)-1.D0) MINI = I END IF 90 CONTINUE C TEMP = LOG10(SIGMIN(MINI)) IF (MINI .EQ. 1) THEN C C--make linear approx to sigmin vs. freq, in dB, using next value C--and find zero crossing C SMINSG=20.D0*LOG10(SIGMIN(MINI)/SIGMIN(MINI+1)) 1 / LOG10(FREQ(MINI)/FREQ(MINI+1)) AWCMIN = FREQ(MINI) * 1 10.D0**(-TEMP / SMINSG) ELSE IF (MINI .EQ. NPTPLT) THEN C C--make linear approx to sigmin vs. freq, in dB, using prev. value C--and find zero crossing C SMINSG=20.D0*LOG10(SIGMIN(MINI)/SIGMIN(MINI-1)) 1 / LOG10(FREQ(MINI)/FREQ(MINI-1)) AWCMIN = FREQ(MINI) * 1 10.D0**(-TEMP / SMINSG) ELSE IF (TEMP .GT. 0.D0) THEN C C--make linear approx to sigmin vs. freq, in dB, using next value C--and find zero crossing C SMINSG=20.D0*LOG10(SIGMIN(MINI)/SIGMIN(MINI+1)) 1 / LOG10(FREQ(MINI)/FREQ(MINI+1)) AWCMIN = FREQ(MINI) * 1 10.D0**(-TEMP / SMINSG) ELSE C C--make linear approx to sigmin vs. freq, in dB, using prev. value C--and find zero crossing C SMINSG=20.D0*LOG10(SIGMIN(MINI)/SIGMIN(MINI-1)) 1 / LOG10(FREQ(MINI)/FREQ(MINI-1)) AWCMIN = FREQ(MINI) * 1 10.D0**(-TEMP / SMINSG) END IF C C--find frequency at which maximum singular value of return ratio C--is nearest one C TEMP = 1.D30 C DO 100, I = 1, NPTPLT IF (ABS(SIGMAX(I)-1.D0) .LT. TEMP) THEN TEMP = ABS(SIGMAX(I)-1.D0) MINI = I END IF 100 CONTINUE C TEMP = LOG10(SIGMAX(MINI)) IF (MINI .EQ. 1) THEN C C--make linear approx to sigmax vs. freq, in dB, using next value C--and find zero crossing C SMAXSG=20.D0*LOG10(SIGMAX(MINI)/SIGMAX(MINI+1)) 1 / LOG10(FREQ(MINI)/FREQ(MINI+1)) AWCMAX = FREQ(MINI) * 1 10.D0**(-TEMP / SMAXSG) ELSE IF (MINI .EQ. NPTPLT) THEN C C--make linear approx to sigmax vs. freq, in dB, using prev. value C--and find zero crossing C SMAXSG=20.D0*LOG10(SIGMAX(MINI)/SIGMAX(MINI-1)) 1 / LOG10(FREQ(MINI)/FREQ(MINI-1)) AWCMAX = FREQ(MINI) * 1 10.D0**(-TEMP / SMAXSG) ELSE IF (TEMP .GT. 0.D0) THEN C C--make linear approx to sigmax vs. freq, in dB, using next value C--and find zero crossing C SMAXSG=20.D0*LOG10(SIGMAX(MINI)/SIGMAX(MINI+1)) 1 / LOG10(FREQ(MINI)/FREQ(MINI+1)) AWCMAX = FREQ(MINI) * 1 10.D0**(-TEMP / SMAXSG) ELSE c C--make linear approx to sigmax vs. freq, in dB, using prev. value C--and find zero crossing C SMAXSG=20.D0*LOG10(SIGMAX(MINI)/SIGMAX(MINI-1)) 1 / LOG10(FREQ(MINI)/FREQ(MINI-1)) AWCMAX = FREQ(MINI) * 1 10.D0**(-TEMP / SMAXSG) END IF C C--write results to output_file C OPEN (UNIT=UNIT1,FILE=OUTFIL,ERR=9993) C WRITE (UNIT1,*,ERR=9993) WRCGML,WWRSCL WRITE (UNIT1,*,ERR=9993) WRCGMH,WWRSCH WRITE (UNIT1,*,ERR=9993) AWCMAX,SMAXSG WRITE (UNIT1,*,ERR=9993) AWCMIN,SMINSG C CLOSE (UNIT=UNIT1) C C--QUIT C RETURN C C--I/O ERRORS COME HERE C 9991 ERRCD = 500 ERRMSG = 'EVLBND: Fatal Error accessing or reading from '// 1 'bounds file '//BNDFIL RETURN 9992 ERRCD = 600 ERRMSG = 'EVLBND: Fatal Error reading or writing to '// 1 'scratch unit.' RETURN 9993 ERRCD = 500 ERRMSG = 'EVLBND: Fatal Error opening or writing to '// 1 'output file '//OUTFIL RETURN END