CFTN77,Q,T,S C$CDS ON C$FILES(0,3) C$EMA //,/FITCOM/ C$MSEG 6 C HP1000 directives precede this comment C INCLUDE 'CALPGM.INC' C PROGRAM CALFIT C********************************************************************** C C Copyright (C) 1989, California Institute of Technology C All rights reserved. U. S. Government Sponsorship under C NASA Contract NAS7-918 is acknowledged. C C Herbert M. Pickett, 20 Mar 1989 C HMP, change code for computing constraints C HMP, change PRCOR C HMP, change for generating binary file of variance C HMP, constained derivatives calculated within HAMX C HMP, IDPAR is now REAL*8 C C C THIS IS A GENERALIZED LINE FITTING PROGRAM C C IT FITS LINES TO PARAMETERS IN A MODEL HAMILTONIAN C C BLENDED LINES ARE TREATED SPECIALLY: C IF THE EXPTL.FREQ. IS THE SAME TO 1 HZ THEN ONLY THE INVERSE ERROR C AVERAGED FREQUENCIES ARE USED IN THE FIT FOR THE BLEND C C*********************************************************************** INTEGER NDPAR,NDEGY,NDFIT,LUPAR,LULIN,LUFIT,LUBAK,LUVAR, : LUBIN,LULBL,LIMPAR,LIMLIN,NDFITM INTEGER*4 HEAPLN,LPTR,EGY,TEIG,EGYDER,VAR,FITBGN,LV,LF,LB C standard dimensions: PARAMETER (NDPAR=999,NDEGY=450,NDFIT=NDPAR+1,HEAPLN=300000) C DOS-PC(16 bit compiler) dimensions: C PARAMETER (NDPAR=180,NDEGY=181,NDFIT=NDPAR+1,HEAPLN=20000) PARAMETER (NDFITM=NDFIT-1) INTEGER NFILE,LPAR,LLIN,LFIT,LBAK,LVAR,LBIN PARAMETER (NFILE=6,LPAR=1,LLIN=2,LFIT=3,LBAK=4,LVAR=5,LBIN=6) INTEGER QNUM(12) CHARACTER*10 PARLBL(NDPAR) CHARACTER*51 PARE CHARACTER*36 AQNUM CHARACTER*80 NAMFIL,CARD CHARACTER*64 CFIL(NFILE) CHARACTER*4 EXT(NFILE) INTEGER RQEXIT,GETPAR,INITL,I0,I1,ITMP,NSIZE,NLINE,NPAR, : NXPAR,IQNFMT,ITD,NQN,NPARX,NDEL,I,II,NMXEGY,ITR, : NBLKPF,K,LU,NITR,LSTF,LBLK,LNEXT,LINE,IBLK, : INDX,IFAC,NCOD,N,NF,NRJ,NFIR,ICNT,LBLND,J,IFLG, : IBLND,HEAD,PCARD,SETOPT,SETBLK,NOPTN,NPARTOT INTEGER*2 CLEN5 PARAMETER (I0=0,I1=1,CLEN5=5) REAL*8 IDPAR(0:NDPAR),PAR(NDPAR),DPAR(NDFIT),ERP(NDPAR), : ERPAR(NDPAR),FIT(NDFIT,NDFITM),FQFAC(4), : ONE,ZERO(2),XERR,XWT,BIGF,BIGD,THRSH,XSQFAC,EPS, : XERRMX,PARFAC,XSQ,XSQL,XSQIR,XSQMW,DDOT,SUM,XSQOLD, : XFRQ,EX,CERR,CALERR,AFRQ,ADIF,RERR,FRQ,DIF,STFAC, : SCALE,DENF,FQFACQ,AVGMW,AVGIR,OLDIF(NDPAR), : PARD(NDPAR),PARBAS PARAMETER (ONE=1.,EPS=1.5E-38) EXTERNAL BRKQR REAL*8 HEAP(0:HEAPLN) COMMON HEAP COMMON /FITCOM/ FIT,PAR,DPAR,ERP,ERPAR,OLDIF,IDPAR SAVE DATA ZERO/0.,0./ C C >>>>>> BEGIN PROGRAM C C following statement is for PRIME computer to enable ^P hot key C CALL MKON$P('QUIT$',CLEN5,BRKQR) BIGF=9999999.9D0 BIGD=999.99999D0 STFAC=1. FQFAC(1)=ONE FQFAC(2)=-1 LULBL=NFILE+1 LUPAR=LPAR LULIN=LLIN LUFIT=LFIT LUBAK=LBAK LUVAR=LVAR LUBIN=LBIN EXT(LPAR)='par' EXT(LLIN)='lin' EXT(LFIT)='fit' EXT(LBAK)='bak' EXT(LVAR)='var' EXT(LBIN)='unf' C C open read and write files C CALL FILGET(CFIL,EXT,NFILE) CALL FILBAK(LUPAR,CFIL(LPAR),LUBAK,CFIL(LBAK)) CALL OPENLU(LUBAK,CFIL(LBAK),'R') REWIND(LUBAK) CALL OPENLU(LUFIT,CFIL(LFIT),'WP') C C read in run parameters C READ (LUBAK,520,END=999) CARD CALL CHTIME(CARD(51:80)) WRITE (LUFIT,520) CARD WRITE (*,521) CARD PAR(1)=NDPAR PAR(2)=32767 PAR(3)=1 PAR(4)=0 PAR(5)=1.E-10 PAR(6)=1.E+06 PAR(7)=1. PAR(8)=1. READ(LUBAK,520,END=999) CARD I=8 IF(PCARD(CARD,PAR,I).EQ.0) GO TO 999 NPARTOT=PAR(1) NLINE= PAR(2) NITR= PAR(3) NXPAR= PAR(4) THRSH= PAR(5) XERRMX= PAR(6) PARFAC= PAR(7) FQFACQ= PAR(8) FQFAC(3)=FQFACQ/29979.2458D0 FQFAC(4)=-FQFAC(3) PARFAC=ABS(PARFAC) IF(ABS(PARFAC-ONE).GT.1.E-10) : WRITE(*,*) 'PARAMETER ERRORS SCALED BY ',PARFAC C .. read in option card C IQNFMT=0 NOPTN=SETOPT(LUBAK,IQNFMT,ITD,NAMFIL) IF(NOPTN.LT.0) GO TO 999 ITMP=10 NQN=MOD(IQNFMT,ITMP)*2 C C read in parameters and variance ( and return true if variance ) C IF (NPARTOT .GT.NDPAR ) THEN WRITE(*,*) ' Warning: number of parameters less' : ,' than requested',NPARTOT,NDPAR NPARTOT =NDPAR ENDIF LIMPAR=NPARTOT LIMLIN=NLINE WRITE (LUFIT,540) NLINE,NPARTOT,NITR,THRSH,XERRMX WRITE (*,540) NLINE,NPARTOT,NITR,THRSH,XERRMX VAR=HEAPLN NPAR=NPARTOT INITL=GETPAR(LUBAK,VAR,HEAP,ERP,NPAR,IDPAR,PAR,I0,LUFIT) NPARTOT=IDPAR(0) NDEL=NPAR+1 IF (NPAR .GT.NDFITM ) THEN WRITE(*,*) ' Warning: number of fitted parameters bigger' : ,' than dimension',NPAR,NDFITM STOP 'big heap' ENDIF C .. fixup supplied values IF (NXPAR.GT.NPARTOT.OR.NXPAR.LT.0) NXPAR=0 IF (NXPAR.GT.0) WRITE(LUFIT,545) NXPAR K=NPARTOT-NXPAR II=0 DO 13 I=1,NPARTOT IF(IDPAR(I).GE.0) THEN II=II+1 PARD(II)=PAR(I) IF(I.LE.K) NPARX=II ENDIF 13 CONTINUE IF (THRSH.LT.0.) THEN XSQFAC=100. IF (THRSH.GT.-EPS) THRSH=-EPS ELSE XSQFAC=1.001 IF (THRSH.LT.EPS) THRSH=EPS ENDIF IF (XERRMX.LT.EPS) XERRMX=1.E+06 IF (ABS(FQFACQ-ONE).GE.1.D-10) WRITE(LUFIT,*) : ' IR Frequencies Scaled by ',FQFACQ LF=VAR-NPAR FITBGN=LF-NPAR IF(FITBGN.LT.0) STOP 'big heap' IF (INITL.NE.0) THEN C ..initialize fit from supplied variance LV=VAR II=0 DO 15 I=1,NPAR IF(ABS( HEAP(LV+II) ).LT.EPS) STOP 'Bad .par' CALL DCOPY(I,HEAP(LV),I1,HEAP(LF+II),NPAR) LV=LV+NPAR 15 II=I C .. invert lower triangle CALL DTRI(HEAP(LF),NDPAR,NPAR,I0) ELSE C ..initialize fit when there is no supplied variance CALL DCOPY(NPAR,ZERO,I0,HEAP(LF),I1) DO 20 I=1,NPAR HEAP(LF)=ONE/HEAP(LF+NPAR) 20 LF=LF+NDEL CALL DCOPY(NPAR,ZERO,I0,HEAP(VAR),NDEL) ENDIF CALL DCOPY(NPAR,ZERO,I0,HEAP(FITBGN),I1) LPTR=FITBGN CALL HEAPNU(NPAR,LPTR) IF(LPTR.GT.FITBGN) THEN LPTR=LPTR-FITBGN WRITE(*,*) ' Warning: HEALPLN should be larger by ',LPTR STOP 'big HEAP' ENDIF C C read lines C CALL OPENLU(LULIN,CFIL(LLIN),'R') CALL LINEIN(LULIN,NLINE,IQNFMT,NBLKPF) CALL OPENLU(LULIN,CFIL(LLIN),'C') C .. initialize block structure NMXEGY=NDEGY K=SETBLK(LUFIT,IDPAR,PAR,NBLKPF,NMXEGY) C .. set parameter labels CALL GETLBL(IDPAR,PARLBL,LULBL,NAMFIL,K) WRITE (LUFIT,550) C .. establish limits on number of lines LPTR=NDEL+NMXEGY LPTR=FITBGN-NMXEGY*LPTR EGY=LPTR CALL HEAPNU(NPAR,LPTR) IF(LPTR.GT.EGY) THEN C .. downsize NMXEGY NMXEGY=0 EGY=LPTR-EGY WRITE(*,*) ' Warning: HEALPLN should be larger by ',EGY LPTR=FITBGN-LPTR EX=0.5*NDEL SUM=EX*EX+LPTR IF(SUM.GT.0) NMXEGY = SQRT(SUM) - EX IF(NMXEGY.LE.0) STOP 'big HEAP' WRITE(*,*) ' Energy Dimension reduced to ',NMXEGY ENDIF K=-1 I=RQEXIT(K) LPTR=NMXEGY TEIG=FITBGN-LPTR*NMXEGY EGY =TEIG -LPTR*NDEL C convert lines LU=LUFIT IF(NITR.GE.0) LU=-LU CALL LINEIX(LU,NLINE,NBLKPF,IQNFMT,HEAD) C******************************************************************** C START ITERATION C ITR=0 NITR=ABS(NITR) XSQOLD=1.E+30 C C find energies and derivatives C 150 LSTF=0 LBLK=0 LNEXT=HEAD C operator interrupt ? 200 IF ( RQEXIT(0).NE.0 ) GO TO 465 CALL GETDBK(LNEXT,LINE,IBLK,INDX,INITL,IFAC) IF(IBLK.NE.LBLK) THEN C .. get size of block CALL GETQN(IBLK,I0,QNUM,NSIZE,NCOD) IF (NSIZE.EQ.0) GO TO 200 LBLK=IBLK IF (NSIZE.GT.NMXEGY) THEN WRITE (LUFIT,620) IBLK,NSIZE STOP ENDIF K=(IBLK-1)/NBLKPF IF (LSTF.NE.K) THEN WRITE(*,*) 'Starting Quantum ',K LSTF=K ENDIF C .. get energies and derivatives EGYDER=EGY+NSIZE CALL HAMX (IBLK,NSIZE,IDPAR,PAR,HEAP(EGY), : HEAP(TEIG),HEAP(EGYDER),I0) ENDIF C .. save energies and derivatives for lines in this block CALL DNUADD(NPAR,NPARX,INITL,INDX,IFAC,HEAP(EGY), : HEAP(EGYDER),NSIZE,LINE,PARD,FQFAC) C repeat until no more lines IF(LNEXT.NE.0) GO TO 200 C C initialize least squares matrix C XSQL=0. XSQ=0. XSQIR=0. XSQMW=0. AVGMW=0. AVGIR=0. LB=FITBGN LF=FITBGN+NPAR DO 240 K=1,NPAR N=NDEL-K CALL DCOPY(N,HEAP(LF),I1,FIT(K,K),I1) FIT(NDEL,K)=DDOT(N,HEAP(LB),I1,FIT(K,K),I1) XSQL=XSQL+FIT(NDEL,K)**2 LB=LB+1 240 LF=LF+NDEL NF=0 NRJ=0 NFIR=0 ICNT=-1 WRITE (LUFIT,630) LINE=1 C C form least squares matrix, LOOP over lines C 330 IF (ICNT.LE.0) THEN WRITE(*,*) 'Fitting Line ',LINE ICNT=50+ICNT ENDIF LBLND=LINE SUM=0. IFLG=1 C ..DO UNTIL all elements of blend are found 260 ICNT=ICNT-1 CALL FRQDAT(LBLND,IBLND,XFRQ,XWT,XERR,QNUM) EX=ABS(XWT/XERR) SUM=SUM+EX C .. accumulate line contributions CALL DNUGET(IFLG,NDEL,EX,LBLND,DPAR) IF(IBLND.EQ.0) THEN LBLND=LBLND+1 IFLG=-1 GO TO 260 ENDIF C ..END DO C C .. calculate errors CERR=CALERR(NPAR,HEAP(VAR),DPAR)/SUM AFRQ=DPAR(NDEL)/SUM ADIF=XFRQ-AFRQ RERR=ADIF*SUM DPAR(NDEL)=RERR IF (ABS(RERR) .LT. XERRMX) THEN NF=NF+1 XSQL=XSQL+RERR*RERR IF (XERR.LT.0.) THEN AVGIR=AVGIR+ADIF XSQIR=XSQIR+ADIF*ADIF NFIR=NFIR+1 ELSE AVGMW=AVGMW+ADIF XSQMW=XSQMW+ADIF*ADIF ENDIF C .. rotate line into FIT matrix CALL JELIM(FIT,DPAR,NDFIT,NDEL) XSQ=XSQ+DPAR(NDEL)**2 ELSE WRITE(LUFIT,*) ' ***** NEXT LINE NOT USED IN FIT' NRJ=NRJ+1 ENDIF IF(ABS(AFRQ).GT.BIGF) AFRQ=SIGN(BIGF,AFRQ) IF(ABS(ADIF).GT.BIGD) ADIF=SIGN(BIGD,ADIF) IF (IBLND.GT.0) THEN CALL QNFMT(NQN,QNUM,AQNUM) WRITE (LUFIT,640) LINE,AQNUM,XFRQ,AFRQ,ADIF,XERR,CERR LBLND=LINE+1 ELSE LBLND=LINE C ..DO UNTIL all elements of blend are printed 290 CALL FRQDAT(LBLND,IBLND,XFRQ,XWT,XERR,QNUM) CALL DNUGET(I0,NPAR,FRQ,LBLND,DPAR) DIF=XFRQ-FRQ IF(ABS(FRQ).GT.BIGF) FRQ=SIGN(BIGF,FRQ) IF(ABS(DIF).GT.BIGD) DIF=SIGN(BIGD,DIF) CALL QNFMT(NQN,QNUM,AQNUM) WRITE (LUFIT,640) LBLND,AQNUM,XFRQ,FRQ,DIF,XERR,CERR, : AFRQ,ADIF,XWT LBLND=LBLND+1 IF (IBLND.EQ.0) GO TO 290 C ..END DO ENDIF LINE=LBLND IF (LINE.LE.NLINE) GO TO 330 C .. END loop over lines IF (NRJ.GT.0) THEN WRITE(LUFIT,700) NRJ WRITE(*,700) NRJ ENDIF IF(XSQL.GT.XSQOLD*1.1) THEN C .. fit diverging STFAC=STFAC*0.5 WRITE(*,'(A,G12.4)') ' fit diverging, factor = ',STFAC LV=VAR K=0 DO 335 I=1,NPARTOT OLDIF(I)=OLDIF(I)*0.5 PAR(I)=PAR(I)-OLDIF(I) IF(IDPAR(I).GE.0) THEN K=K+1 PARD(K)=PAR(I) N=NDEL-K CALL DCOPY(N,HEAP(LV),NPAR,FIT(K,K),I1) LV=LV+NDEL ENDIF 335 CONTINUE GO TO 461 ENDIF XSQOLD=XSQL*XSQFAC C C normalize parameters C DO 350 K=1,NPAR SCALE=ONE/SQRT( DDOT(K,FIT(K,1),NDFIT,FIT(K,1),NDFIT) ) IF(FIT(K,K).LT.0) SCALE=-SCALE DPAR(K)=SCALE CALL DSCAL(K,SCALE,FIT(K,1),NDFIT) IF (FIT(K,K).LT.ABS(THRSH)) THEN WRITE (LUFIT,660) K,(J,FIT(K,J),J=1,K) CALL DCOPY(K,ZERO,I0,FIT(K,1),NDFIT) CALL DCOPY(NDEL,ZERO,I0,ERPAR,I1) ERPAR(K)=ONE DPAR(K)=EPS CALL JELIM(FIT,ERPAR,NDFIT,NDEL) XSQ=XSQ+ERPAR(NDEL)**2 ENDIF 350 CONTINUE WRITE (LUFIT,650) (I,FIT(I,I),I=1,NPAR) C C invert lower triangle C CALL DTRI(FIT,NDFIT,NPAR,I1) C C get estimated errors and print parameters C WRITE (LUFIT,670) LV=VAR LB=FITBGN K=0 DO 460 I=1,NPARTOT IF(IDPAR(I).GE.0) THEN K=K+1 II=I N=NDEL-K DIF=FIT(NDEL,K)*DPAR(K)*STFAC SCALE=DPAR(K)*PARFAC ERPAR(I)=ABS(SCALE)*SQRT(DDOT(N,FIT(K,K),I1,FIT(K,K),I1)) CALL DSCAL(N,SCALE,FIT(K,K),I1) CALL DCOPY(N,FIT(K,K),I1,HEAP(LV),NPAR) HEAP(LB)=HEAP(LB)-DIF PARBAS=PAR(I) PAR(I)=PAR(I)+DIF PARD(K)=PAR(I) LV=LV+NDEL LB=LB+1 ELSE DIF=OLDIF(II) SCALE=PAR(I)/PARBAS DIF=DIF*SCALE ERPAR(I)=ERPAR(II)*ABS(SCALE) PAR(I)=PAR(I)+DIF ENDIF OLDIF(I)=DIF CALL PARER(PAR(I),ERPAR(I),DIF,PARE) WRITE (LUFIT,680) K,IDPAR(I),PARLBL(I),PARE 460 CONTINUE 461 IF (NF.GT.NFIR) THEN AVGMW=AVGMW/(NF-NFIR) XSQMW=SQRT(XSQMW/(NF-NFIR)) ENDIF IF (NFIR.GT.0 ) THEN AVGIR=AVGIR/NFIR XSQIR=SQRT(XSQIR/NFIR) ENDIF WRITE (LUFIT,'(A,G15.5,A,G15.5)') ' MICROWAVE AVG = ',AVGMW, : ' MHz, IR AVG =',AVGIR WRITE (LUFIT,'(A,G15.5,A,G15.5)') ' MICROWAVE RMS = ',XSQMW, : ' MHz, IR RMS =',XSQIR DENF=ONE/SQRT(FLOAT(NF)) XSQ=DENF*SQRT(XSQ) XSQL=DENF*SQRT(XSQL) ITR=ITR+1 WRITE (LUFIT,720) ITR,XSQL,XSQ WRITE (*,720) ITR,XSQL,XSQ IF (ABS(XSQ-XSQL).LT.0.000001*XSQL) ITR=NITR IF (ITR.LT.NITR) GO TO 150 C C end of iteration C C************************************************************************ 465 IF(ITR.EQ.0) THEN WRITE(*,*) ' output files not updated' STOP ENDIF REWIND LUBAK READ (LUBAK,520,END=999) CARD CALL CHTIME(CARD(51:80)) WRITE(*,521) CARD,'FIT COMPLETE' C C compute correlation matrix C CALL PRCORR(LUFIT,NPAR,FIT,NDFIT,ERPAR,IDPAR) WRITE(LUFIT,520) CARD CALL OPENLU(LUFIT,CFIL(LFIT),'C') C save results in output files CALL OPENLU(LUPAR,CFIL(LPAR),'W') CALL OPENLU(LUVAR,CFIL(LVAR),'W') WRITE (LUPAR,520) CARD WRITE (LUVAR,520) CARD READ(LUBAK,520,END=999) CARD WRITE (LUPAR,530) LIMPAR,LIMLIN,NITR,NXPAR,THRSH,XERRMX,PARFAC, : FQFACQ WRITE (LUVAR,530) NPARTOT,NLINE,NITR,NXPAR,THRSH,XERRMX,PARFAC, : FQFACQ DO 490 I=1,NOPTN C save option cards READ (LUBAK,520,END=999) CARD WRITE (LUPAR,520) CARD WRITE (LUVAR,520) CARD 490 CONTINUE CALL OPENLU(LUBAK,CFIL(LBAK),'C') WRITE (LUPAR,730) (IDPAR(I),PAR(I),ERP(I),I=1,NPARTOT) WRITE (LUVAR,730) (IDPAR(I),PAR(I),ERPAR(I),I=1,NPARTOT) CALL OPENLU(LUPAR,CFIL(LPAR),'C') CALL OPENLU(LUBIN,CFIL(LBIN),'WB') CALL PUTVAR(LUVAR,LUBIN,NPAR,IDPAR,HEAP(VAR),ERPAR) CALL OPENLU(LUVAR,CFIL(LVAR),'C') CALL OPENLU(LUBIN,CFIL(LBIN),'C') STOP C 999 WRITE (*,*) ' END FILE ON .PAR BEFORE CARD 3' STOP C 520 FORMAT (A) 521 FORMAT (1X,A) 530 FORMAT (4I5,3E15.4,F13.10) 540 FORMAT (' LINES REQUESTED =',I5,' NUMBER OF PARAMETERS=' :,I4,' NUMBER OF ITERATIONS=',I3/,' THRESHOLD FOR PARAMETER NORM=' :,E11.4,' max (OBS-CALC)/ERROR =',E10.4) 545 FORMAT (' NUMBER OF PARAMETERS EXCLUDED IF INDEX < 0 ',I5) 550 FORMAT (/,24X,' PARAMETERS - A.PRIORI ERROR ') 620 FORMAT (' WARNING .. SIZE OF BLOCK',I5,' IS',I5, : ' AND EXCEEDS DIMENSIONS') 630 FORMAT (/,' ####',40X,' EXP.FREQ. - CALC.FREQ. - DIFF. ' :, '- EXP.ERR.- EST.ERR.-AVG. CALC.FREQ. - DIFF. - WT.') 640 FORMAT (I5,1H:,A,2F14.5,3F10.5,F15.5,F10.5,F7.4) 650 FORMAT (' NORMALIZED DIAGONAL:',11(/,6(I5,E15.5))) 660 FORMAT (' PARAMETER',I3,' NOT DETERMINABLE AND IS FIXED,',10X, : 'COLUMN INFORMATION:',7(/,20X,5(I5,E15.5))) 670 FORMAT (/,31X,' NEW PARAMETER - EST. ERROR - CHANGE THIS ITER.') 680 FORMAT (I4,F16.0,3X,A,1X,A) 700 FORMAT (I5,' Lines rejected from fit') 720 FORMAT (' END OF ITERATION',I3,' OLD, NEW RMS ERROR=',2G15.5) 730 FORMAT (F16.0,E23.15,E15.8) C END SUBROUTINE JELIM(T,VEC,NDM,NZ) C C do orthogonal transformations to rotate vector VEC into matrix T C C NDM = dimensioned column length of T C IBGN = first element of VEC C NZ = length of VEC C REAL*8 T(*),VEC(*),S,C INTEGER NDM,NZ,N,NDM1,K,KD,KI,KDI,I1,NL C EMA T,VEC PARAMETER (I1=1) N=NZ-1 NDM1=NDM+1 C .. KD is vector index of T(K,K) KD=1 NL=NZ DO 320 K=1,N NL=NL-1 KI=K+1 KDI=KD+1 IF (ABS(T(KD)).LT.ABS(VEC(K))) THEN C=T(KD)/VEC(K) S=C*C IF(S.GT.1.E-30) THEN S=SQRT(1 + S) T(KD)=VEC(K)*S S=1/S C=C*S CALL DROT(NL,T(KDI),I1,VEC(KI),I1,C,S) ELSE T(KD)=VEC(K) CALL DSWAP(NL,T(KDI),I1,VEC(KI),I1) ENDIF ELSE S=VEC(K)/T(KD) C=S*S IF(C.GT.1.E-30) THEN C=SQRT(1 + C) T(KD)=T(KD)*C C= 1 / C S=S*C CALL DROT(NL,T(KDI),I1,VEC(KI),I1,C,S) ENDIF ENDIF 320 KD=KD+NDM1 RETURN END C SUBROUTINE QNFMT(NQN,QNUM,AQNUM) INTEGER NQN,QNUM(12) CHARACTER*36 AQNUM C formats quantum numbers for output C NQN = number of quanta C QNUM = vector of quanta C AQNUM = string of quanta in character form C INTEGER I,K WRITE(AQNUM,'(12I3)') QNUM DO 10 I=NQN+1,12 IF(QNUM(I).EQ.0) THEN K=I*3-2 AQNUM(K:K+2)=' ' ENDIF 10 CONTINUE RETURN END SUBROUTINE PARER(PAR,ERRX,DIF,PTMP) REAL*8 PAR,ERRX,DIF CHARACTER*(*) PTMP C C sets up special format for parameters and errors C PAR = parameter value C ERRX = parameter error C DIF = parameter change C PTMP = output string for printing C LOGICAL EFIELD CHARACTER FMT*30,CHEXP*4,CHND*1 REAL*8 APAR,AERR,ADIF,ATEN INTEGER IP,INT,IE,ID,MSDX,LSD,MSDXM,MSDXL C EMA PAR,ERRX,DIF APAR=PAR AERR=ERRX ADIF=DIF C compute exponent fields ID=INT( LOG10(ABS(ADIF)+1.5E-38)-100. ) + 100 IP=INT( LOG10(ABS(APAR)+1.5E-38)-100. ) + 100 ATEN= LOG10(ABS(AERR)+1.5E-38) IE=INT( ATEN-100. ) + 100 MSDX=MAX(IP,IE,ID) LSD=MIN(0,INT( ATEN-102.5 ) + 100 ) C check for too many digits IF(IP.GE.0) IP=IP-MIN(IP,4) LSD=MAX( MAX(IP,IE,ID)-9 ,LSD) C check for too many digits to right of dp or small numbers EFIELD=LSD.LT.-9.OR.LSD.GT.0.OR.MSDX.LE.-4 IF(MSDX.EQ.0) EFIELD=.FALSE. IF(EFIELD) THEN C .. E format ATEN=10.**(-MSDX) IF (MSDX.LT.0) THEN CHEXP='E-00' MSDX=-MSDX ELSE CHEXP='E+00' ENDIF MSDXM=MSDX/10 MSDXL=MSDX IF(MSDXM.GT.0) THEN MSDXL=MSDXL-10*MSDXM CHEXP(3:3)=CHAR(ICHAR('0')+MSDXM) ENDIF CHEXP(4:4)=CHAR(ICHAR('0')+MSDXL) APAR=APAR*ATEN AERR=AERR*ATEN ADIF=ADIF*ATEN ELSE C .. F format MSDX=0 CHEXP=' ' ENDIF CHND(1:1)=CHAR(ICHAR('0')+MIN(MSDX-LSD,9) ) C ..FMT= (F15.n,'E+mm',2(1X,F11.n,'E+mm') ) FMT='(F15.'//CHND//',A,2(1X,F11.'//CHND//',A))' WRITE(PTMP ,FMT,ERR=100) APAR,CHEXP,AERR,CHEXP,ADIF,CHEXP RETURN 100 WRITE(PTMP,'(3E11.5)') PAR,ERRX,DIF RETURN END C SUBROUTINE LINEIN(LUIN,NLINE,IQNFMT,MXQN) INTEGER LUIN,NLINE,IQNFMT,MXQN C get lines from input and stores them C C LUIN= unit for finding lines C NLINE = number of lines C IQNFMT= qunatum number format for line input C MXQN = largest quantum number C REAL*8 XFRQ,XWT,XERR INTEGER IQNUM(12),NQNU,NQNL,NQNT(12),I,K,MXLINE,IPACE,I0,IQF, : GETLIN,DEFLIN I0=0 MXLINE=NLINE NLINE=0 MXQN=1 NQNU=DEFLIN(IQNFMT,NQNT) NQNL=NQNU IF (NQNT(NQNL).EQ.0) NQNU=1 NQNL=NQNL+NQNU IPACE=100 C loop for reading lines C ***************************************************************** DO 140 I=1,MXLINE IF(GETLIN(LUIN,NQNT,IQNUM,XFRQ,XERR,XWT).NE.0) RETURN K=NQNU IQF=IQNUM(NQNU) IF(IQF.EQ.-1) THEN IF(NQNT(K).GT.0) THEN K=NQNT(K) IQF=-IQNUM(K) IF(IQF.EQ.0) IQF=-1 ENDIF IQNUM(NQNU)=IQF ENDIF IF(IQF.LT.0) IQF=-IQF IF(MXQN.LT.IQF) MXQN=IQF K=NQNL IQF=IQNUM(NQNL) IF(IQF.EQ.-1) THEN IF(NQNT(K).GT.0) THEN K=NQNT(K) IQF=-IQNUM(K) IF(IQF.EQ.0) IQF=-1 ENDIF IQNUM(NQNU)=IQF ENDIF IF(IQF.LT.0) IQF=-IQF IF(MXQN.LT.IQF) MXQN=IQF NLINE=I CALL SAVDAT(NLINE,I0,I0,I0,I0,I0,I0, : I0,XFRQ,XWT,XERR,IQNUM) IF (IPACE.LE.NLINE) THEN IPACE=IPACE+100 WRITE(*,*) 'Reading Line ',NLINE ENDIF 140 CONTINUE RETURN END SUBROUTINE LINEIX(LU,NLINE,NBLKPF,IQNFMT,HEAD) INTEGER LU,NLINE,NBLKPF,IQNFMT,HEAD C get lines from input and stores them C C LU = unit for printout of lines ( if > 0 ) C NLINE = number of lines C NBLKPF= number of blocks per F C IQNFMT= qunatum number format for line input C INTEGER MXBLK,NBLK,ORGBLK PARAMETER (MXBLK=2048) INTEGER*2 PRVBLK(MXBLK) LOGICAL BLEND,DONE REAL*8 XFRQ,XFRQEX,XWT,XERR,XWT1,SUM INTEGER NQN,I,IPOS,NBLND,IPACE,MXLINE,NREAD,IBLKU,IBLKL,INDXU, : IQNUM(12),INDXL,L,J,IPREV,IBLND,K,LINKU,LINKL,LINKX,LNLINK MXLINE=NLINE NLINE=0 NBLK=0 C set head of linked list HEAD=0 DO 1 I=1,MXBLK 1 PRVBLK(I)=0 NQN=MOD(IQNFMT,10) I=MOD(IQNFMT/100,5) IF (I.GE.NQN) THEN IPOS=1 ELSE IPOS=NQN ENDIF IF (LU.GT.0) WRITE (LU,570) XFRQEX=1.D38 NBLND=0 IPACE=50 NREAD=0 LINKX=0 ORGBLK=0 C loop for converting lines C ***************************************************************** 140 NREAD=NREAD+1 DONE=NREAD.GT.MXLINE IF (DONE) THEN IPACE=NLINE BLEND=.FALSE. ELSE CALL FRQDAT(NREAD,K,XFRQ,XWT,XERR,IQNUM) IF (XWT.LT.1.E-20) XWT=1 BLEND= ABS(XFRQ-XFRQEX).LT.0.000001 C .. find blocks and index for upper and lower states CALL GETBLK(IBLKU,INDXU,IQNUM ,NBLKPF,IPOS,NQN) CALL GETBLK(IBLKL,INDXL,IQNUM(NQN+1),NBLKPF,IPOS,NQN) IF (IBLKU.EQ.0) THEN C .. print out bad line and try for next L=ABS(LU) WRITE(L,'(A,12I3,F15.5,2F10.5)') ' Bad Line:',IQNUM, : XFRQ,XERR,XWT WRITE(*,'(A,12I3,F15.5,2F10.5)') ' Bad Line:',IQNUM, : XFRQ,XERR,XWT GO TO 140 ENDIF ENDIF C .. let the user know something is happening IF (IPACE.LE.NLINE) THEN IPACE=IPACE+50 WRITE(*,*) 'Converting Line ',NLINE ENDIF IF (BLEND) THEN NBLND=NBLND+1 NLINE=NLINE+1 J=NLINE-NBLND IF (LU.GT.0) WRITE (LU,600) NLINE,IBLKU,INDXU,IBLKL, : INDXL,IQNUM,XFRQ,XERR,XWT, : ' Line Blended with ',J XWT=XWT/XWT1 SUM=SUM+XWT ELSE IF (NBLND.GT.0) THEN C go back and normalize weights for blend C IBLND=1 (no blend) =-1 (last of blend) =0 (blend) C IPREV=NLINE IBLND=-1 DO 10 K=0,NBLND CALL SAVXWT(IPREV,SUM,IBLND) IBLND=0 10 IPREV=IPREV-1 ENDIF IF (DONE) GO TO 200 NBLND=0 NLINE=NLINE+1 IF (LU.GT.0) WRITE (LU,600) NLINE,IBLKU,INDXU,IBLKL, : INDXL,IQNUM,XFRQ,XERR,XWT XWT1=XWT SUM=1. XWT=1. ENDIF IBLND=1 C .. set up links for calculating in order of block IF(IBLKU.LE.IBLKL) THEN LINKU=LNLINK(HEAD,PRVBLK,MXBLK,IBLKU,NLINE,LINKX,0) LINKL=LNLINK(HEAD,PRVBLK,MXBLK,IBLKL,-NLINE,LINKU,0) IF(NBLK.LT.IBLKL) NBLK=IBLKL ELSE LINKL=LNLINK(HEAD,PRVBLK,MXBLK,IBLKL,-NLINE,LINKX,0) LINKU=LNLINK(HEAD,PRVBLK,MXBLK,IBLKU,NLINE,LINKL,0) IF(NBLK.LT.IBLKU) NBLK=IBLKU ENDIF CALL SAVDAT(NLINE,LINKU,LINKL,IBLKU,INDXU,IBLKL,INDXL, : IBLND,XFRQ,XWT,XERR,IQNUM) XFRQEX=XFRQ GO TO 140 C .. while (NBLK.GT.MXBLK) 200 IF(NBLK.LE.MXBLK) RETURN NBLK=MXBLK-1 ORGBLK=ORGBLK+NBLK LINKU=0 DO 201 J=1,NBLK IF(PRVBLK(J).NE.0) THEN LINKU=PRVBLK(J) PRVBLK(J)=0 ENDIF 201 CONTINUE PRVBLK(MXBLK)=0 IF(LINKU.NE.HEAD) THEN ORGBLK=ORGBLK-1 C set PRVBLK(1)=LINKU, LINKU=beginning, old link=0 J=0 IBLKU=1 LINKU=LNLINK(J,PRVBLK,MXBLK,IBLKU,LINKU,LINKU,-1) ELSE HEAD=0 ENDIF NBLK=0 C ..while(LINKU.NE.0) 210 IF(LINKU.NE.0) THEN LINKL=LINKU CALL GETDBK(LINKU,J,IBLKU,INDXU,J,J) IBLKU=IBLKU-ORGBLK IF(NBLK.LT.IBLKU) NBLK=IBLKU LINKX=LNLINK(HEAD,PRVBLK,MXBLK,IBLKU,LINKL,LINKU,-1) GO TO 210 ENDIF GO TO 200 C 570 FORMAT (/,' LINE,BLKU,INDXU,BLKL,INDXL,QUANTUM NUMBERS',19X, : 'ENERGY EXP. ERROR WEIGHTS') 600 FORMAT (1X,5I4,1H:,12I3,F15.4,2F10.4,A,I4) END