c
c     program difjk 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 between individual levels for electron-
c     atom and electron-ion scattering in the jk coupling scheme.
c
      implicit real*8(a-h,o-z)
c     parameter dimensions
c
c     nlvm--maximum number of target levels
c     nlvfm--maximum number of final target levels to which differential
c            cross sections are to be calculated
c     ncm--maximum number of channels in k-matrices
c     nclm--maximum number of channels for a given level
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(nlvm=120,nlvfm=10,ncm=620,nclm=20,nlm=46,njm=2*nlm,
     1 nlbm=2*nlm,nwm=102,nftm=800,ntem=181)
      character name(nlvm)*10,pltunt(5)*12
      dimension title(10),fj(nlvm),enj(nlvm),lpt(nlvm)
      common/blok1/thmn,thmx,einc,theta(ntem),fji,fjf(nlvfm),
     1 sigbri(nlm),sigbrf(nlm,nlvfm),nthta,iprint,icharg,j2mx,nlvf,
     2 ilvi,ilvf(nlvfm),nthmx,limn,limx,lfmx(nlm,nlvfm),
     3 lfmn(nlm,nlvfm),ielas,iphr,nipc,nfpc
      common/blok2/gamma(nftm)
      common/blok3/fjtot(nwm),tr(nwm,ncm,ncm),ti(nwm,ncm,ncm),
     1 flcs(nwm,ncm),fkvs(nwm,ncm),ncci(nwm),nccf(nwm,nlvfm),
     2 icci(nwm,nclm),iccf(nwm,nclm,nlvfm),iwu,imrd
      common/blok4/fmr(nlm,nlvfm,nlm,njm,4),fmi(nlm,nlvfm,nlm,njm,4),
     1 qr(nlm),qi(nlm),alm(nlvfm,nlbm),dsg(ntem,nlvfm),sig(nlvfm),
     2 sigic(nlvfm)
      common/blok5/pl(nlbm),imbmx
      namelist/difjk/iphr,iprint,icharg,j2mx,nlv,ilvi,nlvf,imrd,nthta,
     1 thmn,thmx,ein,iunit
      data zero,threehf,two,onehndrd,one80/0.0d0,1.5d0,2.0d0,1.00d2,
     1 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(unit=5,file='difjin')
      open(unit=6,file='difjout')
      open(unit=7,file='difjplt.dat')
c
c     input from difjin
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 r-matrices were generated with r-matrix code, set
c           iphr equal to 1.  this corrects the differential cross
c           section for the fano phase convention in which the spherical
c           harmonic includes the factor i**l.
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     j2mx--maximum value of 2xjtot to be included in the cross section
c           calculation
c     nlv--number of target levels included in input r-matrices
c     ilvi--level number of the initial level from list of target levels
c           must be in the range 1<=ilvi<=nlv
c     nlvf--number of final levels considered
c     imrd--set equal to zero to read k-matrices from unitarized
c           distorted-wave or close-coupling calculation
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 0<theta<=180 - 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 level
c          must be identical to one of the energies included in the
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
      read(5,difjk)
      fjmx=j2mx
      fjmx=fjmx/two
      if(thmx.eq.zero) thmx=one80 
      if(nlv.gt.nlvm) stop 'increase nlvm to total number of levels'
      if(nlvf.gt.nlvfm) 
     1 stop 'increase nlvfm to the number of final levels'
      if(nthta.gt.ntem-1) stop 'increase ntem to number of angles'
      if(mod(nthta,2).ne.0) nthta=nthta-1
      einev=ein*convev
      write(6,210) j2mx,icharg,nlv,ilvi,nlvf,nthta,iunit,einev
  210 format(2x,'2xJmax=',i3,5x,'charge of target =',i3,5x,
     1 'number of target levels =',i3,/2x,'initial level =',i3,5x,
     2 'number of final levels =',i3,5x,'number of angles =',i3/
     3  2x,'iunit =',i2,5x,'incident energy relative to ground term =',
     4  1pe15.6,' eV')
c
c     ilvf(i)--level numbers of final levels chosen from list of
c              target levels - must all be in the range:
c              ilvi<=ilvf(i)<=nlv
c
      read(5,130) (ilvf(i),i=1,nlvf)
      ielas=0
  130 format(16i5)
      do 15 i=1,nlvf
      if(ilvf(i).lt.ilvi) stop 4
      if(ilvf(i).eq.ilvi) ielas=1
   15 continue
c
c     name(i)--name of ith target term - 8 characters
c     fj(i)--total angular momentum of 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(ilvi),fj(ilvi),lpt(ilvi),enj(ilvi)
  120 format(a,f5.1,i5,f10.5)
      fji=fj(ilvi)
      enji=enj(ilvi)
      write(6,220)
  220 format(////27x,'initial level for this run')
      write(6,230) ilvi,name(ilvi),fji,lpt(ilvi),enj(ilvi)
  230 format(i2,'.',1x,a,1x,'j=',f5.1,1x,'parity=',i2,
     1 1x,'energy(ryd)=',f10.5)
      do i=1,nlvf
        j=ilvf(i)
	read(5,120) name(j),fjf(j),lpt(j),enj(j)
      end do
c
c     calculate incident energy referred to initial level
c
      einc=ein-enji
      write(6,235) einc
  235 format(/5x,'initial energy relative to this level =',f10.5,
     1 ' rydbergs')
      write(6,240)
  240 format(////27x,'final levels for this run')
      do i=1,nlvf
        j=ilvf(i)
	write(6,230) j,name(j),fjf(j),lpt(j),enj(j)
      end do
c
      write(6,260)
  260 format(//)
      fnthta=nthta
      delth=one80/fnthta
      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
      if(icharg.ne.0) then
        nlc=fjmx+fji+threehf
        call coulph(einc,icharg,nlc,sigbri)
      end if
      do 60 i=1,nlvf
      if(icharg.ne.0) then
        j=ilvf(i)  
        efnc=ein-enj(j)
        nlc=fjf(j)+fjmx+threehf
        call coulph(efnc,icharg,nlc,sigbrf(1,i))
      end if
   60 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)=zero
      gamma(2)=zero
      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 95 nf=1,nlvf
      j=ilvf(nf)
      if(j.eq.ilvi) then
        write(6,275) ilvi
  275   format(////21x,'*****',i5,' - target final level - elastic',
     1 ' channel   *****')
      else
        write(6,280) j
  280   format(////21x,'*****',i5,' - target final level   *****')
      end if
      write(6,285)
  285 format(//30x,'angular distribution'//10x,'angle(degrees)',
     1 5x,'differential sigma (au)',5x,'differential sigma (cm2)'//)
      write(7,286) ilvi,j,pltunt(iunit+1)
  286 format('differential cross section from level ',i3,' to level ',
     1 i3,' 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,5x,f10.4,11x,1pe16.8,11x,e16.8)
      write(7,295) theta(nth),dsgplt
  295 format(f6.2,11x,1pe11.5)
   90 continue
      if(icharg.ne.0.and.ilvf(nf).eq.ilvi) go to 95
      sigau=sig(nf)
      sigcm2=convsg*sigau
      write(6,297) sigau,sigcm2
  297 format(//9x,'total cross section = ',1pe16.8,' au',8x,
     1 e16.8,' cm2')
      sigicm=convsg*sigic(nf)
      write(6,300) sigic(nf),sigicm
  300 format(/4x,'integrated cross section = ',1pe16.8,' au',8x,
     1 e16.8,' cm2')
   95 continue
      end
c
c     ******************************************************************
c
      subroutine tmtrix(ein)
c
c     read the r-matrices and calculate and store the t-matrices
c
      implicit real*8(a-h,o-z)
      parameter(nlvm=120,nlvfm=10,ncm=620,nclm=20,nlm=46,njm=2*nlm,
     1 nlbm=2*nlm,nwm=102,nftm=800,ntem=181)
      common/blok1/thmn,thmx,einc,theta(ntem),fji,fjf(nlvfm),
     1 sigbri(nlm),sigbrf(nlm,nlvfm),nthta,iprint,icharg,j2mx,nlvf,
     2 ilvi,ilvf(nlvfm),nthmx,limn,limx,lfmx(nlm,nlvfm),
     3 lfmn(nlm,nlvfm),ielas,iphr,nipc,nfpc
      common/blok2/gamma(nftm)
      common/blok3/fjtot(nwm),tr(nwm,ncm,ncm),ti(nwm,ncm,ncm),
     1 flcs(nwm,ncm),fkvs(nwm,ncm),ncci(nwm),nccf(nwm,nlvfm),
     2 icci(nwm,nclm),iccf(nwm,nclm,nlvfm),iwu,imrd
      common/blok4/fmr(nlm,nlvfm,nlm,njm,4),fmi(nlm,nlvfm,nlm,njm,4),
     1 qr(nlm),qi(nlm),alm(nlvfm,nlbm),dsg(ntem,nlvfm),sig(nlvfm),
     2 sigic(nlvfm)
      dimension fjv(ncm),fkv(ncm),r(ncm,ncm),tnur(ncm,ncm),
     1 tnui(ncm,ncm),a(ncm,ncm),nchf(nlvm),ichi(nclm),ichf(nclm,nlvfm),
     2 ilv(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 nf=1,nlvfm
        do i=1,nlm
          lfmn(i,nf)=1000
          lfmx(i,nf)=0
        end do
      end do
c
c     loop over the partial waves
c
      jtest=0
      iwu=0
      if(imrd.eq.0) then
        open(10,file='KMAT.DAT')
        iread=10
      else
        open(11,file='tmatnu')
        iread=11
      end if
cold    7 read(iread,*,end=110) nchop,istot,j2,eksq
c  200 format(3i5,1pe14.8)                   !independent of stgf format
    7 read(iread,*,end=110) nchop,j2,jptt,eksq
      if(j2.lt.0)stop 'Error: input file contains LS-coupling data'
      jtest=jtest+j2
      if(nchop.gt.ncm) then
        write(6,205) nchop
  205   format('increase ncm to at least ',i5)
        stop 'increase ncm'
      end if
      fjtt=j2
      fjtt=fjtt/two
      lcmx=0
      do ic=1,nchop
        read(iread,210) ilv(ic),lc(ic),kv
  210   format(3i5)
        if(lc(ic).gt.lcmx) lcmx=lc(ic)
        fkv(ic)=kv
        fkv(ic)=fkv(ic)/two
      end do
      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 r-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 j2 for this matrix exceeds j2mx
c
      if(j2.gt.j2mx) go to 7
c
c     next determine whether this matrix contains the levels of
c     interest
c
      k=0
      do 35 ic=1,nchop
      if(ilv(ic).ne.ilvi) go to 35
      k=k+1
      ichi(k)=ic
   35 continue
      if(k.eq.0) go to 7
      nchi=k
      if(nchi.gt.nclm) then
        write(6,235) nchi
  235   format('increase nclm to at least ',i3)
        stop 'increase nclm'
      end if
      ifo=0
      nchfmx=0
      do 45 nf=1,nlvf
      k=0
      do 40 ic=1,nchop
      if(ilv(ic).ne.ilvf(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.nclm) then
        write(6,235) nchfmx
	stop 'increase nclm'
      end if
c
c     store the necessary information
c
      iwu=iwu+1
      if(iwu.gt.nwm)then
        stop 'increase nwm'
      end if
      fjtot(iwu)=fjtt
      ncci(iwu)=nchi
      do i=1,nchi
        icci(iwu,i)=ichi(i)
      end do
      do nf=1,nlvf
        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
        flcs(iwu,ic)=lc(ic)
        fkvs(iwu,ic)=fkv(ic)
      end do
c
c     print the k-matrix
c
      if(iprint.eq.1.and.imrd.eq.0) then
        write(6,300) fjtt,jptt
  300   format(//48x,'r-matrix for j =',f5.1,'  parity =',i3,
     1   /48x,34('*'))
        do i=1,nchop
          write(6,310) i,ilv(i),lc(i),fkv(i),(r(i,ip),ip=1,nchop)
  310     format(1x,i2,2i3,f6.1,1pe12.4,7e13.4/(15x,e12.4,7e13.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) fjtt,jptt
  320   format(//38x,'t-matrix for j =',f5.1,'  parity =',i3,
     1   /38x,34('*'))
        do i=1,nchop
	  flc=lc(i)
          write(6,330) i,ilv(i),fjv(i),flc,fkv(i),(tr(iwu,i,ip),
     1    ti(iwu,i,ip),ip=1,nchop)
  330     format(1x,i2,i3,3f6.1,3(1pe13.4,' +i*',e12.4)/
     1     24x,3(e13.4,' +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,nlvf
      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 if(jtest.eq.0)stop 'Error: old-style stgf header found'
      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(nlvfm=10,ncm=620,nclm=20,nlm=46,njm=2*nlm,
     1 nlbm=2*nlm,nwm=102,nftm=800,ntem=181)
      common/blok1/thmn,thmx,einc,theta(ntem),fji,fjf(nlvfm),
     1 sigbri(nlm),sigbrf(nlm,nlvfm),nthta,iprint,icharg,j2mx,nlvf,
     2 ilvi,ilvf(nlvfm),nthmx,limn,limx,lfmx(nlm,nlvfm),
     3 lfmn(nlm,nlvfm),ielas,iphr,nipc,nfpc
      common/blok2/gamma(nftm)
      common/blok3/fjtot(nwm),tr(nwm,ncm,ncm),ti(nwm,ncm,ncm),
     1 flcs(nwm,ncm),fkvs(nwm,ncm),ncci(nwm),nccf(nwm,nlvfm),
     2 icci(nwm,nclm),iccf(nwm,nclm,nlvfm),iwu,imrd
      common/blok4/fmr(nlm,nlvfm,nlm,njm,4),fmi(nlm,nlvfm,nlm,njm,4),
     1 qr(nlm),qi(nlm),alm(nlvfm,nlbm),dsg(ntem,nlvfm),sig(nlvfm),
     2 sigic(nlvfm)
c
      data zero,half,one,threehf,two/0.0d0,0.5d0,1.0d0,1.5d0,2.0d0/
c
c     initialize the arrays
c
      do ii=1,nlm
        qr(ii)=0.0
        qi(ii)=0.0
      end do
      do ikk=1,4
        do ij=1,njm
          do if=1,nlm
            do il=1,nlvfm
	      do ii=1,nlm 
                fmr(ii,il,if,ij,ikk)=zero
                fmi(ii,il,if,ij,ikk)=zero
              end do
	    end do
	  end do
	end do
      end do
c
c     loop over partial waves
c
      do 130 iw=1,iwu
      fjtt=fjtot(iw)
      fmjtt=two*fjtt+one
      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 nci=1,nchi
          ic=icci(iw,nci)
          li=flcs(iw,ic)
          ii=li+1
          qr(ii)=qr(ii)+fmjtt*(cos(2.*sigbri(ii))*tr(iw,ic,ic)-
     1     sin(2.*sigbri(ii))*ti(iw,ic,ic))
          qi(ii)=qi(ii)+fmjtt*(cos(2.*sigbri(ii))*ti(iw,ic,ic)+
     1     sin(2.*sigbri(ii))*tr(iw,ic,ic))
        end do
      end if
c
c     calculate the contributions to the real and imaginary parts of
c     array m
c
      fkki=fji-threehf
      if(fji.eq.zero) fkki=-half
      ikimx=2
      if(fji.eq.zero) ikimx=1
      do 125 iki=1,ikimx
      fkki=fkki+one
      fmkki=two*fkki+one
      do 120 nci=1,nchi
      ici=icci(iw,nci)
      fli=flcs(iw,ici)
      ii=fli+one
      fmli=two*fli+one
      fki=fkvs(iw,ici)
      fmki=two*fki+one
      coe6j1=s6j(fji,half,fkki,fjtt,fli,fki)
      do 115 nf=1,nlvf
      if(nccf(iw,nf).eq.0) go to 115
      ncfmx=nccf(iw,nf)
      j=ilvf(nf)
      fkkf=fjf(j)-threehf
      if(fjf(j).eq.0.0) fkkf=-half
      ikfmx=2
      if(fjf(j).eq.zero) ikfmx=1
      do 110 ikf=1,ikfmx
      fkkf=fkkf+one
      fmkkf=two*fkkf+one
      ikk=iki+2*ikf-2
      d1=abs(fkkf-fkki)
      s1=fkkf+fkki
      do 105 ncf=1,ncfmx
      icf=iccf(iw,ncf,nf)
      flf=flcs(iw,icf)
      lf=flf
      if=lf-lfmn(ii,nf)+1
      fmlf=two*flf+one
      d2=abs(flf-fli)
      fkf=fkvs(iw,icf)
      fmkf=two*fkf+one
      coe6j2=s6j(fjf(j),half,fkkf,fjtt,flf,fkf)
      iph=fkkf+fki+fkf+fjtt
      coeff=(-1)**iph*fmjtt*sqrt(fmli*fmlf*fmki*fmkf*fmkki*fmkkf)
      fjtmn=d1
      if(d2.gt.d1) fjtmn=d2
      s2=flf+fli
      fjtmx=s1
      if(s2.lt.s1) fjtmx=s2
      inj=0
      fjt=fjtmn-one
  100 fjt=fjt+one
      inj=inj+1
      coe6j3=s6j(fli,flf,fjt,fkkf,fkki,fjtt)
      coe6j=coe6j1*coe6j2*coe6j3
      fmr(ii,nf,if,inj,ikk)=fmr(ii,nf,if,inj,ikk)+coeff*coe6j*
     1 tr(iw,ici,icf)
      fmi(ii,nf,if,inj,ikk)=fmi(ii,nf,if,inj,ikk)+coeff*coe6j*
     1 ti(iw,ici,icf)      
      if(fjt.lt.fjtmx) go to 100
  105 continue
  110 continue
  115 continue
  120 continue
  125 continue
  130 continue
c     print the real and imaginary parts of the m array
      if (iprint.eq.4) then
        write(6,200)
  200   format(/)
        do ikk=1,4
          do inj=1,njm
            do if=1,nlm
              do nf=1,nlvfm
	        do ii=1,nlm
		  if(fmr(ii,nf,if,inj,ikk).ne.zero.or.
     1              fmi(ii,nf,if,inj,ikk).ne.zero) then
                    write(6,210) ii,nf,if,inj,ikk,fmr(ii,nf,if,inj,ikk),
     1                           fmi(ii,nf,if,inj,ikk)
  210               format(2x,'ii =',i2,2x,'nf =',i2,2x,'if =',i2,2x,
     1                     'inj =',i2,2x,'ikk =',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(nlvm=120,nlvfm=10,ncm=620,nclm=20,nlm=46,njm=2*nlm,
     1 nlbm=2*nlm,nwm=102,nftm=800,ntem=181)
      common/blok1/thmn,thmx,einc,theta(ntem),fji,fjf(nlvfm),
     1 sigbri(nlm),sigbrf(nlm,nlvfm),nthta,iprint,icharg,j2mx,nlvf,
     2 ilvi,ilvf(nlvfm),nthmx,limn,limx,lfmx(nlm,nlvfm),
     3 lfmn(nlm,nlvfm),ielas,iphr,nipc,nfpc
      common/blok2/gamma(nftm)
      common/blok3/fjtot(nwm),tr(nwm,ncm,ncm),ti(nwm,ncm,ncm),
     1 flcs(nwm,ncm),fkvs(nwm,ncm),ncci(nwm),nccf(nwm,nlvfm),
     2 icci(nwm,nclm),iccf(nwm,nclm,nlvfm),iwu,imrd
      common/blok4/fmr(nlm,nlvfm,nlm,njm,4),fmi(nlm,nlvfm,nlm,njm,4),
     1 qr(nlm),qi(nlm),alm(nlvfm,nlbm),dsg(ntem,nlvfm),sig(nlvfm),
     2 sigic(nlvfm)
      common/blok5/pl(nlbm),imbmx
      dimension alim(nlvm,nlbm)
c
      data zero,half,one,threehf,two/0.0d0,0.5d0,1.0d0,1.50d0,2.0d0/
c
      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 levels
c
      do 200 nf=1,nlvf
      j=ilvf(nf)
c
c     begin loop over lambda
c
      do 190 imb=1,imbmx
      lmb=imb-1
      flmb=lmb
      fmlmb=two*flmb+one
      alim(nf,imb)=zero
      alm(nf,imb)=zero
      do 180 ii=iimn,iimx
      li=ii-1
      fli=li
      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
      flip=lip
      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
      flf=lf
      inf=lf-lfmn(ii,nf)+1
      id1=iabs(lf-li)
      is1=lf+li
      ifpmn=iabs(lmb-lf)+1
      ifpmn=max0(ifpmn1,ifpmn)
      ifpmx=min0(ifpmx1,(lmb+lf+1))
      do 150 ifp=ifpmn,ifpmx
      lfp=ifp-1
      flfp=lfp
      if(mod(lfp+lf+lmb,2).ne.0) go to 150
      infp=lfp-lfmn(iip,nf)+1
      id1p=iabs(lfp-lip)
      is1p=lfp+lip
c
c     if iphr equal to 1, do not include the phase factor
c     (-1)**((li+lfp-lip-lf)/2) in order to correct for the fano phase
c     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
      fkki=fji-threehf
      if(fji.eq.zero) fkki=-half
      ikimx=2
      if(fji.eq.zero) ikimx=1
      do 140 iki=1,ikimx
      fkki=fkki+one
      fkkf=fjf(j)-threehf
      if(fjf(j).eq.zero) fkkf=-half
      ikfmx=2
      if(fjf(j).eq.zero) ikfmx=1
      do 130 ikf=1,ikfmx
      fkkf=fkkf+one
      ikk=iki+2*ikf-2
      id2=abs(fkkf-fkki)
      is2=fkkf+fkki
      id3=max0(id2,id1)
      is3=min0(is2,is1)
      id3p=max0(id2,id1p)
      ijtmn=max0(id3,id3p)+1
      ijtmx=min0(is3,is1p)+1
      if(ijtmn.gt.ijtmx) go to 130
      do 120 ijt=ijtmn,ijtmx
      jt=ijt-1
      fjt=jt
      inj=ijt-id3
      injp=ijt-id3p
      fmjt=two*fjt+one
      coef3=(-1)**(jt+lmb)*fmjt*fmlmb
      fr=fmr(ii,nf,inf,inj,ikk)
      frp=fmr(iip,nf,infp,injp,ikk)
      fi=fmi(ii,nf,inf,inj,ikk)
      fip=fmi(iip,nf,infp,injp,ikk)
      f1a=fr*frp
      f1b=fi*fip
      f2a=fr*fip
      f2b=frp*fi
      if(f1a.eq.zero.and.f1b.eq.zero.and.f2a.eq.zero.and.f2b.eq.zero)
     1 go to 120
      coef3=coef3*s6j(fli,flf,fjt,flfp,flip,flmb)
      sum1=sum1+coef3*(f1a+f1b)
      if(icharg.ne.0) sum2=sum2+coef3*(f2a-f2b)
  120 continue
  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,nlvf
        write(6,330) nf,lmbmx
  330   format(//31x,'a-lambda coefficient for level 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) nf,lmbmx
  350   format(//31x,'ai-lambda coefficient for level 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 total cross sections
c
      implicit real*8(a-h,o-z)
c
      parameter(nlvfm=10,ncm=620,nclm=20,nlm=46,njm=2*nlm,
     1 nlbm=2*nlm,nwm=102,nftm=800,ntem=181)
      common/blok1/thmn,thmx,einc,theta(ntem),fji,fjf(nlvfm),
     1 sigbri(nlm),sigbrf(nlm,nlvfm),nthta,iprint,icharg,j2mx,nlvf,
     2 ilvi,ilvf(nlvfm),nthmx,limn,limx,lfmx(nlm,nlvfm),
     3 lfmn(nlm,nlvfm),ielas,iphr,nipc,nfpc
      common/blok2/gamma(nftm)
      common/blok3/fjtot(nwm),tr(nwm,ncm,ncm),ti(nwm,ncm,ncm),
     1 flcs(nwm,ncm),fkvs(nwm,ncm),ncci(nwm),nccf(nwm,nlvfm),
     2 icci(nwm,nclm),iccf(nwm,nclm,nlvfm),iwu,imrd
      common/blok4/fmr(nlm,nlvfm,nlm,njm,4),fmi(nlm,nlvfm,nlm,njm,4),
     1 qr(nlm),qi(nlm),alm(nlvfm,nlbm),dsg(ntem,nlvfm),sig(nlvfm),
     2 sigic(nlvfm)
      common/blok5/pl(nlbm),imbmx
c
      data zero,one,two,three,four,eight,fourty5/0.0d0,1.0d0,
     1 2.0d0,3.0d0,4.0d0,8.0d0,4.5d1/
c
c     first calculate the differential cross sections as a function
c     of scattering angle
c
      fmji=two*fji+one
      pi=acos(-one)
      coeff=pi/(two*fmji*einc)
      coef3=one/(eight*fmji*einc)
      cnvrd=atan(one)/fourty5
      if(icharg.ne.0) then
        alpha=float(icharg)/sqrt(einc)
        coef1=alpha**2/(four*einc)
        coef2=-alpha/(four*fmji*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 levels
c
      do 30 nf=1,nlvf
      dsg1=zero
      dsg2=zero
      if(ilvf(nf).eq.ilvi.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
      do 50 nf=1,nlvf
      if(icharg.ne.0.and.ilvf(nf).eq.ilvi) go to 50
      sig(nf)=coeff*alm(nf,1)
c
c     calculate the total cross section by direct integration of the
c     differential cross section as an internal check
c
      h=cnvrd*theta(2)
      angi=cnvrd*theta(1)
      angf=cnvrd*theta(nthmx)
      sum=sin(angi)*dsg(1,nf)+sin(angf)*dsg(nthmx,nf)
      fact=two
      imx=nthmx-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
      sigic(nf)=two*pi*sum*h/three
   50 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=46)
      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=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(fj1,fj2,fj3,fl1,fl2,fl3)
c
c         calculates 6-j symbol(fj1,fj2,fj3,fl1,fl2,fl3) using
c         a 'nested' algorithm
c
      implicit real*8(a-h,o-z)
      data zero,half,one,two/0.0d0,0.5d0,1.0d0,2.0d0/
c
      parameter(nmft=800)
      common /blok2/gamma(nmft)
      s6j=zero
      i1=fj1+fj2+fj3+two
      i2=fj1+fj2-fj3+one
      i3=fj1+fj3-fj2+one
      i4=fj3+fj2-fj1+one
      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=fj1+fl2+fl3+two
      i2=fj1+fl2-fl3+one
      i3=fj1+fl3-fl2+one
      i4=fl3+fl2-fj1+one
      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=fj2+fl1+fl3+two
      i2=fj2+fl1-fl3+one
      i3=fj2+fl3-fl1+one
      i4=fl3+fl1-fj2+one
      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=fl2+fl1+fj3+two
      i2=fl2+fl1-fj3+one
      i3=fl2+fj3-fl1+one
      i4=fj3+fl1-fl2+one
      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)*half
      k1=fj1+fj2+fj3
      k2=fj1+fl2+fl3
      k3=fj2+fl1+fl3
      k4=fl2+fl1+fj3
      k5=fj1+fj2+fl2+fl1
      k6=fj2+fl2+fj3+fl3
      k7=fj1+fl1+fj3+fl3
      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=1.-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=46,nlbm=2*nlm)
      common/blok5/pl(nlbm),imbmx
      x= cos(t)
      a=1.0
      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