program angelmc
c
c MB Chadwick modified angel92 to include photuclear reactions,
c Oct 95, according to the paper to be published in J. Nucl.
c Sci Eng. (Japan), Nov 1995. Changes labeled with "cmbc"
c ANGEL92
c continuum angular distributions
c from empirical systematics based on exponentials of cos theta
c see Phys. Rev. C 37 (1988) 2350
c
c Written by C. Kalbach
c April 1987
c Revised February 1992 to give a smooth transitions
c Single formulae for E1 and E3.
c see Nucl. Sci. Eng. 115 (1993) 43
c
c Third term in 'a' needed for incident N and d but not alpha.
c Tentatively used for incident t and He-3
c Transition energies tentatively set at
c Et1 = 130 MeV (a=n,p,d,t,He3) or 260 MeV (a=alpha)
c Et3 = 35 MeV (a=n,p,d,t,He3)
c (only nucleon values of Et1; N and d values of Et3 known)
c
c *The angel of the Lord encampeth round about them that fear him,
c *and delivereth them. (Psalm 34.7)
c
dimension sigma(19),jang(19),slope(50),xnorm(50),eps(50),fmsd(50)
1,total(50),xcos(19)
character*64 title
character*1 newpg
real*8 arg
common /angels/ xs(5,5)
open(1,file='angel.dat')
open(2,file='angel.out',status='unknown')
newpg=char(12)
do 100 i=1,10
100 sigma(i)=0
do 102 i=1,19
theta=i*10-10
jang(i)=theta
theta=theta*3.1413/180.
102 xcos(i)=cos(theta)
c
c store array of mass excesses
c i=z+1, j=n+1
c
do 104 i=1,5
do 104 j=1,5
104 xs(i,j)=0.
xs(2,2)=2.22
xs(2,3)=8.48
xs(3,2)=7.72
xs(3,3)=28.30
c
c reaction i.d. input
c
irea=1
iwri=2
110 read(irea,112)title
112 format(1a64)
read(irea,116)jcom,jzcom,jpin,jnin,jpout,jnout
116 format(bn,6i5)
if(jcom.le.0)go to 40
read(irea,120)e
120 format(1f6.1)
read(irea,116)neps
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cmbc MB Chadwick, added coding Oct 95, for photonuclear reactions
jflagph=0
if(jpin.eq.0.and.jnin.eq.0)then
c the following is a modification for photonuclear reactions
jflagph=1
jcom=jcom+1
jnin=1
endif
cmbc I follow the prescription in the paper "Photonuclear angular
cmbc distribution systematics in the quasideuteron regime", MB Chadwick,
cmbc PG Young, and S Chiba, to be published in the Nov 1995 issue of
cmbc the Journal of Nucl. Sci and Tech. (Japan). The basic idea is
cmbc to first calculate the Kalbach MSD a-parameter assuming a neutron
cmbc projectile instead of a photon, and the modify a to account for
cmbc (a) reduced forwward peaking since a photon carries less momentum;
cmbc (b) less diffraction/refraction effects since the photon is undistorted.
cmbc end of mbc code (more later)
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c calculate energy independent quantities
c and print headings
c
acom=jcom
azcom=jzcom
jin=jpin+jnin
xin=jin
jout=jpout+jnout
esys=e*(acom-xin)/acom+bind(acom,azcom,jpin,jnin)
bin=bind(acom,azcom,jpout,jnout)
write(iwri,152)newpg,title
152 format(1a1,1a64)
write(iwri,154)esys,bin
154 format(17h sys. excit. en.=,1f6.2,18h; sys. bind. en.=,1f5.2)
er=1.
if (jin.gt.3) er=2.
c er = Et1/130.
esysr = esys / er
e1 = esys
if (esysr.gt.190.) then
e1 = 130.*er
else if (esysr.gt.80.) then
x = (esysr-130.)/12.
x = 1. + exp(x)
e1 = esys / x
x = (136.-esysr)/12.
x = 1. + exp(x)
e1 = e1 + 130.*er/x
end if
if (jin.gt.3) go to 20
e3 = esys
if (esysr.gt.51.) then
e3 = 35.*er
else if (esysr.gt.21.) then
x = (esysr-35.)/3.2
x = 1. + exp(x)
e3 = esys / x
x = (36.6-esysr)/3.2
x = 1. + exp(x)
e3 = e3 + 35.*er/x
end if
xmb=1
if(jout.eq.4)xmb=2.
if(jpout.eq.0)xmb=0.5
c
c energy dependent input
c calculate and print angular distributions
c
20 write(iwri,18)
write(iwri,18)
write(iwri,16)(jang(j),j=1,10)
16 format(5h eps,1i7,9i9)
write(iwri,18)
18 format(1h )
do 28 ne=1,neps
read(irea,14)eps(ne),total(ne),fmsd(ne)
14 format(3f10.4)
if(fmsd(ne).le.0.)fmsd(ne)=1.
epscm=eps(ne)+bin
y = epscm*e1 / (esys*130.)
a = 5.2*y + 4.*y*y*y
if (jin.gt.3) go to 22
y = epscm*e3 / (esys*35.)
a = a + 1.9*y*y*y*y*xmb
22 xnorm(ne) = a*total(ne) / (12.5664*sinh(a))
slope(ne)=a
do 24 i=1,10
arg=a*xcos(i)
sig=fmsd(ne)*dsinh(arg)+dcosh(arg)
cmbc%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cmbc mbc added code for photonuclear reactions. Modify
cmbc a which was calculated for a neutron, to be valid for a
cmbc photon projectile:
if(jflagph.eq.1)then
facmom=sqrt(e/(2.*939))
facrefr=9.3/sqrt(eps(ne))
if(facrefr.lt.1.)facrefr=1.
if(facrefr.gt.4.)facrefr=4.
aphnuc=a*facmom*facrefr
gth=((2.*aphnuc)/(exp(aphnuc)-exp(-aphnuc)))
gth=gth*(1./(12.5664))*exp(aphnuc*xcos(i))
cmbc Now put in MSC as isotropic, giving:
sigma(i)=((1.-fmsd(ne))/(12.5664))+(fmsd(ne)*gth)
sigma(i)=sigma(i)*total(ne)
go to 24
endif
cmbc%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sigma(i)=sig*xnorm(ne)
24 continue
write(iwri,26)eps(ne),(sigma(i),i=1,10)
26 format(1f6.2,10(1pe9.2))
28 continue
write(iwri,18)
write(iwri,18)
write(iwri,30)(jang(j),j=11,19)
30 format(5h eps,1i7,8i9,17h total fmsd)
write(iwri,18)
do 38 ne=1,neps
do 34 j=11,19
arg=slope(ne)*xcos(j)
sig=fmsd(ne)*dsinh(arg)+dcosh(arg)
cmbc%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cmbc mbc added code for photonuclear reactions. Modify
cmbc a which was calculated for a neutron, to be valid for a
cmbc photon projectile:
if(jflagph.eq.1)then
facmom=sqrt(e/(2.*939))
facrefr=9.3/sqrt(eps(ne))
if(facrefr.lt.1.)facrefr=1.
if(facrefr.gt.4.)facrefr=4.
aphnuc=a*facmom*facrefr
gth=((2.*aphnuc)/(exp(aphnuc)-exp(-aphnuc)))
gth=gth*(1./(12.5664))*exp(aphnuc*xcos(j))
cmbc Now put in MSC as isotropic, giving:
sigma(j)=((1.-fmsd(ne))/(12.5664))+(fmsd(ne)*gth)
sigma(j)=sigma(j)*total(ne)
go to 34
endif
cmbc%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sigma(j)=sig*xnorm(ne)
34 continue
fm=0.1*fmsd(ne)
write(iwri,36)eps(ne),(sigma(i),i=11,19),total(ne),fm
36 format(1f6.2,10(1pe9.2),1f6.3)
38 continue
go to 110
40 stop
end
c
c
function bind(acom,azcom,jpout,jnout)
c
c calculate binding energies from semi-empirical masses
c with pairing (and shell) effects not included
c
common /angels/ xs(5,5)
ancom=acom-azcom
xout=jpout+jnout
xpout=jpout
ares=acom-xout
azres=azcom-xpout
anres=ares-azres
acom3=acom**0.33333
ares3=ares**0.33333
bin=15.68*xout
bin=bin-18.56*(acom3*acom3-ares3*ares3)
x=(ancom-azcom)*(ancom-azcom)/acom
y=(anres-azres)*(anres-azres)/ares
bin=bin-28.07*(x-y)
bin=bin+33.228*(x/acom3-y/ares3)
bin=bin-0.717*(azcom*azcom/acom3-azres*azres/ares3)
bin=bin+1.211*(azcom*azcom/acom-azres*azres/ares)
bind=bin-xs(jpout+1,jnout+1)
return
end