program Epax c implicit none Real*4 Ap, Zp, At, Zt, A, Z, result, epaxv2 character*1 cyes c c test case: 58Ni + 9Be -> 49Ni: sigma=9.80E-13 barn c Ap=58. Zp=28. At=9. Zt=4. A=49. Z=28. * C write(6,*) write(6,*) '*****************************************' write(6,*) '* *' write(6,*) '* EPAX Version 2 *' write(6,*) '* *' write(6,*) '* An Empirical Parametrization *' write(6,*) '* of Projectile-Fragmentation *' write(6,*) '* Cross Sections *' write(6,*) '* *' write(6,*) '*****************************************' write(6,*) C C ******* INPUT ******************************************************* C 1 call prso('Enter projectile charge', Zp, 0 ) call prso('Enter projectile mass ', Ap, 0 ) call prso('Enter target charge ', Zt, 0 ) call prso('Enter target mass ', At, 1 ) call prso('Enter fragment charge ', Z, 0 ) call prso('Enter fragment mass ', A, 0 ) C if( Zp .le. 0. .or. Ap .le. 0. .or. Zt .le. 0. .or. At .le. 0. # .or. A .le. 0. .or. Z .le. 0. ) then write(6,*) ': Input parameters incorrect!' write(6,*) goto 999 endif C if( Zp .gt. 83. ) then write(6,*) ': EPAX is not designed to calculate', # ' cross sections for fissile nuclei!' write(6,*) endif C if(Z .gt. Zp .or. A .ge. Ap .or. (A-Z) .gt. (Ap-Zp) ) then Write(6,*) ': EPAX is not designed to calculate', # ' transfer cross sections!' write(6,*) goto 999 endif C C result=epaxv2(Ap, Zp, At, Zt, A, Z) C if(result.lt.0) then write(6,*) ': Input parameters incorrect!' goto 999 endif C write(6,*) write(6,*) '*******************************************' write(6,*) 'Fragmentation cross section:' write(6,10) nint(Ap), nint(Zp), nint(At), nint(Zt), nint(A), # nint(Z) 10 format( # ' Projectile (',I3,',',I2,') + Target (',I3,',',I2, # ')',/,13x,' ---->> Fragment (',I3,',',I2,') :') write(6,11) result 11 format(' Sigma = ', 1PE12.6, ' b') write(6,*) '*******************************************' write(6,*) C 999 continue cyes = 'y' call PYES('Another calculation ', cyes) if( cyes .ne. 'n' .and. cyes .ne. 'N') goto 1 end * ***************************************************************************** * Real*4 function epaxv2(Ap,Zp,At,Zt,A,Z) * * Returns cross section (in barn) for fragment (A,Z) in projectile * fragmentation of projectile (Ap,Zp) on target (At,Zt) * * Version: 2.15 18-JUL-99 K.Suemmerer / B. Blank * ***************************************************************************** C implicit none Real*4 Ap, Zp, At, Zt, A, Z, At13, Ap13 Real*4 S, p, f_mod_y, f_mod, zprob, zbeta, delta, dq, u_n, u_p Real*4 expo, fract, z_exp, dfdz, r, yield_a, zbeta_p real*4 p_Delta(4),p_P(2),p_S(2),p_R(2),p_Un(1),p_Up(3) real*4 p_mn(2),p_mp(2),corr_d(2),corr_r(2),corr_y(2) C epaxv2=-999.0 c Ap13 = Ap**0.33333 At13 = At**0.33333 C C************ Constants of EPAX ******************************* C p_S(1) = -2.38 ! scale factor for xsect in barn p_S(2) = 0.27 * p_P(1) = -2.5840E+00 ! slope of mass yield curve p_P(2) = -7.5700E-03 * p_Delta(1) = -1.0870E+00 ! centroid rel. to beta-stability p_Delta(2) = +3.0470E-02 p_Delta(3) = +2.1353E-04 p_Delta(4) = +7.1350E+01 * p_R(1) = +0.885E+00 ! width parameter R p_R(2) = -9.8160E-03 * p_Un(1) = 1.65 ! slope par. n-rich ride of Z distr. * p_Up(1) = 1.7880 ! slope par. p-rich ride of Z distr. p_Up(2) = +4.7210E-03 p_Up(3) = -1.3030E-05 * p_mn(1) = 0.400 ! memory effect n-rich projectiles p_mn(2) = 0.600 * p_mp(1) = -10.25 ! memory effect p-rich projectiles p_mp(2) = +10.1 * corr_d(1) = -25.0 ! correction close to proj.: centroid dzp corr_d(2) = 0.800 corr_r(1) = +20.0 ! correction close to proj.: width R corr_r(2) = 0.820 corr_y(1) = 200.0 ! correction close to proj.: Yield_a corr_y(2) = 0.90 C C ***************************************************************************** C C *** calculate mass yield S = p_S(2) * (At13 + Ap13 + p_S(1)) p = exp(p_p(2)*Ap + p_p(1)) yield_a = p * S * exp(-p * (Ap - A)) C C ** modification close to projectile f_mod_y=1.0 if (A/Ap.gt.corr_y(2)) then f_mod_y=corr_y(1)*(A/Ap-corr_y(2))**2 + 1.0 endif yield_a= yield_a * f_mod_y * print *,'yield_a=', yield_a C C *** calculate maximum of charge dispersion zprob zbeta = A/(1.98+0.0155*A**(2./3.)) zbeta_p = Ap/(1.98+0.0155*Ap**(2./3.)) if(A.gt.p_delta(4)) then delta = p_delta(1) + p_delta(2)*A else delta = p_delta(3)*A*A endif C C ** modification close to projectile f_mod=1.0 if(A/Ap.gt.corr_d(2)) then F_Mod = corr_d(1)*(A/Ap-corr_d(2))**2 + 1.0 endif Delta = Delta*F_Mod zprob = zbeta+delta C C * correction for proton- and neutron-rich projectiles if((Zp-zbeta_p).gt.0) then !------------- p-rich dq = exp(p_mp(1) + A/Ap*p_mp(2)) else !------------- n-rich dq = p_mn(1)*(A/Ap)**2.0 + p_mn(2)*(A/Ap)**4.0 endif zprob = zprob + dq * (Zp-zbeta_p) C C * small corr. since Xe-129 and Pb-208 are not on Z_beta line zprob = zprob + 0.0020*A * print *,'zprob=',zprob C C *** calculate width parameter R r = exp(p_r(1) + p_r(2)*A) C C ** modification close to projectile f_mod=1.0 if (A/Ap.gt.corr_r(2)) then f_mod = corr_r(1)*Ap*(A/Ap-corr_r(2))**4.0+1.0 endif r = r*f_mod C C ** change width according to dev. from beta-stability if ((Zp-zbeta_p).lt.0.0) then r=r*(1.0-0.0833*abs(Zp-zbeta_p)) endif * print *, 'r=',r C C *** calculate slope parameters u_n, u_p u_n = p_un(1) u_p = p_up(1) + p_up(2)*A + p_up(3)*A**2 * print *,'u:',u_n,u_p C C *** calculate charge dispersion if((zprob-Z).gt.0) then C ** neutron-rich expo = -r*(abs(zprob-Z))**u_n fract = exp(expo)*sqrt(r/3.14159) else C ** proton-rich expo = -r*(abs(zprob-Z))**u_p fract = exp(expo)*sqrt(r/3.14159) C * go to exponential slope dfdz = 1.2 + 0.647*(A/2.)**0.3 z_exp = zprob + dfdz * log(10.) / (2.*r) if( Z .gt. z_exp ) then expo = -r*(abs(zprob-z_exp))**u_p fract = exp(expo)*sqrt(r/3.14159) # / (10**(dfdz))**(Z-z_exp) endif endif C * print *,'fract=',fract epaxv2=fract*yield_a END C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C SUBROUTINE PRSO( C$STR, R$REAL, I) C C PROMPTING - SUBROUTINE C INTEGER I CHARACTER*(*) C$STR CHARACTER*(80) C$REAL REAL*4 R$REAL C C ********************************************************************* 100 IF( I .LT. 0 ) THEN WRITE(6,*) ' : I < 0 IN SUBROUTINE PRSO' RETURN ENDIF IF( I .GT. 6 ) I = 6 C IF( I .EQ. 0 ) THEN WRITE(6, 1) C$STR , R$REAL 1 FORMAT(' ',A,' ',F9.0) GOTO 10 ENDIF IF( I .EQ. 1 ) THEN WRITE(6, 2) C$STR , R$REAL 2 FORMAT(' ',A,' ',F9.1) GOTO 10 ENDIF 10 READ( 5, FMT='(A)', ERR=11, END=11 ) C$REAL IF( C$REAL .EQ. ' ' ) RETURN C IF( C$REAL .NE. ' ' ) THEN IF( C$REAL(1:1) .EQ. ' ' ) CALL NOLBL(C$REAL) ENDIF C IF( C$REAL(1:1) .NE. '-' .AND. C$REAL(1:1) .NE. '+' .AND. # C$REAL(1:1) .NE. '1' .AND. C$REAL(1:1) .NE. '2' .AND. # C$REAL(1:1) .NE. '3' .AND. C$REAL(1:1) .NE. '4' .AND. # C$REAL(1:1) .NE. '5' .AND. C$REAL(1:1) .NE. '6' .AND. # C$REAL(1:1) .NE. '7' .AND. C$REAL(1:1) .NE. '8' .AND. # C$REAL(1:1) .NE. '9' .AND. C$REAL(1:1) .NE. '0' .AND. # C$REAL(1:1) .NE. '.') THEN WRITE(6,*) ': Real number required!!!!!!!' WRITE(6,*) GOTO 100 ENDIF C IF( C$REAL .NE. ' ' ) THEN CALL C$CHRE( C$REAL, R$REAL ) RETURN ENDIF C 11 CONTINUE RETURN END C C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C SUBROUTINE PYES( C$STR, C$ANS ) C C PROMPTING - SUBROUTINE FOR QUESTIONS C CHARACTER*(*) C$STR CHARACTER*(80) C$REAL CHARACTER*(1) C$ANS C C ********************************************************************* 1 WRITE(6, 10) C$STR, C$ANS 10 FORMAT(' ',A,'? --> ',A,'?') READ( 5, FMT='(A)',ERR=100, END=100 ) C$REAL C IF( C$REAL .EQ. ' ' ) C$REAL = C$ANS C IF( C$REAL .NE. ' ' ) THEN CALL NOLBL(C$REAL) C$ANS = C$REAL(1:1) IF( C$ANS .EQ. 'J' .OR. C$ANS .EQ. 'Y' .OR. C$ANS .EQ. 'O' # .OR. C$ANS .EQ. 'j' .OR. C$ANS .EQ. 'y' .OR. C$ANS .EQ. 'o') # THEN C$ANS = 'Y' GOTO 100 ENDIF IF( C$ANS .EQ. 'N' .OR. C$ANS .EQ. 'n') THEN GOTO 100 ENDIF WRITE(6,*) ': Unexpected answer! Answer "y" or "n"' GOTO 1 ENDIF C 100 CONTINUE RETURN END C C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C SUBROUTINE NOLBL(CSTR) C C SUPPRESSES LEADING BLANKS OF STRING CSTR C CHARACTER*(*) CSTR CHARACTER*1 CWORK INTEGER*4 ILEN, I C C ********************************************************************* ILEN = LEN(CSTR) DO 1 I = 1,ILEN CWORK = CSTR(I:I) IF( CWORK .NE. ' ' ) THEN CSTR = CSTR(I:) GOTO 2 ENDIF 1 CONTINUE 2 CONTINUE C END C C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C SUBROUTINE C$RECH( R$REAL, C$STR ) C C CONVERSION: REAL VARIABLE -> CHARACTER - VARIABLE C CHARACTER*8 C$DUMMY CHARACTER*(*) C$STR REAL*4 R$REAL C C ********************************************************************* C$DUMMY = C$STR WRITE( C$DUMMY, FMT='(F8.4)' ) R$REAL CALL NOLBL(C$DUMMY) C$STR = C$DUMMY C RETURN END C C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C SUBROUTINE C$CHRE( C$STR, R$REAL ) C C CONVERSION: CHARACTER VARIABLE -> REAL - VARIABLE C CHARACTER*(*) C$STR REAL*4 R$REAL C C ********************************************************************* READ( C$STR , * ) R$REAL C END C