C C PROGRAM DEIKO.FOR C C---------------------------------------------------------------------- C Elastic and inelastic scattering cross sections for C proton-nucleus or heavy ion scattering at intermediate energies C---------------------------------------------------------------------- C written by C.A. Bertulani / 1994 C GSI, Planckstr. 1, D-64291 Darmstadt, Germany C Last modified version january/1999 C C For a description of the physics in the program, see (ref. 1): C C.A. Bertulani and H. Sagawa, Nucl. Phys. A588 (1995) 667 C C---------------------------------------------------------------------- C Input should be placed in file deiko.in. For a description C of how to make an input see the end of this program C C 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 NTMAX : no. of grid points in the scattering angle C RMAX : maximum value of radial coordinate used in the integrals C TETMAX : maximum value of the scattering angle C C Units are in MeV and fm. Cross sections in mb/sr. C---------------------------------------------------------------------- PROGRAM DEIKO CALL EWA END C SUBROUTINE EWA IMPLICIT NONE INTEGER*4 NRMAX,NTMAX,I,K,IOPT,L,ISO,MU,IFOLD,IFRAC, & ICOUL,IDUM,IW PARAMETER(NRMAX=500,NTMAX=200) INTEGER*4 ISR(NRMAX) REAL*8 PI,ALPHA,HBC,AMU,APMASS,ATMASS,PTMASS,PLAB,GAMMA,ECM, & BETA,ETA,EULER,DELR,DELT,R(NRMAX),B(NRMAX),TET(NTMAX),RP, & Z(NRMAX),V0,W0,RV0,RW0,AV,AW,VS0,RVS0,AVS,AI, & AP,ZP,AT,ZT,ECA,ELAB,VN(NRMAX),CPION,FERMI,EX,BJN, & WN(NRMAX),VS(NRMAX),PHASEC(NRMAX),PCM,ARED,BEL,RMS,RT,RTOT, & QT(NTMAX),SIGRUT(NTMAX),SEC(NTMAX),BESSJ0,BESSJ1,BESSJ, & BETAL,PL,H2M,QL,FAK(0:130),FACT,DELRC, & RRHO,RHO0,ARHO,RGV,RGW,WG0,VG0,SECC(NTMAX),SECN(NTMAX), & BC(NRMAX),RCMAX,AL,BCC(NRMAX),TEMP,DUM,DUM1,DUM2, & RMONO(NRMAX),RSCAL(NRMAX),RTRANS(NRMAX), & BESSK,BKN,ESQ,RHO01,ALPHAL, & RHO02,RRHO1,RRHO2,RHOG,RGV1,RGV2,RHOG1,RHOG2,SIGNN,ANN, & ARHO1,ARHO2,TETMAX,RMAX,CRHO,CRHO1,CRHO2,ACLOSE,BB(NRMAX), & DENS1(NRMAX),DENS2(NRMAX),VGN(NRMAX),WGN(NRMAX),VI0 COMPLEX*16 Y,FC(NTMAX),PHASEN(NRMAX),AMP(NTMAX),AMPS(NTMAX), & DERIV1(NRMAX),UMONO(NRMAX),ZRES(NRMAX),PHASC, & USCAL(NRMAX),UVEC(NRMAX),UTRANS(NRMAX),DERIV2(NRMAX), & AMPN(-5:5,NTMAX),PHASES(NRMAX),COMFAC(0:5),DUM3, & PHASE(NRMAX),AMPC(-5:5,NTMAX),GLM,UOPT(NRMAX) COMMON/A1/HBC,ALPHA,BETA,DELR COMMON/A2/EULER COMMON/A3/ISR COMMON /GFVFAK / FAK C OPEN(10, FILE='DEIKO.IN', STATUS='OLD') OPEN(11, FILE='DEIKO.OUT', STATUS='UNKNOWN') OPEN(12, FILE='DEIKO2.IN', STATUS='UNKNOWN') C CALL GFV(10) C C Inputs C READ(10,*)AP,ZP,AT,ZT,ECA,RMAX,TETMAX,IOPT,IFOLD READ(10,*)V0,RV0,AV,W0,RW0,AW READ(10,*)VG0,WG0,RGV,RGW READ(10,*)RGV1,RGV2,RRHO,ARHO,CRHO READ(10,*)RRHO1,ARHO1,CRHO1,RRHO2,ARHO2,CRHO2 READ(10,*)SIGNN,ANN READ(10,*)VI0,AI,VS0,RVS0,AVS READ(10,*)IW,L,ISO,ICOUL READ(10,*)EX,BEL,ALPHAL WRITE(11,'('' **** This is program DEIKO (1999) ***'')') C WRITE(11,1)AP,ZP,AT,ZT,IOPT,IFOLD,ECA 1 FORMAT(1X,/,1X,' Nuclear mass and charge numbers',/,2X, 1 4F5.1,/,2X, 2 ' Option for the reaction system:',/,2X, 3 I2,/,2X, 4 ' Option for the optical potentials',/,2X, 5 I2,/,2X, 6 ' Laboratory energy/per nucleon',/,2X, 7 F10.3) C C Constants, nuclear sizes, energies, etc. ELAB=ECA*AP Y=DCMPLX(0.,1.) PI=ACOS(-1.) ALPHA=1./137. HBC=197.33 H2M=20.73 AMU=931.5 CPION=1.41 ESQ=1.44 APMASS=AMU*AP ATMASS=AMU*AT PTMASS=(AP+AT)*AMU PLAB=SQRT(ELAB*(ELAB+2.*APMASS)) BETA=PLAB/(ELAB+APMASS) GAMMA=1./SQRT((1.-BETA)*(1.+BETA)) PLAB=PLAB/HBC ECM=SQRT(ATMASS**2+APMASS**2+2.*(ELAB+APMASS)*ATMASS) PCM=PLAB*ATMASS/ECM ETA=ZP*ZT*ALPHA/BETA ARED=APMASS*ATMASS/PTMASS ACLOSE=ETA/PCM EULER=0.5772156649 RP=1.2*AP**(1./3.) RT=1.2*AT**(1./3.) RTOT=RP+RT C grid in R,B,Z and tetha RCMAX=DMAX1(3.D0*GAMMA*HBC*BETA/EX,100.D0) DELR=RMAX/NRMAX DELRC=RCMAX/NRMAX DO I=1,NRMAX R(I)=I*DELR Z(I)=I*DELR B(I)=I*DELR BC(I)=I*DELRC ENDDO C grid in TET DELT=TETMAX*PI/180./NTMAX DO I=1,NTMAX TET(I)=I*DELT ENDDO C coefficients for Simpson's integrations ISR(1)=1 ISR(NRMAX)=1 IDUM=4 DO I=2,NRMAX-1 ISR(I)=IDUM IDUM=6-IDUM ENDDO IDUM=4 C----------------------------------------------------------------- C Calculate optical potential (depends on option IFOLD) C----------------------------------------------------------------- WRITE(11,'('' Parameters of the optical potential '')') GO TO (10,11,12,13,14,15,16,17,18,19),IFOLD+1 10 WRITE(11,110)V0,RV0,AV,W0,RW0,AW 110 FORMAT(1X,' V0,RV0,AV,W0,RW0,AW',/,1X,8(F8.3,X)) DO I=1,NRMAX VN(I)=V0*FERMI(R(I),RV0,AV) WN(I)=W0*FERMI(R(I),RW0,AW) ENDDO GO TO 20 11 WRITE(11,111)VG0,RGV,WG0,RGW,RRHO,ARHO,CRHO 111 FORMAT(' VG0,RGV,WG0,RGW,RRHO,ARHO,CRHO',/,1X,7(F8.3,X)) DO I=1,NRMAX VGN(I)=VG0*EXP(-(R(I)/RGV)**2) WGN(I)=WG0*EXP(-(R(I)/RGW)**2) DENS1(I)=FERMI(R(I),RRHO,ARHO)*(1.+CRHO*R(I)**2/RRHO**2) ENDDO CALL NORM(R,DELR,DENS1,RHO0) RHO0=RHO0*AT CALL TWOFOLD(RHO0,R,RMAX,DELR,VGN,DENS1,VN) CALL TWOFOLD(RHO0,R,RMAX,DELR,WGN,DENS1,WN) GO TO 20 12 WRITE(11,112)SIGNN,ANN,RGV1,RRHO,ARHO,CRHO 112 FORMAT(' SIGNN,ANN,RGV1,RRHO,ARHO,CRHO',/,1X,6(F8.3,X)) RHOG=AP/(PI**1.5*RGV1**3) VG0=-ANN/2.*HBC*BETA*SIGNN*RHOG WG0=-1./2.*HBC*BETA*SIGNN*RHOG DO I=1,NRMAX VGN(I)=VG0*EXP(-(R(I)/RGV1)**2) 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,VGN,DENS1,VN) CALL TWOFOLD(RHO0,R,RMAX,DELR,WGN,DENS1,WN) GO TO 20 13 WRITE(11,113)SIGNN,ANN,RGV1,RGV2 113 FORMAT(' SIGNN,ANN,RGV1,RGV2',/,1X,4(F8.3,X)) 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 VN(I)=-FACT*ANN*EXP(-R(I)**2/(RGV1**2+RGV2**2)) WN(I)=-FACT*EXP(-R(I)**2/(RGV1**2+RGV2**2)) ENDDO GO TO 20 14 WRITE(11,114)SIGNN,ANN,RRHO1,ARHO1,CRHO1,RRHO2,ARHO2,CRHO2 114 FORMAT(' SIGNN,ANN,RRHO1,ARHO1,CRHO1,RRHO2,ARHO2,CRHO2', 1/,1X,8(F8.3,X)) 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*ANN*HBC*BETA/2.*RHO01*RHO02 CALL TWOFOLD(FACT,R,RMAX,DELR,DENS1,DENS2,VN) FACT=-SIGNN*HBC*BETA/2.*RHO01*RHO02 CALL TWOFOLD(FACT,R,RMAX,DELR,DENS1,DENS2,WN) GO TO 20 15 WRITE(11,115)SIGNN,ANN,RGV1,RRHO,ARHO,CRHO 115 FORMAT(' SIGNN,ANN,RGV1,RRHO,ARHO,CRHO', 1/,1X,6(F8.3,X)) RHOG=AP/(PI**1.5*RGV1**3) VG0=-ANN/2.*HBC*BETA*SIGNN*RHOG WG0=-1./2.*HBC*BETA*SIGNN*RHOG DO I=1,NRMAX VGN(I)=VG0*EXP(-(R(I)/RGV1)**2) 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,VGN,DENS1,VN) CALL TWOFOLD(RHO0,R,RMAX,DELR,WGN,DENS1,WN) GO TO 20 16 WRITE(11,116)SIGNN,ANN,RRHO1,ARHO1,CRHO1,RRHO2,ARHO2,CRHO2 116 FORMAT(' SIGNN,ANN,RRHO1,ARHO1,CRHO1,RRHO2,ARHO2,CRHO2', 1/,1X,8(F8.3,X)) 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*ANN*HBC*BETA/2.*RHO01*RHO02 CALL TWOFOLD(FACT,R,RMAX,DELR,DENS1,DENS2,VN) FACT=-SIGNN*HBC*BETA/2.*RHO01*RHO02 CALL TWOFOLD(FACT,R,RMAX,DELR,DENS1,DENS2,WN) GO TO 20 17 WRITE(11,117)SIGNN,ANN,RRHO1,ARHO1,CRHO1,RRHO2,ARHO2,CRHO2 117 FORMAT(' SIGNN,ANN,RRHO1,ARHO1,CRHO1,RRHO2,ARHO2,CRHO2', 1/,1X,8(F8.3,X)) 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*ANN*HBC*BETA/2.*RHO01*RHO02 CALL TWOFOLD(FACT,R,RMAX,DELR,DENS1,DENS2,VN) FACT=-SIGNN*HBC*BETA/2.*RHO01*RHO02 CALL TWOFOLD(FACT,R,RMAX,DELR,DENS1,DENS2,WN) GO TO 20 18 WRITE(11,118)V0,RV0,AV,W0,RW0,AW,VI0,AI,VS0,RVS0,AVS 118 FORMAT(' V0,RV0,AV,W0,RW0,AW,VI0,AI,VS0,RVS0,AVS', 1/,1X,12(F8.3,X)) DO I=1,NRMAX VN(I)=V0*FERMI(R(I),RV0,AV) WN(I)=W0*FERMI(R(I),RW0,AV) & +4.*AI*VI0/AW*EXP((R(I)-RW0)/AW)* & FERMI(R(I),RW0,AW)**2 VS(I)=CPION**2/R(I)*VS0/AVS*EXP((R(I)-RVS0)/AVS)* & FERMI(R(I),RVS0,AVS)**2 ENDDO GO TO 20 19 WRITE(11,118)SIGNN,ANN 119 FORMAT(' SIGNN,ANN',/,1X,2(F8.3,X)) DO I=1,NRMAX READ(12,*)DUM,DENS1(I),DENS2(I) ENDDO CALL NORM(R,DELR,DENS1,RHO01) RHO01=AP*RHO01 CALL NORM(R,DELR,DENS2,RHO02) RHO02=AT*RHO02 FACT=-SIGNN*ANN*HBC*BETA/2.*RHO01*RHO02 CALL TWOFOLD(FACT,R,RMAX,DELR,DENS1,DENS2,VN) FACT=-SIGNN*HBC*BETA/2.*RHO01*RHO02 CALL TWOFOLD(FACT,R,RMAX,DELR,DENS1,DENS2,WN) C----------------------------------------------------------------- C Now, index T is reserved for the excited nucleus (depends on option IW) IF(IW.EQ.1) THEN DUM=AP TEMP=ZP AP=AT ZP=ZT AT=DUM ZT=TEMP ENDIF C---------------------------------------------------------------- C Nuclear and Coulomb phases C---------------------------------------------------------------- C correcting for Coulomb deflection 20 DO I=1,NRMAX BB(I)=ACLOSE+SQRT(ACLOSE**2+B(I)**2) BCC(I)=ACLOSE+SQRT(ACLOSE**2+BC(I)**2) ENDDO c initialize nuclear phase subroutine IF(IFOLD.NE.8) THEN DO I=1,NRMAX PHASES(I)=0. ENDDO CALL PHNUC(Y,VN,WN,BB,Z,PHASEN) ELSE CALL PHNUC(Y,VN,WN,BB,Z,PHASEN) DO I=1,NRMAX WN(I)=CPION**2/R(I)*VS0/AVS*EXP((R(I)-RVS0)/AVS)* & FERMI(R(I),RVS0,AVS)**2 ENDDO CALL PHNUC(Y,VS,WN,BB,Z,PHASES) ENDIF C Coulomb eikonal phase DO I=1,NRMAX PHASEC(I)=2.*ETA*DLOG(PCM*B(I)) C Total phase PHASE(I)=PHASEN(I)+PHASEC(I) ENDDO IF(IOPT.EQ.2)GO TO 50 C---------------------------------------------------------------- C this part is for elastic scattering amplitude and x-sections C---------------------------------------------------------------- C initialize Coulomb amplitude suboutine CALL FCOUL(Y,PCM,ETA,TET,FC) C elastic amplitudes DO I=1,NTMAX QT(I)=2.*PCM*SIN(TET(I)/2.) SIGRUT(I)=ABS(FC(I))**2 AMP(I)=0. AMPS(I)=0. C impact parameter integration - Simpson's rule DO K=1,NRMAX DUM1=QT(I)*B(K) DUM3=PCM*BB(K)*PHASES(K) AMP(I)=AMP(I)+ISR(K)*B(K)*BESSJ0(DUM1)*EXP(Y*PHASEC(K)) & *(1.-EXP(Y*PHASEN(K))*CDCOS(DUM3)) AMPS(I)=AMPS(I)+ISR(K)*B(K)*BESSJ1(DUM1)*EXP(Y*PHASEC(K)) & *EXP(Y*PHASEN(K))*CDSIN(DUM3) ENDDO AMP(I)=Y*PCM*AMP(I)*DELR/3. AMPS(I)=-Y*PCM*AMPS(I)*DELR/3. ENDDO c end integration WRITE(11,'(/,/,'' Angle (degrees) Ratio-to-Rutherford 1 cross section (mb/sr)'')') DO I=1,NTMAX SEC(I)=(ABS(AMP(I)+FC(I))**2+ABS(AMPS(I))**2)/SIGRUT(I) WRITE(11,150)TET(I)*180./PI,SEC(I),SEC(I)*SIGRUT(I)*10 WRITE(*,150)TET(I)*180./PI,SEC(I),SEC(I)*SIGRUT(I)*10 150 FORMAT(1X,3(G8.3,x)) ENDDO GO TO 60 C----------------------------------------------------------------- C this part is for nuclear inelastic amplitudes C----------------------------------------------------------------- C 50 RMS=0.6*RT**2 QL=EX/HBC/BETA PL=PCM-QL c PRINT*,' BEL=',BEL/ESQ,' e^2 fm^2',' ALPHAL=',ALPHAL, ' fm' WRITE(11,'('' Parameters for the inelastic scattering '')') WRITE(11,51)EX,L,ISO,ICOUL,BEL,ALPHAL 51 FORMAT(' EX,L,ISO,ICOUL,BEL,ALPHAL', 1/,1X,F8.3,X,X, 3(I3,X),3(F9.3,X)) C IF(ICOUL.EQ.1)GO TO 200 C------------------------------ C deformed potential model C------------------------------ C 1st AND 2nd derivatives of the nuclear potential - 6 point formulas DO I=1,NRMAX UOPT(I)=VN(I)+Y*WN(I) ENDDO CALL Derivative(UOPT,DERIV1,DERIV2,DELR,NRMAX) C isoscalar and isovector potentials DO I=1,NRMAX UMONO(I)=3.*UOPT(I)+DERIV1(I)*R(I) USCAL(I)=DERIV1(I) UVEC(I)=DERIV1(I)+RT*DERIV2(I)/3. ENDDO DO I=1,NRMAX IF(L.EQ.0.AND.ISO.EQ.0)UTRANS(I)=UMONO(I) IF(L.EQ.0.AND.ISO.EQ.1)UTRANS(I)=(AT-2.*ZT)/AT*UMONO(I) IF(L.GT.0.AND.ISO.EQ.0)UTRANS(I)=USCAL(I) IF(L.EQ.1.AND.ISO.EQ.1)UTRANS(I)=UVEC(I) ENDDO C excitation amplitude - index MU - deformed potential model DO MU=0,L C OBS: Eq. (39) of ref. 1 has a wrong factor sqrt((2*l+1)/(l+1)) COMFAC(MU)=-Y**MU*ARED/HBC/HBC*ALPHAL & *SQRT(FAK(L-MU)/FAK(L+MU)/4./PI) CALL ZINT(L,MU,DELR,QL,B,RMAX,UTRANS,ZRES) DO I=1,NTMAX AMPN(MU,I)=0. QT(I)=2.*SQRT(PCM*PL)*SIN(TET(I)/2.) C integration over impact parameter - Simpson's rule AMP(I)=0. DO K=1,NRMAX DUM1=QT(I)*BB(K) BJN=BESSJ(MU,DUM1) AMP(I)=AMP(I)+ISR(K)*B(K)*BJN*EXP(Y*PHASE(K))*ZRES(K) ENDDO AMPN(MU,I)=COMFAC(MU)*AMP(I)*DELR/3. ENDDO ENDDO C end integration DO MU=-L,-1,1 DO I=1,NTMAX AMPN(MU,I)=(-1)**MU*AMPN(-MU,I) ENDDO ENDDO 105 IF(ICOUL.EQ.0.OR.L.EQ.0)GO TO 210 C--------------------------------------------------- C Coulomb excitation amplitude C--------------------------------------------------- 200 DO MU=0,L COMFAC(MU)=ZP*SQRT(ESQ)*SQRT(PCM*PL)/GAMMA/HBC/BETA*Y**MU* & (EX/HBC)**L*SQRT(BEL)*GLM(L,MU,1./BETA) DO I=1,NTMAX AMPC(MU,I)=0. QT(I)=2.*SQRT(PL*PCM)*SIN(TET(I)/2.) C integration over impact parameter - Simpson's rule AMP(I)=0. DO K=1,NRMAX DUM1=QT(I)*BCC(K) DUM2=EX*BCC(K)/GAMMA/HBC/BETA BJN=BESSJ(MU,DUM1) BKN=BESSK(MU,DUM2) IFRAC=INT(BC(K)/DELR) AL=1.-(BC(K)-DELR*IFRAC)/DELR IF(IFRAC+1.LE.NRMAX.AND.IFRAC.GT.0) THEN PHASC=AL*PHASE(IFRAC)+(1.-AL)*PHASE(IFRAC+1) ELSE PHASC=2.*ETA*DLOG(PCM*BCC(K)) ENDIF AMP(I)=AMP(I)+ISR(K)*BC(K)*BJN*BKN*EXP(Y*PHASC) ENDDO AMPC(MU,I)=AMP(I)*COMFAC(MU)*DELRC/3. ENDDO ENDDO C end integration DO MU=-L,-1,1 DO I=1,NTMAX AMPC(MU,I)=(-1)**MU*AMPC(-MU,I) ENDDO ENDDO C Excitation cross sections 210 WRITE(11,*)' Angle (deg.),',' Inelastic cross section (mb/sr)' & ,' Nuclear, Coulomb, and Total' DO I=1,NTMAX SECN(I)=0. SECC(I)=0. SEC(I)=0. DO MU=-L,L SECN(I)=SECN(I)+PL/PCM*ABS(AMPN(MU,I))**2 SECC(I)=SECC(I)+PL/PCM*ABS(AMPC(MU,I))**2 SEC(I)=SEC(I)+PL/PCM*ABS(AMPN(MU,I)+AMPC(MU,I))**2 ENDDO WRITE(*,125)TET(I)*180/PI,SECN(I)*10,SECC(I)*10,SEC(I)*10 WRITE(11,125)TET(I)*180/PI,SECN(I)*10,SECC(I)*10,SEC(I)*10 125 FORMAT(1X,4(G8.3,X)) ENDDO 60 RETURN END C************ C Fermi function 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 normalizes densities to unity SUBROUTINE NORM(R,DELR,DENS,RHO0) IMPLICIT NONE INTEGER*4 NRMAX,I PARAMETER(NRMAX=500) REAL*8 R(NRMAX),PI,DELR,SUM,DENS(NRMAX),RHO0 INTEGER*4 ISR(NRMAX) COMMON/A3/ISR PI=ACOS(-1.) SUM=0. DO I=1,NRMAX SUM=SUM+ISR(I)*R(I)**2*DENS(I) ENDDO SUM=SUM*DELR/3. RHO0=1./4./PI/SUM RETURN END C************ C nuclear eikonal phase SUBROUTINE PHNUC(Y,VN,WN,B,Z,PHASEN) IMPLICIT NONE INTEGER*4 NRMAX,L,I,K PARAMETER(NRMAX=500) REAL*8 VN(NRMAX),WN(NRMAX),AL,B(NRMAX),Z(NRMAX),RV(NRMAX), & HBC,ALPHA,BETA,DELR COMPLEX*16 Y,PHASEN(NRMAX),VNUCL INTEGER*4 ISR(NRMAX) COMMON/A1/HBC,ALPHA,BETA,DELR COMMON/A3/ISR DO K=1,NRMAX PHASEN(K)=0. C integrate in z coordinate - Simpson's rule DO I=1,NRMAX RV(I)=SQRT(B(K)**2+Z(I)**2) L=INT(RV(I)/DELR) AL=1.-(RV(I)-DELR*L)/DELR IF(L+1.LT.NRMAX.AND.L.GT.0) THEN VNUCL=AL*VN(L)+(1.-AL)*VN(L+1)+Y*(AL*WN(L)+(1.-AL)*WN(L+1)) ELSE VNUCL=0. ENDIF PHASEN(K)=PHASEN(K)+ISR(I)*VNUCL ENDDO PHASEN(K)=-2./HBC/BETA*DELR/3.*PHASEN(K) ENDDO RETURN END C************* C z-integral for the inelastic scattering in the deformed model SUBROUTINE ZINT(L,MU,DELR,QL,B,RMAX,UTRANS,ZRES) IMPLICIT NONE INTEGER*4 NRMAX,NAUX,IFRAC,I,K,L,MU PARAMETER(NRMAX=500,NAUX=200) REAL*8 Z,AL,B(NRMAX),RV,QL,DELR,COSBET, & PLGNDR,SIGN,DELZ,RMAX INTEGER*4 ISR(NRMAX) COMPLEX*16 Y,AUXL,UTRANS(NRMAX),ZRES(NRMAX),FUNC COMMON/A3/ISR Y=DCMPLX(0.,1.) DELZ=RMAX/NAUX DO K=1,NRMAX ZRES(K)=0. C integrate in z coordinate - Simpson's rule DO I=1,NAUX Z=I*DELZ RV=SQRT(B(K)**2+Z**2) IFRAC=INT(RV/DELR) AL=1.-(RV-DELR*IFRAC)/DELR COSBET=Z/RV SIGN=(1.-(-1.)**L)/2. FUNC=((1.-SIGN)*COS(QL*Z)+SIGN*Y*SIN(QL*Z)) & *PLGNDR(L,MU,COSBET) IF(IFRAC+1.LT.NRMAX.AND.IFRAC.GT.0) & AUXL=(AL*UTRANS(IFRAC)+(1.-AL)*UTRANS(IFRAC+1))*FUNC ZRES(K)=ZRES(K)+ISR(I)*AUXL ENDDO ZRES(K)=2.*ZRES(K)*DELZ/3. ENDDO RETURN END C************ C Folding of densities or potentials SUBROUTINE TWOFOLD(FACT,R,RMAX,DELRR,DENS1,DENS2,V) IMPLICIT NONE INTEGER*4 I,K,J,NRMAX,NAUX,IFRAC PARAMETER(NRMAX=500,NAUX=100) REAL*8 FACT,DENS1(NRMAX),DENS2(NRMAX),R(NRMAX),AL,DELX,DELRR, & V(NRMAX),DELR,X(NAUX),RINT(NAUX),SUM1,SUM2,AUX,PI,RMAX INTEGER*4 ISR(NRMAX) COMMON/A3/ISR 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 SUM2=0. DO J=2,NAUX-1 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) AL=1.-(AUX-DELRR*IFRAC)/DELRR IF(IFRAC+1.LT.NRMAX.AND.IFRAC.GT.0) & SUM1=SUM1+ISR(K)*(AL*DENS1(IFRAC)+(1.-AL)*DENS1(IFRAC+1)) ENDDO IFRAC=INT(RINT(J)/DELRR) AL=1.-(RINT(J)-DELRR*IFRAC)/DELRR IF(IFRAC+1.LT.NRMAX.AND.IFRAC.GT.0) & SUM2=SUM2+ISR(J)*RINT(J)**2*(AL*DENS2(IFRAC) & +(1.-AL)*DENS2(IFRAC+1))*DELX/3.*SUM1 ENDDO V(I)=2.*PI*FACT*SUM2/3.*DELR ENDDO RETURN END C************ C elastic Coulomb amplitude SUBROUTINE FCOUL(Y,PCM,ETA,TET,FC) IMPLICIT NONE INTEGER*4 K,NEULER,NTMAX,I PARAMETER (NEULER=200,NTMAX=200) COMPLEX*16 Y, FC(NTMAX) REAL*8 ETA,TET(NTMAX),SUM,ARGAM,FACT,SQRUT,EULER,PCM COMMON/A2/EULER SUM=0. DO K=0,NEULER SUM=SUM+ETA/(1.+DFLOAT(K))-ATAN(ETA/(1.+DFLOAT(K))) ENDDO ARGAM=-EULER*ETA+SUM DO I=1,NTMAX FACT=-2.*ETA*DLOG(SIN(TET(I)/2.))+2.*ARGAM SQRUT=ETA/(2.*PCM*SIN(TET(I)/2.)**2) FC(I)=-EXP(Y*FACT)*SQRUT ENDDO RETURN END C************ C Winther-Alder functions. See appendix of Nucl. Phys. A319 (1971) 518 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 Spherical Bessel function of imaginary argument COMPLEX*16 FUNCTION BESSI(L,X) IMPLICIT NONE INTEGER*4 I,L COMPLEX*16 Y,AUX(0:5) REAL*8 X IF(L.GT.5)PAUSE 'Bad argument L in BESSI' IF(X.GT.20.)GO TO 10 Y=DCMPLX(0.,1.) AUX(0)=SINH(X)/X AUX(1)=Y/X**2*(X*COSH(X)-SINH(X)) AUX(2)=-(3./X**3+1./X)*SINH(X)+3./X**2*COSH(X) IF(L.GT.2) THEN DO I=2,L-1 AUX(I+1)=(2.*I+1.)*AUX(I)/(Y*X)-AUX(I-1) ENDDO ENDIF BESSI=AUX(L) 10 RETURN END C****************** C C Bessel function of zeroth order REAL*8 FUNCTION BESSJ0(X) IMPLICIT NONE REAL*8 Y,P1,P2,P3,P4,P5,Q1,Q2,Q3,Q4,Q5,R1,R2,R3,R4,R5,R6, * S1,S2,S3,S4,S5,S6,AX,Z,X,XX DATA P1,P2,P3,P4,P5/1.D0,-.1098628627D-2,.2734510407D-4, * -.2073370639D-5,.2093887211D-6/, Q1,Q2,Q3,Q4,Q5/-.1562499995D- *1, * .1430488765D-3,-.6911147651D-5,.7621095161D-6,-.934945152D-7/ DATA R1,R2,R3,R4,R5,R6/57568490574.D0,-13362590354.D0,651619640.7D *0, * -11214424.18D0,77392.33017D0,-184.9052456D0/, * S1,S2,S3,S4,S5,S6/57568490411.D0,1029532985.D0, * 9494680.718D0,59272.64853D0,267.8532712D0,1.D0/ IF(ABS(X).LT.8.)THEN Y=X**2 BESSJ0=(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))))) * /(S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))))) ELSE AX=ABS(X) Z=8./AX Y=Z**2 XX=AX-.785398164 BESSJ0=SQRT(.636619772/AX)*(COS(XX)*(P1+Y*(P2+Y*(P3+Y*(P4+Y * *P5))))-Z*SIN(XX)*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5))))) ENDIF RETURN END C************ C Bessel function of first order REAL*8 FUNCTION BESSJ1(X) IMPLICIT NONE REAL*8 Y,P1,P2,P3,P4,P5,Q1,Q2,Q3,Q4,Q5,R1,R2,R3,R4,R5,R6, * S1,S2,S3,S4,S5,S6,X,AX,Z,XX DATA R1,R2,R3,R4,R5,R6/72362614232.D0,-7895059235.D0,242396853.1D0 *, * -2972611.439D0,15704.48260D0,-30.16036606D0/, * S1,S2,S3,S4,S5,S6/144725228442.D0,2300535178.D0, * 18583304.74D0,99447.43394D0,376.9991397D0,1.D0/ DATA P1,P2,P3,P4,P5/1.D0,.183105D-2,-.3516396496D-4,.2457520174D-5 *, * -.240337019D-6/, Q1,Q2,Q3,Q4,Q5/.04687499995D0,-.2002690873D-3 *, * .8449199096D-5,-.88228987D-6,.105787412D-6/ IF(ABS(X).LT.8.)THEN Y=X**2 BESSJ1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))))) * /(S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))))) ELSE AX=ABS(X) Z=8./AX Y=Z**2 XX=AX-2.356194491 BESSJ1=SQRT(.636619772/AX)*(COS(XX)*(P1+Y*(P2+Y*(P3+Y*(P4+Y * *P5))))-Z*SIN(XX)*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5))))) * *SIGN(1.D0,X) ENDIF RETURN END C************** C Bessel function of order n REAL*8 FUNCTION BESSJ(N,X) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (IACC=40,BIGNO=1.E10,BIGNI=1.E-10) IF(N.LE.1)GO TO 100 TOX=2./X IF(X.GT.FLOAT(N))THEN BJM=BESSJ0(X) BJ=BESSJ1(X) DO 11 J=1,N-1 BJP=J*TOX*BJ-BJM BJM=BJ BJ=BJP 11 CONTINUE BESSJ=BJ ELSE M=2*((N+INT(SQRT(FLOAT(IACC*N))))/2) BESSJ=0. JSUM=0 SUM=0. BJP=0. BJ=1. DO 12 J=M,1,-1 BJM=J*TOX*BJ-BJP BJP=BJ BJ=BJM IF(ABS(BJ).GT.BIGNO)THEN BJ=BJ*BIGNI BJP=BJP*BIGNI BESSJ=BESSJ*BIGNI SUM=SUM*BIGNI ENDIF IF(JSUM.NE.0)SUM=SUM+BJ JSUM=1-JSUM IF(J.EQ.N)BESSJ=BJP 12 CONTINUE SUM=2.*SUM-BJ BESSJ=BESSJ/SUM ENDIF GO TO 200 100 IF(N.EQ.0)BESSJ=BESSJ0(X) IF(N.EQ.1)BESSJ=BESSJ1(X) 200 RETURN END C****************** C Legendre Polynomials REAL*8 FUNCTION PLGNDR(L,M,X) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) IF(M.LT.0.OR.M.GT.L.OR.ABS(X).GT.1.)PAUSE 'bad arguments' PMM=1. IF(M.GT.0) THEN SOMX2=SQRT((1.-X)*(1.+X)) FACT=1. DO 11 I=1,M PMM=-PMM*FACT*SOMX2 FACT=FACT+2. 11 CONTINUE ENDIF IF(L.EQ.M) THEN PLGNDR=PMM ELSE PMMP1=X*(2*M+1)*PMM IF(L.EQ.M+1) THEN PLGNDR=PMMP1 ELSE DO 12 LL=M+2,L PLL=(X*(2*LL-1)*PMMP1-(LL+M-1)*PMM)/(LL-M) PMM=PMMP1 PMMP1=PLL 12 CONTINUE PLGNDR=PLL ENDIF ENDIF RETURN END c c========================================================== c Calculates first and second derivatives of a function Wf c========================================================== c SUBROUTINE Derivative(Wf,Wf1,Wf2,Dx,NRMAX) IMPLICIT REAL*8 (a-h,o-z) COMPLEX*16 Wf(NRMAX),Wf1(NRMAX),Wf2(NRMAX) COMMON/coeff/ A(6,6),B(6,6),Ncount DATA Ncount/0/ IF(Ncount.gt.0) go to 1 Ncount=1 NPNTS=NRMAX ds=120.0d0*dx A(1,1)=-274.0d0/ds A(1,2)=600.0d0/ds A(1,3)=-600.0d0/ds A(1,4)=400.0d0/ds A(1,5)=-150.0d0/ds A(1,6)=24.0d0/ds A(2,1)=-24.0d0/ds A(2,2)=-130.0d0/ds A(2,3)=240.0d0/ds A(2,4)=-120.0d0/ds A(2,5)=40.0d0/ds A(2,6)=-6.0d0/ds A(3,1)=6.0d0/ds A(3,2)=-60.0d0/ds A(3,3)=-40.0d0/ds A(3,4)=120.0d0/ds A(3,5)=-30.0d0/ds A(3,6)=4.0d0/ds A(4,1)=-4.0d0/ds A(4,2)=30.0d0/ds A(4,3)=-120.0d0/ds A(4,4)=40.0d0/ds A(4,5)=60.0d0/ds A(4,6)=-6.0d0/ds A(5,1)=6.0d0/ds A(5,2)=-40.0d0/ds A(5,3)=120.0d0/ds A(5,4)=-240.0d0/ds A(5,5)=130.0d0/ds A(5,6)=24.0d0/ds A(6,1)=-24.0d0/ds A(6,2)=150.0d0/ds A(6,3)=-400.0d0/ds A(6,4)=600.0d0/ds A(6,5)=-600.0d0/ds A(6,6)=274.0d0/ds ds=60.0d0*Dx**2 B(1,1)=225.0d0/ds B(1,2)=-770.0d0/ds B(1,3)=1070.0d0/ds B(1,4)=-780.0d0/ds B(1,5)=305.0d0/ds B(1,6)=-50.0d0/ds B(2,1)=50.0d0/ds B(2,2)=-75.0d0/ds B(2,3)=-20.0d0/ds B(2,4)=70.0d0/ds B(2,5)=-30.0d0/ds B(2,6)=5.0d0/ds B(3,1)=-5.0d0/ds B(3,2)=80.0d0/ds B(3,3)=-150.0d0/ds B(3,4)=80.0d0/ds B(3,5)=-5.0d0/ds B(3,6)=0.0d0 B(4,1)=0.0d0 B(4,2)=-5.0d0/ds B(4,3)=80.0d0/ds B(4,4)=-150.0d0/ds B(4,5)=80.0d0/ds B(4,6)=-5.0d0/ds B(5,1)=5.0d0/ds B(5,2)=-30.0d0/ds B(5,3)=70.0d0/ds B(5,4)=-20.0d0/ds B(5,5)=-75.0d0/ds B(5,6)=50.0d0/ds B(6,1)=-50.0d0/ds B(6,2)=305.0d0/ds B(6,3)=-780.0d0/ds B(6,4)=1070.0d0/ds B(6,5)=-770.0d0/ds B(6,6)=225.0d0/ds 1 CONTINUE DO 2 Ir=3,NPNTS-3 WF1(Ir)=A(3,1)*Wf(Ir-2)+A(3,2)*Wf(Ir-1)+A(3,3)*Wf(Ir) & +A(3,4)*Wf(Ir+1)+A(3,5)*Wf(Ir+2)+A(3,6)*Wf(Ir+3) WF2(Ir)=B(3,1)*Wf(Ir-2)+B(3,2)*Wf(Ir-1)+B(3,3)*Wf(Ir) & +B(3,4)*Wf(Ir+1)+B(3,5)*Wf(Ir+2)+B(3,6)*Wf(Ir+3) 2 CONTINUE Ir=1 WF1(IR)=A(1,1)*Wf(Ir)+A(1,2)*Wf(Ir+1)+A(1,3)*Wf(Ir+2) & +A(1,4)*Wf(Ir+3)+A(1,5)*Wf(Ir+4)+A(1,6)*Wf(Ir+5) WF2(IR)=B(1,1)*Wf(Ir)+B(1,2)*Wf(Ir+1)+B(1,3)*Wf(Ir+2) & +B(1,4)*Wf(Ir+3)+B(1,5)*Wf(Ir+4)+B(1,6)*Wf(Ir+5) Ir=2 WF1(IR)=A(2,1)*Wf(Ir-1)+A(2,2)*Wf(Ir)+A(2,3)*Wf(Ir+1) & +A(2,4)*Wf(Ir+2)+A(2,5)*Wf(Ir+3)+A(2,6)*Wf(Ir+4) WF2(IR)=B(2,1)*Wf(Ir-1)+B(2,2)*Wf(Ir)+B(2,3)*Wf(Ir+1) & +B(2,4)*Wf(Ir+2)+B(2,5)*Wf(Ir+3)+B(2,6)*Wf(Ir+4) Ir=NPNTS-2 WF1(IR)=A(4,1)*Wf(Ir-3)+A(4,2)*Wf(Ir-2)+A(4,3)*Wf(Ir-1) & +A(4,4)*Wf(Ir)+A(4,5)*Wf(Ir+1)+A(4,6)*Wf(Ir+2) WF2(IR)=B(4,1)*Wf(Ir-3)+B(4,2)*Wf(Ir-2)+B(4,3)*Wf(Ir-1) & +B(4,4)*Wf(Ir)+B(4,5)*Wf(Ir+1)+B(4,6)*Wf(Ir+2) Ir=NPNTS-1 WF1(IR)=A(5,1)*Wf(Ir-4)+A(5,2)*Wf(Ir-3)+A(5,3)*Wf(Ir-2) & +A(5,4)*Wf(Ir-1)+A(5,5)*Wf(Ir)+A(5,6)*Wf(Ir+1) WF2(IR)=B(5,1)*Wf(Ir-4)+B(5,2)*Wf(Ir-3)+B(5,3)*Wf(Ir-2) & +B(5,4)*Wf(Ir-1)+B(5,5)*Wf(Ir)+B(5,6)*Wf(Ir+1) Ir=NPNTS WF1(IR)=A(6,1)*Wf(Ir-5)+A(6,2)*Wf(Ir-4)+A(6,3)*Wf(Ir-3) & +A(6,4)*Wf(Ir-2)+A(6,5)*Wf(Ir-1)+A(6,6)*Wf(Ir) WF2(IR)=B(6,1)*Wf(Ir-5)+B(6,2)*Wf(Ir-4)+B(6,3)*Wf(Ir-3) & +B(6,4)*Wf(Ir-2)+B(6,5)*Wf(Ir-1)+B(6,6)*Wf(Ir) RETURN END c C**************************** C Calculates sign, sqrt, factorials, etc. of integers and half int. 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 REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) 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 Modified Bessel function REAL*8 FUNCTION BESSK(N,X) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) 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 Modified Bessel function of zeroth order REAL*8 FUNCTION BESSK0(X) REAL*8 Y,P1,P2,P3,P4,P5,P6,P7,X,BESSI0, * 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=(-LOG(X/2.0)*BESSI0(X))+(P1+Y*(P2+Y*(P3+ * Y*(P4+Y*(P5+Y*(P6+Y*P7)))))) ELSE Y=(2.0/X) BESSK0=(EXP(-X)/SQRT(X))*(Q1+Y*(Q2+Y*(Q3+ * Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7)))))) ENDIF RETURN END C**************** C Modified Bessel function of first order REAL*8 FUNCTION BESSK1(X) REAL*8 Y,P1,P2,P3,P4,P5,P6,P7,X,BESSI1, * 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=(LOG(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=(EXP(-X)/SQRT(X))*(Q1+Y*(Q2+Y*(Q3+ * Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7)))))) ENDIF RETURN END C******************* C Bessel function I of order 0 REAL*8 FUNCTION BESSI0(X) REAL*8 Y,P1,P2,P3,P4,P5,P6,P7,X,AX, * 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 (ABS(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=ABS(X) Y=3.75/AX BESSI0=(EXP(AX)/SQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4 * +Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9)))))))) ENDIF RETURN END C*************** C Bessel function I of order 1 REAL*8 FUNCTION BESSI1(X) REAL*8 Y,P1,P2,P3,P4,P5,P6,P7,AX,X, * 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 (ABS(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=ABS(X) Y=3.75/AX BESSI1=(EXP(AX)/SQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+ * Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9)))))))) ENDIF RETURN END C**************************************************************************** C* C************ Explaining the inputs of program DEIKO *************** C ***** The actual input program is denoted DEIKO.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------------------------------------------------------------------------------ C 1st row: AP,ZP,AT,ZT,ECA,RMAX,TETMAX,IOPT,IFOLD C - AP,ZP,AT,ZT are the charge and mass numbers of the nuclei C - ECA is the projectile laboratory energy/per nucleon C - RMAX is the upper limit of the integration over the radial coordinate. C Since the Coulomb part (elastic and inelastic) is obtained analytically C for large r's, RMAX should be given in the range between 20 and 30 fm C - TETMAX is the maximum scatt. angle in degrees C - IOPT is a general option: C - IOPT=1 for nucleus-nucleus or proton-nucleus elastic scattering C - IOPT=2 for inelastic scattering C - IFOLD is an option for constructing the optical potential: C - IFOLD=0 for a Woods-Saxon potential C - IFOLD=1 for a folding of a Gaussian nucleon-nucleus potential C with a modified Fermi density (target) - Very common C for the excitation of a heavy nucleus (lead, gold, etc) C by a light one (alpha's, carbon, etc.) C - IFOLD=2 for a "t-rho-rho" folding potential of a Gaussian C with a modified Fermi density (target) C - IFOLD=3 for a "t-rho-rho" folding potential of two Gaussian C densities C - IFOLD=4 for a "t-rho-rho" folding potential of two modified C Fermi densities C - IFOLD=5 for a "t-rho-rho" folding potential of a mofified C Gaussian (target) with a Gaussian density C - IFOLD=6 for a "t-rho-rho" folding potential of a modified C Gaussian with a modified Fermi density (target) C - IFOLD=7 for a "t-rho-rho" folding potential of two modified C Gaussian densities C - IFOLD=8 for the proton-nucleus potential parameters C - IFOLD=9 to enter your own ground state densities in file DEIKO2.IN. C # of points = NRMAX, 1st column = R, 2n column = DENS1, C 3rd column = DENS2. R must be in steps of RMAX/NRMAX. C 16. 8. 16. 8. 350. 25. 60. 1 1 C------------------------- C 2nd row: V0,RV0,AV,W0,RW0,AW are the parameters of the Woods-Saxon potential. C The potential is of the form: C U(r)=V0/(1+exp((r-RV0)/AV)) + i . W0/(1+exp((r-RW0)/AW)) C -80. 3.8 0.82 -42. 3.9 0.82 C------------------------- C 3rd row: VG0,WG0,RGV,RGW are the strengths and ranges of the nucleon-nucleus C gaussian potential. The potential is of the form: C U(r)=VG0 . exp(-r^2/RGV^2) + i . WG0 . exp(-r^2/RWG^2) C -36.4 -22.5 2. 2. C------------------------- C 4th 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 2. 2. 6.6 0.66 0.545 C-------------------------- C 5th 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 6th row: SIGNN, ANN are the nucleon-nucleon cross sections (in fm^2) and the C real-to-imaginary part of the nucleon-nucleon scattering amplitude C 4.1 0.86 C-------------------------- C 7th row: VI0,AI,VS0,RVS0,AVS are the parameters for the proton-nucleus C potential (with spin-orbit and surface terms): C U(r)=Woods-Saxon with entries of row 2 C + 4 . AI . (d/dr) VI0/(1.+exp((r-RW0)/AW)) (surface term) C + (hbar/2 m r)^2 . (d/dr) VS0/(1.+exp((r-RVS0)/AVS)) (spin-orb.) C -20.06 0.546 -4.26 1.245 1.16 C-------------------------- C 8th row: IW,L,ISO,ICOUL are options for the inelastic scattering 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 ICOUL=0 (1) [2] for nuclear (Coulomb) [nuclear+Coulomb] excitation C 2 2 0 2 C-------------------------- C 9th row:EX,BEL,ALPHAL are parameters for the inelastic excitation C part of the program C EX= Excitation energy of the state C BEL= B-value of the state. This entry is for the reduced C matrix element directly. C ALPHAL= Deformation length of the nuclear excitation. C 13.5 1. 0.1 0.1 C========== C NOTE: According to Sachtler, NPA472 (1987) 215, C U_int = - Alpha0*(3*U_opt+r*dU_opt/dr), (Monopole excitations) C where Alpha0^2 = 2*pi*(hbar^2/m)/(AE_x) C U_int = - Alphal*dU_opt/dr, (isoscalar, l>1, excitations) C where C Alphal^2 = l*(2*l+1)^2/(l+2)^2 * 2*pi*(hbar^2/m)/(AE_x) C * /^2 C U_int = - Alpha1*(dU_opt/dr+(R_0/3)*d^2U_opt/d^2r) C (isovector excitations, l=1), where C Alpha1^2 = pi*(hbar^2/(2*m))*A/(N*Z*E_x) C========== C******************