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