C N. R. BADNELL UoS v1.9 29/06/21 C C BASED-ON MJS VERSION LAST MODIFIED 13.07.05 C C STGFF C ***** C PROGRAM MAIN C IMPLICIT REAL*8 (A-H,O-Y) IMPLICIT COMPLEX*16(Z) C LOGICAL ALLE,ALLO,KPP,QDT1,QDT2,ALL,TOP C CHARACTER*1 NUM(0:9) C SUN REAL*4 TARRY(2) C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) PARAMETER (MXLSD=(MZSLP*(MZSLP+1))/2) C DIMENSION ISE(MZSLP),ILE(MZSLP),ISO(MZSLP),ILO(MZSLP), +KFILE(MZSLP),KFILO(MZSLP) C DIMENSION ND(MXLSD),KPP(MXLSD),KK1(MXLSD),KK2(MXLSD) DIMENSION DR(MZCHF,MZCHF),DI(MZCHF,MZCHF),FT(MZMSH), + SIGMA(MZMSH),OMEGA(MZMSH),GAUNT(MZMSH) C COMMON/LIST/ISEE(MZSLP),ILEE(MZSLP),KFILEE(MZSLP),KEE, + ISOO(MZSLP),ILOO(MZSLP),KFILOO(MZSLP),KOO COMMON/MISC/WBODE(MZPTS),IPERT,MXE,EMESH(MZMSH),BSTO COMMON/CHAN1/ECH1(MZCHF),CC1(MZCHF),LLCH1(MZCHF), X L2P1(MZCHF),KJ1(MZCHF), + NCHF1 COMMON/CHAN2/ECH2(MZCHF),CC2(MZCHF),LLCH2(MZCHF), X L2P2(MZCHF),KJ2(MZCHF), + NCHF2 COMMON/CBUT/IBUT,K10,K20 COMMON/CD/DN(MZCHF,MZCHF),DA(MZCHF,MZCHF),STA(MZMSH),STB(MZMSH), + TOT(MZMSH) COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,TINY,DELTA,LACC COMMON/DEP1/RK1(MZCHF,MZCHF),SE(MZPTS,MZCHF),CE(MZPTS,MZCHF), + NCHOP1,QDT1,F1(MZCHF,MZCHF),FP1(MZCHF,MZCHF),FKNU1(MZCHF) COMMON/DEP2/RK2(MZCHF,MZCHF),SO(MZPTS,MZCHF),CO(MZPTS,MZCHF), + NCHOP2,QDT2,F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FKNU2(MZCHF) COMMON/DIFF1/SEP(MZCHF),CEP(MZCHF) COMMON/DIFF2/SOP(MZCHF),COP(MZCHF) COMMON/CPOINT/RZERO,RTWO,H,KP2 COMMON/CBODE/BDM2(MZPTS),BD0(MZPTS),BD1(MZPTS) COMMON/CALP/ASS(MZCHF,MZCHF),ASC(MZCHF,MZCHF),ACS(MZCHF,MZCHF), + ACC(MZCHF,MZCHF) COMMON/CDKKP/MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, 1 MNP2D2,NCHND2,LRGLD2,NPTYD2,NPAD, 2 DKL(MNPEXT,MNPEXT), CV X DKV(MNPEXT,MNPEXT), 3 AMAT(MZCHF,MZCHF),BMATL(MZCHF,MZCHF) CV X,BMATV(MZCHF,MZCHF) COMMON/SCALE/AZ,AZP,AZ2,AZ4,AZP2 COMMON/CTARG/ENAT(MZTAR),NZED,NELC,NAST,NCONAT(MZTAR), + IST(MZTAR),ILT(MZTAR),IPT(MZTAR) COMMON/NRBZED/TZED C NAMELIST/STGFF/IBUT,ACC,ALLE,ALLO,IPRINT NAMELIST/MESH/M11,M12,M21,M22,NTEMP,TOP C DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ C C FILES FOR INPUT DATA C 0 F00 FOR F INDEX, FROM STGF. C 1 Fnn OUTER REGION DATA, FROM STGF. C 3 D00 DIPOLE INDEX, FROM STG3. C 4 Dmm DIPOLE DATA, FROM STG3 C C DATA READ FROM UNIT 5 C C NAMELIST STGFF: IBUT,AC,ALLE,ALLO C C IBUT= 1 FOR BUTTLE, 0 NOT (DEFAULT) C AC = ACCURACY PARAMETER, 1.E-4 DEFAULT C ALLE C FOR ALLE=TRUE, PROCESS ALL EVEN PARITY SLP CASES. DEFAULT C FOR ALLE=FALSE, READ VALUES OF (ISE, ILE, 0) FOLLOWED C BY (-1 -1 -1) TERMINATOR. C ALLO C FOR ALLO=TRUE, PROCESS ALL ODD PARITY SLP CASES. DEFAULT C FOR ALLO=FALSE, READ VALUES OF (ISO, ILO, 1) FOLLOWED C BY (-1 -1 -1) TERMINATOR. C C NAMELIST MESH: M11,M12,M21,M22,NTEMP,TOP C C ENERGY RANGE (SUBSET OF M1,M2=1,MXE STGF ENERGIES; DEFAULT ALL) C M11, M12, M21, M22 FOR ENERGY RANGES. C M1 IN RANGE M11 TO M12 FOR EVEN. C M2 IN RANGE M21 TO M22 FOR ODD. C C NTEMP C TEMPERATURES IN KELVIN READ (ONE PER RECORD) DEFAULT 0. C TOP-UP IS AUTOMATICALLY INCLUDED, UNLESS C EXPLICITLY SWITCHED-OFF VIA TOP.EQ..FALSE.. C TOP=.TRUE. (DEFAULT) TOPS-UP THERMALLY-AVERAGED XSCNS, NTEMP>0. C C C OUTPUT FILES c ON UNIT 70: C BF.ISEILE_ISOILO FOR TRANSITION (ISE,ILE) TO (ISO, OLO). C VALUES OF: E1,E2 AND X WHERE !mjs C X LINE-STRENGTH FOR TRANSITION, Z-SCALED AND WITH !mjs c ORIGINAL NORMALISATION !mjs C BF FILES WRITTEN ONLY FOR NTEMP.LE.0. C ON UNIT 80: C FF.n WHERE n=1, 2, 3, ...9 FOR THE NTEMP TEMPERATURES READ: C VALUE OF TEMP, THEN VALUES OF OM, SIG, SIG**3, GAUNT: WHERE SIG C HERE IS THERMALLY-AVERAGED FF CROSS-SECTION. C C FACT=ALPHA_fsc/6 DATA FACT/1.216226e-3/ C TFACT=2*PI*SQRT(PI*157894) DATA TFACT/4425.25/ C FACTG=2*(2*PI/3)**(3/2)*ALPHA_fsc DATA FACTG/0.044237/ C GFACT=3**(1.5)/(ALPHA_fsc*16*PI) DATA GFACT/14.1660/ C C OPEN(5,FILE='dstgff') OPEN(6,FILE='routff',STATUS='UNKNOWN') C WRITE(6,200) C IBUT=0 AC=1.D-4 ALLE=.TRUE. ALLO=.TRUE. IPRINT=0 C READ(5,STGFF) C CALL ACSUB C C READ INDEX FILE FOR FREE STATES C CALL RF00 C C NEUTRAL C IF(NZED.EQ.NELC)THEN TZED=0. ELSE TZED=1. ENDIF C C READ TARGET SLPI C CALL TARG C C READ CASES TO BE PROCESSED. C C EVEN PARITY C IF(ALLE)THEN DO K=1,KEE ISE(K)=ISEE(K) ILE(K)=ILEE(K) KFILE(K)=KFILEE(K) ENDDO KE=KEE ELSE K=0 1 K=K+1 READ(5,*)ISE(K),ILE(K) IF(ILE(K).GE.0)THEN DO L=1,KEE IF(ISE(K).EQ.ISEE(L).AND.ILE(K).EQ.ILEE(L))THEN KFILE(K)=KFILEE(L) GOTO 2 ENDIF ENDDO WRITE(6,600)K,ISE(K),ILE(K) STOP 2 GOTO 1 ENDIF KE=K-1 ENDIF C C ODD PARITY C IF(ALLO)THEN DO K=1,KOO ISO(K)=ISOO(K) ILO(K)=ILOO(K) KFILO(K)=KFILOO(K) ENDDO KO=KOO ELSE K=0 3 K=K+1 READ(5,*)ISO(K),ILO(K) IF(ILO(K).GE.0)THEN DO L=1,KOO IF(ISO(K).EQ.ISOO(L).AND.ILO(K).EQ.ILOO(L))THEN KFILO(K)=KFILOO(L) GOTO 4 ENDIF ENDDO WRITE(6,601)K,ISO(K),ILO(K) STOP 4 GOTO 3 ENDIF KO=K-1 ENDIF C C PRINT CASES SELECTED C WRITE(6,610) 610 FORMAT(/' EVEN'/' K IS IL/2J') DO K=1,KE WRITE(6,612)K,ISE(K),ILE(K) ENDDO WRITE(6,611) 611 FORMAT(/' ODD'/' K IS IL/2J') DO K=1,KO WRITE(6,612)K,ISO(K),ILO(K) ENDDO 612 FORMAT(3I5) C C READ REQUIRED ENERGY RANGE C M11=1 M21=1 M12=MXE M22=MXE NTEMP=0 TOP=.TRUE. C READ(5,MESH) C IF(NTEMP.GT.9)THEN WRITE(6,*)'*** CURRENTLY, CAN ONLY PROCESS 9 TEMPS: NTEMP=' X ,NTEMP STOP '***TOO MANY TEMPS SPECIFIED' ENDIF ALL=NTEMP.GT.0 C c m12=min(m12,mxe) c m22=min(m22,mxe) C IF(M11.GT.M12.OR.M12.GT.MXE)THEN WRITE(6,605)M11,M12,MXE STOP '***ERROR ON ENERGY MESH SPECIFICATION' ELSEIF(M21.GT.M22.OR.M22.GT.MXE)THEN WRITE(6,606)M21,M22,MXE STOP '***ERROR ON ENERGY MESH SPECIFICATION' ENDIF C C INITIALISE TOT C IF(ALL)THEN MM=MAX(M12,M22) MM2=2*MM MTM=(MM*(MM+1))/2 IF(MTM.GT.MZMSH)THEN WRITE(6,*)'***DIMENSION ERROR, INCREASE MZMSH TO:',MTM STOP '***DIMENSION ERROR, INCREASE MZMSH' ENDIF DO M=1,MTM TOT(M)=0. ENDDO ENDIF C C READ INDEX FOR D FILES C CALL RD00(ISE,ILE,KE,ISO,ILO,KO,KK1,KK2,ND,KPP,KASE) C WRITE(6,602) DO K=1,KASE K1=KK1(K) K2=KK2(K) IF(K1*K2.GT.0)WRITE(6,603)K,KK1(K),KK2(K),ND(K),KPP(K) ENDDO C C *** START LOOP ON CASES C DO K=1,KASE K1=KK1(K) K2=KK2(K) C IF(K1*K2.EQ.0)GO TO 1000 KT1=K1 KT2=K2 C OPEN(70,FILE='BF.'//NUM(ISE(K1))//NUM(ILE(K1))//'_'// + NUM(ISO(K2))//NUM(ILO(K2)),STATUS='UNKNOWN') C C GET DIPOLE DATA C IF(KPP(K))THEN CALL DKKP(ND(K)) ELSE CALL DKPK(ND(K)) ENDIF C CALL DSCALE IF(NPTYD2.NE.1)THEN PRINT*,' DKKP GIVES NPTYD2=',NPTYD2 STOP ENDIF IF(ISE(K1).NE.NSPND)THEN PRINT*,' K=',K,', ISE(K)=',ISE(K),', NSPND=',NSPND STOP ENDIF IF(LRGLD1.NE.ILE(K1).OR.LRGLD2.NE.ILO(K2))THEN WRITE(6,608)K,K1,K2,LRGLD1,ILE(K1),LRGLD2,ILO(K2) STOP ENDIF C C READ E-INDEPENDENT DATA FOR EVEN PARITY C CALL REIND1(KFILE(K1),ISE(K1),ILE(K1)) C C *START LOOP ON INITIAL ENERGIES C DO M1=1,M12 !nrb E1=EMESH(M1) C C READ E-DEPENDENT DATA FOR EVEN C CALL REDEP1(E1,M1,M11) IF(M1.LT.M11.OR.NCHOP1.EQ.0)GOTO 20 C CALL IJWBK1(E1,NCHOP1) C DO I=1,NCHF1 CALL DIFF(E1-ECH1(I),CC1(I),H,RZERO,SE(1,I),SE(2,I), + CE(1,I),CE(2,I),SEP(I),CEP(I)) ENDDO C C READ E-INDEPENDENT DATA FOR ODD C CALL REIND2(KFILO(K2),ISO(K2),ILO(K2)) C C CHECK AMAT AND BMATL C if(m1.eq.m11)CALL AB(ISE(K1),ILE(K1),ISO(K2),ILO(K2),K) C C GET DV CALL DVECT(E1) C C *START LOOP ON FINAL ENERGIES C DO M2=1,M22 !nrb E2=EMESH(M2) C C READ E-DEPENDENT DATA FOR EVEN C CALL REDEP2(E2,M1,M2,M21) IF(M1.EQ.M2.OR.M2.LT.M21.OR.NCHOP2.EQ.0)GOTO 10 C CALL IJWBK2(E2,NCHOP2) DO I=1,NCHF2 CALL DIFF(E2-ECH2(I),CC2(I),H,RZERO,SO(1,I),SO(2,I), + CO(1,I),CO(2,I),SOP(I),COP(I)) ENDDO C C INNER REGION CONTRIBUTION C CALL DIN(E1,E2) C C OUTER REGION CONTRIBUTION C CALL DOUT(E1,E2) C C SUM INNER AND OUTER C CALL DTOT(DR,DI,E1,E2) C C GET X=LINE-STRENGTH FOR SLPI_SLPI' COMBINATION !mjs C (Z-SCALED AND WITH ORIGINAL NORMALISATION) !mjs X=0. DO J2=1,NCHOP2 DO J1=1,NCHOP1 X=X+DR(J1,J2)**2+DI(J1,J2)**2 ENDDO ENDDO Q=FACT*X/(AZ4*ABS(E2-E1)**3) C C CASE OF NOT ALL IF(.NOT.ALL)THEN C WRITE LINE_STRENGTH TO 70 !mjs WRITE(70,700)E1,E2,X !mjs ELSE C C CASE OF ALL IF(M1.LT.M2)THEN N1=((M1-1)*(MM2-M1))/2+M2 STA(N1)=Q ELSE N2=((M2-1)*(MM2-M2))/2+M1 STB(N2)=Q ENDIF ENDIF C 10 CONTINUE ENDDO ! Ends M2 LOOP !mjs 20 CONTINUE ENDDO ! Ends M1 LOOP !mjs CLOSE(1) CLOSE(2) C C ADD TO TOT C IF(ALL)THEN DO M=1,MTM TOT(M)=TOT(M)+STA(M)+STB(M) ENDDO ENDIF 1000 ENDDO ! Ends Cases LOOP C C *** END LOOP ON CASES C C INCLUDE TOP UP C IF(ALL)THEN C IF(TOP)THEN IF(IPRINT.GT.0)THEN DO M=1,MTM STB(M)=TOT(M) ENDDO ENDIF C CALL TOPUP(ise(kt1),ILE(KT1),ILO(KT2),MM) C IF(IPRINT.GT.0)THEN WRITE(6,701) DO MD=1,MM-1 DO M1=1,MM-MD M2=M1+MD E1=EMESH(M1)*AZ2*0.5 E2=EMESH(M2)*AZ2*0.5 N1=((M1-1)*(MM2-M1))/2+M2 WRITE(6,703)M1,E1,M2,E2,STB(N1),TOT(N1) ENDDO ENDDO ENDIF ELSE IF(IPRINT.GT.0)THEN WRITE(6,702) DO MD=1,MM-1 DO M1=1,MM-MD M2=M1+MD E1=EMESH(M1)*AZ2*0.5 E2=EMESH(M2)*AZ2*0.5 N1=((M1-1)*(MM2-M1))/2+M2 WRITE(6,703)M1,E1,M2,E2,TOT(N1) ENDDO ENDDO ENDIF WRITE(6,*)' ' WRITE(6,*)'****WARNING: NO TOP-UP OF THERMAL AVERAGE XSCTNS' ENDIF C C LOOP OVER TEMPERATURES C DO LN=1,NTEMP READ(5,*,END=30)TEMP !READ TEMP OPEN(80,FILE='FF.'//NUM(LN),STATUS='UNKNOWN') WRITE(80,801)TEMP C CALL THAV(TEMP,MM,FT,MDM) C UXPLUS=PTARG(TEMP) F=TFACT/(PTARG(TEMP)*SQRT(TEMP)) G=.5*AZ2*(EMESH(2)-EMESH(1)) H=SQRT(.5*TEMP/157894.)/FACTG C DO MD=1,MDM SIGMA(MD)=F*FT(MD) OMEGA(MD)=G*MD GAUNT(MD)=GFACT*OMEGA(MD)**3*FT(MD)/(AZ2*UXPLUS) ENDDO C WRITE(80,800)(OMEGA(MD),SIGMA(MD),SIGMA(MD)*OMEGA(MD)**3, + GAUNT(MD),MD=1,MDM) ENDDO C WRITE(6,620) C ENDIF C 30 CONTINUE C C SUN DUM=DTIME(TARRY) TIME=TARRY(1) C CRAY CRAY CALL SECOND(TIME) C TIME=TIME/60.0 WRITE(6,999) TIME 999 FORMAT(//1X,'CPU TIME=',F9.3,' MIN',5X) STOP 'STGFF: NORMAL END' C 200 FORMAT(40X,'STGFF'/40X,'*****'//) 500 FORMAT(L1) 600 FORMAT(' K=',I2,', ISE=',I2,', ILE=',I2, + ' NOT FOUND ON F00 FILE') 601 FORMAT(' K=',I2,', ISO=',I2,', ILO=',I2, + ' NOT FOUND ON F00 FILE') 602 FORMAT(/' TRANSITIONS:'/' K KK1 KK2 ND KPP') 603 FORMAT(4I5,L5) 605 FORMAT(' ERROR IN ENERGY-SELECTION, M11,M12,MXE=',3I5) 606 FORMAT(' ERROR IN ENERGY-SELECTION, M21,M22,MXE=',3I5) 608 FORMAT(' MISMATCH, K=',I2,', K1=',I2,', K2=',I2/ + ' LRGLD1=',I2,', ILE=',I2/ + ' LRGLD2=',I2,', ILO=',I2) 609 FORMAT(' MD=',I5,', FJ=',1P,E12.3) 620 FORMAT(//12X,'THERMALLY AVERAGED XSCTNS WRITTEN TO FILE'/) 700 format(1p,8e11.3) 701 FORMAT(//5X,'M1',9X,'E1',7X,'M2',9X,'E2',7X,'SLINE',6X, X 'INC TOP-UP') 702 FORMAT(//5X,'M1',9X,'E1',7X,'M2',9X,'E2',7X,'SLINE') 703 FORMAT(I7,1P,E13.3,I7,3E13.3) 800 format(1p,4e12.3) 801 FORMAT(1P,E12.3) C END C*********************************************************************** SUBROUTINE AB(ISTOT1,ILTOT1,ISTOT2,ILTOT2,KKK) C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) C PARAMETER (TZERO=0.0) PARAMETER (EPS=1.D-14) C COMMON/CTARG/ENAT(MZTAR),NZED,NELC,NAST,NCONAT(MZTAR), + IST(MZTAR),ILT(MZTAR),IPT(MZTAR) COMMON/CHAN1/ECH1(MZCHF),CC1(MZCHF),LLCH1(MZCHF), X L2P1(MZCHF),KJ1(MZCHF), + NCHF1 COMMON/CHAN2/ECH2(MZCHF),CC2(MZCHF),LLCH2(MZCHF), X L2P2(MZCHF),KJ2(MZCHF), + NCHF2 COMMON/CDKKP/MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, 1 MNP2D2,NCHND2,LRGLD2,NPTYD2,NPAD, 2 DKL(MNPEXT,MNPEXT), CV X DKV(MNPEXT,MNPEXT), 3 AMAT(MZCHF,MZCHF),BMATL(MZCHF,MZCHF) CV X,BMATV(MZCHF,MZCHF) C DIMENSION IT1(MZTAR),IT2(MZTAR),DSTORE(MZTAR,MZTAR) C SAVE DSTORE DATA INIT/0/ C IF(INIT.EQ.0)THEN DO J=1,NAST DO I=1,NAST DSTORE(I,J)=TZERO ENDDO ENDDO INIT=1 ENDIF C DO I1=1,NCHF1 DO K=1,NAST IF(ABS( ECH1(I1)-ENAT(K) ) .LE. 1.E-6*ABS(ENAT(K)) )THEN IT1(I1)=K GOTO 1 ENDIF ENDDO PRINT*,' AB: NOT FOUND, I1=',I1 STOP 1 ENDDO C DO I2=1,NCHF2 DO K=1,NAST IF(ABS(ECH2(I2)-ENAT(K)).LE.1.E-6*ABS(ENAT(K)))THEN IT2(I2)=K GOTO 2 ENDIF ENDDO PRINT*,' AB: NOT FOUND, I2=',I2 STOP 2 ENDDO c c write(6,*)'*******',istot1,iltot1,istot2,iltot2,kkk C K12=0 DO I2=1,NCHF2 ITARG2=IT2(I2) LTARG2=ILT(ITARG2) ISTRG2=IST(ITARG2) LL2=LLCH2(I2) DO I1=1,NCHF1 ITARG1=IT1(I1) LTARG1=ILT(ITARG1) ISTRG1=IST(ITARG1) LL1=LLCH1(I1) IF(NSPND.EQ.0)K12=IABS(KJ1(I1)-KJ2(I2)) IF(ITARG1.EQ.ITARG2.AND.ABS(LL1-LL2).EQ.1.AND.K12.LE.2)THEN FGT=MAX(LL1,LL2) IF(NSPND.NE.0)THEN !LS A=(-1)**(LL1+ILTOT2-1-LTARG2)* + SQRT(DBLE(ISTOT1*(2*ILTOT1+1)*(2*ILTOT2+1))*FGT) + *W1(LL1,LL2,ILTOT1,ILTOT2,LTARG2) ELSE !JK NRB A=(-1)**((KJ2(I2)-LTARG2+2*LL1-2)/2)* X (-1)**((-ILTOT1+KJ2(I2)-3)/2)* X (-1)**((ILTOT1-ILTOT2)/2)* !O-E->E-O NRB X SQRT(DBLE((KJ1(I1)+1)*(KJ2(I2)+1)*(ILTOT1+1)*(ILTOT2+1))) X *(W2(2*LL1,2*LL2,KJ1(I1),KJ2(I2),LTARG1)* X W2(KJ1(I1),KJ2(I2),ILTOT1,ILTOT2,1))*SQRT(FGT) ENDIF c a=sign(a,amat(i1,i2)) !temp, nrb IF(ABS(AMAT(I1,I2)-A).GT.1.E-6*ABS(AMAT(I1,I2)))THEN WRITE(6,600)I1,I2,AMAT(I1,I2),A c write(6,*)'p2',2*ll1,2*ll2,kj1(i1),kj2(i2),ltarg1 c write(6,*)'p1',kj1(i1),kj2(i2),iltot1,iltot2,init STOP '*** ERROR SR.AB: MIS-MATCH ON AMAT' ENDIF ENDIF ENDDO ENDDO C C DO I2=1,NCHF2 ITARG2=IT2(I2) LTARG2=ILT(ITARG2) ISTRG2=IST(ITARG2) LL2=LLCH2(I2) DO I1=1,NCHF1 ITARG1=IT1(I1) LTARG1=ILT(ITARG1) ISTRG1=IST(ITARG1) LL1=LLCH1(I1) IF(ABS(BMATL(I1,I2)).GT.EPS)THEN IF(NSPND.NE.0)THEN !LS Q=(-1)**(-ILTOT1+LTARG2+LL2-1)* + SQRT(DBLE(ISTOT1*(2*ILTOT1+1)*(2*ILTOT2+1)) + /DBLE(ISTRG1))*W1(LTARG1,LTARG2,ILTOT1,ILTOT2,LL2) ELSE !JK NRB Q=(-1)**((-KJ1(I1)+LTARG2+2*LL2-2)/2)* X (-1)**((-ILTOT1+KJ2(I2)-3)/2)* X (-1)**((ILTOT1-ILTOT2)/2)* !O-E->E-O NRB X SQRT(DBLE((KJ1(I1)+1)*(KJ2(I2)+1)*(ILTOT1+1)*(ILTOT2+1))) X *W2(LTARG1,LTARG2,KJ1(I1),KJ2(I2),2*LL2)* X W2(KJ1(I1),KJ2(I2),ILTOT1,ILTOT2,1) ENDIF DD=BMATL(I1,I2)/Q IF(DSTORE(ITARG1,ITARG2).EQ.TZERO)THEN DSTORE(ITARG1,ITARG2)=DD ELSE c dd=sign(dd,dstore(itarg1,itarg2)) !temp IF(ABS(DSTORE(ITARG1,ITARG2)-DD).GT.1.E-6*ABS(DD))THEN WRITE(6,604)ITARG1,ITARG2,I1,I2,DSTORE(ITARG1,ITARG2),DD c write(6,*)'p2',ltarg1,ltarg2,kj1(i1),kj2(i2),2*ll2 c write(6,*)'p1',kj1(i1),kj2(i2),iltot1,iltot2,init STOP '*** ERROR SR.AB: MIS-MATCH ON BMAT' ENDIF ENDIF ENDIF ENDDO ENDDO C 600 FORMAT(' SUBROUTINE AB: I1=',I3,', I2=',I3,', AMAT=', + 1P,E15.6,', A=',E15.6) 604 FORMAT(' SUBROUTINE AB: IT1=',I3,', IT2=',I3,', I1=',I3, +', I2=',I3,', DSTORE=', 1P,E15.6,', DD=',E15.6) C RETURN END C*********************************************************************** SUBROUTINE ACSUB C IMPLICIT REAL*8(A-H,O-Z) C COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,TINY,DELTA,LACC C ACNUM=(24.*AC)**.166666667 ACJWBK=(6.*AC)**.2 ACZP=16.*AC LACC=0 IF(AC.LT.1.E-3)LACC=2 IF(AC.LT.1.E-4)LACC=4 C RETURN END C*********************************************************************** C REAL*8 FUNCTION ARGC(E,L,AC) C C CALCULATES ARG(GAMMA(L+1-I/K)) -1/K -(1/K)*LN(K) - L*PI/2 C NRB: NOT ACTUALLY CALLED BY INJWBK CASE TZED=0. C IMPLICIT REAL*8 (A-H,O-Z) C COMMON/NRBZED/TZED C IF(TZED.EQ.0.)THEN ARGC=-DBLE(L)*1.570796327 RETURN ENDIF C IF(E.GT.0)GOTO 10 ARGC=-(DBLE(L)+.25)*3.141592654 RETURN C 10 FK=SQRT(E) ET=1./FK IP=L+1 P=IP PP=IP*IP C IF(AC.LT.1.E-4)GOTO 100 A1=10.*SQRT(ET)-ET*ET IF(A1.GT.PP)GOTO 20 X=PP*E XP1=X+1. XH=P*FK A=-1.570796327*(P+DBLE(L)-.5) GOTO 200 20 L1=IP IP=1.+SQRT(A1) P=IP PP=IP*IP X=PP*E XP1=X+1. XH=P*FK A=-1.570796327*(P+DBLE(L)-.5) L2=IP-1 DO I=L1,L2 A=A+ATAN(ET/DBLE(I)) ENDDO GOTO 200 C 100 A1=35.*ET**.25-ET*ET IF(A1.GT.PP)GOTO 120 X=PP*E XP1=X+1. XH=P*FK A=-1.570796327*(P+DBLE(L)-.5) GOTO 140 120 L1=IP IP=1.+SQRT(A1) P=IP PP=IP*IP X=PP*E XP1=X+1. XH=P*FK A=-1.570796327*(P+DBLE(L)-.5) L2=IP-1 DO I=L1,L2 A=A+ATAN(ET/DBLE(I)) ENDDO 140 A=A+.000396825540*FK*E*(7.*(1.-3.*X)*XP1*XP1+ C 2.*E*(1.-10.*X+5.*X*X))*XP1**(-5) C 200 A1=FK*X*X*.1667*PP IF(A1.GT.AC)GOTO 210 A=A-FK*(2.-X)*.25*PP GOTO 220 210 A=A-.5*ET*LOG(XP1) 220 A2=(P-.5)*XH A1=A2*X*X IF(A1.GT.AC)GOTO 230 A=A+A2*(1.-X*.33333333) GOTO 240 230 A=A+(P-.5)*ATAN(XH) 240 ARGC=A+FK/(12.*(1.+X)) C RETURN END C C*********************************************************************** C BLOCK DATA C C PROVIDES C BLOCK DATA C C DATA FOR QUADRATURES - C LAGUERRE AND LEGENDRE QUADRATURES WITH NUMBERS OF POINTS C N = 2, 4, 6, 8 AND 10 C IMPLICIT REAL*8(A-H,O-Z) C INCLUDE 'PARAM' C COMMON/CBLK/XLAG(30),WLAG(30),XLEG(15),WLEG(15) COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,TINY,DELTA,LACC C DATA TINY,DELTA/1.E-6,0.1/ !!! was delta=0.1 !!!!! C DATA XLAG/ 1 .58578644,3.4142136, 2 .32254769,1.7457611,4.5366203,9.3950709, 3 .22284660,1.1889321,2.9927363,5.7751436,9.8374674, 4 15.982874, 5 .17027963,.90370178,2.2510866,4.26670017,7.0459054, 6 10.758516,15.7406786,22.8631317, 7 .13779347,.72945455,1.8083429,3.4014337,5.5524961, 8 8.3301527,11.8437858,16.279258,21.996586,29.920697/ DATA WLAG/ 1 .85355339,.14644661, 2 .60315410,.35741869,.38887909E-1,.53929471E-3, 3 .45896467,.41700083,.11337338,.10399197E-1, 4 .26101720E-3,.89854791E-6, 5 .36918859,.41878678,.17579499,3.3343492E-2,2.7945362E-3, 6 9.0765088E-5,8.4857467E-7,1.0480012E-9, 7 .30844112,.40111993,.21806829,6.2087456E-2,9.5015170E-3, 8 7.5300839E-4,2.8259233E-5,4.2493140E-7,1.8395648E-9, 9 9.9118272E-13/ DATA XLEG,WLEG/.577350269, 1 .339981044,.861136312, 2 .238619186,.661209386,.932469514, 3 .183434642,.525532410,.796666477,.960289856, 4 .148874339,.433395394,.679409568,.865063367,.973906529, 5 1., 6 .652145159,.347854845, 7 .467913935,.360761573,.171324492, 8 .362683783,.313706646,.222381034,.101228536, 9 .295524225,.269266719,.219086363,.149451349,.066671344/ C END C C*********************************************************************** SUBROUTINE BODE C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C COMMON/CPOINT/RZERO,RTWO,H,KP2 COMMON/CBODE/BDM2(MZPTS),BD0(MZPTS),BD1(MZPTS) C PARAMETER(B1=14./45.,B2=64./45.,B3=24./45.,B4=28./45.) C RM2(I)=H/(RZERO+H*(I-1))**2 R(I)=H*(RZERO+H*(I-1)) C BDM2(1)=RM2(1)*B1 BDM2(KP2)=RM2(KP2)*B1 BD0(1)=H*B1 BD0(KP2)=H*B1 BD1(1)=R(1)*B1 BD1(KP2)=R(KP2)*B1 C DO I=2,KP2,2 BDM2(I)=RM2(I)*B2 BD0(I)=H*B2 BD1(I)=R(I)*B2 ENDDO DO I=3,KP2,4 BDM2(I)=RM2(I)*B3 BD0(I)=H*B3 BD1(I)=R(I)*B3 ENDDO DO I=5,KP2-1,4 BDM2(I)=RM2(I)*B4 BD0(I)=H*B4 BD1(I)=R(I)*B4 ENDDO C RETURN END C C*********************************************************************** SUBROUTINE DIFF(E,C,H,R1,F1,F2,G1,G2,F1P,G1P) C IMPLICIT REAL*8 (A-H,O-Z) C PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (THREE=3.0) PARAMETER (FOUR=4.0) PARAMETER (SIX=6.0) PARAMETER (SEVEN=7.0) PARAMETER (TWELVE=12.0) PARAMETER (P1=ONE/30.) PARAMETER (P2=ONE/40.) PARAMETER (P3=SEVEN/15.) PARAMETER (P4=TWO/15.) PARAMETER (P5=ONE/360.) PARAMETER (P6=ONE/20.) PARAMETER (P7=ONE/120.) C COMMON/NRBZED/TZED C NRB: C NEUTRAL CASE ADDED C V(X)=EQ+X*(Q2-X*CQ) C Q=H*H EQ=E*Q Q2=TWO*Q CQ=C*Q C X1=ONE/R1 CX=C*X1 HX=H*X1 A1=-TWO*HX*HX*H U1P=A1*(ONE*TZED-CX) A1=-A1*HX U1PP=A1*(TWO*TZED-THREE*CX) A1=-A1*HX*SIX U1PPP=A1*(ONE*TZED-TWO*CX) U1=V(X1) C R2=R1+H X2=ONE/R2 U2=V(X2) C A2=ONE+U2*P1 B1=ONE+U1*(P2*U1-P3)-P4*U1P-P2*U1PP 1 +P5*(FOUR*U1*U1P-U1PPP) C1=H*(ONE+U1*(P5*U1-P4) 1 -P6*U1P-P7*U1PP) C F1P=(A2*F2-B1*F1)/C1 G1P=(A2*G2-B1*G1)/C1 C RETURN END C C*********************************************************************** SUBROUTINE DIN(E1,E2) C IMPLICIT REAL*8 (A-H,O-Z) C LOGICAL QDT1,QDT2 C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) PARAMETER (LWORK=MZCHF*MZCHF) PARAMETER (MWORK=MZMNP) C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) C COMMON/CBUT/IBUT,K10,K20 COMMON/CDV/DV(MNPEXT,MZCHF) COMMON/CHAN2/ECH2(MZCHF),CC2(MZCHF),LLCH2(MZCHF), X L2P2(MZCHF),KJ2(MZCHF), + NCHF2 COMMON/CRMAT2/VALUE2(MZMNP),WMAT2(MZMNP,MZCHF),MNP2 COMMON/DEP1/RK1(MZCHF,MZCHF),SE(MZPTS,MZCHF),CE(MZPTS,MZCHF), + NCHOP1,QDT1,F1(MZCHF,MZCHF),FP1(MZCHF,MZCHF),FKNU1(MZCHF) COMMON/DEP2/RK2(MZCHF,MZCHF),SO(MZPTS,MZCHF),CO(MZPTS,MZCHF), + NCHOP2,QDT2,F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FKNU2(MZCHF) COMMON/CD/DN(MZCHF,MZCHF),DA(MZCHF,MZCHF),STA(MZMSH),STB(MZMSH), + TOT(MZMSH) COMMON/SCALE/AZ,AZP,AZ2,AZ4,AZP2 COMMON/NRBWRK/WORK(LWORK),TEMP1(MWORK),TEMP2(MWORK) C DIMENSION V2(MNPEXT,MZCHF) C IF(E1.EQ.E2)THEN DO I2=1,NCHOP2 DO I1=1,NCHOP1 DN(I1,I2)=TZERO ENDDO ENDDO RETURN ENDIF C DO K2=1,MNP2 TEMP2(K2)=ONE/(VALUE2(K2)-E2) ENDDO DO J2=1,NCHOP2 DO K2=1,MNP2 V2(K2,J2)=TZERO ENDDO DO I=1,NCHF2 DO K2=1,MNP2 V2(K2,J2)=V2(K2,J2)+WMAT2(K2,I)*FP2(I,J2) ENDDO ENDDO DO K2=1,MNP2 V2(K2,J2)=V2(K2,J2)*TEMP2(K2) ENDDO C DO I=1,NCHF2 V2(MNP2+I,J2)=AZ*FP2(I,J2) ENDDO DO J1=1,NCHOP1 DN(J1,J2)=TZERO DO K2=1,K20 DN(J1,J2)=DN(J1,J2)+DV(K2,J1)*V2(K2,J2) ENDDO ENDDO ENDDO C C PUT FACTOR (E2-E1)**2 IN DN C DO J=1,NCHOP2 DO I=1,NCHOP1 DN(I,J)=DN(I,J)*(E2-E1)**2 ENDDO ENDDO C RETURN END C C*********************************************************************** SUBROUTINE DKKP(IND) C IMPLICIT REAL*8 (A-H,O-Z) C CHARACTER*1 NUM(0:9) C C ADAPTED FROM C PROGRAM OF KTT FOR READING DIPOLE MATRIX ELEMENT DATA C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) C COMMON/CDKKP/MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, 1 MNP2D2,NCHND2,LRGLD2,NPTYD2,NPAD, 2 DKL(MNPEXT,MNPEXT), CV X DKV(MNPEXT,MNPEXT), 3 AMAT(MZCHF,MZCHF),BMATL(MZCHF,MZCHF) CV X,BMATV(MZCHF,MZCHF) COMMON/CBUT/IBUT,K10,K20 C CC DIMENSION CGC(MZCHF) C DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ C C OPEN(4,FILE='D'//NUM(IND/10)//NUM(IND-10*(IND/10)), +FORM='UNFORMATTED',STATUS='OLD') READ(4) NOTERM,MNP2D2,NCHND2,LRGLD2,NPTYD2,NSPND, 1 MNP2D1,NCHND1,LRGLD1 C C CALCULATE NPTYD1 NPTYD1=1-NPTYD2+2*(NPTYD2/2) C C IAIN1 = (MNP2D1 - 1) / NOTERM IBIN1 = ( MNP2D2 - 1) / NOTERM MCI = 0 DO 800 IK = 1 , IAIN1 MCH = MCI + 1 MCI = MCI + NOTERM NCI = 0 DO 700 JK = 1 , IBIN1 NCH = NCI + 1 NCI = NCI + NOTERM READ(4) ((DKL(I,J),J=NCH,NCI),I=MCH,MCI) CV 1 ,((DKV(I,J),J=NCH,NCI),I=MCH,MCI) 700 CONTINUE NCH = NCI + 1 NCI = MNP2D2 READ(4) ((DKL(I,J),J=NCH,NCI),I=MCH,MCI) CV 1 ,((DKV(I,J),J=NCH,NCI),I=MCH,MCI) 800 CONTINUE MCH = MCI + 1 MCI = MNP2D1 NCI = 0 DO 900 JK = 1 , IBIN1 NCH = NCI + 1 NCI = NCI + NOTERM READ(4) ((DKL(I,J),J=NCH,NCI),I=MCH,MCI) CV 1 ,((DKV(I,J),J=NCH,NCI),I=MCH,MCI) 900 CONTINUE NCH = NCI + 1 NCI = MNP2D2 READ(4) ((DKL(I,J),J=NCH,NCI),I=MCH,MCI) CV 1 ,((DKV(I,J),J=NCH,NCI),I=MCH,MCI) C C READ BUTTLE PART C IF(IBUT.EQ.-1) GOTO 1000 NCH=1 NCI=MNP2D2 MCH=MNP2D1+1 MCI=MNP2D1+NCHND1 READ(4)((DKL(I,J),I=MCH,MCI),J=NCH,NCI) CV 1 ,((DKV(I,J),I=MCH,MCI),J=NCH,NCI) NCH=MNP2D2+1 NCI=MNP2D2+NCHND2 MCH=1 MCI=MNP2D1 READ(4)((DKL(I,J),J=NCH,NCI),I=MCH,MCI) CV 1 ,((DKV(I,J),J=NCH,NCI),I=MCH,MCI) MCH=MNP2D1+1 MCI=MNP2D1+NCHND1 READ(4)((DKL(I,J),J=NCH,NCI),I=MCH,MCI) CV 1 ,((DKV(I,J),J=NCH,NCI),I=MCH,MCI) C C-----READ ANGULAR COEFFICIENTS 1000 CONTINUE READ(4) MAXM1 CC ,(CGC(M),M=1,MAXM1) READ(4) ((AMAT(J,I),J=1,NCHND1),I=1,NCHND2), 1 ((BMATL(J,I),J=1,NCHND1),I=1,NCHND2) CV 2 ,((BMATV(J,I),J=1,NCHND1),I=1,NCHND2) C CLOSE(4) C RETURN END C C*********************************************************************** C SUBROUTINE DKPK(IND) C IMPLICIT REAL*8 (A-H,O-Z) C CHARACTER*1 NUM(0:9) C C ADAPTED FROM C PROGRAM OF KTT FOR READING DIPOLE MATRIX ELEMENT DATA C AND MODIFIED TO TAKE IN SLP PAIR IN REVERSE ORDER C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) C COMMON/CDKKP/MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, 1 MNP2D2,NCHND2,LRGLD2,NPTYD2,NPAD, 2 DKL(MNPEXT,MNPEXT), CV X DKV(MNPEXT,MNPEXT), 3 AMAT(MZCHF,MZCHF),BMATL(MZCHF,MZCHF) CV X,BMATV(MZCHF,MZCHF) COMMON/CBUT/IBUT,K10,K20 C CC DIMENSION CGC(MZCHF) C DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ C OPEN(4,FILE='D'//NUM(IND/10)//NUM(IND-10*(IND/10)), +FORM='UNFORMATTED',STATUS='OLD') C READ(4) NOTERM,MNP2D2,NCHND2,LRGLD2,NPTYD2,NSPND, 1 MNP2D1,NCHND1,LRGLD1 C C CALCULATE NPTYD1 NPTYD1=1-NPTYD2+2*(NPTYD2/2) C C IAIN1 = (MNP2D1 - 1) / NOTERM IBIN1 = ( MNP2D2 - 1) / NOTERM MCI = 0 DO 800 IK = 1 , IAIN1 MCH = MCI + 1 MCI = MCI + NOTERM NCI = 0 DO 700 JK = 1 , IBIN1 NCH = NCI + 1 NCI = NCI + NOTERM READ(4) ((DKL(J,I),J=NCH,NCI),I=MCH,MCI) CV 1 ,((DKV(J,I),J=NCH,NCI),I=MCH,MCI) 700 CONTINUE NCH = NCI + 1 NCI = MNP2D2 READ(4) ((DKL(J,I),J=NCH,NCI),I=MCH,MCI) CV 1 ,((DKV(J,I),J=NCH,NCI),I=MCH,MCI) 800 CONTINUE MCH = MCI + 1 MCI = MNP2D1 NCI = 0 DO 900 JK = 1 , IBIN1 NCH = NCI + 1 NCI = NCI + NOTERM READ(4) ((DKL(J,I),J=NCH,NCI),I=MCH,MCI) CV 1 ,((DKV(J,I),J=NCH,NCI),I=MCH,MCI) 900 CONTINUE NCH = NCI + 1 NCI = MNP2D2 READ(4) ((DKL(J,I),J=NCH,NCI),I=MCH,MCI) CV 1 ,((DKV(J,I),J=NCH,NCI),I=MCH,MCI) C C READ BUTTLE PART C IF(IBUT.EQ.-1) GOTO 1000 NCH=1 NCI=MNP2D2 MCH=MNP2D1+1 MCI=MNP2D1+NCHND1 READ(4)((DKL(J,I),I=MCH,MCI),J=NCH,NCI) CV 1 ,((DKV(J,I),I=MCH,MCI),J=NCH,NCI) NCH=MNP2D2+1 NCI=MNP2D2+NCHND2 MCH=1 MCI=MNP2D1 READ(4)((DKL(J,I),J=NCH,NCI),I=MCH,MCI) CV 1 ,((DKV(J,I),J=NCH,NCI),I=MCH,MCI) MCH=MNP2D1+1 MCI=MNP2D1+NCHND1 READ(4)((DKL(J,I),J=NCH,NCI),I=MCH,MCI) CV 1 ,((DKV(J,I),J=NCH,NCI),I=MCH,MCI) C C-----READ ANGULAR COEFFICIENTS 1000 READ(4) MAXM1 CC ,(CGC(M),M=1,MAXM1) READ(4) ((AMAT(I,J),J=1,NCHND1),I=1,NCHND2), 1 ((BMATL(I,J),J=1,NCHND1),I=1,NCHND2) CV 2 ,((BMATV(I,J),J=1,NCHND1),I=1,NCHND2) C C EXCHANGE PARAMETERS FOR THE TWO STATES C I=MNP2D1 MNP2D1=MNP2D2 MNP2D2=I I=NCHND1 NCHND1=NCHND2 NCHND2=I I=LRGLD1 LRGLD1=LRGLD2 LRGLD2=I I=NPTYD1 NPTYD1=NPTYD2 NPTYD2=I C CLOSE(4) C RETURN END C C*********************************************************************** SUBROUTINE DOUT(E1,E2) C IMPLICIT REAL*8 (A-H,O-Y) IMPLICIT COMPLEX*16(Z) C LOGICAL QDT1,QDT2 C C COMPUTES OUTER-REGIONS INTEGRALS SFOR FUNCTIONS S AND C. C INTEGRALS MULTIPLIED BY FACTOR (E2-E1)**2 C (GIVING FINITE VALUES FOR E1=E2) C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) C COMMON/CHAN1/ECH1(MZCHF),CC1(MZCHF),LLCH1(MZCHF), X L2P1(MZCHF),KJ1(MZCHF), + NCHF1 COMMON/CHAN2/ECH2(MZCHF),CC2(MZCHF),LLCH2(MZCHF), X L2P2(MZCHF),KJ2(MZCHF), + NCHF2 COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,TINY,DELTA,LACC COMMON/DEP1/RK1(MZCHF,MZCHF),SE(MZPTS,MZCHF),CE(MZPTS,MZCHF), + NCHOP1,QDT1,F1(MZCHF,MZCHF),FP1(MZCHF,MZCHF),FKNU1(MZCHF) COMMON/DEP2/RK2(MZCHF,MZCHF),SO(MZPTS,MZCHF),CO(MZPTS,MZCHF), + NCHOP2,QDT2,F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FKNU2(MZCHF) COMMON/DPP/SE1(MZCHF),SO1(MZCHF),CE1(MZCHF),CO1(MZCHF) COMMON/DIFF1/SEP(MZCHF),CEP(MZCHF) COMMON/DIFF2/SOP(MZCHF),COP(MZCHF) COMMON/CPOINT/RZERO,RTWO,H,KP2 COMMON/CBODE/BDM2(MZPTS),BD0(MZPTS),BD1(MZPTS) COMMON/CDKKP/MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, 1 MNP2D2,NCHND2,LRGLD2,NPTYD2,NPAD, 2 DKL(MNPEXT,MNPEXT), CV X DKV(MNPEXT,MNPEXT), 3 AMAT(MZCHF,MZCHF),BMATL(MZCHF,MZCHF) CV X,BMATV(MZCHF,MZCHF) C COMMON/CALP/ASS(MZCHF,MZCHF),ASC(MZCHF,MZCHF),ACS(MZCHF,MZCHF), + ACC(MZCHF,MZCHF) COMMON/CBET/BSS(MZCHF,MZCHF),BSC(MZCHF,MZCHF),BCS(MZCHF,MZCHF), + BCC(MZCHF,MZCHF) C DE2=(E2-E1)**2 C C INITIALISE ALL OUTER-REGION INTEGRALS DO J=1,NCHF2 DO I=1,NCHF1 ASS(I,J)=0. ASC(I,J)=0. ACS(I,J)=0. ACC(I,J)=0. BSS(I,J)=0. BSC(I,J)=0. BCS(I,J)=0. BCC(I,J)=0. ENDDO ENDDO C C AMAT INTEGRALS C ************** C C OPEN-OPEN DO J=1,NCHOP2 EPS2=E2-ECH2(J) C2=CC2(J) DO I=1,NCHOP1 IF(ABS(AMAT(I,J)).LE.TINY)GOTO 1 EPS1=E1-ECH1(I) C1=CC1(I) C c!!!!!! IF(ABS(FKNU1(I)-FKNU2(J)).GE.DELTA*1.e-10)THEN !!!!!!! 1.e-10 !!!! if(abs(fknu1(i)-fknu2(j)).gt.-1)then !!!!!!!!!!!!!!!!!! C USING ACC FORM C RZERO TO RTWO DO K=1,KP2 ASS(I,J)=ASS(I,J)+SE(K,I)*SO(K,J)*BDM2(K) ASC(I,J)=ASC(I,J)+SE(K,I)*CO(K,J)*BDM2(K) ACS(I,J)=ACS(I,J)+CE(K,I)*SO(K,J)*BDM2(K) ACC(I,J)=ACC(I,J)+CE(K,I)*CO(K,J)*BDM2(K) ENDDO C RWTO TO INFINITY NLAG=8 NLEG=10 !!!!! ZP=ZPLAG(I,J,NLAG,2,EPS1,EPS2) ALP=RTWO*(SQRT(EPS1)-SQRT(EPS2)) IF(ABS(ALP).GT.2)THEN ZQ=ZQLAG(I,J,NLAG,2,EPS1,EPS2) ELSE ZQ=ZQLEG(I,J,NLEG,2,EPS1,EPS2) ENDIF ASS(I,J)=ASS(I,J)+.5*DBLE(ZQ-ZP) ASC(I,J)=ASC(I,J)+.5*DIMAG(ZP+ZQ) ACS(I,J)=ACS(I,J)+.5*DIMAG(ZP-ZQ) ACC(I,J)=ACC(I,J)+.5*DBLE(ZP+ZQ) C ADD SURFACE TERM ASS(I,J)=ASS(I,J)*4.+SFCE(RZERO, + EPS1,C1,EPS2,C2,SE1(I),SEP(I),SO1(J),SOP(J)) ASC(I,J)=ASC(I,J)*4.+SFCE(RZERO, + EPS1,C1,EPS2,C2,SE1(I),SEP(I),CO1(J),COP(J)) ACS(I,J)=ACS(I,J)*4.+SFCE(RZERO, + EPS1,C1,EPS2,C2,CE1(I),CEP(I),SO1(J),SOP(J)) ACC(I,J)=ACC(I,J)*4.+SFCE(RZERO, + EPS1,C1,EPS2,C2,CE1(I),CEP(I),CO1(J),COP(J)) C ELSE C C USING LENGTH FORM DO K=1,KP2 ASS(I,J)=ASS(I,J)+SE(K,I)*SO(K,J)*BD1(K) ASC(I,J)=ASC(I,J)+SE(K,I)*CO(K,J)*BD1(K) ACS(I,J)=ACS(I,J)+CE(K,I)*SO(K,J)*BD1(K) ACC(I,J)=ACC(I,J)+CE(K,I)*CO(K,J)*BD1(K) ENDDO NLAG=8 !!! NLEG=NLAG ZP=ZPLAG(I,J,NLAG,-1,EPS1,EPS2) ALP=RTWO*(SQRT(EPS1)-SQRT(EPS2)) IF(ABS(ALP).GT.2)THEN ZQ=ZQLAG(I,J,NLAG,-1,EPS1,EPS2) ELSE ZQ=ZQLEG(I,J,NLEG,-1,EPS1,EPS2) ENDIF ASS(I,J)=ASS(I,J)+.5*DBLE(ZQ-ZP) ASC(I,J)=ASC(I,J)+.5*DIMAG(ZP+ZQ) ACS(I,J)=ACS(I,J)+.5*DIMAG(ZP-ZQ) ACC(I,J)=ACC(I,J)+.5*DBLE(ZP+ZQ) ENDIF C 1 CONTINUE ENDDO ENDDO C C OPEN-CLOSED DO J=NCHOP2+1,NCHF2 C2=CC2(J) EPS2=E2-ECH2(J) DO I=1,NCHOP1 IF(ABS(AMAT(I,J)).LE.TINY)GOTO 2 C1=CC1(I) EPS1=E1-ECH1(I) DO K=1,KP2 ASC(I,J)=ASC(I,J)+SE(K,I)*CO(K,J)*BDM2(K) ACC(I,J)=ACC(I,J)+CE(K,I)*CO(K,J)*BDM2(K) ENDDO CALL ZSLG12(I,J,NLAG,2,ZS3) ASC(I,J)=ASC(I,J)+DIMAG(ZS3) ACC(I,J)=ACC(I,J)+DBLE(ZS3) ASC(I,J)=ASC(I,J)*4.+SFCE(RZERO, + EPS1,C1,EPS2,C2,SE1(I),SEP(I),CO1(J),COP(J)) ACC(I,J)=ACC(I,J)*4.+SFCE(RZERO, + EPS1,C1,EPS2,C2,CE1(I),CEP(I),CO1(J),COP(J)) 2 CONTINUE ENDDO ENDDO C C CLOSED-OPEN DO J=1,NCHOP2 C2=CC2(J) EPS2=E2-ECH2(J) DO I=NCHOP1+1,NCHF1 IF(ABS(AMAT(I,J)).LE.TINY)GOTO 3 C1=CC1(I) EPS1=E1-ECH1(I) DO K=1,KP2 ACS(I,J)=ACS(I,J)+CE(K,I)*SO(K,J)*BDM2(K) ACC(I,J)=ACC(I,J)+CE(K,I)*CO(K,J)*BDM2(K) ENDDO CALL ZSLG21(I,J,NLAG,2,ZS3) ACS(I,J)=ACS(I,J)+DIMAG(ZS3) ACC(I,J)=ACC(I,J)+DBLE(ZS3) ACS(I,J)=ACS(I,J)*4.+SFCE(RZERO, + EPS1,C1,EPS2,C2,CE1(I),CEP(I),SO1(J),SOP(J)) ACC(I,J)=ACC(I,J)*4.+SFCE(RZERO, + EPS1,C1,EPS2,C2,CE1(I),CEP(I),CO1(J),COP(J)) 3 CONTINUE ENDDO ENDDO C C CLOSED-CLOSED DO J=NCHOP2+1,NCHF2 C2=CC2(J) EPS2=E2-ECH2(J) DO I=NCHOP1+1,NCHF1 IF(ABS(AMAT(I,J)).LE.TINY)GOTO 4 C1=CC1(I) EPS1=E1-ECH1(I) C C IF(ABS(FKNU1(I)-FKNU2(J)).GE.DELTA)THEN C USING ACC FORM DO K=1,KP2 ACC(I,J)=ACC(I,J)+CE(K,I)*CO(K,J)*BDM2(K) ENDDO CALL TLAG(I,J,NLAG,2,T1) ACC(I,J)=ACC(I,J)+T1 ACC(I,J)=ACC(I,J)*4.+SFCE(RZERO, + EPS1,C1,EPS2,C2,CE1(I),CEP(I),CO1(J),COP(J)) C C ELSE C (DELETED USE OF LENGTH FORM FOR CLOSED-CLOSED) !mjs C USING LENGTH FORM C ACC(I,J)=0. C DO K=1,KP2 C ACC(I,J)=ACC(I,J)+CE(K,I)*CO(K,J)*BD1(K) C ENDDO C CALL TLAG(I,J,NLAG,-1,T1) C ACC(I,J)=ACC(I,J)+T1 C ENDIF C 4 CONTINUE ENDDO ENDDO C PUT IN AMAT FACTORS DO J=1,NCHF2 DO I=1,NCHF1 ASS(I,J)=ASS(I,J)*AMAT(I,J) ASC(I,J)=ASC(I,J)*AMAT(I,J) ACS(I,J)=ACS(I,J)*AMAT(I,J) ACC(I,J)=ACC(I,J)*AMAT(I,J) ENDDO ENDDO C C ADD BMAT CONTRIBUTIONS C ********************** C C OPEN-OPEN DO J=1,NCHOP2 EPS2=E2-ECH2(J) DO I=1,NCHOP1 IF(ABS(BMATL(I,J)).GT.TINY)THEN EPS1=E1-ECH1(I) X=DE2/(EPS1-EPS2) BSS(I,J)=BSS(I,J)+BMATL(I,J)* + (SE1(I)*SOP(J)-SEP(I)*SO1(J))*X BSC(I,J)=BSC(I,J)+BMATL(I,J)* + (SE1(I)*COP(J)-SEP(I)*CO1(J))*X BCS(I,J)=BCS(I,J)+BMATL(I,J)* + (CE1(I)*SOP(J)-CEP(I)*SO1(J))*X BCC(I,J)=BCC(I,J)+BMATL(I,J)* + (CE1(I)*COP(J)-CEP(I)*CO1(J))*X ENDIF ENDDO ENDDO C C OPEN=CLOSED DO J=NCHOP2+1,NCHF2 EPS2=E2-ECH2(J) DO I=1,NCHOP1 IF(ABS(BMATL(I,J)).GT.TINY)THEN EPS1=E1-ECH1(I) X=DE2/(EPS1-EPS2) BSC(I,J)=BSC(I,J)+BMATL(I,J)* + (SE1(I)*COP(J)-SEP(I)*CO1(J))*X BCC(I,J)=BCC(I,J)+BMATL(I,J)* + (CE1(I)*COP(J)-CEP(I)*CO1(J))*X ENDIF ENDDO ENDDO C C CLOSED-OPEN DO J=1,NCHOP2 EPS2=E2-ECH2(J) DO I=NCHOP1+1,NCHF2 IF(ABS(BMATL(I,J)).GT.TINY)THEN EPS1=E1-ECH1(I) X=DE2/(EPS1-EPS2) BCS(I,J)=BCS(I,J)+BMATL(I,J)* + (CE1(I)*SOP(J)-CEP(I)*SO1(J))*X BCC(I,J)=BCC(I,J)+BMATL(I,J)* + (CE1(I)*COP(J)-CEP(I)*CO1(J))*X ENDIF ENDDO ENDDO C C CLOSED=CLOSED DO J=NCHOP2+1,NCHF2 EPS2=E2-ECH2(J) DO I=NCHOP1+1,NCHF1 IF(ABS(BMATL(I,J)).GT.TINY)THEN EPS1=E1-ECH1(I) X=DE2/(EPS1-EPS2) BCC(I,J)=BCC(I,J)+BMATL(I,J)* + (CE1(I)*COP(J)-CEP(I)*CO(1,J))*X ENDIF ENDDO ENDDO C RETURN END C C*********************************************************************** SUBROUTINE DSCALE C IMPLICIT REAL*8 (A-H,O-Z) C C SCALING OF DKK DATA C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) C COMMON/CDKKP/MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, 1 MNP2D2,NCHND2,LRGLD2,NPTYD2,NPAD, 2 DKL(MNPEXT,MNPEXT), CV X DKV(MNPEXT,MNPEXT), 3 AMAT(MZCHF,MZCHF),BMATL(MZCHF,MZCHF) CV X,BMATV(MZCHF,MZCHF) COMMON/SCALE/AZ,AZP,AZ2,AZ4,AZP2 COMMON/CBUT/IBUT,K10,K20 C K10=MNP2D1 IF(IBUT.GT.0)K10=K10+NCHND1 K20=MNP2D2 IF(IBUT.GT.0)K20=K20+NCHND2 DO K2=1,K20 DO K1=1,K10 DKL(K1,K2)=AZ*DKL(K1,K2) CV DKV(K1,K2)=AZP*DKV(K1,K2) ENDDO ENDDO C DO I2=1,NCHND2 DO I1=1,NCHND1 BMATL(I1,I2)=AZ*BMATL(I1,I2) CV BMATV(I1,I2)=AZP*BMATV(I1,I2) ENDDO ENDDO C C INCLUDE SPIN FACTOR IN ALL REQUIRED COEFFICIENTS C if(nspnd.eq.0)return C SQRS=SQRT(DBLE(NSPND)) DO K2=1,K20 DO K1=1,K10 DKL(K1,K2)=SQRS*DKL(K1,K2) ENDDO ENDDO DO I2=1,NCHND2 DO I1=1,NCHND1 AMAT(I1,I2)=SQRS*AMAT(I1,I2) BMATL(I1,I2)=SQRS*BMATL(I1,I2) ENDDO ENDDO C RETURN END C C*********************************************************************** SUBROUTINE DTOT(DR,DI, e1,e2) C IMPLICIT REAL*8 (A-H,O-Z) C LOGICAL QDT1,QDT2 C INCLUDE 'PARAM' C PARAMETER (LWORK=MZCHF*MZCHF) PARAMETER (MWORK=MZMNP) PARAMETER (ONE=1.0) PARAMETER (FOUR=4.0) PARAMETER (TZERO=0.0) C COMMON/CD/DN(MZCHF,MZCHF),DA(MZCHF,MZCHF),STA(MZMSH),STB(MZMSH), + TOT(MZMSH) COMMON/CALP/ASS(MZCHF,MZCHF),ASC(MZCHF,MZCHF),ACS(MZCHF,MZCHF), + ACC(MZCHF,MZCHF) COMMON/CBET/BSS(MZCHF,MZCHF),BSC(MZCHF,MZCHF),BCS(MZCHF,MZCHF), + BCC(MZCHF,MZCHF) COMMON/CHAN1/ECH1(MZCHF),CC1(MZCHF),LLCH1(MZCHF), X L2P1(MZCHF),KJ1(MZCHF), + NCHF1 COMMON/CHAN2/ECH2(MZCHF),CC2(MZCHF),LLCH2(MZCHF), X L2P2(MZCHF),KJ2(MZCHF), + NCHF2 COMMON/DEP1/RK1(MZCHF,MZCHF),SE(MZPTS,MZCHF),CE(MZPTS,MZCHF), + NCHOP1,QDT1,F1(MZCHF,MZCHF),FP1(MZCHF,MZCHF),FKNU1(MZCHF) COMMON/DEP2/RK2(MZCHF,MZCHF),SO(MZPTS,MZCHF),CO(MZPTS,MZCHF), + NCHOP2,QDT2,F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FKNU2(MZCHF) COMMON/NRBWRK/WORK(LWORK),TEMP1(MWORK),TEMP2(MWORK) C DIMENSION V(MZCHF),A(MZCHF,MZCHF),B(MZCHF,MZCHF),C(MZCHF,MZCHF), + DR(MZCHF,MZCHF),DI(MZCHF,MZCHF) DIMENSION DNA(MZCHF,MZCHF),DNB(MZCHF,MZCHF) C DO J=1,NCHF2 DO I=1,NCHF1 DNA(I,J)=TZERO DNB(I,J)=TZERO ENDDO ENDDO C C GET DN=TRANSITION AMPLITUDE IN REACTANCE-MATRIX FORM C DNA DO J=1,NCHOP2 DO I=1,NCHOP1 DNA(I,J)=DNA(I,J)+ASS(I,J) DO IP=1,NCHF1 DNA(I,J)=DNA(I,J)+RK1(IP,I)*ACS(IP,J) ENDDO ENDDO ENDDO DO J=1,NCHOP2 DO JP=1,NCHF2 DO I=1,NCHOP1 DNA(I,J)=DNA(I,J)+RK2(JP,J)*ASC(I,JP) ENDDO ENDDO ENDDO DO I=1,NCHOP1 DO JP=1,NCHF2 A(JP,I)=TZERO DO IP=1,NCHF1 A(JP,I)=A(JP,I)+RK1(IP,I)*ACC(IP,JP) ENDDO ENDDO ENDDO DO J=1,NCHOP2 DO I=1,NCHOP1 DO JP=1,NCHF2 DNA(I,J)=DNA(I,J)+A(JP,I)*RK2(JP,J) ENDDO ENDDO ENDDO C DNB DO J=1,NCHOP2 DO I=1,NCHOP1 DNB(I,J)=DNB(I,J)+BSS(I,J) DO IP=1,NCHF1 DNB(I,J)=DNB(I,J)+RK1(IP,I)*BCS(IP,J) ENDDO ENDDO ENDDO DO J=1,NCHOP2 DO JP=1,NCHF2 DO I=1,NCHOP1 DNB(I,J)=DNB(I,J)+RK2(JP,J)*BSC(I,JP) ENDDO ENDDO ENDDO DO I=1,NCHOP1 DO JP=1,NCHF2 B(JP,I)=TZERO DO IP=1,NCHF1 B(JP,I)=B(JP,I)+RK1(IP,I)*BCC(IP,JP) ENDDO ENDDO ENDDO DO J=1,NCHOP2 DO I=1,NCHOP1 DO JP=1,NCHF2 DNB(I,J)=DNB(I,J)+B(JP,I)*RK2(JP,J) ENDDO ENDDO ENDDO c c ALL DO J=1,NCHOP2 DO I=1,NCHOP1 DN(I,J)=DN(I,J)+DNA(I,J)+DNB(I,J) ENDDO ENDDO C C CHANGE TO S-MATRIX FORM. REAL-PART=DR, IMAGINARY-PART=DI C DO J=1,NCHOP1 DO I=1,J A(I,J)=TZERO ENDDO DO K=1,NCHOP1 DO I=1,J A(I,J)=A(I,J)+RK1(I,K)*RK1(K,J) ENDDO ENDDO ENDDO C DO J=1,NCHOP1 DO I=1,J A(J,I)=A(I,J) ENDDO ENDDO DO I=1,NCHOP1 A(I,I)=A(I,I)+ONE ENDDO C CALL VERTS(A,MZCHF,NCHOP1,WORK,IERR) IF (IERR.NE.0) THEN WRITE(6,*)'STOP BECAUSE NO INVERSE FOUND IN SR.DTOT' STOP 'STOP BECAUSE NO INVERSE FOUND IN SR.DTOT' END IF C DO J=1,NCHOP1 DO I=J,NCHOP1 A(J,I)=A(I,J) ENDDO ENDDO C CALL MXMUL(A,DN,B,NCHOP1,NCHOP1,NCHOP2) C DO J=1,NCHOP2 DO I=1,J A(I,J)=TZERO ENDDO DO K=1,NCHOP2 DO I=1,J A(I,J)=A(I,J)+RK2(I,K)*RK2(K,J) ENDDO ENDDO ENDDO DO J=1,NCHOP2 DO I=1,J A(J,I)=A(I,J) ENDDO ENDDO DO I=1,NCHOP2 A(I,I)=A(I,I)+ONE ENDDO C CALL VERTS(A,MZCHF,NCHOP2,WORK,IERR) IF (IERR.NE.0) THEN WRITE(6,*)'STOP BECAUSE NO INVERSE FOUND IN SR.DTOT' STOP 'STOP BECAUSE NO INVERSE FOUND IN SR.DTOT' END IF C DO J=1,NCHOP2 DO I=J,NCHOP2 A(J,I)=A(I,J) ENDDO ENDDO C CALL MXMUL(B,A,C,NCHOP1,NCHOP2,NCHOP2) C DO J=1,NCHOP2 DO I=1,NCHOP1 A(I,J)=FOUR*C(I,J) ! - -> + ENDDO ENDDO C CALL MXMUL(RK1,A,B,NCHOP1,NCHOP1,NCHOP2) CALL MXMUL(B,RK2,DR,NCHOP1,NCHOP2,NCHOP2) C DO J=1,NCHOP2 DO I=1,NCHOP1 DR(I,J)=A(I,J)+DR(I,J) ! -DR -> +DR ENDDO ENDDO C CALL MXMUL(RK1,A,B,NCHOP1,NCHOP1,NCHOP2) CALL MXMUL(A,RK2,C,NCHOP1,NCHOP2,NCHOP2) C DO J=1,NCHOP2 DO I=1,NCHOP1 DI(I,J)=-B(I,J)+C(I,J) ! B -> -B ENDDO ENDDO C RETURN END C C*********************************************************************** SUBROUTINE DVECT(E1) C IMPLICIT REAL*8 (A-H,O-Z) C LOGICAL QDT1 C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) PARAMETER (LWORK=MZCHF*MZCHF) PARAMETER (MWORK=MZMNP) C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) C COMMON/CRMAT1/VALUE1(MZMNP),WMAT1(MZMNP,MZCHF),MNP1 COMMON/CHAN1/ECH1(MZCHF),CC1(MZCHF),LLCH1(MZCHF), X L2P1(MZCHF),KJ1(MZCHF), + NCHF1 COMMON/CDKKP/MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, 1 MNP2D2,NCHND2,LRGLD2,NPTYD2,NPAD, 2 DKL(MNPEXT,MNPEXT), CV X DKV(MNPEXT,MNPEXT), 3 AMAT(MZCHF,MZCHF),BMATL(MZCHF,MZCHF) CV X,BMATV(MZCHF,MZCHF) COMMON/DEP1/RK1(MZCHF,MZCHF),SE(MZPTS,MZCHF),CE(MZPTS,MZCHF), + NCHOP1,QDT1,F1(MZCHF,MZCHF),FP1(MZCHF,MZCHF),FKNU1(MZCHF) COMMON/CDV/DV(MNPEXT,MZCHF) COMMON/CBUT/IBUT,K10,K20 COMMON/SCALE/AZ,AZP,AZ2,AZ4,AZP2 COMMON/NRBWRK/WORK(LWORK),TEMP1(MWORK),TEMP2(MWORK) C DIMENSION V1(MNPEXT,MZCHF) C DO K1=1,MNP1 TEMP1(K1)=ONE/(VALUE1(K1)-E1) ENDDO DO J1=1,NCHOP1 DO K1=1,MNP1 V1(K1,J1)=TZERO ENDDO DO I=1,NCHF1 DO K1=1,MNP1 V1(K1,J1)=V1(K1,J1)+WMAT1(K1,I)*FP1(I,J1) ENDDO ENDDO DO K1=1,MNP1 V1(K1,J1)=V1(K1,J1)*TEMP1(K1) ENDDO C DO I=1,NCHF1 V1(MNP1+I,J1)=AZ*FP1(I,J1) ENDDO DO K2=1,K20 DV(K2,J1)=TZERO DO K1=1,K10 DV(K2,J1)=DV(K2,J1)+V1(K1,J1)*DKL(K1,K2) ENDDO ENDDO ENDDO C RETURN END C C*********************************************************************** SUBROUTINE IJWBK1(E1,NCHOP1) C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C COMMON/CHAN1/ECH1(MZCHF),CC1(MZCHF),LLCH1(MZCHF), X L2P1(MZCHF),KJ1(MZCHF), + NCHF1 C J=0 DO I=1,NCHOP1 EPS=E1-ECH1(I) LL=LLCH1(I) CALL JWBK1(EPS,LL,J) J=J+15 ENDDO RETURN END C C*********************************************************************** SUBROUTINE IJWBK2(E2,NCHOP2) C IMPLICIT REAL*8(A-H,O-Z) C INCLUDE 'PARAM' C COMMON/CHAN2/ECH2(MZCHF),CC2(MZCHF),LLCH2(MZCHF), X L2P2(MZCHF),KJ2(MZCHF), + NCHF2 J=0 DO I=1,NCHOP2 EPS=E2-ECH2(I) LL=LLCH2(I) CALL JWBK2(EPS,LL,J) J=J+15 ENDDO RETURN END C*************************************************************** C SUBROUTINE JWBK1(E,L,J) C C COMPUTES ARRAY D WHICH IS HELD IN COMMON/CJWBK/ AND C USED FOR CALCULATION OF JWBK FUNCTIONS. CNRB C NEUTRAL CASE ADDED C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MX15N=15*MZCHF) C COMMON/CJWBK/D1(MX15N),D2(MX15N) COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,TINY,DELTA,LACC COMMON/NRBZED/TZED C D1(J+1)=E C IF(TZED.EQ.0)THEN FK=SQRT(E) D1(J+2)=FK IF(L.GT.0)THEN D1(J+3)=1./FK C=DBLE(L*(L+1)) D1(J+4)=C SC=SQRT(C) D1(J+5)=SC D1(J+6)=(C+.125)/SC D1(J+7)=E*C D1(J+12)=6.*E*C D1(J+14)=-C*C D1(J+15)=-L*1.5707963 ELSE D1(J+4)=0.0 D1(J+15)=0.0 ENDIF RETURN ENDIF C IF(L.GT.0)GOTO 10 C C CASE OF L.EQ.0 D1(J+4)=0. IF(E.EQ.0)GOTO 30 C CASE OF L.EQ.0 AND E.GT.0 FK=SQRT(E) D1(J+2)=FK D1(J+3)=1./FK GOTO 30 C 10 IF(E.GT.0)GOTO 20 C C CASE OF L.GT.0 AND E.EQ.0 C=DBLE(L*(L+1)) D1(J+4)=C SC=SQRT(C) D1(J+5)=SC D1(J+6)=(C+.125)/SC D1(J+13)=6.*C D1(J+14)=-C*C GOTO 30 C C CASE OF L.GT.0 AND E.GT.0 20 FK=SQRT(E) D1(J+2)=FK D1(J+3)=1./FK C=DBLE(L*(L+1)) D1(J+4)=C SC=SQRT(C) D1(J+5)=SC D1(J+6)=(C+.125)/SC A=1.+E*C D1(J+7)=A A=3.*A D1(J+8)=A-1. D1(J+9)=A+1. D1(J+10)=FK*C D1(J+11)=-4.*E D1(J+12)=-9.+2.*A D1(J+13)=6.*C D1(J+14)=-C*C C C TERM IN ARG GAMMA ETC 30 D1(J+15)=ARGC(E,L,AC) C RETURN END C C*********************************************************************** SUBROUTINE JWBK2(E,L,J) C C COMPUTES ARRAY D WHICH IS HELD IN COMMON/CJWBK/ AND C USED FOR CALCULATION OF JWBK FUNCTIONS. CNRB C NEUTRAL CASE ADDED C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (MX15N=15*MZCHF) C COMMON/CJWBK/D1(MX15N),D2(MX15N) COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,TINY,DELTA,LACC COMMON/NRBZED/TZED C D2(J+1)=E C IF(TZED.EQ.0)THEN FK=SQRT(E) D2(J+2)=FK IF(L.GT.0)THEN D2(J+3)=1./FK C=DBLE(L*(L+1)) D2(J+4)=C SC=SQRT(C) D2(J+5)=SC D2(J+6)=(C+.125)/SC D2(J+7)=E*C D2(J+12)=6.*E*C D2(J+14)=-C*C D2(J+15)=-L*1.5707963 ELSE D2(J+4)=0.0 D2(J+15)=0.0 ENDIF RETURN ENDIF C IF(L.GT.0)GOTO 10 C C CASE OF L.EQ.0 D2(J+4)=0. IF(E.EQ.0)GOTO 30 C CASE OF L.EQ.0 AND E.GT.0 FK=SQRT(E) D2(J+2)=FK D2(J+3)=1./FK GOTO 30 C 10 IF(E.GT.0)GOTO 20 C C CASE OF L.GT.0 AND E.EQ.0 C=DBLE(L*(L+1)) D2(J+4)=C SC=SQRT(C) D2(J+5)=SC D2(J+6)=(C+.125)/SC D2(J+13)=6.*C D2(J+14)=-C*C GOTO 30 C C CASE OF L.GT.0 AND E.GT.0 20 FK=SQRT(E) D2(J+2)=FK D2(J+3)=1./FK C=DBLE(L*(L+1)) D2(J+4)=C SC=SQRT(C) D2(J+5)=SC D2(J+6)=(C+.125)/SC A=1.+E*C D2(J+7)=A A=3.*A D2(J+8)=A-1. D2(J+9)=A+1. D2(J+10)=FK*C D2(J+11)=-4.*E D2(J+12)=-9.+2.*A D2(J+13)=6.*C D2(J+14)=-C*C C C TERM IN ARG GAMMA ETC 30 D2(J+15)=ARGC(E,L,AC) C RETURN END C C*********************************************************************** SUBROUTINE MXMUL(A,B,C,N1,N2,N3) C C CALCULATES C=A*B C C NRB: VECTOR OR SCALAR VERSIONS C IMPLICIT REAL*8(A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (ZERO=0.0) C DIMENSION A(MZCHF,MZCHF),B(MZCHF,MZCHF),C(MZCHF,MZCHF) C CCRAY CRAY DO J=1,N3 CRAY DO I=1,N1 CRAY C(I,J)=ZERO CRAY ENDDO CRAY ENDDO CRAY DO K=1,N2 CRAY DO J=1,N3 CRAY DO I=1,N1 CRAY C(I,J)=C(I,J)+A(I,K)*B(K,J) CRAY ENDDO CRAY ENDDO CRAY ENDDO CCRAY C---- DO J=1,N3 DO I=1,N1 C(I,J)=ZERO ENDDO DO K=1,N2 DO I=1,N1 C(I,J)=C(I,J)+A(I,K)*B(K,J) ENDDO ENDDO ENDDO C---- RETURN END C C*********************************************************************** SUBROUTINE NUMFC1(E1,I,INOUT,C1,C2) C IMPLICIT REAL*8 (A-H,O-Z) C LOGICAL QDT1 C INCLUDE 'PARAM' C PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (TEN=10.0) PARAMETER (TWELVE=12.0) C COMMON/CHAN1/ECH1(MZCHF),CC1(MZCHF),LLCH1(MZCHF), X L2P1(MZCHF),KJ1(MZCHF), + NCHF1 COMMON/DEP1/RK1(MZCHF,MZCHF),SE(MZPTS,MZCHF),CE(MZPTS,MZCHF), + NCHOP1,QDT1,F1(MZCHF,MZCHF),FP1(MZCHF,MZCHF),FKNU1(MZCHF) COMMON/DPP/SE1(MZCHF),SO1(MZCHF),CE1(MZCHF),CO1(MZCHF) COMMON/CPOINT/RZERO,RTWO,H,KP2 COMMON/NRBZED/TZED C NRB: C NEUTRAL CASE ADDED C c V(X)=EQ+X*(Q*TZED-X*CQ) c Q=H*H/TWELVE EQ=Q*(E1-ECH1(I)) CQ=Q*CC1(I) Q=Q*TWO C IF(INOUT.EQ.0)THEN CE(KP2,I)=C1 CE(KP2-1,I)=C2 R=RTWO X=ONE/R U1=V(X) R=R-H X=ONE/R U2=V(X) DO K=KP2-2,1,-1 R=R-H X=ONE/R U3=V(X) CE(K,I)=((TWO-TEN*U2)*CE(K+1,I)-(ONE+U1)*CE(K+2,I))/(ONE+U3) U1=U2 U2=U3 ENDDO ELSE CE(1,I)=C1 CE(2,I)=C2 R=RZERO X=ONE/R U1=V(X) R=R+H X=ONE/R U2=V(X) DO K=3,KP2 R=R+H X=ONE/R U3=V(X) CE(K,I)=((TWO-TEN*U2)*CE(K-1,I)-(ONE+U1)*CE(K-2,I))/(ONE+U3) U1=U2 U2=U3 ENDDO ENDIF CE1(I)=CE(1,I) C CALL THETA1(RTWO,I,T,TP) C RETURN END C C*********************************************************************** SUBROUTINE NUMFC2(E2,I,INOUT,C1,C2) C IMPLICIT REAL*8 (A-H,O-Z) C LOGICAL QDT2 C INCLUDE 'PARAM' C PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (TEN=10.0) PARAMETER (TWELVE=12.0) C COMMON/CHAN2/ECH2(MZCHF),CC2(MZCHF),LLCH2(MZCHF), X L2P2(MZCHF),KJ2(MZCHF), + NCHF2 COMMON/DEP2/RK2(MZCHF,MZCHF),SO(MZPTS,MZCHF),CO(MZPTS,MZCHF), + NCHOP2,QDT2,F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FKNU2(MZCHF) COMMON/DPP/SE1(MZCHF),SO1(MZCHF),CE1(MZCHF),CO1(MZCHF) COMMON/CPOINT/RZERO,RTWO,H,KP2 COMMON/NRBZED/TZED C NRB: C NEUTRAL CASE ADDED C c V(X)=EQ+X*(Q*TZED-X*CQ) c Q=H*H/TWELVE EQ=Q*(E2-ECH2(I)) CQ=Q*CC2(I) Q=Q*TWO C IF(INOUT.EQ.0)THEN CO(KP2,I)=C1 CO(KP2-1,I)=C2 R=RTWO X=ONE/R U1=V(X) R=R-H X=ONE/R U2=V(X) DO K=KP2-2,1,-1 R=R-H X=ONE/R U3=V(X) CO(K,I)=((TWO-TEN*U2)*CO(K+1,I)-(ONE+U1)*CO(K+2,I))/(ONE+U3) U1=U2 U2=U3 ENDDO ELSE CO(1,I)=C1 CO(2,I)=C2 R=RZERO X=ONE/R U1=V(X) R=R+H X=ONE/R U2=V(X) DO K=3,KP2 R=R+H X=ONE/R U3=V(X) CO(K,I)=((TWO-TEN*U2)*CO(K-1,I)-(ONE+U1)*CO(K-2,I))/(ONE+U3) U1=U2 U2=U3 ENDDO ENDIF CO1(I)=CO(1,I) C CALL THETA2(RTWO,I,T,TP) C RETURN END C C*********************************************************************** SUBROUTINE NUMFO1(E1,I,INOUT,S1,S2,C1,C2) C IMPLICIT REAL*8 (A-H,O-Z) C LOGICAL QDT1 C INCLUDE 'PARAM' C PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (TEN=10.0) PARAMETER (TWELVE=12.0) C COMMON/CHAN1/ECH1(MZCHF),CC1(MZCHF),LLCH1(MZCHF), X L2P1(MZCHF),KJ1(MZCHF), + NCHF1 COMMON/DEP1/RK1(MZCHF,MZCHF),SE(MZPTS,MZCHF),CE(MZPTS,MZCHF), + NCHOP1,QDT1,F1(MZCHF,MZCHF),FP1(MZCHF,MZCHF),FKNU1(MZCHF) COMMON/DPP/SE1(MZCHF),SO1(MZCHF),CE1(MZCHF),CO1(MZCHF) COMMON/CPOINT/RZERO,RTWO,H,KP2 COMMON/NRBZED/TZED C NRB: C NEUTRAL CASE ADDED C V(X)=EQ+X*(Q*TZED-X*CQ) C Q=H*H/TWELVE EQ=Q*(E1-ECH1(I)) CQ=Q*CC1(I) Q=Q*TWO C SE(1,I)=S1 SE(2,I)=S2 R=RZERO X=ONE/R U1=V(X) R=R+H X=ONE/R U2=V(X) IF(INOUT.EQ.0)THEN CE(1,I)=C1 CE(2,I)=C2 DO K=3,KP2 R=R+H X=ONE/R U3=V(X) D3=ONE/(ONE+U3) D2=(TWO-TEN*U2)*D3 D1=(ONE+U1)*D3 SE(K,I)=D2*SE(K-1,I)-D1*SE(K-2,I) CE(K,I)=D2*CE(K-1,I)-D1*CE(K-2,I) U1=U2 U2=U3 ENDDO ELSE DO K=3,KP2 R=R+H X=ONE/R U3=V(X) SE(K,I)=((TWO-TEN*U2)*SE(K-1,I)-(ONE+U1)*SE(K-2,I))/(ONE+U3) U1=U2 U2=U3 ENDDO CE(KP2,I)=C1 CE(KP2-1,I)=C2 R=RTWO X=ONE/R U1=V(X) R=R-H X=ONE/R U2=V(X) DO K=KP2-2,1,-1 R=R-H X=ONE/R U3=V(X) CE(K,I)=((TWO-TEN*U2)*CE(K+1,I)-(ONE+U1)*CE(K+2,I))/(ONE+U3) U1=U2 U2=U3 ENDDO ENDIF SE1(I)=SE(1,I) CE1(I)=CE(1,I) C RETURN END C C*********************************************************************** SUBROUTINE NUMFO2(E2,I,INOUT,S1,S2,C1,C2) C IMPLICIT REAL*8 (A-H,O-Z) C LOGICAL QDT2 C INCLUDE 'PARAM' C PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (TEN=10.0) PARAMETER (TWELVE=12.0) C COMMON/CHAN2/ECH2(MZCHF),CC2(MZCHF),LLCH2(MZCHF), X L2P2(MZCHF),KJ2(MZCHF), + NCHF2 COMMON/DEP2/RK2(MZCHF,MZCHF),SO(MZPTS,MZCHF),CO(MZPTS,MZCHF), + NCHOP2,QDT2,F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FKNU2(MZCHF) COMMON/DPP/SE1(MZCHF),SO1(MZCHF),CE1(MZCHF),CO1(MZCHF) COMMON/CPOINT/RZERO,RTWO,H,KP2 COMMON/NRBZED/TZED C NRB: C NEUTRAL CASE ADDED C V(X)=EQ+X*(Q*TZED-X*CQ) C Q=H*H/TWELVE EQ=Q*(E2-ECH2(I)) CQ=Q*CC2(I) Q=Q*TWO C SO(1,I)=S1 SO(2,I)=S2 R=RZERO X=ONE/R U1=V(X) R=R+H X=ONE/R U2=V(X) IF(INOUT.EQ.0)THEN CO(1,I)=C1 CO(2,I)=C2 DO K=3,KP2 R=R+H X=ONE/R U3=V(X) D3=ONE/(ONE+U3) D2=(TWO-TEN*U2)*D3 D1=(ONE+U1)*D3 SO(K,I)=D2*SO(K-1,I)-D1*SO(K-2,I) CO(K,I)=D2*CO(K-1,I)-D1*CO(K-2,I) U1=U2 U2=U3 ENDDO ELSE DO K=3,KP2 R=R+H X=ONE/R U3=V(X) SO(K,I)=((TWO-TEN*U2)*SO(K-1,I)-(ONE+U1)*SO(K-2,I))/(ONE+U3) U1=U2 U2=U3 ENDDO CO(KP2,I)=C1 CO(KP2-1,I)=C2 R=RTWO X=ONE/R U1=V(X) R=R-H X=ONE/R U2=V(X) DO K=KP2-2,1,-1 R=R-H X=ONE/R U3=V(X) CO(K,I)=((TWO-TEN*U2)*CO(K+1,I)-(ONE+U1)*CO(K+2,I))/(ONE+U3) U1=U2 U2=U3 ENDDO ENDIF SO1(I)=SO(1,I) CO1(I)=CO(1,I) C RETURN END C C*********************************************************************** FUNCTION PTARG(TEMP) C C PARTITION FUNCTION FOR TARGET C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C COMMON/CTARG/ENAT(MZTAR),NZED,NELC,NAST,NCONAT(MZTAR), + IST(MZTAR),ILT(MZTAR),IPT(MZTAR) COMMON/SCALE/AZ,AZP,AZ2,AZ4,AZP2 C RYDTZ=TEMP/(AZ2*157894.) PTARG=0. DO N=1,NAST if(ist(n).gt.0)WGHT=IST(N)*(2*ILT(N)+1) if(ist(n).eq.0)WGHT=(ILT(N)+1) !BP PTARG=PTARG+WGHT*EXP(-ENAT(N)/RYDTZ) ENDDO RETURN END C C*********************************************************************** SUBROUTINE RD00(ISE,ILE,KE,ISO,ILO,KO,KK1,KK2,ND,KPP,KASE) C IMPLICIT REAL*8 (A-H,O-Z) C LOGICAL KPP C INCLUDE 'PARAM' C PARAMETER (MXLSD=(MZSLP*(MZSLP+1))/2) C DIMENSION ISE(MZSLP),ILE(MZSLP),ISO(MZSLP),ILO(MZSLP) DIMENSION ND(MXLSD),KPP(MXLSD),KK1(MXLSD),KK2(MXLSD) DIMENSION KS1(MXLSD),KL1(MXLSD),KP1(MXLSD) X ,KS2(MXLSD),KL2(MXLSD),KP2(MXLSD) C OPEN(3,FILE='D00',STATUS='OLD',FORM='UNFORMATTED') C READ(3)KOUNT C IF(KOUNT.GT.MXLSD)THEN WRITE(6,*)'***SR.RD00: INCREASE MZSLP' STOP '***SR.RD00: INCREASE MZSLP' ENDIF C DO K=1,KOUNT READ(3)KS1(K),KL1(K),KP1(K),KS2(K),KL2(K),KP2(K) ENDDO CLOSE(3) C IDIP=1 if(ise(1).eq.0)IDIP=2 !BP holds 2*J C KASE=0 DO K1=1,KE IS1=ISE(K1) IL1=ILE(K1) IP1=0 DO K2=1,KO IS2=ISO(K2) IL2=ILO(K2) IP2=1 IF(IS1.EQ.IS2)THEN IF(IABS(IL1-IL2).LE.IDIP)THEN KASE=KASE+1 KK1(KASE)=K1 KK2(KASE)=K2 ENDIF ENDIF ENDDO ENDDO C DO K=1,KASE K1=KK1(K) K2=KK2(K) IS1=ISE(K1) IL1=ILE(K1) IP1=0 IS2=ISO(K2) IL2=ILO(K2) IP2=1 DO L=1,KOUNT IF(KS1(L).EQ.IS1.AND.KL1(L).EQ.IL1.AND.KP1(L).EQ.IP1.AND. + KS2(L).EQ.IS2.AND.KL2(L).EQ.IL2.AND.KP2(L).EQ.IP2)THEN ND(K)=L KPP(K)=.TRUE. GOTO 1 ENDIF ENDDO DO L=1,KOUNT IF(KS2(L).EQ.IS1.AND.KL2(L).EQ.IL1.AND.KP2(L).EQ.IP1.AND. + KS1(L).EQ.IS2.AND.KL1(L).EQ.IL2.AND.KP1(L).EQ.IP2)THEN ND(K)=L KPP(K)=.FALSE. GOTO 1 ENDIF ENDDO WRITE(6,601)IS1,IL1,IP1,IS2,IL2,IP2 KK1(K)=0 KK2(K)=0 1 CONTINUE ENDDO C RETURN 601 FORMAT('*** SR.RD00: WARNING, KASE NOT FOUND IN FILE D00: ' + ,3I3,2X,3I3) C END C C*********************************************************************** SUBROUTINE REDEP1(E1,M1,M11) C IMPLICIT REAL*8 (A-H,O-Z) C LOGICAL QDT1 C C READS FINAL ENERGY DEPENDENT DATA C INCLUDE 'PARAM' C PARAMETER (ONE=1.0) C COMMON/CHAN1/ECH1(MZCHF),CC1(MZCHF),LLCH1(MZCHF), X L2P1(MZCHF),KJ1(MZCHF), + NCHF1 COMMON/DEP1/RK1(MZCHF,MZCHF),SE(MZPTS,MZCHF),CE(MZPTS,MZCHF), + NCHOP1,QDT1,F1(MZCHF,MZCHF),FP1(MZCHF,MZCHF),FKNU1(MZCHF) COMMON/MISC/WBODE(MZPTS),IPERT,MXE,EMESH(MZMSH),BSTO C DO I=1,NCHF1 EPS=E1-ECH1(I) IF(EPS.GE.0)THEN FKNU1(I)=SQRT(EPS) ELSE FKNU1(I)=ONE/SQRT(-EPS) ENDIF ENDDO C READ(1) QDT1 READ(1) NCHOP1 IF(NCHOP1.EQ.0)RETURN READ(1) ((RK1(I,J),J=1,NCHOP1),I=1,NCHF1) READ(1) ((F1(I,J),J=1,NCHOP1),I=1,NCHF1) READ(1) ((FP1(I,J),J=1,NCHOP1),I=1,NCHF1) READ(1) ((FPP1IJ,J=1,NCHOP1),I=1,NCHF1) DO I=1,NCHOP1 READ(1)INOUT,S1,S2,C1,C2 IF(M1.GE.M11)THEN CALL NUMFO1(E1,I,INOUT,S1,S2,C1,C2) ENDIF ENDDO DO I=NCHOP1+1,NCHF1 READ(1)INOUT,C1,C2 IF(M1.GE.M11)THEN CALL NUMFC1(E1,I,INOUT,C1,C2) ENDIF ENDDO C IF(IPERT.GT.0)THEN READ(1,ERR=92)((DUM,J=1,NCHOP1),I=1,NCHF1) READ(1) READ(1) RETURN 92 BACKSPACE 1 ENDIF C RETURN C END C C*********************************************************************** SUBROUTINE REDEP2(E2,M1,M2,M21) C IMPLICIT REAL*8 (A-H,O-Z) C LOGICAL QDT2 C C READS FINAL ENERGY DEPENDENT DATA C INCLUDE 'PARAM' C PARAMETER (ONE=1.0) C COMMON/CHAN2/ECH2(MZCHF),CC2(MZCHF),LLCH2(MZCHF), X L2P2(MZCHF),KJ2(MZCHF), + NCHF2 COMMON/DEP2/RK2(MZCHF,MZCHF),SO(MZPTS,MZCHF),CO(MZPTS,MZCHF), + NCHOP2,QDT2,F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FKNU2(MZCHF) COMMON/MISC/WBODE(MZPTS),IPERT,MXE,EMESH(MZMSH),BSTO C DO I=1,NCHF2 EPS=E2-ECH2(I) IF(EPS.GE.0)THEN FKNU2(I)=SQRT(EPS) ELSE FKNU2(I)=ONE/SQRT(-EPS) ENDIF ENDDO C READ(2) QDT2 READ(2) NCHOP2 IF(NCHOP2.EQ.0)RETURN READ(2) ((RK2(I,J),J=1,NCHOP2),I=1,NCHF2) READ(2) ((F2(I,J),J=1,NCHOP2),I=1,NCHF2) READ(2) ((FP2(I,J),J=1,NCHOP2),I=1,NCHF2) READ(2) ((FPP2IJ,J=1,NCHOP2),I=1,NCHF2) DO I=1,NCHOP2 READ(2)INOUT,S1,S2,C1,C2 IF(M1.NE.M2.AND.M2.GE.M21)THEN CALL NUMFO2(E2,I,INOUT,S1,S2,C1,C2) ENDIF ENDDO DO I=NCHOP2+1,NCHF2 READ(2)INOUT,C1,C2 IF(M1.NE.M2.AND.M2.GE.M21)THEN CALL NUMFC2(E2,I,INOUT,C1,C2) ENDIF ENDDO C IF(IPERT.GT.0)THEN READ(2,ERR=92)((DUM,J=1,NCHOP2),I=1,NCHF2) READ(2) READ(2) RETURN 92 BACKSPACE 2 ENDIF C RETURN C END C C*********************************************************************** SUBROUTINE REIND1(KFIL,ISE,ILE) C IMPLICIT REAL*8 (A-H,O-Z) C CHARACTER*1 NUM(0:9) C C READS ENERGY INDEPENDENT DATA FOR EVEN PARITY C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) C PARAMETER (HALF=0.5) C COMMON/CHAN1/ECH1(MZCHF),CC1(MZCHF),LLCH1(MZCHF), X L2P1(MZCHF),KJ1(MZCHF), + NCHF1 COMMON/CRMAT1/VALUE1(MZMNP),WMAT1(MZMNP,MZCHF),MNP1 COMMON/CDKKP/MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, 1 MNP2D2,NCHND2,LRGLD2,NPTYD2,NPAD, 2 DKL(MNPEXT,MNPEXT), CV X DKV(MNPEXT,MNPEXT), 3 AMAT(MZCHF,MZCHF),BMATL(MZCHF,MZCHF) CV X,BMATV(MZCHF,MZCHF) COMMON/CBUT/IBUT,K10,K20 COMMON/MISC/WBODE(MZPTS),IPERT,MXE,EMESH(MZMSH),BSTO C DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ C OPEN(1,FILE='F'//NUM(KFIL/10)//NUM(KFIL-10*(KFIL/10)), +STATUS='OLD',FORM='UNFORMATTED') READ(1)IS,IL,IP C C CHECK IS,IL,IP IF(IP.NE.0)THEN PRINT*,' KFIL=',KFIL,', REIND1 HAS IP=',IP STOP ENDIF IF(IS.NE.ISE.OR.IL.NE.ILE)THEN PRINT*,' KFIL=',KFIL,' REIND1 HAS' PRINT*,' ISE, ILE=',ISE,ILE PRINT*,' IS, IL =',IS,IL STOP ENDIF C READ(1)MNP1,NCHF1 IF(MNP1.GT.MZMNP) THEN WRITE(6,613)MNP1 PRINT*,' REIND1 HAS MNP1=',MNP1,' GREATER THAN ', + 'MZMNP=',MZMNP STOP ENDIF IF(NCHF1.GT.MZCHF) THEN WRITE(6,614)NCHF1 PRINT*,' REIND1 HAS NCHF1=',NCHF1,' GREATER THAN ', + 'MZCHF=',MZCHF STOP ENDIF READ(1)(ECH1(I),I=1,NCHF1) IF(IS.NE.0)READ(1)(CC1(I),I=1,NCHF1) if(is.eq.0)read(1)(cc1(i),i=1,nchf1),(l2p1(i),i=1,nchf1) x,(idum,i=1,2*nchf1),(kj1(i),i=1,nchf1) READ(1)(VALUE1(I),I=1,MNP1) READ(1)((WMAT1(J,I),J=1,MNP1),I=1,NCHF1) IF(IPERT.GT.0)READ(1) C DO I=1,NCHF1 LLCH1(I)=NINT(SQRT(CC1(I)+HALF)-HALF) ENDDO C C CHECK NCHF1, MNP1 IF(NCHF1.NE.NCHND1)THEN WRITE(6,600)NCHF1,NCHND1 STOP ENDIF IF(MNP1.NE.MNP2D1)THEN WRITE(6,601)MNP1,MNP2D1 STOP ENDIF C K10=MNP1 IF(IBUT.GT.0)K10=K10+NCHF1 C RETURN 510 FORMAT(//' *** MISMATCH'/4X,'FILE ',A3,' HAS IS,IL,IP,MNP2 = ', 1 4I4//) 613 FORMAT(//' READS MNP2 = ',I4,' WHICH IS LARGER THAN ', 1 'MAXIMUM VALUE OF MZMNP ALLOWED BY DIMENSIONS'//) 614 FORMAT(//' READS NCHF = ',I4,' WHICH IS LARGER THAN MAXIMUM ', 1 'VALUE OF MZCHF (LESS 1 FOR FINAL STATE) ALLOWED BY', 2 ' DIMENSIONS'//) 600 FORMAT(' REIND1, NCHF1=',I2,', MCHND1=',I2) 601 FORMAT(' REIND1, MNP1=',I2,', MNP2D1=',I2) END C C*********************************************************************** SUBROUTINE REIND2(KFIL,ISO,ILO) C IMPLICIT REAL*8 (A-H,O-Z) C CHARACTER*1 NUM(0:9) C C READS ENERGY INDEPENDENT DATA FOR ODD PARITY C INCLUDE 'PARAM' C PARAMETER (MNPEXT=MZMNP+MZCHF) C PARAMETER (HALF=0.5) C COMMON/CHAN2/ECH2(MZCHF),CC2(MZCHF),LLCH2(MZCHF), X L2P2(MZCHF),KJ2(MZCHF), + NCHF2 COMMON/CRMAT2/VALUE2(MZMNP),WMAT2(MZMNP,MZCHF),MNP2 COMMON/CDKKP/MNP2D1,NCHND1,LRGLD1,NPTYD1,NSPND, 1 MNP2D2,NCHND2,LRGLD2,NPTYD2,NPAD, 2 DKL(MNPEXT,MNPEXT), CV X DKV(MNPEXT,MNPEXT), 3 AMAT(MZCHF,MZCHF),BMATL(MZCHF,MZCHF) CV X,BMATV(MZCHF,MZCHF) COMMON/CBUT/IBUT,K10,K20 COMMON/MISC/WBODE(MZPTS),IPERT,MXE,EMESH(MZMSH),BSTO C DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ C OPEN(2,FILE='F'//NUM(KFIL/10)//NUM(KFIL-10*(KFIL/10)), +STATUS='OLD',FORM='UNFORMATTED') REWIND(2) READ(2)IS,IL,IP C C CHECK IS,IL,IP IF(IP.NE.1)THEN PRINT*,' KFIL=',KFIL,', REIND1 HAS IP=',IP STOP ENDIF IF(IS.NE.ISO.OR.IL.NE.ILO)THEN PRINT*,' KFIL=',KFIL,' REIND1 HAS' PRINT*,' ISO, ILO=',ISO,ILO PRINT*,' IS, IL =',IS,IL STOP ENDIF C READ(2)MNP2,NCHF2 IF(MNP2.GT.MZMNP) THEN WRITE(6,613)MNP2 PRINT*,' REIND2 HAS MNP2=',MNP2,' GREATER THAN ', + 'MZMNP=',MZMNP STOP ENDIF IF(NCHF2.GT.MZCHF) THEN WRITE(6,614)NCHF2 PRINT*,' REIND2 HAS NCHF=',NCHF2,' GREATER THAN ', + 'MZCHF=',MZCHF STOP ENDIF READ(2)(ECH2(I),I=1,NCHF2) IF(IS.NE.0)READ(2)(CC2(I),I=1,NCHF2) if(is.eq.0)read(2)(cc2(i),i=1,nchf2),(l2p2(i),i=1,nchf2) x,(idum,i=1,2*nchf2),(kj2(i),i=1,nchf2) READ(2)(VALUE2(I),I=1,MNP2) READ(2)((WMAT2(J,I),J=1,MNP2),I=1,NCHF2) IF(IPERT.GT.0)READ(2) C DO I=1,NCHF2 LLCH2(I)=NINT(SQRT(CC2(I)+HALF)-HALF) ENDDO C C CHECK NCHF1, MNP1 IF(NCHF2.NE.NCHND2)THEN WRITE(6,600)NCHF2,NCHND2 STOP ENDIF IF(MNP2.NE.MNP2D2)THEN WRITE(6,601)MNP2,MNP2D2 STOP ENDIF c K20=MNP2 IF(IBUT.GT.0)K20=K20+NCHF2 C RETURN 510 FORMAT(//' *** MISMATCH'/4X,'FILE ',A3,' HAS IS,IL,IP,MNP2 = ', 1 4I4//) 613 FORMAT(//' READS MNP2 = ',I4,' WHICH IS LARGER THAN ', 1 'MAXIMUM VALUE OF MZMNP ALLOWED BY DIMENSIONS'//) 614 FORMAT(//' READS NCHF = ',I4,' WHICH IS LARGER THAN MAXIMUM ', 1 'VALUE OF MZCHF (LESS 1 FOR FINAL STATE) ALLOWED BY', 2 ' DIMENSIONS'//) 600 FORMAT(' REIND12, NCHF2=',I2,', MCHND2=',I2) 601 FORMAT(' REIND2, MNP2=',I2,', MNP2D2=',I2) END C C*********************************************************************** SUBROUTINE RF00 C IMPLICIT REAL*8 (A-H,O-Z) C C READ DIRECTORY FILE FOR F DATASETS C INCLUDE 'PARAM' C PARAMETER (ONE=1.0) C COMMON/LIST/ISEE(MZSLP),ILEE(MZSLP),KFILEE(MZSLP),KEE, + ISOO(MZSLP),ILOO(MZSLP),KFILOO(MZSLP),KOO COMMON/MISC/WBODE(MZPTS),IPERT,MXE,EMESH(MZMSH),BSTO COMMON/CPOINT/RZERO,RTWO,H,KP2 COMMON/SCALE/AZ,AZP,AZ2,AZ4,AZP2 COMMON/CTARG/ENAT(MZTAR),NZED,NELC,NAST,NCONAT(MZTAR), + IST(MZTAR),ILT(MZTAR),IPT(MZTAR) C iprint=1 !temp ISOO(1)=1 C C FREE STATE DATA OPEN(0,FILE='F00',STATUS='OLD',FORM='UNFORMATTED') C READ(0)NZED,NELC ! Target states READ(0)NAST,(ENAT(I),I=1,NAST),(IST(I),ILT(I),I=1,NAST) READ(0)KP2 ! Range for functions READ(0)RZERO,H ! Boundary and interval READ(0)(WBODE(K),K=1,KP2) READ(0)IPERT READ(0)MXE C WRITE(6,506)NZED,NELC IF(IPERT.EQ.0)WRITE(6,500) IF(IPRINT.GT.0)WRITE(6,503)RZERO,H,KP2 C C CHECK DIMENSION FOR KP2 IF(KP2.GT.MZPTS)THEN WRITE(6,600)KP2,MZPTS STOP ENDIF C C CHECK DIMENSION FOR MXE IF(MXE.GT.MZMSH)THEN WRITE(6,610)MXE,MZMSH STOP 'INCREASE DIMENSION MZMSH' ENDIF C READ(0)(EMESH(I),I=1,MXE) READ(0)BSTO C KSLP=0 C C AND LIST OF SLPI CASES KEE=0 KOO=0 6 READ(0)IS,IL,IP IF(IL.GE.0) THEN KSLP=KSLP+1 IF(IP.EQ.0)THEN KEE=KEE+1 ISEE(KEE)=IS ILEE(KEE)=IL KFILEE(KEE)=KSLP ELSE KOO=KOO+1 ISOO(KOO)=IS ILOO(KOO)=IL KFILOO(KOO)=KSLP ENDIF GOTO 6 ENDIF CLOSE(0) C RTWO=RZERO+H*(KP2-1) C C SCALE PARAMETERS AZ=MAX(NZED-NELC,1) AZP=ONE/AZ AZ2=AZ*AZ AZ4=AZ2*AZ2 AZP2=AZP*AZP C CALL BODE C RETURN C C 500 FORMAT(17X,'(STGF CALCULATION WITHOUT MUTIPOLE POTENTIALS)') 502 FORMAT(10X,I3,4X,I2,3X,I2,3X,I2) 503 FORMAT(/7X,'RZERO =',F9.4/7X,'H =',F8.4/7X,'KP2 =',I5/) 506 FORMAT(//14X,'NUCLEAR CHARGE =',I3, * ', NUMBER OF TARGET ELECTRONS =',I3/14X,52('-')//) 600 FORMAT(///10X,'KP2 =',I6, * ' IS GREATER THAN MZPTS = ',I5///) 610 FORMAT(///10X,'MXE = ',I7,' IS GREATER THAN', + ' MZMSH = ',I7///) C END C C*********************************************************************** SUBROUTINE RTHET1(I,R,TA,TP) C C CALCULATES THETA FOR CHANNEL I AND REAL R C THETA = TA*EXP(TP) C TP = FNUI*LOG(R) - R/FNUI C IMPLICIT REAL*8(A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) C CNRB C NEUTRAL CASE ADDED C COMMON/CTHET1/BB1(MZCHF,MZTET),NSUM1(MZCHF) COMMON/NRBZED/TZED C MI=NSUM1(I) FNUI=BB1(I,1) G=TWO*R/FNUI Y=ONE/G F=ONE S=BB1(I,2) DO M=3,MI F=F*Y S=S+BB1(I,M)*F ENDDO C DLR=LOG(R)*TZED TP=(FNUI*DLR-G/TWO) TA=S C RETURN END C C*********************************************************************** SUBROUTINE RTHET2(I,R,TA,TP) C C CALCULATES THETA FOR CHANNEL I AND REAL R C THETA = TA*EXP(TP) C TP = FNUI*LOG(R) - R/FNUI C IMPLICIT REAL*8(A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) C CNRB C NEUTRAL CASE ADDED C COMMON/CTHET2/BB2(MZCHF,MZTET),NSUM2(MZCHF) COMMON/NRBZED/TZED C MI=NSUM2(I) FNUI=BB2(I,1) G=TWO*R/FNUI Y=ONE/G F=ONE S=BB2(I,2) DO M=3,MI F=F*Y S=S+BB2(I,M)*F ENDDO C DLR=LOG(R)*TZED TP=(FNUI*DLR-G/TWO) TA=S C RETURN END C C*********************************************************************** REAL*8 FUNCTION SFCE(R,E1,C1,E2,C2,F1,FP1,F2,FP2) C IMPLICIT REAL*8 (A-H,O-Z) C PARAMETER (TWO=2.0) C XL=(F1*FP2-FP1*F2)*R XK=(C1-C2)/R**2 XI=F1*F2 FPP2=(C2/R**2-TWO/R-E2)*F2 S1=XL+XI S2=XK*(XL-XI)-TWO*(F1*FPP2-FP1*FP2) SFCE=-(S1*(E1-E2)+S2) C RETURN END C C*********************************************************************** SUBROUTINE TARG C IMPLICIT REAL*8 (A-H,O-Z) C LOGICAL EX C INCLUDE 'PARAM' C COMMON/CTARG/ENAT(MZTAR),NZED,NELC,NAST,NCONAT(MZTAR), + IST(MZTAR),ILT(MZTAR),IPT(MZTAR) C INQUIRE(FILE='T',EXIST=EX) IF(.NOT.EX)RETURN C OPEN(99,FILE='T',STATUS='OLD') READ(99,*)NT IF(NT.NE.NAST)THEN WRITE(6,600)NT,NAST STOP ENDIF DO N=1,NT READ(99,*)IST(N),ILT(N),IPT(N) ENDDO CLOSE(99) 600 FORMAT(' SUBROUTINE TARG READS NT=',I3/ + ' WHICH IS NOT EQUAL TO NAST=',I3) RETURN END C C*********************************************************************** SUBROUTINE THAV(TEMP,MM,FT,MDM) C C GIVES FT(MD) = THERMAL AVERAGE OF TOT FOR ENERGY DIFFERENCES OF C DELTA=DE*MD (DE=EMESH(N+1)-EMESH(N)). c IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) C COMMON/CD/DN(MZCHF,MZCHF),DA(MZCHF,MZCHF),STA(MZMSH),STB(MZMSH), + TOT(MZMSH) COMMON/MISC/WBODE(MZPTS),IPERT,MXE,EMESH(MZMSH),BSTO COMMON/SCALE/AZ,AZP,AZ2,AZ4,AZP2 DIMENSION P(MZMSH),FT(MZMSH) C RYDT=TEMP/(157894.*AZ2) C C LAST EXP=EXP0, LAST MD=MDM EXP0=0.01 F=-LOG(EXP0) DE=EMESH(2)-EMESH(1) MDM=MM-F*RYDT/DE C MM2=2*MM ALP=ONE/RYDT G=RYDT/(EMESH(2)-EMESH(1)) EMESH(MM+1)=TWO*EMESH(MM)-EMESH(MM-1) E2=EXP(-EMESH(2)*ALP) P(1)=E2-ONE E3=E2 E2=1. DO M=2,MM E1=E2 E2=E3 E3=EXP(-EMESH(M+1)*ALP) P(M)=E1-TWO*E2+E3 ENDDO DO MD=1,MDM FJ=TZERO DO M1=1,MM-MD M2=M1+MD N=((M1-1)*(MM2-M1))/2+M2 FJ=FJ+P(M1)*TOT(N) ENDDO M2=1+MD FT(MD)=G*FJ+TOT(M2) ENDDO C RETURN END C C*********************************************************************** SUBROUTINE THETA1(R,I,T,TP) C IMPLICIT REAL*8(A-H,O-Z) C C THETA IS CALLED FOR r=R, FNU=KKNU(I), LL=LLCH(I) C RETURNS T=THETA, TP=D(THETA)/Dr AND BB(I,n) in /CTHET/. C GIVES BB(I,1)=FNU. c FOR ANY r, THETA=r**FNU * EXP(-r/FNU) * S c WHERE S={SUN n=0, NSUM(I)}{BB(i,n+2)*Y**n} WITH Y=FNU/(2*r). C NORMALISATION IS BB(I,2)=R**(-.5*FNU) * EXP(.5*R/FNU). C LOGICAL QDT1 C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) C CNRB C NEUTRAL CASE ADDED C COMMON/CHAN1/ECH1(MZCHF),CC1(MZCHF),LLCH1(MZCHF), X L2P1(MZCHF),KJ1(MZCHF), + NCHF1 COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,TINY,DELTA,LACC COMMON/CTHET1/BB1(MZCHF,MZTET),NSUM1(MZCHF) COMMON/DEP1/RK1(MZCHF,MZCHF),SE(MZPTS,MZCHF),CE(MZPTS,MZCHF), + NCHOP1,QDT1,F1(MZCHF,MZCHF),FP1(MZCHF,MZCHF),FKNU1(MZCHF) COMMON/CPOINT/RZERO,RTWO,H,KP2 COMMON/NRBZED/TZED C LOGICAL FLG C FNU=FKNU1(I) LL=LLCH1(I) M=FNU+LL+12 C IF(M.LE.MZTET)GO TO 10 IF(TZED.EQ.0)THEN M=MZTET !SHOULD CONVERGE AT M=LL GO TO 10 ENDIF PRINT*,' THETA1: M=',M,', MZTET=',MZTET PRINT*,' FNU=',FNU,', LL=',LL PRINT*,' Stopping' STOP C 10 FF=ONE/FNU X=TWO*R*FF Y=ONE/X FL=DBLE(LL) R1=ONE/R C BB1(I,1)=FNU BB1(I,2)=ONE C YN=ONE A1=FL-FNU*TZED A2=FL+FNU*TZED+ONE BN=-TWO*FNU-ONE C BET=ONE C S=ONE U=TZERO C N=0 FLG=.FALSE. DO MN=3,M N=N+1 A1=A1+ONE A2=A2-ONE BN=BN+TWO CN=A1*A2 BET=CN*BET YN=YN*Y U=U+BET*YN AN=ONE/DBLE(N) BET=BET*AN S=S+BET*YN BB1(I,MN)=BET IF(ABS(BET*YN).LT.AC*ABS(S))THEN IF(FLG)THEN GOTO 30 ELSE FLG=.TRUE. ENDIF ENDIF ENDDO C C NOT CONVERGED WRITE(6,610)N STOP C 30 NSUM1(I)=N+2 C C RE-NORMALISE T=S*EXP(-R/FNU) IF(TZED.GT.0)T=T*R**FNU G=CE(KP2,I)/T DO J=2,NSUM1(I) BB1(I,J)=BB1(I,J)*G ENDDO c RETURN C 610 FORMAT(/3X,'SUBROUTINE THETA1'/ 1 ' SUMMATIONS NOT CONVERGED WITH ',I3,' TERMS') END C C*********************************************************************** SUBROUTINE THETA2(R,I,T,TP) C IMPLICIT REAL*8(A-H,O-Z) C C THETA IS CALLED FOR r=R, FNU=KKNU(I), LL=LLCH(I) C RETURNS T=THETA, TP=D(THETA)/Dr AND BB(I,n) in /CTHET/. C GIVES BB(I,1)=FNU. c FOR ANY r, THETA=r**FNU * EXP(-r/FNU) * S c WHERE S={SUN n=0, NSUM(I)}{BB(i,n+2)*Y**n} WITH Y=FNU/(2*r). C NORMALISATION IS BB(I,2)=R**(-.5*FNU) * EXP(.5*R/FNU). C LOGICAL QDT2 C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) C CNRB C NEUTRAL CASE ADDED C COMMON/CHAN2/ECH2(MZCHF),CC2(MZCHF),LLCH2(MZCHF), X L2P2(MZCHF),KJ2(MZCHF), + NCHF2 COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,TINY,DELTA,LACC COMMON/CTHET2/BB2(MZCHF,MZTET),NSUM2(MZCHF) COMMON/DEP2/RK2(MZCHF,MZCHF),SO(MZPTS,MZCHF),CO(MZPTS,MZCHF), + NCHOP2,QDT2,F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FKNU2(MZCHF) COMMON/CPOINT/RZERO,RTWO,H,KP2 COMMON/NRBZED/TZED C LOGICAL FLG C FNU=FKNU2(I) LL=LLCH2(I) M=FNU+LL+12 C IF(M.LE.MZTET)GO TO 10 IF(TZED.EQ.0)THEN M=MZTET !SHOULD CONVERGE AT M=LL GO TO 10 ENDIF PRINT*,' THETA2: M=',M,', MZTET=',MZTET PRINT*,' FNU=',FNU,', LL=',LL PRINT*,' Stopping' STOP C 10 FF=ONE/FNU X=TWO*R*FF Y=ONE/X FL=DBLE(LL) R1=ONE/R C BB2(I,1)=FNU BB2(I,2)=ONE C YN=ONE A1=FL-FNU*TZED A2=FL+FNU*TZED+ONE BN=-TWO*FNU-ONE C BET=ONE C S=ONE U=TZERO C N=0 FLG=.FALSE. DO MN=3,M N=N+1 A1=A1+ONE A2=A2-ONE BN=BN+TWO CN=A1*A2 BET=CN*BET YN=YN*Y U=U+BET*YN AN=ONE/DBLE(N) BET=BET*AN S=S+BET*YN BB2(I,MN)=BET IF(ABS(BET*YN).LT.AC*ABS(S))THEN IF(FLG)THEN GOTO 30 ELSE FLG=.TRUE. ENDIF ENDIF ENDDO C C NOT CONVERGED WRITE(6,610)N STOP C 30 NSUM2(I)=N+2 C C RE-NORMALISE T=S*EXP(-R/FNU) IF(TZED.GT.0)T=T*R**FNU G=CO(KP2,I)/T DO J=2,NSUM2(I) BB2(I,J)=BB2(I,J)*G ENDDO c RETURN C 610 FORMAT(/3X,'SUBROUTINE THETA2'/ 1 ' SUMMATIONS NOT CONVERGED WITH ',I3,' TERMS') END C C*********************************************************************** SUBROUTINE TLAG(I,J,NLAG,LIJ,T1) C C CALCULATES T INTEGRALS FOR CHANNELS I,J USING C LAGUERRE QUADRATURE WITH NLAG POINTS C THE INTEGRALS ARE - C T1 FOR (THETAI,THETAJ) C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) C COMMON/CTHET1/BB1(MZCHF,MZTET),NSUM1(MZCHF) COMMON/CTHET2/BB2(MZCHF,MZTET),NSUM2(MZCHF) COMMON/CBLK/XLAG(30),WLAG(30),XLEG(15),WLEG(15) COMMON/CPOINT/RZERO,RTWO,H,KP2 C C CAP K FNUI=BB1(I,1) FNUJ=BB2(J,1) FK=ONE/FNUI+ONE/FNUJ C INITIALISE FOR QUADRATURE NS=NLAG/2 M=NS*(NS-1) N1=M+1 N2=M+NLAG T1=TZERO C START QUADRATURE DO N=N1,N2 U=XLAG(N) R=RTWO+U/FK C CALCULATE THETA FUNCTIONS CALL RTHET1(I,R,TI,TPI) CALL RTHET2(J,R,TJ,TPJ) C ADD TO SUM C ++ VAX MOD C A1=(R**(-LIJ))*WLAG(N)*EXP(U+TPI+TPJ) A1=(R**(-LIJ))*WLAG(N) U2=U/TWO AI=EXP(U2+TPI) AJ=EXP(U2+TPJ) TI=TI*AI TJ=TJ*AJ C ++ END MOD T1=T1+TI*A1*TJ ENDDO F1=ONE/FK T1=T1*F1 C RETURN END C C*********************************************************************** SUBROUTINE TOPUP(ise,ILE,ILO,MM) C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C COMMON/CD/DN(MZCHF,MZCHF),DA(MZCHF,MZCHF),STA(MZMSH),STB(MZMSH), + TOT(MZMSH) COMMON/MISC/WBODE(MZPTS),IPERT,MXE,EMESH(MZMSH),BSTO C MM2=2*MM FLAM2=MAX(ILE,ILO)**2 if(ise.eq.0)flam2=flam2/4 !BP DO MD=1,MM-1 DO M1=1,MM-MD M2=M1+MD E1=EMESH(M1) E2=EMESH(M2) N1=((M1-1)*(MM2-M1))/2+M2 IF(ILE.GT.ILO)THEN TOT(N1)=TOT(N1)+((1.+FLAM2*E2)*STA(N1) + -(1.+FLAM2*E1)*STB(N1))/(FLAM2*(E1-E2)) ELSE TOT(N1)=TOT(N1)+((1.+FLAM2*E2)*STB(N1) + -(1.+FLAM2*E1)*STA(N1))/(FLAM2*(E1-E2)) ENDIF ENDDO ENDDO RETURN END C C*********************************************************************** SUBROUTINE VERTS(V,LV,N,W,IERR) C C ________________________________________________________ C | | C | INVERT A SYMMETRIC MATRIX WITHOUT PIVOTING | C | | C | INPUT: | C | | C | V --ARRAY CONTAINING MATRIX | C | (ONLY THE LOWER HALF NEED BE DEFINED) | C | | C | LV --LEADING (ROW) DIMENSION OF ARRAY V | C | | C | N --MATRIX DIMENSION | C | | C | W --WORK ARRAY WITH LENGTH AT LEAST N | C | | C | OUTPUT: | C | | C | V --INVERSE (IN LOWER HALF ONLY) | C |________________________________________________________| C REAL*8 V(*),W(*),S,T INTEGER G,H,I,J,K,L,M,N,LV,IERR REAL*8 TZERO,ONE C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) C IERR=0 C C ---------------- CNRB |*** PACK V ***| C ---------------- H=LV-N I=0 M=0 L=N G=(N*(N+1))/2 2 IF(L.EQ.G)GO TO 4 K=L+1 M=M+1 L=L+N-M I=I+H+M DO J=K,L V(J)=V(I+J) ENDDO GO TO 2 C 4 H = N K = 1 10 IF ( H .EQ. 1 ) GOTO 40 C -------------------------- C |*** SAVE PIVOT ENTRY ***| C -------------------------- S = V(K) K = K + H G = K H = H - 1 M = H IF ( S .EQ. TZERO ) GOTO 50 J = 0 20 J = J - M M = M - 1 L = G + M T = V(G+J)/S C --------------------------- C |*** ELIMINATE BY ROWS ***| C --------------------------- DO 30 I = G,L 30 V(I) = V(I) - T*V(I+J) G = L + 1 IF ( M .GT. 0 ) GOTO 20 GOTO 10 40 IF ( V(K) .NE. TZERO ) GOTO 60 IERR=2 RETURN 50 IERR=1 RETURN C ------------------------------------------ C |*** SOLVE FOR ROWS OF INVERSE MATRIX ***| C ------------------------------------------ 60 G = N + N DO 150 M = 1,N L = ((G-M)*(M-1))/2 H = L K = M DO 70 I = M,N 70 W(I) = TZERO W(M) = ONE 80 IF ( K .EQ. N ) GOTO 100 T = W(K)/V(K+L) J = L L = L + N - K K = K + 1 IF ( T .EQ. TZERO ) GOTO 80 DO 90 I = K,N 90 W(I) = W(I) - T*V(I+J) GOTO 80 C ----------------------------------- C |*** BACK SUBSTITUTION BY ROWS ***| C ----------------------------------- 100 W(N) = W(N)/V(K+L) 110 IF ( K .EQ. M ) GOTO 130 J = K K = K - 1 L = L + K - N T = W(K) DO 120 I = J,N 120 T = T - W(I)*V(I+L) W(K) = T/V(K+L) GOTO 110 130 DO 140 I = M,N 140 V(I+H) = W(I) 150 CONTINUE C ------------------ CNRB |*** UNPACK V ***| C ------------------ H=LV-N I=(N-1)*H+(N*(N-1))/2 M=N-1 L=(N*(N+1))/2 K=L 200 IF(I.EQ.0)RETURN DO J=L,K,-1 V(I+J)=V(J) ENDDO I=I-H-M L=K-1 M=M-1 K=K-N+M GO TO 200 C END C C*********************************************************************** REAL*8 FUNCTION W1(A,B,C,D,F) C C CALCULATES RACAH COEFFICIENT W(A,B,C,D,1,F) IMPLICIT INTEGER(A-Q) IMPLICIT REAL*8 (R-Y) C Q=(-1)**(B+D-F) C IF(B-A)10,20,30 C 10 IF(D-C)11,12,13 20 IF(D-C)21,22,23 30 IF(D-C)31,32,33 C 11 X=(B+D+F+2)*(B+D+F+3)*(B+D-F+1)*(B+D-F+2) Y=(B+1)*(2*B+1)*(2*B+3)*(D+1)*(2*D+1)*(2*D+3) GOTO 40 12 X=(B+D+F+2)*(-B+D+F)*(B-D+F+1)*(B+D-F+1) Y=(B+1)*(2*B+1)*(2*B+3)*D*(D+1)*(2*D+1) GOTO 40 13 X=(-B+D+F-1)*(-B+D+F)*(B-D+F+1)*(B-D+F+2) Y=(B+1)*(2*B+1)*(2*B+3)*D*(2*D-1)*(2*D+1) GOTO 40 C 21 X=(B+D+F+2)*(-B+D+F+1)*(B-D+F)*(B+D-F+1) Y=B*(B+1)*(2*B+1)*(D+1)*(2*D+1)*(2*D+3) GOTO 40 22 X=B*(B+1)+D*(D+1)-F*(F+1) Y=B*(B+1)*(2*B+1)*D*(D+1)*(2*D+1) Q=-Q W1=Q*.5*X/SQRT(Y) RETURN 23 X=(B+D+F+1)*(-B+D+F)*(B-D+F+1)*(B+D-F) Y=B*(B+1)*(2*B+1)*D*(2*D-1)*(2*D+1) Q=-Q GOTO 40 C 31 X=(-B+D+F+1)*(-B+D+F+2)*(B-D+F-1)*(B-D+F) Y=B*(2*B-1)*(2*B+1)*(D+1)*(2*D+1)*(2*D+3) GOTO 40 32 X=(B+D+F+1)*(-B+D+F+1)*(B-D+F)*(B+D-F) Y=B*(2*B-1)*(2*B+1)*D*(D+1)*(2*D+1) Q=-Q GOTO 40 33 X=(B+D+F)*(B+D+F+1)*(B+D-F-1)*(B+D-F) Y=B*(2*B-1)*(2*B+1)*D*(2*D-1)*(2*D+1) GOTO 40 C C 40 W1=.5*Q*SQRT(X/Y) C RETURN END C C*********************************************************************** C REAL*8 FUNCTION W2(A,B,C,D,F) C C CALCULATES W(A,B,C,D,1,F) WHERE W IS RACAH COEFFICIENT. C NRB: C*** AND *** A,B,C,D,E,F ARE INPUT TWICE THEIR ACTUAL VALUE. C IMPLICIT INTEGER(A-N) IMPLICIT REAL*8 (W) C C ISIGN=(-1)**((F-A-C)/2) C IF(B-A)10,20,30 C 10 IF(D-C)11,12,13 20 IF(D-C)21,22,23 30 IF(D-C)31,32,33 C 11 WSQ2=.25*DBLE((B+D+F+4)*(B+D+F+6)*(B+D-F+2)*(B+D-F+4))/ + DBLE((B+2)*(B+1)*(B+3)*(D+2)*(D+1)*(D+3))/4. GO TO 40 12 WSQ2=.25*DBLE((B+D+F+4)*(-B+D+F)*(B-D+F+2)*(B+D-F+2))/ + DBLE((B+2)*(B+1)*(B+3)*D*(D+2)*(D+1))/2. ISIGN=-ISIGN GO TO 40 13 WSQ2=.25*DBLE((-B+D+F-2)*(-B+D+F)*(B-D+F+2)*(B-D+F+4))/ + DBLE((B+2)*(B+1)*(B+3)*D*(D-1)*(D+1))/4. GO TO 40 C 21 WSQ2=.25*DBLE((B+D+F+4)*(-B+D+F+2)*(B-D+F)*(B+D-F+2))/ + DBLE(B*(B+2)*(B+1)*(D+2)*(D+1)*(D+3))/2. ISIGN=-ISIGN GO TO 40 22 G=-(B*(B+2)+D*(D+2)-F*(F+2)) IF(G.LT.0)ISIGN=-ISIGN WSQ2=.25*DBLE(G*G)/ + DBLE(B*(B+2)*(B+1)*D*(D+2)*(D+1)) GO TO 40 23 WSQ2=.25*DBLE((B+D+F+2)*(-B+D+F)*(B-D+F+2)*(B+D-F))/ + DBLE(B*(B+2)*(B+1)*D*(D-1)*(D+1))/2. GO TO 40 C 31 WSQ2=.25*DBLE((-B+D+F+2)*(-B+D+F+4)*(B-D+F-2)*(B-D+F))/ + DBLE(B*(B-1)*(B+1)*(D+2)*(D+1)*(D+3))/4 GO TO 40 32 WSQ2=.25*DBLE((B+D+F+2)*(-B+D+F+2)*(B-D+F)*(B+D-F))/ + DBLE(B*(B-1)*(B+1)*D*(D+2)*(D+1))/2. GO TO 40 33 WSQ2=.25*DBLE((B+D+F)*(B+D+F+2)*(B+D-F-2)*(B+D-F))/ + DBLE(B*(B-1)*(B+1)*D*(D-1)*(D+1))/4. GOTO 40 C C 40 W2=ISIGN*SQRT(WSQ2) C END C*********************************************************************** SUBROUTINE ZPHI1(I,ZR,ZAI,ZPI) C C COMPUTES AMPLITUDE ZAI AND PHASE ZPI OF COULOMB FUNCTION ZPHI C FOR COMPLEX RADIAL CO-ORDINATE ZR. C USES DATA IN ARRAY D WHICH IS HELD IN COMMON/CJWBK/ C AND SHOULD HAVE BEEN COMPUTED IN SUBROUTINE INJWBK. C THE STRUCTURE OF ZPHI IS SIMILAR TO THAT OF SUBROUTINE JWBK. CNRB C NEUTRAL CASE ADDED C C IMPLICIT REAL*8(A-H,O-Y) IMPLICIT COMPLEX*16(Z) C INCLUDE 'PARAM' C PARAMETER (MX15N=15*MZCHF) C COMMON/CJWBK/D1(MX15N),D2(MX15N) COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,TINY,DELTA,LACC COMMON/NRBZED/TZED C C J=(I-1)*15 E=D1(J+1) C=D1(J+4) ZX=1./ZR C IF(TZED.EQ.0)THEN FK=D1(J+2) IF(C.EQ.0)THEN ZAI=1.0/SQRT(FK) ZPI=FK*ZR ELSE ZXSQ=ZX*ZX ZW=E-C*ZXSQ ZA1=0.0625*(ZX/ZW)**3 ZCC=ZA1*(D1(J+14)*ZXSQ+D1(J+12))*ZX ZWH=SQRT(ZW) Z=ZR*ZWH ZSQ=Z*Z ZP=(ZSQ-0.125-.2083333*C/ZSQ)/Z+D1(J+15) ZA1=1./(ZR*FK) ZS=D1(J+5)*ZA1 ZG=Z*ZA1 ZP=ZP+D1(J+6)*(0.,-1.)*LOG(ZG+(0.,1.)*ZS) ZPI=ZP ZAI=(1.-ZCC)/SQRT(ZWH) ENDIF RETURN ENDIF C IF(C.EQ.0)GOTO 30 IF(E.EQ.0)GOTO 70 C C CASE OF C.GT.0 AND E.GT.0 ZW=E+ZX*(2.-C*ZX) ZWH=SQRT(ZW) Z=ZR*ZWH FK=D1(J+2) ZRK=ZR*FK ZRMC=ZR-C ZALP=Z+ZRK CK=D1(J+10) C COMPUTE PHASE ZP=Z+D1(J+15) C LOG TERM ZB=FK*ZALP IF(ABS(ZB).GT.ACJWBK)GOTO 10 ZB=-ZB ZP=ZP+ZALP*((((.2*ZB+.25)*ZB+.33333333)*ZB+.5)*ZB+1.) GOTO 20 10 ZP=ZP+D1(J+3)*LOG(1.+ZB) C ARCTAN TERM 20 ZA1=1./(ZR*D1(J+7)) ZS=D1(J+5)*(Z-FK*ZRMC)*ZA1 ZG=(CK*Z+ZRMC)*ZA1 ZP=ZP+D1(J+6)*(0.,-1.)*LOG(ZG+(0.,1.)*ZS) C CAP PHI TERM ZP=ZP+((5.*ZRMC/(Z*Z))-(Z*D1(J+9)+ZRK*D1(J+8)+CK)/ C (ZALP*D1(J+7)))/(24.*Z) C COMPLETE CALCULATION OF ZPHI ZA1=.0625*(ZX/ZW)**3 ZCC=ZA1*(((D1(J+14)*ZX+D1(J+13))*ZX+D1(J+12))*ZX+D1(J+11)) ZPI=ZP ZAI=(1.-ZCC)/SQRT(ZWH) RETURN C 30 IF(E.EQ.0)GOTO 60 C C CASE OF C.EQ.0 AND E.GT.0 ZW=2.*ZX+E ZWH=SQRT(ZW) Z=ZR*ZWH FK=D1(J+2) ZRK=ZR*FK ZALP=Z+ZRK C COMPUTE PHASE ZP=Z+D1(J+15) ZB=FK*ZALP IF(ABS(ZB).GT.ACJWBK)GOTO 40 ZB=-ZB ZP=ZP+ZALP*((((.2*ZB+.25)*ZB+.33333333)*ZB+.5)*ZB+1.) GOTO 50 40 ZP=ZP+D1(J+3)*LOG(1.+ZB) 50 ZP=ZP+1/(4.*ZALP)+(5.*ZR/(Z*Z)-2.*(Z+ZALP)/ZALP)/(24.*Z) C COMPLETE CALCULATION OF ZPHI ZA1=.0625*(ZX/ZW)**3 ZCC=ZA1*(-4.*E-3.*ZX) ZPI=ZP ZAI=(1.-ZCC)/SQRT(ZWH) RETURN C C CASE OF C.EQ.0 AND E.EQ.0 60 ZW=2.*ZX ZWH=SQRT(ZW) Z=ZR*ZWH ZP=2.*Z*(1.+.046875*ZX)+D1(J+15) ZWMQ=1./SQRT(ZWH) ZET=(1.+.0234375*ZX)*ZWMQ ZAI=ZET ZPI=ZP RETURN C C CASE OF E.EQ.0 AND C.GT.0 70 ZW=ZX*(2.-C*ZX) ZWH=SQRT(ZW) Z=ZR*ZWH ZRMC=ZR-C C COMPUTE PHASE ZP=2.*Z+D1(J+15) ZA1=1./ZR ZS=D1(J+5)*Z*ZA1 ZG=ZRMC*ZA1 ZP=ZP+D1(J+6)*(0.,-1.)*LOG(ZG+(0.,1.)*ZS) ZP=ZP-(3.*ZR+C)/(24.*(ZRMC+ZR)*Z) C COMPLETE CALCULATION OF ZPHI ZA1=.0625*(ZX/ZW)**3 ZCC=((D1(J+14)*ZX+D1(J+13))*ZX-3.)*ZX*ZA1 ZAI=(1.-ZCC)/SQRT(ZWH) ZPI=ZP C RETURN END C C*********************************************************************** SUBROUTINE ZPHI2(I,ZR,ZAI,ZPI) C C COMPUTES AMPLITUDE ZAI AND PHASE ZPI OF COULOMB FUNCTION ZPHI C FOR COMPLEX RADIAL CO-ORDINATE ZR. C USES DATA IN ARRAY D WHICH IS HELD IN COMMON/CJWBK/ C AND SHOULD HAVE BEEN COMPUTED IN SUBROUTINE INJWBK. C THE STRUCTURE OF ZPHI IS SIMILAR TO THAT OF SUBROUTINE JWBK. CNRB C NEUTRAL CASE ADDED C C IMPLICIT REAL*8(A-H,O-Y) IMPLICIT COMPLEX*16(Z) C INCLUDE 'PARAM' C PARAMETER (MX15N=15*MZCHF) C COMMON/CJWBK/D1(MX15N),D2(MX15N) COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,TINY,DELTA,LACC COMMON/NRBZED/TZED C C J=(I-1)*15 E=D2(J+1) C=D2(J+4) ZX=1./ZR C IF(TZED.EQ.0)THEN FK=D2(J+2) IF(C.EQ.0)THEN ZAI=1.0/SQRT(FK) ZPI=FK*ZR ELSE ZXSQ=ZX*ZX ZW=E-C*ZXSQ ZA1=0.0625*(ZX/ZW)**3 ZCC=ZA1*(D2(J+14)*ZXSQ+D2(J+12))*ZX ZWH=SQRT(ZW) Z=ZR*ZWH ZSQ=Z*Z ZP=(ZSQ-0.125-.2083333*C/ZSQ)/Z+D2(J+15) ZA1=1./(ZR*FK) ZS=D2(J+5)*ZA1 ZG=Z*ZA1 ZP=ZP+D2(J+6)*(0.,-1.)*LOG(ZG+(0.,1.)*ZS) ZPI=ZP ZAI=(1.-ZCC)/SQRT(ZWH) ENDIF RETURN ENDIF C IF(C.EQ.0)GOTO 30 IF(E.EQ.0)GOTO 70 C C CASE OF C.GT.0 AND E.GT.0 ZW=E+ZX*(2.-C*ZX) ZWH=SQRT(ZW) Z=ZR*ZWH FK=D2(J+2) ZRK=ZR*FK ZRMC=ZR-C ZALP=Z+ZRK CK=D2(J+10) C COMPUTE PHASE ZP=Z+D2(J+15) C LOG TERM ZB=FK*ZALP IF(ABS(ZB).GT.ACJWBK)GOTO 10 ZB=-ZB ZP=ZP+ZALP*((((.2*ZB+.25)*ZB+.33333333)*ZB+.5)*ZB+1.) GOTO 20 10 ZP=ZP+D2(J+3)*LOG(1.+ZB) C ARCTAN TERM 20 ZA1=1./(ZR*D2(J+7)) ZS=D2(J+5)*(Z-FK*ZRMC)*ZA1 ZG=(CK*Z+ZRMC)*ZA1 ZP=ZP+D2(J+6)*(0.,-1.)*LOG(ZG+(0.,1.)*ZS) C CAP PHI TERM ZP=ZP+((5.*ZRMC/(Z*Z))-(Z*D2(J+9)+ZRK*D2(J+8)+CK)/ C (ZALP*D2(J+7)))/(24.*Z) C COMPLETE CALCULATION OF ZPHI ZA1=.0625*(ZX/ZW)**3 ZCC=ZA1*(((D2(J+14)*ZX+D2(J+13))*ZX+D2(J+12))*ZX+D2(J+11)) ZPI=ZP ZAI=(1.-ZCC)/SQRT(ZWH) RETURN C 30 IF(E.EQ.0)GOTO 60 C C CASE OF C.EQ.0 AND E.GT.0 ZW=2.*ZX+E ZWH=SQRT(ZW) Z=ZR*ZWH FK=D2(J+2) ZRK=ZR*FK ZALP=Z+ZRK C COMPUTE PHASE ZP=Z+D2(J+15) ZB=FK*ZALP IF(ABS(ZB).GT.ACJWBK)GOTO 40 ZB=-ZB ZP=ZP+ZALP*((((.2*ZB+.25)*ZB+.33333333)*ZB+.5)*ZB+1.) GOTO 50 40 ZP=ZP+D2(J+3)*LOG(1.+ZB) 50 ZP=ZP+1/(4.*ZALP)+(5.*ZR/(Z*Z)-2.*(Z+ZALP)/ZALP)/(24.*Z) C COMPLETE CALCULATION OF ZPHI ZA1=.0625*(ZX/ZW)**3 ZCC=ZA1*(-4.*E-3.*ZX) ZPI=ZP ZAI=(1.-ZCC)/SQRT(ZWH) RETURN C C CASE OF C.EQ.0 AND E.EQ.0 60 ZW=2.*ZX ZWH=SQRT(ZW) Z=ZR*ZWH ZP=2.*Z*(1.+.046875*ZX)+D2(J+15) ZWMQ=1./SQRT(ZWH) ZET=(1.+.0234375*ZX)*ZWMQ ZAI=ZET ZPI=ZP RETURN C C CASE OF E.EQ.0 AND C.GT.0 70 ZW=ZX*(2.-C*ZX) ZWH=SQRT(ZW) Z=ZR*ZWH ZRMC=ZR-C C COMPUTE PHASE ZP=2.*Z+D2(J+15) ZA1=1./ZR ZS=D2(J+5)*Z*ZA1 ZG=ZRMC*ZA1 ZP=ZP+D2(J+6)*(0.,-1.)*LOG(ZG+(0.,1.)*ZS) ZP=ZP-(3.*ZR+C)/(24.*(ZRMC+ZR)*Z) C COMPLETE CALCULATION OF ZPHI ZA1=.0625*(ZX/ZW)**3 ZCC=((D2(J+14)*ZX+D2(J+13))*ZX-3.)*ZX*ZA1 ZAI=(1.-ZCC)/SQRT(ZWH) ZPI=ZP C RETURN END C C*********************************************************************** COMPLEX*16 FUNCTION ZPLAG(I,J,NLAG,LIJ,EPS1,EPS2) C C CALCULATES P INTEGRALS USING LAGUERRE QUADRATURE C IMPLICIT REAL*8(A-H,O-Y) IMPLICIT COMPLEX*16(Z) C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (EIGHT=8.0) PARAMETER (ZI=(0.0,1.0)) C COMMON/CBLK/XLAG(30),WLAG(30),XLEG(15),WLEG(15) COMMON/CPOINT/RZERO,RTWO,H,KP2 COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,TINY,DELTA,LACC C FK=SQRT(EPS1)+SQRT(EPS2) X=RTWO B=SQRT(EIGHT*X) G=FK*B/EIGHT IF(FK.GT.TZERO)THEN GM=ONE/G G2=ONE+G G2=G2*G2 ENDIF ZB=ZI/B C NS=NLAG/2 M=NS*(NS-1) N1=M+1 N2=M+NLAG C ZPLAG=TZERO DO N=N1,N2 U=XLAG(N) A1=FK*U IF(A1.LE.ACZP)THEN ZET=ONE+ZB*U/TWO ELSE ZET=(SQRT(G2+ZB*G*U)-ONE)*GM ENDIF ZMU=-EIGHT*ZB*(G+ONE/ZET) ZET=ZET*ZET ZR=ZET*X CALL ZPHI1(I,ZR,ZAI,ZPI) CALL ZPHI2(J,ZR,ZAJ,ZPJ) ZPLAG=ZPLAG+ZAI*ZAJ*EXP(ZI*(ZPI+ZPJ)+U)*(ZR**(-LIJ))*WLAG(N)/ZMU ENDDO C RETURN END C C*********************************************************************** COMPLEX*16 FUNCTION ZQLAG(I,J,NLAG,LIJ,EPS1,EPS2) C C CALCULATES Q INTEGRALS USING LAGUERRE QUADRATURE C IMPLICIT REAL*8(A-H,O-Y) IMPLICIT COMPLEX*16(Z) C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (ZI=(0.0,1.0)) C COMMON/CBLK/XLAG(30),WLAG(30),XLEG(15),WLEG(15) COMMON/CPOINT/RZERO,RTWO,H,KP2 C FK=SQRT(EPS1)-SQRT(EPS2) X=RTWO ZMUM=ZI/FK C NS=NLAG/2 M=NS*(NS-1) N1=M+1 N2=M+NLAG C ZQLAG=TZERO DO N=N1,N2 U=XLAG(N) ZR=X+U*ZMUM CALL ZPHI1(I,ZR,ZAI,ZPI) CALL ZPHI2(J,ZR,ZAJ,ZPJ) ZQLAG=ZQLAG+ZAI*ZAJ*EXP(ZI*(ZPI-ZPJ)+U)*(ZR**(-LIJ))*WLAG(N) ENDDO C ZQLAG=ZQLAG*ZMUM C RETURN END C C*********************************************************************** COMPLEX*16 FUNCTION ZQLEG(N,M,NLEG,LIJ,EPS1,EPS2) C C CALCULATES Q INTEGRALS USING LEGENDRE QUADRATURE C IMPLICIT REAL*8(A-H,O-Y) IMPLICIT COMPLEX*16(Z) C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (ZI=(0.0,1.0)) C COMMON/CBLK/XLAG(30),WLAG(30),XLEG(15),WLEG(15) COMMON/CPOINT/RZERO,RTWO,H,KP2 LOGICAL SWAP C FKN=SQRT(EPS1) FKM=SQRT(EPS2) IF(FKN.GE.FKM)THEN I=N J=M FK=FKN-FKM SWAP=.FALSE. ELSE I=M J=N FK=FKM-FKN SWAP=.TRUE. ENDIF C X=RTWO IF(FK.EQ.TZERO)THEN ZA=ONE ELSE ZA=ZI/(ONE+FK*X) ENDIF C NS=NLEG/2 J1=(NS*(NS-1))/2 J2=J1+NS J1=J1+1 C ZQLEG=TZERO DO JJ=J1,J2 V=XLEG(JJ) C ZR1=X*(ONE+ZA*(ONE-V)/(ONE+V)) IF(SWAP)THEN CALL ZPHI2(I,ZR1,ZAI,ZPI) CALL ZPHI1(J,ZR1,ZAJ,ZPJ) ELSE CALL ZPHI1(I,ZR1,ZAI,ZPI) CALL ZPHI2(J,ZR1,ZAJ,ZPJ) ENDIF ZF=ZAI*ZAJ*EXP(ZI*(ZPI-ZPJ))* + (ZR1**(-LIJ))/(ONE+V)**2 C ZR2=X*(ONE+ZA*(ONE+V)/(ONE-V)) IF(SWAP)THEN CALL ZPHI2(I,ZR2,ZAI,ZPI) CALL ZPHI1(J,ZR2,ZAJ,ZPJ) ELSE CALL ZPHI1(I,ZR2,ZAI,ZPI) CALL ZPHI2(J,ZR2,ZAJ,ZPJ) ENDIF ZF=ZF+ZAI*ZAJ*EXP(ZI*(ZPI-ZPJ))* + (ZR2**(-LIJ))/(ONE-V)**2 ZQLEG=ZQLEG+ZF*WLEG(JJ) ENDDO ZQLEG=TWO*X*ZA*ZQLEG C IF(SWAP)ZQLEG=CONJG(ZQLEG) C RETURN END C C*********************************************************************** SUBROUTINE ZSLG12(I,J,NLAG,LIJ,ZS3) C C CALCULATES S INTEGRALS USING LAGUERRE QUADRATURE C IMPLICIT REAL*8(A-H,O-Y) IMPLICIT COMPLEX*16(Z) C LOGICAL QDT1,QDT2 C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (ZI=(0.0,1.0)) C COMMON/CBLK/XLAG(30),WLAG(30),XLEG(15),WLEG(15) COMMON/CPOINT/RZERO,RTWO,H,KP2 COMMON/DEP1/RK1(MZCHF,MZCHF),SE(MZPTS,MZCHF),CE(MZPTS,MZCHF), + NCHOP1,QDT1,F1(MZCHF,MZCHF),FP1(MZCHF,MZCHF),FKNU1(MZCHF) COMMON/DEP2/RK2(MZCHF,MZCHF),SO(MZPTS,MZCHF),CO(MZPTS,MZCHF), + NCHOP2,QDT2,F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FKNU2(MZCHF) C ZFK=DCMPLX(FKNU1(I),ONE/FKNU2(J)) X=RTWO B=SQRT(TWO*X) ZB=ZI/B ZG=ZFK*B/TWO C ZA1=ONE/ZG ZA2=ONE+ZG ZA2=ZA2*ZA2 ZA3=ZB*ZG C NS=NLAG/2 M=NS*(NS-1) N1=M+1 N2=M+NLAG C ZS3=TZERO DO N=N1,N2 U=XLAG(N) ZET=ZA1*(SQRT(ZA2+ZA3*U)-ONE) ZMUM=-ZET/(ZB*(ONE+ZET*ZG)/TWO) ZET=ZET*ZET ZR=ZET*X CALL ZPHI1(I,ZR,ZAI,ZPI) CALL ZTHET2(J,ZR,ZTAJ,ZTPJ) ZB1=(ZR**(-LIJ))*WLAG(N)*ZMUM*EXP(ZI*ZPI+ZTPJ+U)*ZAI ZS3=ZS3+ZB1*ZTAJ ENDDO C RETURN END C C*********************************************************************** SUBROUTINE ZSLG21(I,J,NLAG,LIJ,ZS3) C C CALCULATES S INTEGRALS USING LAGUERRE QUADRATURE C IMPLICIT REAL*8(A-H,O-Y) IMPLICIT COMPLEX*16(Z) C LOGICAL QDT1,QDT2 C INCLUDE 'PARAM' C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) PARAMETER (ZI=(0.0,1.0)) C COMMON/CBLK/XLAG(30),WLAG(30),XLEG(15),WLEG(15) COMMON/CPOINT/RZERO,RTWO,H,KP2 COMMON/DEP1/RK1(MZCHF,MZCHF),SE(MZPTS,MZCHF),CE(MZPTS,MZCHF), + NCHOP1,QDT1,F1(MZCHF,MZCHF),FP1(MZCHF,MZCHF),FKNU1(MZCHF) COMMON/DEP2/RK2(MZCHF,MZCHF),SO(MZPTS,MZCHF),CO(MZPTS,MZCHF), + NCHOP2,QDT2,F2(MZCHF,MZCHF),FP2(MZCHF,MZCHF),FKNU2(MZCHF) C ZFK=DCMPLX(FKNU2(J),ONE/FKNU1(I)) X=RTWO B=SQRT(TWO*X) ZB=ZI/B ZG=ZFK*B/TWO C ZA1=ONE/ZG ZA2=ONE+ZG ZA2=ZA2*ZA2 ZA3=ZB*ZG C NS=NLAG/2 M=NS*(NS-1) N1=M+1 N2=M+NLAG C ZS3=TZERO DO N=N1,N2 U=XLAG(N) ZET=ZA1*(SQRT(ZA2+ZA3*U)-ONE) ZMUM=-ZET/(ZB*(ONE+ZET*ZG)/TWO) ZET=ZET*ZET ZR=ZET*X CALL ZPHI2(J,ZR,ZAI,ZPI) CALL ZTHET1(I,ZR,ZTAJ,ZTPJ) ZB1=(ZR**(-LIJ))*WLAG(N)*ZMUM*EXP(ZI*ZPI+ZTPJ+U)*ZAI ZS3=ZS3+ZB1*ZTAJ ENDDO C RETURN END C C*********************************************************************** SUBROUTINE ZTHET1(I,ZR,ZTA,ZTP) C C CALCULATES THETA FOR CHANNEL I AND COMPLEX ZR C THETA = ZTA*CEXP(ZTP) C ZTP = FNUI*LOG(ZR) - ZR/FNUI C IMPLICIT REAL*8(A-H,O-Y) IMPLICIT COMPLEX*16(Z) C INCLUDE 'PARAM' C PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) C CNRB C NEUTRAL CASE ADDED C COMMON/CTHET1/BB1(MZCHF,MZTET),NSUM1(MZCHF) COMMON/NRBZED/TZED C MI=NSUM1(I) FNUI=BB1(I,1) Z=TWO*ZR/FNUI ZY=ONE/Z ZAS=ONE ZS=BB1(I,2) DO M=3,MI ZAS=ZAS*ZY ZS=ZS+BB1(I,M)*ZAS ENDDO C ZDLR=LOG(ZR)*TZED ZTP=(FNUI*ZDLR-Z/TWO) ZTA=ZS C RETURN END C C*********************************************************************** SUBROUTINE ZTHET2(I,ZR,ZTA,ZTP) C C CALCULATES THETA FOR CHANNEL I AND COMPLEX ZR C THETA = ZTA*CEXP(ZTP) C ZTP = FNUI*LOG(ZR) - ZR/FNUI C IMPLICIT REAL*8(A-H,O-Y) IMPLICIT COMPLEX*16(Z) C INCLUDE 'PARAM' C PARAMETER (ONE=1.0) PARAMETER (TWO=2.0) C CNRB C NEUTRAL CASE ADDED C COMMON/CTHET2/BB2(MZCHF,MZTET),NSUM2(MZCHF) COMMON/NRBZED/TZED C MI=NSUM2(I) FNUI=BB2(I,1) Z=TWO*ZR/FNUI ZY=ONE/Z ZAS=ONE ZS=BB2(I,2) DO M=3,MI ZAS=ZAS*ZY ZS=ZS+BB2(I,M)*ZAS ENDDO C ZDLR=LOG(ZR)*TZED ZTP=(FNUI*ZDLR-Z/TWO) ZTA=ZS C RETURN END