CFTN77,Q,T,S C$CDS ON C INCLUDE 'CALPGM.INC' C PACKAGE FOR DOUBLET PI : DPI.FOR C C Copyright (C) 1989, California Institute of Technology C All rights reserved. U. S. Government Sponsorship under C NASA Contract NAS7-918 is acknowledged C C Herbert M. Pickett 20 Mar 1989 C SUBROUTINE HAMX(IBLK,NSIZD,IDPAR,PAR,EGY,T,DEDP,IFDUMP) C C .. PACKAGE FOR DOUBLET PI C INTEGER IBLK,NSIZD,IFDUMP,NSIZE,ITYMAX,IDM,IERR REAL*8 IDPAR(0:10),PAR(*),EGY(*),T(NSIZD,*),DEDP(NSIZD,*) PARAMETER (ITYMAX=26,IDM=8) INTEGER*2 IP(ITYMAX),ISBLK(IDM+1),IWORK(IDM) REAL*8 W(ITYMAX,IDM,IDM),WORK(IDM),ZV(2) INTEGER IWHOLE,ISDGN,AII,F,AJI,AJF,I0,I1,I2,L,K1,K3,JQ(IDM), & NPAR,ISIG,I,J,K,IDGN,JS,IQ,IPAR,II,KK,ITMP,KQ,JJ,I4, & JF,OFDIAG,KPAR PARAMETER (I0=0,I1=1,I2=2) REAL*8 PTY,XJQ,RX,RY,SQJ,X,Z,AA,BB,AB,RT2,AI,FG, & DJ1,DJ2,G,GG,C6JJ,C3JJ,RT,Q,QQ,DDOT COMMON /SQN/IWHOLE,ISDGN,AII C EMA IDPAR,PAR,EGY,DEDP,T,WORK,W,IWORK SAVE NPAR,ISIG,IP,/SQN/ SAVE W,ZV DATA ZV/0.,0./ NSIZE=NSIZD IF(IBLK.EQ.0) THEN I=MIN(ISDGN+ISDGN,IDM) IF(NSIZE.GT.I) NSIZE=I PAR(1)=NSIZE C UNSCRAMBLE PARAMETER LIST NPAR=IDPAR(0) DO 4 I=1,ITYMAX DO 5 J=1,NPAR K=ABS(IDPAR(J)) IF(K.EQ.I) THEN IP(I)=J GO TO 4 ENDIF 5 CONTINUE IP(I)=0 4 CONTINUE ISIG=1 IF(IP(1).NE.0) ISIG=0 RETURN ENDIF IF(ISIG.EQ.0) THEN K=IP(1) ISIG=1 IF(PAR(K).GT.0) ISIG=-1 ENDIF CALL GETQN(IBLK,I1,JQ,IDGN,I) F=IDGN-1 JS=JQ(1) I =JQ(3) PTY=I IQ=ISIG*I*(1-2*MOD(JS,2)) DO 1 I=1,NSIZE,2 C ....... SET UP J QUANTUM NUMBER JQ(I)=JS*IQ JQ(I+1)=JS IQ=-IQ 1 JS=JS+1 C .. ZERO HAMILTONIAN AND DERIVATIVE W DO 3 I=1,NSIZE DO 3 J=I,NSIZE DO 3 IPAR=1,NPAR 3 W(IPAR,I,J)=0. C ... SET UP FINE STRUCTURE TERMS DO 6 I=1,NSIZE,2 II=I+1 XJQ=JQ(I)+1 RX=JQ(I)*JQ(I) RY=RX-2. SQJ=RX-0.25 X=RX-1. Z=SQRT(X) DO 7 KK=1,10 K=IP(KK) IF(K.EQ.0) GO TO 7 GO TO (10,20,30,40,50,60,70,80,90,100),KK C..A 10 AA=-(1+ISIG)/2 BB=1+AA OFDIAG=0 GO TO 8 C..AJ 20 AA=-RX BB=RY OFDIAG=0 GO TO 8 C..AH 30 AA=RX*RX+X BB=-RY*RY-X AB=0. OFDIAG=1 GO TO 8 C..B+0.5*Q 40 AA=RX BB=RY AB=Z OFDIAG=1 GO TO 8 C..D 50 AA=-RX*RX-X BB=-RY*RY-X AB=-2.*Z*X OFDIAG=1 GO TO 8 C..H 60 AA=(RX*RX+X)*RX+2.*X*X BB=(RY*RY+X)*RY+2.*X*X AB=((3.*X+1.)*X-1.)*Z OFDIAG=1 GO TO 8 C..P 70 AA=0.5*XJQ BB=0. AB=0.25*Z OFDIAG=1 GO TO 8 C..Q 80 AA=0.5+JQ(I) BB=0.5 AB=0.5*Z*JQ(I) OFDIAG=1 GO TO 8 C..PD 90 AA=0.5*XJQ*SQJ BB=0. AB=0.5*Z*SQJ OFDIAG=1 GO TO 8 C..QD 100 AA=0.5*XJQ*XJQ*SQJ BB=0.5*X*SQJ AB=0.5*Z*XJQ*SQJ OFDIAG=1 C 8 IF(OFDIAG.NE.0) W(K,I,II)=AB W(K,I,I)=AA W(K,II,II)=BB 7 CONTINUE 6 CONTINUE C C CHI PARAMETERS (MAGNETIC INTERACTIONS) RT2=SQRT(2.D0) AI=0.5*AII FG=AI*(AI+1.)*(AI+AI+1.) DO 101 I=1,NSIZE,2 AJI=JQ(I+1)*2 -1 DJ1=AJI+1 DO 102 J=I,NSIZE,2 AJF=JQ(J+1)*2-1 DJ2=AJF+1 IF(ABS(AJF-AJI).GT.2) GO TO 102 G=SQRT(DJ1*DJ2*FG)*C6JJ(AJI,I2,AJF,AII,F,AII) ITMP=(AJI+F+AII)/2 IF(MOD(ITMP,I2).NE.0) G=-G GG=G ITMP=AJF/2 IF(MOD(ITMP,I2).NE.0) GG=-GG K=IP(11) K1=-1 IF(K.NE.0)W(K,I,J)=GG*2.*C3JJ(AJI,I2,AJF,K1,I0,I1) K=IP(12) IF(K.NE.0)W(K,I,J)=G*PTY*RT2*C3JJ(AJI,I2,AJF,K1,I2,K1) K=IP(13) K3=3 K1=-K3 IF(K.NE.0)W(K,I+1,J+1)=-C3JJ(AJI,I2,AJF,K1,I0,K3) & *GG/1.5 K=IP(14) IF(K.EQ.0) GO TO 102 W(K,I,J+1)=GG*RT2*C3JJ(AJF,I2,AJI,K1,I2,I1) IF(J.EQ.I) GO TO 102 W(K,I+1,J)=GG*RT2*C3JJ(AJI,I2,AJF,K1,I2,I1) 102 CONTINUE 101 CONTINUE DO 103 I=1,NSIZE,2 Z=JQ(I)*JQ(I)-1 DO 104 KQ=15,19 K=IP(KQ) IF(K.EQ.0) GO TO 104 IF(KQ.NE.19) THEN AA=Z KK=KQ-4 ELSE AA=-JQ(I) KK=14 ENDIF J=IP(KK) IF(J.EQ.0) GO TO 104 JJ=KQ-14 GO TO (106,106,107,108,108),JJ 106 W(K,I,I)=AA*W(J,I,I) IF(JJ.NE.1) GO TO 104 J=IP(13) IF(J.EQ.0) GO TO 104 C 107 W(K,I+1,I+1)=AA*W(J,I+1,I+1) GO TO 104 C 108 W(K,I,I+1)=AA*W(J,I,I+1) 104 CONTINUE 103 CONTINUE C C ETA PARAMETERS (QUADRUPOLE INTERACTIONS) IF(AII.LT.2) GO TO 200 FG=FG*(AI+1.5)/(AI*AI*(AI-0.5)) RT=PTY/SQRT(24.D0) L=4 DO 111 I=1,NSIZE,2 AJI=JQ(I+1)*2-1 DJ1=AJI+1 DO 112 J=I,NSIZE,2 AJF=JQ(J+1)*2-1 DJ2=AJF+1 IF(ABS(AJF-AJI).GT.L) GO TO 112 Q=SQRT(DJ1*DJ2*FG)*C6JJ(AJI,L,AJF,AII,F,AII) ITMP=(AJI+F+AII)/2 IF(MOD(ITMP,I2).NE.0) Q=-Q QQ=Q ITMP=AJF/2 IF(MOD(ITMP,I2).NE.0) QQ=-QQ K=IP(20) K1=-1 IF(K.NE.0)W(K,I,J)=0.25*QQ*C3JJ(AJI,L,AJF,K1,I0,I1) K=IP(21) K3=3 K1=-K3 IF(K.NE.0)W(K,I+1,J+1)=-0.25*QQ*C3JJ(AJI,L,AJF,K1,I0,K3) K=IP(22) IF(K.EQ.0) GO TO 112 K3=-1 W(K,I,J+1)=-RT*Q*C3JJ(AJI,L,AJF,K3,L,K1) IF(I.EQ.J) GO TO 112 I4=4 W(K,I+1,J)=RT*Q*C3JJ(AJI,I4,AJF,K1,L,K3) 112 CONTINUE 111 CONTINUE DO 113 I=1,NSIZE,2 Z=JQ(I)*JQ(I)-1 SQJ=Z+0.75 DO 118 KQ=23,26 K=IP(KQ) IF(K.EQ.0) GO TO 118 KK=KQ-22 GO TO (114,115,116,117),KK 114 AA=SQJ BB=SQJ AB=SQJ/JQ(I) OFDIAG=1 GO TO 119 115 AA=JQ(I) BB=AA*Z/(Z-6.) OFDIAG=0 GO TO 119 116 AA=0. BB=0. AB=-1./JQ(I) OFDIAG=1 GO TO 119 117 AA=0. BB=Z/(6.-Z) OFDIAG=0 GO TO 119 119 IF(OFDIAG.NE.0) THEN KK=IP(22) W(K,I,I+1)=AB*W(KK,I,I+1) ENDIF KK=IP(21) W(K,I+1,I+1)=BB*W(KK,I+1,I+1) KK=IP(20) W(K,I,I)=AA*W(KK,I,I) 118 CONTINUE 113 CONTINUE C C BRING H TOGETHER C 200 DO 201 I=1,NSIZE ISBLK(I)=I DO 201 J=I,NSIZE 201 T(J,I)=DDOT(NPAR,W(1,I,J),I1,PAR,I1) IF(IFDUMP.NE.0) THEN C special to dump Hamiltonian DO 300 I=1,NSIZE EGY(I)=T(I,I) CALL DCOPY(NPAR,ZV,I0,DEDP,NSIZE) 300 CONTINUE RETURN ENDIF C C DIAGONALIZE AND ORDER EIGENVECTORS C ISBLK(NSIZE+1)=NSIZE+1 CALL HDIAG(NSIZD,NSIZE,T,EGY,WORK,IERR) CALL ORDBLK(NSIZD,NSIZE,T,EGY,NSIZE,ISBLK,WORK,IWORK,I0) C C TRANSFORM W TO DEDP KPAR=0 DO 350 IPAR=1,NPAR IF(IDPAR(IPAR).GE.0) THEN KPAR=KPAR+1 KK=IPAR ELSE AB=PAR(IPAR)/PAR(KK) ENDIF DO 351 I=1,NSIZE BB=0. DO 352 K=1,NSIZE AA=0. JF=K-1 DO 354 J=1,JF 354 AA=AA+W(IPAR,J,K)*T(J,I) DO 355 J=K,NSIZE 355 AA=AA+W(IPAR,K,J)*T(J,I) 352 BB=BB+AA*T(K,I) IF(IPAR.EQ.KK) THEN DEDP(I,KPAR)=BB ELSE DEDP(I,KPAR)=DEDP(I,KPAR)+BB*AB ENDIF 351 CONTINUE 350 CONTINUE RETURN END C FUNCTION INTENS(IBLK,ISIZ,JBLK,JSIZ,IDIP,DIP,S,NSIZD) C C ... PACKAGE FOR DOUBLET PI ...ELECTRIC DIPOLE C INTEGER INTENS,IBLK,ISIZ,JBLK,JSIZ,NSIZD INTEGER*4 IDIP(0:*) REAL*8 S(NSIZD,*),DIP(*) C EMA IDIP,DIP,S INTEGER IWHOLE,ISDGN,AII,JJI,JJJ,FI,FJ,I2, & IDGN,JDGN,IX,JX,I,J,K,IM REAL*8 DGN,FAC,AJI,OMEG,C6JJ,AJJ COMMON /SQN/IWHOLE,ISDGN,AII SAVE /SQN/,I2 DATA I2/ 2/ INTENS=0 C START OF INTENSITY CALCULATION IF(IDIP(0).EQ.0) RETURN IF(MOD(IBLK+JBLK,2).EQ.0) RETURN IDGN=((IBLK-1)/2)*2+IWHOLE JDGN=((JBLK-1)/2)*2+IWHOLE IF(ABS(IDGN-JDGN).GT.2) RETURN FI=IDGN-1 FJ=JDGN-1 DGN=IDGN*JDGN IX=ABS(IDGN-ISDGN) JX=ABS(JDGN-ISDGN)-IX+2 FAC=DIP(1)*SQRT(DGN) DO 101 I=1,ISIZ IM=I-1 K=MOD(IM,I2) JJI=IX+IM-K AJI=0.5*JJI IF(K.EQ.0) FAC=-FAC OMEG=K+0.5 DO 102 J=1,JSIZ 102 S(I,J)=0. J=I-JX IF(J.GT.JSIZ) GO TO 101 IF(J.GT.0) THEN JJJ=JJI-2 S(I,J)=FAC*SQRT((AJI-OMEG)*(AJI+OMEG)/AJI)* 1 C6JJ(JJI,I2,JJJ,FJ,AII,FI) ENDIF J=J+2 IF(J.GT.JSIZ) GO TO 101 IF(J.GT.0.AND.JJI.GT.K) THEN S(I,J)=FAC*OMEG*SQRT((AJI+AJI+1.)/(AJI*AJI+AJI))* 1 C6JJ(JJI,I2,JJI,FJ,AII,FI) ENDIF J=J+2 IF(J.GT.0.AND.J.LE.JSIZ) THEN JJJ=JJI+2 AJJ=AJI+1. S(I,J)=-FAC*SQRT((AJJ-OMEG)*(AJJ+OMEG)/AJJ)* 1 C6JJ(JJI,I2,JJJ,FJ,AII,FI) ENDIF 101 CONTINUE INTENS=IDGN RETURN END C SUBROUTINE SETINT(IFDIAG,NSAV,IDIP) C DOES NOTHING TO INITIALIZE INTEGER IFDIAG,NSAV,NDIP INTEGER*4 IDIP(0:*) C EMA IDIP IFDIAG=1 NSAV=1 NDIP=IDIP(0) RETURN END SUBROUTINE GETQN(IBLK,INDX,IQN,IDGN,IQ) C C ..PACKAGE FOR DOUBLET PI C INTEGER IBLK,INDX,IDGN,IQ,IQN(6) INTEGER IWHOLE,ISDGN,AI,KF,MOD COMMON /SQN/IWHOLE,ISDGN,AI SAVE /SQN/ KF=(IBLK-1)/2 IDGN=KF+KF+IWHOLE IF(INDX.EQ.0) THEN IDGN=2*MIN(ISDGN,IDGN) RETURN ENDIF IQ=1 IQN(1)=ABS(IDGN-ISDGN)/2+(INDX+1)/2 IQN(2)=MOD(INDX-1,2)+1 IQN(3)=1-2*MOD(IBLK,2) IQN(4)=KF IF(IQN(1).LT.IQN(2)) IDGN=0 RETURN END C FUNCTION SETOPT(LUIN,IQNFMT,ITD,NAMFIL) C SET OPTIONS INTEGER SETOPT,LUIN,IQNFMT,ITD CHARACTER*80 NAMFIL,CARD INTEGER IWHOLE,ISDGN,AII,I,PCARD REAL*8 VAL(1) COMMON /SQN/IWHOLE,ISDGN,AII SAVE /SQN/ CALL CJJINI READ(LUIN,'(A)',END=999) CARD VAL(1)=1 IWHOLE=1 I=PCARD(CARD,VAL,IWHOLE) ISDGN=VAL(1) IF(ISDGN.LT.1) ISDGN=1 NAMFIL='dpi.nam' ITD=2 AII=ISDGN-1 IWHOLE=MOD(AII,2) IQNFMT=804 IF(IWHOLE.EQ.0) IQNFMT=IQNFMT+10 IF(ISDGN.EQ.1) IQNFMT=803 SETOPT=1 RETURN 999 SETOPT=-1 RETURN END FUNCTION SETBLK(LU,IDPAR,PAR,NBLKPF,NEGY) C SET OPTIONS INTEGER SETBLK,LU,NBLKPF,NEGY REAL*8 IDPAR(*),PAR(*),DMY(2) C EMA IDPAR,DMY REAL*4 AI INTEGER ITMP,IWHOLE,ISDGN,AII COMMON /SQN/IWHOLE,ISDGN,AII SAVE /SQN/ NBLKPF=2 AI=0.5*AII WRITE(LU,900) AI 900 FORMAT(' doublet PI lines for I =',F5.1) ITMP=0 CALL HAMX(ITMP,NEGY,IDPAR,DMY,PAR,PAR,PAR,ITMP) NEGY=DMY(1) SETBLK=1 RETURN END