program linmom
c By M.B. Chadwick, LANL Los Alamos
c Implements Chadwick-Oblozinsky continuum angular distribution theory
c See Phys. Rev. C50, 2490 (1994).
c Sample input corresponds to Figure 1 of this reference.
dimension rat(20),xs(20),g(0:20)
open(unit=8,file='output',status='unknown')
open(unit=7,file='input',status='unknown')
pi=3.1415927
c ***********************************************************************
c Input paramters as follows: (ALL FREE FORMAT)
c ein, eout = incident and emission energies
c bindin and bindout = sparation energies for projectile and ejectile
c from the composite nucleus
c efermi = Fermi energy (35 MeV)
c dsde = angle-integrated preequilibrium cross section at this emiss. energy
c amass = target A value
c nstages = number of preequilibrium stages to be included
c rat(i) = the relative contribution of the preequilibrium stage i
c to the total preequilibrium spectrum, for emission energy eout.
c ***********************************************************************
c READ INPUT PARAMETERS:-
totrat=0.
read(7,*)ein,eout,bindin,bindout,efermi,dsde,amass
read(7,*)nstages
do 8 i=1,nstages
read(7,*)rat(i)
8 totrat=totrat+rat(i)
c determine cross sections for the various preequilibrium stages
c such that they sum to the inputted value dsde:
do 9 i=1,nstages
xs(i)=dsde*rat(i)/totrat
9 continue
c zero arrays
do 20 n=1,20
20 g(n)=0.
c zeta factor to approximate effects of defraction/refraction at low
c energies (only comes into play below 80 MeV)
zeta=9.3/sqrt(eout)
if(zeta.lt.1.)zeta=1.
do 120 nstep=1,nstages
step=float(nstep)
anex=2.*step
write(8,*)'Angular dist. for preequilibrium stage N=',nstep,':'
write(8,*)'Ang. ddxs'
write(8,*)'(deg) (mb/MeV sr)'
write(8,*)'----- -----------'
c establish average energy of the excitons in the res. nucleus:-
np=nstep
nh=nstep
energy=ein+bindin-eout-bindout
glev=amass/13.
c glev eventually cancels
call energyav(glev,np,nh,energy,efermi,eav)
totnorm=0.
do 100 nang=0,18
th=float(nang)
factor=sqrt((ein+bindin+efermi)*(eout+bindout+efermi))/
& (zeta*anex*eav/3.)
ctheta=cos(th*pi/18.)
gdist=exp(factor*ctheta)
fnorm=factor/(2.*pi*( exp(factor) - exp(-factor) ))
gd=gdist*fnorm
gdist=gd*xs(nstep)
write(8,91)th*10.,gdist
totnorm=totnorm+2.*pi*sin(th*pi/18.)*(pi/18.)*gd
91 format(f5.1,3x,1p1e10.3)
g(nang)=g(nang)+gdist
100 continue
120 continue
c Now normalize g(nangle) so int over 4pi=1.
tot=0.
do 200 nth=0,18
th=float(nth)
tot=tot+(2.*pi*sin(th*pi/18.)*(pi/18.)*g(nth))
200 continue
write(8,*)' '
write(8,*)'Angular distribution for sum of all preeq. stages:'
write(8,*)'Ang. ddxs'
write(8,*)'(deg) (mb/MeV sr)'
write(8,*)'----- -----------'
c write out double diff x/s for all preequilibrium stages summed:
do 300 nth=0,18
300 write(8,99)float(nth*10),g(nth)
99 format(f5.1,3x,1p1e10.3)
stop
end
subroutine energyav(glev,np,nh,energy,efermi,eav)
c Calculates average exciton energy relative to bottom of nuclear well,
c within an equidistant-spacing model, see Eq. 5 in Phys Rev C50,2490 (1994).
ap=float(np)
ah=float(nh)
an=ap+ah
t1=2.*ap*(ap+1.)/(an*glev)
call rhorestr(glev,np+1,nh,energy,efermi,om1)
call rhorestr(glev,np,nh,energy,efermi,om2)
t1=t1*om1/om2
t2=energy/an
t3=efermi
eav=t1-t2+t3
return
end
subroutine rhorestr(glev,nnp,nnh,energy,efermi,omega)
c determined restricted p-h state density using Betak-Dobes formula
dimension fact(0:30)
c generate factorials:-
fact(0)=1.
do 30 n=1,20
30 fact(n)=fact(n-1)*float(n)
bp=nnp
bh=nnh
aph=(bp*bp + bh*bh +bp - 3.*bh)/(4.*glev)
tot1=0.
j=-1
10 j=j+1
t1=(-1.)**float(j)
t2=fact(j)/ ( fact(nnh-j)*fact(nnh) )
t3=(energy - float(j)*efermi - aph)
if(t3.le.0.)t3=0.
t3=t3**float(nnp+nnh-1)
tot1=tot1 + (t1*t2*t3)
if(j.le.h)go to 10
con=(glev**float(nnp+nnh) )/(fact(nnp)*fact(nnh)*fact(nnp+nnh-1))
omega=con*tot1
return
end