c     program difls version 2     D. C. Griffin   Last modified: 08/08/08
c
c     program to calculate the angular distribution of the differential
c     excitation cross section and calculate the total excitation cross
c     section for electron-atom and electron-ion scattering in the ls
c     coupling scheme.
c
      implicit real*8(a-h,o-z)
c
c     parameter dimensions
c
c     ntm--maximum number of target terms
c     nfm--maximum number of final target terms to which differential
c           cross sections are to be calculated
c     ncm--maximum number of channels in k-matrices
c     nctm--maximum number of channels for a given term
c     nlm--maximum number of angular momenta for continuum electron
c     njm--maximum number of momentum transfer values
c     nwm--maximum number of partial waves actually used in calculation
c     nftm--maximum number of different factorials used in coupling
c           coefficient calculations - 800 should be enough
c     ntem-maximum number of different angles
c
      parameter(ntm=100,nfm=10,ncm=200,nctm=10,nlm=71,njm=2*nlm,
     1 nlbm=2*nlm,nwm=102,nftm=800,ntem=181)
      character name(ntm)*10,pltunt(5)*12
      dimension title(10),ent(ntm),lt(ntm),lspt(ntm),lpt(ntm)
      common/blok1/einc,theta(ntem),sigbri(nlm),sigbrf(nlm,nfm),
     1 pcf(61,nfm),thpcf(61),thmn,thmx,nthta,iprint,icharg,llmx,ntf,
     2 iti,itf(ntm),nthmx,lti,nspin,ispin(2),ltf(nfm),limn,limx,
     3 lfmx(nlm,nfm),lfmn(nlm,nfm),lspin,ielas,iphr,nipc,nfpc
      common/blok2/gamma(nftm)
      common/blok3/tr(nwm,ncm,ncm),ti(nwm,ncm,ncm),ltot(nwm),lstot(nwm),
     1 ncci(nwm),nccf(nwm,nfm),icci(nwm,nctm),iccf(nwm,nctm,nftm),
     2 lcs(nwm,ncm),iwu,imrd
      common/blok4/fmr(nlm,nfm,nlm,njm,2),fmi(nlm,nfm,nlm,njm,2),
     1 qr(nlm),qi(nlm),alm(nfm,nlbm),dsg(ntem,nfm),sig(nfm),sigic(nfm)
      common/blok5/pl(nlbm),imbmx
      namelist/difls/iphr,iprint,icharg,llmx,nt,iti,ntf,imrd,
     1 nthta,thmn,thmx,ein,iunit
c
      data zero,onehndrd,one80/0.0d0,1.00d2,1.80d2/
      data convev,convsg,convsg1,convsg2,convsg3,convsg4/13.6058d0,
     1 2.80028d-17,1.0d-19,1.0d-18,1.0d-17,1.0d-16/
      data (pltunt(i),i=1,5) /'   a0^2/sr  ','10-19 cm2/sr',
     1 '10-18 cm2/sr','10-17 cm2/sr','10-16 cm2/sr'/
      open(5,file='difin')
      open(6,file='difout')
      open(7,file='difplt.dat')
c
c     input from difin
c
c     title(10)--title for run - 80 characters
c
      read(5,100) title
  100 format(10a8)
      write(6,200) title
  200 format(/80('*')//10a8//80('*')//)
c
c     iphr--if input k-matrices are generated using the fano phase convention
c           set iphr equal to 1.  This corrects the differential cross section
c           for the fact that the spherical harmonics include the factor i**l.
c           For the Condon and Shortley phase convention set iphr equal to 0.
c           (default: 1)
c     iprint--if equal to 1, print k-matrices; if equal to 2, print
c             t-matrices; if equal to 3, print ql and a-lambda 
c             coefficents; if equal to 4, print m array.
c             (default:0)
c     icharg--charge of target
c     llmx--maximum total angular momentum for partial waves to be included
c           in the cross section calculation.
c     nt--number of target terms included in input k-matrices
c     iti--term number of the initial state from list of target terms
c          must be in the range 1<=iti<=nt
c     ntf--number of final states considered
c     imrd--set equal to zero to read k-matrices from unitarized
c           distorted-wave run or R-matrix run from file KMAT.DAT
c           set equal to 1 to read t-matrices from non-unitarized
c           distorted-wave run from file tmatnu
c           (default:0)
c     nthta--number of angles for thmn<theta<=thmx - must be even
c           (default:180)
c     thmn--minimum angle for calculation
c           (default:0.0)
c     thmx--maximum angle for calculation
c           (default:180.)
c     ein--incident energy in rydbergs referred to the ground term
c          must be identical to one of the energies included in the file
c          KMAT.DAT or tmatnu     
c     iunit--unit for differential cross section plot file; set equal to
c            zero for cross sections in units a0**2/sr; 
c            set equal to 1 for cross sections in units of 10-19 cm2/sr; 
c            set equal to 2 for cross sections in units of 10-18 cm2/sr; 
c            set equal to 3 for cross sections in units of 10-17 cm2/sr; 
c            set equal to 4 for cross sections in units of 10-16 cm2/sr.
c            (default:0)
c
      iphr=1
      iprint=0
      imrd=0
      nthta=180
      thmn=zero
      thmx=one80
      iunit=0
c
      read(5,difls)
      if(thmx.eq.0.0) thmx=one80
      if(nt.gt.ntm) stop 'increase ntm to total number of terms'
      if(ntf.gt.nfm) stop 'increase nfm to number of final terms'
      if(nthta.gt.ntem-1) stop 'increase ntem to number of angles'
      if(mod(nthta,2).ne.0) nthta=nthta-1
      write(6,210) iphr,iprint,icharg,llmx,nt,imrd,iunit,iunit,nthta,
     1 thmn,thmx,ein
  210 format(5x,'iphr=',i3,5x,'iprint=',i3,5x,'charge of target =',i3/
     1 5x,'maximum total angular momentum of N+1-electron system =',i3/
     2 5x,'total number of target terms =',i3,5x,'imrd =',i3/
     3 5x,'iunit =',i3,5x,'iunit =',i3,5x,'number of angles considered '
     4 ,i3/5x,'minimum angle =',f6.1,5x,'maximum angle =',f6.1/
     5 5x,'incident energy relative to ground state =',e15.6,
     6 ' rydbergs')
c
c     itf(i)--term numbers of final states chosen from list of
c             target terms - must all be in the range iti<=itf(i)<=nt
c
      read(5,130) (itf(i),i=1,ntf)
  130 format(16i5)
      ielas=0
      do 15 i=1,ntf
      if(itf(i).lt.iti) stop 4
      if(itf(i).eq.iti) ielas=1
   15 continue
c
c     name(i)--name of ith target term - 10 characters
c     lt(i)--total angular momentum of ith target term
c     lspt(i)--2s+1 for the ith target term
c     lpt(i)--parity of ith target term (0-even, 1-odd)
c     ent(i)--energy of ith target term in rydbergs - terms must be
c             energy ordered
c
      read(5,120) name(iti),lt(iti),lspt(iti),lpt(iti),ent(iti)
  120 format(a,3i5,f10.5)
      write(6,220)
  220 format(////27x,'initial term for this run')
      write(6,230) iti,name(iti),lt(iti),lspt(iti),lpt(iti),
     1 ent(iti)
  230 format(i2,'.',1x,a,1x,'l=',i2,1x,'2s+1=',i2,1x,'parity=',
     1 i2,1x,'energy(ryd)=',f10.5)
      do i=1,ntf
        j=itf(i)
        read(5,120) name(j),lt(j),lspt(j),lpt(j),ent(j)
      end do
c
c     calculate incident energy referred to initial state iti
c
      einc=ein-ent(iti)
      write(6,235) einc
  235 format(/5x,'initial energy relative to this term =',f10.5,
     1 ' rydbergs')
      write(6,240)
  240 format(////27x,'final terms for this run')
      do i=1,ntf
        j=itf(i)
        write(6,230) j,name(j),lt(j),lspt(j),lpt(j),ent(j)
      end do
      write(6,260)
  260 format(//)
      delth=(thmx-thmn)/float(nthta)
      theta(1)=thmn
      if(thmn.eq.zero) theta(1)=delth/onehndrd
      theta(2)=thmn+delth
      nthmx=nthta+1
      do 55 i=3,nthmx
      theta(i)=theta(i-1)+delth
   55 continue
      lti=lt(iti)
      if(icharg.ne.0) then
        nlc=llmx+lti+1
        call coulph(einc,icharg,nlc,sigbri)
      end if
      lspin=lspt(iti)
      if(lspin.gt.1) go to 60
      nspin=1
      ispin(1)=2
      go to 65
   60 nspin=2
      ispin(1)=lspin-1
      ispin(2)=lspin+1
   65 do 67 i=1,ntf
      j=itf(i)
      ltf(i)=lt(j)
      if(icharg.ne.0) then
       efnc=ein-ent(j)
       nlc=llmx+ltf(i)+1
       call coulph(efnc,icharg,nlc,sigbrf(1,i))
      end if
   67 continue
c
c     calculate the logarithms of factorials and store them in the
c     array gamma - used in calculation of 3-j and 6-j symbols
c
      gamma(1)=0.0
      gamma(2)=0.0
      do 80 k=3,nftm
      j=k-1
      a=j
      gamma(k)=gamma(j)+log(a)
   80 continue
      call tmtrix(ein)
      call mqaray
      call alcoef
      call cross
c
c     print cross sections
c
      write(6,272)
  272 format(//////30x,30('-')//30x,'*****   cross sections   *****',
     1 //30x,30('-')///)
      do 125 nf=1,ntf
      j=itf(nf)
      if(itf(nf).eq.iti) then
        write(6,275) iti
  275   format(////25x,'*****',i5,' - target final state - elastic',
     1   ' channel   *****')
      else
        write(6,280) j
  280   format(////30x,'*****',i5,' - target final state   *****')
      end if
      write(6,285)
  285 format(//44x,'angular distribution'//12x,'theta(degrees)',
     1 9x,'differential cross section (au/sr)',6x,'differential cross',
     2 ' section (cm2/sr)'//)
      write(7,286) iti,j,pltunt(iunit+1)
  286 format('differential cross section from term ',i3,' to term ',i3,
     1 ' in units of ',a)
      write(7,287)
  287 format('angle (deg)',5x,'  diffsig')
      do 90 nth=1,nthmx
      dsgau=dsg(nth,nf)
      dsgcm=convsg*dsgau
      dsgplt=dsgau
      if(iunit.eq.1) dsgplt=dsgcm/convsg1
      if(iunit.eq.2) dsgplt=dsgcm/convsg2
      if(iunit.eq.3) dsgplt=dsgcm/convsg3
      if(iunit.eq.4) dsgplt=dsgcm/convsg4
      nath=nth-1
      write(6,290) nath,theta(nth),dsgau,dsgcm
  290 format(1x,i4,9x,f10.4,19x,1pe16.8,24x,e16.8)
      write(7,295) theta(nth),dsgplt
  295 format(f6.2,11x,1pe11.5)
   90 continue
      if(icharg.ne.0.and.itf(nf).eq.iti) go to 125
      sigau=sig(nf)
      sigcm2=convsg*sigau
      write(6,297) sigau,sigcm2
  297 format(//21x,'total cross section = ',1pe16.8,' au',21x,e16.8,
     1 ' cm2')
      sigicm=convsg*sigic(nf)
      write(6,300) sigic(nf),sigicm
  300 format(/16x,'integrated cross section = ',1pe16.8,' au',21x,
     1 e16.8,' cm2')
      if(nipc.ne.1) then
      write(6,310)
  310 format(//11x,'partial cross section fraction from 0 to final',
     1 ' angle'//5x,'final angle(degrees)',10x,'partial cross section',
     2 ' fraction (%)'//)
      do 115 ipc=1,nipc
      write(6,320) thpcf(ipc),pcf(ipc,nf)
  320 format(13x,f4.0,37x,f5.1)
  115 continue
      end if
  125 continue
      end
c
c     ******************************************************************
c
      subroutine tmtrix(ein)
c
c     read the k-matrices and calculate and store the t-matrices
c
      implicit real*8(a-h,o-z)
c
      parameter(ntm=100,nfm=10,ncm=200,nctm=10,nlm=71,njm=2*nlm,
     1 nlbm=2*nlm,nwm=102,nftm=800,ntem=181)
      common/blok1/einc,theta(ntem),sigbri(nlm),sigbrf(nlm,nfm),
     1 pcf(61,nfm),thpcf(61),thmn,thmx,nthta,iprint,icharg,llmx,ntf,
     2 iti,itf(ntm),nthmx,lti,nspin,ispin(2),ltf(nfm),limn,limx,
     3 lfmx(nlm,nfm),lfmn(nlm,nfm),lspin,ielas,iphr,nipc,nfpc
      common/blok2/gamma(nftm)
      common/blok3/tr(nwm,ncm,ncm),ti(nwm,ncm,ncm),ltot(nwm),lstot(nwm),
     1 ncci(nwm),nccf(nwm,nfm),icci(nwm,nctm),iccf(nwm,nctm,nftm),
     2 lcs(nwm,ncm),iwu,imrd
      common/blok4/fmr(nlm,nfm,nlm,njm,2),fmi(nlm,nfm,nlm,njm,2),
     1 qr(nlm),qi(nlm),alm(nfm,nlbm),dsg(ntem,nfm),sig(nfm),sigic(nfm)
      common/blok5/pl(nlbm),imbmx
      dimension r(ncm,ncm),tnur(ncm,ncm),tnui(ncm,ncm),a(ncm,ncm),
     1 nchf(nfm),ichi(nctm),ichf(nctm,nfm),itm(ncm),lc(ncm),ipivot(ncm)
      data zero,one,two,etol/0.0d0,1.0d0,2.0d0,1.0d-04/
c
      limn=1000
      limx=0
      do 5 nf=1,nfm
      do 4 i=1,nlm
      lfmn(i,nf)=1000
      lfmx(i,nf)=0
    4 continue
    5 continue
c
c     loop over the partial waves
c
      iwu=0
      if(imrd.eq.0) then
        open(10,file='KMAT.DAT')
        iread=10
      else
        open(11,file='tmatnu')
        iread=11
      end if
    7 read(iread,*,end=110) nchop,istot,ltt,eksq
c  200 format(3i5,1pe14.8)                   !independent of stgf format
      if(nchop.gt.ncm) then
        write(6,205) nchop
  205   format('increase ncm to at least ',i5)
        stop 'increase ncm'
      end if
      lstt=iabs(istot)
      lptot=0
      if(istot.lt.0) lptot=1
      lcmx=0
      do 10 ic=1,nchop
      read(iread,210) itm(ic),lc(ic)
  210 format(2i5)
      if(lc(ic).gt.lcmx) lcmx=lc(ic)
   10 continue
      nlcmx=lcmx+1
      if(nlcmx.gt.nlm) then
        write(6,215) nlcmx
  215   format('increase nlm to at least ',i3)
        stop 'increase nlm'
      end if
c
c     if imrd equal to zero, read the k-matrix for each partial wave
c
      if(imrd.eq.0) then
        do i=1,nchop
          do ip=1,nchop
            read(10,220) r(i,ip)
  220       format(e15.6)
          end do
        end do
      else
c
c     if imrd not equal to zero, read the real and imaginary parts of
c     the t-matrices
c
        do i=1,nchop
          do ip=1,nchop
            read(11,230) tnur(i,ip),tnui(i,ip)
  230       format(2e15.6)
          end do
        end do
      end if
c
c     first determine whether this matrix is for the desired energy
c
      etst=two*abs(ein-eksq)/(ein+eksq)
      if(etst.gt.etol) go to 7
c
c      then determine whether ltt for this matrix exceeds llmx
c
      if(ltt.gt.llmx) go to 7
c
c     next determine whether this matrix contains the terms of
c     interest
c
      do 25 ins=1,nspin
      if(ispin(ins).eq.lstt) go to 30
   25 continue
      go to 7
   30 k=0
      do 35 ic=1,nchop
      if(itm(ic).ne.iti) go to 35
      k=k+1
      ichi(k)=ic
   35 continue
      if(k.eq.0) go to 7
      nchi=k
      if(nchi.gt.nctm) then
        write(6,235) nchi
  235   format('increase nctm to at least ',i3)
        stop 'increase nctm'
      end if
      ifo=0
      nchfmx=0
      do 45 nf=1,ntf
      k=0
      do 40 ic=1,nchop
      if(itm(ic).ne.itf(nf)) go to 40
      ifo=1
      k=k+1
      ichf(k,nf)=ic
   40 continue
      nchf(nf)=k
      if(nchf(nf).gt.nchfmx) nchfmx=nchf(nf)
   45 continue
      if(ifo.eq.0) go to 7
      if(nchfmx.gt.nctm) then
        write(6,235) nchfmx
        stop 'increase nctm'
      end if
c
c     store the necessary information
c
      iwu=iwu+1
      if(iwu.gt.nwm)then
        stop 'increase nwm'
      end if
      ltot(iwu)=ltt
      lstot(iwu)=lstt
      ncci(iwu)=nchi
      do i=1,nchi
        icci(iwu,i)=ichi(i)
      end do
      do nf=1,ntf
        nccf(iwu,nf)=nchf(nf)
        imx=nchf(nf)
        do i=1,imx
          iccf(iwu,i,nf)=ichf(i,nf)
        end do
      end do
      do ic=1,nchop
        lcs(iwu,ic)=lc(ic)
      end do
c
c     print the k-matrix
c
      if(iprint.eq.1.and.imrd.eq.0) then
        write(6,300) ltt,lstt,lptot
  300   format(//44x,'k-matrix for l =',i3,'  2s+1 =',i3,
     1   '  parity =',i3/44x,43('*'))
        do i=1,nchop
          write(6,310) i,itm(i),lc(i),(r(i,ip),ip=1,nchop)
  310     format(1x,i2,2i3,1pe12.4,8e13.4/(9x,e12.4,8e13.4))
        end do
      end if
c
c     if imrd equals zero, calculate the real and imaginary parts of
c     the t-matrix
c
c
c     t = -2i*r/(1-ir)
c     tr = 2*r*(1+r**2)^(-1)*r
c     ti= -2*r*(1+r**2)^(-1)
c
      if(imrd.eq.0) then
        do j=1,nchop
          do i=1,nchop
            a(i,j)=zero  
          end do
          do k=1,nchop
            do i=1,nchop
              a(i,j)=a(i,j)+r(i,k)*r(k,j)
            end do
          end do
        end do
	do i=1,nchop
	  a(i,i)=a(i,i)+one
	end do
        call vert(a,ncm,nchop,ipivot,ierr)
	if(ierr.ne.0) then
	  write(6,400)
  400     format(' the matrix: 1 + k^2 has no inverse in subroutine '
     1           ,'tmtrix - must stop')
          stop 'stop because no inverse found in calculating tmtrix'
        end if
        do j=1,nchop
          do i=1,nchop
            tnui(i,j)=zero
          enddo
          do k=1,nchop
            do i=1,nchop
              tnui(i,j)=tnui(i,j)+r(i,k)*a(k,j)
            end do
          enddo
          do i=1,nchop
            tnui(i,j)=-two*tnui(i,j)
          end do
        end do
        do j=1,nchop
          do i=1,nchop
            tnur(i,j)=-two*a(i,j)
          enddo
        enddo
        do j=1,nchop
          tnur(j,j)=tnur(j,j)+two
        end do      
      end if
      do ip=1,nchop
        do i=1,nchop
          tr(iwu,i,ip)=tnur(i,ip)
          ti(iwu,i,ip)=tnui(i,ip)
        end do
      end do
c
c     print the real and imaginary parts of the t-matrix
c
      if(iprint.eq.2) then
        write(6,320) ltt,lstt,lptot
  320   format(//44x,'t-matrix for l =',i3,'  2s+1 =',i3,
     1   '  parity =',i3/44x,43('*'))
        do i=1,nchop
          write(6,330) i,itm(i),lc(i),(tr(iwu,i,ip),ti(iwu,i,ip),
     1    ip=1,nchop)
  330     format(1x,i2,2i3,4(1pe13.4,' +i*',e12.4)/9x,4(e13.4,
     1     ' +i*',e12.4))
        end do
      end if
c
c     determine the minimum and maximum value of li, and for each final
c     term and each value of li, determine the maximum and minimum
c     values of lf
c
      do 95 nci=1,nchi
      ici=ichi(nci)
      li=lc(ici)
      if(li.lt.limn) limn=li
      if(li.gt.limx) limx=li
      ii=li+1
      do 90 nf=1,ntf
      if(nchf(nf).eq.0) go to 90
      ncfmx=nchf(nf)
      do 87 ncf=1,ncfmx
      icf=ichf(ncf,nf)
      lf=lc(icf)
      if(lf.lt.lfmn(ii,nf)) lfmn(ii,nf)=lf
      if(lf.gt.lfmx(ii,nf)) lfmx(ii,nf)=lf
   87 continue
   90 continue
   95 continue
      go to 7
  110 return
      end
c
c     ******************************************************************
c
      subroutine mqaray
c
c     calculate the real and imaginary parts of the q-array needed for
c     the cross terms for elastic scattering and the the real and
c     imaginary parts of the m-array needed for the calculation of the
c     a-lambda coefficients
c
      implicit real*8(a-h,o-z)
c
      parameter(ntm=100,nfm=10,ncm=200,nctm=10,nlm=71,njm=2*nlm,
     1 nlbm=2*nlm,nwm=102,nftm=800,ntem=181)
      common/blok1/einc,theta(ntem),sigbri(nlm),sigbrf(nlm,nfm),
     1 pcf(61,nfm),thpcf(61),thmn,thmx,nthta,iprint,icharg,llmx,ntf,
     2 iti,itf(ntm),nthmx,lti,nspin,ispin(2),ltf(nfm),limn,limx,
     3 lfmx(nlm,nfm),lfmn(nlm,nfm),lspin,ielas,iphr,nipc,nfpc
      common/blok2/gamma(nftm)
      common/blok3/tr(nwm,ncm,ncm),ti(nwm,ncm,ncm),ltot(nwm),lstot(nwm),
     1 ncci(nwm),nccf(nwm,nfm),icci(nwm,nctm),iccf(nwm,nctm,nftm),
     2 lcs(nwm,ncm),iwu,imrd
      common/blok4/fmr(nlm,nfm,nlm,njm,2),fmi(nlm,nfm,nlm,njm,2),
     1 qr(nlm),qi(nlm),alm(nfm,nlbm),dsg(ntem,nfm),sig(nfm),sigic(nfm)
      common/blok5/pl(nlbm),imbmx
c
      data zero,one,two/0.0d0,1.0d0,2.0d0/
c
c     initialize the arrays
c
      do ii=1,nlm
        qr(ii)=0.0
        qi(ii)=0.0
      end do
      do ins=1,nspin
        do ij=1,njm
	  do il=1,nlm
	    do it=1,nfm
	      do ii=1,nlm
	        fmr(ii,it,il,ij,ins)=zero
                fmi(ii,it,il,ij,ins)=zero
	      end do
	    end do
	  end do
        end do
      end do
c
c     loop over partial waves
c
      do 130 iw=1,iwu
      do 70 ins=1,nspin
      if(ispin(ins).eq.lstot(iw)) go to 80
   70 continue
   80 fmlt=two*float(ltot(iw))+one
      fmst=float(lstot(iw))
      coeff=fmlt*fmst
      coef1=(-1)**ltot(iw)*sqrt(fmst)*fmlt
      nchi=ncci(iw)
c
c     calculate the contribution to the real and imaginary parts of the
c     array q(li)
c
      if(ielas.ne.0.and.icharg.ne.0) then
      do 90 nci=1,nchi
      ic=icci(iw,nci)
      li=lcs(iw,ic)
      ii=li+1
      qr(ii)=qr(ii)+coeff*(cos(2.*sigbri(ii))*tr(iw,ic,ic)-
     1 sin(2.*sigbri(ii))*ti(iw,ic,ic))
      qi(ii)=qi(ii)+coeff*(cos(2.*sigbri(ii))*ti(iw,ic,ic)+
     1 sin(2.*sigbri(ii))*tr(iw,ic,ic))
   90 continue
      end if
c
c     calculate the contributions to the real and imaginary parts of
c     array m
c
      do 120 nci=1,nchi
      ici=icci(iw,nci)
      li=lcs(iw,ici)
      fmli=2.*float(li)+1.
      ii=li+1
      do 115 nf=1,ntf
      if(nccf(iw,nf).eq.0) go to 115
      ncfmx=nccf(iw,nf)
      lft=ltf(nf)
      id1=iabs(lft-lti)
      is1=lft+lti
      do 110 ncf=1,ncfmx
      icf=iccf(iw,ncf,nf)
      lf=lcs(iw,icf)
      fmlf=2.*float(lf)+1.
      il=lf-lfmn(ii,nf)+1
      id2=iabs(lf-li)
      ijtmn=max0(id1,id2)+1
      is2=lf+li
      ijtmx=min0(is1,is2)+1
      coef2=(-1)**(li+lf)*coef1*sqrt(fmli*fmlf)
      do 105 ijt=ijtmn,ijtmx
      jt=ijt-1
      inj=ijt-ijtmn+1
      coe6j=coef2*s6j(lti,lft,jt,lf,li,ltot(iw))
      fmr(ii,nf,il,inj,ins)=fmr(ii,nf,il,inj,ins)+coe6j*tr(iw,ici,icf)
      fmi(ii,nf,il,inj,ins)=fmi(ii,nf,il,inj,ins)+coe6j*ti(iw,ici,icf)
  105 continue
  110 continue
  115 continue
  120 continue
  130 continue
      if (iprint.eq.4) then
        write(6,200)
  200   format(/)
        do ins=1,nspin
          do inj=1,njm
	    do il=1,nlm
	      do nf=1,nfm
	        do ii=1,nlm
		  if(fmr(ii,nf,il,inj,ins).ne.zero.or.
     1             fmi(ii,nf,il,inj,ins).ne.zero) then
                    write(6,210) ii,nf,il,inj,ins,fmr(ii,nf,il,inj,ins),
     1                           fmi(ii,nf,il,inj,ins)
  210               format(2x,'ii =',i2,2x,'nf =',i2,2x,'il =',i2,2x,
     1                     'inj =',i2,2x,'ins =',i2,2x,'fmr =',
     2                      1pe15.7/2x,'fmi =',e15.7)
                  end if  
                end do
	      end do
	    end do
	  end do
        end do
        write(6,200)
      end if 
      return
      end
c
c     ******************************************************************
c
      subroutine alcoef
c
c     calculate the a-lambda coefficients
c
      implicit real*8(a-h,o-z)
c
      parameter(ntm=100,nfm=10,ncm=200,nctm=10,nlm=71,njm=2*nlm,
     1 nlbm=2*nlm,nwm=102,nftm=800,ntem=181)
      common/blok1/einc,theta(ntem),sigbri(nlm),sigbrf(nlm,nfm),
     1 pcf(61,nfm),thpcf(61),thmn,thmx,nthta,iprint,icharg,llmx,ntf,
     2 iti,itf(ntm),nthmx,lti,nspin,ispin(2),ltf(nfm),limn,limx,
     3 lfmx(nlm,nfm),lfmn(nlm,nfm),lspin,ielas,iphr,nipc,nfpc
      common/blok2/gamma(nftm)
      common/blok3/tr(nwm,ncm,ncm),ti(nwm,ncm,ncm),ltot(nwm),lstot(nwm),
     1 ncci(nwm),nccf(nwm,nfm),icci(nwm,nctm),iccf(nwm,nctm,nftm),
     2 lcs(nwm,ncm),iwu,imrd
      common/blok4/fmr(nlm,nfm,nlm,njm,2),fmi(nlm,nfm,nlm,njm,2),
     1 qr(nlm),qi(nlm),alm(nfm,nlbm),dsg(ntem,nfm),sig(nfm),sigic(nfm)
      common/blok5/pl(nlbm),imbmx
      dimension alim(nfm,nlbm)
      data zero,one,two/0.0d0,1.0d0,2.0d0/
      lmbmx=2*limx
      imbmx=lmbmx+1
      iimn=limn+1
      iimx=limx+1
c
c     alm is the true a-lambda coefficient used in the calculation of
c     of the differential cross section and is, of course, real.  alim
c     is the imaginary part of the a-lambda coefficient and should, of
c     course, vanish.  it is calculated only as a test of the code.
c
c     loop over final terms
c
      do 200 nf=1,ntf
      id1=iabs(ltf(nf)-lti)
      is1=ltf(nf)+lti
c
c     begin loop over lambda
c
      do 190 imb=1,imbmx
      lmb=imb-1
      fmlmb=two*float(lmb)+one
      alm(nf,imb)=zero
      alim(nf,imb)=zero
      do 180 ii=iimn,iimx
      li=ii-1
      ifmn=lfmn(ii,nf)+1
      ifmx=lfmx(ii,nf)+1
      iipmn=iabs(lmb-li)+1
      iipmn=max0(iipmn,iimn)
      iipmx=min0(iimx,(lmb+li+1))
      do 170 iip=iipmn,iipmx
      lip=iip-1
      if(mod(lip+li+lmb,2).ne.0) go to 170
      ifpmn1=lfmn(iip,nf)+1
      ifpmx1=lfmx(iip,nf)+1
      do 160 if=ifmn,ifmx
      lf=if-1
      inf=lf-lfmn(ii,nf)+1
      id2=iabs(lf-li)
      id2=max0(id2,id1)
      is2=lf+li
      is2=min0(is2,is1)
      ifpmn=iabs(lmb-lf)+1
      ifpmn=max0(ifpmn1,ifpmn)
      ifpmx=min0(ifpmx1,(lmb+lf+1))
      do 150 ifp=ifpmn,ifpmx
      lfp=ifp-1
      if(mod(lfp+lf+lmb,2).ne.0) go to 150
      infp=lfp-lfmn(iip,nf)+1
      id2p=iabs(lfp-lip)
      id2p=max0(id2p,id1)
      ijtmn=max0(id2,id2p)+1
      is2p=lfp+lip
      ijtmx=min0(is2,is2p)+1
c
c     if iphr equal to 1, do not include phase factor
c     (-1)**((li+lfp-lip-lf)/2) to correct for the fano phase convention
c
      coeff=s3j0(li,lip,lmb)*s3j0(lf,lfp,lmb)
      if(iphr.ne.1) coeff=(-1)**((li+lfp-lip-lf)/2)*coeff
      coef1=coeff
      coef2=zero
      if(icharg.ne.0) then
      pdi=sigbri(ii)-sigbri(iip)
      pdf=sigbrf(if,nf)-sigbrf(ifp,nf)
      cs=cos(pdi+pdf)
      sn=sin(pdi+pdf)
      coef1=coeff*cs
      coef2=coeff*sn
      end if
      sum1=zero
      sum2=zero
      do 140 ijt=ijtmn,ijtmx
      jt=ijt-1
      inj=ijt-id2
      injp=ijt-id2p
      fmjt=two*float(jt)+one
      coef3=(-1)**(jt+lmb)*fmjt*fmlmb
      isa=0
      do 130 is=1,nspin
      fr=fmr(ii,nf,inf,inj,is)
      frp=fmr(iip,nf,infp,injp,is)
      fi=fmi(ii,nf,inf,inj,is)
      fip=fmi(iip,nf,infp,injp,is)
      f1a=fr*frp
      f1b=fi*fip
      f2a=fr*fip
      f2b=frp*fi
      if(f1a.eq.0.0.and.f1b.eq.0.0.and.f2a.eq.0.0.and.f2b.eq.0.0)
     1 go to 130
      isa=isa+1
      if(isa.eq.1) coef3=coef3*s6j(li,lf,jt,lfp,lip,lmb)
      sum1=sum1+coef3*(f1a+f1b)
      if(icharg.ne.0) sum2=sum2+coef3*(f2a-f2b)
  130 continue
  140 continue
      alm(nf,imb)=alm(nf,imb)+coef1*sum1+coef2*sum2
      alim(nf,imb)=alim(nf,imb)+coef2*sum1-coef1*sum2
  150 continue
  160 continue
  170 continue
  180 continue
  190 continue
  200 continue
      if(iprint.eq.3) then
c
c     print the ql and a-lambda coefficients
c
      if(icharg.ne.0.and.ielas.ne.0) then
      write(6,300) limx
  300 format(////51x,'ql coefficients up to li = ',i3/)
      iimx=limx+1
      write(6,310) (qr(i),qi(i),i=1,iimx)
  310 format((8x,4(e13.4,' +i*',e12.4)))
      end if
      write(6,320)
  320 format(//)
      do 210 nf=1,ntf
      write(6,330) itf(nf),lmbmx
  330 format(//31x,'a-lambda coefficient for term number ',i3,
     1 ' up to a maximum lambda of ',i3/)
      write(6,340) (alm(nf,imb),imb=1,imbmx)
  340 format((6x,8e15.6))
      write(6,350) itf(nf),lmbmx
  350 format(//31x,'ai-lambda coefficient for term number ',i3,
     1 ' up to a maximum lambda of ',i3/)
      write(6,340) (alim(nf,imb),imb=1,imbmx)
  210 continue
      end if
      return
      end
c
c     ******************************************************************
c
      subroutine cross
c
c     calculate the excitation cross sections as a function of angle and
c     the partial and total cross sections
c
      implicit real*8(a-h,o-z)
c
      parameter(ntm=100,nfm=10,ncm=200,nctm=10,nlm=71,njm=2*nlm,
     1 nlbm=2*nlm,nwm=102,nftm=800,ntem=181)
      common/blok1/einc,theta(ntem),sigbri(nlm),sigbrf(nlm,nfm),
     1 pcf(61,nfm),thpcf(61),thmn,thmx,nthta,iprint,icharg,llmx,ntf,
     2 iti,itf(ntm),nthmx,lti,nspin,ispin(2),ltf(nfm),limn,limx,
     3 lfmx(nlm,nfm),lfmn(nlm,nfm),lspin,ielas,iphr,nipc,nfpc
      common/blok2/gamma(nftm)
      common/blok3/tr(nwm,ncm,ncm),ti(nwm,ncm,ncm),ltot(nwm),lstot(nwm),
     1 ncci(nwm),nccf(nwm,nfm),icci(nwm,nctm),iccf(nwm,nctm,nftm),
     2 lcs(nwm,ncm),iwu,imrd
      common/blok4/fmr(nlm,nfm,nlm,njm,2),fmi(nlm,nfm,nlm,njm,2),
     1 qr(nlm),qi(nlm),alm(nfm,nlbm),dsg(ntem,nfm),sig(nfm),sigic(nfm)
      common/blok5/pl(nlbm),imbmx
c
      data zero,one,two,three,four,eight,fourty5,onehndrd,one80/
     1 0.0d0,1.0d0,2.0d0,3.0d0,4.0d0,8.0d0,4.5d1,1.00d2,1.80d2/
c
c     first calculate the differential cross sections as a function
c     of scattering angle
c
      fmlti=two*float(lti)+one
      fmsti=float(lspin)
      pi=acos(-one)
      coeff=pi/(two*fmlti*fmsti*einc)
      coef3=one/(eight*fmlti*fmsti*einc)
      cnvrd=atan(one)/fourty5
      if(icharg.ne.0) then
      alpha=float(icharg)/sqrt(einc)
      coef1=alpha**2/(four*einc)
      coef2=-alpha/(four*fmlti*fmsti*einc)
      end if
c
c     loop over angles
c
      do 40 nth=1,nthmx
      ang=cnvrd*theta(nth)
      call legpol(ang)
c
c     loop over the final terms
c
      do 30 nf=1,ntf
      dsg1=zero
      dsg2=zero
      if(itf(nf).eq.iti.and.icharg.ne.0) then
c
c     elastic scattering for a charged target
c
      snth2=sin(ang/two)**2
      dsg1=coef1/(snth2**2)
      tht=alpha*log(snth2)
      cfi=coef2*cos(tht)
      cfr=-coef2*sin(tht)
      sumr=qr(1)
      sumi=qi(1)
      ilmx=limx+1
      do 10 il=2,ilmx
      l=il-1
      sumr=sumr+qr(il)*pl(l)
      sumi=sumi+qi(il)*pl(l)
   10 continue
      dsg2=(cfr*sumr+cfi*sumi)/snth2
      end if
      sum=alm(nf,1)
      do 20 imb=2,imbmx
      lmb=imb-1
      sum=sum+alm(nf,imb)*pl(lmb)
   20 continue
      dsg3=coef3*sum
      dsg(nth,nf)=dsg1+dsg2+dsg3
   30 continue
   40 continue
c
c     calculate total cross section for each inelastic transition
c
      nipc=1
      nfpc=nthmx
      if(thmn.eq.zero.and.thmx.eq.one80.and.nthta.eq.180) then
      nipc=61
      nfpc=61
      nthf=nfpc-2
      do 42 ipc=1,nipc
      nthf=nthf+2
      thpcf(ipc)=theta(nthf)
   42 continue
      end if
      do 70 nf=1,ntf
      if(icharg.ne.0.and.itf(nf).eq.iti) go to 70
      sig(nf)=coeff*alm(nf,1)
c
c     calculate the partial cross section from thmn to thmx by direct
c     integration of the differential cross section
c      when the angular range is from 0 to 180, also calculate the partial
c      cross section fractions from 0 to a final theta starting a 60 and
c      carrying on through 180
c
      h=cnvrd*theta(2)
      angi=cnvrd*theta(1)
      nthf=nfpc-2
      do 50 ipc=1,nipc
      nthf=nthf+2
      angf=cnvrd*theta(nthf)
      sum=sin(angi)*dsg(1,nf)+sin(angf)*dsg(nthf,nf)
      fact=two
      imx=nthf-1
      do 45 i=2,imx
      if(fact.eq.two) then
      fact=four
      else
      fact=two
      end if
      ang=cnvrd*theta(i)
      sum=sum+fact*sin(ang)*dsg(i,nf)
   45 continue
      pcf(ipc,nf)=two*pi*sum*h/three
      if(ipc.eq.nipc) sigic(nf)=pcf(ipc,nf)
   50 continue
      if(nipc.ne.1) then
      do 60 ipc=1,nipc
      pcf(ipc,nf)=(pcf(ipc,nf)/sigic(nf))*onehndrd
   60 continue
      end if
   70 continue
      return
      end
c
c     ******************************************************************
c
      subroutine coulph(e,ic,n,sigbr)
c
c     calculate the coulomb phase shift differences:
c     sigbr(l)=sigma(l)-sigma(0)
c
      implicit real*8(a-h,o-z)
c
      parameter(nlm=71)
      dimension sigbr(nlm)
c
      data zero/0.0d0/
c
      alpha=float(ic)/sqrt(e)
      sigbr(1)=zero
      a=zero
      do 1 l=2,n
      w=float(l-1)
      a=a-atan2(alpha,w)
      sigbr(l)=a
    1 continue
      return
      end
c
c     ******************************************************************
c
      subroutine vert(v,lv,n,w,ierr)
c
c      ________________________________________________________
c     |                                                        |
c     |                invert a general matrix                 |
c     |                                                        |
c     |    input:                                              |
c     |                                                        |
c     |         v     --array containing matrix                |
c     |                                                        |
c     |         lv    --leading (row) dimension of array v     |
c     |                                                        |
c     |         n     --dimension of matrix stored in array v  |
c     |                                                        |
c     |         w     --integer work array with at least n-1   |
c     |                      elements                          |
c     |                                                        |
c     |    output:                                             |
c     |                                                        |
c     |         v     --inverse                                |
c     |                                                        |
c     |    builtin functions: abs                              |
c     |________________________________________________________|
c
      real*8 v(lv,*),s,t
      integer w(*),i,j,k,l,m,n,p
c
      data zero,one/0.0d0,1.0d0/
      ierr=0
c
      if ( n .eq. 1 ) goto 110
      l = 0
      m = 1
10    if ( l .eq. n ) goto 90
      k = l
      l = m
      m = m + 1
c     ---------------------------------------
c     |*** find pivot and start row swap ***|
c     ---------------------------------------
      p = l
      if ( m .gt. n ) goto 30
      s = abs(v(l,l))
      do 20 i = m,n
           t = abs(v(i,l))
           if ( t .le. s ) goto 20
           p = i
           s = t
20    continue
      w(l) = p
30    s = v(p,l)
      v(p,l) = v(l,l)
      if ( s .eq. zero ) goto 120
c     -----------------------------
c     |*** compute multipliers ***|
c     -----------------------------
      v(l,l) = -one
      s = one/s
      do 40 i = 1,n
40         v(i,l) = -s*v(i,l)
      j = l
50    j = j + 1
      if ( j .gt. n ) j = 1
      if ( j .eq. l ) goto 10
      t = v(p,j)
      v(p,j) = v(l,j)
      v(l,j) = t
      if ( t .eq. zero ) goto 50
c     ------------------------------
c     |*** eliminate by columns ***|
c     ------------------------------
      if ( k .eq. 0 ) goto 70
      do 60 i = 1,k
60         v(i,j) = v(i,j) + t*v(i,l)
70    v(l,j) = s*t
      if ( m .gt. n ) goto 50
      do 80 i = m,n
80         v(i,j) = v(i,j) + t*v(i,l)
      goto 50
c     -----------------------
c     |*** pivot columns ***|
c     -----------------------
90    l = w(k)
      do 100 i = 1,n
           t = v(i,l)
           v(i,l) = v(i,k)
100        v(i,k) = t
      k = k - 1
      if ( k .gt. 0 ) goto 90
      return
110   if ( v(1,1) .eq. zero ) goto 120
      v(1,1) = one/v(1,1)
      return
120   ierr=1
      return
      end
c
c     *****************************************************************
c
      real*8 function s6j(j1,j2,j3,l1,l2,l3)
c
c         calculates 6-j symbol(j1,j2,j3,l1,l2,l3) using
c         a 'nested' algorithm
c
      implicit real*8(a-h,o-z)
c  
      parameter(nftm=800)
      common/blok2/gamma(nftm)
      data zero,one/0.0d0,1.0d0/
      s6j=zero
      i1=j1+j2+j3+2
      i2=j1+j2-j3+1
      i3=j1+j3-j2+1
      i4=j3+j2-j1+1
      if(i2.lt.1.or.i3.lt.1.or.i4.lt.1) go to 10
      d1=gamma(i2)+gamma(i3)+gamma(i4)-gamma(i1)
      i1=j1+l2+l3+2
      i2=j1+l2-l3+1
      i3=j1+l3-l2+1
      i4=l3+l2-j1+1
      if(i2.lt.1.or.i3.lt.1.or.i4.lt.1) go to 10
      d2=gamma(i2)+gamma(i3)+gamma(i4)-gamma(i1)
      i1=j2+l1+l3+2
      i2=j2+l1-l3+1
      i3=j2+l3-l1+1
      i4=l3+l1-j2+1
      if(i2.lt.1.or.i3.lt.1.or.i4.lt.1) go to 10
      d3=gamma(i2)+gamma(i3)+gamma(i4)-gamma(i1)
      i1=l2+l1+j3+2
      i2=l2+l1-j3+1
      i3=l2+j3-l1+1
      i4=j3+l1-l2+1
      if(i2.lt.1.or.i3.lt.1.or.i4.lt.1) go to 10
      d4=gamma(i2)+gamma(i3)+gamma(i4)-gamma(i1)
      delt=(d1+d2+d3+d4)*0.5
      k1=j1+j2+j3
      k2=j1+l2+l3
      k3=j2+l1+l3
      k4=l2+l1+j3
      k5=j1+j2+l2+l1
      k6=j2+l2+j3+l3
      k7=j1+l1+j3+l3
      lmin=max0(0,k1,k2,k3,k4)
      lmax=min0(k5,k6,k7)
      l=lmin+1
      lm=lmin-1
      m1=l+1
      m2=l-k1
      m3=l-k2
      m4=l-k3
      m5=l-k4
      m6=k5-lm
      m7=k6-lm
      m8=k7-lm
      delt=delt+gamma(m1)-gamma(m2)-gamma(m3)-gamma(m4)-gamma(m5)
     $-gamma(m6)-gamma(m7)-gamma(m8)
      a=one
      if(lmax.eq.lmin) go to 6
      jj=lmax
      do 1 j=l,lmax
      jk=jj-1
      c=(jj-k1)*(jj-k2)
      a=one-a*(jj+1)*(k5-jk)*(k6-jk)*(k7-jk)/(c*(jj-k3)*(jj-k4))
      jj=jk
    1 continue
    6 ish=iabs(mod(lm+1,2))
      ish=1-ish-ish
      s6j=a*exp(delt)*ish
   10 return
      end
c
c     *****************************************************************
c
      real*8 function s3j0(j1,j2,j3)
c
c         calculation of 3j-symbol (j1,j2,j3,0,0,0)
c
      implicit real*8(a-h,o-z)
c  
      parameter(nftm=800)
      common/blok2/gamma(nftm)
      data zero/0.0d0/
      j=j1+j2+j3
      ie=mod(j,2)
      ia=j-j1-j1
      ib=j-j2-j2
      ic=j-j3-j3
      s3j0=zero
      if(ie.ne.0) go to 10
      if(ia.lt.0.or.ib.lt.0.or.ic.lt.0) go to 10
      i1=j/2
      ish=mod(i1,2)
      ish=1-ish-ish
      i1=i1+1
      i2=j1+j2-j3+1
      i3=j1-j2+j3+1
      i4=-j1+j2+j3+1
      i5=j1+j2+j3+2
      x=(gamma(i2)+gamma(i3)+gamma(i4)-gamma(i5))*0.5
      i2=i1-j3
      i3=i1-j2
      i4=i1-j1
      x=x+gamma(i1)-gamma(i2)-gamma(i3)-gamma(i4)
      s3j0=exp(x)*ish
   10 return
      end
c
c     ******************************************************************
c
      subroutine legpol(t)
c
c         calculation of legendre polynomials by using a standard
c         recurrence relation
c         t = angle in radians
c
      implicit real*8(a-h,o-z)
c     
      parameter(nlm=71,nlbm=2*nlm)
      common/blok5/pl(nlbm),imbmx
      data one/1.0d0/
      x= cos(t)
      a=one
      b=x
      pl(1)=x
      do 1 l=2,imbmx
      l1=l-1
      c=(b*x*(l+l1)-l1*a)/l
      a=b
      b=c
      pl(l)=c
    1 continue
      return
      end