C$large c$debug C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C LFITD - PROGRAM FOR FULLY DIAGONALISING TREATMENT OF K-DOUBLING IN C v=1 E STATE IN A C3v SYMMETRIC TOP C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C The calculation employs factorization of the 2*(2J+1) matrix into an C A matrix and one of two degenerate E matrices, each ca 3J/2 in size. C C C - J UP TO ca 80 C - UP TO 'MAXLIN' TRANSITIONS C - WEIGHTED FIT C - OPTIONAL CORRELATION COEFFICIENTS AND CONTRIBUTIONS OF C INDIVIDUAL CONSTANTS TO TRANSITION FREQUENCIES (EVALUATED C FROM DERIVATIVES) C - FIT WITH ZERO DEGREES OF FREEDOM (ie. N constants to N C frequencies) IS ALLOWED, BUT THE USER IS WARNED BY RINGING C THE BELL C C see: K.M.T.Yamada,M.Bester,M.Tanimoto,G.Winnewisser - J.Mol.Spectrosc. C 126,118(1987) for details on setting up the Hamiltonian C C ver. 13.III.97 ----- Zbigniew KISIEL ----- C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C_________________________/-------------------------------------------------- C C Modification history: C c 9.12.95: minor cosmetics c 29.07.96: more comments and parameterization c 13.03.97: printout also of K and l C_____________________________________________________________________________ C C This program is standard F77 and compiles with MSFortran 5.0. To C transfer to a different Fortran comment out the $LARGE metacommand. C_____________________________________________________________________________ C C DATA INPUT: C C The program provides a conversational way of inputting and editing C the spectroscopic data which is hoped to be self explanatory. The C following assumptions are made: C C - 1/0 convention for YES/NO (usually the answer is read in in fixed C format so a suffices for NO) C - constants are to supplied in the form: 1,constant value C or: 0,constant value C depending on whether the constant is to be fitted or not C C - only DeltaJ=1 transitions are supported which are to be specified C using the two quantum numbers of the lower state: J and Kl-1. C If two different A-type transitions for one value of Kl-1 (such as C the Kl-1=0 doublets) need to be specified then a negative value of C J declares the lower frequency component. The quantum numbers C are to be followed by the frequency and its error eg. C 3,-2,8669.6161,0.01 C C - termination of frequency input loop (including insert and C add modes)is by inputting 0's for every parameter eg. C 0,0,0,0,0,0,0,0 c c - dump of information associated with diagonalization is possible c by following the number of lines in the input file with 1 written c in I5 format. The output is made to file DUMP.K C C It is recommended that on first use of LFITD on a new problem the C data be input in conversational mode and saved with option 7 of the C modification menu prior to any calculations. Data files created C in this way have a self explanatory structure and can then be modifed C with either the modification options of the program or with an external C editor. C C_____________________________________________________________________________ C C The available constants are as follows: C C A, B - rotational constants C DJ, DJK, DK, HJ, HJK, HKJ, HK - centrifugal distortion constants C nJ, nK, nJJ, nJK, nKK - centrifugal distortion on Coriolis resonance C qv, qJ - (2,1) q(+) l-type doubling and resonance constants C rv - (1,-2) r-type interaction constant C q3 - DeltaK=+-3 interaction constant equal to tau.xxxz/4 C h6 - DeltaK=+-6 interaction constant C t - (K=2,l=-1) <-> (K=-2,l=1) interaction term, or Kl-1=-3 splitting C term C C C I - number of transitions C N - number of energy levels to be used in fit (if errors are present C in line declarations this may be less than 2*I) C K1 - number of fitted constants C FR - frequencies of transitions where C (i,1) observed frequencies C (i,2) calculated frequencies C (i,3) obs-calc C ERROR - estimated experimental error in frequency C WEIGHT - weight equal to 1/error**2 C NJA - information on energy levels where C (i,1) J quantum number C (i,2) number of submatrix to which the level belongs =1,2 for C A, E resp. C (i,3) position of the energy level among eigenvalues of the relevant C submatrix in ascending order of magnitude: two eigenvalues for C each value of Kl-1 are possible and are coded into NJA(i,3) C with a+1000*b C (i,4) number of transition associated with the level, this is negative C for the lower level C IA - pointers to which constants are to be fitted (=1, otherwise 0) C LINFIT - pointers to which lines are to be fitted (+1, otherwise 0) C JK - quantum numbers for transitions where C (i,1), (i,2), (i,3) are J, K-1, K+1 for upper level resp. C (i,4), (i,5), (i,6) are J, K-1, K+1 for lower level resp. C IPR - default number of iterations C_____________________________________________________________________________ C C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (MAXLIN=600,MAXTWO=2*MAXLIN,MAXCON=35,MAXDIM=140) COMMON /BLKH/H/BLKR/R COMMON /BLK1/JK,FR,ERROR,WEIGHT,LINFIT,A,IA,I, 1 IOUT/BLK2/NJA,B/BLK3/C,D,DD/TBLOCK/COM,T DIMENSION A(MAXCON),B(MAXCON),C(MAXCON),Z(MAXCON),IA(MAXCON), * H(MAXDIM,MAXDIM),R(MAXDIM,MAXDIM),DD(MAXCON,MAXCON) INTEGER*4 NJA(MAXTWO,4),MNJA,MNJA1 DIMENSION FR(MAXLIN,3),JK(MAXLIN,6),LINFIT(MAXLIN) REAL*4 D(MAXLIN,MAXCON),EE,ERROR(MAXLIN),WEIGHT(MAXLIN) CHARACTER*6 T(MAXCON),TT(MAXCON) CHARACTER*1 COM(79) DATA TT/'A =','B =',' DJ =',' DJK =',' DK =', 1 'HJ =','HJK =','HKJ =','HK =','Azeta=', 1 ' nJ =',' nK =',' nJJ =',' nJK =',' nKK =', 1 'qv =','qj =','qk =','qjj =','qjk =','qkk =', 1 'qjjj =',' rv =',' q3 =',' h6 =', 1 ' t =','LJ =','LJJK =', 1 'LJK =','LKKJ =',' =',' =',' =', 1 ' =',' ='/ NCONST=30 NQNS=2 DO 1234 L=1,NCONST T(L)=TT(L) 1234 CONTINUE IOUT=0 C IFLAG=0 CALL DATAIN(NCONST,NQNS,idump) 400 CALL MODIF(N,IFLAG,NCONST,NQNS) IF(IOUT.EQ.1)GOTO 41 C I1=0 S=1.D50 IPR=3 if(idump.eq.1)then open(2,file='dump.k',status='unknown') do 770 k=1,n-1 m1=iabs(nja(k,4)) write(2,771)nja(k,1),jk(m1,2),(nja(k,l),l=2,4) 771 format(' J,Kl-1: ',2i3,' Nmat,Nposmat:',i3,i7,' line:',i4) 770 continue endif C C C START ITERATION C C 14 WRITE(*,6000) 6000 FORMAT(/' **** CALCULATION IN PROGRESS, CURRENTLY WORKING ON J', 1 ' ='/) I1=I1+1 S1=S DO 15 K=1,I FR(K,2)=0.D0 FR(K,3)=0.D0 DO 15 L=1,NCONST D(K,L)=0.0 15 CONTINUE L1=1 M4=0 WRITE(*,801)NJA(1,1) C C Go through all sorted energy levels C DO 25 K=2,N KK=K-1 IF(NJA(K,1).EQ.NJA(KK,1))GOTO 250 WRITE(*,801)NJA(KK,1) 801 FORMAT(I4,$) 250 IF(NJA(K,1).EQ.NJA(KK,1).AND.NJA(K,2).EQ.NJA(KK,2))GOTO 25 L2=K-1 N1=NJA(KK,2) J1=NJA(KK,1) C C Set up and diagonalize the appropriate submatrix for each group of C common levels C CALL HAMIL(M,J1,N1,A) C--- WRITE(*,1802)J1,N1,M C---1802 FORMAT(' J, MAT.NO., MAT.SIZE: ',3I8) IEGEN=0 if(idump.eq.1)then write(2,1855)j1,n1,m,m 1855 format(' J=',i4,' type ',i1,i6,'x',i3, * ' (first diagonal elements of H, then eigenvalues)') endif if(idump.eq.1)then write(2,1856)(h(l,l),l=1,m) endif c--- CALL HDIAG(M,IEGEN,K1) call jack(m) if(idump.eq.1)then write(2,1856)(h(l,l),l=1,m) 1856 format(8f13.2) endif DO 17 L=L1,L2 M1=NJA(L,4) M2=ABS(M1) MNJA=NJA(L,3) IF(MNJA.GT.1000)THEN MNJA=MNJA/1000 MNJA1=NJA(L,3)-1000*MNJA M=MNJA NJA(L,3)=M IF(JK(M2,1).LT.0)THEN IF(H(MNJA1,MNJA1).LT.H(MNJA,MNJA))M=MNJA1 NJA(L,3)=M ELSE IF(H(MNJA1,MNJA1).GT.H(MNJA,MNJA))M=MNJA1 NJA(L,3)=M ENDIF ELSE M=NJA(L,3) ENDIF FR(M2,2)=FR(M2,2)+SIGN(1,M1)*H(M,M) 17 CONTINUE C C FIND THE DERIVATIVES - since transformation to the space where the Wang C submatrix is diagonal has already been found (eigenvectors), the deri- C vatives for energy levels in this submatrix with respect to a given C constant are obtained by the same similarity transformation of matrix C set up with all other constants set to zero and the required constant C set to 1 (array B) C K1=0 DO 24 J=1,NCONST IF(IA(J).EQ.1)GOTO 18 GOTO 24 18 K1=K1+1 DO 19 L=1,NCONST B(L)=0.D0 19 CONTINUE B(J)=1.D0 CALL HAMIL(M3,J1,N1,B) DO 124 L=L1,L2 M=NJA(L,3) IF(L.NE.L1.AND.M.EQ.M4)GOTO 23 EE=0.0 DO 22 M1=1,M3 DO 122 M2=M1,M3 N2=2 IF(M1.EQ.M2)N2=1 EE=EE+N2*H(M1,M2)*R(M1,M)*R(M2,M) 122 CONTINUE 22 CONTINUE 23 M1=NJA(L,4) M2=ABS(M1) D(M2,K1)=D(M2,K1)+SIGN(1,M1)*EE M4=M 124 CONTINUE 24 CONTINUE L1=K 25 CONTINUE if(idump.eq.1)close(2) C C FIND NEW CONSTANTS - inverse of D'D is evaluated by diagonalization using C inv A = U * inv L * U' where U and L contain the eigenvectors and the C eigenvalues of A. C C To improve conditioning D'D is rescaled to D'D.ij/(D'D.ii D'Djj)**1/2 C (scaling factors in Z) C To save space D'D is placed and diagonalized in the top lh. corner C of H, D'Y is placed in (K1+1)'st column of H and inv D'D is placed C in array DD equivalenced to the rh. side of H C S=0.D0 SW=0.D0 L=0 DO 26 J=1,I FR(J,3)=DSIGN(FR(J,1),FR(J,2))-FR(J,2) IF(LINFIT(J).NE.1)GOTO 26 S=S+FR(J,3)**2 SW=SW+(FR(J,3)/ERROR(J))**2 L=L+1 26 CONTINUE N139=L-K1 IF(N139.LT.0)GOTO 48 IF(N139.GT.0)GOTO 810 WRITE(*,89)L,K1,CHAR(7) SDEV=DSQRT(S) SWDEV=DSQRT(SW) GOTO 811 810 SDEV=DSQRT(S/N139) SWDEV=DSQRT(SW/N139) 811 WRITE(*,78)I1,SDEV,SWDEV 78 FORMAT(//' ITERATION STEP NO.',I2,11X,'Deviation ', 1 'of fit =',F14.6/24X,'Weighted deviation of fit =',F14.6) DO 30 J=1,K1 DO 28 K=J,K1 E=0.D0 DO 27 L=1,I IF(LINFIT(L).NE.1)GOTO 27 E=E+D(L,J)*D(L,K)*WEIGHT(L) 27 CONTINUE H(J,K)=E 28 CONTINUE Z(J)=SQRT(1.D0/H(J,J)) E=0.D0 DO 29 L=1,I IF(LINFIT(L).NE.1)GOTO 29 E=E+D(L,J)*FR(L,3)*WEIGHT(L) 29 CONTINUE H(J,K1+1)=E 30 CONTINUE DO 301 J=1,K1 DO 301 K=J,K1 H(J,K)=H(J,K)*Z(J)*Z(K) H(K,J)=H(J,K) 301 CONTINUE IEGEN=0 CALL HDIAG(K1,IEGEN,L3) DO 35 J=1,K1 B(J)=0.D0 DO 34 K=1,K1 E=0.D0 DO 33 L=1,K1 E=E+R(J,L)*R(K,L)/H(L,L) 33 CONTINUE E=E*Z(J)*Z(K) DD(J,K)=E B(J)=B(J)+E*H(K,K1+1) 34 CONTINUE C(J)=SQRT(DD(J,J))*SWDEV 35 CONTINUE C WRITE(*,81) 81 FORMAT(/' RESULTS - constant, change, standard deviation :'/) K=0 DO 40 J=1,NCONST IF(IA(J).EQ.1)GOTO 39 GOTO 40 39 K=K+1 A(J)=A(J)+B(K) IF(J.LE.8)THEN WRITE(*,82)T(J),A(J),B(K),C(K) ELSE WRITE(*,1082)T(J),A(J),B(K),C(K) ENDIF 82 FORMAT(1X,A6,3F23.12) 1082 FORMAT(1X,A6,F27.16,2F23.16) 40 CONTINUE C C Convergence tests and output C 403 WRITE(*,'(/'' Do you want to interrupt the fit ? '',$)') READ(*,5,ERR=403)IIFLAG IF(IIFLAG.NE.1)GOTO 6 C 404 WRITE(*,401) 401 FORMAT(1X/ * 10X,'=0 OBS-CALC DIFFERENCES =3 NEXT ITERATION'/ * 10X,'=1 MODIFICATIONS =4 EXIT LFITD'/ * 10X,'=2 DISK OUTPUT ......... ',$) READ(*,5,ERR=404)IFLAG 5 FORMAT(I5) IF(IFLAG.LT.0.OR.IFLAG.GT.4)GOTO 404 IF(IFLAG.EQ.1)GOTO 400 IF(IFLAG.EQ.3)GOTO 6 IF(IFLAG.EQ.4)GOTO 99 C if(iflag.eq.0)then WRITE(*,'(1X)') DO 7 J=1,I IF(LINFIT(J).NE.1)GOTO 8 WRITE(*,9)J,(JK(J,L),L=1,nqns),FR(J,1),FR(J,3),ERROR(J) GOTO 7 8 WRITE(*,10)J,(JK(J,L),L=1,nqns),FR(J,1),FR(J,3),ERROR(J) 7 CONTINUE 9 FORMAT(1X,I3,'. ',2I10,F18.4,2F15.4) 10 FORMAT(1X,I3,'. ',2I10,F18.4,F15.4,'EXCL',F11.4) goto 404 endif c if(iflag.eq.2)then IOUT=2 CALL DATOUT(IFLAG,SDEV,SWDEV,NCONST,NQNS,idump) goto 404 endif C 6 iflag=1 E=(S1-S)/S IF(DABS(E).LT.0.001D0)GOTO 41 IF(E.LT.-0.1D0)GOTO 41 IF(I1.LT.IPR)GOTO 14 1900 WRITE(*,450) 450 FORMAT(/' MORE CYCLES ? ',$) READ(*,'(I5)',ERR=1900)L IF(L.NE.1)GOTO 41 I1=0 GOTO 14 C 41 CALL DATOUT(IFLAG,SDEV,SWDEV,NCONST,NQNS,idump) C IF(IFLAG.EQ.1)GOTO 400 GOTO 99 48 WRITE(*,89)L,K1,' - ITERATION STOPPPED' 89 FORMAT(//' ONLY',I3,' MEASURED TRANSITIONS FOR',I3, 1 ' CONSTANTS',A) GOTO 41 C 99 STOP END C C_____________________________________________________________________________ C SUBROUTINE DATAIN(NCONST,NQNS,idump) C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (MAXLIN=600,MAXCON=35) COMMON /BLK1/JK,FR,ERROR,WEIGHT,LINFIT,A,IA,I,IOUT COMMON /TBLOCK/COM,T DIMENSION JK(MAXLIN,6),FR(MAXLIN,3),A(MAXCON),IA(MAXCON), * LINFIT(MAXLIN) CHARACTER*6 T(MAXCON),cdummy CHARACTER*1 COM(79) CHARACTER*25 FILNAM REAL*4 ERROR(MAXLIN),WEIGHT(MAXLIN) C WRITE(*,3344) 3344 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | LFITD - FITTING OF ROTATIONAL SPECTRUM', * ' IN v=1 E STATE OF A C3v SYMMETRIC',T79,'|'/ * ' | TOP (diagonalization of E,A factorised ', * 'complete eigenmatrix)',T79,'|'/ * ' |',76(1H_),'|'/' version 13.III.1997',T64,'Zbigniew KISIEL'/) 1900 WRITE(*,12) 12 FORMAT(//' DATA TO BE CREATED (1), OR READ FROM DISK (2)', 1 ': ',$) READ(*,14,ERR=1900)IFLAG 14 FORMAT(2I5) C 31 IF(IFLAG.EQ.2)GOTO 15 1901 WRITE(*,16) 16 FORMAT(/' COMMENT :') READ(*,17,ERR=1901)COM 17 FORMAT(79A1) WRITE(*,18) 18 FORMAT(/' FITTING PARAMETER, INITIAL VALUE OF CONSTANT : ' 1 /) NFIN=NCONST DO 19 L=1,NFIN 1903 WRITE(*,20)T(L) 20 FORMAT(1X,A6,2X,$) READ(*,*,ERR=1902)IA(L),A(L) GOTO 19 1902 WRITE(*,'(1H+,A1)')CHAR(7) GOTO 1903 19 CONTINUE I=0 WRITE(*,22) 22 FORMAT(/' LINE ASSIGNMENTS, FREQUENCIES AND ERRORS ', 1 '(Separated by commas) :'/) 24 I=I+1 1906 READ(*,*,ERR=1904)(JK(I,K),K=1,6),FR(I,1),ERROR(I) IF(ERROR(I).EQ.0.)ERROR(I)=1. WEIGHT(I)=1./ERROR(I)**2 GOTO 1905 1904 WRITE(*,'(1H+,A1)')CHAR(7) GOTO 1906 1905 IF(JK(I,1).GT.0)GOTO 24 I=I-1 DO 70 L=1,I LINFIT(L)=1 70 CONTINUE GOTO 41 C C...input from disk C 15 WRITE(*,35) 35 FORMAT(/' NAME OF DATA FILE : ',$) READ(*,'(A)',ERR=15)FILNAM OPEN(2,FILE=FILNAM,ERR=15) READ(2,3500,ERR=15)COM READ(2,14,ERR=15)I,idump READ(2,150,ERR=15)(cdummy,IA(K),A(K),K=1,NCONST) READ(2,140)((JK(K,M),M=1,NQNS),FR(K,1),ERROR(K),LINFIT(K),K=1,I) CLOSE(2) 3500 FORMAT(1X,79A1) 150 FORMAT(A6,I5,F30.19) 140 FORMAT(2I5,F20.6,F15.6,I5) DO 1020 K=1,I IF(ERROR(K).EQ.0.)ERROR(K)=1.0 WEIGHT(K)=1./ERROR(K)**2 1020 CONTINUE C 41 RETURN END C C_____________________________________________________________________________ C SUBROUTINE MODIF(N,IFLAG,NCONST,NQNS) C C Modification of data for the fit and identification and sorting of C energy levels: C C N - on output contains number of energy levels identified for fitting C IFLAG - if set to 1 then check of data, listing and reevaluation of NJA C are skipped CNCONST - total number of constants that can be fitted with current version C of the program C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (MAXLIN=600,MAXTWO=2*MAXLIN,MAXCON=35) COMMON /BLK1/JK,FR,ERROR,WEIGHT,LINFIT,A,IA,I, 1 IOUT/BLK2/NJA,B/TBLOCK/COM,T DIMENSION JK(MAXLIN,6),FR(MAXLIN,3),A(MAXCON),IA(MAXCON), * B(MAXCON),LINFIT(MAXLIN) CHARACTER*6 T(MAXCON) INTEGER*4 NJA(MAXTWO,4),MNJA INTEGER*2 IINT(6) CHARACTER*1 COM(79) REAL*4 ERROR(MAXLIN),WEIGHT(MAXLIN) IONE=1 C IF(IFLAG.EQ.1)GOTO 11 C C Check of data, evaluation of NJA (identification of energy levels) and C listing C 29 WRITE(*,30)COM WRITE(*,31) DO 32 K=1,NCONST WRITE(*,33)T(K),A(K),IA(K) 32 CONTINUE WRITE(*,34)I C N=0 DO 10 KK=1,I C C...is Kl-1 in the range -(J+1) to (J-1) ? IF(JK(KK,2).LT.-(IABS(JK(KK,1))+1).OR.JK(KK,2).GT. 1 IABS(JK(KK,1))-1)GOTO 9 C C...is Kl-1 of the lower A-type component specified by negative J a multiple c of 3 ? IF(JK(KK,1).LT.0)THEN J=IABS(JK(KK,2)) IF(J-3*(J/3).NE.0)GOTO 9 ENDIF C IF(JK(KK,2)-3*(JK(KK,2)/3).EQ.0)THEN NMAT=1 ELSE NMAT=2 ENDIF C DO 611 KKK=1,2 N=N+1 NJA(N,2)=NMAT IF(KKK.EQ.1)THEN NJA(N,1)=IABS(JK(KK,1)) NJA(N,4)=-KK ELSE NJA(N,1)=IABS(JK(KK,1))+1 NJA(N,4)=KK ENDIF C C...position of level among eigenvalues of the matrix NJA(N,3)=0 J=NJA(N,1) C C...matrix A: two eigenvalues for each value of Kl-1 are possible and C are coded into NJA(i,3) with a+1000*b if(NMAT.eq.1)then L=1 KE=J IF(J-3*(J/3).EQ.0)KE=KE-1 KA=-KE KLM1=KA*L-1 IF(KLM1-3*(KLM1/3).NE.0)L=-1 C--- WRITE(*,*)J,KA,L C II=0 DO 612 K=KA,KE IF(K-3*(K/3).EQ.0)GOTO 612 II=II+1 klm1=K*L-1 IF(KLM1.EQ.JK(KK,2))THEN IF(NJA(N,3).NE.0)THEN NJA(N,3)=NJA(N,3)+1000*II C--- WRITE(*,*)NJA(N,1),NJA(N,2),NJA(N,3),NJA(N,4) GOTO 611 ELSE NJA(N,3)=II ENDIF ENDIF L=-L 612 CONTINUE endif C C...matrix E if(NMAT.eq.2)then KA=-J KE=J JJ=J-3*(J/3) IF(JJ.EQ.1)KE=J-1 IF(JJ.EQ.2)KA=KA+1 L=1 IF(JJ.EQ.0)L=-1 C--- WRITE(*,*)J,KA,L C II=0 DO 614 K=KA,KE KLM1=K*L-1 IF(KLM1-3*(KLM1/3).EQ.0)GOTO 614 II=II+1 KLM1=K*L-1 IF(KLM1.EQ.JK(KK,2))THEN NJA(N,3)=II C--- WRITE(*,*)NJA(N,1),NJA(N,2),NJA(N,3),NJA(N,4) GOTO 611 ENDIF L=-L 614 CONTINUE endif 611 CONTINUE C IF(LINFIT(KK).NE.1)GOTO 99 WRITE(*,96)KK,(JK(KK,M),M=1,NQNS),FR(KK,1),ERROR(KK) GOTO 10 99 WRITE(*,98)KK,(JK(KK,M),M=1,NQNS),FR(KK,1),ERROR(KK) GOTO 10 C C...error in quantum numbers 9 WRITE(*,97)KK,(JK(KK,M),M=1,NQNS),FR(KK,1),ERROR(KK) LINFIT(KK)=0 10 CONTINUE C 30 FORMAT(/1X,79A1) 31 FORMAT(/' CURRENT CONSTANTS :'/) 33 FORMAT(1X,A6,F22.12,10X,I1) 34 FORMAT(/' MEASURED TRANSITIONS : ',I3/) 96 FORMAT(1X,I3,'. ',2I10,F19.4,F14.4) 97 FORMAT(1X,I3,'. ',2I10,F19.4,F14.4,'***INPUT ERROR') 98 FORMAT(1X,I3,'. ',2I10,F19.4,F14.4,' EXCL.') C C Selection and execution of required modification C 1900 WRITE(*,40) 40 FORMAT(/' MODIFICATIONS ? ',$) READ(*,14,ERR=1900)IFLAG 14 FORMAT(I5) IF(IFLAG.NE.1)GOTO 41 11 WRITE(*,70) 70 FORMAT(/' MODIFICATIONS ARE SELECTED AS FOLLOWS :' 1 //' = 0 LISTING AND EXIT TO CALCULATION',6X, 1 '= 5 INSERT LINES'/' = 1 CHANGE CONSTANTS',21X, 1 '= 6 CHANGE THE COMMENT'/' = 2 DELETE LINES',25X, 1 '= 7 STORE DATA'/' = 3 CHANGE LINES',25X, 1 '= 8 EXCLUDE LINES FROM FIT'/' = 4 ADD LINES',28X, 1 '= 9 INCLUDE LINES IN FIT'/43X,'= 10 CHANGE ERROR'/) 72 WRITE(*,71) 71 FORMAT(/' CONTROL PARAMETER = ',$) READ(*,14,ERR=72)ICONTR ICONTR=ICONTR+1 GOTO(29,74,75,76,77,500,178,179,210,220),ICONTR C C Change of constant C 74 WRITE(*,31) LL=NCONST/2 DO 37 L=1,LL WRITE(*,35)(K,T(K),A(K),IA(K),K=L,L+LL,LL) 37 CONTINUE 35 FORMAT(1X,I2,': ',A6,F26.16,I2,2X,I2,': ',A6,F26.16,I2) IF(2*LL.LT.NCONST)WRITE(*,36)NCONST,T(NCONST),A(NCONST), 1 IA(NCONST) 36 FORMAT(39X,2X,I2,': ',A6,F26.16,I2) 1910 WRITE(*,78) 78 FORMAT(/' NUMBER OF CONSTANT TO BE CHANGED ? ',$) READ(*,14,ERR=1910)NN IF(NN.LT.1.OR.NN.GT.NCONST)GOTO 72 1911 WRITE(*,79)T(NN) 79 FORMAT(/' FITTING PARAMETER, NEW VALUE OF ',A6,' ',$) READ(*,*,ERR=1911)IA(NN),A(NN) GOTO 74 C C Line deletion C 75 WRITE(*,81) 81 FORMAT(/' NUMBER OF LINE TO BE DELETED ? ',$) READ(*,14,ERR=75)NN IF(NN.LT.1.OR.NN.GT.I)GOTO 72 I=I-1 DO 1 L=NN,I FR(L,1)=FR(L+1,1) LINFIT(L)=LINFIT(L+1) ERROR(L)=ERROR(L+1) WEIGHT(L)=WEIGHT(L+1) DO 1 LL=1,NQNS JK(L,LL)=JK(L+1,LL) 1 CONTINUE GOTO 75 C C Change of lines C 76 WRITE(*,82) 82 FORMAT(/' NUMBER OF LINE TO BE CHANGED ? ',$) READ(*,14,ERR=76)NN IF(NN.LT.1.OR.NN.GT.I)GOTO 72 1912 WRITE(*,90)NN,(JK(NN,L),L=1,NQNS),FR(NN,1),ERROR(NN) 90 FORMAT(/' LINE NO.',I3,' :',2(I5,2I3),F18.4,F14.4/ 1 ' NEW LINE : ',$) READ(*,*,ERR=1912)(JK(NN,L),L=1,NQNS),FR(NN,1),ERROR(NN) IF(ERROR(NN).EQ.0.)ERROR(NN)=1. WEIGHT(NN)=1./ERROR(NN)**2 GOTO 76 C C Line addition C 77 I=I+1 WRITE(*,190) 190 FORMAT(/' LINE (LINES) TO BE ADDED :'/) 191 READ(*,*,ERR=1915)(JK(I,L),L=1,NQNS),FR(I,1),ERROR(I) GOTO 1916 1915 WRITE(*,'(1H+,A1)')CHAR(7) GOTO 191 1916 IF(JK(I,1).EQ.0.AND.JK(I,4).EQ.0)GOTO 192 LINFIT(I)=1 IF(ERROR(I).EQ.0.)ERROR(I)=1. WEIGHT(I)=1./ERROR(I)**2 I=I+1 GOTO 191 192 I=I-1 GOTO 72 C C Line insertion C 500 WRITE(*,501) 501 FORMAT(/' INSERTIONS ARE TO PRECEDE LINE NO. : ',$) READ(*,14,ERR=500)NINS WRITE(*,502) 502 FORMAT(/' LINE (LINES) TO BE INSERTED :'/) 503 READ(*,*,ERR=1918)(IINT(L),L=1,NQNS),FREQ,TEMP GOTO 1919 1918 WRITE(*,'(1H+,A1)')CHAR(7) GOTO 503 1919 IF(IINT(1).EQ.0.AND.IINT(4).EQ.0)GOTO 72 I=I+1 DO 504 L=I,NINS,-1 FR(L,1)=FR(L-1,1) LINFIT(L)=LINFIT(L-1) ERROR(L)=ERROR(L-1) WEIGHT(L)=WEIGHT(L-1) DO 504 LL=1,NQNS JK(L,LL)=JK(L-1,LL) 504 CONTINUE FR(NINS,1)=FREQ LINFIT(NINS)=1 IF(TEMP.EQ.0.)TEMP=1. ERROR(NINS)=TEMP WEIGHT(NINS)=1./TEMP**2 DO 505 LL=1,NQNS JK(NINS,LL)=IINT(LL) 505 CONTINUE NINS=NINS+1 GOTO 503 C C Change of comment C 178 WRITE(*,180) 180 FORMAT(/' NEW COMMENT :') READ(*,181,ERR=178)COM 181 FORMAT(79A1) GOTO 72 C C Saving of data C 179 IOUT=1 GOTO 200 C C Line exclusion C 210 WRITE(*,211) 211 FORMAT(/' NUMBER OF LINE TO BE EXCLUDED FROM FIT ? ',$) READ(*,14,ERR=210)NN IF(NN.LT.1.OR.NN.GT.I)GOTO 72 LINFIT(NN)=0 GOTO 210 C C Line inclusion C 220 WRITE(*,221) 221 FORMAT(/' NUMBER OF LINE TO BE INCLUDED IN FIT ? ',$) READ(*,14,ERR=220)NN IF(NN.LT.1.OR.NN.GT.I)GOTO 72 LINFIT(NN)=1 GOTO 220 C C Change of error C 300 WRITE(*,301) 301 FORMAT(/' ERROR TO BE CHANGED ON LINE NUMBER ? ',$) READ(*,*,ERR=300)NN IF(NN.LT.1.OR.NN.GT.I)GOTO 72 WRITE(*,302)NN,(JK(NN,K),K=1,NQNS),FR(NN,1),ERROR(NN) 302 FORMAT(/' LINE NO.',I3,' :',2(I5,2I3),F18.4,' +-',F14.4/ 1 ' NEW ERROR : ',$) READ(*,*,ERR=300)TEMP IF(TEMP.EQ.0.)TEMP=1. ERROR(NN)=TEMP WEIGHT(NN)=1./TEMP**2 GOTO 300 C C Sort of energy levels in NJA according to J, submatrix number and order C of eigenvalue in the submatrix C 41 DO 13 J=1,N-1 K1=J+1 L1=J DO 12 K=K1,N IF(NJA(K,1)-NJA(L1,1))411,401,12 401 IF(NJA(K,2)-NJA(L1,2))411,402,12 402 IF(NJA(K,3)-NJA(L1,3))411,12,12 411 L1=K 12 CONTINUE IF(L1.EQ.J)GOTO 13 DO 400 L=1,4 MNJA=NJA(J,L) NJA(J,L)=NJA(L1,L) NJA(L1,L)=MNJA 400 CONTINUE 13 CONTINUE N=N+1 NJA(N,1)=0 NJA(N,2)=0 NJA(N,3)=0 C 200 RETURN END C C_____________________________________________________________________________ C SUBROUTINE DATOUT(IFLAG,SDEV,SWDEV,NCONST,NQNS,idump) C C OUTPUT OF DATA AND RESULTS - BOTH TO SCREEN AND TO DISK C C IFLAG - on output signals whether further modifications are required C SDEV - standard deviation of the fit C NCONST - total no of fittable constants C NQNS - number of quantum numbers C also: C IOUT - if set to 1 on input, disk output of data is carried out only C if set to 2, disk output of both data and intermediate results C is enabled C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (MAXLIN=600,MAXDIM=140,MAXCON=35) COMMON /BLKH/H/BLKR/R COMMON /BLK1/JK,FR,ERROR,WEIGHT,LINFIT,A,IA,I, 1 IOUT/BLK3/C,D,DD/TBLOCK/COM,T INTEGER*2 JK(MAXLIN,6),IA(MAXCON),TPT(MAXCON),LINFIT(MAXLIN), * NCC(20),UNITS(MAXCON) REAL*8 FR(MAXLIN,3),A(MAXCON),C(MAXCON),H(MAXDIM,MAXDIM), * R(MAXDIM,MAXDIM),DD(MAXCON,MAXCON) REAL*4 D(MAXLIN,MAXCON),DER,DIV,ERROR(MAXLIN),WEIGHT(MAXLIN) CHARACTER*6 T(MAXCON) CHARACTER*25 FILNAM,CONVAL,ERVAL CHARACTER COM(79),FORMAR(54),symmet*2 CHARACTER*3 UNIT,NAMUN(4) CHARACTER*38 OUTCON(MAXCON) CHARACTER*54 FORMAL,F1,F2 EQUIVALENCE (FORMAL,FORMAR(1)) NOUT=0 C C...UNITS defines the units in which various constants are to be printed C out, NAMUN contains the names of these units DATA NAMUN/'MHz','kHz',' Hz','mHz'/ DATA UNITS/1,1, 2,2,2, 3,3,3,3, 1, 2,2, 3,3,3, 1,2,2, 3,3,3, 4, * 1,2,4,3, 4,4,4,4, 5*1/ C C...F1 - format for first line of printout of contributions, F2 - format C for all following lines for the same transition F1='(1X,I3,3H. ,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0)' F2='(1X,6H ,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0)' IF(IOUT.EQ.1)GOTO 100 C K1=0 DO 9 N=1,NCONST IF(IA(N).NE.1)GOTO 9 K1=K1+1 TPT(K1)=N 9 CONTINUE C IF(IOUT.EQ.2)GOTO 300 C C...Output of results (standard listing) which can be routed to screen C (NOUT=0) or to disk (NOUT=3) C 25 WRITE(NOUT,1)COM 1 FORMAT(1X,79(1H_)//1X,79A1/1X,79(1H_)//' Fit of K-doubling in v=1 1 E state of a C3v symmetric top '//) L=0 DO 2 N=1,NCONST IF(IA(N).EQ.0)GOTO 3 L=L+1 WRITE(NOUT,4)T(N),A(N),C(L) GOTO 2 3 WRITE(NOUT,5)T(N),A(N) 2 CONTINUE 4 FORMAT(1X,A6,F30.19,' +-',F25.19) 5 FORMAT(1X,A6,F30.19) WRITE(NOUT,405)SDEV,SWDEV 405 FORMAT(/' Deviation of fit = ',F14.6/ 1 ' Weighted deviation of fit = ',F14.6) c c...Write out lines with assigned symmetry: c c 1/ components with Kl-1 not divisible by 3 are E c 2/ A+ and A- denote upper and lower frequency components resp. c 3/ according to conventions of Cartwright and Mills (JMS 34,415-439,1970) c A components are such that, for positive q: c A+ = A1 A- = A2 for even J c A+ = A2 A- = A1 for odd J c WRITE(NOUT,6) 6 FORMAT(/10X,'J Kl-1 K l Observed ', * ' Obs-Calc Error'/) NINFIT=0 DO 7 N=1,I kl=jk(n,2)+1 kq=iabs(kl) lq=isign(1,kl) if(mod(iabs(jk(n,2)),3).eq.0)then symmet='A+' else symmet='E ' endif if(jk(n,1).lt.0)then kq=-kq lq=-lq symmet='A-' endif IF(LINFIT(N).NE.1)GOTO 70 WRITE(NOUT,8)N,(JK(N,L),L=1,NQNS),kq,lq,symmet,FR(N,1), * FR(N,3),ERROR(N) NINFIT=NINFIT+1 GOTO 7 70 WRITE(NOUT,71)N,(JK(N,L),L=1,NQNS),kq,lq,symmet,FR(N,1), * FR(N,3),ERROR(N) 7 CONTINUE c 8 FORMAT(1X,I3,'. ',I5,i4,3x,i4,',',i2,2x,a, * 2F15.4,f11.4) 71 FORMAT(1X,I3,'. ',I5,i4,3x,i4,',',i2,2x,a, * 2F15.4,f11.4,'EXCL') c IF(NOUT.NE.0)WRITE(NOUT,410) * SDEV,NINFIT,SWDEV,NINFIT-K1 410 FORMAT(/' Sigma = ',F14.6,28X,I4,' Transitions in fit'/ 1 ' Sigma.w = ',F14.6,28X,I4,' Degrees of freedom') C C...On output to disk additional conversion of constants into customary units C and writing in the convention of error given in units of last quoted C digit in constant value (three digit error) C IF(NOUT.EQ.3)THEN WRITE(NOUT,'(1X)') L=0 DO 1002 N=1,NCONST UNIT=NAMUN(UNITS(N)) ASCLD=A(N)*10.D0**(3*(UNITS(N)-1)) IF(IA(N).EQ.0)GOTO 1003 C C...value with error L=L+1 ESCLD=C(L)*10.D0**(3*(UNITS(N)-1)) WRITE(CONVAL,'(F25.16)')ASCLD WRITE(ERVAL,'(F25.16)')ESCLD CALL CONFOR(CONVAL,ERVAL,NDCON,NDEROR) WRITE(OUTCON(N),411)T(N)(1:5),UNIT,CONVAL(1:NDCON), * ERVAL(1:NDEROR) GOTO 1002 C C...value without error 1003 WRITE(CONVAL,'(F25.16)')ASCLD IF(ASCLD.EQ.0.0D0)CONVAL=' 0.0 ' IF(CONVAL(8:8).EQ.' ')CONVAL(8:8)='0' IF(CONVAL(8:8).EQ.'-')THEN CONVAL(8:8)='0' CONVAL(7:7)='-' ENDIF DO 1104 M=25,1,-1 IF(CONVAL(M:M).NE.'0'.AND.CONVAL(M:M).NE.' ')GOTO 1105 IF(CONVAL(M:M).EQ.'0')CONVAL(M:M)=' ' 1104 CONTINUE 1105 WRITE(OUTCON(N),1005)T(N)(1:5),UNIT,CONVAL 1002 CONTINUE LLLL=NCONST/2 DO 1107 N=1,LLLL WRITE(NOUT,1108)OUTCON(N),OUTCON(N+LLLL) 1107 CONTINUE IF(2*LLLL.NE.NCONST)WRITE(NOUT,1108)' ',OUTCON(NCONST) ENDIF 1108 FORMAT(1X,A38,1X,A38) 411 FORMAT(A,'/',A,A,'(',A,') ') 1005 FORMAT(A,'/',A,A) C IF(NOUT.EQ.0)GOTO 24 C C Variance covariance matrix (not printed) and correlation coefficients C (lower triangle only). Also mean values of correlation coefficients. C IF(K-1.LT.1)GOTO 31 SDEVSQ=SWDEV*SWDEV DO 200 N=1,K1 DO 200 L=N,K1 H(N,L)=DD(N,L)*SDEVSQ 200 CONTINUE DO 201 N=1,K1 DO 201 L=1,N H(N,L)=H(L,N)/(C(L)*C(N)) 201 CONTINUE C WRITE(NOUT,202) 202 FORMAT(/' CORRELATION COEFFICIENTS, C.ij:') DO 3035 N=1,20 NCC(N)=0 3035 CONTINUE NCORCO=0 CORCAV=0.D0 CCSUM=0.0D0 M=1 305 MM=M+7 IF(MM.GT.K1)MM=K1 WRITE(NOUT,301) WRITE(NOUT,230)(T(TPT(L)),L=M,MM) 230 FORMAT(7X,8(3X,A5,1X)) WRITE(NOUT,301) DO 302 N=M,K1 LLLL=MM IF(LLLL.GT.N)LLLL=N WRITE(NOUT,303)T(TPT(N)),(H(N,L),L=M,LLLL) DO 3030 L=M,LLLL IF(L.NE.N)THEN NCORCO=NCORCO+1 CORCAV=CORCAV+DABS(H(N,L)) CCSUM=CCSUM+H(N,L) NNN=DABS(H(N,L))/0.05+1 IF(NNN.EQ.21)NNN=20 NCC(NNN)=NCC(NNN)+1 IF(NCC(NNN).GT.67)NCC(NNN)=67 ENDIF 3030 CONTINUE 302 CONTINUE 303 FORMAT(1X,A5,1X,8(F9.4)) IF(MM.GE.K1)THEN GOTO 3031 ELSE M=M+8 GOTO 305 ENDIF 3031 IF(NCORCO.NE.0)THEN WRITE(NOUT,3032)CORCAV/NCORCO,CCSUM/NCORCO 3032 FORMAT(1X/' Mean value of |C.ij|, i.ne.j =', * F9.4/ * ' Mean value of C.ij, i.ne.j =', * F9.4/) ENDIF C C...bar graph of distribution of correlation coefficients IF(K.EQ.3)THEN WRITE(NOUT,'(1X/'' Distribution of values of |C.ij|, i.ne.j:'' * //10X,6(8X,I2)/10X,6(''....,....|''),''....,..'' * )')(N,N=10,60,10) DO 3036 N=1,20 WRITE(NOUT,3037)NCC(N),N*0.05,('*',LLLL=1,NCC(N)) 3036 CONTINUE 3037 FORMAT(1X,I2,' <',F4.2,1x,67A1) ENDIF C C Contributions to frequencies (disk only), F descriptor in output format C is chosen individually for each derivative according to its magnitude C C...values of contributions 304 IF(K-2.LT.1)GOTO 31 DO 203 N=1,I DO 203 L=1,K1 D(N,L)=D(N,L)*A(TPT(L)) 203 CONTINUE C C...heading WRITE(NOUT,204) 204 FORMAT(/' CONTRIBUTIONS OF INDIVIDUAL CONSTANTS TO FREQUENCY:'/) IST=1 ISP=8 390 IF(ISP.GT.K1)ISP=K1 WRITE(NOUT,230)(T(TPT(L)),L=IST,ISP) IF(ISP.LT.K1)THEN IST=IST+8 ISP=ISP+8 GOTO 390 ENDIF WRITE(NOUT,301) 301 FORMAT(1X) C DO 320 N=1,I C FORMAL=F1 IST=1 315 ISP=IST+7 IF(ISP.GT.K1)ISP=K1 314 DO 310 L=IST,ISP C C...determine LLL in format descriptor F9.LLL such that largest number C of significant digits are displayed for the current contribution LLL=5 DER=ABS(D(N,L)) DIV=10. 312 IF(DER/DIV.LT.1.)GOTO 311 DIV=DIV*10. LLL=LLL-1 GOTO 312 311 IF(LLL.LT.0)LLL=0 LLLL=L IF(LLLL.GT.8)LLLL=L-IST+1 FORMAR(13+LLLL*5)=CHAR(LLL+ICHAR('0')) 310 CONTINUE IF(ISP.LE.8)THEN WRITE(NOUT,FORMAL)N,(D(N,L),L=IST,ISP) ELSE WRITE(NOUT,FORMAL)(D(N,L),L=IST,ISP) ENDIF IF(ISP.EQ.K1)GOTO 320 FORMAL=F2 IST=IST+8 GOTO 315 C 320 CONTINUE C C 31 NOUT=0 CLOSE(3) GOTO 23 30 WRITE(NOUT,405)SDEV,SWDEV 24 WRITE(*,20) IFLAG=0 20 FORMAT(/' PRESS TO CONTINUE, 1 TO EXIT PROGRAM: ',$) READ(*,15,ERR=24)IIFLAG IF(IIFLAG.EQ.1)GOTO 500 IFLAG=1 IF(IOUT.NE.1.AND.IOUT.NE.2)GOTO 300 500 IOUT=0 GOTO 16 C 300 WRITE(*,22) 22 FORMAT(/' HARD COPY OUTPUT - CHOOSE ONE OF:'//20X,'0 = no hardcopy 1'/20X,'1 = standard output'/20X,'2 = so. + correlation coeff.'/ 1 20X,'3 = so. + cc. + distribution + contributions ... ',$) READ(*,15,ERR=300)K IF(K.LT.1.OR.K.GT.3)GOTO 23 NOUT=3 1901 WRITE(*,17) READ(*,'(A)',ERR=1901)FILNAM OPEN(3,FILE=FILNAM,STATUS='UNKNOWN',ERR=1901) GOTO 25 C 23 WRITE(*,14) 14 FORMAT(/' DATA TO BE WRITTEN ON DISK ? ',$) READ(*,15,ERR=23)L 15 FORMAT(I5) IOUT=1 IF(L.NE.1)GOTO 24 100 WRITE(*,17) 17 FORMAT(/' NAME TO BE USED FOR THE OUTPUT FILE : ',$) READ(*,'(A)',ERR=100)FILNAM OPEN(2,FILE=FILNAM,STATUS='UNKNOWN',ERR=100) WRITE(2,51)COM WRITE(2,52)I,idump WRITE(2,150)(T(K),IA(K),A(K),K=1,NCONST) WRITE(2,140)((JK(K,M),M=1,NQNS),FR(K,1),ERROR(K),LINFIT(K),K=1,I) CLOSE(2) 51 FORMAT(1X,79A1) 52 FORMAT(2I5) 140 FORMAT(2I5,F20.6,F15.6,I5) 150 FORMAT(A6,I5,F30.19) C IOUT=1 GOTO 24 C 16 RETURN END C C_____________________________________________________________________________ C SUBROUTINE CONFOR(CONVAL,ERVAL,NDCON,NDEROR) C C Constant and error formating for output C IMPLICIT INTEGER*2 (I-N) CHARACTER*25 CONVAL,ERVAL,CONOUT,EROUT C NERD=2 NDCON=0 NDEROR=0 NDNOTZ=0 ICZERO=0 C C...Go through digits of constant value adding those to output buffer while C checking at the same time digits of the error value, reacting as C necessary. Terminate when either three digits in the error value are C transferred or, if error has more than three digits before the decimal C point, the decimal point is reached. DO 1 N=1,25 NDCON=NDCON+1 CONOUT(NDCON:NDCON)=CONVAL(N:N) C C...ensure that zero precedes the decimal point and use ICZERO to monitor C whether decimal point has been reached IF(CONVAL(N:N).EQ.'.')ICZERO=1 IF(CONVAL(N:N).EQ.'.'.AND.CONVAL(N-1:N-1).EQ.' ') * CONOUT(NDCON-1:NDCON-1)='0' IF(CONVAL(N:N).EQ.'.'.AND.CONVAL(N-1:N-1).EQ.'-')THEN CONOUT(NDCON-2:NDCON-2)='-' CONOUT(NDCON-1:NDCON-1)='0' ENDIF C C...use NDNOTZ (number of digits not zero) to monitor whether significant C digits in constant value have been reached IF(CONVAL(N:N).GE.'1'.AND.CONVAL(N:N).LE.'9') * NDNOTZ=NDNOTZ+1 C C...do not transfer error digit if it is a leading zero, dot or space IF(NDEROR.EQ.0 .AND. (ERVAL(N:N).EQ.' '.OR. * ERVAL(N:N).EQ.'0'.OR.ERVAL(N:N).EQ.'.') )GOTO 1 C C...exit if error larger than value and decimal point reached IF(NDEROR.GE.NERD .AND. NDNOTZ.GT.0 .AND. ICZERO.EQ.1)GOTO 2 C C...do not transfer the dot in error value IF(ERVAL(N:N).EQ.'.')GOTO 1 C C...transfer error digits until three or, if more, the first significant C digit in value is reached NDEROR=NDEROR+1 EROUT(NDEROR:NDEROR)=ERVAL(N:N) IF(NDEROR.GE.NERD .AND. NDNOTZ.GT.0 .AND. ICZERO.EQ.1)GOTO 2 1 CONTINUE 2 CONVAL(1:NDCON)=CONOUT(1:NDCON) ERVAL(1:NDEROR)=EROUT(1:NDEROR) C RETURN END C C_____________________________________________________________________________ C SUBROUTINE HAMIL(I,J,N,A) C C Evaluation of a A or E v=1 symmetric top submatrix where: C C I - on output contains the dimensionality of the calculated subm. C J - J quantum number for the subm. C N - number of the required subm. = 1,2 for A,E resp. C A - array containing the spectroscopic constants C C Evaluated submatrix is placed in H C C Constant numbers: C C A B DJ DJK DK HJ HJK HKJ HK Azeta n.J n.K n.JJ n.JK n.KK Q.v Q.J C 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 C C qk qjj qjk qkk qjjj r.v q.3 h.6 t LJ LJJK LJK LKKJ C 18 19 20 21 22 23 24 25 26 27 28 29 30 C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (MAXDIM=140,MAXCON=35) COMMON /BLKH/H/BLKR/U DIMENSION A(MAXCON),H(MAXDIM,MAXDIM),U(MAXDIM,MAXDIM) C V=1 JJ1=(J+1)*J Q=JJ1 QQ=Q*Q C M=3*J/2+1 DO 2 I=1,M DO 2 K=1,M H(I,K)=0.D0 2 CONTINUE C C SET UP MATRICES. The A matrix contains only basis functions |K,l> such C that (Kl-1) is a multiple of three, E matrices contain all other C functions. (Kl-1) ranges from -(J+1) to J-1. C KA,KE are starting and terminal values for K, L is the C l quantum number (+1 or -1). C C Specimen matrix sizes J.....3..4..5..6..7.. ...3....4.....5.....6.....7.. C A 4 6 8 8 10 -2,1 -4,-1 -5,1 -5,1 -7,-1 C E 5 6 7 9 10 -3,-1 -4,1 -4,1 -6,-1 -7,1 C 2(2J+1) 14 18 22 26 30 C C...A-matrix: 1. Values for K omit 0 and multiples of three and range over C -(J-1) to (J-1) for J=3n C -J to J for J=3n+1 and J=3n+2 C 2. The starting value for l is chosen so that Kl-1 is a C multiple of three, l alternates in sign on iterating K C if(n.eq.1)then L=1 KE=J IF(J-3*(J/3).EQ.0)KE=KE-1 KA=-KE KLM1=KA*L-1 IF(KLM1-3*(KLM1/3).NE.0)L=-1 C i=0 DO 1 K=KA,KE IF(K-3*(K/3).EQ.0)GOTO 1 I=I+1 RK=K RK2=RK*RK RK4=RK2*RK2 RKL=K*L c H(I,I)=(A(1)-A(2))*RK2+A(2)*Q-2.D0*RKL*A(10) 1 -A(3)*QQ-A(4)*Q*RK2-A(5)*RK4 1 +A(6)*QQ*Q+A(7)*QQ*RK2+A(8)*Q*RK4+A(9)*RK4*RK2 1 +A(11)*Q*RKL+A(12)*RKL*RK2+A(13)*QQ*RKL 1 +A(14)*Q*RK2*RKL+A(15)*RK4*RKL 1 +A(27)*QQ*QQ+A(28)*QQ*Q*rk2+A(29)*QQ*rk2**2 1 +A(30)*Q*rk2**3 c if(k.eq.ke)goto 1 c c...(2,2) term IF(L.EQ.-1) 1 H(I,I+1)=0.25D0*(A(16)-A(17)*Q-A(18)*(rk+1)**2 1 +A(19)*QQ+A(20)*Q*(rk+1)**2+A(21)*(rk+1)**4 1 +A(22)*QQ*Q ) 1 *DSQRT((Q-RK*(RK+1.D0))*(Q-(RK+1.D0)*(RK+2.D0)) 1 *((V+1.D0)**2-(L+1)**2)) c c...(1,-2) term IF(L.EQ.1) 1 H(I,I+1)=A(23)*(2.D0*RK+1.D0) 1 *DSQRT((Q-RK*(RK+1.D0))*((V+1.D0)**2-(L-1)**2)) c c...(3,0) term IF(A(24).NE.0.0D0) 1 H(I,I+2)=A(24)*(2.D0*RK+3.D0)*DSQRT( 1 (Q- RK *(RK+1.D0)) * (Q-(RK+1.D0)*(RK+2.D0)) 1 * (Q-(RK+2.D0)*(RK+3.D0)) ) c c C...(6,0) term IF(A(25).NE.0.0D0)H(I,I+4)=A(25)*DSQRT( 1 (Q- RK *(RK+1.D0)) * (Q-(RK+1.D0)*(RK+2.D0)) 1 * (Q-(RK+2.D0)*(RK+3.D0)) * (Q-(RK+3.D0)*(RK+4.D0)) 1 * (Q-(RK+4.D0)*(RK+5.D0)) * (Q-(RK+5.D0)*(RK+6.D0)) ) c c...Kl-1=-3 splitting term IF(K.EQ.-2) 1 H(I,I+3)=A(26)*Q*(Q-2.D0) H(I+1,I)=H(I,I+1) H(I+2,I)=H(I,I+2) H(I+3,I)=H(I,I+3) H(I+4,I)=H(I,I+4) c L=-L 1 CONTINUE endif C C...E-matrix: 1. Values for K range over C -J to J for J=3n C -J to (J-1) for J=3n+1 C -(J-1) to J for J=3n+2 C and on iteration values which lead to multiple of 3 (Kl-1) C values are omitted C 2. The starting value of l is -1 for J=3n and 3n+2 and C 1 for J=3n+1, l alternates in sign on iterating K C if(n.eq.2)then KA=-J KE=J JJ=J-3*(J/3) IF(JJ.EQ.1)KE=J-1 IF(JJ.EQ.2)KA=KA+1 L=1 IF(JJ.EQ.0)L=-1 C i=0 DO 4 K=KA,KE KLM1=K*L-1 IF(KLM1-3*(KLM1/3).EQ.0)GOTO 4 I=I+1 RK=K RK2=RK*RK RK4=RK2*RK2 RKL=K*L c H(I,I)=(A(1)-A(2))*RK2+A(2)*Q-2.D0*RKL*A(10) 1 -A(3)*QQ-A(4)*Q*RK2-A(5)*RK4 1 +A(6)*QQ*Q+A(7)*QQ*RK2+A(8)*Q*RK4+A(9)*RK4*RK2 1 +A(11)*Q*RKL+A(12)*RKL*RK2+A(13)*QQ*RKL 1 +A(14)*Q*RK2*RKL+A(15)*RK4*RKL 1 +A(27)*QQ*QQ+A(28)*QQ*Q*rk2+A(29)*QQ*rk2**2 1 +A(30)*Q*rk2**3 c if(k.eq.ke)goto 4 c c...(2,2) term IF(L.EQ.-1) 1 H(I,I+1)=0.25D0*(A(16)-A(17)*Q-A(18)*(rk+1)**2 1 +A(19)*QQ+A(20)*Q*(rk+1)**2+A(21)*(rk+1)**4 1 +A(22)*QQ*Q ) 1 *DSQRT((Q-RK*(RK+1.D0))*(Q-(RK+1.D0)*(RK+2.D0)) 1 *((V+1.D0)**2-(L+1)**2)) c c...(1,-2) term IF(L.EQ.1) 1 H(I,I+1)=A(23)*(2.D0*RK+1.D0) 1 *DSQRT((Q-RK*(RK+1.D0))*((V+1.D0)**2-(L-1)**2)) c c...(3,0) term IF(A(24).NE.0.0D0)H(I,I+2)=A(24)*(2.D0*RK+3.D0)* 1 DSQRT( (Q-RK*(RK+1.D0))* (Q-(RK+1.D0)*(RK+2.D0)) 1 * (Q-(RK+2.D0)*(RK+3.D0)) ) c c C...(6,0) term IF(A(25).NE.0.D0)H(I,I+4)=A(25)*DSQRT( 1 (Q- RK *(RK+1.D0)) * (Q-(RK+1.D0)*(RK+2.D0)) 1 * (Q-(RK+2.D0)*(RK+3.D0)) * (Q-(RK+3.D0)*(RK+4.D0)) 1 * (Q-(RK+4.D0)*(RK+5.D0)) * (Q-(RK+5.D0)*(RK+6.D0)) ) c H(I+1,I)=H(I,I+1) H(I+2,I)=H(I,I+2) H(I+4,I)=H(I,I+4) c L=-L 4 CONTINUE endif C C--- write(*,'(1x/'' j,n,ka,ke,i: '',5i5)')j,n,ka,ke,i C--- IF(J.GE.39)WRITE(*,'(1X/I5)')I C RETURN END C C_____________________________________________________________________________ C SUBROUTINE HDIAG(N,IEGEN,NR) C C EIGENVALUES AND EIGENVECTORS OF A SQUARE, SYMMETRIC ARRAY C C N - dimensionality of array to be diagonalized C IEGEN - if =0 then calculate eigenvectors, if not then do not C NR - on output contains the number of iterations that were carried out C C Matrix to be diagonalized is placed in H C Eigenvalues are output on the diagonal of H, in ascending order C Eigenvectors are output in U in order of eigenvalues C C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (MAXDIM=140) COMMON /BLKH/H/BLKR/U DIMENSION H(MAXDIM,MAXDIM),U(MAXDIM,MAXDIM),X(MAXDIM), * IQ(MAXDIM) C IF(IEGEN.NE.0)GOTO 15 DO 14 I=1,N DO 14 J=1,N IF(I-J.NE.0)GOTO 12 U(I,J)=1.D0 GOTO 14 12 U(I,J)=0.D0 14 CONTINUE 15 NR=0 IF(N-1.LE.0)GOTO 1000 NMI1=N-1 DO 30 I=1,NMI1 X(I)=0.D0 IPL1=I+1 DO 30 J=IPL1,N IF(X(I)-ABS(H(I,J)).GT.0.D0)GOTO 30 X(I)=ABS(H(I,J)) IQ(I)=J 30 CONTINUE 40 DO 70 I=1,NMI1 IF(I-1.LE.0)GOTO 60 IF(XMAX-X(I))60,70,70 60 XMAX=X(I) IPIV=I JPIV=IQ(I) 70 CONTINUE EPS=1.0D-12 IF(XMAX-EPS)1000,1000,148 148 CONTINUE NR=NR+1 130 IF(H(IPIV,IPIV)-H(JPIV,JPIV))140,141,142 140 TANG=2.D0*H(IPIV,JPIV)/(H(IPIV,IPIV)-H(JPIV,JPIV) 1 -SQRT((H(IPIV,IPIV)-H(JPIV,JPIV))**2+4.D0*H(IPIV,JPIV)**2)) GOTO 150 141 TANG=1.D0 GOTO 150 142 TANG=2.D0*H(IPIV,JPIV)/(H(IPIV,IPIV)-H(JPIV,JPIV) 1 +SQRT((H(IPIV,IPIV)-H(JPIV,JPIV))**2+4.D0*H(IPIV,JPIV)**2)) 150 COSINE=1.D0/SQRT(1.D0+TANG**2) SINE=TANG*COSINE HII=H(IPIV,IPIV) H(IPIV,IPIV)=COSINE**2*(HII+TANG*(2.D0*H(IPIV,JPIV)+TANG* 1 H(JPIV,JPIV))) H(JPIV,JPIV)=COSINE**2*(H(JPIV,JPIV)-TANG*(2.D0*H(IPIV,JPIV)- 1 TANG*HII)) H(IPIV,JPIV)=0.D0 IF(H(IPIV,IPIV)-H(JPIV,JPIV))153,153,152 152 HTEMP=H(IPIV,IPIV) H(IPIV,IPIV)=H(JPIV,JPIV) H(JPIV,JPIV)=HTEMP HTEMP=SIGN(1.D0,-SINE)*COSINE COSINE=ABS(SINE) SINE=HTEMP 153 CONTINUE DO 350 I=1,NMI1 IF(I-IPIV)210,350,200 200 IF(I-JPIV)210,350,210 210 IF(IQ(I)-IPIV.EQ.0)GOTO 240 IF(IQ(I)-JPIV)350,240,350 240 K=IQ(I) 250 HTEMP=H(I,K) H(I,K)=0.D0 IPL1=I+1 X(I)=0.D0 DO 320 J=IPL1,N IF(X(I)-ABS(H(I,J)).GT.0.D0)GOTO 320 X(I)=ABS(H(I,J)) IQ(I)=J 320 CONTINUE H(I,K)=HTEMP 350 CONTINUE X(IPIV)=0.D0 X(JPIV)=0.D0 DO 530 I=1,N IF(I-IPIV)370,530,420 370 HTEMP=H(I,IPIV) H(I,IPIV)=COSINE*HTEMP+SINE*H(I,JPIV) IF(X(I)-ABS(H(I,IPIV)).GE.0.D0)GOTO 390 X(I)=ABS(H(I,IPIV)) IQ(I)=IPIV 390 H(I,JPIV)=-SINE*HTEMP+COSINE*H(I,JPIV) IF(X(I)-ABS(H(I,JPIV)).GE.0.D0)GOTO 530 X(I)=ABS(H(I,JPIV)) IQ(I)=JPIV GOTO 530 420 IF(I-JPIV)430,530,480 430 HTEMP=H(IPIV,I) H(IPIV,I)=COSINE*HTEMP+SINE*H(I,JPIV) IF(X(IPIV)-ABS(H(IPIV,I)).GE.0.D0)GOTO 450 X(IPIV)=ABS(H(IPIV,I)) IQ(IPIV)=I 450 H(I,JPIV)=-SINE*HTEMP+COSINE*H(I,JPIV) IF(X(IPIV)-ABS(H(I,JPIV)).GE.0.D0)GOTO 530 X(I)=ABS(H(I,JPIV)) IQ(I)=JPIV GOTO 530 480 HTEMP=H(IPIV,I) H(IPIV,I)=COSINE*HTEMP+SINE*H(JPIV,I) IF(X(IPIV)-ABS(H(IPIV,I)).GE.0.D0)GOTO 500 X(IPIV)=ABS(H(IPIV,I)) IQ(IPIV)=I 500 H(JPIV,I)=-SINE*HTEMP+COSINE*H(JPIV,I) IF(X(JPIV)-ABS(H(JPIV,I)).GE.0.D0)GOTO 530 510 X(JPIV)=ABS(H(JPIV,I)) IQ(JPIV)=I 530 CONTINUE IF(IEGEN.NE.0)GOTO 40 DO 550 I=1,N HTEMP=U(I,IPIV) U(I,IPIV)=COSINE*HTEMP+SINE*U(I,JPIV) U(I,JPIV)=-SINE*HTEMP+COSINE*U(I,JPIV) 550 CONTINUE GOTO 40 1000 CONTINUE C RETURN END C_____________________________________________________________________________ C SUBROUTINE JACK(N) c c JACOBI type matrix diagonalization routine for real, symmetric matrices c (optimised version of AUTOValori from D.Lister's LINSHP program) c c H - array to be diagonalised (contains the eigenvalues on the diagonal c on output c S - contains the eigenvectors on output (in order of eigenvalues) c N - dimension of the array to be diagonalised contained in the top c left hand corner of H c c IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (MAXDIM=140) COMMON /BLKH/H/BLKR/S REAL*8 H(MAXDIM,MAXDIM),S(MAXDIM,MAXDIM) c c ... initialization and setup of limits c ss=0.7071067811865475d0 n1=n-1 sum=0.d0 do 20 i=1,n1 s(i,i)=1.d0 k=i+1 do 20 m=k,n sum=sum+h(i,m)**2 s(i,m)=0.d0 s(m,i)=0.d0 20 continue s(n,n)=1.d0 unu=dsqrt(2.d0*sum) ene=n c c ... convergence criterion c unufi=1.d-12 c c c ... iteration sequence - determination of SIN and COS for rotation c 54 unu=unu/ene 55 ind=0 do 89 m=2,n l=m-1 do 89 i=1,l if(h(i,m))60,89,60 60 q=dabs(h(i,m)) if(q-unu)89,62,62 62 ind=1 if(h(i,i)-h(m,m))68,64,68 64 si=dsign(ss,h(i,m)) co=ss goto 77 68 x1=-h(i,m)*2.d0/(h(i,i)-h(m,m)) x1a=dabs(x1) if(1.d30-x1a)66,66,71 66 si=dsign(ss,x1) co=ss goto 77 71 if(1.d-30-x1a)72,89,89 72 t=datan(x1)*0.5d0 si=dsin(t) co=dcos(t) c c ... transformation loop c 77 do 79 j=1,n hji=h(j,i) sji=s(j,i) h(j,i)=hji*co-h(j,m)*si h(j,m)=hji*si+h(j,m)*co s(j,i)=sji*co-s(j,m)*si s(j,m)=sji*si+s(j,m)*co if(j.eq.i.or.j.eq.m)goto 79 h(i,j)=h(j,i) h(m,j)=h(j,m) 79 continue c h(i,i)=co*h(i,i)-si*h(m,i) h(m,m)=co*h(m,m)+si*h(i,m) h(i,m)=0.d0 h(m,i)=0.d0 89 continue if(ind)92,91,55 91 if(unu-unufi)92,92,54 c 92 return end C C_____________________________________________________________________________ C_____________________________________________________________________________