C C PROGRAM COULEIK C------------------------------------------------------------------------- C C Coulomb fragmentation cross sections. Use as input photonuclear cross C sections parameterized as Lorentzians. Sequential excitation of C one- and two-phonons are considered. C Valid for energies greater than 30 MeV/nucleon. C Based on C Eqs. (13)-(15) of C.A. Bertulani and A.M. Nathan, Nucl. Phys. A554 C (1993) 158. C The excitation probabilities are calculated via the harmonic C oscillator model. C C------------------------------------------------------------------------- C written by C.A. Bertulani C GSI, Planckstr. 1, D-64291 Darmstadt, Germany C Last modified version Jan/15/95 C------------------------------------------------------------------------- C Input should be placed in file couleik.in. An explanation of how C to make an input file is given at the end of the program C C Inportant variables: C NRMAX : no. of grid points in radial coord. of optical pot. C also used as no. of grid points for Z-integration and C b-integration in the eikonal approximation. C An optimal value is 200. C NEX : no. of grid points for the excitation energy. C RMAX : maximum value of radial coordinate used in the integrals C An optimal value is 25. C C See the end of this program for more details C C Units are in MeV and fm. Output cross sections are in barns/sr and C barns. C---------------------------------------------------------------------- PROGRAM COULEIK CALL EWA END C SUBROUTINE EWA IMPLICIT NONE INTEGER*4 NRMAX,NEX,IW,I,K,L,ISO,MU,IFOLD,N,J,NINT PARAMETER(NRMAX=201,NEX=401,NINT=1001) REAL*8 PI,ALPHA,HBC,AMU,APMASS,ATMASS,PTMASS,PLAB,GAMMA,G,WID, & BETA,ETA,EULER,DELR,DELEX,R(NRMAX),B(NRMAX),RP,EMAX,STR,E0, & Z(NRMAX),AP,ZP,AT,ZT,ELAB,FERMI,EX(NEX),FRAC,GM,SEL, & WN(NRMAX),PCM,ARED,RT,RTOT,DELDUM,WG0,ECA,SUM1,SUM2, & H2M,FAK(0:130),FACT,RRHO,RHO0,ARHO,BKNP,BKNM,COMFAC, & AUX,BESSK,BKN,ESQ,RHO01,PHOTO(NEX),WGN(NRMAX),EMIN, & RHO02,RRHO1,RRHO2,RHOG,RGV1,RGV2,RHOG1,RHOG2,SIGNN,SUM,TEMP, & ARHO1,ARHO2,RMAX,CRHO,CRHO1,CRHO2,ACLOSE,AN,EE,PB(NRMAX), & DENS1(NRMAX),DENS2(NRMAX),DUM,FAD(0:130),DELBM, & VIRT(NRMAX,NEX),PROB(NRMAX,NEX),PROBM(NRMAX,NEX),SEC(NEX), & PROB2M(NRMAX,NEX),BM(NRMAX),SEC1(NEX),SEC2(NEX),PEE,PEEM, & VIRTM(NRMAX,NEX),SEC2M(NEX),PROB2(NRMAX,NEX),PBM(NRMAX), & VIRTSC(NRMAX,NEX),BMIN,PROBSC(NRMAX,NEX),PBSC(NRMAX), & SECSC1(NEX),SECSC2(NEX),PROSC2(NRMAX,NEX), & PEESC,SECSC(NEX),SUMSC1,SUMSC2,SUMSC,INTR(NRMAX),INTE(NEX), & INTT(NINT),RHONP,RHOPP COMPLEX*16 Y,PHASEN(NRMAX),GLM COMMON/A1/HBC,ALPHA,BETA,DELR COMMON/A2/EULER COMMON /GFVFAK / FAK COMMON /GFVFAD / FAD COMMON /INT/ INTR,INTE,INTT C OPEN(10, FILE='COULEIK.IN', STATUS='OLD') OPEN(11, FILE='COULEIK.OUT', STATUS='UNKNOWN') OPEN(12, FILE='GARB1.OUT', STATUS='UNKNOWN') OPEN(13, FILE='GARB2.OUT', STATUS='UNKNOWN') OPEN(14, FILE='GARB4.OUT', STATUS='UNKNOWN') C CALL GFV(10) C ***** C Inputs C ***** READ(10,*)AP,ZP,AT,ZT,ECA,RMAX,IFOLD,BMIN READ(10,*)RGV1,RGV2,RRHO,ARHO,CRHO READ(10,*)RRHO1,ARHO1,CRHO1,RRHO2,ARHO2,CRHO2 READ(10,*)SIGNN READ(10,*)IW,L,ISO READ(10,*)E0,WID,FRAC,EMIN,EMAX WRITE(11,'('' **** This is program COULEIK ***'')') WRITE(11,1)AP,ZP,AT,ZT,IFOLD,ECA,L,WID,E0,ISO,FRAC 1 FORMAT(1X,/,1X,' Nuclear mass and charge numbers',/,2X, 1 4(F5.1,2X),/,2X, 2 ' Option for the optical potentials',/,2X, 3 I2,/,2X, 4 ' Laboratory energy/per nucleon',/,2X, 5 F10.3,/,2X, 6 ' Multipolarity, width and energy of the state',/,2X, 7 I2,3X,F10.3,3X,F10.3,/,2X, 8 ' Isoscalar (0) or Isovector (1), Fraction of the sum rule' 9 ,/,2X,I2,3X,F10.3) C ***** C Integers used in Simpson's integrations C ***** INTR(1)=1. INTE(1)=1. INTT(1)=1. INTR(NRMAX)=1. INTE(NEX)=1. INTT(NINT)=1. G=4. DO I=2,NRMAX-1 INTR(I)=G G=6.-G ENDDO DO I=2,NEX-1 INTE(I)=G G=6.-G ENDDO DO I=2,NINT-1 INTT(I)=G G=6.-G ENDDO C ***** C Constants, nuclear sizes, energies, etc. C ***** ELAB=AP*ECA Y=DCMPLX(0.,1.) PI=ACOS(-1.) ALPHA=1./137. HBC=197.33 H2M=20.73 AMU=931.5 ESQ=1.44 APMASS=AMU*AP ATMASS=AMU*AT PTMASS=(AP+AT)*AMU PLAB=SQRT(2.*APMASS*ELAB)/HBC GAMMA=1.+ELAB/APMASS BETA=SQRT(GAMMA**2-1.)/GAMMA ETA=ZP*ZT*ALPHA/BETA ARED=APMASS*ATMASS/PTMASS PCM=ARED/APMASS*PLAB ACLOSE=ETA/PCM EULER=0.5772156649 C ***** C Now, index T is reserved for the excited nucleus (depends on option IW) C ***** IF(IW.EQ.1) THEN DUM=AP AUX=ZP AP=AT ZP=ZT AT=DUM ZT=AUX ENDIF RP=1.2*AP**(1./3.) RT=1.2*AT**(1./3.) RTOT=RP+RT C ***** C grid in R,B,Z (0 - rmax) C ***** DELR=RMAX/NRMAX DO I=1,NRMAX R(I)=I*DELR Z(I)=I*DELR B(I)=I*DELR ENDDO C ***** C grid for b > rmax (up to 300 fm) C ***** DELBM=300./NRMAX DO I=1,NRMAX BM(I)=RMAX+I*DELBM ENDDO C ***** C grid in excitation energy C ***** DELEX=EMAX/(NEX-1) DO I=1,NEX EX(I)=I*DELEX ENDDO C C------------------------------------------------------------------------ C Gamma-neutron cross-sections normalized to a fraction of the sum rule C------------------------------------------------------------------------ C sum rules for the reduced matrix elements C IF(L.EQ.1) THEN SEL=14.8*(AT-ZT)*ZT/AT*ESQ ELSE SEL=5.*L*(2.*L+1)*ZT**2/AT*RT**(2.*L-2.)*ESQ ENDIF IF(L.GT.1.AND.ISO.EQ.1)SEL=(AT-ZT)/ZT*SEL C ***** C photo nuclear cross sections C ***** AUX=HBC**(2.*L-1)*L*FAD(2*L+1)**2/(2.*PI)**3/(L+1) SUM=0. DELDUM=5.*E0/NINT DO K=1,NINT DUM=K*DELDUM IF(DUM.GE.EMIN) THEN SUM=SUM+INTT(K)/(1.+(DUM**2-E0**2)**2/DUM**2/WID**2) & /DUM**(2.*L-2.) ENDIF ENDDO SUM=AUX*SUM/3.*DELDUM STR=FRAC*SEL/SUM DO I=1,NEX PHOTO(I)=0. IF(EX(I).GE.EMIN) THEN PHOTO(I)=STR/(1.+(EX(I)**2-E0**2)**2/EX(I)**2/WID**2) ENDIF ENDDO C C----------------------------------------------------------------- C Calculate absorptive potential (depends on option IFOLD) C----------------------------------------------------------------- C WRITE(11,'('' Parameters of the optical potential '')') C GO TO (12,13,14,15,16,17,18),IFOLD-1 C 12 WRITE(11,112)SIGNN,RGV1,RRHO,ARHO,CRHO 112 FORMAT(' SIGNN,RGV1,RRHO,ARHO,CRHO',/,1X,5(F8.3,1x)) RHOG=AP/(PI**1.5*RGV1**3) WG0=-1./2.*HBC*BETA*SIGNN*RHOG DO I=1,NRMAX WGN(I)=WG0*EXP(-(R(I)/RGV1)**2) DENS1(I)=FERMI(R(I),RRHO,ARHO)*(1.+CRHO*R(I)**2/RRHO**2) ENDDO CALL NORM(R,DELR,DENS1,RHO0) RHO0=AT*RHO0 CALL TWOFOLD(RHO0,R,RMAX,DELR,WGN,DENS1,WN) GO TO 20 C 13 WRITE(11,113)SIGNN,RGV1,RGV2 113 FORMAT(' SIGNN,RGV1,RGV2',/,1X,3(F8.3,1x)) RHOG1=AP/(PI**1.5*RGV1**3) RHOG2=AT/(PI**1.5*RGV2**3) FACT=0.5*PI**1.5*SIGNN*HBC*BETA*RHOG1*RHOG2* & (RGV1*RGV2)**3/(RGV1**2+RGV2**2)**1.5 DO I=1,NRMAX WN(I)=-FACT*EXP(-R(I)**2/(RGV1**2+RGV2**2)) ENDDO GO TO 20 C 14 WRITE(11,114)SIGNN,RRHO1,ARHO1,CRHO1,RRHO2,ARHO2,CRHO2 114 FORMAT(' SIGNN,RRHO1,ARHO1,CRHO1,RRHO2,ARHO2,CRHO2', 1/,1X,7(F8.3,1x)) DO I=1,NRMAX DENS1(I)=FERMI(R(I),RRHO1,ARHO1)*(1.+CRHO1*R(I)**2 & /RRHO1**2) DENS2(I)=FERMI(R(I),RRHO2,ARHO2)*(1.+CRHO2*R(I)**2 & /RRHO2**2) ENDDO CALL NORM(R,DELR,DENS1,RHO01) RHO01=AP*RHO01 CALL NORM(R,DELR,DENS2,RHO02) RHO02=AT*RHO02 FACT=-SIGNN*HBC*BETA/2.*RHO01*RHO02 CALL TWOFOLD(FACT,R,RMAX,DELR,DENS1,DENS2,WN) GO TO 20 C 15 WRITE(11,115)SIGNN,RGV1,RRHO,ARHO,CRHO 115 FORMAT(' SIGNN,RGV1,RRHO,ARHO,CRHO', 1/,1X,5(F8.3,1x)) RHOG=AP/(PI**1.5*RGV1**3) WG0=-1./2.*HBC*BETA*SIGNN*RHOG DO I=1,NRMAX WGN(I)=WG0*EXP(-(R(I)/RGV1)**2) TEMP=(R(I)**2-RRHO**2)/ARHO**2 IF(TEMP.LT.20.) THEN DENS1(I)=1./(1.+EXP(TEMP))*(1.+CRHO*R(I)**2/RRHO**2) ELSE DENS1(I)=0. ENDIF ENDDO CALL NORM(R,DELR,DENS1,RHO0) RHO0=AT*RHO0 CALL TWOFOLD(RHO0,R,RMAX,DELR,WGN,DENS1,WN) GO TO 20 C 16 WRITE(11,116)SIGNN,RRHO1,ARHO1,CRHO1,RRHO2,ARHO2,CRHO2 116 FORMAT(' SIGNN,RRHO1,ARHO1,CRHO1,RRHO2,ARHO2,CRHO2', 1/,1X,7(F8.3,1x)) DO I=1,NRMAX TEMP=(R(I)**2-RRHO1**2)/ARHO1**2 IF(TEMP.LT.20.) THEN DENS1(I)=1./(1.+EXP(TEMP))*(1.+CRHO1*R(I)**2/RRHO1**2) ELSE DENS1(I)=0. ENDIF DENS2(I)=FERMI(R(I),RRHO2,ARHO2) & *(1.+CRHO2*R(I)**2/RRHO2**2) ENDDO CALL NORM(R,DELR,DENS1,RHO01) RHO01=AP*RHO01 CALL NORM(R,DELR,DENS2,RHO02) RHO02=AT*RHO02 FACT=-SIGNN*HBC*BETA/2.*RHO01*RHO02 CALL TWOFOLD(FACT,R,RMAX,DELR,DENS1,DENS2,WN) GO TO 20 C 17 WRITE(11,117)SIGNN,RRHO1,ARHO1,CRHO1,RRHO2,ARHO2,CRHO2 117 FORMAT(' SIGNN,RRHO1,ARHO1,CRHO1,RRHO2,ARHO2,CRHO2', 1/,1X,7(F8.3,1x)) DO I=1,NRMAX TEMP=(R(I)**2-RRHO1**2)/ARHO1**2 IF(TEMP.LT.20.) THEN DENS1(I)=1./(1.+EXP(TEMP))*(1.+CRHO1*R(I)**2/RRHO1**2) ELSE DENS1(I)=0. ENDIF TEMP=(R(I)**2-RRHO2**2)/ARHO2**2 IF(TEMP.LT.20.) THEN DENS2(I)=1./(1.+EXP(TEMP))*(1.+CRHO2*R(I)**2/RRHO2**2) ELSE DENS2(I)=0. ENDIF ENDDO CALL NORM(R,DELR,DENS1,RHO01) RHO01=AP*RHO01 CALL NORM(R,DELR,DENS2,RHO02) RHO02=AT*RHO02 FACT=-SIGNN*HBC*BETA/2.*RHO01*RHO02 CALL TWOFOLD(FACT,R,RMAX,DELR,DENS1,DENS2,WN) C 18 WRITE(11,118) 118 FORMAT(' ---- IFOLD=8 ---- DROPLET DENSITIES WERE USED') DO I=1,NRMAX DENS1(I)=RHONP(R(I),AP,ZP)+RHOPP(R(I),AP,ZP) DENS2(I)=RHONP(R(I),AT,ZT)+RHOPP(R(I),AT,ZT) ENDDO CALL NORM(R,DELR,DENS1,RHO01) RHO01=AP*RHO01 CALL NORM(R,DELR,DENS2,RHO02) RHO02=AT*RHO02 FACT=-SIGNN*HBC*BETA/2.*RHO01*RHO02 CALL TWOFOLD(FACT,R,RMAX,DELR,DENS1,DENS2,WN) C C ***** C initialize nuclear phase subroutine C ***** 20 CALL PHNUC(Y,WN,B,Z,PHASEN) C C------------------------------------------------------------------------ C Coulomb fragmentation cross sections C------------------------------------------------------------------------ C C Equivalent photon number as a function of the excitation energy (I) C and of impact parameter (K) C COMFAC=ZP**2*ALPHA*L*FAD(2*L+1)**2/(2.*PI)**3/(L+1) DO I=1,NEX DO K=1,NRMAX VIRT(K,I)=0. VIRTM(K,I)=0. VIRTSC(K,I)=0. DO MU=0,L IF(MU.EQ.0) THEN J=1 ELSE J=2 ENDIF DUM=EX(I)*B(K)/GAMMA/HBC/BETA BKN=BESSK(MU,DUM) VIRT(K,I)=VIRT(K,I)+J*COMFAC*ABS(GLM(L,MU,1./BETA))**2* & (DUM/B(K))**2*BKN**2*EXP(-2*DIMAG(PHASEN(K))) C ***** C to compare with sharp cutoff approximation C ***** IF(B(K).GE.BMIN-DELR/2.) THEN VIRTSC(K,I)=VIRTSC(K,I)+J*COMFAC* & ABS(GLM(L,MU,1./BETA))**2*(DUM/B(K))**2*BKN**2 ENDIF C ***** C analytical part for b > RMAX C ***** DUM=EX(I)*BM(K)/GAMMA/HBC/BETA BKN=BESSK(MU,DUM) VIRTM(K,I)=VIRTM(K,I)+J*COMFAC*ABS(GLM(L,MU,1./BETA))**2* & (DUM/BM(K))**2*BKN*BKN ENDDO ENDDO ENDDO C C-------------------------------------------------------------- C cross section (b/MeV) as a function of the excitation energy C-------------------------------------------------------------- C contribution from one-phonon excitation C-------------------------------------------------------------- C C excitation probability (1/MeV) as a function of energy (I) and C impact parameter (K) - first-order perturbation theory C DO K=1,NRMAX DO I=1,NEX PROB(K,I)=VIRT(K,I)/EX(I)*PHOTO(I) PROBM(K,I)=VIRTM(K,I)/EX(I)*PHOTO(I) PROBSC(K,I)=VIRTSC(K,I)/EX(I)*PHOTO(I) ENDDO ENDDO C ***** C integration over excitation energy - One-phonon excitation prob. C as a function of impact parameter - harmonic oscillator probability C ***** DO K=1,NRMAX DO I=1,NEX PB(K)=PB(K)+INTE(I)*PROB(K,I)*DELEX/3. PBM(K)=PBM(K)+INTE(I)*PROBM(K,I)*DELEX/3. PBSC(K)=PBSC(K)+INTE(I)*PROBSC(K,I)*DELEX/3. ENDDO ENDDO DO K=1,NRMAX DO I=1,NEX PROB(K,I)=PROB(K,I)*EXP(-PB(K)) PROBM(K,I)=PROBM(K,I)*EXP(-PBM(K)) PROBSC(K,I)=PROBSC(K,I)*EXP(-PBSC(K)) ENDDO ENDDO DO K=1,NRMAX WRITE(12,*)B(K),PB(K)*EXP(-PB(K)) ENDDO DO K=1,NRMAX WRITE(12,*)BM(K),PBM(K)*EXP(-PBM(K)) ENDDO C ***** C excitation spectrum of one-phonon resonance C ***** DO I=1,NEX SEC1(I)=0. SECSC1(I)=0. DO K=1,NRMAX SEC1(I)=SEC1(I)+INTR(K)*2.*PI*B(K)*PROB(K,I) SECSC1(I)=SECSC1(I)+INTR(K)*2.*PI*B(K)*PROBSC(K,I) G=6.-G ENDDO SEC1(I)=SEC1(I)*DELR/3./100 SECSC1(I)=SECSC1(I)*DELR/3./100 C ***** C end integration - add analytical part C ***** DO MU=0,L IF(MU.EQ.0) THEN J=1 ELSE J=2 ENDIF DUM=EX(I)*RMAX/GAMMA/HBC/BETA BKN=BESSK(ABS(MU),DUM) BKNP=BESSK(ABS(MU+1),DUM) BKNM=BESSK(ABS(MU-1),DUM) GM=PI*DUM**2*(BKNP*BKNM-BKN**2) SEC1(I)=SEC1(I)+J*COMFAC*ABS(GLM(L,MU,1./BETA))**2*GM* & EXP(-PB(NRMAX))*PHOTO(I)/EX(I)/100 SECSC1(I)=SECSC1(I)+J*COMFAC*ABS(GLM(L,MU,1./BETA))**2*GM* & EXP(-PBSC(NRMAX))*PHOTO(I)/EX(I)/100 ENDDO ENDDO C C-------------------------------------------- C contribution from double-phonon excitation C-------------------------------------------- C C excitation probability (1/MeV) as a function of energy (I) and C impact parameter (K) - harmonic oscillator C DO K=1,NRMAX DO I=2,NEX PROB2(K,I)=0. PROB2M(K,I)=0. PROSC2(K,I)=0. DO J=1,I-1 EE=EX(I)-EX(J) N=INT(EE/DELEX+1.) AN=1.-(EE-DELEX*N)/DELEX IF(N+1.LT.NEX.AND.N.GT.0) THEN PEE=AN*PROB(K,N)+(1.-AN)*PROB(K,N+1) PEEM=AN*PROBM(K,N)+(1.-AN)*PROBM(K,N+1) PEESC=AN*PROBSC(K,N)+(1.-AN)*PROBSC(K,N+1) ELSE PEE=0. PEEM=0. PEESC=0. ENDIF PROB2(K,I)=PROB2(K,I)+INTE(J)*PROB(K,J)*PEE PROB2M(K,I)=PROB2M(K,I)+INTE(J)*PROBM(K,J)*PEEM PROSC2(K,I)=PROSC2(K,I)+INTE(J)*PROBSC(K,J)*PEESC ENDDO PROB2(K,I)=0.5*PROB2(K,I)*EXP(PB(K))*DELEX/3. PROB2M(K,I)=0.5*PROB2M(K,I)*EXP(PBM(K))*DELEX/3. PROSC2(K,I)=0.5*PROSC2(K,I)*EXP(PBSC(K))*DELEX/3. ENDDO ENDDO C ***** C integration over impact parameter - Simpson's rule C ***** DO I=1,NEX SEC2(I)=0. SEC2M(I)=0. SECSC2(I)=0. DO K=1,NRMAX SEC2(I)=SEC2(I)+INTR(K)*2.*PI*B(K)*PROB2(K,I) SEC2M(I)=SEC2M(I)+INTR(K)*2.*PI*BM(K)*PROB2M(K,I) SECSC2(I)=SECSC2(I)+INTR(K)*2.*PI*B(K)*PROSC2(K,I) ENDDO SEC2(I)=(SEC2(I)+SEC2M(I))*DELR/3./100 SECSC2(I)=(SECSC2(I)+SEC2M(I))*DELR/3./100 SEC(I)=SEC1(I)+SEC2(I) SECSC(I)=SECSC1(I)+SECSC2(I) ENDDO C C-------------------------------------------------------------------- C Excitation spectrum and total Coulomb fragmentation cross section C-------------------------------------------------------------------- C C Total cross section C SUM1=0. SUM2=0. SUM=0. SUMSC1=0. SUMSC2=0. SUMSC=0. DO I=1,NEX SUM1=SUM1+INTE(I)*SEC1(I) SUM2=SUM2+INTE(I)*SEC2(I) SUMSC1=SUMSC1+INTE(I)*SECSC1(I) SUMSC2=SUMSC2+INTE(I)*SECSC2(I) ENDDO SUM1=DELEX*SUM1/3. SUM2=DELEX*SUM2/3. SUM=SUM1+SUM2 SUMSC1=DELEX*SUMSC1/3. SUMSC2=DELEX*SUMSC2/3. SUMSC=SUMSC1+SUMSC2 C ***** C Output C ***** C WRITE(11,*)' Energy (MeV),',' Excitation spectrum (b/MeV),', &' 1-phonon, 2-phonon, total' DO I=1,NEX WRITE(*,125)EX(I),SEC1(I),SEC2(I),SEC(I) WRITE(11,125)EX(I),SEC1(I),SEC2(I),SEC(I) 125 FORMAT(1X,4(G10.4,1x)) C WRITE(12,*)EX(I) C WRITE(13,*)SEC2(I) WRITE(14,125)EX(I),SECSC1(I),SECSC2(I),SECSC(I) C WRITE(12,*)EX(I) C WRITE(13,*)SECSC2(I) ENDDO C WRITE(11,*)' Total cross section (b), 1-phonon, 2-phonon, total' WRITE(*,126)SUM1,SUM2,SUM WRITE(11,126)SUM1,SUM2,SUM 126 FORMAT(1X,G8.3) PRINT*,' sharp-cutoff result' WRITE(*,126)SUMSC1,SUMSC2,SUMSC WRITE(14,126)SUMSC1,SUMSC2,SUMSC RETURN END C C********************* End of Main Program **************************** C C*********************************************************************** C C Fermi function C C*********************************************************************** C REAL*8 FUNCTION FERMI(R,R0,A) IMPLICIT NONE REAL*8 R,R0,A,TEMP TEMP=(R-R0)/A IF(TEMP.LT.30.) THEN FERMI=1./(1.+EXP(TEMP)) ELSE FERMI=0. ENDIF RETURN END C C*********************************************************************** C Droplet densities for neutrons - From Myers and Swiatecki, Ann. of C Phys. A55 (1969) 395 and C Ann. of Phys. A84 (1974) 186 C C*********************************************************************** C REAL*8 FUNCTION RHONP(RG,AP,ZP) REAL*8 RG, TNP, BNP, EPSP, DELTAP, RP, DP, RNP, CNP, RHO0NP REAL*8 PI/3.141597265/ REAL*8 AP, ZP TNP = 2.4 BNP = 0.413 * TNP DELTAP = ((AP-2.*ZP)/AP + 8.076E-3*ZP*AP**(-.66666)) # / ( 1. + 4.871 * AP**(-.33333) ) EPSP = -0.1724*AP**(-.33333) + 3.051E-3 * ZP**2 * AP**(-1.33333) # + .4166666 * DELTAP**2 RP = 1.18 * AP**0.33333 * ( 1. + EPSP ) DP = 0.66666 * ((AP-2.*ZP)/AP - DELTAP ) * RP RNP = RP + (ZP/AP) * DP CNP = RNP * ( 1. - (BNP/RNP)**2 ) RHO0NP = 3. * ( AP - ZP )/ ( 12.5664 * CNP**3 * # ( 1. + PI**2 * TNP**2 / (19.36 * CNP**2))) RHONP = RHO0NP / ( 1. + EXP((RG-CNP)/(TNP/4.4)) ) RETURN END C C*********************************************************************** C Droplet densities for protons - See Myers and Swiatecki, Ann. of C Phys. A55 (1969) 395 and C Ann. of Phys. A84 (1974) 186 C C*********************************************************************** C REAL*8 FUNCTION RHOPP(RG,AP,ZP) REAL*8 RG, TPP, BPP, EPSP, DELTAP, RP, DP, RPP, CPP, RHO0PP REAL*8 PI/3.141597265/ REAL*8 AP, ZP TPP = 2.4 BPP = 0.413 * TPP DELTAP = ((AP-2.*ZP)/AP + 8.076E-3*ZP*AP**(-.66666)) # / ( 1. + 4.871 * AP**(-.33333) ) EPSP = -0.1724*AP**(-.33333) + 3.051E-3 * ZP**2 * AP**(-1.33333) # + .4166666 * DELTAP**2 RP = 1.18 * AP**0.33333 * ( 1. + EPSP ) DP = 0.66666 * ((AP-2.*ZP)/AP - DELTAP ) * RP RPP = RP - ((AP-ZP)/AP) * DP CPP = RPP * ( 1. - (BPP/RPP)**2 ) RHO0PP = 3. * ZP / ( 12.5664 * CPP**3 * # ( 1. + PI**2 * TPP**2 / (19.36 * CPP**2))) RHOPP = RHO0PP / ( 1. + EXP((RG-CPP)/(TPP/4.4)) ) RETURN END C C*********************************************************************** C C normalizes densities to unity C C*********************************************************************** C SUBROUTINE NORM(R,DELR,DENS,RHO0) IMPLICIT NONE INTEGER NRMAX,I,NEX,NINT PARAMETER(NRMAX=201,NEX=401,NINT=1001) REAL*8 R(NRMAX),PI,DELR,SUM,DENS(NRMAX),RHO0, & INTR(NRMAX),INTE(NEX),INTT(NINT) COMMON /INT/ INTR,INTE,INTT PI=ACOS(-1.) SUM=0. DO I=2,NRMAX-1 SUM=SUM+INTR(I)*R(I)**2*DENS(I) ENDDO SUM=SUM*DELR/3. RHO0=1./4./PI/SUM RETURN END C C********************************************************************** C C nuclear eikonal phase C C********************************************************************** C SUBROUTINE PHNUC(Y,WN,B,Z,PHASEN) IMPLICIT NONE INTEGER*4 NRMAX,L,I,K,NEX,NINT PARAMETER(NRMAX=201,NINT=1001,NEX=401) REAL*8 WN(NRMAX),AL,B(NRMAX),Z(NRMAX),RV(NRMAX), & HBC,ALPHA,BETA,DELR,INTR(NRMAX),INTE(NEX), & INTT(NINT) COMPLEX*16 Y,PHASEN(NRMAX),VNUCL COMMON/A1/HBC,ALPHA,BETA,DELR COMMON /INT/ INTR,INTE,INTT DO K=1,NRMAX PHASEN(K)=0. C ***** C integrate in z coordinate - Simpson's rule C ***** DO I=2,NRMAX-1 RV(I)=SQRT(B(K)**2+Z(I)**2) L=INT(RV(I)/DELR+1.) AL=1.-(RV(I)-DELR*L)/DELR IF(L+1.LT.NRMAX.AND.L.GT.0) THEN VNUCL=Y*(AL*WN(L)+(1.-AL)*WN(L+1)) ELSE VNUCL=0. ENDIF PHASEN(K)=PHASEN(K)+INTR(I)*VNUCL ENDDO PHASEN(K)=-2./HBC/BETA*DELR/3.*PHASEN(K) ENDDO RETURN END C C********************************************************************** C C Folding of densities or potentials C C********************************************************************** C SUBROUTINE TWOFOLD(FACT,R,RMAX,DELRR,DENS1,DENS2,V) IMPLICIT NONE INTEGER*4 I,K,J,NRMAX,NAUX,IFRAC PARAMETER(NRMAX=201,NAUX=101) REAL*8 FACT,DENS1(NRMAX),DENS2(NRMAX),R(NRMAX),AL,DELX,DELRR, & V(NRMAX),DELR,X(NAUX),H,RINT(NAUX),G,SUM1,SUM2,AUX,PI,RMAX PI=ACOS(-1.) DELR=RMAX/NAUX DELX=2./NAUX DO I=1,NAUX X(I)=-1.+I*DELX RINT(I)=I*DELR ENDDO DO I=1,NRMAX H=4. SUM2=0. DO J=2,NAUX-1 G=4. SUM1=0. DO K=2,NAUX-1 AUX=SQRT(ABS(R(I)**2+RINT(J)**2-2.*R(I)*RINT(J)*X(K))) IFRAC=INT(AUX/DELRR+1.) AL=1.-(AUX-DELRR*IFRAC)/DELRR IF(IFRAC+1.LT.NRMAX.AND.IFRAC.GT.0) & SUM1=SUM1+G*(AL*DENS1(IFRAC)+(1.-AL)*DENS1(IFRAC+1)) G=6.-G ENDDO IFRAC=INT(RINT(J)/DELRR) AL=1.-(RINT(J)-DELRR*IFRAC)/DELRR IF(IFRAC+1.LT.NRMAX.AND.IFRAC.GT.0) & SUM2=SUM2+H*RINT(J)**2*(AL*DENS2(IFRAC) & +(1.-AL)*DENS2(IFRAC+1))*DELX/3.*SUM1 H=6.-H ENDDO V(I)=2.*PI*FACT*SUM2/3.*DELR ENDDO RETURN END C C********************************************************************** C C Winther-Alder functions - See Appendix B of A319 (1979) 518 C C********************************************************************** C COMPLEX*16 FUNCTION GLM(L,M,X) INTEGER*4 L,M REAL*8 X,PI COMPLEX*16 Y,VEC(3,0:3) IF(L.GT.3)PAUSE 'Bad argument in function GLM' Y=DCMPLX(0.,1.) PI=ACOS(-1.) VEC(1,1)=1./3.*SQRT(8.*PI)*X VEC(1,0)=Y*4./3.*SQRT(PI*(X**2-1.)) VEC(2,2)=-2./5.*X*SQRT(PI*(X**2-1.)/6.) VEC(2,1)=Y*2./5.*SQRT(PI/6.)*(2.*X**2-1.) VEC(2,0)=2./5.*X*SQRT(PI*(X**2-1.)) VEC(3,3)=1./21.*SQRT(PI/5.)*X*(X**2-1.) VEC(3,2)=-Y/21.*SQRT(2.*PI*(X**2-1.)/15.)*(3.*X**2-1.) VEC(3,1)=-1./105.*X*SQRT(PI/3.)*(15.*X**2-11.) VEC(3,0)=Y*2./105.*SQRT(PI*(X**2-1.))*(5.*X**2-1.) GLM=VEC(L,M) RETURN END C C********************************************************************** C C Calculates sign, sqrt, factorials, etc. of integers and half int. C C********************************************************************** C SUBROUTINE GFV(IE) C C IV(N) = (-1)**N C SQ(N) = SQRT(N) C SQI(N) = 1./SQRT(N) C SQH(N) = SQRT(N+1/2) C SHI(N) = 1./SQRT(N+1/2) C FAK(N) = N! C FAD(N) = N!! C FI(N) = 1./N! C WF(N) = SQRT(N!) C WFI(N) = 1./SQRT(N!) C GM2(N) = GAMMA(N+1/2) C WG(N) = SQRT(GAMMA(N+1/2)) C WGI(N) = 1./SQRT(GAMMA(N+1/2)) C C----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) C REAL*8 FAK(0:130) COMMON /GFVMAX/ MAXGFV COMMON /GFVIV / IV(0:130) COMMON /GFVSQ / SQ(0:130) COMMON /GFVSQI/ SQI(0:130) COMMON /GFVSQH/ SQH(0:130) COMMON /GFVSHI/ SHI(0:130) COMMON /GFVFAK/ FAK COMMON /GFVFAD/ FAD(0:130) COMMON /GFVFI / FI(0:130) COMMON /GFVWF / WF(0:130) COMMON /GFVWFI/ WFI(0:130) COMMON /GFVGM2/ GM2(0:130) COMMON /GFVWG / WG(0:130) COMMON /GFVWGI/ WGI(0:130) COMMON /SWITCH/ RSW(100),ISW(100),IPR(100) C IF (IE.GT.130) STOP 'STOP IN GFV: IE LARGER 130 ' MAXGFV = IE IV(0) = +1 SQ(0) = 0. SQI(0) = 1.D30 SQH(0) = SQRT(0.5D0) SHI(0) = 1./SQH(0) FAK(0) = 1.D0 FI(0) = 1.D0 WF(0) = 1.D0 WFI(0) = 1.D0 GM2(0) = 1.7724538509 WG(0) = SQRT(GM2(0)) WGI(0) = 1./WG(0) DO 10 I=1,130 IV(I) = -IV(I-1) SQ(I) = SQRT(FLOAT(I)) SQH(I) = SQRT(I+0.5D0) SQI(I) = 1./SQ(I) 10 SHI(I) = 1./SQH(I) DO 11 I=1,IE FAK(I) = I*FAK(I-1) WF(I) = SQ(I)*WF(I-1) GM2(I) = (I-0.5D0)*GM2(I-1) WG(I) = SQH(I-1)*WG(I-1) FI(I) = 1.D0/FAK(I) WFI(I) = 1.D0/WF(I) 11 WGI(I) = 1.D0/WG(I) FAD(1) = 1.D0 DO 12 I=3,IE,2 FAD(I) = I*FAD(I-2) 12 CONTINUE C IF (IPR(4).EQ.1) &PRINT *,'*** GFV ************************************************' RETURN C-END-GFV END C C*********************************************************************** C C Modified Bessel function of order N C C*********************************************************************** C REAL*8 FUNCTION BESSK(N,X) IMPLICIT REAL*8 (A-H,O-Z) IF(N.LE.1)GO TO 100 TOX=2.0/X BKM=BESSK0(X) BK=BESSK1(X) DO 11 J=1,N-1 BKP=BKM+J*TOX*BK BKM=BK BK=BKP 11 CONTINUE BESSK=BK GO TO 200 100 IF(N.EQ.0)BESSK=BESSK0(X) IF(N.EQ.1)BESSK=BESSK1(X) 200 RETURN END C C*********************************************************************** C C Modified Bessel function of zeroth order C C*********************************************************************** C REAL*8 FUNCTION BESSK0(X) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 Y,P1,P2,P3,P4,P5,P6,P7, * Q1,Q2,Q3,Q4,Q5,Q6,Q7 DATA P1,P2,P3,P4,P5,P6,P7/-0.57721566D0,0.42278420D0,0.23069756D0, * 0.3488590D-1,0.262698D-2,0.10750D-3,0.74D-5/ DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7/1.25331414D0,-0.7832358D-1,0.2189568D-1, * -0.1062446D-1,0.587872D-2,-0.251540D-2,0.53208D-3/ IF (X.LE.2.0) THEN Y=X*X/4.0 BESSK0=(-DLOG(X/2.0)*BESSI0(X))+(P1+Y*(P2+Y*(P3+ * Y*(P4+Y*(P5+Y*(P6+Y*P7)))))) ELSE Y=(2.0/X) BESSK0=(DEXP(-X)/DSQRT(X))*(Q1+Y*(Q2+Y*(Q3+ * Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7)))))) ENDIF RETURN END C C********************************************************************** C C Modified Bessel function of first order C C********************************************************************** C REAL*8 FUNCTION BESSK1(X) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 Y,P1,P2,P3,P4,P5,P6,P7, * Q1,Q2,Q3,Q4,Q5,Q6,Q7 DATA P1,P2,P3,P4,P5,P6,P7/1.0D0,0.15443144D0,-0.67278579D0, * -0.18156897D0,-0.1919402D-1,-0.110404D-2,-0.4686D-4/ DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7/1.25331414D0,0.23498619D0,-0.3655620D-1, * 0.1504268D-1,-0.780353D-2,0.325614D-2,-0.68245D-3/ IF (X.LE.2.0) THEN Y=X*X/4.0 BESSK1=(DLOG(X/2.0)*BESSI1(X))+(1.0/X)*(P1+Y*(P2+ * Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))) ELSE Y=2.0/X BESSK1=(DEXP(-X)/DSQRT(X))*(Q1+Y*(Q2+Y*(Q3+ * Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7)))))) ENDIF RETURN END C C*********************************************************************** C C Bessel function I of order 0 C C*********************************************************************** C REAL*8 FUNCTION BESSI0(X) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 Y,P1,P2,P3,P4,P5,P6,P7, * Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9 DATA P1,P2,P3,P4,P5,P6,P7/1.0D0,3.5156229D0,3.0899424D0,1.2067492D *0, * 0.2659732D0,0.360768D-1,0.45813D-2/ DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,0.1328592D-1, * 0.225319D-2,-0.157565D-2,0.916281D-2,-0.2057706D-1, * 0.2635537D-1,-0.1647633D-1,0.392377D-2/ IF (DABS(X).LT.3.75) THEN Y=(X/3.75)**2 BESSI0=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))) ELSE AX=DABS(X) Y=3.75/AX BESSI0=(DEXP(AX)/DSQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4 * +Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9)))))))) ENDIF RETURN END C C*********************************************************************** C C Bessel function I of order 1 C C*********************************************************************** C REAL*8 FUNCTION BESSI1(X) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 Y,P1,P2,P3,P4,P5,P6,P7, * Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9 DATA P1,P2,P3,P4,P5,P6,P7/0.5D0,0.87890594D0,0.51498869D0, * 0.15084934D0,0.2658733D-1,0.301532D-2,0.32411D-3/ DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,-0.3988024D-1, * -0.362018D-2,0.163801D-2,-0.1031555D-1,0.2282967D-1, * -0.2895312D-1,0.1787654D-1,-0.420059D-2/ IF (DABS(X).LT.3.75) THEN Y=(X/3.75)**2 BESSI1=X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))) ELSE AX=DABS(X) Y=3.75/AX BESSI1=(DEXP(AX)/DSQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+ * Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9)))))))) ENDIF RETURN END C C C**************************************************************************** C************ Explaining the inputs of program COULEIK *************** C ***** The input program is denoted COULEIK.IN ******* C ***** Units are in MeV and fm ***** C Dummy entries should be given in all rows, even if they will not be used. C An example is given for bismuth+gold at 1 GeV/nucleon C------------------------------------------------------------------------------ C C 1st row: AP,ZP,AT,ZT,ECA,RMAX,IFOLD,BMIN C - AP,ZP,AT,ZT are the charge and mass numbers of the nuclei C - ECA is the laboratory energy per nucleon C - RMAX is the upper limit of the integration over the radial coordinate. C Since a part of the cross section is obtained analytically for C large r's, RMAX should be given in the range between 20 and 30 fm C - IFOLD is an option for constructing the optical potential: C - IFOLD=2 for a "rho-rho" folding potential of a Gaussian C with a modified Fermi density (target) C - IFOLD=3 for a "rho-rho" folding potential of two Gaussian C densities C - IFOLD=4 for a "rho-rho" folding potential of two modified C Fermi densities C - IFOLD=5 for a "rho-rho" folding potential of a mofified C Gaussian (target) with a Gaussian density C - IFOLD=6 for a "rho-rho" folding potential of a modified C Gaussian with a modified Fermi density (target) C - IFOLD=7 for a "rho-rho" folding potential of two modified C Gaussian densities C - IFOLD=8 Use droplet densities for neutrons and protons C (From Myers and Swiatecki, ref. given in prog.) C C 136. 54. 136. 54. 600. 4 2 C - BMIN= minimum impact parameter to be used with a sharp-cutoff C------------------------ C 2th row: RGV1,RGV2,RRHO,ARHO,CRHO. The first two parameters, RGV1 and RGV2, C are used for the Gaussian densities C Rho1(r)=Rho0 . exp(-r^2/RGV1^2) C Rho2(r)=Rho0 . exp(-r^2/RGV2^2) C In the folding of a single gaussian density with another density C the function Rho1(r) is used. C The last three parameters are needed for the modified Fermi density C Rho(r)= Rho0/(1+exp((r-RRHO)/ARHO)) . (1+CRHO . r^2/RRHO^2) C or the modified Gaussian density C Rho(r)=Rho0/(1+exp((r^2-RRHO^2)/ARHO^2)) . (1+CRHO . r^2/RRHO^2) C The quantity Rho0 is not needed since the density is normalized C to the particle number. C 1.94 1.94 6.6 0.66 0.545 C----------------------- C 3th row: RRHO1,ARHO1,CRHO1,RRHO2,ARHO2,CRHO2 are needed for the "t-rho-rho" C folding of two modified Fermi or two modified Gaussian densities. C 5. 0.7 0.1 5. 0.7 0.1 C----------------------- C 4th row: SIGNN is the nucleon-nucleon cross sections (in fm^2) C 4.1 C----------------------- C 5th row: IW,L,ISO are options for the excitation mode C IW= 1 (2) for projectile (target) excitation C L= Angular momentum of the excitation C ISO= 0 (1) for a isoscalar (isovector) transition C 1 2 0 C----------------------- C 6th row:E0,WID,FRAC,EMAX are parameters for the (gamma,n) X-section, C in terms of a Lorentzian C E0= Centroid of the Lorentzian fit to (gamma, n) X-section C WID= Width of the Lorentzian fit C FRAC= Fraction of the sum rule exhausted by the (gamma, n) X-sect. C EMIN= Threshold for neutron emission C EMAX= Maximun photon energy. Use EMAX > 50 (MeV) C 15.3 4.8 1. 80. C*************************************************************************** C END OF PROGRAM C***************************************************************************