C--------------------- DOUBLE PRECISION VERSION 2B ---------------------MAI00100 C DISCRETE - VERSION 2B (December 1990). FOR THE AUTOMATIC ANALYSIS OFMAI00200 C MULTICOMPONENT EXPONENTIAL DECAY DATA FOR UP TO 9 COMPONENTS. MAI00300 C MUST BE RUN ON A COMPUTER WITH AT LEAST 5 SIGNIFICANT FIGURES IN MAI00400 C ITS SINGLE PRECISION FORTRAN ARITHMETIC AND 9 SIG. FIGS. IN ITS MAI00500 C DOUBLE PRECISION PRECISION FORTRAN ARITHMETIC. MAI00600 C REFERENCES - S. W. PROVENCHER (1976) WRITE UP OF DISCRETE, MAI00700 C (1980) Math. Biosci., Vol. 50, 251-262,MAI00800 C (1976) J. CHEM. PHYS., VOL. 64, 2772-7,MAI00900 C (1976) BIOPHYS. J., VOL. 16, 27-41. MAI01000 C-----------------------------------------------------------------------MAI01100 C CALLS SUBPROGRAMS - BLOCK DATA, DATAIO, WEIGHT, FANLYZ, YANLYZ MAI01200 C WHICH IN TURN CALL - LSTSQR, EVAR, VARF, ETHEOR, PIVOT, PIVOT1, MAI01300 C ANLERR, FISHER, RESIDU, PLPRES MAI01400 C-----------------------------------------------------------------------MAI01500 LOGICAL LAST, REGINT, LBLOKA MAI01600 DOUBLE PRECISION GSE, YLYFIT, DELTAT, E, TSTART, DBLOKA, T MAI01700 COMMON /BLOK/ DBLOKA(5933), RBLOKA(471), IBLOKA(301), NMAX, MAI01800 1 NINTMX, IBLOKB(12), REGINT, LBLOKA(45) MAI01900 C MAI02000 C***********************************************************************MAI02100 C***********************************************************************MAI02200 C MAI02300 C THE FOLLOWING INSTRUCTIONS DESCRIBE ALL POSSIBLE NECESSARY CHANGES MAI02400 C THAT YOU MAY HAVE TO MAKE IN THE MAIN PROGRAM. (SEE ALSO THE MAI02500 C BLOCK DATA SUBPROGRAM.) MAI02600 C MAI02700 C***********************************************************************MAI02800 C IF YOU CHANGE NMAX IN THE BLOCK DATA SUBPROGRAM, YOU MUST READJUST MAI02900 C THE DIMENSIONS IN THE FOLLOWING STATEMENT TO NMAX - MAI03000 C MAI03100 DIMENSION T(200), Y(200), SQRTW(200), YLYFIT(200) MAI03200 C MAI03300 C***********************************************************************MAI03400 C ORDINARILY E SHOULD BE DIMENSIONED E(NMAX,10) IN THE FOLLOWING MAI03500 C STATEMENT. HOWEVER, IF YOU WILL ALWAYS INPUT REGINT=.TRUE., AND MAI03600 C IF NMAX IS GREATER THAN 200, YOU CAN SAVE STORAGE BY DIMENSIONING MAI03700 C IT AS E(1,10) - MAI03800 C MAI03900 DIMENSION E(200,10) MAI04000 C MAI04100 C***********************************************************************MAI04200 C ORDINARILY GSE SHOULD BE DIMENSIONED GSE(500,4). HOWEVER, IF YOU MAI04300 C WILL ALWAYS INPUT REGINT=.FALSE., AND IF NMAX IS LESS THAN 200, MAI04400 C YOU CAN SAVE A LITTLE STORAGE BY DIMENSIONING IT AS GSE(1,4) IN MAI04500 C THE FOLLOWING STATEMENT - MAI04600 C MAI04700 DIMENSION GSE(500,4) MAI04800 C MAI04900 C***********************************************************************MAI05000 C IF YOU CHANGE NINTMX IN THE BLOCK DATA SUBPROGRAM, YOU MUST MAI05100 C READJUST THE DIMENSIONS IN THE FOLLOWING STATEMENT TO NINTMX - MAI05200 C MAI05300 DIMENSION TSTART(10), NT(10), DELTAT(10) MAI05400 C MAI05500 C THIS IS THE END OF ALL POSSIBLE CHANGES YOU MAY HAVE TO MAKE IN MAI05600 C THE MAIN PROGRAM. MAI05700 C MAI05800 C***********************************************************************MAI05900 C***********************************************************************MAI06000 C MAI06100 EQUIVALENCE (E(1,1),GSE(1,1)) MAI06200 C-----------------------------------------------------------------------MAI06300 C DATAIO IS FOR INPUT AND OUTPUT OF INPUT DATA MAI06400 100 CALL DATAIO (T,Y,SQRTW,NMAX,TSTART,DELTAT,NT,NINTMX,LAST) MAI06500 IF (REGINT) GO TO 110 MAI06600 NDIME=NMAX MAI06700 NDIMG=1 MAI06800 GO TO 200 MAI06900 110 NDIME=1 MAI07000 NDIMG=500 MAI07100 C-----------------------------------------------------------------------MAI07200 C WEIGHT GENERATES CRUDE STARTING VALUES FOR CONSTRAINED LEAST SQUARES MAI07300 C FITS TO RAW DATA TO GENERATE WEIGHTS AND DOES INITIAL MAI07400 C COMPUTATIONS OF QUANTITIES FOR LATER USE. MAI07500 200 CALL WEIGHT (.FALSE.,Y,T,SQRTW,YLYFIT,E,NMAX,NDIME, MAI07600 1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG) MAI07700 C-----------------------------------------------------------------------MAI07800 C FANLYZ IS FOR THE CONSTRAINED STEPWISE LEAST SQUARES ANALYSIS OF MAI07900 C THE RAW DATA AND THE TRANSFORMS TO GET STARTING ESTIMATES MAI08000 C FOR THE FINAL LEAST SQUARES ANALYSES OF THE RAW DATA. MAI08100 CALL FANLYZ (T,Y,SQRTW,YLYFIT,E,NMAX,NDIME, MAI08200 1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG) MAI08300 C-----------------------------------------------------------------------MAI08400 C WEIGHT USES STARTING VALUES FROM FANLYZ FOR CONSTRAINED LEAST MAI08500 C SQUARES FITS TO RAW DATA TO GENERATE WEIGHTS FOR THE FINAL MAI08600 C LEAST SQUARES FITS TO THE RAW DATA. MAI08700 CALL WEIGHT (.TRUE.,Y,T,SQRTW,YLYFIT,E,NMAX,NDIME, MAI08800 1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG) MAI08900 C-----------------------------------------------------------------------MAI09000 C YANLYZ IS FOR THE CONSTRAINED STEPWISE LEAST SQUARES ANALYSIS OF MAI09100 C THE RAW DATA (USING THE STARTING VALUES OBTAINED FROM FANLYZ MAI09200 C AND THE WEIGHTS FROM WEIGHT). MAI09300 CALL YANLYZ (T,Y,SQRTW,YLYFIT,E,NMAX,NDIME, MAI09400 1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG) MAI09500 C-----------------------------------------------------------------------MAI09600 IF (.NOT.LAST) GO TO 100 MAI09700 STOP MAI09800 END MAI09900 C--------------------- DOUBLE PRECISION VERSION 2B ---------------------BLK00100 BLOCK DATA BLK00200 LOGICAL REGINT, FAILED, PRPREL, NOBASE, WTED, REPEAT, PLOTRS, BLK00300 1 NONNEG, FFAIL, PIVALP, PRFINL BLK00400 REAL LAMST, LAMF, LAMNMX BLK00500 DOUBLE PRECISION TRBASE, GF, DGFL, DGRIDF, GFSTL, FHAT, VARMIN, BLK00600 1 SUMWTY, SUMWT, ZERO, DGRIDE, GESTL, VAR, PALPHA, PLAM, BLK00700 2 SE, PLMTRY, DELTAP, DFCAPL, SFDE, ADUB, DDSE, DTOLER, F, DF BLK00800 COMMON /BLOK/ TRBASE(19), GF(19,126), DGFL(19,126), DGRIDF, BLK00900 1 GFSTL, FHAT(19), VARMIN, SUMWTY, SUMWT, ZERO(10), BLK01000 2 DGRIDE, GESTL, VAR, PALPHA(10), PLAM(10), SE(10,10), BLK01100 3 PLMTRY(10), DELTAP(19), DFCAPL(10,9), SFDE(9), ADUB(19,19), BLK01200 4 DDSE(9,9), DTOLER(19), F(19,10), DF(19,10), DELTA, FHATMX, BLK01300 5 CONVRG, LAMST(9), LAMF(9,9,2), PRECIS, LAMNMX(2,2), SIGYY BLK01400 COMMON /BLOK/ ENPHI, SIGMAP(19), PCERR(19), VGRID2(252), JGEADD, BLK01500 1 NBINOM(20,9), NVGRID(9), NPERG2(9), MDIMA, IBASE, IWTSGN, BLK01600 2 MINTER, MLAMMX, NGRID2(9), NLAMMX, IWT, LABEL(80), NINT, NF, BLK01700 3 N, IWRITE, IREAD, LINEPG, NMAX, NINTMX, MTRY, ITER, JLAMP1, BLK01800 4 NFTRY(9), REGINT, FAILED, PRPREL, NOBASE, WTED, REPEAT, BLK01900 5 PLOTRS, PRFINL, NONNEG, FFAIL(9,2), PIVALP(19) BLK02000 C BLK02100 C***********************************************************************BLK02200 C***********************************************************************BLK02300 C BLK02400 C THE FOLLOWING INSTRUCTIONS DESCRIBE ALL POSSIBLE NECESSARY CHANGES BLK02500 C THAT YOU MAY HAVE TO MAKE IN THE BLOCKDATA SUBPROGRAM. (SEE THE BLK02600 C MAIN PROGRAM ALSO.) BLK02700 C BLK02800 C***********************************************************************BLK02900 C ALL STATEMENTS FOR READING CARDS AND FOR PRINTING OUTPUT ARE OF BLK03000 C THE FORM READ (IREAD,... AND WRITE(IWRITE,.... THUS, IREAD AND BLK03100 C IWRITE MAY HAVE TO BE CHANGED TO THE VALUES APPROPRIATE FOR YOUR BLK03200 C COMPUTER CENTER IN THE FOLLOWING STATEMENT - BLK03300 C BLK03400 DATA IREAD/5/,IWRITE/6/ BLK03500 C BLK03600 C***********************************************************************BLK03700 C IN THE FOLLOWING STATEMENT, LINEPG MUST BE SET TO THE NUMBER OF BLK03800 C LINES PER PAGE AVAILABLE FOR PRINTED OUTPUT ON YOUR PRINTER - BLK03900 C BLK04000 DATA LINEPG/62/ BLK04100 C BLK04200 C***********************************************************************BLK04300 C NMAX MUST BE GREATER THAN OR EQUAL TO THE MAXIMUM NUMBER OF DATA BLK04400 C POINTS YOU WILL EVER USE. CORE STORAGE CAN BE SAVED BY KEEPING BLK04500 C NMAX AS LITTLE ABOVE THIS MAXIMUM AS POSSIBLE. IF YOU CHANGE NMAX BLK04600 C IN THE FOLLOWING STATEMENT, YOU MUST ALSO RESET THE DIMENSIONS AS BLK04700 C DIRECTED IN THE MAIN PROGRAM - BLK04800 C BLK04900 DATA NMAX/200/ BLK05000 C BLK05100 C***********************************************************************BLK05200 C NINTMX MUST BE GREATER THAN OR EQUAL TO THE MAXIMUM VALUE OF NINT BLK05300 C YOU WILL EVER USE WITH REGINT=.TRUE.. IF YOU CHANGE BLK05400 C NINTMX IN THE FOLLOWING STATEMENT, YOU MUST ALSO RESET THE BLK05500 C DIMENSIONS AS DIRECTED IN THE MAIN PROGRAM - BLK05600 C BLK05700 DATA NINTMX /10/ BLK05800 C BLK05900 C THIS IS THE END OF ALL POSSIBLE NECESSARY CHANGES YOU MAY HAVE TO BLK06000 C MAKE IN THE BLOCK DATA SUBPROGRAM BLK06100 C BLK06200 C***********************************************************************BLK06300 C***********************************************************************BLK06400 C BLK06500 DATA NGRID2/0, 20, 12, 5*10, 12/,MINTER/7/, BLK06600 1 CONVRG/5.E-5/,MLAMMX/9/, SIGMAP/19*0./, PRECIS/1.E-7/ BLK06700 END BLK06800 C--------------------- DOUBLE PRECISION VERSION 2B ---------------------DAT00100 C SUBROUTINE DATAIO. FOR INPUT AND OUTPUT OF INITIAL INPUT DATA DAT00200 C-----------------------------------------------------------------------DAT00300 C CALLS NO EXTERNAL SUBPROGRAMS. DAT00400 C-----------------------------------------------------------------------DAT00500 SUBROUTINE DATAIO (T,Y,SQRTW,NMAX,TSTART,DELTAT,NT,NINTMX,LAST) DAT00600 LOGICAL REGINT, LAST, NONNEG, NOBASE, PRY, PRPREL, PLOTRS, REPEAT,DAT00700 1 LBLOKA, LBLOKB, LBLOKC, ABORT, PRFINL DAT00800 DOUBLE PRECISION TEND, DUB, T, TSTART, DELTAT, DBLOKA DAT00900 COMMON /BLOK/ DBLOKA(5933), RBLOKA(471), IBLOKA(201), IWTSGN, DAT01000 1 IBLOKB, MLAMMX, IBLOKC(9), NLAMMX, IWT, LABEL(80), NINT, DAT01100 2 IBLOKD, N, IWRITE, IREAD, IBLOKE(3), MTRY, IBLOKF(11), DAT01200 3 REGINT, LBLOKA, PRPREL, NOBASE, LBLOKB, REPEAT, PLOTRS, DAT01300 4 PRFINL, NONNEG, LBLOKC(37) DAT01400 DIMENSION T(NMAX), Y(NMAX), SQRTW(NMAX), TSTART(NINTMX), DAT01500 1 NT(NINTMX), DELTAT(NINTMX), LA(9), LATF(2) DAT01600 DATA LATF/1HT, 1HF/ DAT01700 5010 FORMAT (80A1) DAT01800 READ (IREAD,5010) LABEL DAT01900 5020 FORMAT (38H1DISCRETE - Version 2B (December 1990),10X,80A1) DAT02000 WRITE (IWRITE,5020) LABEL DAT02100 ABORT=.FALSE. DAT02200 5100 FORMAT (9A1,3I3) DAT02300 READ (IREAD,5100) LA,NLAMMX,IWT,MTRY DAT02400 DO 101 J=1,9 DAT02500 IF (LA(J).EQ.LATF(1) .OR. LA(J).EQ.LATF(2)) GO TO 101 DAT02600 5101 FORMAT (/43H THE FOLLOWING CARD HAS NO T OR F IN COLUMN,I2) DAT02700 WRITE (IWRITE,5101) J DAT02800 ABORT=.TRUE. DAT02900 101 CONTINUE DAT03000 IF (NLAMMX.GE.1 .AND. NLAMMX.LE.MLAMMX) GO TO 102 DAT03100 5102 FORMAT (/9H NLAMMX =,I4,21H IS NOT BETWEEN 1 AND,I2) DAT03200 WRITE (IWRITE,5102) NLAMMX, MLAMMX DAT03300 ABORT=.TRUE. DAT03400 102 IF (IWT.EQ.1 .OR. IABS(IWT).EQ.2 .OR. IABS(IWT).EQ.3 DAT03500 1 .OR. IWT.EQ.4) GO TO 103 DAT03600 5103 FORMAT (/1X,I4,29H IS AN IMPROPER VALUE FOR IWT) DAT03700 WRITE (IWRITE,5103) IWT DAT03800 ABORT=.TRUE. DAT03900 103 IF (MTRY.GE.1 .AND. MTRY.LE.45) GO TO 104 DAT04000 5104 FORMAT (/7H MTRY =,I4,24H IS NOT BETWEEN 1 AND 45) DAT04100 WRITE (IWRITE,5104) MTRY DAT04200 ABORT=.TRUE. DAT04300 104 IF (.NOT.ABORT) GO TO 105 DAT04400 5105 FORMAT (/21H ERRONEOUS CARD IMAGE,10X,9A1,3I3) DAT04500 WRITE (IWRITE,5105) LA,NLAMMX,IWT,MTRY DAT04600 STOP DAT04700 105 LAST=LA(1) .EQ. LATF(1) DAT04800 REGINT=LA(2) .EQ. LATF(1) DAT04900 NOBASE=LA(3) .EQ. LATF(1) DAT05000 NONNEG=LA(4) .EQ. LATF(1) DAT05100 PRY=LA(5) .EQ. LATF(1) DAT05200 PRPREL=LA(6) .EQ. LATF(1) DAT05300 PRFINL=LA(7) .EQ. LATF(1) DAT05400 PLOTRS=LA(8) .EQ. LATF(1) DAT05500 REPEAT=LA(9) .EQ. LATF(1) DAT05600 5106 FORMAT (/9H LAST =,L2,4X,8HREGINT =,L2,5X,8HNOBASE =,L2,5X, DAT05700 1 8HNONNEG =,L2,5X,5HPRY =,L2,8X,8HPRPREL =,L2,5X,8HPRFINL =,L2,5X,DAT05800 2 8HPLOTRS =,L2,5X,8HREPEAT =,L2/9H NLAMMX =,I2,4X,8HIWT =, DAT05900 3 I2,5X,8HMTRY =,I3) DAT06000 WRITE (IWRITE,5106) LAST,REGINT,NOBASE,NONNEG,PRY,PRPREL,PRFINL, DAT06100 1 PLOTRS,REPEAT,NLAMMX,IWT,MTRY DAT06200 IWTSGN=IWT DAT06300 IWT=IABS(IWT) DAT06400 IF (REGINT) GO TO 150 DAT06500 C-----------------------------------------------------------------------DAT06600 C READS IN T VALUES (REGINT=.FALSE.). DAT06700 C-----------------------------------------------------------------------DAT06800 READ (IREAD,5110) N DAT06900 5110 FORMAT (I5) DAT07000 READ (IREAD,5120) (T(J),J=1,N) DAT07100 5120 FORMAT (5D15.6) DAT07200 NINT=1 DAT07300 TSTART(1)=0.D0 DAT07400 DELTAT(1)=0.D0 DAT07500 NT(1)=N DAT07600 GO TO 200 DAT07700 C-----------------------------------------------------------------------DAT07800 C CALCULATES T VALUES FROM TSTART, TEND, AND NT (REGINT=.TRUE.). DAT07900 C-----------------------------------------------------------------------DAT08000 150 READ (IREAD,5110) NINT DAT08100 IF (NINT.LE.NINTMX .AND. NINT.GE.1) GO TO 160 DAT08200 5160 FORMAT (///7H NINT =,I5,30H IS NOT BETWEEN 1 AND NINTMX =,I2) DAT08300 WRITE (IWRITE,5160) NINT,NINTMX DAT08400 STOP DAT08500 5165 FORMAT (///9X,6HTSTART,11X,4HTEND,3X,2HNT) DAT08600 160 WRITE (IWRITE,5165) DAT08700 N=0 DAT08800 DO 170 J=1,NINT DAT08900 READ (IREAD,5170) TSTART(J),TEND,NT(J) DAT09000 5170 FORMAT (2D15.6,I5) DAT09100 WRITE (IWRITE,5170) TSTART(J),TEND,NT(J) DAT09200 IF (NT(J) .GE. 2) GO TO 175 DAT09300 5175 FORMAT (/24H ABOVE NT IS LESS THAN 2) DAT09400 WRITE (IWRITE,5175) DAT09500 STOP DAT09600 175 DUB=(TEND-TSTART(J))/DBLE(FLOAT(NT(J)-1)) DAT09700 N=N+1 DAT09800 T(N)=TSTART(J) DAT09900 L=NT(J) DAT10000 DELTAT(J)=DUB DAT10100 DO 180 K=2,L DAT10200 N=N+1 DAT10300 180 T(N)=T(N-1)+DUB DAT10400 170 CONTINUE DAT10500 C-----------------------------------------------------------------------DAT10600 C READS IN Y VALUES. DAT10700 C-----------------------------------------------------------------------DAT10800 200 IF (N.GT.2*NLAMMX+3 .AND. N.LE.NMAX) GO TO 205 DAT10900 5205 FORMAT (///4H N =,I5,37H IS NOT BETWEEN 2*NLAMMX+3 AND NMAX =,I5) DAT11000 WRITE (IWRITE,5205) N,NMAX DAT11100 STOP DAT11200 205 IF (PRY) WRITE (IWRITE,5004) DAT11300 5004 FORMAT (///1X) DAT11400 5207 FORMAT (5E15.6) DAT11500 READ (IREAD,5207) (Y(J),J=1,N) DAT11600 IF (IWT .EQ. 4) GO TO 220 DAT11700 IF (.NOT.PRY) GO TO 500 DAT11800 5210 FORMAT (4(19X,1HT,12X,1HY)/) DAT11900 WRITE (IWRITE,5210) DAT12000 5220 FORMAT (4(1PD20.5,E13.5)) DAT12100 WRITE (IWRITE,5220) (T(J),Y(J),J=1,N) DAT12200 GO TO 500 DAT12300 C-----------------------------------------------------------------------DAT12400 C READS IN SPECIAL WEIGHTS (IWT=4). DAT12500 C-----------------------------------------------------------------------DAT12600 220 READ (IREAD,5207) (SQRTW(J),J=1,N) DAT12700 DO 227 J=1,N DAT12800 227 SQRTW(J)=SQRT(SQRTW(J)) DAT12900 IF (.NOT.PRY) GO TO 500 DAT13000 5230 FORMAT (3(17X,1HT,12X,1HY,5X,8HSQRT(WT))/) DAT13100 WRITE (IWRITE,5230) DAT13200 5240 FORMAT (3(1PD18.5,2E13.5)) DAT13300 WRITE (IWRITE,5240) (T(J),Y(J),SQRTW(J),J=1,N) DAT13400 500 RETURN DAT13500 END DAT13600 C--------------------- DOUBLE PRECISION VERSION 2B ---------------------WT000100 C SUBROUTINE WEIGHT. DOES INITIAL LEAST SQUARES FITS TO RAW DATA TO WT000200 C GENERATE WEIGHTS, EVALUATES DATA TRANSFORMS, TRANSFORM AND WT000300 C EXPONENTIAL SUMS AT GRID POINTS FOR INTERPOLATION, AND OTHER WT000400 C QUANTITIES FOR LATER USE. WT000500 C IF FINAL=.TRUE. USES LAMF (VALUES FROM FANLYZ). WT000600 C IF FINAL=.FALSE., GENERATES ITS OWN CRUDE STARTING LAMBDAS. WT000700 C-----------------------------------------------------------------------WT000800 C CALLS SUBPROGRAMS - LSTSQR, ETHEOR WT000900 C WHICH IN TURN CALL - EVAR, PIVOT, PIVOT1, ANLERR WT001000 C-----------------------------------------------------------------------WT001100 SUBROUTINE WEIGHT (FINAL,Y,T,SQRTW,TIMU,E,NMAX,NDIME, WT001200 1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG) WT001300 EXTERNAL EVAR WT001400 LOGICAL FINAL, REGINT, NOBASE, WTED, FAILED, LDUM, LBLOKA, PRPREL,WT001500 1 LBLOKB, FFAIL WT001600 REAL LAMNMX, LAMST, LAMF WT001700 DOUBLE PRECISION VARBES, RGRIDE, GEST, GFST, MU, DUB, DDUB, WT001800 1 DDDUB, RGRIDF, DS, DSS, DSSS, TRBASE, GF, DGFL, DGRIDF, WT001900 2 GFSTL, FHAT, VARMIN, SUMWTY, GSE, SUMWT, ZERO, DGRIDE, GESTL, WT002000 3 VAR, PALPHA, PLAM, TIMU, E, TSTART, DELTAT, DBLOKA, T WT002100 COMMON /BLOK/ TRBASE(19), GF(19,126), DGFL(19,126), DGRIDF, WT002200 1 GFSTL, FHAT(19), VARMIN, SUMWTY, SUMWT, ZERO(10), DGRIDE, WT002300 2 GESTL, VAR, PALPHA(10), PLAM(10), DBLOKA(1069), DELTA, FHATMX, WT002400 3 CONVRG, LAMST(9), LAMF(9,9,2), RBLOKA, LAMNMX(2,2), WT002500 4 RBLOKB(292), JGEADD, NBINOM(20,9), NVGRID(9), NPERG2(9), WT002600 5 MDIMA, IBASE, IWTSGN, MINTER, MLAMMX, NGRID2(9), NLAMMX WT002700 COMMON /BLOK/ IWT, LABEL(80), NINT, NF, N, IWRITE, IBLOKA(16), WT002800 1 REGINT, FAILED, PRPREL, NOBASE, WTED, LBLOKA(4), FFAIL(9,2), WT002900 2 LBLOKB(19) WT003000 DIMENSION E(NDIME,10), WTALP(10), WTLAM(10), MU(19) WT003100 DIMENSION Y(NMAX), T(NMAX), SQRTW(NMAX), TIMU(NMAX), WT003200 1 TSTART(NINTMX), DELTAT(NINTMX), NT(NINTMX), GSE(NDIMG,4), WT003300 2 ALPH(10,2), TNLAMN(2), GAMMA(21) WT003400 DATA ALPH /1HT, 1HR, 1HA, 1HN, 1HS, 1HF, 1HO, 1HR, 1HM, 1HS, WT003500 1 1HR, 1HA, 1HW, 1H , 1HD, 1HA, 1HT, 1HA, 1H , 1H /, WT003600 2 TNLAMN/.2,.02/,T1LAMX/2.08/,NGRID1/240/,NG2MAX/20/ WT003700 IF (FINAL) GO TO 120 WT003800 FHATMX=0. WT003900 MDIMA=2*MLAMMX+1 WT004000 K=MLAMMX+1 WT004100 DO 105 J=1,K WT004200 105 ZERO(J)=0.D0 WT004300 TN=-1.E+20 WT004400 DDUM=-1.E+20 WT004500 DO 106 K=1,5 WT004600 DUM=1.E+20 WT004700 DO 107 J=1,N WT004800 DDDUM=SNGL(T(J)) WT004900 IF (K .EQ. 1) TN=AMAX1(TN,DDDUM) WT005000 107 IF (DDDUM.GT.DDUM .AND. DDDUM.LT.DUM) DUM=DDDUM WT005100 DDUM=DUM+1.E-5*ABS(DUM) WT005200 IF (K .EQ. 1) T1=DUM WT005300 IF (K .EQ. 2) T2=DUM WT005400 106 CONTINUE WT005500 T5=DUM WT005600 TNL1=-1.E+20 WT005700 DUM=TN-1.E-5*ABS(TN) WT005800 DO 108 J=1,N WT005900 108 IF (SNGL(T(J)) .LT. DUM) TNL1=AMAX1(TNL1,SNGL(T(J))) WT006000 DUM=.25*(T5-T1) WT006100 DELTA=DUM-T1 WT006200 C-----------------------------------------------------------------------WT006300 C LAMNMX(J,K) = MIN. (J=1) AND MAX. (J=2) ALLOWED VALUES OF LAMBDA WT006400 C FOR LEAST SQUARES PARAMETERS (K=1) AND FOR STARTING VALUES (K=2).WT006500 C-----------------------------------------------------------------------WT006600 LAMNMX(2,1)=T1LAMX/(T1+DELTA) WT006700 LAMNMX(2,2)=AMIN1(LAMNMX(2,1),.693/(T1+DELTA)) WT006800 LAMNMX(1,2)=TNLAMN(1)/(TN+DELTA) WT006900 IF (NOBASE) GO TO 110 WT007000 NPARMX=2*NLAMMX+1 WT007100 IBASE=1 WT007200 LAMNMX(1,1)=LAMNMX(1,2) WT007300 GO TO 115 WT007400 110 IBASE=0 WT007500 NPARMX=2*NLAMMX WT007600 FHAT(NPARMX+1)=0.D0 WT007700 LAMNMX(1,1)=TNLAMN(2)/(TN+DELTA) WT007800 C-----------------------------------------------------------------------WT007900 C IF NECESSARY, INCREASE DELTA SO THAT DISTANCE BETWEEN FIRST TWO WT008000 C POINTS ON X=ALOG(T+DELTA) AXIS IS LESS THAN A HALF-PERIOD OF WT008100 C THE HIGHEST FREQUENCY USED IN THE TRANSFORMS WT008200 C-----------------------------------------------------------------------WT008300 115 XRANGE=ALOG((TN+DELTA)**2/((T1+DELTA)*(TNL1+DELTA))) WT008400 IF (ALOG((T2+DELTA)/(T1+DELTA))*FLOAT(2*NLAMMX) .LT. XRANGE) WT008500 1 GO TO 150 WT008600 DELTA=DELTA+DUM WT008700 GO TO 115 WT008800 120 IWT=IABS(IWTSGN) WT008900 IF (IWT.EQ.1 .OR. IWT.EQ.4) RETURN WT009000 DO 130 J=1,N WT009100 130 Y(J)=Y(J)/SQRTW(J) WT009200 150 IF (IWT .EQ. 4) GO TO 400 WT009300 SUMWT=DBLE(FLOAT(N)) WT009400 SUMWTY=0.D0 WT009500 DUM=0. WT009600 DO 160 J=1,N WT009700 SUMWTY=SUMWTY+DBLE(Y(J)) WT009800 DUM=DUM+Y(J)*Y(J) WT009900 160 SQRTW(J)=1. WT010000 VARMIN=DBLE(DUM*CONVRG*CONVRG) WT010100 IF (IWT.EQ.1) GO TO 400 WT010200 WTED=.FALSE. WT010300 VARBES=1.D+20 WT010400 RATIOL=LAMNMX(2,2)/LAMNMX(1,2) WT010500 C-----------------------------------------------------------------------WT010600 C START OF MAIN LOOP FOR STEPWISE ANALYSIS WT010700 C-----------------------------------------------------------------------WT010800 DO 200 JLAM=1,NLAMMX WT010900 NF=2*JLAM+IBASE WT011000 5200 FORMAT (1H1,25X,80A1//30X, WT011100 1 41HPRELIMINARY ANALYSIS TO DETERMINE WEIGHTS/30X, WT011200 2 17HANALYSIS ASSUMING,I2,11H COMPONENTS) WT011300 IF (PRPREL) WRITE (IWRITE,5200) LABEL,JLAM WT011400 IF (FINAL) GO TO 220 WT011500 FFAIL(JLAM,1)=.TRUE. WT011600 DUM=RATIOL**(1./FLOAT(JLAM+1)) WT011700 LAMST(1)=LAMNMX(1,2)*DUM WT011800 IF (JLAM .EQ. 1) GO TO 240 WT011900 DO 210 J=2,JLAM WT012000 210 LAMST(J)=LAMST(J-1)*DUM WT012100 GO TO 240 WT012200 220 DO 230 J=1,JLAM WT012300 230 LAMST(J)=LAMF(J,JLAM,1) WT012400 240 CALL LSTSQR (LAMST,JLAM,.TRUE.,N,Y,T,JLAM+IBASE,EVAR,GSE,NDIMG, WT012500 1 SQRTW,TIMU,E,NMAX,TSTART,DELTAT,NT,NINT,1,FFAIL(JLAM,1),NDIME, WT012600 2 .FALSE.,PRPREL) WT012700 IF (FAILED .OR. VAR.GE.VARBES) GO TO 300 WT012800 VARBES=VAR WT012900 NLINWT=JLAM+IBASE WT013000 DO 260 J=1,NLINWT WT013100 WTALP(J)=SNGL(PALPHA(J)) WT013200 260 WTLAM(J)=SNGL(PLAM(J)) WT013300 IF (IWTSGN .LT. 0) NLINWT=JLAM WT013400 200 CONTINUE WT013500 300 IF (VARBES .LT. 1.D+20) GO TO 310 WT013600 5300 FORMAT (////34H 1-COMPONENT ANALYSIS TO DETERMINE, WT013700 1 54H WEIGHTS SOMEHOW FAILED - UNIT WEIGHTS WILL BE USED. , WT013800 2 22(2H**)) WT013900 WRITE (IWRITE,5300) WT014000 IWT=1 WT014100 GO TO 400 WT014200 C-----------------------------------------------------------------------WT014300 C CALCULATE WEIGHTS INCLUDING ERRFIT (STANDARD ERROR IN REGION OF WT014400 C THE THEORETICAL CURVE WHERE IT IS CLOSEST TO Y=0., WHERE A WT014500 C WEIGHT COULD OTHERWISE BE DISASTROUSLY LARGE IF ERRFIT WERE WT014600 C NOT ACCOUNTED FOR). WT014700 C-----------------------------------------------------------------------WT014800 310 DUM=1.E+20 WT014900 DO 320 K=1,N WT015000 S=0. WT015100 DDUM=SNGL(T(K)) WT015200 DO 330 J=1,NLINWT WT015300 330 S=S+WTALP(J)*EXP(-DDUM*WTLAM(J)) WT015400 SQRTW(K)=S WT015500 IF (ABS(S) .GE. DUM) GO TO 320 WT015600 DUM=ABS(S) WT015700 L=K WT015800 320 CONTINUE WT015900 K=MIN0(N,L+5) WT016000 I=MAX0(1,K-9) WT016100 DUM=0. WT016200 DDUM=0. WT016300 IF (IWTSGN.LT.0 .AND. .NOT.NOBASE) DDUM=WTALP(NLINWT+1) WT016400 DO 340 J=I,K WT016500 340 DUM=DUM+(Y(J)-SQRTW(J)-DDUM)**2 WT016600 ERRFIT=SQRT(DUM/FLOAT(K-I+1)) WT016700 DO 350 K=1,N WT016800 SQRTW(K)=1./(ABS(SQRTW(K))+ERRFIT) WT016900 IF (IWT .EQ. 2) SQRTW(K)=SQRT(SQRTW(K)) WT017000 350 CONTINUE WT017100 L=1 WT017200 IF (FINAL) L=2 WT017300 5350 FORMAT (/////41H PARAMETERS USED TO GENERATE WEIGHTS FOR , WT017400 1 10A1//5X,6HLAMBDA,9X,5HALPHA,10X, WT017500 2 49HERRFIT (UNCERTAINTY TERM ADDED TO ABSOLUTE VALUES, WT017600 3 24H OF THEORETICAL CURVE) =,1PE9.2/(1X,1PE10.3,E14.3)) WT017700 WRITE (IWRITE,5350) (ALPH(J,L),J=1,10),ERRFIT,(WTLAM(K),WTALP(K), WT017800 1 K=1,NLINWT) WT017900 400 WTED=IWT.NE.1 WT018000 IF (FINAL) GO TO 600 WT018100 C-----------------------------------------------------------------------WT018200 C INITIAL COMPUTATION OF QUANTITIES TO BE USED LATER WT018300 C-----------------------------------------------------------------------WT018400 DGRIDE=DBLE(ALOG(LAMNMX(2,2)/LAMNMX(1,2))/FLOAT(NGRID1-1)) WT018500 RGRIDE=DEXP(DGRIDE) WT018600 JGEADD=MINTER/2 WT018700 IF (NOBASE) JGEADD=JGEADD+INT(ALOG(LAMNMX(1,2)/LAMNMX(1,1))/ WT018800 1 SNGL(DGRIDE)+.5) WT018900 GEST=DBLE(LAMNMX(1,2)/SNGL(RGRIDE)**(JGEADD+1)) WT019000 GESTL=DLOG(GEST) WT019100 NGE1=NGRID1+JGEADD+INT(ALOG(2.*LAMNMX(2,1)/LAMNMX(2,2))/ WT019200 1 SNGL(DGRIDE)+.5)+MINTER/2 WT019300 IF (NGE1 .LE. 500) GO TO 401 WT019400 5401 FORMAT (////20H FIX LAMNMX. NGE1 =,I9,1X,51(2H**)) WT019500 WRITE (IWRITE,5401) NGE1 WT019600 STOP WT019700 401 JGEADD=JGEADD+1 WT019800 IF (NLAMMX .EQ. 1) GO TO 403 WT019900 DO 402 J=2,NLAMMX WT020000 402 NPERG2(J)=NGRID1/NGRID2(J) WT020100 403 GAMMA(1)=1. WT020200 DO 404 J=1,NG2MAX WT020300 404 GAMMA(J+1)=GAMMA(J)*FLOAT(J) WT020400 DO 405 J=1,NLAMMX WT020500 DO 406 K=J,NG2MAX WT020600 L=K-J+1 WT020700 406 NBINOM(K,J)=INT(GAMMA(K)/(GAMMA(L)*GAMMA(J))+.5) WT020800 IF (J .EQ. 1) GO TO 405 WT020900 L=NGRID2(J)+1 WT021000 K=L-J WT021100 NVGRID(J)=INT(GAMMA(L)/(GAMMA(J+1)*GAMMA(K))+.5) WT021200 405 CONTINUE WT021300 NGF=NGRID1/2+2*(MINTER/2) WT021400 DGRIDF=DBLE(ALOG(LAMNMX(2,1)/LAMNMX(1,1))/FLOAT(NGRID1/2-1)) WT021500 RGRIDF=DEXP(DGRIDF) WT021600 GFST=DBLE(LAMNMX(1,1)/SNGL(RGRIDF)**(MINTER/2+1)) WT021700 GFSTL=DLOG(GFST) WT021800 C-----------------------------------------------------------------------WT021900 C COMPUTE FHAT (TRANSFORMS OF DATA) AND GF AND DGFL (INTERPOLATION WT022000 C GRID POINTS FOR TRANSFORMS AND THEIR DERIVATIVES W.R.T. WT022100 C ALOG(LAMBDA)) WT022200 C-----------------------------------------------------------------------WT022300 IF (.NOT.REGINT) GO TO 420 WT022400 DO 410 J=1,NINT WT022500 410 TSTART(J)=TSTART(J)+DBLE(DELTA) WT022600 420 DUB=DBLE(DELTA) WT022700 DO 430 J=1,N WT022800 430 T(J)=T(J)+DUB WT022900 MU(1)=DBLE(6.283185/XRANGE) WT023000 DO 440 K=1,NPARMX WT023100 IF (K .GT. 1) GO TO 450 WT023200 DO 445 J=1,N WT023300 FHATMX=FHATMX+SQRTW(J)*ABS(Y(J)) WT023400 445 TIMU(J)=DBLE(SQRTW(J)) WT023500 GO TO 470 WT023600 450 L=K/2 WT023700 IF (L+L .LT. K) GO TO 460 WT023800 DUB=MU(1) WT023900 IF (L .GT. 1) DUB=MU(L-1)+DUB WT024000 MU(L)=DUB WT024100 DO 455 J=1,N WT024200 455 TIMU(J)=DCOS(DLOG(T(J))*DUB)*DBLE(SQRTW(J)) WT024300 GO TO 470 WT024400 460 DUB=MU(L) WT024500 DO 465 J=1,N WT024600 465 TIMU(J)=DSIN(DLOG(T(J))*DUB)*DBLE(SQRTW(J)) WT024700 470 DUB=0.D0 WT024800 DDUB=0.D0 WT024900 DO 475 J=1,N WT025000 DUB=DUB+Y(J)*TIMU(J) WT025100 475 DDUB=DDUB+TIMU(J) WT025200 FHAT(K)=DUB WT025300 TRBASE(K)=DDUB WT025400 DUB=GFST WT025500 DO 480 L=1,NGF WT025600 DUB=DUB*RGRIDF WT025700 DS=0.D0 WT025800 DSS=0.D0 WT025900 IF (REGINT) GO TO 500 WT026000 DO 490 J=1,N WT026100 DDUB=DEXP(-DUB*T(J))*TIMU(J) WT026200 DS=DS+DDUB WT026300 490 DSS=DSS-DDUB*T(J) WT026400 GO TO 550 WT026500 C-----------------------------------------------------------------------WT026600 C USE RECURSION RELATIONS IF REGINT=.TRUE. WT026700 C-----------------------------------------------------------------------WT026800 500 J=0 WT026900 DO 510 JJ=1,NINT WT027000 DDUB=DEXP(-DUB*DELTAT(JJ)) WT027100 DDDUB=DEXP(-DUB*TSTART(JJ)) WT027200 II=NT(JJ) WT027300 DO 520 I=1,II WT027400 J=J+1 WT027500 DSSS=DDDUB*TIMU(J) WT027600 DS=DS+DSSS WT027700 DSS=DSS-T(J)*DSSS WT027800 DDDUB=DDDUB*DDUB WT027900 520 CONTINUE WT028000 510 CONTINUE WT028100 550 GF(K,L)=DS WT028200 DGFL(K,L)=DSS*DUB WT028300 480 CONTINUE WT028400 440 CONTINUE WT028500 IF (PRPREL) WRITE (IWRITE,5560) DELTA,FHAT(1),(MU(K),FHAT(2*K), WT028600 1 FHAT(2*K+1),K=1,NLAMMX) WT028700 5560 FORMAT (////6X,2HMU,13X, WT028800 1 38HREAL AND IMAGINARY PARTS OF TRANSFORMS,15X,7HDELTA =,1PE9.2 WT028900 2 /4X,4H.000,1PD17.3/(1X,0PF7.3,1PD17.3,D14.3)) WT029000 DUB=DBLE(DELTA) WT029100 DO 560 J=1,N WT029200 560 T(J)=T(J)-DUB WT029300 IF (.NOT.REGINT) GO TO 600 WT029400 DO 570 J=1,NINT WT029500 570 TSTART(J)=TSTART(J)-DBLE(DELTA) WT029600 600 IF (.NOT.WTED) GO TO 620 WT029700 DUM=0. WT029800 SUMWT=0.D0 WT029900 SUMWTY=0.D0 WT030000 DO 610 J=1,N WT030100 Y(J)=Y(J)*SQRTW(J) WT030200 DUM=DUM+Y(J)*Y(J) WT030300 IF (NOBASE) GO TO 610 WT030400 DUB=DBLE(SQRTW(J)) WT030500 SUMWT=SUMWT+DUB*DUB WT030600 SUMWTY=SUMWTY+DUB*DBLE(Y(J)) WT030700 610 CONTINUE WT030800 VARMIN=DBLE(DUM*CONVRG*CONVRG) WT030900 620 IF (FINAL .OR. .NOT.REGINT) GO TO 900 WT031000 C-----------------------------------------------------------------------WT031100 C COMPUTE GSE (INTERPOLATION GRID POINTS FOR EXPONENTIAL SUMS) WT031200 C-----------------------------------------------------------------------WT031300 DUB=GEST WT031400 L=NGE1-INT(.69315/SNGL(DGRIDE)) WT031500 DO 630 K=1,NGE1 WT031600 DUB=DUB*RGRIDE WT031700 LDUM=K .LE. L WT031800 CALL ETHEOR (DUB,.TRUE.,WTED,T,SQRTW,Y,NMAX,TSTART,DELTAT,NT, WT031900 1 NINT,DS,DSS,DSSS,DDDUB,.TRUE.,LDUM) WT032000 GSE(K,1)=DS WT032100 GSE(K,2)=DSS WT032200 GSE(K,3)=DSSS WT032300 GSE(K,4)=DDDUB WT032400 630 CONTINUE WT032500 900 RETURN WT032600 END WT032700 C--------------------- DOUBLE PRECISION VERSION 2B ---------------------FAN00100 C SUBROUTINE FANLYZ. DOES CONSTRAINED STEPWISE LEAST SQUARES ANALYSIS FAN00200 C OF RAW DATA (USING GRID SEARCH IF NECESSARY) TO GET STARTING FAN00300 C VALUES FOR ANALYSIS OF TRANSFORMS, WHICH YEILD LAMF (STARTING FAN00400 C VALUES FOR FINAL ANALYSIS OF RAW DATA IN YANLYZ). FAN00500 C-----------------------------------------------------------------------FAN00600 C CALLS SUBPROGRAMS - LSTSQR, EVAR FAN00700 C WHICH IN TURN CALL - PIVOT, PIVOT1, ETHEOR, VARF FAN00800 C-----------------------------------------------------------------------FAN00900 SUBROUTINE FANLYZ (T,Y,SQRTW,YLYFIT,E,NMAX,NDIME, FAN01000 1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG) FAN01100 EXTERNAL VARF, EVAR FAN01200 LOGICAL PRPREL, WTED, FAILED, SEARCH, LDUM, FFAIL, LBLOKA, FAN01300 1 LBLOKB, LBLOKC, LBLOKD, SAVEPK FAN01400 REAL LAMNMX, LAMST, LAMF FAN01500 DOUBLE PRECISION YLYFIT, E, TSTART, DELTAT, GSE, DUB, DDUB, T, FAN01600 1 DGRIDE, GESTL, VAR, PALPHA, PLAM, VARE, VARTRY, DBLOKA, DBLOKB FAN01700 COMMON /BLOK/ DBLOKA(4841), DGRIDE, GESTL, VAR, PALPHA(10), FAN01800 1 PLAM(10), DBLOKB(1069), RBLOKA(3), LAMST(9), LAMF(9,9,2), FAN01900 2 RBLOKB, LAMNMX(2,2), RBLOKC(40), VGRID2(252), JGEADD, FAN02000 3 IBLOKA(180), NVGRID(9), NPERG2(9), IBLOKB, IBASE, IBLOKC(2), FAN02100 4 MLAMMX, NGRID2(9), NLAMMX, IWT, LABEL(80), NINT, NF, N, FAN02200 5 IWRITE, IBLOKD(4), MTRY, IBLOKE(2), NFTRY(9), LBLOKA, FAILED FAN02300 COMMON /BLOK/ PRPREL, LBLOKB, WTED, LBLOKC(4), FFAIL(9,2), FAN02400 1 LBLOKD(19) FAN02500 DIMENSION E(NDIME,10), LENDEQ(9), LEQ(9), JCENTR(20), LSTART(9,45)FAN02600 DIMENSION T(NMAX), Y(NMAX), SQRTW(NMAX), GSE(NDIMG,4), FAN02700 1 YLYFIT(NMAX), TSTART(NINTMX), DELTAT(NINTMX), NT(NINTMX) FAN02800 EQUIVALENCE (LENDEQ(1),LEND1), (LENDEQ(2),LEND2), (LENDEQ(3), FAN02900 1 LEND3), (LENDEQ(4),LEND4), (LENDEQ(5),LEND5), (LENDEQ(6),LEND6), FAN03000 2 (LENDEQ(7),LEND7), (LENDEQ(8),LEND8), (LENDEQ(9),LEND9), FAN03100 3 (LEQ(1),L1), (LEQ(2),L2), (LEQ(3),L3), (LEQ(4),L4), (LEQ(5),L5), FAN03200 4 (LEQ(6),L6), (LEQ(7),L7), (LEQ(8),L8), (LEQ(9),L9) FAN03300 DATA LLPEAK/1/ FAN03400 WTED=IWT .NE. 1 FAN03500 LAMST(1)=SQRT(LAMNMX(2,2)*LAMNMX(1,2)) FAN03600 C-----------------------------------------------------------------------FAN03700 C START OF MAIN LOOP FOR STEPWISE ANALYSIS FAN03800 C-----------------------------------------------------------------------FAN03900 DO 200 JLAM=1,NLAMMX FAN04000 C-----------------------------------------------------------------------FAN04100 C INITIALIZE FOR FITS ASSUMING JLAM COMPONENTS FAN04200 C-----------------------------------------------------------------------FAN04300 5150 FORMAT (1H1,25X,80A1/) FAN04400 IF (PRPREL) WRITE (IWRITE,5150) LABEL FAN04500 VARE=1.D+20 FAN04600 VARTRY=1.D+20 FAN04700 NF=2*JLAM+IBASE FAN04800 IF (JLAM .EQ. 1) GO TO 220 FAN04900 LVGRID=NVGRID(JLAM) FAN05000 DO 210 J=1,LVGRID FAN05100 210 VGRID2(J)=0. FAN05200 220 JLAML1=JLAM-1 FAN05300 JLAMP1=JLAM+1 FAN05400 IF (PRPREL) WRITE (IWRITE,5210) JLAM FAN05500 5210 FORMAT (30X,45HINITIAL ANALYSIS OF THE RAW DATA AND THEN THE, FAN05600 2 22H TRANSFORMS - ASSUMING,I2,11H COMPONENTS) FAN05700 C-----------------------------------------------------------------------FAN05800 C START OF LOOP FOR UP TO MTRY JLAM-COMPONENT FITS FAN05900 C-----------------------------------------------------------------------FAN06000 DO 300 JTRY=1,MTRY FAN06100 5310 FORMAT (////30X,22HFIT USING STARTING SET,I3) FAN06200 IF (PRPREL) WRITE (IWRITE,5310) JTRY FAN06300 IF (JLAM.EQ.1) GO TO 600 FAN06400 IF (JTRY-2) 315,350,500 FAN06500 C-----------------------------------------------------------------------FAN06600 C SELECT FIRST SET OF STARTING VALUES FAN06700 C-----------------------------------------------------------------------FAN06800 315 DO 320 J=1,JLAML1 FAN06900 320 LAMST(J)=LAMF(J,JLAML1,2) FAN07000 IF (LPEAK.EQ.1 .OR. LPEAK.EQ.JLAM) GO TO 345 FAN07100 LAMST(JLAM)=SQRT(LAMF(LPEAK-1,JLAML1,2)*LAMF(LPEAK,JLAML1,2)) FAN07200 GO TO 600 FAN07300 345 IF (LPEAK .EQ. 1) LAMST(JLAM)=SQRT(LAMNMX(1,2)*LAMF(1,JLAML1,2)) FAN07400 IF (LPEAK .EQ. JLAM) LAMST(JLAM)=SQRT(LAMF(JLAML1,JLAML1,2)* FAN07500 1 LAMNMX(2,2)) FAN07600 GO TO 600 FAN07700 C-----------------------------------------------------------------------FAN07800 C DO COARSE GRID SEARCH FAN07900 C-----------------------------------------------------------------------FAN08000 350 L=NGRID2(JLAM)-JLAM FAN08100 DO 360 J=1,JLAM FAN08200 L=L+1 FAN08300 360 LENDEQ(J)=L FAN08400 IF (JLAM .EQ. MLAMMX) GO TO 368 FAN08500 DO 365 J=JLAMP1,MLAMMX FAN08600 365 LENDEQ(J)=LENDEQ(JLAM) FAN08700 368 JCENTR(1)=JGEADD-1+(NPERG2(JLAM)+1)/2 FAN08800 L=NGRID2(JLAM) FAN08900 DO 370 J=2,L FAN09000 370 JCENTR(J)=JCENTR(J-1)+NPERG2(JLAM) FAN09100 SEARCH=.TRUE. FAN09200 C-----------------------------------------------------------------------FAN09300 C ENTRY POINT TO FIND STARTING LAMBDA VALUES FROM GRID FAN09400 C-----------------------------------------------------------------------FAN09500 400 LOC=0 FAN09600 DO 410 L1=1,LEND1 FAN09700 LST2=L1+1 FAN09800 DO 420 L2=LST2,LEND2 FAN09900 LST3=L2+1 FAN10000 IF (JLAM .LT. 3) LST3=LEND3 FAN10100 DO 430 L3=LST3,LEND3 FAN10200 LST4=L3+1 FAN10300 IF (JLAM .LT. 4) LST4=LEND4 FAN10400 DO 440 L4=LST4,LEND4 FAN10500 LST5=L4+1 FAN10600 IF (JLAM .LT. 5) LST5=LEND5 FAN10700 DO 450 L5=LST5,LEND5 FAN10800 LST6=L5+1 FAN10900 IF (JLAM .LT. 6) LST6=LEND6 FAN11000 DO 460 L6=LST6,LEND6 FAN11100 LST7=L6+1 FAN11200 IF (JLAM .LT. 7) LST7=LEND7 FAN11300 DO 470 L7=LST7,LEND7 FAN11400 LST8=L7+1 FAN11500 IF (JLAM .LT. 8) LST8=LEND8 FAN11600 DO 480 L8=LST8,LEND8 FAN11700 LST9=L8+1 FAN11800 IF (JLAM .LT. 9) LST9=LEND9 FAN11900 DO 490 L9=LST9,LEND9 FAN12000 LOC=LOC+1 FAN12100 IF (SEARCH) GO TO 498 FAN12200 IF (VGRID2(LOC) .GT. VGMAX) GO TO 496 FAN12300 MINDIF=1000 FAN12400 IF (JTRY .LE. 2) GO TO 493 FAN12500 L=JTRY-1 FAN12600 DO 491 J=2,L FAN12700 I=0 FAN12800 DO 492 K=1,JLAM FAN12900 492 I=I+IABS(LEQ(K)-LSTART(K,J)) FAN13000 MINDIF=MIN0(MINDIF,I) FAN13100 491 CONTINUE FAN13200 IF (MINDIF-MAXDIF) 496,493,494 FAN13300 493 IF (VGRID2(LOC) .GE. VGBEST) GO TO 496 FAN13400 494 VGBEST=VGRID2(LOC) FAN13500 MAXDIF=MINDIF FAN13600 DO 495 J=1,JLAM FAN13700 495 LSTART(J,JTRY)=LEQ(J) FAN13800 496 IF (LOC .LT. LVGRID) GO TO 490 FAN13900 DO 497 J=1,JLAM FAN14000 L=LSTART(J,JTRY) FAN14100 497 LAMST(J)=EXP(SNGL(GESTL+FLOAT(JCENTR(L))*DGRIDE)) FAN14200 GO TO 600 FAN14300 498 DO 499 J=1,JLAM FAN14400 L=LEQ(J) FAN14500 499 PLAM(J)=DEXP(GESTL+FLOAT(JCENTR(L))*DGRIDE) FAN14600 CALL EVAR (0.,JLAM,DUB,T,SQRTW,NMAX,TSTART,DELTAT,NT, FAN14700 1 NINT,.FALSE.,LDUM,E,Y,YLYFIT,JLAM+IBASE,1.D+20,0,NDIME,.TRUE., FAN14800 2 GSE,NDIMG,J) FAN14900 490 CONTINUE FAN15000 480 CONTINUE FAN15100 470 CONTINUE FAN15200 460 CONTINUE FAN15300 450 CONTINUE FAN15400 440 CONTINUE FAN15500 430 CONTINUE FAN15600 420 CONTINUE FAN15700 410 CONTINUE FAN15800 C-----------------------------------------------------------------------FAN15900 C FIND LOWEST VARIANCE ON COARSE GRID FAN16000 C-----------------------------------------------------------------------FAN16100 500 VGMAX=1.E+20 FAN16200 DO 510 J=1,LVGRID FAN16300 510 VGMAX=AMIN1(VGMAX,VGRID2(J)) FAN16400 IF (JTRY .GT. 2) VGMAX=10.*VGMAX FAN16500 MAXDIF=-1 FAN16600 SEARCH=.FALSE. FAN16700 VGBEST=1.E+20 FAN16800 GO TO 400 FAN16900 C-----------------------------------------------------------------------FAN17000 C DO LEAST SQUARES FITS FAN17100 C-----------------------------------------------------------------------FAN17200 600 LDUM=JTRY.EQ.1 FAN17300 CALL LSTSQR (LAMST,JLAM,.TRUE.,N,Y,T,JLAM+IBASE, FAN17400 1 EVAR,GSE,NDIMG,SQRTW,YLYFIT,E,NMAX,TSTART,DELTAT,NT,NINT, FAN17500 2 2, LDUM ,NDIME,.TRUE.,PRPREL) FAN17600 FAILED=FAILED .OR. VAR.GT.VARE FAN17700 IF (VAR .GT. VARE) GO TO 300 FAN17800 VARE=VAR FAN17900 DO 605 J=1,JLAM FAN18000 LAMST(J)=SNGL(PLAM(J)) FAN18100 605 LAMF(J,JLAM,1)=LAMST(J) FAN18200 SAVEPK=VARTRY .GE. 1.D+20 FAN18300 IF (FAILED) GO TO 610 FAN18400 5600 FORMAT (/51H END OF FIT TO RAW DATA, START OF FIT TO TRANSFORMS) FAN18500 IF (PRPREL) WRITE (IWRITE,5600) FAN18600 CALL LSTSQR (LAMST,JLAM,.FALSE.,NF,Y,T,JLAM+IBASE, FAN18700 1 VARF,GSE,NDIMG,SQRTW,YLYFIT,E,NMAX,TSTART,DELTAT,NT,NINT, FAN18800 2 3,.FALSE.,NDIME,.TRUE.,PRPREL) FAN18900 IF (FAILED .AND. VAR.GE.VARTRY) GO TO 300 FAN19000 SAVEPK=.TRUE. FAN19100 VARTRY=VAR FAN19200 DO 608 J=1,JLAM FAN19300 608 LAMF(J,JLAM,2)=SNGL(PLAM(J)) FAN19400 610 IF (.NOT.SAVEPK .AND. JLAM.NE.1) GO TO 300 FAN19500 IF (JLAM .EQ. NLAMMX) GO TO 670 FAN19600 C-----------------------------------------------------------------------FAN19700 C FINDS MAX SUM OF ABSOLUTE ADJACENT AMPLITUDES FOR LATER DETERMINATIONFAN19800 C OF STARTING VALUE FOR LAMST(JLAMP1) FAN19900 C-----------------------------------------------------------------------FAN20000 DUB=DABS(PALPHA(1))+DABS(PALPHA(JLAMP1)) FAN20100 LLPEAK=1 FAN20200 IF (DUB .GT. DABS(PALPHA(JLAM))) GO TO 620 FAN20300 LLPEAK=JLAMP1 FAN20400 DUB=DABS(PALPHA(JLAM)) FAN20500 620 IF (JLAM .EQ. 1) GO TO 690 FAN20600 DO 630 J=2,JLAM FAN20700 DDUB=DABS(PALPHA(J-1))+DABS(PALPHA(J)) FAN20800 IF (DDUB .LE. DUB) GO TO 630 FAN20900 DUB=DDUB FAN21000 LLPEAK=J FAN21100 630 CONTINUE FAN21200 670 IF (.NOT.FAILED .OR. JLAM.EQ.1) GO TO 690 FAN21300 300 CONTINUE FAN21400 JTRY=MTRY FAN21500 690 LPEAK=LLPEAK FAN21600 NFTRY(JLAM)=JTRY FAN21700 FFAIL(JLAM,2)=FAILED FAN21800 FFAIL(JLAM,1)=VARTRY .GE. 1.D+20 FAN21900 IF (.NOT.FFAIL(JLAM,1)) GO TO 200 FAN22000 DO 695 J=1,JLAM FAN22100 695 LAMF(J,JLAM,2)=LAMF(J,JLAM,1) FAN22200 200 CONTINUE FAN22300 RETURN FAN22400 END FAN22500 C--------------------- DOUBLE PRECISION VERSION 2B ---------------------YAN00100 C SUBROUTINE YANLYZ. FOR CONSTRAINED STEPWISE LEAST SQUARES ANALYSIS YAN00200 C OF RAW DATA USING LAMF (STARTING VALUES OBTAINED FROM FANLYZ). YAN00300 C-----------------------------------------------------------------------YAN00400 C CALLS SUBPROGRAMS - LSTSQR, FISHER, RESIDU YAN00500 C WHICH IN TURN CALL - EVAR, ETHEOR, PIVOT, PIVOT1, PLPRES YAN00600 C-----------------------------------------------------------------------YAN00700 SUBROUTINE YANLYZ (T,Y,SQRTW,YLYFIT,E,NMAX,NDIME, YAN00800 1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG) YAN00900 EXTERNAL EVAR YAN01000 LOGICAL NOBASE, PRFINL, PLOTRS, REPEAT, WTED, FAILED, YFAIL, DONE,YAN01100 1 REGINT, FFAIL, LBLOKA, LBLOKB, LBLOKC YAN01200 REAL LAMNMX, LAMST, LAMF YAN01300 DOUBLE PRECISION VARBES, VARMIN, VAR, PALPHA, PLAM, YLYFIT, E, YAN01400 1 TSTART, DELTAT, GSE, DBLOKA, DBLOKB, T YAN01500 COMMON /BLOK/ DBLOKA(4843), VAR, PALPHA(10), PLAM(10), YAN01600 1 DBLOKB(1069), RBLOKA(3), LAMST(9), LAMF(9,9,2), RBLOKB, YAN01700 2 LAMNMX(2,2), SIGYY, ENPHI, SIGMAP(19), PCERR(19), RBLOKC(252), YAN01800 3 IBLOKA(200), IBASE, IBLOKB(12), NLAMMX, IWT, LABEL(80), NINT, YAN01900 4 NF, N, IWRITE, IBLOKC, LINEPG, IBLOKD(3), ITER, JLAMP1, YAN02000 5 NFTRY(9), REGINT, FAILED, LBLOKA, NOBASE, WTED, REPEAT YAN02100 COMMON /BLOK/ PLOTRS, PRFINL, LBLOKB, FFAIL(9,2), LBLOKC(19) YAN02200 DIMENSION E(NDIME,10), YFAIL(9), ASAVE(10,9,6), YAN02300 1 PROBSV(9), NEXBES(6,5), VARSV(9), SIGYSV(9), ITERSV(9), YAN02400 2 ENPHSV(9), PROBLN(9), DONE(9) YAN02500 DIMENSION T(NMAX), Y(NMAX), SQRTW(NMAX), GSE(NDIMG,4), YAN02600 1 YLYFIT(NMAX), TSTART(NINTMX), DELTAT(NINTMX), YAN02700 2 NT(NINTMX), PUNCOR(6) YAN02800 DATA NEXBES/1H , 1H , 1H , 1H , 1H , 1H , 1HS, 1HE, 1HC, 1HO, YAN02900 1 1HN, 1HD, 1H , 1HT, 1HH, 1HI, 1HR, 1HD, 1HF, 1HO, 1HU, 1HR, 1HT, YAN03000 2 1HH, 1H , 1HF, 1HI, 1HF, 1HT, 1HH/ YAN03100 WTED=IWT .NE. 1 YAN03200 JLAMBS=1 YAN03300 VARBES=1.D+20 YAN03400 VARMIN=1.D+20 YAN03500 ABSYMX=0. YAN03600 DO 180 J=1,N YAN03700 180 ABSYMX=AMAX1(ABSYMX,ABS(Y(J))) YAN03800 C-----------------------------------------------------------------------YAN03900 C START OF MAIN LOOP FOR STEPWISE ANALYSIS YAN04000 C-----------------------------------------------------------------------YAN04100 DO 200 JLAM=1,NLAMMX YAN04200 JLAMP1=JLAM+1 YAN04300 NF=2*JLAM+IBASE YAN04400 IF (.NOT.PRFINL) GO TO 210 YAN04500 5200 FORMAT (1H1,25X,80A1/) YAN04600 WRITE (IWRITE,5200) LABEL YAN04700 5240 FORMAT (40X,23HFINAL ANALYSIS ASSUMING,I2,11H COMPONENTS) YAN04800 WRITE (IWRITE,5240) JLAM YAN04900 210 DO 220 J=1,JLAM YAN05000 220 LAMST(J)=LAMF(J,JLAM,2) YAN05100 CALL LSTSQR (LAMST,JLAM,.TRUE.,N,Y,T,JLAM+IBASE,EVAR, YAN05200 1 GSE,NDIMG,SQRTW,YLYFIT,E,NMAX,TSTART,DELTAT,NT,NINT, YAN05300 2 4,FFAIL(JLAM,2),NDIME,.FALSE.,PRFINL) YAN05400 YFAIL(JLAM)=FAILED YAN05500 IF (.NOT.FAILED) GO TO 225 YAN05600 IF (FFAIL(JLAM,1)) GO TO 250 YAN05700 5222 FORMAT (////43H SECOND ATTEMPT USING STARTING LAMBDAS FROM, YAN05800 1 29H INITIAL ANALYSIS OF RAW DATA) YAN05900 IF (PRFINL) WRITE (IWRITE,5222) YAN06000 DO 222 J=1,JLAM YAN06100 222 LAMST(J)=LAMF(J,JLAM,1) YAN06200 CALL LSTSQR (LAMST,JLAM,.TRUE.,N,Y,T,JLAM+IBASE,EVAR, YAN06300 1 GSE,NDIMG,SQRTW,YLYFIT,E,NMAX,TSTART,DELTAT,NT,NINT, YAN06400 2 4,.FALSE.,NDIME,.FALSE.,PRFINL) YAN06500 YFAIL(JLAM)=FAILED YAN06600 IF (FAILED) GO TO 250 YAN06700 C-----------------------------------------------------------------------YAN06800 C SAVE VALUES FOR LATER USE AND FOR FINAL OUTPUT YAN06900 C-----------------------------------------------------------------------YAN07000 225 SIGYSV(JLAM)=SIGYY YAN07100 VARSV(JLAM)=SNGL(VAR) YAN07200 ITERSV(JLAM)=ITER YAN07300 ENPHSV(JLAM)=ENPHI YAN07400 I=JLAM YAN07500 DO 230 K=1,JLAM YAN07600 ASAVE(K,JLAM,1)=SNGL(PALPHA(K)) YAN07700 ASAVE(K,JLAM,2)=SNGL(PLAM(K)) YAN07800 ASAVE(K,JLAM,3)=SIGMAP(K) YAN07900 I=I+1 YAN08000 ASAVE(K,JLAM,4)=SIGMAP(I) YAN08100 ASAVE(K,JLAM,5)=PCERR(K) YAN08200 ASAVE(K,JLAM,6)=PCERR(I) YAN08300 230 CONTINUE YAN08400 IF (NOBASE) GO TO 240 YAN08500 ASAVE(JLAMP1,JLAM,2)=0. YAN08600 ASAVE(JLAMP1,JLAM,3)=0. YAN08700 ASAVE(JLAMP1,JLAM,5)=0. YAN08800 ASAVE(JLAMP1,JLAM,1)=SNGL(PALPHA(JLAMP1)) YAN08900 ASAVE(JLAMP1,JLAM,4)=SIGMAP(NF) YAN09000 ASAVE(JLAMP1,JLAM,6)=PCERR(NF) YAN09100 C-----------------------------------------------------------------------YAN09200 C TEST IF THIS SOLUTION IS THE BEST SO FAR YAN09300 C-----------------------------------------------------------------------YAN09400 240 IF (VARBES .GE. 1.D+20) GO TO 248 YAN09500 IF (VAR .GE. VARMIN) GO TO 200 YAN09600 K=N-NF YAN09700 DUM=SNGL(VARBES/VAR-1.D0) YAN09800 L=2*(JLAM-JLAMBS) YAN09900 IF (FISHER(FLOAT(K)*DUM/(FLOAT(L)*(1.+ENPHI*FLOAT((NF+2)*N)/ YAN10000 1 FLOAT(NF*K))),L,K) .LE. .5) GO TO 250 YAN10100 248 JLAMBS=JLAM YAN10200 VARBES=VAR YAN10300 250 VARMIN=DMIN1(VARMIN,VAR) YAN10400 200 CONTINUE YAN10500 C-----------------------------------------------------------------------YAN10600 C END OF MAIN LOOP FOR STEPWISE ANALYSIS YAN10700 C-----------------------------------------------------------------------YAN10800 IF (VARBES .LT. 1.D+20) GO TO 300 YAN10900 5300 FORMAT (///50H ALL ANALYSES OF THE RAW DATA FAILED - CHECK INPUT, YAN11000 1 22H DATA FOR GROSS ERRORS ) YAN11100 WRITE (IWRITE,5300) YAN11200 RETURN YAN11300 C-----------------------------------------------------------------------YAN11400 C CALCULATE PROBSV AND PROBLN (PROBABILITY AND UNCORRECTED PROBABILITY YAN11500 C (I.E., USING ENPHI=0.) THAT THE BEST SOLUTION IS ACTUALLY YAN11600 C BETTER THAN THE OTHERS) YAN11700 C-----------------------------------------------------------------------YAN11800 300 PNXBES=1.E+20 YAN11900 JLAMNX=JLAMBS YAN12000 DO 310 J=1,NLAMMX YAN12100 PROBSV(J)=2. YAN12200 IF (J.EQ.JLAMBS .OR. YFAIL(J)) GO TO 310 YAN12300 I=2*MAX0(JLAMBS,J)+IBASE YAN12400 K=N-I YAN12500 L=2*(JLAMBS-J) YAN12600 DUM=VARSV(J)/SNGL(VARBES) YAN12700 IF (L .LT. 0) GO TO 315 YAN12800 DDUM=FLOAT(K)*(DUM-1.)/FLOAT(L) YAN12900 DDDUM=1.+ENPHSV(JLAMBS)*FLOAT((I+2)*N)/FLOAT(I*K) YAN13000 PROBSV(J)=FISHER(DDUM/DDDUM,L,K) YAN13100 PROBLN(J)=FISHER(DDUM,L,K) YAN13200 GO TO 325 YAN13300 315 IF (DUM .LT. 1.) GO TO 320 YAN13400 PROBSV(J)=1. YAN13500 PROBLN(J)=1. YAN13600 GO TO 325 YAN13700 320 DDUM=FLOAT(K)*(1.-1./DUM)/FLOAT(L) YAN13800 DDDUM=1.+ENPHSV(J)*FLOAT((I+2)*N)/FLOAT(I*K) YAN13900 PROBLN(J)=1.-FISHER(DDUM,-L,K) YAN14000 PROBSV(J)=1.-FISHER(DDUM/DDDUM,-L,K) YAN14100 325 IF (PROBSV(J) .GE. PNXBES) GO TO 310 YAN14200 PNXBES=PROBSV(J) YAN14300 PNXLIN=PROBLN(J) YAN14400 JLAMNX=J YAN14500 310 CONTINUE YAN14600 PROBSV(JLAMBS)=0. YAN14700 IF (PNXBES .LT. 1.E+20) GO TO 330 YAN14800 PNXBES=1. YAN14900 PNXLIN=1. YAN15000 C-----------------------------------------------------------------------YAN15100 C FINAL SUMMARY OF RESULTS - SUMMARIZES UP TO BEST 5 SOLUTIONS. YAN15200 C-----------------------------------------------------------------------YAN15300 330 CALL RESIDU (JLAMBS,PLOTRS,T,Y,SQRTW,YLYFIT,NMAX,PUNCOR,ASAVE) YAN15400 5330 FORMAT (38H1DISCRETE - Version 2B (December 1990),10X,80A1/) YAN15500 340 WRITE (IWRITE,5330) LABEL YAN15600 DO 350 J=1,NLAMMX YAN15700 350 DONE(J)=YFAIL(J) YAN15800 L=17+JLAMBS+IBASE YAN15900 5350 FORMAT (30X,6A1,14H BEST SOLUTION,I13,11H COMPONENTS) YAN16000 WRITE (IWRITE,5350) (NEXBES(J,1),J=1,6),JLAMBS YAN16100 IF (PNXBES .GT. .95) GO TO 360 YAN16200 5360 FORMAT (48H ++++++ NOTE - THE APPROX. PROBABILITY THAT THIS, YAN16300 1 34H SOLUTION IS ACTUALLY BEST IS ONLY,F6.3, YAN16400 2 38H. LOOK AT SECOND BEST SOLUTION ALSO. ,6(1H+)) YAN16500 WRITE (IWRITE,5360) PNXBES YAN16600 L=L-1 YAN16700 360 DO 370 J=1,NLAMMX YAN16800 IF (J .GT. 1) GO TO 375 YAN16900 DONE(JLAMBS)=.TRUE. YAN17000 I=JLAMBS YAN17100 GO TO 400 YAN17200 375 DUM=1.E+20 YAN17300 DO 380 K=1,NLAMMX YAN17400 IF (PROBSV(K).GT.DUM .OR. DONE(K)) GO TO 380 YAN17500 DUM=PROBSV(K) YAN17600 I=K YAN17700 380 CONTINUE YAN17800 IF (DUM.GE.1.E+20 .OR. J.GT.5) GO TO 500 YAN17900 DONE(I)=.TRUE. YAN18000 CALL RESIDU (I,.FALSE.,T,Y,SQRTW,YLYFIT,NMAX,PUNCOR,ASAVE) YAN18100 L=L+15+I+IBASE YAN18200 IF (L-2 .LE. LINEPG) GO TO 385 YAN18300 L=14+I+IBASE YAN18400 WRITE (IWRITE,5999) YAN18500 5999 FORMAT (1H1) YAN18600 GO TO 390 YAN18700 385 WRITE (IWRITE,5001) YAN18800 WRITE (IWRITE,5001) YAN18900 5001 FORMAT (1X) YAN19000 390 IF (DUM .GT. .95) WRITE (IWRITE,5350) (NEXBES(K,J),K=1,6),I YAN19100 IF (DUM .LE. .95) WRITE (IWRITE,5385) (NEXBES(K,J),K=1,6),I YAN19200 5385 FORMAT (30X,6A1,14H BEST SOLUTION,I13,11H COMPONENTS,3X, YAN19300 1 30H- A SIGNIFICANT POSSIBILITY ,25(1H+)) YAN19400 400 WRITE (IWRITE,5001) YAN19500 WRITE (IWRITE,5400) YAN19600 5400 FORMAT (1X,70(1H*)) YAN19700 WRITE (IWRITE,5402) NFTRY(I) YAN19800 5402 FORMAT (2H *,5X,25HALPHA +- STD ERR PERCENT,11X, YAN19900 1 31HLAMBDA +- STD ERR PERCENT * , YAN20000 2 41HSTARTING LAMBDA (FROM FIT TO TRANSFORMS -, YAN20100 3 I4,7H TRIES)) YAN20200 5404 FORMAT (2H *,1PE10.3,3H +-,E8.1,0PF9.3,1PE17.3,3H +-, YAN20300 1 E8.1,0PF9.3,2H *,1PE18.3, YAN20400 2 39H (NO EXACT FIT TO THE TRANSFORMS FOUND)) YAN20500 IF (FFAIL(I,2)) WRITE (IWRITE,5404) (ASAVE(K,I,1),ASAVE(K,I,4), YAN20600 1 ASAVE(K,I,6),ASAVE(K,I,2),ASAVE(K,I,3),ASAVE(K,I,5),LAMF(K,I,2), YAN20700 2 K=1,I) YAN20800 5405 FORMAT (2H *,1PE10.3,3H +-,E8.1,0PF9.3,1PE17.3,3H +-, YAN20900 1 E8.1,0PF9.3,2H *,1PE18.3) YAN21000 IF (.NOT.FFAIL(I,2)) WRITE (IWRITE,5405) (ASAVE(K,I,1), YAN21100 1 ASAVE(K,I,4),ASAVE(K,I,6),ASAVE(K,I,2),ASAVE(K,I,3),ASAVE(K,I,5),YAN21200 2 LAMF(K,I,2),K=1,I) YAN21300 IF (NOBASE) GO TO 410 YAN21400 K=I+1 YAN21500 WRITE (IWRITE,5404) ASAVE(K,I,1),ASAVE(K,I,4),ASAVE(K,I,6), YAN21600 1 ASAVE(K,I,2),ASAVE(K,I,3),ASAVE(K,I,5) YAN21700 410 WRITE (IWRITE,5400) YAN21800 IF (J.EQ.1 .AND. PNXBES.GT..95) WRITE (IWRITE,5410) JLAMNX, YAN21900 1 JLAMBS, PNXBES YAN22000 5410 FORMAT (93X,20(1H*)/34H APPROXIMATE PROBABILITY THAT THIS, YAN22100 1 56H SOLUTION IS REALLY BETTER THAN THE SECOND BEST SOLUTION, YAN22200 2 9H = * PNG(,I1,1H/,I1,3H) =,F6.3,2H */93X,20(1H*)) YAN22300 5412 FORMAT (/5H PNG(,I1,1H/,I1,1H),12X,1H=,F21.3,13X,4HNPHI,13X, YAN22400 1 1H=,1PE11.2,13X,25H(UNCORRECTED PNG WOULD BE,0PF8.3,1H)) YAN22500 IF (J .GT. 1) WRITE (IWRITE,5412) I,JLAMBS,PROBSV(I),ENPHSV(I), YAN22600 1 PROBLN(I) YAN22700 5414 FORMAT (22H ITERATIONS IN FIT =,I21,13X, YAN22800 1 18HSTD. DEV. OF FIT =,1PE11.4,13X, YAN22900 2 27HSIGNAL/NOISE RATIO OF FIT =,0PF7.0) YAN23000 DDUM=ABSYMX/SIGYSV(I) YAN23100 IF (WTED) WRITE (IWRITE,5414) ITERSV(I),SIGYSV(I) YAN23200 IF (.NOT.WTED) WRITE (IWRITE,5414) ITERSV(I),SIGYSV(I),DDUM YAN23300 5416 FORMAT (21H LAMBDA HELD BETEWEEN,1PE9.2,4H AND,E9.2,13X, YAN23400 1 4HNPHI,13X,1H=,E11.2) YAN23500 IF (J .EQ. 1) WRITE (IWRITE,5416) LAMNMX(1,1),LAMNMX(2,1), YAN23600 1 ENPHSV(JLAMBS) YAN23700 5418 FORMAT (/4H LAG,25X,5(4X,3HK =,I2),23X,16(1H*)/ YAN23800 1 29H PROB. RESIDUALS UNCORRELATED,5F9.3,4X, YAN23900 2 27HWEIGHTED AVERAGE = * PUNC =,F6.3,2H *) YAN24000 WRITE (IWRITE,5418) (K,K=1,5),PUNCOR YAN24100 5422 FORMAT (1X,131(1H=)) YAN24200 IF (L-1 .LE. LINEPG) WRITE (IWRITE,5422) YAN24300 IF (L .LE. LINEPG) WRITE (IWRITE,5422) YAN24400 370 CONTINUE YAN24500 500 IF (.NOT.REPEAT) RETURN YAN24600 REPEAT=.FALSE. YAN24700 CALL RESIDU (JLAMBS,.FALSE.,T,Y,SQRTW,YLYFIT,NMAX,PUNCOR,ASAVE) YAN24800 GO TO 340 YAN24900 END YAN25000 C--------------------- DOUBLE PRECISION VERSION 2B ---------------------LSQ00100 C SUBROUTINE LSTSQR. CONSTRAINED LEAST SQUARES FIT USING STEPWISE LSQ00200 C REGRESSION APPROACH. LSQ00300 C N = NO. OF DATA POINTS TO BE FIT, I.E., =N (RAWDAT=.TRUE.) OR LSQ00400 C =2*JLAM+IBASE (RAWDAT=.FALSE.). LSQ00500 C Y,T = RAW DATA, REGARDLESS OF RAWDAT. LSQ00600 C ITYPE=1,2,3,4 FOR PRELIMINARY DETERMINATION OF WEIGHTS, FIT TO RAW LSQ00700 C DATA WITH INTERPOLATION, FIT TO TRANSFORMS, AND FINAL FIT TO LSQ00800 C RAW DATA, RESPECTIVELY. LSQ00900 C-----------------------------------------------------------------------LSQ01000 C CALLS SUBPROGRAMS - VARDUM (A DUMMY NAME), PIVOT, ANLERR LSQ01100 C WHICH IN TURN CALL - PIVOT1, ETHEOR LSQ01200 C-----------------------------------------------------------------------LSQ01300 SUBROUTINE LSTSQR (LAMST,JLAM,RAWDAT,N,Y,T,NLIN,VARDUM, LSQ01400 1 GSE,NDIMG,SQRTW,YLYFIT,E,NMAX,TSTART,DELTAT,NT,NINT, LSQ01500 2 ITYPE,SPLIT,NDIME,INTER,PRLSQ) LSQ01600 EXTERNAL VARDUM LSQ01700 REAL LAMNMX, LAMST LSQ01800 LOGICAL RAWDAT, NOBASE, PRLSQ, FAILED, ALLPIV(2), PIVLAM, LSQ01900 1 LDUM, INTER, PIVALP, LBLOKA, LBLOKB, LBLOKC, SPLIT LSQ02000 DOUBLE PRECISION GSE, YLYFIT, E, TSTART, DELTAT, VARMIN, VAR, LSQ02100 1 PALPHA, PLAM, PLMTRY, DELTAP, DFCAPL, SFDE, ADUB, DDSE, LSQ02200 2 DTOLER, F, DF, VAROLD, VARG, R, DFLAT, DALPL, DUB, DDUB, VT, LSQ02300 3 DBLOKA, DBLOKB, DBLOKC, T LSQ02400 COMMON /BLOK/ DBLOKA(4828), VARMIN, DBLOKB(14), VAR, PALPHA(10), LSQ02500 1 PLAM(10), DBLOKC(100), PLMTRY(10), DELTAP(19), DFCAPL(10,9), LSQ02600 2 SFDE(9), ADUB(19,19), DDSE(9,9), DTOLER(19), F(19,10), LSQ02700 3 DF(19,10), DELTA, FHATMX, CONVRG, RBLOKA(171), PRECIS, LSQ02800 4 LAMNMX(2,2), RBLOKB(2), SIGMAP(19), PCERR(19), RBLOKC(252), LSQ02900 5 IBLOKA(199), MDIMA, IBLOKB(96), NF, IBLOKC, IWRITE, IBLOKD(5) LSQ03000 COMMON /BLOK/ ITER, JLAMP1, IBLOKE(9), LBLOKA, FAILED, LBLOKB, LSQ03100 1 NOBASE, LBLOKC(23), PIVALP(19) LSQ03200 DIMENSION PIVLAM(9), R(10), DALPL(10,9) LSQ03300 DIMENSION LAMST(JLAM), Y(NMAX), T(NMAX), LSQ03400 1 HOLRTH(2), SQRTW(NMAX), YLYFIT(NMAX), LSQ03500 2 E(NDIME,NLIN), TSTART(NINT), QT(3), VT(3), GSE(NDIMG,4), LSQ03600 3 DELTAT(NINT), NT(NINT), MXITER(4), RLAMMN(2) LSQ03700 DATA QMIN/.001/, HOLRTH/1H , 1H*/, MCONV/3/, NABORT/4/, LSQ03800 1 MXITER/30, 30, 40, 40/,RLAMMN/1.1, 2./,DFLAT/-1.D-2/ LSQ03900 ISPLIT=1 LSQ04000 IF (SPLIT) ISPLIT=2 LSQ04100 JLAMP1=JLAM+1 LSQ04200 NCONV=0 LSQ04300 NVARUP=0 LSQ04400 NFLAT=0 LSQ04500 ALLPIV(2)=.TRUE. LSQ04600 IF (NOBASE) PALPHA(JLAMP1)=0.D0 LSQ04700 PLAM(JLAMP1)=0.D0 LSQ04800 PLMTRY(JLAMP1)=0.D0 LSQ04900 PIVALP(JLAMP1)=.TRUE. LSQ05000 C-----------------------------------------------------------------------LSQ05100 C PUT PLAM IN ASCENDING ORDER AND SEPARATE THEM IF NECESSARY LSQ05200 C-----------------------------------------------------------------------LSQ05300 DO 150 J=1,JLAM LSQ05400 PIVLAM(J)=.FALSE. LSQ05500 150 DELTAP(J)=0.D0 LSQ05600 DO 160 J=1,JLAM LSQ05700 DUM=1.E+20 LSQ05800 DO 165 K=1,JLAM LSQ05900 IF (PIVLAM(K) .OR. LAMST(K).GT.DUM) GO TO 165 LSQ06000 L=K LSQ06100 DUM=LAMST(K) LSQ06200 165 CONTINUE LSQ06300 PIVLAM(L)=.TRUE. LSQ06400 IF (J .EQ. 1) PLAM(J)=DBLE(AMAX1(DUM,LAMNMX(1,ISPLIT))) LSQ06500 IF (J .GT. 1) PLAM(J)=DBLE(AMIN1(AMAX1(DUM,RLAMMN(ISPLIT)* LSQ06600 1 SNGL(PLAM(J-1))),LAMNMX(2,1))) LSQ06700 160 CONTINUE LSQ06800 VAROLD=1.D+20 LSQ06900 VARG=1.D+20 LSQ07000 ITER=-1 LSQ07100 Q=0. LSQ07200 IF (PRLSQ) WRITE (IWRITE,5205) LSQ07300 5205 FORMAT (/4H ITR,2X,8HVARIANCE,2X,9HDAMPING Q,3X,8HBASELINE, LSQ07400 1 5(5X,5HALPHA,3X,6HLAMBDA)) LSQ07500 C-----------------------------------------------------------------------LSQ07600 C START OF MAIN LOOP FOR THE LEAST SQUARES FIT LSQ07700 C-----------------------------------------------------------------------LSQ07800 210 NTRY=1 LSQ07900 215 ITER=ITER+1 LSQ08000 220 CALL VARDUM (Q,JLAM,VAR,T,SQRTW,NMAX,TSTART,DELTAT,NT, LSQ08100 1 NINT,.TRUE.,LDUM,E,Y,YLYFIT,NLIN,VARG,NTRY,NDIME,INTER, LSQ08200 2 GSE,NDIMG,IERROR) LSQ08300 IF (IERROR .EQ. 5) GO TO 800 LSQ08400 ALLPIV(1)=LDUM LSQ08500 IF (VAR.LE.VARG .OR. NTRY.EQ.2) GO TO 240 LSQ08600 NTRY=2 LSQ08700 Q=QG LSQ08800 GO TO 220 LSQ08900 240 DO 242 J=1,JLAM LSQ09000 R(J)=PALPHA(J) LSQ09100 IF (.NOT.RAWDAT) R(J)=PALPHA(J)*EXP(-SNGL(PLMTRY(J))*DELTA) LSQ09200 242 PLAM(J)=PLMTRY(J) LSQ09300 DUB=VAR-VAROLD LSQ09400 IF (.NOT.PRLSQ) GO TO 250 LSQ09500 NU=MIN0(5,JLAM) LSQ09600 DDUM=HOLRTH(1) LSQ09700 IF (DUB .GT. 0.D0) DDUM=HOLRTH(2) LSQ09800 DO 244 J=1,JLAM LSQ09900 PCERR(J)=HOLRTH(1) LSQ10000 IF (.NOT.PIVLAM(J)) PCERR(J)=HOLRTH(2) LSQ10100 244 CONTINUE LSQ10200 5240 FORMAT (1X,I2,1PD11.4,A1,E10.2,A1,D10.2,1X,5(1P2D9.2,A1)) LSQ10300 WRITE (IWRITE,5240) ITER,VAR,DDUM,Q,HOLRTH(NTRY),PALPHA(JLAMP1), LSQ10400 1 (R(J),PLAM(J),PCERR(J),J=1,NU) LSQ10500 IF (JLAM .LE. 5) GO TO 250 LSQ10600 5250 FORMAT (37X,4(1P2D9.2,A1)) LSQ10700 WRITE (IWRITE,5250) (R(J),PLAM(J),PCERR(J),J=6,JLAM) LSQ10800 5001 FORMAT (1X) LSQ10900 WRITE (IWRITE,5001) LSQ11000 C-----------------------------------------------------------------------LSQ11100 C TESTS FOR CONVERGENCE AND DIVERGENCE LSQ11200 C-----------------------------------------------------------------------LSQ11300 250 IF ((DABS(DUB).GT.DMAX1(VARMIN,VAROLD*CONVRG) .AND. RAWDAT) .OR. LSQ11400 1 (VAR.GT.DBLE(FLOAT(NF)*(CONVRG*FHATMX)**2) .AND. .NOT.RAWDAT)) LSQ11500 2 GO TO 270 LSQ11600 NCONV=NCONV+1 LSQ11700 IF (NCONV.GE.MCONV .OR. .NOT.RAWDAT) GO TO 800 LSQ11800 NVARUP=0 LSQ11900 GO TO 300 LSQ12000 270 NCONV=0 LSQ12100 IF (DUB .LE. 0.D0) GO TO 275 LSQ12200 NVARUP=NVARUP+1 LSQ12300 IF (NVARUP.LT.NABORT .OR. ITER.LT.JLAM) GO TO 280 LSQ12400 IERROR=1 LSQ12500 GO TO 800 LSQ12600 275 NVARUP=0 LSQ12700 280 IF (ITER.LT.MXITER(ITYPE)) GO TO 290 LSQ12800 IERROR=2 LSQ12900 GO TO 800 LSQ13000 290 IF (DUB/VAROLD.LT.DFLAT .OR. RAWDAT) GO TO 295 LSQ13100 NFLAT=NFLAT+1 LSQ13200 IF (NFLAT .LT. NABORT) GO TO 300 LSQ13300 IERROR=3 LSQ13400 GO TO 800 LSQ13500 295 NFLAT=0 LSQ13600 C-----------------------------------------------------------------------LSQ13700 C FINISH COMPUTATION OF DFCAPL (PARTIAL OF FCAP W.R.T. LAMBDA) LSQ13800 C-----------------------------------------------------------------------LSQ13900 300 DO 310 K=1,JLAM LSQ14000 DO 320 J=1,NLIN LSQ14100 DFCAPL(J,K)=DFCAPL(J,K)*PALPHA(K) LSQ14200 320 CONTINUE LSQ14300 DFCAPL(K,K)=DFCAPL(K,K)+SFDE(K) LSQ14400 310 CONTINUE LSQ14500 C-----------------------------------------------------------------------LSQ14600 C CALCULATE DALPL (PARTIAL OF PALPHA W.R.T. LAMBDA) FROM MATRIX LSQ14700 C PRODUCT ADUB*DFCAPL LSQ14800 C-----------------------------------------------------------------------LSQ14900 DO 400 J=1,NLIN LSQ15000 DO 400 K=1,JLAM LSQ15100 400 DALPL(J,K)=0.D0 LSQ15200 DO 410 K=1,JLAM LSQ15300 IF (.NOT.PIVALP(K)) GO TO 410 LSQ15400 DO 420 J=1,NLIN LSQ15500 IF (.NOT.PIVALP(J)) GO TO 420 LSQ15600 DUB=0.D0 LSQ15700 DO 430 I=1,NLIN LSQ15800 IF (PIVALP(I)) DUB=DUB+ADUB(J,I)*DFCAPL(I,K) LSQ15900 430 CONTINUE LSQ16000 DALPL(J,K)=DUB LSQ16100 420 CONTINUE LSQ16200 410 CONTINUE LSQ16300 C-----------------------------------------------------------------------LSQ16400 C CALCULATE ADUB (NORMAL EQUATION MATRIX FOR NONLINEAR LEAST SQUARES) LSQ16500 C AND INVERT IT LSQ16600 C-----------------------------------------------------------------------LSQ16700 IF (.NOT.RAWDAT) GO TO 550 LSQ16800 DO 510 K=1,JLAM LSQ16900 DO 520 J=1,K LSQ17000 DUB=0.D0 LSQ17100 DO 530 L=1,NLIN LSQ17200 530 DUB=DUB-DALPL(L,J)*DFCAPL(L,K) LSQ17300 IF (DABS(PALPHA(K)) .GT. DABS(PALPHA(J))) GO TO 540 LSQ17400 DUB=DUB+DALPL(J,K)*SFDE(J) LSQ17500 DUB=DUB+DALPL(K,J)*SFDE(K) LSQ17600 GO TO 545 LSQ17700 540 DUB=DUB+DALPL(K,J)*SFDE(K) LSQ17800 DUB=DUB+DALPL(J,K)*SFDE(J) LSQ17900 545 DDUB=PALPHA(K)*PALPHA(J)*DDSE(J,K) LSQ18000 ADUB(J,K)=DUB+DDUB LSQ18100 520 CONTINUE LSQ18200 DTOLER(K)=DDUB LSQ18300 ADUB(K,JLAMP1)=PALPHA(K)*SFDE(K) LSQ18400 510 CONTINUE LSQ18500 GO TO 590 LSQ18600 550 DO 560 J=1,JLAM LSQ18700 DTOLER(J)=0.D0 LSQ18800 R(J)=0.D0 LSQ18900 DO 565 K=J,JLAMP1 LSQ19000 565 ADUB(J,K)=0.D0 LSQ19100 560 CONTINUE LSQ19200 DO 570 NU=1,N LSQ19300 DO 575 K=1,JLAM LSQ19400 IF (.NOT.PIVALP(K)) GO TO 575 LSQ19500 DUB=0.D0 LSQ19600 DO 578 J=1,NLIN LSQ19700 578 DUB=DUB+F(NU,J)*DALPL(J,K) LSQ19800 DDUB=DF(NU,K)*PALPHA(K) LSQ19900 DTOLER(K)=DTOLER(K)+DDUB*DDUB LSQ20000 R(K)=DUB+DDUB LSQ20100 575 CONTINUE LSQ20200 DDUB=YLYFIT(NU) LSQ20300 DO 580 K=1,JLAM LSQ20400 IF (.NOT.PIVALP(K)) GO TO 580 LSQ20500 DUB=R(K) LSQ20600 DO 585 J=1,K LSQ20700 585 ADUB(J,K)=ADUB(J,K)+R(J)*DUB LSQ20800 ADUB(K,JLAMP1)=ADUB(K,JLAMP1)+R(K)*DDUB LSQ20900 580 CONTINUE LSQ21000 570 CONTINUE LSQ21100 590 DO 595 J=1,JLAM LSQ21200 DTOLER(J)=PRECIS*DMAX1(DTOLER(J),ADUB(J,J)) LSQ21300 595 R(J)=ADUB(J,JLAMP1) LSQ21400 ADUB(JLAMP1,JLAMP1)=0.D0 LSQ21500 CALL PIVOT (ADUB,MDIMA,JLAM,.FALSE.,.FALSE.,.TRUE., LSQ21600 1 DBLE(LAMNMX(1,1)),DBLE(LAMNMX(2,1)),PLAM,DTOLER,DELTAP,LDUM, LSQ21700 2 QG,PIVLAM,JLAM,QMAX,IERROR) LSQ21800 IF (IERROR .EQ. 5) GO TO 800 LSQ21900 ALLPIV(2)=LDUM LSQ22000 C-----------------------------------------------------------------------LSQ22100 C USE MODIFICATION-A OF G. E. P. BOX TO ESTIMATE BEST Q (FRACTIONAL LSQ22200 C STEP SIZE) LSQ22300 C-----------------------------------------------------------------------LSQ22400 VAROLD=VAR LSQ22500 DDUM=0. LSQ22600 DO 600 J=1,JLAM LSQ22700 600 DDUM=DDUM+SNGL(DELTAP(J)*R(J)) LSQ22800 IF (DDUM .GT. 0.) GO TO 610 LSQ22900 5610 FORMAT (42H GAUSS VECTOR POINTING IN WRONG DIRECTION.,2X, LSQ23000 1 44(2H**)) LSQ23100 IF (PRLSQ) WRITE (IWRITE,5610) LSQ23200 NTRY=2 LSQ23300 Q=AMIN1(QMIN,QMAX) LSQ23400 GO TO 215 LSQ23500 610 CALL VARDUM (QG,JLAM,VARG,T,SQRTW,NMAX,TSTART,DELTAT,NT, LSQ23600 1 NINT,.FALSE.,LDUM,E,Y,YLYFIT,NLIN,VAR,NTRY,NDIME,INTER, LSQ23700 2 GSE,NDIMG,IERROR) LSQ23800 IF (IERROR .EQ. 5) GO TO 800 LSQ23900 ALLPIV(1)=LDUM LSQ24000 DUM=SNGL(VARG-VAROLD)+2.*DDUM*QG LSQ24100 IF (DUM .LE. 0.) DUM=.25*DDUM*QG*QG LSQ24200 Q=AMIN1(DDUM*QG*QG/DUM,4.,QMAX) LSQ24300 IF (VARG.LE.VAROLD .OR. (Q.GE..125*QG .AND. RAWDAT)) GO TO 210 LSQ24400 C-----------------------------------------------------------------------LSQ24500 C PUT QUADRATIC PARABOLA THRU ACTUAL VARIANCE SURFACE WHEN A LSQ24600 C STRONG NONLINEARITY IS INDICATED BY A SMALL Q LSQ24700 C-----------------------------------------------------------------------LSQ24800 QT(1)=AMIN1(.125*QG,AMAX1(Q,.02)) LSQ24900 IF (Q .GT. .125*QG) QT(1)=Q LSQ25000 CALL VARDUM (QT(1),JLAM,DUB,T,SQRTW,NMAX,TSTART,DELTAT,NT, LSQ25100 1 NINT,.FALSE.,LDUM,E,Y,YLYFIT,NLIN,VAR,NTRY,NDIME,INTER, LSQ25200 2 GSE,NDIMG,IERROR) LSQ25300 IF (IERROR .EQ. 5) GO TO 800 LSQ25400 VT(1)=DUB LSQ25500 IF (VT(1) .GE. VAROLD) GO TO 615 LSQ25600 QT(2)=4.*QT(1) LSQ25700 IF (QT(1) .LE. .125*QG) GO TO 617 LSQ25800 QT(2)=0. LSQ25900 VT(2)=VAROLD LSQ26000 QT(3)=QG LSQ26100 VT(3)=VARG LSQ26200 GO TO 630 LSQ26300 615 QT(2)=.5*QT(1) LSQ26400 QT(2)=AMIN1(QT(2),.05) LSQ26500 617 CALL VARDUM (QT(2),JLAM,DUB,T,SQRTW,NMAX,TSTART,DELTAT,NT, LSQ26600 1 NINT,.FALSE.,LDUM,E,Y,YLYFIT,NLIN,VAR,NTRY,NDIME,INTER, LSQ26700 2 GSE,NDIMG,IERROR) LSQ26800 IF (IERROR .EQ. 5) GO TO 800 LSQ26900 VT(2)=DUB LSQ27000 IF (VT(1).GE.VAROLD .OR. VT(2).GE.VT(1)) GO TO 620 LSQ27100 QT(3)=QG LSQ27200 VT(3)=VARG LSQ27300 GO TO 630 LSQ27400 620 QT(3)=0. LSQ27500 VT(3)=VAROLD LSQ27600 630 VARG=DMIN1(VT(1),VT(2),VARG) LSQ27700 IF (VT(1) .LE. VARG) QG=QT(1) LSQ27800 IF (VT(2) .LE. VARG) QG=QT(2) LSQ27900 DUM=SNGL((VT(1)-VT(3)))/(QT(1)-QT(3)) LSQ28000 DDUM=(DUM-SNGL((VT(2)-VT(3)))/(QT(2)-QT(3)))/(QT(1)-QT(2)) LSQ28100 IF (DDUM .GT. 0.) GO TO 650 LSQ28200 640 Q=AMIN1(Q,.125*QT(2),.01) LSQ28300 GO TO 210 LSQ28400 650 DUM=.5*((QT(1)+QT(3))-DUM/DDUM) LSQ28500 IF (DUM .LE. 0.) GO TO 640 LSQ28600 Q=AMIN1(DUM,8.,QMAX) LSQ28700 GO TO 210 LSQ28800 C-----------------------------------------------------------------------LSQ28900 C CHECK FOR ABNORMAL EXIT AND CALL ANLERR LSQ29000 C-----------------------------------------------------------------------LSQ29100 800 FAILED=IERROR.NE.0 .OR. .NOT.ALLPIV(1) LSQ29200 IF (.NOT.RAWDAT) GO TO 890 LSQ29300 FAILED=FAILED .OR. .NOT.ALLPIV(2) LSQ29400 IF (.NOT.FAILED .OR. ITYPE.EQ.4) CALL ANLERR (JLAM,T,SQRTW,NMAX,E,LSQ29500 1 NLIN,NDIME,TSTART,DELTAT,NT,NINT,N,ITYPE,INTER,PRLSQ) LSQ29600 890 IF (.NOT.(PRLSQ .AND. FAILED)) GO TO 950 LSQ29700 5889 FORMAT (/47H VARIANCE DID NOT SIGNIFICANTLY DECREASE IN THE, LSQ29800 1 5H LAST,I2,7H TIMES.,5X,6(11H***TOO SLOW)) LSQ29900 IF (IERROR .EQ. 3) WRITE (IWRITE,5889) NABORT LSQ30000 5890 FORMAT (/41H MAXIMUM ITERATIONS WITHOUT CONVERGENCE. , LSQ30100 1 7(13H***ITERATIONS)) LSQ30200 IF (IERROR .EQ. 2) WRITE (IWRITE,5890) LSQ30300 5892 FORMAT (/28H VARIANCE INCREASED THE LAST,I2,8H TIMES. , LSQ30400 1 3X,7(13H***DIVERGENCE)) LSQ30500 IF (IERROR .EQ. 1) WRITE (IWRITE,5892) NVARUP LSQ30600 5894 FORMAT (/42H NEARLY SINGULAR NORMAL EQUATIONS MATRIX. , LSQ30700 1 6(15H****SINGULARITY)) LSQ30800 IF (IERROR .EQ. 5) WRITE (IWRITE,5894) LSQ30900 C-----------------------------------------------------------------------LSQ31000 C PUT PLAM IN ASCENDING ORDER LSQ31100 C-----------------------------------------------------------------------LSQ31200 950 DO 960 J=1,JLAM LSQ31300 DUB=1.E+20 LSQ31400 DO 985 K=J,JLAM LSQ31500 IF (PLAM(K) .GT. DUB) GO TO 985 LSQ31600 L=K LSQ31700 DUB=PLAM(K) LSQ31800 985 CONTINUE LSQ31900 PLAM(L)=PLAM(J) LSQ32000 PLAM(J)=DUB LSQ32100 DUB=PALPHA(L) LSQ32200 PALPHA(L)=PALPHA(J) LSQ32300 PALPHA(J)=DUB LSQ32400 IF (.NOT.RAWDAT) GO TO 960 LSQ32500 DUM=SIGMAP(L) LSQ32600 SIGMAP(L)=SIGMAP(J) LSQ32700 SIGMAP(J)=DUM LSQ32800 K=L+JLAM LSQ32900 JK=J+JLAM LSQ33000 DUM=SIGMAP(K) LSQ33100 SIGMAP(K)=SIGMAP(JK) LSQ33200 SIGMAP(JK)=DUM LSQ33300 PCERR(J)=100.*SIGMAP(J)/SNGL(PLAM(J)) LSQ33400 IF (DABS(PALPHA(J)) .GT. 0.D0) PCERR(JK)=100.*SIGMAP(JK)/ABS( LSQ33500 1 SNGL(PALPHA(J))) LSQ33600 960 CONTINUE LSQ33700 RETURN LSQ33800 END LSQ33900 C--------------------- DOUBLE PRECISION VERSION 2B ---------------------EVR00100 C SUBROUTINE EVAR. FOR FRACTIONAL STEP SIZE Q, EVALUATES PALPHA EVR00200 C (LINEAR LEAST SQUARES PARAMETERS), VAR (WEIGHTED VARIANCE) AND EVR00300 C ACCUMULATES IT IN VGRID2 (COARSE GRID SUM), AND, IF EVR00400 C INVERT=.TRUE., EVALUATES ARRAYS SFDE, DFCAPL, AND FULL UPPER EVR00500 C TRIANGLES OF DDSE AND SE FOR LATER USE IN NONLINEAR LEAST EVR00600 C SQUARES. EVR00700 C-----------------------------------------------------------------------EVR00800 C CALLS SUBPROGRAMS - ETHEOR, PIVOT EVR00900 C WHICH IN TURN CALL - PIVOT1 EVR01000 C-----------------------------------------------------------------------EVR01100 SUBROUTINE EVAR (Q,JLAM,VAR,T,SQRTW,NMAX,TSTART,DELTAT,NT,NINT, EVR01200 1 INVERT,ALLPIV,E,Y,YLYDUM,NLIN,VARG,NTRY,NDIME,INTER,GSE,NDIMG, EVR01300 2 IERROR) EVR01400 LOGICAL INVERT, ALLPIV, LDDUM, INTER, LBLOKA, LBLOKB, LBLOKC, EVR01500 1 NOTINV, NOBASE, WTED, REGINT, NONNEG, PIVALP EVR01600 REAL LAMNMX EVR01700 DOUBLE PRECISION TSTART, DELTAT, E, YLYDUM, GSE, SUMWTY, SUMWT, EVR01800 1 ZERO, DGRIDE, GESTL, PALPHA, PLAM, SE, PLMTRY, DELTAP, DFCAPL, EVR01900 2 SFDE, ADUB, DDSE, DTOLER, DBLOKA, DBLOKB, DBLOKC, ADUM, DUM, EVR02000 3 DDUM, DDDUM, XLAM, HZ, H, HH, DX, DDX, S, FACTOR, XEXP, EVR02100 4 YLYFIT, VAR, VARG, T EVR02200 COMMON /BLOK/ DBLOKA(4829), SUMWTY, SUMWT, ZERO(10), DGRIDE, EVR02300 1 GESTL, DBLOKB, PALPHA(10), PLAM(10), SE(10,10), PLMTRY(10), EVR02400 2 DELTAP(19), DFCAPL(10,9), SFDE(9), ADUB(19,19), DDSE(9,9), EVR02500 3 DTOLER(19), DBLOKC(380), RBLOKA(174), PRECIS, LAMNMX(2,2), EVR02600 4 RBLOKB(40), VGRID2(252), JGEADD, NBINOM(20,9), IBLOKA(9), EVR02700 5 NPERG2(9), MDIMA, IBLOKB(4), NGRID2(9), IBLOKC(84), N EVR02800 COMMON /BLOK/ IBLOKD(17), REGINT, LBLOKA(2), NOBASE, WTED, EVR02900 1 LBLOKB(3), NONNEG, LBLOKC(18), PIVALP(19) EVR03000 DIMENSION ADUM(10), FACTOR(9), XEXP(9), JG2(9) EVR03100 DIMENSION T(NMAX), SQRTW(NMAX), TSTART(NINT), DELTAT(NINT), EVR03200 1 NT(NINT), E(NDIME,NLIN), Y(NMAX), YLYDUM(1), GSE(NDIMG,4) EVR03300 DO 110 J=1,JLAM EVR03400 110 PLMTRY(J)=DMAX1(DBLE(LAMNMX(1,1)), EVR03500 1 DMIN1(DBLE(LAMNMX(2,1)),PLAM(J)+Q*DELTAP(J))) EVR03501 NLINP=NLIN+1 EVR03600 DO 160 K=1,NLINP EVR03700 160 ADUB(K,NLINP)=0.D0 EVR03800 IF (NOBASE) GO TO 170 EVR03900 PLMTRY(NLIN)=0.D0 EVR04000 ADUB(NLIN,NLIN)=SUMWT EVR04100 ADUB(NLIN,NLINP)=SUMWTY EVR04200 170 NOTINV=.NOT.INVERT EVR04300 IF (REGINT) GO TO 300 EVR04400 C-----------------------------------------------------------------------EVR04500 C ORDINARY COMPUTATION OF E (EXPONENTIALS), ADUB (NORMAL EQUATIONS EVR04600 C MATRIX FOR LINEAR LEAST SQUARES) AND, IF INVERT=.TRUE., PART EVR04700 C OF DFCAPL (PARTIAL OF FCAP W.R.T. LAMBDA) WHEN REGINT=.FALSE. EVR04800 C-----------------------------------------------------------------------EVR04900 DO 190 K=1,NLIN EVR05000 DO 195 J=1,K EVR05100 IF (J .GT. JLAM) GO TO 190 EVR05200 ADUB(J,K)=0.D0 EVR05300 DFCAPL(K,J)=0.D0 EVR05400 195 IF (K .LE. JLAM) DDSE(J,K)=0.D0 EVR05500 190 CONTINUE EVR05600 DO 200 J=1,N EVR05700 DUM=-T(J) EVR05800 DO 210 K=1,JLAM EVR05900 ADUM(K)=DBLE(SQRTW(J))*DEXP(PLMTRY(K)*DUM) EVR06000 210 E(J,K)=ADUM(K) EVR06100 IF (NOBASE) GO TO 215 EVR06200 ADUM(NLIN)=DBLE(SQRTW(J)) EVR06300 E(J,NLIN)=ADUM(NLIN) EVR06400 215 DO 220 K=1,NLIN EVR06500 DDUM=ADUM(K) EVR06600 LDDUM=K .LE. JLAM EVR06700 DO 230 L=1,K EVR06800 IF (L .GT. JLAM) GO TO 220 EVR06900 DDDUM=ADUM(L)*DDUM EVR07000 ADUB(L,K)=ADUB(L,K)+DDDUM EVR07100 IF (NOTINV) GO TO 230 EVR07200 DDDUM=DDDUM*DUM EVR07300 DFCAPL(K,L)=DFCAPL(K,L)-DDDUM EVR07400 IF (LDDUM) DDSE(L,K)=DDSE(L,K)+DDDUM*DUM EVR07500 230 CONTINUE EVR07600 ADUB(K,NLINP)=ADUB(K,NLINP)+DDUM*DBLE(Y(J)) EVR07700 220 CONTINUE EVR07800 200 CONTINUE EVR07900 GO TO 500 EVR08000 C-----------------------------------------------------------------------EVR08100 C EVALUATE ADUB AND, IF INVERT=.TRUE., PART OF DFCAPL AND DDSE EVR08200 C (SECOND DERIVATIVES OF EXPONENTIAL SUMS W.R.T. LAMBDA) EVR08300 C BY HERMITE AND LAGRANGE INTERPOLATION, IF REGINT=INTER=.TRUE. EVR08400 C-----------------------------------------------------------------------EVR08500 300 IF (.NOT.INTER) GO TO 400 EVR08600 DO 310 K=1,NLIN EVR08700 LDDUM=K .LE. JLAM EVR08800 DO 320 J=1,K EVR08900 IF (J .GT. JLAM) GO TO 310 EVR09000 XLAM=PLMTRY(K)+PLMTRY(J) EVR09100 DUM=(DLOG(XLAM)-GESTL)/DGRIDE EVR09200 JZ=IDINT(DUM+.5D0) EVR09300 HZ=DUM-DBLE(FLOAT(JZ)) EVR09400 IF (DABS(HZ) .GT. 1.D-5) GO TO 330 EVR09500 ADUB(J,K)=GSE(JZ,1) EVR09600 IF (NOTINV) GO TO 320 EVR09700 DFCAPL(K,J)=-GSE(JZ,2)/XLAM EVR09800 IF (LDDUM) DDSE(J,K)=(GSE(JZ,3)-GSE(JZ,2))/(XLAM*XLAM) EVR09900 GO TO 320 EVR10000 330 JM=IDINT(DUM) EVR10100 JP=JM+1 EVR10200 H=DUM-DBLE(FLOAT(JM)) EVR10300 HH=1.D0-H EVR10400 DDUM=1.D0+3.D0*H EVR10500 DDDUM=1.D0+3.D0*HH EVR10600 DX=H*DGRIDE EVR10700 DDX=-HH*DGRIDE EVR10800 ADUB(J,K)=HH*HH*HH*((DDUM+6.D0*H*H)*GSE(JM,1)+DX*(DDUM*GSE(JM,2)+ EVR10900 1 .5D0*DX*GSE(JM,3)))+H*H*H*((DDDUM+6.D0*HH*HH)*GSE(JP,1)+ EVR11000 2 DDX*(DDDUM*GSE(JP,2)+.5D0*DDX*GSE(JP,3))) EVR11100 IF (NOTINV) GO TO 320 EVR11200 HH=HZ*HZ EVR11300 DX=(HH-1.D0)*HZ*.5D0 EVR11400 DUM=1.D0/(HZ+1.D0) EVR11500 DDUM=1.D0/HZ EVR11600 DDDUM=1.D0/(HZ-1.D0) EVR11700 H=DX*DX*(((3.D0*HZ+4.D0)*GSE(JZ-1,2)*DUM+DGRIDE*GSE(JZ-1,3))*DUM+ EVR11800 1 4.D0*DDUM*(DDUM*GSE(JZ,2)+DGRIDE*GSE(JZ,3))+DDDUM*((4.D0-3.D0* EVR11900 2 HZ)*GSE(JZ+1,2)*DDDUM+DGRIDE*GSE(JZ+1,3))) EVR12000 DFCAPL(K,J)=-H/XLAM EVR12100 IF (.NOT.LDDUM) GO TO 320 EVR12200 DUM=DX*(HH-4.D0)*(HH-9.D0)*(GSE(JZ+3,3)/(HZ-3.D0)+GSE(JZ-3,3)/ EVR12300 1 (HZ+3.D0)-6.D0*(GSE(JZ+2,3)/(HZ-2.D0)+GSE(JZ-2,3)/(HZ+2.D0))+ EVR12400 2 15.D0*(GSE(JZ+1,3)*DDDUM+GSE(JZ-1,3)*DUM)-20.D0*GSE(JZ,3)*DDUM)/ EVR12500 3 360.D0 EVR12600 DDSE(J,K)=(DUM-H)/(XLAM*XLAM) EVR12700 320 CONTINUE EVR12800 DUM=(DLOG(PLMTRY(K))-GESTL)/DGRIDE EVR12900 JZ=IDINT(DUM+.5D0) EVR13000 HZ=DUM-DBLE(FLOAT(JZ)) EVR13100 IF (DABS(HZ) .GT. 1.D-5) GO TO 340 EVR13200 ADUB(K,NLINP)=GSE(JZ,4) EVR13300 GO TO 310 EVR13400 340 HH=HZ*HZ EVR13500 ADUB(K,NLINP)=HZ*(HH-1.D0)*(HH-4.D0)*(HH-9.D0)*(GSE(JZ+3,4)/(HZ- EVR13600 1 3.D0)+GSE(JZ-3,4)/(HZ+3.D0)-6.D0*(GSE(JZ+2,4)/(HZ-2.D0)+ EVR13700 2 GSE(JZ-2,4)/(HZ+2.D0))+15.D0*(GSE(JZ+1,4)/(HZ-1.D0)+GSE(JZ-1,4)/ EVR13800 3 (HZ+1.D0))-20.D0*GSE(JZ,4)/HZ)/720.D0 EVR13900 310 CONTINUE EVR14000 GO TO 500 EVR14100 C-----------------------------------------------------------------------EVR14200 C EVALUATE ADUB, PART OF DFCAPL, AND DDSE BY CALLING ETHEOR EVR14300 C IF REGINT=.TRUE. AND INTER=.FALSE. EVR14400 C-----------------------------------------------------------------------EVR14500 400 DO 410 K=1,NLIN EVR14600 LDDUM=K .LE. JLAM EVR14700 DO 420 J=1,K EVR14800 IF (J .GT. JLAM) GO TO 410 EVR14900 XLAM=PLMTRY(K)+PLMTRY(J) EVR15000 CALL ETHEOR (XLAM,INVERT,WTED,T,SQRTW,Y,NMAX,TSTART,DELTAT, EVR15100 1 NT,NINT,DUM,DDUM,DDDUM,S,.TRUE.,.FALSE.) EVR15200 ADUB(J,K)=DUM EVR15300 IF (NOTINV) GO TO 420 EVR15400 DFCAPL(K,J)=-DDUM/XLAM EVR15500 IF (LDDUM) DDSE(J,K)=(DDDUM-DDUM)/(XLAM*XLAM) EVR15600 420 CONTINUE EVR15700 CALL ETHEOR (PLMTRY(K),.FALSE.,WTED,T,SQRTW,Y,NMAX,TSTART,DELTAT, EVR15800 1 NT,NINT,DUM,DDUM,DDDUM,S,.FALSE.,.TRUE.) EVR15900 ADUB(K,NLINP)=S EVR16000 410 CONTINUE EVR16100 500 DO 510 K=1,NLIN EVR16200 DTOLER(K)=10.*PRECIS*ADUB(K,K) EVR16300 DO 515 J=1,K EVR16400 515 SE(J,K)=ADUB(J,K) EVR16500 510 CONTINUE EVR16600 C-----------------------------------------------------------------------EVR16700 C SOLVE NORMAL EQUATIONS FOR PALPHA USING STEPWISE REGRESSION EVR16800 C-----------------------------------------------------------------------EVR16900 CALL PIVOT (ADUB,MDIMA,NLIN,.FALSE.,INVERT,NONNEG,0.D0, EVR17