C$DEBUG C$LARGE C------------------------------------------------------------------------- C C PROGRAM STARK C C MODIFIED RIBEAUD'S PROGRAM FOR THE CALCULATION OF THE STARK EFFECT C FOR AN ASYMMETRIC ROTATOR. C C (This modification for PC operation is based on copy of STARK C from group of A.P.Cox) C C ver 4/1992 ----- ZBIGNIEW KISIEL ----- C C------------------------------------------------------------------------- C C STARK reads data from file STARK.DAT and outputs to STARK.RES. C Structure of STARK.DAT (repeatable): C C LINE.1: PROG,OB,JMIN,JMAX,(ITXT(I),I=1,18) (4I2,18A4) C LINE.2: AA(1),AA(2),AA(3) (5F15.0) C LINE.3: AA(4),AA(5),AA(6),AA(7),AA(8) (5F15.0) C LINE.4: EFELD,EPS,EP,FMIN,FMAX (5F15.0) C LINE.5: (DM(I),I=1,3),(DS(I),I=1,3) (6F10.0) C C where: C C PROG - program control parameter to be set to 1,2 or 3 C 1 = only frequencies C 2 = frequencies and energy levels C 3 = only energy levels C OB - to be set to -1,0 or +1 (delta M? - in any case standard DeltaM=0 C results are obtained with OB=0) C JMIN - minimum J for the calculation C JMAX - maximum J for the calculation, up to 60 C ITXT - alphanumeric comment, up to 72 characters C AA - rotational constants A,B,C and quartic cd. constants C EFELD - electric field (Volts/cm) C EPS - threshold criterion for modified calculation for degeneracies? C EP - intensity cutoff for transitions C FMIN - lower frequency cutoff for calculated transitions C FMAX - upper frequency cutoff for calculated transitions C DM - effective dipole moment components (Debyes) C DS - stark? dipole moment components (Debyes, to be set equal to DM - C a difference between DM and DS presumably allows to account C for inducerd dipoles) C C SPECIMEN DATA SET FOR THE SO2 EXAMPLE IN GORDY AND COOK (pure second C order effect, mu.b only, frequencies only, cd. neglected): C C---- 1 0 0 3 SO2 - example from Gordy and Cook C---- 60778.79 10318.10 8799.96 C---- 0.0 0.0 0.0 0.0 0.0 C---- 2005.1 0.001 0.1000 50000.0 80000. C---- 0.0 1.58 0.0 0.0 1.58 0.0 C C ALL CONSTANTS WITH UNITS OF FREQUENCY ARE TO BE EXPRESSED IN MHz. C Watson's A-type centrifugal distortion constants and Wang C factorisation are employed. C C------------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION H(32,32),E(12,32),A(7,12,32),D(36,32,32),T(8,32,32), 1F(32,32),DM(3),DS(3),DSS(3),AA(8),ITXT(18),KK1(18),KK2(18), 2LL1(18),LL2(18),LL4(18),MM1(18),MM2(18) INTEGER B(12,34),V,VV,OB,PROG,SS1(18),SS2(18),SN,SN1,SN2 C DATA KK1/0,1,0,2,0,2,0,1,1,2,0,2,1,1,0,2,1,1/ DATA KK2/2,1,1,1,1,1,0,1,1,2,1,1,0,2,1,1,2,0/ DATA LL1/0,1,0,1,0,1,0,1,0,1,0,1,1,0,0,1,1,0/ DATA LL2/1,0,0,1,1,0,0,1,0,1,1,0,0,1,0,1,1,0/ DATA LL4/1,1,2,2,3,3,1,1,1,1,2,2,2,2,3,3,3,3/ DATA SS1/1,4,1,2,1,2,5,8,7,6,5,6,8,7,5,6,8,7/ DATA SS2/2,3,3,4,4,3,1,4,3,2,4,3,1,2,3,4,2,1/ DATA MM1/2,1,2,0,2,0,1,0,0,-1,1,-1,0,0,1,-1,0,0/ DATA MM2/0,1,1,1,1,1,2,1,1,0,1,1,2,0,1,1,0,2/ DATA MSA/64/,MSB/63/ C WRITE(*,121) WRITE(*,200) 200 FORMAT(1X,'STARK - PREDICTION OF STARK FREQUENCY SHIFTS F', 1 'OR AN ASYMMETRIC ROTATOR'/1X,70(1H=)//) open(1,file='stark.dat') open(2,file='stark.res',status='unknown',access='append') 80 READ(1,110,ERR=90,END=90) PROG,OB,JMIN,JMAX,(ITXT(I),I=1,18) IF (PROG.LT.1.OR.PROG.GT.3) GO TO 90 READ(1,111) (AA(I),I=1,3) READ(1,111) (AA(I),I=4,8) READ(1,111) EFELD,EPS,EP,FMIN,FMAX READ(1,112) (DM(I),I=1,3),(DS(I),I=1,3) IF(JMIN.LT.0.OR.JMAX.GT.61) GO TO 70 FM=FMAX-FMIN C IF (PROG.NE.3.AND.FM.LE.1.0D-8) GO TO 70 UTBYTT IF (PROG.NE.3.AND.FM.LE.1.0D-16)GO TO 70 NYTT IG=(JMAX+2)/2 DO 6 I=1,IG DO 5 K=1,IG H(I,K)=0.0D0 DO 4 IK=1,8 F(I,K)=H(I,K) T(IK,I,K)=0.0D0 4 CONTINUE DO 5 KI=1,36 D(KI,I,K)=0.0D0 5 CONTINUE DO 6 L=1,12 E(L,I)=0.0D0 B(L,I)=0 DO 6 LI=1,7 A(LI,L,I)=0.0D0 6 CONTINUE IF(EFELD.EQ.0) GO TO 12 GO TO 13 12 EFELD=100.0D0 C C...note the conversion factor 0.503447 MHz (Debye Volt/cm)-1 for use C of customary units C 13 C3=EFELD*0.503447D0 NYTT C 13 C3=EFELD/1.98768D0 UTBYTT C IF (DABS(EPS).LT.1.0D-8) GO TO 14 UTBYTT IF (DABS(EPS).LT.1.0D-16)GO TO 14 NYTT GO TO 15 14 EPS=10.0D0 15 DO 19 L=1,3 C IF (DABS(DS(L)).LT.1.0D-8) GO TO 17 UTBYTT IF (DABS(DS(L)).LT.1.0D-16)GO TO 17 NYTT GO TO 18 17 DSS(L)=0.0D0 GO TO 19 18 DSS(L)=(DM(L)/DS(L))**2 19 CONTINUE WRITE(2,100) (ITXT(I),I=1,18) WRITE(2,106) WRITE(2,108) WRITE(2,131) 'A',AA(1),DM(1),DS(1) WRITE(2,131) 'B',AA(2),DM(2),DS(2) WRITE(2,131) 'C',AA(3),DM(3),DS(3) WRITE(2,107) WRITE(2,120) (AA(LU),LU=4,8) WRITE(2,102) JMIN,JMAX WRITE(2,105) EFELD NULL=1 GO TO (21,21,20), PROG 20 NULL=0 GO TO 40 21 WRITE(2,103) FMIN,FMAX WRITE(2,104) EP,EPS,OB WRITE(2,109) C C...MAIN LOOP OVER J C WRITE(*,'('' Now working on J =''/)') 40 DO 60 J=JMIN,JMAX WRITE(*,'(I4,$)')J WRITE(2,99) JOB=15 CALL ASROHA(J,JOB,AA,E,H,F,T) IF (J.EQ.JMIN) GO TO 44 GO TO 45 44 MIN=1 IF (J.EQ.0) GO TO 60 45 V=1 VV=87108 CALL DIPOL(J,MIN,V,VV,NULL,MSA,MSB,C3,EPS,DSS,DS, * LL1,LL2,LL4,KK1,KK2,MM1,MM2,SS1,SS2,E,D,T,A,B,OB) MIN=0 JV=J-1 GO TO (47,46,46), PROG 46 CALL EDRUCK(JV,JMIN,E,A,DM,DSS) GO TO (47,47,55), PROG 47 IF (J.GT.JMIN) GO TO 48 GO TO 60 48 IF (JV.EQ.JMIN.AND.JMIN.NE.0) GO TO 49 GO TO 50 49 MIN=1 50 IF (JV.EQ.0) GO TO 60 CALL UBERG(JV,MIN,OB,EP,FMIN,FMAX,E,D,A,B,DM,DSS,LL1,LL2,LL4, * KK1,KK2,MM1,MM2,SS1,SS2) 55 JV=J-2 LIM=(JV+2)/2 SN1=(JV-JV/3*3)*4+1 SN2=SN1+3 DO 56 SN=SN1,SN2 DO 56 L=1,LIM E(SN,L)=0.0D0 B(SN,L)=0 DO 56 I=1,7 A(I,SN,L)=0.0D0 56 CONTINUE MIN=0 60 CONTINUE C 99 FORMAT(1X) 100 FORMAT(1x,72(1H-)/1x,18A4/1x,72(1H-)/) 101 FORMAT(4X,F15.7,3X,2F10.5,5F12.7) 131 FORMAT(A1,2X,F16.7,3X,2F10.5,5F12.7) 120 FORMAT(1X,5F15.8) 102 FORMAT(/' J LIMITS:', 1 I4,' to',I4) 103 FORMAT(' FREQUENCY LIMITS:',F12.2,' to',F12.2, 1 ' MHz') 104 FORMAT(' MINIMUM INTENSITY:',F12.8/ 1 ' EXPANSION IN TERMS OF EPS:',F8.4/ 1 ' OB:',I4) 105 FORMAT(' ELECTRIC FIELD:',F11.5,' Volts/cm') 106 FORMAT(' --- ROTATION AND CD. CONSTANTS MUST BE', * ' GIVEN IN MHz,'/' ALL RESULTS ARE IN MHz') 107 FORMAT(/' QUARTIC CENTRIFUGAL DISTORTION CONSTANTS:'/) 108 FORMAT(/5X,'ROTATION CONST',4X,'EFF.DIPOLE',2X,'ST.DIPOLE'/) 109 FORMAT(/' For second order Stark effect: dF = A + M**2 B'// 1 1x,72(1H-)/ 1 ' J K K <- J K K Dipole Frequency Strength'/ 1 ' -1 +1 -1 +1 ',35X,'A B') 110 FORMAT(4I2,18A4) 111 FORMAT(5F15.0) 112 FORMAT(6F10.0) 121 FORMAT(/1X,70(1H=)) 70 GO TO 80 90 CONTINUE WRITE(2,121) CLOSE(1) CLOSE(2) WRITE(*,121) C STOP END C C------------------------------------------------------------------------- C SUBROUTINE DIPOL(J,MIN,V,VV,NULL,MSA,MSB,C3,EPS,DSS, 1DS,LL1,LL2,LL4,KK1,KK2,MM1,MM2,SS1,SS2,E,D,T,A,B,OB) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DS(3),E(12,32),D(36,32,32),T(8,32,32),A(7,12,32), 1F(32,32),DSS(3) INTEGER S1,S2,S3,V,VV,ST1,ST2,OB,B(12,34) INTEGER LL1(18),LL2(18),LL4(18),KK1(18),KK2(18),MM1(18), 1MM2(18),SS1(18),SS2(18) C NNUM=1+(J-J/2*2)*18 NS=(J-J/2*2)*4 NSD=NS IF (MIN.EQ.0) GO TO 5 MAX=NNUM+5 GO TO 6 5 MAX=NNUM+17 6 DO 250 L=NNUM,MAX LLL=L/19 LL=L-LLL*18 L4=LL4(LL) IF (LL.GT.6) GO TO 10 S3=0 GO TO 11 10 S3=1 NSD=-NS 11 K1=KK1(LL) K2=KK2(LL) L1=LL1(LL) L2=LL2(LL) K3=(K1*K2+1)/2 L3=(-1)**(L1-L2) M1=(J+MM1(LL))/2 M2=(J+MM2(LL))/2 JM=J-S3 ST1=SS1(LL)+NSD S1=ST1+(JM/2*2-JM/3*3)*4 ST2=SS2(LL)+NS S2=ST2+(J/2*2-J/3*3)*4 IZ=VV-(VV/2)*2 V=V+IZ VV=VV/2 C IF (DABS(DS(L4)).LE.1.0D-8) GO TO 250 UTBYTT IF (DABS(DS(L4)).LE.1.0D-16)GO TO 250 NYTT GO TO (50,40,160,140,150,140,150),V 40 J1=J-1 IF (K3.EQ.0) GO TO 41 GO TO 42 41 D(L,1,1)=DSQRT(2.0D0*J*(J+1)) 42 IF (J1.GT.0) GO TO 43 GO TO 45 43 DO 44 I=1,J1 I1=(I-2*K3+3)/2 I2=I/2+1 I3=2*I-(I+1)/2*4+L3+2 I3=ISIGN(1,I3) G=J*(J+1)-I*(I+1) 44 D(L,I1,I2)=I3*DSQRT(G) 45 GO TO 199 50 J1=(J-K2)/2 IF (J1.GE.0) GO TO 51 GO TO 55 51 JJ1=J1+1 DO 52 IPI=1,JJ1 I=IPI-1 I1=I-K3+2 I2=I+1 52 D(L,I1,I2)=4*I+2*K2 55 GO TO 199 140 IF (K3.EQ.0) GO TO 141 GO TO 142 141 D(L,1,1)=L3*DSQRT(2.0D0*J*(J+1)) 142 JJ=(J-1)/2*2 IF (JJ.GE.1) GO TO 143 GO TO 145 143 DO 144 I=1,JJ I1=(I-2*K3+3)/2 I2=I/2+1 I3=(I+1)/2 I3=2*I-4*I3+1 I4=L3-I3+1 I4=ISIGN(1,I4) G=(J+I*I3)*(J+(I+1)*I3) 144 D(L,I1,I2)=I4*DSQRT(G) 145 GO TO 199 150 IF (K3.EQ.0) GO TO 151 GO TO 152 151 D(L,1,1)=DSQRT(2.0D0*J*(J-1)) 152 JJ=J/2*2-1 IF (JJ.GE.1) GO TO 153 GO TO 155 153 DO 154 I=1,JJ I1=I/2+1 I2=(I-2*K3+3)/2 I3=(I+1)/2 I3=4*I3-2*I-1 I4=L3-I3+1 I4=ISIGN(1,I4) G=(J+I*I3)*(J+(I+1)*I3) 154 D(L,I1,I2)=I4*DSQRT(G) 155 GO TO 199 160 JJ=(J-K1-1)/2 IF (JJ.GE.0) GO TO 161 GO TO 199 161 JJ11=JJ+1 DO 162 IPI=1,JJ11 I=IPI-1 I1=I+1 G=J**2-(K1+2*I)**2 162 D(L,I1,I1)=2*DSQRT(G) 199 IF (M1.GE.1.AND.M2.GE.1) GO TO 200 GO TO 250 200 DO 201 I=1,M2 DO 201 K=1,M1 F(I,K)=0.0D0 DO 201 LI=1,M1 201 F(I,K)=F(I,K)+T(ST1,LI,K)*D(L,LI,I) DO 202 I=1,M1 DO 202 K=1,M2 D(L,I,K)=0.0D0 DO 202 LI=1,M2 202 D(L,I,K)=D(L,I,K)+F(LI,I)*T(ST2,LI,K) CALL ENTART(J,JM,L,NULL,K1,K2,L1,L2,L4,S1,S2,S3,M1,M2,MSA,MSB, 1C3,EPS,DS,E,D,A,B,F,DSS,OB) IF (LL.GT.6) GO TO 220 CALL Q2STRK(J,F,L4,M1,M2,S1,S2,A) GO TO 250 220 CALL R2STRK(J,F,L4,M1,M2,S1,S2,A) 250 CONTINUE C RETURN END C C------------------------------------------------------------------------- C SUBROUTINE EDRUCK(J,JMIN,E,A,DM,DSS) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION E(12,32),A(7,12,32),DM(3),DSS(3),ST(9) C WRITE(2,102) J DO 21 N=1,4 NS=(J-J/3*3)*4+N L=N-1 M=N/2*2 K=L/2+(4-N)/2*M L=J-K+M-L M=(J-K+2)/2 IF (M.GE.1) GO TO 11 GO TO 21 11 DO 20 I=1,M DO 12 II=1,7 12 ST(II)=A(II,NS,I) ST(8)=ST(1)*DSS(1)+ST(3)*DSS(2)+ST(5)*DSS(3) ST(9)=ST(2)*DSS(1)+ST(4)*DSS(2)+ST(6)*DSS(3) WRITE(2,100) J,K,L,E(NS,I),(ST(LA),LA=7,9) WRITE(2,101) (ST(LA),LA=1,6) K=K+2 L=L-2 20 CONTINUE 21 CONTINUE C 100 FORMAT(2X,3I3,8X,F25.10,34X,3F12.7) 101 FORMAT(46X,6F12.7) 102 FORMAT(' ENERGY LEVELS J =',I3/) C RETURN END C C------------------------------------------------------------------------- C SUBROUTINE ASROHA(J,L,AA,E,H,R,T) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(8),E(12,32),R(32,32),H(32,32),T(8,32,32) C Q=(J+1)*J C=(AA(2)+AA(3))/2.0D0 B=(-AA(4)*Q+C)*Q C=-AA(5)*Q+AA(1)-C D=(AA(2)-AA(3))/4.0-AA(7)*Q F=-AA(6) G=-AA(8) DO 12 N=1,4 LN=L-2*(L/2) L=L/2 IF (LN.EQ.0) GO TO 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) GO TO 12 I=0 DO 7 K=KA,KE,2 I=I+1 M=K*K H(I,I)=(M*F+C)*M+B IF (K.LT.KE) GO TO 6 GO TO 7 6 M=M+K+K+1 H(I,I+1)=(((K+2)*K+2)*G+D)*DSQRT((-2*M+Q)*Q+(K+2)*K*M) H(I+1,I)=H(I,I+1) 7 CONTINUE IF (2-N) 8,10,9 8 H(1,1)=(-2*N+7)*(G+D)*Q+H(1,1) GO TO 10 9 H(1,2)=H(1,2)*DSQRT(2.0D0) H(2,1)=H(1,2) 10 CALL HDIAG(H,I,0,R,LN) LN=(J-J/3*3)*4+N LNC=(J-J/2*2)*4+N DO 11 K=1,I E(LN,K)=H(K,K) DO 11 M=1,I T(LNC,K,M)=R(K,M) 11 H(K,M)=0 12 CONTINUE C RETURN END C C------------------------------------------------------------------------- C SUBROUTINE UBERG(J,MIN,OB,EP,FMIN,FMAX,E,D,A,B,DM, 1DSS,LL1,LL2,LL4,KK1,KK2,MM1,MM2,SS1,SS2) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION E(12,32),D(36,32,32),A(7,12,32),DM(3),DSS(3),LL1(18), 1LL2(18),LL4(18),KK1(18),KK2(18),MM1(18),MM2(18),ST(9) INTEGER B(12,34),SS1(18),SS2(18),SIG,S3,S2,S1 INTEGER OB NYTT CHARACTER DIPCOM(3) DIPCOM(1)='a' DIPCOM(2)='b' DIPCOM(3)='c' C 104 FORMAT(1X) 103 FORMAT(' TRANSITIONS J =',I3/) 102 FORMAT(45X,6F12.7) 101 FORMAT(3I3,1X,3I3,3X,A1,F14.4,F13.6/37X,3F13.6,2H *) 100 FORMAT(3I3,1X,3I3,3X,A1,F14.4,F13.6/37X,3F13.6) 106 FORMAT(3I3,1X,3I3,3X,A1,F14.4,F13.6,2F13.6,2H *) 105 FORMAT(3I3,1X,3I3,3X,A1,F14.4,F13.6,2F13.6) C WRITE(2,103) J NNUM=1+(J-J/2*2)*18 NS=(J-J/2*2)*4 NSD=NS IF (MIN.EQ.0) GO TO 2 MAX=NNUM+5 GO TO 3 2 MAX=NNUM+17 3 DO 22 L=NNUM,MAX LLL=L/19 LL=L-LLL*18 L4=LL4(LL) DK=DM(L4)**2 M1=(J+MM1(LL))/2 M2=(J+MM2(LL))/2 IF (DK)40,40,4 4 K1=KK1(LL) K2=KK2(LL) L1=LL1(LL) L2=LL2(LL) K3=(K1*K2+1)/2 L3=(-1)**(L1-L2) IF (LL.GT.6) GO TO 6 G=(2.0D0*J+1)/(4*J*(J+1)) S3=0 GO TO 7 6 G=1.0D0/(4*J) S3=1 NSD=-NS 7 JM=J-S3 S1=SS1(LL)+NSD+(JM/2*2-JM/3*3)*4 S2=SS2(LL)+NS+(J/2*2-J/3*3)*4 IF (M1.LT.1.OR.M2.LT.1) GO TO 22 DO 18 K=1,M2 DO 18 I=1,M1 DIK=G*D(L,I,K)**2 DD=DIK*DK IF (DD-EP)18,18,8 8 DE=E(S2,K)-E(S1,I) DF=DABS(DE) IF(DF-FMIN)18,9,9 9 IF (DF-FMAX)10,10,18 10 I1=2*(I-1) I2=2*(K-1) KM1=I1+K1 KP1=JM-I1-K1+L1 KM2=I2+K2 KP2=J-I2-K2+L2 SIG=DSIGN(1.0D0,DE) DO 11 LA=1,6 11 ST(LA)=SIG*(A(LA,S2,K)-A(LA,S1,I)) ST(7)=A(7,S2,K)-A(7,S1,I) ST(7)=DABS(ST(7)) ST(8)=ST(1)*DSS(1)+ST(3)*DSS(2)+ST(5)*DSS(3) ST(9)=ST(2)*DSS(1)+ST(4)*DSS(2)+ST(6)*DSS(3) IF (B(S1,I).NE.0.OR.B(S2,K).NE.0) GO TO 16 C IF(ST(7).NE.0)THEN WRITE(2,100) J,KM2,KP2,JM,KM1,KP1,DIPCOM(L4), * DE,DD,(ST(LA),LA=7,9) ELSE WRITE(2,105) J,KM2,KP2,JM,KM1,KP1,DIPCOM(L4), * DE,DD,(ST(LA),LA=8,9) ENDIF GO TO 17 16 IF(ST(7).NE.0)THEN WRITE(2,101) J,KM2,KP2,JM,KM1,KP1,DIPCOM(L4), * DE,DD,(ST(LA),LA=7,9) ELSE WRITE(2,106) J,KM2,KP2,JM,KM1,KP1,DIPCOM(L4), * DE,DD,(ST(LA),LA=8,9) ENDIF C 17 IF (OB)18,18,12 12 WRITE(2,102) (ST(LA),LA=1,6) 18 D(L,I,K)=0.0D0 GO TO 50 40 IF (M1.LT.1.OR.M2.LT.1) GO TO 22 DO 20 K=1,M2 DO 20 I=1,M1 20 D(L,I,K)=0.0D0 50 IF (LL.EQ.6) GO TO 21 GO TO 22 21 WRITE(2,104) 22 CONTINUE C RETURN END C C------------------------------------------------------------------------- C SUBROUTINE Q2STRK(J,F,L4,M1,M2,S1,S2,A) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION F(32,32),A(7,12,32) INTEGER S1,S2 C L6=2*L4 J2=J**2 DO 2 I=1,M1 X=0.0D0 DO 1 K=1,M2 1 X=X+F(I,K) 2 A(L6,S1,I)=A(L6,S1,I)+X/J2 DO 4 K=1,M2 X=0.0D0 DO 3 I=1,M1 3 X=X-F(I,K) 4 A(L6,S2,K)=A(L6,S2,K)+X/J2 C RETURN END C C------------------------------------------------------------------------- C SUBROUTINE R2STRK(J,F,L4,M1,M2,S1,S2,A) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION F(32,32),A(7,12,32) INTEGER S1,S2 C J2=J**2 L6=2*L4 L5=L6-1 DO 2 I=1,M1 X=0.0D0 DO 1 K=1,M2 1 X=X+F(I,K) A(L5,S1,I)=A(L5,S1,I)+X 2 A(L6,S1,I)=A(L6,S1,I)-X/J2 DO 4 K=1,M2 X=0.0D0 DO 3 I=1,M1 3 X=X-F(I,K) A(L5,S2,K)=A(L5,S2,K)+X 4 A(L6,S2,K)=A(L6,S2,K)-X/J2 C RETURN END C C------------------------------------------------------------------------- C SUBROUTINE ENTART(J,JM,L,NULL,K1,K2,L1,L2,L4,S1,S2,S3, 1M1,M2,MSA,MSB,C3,EPS,DS,E,D,A,B,F,DSS,OB) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DS(3),E(12,32),D(36,32,32),A(7,12,32),DSS(3), 1IND(2,5),F(32,32),SP(64) INTEGER S1,S2,S3,ZBB,B(12,34),ZB(2),OB C 100 FORMAT(' THREE-OR FOURFOLD DEGENERATION') 101 FORMAT(2X,3I3,3X,F20.10,3I3,3X,F20.10,F20.10,12X,I4) 102 FORMAT(2X,3I3,3X,F20.10,3I3,3X,F20.10) 103 FORMAT(' STARK EFFECT, MIXED ORDER') 104 FORMAT(2X,3I3,3X,F20.10,3I3,3X,F20.10,F20.10,F12.6,3X,I4) 105 FORMAT(' 1ST ORDER STARK EFFECT') 106 FORMAT(2X,3I3,3X,F20.10,3I3,3X,F20.10,F20.10,15X,I4) 107 FORMAT(2X,3(11F10.5,/)) C IF (S3.EQ.0) GO TO 2 CJM=1.0D0/(4*(4*J**2-1)) GO TO 3 2 CJM=1.0D0/(4*((J+1)**2)) 3 CJM=DSQRT(CJM) C5=DS(L4) C4=C3*C5 J2=J**2 JT=(-1)**S3 DO 50 I=1,M1 DO 50 K=1,M2 H12=C4*D(L,I,K) H12=DABS(H12)*CJM D(L,I,K)=D(L,I,K)*NULL IF (DABS(H12).LT.1.0D-28) GO TO 50 E1=E(S1,I) E2=E(S2,K) DE=E1-E2 CR=DABS(DE)-H12*EPS IF (CR)4,4,49 4 F(I,K)=0.0D0 KM1=K1+2*(I-1) KP1=JM-2*(I-1)-K1+L1 KM2=2*(K-1)+K2 KP2=J-2*(K-1)-K2+L2 ZB(1)=B(S1,I) ZB(2)=B(S2,K) ZBB=ZB(1)+ZB(2) IF (ZBB.EQ.0) GO TO 20 DO 7 II=1,2 DO 7 KK=1,5 IND(II,KK)=ZB(II)-MSA*(ZB(II)/MSA) 7 ZB(II)=ZB(II)/MSA ID1=IND(1,4) ID2=IND(1,5) ID3=IND(2,4) ID4=IND(2,5) IF (ID1.EQ.0) GO TO 8 GO TO 9 8 EN1=0.0D0 GO TO 10 9 EN1=E(ID1,ID2) 10 IF (ID3.EQ.0) GO TO 11 GO TO 12 11 EN2=0.0D0 GO TO 13 12 EN2=E(ID3,ID4) 13 WRITE(2,100) WRITE(2,101) JM,KM1,KP1,E1,J,KM2,KP2,E2,H12,L4 WRITE(2,102) IND(1,3),IND(1,2),IND(1,1),EN1,IND(2,3),IND(2,2), 1IND(2,1),EN2 GO TO 40 20 Y1=H12/J DELTA=DE/2 Y2=DABS(DELTA)*EPS IF (Y1.GE.Y2) GO TO 30 COEF=(Y1/DELTA)**2 WRITE(2,103) WRITE(2,104) JM,KM1,KP1,E1,J,KM2,KP2,E2,DELTA,COEF,L4 IF ((OB.EQ.1).OR.(OB.EQ.-1)) GO TO 22 GO TO 29 22 DO 23 IJ=1,J JMM=(IJ-S3)**2 JMM=S3*J2+JT*JMM CO=1.0D0+JMM*COEF 23 SP(IJ)=DELTA*(DSQRT(CO)-1.0D0) WRITE(2,107) (SP(IJ),IJ=1,J) 29 GO TO 40 30 WRITE(2,105) WRITE(2,106) JM,KM1,KP1,E1,J,KM2,KP2,E2,Y1,L4 Y1=Y1*DSQRT(DSS(L4)) A(7,S1,I)=Y1 A(7,S2,K)=Y1 40 B(S1,I)=(((MSA*K+S2)*MSA+J)*MSA+KM2)*MSA+KP2 B(S2,K)=(((MSA*I+S1)*MSA+JM)*MSA+KM1)*MSA+KP1 GO TO 50 49 F(I,K)=H12**2/DE 50 CONTINUE C RETURN END C C------------------------------------------------------------------------- C SUBROUTINE HDIAG(H,N,IEGEN,U,NR) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION H(32,32),U(32,32),X(32),IQ(32) C IF (IEGEN) 15,10,15 10 DO 14 I=1,N DO 14 J=1,N IF (I-J)12,11,12 11 U(I,J)=1.0D0 GO TO 14 12 U(I,J)=0.0D0 14 CONTINUE 15 NR=0 IF (N-1)1000,1000,17 17 NMI1=N-1 DO 30 I=1,NMI1 X(I)=0.0D0 IPL1=I+1 DO 30 J=IPL1,N IF (X(I)-DABS(H(I,J)))20,20,30 20 X(I)=DABS(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,148 148 CONTINUE 120 NR=NR+1 130 IF (H(IPIV,IPIV)-H(JPIV,JPIV))140,141,142 140 TANG=2.0D0*H(IPIV,JPIV)/(H(IPIV,IPIV)-H(JPIV,JPIV) 1-DSQRT((H(IPIV,IPIV)-H(JPIV,JPIV))**2+4.0D0*H(IPIV,JPIV)**2)) GO TO 150 141 TANG=1.0D0 GO TO 150 142 TANG=2.0D0*H(IPIV,JPIV)/(H(IPIV,IPIV)-H(JPIV,JPIV) 1+DSQRT((H(IPIV,IPIV)-H(JPIV,JPIV))**2+4.0D0*H(IPIV,JPIV)**2)) 150 COSINE=1.0D0/DSQRT(1.0D0+TANG**2) SINE=TANG*COSINE HII=H(IPIV,IPIV) H(IPIV,IPIV)=COSINE**2*(HII+TANG*(2.0D0*H(IPIV,JPIV)+TANG*H(JPIV, *JPIV))) H(JPIV,JPIV)=COSINE**2*(H(JPIV,JPIV)-TANG*(2.0D0*H(IPIV,JPIV)-TANG 1*HII)) H(IPIV,JPIV)=0.0D0 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=DSIGN(1.0D0,-SINE)*COSINE COSINE=DABS(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)230,240,230 230 IF (IQ(I)-JPIV)350,240,350 240 K=IQ(I) 250 HTEMP=H(I,K) H(I,K)=0.0D0 IPL1=I+1 X(I)=0.0D0 DO 320 J=IPL1,N IF (X(I)-DABS(H(I,J)))300,300,320 300 X(I)=DABS(H(I,J)) IQ(I)=J 320 CONTINUE H(I,K)=HTEMP 350 CONTINUE X(IPIV)=0.0D0 X(JPIV)=0.0D0 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)-DABS(H(I,IPIV)))380,390,390 380 X(I)=DABS(H(I,IPIV)) IQ(I)=IPIV 390 H(I,JPIV)=-SINE*HTEMP+COSINE*H(I,JPIV) IF (X(I)-DABS(H(I,JPIV)))400,530,530 400 X(I)=DABS(H(I,JPIV)) IQ(I)=JPIV GO TO 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)-DABS(H(IPIV,I)))440,450,450 440 X(IPIV)=DABS(H(IPIV,I)) IQ(IPIV)=I 450 H(I,JPIV)=-SINE*HTEMP+COSINE*H(I,JPIV) IF (X(IPIV)-DABS(H(I,JPIV)))460,530,530 460 X(I)=DABS(H(I,JPIV)) IQ(I)=JPIV GO TO 530 480 HTEMP=H(IPIV,I) H(IPIV,I)=COSINE*HTEMP+SINE*H(JPIV,I) IF (X(IPIV)-DABS(H(IPIV,I)))490,500,500 490 X(IPIV)=DABS(H(IPIV,I)) IQ(IPIV)=I 500 H(JPIV,I)=-SINE*HTEMP+COSINE*H(JPIV,I) IF (X(JPIV)-DABS(H(JPIV,I)))510,530,530 510 X(JPIV)=DABS(H(JPIV,I)) IQ(JPIV)=I 530 CONTINUE IF (IEGEN)40,540,40 540 DO 550 I=1,N HTEMP=U(I,IPIV) U(I,IPIV)=COSINE*HTEMP+SINE*U(I,JPIV) 550 U(I,JPIV)=-SINE*HTEMP+COSINE*U(I,JPIV) GO TO 40 C 1000 RETURN END C C------------------------------------------------------------------------- C-------------------------------------------------------------------------