$DEBUG C$LARGE C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C ASROTL - ASYMMETRIC ROTOR ENERGY LEVELS, TRANSITIONS, LINE STRENGTHS AND C INTENSITIES C C NOTE: in this version the intensities for representation III are C only approximate C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C - J UP TO 200 (MAXJ) C - WATSON'S REDUCED HAMILTONIAN (ANY COMBINATION OF REDUCTION 'A' OR C 'S' AND REPRESENTATIONS I OR III) C - UP TO AND INCLUDING DECADIC TERMS C - UP TO 8150 LINES C - ALL TYPES OF TRANSITIONS AND BROAD RANGE OF SELECTION RULES ARE C SUPPORTED, CLEAR INDICATION OF TRANSITION TYPE IS GIVEN IN THE C PRINTOUT C - LINE INTENSITIES (MAXIMUM ABSORPTION COEFFICIENTS) ARE CALCULATED C FOR A SPECIFIED TEMPERATURE, VALUES ASSUMED FOR OTHER KEY C PARAMETERS ARE PRINTED AT THE END OF THE FREQUENCY LISTING C C C This program, though much developed, is still recognisably descended C from ASROT of P.J.Mjoberg et al (ca 1972) - handling of Wang C factorisation is common with that used in ASFIT. C The program is standard F77 and compiles with many FORTRANs including C MSF, SUN, SG, NDP. C_____________________________________________________________________________ C C ver. 3/94 ----- Zbigniew KISIEL ----- C C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, 02-668 Warszawa, POLAND | C | kisiel@ifpan.edu.pl kisiel@planif61.bitnet | C_________________________/-------------------------------------------------- C C C ASROT can be used completely interactively, questions being asked by the C program being (hopefully) self explanatory. C C In cases where many constants are to be put in, online typing-in C may not be preferred and it is possible to read these from an ASFITL C data file or from a constants file having following format: C C Line.1: NCON = 8 C Line.2: A = 48553.1122 C . B = 5010.04024 C . C = 4529.89707 C DJ = .00230008 C DJK = -.0691118 C DK = 1.74704 C dJ = .000477929 C dK = .029516 C C This file is read in under a fixed (A7,F30) FORMAT, the alphanumeric C field being a dummy field for constant identification. NCON can be set C as required, but no constants can be skipped ie. these have to be in C the following sequence: C C 'A' reduction: C C A, B, C C DJ, DJK, DK, dJ, dK C HJ, HJK, HKJ, HK, hJ, hJK, hK C LJ, LJJK, LJK, LKKJ, LK, lJ, lJK, lKJ, lK C PJ, PJJK, PJK, PKJ, PKKJ, PK, pJ, pJJK, pJK, pKKJ, pK C C 'S' reduction: C C A, B, C C DJ, DJK, DK, d1, d2 C HJ, HJK, HKJ, HK, h1, h2, h3 C LJ, LJJK, LJK, LKKJ, LK, l1, l2, l3, l4 C PJ, PJJK, PJK, PKJ, PKKJ, PK, p1, p2, p3, p4, p5 C C C - On interactive input the cutoff in constants is specified by means C of the order of the Hamiltonian eg. 4 for up to and including quartic, C 6 - sextic etc. Any other values are truncated to their nearest lower C multiple of 2. C C - Various printout combinations are possible (see appropriate question C in the interactive mode), options 1 (frequencies only)+sorting+discard, C and option 3 (energy levels only) being recommended for routine use C C - ASROT has built in detection of close K-1, or K+1 doublets on sorting; C when these are closer than 0.001 in frequency (FRESOL in SORT0) an C entry for a doublet is printed - this has a D after the frequency and C only the common K quantum number is printed. C C - ASROT writes identified transitions to a temporary file, each line C taking 24 bytes. At the present line maximum of 8150 lines this C implies a maximum of ca 200Kbytes for this file. The ASCII output C file may be considerably larger. C C - if ASROT is run under options: 'frequencies only, lines to be sorted, C unsorted lines to be discarded' then the results can be viewed C graphically with program ASPEC/ASC C C - On repetitive use of ASROT for similar predictions interactive use C becomes tedious and many systems allow a pipeline mechanism such as C ASROTL reduction 'A', IRED=+-1 -> reduction 'S' C +ve values -> representation I, -ve values -> representation III C IRED=2 IF(IP.LT.0)THEN IRED=1 DO 3400 I=1,NCONST ISTP(I)=ISTPS(I) 3400 CONTINUE ENDIF IF(IABS(IP).GT.10)THEN IRED=-IRED IP=IABS(IP)-10 namrep='III' ENDIF IP=IABS(IP) C 2901 WRITE(*,4301) 4301 FORMAT(/' JMIN = ',$) READ(*,*,ERR=2901)JMIN 2903 WRITE(*,4302) 4302 FORMAT(' JRMAX = ',$) READ(*,*,ERR=2903)JRMAX 2904 WRITE(*,4303) 4303 FORMAT(' JQMAX = ',$) READ(*,*,ERR=2904)JQMAX 200 WRITE(*,201) 201 FORMAT(' K-1.MAX = ',$) READ(*,*,ERR=200)KMMAX C 2905 WRITE(*,4304) 4304 FORMAT(/' LINES TO BE SORTED (1=YES, 0=NO) ? ',$) READ(*,40,ERR=2905)IS IF(IS.LT.0.OR.IS.GT.1)GOTO 2905 ISON=0 IF(IP.NE.1.OR.IS.NE.1)GOTO 4390 2906 WRITE(*,4391) 4391 FORMAT(/' ARE UNSORTED FREQUENCIES TO BE DISCARDED ? ',$) READ(*,40,ERR=2906)ISON IF(ISON.LT.0.OR.ISON.GT.1)GOTO 2906 4390 WRITE(*,4305) 4305 FORMAT(/' COMMENT :') READ(*,4000,ERR=4390)COMENT 4000 FORMAT(72A1) WRITE(4,42)COMENT WRITE(4,43)namred(iabs(ired)),namrep DO 5848 I=1,NCONST A(I)=0.D0 5848 CONTINUE C 2907 WRITE(*,6853) 6853 FORMAT(/' ARE ROTATIONAL CONSTANTS AVAILABLE ON DISK ? ',$) READ(*,40,ERR=2907)I IF(I.NE.1)GOTO 6854 2908 WRITE(*,6855) 6855 FORMAT(/' NAME OF FILE CONTAINING ROTATIONAL CONSTANTS : ',$) 2909 READ(*,'(A)',ERR=2908)FILROT OPEN(3,FILE=FILROT,STATUS='OLD',ERR=2908) READ(3,'(A)')LINE READ(LINE,6866,ERR=2950)BLUFF,NCON 6866 FORMAT(A7,I5) GOTO 2951 C C...Constants from a .ASL file C 2950 READ(3,'(A)')LINE NCON=0 DO 2953 I=1,NCONST READ(3,'(I5,F30.19)',ERR=6851)K,A(I) NCON=NCON+1 2953 CONTINUE GOTO 6851 C C...Constants from a .CON file C 2951 DO 6856 I=1,NCON 2911 READ(3,6867,ERR=2910)BLUFF,A(I) 6867 FORMAT(A7,F30.19) GOTO 6856 2910 WRITE(*,'(1H+,A)')CHAR(7) GOTO 2911 6856 CONTINUE GOTO 6851 C C...Keyboard input of constants C 6854 WRITE(*,6849) 6849 FORMAT(/' Order of Hamiltonian (eg. 2, 4, 6 etc.) = ',$) READ(*,40,ERR=6854)IORDER IORDER=(IORDER/2)*2 NCON=0 DO 6850 I=IORDER,2,-2 NCON=NCON+I+1 6850 CONTINUE WRITE(*,4306) 4306 FORMAT(/' ROTATIONAL AND CENTRIFUGAL DISTORTION CONSTANTS /MHz', 1 ' :'/) DO 5000 I=1,NCONST WRITE(*,5001)ISTP(I) 2916 READ(*,41,ERR=2915)A(I) IF(I.GE.NCON)GOTO 6851 GOTO 5000 2915 WRITE(*,'(1H+,A1)')CHAR(7) GOTO 2916 5000 CONTINUE 5001 FORMAT(20X,A6,' ',$) C 6851 WRITE(*,7005) 7005 FORMAT(//' CALCULATION WILL BE CARRIED OUT WITH FOLLOWING ', 1 'CONSTANTS:'/) DO 6870 I=1,NCON IF(A(I).EQ.0.D0)GOTO 6870 WRITE(*,44)ISTP(I),A(I) WRITE(4,44)ISTP(I),A(I) 6870 CONTINUE C IF(IP.EQ.3)GOTO 5 WRITE(*,4307) 4307 FORMAT(/' DIPOLE MOMENT COMPONENTS /D :'/) DO 5002 I=1,3 2918 WRITE(*,5001)ISTP(I+NCONST) READ(*,41,ERR=2917)D(I) GOTO 5002 2917 WRITE(*,'(1H+,A1)')CHAR(7) GOTO 2918 5002 CONTINUE 2919 WRITE(*,5003) 5003 FORMAT(/' CUTOFF LINE STRENGTH = ',$) READ(*,41,ERR=2919)EPS 3900 WRITE(*,3901) 3901 FORMAT(' ROTATIONAL TEMPERATURE /K = ',$) READ(*,41,ERR=3900)TEMP 2920 WRITE(*,4308) 4308 FORMAT(' LOW FREQUENCY LIMIT /MHz = ',$) READ(*,41,ERR=2920)FMIN 2921 WRITE(*,4309) 4309 FORMAT(' HIGH FREQUENCY LIMIT /MHz = ',$) READ(*,41,ERR=2921)FMAX WRITE(4,50) WRITE(4,44)(ISTP(I+NCONST),D(I),I=1,3) C 5 JMIN=ABS(JMIN) JRMAX=MIN(ABS(JRMAX),200) JQMAX=MIN(ABS(JQMAX),200) JMAX=MAX(JRMAX,JQMAX) IF(KMMAX.EQ.0)KMMAX=JMAX K=1 IF(JMIN.GT.JMAX)GOTO 28 K=2 IF(A(1).LT.A(2).OR.A(2).LT.A(3))GOTO 28 IF(IP.NE.3)GOTO 9 WRITE(4,43)namred(iabs(ired)),namrep GOTO 14 9 EPS=DABS(EPS) FMIN=DABS(FMIN) FMAX=DABS(FMAX) K=4 IF(FMAX.LT.FMIN)GOTO 29 RJR=0.D0 RJQ=0.D0 DO 13 I=1,3 RJR=RJR*16.D0 RJQ=RJQ*4.D0 IF(D(4-I).EQ.0.D0)GOTO 13 RJR=RJR+15.D0 RJQ=RJQ+12288.D0 13 CONTINUE 14 M=JMAX/2+1 DO 15 I=1,M DO 15 K=1,M H(I,K)=0.D0 15 CONTINUE WRITE(*,2300) 2300 FORMAT(/' ***** CALCULATION IN PROGRESS, NOW WORKING ON J = '/) JQMIN=JMIN IF(JQMAX.GT.JMIN)GOTO 7006 JQMAX=0 JQMIN=0 7006 WRITE(4,46)JMIN,JRMAX,FMIN,FMAX,JQMIN,JQMAX,EPS,KMMAX IF(ISON.EQ.0)WRITE(4,4392) JJMIN=JMIN+1 JJMAX=JMAX+1 CLOSE(4) OPEN(4,FILE=FILNAM,ACCESS='APPEND') C C---------------------------------- C Main calculation loop - calculate energy levels for the specified C range of J C DO 24 JJ=JJMIN,JJMAX J=JJ-1 WRITE(*,2400)J 2400 FORMAT(I5,$) L=15 IF(IABS(IRED).NE.1)THEN CALL ASROHA(NCON) ELSE CALL ASROHS ENDIF IF(IP.LE.1)GOTO 19 IF(ISON.EQ.0)WRITE(4,50) DO 118 N=1,4 L=N-1 M=N/2*2 K=L/2+(4-N)/2*M L=J-K+M-L M=(J-K+2)/2 NN=(J-J/2*2)*4+N IF(M.LT.1)GOTO 118 DO 18 I=1,M WRITE(4,47)J,K,L,E(NN,I) K=K+2 L=L-2 18 CONTINUE 118 CONTINUE IF(IP.EQ.3)GOTO 24 19 RJS=1.D-6 IF(J.EQ.JMIN)GOTO 21 IF(J.GT.JRMAX)GOTO 21 RJS=RJS+RJR 21 IF(J.GT.JQMAX)GOTO 23 RJS=RJS+RJQ 23 CALL DIPTRA(DIDUTY,RJS,ISON,KMMAX) 24 CONTINUE C C---------------------------------- C C...IP = 1/2/3 calculate freq/freq+el/el IF(IP.EQ.3)GOTO 2 C C...IS = 0/1 do not sort/sort frequencies 25 IF(IS.NE.0)CALL SORT0(FV,SIGMA,GI,HALFW,ABUND,TEMP) GOTO 31 C C...Various error reports C 28 WRITE(*,48)K,JMIN,JMAX IF(IP.EQ.3)GOTO 2 GOTO 30 29 WRITE(*,49)K,EPS,FMIN,FMAX 30 IF(IS.NE.1.OR.NT.GT.1)GOTO 25 C 31 WRITE(*,50) C C 2 STOP END C C_____________________________________________________________________________ C SUBROUTINE ASROHA(NCON) C C ASYMMETRIC ROTOR ENERGY LEVELS FOR A GIVEN J - Reduction 'A', (IRED=+-2), C the sign of IRED determines whether representation I (+2) or III (-2) C (all four Wang matrices for a given J are calculated in one call of C ASROHA, the value of J is passed over in common block ASRO) C C The coding is set up to confer a slight speed advantage when constants C do not go beyond quartic. C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NCONST=35, MAXJ=200, MAXDIM=MAXJ/2+1) C COMMON /ASRO/I,J,L,NT,EA,EB,EE,EPS,FMIN,FMAX, 1 A(NCONST),S(3),E(8,MAXDIM),V(MAXDIM), 1 ID(MAXDIM),IP,MAP(8),IRED COMMON /BLKH/H(MAXDIM,MAXDIM)/BLKR/R COMMON /BIG/T C REAL*4 T(8,MAXDIM,MAXDIM),R(MAXDIM,MAXDIM) INTEGER*4 MM C C...various ancillary coefficients Q=J+1 Q=Q*J if(ired.eq.2)then C=(A(2)+A(3))*0.5D0 else C=(A(1)+A(2))*0.5D0 endif B=(-A(4)*Q+C)*Q if(ired.eq.2)then C=-A(5)*Q+A(1)-C D=(A(2)-A(3))*0.25D0-A(7)*Q else C=-A(5)*Q+A(3)-C D=(A(1)-A(2))*0.25D0-A(7)*Q endif F=-A(6) G=-A(8) IF(NCON.LE.8)GOTO 50 QQ=Q*Q QQQ=QQ*Q QQQQ=QQQ*Q QQQQQ=QQQQ*Q B=B+QQQ*A(9)+QQQQ*A(16)+QQQQQ*A(25) C=C+A(10)*QQ+A(17)*QQQ+A(26)*QQQQ F=F+A(11)*Q+A(18)*QQ+A(27)*QQQ FK6=A(12)+A(19)*Q+A(28)*QQ FK8=A(20)+A(29)*Q FK10=A(30) D=D+A(13)*QQ+A(21)*QQQ+A(31)*QQQQ G=G+A(14)*Q+A(22)*QQ+A(32)*QQQ FKK4=0.5D0*(A(15)+A(23)*Q+A(33)*QQ) FKK6=0.5D0*(A(24)+A(34)*Q) FKK8=0.5D0*A(35) C C...Process the four Wang matrices for a given J; the dimensionality of C each matrix is I C 50 DO 12 N=1,4 LN=L-2*(L/2) L=L/2 IF(LN.EQ.0)GOTO 12 K=(N-1)/2 M=N/2*2 KA=(4-N)/2*M+K KE=(J+K)/2*2-K IF(KE.LT.KA)GOTO 12 I=0 KA=KA+1 KE=KE+1 C C...go over K (the loop variable KK is K+1 to avoid problems for K=0) DO 7 KK=KA,KE,2 K=KK-1 I=I+1 MM=K MM=MM*K FLTA=MM IF(NCON.GT.8)GOTO 51 H(I,I)=(FLTA*F+C)*FLTA+B GOTO 52 51 H(I,I)=((((FK10*FLTA+FK8)*FLTA+FK6)*FLTA+F)*FLTA+C)*FLTA+B RMP2=(K+2)**2 SRMP2=RMP2*RMP2 SKSQ=FLTA*FLTA FLT4=SKSQ+SRMP2 FLT6=SKSQ*FLTA+SRMP2*RMP2 FLT8=SKSQ*SKSQ+SRMP2*SRMP2 52 IF(KK.GE.KE)GOTO 7 MM=MM+K+K+1 FLTA=MM+1 FLTB=-2*MM FLTC=MM FLTC=(FLTC-1.D0)*FLTC IF(NCON.GT.8)GOTO 60 H(I,I+1)=(D+FLTA*G)*SQRT((FLTB+Q)*Q+FLTC) GOTO 61 60 H(I,I+1)=(D+FLTA*G+FLT4*FKK4+FLT6*FKK6+FLT8*FKK8)* 1 SQRT((FLTB+Q)*Q+FLTC) 61 H(I+1,I)=H(I,I+1) 7 CONTINUE C C...touch up the top lh. corner elements IF(2-N)8,10,9 8 FLTA=-2*N+7 IF(NCON.LE.8)GOTO 55 H(1,1)=FLTA*(D+G+2.D0*(FKK4+FKK6+FKK8))*Q+H(1,1) GOTO 56 55 H(1,1)=FLTA*(D+G)*Q+H(1,1) 56 GOTO 10 9 H(1,2)=H(1,2)*1.414213562373095D0 H(2,1)=H(1,2) C C...diagonalise the Wang matrix, distribute eigenvalues and eigenvectors. C If representation III then the Wang matrices are loaded into E so that C the eigenvalues are in the same places as for representation I C 10 CALL HDIAG0 CALL HSORT(I) if(ired.eq.2)then LN=(J-J/2*2)*4+N else IMAP=(J-J/2*2)*4+N LN=(J-J/2*2)*4+MAP(IMAP) endif DO 111 K=1,I E(LN,K)=H(K,K) c WRITE(4,1111)H(K,K) c WRITE(4,1112)(R(K,M),M=1,I) c1111 FORMAT(' Eigenvalue = ',F12.3) c1112 FORMAT(' Eigenvector = ',12F10.3) DO 11 M=1,I T(LN,K,M)=R(K,M) H(K,M)=0.0D0 11 CONTINUE 111 CONTINUE 12 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE ASROHS C C ASYMMETRIC ROTOR ENERGY LEVELS FOR A GIVEN J - Reduction 'S', (IRED=+-1), C the sign of IRED determines whether representation I (+1) or III (-1) 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,3,4 for E+,E-,O+,O- resp. C A - array containing rotational and cd. constants in usual order C C Evaluated submatrix is placed in H, then diagonalised and eigenvalues C transferred to E, eigenvectors to T C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NCONST=35, MAXJ=200, MAXDIM=MAXJ/2+1) C COMMON /ASRO/I,J,L,NT,EA,EB,EE,EPS,FMIN,FMAX, 1 A(NCONST),S(3),E(8,MAXDIM),V(MAXDIM), 1 ID(MAXDIM),IP,MAP(8),IRED COMMON /BLKH/H(MAXDIM,MAXDIM)/BLKR/R COMMON /BIG/T C REAL*4 T(8,MAXDIM,MAXDIM),R(MAXDIM,MAXDIM) INTEGER*4 MM C C...ancillary coefficients Q=J+1 Q=Q*J QQ=Q*Q QQQ=QQ*Q QQQQ=QQQ*Q QQQQQ=QQQQ*Q C if(ired.eq.1)then C=0.5D0*(A(2)+A(3)) else C=0.5D0*(A(1)+A(2)) endif DK0=C*Q-A(4)*QQ+A(9)*QQQ+A(16)*QQQQ+A(25)*QQQQQ if(ired.eq.1)then DK2=A(1)-C-A(5)*Q+A(10)*QQ+A(17)*QQQ+A(26)*QQQQ else DK2=A(3)-C-A(5)*Q+A(10)*QQ+A(17)*QQQ+A(26)*QQQQ endif DK4=-A(6)+A(11)*Q+A(18)*QQ+A(27)*QQQ DK6=A(12)+A(19)*Q+A(28)*QQ DK8=A(20)+A(29)*Q DK10=A(30) C if(ired.eq.1)then OKK2=0.25D0*(A(2)-A(3))+A(7)*Q+A(13)*QQ+A(21)*QQQ+A(31)*QQQQ else OKK2=0.25D0*(A(1)-A(2))+A(7)*Q+A(13)*QQ+A(21)*QQQ+A(31)*QQQQ endif OKK4=A(8)+A(14)*Q+A(22)*QQ+A(32)*QQQ OKK6=A(15)+A(23)*Q+A(33)*QQ OKK8=A(24)+A(34)*Q OKK10=A(35) C C...Process the four Wang matrices for the specified J C C (KA and KE are the lowest and highest K in the submatrix: C KA = 0, 2, 1, 1 for E+, E-, O+, O- resp.) C DO 120 N=1,4 K=(N-1)/2 M=N/2*2 KA=(4-N)/2*M+K KE=(J+K)/2*2-K M=(KE-KA)/2+1 C DO 1 I=1,M DO 1 K=1,M H(I,K)=0.D0 1 CONTINUE H(1,1)=0.D0 H(1,2)=0.D0 C C For summing over K loop variable KK equal to K+1 is used to avoid C loop count from zero for E+ C Variable I is a pointer to elements in the Wang submatrix C I=0 KA1=KA+1 KE1=KE+1 DO 3 KK=KA1,KE1,2 K=KK-1 I=I+1 MM=K MM=MM*K FK2=MM c c... H(I,I)=((((DK10*FK2+DK8)*FK2+DK6)*FK2+DK4)*FK2+DK2)*FK2 * +DK0 c c... IF(K.GE.KE)GOTO 3 H(I,I+1)=OKK2*DSQRT(FPLUS(J,K,2)) H(I+1,I)=H(I,I+1) c c... IF(K.GE.KE-2)GOTO 3 H(I,I+2)=OKK4*DSQRT(FPLUS(J,K,4)) H(I+2,I)=H(I,I+2) c c... IF(K.GE.KE-4)GOTO 3 H(I,I+3)=OKK6*DSQRT(FPLUS(J,K,6)) H(I+3,I)=H(I,I+3) c c... IF(K.GE.KE-6)GOTO 3 H(I,I+4)=OKK8*DSQRT(FPLUS(J,K,8)) H(I+4,I)=H(I,I+4) c c... IF(K.GE.KE-8)GOTO 3 H(I,I+5)=OKK10*DSQRT(FPLUS(J,K,10)) H(I+5,I)=H(I,I+5) c 3 CONTINUE C C Final touch up: 4 -> O+ and O-, 12 -> E-, 5 -> E+ C IF(2-N)4,12,5 c c...O+ and O-, FLTA = +1 for O+, -1 for O- 4 FLTA=-2*N+7 H(1,1)=FLTA*OKK2*DSQRT(FMINUS(J,1,2))+H(1,1) IF(KE.LT.3)GOTO 6 H(1,2)=FLTA*OKK4*DSQRT(FMINUS(J,3,4))+H(1,2) H(2,1)=H(1,2) IF(KE.LT.5)GOTO 6 H(2,2)=FLTA*OKK6*DSQRT(FMINUS(J,3,6))+H(2,2) H(1,3)=FLTA*OKK6*DSQRT(FMINUS(J,5,6))+H(1,3) H(3,1)=H(1,3) IF(KE.LT.7)GOTO 6 H(1,4)=FLTA*OKK8*DSQRT(FMINUS(J,7,8))+H(1,4) H(4,1)=H(1,4) H(2,3)=FLTA*OKK8*DSQRT(FMINUS(J,5,8))+H(2,3) H(3,2)=H(2,3) IF(KE.LT.9)GOTO 6 H(1,5)=FLTA*OKK10*DSQRT(FMINUS(J,9,10))+H(1,5) H(5,1)=H(1,5) H(2,4)=FLTA*OKK10*DSQRT(FMINUS(J,7,10))+H(2,4) H(4,2)=H(2,4) H(3,3)=FLTA*OKK10*DSQRT(FMINUS(J,5,10))+H(3,3) GOTO 6 c c...root(2) elements of E+ 5 IF(KE.LT.2)GOTO 6 K=4 IF(KE.LT.6)K=1+KE/2 DO 10 KK=2,K H(1,KK)=H(1,KK)*1.414213562373095D0 H(KK,1)=H(1,KK) 10 CONTINUE c c...E- and remaining elements of E+, FLTA = +1 for E+ and -1 for E- 12 IF(KE.LT.2)GOTO 6 FLTA=-2*N+3 K=2 IF(N.EQ.2)K=1 H(K,K)= FLTA*OKK4*DSQRT(FMINUS(J,2,4))+H(K,K) IF(KE.LT.4)GOTO 6 H(K,K+1)=FLTA*OKK6*DSQRT(FMINUS(J,4,6))+H(K,K+1) H(K+1,K)=H(K,K+1) IF(KE.LT.6)GOTO 6 H(K,K+2)=FLTA*OKK8*DSQRT(FMINUS(J,6,8))+H(K,K+2) H(K+2,K)=H(K,K+2) H(K+1,K+1)=FLTA*OKK8*DSQRT(FMINUS(J,4,8))+H(K+1,K+1) IF(KE.LT.8)GOTO 6 H(K,K+3)=FLTA*OKK10*DSQRT(FMINUS(J,8,10))+H(K,K+3) H(K+3,K)=H(K,K+3) H(K+1,K+2)=FLTA*OKK10*DSQRT(FMINUS(J,6,10))+H(K+1,K+2) H(K+2,K+1)=H(K+1,K+2) C C...diagonalise the Wang matrix, distribute eigenvalues and eigenvectors C 6 CALL HDIAG0 CALL HSORT(I) if(ired.eq.1)then LN=(J-J/2*2)*4+N else IMAP=(J-J/2*2)*4+N LN=(J-J/2*2)*4+MAP(IMAP) endif DO 111 K=1,I E(LN,K)=H(K,K) DO 112 M=1,I T(LN,K,M)=R(K,M) H(K,M)=0.0D0 112 CONTINUE 111 CONTINUE 120 CONTINUE C RETURN END C C_____________________________________________________________________________ C FUNCTION FPLUS(J,K,KDELTA) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) INTEGER*4 KDELTA C C Values of f for matrix elements : C C dK-1 dK dK C + ___ ___ ___ C f (J,K) = | | (J-K-l) | | (J+K+l) = | | [J(J+1)-(K+l)(K+l-1)] C dK l=0 l=1 l=1 C FPLUS=1.0D0 RJ=J RK=K C DO 1 L=0,KDELTA-1 FPLUS=FPLUS*(RJ-RK-L) 1 CONTINUE DO 2 L=1,KDELTA FPLUS=FPLUS*(RJ+RK+L) 2 CONTINUE C RETURN END C C_____________________________________________________________________________ C FUNCTION FMINUS(J,K,KDELTA) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) INTEGER*4 K,KDELTA C C Values of f for matrix elements : C C dK-1 dK dK C - ___ ___ ___ C f (J,K) = | | (J+K-l) | | (J-K+l) = | | [J(J+1)-(K-l)(K-l+1)] C dK l=0 l=1 l=1 C FMINUS=1.0D0 RJ=J RK=K C DO 1 L=0,KDELTA-1 FMINUS=FMINUS*(RJ+RK-L) 1 CONTINUE DO 2 L=1,KDELTA FMINUS=FMINUS*(RJ-RK+L) 2 CONTINUE C RETURN END C_____________________________________________________________________________ C SUBROUTINE DIPTRA(DIDUTY,RJS,ISON,KMMAX) C C Calculate and determine which are acceptable transitions between levels C currently being worked on. Write these out to an intermediate file for C later sorting. C Space saving coding of quantum numbers from the original ASROT C is still retained (though it could now probably be dispensed with) C C Calculation of transitions is governed by numbers encoded into DIDUTY which C decode as follows: C ID( 3 4 5 6 7 8 9 10 11) XXX C C diduty( 1) =402657856. 0 0 1 1 1 0 0 0 0 3.0 IRQ=1 C diduty( 2) =421828864. 0 0 4 4 1 1 1 1 1 3.0 (delta J=+-1) C diduty( 3) =402953920. 0 0 3 3 1 1 1 0 0 3.0 C diduty( 4) =422122624. 0 0 2 2 1 2 2 1 1 3.0 C diduty( 5) =151267392. 0 0 1 4 2 0 1 0 1 1.0 C diduty( 6) =136652417. 1 0 2 3 2 2 1 1 0 1.0 C diduty( 7) =270574336. 0 0 4 1 2 1 0 1 0 2.0 C diduty( 8) =285779137. 1 0 3 2 2 1 2 0 1 2.0 C diduty( 9) =134493776. 0 2 1 3 3 0 1 0 0 1.0 C diduty(10) =153434257. 1 2 2 4 3 2 1 1 1 1.0 C diduty(11) =287880465. 1 2 4 2 3 1 2 1 1 2.0 C diduty(12) =268481232. 0 2 3 1 3 1 0 0 0 2.0 C --------------------------------------------------------------------- C diduty(13) =688395328. 0 0 1 2 1 0 2 0 1 5.0 IRQ=0 C diduty(14) =673486593. 1 0 4 3 1 1 1 1 0 5.0 (delta J=0) C diduty(15) =537142864. 0 2 1 3 2 0 1 0 0 4.0 C diduty(16) =556083345. 1 2 2 4 2 2 1 1 1 4.0 C diduty(17) =553924672. 0 0 1 4 3 0 1 0 1 4.0 C diduty(18) =539309697. 1 0 2 3 3 2 1 1 0 4.0 C | | | |__| C | | | used in calc. of no. of eigenv. C | | | in Wang matrices for J'' and J' C | | dipole component C | Wang submatrix for J' C Wang submatrix for J'' C C ID(10) and ID(11) are used in calculation of K+1 for energy levels for C J'' and J' resp. from K-1 and J, and define whether K-1 + K+1 = J (0) C or J+1 (1). C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NCONST=35, MAXJ=200, MAXDIM=MAXJ/2+1) C COMMON /ASRO/I,J,JS,NT,EA,EB,EE,EPS,FMIN,FMAX, 1 A(NCONST),D(3),E(8,MAXDIM),V(MAXDIM), 1 ID(MAXDIM),IP,MAP(8),IRED COMMON /BLKH/H(MAXDIM,MAXDIM)/BLKR/U COMMON /DEF/NREC COMMON /BIG/T C REAL*8 DIDUTY(18) REAL*4 T(8,MAXDIM,MAXDIM),U(MAXDIM,MAXDIM),EASP,ELOW C GF(K)=SIGN(1,L-N3+1)*SQRT(DBLE(REAL(((K+1)*N3+J)*(K*N3+J) 1 *(1+1/(K+1))))) NF(K)=2*K-(K+1)/2*4+1 C C...The variable J is J'. If this is even then the four sets of Wang c eigenvalues for this J are in E(1-4,..), those for J-1 being in c E(5-8,..). If J is odd then the eiegnvalues for this J are in E(5-8,..) c and for J-1 in E(1-4,..) c JP=0 IRQ=1 ID(1)=J c c...J even: I=+1 ID(12)=4 ID(13)=0 c J odd : -1 0 4 c I=(J/2*2-J)*2+1 ID(12)=(I+1)*2 ID(13)=(1-I)*2 DO 80 I=1,18 IF(I.NE.13)GOTO 3 JP=0 IRQ=0 3 K=MOD(RJS,2.D0) RJS=(RJS-K)*0.5D0 IF(K.EQ.0)GOTO 80 XXX=DIDUTY(I) DO 5 K=3,11 XXXX=MOD(XXX,8.D0) ID(K)=XXXX+1.D-6 XXX=(XXX-XXXX)*.125D0 5 CONTINUE IS=XXX+1.D-6 c c...I1,I2 are numbers of Wang matrices for J'' and J' resp. c M1,M2 are numbers of eigenvalues in Wang matrices E(I1,..) and E(I2,..) M1=(J+2-IRQ-ID(8))/2 IF(M1.EQ.0)GOTO 80 M2=(J+2-ID(9))/2 IF(M2.EQ.0)GOTO 80 I1=ID(5)+ID(13-IRQ) I2=ID(6)+ID(13) c c...L is either 1 or -1, N is either 0 or 1 c (the code for representation III is a quick hack to get the intensities c for c.R type transitions approximately right - other transitions types c NOT YET SUPPORTED for this representation) c ================= if(ired.lt.0)then if(id(7).eq.3.and.i.le.12)then id(20)=id(8) if(i.eq.9)id(20)=1 if(i.eq.12)id(20)=0 id(20)=id(20)+(j-(j/2)*2) endif endif L=ID(4)-1 N=ID(3) ID(4)=J-IRQ c c...SET UP THE DIRECTION COSINE MATRIX c if(ired.lt.0)then if(i.eq.1.or.i.eq.2)is=2 if(i.eq.3.or.i.eq.4)is=1 if(i.eq.6.or.i.eq.5)is=1 if(i.eq.7.or.i.eq.8)is=2 if(i.ge.9.and.i.le.12)is=3 if(i.eq.13.or.i.eq.14)is=4 if(i.eq.17.or.i.eq.18)is=5 endif GOTO(10,20,30,40,50),IS c c...IS=1, x.R and y.R transitions c 10 M=(J-1)/2*2 NN=N+1 MM=M+1 DO 11 KKK=NN,MM K=KKK-1 N1=(K-2*N+3)/2 N2=K/2+1 N3=NF(K) H(N1,N2)=GF(K) 11 CONTINUE GOTO 60 c c...IS=2, x.R and y.R transitions c 20 M=J/2*2-1 NN=N+1 MM=M+1 DO 21 KKK=NN,MM K=KKK-1 N1=(K-2*N+3)/2 N2=K/2+1 N3=-NF(K) H(N2,N1)=GF(K) 21 CONTINUE GOTO 60 c c...IS=3, z.R transitions c 30 DO 31 K=1,M1 if(ired.gt.0)then G=J**2-((K-1)*2+ID(8))**2 else G=J**2-((K-1)*2+ID(20))**2 endif H(K,K)=2.D0*SQRT(DABS(G)) 31 CONTINUE GOTO 60 c c...IS=4, x.Q and y.Q transitions c 40 M=J-1 NN=N+1 MM=M+1 DO 41 KKK=NN,MM K=KKK-1 N1=(K-2*N+3)/2 N2=K/2+1 N3=NF(K)+L+1 DBLA=((J+1)*J-(K+1)*K)*(1+1/(K+1)) H(N1,N2)=SIGN(1,N3)*SQRT(DBLA) 41 CONTINUE GOTO 60 c c...IS=5, z.Q transitions c 50 DO 51 K=1,M2 N1=K-N+1 H(N1,K)=((K-1)*2+ID(9))*2 51 CONTINUE c c...**2 c 60 G=((1-IRQ)*(2*J+1)/(J+1.D0)+IRQ)*0.25D0/J c c...loop over energy levels for J'' and J' c DO 169 K=1,M2 N1=0 DO 69 L=1,M1 EE=E(I2,K)-E(I1,L) EB=ABS(EE) IF(EB.LT.FMIN.OR.EB.GT.FMAX)GOTO 69 ELOW=E(I2,K) IF(E(I1,L).LT.E(I2,K))ELOW=E(I1,L) ELOW=ELOW*3.3356409E-5 IF(N1.NE.0)GOTO 64 N1=1 DO 63 M=1,M1 V(M)=0.D0 DO 63 N=1,M2 V(M)=T(I2,N,K)*H(M,N)+V(M) 63 CONTINUE 64 EA=0.D0 DO 65 M=1,M1 EA=T(I1,M,L)*V(M)+EA 65 CONTINUE N2=ID(7) EA=EA**2*G IF(EA.LT.EPS)GOTO 69 N2=(L-1)*2 N3=N2+ID(8) ID(5)=N3 IF(ID(5).GT.KMMAX)GOTO 69 ID(6)=ID(4)-N3+ID(10) N2=(K-1)*2 N3=N2+ID(9) ID(2)=N3 IF(ID(2).GT.KMMAX)GOTO 69 ID(3)=J-N3+ID(11) IF(JP.NE.0)GOTO 68 IF(ISON.EQ.0)WRITE(4,90) 90 FORMAT(/1X) C C...printout of quantum numbers (including dipole type), frequency and C line strength of identified, unsorted transitions 68 JP=1 IF(ISON.EQ.0)WRITE(4,91)(ID(KSQ),KSQ=1,7),EE,EA 91 FORMAT(1X,3I4,4X,3I4,I9,F21.6,F20.6) EASP=EA IF(NT.GE.8150)GOTO 269 C C...compression of quantum numbers into XXX and output of spectral line to C temporary file for later sorting with SORT0 C ID(7) - dipole component for transition C ID(1),ID(2),ID(3) - J,K-1,K+1 for upper state C ID(4),ID(5),ID(6) - J,K-1,K+1 for lower state C Five integer quantities are compressed: C Jupper,(Tau+J+1)upper,Jlower,(Tau+J+1)lower,dipole C These values are compressed (in mod.410) in reverse order so that on C unpacking they will come out as above C NT=NT+1 XXX=ID(7) IF(EE.GT.0.0d0)THEN ID(21)=ID(1) ID(22)=ID(2)-ID(3)+ID(1)+1 ID(23)=ID(4) ID(24)=ID(5)-ID(6)+ID(4)+1 ELSE ID(21)=ID(4) ID(22)=ID(5)-ID(6)+ID(4)+1 ID(23)=ID(1) ID(24)=ID(2)-ID(3)+ID(1)+1 ENDIF DO 13 KKK=24,21,-1 XXXX=ID(KKK) XXX=DNINT(XXX*410.D0+XXXX) 13 CONTINUE NREC=NT+1 WRITE(2,REC=NREC)XXX,EB,EASP,ELOW 69 CONTINUE WRITE(2,REC=1)NT 169 CONTINUE 269 WRITE(2,REC=1)NT c c...clear direction cosine matrix for next transition type H(1,1)=0.D0 N3=J/2+1 DO 70 K=2,N3 H(K,K)=0.D0 H(K,K-1)=0.D0 H(K-1,K)=0.D0 70 CONTINUE 80 CONTINUE C C RETURN END C C_____________________________________________________________________________ C SUBROUTINE SORT0(FV,SIGMA,GI,HALFW,ABUND,TEMP) C C Read and decode transitions stored in an intermediate file and sort C these in ascending order of frequency. C C NQ(I,1),NQ(I,2),NQ(I,3) - J,K-1,K+1 of the upper state C NQ(I,4),NQ(I,5),NQ(I,6) - J,K-1,K+1 of the lower state C NQ(I,7) - type of transition flag = 1,2,3 for a,b,c C NQ(I,8) - Delta K-1 C NQ(I,9) - Delta K+1 C NQ(I,10) - (Delta J) + 2 C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NCONST=35, MAXJ=200, MAXDIM=MAXJ/2+1) C COMMON /ASRO/B,C,LL,N,EA,EB,EE,EPS,FMIN,FMAX, 1 A(NCONST),D(3),E(8,MAXDIM),V(MAXDIM), 1 ID(MAXDIM),IP,MAP(8),IRED COMMON /BLKH/S/BLKR/M COMMON /DEF/NREC COMMON /BIG/T C REAL*8 L,JP(8150),F(8150),T(40804) INTEGER*2 B,C,NQ(8150,10),M(20402) INTEGER*4 ITRANS(8150) CHARACTER DIPOLE(3),TYPETR(3) C C...To save space JP, F and EL are equivalenced with T. Other large C arrays NQ, ITRANS and S are not equivalenced although S is read C into the eigenmatrix array in the /ASRO/ block. The eigenvector C array in /ASRO/ (here named M) is used to hold pointers for sorting C REAL*4 S(20402),EL(8150),ALFA,FV,SIGMA,GI,HALFW,ABUND,TEMP, 1 CONST,CONST1,CONST3,FREQ EQUIVALENCE (T(1),JP(1)),(T(8151),F(1)),(T(16301),EL(1)), * (S(8151),ITRANS(1)),(T(20376),NQ(1,1)) DATA DIPOLE/'a','b','c'/TYPETR/'P','Q','R'/ FRESOL=0.001D0 C 20 FORMAT(/1x,78(1H-)/' LINES IN ORDER OF ASCENDING FREQUENCY:', 1 10X,I5,' lines in temp. file' 1 //' FREQUENCY', 1 4X,'Alpha max',4X,'J K K J K K',4X,'Branch',4X, 1 'Line',8X,'Lower'/34X,'+1 -1 +1 -1',12X,'Strength', 1 4X,'Level'/4X,'/MHz',9X,'/cm-1',52X,'cm-1'/) 2500 FORMAT(//' *** LINES ARE BEING SORTED') 2501 FORMAT(1X) C WRITE(*,2500) C C Read unsorted frequencies: JP contains packed quantum numbers, C F - frequencies, S - line strengths, EL - energy of the lower level C READ(2,REC=1)N DO 200 I=1,N NREC=I+1 READ(2,REC=NREC)JP(I),F(I),S(I),EL(I) 200 CONTINUE C C Sorting loop - numbers of transitions in ascending order of frequency C are placed in M, F is sorted, JP,S and EL remain unsorted C WRITE(4,20)N DO 1 I=1,N M(I)=I 1 CONTINUE N1=N-1 DO 4 I=1,N1 EE=F(I) J=I LX=I+1 DO 3 K=LX,N IF(F(K).GE.EE)GOTO 3 EE=F(K) J=K 3 CONTINUE F(J)=F(I) F(I)=EE K=M(J) M(J)=M(I) M(I)=K 4 CONTINUE C C Unpack quantum numbers and generate transition type index ITRANS C - this is used for identification of unresolved K+1 or K-1 doublets C (note that Tau+J+1 is compressed) C K+1 = (J + 1 - Tau) /2 , K-1 = Tau + K+1 C DO 10 I=1,N J=M(I) L=JP(J) DO 5 K=1,5 EE=DMOD(L,410.D0) NQ(I,K)=NINT(EE) L=DNINT((L-EE)/410.D0) 5 CONTINUE NQ(I,7)=NQ(I,5) NQ(I,4)=NQ(I,4)-NQ(I,3)-1 NQ(I,2)=NQ(I,2)-NQ(I,1)-1 NQ(I,6)=(NQ(I,3)+1-NQ(I,4))/2 NQ(I,5)=NQ(I,4)+NQ(I,6) NQ(I,4)=NQ(I,3) NQ(I,3)=(NQ(I,1)+1-NQ(I,2))/2 NQ(I,2)=NQ(I,2)+NQ(I,3) NQ(I,8)=NQ(I,2)-NQ(I,5) NQ(I,9)=NQ(I,3)-NQ(I,6) NQ(I,10)=NQ(I,1)-NQ(I,4)+2 ITRANS(I)=(NQ(I,7)*10+NQ(I,10))*10 ITRANS(I)=(ITRANS(I)+(NQ(I,8)+50))*100+NQ(I,9)+50 10 CONTINUE C C Calculate intensities, establish types of transitions and C write out sorted transitions C CONST=3.85E-14*DSQRT(A(1)*A(2)*A(3))*FV*GI*SIGMA/(HALFW 1 *TEMP**2.5)*ABUND CONST1=0.5/(29979*0.69497) IOMIT=0 NOUTP=1 C C DO 6 I=1,N IF(IOMIT.EQ.1)GOTO 11 J=M(I) L=JP(J) C C...absorption coefficient CONST3=-EL(J)/(0.69497*TEMP) IF(CONST3.LT.-75.)GOTO 6 FREQ=F(I) ALFA=CONST*SNGL(D(NQ(I,7)))**2*S(J)*FREQ**2*EXP(CONST3) 1 *(1.-CONST1*FREQ/TEMP) IF(ALFA.LT.1.E-12)GOTO 6 C C...check for unresolved doublets IF(I.EQ.N)GOTO 12 IF(DABS(F(I)-F(I+1)).GT.FRESOL)GOTO 12 IF(ITRANS(I).NE.ITRANS(I+1))GOTO 12 IF((NQ(I,1).NE.NQ(I+1,1)).OR.(NQ(I,4).NE.NQ(I+1,4)))GOTO 12 IF((NQ(I,2).EQ.NQ(I+1,2)).AND.(NQ(I,5).EQ.NQ(I+1,5)))GOTO 14 IF((NQ(I,3).EQ.NQ(I+1,3)).AND.(NQ(I,6).EQ.NQ(I+1,6)))GOTO 15 GOTO 12 C C...doublet with K-1's the same 14 ALFA=ALFA*2. F(I)=0.5D0*(F(I)+F(I+1)) IOMIT=1 WRITE(4,17)F(I),ALFA,NQ(I,1),NQ(I,2),NQ(I,4),NQ(I,5), 1 DIPOLE(NQ(I,7)), 1 TYPETR(NQ(I,10)),NQ(I,8),NQ(I,9),S(J),EL(J) 17 FORMAT(F14.6,1PE10.2,' D',I3,',',I3,',',I8,',',I3,',', 1 4X,A1,'.',A1,I2,',',I2,0PF10.6,F10.3) IF((NOUTP/10)*10.EQ.NOUTP)WRITE(4,2501) NOUTP=NOUTP+1 GOTO 6 C C...doublet with K+1's the same 15 ALFA=ALFA*2. F(I)=0.5D0*(F(I)+F(I+1)) IOMIT=1 WRITE(4,16)F(I),ALFA,NQ(I,1),NQ(I,3),NQ(I,4),NQ(I,6), 1 DIPOLE(NQ(I,7)), 1 TYPETR(NQ(I,10)),NQ(I,8),NQ(I,9),S(J),EL(J) 16 FORMAT(F14.6,1PE10.2,' D',I3,', ,',I3,I5,', ,',I3, 1 1X,A1,'.',A1,I2,',',I2,0PF10.6,F10.3) IF((NOUTP/10)*10.EQ.NOUTP)WRITE(4,2501) NOUTP=NOUTP+1 GOTO 6 C C...single transition 12 WRITE(4,21)F(I),ALFA,(NQ(I,K),K=1,6),DIPOLE(NQ(I,7)), 1 TYPETR(NQ(I,10)),NQ(I,8),NQ(I,9),S(J),EL(J) 21 FORMAT(F14.6,1PE10.2,2X,2(I3,',')I3,I5,2(',',I3), 1 1X,A1,'.',A1,I2,',',I2,0PF10.6,F10.3) IF((NOUTP/10)*10.EQ.NOUTP)WRITE(4,2501) NOUTP=NOUTP+1 GOTO 6 C 11 IOMIT=0 C 6 CONTINUE C C C Concluding information C WRITE(4,666)NOUTP,FV,SIGMA,GI,HALFW,ABUND,TEMP 666 FORMAT(1X/49X,I5,' lines printed' 1 //' ABSORPTION COEFFICIENTS WERE CALCULATED WITH:'// 1 7X,' Vibrational population factor, F.v = ',F9.4/ 1 19X,' Symmetry number, sigma = ',F9.4/ 1 2X,' Reduced nuclear statistical weight, g.I = ',F9.4/ 1 ' Pressure broadening parameter /MHz Torr-1 = ',F9.4/ 1 11X,' Isotopic abundance factor, i.c = ',F9.4/ 1 25X,' Temperature, T/K = ',F9.4/1x,78(1H-)) C RETURN END C C_____________________________________________________________________________ C SUBROUTINE HDIAG0 C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NCONST=35, MAXJ=200, MAXDIM=MAXJ/2+1) C COMMON /ASRO/N,M,L,NT,EA,EB,EE,SPE,FMIN,FMAX, 1 A(NCONST),D(3),E(8,MAXDIM),X(MAXDIM), 1 IQ(MAXDIM),IP,MAP(8),IRED COMMON /BLKH/H(MAXDIM,MAXDIM)/BLKR/U COMMON /BIG/T C REAL*4 T(8,MAXDIM,MAXDIM),U(MAXDIM,MAXDIM),RSIN,RCOS,RHTEMP C IF(IP.EQ.3)GOTO 15 DO 14 I=1,N DO 14 J=1,N IF(I-J)12,11,12 11 U(I,J)=1. GOTO 14 12 U(I,J)=0. 14 CONTINUE 15 NR=0 IF(N-1)1000,1000,17 17 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)))20,20,30 20 X(I)=ABS(H(I,J)) IQ(I)=J 30 CONTINUE 40 DO 70 I=1,NMI1 IF(I-1)60,60,45 45 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,120 120 NR=NR+1 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 SUNO=-SINE HTEMP=SIGN(1.D0,SUNO)*COSINE COSINE=ABS(SINE) SINE=HTEMP 153 DO 350 I=1,NMI1 IF(I-IPIV)210,350,200 200 IF(I-JPIV)210,350,210 210 IF(IQ(I)-IPIV)230,240,230 230 IF(IQ(I)-JPIV)350,240,350 240 K=IQ(I) 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)))300,300,320 300 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)))380,390,390 380 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)))400,530,530 400 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)))440,450,450 440 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)))460,530,530 460 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)))490,500,500 490 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)))510,530,530 510 X(JPIV)=ABS(H(JPIV,I)) IQ(JPIV)=I 530 CONTINUE IF(IP.EQ.3)GOTO 40 RSIN=SINE RCOS=COSINE DO 550 I=1,N RHTEMP=U(I,IPIV) U(I,IPIV)=RCOS*RHTEMP+RSIN*U(I,JPIV) U(I,JPIV)=-RSIN*RHTEMP+RCOS*U(I,JPIV) 550 CONTINUE GOTO 40 1000 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE HSORT(M) C C SORTING OF EIGENVALUES AND EIGENVECTORS COMING OUT OF HDIAG C - introduced after it was found that for representation III some C HDIAG eigenvalues did not keep ascending order. C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (MAXJ=200, MAXDIM=MAXJ/2+1) C COMMON /BLKH/A/BLKR/U/SORTCC/WK,IPT REAL*8 A(MAXDIM,MAXDIM),WK(MAXDIM),D(MAXDIM) REAL*4 U(MAXDIM,MAXDIM) INTEGER*2 IPT(MAXDIM) isort=1 C C...Check whether additional sorting necessary C DO 150 I=1,M-1 IF(A(I+1,I+1).LT.A(I,I))GOTO 151 150 CONTINUE RETURN C C...Transfer eigenvalues from the diagonal of matrix A to D, and eigenvectors C from U to A C 151 DO 115 I=1,M D(I)=A(I,I) DO 115 J=1,M A(I,J)=U(I,J) 115 CONTINUE C C...Straightforward magnitude sorting C DO 122 I=1,M IPT(I)=I WK(I)=D(I) 122 CONTINUE CALL SORTC(isort,M) C C...Return of sorted eigenvectors in U and eigenvalues on the diagonal of A C 140 DO 117 II=1,M I=IPT(II) DO 117 J=1,M U(J,II)=A(J,I) 117 CONTINUE DO 118 II=1,M I=IPT(II) A(II,II)=D(I) 118 CONTINUE C RETURN END C C____________________________________________________________________________ C SUBROUTINE SORTC(N,M) IMPLICIT INTEGER*2 (I-N) PARAMETER (JMAX=200,MAXDIM=JMAX/2+1) C COMMON /SORTCC/WK,IPT INTEGER*2 IPT(MAXDIM) REAL*8 WK(MAXDIM),EE C C ... This routine sorts the quantities in vector WK in ascending order C of magnitude and also accordingly rearranges vector IPT of pointers C 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_____________________________________________________________________________