*module cstes
      MODULE cstes
c Physical constants: hc=hbar*c; e2=e^2/4 pi epsilon_0;hm=hbar^2/2m_N; ci=i
*      real*8, parameter :: pi=3.141592653589793238d0
      real*8, parameter :: pi=acos(-1d0)
      real*8 :: hc,e2,hm
      complex*16, parameter :: ci=(0,1)
      END MODULE

*module reseaufo
      MODULE reseaufo
c Wave function meshes
      integer :: Nr,Nt,Np
      real*8, parameter :: eps=1e-8
      real*8 :: rm,hr
      END MODULE

*module reseaueik
      MODULE reseaueik
c Meshes for eikonal phase calculation
      integer :: Nb,NZ,Ns
      real*8 :: bmax,hbpol,v0,hZ,hs
      real*8, allocatable, dimension (:) :: xb
      END MODULE

*module collision
      MODULE collision
c Data about the collision: spin of the fragment, mass and charge of the
c core(c), fragment(f), and target (T), reduced c-f mass,...
      integer :: iS,iSc,kZf,kZc,kZT
      real*8 :: e2ct,e2ft,Ac,Af,Apr,cc,cf,qc,qf,ATa,deuxm
      real*8 :: EAP,EKPT,etaPT
      END MODULE

*module potentials
      MODULE potentials
c Parameters of the two-body potentials: Vcf,VcT,VfT
      integer  :: ntypo,Njpo,npotN,lamax,lamin,jelj
      integer, allocatable, dimension (:) :: Jpo,lpo
      real*8, allocatable, dimension (:) :: Vp,rp,ap,Vls,rls,als,rC
      real*8, allocatable, dimension (:,:) :: Vcf
      real*8 :: VcT,WcT,WDcT,rRcT,rIcT,rDcT,aRcT,aIcT,aDcT,rCcT
      real*8 :: VfT,WfT,WDfT,rRfT,rIfT,rDfT,aRfT,aIfT,aDfT,rCfT
      END MODULE

*module boundstates
      MODULE boundstates
c Physical bound states of the projectile at the mesh points
      integer :: Nel
      integer, allocatable, dimension (:) :: NRel,Jel,lel
      real*8, allocatable, dimension (:) :: Eel
      real*8, allocatable, dimension (:,:) :: foel
      END MODULE

*module eli
      MODULE eli
c Bound states of the c-f potential to be removed by SuSy
      integer :: Neli
      integer, allocatable, dimension (:) :: NReli,Jeli,leli
      real*8, allocatable, dimension (:) :: Eeli
      END MODULE

*module phaseik
      MODULE phaseik
c Eikonal phase, and parameters for their calculation:
c lec determines the type of approximation
c     relcorr = 1 imply the use of gamma in the first order E1 phase
c     gaminv is the inverse of gamma => gaminv=sqrt(1-(v/c)^2)
      integer :: lec,Nl,relcorr
      real*8 :: gaminv
      complex*16, allocatable, dimension (:,:,:,:):: feik
      END MODULE

*module Sbu
      MODULE Sbu
c Breakup amplitude and angular amplitude
c Dimension of these amplitudes
      integer :: NlJM,Nsed,Nq,NE
      real*8, allocatable, dimension (:) :: theta,Enrj
      complex*16, allocatable, dimension(:,:,:,:) :: SklJM,fbu
      END MODULE

*pgm principal Chaconne
      program Chaconne
c     Pgm f90
c     Programme that computes reactions involving one-nucleon halo nuclei
c     (2-body projectile) within the Coulomb-Corrected Eikonal approximation.
c     It can compute the elastic and inelastic scattering cross sections
c     (angular distributions) and the breakup cross sections
c     (angular distribution over a range of continuum energies,
c     energy distribution, parallel-momentum distribution)
c     Coulomb breakup reactions including both
c     The optical potentials between projectile and target
c     include  both nuclear and Coulomb P-T interactions.
c     The Coulomb eikonal phase is renormalised as proposed by
c     Abu-Ibrahim and Suzuki [Prog. Theor. Phys. 112, 1013 (2004)],
c     following an idea of Margueron, Bonnaccorso, and Brink
c     [NPA 720, 337 (2003)].
c     The programme is based on a subroutine written by D. Baye 'multeik'.

      USE cstes
      USE reseaufo
      USE reseaueik
      USE collision
      USE potentials
      USE boundstates
      USE eli
      USE phaseik
      USE Sbu

      implicit real*8 (a-h,o-z)
      parameter (lnom=80,nbes=2,epsrel=5d-2)
      character (len=2) typOL,lnb,c2l,Rel
      character (lnom) fich,fichres,fichs0,fichsk,fichcT,fichfT
      real*8 :: fp(30),yp(30),Pl(30),Ylm(466)
      complex*16 rint,csljp,csljm,csljpE,csljmE,cgamma,cg,cres,ctmp,
     a     feik1,feik2
      integer, allocatable, dimension(:):: Nob
      real*8, allocatable, dimension(:):: bmin,hb
      real*8, allocatable, dimension(:):: fBessJ,fBessY
      real*8, allocatable, dimension(:):: sruth,fibb
      real*8, allocatable, dimension(:):: S,sdE,sedbdl,selast
      real*8, allocatable, dimension(:,:):: sinel
      real*8, allocatable, dimension(:):: Si02
      real*8, allocatable, dimension(:):: sedp
      real*8, allocatable, dimension (:,:,:):: foli
      real*8, allocatable, dimension(:,:):: ADlj
      real*8, allocatable, dimension(:,:,:,:):: ADfull
      complex*16, allocatable, dimension(:):: fcoul,ftild
      complex*16, allocatable, dimension(:,:):: phasOL
      complex*16, allocatable, dimension(:):: S20,Pstr
      complex*16, allocatable, dimension(:,:,:):: S0M,S0iM


c     Input reading
c Physical constants: hc=hbar*c; e2=e^2/4 pi epsilon_0;hm=hbar^2/2m_N
      read(*,*)hc,e2,hm
      amn=hc**2/(2*hm)
*      write(*,*)' hbar^2/2m_N',hm,' m_N = ',amn
c Mass, charge, and spin of core(c) and the fragment(f)
      read(*,*)Ac,Af,kZc,kZf,iSc,iS
      Apr=Af+Ac
      deuxm=(Af*Ac)/Apr/hm      !=2\mu/\hbar^2
*      write(*,*)' hbar^2/2mu',1/deuxm
c Spin of f and c (up to now Sc is set equal to 0)
      nAc=Ac+0.5d0
      nAf=Af+0.5d0
      lzc=2-mod(nAc,2)
      if(iS==0)then
c     fragment spin is neglected
         lzf=2
      else
         lzf=2-mod(nAf,2)
      endif

      iSc=iSc*lzc
      iS=iS*lzf

c c-f potential
      read(*,*)ntypo,Njpo
c ntypo=
c     1,11: different potentials for Njpo-1 partial waves +1 for all the others
c     2,12: different potentials for l-even and l-odd partial waves

c     1,2: Woods Saxon potentials
c     11,12: Gaussian potentials

      if(mod(ntypo,10)==2)then
         Njpo=2
      elseif(mod(ntypo,10)/=1)then
         write(*,*)'Unprogrammed potential'
         goto 9999
      endif
      allocate(Jpo(Njpo),lpo(Njpo))
      allocate(Vp(Njpo),rp(Njpo),ap(Njpo),rC(Njpo))
      allocate(Vls(Njpo),rls(Njpo),als(Njpo))

      If(mod(ntypo,10)==1)then
         do ipo=1,Njpo-1
            read(*,*)Jpo(ipo),lpo(ipo),Vp(ipo),rp(ipo),ap(ipo),
     a           Vls(ipo),rls(ipo),als(ipo),rC(ipo)
         enddo
         Jpo=Jpo*lzf
         read(*,*)Vp(Njpo),rp(Njpo),ap(Njpo),
     a        Vls(Njpo),rls(Njpo),als(Njpo),rC(Njpo)
      ElseIf(mod(ntypo,10)==2)then
            read(*,*)Vp(1),rp(1),ap(1),Vls(1),rls(1),als(1),rC(1)!1=odd l
            read(*,*)Vp(2),rp(2),ap(2),Vls(2),rls(2),als(2),rC(2)!2=even l
      EndIf

c Bound states (etats lie: el)
      read(*,*)Nel
      allocate(NRel(Nel),Jel(Nel),lel(Nel),Eel(Nel))
      read(*,*)(NRel(iel),Jel(iel),lel(iel),Eel(iel),iel=1,Nel)
      iJ0=Jel(1)
      l0=lel(1)
      Nef=iJ0+1
      Nsed=(Nef+1)/2
      Nelm=Nef
      do iel=2,Nel
         Nelm=Nelm+Jel(iel)+1
      enddo

**********************************************************************
c Forbidden states to be removed by SuSy
*      read(*,*)Neli
      Neli=0
*     allocate(NReli(Neli),Jeli(Neli),leli(Neli),Eeli(Neli))
*      read(*,*)(NReli(ieli),Jeli(ieli),leli(ieli),Eeli(ieli),
*     a     ieli=1,Neli)
*      Jeli=Jeli*lzf
**********************************************************************

c Collision data (target mass and charge, P energy/nucleon, trajectory type)
      read(*,*)ATa,kZT,EAP !,Rel

c Rel=
c     RV: relativistic velocity
c     NR: non-relativistic velocity
      Rel='RV'
      select case (Rel)
      case('RV','Rv','rV','rv') !relativistic velocity
         Rel='RV'
         v0=1+2*EAP*hm/hc**2
         v0=hc*sqrt(1-1/v0**2) !relativistic relative velocity*hbar
      case('NR','Nr','nR','nr')
         Rel='NR'
         v0=2*sqrt(EAP*hm) !non-relativistic velocity*hbar
      end select

      EKPT=abs(Apr*ATa/(Apr+ATa)/2/hm*v0) !incident wave vector

      e2ct=e2*kZc*kZT
      e2ft=e2*kZf*kZT

      cc=-Af/Apr
      cf=Ac/Apr

      etaPT=(e2ct+e2ft)/abs(v0)

      qc=2*e2ct*cc/v0
      qf=2*e2ft*cf/v0
c Relativistic corrections  (only if realcorr = 1)    
*      read(*,*) relcorr
      relcorr=0
      if(relcorr==1)then
         gaminv = 1.0/(1+2*EAP*hm/hc**2) !sqrt(1.0-(v0/hc)**2) !
      else
         gaminv = 1.0
      endif
c Projectile-target potentials: npotN=V_PT type; bmaxN=max b with nuclear V_PT
      read(*,*)npotN,bmaxN,lamax
c npotN=
c     0: purely Coulomb (point-point)
c     10: Woods-Saxon optical potential with point-spere Coulomb
c     20: Gaussian potential (no Coulomb)
c     100: chi phaseshift function given in a file
      If(lamax>=0)then        !sum multipoles from 0 up to lamax
         lamax=lamax
         lamin=0
      Else                      !only lamax term
         lamax=-lamax
         lamin=lamax
      Endif
      If(npotN<10)then !no potential to be read
         read(*,*)
         read(*,*)
         VcT=0
         WcT=0
         WDcT=0
         VfT=0
         WfT=0
         WDfT=0
      ElseIf(npotN==10)then !optical potentials
c Potentiel c-T
         read(*,*)VcT,WcT,WDcT,rRcT,rIcT,rDcT,aRcT,aIcT,aDcT,rCcT
c Potentiel f-T
         read(*,*)VfT,WfT,WDfT,rRfT,rIfT,rDfT,aRfT,aIfT,aDfT,rCfT
      ElseIf(npotN<30)then !Gaussian potentials
         WcT=0
         WDcT=0
         VfT=0
         WfT=0
         WDfT=0
c Potentiel c-T
         read(*,*)VcT,aRcT
c Potentiel f-T
         read(*,*)VfT,aRfT
      ElseIf(npotN==100)then
         read(*,*)fichcT
         fichcT=trim(fichcT)
         read(*,*)fichfT
         fichfT=trim(fichfT)
         open(unit=101,status='unknown',file=fichcT)
         open(unit=102,status='unknown',file=fichfT)
      EndIf
c Discretisation meshes
c Radial and angular meshes:
      read(*,*)Nr,rm,Nt,Np
      if(mod(Nt,2)/=0)then
         write(*,*)'Nt must be even: Nt=',Nt
         STOP
      endif
      if(mod(Np,2)/=0)then
         write(*,*)'Np must be even: Np=',Np
         STOP
      endif
      hr=rm/Nr
c Z and d meshes
      read(*,*)NZ,hZ,Ns,hs

c b mesh (impact parameter)
c Discretisation meshes b dependent: Nzb=Nb of b zones, bmax= max b
      read(*,*)Nzb,bmax,hbpol
      allocate(bmin(Nzb),hb(Nzb),Nob(Nzb))

      Nb=0      !total Nb of b
      izb=1
      DO
         read(*,*)bmin(izb),hb(izb)
         IF(bmin(izb)>=bmax)then
            if(izb==1)then
               Nob(izb)=1
               hb=1
            else
               Nob(izb)=0
               hb(izb)=0
            endif
            do izb2=izb+1,Nzb !skip next (useless) lines
               read(*,*)bmin(izb2),hb(izb2)
               Nob(izb2)=0
               hb(izb2)=0
            enddo
            Nzb=max(1,izb-1)
            EXIT
         ELSE
            If(izb>1)then
               Nob(izb-1)=(bmin(izb)-bmin(izb-1))/hb(izb-1)
               Nb=Nb+Nob(izb-1)
               if(Nob(izb-1)<0)then
                  write(*,'(''Invalid bmin in b zone'',i3)')izb
                  GOTO 9999
               elseif(abs(bmin(izb-1)+Nob(izb-1)*hb(izb-1)
     a                 -bmin(izb))>1d-5)then
                  write(*,'(''Invalid hb in b zone'',i3)')izb-1
                  GOTO 9999
               endif
            EndIf
         ENDIF
         if(izb==Nzb) EXIT
         izb=izb+1
      ENDDO
      If(bmin(Nzb)<bmax)then
         Nob(Nzb)=(bmax-bmin(Nzb))/hb(Nzb)+1
         Nb=Nb+Nob(Nzb)
         bmax=bmin(Nzb)+(Nob(Nzb)-1)*hb(Nzb)
      Else ! 1!b
         Nb=1
         bmax=bmin(1)
      EndIf

      allocate(xb(Nb))

      ib=0
      Do izb=1,Nzb
         do izb2=0,Nob(izb)-1
            ib=ib+1
            xb(ib)=bmin(izb)+izb2*hb(izb)
         enddo
      EndDo

c energies for outputs
c     for sdE: E=E0+(iE-1)*hE
c     for sdp: p=+/-(E0+ip*hE)
      read(*,*)NE,E0,hE
c angles for the angular-distribution output
      read(*,*)Nq,th0,hth
c Output selection variable
      read(*,*)lec
c lec=
c     0: sdb (breakup probability as function of b; only for the 1st energy)
c     1: sdE (energy distribution)
c     2: sdp (parallel-momentum distribution)
c     3: sdt (angular distribution integrated over an energy range)
c     7: str (one-neutron stripping cross section)
c     8: but (total breakup cross section)
c     9: sel (elastic scattering)(not programmed: issue with Hankel functions)
c     <10: usual eikonal
c     <20: modified eikonal
c     <30: modified eikonal with corection for scT>str
c     100<=lec<110: only first-order Coulomb FF approximation
      if(lec>=100.and.lec<110)npotN=0
      ilec=mod(lec,10)

c Name of output files
c     wave function in 'fichfo' & additional calculations in 'fichres'
      read(*,'(a80)')fich
      if(ilec<7)then
         fichsk=trim(fich)//'.sk'
      else
         fichs0=trim(fich)//'.s0'
      endif
      select case (ilec)
      case(0)
         fichres=trim(fich)//'.sdb'
         NE=1
         hE=1
      case(1)
         fichres=trim(fich)//'.sdE'
      case(2)
         fichres=trim(fich)//'.sdp'
      case(3)
         fichres=trim(fich)//'.sdt'
      case(7)
         fichres=trim(fich)//'.str'
         NE=0
         hE=1
      case(8)
         fichres=trim(fich)//'.but'
         NE=0
         hE=1
      case(9)
         fichres=trim(fich)//'.sel'
         NE=0
         hE=1
      end select

c     Printing information about the calculation inputs
c     (value of the constants, potentials used, numerical parameters,...)
      write(*,7001)hc,e2,hm
 7001 format(' hc=',3pd22.14,' e2=',1pd22.14,' hm=',2pd22.14)
      write(*,7002)Ac,Af,deuxm
 7002 format('Ac=',f5.2,' Af=',f5.2,' deuxm=',1pd22.14)
      write(*,7003)kZc,kZf
 7003 format('Zc=',i2,' Zf=',i2)
      write(*,7004)iSc,iS
 7004 format('Sc=',i2,' Sf=',i2)

      if(ntypo<10)then
         write(*,*)' Woods-Saxon Potential'
         select case (ntypo)
         case(1)
            do ipo=1,Njpo-1
               write(*,7006)Jpo(ipo),lpo(ipo),
     b              Vp(ipo),rp(ipo),ap(ipo),
     a              Vls(ipo),rls(ipo),als(ipo),rC(ipo)
            enddo
 7006       format(' J=',i3,'l=',i3,' Vp=',f8.3,' rp=',f7.3,' ap=',f7.3,
     a           ' Vls=',f7.3,' rls=',f7.3,' als=',f7.3,' rC=',f7.3)
            write(*,*)'Others:'
            write(*,7005)Vp(Njpo),rp(Njpo),ap(Njpo),
     a           Vls(Njpo),rls(Njpo),als(Njpo),rC(Njpo)
 7005       format(' Vp=',f8.3,' rp=',f7.3,' ap=',f7.3,' Vls=',f7.3,
     a           ' rls=',f7.3,' als=',f7.3,' rC=',f7.3)
         case(2)
            write(*,7035)Vp(1),rp(1),ap(1),Vls(1),rls(1),als(1),rC(1)
 7035       format(' Odd-l: Vp=',f8.3,' rp=',f7.3,' ap=',f7.3,
     a           ' Vls=',f7.3,' rls=',f7.3,' als=',f7.3,' rC=',f7.3)
            write(*,7025)Vp(2),rp(2),ap(2),Vls(2),rls(2),als(2),rC(2)
 7025       format('Even-l: Vp=',f8.3,' rp=',f7.3,' ap=',f7.3,
     a           ' Vls=',f7.3,' rls=',f7.3,' als=',f7.3,' rC=',f7.3)
         end select
      elseif(ntypo<20)then
         write(*,*)' Gaussian Potential'
         select case (ntypo)
         case(11)
            do ipo=1,Njpo-1
               write(*,7206)Jpo(ipo),lpo(ipo),
     b              Vp(ipo),rp(ipo),ap(ipo),
     a              Vls(ipo),rls(ipo),als(ipo),rC(ipo)
            enddo
 7206       format(' J=',i3,'l=',i3,' Vp=',f8.3,' rp=',f7.3,' ap=',f7.3,
     a           ' Vpp=',f7.3,' rpp=',f7.3,' app=',f7.3,' rC=',f7.3)
            write(*,*)'Others:'
            write(*,7205)Vp(Njpo),rp(Njpo),ap(Njpo),
     a           Vls(Njpo),rls(Njpo),als(Njpo),rC(Njpo)
 7205       format(' Vp=',f8.3,' rp=',f7.3,' ap=',f7.3,' Vpp=',f7.3,
     a           ' rpp=',f7.3,' app=',f7.3,' rC=',f7.3)

         case(12)
            write(*,7235)Vp(1),rp(1),ap(1),Vls(1),rls(1),als(1),rC(1)
 7235       format(' Odd-l: Vp=',f8.3,' rp=',f7.3,' ap=',f7.3,
     a           ' Vpp=',f7.3,' rpp=',f7.3,' app=',f7.3,' rC=',f7.3)
            write(*,7225)Vp(2),rp(2),ap(2),Vls(2),rls(2),als(2),rC(2)
 7225       format('Even-l: Vp=',f8.3,' rp=',f7.3,' ap=',f7.3,
     a           ' Vpp=',f7.3,' rpp=',f7.3,' app=',f7.3,' rC=',f7.3)
         end select
      endif

      write(*,7007)Nel
 7007 format('Nel=',i2)
      write(*,7008)(iel,NRel(iel),Jel(iel),lel(iel),Eel(iel),iel=1,Nel)
 7008 format('state',i2,': NR=',i2,' J= ',i2,' l= ',i2,' E0= ',f8.4)
      write(*,7011)ATa,kZT,EAP,v0/hc,VCK
 7011 format('ATa=',f6.2,' ZT=',i2,' E/A=',f7.3,' v0/c=',f6.3,
     a     ' VCK=',1pd13.6)
      write(*,7042)npotN
*      if(npotN==10.or.npotN==11.or.npotN==12)then
      if(npotN==10)then
         write(*,7040)VcT,WcT,WDcT,rRcT,rIcT,rDcT,aRcT,aIcT,aDcT,rCcT
         write(*,7041)VfT,WfT,WDfT,rRfT,rIfT,rDfT,aRfT,aIfT,aDfT,rCfT
*         if(npotN==11) write(*,*)'No f-T Coulomb potential'
*         if(npotN==12) write(*,*)'No c-T Coulomb potential'
 7042 format('npotN= ',i2)
 7040 format('VcT=',f9.4,' WcT=',f9.4,' WDcT=',f9.4,' rRcT=',f7.4,
     a     ' rIcT=',f7.4,' rDcT=',f7.4,' aRcT=',f7.4,' aIcT=',f7.4,
     b     ' aDcT=',f7.4,' rCcT=',f7.4)
 7041 format('VfT=',f9.4,' WfT=',f9.4,' WDfT=',f9.4,' rRfT=',f7.4,
     a     ' rIfT=',f7.4,' rDfT=',f7.4,' aRfT=',f7.4,' aIfT=',f7.4,
     b     ' aDfT=',f7.4,' rCfT=',f7.4)
      elseif(npotN==20)then
         write(2,7044)VcT,aRcT
         write(2,7045)VfT,aRfT
 7044 format('VcT=',f9.4,' aRcT=',f7.4)
 7045 format('VfT=',f9.4,' aRfT=',f7.4)
      elseif(npotN==100)then
         write(*,*)' Numerical eikonal phases read externally'
      endif
      write(*,7012)Nr,rm
 7012 format('Nr=',i7,' rm=',f9.2)
      write(*,7212)lamax,Nt,Np
 7212 format('lamax=',i3,' Nt=',i3,' Np=',i3)
      write(*,7013)NZ,hZ,Ns,hs
 7013 format('NZ= ',i5,' hZ=',f6.3,' Ns= ',i5,' hs=',f6.3)
      write(*,7014)bmin(1),bmax,Nb
 7014 format('bmin=',f6.2,' bmax=',f6.2,' Nbtot=',i8)
      do izb=1,Nzb
         write(*,7015)izb,bmin(izb),Nob(izb),hb(izb)
      enddo
 7015 format('b zone ',i3,' bmin=',f6.2,' Nb=',i4,' hb=',f7.4)
      write(*,7016)lec
 7016 format('lec=',i3)

      Nl=lel(1)+lamax
      Nlj=(iS+1)*(Nl+1)
      NlJM=(iS+1)*(Nl+1)**2
      ALLOCATE(foel(Nr,Nel),Vcf(Nr,Nlj),S(Nr),fibb(Nb))

      fibb(1)=0
      ib=0
      Do izb=1,Nzb
         ib=ib+1
         If(mod(Nob(izb),2)==1)then   !Simpson
            fibb(ib)=fibb(ib)+xb(ib)*hb(izb)/300
            fibb(ib+1)=4*xb(ib+1)*hb(izb)/300
            do ib2=3,Nob(izb)-1,2
               ib=ib+2
               fibb(ib)=2*xb(ib)*hb(izb)/300
               fibb(ib+1)=4*xb(ib+1)*hb(izb)/300
            enddo
            ib=ib+1
            fibb(ib)=xb(ib)*hb(izb)/300
         Else                   !Trapezes
            fibb(ib)=fibb(ib)+xb(ib)*hb(izb)/200
            do ib2=2,Nob(izb)-1
               ib=ib+1
               fibb(ib)=xb(ib)*hb(izb)/100
            enddo
            ib=ib+1
            fibb(ib)=xb(ib)*hb(izb)/200
         EndIf
      EndDo

      SELECT CASE(ilec)
      CASE(0:3)                  !breakup cross sections
         ALLOCATE(Enrj(NE),phasOL(NE,Nlj),foli(Nr,Nlj,NE))
         ALLOCATE(SklJM(NlJM,NE,Nsed,Nb))
         if(ilec==3)then        !angular distribution
            ALLOCATE(fbu(Nq,Nsed,NlJM,NE))
            ALLOCATE(theta(Nq))
            Do iq=1,Nq
               theta(iq)=pi/180*(th0+(iq-1)*hth)
            EndDo
            allocate(ADlj(Nq,0:Nlj))
         endif

         Do iE=1,NE             !defining energies at which SklJM are computed
            If(ilec==2)then     !parallel-momentum distribution
               Enrj(iE)=(E0+(iE-1)*hE)**2/deuxm
            Else                !energy or angular distributions
               Enrj(iE)=E0+(iE-1)*hE
            EndIf
         EndDo
      CASE(7)                   !total one-neutron stripping cross section
         ALLOCATE(Pstr(Nb))
      CASE(8)                   !total breakup cross section
         ALLOCATE(S20(Nb),Si02(Nb))
      CASE(9)                   !elastic scattering cross section
         ALLOCATE(theta(Nq),sruth(Nq),fcoul(Nq),ftild(Nq))
         ALLOCATE(S0M(Nb,Nef,Nsed),selast(Nq))
         ALLOCATE(S0iM(Nb,Nelm,Nsed),sinel(Nel,Nq))
         Do iq=1,Nq
            theta(iq)=pi/180*(th0+(iq-1)*hth)
         EndDo
      END SELECT

c Computes the c-f potential and the bound-state, and continuum wave functions
      jlj=0
      Do l=0,Nl
         Do iJ=abs(2*l-iS),2*l+iS,2
            jlj=jlj+1
            do ir=1,Nr          !computes Vcf=the core-fragment potential
               r=ir*hr
               Vcf(ir,jlj)=vJ(r,l,iJ)*deuxm
            enddo

******************************************************
c elimination of states by SuSy
*            U=Vcf(:,jlj)

            Do ieli=1,Neli
               if(l==leli(ieli).AND.iJ==Jeli(ieli))then
                  ins=l+1
                  call super(Nr,hr,ins,Eeli(ieli),NReli(ieli),
     a                 Vcf(:,jlj),S,eps)
               endif
            EndDo  !ieli
******************************************************

            do iel=1,Nel        !computes bound states
               If(l==lel(iel).and.iJ==Jel(iel))then
                  S2=0
                  Ek=0
                  S=Vcf(:,jlj)
                  call num1l(Nr,hr,Ek,S2,S,foel(:,iel),NRel(iel),eps)
                  if(abs(1-Ek/deuxm/Eel(iel))>epsrel)then
                     write(*,8001)iel,Eel(iel),Ek/deuxm
                     goto 9998
                  else
                     Eel(iel)=Ek/deuxm
                     write(*,8002)iel,Eel(iel)
                  endif
               EndIf
            enddo

           do iE=1,NE           !computes continuum states
              Ek=sqrt(Enrj(iE)*deuxm)
              eta=kZc*kZf*e2*deuxm/2/Ek
c     sigma_l
              If(eta/=0)then
                 if(l==0)then
                    cg=1+ci*eta
                    cg=cgamma(cg)
                    phasOL(iE,:)=atan2(aimag(cg),real(cg))
                 elseif(iJ==abs(2*l-iS))then
                    phasOL(iE,jlj:NlJ)=phasOL(iE,jlj:NlJ)+atan(eta/l)
                 endif
              Else
                 phasOL(iE,jlj)=0
              EndIf
c     Calcul de l'onde distordue sur le reseau uniforme. 'dephase' normalise 
c     la fct d'onde a cos(dl)*Fl(eta,k*r)+sin(dl)*Gl(eta,k*r) pour r->infini
              call dephase0(Nr,hr,Vcf(:,jlj),foli(:,jlj,iE),eps,dep,
     a              l,eta,Ek,rm)
c     Ajout du dephasage nucleaire a la phase et passage en exp.
              phasOL(iE,jlj)=phasOL(iE,jlj)+dep-mod(l,4)*pi/2
              phasOL(iE,jlj)=exp(ci*phasOL(iE,jlj))
           enddo !iE
         EndDo !iJ
      EndDo !l

 8001 Format('Error in binding energy of state',i2,2d16.8)
 8002 Format('Binding energy of state',i2,'=',1pd16.8,'MeV')
c Computes the eikonal phase with multipole expansion of the PT interaction
      call multeik

      if(npotN==100)then
         close(unit=101)
         close(unit=102)
      endif

      open(unit=2,status='unknown',file=fichres)

      SELECT CASE (ilec)
      CASE(0:3)
c Computes the Skljm
      SklJM=0

      cE=sqrt((1d0*iJ0+1d0)*(2d0*l0+1d0)/4/pi)
      DO 1 ib=1,Nb
         DO 2 ief=1,Nsed
            iM0=-iJ0+2*ief-2
            DO 3 iE=1,NE
               Egv=(Enrj(iE)-Eel(1))*gaminv/v0 !relativistic correction
               x=Egv*xb(ib)
               if(x==0)then
                  fp = 0
                  yp = 0
               else
                  call beskn(-x,nbes,fp,yp)
               endif
               If(lec<10)then   !Usual Eikonal => no correction
                  cor2=0
                  cor3=0
               Else             !Coulomb Corrected Eikonal
                  cor2=(qc+qf)*Egv*yp(2)
                  cor3=(qc+qf)*Egv*yp(1)
               EndIf

               jlJ=0
               jlJM=0
         Do l=0,Nl
            Do iJ=abs(2*l-iS),2*l+iS,2
               jlJ=jlJ+1
               Do iM=-iJ,iJ,2
                  jlJM=jlJM+1

                 iye=0
                 iyo=0
                 mu=(iM-iM0)/2  !fixed by Q3m
                 if(mu<0.and.mod(mu,2)==1)then
                    mfac=-1
                 else
                    mfac=1
                 endif
                 do la=0,lamax
                    if(abs(mu)>la)goto 421
                    call troisj(2*l,2*la,2*l0,0,0,0,Q3)
                    if(Q3==0)goto 421
                    call sixj(iJ,2*l,iS,2*l0,iJ0,2*la,Q6)
                    if(Q6==0)goto 421
                    call troisj(iJ0,2*la,iJ,iM0,2*mu,-iM,Q3m)
                    if(Q3m==0)goto 421

                    rint=0
                    If(mod(la+mu,2)==0)then !only feik 1 and 2
                       jye=iye+abs(mu)/2+1 !verif OK
                       do jr=1,Nr !radial integral
                          rint=rint+foel(jr,1)*foli(jr,jlj,iE)*
     a                         (feik(jye,jr,ib,1)
     b                         +ci*cor2*feik(jye,jr,ib,2))
                       enddo
                    Else        !only feik 3
                       jyo=iyo+abs(mu)/2+1
                       do jr=1,Nr !radial integral
                          rint=rint-foel(jr,1)*foli(jr,jlj,iE)*
     a                         cor3*feik(jyo,jr,ib,3)
                       enddo
                    EndIf
                    SklJM(jlJM,iE,ief,ib)=SklJM(jlJM,iE,ief,ib)
     a                   +rint*hr*Q3*Q6*Q3m*sqrt(2d0*la+1d0)
 421                continue
                    iye=iye+la/2+1
                    iyo=iyo+(la+1)/2
                 enddo !la
                  if(mod(iS-iM,4)==0)then
                     SklJM(jlJM,iE,ief,ib)=SklJM(jlJM,iE,ief,ib)
     a                    *mfac*cE*phasOL(iE,jlJ)
     b                    *sqrt((1d0*iJ+1d0)*(2d0*l+1d0))
                  else
                     SklJM(jlJM,iE,ief,ib)=-SklJM(jlJM,iE,ief,ib)
     a                    *mfac*cE*phasOL(iE,jlJ)
     b                    *sqrt((1d0*iJ+1d0)*(2d0*l+1d0))
                  endif
              EndDo !iM
           EndDo !iJ
        EndDo !il
 3    continue  !iE
 2    continue  !m0
 1    continue  !b

*ERSk
c     Printing of SklJM
*      open(unit=21,status='unknown',file=fichsk)
*      do ib=1,Nb
*         b=b0+(ib-1)*hb
*         write(21,*)b,(abs(SklJM(icomp,1,1,ib)),
*     a        phase(SklJM(icomp,1,1,ib)),icomp=1,NlJM)
*      enddo
*      close(21)

      Select Case(ilec)
      Case(0,1)
         ALLOCATE(sedbdl(0:Nlj),sdE(0:Nlj))
c     Energy distribution or Pbu(b) at 1 energy
         DO iE=1,NE
            Ek=sqrt(deuxm*Enrj(iE))
            sdE=0
            Do ib=1,Nb
               sedbdl=0
               do ief=1,Nsed
                  if(ief==Nsed.and.lzf==2)then !m0=0
c     factor for the summation over m0
                     m0f=1
                  else
                     m0f=2
                  endif
                  jlJM=0
                  jlJ=0
                  Do l=0,Nl
                     Do iJ=abs(2*l-iS),2*l+iS,2
                        jlj=jlj+1
                        slj=0
                        do iM=-iJ,iJ,2
                           jlJM=jlJM+1
                           sedbdl(jlj)=sedbdl(jlj)+
     a                          m0f*abs(SklJM(jlJM,iE,ief,ib))**2
                           sedbdl(0)=sedbdl(0)+
     a                          m0f*abs(SklJM(jlJM,iE,ief,ib))**2
                        enddo   !iM
                     EndDo      !iJ
                  EndDo         !l
               enddo            !m0
               sedbdl=sedbdl*deuxm/pi/Ek/Nef
*ERb
               If(ilec==0)then  !thus NE=1
                  write(2,3003)xb(ib),(sedbdl(il),il=0,Nlj)
               Else
c     Integration over b
                  sdE=sdE+fibb(ib)*sedbdl
               EndIf
            EndDo !ib
*ERE
            If(ilec==1)then
               sdE=2*pi*sdE
               write(2,3003)Enrj(iE),(sdE(il),il=0,Nlj)
            EndIf
         ENDDO !iE
         DEALLOCATE(sedbdl,sdE)

         Case(2)
c     Parallel momentum distribution in the projectile CM rest frame
            ALLOCATE(sedp(-NE:NE))
            sedp=0
      do ib=1,Nb
c     Integration over b
         do ief=1,Nsed
            if(ief==Nsed.and.lzf==2)then
c     factor for the summation over m0
               m0f=1
            else
               m0f=2
            endif
            DO ikp=0,NE
               if (ikp==0)then
                  Ekp=0
               else
                  Ekp=E0+(ikp-1)*hE
               endif
               DO ik=max(ikp,1),NE
                  Ek=E0+(ik-1)*hE
                  costk=Ekp/Ek
                  call HARSPH(costk,-Nl-1,Pl,Ylm)
                  smmsp=0
                  smmsm=0
                  Do iM=-2*Nl-iS,2*Nl+iS,2
                     Do ims=-iS,iS,2
                        iml=iM-ims
                        csljp=0
                        csljm=0
                        do il=abs(iml),2*Nl,2
                           l=il/2
                           do iJ=abs(il-iS),il+iS,2
                              If(iJ>=abs(iM))then
                                 jlJM=jfromlSJM(l,iS,iJ,iM)
                                 call clebs(il,iS,iJ,iml,ims,iM,CGO)
                                 jYlm=l*(l+1)
                                 jYlm=jYlm/2+abs(iml)/2+1
c     NB: HARSPH gives the Ylm for m>=0 only, and phi=0
c     values for m<0 are the same with a phase (-)^m
c     this phase is irrelevant here because of the |...|^2 in the sum over m
                                 csljp=csljp+CGO*Ylm(jYlm)
     a                                *SklJM(jlJM,ik,ief,ib)
                                 if(mod(il-iml,4)==0)then
                                    csljm=csljm+CGO*Ylm(jYlm)
     a                                   *SklJM(jlJM,ik,ief,ib)
                                 else
                                    csljm=csljm-CGO*Ylm(jYlm)
     a                                   *SklJM(jlJM,ik,ief,ib)
                                 endif
                              EndIf
                           enddo !iJ
                        enddo !l
                        smmsp=smmsp+abs(csljp)**2
                        smmsm=smmsm+abs(csljm)**2
                     EndDo !ims
                  EndDo !iM
c     Integral over k by Simpson's rule
c     (breakup probabilties above Nk are assumed to be negligible)
c     Sum over m0 and integral over b (fibb) are included
                  if(ikp==0)then
                     if(ik==1)then !first step by trapezes
                        sedp(0)=sedp(0)
     a                       +fibb(ib)*m0f*smmsp/Ek*(E0/2+hE/3)
                     elseif(mod(ik,2)==0)then !rest by Simpson
                        sedp(0)=sedp(0)
     a                       +fibb(ib)*m0f*4*smmsp*hE/Ek/3
                     else
                        sedp(0)=sedp(0)
     a                       +fibb(ib)*m0f*2*smmsp*hE/Ek/3
                     endif
                  elseif(ik==ikp)then
                     sedp(ikp)=sedp(ikp)
     a                    +fibb(ib)*m0f*smmsp*hE/Ek/3
                     sedp(-ikp)=sedp(-ikp)
     a                    +fibb(ib)*m0f*smmsm*hE/Ek/3
                  elseif(mod(ik-ikp,2)==1)then
                     sedp(ikp)=sedp(ikp)
     a                    +fibb(ib)*m0f*4*smmsp*hE/Ek/3
                     sedp(-ikp)=sedp(-ikp)
     a                    +fibb(ib)*m0f*4*smmsm*hE/Ek/3
                  else
                     sedp(ikp)=sedp(ikp)
     a                    +fibb(ib)*m0f*2*smmsp*hE/Ek/3
                     sedp(-ikp)=sedp(-ikp)
     a                    +fibb(ib)*m0f*2*smmsm*hE/Ek/3
                  endif
               ENDDO !ik
            ENDDO !ikp
         enddo !m0
      enddo !ib

      sedp=8*pi*sedp/Nef
*ERp
      Do ikp=-NE,NE
         if(ikp<0)then
            write(2,'(1f8.4,1pd18.8)')-E0+(ikp+1)*hE,sedp(ikp)
         elseif(ikp==0)then
            write(2,'(1f8.4,1pd18.8)')ikp*hE,sedp(ikp)
         else
            write(2,'(1f8.4,1pd18.8)')E0+(ikp-1)*hE,sedp(ikp)
         endif
      EndDo
      DEALLOCATE (sedp)
      Case(3)
c Angular distribution including integration over energy

         call calculFBU

         ADlj=0
         do iq=1,Nq
            Do ief=1,Nsed
               if(ief==Nsed.and.lzf==2)then !m0=0
c     factor for the summation over m0
                  m0f=1
               else
                  m0f=2
               endif
               icomp=0
               jlJ=0
               Do l=0,Nl
                  Do iJ=abs(2*l-iS),2*l+iS,2
                     jlj=jlj+1
                     do iM=-iJ,iJ,2
                        icomp=icomp+1
*                  if(NE==1)hE=1.0d0!si rectangles
                  if(NE==1)hE=2 !si trapezes
                  do iE=1,NE
                     cbu=deuxm/pi/sqrt(deuxm*Enrj(iE))*EKPT**2/100 ! in barn
!     Integration par trapezes
                     if(iE==1.or.iE==NE)then
                        ADlj(iq,jlj)=ADlj(iq,jlj)+m0f*hE/2*cbu
     a                       *(abs(fbu(iq,ief,icomp,iE)))**2
                        ADlj(iq,0)=ADlj(iq,0)+m0f*hE/2*cbu
     a                       *(abs(fbu(iq,ief,icomp,iE)))**2
                     else
                        ADlj(iq,jlj)=ADlj(iq,jlj)+m0f*hE*cbu
     a                       *(abs(fbu(iq,ief,icomp,iE)))**2
                        ADlj(iq,0)=ADlj(iq,0)+m0f*hE*cbu
     a                       *(abs(fbu(iq,ief,icomp,iE)))**2
                     endif
                  enddo !iE
                     enddo !iM
                  EndDo !iJ
               EndDo !l
            enddo !ief
            ADlj(iq,:)=ADlj(iq,:)/Nef
*ERt
*            write(2,'(1pd10.4,60pd18.8)')180/pi*theta(iq),AD(iq),
            write(2,3003)180/pi*theta(iq),(ADlj(iq,jlj),jlj=0,Nlj)
         enddo

      DEALLOCATE (fbu,theta,ADlj)
      DEALLOCATE (Enrj,phasOL,foli,SklJM)

      End Select                !ilec for breakup distributions

      CASE(7)
c Computes the total one-neutron stripping cross section
c ie total inelastic one-neutron removal cross section
      Pstr=0

      il0=2*lel(1)
      iJ0=Jel(1)
      cstr=(1d0*il0+1d0)/2/sqrt(pi)
      if(mod(iJ0,2)/=0)cstr=-cstr
      iye=1
      Do la=0,min(lamax,il0,iJ0),2 !only even la because of Q3
         ila=2*la
         call troisj(il0,ila,il0,0,0,0,Q3)
         if(Q3==0)goto 827
         call sixj(iJ0,il0,iS,il0,iJ0,ila,Q6)
         if(Q6==0)goto 827

         sum0=0
         do im0=-iJ0,iJ0,2
            call troisj(iJ0,ila,iJ0,iM0,0,-iM0,Q3j)
            if(mod(iS+iM0,2)==0)then
               sum0=sum0+Q3j
            else
               sum0=sum0-Q3j
            endif
         enddo
         DO ib=1,Nb
            rint=0
!only feik 1 because no correction E1
            do jr=1,Nr          !radial integral
               if(la==0)then
                  feik1=feik(iye,jr,ib,1)+sqrt(4*pi)
               else
                  feik1=feik(iye,jr,ib,1)
               endif
               rint=rint+(foel(jr,1))**2*feik1
            enddo !jr
            Pstr(ib)=Pstr(ib)+rint*hr*Q3*Q6*Q3j*cstr
         ENDDO   !b
 827     continue
         iye=iye+la+2
      EndDo                     !la

      open(unit=20,status='unknown',file=fichs0)
      do ib=1,Nb
         write(20,'(1f10.4,2(1pd18.8))')xb(ib),
     a        real(Pstr(ib)),aimag(Pstr(ib))
      enddo
      close(unit=20)
c Computes total one-neutron stripping cross section (integral over b)
      Str=0
      do ib=2,Nb-1
         Str=Str+fibb(ib)*real(Pstr(ib))
      enddo
      Str=2*pi*Str
      write(*,'(''Total stripping cross section='',1pd18.8)')Str
      deallocate(Pstr)
*ERbut
      write(2,'(1pd18.8)')Str

      CASE(8)
c Computes the total breakup cross section through closing relation.
c 1. Computes the norm of the wave function
      S20=0
      il0=2*l0
      CS2=sqrt((1d0*il0+1d0)/(1d0*iJ0+1d0))/2/pi
      if(mod(il0+iS+iJ0,4)/=0)Cs2=-Cs2
      iLAmax=2*min(il0,iJ0)
      DO 800 iLA=0,iLAmax,4     !LA=2*Lambda, and Lambda is even by Q3l
         if(lzf==2)then         !case iM0=0
            call clebs(iJ0,iLA,iJ0,0,0,0,Q3jm)
            Q3jm=Q3jm/2
         else                   !case iM0=1
            call clebs(iJ0,iLA,iJ0,1,0,1,Q3jm)
         endif
         do iM0=iJ0,2,-2
            call clebs(iJ0,iLA,iJ0,iM0,0,iM0,Q3)
            Q3jm=Q3jm+Q3
         enddo
         call clebs(il0,iLA,il0,0,0,0,Q3l)
         if(Q3l==0)goto 800
         call sixj(iJ0,il0,iS,il0,iJ0,iLA,Q6)
         if(Q6==0)goto 800
         iye1=0
         Do la1=0,lamax !lambda
            ila1=2*la1
            iye2=0
            Do la2=0,lamax !lambda prime
               ila2=2*la2
               call clebs(ila2,ila1,iLA,0,0,0,Q3lala)
               if(Q3lala==0)goto 802
               Q3lala=Q3lala*sqrt((1d0*ila1+1d0)*(1d0*ila1+1d0))
               do 803 mu=-la1,la1
                  if(abs(mu)>la2)goto 803
                  call clebs(ila2,ila1,iLA,-2*mu,2*mu,0,Q3lamu)
                  if(Q3lamu==0)goto 803
                  DO ib=1,Nb
                     rint=0
                     If(mod(la1+mu,2)==0.and.mod(la2+mu,2)==0)then
!only feik 1 because the E1 correction is neglected to compute norm of wf
                        jye1=iye1+abs(mu)/2+1
                        jye2=iye2+abs(mu)/2+1
                        do jr=1,Nr
                           if(la1==0.and.mu==0)then
                              feik1=feik(jye1,jr,ib,1)+sqrt(4*pi)
                           else
                              feik1=feik(jye1,jr,ib,1)
                           endif
                           if(la2==0.and.mu==0)then
                              feik2=conjg(feik(jye2,jr,ib,1)+sqrt(4*pi))
                           else
                              feik2=conjg(feik(jye2,jr,ib,1))
                           endif
                           rint=rint+foel(jr,1)**2*feik1*feik2
                        enddo
                        if(mod(mu,2)==0)then
                           S20(ib)=S20(ib)
     a                          +rint*hr*Q3lamu*Q3lala*Q6*Q3l*Q3jm
                        else
                           S20(ib)=S20(ib)
     a                          -rint*hr*Q3lamu*Q3lala*Q6*Q3l*Q3jm
                        endif
                     EndIf
                  ENDDO !ib
 803           continue         !mu
 802           Continue
               iye2=iye2+la2/2+1
            EndDo               !la2
            iye1=iye1+la1/2+1
         EndDo                  !la1
 800  CONTINUE                  !LA
      S20=S20*CS2

c 2. Computes the S0 for the P-T elastic scattering (simple modif of CASE(9))
      Si02=0

      DO iel=1,Nel
         il=2*lel(iel)
         iJ=Jel(iel)
         cSi0=(1d0*il0+1d0)*(1d0*il+1d0)*(1d0*iJ+1d0)/4/pi
         iye=0
         Do la=0,lamax
            ila=2*la
            if(ila<abs(il-il0).or.ila>il+il0)goto 820
            call troisj(il,ila,il0,0,0,0,Q3)
            if(Q3==0)goto 820
            call sixj(iJ,il,iS,il0,iJ0,ila,Q6)
            if(Q6==0)goto 820
            do mu=-la,la
               DO ib=1,Nb
                  rint=0
                  If(mod(la+mu,2)==0)then !only feik 1 because no correction E1
                     jye=iye+abs(mu)/2+1
                     do jr=1,Nr !radial integral
                        if(la==0.and.mu==0)then
                           feik1=feik(jye,jr,ib,1)+sqrt(4*pi)
                        else
                           feik1=feik(jye,jr,ib,1)
                        endif
                        rint=rint+foel(jr,iel)*foel(jr,1)*feik1
                     enddo
                  EndIf
                  Si02(ib)=Si02(ib)+(abs(rint*hr*Q3*Q6))**2*cSi0
               ENDDO !ib
            enddo
 820        continue
            iye=iye+la/2+1
         EndDo !la
      ENDDO !iel

      open(unit=20,status='unknown',file=fichs0)
      do ib=1,Nb
         write(20,'(1f10.4,3(1pd18.8))')xb(ib),real(S20(ib)),
     a        aimag(S20(ib)),Si02(ib)
      enddo
      close(unit=20)
c Computes total breakup cross section (integral over b)
      Sbut=0
      do ib=1,Nb
         Sbut=Sbut+fibb(ib)*(real(S20(ib))-Si02(ib))
      enddo
      Sbut=2*pi*Sbut
      write(*,*)"computing bu"
      write(*,'(''Total breakup cross section='',1pd18.8)')Sbut
      deallocate(S20,Si02)
*ERbut
      write(2,'(1pd18.8)')Sbut

      CASE (9)
c Computes the S0 for the P-T elastic scattering
      S0M=0

      cE=(1d0*iJ0+1d0)*(2d0*l0+1d0)/2/sqrt(pi)

      Do ief=1,Nsed
         iM0=-iJ0+2*ief-2
         do ief2=1,Nef
            iM=-iJ0+2*ief2-2
            mu=(iM-iM0)/2       !fixed by Q3m
            iye=0
            if(mu<0.and.mod(mu,2)==1)then
               mfac=-1
            else
               mfac=1
            endif
            do la=0,lamax
               if(abs(mu)>la)goto 429
               call troisj(2*l0,2*la,2*l0,0,0,0,Q3)
               if(Q3==0)goto 429
               call sixj(iJ0,2*l0,iS,2*l0,iJ0,2*la,Q6)
               if(Q6==0)goto 429
               call troisj(iJ0,2*la,iJ0,iM0,2*mu,-iM,Q3m)
               if(Q3m==0)goto 429

               DO ib=1,Nb
                  rint=0
                  If(mod(la+mu,2)==0)then !only feik 1 because no correction E1
                     jye=iye+abs(mu)/2+1
                     do jr=1,Nr !radial integral
c     No correction E1 at first order since elastic (l=l0, and J=J0)
                        rint=rint+foel(jr,1)**2*feik(jye,jr,ib,1)
                     enddo
                  EndIf
                  S0M(ib,ief2,ief)=S0M(ib,ief2,ief)
     a                 +rint*hr*Q3*Q6*Q3m*sqrt(2d0*la+1d0)
               ENDDO !b
 429           continue
               iye=iye+la/2+1
            enddo !la
            if(mod(iS-iM,4)==0)then
               S0M(:,ief2,ief)=S0M(:,ief2,ief)*mfac*cE
            else
               S0M(:,ief2,ief)=-S0M(:,ief2,ief)*mfac*cE
            endif
         enddo !ief2 (ie m)
      EndDo !ief (ie m0)

*ERS0
c     Printing of S0
      open(unit=20,status='unknown',file=fichs0)
      do ib=1,Nb
         write(20,'(1f10.4,8(1pd18.8))')xb(ib),(abs(S0M(ib,icomp,1)),
     a        phase(S0M(ib,icomp,1)),icomp=1,Nef)
      enddo !ib
      close(20)

c     Computes sigma_0 for PT motion
      cg=1+ci*etaPT
      cg=cgamma(cg)
      s0PT=atan2(aimag(cg),real(cg)) !OK

c Computes Rutherford cross section 'sruth', 
c and Coulomb elastic scattering amplitude 'fcoul'
      Do iq=1,Nq
         sruth(iq)=(etaPT/EKPT)**2/sin(theta(iq)/2)**4/400
         fcoul(iq)=-ci*etaPT/2.0d0/EKPT**2/sin(theta(iq)/2)**2
     a        *exp(-2*ci*etaPT*log(2*EKPT*sin(theta(iq)/2))
     b        +2*ci*S0PT)
      EndDo

c Computes elastic scattering cross section
      open(unit=20,status='unknown',file=fichs0)

      nla=20 !nb of points for the Gauss integration in 'intbesm'

      NJmax=iJ0
      Allocate(fBessJ(NJmax+1),fBessY(NJmax+1))
      Do ief=1,Nsed
         iM0=-iJ0+(ief-1)*2  !2m0
         Do ief2=1,Nef
            iM02=-iJ0+(ief2-1)*2 !2m0'
            ibess=abs(iM0-iM02)/2+1
            ftild=0
            DO ib=1,Nb-1
               Npol=nint((xb(ib+1)-xb(ib))/hbpol)
               if(abs(Npol*1.0d0-(xb(ib+1)-xb(ib))/hbpol)>1d-3)then
                  write(*,*)'Error in interpolation at ib=',ib
                  write(*,*)'hb=',xb(ib+1)-xb(ib)
                  write(*,*)'hbpol=',hbpol
                  STOP
               endif

               do ipol=1,Npol
                  bpol=xb(ib)+(ipol-1)*hbpol
                  if(bpol<1d-3)then
                     chic2=0
                  else
*                     chic2=2*etaPT*log(bpol/2)
c     The Coulomb scattering phase is chosen equal to the
c     Coulomb eikonal phase that reproduces the exact Coulomb amplitude
c     cf. Bertulani & Danielewicz pg 423 Eq.(8.77).
                  chic2=2*etaPT*log(EKPT*bpol)
                  endif
                  ctmp=S0M(ib,ief2,ief)+(S0M(ib+1,ief2,ief)
     a                 -S0M(ib,ief2,ief))*(ipol-1)/Npol
                  do iq=1,Nq
                     q=2*EKPT*sin(theta(iq)/2)
*                     z=q*bpol !version w/ S17*
                     x=q*bpol
                     If(ibess==1)then !iM02==iM0
*                        BessJ0=S17AEF(q*bpol,IFAIL)
*                        call BESSJ(NJmax,q*bpol,fBessJ,fBessY)
                        call bessj0y0(x,fBessJ(1),fBessY(1))
                        ftild(iq)=ftild(iq)+fBessJ(ibess)*ctmp
     a                       *bpol*exp(ci*chic2)*hbpol
                     ElseIf(ibess==2)then !iM02=iM0+-1
                        fBessJ(2)=bessj1(x)
                        ftild(iq)=ftild(iq)+fBessJ(ibess)*ctmp
     a                       *bpol*exp(ci*chic2)*hbpol
                     Else
*                        z=q*bpol
*                        call S17DEF(0.0d0,z,NJmax,'U',cBessJ,NZ,IFAIL)
                        call BESSJ(NJmax,x,fBessJ,fBessY)
*                        ftild(iq)=ftild(iq)+cBessJ(ibess+1)*ctmp
*     a                       *bpol*exp(ci*chic2)*hbpol
                        ftild(iq)=ftild(iq)+fBessJ(ibess)*ctmp
     a                       *bpol*exp(ci*chic2)*hbpol
                     Endif
                  enddo !iq
               enddo !bpol
            ENDDO !b

            Do iq=1,Nq
               q=2*EKPT*sin(theta(iq)/2)
               call intbesm(ibess-1,bmax,q,etaPT,nla,cres)
               ftild(iq)=ftild(iq)+S0M(Nb,ief2,ief)*bmax**2*cres
               If(iM0==iM02)then
                  selast(iq)=selast(iq)+2*abs(fcoul(iq)+ftild(iq))**2
               Else
                  selast(iq)=selast(iq)+2*abs(ftild(iq))**2
               Endif
            EndDo
         EndDo !ief2 (ie m)
      EndDo !ief(ie m0)

      Deallocate(fBessJ,fBessY)

      close(unit=20)

      do iq=1,Nq
         selast(iq)=EKPT**2*selast(iq)/Nef/100
      enddo

*********************************************
c Computes the S0iM
      S0iM=0

      cE=sqrt((1d0*iJ0+1d0)*(2d0*l0+1d0)/4/pi)
      DO ib=1,Nb
         DO ief=1,Nsed
            iM0=-iJ0+2*ief-2
            ielm=Nef
            DO iel=2,Nel
               Egv=(Eel(iel)-Eel(1))*gaminv/v0 !relativistic correction
               x=Egv*xb(ib) 
               if(x==0)then
                  fp = 0
                  yp = 0
               else
                  call beskn(-x,nbes,fp,yp)
               endif
               If(lec<10)then   !Usual Eikonal => no correction
                  cor2=0
                  cor3=0
               Else             !Coulomb Corrected Eikonal
                  cor2=(qc+qf)*Egv*yp(2)
                  cor3=(qc+qf)*Egv*yp(1)
               EndIf

               l=lel(iel)
               iJ=Jel(iel)

               Do iM=-iJ,iJ,2
                  ielm=ielm+1
                  iye=0
                  iyo=0
                  mu=(iM-iM0)/2 !fixed by Q3m
                  if(mu<0.and.mod(mu,2)==1)then
                     mfac=-1
                  else
                     mfac=1
                  endif
                  do la=0,lamax
                     if(abs(mu)>la)goto 422
                     call troisj(2*l,2*la,2*l0,0,0,0,Q3)
                     if(Q3==0)goto 422
                     call sixj(iJ,2*l,iS,2*l0,iJ0,2*la,Q6)
                     if(Q6==0)goto 422
                     call troisj(iJ0,2*la,iJ,iM0,2*mu,-iM,Q3m)
                     if(Q3m==0)goto 422

                     rint=0
                     If(mod(la+mu,2)==0)then !only feik 1 and 2
                        jye=iye+abs(mu)/2+1 !verif OK
                        do jr=1,Nr !radial integral
                           rint=rint+foel(jr,1)*foel(jr,iel)*
     a                          (feik(jye,jr,ib,1)
     b                          +ci*cor2*feik(jye,jr,ib,2))
                        enddo
                     Else       !only feik 3
                        jyo=iyo+abs(mu)/2+1 !verif OK
                        do jr=1,Nr !radial integral
                           rint=rint-foel(jr,1)*foel(jr,iel)*
     a                          cor3*feik(jyo,jr,ib,3)
                        enddo
                     EndIf
                     S0iM(ib,ielm,ief)=S0iM(ib,ielm,ief)
     a                   +rint*hr*Q3*Q6*Q3m*sqrt(2d0*la+1d0)
 422                 continue
                     iye=iye+la/2+1
                     iyo=iyo+(la+1)/2
                  enddo         !la
                  if(mod(iS-iM,4)==0)then
                     S0iM(ib,ielm,ief)=S0iM(ib,ielm,ief)
     a                    *mfac*cE
     b                    *sqrt((1d0*iJ+1d0)*(2d0*l+1d0))
                  else
                     S0iM(ib,ielm,ief)=-S0iM(ib,ielm,ief)
     a                    *mfac*cE
     b                    *sqrt((1d0*iJ+1d0)*(2d0*l+1d0))
                  endif
               EndDo !iM
            ENDDO !iel
         ENDDO !iM0
      ENDDO !ib

c Calcul de Sinel
         write(*,*)'Calcul de la SINEL'
         sinel=0

         hbpol=0.01d0
         NJmax=Jel(1)+maxval(Jel)+1
         Allocate(fBessJ(NJmax+1),fBessY(NJmax+1))

         Do ief=1,Nsed
            iM0=-Jel(1)+(ief-1)*2 !2m0
            ielm=Nef
            Do iel=2,Nel
            Do iM02=-Jel(iel),Jel(iel),2 !2m
               ielm=ielm+1
               ibess=abs(iM0-iM02)/2+1
               ftild=0
               do ib=1,Nb-1
                  Npol=nint((xb(ib+1)-xb(ib))/hbpol)
                  if(abs(Npol*1.0d0-(xb(ib+1)-xb(ib))/hbpol)>1d-3)then
                     write(*,*)'Error in interpolation at ib=',ib
                     write(*,*)'hb=',xb(ib+1)-xb(ib)
                     write(*,*)'hbpol=',hbpol
                     STOP
                  endif
                  do ipol=0,Npol-1
                     bpol=xb(ib)+ipol*hbpol

                     ctmp=(S0iM(ib+1,ielm,ief)-S0iM(ib,ielm,ief))
     a                    *ipol/Npol + S0iM(ib,ielm,ief)

!phase menant a l'amplitude de diffusion Coulombienne exacte
                     chic2=2*etaPT*log(EKPT*bpol)
                     if(abs(bpol)<1e-3)chic2=0
                     do iq=1,Nq
                        q=2*EKPT*sin(theta(iq)/2)
*                        z=q*bpol
*                        call S17DEF(0.0d0,z,NJmax,'U',BessJ,NZ,IFAIL)
*                        ftild(iq)=ftild(iq)+BessJ(ibess)*ctmp
*     a                       *bpol*exp(ci*chic2)*hbpol

*                        z=q*bpol !version w/ S17DEF
*                        if(z==0)z=1d-15
                        x=q*bpol
                        if(x==0)x=1d-15

                        If(ibess==1)then !iM02==iM0
                           call bessj0y0(x,fBessJ(1),fBessY(1))
                           ftild(iq)=ftild(iq)+fBessJ(ibess)*ctmp
     a                          *bpol*exp(ci*chic2)*hbpol
                        ElseIf(ibess==2)then !iM02=iM0+-1
                           fBessJ(2)=bessj1(x)
                           ftild(iq)=ftild(iq)+fBessJ(ibess)*ctmp
     a                          *bpol*exp(ci*chic2)*hbpol
                        Else
                           call BESSJ(NJmax,x,fBessJ,fBessY)
                           ftild(iq)=ftild(iq)+fBessJ(ibess)*ctmp
     a                          *bpol*exp(ci*chic2)*hbpol
                        Endif
                     enddo      !iq
                  enddo !bpol
               enddo !ib
               do iq=1,Nq
                  sinel(iel,iq)=sinel(iel,iq)+2*abs(ftild(iq))**2
               enddo !iq
            EndDo !iM02
            Enddo !iel
         Enddo !ief

         Deallocate(fBessJ,fBessY)

         Do iq=1,Nq
            sinel(:,iq)=EKPT**2*sinel(:,iq)/Nef/100
         Enddo

*********************************************

*ERel
c Printing of elastic scattering cross section
      Do iq=1,Nq
         if(theta(iq).gt.1e-5)then
            write(2,'(1f8.4,10(1pd16.8))')180/pi*theta(iq),
     a           sruth(iq),selast(iq),(sinel(iel,iq),iel=2,Nel)
         endif
      enddo

      DEALLOCATE(theta,sruth,fcoul,ftild)
      DEALLOCATE(S0M,S0iM,selast,sinel)

      END SELECT !ilect for reaction channel (el.scatt., total or distr. bu)

 1001 format(/,' Breakup: ',a2,i2,' -> ',a2,i2,' + ',a2,i2,/) 
 1003 format(1h ,10f7.3)
 1005 format(f11.6,4x,6(1pe14.5,2x))
 1006 format(1h ,f8.2,1p,5e13.4)
 3000 format(2f10.6,2i3)
 3001 format(3i3)
 3003 format(f8.3,1p,60d16.8)

 9998 continue
      DEALLOCATE(NRel,Jel,lel,Eel)
      DEALLOCATE(Jpo,lpo,Vp,rp,ap,rC,Vls,rls,als)
      DEALLOCATE(Vcf,foel,S,xb,fibb)

      close(1)
      close(2)


 9999 continue
*END
      END

*multeik
      subroutine multeik
c Computes the eikonal phase (including Coulomb correction)
c Note the change in the definition of Vcs and Vfs from 'reneik12.f':
c here they include i/hbar/v. In this way we can use directly the
c chi phase-shift functions of Suzuki in the variable chi
      USE cstes
      USE collision
      USE potentials
      USE reseaufo
      USE reseaueik
      USE phaseik
      implicit real*8(a-h,o-z)
      real*8, parameter :: epslu=1d-4
      real*8 :: xt(Nt),wt(Nt)
      complex*16 :: voptN,vopt,cpotct,cpotft,echin,echic,echi(3)
      complex*16 :: chincT,chinfT,chinPT
      complex*16 :: Vcs(0:Ns),Vfs(0:Ns),Vsn,Vsp
      real*8, allocatable, dimension (:):: yl,p
      real*8, allocatable, dimension (:,:):: ylm
      complex*16, allocatable, dimension (:,:):: f
***
c     Evaluate the dimension of the arrays, and allocate them
      if(mod(lamax,2)==0)then
         ny=(lamax/2+1)**2
      else
         ny=(lamax/2+1)*(lamax/2+2)
      endif
*      print*,' ny =',ny
      nydim=(Nl+1)*(Nl+2)/2

      allocate (ylm(nydim,Nt/2),feik(ny,Nr,Nb,3),f(ny,3))
      allocate (p(Nl+1),yl(nydim))
***
c     Computes the points and weights of the Gauss quadrature on theta
      ifail=0
      zero=0
!     old: call d01bcf (0,-one,one,zero,zero,Nt,wt,xt,ifail)
       x1= -1.d0
       x2= 1.d0
       call gauleg(x1,x2,xt,wt,Nt)
       if(ifail.ne.0)print*,Nt,ifail
c     Computes the Ylm needed because of symmetries
c     (ie only phi=0, mu>=0, and la+mu even)
      Do it=1,Nt/2 !Nt is even to simplify the calculation
         x=xt(it)
         call harsph (x,-(lamax+1),p,yl)
         iy=0
         do la=0,lamax
            do mu=0,la
               iy=iy+1
               ylm(iy,it)=yl(iy)
            enddo !mu
         enddo !la
      enddo ! it

      IF(npotN>=10.and.npotN<100)then
c     Precalculation of the integral over Z of the P-T nuclear interaction
c     as a function of the perpendicular distance s
c     NB integrand symmetric in Z => *2 except for Z=0
c     and invariant on the origin => no zz-z
c     Tested and OK for 11Be @70AMeV on lead for b=10-100fm, E=0-3MeV
c     with Ns=600 hs=0.05fm. Another interpolation could be better
c     eg 3pt instead of 2 to have at least the curve in the interpolation.
         Do js=0,Ns
            if(js==0)then
               s=1d-15           !at s=0 voptN diverges
            else
               s=js*hs
            endif

            if(npotN==10)then   !optical potentials
*               Vcs(js)=vopt(VcT,rRcT,aRcT,WcT,rIcT,aIcT,
*     a              WDcT,rDcT,aDcT,kZc,kZT,rCcT,s)
*               Vfs(js)=vopt(VfT,rRfT,aRfT,WfT,rIfT,aIfT,
*     a              WDfT,rDfT,aDfT,kZf,kZT,rCfT,s)

               Vcs(js)=voptN(VcT,rRcT,aRcT,WcT,rIcT,aIcT,
     a              WDcT,rDcT,aDcT,kZc,kZT,rCcT,s)
               Vfs(js)=voptN(VfT,rRfT,aRfT,WfT,rIfT,aIfT,
     a              WDfT,rDfT,aDfT,kZf,kZT,rCfT,s)
            elseif(npotN==20)then  !gaussian potentials
               Vcs(js)=vgauss(VcT,aRcT,s)
               Vfs(js)=vgauss(VfT,aRfT,s)
            endif

            do iz=1,Nz
               zz=iz*hz
               rxT=sqrt(s**2+zz**2)
               if(npotN==10)then !optical potentials
*                  Vcs(js)=Vcs(js)+2*vopt(VcT,rRcT,aRcT,WcT,rIcT,aIcT,
*     a                 WDcT,rDcT,aDcT,kZc,kZT,rCcT,rxT)
*                  Vfs(js)=Vfs(js)+2*vopt(VfT,rRfT,aRfT,WfT,rIfT,aIfT,
*     a                 WDfT,rDfT,aDfT,kZf,kZT,rCfT,rxT)

                  Vcs(js)=Vcs(js)+2*voptN(VcT,rRcT,aRcT,WcT,rIcT,aIcT,
     a                 WDcT,rDcT,aDcT,kZc,kZT,rCcT,rxT)
                  Vfs(js)=Vfs(js)+2*voptN(VfT,rRfT,aRfT,WfT,rIfT,aIfT,
     a                 WDfT,rDfT,aDfT,kZf,kZT,rCfT,rxT)
               elseif(npotN==20)then !gaussian potentials
                  Vcs(js)=Vcs(js)+2*vgauss(VcT,aRcT,rxT)
                  Vfs(js)=Vfs(js)+2*vgauss(VfT,aRfT,rxT)
               endif
            enddo !iz
         EndDo !js
! Old version where Vcs and Vfs contain only the integral of the potential
*         Vcs=Vcs*hz
*         Vfs=Vfs*hz
! New version where Vcs and Vfs contain also the -i/hbar/v factor
         Vcs=-ci*Vcs*hz/v0
         Vfs=-ci*Vfs*hz/v0
*         open(unit=27,file='Vopt.chi')
*         do js=0,Ns
*            write(27,*)js*hs,real(Vcs(js)),aimag(Vcs(is)),
*     a           real(Vfs(js)),aimag(Vfs(is))
*         enddo
*         close(unit=27)

      ELSEIF(npotN==100)then
c     The Chi phase-shift functions are precalculated from a GammaNN profile
c     function and read in separate files 101='fichcT' and 102='fichfT'.
         Do js=0,Ns
            s=js*hs
*            read(101,*)slu,Vsn,Vsp
*            write(*,*)s,Vsn,Vsp
            read(101,*)slu,VsnR,VsnI,VspR,VspI
*            write(*,*)s,VsnR,VsnI,VspR,VspI
            if(abs(slu-s)>epslu)then
               write(*,*)'Error in reading s in fichcT: ',s,slu
               STOP
            endif
*            Vcs(js)=-(Vsn+Vsp)    !old version when there was a mistake in Suzuki's phaseshift functions
*            Vcs(js)=Vsn+Vsp
            Vcs(js)=VsnR+VspR+ci*(VsnI+VspI)

*            read(102,*)slu,Vsn,Vsp
*            write(*,*)s,Vsn,Vsp
            read(102,*)slu,VsnR,VsnI,VspR,VspI
            if(abs(slu-s)>epslu)then
               write(*,*)'Error in reading s in fichfT: ',s,slu
               STOP
            endif
*            Vfs(js)=-(Vsn+Vsp)    !old version when there was a mistake in Suzuki's phaseshift functions
*            Vfs(js)=Vsn+Vsp
            Vfs(js)=VsnR+VspR+ci*(VsnI+VspI)
         EndDo
      ENDIF

      feik=0

c     calculation of the eikonal phase corrected for Coulomb at first-order
      do 3 ib=1,nb              !for each impact parameter
         if(xb(ib)==0)then
            b=1d-15             !at b=0 the Coulomb phase diverges
         else
*            b=b0+(ib-1)*hb
            b=xb(ib)
         endif
         do 4 ir=1,nr           !for each radius
            r=ir*hr
c            f=0
            do 5 it=1,Nt/2      !we only need 1/2 of the theta NB => Nt even?
               z=r*xt(it)
               sp=r*sqrt(1-xt(it)**2)
               do 6 ip=1,Np/2   !we only need 1/2 of the phi NB => Np even?
                  phi=(2*ip-1)*pi/Np
                  x=sp*cos(phi)
                  y=sp*sin(phi)
c     Calculation of the nuclear eikonal phases
         IF(npotN==0)then !ie purely point-point Coulomb potential
            echin=1             !ie not nuclear potential
*            if(lec>=20)then
*               scT=sqrt((b+cc*x)**2+(cc*y)**2)
*               sfT=sqrt((b+cf*x)**2+(cf*y)**2)
*            endif
         ELSE
            scT=sqrt((b+cc*x)**2+(cc*y)**2)
            sfT=sqrt((b+cf*x)**2+(cf*y)**2)
            ics=scT/hs
            ifs=sfT/hs

c     2-point interpolation
*            if(ics>=Ns)then
*               chi=0
*            else !interpolation
*               chi=(ics+1-scT/hs)*Vcs(ics)+(scT/hs-ics)*Vcs(ics+1)
*            endif
*            if(ifs<Ns)then
*               chi=chi+(ifs+1-sfT/hs)*Vfs(ifs)+(sfT/hs-ifs)*Vfs(ifs+1)
*            endif
c     3-point interpolation
c     Interpolation of cT phase
            If(ics>=Ns)then
               chincT=0         !beyond maximum impact parameter
            Else
               if(ics==0.or.scT/hs-ics>0.5d0) ics=ics+1
               ps=scT/hs-ics
               if(ics==Ns)then
                  chincT=ps*(ps-1)/2*Vcs(ics-1)+(1-ps**2)*Vcs(ics)
               else
                  chincT=ps*(ps-1)/2*Vcs(ics-1)+(1-ps**2)*Vcs(ics)
     a                 +ps*(ps+1)/2*Vcs(ics+1)
               endif
            Endif
c Interpolation of fT phase
            If(ifs>=Ns)then
               chinfT=0         !beyond maximum impact parameter
            Else
               if(ifs==0.or.sfT/hs-ifs>0.5d0) ifs=ifs+1
               ps=sfT/hs-ifs
               if(ifs==Ns)then
                  chinfT=ps*(ps-1)/2*Vfs(ifs-1)+(1-ps**2)*Vfs(ifs)
               else
                  chinfT=ps*(ps-1)/2*Vfs(ifs-1)+(1-ps**2)*Vfs(ifs)
     a                 +ps*(ps+1)/2*Vfs(ifs+1)
               endif
            Endif
! Old version where Vcs and Vfs contain only the integral of the potential
*            chi=-chi/v0
*            echin=exp(ci*chi)   !nulear part of the phase
! New version where Vcs and Vfs contain also the -i/hbar/v factor
! this part is now below, where the Coulomb phase is evaluated
*            echin=exp(chi)
*            if(ib==200)write(*,*)b,chi,echin
! to include computation of stripping cross section
            if(lec==7)then      !computes stripping cross section
               echin=abs(exp(chincT))**2*(1-abs(exp(chinfT))**2)
            else
               echin=exp(chincT+chinfT)
            endif
         ENDIF !npotN

c     Calculation of Coulomb eikonal phase
c     Exact Coulomb phase (pt-pt)
            chicct=1+2*cc*x/b+cc**2*(x**2+y**2)/b**2
            chicct=e2ct/v0*log(chicct)
            chicft=1+2*cf*x/b+cf**2*(x**2+y**2)/b**2
            chicft=e2ft/v0*log(chicft)

c     Coulomb with approximation of the logarithm
*            chicct=qc*x/b
*            chicft=qf*x/b

c     Coulomb (analytical integration up to Zm) (a=-monopole term)
*         Zm=10*Nz*hz
*         chicct=Zm/sqrt((b+cc*x)**2+cc**2*y**2)
*         chicct=-2*e2ct/v0*log(chicct+sqrt(chicct**2+1))
*     a        +2*e2ct/v0*log(Zm/b+sqrt((Zm/b)**2+1))
*         chicft=Zm/sqrt((b+cf*x)**2+cf**2*y**2)
*         chicft=-2*e2ft/v0*log(chicft+sqrt(chicft**2+1))
*     a        +2*e2ft/v0*log(Zm/b+sqrt((Zm/b)**2+1))

c     Coulomb V(RxT)-V(R) integrated up to Zm (in a different way)
*         Zm=10*Nz*hz
*         chicct=-2*e2ct/v0*log(b/sqrt((b+cc*x)**2+cc**2*y**2))
*     a        -2*e2ct/v0*log((Zm+sqrt(Zm**2+(b+cc*x)**2+cc**2*y**2))
*     b        /(Zm+sqrt(Zm**2+b**2)))
*         chicft=-2*e2ft/v0*log(b/sqrt((b+cf*x)**2+cf**2*y**2))
*     a        -2*e2ft/v0*log((Zm+sqrt(Zm**2+(b+cf*x)**2+cf**2*y**2))
*     b        /(Zm+sqrt(Zm**2+b**2)))

            chicPT=chicft+chicct
***      
c     echic=exp(ci*(chicct+chicft))
**                  echic=exp(ci*chicct)-ci*chicct
*                  echic=ci*(chicct-chicft) !Coulomb phase
*                  echic=exp(echic)-echic !Subtraction of the 1st-order
c Coulomb phase and subtraction of the 1st-order of the eikonal
            If(lec<10)then      !Usual Eikonal => no correction
               if(lec==7)then   !computes stripping cross section
                  echic=1
               else
                  echic=exp(ci*chicPT)
               endif
c     test without nuclear and/or Coulomb interaction
*                  echin=1
*                  echic=1
               echi(1)=echin*echic !total phase
               echi(2)=0
               echi(3)=0
            ElseIf(lec<20)then  !Modified Eikonal
               echic=exp(ci*chicPT)-ci*chicPT
*               echic=exp(ci*chicPT) !sans correction !supprimer cor2 et cor3
*               echic=exp(-ci*chicPT)+ci*chicPT !signes de Suzuki (ntypo=0)
*               echic=exp(ci*chicPT) !test pr verifier signes
c     test without nuclear and/or Coulomb interaction
*                  echin=1
*                  echic=1
               echi(1)=echin*echic !total phase
               echi(2)=echin*x  !part to be corrected with K1
               echi(3)=echin*z  !part to be corrected with K0
*            ElseIf(lec<30)then  !Modified Eikonal with correction scT>str
*               if(scT<=str)then ! no correction
*                  echic=exp(ci*chicPT)
*c     test without nuclear and/or Coulomb interaction
**                  echin=1
**                  echic=1
*                  echi(1)=echin*echic !total phase
*                  echi(2)=0
*                  echi(3)=0
*               else
*                  echic=exp(ci*chicPT)-ci*chicPT
**               echic=exp(ci*chicPT) !sans correction !supprimer cor2 et cor3
**               echic=exp(-ci*chicPT)+ci*chicPT !signes de Suzuki (ntypo=0)
**               echic=exp(ci*chicPT) !test pr verifier signes
*c     test without nuclear and/or Coulomb interaction
**                  echin=1
**                  echic=1
*                  echi(1)=echin*echic !total phase
*                  echi(2)=echin*x !part to be corrected with K1
*                  echi(3)=echin*z !part to be corrected with K0
*               endif
            ElseIf(lec<110)then  !First-order approximation
c     ie calculation without nuclear and/or Coulomb interaction
               echi(1)=1 !total phase
               echi(2)=x  !part to be corrected with K1
               echi(3)=z  !part to be corrected with K0
            EndIf
*                  write(*,*)echi

c     test without any interaction
*                  echi(1)=1
*                  echi(2)=0
*                  echi(3)=0
c     Multipolar expansion of the eikonal phase
            iy=0
            iye=0
            iyo=0
            do la=0,lamax       !projection on spherical harmonics
               do mu=0,la
                  iy=iy+1       !=(la,mu)
                  if(mod(la+mu,2)==0)then
                     iye=iye+1
*      if(ib==1.and.ip==1.and.it==1)write(*,*)la,mu,iye
                     feik(iye,ir,ib,1:2)=feik(iye,ir,ib,1:2)
     a                    +wt(it)*ylm(iy,it)
     b                    *cos(mu*phi)*echi(1:2)*8*pi/Np
                  else
                     iyo=iyo+1
*      if(ib==1.and.ip==1.and.it==1)write(*,*)la,mu,iyo
                     feik(iyo,ir,ib,3)=feik(iyo,ir,ib,3)
     a                    +wt(it)*ylm(iy,it)
     b                    *cos(mu*phi)*echi(3)*8*pi/Np
                  endif
               enddo
            enddo
 6       continue               !ip
 5    continue                  !it
c      feik(:,ir,ib,:)=f(:,:)*8*pi/Np !integral on the whole sphere
      feik(1,ir,ib,1)=feik(1,ir,ib,1)-sqrt(4*pi)
c     subtract the (unsignificant) part which'd correspond to echi=1
 4    continue                  !ir
 3    continue !ib
      deallocate(ylm,f,p,yl)
      return
 1000 format(1h ,10f12.6)
      end

*calculFBU
      subroutine calculFBU

      USE cstes
      USE reseaueik
      USE collision
      USE boundstates
      USE phaseik
      USE Sbu

      implicit real*8 (a-h,o-z)
*      real*8 :: fBessJ(30),fBessY(30)
      real*8, allocatable, dimension(:):: fBessJ,fBessY
*     complex*16 :: ctmp,z,cBessJ(30) !old version using S17DEF
      complex*16 :: ctmp

      fbu=0
*      hbpol=0.05d0
*      Nbesmax=(2*Nl+iS+Jel(1))/2 !old version with S17DEF
      Nbesmax=(2*Nl+iS+Jel(1))/2+1 !version using BESSJ from PHYSTM

      Allocate(fBessJ(Nbesmax+1),fBessY(Nbesmax+1))
      
      do iq=1,Nq
         q=2*EKPT*sin(theta(iq)/2)
         do ib=1,Nb-1
*            bleft=b0+(ib-1)*hb
            bleft=xb(ib)

*            Npol=hb/hbpol
            Npol=nint((xb(ib+1)-xb(ib))/hbpol)
            if(abs(Npol*1.0d0-(xb(ib+1)-xb(ib))/hbpol)>1d-3)then
               write(*,*)'Error in interpolation at ib=',ib
               write(*,*)'hb=',xb(ib+1)-xb(ib)
               write(*,*)'hbpol=',hbpol
               STOP
            endif

            do ipol=1,Npol
               bpol=bleft+(ipol-1)*hbpol
               if(abs(bpol)<1d-3)then
                  chic2=0
               else
*                  tmpchic=2.0d0*bpol/abs(v0)/(Tout-Tin)
*                  chic2=2*etaa*log((1+sqrt(1+tmpchi**2))/tmpchi)
c     The Coulomb scattering phase is chosen equal to the the
c     Coulomb eikonal phase that reproduces the exact Coulomb amplitude
c     cf. Bertulani & Danielewicz pg 423 Eq.(8.77).
                  chic2=2*etaPT*log(EKPT*bpol)
               endif
*               z=q*bpol
*               call S17DEF(0,z,Nbesmax,'U',cBessJ,NZ,IFAIL)
               x=q*bpol
               call BESSJ(Nbesmax,x,fBessJ,fBessY)
               DO ief=1,Nsed 
                  iM0=-Jel(1)+(ief-1)*2
                  jlJM=0
                  Do l=0,Nl
                     Do iJ=abs(2*l-iS),2*l+iS,2
                        Do iM=-iJ,iJ,2
                           jlJM=jlJM+1
                           ibess=abs(iM-iM0)/2+1
*                           Bessel=real(cBessJ(ibess)) !version w/S17DEF
                           do iE=1,NE
c     Two point interpolation of Skljm
                              ctmp=(SklJM(jlJM,iE,ief,ib+1)
     a                             -SklJM(jlJM,iE,ief,ib))*(ipol-1)/Npol
     b                             +SklJM(jlJM,iE,ief,ib)
c     Integral over b
*                              fbu(iq,ief,jlJM,iE)=
*     a                             fbu(iq,ief,jlJM,iE)
*     b                             +Bessel*ctmp*bpol*hbpol !version w/S17DEF
*     c                             *exp(ci*chic2)
                              fbu(iq,ief,jlJM,iE)=
     a                             fbu(iq,ief,jlJM,iE)
     b                             +fBessJ(ibess)*ctmp*bpol*hbpol
     c                             *exp(ci*chic2)
                           enddo !iE
                        EndDo !iM
                     EndDo !iJ
                  EndDo !l
               ENDDO !ief
            enddo !ipol
         enddo !ib
      enddo !iq
      Deallocate(fBessJ,fBessY)
      end

*vJ
      function vJ(r,l,iJ)
c Effective potential between the core and the fragment in the l J channel
c It contains nuclear (Woods-Saxon or Gaussian), Coulomb (point-sphere) and centrifugal
      USE cstes
      USE collision
      USE potentials
      implicit real*8 (a-h,o-z)

      jls=iJ*(iJ+2)-4*l*(l+1)-iS*(iS+2)

      IF(mod(ntypo,10)==1)then      !Different potentials for various partial waves
         ipo=1
         do
            if(l==lpo(ipo).and.iJ==Jpo(ipo).or.ipo==Njpo)then
c     special channels or all others
               exit
            else
               ipo=ipo+1
            endif
         enddo
      ELSEIF(mod(ntypo,10)==2)then      !Different potentials for odd and even waves
         ipo=2-mod(l,2)
      ELSE
         write(*,*)'Uprogrammed potential'
         STOP
      ENDIF

      IF(ntypo<10)then !Woods-Saxon form
         eve=exp(-(r-rp(ipo))/ap(ipo))
         vJ=Vp(ipo)*eve/(1+eve)+l*(l+1)/r**2/deuxm
         if(Vls(ipo)/=0)then
            eve=exp(-(r-rls(ipo))/als(ipo))
            vJ=vJ-jls*Vls(ipo)*eve/(1+eve)**2/als(ipo)/r/8
         endif
      ELSEIF(ntypo<20)then !Gaussian form
         eve=(r-rp(ipo))/ap(ipo)
         eve=exp(-eve**2/2)
         vJ=Vp(ipo)*eve+l*(l+1)/r**2/deuxm
         if(Vls(ipo)/=0)then !the ls term is used for the derivative
            eve=(r-rls(ipo))/als(ipo)
            eve=(r-rls(ipo))**2*exp(-eve**2/2) !(r-r0)^2*Gaussian form
            vJ=vJ+Vls(ipo)*eve
         endif
      ELSE
         write(*,*)'Unprogrammed potential'
         STOP
      ENDIF

      If (kZf/=0) then
c Point-sphere Coulomb potential
         if(r>=rC(ipo))then
            c1=kZf*kZc*e2/r
            vJ=vJ+c1
         else
            c1=kZf*kZc*e2/2/rC(ipo)
            c2=3-(r/rC(ipo))**2
            vJ=vJ+c1*c2
         endif
      Endif
      return
      end

*vSW
      function vSW(V,r0,a,r)
c Computes the value of a Woods-Saxon potential in r
      implicit none
      real*8 vSW,V,r0,a,r,eve
      eve=exp((r0-r)/a)
      vSW=V*eve/(1+eve)
      return
      end

*dvSW
      function dvSW(V,r0,a,r)
c Computes - the first derivative of a Woods-Saxon potential in r
      implicit none
      real*8 dvSW,V,r0,a,r,eve
      eve=exp((r0-r)/a)
      dvSW=V*eve/(1+eve)**2
      return
      end

*vopt
      function vopt(V,r0R,aR,W,r0I,aI,WD,r0D,aD,kZ1,kZ2,rC,r)
c Computes the optical (complex) potential at r.
c Woods-Saxon form factor is assumed, and the Coulomb potential is point-sphere
      USE cstes
      USE collision
      implicit real*8 (a-h,o-z)
      complex*16 vopt

      vopt=vSW(V,r0R,aR,r)+ci*vSW(W,r0I,aI,r)
      If(WD/=0) vopt=vopt+ci*dvSW(WD,r0D,aD,r)
      If(kZ1/=0.and.kZ2/=0)then
c     Coulomb point-sphere potential
         if(r>=rC)then
            c1=kZ1*kZ2*e2/r
            vopt=vopt+c1
         else
            c1=kZ1*kZ2*e2/2/rC
            c2=3-(r/rC)**2
            vopt=vopt+c1*c2
         endif         
      Endif
      return
      end

*voptN
      function voptN(V,r0R,aR,W,r0I,aI,WD,r0D,aD,kZ1,kZ2,rC,r)
c Computes the nuclear part of the optical (complex) potential at r.
c Woods-Saxon form factor is assumed, and the Coulomb potential is point-sphere
c The pure Coulomb is subtracted
      USE cstes
      USE collision
      implicit real*8 (a-h,o-z)
      complex*16 voptN

      voptN=vSW(V,r0R,aR,r)+ci*vSW(W,r0I,aI,r)
      If(WD/=0) voptN=voptN+ci*dvSW(WD,r0D,aD,r)
      If(kZ1/=0.and.kZ2/=0)then
c     Coulomb point-sphere potential minus pure Coulomb
         if(r<rC)then
            c2=3-(r/rC)**2
            voptN=voptN+kZ1*kZ2*e2*(c2/2/rC-1/r)
         endif         
      Endif
      return
      end

*vgauss
      function vgauss(V,a,r)
c Computes at r a Gaussian potential centered in 0 with diffuseness a,
c and amplitude V
      implicit none
      real*8 vgauss,V,a,r
      vgauss=V*exp(-a*r**2)
      return
      end

*super
      subroutine super (n,h,ins,e0,nrad,u,s,eps)
*
*     n: number of discretization points
*     h: step
*     ins: singularity modification at the origin
*     e0: approxaimation of the bound-state energy
*     nrad: number of nodes
*     u: potential multiplied by 2m/hbar^2
*     s: auxiliary vector (wave function, potential modification)
*     eps: accuracy
*     
      implicit real*8 (a-h,o-z)
      dimension u(n),v(n),s(n)

      e=e0
      v=u                       !utilisation d'un vecteur tampon
      s2=0 
*
*     calculation of bound-state wave function
*     s: wave function
*
      call num1l(n,h,e,s2,v,s,nrad,eps)
*      write(*,*)'Super:Energie de l''etat lie:',nrad,e
*
*     integration of psi^2 by Adams formula
*     (Abramowicz and Stegun: 25.5.5)
*
      s1=s(3) 
      s2=s(2) 
      s3=s(1) 
      k=2*ins+1
      s(1)=log(h*s(1)**2/k)
      ss=2*h*s(2)**2/k
      s(2)=log(ss) 
      ss=ss+h/24*(9*s(3)**2+19*s2**2-5*s3**2) 
      s(3)=log(ss) 
      do 3 j=4,n
         ss=ss+h/24*(9*s(j)**2+19*s1**2-5*s2**2+s3**2) 
         s3=s2 
         s2=s1 
         s1=s(j) 
         s(j)=log(ss) 
 3    continue
*
*     5-point second derivative of the logarithm of the integral
*     (Abramowicz and Stegun: 25.3.24)
*
      s2=s(1)
      s1=s(2)
      s(1)=-k/h**2
      s(2)=s(1)/4
      do 4 j=5,n
         s3=s2
         s2=s1
         s1=s(j-2)
         s(j-2)=(-(s(j)+s3)+16*(s(j-1)+s2)-30*s1)/(12*h**2)
 4    continue
      s(n-1)=0
      s(n)=0
      do 5 j=1,n
         u(j)=u(j)-s(j)*2
 5    continue
      ins=ins+2
      return
      end

*superpotu
      subroutine superpotu (n,h,ins,e0,u,eps)
c Computes the SuSy partner of a potential u by removing the bound state
c of approximate energy e0.
c Using the potential u given in imputs, it finds the state to be removed
c using 'num1l'.
c The forbidden state is then removed using the pair of SuSy transformations
c given in D.Baye, PRL 58, 2738 (1987), and D.Baye, J.Phys.A 20, 5529 (1987).
c The new potential is sent in output as the new u.
*     n: number of discretization points
*     h: step
*     ins: singularity modification at the origin
*     e0: approximation of the bound-state energy
*     u: potential multiplied by 2m/hbar^2
*     eps: accuracy
*     
      implicit real*8 (a-h,o-z)
      parameter (Nradmax=10,epsrel=5d-2)
      dimension u(n),v(n),s(n)

*
*     calculation of bound-state wave function
*     s: wave function
      nrad=0
      do
         e=e0
         v=u                    !use of a buffer for v is modified by num1l
         s2=0 
         call num1l(n,h,e,s2,v,s,nrad,eps)
         if(nrad>=0.and.abs(1-e/e0)<epsrel)then
*            write(*,*)'Superpotu: E_eli=',nrad,e
            exit
         elseif(nrad>Nradmax)then
            write(*,*)'Forbidden state not found: E0i=',e0
            STOP
         else
            nrad=nrad+1
         endif
      enddo
*
*     integration of psi^2 by Adams formula
*     (Abramowicz and Stegun: 25.5.5)
*
      s1=s(3) 
      s2=s(2) 
      s3=s(1) 
      k=2*ins+1
      s(1)=log(h*s(1)**2/k)
      ss=2*h*s(2)**2/k
      s(2)=log(ss) 
      ss=ss+h/24*(9*s(3)**2+19*s2**2-5*s3**2) 
      s(3)=log(ss) 
      do 3 j=4,n
         ss=ss+h/24*(9*s(j)**2+19*s1**2-5*s2**2+s3**2) 
         s3=s2 
         s2=s1 
         s1=s(j) 
         s(j)=log(ss) 
 3    continue
*
*     5-point second derivative of the logarithm of the integral
*     (Abramowicz and Stegun: 25.3.24)
*
      s2=s(1)
      s1=s(2)
      s(1)=-k/h**2
      s(2)=s(1)/4
      do 4 j=5,n
         s3=s2
         s2=s1
         s1=s(j-2)
         s(j-2)=(-(s(j)+s3)+16*(s(j-1)+s2)-30*s1)/(12*h**2)
 4    continue
      s(n-1)=0
      s(n)=0
      do 5 j=1,n
         u(j)=u(j)-s(j)*2
 5    continue
      ins=ins+2
      return
      end

*jfromlSJM
      integer function jfromlSJM(l,iS,iJ,iM)
c     Donne le compteur jlJM a partir de l,iS,iJ et iM
      implicit none
      integer :: l,iS,iJ,iM
      integer :: it,jl,jJ,jM
      jfromlSJM=l*(l+1)*iS+(iJ**2-iS**2)/4+(iJ+iM)/2+1
c      verification du calcul
*      it=0
*      do jl=0,2*l-2,2
*         do jJ=abs(jl-iS),jl+iS,2
*            do jM=-jJ,jJ,2
*               it=it+1
*            enddo
*         enddo
*      enddo
*      do jJ=abs(2*l-iS),iJ-2,2
*         do jM=-jJ,jJ,2
*            it=it+1
*         enddo
*      enddo
*      it=it+(iJ+iM)/2+1
*      if(it/=jfromlSJM)write(*,*)it,jfromlSJM
      return
      end

*th
      function th(x)
c Computes tanh(x) for large positive x
      implicit none
      real*8 th,x
      th=(1-exp(-2*x))/(1+exp(-2*x))
      end

*phase
      function phase(z)
c     Calcule la phase du nombre complexe z comprise entre -pi et pi
      implicit none
*      real*8, parameter :: pi=3.141592653589793238d0
      real*8 :: pi,phase,ph
      complex*16 :: z
      pi=acos(-1d0)
      ph=atan(aimag(z)/real(z))
      If(real(z)<0.and.aimag(z).ge.0) then
         ph=ph+pi
      Elseif(real(z)<0.and.aimag(z)<0) then
         ph=ph-pi
      Endif
      phase=ph
      return
      end

*c2l
      character(len=2) function c2l(n)
c Gives the 2 letters equivalent of a positive integer<100
      implicit NONE
      integer n,nu,nd
      character lu,ld
      nu=48+mod(n,10)
      lu=char(nu)
      nd=n/10
      nd=48+mod(nd,10)
      ld=char(nd)
      c2l=ld//lu
      return
      end
*BESKN
      SUBROUTINE BESKN(XX,N,FP,YP) 
      IMPLICIT real*8 (A-H,O-Z) 
      DIMENSION FP(30),YP(30),C1(5,2),C2(5,2)
      LOGICAL MAV,XL0,B1E,B2E
      SAVE
      DATA PIDEMI/1.57079633D0/,PI/3.14159265D0/
      DATA C1/106.88679D0,5.8486996D0,.54284668D0,.10351563D0,.0625D0
     1 ,-132.76182D0,-7.7399540D0,-.80529785D0,-.19335938D0,-.1875D0/ 
      DATA C2/ 23.75128D0,1.6380462D0,.20957031D0,.651041666D-1,.125D0
     1 ,-30.624012D0,-2.3709870D0,-.37089844D0,-.16406250D0,-.375D0/  
C
      LL=1  
      NN=ABS(N)
      XL0=SIGN(1d0,XX).LT.0
      X=ABS(XX) 
      IF(X.EQ.0)GOTO16 
      Y=2/X
      M=X 
      MAV=X.GT.10.0D0.AND.NN.LE.(M/8+4)
      B1E=X.GT.10.0D0.OR..NOT.XL0
      B2E=X.GT.1.2D0
      IF(MAV)GOTO1
      M=MAX(NN+10,M+25)
   15 M=M/2 
      YNU=0
      I=M+M 
      FLI=I 
      FL2=FLI/2
      LIM=MAX(20,NN)/2*2 
      LM=(I-1)/LIM+1  
      K=MOD(I-1,LIM)+1
      FP(1)=1.D-250 
      FP(2)=0
      ALF=0
      DO5J=1,LM 
      I=K
      FP(K+1)=FP(1) 
      FP(K+2)=FP(2) 
      DO4L=1,K
      FP(I)=FLI*FP(I+1)*Y+FP(I+2) 
      ALF=ALF+FP(I) 
      I=I-1
    4 FLI=FLI-1
      IF(B1E.OR.B2E)GOTO5
      M=K/2
      I=M
      DO7L=1,M
      YNU=YNU+FP(2*I+1)/FL2
      I=I-1
    7 FL2=FL2-1
    5 K=LIM 
      IF(X.GT.10.D0)GOTO1 
      E=EXP(X)
      ALF=E/(ALF+ALF-FP(1))
   11 DO6I=1,NN 
    6 FP(I)=FP(I)*ALF 
      GOTO18
    1 UX2=1/X**2 
    8 ALF=0
      PHI=0
      DO9I=1,5
      ALF=(ALF+C1(I,LL))*UX2
    9 PHI=(PHI+C2(I,LL))*UX2
      PHI=(PHI+1)*X
      ALF=(ALF+1)/SQRT(X*PIDEMI)
      SH=EXP(PHI)
      IF(XL0)YP(LL)=ALF/SH*PIDEMI
      ALF=ALF*(SH-1/SH)/2
      LL=LL+1 
      IF(MAV)GOTO12 
      ALF=ALF/FP(LL-1)
      IF(XL0.AND.LL.EQ.2)GOTO8
      GOTO11
   12 FP(LL-1)=ALF
      IF(LL.EQ.2)GOTO8
      IF(NN.LE.2)RETURN 
      FLI=0
      DO13I=3,NN
      FLI=FLI+Y 
   13 FP(I)=-FLI*FP(I-1)+FP(I-2) 
   18 IF(.NOT.XL0)RETURN
      IF(B1E)GOTO20
      IF(B2E)GOTO45
      YNU=(0.115931515D0-LOG(X))*FP(1)+2*YNU*ALF
      YP(1)=YNU
   44 YP(2)=(1/X-YP(1)*FP(2))/FP(1)
   20 FLI=0
      IF(NN.LE.2)RETURN 
      DO23I=3,NN
      FLI=FLI+1
   23 YP(I)=Y*FLI*YP(I-1)+YP(I-2) 
      RETURN
   45 Z=1/X
      ANUM=1+Z*(1.04375D1+Z*(3.50585937D1+Z*(4.51474609D1+Z*(
     1 2.07832947D1+Z*(2.52649784D0+Z*3.08763385D-2)))))
      ADEN=1+Z*(1.05625D1+Z*(3.63085937D1+Z*(4.90166016D1+Z*(
     1 2.50188904D1+Z*(3.90920162D0+Z*1.04710758D-1)))))
      YP(1)=1.25331414D0*ANUM/ADEN*SQRT(Z)/E
      GOTO44
   16 ZERO=0
      DO17I=1,NN
      IF(XL0)YP(I)=1/ZERO
   17 FP(I)=0
      FP(1)=1
      IF(XL0)PRINT1001
      RETURN
 1001 FORMAT(20H KN(0) N''EXISTE PAS ) 
      END
*intbes
      subroutine intbes(b0,q,eta,xk,nla,cres)
      implicit real*8(a,b,d-h,o-z)
      implicit complex*16(c)
      dimension ug(100),wg(100),cy(10),cy2(10)
      data ci/(0,1)/
c     old: call d01bcf(-3,0.0d0,1.0d0,0.0d0,0.0d0,nla,wg,ug,ifail)
       x1= 0.d0
       x2= 1.d0
       call gauleg(x1,x2,ug,wg,nla)
      cres1=0
      cres2=0
      do i=1,nla
         cb1=cmplx(b0,ug(i),16)
         cb2=cmplx(b0,-ug(i),16)
         cz1=exp(2*ci*eta*log(xk*cb1))
         cz2=exp(2*ci*eta*log(xk*cb2))
c       call s17dlf(1,0.0d0,q*cb1,1,'U',cy,nz,ifail)
c       call s17dlf(2,0.0d0,q*cb2,1,'U',cy2,nz,ifail)
         cres1=cres1+wg(i)*cy(1)*cz1/cb1
         cres2=cres2+wg(i)*cy2(1)*cz2/cb2
      end do
      cres=ci*(cres1-cres2)/2
      return
      end
*intbesm
      subroutine intbesm(mu,b0,q,eta,nla,cres)
      USE cstes
c     computes the integral over b beyond b_max needed to converge
c     when computing the elastic-scattering cross section.
c     Idea: replace the Hankel functions by their long-range expansion?
c     Is that valid for very small angle?
c     Do we need it if only breakup is analysed?
      implicit real*8(a,b,d-h,o-z)
      implicit complex*16(c)
      dimension ug(100),wg(100),cy1(0:10),cy2(0:10)
!     old: call d01bcf(-3,0.0d0,1.0d0,0.0d0,0.0d0,nla,wg,ug,ifail)
       x1= 0.d0
       x2= 1.d0
       call gauleg(x1,x2,ug,wg,nla)
      cres1=0
      cres2=0
      do i=1,nla
         cb1=cmplx(b0,ug(i),16)
         cb2=cmplx(b0,-ug(i),16)
         cz1=exp(2*ci*eta*log(cb1/2))
         cz2=exp(2*ci*eta*log(cb2/2))
!     Asymptotic expression of Hankel functions (in progress)
c         cy1(1)=sqrt(2/pi/q*cb1)*exp(-ci*(q*cb1-pi/4))/2
c         do i=2,Nbesmax
c         cy1(i)=ci*BessJ(i-1)
c         enddo         
!      old: call s17dlf(1,0.0d0,q*cb1,mu+1,'U',cy1,nz,ifail)
!      old: call s17dlf(2,0.0d0,q*cb2,mu+1,'U',cy2,nz,ifail)
         cres1=cres1+wg(i)*cy1(mu)*cz1/cb1
         cres2=cres2+wg(i)*cy2(mu)*cz2/cb2
      end do
      cres=ci*(cres1-cres2)/2
      return
      end
*DEPHASE
      SUBROUTINE DEPHASE0(MAX,H,W,Y,EPS,DELTA,L,ETA,QK,RFIN) 
c my version of dephase with a call of coufra(*,*,0,L etc)
c this helps getting a better convergence at low energy
      IMPLICIT REAL*8 (A-H,O-Z) 
      DIMENSION FG(2),DFC(101),GC(101),DGC(101),FC(101),W(1),Y(1) 
      DD=0 
      PAS=H*EPS 
      QQ=QK*QK
      LL=L+1
      L1=L*LL 
      R=2*ETA*QK
      RMAX=MAX*H
      DO20 K=2,MAX
      RMAX=RMAX-H 
      IF(ABS(W(MAX-K+1)-(L1/RMAX+R)/RMAX).LT.EPS)GOTO20 
      IF(K.GT.3)GOTO21
      PRINT2001 
      STOP  
   20 CONTINUE
      RMAX=H*(MAX-2)  
   21 H2 = H**2 
      H212 = H2/12 
      AA=H212*(QQ-W(1)) 
      Y(1)=H**LL
      B0=0  
      B1=Y(1)*(1+AA)  
      DO22 K=2,MAX
      AA=H212*(QQ-W(K)) 
      B2=12*Y(K-1)-10*B1-B0 
      Y(K)=B2/(1+AA)  
      B0=B1 
      B1=B2 
   22 CONTINUE
      R=RMAX
      N1=INT((RMAX+1d-5)/H)+1 
      CALL COUFRA(QK*R,ETA,0,L,FC,DFC,GC,DGC) 
      FG(1)=FC(LL)
      FG(2)=GC(LL)
      DO23 N=N1,MAX 
      R=R+H 
      CALL COUFRA(QK*R,ETA,0,L,FC,DFC,GC,DGC) 
      C=FG(2)*Y(N)-GC(LL)*Y(N-1)
      D=FC(LL)*Y(N-1)-FG(1)*Y(N)
      Z=FG(1)*GC(LL)-FG(2)*FC(LL) 
      DEL=ABS(Z)/SQRT(C*C+D*D)  
      IF(ABS(1-DD/DEL).LT.PAS) GO TO 24 
      DD = DEL
      FG(1)=FC(LL)
      FG(2)=GC(LL)
   23 CONTINUE
      PRINT2000,L,QK  
   24 RFIN=R
      D=D/C 
      IF(Y(N)*(FC(LL)+GC(LL)*D).LT.0.)DEL=-DEL
      DELTA=ATAN(D) 
      DO25 K=1,MAX
   25 Y(K)=Y(K)*DEL 
      RETURN
 2000 FORMAT(/27H PAS DE CONVERGENCE POUR L=,I3,5H   K=,F10.6/) 
 2001 FORMAT(/46H LE POTENTIEL N'EST JAMAIS PUREMENT COULOMBIEN)
      END
*CLEBS
      SUBROUTINE CLEBS (L1,L2,L3,M1,M2,M3,Q)
      IMPLICIT REAL*8 (A-H,O-Z) 
      DIMENSION FT(0:1000)
      SAVE  
      DATALMEM,(FT(I),I=0,10)/9,2*1.0d0,2.0d0,6.0d0,24.0d0,
     1 120.0d0,720.0d0,5040.0d0,40320.0d0,362880.0d0,3628800.0d0/
      IS=(L1+L2+M1-M2)/2
      K3=-M3
      Q1=L3+1 
   12 Q=0  
      I1=L1 
      I2=L2 
      I3=L3 
      K1=M1 
      K2=M2 
      IF(K1+K2+K3.NE.0)RETURN 
      L=I1+I2+I3
      IF(MOD(L,2).EQ.1)GOTO8
      L=L/2 
      IF(L.LE.LMEM)GOTO6
      DO10I=LMEM,L
   10 FT(I+1)=(I+1)*FT(I) 
      LMEM=L
    6 J1=ABS(K1) 
      J2=ABS(K2) 
      J3=ABS(K3) 
      IF(I1.LT.ABS(I2-I3).OR.I1.GT.I2+I3)RETURN
      IF(J1+J2.EQ.0)GOTO11
      J1=I1-J1
      J2=I2-J2
      J3=I3-J3
      IF(J1)8,2,1 
    2 IF(I1.NE.0)GOTO1
      IF(J2.LT.0)GOTO8
   13 IF(J3.LT.0)GOTO8
    4 Q=SQRT(Q1/(I2+1)) 
      IS=IS+(I2-K2)/2 
      IF(MOD(IS,2).EQ.1)Q=-Q
      RETURN
    1 IF(J2.GT.J1)GOTO3 
      IF(J2.LT.0)GOTO8
      IS=IS+L 
      J1=J2 
      K1=K2 
      K2=M1 
      I1=I2 
      I2=L1 
      IF(I1.EQ.0)GOTO13 
    3 IF(J3.GT.J1)GOTO5 
      IF(J3.LT.0)GOTO8
      IS=IS+L 
      J1=K3 
      K3=K1 
      K1=J1 
      I3=I1 
      I1=L3 
      J1=J3 
      IF(I1.EQ.0)GOTO4
    5 IF(K1.GE.0)GOTO9
      K1=-K1
      K2=-K2
      K3=-K3
      IS=IS+L 
    9 CONTINUE
      Q1=Q1*FT(L-I3)/FT(L-I1)/FT(L-I2)/FT(L+1)
      I2=(I2+K2)/2
      I3=(I3+K3)/2
      K2=I2-K2
      K3=I3-K3
      J1=J1/2 
      I1=J1+K1
      J2=I3-K2
      J3=MAX(J2,0) 
      IS=IS+I1+K2 
      X=0  
      DO7I=J3,J1
    7 X=-X+FT(I1+I)*FT(I2+I3-I)/FT(J1-I)/FT(I3-I)/FT(I-J2)/FT(I)
      Q=X*       SQRT(Q1*FT(J1)*FT(K2)*FT(K3)*FT(I3)/FT(I1)/FT(I2)) 
      IF(MOD(IS,2).EQ.1)Q=-Q
      RETURN
    8 Q=0  
      PRINT1010,L1,L2,L3,M1,M2,M3 
 1010 FORMAT(10H ERREUR 3J,2(3X,3I3)) 
      RETURN
   11 IF(MOD(L,2).EQ.1)RETURN 
      I1=L-I1 
      I2=L-L2 
      I3=L-L3 
      Q=SQRT(FT(I1)*FT(I2)*FT(I3)/FT(L+1)*Q1) 
      I1=I1/2 
      I2=I2/2 
      I3=I3/2 
      L =L/2
      Q=Q*FT(L )/FT(I1)/FT(I2)/FT(I3) 
      IF(MOD(L +IS,2).EQ.1)Q=-Q 
      RETURN
      ENTRY TROISJI (L1,L2,L3,M1,M2,M3,Q) 
      ENTRY TROISJ (L1,L2,L3,M1,M2,M3,Q)  
      K3=M3 
      Q1=1 
      IS=0  
      GOTO12
      END 
*SIXJ 
      SUBROUTINE SIXJ(J1,J2,J3,L1,L2,L3,Q)
      IMPLICIT REAL*8 (A-H,O-Z) 
      DIMENSION M(7),M1(4),M2(4),M3(4),FT(0:1000)
      LOGICAL ICE 
      COMMON/TEUFJ/ICE
      SAVE  
      DATA LMEM,(FT(I),I=0,10)/9,2*1.0d0,2.0d0,6.0d0,24.0d0,
     1 120.0d0,720.0d0,5040.0d0,40320.0d0,362880.0d0,3628800.0d0/  
      ENTRYSIXJI(J1,J2,J3,L1,L2,L3,Q) 
      I1=J1 
      I2=J2 
      I3=J3 
      K1=L1 
      K2=L2 
      K3=L3 
      IS=0  
   24 Q=0  
      M(1)=I1+I2+I3 
      M(2)=I1+K2+K3 
      M(3)=K1+I2+K3 
      M(4)=K1+K2+I3 
      DO 17 I=1,4 
      IF(MOD(M(I),2).EQ.1) GO TO 8
   17 CONTINUE
      L=MAX(I1+I2+K1+K2,I1+I3+K1+K3,I2+I3+K2+K3) 
      L=L/2 
      IF(L.LE.LMEM) GO TO 6 
      DO 10 I=LMEM,L  
   10 FT(I+1)=FT(I)*(I+1) 
      LMEM=L
    6 IF(I1.LT.ABS(I2-I3).OR.I1.GT.I2+I3) RETURN 
      IF(I1.LT.ABS(K2-K3).OR.I1.GT.K2+K3) RETURN 
      IF(K1.LT.ABS(I2-K3).OR.K1.GT.I2+K3) RETURN 
      IF(K1.LT.ABS(K2-I3).OR.K1.GT.K2+I3) RETURN 
      IF(I1) 8,2,1
    2 IF(I2.LT.0) GO TO 8 
    9 IF(I3.LT.0) GO TO 8 
   14 IF(K1.LT.0) GO TO 8 
   19 IF(K2.LT.0) GO TO 8 
   23 IF(K3.LT.0) GO TO 8 
   27 Q=SQRT(1.0d0/(I2+1)/(K2+1))  
      IS=(I2+K2+K1)/2+IS
      IF(MOD(IS,2).EQ.1) Q=-Q 
      RETURN
    1 IF(I1.GT.1) GO TO 3 
      IF(I2.LT.0) RETURN
   12 IF(I3.LT.0) RETURN
   16 IF(K1.LT.0) RETURN
   21 IF(K2.LT.0) RETURN
   25 IF(K3.LT.0) RETURN
   28 IF(I2.LT.I3) GO TO 4
      IC=I2 
      I2=I3 
      I3=IC 
      IC=K2 
      K2=K3 
      K3=IC 
    4 IF(K2.GT.K3) GO TO 5
      I11=I1+K1+I2-K2 
      I11=I11/2 
      I12=I11-I2+K2 
      Q=SQRT(I11*I12*1./I3/(I3+1)/K3/(K3+1))
      IS =I11+K2+IS 
      IF(MOD(IS,2).EQ.1) Q=-Q 
      RETURN
    5 I11=K3-K1+I2
      I11=I11/2+1 
      I12=I11+K1+1
      Q=SQRT(I11*I12*1.0d0/I3/(I3+1)/K2/(K2+1))
      IS =I12-1+IS
      IF(MOD(IS ,2).EQ.1) Q=-Q  
      RETURN
    3 IF(I2.GE.I1) GO TO 7
      IF(I2.LT.0) GO TO 8 
      IC=I2 
      I2=I1 
      I1=IC 
      IC=K1 
      K1=K2 
      K2=IC 
      IF(I1.EQ.0) GO TO 9 
      IF(I1.EQ.1) GO TO 12
    7 IF(I3.GE.I1) GO TO 13 
      IF(I3.LT.0) GO TO 8 
      IC=I3 
      I3=I1 
      I1=IC 
      IC=K3 
      K3=K1 
      K1=IC 
      IF(I1.EQ.0) GO TO 14
      IF(I1.EQ.1) GO TO 16
   13 IF(K1.GE.I1) GO TO 18 
      IF(K1.LT.0) GO TO 8 
      IC=K1 
      K1=I1 
      I1=IC 
      IC=K2 
      K2=I2 
      I2=IC 
      IF(I1.EQ.0) GO TO 19
      IF(I1.EQ.1) GO TO 21
   18 IF(K2.GE.I1) GO TO 22 
      IF(K2.LT.0) GO TO 8 
      IC=K2 
      K2=I1 
      I1=IC 
      IC=K1 
      K1=I2 
      I2=IC 
      IF (I1.EQ.0) GO TO 23 
      IF(I1.EQ.1) GO TO 25
   22 IF(K3.GE.I1) GO TO 26 
      IF(K3.LT.0) GO TO 8 
      IC=K3 
      K3=I1 
      I1=IC 
      IC=K1 
      K1=I3 
      I3=IC 
      IF(I1.EQ.0) GO TO 27
      IF(I1.EQ.1) GO TO 28
   26 M1(4)=I3
      M1(1)=I3
      M1(3)=K3
      M1(2)=K3
      M2(2)=I1
      M2(1)=I1
      M2(4)=K1
      M2(3)=K1
      M3(3)=I2
      M3(1)=I2
      M3(4)=K2
      M3(2)=K2
      M(1)=I1+I2+I3 
      M(2)=I1+K2+K3 
      M(3)=K1+I2+K3 
      M(4)=K1+K2+I3 
      Q1=1 
      DO 11 I=1,4 
      M(I)=M(I)/2 
   11 Q1=FT(M(I)-M1(I))*FT(M(I)-M2(I))*FT(M(I)-M3(I))*Q1/FT(M(I)+1) 
      Q1=SQRT(Q1) 
      M1(1)=I1+K1 
      M1(2)=I2+K2 
      M1(3)=I3+K3 
      IC=M1(1)+M1(2)  
      M(5)=IC/2 
      IC=M1(2)+M1(3)  
      M(6)=IC/2 
      IC=M1(1)+M1(3)  
      M(7)=IC/2 
      MAXZ=MIN(M(5),M(6),M(7)) 
      MINZ=MAX(M(1),M(2),M(3),M(4))
      X=0  
      DO 15 I=MINZ,MAXZ 
      Q2=1 
      DO 20 J=1,7 
      IJ=I-M(J) 
      IF(J.GT.4) IJ=-IJ 
   20 Q2=Q2*FT(IJ)
      Q2=FT(I+1)/Q2 
   15 X=-X+Q2 
      Q=X*Q1
      IS=MAXZ+IS
      IF(MOD(IS,2).EQ.1) Q=-Q 
      RETURN
    8 PRINT 1010,J1,J2,J3,L1,L2,L3
      ICE=.FALSE. 
 1010 FORMAT(10H ERREUR 6J,2(3X,3I3)) 
      RETURN
      ENTRY RACAH(J1,J2,J3,L1,L2,L3,Q)
      IS=(J1+J2+J3+L1)/2
      I1=J1 
      I2=J2 
      I3=L2 
      K1=L1 
      K2=J3 
      K3=L3 
      GO TO 24
      END 
*NEUFJ
      SUBROUTINE NEUFJ (J1,J2,J12,J3,J4,J34,J13,J24,J,Q)
      IMPLICIT REAL*8 (A-H,O-Z) 
      LOGICAL ICE 
      COMMON/TEUFJ/ICE
      ENTRYNEUFJI(J1,J2,J12,J3,J4,J34,J13,J24,J,Q)  
      I1=ABS(J1-J) 
      I2=ABS(J34-J2) 
      I3=ABS(J3-J24) 
      MING=MAX(I1,I2,I3)+2 
      I1=J1+J 
      I2=J34+J2 
      I3=J3+J24 
      MAXG=MIN(I1,I2,I3)+2 
      Q=0  
      IF(MING.GT.MAXG)RETURN
      ICE=.TRUE.
      DO 1 I=MING,MAXG,2
      IG=I-2
      CALL SIXJ(J1,J2,J12,J34,J,IG,Q1)
      CALL SIXJ(J3,J4,J34,J2,IG,J24,Q2)
      CALL SIXJ(J13,J24,J,IG,J1,J3,Q3)
      Q1=Q1*Q2*Q3*(IG+1)
    1 Q=Q+Q1
      IF(MOD(MAXG,2).EQ.1) Q=-Q 
      IF(ICE)   RETURN
      PRINT 10,J1,J2,J12,J3,J4,J34,J13,J24,J
   10 FORMAT(10H ERREUR 9J,3(3X,3I3)) 
      RETURN
      END 
*ERF
      FUNCTION ERF(X)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION C(18),D(17)
      DATA NCFC,NCFD/18,17/,XUP/6.25D0/,SQRTPI/1.7724538509055160D0/
     A,C(1),C(2),C(3),C(4),C(5),C(6),C(7),C(8),C(9),C(10),C(11)
     A,C(12),C(13),C(14),C(15),C(16),C(17),C(18)
     A/1.9449071068178803D0,4.20186582324414D-2,-1.86866103976769D-2
     A,5.1281061839107D-3,-1.0683107461726D-3,1.744737872522D-4
     A,-2.15642065714D-5,1.7282657974D-6,-2.00479241D-8
     A,-1.64782105D-8,2.0008475D-9,2.57716D-11,-3.06343D-11
     A,1.9158D-12,3.703D-13,-5.43D-14,-4.0D-15,1.2D-15/
     A,D(1),D(2),D(3),D(4),D(5),D(6),D(7),D(8),D(9),D(10),D(11)
     A,D(12),D(13),D(14),D(15),D(16),D(17)
     A/1.4831105640848036D0,-3.010710733865950D-1,6.89948306898316D-2
     A,-1.39162712647222D-2,2.4207995224335D-3,-3.658639685849D-4
     A,4.86209844323D-5,-5.7492565580D-6,6.113243578D-7
     A,-5.89910153D-8,5.2070091D-9,-4.232976D-10,3.18811D-11
     A,-2.2361D-12,1.467D-13,-9.0D-15,5.0D-16/
C
      DATA                             ZERO, ONE, TWO, THREE, TWENTY,
     *                                 HALF/0.0D0, 1.0D0, 2.0D0, 3.0D0,
     *                                 20.0D0, 0.5D0/
C     .. Executable Statements ..
C
C     NO FAILURE EXITS
      IFAIL = 0
      XV = ABS(X)
      IF (XV.GE.XUP) GO TO 120
      IF (XV.LE.TWO) GO TO 60
      X2 = TWO - TWENTY/(XV+THREE)
C
C     SUMMATION
      BJP2 = ZERO
      BJP1 = C(NCFC)
      J = NCFC - 1
   20 BJ = X2*BJP1 - BJP2 + C(J)
      IF (J.EQ.1) GO TO 40
      BJP2 = BJP1
      BJP1 = BJ
      J = J - 1
      GO TO 20
   40 X2 = HALF*(BJ-BJP2)/XV*EXP(-X*X)/SQRTPI
      ERF = (ONE-X2)*SIGN(ONE,X)
      GO TO 140
C
   60 X2 = X*X - TWO
C     SUMMATION
      BJP2 = ZERO
      BJP1 = D(NCFD)
      J = NCFD - 1
   80 BJ = X2*BJP1 - BJP2 + D(J)
      IF (J.EQ.1) GO TO 100
      BJP2 = BJP1
      BJP1 = BJ
      J = J - 1
      GO TO 80
  100 ERF = HALF*(BJ-BJP2)*X
      GO TO 140
C
  120 ERF = SIGN(ONE,X)
  140 RETURN
      END
*HARSPH 
      SUBROUTINE HARSPH (X,NN,P,Y)
      IMPLICIT REAL*8 (A-H,O-Z) 
      DIMENSIONP(2),Y(2),S(0:1000)
      SAVE  
      DATANMEM,S(0),S(1),S(2)/1,0.0d0,1.0d0,1.0d0/ 
      DATAC/0.28209479177388D0/ 
      N=ABS(NN) 
      IF(N+N-1.LE.NMEM)GOTO1
      N2=N+N-1
      NMEM=NMEM+1 
      DO2L=NMEM,N2
      Z=L
    2 S(L)=SQRT(Z) 
      NMEM=N2 
    1 P(1)=1.0d0 
      P(2)=X
      IF(N.LT.3)GOTO3 
      DO4L=3,N
    4 P(L)=((L+L-3)*X*P(L-1)-(L-2)*P(L-2))/(L-1)
    3 IF(NN.GT.0)RETURN 
      Q=SQRT(1-X**2) 
      K=2 
      Y(1)=1 
      DO5L=2,N
      Y(K)=P(L) 
      Y(K+L-1)=Y(K-1)*Q*S(L+L-3)/S(L+L-2) 
    5 K=K+L 
      IF(N.LT.3)GOTO6 
      N2=N-1
      K=1 
      DO7M=2,N2 
      K=K+M 
    7 Y(K+M)=S(M+M-1)*X*Y(K)
      N2=N-2
      IF(N.LT.4)GOTO6 
      KK=1  
      DO8M=2,N2 
      KK=KK+M 
      K=KK  
      DO9L=M,N2 
      K=K+L 
    9 Y(K+L+1)=((L+L+1)*X*Y(K)-S(L+M-1)*S(L-M+1)*Y(K-L))
     2 /(S(L+M)*S(L-M+2)) 
    8 CONTINUE
    6 K=0 
      DO10L=1,N 
      Q=S(L+L-1)*C
      DO10M=1,L 
      K=K+1 
      Y(K)=Q*Y(K) 
   10 Q=-Q  
      RETURN
      END 
*HANSPH 
      SUBROUTINE HANSPH (XX,N,BP,BM)
      IMPLICIT REAL*8 (A-H,O-Z) 
      DIMENSIONBP(2),BM(2)
      NN=ABS(N)
      X=ABS(XX) 
      IF(X.EQ.0)GOTO8
      E=1
      IF(XX.GT.0)E=EXP(-X) 
      I=X 
      IF(NN.LE.(I/5+1))GOTO1
      I=MAX(NN,I)+20 
      K=I/2*2 
      FLI=K+K+1 
      LIM=MAX(20,NN)/2*2 
      LM=(K-1)/LIM+1  
      K=MOD(K-1,LIM)+1
      BP(2)=0
      BP(1)=1/E 
      DO3J=1,LM 
      BP(K+1)=BP(1) 
      BP(K+2)=BP(2) 
      I=K 
      DO4L=1,K
      BP(I)=FLI/X*BP(I+1)+BP(I+2) 
      I=I-1 
    4 FLI=FLI-2
    3 K=LIM 
      FLI=((BP(1)+BP(2))*X+BP(1))*E 
      DO5L=1,NN 
    5 BP(L)=BP(L)/FLI 
    6 IF(N.GE.0)RETURN
      BM(1)=E/X*1.5707963267949D0 
      BM(2)=BM(1)*(1+1/X) 
      FLI=1
      DO2L=3,NN 
      FLI=FLI+2
    2 BM(L)=BM(L-2)+BM(L-1)*FLI/X 
      RETURN
    1 FLI=.5d0/E
      E2=E/2 
      IF(E.EQ.1.d0)E2=EXP(-X-X)/2
      BP(1)=(FLI-E2)/X
      BP(2)=(FLI+E2-BP(1))/X
      FLI=1
      DO7I=3,NN 
      FLI=FLI+2
    7 BP(I)=BP(I-2)-BP(I-1)*FLI/X 
      GOTO6 
    8 DO9I=2,NN 
    9 BP(I)=0
      BP(1)=1
      IF(N.LT.0)PRINT1001 
      RETURN
 1001 FORMAT(37H ISPH(0) N'EXISTE PAS POUR N NEGATIF )
      END 
*COUFRA 
      SUBROUTINE COUFRA(RHO,ETA,MINL,MAXL,FC,FCP,GC,GCP)
C*** FONCTIONS COULOMBIENNES CALCULEES EN R = RHO PAR LA METHODE DES FRA
C*** CONTINUES DE STEED. MINL ET MAXL CORRESPONDENT AUX VRAIES VALEURS D
C*** VOIR BARNETT, FENG, STEED ET GOLDFARB, CPC 1974 *******************
      IMPLICIT REAL*8 (A-H,O-Z) 
      REAL*8 K,K1,K2,K3,K4,M1,M2,M3,M4
      DIMENSION FC(MAXL),FCP(MAXL),GC(MAXL),GCP(MAXL)
      SAVE  
      DATA ACCUR,STEP/1.D-7,100.0D0/
      PACE = STEP 
      ACC = ACCUR 
      R = RHO 
      KTR = 1 
      LMAX = MAXL 
      LMIN1 = MINL+1  
      XLL1 = MINL*LMIN1 
      ETA2 = ETA**2 
      TURN = ETA+SQRT(ETA2+XLL1)
      IF(R.LT.TURN.AND.ABS(ETA).GE.1.D-6) KTR = -1  
      KTRP = KTR
      GO TO 2 
    1 R = TURN
      TF = F
      TFP = FP
      LMAX = MINL 
      KTRP = 1
    2 ETAR = ETA*R
   21 RHO2=R*R
      PL = LMAX+1 
      PMX = PL+0.5D0  
C** FRACTION CONTINUE POUR FP(MAXL)/F(MAXL) ; XL=F ; XLPRIME=FP ********
      FP = ETA/PL+PL/R
      DK = ETAR+ETAR  
      DEL = 0
      D = 0
      F = 1 
      K = (PL*PL-PL+ETAR)*(PL+PL-1) 
      IF(PL*PL+PL+ETAR.NE.0.) GO TO 3 
      R = R*1.0000001D0 
      GO TO 2 
    3 H = (PL*PL+ETA2)*(1-PL*PL)*RHO2 
      K = K+DK+PL*PL*6
      D = 1/(D*H+K) 
      DEL = DEL*(D*K-1) 
      IF(PL.LT.PMX) DEL = -R*(PL*PL+ETA2)*(PL+1)*D/PL 
      PL = PL+1 
      FP = FP+DEL 
      IF(D.LT.0) F = -F
      IF(PL.GT.20000.0D0) GO TO 11
      IF(ABS(DEL/FP).GE.ACC) GO TO 3
      FP = F*FP 
      IF(LMAX.EQ.MINL) GO TO 5  
      FC(LMAX+1) = F  
      FCP(LMAX+1) = FP
C*** RECURRENCE ARRIERE POUR F ET FP ; GC,GCP UTILISES POUR STOCKAGE ***
      L = LMAX
      DO 4 LP=LMIN1,LMAX
      PL = L
      GC(L+1) = ETA/PL+PL/R 
      GCP(L+1) = SQRT(ETA2+PL*PL)/PL
      FC(L) =(GC(L+1)*FC(L+1)+FCP(L+1))/GCP(L+1)
      FCP(L) = GC(L+1)*FC(L)-GCP(L+1)*FC(L+1) 
    4 L = L-1 
      F = FC(LMIN1) 
      FP = FCP(LMIN1) 
    5 IF(KTRP.EQ.-1) GO TO 1
C*** MEME CALCUL POUR R = TURN SI RHO.LT.TURN 
C*** P + I.Q CALCULE EN MINL , EQUATION (32)
      P = 0
      Q = R-ETA 
      PL = 0 
      AR = -(ETA2+XLL1) 
      AI = ETA
      BR = Q+Q
      BI = 2
      WI = ETA+ETA
      DR = BR/(BR*BR+BI*BI) 
      DI = -BI/(BR*BR+BI*BI)
      DP = -(AR*DI+AI*DR) 
      DQ = AR*DR-AI*DI
    6 P = P+DP
      Q = Q+DQ
      PL = PL+2 
      AR = AR+PL
      AI = AI+WI
      BI = BI+2 
      D = AR*DR-AI*DI+BR
      DI = AI*DR+AR*DI+BI 
      T = 1/(D*D+DI*DI) 
      DR = T*D
      DI = -T*DI
      H = BR*DR-BI*DI-1
      K = BI*DR+BR*DI 
      T = DP*H-DQ*K 
      DQ = DP*K+DQ*H  
      DP = T
      IF(PL.GT.46000.0D0) GO TO 11
      IF(ABS(DP)+ABS(DQ).GE.(ABS(P)+ABS(Q))*ACC) GO TO 6
      P = P/R 
      Q = Q/R 
C*** CALCUL DE FP,G,GP, ET NORMALISATION DE F EN L = MINL **************
      G = (FP-P*F)/Q  
      GP = P*G-Q*F
      W = 1/SQRT(FP*G-F*GP) 
      G = W*G 
      GP = W*GP 
      IF(KTR.EQ.1) GO TO 8
      F = TF
      FP = TFP
      LMAX = MAXL 
C*** CALCUL DE G(MINL) ET GP(MINL) PAR INTEGRATION RUNGE-KUTTA A PARTIR 
C***         VOIR FOX ET MAYERS(1968) PG 202
      IF(RHO.LT.0.2D0*TURN) PACE = 999.0D0
      R3=1.0D0/3.0D0  
      H = (RHO-TURN)/(PACE+1) 
      H2 = H/2
      I2 = INT(PACE+0.001D0)
      ETAH = ETA*H
      H2LL = H2*XLL1  
      S = (ETAH+H2LL/R)/R-H2
    7 RH2 = R+H2
      T = (ETAH+H2LL/RH2)/RH2-H2
      K1 = H2*GP
      M1 = S*G
      K2 = H2*(GP+M1) 
      M2 = T*(G+K1) 
      K3 = H*(GP+M2)  
      M3 = T*(G+K2) 
      M3 = M3+M3
      K4 = H2*(GP+M3) 
      RH = R+H
      S = (ETAH+H2LL/RH)/RH-H2  
      M4 = S*(G+K3) 
      G = G+(K1+K2+K2+K3+K4)*R3 
      GP = GP+(M1+M2+M2+M3+M4)*R3 
      R = RH
      I2 = I2-1 
      IF(ABS(GP).GT.1.D300) GO TO 11
      IF(I2.GE.0) GO TO 7 
      W = 1/(FP*G-F*GP) 
C*** RECURRENCE AVANT A PARTIR DE GC(MINL) ET GCP(MINL) 
C*** RENORMALISATION DE FC ET FCP POUR CHAQUE VALEUR DE L **************
    8 GC(LMIN1) = G 
      GCP(LMIN1) = GP 
      IF(LMAX.EQ.MINL) GO TO 10 
      DO 9 L=LMIN1,LMAX 
      T = GC(L+1) 
      GC(L+1) = (GC(L)*GC(L+1)-GCP(L))/GCP(L+1) 
      GCP(L+1) = GC(L)*GCP(L+1)-GC(L+1)*T 
      FC(L+1) = W*FC(L+1) 
    9 FCP(L+1) = W*FCP(L+1) 
      FC(LMIN1) = W*FC(LMIN1) 
      FCP(LMIN1) = W*FCP(LMIN1) 
      RETURN
   10 FC(LMIN1) = W*F 
      FCP(LMIN1) = W*FP 
      RETURN
   11 W = 0
      G = 0
      GP = 0 
      GO TO 8 
      END 
      SUBROUTINE WITT(AL,AM,X,W,WP,EPS)
      IMPLICIT REAL*8 (A-H,O-Z) 
      CALL WIT1(AL-1,AM,X,W1,EPS)
      CALL WIT1(AL-2,AM,X,W2,EPS)
      W=(X-2*(AL-1))*W1+(AM+1.5d0-AL)*(AM-1.5d0+AL)*W2
      WP=(AL/X-0.5d0)*W-(AM-AL+0.5d0)*(AM+AL-0.5d0)*W1/X
      RETURN
      END
      SUBROUTINE WIT1(AL,AM,X,W,EPS)
      IMPLICIT REAL*8 (A-H,O-Z) 
      AP=0.5d0+AM-AL
      BP=1+AM+AM
      W1=0
      DO1 L=1,7
      CALL WIT(L,AP,BP,X,W2)
      IF(ABS((W1-W2)/W2).LT.EPS)GOTO2
      W1=W2
    1 CONTINUE
      PRINT1000,AL,AM,X,EPS
    2 W=W2*X**AL*EXP(-X/2)/GAMA(AP)
      RETURN
 1000 FORMAT('PAS DE CONVERGENCE COWITT,AL,AM,X,EPS=',1P,4d13.5)
      END
      SUBROUTINE WIT(L,AP,BP,X,W)
      IMPLICIT REAL*8 (A-H,O-Z) 
      dimension KLA(8),AI(136),R(136)
      DATA KLA/ 0, 4,12,24,40,64,96,136/
      DATA (AI(I),I=  1,52 )/
     1  3.2254768961939D-01, 1.7457611011584D+00, 4.5366202969211D+00,
     1  9.3950709123010D+00,
     1  1.7027963230510D-01, 9.0370177679937D-01, 2.2510866298660D+00,
     1  4.2667001702877D+00, 7.0459054023936D+00, 1.0758516010181D+01,
     1  1.5740678641278D+01, 2.2863131736889D+01,
     1  1.1572211735804D-01, 6.1175748451513D-01, 1.5126102697764D+00,
     1  2.8337513377435D+00, 4.5992276394182D+00, 6.8445254531150D+00,
     1  9.6213168424567D+00, 1.3006054993306D+01, 1.7116855187462D+01,
     1  2.2151090379396D+01, 2.8487967250984D+01, 3.7099121044467D+01,
     1  8.7649410478933D-02, 4.6269632891507D-01, 1.1410577748312D+00,
     1  2.1292836450983D+00, 3.4370866338932D+00, 5.0780186145497D+00,
     1  7.0703385350481D+00, 9.4383143363917D+00, 1.2214223368866D+01,
     1  1.5441527368781D+01, 1.9180156856753D+01, 2.3515905693991D+01,
     1  2.8578729742881D+01, 3.4583398702284D+01, 4.1940452647686D+01,
     1  5.1701160339542D+01,
     1  5.9019852181480D-02, 3.1123914619844D-01, 7.6609690554590D-01,
     1  1.4255975908036D+00, 2.2925620586321D+00, 3.3707742642089D+00,
     1  4.6650837034671D+00, 6.1815351187366D+00, 7.9275392471720D+00,
     1  9.9120980150776D+00, 1.2146102711729D+01, 1.4642732289596D+01/
      DATA (AI(I),I= 53,108)/
     1  1.7417992646508D+01, 2.0491460082615D+01, 2.3887329848169D+01,
     1  2.7635937174332D+01, 3.1776041352373D+01, 3.6358405801649D+01,
     1  4.1451720484868D+01, 4.7153106445154D+01, 5.3608574544692D+01,
     1  6.1058531447214D+01, 6.9962240035101D+01, 8.1498279233948D+01,
     1  4.4489365833263D-02, 2.3452610951958D-01, 5.7688462930186D-01,
     1  1.0724487538177D+00, 1.7224087764446D+00, 2.5283367064257D+00,
     1  3.4922132730219D+00, 4.6164567697496D+00, 5.9039585041741D+00,
     1  7.3581267331861D+00, 8.9829409242124D+00, 1.0783018632539D+01,
     1  1.2763697986742D+01, 1.4931139755522D+01, 1.7292454336715D+01,
     1  1.9855860940335D+01, 2.2630889013196D+01, 2.5628636022458D+01,
     1  2.8862101816323D+01, 3.2346629153963D+01, 3.6100494805750D+01,
     1  4.0145719771537D+01, 4.4509207995752D+01, 4.9224394987306D+01,
     1  5.4333721333392D+01, 5.9892509162130D+01, 6.5975377287930D+01,
     1  7.2687628090655D+01, 8.0187446977905D+01, 8.8735340417883D+01,
     1  9.8829542868275D+01, 1.1175139809793D+02,
     1  3.5700394308758D-02, 1.8816228315858D-01, 4.6269428131451D-01,
     1  8.5977296397281D-01, 1.3800108205271D+00, 2.0242091359225D+00,
     1  2.7933693535066D+00, 3.6887026779080D+00, 4.7116411465546D+00,
     1  5.8638508783433D+00, 7.1472479081018D+00, 8.5640170175858D+00/
      DATA (AI(I),I=109,136)/
     1  1.0116634048452D+01, 1.1807892294004D+01, 1.3640933712536D+01,
     1  1.5619285893338D+01, 1.7746905950095D+01, 2.0028232834574D+01,
     1  2.2468249983498D+01, 2.5072560772425D+01, 2.7847480009168D+01,
     1  3.0800145739444D+01, 3.3938657084912D+01, 3.7272245880475D+01,
     1  4.0811492823884D+01, 4.4568603175330D+01, 4.8557763533054D+01,
     1  5.2795611187213D+01, 5.7301863323389D+01, 6.2100179072770D+01,
     1  6.7219370927125D+01, 7.2695158847608D+01, 7.8572802911565D+01,
     1  8.4911231135695D+01, 9.1789874671229D+01, 9.9320808717441D+01,
     1  1.0767244063938D+02, 1.1712230951268D+02, 1.2820184198825D+02,
     1  1.4228004446916D+02/
      DATA (R(I),I=  1,52 )/
     1  8.3273912383786D-01, 2.0481024384543D+00, 3.6311463058215D+00,
     1  6.4871450844076D+00,
     1  4.3772341049290D-01, 1.0338693476656D+00, 1.6697097656587D+00,
     1  2.3769247017586D+00, 3.2085409133480D+00, 4.2685755108251D+00,
     1  5.8180833686718D+00, 8.9062262152922D+00,
     1  2.9720963604447D-01, 6.9646298043062D-01, 1.1077813946158D+00,
     1  1.5384642390428D+00, 1.9983276062742D+00, 2.5007457691008D+00,
     1  3.0653215182824D+00, 3.7232891107827D+00, 4.5298140299818D+00,
     1  5.5972584618353D+00, 7.2129954609262D+00, 1.0543837461910D+01,
     1  2.2503631486429D-01, 5.2583605276235D-01, 8.3196139168709D-01,
     1  1.1460992409637D+00, 1.4717513169668D+00, 1.8131346873813D+00,
     1  2.1755175196946D+00, 2.5657627501650D+00, 2.9932150863713D+00,
     1  3.4712344831021D+00, 4.0200440864446D+00, 4.6725166077329D+00,
     1  5.4874206579862D+00, 6.5853612332889D+00, 8.2763579843641D+00,
     1  1.1824277551659D+01,
     1  1.5149441285947D-01, 3.5325658252990D-01, 5.5678456328812D-01,
     1  7.6268531769732D-01, 9.7187263224651D-01, 1.1853578930378D+00,
     1  1.4042656272844D+00, 1.6298686157570D+00, 1.8636350553320D+00,
     1  2.1072911510814D+00, 2.3629058910419D+00, 2.6330087531638D+00/
      DATA (R(I),I= 53,108)/
     1  2.9207575797277D+00, 3.2301851334923D+00, 3.5665733773686D+00,
     1  3.9370437554551D+00, 4.3515311888634D+00, 4.8244818548978D+00,
     1  5.3780220797891D+00, 6.0484178126197D+00, 6.9008983521807D+00,
     1  8.0699651561471D+00, 9.9027933194843D+00, 1.3820532094793D+01,
     1  1.1418710576812D-01, 2.6606521689761D-01, 4.1879313732490D-01,
     1  5.7253284649977D-01, 7.2764878838098D-01, 8.8453671934022D-01,
     1  1.0436188758921D+00, 1.2053492741523D+00, 1.3702213385218D+00,
     1  1.5387772564687D+00, 1.7116193526864D+00, 1.8894240634494D+00,
     1  2.0729593402465D+00, 2.2631066339970D+00, 2.4608890724883D+00,
     1  2.6675081263971D+00, 2.8843920929220D+00, 3.1132613270395D+00,
     1  3.3562176925959D+00, 3.6158698564843D+00, 3.8955130449484D+00,
     1  4.1993941047115D+00, 4.5331149785341D+00, 4.9042702876116D+00,
     1  5.3235009720236D+00, 5.8063332142336D+00, 6.3766146741598D+00,
     1  7.0735265807070D+00, 7.9676935092955D+00, 9.2050403312780D+00,
     1  1.1163013090768D+01, 1.5390180415262D+01,
     1  9.1625471157176D-02, 2.1342058490497D-01, 3.3571811668027D-01,
     1  4.5854093503343D-01, 5.8206816577905D-01, 7.0649521636717D-01,
     1  8.3202690300345D-01, 9.5887819879439D-01, 1.0872761620305D+00,
     1  1.2174623279777D+00, 1.3496954913567D+00, 1.4842549297768D+00/
      DATA (R(I),I=109,136)/
     1  1.6214441628118D+00, 1.7615953746767D+00, 1.9050746658947D+00,
     1  2.0522883472617D+00, 2.2036905532450D+00, 2.3597925385232D+00,
     1  2.5211741403764D+00, 2.6884980554088D+00, 2.8625278132105D+00,
     1  3.0441506653115D+00, 3.2344070972636D+00, 3.4345293984277D+00,
     1  3.6459928249939D+00, 3.8705845972164D+00, 4.1104986804328D+00,
     1  4.3684687232542D+00, 4.6479589840745D+00, 4.9534461124097D+00,
     1  5.2908484059009D+00, 5.6682046090332D+00, 6.0967964147430D+00,
     1  6.5931088610398D+00, 7.1824959955372D+00, 7.9066663113853D+00,
     1  8.8408924928109D+00, 1.0140899265621D+01, 1.2210021299206D+01,
     1  1.6705520642027D+01/
      save
      NW1=KLA(L)+1
      NW2=KLA(L+1)
      W=0
      DO201 I=NW1,NW2
      XX=AI(I)
  201 W=W+R(I)*XX**(AP-1)*(1+XX/X)**(BP-AP-1)*EXP(-XX)
      RETURN
      END
*NUM1L
      SUBROUTINE NUM1L(N,H,E,S2,U,S,NO,EPS) 
CC  VERSION CORRIGEE LE 21 NOV 72 
C*****INTEGRATION DE L"EQUATION DE SCHROEDINGER PAR LA METHODE DE NUMERO
C*****POUR E NEGATIF  
C*****RECHERCHE DE L"ENERGIE PROPRE PAR LA METHODE DE RAPHSON-NEWTON
      IMPLICIT REAL*8 (A-H,O-Z) 
      DIMENSION U(1),S(1) 
      DATA RAP1,RAP2/0,0/ 
      H12=H*H/12 
C*****CONTROLE DES CONDITIONS ASYMPTOTIQUES 
      IF(E.GT.0) E=0
      DEI=0
      EPSS=.1d-10 
      IF(U(N-1).GT.EPSS) GO TO 10 
      DEI=U(N-1)-EPSS 
      DO 8 K=1,N
    8 U(K)=U(K)-DEI 
   10 U(N)=U(N-1) 
C*****CALCUL DU NOMBRE D"ETATS LIES PAR INTEGRATION A ENERGIE NULLE 
      S(1)=1.d-10 
      B0=0 
      AA=H12*U(1) 
      IF (S2) 16,18,16
   16 B0=-S(1)*AA 
   18 B1=S(1)*(1-AA) 
      DO 38 K=2,N 
      B2=12*S(K-1)-10*B1-B0 
      IF (ABS(B2).LT.1.d+10) GO TO 22 
      B2=B2*1.d-20
      B1=B1*1.d-20
   22 AA=H12*U(K) 
      S(K)=B2/(1-AA) 
      B0=B1 
   38 B1=B2 
      DO 42 K=5,N 
      N0=K  
      IF(U(K).LT.0) GO TO 44 
   42 CONTINUE
   44 NEL=0 
      DO 52 K=N0,N
      IF (S(K-1)*S(K)) 46,50,52 
   46 NEL=NEL+2 
      GO TO 52
   50 NEL=NEL+1 
   52 CONTINUE
      NEL=NEL/2 
      IF(NEL.GT.NO) GO TO 64
      IF(NEL.EQ.NO) GO TO 60
   62 NO=-1 
      RETURN
   60 RAP1=S(N-1)/S(N)
      RAP2=EXP(H*SQRT(U(N-1)-E))
      IF(RAP1.LT.RAP2) GO TO 62 
C*****CALCUL DE EMIN ET EMAX ENTRE LESQUELLES SE TROUVE L"ENERGIE PROPRE
   64 UMIN=U(1) 
      DO 70 K=2,N 
      IF(U(K).LT.UMIN) UMIN=U(K)
   70 CONTINUE
      EMIN=UMIN 
      EMAX=0 
C*****DEBUT DE LA RECHERCHE DE L"ENERGIE PROPRE DANS L"INTERVALLE MAXIMU
      TE=EMAX-EMIN
C*****REJET DE L"ENERGIE D"ESSAI E PROPOSEE SI ELLE EST A L"EXTERIEUR DE
C*****BORNES (EMIN,EMAX)
      IF((E.LT.EMIN).OR.(E.GT.EMAX)) E=EMIN+TE/2 
      E1=EMIN 
      E2=EMAX 
      J=2 
      I=1 
      GO TO 102 
C*****REDUCTION DES BORNES EMIN ET EMAX 
   90 EMIN=E1 
      EMAX=E2 
      TE=EMAX-EMIN
      J=2 
   98 I=1 
  100 E=EMIN+TE*I/J 
  102 DE=0 
  104 E=E+DE
      IF(E.GT.0) GO TO 204 
      S(N)=1.d-10 
      N1=N-1
      RAP2=EXP(H*SQRT((U(N-1)+U(N))/2-E))
      S(N1)=S(N)*RAP2 
      AA=H12*(U(N1)-E)
      B0=S(N)*(1-AA) 
      B1=S(N1)*(1-AA)
      N1=N-2
      DO 138 KAUX=1,N1
      K=N1-KAUX+1 
      B2=12*S(K+1)-10*B1-B0 
      AA=H12*(U(K)-E) 
      S(K)=B2/(1-AA) 
      B0=B1 
      B1=B2 
      IF(U(K).LT.E) GO TO 140 
  138 CONTINUE
  140 N1=K  
C*****NORMALISATION DE LA FONCTION D"ONDE A S(N1) 
      DO 146 KAUX=N1,N
      K=N-KAUX+N1 
  146 S(K)=S(K)/S(N1) 
C*****DEBUT DE L"INTEGRATION VERS L"EXTERIEUR JUSQU"A N1
      S(1)=1.d-10 
      B0=0 
      AA=H12*(U(1)-E) 
      IF(S2) 156,158,156
  156 B0=-S(1)*AA 
  158 B1=S(1)*(1-AA) 
      DO 170 K=2,N1 
      B2=12*S(K-1)-10*B1-B0 
      AA=H12*(U(K)-E) 
      S(K)=B2/(1-AA) 
      B0=B1 
  170 B1=B2 
C*****NORMALISATION DE LA FONCTION A S(N1)
      DO 174 K=1,N1 
  174 S(K)=S(K)/S(N1) 
C*****CALCUL DE LA CORECTION D"ENERGIE
      SOM=0
      DO 180 K=1,N
  180 SOM=SOM+S(K)*S(K) 
      DE=((-S(N1-1)+2-S(N1+1))/(H*H)+U(N1)-E)/SOM  
      IF(ABS(DE).GT.EPS) GO TO 104
C*****CALCUL DU NOMBRE DE NOEUDS DE L"ETAT PROPRE TROUVE
      DO 182 K=5,N
      IF(U(K).LT.E) GO TO 184 
  182 CONTINUE
  184 N0=K  
      NEL=0 
      DO 192 K=N0,N1  
      IF(S(K-1)*S(K)) 186,190,192 
  186 NEL=NEL+2 
      GO TO 192 
  190 NEL=NEL+1 
  192 CONTINUE
      NEL=NEL/2 
C*****L"ETAT PROPRE TROUVE EST-IL LE BON  
      IF(NEL-NO) 198,214,202
  198 IF(E.GT.E1) E1=E
      GO TO 204 
  202 IF(E.LT.E2) E2=E
  204 I=I+2 
      IF (I.LE.J)  GO TO 100
      J=2*J 
      IF(ABS(E1-EMIN).GT.EPS.OR.ABS(EMAX-E2).GT.EPS) GO TO 90 
      GO TO 98
C*****NORMALISATION DE LA FONCTION PROPRE 
  214 SOM=1/SQRT(SOM*H)
      DO 218 K=1,N
  218 S(K)=S(K)*SOM 
      E=E+DEI 
      RETURN
C*****DEBUT FORMATS 
 2000 FORMAT(/,35X,56HL"ETAT DEMANDE N"EST PAS LIE. RETOUR DE NUM1L AVEC
     1 NO=-1,/) 
C*****FIN FORMATS 
      END 
*GAMMA
      FUNCTION GAMA (XX)
C 
C     GAMMA FUNCTION  
C---                              BY CHEBYSHEV POLYNOMIALS
C     LOGARITHM OF GAMMA FUNCTION 
C---                                             - CDC 6400 - 
C 
      IMPLICIT REAL*8 (A-H,O-Z) 
      DIMENSION CT(20),B(5)
      SAVE  
      DATA B / 12,-360,1260,-1680,1188/
      DATA ALPID/ 0.91893853320467d0/ 
      DATA CK / 0.57721566490153d0/ 
      DATA PI /3.1415926535898d0/ 
C 
      DATA CT/1.8835711955910D0, 0.44153813248410D-2, 
     1        0.56850436815994  D-1,-0.42198353964186  D-2, 
     2        0.13268081812125  D-2,-0.18930245297989  D-3, 
     3        0.36069253274412  D-4,-0.60567619044609  D-5, 
     4        0.10558295463023  D-5,-0.18119673655426  D-6, 
     5        0.3117724964715   D-7,-0.535421963902    D-8, 
     6 0.91932755199D- 9,-0.15779412803D- 9, 0.2707980623 D-10, 
     7-0.464681865  D-11, 0.79733502   D-12,-0.13680782   D-12, 
     8 0.2347319    D-13,-0.402743     D-14/
C 
C *************************** GAMMA (X) *************************** 
C 
      ENTRY GAMMA(XX) 
      X=XX  
      LY   = 1
      IF (X.LE.0) GO TO 99 
  100 IF (X.LE.1.0d-6) GO TO 300
      IF (X.EQ.2) GO TO 10
      Y    = X - 1  
      LG    = 1 
      T     = Y 
      IF (Y) 1,10,2 
    1 IY    = ABS(Y)+1
      T     = Y + IY 
      LG    = 2 
      GO TO 3 
    2 IF (Y.LE.1) GO TO 3 
      IY    = Y 
      T     = Y - IY 
      LG    = 3 
    3 F     = 4*T - 2
      A2    = CT(20)  
      A1    = F * A2      + CT(19)
      DO 4 I=3,18,2 
      J     = 20 - I  
      A2    = F * A1 - A2 + CT(J+1) 
      A1    = F * A2 - A1 + CT(J) 
    4 CONTINUE
      A2    = F * A1 - A2 + CT(2) 
      A2    = F * A2 - A1 + CT(1) 
      A     = (A2-A1)/2     
      IF(LG-2)5,6,8 
    5 GAMA = A
      IF (LY.EQ.2) GO TO 200
      GOTO11
    6 D     = 1 
      DO 7 I=1,IY 
      D     = D * (T - (I-1))
    7 CONTINUE
      GAMA = A / D
      IF (LY.EQ.2) GO TO 200
      GOTO11
    8 D     = 1
      DO 9 I=1,IY 
      D     = D * (T + I)
    9 CONTINUE
      GAMA = A * D
      IF (LY.EQ.1) GOTO11 
  200 GAMA = LOG(GAMA)
      GOTO11
   10 GAMA = 1
      IF (LY.EQ.1) GOTO11 
      GAMA = 0
      GOTO11
C 
C ************************* LOG GAMMA (X) ************************* 
C 
      ENTRY ALGAMA(XX)
      X=XX  
      LY   = 2
      IF (X.LE.0) GO TO 99 
   20 IF (X.LT.15) GO TO 100 
      XL   = LOG(X) 
      GAMA = X*(XL-1)-XL/2+ALPID
      IF (X.GT.2**12) GOTO11
      XL   = X * X
      F    = X
      G    = 1/(X*B(1))
      DO 21 I=2,5 
      F    = XL * F 
      G    = G + 1/(F*B(I))
   21 CONTINUE
      GAMA = GAMA + G 
      GOTO11
  300 ZZ   = (CK * X + 1) * X 
      GAMA = 1/ZZ 
      IF(LY.EQ.2)GOTO200
      GOTO11
   99 X=-X+1 
      IF(X.EQ.INT(X))GOTO12
      IF(LY.EQ.2)GOTO20 
      GOTO100 
   11 IF(XX.GT.0)RETURN
      IF(LY.EQ.1)GAMA=PI/(GAMA*SIN(PI*X)) 
      IF(LY.EQ.2)GAMA=GAMA+LOG(PI/ABS(SIN(PI*X))) 
      RETURN
   12 PRINT1000,INT(XX) 
      GAMA=0 
      RETURN
 1000 FORMAT(8H GAMMA (,I4,15H ) N EXISTE PAS)
      END   
   

*AMBLE
      SUBROUTINE AMBLE(A,B,E,S,M,NM)
      IMPLICIT real*8 (A-H,O-Z)
      DIMENSIONA(NM,NM),B(NM,NM),S(NM,NM),E(NM,2) 
      LOGICALTEST 
      SAVE
      DATA NZ/1/
      N=ABS(M) 
      CALLREDUC(M,NM,A,B,E(1,2),TEST) 
      IF(.NOT.TEST)GOTO7
      PRINT1000 
      RETURN
    7 DO1I=1,N
      DO1J=1,N
    1 A(I,J)=A(J,I) 
      CALL MARDI(A,E,S,N,0,NM,1.d-14)
      DO3J=1,N
      S(N,J)=S(N,J)/E(N,2)
      IF(N.EQ.1)RETURN
      DO3L=2,N
      I=N+2-L 
      I1=I-1
      X=S(I1,J) 
      DO2K=I,N
    2 X=X-B(K,I1)*S(K,J)
    3 S(I1,J)=X/E(I1,2) 
      IF(NZ.EQ.0.OR.N.EQ.1) RETURN
      IM=N-1
      DO 4 I=1,IM 
      JM=I+1
      DO 5 J=JM,N 
      IF(E(I,1).LT.E(J,1)) GOTO 5 
      Z=E(J,1)
      E(J,1)=E(I,1) 
      E(I,1)=Z
      DO 6 K=1,N
      ZZ=S(K,J) 
      S(K,J)=S(K,I) 
      S(K,I)=ZZ 
    6 CONTINUE
    5 CONTINUE
    4 CONTINUE
      RETURN
 1000 FORMAT(40H LA MATRICE B N EST PAS DEFINIE POSITIVE) 
      END 
*MARDIF
      SUBROUTINE MARDI(A,E,S,NLIM,NT,N,RHO) 
      IMPLICIT real*8 (A-H,O-Z)
      DIMENSION A(1),S(1),E(1)  
      M = N*NLIM
      KII = NT*N
      DO 3 I=1,M
      KII = KII+1 
    3 S(KII) = 0
      K = N+1 
      NS = NT*K 
      IF(NLIM.GT.1) GO TO 20
      E(NT+1) = A(NS+1) 
      S(NS+1) = 1
      RETURN
   20 XNU = 0
      ALIM = NLIM 
      ALIM = 1/ALIM
      KIJ = N+NS
      KNS = NS-K
      NS = NS-N 
      DO 1 J=2,NLIM 
      DO 2 I=2,J
      KIJ = KIJ+1 
    2 XNU = XNU + A(KIJ)**2 
    1 KIJ = KIJ+K-J 
      XNU = SQRT(2*XNU)
      IF(RHO.LT.0) PRINT 1000,XNU
 1000 FORMAT(1H ,10X,17H(NORME INITIALE = ,1PE10.2,1H)) 
      XNUF = XNU* ABS(RHO)*ALIM 
      KII = NS
      DO 4 I=1,NLIM 
      KII = KII+K 
      J = I+NT
      E(J) = A(KII) 
    4 S(KII) = 1
    5 XNU = XNU*ALIM  
    6 M = 0 
      KJJ = NS+K
      DO 7 J=2,NLIM 
      KJJ = KJJ+K 
      KIJ = KJJ-J 
      DO 7 I=2,J
      KIJ = KIJ+1 
      AB = ABS(A(KIJ))
      IF(AB.LE.XNU) GO TO 7 
      M = I-1 
      KIL = I+KNS 
      KII = N*M+KIL 
      KLI = KII-M 
      KLJ = KJJ-J 
      KJL = J+NS
      AII = A(KII)
      AJJ = A(KJJ)
      AIJ = A(KIJ)
      B = AII-AJJ 
      IF(AB.LT.(.000175d0*ABS(B))) GOTO 9
      C = - SIGN(1.00d0,B)*AIJ/ SQRT(AIJ*AIJ+B*B/4)
      B = (SQRT(1-C*C)+1)*2
      SI = C/ SQRT(B) 
      S2 = SI*SI
      C2 = 1-S2
      C = SQRT(C2)
      GO TO 10
    9 C = -AIJ/B
      B = C*C/2
      SI = C*(1-3*B)
      C = 1-B
      S2 = SI*SI
      C2 = 1-S2
   10 ASC = 2*SI*C*AIJ
      DO 11 L=1,NLIM  
      KLI = KLI+1 
      KLJ = KLJ+1 
      ALI = S(KLI)
      S(KLI) = ALI*C-S(KLJ)*SI  
      S(KLJ) = ALI*SI+S(KLJ)*C  
      KIL = KIL+N 
      KJL = KJL+N 
      KL1 = MAX(KIL,KLI) 
      KL2 = MAX(KJL,KLJ) 
      ALI = A(KL1)
      A(KL1) = ALI*C-A(KL2)*SI  
   11 A(KL2) = ALI*SI+A(KL2)*C  
      A(KII) = AII*C2+AJJ*S2-ASC
      A(KJJ) = AII*S2+AJJ*C2+ASC
      A(KIJ) = 0
    7 CONTINUE
      IF(M.NE.0) GO TO 6
      IF(XNU.GT.XNUF) GO TO 5 
      KII = NS
      DO 12 I=1,NLIM  
      J = I+NT
      KII = KII+K 
      AII = A(KII)
      A(KII) = E(J) 
   12 E(J) = AII
      KIJ = NS+2*N
      DO 13 J=2,NLIM  
      KJI = J+NS
      DO 14 I=2,J 
      KIJ = KIJ+1 
      KJI = KJI+N 
   14 A(KIJ) = A(KJI) 
   13 KIJ = KIJ+K-J 
      call order(e,s,nlim,n) !ajout le 25/10/07
      RETURN
      END 
      SUBROUTINE REDUC(N,NM,A,B,DL,FAIL)
      IMPLICIT real*8(A-H,O-Z)
      DIMENSIONA(NM,NM),B(NM,NM),DL(NM) 
      LOGICALFAIL 
      FAIL=.TRUE. 
      IF(N.LT.0)GOTO1 
      DO2I=1,N
      DO2J=I,N
      X=B(I,J)
      I1=I-1
      IF(I1.EQ.0)GOTO10 
      DO3K=1,I1 
    3 X=X-B(I,K)*B(J,K) 
   10 IF(I.NE.J)GOTO4 
      FAIL=X.LE.0.d0
      IF(FAIL)RETURN  
      DL(I)=SQRT(X) 
      Y=DL(I) 
      GOTO2 
    4 B(J,I)=X/Y
    2 CONTINUE
    1 DO5I=1,N
      I1=I-1
      Y=DL(I) 
      DO5J=I,N
      X=A(I,J)
      IF(I1.EQ.0)GOTO5
      DO6K=1,I1 
    6 X=X-B(I,K)*A(J,K) 
    5 A(J,I)=X/Y
      DO7J=1,N
      J1=J-1
      DO7I=J,N
      I1=I-1
      X=A(I,J)
      IF(I.EQ.J)GOTO11
      DO8K=J,I1 
    8 X=X-A(K,J)*B(I,K) 
   11 IF(J.EQ.1)GOTO7 
      DO9K=1,J1 
    9 X=X-A(J,K)*B(I,K) 
    7 A(I,J)=X/DL(I)  
      RETURN
      END 
*BESSPH 
      SUBROUTINEBESSPH(X,N,BP,BM) 
      IMPLICIT REAL*8 (A-H,O-Z) 
      DIMENSIONBP(N),BM(N)
      SAVE
      NN=ABS(N)
      XX=ABS(X) 
      IF(NN.EQ.0)NN=1 
      IF(X.EQ.0)GOTO8
      I=X 
      IF(NN.LE.I)GOTO1
      I=MAX(NN+10,I+25)
   10 K=I/2*2 
      FLI=K+K+1 
      LIM=MAX(20,NN)/2*2 
      LM=(K-1)/LIM+1  
      K=MOD(K-1,LIM)+1
      BP(2)=0
      BP(1)=1.d-293 
      DO13J=1,LM
      BP(K+1)=BP(1) 
      BP(K+2)=BP(2) 
      I=K 
      DO4L=1,K
      BP(I)=FLI/X*BP(I+1)-BP(I+2) 
      I=I-1 
    4 FLI=FLI-2
   13 K=LIM 
      BM1=BP(2)-BP(1)/X 
      IF(ABS(BP(1)).LT.ABS(BM1))GOTO11  
      FLI=BP(1)/SIN(X)*X
      GOTO12
   11 FLI=-BM1/COS(X)*X 
   12 DO5L=1,NN 
    5 BP(L)=BP(L)/FLI 
      BM1=BM1/FLI 
    6 IF(N.GT.0)RETURN
      BM(1)=BM1 
      BM(2)=BM(1)/X-BP(1) 
      IF(NN.LE.2)RETURN 
      FLI=1
      DO7I=3,NN 
      FLI=FLI+2
    7 BM(I)=FLI/X*BM(I-1)-BM(I-2) 
      RETURN
    1 BP(1)=SIN(X)/X  
      BM1=-COS(X)/X 
      BP(2)=BP(1)/X+BM1 
      FLI=1
      DO2I=3,NN 
      FLI=FLI+2
    2 BP(I)=FLI/X*BP(I-1)-BP(I-2) 
      GOTO6 
    8 DO9I=2,NN 
    9 BP(I)=0
      BP(1)=1
      IF(N.LT.0)PRINT1001 
      RETURN
 1001 FORMAT(37H JSPH(0) N'EXISTE PAS POUR N NEGATIF )
      END 
*MATIN
      SUBROUTINE MATIN (A,IDIM1,N1,N2,INDEX,NERROR,DETERM)
      IMPLICIT REAL*8(A-H,O-Z)
C
C         MATRIX INVERSION WITH ACCOMPANYING SOLUTION OF LINEAR EQUATION
C
      DIMENSION A(IDIM1*IDIM1),INDEX(IDIM1)
      DETER=1
      N=N1
      IEMAT=N+N2
      IDIM=IDIM1
      NMIN1=N-1
C        THE ROUTINE DOES ITS OWN EVALUATION FOR DOUBLE SUBSCRIPTING OF
C        ARRAY A.
      IPIVC=1-IDIM
C        MAIN LOOP TO INVERT THE MATRIX
      DO 11 MAIN=1,N
      PIVOT=0
      IPIVC=IPIVC+IDIM
C        SEARCH FOR NEXT PIVOT IN COLUMN MAIN.
      IPIVC1=IPIVC+MAIN-1
      IPIVC2=IPIVC +NMIN1
      DO 2 I1=IPIVC1,IPIVC2
      IF(ABS(A(I1))-ABS(PIVOT)) 2,2,1
    1 PIVOT=A(I1)
      LPIV=I1
    2 CONTINUE
C        IS PIVOT DIFFERENT FROM ZERO
      IF(PIVOT) 3,15,3
C        GET THE PIVOT-LINE INDICATOR AND SWAP LINES IF NECESSARY
    3 ICOL=LPIV-IPIVC+1
      INDEX(MAIN)=ICOL
      IF(ICOL-MAIN) 6,6,4
C        COMPLEMENT THE DETERMINANT
    4 DETER=-DETER
C        POINTER TO LINE PIVOT FOUND
      ICOL=ICOL-IDIM
C        POINTER TO EXACT PIVOT LINE
      I3=MAIN-IDIM
      DO 5 I=1,IEMAT
      ICOL=ICOL+IDIM
      I3=I3+IDIM
      SWAP=A(I3)
      A(I3)=A(ICOL)
    5 A(ICOL)=SWAP
C        COMPUTE DETERMINANT
    6 DETER=DETER*PIVOT
      PIVOT=1/PIVOT
C        TRANSFORM PIVOT COLUMN
      I3=IPIVC+NMIN1
      DO 7 I=IPIVC,I3
    7 A(I)=-A(I)*PIVOT
      A(IPIVC1)=PIVOT
C        PIVOT ELEMENT TRANSFORMED
C
C        NOW CONVERT REST OF THE MATRIX
      I1=MAIN-IDIM
C        POINTER TO PIVOT LINE ELEMENTS
      ICOL=1-IDIM
C        GENERAL COLUMN POINTER
      DO 10 I=1,IEMAT
      ICOL=ICOL+IDIM
      I1=I1+IDIM
C        POINTERS MOVED
      IF(I-MAIN) 8,10,8
C        PIVOT COLUMN EXCLUDED
    8 JCOL=ICOL+NMIN1
      SWAP=A(I1)
      I3=IPIVC-1
      DO 9 I2=ICOL,JCOL
      I3=I3+1
    9 A(I2)=A(I2)+SWAP*A(I3)
      A(I1)=SWAP*PIVOT
   10 CONTINUE
   11 CONTINUE
C        NOW REARRANGE THE MATRIX TO GET RIGHT INVERS
      DO 14 I1=1,N
      MAIN=N+1-I1
      LPIV=INDEX(MAIN)
      IF(LPIV-MAIN) 12,14,12
   12 ICOL=(LPIV-1)*IDIM+1
      JCOL=ICOL+NMIN1
      IPIVC=(MAIN-1)*IDIM+1-ICOL
      DO 13 I2=ICOL,JCOL
      I3=I2+IPIVC
      SWAP=A(I2)
      A(I2)=A(I3)
   13 A(I3)=SWAP
   14 CONTINUE
      DETERM=DETER
      NERROR=0
      RETURN
   15 NERROR=MAIN
      DETERM=DETER
      RETURN
      END
      SUBROUTINE INSPL5 (X,CF,H,N)
C 
C  J.-P. HENNART,P. HALLET AND E. MUND
C    U.L.B. COMPUTER CENTER - METROLOGIE NUCLEAIRE  
C    AVENUE AD. BUYL  91
C    1150  BRUSSELS 
C 
C  CAT. A - INTERPOLATION - APPROXIMATION 
C 
C 
C  CALCULATION OF THE PARAMETERS OF A NATURAL QUINTIC SPLINE  
C  OF INTERPOLATION 
C 
C 
C  N            NUMBER OF INTERPOLATION POINTS
C 
C  X  (N)       VECTOR OF ABSCISSAE 
C  CF(N,1)      ORDINATES 
C  CF(N,2)      FIRST DERIVATIVES 
C  CF(N,3)      SECOND DERIVATIVES
C  H(6*N-3)     AUXILIARY VECTOR
C 
      IMPLICIT REAL*8 (A-H,O-Z) 
      DIMENSION X(N),CF(N,3),H(6*N-3) 
      J=N 
      DO 90 I=2,J 
        K=I-1 
        IF(X(I).GT.X(K))GOTO90  
        PRINT1000,K,X(K),I,X(I) 
        STOP
   90   CONTINUE
      IA=4*J-4
      DP1=0
      DP2=0
      DP3=0
      SF1=0
      SF2=0
C 
C  COMPUTATION OF MATRIX ELEMENTS AND SIMULTANEOUS  
C  FORWARD ELIMINATION
C 
      DO 60 I=1,J 
        IB=4*I-3
        ID=IA+I+I-1 
        IF (I.LE.1) GO TO 10
        IC=IB-4 
        IE=ID-2 
        DI1=-DS1
        DI2=+DS2
        DI3=+DS3
        DI4=-DS4
        DP1=- 12*P2  
        DP2=- 64*P3  
        DP3=+  3*P1  
        DP4=-DP1
        SF1=FP3 
        SF2=FP4 
   10   IF (I.GE.J) GO TO 20
        P1=1/(X(I+1)-X(I)) 
        P2=P1*P1
        P3=P2*P1
        FP =CF(I,1)-CF(I+1,1) 
        FP3=20*P3*FP 
        FP4= 6*P1*FP3
        DP1=DP1+ 12*P2 
        DP2=DP2- 64*P3 
        DP3=DP3+  3*P1 
        DP4=-DP1
        DS1=-  8*P2  
        DS2=+ 56*P3  
        DS3=+     P1  
        DS4=DS1 
        SF1=SF1-FP3 
        SF2=SF2+FP4 
   20   T1=DP1
        T2=DP2
        T3=DP3
        T4=DP4
        IF (I.LE.1) GO TO 30
        T1=T1-DI1*H(IC+1)-DI3*H(IC+2) 
        T2=T2-DI2*H(IC+1)-DI4*H(IC+2) 
        T3=T3-DI1*H(IC+3)-DI3*H(IC+4) 
        T4=T4-DI2*H(IC+3)-DI4*H(IC+4) 
   30   DET=T1*T4-T2*T3 
        IF (I.GE.J) GO TO 40
        H(IB+1)=(+T4*DS1-T3*DS2)/DET
        H(IB+2)=(-T2*DS1+T1*DS2)/DET
        H(IB+3)=(+T4*DS3-T3*DS4)/DET
        H(IB+4)=(-T2*DS3+T1*DS4)/DET
   40   U1=SF1
        U2=SF2
        IF (I.LE.1) GO TO 50
        U1=U1+DI1*H(IE+1)+DI3*H(IE+2) 
        U2=U2+DI2*H(IE+1)+DI4*H(IE+2) 
   50   H(ID+1)=(+T4*U1-T3*U2)/DET
        H(ID+2)=(-T2*U1+T1*U2)/DET
   60 CONTINUE
C 
C  BACKWARD SUBSTITUTION
C 
      DO 80 I=1,J 
        K=J+1-I 
        IF (I.LE.1) GOTO 70 
        IB=4*K-3
   70   ID=IA+K+K-1 
        CF(K,2)=H(ID+1) 
        CF(K,3)=H(ID+2) 
        IF (I.LE.1) GO TO 80
        CF(K,2)=CF(K,2)+H(IB+1)*CF(K+1,2)+H(IB+3)*CF(K+1,3) 
        CF(K,3)=CF(K,3)+H(IB+2)*CF(K+1,2)+H(IB+4)*CF(K+1,3) 
   80 CONTINUE
      RETURN
 1000 FORMAT(45H    ORDRE DES ABSCISSES INCORRECT DANS SPLIN5/
     1 8H      X(,I3,4H ) =,G9.4,8H  GE  X(,I3,4H ) =,G9.4) 
      END 
*SPLIN5 
      SUBROUTINE SPLIN5 (XX,F,X,CF,N,NF)  
C MOD. 12/2/76
C 
C  J.-P. HENNART,P. HALLET AND E. MUND
C    U.L.B. COMPUTER CENTER - METROLOGIE NUCLEAIRE  
C    AVENUE AD. BUYL  91
C    1150  BRUSSELS 
C 
C  CAT. A - INTERPOLATION - APPROXIMATION 
C 
C 
C  CALCULATION OF THE VALUES OF A NATURAL QUINTIC SPLINE OF 
C  INTERPOLATION AND ( OR ) ITS FIRST AND SECOND DERIVATIVES  
C  AT ABSCISSA XX 
C 
C 
C  N            NUMBER OF INTERPOLATION POINTS
C 
C  X  (N)       VECTOR OF ABSCISSAE 
C  CF(N,1)      ORDINATES 
C  CF(N,2)      FIRST DERIVATIVES 
C  CF(N,3)      SECOND DERIVATIVES
C 
C 
C  XX           ABSCISSA WHERE VALUES OF THE FUNCTION AND 
C               (OR) ITS FIRST AND SECOND DERIVATIVES ARE 
C               COMPUTED FOLLOWING INDEX NF VALUE 
C 
C  NF=1         THE FUNCTION IS CALCULATED
C  NF=2         THE FUNCTION AND ITS FIRST DERIVATIVE ARE 
C               CALCULATED
C  NF=3         THE FUNCTION AND ITS FIRST AND SECOND 
C               DERIVATIVES ARE CALCULATED
C 
C  F(1)         VALUE OF THE FUNCTION 
C  F(2)         VALUE OF THE FIRST DERIVATIVE 
C  F(3)         VALUE OF THE SECOND DERIVATIVE
C 
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION F(3),X(N),CF(N,3)
      SAVE
      DATA I/2/
      IF(I.GT.N)I=2 
   12 XN = XX-X(I)
      IF(XN)11,40,20  
   11 XN = XX-X(I-1)  
      IF(XN)70,51,51  
   50 XN = XX-X(I-1)  
   51 H=X(I)-X(I-1) 
      CC=CF(I-1,1)-CF(I,1)
      XR1=XN/H
      XR2=XR1*XR1 
      XR3=XR2*XR1 
      Y=CF(I-1,1)+CC*XR3*(-10+XR1*(+15- 6*XR1))  
      W=CF(I-1,2)*XR1 
     1    +CF(I-1,2)*XR3*(-6+XR1*(8-3*XR1))  
     2    -CF(I,2  )*XR3*(4+XR1*(-7+3*XR1))  
      U=CF(I-1,3)*XR2/2
     1    +CF(I-1,3)*XR3*(-1.5d0+XR1*(+1.5d0- .5d0*XR1))  
     2    +CF(I,3  )*XR3*(+ .5d0+XR1*(-1+ .5d0*XR1))  
      F(1)=Y+H*(W+H*U)
      IF (NF.EQ.1) GO TO 60 
      Y=      30*CC*XR2*(-1+XR1*(2-XR1))  
      W=CF(I-1,2) 
     1    +CF(I-1,2)*XR2*(-18+XR1*(32-15*XR1))  
     2    -CF(I,2  )*XR2*(+12+XR1*(-28+15*XR1))  
      U=CF(I-1,3)*XR1 
     1    +CF(I-1,3)*XR2*(-4.5d0+XR1*(6-2.5d0*XR1))  
     2    +CF(I,3  )*XR2*(+1.5d0+XR1*(-4+2.5d0*XR1))  
      F(2)=Y/H+W+H*U  
      IF (NF.EQ.2) GO TO 60 
      Y=      60*CC*XR1*(- 1+XR1*(+ 3- 2*XR1))  
      W=   CF(I-1,2)*XR1*(-36+XR1*(+96-60*XR1))  
     1    -CF(I,2  )*XR1*(+24+XR1*(-84+60*XR1))  
      U=CF(I-1,3) 
     1    +CF(I-1,3)*XR1*(- 9+XR1*(+18-10*XR1))  
     2    +CF(I,3  )*XR1*(+ 3+XR1*(-12+10*XR1))  
      F(3)=(Y/H+W)/H+U
   60 RETURN
   40 F(1)=CF(I,1)
      IF (NF.EQ.1) GO TO 60 
      F(2)=CF(I,2)
      IF (NF.EQ.2) GO TO 60 
      F(3)=CF(I,3)
      GO TO 60
   20 IF(I.EQ.N) GO TO 13 
      IMIN = I+1
      DO 30 I=IMIN,N  
        XN=XX-X(I)
        IF (XN) 50,40,30
   30 CONTINUE
      I = N 
   13 J = N 
   10 CONTINUE
      F(1) = CF(J,1)+XN*(CF(J,2)+XN*CF(J,3)/2)
      IF (NF.EQ.1) GO TO 60 
      F(2)=CF(J,2)+XN*CF(J,3) 
      IF (NF.EQ.2) GO TO 60 
      F(3)=CF(J,3)
      GO TO 60
   70 IF(I.EQ.2) GO TO 80 
      I = 2 
      GO TO 12
   80 J = 1 
      GO TO 10
C 
      END 
      SUBROUTINE INIGAUS
      implicit real*8(a-h,o-z)
      COMMON/GAU/KTAB(10),U(152),W(152)
      return
      end
      blockdata legendre
      implicit real*8(a-h,o-z)
      COMMON/GAU/KTAB(10),U(152),W(152)
      DATA KTAB/0,2,6,12,20,32,48,72,104,152/
      DATA (U(I) ,I=1,38)/
     1    3.3998104358485 D-01,    8.6113631159405 D-01,
     1    1.8343464249565 D-01,    5.2553240991633 D-01,
     1    7.9666647741362 D-01,    9.6028985649753 D-01,
     1    1.2523340851146 D-01,    3.6783149899818 D-01,
     1    5.8731795428661 D-01,    7.6990267419430 D-01,
     1    9.0411725637047 D-01,    9.8156063424671 D-01,
     1    9.5012509837637 D-02,    2.8160355077925 D-01,
     1    4.5801677765722 D-01,    6.1787624440264 D-01,
     1    7.5540440835500 D-01,    8.6563120238783 D-01,
     1    9.4457502307323 D-01,    9.8940093499165 D-01,
     1    6.4056892862605 D-02,    1.9111886747361 D-01,
     1    3.1504267969616 D-01,    4.3379350762604 D-01,
     1    5.4542147138883 D-01,    6.4809365193697 D-01,
     1    7.4012419157855 D-01,    8.2000198597390 D-01,
     1    8.8641552700440 D-01,    9.3827455200273 D-01,
     1    9.7472855597131 D-01,    9.9518721999702 D-01,
     1    4.8307665687737 D-02,    1.4447196158279 D-01,
     1    2.3928736225213 D-01,    3.3186860228212 D-01,
     1    4.2135127613063 D-01,    5.0689990893223 D-01/
      DATA (U(I),I=39,76)/
     1    5.8771575724076 D-01,    6.6304426693021 D-01,
     1    7.3218211874029 D-01,    7.9448379596794 D-01,
     1    8.4936761373256 D-01,    8.9632115576605 D-01,
     1    9.3490607593773 D-01,    9.6476225558750 D-01,
     1    9.8561151154526 D-01,    9.9726386184948 D-01,
     1    3.2380170962869 D-02,    9.7004699209462 D-02,
     1    1.6122235606889 D-01,    2.2476379039468 D-01,
     1    2.8736248735545 D-01,    3.4875588629216 D-01,
     1    4.0868648199071 D-01,    4.6690290475095 D-01,
     1    5.2316097472223 D-01,    5.7722472608397 D-01,
     1    6.2886739677651 D-01,    6.7787237963266 D-01,
     1    7.2403413092381 D-01,    7.6715903251573 D-01,
     1    8.0706620402944 D-01,    8.4358826162439 D-01,
     1    8.7657202027424 D-01,    9.0587913671556 D-01,
     1    9.3138669070655 D-01,    9.5298770316043 D-01,
     1    9.7059159254624 D-01,    9.8412458372282 D-01,
     1    9.9353017226635 D-01,    9.9877100725242 D-01,
     1    2.4350292663424 D-02,    7.2993121787798 D-02,
     1    1.2146281929612 D-01,    1.6964442042399 D-01/
      DATA (U(I),I=77,114)/
     1    2.1742364374000 D-01,    2.6468716220876 D-01,
     1    3.1132287199021 D-01,    3.5722015833766 D-01,
     1    4.0227015796399 D-01,    4.4636601725346 D-01,
     1    4.8940314570705 D-01,    5.3127946401989 D-01,
     1    5.7189564620263 D-01,    6.1115535517239 D-01,
     1    6.4896547125465 D-01,    6.8523631305423 D-01,
     1    7.1988185017161 D-01,    7.5281990726053 D-01,
     1    7.8397235894334 D-01,    8.1326531512279 D-01,
     1    8.4062929625258 D-01,    8.6599939815409 D-01,
     1    8.8931544599511 D-01,    9.1052213707850 D-01,
     1    9.2956917213193 D-01,    9.4641137485840 D-01,
     1    9.6100879965205 D-01,    9.7332682778991 D-01,
     1    9.8333625388462 D-01,    9.9101337147674 D-01,
     1    9.9634011677195 D-01,    9.9930504173577 D-01,
     1    1.6276744849602 D-02,    4.8812985136049 D-02,
     1    8.1297495464425 D-02,    1.1369585011066 D-01,
     1    1.4597371465489 D-01,    1.7809688236761 D-01,
     1    2.1003131046056 D-01,    2.4174315616384 D-01,
     1    2.7319881259104 D-01,    3.0436494435449 D-01/
      DATA (U(I),I=115,152)/
     1    3.3520852289262 D-01,    3.6569686147231 D-01,
     1    3.9579764982890 D-01,    4.2547898840729 D-01,
     1    4.5470942216774 D-01,    4.8345797392059 D-01,
     1    5.1169417715466 D-01,    5.3938810832435 D-01,
     1    5.6651041856139 D-01,    5.9303236477757 D-01,
     1    6.1892584012546 D-01,    6.4416340378496 D-01,
     1    6.6871831004391 D-01,    6.9256453664217 D-01,
     1    7.1567681234896 D-01,    7.3803064374439 D-01,
     1    7.5960234117664 D-01,    7.8036904386743 D-01,
     1    8.0030874413914 D-01,    8.1940031073793 D-01,
     1    8.3762351122818 D-01,    8.5495903343460 D-01,
     1    8.7138850590929 D-01,    8.8689451740242 D-01,
     1    9.0146063531585 D-01,    9.1507142312089 D-01,
     1    9.2771245672230 D-01,    9.3937033975275 D-01,
     1    9.5003271778443 D-01,    9.5968829144874 D-01,
     1    9.6832682846326 D-01,    9.7593917458513 D-01,
     1    9.8251726356301 D-01,    9.8805412632962 D-01,
     1    9.9254390032376 D-01,    9.9598184298720 D-01,
     1    9.9836437586317 D-01,    9.9968950388322 D-01/
      DATA (W(I),I=1,38)/
     1    6.5214515486254 D-01,    3.4785484513745 D-01,
     1    3.6268378337836 D-01,    3.1370664587788 D-01,
     1    2.2238103445337 D-01,    1.0122853629037 D-01,
     1    2.4914704581340 D-01,    2.3349253653835 D-01,
     1    2.0316742672306 D-01,    1.6007832854334 D-01,
     1    1.0693932599531 D-01,    4.7175336386511 D-02,
     1    1.8945061045506 D-01,    1.8260341504492 D-01,
     1    1.6915651939500 D-01,    1.4959598881657 D-01,
     1    1.2462897125553 D-01,    9.5158511682492 D-02,
     1    6.2253523938647 D-02,    2.7152459411754 D-02,
     1    1.2793819534675 D-01,    1.2583745634682 D-01,
     1    1.2167047292780 D-01,    1.1550566805372 D-01,
     1    1.0744427011596 D-01,    9.7618652104114 D-02,
     1    8.6190161531953 D-02,    7.3346481411080 D-02,
     1    5.9298584915436 D-02,    4.4277438817419 D-02,
     1    2.8531388628933 D-02,    1.2341229799987 D-02,
     1    9.6540088514727 D-02,    9.5638720079274 D-02,
     1    9.3844399080804 D-02,    9.1173878695763 D-02,
     1    8.7652093004403 D-02,    8.3311924226946 D-02/
      DATA (W(I),I=39,76)/
     1    7.8193895787070 D-02,    7.2345794108847 D-02,
     1    6.5822222776361 D-02,    5.8684093478535 D-02,
     1    5.0998059262376 D-02,    4.2835898022226 D-02,
     1    3.4273862913021 D-02,    2.5392065309262 D-02,
     1    1.6274394730905 D-02,    7.0186100094699 D-03,
     1    6.4737696812683 D-02,    6.4466164435950 D-02,
     1    6.3924238584648 D-02,    6.3114192286254 D-02,
     1    6.2039423159892 D-02,    6.0704439165893 D-02,
     1    5.9114839698394 D-02,    5.7277292100403 D-02,
     1    5.5199503699983 D-02,    5.2890189485193 D-02,
     1    5.0359035553853 D-02,    4.7616658492490 D-02,
     1    4.4674560856694 D-02,    4.1545082943464 D-02,
     1    3.8241351065829 D-02,    3.4777222564770 D-02,
     1    3.1167227832798 D-02,    2.7426509708356 D-02,
     1    2.3570760839324 D-02,    1.9616160457355 D-02,
     1    1.5579315722943 D-02,    1.1477234579234 D-02,
     1    7.3275539012760 D-03,    3.1533460523049 D-03,
     1    4.8690957009138 D-02,    4.8575467441503 D-02,
     1    4.8344762234802 D-02,    4.7999388596458 D-02/
      DATA (W(I),I=77,114)/
     1    4.7540165714830 D-02,    4.6968182816210 D-02,
     1    4.6284796581314 D-02,    4.5491627927418 D-02,
     1    4.4590558163756 D-02,    4.3583724529322 D-02,
     1    4.2473515123653 D-02,    4.1262563242622 D-02,
     1    3.9953741132720 D-02,    3.8550153178615 D-02,
     1    3.7055128540240 D-02,    3.5472213256882 D-02,
     1    3.3805161837141 D-02,    3.2057928354851 D-02,
     1    3.0234657072402 D-02,    2.8339672614259 D-02,
     1    2.6377469715053 D-02,    2.4352702568710 D-02,
     1    2.2270173808383 D-02,    2.0134823153530 D-02,
     1    1.7951715775697 D-02,    1.5726030476024 D-02,
     1    1.3463047896718 D-02,    1.1168139460131 D-02,
     1    8.8467598263630 D-03,    6.5044579689780 D-03,
     1    4.1470332605620 D-03,    1.7832807216960 D-03,
     1    3.2550614492363 D-02,    3.2516118713868 D-02,
     1    3.2447163714064 D-02,    3.2343822568575 D-02,
     1    3.2206204794029 D-02,    3.2034456231992 D-02,
     1    3.1828758894411 D-02,    3.1589330770726 D-02,
     1    3.1316425596861 D-02,    3.1010332586313 D-02/
      DATA (W(I),I=115,152)/
     1    3.0671376123669 D-02,    3.0299915420827 D-02,
     1    2.9896344136328 D-02,    2.9461089958167 D-02,
     1    2.8994614150555 D-02,    2.8497411065085 D-02,
     1    2.7970007616848 D-02,    2.7412962726029 D-02,
     1    2.6826866725591 D-02,    2.6212340735672 D-02,
     1    2.5570036005349 D-02,    2.4900633222483 D-02,
     1    2.4204841792364 D-02,    2.3483399085926 D-02,
     1    2.2737069658329 D-02,    2.1966644438744 D-02,
     1    2.1172939892191 D-02,    2.0356797154333 D-02,
     1    1.9519081140145 D-02,    1.8660679627411 D-02,
     1    1.7782502316044 D-02,    1.6885479864245 D-02,
     1    1.5970562902562 D-02,    1.5038721026994 D-02,
     1    1.4090941772314 D-02,    1.3128229566961 D-02,
     1    1.2151604671088 D-02,    1.1162102099838 D-02,
     1    1.0160770535008 D-02,    9.1486712307830 D-03,
     1    8.1268769256980 D-03,    7.0964707911530 D-03,
     1    6.0585455042350 D-03,    5.0142027429269 D-03,
     1    3.9645543384440 D-03,    2.9107318179340 D-03,
     1    1.8539607889460 D-03,    7.9679206555199 D-04/
      END
*DJKM 
      FUNCTION DJKM (JJ,KK,MM,X)
      IMPLICIT real*8 (A-H,O-Z) 
      J=JJ  
      K=KK  
      M=MM  
      IK=ABS(K)
      IM=ABS(M)
      IF(MOD(J,2).NE.MOD(IK,2).OR.MOD(J,2).NE.MOD(IM,2))GOTO7 
      IS=1  
      IF(IK.GE.IM)GOTO1 
      L=K 
      K=M 
      M=L 
      IM=IK 
      IF(2*MOD(ABS(K-M)/2,2).NE.0)IS=-1  
    1 IF(K.GE.0)GOTO2 
      K=-K  
      M=-M  
      IF(2*MOD(ABS(K-M)/2,2).NE.0)IS=-IS 
    2 IF(K.GT.J)GOTO8 
      KM=(K+IM)/2 
      D=1
      L1=(J+IM)/2+1 
      L2=(J+K)/2
      IF(L1.GT.L2)GOTO3 
      DO4L=L1,L2
    4 D=L*D/(L-KM)
    3 S=SIN(x/2) 
      C=COS(x/2) 
      CS=C**2-S**2
      D1=SQRT(D)
      IF(M.NE.-K)D1=D1*C**((K+M)/2) 
      IF(M.NE.K)D1=D1*S**((K-M)/2)
      J=J-1 
      D2=0
      IF(K+1.GT.J)GOTO6 
      D2=D1 
      D1=((K+2)*CS-M)*D1/2
      K3=K+3
      IF(K3.GT.J)GOTO6
      KM=K*M
      DO5LL=K3,J,2
      L=LL-1
      L2=L+2
      D=((L*L2*CS-KM)*D1*(L+L2)-L2*(L-M)*(L+M)*D2)  
     1 /(L*(L2+K)*(L2-K)) 
      D2=D1 
    5 D1=D  
    6 DJKM=D1*IS
      RETURN
    7 PRINT1000,J,KK,MM 
      RETURN
    8 PRINT1001,KK,MM,J 
      RETURN
 1000 FORMAT(/4H J =,I3,5H ,K =,I3,5H ,M =,I3,13H NE SONT PAS 
     1 ,37HSIMULTANEMENT ENTIERS OU DEMI-ENTIERS/)  
 1001 FORMAT(/4H K =,I3,7H OU M =,I3,19H PLUS GRAND QUE J =,I3/)
C**  THIS PROGRAM VALID ON FTN4 AND FTN5 ** 
      END 
      SUBROUTINE dede (A,IDIM1,N1,N2,INDEX,NERROR,DETERM)
      IMPLICIT REAL*8(A-H,O-Z)
C
C         MATRIX INVERSION WITH ACCOMPANYING SOLUTION OF LINEAR EQUATION
C
      DIMENSION A(IDIM1*IDIM1),INDEX(IDIM1)
      DETER=1
      N=N1
      IEMAT=N+N2
      IDIM=IDIM1
      NMIN1=N-1
C        THE ROUTINE DOES ITS OWN EVALUATION FOR DOUBLE SUBSCRIPTING OF
C        ARRAY A.
      IPIVC=1-IDIM
C        MAIN LOOP TO INVERT THE MATRIX
      DO 11 MAIN=1,N
      PIVOT=0
      IPIVC=IPIVC+IDIM
C        SEARCH FOR NEXT PIVOT IN COLUMN MAIN.
      IPIVC1=IPIVC+MAIN-1
      IPIVC2=IPIVC +NMIN1
      DO 2 I1=IPIVC1,IPIVC2
      IF(ABS(A(I1))-ABS(PIVOT)) 2,2,1
    1 PIVOT=A(I1)
      LPIV=I1
    2 CONTINUE
C        IS PIVOT DIFFERENT FROM ZERO
      IF(PIVOT) 3,15,3
C        GET THE PIVOT-LINE INDICATOR AND SWAP LINES IF NECESSARY
    3 ICOL=LPIV-IPIVC+1
      INDEX(MAIN)=ICOL
      IF(ICOL-MAIN) 6,6,4
C        COMPLEMENT THE DETERMINANT
    4 DETER=-DETER
C        POINTER TO LINE PIVOT FOUND
      ICOL=ICOL-IDIM
C        POINTER TO EXACT PIVOT LINE
      I3=MAIN-IDIM
      DO 5 I=1,IEMAT
      ICOL=ICOL+IDIM
      I3=I3+IDIM
      SWAP=A(I3)
      A(I3)=A(ICOL)
    5 A(ICOL)=SWAP
C        COMPUTE DETERMINANT
    6 DETER=DETER*PIVOT
      pivot=1/pivot
C        TRANSFORM PIVOT COLUMN
      I3=IPIVC+NMIN1
      DO 7 I=IPIVC,I3
    7 A(I)=-A(I)*PIVOT
      A(IPIVC1)=PIVOT
C        PIVOT ELEMENT TRANSFORMED
C
C        NOW CONVERT REST OF THE MATRIX
      I1=MAIN-IDIM
C        POINTER TO PIVOT LINE ELEMENTS
      ICOL=1-IDIM
C        GENERAL COLUMN POINTER
      DO 10 I=1,IEMAT
      ICOL=ICOL+IDIM
      I1=I1+IDIM
C        POINTERS MOVED
      IF(I-MAIN) 8,10,8
C        PIVOT COLUMN EXCLUDED
    8 JCOL=ICOL+NMIN1
      SWAP=A(I1)
      I3=IPIVC-1
      DO 9 I2=ICOL,JCOL
      I3=I3+1
    9 A(I2)=A(I2)+SWAP*A(I3)
      A(I1)=SWAP*PIVOT
   10 CONTINUE
   11 CONTINUE
      DETERM=DETER
      NERROR=0
      RETURN
   15 NERROR=MAIN
c     DETERM=DETER
      determ=0
      RETURN
      END
*ORDER
      SUBROUTINE ORDER(E,V,N,NN)
      IMPLICIT real*8 (A-H,O-Z) 
      DIMENSION E(NN),V(NN,NN)  
      IF(N.EQ.1) RETURN 
      IM=N-1
      DO 1 I=1,IM 
      JM=I+1
      DO 2 J=JM,N 
      IF(E(I).LT.E(J)) GOTO 2 
      Z=E(J)
      E(J)=E(I) 
      E(I)=Z
      DO 3 K=1,N
      ZZ   =V(K,J)
      V(K,J)=V(K,I) 
      V(K,I)=ZZ 
    3 CONTINUE
    2 CONTINUE
    1 CONTINUE
      RETURN
      END 
c     function second(a)
c     real*8 second
c     real*4time0(2),etime
c     second=etime(time0)
c     return
c     end
*CGAMMA 
      FUNCTION CGAMMA(Z)
      IMPLICIT real*8 (A-H,O-Z) 
      COMPLEX*16 Z,ZZ,Y,A5,CGAMMA,X  
      DIMENSION C(6),CZ(2),CY(2),CA5(2),CX(2) 
      EQUIVALENCE (CZ,ZZ),(CY,Y),(CA5,A5),(CX,X),(CX,ZR),(CX(2),ZIM)
      SAVE
      DATA PI/3.14159265358979/ 
      DATA C/76.180091729406,-86.505320327112,24.014098222230,
     1 -1.231739516140,0.001208580030,-0.000005363820/
      ZZ = Z
      CZ(1) = CZ(1)-1.
      REZ=CZ(1) 
      IF(REZ) 3,2,2 
    3 ZZ = -ZZ
    2 ZR = CZ(1)
      ZIM = -CZ(2)
      Z2 = ZIM**2 
      A5 = (1.000000000178,0.)  
      DO 1 J=1,6
      ZR = ZR+1.
      CY(2) = C(J)/(ZR**2+Z2) 
      CY(1) = ZR*CY(2)
      CY(2) = ZIM*CY(2) 
    1 A5 = A5+Y 
      CY(1) = CZ(1)+5.5 
      CY(2) = CZ(2) 
      CX(2)=CY(2) 
      CX(1) = CZ(1)+0.5 
      CA5(1) = CA5(1)*2.506628274631001 
      CA5(2) = CA5(2)*2.506628274631001 
      CGAMMA = EXP(X*LOG(Y)-Y)*A5 
      IF(REZ) 4,5,5 
    4 CZ(1) = PI*CZ(1)
      CZ(2) = PI*CZ(2)
      CGAMMA = ZZ/(SIN(ZZ)*CGAMMA) 
    5 RETURN
      END 

      subroutine bessj(n,x,y,z)
      implicit real*8(a-h,o-z)
      PARAMETER (IACC=40,BIGNO=1.d10,BIGNI=1.d-10)
      data pi/3.14159265358979d0/
      dimension y(n+1),z(n+1)
      allocatable yy(:)
      y=0
      z=0
      if(x.eq.0)then
      y(1)=1
      return
      endif
      call bessj0y0(x,y(1),z(1))
      if(n.eq.0)return
      y(2)=bessj1(x)
      z(2)=(y(2)*z(1)-2/pi/x)/y(1)
      if(n.eq.1)return
      ax=abs(x)
c     y(2)=abs(y(2))
      nlim=ax
      nlim=min(nlim,n)
      if(nlim.eq.0)nlim=1
      tox=2/ax
        do 11 j=1,nlim-1
          y(j+2)=j*tox*y(j+1)-y(j)
11      continue
        do 10 j=1,n-1
          z(j+2)=j*tox*z(j+1)-z(j)
10      continue
      if(nlim.eq.n)goto13
      xnor=y(nlim+1)
c       m=2*((n+int(sqrt(float(IACC*n))))/2)
        m=n+10
      allocate (yy(m+2))
        yy(m+2)=0
        yy(m+1)=tiny(x)
        do 12 j=m,nlim+1,-1
          yy(j)=j*yy(j+1)*tox-yy(j+2)
12      continue
        y(nlim+1:n+1)=yy(nlim+1:n+1)*xnor/yy(nlim+1)
  13  if(x.lt.0)then
      y(2:n+1:2)=-y(2:n+1:2)
      z=0
      endif
      return
      END

*BESSJ0
      subroutine bessj0y0(X,bessj0,bessy0)
      implicit real*8(a-h,o-z)
      DATA XBIG,XVSMAL,TBPI/1.0D+16,3.2D-9,6.36619772367581343D-1/
      T = ABS(X)
      IF (T.GT.8.0D0) GO TO 40
C     SMALL X
   20 T = 3.125D-2*T*T - 1
      T2 = 2*T
C     EXPANSION (0027) EVALUATED AS Y(T)  --PRECISION 17E.18
      A = +1.22200000000000000D-17
      B = T2*A - 7.58850000000000000D-16
      C = T2*B - A + 4.12532100000000000D-14
      A = T2*C - B - 1.94383469000000000D-12
      B = T2*A - C + 7.84869631400000000D-11
      C = T2*B - A - 2.67925353056000000D-9
      A = T2*C - B + 7.60816359241900000D-8
      B = T2*A - C - 1.76194690776215000D-6
      C = T2*B - A + 3.24603288210050800D-5
      A = T2*C - B - 4.60626166206275050D-4
      B = T2*A - C + 4.81918006946760450D-3
      C = T2*B - A - 3.48937694114088852D-2
      A = T2*C - B + 1.58067102332097261D-1
      B = T2*A - C - 3.70094993872649779D-1
      C = T2*B - A + 2.65178613203336810D-1
      A = T2*C - B - 8.72344235285222129D-3
      Y = T*A - C + 1.57727971474890120D-1
      bessj0 = Y
            A = +3.90000000000000000D-19
            B = T2*A - 2.69800000000000000D-17
            C = T2*B - A + 1.64349000000000000D-15
            A = T2*C - B - 8.74734100000000000D-14
            B = T2*A - C + 4.02633082000000000D-12
            C = T2*B - A - 1.58375525420000000D-10
            A = T2*C - B + 5.24879478733000000D-9
            B = T2*A - C - 1.44072332740190000D-7
            C = T2*B - A + 3.20653253765480000D-6
            A = T2*C - B - 5.63207914105698700D-5
            B = T2*A - C + 7.53113593257774230D-4
            C = T2*B - A - 7.28796247955207918D-3
            A = T2*C - B + 4.71966895957633869D-2
            B = T2*A - C - 1.77302012781143582D-1
            C = T2*B - A + 2.61567346255046637D-1
            A = T2*C - B + 1.79034314077182663D-1
            B = T2*A - C - 2.74474305529745265D-1
            G = T*B - A - 3.31461132032849417D-2
            bessy0 = TBPI*Y*LOG(X) + G
      GO TO 80
C
C     LARGE X
   40 G = T - 0.5D0/TBPI
      Y = SQRT(TBPI/T)
      CX = COS(G)*Y
      SX = -SIN(G)*Y*8/T
      CX2 = COS(G)*Y*8.0D0/X
      SX2 = SIN(G)*Y
      T = 128/(T*T) - 1
C     EXPANSION (0029) EVALUATED AS Y(T)  --PRECISION 17E.18
      Y = (((((((((((((-1.17964800000000000D-14)
     *    *T+3.34028800000000000D-14)*T-6.41843200000000000D-14)
     *    *T+2.45294080000000000D-13)*T-1.06365696000000000D-12)
     *    *T+4.78702080000000000D-12)*T-2.48827276800000000D-11)
     *    *T+1.55005642880000000D-10)*T-1.21211819632000000D-9)
     *    *T+1.28037614434400000D-8)*T-2.05274481565160000D-7)
     *    *T+6.13741608010926000D-6)*T-5.36367319213004570D-4)*T +
     *    9.99457275788251954D-1
C     EXPANSION (0030) EVALUATED AS G(T)  --PRECISION 17E.18
      G = ((((((((((((((((-9.83040000000000000D-16)
     *    *T+2.12992000000000000D-15)*T-1.14688000000000000D-15)
     *    *T+4.75136000000000000D-15)*T-2.27942400000000000D-14)
     *    *T+6.95193600000000000D-14)*T-2.28807680000000000D-13)
     *    *T+8.59855360000000000D-13)*T-3.59094272000000000D-12)
     *    *T+1.70330918400000000D-11)*T-9.47034774400000000D-11)
     *    *T+6.43278173280000000D-10)*T-5.66871613024000000D-9)
     *    *T+7.10621485930000000D-8)*T-1.47713883264594000D-6)
     *    *T+6.83314909934390000D-5)*T - 1.55551138795135187D-2
      bessj0 = Y*CX + G*SX
      bessy0 = G*CX2 + Y*SX2
   80 RETURN
      END

*BESSJ1
      FUNCTION bessj1(X)
      implicit Real*8(a-h,o-z)         
      DATA XBIG,XVSMAL,TBPI/1.0D+16,3.2D-9,6.36619772367581343D-1/
      T = ABS(X)
      IF (T.GT.8.0D0) GO TO 40
C     SMALL X
      Y = 4
C      TEST FOR VERY SMALL X
      IF (T.LE.XVSMAL) GO TO 20
      T = 3.125D-2*T*T - 1
      T2 = 2*T
      A = +2.95000000000000000D-18
      B = T2*A - 1.95540000000000000D-16
      C = T2*B - A + 1.13857200000000000D-14
      A = T2*C - B - 5.77740420000000000D-13
      B = T2*A - C + 2.52812366400000000D-11
      C = T2*B - A - 9.42421298160000000D-10
      A = T2*C - B + 2.94970700727800000D-8
      B = T2*A - C - 7.61758780540030000D-7
      C = T2*B - A + 1.58870192399321300D-5
      A = T2*C - B - 2.60444389348580680D-4
      B = T2*A - C + 3.24027018268385747D-3
      C = T2*B - A - 2.91755248061542077D-2
      A = T2*C - B + 1.77709117239728283D-1
      B = T2*A - C - 6.61443934134543253D-1
      C = T2*B - A + 1.28799409885767762D+0
      A = T2*C - B - 1.19180116054121687D+0
      Y = T*A - C + 6.48358770605264921D-1
   20 bessj1 = Y*X/8
      GO TO 80
C
C     LARGE X
   40 G = T - 1.5D0/TBPI
      Y = SIGN(SQRT(TBPI/T),X)
      CX = COS(G)*Y
      SX = -SIN(G)*Y*8.0D0/T
      T = 128.0D0/(T*T) - 1.0D0
      Y = (((((((((((((+1.24928000000000000D-14)
     *    *T-3.54508800000000000D-14)*T+6.86387200000000000D-14)
     *    *T-2.63372800000000000D-13)*T+1.14622464000000000D-12)
     *    *T-5.19113984000000000D-12)*T+2.72040729600000000D-11)
     *    *T-1.71310758400000000D-10)*T+1.36030580128000000D-9)
     *    *T-1.47085129889600000D-8)*T+2.45367662227560000D-7)
     *    *T-7.95969469843846000D-6)*T+8.98804941670557880D-4)*T +
     *    1.00090702627808217D+0
      G = (((((((((((((((-2.29376000000000000D-15)
     *    *T+5.32480000000000000D-15)*T-4.83328000000000000D-15)
     *    *T+1.75718400000000000D-14)*T-7.43936000000000000D-14)
     *    *T+2.50224640000000000D-13)*T-9.23228160000000000D-13)
     *    *T+3.87554432000000000D-12)*T-1.85248000000000000D-11)
     *    *T+1.03950931840000000D-10)*T-7.15110504800000000D-10)
     *    *T+6.42013250344000000D-9)*T-8.29196070929200000D-8)
     *    *T+1.82120185123076000D-6)*T-9.62145882205441600D-5)*T +
     *    4.67768740274489776D-2
      bessj1 = Y*CX + G*SX
   80 RETURN
      END
      SUBROUTINE cMATIN (A,IDIM1,N1,N2,INDEX,NERROR,DETERM)
      IMPLICIT complex*16(A-H,O-Z)
      real*8 pivot2
C
C         MATRIX INVERSION WITH ACCOMPANYING SOLUTION OF LINEAR EQUATION
C
      DIMENSION A(IDIM1*IDIM1),INDEX(IDIM1)
      DETER=1
      N=N1
      IEMAT=N+N2
      IDIM=IDIM1
      NMIN1=N-1
C        THE ROUTINE DOES ITS OWN EVALUATION FOR DOUBLE SUBSCRIPTING OF
C        ARRAY A.
      IPIVC=1-IDIM
C        MAIN LOOP TO INVERT THE MATRIX
      DO 11 MAIN=1,N
      PIVOT=0
      IPIVC=IPIVC+IDIM
C        SEARCH FOR NEXT PIVOT IN COLUMN MAIN.
      IPIVC1=IPIVC+MAIN-1
      IPIVC2=IPIVC +NMIN1
      DO 2 I1=IPIVC1,IPIVC2
      IF(ABS(A(I1))-ABS(PIVOT)) 2,2,1
    1 PIVOT=A(I1)
      LPIV=I1
    2 CONTINUE
C        IS PIVOT DIFFERENT FROM ZERO
      pivot2=pivot
      IF(PIVOT2) 3,15,3
C        GET THE PIVOT-LINE INDICATOR AND SWAP LINES IF NECESSARY
    3 ICOL=LPIV-IPIVC+1
      INDEX(MAIN)=ICOL
      IF(ICOL-MAIN) 6,6,4
C        COMPLEMENT THE DETERMINANT
    4 DETER=-DETER
C        POINTER TO LINE PIVOT FOUND
      ICOL=ICOL-IDIM
C        POINTER TO EXACT PIVOT LINE
      I3=MAIN-IDIM
      DO 5 I=1,IEMAT
      ICOL=ICOL+IDIM
      I3=I3+IDIM
      SWAP=A(I3)
      A(I3)=A(ICOL)
    5 A(ICOL)=SWAP
C        COMPUTE DETERMINANT
    6 DETER=DETER*PIVOT
      PIVOT=1/PIVOT
C        TRANSFORM PIVOT COLUMN
      I3=IPIVC+NMIN1
      DO 7 I=IPIVC,I3
    7 A(I)=-A(I)*PIVOT
      A(IPIVC1)=PIVOT
C        PIVOT ELEMENT TRANSFORMED
C
C        NOW CONVERT REST OF THE MATRIX
      I1=MAIN-IDIM
C        POINTER TO PIVOT LINE ELEMENTS
      ICOL=1-IDIM
C        GENERAL COLUMN POINTER
      DO 10 I=1,IEMAT
      ICOL=ICOL+IDIM
      I1=I1+IDIM
C        POINTERS MOVED
      IF(I-MAIN) 8,10,8
C        PIVOT COLUMN EXCLUDED
    8 JCOL=ICOL+NMIN1
      SWAP=A(I1)
      I3=IPIVC-1
      DO 9 I2=ICOL,JCOL
      I3=I3+1
    9 A(I2)=A(I2)+SWAP*A(I3)
      A(I1)=SWAP*PIVOT
   10 CONTINUE
   11 CONTINUE
C        NOW REARRANGE THE MATRIX TO GET RIGHT INVERS
      DO 14 I1=1,N
      MAIN=N+1-I1
      LPIV=INDEX(MAIN)
      IF(LPIV-MAIN) 12,14,12
   12 ICOL=(LPIV-1)*IDIM+1
      JCOL=ICOL+NMIN1
      IPIVC=(MAIN-1)*IDIM+1-ICOL
      DO 13 I2=ICOL,JCOL
      I3=I2+IPIVC
      SWAP=A(I2)
      A(I2)=A(I3)
   13 A(I3)=SWAP
   14 CONTINUE
      DETERM=DETER
      NERROR=0
      RETURN
   15 NERROR=MAIN
      DETERM=DETER
      RETURN
      END

      SUBROUTINE REGSOL(N,H,NV,POT,E,CI,SOL)
      implicit real*8 (a-h,o-z)
c     parameter(nv=2)
      DIMENSION X1(NV),X2(NV),AU(NV),CI(NV),POT(NV,NV,N),SOL(NV,N)
      H2=H*H
      AP=12/H2
* INITIALISATION DE LA SOLUTION
      DO2 J=1,NV
      SOL(J,1)=CI(J)
    2 CONTINUE
      DO3 J=1,NV
        A=0
        DO4 K=1,NV
          B=POT(J,K,1)
          IF(K.EQ.J)B=B-AP-E
    4     A=A+SOL(K,1)*B
        X1(J)=0
        X2(J)=-A/AP
    3 CONTINUE
      DO8 J=1,NV
        A1=0
        DO5 K=1,NV
          B=POT(J,K,1)
          IF(J.EQ.K)B=B-E
    5     A1=A1+B*X2(K)
        A2=0
        DO6 K=1,NV
          B1=POT(J,K,1)
          IF(J.EQ.K)B1=B1-E
          B2=0
          DO7 M=1,NV
            B=POT(K,M,1)
            IF(K.EQ.M)B=B-E
    7       B2=B2+B*X2(M)
    6     A2=A2+B1*B2
        AU(J)=H2*(A1+A2/AP)
    8   CONTINUE
* ITERATION
      DO10 I=2,N
        DO11 K=1,NV
          Z=2*X2(K)-X1(K)+AU(K)
          X1(K)=X2(K)
          X2(K)=Z
   11     CONTINUE
        DO13 J=1,NV
          A1=0
          DO12 K=1,NV
            B=POT(J,K,I)
            IF(J.EQ.K)B=B-E
   12       A1=A1+B*X2(K)
          A2=0
          DO14 K=1,NV
            B1=POT(J,K,I)
            IF(J.EQ.K)B1=B1-E
            B2=0
            DO15 M=1,NV
              B=POT(K,M,I)
              IF(K.EQ.M)B=B-E
   15         B2=B2+B*X2(M)
   14       A2=A2+B1*B2
          AU(J)=H2*(A1+A2/AP)
          SOL(J,I)=X2(J)+AU(J)/12
   13     CONTINUE
   10   CONTINUE
      RETURN
      END

***********************************************************************
      SUBROUTINE REGSOLI(N,H,NV,POT,E,CI,SOL)
      implicit real*8 (a-h,o-z)
c     parameter(nv=2)
      DIMENSION X1(NV),X2(NV),AU(NV),CI(NV),POT(NV,NV,N),SOL(NV,N)
      H2=H*H
      AP=12/H2
* INITIALISATION DE LA SOLUTION
      DO  J=1,NV
        SOL(J,N)=CI(J)
        SOL(J,N-1)=CI(J)*exp(sqrt(-e)*h)
      ENDDO
      DO3 J=1,NV
        A=0
        DO4 K=1,NV
          B=POT(J,K,N)
          IF(K.EQ.J)B=B-AP-E
    4     A=A+SOL(K,N)*B
        X1(J)=-A/AP
    3 CONTINUE
      DO33 J=1,NV
        A=0
        DO44 K=1,NV
          B=POT(J,K,N-1)
          IF(K.EQ.J)B=B-AP-E
   44     A=A+SOL(K,N-1)*B
        X2(J)=-A/AP
   33 CONTINUE
      DO8 J=1,NV
        A1=0
        DO5 K=1,NV
          B=POT(J,K,N-1)
          IF(J.EQ.K)B=B-E
    5     A1=A1+B*X2(K)
        A2=0
        DO6 K=1,NV
          B1=POT(J,K,N-1)
          IF(J.EQ.K)B1=B1-E
          B2=0
          DO7 M=1,NV
            B=POT(K,M,N-1)
            IF(K.EQ.M)B=B-E
    7       B2=B2+B*X2(M)
    6     A2=A2+B1*B2
        AU(J)=H2*(A1+A2/AP)
    8   CONTINUE
* ITERATION
      DO10 I=N-2,1,-1
        DO11 K=1,NV
          Z=2*X2(K)-X1(K)+AU(K)
          X1(K)=X2(K)
          X2(K)=Z
   11     CONTINUE
        DO13 J=1,NV
          A1=0
          DO12 K=1,NV
            B=POT(J,K,I)
            IF(J.EQ.K)B=B-E
   12       A1=A1+B*X2(K)
          A2=0
          DO14 K=1,NV
            B1=POT(J,K,I)
            IF(J.EQ.K)B1=B1-E
            B2=0
            DO15 M=1,NV
              B=POT(K,M,I)
              IF(K.EQ.M)B=B-E
   15         B2=B2+B*X2(M)
   14       A2=A2+B1*B2
          AU(J)=H2*(A1+A2/AP)
          SOL(J,I)=X2(J)+AU(J)/12
   13     CONTINUE
   10   CONTINUE
      RETURN
      END
      SUBROUTINE gauleg(x1,x2,x,w,n)
      INTEGER n
      REAL*8 x1,x2,x(n),w(n)
      DOUBLE PRECISION EPS
      PARAMETER (EPS=3.d-14)
! EPS is the relative precision.
! Given the lower and upper limits of integration x1 and x2, and given n, this routine returns
! arrays x(1:n) and w(1:n) of length n, containing the abscissas and weights of the Gauss-Legendre n-point quadrature formula.
      INTEGER i,j,m
      DOUBLE PRECISION p1,p2,p3,pp,xl,xm,z,z1
!High precision is a good idea for this routine.
      m=(n+1)/2 !The roots are symmetric in the interval, so we only have to find half of them.
      xm=0.5d0*(x2+x1)
      xl=0.5d0*(x2-x1)

      do i=1,m
!Loop over the desired roots.
      z=cos(3.141592654d0*(i-.25d0)/(n+.5d0))

! Starting with the above approximation to the ith root, we enter the main loop of refinement by Newton’s method.
    1 continue
      p1=1.d0
      p2=0.d0
      do j=1,n
! Loop up the recurrence relation to get the Legendre polynomial evaluated at z.
      p3=p2
      p2=p1
      p1=((2.d0*j-1.d0)*z*p2-(j-1.d0)*p3)/j
      enddo
! p1 is now the desired Legendre polynomial. We next compute pp, its derivative, by a standard relation involving also p2, the polynomial of one lower order.
      pp=n*(z*p1-p2)/(z*z-1.d0)
      z1=z
      z=z1-p1/pp
! Newton’s method.
      if(abs(z-z1).gt.EPS)goto 1
      x(i)=xm-xl*z    !Scale the root to the desired interval,
      x(n+1-i)=xm+xl*z   !and put in its symmetric counterpart.
!      write(*,*) x(i)
      w(i)=2.d0*xl/((1.d0-z*z)*pp*pp)  !Compute the weight
      w(n+1-i)=w(i)                     !and its symmetric counterpart.
      enddo
      return
      END     
