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