C C Program for coupled-channels equations for heavy-ion C fusion reactions C C The couplings are always treated to all orders. C C last modification: April 29, 1999 C IMPLICIT REAL*8 (A-H,O-Z) C C Parameters: C NLEVELMAX - maximum number of intrinsic levels to be included C PARAMETER (NLEVELMAX=25) C C Principal global variables C C P - penetrability C SIGMA - fusion cross section, unit mb C SPIN - mean angular momentum C C AP,ZP,RP - atomic #, proton # and radius of the projectile C AT,ZT,RT - those of the target C RMASS - reduced mass, unit MeV/c**2 C ANMASS - nucleon mass, unit MeV/c**2 C E - bombarding energy in the center of mass frame, unit MeV C V0,R0,A0 - depth, range, and dissuseness parameters of uncoupled C nuclear potential, which is assumed to be a Woods-Saxon form C IVIBROT - option for intrinsic degree of freedom C = -1; no excitation (inert)/ =0 ; vibrational coupling/ C =1; rotational coupling C BETA - defeormation parameter C OMEGA - excitation energy of the oscillator C LAMBDA - multipolarity C NPHONON - the number of phonon to be included C BETA2 - static quadrupole deformation parameter C BETA4 - static hexadecapole deformation parameter C E2 - excitation energy of 2+ state in a rotational band C NROT - the number of levels in the rotational band to be included C (up to I^pi=2*NROT+ states are included) C L - angular momentum of the relative motion C CPOT - coupling matrix C DIMENSION CPOT0(NLEVELMAX,NLEVELMAX) COMMON/HION/AP,ZP,RP,AT,ZT,RT COMMON/CONST/HBAR,PI COMMON/INTR/IVIBROTP,IVIBROTT,NPROJ,NTARG COMMON/OSCP/BETAP,BETAPN,OMEGAP,LAMBDAP,NPHONONP COMMON/ROTP/BETA2P,BETA4P,E2P,NROTP COMMON/OSCT/BETAT,BETATN,OMEGAT,LAMBDAT,NPHONONT COMMON/OSCT2/BETAT2,BETAT2N,OMEGAT2,LAMBDAT2,NPHONONT2 COMMON/ROTT/BETA2T,BETA4T,E2T,NROTT COMMON/DIMEN/NLEVEL COMMON/POT/V0,R0,A0 COMMON/DYN/E,RMASS COMMON/ANGULAR/L COMMON/SHAPE/RB,VB,CURV COMMON/GRID/RMIN,RMAX,DR COMMON/COUP/CPOT(NLEVELMAX,NLEVELMAX,0:3000) COMMON/COUPH/CPOTH(NLEVELMAX,NLEVELMAX) COMMON/TRANS/FTR,QTRANS,NTRANS COMMON/MUT/IMUTUAL(0:NLEVELMAX,0:NLEVELMAX) COMMON/EIGEN/EV(NLEVELMAX),EVEC(NLEVELMAX,NLEVELMAX) C C Output files: C OUTPUT: detailed information of the calculation C cross.dat : fusion cross sections C spin.dat : averaged angular momenta C partial.dat: partial cross sections (being calculated only when a single C value of the energy is entered in the input file) C OPEN(7,FILE='cross.dat',STATUS='UNKNOWN') OPEN(8,FILE='spin.dat',STATUS='UNKNOWN') OPEN(9,FILE='OUTPUT',STATUS='UNKNOWN') OPEN(10,FILE='ccfull.inp',STATUS='UNKNOWN') C C Define two constants used in various places C HBAR=197.329D0 PI=3.141592653D0 C C Input phase C ANMASS=938.D0 AP=16. ZP=8. AT=148. ZT=62. READ(10,*)AP,ZP,AT,ZT C r_coup, unit fm RP=1.2D0*AP**(1.D0/3.D0) RT=1.06D0*AT**(1.D0/3.D0) READ(10,*)RP,IVIBROTP,RT,IVIBROTT RP=RP*AP**(1.D0/3.D0) RT=RT*AT**(1.D0/3.D0) RMASS=AP*AT/(AP+AT)*ANMASS C----- c IVIBROTT=0 C = -1 for inert, =0 for vibration, =1 for rotation IF(IVIBROTT.EQ.-1) THEN NTARG=0 READ(10,*)OMEGA,BETA,ALAMBDA,NPHONON ELSEIF(IVIBROTT.EQ.0) THEN C (input for target phonon) OMEGAT=0.55D0 BETAT=0.182D0 LAMBDAT=2 NPHONONT=3 READ(10,*)OMEGAT,BETAT,LAMBDAT,NPHONONT NTARG=NPHONONT ELSE C (input for target rotation) E2T=0.D0 BETA2T=0.D0 BETA4T=0.D0 NROTT=0 LAMBDAT=2 READ(10,*)E2T,BETA2T,BETA4T,NROTT NTARG=NROTT ENDIF C----- C (input for target phonon (the second mode of excitation)) NPHONONT2=3 C =0; no second mode in the target OMEGAT2=1.16D0 BETAT2=0.236D0 LAMBDAT2=3 READ(10,*)OMEGAT2,BETAT2,LAMBDAT2,NPHONONT2 IF(IVIBROTT.EQ.-1) NPHONONT2=0 C----- C IVIBROTP=-1 C = -1 for inert, =0 for vibration, =1 for rotation IF(IVIBROTP.EQ.-1) THEN NPROJ=0 READ(10,*)OMEGAP,BETAP,ALAMBDAP,NPHONONP ELSE IF(IVIBROTP.EQ.0) THEN C (input for projectile phonon) OMEGAP=3.74D0 BETAP=0.339D0 LAMBDAP=3 NPHONONP=1 READ(10,*)OMEGAP,BETAP,LAMBDAP,NPHONONP NPROJ=NPHONONP ELSE C (input for projectile rotation) E2P=0.D0 BETA2P=0.D0 BETA4P=0.D0 NROTP=0 READ(10,*)E2P,BETA2P,BETA4P,NROTP NPROJ=NROTP ENDIF NTRANS=0 C NTRANS= 0 ; no transfer channel/ =1 ; with transfer channel C FTRANS(R)=FTR*dVN(R)/dR C FTR=7.D0 QTRANS=5.21D0 READ(10,*)NTRANS,QTRANS,FTR C--------------- parameters for the nuclear potential R1=1.233D0*AP**(1.D0/3.D0)-0.98D0*AP**(-1.D0/3.D0) R2=1.233D0*AT**(1.D0/3.D0)-0.98D0*AT**(-1.D0/3.D0) R0=R1+R2+0.29D0 R12=R1*R2/(R1+R2) A0=0.63D0 GAMMA=0.95D0*(1.D0-1.8D0*(AP-2.D0*ZP)/AP*(AT-2.D0*ZT)/AT) V0=16.D0*PI*GAMMA*R12*A0 v0=31.67*r1*r2/(r1+r2) V0=1551D0 R0=0.95D0*(AP**(1.D0/3.D0)+AT**(1.D0/3.D0)) A0=1.05D0 READ(10,*)V0,R0,A0 R0=R0*(AP**(1.D0/3.D0)+AT**(1.D0/3.D0)) C the minimum energy, the maximum energy, and the energy step EMIN=53.D0 EMAX=72.D0 DE=2.D0*AT/(AP+AT)/6.D0 READ(10,*)EMIN,EMAX,DE RMAX=50.D0 DR=0.05D0 READ(10,*)RMAX,DR C========================================== end of input phase IF(IVIBROTT.NE.-1.AND.NTARG.EQ.0.AND.NPHONONT2.NE.0) THEN WRITE(6,91) 91 FORMAT(1H ,'The target excitation should be input' & 1H 'in the first mode :^(') STOP ENDIF IF(IVIBROTT.EQ.1.AND.IVIBROTP.EQ.1) THEN WRITE(6,92) 92 FORMAT(1H ,'This code cannot handle fusion between' & 1H 'two deformed nuclei :^(') STOP ENDIF IF(IVIBROTT.EQ.1.AND.NPHONONT2.NE.0) THEN WRITE(6,*)'!!!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!!!' WRITE(6,*)' ' WRITE(6,*)'The vibrational phonon in the target is NOT' WRITE(6,*)'a deformed phonon, but a *spherical* phonon,' WRITE(6,*)'i.e. the rotation-vibration model is NOT used.' WRITE(6,*)' ' WRITE(6,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' WRITE(6,*)' ' ENDIF CALL NUCLEUS NLEVEL=NLEVEL+NTRANS IF(NLEVEL.GT.NLEVELMAX) THEN WRITE(6,*)'Too many channels :^(' STOP ENDIF 100 WRITE(6,*) 'NLEVEL=',NLEVEL C Evaluation coupling matrix elements C------------------------------------------------------ L=0 C search a minimum of the Coulomb pocket RMIN CALL POTSHAPE ITERAT=INT((RMAX-RMIN)/DR) IF(ITERAT.GT.2999) THEN WRITE(6,*)'Too large grid !! :^(' STOP ENDIF RMAX=RMIN+DR*ITERAT CALL CMAT0 DO 29 IR=0,ITERAT+1 R=RMIN+DR*IR CALL CMAT(R,CPOT0) DO 31 I=1,NLEVEL DO 32 J=1,NLEVEL CPOT(I,J,IR)=CPOT0(I,J) 32 CONTINUE 31 CONTINUE 29 CONTINUE RH=RMIN+DR/2.D0 CALL CMAT(RH,CPOTH) WRITE(9,99) WRITE(9,97) 99 FORMAT(/,' Ecm (MeV) sigma (mb) ') 97 FORMAT(' -------------------------------------') C start the energy roop C------------------------------------------------------------- IF(DE.EQ.0.D0) THEN IESTEP=0 ELSE IESTEP=(EMAX-EMIN)/DE ENDIF IF(IESTEP.EQ.0) THEN OPEN(11,FILE='partial.dat',STATUS='UNKNOWN') ENDIF DO 10 I=0,IESTEP E=EMIN+DE*I SIGMA=0.D0 SPIN=0.D0 C angular momentum roop C----------------------- DO 20 L=0,200 S0=SIGMA RMIN0=RMIN CALL POTSHAPE RMIN00=RMIN RMIN=RMIN0 IF(V(RMIN00).GT.E.OR.RMIN00.LT.0.D0) THEN P=0.D0 GO TO 15 ENDIF C C integration of the c.c. equations CALL MNUMEROV(P) C IF(L.EQ.0) WRITE(9,*)E,P IF(NLEVEL.GT.5) WRITE(6,*)E,L*1.,P IF(IESTEP.EQ.0) THEN WRITE(11,*)E,L*1.,(2.D0*L+1.D0)*P*PI*HBAR**2/2.D0/RMASS/E*10.D0 ENDIF SIGMA=SIGMA+(2.D0*L+1.D0)*P*PI*HBAR**2/2.D0/RMASS/E*10.D0 SPIN=SPIN+L*(2.D0*L+1.D0)*P*PI*HBAR**2/2.D0/RMASS/E*10.D0 IF(SIGMA-S0.LT.S0*1.D-4) GOTO 15 20 CONTINUE 15 IF(SIGMA.NE.0.D0) SPIN=SPIN/SIGMA WRITE(7,*)E,SIGMA WRITE(8,*)E,SPIN IF(SIGMA.LT.1.D-2) THEN WRITE(6,81)E,SIGMA,SPIN WRITE(9,81)E,SIGMA,SPIN ELSE WRITE(6,82)E,SIGMA,SPIN WRITE(9,82)E,SIGMA,SPIN ENDIF 10 CONTINUE 81 FORMAT(F15.5,E15.5,F15.5) 82 FORMAT(3F15.5) WRITE(6,*)'The program terminated normally :^)' STOP END C************************************************************* SUBROUTINE NUCLEUS C C Subroutine to record several input parameters in the C 'OUTPUT' file C C************************************************************* IMPLICIT REAL*8(A-H,O-Z) PARAMETER (NLEVELMAX=25) CHARACTER*1 ANS CHARACTER*2 TEXT,PRO,TARG DIMENSION TEXT(109) COMMON/HION/AP,ZP,RP,AT,ZT,RT COMMON/CONST/HBAR,PI COMMON/INTR/IVIBROTP,IVIBROTT,NPROJ,NTARG COMMON/OSCP/BETAP,BETAPN,OMEGAP,LAMBDAP,NPHONONP COMMON/ROTP/BETA2P,BETA4P,E2P,NROTP COMMON/OSCT/BETAT,BETATN,OMEGAT,LAMBDAT,NPHONONT COMMON/OSCT2/BETAT2,BETAT2N,OMEGAT2,LAMBDAT2,NPHONONT2 COMMON/ROTT/BETA2T,BETA4T,E2T,NROTT COMMON/PHONON/NPHONON COMMON/DIMEN/NLEVEL COMMON/POT/V0,R0,A0 COMMON/DYN/E,RMASS COMMON/SHAPE/RB,VB,CURV COMMON/GRID/RMIN,RMAX,DR COMMON/TRANS/FTR,QTRANS,NTRANS COMMON/MUT/IMUTUAL(0:NLEVELMAX,0:NLEVELMAX) DATA TEXT/'H ','He','Li','Be','B ','C ','N ','O ','F ','Ne', &'Na','Mg','Al','Si','P ','S ','Cl','Ar','K ','Ca','Sc','Ti','V ', &'Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr', &'Rb','Sr','Y ','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In', &'Sn','Sb','Te','I ','Xe','Cs','Ba','La','Ce','Pr','Nd','Pm','Sm', &'Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb','Lu','Hf','Ta','W ','Re', &'Os','Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn','Fr','Ra', &'Ac','Th','Pa','U ','Np','Pu','Am','Cm','Bk','Cf','Es','Fm','Md', &'No','Lr','XX','X1','X2','X3','X4','04'/ NLEVEL=(NPROJ+1)*(NTARG+1) NP=ZP NPP=AP PRO=TEXT(NP) NT=ZT NTT=AT TARG=TEXT(NT) WRITE(9,*)NPP,PRO,' + ',NTT,TARG, ' Fusion reaction' WRITE(6,*)NPP,PRO,' + ',NTT,TARG, ' Fusion reaction' WRITE(9,30) WRITE(6,30) R0P=RP/AP**(1.D0/3.D0) R0T=RT/AT**(1.D0/3.D0) IF(NTARG.NE.0) THEN IF(IVIBROTT.EQ.0) THEN WRITE(6,21)BETAT,OMEGAT,LAMBDAT,NPHONONT 21 FORMAT('Phonon Excitation in the targ.: beta=',F6.3, & ', omega=',F5.2, '(MeV), Lambda=',I2,', Nph=',I2) BETATN=BETAT WRITE(6,*)' ' WRITE(6,*)' Different beta_N from beta_C for this mode(n/y)?' READ(5,1)ANS IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN WRITE(6,*)'beta_N=?' READ(5,*)BETATN ENDIF WRITE(9,211)BETATN,BETAT,r0T,OMEGAT,LAMBDAT,NPHONONT 211 FORMAT('Phonon Excitation in the targ.: beta_N=',F6.3, & ', beta_C=',F6.3,', r0=',F5.2,'(fm)', & 40X, 'omega=',F5.2, '(MeV), Lambda=',I2,', Nph=',I2) ENDIF IF(IVIBROTT.EQ.1) THEN WRITE(9,22)BETA2T,BETA4T,R0T,E2T,NROTT WRITE(6,22)BETA2T,BETA4T,R0T,E2T,NROTT 22 FORMAT('Rotational Excitation in the targ.: beta2=',F6.3, & ', beta4=',F6.3,', r0=',F5.2,'(fm)', & 42X,'E2=',F5.2,'(MeV), Nrot=',I2) ENDIF ENDIF IF(NPHONONT2.NE.0) THEN WRITE(6,23)BETAT2,OMEGAT2,LAMBDAT2,NPHONONT2 23 FORMAT('Phonon Excitation in the targ.: beta=',F6.3, & ', omega=',F5.2, '(MeV), Lambda=',I2,', Nph=',I2) BETAT2N=BETAT2 WRITE(6,*)' ' WRITE(6,*)' Different beta_N from beta_C for this mode(n/y)?' READ(5,1)ANS IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN WRITE(6,*)'beta_N=?' READ(5,*)BETAT2N ENDIF WRITE(9,233)BETAT2N,BETAT2,R0T,OMEGAT2,LAMBDAT2,NPHONONT2 233 FORMAT('Phonon Excitation in the targ.: beta_N=',F6.3, & ', beta_C=',F6.3,', r0=',F5.2,'(fm)', & 40X, 'omega=',F5.2, '(MeV), Lambda=',I2,', Nph=',I2) CALL MUTUAL ENDIF IF(NPROJ.NE.0) THEN IF(IVIBROTP.EQ.0) THEN WRITE(6,24)BETAP,OMEGAP,LAMBDAP,NPHONONP 24 FORMAT('Phonon Excitation in the proj.: beta=',F6.3, & ', omega=',F5.2, '(MeV), Lambda=',I2,', Nph=',I2) BETAPN=BETAP WRITE(6,*)' ' WRITE(6,*)' Different beta_N from beta_C for this mode(n/y)?' READ(5,1)ANS IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN WRITE(6,*)'beta_N=?' READ(5,*)BETAPN ENDIF WRITE(9,244)BETAPN,BETAP,R0P,OMEGAP,LAMBDAP,NPHONONP 244 FORMAT('Phonon Excitation in the proj.: beta_N=',F6.3, & ', beta_C=',F6.3,', r0=',F5.2,'(fm)', & 40X, 'omega=',F5.2, '(MeV), Lambda=',I2,', Nph=',I2) ENDIF IF(IVIBROTP.EQ.1) THEN WRITE(9,25)BETA2P,BETA4P,R0P,E2P,NROTP WRITE(6,25)BETA2P,BETA4P,R0P,E2P,NROTP 25 FORMAT('Rotational Excitation in the proj.: beta2=',F6.3, & ', beta4=',F6.3,', r0=',F5.2,'(fm)', & 42X,'E2=',F5.2,'(MeV), Nrot=',I2) ENDIF ENDIF IF(NTRANS.NE.0) THEN WRITE(9,30) WRITE(9,26)FTR, QTRANS WRITE(6,26)FTR, QTRANS 26 FORMAT('Transfer channel: Strength=',F5.2, & ', Q=',F5.2, '(MeV)') ENDIF WRITE(9,30) R00=R0/(AP**(1.D0/3.D0)+AT**(1.D0/3.D0)) WRITE(9,27)V0,R00,A0 27 FORMAT('Potential parameters: V0=', F8.2, '(MeV), r0=', & F5.2, '(fm), a=', F5.2, '(fm)') CALL POTSHAPE WRITE(9,28)RB,VB,CURV 28 FORMAT(' Uncoupled barrier: Rb=', F5.2, '(fm), Vb=',F8.2, & '(MeV), Curv=',F5.2, '(MeV)') WRITE(9,30) 30 FORMAT('-------------------------------------------------') 1 FORMAT(A1) RETURN END C***************************************************************************** SUBROUTINE MUTUAL C C Subroutine to select levels to be included in the c.c. calculations C C***************************************************************************** IMPLICIT REAL*8(A-H,M,O-Z) PARAMETER (NLEVELMAX=25) CHARACTER*1 ANS,SIGN1,SIGN2 COMMON/INTR/IVIBROTP,IVIBROTT,NPROJ,NTARG COMMON/OSCT/BETAT,BETATN,OMEGAT,LAMBDAT,NPHONONT COMMON/OSCT2/BETAT2,BETAT2N,OMEGAT2,LAMBDAT2,NPHONONT2 COMMON/ROTT/BETA2T,BETA4T,E2T,NROTT COMMON/DIMEN/NLEVEL COMMON/TRANS/FTR,QTRANS,NTRANS COMMON/MUT/IMUTUAL(0:NLEVELMAX,0:NLEVELMAX) DO 10 I=0,NLEVELMAX DO 20 J=0,NLEVELMAX IMUTUAL(I,J)=0 20 CONTINUE 10 CONTINUE WRITE(6,*)' ' WRITE(6,*)'Mutual excitations in the *target* nucleus' WRITE(6,*)' Include the mutual excitations (y/n)?' READ(5,1)ANS 1 FORMAT(A1) IF(ANS.EQ.'n' .OR. ANS.EQ.'N') THEN IMUT=0 NLEVEL=(NTARG+NPHONONT2+1)*(NPROJ+1) c WRITE(9,99) c 99 FORMAT(1H ,'No mutual excitations in the target are included.') DO 80 I=0,NTARG DO 81 J=0,NPHONONT2 IF(I.EQ.0.OR.J.EQ.0) IMUTUAL(I,J)=1 81 CONTINUE 80 CONTINUE GOTO 70 ENDIF WRITE(6,*)' All the possible mutual excitation channels (n/y)?' READ(5,1)ANS IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN IMUT=1 NLEVEL=(NTARG+1)*(NPHONONT2+1)*(NPROJ+1) c WRITE(9,90) c 90 FORMAT(1H ,'All the possible mutual excitation channels' c & 1H ,'in the target are included.') DO 82 I=0,NTARG DO 83 J=0,NPHONONT2 IMUTUAL(I,J)=1 83 CONTINUE 82 CONTINUE GOTO 70 ENDIF IMUT=2 IF((-1.D0)**LAMBDAT.EQ.1.D0) SIGN1='+' IF((-1.D0)**LAMBDAT.EQ.-1.D0) SIGN1='-' IF((-1.D0)**LAMBDAT2.EQ.1.D0) SIGN2='+' IF((-1.D0)**LAMBDAT2.EQ.-1.D0) SIGN2='-' DO 30 I=0,NTARG DO 40 J=0,NPHONONT2 IF(I.EQ.0.OR.J.EQ.0) THEN IMUTUAL(I,J)=1 GOTO 40 ENDIF WRITE(6,94)' Include (',LAMBDAT,SIGN1,'^',I,',',LAMBDAT2,SIGN2, & '^',J,') state ? (y/n)' READ(5,1)ANS IF(ANS.EQ.'n' .OR. ANS.EQ.'N') THEN IMUTUAL(I,J)=0 GOTO 40 ENDIF IMUTUAL(I,J)=1 40 CONTINUE 30 CONTINUE 70 WRITE(9,*)' ' WRITE(9,93) WRITE(6,93) 93 FORMAT(1H ,' Excited states in the target to be included: ') NLEVEL=0 IF((-1.D0)**LAMBDAT.EQ.1.D0) SIGN1='+' IF((-1.D0)**LAMBDAT.EQ.-1.D0) SIGN1='-' IF((-1.D0)**LAMBDAT2.EQ.1.D0) SIGN2='+' IF((-1.D0)**LAMBDAT2.EQ.-1.D0) SIGN2='-' DO 50 I=0,NTARG DO 60 J=0,NPHONONT2 IF(IMUTUAL(I,J).EQ.0) GOTO 60 WRITE(9,94)' (',LAMBDAT,SIGN1,'^',I,',',LAMBDAT2,SIGN2, & '^',J,') state' WRITE(6,94)' (',LAMBDAT,SIGN1,'^',I,',',LAMBDAT2,SIGN2, & '^',J,') state' 94 FORMAT(1H ,A,I2,2A,I2,A,I2,2A,I2,A) NLEVEL=NLEVEL+1 60 CONTINUE 50 CONTINUE NLEVEL=NLEVEL*(NPROJ+1) RETURN END C***************************************************************************** SUBROUTINE POTSHAPE C C shape of the bare Coulomb barrier, i.e. the barrier position, the curvature, C and the position of the Coulomb pocket C C***************************************************************************** IMPLICIT REAL*8(A-H,M,O-Z) COMMON/SHAPE/RB,VB,CURV COMMON/GRID/RMIN,RMAX,DR COMMON/DYN/E,RMASS COMMON/CONST/HBAR,PI EXTERNAL DV,V R=50.5D0 U0=DV(R) 10 R=R-1.D0 IF(R.LT.0.D0) THEN RB=-5.D0 RMIN=-5.D0 RETURN ENDIF U1=DV(R) IF(U0*U1.GT.0.D0) GOTO 10 RA=R+1.D0 RB=R TOLK=1.D-6 N=LOG10(ABS(RB-RA)/TOLK)/LOG10(2.D0)+0.5D0 DO 20 I=1,N R=(RA+RB)/2.D0 U=DV(R) IF(U0*U.LT.0.D0) THEN RB=R ELSE RA=R ENDIF 20 CONTINUE RB=R DDV00=(DV(RB+1.D-5)-DV(RB-1.D-5))/2.D-5 IF(DDV00.GT.0.D0) THEN WRITE(6,*)'SOMETHING IS STRANGE :^(' WRITE(8,*)'SOMETHING IS STRANGE :^(' STOP ENDIF CURV=HBAR*SQRT(ABS(DDV00)/RMASS) VB=V(RB) RA=RB-0.5D0 U0=DV(RA) RBB=0.5D0 U1=DV(RBB) TOLK=1.D-6 N=LOG10(ABS(RB-RA)/TOLK)/LOG10(2.D0)+0.5D0 DO 21 I=1,N R=(RA+RBB)/2.D0 U=DV(R) IF(U0*U.LT.0.D0) THEN RBB=R ELSE RA=R ENDIF 21 CONTINUE RMIN=R C WRITE(6,*)'RMIN=',RMIN,V(RMIN) C WRITE(6,*)'RB=',RB,V(RB) C WRITE(6,*)'CURV=',CURV00 RETURN END C************************************************************** SUBROUTINE MNUMEROV(P) C C Subroutine for integration of the c.c. eqs. by modified C Numerov method C C************************************************************** IMPLICIT REAL*8 (A-H,O-Z) EXTERNAL V PARAMETER (NLEVELMAX=25) DIMENSION PSI(NLEVELMAX),PSI0(NLEVELMAX),PSI1(NLEVELMAX) DIMENSION XI(NLEVELMAX),XI0(NLEVELMAX),XI1(NLEVELMAX) DIMENSION PHI0(NLEVELMAX) DIMENSION BB(NLEVELMAX,NLEVELMAX), & BIN(NLEVELMAX,NLEVELMAX) DIMENSION CC(NLEVELMAX,NLEVELMAX), & CIN(NLEVELMAX,NLEVELMAX) DIMENSION DD0(NLEVELMAX,NLEVELMAX),DD1(NLEVELMAX,NLEVELMAX) DIMENSION DD(NLEVELMAX) DIMENSION FCW(0:200),GCW(0:200),FPCW(0:200),GPCW(0:200) DIMENSION SIGMAD(0:200),IEXP(0:200) DIMENSION ECH(NLEVELMAX) COMMON/HION/AP,ZP,RP,AT,ZT,RT COMMON/CONST/HBAR,PI COMMON/DIMEN/NLEVEL COMMON/DYN/E,RMASS COMMON/ANGULAR/L COMMON/GRID/RMIN,RMAX,DR COMMON/COUP/CPOT(NLEVELMAX,NLEVELMAX,0:3000) COMPLEX*16 K,KK COMPLEX*16 PSI,PSI0,PSI1,PHI0 COMPLEX*16 XI,XI0,XI1 COMPLEX*16 AI COMPLEX*16 CWUP0,CWDOWN0,CWUP1,CWDOWN1 COMPLEX*16 BB,BIN,CC,CIN,DD AI=(0.D0,1.D0) C define two constants used in this subroutine FAC=DR**2*(2.D0*RMASS/HBAR**2) ITERAT=INT((RMAX-RMIN)/DR) RMAX=RMIN+ITERAT*DR DO 22 LC=0,200 FCW(LC)=0.D0 GCW(LC)=0.D0 FPCW(LC)=0.D0 GPCW(LC)=0.D0 SIGMAD(LC)=0.D0 IEXP(LC)=0D0 22 CONTINUE DO 51 I=1,NLEVELMAX DO 52 J=1,NLEVELMAX BB(I,J)=0.D0 BIN(I,J)=0.D0 CC(I,J)=0.D0 CIN(I,J)=0.D0 52 CONTINUE DD(I)=0.D0 51 CONTINUE C integration of the io-th channel wave function from RMIN C========================================================== DO 15 IO=1,NLEVEL c write(6,*)'IO=',IO*1. DO 200 J1=1,NLEVEL PSI(J1)=0.D0 PSI0(J1)=0.D0 PSI1(J1)=0.D0 PHI0(J1)=0.D0 XI0(J1)=0.D0 XI1(J1)=0.D0 XI(J1)=0.D0 200 CONTINUE C initial conditions C---------- ECH(IO)=E-V(RMIN)-CPOT(IO,IO,0) C -0.25D0*HBAR**2/2.D0/RMASS/RMIN**2 IF(ECH(IO).GT.0.D0) THEN K=SQRT(2.D0*RMASS/HBAR**2*ECH(IO)) PSI0(IO)=EXP(-AI*K*RMIN) PHI0(IO)=-AI*K*PSI0(IO) ELSE K=SQRT(2.D0*RMASS/HBAR**2*ABS(ECH(IO))) PSI0(IO)=EXP(K*RMIN) PHI0(IO)=K*PSI0(IO) ENDIF CALL RKUTTA00(PSI0,PHI0,PSI1) DO 91 I0=1,NLEVEL XI0(I0)=(1.D0-FAC/12.D0*(V(RMIN)-E))*PSI0(I0) XI1(I0)=(1.D0-FAC/12.D0*(V(RMIN+DR)-E))*PSI1(I0) DO 92 IC=1,NLEVEL XI0(I0)=XI0(I0)-FAC/12.D0*CPOT(I0,IC,0)*PSI0(IC) XI1(I0)=XI1(I0)-FAC/12.D0*CPOT(I0,IC,1)*PSI1(IC) 92 CONTINUE 91 CONTINUE DO 29 IR=2,ITERAT+1 C---------------------------------------------- ITERATIONS START R=RMIN+DR*IR R0=RMIN+DR*(IR-2) R1=RMIN+DR*(IR-1) DO 71 I0=1,NLEVEL DO 72 IC=1,NLEVEL DD0(I0,IC)=FAC/SQRT(12.D0)*CPOT(I0,IC,IR-1) IF(I0.EQ.IC) DD0(I0,IC)=DD0(I0,IC)+FAC/SQRT(12.D0)*(V(R1)-E) & +SQRT(3.D0) 72 CONTINUE 71 CONTINUE DO 73 I0=1,NLEVEL DO 74 IC=1,NLEVEL DD1(I0,IC)=0.D0 IF(I0.EQ.IC) DD1(I0,IC)=DD1(I0,IC)-1.D0 DO 75 IK=1,NLEVEL DD1(I0,IC)=DD1(I0,IC)+DD0(I0,IK)*DD0(IK,IC) 75 CONTINUE 74 CONTINUE 73 CONTINUE DO 54 I0=1,NLEVEL XI(I0)=-XI0(I0) DO 53 IC=1,NLEVEL XI(I0)=XI(I0)+DD1(I0,IC)*XI1(IC) 53 CONTINUE 54 CONTINUE IF(IR.EQ.ITERAT+1) GOTO 66 DO 42 I0=1,NLEVEL XI0(I0)=XI1(I0) XI1(I0)=XI(I0) 42 CONTINUE 29 CONTINUE C-------------------------------------------------------------- C Matching to the Coulomb wave function at RMAX 66 DO 56 I0=1,NLEVEL DO 55 IC=1,NLEVEL CC(I0,IC)=-FAC/12.D0*CPOT(I0,IC,ITERAT-1) IF(I0.EQ.IC) CC(I0,IC)=CC(I0,IC)-FAC/12.D0*(V(RMAX-DR)-E)+1.D0 55 CONTINUE 56 CONTINUE CALL MATINV(NLEVEL,CC,CIN) DO 94 I0=1,NLEVEL PSI0(I0)=0.D0 DO 93 IC=1,NLEVEL PSI0(I0)=PSI0(I0)+CIN(I0,IC)*XI0(IC) 93 CONTINUE 94 CONTINUE DO 556 I0=1,NLEVEL DO 555 IC=1,NLEVEL CC(I0,IC)=-FAC/12.D0*CPOT(I0,IC,ITERAT+1) IF(I0.EQ.IC) CC(I0,IC)=CC(I0,IC)-FAC/12.D0*(V(RMAX+DR)-E)+1.D0 555 CONTINUE 556 CONTINUE CALL MATINV(NLEVEL,CC,CIN) DO 944 I0=1,NLEVEL PSI(I0)=0.D0 DO 933 IC=1,NLEVEL PSI(I0)=PSI(I0)+CIN(I0,IC)*XI(IC) 933 CONTINUE 944 CONTINUE DO 60 II=1,NLEVEL C Coulomb wave function EC=E-CPOT(II,II,ITERAT-1) RHO=SQRT(2.D0*RMASS*EC)/HBAR*(RMAX-DR) ETA=(ZP*ZT/137.D0)*SQRT(RMASS/2.D0/EC) CALL DFCOUL(ETA,RHO,FCW,FPCW,GCW,GPCW,SIGMAD,L,IEXP) CWUP0=(GCW(L)+AI*FCW(L)) CWDOWN0=GCW(L)-AI*FCW(L) EC=E-CPOT(II,II,ITERAT+1) RHO=SQRT(2.D0*RMASS*EC)/HBAR*(RMAX+DR) ETA=(ZP*ZT/137.D0)*SQRT(RMASS/2.D0/EC) CALL DFCOUL(ETA,RHO,FCW,FPCW,GCW,GPCW,SIGMAD,L,IEXP) CWUP1=(GCW(L)+AI*FCW(L)) CWDOWN1=GCW(L)-AI*FCW(L) BB(IO,II)=(CWUP0*PSI(II)-CWUP1*PSI0(II)) & /(CWUP0*CWDOWN1-CWUP1*CWDOWN0) 60 CONTINUE C--------------------------------------------------------------- 15 CONTINUE C=============================================================== C PENETRATION PROBABILITY C calculate the inversion of the matrix BB CALL MATINV(NLEVEL,BB,BIN) C P=0.D0 DO 85 IO=1,NLEVEL IF(ECH(IO).LT.0.D0) GOTO 85 K=SQRT((2.D0*RMASS/HBAR**2*ECH(IO))) KK=SQRT(2.D0*RMASS/HBAR**2*E) P=P+(ABS(BIN(1,IO)))**2*K/KK 85 CONTINUE IF(P.GT.1.1D0) THEN WRITE(6,*)'Penetrability > 1 !!!!!' WRITE(6,*)'Something must be strange :^(' STOP ENDIF RETURN END C************************************************************** SUBROUTINE RKUTTA00(PSI0,PHI0,PSI1) C C Wave functions at R=RMIN+DR C C************************************************************** IMPLICIT REAL*8 (A-H,O-Z) EXTERNAL V PARAMETER (NLEVELMAX=25) DIMENSION PSI0(NLEVELMAX),PHI0(NLEVELMAX),PSI1(NLEVELMAX) DIMENSION AK1(NLEVELMAX),AK2(NLEVELMAX), & AK3(NLEVELMAX),AK4(NLEVELMAX) DIMENSION BK1(NLEVELMAX),BK2(NLEVELMAX), & BK3(NLEVELMAX),BK4(NLEVELMAX) COMMON/HION/AP,ZP,RP,AT,ZT,RT COMMON/CONST/HBAR,PI COMMON/DIMEN/NLEVEL COMMON/DYN/E,RMASS COMMON/ANGULAR/L COMMON/GRID/RMIN,RMAX,DR COMMON/COUP/CPOT(NLEVELMAX,NLEVELMAX,0:3000) COMMON/COUPH/CPOTH(NLEVELMAX,NLEVELMAX) COMPLEX*16 PSI0,PHI0,PSI1 COMPLEX*16 AI COMPLEX*16 AK1,AK2,AK3,AK4 COMPLEX*16 BK1,BK2,BK3,BK4 AI=(0.D0,1.D0) C define two constants used in this subroutine FAC=DR*(2.D0*RMASS/HBAR**2) R=RMIN RPP=RMIN+DR RH=RMIN+DR/2.D0 DO 19 J1=1,NLEVEL AK1(J1)=0.D0 AK2(J1)=0.D0 AK3(J1)=0.D0 AK4(J1)=0.D0 BK1(J1)=0.D0 BK2(J1)=0.D0 BK3(J1)=0.D0 BK4(J1)=0.D0 19 CONTINUE DO 58 I0=1,NLEVEL DO 57 IC=1,NLEVEL AK1(I0)=AK1(I0)+FAC*CPOT(I0,IC,0)*PSI0(IC) 57 CONTINUE AK1(I0)=AK1(I0)-FAC*(E-V(R))*PSI0(I0) BK1(I0)=DR*PHI0(I0) 58 CONTINUE DO 56 I0=1,NLEVEL DO 55 IC=1,NLEVEL AK2(I0)=AK2(I0)+FAC*CPOTH(I0,IC)*(PSI0(IC)+1.D0/2.D0*BK1(IC)) 55 CONTINUE AK2(I0)=AK2(I0)-FAC*(E-V(RH))*(PSI0(I0)+1.D0/2.D0*BK1(I0)) BK2(I0)=DR*(PHI0(I0)+1./2.*AK1(I0)) 56 CONTINUE DO 54 I0=1,NLEVEL DO 53 IC=1,NLEVEL AK3(I0)=AK3(I0)+FAC*CPOTH(I0,IC)*(PSI0(IC)+1.D0/2.D0*BK2(IC)) 53 CONTINUE AK3(I0)=AK3(I0)-FAC*(E-V(RH))*(PSI0(I0)+1.D0/2.D0*BK2(I0)) BK3(I0)=DR*(PHI0(I0)+1.D0/2.D0*AK2(I0)) 54 CONTINUE DO 42 I0=1,NLEVEL DO 41 IC=1,NLEVEL AK4(I0)=AK4(I0)+FAC*CPOT(I0,IC,1)*(PSI0(IC)+BK3(IC)) 41 CONTINUE AK4(I0)=AK4(I0)-FAC*(E-V(RPP))*(PSI0(I0)+BK3(I0)) BK4(I0)=DR*(PHI0(I0)+AK3(I0)) 42 CONTINUE C wave functions at RMIN+DR DO 59 IS=1,NLEVEL PSI1(IS)=PSI0(IS)+1.D0/6.D0*(BK1(IS)+2.D0*BK2(IS) & +2.D0*BK3(IS)+BK4(IS)) 59 CONTINUE RETURN END C****************************************************************** SUBROUTINE CMAT0 C C Preparation for the nuclear coupling matrix C C***************************************************************** IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NLEVELMAX=25) DIMENSION A(NLEVELMAX,NLEVELMAX) COMMON/HION/AP,ZP,RP,AT,ZT,RT COMMON/CONST/HBAR,PI COMMON/INTR/IVIBROTP,IVIBROTT,NPROJ,NTARG COMMON/OSCP/BETAP,BETAPN,OMEGAP,LAMBDAP,NPHONONP COMMON/OSCT/BETAT,BETATN,OMEGAT,LAMBDAT,NPHONONT COMMON/OSCT2/BETAT2,BETAT2N,OMEGAT2,LAMBDAT2,NPHONONT2 COMMON/ROTP/BETA2P,BETA4P,E2P,NROTP COMMON/ROTT/BETA2T,BETA4T,E2T,NROTT COMMON/DIMEN/NLEVEL COMMON/MUT/IMUTUAL(0:NLEVELMAX,0:NLEVELMAX) COMMON/EIGEN/EV(NLEVELMAX),EVEC(NLEVELMAX,NLEVELMAX) COMMON/TRANS/FTR,QTRANS,NTRANS EXTERNAL CG FVIBT=RT*BETATN/SQRT(4.D0*PI) FROT2T=RT*BETA2T FROT4T=RT*BETA4T FVIBT2=RT*BETAT2N/SQRT(4.D0*PI) FVIBP=RP*BETAPN/SQRT(4.D0*PI) FROT2P=RP*BETA2P FROT4P=RP*BETA4P DO 10 I=1,NLEVELMAX DO 20 J=1,NLEVELMAX A(I,J)=0.D0 20 CONTINUE 10 CONTINUE NDIM=NLEVEL-NTRANS IF(NDIM.EQ.1) RETURN I=0 DO 31 IP=0,NPROJ DO 32 IT=0,NTARG DO 33 IT2=0,NPHONONT2 IF(NPHONONT2.NE.0) THEN IF(IMUTUAL(IT,IT2).EQ.0) GOTO 33 ENDIF I=I+1 J=0 DO 41 JP=0,NPROJ DO 42 JT=0,NTARG DO 43 JT2=0,NPHONONT2 IF(NPHONONT2.NE.0) THEN IF(IMUTUAL(JT,JT2).EQ.0) GOTO 43 ENDIF J=J+1 IF(I.GT.J) THEN A(I,J)=A(J,I) GOTO 43 ENDIF C=0.D0 IF(IP.EQ.JP.AND.IT2.EQ.JT2) THEN IF(IVIBROTT.EQ.0) THEN IF(IT.EQ.JT-1) C=SQRT(JT*1.D0)*FVIBT IF(JT.EQ.IT-1) C=SQRT(IT*1.D0)*FVIBT ELSE C=FROT2T*SQRT((2.D0*2*IT+1.D0)*5.D0*(2.D0*2*JT+1.D0) & /4.D0/PI)*CG(2*IT,0,2,0,2*JT,0)**2/(2.D0*2*JT+1.D0) C=C+FROT4T*SQRT((2.D0*2*IT+1.D0)*9.D0*(2.D0*2*JT+1.D0) & /4.D0/PI)*CG(2*IT,0,4,0,2*JT,0)**2/(2.D0*2*JT+1.D0) ENDIF ENDIF IF(IP.EQ.JP.AND.IT.EQ.JT) THEN IF(IT2.EQ.JT2-1) C=C+SQRT(JT2*1.D0)*FVIBT2 IF(JT2.EQ.IT2-1) C=C+SQRT(IT2*1.D0)*FVIBT2 ENDIF IF(IT.EQ.JT.AND.IT2.EQ.JT2) THEN IF(IVIBROTP.EQ.0) THEN IF(IP.EQ.JP-1) C=C+SQRT(JP*1.D0)*FVIBP IF(JP.EQ.IP-1) C=C+SQRT(IP*1.D0)*FVIBP ELSE C=C+FROT2P*SQRT((2.D0*2*IP+1.D0)*5.D0*(2.D0*2*JP+1.D0) & /4.D0/PI)*CG(2*IP,0,2,0,2*JP,0)**2/(2.D0*2*JP+1.D0) C=C+FROT4P*SQRT((2.D0*2*IP+1.D0)*9.D0*(2.D0*2*JP+1.D0) & /4.D0/PI)*CG(2*IP,0,4,0,2*JP,0)**2/(2.D0*2*JP+1.D0) ENDIF ENDIF A(I,J)=C 43 CONTINUE 42 CONTINUE 41 CONTINUE 33 CONTINUE 32 CONTINUE 31 CONTINUE CALL MDIAG(A,NDIM) RETURN END C****************************************************************** SUBROUTINE CMAT(R,CPOT) C C Coupling matrix C C***************************************************************** IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NLEVELMAX=25) DIMENSION CPOT(NLEVELMAX,NLEVELMAX) COMMON/HION/AP,ZP,RP,AT,ZT,RT COMMON/CONST/HBAR,PI COMMON/INTR/IVIBROTP,IVIBROTT,NPROJ,NTARG COMMON/OSCP/BETAP,BETAPN,OMEGAP,LAMBDAP,NPHONONP COMMON/OSCT/BETAT,BETATN,OMEGAT,LAMBDAT,NPHONONT COMMON/OSCT2/BETAT2,BETAT2N,OMEGAT2,LAMBDAT2,NPHONONT2 COMMON/ROTP/BETA2P,BETA4P,E2P,NROTP COMMON/ROTT/BETA2T,BETA4T,E2T,NROTT COMMON/DIMEN/NLEVEL COMMON/MUT/IMUTUAL(0:NLEVELMAX,0:NLEVELMAX) COMMON/EIGEN/EV(NLEVELMAX),EVEC(NLEVELMAX,NLEVELMAX) COMMON/TRANS/FTR,QTRANS,NTRANS EXTERNAL FCT,FCP,FACT,VNCC,CG,FCT2,FCT4,FCP2,FCP4,FCTT EXTERNAL FTRANS DO 10 I=1,NLEVELMAX DO 20 J=1,NLEVELMAX CPOT(I,J)=0.D0 20 CONTINUE 10 CONTINUE NDIM=NLEVEL-NTRANS IF(NLEVEL.EQ.1) RETURN IF(NDIM.EQ.1) GOTO 55 I=0 DO 31 IP=0,NPROJ DO 32 IT=0,NTARG DO 33 IT2=0,NPHONONT2 IF(NPHONONT2.NE.0) THEN IF(IMUTUAL(IT,IT2).EQ.0) GOTO 33 ENDIF I=I+1 J=0 DO 41 JP=0,NPROJ DO 42 JT=0,NTARG DO 43 JT2=0,NPHONONT2 IF(NPHONONT2.NE.0) THEN IF(IMUTUAL(JT,JT2).EQ.0) GOTO 43 ENDIF J=J+1 IF(I.GT.J) THEN CPOT(I,J)=CPOT(J,I) GOTO 43 ENDIF C=0.D0 C Nuclear coupling DO 50 K=1,NDIM C=C+VNCC(R,EV(K))*EVEC(I,K)*EVEC(J,K) 50 CONTINUE C Coulomb coupling IF(IP.EQ.JP.AND.IT2.EQ.JT2) THEN IF(IVIBROTT.EQ.0) THEN IF(IT.EQ.JT-1) C=C+SQRT(JT*1.D0)*FCT(R) IF(JT.EQ.IT-1) C=C+SQRT(IT*1.D0)*FCT(R) ELSE C=C+SQRT((2.D0*2*IT+1.D0)*5.D0*(2.D0*2*JT+1.D0) & /4.D0/PI)*CG(2*IT,0,2,0,2*JT,0)**2/(2.D0*2*JT+1.D0)*FCT2(R) & +SQRT((2.D0*2*IT+1.D0)*9.D0*(2.D0*2*JT+1.D0) & /4.D0/PI)*CG(2*IT,0,4,0,2*JT,0)**2/(2.D0*2*JT+1.D0)*FCT4(R) ENDIF ENDIF IF(IT.EQ.JT.AND.IT2.EQ.JT2) THEN IF(IVIBROTP.EQ.0) THEN IF(IP.EQ.JP-1) C=C+SQRT(JP*1.D0)*FCP(R) IF(JP.EQ.IP-1) C=C+SQRT(IP*1.D0)*FCP(R) ELSE C=C+SQRT((2.D0*2*IP+1.D0)*5.D0*(2.D0*2*JP+1.D0) & /4.D0/PI)*CG(2*IP,0,2,0,2*JP,0)**2/(2.D0*2*JP+1.D0)*FCP2(R) & +SQRT((2.D0*2*IP+1.D0)*9.D0*(2.D0*2*JP+1.D0) & /4.D0/PI)*CG(2*IP,0,4,0,2*JP,0)**2/(2.D0*2*JP+1.D0)*FCP4(R) ENDIF ENDIF IF(IP.EQ.JP.AND.IT.EQ.JT) THEN IF(IT2.EQ.JT2-1) C=C+SQRT(JT2*1.D0)*FCTT(R) IF(JT2.EQ.IT2-1) C=C+SQRT(IT2*1.D0)*FCTT(R) ENDIF C Excitation Energy IF(IT.EQ.JT.AND.IP.EQ.JP.AND.IT2.EQ.JT2) THEN IF(IVIBROTT.EQ.0) THEN C=C+IT*OMEGAT ELSE C=C+2.D0*IT*(2.D0*IT+1.D0)/6.D0*E2T ENDIF IF(IVIBROTP.EQ.0) THEN C=C+IP*OMEGAP ELSE C=C+2.D0*IP*(2.D0*IP+1.D0)/6.D0*E2P ENDIF C=C+IT2*OMEGAT2 ENDIF CPOT(I,J)=C 43 CONTINUE 42 CONTINUE 41 CONTINUE 33 CONTINUE 32 CONTINUE 31 CONTINUE C0=CPOT(1,1) DO 51 I=1,NDIM CPOT(I,I)=CPOT(I,I)-VN(R) c CPOT(I,I)=CPOT(I,I)-C0 51 CONTINUE c transfer channel 55 IF(NTRANS.EQ.1) THEN CPOT(1,NLEVEL)=FTRANS(R) CPOT(NLEVEL,1)=FTRANS(R) CPOT(NLEVEL,NLEVEL)=CPOT(1,1)-QTRANS ENDIF RETURN END C***************************************************************** FUNCTION V(R) C C potential for the relative motion C C***************************************************************** IMPLICIT REAL*8(A-H,O-Z) EXTERNAL VC,VN,VCENT V=VC(R)+VN(R)+VCENT(R) RETURN END C***************************************************************** FUNCTION VC(R) C C Coulomb potential C C***************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/HION/AP,ZP,RP,AT,ZT,RT COMMON/CONST/HBAR,PI C Coulomb radius C RC=1.2D0*(AP**(1.D0/3.D0)+AT**(1.D0/3.D0)) C IF(R.GT.RC) THEN VC=ZP*ZT/R*HBAR/137.D0 C ELSE C VC=ZP*ZT*HBAR/137.D0*(3.D0*RC**2-R**2)/2.D0/RC**3 C ENDIF RETURN END C***************************************************************** FUNCTION DVC(R) C C first derivative of the Coulomb potential C C***************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/HION/AP,ZP,RP,AT,ZT,RT COMMON/CONST/HBAR,PI C Coulomb radius C RC=1.2D0*(AP**(1.D0/3.D0)+AT**(1.D0/3.D0)) C IF(R.GT.RC) THEN DVC=-ZP*ZT/R**2*HBAR/137.D0 C ELSE C VC=ZP*ZT*HBAR/137.D0*(3.D0*RC**2-R**2)/2.D0/RC**3 C ENDIF RETURN END C**************************************************************** FUNCTION VN(R) C C Woods-Saxon potential for nuclear potential C C**************************************************************** IMPLICIT REAL*8 (A-H,O-Z) COMMON/HION/AP,ZP,RP,AT,ZT,RT COMMON/CONST/HBAR,PI COMMON/POT/V0,R0,A0 IF(R.LT.50.D0) THEN VN=-V0/(1.D0+EXP((R-R0)/A0)) ELSE VN=0.D0 ENDIF RETURN END C**************************************************************** FUNCTION VNCC(R,XT) C C Woods-Saxon potential for nuclear potential C C**************************************************************** IMPLICIT REAL*8 (A-H,O-Z) COMMON/HION/AP,ZP,RP,AT,ZT,RT COMMON/CONST/HBAR,PI COMMON/POT/V0,R0,A0 IF(R.LT.50.D0) THEN RR=R-R0-XT VNCC=-V0/(1.D0+EXP(RR/A0)) ELSE VNCC=0.D0 ENDIF RETURN END C*************************************************************** FUNCTION DVN(R) C C The first derivative of VN(R) C C************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/POT/V0,R0,A IF(R.LT.50.D0) THEN DVN=V0/A*EXP((R-R0)/A)/(1.D0+EXP((R-R0)/A))**2 ELSE DVN=0.D0 ENDIF RETURN END C************************************************************ FUNCTION VCENT(R) C C Centrifugal potential C C************************************************************ IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HBAR,PI COMMON/ANGULAR/L COMMON/DYN/E,RMASS VCENT=L*(L+1.D0)*HBAR**2/2.D0/RMASS/R**2 RETURN END C************************************************************ FUNCTION DVCENT(R) C C The first derivative of the centrifugal potential C C************************************************************ IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HBAR,PI COMMON/ANGULAR/L COMMON/DYN/E,RMASS DVCENT=-2.D0*L*(L+1.D0)*HBAR**2/2.D0/RMASS/R**3 RETURN END C*************************************************************** FUNCTION DV(R) C C The first derivative of V(R) C C*************************************************************** IMPLICIT REAL*8(A-H,O-Z) EXTERNAL DVN,DVC,DVCENT DV=DVN(R)+DVC(R)+DVCENT(R) RETURN END C***************************************************************** FUNCTION FCP(R) C C Coulomb coupling form factor for projectile phonon excitation C C***************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/HION/AP,ZP,RP,AT,ZT,RT COMMON/CONST/HBAR,PI COMMON/OSCP/BETAP,BETAPN,OMEGA,LAMBDA,NPHONON RC=RP IF(R.GT.RC) THEN FCP=3.D0/(2.D0*LAMBDA+1.D0)*ZP*ZT/137.D0*HBAR/R & *(RP/R)**LAMBDA ELSE FCP=3.D0/(2.D0*LAMBDA+1.D0)*ZP*ZT/137.D0*HBAR/RC & *(R/RP)**LAMBDA ENDIF FCP=FCP*BETAP/SQRT(4.D0*PI) RETURN END C***************************************************************** FUNCTION FCP2(R) C C Coulomb coupling form factor for projectile excitation C (rotational E2 coupling) C C***************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/HION/AP,ZP,RP,AT,ZT,RT COMMON/CONST/HBAR,PI COMMON/ROTP/BETA2P,BETA4P,E2P,NROTP LAMBDA=2 RC=RP IF(R.GT.RC) THEN FCP2=3.D0/(2.D0*LAMBDA+1.D0)*ZP*ZT/137.D0*HBAR/R & *(RP/R)**LAMBDA ELSE FCP2=3.D0/(2.D0*LAMBDA+1.D0)*ZP*ZT/137.D0*HBAR/RC & *(R/RP)**LAMBDA ENDIF FCP2=FCP2*(BETA2P+2.D0*SQRT(5.D0/PI)*BETA2P**2/7.D0) RETURN END C***************************************************************** FUNCTION FCP4(R) C C Coulomb coupling form factor for projectile excitation C (rotational E4 coupling) C C***************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/HION/AP,ZP,RP,AT,ZT,RT COMMON/CONST/HBAR,PI COMMON/ROTP/BETA2P,BETA4P,E2P,NROTP LAMBDA=4 RC=RP IF(R.GT.RC) THEN FCP4=3.D0/(2.D0*LAMBDA+1.D0)*ZP*ZT/137.D0*HBAR/R & *(RP/R)**LAMBDA ELSE FCP4=3.D0/(2.D0*LAMBDA+1.D0)*ZP*ZT/137.D0*HBAR/RC & *(R/RP)**LAMBDA ENDIF FCP4=FCP4*(BETA4P+9.D0*BETA2P**2/7.D0/SQRT(PI)) RETURN END C***************************************************************** FUNCTION FCT(R) C C Coulomb coupling form factor for the target phonon excitation C C***************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/HION/AP,ZP,RP,AT,ZT,RT COMMON/CONST/HBAR,PI COMMON/OSCT/BETAT,BETATN,OMEGA,LAMBDA,NPHONON RC=RT IF(R.GT.RC) THEN FCT=3.D0/(2.D0*LAMBDA+1.D0)*ZP*ZT/137.D0*HBAR/R & *(RT/R)**LAMBDA ELSE FCT=3.D0/(2.D0*LAMBDA+1.D0)*ZP*ZT/137.D0*HBAR/RC & *(R/RT)**LAMBDA ENDIF FCT=FCT*BETAT/SQRT(4.D0*PI) RETURN END C***************************************************************** FUNCTION FCTT(R) C C Coulomb coupling form factor for the second target phonon excitation C C***************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/HION/AP,ZP,RP,AT,ZT,RT COMMON/CONST/HBAR,PI COMMON/OSCT2/BETAT2,BETAT2N,OMEGA2,LAMBDA2,NPHONON2 RC=RT IF(R.GT.RC) THEN FCTT=3.D0/(2.D0*LAMBDA2+1.D0)*ZP*ZT/137.D0*HBAR/R & *(RT/R)**LAMBDA2 ELSE FCTT=3.D0/(2.D0*LAMBDA2+1.D0)*ZP*ZT/137.D0*HBAR/RC & *(R/RT)**LAMBDA2 ENDIF FCTT=FCTT*BETAT2/SQRT(4.D0*PI) RETURN END C***************************************************************** FUNCTION FCT2(R) C C Coulomb coupling form factor for target excitation C (rotational E2 coupling) C C***************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/HION/AP,ZP,RP,AT,ZT,RT COMMON/CONST/HBAR,PI COMMON/ROTT/BETA2T,BETA4T,E2T,NROTT LAMBDA=2 RC=RT IF(R.GT.RC) THEN FCT2=3.D0/(2.D0*LAMBDA+1.D0)*ZP*ZT/137.D0*HBAR/R & *(RT/R)**LAMBDA ELSE FCT2=3.D0/(2.D0*LAMBDA+1.D0)*ZP*ZT/137.D0*HBAR/RC & *(R/RT)**LAMBDA ENDIF FCT2=FCT2*(BETA2T+2.D0*SQRT(5.D0/PI)*BETA2T**2/7.D0) RETURN END C***************************************************************** FUNCTION FCT4(R) C C Coulomb coupling form factor for target excitation C (rotational E4 coupling) C C***************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/HION/AP,ZP,RP,AT,ZT,RT COMMON/CONST/HBAR,PI COMMON/ROTT/BETA2T,BETA4T,E2T,NROTT LAMBDA=4 RC=RT IF(R.GT.RC) THEN FCT4=3.D0/(2.D0*LAMBDA+1.D0)*ZP*ZT/137.D0*HBAR/R & *(RT/R)**LAMBDA ELSE FCT4=3.D0/(2.D0*LAMBDA+1.D0)*ZP*ZT/137.D0*HBAR/RC & *(R/RT)**LAMBDA ENDIF FCT4=FCT4*(BETA4T+9.D0*BETA2T**2/7.D0/SQRT(PI)) RETURN END C*************************************************************** FUNCTION DDVN(R) C C The second derivative of VN(R) C C************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/POT/V0,R0,A IF(R.LT.50.D0) THEN DDVN=V0/A**2*EXP((R-R0)/A)/(1.D0+EXP((R-R0)/A))**2 & -2.D0*V0/A**2*EXP(2.D0*(R-R0)/A) & /(1.D0+EXP((R-R0)/A))**3 ELSE DDVN=0.D0 ENDIF RETURN END C**************************************************************** FUNCTION FTRANS(R) C C Coupling form factor for transfer reactions C C*************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HBAR,PI COMMON/TRANS/FTR,QTRANS,NTRANS EXTERNAL DVN FTRANS=FTR*DVN(R) RETURN END C************************************************************* SUBROUTINE MDIAG(A,N) C C The Jacobi method to compute all eigenvalues and eigenvectors C of a real symmetric matrix A, which is of size N by N, stored C in a physical NP by NP array. On output, the elements of A C above the diagonal are destroyed. D returns the eigenvalues of C A in its first N elements. V is a matrix with the same logical C and physical dimensions as A whose columns contain, on output, C the normalized eigenvectors of A. NROT returns the number of C Jacobi rotations which were required. C C N.B.; The I-th component of the eigenvector for the K-th C eigenvalue is given by V(I,K). C C************************************************************ IMPLICIT REAL*8(A-H,O-Z) PARAMETER(NMAX=100) PARAMETER (NLEVELMAX=25) DIMENSION A(NLEVELMAX,NLEVELMAX) DIMENSION B(NMAX),Z(NMAX) COMMON/EIGEN/D(NLEVELMAX),V(NLEVELMAX,NLEVELMAX) DO 12 IP=1,N DO 11 IQ=1,N V(IP,IQ)=0.D0 11 CONTINUE V(IP,IP)=1.D0 12 CONTINUE DO 13 IP=1,N B(IP)=A(IP,IP) D(IP)=B(IP) Z(IP)=0.D0 13 CONTINUE NROT=0 DO 24 I=1,50 SM=0.D0 DO 15 IP=1,N-1 DO 14 IQ=IP+1,N SM=SM+ABS(A(IP,IQ)) 14 CONTINUE 15 CONTINUE IF(SM.EQ.0.D0) RETURN IF(I.LT.4) THEN TRESH=0.2D0*SM/N**2 ELSE TRESH=0.D0 ENDIF DO 22 IP=1,N-1 DO 21 IQ=IP+1,N G=100.D0*ABS(A(IP,IQ)) IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP))) & .AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ)))) THEN A(IP,IQ)=0.D0 ELSEIF(ABS(A(IP,IQ)).GT.TRESH) THEN H=D(IQ)-D(IP) IF(ABS(H)+G.EQ.ABS(H)) THEN T=A(IP,IQ)/H ELSE THETA=0.5D0*H/A(IP,IQ) T=1.D0/(ABS(THETA)+SQRT(1.D0+THETA**2)) IF(THETA.LT.0.D0) T=-T ENDIF C=1.D0/SQRT(1.D0+T**2) S=T*C TAU=S/(1.D0+C) H=T*A(IP,IQ) Z(IP)=Z(IP)-H Z(IQ)=Z(IQ)+H D(IP)=D(IP)-H D(IQ)=D(IQ)+H A(IP,IQ)=0.D0 DO 16 J=1,IP-1 G=A(J,IP) H=A(J,IQ) A(J,IP)=G-S*(H+G*TAU) A(J,IQ)=H+S*(G-H*TAU) 16 CONTINUE DO 17 J=IP+1,IQ-1 G=A(IP,J) H=A(J,IQ) A(IP,J)=G-S*(H+G*TAU) A(J,IQ)=H+S*(G-H*TAU) 17 CONTINUE DO 18 J=IQ+1,N G=A(IP,J) H=A(IQ,J) A(IP,J)=G-S*(H+G*TAU) A(IQ,J)=H+S*(G-H*TAU) 18 CONTINUE DO 19 J=1,N G=V(J,IP) H=V(J,IQ) V(J,IP)=G-S*(H+G*TAU) V(J,IQ)=H+S*(G-H*TAU) 19 CONTINUE NROT=NROT+1 ENDIF 21 CONTINUE 22 CONTINUE DO 23 IP=1,N B(IP)=B(IP)+Z(IP) D(IP)=B(IP) Z(IP)=0.D0 23 CONTINUE 24 CONTINUE PAUSE '50 iterations should never happen' RETURN END C**************************************************************** function cg(j1,m1,j2,m2,j3,m3) C C Clebsh-Gordan Coefficient C**************************************************************** implicit real*8(a-h,o-z) external fact if(m1+m2.ne.m3) then cg=0.d0 return endif if(j3.lt.abs(j1-j2)) then cg=0.d0 return endif if(j3.gt.j1+j2) then cg=0.d0 return endif ka=j1+j2-j3 kb=j3+j1-j2 kc=j2+j3-j1 kd=j1+j2+j3+1 del=sqrt(fact(ka)*fact(kb)*fact(kc)/fact(kd)) cg=0.d0 do 10 n=0,max(j1+j2-j3,j1-m1,j2+m2) ka1=j1+j2-j3-n if(ka1.lt.0.d0) go to 10 ka2=j3-j2+m1+n if(ka2.lt.0.d0) go to 10 ka3=j3-j1-m2+n if(ka3.lt.0.d0) go to 10 ka4=j1-m1-n if(ka4.lt.0.d0) go to 10 ka5=j2+m2-n if(ka5.lt.0.d0) go to 10 an=n*1.d0 cg=cg+(-1.d0)**n & /(fact(n)*fact(ka1)*fact(ka2)*fact(ka3)*fact(ka4)*fact(ka5)) 10 continue 20 cg=cg*sqrt(fact(j1+m1)*fact(j1-m1)) cg=cg*sqrt(fact(j2+m2)*fact(j2-m2)) cg=cg*sqrt(fact(j3+m3)*fact(j3-m3)) cg=cg*sqrt(2.d0*j3+1.d0)*del return end C*************************************************************** SUBROUTINE MATINV(NMAX,C,D) C C subroutine for calculating the inversion of a matrix C, C whose dimension is NMAX C C*************************************************************** IMPLICIT COMPLEX*16 (A-H,O-Z) PARAMETER (NLEVELMAX=25) DIMENSION U(NLEVELMAX),V(NLEVELMAX), & C(NLEVELMAX,NLEVELMAX),D(NLEVELMAX,NLEVELMAX) DETER=1.0 DO 1 M=1,NMAX DO 1 N=1,NMAX D(N,M)=0.0 IF(N.EQ.M) D(N,M)=1.0 1 CONTINUE DO 10 N=1,NMAX T=C(N,N) IF (ABS(T).LT.1.E-10) GO TO3 GO TO 6 3 J=N 4 J=J+1 IF (J.GT.NMAX) GO TO 11 T=C(N,J) DETER=-DETER IF (ABS(T).LE.1.E-10) GO TO 4 DO 5 K=1,NMAX U(K)=C(N,K) V(K)=D(N,K) C(N,K)=C(J,K) D(N,K)=D(J,K) C(J,K)=U(K) 5 D(J,K)=V(K) 6 DO 7 K=1,NMAX IF (K.EQ.N) GO TO 7 A=C(K,N)/C(N,N) DO 8 L=1,NMAX C(K,L)=C(K,L)-A*C(N,L) D(K,L)=D(K,L)-A*D(N,L) 8 CONTINUE 7 CONTINUE B=C(N,N) DETER=B*DETER DO 10 M=1,NMAX C(N,M)=C(N,M)/B D(N,M)=D(N,M)/B 10 CONTINUE RETURN 11 PRINT 12 12 FORMAT('MATRIX NOT INVERTIBLE') RETURN END function fact(n) implicit real*8(a-h,o-z) if(n.lt.0) then fact=0.d0 return endif if (n.eq.0) then fact=1.d0 go to 10 endif fact=1.d0 do 20 i=1,n fact=fact*i*1.d0 20 continue 10 return end C***************************************************************** C Subroutine for the Coulomb wave function C***************************************************************** SUBROUTINE JFLGAM(XD,YD,PAR,PAI,NBCHIF) DOUBLE PRECISION XD,YD,PAR,PAI,TEST,C,HLO2PI,PI,PIS2,PIS4 DOUBLE PRECISION X,Y,U,V,TRA,TRA1,ALO2PI,RAC2,COSI,SINI DOUBLE PRECISION COS2I,SIN2I,ZMOD,DEPI DOUBLE PRECISION XX DOUBLE PRECISION ALOPI DOUBLE PRECISION SUPINT DIMENSION TEST(7),C(6) DATA TEST/2.9152D+7,2.2958D+3,1.4124D+2,3.9522D+1,19.6611D0, &12.791D0,-10.D0/ DATA C/8.333333333333333D-2,-2.777777777777777D-3, &7.936507936507937D-4,-5.952380952380952D-4,8.417508417508418D-4, &-1.917526917526918D-3/ DATA HLO2PI/0.918938533204672/,PI/3.141592653589793/ DATA PIS2/1.570796326794897/,PIS4/0.785398163397448/ DATA ALO2PI/1.837877066409345/,RAC2/0.3465735902799726/ DATA DEPI/6.283185307179586/,ALOPI/1.1447298858494002/ DATA SUPINT/2147483647.D0/ NBCHIF=15 X=DABS(XD) XX=X IF(YD)1,2,1 1 Y=DABS(YD) KR=1 I=DMOD(10.99D0-X,SUPINT) C TRANSLATION IF(I)3,3,4 4 TRA=I X=X+TRA C LOGARITHME(X+IY) (X,Y POSITIFS) 3 IF(X-Y)5,6,7 5 TRA1=X/Y IF(TRA1)8,8,9 8 U=DLOG(Y) V=PIS2 GO TO 10 6 U=RAC2+DLOG(X) V=PIS4 GO TO 10 9 TRA=Y*DSQRT(1.D0+TRA1*TRA1) TRA1=Y/X 11 U=DLOG(TRA) V=DATAN(TRA1) 10 GO TO (12,19,23),KR 7 TRA1=Y/X TRA=X*DSQRT(1.D0+TRA1*TRA1) GO TO 11 12 KR=2 C DEVELOPPEMENT ASYMPTOTIQUE ( X SUPERIEUR A 10 ) TRA=X-0.5D0 PAI=V*TRA+Y*(U-1.D0) PAR=-X+HLO2PI+U*TRA-V*Y ZMOD=X+Y IF(ZMOD-TEST(1))13,13,14 13 TRA=X*X+Y*Y COSI=X/TRA SINI=Y/TRA SIN2I=(SINI*COSI)+(SINI*COSI) COS2I=(COSI+SINI)*(COSI-SINI) K=1 GO TO 15 16 TRA=COSI*COS2I-SINI*SIN2I SINI=SINI*COS2I+COSI*SIN2I COSI=TRA 15 PAR=PAR+C(K)*COSI PAI=PAI-C(K)*SINI K=K+1 IF(ZMOD-TEST(K))16,16,14 C TRANSLATION INVERSE 17 I=I-1 X=I X=XX+X GO TO 3 19 PAR=PAR-U PAI=PAI-V 14 IF(I-1)18,60,17 60 IF(XD)17,61,17 C CONTROLE DU QUADRANT 18 IF(XD)20,61,21 61 TRA=PI*Y IF(TRA-1.D-2)300,300,301 300 PAR= TRA*(2.D0+TRA*(-2.D0+TRA*(1.333333333333333D0+TRA*( &-0.6666666666666666D0+TRA*(0.2666666666666666D0+TRA*( &-0.08888888888888888D0+TRA*0.02539682539682540D0)))))) TRA1=-DLOG(Y)-DLOG(PAR) GO TO 302 301 PAR=1.D0-DEXP(-TRA-TRA) TRA1=-DLOG(Y*PAR) 302 PAR=0.5D0*(ALO2PI-TRA+TRA1) PAI=PAI-PIS2 21 IF(YD)28,100,100 C X+IY CHANGE EN -X-IY 20 TRA=PI*Y PAR=ALO2PI-U-PAR-TRA PAI=PI-V-PAI TRA=DEXP(-TRA-TRA) X=PI*DMOD(X,2.D0) SINI=(1.D0-TRA)*DCOS(X) COSI=(1.D0+TRA)*DSIN(X) KR=3 X=DABS(COSI) Y=DABS(SINI) GO TO 3 23 IF(COSI)24,25,25 24 V=PI-DSIGN(V,SINI) GO TO 26 25 IF(SINI)27,26,26 27 V=-V 26 PAR=PAR-U PAI=PAI-V IF(YD)100,100,28 28 PAI=-PAI C ARGUMENT DANS -PI,PI 100 TRA=DABS(PAI/DEPI) IF(TRA-1.D+15)203,204,204 204 NBCHIF=0 PAI=0.D0 GO TO 29 203 IF(TRA-1.D0)205,205,206 206 NBCHIF=DLOG10(TRA) NBCHIF=14-NBCHIF TRA=DMOD(TRA,SUPINT) PAI=DMOD(TRA,1.D0)*DSIGN(DEPI,PAI) 205 IF(DABS(PAI)-PI)29,29,207 207 PAI=PAI-DSIGN(DEPI,PAI) GO TO 29 C JFLGAM REEL 2 PAI=0.D0 IF(XD)31,32,33 C CONDITIONS D EXISTENCE 32 WRITE (6,1000) 1000 FORMAT (21H JFLGAM(0) EST INFINI) GO TO 50 31 IF(X-4503599627370496.D0)34,35,35 35 WRITE (6,1001) 1001 FORMAT (30H ARGUMENT DE JFLGAM TROP GRAND) GO TO 50 34 Y=DMOD(X,SUPINT) IF(Y)400,36,400 400 IF(Y-0.99D0)33,33,405 405 TRA=IDINT(Y+0.1D0) IF(DABS(Y-TRA)-5.D-15)36,36,33 36 WRITE (6,1002) 1002 FORMAT (28H JFLGAM (-ENTIER) EST INFINI) 50 PAR=1.D+74 NBCHIF=0 GO TO 29 C TRANSLATION 33 I=DMOD(10.99D0-X,SUPINT) IF(I)37,37,38 38 TRA=I X=X+TRA C DEVELOPPEMENT ASYMPTOTIQUE 37 Y=DLOG(X) PAR=-X+HLO2PI+Y*(X-0.5D0) IF(X-TEST(1))39,39,43 39 COSI=1.D0/X COS2I=COSI*COSI K=1 GO TO 41 42 COSI=COSI*COS2I 41 PAR=PAR+C(K)*COSI K=K+1 IF(X-TEST(K))42,42,40 C TRANSLATION INVERSE 44 X=X-1.D0 48 Y=DLOG(X) PAR=PAR-Y I=I-1 40 IF(I-1)43,49,44 49 X=DABS(XD) GO TO 48 C X NEGATIF 43 IF(XD)45,29,29 45 PAR=ALOPI-PAR-Y Y=PI*DMOD(X,2.D0) Y=-DSIN(Y) IF(Y)46,36,47 46 Y=-Y PAI=PI 47 PAR=PAR-DLOG(Y) ENTRY JFLGV1 29 RETURN END SUBROUTINE YFCLEN(ETA,RO,U,UP,V,VP,SIGMA0,IDIV,NN) IMPLICIT COMPLEX*16(A-D,F-H),REAL*8(E,P-Z) C IF(NN.EQ.1)GO TO 20 C ETA2=ETA*ETA FA=DCMPLX(1.D0,ETA) M=0.25*ETA+4.D1 C C POLYNOMES DE TCHEBICHEV JUSQU'AU RANG M C K=M+2 X=(ETA+ETA)/RO XX=X+X-1.D0 T0=1.D0 T1=XX XX=XX+XX DO 6 J=2,K TJ=XX*T1-T0 T0=T1 T1=TJ 6 CONTINUE TM=T1 TL=T0 C C INITIALISATION C AM=(0.D0,0.D0) AL=(1.D-40,1.D-40) BN=(0.D0,0.D0) BM=(1.D-40,1.D-40) BL=(0.D0,0.D0) BK=DCMPLX(4.D0*DFLOAT(M+1),0.D0)*AL+BM F=(0.D0,0.D0) G =(0.D0,0.D0) GP=(0.D0,0.D0) C C RECURRENCE DESCENDANTE C K=M 10 R=K TK=XX*TL-TM TM=TL TL=TK HK=DCMPLX(TK,0.D0) C1=DCMPLX(R*(R+1.D0)-ETA2,ETA*(R+R+1.D0)) C2=(4.D0,0.D0)*DCMPLX(R+1.D0,0.D0) C2=C2*DCMPLX(-R-1.D0,ETA*3.D0) C3=FA*DCMPLX(-R-R-4.D0,ETA) C4=DCMPLX((7.D0*R+5.D0)/4.D0,0.D0) C5=DCMPLX(R+R+2.D0,0.D0) C6=DCMPLX((R+3.D0)/4.D0,0.D0) AK=(C2*AL+C3*AM-C4*BL-C5*BM-C6*BN)/C1 J=K/2 J=K-J-J IF(J)1,2,1 1 F=F-AK GO TO 3 2 F=F+AK 3 Z=ABS(AK) G=G+HK*AK GP=GP+HK*BK C C F=A0/2-A1+A2-A3+A4-A5+... C C CONGRUENCE MODULO 10**60 C IF(Z-1.D60)4,5,5 5 D=(1.D60,0.D0) AK=AK/D AL=AL/D AM=AM/D BK=BK/D BL=BL/D BM=BM/D BN=BN/D F=F/D G=G/D GP=GP/D 4 IF(K)8,8,9 9 D=DCMPLX(4.D0*R,0.D0) BJ=D*AK+BL AM=AL AL=AK BN=BM BM=BL BL=BK BK=BJ K=K-1 GO TO 10 C C NORMALISATION ET CALCUL DE Z(RO) C 8 D=(0.5D0,0.D0)*AK F=F-D G=G-D GP=GP-(0.5D0,0.D0)*BK D=DCMPLX(0.D0,-ETA*DLOG(2.D0)+SIGMA0) AXPO=EXP(D) F=F/AXPO G=G/F GP=GP/F C C CALCUL DE F ET G C D=DCMPLX(0.D0,RO-ETA*DLOG(RO)) AXPO=EXP(D) D=G*AXPO GP=AXPO*(DCMPLX(0.D0,1.D0-ETA/RO)*G-DCMPLX(X/RO,0.D0)*GP) V=D D=(0.D0,-1.D0)*D U=D VP=GP GP=(0.D0,-1.D0)*GP UP=GP IDIV=0 RETURN C C SERIE ORIGINE C 20 PI=3.141592653589793D0 XA=0.577215664901533D0 RO2=RO*RO ETAP=ETA+ETA PIETA=PI*ETA Z=138.15510557964276D0 IDIV=0 IF(DABS(PIETA)-Z)21,21,22 22 INDG=IDINT(PIETA/Z) IDIV=60*INDG IF(ETA.LT.0) IDIV=0 PIETA=PIETA-Z*DFLOAT(INDG) 21 PIETA2=0.5D0*PIETA P=DEXP(PIETA2)*DSQRT(DSINH(PIETA)/PIETA) CALL JFDELG(1.D0,ETA,PAR,PAI,NB) Z1=ETAP*(XA+XA+DLOG(2.D0)-1.D0+PAR) U0=0.D0 U1=RO V0=1.D0 V1=Z1*RO U=U0+U1 V=V0+V1 UP=1.D0 VP=Z1 XN=2.D0 DO 104N=2,10000 XN1=XN*(XN-1.D0) U2=(ETAP*RO*U1-RO2*U0)/XN1 U=U+U2 V2=(ETAP*RO*V1-RO2*V0-ETAP*(XN+XN-1.D0)*U2)/XN1 V=V+V2 UP=UP+XN*U2/RO VP=VP+XN*V2/RO IF(DABS(U2/U).GT.1.D-14)GOTO18 IF(DABS(V2/V).LE.1.D-14)GOTO19 18 U0=U1 U1=U2 V0=V1 V1=V2 XN=XN+1.D0 104 CONTINUE 19 PP=V+ETAP*U*DLOG(RO) W=U/P WP=UP/P V=P*PP VP=P*(VP+ETAP*(UP*DLOG(RO)+U/RO)) U=W UP=WP RETURN END SUBROUTINE YFASYM(ETA,RAU,FO,FPO,GO,GPO,SIGO,IEXP) IMPLICIT REAL*8 (A-H,O-Z) IEXP=0 TRB=0.D0 RAU2=RAU+RAU ETAC=ETA*ETA CALL JFLGAM(1.D0,ETA,TRA,SIGO,NTRUC) 40 N=0 PS=1.D0 GS=0.D0 PT=0.D0 GT=1.D0-ETA/RAU SF=PS SG=GS SPF=PT SPG=GT 45 DENOM= DFLOAT(N+1)*RAU2 AN= DFLOAT(N+N+1)*ETA/DENOM BN= (ETAC-DFLOAT(N*(N+1)))/DENOM PS1=AN*PS-BN*PT GS1=AN*GS-BN*GT-PS1/RAU PT1=AN*PT+BN*PS GT1=AN*GT+BN*GS-PT1/RAU 42 SF=SF+PS1 SG=SG+GS1 SPF=SPF+PT1 SPG=SPG+GT1 N=N+1 IF(N-17)46,48,44 48 TRA=PS*PS+PT*PT 44 TRB=PS1*PS1+PT1*PT1 TEST=TRA-TRB IF(TEST)47,47,46 46 PS=PS1 GS=GS1 PT=PT1 GT=GT1 TRA=TRB GOTO 45 47 TETAO= RAU-ETA*DLOG (RAU2)+SIGO TRA= DSIN(TETAO) TRB=DCOS(TETAO) GO=SF*TRB-SPF*TRA GPO=SG*TRB-SPG*TRA FO=SPF*TRB+SF*TRA FPO=SPG*TRB+SG*TRA RETURN END SUBROUTINE DFCOUL(ETA,RO,F,FP,G,GP,SIGMA,L,IEXP) IMPLICIT REAL*8(A-H,O-Z) DIMENSION F(*),FP(*),G(*),GP(*),IEXP(*),SIGMA(*) DATA DEPI/6.283185307179586D0/ ETAC=ETA*ETA L1=L+1 CALL DFCZ0(ETA,RO,F0,FP0,G0,GP0,S,I) F(1)=F0 FP(1)=FP0 G(1)=G0 GP(1)=GP0 IEXP(1)=I SIGMA(1)=S IF(L)1,1,2 1 RETURN 2 LINF=0 IND=0 IF((ETA.GT.0).AND.(RO.LT.(ETA+ETA)))GO TO 21 Z=ETA+DSQRT(ETAC+DFLOAT(L*(L+1))) IF(RO.GE.Z)GO TO 20 7 ROINF=ETA+DSQRT(ETAC+DFLOAT(LINF*(LINF+1))) IF(RO-ROINF)3,4,4 4 IF(LINF-L)5,6,6 5 LINF=LINF+1 GO TO 7 3 IND=1 6 LIN=LINF+1 20 XM=1.D0 IF(IND.EQ.0)LIN=L1 DO 8 J=2,LIN ZIG=(DSQRT(ETAC+XM*XM))/XM ZAG=ETA/XM+XM/RO F(J)=(ZAG*F(J-1)-FP(J-1))/ZIG FP(J)=ZIG*F(J-1)-ZAG*F(J) G(J)=(ZAG*G(J-1)-GP(J-1))/ZIG GP(J)=ZIG*G(J-1)-ZAG*G(J) IEXP(J)=I SIG=SIGMA(J-1)+DATAN(ETA/(J-1)) IPI=SIG/DEPI SIG=SIG-IPI*DEPI IF(SIG)60,50,70 60 IF(SIG.LT.(-DEPI/2.D0))SIG=SIG+DEPI GO TO 50 70 IF(SIG.GT.(DEPI/2.D0))SIG=SIG-DEPI 50 SIGMA(J)=SIG 8 XM=XM+1.D0 IF(IND.EQ.0)RETURN GO TO 22 21 LIN=1 22 FTEST=F(LIN) FPTEST=FP(LIN) LMAX=LINF+25+IDINT(5.D0*DABS(ETA)) IF(LMAX-L)9,10,10 9 LMAX=L 10 FI=1.D0 FPI=1.D0 13 XM=LMAX+1 ZIG=(DSQRT(ETAC+XM*XM))/XM ZAG=ETA/XM+XM/RO FL=(ZAG*FI+FPI)/ZIG FPL=ZAG*FL-ZIG*FI IF(DABS(FL)-1.D15)26,27,27 26 IF(DABS(FPL)-1.D15)28,27,27 27 FL=FL*1.D-15 FPL=FPL*1.D-15 28 FI=FL FPI=FPL IF(LMAX-L)11,11,12 12 LMAX=LMAX-1 GO TO 13 11 F(LMAX+1)=FL FP(LMAX+1)=FPL IF(LMAX-LINF)15,15,14 14 GO TO 12 15 FACT=FTEST/F(LIN) FACTP=FPTEST/FP(LIN) INDICE=I/60 XM=LINF DO 16 J=LIN,L1 F(J)=F(J)*FACT FP(J)=FP(J)*FACTP 25 IF(J.EQ.1)GO TO 16 ZIG=(DSQRT(ETAC+XM*XM))/XM ZAG=ETA/XM+XM/RO G(J)=(ZAG*G(J-1)-GP(J-1))/ZIG GP(J)=ZIG*G(J-1)-ZAG*G(J) IF(DABS(G(J))-1.D60)17,18,18 17 IF(DABS(GP(J))-1.D60)19,18,18 18 G(J)=G(J)/1.D60 GP(J)=GP(J)/1.D60 INDICE=INDICE+1 19 IEXP(J)=INDICE*60 A=FP(J)*G(J) B=-F(J)*GP(J) IF(A-1.D0)29,30,30 29 I1=IDINT(DLOG10(A)) I2=IDINT(DLOG10(B)) GO TO 31 30 I1=IDINT(DLOG10(A))+1 I2=IDINT(DLOG10(B))+1 31 F(J)=F(J)*1.D1**(-I2) FP(J)=FP(J)*1.D1**(-I1) SIG=SIGMA(J-1)+DATAN(ETA/(J-1)) IPI=SIG/DEPI SIG=SIG-IPI*DEPI IF(SIG)61,51,71 61 IF(SIG.LT.(-DEPI/2.D0))SIG=SIG+DEPI GO TO 51 71 IF(SIG.GT.(DEPI/2.D0))SIG=SIG-DEPI 51 SIGMA(J)=SIG 16 XM=XM+1.D0 RETURN END SUBROUTINE YFIREG(ETA,RO,G0,GP0) IMPLICIT REAL*8(A-H,O-Z) IF(ETA.LE.0.D0)GOTO250 IF(ETA.LE.3.D0)GOTO251 IF(ETA.LE.1.D1)GOTO252 IF(ETA.LE.18.D0)GOTO253 IF(ETA.LE.22.D0)GOTO254 IF(RO.LE.0.3D0+(3.D1-ETA)/8.D1)GOTO200 C SERIE DE TAYLOR DEPART RAU0 300 CONTINUE RAU0=1.666666666666667D0*DABS(ETA)+7.5D0 CALL YFASYM(ETA,RAU0,F0,FP0,G0,GP0,SIGMA0,IEXP) X=RAU0-RO X2=X*X X3=X*X2 UNR=1.D0/RAU0 ETR0=1.D0-2.D0*ETA*UNR U0=G0 U1=-X*GP0 U2=-0.5D0*ETR0*X2*U0 S=U0+U1+U2 V1=U1/X V2=2.D0*U2/X T=V1+V2 XN=3.D0 DO10N=3,10000 C N=N XN1=XN-1.D0 XN1=XN*XN1 U3=X*U2*UNR*(1.D0-2.D0/XN)-ETR0*U1*X2/XN1+X3*U0*UNR/XN1 S=S+U3 V3=XN*U3/X T=T+V3 16 IF(DABS(U3/S).GT.1.D-11)GO TO 11 IF(DABS(V3/T).LE.1.D-11)GO TO 20 11 U0=U1 U1=U2 U2=U3 XN=XN+1.D0 10 CONTINUE 20 G0=S GP0=-T RETURN C SERIE ORIGINE 200 CONTINUE PI=3.141592653589793D0 GA=0.577215664901533D0 ETA2=ETA*ETA RO2=RO*RO ETAP=ETA+ETA PIETA=PI*ETA PIETA2=0.5D0*PIETA B=DEXP(PIETA2)*DSQRT(DSINH(PIETA)/PIETA) CALL JFDELG(1.D0,ETA,PAR,PAI,NB) C1=ETAP*(GA+GA+DLOG(2.D0)-1.D0+PAR) U0=0.D0 U1=RO V0=1.D0 V1=C1*RO U=U0+U1 V=V0+V1 UP=1.D0 VP=C1 XN=2.D0 DO 104N=2,10000 XN1=XN*(XN-1.D0) U2=(ETAP*RO*U1-RO2*U0)/XN1 U=U+U2 V2=(ETAP*RO*V1-RO2*V0-ETAP*(XN+XN-1.D0)*U2)/XN1 V=V+V2 UP=UP+XN*U2/RO VP=VP+XN*V2/RO 17 IF(DABS(U2/U).GT.1.D-14)GOTO18 IF(DABS(V2/V).LE.1.D-14)GOTO19 18 U0=U1 U1=U2 V0=V1 V1=V2 XN=XN+1.D0 104 CONTINUE 19 GP=V+ETAP*U*DLOG(RO) G0=B*GP GP0=B*(VP+ETAP*(UP*DLOG(RO)+U/RO)) RETURN 250 IF(RO.LE.0.5D0*ETA+9.D0)GOTO200 GOTO 300 251 IF(RO.LE.2.25D0+7.35D0*(3.D0-ETA))GOTO200 GOTO 300 252 IF(RO.LE.1.2D0+1.5D-1*(1.D1-ETA))GOTO200 GOTO 300 253 IF(RO.LE.0.6D0+0.75D-1*(18.D0-ETA))GOTO200 GOTO 300 254 IF(RO.LE.0.4D0+0.5D-1*(22.D0-ETA))GOTO200 GOTO 300 END SUBROUTINE YFRICA(ETA,RAU,FO,FPO,GO,GPO,SIGMA0,IDIV) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(5),QP(5) C C COEFFICIENTS RICCATI C DATA G61,G62,G63,G64,G65,G66,G67,G68,G69,G610, &G611/ 0.1159057617187498D-01,0.3863525390624998D-01, &0.4660034179687498D-01,0.4858398437499998D-01, &0.1156514485677080D 01,0.5687475585937496D 01, &0.1323888288225445D 02,0.1713083224826384D 02, &0.1269003295898436D 02,0.5055236816406248D 01, &0.8425394694010415D 00/ DATA G81,G82,G83,G84,G85,G86,G87,G88,G89,G810,G811,G812,G813,G814, &G815/ 0.1851092066083633D-01,0.8638429641723630D-01, &0.1564757823944092D 00,0.1430139541625977D 00, &0.1924622058868408D 00,0.8500803152720129D 01, &0.7265429720878595D 02,0.3057942376817972D 03, &0.7699689544836672D 03,0.1254157054424285D 04, &0.1361719536066055D 04,0.9831831171035763D 03, &0.4547869927883148D 03,0.1222640538215636D 03, &0.1455524450256709D 02/ DATA GP61,GP62,GP63,GP64,GP65,GP66/ 0.2897644042968748D-01, &0.2318115234375000D 00,0.8056640625000000D 00, &0.1601562499999998D 01,0.3046875000000000D 00, &0.5624999999999998D 01/ DATA GP81,GP82,GP83,GP84,GP85,GP86,GP87, &GP88/ 0.6478822231292720D-01,0.6910743713378906D 00, &0.3322952270507811D 01,0.9483032226562498D 01, &0.1769653320312499D 02,0.3478710937499998D 02, &0.5020312499999999D 02,0.7874999999999999D 02/ DATA Q /0.4959570165D-1,0.8888888889D-2,0.2455199181D-2, &0.9108958061D-3,0.2534684115D-3/ DATA QP /0.1728260369D0,0.3174603174D-3,0.3581214850D-2, &0.3117824680D-3,0.9073966427D-3/ CALL JFLGAM(1.D0,ETA,TRA,SIGMA0,IND) RAU2=RAU+RAU RAUC=RAU*RAU ETAC=ETA*ETA ETA2=ETA+ETA ETARO=ETA*RAU ETARO2=ETARO+ETARO PIETA=3.141592653589793*ETA IND=0 JND=0 IG=0 IF(ETA)20,20,21 20 IF(-ETARO-14.0625D0)32,22,22 22 INDICE=1 C RICCATI 3 IDIV=0 GO TO 2 21 IF(DABS(RAU-ETA2).LE.1.D-9)GO TO 18 IF(RAU-ETA2)30,18,31 31 IF(RAU-ETA2-2.D1*(ETA**0.25D0))34,33,33 33 INDICE=0 C RICCATI 2 IDIV=0 GO TO 2 32 NN=1 GO TO 35 34 NN=0 35 CALL YFCLEN(ETA,RAU,FO,FPO,GO,GPO,SIGMA0,IDIV,NN) RETURN 30 IF(ETARO-12.D0)32,32,23 23 TRA=ETA2-6.75D0*(ETA**0.4D0) IF(RAU-TRA)6,6,24 24 IND=1 JND=1 RO=RAU RAU=TRA RAU0=TRA C RICCATI 1 6 X=RAU/ETA2 U= (1.D0-X)/X X2=X*X RU= DSQRT(U) RX= DSQRT(X) TRE= 1.D0/(U*RU*ETA2) TRB=TRE*TRE FI= (DSQRT((1.D0-X)*X)+DASIN(RX)-1.570796326794897D0)*ETA2 TR1= -0.25D0*DLOG(U) 602 TR2= -((9.D0*U+6.D0)*U+5.D0)/48.D0 603 TR3= ((((-3.D0*U-4.D0)*U+6.D0)*U+12.D0)*U+5.D0)/64.D0 604 TR4=- ((((((U+2.D0)*945.D0*U+1395.D0)*U+12300.D0)*U+25191.D0) &*U+19890.D0)*U+5525.D0)/46080.D0 605 TR5= ((((((((-27.D0*U-72.D0)*U-68.D0)*U+360.D0)*U+2190.D0) &*U+4808.D0)*U+5148.D0)*U+2712.D0)*U+565.D0)/2048.D0 606 TR6=- (((((((((G61*U+G62)*U+G63)*U+G64)*U+G65)*U+G66)*U+G67) &*U+G68)*U+G69)*U+G610)*U+G611 607 TR7= ((((((((((((-81.*U-324.)*U-486.)*U-404.)*U+4509.)*U+52344.) &*U+233436.)*U+567864.)*U+838521.)*U+775884.)*U+441450.) &*U+141660.)*U+19675.) /6144. 608 TR8= (((((((((((((G81*U+G82)*U+G83)*U+G84)*U+G85)*U+G86)*U+G87) &*U+G88)*U+G89)*U+G810)*U+G811)*U+G812)*U+G813)*U+G814)*U+G815 PSIP=PSIP+TRA XXX=138.1551055796428D0 FI= FI+TRE*(TR2+TRB*(TR4+TRB*(TR6+TRB*TR8))) PSI=-FI INDG=IDINT(PSI/XXX) IDIV=60*INDG TRA= TR1+TRB*(TR3+TRB*(TR5+TRB*TR7)) FI=FI+TRA PSI=PSI+TRA FIP=RU*ETA2 TRA=1.D0/(X2*U) TR1=0.25D0 TRE= TRE/(X2*X2*U) TRB=TRB/(X2*X2) 702 TR2= -(8.D0*X-3.D0)/32.D0 703 TR3= ((24.D0*X-12.D0)*X+3.D0)/64.D0 704 TR4= (((-1536.D0*X+704.D0)*X-336.D0)*X+63.D0)/2048.D0 705 TR5= ((((1920.D0*X-576.D0)*X+504.D0)*X-180.D0)*X+27.D0)/1024.D0 706 TR6 = ((((-GP66*X+GP65)*X-GP64)*X+GP63)*X-GP62)*X+GP61 707 TR7= - ((((((-40320.D0*X-10560.D0)*X-13248.D0)*X+7560.D0) &*X-3132.D0)*X+756.D0)*X-81.D0) / 2048.D0 708 TR8 =- ((((((GP88*X+GP87)*X+GP86)*X-GP85)*X+GP84)*X-GP83)*X+GP82) &*X-GP81 FIP=FIP +TRE*(TR2+TRB*(TR4+TRB*(TR6+TRB*TR8))) TRA= TRA*(TR1+TRB*(TR3+TRB*(TR5+TRB*TR7))) FIP=FIP+TRA PSIP=-FIP IF(INDG.EQ.0)GO TO 8 PSI=PSI-XXX*DFLOAT(INDG) FI =FI +XXX*DFLOAT(INDG) 8 FO=0.5D0*DEXP(FI) GO= DEXP(PSI) FPO= FO*FIP/ETA2 GPO=GO*PSIP/ETA2 IF(JND.EQ.0)RETURN RAU=RO GO=FO GPO=FPO 27 X=RAU0-RO X2=X*X X3=X*X2 UNR=1.D0/RAU0 ETR0=1.D0-2.D0*ETA*UNR U0=GO U1=-X*GPO U2=-0.5D0*ETR0*X2*U0 S=U0+U1+U2 V1=U1/X V2=2.D0*U2/X T=V1+V2 XN=3.D0 DO10N=3,10000 C N=N XN1=XN-1.D0 XN1=XN*XN1 U3=X*U2*UNR*(1.D0-2.D0/XN)-ETR0*U1*X2/XN1+X3*U0*UNR/XN1 S=S+U3 V3=XN*U3/X T=T+V3 16 IF(DABS(U3/S).GT.1.D-10)GO TO 11 IF(DABS(V3/T).LE.1.D-10)GO TO 17 11 U0=U1 U1=U2 U2=U3 XN=XN+1.D0 10 CONTINUE 17 IF(IG)25,26,25 25 GO=S GPO=-T FO=HO FPO=HPO RETURN 26 HO=S HPO=-T 18 ET0=ETA**(0.166666666666667) ETAD=ETAC*ETAC ET=ETA**(0.6666666666666667) ET1=ET*ET ET2=ET1*ET1 ET3=ET2*ET ET4=ETAD*ET ET5=ET4*ET FO=1.D0-Q(1)/ET1-Q(2)/ETAC-Q(3)/ET3-Q(4)/ETAD-Q(5)/ET5 GO=1.D0+Q(1)/ET1-Q(2)/ETAC+Q(3)/ET3-Q(4)/ETAD+Q(5)/ET5 FPO=1.D0+QP(1)/ET+QP(2)/ETAC+QP(3)/ET2+QP(4)/ETAD+QP(5)/ET4 GPO=1.D0-QP(1)/ET+QP(2)/ETAC-QP(3)/ET2+QP(4)/ETAD-QP(5)/ET4 FO=0.7063326373D0*ET0*FO GO=1.223404016D0*ET0*GO FPO=0.4086957323D0*FPO/ET0 GPO=-0.7078817734D0*GPO/ET0 IDIV=0 IF(IND.EQ.0)RETURN IG=1 RAU0=ETA2 GO TO 27 2 X=ETA2/RAU X2=X*X U=1.D0-X RU= DSQRT(U) U3=U*U*U TRD= 1.D0/(U3*ETA2*ETA2) TRC=X2*TRD TRE=1.D0/(U*RU*ETA2) FI= -0.25D0*DLOG(U) TRB=TRD/64.D0 TR3= (((3.D0*U-4.D0)*U-6.D0)*U+12.D0)*U-5.D0 501 TR5= ((((((((-27.D0*U+72.D0)*U-68.D0)*U-360.D0)*U+2190.D0) &*U-4808.D0)*U+5148.D0)*U-2712.D0)*U+565.D0)/32.D0 502 TR7= ((((((((((((81.D0*U-324.D0)*U+486.D0)*U-404.D0)*U-4509.D0) &*U+52344.D0)*U-233436.D0)*U+567864.D0)*U-838521.D0)*U+775884.D0) &*U-441450.D0)*U+141660.D0)*U-19675.D0)/96.D0 FI= FI+TRB*(TR3+TRD*(TR5+TRD*TR7)) FIP=0.25D0/U TRB=3.D0*TRC/(64.D0*U) TR3= (X-4.D0)*X+8.D0 511 TR5= ((((9.D0*X-60.D0)*X+168.D0)*X-192.D0)*X+640.D0)/16.D0 512 TR7= ((((((-27.D0*X+252.D0)*X-1044.D0)*X+2520.D0)*X-4416.D0) &*X-3520.D0)*X-13440.D0)/32.D0 FIP =FIP+TRB*(TR3+TRC*(TR5+TRC*TR7)) TRA= DABS((RU-1.D0)/(RU+1.D0)) PSI= (0.5D0*DLOG(TRA)+RU/X)*ETA2+0.785398163397448D0 TR2= -((9.D0*U-6.D0)*U+5.D0)/48.D0 521 TR4= ((((((U-2.D0)*945.D0*U+1395.D0)*U-12300.D0)*U+25191.D0) &*U-19890.D0)*U+5525.D0)/46080.D0 522 TR6 = (((((((((-G61*U+G62)*U-G63)*U+G64)*U-G65)*U+G66)*U-G67) &*U+G68)*U-G69)*U+G610)*U-G611 523 TR8= (((((((((((((G81*U-G82)*U+G83)*U-G84)*U+G85)*U-G86)*U+G87) &*U-G88)*U+G89)*U-G810)*U+G811)*U-G812)*U+G813)*U-G814)*U+G815 PSI= PSI+TRE*(TR2+TRD*(TR4+TRD*(TR6+TRD*TR8))) PSIP = -RU*ETA2/X2 TRB=TRE*X/U TR2= (3.D0*X-8.D0)/32.D0 531 TR4= - (((63.D0*X-336.D0)*X+704.D0)*X-1536.D0)/2048.D0 532 TR6 = ((((GP61*X-GP62)*X+GP63)*X-GP64)*X+GP65)*X-GP66 533 TR8 = ((((((-GP81*X+GP82)*X-GP83)*X+GP84)*X-GP85)*X+GP86)*X+GP87) &*X+GP88 PSIP =PSIP+ TRB*(TR2+TRC*(TR4+TRC*(TR6+TRC*TR8))) TRA= DEXP(FI) FO= TRA*DSIN(PSI) GO= TRA*DCOS(PSI) IF(INDICE)535,536,535 535 TRA=FO FO=GO GO=-TRA 536 TRA=-ETA2/RAUC FPO=(FIP*FO+PSIP*GO)*TRA GPO=(FIP*GO-PSIP*FO)*TRA RETURN END SUBROUTINE DFCZ0(ETA,RO,F0,FP0,G0,GP0,SIGMA0,IEXP) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A1(110),A2(110),B1(110),B2(110) IF(RO)2,2,1 2 WRITE (6,1000) 1000 FORMAT (21H RO NEGATIF OU NUL **) RETURN 1 IF(ETA-30.D0)3,3,4 4 IF(DABS(ETA)-5.D2)28,28,29 28 CALL YFRICA(ETA,RO,F0,FP0,G0,GP0,SIGMA0,IEXP) RETURN 29 WRITE (6,1001) 1001 FORMAT (42H VALEUR ABSOLUE DE ETA SUPE-&EU-E A 500 **) RETURN 3 IF(ETA+8.D0)4,5,5 5 IF(ETA)6,7,6 7 F0=DSIN(RO) G0=DCOS(RO) FP0=G0 GP0=-F0 IEXP=0 SIGMA0=0.D0 RETURN 6 BORNE=1.666666666666667D0*DABS(ETA)+7.5D0 IF(RO-BORNE)8,9,9 9 CALL YFASYM(ETA,RO,F0,FP0,G0,GP0,SIGMA0,IEXP) RETURN 8 IF(ETA-10.D0)10,11,11 10 IF(ETA)12,12,13 13 IF(RO-2.D0)14,12,12 11 IF(ETA-(5.D0*RO+6.D1)/7.D0)12,12,14 12 CALL YFASYM(ETA,BORNE,F0,FP0,G0,GP0,SIGMA0,IEXP) H=BORNE DH=F0/H IF(ETA)20,21,21 20 N=-0.5D0*ETA+5.D0 GO TO 22 21 N=ETA/5.D0+5.D0 22 N=5*(N+1) Z=4.D0/H Y=1.D0-(ETA+ETA)*Z A1(N+2)=1.D-55 A1(N+3)=0.D0 A1(N+4)=1.D-64 B1(N+3)=1.D-50 B1(N+4)=1.D-68 A2(N+2)=0.D0 A2(N+3)=1.D-74 A2(N+4)=1.D-53 B2(N+3)=0.D0 B2(N+4)=1.D-66 M=N+2 DI=N DO 23 II=2,M I=M-II+2 B1(I)=B1(I+2)+Z*(DI+1.D0)*A1(I+1) S=A1(I+2)+Y*(A1(I+1)-A1(I)) Q=(DI+2.D0)*B1(I)+(DI-1.D0)*B1(I+1) A1(I-1)=S-Z*Q B2(I)=B2(I+2)+Z*(DI+1.D0)*A2(I+1) S=A2(I+2)+Y*(A2(I+1)-A2(I)) Q=(DI+2.D0)*B2(I)+(DI-1.D0)*B2(I+1) A2(I-1)=S-Z*Q IF(I.GE.N)GO TO 23 D=-(B2(I+2)+B2(I))/(B1(I+2)+B1(I)) DO 24 J=I,M A2(J)=A2(J)+D*A1(J) B2(J)=B2(J)+D*B1(J) 24 CONTINUE A2(I-1)=A2(I-1)+D*A1(I-1) 23 DI=DI-1.D0 Q=A1(3)-A1(1) C=A2(3)-A2(1) C=Q/C X1=0.5D0*(A1(2)-C*A2(2)) DO 25 I=3,M X1=X1+A1(I)-C*A2(I) 25 CONTINUE X1=DH/X1 X2=-C*X1 DO 26 I=2,M B1(I)=X1*B1(I)+X2*B2(I) A1(I)=X1*A1(I)+X2*A2(I) 26 CONTINUE A1(1)=X1*A1(1)+X2*A2(1) B1(1)=0.D0 X=RO/H Y=2.D0*X-1.D0 T1=1.D0 T2=Y RESULT=0.5D0*A1(2)+Y*A1(3) DERIVE=0.5D0*B1(2)+Y*B1(3) DO 27 I=2,N TI=2.D0*Y*T2-T1 T1=T2 T2=TI RESULT=RESULT+TI*A1(I+2) DERIVE=DERIVE+TI*B1(I+2) 27 CONTINUE F0=RESULT*RO FP0=DERIVE*RO+RESULT GO TO 30 C SERIE ORIGINE REGULIERE 14 PI=3.141592653589793D0 CALL JFLGAM(1.D0,ETA,TRA,SIGMA0,NTRUC) IEXP=0 RO2=RO*RO ETAP=ETA+ETA PIETA=PI*ETA PIETA2=0.5D0*PIETA B=DEXP(PIETA2)*DSQRT(DSINH(PIETA)/PIETA) U0=0.D0 U1=RO U=U0+U1 UP=1.D0 XN=2.D0 DO 15 N=2,10000 XN1=XN*(XN-1.D0) U2=(ETAP*RO*U1-RO2*U0)/XN1 U=U+U2 UP=UP+XN*U2/RO 17 IF(DABS(U2/U).LT.1.D-10)GO TO 19 18 U0=U1 U1=U2 XN=XN+1.D0 15 CONTINUE 19 F0=U/B FP0=UP/B 30 CALL YFIREG(ETA,RO,G0,GP0) RETURN END SUBROUTINE JFDELG (XD,YD,PAR,PAI,NBCHIF) DOUBLE PRECISION XD,YD,PAR,PAI,TEST,C,PI DOUBLE PRECISION X,Y,U,V,TRA,TRA1,COSI,SINI DOUBLE PRECISION COS2I,SIN2I,ZMOD,DEPI DOUBLE PRECISION TRB,XX DOUBLE PRECISION RAC2,PIS4 DOUBLE PRECISION SUPINT DIMENSION TEST(7),C(6) DATA TEST/2.9152D+7,2.2958D+3,1.4124D+2,3.9522D+1,19.6611D0, &12.791D0,-10.D0/ DATA RAC2/0.3465735902799726/,PIS4/0.785398163397448/ DATA C/8.333333333333333D-2,-8.33333333333333D-3, &3.968253968253968D-3,-4.166666666666667D-3,7.575757575757576D-3, &-2.109279609279609D-2/ DATA PI/3.141592653589793/ DATA DEPI/6.283185307179586/ DATA SUPINT/2147483647.D0/ X=DABS(XD) XX=X NBCHIF=15 IF(YD)1,2,1 1 Y=DABS(YD) KR=1 I=DMOD(10.99D0-X,SUPINT) C TRANSLATION IF(I)3,3,4 4 TRA=I X=X+TRA C LOGARITHME(X+IY) (X,Y POSITIFS) 3 IF(X-Y)5,6,7 5 TRA1=X/Y TRB=1.D0+TRA1*TRA1 TRA=Y*DSQRT(TRB) SINI=1./(TRB*Y) COSI=SINI*TRA1 TRA1=Y/X GO TO 11 6 U=RAC2+DLOG(X) V=PIS4 SINI=0.5D0/X COSI=SINI GO TO 10 7 TRA1=Y/X TRB=1.D0+TRA1*TRA1 TRA=X*DSQRT(TRB) COSI=1./(TRB*X) SINI=COSI*TRA1 11 U=DLOG(TRA) V=DATAN(TRA1) C DEVELOPPEMENT ASYMPTOTIQUE ( X SUPERIEUR A 10 ) 10 PAR=U-0.5*COSI PAI=V+0.5*SINI ZMOD=X+Y IF(ZMOD-TEST(1))13,13,14 13 SIN2I=(SINI*COSI)+(SINI*COSI) COS2I=(COSI+SINI)*(COSI-SINI) SINI=SIN2I COSI=COS2I K=1 GO TO 15 16 TRA=COSI*COS2I-SINI*SIN2I SINI=SINI*COS2I+COSI*SIN2I COSI=TRA 15 PAR=PAR-C(K)*COSI PAI=PAI+C(K)*SINI K=K+1 IF(ZMOD-TEST(K))16,16,14 C TRANSLATION INVERSE 17 I=I-1 X=I X=XX+X 56 IF(X-Y)55,55,57 55 TRA1=X/Y TRB=X*TRA1+Y SINI=1.D0/TRB COSI=TRA1/TRB GO TO 19 57 TRA1=Y/X TRB=X+Y*TRA1 COSI=1.D0/TRB SINI=TRA1/TRB 19 PAR=PAR-COSI PAI=PAI+SINI 14 IF(I)18,18,17 C CONTROLE DU QUADRANT 18 IF(XD)20,61,21 61 TRA=PI*Y IF(TRA-1.D-2)300,300,301 300 TRB= TRA*(2.D0+TRA*(-2.D0+TRA*(1.333333333333333D0+TRA*( &-0.6666666666666666D0+TRA*(0.2666666666666666D0+TRA*( &-0.08888888888888888D0+TRA*0.02539682539682540D0)))))) TRB=(2.D0-TRB)/TRB GO TO 302 301 TRB= DEXP(-TRA-TRA) TRB=(1.D0+TRB)/(1.D0-TRB) 302 PAI=0.5D0*(1.D0/Y+PI*TRB) 21 IF(YD)28,100,100 C X+IY CHANGE EN -X-IY 20 TRA=DEXP(-DEPI*Y) TRB=TRA*TRA COS2I=DEPI*DMOD(X,1.D0) SIN2I=-2.D0*TRA*DCOS(COS2I)+1.D0+TRB PAR=PAR+COSI+DEPI*TRA*DSIN(COS2I)/SIN2I PAI=PAI-SINI+PI*(TRB-1.D0)/SIN2I IF(YD)100,100,28 28 PAI=-PAI 35 WRITE (6,1001) 1001 FORMAT (30H ARGUMENT DE JFDELG TROP GRAND) GO TO 50 34 Y=DMOD(X,SUPINT) IF(Y) 400,36,400 400 IF(Y-0.99D0) 33,33,405 405 TRA= IDINT(Y+0.1D0) IF(DABS(Y-TRA)-5.D-15)36,36,33 31 IF(X-4503599627370496.D0)34,35,35 C ARGUMENT DANS -PI,PI 100 TRA=DABS(PAI/DEPI) IF(TRA-1.D+15)203,204,204 204 NBCHIF=0 PAI=0.D0 GO TO 29 203 IF(TRA-1.D0)205,205,206 206 NBCHIF=DLOG10(TRA) NBCHIF=14-NBCHIF TRA=DMOD(TRA,SUPINT) PAI=DMOD(TRA,1.D0)*DSIGN(DEPI,PAI) 205 IF(DABS(PAI)-PI)29,29,207 207 PAI=PAI-DSIGN(DEPI,PAI) GO TO 29 C DELGAMMA REEL 2 PAI=0.D0 IF(XD)31,32,33 C CONDITIONS D EXISTENCE 32 WRITE (6,1000) 1000 FORMAT (21H JFDELG(0) EST INFINI) GO TO 50 36 WRITE (6,1002) 1002 FORMAT (28H JFDELG (-ENTIER) EST INFINI) 50 PAR=1.D+74 NBCHIF=0 GO TO 29 C TRANSLATION 33 I=DMOD(10.99D0-X,SUPINT) IF(I)37,37,38 38 TRA=I X=X+TRA C DEVELOPPEMENT ASYMPTOTIQUE 37 Y=DLOG(X) PAR=Y-0.5D0/X IF(X-TEST(1))39,39,43 39 COS2I=1.D0/(X*X) COSI=COS2I K=1 GO TO 41 42 COSI=COSI*COS2I 41 PAR=PAR-C(K)*COSI K=K+1 IF(X-TEST(K))42,42,40 C TRANSLATION INVERSE 44 I=I-1 X=I X=XX+X PAR=PAR-1.D0/X 40 IF(I)43,43,44 C X NEGATIF 43 IF(XD)45,29,29 45 PAR=PAR+1.D0/X Y=PI*DMOD(X,2.D0) PAR=PAR+PI*DCOS(Y)/DSIN(Y) ENTRY JFDEV1 29 RETURN END