c$debug MS-Fort C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C Q2FIT - PROGRAM TO FIT THE SPECTRUM OF AN ASYMMETRIC ROTOR WITH C ONE OR TWO QUADRUPOLAR NUCLEI AND SPIN-SPIN AND SPIN-ROTATION C INTERACTIONS C c This program was originally developed to deal with strong c accidental quadrupole-mediated perturbation observed in c the hydrogen bonded dimer vinyl fluoride...HCl C C Z.Kisiel,P.W.Fowler,A.C.Legon - J.Chem.Phys. 93,3054-3062(1990) C C C Ver 8.II.2001 ------- Zbigniew KISIEL ------- C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, 02-668 Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C | http://info.ifpan.edu.pl/~kisiel/prospe.htm | C_________________________/-------------------------------------------------- c C Modification history: C c ca. 1990: creation from program READ c ca. 1994: addition of spin-rotation corrections from Claudio di Esposti c 9.II.01: elimination of longstanding bug in LSQ which affected deviations c on constants when there was a large number of unfitted lines C C_____________________________________________________________________________ C C COUPLING SCHEME: I = I1 + I2 C F = J + I C C FOR MATRIX ELEMENTS SEE: C 1. Keenan,Wozniak,Flygare - J.Chem.Phys. 75,631(1981) C 2. Read,Flygare - J.Chem.Phys. 76,2238(1982) C C - MAXIMUM NUMBER OF TRANSITIONS IS DETERMINED BY VARIABLE 'NLIN' IN C PARAMETER STATEMENTS C - QUANTUM NUMBER RANGE IS LIMITED BY THE DOUBLE SUM OF (2J+1) C OVER I AND J FOR A GIVEN F BEING EQUAL TO 'MDIM', DECLARED IN C RELEVANT PARAMETER STATEMENTS. SETTING MDIM=168 AND NLIN=200 C ALLOWS THE PROGRAM TO RUN IN THE LOW MEMORY IN A PC (see comments C in routine MHIGH on how to determine the required value of MDIM) C - ALL SEXTIC CENTRIFUGAL DISTORTION CONSTANTS AND LKKJ C - THE FIT IS NOT WEIGHTED BUT LINES CAN BE EXCLUDED FROM THE FIT C (for predictive purposes, for example) C C C This program is a fitting version of modified program READ written C by G.W.Read ca 1981, containing angular momentum routines written C by M.P.Keenan ca 1980. C C NOTE: Q2FIT has been exhaustively tested for the single nucleus C case including X.ab. The calculation for two nuclei has been tested C with Ar...ClCN. Spin rotation has only been tested in the C symmetric limit, and spin-spin has not yet been tested at all. C Apart from the -q error in 3j symbol eq.11 of paper [1], the last two C spherical tensor components in eq.12 should have signs -+ and + resp. C (instead of +- and +-). The off-diagonal quadrupole constants are C actually X.ab+-iX.ac and X.bb-X.cc+-2iX.bc, the use of a single C constant X.ab here, in real formulation means that only X.ab or X.ac C can be fitted under conditions of the other two being zero (a version C of program for all three has also been written) C In paper 2 eq.A3 the closing square bracket after delta(J'J)delta(K'K) C is missing. C C COMPILATION NOTES: C C Compile with any 32-bit compiler, remembering to use the appropriate c option for static allocation of variables (eg. -static with f77) C_____________________________________________________________________________ C C C DATA FILE IS TO HAVE THE FOLLOWING FORM: C ======================================== C C Line.1: COMENT C Line.2: I1 I2 C Line.3: NIT IDUMP ICONTR C Line.4: NT C C LINE.5: DUMMY IA(1) CON(1) :FORMAT (A14,I2,F30), DUMMY is a C . . :14 character dummy alphan. varia- C Line.28: DUMMY IA(31) CON(31) :ble allowing for use of constant C :names in the data file C C Line.29: J' K-1' K+1' J K-1 K+1 I' F' I F Frequency LINFIT(1) C . . C Line.28+NT: J' K-1' K+1' J K-1 K+1 I' F' I F Frequency LINFIT(NT) C C C The meaning of the various parameters is explained in the C section below. C All lines, with the exception of the constant specification C lines are read in free format. Primed quantum numbers refer C to upper level in the transition. C C This version of Q2FIT is designed for semi-interactive running C - the user is asked to provide the input and output file names C from the console and can follow the progress of the fit by C means of intermediate results displayed on the screen. C C_____________________________________________________________________________ C C COMENT - up to 78 character alphanumeric comment C I1,I2 - spins of the two quadrupolar nuclei (in units of 1/2) C FH,FL - highest and lowest values of the quantum number F (in units C of 1/2) C C IA(i) - fitting parameters for the spectroscopic constansts C (=1 constant is to be fitted, =0 unchanged) C CON(i) - spectroscopic constants C NCON - number of constants to be fitted C ICFIT(i) - pointers to constants to be fitted (ie. only the first NCON C elements are set C PO(i) - names of constants used in the printout C POSHRT(i) - short names of constants for correlation coefficients and C contributions parts of printout C X(i) - changes in values of constants after a least squares iteration C SD(i) - standard deviations on fitted constants C C NT - number of transitions C LABO(i,j,k) - quantum numbers for transitions: C i is the number of the transition C j defines the level involved in the transition and is 1 for the C lower and 2 for the upper C k is a quantum number index for a level in the transition: C =1 I C =2 J C =3 Tau+J+1 which runs from 1 to 2J+1 for a given J, and counts C the postion of the energy level in sorted levels in a given C J block C =4 F C FREQ(i,j) - frequencies of transitions for i from 1 to NT where: C (i,1) observed frequencies C (i,2) calculated frequencies C (I,3) obs-calc C LINFIT(i) - pointer to whether line i is to be fitted (=1) or omitted C from the fit, ie predicted only (=0) C INFOE(i,j) - information on energy levels to be calculated where: C (i,1) points to the transition the level is involved in and C is negative for the lower, and positive for the upper level C in the transition C (i,2) is the value of F for the level (in units of 1/2) C C NIT - number of iterations C IDUMP - a dump of the eigenmatrix of H and the resulting eigenvalues C is produced when this is set to 2, eigenvalues only when 1 and C no dump when 0. The dump is to a file called DUMP. C ICONTR - contributions of fitted constants to frequencies are produced C when this is set to 1, quadrupole removed frequencies are C calculated when set to 2 (otherwise set to 0) C C_____________________________________________________________________________ C C IMPLICIT REAL*8(A-E,G-H,O-Z),INTEGER(F) PARAMETER (NCONST=31,NLIN=200,NLIN2=2*NLIN) COMMON /BLKEN/TFREQ,DER,LABO,INFOE,NELEV,ICFIT COMMON /BLKFIT/LINFIT,CORELC(NCONST,NCONST) DIMENSION CON(NCONST),IA(NCONST) DIMENSION LABO(NLIN,2,4),DER(NLIN,NCONST),LINFIT(NLIN) DIMENSION ICFIT(NCONST),TFREQ(NLIN,3),INFOE(NLIN2,2) DIMENSION SD(NCONST),X(NCONST),INQ(NLIN,10),UNITS(NCONST) CHARACTER*78 COMENT CHARACTER*30 FILNAM CHARACTER*14 DUMMY CHARACTER*12 PO(NCONST) CHARACTER*6 POSHRT(NCONST) CHARACTER FORMAR(54),CONVAL*25,ERVAL*25 CHARACTER*39 OUTCON(NCONST) CHARACTER*7 UNIT,NAMUN(6) CHARACTER*54 FORMAL,F1,F2 EQUIVALENCE (FORMAL,FORMAR(1)) 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 ','microHz', * 'nHz '/ DATA UNITS/1,1,1, 2,2,2,2,2, 3,3,3,3,3,3,3, 4, 1,1,1, * 1,1,1, 1,1,1, 1,1,1, 2,2,2/ C DATA PO/' A ',' B ',' C ', 1 ' DJ ',' DJK ',' DK ', 1 ' dJ ',' dK ',' HJ ', 1 ' HJK ',' HKJ ',' HK ', 1 ' hJ ',' hJK ',' hK ', 1 ' LKKJ ', 1 'I1: X.aa ','I1: Xbb-Xcc ','I1: Xab ', 1 'I2: X.aa ','I2: Xbb-Xcc ','I2: Xab ', 1 'I1: Maa ','I1: Mbb-Mcc ','I1: Trace(M)', 1 'I2: Maa ','I2: Mbb-Mcc ','I2: Trace(M)', 1 ' Daa ',' Dbb-Dcc ',' Dba '/ DATA POSHRT/'A ','B ','C ','DJ ','DJK ', 1 'DK ','dJ ','dK ','HJ ', 1 'HJK ','HKJ ','HK ','hJ ', 1 'hJK ','hK ','LKKJ ', 1 '1:Xa ','1:Xb-c','1:Xab ', 1 '2:Xa ','2:Xb-c','2:Xab ','1:Ma ','1:Mb-c','1:Tr ', 1 '2:Ma ','2:Mb-c','2:Tr ','Da ','Db-c ','Dba '/ 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)' C C C...Input C WRITE(*,200) 200 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | Q2FIT - FIT OF SPECTRUM OF ASYMMETRIC TOP WITH TWO QUAD', * 'RUPOLAR NUCLEI',T79,'|'/ * ' | (cd. in reduction A + ss. + sr.)', * T79,'|'/ * ' |',76(1H_),'|'/' version 8.II.2001',T64,'Zbigniew KISIEL'/) 202 WRITE(*,201)'Name of data file: ' 201 FORMAT(1X/1X,A,$) READ(*,'(A)',ERR=202)FILNAM OPEN(2,FILE=FILNAM,status='old',ERR=202) READ(2,'(A)')COMENT READ(2,*)I1,I2 READ(2,*)NIT,IDUMP,ICONTR READ(2,*) NT C NCON=0 DO 350 I=1,NCONST READ(2,370)DUMMY,IA(I),CON(I) 370 FORMAT(A14,I2,F30.0) IF(IA(I).EQ.0)GOTO 350 NCON=NCON+1 ICFIT(NCON)=I 350 CONTINUE C NLFIT=0 DO 2 I=1,NT READ(2,*,ERR=351)(INQ(I,J),J=1,10),TFREQ(I,1),LINFIT(I) IF(LINFIT(I).EQ.1)NLFIT=NLFIT+1 LABO(I,1,1)=INQ(I,9) LABO(I,1,2)=INQ(I,4) LABO(I,1,3)=INQ(I,5)-INQ(I,6)+INQ(I,4)+1 LABO(I,1,4)=INQ(I,10) LABO(I,2,1)=INQ(I,7) LABO(I,2,2)=INQ(I,1) LABO(I,2,3)=INQ(I,2)-INQ(I,3)+INQ(I,1)+1 LABO(I,2,4)=INQ(I,8) 2 CONTINUE GOTO 352 351 WRITE(*,'(1X/1X,A,I5,a/10x,10i4,f15.3/)') * '**** INPUT ERROR ON LINE =',I,' following:',(inq(i-1,j), * j=1,10),tfreq(i-1,1) STOP 352 CLOSE(2) C WRITE(*,257)NT,NLFIT,NCON 257 FORMAT(1X/1X,I3,' Lines read in'/1X,I3,' Lines are to be used ', 1 'in the fit'/1X,I3,' Constants are to be fitted') NDEGF=NLFIT-NCON IF(NDEGF.LT.1)NDEGF=1 C C...Initialise output C 203 WRITE(*,201)'Name of output file: ' READ(*,'(A)',ERR=203)FILNAM OPEN(3,FILE=FILNAM,ERR=203,STATUS='UNKNOWN') WRITE(3,1)COMENT 1 FORMAT(1X,78(1H-)/1X,A/1X,78(1H-)/) C C...Extract the energy levels to be calculated, check the validity of quantum C numbers and then sort the energy levels in ascending order of F C NELEV=0 DO 230 I=1,NT DO 231 J=1,2 NELEV=NELEV+1 INFOE(NELEV,1)=-I IF(J.EQ.2)INFOE(NELEV,1)=I INFOE(NELEV,2)=LABO(I,J,4) 231 CONTINUE 230 CONTINUE C DO 232 J=1,NELEV-1 K1=J+1 L1=J DO 233 K=K1,NELEV IF(INFOE(K,2).GE.INFOE(L1,2))GOTO 233 L1=K 233 CONTINUE IF(L1.EQ.J)GOTO 232 DO 234 L=1,2 M=INFOE(J,L) INFOE(J,L)=INFOE(L1,L) INFOE(L1,L)=M 234 CONTINUE 232 CONTINUE FH=INFOE(NELEV,2) c if(idump.eq.1.or.idump.eq.2)open(4,file='dump',status='unknown') dump C C C------------------ MAIN FITTING LOOP -------------------- C C N=0 DO 15 LP=1,NIT C if(idump.ne.1.and.idump.ne.2)goto 1302 dump do 1300 i=1,NCONST dump if(con(i).eq.0.d0)goto 1300 dump write(4,1301)po(i),con(i) dump 1301 format (1X,A,' =',F20.9) dump 1300 continue dump c 1302 DO 3 J=1,NT TFREQ(J,2)=0.D0 DO 3 M=1,NCON DER(J,M)=0.0D0 3 CONTINUE C WRITE(*,210)'CALL TO READ - Iteration',LP 210 FORMAT(1X/1X,A,I3) C CALL READ (CON,NCON,I1,I2,FH,IDUMP) C WRITE(3,255) 255 FORMAT(1X/' TRANSITIONS:'// 1 7X,'J K K <- J K K I F <- I F',7X, 2 'Obs',13X,'Obs-Calc'/ 3 10X, '-1 +1 -1 +1'/) SUMSQ=0.D0 DO 7 I=1,NT TFREQ(I,3)=TFREQ(I,1)-TFREQ(I,2) IF(LINFIT(I).NE.1)GOTO 251 SUMSQ=SUMSQ+TFREQ(I,3)**2 WRITE(3,211)I,(INQ(I,J),J=1,10),TFREQ(I,1),TFREQ(I,3) 211 FORMAT(1X,I3,'. ',3I3,' ',3I3,2X,2(I3,'/2'),' ',2(I3,'/2'), 1 F15.5,F14.5) GOTO 7 251 WRITE(3,256)I,(INQ(I,J),J=1,10),TFREQ(I,1),TFREQ(I,3) 256 FORMAT(1X,I3,'. ',3I3,' ',3I3,2X,2(I3,'/2'),' ',2(I3,'/2'), 1 F15.5,F14.5,'--') 7 CONTINUE SDEVSQ=DSQRT(SUMSQ/NDEGF) WRITE(3,212)SDEVSQ WRITE(*,212)SDEVSQ 212 FORMAT(1X/' Standard deviation = ',F14.6/) C C...Determine increments to constants by linear least squares C CALL LSQR(NT,NCON,SD,X) DO 10 I=1,NCON J=ICFIT(I) CON(J)=CON(J)+X(I) 10 CONTINUE C C...Write constants predicted by each iteration C WRITE(3,11)LP WRITE(*,11)LP 11 FORMAT (1X/' ITERATION NO =',I2,' CONSTANTS, deviations and', 1 ' changes:'/) DO 13 I=1,NCON J=ICFIT(I) WRITE(*,12)PO(J),CON(J),SD(I),X(I) WRITE(3,12)PO(J),CON(J),SD(I),X(I) 12 FORMAT (1X,A,' =',F20.9,' +-',F17.9,F22.9) 13 CONTINUE 14 CONTINUE 15 CONTINUE c C--------------------------------------------------------- C C Write final results C C...Constants, with additional conversion of values 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 WRITE(3,16) 16 FORMAT (1X///' FINAL RESULTS OF LEAST SQUARES FITTING PROCEDURE' 1 /1X,48(1H=)//' FITTED CONSTANTS:'/) C L=0 DO 1002 N=1,NCONST UNIT=NAMUN(UNITS(N)) ASCLD=CON(N)*10.D0**(3*(UNITS(N)-1)) IF(IA(N).EQ.0)GOTO 1003 C C...value with error L=L+1 ESCLD=SD(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)POSHRT(N),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)POSHRT(N),UNIT,CONVAL 1002 CONTINUE LLLL=NCONST/2 DO 1107 N=1,LLLL WRITE(3,1108)OUTCON(N),OUTCON(N+LLLL) 1107 CONTINUE IF(2*LLLL.NE.NCONST)WRITE(3,1108)' ',OUTCON(NCONST) 1108 FORMAT(1X,A38,1X,A38) 411 FORMAT(A,'/',A,A,'(',A,') ') 1005 FORMAT(A,'/',A,A) C write(3,'(1x)') IF(CON(17).NE.0.D0.OR.CON(18).NE.0.0D0)THEN XCC=-0.5D0*(CON(17)+CON(18)) XBB=CON(18)+XCC WRITE(3,115)'1:',XBB,'1:',XCC ENDIF IF(CON(20).NE.0.D0.OR.CON(21).NE.0.0D0)THEN XCC=-0.5D0*(CON(20)+CON(21)) XBB=CON(21)+XCC WRITE(3,115)'2:',XBB,'2:',XCC ENDIF 115 FORMAT(1X,A,'X.bb /MHz ',F23.12 1 /1X,A,'X.cc /MHz ',F23.12) C C...Correlation coefficients C WRITE(3,121) 121 FORMAT(/' CORRELATION COEFFICIENTS:') NMIN=1 NMAX=8 124 IF(NMAX.GT.NCON)NMAX=NCON WRITE(3,'(1X)') WRITE(3,120)(POSHRT(ICFIT(I)),I=NMIN,NMAX) 120 FORMAT(7X,8(3X,A6)) WRITE(3,'(1X)') DO 122 I=1,NCON LMAX=NMAX IF(I.LT.LMAX)LMAX=I IF(LMAX.LT.NMIN)GOTO 122 WRITE(3,123)POSHRT(ICFIT(I)),(CORELC(I,J),J=NMIN,LMAX) 122 CONTINUE 123 FORMAT(1X,A6,8(F9.4)) IF(NMAX.EQ.NCON)GOTO 125 NMIN=NMIN+8 NMAX=NMAX+8 GOTO 124 C C...Contributions to frequencies (optional) C 125 IF(ICONTR.LT.1)GOTO 126 DO 130 J=1,NCON K=ICFIT(J) DO 130 I=1,NT DER(I,J)=DER(I,J)*CON(K) 130 CONTINUE C IF(ICONTR.EQ.2)THEN WRITE(3,153) 153 FORMAT(1X/' TRANSITIONS:'// 1 7X,'J K K <- J K K I F <- I F',7X, 2 'Obs',10X,'Quadr. free'/ 3 10X, '-1 +1 -1 +1'/) DO 151 I=1,NT TFREQ(I,3)=TFREQ(I,1) DO 152 J=1,NCON K=ICFIT(J) IF(K.GE.17)THEN TFREQ(I,3)=TFREQ(I,3)-DER(I,J) if(k.eq.19.or.k.eq.22) * tfreq(i,3)=tfreq(i,3)+0.5d0*der(i,j) ENDIF 152 CONTINUE WRITE(3,211)I,(INQ(I,J),J=1,10),TFREQ(I,1),TFREQ(I,3) 151 CONTINUE GOTO 126 ENDIF C WRITE(3,131) 131 FORMAT(//' CONTRIBUTIONS OF INDIVIDUAL CONSTANTS TO FREQUENCY:'/) NMIN=1 NMAX=8 132 IF(NMAX.GT.NCON)NMAX=NCON WRITE(3,120)(POSHRT(ICFIT(I)),I=NMIN,NMAX) IF(NMAX.EQ.NCON)GOTO 133 NMIN=NMIN+8 NMAX=NMAX+8 GOTO 132 133 WRITE(3,'(1X)') C DO 134 N=1,NT FORMAL=F1 NMIN=1 NMAX=8 138 IF(NMAX.GT.NCON)NMAX=NCON DO 135 L=NMIN,NMAX LLL=5 DERIV=DABS(DER(N,L)) DIV=10.D0 137 IF(DERIV/DIV.LT.1.D0)GOTO 136 DIV=DIV*10.D0 LLL=LLL-1 GOTO 137 136 IF(LLL.LT.0)LLL=0 LLLL=L IF(LLLL.GT.8)LLLL=LLLL-8 FORMAR(13+LLLL*5)=CHAR(LLL+ICHAR('0')) 135 CONTINUE IF(NMIN.GT.1)GOTO 139 WRITE(3,FORMAL)N,(DER(N,L),L=NMIN,NMAX) GOTO 140 139 WRITE(3,FORMAL)(DER(N,L),L=NMIN,NMAX) 140 IF(NMAX.EQ.NCON)GOTO 134 FORMAL=F2 NMIN=NMIN+8 NMAX=NMAX+8 GOTO 138 134 CONTINUE C C 126 CLOSE(3) C STOP END C C_____________________________________________________________________________ C SUBROUTINE CONFOR(CONVAL,ERVAL,NDCON,NDEROR) C C Constant and error formating for output C CHARACTER*25 CONVAL,ERVAL,CONOUT,EROUT C NERD=3 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 READ (CON,NCON,I1,I2,FH,IDUMP) C C Routine to: C C - set up the Hamiltonian for given quantum numbers I1,I2 and F, C and spectroscopic constants in CON C - find the energy levels required for calculated frequencies C - evaluate the derivatives of the energy levels, and hence the C transitions with respect to the NCON constants to be fitted C C NELINF(i,j) - information on eigenvalues of H required by the spectral data: C (i,1) is the number of transition the level belongs to (-ve for C the lower level) C (i,2) is the number of this eigenvalue among the eigenvalues of H C NEXTR - number of elements of NELINF that set for each value of F C C NCOL - the number of columns of H that the printer being used can C print on one line C 17 is suitable for a 15" printer in condensed mode C 12 is suitable for 10" in condensed elite C IMPLICIT REAL*8(A-E,G-H,O-Z),INTEGER(F) c PARAMETER (MDIM=1600,NCONST=31,NLIN=200,NLIN2=2*NLIN) PARAMETER (MDIM=300,NCONST=31,NLIN=200,NLIN2=2*NLIN) COMMON /BLKEN/TFREQ,DER,LABO,INFOE,NELEV,ICFIT COMMON /HBLK/H1/SBLK/T COMMON /SORTCC/WK(MDIM),IPT INTEGER*2 IPT(MDIM) DIMENSION H1(MDIM,MDIM),NELINF(MDIM,2),LABO(NLIN,2,4), * ICFIT(NCONST),CON(NCONST),CONST(NCONST),TFREQ(NLIN,3), * INFOE(NLIN2,2),DER(NLIN,NCONST) REAL*4 T(MDIM,MDIM) c NCOLS=12 C C C...Set up and diagonalize the Hamiltonian for each required value of F, in C ascending order from the lowest in the data set to the highest. C The necessary eigenvalues and their derivatives are extracted C following the diagonalisation for each F. C FMIN=INFOE(1,2) NEXTF=FMIN C DO 100 F=FMIN,FH,2 IF(F.NE.NEXTF)GOTO 100 IF(FMIN.EQ.0)FM=(F+2)/2 IF(FMIN.EQ.1)FM=(F+1)/2 MS=MHIGH(I1,I2,F) WRITE(*,220)F,MS,MS 220 FORMAT(10X,'F = ',I3,'/2',5X,'H size = ',I3,' x',I3) IF(MS.GT.MDIM)GOTO 103 C C...set up and diagonalise H C DO 105 I=1,NCONST CONST(I)=CON(I) 105 CONTINUE CALL SETHAM(MS,I1,I2,F,CONST) IF(IDUMP.NE.1.AND.IDUMP.NE.2)GOTO 1014 c write(4,'(1x)') dump if(idump.eq.1)goto 1014 dump nmin=1 dump nmax=ms dump if(ms.gt.ncols)nmax=ncols dump 1011 do 1001 i=1,ms dump write(4,1002)i,(h1(i,j),j=nmin,nmax) dump 1001 continue dump 1002 format(1x,i2,':',60f13.4) dump if(nmax.eq.ms)goto 1010 dump nmin=nmin+ncols dump nmax=nmax+ncols dump if(nmax.gt.ms)nmax=ms dump write(4,'(1x)') dump goto 1011 dump 1010 continue dump c 1014 IEGEN=0 IORD=0 CALL HDIAG(MS,MDIM,IEGEN,IORD) write(*,'('' finished '',i3,''x'',i3)')ms,ms C DO 1020 I=1,MS WK(I)=H1(I,I) IPT(I)=I 1020 CONTINUE c mm=0 ih=I1+I2 il=iabs(I1-I2) do 1003 i=il,ih,2 jh=F+i jl=iabs(i-F) DO 1003 j=jl,jh,2 IF(IDUMP.NE.0) dump * write(4,1005)i,F,j/2,(h1(nn,nn),nn=mm+1,mm+(j+1)) dump CALL SORTC(mm+1,mm+j+1) IF(IDUMP.NE.0) dump * write(4,1025)i,F,j/2,(h1(IPT(nn),IPT(nn)), dump * nn=mm+1,mm+j+1) dump mm=(j+1)+mm 1003 continue 1005 format(1X,'I=',I2,'/2 F=',I2,'/2 J=',I2,' RAW E-levels:', dump * 20(/1X,11F13.4)) dump 1025 format(1X,'I=',I2,'/2 F=',I2,'/2 J=',I2,' SORTED LEVELS:', dump * 20(/1X,11F13.4)) dump C C C...determine which eigenvalues will be necessary and load the relevant C eigenvalues into TFREQ(i,2) with the appropriate sign C 1015 NEXTR=0 DO 102 LL=1,NELEV IF(INFOE(LL,2).LT.F)GOTO 102 IF(INFOE(LL,2).GT.F)GOTO 106 NLEV=2 IF(INFOE(LL,1).LT.0)NLEV=1 NLINE=IABS(INFOE(LL,1)) I=LABO(NLINE,NLEV,1) J=LABO(NLINE,NLEV,2) K=LABO(NLINE,NLEV,3) c CALL NUMENL(I1,I2,F,I,J,K,M) M=IPT(M) c 120 NEXTR=NEXTR+1 NELINF(NEXTR,1)=INFOE(LL,1) NELINF(NEXTR,2)=M TFREQ(NLINE,2)=TFREQ(NLINE,2)+H1(M,M)*ISIGN(1,INFOE(LL,1)) 102 CONTINUE C C...evaluate the derivatives for the levels extracted above C 106 NEXTF=INFOE(LL,2) DO 107 LL=1,NCON DO 108 I=1,NCONST CONST(I)=0.0D0 108 CONTINUE CONST(ICFIT(LL))=1.D0 C CALL SETHAM(MS,I1,I2,F,CONST) C DO 109 I=1,NEXTR NLINE=IABS(NELINF(I,1)) M=NELINF(I,2) EE=0.0 DO 110 M1=1,MS DO 111 M2=M1,MS RMULT=2.D0 IF(M1.EQ.M2)RMULT=1.D0 EE=EE+RMULT*H1(M1,M2)*T(M1,M)*T(M2,M) 111 CONTINUE 110 CONTINUE DER(NLINE,LL)=DER(NLINE,LL)+ISIGN(1,NELINF(I,1))*EE 109 CONTINUE 107 CONTINUE C 100 CONTINUE GOTO 104 C 103 WRITE(*,'(1X//1X,3A//)')CHAR(7),CHAR(7), 1 '**** H matrix too large - !!!!CRASH!!!!' C 104 RETURN END c c_____________________________________________________________________________ c FUNCTION MHIGH(I1,I2,F) C C MHIGH returns the dimension of the eigenvalue matrix H for given I1,I2 C and F where: C C I1+I2 I+F C \--- \--- C MHIGH = > > (2J+1) C /--- /--- C I=|I1-I2| J=|I-F| C C eg. - for I1=1, I2=3/2 and Fmax=13/2 there will be: C C I = 1/2: J = 6,7 C I = 3/2: J = 5,6,7,8 C I = 5/2: J = 4,5,6,7,8,9 ---> MHIGH = 168 C C - for I1=3/2, I2=0, Fmax=41/2 ---> MHIGH = 168 C - for I1=3/2, I2=0, Fmax=165/2 (J=80) ---> MHIGH = 664 C..................................................... C Convenient program size estimator: C C PROGRAM SIZE (bytes) = nmax*nmax*12 + (200000 to 400000) C C where: nmax=(2I1+1)(2I2+1)(2F+1), and the 0.2-0.4Mb element C depends on system and FORTRAN used C.................................................... C IMPLICIT REAL*8(A-E,G-H,O-Z),INTEGER(F) c MS=0 IH=I1+I2 IL=IABS(I1-I2) DO 1 I=IL,IH,2 JH=F+I JL=IABS(I-F) DO 1 J=JL,JH,2 MS=(J+1)+MS 1 CONTINUE MHIGH=MS c RETURN END c c_____________________________________________________________________________ c SUBROUTINE QUANTN(I1,I2,F,M,I,J,K) C C This routine returns I,J,K for the M'th row or column of the eigenvalue C matrix for given I1,I2 and F C C NOTE: both I and J are returned in units of 1/2 C C IMPLICIT REAL*8(A-E,G-H,O-Z),INTEGER(F) c M1=M IH=I1+I2 IL=IABS(I1-I2) DO 1 I=IL,IH,2 JH=F+I JL=IABS(I-F) DO 1 J=JL,JH,2 M2=M1-J-1 IF(M2.LE.0)GOTO 2 M1=M2 1 CONTINUE 2 K=M1-1-J/2 C RETURN END c c_____________________________________________________________________________ c SUBROUTINE NUMENL(I1,I2,F,I,J,K,M) C C This routine returns M, the number of eigenvalue among the eigenvalues of C H, for given I,J,K=(Tau+J+1) and I1,I2,F C C NOTE: I1,I2,F and I are expected to be in units of 1/2, whereas C J and K are to be in units of 1 C C IMPLICIT REAL*8(A-E,G-H,O-Z),INTEGER(F) C J2=J*2 IH=I1+I2 IL=IABS(I1-I2) M=0 DO 1 II1=IL,IH,2 JH=II1+F JL=IABS(II1-F) DO 2 JJ1=JL,JH,2 IF(II1.LT.I)GOTO 3 IF(JJ1.LT.J2)GOTO 3 M=M+K GOTO 5 3 M=M+JJ1+1 2 CONTINUE 1 CONTINUE C 5 RETURN END c c_____________________________________________________________________________ c SUBROUTINE SORTC(N,M) c PARAMETER (MDIM=1600) PARAMETER (MDIM=300) COMMON /SORTCC/WK,IPT INTEGER*2 IPT(MDIM) REAL*8 WK(MDIM),EE C C ... This routine sorts the quantities part of vector WK from N to M in C ascending order of magnitude and also accordingly rearranges vector C IPT of pointers to original positions of sorted quantities C DO 101 I=N,M-1 J=I 106 J=J+1 IF(WK(J)-WK(I))103,104,104 103 EE=WK(I) WK(I)=WK(J) WK(J)=EE K=IPT(I) IPT(I)=IPT(J) IPT(J)=K 104 IF(J.EQ.M)GOTO 101 GOTO 106 101 CONTINUE C RETURN END C c_____________________________________________________________________________ c SUBROUTINE SETHAM(MS,I1,I2,F,CON) C C Set up the H matrix of size MSxMS for values of F,I1 and I2 provided, C and constants in CON C C IMPLICIT REAL*8(A-E,G-H,O-Z),INTEGER(F) c PARAMETER (MDIM=1600,NCONST=31) PARAMETER (MDIM= 300,NCONST=31) COMMON /HBLK/H1/HAMBLK/D,G,V,SR1,SR2,CONST DIMENSION CON(NCONST),CONST(NCONST) DIMENSION H1(MDIM,MDIM) DIMENSION D(5),G(5),V(5),SR1(5),SR2(5) ROOT6=DSQRT(6.D0) ROOT23=DSQRT(2.D0/3.D0) C C...Use the constants to set up the following tensors: C C D(5) - spherical spin-spin tensor C G(5) - spherical quadrupole tensor for nucleus with spin I1 C V(5) - spherical quadrupole tensor for nucleus with spin I2 C SR1(5) - spherical spin-rotation tensor for nucleus with spin I1 C SR2(5) - spherical spin rotation tensor for nucleus with spin I2 C C Order of components in all tensors: C(-2),C(-1),C(O),C(1),C(2) C DO 30 I=1,5 D(I)=0.0D0 G(I)=0.0D0 V(I)=0.0D0 SR1(I)=0.0D0 SR2(I)=0.0D0 30 CONTINUE C G(3)=CON(17) G(1)=CON(18)/ROOT6 G(5)=CON(18)/ROOT6 G(2)=ROOT23*CON(19) G(4)=-ROOT23*CON(19) C V(3)=CON(20) V(1)=CON(21)/ROOT6 V(5)=CON(21)/ROOT6 V(2)=ROOT23*CON(22) V(4)=-ROOT23*CON(22) C SR1(3)=CON(23)-CON(25)/3.D0 SR1(1)=CON(24)/ROOT6 SR1(5)=CON(24)/ROOT6 C SR2(3)=CON(26)-CON(28)/3.D0 SR2(1)=CON(27)/ROOT6 SR2(5)=CON(27)/ROOT6 C D(3)=CON(29) D(1)=CON(30)/ROOT6 D(5)=CON(30)/ROOT6 D(2)=ROOT23*CON(31) D(4)=-ROOT23*CON(31) C DO 102 L=1,NCONST CONST(L)=CON(L) 102 CONTINUE C C...Set up H, assuming it is Hermitian C DO 101 L=1,MS DO 101 M=L,MS CALL QUANTN(I1,I2,F,L,IR,JR,KR) CALL QUANTN(I1,I2,F,M,IC,JC,KC) H1(L,M)=HAMIL(I1,I2,F,IR,IC,JR,JC,KR,KC) IF(L.NE.M)H1(M,L)=H1(L,M) 101 CONTINUE C RETURN END c c_____________________________________________________________________________ c FUNCTION HAMIL(I1,I2,F,IR,IC,JR,JC,KR,KC) C C Evaluation of a Hamiltonian element for given I1,I2,F and: C C row I,J,K quantum numbers: IR,JR,KR C column I,J,K quantum numbers: IC,JC,KC C C I1,I2,F,IR,IC,JR,JC are all in units of 1/2 C C IMPLICIT REAL*8(A-E,G-H,O-Z),INTEGER(F) PARAMETER (NCONST=31) COMMON /HAMBLK/D,G,H,SR1,SR2,CON DIMENSION D(5),G(5),H(5),SR1(5),SR2(5),CON(NCONST) REAL*8 J,JPR,K,KPR,IA,IB,Q,I,IPR,FF,JJP1,KSQ c S=0.0D0 ZERO=0.d0 ONE=1.d0 TWO=2.d0 C IF(IR.NE.IC.OR.JR.NE.JC)GOTO 1 JR=JR/2 JC=JC/2 JJP1=JR*(JR+1) KSQ=KR*KR IF(KR.EQ.KC)GOTO 2 IF(KR.EQ.KC+2.OR.KR.EQ.KC-2)GOTO 3 GOTO 7 C C C 2 HBPC=0.5D0*(CON(2)+CON(3)) S=S+HBPC*JJP1+(CON(1)-HBPC)*KSQ-CON(4)*(JJP1*JJP1) 1 -CON(5)*(JJP1*KSQ)-CON(6)*(KSQ*KSQ) 1 +CON(9)*JJP1**3+CON(10)*JJP1*JJP1*KSQ 1 +CON(11)*(JJP1*KSQ*KSQ)+CON(12)*KSQ*KSQ*KSQ GOTO 7 C C C 3 S=S+(0.25D0*(CON(2)-CON(3))-CON(7)*JJP1-0.5D0*CON(8) 1 *((KR+2)**2+KSQ)+CON(13)*JJP1*JJP1 1 +0.5D0*CON(14)*JJP1*((KR+2)**2+KSQ) 1 +0.5D0*CON(15)*((KR+2)**4+KSQ*KSQ)) 2 *DSQRT(DBLE((JJP1-KR*(KR+1))*(JJP1-(KR+1)*(KR+2)))) C C 7 JR=2*JR JC=2*JC C 1 M=KR-KC J=DBLE(JR)*0.5D0 JPR=DBLE(JC)*0.5D0 K=KR KPR=KC IA=DBLE(I1)*0.5D0 IB=DBLE(I2)*0.5D0 I=DBLE(IR)*0.5D0 IPR=DBLE(IC)*0.5D0 FF=DBLE(F)*0.5D0 Q=KC-KR C C M=M+3 C IF(M.GT.5.OR.M.LT.1)S=0.0D0 IF(M.GT.5.OR.M.LT.1)GOTO 9 C C C IF(I1.LE.1)GOTO 5 IF(G(M).EQ.0.0D0)GOTO 5 C N=IC+(F+I1+I2)/2+KR S=S+(-1)**N*0.25D0*DSQRT(DBLE((IR+1)*(IC+1)*(JR+1)*(JC+1))) 1 *SIXJ(I,IPR,TWO,JPR,J,FF)*THREEJ(J,TWO,JPR,-K,-Q,KPR) 3 *(SIXJ(I,IPR,TWO,IA,IA,IB)/THREEJ(IA,TWO,IA,-IA,ZERO,IA)) 4 *G(M) C C C 5 IF(I2.LE.1)GOTO 10 IF(H(M).EQ.0.0D0)GOTO 10 C N=KR+(F+IC+IR+I1+I2)/2 S=S+(-1)**N*0.25D0*DSQRT(DBLE((IR+1)*(IC+1)*(JR+1)*(JC+1))) 1 *SIXJ(I,IPR,TWO,JPR,J,FF)*THREEJ(J,TWO,JPR,-K,-Q,KPR) 3 *(SIXJ(I,IPR,TWO,IB,IB,IA)/THREEJ(IB,TWO,IB,-IB,ZERO,IB)) 4 *H(M) C C C 10 IF(D(M).EQ.0.0D0)GOTO 6 C N=KR+(IC+F)/2 S=S+(-1)**N*DSQRT((2.0D0*I+1.0D0)*(2.0D0*IPR+1.0D0)*(2.0D0*J 1+1.0D0)*30.D0*(2.0D0*JPR+1.0D0)*(2.0D0*IA+1.0D0)*(IA+1.0D0)*IA 2*(2.0D0*IB+1.0D0)*(IB+1.0)*IB) 3*THREEJ(J,TWO,JPR,-K,-Q,KPR)*SIXJ(I,IPR,TWO,JPR,J,FF) 4*X9J(IA,IA,ONE,IB,IB,ONE,I,IPR,TWO)*D(M)/2.0D0 C C C 6 IF(SR1(M).EQ.0.0D0)GOTO 11 C N=IC+KR+(I1+I2+JR+JC+F)/2 S=S+(-1)**N*DSQRT(30.0D0*(2.0D0*J+1.0D0)*(2.0D0*JPR+1.0D0) 1*(2.0D0*J+1.0D0)*(J+1.0D0)*J*(2.0D0*IPR+1.0D0)*(2.0D0*I+1.0D0) 2*(2.0D0*IA+1.0D0)*(IA+1.0D0)*IA)*SIXJ(IA,I,IB,IPR,IA,ONE) 3*SIXJ(I,J,FF,JPR,IPR,ONE)*SIXJ(ONE,TWO,ONE,JPR,J,J) 4*THREEJ(J,TWO,JPR,-K,-Q,KPR)*SR1(M)/2.0D0 C 11 IF(JR.NE.JC.OR.KR.NE.KC)GOTO 8 C.D.Esp IF(CON(25).EQ.0.0D0)GOTO 8 C N=(I1+I2+F+JR)/2+IC+1 S=S-(-1)**N*DSQRT((2.0D0*IPR+1.0D0)*(2.0D0*I+1.0D0) 1*(2.0D0*IA+1.0)*(IA+1.0D0)*IA*(2.0D0*J+1.0D0)*(J+1.0D0) 2*J)*SIXJ(I,J,FF,J,IPR,ONE)*SIXJ(IA,I,IB,IPR,IA,ONE) 3*CON(25)/3.D0 C.D.Esp C C C 8 IF(SR2(M).EQ.0.0D0)GOTO 12 C N=(IC+JR+F+I1+I2+IR+JC)/2+KR S=S+(-1)**N*DSQRT(30.D0*(2.0D0*J+1.0D0)*(2.0D0*JPR+1.0D0) 1 *(2.0D0*J+1.0D0)*(J+1.0D0)*J*(2.0D0*IPR+1.0D0)*(2.0D0*I+1.0D0) 2 *(2.0D0*IB+1.0D0)*(IB+1.0D0)*IB)*SIXJ(I,ONE,IPR,IB,IA,IB) 3 *SIXJ(I,J,FF,JPR,IPR,ONE)*SIXJ(ONE,TWO,ONE,JPR,J,J) 4 *THREEJ(J,TWO,JPR,-K,-Q,KPR)*SR2(M)/2.0D0 C 12 IF(JR.NE.JC.OR.KR.NE.KC)GOTO 9 C.D.Esp IF(CON(28).EQ.0.0D0)GOTO 9 C N=(I1+I2+IR+IC+JR+F)/2+1 S=S-(-1)**N*DSQRT((2.0D0*IPR+1.0D0)*(2.0D0*I+1.0D0) 1 *(2.0D0*IB+1.0D0)*(IB+1.0D0)*IB*(2.0D0*J+1.0D0)*(J+1.0D0) 2 *J)*SIXJ(I,J,FF,J,IPR,ONE)*SIXJ(I,ONE,IPR,IB,IA,IB) 3 *CON(28)/3.D0 C.D.Esp C C 9 IF(ABS(S).LE.1.D-300)S=0.D0 HAMIL=S c RETURN END c c_____________________________________________________________________________ c C C ANGULAR MOMENTUM MATRIX ELEMENT CALCULATIONS C c c_____________________________________________________________________________ c c FUNCTION THREEJ(X1,X2,X3,Y1,Y2,Y3) IMPLICIT REAL*8(A-H,O-Z) c IF(ABS(Y1+Y2+Y3).GT.0.1D0)GOTO 50 TR=TRIANG(X1,X2,X3) IF(NINT(TR).NE.0)GOTO 60 50 THREEJ=0.0D0 RETURN 60 IR=X2+Y3-X1 PP=(-1)**IR QQ=TRI(X1,X2,X3) RR=DSQRT(FAC(X1-Y1))*DSQRT(FAC(X1+Y1))*DSQRT(FAC(X2-Y2)) RR=RR*DSQRT(FAC(X2+Y2))*DSQRT(FAC(X3-Y3))*DSQRT(FAC(X3+Y3)) SS=SUML(X1,X2,X3,Y1,Y2,Y3) THREEJ=QQ*RR*SS*PP c RETURN END c c_____________________________________________________________________________ c FUNCTION SIXJ(X1,X2,X3,Y1,Y2,Y3) IMPLICIT REAL*8(A-H,O-Z) c SIXK=TRIANG(X1,X2,X3)*TRIANG(Y1,X2,Y3)*TRIANG(X1,Y2,Y3)*TRIANG(Y1, 1Y2,X3) IF(NINT(SIXK).EQ.0)GOTO 2 IR=X1+X2+Y1+Y2 PP=(-1)**IR SIXK=PP*TRI(X1,X2,X3)*TRI(Y1,Y2,X3)*TRI(Y1,X2,Y3)*TRI(X1,Y2,Y3) 1 *SUMK(X1,X2,X3,Y1,Y2,Y3) 2 CONTINUE SIXJ=SIXK c RETURN END c c_____________________________________________________________________________ c FUNCTION X9J(X1,X2,X3,X4,X5,X6,X7,X8,X9) IMPLICIT REAL*8(A-H,O-Z) c S=0.0D0 Y1=TRIANG(X1,X2,X3) Y2=TRIANG(X4,X5,X6) Y3=TRIANG(X7,X8,X9) Y4=TRIANG(X1,X4,X7) Y5=TRIANG(X2,X5,X8) Y6=TRIANG(X3,X6,X9) T=Y1*Y2*Y3*Y4*Y5*Y6 IF(T.LT.0.1D0)GOTO 1 Z1=X1+X9 Z2=X6+X2 Z3=X4+X8 N1=MINT(Z1) N2=MINT(Z2) N3=MINT(Z3) IF(N1.NE.N2)GOTO 1 IF(N2.NE.N3)GOTO 1 IF(N1.EQ.0)GOTO 2 DO 100 I=0,10 X0=I S=S+(2.0D0*X0+1.0D0)*SIXJ(X1,X2,X3,X6,X9,X0) 1*SIXJ(X4,X5,X6,X2,X0,X8)*SIXJ(X7,X8,X9,X0,X1,X4) 100 CONTINUE GOTO 1 2 DO 101 I=1,19,2 X0=I/2.0D0 S=S-(2.0D0*X0+1.0D0)*SIXJ(X1,X2,X3,X6,X9,X0) 1*SIXJ(X4,X5,X6,X2,X0,X8)*SIXJ(X7,X8,X9,X0,X1,X4) 101 CONTINUE 1 X9J=S c RETURN END c c_____________________________________________________________________________ c FUNCTION TRIANG(O,P,Q) C C Test for the triangular condition: this is not satisfied if sum of any C two of the three numbers is smaller than the third number C IMPLICIT REAL*8(A-H,O-Z) c Z0=0.0D0 IF((O+P-Q).LT.Z0.OR.(O-P+Q).LT.Z0.OR.(-O+P+Q).LT.Z0)GOTO 1 GOTO 2 1 STRIA=Z0 TRIANG=STRIA RETURN 2 STRIA=1.0D0 TRIANG=STRIA c RETURN END c c_____________________________________________________________________________ c FUNCTION TRI(U,V,W) IMPLICIT REAL*8(A-H,O-Z) C Z1=1.0D0 TRK=DSQRT((FAC(U+V-W)*FAC(U-V+W)*FAC(-U+V+W))/(FAC(U+V+W+Z1))) TRI=TRK C RETURN END c c_____________________________________________________________________________ c FUNCTION FAC(Z) IMPLICIT REAL*8(A-H,O-Z) c X=Z Z1=1.0D0 FAC=Z1 1 IF(X.LE.Z1)GOTO 3 2 FAC=FAC*X X=X-Z1 GOTO 1 3 CONTINUE c RETURN END c c_____________________________________________________________________________ c FUNCTION SUML(U,V,W,X,Y,Z) IMPLICIT REAL*8(A-H,O-Z) c DIMENSION FF(5) FF(1)=U-X FF(2)=V+Y FF(3)=U+V-W FF(4)=W-V+X FF(5)=W-U-Y DO 75 IK=2,3 IF(FF(IK).GE.FF(1))GOTO 75 FF(1)=FF(IK) 75 CONTINUE KMAX=FF(1) IF(FF(4).LE.FF(5))GOTO 80 FF(4)=FF(5) 80 IF(FF(4).LE.0.0D0)GOTO 85 KMIN=0 GOTO 95 85 KMIN=-FF(4) 95 SUMQ=0.0D0 DO 97 KF=KMIN,KMAX XY=(-1)**KF RK=KF SUMQ=SUMQ+XY/(FAC(U-X-RK)*FAC(W-V+X+RK)*FAC(V+Y-RK)* 1FAC(W-U-Y+RK)*FAC(RK)*FAC(U+V-W-RK)) 97 CONTINUE SUML=SUMQ c RETURN END c c_____________________________________________________________________________ c FUNCTION SUMK(S1,S2,S3,T1,T2,T3) IMPLICIT REAL*8(A-H,O-Z) DIMENSION FA(7) c Z1=1.0D0 FA(1)=S1+S2+T1+T2+Z1 FA(2)=S1+S2-S3 FA(3)=T1+T2-S3 FA(4)=S1+T2-T3 FA(5)=T1+S2-T3 FA(6)=-S1-T1+S3+T3 FA(7)=-S2-T2+S3+T3 DO 1 IK=2,5 IF(FA(IK).LT.FA(1))GOTO 2 GOTO 1 2 FA(1)=FA(IK) 1 CONTINUE KFMAX=FA(1) IF(FA(6).LE.FA(7))GOTO 3 4 FA(6)=FA(7) Z0=0.0D0 3 IF(FA(6).LE.Z0)GOTO 5 GOTO 6 5 KFMIN=-FA(6) GOTO 8 6 KFMIN=0 8 SUMP=Z0 DO 7 KF=KFMIN,KFMAX XY=(-1)**KF RK=KF Z1=1.0D0 7 SUMP=SUMP+XY*FAC(S1+S2+T1+T2+Z1-RK)/(FAC(S1+S2-S3-RK)*FAC(T1+T2 1 -S3-RK)*FAC(S1+T2-T3-RK)*FAC(T1+S2-T3-RK)*FAC(-S1-T1+S3+T3+RK) 1 *FAC(-S2-T2+S3+T3+RK)*FAC(RK)) SUMK=SUMP c RETURN END c c_____________________________________________________________________________ c FUNCTION MINT(A) IMPLICIT REAL*8(A-H,O-Z) C C This is apparently a test for whether A is close to an integer (MINT=1) C or half integer (MINT=0) C I=2.00001*A N=(1+(-1)**I)/2 MINT=N C RETURN END C c_____________________________________________________________________________ C C c_____________________________________________________________________________ c SUBROUTINE LSQR(IN,LX,SD,X) C C Least squares routine from A.S.G.'s programs CDFIT C (ZK: note that the NINFIT variable was added to eliminate the problem c of counting unfitted lines in calculating the standard deviation) C C IN - number of equations C LX - number of parameters C SD - standard deviations on the fitted parameters C X - calculated increments to parameters C R - correlation matrix C IMPLICIT REAL*8(A-E,G-H,O-Z),INTEGER(F) PARAMETER (NCONST=31,NLIN=200,NLIN2=2*NLIN) COMMON /BLKEN/TFREQ,DER,LABO,INFOE,NELEV,ICFIT COMMON /BLKINV/A/BLKFIT/LINFIT,R DIMENSION TFREQ(NLIN,3),DER(NLIN,NCONST),LABO(NLIN,2,4), * ICFIT(NCONST),INFOE(NLIN2,2),LINFIT(NLIN) DIMENSION A(NCONST,NCONST),R(NCONST,NCONST),C(NCONST), * SD(NCONST),X(NCONST),BC(NLIN),Z(NCONST) c SIGMA=0.0D0 ninfit=0 DO 110 I=1,IN BC(I)=TFREQ(I,3) IF(LINFIT(I).EQ.1)then SIGMA=SIGMA+BC(I)*BC(I) ninfit=ninfit+1 endif 110 CONTINUE if(ninfit-LX.gt.0)then SIGMA=DSQRT(SIGMA/(ninfit-LX)) else SIGMA=DSQRT(SIGMA) endif C IFAIL=0 DO 101 I=1,NCONST DO 101 J=1,NCONST R(I,J)=0.0D0 101 CONTINUE C C...Form: A = DER' * DER and rescale it to improve conditioning C DO 50 I1=1,LX DO 3 I2=I1,LX SUM=0.0D0 DO 5 IS=1,IN IF(LINFIT(IS).NE.1)GOTO 5 SUM=SUM+DER(IS,I1)*DER(IS,I2) 5 CONTINUE A(I1,I2)=SUM 3 CONTINUE Z(I1)=DSQRT(1.D0/A(I1,I1)) 50 CONTINUE C DO 51 I1=1,LX DO 51 I2=I1,LX A(I1,I2)=A(I1,I2)*Z(I1)*Z(I2) A(I2,I1)=A(I1,I2) 51 CONTINUE C C...Form: C = DER' * BC C DO 6 I1=1,LX SUM=0.0D0 DO 8 IS=1,IN IF(LINFIT(IS).NE.1)GOTO 8 SUM=SUM+DER(IS,I1)*BC(IS) 8 CONTINUE C(I1)=SUM 6 CONTINUE C C...Invert A ie. (DER'*DER) and 'unscale' C CALL INV(LX) C DO 52 I1=1,LX DO 52 I2=1,LX A(I1,I2)=A(I1,I2)*Z(I1)*Z(I2) 52 CONTINUE C C...Increments to constants X = A * C, standard errors SD = SIGMA * A**1/2 C and the correlation matrix R C DO 9 I1=1,LX SUM=0.0D0 DO 10 IS=1,LX SUM=SUM+A(I1,IS)*C(IS) 10 CONTINUE X(I1)=SUM SD(I1)=SIGMA*DSQRT(A(I1,I1)) 9 CONTINUE C SDEVSQ=SIGMA*SIGMA DO 15 I1=1,LX DO 15 I2=1,I1 R(I1,I2)=SDEVSQ*A(I1,I2)/(SD(I1)*SD(I2)) 15 CONTINUE c RETURN END c c_____________________________________________________________________________ c SUBROUTINE HDIAG(N,NDIM,IEGEN,IORD) C C A=MATRIX TO BE DIAGONALIZED C NDIM=DIMENSION OF A (Dummy parameter here since COMMON's are used) C N=DIMENSION OF SUBMATRIX TO BE DIAGONALIZED C IEGEN=0 EIGENVALUES AND EIGENVECTORS C =1 EIGENVALUES ONLY C IORD=0 NO ORDERING OF EIGENVALUES OR EIGENVECTORS (ORDER IN = ORDER OUT) C =1 ORDERING BY SIZE OF EIGENVALUES C EIVR=TRANSFORMATION MATRIX (MATRIX OF EIGENVECTORS) C C THIS ROUTINE USES A VARIABLE THRESHOLD JACOBI METHOD C IMPLICIT REAL*8(A-H,O-Z) c PARAMETER (MDIM=1600) PARAMETER (MDIM= 300) COMMON /HBLK/A/SBLK/EIVR DIMENSION A(MDIM,MDIM) REAL*4 EIVR(MDIM,MDIM) c write(*,'('' starting '',i3,''x'',i3)')n,n IF(N.GT.0)GOTO 1 EIVR(1,1)=1.0 RETURN 1 IF(IEGEN.GT.0)GOTO 102 DO 101 J=1,N DO 100 I=1,N 100 EIVR(I,J)=0.0 101 EIVR(J,J)=1.0 C C FIND THE ABSOLUTELY LARGEST ELEMENT OF A C 102 ATOP=0.0D0 DO 111 I=1,N DO 111 J=1,N IF(ATOP.GE.DABS(A(I,J)))GOTO 111 ATOP=DABS(A(I,J)) 111 CONTINUE IF(ATOP)109,109,113 109 RETURN C C CALCULATE THE STOPPING CRITERION -- DSTOP C 113 AVGF=(N*(N-1))*0.55D0 D=0.0D0 DO 114 JJ=2,N DO 114 II=2,JJ S=A(II-1,JJ)/ATOP 114 D=S*S+D DSTOP=(1.0D-06)*D C C CALCULATE THE THRESHOLD -- THRSH C THRSH=DSQRT(D/AVGF)*ATOP C C START A SWEEP C 115 IFLAG=0 DO 130 JCOL=2,N JCOL1=JCOL-1 DO 130 IROW=1,JCOL1 AIJ=A(IROW,JCOL) C C COMPARE THE OFF DIAGONAL ELEMENT WITH THRSH C IF(DABS(AIJ).LE.THRSH)GOTO 130 AII=A(IROW,IROW) AJJ=A(JCOL,JCOL) S=AJJ-AII C C CHECK TO SEE IF THE CHOSEN ROTATION IS LESS THAN THE ROUNDING ERROR C IF SO, THEN DO NOT ROTATE C IF(DABS(AIJ).LE.(1.0D-09*DABS(S)))GOTO 130 IFLAG=1 C C IF THE ROTATION IS VERY CLOSE TO 45 DEGREES, SET SIN AND COS C TO 1/(ROOT 2) C IF((1.0D-10*DABS(AIJ)).LT.DABS(S))GOTO 116 S=0.70710678118654752D0 C=S GOTO 120 C C CALCULATE SIN AND COS FOR ROTATION THAT IS NOT VERY CLOSE C TO 45 DEGREES C 116 T=AIJ/S S=0.25D0/DSQRT(0.25D0+T*T) C C COS=C, SIN=S C C=DSQRT(0.5D0+S) S=2.0D0*T*S/C C C CALCULATION OF THE NEW ELEMENTS OF MATRIX A C 120 DO 121 I=1,IROW T=A(I,IROW) U=A(I,JCOL) A(I,IROW)=C*T-S*U 121 A(I,JCOL)=S*T+C*U I2=IROW+2 IF(I2.GT.JCOL)GOTO 123 CONTINUE DO 122 I=I2,JCOL T=A(I-1,JCOL) U=A(IROW,I-1) A(I-1,JCOL)=S*U+C*T 122 A(IROW,I-1)=C*U-S*T 123 A(JCOL,JCOL)=S*AIJ+C*AJJ A(IROW,IROW)=C*A(IROW,IROW)-S*(C*AIJ-S*AJJ) DO 124 J=JCOL,N T=A(IROW,J) U=A(JCOL,J) A(IROW,J)=C*T-S*U 124 A(JCOL,J)=S*T+C*U C C ROTATION COMPLETED C SEE IF EIGENVECTORS REQUIRED C IF(IEGEN.GT.0)GOTO 126 DO 125 I=1,N T=EIVR(I,IROW) EIVR(I,IROW)=SNGL(C*T-EIVR(I,JCOL)*S) 125 EIVR(I,JCOL)=SNGL(S*T+EIVR(I,JCOL)*C) C C CALCULATE THE NEW NORM D AND COMPARE WITH DSTOP C 126 CONTINUE S=AIJ/ATOP D=D-S*S IF(D.GE.DSTOP)GOTO 129 C C RECALCULATE DSTOP AND THRSH TO DISCARD ROUNDING ERRORS C D=0.0D0 DO 128 JJ=2,N DO 128 II=2,JJ S=A(II-1,JJ)/ATOP 128 D=S*S+D DSTOP=(1.0D-06)*D 129 THRSH=DSQRT(D/AVGF)*ATOP 130 CONTINUE IF(IFLAG.NE.0)GOTO 115 C C ARRANGE THE EIGENVALUES IN THE ORDER OF INCREASING ENERGY C ARRANGE THE EIGENVECTORS IN THE SAME ORDER C IF(IORD.EQ.0)RETURN NU=N DO 11 I=1,N IF(I.GE.NU)RETURN AMIN=A(I,I) DO 10 J=I,NU IF(A(J,J).GE.AMIN)GOTO 10 C C IF IEGEN IS -1 EXCLUDE UNCONVERGED EIGENVALUES FROM THIS ORDERING C TE=ABS(EIVR(N,J))+ABS(EIVR(N-1,J)) IF((TE.GT.0.05D0).AND.(IEGEN.EQ.-1))GOTO 15 II=I AMIN=A(J,J) A(J,J)=A(I,I) A(I,I)=AMIN 16 DO 12 K=1,N TEMP=EIVR(K,II) EIVR(K,II)=EIVR(K,J) 12 EIVR(K,J)=SNGL(TEMP) GOTO 10 15 AM=A(J,J) A(J,J)=A(NU,NU) A(NU,NU)=AM II=NU NU=NU-1 GOTO 16 10 CONTINUE 11 CONTINUE c RETURN END c c_____________________________________________________________________________ C SUBROUTINE INV(N) C C INVERSION OF A MATRIX (Basically SSP vintage) C PARAMETER (NCONST=31,NSQCON=NCONST*NCONST) COMMON /BLKINV/AA REAL*8 AA(NCONST,NCONST),A(NSQCON) INTEGER*2 L(NCONST),M(NCONST) REAL*8 D,BIGA,HOLD C IJ=0 DO 5000 J=1,N DO 5000 I=1,N IJ=IJ+1 5000 A(IJ)=AA(I,J) C D=1.0D0 NK=-N DO 80 K=1,N NK=NK+N L(K)=K M(K)=K KK=NK+K BIGA=A(KK) DO 20 J=K,N IZ=N*(J-1) DO 20 I=K,N IJ=IZ+I 10 IF(DABS(BIGA)-DABS(A(IJ))) 15,20,20 15 BIGA=A(IJ) L(K)=I M(K)=J 20 CONTINUE C J=L(K) IF(J-K) 35,35,25 25 KI=K-N DO 30 I=1,N KI=KI+N HOLD=-A(KI) JI=KI-K+J A(KI)=A(JI) 30 A(JI)=HOLD C 35 I=M(K) IF(I-K) 45,45,38 38 JP=N*(I-1) DO 40 J=1,N JK=NK+J JI=JP+J HOLD=-A(JK) A(JK)=A(JI) 40 A(JI)=HOLD C 45 IF(BIGA) 48,46,48 46 D=0.0D0 WRITE(*,'(1X/'' ***** ZERO DETERMINANT *****''/)') RETURN 48 DO 55 I=1,N IF(I-K) 50,55,50 50 IK=NK+I A(IK)=A(IK)/(-BIGA) 55 CONTINUE C DO 65 I=1,N IK=NK+I HOLD=A(IK) IJ=I-N DO 65 J=1,N IJ=IJ+N IF(I-K) 60,65,60 60 IF(J-K) 62,65,62 62 KJ=IJ-I+K A(IJ)=HOLD*A(KJ)+A(IJ) 65 CONTINUE C KJ=K-N DO 75 J=1,N KJ=KJ+N IF(J-K)70,75,70 70 A(KJ)=A(KJ)/BIGA 75 CONTINUE C D=D*BIGA C A(KK)=1./BIGA 80 CONTINUE C K=N 100 K=K-1 IF(K) 150,150,105 105 I=L(K) IF(I-K) 120,120,108 108 JQ=N*(K-1) JR=N*(I-1) DO 110 J=1,N JK=JQ+J HOLD=A(JK) JI=JR+J A(JK)=-A(JI) 110 A(JI)=HOLD 120 J=M(K) IF(J-K) 100,100,125 125 KI=K-N DO 130 I=1,N KI=KI+N HOLD=A(KI) JI=KI-K+J A(KI)=-A(JI) 130 A(JI)=HOLD GOTO 100 C 150 DO 5001 I=1,N DO 5001 J=1,N IJ=N*(J-1)+I 5001 AA(I,J)=A(IJ) C RETURN END c c_____________________________________________________________________________ c_____________________________________________________________________________