SUBROUTINE BODFIX CHARACTER*80 Label DIMENSION F(19,10,10),UL(19,10,10),SSQ(19,10,10),SRE(19,10,10), 1 SIM(19,10,10), CO(19,10),SI(19,10),ETA(19,10),RT(19,10),C(10,10), 2 U(10,10),EIGEN(10) Logical Dbug COMMON/Blank/ E,B,PK,PL,QSQ(19,10),QSS(19,10) COMMON/B2/EIG(19,10) COMMON/B3/A12,A6 EQUIVALENCE(F,SRE) DBug=.false. READ(5,'(A)')Label WRITE(6,'(A)')Label C------------------------------------------------------- C read in data C read parameters controlling accuracy of jacobi etc C------------------------------------------------------- READ(5,10) ITMAX,EPS1,EPS2,EPS3 10 FORMAT(I5,3E10.3) C------------------------------------------------------- C nx is the number of points to be taken in the gauss mehler integration C------------------------------------------------------- READ(5,20) NX 20 FORMAT(I5) INDEX = 0 C------------------------------------------------------- C b and e are the well depth and energy in reduced units C------------------------------------------------------- READ(5,30) B,E 30 FORMAT(2F10.7) C------------------------------------------------------- C asymmetry parameters of the repulsive and attractive potentials C------------------------------------------------------- READ(5,40) A12,A6 40 FORMAT(2F10.4) C------------------------------------------------------- C jmin is smallest rotor state. jmax is largest. lam is component C of angular momenta along body fixed axis. lamax is less than or C equal to jmax C lamax need only be as large as the largest lower j of interest C for example if only cross sections from j=0 to j=n are wanted C set lamax =0 for any n C------------------------------------------------------- READ(5,50) JMIN,JMAX,LAMAX 50 FORMAT(3I5) C------------------------------------------------------- C jc is the total angular momentum C read in the maximum and minimum values of jc desired and the grid C spacing jcjump desired C------------------------------------------------------- READ(5,50) JCMIN,JCMAX,JCJUMP WRITE(6,60) 60 FORMAT(1x,'BODY FIXED INFINITE ORDER SUDDEN APPROXIMATION') WRITE(6,70) ITMAX,EPS1,EPS2,EPS3 70 FORMAT(1x,' ITMAX=',I4,' EPS1=',E12.3,' EPS2=',E12.3, > ' EPS3=',E12.3) PK = SQRT(E ) WRITE(6,80) B,E,PK 80 FORMAT(1x,'B=',E15.7,' E=',E15.7,' PK=',E15.7) WRITE(6,90) A12,A6 90 FORMAT(1x, ' A12= ',F10.5,5X, ' A6= ',F10.5 / 1H ) WRITE(6,100) JMIN,JMAX,LAMAX 100 FORMAT(1x,'JMIN=',I5,' JMAX=',I5,' LAMAX=',I5) C------------------------------------------------------- C set up initial values and estimates C------------------------------------------------------- XR = (1.+SQRT(1.+E/B))**(-1./6.) PI=3.1415927 N=(JMAX-JMIN)/2 +1 LM = LAMAX+1 K=2 DO 140 L=1,LM LAM = L-1 C------------------------------------------------------- C note the integer division C------------------------------------------------------- IM= 1+L/2 C------------------------------------------------------- C set up f matrices C construc matrices from legendre poly of order k in potential C------------------------------------------------------- 110 CALL POT(F,K,L,N,JMIN,IM) C------------------------------------------------------- C diagonalize the f matrices C------------------------------------------------------- DO 120 I=IM,N RT(L,I) = XR EIG(L,I)=0.0 DO 120 J=IM,N UL(L,I,J) =0. 120 C(I,J) = F(L,I,J) CALL JACOBI(C,U,IM,N,EIGEN,EPS1,EPS2,EPS3,ITMAX) DO 130 I=IM,N EIG(L,I)=EIGEN(I) DO 130 J=IM,N 130 UL(L,I,J)=U(I,J) 140 CONTINUE C------------------------------------------------------- C start calculation for a given total j C------------------------------------------------------- SIGG=0.0 SIGU=0.0 Sigma_Sum0=0.0 Sigma_Sum2=0.0 Sigma_Sum4=0.0 Sigma_Sum6=0.0 DO 290 JC=JCMIN,JCMAX,JCJUMP PL=JC WRITE(6,*) WRITE(6,*) WRITE(6,150) JC 150 FORMAT(1x,'JTOT=',I5) C------------------------------------------------------- C calculate turning points C------------------------------------------------------- XJ = JC**2 DO 260 L=1,LM LAM = L-1 IM= 1+L/2 DO 200 I=IM,N XL = RT(L,I) CALL POTS(XL,L,I) YL=QSQ(L,I) XR = RT(L,I) +0.001 +(0.6E-06)*XJ CALL POTS(XR,L,I) YR=QSQ(L,I) DO 170 J=1,ITMAX RR=(XL*YR-XR*YL)/(YR-YL) CALL POTS(RR,L,I) FR=QSQ(L,I) IF(ABS(FR).LE.0.1E-04) GO TO 180 IF( ABS(XR-RR) .LT. 1.5E-7 ) GOTO 180 IF(ABS(YL).LT.ABS(YR)) GO TO 160 XL=XR YL=YR 160 CONTINUE XR=RR 170 YR=FR 180 RT(L,I)=RR IF(Dbug)WRITE(6,190) L,I,RT(L,I),J,FR 190 FORMAT(1H ,'LAM +1=',I5,' I=',I5,' RT=',E14.7,' ITER=',I5, > ' FR=',E14.7) 200 CONTINUE C------------------------------------------------------- C calculate phase shifts C------------------------------------------------------- DO 210 I=IM,N CALL GAUMEH(L,I,RT(L,I),ETA(L,I),NX,INDEX) 210 IF(Dbug)WRITE(6,220) NX,L,I,ETA(L,I) 220 FORMAT(1H ,'NX=',I5,' L=',I5,' I=',I5,' ETA=',E14.7) C------------------------------------------------------- C calculate scattering matrix C------------------------------------------------------- DO 230 I=IM,N PH=2.*ETA(L,I) CO(L,I)=COS(PH) SI(L,I) = SIN(PH) DO 230 J=IM,N SSQ(L,I,J)=0. SRE(L,I,J)=0. 230 SIM(L,I,J)=0. DO 250 I=IM,N DO 250 J=I,N DO 240 M=IM,N UU = UL(L,I,M)*UL(L,J,M) SRE(L,I,J)=SRE(L,I,J)+UU*CO(L,M) 240 SIM(L,I,J)=SIM(L,I,J)+UU*SI(L,M) 250 SSQ(L,I,J)=SRE(L,I,J)**2+SIM(L,I,J)**2 260 CONTINUE OP0 = 1.0+SSQ(1,1,1)-2.*SRE(1,1,1) OP2 = SSQ(1,1,2) OP4 = SSQ(1,1,3) OP6 = SSQ(1,1,4) TOP0 = (2.*FLOAT(JC)+1.)*OP0 TOP2 = (2.*FLOAT(JC)+1.)*OP2 TOP4=(2.*FLOAT(JC)+1.)*OP4 TOP6=(2.*FLOAT(JC)+1.)*OP6 WRITE(6,270) JC,OP0 ,TOP0,OP2,TOP2,OP4,TOP4,OP6,TOP6 WRITE(6,'(i5,",",8(1pe15.7,","))')JC, OP0, OP2, OP4, OP6, > TOP0,TOP2,TOP4,TOP6 WRITE(8,'(i5,",",8(1pe15.7,","))')JC, OP0, OP2, OP4, OP6, > TOP0,TOP2,TOP4,TOP6 270 FORMAT(1x,' J=',I5,' TSQ=',2F12.5,' OP2=',2F10.5,' OP4=', 1 2F10.5,' OP6=',2F10.5) IF(MOD(JC,2).EQ.0)SIGG=SIGG+TOP0 IF(MOD(JC,2).EQ.1)SIGU=SIGU+TOP0 WRITE(6,280)SIGG,SIGU 280 FORMAT(1x,' SIGG=',E14.7,' SIGU=',E14.7) WRITE(9,'(i5,",",4(1pe15.7,","))')JC, TOP0, TOP2, TOP4, TOP6 !IF(MOD(JC,2).EQ.0)THEN Sigma_Sum0=Sigma_Sum0+TOP0 Sigma_Sum2=Sigma_Sum2+TOP2 Sigma_Sum4=Sigma_Sum4+TOP4 Sigma_Sum6=Sigma_Sum6+TOP6 !ENDIF 290 CONTINUE WRITE(6,'(" Sigma_Sum0=",1PE15.7)')Sigma_Sum0*PI/E/0.529**2 WRITE(6,'(" Sigma_Sum2=",1PE15.7)')Sigma_Sum2*PI/E/0.529**2 WRITE(6,'(" Sigma_Sum4=",1PE15.7)')Sigma_Sum4*PI/E/0.529**2 WRITE(6,'(" Sigma_Sum6=",1PE15.7)')Sigma_Sum6*PI/E/0.529**2 RETURN END