$DEBUG C----------------------------------------------------------------------------- C C SYMTOP - This is a simple calculator of tables of symmetric top C transition frequencies and of their relative intensities C c A useful feature is calculation of uncertainties in predicted c frequencies by using both uncertainties in spectroscopic c constants and intercorrelations between them. C C Ver 21.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 Nov 1992: rescued from oblivion by porting from PDP-11 and 8 inch floppies C Feb 2001: dragged into the XXI-st century C C----------------------------------------------------------------------------- C C Input is from a disk data file - reedit a previous version, which should c be self-explanatory. c c - do not not cut out any of the dividing lines c - dots in the dividing lines serve to indicate which of the columns c beneath are read on input c - if correlation coefficients are not available the relevant block will C C----------------------------------------------------------------------------- C C IMPLICIT REAL*8 (A-H,O-Z) c parameter (nconst=6) c REAL*8 JOPT,IJOPT,CONST(nconst),ECONST(nconst), * VARCOV(nconst,nconst),intens CHARACTER*8 FORMAR(8),FDESCR,PEDESC,EDESCR CHARACTER FOROUT*64,unit*3 CHARACTER FILout*40,FILinp*40,cdummy*20,errmes*60,coment*78 REAL*4 CONNAM(nconst),INTERM(nconst) EQUIVALENCE (B,CONST(1)) EQUIVALENCE(FOROUT,FORMAR(1)) COMMON CONST,ECONST,VARCOV,A,Temp,SPIN C NLINES=72 DATA CONNAM/'B ','DJ ','DJK ','HJ ','HJK ','HKJ '/ DATA FORMAR(1),FORMAR(8),FDESCR,PEDESC,EDESCR/'(1X,A4',')', 1 ',F12.4',',1PE12.4',',E12.4'/ C WRITE(*,3344) 3344 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | S Y M T O P - Calculation of symmetric top transition', * ' frequencies',T79,'|'/ * ' | and their relative intensities', * T79,'|'/ * ' |',76(1H_),'|'/' version 21.II.2001',T64,'Zbigniew KISIEL'/) c C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C INPUT C 752 WRITE(*,751) 751 FORMAT(1X/' Data file name: ',$) READ(*,'(A)',ERR=752)FILinp OPEN(4,FILE=FILinp,ERR=752,STATUS='OLD') NIN=4 c errmes='Cannot read the main comment' READ(NIN,'(a)',ERR=50)COMENT write(*,'(1x/1x,a/)')coment C DO 570 I=1,NCONST DO 570 J=1,NCONST VARCOV(I,J)=0.0D0 570 CONTINUE C C...constants and errors (for historical reasons internal calculations C are in kHz) C errmes='Cannot read one of the constants' read(nin,'(1a)')cdummy DO 2 N=1,NCONST READ(NIN,3,ERR=50)cdummy,CONST(N),econst(n) WRITE(*,4)cdummy,CONST(n),econst(n) if(n.eq.1)then const(n)=const(n)*1.d+3 econst(n)=econst(n)*1.d+3 endif if(n.ge.4.and.n.le.6)then const(n)=const(n)*1.d-6 econst(n)=econst(n)*1.d-6 endif VARCOV(N,N)=ECONST(N)**2 2 CONTINUE 3 FORMAT(a,F25.0,5x,f25.0) 4 FORMAT(1x,a,F20.12,' +-',F20.12) C C...correlation coefficients C errmes='Bad correlation coefficients flag' read(nin,'(1a)')cdummy READ(NIN,1002,ERR=50)IFLAG 1002 format(35x,i15) read(nin,'(1a)')cdummy IF(IFLAG.NE.1)then do 5 i=1,6 read(nin,'(1a)')cdummy 5 continue GO TO 501 endif c errmes='Problems reading correlation coefficients' read(nin,'(1a)')cdummy DO 1502 I=2,NCONST READ(NIN,1503,ERR=50)(varcov(I,J),J=1,I-1) do 1504 j=1,i-1 VARCOV(J,I)=VARCOV(I,J)*DSQRT(VARCOV(I,I)*VARCOV(J,J)) 1504 continue 1502 CONTINUE 1503 format(5x,10F8.0) C C...A, temperature, frequency limits and K limit C 501 read(nin,'(1a)')cdummy errmes='Reading A' READ(NIN,1003,ERR=50)A 1003 format(35x,F15.0) A=A*1.D-3 c errmes='Reading temperature' READ(NIN,1003,ERR=50)Temp C errmes='Reading lower frequency limit' READ(NIN,1003,ERR=50)FMIN C errmes='Reading upper frequency limit' READ(NIN,1003,ERR=50)FMAX c FMAX=FMAX*1.D6 FMIN=FMIN*1.D6 JMIN=(FMIN*0.5D0/B)-1.D0 IF(JMIN.LT.0)JMIN=0 JMAX=(FMAX*0.5D0/B)+1.D0 c errmes='Reading upper limit on K' READ(NIN,1002,ERR=50)KMAX C errmes='Reading spin for equivalent nuclei' READ(NIN,1003,ERR=50)SPIN C CLOSE(NIN) C C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C OUTPUT C 2000 WRITE(*,668) 668 FORMAT(1X/' File name for output: ',$) READ(*,'(A)',ERR=2000)FILout OPEN(3,FILE=FILout,STATUS='UNKNOWN',err=2000) C WRITE(3,11)COMENT 11 FORMAT(1X,78(1H-)/1X,a/1X,78(1H-)/) DO 8 N=1,NCONST cmult=1.d0 unit='kHz' if(n.eq.1)then cmult=1.D-3 unit='MHz' endif if(n.ge.4.and.n.le.6)then cmult=1.D+6 unit='mHz' endif WRITE(3,12)CONNAM(N),CONST(N)*cmult,ECONST(N)*cmult,unit 8 CONTINUE 12 FORMAT(1X,A4,'=',F21.10,' +- ',F16.10,1x,a) c c...Zeroth order approximation for JOPT: Gordy&Cook 2nd ed., eq.5.41 c JOPT=5.5D0*DSQRT(1.D6*Temp/B) c c...crude iterative refinement of JOPT c step=0.1d0 3000 ijopt=intens(b,temp,jopt) djopt=intens(b,temp,jopt+step) if(djopt-ijopt.gt.0.d0)then jopt=jopt+step goto 3000 else if(abs((ijopt-djopt)/ijopt).gt.1.d-8)then step=-step/5. goto 3000 endif c WRITE(3,9)A*1.E+3,temp,jopt 9 FORMAT(/' A = ',F14.4,' MHz'/ * ' Temp= ',F14.4,' K'/ * ' Jopt= ',F14.4) if(spin.lt.0.)then write(3,'(22x,''Not a simple C3v symmetric top'')') else write(3,'('' Spin= '',F14.4, * '' for off-axis nuclei in a C3v symmetric top'')')spin endif c IF(IFLAG.NE.1)GO TO 505 WRITE(3,506) 506 FORMAT(/' Variance covariance matrix (kHz**2) and ', 1 'correlation coefficients :'/) NCON=0 DO 601 N=1,NCONST NCON=NCON+1 INTERM(NCON)=CONNAM(N) 601 CONTINUE WRITE(3,507)(INTERM(I),I=1,NCON) 507 FORMAT(3X,6(8x,A4)) WRITE(3,508) 508 FORMAT(1X) FORMAR(2)=PEDESC DO 3700 N=3,7 FORMAR(N)=EDESCR 3700 CONTINUE c DO 509 I=1,NCONST NCON=0 DO 600 J=1,NCONST NCON=NCON+1 INTERM(NCON)=VARCOV(I,J) 600 CONTINUE WRITE(3,FOROUT)CONNAM(I),(INTERM(J),J=1,NCON) IF(I.GE.NCONST)GO TO 509 FORMAR(I+1)=FDESCR FORMAR(I+2)=PEDESC 509 CONTINUE C 505 WRITE(3,60) 60 FORMAT(//' Symmetric top transition frequencies /MHz, their', 1 ' uncertainties /MHz,'/' and relative intensities:'// 2 ' _ intensity relative to intensity at', 3 ' J.opt'/' |') c C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C CALCULATION C WRITE(*,41) 41 FORMAT(//' **** CALCULATION IN PROGRESS, now working ', 1 'on J ='//) C C...Main loop for TRANSITION FREQUENCIES AND RELATIVE INTENSITIES C DO 29 J=JMIN,JMAX,3 C WRITE(*,777)J,J+1,J+2 777 FORMAT(1H+,3I4) if(j.ne.jmin)write(3,'(1x/)') C RIKZ =intens(b,temp,dble(j)) RIKZA=intens(b,temp,dble(j+1)) RIKZB=intens(b,temp,dble(j+2)) RIJ =RIKZ /IJOPT RIJA=RIKZA/IJOPT RIJB=RIKZB/IJOPT C 666 WRITE(3,45)J+1,J,RIJ,J+2,J+1,RIJA,J+3,J+2,RIJB 45 FORMAT(' J =',I3,' <-',I3,3X,'(',F5.3,')',T29, 1 'J =',I3,' <-',I3,3X,'(',F5.3,')',T56,'J =', 1 I3,' <-',I3,3X,'(',F5.3,')'/) C IF(J.EQ.JMIN)THEN WRITE(3,755) 755 FORMAT(' _ K quantum number'/ 1' | _ calculated frequency'/ 1' | _ estimated error in prediction'/ 1' | _ intensity relative to 2*I(K=0) ', 1 'for this J'/ 1' |') ENDIF C CALL FREINT(J,0,F1,ERF1,RIk01) CALL FREINT(J+1,0,F2,ERF2,RIk02) CALL FREINT(J+2,0,F3,ERF3,RIk03) DO 20 K=0,J+2 IF(K.EQ.KMAX+1)GO TO 29 CALL FREINT(J,K,F1,ERF1,RI1) CALL FREINT(J+1,K,F2,ERF2,RI2) CALL FREINT(J+2,K,F3,ERF3,RI3) RI1=RI1/(2.*RIK01) RI2=RI2/(2.*RIK02) RI3=RI3/(2.*RIK03) IF(K.EQ.J+1)THEN WRITE(3,27)K,F2,ERF2,RI2,K,F3,ERF3,RI3 GO TO 20 ENDIF IF(K.EQ.J+2)THEN WRITE(3,28)K,F3,ERF3,RI3 GO TO 20 ENDIF kk=k if(k.ge.100)kk=k-100 WRITE(3,25)KK,F1,ERF1,RI1,K,F2,ERF2,RI2,K,F3,ERF3,RI3 20 CONTINUE C 29 CONTINUE C 25 FORMAT(I2,F12.3,F6.3,F5.2,2(1X,I3,F12.3, 1 F6.3,F5.2)) 27 FORMAT(25X,2(1X,I3,F12.3,F6.3,F5.2)) 28 FORMAT(53X,I3,F12.3,F6.3,F5.2) C 30 write(3,'(/1x,78(1H-))') CLOSE(3) write(*,'(1x)') C 51 STOP c 50 write(*,'(1x//'' **** INPUT ERROR: '',a//)')errmes stop c END C C----------------------------------------------------------------------------- C SUBROUTINE FREINT(J,K,F,ERF,RI) c c...Frequency, error and intensity of a symmetric top transition C IMPLICIT REAL*8(A-H,O-Z) parameter (nconst=6) c COMMON CONST,ECONST,VARCOV,A,Temp,SPIN EQUIVALENCE (B,CONST(1)) REAL*8 CONST(nconst),ECONST(nconst),VARCOV(nconst,nconst), * COEF(nconst) C X=J+1 RJ=J ZK=K**2 Y=2.D0*X**2 COEF(1)=1.D0 COEF(2)=-Y COEF(3)=-ZK COEF(4)=Y*0.25D0*((RJ+2.D0)**3-RJ**3) COEF(5)=Y*ZK COEF(6)=ZK*ZK F=0.D0 ERF=0.D0 DO 2 LI=1,NCONST F=F+CONST(LI)*COEF(LI) DO 3 LJ=1,NCONST IF(LI.GT.LJ)GO TO 3 TERM=VARCOV(LI,LJ)*COEF(LI)*COEF(LJ) IF(LI.NE.LJ)TERM=2.D0*TERM ERF=ERF+TERM 3 CONTINUE 2 CONTINUE F=2.D-3*X*F ERF=2.D-3*DSQRT(ERF)*X RI=(1.D0-ZK/X**2)*DEXP(-4.8D-2*(A-B*1.D-6)*ZK/Temp) c IF(K.ne.0)RI=RI*2.d0 IF(SPIN.lt.0)goto 1 c IF((K/3)*3.EQ.K)then swt=(2.*spin+1)*(4.*spin*spin+4.*spin+3.)/3. else swt=(2.*spin+1)*(4.*spin*spin+4.*spin)/3. endif RI=RI*swt C 1 RETURN END C C----------------------------------------------------------------------------- c FUNCTION intens(b,temp,j) real*8 b,temp,j,intens,nuzero c c...Intensity of K=0 transition Gordy&Cook 2nd ed., eq.6.61 with X set to 1 c (note that this will begin failing for very low T, giving rise to -ve c intensities) c nuzero=2.d0*B*(J+1.D0)*1.D-6 intens=nuzero**3*(1.D0-0.024D0*nuzero/Temp)* 1 DEXP(-0.024D0*J*nuzero/Temp) if(intens.eq.0.d0)intens=1.d-10 return end c C-----------------------------------------------------------------------------