      program accba
c     Adiabatic Coupled-Channels Born Approx. valid only for (d,2He) Reaction
      implicit real*8(a-h,o-z)
      include 'parameter.f'
      integer*4  sa,spa,sb,sp,t,jtr(njtrmx)
      real*8     IA,IB,E12(nE12mx),wE12(nE12mx),mp,mu
      parameter (sa=1, spa=1, sb=0, sp=1)
      parameter (mp=2., Z1=1., nua=0, nub=-1, t=1, Ast=-1.) ! 0.5) ???
      parameter (mu=931.494*(1.007275+1.008665)/2.)
      parameter (hc=197.32705, e2=hc/137.036)
      real*8     ff(0:nRmax,2,2,njtrmx),wrs(0:nRmax)
      complex*16 H(0:nRmax,0:L1max)
      real*8     chi0(0:nRmax,nE12mx),sigma1(0:Lmax0),sigma2(0:Lmax0)
      complex*16 chi1(0:nRmax,-1:1,0:Lmax0)
      complex*16 chi2(0:nRmax,0:Lmax0),pUopt0
      complex*16 Y(0:(Lmax0+1)*(Lmax0+2))
      complex*16 beta(0:Lmax0,-sa:sa,0:jtrmax,njtrmx,nE12mx)
      complex*16 Samp(-sa:sa),S
      parameter (nobs=5)  ! cs,Ay,Ayy,Axx,Axz
      real*8     obs(nobs),obs0(nobs)
      logical*4  DWBA,pntprj,tridif,decay
      character  title*80
      pi =atan(1.D0)*4.D0
      pi4=pi*4.D0
      call inirac	! initialize faclog for cleb, rac, u9
c     basic configuration
      read(5,'(A)') title
      read(5,*) DWBA,pntprj,tridif,decay
      tridif=(.not.pntprj.and.tridif)
      read(5,*) thmin,thmax,dth
      read(5,*) A,ZA,IA,IB,Elab,Qval,Ex
      B =A
      ZB=ZA-nua+nub
      Z2=Z1+nua-nub
      read(5,*) iptr,njtr,(jtr(i),i=1,njtr)
      iptr =sign(1,iptr)  ! target parity transfer
      if (njtr.gt.njtrmx) stop 'too large njtr'
      mxjtr=0
      do ijtr=1,njtr
         if (jtr(ijtr).gt.jtrmax) stop 'too large jtr'
         mxjtr=max(mxjtr,jtr(ijtr))
      enddo
      read(5,*) nRL,dRL,Lmax
      if (nRL.gt.nRmax) stop 'too large nRL'
      if (mod(nRL,2).ne.0) stop 'nRL must be even'
      if (Lmax+mxjtr+sa.gt.Lmax0) stop 'too large Lmax'
      read(5,*) nrs,drs,mrs,lamax,lbmax
      if (nrs*mrs.gt.nRmax) stop 'too large nrs*mrs'
      if (lamax.ne.0) lamax=2 ! include deuteron D-state
      if (lbmax.ne.0) lbmax=2 ! allow coupling up to singlet-D
      read(5,*) nE12,E12min,E12max
      if (nE12.gt.nE12mx) stop 'too large nE12'
      if (pntprj) then
         nrs=1
         drs=0.
         nE12=1
         E12min=0.
         E12max=0.
         E12(1)=0.
         wE12(1)=1.
         lamax=0
         lbmax=0
      else
         if (mod(nrs,2).ne.0) stop 'nrs must be even'
         call gauleg(E12min,E12max,E12,wE12,nE12)
         call rsck(nE12,E12,nrs,drs,mrs,chi0) ! prepare 2He wave-function
         wrs( 0 )=drs/3.D0  ! not used as chi0(0)=0
         wrs(nrs)=drs/3.D0
         do irs=1,nrs-1     ! weight for Simpson method
            wrs(irs)=drs*(mod(irs,2)+1)*2.D0/3.D0
         enddo
      endif
      ltitle=len(title)
      do while(ltitle.gt.1.and.title(ltitle:ltitle).eq.' ')
         ltitle=ltitle-1
      enddo
      write(6,'(3A)') '"',title(:ltitle),'"'
      write(6,'(/A,4(2X,A,L2))')
     &' calc. options:','DWBA=',DWBA,'point-proj.=',pntprj,
     &'triple-diff.c.s.=',tridif,'particle-decay=',decay
      write(6,'(/A,SP,I3,A,S,10(I3,:1H,))')
     &' parity-transfer=',iptr,'  j-transfer(s)=',(jtr(i),i=1,njtr)
      write(6,'(/A/)') ' integration & summation parameters.'
      write(6,'(A,I4,A,F5.2,9X,A,I4)')
     &'  nRL=',nRL,'  dRL=',dRL,'  Lfmax=',Lmax
      write(6,'(A,I4,A,F5.2,A,I3,2(A,I2))')
     &'  nrs=',nrs,'  drs=',drs,'  mrs=',mrs,
     &'  lamax=',lamax,'  lbmax=',lbmax
      write(6,'(A,I3,2(A,F6.2))')
     &'  nE12=',nE12,'  E12min=',E12min,'  E12max=',E12max
      write(6,'(/A/)') ' angular dist. parameters.'
      write(6,'(3(A,F6.1))')
     &'  thmin=',thmin,'  thmax=',thmax,'  dth=',dth
c     prepare incident-channel distorted wave-func.
      Ecm1=Elab*A/(mp+A)
      amu =mp*A/(mp+A)
      fk1 =sqrt(2.D0*mu*amu*Ecm1)/hc
      eta1=Z1*ZA*e2*mu*amu/(hc**2*fk1)
      write(6,'(/A/)') ' incident channel parameters.'
      write(6,'(A,F7.2,A,F5.1,5(2X,A,F4.1))') '  Elab=',Elab,
     &'  A=',A,'ZA=',ZA,'IA=',IA,'a=',mp,'Za=',Z1,'sa=',sa*1.
      write(6,'(A,F8.3,2(2X,A,F8.5))')'  Ecm=',Ecm1,'k=',fk1,'eta=',eta1
      call cphase(Lmax,eta1,sigma1)
      call distwf(A,Z1*ZA,amu,fk1,eta1,nRL,dRL,Lmax+mxjtr+sa,sa,chi1)!read pot.
c     prepare exit-channel scattering wave-func.
      Ecm2=Ecm1+Qval-Ex
      amu =mp*B/(mp+B)
      fk2 =sqrt(2.D0*mu*amu*Ecm2)/hc
      eta2=Z2*ZB*e2*mu*amu/(hc**2*fk2)
      write(6,'(/A/)') ' exit channel parameters.'
      write(6,'(A,F7.2,A,F5.1,5(2X,A,F4.1),2X,A,F6.2)') '  Qval=',Qval,
     &'  B=',B,'ZB=',ZB,'IB=',IB,'b=',mp,'Zb=',Z2,'sb=',sb*1.,'Ex=',Ex
      write(6,'(A,F8.3,2(2X,A,F8.5))')'  Ecm=',Ecm2,'k=',fk2,'eta=',eta2
      call cphase(Lmax,eta2,sigma2)
      if (DWBA) then
         call distwf(B,Z2*ZB,amu,fk2,eta2,nRL,dRL,Lmax,sb,chi2) ! read pot.
      else
         dummy=dble(pUopt0(B,ZB,mp,Z2/2.D0,1)) ! read-in p-B optical pot.
      endif
c     prepare "ordinary" form-factor
      call ffint0                    ! read-in effective interaction param.
      call shell0(A,ZA)              ! read-in single-particle potential param.
      call formf0(A,nRL,dRL,0,*90000)! read-in & prepare single-part. orbits.
      do ijtr=1,njtr
         j=jtr(ijtr)
         iptr2=(1-iptr*(-1)**j)/2 ! 0 for natural-par., 1 for unnatural-par.
         ilp=0
         call formf1(j,*90000)       ! read-in p-h configurations
         write(6,'(/A/)') ' "ordinary" form-factor calculated for ...'
         do lp=abs(j-sp*iptr2),j+sp*iptr2,2
            ilp=ilp+1
            call formf(lp,sp,j,t,IB,nua,nub,Ast,
     &                 nRL,dRL,ff(0,1,ilp,ijtr))
            write(6,'(4(A,I1))') '  j=',j,', l=',lp,', s=',sp,', t=',t
         enddo
      enddo
c     initialize scattering amplitude
      do Lf=0,Lmax
         do ma=-sa,sa
            do m=0,mxjtr
               do ijtr=1,njtr
                  do iE12=1,nE12
                     beta(Lf,ma,m,ijtr,iE12)=(0.D0,0.D0)
                  enddo
               enddo
            enddo
         enddo
      enddo
c     calculate scattering amplitude for each rs
      do irs=1,nrs
         rs=irs*drs
         write(0,'(A,F6.2,2A)') ' rs=',rs,char(27),'M'
         if (.not.DWBA) call adiawf(rs,lbmax,nRL,dRL,Lmax,fk2,eta2,chi2)
         do ijtr=1,njtr
            j=jtr(ijtr)
            iptr2=(1-iptr*(-1)**j)/2 ! 0 for natural-par., 1 for unnatural-par.
            ilp=0
            do lp=abs(j-sp*iptr2),j+sp*iptr2,2
               ilp=ilp+1
               do la=0,lamax,2
                  call ffexpd(rs,lp,la,nRL,dRL,ff(0,1,ilp,ijtr),H)
                  do L1=max(abs(lp-la),abs(j-sa)),
     &                  min(    lp+la ,    j+sa ),2
                     fac1=sqrt(dble((2*L1+1)*(2*lp+1)*(2*sa+1)))*
     &                    rac(2*L1,2*la,2*j,2*spa,2*lp,2*sa)
                     do Lf=0,Lmax
                        do Li=abs(Lf-L1),Lf+L1,2
                           fac2=sqrt(dble((2*Lf+1)*(2*Li+1)))*
     &                          clebz(2*Lf,2*L1,2*Li)*
     &                          (-1)**((Li-Lf-L1)/2)
                           do Ji=max(abs(Li-sa),abs(Lf-j)),
     &                           min(    Li+sa ,    Lf+j )
                              fac3=rac(2*Lf,2*L1,2*Ji,2*sa,2*Li,2*j)
                              fac=fac1*fac2*fac3
                              if (fac.eq.0.D0) write(6,*)
     &                     'j,lp,la,L1,Lf,Li,Ji=',j,lp,la,L1,Lf,Li,Ji
                              JS=Ji-Li
c                             integrate over RL with Simpson method
                              S=chi1(nRL,JS,Li)*H(nRL,L1)*chi2(nRL,Lf)!chi(0)=0
                              do iR=1,nRL-1
                                 S=S+chi1(iR,JS,Li)*H(iR,L1)*chi2(iR,Lf)
     &                              *(mod(iR,2)+1)*2
                              enddo
                              S=S/3.D0*dRL*fac1*fac2*fac3
     &                         *exp(dcmplx(0.D0,sigma1(Li)+sigma2(Lf)))
                              if (pntprj) then
                                 S=S*sqrt(pi4)/Ast  ! w/o 4*pi*Y00(k12)
                              else
                                 S=S*pi4*rscd(la,rs)*wrs(irs)
                              endif
                              do ma=-min(sa,Ji),min(sa,Ji)
                                 fac4=cleb(2*Li,0,2*sa,2*ma,2*Ji,2*ma)
                                 do m=0,j
                                    fac5=cleb(2*Lf,2*(ma-m),2*j,2*m,
     &                                        2*Ji,2*ma)
                                    if (pntprj) then
                                       beta(Lf,ma,m,ijtr,1)=
     &                                 beta(Lf,ma,m,ijtr,1)+fac4*fac5*S
                                    else
                                       do iE12=1,nE12
                                       beta(Lf,ma,m,ijtr,iE12)=
     &                                 beta(Lf,ma,m,ijtr,iE12)+
     &                                    fac4*fac5*S*chi0(irs,iE12)
                                       enddo
                                    endif
                                 enddo
                              enddo
                           enddo
                        enddo
                     enddo
                  enddo
               enddo
            enddo
         enddo
      enddo
c     readin particle-decay parameter
70000 if (decay) call decay0(IB,*70000,*80000)
c     sum-up scattering amplitudes at each angle, using reflection-invariance:
c     beta(Lf,ma,m)*Y(Lf,ma-m)=iptr*(-)**(j-m+sa-ma)*beta(Lf,-ma,-m)*Y(Lf,m-ma)
      write(6,'(/2A/)') '# ',title(:ltitle)
      write(6,'(1H#,A6,A10,4A9)') 'angle','xsec','Ay','Ayy','Axx','Axz'
      fac=10./Ecm1/Ecm2*fk2/fk1*(2*IB+1)/(2*IA+1)/(2*sa+1) ! 10 for fm^2 => mb
      if (pntprj) then
         fac=fac/pi4**2
      else
         fac2=fac*4./pi4**5 ! for triple-diff. cross section
         fac=fac2*pi4        ! 4*pi for int. over p-p rel. angle
      endif
      do thdeg=thmin,thmax+0.1*dth,dth
         if (tridif) write(1,'(A,F7.2)') '# th=',thdeg
         th=thdeg/180.D0*pi
         call spherical(th,0.D0,Lmax,Y) ! phi=0 i.e. s//y-axis in Madison conv.
         do i=1,nobs
            obs(i)=0.D0
         enddo
         if (decay) call decay1(njtr,jtr) ! clear storage
         do iE12=1,nE12
            do i=1,nobs
               obs0(i)=0.D0
            enddo
            do ijtr=1,njtr
               j=jtr(ijtr)
               do m=-j,j
                  do ma=-sa,sa  ! j,m,ma are incoherent sum
                     Samp(ma)=(0.D0,0.D0)
                     do Lf=abs(ma-m),Lmax  ! coherent sum
                        if (m.ge.0) then
                           Samp(ma)=Samp(ma)+
     &                       beta(Lf, ma, m,ijtr,iE12)*Y(Lf*(Lf+1)+ma-m)
                        else
                           Samp(ma)=Samp(ma)+iptr*(-1)**(j-m+sa-ma)*
     &                       beta(Lf,-ma,-m,ijtr,iE12)*Y(Lf*(Lf+1)-ma+m)
                        endif
                     enddo
                  enddo
                  call amp2cs(Samp,obs0)
                  if (decay) call decay2(Samp,m,ijtr) ! store Samp
               enddo
            enddo
            if (decay) call decay3(njtr,jtr,IA,IB,wE12(iE12)) ! calc. angcor
            fk=sqrt(mu*E12(iE12))/hc
            do i=1,nobs
               if (pntprj) then
                  obs(i)=obs(i)+obs0(i)*wE12(iE12) ! weighted mean
               else
                  obs(i)=obs(i)+obs0(i)*wE12(iE12)*fk/E12(iE12)
               endif
            enddo
            if (tridif) write(1,'(0pF8.5,1pE12.4)')
     &           E12(iE12),obs0(1)*fac2*fk/E12(iE12)
         enddo
         if (tridif) write(1,'(A)') 'join'
         do i=2,nobs
            obs(i)=obs(i)/obs(1)
         enddo
         obs(1)=obs(1)*fac
         write(6,'(0pF7.2,1pE12.3,0p4F9.4)') thdeg,(obs(i),i=1,nobs)
         if (decay) call decay4(E12max-E12min) ! output angular correlation
      enddo
      if (decay) goto 70000
80000 stop
90000 stop 1024
      end

      include 'decay.f'
      include 'denmat.f'
