C N. R. BADNELL UoS v1.5 21/10/99 C C PROGRAM IMPACT C ************** C C CREES, SEATON, STOREY, TAYLOR & WILSON (1981) C C THIS VERSION (BADNELL 1999) REQUIRES NO PREPROCESSING AND MAKES C USE OF BLAS/LAPACK LIBRARIES WHICH MUST BE LINKED BY THE USER. C C--------------------------------------------------------------------- C C SOLVES COUPLED INTEGRO-DIFFERENTIAL EQUATIONS OF ELECTRON-ATOM C COLLISION THEORY C C COMPUTES - ELECTRON-ATOM COLLISION STRENGTHS C - ATOMIC STRUCTURES IN FROZEN CORES APPROXIMATION C C REFERENCES - GENERAL THEORY, EISSNER AND SEATON,J. PHYS. B.,5,2187, C 1972 C - FROZEN CORES APPROXIMATION, SEATON AND WILSON, C J. PHYS. B., 5,L1,1972 C - NUMERICAL METHODS, SEATON, J. PHYS. B, 7, 1817, 1974 C - FOR ORIGINAL VERSIONS OF PROGRAMS IMPACT AND IMPPRO C CREES, SEATON AND WILSON, COMPUT. PHYS. COMMUN. C 15, (1978) 23. REFERRED TO AS CPC ON COMMENT CARDS. C - PROGRAM ASYPCK, CREES COMPUT. PHYS. COMMUN. 19 C (1980) 103. ASYPCK2, 23 (1981) 181. C C--------------------------------------------------------------------- C C THE BASIC DIMENSIONS ARE INSERTED BY AN INCLUDE STATEMENT THAT C READS A FILE CALLED PARAM - NRB. C (SEE CPC TABLE 3, PLUS MFG) C C DIMENSIONS - MAXIMUM VALUE OF C KM = MZ01 C NTERM = MZ02 C NRAD = MZ03 C JMAX = MZ04 C ITC = MZ05 C N = MZ06 C NCHF = MZ07 C NTOT = MZ08 C MFG+1 = MZ22 C C TOTAL SIZE OF LARGE ARRAY IS: MZ08*(MZ08+MZ07) C C--------------------------------------------------------------------- C C ADDITIONAL FEATURES NOT DESCRIBED IN CPC C ======================================== C C 1. ENERGY DEPENDENT MFG C C IN THE CPC VERSION ONLY ONE VALUE OF MFG IS ALLOWED FOR A RUN C INVOLVING SEVERAL VALUES OF THE ELECTRON ENERGY. IT HAS BEEN C FOUND USEFUL TO PERMIT THE SPECIFICATION OF A DIFFERENT VALUE C OF MFG FOR EACH ENERGY, BY USING AN ARRAY MFGK(K). THE PROGRAM C WILL READ IN A VALUE FOR MFG IN THE STANDARD (I.E. PUBLISHED) C MANNER. IMPACT WILL USE THIS VALUE FOR ALL ENERGIES UNLESS A C VALUE OF MFG FOR A GIVEN ENERGY HAS BEEN READ IN TO OVERWRITE C THIS, UNDER STATEMENT 2100 (SEE ALSO COMMENTS FOLLOWING CARD). C C 2. CHANGE IN CRITERIA FOR CALLING ITERA METHOD FOR SOLUTIONS C IN THE ASYMPTOTIC REGION C C FOR THE CASES OF POSITIVE IONS AND ENERGIES JUST ABOVE A C THRESHOLD, THERE IS A CHANGE IN THE CRITERIA FOR CHOICE OF C USING THE EXPAN OR ITERA METHOD FOR OBTAINING ASYMPTOTIC C SOLUTIONS. DETAILS ARE GIVEN IN COMMENT CARDS IN SUBROUTINE C ASYSM AND IN SECTION 9 OF THE WRITE-UP OF THE PROGRAM ASYPCK C (CREES, COMPUT. PHYS. COMMUN., 19 (1980) 103) C C 3. BLAS/LAPACK C C THIS VERSION CONTAINS ROUTINES FROM THE BLAS/LAPACK LIBRARIES C FOR SOLVING THE SYSTEM OF LINEAR ALGEBRAIC EQUATIONS. C THIS VERSION IS FASTER THAN THE AMENDED STANDARD VERSION. C IT CONTAINS ALL THE ADDITIONAL FEATURES IN THE STANDARD VERSION. C C----------------------------------------------------------------------- C PROGRAM MAIN C INCLUDE 'PARAM' C INTEGER HEAD DIMENSION MFGK(MZ01) COMMON C(MZ08,MZ08),G(MZ08,MZ07) C COMMON /HOLD/DEL1,DEL2,DELNU,ELOW,EMAX,FINT,FRAT,HEAD,IDM,SZ,Z, C BLAM(MZ07,MZ07,5),EKLSQ(MZ01),ET(MZ07),ETERM(MZ02), C F1(2,MZ07,MZ07),F2(2,MZ07,MZ07),FI(MZ05,MZ05),FLAG(MZ09,3,4), C FX(2,MZ07,MZ07),P(MZ06,MZ03),Q(MZ06,MZ03),R(MZ06),S(MZ10,3), C SD1(MZ05,MZ05),SD2(MZ05,MZ05),SH(4,5),T(MZ10,3), C IFILE,IPRINT,ISTOT,ITC,ITCP1,ITMX,JFILE,KCASE,LAMMX,LCUT,LTOT, C MCH,MDG,MODE,N,N1,N2,NCHB,NCHCL,NCHF,NCHFMX,NCHOP,NLCTRN,NLOW, C NMX,NRAD,NTERM,NTHROW,NTOT,NTOTMX, C ISTERM(MZ02),ITERM(MZ07),LL(MZ07),LQN(MZ03),LTERM(MZ02), C MDEGEN(MZ01),MTCHCL(MZ01),NQN(MZ03) COMMON /FGCOM/R1,R2,KJ(MZ07),METHOD,MFG COMMON /OUT/IPUNCH,ICONV COMMON /FORM/JFORM COMMON /OV/SQK,IERROR,NT,NC COMMON /POL1/ALPHA,RDICUT,NPOL,NDIEL,IGC COMMON /SCRATCH/KFILE COMMON /TXT/TEXT(8) COMMON /WORK1/PP(MZ06),PP1(MZ06),PP2(MZ06),RHO(5) COMMON /WORK2/FDR(4,MZ07,MZ07) COMMON /WORK3/EL(MZ07),WF(MZ07),WF1(MZ07),PHI(MZ05) C C INITIALIZATION C ************** C REWIND 10 C WRITE DIMENSIONS WRITE(6,600)MZ01,MZ02,MZ03,MZ04,MZ05,MZ06,MZ07,MZ08,MZ22 C TIME INITIALISATION CALL TIMES(0) C MAXIMUM VALUES FOR INITIALIZATION OF ARRAYS NMX= MZ06 NTOTMX= MZ08 NCHFMX=MZ07 C COMPUTE SH (SMALL H) USING INFORMATION FROM BLOCK DATA A1=1./720. DO 101 J=1,4 DO 101 I=1,5 101 SH(J,I)=SH(J,I)*A1 C C READ - C ====== C (SEE CPC SECTION 2.1(I)) C NCASE = NUMBER OF CASES OF DIFFERENT (S,L,PARITY) C KM = NUMBER OF ENERGIES FOR EACH CASE C IFILE = UNIT FOR READING DATA IN SUBROUTINES ATOM AND C EQNS (DEFAULT - UNIT = 5) C JFILE = UNIT FOR PUNCHED OUTPUT (DEFAULT - NO PUNCHING) C USE JFILE POSITIVE IF INFORMATION ON R MATRICES C IS ONLY PUNCHED OUTPUT REQUIRED. C IF JFILE NEGATIVE THE RADIAL FUNCTIONS ARE C PUNCHED ON IABS(JFILE) C MODE = 0 FOR USE OF UNIT 10 AS SCRATCH FILE OR FOR C CREATING PERMANENT FILE FOR ONE SLPI CASE C 1 FOR READING FROM PERMANENT FILE PREVIOUSLY CREATED C UNDER MODE 0 OR 10 C 2 FOR SOLUTION OF COULOMB PROBLEM C 10 FOR CREATING PERMANENT FILE FOR MANY SLPI CASES C (NOTE THAT MODE 10 CAN BE USED WITH SCRATCH FILE C BUT THIS IS EXTRAVAGANT IN USE OF FILE SPACE) C MFG = PARAMETER FOR FOX-GOODWIN METHOD (SEE COMMENTS IN C SUBROUTINE FG). FG METHOD NOT USED FOR MFG BLANK C OR ZERO (SEE ALSO COMMENTS UNDER NOTE 1 OF ADDITIONAL C FEATURES ABOVE AND AFTER CARD IMPM 187 BELOW C NPOL = 0 FOR NO POLARIZATION C = 1 FOR CORE POLARIZATION C = 2 FOR AS 1 + DI-ELECTRONIC POLARIZATION C = 3 FOR AS 1 + EXTRA CHANNEL POLARIZATION C = 4 FOR AS 3 + DI-ELECTRONIC POLARIZATION C TEXT = INFORMATION FOR SUBSEQUENT FILE IDENTIFICATION C 2000 READ(5,501) NCASE,KM,IFILE,JFILE,JFORM,MODE,MFG,NPOL,TEXT KR=0 IF(KM.LT.0) KR=-KM KM=IABS(KM) C C WHEN ALL CASES PROCESSED PROGRAM RETURNS TO STATEMENT 2000 C TERMINATION BY BLANK CARD IF (KM*NCASE.EQ.0) STOP C C DEFAULT OPTIONS AND DIMENSION CHECKS C CHECK KM IF(KM.LE.MZ01) GO TO 110 WRITE(6,605) KM STOP 110 IF(KM.GT.0) GO TO 120 WRITE(6,625) KM STOP C C CHECK NCASE 120 IF(NCASE.GT.0) GO TO 130 WRITE(6,626) NCASE STOP C C CHECK MFG 130 IF(MFG.GE.0) GO TO 140 WRITE(6,629) MFG=0 C C CHECK MODE 140 IF(MODE.EQ.10) GO TO 150 IF(MODE.GE.0.AND.MODE.LE.2) GO TO 150 WRITE(6,630) MODE MODE=0 C C CHECK NPOL 150 IF(NPOL.GE.0.AND.NPOL.LE.4) GO TO 160 WRITE(6,631) NPOL NPOL=0 C C CHECK IFILE 160 IF(IFILE.LE.0) IFILE=5 C C IF JFILE NEGATIVE, SET IPUNCH=1 AND JFILE=IABS(JFILE) IPUNCH = 0 IF(JFILE.GE.0) GO TO 4 IPUNCH=1 JFILE=-JFILE IF(MFG.LE.0) MFG=1 L22M1 = MZ22- 1 IF(MFG.LT.L22M1) GO TO 4 WRITE (6,607) MFG,L22M1 MFG=L22M1 C UNITS 10 AND 11 ARE NOT ALLOWED FOR IJILE OR JFILE 4 IF(IFILE.EQ.10.OR.JFILE.EQ.10) GO TO 5 IF(IFILE.EQ.11.OR.JFILE.EQ.11) GO TO 5 GO TO 6 5 WRITE(6,606) STOP 6 MFGS=MFG C C ---- INSERT FOR DI-ELECTRONIC POLARIZATION C IF REQUIRED, READ DATA FOR DI-ELECTRONIC POLARIZATION C C READ - C ====== C (SEE CPC SECTION 13.4(II)) C ALPHA = CORE POLARIZABILITY C RDICUT = CUT-OFF RADIUS FOR DIELECTRONIC POLARIZATION C IGC IS SUCH THAT ORBITALS WITH (IG.LE.IGC) ARE C TREATED AS CORE ORBITALS C NDIEL=2*(NPOL/2)-NPOL+1 IF(NPOL.EQ.0) NDIEL=0 IF(NDIEL.GT.0) READ(5,503) ALPHA,RDICUT,IGC C ---- END INSERT C C READ - C ====== C (SEE CPC SECTION 2.1(II)) C FOR EACH ENERGY CASE, K=1,KM C EKLSQ(K) = KINETIC ENERGY OF ELECTRON (RYDBERGS) C WHEN TARGET SYSTEM IS IN GROUND STATE C MTCHCL(K)= CHANNEL IN WHICH CONVERGENCE TEST IS APPLIED C FOR ALL CHANNELS CLOSED C MDEGEN(K)= CHANNEL IN WHICH EIGENVALUE IS TO BE FOUND C IN THE DEGENERATE COULOMB PROBLEM C MFGK(K) = MFG FOR EACH ENERGY CASE IF REQUIRED TO BE C DIFFERENT FROM THAT READ UNDER STATEMENT C 2000 ABOVE (SEE ALSO NOTE 1 OF ADDITIONAL C FEATURES) C IF(KR.EQ.0) GO TO 990 KM=0 WRITE(6,940) 940 FORMAT(/,59H ***WARNING*** ENERGIES INPUT USING MENDOZA MODIFICA CTION. ,/,17X,14HENERGIES ARE - ,/) DO 980 IJK=1,KR READ(5,950) EKL,DELTE,NENER 950 FORMAT(2E14.6,I5) IF(NENER.LE.0) NENER=1 DO 970 LMN=1,NENER KM=KM+1 IF(KM.LE.MZ01) GO TO 960 WRITE(6,605) KM STOP 960 EKLSQ(KM)=EKL+(LMN-1)*DELTE MFGK(KM)=0 MTCHCL(KM)=0 970 MDEGEN(KM)=0 980 CONTINUE WRITE(6,945) (EKLSQ(IJK),IJK=1,KM) 945 FORMAT(8E14.6) GO TO 995 990 CONTINUE READ(5,502) (EKLSQ(K),MTCHCL(K),MDEGEN(K),MFGK(K),K=1,KM) C C READ - C ====== C (SEE CPC SECTION 2.1(III)) C IPRINT = PRINT LEVEL, DESCRIBED IN CPC APPENDIX E AND C NOTE 2 OF ADDITIONAL FEATURES C ITC = NUMBER OF POINTS IN THE STARTING REGION C FINT,EMAX,FRAT,DEL1 AND DEL2 = PARAMETERS FOR DETERMINING C THE TABULAR POINTS C DELNU = ACCURACY CRITERION FOR BOUND STATES C ITMX = MAXIMUM NUMBER OF ITERATIONS FOR BOUND STATES C 995 CONTINUE READ(5,500) IPRINT,ITC,FINT,EMAX,FRAT,DEL1,DEL2,DELNU,ITMX C C DEFAULT OPTIONS AND DIMENSION CHECKS C CHECK PAGE THROW SUPPRESSION OPTION NTHROW=0 IF(IPRINT) 1,2,7 1 NTHROW=1 IPRINT=-IPRINT GO TO 7 2 WRITE(6,616) IPRINT=-1 NTHROW=1 GO TO 8 7 IF(IPRINT.LE.5) GO TO 8 WRITE(6,617) IPRINT=5 8 IF(ITC.LE.MZ05) GO TO 10 WRITE(6,601) ITC ITC=MZ05 10 IF (ITC.GE.3) GO TO 20 WRITE(6,602) ITC ITC=3 20 ITCP1=ITC+1 IF(FINT.GE.3.) GO TO 22 IF(FINT.GT.0.) GO TO 21 WRITE(6,609) FINT=3. GO TO 22 21 WRITE(6,615) FINT 22 IF(EMAX.GE.0.) GO TO 23 EMAX=ABS(EMAX) WRITE(6,618) 23 IF(FRAT.GE.1.) GO TO 24 FRAT=1. WRITE(6,619) GO TO 25 24 IF(FRAT.GT.1.3) WRITE(6,604) FRAT 25 IF(DEL1.GE.0.03) GO TO 26 WRITE(6,620) DEL1=0.15 GO TO 27 26 IF(DEL1.LE.0.3) GO TO 27 WRITE(6,621) DEL1 27 IF(DEL2.GT.0..AND.DEL2.LE.0.1) GO TO 28 WRITE(6,603) DEL2 DEL2=0.01 28 IF(DELNU.GT.0.) GO TO 30 DELNU=DEL2 WRITE(6,622) DELNU 30 IF(ITMX.LE.0) ITMX=10 C C CHECK VALUES OF MFG FOR EACH ENERGY CASE DO 170 K=1,KM IF(MFGK(K).GE.0) GO TO 166 WRITE(6,632) K MFGK(K)=0 166 IF(IPUNCH.EQ.0) GO TO 170 IF(MFGK(K).LE.L22M1) GO TO 170 WRITE(6,634) K,MFGK(K),L22M1 MFGK(K)=L22M1 170 CONTINUE C C SUB1 - C ****** C C READ DATA FOR TARGET STATES, COMPUTE TABULAR POINTS AND TARGET C RADIAL FUNCTIONS AT THESE POINTS. COMPUTE ARRAYS C FLAG,FI,SD1,SD2,S AND T (SEE CPC SECTION 5) CALL SUB1 CALL TIMES(1) C C FOR MODE.NE.1, WRITE FIRST DATA SET ON UNIT 10 NHD=HEAD IF(MODE.EQ.1) GO TO 31 MCASE=NCASE IF(MODE.EQ.0) MCASE=1 WRITE(10) NHD,MCASE,NPOL,ITC,FINT,EMAX,FRAT,DEL1,DEL2, C IPUNCH,NCHFMX,NTOTMX GO TO 32 C FOR MODE.EQ.1 , CHECK FIRST DATA SET ON UNIT 10 31 READ(10)NHDP,NCASEP,NPOLP,ITCP,FINTP,EMAXP,FRATP,DEL1P, C DEL2P,IPNCHP,NCHFP,NTOTP IF(NHD.NE.NHDP) GO TO 320 IF(NCASE.EQ.NCASEP) GO TO 302 WRITE(6,608) NCASE,NCASEP IF(NCASE.GT.NCASEP) NCASE=NCASEP 302 IF(NPOL.NE.NPOLP) GO TO 320 IF(ITC.NE.ITCP) GO TO 320 IF(FINT.NE.FINTP) GO TO 320 IF(EMAX.NE.EMAXP) GO TO 320 IF(FRAT.NE.FRATP) GO TO 320 IF(DEL1.NE.DEL1P) GO TO 320 IF(DEL2.NE.DEL2P) GO TO 320 IF(IPUNCH.EQ.IPNCHP) GO TO 304 IF(IPUNCH.EQ.1) GO TO 320 304 IF(NCHFMX.NE.NCHFP) GO TO 320 IF(NTOTMX.NE.NTOTP) GO TO 320 GO TO 32 C CATASTROPHIC DISCREPANCY 320 WRITE(6,612) WRITE(6,613) WRITE(6,611) NHD,NCASE,NPOL,ITC,FINT,EMAX,FRAT,DEL1,DEL2, C IPUNCH,NCHFMX,NTOTMX WRITE(6,614) WRITE(6,611) NHDP,NCASEP,NPOLP,ITCP,FINTP,EMAXP,FRATP, C DEL1P,DEL2P,IPNCHP,NCHFP,NTOTP STOP C C PUNCH ON JFILE 32 IF(JFILE.EQ.0) GO TO 34 WRITE(JFILE,700) HEAD,TEXT WRITE(JFILE,701) ITC,FINT,EMAX,FRAT,DEL1,DEL2,SZ,Z,NPOL,IPUNCH WRITE(JFILE,702) NCASE,KM,NTERM WRITE(JFILE,704) (EKLSQ(K),K=1,KM) WRITE(JFILE,703) (JTERM,ISTERM(JTERM),LTERM(JTERM),ETERM(JTERM), 1 JTERM=1,NTERM ) IF(IPUNCH.EQ.0) GO TO 34 WRITE(JFILE,702) N,N2 WRITE(JFILE,710) (R(I),I=1,N2) DO 33 I=1,ITC 33 WRITE(JFILE,710)(FI(I,J),J=1,ITC) C C BEGIN DO LOOP ON CASES (S,L,PARITY) C ----------------------------------- C 34 DO 1000 ICASE=1,NCASE KCASE=ICASE IF (EKLSQ(1).GE.0.0.OR.NCASE.EQ.1) GO TO 1239 READ(5,501) KKASE,KM IF (KKASE.EQ.ICASE) GO TO 1237 PRINT 1238 PRINT 1240 STOP 1238 FORMAT (32H MISMATCH OF ICASE DATA ) 1237 CONTINUE DO 1235 K=1,KM READ (5,502) EKLSQ(K), MTCHCL(K), MDEGEN(K), MFGK(K) 1235 CONTINUE 1239 CONTINUE 1240 FORMAT (//,10X, 33H*** NOTICE OF PROGRAM CHANGE *** //, 1 5X, 38H(FOR NCASE EQ ONE OR EKLSQ(1) GT ZERO /, 2 5X, 29HTHE PROGRAM WORKS NORMALLY.) //, 3 5X, 44HTO USE DIFFERENT ENERGIES FOR EACH SLP CASE /, 4 5X, 56HDESIGNATE NCASE GT 1 AND KM = 1 ON THE FIRST IMPUT CARD;/, 5 5X, 45HTHEN USE EKLSQ(1) LT 0.0 ON THE SECOND CARD. /, 6 5X, 55HAFTER FINT CARD,INPUT A SET OF CARDS FOR EACH SLP CASE./, 7 5X, 50HTHE FIRST OF A SET GIVES ICASE, KM IN 2I5 FORMAT. /, 8 5X, 43HTHE NEXT CARDS ARE THE USUAL ENERGY CARDS. //) C C SUB2 - C ****** C C SET UP EQUATIONS WITHOUT ENERGY TERMS (CPC SECTION 6) PRINT*,EKLSQ(1) CALL SUB2 PRINT*,EKLSQ(1) CALL TIMES(2) NT=NTOT C C FOR MODE.NE.1, WRITE SECOND DATA SET ON UNIT 10 IF(MODE.EQ.1) GO TO 35 WRITE(10) ISTOT,LTOT,NTOT GO TO 36 C FOR MODE.EQ.1, CHECK SECOND DATA SET ON UNIT 10 35 READ(10) ISTOTP,LTOTP,NTOTP IF(ISTOT.EQ.ISTOTP.AND.LTOT.EQ.LTOTP.AND.NTOT.EQ.NTOTP) GO TO 36 C CATASTROPHIC DISCREPANCY WRITE(6,610) ISTOT,LTOT,NTOT,ISTOTP,LTOTP,NTOTP STOP C C SKIP IF NCHF.EQ.0 36 IF(NCHF.EQ.0) GO TO 1000 C C DEFINE SCRATCH FILE KFILE=10 IF(MODE.NE.0.AND.MODE.NE.2.AND.ICASE.GT.1) KFILE=12 C C WRITE DATA SET TO UNIT10 FOR MODE.NE.1 IF(MODE.EQ.1) GO TO 37 WRITE(10) LAMMX,BLAM CNOAA WRITE(10) C,G CALL OUTMAT(C,NTOTMX,NTOT,NTOT,10) CALL OUTMAT(G,NTOTMX,NTOT,NCHF,10) IF(KFILE.EQ.12) GO TO 371 GO TO 372 C C READ DATA SET FROM UNIT10 IF KFILE.EQ.12 37 IF(KFILE.EQ.10) GO TO 39 READ(10) LAMMX,BLAM CNOAA READ(10) C,G CALL INMAT(C,NTOTMX,NTOT,NTOT,10) CALL INMAT(G,NTOTMX,NTOT,NCHF,10) C C WRITE DATA SET TO UNIT12 IF KFILE.EQ.12 371 REWIND 12 WRITE(12) LAMMX,BLAM CNOAA WRITE(12) C,G CALL OUTMAT(C,NTOTMX,NTOT,NTOT,12) CALL OUTMAT(G,NTOTMX,NTOT,NCHF,12) C C WRITE FIRST FORM OF C TO FILES 372 IF(IPUNCH.NE.1) GO TO 375 C C READ IT FROM UNIT11 IF MODE.NE.1 IF(MODE.EQ.1) GO TO 373 REWIND 11 READ(11) C CNOAA WRITE(10) C CALL OUTMAT(C,NTOTMX,NTOT,NTOT,10) IF(KFILE.EQ.12) GO TO 374 GO TO 375 C C READ IT FROM UNIT10 IF KFILE.EQ.12 373 IF(KFILE.EQ.10) GO TO 375 CNOAA READ(10) C CALL INMAT(C,NTOTMX,NTOT,NTOT,10) C C WRITE IT TO UNIT12 IF KFILE.EQ.12 374 CONTINUE CNOAA WRITE(12) C CALL OUTMAT(C,NTOTMX,NTOT,NTOT,12) C C POSITION SCRATCH FILE AT LAMMX,BLAM 375 REWIND KFILE IF(KFILE.EQ.12) GO TO 39 READ(10) READ(10) 39 DO 1002 K=1,KM KCASE=K SQK=EKLSQ(KCASE) PRINT*, SQK,EKLSQ(1) MFG=MFGS IF(MFGK(KCASE).NE.0) MFG=MFGK(KCASE) C C CALCULATE ASYMPTOTIC FORM OF SOLUTIONS IF(IPUNCH.EQ.1) REWIND 11 READ(KFILE) LAMMX,BLAM C C FOR IPUNCH = 1 AND MODE = 1, WRITE LAMMX AND BLAM ONTO JFILE IF(K.GT.1) GO TO 44 IF(IPUNCH.NE.1) GO TO 44 IF(MODE.NE.1) GO TO 44 WRITE(JFILE,702) LAMMX DO 42 I1=1,LAMMX DO 42 I2=1,NCHF 42 WRITE(JFILE,710) (BLAM(I3,I2,I1),I3=1,NCHF) 44 CONTINUE C C SUB3 - C ****** C C CALCULATE ASYMPTOTIC SOLUTIONS (SEE CPC SECTION 7) CALL SUB3 CALL TIMES(3) NC=NCHOP C IF (IERROR.EQ.1) GO TO 1002 IF(JFILE.GT.0) WRITE(JFILE,720) EKLSQ(KCASE),KCASE,ISTOT,LTOT,HEAD C C READ C AND G MATRICES FROM DISK CNOAA READ(KFILE) C,G CALL INMAT(C,NTOTMX,NTOT,NTOT,KFILE) CALL INMAT(G,NTOTMX,NTOT,NCHF,KFILE) C C SUB4 - C ****** C C PUT ENERGY TERMS IN C AND G, SOLVE THE EQUATIONS, OBTAIN AND C PRINT FINAL RESULTS (SEE CPC SECTION 8) CALL SUB4 CALL TIMES(4) C C REPOSITION SCRATCH FILE AT LAMMX,BLAM WHEN LOOP C EXHAUSTED AND MODE.EQ.1.OR.MODE.EQ.10 IF(MODE.NE.0.AND.MODE.NE.2.AND.KCASE.EQ.KM) GO TO 1002 REWIND KFILE IF(KFILE.EQ.12) GO TO 1002 READ(10) READ(10) 1002 CONTINUE C C SKIP FIRST FORM OF C IF NOT REQUIRED CNOAA IF(MODE.EQ.1.AND.IPUNCH.NE.IPNCHP) READ(10) C IF(MODE.EQ.1.AND.IPUNCH.NE.IPNCHP) CALL INMAT(C,NTOTMX, 1NTOT,NTOT,10) C C REUSE SPACE ON UNIT10 FOR MODE.EQ.10.OR.MODE.EQ.2 IF(MODE.NE.0.AND.MODE.NE.2) GO TO 1000 REWIND 10 READ(10) 1000 CONTINUE IF(MODE.EQ.0.OR.MODE.EQ.2) REWIND 10 CALL TIMES(5) GO TO 2000 C C READ FORMATS 500 FORMAT (2I5,3F10.6,3E10.3,I2) 501 FORMAT(4I5,I2,I3,2I5,5X,8A4) 502 FORMAT(E14.6,3I5) 503 FORMAT(2F10.6,I5) C C PUNCH FORMATS 700 FORMAT(1X,24HIMPACT OF ELECTRON WITH ,A4,4X,7HTEXT = ,8A4) 701 FORMAT(I3,F5.2,F5.1,F5.2,2E9.2,2F3.0,2I3) 702 FORMAT (5I3) 703 FORMAT (3I3,E14.6) 704 FORMAT (E14.6) 710 FORMAT(5E14.7) 720 FORMAT(1X,7HEKSQ = ,E14.6,13H(INPUT EKLSQ(,I2,2H)), 1 10H FOR SLPI=,I3,1H,,I3,8H HEAD = ,A4) C C MESSAGE FORMATS 600 FORMAT(30X,'************'/30X,'** IMPACT **'/30X,'************'// X/'DIMENSIONS:'/ X' KM =',I5/ X' NTERM =',I5/ X' NRAD =',I5/ X' JMAX =',I5/ X' ITC =',I5/ X' N =',I5/ X' NCHF =',I5/ X' NTOT =',I5/ X' MFG+1 =',I5/) 601 FORMAT(30H *** WARNING *** ITC READ AS ,I3,45H WHICH IS GREATER T CHAN ALLOWED BY DIMENSIONS./,18X,29HITC PUT EQUAL TO ITCMX = MZ05) 602 FORMAT(30H *** WARNING *** ITC READ AS ,I3,50H WHICH IS LESS THAN C 3 (MINIMUM ALLOWED BY METHOD)./,18X,18HITC SET EQUAL TO 3) 603 FORMAT(31H *** WARNING *** DEL2 READ AS ,E10.3,44H WHICH IS UNSUI CTABLE. DEL2 SET EQUAL TO 0.01) 604 FORMAT(' *** WARNING *** FRAT = ',F10.6,' . FRAT LARGER THAN ', C'1.3 IS NOT RECOMMENDED') 605 FORMAT ( 31H0NUMBER OF ENERGIES TO BE RUN =,I3,52H WHICH IS GREA 1TER THAN MAXIMUM ALLOWED BY DIMENSIONS ) 606 FORMAT(47H UNITS 10 AND 11 NOT ALLOWED FOR IFILE OR JFILE) 607 FORMAT(12H MFG READ AS ,I5,11H REDUCED TO ,I5, C 55H WHICH IS MAXIMUM ALLOWED BY DIMENSIONS FOR IPUNCH.GT.0 ) 608 FORMAT(8H WARNING/1X,7(1H*)/5X, C 49HNUMBER OF SLPI CASES SPECIFIED ON UNIT 5, NCASE =,I3/5X, C 46HDIFFERS FROM NUMBER WRITTEN ON IFILE, NCASEP =,I3/2X, C 47HIF NCASE.GT.NCASEP, NCASE TAKEN EQUAL TO NCASEP/) 609 FORMAT(70H *** WARNING *** FINT READ AS NEGATIVE OR ZERO. FINT SE CT EQUAL TO 3.0) 610 FORMAT(//33H CATASTROPHIC DISCREPANCY BETWEEN/8H ISTOT =,I5, C 8H, LTOT =,I5,8H, NTOT =,I5,13H FROM INFILE/4H AND/ C 9H ISTOTP =,I5,9H, LTOTP =,I5,9H, NTOTP =,I5, C 14H FROM UNIT 10) 611 FORMAT(5X,A4,/5X,3I5/5X,3F10.6,2E14.6/5X,I5,2I10) 612 FORMAT(//43H CATASTROPHIC DISCREPANCY BETWEEN VALUES OF// C 5X,4HHEAD/5X,16HNCASE, NPOL, ITC, / C 5X,28HFINT, EMAX, FRAT, DEL1, DEL2/5X, C 22HIPUNCH, NCHFMX, NTOTMX /) 613 FORMAT(35H VALUES FROM UNITS 5 AND INFILE ARE/) 614 FORMAT(24H VALUES FROM UNIT 10 ARE /) 615 FORMAT(25H *** WARNING *** FINT = ,F10.6,41H . FINT LESS THAN 3. C0 IS NOT RECOMMENDED) 616 FORMAT(' *** WARNING *** IPRINT READ AS ZERO. IPRINT SET EQUAL ', C'TO 1 (LOWEST PRINT LEVEL)') 617 FORMAT(98H *** WARNING *** IABS(IPRINT) READ AS GREATER THAN 5. I CPRINT SET EQUAL TO 5 (HIGHEST PRINT LEVEL)) 618 FORMAT(66H *** WARNING *** EMAX READ AS NEGATIVE. EMAX CHANGED TO C ABS(EMAX)) 619 FORMAT(63H *** WARNING *** FRAT READ AS LESS THAN 1. FRAT SET EQU CAL TO 1) 620 FORMAT(' *** WARNING *** DEL1 READ AS LESS THAN 0.03. DEL1 SET ', C'EQUAL TO 0.15') 621 FORMAT(25H *** WARNING *** DEL1 = ,E10.3,44H . DEL1 GREATER THAN C 0.3 IS NOT RECOMMENDED) 622 FORMAT(75H *** WARNING *** DELNU READ AS NEGATIVE OR ZERO. DELNU CSET EQUAL TO DEL2 = ,E10.3) 625 FORMAT(/,63H ***ERROR*** NUMBER OF ENERGIES READ IN (KM) MUST BE C POSITIVE,/,15X,15HKM SPECIFIED AS,I3) 626 FORMAT(/,63H ***ERROR*** NUMBER OF CASES READ IN (NCASE) MUST BE C POSITIVE,/,15X,18HNCASE SPECIFIED AS,I3) 629 FORMAT(/,57H ***WARNING*** MFG SPECIFIED AS NEGATIVE. RESET TO Z CERO) 630 FORMAT(/,50H ***WARNING*** ILLEGAL VALUE OF MODE SPECIFIED (,I5, C17H). RESET TO ZERO. ) 631 FORMAT(/,50H ***WARNING*** ILLEGAL VALUE OF NPOL SPECIFIED (,I5, C 17H). RESET TO ZERO. ) 632 FORMAT(/,22H ***WARNING*** MFGK(,I2,39H) SPECIFIED AS NEGATIVE. CRESET TO ZERO. ) 634 FORMAT(/,22H ***WARNING*** MFGK(,I2,14H) SPECIFIED AS,I5, C11H REDUCED TO ,I5,9H WHICH IS,/,17X, C56HTHE MAXIMUM ALLOWED FOR BY DIMENSIONS FOR JFILE NEGATIVE ) END C*********************************************************************** SUBROUTINE ASYMPT(SQK,IERROR) C INCLUDE 'PARAM' C C SUBROUTINE ASYMPT ORGANISES CALCULATION OF FUNCTIONS IN THE C ASYMPTOTIC REGION, AND PUTS RESULTS OBTAINED IN FORM C REQUIRED BY IMPACT. C IT CALLS A SUBROUTINE ASYSM, WHICH IS A MODIFIED VERSION OF THE C PROGRAM ASYM OF NORCROSS, CPC,1,88,1969. C COMMON /WORK2/FX1(4,MZ07,MZ07) COMMON /WORK3/EL(MZ07),WF(MZ07),WF1(MZ07),PHI(MZ05) C C COMMON NCHD1,IJNZ(MZ25),KNZ(5),BLIN(MZ25) C ,DUMFG(MZ22,MZ07,MZ07) C INTEGER HEAD COMMON /HOLD/DEL1,DEL2,DELNU,ELOW,EMAX,FINT,FRAT,HEAD,IDM,SZ,Z, C BLAM(MZ07,MZ07,5),EKLSQ(MZ01),ET(MZ07),ETERM(MZ02), C F1(2,MZ07,MZ07),F2(2,MZ07,MZ07),FI(MZ05,MZ05),FLAG(MZ09,3,4), C FX(2,MZ07,MZ07),P(MZ06,MZ03),Q(MZ06,MZ03),R(MZ06),S(MZ10,3), C SD1(MZ05,MZ05),SD2(MZ05,MZ05),SH(4,5),T(MZ10,3), C IFILE,IPRINT,ISTOT,ITC,ITCP1,ITMX,JFILE,KCASE,LAMMX,LCUT,LTOT, C MCH,MDG,MODE,N,N1,N2,NCHB,NCHCL,NCHF,NCHFMX,NCHOP,NLCTRN,NLOW, C NMX,NRAD,NTERM,NTHROW,NTOT,NTOTMX, C ISTERM(MZ02),ITERM(MZ07),LL(MZ07),LQN(MZ03),LTERM(MZ02), C MDEGEN(MZ01),MTCHCL(MZ01),NQN(MZ03) COMMON /CMN02/ NCNL,NOP,LAMAX COMMON /CMN03/ NDL(MZ07),NDU(MZ07) COMMON /FGCOM/R1,R2,KJ(MZ07),METHOD,MFG COMMON /OUT/IPUNCH,ICONV COMMON /ASYMT/ R3,DDR C C WRITE HEADING FOR THIS ENERGY CASE IF(NTHROW.EQ.0) WRITE(6,603) SQK IF(NTHROW.EQ.1) WRITE(6,1603) SQK C C SET UP ARRAYS KNZ AND IJNZ FOR FAST ADDRESSING TECHNIQUE C (SEE CPC APPENDIX C) NCHD1=MZ07-NCHF L=0 DO 7 K=1,LAMMX IJ=0 DO 6 J=1,NCHF DO 5 I=J,NCHF IJ=IJ+1 IF(BLAM(I,J,K).EQ.0.) GO TO 5 L=L+1 IJNZ(L)=IJ 5 CONTINUE 6 IJ=IJ+NCHD1+J 7 KNZ(K)=L C C SETUP IERROR,NCHOP,DEL3,DEL4,R1,R2 AND DDR IERROR=0 NCHOP=NCHF DEL4=-.0025*SZ*SZ DEL3=0. R1=R(N1) R2=R(N2) DDR=(R2-R1)*0.25 C C INITIALISE F1 AND F2 DO 10 K=1,2 DO 10 I=1,NCHF DO 10 J=1,NCHF F1(K,I,J)=0. 10 F2(K,I,J)=0. C C CALCULATE NUMBER OF OPEN CHANNELS, NUMBER OF CLOSED CHANNELS, C AND FOR I=1,NCHF SET ET(I)= SQK-ETERM(ITERM(I)) AND EL(I)=LL(I) C (NOTING DEFINITIONS OF OPEN AND CLOSED CHANNELS GIVEN IN C CPC SECTION 7.1(II)) DO 20 I=1,NCHF EL(I)=LL(I) K=ITERM(I) E1=SQK-ETERM(K) IF(E1.GE.0.) GO TO 19 IF(E1.GE.DEL4) GO TO 15 NCHOP=NCHOP-1 GO TO 19 15 IF(E1.GE.DEL3) GO TO 19 DEL3=E1 J=I 19 IF(ABS(E1).GT.EMAX)WRITE(6,617)I,E1,EMAX 20 ET(I)=E1 C CHECK NO ET(I).EQ.0.0 FOR SZ.EQ.0 (SEE CPC SECTION 7.1 (II)) IF(SZ.GT.0.) GO TO 40 IF(NCHOP.EQ.0) GO TO 50 DO 35 I=1,NCHOP IF (ET(I).GT.0.) GO TO 35 IF(NTHROW.EQ.0) WRITE(6,601) SQK,I IF(NTHROW.EQ.1) WRITE(6,1601) SQK,I IERROR=1 RETURN 35 CONTINUE GO TO 50 C INCREMENT ET(I) FOR THRESHOLD CASE WITH SZ.GT.0 C (SEE CPC SECTION 7.1 (II)) 40 IF(DEL3.EQ.0.) GO TO 50 WRITE(6,602) DEL3,J DO 48 I=1,NCHF 48 ET(I)=ET(I)-DEL3 50 NCHCL=NCHF-NCHOP C WRITE NUMBERS OF OPEN AND CLOSED CHANNELS WRITE(6,605) NCHOP,NCHCL M=NCHOP+1 c c nrb - skip asympt but put something meaningful (sin(kr)) c also good for high energy ions. c c pih=0.5*acos(-1.0) c do i=1,nchop c ee2=sqrt(et(i)) c ee4=1.0/sqrt(ee2) c f1(1,i,i)=ee4*sin(ee2*r1-el(i)*pih) c f1(2,i,i)=ee4*sin(ee2*r2-el(i)*pih) c f2(1,i,i)=ee4*cos(ee2*r1-el(i)*pih) c f2(2,i,i)=ee4*cos(ee2*r2-el(i)*pih) c enddo c if(nchop.eq.nchf)then c return c else c stop ' can only use this test when nchop.eq.nchf' c need closed channel equivalent for use here c endif c C C R1 AND R2 FOR MFG.NE.0 IF(MFG.EQ.0.AND.NCHOP.EQ.0.AND.SZ.GT.0) MFG=1 IF(MFG.EQ.0) GO TO 51 DR=(R2-R1)*0.5 R1=R1+MFG*DR R2=R1+DR C C INTERVAL FOR INTEGRATION FROM R2 TO R1 51 MD=4 H=.25*(R1-R2) C WRITE VALUE OF MFG USED FOR THIS ENERGY IF(MFG.EQ.0) WRITE(6,618) IF(MFG.GT.0) WRITE(6,615) MFG C IF (NCHOP.EQ.0) GO TO 115 C C FIND ASYMPTOTIC SOLUTIONS AT R2 FOR NCHOP GREATER THAN ZERO R3=R2 CALL ASYSM DO 60 I=1,NCHF DO 70 J=1,NCHOP DO 70 K=1,2 70 F2(K,I,J)=FX(K,I,J) IF (NCHCL.EQ.0) GO TO 60 DO 80 J=M,NCHF 80 F2(1,I,J)=FX(1,I,J) 60 CONTINUE C FIND SOLUTIONS AT R1 BY INTEGRATING INWARDS FROM R2 MM=2 DO 90 J=1,NCHF IF (J.GT.NCHOP) MM=1 DO 90 K=1,MM DO 100 I=1,NCHF WF(I)=FX(K,I,J) 100 WF1(I)=FX1(K,I,J) RS=R2 CALL SOLV(RS,H,MD,WF,WF1) DO 110 I=1,NCHF 110 F1(K,I,J)=WF(I) 90 CONTINUE GO TO 190 C C CASE OF ALL CHANNELS CLOSED (SEE CPC SECTION 7.4) C CALCULATE DEL1, THE INCREMENT IN ENERGY 115 C2=5.*DELNU IF(C2.GT.0.01) C2=0.01 IF(C2.LT.0.001) C2=0.001 A1=-ET(1) IF(SZ.NE.0) GO TO 116 DEL1=A1*C2 GO TO 117 116 DEL1=2.*A1*C2*SQRT(A1)/SZ C FIND ASYMPTOTIC SOLUTIONS FOR ENERGIES SQK AND SQK + DEL1 117 DO 170 K=1,2 IF(K.EQ.1) GO TO 130 DO 120 I=1,NCHF 120 ET(I)=ET(I)+DEL1 C SOLUTIONS AT R2 130 R3=R2 CALL ASYSM DO 140 I=1,NCHF DO 140 J=1,NCHF 140 F2(K,I,J)=FX(1,I,J) C SOLUTIONS AT R1 DO 160 J=1,NCHF DO 150 I=1,NCHF WF(I)=FX(1,I,J) 150 WF1(I)=FX1(1,I,J) RS=R2 CALL SOLV(RS,H,MD,WF,WF1) DO 160 I=1,NCHF 160 F1(K,I,J)=WF(I) 170 CONTINUE DO 180 I=1,NCHF 180 ET(I)=ET(I)-DEL1 C C PARAMETER NFG DETERMINES NORMALISATION FOR CLOSED CHANNELS 190 NFG=0 IF(NCHOP.GT.0.OR.SZ.EQ.0) NFG=1 C C USE OF FOX-GOODWIN INTEGRATION TECHNIQUE IF(MFG.EQ.0) GO TO 193 C COMPUTE ARRAY KJ(J) WHICH DETERMINES METHOD OF NORMALISATION C FOR ALL CHANNELS CLOSED AND SZ GREATER THAN 0 C (SEE CPC SECTION 10.3(III)) C C INITIALISE KJ(J) = 1 DO 183 J=1,NCHF 183 KJ(J)=1 IF(NFG.EQ.1) GO TO 192 A1=1./R(N1) C1=2.*SZ*A1 C2=A1*A1 DO 191 J=1,NCHF FL=LL(J) C COMPUTE WJ=2*SZ/RN1-FL*(FL+1)/RN1**2+ET(J) WJ=C1-FL*(FL+1.)*C2+ET(J) C IF WJ.GT.0 PUT KJ(J)=2 191 IF(WJ.GT.0) KJ(J)=2 C 192 CONTINUE CALL FG C C RENORMALIZE CLOSED CHANNEL FUNCTIONS 193 IF(NCHCL.EQ.0) GO TO 440 M=NCHOP+1 A5=SQRT(SZ/2.) DO 290 J=M,NCHF J1=NDL(J) J2=NDU(J) E1=-ET(J) FL=LL(J) DO 280 K=1,2 IF(NFG.EQ.1) GO TO 240 C SCAN FOR DEGENERATE GROUP OF CLOSED CHANNELS. DO 200 L=J1,J2 IF(KJ(L).EQ.2) GO TO 210 200 CONTINUE GO TO 240 C FOR KJ(J)=2 COMPUTE NORMALISING FACTOR C - A1=SQRT(SZ/2.)*(FNU**(FL+.5))/(GAMMA(FNU+FL+1.) 210 A1=SZ/SQRT(E1) A2=A1**(FL+.5) A1=A1+FL NX=A1 X=NX X=A1-X B1=1. IF(NX.EQ.0) GO TO 230 DO 220 L=1,NX B2=L-1 220 B1=B1*(A1-B2) C COMPUTE B2=GAMMA(1+X) FOR 0.LE.X.LE.1 230 B2=1.+X*(-.5748646+X*(.9512363+X*(-.6998588+X*(.4245549 C -.1010678*X)))) A1=B1*B2 A1=A5*A2/A1 GO TO 260 C FOR KJ(J)=1 COMPUTE NORMALISING FACTOR C - A1=1./SQRT(SUM OVER I OF (F1(K,I,J)**2+F2(K,I,J)**2)) 240 A1=0. DO 250 I=1,NCHF 250 A1=A1+F1(K,I,J)**2+F2(K,I,J)**2 A1=1./SQRT(A1) C MULTIPLY FUNCTIONS BY NORMALISING FACTOR A1 260 DO 270 I=1,NCHF F1(K,I,J)=A1*F1(K,I,J) 270 F2(K,I,J)=A1*F2(K,I,J) IF(NCHOP.GT.0) GO TO 290 280 E1=E1-DEL1 290 CONTINUE C C RENORMALIZE OPEN CHANNEL FUNCTIONS FOR SZ.GT.1 440 IF(SZ.LE.1) GO TO 435 IF(NCHOP.EQ.0) GO TO 435 A1=SQRT(SZ) DO 430 K=1,2 DO 430 I=1,NCHF DO 430 J=1,NCHOP F1(K,I,J)=F1(K,I,J)*A1 430 F2(K,I,J)=F2(K,I,J)*A1 C C PRINT FUNCTIONS AT MATCHING POINTS IF IPRINT.GE.4 435 IF(IPRINT.LT.4) RETURN IF(NTHROW.EQ.0) WRITE(6,610) IF(NTHROW.EQ.1) WRITE(6,1610) IF(MFG.GT.0.AND.NCHOP.EQ.0) WRITE(6,616) METHOD DO 442 NPI=1,2 WRITE(6,611) NPI M3=NCHF IF(NCHOP.GT.0) GO TO 436 WRITE(6,620) GO TO 437 436 IF(NCHCL.EQ.0) WRITE(6,612) IF(NCHCL.NE.0) WRITE(6,622) 437 M1=(M3-1)/10+1 M2=1 DO 4438 I1=1,M1 IF(I1.GT.1) WRITE(6,619) M=M2+9 IF(M.GT.NCHF) M= NCHF DO 438 I=1,NCHF IF(NPI.EQ.1) WRITE(6,614) (F1(1,I,J),J=M2,M) 438 IF(NPI.EQ.2) WRITE(6,614) (F2(1,I,J),J=M2,M) 4438 M2=M2+10 IF(NCHOP.GT.0) GO TO 439 WRITE(6,621) GO TO 443 439 M1=(NCHOP-1)/10+1 M3=NCHOP IF(NCHCL.EQ.0) WRITE(6,613) IF(NCHCL.NE.0) WRITE(6,623) 443 M2=1 DO 4441 I1=1,M1 M=M2+9 IF(I1.GT.1) WRITE(6,619) IF(M.GT.M3) M=M3 DO 441 I=1,NCHF IF(NPI.EQ.1) WRITE(6,614) (F1(2,I,J),J=M2,M) 441 IF(NPI.EQ.2) WRITE(6,614) (F2(2,I,J),J=M2,M) 4441 M2=M2+10 442 CONTINUE RETURN C C FORMATS FOR MESSAGES 601 FORMAT(14H1 FOR EKLSQ = ,E14.6,6H , ET(,I2,6H) =0.0/ C 46H THIS CASE NOT ALLOWED FOR SZ=0, RUN BYPASSED ) 1601 FORMAT(/////,14H FOR EKLSQ = ,E14.6,6H , ET(,I2,6H) =0.0/, C46H THIS CASE NOT ALLOWED FOR SZ=0, RUN BYPASSED ) 602 FORMAT(30H ENERGIES INCREMENTED BY DEL3=,E14.6,12H TO GIVE ET(,I2, C 5H)=0.0,/) 617 FORMAT( 54H *** WARNING *** IT IS RECOMMENDED TO HAVE ABS(ET(I)), C 18H.LE.EMAX FOR ALL I/,17X,32H THIS CONDITION IS VIOLATED (ET(, C I2,2H)=,E12.4,10H AND EMAX=,E12.4, C 38H) AND NUMERICAL ERRORS MAY BE PRODUCED,/) C PRINT FORMATS 603 FORMAT(1H1,23HRESULTS FOR EKLSQ = ,F12.6,53H (ENERGY IN RYDBER CGS RELATIVE TO LOWEST TARGET TERM),/,1X,35(1H*)//) 1603 FORMAT(/////,24H RESULTS FOR EKLSQ = ,F12.6,53H (ENERGY IN RYD CBERGS RELATIVE TO LOWEST TARGET TERM),/,1X,35(1H*)//) 605 FORMAT(28H NUMBER OF OPEN CHANNELS =,I2/ 1 28H NUMBER OF CLOSED CHANNELS =,I2/) 610 FORMAT(1H1,28HRADIAL FUNCTIONS FROM ASYMPT,16H (COLUMN J FOR , C 42HBOUNDARY CONDITION J, ROW I FOR CHANNEL I)/1X,28(1H*)) 1610 FORMAT(/////,29H RADIAL FUNCTIONS FROM ASYMPT,16H (COLUMN J FOR , C42HBOUNDARY CONDITION J, ROW I FOR CHANNEL I)/1X,28(1H*)) 611 FORMAT(/5X,17HSOLUTIONS AT R(N+,I1,1H)/5X,19(1H-)) 612 FORMAT(/,10X,13HSIN SOLUTIONS,/) 613 FORMAT(/,10X,13HCOS SOLUTIONS,/) 619 FORMAT(1H0) 620 FORMAT(/,10X,41HDECAYING SOLUTIONS AT ENERGY EKLSQ(KCASE),/) 621 FORMAT(/,10X,40HDECAYING SOLUTIONS AT INCREMENTED ENERGY,/) 622 FORMAT(/,10X,71HSIN SOLUTIONS FOR OPEN CHANNELS, DECAYING SOLUTION CS FOR CLOSED CHANNELS,/) 623 FORMAT(/,10X,31HCOS SOLUTIONS FOR OPEN CHANNELS,/) 614 FORMAT(10E12.4) 615 FORMAT(40H FOX-GOODWIN INTEGRATION USED WITH MFG = ,I5) 616 FORMAT(48H FUNCTIONS NORMALISED USING FOX-GOODWIN METHOD = ,I2) 618 FORMAT(40H FOX-GOODWIN INTEGRATION METHOD NOT USED) END C********************************************************************** SUBROUTINE ASYSM C INCLUDE 'PARAM' C C THIS IS A MODIFIED VERSION OF THE PROGRAM ASYM OF NORCROSS, COMP. C PHYS. COMM. 1 (1969) 88. MODIFICATIONS ARE DESCRIBED IN IMPACT C CPC PAPER, SECTION 7 AND APPENDIX C, AND ALSO ON COMMENT CARDS C BELOW AND IN SECTION 9 OF WRITE-UP OF PROGRAM ASYPCK (CREES, C COMPUT. PHYS. COMMUN., SUBMITTED JUNE 1979) C C --------------------------------------------------------------------- C C INPUT DATA - C ************ C C ET(N)=ENERGIES (RYDBERG UNITS) C ALL OPEN CHANNELS MUST PRECEED ALL CLOSED CHANNELS C DEGENERATE CHANNELS MUST BE GROUPED TOGETHER. C DEGENERACEY ASSUMED IF ENERGIES AGREE FRACTIONALLY TO 1.D-4. C EL(N)=ANGULAR MOMENTUM VALUES (ONLY USED FOR INPHS=.TRUE.) C MUST ALWAYS BE INCLUDED IN ARRAY CN, I.E. CN(I,I,1)=-L*(L+1) C N=NUMBER OF CHANNELS (LESS THAN OR EQUAL TO NCHFMX) C SZ=RESIDUAL CHARGE (GREATER THAN OR EQUAL TO 0) C CN=COEFFICIENTS OF (1/R)**(K+1). A SYMMETRIC MATRIX C POSITIVE/NEGATIVE FOR ATTRACTIVE/REPULSIVE POTENTIALS C M=MAXIMUM VALUE OF K IN ABOVE SERIES (LESS THAN OR EQUAL TO 5) C RO=VALUE OF R (ATOMIC UNITS) AT WHICH SOLUTION DESIRED) C DR=INITIAL ESTIMATE OF STEP-LENGTH FOR INWARD INTEGRATIONS C C --------------------------------------------------------------------- C C PARAMETERS SET-UP INTERNALLY - C ****************************** C C NZD=RESIDUAL CHARGE (INTEGER SZ) C DELE=SOME SMALL POSITIVE ENERGY. TYPICALLY 0.07*SZ*SZ C FOR DETAILS SEE NORCROSS,CPC,1,88,1969 C EROR=CRITERION FOR CONVERGENCE,ENERGY DEGENERACY AND INTEGRATION C MUST BE GREATER THAN OR EQUAL TO 1.E-5 C IBUG=0 FOR NO PRINT OUT =1 FOR SOME PRINT OUT C SET EQUAL TO 0 IF IPRINT =1,2,3 OR 4 C SET EQUAL TO 1 IF IPRINT = 5 C WARNING MESSAGES PRINTED WITH IBUG = 0 OR 1 C IOUT=SYSTEM I/O DEVICE NUMBER FOR WRITE STATEMENTS C C --------------------------------------------------------------------- C C OUTPUT RESULTS - C **************** C C F(K,I,J) IS I-TH MEMBER OF J-TH BOUNDARY CONDITION C (ORDERED AS INPUT ET(J)) C C F(1,I,J) = REGULAR SOLUTIONS C F(2,I,J) = IRREGULAR SOLUTIONS C FDR(1,I,J) = DERIVATIVES OF REGULAR SOLUTIONS C FDR(2,I,J) = DERIVATIVES OF IRREGULAR SOLUTIONS C C FOR OPEN CHANNELS C F(1,I,I) ASYMPTOTIC TO (ET(I)**(-1/4))*SIN(PHI) C F(2,I,I) ASYMPTOTIC TO (ET(I)**(-1/4))*COS(PHI) C C FOR CLOSED CHANNELS C F(1,I,I) ASYMPTOTIC TO (ET(I)**(-1/4)*EXP(-PHI) C FOR ET(I) = 0.0, ET(I) ARE REPLACED IN C ASYMPTOTIC AMPLITUDE BY 2.0*SZ/R C F(2,I,I) ASYMPTOTIC SOLUTIONS NOT CALCULATED C C C --------------------------------------------------------------------- C C INTEGER HEAD COMMON NCHD1,IJNZ(MZ25),KNZ(5),BLIN(MZ25),NZD,NZDM, C AB1(MZ07,13,10),AB2(MZ07,13,10),ALP(13,13),CC(MZ07,MZ07,5), C E2(MZ07),PP(12,10),T1(MZ07),V(MZ11,11,10),ZET(12,10),ZTIN(12,10) C ,GAM1(MZ07,15),GAM2(MZ07,15) C , VF1(MZ07),VF2(MZ07),VFDR1(MZ07),VFDR2(MZ07),Z1,Z2 COMMON /CMN02/ NCNL,NOP,LAMAX COMMON /CMN03/ NDL(MZ07), NDU(MZ07) COMMON /CMN04/ HY,NHY,JJ COMMON /ASYMT/ RO,DR COMMON /WORK2/FDR(4,MZ07,MZ07) COMMON /WORK3/EL(MZ07),WF(MZ07),WF1(MZ07),PHI(MZ05) COMMON /HOLD/DEL1,DEL2,DELNU,ELOW,EMAX,FINT,FRAT,HEAD,IDM,SZ,Z, C CN (MZ07,MZ07,5),EKLSQ(MZ01),ET(MZ07),ETERM(MZ02), C FC(2,MZ07,MZ07),F2(2,MZ07,MZ07),FI(MZ05,MZ05),FLAG(MZ09,3,4), C F (2,MZ07,MZ07),P(MZ06,MZ03),Q(MZ06,MZ03),R(MZ06),S(MZ10,3), C SD1(MZ05,MZ05),SD2(MZ05,MZ05),SH(4,5),T(MZ10,3), C IFILE,IPRINT,ISTOT,ITC,ITCP1,ITMX,JFILE,KCASE,M ,LCUT,LTOT, C MCH,MDG,MODE,NNN,N1,N2,NCHB,NCHCL,N,NCHFMX,NCHOP,NLCTRN,NLOW, C NMX,NRAD,NTERM,NTHROW,NTOT,NTOTMX, C ISTERM(MZ02),ITERM(MZ07),LZ(MZ07),LQN(MZ03),LTERM(MZ02), C MDEGEN(MZ01),MTCHCL(MZ01),NQN(MZ03) DIMENSION WFD(MZ07),WFD1(MZ07),ISF(MZ07) C C --------------------------------------------------------------------- C C COMMENTS - C ********** C C NOT APPLICABLE FOR C 1) SZ = 0 IF ANY ET(J) = 0 C 2) SZ LESS THAN 0 C RO INCREASED TO RF UNTIL SOLUTION ACCURATE TO 1.E-4 C OBTAINED. RESULT THEN INTEGRATED BACK TO RO C LARGE RF MAY BE REQUIRED (WARNING PRINTED) FOR CASES C 1) LARGE POSITIVE ENERGIES JUST BELOW THRESHOLD C 2) SZ = 0 JUST ABOVE THRESHOLD C 3) NEARLY DEGENERATE CHANNELS C ISF(I) IS THE NUMBER OF TERMS TAKEN IN ASYMPTOTIC C SERIES FOR CHANNEL I C USES INTERNAL SUBROUTINE EXPAN,ITERA,POTS,ZETA,PPFS,QROP, C POP,DCHAIN,PHAS,BOOLE,COULGM C AND SUBROUTINES SOLV,WOUTER C C --------------------------------------------------------------------- C C CHANGES MADE FROM ORIGINAL ASYM - C ********************************* C C THE SUBROUTINE ASYSM (AND SUBROUTINES CALLED BY IT) DIFFERS C FROM ASYM OF NORCROSS (AND SUBROUTINES CALLED BY IT) AS C FOLLOWS (SEE CPC APPENDIX C) - C 1) MAJOR RECODING TO IMPROVE EFFICIENCY C 2) INSERTION OF COMMON BLOCKS HOLD AND WORK C 3) MORE USE OF BLANK COMMON AS WORKING SPACE C 4) SUBROUTINES VMUL AND FACT CALLED BY ASYM HAVE BEEN DELETED C MZ05) THE INPHS OPTION OF ASYM HAS BEEN DELETED C MZ01) DIFFERENT METHOD USED FOR CALCULATION OF PHASE C OF FUNCTIONS C C NO CHANGES HAVE BEEN MADE IN NUMERICAL METHODS, EXCEPT AS C MENTIONED IN CPC SECTION 7.1, APPENDIX C(X), IN SUBROUTINES PHAS C AND COULGM, AND ON COMMENT CARDS BELOW. C C --------------------------------------------------------------------- C C CHANGE FROM THE PUBLISHED PROGRAM IN CRITERIA FOR CALLING ITERA C *************************************************************** C C FOR THE CASE OF Z.GT.0 (POSITIVE IONS) AND WHERE THE CHANNEL ENERGY, C E2J, LIES JUST ABOVE A THRESHOLD (I.E. E2J LIES IN THE RANGE C (0.LE.E2J.LE.DELE), WHERE DELE = 0.07 * Z * Z), IN THE PUBLISHED C VERSION OF THE PROGRAM, IMPACT TRIES BOTH THE EXPAN AND ITERA C METHODS OF SOLUTION AND USES SOLUTIONS FROM THE METHOD GIVING THE C VALUE OF RF CLOSEST TO RO (I.E. WHERE THE SOLUTIONS ARE REQUIRED) C SO THAT EXPAN IS TRIED FIRST AND IF RFE IS NOT GREATER THAN C -- SEE NORCROSS, CPC, 1, (1969), 88. THE PROCEDURE HAS BEEN CHANGED C RO * 1.15, AND PROVIDING THAT TAKING THIS VALUE FOR RF IMPLIES C THAT SOLV WOULD REQUIRE LESS THAN 250 INTEGRATION POINTS TO C INTEGRATE INWARDS FROM RF TO RO, IMPACT WILL NOW TAKE THE EXPAN C SOLUTIONS WITHOUT TRYING THE ITERA METHOD. THIS CAN AVOID C UNNECESSARY CALLS TO ITERA, AND CONSEQUENTLY CAN LEAD TO A SUB- C STANTIAL REDUCTION IN CPU TIME, WITHOUT AFFECTING THE ACCURACY OF C THE SOLUTIONS. FULLER DETAILS ARE GIVEN IN SECTION 9 OF C CREES, COMPUT. PHYS. COMMUN., SUBMITTED JUNE 1979 C C ---------------------------------------------------------------------- C C INITIALIZATION C -------------- C DATA PI,IAS,NIT,ALZ/3.1415926535898,15,6,1.0/ NHY=10 RST=0.0 C C SET UP IBUG,IOUT,NZD,DELE AND EROR AND CHECK INPUT DATA C PUT NZD = INTEGER SZ (ALLOWING FOR POSSIBLE ROUNDING ERRORS) NZD=SZ+0.5 IOUT=6 IF(N.GT.MZ07) GO TO 2 IF(M.GT.5) GO TO 2 IF(NZD.GE.0) GO TO 3 2 WRITE(IOUT,290) N,M,NZD CALL EXIT(0) 3 NCNL=N LAMAX=M IBUG=0 IF(IPRINT.EQ.5) IBUG=1 DELE=0.07*SZ*SZ EROR=DEL2*0.2 IF(EROR.LE.1.E-5) EROR=1.E-5 C C SETUP RFT AND RFU C ================= RFT=RO*1.15 RFU=RO+250*DR IF(IBUG.EQ.0) GO TO 4 WRITE(IOUT,300) DO 303 K=1,LAMAX WRITE(IOUT,301) K DO 303 I=1,NCNL WRITE(IOUT,801) (CN(I,J,K),J=1,NCNL) 303 CONTINUE C C CHANGE TO KAPPA-RHO UNITS C DEFINE NOP=NUMBER OF OPEN CHANNELS C DEFINE UPPER AND LOWER CHANNEL NUMBERS FOR ENERGY DEGENERACY C 4 NOP=0 IF (NZD.NE.0 ) ALZ=NZD ALZD=1./ALZ ALZ2=ALZD*ALZD IF(NCHOP.NE.NCNL) AETN=LOG(2./EROR)/SQRT(ABS(ET(NCNL))) DO 10 I=1,NCNL E2(I)=ET(I)*ALZ2 IF(E2(I).LT.0.0) GO TO 5 NOP=I 5 DO 10 J=1,NCNL AB12=ALZD DO 10 K=1,LAMAX AB12=AB12*ALZ 10 CC(I,J,K)=CN(I,J,K)*AB12 DO 70 J=1,NCNL E2J=E2(J) DO 40 I=1,J IF(ABS(ET(I)-ET(J)).NE.0.0) GO TO 40 NDL(J)=I GO TO 50 40 CONTINUE 50 DO 60 I=J,NCNL IF(ABS(ET(I)-ET(J)).NE.0.0) GO TO 70 60 NDU(J)=I 70 CONTINUE C C SOLUTIONS FOR EACH BOUNDARY CONDITION J C --------------------------------------- C DO 250 J=1,NCNL C C INITIALIZATION E2J=E2(J) ETJ=ET(J) IF(IBUG.EQ.1) WRITE(IOUT,600) JJ=J RFI=1.0E+8 RFE=1.0E+8 C C TEST FOR CASE OF ZERO ENERGY C ============================ IF(E2J.EQ.0.0) GO TO 76 C C TRY ASYMPTOTIC EXPANSION METHOD C =============================== ZA1=NZD RS=AMAX1(RO,ABS(2.0*ZA1/ETJ)) IF(ETJ.LT.0) RS=RO CALL EXPAN(IAS,EROR,RS,ISF,RFE,IBUG,IOUT) C C CHECK FOR POSSIBLE USE OF ITERATIVE WBK METHOD C ============================================== IF(NZD.EQ.0) GO TO 86 IF(J.GT.NOP) GO TO 86 IF(ETJ.GT.DELE) GO TO 86 IF(RFE.LT.RFT.AND.RFE.LE.RFU) GO TO 86 C C TRY WBK-ITERATION SOLUTION C -------------------------- C COMPUTE RG, APPROXIMATELY EQUAL TO TWICE VALUE OF R FOR C FIRST POINT OF INFLECTION IN COULOMB FUNCTION 76 CNJJ=-CN(J,J,1) E2CN=E2J*CNJJ IF(E2CN.LT.0.01) GO TO 77 RG=(SQRT(1.+E2CN)-1.)*2.*ALZ/ETJ GO TO 78 77 RG=CNJJ*ALZD C DEFINE RS AS LARGER OF RG, R0 78 RS=AMAX1(RO,RG) IF(RS.EQ.RST) GO TO 79 C SET UP MESH AND POTENTIAL MATRIX FOR USE IN SUBROUTINE ITERA HY=1.0/(FLOAT(NHY)*SQRT(RS*ALZ)) CALL POTS(2*NIT-1) 79 RST=RS CALL ITERA(NIT,EROR/10.0,RS,RFI,IBUG,IOUT) C C DETERMINATION OF METHOD BY CRITERION OF SMALLEST RF C --------------------------------------------------- C 85 RF=MIN(RFI,RFE) IF(RFI.LE.RFE) GO TO 130 C C EXPANSION METHOD USED C --------------------- C COMPUTE AMPLITUDE FACTORS 86 RF=RFE IF(IBUG.EQ.1) WRITE(IOUT,400) J RHO=RF*ALZ IF(J.GT.NOP) GO TO 93 RZZ=1./RHO DO 92 I=1,NCNL ISG=ISF(I) ISG1=ISG-1 AB11=GAM1(I,1) AB12=0. AB21=GAM2(I,1) AB22=0. A1=1. F1=0. DO 90 LL=2,ISG1 F1=F1+RZZ A1=A1*RZZ TERM1=GAM1(I,LL)*A1 TERM2=GAM2(I,LL)*A1 AB11=AB11+TERM1 AB12=AB12-TERM1*F1 AB21=AB21+TERM2 90 AB22=AB22-TERM2*F1 F1=F1+RZZ A1=A1*RZZ*0.5 TERM1=GAM1(I,ISG)*A1 TERM2=GAM2(I,ISG)*A1 VF1(I)=AB11+TERM1 VF2(I)=AB21+TERM2 VFDR1(I)=AB12-TERM1*F1 92 VFDR2(I)=AB22-TERM2*F1 93 XKEY=0. EKAP=SQRT(ABS(E2J)) SGN=NZD IF (J.GT.NOP) SGN=-NZD Z1=EKAP+SGN/(SQRT(ABS(ETJ))*RHO) Z2=0. SQZ=SQRT(EKAP*ALZ) GO TO 135 C C ITERATION METHOD USED C --------------------- 130 XKEY=1.0 IF(IBUG.EQ.1) WRITE(IOUT,500) J SQZ=SQRT(Z1*ALZ) C C WAVE FUNCTIONS AT RF C -------------------- C COMPUTE PHASE (SEE CPC APPANDIX C SECTION (X)) 135 IF(NZD.GT.0) GO TO 136 C CASE OF Z = 0 PH=RF*ALZ*SQRT(ABS(E2J)) IF(J.GT.NOP) GO TO 141 PH=PH-EL(J)*0.5*PI GO TO 139 136 IF(E2J.NE.0.) GO TO 137 C CASE OF Z .GT. 0 , K**2 = 0 PH=SQRT(8.*RF*ALZ)-(EL(J)+0.25)*PI GO TO 139 137 IF(E2J.GT.0) GO TO 138 C CASE OF Z .GT. 0 , K**2 .LT. 0 AG=SQRT(ABS(E2J)) RH=RF*ALZ*AG PH=RH-LOG(RH+RH)/AG GO TO 141 C CASE OF Z .GT. 0 , K**2 .GT. 0 138 PH=PHAS(RF*ALZ,E2J,XKEY) AG=EL(J) CALL COULGM(AG,E2J,PI,DG) PH=PH+DG 139 CONTINUE C C COMPLETE COMPUATION OF FUNCTIONS AT RF SI=SIN(PH)/SQZ CO=COS(PH)/SQZ DO 140 I=1,NCNL AB11=VF1(I) AB12=VFDR1(I) AB21=VF2(I) AB22=VFDR2(I) WF(I)=SI*AB11+AB21*CO WF1(I)=CO*AB11-AB21*SI TA=AB12-AB21*Z1-Z2*AB11 TB=AB22+AB11*Z1-Z2*AB21 WFD(I)=ALZ*(TA*SI+TB*CO) 140 WFD1(I)=ALZ*(TA*CO-TB*SI) GO TO 150 141 EXD=EXP(-PH)/SQZ AB11=Z1+Z2 DO 145 I=1,NCNL WF(I)=VF1(I)*EXD 145 WFD(I)=((VFDR1(I)*EXD-WF(I)*AB11))*ALZ C C INTEGRATION FROM RF TO DESIRED R0 C --------------------------------- 150 IF(RF.GT.RO) GO TO 149 DO 146 I=1,NCNL F(1,I,J)=WF(I) 146 FDR(1,I,J)=WFD(I) IF(J.GT.NOP) GO TO 250 DO 147 I=1,NCNL F(2,I,J)=WF1(I) 147 FDR(2,I,J)=WFD1(I) GO TO 250 149 NHR=1+ABS((RF-RO))/DR ! ABS - NRB C IF(RO-RF.GT.DR)WRITE(6,731)RF,RO,J ! NRB C C CHECK FOR LARGE NUMBER OF INTEGRATION POINTS IF (FLOAT(NHR).LE.EROR*2.5E6) GO TO 720 WRITE(IOUT,710) NHR,J 720 IF(NOP.EQ.NCNL) GO TO 750 IF(J.GT.NOP) GO TO 750 A1=RO+AETN C C CHECK FOR LARGE VALUE OF RF IF (RF.GT.A1) WRITE(IOUT,730) RF,J 750 HR=(RO-RF)/FLOAT(NHR) IF(IBUG.EQ.1) WRITE(IOUT,700) RF,NHR,HR,J C C CALL SOLV FOR INWARD INTEGRATIONS RA=RF CALL SOLV(RA,HR,NHR,WF,WFD) DO 235 I=1,NCNL F(1,I,J)=WF(I) 235 FDR(1,I,J)=WFD(I) IF(J.GT.NOP) GO TO 240 RA=RF CALL SOLV(RA,HR,NHR,WF1,WFD1) 240 DO 245 I=1,NCNL F(2,I,J)=WF1(I) 245 FDR(2,I,J)=WFD1(I) C 250 CONTINUE C COMPLETES LOOP ON SOLUTIONS FOR DIFFERENT BOUNDARY CONDITIONS C ------------------------------------------------------------- C C DE-BUG PRINTING C --------------- IF(IBUG.EQ.0) RETURN WRITE(IOUT,800) RO,(I,I=1,NCNL) WRITE(IOUT,798) M3=NCNL DO 810 K=1,2 M1=(M3-1)/10+1 DO 808 KK=1,2 IF(KK.EQ.2) WRITE(IOUT,797) M2=1 DO 806 I1=1,M1 MCNL=M2+9 IF(MCNL.GT.M3) MCNL=M3 IF(I1.GT.1) WRITE(IOUT,796) DO 804 I=1,NCNL IF(KK.EQ.1) WRITE(IOUT,801) (F(K,I,J),J=M2,MCNL) 804 IF(KK.EQ.2) WRITE(IOUT,801) (FDR(K,I,J),J=M2,MCNL) 806 M2=M2+10 808 CONTINUE IF(K.EQ.2) GO TO 260 IF(NOP.EQ.0) GO TO 260 WRITE(IOUT,799) 810 M3=NOP C 260 WRITE(IOUT,900) RETURN C C PRINT FORMATS 290 FORMAT(/,26H INPUT DATA ERROR IN ASYSM,3I5) 300 FORMAT(/38H *** ASYSM PRINT-OUT WITH IBUG = 1 ***/) 301 FORMAT(27H0COEFFICIENTS OF 1.0/R**(1+,I1,1H)) 400 FORMAT(46H USED EXPANSION METHOD FOR SOLUTION IN CHANNEL ,I3) 500 FORMAT(46H USED ITERATION METHOD FOR SOLUTION IN CHANNEL ,I3) 600 FORMAT(//) 700 FORMAT(20H INTEGRALS FROM RF = ,F9.3,8H USING A ,I6, C 30H POINT MESH WITH STEP LENGTH =,F8.5,12H IN CHANNEL ,I2) 710 FORMAT(38H *** WARNING *** INTEGRATION REQUIRED,I5,19H POINTS IN 1CHANNEL ,I3,36H , RESULTS PRODUCED MAY BE INCORRECT /17X,67H DUE T 2O INACCURACIES ARISING FROM LARGE NUMBER OF INTEGRATION STEPS ) 730 FORMAT(/,78H *** WARNING *** POSSIBLE CLOSED CHANNEL INSTABILITY 1ON INTEGRATION FROM RF =,F10.5,11H IN CHANNEL,I2,/,18X,86H(SEE NOR 2CROSS, J.PHYS.B,4,1458,1971 , AND NORCROSS AND SEATON, J.PHYS.B, 36,614,1973 )/,18X,61HUSE OF FOX-GOODWIN METHOD IS RECOMMENDED (SEE 4 CPC SECTION 10)) 731 FORMAT(' ***WARNING: OUTWARD INTEGRATION FROM',F9.3,' TO',F9.3, X' FOR CHANNEL',I3,' TRY INCREASING EMAX?') 796 FORMAT(/) 797 FORMAT(12H0DERIVATIVES) 798 FORMAT(18H0REGULAR SOLUTIONS) 799 FORMAT(20H0IRREGULAR SOLUTIONS) 800 FORMAT(//,23H WAVE FUNCTIONS AT RO =,F12.5,24H FOR BOUNDARY CONDIT 1IONS,/,(I8,9I12)) 801 FORMAT(10E12.4) 900 FORMAT(/31H *** END OF ASYSM PRINT-OUT ***) END C*********************************************************************** SUBROUTINE ATOM C INCLUDE 'PARAM' C DESCRIBED IN CPC SECTIONS 2.2 AND 5.1 C C INTEGER HEAD COMMON D(4),FL(4),FL1(4),NCARDM(MZ03),PIN(MZ04,MZ03),PHI(MZ04), CQIN(MZ04,MZ03),RIN(MZ04),ZEFF(MZ04) COMMON /HOLD/DEL1,DEL2,DELNU,ELOW,EMAX,FINT,FRAT,HEAD,IDM,SZ,Z, C BLAM(MZ07,MZ07,5),EKLSQ(MZ01),ET(MZ07),ETERM(MZ02), C F1(2,MZ07,MZ07),F2(2,MZ07,MZ07),FI(MZ05,MZ05),FLAG(MZ09,3,4), C FX(2,MZ07,MZ07),P(MZ06,MZ03),Q(MZ06,MZ03),R(MZ06),S(MZ10,3), C SD1(MZ05,MZ05),SD2(MZ05,MZ05),SH(4,5),T(MZ10,3), C IFILE,IPRINT,ISTOT,ITC,ITCP1,ITMX,JFILE,KCASE,LAMMX,LCUT,LTOT, C MCH,MDG,MODE,N,N1,N2,NCHB,NCHCL,NCHF,NCHFMX,NCHOP,NLCTRN,NLOW, C NMX,NRAD,NTERM,NTHROW,NTOT,NTOTMX, C ISTERM(MZ02),ITERM(MZ07),LL(MZ07),LQN(MZ03),LTERM(MZ02), C MDEGEN(MZ01),MTCHCL(MZ01),NQN(MZ03) COMMON /POL1/ALPHA,RDICUT,NPOL,NDIEL,IGC COMMON /TXT/TEXT(8) DATA HED/4HP(I,/ C C READ ATOM DATA READ(IFILE,501) NRAD,JMAX,NLCTRN,Z,HEAD IF (JMAX.LE.MZ04) GO TO 202 WRITE(6,626) JMAX STOP 202 IF (NRAD.LE.MZ03) GO TO 203 WRITE(6,627) NRAD STOP 203 READ(IFILE,502) (I,RIN(I),ZEFF(I),J=1,JMAX) DO 200 IRAD=1,NRAD READ(IFILE,503) IG,NQN(IG),LQN(IG),NCARDS IF(IG.LE.NRAD) GO TO 199 WRITE(6,629) IG,NRAD STOP 199 IMAX=2*NCARDS READ(IFILE,502) (J,PIN(J,IG),QIN(J,IG),I=1,IMAX) IF (IMAX.EQ.JMAX) GO TO 200 IMP1=IMAX+1 DO 201 I=IMP1,JMAX PIN(I,IG)=0. 201 QIN(I,IG)=0. 200 CONTINUE A1=NLCTRN SZ=Z-A1 READ(IFILE,504) NTERM,NLOW,ELOW IF (NTERM.LE.MZ02) GO TO 204 WRITE(6,628) NTERM STOP 204 READ(IFILE,505) (JTERM,ISTERM(JTERM),LTERM(JTERM),ETERM(JTERM), 1I=1,NTERM) WRITE(6,600) HEAD,Z,NLCTRN,NRAD,NTERM,ELOW,NLOW WRITE(6,612) TEXT IF(MODE.EQ.2) WRITE(6,611) SZ C C CHECK TERM ENERGIES. C TWO TERMS ARE ASSUMED DEGENERATE IF DIFFERENCE IN C ENERGIES READ IS LESS THAN EROR. C COMPUTE EROR EROR=DEL2*0.2 IF(EROR.LE.1.E-5) EROR=1.E-5 IF(SZ.GT.1.5) EROR=EROR*SZ*SZ C CHECK ETERM(1) A1=ETERM(1) IF(A1.EQ.0) GO TO 211 IF(ABS(A1).LE.EROR) GO TO 210 WRITE(6,630) A1 STOP 210 WRITE(6,631) A1 ETERM(1)=0. C SCAN ENERGIES FOR DEGENERACIES A1=0. 211 KB=2 KA=1 212 IF(KB.GT.NTERM) GO TO 216 KC=0 A1=ETERM(KA) A3=A1-EROR DO 214 J=KB,NTERM A2=ETERM(J) IF(A2.EQ.A1) GO TO 214 IF(A2.GE.A3) GO TO 213 WRITE(6,632) J,A2,KA,A1 STOP 213 IF(ABS(A2-A1).GT.EROR) GO TO 215 WRITE(6,633) J,A2,KA,A1 ETERM(J)=A1 214 KC=KC+1 215 KB=KB+KC+1 KA=KB-1 GO TO 212 216 CONTINUE WRITE(6,607) C C WRITE OUT INFORMATION ON TERMS WRITE(6,601) (I,ISTERM(I),LTERM(I),ETERM(I),I=1,NTERM) DEL1P=DEL1 C C FIND PMAX, AND MODIFY DEL2 DEL2P=DEL2 PMAX=1.E10 DO 30 IG=1,NRAD A1=0. DO 35 IT=2,JMAX A2=ABS(PIN(IT,IG)) IF (A2.GT.A1) A1=A2 35 CONTINUE IF (A1.LT.PMAX) PMAX=A1 30 CONTINUE DEL2=PMAX*DEL2 C C FIND LARGEST VALUE OF ABS(PIN(JMAX,IG)) AND CHECK .LT. DEL2 A1=0. DO 40 IG=1,NRAD A2=ABS(PIN(JMAX,IG)) IF (A2.GT.A1) A1=A2 40 CONTINUE IF (DEL2.GE.A1) GO TO 41 WRITE(6,624)A1,DEL2,A1 DEL2=A1 C C TO FIND MAXJ, THE SMALLEST J SUCH THAT ABS(PIN(J,IG)).LT.DEL2 FOR C IG=1,NRAD 41 IT=JMAX+1 DO 1 I=1,JMAX IT=IT-1 DO 2 IG=1,NRAD IF (ABS(PIN(IT,IG)) .LT. DEL2) GO TO 2 MAXJ=IT+1 IF (I.EQ.1) GO TO 25 GO TO 3 2 CONTINUE 1 CONTINUE 25 MAXJ=MAXJ-1 C C TO FIND PHI(I), I=1,JMAX C PHI(I)=INTEGRAL FROM 0.0 TO RIN(I) OF SQRT(2.*ZEFF(R)/R+EMAX) 3 PHI(1)=0. A1=RIN(1) A3=ZEFF(1) DO 4 I=2,JMAX A2=RIN(I) A4=ZEFF(I) DEL=A2-A1 A=(A3*A2-A4*A1)/DEL B=(A4-A3)/DEL A=2.*A B=2.*B+EMAX X=0.5*(A1+A2) RINT=A+B*X RINT=RINT*RINT RINT=(3.*A/4.+B*X)/RINT A1=24.*X*X A1=DEL*DEL*A/A1 RINT=1.+A1*RINT A1=(A/X+B) A1=DEL*SQRT(A1) RINT=A1*RINT A1=A2 A3=A4 4 PHI(I)=PHI(I-1)+RINT C C SET UP R(1) AND R(2) AND COMPUTE VALUES OF P(1,IG) AND Q(1,IG) DEL=3.141593/FINT DEL1=DEL1/Z R(1)=0. R(2)=DEL1 DO 23 IG=1,NRAD P(1,IG)=PIN(1,IG) 23 Q(1,IG)=QIN(1,IG) JM1=1 JMXM1=JMAX-1 RJMXM1=RIN(JMXM1) C C BEGIN DO LOOP WHICH COMPUTES INTERVALS FOR IT=3,N AND FUNCTIONS C P AND Q FOR IT=2,N RM1=0. ICON=0 DO 5 IT=2,NMX IF (IT.EQ.2) GO TO 12 C C CALCULATES TABULATION POINTS R(IT) RM2=RM1 RM1=R(IT-1) DEL1=(RM1-RM2)*FRAT A1=RM1-RIN(JM1) A2=RIN(J)-RM1 PH=(A1*PHI(J)+A2*PHI(JM1))/(A1+A2) DO 8 IJ=JM1,JMXM1 A3=PHI(IJ)-PH IF (A3 .LT. DEL) IM1=IJ IF (IM1 .EQ. IJ) GO TO 8 I=IJ GO TO 9 8 CONTINUE R(IT)=RM1+DEL1 IF (R(IT).LT.RJMXM1) GO TO 12 N=IT-1 GO TO 10 9 A1=PHI(IM1)-PH A2=PHI(I)-PH A3=RIN(IM1) A4=RIN(I) A=(A4-A3)*DEL+A2*A3-A1*A4 A=A/(A2-A1) DEL2=A-RM1 IF (DEL1 .LE. DEL2) GO TO 11 R(IT)=A GO TO 7 11 R(IT)=RM1+DEL1 12 CONTINUE C R(IT) NOW CALCULATED C DO 6 IJ=JM1,JMXM1 IF (RIN(IJ).LE.R(IT)) IM1=IJ IF (IM1.EQ.IJ) GO TO 6 I=IJ GO TO 7 6 CONTINUE C C CALCULATES P(IT,IG) AND Q(IT,IG) FOR IG=1,NRAD 7 IMIN=1 IF (I.EQ.2) IMIN=2 DO 13 IJ=IMIN,4 13 D(IJ)=RIN(I+IJ-3) A1=R(IT) DO 16 IJ=IMIN,4 A2=1. DO 14 K=IMIN,4 IF (K.EQ.IJ) GO TO 14 A2=A2*(D(IJ)-D(K)) 14 CONTINUE A2=1./A2 DO 15 K=IMIN,4 IF (K.EQ.IJ) GO TO 15 A2=A2*(A1-D(K)) 15 CONTINUE FL(IJ)=A2 16 FL1(IJ)=A2 IJMIN=1 IF(I.EQ.2) IJMIN=3 IF(I.EQ.3) IJMIN=2 L0=I+IMIN-4 DO 17 IG=1,NRAD IF (IT.GT.ITC) GO TO 18 L1=-LQN(IG)-1 DO 19 IJ=IJMIN,4 19 FL1(IJ)=FL(IJ)*D(IJ)**L1 18 A1=0. A2=0. L1=L0 DO 20 K=IMIN,4 L1=L1+1 A1=A1+FL1(K)*PIN(L1,IG) 20 A2=A2+FL1(K)*QIN(L1,IG) P(IT,IG)=A1 17 Q(IT,IG)=A2 J=I JM1=IM1 IF (J.LE.MAXJ) GO TO 5 N=IT GO TO 10 5 CONTINUE N=NMX WRITE(6,621) ICON=1 C 10 IF(N.GT.ITC) GO TO 65 WRITE(6,625) N,ITC STOP 65 N1=N+1 N2=N+2 A1=NLCTRN SZ=Z-A1 IF (MODE.NE.2) GO TO 66 Z=SZ NLCTRN=0 66 DEL2=DEL2P DEL1=DEL1P C C EXTEND RANGE OF R RIT=R(N) A1=RIT-R(N-1) R(N+1)=RIT+A1 R(N+2)=RIT+A1+A1 C C PRINT INTERVALS AND ATOMIC RADIAL FUNCTIONS WRITE(6,610) WRITE(6,602) FINT,FRAT,EMAX,ITC,DEL1,DEL2,N JPRINT=IPRINT IF(NTHROW.EQ.1) JPRINT=-JPRINT WRITE(6,613) MODE,NPOL,JPRINT IF(ICON.EQ.1) GO TO 72 IF(IPRINT.EQ.1.OR.IPRINT.EQ.3) RETURN 72 IF(NTHROW.EQ.0) WRITE(6,603) IF(NTHROW.EQ.1) WRITE(6,1603) WRITE(6,604)(I,NQN(I),LQN(I),I=1,NRAD) I1=NRAD/9+1 I2=1 DO 80 K=1,I1 I3=I2+7 IF(I3.GT.NRAD) I3=NRAD WRITE(6,605)(HED,J,J=I2,I3) WRITE(6,607) DO 78 I=1,ITC 78 WRITE(6,606) I,R(I),(P(I,J),J=I2,I3) WRITE(6,607) DO 79 I=ITCP1,N 79 WRITE(6,606) I,R(I),(P(I,J),J=I2,I3) 80 I2=I2+8 WRITE(6,607) N1=N+1 DO 81 IT=N1,N2 81 WRITE(6,606) IT,R(IT) RETURN C 501 FORMAT (5X,I5,15X,I3,I4,F4.0,36X,A4) 502 FORMAT(5X,I4,2E14.7,I4,2E14.7) 503 FORMAT (5X,2I5,I3,19X,I6) 504 FORMAT (5X,I5,5X,I3,F20.8) 505 FORMAT (5X,2I5,I3,F20.8) 600 FORMAT(//,29H ELECTRON COLLISIONS WITH ,A4,/,1H ,32(1H*), 1 //36H SPECIFICATION OF TARGET SYSTEM - / 19H NUCLEAR 2 CHARGE =,F4.0/ 25H NUMBER OF ELECTRONS =,I3/ 41H NUMBER OF 3 ATOMIC RADIAL FUNCTIONS =,I3/ 21H NUMBER OF TERMS =,I3/ 428H ENERGY OF LOWEST TERM =,F20.8,11H (RYDBERGS),/, 5 29H ( LOWEST TERM IS NUMBER ,I3,2H )/) 601 FORMAT(46H TERM PARITY*(2S+1) L ENERGY (RYDBERGS) 1/,(I3,7X,I4,8X,I3,F12.6)) 602 FORMAT(6H FINT=,F4.2,7H, FRAT=,F4.2,7H, EMAX=,F7.2,6H, ITC=, C I2,7H, DEL1=,E10.3,7H, DEL2=,E10.3,/,14H THESE GIVE N=,I4) 603 FORMAT(1H1,23HTARGET RADIAL FUNCTIONS/1X,23(1H*)//1X,4HIGAM,4X, C 3HNQN,5X,3HLQN) 1603 FORMAT(/////,24H TARGET RADIAL FUNCTIONS/1X,23(1H*)//1X,4HIGAM,4X, C3HNQN,5X,3HLQN) 604 FORMAT(I4,6X,I2,6X,I2) 605 FORMAT(//3X,1HI,6X,4HR(I),5X,8(4X,A4,I2,2H) )) 606 FORMAT(I4,F12.6,3X,8F12.6) 607 FORMAT(/) 610 FORMAT(//35H PARAMETERS DETERMINING INTERVALS - ) 611 FORMAT(50H OPERATING IN MODE 2 - SOLUTION OF COULOMB PROBLEM , C 14H FOR CHARGE = ,F5.1) 612 FORMAT(/,24H TEXT READ ON FIRST CARD,10X,10(1H*),8A4,10(1H*)/) 613 FORMAT(/18H RUN WITH MODE =,I2,/,10X,8HNPOL =,I2,/,10X,8HIPRIN CT =,I2) C C FORMATS FOR MESSAGES 621 FORMAT ( 110H0NUMBER OF POINTS N NEEDED FOR P(IT,IG) TO BE .LT. 1DEL2 FOR ALL R .GT. R(N) EXCEEDS MAX ALLOWED BY DIMENSIONS. / 229H CHECK P(N,IG) FOR MAGNITUDE. ) 624 FORMAT ( 38H0MAX VALUE OF PIN(JMAX,IG),IG=1,NRAD = ,E10.4, 1 29H ,WHICH IS GREATER THAN DEL2= ,E10.4,21H . DEL2 PUT EQUAL TO , 2 E10.4 ) 625 FORMAT ( 28H0CALCULATED VALUE OF N =,I2,37H, IS .LE. INPU 1T VALUE OF ITC =,I2,24H RUN STOPPED IN ATOM ) 626 FORMAT ( 15H0VALUE OF JMAX=,I3,52H WHICH IS GREATER THAN MAXIMUM 1 ALLOWED BY DIMENSIONS ) 627 FORMAT ( 15H0VALUE OF NRAD=,I3,52H WHICH IS GREATER THAN MAXIMUM 1 ALLOWED BY DIMENSIONS ) 628 FORMAT ( 18H0NUMBER OF TERMS =,I3,52H WHICH IS GREATER THAN MAXI 1MUM ALLOWED BY DIMENSIONS ) 629 FORMAT(//,34H ***ERROR*** ORBITAL INDEX IG = ,I3,32H ON COLALG I 1NPUT EXCEEDS NRAD = ,I3,/,15X,30H.THIS IS NOT A DIMENSION FAULT,/, 215X,51HPROBABLE CAUSE IS THAT THE TARGET ORBITALS WERE NOT,/,56HNU 3MBERED CONSECUTIVELY DURING COLALG (OR EQUIVALENT) RUN ) 630 FORMAT(/,25H *** ERROR *** ETERM(1)=,E15.9,18H (SHOULD BE ZERO)) 631 FORMAT(/,23H *** ETERM(1) READ AS ,E15.9,19H. SET EQUAL TO 0.0) 632 FORMAT(/,22H *** ERROR *** ETERM(,I2,2H)=,E15.9,20H IS LESS THAN CETERM(,I2,2H)=,E15.9,/,16X,45HTERMS SHOULD BE IN ORDER OF INCREASI CNG ENERGY) 633 FORMAT(/,12H *** ETERM(,I2,9H) READ AS,E15.9,22H. SET EQUAL TO E CTERM(,I2,2H)=,E15.9) END C*********************************************************************** BLOCK DATA C INCLUDE 'PARAM' C C USED TO INITIALISE ARRAY SH C INTEGER HEAD COMMON /HOLD/DEL1,DEL2,DELNU,ELOW,EMAX,FINT,FRAT,HEAD,IDM,SZ,Z, C BLAM(MZ07,MZ07,5),EKLSQ(MZ01),ET(MZ07),ETERM(MZ02), C F1(2,MZ07,MZ07),F2(2,MZ07,MZ07),FI(MZ05,MZ05),FLAG(MZ09,3,4), C FX(2,MZ07,MZ07),P(MZ06,MZ03),Q(MZ06,MZ03),R(MZ06),S(MZ10,3), C SD1(MZ05,MZ05),SD2(MZ05,MZ05),SH(4,5),T(MZ10,3), C IFILE,IPRINT,ISTOT,ITC,ITCP1,ITMX,JFILE,KCASE,LAMMX,LCUT,LTOT, C MCH,MDG,MODE,N,N1,N2,NCHB,NCHCL,NCHF,NCHFMX,NCHOP,NLCTRN,NLOW, C NMX,NRAD,NTERM,NTHROW,NTOT,NTOTMX, C ISTERM(MZ02),ITERM(MZ07),LL(MZ07),LQN(MZ03),LTERM(MZ02), C MDEGEN(MZ01),MTCHCL(MZ01),NQN(MZ03) DATA SH /251.,232.,243.,224.,646.,992.,918.,1024.,-264., 1192.,648.,384.,106.,32.,378.,1024.,-19.,-8.,-27.,224./ END C*********************************************************************** SUBROUTINE BOOLE(BS) C INCLUDE 'PARAM' C C BOOLE"S RULE INTEGRATION OF N POINTS ON MESH H C INTEGRATION FROM ORIGIN TO EACH MESH POINT C BI IS INPUT ARRAY C BO IS OUTPUT ARRAY C BS=VALUE OF INTEGRAND AT ORIGIN C COMMON /CMN14/ CHN(10),BI(10),BO(10),P1(10),Z1(10) COMMON /CMN04/ H,N,JJ H2=(H+H)/45. B1=BI(1) B2=BI(2) B3=BI(3) B4=BI(4) BO(1)=H*(646.*B1-264.*B2+106.*B3-19.*B4+251.*BS)/720. BO(2)=H*(124.*B1+24.*B2+4.*B3-B4+29.*BS)/90. BO(3)=3.*H*(34.*B1+24.*B2+14.*B3-B4+9.*BS)/80. BO(4)=H2*(7.*(B4+BS)+32.*(B1+B3)+12.*B2) DO 90 J=5,N B5=BI(J) BO(J)=BO(J-4)+H2*(7.*(B1+B5)+32.*(B2+B4)+12.*B3) B1=B2 B2=B3 B3=B4 90 B4=B5 RETURN END C********************************************************************** SUBROUTINE BOUND C INCLUDE 'PARAM' C DESCRIBED IN CPC SECTION 8.3 C C FOR THE CASE OF ALL CHANNELS CLOSED, CALCULATES THE ENERGY EIGEN- C VALUE (E1) CLOSEST TO THE INPUT ENERGY (SQK=EKLSQ(KCASE)), AND THE C CORRESPONDING EIGENVECTORS (G). C PUT E1=SQK + FL . THE EQUATIONS TO BE SOLVED ARE (C+FL*C1)*G=0 . C THE ARRAY C HAS BEEN CALCULATED FOR ENERGY SQK. C THE ELEMENTS OF THE ARRAY C1 ARE COMPUTED WHEN REQUIRED (THE ARRAY C C1 IS NOT STORED). C THE EQUATIONS ARE SOLVED USING AN ITERATIVE NUMERICAL METHOD C DESCRIBED IN "MODERN COMPUTING METHODS", 1961, HM STATIONARY C OFFICE, PAGE 25. C C INTEGER HEAD COMMON C(MZ08,MZ08),G(MZ08),GP(MZ08),GPP(MZ08),CNORM(MZ08) COMMON /WORK3/C1(MZ07),C2(MZ07),FNU(MZ07),PHI(MZ05) COMMON /ORDER/IROW(MZ08),ICOL(MZ12) COMMON /OUT/IPUNCH,ICONV COMMON /HOLD/DEL1,DEL2,DELNU,ELOW,EMAX,FINT,FRAT,HEAD,IDM,SZ,Z, C BLAM(MZ07,MZ07,5),EKLSQ(MZ01),ET(MZ07),ETERM(MZ02), C F1(2,MZ07,MZ07),F2(2,MZ07,MZ07),FI(MZ05,MZ05),FLAG(MZ09,3,4), C FX(2,MZ07,MZ07),P(MZ06,MZ03),Q(MZ06,MZ03),R(MZ06),S(MZ10,3), C SD1(MZ05,MZ05),SD2(MZ05,MZ05),SH(4,5),T(MZ10,3), C IFILE,IPRINT,ISTOT,ITC,ITCP1,ITMX,JFILE,KCASE,LAMMX,LCUT,LTOT, C MCH,MDG,MODE,N,N1,N2,NCHB,NCHCL,NCHF,NCHFMX,NCHOP,NLCTRN,NLOW, C NMX,NRAD,NTERM,NTHROW,NTOT,NTOTMX, C ISTERM(MZ02),ITERM(MZ07),LL(MZ07),LQN(MZ03),LTERM(MZ02), C MDEGEN(MZ01),MTCHCL(MZ01),NQN(MZ03) COMMON /SCALE/CSCALE(MZ12),RSCALE(MZ08) C C INITIALIZATION DO 4 I=1,NTOT 4 GP(I)=0. IF(NTHROW.EQ.0) WRITE(6,630) IF(NTHROW.EQ.1) WRITE(6,1630) SQK=EKLSQ(KCASE) FL=0. NITT=0 IJKM=ITMX NTM1=NTOT-1 IF(MODE.EQ.2) NTM1=MCH*N2-1 NTM=NTOT IF(MODE.EQ.2) NTM=NTM1+1 ICONV=0 C C CONVERGENCE OF EFFECTIVE QUANTUM NUMBER IS TESTED IN CHANNEL MCH ETMCH=ET(MCH) C INPUT VALUE OF EFFECTIVE QUANTUM NUMBER IN CHANNEL MCH IS FN FN=-ETMCH IF(SZ.NE.0) FN=SZ/SQRT(FN) WRITE(6,631)SQK,MCH,FN WRITE(6,632)NITT,SQK,FN FNP=FN C C OBTAIN INITIAL ESTIMATE OF EIGENVECTOR. THIS ESTIMATE IS GP IN C PRIMED REPRESENTATION. PUT GP(NTOT)=1 AND SOLVE C*GP=0 FOR ROWS C 1 TO (NTOT-1) I=NTM1 A1=-C(I,NTM)/C(I,I) A2=1+A1*A1 GP(I)=A1 DO 10 K=3,NTM IP1=I I=I-1 A1=C(I,NTM) DO 5 J=IP1,NTM1 5 A1=A1+C(I,J)*GP(J) A1=-A1/C(I,I) A2=A2+A1*A1 10 GP(I)=A1 A2=1./SQRT(A2) C NORMALISE GP TO (SUM OVER I OF GP(I)**2)=1 DO 15 I=1,NTM1 15 GP(I)=GP(I)*A2 GP(NTM)=A2 NTM1=NTOT-1 IF(IPRINT.EQ.5) CALL PRINTB(FL) C C BEGIN DO LOOP ON ITERATIONS C --------------------------- DO 20 IJK=1,ITMX C C STORE PREVIOUS VECTOR GP IN GPP DO 21 I=1,NTOT 21 GPP(I)=GP(I) C C COMPUTE G=(GP IN UNPRIMED REPRESENTATION) DO 22 I=1,NTOT K=ICOL(I) 22 G(K)=GP(I) C C RETURN TO ORIGINAL NORMALISATION IN STARTING REGION, DO 35 I=1,NTOT 35 G(I)=G(I)*CSCALE(I) C C REPLACE G BY -C1*G C ------------------ K=0 DO 40 I=1,NCHF K=K+N2 C1(I)=G(K-1) 40 C2(I)=G(K) C C ......FOR FREE CHANNEL EQUATIONS ...... LIP=-N2 DO 50 I=1,NCHF LIP=LIP+N2 L=LL(I) FL=L B1=C2(I) A1=Z/(FL+1.) C DO 60 IT=1,ITC 60 PHI(IT)=1.-A1*R(IT) C C ......FOR IT=ITCP1,N2 DO 70 LI=ITCP1,N2 IT=N2+ITCP1-LI M=3 IF (IT.EQ.N1) M=2 IF (IT.EQ.N2) M=1 C A1=0. DO 80 J=1,M IS=IT-3+J IP=LIP+IS A2=1. IF (IS.LE.ITC) A2=R(IS)**(L+3) 80 A1=A1-T(IT,J)*A2*G(IP) 70 G(LIP+IT)=A1 C LI=LIP+ITCP1 A1=PHI(ITC-1)*R(ITC-1)**(L+1) A2=PHI(ITC)*R(ITC)**(L+1) G(LI)=G(LI)-(T(ITCP1,1)*A1+T(ITCP1,2)*A2)*B1 LI=LI+1 G(LI)=G(LI)-T(ITC+2,1)*A2*B1 C A1=0. A2=0. DO 90 J=1,NCHF A3=C1(J) A1=A1+FX(1,I,J)*A3 A2=A2+FX(2,I,J)*A3 90 CONTINUE C G(LIP+N1)=G(LIP+N1)-A1 G(LIP+N2)=G(LIP+N2)-A2 C C ......FOR IT=1,ITC G(LIP+1)=PHI(1)*B1 DO 100 IT=2,ITC RIT=R(IT) IS=LIP+IT A2=RIT**(L+1) A1=PHI(IT)*B1*A2 100 G(IS)=G(IS)*A2*RIT*RIT+A1 50 CONTINUE C C ......FOR LAM EQUATIONS LIP=N2*NCHF+NCHB IF (LIP.EQ.NTOT) GO TO 110 LI=LIP+1 DO 120 IT=LI,NTOT 120 G(IT)=0. 110 CONTINUE C C REPLACEMENT COMPLETED C --------------------- C C NORMALIZE G FOR ROW NORMALIZATION OF MATRIX DO 125 I=1,NTOT 125 G(I)=G(I)*RSCALE(I) C DO 126 I=1,NTOT 126 GP(I)=G(I) C +++++ DELETED SECTION OF PETES VERSION FOR ROW INTERCHANGES ++++ C C REPLACE GP BY (C**(-1))*GP C INVERT LOWER TRIANGULAR MATRIX K=0 DO 127 J=2,NTOT K=K+1 A1=GP(J) DO 128 M=1,K 128 A1=A1-C(J,M)*GP(M) 127 GP(J)=A1 C INVERT UPPER TRIANGULAR MATRIX K=NTOT GP(K)=GP(K)/C(K,K) DO 129 I=1,NTM1 K=K-1 A1=GP(K) M=NTOT+1 DO 131 J=1,I M=M-1 131 A1=A1-C(K,M)*GP(M) 129 GP(K)=A1/C(K,K) C C COMPUTE NEW ESTIMATE OF EIGENVALUE FL A1=0. A2=0. DO 400 I=1,NTOT A3=GP(I) A1=A1+GPP(I)*A3 400 A2=A2+A3*A3 A2=1./A2 FL=A1*A2 C C CHECK THAT NEW ENERGY GIVES ALL CHANNELS CLOSED DO 405 I=1,NCHF A1=ET(I)+FL IF(A1.LT.0) GO TO 405 WRITE(6,613) IJK,I,A1 RETURN 405 CONTINUE C C NEW ESTIMATE FOR ENERGY IS E1=SQK+FL C NEW EFFECTIVE QUANTUM NUMBER IN CHANNEL MCH IS FNP=-ETMCH-FL IF(SZ.NE.0) FNP=SZ/SQRT(FNP) WRITE(6,632)IJK,E1,FNP C C NORMALIZE GP A2=SQRT(A2) DO 410 I=1,NTOT 410 GP(I)=GP(I)*A2 C C TEST CONVERGENCE IF(ABS(FN-FNP).GT.DELNU) GO TO 415 C C EXIT FROM DO LOOP IF CONVERGENCE TEST SATISFIED NITT=IJK GO TO 420 C C CONTINUE IF CONVERGENCE TEST NOT SATISFIED 415 FN=FNP IF(IPRINT.EQ.5) CALL PRINTB(FL) C 20 CONTINUE C DO LOOP ON ITERATIONS COMPLETED C ------------------------------- C C FAILURE TO CONVERGE AFTER ITMX ITERATIONS WRITE(6,611) ITMX,DELNU RETURN C C CONVERGENCE CRITERION SATISFIED 420 ICONV=1 CALL PRINTB(FL) WRITE(6,610) WRITE(6,612)NITT,MCH,DELNU WRITE(6,633)E1 C C COMPUTE AND PRINT EFFECTIVE QUANTUM NUMBERS IF(JFILE.GT.0) FL=0. DO 240 I=1,NCHF 240 FNU(I)=SZ/SQRT(-ET(I)-FL) WRITE(6,603)(I,I=1,NCHF) WRITE(6,604)(FNU(I),I=1,NCHF) C 603 FORMAT(/12H CHANNEL, I= ,8X,I2,9(10X,I2)/ C (20X,I2,9(10X,I2))) 604 FORMAT(12H FNU(I) =,10F12.6,/(12X,10F12.6)) 610 FORMAT (//) 630 FORMAT(1H1,18H SUBROUTINE BOUND/19H ****************// C55H CALCULATES ENERGY E1 (RYDBERG UNITS) OF ADDED ELECTRON , C37H WHEN ION CORE IS IN THE GROUND STATE//) 1630 FORMAT(/////,19H SUBROUTINE BOUND/3X,16(1H*)//, C55H CALCULATES ENERGY E1 (RYDBERG UNITS) OF ADDED ELECTRON , C37H WHEN ION CORE IS IN THE GROUND STATE //) 631 FORMAT(28H INITIAL ESTIMATE OF E1=SQK=,E13.5/ C56H INITIAL ESTIMATE OF EFFECTIVE QUANTUM NUMBER IN CHANNEL,I3, C6H , FN=,F10.6/) 632 FORMAT(10H ITERATION,I3, 5H ,E1=,E13.5, 5H ,FN=,F10.6) 633 FORMAT( 4H E1=,E13.5//) C FORMATS FOR MESSAGES 611 FORMAT (1H ,6HAFTER ,I3,46H ITERATIONS , CONVERGENCE TO ACCURA 1CY OF ,E13.4,19H HAD NOT OCCURED ) 612 FORMAT(1H ,6HAFTER ,I2,48H ITERATIONS EFFECTIVE QUANTUM NUMBER IN 1CHANNEL ,I2,26H CONVERGED TO ACCURACY OF ,E13.4) 613 FORMAT (1H0,20HON ITERATION NUMBER ,I3,35H VALUE OF ENERGY IN CHAN 1NEL NUMBER ,I3,5H WAS ,E12.4,26H, I.E. NOT LESS THAN ZERO. / 286H THE PROGRAM THEREFORE BYPASSED THE EIGENVALUE PROBLEM AND WENT 3 ON TO NEXT ENERGY/CASE ) RETURN END C ********************************************************************** SUBROUTINE COULGM(AG,E2,PI,DG) C INCLUDE 'PARAM' C SEE CPC APPENDIX C SECTION (X) C C AG = L(J) ANGULAR MOMENTUM OF CHANNEL J C E2 = SCALED ENERGY (K**2)/(Z**2) C PI = 3.1415926 C C USED ONLY FOR Z .GT. 0 AND E2 .GT. 0 C C CALCULATES DG AS DEFINED BY EQUATIONS (C24) AND (C30) C L=IFIX(AG) B=1./SQRT(E2) C C INITIALISE P = L+1 IP=L+1 A=FLOAT(IP) C C INITIALISE DG TO MULTIPLE OF PI DG=-(AG+0.25)*PI C C SET UP CONDITION (C31) L1=IP A1=1./3. A1=25.*(B**A1)-B*B C C CHECK FOR LARGE ENOUGH VALUE OF P (C31) IF(A1.LE.0.) GO TO 2 1 A2=A*A IF(A2.GT.A1) GO TO 2 C INCREMENT P BY 1 IP=IP+1 A=A+1. GO TO 1 C C ADD ON TERMS TO OBTAIN DG AS DEFINED BY EQUATION (C30) 2 A1=A/B A2=A1*A1 A3=A2+1. A4=B*A3 DG=DG-B*0.5*LOG(A3)+(A-.5)*ATAN(A1)+1./(12.*A4) C +(1.-3.*A2)/(360.*A4*A4*A4) C C ADD ON TERMS OF EQUATION (C32) TO DG C (NOTING THAT CORRECT MULTIPLE OF PI HAS BEEN ACCOUNTED FOR ABOVE) IF(IP.EQ.L1) RETURN L2=IP-1 DO 3 I=L1,L2 A=FLOAT(I)/B 3 DG=DG-ATAN(A) RETURN C END C*********************************************************************** SUBROUTINE DBLU(W,X) C INCLUDE 'PARAM' C C COMPUTES POTENTIAL W AT R=X C C INTEGER HEAD COMMON NCHD1,IJNZ(MZ25),KNZ(5),BLIN(MZ25) C , RFGST(MZ22,MZ07,MZ07) COMMON /HOLD/DEL1,DEL2,DELNU,ELOW,EMAX,FINT,FRAT,HEAD,IDM,SZ,Z, C BLAM(MZ07,MZ07,5),EKLSQ(MZ01),ET(MZ07),ETERM(MZ02), C F1(2,MZ07,MZ07),F2(2,MZ07,MZ07),FI(MZ05,MZ05),FLAG(MZ09,3,4), C FX(2,MZ07,MZ07),P(MZ06,MZ03),Q(MZ06,MZ03),R(MZ06),S(MZ10,3), C SD1(MZ05,MZ05),SD2(MZ05,MZ05),SH(4,5),T(MZ10,3), C IFILE,IPRINT,ISTOT,ITC,ITCP1,ITMX,JFILE,KCASE,LAMMX,LCUT,LTOT, C MCH,MDG,MODE,N,N1,N2,NCHB,NCHCL,NCHF,NCHFMX,NCHOP,NLCTRN,NLOW, C NMX,NRAD,NTERM,NTHROW,NTOT,NTOTMX, C ISTERM(MZ02),ITERM(MZ07),LL(MZ07),LQN(MZ03),LTERM(MZ02), C MDEGEN(MZ01),MTCHCL(MZ01),NQN(MZ03) DIMENSION W(MZ24) C A1=1./X A2=2.*SZ*A1 C C INITIALIZE W ON AND BELOW DIAGONAL L=1 IF(NCHF.EQ.1) GO TO 25 DO 20 J=2,NCHF W(L)=ET(J-1)+A2 DO 10 I=J,NCHF L=L+1 10 W(L)=0. 20 L=L+NCHD1+J 25 W(L)=ET(NCHF)+A2 C C COMPLETE CALCULATION OF W ON AND BELOW DIAGONAL A2=A1 K2=0 DO 40 K=1,LAMMX A2=A2*A1 K1=K2+1 K2=KNZ(K) IF(K2.LT.K1) GO TO 40 DO 30 L=K1,K2 IJ=IJNZ(L) 30 W(IJ)=W(IJ)+A2*BLIN(L) 40 CONTINUE C C W ABOVE DIAGONAL IF(NCHF.EQ.1) RETURN L=1 DO 50 I=2,NCHF IJ=L DO 45 J=I,NCHF L=L+1 IJ=IJ+MZ07 45 W(IJ)=W(L) 50 L=L+NCHD1+I C RETURN END C*********************************************************************** SUBROUTINE DCHAIN( CHN,I,NP,NV,NAB,ND,C1,N1,C2,N2) C INCLUDE 'PARAM' C C CHAIN RULE DERIVATIVE OF GENERAL EXPRESSION C I=POSITION IN VECTOR ARRAYS C NP=OVERALL DERIVATIVE NP-1, MUST BE G.E.1 C NV=DERIVATE NV-1 OF POTENTIAL MATRIX C NV=0 IMPLIES SUBSTITUTE DELTA(I,K) FOR POTENTIAL MATRIX DERIVATIVE C NAB=1 FOR A VECTOR, =2 FOR B VECTOR C ND=DERIVATIVE ND-1 OF AB VECTOR, MUST BE G.E.1 C C1=SCALAR FUNCTION C1(J) C C2=SCALAR FUNCTION C2(J) C ARRAY CHN OUTPUTS RESULT AT EACH MESH POINT C DIMENSION C1(12,10),C2(12,10),CHN(10) C C COMMON NCHD1,IJNZ(MZ25),KNZ(5),BLIN(MZ25),NZD,NZDM, C AB1(MZ07,13,10),AB2(MZ07,13,10),ALP(13,13),CC(MZ07,MZ07,5), C E2(MZ07),PP(12,10),T1(MZ07),V(MZ11,11,10),ZET(12,10),ZTIN(12,10) C ,GAM(2,MZ07,15) COMMON /CMN02/ NCNL,NOP,LAMAX COMMON /CMN04/ HY,NHY,JJ DIMENSION NJ(MZ07) DIMENSION BLP(13),D1(12),D2(12) NU=NP K=NU-1 IF(K.EQ.0) K=1 IF(NV.NE.0) GO TO 101 IF(N1.NE.0) GO TO 104 IF(N2.NE.0) GO TO 103 R=ALP(1,K) IND=NU+ND-1 IF(NAB.GT.1) GO TO 368 DO 110 IPT=1,NHY 110 CHN(IPT)=AB1(I,IND,IPT)*R RETURN 368 DO 112 IPT=1,NHY 112 CHN(IPT)=AB2(I,IND,IPT)*R RETURN 101 DO 10 KVD=I,NCNL 10 NJ(KVD)=I+(KVD*KVD-KVD)/2 IF(I.EQ.1) GO TO 102 I2=I-1 I1=(I*I2)/2 DO 11 KVD=1,I2 11 NJ(KVD)=KVD+I1 C 102 CONTINUE DO 19 KP=1,NU 19 BLP(KP)=ALP(KP,K) ND1=ND-1 NUP1=NU+1 IF(NAB.GT.1) GO TO 187 DO 400 IPT=1,NHY DUM=0. IN1=0 NU1=NUP1 DO 390 KP=1,NU INV=NU1 NU1=NU1-1 K1=NU1-1 IF(K1.EQ.0) K1=1 CHA1=0. IND=ND1 DO 190 KQ1=1,NU1 INV=INV-1 IND=IND+1 S=0. DO 180 KD=1,NCNL IJ=NJ(KD) 180 S=S+V(IJ,INV,IPT)*AB1(KD,IND,IPT) 190 CHA1=CHA1+S*ALP(KQ1,K1) IF(N1.NE.0) GO TO 260 DUM=BLP(1)*CHA1 GO TO 409 260 IN1=IN1+1 CHA2=C1(IN1,IPT) DUM=DUM+BLP(KP)*CHA1*CHA2 390 CONTINUE 409 IF(DUM.EQ.0.) GO TO 600 400 CHN(IPT)=DUM RETURN C 187 DO 300 IPT=1,NHY DUM=0. IN1=0 NU1=NUP1 DO 290 KP=1,NU INV=NU1 NU1=NU1-1 K1=NU1-1 IF(K1.EQ.0) K1=1 CHA1=0. IND=ND1 DO 191 KQ1=1,NU1 INV=INV-1 IND=IND+1 S=0. DO 181 KD=1,NCNL IJ=NJ(KD) 181 S=S+V(IJ,INV,IPT)*AB2(KD,IND,IPT) 191 CHA1=CHA1+S*ALP(KQ1,K1) IF(N1.NE.0) GO TO 360 DUM=BLP(1)*CHA1 GO TO 509 360 IN1=IN1+1 CHA2=C1(IN1,IPT) DUM=DUM+BLP(KP)*CHA1*CHA2 290 CONTINUE 509 IF(DUM.EQ.0.) GO TO 600 300 CHN(IPT)=DUM RETURN 600 DO 610 IPT=1,NHY 610 CHN(IPT)=0. RETURN 103 DO 29 KP=1,NU 29 BLP(KP)=ALP(KP,K) NUD=NU+ND IF(NAB.EQ.2) GO TO 131 DO 500 IPT=1,NHY DUM=0. IN2=0 IND=NUD DO 490 KP=1,NU IN2=IN2+1 IND=IND-1 490 DUM=DUM+AB1(I,IND,IPT)*C2(IN2,IPT)*BLP(KP) IF(DUM.EQ.0.) GO TO 600 500 CHN(IPT)=DUM RETURN 131 DO 700 IPT=1,NHY DUM=0. IND=NUD IN2=0 DO 590 KP=1,NU IN2=IN2+1 IND=IND-1 590 DUM=DUM+AB2(I,IND,IPT)*C2(IN2,IPT)*BLP(KP) IF(DUM.EQ.0.) GO TO 600 700 CHN(IPT)=DUM RETURN 104 DO 39 KP=1,NU 39 BLP(KP)=ALP(KP,K) NUD=NU+ND N21=N2-1 IF(NAB.EQ.2) GO TO 850 DO 800 IPT=1,NHY DO 781 I1=1,12 D1(I1)=C1(I1,IPT) 781 D2(I1)=C2(I1,IPT) DUM=0. IND=NUD DO 790 KP=1,NU NU2=KP K2=NU2-1 IF(K2.EQ.0) K2=1 IN1=NU2+N1 IN2=N21 CHA2=0. DO 780 KQ2=1,NU2 IN1=IN1-1 IN2=IN2+1 780 CHA2=CHA2+ALP(KQ2,K2)*D1(IN1)*D2(IN2) IND=IND-1 DUM=DUM+BLP(KP)*CHA2*AB1(I,IND,IPT) 790 CONTINUE IF(DUM.EQ.0.0) GO TO 600 800 CHN(IPT)=DUM RETURN 850 DO 900 IPT=1,NHY DO 881 I1=1,12 D1(I1)=C1(I1,IPT) 881 D2(I1)=C2(I1,IPT) DUM=0. IND=NUD DO 890 KP=1,NU NU2=KP K2=NU2-1 IF(K2.EQ.0) K2=1 CHA2=0. IN1=NU2+N1 IN2=N21 DO 880 KQ2=1,NU2 IN1=IN1-1 IN2=IN2+1 880 CHA2=CHA2+ALP(KQ2,K2)*D1(IN1)*D2(IN2) IND=IND-1 DUM=DUM+BLP(KP)*CHA2*AB2(I,IND,IPT) 890 CONTINUE IF(DUM.EQ.0.0) GO TO 600 900 CHN(IPT)=DUM RETURN END C ********************************************************************* SUBROUTINE DIPOT C INCLUDE 'PARAM' C C REQUIRED FOR DI-ELECTRONIC POLARISATION C (SEE CPC SECTION 13.4, NPOL = 2) C COMPUTES DI(IT) SUCH THAT (P(IG)/SF/P(IGP))=SUM(IT=1,N OF C DI(IT)*P(IT,IGP)), WHERE L=LQN(IGP) C INTEGER HEAD COMMON /POL1/ALPHA,RDICUT,NPOL,NDIEL,IGC COMMON /POL2/DI(MZ06),SF(MZ06,4),SFP(MZ05) COMMON /OVSUB2/ IG1,IGP,IGPP,IG,L,IPQ,MOAD,LAM COMMON /HOLD/DEL1,DEL2,DELNU,ELOW,EMAX,FINT,FRAT,HEAD,IDM,SZ,Z, C BLAM(MZ07,MZ07,5),EKLSQ(MZ01),ET(MZ07),ETERM(MZ02), C F1(2,MZ07,MZ07),F2(2,MZ07,MZ07),FI(MZ05,MZ05),FLAG(MZ09,3,4), C FX(2,MZ07,MZ07),P(MZ06,MZ03),Q(MZ06,MZ03),R(MZ06),S(MZ10,3), C SD1(MZ05,MZ05),SD2(MZ05,MZ05),SH(4,5),T(MZ10,3), C IFILE,IPRINT,ISTOT,ITC,ITCP1,ITMX,JFILE,KCASE,LAMMX,LCUT,LTOT, C MCH,MDG,MODE,N,N1,N2,NCHB,NCHCL,NCHF,NCHFMX,NCHOP,NLCTRN,NLOW, C NMX,NRAD,NTERM,NTHROW,NTOT,NTOTMX, C ISTERM(MZ02),ITERM(MZ07),LL(MZ07),LQN(MZ03),LTERM(MZ02), C MDEGEN(MZ01),MTCHCL(MZ01),NQN(MZ03) COMMON /WORK1/PP(MZ06),PP1(MZ06),PP2(MZ06),RHO(5) C LG=LQN(IG)+1 DO 10 IT=1,N DI(IT)=0. 10 PP(IT)=P(IT,IG) C ITMC=0 ITM1=ITC-1 A5=R(ITM1)**LG PP(ITM1)=PP(ITM1)*A5 RIT=R(ITC) A6=RIT**LG PP(ITC)=PP(ITC)*A6 JMAX=4 A3=PP(ITC)*SF(ITC,1) IT3=ITC-3 DO 30 IT=ITCP1,N IF(IT.EQ.N) JMAX=3 ITM1=ITM1+1 ITMC=ITMC+1 IT3=IT3+1 RITM1=RIT RIT=R(IT) A4=A3 A3=PP(IT)*SF(IT,1) DR=0.25*(RIT-RITM1) DI(ITM1)=DI(ITM1)+SH(4,1)*DR*A4 DI(IT)=DI(IT)+SH(4,5)*DR*A3 DO 30 MU=2,4 MU1=MU-1 A1=0. DO 20 J=1,JMAX 20 A1=A1+FLAG(ITMC,MU1,J)*PP(IT3+J) A1=A1*SH(4,MU)*DR*SF(ITM1,MU) IS=IT3 DO 30 J=1,JMAX IS=IS+1 30 DI(IS)=DI(IS)+A1*FLAG(ITMC,MU1,J) DI(ITC-1)=DI(ITC-1)*A5 DI(ITC)=DI(ITC)*A6 C MU=LG+L DO 50 I=1,ITC A1=0. DO 40 IT=1,ITC 40 A1=A1+FI(I,IT)*P(IT,IG) 50 PP(I)=A1 C A1=R(ITC) A2=A1**MU DO 70 I=1,ITC A2=A2*A1 A3=A2*PP(I) MU=MU+1 I1=MU DO 70 IP=1,ITC A3=A3*A1 A4=A3*SFP(IP) I1=I1+1 A5=I1 DO 70 IPP=1,ITC A4=A4*A1 A5=A5+1. A6=A4/A5 DO 60 IT=1,ITC 60 DI(IT)=DI(IT)+A6*FI(IPP,IT) 70 CONTINUE RETURN END C*********************************************************************** SUBROUTINE ENERGY C INCLUDE 'PARAM' C DESCRIBED IN CPC SECTION 8.1 C C INSERTS ENERGY-DEPENDENT TERMS INTO MATRICES C COMMON /WORK3/C1(MZ07),C2(MZ07),FNU(MZ07),PHI(MZ05) C C INTEGER HEAD COMMON C(MZ08,MZ08),G(MZ08,MZ07) COMMON /HOLD/DEL1,DEL2,DELNU,ELOW,EMAX,FINT,FRAT,HEAD,IDM,SZ,Z, C BLAM(MZ07,MZ07,5),EKLSQ(MZ01),ET(MZ07),ETERM(MZ02), C F1(2,MZ07,MZ07),F2(2,MZ07,MZ07),FI(MZ05,MZ05),FLAG(MZ09,3,4), C FX(2,MZ07,MZ07),P(MZ06,MZ03),Q(MZ06,MZ03),R(MZ06),S(MZ10,3), C SD1(MZ05,MZ05),SD2(MZ05,MZ05),SH(4,5),T(MZ10,3), C IFILE,IPRINT,ISTOT,ITC,ITCP1,ITMX,JFILE,KCASE,LAMMX,LCUT,LTOT, C MCH,MDG,MODE,N,N1,N2,NCHB,NCHCL,NCHF,NCHFMX,NCHOP,NLCTRN,NLOW, C NMX,NRAD,NTERM,NTHROW,NTOT,NTOTMX, C ISTERM(MZ02),ITERM(MZ07),LL(MZ07),LQN(MZ03),LTERM(MZ02), C MDEGEN(MZ01),MTCHCL(MZ01),NQN(MZ03) C SQK=EKLSQ(KCASE) C C CALCULATE ARRAY FX USED FOR LINEARISING C-MATRIX IN SUBROUTINE BOUND IF (NCHOP.GT.0) GO TO 10 LI=-1 T2=T(N2,2) T3=T(N1,3) T4=T(N2,3) DO 200 I=1,NCHF LI=LI+N2 LIP=LI+1 E1=ET(I) E2=E1+DEL1 DO 200 J=1,NCHF A1=E2*F1(2,I,J)-E1*F1(1,I,J) A2=T2*A1 A1=T3*A1 A3=E2*F2(2,I,J)-E1*F2(1,I,J) A2=A2+T4*A3 K=0 DO 220 IP=1,NCHF K=K+N2 A3=F1(2,IP,J)-F1(1,IP,J) A4=F2(2,IP,J)-F2(1,IP,J) A1=A1+C(LI,K-1)*A3 220 A2=A2+C(LIP,K-1)*A3+C(LIP,K)*A4 FX(1,I,J)=A1/DEL1 200 FX(2,I,J)=A2/DEL1 C C INCLUDE ENERGY TERMS IN BOUND CHANNEL PART OF MATRIX 10 IF (NCHB.EQ.0) GO TO 20 LI=NCHF*N2 DO 30 I=1,NCHB LI=LI+1 30 C(LI,LI)=C(LI,LI)-SQK C C ADD TERMS INVOLVING K(I)**2,I=1,NCHF 20 DO 40 I=1,NCHF L=LL(I) FL=L E1=ET(I) A1=Z/(FL+1.) DO 50 IT=1,ITC 50 PHI(IT)=1.-A1*R(IT) LI=(I-1)*N2+1 G(LI,I)=G(LI,I)+E1*PHI(1) DO 60 IT=2,ITC RIT=R(IT) LI=LI+1 G(LI,I)=G(LI,I)+E1*PHI(IT)*RIT**(L+1) 60 C(LI,LI)=C(LI,LI)-E1*RIT**(L+3) A1=PHI(ITC-1)*R(ITC-1)**(L+1) A2=PHI(ITC)*R(ITC)**(L+1) A1=T(ITCP1,1)*A1+T(ITCP1,2)*A2 A2=T(ITC+2,1)*A2 IP=LI+1 G(IP,I)=G(IP,I)-A1*E1 IP=IP+1 G(IP,I)=G(IP,I)-A2*E1 DO 70 IT=ITCP1,N2 LI=LI+1 DO 80 J=1,3 IS=IT-3+J LIP=LI-3+J A1=T(IT,J)*E1 IF (IS.LE.ITC) A1=A1*R(IS)**(L+3) 80 C(LI,LIP)=C(LI,LIP)+A1 70 CONTINUE 40 CONTINUE C C INCORPORATE ASYMPTOTIC SOLUTIONS INTO MATRIX L1=-2 DO 90 I=1,NCHF L1=L1+N2 LI=L1 DO 90 M=1,2 LI=LI+1 LIP=0 DO 100 IP=1,NCHF LIP=LIP+N2 C1(IP)=C(LI,LIP-1) 100 C2(IP)=C(LI,LIP) L2=-2 DO 90 IP=1,NCHF L2=L2+N2 LIP=L2 DO 90 K=1,2 LIP=LIP+1 A1=0. DO 110 J=1,NCHF 110 A1=A1+C1(J)*F1(K,J,IP)+C2(J)*F2(K,J,IP) 90 C(LI,LIP)=A1 IF(NCHCL.EQ.0) RETURN DO 130 I=1,NCHCL NCL=NCHOP+I LI=NCL*N2 DO 130 IT=1,NTOT 130 C(IT,LI)=-G(IT,NCL) C RETURN END C*********************************************************************** SUBROUTINE EQNS C INCLUDE 'PARAM' C DESCRIBED IN CPC SECTION 6.1 C INTEGER HEAD COMMON C(MZ08,MZ08),G(MZ08,MZ07) COMMON /HOLD/DEL1,DEL2,DELNU,ELOW,EMAX,FINT,FRAT,HEAD,IDM,SZ,Z, C BLAM(MZ07,MZ07,5),EKLSQ(MZ01),ET(MZ07),ETERM(MZ02), C F1(2,MZ07,MZ07),F2(2,MZ07,MZ07),FI(MZ05,MZ05),FLAG(MZ09,3,4), C FX(2,MZ07,MZ07),P(MZ06,MZ03),Q(MZ06,MZ03),R(MZ06),S(MZ10,3), C SD1(MZ05,MZ05),SD2(MZ05,MZ05),SH(4,5),T(MZ10,3), C IFILE,IPRINT,ISTOT,ITC,ITCP1,ITMX,JFILE,KCASE,LAMMX,LCUT,LTOT, C MCH,MDG,MODE,N,N1,N2,NCHB,NCHCL,NCHF,NCHFMX,NCHOP,NLCTRN,NLOW, C NMX,NRAD,NTERM,NTHROW,NTOT,NTOTMX, C ISTERM(MZ02),ITERM(MZ07),LL(MZ07),LQN(MZ03),LTERM(MZ02), C MDEGEN(MZ01),MTCHCL(MZ01),NQN(MZ03) COMMON /OPS/Y(MZ06,MZ06),YPOT(MZ06) COMMON /OUT/IPUNCH,ICONV COMMON /OVSUB2/ IG,IGP,IGPP,IG1,L,IPQ,MOAD,LAM DIMENSION RL(MZ06),FLAP(MZ06),JB(MZ03),JORTHO(MZ03),PHI(MZ05) DIMENSION RITN1(5),RITN2(5),MDIAG(2) EQUIVALENCE (YPOT(1),RL(1),FLAP(1)),(JB(1),JORTHO(1),PHI(1)) EQUIVALENCE (FF,GG,AA,BB) COMMON /POL1/ ALPHA,RDICUT,NPOL,NDIEL,IGC C C ----INSERT FOR DI-ELECTRONIC POLARIZATION COMMON /POL2/DI(MZ06),SF(MZ06,4),SFP(MZ05) C FUNSF(SR) IS DEFINED IN FUNCTION SUBROUTINE FUNSF (NOTE THAT C THIS IS THE SQUARE ROOT OF THE POLARISATION POTENTIAL WHEN FUNSF C IS CALLED WITH NDIEL = 1) C COMPUTE SF(IT,J) = FUNSF(R(IT)+0.25*(J-1)*(R(IT+1)-R(IT))), C FUNSF = SUM I =1,2 ... ITC OF SFP(I)*R**I C SUCH THAT FUNSF(R)=(SUM I=1,ITC OF SFP(I)*R**I) IF(NDIEL.EQ.0) GO TO 9065 SF(1,1)=SQRT(ALPHA)/(RDICUT**3) DO 9010 I=2,ITC A1=R(I) 9010 SF(I,1)=FUNSF(A1)/A1 DO 9030 I=1,ITC A1=0. DO 9020 IT=1,ITC 9020 A1=A1+FI(I,IT)*SF(IT,1) 9030 SFP(I)=A1 DO 9040 I=2,ITC 9040 SF(I,1)=SF(I,1)*R(I) A1=R(ITC) I=ITC-1 DO 9060 IT=ITCP1,N I=I+1 A2=0.25*(R(IT)-A1) DO 9050 J=1,4 SF(I,J)=FUNSF(A1) 9050 A1=A1+A2 9060 CONTINUE SF(N,1)=FUNSF(A1) 9065 CONTINUE C ----END INSERT C C INITIALISE C AND BLAM TO ZERO - NRB DO LIP=1,NTOTMX DO LI=1,NTOTMX C(LI,LIP)=0. ENDDO ENDDO DO JLAM=1,5 DO IP=1,NCHFMX DO I=1,NCHFMX BLAM(I,IP,JLAM)=0. ENDDO ENDDO ENDDO C C READ AND PRINT GENERAL INFORMATION ON SYSTEM INCLUDING (2.*S+1.), C TOTAL L, NO OF FREE CHANNELS AND TOTAL NO OF CHANNELS READ(IFILE,501) ISTOT,LTOT,NCHF,NCHT,LCUT,ITMCUT,MDIAG IF (NCHF.LE.NCHFMX) GO TO 8 WRITE(6,611) NCHF GO TO 400 8 NCHB=NCHT-NCHF ICASE=KCASE IP=IABS(ITMCUT) J=IABS(ISTOT) LAMMX=1 IF(NTHROW.EQ.0) WRITE(6,600) IF(NTHROW.EQ.1) WRITE(6,1600) IF(ISTOT.GT.0) WRITE(6,606)J,LTOT IF(ISTOT.LT.0) WRITE(6,607)J,LTOT IF (LCUT.GT.0) WRITE(6,608) LCUT IF(IP.GT.0) WRITE(6,605) ITMCUT,IP IF(ITMCUT.LT.0) WRITE(6,612) IF(JFILE.GT.0) WRITE(6,615) MDIAG WRITE(6,601) NCHF,NCHB NTOT=NCHF*N2 IF (MODE.NE.2) NTOT=NTOT+NCHB IF(NTOT.LE.NTOTMX) GO TO 14 WRITE(6,610) NTOT GO TO 400 C C READ AND PRINT FREE CHANNEL DATA 14 DO 10 ICHF=1,NCHF READ(IFILE,502) I,ITERM(I),LL(I),(JORTHO(K),K=1,NRAD) WRITE(6,602) I,ITERM(I),LL(I),(JORTHO(K),K=1,NRAD) IP=N2*(I-1) C C INSERT ORTHOGONALITY CONDITIONS IF (MODE.EQ.2) GO TO 10 DO 9 K=1,NRAD IG=JORTHO(K) IF(IG.EQ.0) GO TO 9 NTOT=NTOT+1 IF(NTOT.LE.NTOTMX) GO TO 13 WRITE(6,610) NTOT GO TO 400 13 IF(MODE.EQ.1) GO TO 9 ILQN=LQN(IG)+1 IT=IP+1 C(IT,NTOT)=P(1,IG) DO 15 LLL=2,ITC IT=IT+1 15 C(IT,NTOT)=P(LLL,IG)*(R(LLL)**ILQN) DO 20 LLL=ITCP1,N IT=IT+1 20 C(IT,NTOT)=P(LLL,IG) IPQ=1 C CALL LAP IT=IP DO 30 LLL=1,N IT=IT+1 30 C(NTOT,IT)=FLAP(LLL) 9 CONTINUE 10 CONTINUE C C READ AND PRINT BOUND CHANNEL DATA IF(NCHB.EQ.0) GO TO 41 IF(MODE.NE.2) WRITE(6,603) (LIG,LIG=1,NRAD) DO 40 ICHB=1,NCHB READ(IFILE,505) J,(JB(K),K=1,NRAD) J=J-NCHF IF(MODE.NE.2) WRITE(6,604)J,(JB(K),K=1,NRAD) 40 CONTINUE IF (MODE.EQ.2) NCHB=0 C C PUNCH ON JFILE 41 IF(JFILE.EQ.0) GO TO 42 WRITE(JFILE,702) ISTOT,LTOT,HEAD,MDIAG WRITE(JFILE,701)ICASE,ISTOT,LTOT,NCHF,LCUT,NCHB WRITE(JFILE,701)(LL(I),I=1,NCHF) WRITE(JFILE,701)(ITERM(I),I=1,NCHF) C C READ ALL ALGEBRAIC COEFFICIENTS REQUIRED FOR SETTING UP EQUATIONS 42 MDE=(MODE+1)/2 IF(MDE.EQ.1) GO TO 51 50 READ(IFILE,504) KEY,J,I,IG,IP,IGP,LAM,FF GO TO (60,70,80,90,95,100),KEY C FOR MODE = 1 AND MODE = 2, SKIP PROCESSING OF CARDS C WITH KEY = 1 TO 6 51 READ(IFILE,506) KEY,LAM IF(KEY.EQ.6) GO TO 100 GO TO 51 C C COMPUTE FIRST FORM OF C AS DEFINED BY CPC EQNS. (6.14) TO (6.17) C C F COEFFICIENT 60 MOAD=1 CALL POT C ----INSERT FOR DI-ELECTRONIC POLARIZATION IF(NDIEL.EQ.0) GO TO 9090 IF(LAM.NE.1) GO TO 9090 IF(IG.LE.IGC.OR.IGP.LE.IGC) GO TO 9090 IG1=IG L=LQN(IGP) CALL DIPOT A1=0. DO 9070 IT=1,N 9070 A1=A1+DI(IT)*P(IT,IGP) DO 9080 IT=1,N 9080 YPOT(IT)=YPOT(IT)-A1*SF(IT,1) 9090 CONTINUE C ----END INSERT DO 110 IT=1,N 110 YPOT(IT)=2.*FF*YPOT(IT) LI=(I-1)*N2+1 LIP=(IP-1)*N2+1 I1=LAM+LL(IP)-LL(I) IF (I1.EQ.0) C(LI,LIP)=C(LI,LIP)+YPOT(1) IF (I.EQ.IP) GO TO 115 I2=LAM+LL(I)-LL(IP) IF (I2.EQ.0) C(LIP,LI)=C(LIP,LI)+YPOT(1) 115 CONTINUE DO 120 IT=2,ITC LI=LI+1 LIP=LIP+1 A2=R(IT)**(LL(I)+1) IF (I.EQ.IP) GO TO 120 A1=R(IT)**(LL(IP)+1) C(LI,LIP)=C(LI,LIP)+YPOT(IT)*A1 120 C(LIP,LI)=C(LIP,LI)+YPOT(IT)*A2 DO 130 IT=ITCP1,N LI=LI+1 LIP=LIP+1 IF (I.EQ.IP) GO TO 130 C(LI,LIP)=C(LI,LIP)+YPOT(IT) 130 C(LIP,LI)=C(LIP,LI)+YPOT(IT) IF(LAM.LT.1.OR.LAM.GT.5) GO TO 50 IF (LAM.GT.LAMMX) LAMMX=LAM A1=YPOT(N)*R(N)**(LAM+1) BLAM(I,IP,LAM)=BLAM(I,IP,LAM)+A1 IF (I.EQ.IP) GO TO 50 BLAM(IP,I,LAM)=BLAM(IP,I,LAM)+A1 GO TO 50 C C G COEFFICIENT 70 I1=I I2=IP IG1=IG IG2=IGP DO 140 IJK=1,2 L=LL(I2) MOAD=2 CALL POT C ----INSERT FOR DI-ELECTRONIC POLARIZATION IF(NDIEL.EQ.0) GO TO 9110 IF(LAM.NE.1) GO TO 9110 IF(IG1.LE.IGC.OR.IG2.LE.IGC) GO TO 9110 CALL DIPOT DO 9100 IT=1,N A1=SF(IT,1) DO 9100 IS=1,N 9100 Y(IT,IS)=Y(IT,IS)-A1*DI(IS) 9110 CONTINUE C ----END INSERT A2=-2.*GG DO 150 IT=1,N A1=A2*P(IT,IG2) DO 150 IS=1,N 150 Y(IT,IS)=Y(IT,IS)*A1 J=1 IF ((LAM+LQN(IG2)-LL(I1)).NE.0) J=2 LI=(I1-1)*N2-1+J LJ=(I2-1)*N2 IF(J.EQ.2) GO TO 158 LI=LI+1 LIP=LJ DO 155 IS=1,N LIP=LIP+1 155 C(LI,LIP)=C(LI,LIP)+Y(1,IS) 158 DO 160 IT=2,ITC LI=LI+1 LIP=LJ A1=R(IT)**(LQN(IG2)+1) DO 160 IS=1,N LIP=LIP+1 160 C(LI,LIP)=C(LI,LIP)+A1*Y(IT,IS) DO 163 IT=ITCP1,N LI=LI+1 LIP=LJ DO 163 IS=1,N LIP=LIP+1 163 C(LI,LIP)=C(LI,LIP)+Y(IT,IS) IF (I.EQ.IP) GO TO 50 I1=IP I2=I IG1=IGP 140 IG2=IG GO TO 50 C C A COEFFICIENT 80 J=J-NCHF IG=IP LI=(I-1)*N2+1 LJ=LQN(IG)+1 LIP=NCHF*N2+J C(LI,LIP)=C(LI,LIP)+AA*Q(1,IG) DO 170 IT=2,ITC LI=LI+1 170 C(LI,LIP)=C(LI,LIP)+AA*Q(IT,IG)*R(IT)**LJ DO 180 IT=ITCP1,N LI=LI+1 180 C(LI,LIP)=C(LI,LIP)+AA*Q(IT,IG) LI=(I-1)*N2 IPQ=-1 CALL LAP DO 190 IS=1,N LI=LI+1 190 C(LIP,LI)=C(LIP,LI)+AA*FLAP(IS) GO TO 50 C C B COEFFICIENT 90 J=J-NCHF MOAD=1 CALL POT IGPP=IGP IGP=IP C ----INSERT FOR DI-ELECTRONIC POLARIZATION IJK=0 IF(NDIEL.EQ.0) GO TO 9140 IF(LAM.NE.1) GO TO 9140 IF(IG.LE.IGC.OR.IGP.LE.IGC.OR.IGPP.LE.IGC ) GO TO 9140 IJK=1 IG1=IG L=LQN(IGPP) CALL DIPOT A1=0. DO 9120 IT=1,N 9120 A1=A1+DI(IT)*P(IT,IGPP) DO 9130 IT=1,N 9130 YPOT(IT)=YPOT(IT)-A1*SF(IT,1) 9140 CONTINUE C ----END INSERT CB=BB+BB DO 200 IT=1,N 200 YPOT(IT)=CB*YPOT(IT)*P(IT,IGP) LI=(I-1)*N2+1 LIP=NCHF*N2+J I1=LAM+LQN(IGP)-LL(I) IF (I1.EQ.0) C(LI,LIP)=C(LI,LIP)+YPOT(1) DO 210 IT=2,ITC LI=LI+1 210 C(LI,LIP)=C(LI,LIP)+YPOT(IT)*R(IT)**(LQN(IGP)+1) DO 220 IT=ITCP1,N LI=LI+1 220 C(LI,LIP)=C(LI,LIP)+YPOT(IT) L=LL(I) CALL RLAM C ----INSERT FOR DI-ELECTRONIC POLARIZATION IF(IJK.NE.1) GO TO 9160 IG1=IGP CALL DIPOT DO 9150 IT=1,N 9150 RL(IT)=RL(IT)-A1*DI(IT) 9160 CONTINUE C ----END INSERT LI=(I-1)*N2 CB=BB+BB DO 230 IT=1,N LI=LI+1 230 C(LIP,LI)=C(LIP,LI)+CB*RL(IT) GO TO 50 C C H(J,JP) COEFFICIENT 95 LI=NCHF*N2+I LIP=NCHF*N2+IP C(LI,LIP)=C(LI,LIP)+FF IF (I.EQ.IP) GO TO 50 C(LIP,LI)=C(LIP,LI)+FF GO TO 50 C C KEY = 6 CARD READ 100 CONTINUE C C ----INSERT FOR DI-ELECTRONIC POLARIZATION C MODIFY MATRIX ELEMENT SCRIPT H(J,JP) C 2.*AA IS COEFFICIENT OF R(IG,IGP,IGPP,IGPPP,LAM=1) C TERMINATE WITH KEY=8 IF(NDIEL.EQ.0) GO TO 101 C CHECK VALUE OF LAM READ ON KEY = 6 CARD IF(LAM.EQ.0) WRITE(6,613) 9165 READ(IFILE,9510)KEY,J,JP,IG,IGP,IGPP,IGPPP,AA IF(KEY.EQ.8) GO TO 101 IF(MDE.EQ.1) GO TO 9165 IG1=IG L=LQN(IGPP) CALL DIPOT A1=0. DO 9170 IT=1,N 9170 A1=A1+DI(IT)*P(IT,IGPP) IG1=IGP L=LQN(IGPPP) CALL DIPOT A2=0. DO 9180 IT=1,N 9180 A2=A2+DI(IT)*P(IT,IGPPP) A1=2.*A1*A2*AA LI=NCHF*N2+J LIP=NCHF*N2+JP C(LI,LIP)=C(LI,LIP)-A1 IF(J.EQ.JP) GO TO 9165 C(LIP,LI)=C(LIP,LI)-A1 GO TO 9165 C ----END INSERT C C FOR NDIEL = 0 ,CHECK VALUE OF LAM READ ON KEY = 6 CARD 101 IF(NDIEL.EQ.1) GO TO 103 IF(LAM.EQ.0) GO TO 103 WRITE(6,614) LAM 102 READ(IFILE,506) KEY IF(KEY.EQ.8) GO TO 103 GO TO 102 C INSERT TERMS FOR IT=N1,N2 103 RTN1=1./R(N1) RTN2=1./R(N2) A1=RTN1 A2=RTN2 DO 235 JLAM=1,5 LAM=JLAM A1=A1*RTN1 A2=A2*RTN2 RITN1(LAM)=A1 235 RITN2(LAM)=A2 RTN1=2.*RTN1*NLCTRN RTN2=2.*RTN2*NLCTRN N12=N1-N2 I1=N12 DO 240 I=1,NCHF I1=I1+N2 I2=N12 DO 240 IP=1,NCHF I2=I2+N2 IF(I.NE.IP) GO TO 248 C(I1,I2)=RTN1 C(I1+1,I2+1)=RTN2 248 A1=0. A2=0. DO 250 JLAM=1,5 LAM=JLAM A3=BLAM(I,IP,LAM) A1=A1+A3*RITN1(LAM) 250 A2=A2+A3*RITN2(LAM) C(I1,I2)=C(I1,I2)+A1 240 C(I1+1,I2+1)=C(I1+1,I2+1)+A2 C C ----INSERT FOR DI-ELECTRONIC POLARISATION C WRITE AND PUNCH INFORMATION IF(NDIEL.EQ.0) GO TO 9190 WRITE(6,9650) IGC,ALPHA,RDICUT NP=0 IF(JFILE.NE.0) WRITE(JFILE,700) NP,IGC,ALPHA,RDICUT 9190 CONTINUE C ----END INSERT C C ----INSERT FOR CORE POLARIZATION AND EXTRA CHANNEL POLARIZATION IF(NPOL.EQ.0) GO TO 9250 NDI=NDIEL NDIEL=0 C WITH NDIEL=0, FUNSF(SR)=ALPHA*(1.-EXP(-SR/RDICUT)**6)/(SR**4) NP=(NPOL+1)/2 IF(LAMMX.LT.3) LAMMX=3 DO 9240 IJM=1,NP IF(IJM.EQ.1) WRITE(6,9620) IF(IJM.EQ.2) WRITE(6,9630) WRITE(6,9635) DO 9230 IP=1,NCHF C READ INFORMATION ON POLARIZATION C C READ - C ====== C (SEE CPC SECTION 13.4(I)) C I = CHANNEL INDEX C ALPHA = CORE POLARIZABILITY IN CHANNEL I C RCUT = CUT-OFF PARAMETER FOR CORE POLARIZATION C IN CHANNEL I C READ(5,9500) I,ALPHA,RCUT C PUNCH INFORMATION ON POLARIZATION IF(JFILE.NE.0) WRITE(JFILE,700) IJM,I,ALPHA,RCUT WRITE(6,9640) I,ALPHA,RCUT RDICUT=RCUT LLIP1=LL(I)+1 LI=(I-1)*N2+1 DO 9210 IT=2,ITC LI=LI+1 RIT=R(IT) 9210 C(LI,LI)=C(LI,LI)-FUNSF(RIT)*RIT**LLIP1 DO 9220 IT=ITCP1,N2 LI=LI+1 9220 C(LI,LI)=C(LI,LI)-FUNSF(R(IT)) 9230 BLAM(I,I,3)=BLAM(I,I,3)-ALPHA 9240 CONTINUE NDIEL=NDI 9250 CONTINUE C ----END INSERT C C COMPUTATION OF FIRST FORM OF C NOW COMPLETE C IF REQUIRED, WRITE C ON UNIT 11 IF(IPUNCH.NE.1) GO TO 9255 IF(MODE.EQ.1) GO TO 9255 REWIND 11 WRITE(11) C 9255 CONTINUE C C INSERT DERIVATIVE TERMS, TO OBTAIN C AS DEFINED AT THE END OF C SECTION 6 OF CPC C C INCORPORATE S AND T DO 260 I=1,NCHF FL=LL(I) LI=I*N2 DO 260 K=ITCP1,N2 IT=N2+ITCP1-K DO 270 LIP=1,NTOT A1=0. DO 280 J=1,3 IP=LI-3+J 280 A1=A1+C(IP,LIP)*T(IT,J) 270 C(LI,LIP)=-A1 DO 290 J=1,3 LIP=LI-3+J IS=IT-3+J A1=R(IS)*R(IS) A1=S(IT,J)-FL*(FL+1.)*T(IT,J)/A1 IF (IS.LE.ITC) A1=A1*R(IS)**(LL(I)+1) 290 C(LI,LIP)=C(LI,LIP)+A1 LI=LI-1 260 CONTINUE C C MODIFY MATRIX TO ACT ON GS INSTEAD OF F , WHERE , C F(IT)=PHI(IT)+R(IT)**2*GS(IT) FOR IT=1,ITC C F(IT)=GS(IT) FOR IT=ITCP1,N2 LI=-N2 DO 300 I=1,NCHF FL=LL(I) A1=Z/(FL+1.) DO 310 IT=1,ITC 310 PHI(IT)=1.-A1*R(IT) LI=LI+N2 DO 300 IT=1,NTOT A1=0. LIP=LI DO 320 IS=1,ITC LIP=LIP+1 A1=A1+C(IT,LIP)*PHI(IS) 320 C(IT,LIP)=C(IT,LIP)*R(IS)*R(IS) 300 G(IT,I)=-A1 C C INSERT DERIVATIVE TERMS FOR IT=1,ITC A4=2.*Z*Z LI=-N2 DO 330 I=1,NCHF LI=LI+N2 L=LL(I) FL=L LIP=LI+1 C(LIP,LIP)=C(LIP,LIP)-(4.*FL+6.) A1=A4/(FL+1.) G(LIP,I)=G(LIP,I)-A1 DO 330 IT=2,ITC LIP=LI+IT A2=R(IT)**(L+1) G(LIP,I)=G(LIP,I)-A2*A1 DO 330 IS=1,ITC IP=LI+IS A3=FL*SD1(IT,IS)+SD2(IT,IS) 330 C(LIP,IP)=C(LIP,IP)-A2*A3 C C MODIFY BLAM FOR INPUT INTO ASYMPT DO 340 I=1,NCHF FL=LL(I) DO 350 IP=1,NCHF DO 350 JLAM=1,LAMMX 350 BLAM(I,IP,JLAM)=-BLAM(I,IP,JLAM) 340 BLAM(I,I,1)=BLAM(I,I,1)-FL*(FL+1.) C PUNCH BLAM IF(MODE.EQ.1) RETURN IF(IPUNCH.NE.1) RETURN WRITE(JFILE,701) LAMMX DO 360 LLL=1,LAMMX DO 360 J=1,NCHF 360 WRITE(JFILE,703) (BLAM(I,J,LLL),I=1,NCHF) RETURN C C PROCEDURE FOR SKIPPING CASE IF DIMENSIONS EXCEEDED 400 READ(IFILE,506) KEY,LAM IF (KEY.EQ.6) GO TO 405 GO TO 400 405 IF(LAM.EQ.0) GO TO 410 407 READ(IFILE,504) KEY IF(KEY.EQ.8) GO TO 410 GO TO 407 410 IF(NPOL.EQ.0) GO TO 412 NP=(NPOL+1)/2 DO 408 IP=1,NP DO 408 J=1,NCHF 408 READ(5,9500) I,ALPHA,RCUT 412 NCHF=0 RETURN C C READ FORMATS 501 FORMAT (5X,2I3,I5,I3, 3X,2I3,6X,2A3) 502 FORMAT (13X,I7,2I3,9X,15I2) 504 FORMAT (I5,6X,I5,6X,4I3,I4,E14.6) 505 FORMAT (13X,I7,15X,15I2) 506 FORMAT(I5,29X,I4) 9500 FORMAT(I10,2F10.4) 9510 FORMAT(I5,5X,2I5,10X,4I5,E20.6) C C PRINT FORMATS 600 FORMAT(1H1,32HSTATE OF ELECTRON-ATOM SYSTEM /1X,32(1H*)//) 1600 FORMAT(/////,33H STATE OF ELECTRON-ATOM SYSTEM /1X,32(1H*)//) 601 FORMAT(/26H NUMBER OF FREE CHANNELS =,I3/ C 27H NUMBER OF BOUND CHANNELS =,I3// C 13H CHANNEL LIST/1X,12(1H-)// C 14H FREE CHANNELS//2X,1HI,4X,4HTERM,4X,7HSMALL L,4X, C 59HCHANNEL RADIAL FUNCTION ORTHOGONAL TO P(IT,IGAM) FOR IGAM = ) 602 FORMAT(I3,5X,I2,6X,I3,6X,25I4/(24X,25I4)) 603 FORMAT(//15H BOUND CHANNELS//3H J,8X,14HCONFIGURATION,/ C 11X,5HIGAM=,I2,25I4/(16X,I2,25I4)) 604 FORMAT(I3,13X,I2,25I4/(15X,I2,25I4)) 606 FORMAT(10H (2*S+1) =,I2,6H, L =,I2,13H, EVEN PARITY) 607 FORMAT(10H (2*S+1) =,I2,6H, L =,I2,12H, ODD PARITY) 9620 FORMAT(//18H CORE POLARIZATION/1X,17(1H-)/) 9630 FORMAT(//27H EXTRA CHANNEL POLARIZATION/1X,26(1H-)/) 9635 FORMAT(9H CHANNEL,5X,14HPOLARIZABILITY,5X,4HRCUT) 9640 FORMAT(I3,9X,F8.4,13X,F6.4) 9650 FORMAT(//27H DI-ELECTRONIC POLARISATION/1X,26(1H-)//6H IGC=,I2, C8H, ALPHA=,F8.4,9H, RDICUT=,F6.4) C C PUNCH FORMATS 700 FORMAT(I2,I3,2E14.7) 701 FORMAT(24I3) 702 FORMAT(1X,7HSLPI = ,I3,1H,,I3,15H FOR IMPACT ON ,A4,2A3) 703 FORMAT(5E14.7) C C FORMATS FOR MESSAGES 605 FORMAT(/8H ITMCUT=,I3,5H THE,I3, C 26H LOWEST TERMS ARE RETAINED) 608 FORMAT(47H CHANNELS WITH SMALL L GREATER THAN OR EQUAL TO,I2, 115H ARE SUPPRESSED ) 610 FORMAT (1H0,34HERROR - DIMENSIONS EXCEEDED. NTOT= ,I5, 134H GREATER THAN NTOTMX. CASE NOT RUN ) 611 FORMAT(26H0NUMBER OF FREE CHANNELS =,I3, C 52H WHICH IS GREATER THAN MAXIMUM ALLOWED BY DIMENSIONS) 612 FORMAT(43H ITMCUT NEGATIVE MEANS ADDITIONAL CHANNELS, C 25H MAY HAVE BEEN SUPPRESSED) 613 FORMAT(/,54H *** WARNING *** FOR NDIEL = 1, NON-ZERO VALUE OF LAM C,31H SHOULD BE READ ON KEY = 6 CARD,/,18X,19HVALUE READ IS LAM = C ,51H0. ALGEBRA FILE MAY NOT CONTAIN KEY = 7 AND 8 CARDS) 614 FORMAT(/,24H FOR NDIEL = 0 AND LAM =,I4, C60H ON KEY=6 CARD, KEY=7 AND 8 CARDS ARE READ BUT NOT PROCESSED/) 615 FORMAT(/,33H COLLISION ALGEBRA CREATED USING ,2A3) END C*********************************************************************** SUBROUTINE EXPAN(IAS,EROR,RO,ISF,RF,IBUG,IOUT) C INCLUDE 'PARAM' C C PHASE FROM DIAGONAL CHANNEL AND AMPLITUDE BY ASYMPTOTIC EXPANSION C IAS=MAXIMUM NUMBER OF TERMS IN EXPANSION C EROR=MAXIMUM VALUE OF LAST TERM IN EXPANSION C RO=DESIRED VALUE OF R C ISF,IAS OR POINT WHERE EXPANSION DIVERGES C RF=MINIMUM VALUE OF R FOR WHICH EXPANSION VALID FOR GIVEN ERROR C DIMENSION ISF(MZ07) C C COMMON NCHD1,IJNZ(MZ25),KNZ(5),BLIN(MZ25),NZD,NZDM, C AB1(MZ07,13,10),AB2(MZ07,13,10),ALP(13,13),CC(MZ07,MZ07,5), C E2(MZ07),PP(12,10),T1(MZ07),V(MZ11,11,10),ZET(12,10),ZTIN(12,10) C ,GAM1(MZ07,15),GAM2(MZ07,15) C , VF1(MZ07),VF2(MZ07),VFDR1(MZ07),VFDR2(MZ07) COMMON /CMN02/ NCNL,NOP,LAMAX COMMON /CMN03/ NDL(MZ07), NDU(MZ07) COMMON /CMN04/ HY,NHY,JJ NDLJJ=NDL(JJ) NDUJJ=NDU(JJ) E2JJ=E2(JJ) EKAP=SQRT(ABS(E2JJ)) RA=RO IF(NZD.EQ.0) GO TO 12 ALZ=NZD ETA=1./EKAP A5=ETA*ETA ETA2=A5*0.5 GO TO 13 12 ALZ=1. ETA=0. A5=0. ETA2=0. 13 RP1S=RO*ALZ RP2S=1.0E+08*ALZ EKAP2=EKAP+EKAP IF(JJ.GT.NOP) GO TO 1000 C C OPEN CHANNELS FOR Z GREATER THAN OR EQUAL TO ZERO C DO 20 J=1,NCNL GAM1(J,1)=0. GAM1(J,2)=0. 20 GAM2(J,1)=0. GAM1(JJ,1)=1. GAM1(JJ,2)=-ETA2 C DO 35 J=1,NCNL G3=0. IF(J.LT.NDLJJ) GO TO 35 IF(J.GT.NDUJJ) GO TO 35 G3=-CC(J,JJ,1)*0.5/EKAP IF(J.NE.JJ) GO TO 35 G3=G3+ETA2/EKAP 35 GAM2(J,2)=G3 C A2=0. A3=A5 A6=ETA ETA2=ETA+ETA M1=1 M3=0 DO 100 I=3,IAS A4=A3 A5=A6 M2=M1 M1=M1+1 M3=M3+2 A1=1./(FLOAT(M1+M1)*EKAP) A2=A2+EKAP2 A3=A3-M3 A6=A6+ETA2 MDG=MIN(M1,LAMAX) MNDG=MIN(M2,LAMAX) DO 100 J=1,NCNL G1=GAM1(J,M1) G2=GAM2(J,M1) IF(J.LT.NDLJJ) GO TO 50 IF(J.GT.NDUJJ) GO TO 50 C C DEGENERATE ENERGIES DUMA=A3*G2+A6*G1 DUMB=A3*G1-A6*G2 DO 40 L=1,NCNL IMIN=I DO 40 M=1,MDG IMIN=IMIN-1 G3=CC(J,L,M) DUMA=DUMA-G3*GAM2(L,IMIN) 40 DUMB=DUMB-G3*GAM1(L,IMIN) GAM1(J,I)=-DUMA*A1 GAM2(J,I)=DUMB*A1 GO TO 100 C C NON-DEGENERATE ENERGIES 50 G3=GAM1(J,M2) G4=GAM2(J,M2) DUMA=A4*G3-A5*G4-A2*G2 DUMB=A4*G4+A5*G3+A2*G1 DO 70 L=1,NCNL IMIN=M1 DO 70 M=1,MNDG IMIN=IMIN-1 G3=CC(J,L,M) DUMA=DUMA-G3*GAM1(L,IMIN) 70 DUMB=DUMB-G3*GAM2(L,IMIN) G3=1./(E2(J)-E2JJ) GAM1(J,I)=DUMA*G3 GAM2(J,I)=DUMB*G3 100 CONTINUE C DO 150 J=1,NCNL RP1=RP1S RP2=RP2S DO 130 I=4,IAS IT=I-1 ASM=AMAX1(ABS(GAM1(J,I)),ABS(GAM2(J,I))) IF(ASM.LT.0.1) GO TO 130 RP1=(ASM/EROR)**(1./FLOAT(IT)) IF(RP1.GT.RP2) GO TO 140 RP2=RP1 130 CONTINUE ISF(J)=IAS RF=RP1/ALZ GO TO 150 140 ISF(J)=IT RF=RP2/ALZ 150 RA=AMAX1(RF,RA) GO TO 250 C C CLOSED CHANNELS C 1000 DO 920 J=1,NCNL 920 GAM1(J,1)=0. GAM1(JJ,1)=1. C DO 835 J=1,NCNL G3=0. IF(J.LT.NDLJJ) GO TO 835 IF(J.GT.NDUJJ) GO TO 835 G3=-CC(J,JJ,1)*0.5/EKAP IF(J.NE.JJ) GO TO 835 G3=G3+(ETA*0.5-ETA2)/EKAP 835 GAM1(J,2)=G3 ETAP2=ETA+ETA C M1=1 A2=0. M3=0 A6=ETA-A5 DO 900 I=3,IAS M2=M1 A5=A6 M1=M1+1 M3=M3+2 A1=1./(FLOAT(M1+M1)*EKAP) A6=A6+ETAP2-M3 A2=A2-EKAP2 MDG=MIN(M1,LAMAX) MNDG=MIN(M1-1,LAMAX) DO 900 J=1,NCNL IF(J.LT.NDLJJ) GO TO 950 IF(J.GT.NDUJJ) GO TO 950 C C TERMS WITH DEGENERATE ENERGIES DUMA=GAM1(J,M1)*A6 DO 840 L=1,NCNL IMIN=I DO 840 M=1,MDG IMIN=IMIN-1 840 DUMA=DUMA-CC(J,L,M)*GAM1(L,IMIN) GAM1(J,I)=A1*DUMA GO TO 900 C C TERMS WITH NON-DEGENERATE ENERGIES 950 DUMA=A2*GAM1(J,M1)+A5*GAM1(J,M2) DO 870 L=1,NCNL IMIN=M1 DO 870 M=1,MNDG IMIN=IMIN-1 870 DUMA=DUMA-CC(J,L,M)*GAM1(L,IMIN) GAM1(J,I)=DUMA/(E2(J)-E2JJ) 900 CONTINUE C DO 1050 J=1,NCNL RP1=RP1S RP2=RP2S DO 1030 I=4,IAS IT=I-1 ASM=ABS(GAM1(J,I)) IF(ASM.LT.0.1) GO TO 1030 RP1=(ASM/EROR)**(1./FLOAT(IT)) IF(RP1.GT.RP2) GO TO 1040 RP2=RP1 1030 CONTINUE ISF(J)=IAS RF=RP1/ALZ GO TO 1050 1040 ISF(J)=IT RF=RP2/ALZ 1050 RA=AMAX1(RF,RA) C C CALCULATE AMPLITUDE PARAMETERS IN CLOSED CHANNEL CASE RHO=RA*ALZ RZZ=1./RHO DO 392 I=1,NCNL ISG=ISF(I) ISG1=ISG-1 A4=GAM1(I,1) A5=0. A1=1. A2=0. DO 390 LL=2,ISG1 A1=A1*RZZ A2=A2+RZZ A3=GAM1(I,LL)*A1 A4=A4+A3 390 A5=A5-A3*A2 A1=A1*RZZ*0.5 A2=A2+RZZ A3=GAM1(I,ISG)*A1 VF1(I)=A4+A3 392 VFDR1(I)=A5-A3*A2 C 250 RF=RA IF(IBUG.EQ.1) WRITE(IOUT,801) JJ,RF,(ISF(I),I=1,NCNL) 801 FORMAT(34H TRIED EXPANSION METHOD IN CHANNEL ,I3,13H, CONVERGENCE C ,E12.5,/,9H ISF(I) =,20I5) RETURN END C********************************************************************** SUBROUTINE FG C INCLUDE 'PARAM' C DESCRIBED IN CPC SECTION 10.3 C C FOR MFG NON-ZERO, ASYMPTOTIC POINTS HAVE BEEN REDEFINED AS C R1=R(N1)+MFG*A1, C R2=R1+A1, WHERE A1=0.5*(R(N2)-R(N1)) C ASYSM IS USED TO CALCULATE FUNCTIONS AT R1,R2. FOX-GOODWIN C TECHNIQUE IS USED TO INTEGRATE FROM R1,R2 TO R(N1),R(N2) C REFERENCE - NORCROSS AND SEATON, J. PHYS. B.,6,614,1973 C COMMENTS REFER TO THEIR NOTATION - C P AND Q ARE DEFINED BY THEIR EQN. (3.16). RFG IS MATRIX R IN THEIR C EQN. (3.7) AND SR IS MATRIX SCRIPT R IN THEIR EQN. (3.11). C C INTEGER HEAD COMMON NCHD1,IJNZ(MZ25),KNZ(5),BLIN(MZ25) C , RFGST(MZ22,MZ07,MZ07) COMMON /HOLD/DEL1,DEL2,DELNU,ELOW,EMAX,FINT,FRAT,HEAD,IDM,SZ,Z, C BLAM(MZ26),EKLSQ(MZ01),ET(MZ07),ETERM(MZ02), C F1(2,MZ07,MZ07),F2(2,MZ07,MZ07),FI(MZ05,MZ05),FLAG(MZ09,3,4), C FX(2,MZ07,MZ07),P(MZ06,MZ03),Q(MZ06,MZ03),R(MZ06),S(MZ10,3), C SD1(MZ05,MZ05),SD2(MZ05,MZ05),SH(4,5),T(MZ10,3), C IFILE,IPRINT,ISTOT,ITC,ITCP1,ITMX,JFILE,KCASE,LAMMX,LCUT,LTOT, C MCH,MDG,MODE,N,N1,N2,NCHB,NCHCL,NCHF,NCHFMX,NCHOP,NLCTRN,NLOW, C NMX,NRAD,NTERM,NTHROW,NTOT,NTOTMX, C ISTERM(MZ02),ITERM(MZ07),LL(MZ07),LQN(MZ03),LTERM(MZ02), C MDEGEN(MZ01),MTCHCL(MZ01),NQN(MZ03) COMMON /WORK2/A(MZ07,MZ07),RFG(MZ07,MZ07),RFGP(MZ07,MZ07) X ,SR(MZ07,MZ07) DIMENSION W1(MZ07,MZ07),W2(MZ07,MZ07),DD(MZ07,MZ07) COMMON /FGCOM/R1,R2,KJ(MZ07),METHOD,MFG COMMON /OUT/IPUNCH,ICONV DIMENSION X(MZ07) COMMON /CMN03/ NDL(MZ07), NDU(MZ07) EQUIVALENCE (DD(1,1),X(1)) C LK=0 K2=0 DO 22 K=1,LAMMX K1=K2+1 K2=KNZ(K) IF(K2.LT.K1) GO TO22 DO 21 L=K1,K2 IJ=IJNZ(L) 21 BLIN(L)=BLAM(IJ+LK) 22 LK=LK+MZ24 C DR=R2-R1 CFG=12./(DR*DR) MFGP1 = MFG + 1 CRYA ICRAY=1 CRYA ND= 100 C C INCLUDE CLOSED CHANNEL FUNCTIONS IN IRREGULAR SOLUTIONS IF (NCHOP.EQ.0.OR.NCHCL.EQ.0) GO TO 10 M=NCHOP+1 DO 5 I=1,NCHF DO 5 J=M,NCHF F1(2,I,J)=F1(1,I,J) 5 F2(2,I,J)=F2(1,I,J) C 10 DO 200 K=1,2 C INCREMENT ENERGY FOR NCHOP=0 AND K=2 IF(K.EQ.1) GO TO 20 IF(NCHOP.GT.0) GO TO 20 DO 15 I=1,NCHF 15 ET(I)=ET(I)+DEL1 C INITIALIZE SR AND RFG 20 DO 25 I=1,NCHF DO 25 J=1,NCHF 25 A(I,J)=F1(K,I,J) CALL INV(NCHF) DO 35 I=1,NCHF DO 35 J=1,NCHF A1=0. DO 30 IP=1,NCHF 30 A1=A1+F2(K,I,IP)*A(IP,J) SR(I,J)=A1 35 RFG(I,J)=A1 C STORE RFG IF IPUNCH NON-ZERO IF (IPUNCH.EQ.0) GO TO 352 MST = MFGP1 DO 351 I=1,NCHF DO 351 J=1,NCHF 351 RFGST(MST,I,J) = RFG(I,J) C INITIALIZE W 352 CALL DBLU(W1,R1) CALL DBLU(W2,R2) C GENERATE RFG AND SR DO 90 IS=1,MFG C IF IS=MFG, STORE OLD RFG AS RFGP IF(IS.NE.MFG) GO TO 37 DO 36 I=1,NCHF DO 36 J=1,NCHF 36 RFGP(I,J)=RFG(I,J) C C COMPUTE NEW RFG 37 DO 3 I=1,NCHF DO 3 J=1,NCHF A1=-10.*W1(I,J)-CFG*RFG(I,J) IF(I.EQ.J) A1=A1+2.*CFG DO 2 IP=1,NCHF 2 A1=A1-W2(I,IP)*RFG(IP,J) 3 A(I,J)=A1 CRAY DO 2 J=1,NCHF CRAY2 A(I,J)=-10.*W1(I,J)-CFG*RFG(I,J) CRAY A(I,I)=A(I,I)+2.*CFG CRAY DO 3 IP=1,NCHF CRAY DO 3 J=1,NCHF CRAY3 A(I,J)=A(I,J)-W2(I,IP)*RFG(IP,J) CALL INV(NCHF) C COMPUTE NEW W1,W2 AND R1 CRYA IF(ICRAY.EQ.1) GO TO 8 DO 55 I=1,NCHF DO 55 J=1,NCHF 55 W2(I,J)=W1(I,J) CRYA8 CALL SCOPY(ND,W1,1,W2,1) R1=R1-DR CALL DBLU(W1,R1) C COMPLETE COMPUTATION OF RFG DO 4 I=1,NCHF DO 4 J=1,NCHF A1=A(I,J)*CFG DO 1 IP=1,NCHF 1 A1=A1+A(I,IP)*W1(IP,J) 4 RFG(I,J)=A1 CRAY DO 1 J=1,NCHF CRAY1 RFG(I,J)=A(I,J)*CFG CRAY DO 4 IP=1,NCHF CRAY DO 4 J=1,NCHF CRAY4 RFG(I,J)=RFG(I,J)+A(I,IP)*W1(IP,J) C NEW RFG NOW COMPUTED C C STORE RFG IF IPUNCH NON - ZERO IF (IPUNCH.EQ.0) GO TO 67 MST=MST-1 DO 66 I=1,NCHF DO 66 J=1,NCHF 66 RFGST(MST,I,J) = RFG(I,J) C C COMPUTE NEW SR 67 DO 7 I=1,NCHF DO 7 J=1,NCHF A1=0. DO 6 IP=1,NCHF 6 A1=A1+SR(I,IP)*RFG(IP,J) 7 DD(I,J)=A1 CRAY DO 6 J=1,NCHF CRAY6 DD(I,J)=0. CRAY DO 7 IP=1,NCHF CRAY DO 7 J=1,NCHF CRAY7 DD(I,J)=DD(I,J)+SR(I,IP)*RFG(IP,J) CRYA IF(ICRAY.EQ.1) GO TO 9 DO 80 I=1,NCHF DO 80 J=1,NCHF 80 SR(I,J)=DD(I,J) CRYA9 CALL SCOPY(ND,DD,1,SR,1) C NEW SR NOW COMPUTED C 90 CONTINUE C INWARD INTEGRATION NOW COMPLETED C C COMPUTE NEW F1 DO 100 I=1,NCHF DO 100 J=1,NCHF A1=0. IF(I.EQ.J) A1=1. 100 F1(K,I,J)=A1 C F1 NOW INITIALIZED TO UNIT MATRIX IF(NCHOP.EQ.0) GO TO 170 IF (NCHCL.EQ.0) GO TO 118 C COMPUTE Q DO 105 I=1,NCHCL DO 105 J=1,NCHCL 105 A(I,J)=F2(K,I+NCHOP,J+NCHOP) CALL INV(NCHCL) DO 115 I=1,NCHOP DO 115 J=1,NCHCL A1=0. DO 110 IP=1,NCHCL 110 A1=A1+F2(K,I,IP+NCHOP)*A(IP,J) 115 W1(I,J)=A1 C Q=W1 NOW COMPUTED C COMPUTE P 118 DO 125 I=1,NCHOP DO 125 J=1,NCHOP A1=-SR(I,J) IF (NCHCL.EQ.0 ) GO TO 125 DO 120 IP=1,NCHCL 120 A1=A1+W1(I,IP)*SR(IP+NCHOP,J) 125 A(I,J)=-A1 CALL INV(NCHOP) C P=A NOW COMPUTED C COMPUTE OPEN-OPEN AND OPEN-CLOSED ELEMENTS OF F1 DO 135 I=1,NCHOP DO 135 J=1,NCHOP A1=F2(K,I,J) IF (NCHCL.EQ.0 ) GO TO 135 DO 130 IP=1,NCHCL 130 A1=A1-W1(I,IP)*F2(K,NCHOP+IP,J) 135 W2(I,J)=A1 C F(OPEN-OPEN)-Q*F(CLOSED-OPEN)=W2 NOW COMPUTED DO 145 I=1,NCHOP DO 145 J=1,NCHOP A1=0. DO 140 IP=1,NCHOP 140 A1=A1+A(I,IP)*W2(IP,J) 145 F1(K,I,J)=A1 C NEW F1(OPEN-OPEN) NOW COMPUTED IF (NCHCL.EQ.0) GO TO 170 DO 155 I=1,NCHOP DO 155 J=1,NCHCL A1=SR(I,NCHOP+J) DO 150 IP=1,NCHCL 150 A1=A1-W1(I,IP)*SR(IP+NCHOP,J+NCHOP) 155 W2(I,J)=-A1 C Q*SR(CLOSED-CLOSED) - SR(OPEN-CLOSED) = W2 NOW COMPUTED DO 165 I=1,NCHOP DO 165 J=1,NCHCL A1=0. DO 160 IP=1,NCHOP 160 A1=A1+A(I,IP)*W2(IP,J) 165 F1(K,I,NCHOP+J)=A1 C NEW F1(OPEN-CLOSED) NOW COMPUTED 170 CONTINUE C NEW F1 NOW COMPUTED C C TRANSFORMATION OF F1 FOR NCHOP.EQ.0 C THE METHOD USED (SEE CPC SECTION 10.3(III)) DEPENDS ON C ARRAY KJ(J) COMPUTED IN ASYMPT. C NOTE THAT FINAL NORMALIZATION (CPC EQUATIONS (10.14) OR (7.7)) C IS PERFORMED IN ASYMPT AFTER RETURN FROM SUBROUTINE FG. IF(NCHOP.NE.0) GO TO 1790 CARD DELETED-SEE CPC NOTE IN PROOF C IF KJ(J) = 2 FOR ANY J, TRANSFORM F1 FOR K=1 AND K=2 IF(K.EQ.2) GO TO 1715 METHOD=1 DO 1705 J=1,NCHF IF(KJ(J).EQ.2) GO TO 1710 1705 CONTINUE C METHOD=1 IS USED. NO CHANGE REQUIRED AT THIS STAGE GO TO 1790 C METHOD=2 IS USED. 1710 METHOD=2 1715 IF(METHOD.EQ.1) GO TO 1790 C PROCEDURE FOR METHOD=2 IS DESCRIBED IN CPC SECTION 10.3 (III). C STORE F2 IN A AND INVERT DO 1720 I=1,NCHF DO 1720 J=1,NCHF 1720 A(I,J)=F2(K,I,J) CALL INV(NCHF) C COMPUTE A*SR AND STORE IN DD DO 1730 I=1,NCHF DO 1730 J=1,NCHF A1=0. DO 1725 L=1,NCHF 1725 A1=A1+A(I,L)*SR(L,J) 1730 DD(I,J)=A1 C SOLVE THE EQUATIONS DD*UPPERT=UNITLT WHERE - C UPPERT IS UPPER TRIANGULAR C UNITLT IS UNIT LOWER TRIANGULAR. C STORE UPPERT AND LOWER PART OF UNITLT IN A. DO 1765 J=1,NCHF IF(J.EQ.1) GO TO 1740 J1=J-1 DO 1735 I=J,NCHF A1=DD(J1,I) DO 1735 L=J,NCHF 1735 DD(L,I)=DD(L,I)-DD(L,J1)*A1 1740 A1=1./DD(J,J) A(J,J)=A1 IF(J.EQ.NCHF) GO TO 1750 J2=J+1 DO 1745 I=J2,NCHF DD(J,I)=A1*DD(J,I) 1745 A(I,J)=DD(I,J)*A1 1750 IF(J.EQ.1) GO TO 1765 I=J DO 1760 L=1,J1 I2=I I=I-1 A1=0. DO 1755 L2=I2,J 1755 A1=A1-DD(I,L2)*A(L2,J) 1760 A(I,J)=A1 1765 CONTINUE C SCAN FOR DEGENERATE GROUPS OF CHANNELS. J1=1 2000 J2=NDU(J1) C IF KJ(J) = 2 FOR ANY MEMBER OF GROUP, USE CASE (A) OF CPC C SECTION 10.3 (III). C OTHERWISE USE CASE (B), FOR WHICH NO CHANGE REQUIRED AT THIS STAGE. DO 2010 J=J1,J2 IF(KJ(J).EQ.2) GO TO 2020 2010 CONTINUE C CASE (B) IS USED. GO TO 1790 C CASE (A) IS USED. MAKE TRANSFORMATION REQUIRED IN ORDER TO C SATISFY CPC EQUATIONS (10.17), (10.18). 2020 IF(J1.EQ.J2) GO TO 2090 J2M=J2-1 DO 2080 J=J1,J2M JP=J+1 X(J)=1. DO 2040 I=JP,J2 IM=I-1 A1=0. DO 2030 L=J,IM 2030 A1=A1+A(I,L)*X(L) 2040 X(I)=-A1 DO 2060 I=1,J A1=A(I,J) DO 2050 L=JP,J2 2050 A1=A1+A(I,L)*X(L) 2060 F1(K,I,J)=A1 DO 2070 I=JP,J2 A1=0. DO 2065 L=I,J2 2065 A1=A1+A(I,L)*X(L) 2070 F1(K,I,J)=A1 2080 CONTINUE 2090 DO 2095 I=1,J2 2095 F1(K,I,J2)=A(I,J2) IF(J2.EQ.NCHF) GO TO 1790 J1=J2+1 GO TO 2000 1790 CONTINUE C C COMPUTE NEW F2 C COMPUTE RFGP*RFG=DD DO 180 I=1,NCHF DO 180 J=1,NCHF A1=0. DO 175 IP=1,NCHF 175 A1=A1+RFGP(I,IP)*RFG(IP,J) 180 DD(I,J)=A1 C COMPUTE NEW F2=DD*F1 DO 190 I=1,NCHF DO 190 J=1,NCHF A1=0. DO 185 IP=1,NCHF 185 A1=A1+DD(I,IP)*F1(K,IP,J) 190 F2(K,I,J)=A1 C NEW F2 NOW COMPUTED C R1=R1+MFG*DR C WRITE STORED RFG ON UNIT 11 IF ( IPUNCH.EQ.0 ) GO TO 200 DO 192 MST=1,MFGP1 DO 191 I=1,NCHF DO 191 J=1,NCHF 191 RFG(I,J) = RFGST(MST,I,J) 192 WRITE (11) RFG 200 CONTINUE C IF(NCHOP.GT.0) RETURN DO 210 I=1,NCHF 210 ET(I)=ET(I)-DEL1 RETURN END C*********************************************************************** FUNCTION FUNSF(SR) C INCLUDE 'PARAM' C DESCRIBED IN CPC SECTION 13.4 C C THE CORE POLARIZATION POTENTIAL IS TAKEN TO BE C FUNSF(SR)=(1.-EXP(-(SR/RCUT)**6))*ALPHA/(SR**4) C WHERE SR IS THE RADIAL VARIABLE. C NOTE THAT - C (1) THIS FUNCTION SUBROUTINE CAN BE CHANGED TO COMPUTE C ANY OTHER POLARIZATION POTENTIAL THE USER MAY PREFER C (2) FOR DI-ELECTRONIC POLARIZATION, FUNSF IS COMPUTED AS C SQRT OF ABOVE EXPRESSION C (3) FOR POLARIZATION TERMS OTHER THAN DI-ELECTRONIC, C THIS ROUTINE IS CALLED WITH C ALPHA=CHANNEL POLARIZATION C RCUT=CHANNEL CUT-OFF PARAMETER C IF NDIEL NON-ZERO, NDIEL TEMPORARILY SET TO ZERO FOR C THESE CASES C COMMON /POL1/ ALPHA,RCUT,NPOL,NDIEL,IGC A2=(SR/RCUT)**6 A6=ALPHA/(SR**4) IF(A2.GT..5) GO TO 20 A3=A2 A4=1. A5=A3 10 A4=A4+1. A3=-A3*A2/A4 A5=A5+A3 IF(ABS(A3).GT.1.E-7) GO TO 10 A6=A6*A5 GO TO 30 20 IF(A2.GT.16.) GO TO 30 A6=(1.-EXP(-A2))*A6 30 IF(NDIEL.EQ.1) A6=SQRT(A6) FUNSF=A6 RETURN END C*********************************************************************** SUBROUTINE INMAT(A,NDIMA,NRR,NCR,NTP) C INCLUDE 'PARAM' DIMENSION A(NDIMA,*) DO 10 I=1,NCR 10 READ(NTP) (A(J,I),J=1,NRR) RETURN END C*********************************************************************** SUBROUTINE INV(N) C INCLUDE 'PARAM' C C MATRIX INVERSION ROUTINE WITH FULL PIVOTING C CALLED BY SUBROUTINES FG AND OMEGA C COMMON /WORK2/A(MZ07,MZ07),CC(MZ07,MZ07),B(MZ07,MZ07) X ,CS(MZ07,MZ07) DIMENSION IPIVOT(MZ07),IND1(MZ07),IND2(MZ07) C IF(N.EQ.1) GO TO 741 C INITIALISATION DO 20 J=1,N 20 IPIVOT(J)=0 DO 550 I=1,N C SEARCH FOR PIVOT C AMAX=0. DO 105 J=1,N IF(IPIVOT(J)-1) 60,105,60 60 DO 100 K=1,N IF(IPIVOT(K)-1) 80,100,740 80 IF(ABS(AMAX)-ABS(A(J,K))) 85,100,100 85 IR=J IC=K AMAX =A(J,K) 100 CONTINUE 105 CONTINUE IPIVOT(IC)=IPIVOT(IC)+1 C IF(IR-IC) 150,260,150 150 DO 200 L=1,N SWAP=A(IR,L) A(IR,L)=A(IC,L) 200 A(IC,L)=SWAP 260 IND1(I)=IR IND2(I)=IC PIVI=1./A(IC,IC) A(IC,IC)=1. DO 350 L=1,N 350 A(IC,L)=A(IC,L)*PIVI DO 550 L1=1,N IF(L1-IC) 400,550,400 400 T=A(L1,IC) A(L1,IC)=0. DO 450 L=1,N 450 A(L1,L)=A(L1,L)-A(IC,L)*T 550 CONTINUE L=N+1 DO 710 I=1,N L=L-1 IF(IND1(L)-IND2(L)) 630,710,630 630 JR=IND1(L) JC=IND2(L) DO 705 K=1,N SWAP=A(K,JR) A(K,JR)=A(K,JC) 705 A(K,JC)=SWAP 710 CONTINUE 740 RETURN C 741 A(1,1)=1./A(1,1) RETURN END C*********************************************************************** SUBROUTINE ITERA(NIT,EROR,RO,RF,IBUG,IOUT) C INCLUDE 'PARAM' C C PHASE BY WBK AND AMPLITUDE BY ITERATION OF DIFFERENTIAL EQUATION C NOT APPLICABLE FOR NZD = 0 C NIT=MAXIMUM NUMBER OF ITERATIONS TO BE ATTEMPTED C EROR=CONVERGENCE CRITERION C RO=DESIRED VALUE OF R, USED FOR LAST MESH POINT C RF=MINIMUM VALUE OF R FOR WHICH CONVERGENCE OBTAINED C C LOGICAL CNVG COMMON NCHD1,IJNZ(MZ25),KNZ(5),BLIN(MZ25),NZD,NZDM, C AB1(MZ07,13,10),AB2(MZ07,13,10),ALP(13,13),CC(MZ07,MZ07,5), C E2(MZ07),PP(12,10),T1(MZ07),V(MZ11,11,10),ZET(12,10),ZTIN(12,10) C ,GAM(2,15,MZ07) C , VF1(MZ07),VF2(MZ07),VFDR1(MZ07),VFDR2(MZ07),ZA,ZB DIMENSION DUMTM2(14,MZ07),DUMTM1(14,MZ07) DIMENSION TM1N(MZ07),TM2N(MZ07),TM1(14*MZ07),TM2(14*MZ07) EQUIVALENCE (DUMTM1(1,1),TM1(1)),(DUMTM2(1,1),TM2(1)) COMMON /CMN02/ NCNL,NOP,LAMAX COMMON /CMN04/ HY,NHY,J C C MESH FOR INTEGRATION OVER VARIABLE U=SQRT(1./R) IN EQUATION (24) C OF NORCROSS, CPC 1 (1969) 88. INTEGRATIONS FROM U=0. TO C U=SQRT(1./R0) ARE PERFORMED WITH NHY EQUAL INTERVALS IN U, C WHERE NHY IS SET EQUAL TO 10 IN ASYSM. IF CONVERGENCE IS NOT C SATISFACTORY AT R = R0, ITERA USES A LARGER VALUE , RF , OF R C AND CORRESPONDING SMALLER NUMBER , IMAX , OF MESH INTERVALS. C THE SMALLEST VALUE OF IMAX CONSIDERED TO BE SATISFACTORY IS C (NGY + 1) WHERE NGY IS SET EQUAL TO NHY/3 . C C INITIALIZATION C NGY=NHY/3 ZA1= NZD RF=RO NHY2=NHY-NGY NHY1=NHY-1 NITM1=NIT-1 NDRV=NIT+NIT NDAB=NDRV+1 NPP=NDRV-1 C CALL ZETA(NDRV) CALL PPFS(NPP) C C DEFINE CHAIN RULE COEFFICIENTS FOR USE IN SUBROUTINE DCHAIN C ALP(1,1)=1.0 ALP(2,1)=1.0 ALP(1,2)=1. ALP(2,2)=2. ALP(3,2)=1. K=2 KM=1 DO 50 KK=4,NPP K=K+1 KM=KM+1 ALP(1,K)=1. ALP(2,K)=K ALP(KK,K)=1. DO 50 KP=3,K 50 ALP(KP,K)=ALP(KP,KM)+ALP(KP-1,KM) C C FIRST APPROXIMATION C DO 81 I=1,NCNL DO 81 ND=1,NDAB DO 81 IM=1,NHY AB1(I,ND,IM)=0. 81 AB2(I,ND,IM)=0. DO 82 IM=1,NHY 82 AB1(J,1,IM)=1.0 CNVG=.TRUE. C C START ITERATING C DO 120 NT=2,NIT NDRV=NDRV-2 NDG=NDRV+1 CALL POP(1,1,2) CALL QROP(1,1,1,2) DO 104 NDR=2,NDG ND=NDR CALL POP(ND,1,2) CALL QROP(ND,1,1,2) ND=ND-1 CALL POP(ND,2,1) CALL QROP(ND,2,2,1) 104 CONTINUE IF(NT.EQ.2) GO TO 110 C C CONVERGENCE CHECK C DO 105 I=1,NCNL IF(ABS(AB1(I,1,NHY)-TM1N(I)).GT.EROR) GO TO 153 IF(ABS(AB2(I,1,NHY)-TM2N(I)).GT.EROR) GO TO 153 105 CONTINUE GO TO 130 153 IF(NT.EQ.NIT) GO TO 129 110 DO 115 I=1,NCNL TM2N(I)=AB2(I,1,NHY) 115 TM1N(I)=AB1(I,1,NHY) IF(NT.NE.NITM1) GO TO 120 ICT=0 DO 182 IM=NGY,NHY1 DO 182 I=1,NCNL ICT=ICT+1 TM2(ICT)=AB2(I,1,IM) 182 TM1(ICT)=AB1(I,1,IM) DO 183 I=1,NCNL ICT=ICT+1 TM1(ICT)=TM1N(I) 183 TM2(ICT)=TM2N(I) 120 CONTINUE C 129 CNVG=.FALSE. 130 IMAX=NHY C C CONVERGENCE FOR RO C IF(CNVG) GO TO 140 C C SEARCH FOR MINIMUM VALUE OF R FOR WHICH CONVERGENCE HOLDS C ICT=0 DO 107 IM=NGY,NHY IN=IM DO 107 I=1,NCNL ICT=ICT+1 IF(ABS(AB2(I,1,IM)-TM2(ICT)).LE.EROR) GO TO 1029 GO TO 1030 1029 IF(ABS(AB1(I,1,IM)-TM1(ICT)).LE.EROR) GO TO 107 1030 IMAX=IN-1 GO TO 119 107 CONTINUE 119 ZA1=FLOAT(IMAX)*HY RF=1./(ZA1*ZA1*FLOAT(NZD)) IF(IMAX.GT.NGY) GO TO 200 C C IMAX TOO SMALL WRITE(IOUT, 601) IMAX,RF RF=1.E+8 RETURN C C CONVEGENCE OBTAINED C DEFINITION OF AMPLITUDE PARAMETERS 200 CONTINUE 140 DO 125 I=1,NCNL VF1(I)=AB1(I,1,IMAX) VFDR1(I)=AB1(I,2,IMAX) VF2(I)=AB2(I,1,IMAX) 125 VFDR2(I)=AB2(I,2,IMAX) ZA=ZET(1,IMAX) ZB=ZET(2,IMAX)*ZTIN(1,IMAX)*0.5 IF(IBUG.EQ.1) WRITE(IOUT,600) J,IMAX,RF RETURN C 600 FORMAT(34H TRIED ITERATION METHOD IN CHANNEL,I3,27H, CONVERGENCE A CT MESH POINT,I3,7H (RFI =,E13.5,1H)) 601 FORMAT(' ITERATIVE METHOD FAILED. OBTAINED IMAX =',I5, C ' WHICH IS TOO SMALL.'/' IF ITERATIVE METHOD REQUIRED ', C 'USE VALUE OF MFG SUCH THAT ASYMPTOTIC ROUTINES CALLED FOR ', C 'R GREATER THAN ',E13.5) C END C*********************************************************************** SUBROUTINE LAP C INCLUDE 'PARAM' C C COMPUTES MATRIX OPERATORS FOR OVERLAP INTEGRALS WITH P(IT,IG) FOR C (IPQ.NE.-1), AND WITH Q(IT,IG) FOR (IPQ.EQ.-1) C THE METHOD USED IS DESCRIBED IN CPC APPENDIX B(III,A) C INTEGER HEAD COMMON /HOLD/DEL1,DEL2,DELNU,ELOW,EMAX,FINT,FRAT,HEAD,IDM,SZ,Z, C BLAM(MZ07,MZ07,5),EKLSQ(MZ01),ET(MZ07),ETERM(MZ02), C F1(2,MZ07,MZ07),F2(2,MZ07,MZ07),FI(MZ05,MZ05),FLAG(MZ09,3,4), C FX(2,MZ07,MZ07),P(MZ06,MZ03),Q(MZ06,MZ03),R(MZ06),S(MZ10,3), C SD1(MZ05,MZ05),SD2(MZ05,MZ05),SH(4,5),T(MZ10,3), C IFILE,IPRINT,ISTOT,ITC,ITCP1,ITMX,JFILE,KCASE,LAMMX,LCUT,LTOT, C MCH,MDG,MODE,N,N1,N2,NCHB,NCHCL,NCHF,NCHFMX,NCHOP,NLCTRN,NLOW, C NMX,NRAD,NTERM,NTHROW,NTOT,NTOTMX, C ISTERM(MZ02),ITERM(MZ07),LL(MZ07),LQN(MZ03),LTERM(MZ02), C MDEGEN(MZ01),MTCHCL(MZ01),NQN(MZ03) COMMON /WORK1/PP(MZ06),PP1(MZ06),PP2(MZ06),RHO(5) COMMON /OPS/Y(MZ06,MZ06),FLAP(MZ06) COMMON /OVSUB2/ IG,IGP,IGPP,IG1,L,IPQ,MOAD,LAM C L1=LQN(IG) L=L1 RITC=R(ITC) M1=L+L1+1 C C INITIALISE FLAP TO ZERO DO 1 IT=ITCP1,N 1 FLAP(IT)=0. C C FOR IPQ.NE.-1 , STORE P(J,IG) IN PP(J) IF(IPQ.EQ.-1) GO TO 110 DO 100 J=1,N 100 PP(J)=P(J,IG) GO TO 130 C FOR IPQ.EQ.-1 , STORE Q(J,IG) IN PP(J) 110 DO 120 J=1,N 120 PP(J)=Q(J,IG) C C COMPUTE COEFFICIENTS FOR EXPANSION OF PP AND STORE IN PP1 130 DO 2 I=1,ITC A1=0. DO 3 J=1,ITC 3 A1=A1+FI(I,J)*PP(J) 2 PP1(I)=A1 C C INITIALISE FLAP(K),K=1,ITC A2=RITC**M1 DO 4 K=1,ITC A1=0. A3=A2 M2=M1 DO 5 J=1,ITC B1=FI(J,K) A3=A3*RITC B3=A3 M2=M2+1 B2=M2 DO 5 I=1,ITC B2=B2+1. B3=B3*RITC 5 A1=A1+B1*B3*PP1(I)/B2 4 FLAP(K)=A1 C C REPLACE PP(I) BY PP(I)*R(I)**(L+1) FOR I=ITC-1, ITC I=ITC-1 PP(I)=PP(I)*R(I)**(L+1) PP(ITC)=PP(ITC)*RITC**(L+1) C C COMPUTE AND STORE FLAP(IS),IS=1,N RIT=RITC IT3=ITC-3 DO 8 IT=ITCP1,N IT3=IT3+1 JMAX=4 IF (IT.EQ.N) JMAX=3 ITMC=IT-ITC RITM1=RIT RIT=R(IT) DR=0.25*(RIT-RITM1) C A1=SH(4,1)*DR*PP(IT-1) IF (IT.EQ.ITCP1) A1=A1*RITM1**(L+1) FLAP(IT-1)=FLAP(IT-1)+A1 C A1=SH(4,5)*DR*PP(IT) FLAP(IT)=FLAP(IT)+A1 C DO 8 MU=2,4 A1=0. MU1=MU-1 DO 11 J=1,JMAX 11 A1=A1+FLAG(ITMC,MU1,J)*PP(IT3+J) A1=A1*SH(4,MU)*DR IS=IT3 DO 8 J=1,JMAX IS=IS+1 B1=FLAG(ITMC,MU1,J) IF (IS.LE.ITC) B1=B1*R(IS)**(L+1) 8 FLAP(IS)=FLAP(IS)+A1*B1 C RETURN END C*********************************************************************** SUBROUTINE MATS C INCLUDE 'PARAM' C DESCRIBED IN CPC SECTION 5.2 AND APPENDIX B C C INTEGER HEAD C C COMMON /HOLD/DEL1,DEL2,DELNU,ELOW,EMAX,FINT,FRAT,HEAD,IDM,SZ,Z, C BLAM(MZ07,MZ07,5),EKLSQ(MZ01),ET(MZ07),ETERM(MZ02), C F1(2,MZ07,MZ07),F2(2,MZ07,MZ07),FI(MZ05,MZ05),FLAG(MZ09,3,4), C FX(2,MZ07,MZ07),P(MZ06,MZ03),Q(MZ06,MZ03),R(MZ06),S(MZ10,3), C SD1(MZ05,MZ05),SD2(MZ05,MZ05),SH(4,5),T(MZ10,3), C IFILE,IPRINT,ISTOT,ITC,ITCP1,ITMX,JFILE,KCASE,LAMMX,LCUT,LTOT, C MCH,MDG,MODE,N,N1,N2,NCHB,NCHCL,NCHF,NCHFMX,NCHOP,NLCTRN,NLOW, C NMX,NRAD,NTERM,NTHROW,NTOT,NTOTMX, C ISTERM(MZ02),ITERM(MZ07),LL(MZ07),LQN(MZ03),LTERM(MZ02), C MDEGEN(MZ01),MTCHCL(MZ01),NQN(MZ03) COMMON C(MZ06),RHO(5),D(4) Z2=Z+Z C C COMPUTE AND STORE FLAG ITMC=0 IT3=ITC-3 DO 1 IT=ITCP1,N ITMC=ITMC+1 IT3=IT3+1 JMAX=4 IF (IT.EQ.N) JMAX=3 I=IT3 DO 2 J=1,JMAX I=I+1 2 D(J)=R(I) A=0.25*(D(3)-D(2)) RHO(1)=D(2)+A RHO(2)=RHO(1)+A RHO(3)=RHO(2)+A DO 1 I=1,JMAX A=1. DO 4 J=1,JMAX IF (J.EQ.I) GO TO 4 A=A*(D(I)-D(J)) 4 CONTINUE A=1./A DO 1 IS=1,3 A1=A DO 6 J=1,JMAX IF (J.EQ.I) GO TO 6 A1=A1*(RHO(IS)-D(J)) 6 CONTINUE 1 FLAG(ITMC,IS,I)=A1 C C COMPUTE AND STORE FI MS =ITC-1 DO 11 IT=1,ITC A1=1. R0=R(IT) J=0 DO 7 I=1,ITC IF (I.EQ.IT) GO TO 7 A1=A1*(R0-R(I)) J=J+1 C(J)=R(I) 7 CONTINUE A1=1./A1 FI(2,IT)=A1 FI(1,IT)=-A1*C(1) DO 9 NS=2,MS NS2=NS+1 FI(NS2,IT)=FI(NS,IT) A2=C(NS) DO 8 NS1=2,NS NS2=NS2-1 8 FI(NS2,IT)=FI(NS2-1,IT)-A2*FI(NS2,IT) 9 FI(1,IT)=-A2*FI(1,IT) 11 CONTINUE C C COMPUTE AND STORE SD1 AND SD2 DO 13 IT=1,ITC RIT=R(IT) DO 14 IS=1,ITC A1=0. A2=0. A4=1. DO 15 K=2,ITC RK=K A4=A4*RIT A3=FI(K,IS)*(RK-1.)*A4 A1=A1+A3 15 A2=A2+(RK+4.)*A3 SD1(IT,IS)=2.*A1 14 SD2(IT,IS)=A2 SD1(IT,IT)=SD1(IT,IT)+4. 13 SD2(IT,IT)=SD2(IT,IT)+6.+Z2*RIT C C COMPUTE AND STORE S AND T R0=R(ITC) R1=R(ITC-1) A=R0-R1 A2=A*A DO 10 I=ITCP1,N2 A1=A A3=A2 R2=R1 R1=R0 R0=R(I) A=R0-R1 A2=A*A A4=A+A1 RK=A*(A1*A4-A2) T(I,1)=RK S(I,1)=12.*A+Z2*RK/R2 RK=A4*(3.*A*A1+A3+A2) T(I,2)=RK S(I,2)=-12.*A4+Z2*RK/R1 RK=A1*(A*A4-A3) T(I,3)=RK 10 S(I,3)=12.*A1+Z2*RK/R0 C RETURN END C********************************************************************** SUBROUTINE OMEGA C INCLUDE 'PARAM' C DESCRIBED IN CPC SECTION 8.3 C C CALCULATES COLLISION STRENGTHS OMEGA AND OTHER OUTPUT FOR THE C COLLISION PROBLEM C C INTEGER HEAD COMMON C(MZ08,MZ08),H(MZ08,MZ07) COMMON /HOLD/DEL1,DEL2,DELNU,ELOW,EMAX,FINT,FRAT,HEAD,IDM,SZ,Z, C BLAM(MZ07,MZ07,5),EKLSQ(MZ01),ET(MZ07),ETERM(MZ02), C F1(2,MZ07,MZ07),F2(2,MZ07,MZ07),FI(MZ05,MZ05),FLAG(MZ09,3,4), C FX(2,MZ07,MZ07),P(MZ06,MZ03),Q(MZ06,MZ03),R(MZ06),S(MZ10,3), C SD1(MZ05,MZ05),SD2(MZ05,MZ05),SH(4,5),T(MZ10,3), C IFILE,IPRINT,ISTOT,ITC,ITCP1,ITMX,JFILE,KCASE,LAMMX,LCUT,LTOT, C MCH,MDG,MODE,N,N1,N2,NCHB,NCHCL,NCHF,NCHFMX,NCHOP,NLCTRN,NLOW, C NMX,NRAD,NTERM,NTHROW,NTOT,NTOTMX, C ISTERM(MZ02),ITERM(MZ07),LL(MZ07),LQN(MZ03),LTERM(MZ02), C MDEGEN(MZ01),MTCHCL(MZ01),NQN(MZ03) COMMON /FGCOM/R1,R2,KJ(MZ07),METHOD,MFG COMMON /FORM/JFORM COMMON /ORDER/IROW(MZ08),ICOL(MZ12) COMMON /OUT/IPUNCH,ICONV COMMON /SCRATCH/ KFILE COMMON /WORK2/CC(MZ07,MZ07),A(MZ07,MZ07),B(MZ07,MZ07) X ,CS(MZ07,MZ07) C DIMENSION F(MZ08),FAS(MZ22,MZ07,MZ07),FBAR(MZ08),G(MZ12,MZ07), C GQ(MZ12),RFG(MZ07,MZ07),RFGST(MZ22,MZ07,MZ07) DIMENSION AINT(MZ07,MZ07),BINT(MZ07,MZ07),CINT(MZ07,MZ07) C EQUIVALENCE (C(1,1),FAS(1,1,1),G(1,1),RFGST(1,1,1)),(FX(1,1,1) C,RFG(1,1)) C EQUIVALENCE (FBAR(1),IROW(1)),(GQ(1),ICOL(1)) C C SQK=EKLSQ(KCASE) C C PUT SOLUTION MATRIX IN ARRAY G C PUT CAPITAL A FOR OPEN CHANNELS IN RFG DO 901 I=1,NCHOP DO 901 J=1,NCHOP 901 RFG(I,J)=C(I,J) C PUT H IN FIRST NTOT ROWS OF G DO 902 I=1,NTOT DO 902 J=1,NCHOP 902 G(I,J)=H(I,J) C PUT RFG IN LAST NCHOP ROWS OF G I=NTOT DO 903 K=1,NCHOP I=I+1 DO 903 J=1,NCHOP 903 G(I,J)=RFG(K,J) C C STORES MATRICES A AND B SUCH THAT (FOR ALL J=1,NCHOP) C A(I,J)= LOWER CASE A FOR I=1,NCHOP C B(I,J)= LOWER CASE B FOR I=1,NCHOP C A(I,J)= LOWER CASE D FOR I.GT.NCHOP C B(I,J)= CAPITAL A FOR I.GT.NCHOP LI=0 DO 15 I=1,NCHF LI=LI+N1 DO 12 IP=1,NCHOP 12 A(I,IP)=G(LI,IP) LI=LI+1 DO 14 IP=1,NCHOP 14 B(I,IP)=G(LI,IP) 15 CONTINUE C C COMPUTE RADIAL FUNCTIONS AT POINTS (N+1) AND (N+2) J1=-1 J2=0 DO 16 I=1,NCHF J1=J1+N2 J2=J2+N2 DO 16 IP=1,NCHOP A1=0. A2=0. DO 17 K=1,NCHF A3=A(K,IP) A1=A1+F1(1,I,K)*A3 17 A2=A2+F2(1,I,K)*A3 DO 18 K=1,NCHOP A3=B(K,IP) A1=A1+F1(2,I,K)*A3 18 A2=A2+F2(2,I,K)*A3 G(J1,IP)=A1 16 G(J2,IP)=A2 C C PUT ARRAY CAPITAL A (DETERMINING BEHAVIOUR AT THE ORIGIN) C IN LOCATION H LI=NCHOP+1 DO 21 J=1,NCHOP DO 22 I=1,NCHOP 22 H(I,J)=G(NTOT+I,J) IF(NCHCL.EQ.0) GO TO 21 DO 23 I=LI,NCHF 23 H(I,J)=B(I,J) 21 CONTINUE C C PRINT THE RADIAL FUNCTIONS IF (IPRINT.NE.4.AND.IPRINT.NE.5) GO TO 25 M1=(NCHF-1)/10+1 DO 20 IP=1,NCHOP IF(NTHROW.EQ.0) WRITE(6,600) IP IF(NTHROW.EQ.1) WRITE(6,1600) IP M2=1 DO 30 I1=1,M1 M3=M2+9 IF(M3.GT.NCHF) M3=NCHF WRITE(6,601) (L,L=M2,M3) WRITE(6,602) LI=(M3-1)*N2 M4=(M2-1)*N2 DO 27 IT=1,N2 M4=M4+1 LI=LI+1 WRITE(6,603) IT,(G(LIP,IP),LIP=M4,LI,N2) 27 IF(IT.EQ.ITC.OR.IT.EQ.N) WRITE(6,602) 30 M2=M2+10 WRITE(6,624)(H(I,IP),I=1,NCHF) IF (NCHB.EQ.0) GO TO 40 WRITE(6,604) DO 39 J=1,NCHB I=LI+J 39 WRITE(6,612)J,G(I,IP) 40 LI=LI+NCHB LIP=NTOT-LI IF (LIP.EQ.0) GO TO 41 WRITE(6,605) DO 19 J=1,LIP I=LI+J 19 WRITE(6,612)J,G(I,IP) 41 CONTINUE C C PRINT THE MATRICES SMALL A AND SMALL B (LOWER CASE A AND B) WRITE(6,625)(A(I,IP),I=1,NCHF) WRITE(6,626)(B(I,IP),I=1,NCHOP) 20 CONTINUE 25 IF(IPRINT.NE.3) GO TO 60 IF(NTHROW.EQ.0) WRITE(6,606) IF(NTHROW.EQ.1) WRITE(6,1606) DO 70 I=1,NCHOP 70 WRITE(6,607)I,(A(I,IP),IP=1,NCHOP) WRITE(6,608) DO 80 I=1,NCHOP 80 WRITE(6,607)I,(B(I,IP),IP=1,NCHOP) C C PUNCH OUT SOME QUANTITIES C WRITE(7,1003) NCHF,NCHOP WRITE(7,1001) (ET(I),I=1,NCHF) WRITE(7,1003) (LL(I),I=1,NCHF) DO 800 I=1,NCHOP DO 790 JJJ=1,NCHOP BINT(I,JJJ)=B(I,JJJ) CINT(I,JJJ)=CC(I,JJJ) CC(I,JJJ)=A(I,JJJ) 790 CONTINUE 800 CONTINUE C C INVERT THE A MATRIX C CALL INV(NCHOP) C C FORM THE REACTANCE MATRIX C WRITE(6,1000) DO 830 I=1,NCHOP DO 820 JJJ=1,NCHOP AINT(I,JJJ)=0.0 DO 810 KKK=1,NCHOP AINT(I,JJJ)=AINT(I,JJJ)+BINT(I,KKK)*CC(KKK,JJJ) 810 CONTINUE 820 CONTINUE WRITE(6,1001) (AINT(I,JJJ),JJJ=1,NCHOP) WRITE(7,1001) (AINT(I,JJJ),JJJ=1,NCHOP) 830 CONTINUE C C REPLACE THE CONTENTS OF CC DO 850 I=1,NCHOP DO 840 JJJ=1,NCHOP CC(I,JJJ)=CINT(I,JJJ) 840 CONTINUE 850 CONTINUE C C FORM THE ANALYTIC R MATRIX C IF(SZ.NE.0.0) GO TO 930 DO 880 I=1,NCHOP AK=SQRT(ET(I)) AAK=SQRT(AK) KL=2*LL(I)+1 PROD=1.0 DO 860 KKK=1,KL PROD=PROD*AAK 860 CONTINUE PROD=1.0/PROD DO 870 JJJ=1,NCHOP AINT(I,JJJ)=AINT(I,JJJ)*PROD 870 CONTINUE 880 CONTINUE DO 910 JJJ=1,NCHOP AK=SQRT(ET(JJJ)) AAK=SQRT(AK) KL=2*LL(JJJ)+1 PROD=1.0 DO 890 KKK=1,KL PROD=PROD*AAK 890 CONTINUE PROD=1.0/PROD DO 900 I=1,NCHOP AINT(I,JJJ)=AINT(I,JJJ)*PROD 900 CONTINUE 910 CONTINUE C C WRITE AND PUNCH OUT THE ANALYTIC R MATRIX C WRITE(6,1002) DO 920 I=1,NCHOP WRITE(6,1001) (AINT(I,JJJ),JJJ=1,NCHOP) WRITE(7,1001) (AINT(I,JJJ),JJJ=1,NCHOP) 920 CONTINUE 930 CONTINUE C C CALCULATE AND PRINT THE R MATRIX 60 IF(NTHROW.EQ.0) WRITE(6,609) ISTOT,LTOT,EKLSQ(KCASE) IF(NTHROW.EQ.1) WRITE(6,1609) ISTOT,LTOT,EKLSQ(KCASE) DO 100 I=1,NCHOP DO 100 J=1,NCHOP 100 CC(I,J) = A(I,J) CALL INV(NCHOP) DO 110 I=1,NCHOP DO 110 J=1,NCHOP A1=0. DO 115 L=1,NCHOP 115 A1 = A1 + B(I,L)*CC(L,J) 110 CS(I,J)=A1 DO 120 I=1,NCHOP 120 WRITE(6,607) I,(CS(I,IP),IP=1,NCHOP) 90 CONTINUE C C CALCULATE AND PRINT THE COLLISION STRENGTHS DO 200 I=1,NCHOP DO 200 J=1,NCHOP A1=0. DO 210 K=1,NCHOP 210 A1=A1+A(K,I)*A(K,J)+B(K,I)*B(K,J) 200 CC(I,J) = A1 CALL INV(NCHOP) A4=(2*LTOT+1)*IABS(ISTOT) DO 220 I=1,NCHOP DO 220 J=I,NCHOP A1=0. A2=0. DO 230 K=1,NCHOP A0=A(J,K) B0=B(J,K) DO 230 L=1,NCHOP A3 = CC(L,K) A1=A1+A3*(A(I,L)*A0-B(I,L)*B0) 230 A2=A2+A3*(A(I,L)*B0+B(I,L)*A0) IF (I.EQ.J) A1=A1-1. A3 = A4*(A1*A1+A2*A2)*.5 CS(I,J) = A3 220 CS(J,I) = A3 IF(NTHROW.EQ.0) WRITE(6,610) ISTOT,LTOT,EKLSQ(KCASE) IF(NTHROW.EQ.1) WRITE(6,1610) ISTOT,LTOT,EKLSQ(KCASE) DO 130 I=1,NCHOP DO 130 IP=I,NCHOP 130 WRITE(6,611) I,ITERM(I),LL(I),IP,ITERM(IP),LL(IP),CS(I,IP) C C PUNCHING ON JFILE IF(JFILE.EQ.0) RETURN WRITE(JFILE,730)KCASE,NCHOP,SQK C PUNCH SMALL A AND SMALL B DO 300 I=1,NCHOP DO 300 IP=1,NCHOP IF(JFORM.NE.0) GO TO 290 WRITE(JFILE,750) ISTOT,LTOT,I,IP,A(I,IP),B(I,IP),SQK,HEAD GO TO 300 290 WRITE(JFILE,760) ISTOT,LTOT,I,IP,A(I,IP),B(I,IP),SQK,HEAD 300 CONTINUE WRITE(6,630) C PUNCHING OF FUNCTIONS IF(IPUNCH.EQ.0) RETURN C - STORE FUNCTIONS AT R(N+1) IN F2(K,I,J),K=1,2 AND I=1,NCHF C - AND J=1,NCHOP WHERE C -- FOR K=1, F2=(LINEAR COMBINATION OF SIN SOLUTIONS AND C DECAYING SOLUTIONS) C -- FOR K=2, F2=(LINEAR COMBINATION OF COS SOLUTIONS) DO 405 I=1,NCHF DO 405 J=1,NCHOP A1=0. DO 400 L=1,NCHF 400 A1=A1+F1(1,I,L)*A(L,J) 405 F2(1,I,J)=A1 DO 415 I=1,NCHF DO 415 J=1,NCHOP A1=0. DO 410 L=1,NCHOP 410 A1=A1+F1(2,I,L)*B(L,J) 415 F2(2,I,J)=A1 C - RETURN TO ORIGINAL NORMALISATION IN STARTING REGION C -- RE-NORMALIZE CAP A AND STORE IN A DO 305 I=1,NCHF DO 305 J=1,NCHOP 305 A(I,J)=H(I,J) C STORE FUNCTIONS IN H L=0 DO 320 I=1,NCHF DO 310 IT=1,ITC L=L+1 DO 310 J=1,NCHOP 310 H(L,J)=G(L,J) DO 320 IT=ITCP1,N2 L=L+1 DO 320 J=1,NCHOP 320 H(L,J)=G(L,J) IF(L.EQ.NTOT) GO TO 324 L=L+1 DO 321 J=1,NCHOP DO 321 K=L,NTOT 321 H(K,J)=G(K,J) C - PUNCH CAP A 324 DO 325 J=1,NCHOP 325 WRITE(JFILE,700)(A(I,J),I=1,NCHF) C - PUNCH SMALL G FOR IT.LE.ITC,AND F FOR IT.GT.ITC DO 335 J=1,NCHOP L=1 DO 330 I=1,NCHF M=L+N2-1 WRITE(JFILE,700)(H(K,J),K=L,M) 330 L=L+N2 335 CONTINUE C - PUNCH SMALL C IF(NCHB.EQ.0) GO TO 345 L=N2*NCHF M=L+NCHB L=L+1 DO 340 J=1,NCHOP 340 WRITE(JFILE,700)(H(K,J),K=L,M) C COMPUTE AND PUNCH Q 345 REWIND 11 CNOAA READ(KFILE) C CALL INMAT(C,NTOTMX,NTOT,NTOT,KFILE) NCN2=NCHF*N2 C ... BEGIN LOOP 380 ON LINEARLY INDEPENDENT SOLUTIONS ... DO 380 J=1,NCHOP C COMPUTE F=(RADIAL FUNCTION)*ET C COMPUTE FBAR AND STORE IN FBAR L=1 DO 360 I=1,NCHF A1=A(I,J) A3=ET(I) F(L)=A1*A3 FBAR(L)=A1 K=LL(I)+1 A4=Z/FLOAT(K) DO 350 IT=2,ITC L=L+1 A2=R(IT) A5=A1*(1.-A2*A4) +A2*A2*H(L,J) FBAR(L)=A5 350 F(L)=A5*(A2**K)*A3 DO 355 IT=ITCP1,N2 L=L+1 A5=H( L,J) F(L)=A5*A3 355 FBAR(L)=A5 360 L=L+1 IF(L.GT.NTOT) GO TO 364 DO 361 I=L,NTOT 361 FBAR(I)=H(I,J) C COMPUTE Q AND STORE IN GQ 364 DO 370 L=1,NCN2 A4=-F(L) DO 365 IP=1,NTOT 365 A4=A4+C(L,IP)*FBAR(IP) 370 GQ(L)=-A4 C - PUNCH Q L=1 DO 375 I=1,NCHF M=L+N2-1 WRITE(JFILE,700)(GQ(K),K=L,M) 375 L=L+N2 C ... END LOOP 380 ... 380 CONTINUE C FUNCTIONS IN ASYMPTOTIC REGION COMPUTED AT R=R(N+1)+M*A1 C WHERE M=1,(MFG+1) AND A1=(R(N+2)-R(N+1))*0.5 . C FUNCTIONS STORED IN FAS(M,I,J),I=1,NCHF AND J=1,NCHOP C - DEFINE MFGP1=MFG+1 WRITE(JFILE,710) MFGP1 C - READ RFG FOR SIN SOLUTIONS AND STORE IN RFGST DO 420 M=1,MFGP1 READ(11) RFG DO 420 I=1,NCHF DO 420 J=1,NCHF 420 RFGST(M,I,J)=RFG(I,J) C ... BEGIN LOOP 450 ON M ... DO 450 M=1,MFGP1 C - PUT OLD SOLUTIONS IN F1 DO 425 K=1,2 DO 425 I=1,NCHF DO 425 J=1,NCHOP 425 F1(K,I,J)=F2(K,I,J) C - COMPUTE NEW SIN SOLUTIONS DO 435 I=1,NCHF DO 435 J=1,NCHOP A1=0. DO 430 L=1,NCHF 430 A1=A1+RFGST(M,I,L)*F1(1,L,J) 435 F2(1,I,J)=A1 C - READ RFG FOR COS SOLUTIONS READ (11) RFG C - COMPUTE NEW COS SOLUTIONS C - AND LINEAR COMBINATIONS OF NEW SIN AND COS SOLUTIONS DO 445 I=1,NCHF DO 445 J=1,NCHOP A1=0. DO 440 L=1,NCHF 440 A1=A1+RFG(I,L)*F1(2,L,J) F2(2,I,J)=A1 445 FAS(M,I,J)=F2(1,I,J)+A1 C ... END LOOP 450 ... 450 CONTINUE C PUNCH ASYMPTOTIC FUNCTIONS DO 460 J=1,NCHOP DO 460 I=1,NCHF 460 WRITE(JFILE,700)(FAS(M,I,J),M=1,MFGP1) WRITE(6,640) RETURN C 600 FORMAT(1H1,17HSUBROUTINE OMEGA /,1X,17(1H*)//1X, C38HSOLUTIONS OF THE EQUATIONS - SOLUTION ,I2) 1600 FORMAT(/////,18H SUBROUTINE OMEGA /,2X,17(1H*)//,1X, C38HSOLUTIONS OF THE EQUATIONS - SOLUTION ,I2) 601 FORMAT(//,20X,9H CHANNELS/(9X,9(I2,11X),I2)) 602 FORMAT (1H0) 603 FORMAT(1H ,I3,E12.4,9E13.4) 604 FORMAT(/14H C PARAMETERS/) 605 FORMAT(/16H LAM PARAMETERS/) 1606 FORMAT(/////,29H MATRICES SMALL A AND SMALL B/1H ,28(1H*)// C 4H A ,/,4H *** ) 606 FORMAT(1H1,29H MATRICES SMALL A AND SMALL B/1H ,28(1H*)//4H A /4H C ***) 607 FORMAT(1H0,I2,E12.4,9E13.4/(3X,E12.4,9E13.4)) 608 FORMAT (/4H B /4H ***) 609 FORMAT(1H1,24HR MATRIX FOR SLPI = ,2I3,9H EKLSQ =,F12.6/10H * Z********//) 1609 FORMAT(/////,25H R MATRIX FOR SLPI = ,2I3,9H EKLSQ = ,F12.6,/ C,1H ,9(1H*)) 610 FORMAT(1H1,34H COLLISION STRENGTHS FOR SLPI = ,2I3,9H EKLSQ =, YF12.6/1H ,19(1H*)//1H0,60H CHANNEL TERM LL CHANNELP TER ZMP LLP OMEGA ) 1610 FORMAT(/////,35H COLLISION STRENGTHS FOR SLPI = ,2I3,9H EKLSQ C= ,F12.6,/,1H ,19(1H*)//1H0,60H CHANNEL TERM LL CHANNELP C TERM LLP OMEGA ) 611 FORMAT (1H0,2X,I3,6X,I3,3X,I3,9X,I3,5X,I3,3X,I3,6X,E11.4) 612 FORMAT(I4,E13.4) 624 FORMAT(/,12H CAPITAL A //(4X,E12.4,9E13.4)) 625 FORMAT(/39H VALUES OF SMALL A FOR OPEN CHANNELS, , C28H SMALL D FOR CLOSED CHANNELS //(4X,E12.4,9E13.4)) 626 FORMAT(/39H VALUES OF SMALL B FOR OPEN CHANNELS //,(4X,E12.4,9E C13.4)) 630 FORMAT(//40H MATRICES SMALL A AND B WRITTEN ON JFILE) 640 FORMAT(/34H RADIAL FUNCTIONS WRITTEN ON JFILE) C 700 FORMAT(5E14.7) 710 FORMAT(24I3) 730 FORMAT(2I3,E14.6) 750 FORMAT(2I3,16X,2I3,2X,2E12.5,E14.6,A4) 1000 FORMAT(//24H THE REACTANCE MATRIX IS) 1001 FORMAT(5E14.6) 1002 FORMAT(//24H THE ANALYTIC MATRIX IS) 1003 FORMAT(12I5) 760 FORMAT(2I3,16X,2I3,2X,2E14.7,E14.6,A4) END C*********************************************************************** SUBROUTINE OUTMAT(A,NDIMA,NRW,NCW,NTP) C INCLUDE 'PARAM' DIMENSION A(NDIMA,*) DO 10 I=1,NCW 10 WRITE(NTP) (A(J,I),J=1,NRW) RETURN END C*********************************************************************** FUNCTION PHAS(RHO,E2,XKEY) C INCLUDE 'PARAM' C SEE CPC APPENDIX C SECTION (X) C C RHO = SCALED RADIAL VARIABLE = Z * R C E2 = SCALED ENERGY = (K**2)/(Z**2) C XKEY = 0 FOR EXPAN METHOD OF SOLUTION C = 1 FOR ITERA METHOD OF SOLUTION C C ONLY CALLED FOR CASE Z .GT. 0 , E2 .GT. 0 C CALCULATES PH PART OF PHASE, EQUATIONS (C25) AND (C26) C ALF=SQRT(E2) IF(XKEY.EQ.1.0) GO TO 30 C C EXPAN METHOD C ------------ C C USE FORM OF PHAS GIVEN BY EQUATION (C25) RH=ALF*RHO A1=1./ALF PHAS=RH+A1*(1.+LOG(2.*RH*ALF)) RETURN C C ITERA METHOD C ------------ C C PHAS DEFINED BY EQUATION (C26) 30 RX=RHO*SQRT(E2+2./RHO) A1=1.+RHO*E2 A2=ALF*RX/A1 IF(A2.GE.0.1) GO TO 60 C USE FORM OF PHAS GIVEN BY EQUATION (C28) A3=A2*A2 A2=A1+A1 A4=1. A6=1. 40 A5=A4/A1 A6=A6+A5 IF(A5.LT.1.E-6) GO TO 50 A1=A1+A2 A4=A4*A3 GO TO 40 50 PHAS=RX*A6 RETURN C USE FORM OF PHAS GIVEN BY EQUATION (C26) 60 PH=1.+RHO*E2+RX*ALF PHAS=RX+LOG(PH)/ALF RETURN C END C*********************************************************************** SUBROUTINE POP(NDRV,NBA,NAP) C INCLUDE 'PARAM' C C ITERATION OPERATOR FOR DEGENERATE ENERGIES C OPERATOR P ON A(NAP=1) OR B(NAP=2) C RESULT STORED IN A(NBA=1) OR B(NBA=2) C RESULT TIMES MINUS 1 IF NBA=1 C DERIVATIVE NDRV-1 OF P IN STATE I C RESULT OBTAINED FOR EACH MESH POINT C C COMMON NCHD1,IJNZ(MZ25),KNZ(5),BLIN(MZ25),NZD,NZDM, C AB1(MZ07,13,10),AB2(MZ07,13,10),ALP(13,13),CC(MZ07,MZ07,5), C E2(MZ07),PP(12,10),T1(MZ07),V(MZ11,11,10),ZET(12,10),ZTIN(12,10) C ,GAM(2,MZ07,15) COMMON /CMN02/ NCNL,NOP,LAMAX COMMON /CMN03/ NDL(MZ07), NDU(MZ07) COMMON /CMN04/ HY,NHY,J COMMON /CMN14/ CH(10),BI(10),BO(10),P1(10),Z1(10) DIMENSION AB(10,MZ07),NJ(MZ07) NDLJ=NDL(J) NDUJ=NDU(J) SGN=0.5 IF(NBA.EQ.1) SGN=-SGN IF(J.GT.NOP) SGN=-SGN IF(NDRV.GT.1) GO TO 45 IF(NAP.EQ.2) GO TO 94 DO 93 IM=1,NHY Z=ZTIN(1,IM) Z1(IM)=Z P1(IM)=PP(1,IM) DO 93 IK=1,NCNL 93 AB(IM,IK)=AB1(IK,1,IM)*Z GO TO 96 94 DO 95 IM=1,NHY Z=ZTIN(1,IM) Z1(IM)=Z P1(IM)=PP(1,IM) DO 95 IK=1,NCNL 95 AB(IM,IK)=AB2(IK,1,IM)*Z 96 CONTINUE DO 80 I=1,NCNL IF(I.LT.NDLJ) GO TO 80 IF(I.GT.NDUJ) GO TO 80 DO 15 IK=I,NCNL 15 NJ(IK)=I+(IK*IK-IK)/2 IF(I.EQ.1) GO TO 17 I1=I-1 I2=(I*I1)/2 DO 16 IK=1,I1 16 NJ(IK)=IK+I2 17 CONTINUE Y=0. DO 30 IM=1,NHY SUM=0. Y=Y+HY DO 20 IK=1,NCNL IJ=NJ(IK) 20 SUM=SUM+AB(IM,IK)*V(IJ,1,IM) 30 BI(IM)=2.*(SUM+AB(IM,I)*P1(IM))/(Y*Y*Y) BOR=0.0 IF(I.NE.J) GO TO 135 IF(NBA.EQ.1) GO TO 135 IF(E2(J).NE.0.0) GO TO 135 BOR=SQRT(2.0)*(CC(J,J,1)-0.1875) 135 CALL BOOLE(BOR) CON=0. IF(I.NE.J) GO TO 137 IF(NBA.EQ.1) GO TO 136 IF(J.LE.NOP) GO TO 137 136 CON=1. 137 IF(NBA.EQ.2) GO TO 37 IF(NAP.EQ.2) GO TO 35 DO 34 IM=1,NHY 34 AB1(I,1,IM)=SGN*(AB1(I,2,IM)*Z1(IM)-BO(IM))+CON GO TO 80 35 DO 36 IM=1,NHY 36 AB1(I,1,IM)=SGN*(AB2(I,2,IM)*Z1(IM)-BO(IM))+CON GO TO 80 37 IF(NAP.EQ.2) GO TO 39 DO 38 IM=1,NHY 38 AB2(I,1,IM)=SGN*(AB1(I,2,IM)*Z1(IM)-BO(IM))+CON GO TO 80 39 DO 40 IM=1,NHY 40 AB2(I,1,IM)=SGN*(AB2(I,2,IM)*Z1(IM)-BO(IM))+CON 80 CONTINUE RETURN C 45 NP1=NDRV NP2=NDRV-1 DO 180 I=1,NCNL IF(I.LT.NDLJ) GO TO 180 IF(I.GT.NDUJ) GO TO 180 II=I CALL DCHAIN(P1,II,NP1,0,NAP,2,PP,0,ZTIN,1) CALL DCHAIN(Z1,II,NP2,1,NAP,1,ZTIN,1,PP,0) CALL DCHAIN(CH,II,NP2,0,NAP,1,PP,1,ZTIN,1) IF(NBA.EQ.2) GO TO 178 DO 60 IM=1,NHY 60 AB1(I,NDRV,IM)=SGN*(P1(IM)+Z1(IM)+CH(IM)) GO TO 180 178 DO 179 IM=1,NHY 179 AB2(I,NDRV,IM)=SGN*(P1(IM)+Z1(IM)+CH(IM)) 180 CONTINUE RETURN END C ********************************************************************* SUBROUTINE POT C INCLUDE 'PARAM' C C COMPUTES POTENTIAL OPERATORS YPOT(IT) C (SEE CPC SECTION 6.2(II) AND APPENDIX B(III,C) C INTEGER HEAD COMMON /HOLD/DEL1,DEL2,DELNU,ELOW,EMAX,FINT,FRAT,HEAD,IDM,SZ,Z, C BLAM(MZ07,MZ07,5),EKLSQ(MZ01),ET(MZ07),ETERM(MZ02), C F1(2,MZ07,MZ07),F2(2,MZ07,MZ07),FI(MZ05,MZ05),FLAG(MZ09,3,4), C FX(2,MZ07,MZ07),P(MZ06,MZ03),Q(MZ06,MZ03),R(MZ06),S(MZ10,3), C SD1(MZ05,MZ05),SD2(MZ05,MZ05),SH(4,5),T(MZ10,3), C IFILE,IPRINT,ISTOT,ITC,ITCP1,ITMX,JFILE,KCASE,LAMMX,LCUT,LTOT, C MCH,MDG,MODE,N,N1,N2,NCHB,NCHCL,NCHF,NCHFMX,NCHOP,NLCTRN,NLOW, C NMX,NRAD,NTERM,NTHROW,NTOT,NTOTMX, C ISTERM(MZ02),ITERM(MZ07),LL(MZ07),LQN(MZ03),LTERM(MZ02), C MDEGEN(MZ01),MTCHCL(MZ01),NQN(MZ03) COMMON /OPS/Y(MZ06,MZ06),YPOT(MZ06) COMMON /OVSUB2/ IG,IGP,IGPP,IG1,L,IPQ,MOAD,LAM C IF(MOAD.EQ.1) GO TO 5 CALL YLAM(IG1,L,LAM) RETURN C 5 L=LQN(IGP) CALL YLAM(IG,L,LAM) C DO 10 IT=1,N YP=0. DO 11 IS=1,N 11 YP=YP+Y(IT,IS)*P(IS,IGP) 10 YPOT(IT)=YP C C MODIFICATION FOR LAM=0 (CPC EQUATION (6.4)) IF (LAM.NE.0) RETURN L1=-2*(L+1) A1=-YPOT(N)*R(N) IF(IG.EQ.IGP)A1=1.+A1 A2=1./(P(1,IG)*P(1,IGP)) A1=A1/(1.-A2*P(N,IG)*P(N,IGP)*R(N)**L1) DO 12 IT=2,ITC 12 YPOT(IT)=YPOT(IT)+A1*(1.-A2*P(IT,IG)*P(IT,IGP))/R(IT) DO 13 IT=ITCP1,N 13 YPOT(IT)=YPOT(IT)+A1*(1.-A2*P(IT,IG)*P(IT,IGP)*R(IT)**L1)/R(IT) RETURN END C*********************************************************************** SUBROUTINE POTS(NDRV) C INCLUDE 'PARAM' C C POTENTIALS FOR USE IN ITERATION METHOD (SUBROUTINE ITERA) C POTENTIAL AND DERIVATIVES UP TO NDRV-1 FOR ALL MESH POINTS C C COMMON NCHD1,IJNZ(MZ25),KNZ(5),BLIN(MZ25),NZD,NZDM, C AB1(MZ07,13,10),AB2(MZ07,13,10),ALP(13,13),CC(MZ07,MZ07,5), C E2(MZ07),PP(12,10),T1(MZ07),V(MZ11,11,10),ZET(12,10),ZTIN(12,10) C ,GAM(2,MZ07,15) COMMON /CMN02/ NCNL,NOP,LAMAX COMMON /CMN04/ HY,NHY,JJ C IF(NDRV.GT.1) GO TO 4 DO 2 M=1,LAMAX 2 ALP(M,1)=1. GO TO 12 4 DO 10 M=1,LAMAX ALP(M,1)=1. A1=1. FF=M DO 5 ID=2,NDRV FF=FF+1. A1=A1*FF 5 ALP(M,ID)=A1 10 CONTINUE C 12 HY2=HY*HY DO 20 IM=1,NHY Y1=FLOAT(IM*IM)*HY2 YY=-1. DO 20 ID=1,NDRV YY=-Y1*YY DO 20 I=1,NCNL DO 20 J=I,NCNL IJ=I+(J*(J-1))/2 A1=0. YZ=YY DO 15 M=1,LAMAX YZ=Y1*YZ 15 A1=A1+CC(I,J,M)*YZ*ALP(M,ID) 20 V(IJ,ID,IM)=A1 C RETURN END C*********************************************************************** SUBROUTINE PPFS(NDRV) C INCLUDE 'PARAM' C C CORRECTION TERM IN WBK AMPLITUDE DIFFERENTIAL EQUATION C CORRECTION FACTOR Q AND DERIVATIVES UP TO NDRV-1 FOR ALL MESH C POINTS FOR SUBROUTINE ITERA C C COMMON NCHD1,IJNZ(MZ25),KNZ(5),BLIN(MZ25),NZD,NZDM, C AB1(MZ07,13,10),AB2(MZ07,13,10),ALP(13,13),CC(MZ07,MZ07,5), C E2(MZ07),PP(12,10),T1(MZ07),V(MZ11,11,10),ZET(12,10),ZTIN(12,10) C ,GAM(2,MZ07,15) COMMON /CMN04/ HY,NHY,J SGN=1.0 IF(E2(J).LT.0.0) SGN=-1.0 Y1=0. DO 20 IPT=1,NHY Y1=Y1+HY Y2=Y1*Y1 Y3=ZET(1,IPT) Y3=Y2/(Y3*Y3) 20 PP(1,IPT)=(-SGN*Y2*Y2*Y3)*(1.-SGN*1.25*Y3) IF(NDRV.EQ.1) RETURN C S2=SGN+SGN ALP(1,1)=S2+SGN ALP(2,1)=-7. ALP(3,1)=5*SGN DO 360 KK=3,NDRV K=KK-1 KM1=K-1 KU=KK+1 ALP(1,K)=-KU*ALP(1,KM1) ALP(KU,K)=S2*KK*ALP(KK,KM1) DO 360 KP=2,KK KS=KP-1 360 ALP(KP,K)=S2*KS*ALP(KS,KM1)-(KP+KK)*ALP(KP,KM1) C 70 Y0=0. DO 81 IPT=1,NHY Y0=Y0+HY Y1=Y0*Y0 Y2=Y0/ZET(1,IPT) TEMP=Y2*Y2 Y2=Y1*Y1 YMIN=1.0D-76/TEMP K=0 KU=2 DO 81 KK=2,NDRV K=K+1 KU=KU+1 Y2=Y2*Y1 Y3=Y2 S=0. DO 79 KP=1,KU IF(Y3.LT.YMIN) GO TO 80 Y3=TEMP*Y3 79 S=S+ALP(KP,K)*Y3 80 PP(KK,IPT)=S 81 CONTINUE RETURN END C*********************************************************************** SUBROUTINE PRINTB(FL) C INCLUDE 'PARAM' C DESCRIBED IN CPC SECTION 8.3 C C PRINTS THE OUTPUT FROM SUBROUTINE BOUND C C INTEGER HEAD COMMON C(MZ08,MZ08),G(MZ08),GP(MZ08),GPP(MZ08),CNORM(MZ08) COMMON /FGCOM/R1,R2,KJ(MZ07),METHOD,MFG DIMENSION RFG(MZ07,MZ07),RFGST(MZ22,MZ07,MZ07),FAS(MZ22,MZ07) DIMENSION G1(2,MZ07),G2(2,MZ07),FLL(MZ07) COMMON /HOLD/DEL1,DEL2,DELNU,ELOW,EMAX,FINT,FRAT,HEAD,IDM,SZ,Z, C BLAM(MZ07,MZ07,5),EKLSQ(MZ01),ET(MZ07),ETERM(MZ02), C F1(2,MZ07,MZ07),F2(2,MZ07,MZ07),FI(MZ05,MZ05),FLAG(MZ09,3,4), C FX(2,MZ07,MZ07),P(MZ06,MZ03),Q(MZ06,MZ03),R(MZ06),S(MZ10,3), C SD1(MZ05,MZ05),SD2(MZ05,MZ05),SH(4,5),T(MZ10,3), C IFILE,IPRINT,ISTOT,ITC,ITCP1,ITMX,JFILE,KCASE,LAMMX,LCUT,LTOT, C MCH,MDG,MODE,N,N1,N2,NCHB,NCHCL,NCHF,NCHFMX,NCHOP,NLCTRN,NLOW, C NMX,NRAD,NTERM,NTHROW,NTOT,NTOTMX, C ISTERM(MZ02),ITERM(MZ07),LL(MZ07),LQN(MZ03),LTERM(MZ02), C MDEGEN(MZ01),MTCHCL(MZ01),NQN(MZ03) COMMON /WORK3/A(MZ07), H(MZ07),FNU(MZ07),PHI(MZ05) COMMON /ORDER/IROW(MZ08),ICOL(MZ12) COMMON /OUT/IPUNCH,ICONV COMMON /SCALE/CSCALE(MZ12),RSCALE(MZ08) COMMON /SCRATCH/ KFILE EQUIVALENCE (C(1,1),RFGST(1,1,1),FAS(1,1)),(FX(1,1,1),RFG(1,1)) X, (G1(1,1),F1(1,1,1)),(G2(1,1),F2(1,1,1)),(FX(1,1,1),FLL(1)) C C COMPUTE G=(GP IN UNPRIMED REPRESENTATION) C AND RETURN TO ORIGINAL NORMALIZATION IF(IPRINT.EQ.5) GO TO 5 IF(IPUNCH.EQ.0.AND.IPRINT.LE.3) GO TO 90 5 DO 10 I=1,NTOT K=ICOL(I) 10 G(K)=GP(I)*CSCALE(K) C C STORES MATRICES A AND H SUCH THAT C A(I)= LOWER CASE D FOR I=1,NCHF C H(I)= CAPITAL A FOR I=1,NCHF LI=0 DO 20 I=1,NCHF LI=LI+N1 A(I)=G(LI) LI=LI+1 H(I)=G(LI) 20 CONTINUE C C COMPUTE RADIAL FUNCTIONS AT POINTS (N+1) AND (N+2) J1=-1 J2=0 B2=FL/DEL1 B1=1.-B2 DO 30 I=1,NCHF J1=J1+N2 J2=J2+N2 A1=0. A2=0. DO 35 K=1,NCHF A3=A(K) A4=A3*B2 A3=A3*B1 A1=A1+F1(1,I,K)*A3+F1(2,I,K)*A4 35 A2=A2+F2(1,I,K)*A3+F2(2,I,K)*A4 G(J1)=A1 30 G(J2)=A2 C C PRINT THE RADIAL FUNCTIONS IF(IPRINT.LT.4) GO TO 90 M1=(NCHF-1)/10+1 WRITE(6,599) M2=1 DO 42 I1=1,M1 M3=M2+9 IF(M3.GT.NCHF) M3=NCHF WRITE(6,601) (L,L=M2,M3) WRITE(6,602) LI=(M3-1)*N2 M4=(M2-1)*N2 DO 40 IT=1,N2 LI=LI+1 M4=M4+1 WRITE(6,603) IT,(G(LIP),LIP=M4,LI,N2) IF (IT.EQ.ITC.OR.IT.EQ.N) WRITE(6,602) 40 CONTINUE 42 M2=M2+10 WRITE(6,624) (H(I),I=1,NCHF) LI=NCHF*N2 IF (NCHB.EQ.0) GO TO 45 WRITE(6,604) DO 50 J=1,NCHB I=LI+J 50 WRITE(6,612) J,G(I) 45 LI=LI+NCHB LIP=NTOT-LI IF (LIP.EQ.0) GO TO 60 WRITE(6,605) DO 70 J=1,LIP I=LI+J 70 WRITE(6,612) J,G(I) C C PRINT THE MATRIX SMALL D 60 WRITE(6,625) (A(I),I=1,NCHF) IF(ICONV.NE.0) GO TO 91 IF(NTHROW.EQ.0) WRITE(6,600) IF(NTHROW.EQ.1) WRITE(6,1600) RETURN C C OUTPUT ON JFILE 90 IF(ICONV.EQ.0) RETURN 91 IF(JFILE.EQ.0) RETURN C PUT CONVERGED ENERGY IN ET(I) DO 1000 I=1,NCHF 1000 ET(I)=ET(I)+FL C PUNCH KCASE,NCHOP AND CONVERGED ENERGY SQK=EKLSQ(KCASE)+FL WRITE(JFILE,730) KCASE,NCHOP,SQK IF(IPUNCH.NE.1) RETURN WRITE(6,640) C PUNCH CAP A WRITE(JFILE,700)(H(I),I=1,NCHF) C PUNCH SMALL G FOR IT=1,ITC AND F FOR IT GREATER THAN ITC L=1 DO 1030 I=1,NCHF M=L+N2-1 WRITE(JFILE,700)(G(K),K=L,M) 1030 L=L+N2 C PUNCH SMALL C IF(NCHB.EQ.0) GO TO 1035 L=M+1 M=L+NCHB -1 WRITE(JFILE,700)(G(K),K=L,M) C COMPUTE F AND STORE IN GP C COMPUTE FBAR AND STORE IN G 1035 L=1 DO 1050 I=1,NCHF A1=H(I) G(L)=A1 GP(L)=A1 DO 1040 IT=ITCP1,N2 K=L+IT-1 1040 GP(K)=G(K) 1050 L=L+N2 DO 1060 IT=2,ITC A1=R(IT) A2=A1*A1 A3=Z*A1 L=IT DO 1060 I=1,NCHF K=LL(I)+1 A4=K A5=H(I)*(1.-A3/A4)+A2*G(L) G(L)=A5 GP(L)=A5*(A1**K) 1060 L=L+N2 C COMPUTE Q AND STORE IN GPP REWIND 11 CNOAA READ(KFILE) C CALL INMAT(C,NTOTMX,NTOT,NTOT,KFILE) L=0 DO 1080 I=1,NCHF A1=ET(I) DO 1080 IT=1,N2 L=L+1 A4=-A1*GP(L) DO 1070 IP=1,NTOT 1070 A4=A4+C(L,IP)*G(IP) 1080 GPP(L)=-A4 C PUNCH Q L=1 DO 1090 I=1,NCHF M=L+N2-1 WRITE(JFILE,700)(GPP(K),K=L,M) 1090 L=L+N2 C C FUNCTIONS IN ASYMPTOTIC REGION COMPUTED AT R=R(N+1)+M*DR C WHERE M=1,(MFG+1) AND DR=.5*(R(N+2)-R(N+1)) C FUNCTIONS STORED IN FAS(M,I),I=1,NCHF C - DEFINE MFGP1=MFG+1 WRITE (JFILE,710) MFGP1 C - READ RFG FOR INPUT ENERGY AND STORE IN RFGST DO 1110 M=1,MFGP1 READ (11) RFG DO 1110 I=1,NCHF DO 1110 J=1,NCHF 1110 RFGST(M,I,J) = RFG(I,J) C TWO METHODS ARE USED C FOR EACH CHANNEL AND EACH VALUE OF R DEFINE C - WJ=2*SZ/R-LL*(LL+1)/R**2+ET(J) C METHOD = 1 IS USED IF WJ.LT.0 FOR ALL J C - INTERPOLATE BETWEEN RFG(INPUT ENERGY) AND C - RFG(INCREMENTED ENERGY). C - COMPUTE RADIAL FUNCTIONS FOR CONVERGED ENERGY C METHOD = 2 IS USED IF WJ.GT.0 FOR ANY J C - COMPUTE RADIAL FUNCTIONS FOR INPUT ENERGY C - AND FOR INCREMENTED ENERGY, C - AND THEN INTERPOLATE FOR FUNCTIONS AT CONVERGED ENERGY C METHOD FOR R(N1) HAS BEEN COMPUTED IN SUBROUTINE FG IF(METHOD.EQ.2) GO TO 1141 C COMPUTE FOR METHOD = 1 C INITIALIZE A, FUNCTIONS AT R(N+1) J=-1 DO 1115 I=1,NCHF J=J+N2 1115 A(I)=G(J) MM=1 C DO LOOP ON M FOR METHOD 1 1116 DO 1140 M=MM,MFGP1 C - READ RFG AT INCREMENTED ENERGY C - AND COMPUTE RFG AT CONVERGED ENERGY READ (11) RFG DO 1120 I=1,NCHF DO 1120 J=1,NCHF 1120 RFG(I,J) = RFGST(M,I,J)*B1 + RFG(I,J)*B2 DO 1130 I=1,NCHF A2 = 0. DO 1125 J=1,NCHF 1125 A2 = A2 + RFG(I,J)*A(J) 1130 FAS(M,I) = A2 DO 1135 I=1,NCHF 1135 A(I) = FAS(M,I) 1140 CONTINUE GO TO 1160 C COMPUTE FOR METHOD = 2 1141 RR=R(N1) DR=.5*(R(N2)-RR) C1=2.*SZ C INITIALIZE FUNCTIONS AT OLD R C G1(1,J) = FUNCTION AT INPUT ENERGY C G1(2,J) = FUNCTION AT INCREMENTED ENERGY C - INITIALIZE AT OLD R = R(N1) DO 1143 J=1,NCHF A1=0. A2=0. DO 1142 K=1,NCHF A3=A(K) A1=A1+F1(1,J,K)*A3 1142 A2=A2+F1(2,J,K)*A3 G1(1,J)=A1 1143 G1(2,J)=A2 C COMPUTE ARRAY FLL(J)=LL(J)*(LL(J)+1) DO 1144 J=1,NCHF A1=LL(J) 1144 FLL(J)=A1*(A1+1.) C DO LOOP ON M FOR METHOD 2 DO 1159 M=1,MFGP1 C COMPUTE NEW R RR=RR+DR C COMPUTE FUNCTIONS AT NEW R C NEW FUNCTIONS ARE C - G2(1,J) AT INPUT ENERGY C - G2(2,J) AT INCREMENTED ENERGY READ(11) RFG DO 1146 J=1,NCHF A1=0. A2=0. DO 1145 K=1,NCHF A1=A1+RFGST(M,J,K)*G1(1,K) 1145 A2=A2+RFG(J,K)*G1(2,K) G2(1,J)=A1 1146 G2(2,J)=A2 C INTERPOLATE FUNCTIONS DO 1147 J=1,NCHF 1147 FAS(M,J)=G2(1,J)*B1+G2(2,J)*B2 C CHECK METHOD IF(M.EQ.MFGP1) GO TO 1159 A1=1./RR A2=A1*A1 A1=C1*A1 METHOD=1 DO 1149 J=1,NCHF IF(KJ(J).EQ.1) GO TO 1149 WJ=A1-FLL(J)*A2+ET(J) IF(WJ.LT.0.) GO TO 1148 METHOD=2 GO TO 1149 1148 KJ(J)=1 1149 CONTINUE IF(METHOD.EQ.2) GO TO 1151 C CHANGE TO METHOD = 1 MM=M+1 DO 1150 J=1,NCHF 1150 A(J)=FAS(M,J) WRITE(6,630) M GO TO 1116 C COMPUTE FUNCTIONS AT NEW VALUE OF OLD R 1151 DO 1152 J=1,NCHF G1(1,J)=G2(1,J) 1152 G1(2,J)=G2(2,J) 1159 CONTINUE C C PUNCH FUNCTIONS IN ASYMPTOTIC REGION 1160 DO 1170 J=1,NCHF 1170 WRITE(JFILE,700)(FAS(M,J),M=1,MFGP1) C IF(NTHROW.EQ.0) WRITE(6,600) IF(NTHROW.EQ.1) WRITE(6,1600) RETURN C C PRINT FORMATS 599 FORMAT(1H0,31HEIGENSOLUTIONS OF THE EQUATIONS ) 600 FORMAT (1H1) 1600 FORMAT(/////) 601 FORMAT(//,20X,8HCHANNELS/(9X,9(I2,11X),I2)) 602 FORMAT (1H0) 603 FORMAT(1H ,I3,E12.4,9E13.4) 604 FORMAT(/14H C PARAMETERS/) 605 FORMAT(/16H LAM PARAMETERS/) 612 FORMAT(I3,E13.4) 624 FORMAT(/,12H CAPITAL A //(4X,E12.4,9E13.4)) 625 FORMAT(/,21H VALUES OF SMALL D ,//(4X,E12.4,9E13.4)) 630 FORMAT(//51H FOR CALCULATION OF FUNCTIONS IN ASYMPTOTIC REGION, C /32H CHANGE TO METHOD = 1 AFTER M = ,I4) 640 FORMAT(///51H CONVERGENCE CRITERION SATISFIED, FUNCTIONS WRITTEN C ,9H ON JFILE) C PUNCH FORMATS 700 FORMAT(5E14.7) 710 FORMAT(24I3) 730 FORMAT(2I3,E14.6,19H (=CONVERGED EKLSQ) ) END C*********************************************************************** SUBROUTINE QROP(NDRV,NBA,NAR,NAQ) C INCLUDE 'PARAM' C C ITERATION OPERATOR FOR NON-DEGENARATE ENERGIES C OPERATORS Q AND R ON A OR B FOR NAR OR NAQ 1 OR 2 RESP. C COMPUTES Q+R FOR NBA=1, Q-R FOR NBA=2 C DERIVATIVE NDRV-1 OF QR IN STATE I C RESULT OBTAINED FOR EACH MESH POINT C COMMON NCHD1,IJNZ(MZ25),KNZ(5),BLIN(MZ25),NZD,NZDM, C AB1(MZ07,13,10),AB2(MZ07,13,10),ALP(13,13),CC(MZ07,MZ07,5), C E2(MZ07),PP(12,10),T1(MZ07),V(MZ11,11,10),ZET(12,10),ZTIN(12,10) C ,GAM(2,MZ07,15) COMMON /CMN02/ NCNL,NOP,LAMAX COMMON /CMN03/ NDL(MZ07), NDU(MZ07) COMMON /CMN04/ HY,NHY,J COMMON /CMN14/ CH(10),CHN1(10),CHN2(10),CHN3(10),CHN4(10) SGN=2.0 IF(NBA.EQ.1) SGN=-2.0 NP=NDRV NDLJ=NDL(J) NDUJ=NDU(J) E2J=E2(J) DO 80 I=1,NCNL IF(I.LT.NDLJ) GO TO 25 IF(I.GT.NDUJ) GO TO 25 GO TO 80 25 II=I DE=E2(I)-E2J CALL DCHAIN(CHN1,II,NP,0,NAR,3,ZTIN,0,ZET,0) CALL DCHAIN(CHN2,II,NP,0,NAR,2,ZTIN,2,ZET,1) CALL DCHAIN(CHN3,II,NP,1,NAR,1,ZTIN,0,ZET,0) CALL DCHAIN(CHN4,II,NP,0,NAR,1,ZET,0,PP,1) CALL DCHAIN(CH,II,NP,0,NAQ,2,PP,0,ZET,1) IF(NBA.GT.1) GO TO 68 DO 70 IM=1,NHY 70 AB1(I,NP,IM)=-(CHN1(IM)+CHN2(IM)+CHN3(IM)+CHN4(IM)+SGN*CH(IM))/DE GO TO 80 68 DO 71 IM=1,NHY 71 AB2(I,NP,IM)=-(CHN1(IM)+CHN2(IM)+CHN3(IM)+CHN4(IM)+SGN*CH(IM))/DE 80 CONTINUE RETURN END C ********************************************************************* SUBROUTINE RLAM C INCLUDE 'PARAM' C C COMPUTES OPERATOR RL(IT) (SEE CPC APPENDIX B(III,D) C INTEGER HEAD COMMON /HOLD/DEL1,DEL2,DELNU,ELOW,EMAX,FINT,FRAT,HEAD,IDM,SZ,Z, C BLAM(MZ07,MZ07,5),EKLSQ(MZ01),ET(MZ07),ETERM(MZ02), C F1(2,MZ07,MZ07),F2(2,MZ07,MZ07),FI(MZ05,MZ05),FLAG(MZ09,3,4), C FX(2,MZ07,MZ07),P(MZ06,MZ03),Q(MZ06,MZ03),R(MZ06),S(MZ10,3), C SD1(MZ05,MZ05),SD2(MZ05,MZ05),SH(4,5),T(MZ10,3), C IFILE,IPRINT,ISTOT,ITC,ITCP1,ITMX,JFILE,KCASE,LAMMX,LCUT,LTOT, C MCH,MDG,MODE,N,N1,N2,NCHB,NCHCL,NCHF,NCHFMX,NCHOP,NLCTRN,NLOW, C NMX,NRAD,NTERM,NTHROW,NTOT,NTOTMX, C ISTERM(MZ02),ITERM(MZ07),LL(MZ07),LQN(MZ03),LTERM(MZ02), C MDEGEN(MZ01),MTCHCL(MZ01),NQN(MZ03) COMMON /WORK1/PP(MZ06),PP1(MZ06),PP2(MZ06),RHO(5) COMMON /OPS/Y(MZ06,MZ06),RL(MZ06) COMMON /OVSUB2/ IG,IGP,IGPP,IG1,L,IPQ,MOAD,LAM DIMENSION B(MZ06),AA(5),AB(4),AC(4) DIMENSION RHOL(5),RHOL1(5) C LAM1=-LAM-1 C L1=LQN(IG) L2=LQN(IGP) L3=LQN(IGPP) RITC=R(ITC) C C COMPUTE PP(IT),PP1(IT),PP2(IT) AS COEFFICIENTS FOR EXPANSION OF C P(IT,IG),P(IT,IGP),P(IT,IGPP) FOR IT=1,ITC DO 1 NN=1,ITC A1=0. A2=0. A3=0. DO 2 M=1,ITC A4=FI(NN,M) A1=A1+A4*P(M,IG) A2=A2+A4*P(M,IGP) 2 A3=A3+A4*P(M,IGPP) PP(NN)=A1 PP1(NN)=A2 1 PP2(NN)=A3 C C COMPUTE RL,B,CC, AND BB IN STARTING REGION, AS DEFINED C BY CPC EQN. (B29) M1=L+L1+L2+L3+1 RITCM1=RITC**M1 D1=L1+L3-LAM D2=L1+L3+LAM+1 DO 3 K=1,ITC A4=0. B2=1. DO 4 M=1,ITC B1=PP2(M) M0=M B2=B2*RITC E1=B2 DO 4 NN=1,ITC E1=E1*RITC M0=M0+1 A1=D1+M0 A2=D2+M0 A1=1./A1 A2=1./A2 B3=B1*PP(NN)*(A2-A1) E2=E1 DO 4 J=1,ITC M3=M1+M0+J E2=E2*RITC B4=B3*FI(J,K) A3=E2 DO 4 I=1,ITC A1=M3+I A1=B4/A1 A3=A3*RITC 4 A4=A4+A3*A1*PP1(I) 3 RL(K)=A4*RITCM1 C COMPUTE AND STORE B(K) , K=1,ITC M1=LAM+L+L2+1 RITCM1=RITC**M1 DO 6 K=1,ITC A1=0. A2=1. M0=M1 DO 7 J=1,ITC A2=A2*RITC B3=A2 M0=M0+1 M2=M0 B1=FI(J,K) DO 7 I=1,ITC B3=B3*RITC M2=M2+1 E2=M2 7 A1=A1+PP1(I)*B3*B1/E2 6 B(K)=A1*RITCM1 C COMPUTE AND STORE BB AND CC AT ITC BB=0. CC=0. A3=1. M1=0 DO 8 J=1,ITC B1=PP2(J) A3=A3*RITC M1=M1+1 B3=A3 M3=M1 DO 8 I=1,ITC M3=M3+1 B3=B3*RITC A1=D1+M3 A2=D2+M3 B4=B1*B3*PP(I) BB=BB+B4/A2 8 CC=CC+B4/A1 BB=BB*RITC**D2 CC=CC*RITC**D1 C C COMPUTE PP(IT),PP1(IT),PP(2(IT) FOR IT=(ITC-1),N AS RADIAL C FUNCTIONS P(IT,IG),P(IT,IGP),P(IT,IGPP) (NOT BAR FUNCTIONS) ITCM1=ITC-1 DO 9 I=ITCM1,ITC RIT=R(I) PP(I)=P(I,IG)*RIT**(L1+1) PP1(I)=P(I,IGP)*RIT**(L2+1) 9 PP2(I)=P(I,IGPP)*RIT**(L3+1) DO 10 I=ITCP1,N PP(I)=P(I,IG) PP1(I)=P(I,IGP) 10 PP2(I)=P(I,IGPP) C C COMPUTE RL,B AND CC AS DEFINED BY CPC EQN. (B30) DO 5 K=ITCP1,N B(K)=0. 5 RL(K)=0. C RIT=R(ITC) C1=RIT**LAM C2=1./(C1*RIT) C6=PP(ITC)*PP2(ITC) DO 12 IT=ITCP1,N IT3=IT-3 ITM1=IT-1 ITMC=IT-ITC JMAX=4 IF (IT.EQ.N) JMAX=3 RITM1=RIT RIT=R(IT) C3=C1 C4=C2 C5=C6 C1=RIT**LAM C2=1./(C1*RIT) C6=PP(IT)*PP2(IT) DR=0.25*(RIT-RITM1) RHO(1)=RITM1 DO 13 MU=2,5 13 RHO(MU)=RHO(MU-1)+DR DO 25 MU=1,5 B2=RHO(MU) B3=B2**LAM RHOL(MU)=B3 25 RHOL1(MU)=1./(B2*B3) DO 14 IS=1,4 B3=SH(IS,1)*C5 B4=SH(IS,5)*C6 B1=B3*C3+B4*C1 B2=B3*C4+B4*C2 DO 15 MU=2,4 MU1=MU-1 E1=RHOL(MU) E2=RHOL1(MU) D1=SH(IS,MU) DO 15 I=1,JMAX D2=D1*FLAG(ITMC,MU1,I)*PP(IT3+I) DO 15 J=1,JMAX D3=D2*FLAG(ITMC,MU1,J)*PP2(IT3+J) B1=B1+D3*E1 15 B2=B2+D3*E2 AB(IS)=BB+B1*DR 14 AC(IS)=CC+B2*DR AA(1)=C4*BB-C3*CC BB=AB(4) CC=AC(4) DO 16 IS=2,5 16 AA(IS)=RHOL1(IS)*AB(IS-1)-RHOL(IS)*AC(IS-1) A3=SH(4,1)*DR*PP1(ITM1) IF (IT.EQ.ITCP1) A3=A3*RITM1**(L+1) B(ITM1)=B(ITM1)+A3*C3 RL(ITM1)=RL(ITM1)+A3*AA(1) A3=SH(4,5)*DR*PP1(IT) B(IT)=B(IT)+A3*C1 RL(IT)=RL(IT)+A3*AA(5) DO 12 MU=2,4 MU1=MU-1 A1=RHOL(MU) A3=0. DO 17 J=1,JMAX 17 A3=A3+FLAG(ITMC,MU1,J)*PP1(IT3+J) A3=A3*SH(4,MU)*DR A1=A1*A3 A2=AA(MU)*A3 DO 12 J=1,JMAX IS=IT3+J A3=FLAG(ITMC,MU1,J) IF (IS.LE.ITC) A3=A3*R(IS)**(L+1) B(IS)=B(IS)+A1*A3 12 RL(IS)=RL(IS)+A2*A3 C C COMPLETE COMPUTATION OF RL , CPC EQN.7B31) DO 18 K=1,N 18 RL(K)=RL(K)+CC*B(K) RETURN END C*********************************************************************** SUBROUTINE SOLN C INCLUDE 'PARAM' C DESCRIBED IN CPC SECTION 8.2 C C CALLS ROUTINES FOR SOLVING THE EQUATIONS AND FOR CALCULATING AND C PRINTING FINAL RESULTS C INTEGER HEAD COMMON /HOLD/DEL1,DEL2,DELNU,ELOW,EMAX,FINT,FRAT,HEAD,IDM,SZ,Z, C BLAM(MZ07,MZ07,5),EKLSQ(MZ01),ET(MZ07),ETERM(MZ02), C F1(2,MZ07,MZ07),F2(2,MZ07,MZ07),FI(MZ05,MZ05),FLAG(MZ09,3,4), C FX(2,MZ07,MZ07),P(MZ06,MZ03),Q(MZ06,MZ03),R(MZ06),S(MZ10,3), C SD1(MZ05,MZ05),SD2(MZ05,MZ05),SH(4,5),T(MZ10,3), C IFILE,IPRINT,ISTOT,ITC,ITCP1,ITMX,JFILE,KCASE,LAMMX,LCUT,LTOT, C MCH,MDG,MODE,N,N1,N2,NCHB,NCHCL,NCHF,NCHFMX,NCHOP,NLCTRN,NLOW, C NMX,NRAD,NTERM,NTHROW,NTOT,NTOTMX, C ISTERM(MZ02),ITERM(MZ07),LL(MZ07),LQN(MZ03),LTERM(MZ02), C MDEGEN(MZ01),MTCHCL(MZ01),NQN(MZ03) COMMON CP(MZ08,MZ12) COMMON /OVSUB4/ NOROWS,NOCOLS COMMON /ORDER/ IROW(MZ08),ICOL(MZ12) COMMON /SCALE/CSCALE(MZ12),RSCALE(MZ08) C C ARRAY C CONTAINS MATRIX WITH NTOT ROWS AND NC COLUMNS WHERE NC=NTOT+NCHOP MCH=MTCHCL(KCASE) IF (MCH.LE.0) MCH=1 IF (MCH.GT.NCHF) MCH=1 C C CLOSE-UP CP AND G AND CHANGE SIGN OF G K1=NTOTMX K2=NTOT IF(NCHOP.EQ.0) GO TO 905 DO 91 J=1,NCHOP K1=K1+1 K2=K2+1 DO 91 I=1,NTOT 91 CP(I,K2)=-CP(I,K1) 905 CONTINUE C C FOR DEGENERATE COULOMB BOUND STATES, SET ALL CHANNEL MATRICES C EQUAL TO UNIT MATRICES EXCEPT FOR CHANNEL MDG IF (MODE.NE.2) GO TO 9 IF(NCHOP.NE.0) GO TO 9 MDG=MDEGEN(KCASE) IF(MDG.LE.0) GO TO 9 IF(MDG.GT.NCHF) GO TO 9 K2=0. DO 8 I=1,NCHF K1=K2+1 K2=K2+N2 IF(I.EQ.MDG) GO TO 8 DO 7 IS=K1,K2 DO 6 IT=K1,K2 6 CP(IS,IT)=0. 7 CP(IS,IS)=1. 8 CONTINUE 9 CONTINUE C NOROWS=NTOT NOCOLS=NC C C C NRB: CALL VSOLVE DRIVER FOR SOLUTION OF LINEAR ALGEBRAIC EQUATIONS C USING BLAS/LAPACK LIBRARIES. C CALL VSOLVE(CP,CSCALE,RSCALE,IROW,NTOTMX,NTOT,NC,IER) C C IF(IER.EQ.0) GO TO 930 WRITE(6,693) IER 693 FORMAT(//37H RETURN FROM MATRIX ROUTINE WITH IER=,I3) STOP 930 CONTINUE C C OPEN-UP C AND G IF(NCHOP.EQ.0) GO TO 60 K1=NCHOP+NTOTMX+1 K2=NCHOP+NTOT+1 DO 50 I=1,NCHOP K1=K1-1 K2=K2-1 DO 50 J=1,NTOT 50 CP(J,K1)=CP(J,K2) 60 CONTINUE C C CALLS OMEGA IF NCHOP.GT.0 - CALLS BOUND IF NCHOP.EQ.0 IF(NCHOP.EQ.0) GO TO 80 C -- CASE OF NCHOP.GT.0 CALL OMEGA RETURN 80 CONTINUE C -- CASE OF NCHOP.EQ.0 C C SUBROUTINE VSOLVB RETURNS IROW WHICH GIVES INFORMATION ON COLUMN SWAPPING C IN MATRIX COMPUTE ICOL WHICH GIVES INFORMATION FOR COLUMN INTERCHANGES C K1=NOROWS-1 DO 300 J=1,NOCOLS 300 ICOL(J)=J DO 310 I=1,K1 J=IROW(I) K=ICOL(J) ICOL(J)=ICOL(I) 310 ICOL(I)=K CALL BOUND RETURN END C*********************************************************************** SUBROUTINE SOLV(R,H,MD,WF,WF1) C INCLUDE 'PARAM' C CALLED BY SUBROUTINES ASYMPT AND ASYSM C C FOURTH ORDER RUNGE-KUTTA INTEGRATION OF COUPLED EQUATIONS C WF AND WF1 ARE STARTINE VALUES OF SOLUTION VECTOR AND DERIVATIVE C H = SIZE OF MESH POINTS, NEGATIVE FOR INWARD INTEGRATION C MD , NUMBER OF MESH POINTS TO MATCHING POINT RF=RS+MD*H C RETURNS WF AND WF1 VIA ARGUMENT LIST C FINAL VALUE OF R IS ALSO RETURNED C DIMENSION WF(MZ07),WF1(MZ07) C COMMON NCHD1,IJNZ(MZ25),KNZ(5),BLIN(MZ25),NZD,NZDM, C AB1(MZ07,13,10),AB2(MZ07,13,10),ALP(13,13),CC(MZ26), C E2(MZ07),PP(12,10),T1(MZ07),V(MZ11,11,10),ZET(12,10),ZTIN(12,10) C , W(MZ07,MZ07),C1(MZ07),C2(MZ07),D1(MZ07),D2(MZ07),D3(MZ07) COMMON /CMN02/ N,NOP,M COMMON /WTR/B1,B2 B1=0. B2=0. HH=.5*H C C SCALE E2,WF1 Z2=H*H IF(NZD.EQ.0) GO TO 105 A1=NZD*NZD Z2=A1*Z2 B2=2.*Z2 105 DO 10 I=1,N E2(I)=E2(I)*Z2 10 WF1(I)=WF1(I)*H C C SCALE CC AND PUT INTO BLIN FOR FAST ADDRESSING LK=0 K2=0 DO 15 K=1,M K1=K2+1 K2=KNZ(K) IF(K2.LT.K1) GO TO 15 DO 12 L=K1,K2 IJ=IJNZ(L) 12 BLIN(L)=CC(IJ+LK)*Z2 15 LK=LK+MZ24 C CALL WOUTER(R) C C START INWARD INTEGRATION DO 100 IJ=1,MD CRAY DO 1 I=1,N CRAY D1(I)=0. CRAY D2(I)=0. CRAY D3(I)=0. CRAY1 C2(I)=WF(I)+0.5*WF1(I) DO 3 I=1,N A1=0. DO 2 J=1,N 2 A1=A1+W(I,J)*WF(J) C2(I)=WF(I)+0.5*WF1(I) 3 D1(I)=A1 CRAY DO 3 J=1,N CRAY3 D1(J)=D1(J)+W(J,I)*WF(I) R=R+HH CALL WOUTER(R) CRAY DO 2 I=1,N CRAY2 C1(I)=C2(I)-0.25*D1(I) DO 5 I=1,N A2=D1(I) A1=0. DO 4 J=1,N 4 A1=A1+W(I,J)*C2(J) D2(I)=A1 C1(I)=C2(I)-0.25*A2 5 D1(I)=A1+A2 CRAY DO 5 J=1,N CRAY5 D2(J)=D2(J)+W(J,I)*C2(I) CRAY DO 4 I=1,N CRAY D1(I)=D1(I)+D2(I) CRAY4 C2(I)=WF(I)+WF1(I)-0.5*D2(I) DO 7 I=1,N A1=0. A2=D2(I) DO 6 J=1,N 6 A1=A1+W(I,J)*C1(J) D3(I)=A1 C2(I)=WF(I)+WF1(I)-0.5*A2 7 D2(I)=A2+2.*A1 CRAY DO 7 J=1,N CRAY7 D3(J)=D3(J)+W(J,I)*C1(I) CRAY DO 6 I=1,N CRAY6 D2(I)=D2(I)+2.*D3(I) R=R+HH CALL WOUTER(R) A3=1./6. DO 9 I=1,N A2=WF1(I) A1=0. DO 8 J=1,N 8 A1=A1+W(I,J)*C2(J) WF(I)=WF(I)+A2-(D1(I)+D3(I))*A3 9 WF1(I)=A2-(D1(I)+D2(I)+A1)*A3 CRAY WF(I)=WF(I)+WF1(I)-(D1(I)+D3(I))*A3 CRAY DO 9 J=1,N CRAY9 D2(J)=D2(J)+W(J,I)*C2(I) CRAY DO 8 I=1,N CRAY8 WF1(I)=WF1(I)-A3*(D1(I)+D2(I)) 100 CONTINUE C END INWARD INTEGRATION C C UNSCALE E2,WF1 A1=1./H Z2=1./Z2 DO 110 I=1,N WF1(I)=WF1(I)*A1 110 E2(I)=E2(I)*Z2 RETURN END C*********************************************************************** SUBROUTINE SUB1 INCLUDE 'PARAM' CALL ATOM CALL MATS RETURN END C*********************************************************************** SUBROUTINE SUB2 INCLUDE 'PARAM' CALL EQNS RETURN END C*********************************************************************** SUBROUTINE SUB3 INCLUDE 'PARAM' COMMON /OV/SQK,IERROR,NT,NC CALL ASYMPT(SQK,IERROR) RETURN END C*********************************************************************** SUBROUTINE SUB4 INCLUDE 'PARAM' CALL ENERGY CALL SOLN RETURN END C*********************************************************************** SUBROUTINE TIMES(IPOS) C INCLUDE 'PARAM' C CNRB SAVE TIM0 C C WS AND PC - TIME SINCE LAST CALL REAL*4 TARRY(2) C DUM=DTIME(TARRY) FT=TARRY(1) FT=FT/60. IF(IPOS.EQ.0)THEN TIM0=FT RETURN ELSE TIM0=TIM0+FT ENDIF C WS AND PC C C CRAY - AGGREGATE TIME C IF(IPOS.EQ.0)THEN C CALL SECOND(TIM0) C TIM0=TIM0/60. C RETURN C ENDIF C CALL SECOND(TIM1) C TIM1=TIM1/60. C FT=TIM1-TIM0 C TIM0=TIM1 C CRAY C IF(IPOS.EQ.5)THEN WRITE(6,2) TIM0 RETURN ENDIF WRITE(6,1) IPOS,FT 1 FORMAT(//90X,'TIME FOR SUB',I1,' =',F10.3,' MIN'//) 2 FORMAT(//90X,'TOTAL RUN TIME =',F10.3,' MIN'//) RETURN END C*********************************************************************** SUBROUTINE VSOLVB(A,SC,SR,IP,IA,M,N,IER) C INCLUDE 'PARAM' C CNRB C SUBROUTINE TO LU DECOMPOSE MATRIX LINEAR SYSTEM C*G = 0 C C PARAMETERS...A AN ARRAY DIMENSIONED IA BY AT LEAST N. C COLUMNS 1 THROUGH M CONTAIN C, C C SC ARRAY OF LENGTH N TO HOLD THE C COLUMN SCALING FACTORS C C SR ARRAY OF LENGTH M TO HOLD THE C ROW SCALING FACTORS C C IP ARRAY OF LENGTH M TO HOLD C COLUMN PIVOT INFORMATION C C IA ROW DIMENSION OF ARRAY A IN C CALLING ROUTINE C IA MUST BE .GE. M C C M NUMBER OF ROWS AND COLUMNS IN C, C AND NUMBER OF ROWS IN G C C N M+NRHS, WHERE NRHS IS THE NUMBER C OF RIGHT HAND SIDES (NUMBER OF C COLUMNS IN G) = 0 HERE C C IER ERROR INDICATION PARAMETER C C ON ENTRY, PARAMETERS A,IA,M,N SHOULD BE SET C BY CALLING ROUTINE. C ON EXIT, PARAMETERS A,SC,SR,IP,IER WILL BE SET C (OR RESET) BY THIS ROUTINE. C C THIS ROUTINE PERFORMS AN L-U DECOMPOSITION OF THE C FIRST M COLUMNS OF A COLUMNWISE PERMUTATION OF THE C MATRIX A, USING THE CROUT ALGORITHM. C C SEE VSOLF FOR CASE OF N > M. C C VSOLVE CALLS THE BLAS ROUTINES C ISAMAX, SSCAL, SSWAP, AND SGEMV - NRB 05OCT99 C C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) C DIMENSION A(IA,N),SC(N),SR(M),IP(M) DATA EPSP/1.0E-10/ C IER = 0 C C CHECK CONSISTENCY OF INPUT PARAMETERS IF(N .LT. M) IER = -1 IF(IA .LT. M) IER = -2 IF(M .LT. (N-M)) IER = -3 IF(N-M .GT. 0)IER=-4 ! CASE NCHOP>0 IF(IER .NE. 0) RETURN C ALF=-ONE BET=ONE C C SCALE COLUMNS OF A DO J = 1,N CS MC = ISAMAX(M,A(1,J),1) MC = IDAMAX(M,A(1,J),1) CMAX = A(MC,J) IF(CMAX .EQ. TZERO) CMAX = ONE SC(J) = ONE/CMAX CS CALL SSCAL(M,ONE/CMAX,A(1,J),1) CALL DSCAL(M,ONE/CMAX,A(1,J),1) ENDDO C C SCALE ROWS OF A DO I = 1,M CS MR = ISAMAX(N,A(I,1),IA) MR = IDAMAX(N,A(I,1),IA) RMAX = A(I,MR) IF(RMAX .EQ. TZERO) RMAX = ONE SR(I) = ONE/RMAX CS CALL SSCAL(N,ONE/RMAX,A(I,1),IA) CALL DSCAL(N,ONE/RMAX,A(I,1),IA) ENDDO C C FOR EACH COLUMN, PERFORM THE REDUCTION APPROPRIATE TO LU C FACTORIZATION DO J = 2,M JM1 = J - 1 CS MR = ISAMAX(N+1-JM1,A(JM1,JM1),IA) + J - 2 MR = IDAMAX(N+1-JM1,A(JM1,JM1),IA) + J - 2 IP(JM1) = MR IF(MR .NE. JM1) THEN C C SWAP COMPLETE COLUMNS CS CALL SSWAP(M,A(1,JM1),1,A(1,MR),1) CALL DSWAP(M,A(1,JM1),1,A(1,MR),1) ENDIF C C OBTAIN COLUMN J-1 OF MATRIX L IF(JM1 .GT. 1) CS X CALL SGEMV('N',M-JM1,J-2,ALF,A(J,1),IA,A(1,JM1),1,BET,A(J,JM1),1) X CALL DGEMV('N',M-JM1,J-2,ALF,A(J,1),IA,A(1,JM1),1,BET,A(J,JM1),1) C C DIVIDE PRECEEDING COLUMN OF L BY A(JM1,JM1) IF(ABS(A(JM1,JM1)) .LE. EPSP) THEN IER = JM1 RETURN ENDIF CS CALL SSCAL(M-JM1,ONE/A(JM1,JM1),A(J,JM1),1) CALL DSCAL(M-JM1,ONE/A(JM1,JM1),A(J,JM1),1) C C OBTAIN ROW J OF MATRIX U CS CALL SGEMV('T',JM1,N-JM1,ALF,A(1,J),IA,A(J,1),IA,BET,A(J,J),IA) CALL DGEMV('T',JM1,N-JM1,ALF,A(1,J),IA,A(J,1),IA,BET,A(J,J),IA) ENDDO C IF(ABS(A(M,M)) .LT. EPSP) IER = M C RETURN END C*********************************************************************** SUBROUTINE VSOLVE(A,SC,SR,IP,IA,M,N,IER) C INCLUDE 'PARAM' C CNRB C DRIVER ROUTINE TO PICK APPROPRIATE VSOLV: C VSOLVF N-M>0 C VSOLVB N-M=0 C DIMENSION A(IA,N),SC(N),SR(M),IP(M) C IF(N.GT.M)CALL VSOLVF(A,SC,SR,IP,IA,M,N,IER) IF(N.EQ.M)CALL VSOLVB(A,SC,SR,IP,IA,M,N,IER) IF(N.LT.M)IER=-1 C RETURN END C*********************************************************************** SUBROUTINE VSOLVF(A,SC,SR,IP,IA,M,N,IER) C INCLUDE 'PARAM' C CNRB C SUBROUTINE FOR THE SOLUTION OF THE C MATRIX LINEAR SYSTEM C*G = X*B (B=1 DEFAULT) C (USES LAPACK/BLAS LIBRARIES) C C PARAMETERS...A AN ARRAY DIMENSIONED IA BY AT LEAST N. C COLUMNS 1 THROUGH M CONTAIN C, C COLUMNS M+1 THROUGH N CONTAIN -X C C SC ARRAY OF LENGTH N TO HOLD THE C COLUMN SCALING FACTORS C C SR ARRAY OF LENGTH M TO HOLD THE C ROW SCALING FACTORS C C IP ARRAY OF LENGTH M TO HOLD C ROW PIVOT INFORMATION C C IA ROW DIMENSION OF ARRAY A IN C CALLING ROUTINE C IA MUST BE .GE. M C C M NUMBER OF ROWS AND COLUMNS IN C, C AND NUMBER OF ROWS IN G AND X C C N M+NRHS, WHERE NRHS IS THE NUMBER C OF RIGHT HAND SIDES (NUMBER OF C COLUMNS IN G AND X) C C IER ERROR INDICATION PARAMETER C C ON ENTRY, PARAMETERS A,IA,M,N SHOULD BE SET C BY CALLING ROUTINE. C ON EXIT, PARAMETERS A,SC,SR,IP,IER WILL BE SET C (OR RESET) BY THIS ROUTINE. C C THIS ROUTINE PERFORMS AN L-U DECOMPOSITION OF THE C FIRST M COLUMNS OF MATRIX A USING PARTIAL PIVOTING. C THE FREE VARIABLES ARE CHOSEN TO BE THE C UNIT MATRIX, AND THE N-M>0 LINEARLY INDEPENDENT C SOLUTIONS OF THE SYSTEM (THE COLUMNS OF G) ARE C DETERMINED AND STORED IN COLUMNS M+1 THROUGH N C IN THE MATRIX A. C THE SQUARE MATRIX B IS RETURNED IN COLUMNS 1 C THROUGH NRHS IN ARRAY A. C C SEE VSOLB FOR CASE OF N = M. C C VSOLVF CALLS THE LAPACK ROUTINES SGETRF, SGETRS C AND THE BLAS ROUTINES ISAMAX, SSCAL - NRB 05OCT99 C C PARAMETER (TZERO=0.0) PARAMETER (ONE=1.0) C DIMENSION A(IA,N),SC(N),SR(M),IP(M) C IER = 0 C C CHECK CONSISTENCY OF INPUT PARAMETERS IF(N .LT. M) IER = -1 IF(IA .LT. M) IER = -2 IF(M .LT. (N-M)) IER = -3 IF(N-M .EQ. 0)IER=-4 ! CASE NCHOP=0 IF(IER .NE. 0) RETURN C C INITIAL SET-UP MPLUS1 = M + 1 NRHS = N - M C C SCALE MATRIX X ALF=-ONE BET=ONE TAB=ALF*BET DO J=MPLUS1,N CS CALL SSCAL(M,TAB,A(1,J),1) CALL DSCAL(M,TAB,A(1,J),1) ENDDO C C SCALE COLUMNS OF A DO J = 1,N CS MC = ISAMAX(M,A(1,J),1) MC = IDAMAX(M,A(1,J),1) CMAX = A(MC,J) IF(CMAX .EQ. TZERO) CMAX = ONE SC(J) = ONE/CMAX CS CALL SSCAL(M,ONE/CMAX,A(1,J),1) CALL DSCAL(M,ONE/CMAX,A(1,J),1) ENDDO C C SCALE ROWS OF A DO I = 1,M CS MR = ISAMAX(N,A(I,1),IA) MR = IDAMAX(N,A(I,1),IA) RMAX = A(I,MR) IF(RMAX .EQ. TZERO) RMAX = ONE SR(I) = ONE/RMAX CS CALL SSCAL(N,ONE/RMAX,A(I,1),IA) CALL DSCAL(N,ONE/RMAX,A(I,1),IA) ENDDO C C PERFORM THE REDUCTION APPROPRIATE TO LU FACTORIZATION CS CALL SGETRF(M,M,A,IA,IP,IER) CALL DGETRF(M,M,A,IA,IP,IER) IF(IER .NE. 0) RETURN C C SOLVE THE LU-DECOMPOSED SYSTEM C CHOOSE THE FREE VARIABLES WITH THE UNIT MATRIX CS CALL SGETRS('N',M,NRHS,A,IA,IP,A(1,MPLUS1),IA,IER) CALL DGETRS('N',M,NRHS,A,IA,IP,A(1,MPLUS1),IA,IER) IF(IER .NE. 0) RETURN C C ENTER AN NRHS BY NRHS UNIT MATRIX INTO A DO J = 1,NRHS DO I = 1,NRHS A(I,J) = TZERO ENDDO A(J,J) = ONE ENDDO C C RESCALE SOLUTION VECTORS BY COLUMN SCALING FACTORS DO JJ = 1,N J = JJ ICOL = MPLUS1 IF(JJ .GT. M) THEN ICOL = 1 J = J - M ENDIF CS CALL SSCAL(NRHS,SC(JJ),A(J,ICOL),IA) CALL DSCAL(NRHS,SC(JJ),A(J,ICOL),IA) ENDDO C RETURN END C*********************************************************************** SUBROUTINE WOUTER(R) C INCLUDE 'PARAM' C C ENERGY MATRIX IN THE ASYMPTOTIC REGION C COMMON NCHD1,IJNZ(MZ25),KNZ(5),BLIN(MZ25),NZD,NZDM, C AB1(MZ07,13,10),AB2(MZ07,13,10),ALP(13,13),CC(MZ07,MZ07,5), C E2(MZ07),PP(12,10),T1(MZ07),V(MZ11,11,10),ZET(12,10),ZTIN(12,10) C , W(MZ24),C1(MZ07),C2(MZ07),D1(MZ07),D2(MZ07),D3(MZ07) COMMON /CMN02/ N,NOP,M COMMON /WTR/B1,B2 X=1./R IF(NZD.EQ.0) GO TO 1 X=X/FLOAT(NZD) B1=B2*X C C INITIALISE W ON AND BELOW DIAGONAL 1 L=1 IF(N.EQ.1) GO TO 20 DO 15 J=2,N W(L)=E2(J-1)+B1 DO 10 I=J,N L=L+1 10 W(L)=0. 15 L=L+NCHD1+J 20 W(L)=E2(N)+B1 C C CONTINUE CALCULATION OF W ON AND BELOW DIAGONAL A2=X K2=0 DO 40 K=1,M A2=A2*X K1=K2+1 K2=KNZ(K) IF(K2.LT.K1) GO TO 40 DO 30 L=K1,K2 IJ=IJNZ(L) 30 W(IJ)=W(IJ)+A2*BLIN(L) 40 CONTINUE C C W ABOVE DIAGONAL IF(N.EQ.1) RETURN L=1 DO 50 I=2,N IJ=L DO 45 J=I,N L=L+1 IJ=IJ+MZ07 45 W(IJ)=W(L) 50 L=L+NCHD1+I C RETURN END C*********************************************************************** SUBROUTINE YLAM(IG,L,LAM) C INCLUDE 'PARAM' C C COMPUTES OPERATOR Y(IT,IS) USING METHOD DESCRIBED IN CPC, C APPENDIX B(III,B) C INTEGER HEAD COMMON /HOLD/DEL1,DEL2,DELNU,ELOW,EMAX,FINT,FRAT,HEAD,IDM,SZ,Z, C BLAM(MZ07,MZ07,5),EKLSQ(MZ01),ET(MZ07),ETERM(MZ02), C F1(2,MZ07,MZ07),F2(2,MZ07,MZ07),FI(MZ05,MZ05),FLAG(MZ09,3,4), C FX(2,MZ07,MZ07),P(MZ06,MZ03),Q(MZ06,MZ03),R(MZ06),S(MZ10,3), C SD1(MZ05,MZ05),SD2(MZ05,MZ05),SH(4,5),T(MZ10,3), C IFILE,IPRINT,ISTOT,ITC,ITCP1,ITMX,JFILE,KCASE,LAMMX,LCUT,LTOT, C MCH,MDG,MODE,N,N1,N2,NCHB,NCHCL,NCHF,NCHFMX,NCHOP,NLCTRN,NLOW, C NMX,NRAD,NTERM,NTHROW,NTOT,NTOTMX, C ISTERM(MZ02),ITERM(MZ07),LL(MZ07),LQN(MZ03),LTERM(MZ02), C MDEGEN(MZ01),MTCHCL(MZ01),NQN(MZ03) C COMMON /WORK1/B(MZ06),C(MZ06),PP(MZ06),RHO(5) COMMON /OPS/Y(MZ06,MZ06),YPOT(MZ06) C L1=LQN(IG) M=L1+L+2 LAM1=LAM+M-1 LAM2=-LAM+M-2 DO 50 IT=2,ITC DO 50 IS=ITCP1,N 50 Y(IT,IS)=0. C C COMPUTE PP AS POWER SERIES COEFFICIENTS FOR P(IT,IG) DO 1 I=1,ITC A1=0. DO 2 J=1,ITC 2 A1=A1+FI(I,J)*P(J,IG) 1 PP(I)=A1 C C INITIALISE ( B(IS), C(IS), Y(IT,IS),IT=2,ITC ) ,IS=1,ITC DO 3 IT=2,ITC RIT=R(IT) DO 4 IS=1,ITC B(IS)=0. 4 C(IS)=0. RI=1. DO 5 I=1,ITC RI=RI*RIT RJ=RI*PP(I) D1=LAM1+I D2=LAM2+I DO 5 J=1,ITC RJ=RJ*RIT D1=D1+1. D2=D2+1. A1=RJ/D1 A2=RJ/D2 DO 5 IS=1,ITC A=FI(J,IS) B(IS)=B(IS)+A*A1 5 C(IS)=C(IS)+A*A2 RIT=RIT**(M-2) DO 3 IS=1,ITC 3 Y(IT,IS)=RIT*(B(IS)-C(IS)) RIT=R(ITC) A1=RIT**LAM1 A2=RIT**LAM2 DO 6 IS=1,ITC B(IS)=B(IS)*A1 6 C(IS)=C(IS)*A2 DO 7 IS=ITCP1,N B(IS)=0. 7 C(IS)=0. C C COMPUTE PP AS RADIAL FUNCTION P(R,IG) (NOT BAR FUNCTION) DO 8 I=2,ITC 8 PP(I)=P(I,IG)*R(I)**(L1+1) DO 9 I=ITCP1,N 9 PP(I)=P(I,IG) C C COMPUTE AND STORE (Y(IT,IS),IT=2,N),IS=1,N LAM1=-LAM-1 A=RIT**(L+1) C1=RIT**LAM C2=A/(C1*RIT) C1=A*C1 DO 10 IT=ITCP1,N IT1=IT-1 IT3=IT-3 JMAX=4 IF (IT.EQ.N) JMAX=3 ITMC=IT-ITC RITM1=RIT D1=C1 D2=C2 RIT=R(IT) C1=RIT**LAM C2=1./(C1*RIT) DR=0.25*(RIT-RITM1) C COMPUTE RHO RHO(1)=RITM1 DO 11 MU=2,5 11 RHO(MU)=RHO(MU-1)+DR C A=SH(4,1)*DR*PP(IT1) B(IT1)=B(IT1)+A*D1 C(IT1)=C(IT1)+A*D2 C A=SH(4,5)*DR*PP(IT) B(IT)=B(IT)+A*C1 C(IT)=C(IT)+A*C2 C DO 13 MU=2,4 MU1=MU-1 RHOMU=RHO(MU) A=0. DO 12 J=1,JMAX 12 A=A+FLAG(ITMC,MU1,J)*PP(IT3+J) A=A*SH(4,MU)*DR A1=RHOMU**LAM A2=A/(A1*RHOMU) A1=A*A1 DO 13 J=1,JMAX IS=IT3+J A=FLAG(ITMC,MU1,J) IF (IS.LE.ITC) A=A*R(IS)**(L+1) B(IS)=B(IS)+A1*A 13 C(IS)=C(IS)+A2*A DO 10 IS=1,N 10 Y(IT,IS)=B(IS)*C2-C(IS)*C1 C C COMPLETE CALCULATION OF Y(IT,IS) FOR IT=2,N 31 DO 32 IT=2,N A1=R(IT)**LAM DO 32 IS=1,N 32 Y(IT,IS)=Y(IT,IS)+C(IS)*A1 C C COMPUTE Y(1,IS) DO 33 IS=1,N 33 Y(1,IS)=C(IS) RETURN END C*********************************************************************** SUBROUTINE ZETA(N) C INCLUDE 'PARAM' C C WBK AMPLITUDE ZET AND ITS INVERSE FOR ALL MESH POINTS FOR ITERA C COMPUTES AND STORES ZET(J) AND 1.0/ZET(J) AND DERIVS. UP TO N-1 C COMMON NCHD1,IJNZ(MZ25),KNZ(5),BLIN(MZ25),NZD,NZDM, C AB1(MZ07,13,10),AB2(MZ07,13,10),ALP(13,13),CC(MZ07,MZ07,5), C E2(MZ07),PP(12,10),T1(MZ07),V(MZ11,11,10),ZET(12,10),ZTIN(12,10) C , GAM(2,MZ07,15) DIMENSION BLP(12,11) COMMON /CMN04/ HY,NHY,J SGN=1. E2J=E2(J) AE2J=ABS(E2J) IF(E2J.LT.0.) SGN=-1. S=SGN+SGN C Y1=0. DO 40 IPT=1,NHY Y1=Y1+HY Y2=SQRT(AE2J+S*Y1*Y1) ZET(1,IPT)=Y2 40 ZTIN(1,IPT)=1./Y2 C IF(N.EQ.1) RETURN C S=SGN DO 100 NZ=1,2 S=-S ALP(1,1)=-S IF(N.EQ.2) GO TO 85 ALP(1,2)=S+S ST=SGN*(1.-NZ-NZ) ALP(2,2)=-(ST+4.*SGN)*S IF(N.EQ.3) GO TO 85 DO 170 KK=4,N K=KK-1 KM1=K-1 ALP(1,K)=-(1+KM1)*ALP(1,KM1) ALP(K,K)=(ST+SGN*(K+K))*ALP(KM1,KM1) DO 170 KP=2,KM1 170 ALP(KP,K)=-(KP+KM1)*ALP(KP,KM1)+(ST+SGN*(KP+KP))*ALP(KP-1,KM1) C 85 IF(NZ.EQ.2) GO TO 105 N1=N-1 DO 90 KK=1,N1 DO 90 KP=1,KK 90 BLP(KP,KK)=ALP(KP,KK) 100 CONTINUE C 105 Y2=0. DO 120 IPT=1,NHY Z1=ZET(1,IPT) Y2=Y2+HY Y1=Y2*Y2 YX=1. TEMP=Y1/(Z1*Z1) YMIN=1.0D-76/TEMP K=0 DO 120 KK=2,N K=K+1 YX=Y1*YX YZ=YX Y3=0. Y4=0. DO 88 KP=1,K IF(YZ.LT.YMIN) GO TO 89 YZ=YZ*TEMP Y3=Y3+BLP(KP,K)*YZ 88 Y4=Y4+ALP(KP,K)*YZ 89 ZTIN(KK,IPT)=Y3/Z1 120 ZET(KK,IPT)=Y4*Z1 C RETURN END