                                                C
C     Program for coupled-channels equations for heavy-ion 
C     fusion reactions
C
C     The couplings are always treated to all orders.
C
C     Anharmonic couplings and generalised rotational couplings (but for E2 
C     only: E4 couplings are yet to be worked out) are included in this 
C     version of CCFULL. 
C
C     last modification: December 14th, 2010
c                        June 29th, 2010: closed channels
c                        December 14th, 2010: delete the option "P=1-R"
C
      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER*1 ANS
C
C  Parameters:
C    NLEVELMAX - maximum number of intrinsic levels to be included
C
      PARAMETER (NLEVELMAX=30,NRMAX=5000)
C
C Principal global variables
C
C    P         - penetrability 
C    SIGMA     - fusion cross section, unit mb
C    SPIN      - mean angular momentum
C
C    AP,ZP,RP  - atomic #, proton # and radius of the projectile
C    AT,ZT,RT  - those of the target
C    RMASS     - reduced mass, unit MeV/c**2
C    ANMASS    - nucleon mass, unit MeV/c**2
C    E         - bombarding energy in the center of mass frame, unit MeV
C    V0,R0,A0  - depth, range, and dissuseness parameters of uncoupled 
C                nuclear potential, which is assumed to be a Woods-Saxon form
C    IVIBROT   - option for intrinsic degree of freedom 
C                = -1; no excitation (inert)/ =0 ; vibrational coupling/ 
C                =1; rotational coupling
C    BETA      - defeormation parameter
C    OMEGA     - excitation energy of the oscillator
C    LAMBDA    - multipolarity
C    NPHONON   - the number of phonon to be included
C    BETA2     - static quadrupole deformation parameter
C    BETA4     - static hexadecapole deformation parameter
C    E2        - excitation energy of 2+ state in a rotational band
C    NROT      - the number of levels in the rotational band to be included 
C                (up to I^pi=2*NROT+ states are included)
C    L         - angular momentum of the relative motion
C    CPOT      - coupling matrix
C    

        dimension cpot0(nlevelmax,nlevelmax)

	common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi
	common/intr/ivibrotp,ivibrott,nproj,ntarg
        common/oscp/betap,betapn,omegap,lambdap,nphononp
	common/rotp/beta2p,beta4p,e2p,nrotp
        common/osct/betat,betatn,omegat,lambdat,nphonont
        common/osct2/betat2,betat2n,omegat2,lambdat2,nphonont2
	common/rott/beta2t,beta4t,e2t,nrott
	common/dimen/nlevel
        common/pot/v0,r0,a0,pow
	common/dyn/e,rmass
        common/angular/l
	common/shape/rb,vb,curv
	common/grid/rmin,rmax,dr
        common/coup/cpot(nlevelmax,nlevelmax,0:nrmax)
        common/couph/cpoth(nlevelmax,nlevelmax)                  
        common/trans/ftr,qtrans,iq,ntrans
	common/mut/imutual(0:nlevelmax,0:nlevelmax)
        common/eigen/ev(nlevelmax),evec(nlevelmax,nlevelmax)

C
C Output files: 
C   OUTPUT: detailed information of the calculation
C   cross.dat : fusion cross sections
C   spin.dat  : averaged angular momenta
C   partial.dat: partial cross sections (being calculated only when a single 
C                          value of the energy is entered in the input file)
C
	OPEN(7,FILE='cross.dat',STATUS='UNKNOWN')
	OPEN(8,FILE='spin.dat',STATUS='UNKNOWN')
	OPEN(9,FILE='OUTPUT',STATUS='UNKNOWN')
	OPEN(10,FILE='ccfull.inp',STATUS='UNKNOWN')

C
C Define two constants used in various places
C
      hbar=197.329d0
      pi=3.141592653d0
C
C Input phase
C
      anmass=938.d0

      ap=16.
      zp=8.
      at=148.
      zt=62.
	read(10,*)ap,zp,at,zt
C  r_coup, unit fm
      rp=1.2d0*ap**(1.d0/3.d0)
      rt=1.06d0*at**(1.d0/3.d0)
	read(10,*)rp,ivibrotp,rt,ivibrott
	rp=rp*ap**(1.d0/3.d0)
	rt=rt*at**(1.d0/3.d0)

      rmass=ap*at/(ap+at)*anmass

C-----
c      IVIBROTT=0
C             = -1 for inert, =0 for vibration, =1 for rotation
	if(ivibrott.eq.-1) then
	      ntarg=0
              read(10,*)omega,beta,alambda,nphonon
	elseif(ivibrott.eq.0) then
c (input for target phonon)
              omegat=0.55d0
              betat=0.182d0
              lambdat=2
              nphonont=3
              read(10,*)omegat,betat,lambdat,nphonont
	      ntarg=nphonont
	else
c (input for target rotation)
	      e2t=0.d0
	      beta2t=0.d0
	      beta4t=0.d0
              nrott=0
              lambdat=2
              read(10,*)e2t,beta2t,beta4t,nrott
	      ntarg=nrott
	endif
C-----
C (input for target phonon (the second mode of excitation))
              nphonont2=3
c                     =0; no second mode in the target
              omegat2=1.16d0
              betat2=0.236d0
              lambdat2=3
              read(10,*)omegat2,betat2,lambdat2,nphonont2
	      if(ivibrott.eq.-1) nphonont2=0
c-----
c      ivibrotp=-1
c             = -1 for inert, =0 for vibration, =1 for rotation
	if(ivibrotp.eq.-1) then
	      nproj=0
              read(10,*)omegap,betap,alambdap,nphononp
	else if(ivibrotp.eq.0) then
c (input for projectile phonon)
              omegap=3.74d0
              betap=0.339d0
              lambdap=3
              nphononp=1
              read(10,*)omegap,betap,lambdap,nphononp
              nproj=nphononp
	else
c (input for projectile rotation)
              e2p=0.d0
	      beta2p=0.d0
	      beta4p=0.d0
	      nrotp=0
              read(10,*)e2p,beta2p,beta4p,nrotp
	      nproj=nrotp
	endif

      ntrans=0
c ntrans= 0 ; no transfer channel/ =1 ; with transfer channel
c ftrans(r)=ftr*dvn(r)/dr
c
      ftr=7.d0
      qtrans=5.21d0
        	read(10,*)ntrans,qtrans,ftr

c--------------- parameters for the nuclear potential
              r1=1.233d0*ap**(1.d0/3.d0)-0.98d0*ap**(-1.d0/3.d0)
              r2=1.233d0*at**(1.d0/3.d0)-0.98d0*at**(-1.d0/3.d0)
              r0=r1+r2+0.29d0
              r12=r1*r2/(r1+r2)
              a0=0.63d0
              gamma=0.95d0*(1.d0-1.8d0*(ap-2.d0*zp)/ap*(at-2.d0*zt)/at)
              v0=16.d0*pi*gamma*r12*a0
              v0=31.67*r1*r2/(r1+r2)

                v0=1551d0
                r0=0.95d0*(ap**(1.d0/3.d0)+at**(1.d0/3.d0))
                a0=1.05d0
                        read(10,*)v0,r0,a0
	                r0=r0*(ap**(1.d0/3.d0)+at**(1.d0/3.d0))

        pow=1.d0
	WRITE(6,*)' Modified Woods-Saxon (power of WS) (n/y)?'
	read(5,1)ans
 1	format(a1)
	     IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
	     WRITE(6,*)'power=?'
	     read(5,*)pow
                if(pow.eq.0.d0) then
                WRITE(6,*)'Power cannot be zero :^('
                stop
                endif
	     endif

c  the minimum energy, the maximum energy, and the energy step

	emin=53.d0
	emax=72.d0
	de=2.d0*at/(ap+at)/6.d0
        	read(10,*)emin,emax,de

      rmax=50.d0
      dr=0.05d0
        	read(10,*)rmax,dr

c========================================== end of input phase

	if(ivibrott.ne.-1.and.ntarg.eq.0.and.nphonont2.ne.0) then
	write(6,91)
 91	FORMAT(1H ,'The target excitation should be input'
     &  /' in the first mode :^(')
	stop
	endif

c	if(ivibrott.eq.1.and.ivibrotp.eq.1) then
c	write(6,92)
c 92	format(1h ,'this code cannot handle fusion between'
c     &  / 'two deformed nuclei :^(')
c	stop
c	endif

c	if(ivibrott.eq.1.and.nphonont2.ne.0) then
c 	write(6,*)'!!!!!!!!!!!!!!! warning !!!!!!!!!!!!!!!!!!'
c	write(6,*)'    '
c	write(6,*)'the vibrational phonon in the target is not'
c        write(6,*)'a deformed phonon, but a *spherical* phonon,'
c        write(6,*)'i.e. the rotation-vibration model is not used.'
c	write(6,*)'    '
c 	write(6,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
c	write(6,*)'    '
c	endif

         call nucleus
	 nlevel=nlevel+ntrans
	 if(nlevel.gt.nlevelmax) then
	 write(6,*)'too many channels :^('
	 stop
	 endif
 100	write(6,*) 'nlevel=',nlevel

C Evaluation coupling matrix elements
C------------------------------------------------------
      l=0
C   search a minimum of the Coulomb pocket RMIN
      call potshape
      iterat=int((rmax-rmin)/dr)+1.d-6
             if(iterat+1.gt.nrmax) then
	     write(6,*)'too large grid !!  :^('
	     stop
	     endif
      rmax=rmin+dr*iterat

      call cmat0
      do 29 ir=0,iterat+1
      r=rmin+dr*ir
      call cmat(r,cpot0)
      do 31 i=1,nlevel
      do 32 j=1,nlevel
         cpot(i,j,ir)=cpot0(i,j)
 32      continue
 31      continue
 29      continue

      rh=rmin+dr/2.d0                
      call cmat(rh,cpoth)            

	write(9,99)
	write(9,97)
 99	FORMAT(/,'       Ecm (MeV)    sigma (mb)          <l> ')
 97	format('       -------------------------------------')

        iwarn=0

c  start the energy loop
c-------------------------------------------------------------
	if(de.eq.0.d0) then
	iestep=0
	else
        iestep=(emax-emin)/de+1.d-6
	endif

	if(iestep.eq.0) then
	open(11,file='partial.dat',status='unknown')
	endif

      do 10 i=0,iestep
      e=emin+de*i

      sigma=0.d0
      spin=0.d0
c  angular momentum loop
c-----------------------
      do 20 l=0,200
      s0=sigma

      rmin0=rmin
      call potshape
      rmin00=rmin
      rmin=rmin0
      if(v(rmin00).gt.e.or.rmin00.lt.0.d0) then
      p=0.d0
        if(l.ne.0.and.iwarn.eq.0) then
        write(6,*)'    '
        WRITE(6,*)'!!!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!!!!!!!'
        WRITE(6,*)'Shallow potential!!!!'
        WRITE(6,*)'You may get an oscillatioin in BD at E > ', E
        WRITE(6,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
        WRITE(6,*)'    '
        WRITE(9,*)'    '
        WRITE(9,*)'!!!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!!!!!!'
        WRITE(9,*)'Shallow potential!!!!'
        WRITE(9,*)'You may get an oscillatioin in BD at E > ', E
        WRITE(9,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
        write(9,*)'    '
        iwarn=1
        endif
      go to 15
      endif

c
c integration of the c.c. equations
      call mnumerov(p)

c	if(l.eq.0) write(9,*)e,p
	if(nlevel.gt.10) write(6,*)e,l*1.,p

      if(iestep.eq.0) then
      write(11,*)e,l*1.,(2.d0*l+1.d0)*p*pi*hbar**2/2.d0/rmass/e*10.d0
      endif

      sigma=sigma+(2.d0*l+1.d0)*p*pi*hbar**2/2.d0/rmass/e*10.d0
      spin=spin+l*(2.d0*l+1.d0)*p*pi*hbar**2/2.d0/rmass/e*10.d0
      if(sigma-s0.lt.s0*1.d-4) goto 15
 20   continue

 15   if(sigma.ne.0.d0) spin=spin/sigma

      write(7,*)e,sigma
      write(8,*)e,spin
       if(sigma.lt.1.d-2) then
       write(6,81)e,sigma,spin
       write(9,81)e,sigma,spin
       else
       write(6,82)e,sigma,spin
       write(9,82)e,sigma,spin
       endif

 10   continue
 81   format(f15.5,e15.5,f15.5)
 82   format(3f15.5)

      WRITE(6,*)'The program terminated normally :^)'	
      if(iwarn.eq.1) then
      write(6,*)'    '
      WRITE(6,*)'.....but you may get an oscillation in BD'
      WRITE(6,*)' (see OUTPUT file).'
      endif

      stop
      end

c*************************************************************
      subroutine nucleus
c
c subroutine to record several input parameters in the 
c 'output' file
c
c*************************************************************
      implicit real*8(a-h,o-z)
      parameter (nlevelmax=30)
      character*1 ans,sgn
      character*2 text,pro,targ
      dimension text(109)
	common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi
	common/intr/ivibrotp,ivibrott,nproj,ntarg
        common/oscp/betap,betapn,omegap,lambdap,nphononp
	common/rotp/beta2p,beta4p,e2p,nrotp
        common/osct/betat,betatn,omegat,lambdat,nphonont
        common/osct2/betat2,betat2n,omegat2,lambdat2,nphonont2
	common/rott/beta2t,beta4t,e2t,nrott
	common/phonon/nphonon
	common/dimen/nlevel
        common/pot/v0,r0,a0,pow
	common/dyn/e,rmass
	common/shape/rb,vb,curv
	common/grid/rmin,rmax,dr
        common/trans/ftr,qtrans,iq,ntrans
	common/mut/imutual(0:nlevelmax,0:nlevelmax)
        common/ahv/betnahv(0:nlevelmax,0:nlevelmax),
     &             betcahv(0:nlevelmax,0:nlevelmax),
     &             omeahv(0:nlevelmax)
        common/ahv2/betnahv2(0:nlevelmax,0:nlevelmax),
     &             betcahv2(0:nlevelmax,0:nlevelmax),
     &             omeahv2(0:nlevelmax)
        common/ahvp/betnahvp(0:nlevelmax,0:nlevelmax),
     &             betcahvp(0:nlevelmax,0:nlevelmax),
     &             omeahvp(0:nlevelmax)
        common/grott/bett(0:nlevelmax,0:nlevelmax),erott(0:nlevelmax)
        common/grotp/betp(0:nlevelmax,0:nlevelmax),erotp(0:nlevelmax)

      DATA TEXT/'H ','He','Li','Be','B ','C ','N ','O ','F ','Ne',
     &'Na','Mg','Al','Si','P ','S ','Cl','Ar','K ','Ca','Sc','Ti','V ',
     &'Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr',
     &'Rb','Sr','Y ','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In',
     &'Sn','Sb','Te','I ','Xe','Cs','Ba','La','Ce','Pr','Nd','Pm','Sm',
     &'Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb','Lu','Hf','Ta','W ','Re',
     &'Os','Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn','Fr','Ra',
     &'Ac','Th','Pa','U ','Np','Pu','Am','Cm','Bk','Cf','Es','Fm','Md',
     &'No','Lr','XX','X1','X2','X3','X4','04'/
 
        nlevel=(nproj+1)*(ntarg+1)
	np=zp
	npp=ap
	pro=text(np)
	nt=zt
	ntt=at
	targ=text(nt)
	WRITE(9,*)NPP,PRO,'  +  ',NTT,TARG, '     Fusion reaction'
	WRITE(6,*)NPP,PRO,'  +  ',NTT,TARG, '     Fusion reaction'

	write(9,30)
	write(6,30)

	r0p=rp/ap**(1.d0/3.d0)
	r0t=rt/at**(1.d0/3.d0)

	if(ntarg.ne.0) then
	if(ivibrott.eq.0) then
	write(6,21)betat,omegat,lambdat,nphonont 
 21	FORMAT('Phonon Excitation in the targ.: beta=',F6.3,
     &    ', omega=',F5.2, '(MeV), Lambda=',I2,', Nph=',I2)

	betatn=betat
	write(6,*)'    '
	WRITE(6,*)' Different beta_N from beta_C for this mode(n/y)?'
	READ(5,1)ANS
	     IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
	     WRITE(6,*)'beta_N=?'
	     read(5,*)betatn
	     endif
	write(9,211)betatn,betat,r0t,omegat,lambdat,nphonont 
 211	FORMAT('Phonon Excitation in the targ.: beta_N=',F6.3,
     &    ', beta_C=',F6.3,', r0=',F5.2,'(fm),',
     &     /32X,'omega=',F5.2, '(MeV), Lambda=',I2,', Nph=',I2)
	endif

	if(ivibrott.eq.1) then
	write(9,22)beta2t,beta4t,r0t,e2t,nrott
	write(6,22)beta2t,beta4t,r0t,e2t,nrott
 22     FORMAT('Rotational Excitation in the targ.: beta2=',F6.3,
     &    ', beta4=',F6.3,', r0=',F5.2,'(fm),',
     &      /36X,'E2=',F5.2,'(MeV), Nrot=',I2)
	endif
	endif

	if(nphonont2.ne.0) then
	write(6,23)betat2,omegat2,lambdat2,nphonont2 
 23	FORMAT('Phonon Excitation in the targ.: beta=',F6.3,
     &    ', omega=',F5.2, '(MeV), Lambda=',I2,', Nph=',I2)

	betat2n=betat2
	write(6,*)'    '
	WRITE(6,*)' Different beta_N from beta_C for this mode(n/y)?'
	READ(5,1)ANS
	     IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
	     WRITE(6,*)'beta_N=?'
	     read(5,*)betat2n
	     endif

	write(9,233)betat2n,betat2,r0t,omegat2,lambdat2,nphonont2 
 233	FORMAT('Phonon Excitation in the targ.: beta_N=',F6.3,
     &    ', beta_C=',F6.3,', r0=',F5.2,'(fm),',
     &    /32X, 'omega=',F5.2, '(MeV), Lambda=',I2,', Nph=',I2)

	call mutual
	endif

	if(nproj.ne.0) then
	if(ivibrotp.eq.0) then
	write(6,24)betap,omegap,lambdap,nphononp 
 24	FORMAT('Phonon Excitation in the proj.: beta=',F6.3,
     &    ', omega=',F5.2, '(MeV), Lambda=',I2,', Nph=',I2)

	betapn=betap
	write(6,*)'    '
	WRITE(6,*)' Different beta_N from beta_C for this mode(n/y)?'
	READ(5,1)ANS
	     IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
	     WRITE(6,*)'beta_N=?'
	     read(5,*)betapn
	     endif
	write(9,244)betapn,betap,r0p,omegap,lambdap,nphononp 
 244	FORMAT('Phonon Excitation in the proj.: beta_N=',F6.3,
     &    ', beta_C=',F6.3,', r0=',F5.2,'(fm),',
     &    /32X, 'omega=',F5.2, '(MeV), Lambda=',I2,', Nph=',I2)
	endif

	if(ivibrotp.eq.1) then
	write(9,25)beta2p,beta4p,r0p,e2p,nrotp
	write(6,25)beta2p,beta4p,r0p,e2p,nrotp
 25     FORMAT('Rotational Excitation in the proj.: beta2=',F6.3,
     &    ', beta4=',F6.3,', r0=',F5.2,'(fm),',
     &      /36X,'E2=',F5.2,'(MeV), Nrot=',I2)
	endif
	endif

        call anharmonicity
        call grotation

	if(ntrans.ne.0) then
	write(9,30)
	write(9,26)ftr, qtrans
	write(6,26)ftr, qtrans
 26	FORMAT('Transfer channel: Strength=',F5.2,
     &    ', Q=',F5.2, '(MeV)')
        IQ=0
        WRITE(6,*)'     '
        WRITE(6,*)'Which type (charge change)?'
        WRITE(6,*)'         [0 for n-transfer, 1 for p-pick up,', 
     &              ' -1 for p-stripping, etc.]'
        read(5,*)iq
        if(iq.gt.0) then
        sgn='+'
        write(6,291)sgn,iq
        write(9,291)sgn,iq
        else
        write(6,29)iq
        write(9,29)iq
        endif
 29     FORMAT('   Charge difference of the transfer partition:',I5)
 291    FORMAT('   Charge difference of the transfer partition:  ',
     &          a,i2)
        write(6,*)'     '
	endif

	write(9,30)

	r00=r0/(ap**(1.d0/3.d0)+at**(1.d0/3.d0))
	write(9,27)v0,r00,a0,pow
 27	FORMAT('Potential parameters: V0=', F8.2, '(MeV), r0=', 
     &        F5.2, '(fm), a=', F5.2, '(fm), power=',F5.2)
		CALL POTSHAPE
		WRITE(9,28)RB,VB,CURV
 28 	FORMAT('   Uncoupled barrier: Rb=', F5.2, '(fm), Vb=',F8.2, 
     &        '(MeV), Curv=',F5.2, '(MeV)')
	WRITE(9,30)
 30	FORMAT('-------------------------------------------------')

 1	FORMAT(A1)
      RETURN
      END

c*************************************************************
      subroutine anharmonicity
c
c anharmonic couplings for vibrations
c
c*************************************************************
      implicit real*8(a-h,o-z)
      parameter (nlevelmax=30)
      character*1 ans,sign
	common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi
	common/intr/ivibrotp,ivibrott,nproj,ntarg
        common/oscp/betap,betapn,omegap,lambdap,nphononp
	common/rotp/beta2p,beta4p,e2p,nrotp
        common/osct/betat,betatn,omegat,lambdat,nphonont
        common/osct2/betat2,betat2n,omegat2,lambdat2,nphonont2
	common/rott/beta2t,beta4t,e2t,nrott
	common/phonon/nphonon
	common/dimen/nlevel
        common/pot/v0,r0,a0,pow
	common/dyn/e,rmass
	common/shape/rb,vb,curv
	common/grid/rmin,rmax,dr
        common/trans/ftr,qtrans,iq,ntrans
	common/mut/imutual(0:nlevelmax,0:nlevelmax)
        common/ahv/betnahv(0:nlevelmax,0:nlevelmax),
     &             betcahv(0:nlevelmax,0:nlevelmax),
     &             omeahv(0:nlevelmax)
        common/ahv2/betnahv2(0:nlevelmax,0:nlevelmax),
     &             betcahv2(0:nlevelmax,0:nlevelmax),
     &             omeahv2(0:nlevelmax)
        common/ahvp/betnahvp(0:nlevelmax,0:nlevelmax),
     &             betcahvp(0:nlevelmax,0:nlevelmax),
     &             omeahvp(0:nlevelmax)

c-----------------------------
c ahv coupling for the first mode in the target
           do it=0,ntarg
           do jt=0,ntarg

           betnahv(it,jt)=0.d0
           betcahv(it,jt)=0.d0

           if(it.eq.jt-1) then
              betnahv(it,jt)=betatn*sqrt(jt*1.d0)
              betcahv(it,jt)=betat*sqrt(jt*1.d0)
           else if(jt.eq.it-1) then
              betnahv(it,jt)=betatn*sqrt(it*1.d0)
              betcahv(it,jt)=betat*sqrt(it*1.d0)
           endif

           enddo
              omeahv(it)=it*omegat
           enddo

        if(ivibrott.eq.0.and.ntarg.gt.1) then

	WRITE(6,90)'AHV couplings for the first mode ',
     &             'in the target phonon (n/y)?'
 90	format(2a)
	read(5,1)ans
        IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
           WRITE(9,*)'            '
           WRITE(9,90)'**** AHV Couplings in the target',
     &                ' (the 1st mode) ****'
           write(9,*)'             '
	if((-1.d0)**lambdat.eq.1.d0) sign='+'
	if((-1.d0)**lambdat.eq.-1.d0) sign='-'

           do it=0,ntarg
           do 44 jt=0,ntarg
           if(it.gt.jt) goto 44
           
           if(it.eq.0.and.jt.eq.0) goto 44
           if(it.eq.0.and.jt.eq.1) goto 44

        write(6,*)'             '
	WRITE(6,94)'Transition from the',LAMBDAT,SIGN,'^',IT,
     &                  ' to the',LAMBDAT,SIGN,'^',JT,'state:'
        WRITE(6,95)'   beta_N and beta_C in the HO limit=',
     &             BETNAHV(IT,JT),BETCAHV(IT,JT)
        WRITE(6,90)'    Modify these beta_N a/o beta_C (n/y)?'
	READ(5,1)ANS

        IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
	WRITE(9,94)'Transition from the',LAMBDAT,SIGN,'^',IT,
     &                  ' to the',LAMBDAT,SIGN,'^',JT,'state:'
        WRITE(9,95)'   beta_N and beta_C in the HO limit=',
     &             BETNAHV(IT,JT),BETCAHV(IT,JT)
        WRITE(6,*)'   beta_N and beta_C =?'
        READ(5,*)BETNAHV(IT,JT),BETCAHV(IT,JT)
        WRITE(9,97)BETNAHV(IT,JT),BETCAHV(IT,JT)
        ENDIF

 94     format(1h ,a,i2,2a,i2,a,i2,2a,i2,a)
 95     format(1h ,a,2g15.5)
 97     format('    beta_n=',g15.5,',   beta_c=',g15.5)

 44        enddo
           enddo

              WRITE(6,*)'           '
              WRITE(6,*)'Excitation energy for the first mode:'
              DO I=2,NTARG
              WRITE(6,98)'    Energy of the',I,
     &              '-phonon state in the HO=',I*OMEGAT
              WRITE(6,*)'   Modify this energy(n/y)?'
              READ(5,1)ANS
              IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
              WRITE(6,*)'   Energy=?'
              READ(5,*)OMEAHV(I)
              WRITE(6,*)'       '
              WRITE(9,61)'Energy of the',I,'-phonon state=',OMEAHV(I),
     &                         '(HO: ',I*OMEGAT,')'
              endif
              enddo
              write(9,*)'          '              
              write(6,*)'          '              
 98           format(a,i2,a,g15.5)
 61           format(1h ,a,i2,a,g15.5,a,g12.5,a)
            endif
            endif
c-----------------------------
c ahv coupling for the second mode in the target
           do it=0,nphonont2
           do jt=0,nphonont2

           betnahv2(it,jt)=0.d0
           betcahv2(it,jt)=0.d0

           if(it.eq.jt-1) then
              betnahv2(it,jt)=betat2n*sqrt(jt*1.d0)
              betcahv2(it,jt)=betat2*sqrt(jt*1.d0)
           else if(jt.eq.it-1) then
              betnahv2(it,jt)=betat2n*sqrt(it*1.d0)
              betcahv2(it,jt)=betat2*sqrt(it*1.d0)
           endif

           enddo
              omeahv2(it)=it*omegat2
           enddo

        if(nphonont2.gt.1) then

	WRITE(6,90)'AHV couplings for the second mode ',
     &             'in the target phonon (n/y)?'
	READ(5,1)ANS
        IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
           WRITE(9,*)'            '
           WRITE(9,90)'**** AHV Couplings in the target',
     &                ' (the 2nd mode) ****'
           write(9,*)'             '

	if((-1.d0)**lambdat2.eq.1.d0) sign='+'
	if((-1.d0)**lambdat2.eq.-1.d0) sign='-'

           do it=0,nphonont2
           do 54 jt=0,nphonont2
           if(it.gt.jt) goto 54
           
           if(it.eq.0.and.jt.eq.0) goto 54
           if(it.eq.0.and.jt.eq.1) goto 54

        write(6,*)'             '
	WRITE(6,94)'Transition from the',LAMBDAT2,SIGN,'^',IT,
     &                  ' to the',LAMBDAT2,SIGN,'^',JT,'state:'
        WRITE(6,95)'    beta_N and beta_C in the HO limit=',
     &             BETNAHV2(IT,JT),BETCAHV2(IT,JT)
        WRITE(6,90)'     Modify these beta_N a/o beta_C (n/y)?'
	READ(5,1)ANS

        IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
	WRITE(9,94)'Transition from the',LAMBDAT2,SIGN,'^',IT,
     &                  ' to the',LAMBDAT2,SIGN,'^',JT,'state:'
        WRITE(9,95)'   beta_N and beta_C in the HO limit=',
     &             BETNAHV(IT,JT),BETCAHV(IT,JT)
        WRITE(6,*)'    beta_N and beta_C =?'
        READ(5,*)BETNAHV2(IT,JT),BETCAHV2(IT,JT)
        WRITE(9,97)BETNAHV2(IT,JT),BETCAHV2(IT,JT)
        ENDIF

 54        ENDDO
           ENDDO

              WRITE(6,*)'           '
              WRITE(6,*)'Excitation energy for the second mode:'
              DO I=2,NPHONONT2
              WRITE(6,98)'    Energy of the',I,
     &              '-phonon state in the HO=',I*OMEGAT2
            WRITE(6,*)'   Modify this energy(n/y)?'
              READ(5,1)ANS
              IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
              WRITE(6,*)'   Energy=?'
              READ(5,*)OMEAHV2(I)
              WRITE(9,61)'Energy of the',I,'-phonon state=',OMEAHV2(I),
     &                         '(HO: ',I*OMEGAT2,')'
              ENDIF
              ENDDO
              WRITE(9,*)'          '              
              WRITE(6,*)'          '              
            ENDIF
            ENDIF

c-----------------------------
c ahv coupling for the projectile phonon 
           do it=0,nproj
           do jt=0,nproj

           betnahvp(it,jt)=0.d0
           betcahvp(it,jt)=0.d0

           if(it.eq.jt-1) then
              betnahvp(it,jt)=betapn*sqrt(jt*1.d0)
              betcahvp(it,jt)=betap*sqrt(jt*1.d0)
           else if(jt.eq.it-1) then
              betnahvp(it,jt)=betapn*sqrt(it*1.d0)
              betcahvp(it,jt)=betap*sqrt(it*1.d0)
           endif

           enddo
              omeahvp(it)=it*omegap
           enddo

        if(ivibrotp.eq.0.and.nproj.gt.1) then

	write(6,72)
 72     FORMAT('AHV couplings for the projectile phonon (n/y)?')
	READ(5,1)ANS
        IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
           WRITE(9,*)'            '
           WRITE(9,73)
 73        FORMAT(1H ,'**** AHV Couplings in the projectile ****') 
           write(9,*)'            '
	if((-1.d0)**lambdap.eq.1.d0) sign='+'
	if((-1.d0)**lambdap.eq.-1.d0) sign='-'

           do it=0,nproj
           do 64 jt=0,nproj
           if(it.gt.jt) goto 64
           
           if(it.eq.0.and.jt.eq.0) goto 64
           if(it.eq.0.and.jt.eq.1) goto 64

        write(6,*)'             '
	WRITE(6,94)'Transition from the',LAMBDAP,SIGN,'^',IT,
     &                  ' to the',LAMBDAP,SIGN,'^',JT,'state:'
        WRITE(6,95)'  beta_N and beta_C in the HO limit=',
     &             BETNAHVP(IT,JT),BETCAHVP(IT,JT)
        WRITE(6,90)'   Modify these beta_N a/o beta_C (n/y)?'
	read(5,1)ans

        IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
	WRITE(9,94)'Transition from the',LAMBDAP,SIGN,'^',IT,
     &                  ' to the',LAMBDAP,SIGN,'^',JT,'state:'
        WRITE(9,95)'   beta_N and beta_C in the HO limit=',
     &             BETNAHVP(IT,JT),BETCAHVP(IT,JT)
        WRITE(6,*)'  beta_N and beta_C =?'
        READ(5,*)BETNAHVP(IT,JT),BETCAHVP(IT,JT)
        WRITE(9,97)BETNAHVP(IT,JT),BETCAHVP(IT,JT)
        ENDIF

 64        ENDDO
           ENDDO

              WRITE(6,*)'           '
           WRITE(6,*)'Excitation energy for the projectile phonon:'
           DO I=2,NPROJ
              WRITE(6,98)'    Energy of the',I,
     &              '-phonon state in the HO=',I*OMEGAP
            WRITE(6,*)'   Modify this energy(n/y)?'
              READ(5,1)ANS
              IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
              WRITE(6,*)'   Energy=?'
              READ(5,*)OMEAHVP(I)
              WRITE(9,61)'Energy of the',I,'-phonon state=',OMEAHVP(I),
     &                         '(HO: ',I*OMEGAP,')'
              ENDIF
              ENDDO
              WRITE(9,*)'          '              
              WRITE(6,*)'          '              
            ENDIF
            ENDIF

 1	FORMAT(A1)

      RETURN
      END

c*************************************************************
      subroutine grotation
c
c generalised e2 rotational couplings
c
c*************************************************************
      implicit real*8(a-h,o-z)
      parameter (nlevelmax=30)
      character*1 ans
	common/hion/ap,zp,rp,at,zt,rt
	common/intr/ivibrotp,ivibrott,nproj,ntarg
	common/rotp/beta2p,beta4p,e2p,nrotp
	common/rott/beta2t,beta4t,e2t,nrott
        common/grott/bett(0:nlevelmax,0:nlevelmax),erott(0:nlevelmax)
        common/grotp/betp(0:nlevelmax,0:nlevelmax),erotp(0:nlevelmax)

c-----------------------------
c the generalised couplings for the target rotor
           do it=0,ntarg
           do jt=0,ntarg
           bett(it,jt)=0.d0
           if(it.eq.jt-1.or.it.eq.jt+1.or.it.eq.jt) bett(it,jt)=beta2t
           enddo
           erott(it)=2.d0*it*(2.d0*it+1.d0)/6.d0*e2t
           enddo

        if(ivibrott.eq.1.and.ntarg.gt.1) then

        WRITE(6,*)'      '
	WRITE(6,*)' Generalised E2 couplings for the target rotor (n/y)?'
	READ(5,1)ANS
        IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
           WRITE(9,*)'            '
           WRITE(9,90)'**** Generalised couplings for the ',
     &                ' target rotor ****'
           WRITE(9,*)'             '
           WRITE(6,*)'             '

           WRITE(6,*)'Beta2_targ=',BETA2T

           DO IT=1,NTARG-1
           JT=IT+1

           WRITE(6,95)'Modify beta2 for the transition from',
     &                2*IT,'+ to',2*JT,'+ ? (n/y)'
   	   READ(5,1)ANS

           IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
           WRITE(6,*)'   BETA2=?'           
           READ(5,*)BETT(IT,JT)
               BETT(JT,IT)=BETT(IT,JT)
           WRITE(9,94)'Beta2 for the transition from',2*IT,'+ to',
     &                2*JT,'+ =: ',BETT(IT,JT)
           ENDIF
           ENDDO

              WRITE(6,*)'           '
              WRITE(6,*)'Reorientation terms:'

           DO IT=1,NTARG
           JT=IT

           WRITE(6,95)'Modify beta2 for the transition from',
     &                2*IT,'+ to',2*JT,'+ ? (n/y)'
   	   READ(5,1)ANS

           IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
           WRITE(6,*)'   BETA2=?'           
           READ(5,*)BETT(IT,JT)
               BETT(JT,IT)=BETT(IT,JT)
           WRITE(9,94)'Beta2 for the transition from',2*IT,'+ to',
     &                2*JT,'+ =: ',BETT(IT,JT)
           ENDIF
           ENDDO


              WRITE(6,*)'           '
              WRITE(6,*)'Excitation energies:'
              DO I=2,NTARG
              WRITE(6,98)'    Energy of the',2*I,
     &              '+ state for the pure rotor=',EROTT(I)
              WRITE(6,*)'   Modify this energy(n/y)?'
              READ(5,1)ANS
              IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
                 EROTT0=EROTT(I)
              WRITE(6,*)'   Energy=?'
              READ(5,*)EROTT(I)
              WRITE(6,*)'       '
              WRITE(9,61)'Energy of the',2*I,'+ state=',EROTT(I),
     &                         '(rotor: ',EROTT0,')'
              ENDIF
              ENDDO
              WRITE(9,*)'          '              
              WRITE(6,*)'          '              
 98           FORMAT(A,I2,A,G15.5)
 61           FORMAT(1H ,A,I2,A,G15.5,A,G12.5,A)
            ENDIF
            ENDIF

 94     FORMAT(1H ,A,I2,A,I2,A,G15.5)
 95     FORMAT(1H ,A,I2,A,I2,A)
 90	FORMAT(2A)

C-----------------------------
c The generalised couplings for the projectile rotor
           DO IP=0,NPROJ
           DO JP=0,NPROJ
           BETP(IP,JP)=0.D0
           IF(IP.EQ.JP-1.OR.IP.EQ.JP+1.OR.IP.EQ.JP) BETP(IP,JP)=BETA2P
           ENDDO
           EROTP(IP)=2.D0*IP*(2.D0*IP+1.D0)/6.D0*E2P
           ENDDO

        IF(IVIBROTP.EQ.1.AND.NPROJ.GT.1) THEN

        WRITE(6,*)'      '
	WRITE(6,90)' Generalised E2 couplings for the projectile', 
     &             ' rotor (n/y)?'
	READ(5,1)ANS
        IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
           WRITE(9,*)'            '
           WRITE(9,90)'**** Generalised couplings for the ',
     &                ' projectile rotor ****'
           WRITE(9,*)'             '
           WRITE(6,*)'             '

           WRITE(6,*)'Beta2_proj=',BETA2P

           DO IP=1,NPROJ-1
           JP=IP+1

           WRITE(6,95)'Modify beta2 for the transition from',
     &                2*IP,'+ to',2*JP,'+ ? (n/y)'
   	   READ(5,1)ANS

           IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
           WRITE(6,*)'   BETA2=?'           
           READ(5,*)BETP(IP,JP)
               BETP(JP,IP)=BETP(IP,JP)
           WRITE(9,94)'Beta2 for the transition from',2*IP,'+ to',
     &                2*JP,'+ =: ',BETP(IP,JP)
           ENDIF
           ENDDO

              WRITE(6,*)'           '
              WRITE(6,*)'Reorientation terms:'

           DO IP=1,NPROJ
           JP=IP

           WRITE(6,95)'Modify beta2 for the transition from',
     &                2*IP,'+ to',2*JP,'+ ? (n/y)'
   	   READ(5,1)ANS

           IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
           WRITE(6,*)'   BETA2=?'           
           READ(5,*)BETP(IP,JP)
               BETP(JP,IP)=BETP(IP,JP)
           WRITE(9,94)'Beta2 for the transition from',2*IP,'+ to',
     &                2*JP,'+ =: ',BETP(IP,JP)
           ENDIF
           ENDDO


              WRITE(6,*)'           '
              WRITE(6,*)'Excitation energies:'
              DO I=2,NPROJ
              WRITE(6,98)'    Energy of the',2*I,
     &              '+ state for the pure rotor=',EROTP(I)
              WRITE(6,*)'   Modify this energy(n/y)?'
              READ(5,1)ANS
              IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
                 EROTP0=EROTP(I)
              WRITE(6,*)'   Energy=?'
              READ(5,*)EROTP(I)
              WRITE(6,*)'       '
              WRITE(9,61)'Energy of the',2*I,'+ state=',EROTP(I),
     &                         '(rotor: ',EROTP0,')'
              ENDIF
              ENDDO
              WRITE(9,*)'          '              
              WRITE(6,*)'          '              
            ENDIF
            ENDIF

 1	FORMAT(A1)

        RETURN
        END

c*****************************************************************************
      subroutine mutual
c
c subroutine to select levels to be included in the c.c. calculations
c
c*****************************************************************************
        implicit real*8(a-h,m,o-z)
        parameter (nlevelmax=30)
	character*1 ans,sign1,sign2
	common/intr/ivibrotp,ivibrott,nproj,ntarg
        common/osct/betat,betatn,omegat,lambdat,nphonont
        common/osct2/betat2,betat2n,omegat2,lambdat2,nphonont2
	common/rott/beta2t,beta4t,e2t,nrott
	common/dimen/nlevel
        common/trans/ftr,qtrans,iq,ntrans
	common/mut/imutual(0:nlevelmax,0:nlevelmax)

	do 10 i=0,nlevelmax
	do 20 j=0,nlevelmax
	imutual(i,j)=0
 20	continue
 10	continue

	write(6,*)'    '
	WRITE(6,*)'Mutual excitations in the *target* nucleus'

	WRITE(6,*)'  Include the mutual excitations (y/n)?'
	READ(5,1)ANS
 1	FORMAT(A1)
	     IF(ANS.EQ.'n' .OR. ANS.EQ.'N') THEN
	     IMUT=0
	     NLEVEL=(NTARG+NPHONONT2+1)*(NPROJ+1)
c	     WRITE(9,99)
c 99	     FORMAT(1H ,'No mutual excitations in the target are included.')
	     DO 80 I=0,NTARG
	     DO 81 J=0,NPHONONT2
	     IF(I.EQ.0.OR.J.EQ.0) IMUTUAL(I,J)=1
 81          CONTINUE
 80	     CONTINUE
	     GOTO 70
	     ENDIF

	WRITE(6,*)'  All the possible mutual excitation channels (n/y)?'
	READ(5,1)ANS
	     IF(ANS.EQ.'Y' .OR. ANS.EQ.'y') THEN
	     IMUT=1
	     NLEVEL=(NTARG+1)*(NPHONONT2+1)*(NPROJ+1)
c	     WRITE(9,90)
c 90	     FORMAT(1H ,'All the possible mutual excitation channels'
c     &                  1H ,'in the target are included.')
	     DO 82 I=0,NTARG
	     DO 83 J=0,NPHONONT2
	     IMUTUAL(I,J)=1
 83          CONTINUE
 82	     CONTINUE
	     GOTO 70
	     ENDIF

	IMUT=2
	IF((-1.D0)**LAMBDAT.EQ.1.D0) SIGN1='+'
	IF((-1.D0)**LAMBDAT.EQ.-1.D0) SIGN1='-'
	IF((-1.D0)**LAMBDAT2.EQ.1.D0) SIGN2='+'
	IF((-1.D0)**LAMBDAT2.EQ.-1.D0) SIGN2='-'
	DO 30 I=0,NTARG
	DO 40 J=0,NPHONONT2
		IF(I.EQ.0.OR.J.EQ.0) THEN 
		IMUTUAL(I,J)=1
		GOTO 40
		ENDIF
	WRITE(6,94)'  Include (',LAMBDAT,SIGN1,'^',I,',',LAMBDAT2,SIGN2,
     &            '^',J,') state ? (y/n)'
	READ(5,1)ANS
	     IF(ANS.EQ.'n' .OR. ANS.EQ.'N') THEN
	     IMUTUAL(I,J)=0
	     GOTO 40
	     ENDIF
	IMUTUAL(I,J)=1
 40	CONTINUE
 30	CONTINUE

 70	WRITE(9,*)'    '
	WRITE(9,93)
	WRITE(6,93)
 93	FORMAT(1H ,'  Excited states in the target to be included: ')

	NLEVEL=0
	IF((-1.D0)**LAMBDAT.EQ.1.D0) SIGN1='+'
	IF((-1.D0)**LAMBDAT.EQ.-1.D0) SIGN1='-'
	IF((-1.D0)**LAMBDAT2.EQ.1.D0) SIGN2='+'
	IF((-1.D0)**LAMBDAT2.EQ.-1.D0) SIGN2='-'
	DO 50 I=0,NTARG
	DO 60 J=0,NPHONONT2
		IF(IMUTUAL(I,J).EQ.0) GOTO 60
	WRITE(9,94)'    (',LAMBDAT,SIGN1,'^',I,',',LAMBDAT2,SIGN2,
     &                   '^',J,') state'
	WRITE(6,94)'    (',LAMBDAT,SIGN1,'^',I,',',LAMBDAT2,SIGN2,
     &                   '^',J,') state'
 94	FORMAT(1H ,A,I2,2A,I2,A,I2,2A,I2,A)
	NLEVEL=NLEVEL+1
 60	CONTINUE
 50	CONTINUE
        NLEVEL=NLEVEL*(NPROJ+1)

	RETURN
	END
c*****************************************************************************
      subroutine potshape
c
c shape of the bare coulomb barrier, i.e. the barrier position, the curvature,
c and the position of the coulomb pocket
c
c*****************************************************************************
      implicit real*8(a-h,m,o-z)
      common/shape/rb,vb,curv
      common/grid/rmin,rmax,dr
      common/dyn/e,rmass
      common/const/hbar,pi
      external dv,v

        r=50.5d0
        u0=dv(r)
 10	r=r-1.d0

           if(r.lt.0.d0) then
              rb=-5.d0
	      rmin=-5.d0
              return
              endif

	u1=dv(r)
	if(u0*u1.gt.0.d0) goto 10
		ra=r+1.d0
		rb=r
	        tolk=1.d-6
                n=log10(abs(rb-ra)/tolk)/log10(2.d0)+0.5d0
                do 20 i=1,n
                r=(ra+rb)/2.d0
                u=dv(r)
                if(u0*u.lt.0.d0) then
                rb=r
                else
                ra=r
                endif
 20             continue

	rb=r

        ddv00=(dv(rb+1.d-5)-dv(rb-1.d-5))/2.d-5

	if(ddv00.gt.0.d0) then
	write(6,*)'something is strange :^('
	write(8,*)'something is strange :^('
	stop
	endif

	curv=hbar*sqrt(abs(ddv00)/rmass)
	vb=v(rb)

		ra=rb-0.5d0
		u0=dv(ra)
		rbb=0.5d0
		u1=dv(rbb)
	        tolk=1.d-6
                n=log10(abs(rb-ra)/tolk)/log10(2.d0)+0.5d0
                do 21 i=1,n
                r=(ra+rbb)/2.d0
                u=dv(r)
                if(u0*u.lt.0.d0) then
                rbb=r
                else
                ra=r
                endif
 21             continue

	rmin=r

c	write(6,*)'rmin=',rmin,v(rmin)
c       write(6,*)'rb=',rb,v(rb)
c       write(6,*)'curv=',curv00

	return
	end

c**************************************************************
      subroutine mnumerov00(p)
c
c subroutine for integration of the c.c. eqs. by modified 
c numerov method
c
c**************************************************************
      implicit real*8 (a-h,o-z)
      external v
      parameter (nlevelmax=30,nrmax=5000)

      dimension psi(nlevelmax),psi0(nlevelmax),psi1(nlevelmax)
      dimension xi(nlevelmax),xi0(nlevelmax),xi1(nlevelmax)
      dimension phi0(nlevelmax)
      dimension bb(nlevelmax,nlevelmax),
     &          bin(nlevelmax,nlevelmax)
      dimension bb2(nlevelmax,nlevelmax)
      dimension cc(nlevelmax,nlevelmax),
     &          cin(nlevelmax,nlevelmax)
      dimension dd0(nlevelmax,nlevelmax),dd1(nlevelmax,nlevelmax)
      dimension dd(nlevelmax)
      dimension fcw(0:200),gcw(0:200),fpcw(0:200),gpcw(0:200)
      dimension sigmad(0:200),iexp(0:200)
      dimension ech(nlevelmax),ech2(nlevelmax)
 
	common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi
	common/dimen/nlevel
	common/dyn/e,rmass
        common/angular/l
	common/grid/rmin,rmax,dr
        common/coup/cpot(nlevelmax,nlevelmax,0:nrmax)

      complex*16 k,kk,k2
      complex*16 psi,psi0,psi1,phi0
      complex*16 xi,xi0,xi1
      complex*16 ai
      complex*16 cwup0,cwdown0,cwup1,cwdown1
      complex*16 bb,bb2,bin,cc,cin,dd
      complex*16 dummy,dummy2

      ai=(0.d0,1.d0)
c define two constants used in this subroutine
      fac=dr**2*(2.d0*rmass/hbar**2) 
      iterat=int((rmax-rmin)/dr)
      rmax=rmin+iterat*dr

	do 22 lc=0,200
	fcw(lc)=0.d0
	gcw(lc)=0.d0
	fpcw(lc)=0.d0
	gpcw(lc)=0.d0
	sigmad(lc)=0.d0
	iexp(lc)=0d0
 22	continue

      do 51 i=1,nlevelmax
      do 52 j=1,nlevelmax
      bb(i,j)=0.d0
      bin(i,j)=0.d0
      cc(i,j)=0.d0
      cin(i,j)=0.d0
 52   continue
      dd(i)=0.d0
 51   continue

c integration of the io-th channel wave function from rmin 
c==========================================================
      do 15 io=1,nlevel
c	write(6,*)'io=',io*1.

      do 200 j1=1,nlevel
      psi(j1)=0.d0
      psi0(j1)=0.d0
      psi1(j1)=0.d0
      phi0(j1)=0.d0
      xi0(j1)=0.d0
      xi1(j1)=0.d0
      xi(j1)=0.d0
 200  continue

        if(io.eq.1) then
        do io2=1,nlevel
        ech(io2)=e-v(rmin)-cpot(io2,io2,0)
        ech2(io2)=e-cpot(io2,io2,iterat)
        enddo
        endif

c  initial conditions
c----------
	if(ech(io).gt.0.d0) then
        k=sqrt(2.d0*rmass/hbar**2*ech(io))
        psi0(io)=exp(-ai*k*rmin)
        phi0(io)=-ai*k*psi0(io)
	else
        k=sqrt(2.d0*rmass/hbar**2*abs(ech(io)))
        psi0(io)=exp(k*rmin)
        phi0(io)=k*psi0(io)
	endif

      call rkutta00(psi0,phi0,psi1)

      do 91 i0=1,nlevel
	xi0(i0)=(1.d0-fac/12.d0*(v(rmin)-e))*psi0(i0)
	xi1(i0)=(1.d0-fac/12.d0*(v(rmin+dr)-e))*psi1(i0)
      do 92 ic=1,nlevel
      xi0(i0)=xi0(i0)-fac/12.d0*cpot(i0,ic,0)*psi0(ic)
      xi1(i0)=xi1(i0)-fac/12.d0*cpot(i0,ic,1)*psi1(ic)
 92   continue
 91   continue

      do 29 ir=2,iterat+1
c----------------------------------------------  iterations start
      r=rmin+dr*ir
      r0=rmin+dr*(ir-2)
      r1=rmin+dr*(ir-1)

      do 71 i0=1,nlevel
      do 72 ic=1,nlevel
         dd0(i0,ic)=fac/sqrt(12.d0)*cpot(i0,ic,ir-1)
      if(i0.eq.ic) dd0(i0,ic)=dd0(i0,ic)+fac/sqrt(12.d0)*(v(r1)-e)
     &                        +sqrt(3.d0)
 72   continue
 71   continue

      do 73 i0=1,nlevel
      do 74 ic=1,nlevel
      dd1(i0,ic)=0.d0
      if(i0.eq.ic) dd1(i0,ic)=dd1(i0,ic)-1.d0
      do 75 ik=1,nlevel
         dd1(i0,ic)=dd1(i0,ic)+dd0(i0,ik)*dd0(ik,ic)
 75   continue
 74   continue
 73   continue
	
      do 54 i0=1,nlevel
         xi(i0)=-xi0(i0)
      do 53 ic=1,nlevel
         xi(i0)=xi(i0)+dd1(i0,ic)*xi1(ic)
 53   continue
 54   continue

      if(ir.eq.iterat+1) goto 66
      do 42 i0=1,nlevel
         xi0(i0)=xi1(i0)
         xi1(i0)=xi(i0)
 42   continue
 29   continue
c--------------------------------------------------------------
c  matching to the coulomb wave function at rmax

 66   do 56 i0=1,nlevel
      do 55 ic=1,nlevel
      cc(i0,ic)=-fac/12.d0*cpot(i0,ic,iterat-1)
      if(i0.eq.ic) cc(i0,ic)=cc(i0,ic)-fac/12.d0*(v(rmax-dr)-e)+1.d0
 55   continue
 56   continue
      call matinv(nlevel,cc,cin)
      do 94 i0=1,nlevel
         psi0(i0)=0.d0
      do 93 ic=1,nlevel
         psi0(i0)=psi0(i0)+cin(i0,ic)*xi0(ic)
 93   continue
 94   continue

      do 556 i0=1,nlevel
      do 555 ic=1,nlevel
      cc(i0,ic)=-fac/12.d0*cpot(i0,ic,iterat+1)
      if(i0.eq.ic) cc(i0,ic)=cc(i0,ic)-fac/12.d0*(v(rmax+dr)-e)+1.d0
 555   continue
 556   continue
      call matinv(nlevel,cc,cin)
      do 944 i0=1,nlevel
         psi(i0)=0.d0
      do 933 ic=1,nlevel
         psi(i0)=psi(i0)+cin(i0,ic)*xi(ic)
 933   continue
 944   continue

      do 60 ii=1,nlevel
c   coulomb wave function
      ec=e-cpot(ii,ii,iterat-1)
      rho=sqrt(2.d0*rmass*ec)/hbar*(rmax-dr)
      eta=(zp*zt/137.d0)*sqrt(rmass/2.d0/ec)
      call dfcoul(eta,rho,fcw,fpcw,gcw,gpcw,sigmad,l,iexp)
      cwup0=(gcw(l)+ai*fcw(l))
      cwdown0=gcw(l)-ai*fcw(l)

      ec=e-cpot(ii,ii,iterat+1)
      rho=sqrt(2.d0*rmass*ec)/hbar*(rmax+dr)
      eta=(zp*zt/137.d0)*sqrt(rmass/2.d0/ec)
      call dfcoul(eta,rho,fcw,fpcw,gcw,gpcw,sigmad,l,iexp)
      cwup1=(gcw(l)+ai*fcw(l))
      cwdown1=gcw(l)-ai*fcw(l)

      bb(io,ii)=(cwup0*psi(ii)-cwup1*psi0(ii))
     &                         /(cwup0*cwdown1-cwup1*cwdown0)
      bb2(io,ii)=(cwdown1*psi0(ii)-cwdown0*psi(ii))
     &                         /(cwup0*cwdown1-cwup1*cwdown0)
 60   continue
c---------------------------------------------------------------
 15   continue

c===============================================================
c                                        penetration probability
c  calculate the inversion of the matrix bb
      call matinv(nlevel,bb,bin)
c
      p=0.d0
      ref=0.d0

      do 85 io=1,nlevel
      if(ech(io).lt.0.d0) goto 85
      k=sqrt((2.d0*rmass/hbar**2*ech(io)))
      k2=sqrt((2.d0*rmass/hbar**2*ech2(io)))
      kk=sqrt(2.d0*rmass/hbar**2*e)
      
             dummy=0.d0
             dummy2=0.d0
             do io2=1,nlevel
                dummy2=dummy2+bin(1,io2)*bb2(io2,io)
                enddo

      p=p+(abs(bin(1,io)))**2*k/kk
      ref=ref+(abs(dummy2))**2*k2/kk
 85   continue

c      if(p+ref.gt.1.1d0) then
c         write(6,*)'penetrability > 1.1 !!!!!'
c         write(6,*)'something must be strange (accuracy problem) :^('
c         stop
c      endif

      if(p.gt.0.99d0) p=1.d0-ref

      return
      end

c**************************************************************
      subroutine mnumerov0(p)
c
c subroutine for integration of the c.c. eqs. by modified 
c numerov method
c
c a similar stabilization method as that given in 
c   z.h. levine, pra30('84) 1120
c   t.n. rescigno and a.e. orel, pra25('82)2402
c is used in the routine. 
c
c the linear independence condition is imposed at the barrier position.
c
c**************************************************************
      implicit real*8 (a-h,o-z)
      external v
      parameter (nlevelmax=30,istab0=5,nrmax=5000)

      dimension psi(nlevelmax),psi0(nlevelmax),psi1(nlevelmax)
      dimension xi(nlevelmax,nlevelmax),xi0(nlevelmax,nlevelmax),
     &          xi1(nlevelmax,nlevelmax)
      dimension phi0(nlevelmax)
      dimension bb(nlevelmax,nlevelmax),
     &          bin(nlevelmax,nlevelmax)
      dimension bb2(nlevelmax,nlevelmax)
      dimension cc(nlevelmax,nlevelmax),
     &          cin(nlevelmax,nlevelmax)
      dimension dd0(nlevelmax,nlevelmax),dd1(nlevelmax,nlevelmax)
      dimension dd(nlevelmax)
      dimension fcw(0:200),gcw(0:200),fpcw(0:200),gpcw(0:200)
      dimension sigmad(0:200),iexp(0:200)
      dimension ech(nlevelmax),ech2(nlevelmax)
 
	common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi
	common/dimen/nlevel
	common/dyn/e,rmass
        common/angular/l
	common/grid/rmin,rmax,dr
        common/coup/cpot(nlevelmax,nlevelmax,0:nrmax)
	common/shape/rb,vb,curv
        
      complex*16 k,kk,k2
      complex*16 psi,psi0,psi1,phi0
      complex*16 xi,xi0,xi1
      complex*16 ai
      complex*16 cwup0,cwdown0,cwup1,cwdown1
      complex*16 bb,bb2,bin,cc,cin,dd
      complex*16 dummy,dummy2
      complex*16 xi1d(nlevelmax,nlevelmax)
      complex*16 xi0d(nlevelmax,nlevelmax)
      complex*16 aa(nlevelmax,nlevelmax)
      complex*16 bb0(nlevelmax,nlevelmax)
      complex*16 bb20(nlevelmax,nlevelmax)

      ai=(0.d0,1.d0)
c define two constants used in this subroutine
      fac=dr**2*(2.d0*rmass/hbar**2) 
      iterat=int((rmax-rmin)/dr)
      rmax=rmin+iterat*dr

      ibarrier=(rb-rmin)/dr+1.d-6

	do 22 lc=0,200
	fcw(lc)=0.d0
	gcw(lc)=0.d0
	fpcw(lc)=0.d0
	gpcw(lc)=0.d0
	sigmad(lc)=0.d0
	iexp(lc)=0d0
 22	continue

      do 51 i=1,nlevelmax
      do 52 j=1,nlevelmax
      bb(i,j)=0.d0
      bin(i,j)=0.d0
      cc(i,j)=0.d0
      cin(i,j)=0.d0
      aa(i,j)=0.d0
 52   continue
      dd(i)=0.d0
 51   continue

      do i=1,nlevel
         aa(i,i)=1.d0
         enddo

c  initial conditions
c==========================================================
      do 15 io=1,nlevel
c	write(6,*)'io=',io*1.

      do 200 j1=1,nlevel
      psi(j1)=0.d0
      psi0(j1)=0.d0
      psi1(j1)=0.d0
      phi0(j1)=0.d0
 200  continue

        if(io.eq.1) then
        do io2=1,nlevel
        ech(io2)=e-v(rmin)-cpot(io2,io2,0)
        ech2(io2)=e-cpot(io2,io2,iterat)
        enddo
        endif
        
	if(ech(io).gt.0.d0) then
        k=sqrt(2.d0*rmass/hbar**2*ech(io))
        psi0(io)=exp(-ai*k*rmin)
        phi0(io)=-ai*k*psi0(io)
	else
        k=sqrt(2.d0*rmass/hbar**2*abs(ech(io)))
        psi0(io)=exp(k*rmin)
        phi0(io)=k*psi0(io)
        endif

      call rkutta00(psi0,phi0,psi1)

      do 91 i0=1,nlevel
	xi0(i0,io)=(1.d0-fac/12.d0*(v(rmin)-e))*psi0(i0)
	xi1(i0,io)=(1.d0-fac/12.d0*(v(rmin+dr)-e))*psi1(i0)
      do 92 ic=1,nlevel
      xi0(i0,io)=xi0(i0,io)-fac/12.d0*cpot(i0,ic,0)*psi0(ic)
      xi1(i0,io)=xi1(i0,io)-fac/12.d0*cpot(i0,ic,1)*psi1(ic)
 92   continue
 91   continue
 15   continue

c===============
      istab=0
      do 29 ir=2,iterat+1
c----------------------------------------------  iterations start
      r=rmin+dr*ir
      r0=rmin+dr*(ir-2)
      r1=rmin+dr*(ir-1)
      istab=istab+1

      do 71 i0=1,nlevel
      do 72 ic=1,nlevel
         dd0(i0,ic)=fac/sqrt(12.d0)*cpot(i0,ic,ir-1)
      if(i0.eq.ic) dd0(i0,ic)=dd0(i0,ic)+fac/sqrt(12.d0)*(v(r1)-e)
     &                        +sqrt(3.d0)
 72   continue
 71   continue

      do 73 i0=1,nlevel
      do 74 ic=1,nlevel
      dd1(i0,ic)=0.d0
      if(i0.eq.ic) dd1(i0,ic)=dd1(i0,ic)-1.d0
      do 75 ik=1,nlevel
         dd1(i0,ic)=dd1(i0,ic)+dd0(i0,ik)*dd0(ik,ic)
 75   continue
 74   continue
 73   continue
	
      do 16 ich=1,nlevel
      do 54 i0=1,nlevel
         xi(i0,ich)=-xi0(i0,ich)
      do 53 ic=1,nlevel
         xi(i0,ich)=xi(i0,ich)+dd1(i0,ic)*xi1(ic,ich)
 53   continue
 54   continue
 16   continue

      if(ir.eq.iterat+1) goto 66
      if(r.lt.12.d0.and.istab.eq.istab0) then
        call stabilize(xi1,xi,aa,ir)
        istab=0
        endif
c      if(ir.eq.ibarrier) call stabilize(xi1,xi,aa,ir)

      do ich=1,nlevel
      do i0=1,nlevel
         xi0(i0,ich)=xi1(i0,ich)
         xi1(i0,ich)=xi(i0,ich)
      enddo
      enddo

 29   continue
c--------------------------------------------------------------
c  matching to the coulomb wave function at rmax

 66   do 17 io=1,nlevel

      do 56 i0=1,nlevel
      do 55 ic=1,nlevel
      cc(i0,ic)=-fac/12.d0*cpot(i0,ic,iterat-1)
      if(i0.eq.ic) cc(i0,ic)=cc(i0,ic)-fac/12.d0*(v(rmax-dr)-e)+1.d0
 55   continue
 56   continue
      call matinv(nlevel,cc,cin)
      do 94 i0=1,nlevel
         psi0(i0)=0.d0
      do 93 ic=1,nlevel
         psi0(i0)=psi0(i0)+cin(i0,ic)*xi0(ic,io)
 93   continue
 94   continue

      do 556 i0=1,nlevel
      do 555 ic=1,nlevel
      cc(i0,ic)=-fac/12.d0*cpot(i0,ic,iterat+1)
      if(i0.eq.ic) cc(i0,ic)=cc(i0,ic)-fac/12.d0*(v(rmax+dr)-e)+1.d0
 555   continue
 556   continue
      call matinv(nlevel,cc,cin)
      do 944 i0=1,nlevel
         psi(i0)=0.d0
      do 933 ic=1,nlevel
         psi(i0)=psi(i0)+cin(i0,ic)*xi(ic,io)
 933   continue
 944   continue

      do 60 ii=1,nlevel
c   coulomb wave function
      ec=e-cpot(ii,ii,iterat-1)
      rho=sqrt(2.d0*rmass*ec)/hbar*(rmax-dr)
      eta=(zp*zt/137.d0)*sqrt(rmass/2.d0/ec)
      call dfcoul(eta,rho,fcw,fpcw,gcw,gpcw,sigmad,l,iexp)
      cwup0=(gcw(l)+ai*fcw(l))
      cwdown0=gcw(l)-ai*fcw(l)

      ec=e-cpot(ii,ii,iterat+1)
      rho=sqrt(2.d0*rmass*ec)/hbar*(rmax+dr)
      eta=(zp*zt/137.d0)*sqrt(rmass/2.d0/ec)
      call dfcoul(eta,rho,fcw,fpcw,gcw,gpcw,sigmad,l,iexp)
      cwup1=(gcw(l)+ai*fcw(l))
      cwdown1=gcw(l)-ai*fcw(l)

      bb0(ii,io)=(cwup0*psi(ii)-cwup1*psi0(ii))
     &                         /(cwup0*cwdown1-cwup1*cwdown0)
      bb20(ii,io)=(cwdown1*psi0(ii)-cwdown0*psi(ii))
     &                         /(cwup0*cwdown1-cwup1*cwdown0)
 60   continue
c---------------------------------------------------------------
 17   continue

c===============================================================
c                                        penetration probability
      do i=1,nlevel
      do j=1,nlevel
         bb(i,j)=0.d0
         bb2(i,j)=0.d0
         do j2=1,nlevel
            bb(i,j)=bb(i,j)+bb0(i,j2)*aa(j2,j)
            bb2(i,j)=bb2(i,j)+bb20(i,j2)*aa(j2,j)
            enddo
            enddo
            enddo

c  calculate the inversion of the matrix bb
      call matinv(nlevel,bb,bin)
c
      p=0.d0
      ref=0.d0

      do 85 io=1,nlevel
      if(ech(io).lt.0.d0) goto 85
      k=sqrt((2.d0*rmass/hbar**2*ech(io)))
      k2=sqrt((2.d0*rmass/hbar**2*ech2(io)))
      kk=sqrt(2.d0*rmass/hbar**2*e)
      
             dummy=0.d0
             dummy2=0.d0
             do io2=1,nlevel
                dummy2=dummy2+bb2(io,io2)*bin(io2,1)
                enddo

      p=p+(abs(bin(io,1)))**2*k/kk
      ref=ref+(abs(dummy2))**2*k2/kk
 85   continue

c      if(p+ref.gt.1.1d0) then
c         write(6,*)'penetrability > 1.1 !!!!!'
c         write(6,*)'something must be strange (accuracy problem) :^('
c         stop
c      endif

      if(p.gt.0.99d0) p=1.d0-ref

      return
      end

c**************************************************************
      subroutine mnumerov(p)
c
c subroutine for integration of the c.c. eqs. by modified 
c numerov method
c
c a similar stabilization method as that given in 
c   z.h. levine, pra30('84) 1120
c   t.n. rescigno and a.e. orel, pra25('82)2402
c is used in the routine. 
c
c the linear independence condition is imposed at the barrier position.
c
c**************************************************************
      implicit real*8 (a-h,o-z)
      external v
      parameter (nlevelmax=30,nrmax=5000)

      dimension psi(nlevelmax),psi0(nlevelmax),psi1(nlevelmax)
      dimension xi(nlevelmax,nlevelmax),xi0(nlevelmax,nlevelmax),
     &          xi1(nlevelmax,nlevelmax)
      dimension phi0(nlevelmax)
      dimension bb(nlevelmax,nlevelmax),
     &          bin(nlevelmax,nlevelmax)
      dimension bb2(nlevelmax,nlevelmax)
      dimension cc(nlevelmax,nlevelmax),
     &          cin(nlevelmax,nlevelmax)
      dimension dd0(nlevelmax,nlevelmax),dd1(nlevelmax,nlevelmax)
      dimension dd(nlevelmax)
      dimension fcw(0:200),gcw(0:200),fpcw(0:200),gpcw(0:200)
      dimension sigmad(0:200),iexp(0:200)
      dimension ech(nlevelmax),ech2(nlevelmax)
 
	common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi
	common/dimen/nlevel
	common/dyn/e,rmass
        common/angular/l
	common/grid/rmin,rmax,dr
        common/coup/cpot(nlevelmax,nlevelmax,0:nrmax)
	common/shape/rb,vb,curv
        
      complex*16 k,kk,k2
      complex*16 psi,psi0,psi1,phi0
      complex*16 xi,xi0,xi1
      complex*16 ai
      complex*16 cwup0,cwdown0,cwup1,cwdown1
      complex*16 bb,bb2,bin,cc,cin,dd
      complex*16 dummy,dummy2
      complex*16 xi1d(nlevelmax,nlevelmax)
      complex*16 xi0d(nlevelmax,nlevelmax)
      complex*16 aa(nlevelmax,nlevelmax)
      complex*16 bb0(nlevelmax,nlevelmax)
      complex*16 bb20(nlevelmax,nlevelmax)

      ai=(0.d0,1.d0)
c define two constants used in this subroutine
      fac=dr**2*(2.d0*rmass/hbar**2) 
      iterat=int((rmax-rmin)/dr)
      rmax=rmin+iterat*dr

      ibarrier=(rb-rmin)/dr+1.d-6

	do lc=0,200
	fcw(lc)=0.d0
	gcw(lc)=0.d0
	fpcw(lc)=0.d0
	gpcw(lc)=0.d0
	sigmad(lc)=0.d0
	iexp(lc)=0d0
        enddo

      do i=1,nlevelmax
      do j=1,nlevelmax
      bb(i,j)=0.d0
      bin(i,j)=0.d0
      cc(i,j)=0.d0
      cin(i,j)=0.d0
      aa(i,j)=0.d0
      enddo
      dd(i)=0.d0
      enddo

      do i=1,nlevel
         aa(i,i)=1.d0
         enddo

c  initial conditions
c==========================================================
      do 15 io=1,nlevel
c	write(6,*)'io=',io*1.

      do j1=1,nlevel
      psi(j1)=0.d0
      psi0(j1)=0.d0
      psi1(j1)=0.d0
      phi0(j1)=0.d0
      enddo

        if(io.eq.1) then
        do io2=1,nlevel
        ech(io2)=e-v(rmin)-cpot(io2,io2,0)
        ech2(io2)=e-cpot(io2,io2,iterat)
        enddo
        endif
        
	if(ech(io).gt.0.d0) then
        k=sqrt(2.d0*rmass/hbar**2*ech(io))
        psi0(io)=exp(-ai*k*rmin)
        phi0(io)=-ai*k*psi0(io)
	else
        k=sqrt(2.d0*rmass/hbar**2*abs(ech(io)))
        psi0(io)=exp(k*rmin)
        phi0(io)=k*psi0(io)
        endif

      call rkutta00(psi0,phi0,psi1)

      do i0=1,nlevel
	xi0(i0,io)=(1.d0-fac/12.d0*(v(rmin)-e))*psi0(i0)
	xi1(i0,io)=(1.d0-fac/12.d0*(v(rmin+dr)-e))*psi1(i0)
      do ic=1,nlevel
      xi0(i0,io)=xi0(i0,io)-fac/12.d0*cpot(i0,ic,0)*psi0(ic)
      xi1(i0,io)=xi1(i0,io)-fac/12.d0*cpot(i0,ic,1)*psi1(ic)
      enddo
      enddo

 15   continue

c===============
      do 29 ir=2,iterat+1
c----------------------------------------------  iterations start
      r=rmin+dr*ir
      r0=rmin+dr*(ir-2)
      r1=rmin+dr*(ir-1)

      do i0=1,nlevel
      do ic=1,nlevel
         dd0(i0,ic)=fac/sqrt(12.d0)*cpot(i0,ic,ir-1)
      if(i0.eq.ic) dd0(i0,ic)=dd0(i0,ic)+fac/sqrt(12.d0)*(v(r1)-e)
     &                        +sqrt(3.d0)
      enddo
      enddo

      do i0=1,nlevel
      do ic=1,nlevel
      dd1(i0,ic)=0.d0
      if(i0.eq.ic) dd1(i0,ic)=dd1(i0,ic)-1.d0
      do ik=1,nlevel
         dd1(i0,ic)=dd1(i0,ic)+dd0(i0,ik)*dd0(ik,ic)
         enddo
         enddo
         enddo
	
      do ich=1,nlevel
      do i0=1,nlevel
         xi(i0,ich)=-xi0(i0,ich)
      do ic=1,nlevel
         xi(i0,ich)=xi(i0,ich)+dd1(i0,ic)*xi1(ic,ich)
         enddo
         enddo
         enddo

      if(ir.eq.iterat+1) goto 66
      if(ir.eq.ibarrier) call stabilize(xi1,xi,aa,ir)

      do ich=1,nlevel
      do i0=1,nlevel
         xi0(i0,ich)=xi1(i0,ich)
         xi1(i0,ich)=xi(i0,ich)
      enddo
      enddo

 29   continue
c--------------------------------------------------------------
c  matching to the coulomb wave function at rmax

 66   do 17 io=1,nlevel

      do i0=1,nlevel
      do ic=1,nlevel
      cc(i0,ic)=-fac/12.d0*cpot(i0,ic,iterat-1)
      if(i0.eq.ic) cc(i0,ic)=cc(i0,ic)-fac/12.d0*(v(rmax-dr)-e)+1.d0
      enddo
      enddo

      call matinv(nlevel,cc,cin)
      do i0=1,nlevel
         psi0(i0)=0.d0
      do ic=1,nlevel
         psi0(i0)=psi0(i0)+cin(i0,ic)*xi0(ic,io)
      enddo
      enddo

      do i0=1,nlevel
      do ic=1,nlevel
      cc(i0,ic)=-fac/12.d0*cpot(i0,ic,iterat+1)
      if(i0.eq.ic) cc(i0,ic)=cc(i0,ic)-fac/12.d0*(v(rmax+dr)-e)+1.d0
      enddo
      enddo

      call matinv(nlevel,cc,cin)
      do i0=1,nlevel
         psi(i0)=0.d0
      do ic=1,nlevel
         psi(i0)=psi(i0)+cin(i0,ic)*xi(ic,io)
      enddo
      enddo

      do 60 ii=1,nlevel
c   coulomb wave function
      ec=e-cpot(ii,ii,iterat-1)

      if(ec.lt.0.) then
         r1=rmax-dr
         r2=rmax+dr
         ak=sqrt(2.d0*rmass*abs(ec)/hbar**2)
         
         bb0(ii,io)=(exp(-ak*r2)*psi0(ii)-exp(-ak*r1)*psi(ii))
     &               /(exp(ak*(r1-r2))-exp(-ak*(r1-r2)))
         bb20(ii,io)=-(exp(ak*r2)*psi0(ii)-exp(ak*r1)*psi(ii))
     &               /(exp(ak*(r1-r2))-exp(-ak*(r1-r2)))
         else
      rho=sqrt(2.d0*rmass*ec)/hbar*(rmax-dr)
      eta=(zp*zt/137.d0)*sqrt(rmass/2.d0/ec)
      call dfcoul(eta,rho,fcw,fpcw,gcw,gpcw,sigmad,l,iexp)
      cwup0=(gcw(l)+ai*fcw(l))
      cwdown0=gcw(l)-ai*fcw(l)

      ec=e-cpot(ii,ii,iterat+1)
      rho=sqrt(2.d0*rmass*ec)/hbar*(rmax+dr)
      eta=(zp*zt/137.d0)*sqrt(rmass/2.d0/ec)
      call dfcoul(eta,rho,fcw,fpcw,gcw,gpcw,sigmad,l,iexp)
      cwup1=(gcw(l)+ai*fcw(l))
      cwdown1=gcw(l)-ai*fcw(l)

      bb0(ii,io)=(cwup0*psi(ii)-cwup1*psi0(ii))
     &                         /(cwup0*cwdown1-cwup1*cwdown0)
      bb20(ii,io)=(cwdown1*psi0(ii)-cwdown0*psi(ii))
     &                         /(cwup0*cwdown1-cwup1*cwdown0)
      endif

 60   continue
c---------------------------------------------------------------
 17   continue

c===============================================================
c                                        penetration probability
      do i=1,nlevel
      do j=1,nlevel
         bb(i,j)=0.d0
         bb2(i,j)=0.d0
         do j2=1,nlevel
            bb(i,j)=bb(i,j)+bb0(i,j2)*aa(j2,j)
            bb2(i,j)=bb2(i,j)+bb20(i,j2)*aa(j2,j)
            enddo
            enddo
            enddo

c  calculate the inversion of the matrix bb
      call matinv(nlevel,bb,bin)
c
      p=0.d0
      ref=0.d0

      do 85 io=1,nlevel
      if(ech(io).lt.0.d0) goto 85
      k=sqrt((2.d0*rmass/hbar**2*ech(io)))
      kk=sqrt(2.d0*rmass/hbar**2*e)
      p=p+(abs(bin(io,1)))**2*k/kk
 85   continue

      do 86 io=1,nlevel
      if(ech2(io).lt.0.d0) goto 86
      k2=sqrt((2.d0*rmass/hbar**2*ech2(io)))
      kk=sqrt(2.d0*rmass/hbar**2*e)
             dummy=0.d0
             dummy2=0.d0
             do io2=1,nlevel
                dummy2=dummy2+bb2(io,io2)*bin(io2,1)
                enddo
      ref=ref+(abs(dummy2))**2*k2/kk
 86   continue

c      if(p+ref.gt.1.1d0) then
c         write(6,*)'penetrability > 1.1 !!!!!'
c         write(6,*)'something must be strange (accuracy problem) :^('
c         stop
c      endif

c      if(p.gt.0.99d0) p=1.d0-ref

      return
      end

c**************************************************************
      subroutine stabilize(xi1,xi,aa,ir1)
c
c a subroutine to stabilize the solutions by imposing the linear 
c independence. 
c
c refs.
c   z.h. levine, pra30('84) 1120
c   t.n. rescigno and a.e. orel, pra25('82)2402
c
c**************************************************************
      implicit real*8 (a-h,o-z)
      parameter (nlevelmax=30,nrmax=5000)

      dimension psi(nlevelmax,nlevelmax)
      dimension xi(nlevelmax,nlevelmax),xi1(nlevelmax,nlevelmax)
      dimension cc(nlevelmax,nlevelmax),
     &          cin(nlevelmax,nlevelmax)

	common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi
	common/dimen/nlevel
	common/dyn/e,rmass
        common/angular/l
	common/grid/rmin,rmax,dr
        common/coup/cpot(nlevelmax,nlevelmax,0:nrmax)

      complex*16 psi,xi,xi1
      complex*16 ai
      complex*16 cc,cin
      complex*16 aa(nlevelmax,nlevelmax)
      complex*16 aa0(nlevelmax,nlevelmax)
      complex*16 xid(nlevelmax,nlevelmax),xid1(nlevelmax,nlevelmax)

      do i=1,nlevel
      do j=1,nlevel
         aa0(i,j)=aa(i,j)
         xid(i,j)=xi(i,j)
         xid1(i,j)=xi1(i,j)
         enddo
         enddo

      ai=(0.d0,1.d0)
      fac=dr**2*(2.d0*rmass/hbar**2) 

c     transforme to the physical wave functions psi from xi 
      r=rmin+ir1*dr
         do i0=1,nlevel
         do ic=1,nlevel
          cc(i0,ic)=-fac/12.d0*cpot(i0,ic,ir1)
          if(i0.eq.ic) cc(i0,ic)=cc(i0,ic)-fac/12.d0*(v(r)-e)+1.d0
          enddo
          enddo
        call matinv(nlevel,cc,cin)
      do ich=1,nlevel
        do i0=1,nlevel
          psi(i0,ich)=0.d0
          do ic=1,nlevel
           psi(i0,ich)=psi(i0,ich)+cin(i0,ic)*xi(ic,ich)
           enddo
           enddo
           enddo

         do i=1,nlevel
         do j=1,nlevel
              aa(i,j)=0.d0
              do k=1,nlevel
                 aa(i,j)=aa(i,j)+psi(i,k)*aa0(k,j)
                 enddo
                 enddo
                 enddo

        call matinv(nlevel,psi,cin)
        
        do i=1,nlevel
        do j=1,nlevel
           xi(i,j)=0.d0
           xi1(i,j)=0.d0
           do k=1,nlevel
              xi(i,j)=xi(i,j)+xid(i,k)*cin(k,j)
              xi1(i,j)=xi1(i,j)+xid1(i,k)*cin(k,j)
              enddo
              enddo
              enddo

      return
      end

c**************************************************************
      subroutine rkutta00(psi0,phi0,psi1)
c
c wave functions at r=rmin+dr
c
c**************************************************************
      implicit real*8 (a-h,o-z)
      external v
      parameter (nlevelmax=30,nrmax=5000)

      dimension psi0(nlevelmax),phi0(nlevelmax),psi1(nlevelmax)
      dimension ak1(nlevelmax),ak2(nlevelmax),
     &          ak3(nlevelmax),ak4(nlevelmax)
      dimension bk1(nlevelmax),bk2(nlevelmax),
     &          bk3(nlevelmax),bk4(nlevelmax)
 
	common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi
	common/dimen/nlevel
	common/dyn/e,rmass
        common/angular/l
	common/grid/rmin,rmax,dr
        common/coup/cpot(nlevelmax,nlevelmax,0:nrmax)
        common/couph/cpoth(nlevelmax,nlevelmax)                  

      complex*16 psi0,phi0,psi1
      complex*16 ai
      complex*16 ak1,ak2,ak3,ak4
      complex*16 bk1,bk2,bk3,bk4

      ai=(0.d0,1.d0)
c define two constants used in this subroutine
      fac=dr*(2.d0*rmass/hbar**2) 

      r=rmin
      rpp=rmin+dr
      rh=rmin+dr/2.d0

      do 19 j1=1,nlevel
      ak1(j1)=0.d0
      ak2(j1)=0.d0
      ak3(j1)=0.d0
      ak4(j1)=0.d0
      bk1(j1)=0.d0
      bk2(j1)=0.d0
      bk3(j1)=0.d0
      bk4(j1)=0.d0
 19   continue

      do 58 i0=1,nlevel
      do 57 ic=1,nlevel
      ak1(i0)=ak1(i0)+fac*cpot(i0,ic,0)*psi0(ic)
 57   continue
      ak1(i0)=ak1(i0)-fac*(e-v(r))*psi0(i0)
      bk1(i0)=dr*phi0(i0)
 58   continue

      do 56 i0=1,nlevel
      do 55 ic=1,nlevel
      ak2(i0)=ak2(i0)+fac*cpoth(i0,ic)*(psi0(ic)+1.d0/2.d0*bk1(ic))
 55   continue
      ak2(i0)=ak2(i0)-fac*(e-v(rh))*(psi0(i0)+1.d0/2.d0*bk1(i0))
      bk2(i0)=dr*(phi0(i0)+1./2.*ak1(i0))
 56   continue

      do 54 i0=1,nlevel
      do 53 ic=1,nlevel
      ak3(i0)=ak3(i0)+fac*cpoth(i0,ic)*(psi0(ic)+1.d0/2.d0*bk2(ic))
 53   continue
      ak3(i0)=ak3(i0)-fac*(e-v(rh))*(psi0(i0)+1.d0/2.d0*bk2(i0))
      bk3(i0)=dr*(phi0(i0)+1.d0/2.d0*ak2(i0))
 54   continue

      do 42 i0=1,nlevel
      do 41 ic=1,nlevel
      ak4(i0)=ak4(i0)+fac*cpot(i0,ic,1)*(psi0(ic)+bk3(ic))
 41   continue
      ak4(i0)=ak4(i0)-fac*(e-v(rpp))*(psi0(i0)+bk3(i0))
      bk4(i0)=dr*(phi0(i0)+ak3(i0))
 42   continue

c wave functions at rmin+dr
      do 59 is=1,nlevel
      psi1(is)=psi0(is)+1.d0/6.d0*(bk1(is)+2.d0*bk2(is)
     &        +2.d0*bk3(is)+bk4(is))  
 59   continue
      return
      end

c******************************************************************
	subroutine cmat0
c
c preparation for the nuclear coupling matrix
c
c*****************************************************************
      implicit real*8 (a-h,o-z)
      parameter (nlevelmax=30)
      dimension a(nlevelmax,nlevelmax)
      common/hion/ap,zp,rp,at,zt,rt
      common/const/hbar,pi
      common/intr/ivibrotp,ivibrott,nproj,ntarg
      common/oscp/betap,betapn,omegap,lambdap,nphononp
      common/osct/betat,betatn,omegat,lambdat,nphonont
      common/osct2/betat2,betat2n,omegat2,lambdat2,nphonont2
      common/rotp/beta2p,beta4p,e2p,nrotp
      common/rott/beta2t,beta4t,e2t,nrott
      common/dimen/nlevel
      common/mut/imutual(0:nlevelmax,0:nlevelmax)
      common/eigen/ev(nlevelmax),evec(nlevelmax,nlevelmax)
      common/trans/ftr,qtrans,iq,ntrans
        common/ahv/betnahv(0:nlevelmax,0:nlevelmax),
     &             betcahv(0:nlevelmax,0:nlevelmax),
     &             omeahv(0:nlevelmax)
        common/ahv2/betnahv2(0:nlevelmax,0:nlevelmax),
     &             betcahv2(0:nlevelmax,0:nlevelmax),
     &             omeahv2(0:nlevelmax)
        common/ahvp/betnahvp(0:nlevelmax,0:nlevelmax),
     &             betcahvp(0:nlevelmax,0:nlevelmax),
     &             omeahvp(0:nlevelmax)
        common/grott/bett(0:nlevelmax,0:nlevelmax),erott(0:nlevelmax)
        common/grotp/betp(0:nlevelmax,0:nlevelmax),erotp(0:nlevelmax)

      external cg

        fvibt=rt*betatn/sqrt(4.d0*pi)
	frot2t=rt*beta2t
	frot4t=rt*beta4t
        fvibt2=rt*betat2n/sqrt(4.d0*pi)
        fvibp=rp*betapn/sqrt(4.d0*pi)
	frot2p=rp*beta2p
	frot4p=rp*beta4p

	do 10 i=1,nlevelmax
	do 20 j=1,nlevelmax
	a(i,j)=0.d0
 20	continue
 10	continue

	ndim=nlevel-ntrans
	if(ndim.eq.1) return
	
	   i=0
           do 31 ip=0,nproj
           do 32 it=0,ntarg
           do 33 it2=0,nphonont2
			if(nphonont2.ne.0) then
			if(imutual(it,it2).eq.0) goto 33
			endif
	   i=i+1
           j=0
           do 41 jp=0,nproj
           do 42 jt=0,ntarg
           do 43 jt2=0,nphonont2
			if(nphonont2.ne.0) then
			if(imutual(jt,jt2).eq.0) goto 43
			endif
	   j=j+1

              if(i.gt.j) then
                 a(i,j)=a(j,i)
                 goto 43
                 endif

	   c=0.d0
	   if(ip.eq.jp.and.it2.eq.jt2) then
             if(ivibrott.eq.0) then
c               if(it.eq.jt-1) c=sqrt(jt*1.d0)*fvibt
c               if(jt.eq.it-1) c=sqrt(it*1.d0)*fvibt

               c=rt*betnahv(it,jt)/sqrt(4.d0*pi)

	     else
c		c=frot2t*sqrt((2.d0*2*it+1.d0)*5.d0*(2.d0*2*jt+1.d0)
c     &    /4.d0/pi)*cg(2*it,0,2,0,2*jt,0)**2/(2.d0*2*jt+1.d0)
		c=rt*bett(it,jt)
     &             *sqrt((2.d0*2*it+1.d0)*5.d0*(2.d0*2*jt+1.d0)
     &             /4.d0/pi)*cg(2*it,0,2,0,2*jt,0)**2/(2.d0*2*jt+1.d0)
		c=c+frot4t*sqrt((2.d0*2*it+1.d0)*9.d0*(2.d0*2*jt+1.d0)
     &             /4.d0/pi)*cg(2*it,0,4,0,2*jt,0)**2/(2.d0*2*jt+1.d0)
            endif
	   endif
	   if(ip.eq.jp.and.it.eq.jt) then
c               if(it2.eq.jt2-1) c=c+sqrt(jt2*1.d0)*fvibt2
c               if(jt2.eq.it2-1) c=c+sqrt(it2*1.d0)*fvibt2

               c=c+rt*betnahv2(it2,jt2)/sqrt(4.d0*pi)

	   endif

c-------
	   if(it.eq.jt.and.it2.eq.jt2) then
             if(ivibrotp.eq.0) then
c               if(ip.eq.jp-1) c=c+sqrt(jp*1.d0)*fvibp
c               if(jp.eq.ip-1) c=c+sqrt(ip*1.d0)*fvibp

               c=c+rp*betnahvp(ip,jp)/sqrt(4.d0*pi)

	     else
c		c=c+frot2p*sqrt((2.d0*2*ip+1.d0)*5.d0*(2.d0*2*jp+1.d0)
c     &     /4.d0/pi)*cg(2*ip,0,2,0,2*jp,0)**2/(2.d0*2*jp+1.d0)
		c=c+rp*betp(ip,jp)
     &              *sqrt((2.d0*2*ip+1.d0)*5.d0*(2.d0*2*jp+1.d0)
     &              /4.d0/pi)*cg(2*ip,0,2,0,2*jp,0)**2/(2.d0*2*jp+1.d0)
		c=c+frot4p*sqrt((2.d0*2*ip+1.d0)*9.d0*(2.d0*2*jp+1.d0)
     &              /4.d0/pi)*cg(2*ip,0,4,0,2*jp,0)**2/(2.d0*2*jp+1.d0)
             endif
            endif


	a(i,j)=c

 43              continue
 42              continue
 41              continue
 33              continue
 32              continue
 31              continue

	call mdiag(a,ndim)

      return
      end
c******************************************************************
	subroutine cmat(r,cpot)
c
c coupling matrix
c
c*****************************************************************
      implicit real*8 (a-h,o-z)
      parameter (nlevelmax=30)
      dimension cpot(nlevelmax,nlevelmax)
      common/hion/ap,zp,rp,at,zt,rt
      common/const/hbar,pi
      common/intr/ivibrotp,ivibrott,nproj,ntarg
      common/oscp/betap,betapn,omegap,lambdap,nphononp
      common/osct/betat,betatn,omegat,lambdat,nphonont
      common/osct2/betat2,betat2n,omegat2,lambdat2,nphonont2
      common/rotp/beta2p,beta4p,e2p,nrotp
      common/rott/beta2t,beta4t,e2t,nrott
      common/dimen/nlevel
      common/mut/imutual(0:nlevelmax,0:nlevelmax)
      common/eigen/ev(nlevelmax),evec(nlevelmax,nlevelmax)
      common/trans/ftr,qtrans,iq,ntrans
        common/ahv/betnahv(0:nlevelmax,0:nlevelmax),
     &             betcahv(0:nlevelmax,0:nlevelmax),
     &             omeahv(0:nlevelmax)
        common/ahv2/betnahv2(0:nlevelmax,0:nlevelmax),
     &             betcahv2(0:nlevelmax,0:nlevelmax),
     &             omeahv2(0:nlevelmax)
        common/ahvp/betnahvp(0:nlevelmax,0:nlevelmax),
     &             betcahvp(0:nlevelmax,0:nlevelmax),
     &             omeahvp(0:nlevelmax)
        common/grott/bett(0:nlevelmax,0:nlevelmax),erott(0:nlevelmax)
        common/grotp/betp(0:nlevelmax,0:nlevelmax),erotp(0:nlevelmax)
      external fct,fcp,fact,vncc,cg,fct2,fct4,fcp2,fcp4,fctt
      external fct2v,fcp2v
      external ftrans

	do 10 i=1,nlevelmax
	do 20 j=1,nlevelmax
	cpot(i,j)=0.d0
 20	continue
 10	continue

	ndim=nlevel-ntrans
	if(nlevel.eq.1) return
	if(ndim.eq.1) goto 55
	
	   i=0
           do 31 ip=0,nproj
           do 32 it=0,ntarg
           do 33 it2=0,nphonont2
			if(nphonont2.ne.0) then
			if(imutual(it,it2).eq.0) goto 33
			endif
	   i=i+1
           j=0
           do 41 jp=0,nproj
           do 42 jt=0,ntarg
           do 43 jt2=0,nphonont2
			if(nphonont2.ne.0) then
			if(imutual(jt,jt2).eq.0) goto 43
			endif
	   j=j+1

              if(i.gt.j) then
                 cpot(i,j)=cpot(j,i)
                 goto 43
                 endif
	c=0.d0
c nuclear coupling

	do 50 k=1,ndim
	c=c+vncc(r,ev(k))*evec(i,k)*evec(j,k)
 50	continue

c coulomb coupling 

	   if(ip.eq.jp.and.it2.eq.jt2) then
             if(ivibrott.eq.0) then
c               if(it.eq.jt-1) c=c+sqrt(jt*1.d0)*fct(r)
c               if(jt.eq.it-1) c=c+sqrt(it*1.d0)*fct(r)

               if(it.ne.jt) then
               cc=betcahv(it,jt)*fct(r)
               if(betat.ne.0.d0) cc=cc/betat
               else
               cc=betcahv(it,jt)*fct2v(r)/sqrt(4.d0*pi)
               endif

               c=c+cc

	     else

        	bet=bett(it,jt)+2.d0*sqrt(5.d0/pi)*bett(it,jt)**2/7.d0
c		c=c+sqrt((2.d0*2*it+1.d0)*5.d0*(2.d0*2*jt+1.d0)
c     &     /4.d0/pi)*cg(2*it,0,2,0,2*jt,0)**2/(2.d0*2*jt+1.d0)*fct2(r)
		c=c+sqrt((2.d0*2*it+1.d0)*5.d0*(2.d0*2*jt+1.d0)
     &     /4.d0/pi)*cg(2*it,0,2,0,2*jt,0)**2/(2.d0*2*jt+1.d0)
     &              *fct2(r)*bet
     &		   +sqrt((2.d0*2*it+1.d0)*9.d0*(2.d0*2*jt+1.d0)
     &     /4.d0/pi)*cg(2*it,0,4,0,2*jt,0)**2/(2.d0*2*jt+1.d0)*fct4(r)
            endif
	   endif

	   if(it.eq.jt.and.it2.eq.jt2) then
             if(ivibrotp.eq.0) then
c               if(ip.eq.jp-1) c=c+sqrt(jp*1.d0)*fcp(r)
c               if(jp.eq.ip-1) c=c+sqrt(ip*1.d0)*fcp(r)

               if(ip.ne.jp) then
               cc=betcahvp(ip,jp)*fcp(r)
               if(betap.ne.0.d0) cc=cc/betap
               else
               cc=betcahvp(ip,jp)*fcp2v(r)/sqrt(4.d0*pi)
               endif
               c=c+cc

	     else
        	bet=betp(ip,jp)+2.d0*sqrt(5.d0/pi)*betp(ip,jp)**2/7.d0
c		c=c+sqrt((2.d0*2*ip+1.d0)*5.d0*(2.d0*2*jp+1.d0)
c     &     /4.d0/pi)*cg(2*ip,0,2,0,2*jp,0)**2/(2.d0*2*jp+1.d0)*fcp2(r)
		c=c+sqrt((2.d0*2*ip+1.d0)*5.d0*(2.d0*2*jp+1.d0)
     &     /4.d0/pi)*cg(2*ip,0,2,0,2*jp,0)**2/(2.d0*2*jp+1.d0)
     &          *fcp2(r)*bet
     &            +sqrt((2.d0*2*ip+1.d0)*9.d0*(2.d0*2*jp+1.d0)
     &     /4.d0/pi)*cg(2*ip,0,4,0,2*jp,0)**2/(2.d0*2*jp+1.d0)*fcp4(r)
             endif
            endif

	   if(ip.eq.jp.and.it.eq.jt) then
c               if(it2.eq.jt2-1) c=c+sqrt(jt2*1.d0)*fctt(r)
c               if(jt2.eq.it2-1) c=c+sqrt(it2*1.d0)*fctt(r)

               if(it2.ne.jt2) then
               cc=betcahv2(it2,jt2)*fctt(r)
               if(betat2.ne.0.d0) cc=cc/betat2
               else
               cc=betcahv2(it2,jt2)*fct2v(r)/sqrt(4.d0*pi)
               endif
               c=c+cc

	   endif

c excitation energy
           if(it.eq.jt.and.ip.eq.jp.and.it2.eq.jt2) then
		if(ivibrott.eq.0) then
c		c=c+it*omegat
		c=c+omeahv(it)
		else
c		c=c+2.d0*it*(2.d0*it+1.d0)/6.d0*e2t
		c=c+erott(it)
		endif
		if(ivibrotp.eq.0) then
c		c=c+ip*omegap
		c=c+omeahvp(ip)
		else
c		c=c+2.d0*ip*(2.d0*ip+1.d0)/6.d0*e2p
		c=c+erotp(ip)
		endif
c	        c=c+it2*omegat2
	        c=c+omeahv2(it2)
 	   endif

	cpot(i,j)=c

 43              continue
 42              continue
 41              continue
 33              continue
 32              continue
 31              continue

      c0=cpot(1,1)
      do 51 i=1,ndim
          cpot(i,i)=cpot(i,i)-vn(r)
c         cpot(i,i)=cpot(i,i)-c0
 51      continue

c    transfer channel
 55      if(ntrans.eq.1) then
            cpot(1,nlevel)=ftrans(r)
            cpot(nlevel,1)=ftrans(r)
	    cpot(nlevel,nlevel)=cpot(1,1)-qtrans
     &          +(zp+iq)*(zt-iq)/r*hbar/137.d0-zp*zt/r*hbar/137.d0
            endif

      return
      end

c*****************************************************************
	function v(r)
c
c  potential for the relative motion
c
c*****************************************************************
	implicit real*8(a-h,o-z)
	external vc,vn,vcent
	v=vc(r)+vn(r)+vcent(r)
	return
	end
c*****************************************************************
	function vc(r)
c
c coulomb potential
c
c*****************************************************************
	implicit real*8(a-h,o-z)
	common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi
c  coulomb radius
c	rc=1.2d0*(ap**(1.d0/3.d0)+at**(1.d0/3.d0))
c	if(r.gt.rc) then
	vc=zp*zt/r*hbar/137.d0
c	else
c	vc=zp*zt*hbar/137.d0*(3.d0*rc**2-r**2)/2.d0/rc**3
c	endif
	return
	end
c*****************************************************************
	function dvc(r)
c
c first derivative of the coulomb potential
c
c*****************************************************************
	implicit real*8(a-h,o-z)
	common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi
c  coulomb radius
c	rc=1.2d0*(ap**(1.d0/3.d0)+at**(1.d0/3.d0))
c	if(r.gt.rc) then
	dvc=-zp*zt/r**2*hbar/137.d0
c	else
c	vc=zp*zt*hbar/137.d0*(3.d0*rc**2-r**2)/2.d0/rc**3
c	endif
	return
	end
c****************************************************************
	function vn(r)
c
c (power of) woods-saxon potential for nuclear potential
c
c****************************************************************
	implicit real*8 (a-h,o-z)
	common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi
        common/pot/v0,r0,a0,pow

	if(r.lt.50.d0) then
        vn=-v0/(1.d0+exp((r-r0)/a0/pow))**pow
	else
	vn=0.d0
        endif

	return
	end
c****************************************************************
	function vncc(r,xt)
c
c (power of) woods-saxon potential for nuclear potential 
c
c****************************************************************
	implicit real*8 (a-h,o-z)

	common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi
        common/pot/v0,r0,a0,pow

	if(r.lt.50.d0) then
           rr=r-r0-xt
        vncc=-v0/(1.d0+exp(rr/a0/pow))**pow
	else
	vncc=0.d0
        endif

	return
	end

c***************************************************************
      function dvn(r)
c     
c the first derivative of vn(r)
c
c**************************************************************
      implicit real*8(a-h,o-z)
      common/pot/v0,r0,a,pow
      external vn

      if(r.lt.50.d0) then
c      dvn=v0/a*exp((r-r0)/a)/(1.d0+exp((r-r0)/a))**2
      dvn=(vn(r+1.d-5)-vn(r-1.d-5))/2.d-5
      else
      dvn=0.d0
      endif

      return
      end
c************************************************************
	function vcent(r)
c
c centrifugal potential
c
c************************************************************
	implicit real*8(a-h,o-z)
	common/const/hbar,pi
	common/angular/l
	common/dyn/e,rmass
	vcent=l*(l+1.d0)*hbar**2/2.d0/rmass/r**2
	return
	end
c************************************************************
	function dvcent(r)
c
c the first derivative of the centrifugal potential
c
c************************************************************
	implicit real*8(a-h,o-z)
	common/const/hbar,pi
	common/angular/l
	common/dyn/e,rmass
	dvcent=-2.d0*l*(l+1.d0)*hbar**2/2.d0/rmass/r**3
	return
	end

c***************************************************************
	function dv(r)
c
c the first derivative of v(r)
c
c***************************************************************
	implicit real*8(a-h,o-z)
	external dvn,dvc,dvcent

	dv=dvn(r)+dvc(r)+dvcent(r)

	return
	end
c*****************************************************************
	function fcp(r)
c
c coulomb coupling form factor for projectile phonon excitation
c
c*****************************************************************
	implicit real*8(a-h,o-z)
        common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi
        common/oscp/betap,betapn,omega,lambda,nphonon

	rc=rp
	if(r.gt.rc) then
	fcp=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/r
     &                                         *(rp/r)**lambda
	else
	fcp=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/rc
     &                                         *(r/rp)**lambda
	endif

        fcp=fcp*betap/sqrt(4.d0*pi)

	return
	end

c*****************************************************************
	function fcp2v(r)
c
c coulomb coupling form factor for projectile excitation
c (for anharmonic vibration)
c
c*****************************************************************
	implicit real*8(a-h,o-z)
        common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi

	lambda=2
	rc=rp
	if(r.gt.rc) then
	fcp2v=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/r
     &                                         *(rp/r)**lambda
	else
	fcp2v=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/rc
     &                                         *(r/rp)**lambda
	endif

	return
	end
c*****************************************************************
	function fcp2(r)
c
c coulomb coupling form factor for projectile excitation
c (rotational e2 coupling)
c
c*****************************************************************
	implicit real*8(a-h,o-z)
        common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi
	common/rotp/beta2p,beta4p,e2p,nrotp

	lambda=2
	rc=rp
	if(r.gt.rc) then
	fcp2=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/r
     &                                         *(rp/r)**lambda
	else
	fcp2=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/rc
     &                                         *(r/rp)**lambda
	endif

c	fcp2=fcp2*(beta2p+2.d0*sqrt(5.d0/pi)*beta2p**2/7.d0)

	return
	end
c*****************************************************************
	function fcp4(r)
c
c coulomb coupling form factor for projectile excitation
c (rotational e4 coupling)
c
c*****************************************************************
	implicit real*8(a-h,o-z)
        common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi
	common/rotp/beta2p,beta4p,e2p,nrotp

	lambda=4
	rc=rp
	if(r.gt.rc) then
	fcp4=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/r
     &                                         *(rp/r)**lambda
	else
	fcp4=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/rc
     &                                         *(r/rp)**lambda
	endif

	fcp4=fcp4*(beta4p+9.d0*beta2p**2/7.d0/sqrt(pi))

	return
	end
c*****************************************************************
	function fct(r)
c
c coulomb coupling form factor for the target phonon excitation
c
c*****************************************************************
	implicit real*8(a-h,o-z)
        common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi
        common/osct/betat,betatn,omega,lambda,nphonon

	rc=rt
	if(r.gt.rc) then
	fct=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/r
     &                                         *(rt/r)**lambda
	else
	fct=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/rc
     &                                         *(r/rt)**lambda
	endif

        fct=fct*betat/sqrt(4.d0*pi)

	return
	end

c*****************************************************************
	function fctt(r)
c
c coulomb coupling form factor for the second target phonon excitation
c
c*****************************************************************
	implicit real*8(a-h,o-z)
        common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi
        common/osct2/betat2,betat2n,omega2,lambda2,nphonon2

	rc=rt
	if(r.gt.rc) then
	fctt=3.d0/(2.d0*lambda2+1.d0)*zp*zt/137.d0*hbar/r
     &                                         *(rt/r)**lambda2
	else
	fctt=3.d0/(2.d0*lambda2+1.d0)*zp*zt/137.d0*hbar/rc
     &                                         *(r/rt)**lambda2
	endif

        fctt=fctt*betat2/sqrt(4.d0*pi)

	return
	end
c*****************************************************************
	function fct2v(r)
c
c coulomb coupling form factor for target excitation
c (for anharmonic vibration)
c
c*****************************************************************
	implicit real*8(a-h,o-z)
        common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi

	lambda=2
	rc=rt
	if(r.gt.rc) then
	fct2v=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/r
     &                                         *(rt/r)**lambda
	else
	fct2v=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/rc
     &                                         *(r/rt)**lambda
	endif

	return
	end


c*****************************************************************
	function fct2(r)
c
c coulomb coupling form factor for target excitation
c (rotational e2 coupling)
c
c*****************************************************************
	implicit real*8(a-h,o-z)
        common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi
	common/rott/beta2t,beta4t,e2t,nrott

	lambda=2
	rc=rt
	if(r.gt.rc) then
	fct2=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/r
     &                                         *(rt/r)**lambda
	else
	fct2=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/rc
     &                                         *(r/rt)**lambda
	endif

c	fct2=fct2*(beta2t+2.d0*sqrt(5.d0/pi)*beta2t**2/7.d0)

	return
	end

c*****************************************************************
	function fct4(r)
c
c coulomb coupling form factor for target excitation
c (rotational e4 coupling)
c
c*****************************************************************
	implicit real*8(a-h,o-z)
        common/hion/ap,zp,rp,at,zt,rt
	common/const/hbar,pi
	common/rott/beta2t,beta4t,e2t,nrott

	lambda=4
	rc=rt
	if(r.gt.rc) then
	fct4=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/r
     &                                         *(rt/r)**lambda
	else
	fct4=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/rc
     &                                         *(r/rt)**lambda
	endif

	fct4=fct4*(beta4t+9.d0*beta2t**2/7.d0/sqrt(pi))

	return
	end

c***************************************************************
      function ddvn(r)
c     
c the second derivative of vn(r)
c
c**************************************************************
      implicit real*8(a-h,o-z)
      common/pot/v0,r0,a,pow
      external dvn

      if(r.lt.50.d0) then
c      ddvn=v0/a**2*exp((r-r0)/a)/(1.d0+exp((r-r0)/a))**2
c     &    -2.d0*v0/a**2*exp(2.d0*(r-r0)/a)
c     &                           /(1.d0+exp((r-r0)/a))**3
      ddvn=(dvn(r+1.d-5)-dvn(r-1.d-5))/2.d-5
      else
      ddvn=0.d0
      endif

      return
      end
c****************************************************************
	function ftrans(r)
c
c coupling form factor for transfer reactions
c
c***************************************************************
	implicit real*8(a-h,o-z)
	common/const/hbar,pi
        common/trans/ftr,qtrans,iq,ntrans
	external dvn

	ftrans=ftr*dvn(r)
	return
	end

c*************************************************************
	subroutine mdiag(a,n)
c
c the jacobi method to compute all eigenvalues and eigenvectors 
c of a real symmetric matrix a, which is of size n by n, stored 
c in a physical np by np array. on output, the elements of a 
c above the diagonal are destroyed. d returns the eigenvalues of 
c a in its first n elements. v is a matrix with the same logical 
c and physical dimensions as a whose columns contain, on output, 
c the normalized eigenvectors of a. nrot returns the number of 
c jacobi rotations which were required. 
c
c n.b.; the i-th component of the eigenvector for the k-th 
c       eigenvalue is given by v(i,k). 
c
c************************************************************
	implicit real*8(a-h,o-z)
	parameter(nmax=100)
        parameter (nlevelmax=30)
        dimension a(nlevelmax,nlevelmax)
	dimension b(nmax),z(nmax)
        common/eigen/d(nlevelmax),v(nlevelmax,nlevelmax)

	do 12 ip=1,n
	do 11 iq=1,n
	v(ip,iq)=0.d0
 11	continue
	v(ip,ip)=1.d0
 12	continue

	do 13 ip=1,n
	b(ip)=a(ip,ip)
	d(ip)=b(ip)
	z(ip)=0.d0
 13	continue

	nrot=0

	do 24 i=1,50
	sm=0.d0

	do 15 ip=1,n-1
	do 14 iq=ip+1,n
	sm=sm+abs(a(ip,iq))
 14	continue
 15     continue

	if(sm.eq.0.d0) return

	if(i.lt.4) then
	tresh=0.2d0*sm/n**2
	else
	tresh=0.d0
	endif

	do 22 ip=1,n-1
	do 21 iq=ip+1,n
	g=100.d0*abs(a(ip,iq))
	
	if((i.gt.4).and.(abs(d(ip))+g.eq.abs(d(ip)))
     &     .and.(abs(d(iq))+g.eq.abs(d(iq)))) then
	   a(ip,iq)=0.d0
	elseif(abs(a(ip,iq)).gt.tresh) then
	   h=d(iq)-d(ip)

	   if(abs(h)+g.eq.abs(h)) then
	     t=a(ip,iq)/h
	   else
	     theta=0.5d0*h/a(ip,iq)
	     t=1.d0/(abs(theta)+sqrt(1.d0+theta**2))
	     if(theta.lt.0.d0) t=-t
	   endif

	c=1.d0/sqrt(1.d0+t**2)
	s=t*c
	tau=s/(1.d0+c)
	h=t*a(ip,iq)
	z(ip)=z(ip)-h
	z(iq)=z(iq)+h
	d(ip)=d(ip)-h
	d(iq)=d(iq)+h
	a(ip,iq)=0.d0

	do 16 j=1,ip-1
	g=a(j,ip)
	h=a(j,iq)
	a(j,ip)=g-s*(h+g*tau)
	a(j,iq)=h+s*(g-h*tau)
 16	continue

	do 17 j=ip+1,iq-1
	g=a(ip,j)
	h=a(j,iq)
	a(ip,j)=g-s*(h+g*tau)
	a(j,iq)=h+s*(g-h*tau)
 17	continue

	do 18 j=iq+1,n
	g=a(ip,j)
	h=a(iq,j)
	a(ip,j)=g-s*(h+g*tau)
	a(iq,j)=h+s*(g-h*tau)
 18	continue

	do 19 j=1,n
	g=v(j,ip)
	h=v(j,iq)
	v(j,ip)=g-s*(h+g*tau)
	v(j,iq)=h+s*(g-h*tau)
 19	continue
	
	nrot=nrot+1
	endif

 21	continue
 22	continue

	do 23 ip=1,n
	b(ip)=b(ip)+z(ip)
	d(ip)=b(ip)
	z(ip)=0.d0

 23	continue
 24	continue

	pause '50 iterations should never happen'

	return
	end

c****************************************************************
      function cg(j1,m1,j2,m2,j3,m3)
c
c  clebsh-gordan coefficient <j1 m1 j2 m2 | j3 m3>
c****************************************************************
      implicit real*8(a-h,o-z)
	external fact

	if(m1+m2.ne.m3) then
	cg=0.d0
	return
	endif

	if(j3.lt.abs(j1-j2)) then
	cg=0.d0
	return
	endif

	if(j3.gt.j1+j2) then
	cg=0.d0
	return
	endif

      ka=j1+j2-j3
      kb=j3+j1-j2
      kc=j2+j3-j1
      kd=j1+j2+j3+1
      del=sqrt(fact(ka)*fact(kb)*fact(kc)/fact(kd))

      cg=0.d0

      do 10 n=0,max(j1+j2-j3,j1-m1,j2+m2)

      ka1=j1+j2-j3-n
      if(ka1.lt.0.d0) go to 10
      ka2=j3-j2+m1+n
      if(ka2.lt.0.d0) go to 10
      ka3=j3-j1-m2+n
      if(ka3.lt.0.d0) go to 10
      ka4=j1-m1-n
      if(ka4.lt.0.d0) go to 10
      ka5=j2+m2-n
      if(ka5.lt.0.d0) go to 10
      an=n*1.d0

      cg=cg+(-1.d0)**n
     &  /(fact(n)*fact(ka1)*fact(ka2)*fact(ka3)*fact(ka4)*fact(ka5))

 10   continue

 20   cg=cg*sqrt(fact(j1+m1)*fact(j1-m1))
      cg=cg*sqrt(fact(j2+m2)*fact(j2-m2))
      cg=cg*sqrt(fact(j3+m3)*fact(j3-m3))

      cg=cg*sqrt(2.d0*j3+1.d0)*del

      return
      end
c***************************************************************
      subroutine matinv(nmax,c,d)
c     
c subroutine for calculating the inversion of a matrix c, 
c whose dimension is nmax
c
c***************************************************************
      implicit complex*16 (a-h,o-z)
      parameter (nlevelmax=30)
      dimension u(nlevelmax),v(nlevelmax),
     &          c(nlevelmax,nlevelmax),d(nlevelmax,nlevelmax)
      deter=1.0
      do 1 m=1,nmax
      do 1 n=1,nmax
      d(n,m)=0.0
      if(n.eq.m) d(n,m)=1.0
 1    continue
      do 10 n=1,nmax
      t=c(n,n)
      if (abs(t).lt.1.e-10) go to3
      go to 6
 3    j=n
 4    j=j+1
      if (j.gt.nmax) go to 11
      t=c(n,j)
      deter=-deter
      if (abs(t).le.1.e-10) go to 4
      do 5 k=1,nmax
      u(k)=c(n,k)
      v(k)=d(n,k)
      c(n,k)=c(j,k)
      d(n,k)=d(j,k)
      c(j,k)=u(k)
 5    d(j,k)=v(k)
 6    do 7 k=1,nmax
      if (k.eq.n) go to 7
      a=c(k,n)/c(n,n)
      do 8 l=1,nmax
      c(k,l)=c(k,l)-a*c(n,l)                              
      d(k,l)=d(k,l)-a*d(n,l)
 8    continue
 7    continue
      b=c(n,n)
      deter=b*deter
      do 10 m=1,nmax
      c(n,m)=c(n,m)/b
      d(n,m)=d(n,m)/b
 10   continue
      return
 11   print 12
 12   format('matrix not invertible')
      return
      end

                         function fact(n)
                         implicit real*8(a-h,o-z)
			 if(n.lt.0) then
			 fact=0.d0
			 return
			 endif
                         if (n.eq.0) then
                            fact=1.d0
                            go to 10
                            endif
                            fact=1.d0
                            do 20 i=1,n
                               fact=fact*i*1.d0
 20                            continue
 10                            return
                               end

C*****************************************************************
C Subroutine for the Coulomb wave function
C*****************************************************************
      SUBROUTINE JFLGAM(XD,YD,PAR,PAI,NBCHIF)
      DOUBLE PRECISION XD,YD,PAR,PAI,TEST,C,HLO2PI,PI,PIS2,PIS4
      DOUBLE PRECISION X,Y,U,V,TRA,TRA1,ALO2PI,RAC2,COSI,SINI
      DOUBLE PRECISION COS2I,SIN2I,ZMOD,DEPI
      DOUBLE PRECISION XX
      DOUBLE PRECISION ALOPI
      DOUBLE PRECISION SUPINT
      DIMENSION TEST(7),C(6)
      DATA TEST/2.9152D+7,2.2958D+3,1.4124D+2,3.9522D+1,19.6611D0,
     &12.791D0,-10.D0/
      DATA C/8.333333333333333D-2,-2.777777777777777D-3,
     &7.936507936507937D-4,-5.952380952380952D-4,8.417508417508418D-4,
     &-1.917526917526918D-3/
      DATA HLO2PI/0.918938533204672/,PI/3.141592653589793/
      DATA PIS2/1.570796326794897/,PIS4/0.785398163397448/
      DATA ALO2PI/1.837877066409345/,RAC2/0.3465735902799726/
      DATA DEPI/6.283185307179586/,ALOPI/1.1447298858494002/
      DATA SUPINT/2147483647.D0/
      NBCHIF=15
      X=DABS(XD)
      XX=X
      IF(YD)1,2,1
    1 Y=DABS(YD)
      KR=1
      I=DMOD(10.99D0-X,SUPINT)
C     TRANSLATION
      IF(I)3,3,4
    4 TRA=I
      X=X+TRA
C     LOGARITHME(X+IY) (X,Y POSITIFS)
    3 IF(X-Y)5,6,7
    5 TRA1=X/Y
      IF(TRA1)8,8,9
    8 U=DLOG(Y)
      V=PIS2
      GO TO 10
    6 U=RAC2+DLOG(X)
      V=PIS4
      GO TO 10
    9 TRA=Y*DSQRT(1.D0+TRA1*TRA1)
      TRA1=Y/X
   11 U=DLOG(TRA)
      V=DATAN(TRA1)
   10 GO TO (12,19,23),KR
    7 TRA1=Y/X
      TRA=X*DSQRT(1.D0+TRA1*TRA1)
      GO TO 11
   12 KR=2
C     DEVELOPPEMENT ASYMPTOTIQUE ( X SUPERIEUR A 10 )
      TRA=X-0.5D0
      PAI=V*TRA+Y*(U-1.D0)
      PAR=-X+HLO2PI+U*TRA-V*Y
      ZMOD=X+Y
      IF(ZMOD-TEST(1))13,13,14
   13 TRA=X*X+Y*Y
      COSI=X/TRA
      SINI=Y/TRA
      SIN2I=(SINI*COSI)+(SINI*COSI)
      COS2I=(COSI+SINI)*(COSI-SINI)
      K=1
      GO TO 15
   16 TRA=COSI*COS2I-SINI*SIN2I
      SINI=SINI*COS2I+COSI*SIN2I
      COSI=TRA
   15 PAR=PAR+C(K)*COSI
      PAI=PAI-C(K)*SINI
      K=K+1
      IF(ZMOD-TEST(K))16,16,14
C     TRANSLATION INVERSE
   17 I=I-1
      X=I
      X=XX+X
      GO TO 3
   19 PAR=PAR-U
      PAI=PAI-V
   14 IF(I-1)18,60,17
   60 IF(XD)17,61,17
C     CONTROLE DU QUADRANT
   18 IF(XD)20,61,21
   61 TRA=PI*Y
      IF(TRA-1.D-2)300,300,301
  300 PAR= TRA*(2.D0+TRA*(-2.D0+TRA*(1.333333333333333D0+TRA*(
     &-0.6666666666666666D0+TRA*(0.2666666666666666D0+TRA*(
     &-0.08888888888888888D0+TRA*0.02539682539682540D0))))))
      TRA1=-DLOG(Y)-DLOG(PAR)
      GO TO 302
  301 PAR=1.D0-DEXP(-TRA-TRA)
      TRA1=-DLOG(Y*PAR)
  302 PAR=0.5D0*(ALO2PI-TRA+TRA1)
      PAI=PAI-PIS2
   21 IF(YD)28,100,100
C     X+IY CHANGE EN -X-IY
   20 TRA=PI*Y
      PAR=ALO2PI-U-PAR-TRA
      PAI=PI-V-PAI
      TRA=DEXP(-TRA-TRA)
      X=PI*DMOD(X,2.D0)
      SINI=(1.D0-TRA)*DCOS(X)
      COSI=(1.D0+TRA)*DSIN(X)
      KR=3
      X=DABS(COSI)
      Y=DABS(SINI)
      GO TO 3
   23 IF(COSI)24,25,25
   24 V=PI-DSIGN(V,SINI)
      GO TO 26
   25 IF(SINI)27,26,26
   27 V=-V
   26 PAR=PAR-U
      PAI=PAI-V
      IF(YD)100,100,28
   28 PAI=-PAI
C     ARGUMENT DANS -PI,PI
  100 TRA=DABS(PAI/DEPI)
      IF(TRA-1.D+15)203,204,204
  204 NBCHIF=0
      PAI=0.D0
      GO TO 29
  203 IF(TRA-1.D0)205,205,206
  206 NBCHIF=DLOG10(TRA)
      NBCHIF=14-NBCHIF
      TRA=DMOD(TRA,SUPINT)
      PAI=DMOD(TRA,1.D0)*DSIGN(DEPI,PAI)
  205 IF(DABS(PAI)-PI)29,29,207
  207 PAI=PAI-DSIGN(DEPI,PAI)
      GO TO 29
C     JFLGAM REEL
    2 PAI=0.D0
      IF(XD)31,32,33
C     CONDITIONS D EXISTENCE
   32 WRITE (6,1000)
 1000 FORMAT (21H JFLGAM(0) EST INFINI)
      GO TO 50
   31 IF(X-4503599627370496.D0)34,35,35
   35 WRITE (6,1001)
 1001 FORMAT (30H ARGUMENT DE JFLGAM TROP GRAND)
      GO TO 50
   34 Y=DMOD(X,SUPINT)
      IF(Y)400,36,400
  400 IF(Y-0.99D0)33,33,405
  405 TRA=IDINT(Y+0.1D0)
      IF(DABS(Y-TRA)-5.D-15)36,36,33
   36 WRITE (6,1002)
 1002 FORMAT (28H JFLGAM (-ENTIER) EST INFINI)
   50 PAR=1.D+74
      NBCHIF=0
      GO TO 29
C     TRANSLATION
   33 I=DMOD(10.99D0-X,SUPINT)
      IF(I)37,37,38
   38 TRA=I
      X=X+TRA
C     DEVELOPPEMENT ASYMPTOTIQUE
   37 Y=DLOG(X)
      PAR=-X+HLO2PI+Y*(X-0.5D0)
      IF(X-TEST(1))39,39,43
   39 COSI=1.D0/X
      COS2I=COSI*COSI
      K=1
      GO TO 41
   42 COSI=COSI*COS2I
   41 PAR=PAR+C(K)*COSI
      K=K+1
      IF(X-TEST(K))42,42,40
C     TRANSLATION INVERSE
   44 X=X-1.D0
   48 Y=DLOG(X)
      PAR=PAR-Y
      I=I-1
   40 IF(I-1)43,49,44
   49 X=DABS(XD)
      GO TO 48
C     X NEGATIF
   43 IF(XD)45,29,29
   45 PAR=ALOPI-PAR-Y
      Y=PI*DMOD(X,2.D0)
      Y=-DSIN(Y)
      IF(Y)46,36,47
   46 Y=-Y
      PAI=PI
   47 PAR=PAR-DLOG(Y)
      ENTRY JFLGV1
   29 RETURN
      END

      SUBROUTINE YFCLEN(ETA,RO,U,UP,V,VP,SIGMA0,IDIV,NN)
      IMPLICIT COMPLEX*16(A-D,F-H),REAL*8(E,P-Z)
C
      IF(NN.EQ.1)GO TO 20
C
      ETA2=ETA*ETA
      FA=DCMPLX(1.D0,ETA)
      M=0.25*ETA+4.D1
C
C          POLYNOMES DE TCHEBICHEV JUSQU'AU RANG M
C
      K=M+2
      X=(ETA+ETA)/RO
      XX=X+X-1.D0
      T0=1.D0
      T1=XX
      XX=XX+XX
      DO 6 J=2,K
      TJ=XX*T1-T0
      T0=T1
      T1=TJ
    6 CONTINUE
      TM=T1
      TL=T0
C
C          INITIALISATION
C
      AM=(0.D0,0.D0)
      AL=(1.D-40,1.D-40)
      BN=(0.D0,0.D0)
      BM=(1.D-40,1.D-40)
      BL=(0.D0,0.D0)
      BK=DCMPLX(4.D0*DFLOAT(M+1),0.D0)*AL+BM
      F=(0.D0,0.D0)
      G =(0.D0,0.D0)
      GP=(0.D0,0.D0)
C
C          RECURRENCE DESCENDANTE
C
      K=M
   10 R=K
      TK=XX*TL-TM
      TM=TL
      TL=TK
      HK=DCMPLX(TK,0.D0)
      C1=DCMPLX(R*(R+1.D0)-ETA2,ETA*(R+R+1.D0))
      C2=(4.D0,0.D0)*DCMPLX(R+1.D0,0.D0)
      C2=C2*DCMPLX(-R-1.D0,ETA*3.D0)
      C3=FA*DCMPLX(-R-R-4.D0,ETA)
      C4=DCMPLX((7.D0*R+5.D0)/4.D0,0.D0)
      C5=DCMPLX(R+R+2.D0,0.D0)
      C6=DCMPLX((R+3.D0)/4.D0,0.D0)
      AK=(C2*AL+C3*AM-C4*BL-C5*BM-C6*BN)/C1
      J=K/2
      J=K-J-J
      IF(J)1,2,1
    1 F=F-AK
      GO TO 3
    2 F=F+AK
    3 Z=ABS(AK)
      G=G+HK*AK
      GP=GP+HK*BK
C
C          F=A0/2-A1+A2-A3+A4-A5+...
C
C          CONGRUENCE MODULO 10**60
C
      IF(Z-1.D60)4,5,5
    5 D=(1.D60,0.D0)
      AK=AK/D
      AL=AL/D
      AM=AM/D
      BK=BK/D
      BL=BL/D
      BM=BM/D
      BN=BN/D
      F=F/D
      G=G/D
      GP=GP/D
    4 IF(K)8,8,9
    9 D=DCMPLX(4.D0*R,0.D0)
      BJ=D*AK+BL
      AM=AL
      AL=AK
      BN=BM
      BM=BL
      BL=BK
      BK=BJ
      K=K-1
      GO TO 10
C
C          NORMALISATION ET CALCUL DE Z(RO)
C
    8 D=(0.5D0,0.D0)*AK
      F=F-D
      G=G-D
      GP=GP-(0.5D0,0.D0)*BK
      D=DCMPLX(0.D0,-ETA*DLOG(2.D0)+SIGMA0)
      AXPO=EXP(D)
      F=F/AXPO
      G=G/F
      GP=GP/F
C
C          CALCUL DE F ET G
C
      D=DCMPLX(0.D0,RO-ETA*DLOG(RO))
      AXPO=EXP(D)
      D=G*AXPO
      GP=AXPO*(DCMPLX(0.D0,1.D0-ETA/RO)*G-DCMPLX(X/RO,0.D0)*GP)
      V=D
      D=(0.D0,-1.D0)*D
      U=D
      VP=GP
      GP=(0.D0,-1.D0)*GP
      UP=GP
      IDIV=0
      RETURN
C
C          SERIE ORIGINE
C
   20 PI=3.141592653589793D0
      XA=0.577215664901533D0
      RO2=RO*RO
      ETAP=ETA+ETA
      PIETA=PI*ETA
      Z=138.15510557964276D0
      IDIV=0
      IF(DABS(PIETA)-Z)21,21,22
   22 INDG=IDINT(PIETA/Z)
      IDIV=60*INDG
      IF(ETA.LT.0) IDIV=0
      PIETA=PIETA-Z*DFLOAT(INDG)
   21 PIETA2=0.5D0*PIETA
      P=DEXP(PIETA2)*DSQRT(DSINH(PIETA)/PIETA)
      CALL JFDELG(1.D0,ETA,PAR,PAI,NB)
      Z1=ETAP*(XA+XA+DLOG(2.D0)-1.D0+PAR)
      U0=0.D0
      U1=RO
      V0=1.D0
      V1=Z1*RO
      U=U0+U1
      V=V0+V1
      UP=1.D0
      VP=Z1
      XN=2.D0
      DO 104N=2,10000
      XN1=XN*(XN-1.D0)
      U2=(ETAP*RO*U1-RO2*U0)/XN1
      U=U+U2
      V2=(ETAP*RO*V1-RO2*V0-ETAP*(XN+XN-1.D0)*U2)/XN1
      V=V+V2
      UP=UP+XN*U2/RO
      VP=VP+XN*V2/RO
      IF(DABS(U2/U).GT.1.D-14)GOTO18
      IF(DABS(V2/V).LE.1.D-14)GOTO19
   18 U0=U1
      U1=U2
      V0=V1
      V1=V2
      XN=XN+1.D0
  104 CONTINUE
   19 PP=V+ETAP*U*DLOG(RO)
      W=U/P
      WP=UP/P
      V=P*PP
      VP=P*(VP+ETAP*(UP*DLOG(RO)+U/RO))
      U=W
      UP=WP
      RETURN
      END

      SUBROUTINE YFASYM(ETA,RAU,FO,FPO,GO,GPO,SIGO,IEXP)
      IMPLICIT REAL*8 (A-H,O-Z)
      IEXP=0
      TRB=0.D0
      RAU2=RAU+RAU
      ETAC=ETA*ETA
      CALL JFLGAM(1.D0,ETA,TRA,SIGO,NTRUC)
   40 N=0
      PS=1.D0
      GS=0.D0
      PT=0.D0
      GT=1.D0-ETA/RAU
      SF=PS
      SG=GS
      SPF=PT
      SPG=GT
   45 DENOM= DFLOAT(N+1)*RAU2
      AN= DFLOAT(N+N+1)*ETA/DENOM
      BN= (ETAC-DFLOAT(N*(N+1)))/DENOM
      PS1=AN*PS-BN*PT
      GS1=AN*GS-BN*GT-PS1/RAU
      PT1=AN*PT+BN*PS
      GT1=AN*GT+BN*GS-PT1/RAU
   42 SF=SF+PS1
      SG=SG+GS1
      SPF=SPF+PT1
      SPG=SPG+GT1
      N=N+1
      IF(N-17)46,48,44
   48 TRA=PS*PS+PT*PT
   44 TRB=PS1*PS1+PT1*PT1
      TEST=TRA-TRB
      IF(TEST)47,47,46
   46 PS=PS1
      GS=GS1
      PT=PT1
      GT=GT1
      TRA=TRB
      GOTO 45
   47 TETAO= RAU-ETA*DLOG (RAU2)+SIGO
      TRA= DSIN(TETAO)
      TRB=DCOS(TETAO)
      GO=SF*TRB-SPF*TRA
      GPO=SG*TRB-SPG*TRA
      FO=SPF*TRB+SF*TRA
      FPO=SPG*TRB+SG*TRA
      RETURN
      END

      SUBROUTINE DFCOUL(ETA,RO,F,FP,G,GP,SIGMA,L,IEXP)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION F(*),FP(*),G(*),GP(*),IEXP(*),SIGMA(*)
      DATA DEPI/6.283185307179586D0/
      ETAC=ETA*ETA
      L1=L+1
      CALL DFCZ0(ETA,RO,F0,FP0,G0,GP0,S,I)
      F(1)=F0
      FP(1)=FP0
      G(1)=G0
      GP(1)=GP0
      IEXP(1)=I
      SIGMA(1)=S
      IF(L)1,1,2
    1 RETURN
    2 LINF=0
      IND=0
      IF((ETA.GT.0).AND.(RO.LT.(ETA+ETA)))GO TO 21
      Z=ETA+DSQRT(ETAC+DFLOAT(L*(L+1)))
      IF(RO.GE.Z)GO TO 20
    7 ROINF=ETA+DSQRT(ETAC+DFLOAT(LINF*(LINF+1)))
      IF(RO-ROINF)3,4,4
    4 IF(LINF-L)5,6,6
    5 LINF=LINF+1
      GO TO 7
    3 IND=1
    6 LIN=LINF+1
   20 XM=1.D0
      IF(IND.EQ.0)LIN=L1
      DO 8 J=2,LIN
      ZIG=(DSQRT(ETAC+XM*XM))/XM
      ZAG=ETA/XM+XM/RO
      F(J)=(ZAG*F(J-1)-FP(J-1))/ZIG
      FP(J)=ZIG*F(J-1)-ZAG*F(J)
      G(J)=(ZAG*G(J-1)-GP(J-1))/ZIG
      GP(J)=ZIG*G(J-1)-ZAG*G(J)
      IEXP(J)=I
      SIG=SIGMA(J-1)+DATAN(ETA/(J-1))
      IPI=SIG/DEPI
      SIG=SIG-IPI*DEPI
      IF(SIG)60,50,70
   60 IF(SIG.LT.(-DEPI/2.D0))SIG=SIG+DEPI
      GO TO 50
   70 IF(SIG.GT.(DEPI/2.D0))SIG=SIG-DEPI
   50 SIGMA(J)=SIG
    8 XM=XM+1.D0
      IF(IND.EQ.0)RETURN
      GO TO 22
   21 LIN=1
   22 FTEST=F(LIN)
      FPTEST=FP(LIN)
      LMAX=LINF+25+IDINT(5.D0*DABS(ETA))
      IF(LMAX-L)9,10,10
    9 LMAX=L
   10 FI=1.D0
      FPI=1.D0
   13 XM=LMAX+1
      ZIG=(DSQRT(ETAC+XM*XM))/XM
      ZAG=ETA/XM+XM/RO
      FL=(ZAG*FI+FPI)/ZIG
      FPL=ZAG*FL-ZIG*FI
      IF(DABS(FL)-1.D15)26,27,27
   26 IF(DABS(FPL)-1.D15)28,27,27
   27 FL=FL*1.D-15
      FPL=FPL*1.D-15
   28 FI=FL
      FPI=FPL
      IF(LMAX-L)11,11,12
   12 LMAX=LMAX-1
      GO TO 13
   11 F(LMAX+1)=FL
      FP(LMAX+1)=FPL
      IF(LMAX-LINF)15,15,14
   14 GO TO 12
   15 FACT=FTEST/F(LIN)
      FACTP=FPTEST/FP(LIN)
      INDICE=I/60
      XM=LINF
      DO 16 J=LIN,L1
      F(J)=F(J)*FACT
      FP(J)=FP(J)*FACTP
   25 IF(J.EQ.1)GO TO 16
      ZIG=(DSQRT(ETAC+XM*XM))/XM
      ZAG=ETA/XM+XM/RO
      G(J)=(ZAG*G(J-1)-GP(J-1))/ZIG
      GP(J)=ZIG*G(J-1)-ZAG*G(J)
      IF(DABS(G(J))-1.D60)17,18,18
   17 IF(DABS(GP(J))-1.D60)19,18,18
   18 G(J)=G(J)/1.D60
      GP(J)=GP(J)/1.D60
      INDICE=INDICE+1
   19 IEXP(J)=INDICE*60
      A=FP(J)*G(J)
      B=-F(J)*GP(J)
      IF(A-1.D0)29,30,30
   29 I1=IDINT(DLOG10(A))
      I2=IDINT(DLOG10(B))
      GO TO 31
   30 I1=IDINT(DLOG10(A))+1
      I2=IDINT(DLOG10(B))+1
   31 F(J)=F(J)*1.D1**(-I2)
      FP(J)=FP(J)*1.D1**(-I1)
      SIG=SIGMA(J-1)+DATAN(ETA/(J-1))
      IPI=SIG/DEPI
      SIG=SIG-IPI*DEPI
      IF(SIG)61,51,71
   61 IF(SIG.LT.(-DEPI/2.D0))SIG=SIG+DEPI
      GO TO 51
   71 IF(SIG.GT.(DEPI/2.D0))SIG=SIG-DEPI
   51 SIGMA(J)=SIG
   16 XM=XM+1.D0
      RETURN
      END

      SUBROUTINE YFIREG(ETA,RO,G0,GP0)
      IMPLICIT REAL*8(A-H,O-Z)
      IF(ETA.LE.0.D0)GOTO250
      IF(ETA.LE.3.D0)GOTO251
      IF(ETA.LE.1.D1)GOTO252
      IF(ETA.LE.18.D0)GOTO253
      IF(ETA.LE.22.D0)GOTO254
      IF(RO.LE.0.3D0+(3.D1-ETA)/8.D1)GOTO200
C   SERIE DE TAYLOR DEPART RAU0
  300 CONTINUE
      RAU0=1.666666666666667D0*DABS(ETA)+7.5D0
      CALL YFASYM(ETA,RAU0,F0,FP0,G0,GP0,SIGMA0,IEXP)
      X=RAU0-RO
      X2=X*X
      X3=X*X2
      UNR=1.D0/RAU0
      ETR0=1.D0-2.D0*ETA*UNR
      U0=G0
      U1=-X*GP0
      U2=-0.5D0*ETR0*X2*U0
      S=U0+U1+U2
      V1=U1/X
      V2=2.D0*U2/X
      T=V1+V2
      XN=3.D0
      DO10N=3,10000
C     N=N
      XN1=XN-1.D0
      XN1=XN*XN1
      U3=X*U2*UNR*(1.D0-2.D0/XN)-ETR0*U1*X2/XN1+X3*U0*UNR/XN1
      S=S+U3
      V3=XN*U3/X
      T=T+V3
   16 IF(DABS(U3/S).GT.1.D-11)GO TO 11
      IF(DABS(V3/T).LE.1.D-11)GO TO 20
   11 U0=U1
      U1=U2
      U2=U3
      XN=XN+1.D0
   10 CONTINUE
   20 G0=S
      GP0=-T
      RETURN
C   SERIE  ORIGINE
  200 CONTINUE
      PI=3.141592653589793D0
      GA=0.577215664901533D0
      ETA2=ETA*ETA
      RO2=RO*RO
      ETAP=ETA+ETA
      PIETA=PI*ETA
      PIETA2=0.5D0*PIETA
      B=DEXP(PIETA2)*DSQRT(DSINH(PIETA)/PIETA)
      CALL JFDELG(1.D0,ETA,PAR,PAI,NB)
      C1=ETAP*(GA+GA+DLOG(2.D0)-1.D0+PAR)
      U0=0.D0
      U1=RO
      V0=1.D0
      V1=C1*RO
      U=U0+U1
      V=V0+V1
      UP=1.D0
      VP=C1
      XN=2.D0
      DO 104N=2,10000
      XN1=XN*(XN-1.D0)
      U2=(ETAP*RO*U1-RO2*U0)/XN1
      U=U+U2
      V2=(ETAP*RO*V1-RO2*V0-ETAP*(XN+XN-1.D0)*U2)/XN1
      V=V+V2
      UP=UP+XN*U2/RO
      VP=VP+XN*V2/RO
   17 IF(DABS(U2/U).GT.1.D-14)GOTO18
      IF(DABS(V2/V).LE.1.D-14)GOTO19
   18 U0=U1
      U1=U2
      V0=V1
      V1=V2
      XN=XN+1.D0
  104 CONTINUE
   19 GP=V+ETAP*U*DLOG(RO)
      G0=B*GP
      GP0=B*(VP+ETAP*(UP*DLOG(RO)+U/RO))
      RETURN
  250 IF(RO.LE.0.5D0*ETA+9.D0)GOTO200
      GOTO 300
  251 IF(RO.LE.2.25D0+7.35D0*(3.D0-ETA))GOTO200
      GOTO 300
  252 IF(RO.LE.1.2D0+1.5D-1*(1.D1-ETA))GOTO200
      GOTO 300
  253 IF(RO.LE.0.6D0+0.75D-1*(18.D0-ETA))GOTO200
      GOTO 300
  254 IF(RO.LE.0.4D0+0.5D-1*(22.D0-ETA))GOTO200
      GOTO 300
      END

      SUBROUTINE YFRICA(ETA,RAU,FO,FPO,GO,GPO,SIGMA0,IDIV)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION Q(5),QP(5)
C
C        COEFFICIENTS RICCATI
C
      DATA G61,G62,G63,G64,G65,G66,G67,G68,G69,G610,
     &G611/ 0.1159057617187498D-01,0.3863525390624998D-01,
     &0.4660034179687498D-01,0.4858398437499998D-01,
     &0.1156514485677080D 01,0.5687475585937496D 01,
     &0.1323888288225445D 02,0.1713083224826384D 02,
     &0.1269003295898436D 02,0.5055236816406248D 01,
     &0.8425394694010415D 00/
      DATA G81,G82,G83,G84,G85,G86,G87,G88,G89,G810,G811,G812,G813,G814,
     &G815/ 0.1851092066083633D-01,0.8638429641723630D-01,
     &0.1564757823944092D 00,0.1430139541625977D 00,
     &0.1924622058868408D 00,0.8500803152720129D 01,
     &0.7265429720878595D 02,0.3057942376817972D 03,
     &0.7699689544836672D 03,0.1254157054424285D 04,
     &0.1361719536066055D 04,0.9831831171035763D 03,
     &0.4547869927883148D 03,0.1222640538215636D 03,
     &0.1455524450256709D 02/
      DATA GP61,GP62,GP63,GP64,GP65,GP66/ 0.2897644042968748D-01,
     &0.2318115234375000D 00,0.8056640625000000D 00,
     &0.1601562499999998D 01,0.3046875000000000D 00,
     &0.5624999999999998D 01/
      DATA GP81,GP82,GP83,GP84,GP85,GP86,GP87,
     &GP88/ 0.6478822231292720D-01,0.6910743713378906D 00,
     &0.3322952270507811D 01,0.9483032226562498D 01,
     &0.1769653320312499D 02,0.3478710937499998D 02,
     &0.5020312499999999D 02,0.7874999999999999D 02/
      DATA Q /0.4959570165D-1,0.8888888889D-2,0.2455199181D-2,
     &0.9108958061D-3,0.2534684115D-3/
      DATA QP /0.1728260369D0,0.3174603174D-3,0.3581214850D-2,
     &0.3117824680D-3,0.9073966427D-3/
      CALL JFLGAM(1.D0,ETA,TRA,SIGMA0,IND)
      RAU2=RAU+RAU
      RAUC=RAU*RAU
      ETAC=ETA*ETA
      ETA2=ETA+ETA
      ETARO=ETA*RAU
      ETARO2=ETARO+ETARO
      PIETA=3.141592653589793*ETA
      IND=0
      JND=0
      IG=0
      IF(ETA)20,20,21
   20 IF(-ETARO-14.0625D0)32,22,22
   22 INDICE=1
 
C             RICCATI 3
      IDIV=0
      GO TO 2
   21 IF(DABS(RAU-ETA2).LE.1.D-9)GO TO 18
      IF(RAU-ETA2)30,18,31
   31 IF(RAU-ETA2-2.D1*(ETA**0.25D0))34,33,33
   33 INDICE=0
C             RICCATI  2
      IDIV=0
      GO TO 2
   32 NN=1
      GO TO 35
   34 NN=0
   35 CALL YFCLEN(ETA,RAU,FO,FPO,GO,GPO,SIGMA0,IDIV,NN)
      RETURN
   30 IF(ETARO-12.D0)32,32,23
   23 TRA=ETA2-6.75D0*(ETA**0.4D0)
      IF(RAU-TRA)6,6,24
   24 IND=1
      JND=1
      RO=RAU
      RAU=TRA
      RAU0=TRA
C             RICCATI  1
 
    6 X=RAU/ETA2
      U= (1.D0-X)/X
      X2=X*X
      RU= DSQRT(U)
      RX= DSQRT(X)
      TRE= 1.D0/(U*RU*ETA2)
      TRB=TRE*TRE
      FI= (DSQRT((1.D0-X)*X)+DASIN(RX)-1.570796326794897D0)*ETA2
      TR1= -0.25D0*DLOG(U)
  602 TR2= -((9.D0*U+6.D0)*U+5.D0)/48.D0
  603 TR3= ((((-3.D0*U-4.D0)*U+6.D0)*U+12.D0)*U+5.D0)/64.D0
  604 TR4=- ((((((U+2.D0)*945.D0*U+1395.D0)*U+12300.D0)*U+25191.D0)
     &*U+19890.D0)*U+5525.D0)/46080.D0
  605 TR5= ((((((((-27.D0*U-72.D0)*U-68.D0)*U+360.D0)*U+2190.D0)
     &*U+4808.D0)*U+5148.D0)*U+2712.D0)*U+565.D0)/2048.D0
  606 TR6=- (((((((((G61*U+G62)*U+G63)*U+G64)*U+G65)*U+G66)*U+G67)
     &*U+G68)*U+G69)*U+G610)*U+G611
  607 TR7= ((((((((((((-81.*U-324.)*U-486.)*U-404.)*U+4509.)*U+52344.)
     &*U+233436.)*U+567864.)*U+838521.)*U+775884.)*U+441450.)
     &*U+141660.)*U+19675.) /6144.
  608 TR8= (((((((((((((G81*U+G82)*U+G83)*U+G84)*U+G85)*U+G86)*U+G87)
     &*U+G88)*U+G89)*U+G810)*U+G811)*U+G812)*U+G813)*U+G814)*U+G815
      PSIP=PSIP+TRA
      XXX=138.1551055796428D0
      FI= FI+TRE*(TR2+TRB*(TR4+TRB*(TR6+TRB*TR8)))
      PSI=-FI
      INDG=IDINT(PSI/XXX)
      IDIV=60*INDG
      TRA= TR1+TRB*(TR3+TRB*(TR5+TRB*TR7))
      FI=FI+TRA
      PSI=PSI+TRA
 
      FIP=RU*ETA2
      TRA=1.D0/(X2*U)
      TR1=0.25D0
      TRE= TRE/(X2*X2*U)
      TRB=TRB/(X2*X2)
  702 TR2= -(8.D0*X-3.D0)/32.D0
  703 TR3= ((24.D0*X-12.D0)*X+3.D0)/64.D0
  704 TR4= (((-1536.D0*X+704.D0)*X-336.D0)*X+63.D0)/2048.D0
  705 TR5= ((((1920.D0*X-576.D0)*X+504.D0)*X-180.D0)*X+27.D0)/1024.D0
  706 TR6 = ((((-GP66*X+GP65)*X-GP64)*X+GP63)*X-GP62)*X+GP61
  707 TR7= - ((((((-40320.D0*X-10560.D0)*X-13248.D0)*X+7560.D0)
     &*X-3132.D0)*X+756.D0)*X-81.D0) / 2048.D0
  708 TR8 =- ((((((GP88*X+GP87)*X+GP86)*X-GP85)*X+GP84)*X-GP83)*X+GP82)
     &*X-GP81
      FIP=FIP +TRE*(TR2+TRB*(TR4+TRB*(TR6+TRB*TR8)))
      TRA= TRA*(TR1+TRB*(TR3+TRB*(TR5+TRB*TR7)))
      FIP=FIP+TRA
      PSIP=-FIP
      IF(INDG.EQ.0)GO TO 8
      PSI=PSI-XXX*DFLOAT(INDG)
      FI =FI +XXX*DFLOAT(INDG)
    8 FO=0.5D0*DEXP(FI)
      GO= DEXP(PSI)
      FPO= FO*FIP/ETA2
      GPO=GO*PSIP/ETA2
      IF(JND.EQ.0)RETURN
      RAU=RO
      GO=FO
      GPO=FPO
   27 X=RAU0-RO
      X2=X*X
      X3=X*X2
      UNR=1.D0/RAU0
      ETR0=1.D0-2.D0*ETA*UNR
      U0=GO
      U1=-X*GPO
      U2=-0.5D0*ETR0*X2*U0
      S=U0+U1+U2
      V1=U1/X
      V2=2.D0*U2/X
      T=V1+V2
      XN=3.D0
 
      DO10N=3,10000
C     N=N
      XN1=XN-1.D0
      XN1=XN*XN1
      U3=X*U2*UNR*(1.D0-2.D0/XN)-ETR0*U1*X2/XN1+X3*U0*UNR/XN1
      S=S+U3
      V3=XN*U3/X
      T=T+V3
   16 IF(DABS(U3/S).GT.1.D-10)GO TO 11
      IF(DABS(V3/T).LE.1.D-10)GO TO 17
   11 U0=U1
      U1=U2
      U2=U3
      XN=XN+1.D0
   10 CONTINUE
   17 IF(IG)25,26,25
   25 GO=S
      GPO=-T
      FO=HO
      FPO=HPO
      RETURN
   26 HO=S
      HPO=-T
   18 ET0=ETA**(0.166666666666667)
      ETAD=ETAC*ETAC
      ET=ETA**(0.6666666666666667)
      ET1=ET*ET
      ET2=ET1*ET1
      ET3=ET2*ET
      ET4=ETAD*ET
      ET5=ET4*ET
      FO=1.D0-Q(1)/ET1-Q(2)/ETAC-Q(3)/ET3-Q(4)/ETAD-Q(5)/ET5
      GO=1.D0+Q(1)/ET1-Q(2)/ETAC+Q(3)/ET3-Q(4)/ETAD+Q(5)/ET5
      FPO=1.D0+QP(1)/ET+QP(2)/ETAC+QP(3)/ET2+QP(4)/ETAD+QP(5)/ET4
      GPO=1.D0-QP(1)/ET+QP(2)/ETAC-QP(3)/ET2+QP(4)/ETAD-QP(5)/ET4
      FO=0.7063326373D0*ET0*FO
      GO=1.223404016D0*ET0*GO
      FPO=0.4086957323D0*FPO/ET0
      GPO=-0.7078817734D0*GPO/ET0
      IDIV=0
      IF(IND.EQ.0)RETURN
      IG=1
      RAU0=ETA2
      GO TO 27
    2 X=ETA2/RAU
      X2=X*X
      U=1.D0-X
      RU= DSQRT(U)
      U3=U*U*U
      TRD= 1.D0/(U3*ETA2*ETA2)
      TRC=X2*TRD
      TRE=1.D0/(U*RU*ETA2)
      FI= -0.25D0*DLOG(U)
      TRB=TRD/64.D0
      TR3= (((3.D0*U-4.D0)*U-6.D0)*U+12.D0)*U-5.D0
  501 TR5= ((((((((-27.D0*U+72.D0)*U-68.D0)*U-360.D0)*U+2190.D0)
     &*U-4808.D0)*U+5148.D0)*U-2712.D0)*U+565.D0)/32.D0
  502 TR7= ((((((((((((81.D0*U-324.D0)*U+486.D0)*U-404.D0)*U-4509.D0)
     &*U+52344.D0)*U-233436.D0)*U+567864.D0)*U-838521.D0)*U+775884.D0)
     &*U-441450.D0)*U+141660.D0)*U-19675.D0)/96.D0
      FI= FI+TRB*(TR3+TRD*(TR5+TRD*TR7))
 
      FIP=0.25D0/U
      TRB=3.D0*TRC/(64.D0*U)
      TR3= (X-4.D0)*X+8.D0
  511 TR5= ((((9.D0*X-60.D0)*X+168.D0)*X-192.D0)*X+640.D0)/16.D0
  512 TR7= ((((((-27.D0*X+252.D0)*X-1044.D0)*X+2520.D0)*X-4416.D0)
     &*X-3520.D0)*X-13440.D0)/32.D0
      FIP =FIP+TRB*(TR3+TRC*(TR5+TRC*TR7))
      TRA= DABS((RU-1.D0)/(RU+1.D0))
      PSI= (0.5D0*DLOG(TRA)+RU/X)*ETA2+0.785398163397448D0
      TR2= -((9.D0*U-6.D0)*U+5.D0)/48.D0
  521 TR4= ((((((U-2.D0)*945.D0*U+1395.D0)*U-12300.D0)*U+25191.D0)
     &*U-19890.D0)*U+5525.D0)/46080.D0
  522 TR6 = (((((((((-G61*U+G62)*U-G63)*U+G64)*U-G65)*U+G66)*U-G67)
     &*U+G68)*U-G69)*U+G610)*U-G611
  523 TR8= (((((((((((((G81*U-G82)*U+G83)*U-G84)*U+G85)*U-G86)*U+G87)
     &*U-G88)*U+G89)*U-G810)*U+G811)*U-G812)*U+G813)*U-G814)*U+G815
      PSI= PSI+TRE*(TR2+TRD*(TR4+TRD*(TR6+TRD*TR8)))
      PSIP = -RU*ETA2/X2
      TRB=TRE*X/U
      TR2= (3.D0*X-8.D0)/32.D0
  531 TR4= - (((63.D0*X-336.D0)*X+704.D0)*X-1536.D0)/2048.D0
  532 TR6 = ((((GP61*X-GP62)*X+GP63)*X-GP64)*X+GP65)*X-GP66
  533 TR8 = ((((((-GP81*X+GP82)*X-GP83)*X+GP84)*X-GP85)*X+GP86)*X+GP87)
     &*X+GP88
      PSIP =PSIP+ TRB*(TR2+TRC*(TR4+TRC*(TR6+TRC*TR8)))
      TRA= DEXP(FI)
      FO= TRA*DSIN(PSI)
      GO= TRA*DCOS(PSI)
      IF(INDICE)535,536,535
  535 TRA=FO
      FO=GO
      GO=-TRA
  536 TRA=-ETA2/RAUC
      FPO=(FIP*FO+PSIP*GO)*TRA
      GPO=(FIP*GO-PSIP*FO)*TRA
      RETURN
      END

      SUBROUTINE DFCZ0(ETA,RO,F0,FP0,G0,GP0,SIGMA0,IEXP)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A1(110),A2(110),B1(110),B2(110)
      IF(RO)2,2,1
    2 WRITE (6,1000)
 1000 FORMAT (21H RO NEGATIF OU NUL **)
      RETURN
    1 IF(ETA-30.D0)3,3,4
    4 IF(DABS(ETA)-5.D2)28,28,29
   28 CALL YFRICA(ETA,RO,F0,FP0,G0,GP0,SIGMA0,IEXP)
      RETURN
   29 WRITE (6,1001)
 1001 FORMAT (42H VALEUR ABSOLUE DE ETA SUPE-&EU-E A 500 **)
      RETURN
    3 IF(ETA+8.D0)4,5,5
    5 IF(ETA)6,7,6
    7 F0=DSIN(RO)
      G0=DCOS(RO)
      FP0=G0
      GP0=-F0
      IEXP=0
      SIGMA0=0.D0
      RETURN
    6 BORNE=1.666666666666667D0*DABS(ETA)+7.5D0
      IF(RO-BORNE)8,9,9
    9 CALL YFASYM(ETA,RO,F0,FP0,G0,GP0,SIGMA0,IEXP)
      RETURN
    8 IF(ETA-10.D0)10,11,11
   10 IF(ETA)12,12,13
   13 IF(RO-2.D0)14,12,12
   11 IF(ETA-(5.D0*RO+6.D1)/7.D0)12,12,14
   12 CALL YFASYM(ETA,BORNE,F0,FP0,G0,GP0,SIGMA0,IEXP)
      H=BORNE
      DH=F0/H
      IF(ETA)20,21,21
   20 N=-0.5D0*ETA+5.D0
      GO TO 22
   21 N=ETA/5.D0+5.D0
   22 N=5*(N+1)
      Z=4.D0/H
      Y=1.D0-(ETA+ETA)*Z
      A1(N+2)=1.D-55
      A1(N+3)=0.D0
      A1(N+4)=1.D-64
      B1(N+3)=1.D-50
      B1(N+4)=1.D-68
      A2(N+2)=0.D0
      A2(N+3)=1.D-74
      A2(N+4)=1.D-53
      B2(N+3)=0.D0
      B2(N+4)=1.D-66
      M=N+2
      DI=N
      DO 23 II=2,M
      I=M-II+2
      B1(I)=B1(I+2)+Z*(DI+1.D0)*A1(I+1)
      S=A1(I+2)+Y*(A1(I+1)-A1(I))
      Q=(DI+2.D0)*B1(I)+(DI-1.D0)*B1(I+1)
      A1(I-1)=S-Z*Q
      B2(I)=B2(I+2)+Z*(DI+1.D0)*A2(I+1)
      S=A2(I+2)+Y*(A2(I+1)-A2(I))
      Q=(DI+2.D0)*B2(I)+(DI-1.D0)*B2(I+1)
      A2(I-1)=S-Z*Q
      IF(I.GE.N)GO TO 23
      D=-(B2(I+2)+B2(I))/(B1(I+2)+B1(I))
      DO 24 J=I,M
      A2(J)=A2(J)+D*A1(J)
      B2(J)=B2(J)+D*B1(J)
   24 CONTINUE
      A2(I-1)=A2(I-1)+D*A1(I-1)
   23 DI=DI-1.D0
      Q=A1(3)-A1(1)
      C=A2(3)-A2(1)
      C=Q/C
      X1=0.5D0*(A1(2)-C*A2(2))
      DO 25 I=3,M
      X1=X1+A1(I)-C*A2(I)
   25 CONTINUE
      X1=DH/X1
      X2=-C*X1
      DO 26 I=2,M
      B1(I)=X1*B1(I)+X2*B2(I)
      A1(I)=X1*A1(I)+X2*A2(I)
   26 CONTINUE
      A1(1)=X1*A1(1)+X2*A2(1)
      B1(1)=0.D0
      X=RO/H
      Y=2.D0*X-1.D0
      T1=1.D0
      T2=Y
      RESULT=0.5D0*A1(2)+Y*A1(3)
      DERIVE=0.5D0*B1(2)+Y*B1(3)
      DO 27 I=2,N
      TI=2.D0*Y*T2-T1
      T1=T2
      T2=TI
      RESULT=RESULT+TI*A1(I+2)
      DERIVE=DERIVE+TI*B1(I+2)
   27 CONTINUE
      F0=RESULT*RO
      FP0=DERIVE*RO+RESULT
      GO TO 30
C   SERIE ORIGINE REGULIERE
   14 PI=3.141592653589793D0
      CALL JFLGAM(1.D0,ETA,TRA,SIGMA0,NTRUC)
      IEXP=0
      RO2=RO*RO
      ETAP=ETA+ETA
      PIETA=PI*ETA
      PIETA2=0.5D0*PIETA
      B=DEXP(PIETA2)*DSQRT(DSINH(PIETA)/PIETA)
      U0=0.D0
      U1=RO
      U=U0+U1
      UP=1.D0
      XN=2.D0
      DO 15 N=2,10000
      XN1=XN*(XN-1.D0)
      U2=(ETAP*RO*U1-RO2*U0)/XN1
      U=U+U2
      UP=UP+XN*U2/RO
   17 IF(DABS(U2/U).LT.1.D-10)GO TO 19
   18 U0=U1
      U1=U2
      XN=XN+1.D0
   15 CONTINUE
   19 F0=U/B
      FP0=UP/B
   30 CALL YFIREG(ETA,RO,G0,GP0)
      RETURN
      END

      SUBROUTINE JFDELG (XD,YD,PAR,PAI,NBCHIF)
      DOUBLE PRECISION XD,YD,PAR,PAI,TEST,C,PI
      DOUBLE PRECISION X,Y,U,V,TRA,TRA1,COSI,SINI
      DOUBLE PRECISION COS2I,SIN2I,ZMOD,DEPI
      DOUBLE PRECISION TRB,XX
      DOUBLE PRECISION RAC2,PIS4
      DOUBLE PRECISION SUPINT
      DIMENSION TEST(7),C(6)
      DATA TEST/2.9152D+7,2.2958D+3,1.4124D+2,3.9522D+1,19.6611D0,
     &12.791D0,-10.D0/
      DATA RAC2/0.3465735902799726/,PIS4/0.785398163397448/
      DATA C/8.333333333333333D-2,-8.33333333333333D-3,
     &3.968253968253968D-3,-4.166666666666667D-3,7.575757575757576D-3,
     &-2.109279609279609D-2/
      DATA PI/3.141592653589793/
      DATA DEPI/6.283185307179586/
      DATA SUPINT/2147483647.D0/
      X=DABS(XD)
      XX=X
      NBCHIF=15
      IF(YD)1,2,1
    1 Y=DABS(YD)
      KR=1
      I=DMOD(10.99D0-X,SUPINT)
C     TRANSLATION
      IF(I)3,3,4
    4 TRA=I
      X=X+TRA
C     LOGARITHME(X+IY) (X,Y POSITIFS)
    3 IF(X-Y)5,6,7
    5 TRA1=X/Y
      TRB=1.D0+TRA1*TRA1
      TRA=Y*DSQRT(TRB)
      SINI=1./(TRB*Y)
      COSI=SINI*TRA1
      TRA1=Y/X
      GO TO 11
    6 U=RAC2+DLOG(X)
      V=PIS4
      SINI=0.5D0/X
      COSI=SINI
      GO TO 10
    7 TRA1=Y/X
      TRB=1.D0+TRA1*TRA1
      TRA=X*DSQRT(TRB)
      COSI=1./(TRB*X)
      SINI=COSI*TRA1
   11 U=DLOG(TRA)
      V=DATAN(TRA1)
C     DEVELOPPEMENT ASYMPTOTIQUE ( X SUPERIEUR A 10 )
   10 PAR=U-0.5*COSI
      PAI=V+0.5*SINI
      ZMOD=X+Y
      IF(ZMOD-TEST(1))13,13,14
   13 SIN2I=(SINI*COSI)+(SINI*COSI)
      COS2I=(COSI+SINI)*(COSI-SINI)
      SINI=SIN2I
      COSI=COS2I
      K=1
      GO TO 15
   16 TRA=COSI*COS2I-SINI*SIN2I
      SINI=SINI*COS2I+COSI*SIN2I
      COSI=TRA
   15 PAR=PAR-C(K)*COSI
      PAI=PAI+C(K)*SINI
      K=K+1
      IF(ZMOD-TEST(K))16,16,14
C     TRANSLATION INVERSE
   17 I=I-1
      X=I
      X=XX+X
   56 IF(X-Y)55,55,57
   55 TRA1=X/Y
      TRB=X*TRA1+Y
      SINI=1.D0/TRB
      COSI=TRA1/TRB
      GO TO 19
   57 TRA1=Y/X
      TRB=X+Y*TRA1
      COSI=1.D0/TRB
      SINI=TRA1/TRB
   19 PAR=PAR-COSI
      PAI=PAI+SINI
   14 IF(I)18,18,17
 
C     CONTROLE DU QUADRANT
   18 IF(XD)20,61,21
   61 TRA=PI*Y
      IF(TRA-1.D-2)300,300,301
  300 TRB= TRA*(2.D0+TRA*(-2.D0+TRA*(1.333333333333333D0+TRA*(
     &-0.6666666666666666D0+TRA*(0.2666666666666666D0+TRA*(
     &-0.08888888888888888D0+TRA*0.02539682539682540D0))))))
      TRB=(2.D0-TRB)/TRB
      GO TO 302
  301 TRB= DEXP(-TRA-TRA)
      TRB=(1.D0+TRB)/(1.D0-TRB)
  302 PAI=0.5D0*(1.D0/Y+PI*TRB)
   21 IF(YD)28,100,100
C     X+IY CHANGE EN -X-IY
   20 TRA=DEXP(-DEPI*Y)
      TRB=TRA*TRA
      COS2I=DEPI*DMOD(X,1.D0)
      SIN2I=-2.D0*TRA*DCOS(COS2I)+1.D0+TRB
      PAR=PAR+COSI+DEPI*TRA*DSIN(COS2I)/SIN2I
      PAI=PAI-SINI+PI*(TRB-1.D0)/SIN2I
      IF(YD)100,100,28
   28 PAI=-PAI
   35 WRITE (6,1001)
 1001 FORMAT (30H ARGUMENT DE JFDELG TROP GRAND)
      GO TO 50
   34 Y=DMOD(X,SUPINT)
      IF(Y) 400,36,400
  400 IF(Y-0.99D0) 33,33,405
  405 TRA= IDINT(Y+0.1D0)
      IF(DABS(Y-TRA)-5.D-15)36,36,33
   31 IF(X-4503599627370496.D0)34,35,35
 
C     ARGUMENT DANS -PI,PI
  100 TRA=DABS(PAI/DEPI)
      IF(TRA-1.D+15)203,204,204
  204 NBCHIF=0
      PAI=0.D0
      GO TO 29
  203 IF(TRA-1.D0)205,205,206
  206 NBCHIF=DLOG10(TRA)
      NBCHIF=14-NBCHIF
      TRA=DMOD(TRA,SUPINT)
      PAI=DMOD(TRA,1.D0)*DSIGN(DEPI,PAI)
  205 IF(DABS(PAI)-PI)29,29,207
  207 PAI=PAI-DSIGN(DEPI,PAI)
      GO TO 29
C        DELGAMMA REEL
    2 PAI=0.D0
      IF(XD)31,32,33
C     CONDITIONS D EXISTENCE
   32 WRITE (6,1000)
 1000 FORMAT (21H JFDELG(0) EST INFINI)
      GO TO 50
   36 WRITE (6,1002)
 1002 FORMAT (28H JFDELG (-ENTIER) EST INFINI)
   50 PAR=1.D+74
      NBCHIF=0
      GO TO 29
C     TRANSLATION
   33 I=DMOD(10.99D0-X,SUPINT)
      IF(I)37,37,38
   38 TRA=I
      X=X+TRA
C     DEVELOPPEMENT ASYMPTOTIQUE
   37 Y=DLOG(X)
      PAR=Y-0.5D0/X
      IF(X-TEST(1))39,39,43
   39 COS2I=1.D0/(X*X)
      COSI=COS2I
      K=1
      GO TO 41
   42 COSI=COSI*COS2I
   41 PAR=PAR-C(K)*COSI
      K=K+1
      IF(X-TEST(K))42,42,40
C     TRANSLATION INVERSE
   44 I=I-1
      X=I
      X=XX+X
      PAR=PAR-1.D0/X
   40 IF(I)43,43,44
C     X NEGATIF
   43 IF(XD)45,29,29
   45 PAR=PAR+1.D0/X
      Y=PI*DMOD(X,2.D0)
      PAR=PAR+PI*DCOS(Y)/DSIN(Y)
      ENTRY JFDEV1
   29 RETURN
      END




