c======================================================================
      PROGRAM ANGDIST
c======================================================================
c calculates angular distribution in nucleus-nucleus collisions, having
c as input the system and the nuclear S-matrix at each partial wave.
c C.A. Bertulani, 1990
c----------------------------------------------------------------------
c
      parameter(lmax=200)
      implicit real*8(a-h,o-z)
      character*12 file_in,file_out,file_dat
      complex*16 sn(0:lmax),fcou,fn,ci
      dimension cph(0:lmax),pl(0:lmax)
c
      open(11,file='angdist.in',status='old')
      read(11,*)file_in,file_dat,file_out
      open(12,file=file_in,status='old')
      open(13,file=file_dat,status='unknown')
      open(14,file=file_out,status='unknown')
c 
      read(11,*)zp,ap,zt,at,ecm,tmin,tmax,dt
      nt=int((tmax-tmin)/dt)+1
      pi=4.d0*atan(1.d0)
      ci=cmplx(0.d0,1.d0)        
      am=931.5*ap*at/(ap+at)
      am2=2.*am
      h2=197.33**2
      aux1=am2*ecm/h2
      xk=dsqrt(aux1)
      eta=1.44*zp*zt*am/xk/h2
c
      write(13,100)zp,ap,zt,at,ecm,lmax
      call coul(lmax,eta,cph)
c
      do l=0,lmax
         read(12,*)ll,sn(l)
         cps=cph(l)/pi
         snr=real(sn(l))
         sni=dimag(sn(l))
         psn=atan2(sni,snr)/pi/2.
         write(13,150)l,abs(sn(l)),psn,cps
      end do
c
      write(13,200)
      do it=1,nt
         tdeg=tmin+dble(it-1)*dt
         t=tdeg*pi/180.
         cost=cos(t)
         call legpol(lmax,cost,pl)         
         fn=cmplx(0.d0,0.d0)
         do l=0,lmax
            fn=fn+dble(2*l+1)*exp(2.*ci*cph(l))*
     &         (sn(l)-cmplx(1.d0,0.d0))*pl(l)/(2.*ci*xk)
         end do
         s2=sin(t/2.)**2
         au=-eta*log(s2)+pi+2.*cph(0)
         fcou=eta*exp(ci*au)/(2.*xk*s2)
         sigma=10.*abs(fn+fcou)**2
c multiplication by 10 transforms the cross section in mb
c
         a=eta/xk
         sigruth=10.*(a/(2.*s2))**2
         ratio=sigma/sigruth
         write(14,250)tdeg,ratio,sigma
         write(13,250)tdeg,ratio,sigma
      end do
c
100   format(//, ' System: zp = ',f4.0,'  ap = ',f4.0,' zt = ',f4.0,
     &       '  at = ',f4.0,/,
     &       ' Ecm = ',g10.2,4x,'lmax = ',i3,//,
     &       '   l  ',4x,'  | S_l |',4x,' N p-shift/pi ',4x,
     &       ' C p-shift/pi',/)
200   format(//,' Theta(deg)',5x,'sig/sruth ',2x,'sigma (mb) ',/)
150   format(1x,i3,4x,g12.4,4x,g12.4,4x,g12.4,4x,g12.4)
250   format(3(1x,g12.4))
      stop
      end
c
c
c=====================================================================
      subroutine legpol(maxl,cosbet,pl)
c=====================================================================
c
c     Calculates Legenrde Polynomials  PL(L) = P_L (COSBET)
c
c---------------------------------------------------------------------
c
      implicit real*8 (a-h,o-z) 
      dimension pl(0:maxl)
c
      pl(0) = 1.d0
      pl(1) = cosbet
      do 1 l = 2,maxl
    1 pl(l) = ((2*l-1)*cosbet*pl(l-1)-(l-1)*pl(l-2))/l
c
      return
c-end-legpol
      end
c
c======================================================================
      subroutine coul(lmax,eta,cph)
c----------------------------------------------------------------------
c calculates Coulomb phase-shifts
c----------------------------------------------------------------------
      implicit real*8(a-h,o-z)
      dimension cph(0:lmax)
c
      sto=16.d0+eta*eta
      cph(0)=-eta+(eta/2.d0)*dlog(sto)+3.5d0*datan(eta/4.d0)-(
     1 datan(eta)+datan(eta/2.d0)+datan(eta/3.d0))-(eta/(12.d0*
     2 sto))*(1.d0+(1.d0/30.d0)*(eta**2-48.d0)/sto**2+(1.d0/105.d0)
     3 *(eta**4-160.d0*eta**2+1280.d0)/sto**4)
c
      do 1 ii=1,lmax
         fi=dble(ii)
         cph(ii)=cph(ii-1)+datan(eta/fi)
1     continue
      return
      end