C************************************************************ C***** ***** C***** MAIN PROGRAM ***** C***** ***** C***** SCATTERING WAVE FUNCTION CALCULATED WITH ***** C***** EIKONAL APPROXIMATION ***** C***** ***** C************************************************************ C C======================================================================= PROGRAM WFEIKONAL02 IMPLICIT REAL*8(A-H,O-Z) C======================================================================= C PARAMETER (HC0=197.3269601D0,AC0=931.494013D0,EC0=137.03599976D0) PARAMETER(NBDIM=1001) C COMPLEX*16 SEK(NBDIM) C COMMON /PHYSCONST/HC,AC,EC COMMON /MACHINE/EPSLIM,EXPMAX CCC--------------------------------------------------------------------- HC=HC0 AC=AC0 EC=EC0 C C------------------------------------- C--- LIMIT VALUE FOR COMPUTATION --- C------------------------------------- C EXPMAX=300.0D0 EPSLIM=1.0D-14 C C C------------------- C--- FILE OPEN --- C------------------- C CALL FOPEN C C C------------------------- C--- PARAMETER INPUT --- C------------------------- C CALL INPUT(AP,AT,ZP,ZT,ELAB,RMAX,DR,NRMAX,V0,RV,AV 1 ,W0,RW,AW,WD,RD,AD,BMIN,BMAX,DB,NBMAX,KIBCS 2 ,THMIN,THMAX,DTH,NTHMAX,ICM,IREL) C C C----------------------------------- C--- SCATTERING WAVE FUNCTION --- C----------------------------------- C CALL WF(AP,AT,ZP,ZT,ELAB,RMAX,DR,NRMAX,V0,RV,AV 1 ,W0,RW,AW,WD,RD,AD,BMIN,BMAX,DB,NBMAX,NBDIM,SEK,FK 2 ,ICM,IREL) C C C----------------------- C--- CROSS SECTION --- C----------------------- C CALL CS(FK,SEK,BMIN,BMAX,DB,NBMAX,NBDIM,KIBCS 1 ,THMIN,THMAX,DTH,NTHMAX) C STOP 0 END C C C Cver0.2 C======================================================================= SUBROUTINE CS(FK,SEK,BMIN,BMAX,DB,NBMAX,NBDIM,KIBCS 1 ,THMIN,THMAX,DTH,NTHMAX) IMPLICIT REAL*8 (A-H,O-Z) C======================================================================= C COMPLEX*16 CSUM,SEK(NBDIM) CCC--------------------------------------------------------------------- C PI=ACOS(-1.0D0) PIAD=PI/180.0D0 C WRITE(KIBCS,*) ' THETA CS[mb/sr]' C DO ITH=1,NTHMAX THT=(ITH-1)*DTH+THMIN RAD=THT*PIAD C Q=2.0D0*FK*SIN(RAD/2.0D0) C CSUM=0.0D0 SUME=0.0D0 SUMR=0.0D0 SUMT=0.0D0 DO IB=1,NBMAX B=(IB-1)*DB+BMIN ARG=B*FK*RAD CSUM=CSUM+BESSJ0(ARG)*(SEK(IB)-1.0D0)*B SUME=SUME+ABS(1.0D0-SEK(IB))**2*B SUMR=SUMR+(1.0D0-ABS(SEK(IB))**2)*B SUMT=SUMT+2.0D0*(1.0D0-DBLE(SEK(IB)))*B END DO CSUM=CSUM*DB C SIG=FK**2*ABS(CSUM)**2*10.0D0 C WRITE(KIBCS,600) THT,SIG 600 FORMAT(F8.4,E13.5) C END DO C SIGE=SUME*2.0D0*PI*DB*10.0D0 SIGR=SUMR*2.0D0*PI*DB*10.0D0 SIGT=SUMT*2.0D0*PI*DB*10.0D0 C C-----OUTPUT C WRITE(6,601) SIGE,SIGR,SIGT 601 FORMAT(/3X,'-- INTEGRATED CROSS SECTION --' 1 /8X,'NUCLEAR ELASTIC :',F10.3,' [mb]' 2 /8X,'REACTION :',F10.3,' [mb]' 3 /8X,'TOTAL :',F10.3,' [mb]') C RETURN C END Cver0.2 C C C C======================================================================= SUBROUTINE WF(AP,AT,ZP,ZT,ELAB,RMAX,DR,NRMAX,V0,RV,AV 1 ,W0,RW,AW,WD,RD,AD,BMIN,BMAX,DB,NBMAX,NBDIM,SEK,FK 2 ,ICM,IREL) IMPLICIT REAL*8 (A-H,O-Z) C======================================================================= C COMPLEX*16 CSUM,WI COMPLEX*16 SEK(NBDIM) C COMMON /PHYSCONST/HC,AC,EC COMMON /MACHINE/EPSLIM,EXPMAX CCC--------------------------------------------------------------------- C C-----DIMENSION CHECK FOR B C IF(NBMAX.GT.NBDIM) GO TO 910 C C-----INITIALIZATION C SEK(:)=0.0D0 C WI=(0.0D0,1.0D0) RMAX=(NRMAX-1)*DR FM=AP*AC FK=SQRT(2.0D0*FM*ELAB)/HC HVEL=FK*HC**2/FM C C-----KINEMATICAL CORRECTION C IF(IREL.EQ.1) THEN FMP=AP*AC FMT=AT*AC EPL=ELAB+FMP ETL=FMT FKPL=SQRT(ELAB*(ELAB+2.0D0*FMP))/HC IF(ICM.EQ.0) THEN FK=FKPL HVEL=FKPL*HC**2/EPL ! invariant ELSE FK=FKPL*FMT/SQRT(FMP**2+FMT**2+2.0D0*FMT*EPL) HVEL=FKPL*HC**2/EPL ! invariant END IF ELSE IF(ICM.EQ.1) THEN ECM=ELAB*AT/(AP+AT) FM=AP*AT/(AP+AT)*AC FK=SQRT(2.0D0*FM*ECM)/HC HVEL=FK*HC**2/FM END IF C C-----OUTPUT C WRITE(6,600) FK,HVEL 600 FORMAT(/3X,'-- KINEMATICS --' 1 /8X,'FK (WAVE NUMBER) :',F10.4,' [1/fm]' 2 /8X,'HBAR*C :',F10.4,' [MeV fm]') C WRITE(6,601) 601 FORMAT(/3X,'-- S-MATRIX --' 1 /7X,'b realS imagS') C C DO IB=1,NBMAX B=(IB-1)*DB+BMIN IF(B.GT.RMAX) THEN NBMAX=IB-1 RETURN END IF C DZ=DR ZMAX=SQRT(RMAX**2-B**2) NZMAX0=ZMAX/DZ+1.01 NZMAX=2*(NZMAX0-1)+1 ZMAX=(NZMAX0-1)*DZ C CSUM=0.0D0 DO IZ=1,NZMAX Z=(IZ-1)*DZ-ZMAX R=SQRT(Z**2+B**2) ARGV=(R-RV)/AV IF(ARGV.GT.EXPMAX) THEN VOPT=0.0D0 ELSE VOPT=-V0/(1.0D0+EXP(ARGV)) END IF ARGW=(R-RW)/AW IF(ARGW.GT.EXPMAX) THEN WOPTV=0.0D0 ELSE WOPTV=-W0/(1.0D0+EXP(ARGW)) END IF ARGD=(R-RD)/AD IF(2.0D0*ARGD.GT.EXPMAX) THEN WOPTS=0.0D0 ELSE WOPTS=-4.0D0*WD*EXP(ARGD)/(1.0D0+EXP(ARGD))**2 END IF WOPT=WOPTV+WOPTS CSUM=CSUM+VOPT+WI*WOPT END DO C SEK(IB)=EXP(CSUM*DZ/(HVEL*WI)) C WRITE(6,602) B,SEK(IB) 602 FORMAT(F10.4,2E13.5) C END DO C RETURN C 910 WRITE(6,6910) NBMAX,NBDIM 6910 FORMAT(/3X,'-- DIMENSION ERROR IN WF --' 1 /8X,'NBMAX > NBDIM: NBMAX,NBDIM=',2I8) STOP C END C C C C======================================================================= SUBROUTINE INPUT(AP,AT,ZP,ZT,ELAB,RMAX,DR,NRMAX,V0,RV,AV 1 ,W0,RW,AW,WD,RD,AD,BMIN,BMAX,DB,NBMAX,KIBCS 2 ,THMIN,THMAX,DTH,NTHMAX,ICM,IREL) IMPLICIT REAL*8 (A-H,O-Z) C======================================================================= C CHARACTER*2 NAMEP,NAMET CHARACTER*50 COMMENT C COMMON /MACHINE/EPSLIM,EXPMAX CCC--------------------------------------------------------------------- C KIBAN=5 C AV=0.6D0 AW=0.6D0 AD=0.6D0 C READ(KIBAN,500) COMMENT 500 FORMAT(A50) READ(KIBAN,501) AP,ZP,AT,ZT 501 FORMAT(5F10.0) READ(KIBAN,501) ELAB,RMAX,DR READ(KIBAN,501) BMIN,BMAX,DB READ(KIBAN,501) FICM,FIREL READ(KIBAN,501) V0,RV,AV READ(KIBAN,501) W0,RW,AW READ(KIBAN,501) WD,RD,AD READ(KIBAN,500) DUMMY READ(KIBAN,501) FKIBCS,THMIN,THMAX,DTH C NRMAX=RMAX/DR+1.01 NBMAX=(BMAX-BMIN)/DB+1.01 BMAX=BMIN+(NBMAX-1)*DB C ICM=FICM IREL=FIREL KIBCS=FKIBCS NTHMAX=(THMAX-THMIN)/DTH+1.01 C IF(RV.LT.-EPSLIM) THEN RV=-RV ELSE RV=RV*AT**(1.0D0/3.0D0) END IF IF(RW.LT.-EPSLIM) THEN RW=-RW ELSE RW=RW*AT**(1.0D0/3.0D0) END IF IF(RD.LT.-EPSLIM) THEN RD=-RD ELSE RD=RD*AT**(1.0D0/3.0D0) END IF C C---------------------------------------- C--- NAMES OF PROJECTILE AND TARGET --- C---------------------------------------- C CALL ELEMNT(AP,ZP,MASSP,NAMEP) CALL ELEMNT(AT,ZT,MASST,NAMET) C C C------------------------------- C--- OUTPUT FOR THE INPUTS --- C------------------------------- C WRITE(6,600) 600 FORMAT(5X,'***************************************************' 1 ,'**************' 2 /5X,'***** ' 3 ,' *****' 4 /5X,'***** PROGRAM WFEK ' 5 ,' *****' 6 /5X,'***** ' 7 ,' *****' 8 /5X,'***************************************************' 9 ,'**************') C WRITE(6,601) MASSP,NAMEP,MASST,NAMET,NINT(ELAB) 601 FORMAT(//3X,'===================================================' 1 /8X,'REACTION:',3X,I3,A2,' + ',I3,A2,5X,'at ',I5,'MeV' 2 /3X,'===================================================') C WRITE(6,602) COMMENT 602 FORMAT(/5X,'USER''S COMMENT:'/10X,A50) C WRITE(6,603) AP,ZP,AT,ZT, ELAB,RMAX,DR, BMIN,BMAX,DB ,ICM,IREL 603 FORMAT(/3X,'-- GENERAL INPUTS --' 1 /8X,'AP ZP AT ZT :',2F10.5,2F10.2 2 /8X,'ELAB RMAX DR :',3F10.5 3 /8X,'BMIN BMAX DB :',3F10.5 4 /8X,'ICM IREL :',2I5) C IF(ABS(ZP*ZT).GT.EPSLIM) WRITE(6,604) 604 FORMAT(/8X,'NOTE: BOTH PARTICLES ARE CHARGED BUT NO COULOMB' 1 ,' INTERACTION IS INCLUDED') C WRITE(6,605) 605 FORMAT(/3X,'-- OPTICAL POTENTIAL PARAMETERS --') C IF(ABS(V0).GT.EPSLIM) THEN WRITE(6,606) V0,RV,AV ELSE WRITE(6,607) END IF 606 FORMAT(8X,'V0 RV AV :',3F10.5) 607 FORMAT(8X,'V0 RV AV : NO REAL POT.') IF(ABS(W0).GT.EPSLIM) THEN WRITE(6,608) W0,RW,AW ELSE WRITE(6,609) END IF 608 FORMAT(8X,'W0 RW AW :',3F10.5) 609 FORMAT(8X,'W0 RW AW : NO VOLUME IMAGINARY POT.') IF(ABS(W0).GT.EPSLIM) THEN WRITE(6,610) WD,RD,AD ELSE WRITE(6,611) END IF 610 FORMAT(8X,'WD RD AD :',3F10.5) 611 FORMAT(8X,'WD RD AD : NO SURFACE IMAGINARY POT.') C WRITE(6,612) KIBCS,THMIN,THMAX,DTH 612 FORMAT(/3X,'-- CONTROL OF OUTPUT --' 1 /8X,'KIBCS THMIN THMAX DTH :',I5,3F10.5) C RETURN C END C C C C======================================================================= FUNCTION BESSJ0(Z) IMPLICIT REAL*8 (A-H,O-Z) C======================================================================= C CCC--------------------------------------------------------------------- C IF(ABS(Z).LE.8.0D0) THEN P00= 0.5756849057376692D11 P01=-0.1336259035395191D11 P02= 0.6516196407500873D09 P03=-0.1121442418227296D08 P04= 0.7739233017240460D05 P05=-0.1849052456430540D03 C Q00= 0.5756849041067652D11 Q01= 0.1029532985766340D10 Q02= 0.9494680718778989D07 Q03= 0.5927264853423237D05 Q04= 0.2678532712126684D03 Q05= 1.0D0 C Y=Z**2 BESSJ0=(P00+Y*(P01+Y*(P02+Y*(P03+Y*(P04+Y*P05))))) 1 /(Q00+Y*(Q01+Y*(Q02+Y*(Q03+Y*(Q04+Y*Q05))))) C ELSE Y=(8.0D0/Z)**2 C P00= 0.9999999999947955D0 P01=-0.109863241367192D-2 P02= 0.2737589862975D-4 P03=-0.216060245228D-5 P04= 0.30994056326D-6 P05=-0.4048673907D-7 C P0=P00+Y*(P01+Y*(P02+Y*(P03+Y*(P04+Y*P05)))) C P00=-0.1562499999676469D-1 P01= 0.14305089740097D-3 P02=-0.692763042016D-5 P03= 0.80888732887D-6 P04=-0.14748826511D-6 P05= 0.2176282453D-7 C Q0=8.0D0/ABS(Z)*(P00+Y*(P01+Y*(P02+Y*(P03+Y*(P04+Y*P05))))) C PI=ACOS(-1.0D0) ARG=ABS(Z)-PI/4.0D0 BESSJ0=SQRT(2.0D0/PI/ABS(Z))*(P0*COS(ARG)-Q0*SIN(ARG)) END IF C RETURN END C C C C----------------------------------------------------------------------- C--- Subroutines below are made by Y. Iseri (Chiba-Keizai College) --- C----------------------------------------------------------------------- CCC CCC ***** FOPEN ***** CCC SUBROUTINE FOPEN IMPLICIT REAL*8(A-H,O-Z) CHARACTER FNAME*50,STA*8,COMM*60,OFF*1 C ========================================== C --- UNIT-5 --- C WRITE(*,100) C 100 FORMAT(//1x,'>>> Input File-Name of Unit 5 ==> ',$) C READ(*,'(A)') FNAME C OPEN(5,FILE=FNAME,STATUS='OLD') C --- OTHER UNITS --- READ(5,'(A)') COMM WRITE(*,200) COMM 200 FORMAT(/1x,'---- (Comment in File) ----' 1 /1x,' << ',A,' >>' 2 //1x,'---- ( Open Files ) ----' 3 //1x,' Unit Status File-Name') 1 READ(5,300) OFF,IUNIT,STA,FNAME 300 FORMAT(A1,I3,1X,A8,2X,A50) IF(OFF.NE.' ') GO TO 1 IF(IUNIT.GE.100 .OR. IUNIT.LT.0) GO TO 900 WRITE(*,310) IUNIT,STA,FNAME 310 FORMAT(1x,3X,I2,3X,A,3X,A) OPEN(IUNIT,FILE=FNAME,STATUS=STA) GO TO 1 C --- ONE MORE LINE FOR DISCRIMINATION --- 900 READ(5,'(A)') COMM RETURN END CCCCC CCCCC ***** ELEMNT ***** CCCCC SUBROUTINE ELEMNT(A,Z,MASS,NAME) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (MMAX=109) CHARACTER*2 NAME CHARACTER*2 NAT(MMAX),NAP(3),NAA(2),NEU,NON DATA NEU,NAT / 'n ', 1 'H ','He','Li','Be','B ','C ','N ','O ','F ','Ne','Na', 2 'Mg','Al','Si','P ','S ','Cl','Ar','K ','Ca','Sc','Ti', 3 'V ','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As', 4 'Se','Br','Kr','Rb','Sr','Y ','Zr','Nb','Mo','Tc','Ru', 5 'Rh','Pd','Ag','Cd','In','Sn','Sb','Te','I ','Xe','Cs', 6 'Ba','La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy', 7 'Ho','Er','Tm','Yb','Lu','Hf','Ta','W ','Re','Os','Ir', 8 'Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn','Fr','Ra', 9 'Ac','Th','Pa','U ','Np','Pu','Am','Cm','Bk','Cf','Es', A 'Fm','Md','No','Lr','Rf','Db','Sg','Bh','Hs','Mt'/ DATA NAP /'p ','d ','t '/, NAA/'-p','-n'/, NON/'??'/ CCC ============================================================ MASS=NINT(A) IZ=NINT(Z) C IF(IZ.EQ.-1) THEN NAME= NAA(1) ELSE IF(IZ.EQ. 0) THEN NAME= NEU ELSE IF(IZ.EQ. 1 .AND. MASS.LE.3) THEN NAME= NAP(MASS) ELSE IF(IZ.GE.1 .AND. IZ.LE.MMAX) THEN NAME= NAT(IZ) ELSE NAME= NON END IF C RETURN END