
      SUBROUTINE CATHEN(LOPT,IOPT,JOPT,K ,KMAX
     1                 ,UB,VB,SI,DRT,IREL,ICON19,IPRT,IERR,CONST)
c
      logical itnt
c
      DIMENSION UB(2*kmax,2),VB(4*kmax),SI(2*kmax),C(10),QNUM(4,2),G(4)
     1,IQN(4),BKIN(32),JI(2),JJ(2),FJT(2),CON(32),CONST(3),ik9(2)
c
      EQUIVALENCE (EN,BKIN(1)),(FM,BKIN(2)),(Z,BKIN(3)),(FMA,BKIN(4))
     1,(ZA,BKIN(5)),(RY,BKIN(6)),(FS,BKIN(7)),(ECM,BKIN(8))
     2,(FMU,BKIN(9)),(RC,BKIN(10)),(FK2,BKIN(11)),(FK,BKIN(12))
     3,(ETA,BKIN(13)),(AFACT,BKIN(14)),(PFACT,BKIN(15)),(DRF,BKIN(16))
     4,(FN,BKIN(17)),(PNLOC,BKIN(18)),(amass,bkin(19)),(vfact,BKIN(20))
     5,(sfact,bkin(21)),(fk0,bkin(22))
     6,(BIF,BKIN(28)),(G,BKIN(29))
     7,(C(1),CNTROL),(C(2),QCODE),(C(3),FMUV),(C(4),OPT)
c
	common/ch3files/input, ioutp
      DATA SQR4PI,SQRTEN,PI/3.54490780,3.162277660,3.1415926/
      DATA ECMT/0.0/
C
      iz=0
      itnt=.false.
      zero=0.0 
      IQFLG=0
      KMFLG=0
      km=0
      amass = 0.0
      IQN(1)=LOPT
      IQN(2)=IOPT
      IQN(3)=JOPT
      KT=KMAX
  100 CONTINUE
      READ (input,9001)C
      CNTRL=ABS(CNTROL)
      IF(CNTRL.EQ.0.0) GO TO 3010
      if(icon19  .eq.0) then
      WRITE(ioutp,9051)QCODE,FMUV,OPT
      endif
      IF(QCODE)2020,3000,1100
 1100 CONTINUE
C
C     LOOP FOR TWO ORBITALS
C
      iprint = 0
      if(iprt.eq.4 .or. icon19.ne.0) then iprint = 1
      ML=1
      DO 2015 I=1,2
      FJT(1)=C(5)
      FJT(2)=C(I+5)
      IF(I.EQ.1 .or. CNTRL.EQ.2.0) then
      itemp=0
      CALL KINEX(BKIN,UB(1,2),ECMT,DRT,itemp,KMAX,IREL,iprint,CONST)
      DO 2003 M=1,KMAX
      MK=M+M-1
      VB(MK  )=0.0
      VB(MK+1)=0.0
 2003 CONTINUE
C
C
      CALL POTENT(BKIN,UB(1,2),VB,K,KMAX,2,iz,iz,zero,IERR,iprint)
      FJ2=G(3)
      if(kmflg.eq.0) then
      km=bif
      KMFLG=1
      endif
      KT=MIN0(KT,K)
 2005 CONTINUE
C      IF(FJT(1).EQ.0.0) FJT(2)=FJ2
C
C     PRINT OUT SINGLE PARTICLE INFORMATION
C
      if(icon19  .eq.0) then
      WRITE(ioutp,9501)
      if(icode.le.3 .or. icode.eq.6) WRITE(ioutp,9100)FJ2,FJT
      IF(IPRT.ne.2) then
        R=DRF
        K2=K+K
        DO 2006 M=1,K2,20
        MK4=MIN0(M+19,K2)
        WRITE(ioutp,9052)R,(UB(N,2),N=M,MK4,2)
        R=R+DRF*10.0
 2006   CONTINUE
      endif
      endif
      endif
c
      ik9(i)=k
      DO 2008 M=1,kmax
      MK=M+M-1
      UB(ML,1)=UB(MK,2)
      ML=ML+1
 2008 CONTINUE
      DO 2009 M=1,4
      QNUM(M,I)=G(M)
 2009 CONTINUE
      JI(I)=FJT(1)
      JJ(I)=FJT(2)
 2015 CONTINUE
 2020 ICODE=ABS(QCODE)
      LTR1=IQN(1)+1
C     ENTER TIME REVERSAL PHASE
      LVR=QNUM(2,1)+QNUM(2,2)
      OPT=OPT*PHASEF((LVR+IQN(1))/2)
      if    (icode.le.10.or. icode.gt.20) then
        iflag=1
      elseif(icode.gt.10.and.icode.le.20) then
        iflag=2
        icode=icode-10
      endif
C
C     ICODE= 1  YUKAWA      POTENTIAL
C     ICODE= 2  COULOMB     POTENTIAL
C     ICODE= 3  OPEP TENSOR POTENTIAL
C     ICODE= 4  TNT off-shell FR correction
C     ICODE= 5  TWO PARTICLE TRANSFER
C     ICODE= 6  ZERO RANGE KNOCK OUT
C     Icode= 9  (e,e'D) Knock-out. laguerre integration, with    exp(ik_ir/2)
C     Icode=10  (e,e'D) Knock-out. hermite  integration, with    exp(ik_ir/2)
C
c     Icode=11-20  gives imaginary form factor
c     Icode=29  (e,e'D) Knock-out. laguerre integration, without exp(ik_ir/2)
c     Icode=30  (e,e'D) Knock-out. hermite  integration, without exp(ik_ir/2)
c
      GO TO (2100,2100,2300,2400,2500,2100,2990,2990,2900,2900),ICODE
      go to (2990,2990,2990,2990,2990,2990,2990,2990,2900,2900),icode-20
      write(ioutp,'(a,i3)')' Invalid ICODE parameter in CATHEN =',icode
      GO TO 2990
c
 2100 LVR=0
      TOPT=1.0
      GO TO 2350
 2300 ICODE=1
      LVR=2
      TOPT=-SQRTEN
 2350 MNNL=IABS(LTR1-LVR-1)+1
      MXXL=LTR1+LVR
      LTR2=LTR1+LTR1-2
      LVR2=LVR+LVR
      LLX=QNUM(2,1)+QNUM(2,1)
      LPX=QNUM(2,2)+QNUM(2,2)
      JLX=QNUM(3,1)
      JPX=QNUM(3,2)
      IS1=QNUM(4,1)
      IS2=QNUM(4,2)
      OPT=OPT*PHASEF(LPX/2)
      DRX=DRT
      DO 2390 LAM =MNNL,MXXL,2
      LAM2 =LAM +LAM -2
      VOPT=TOPT*OPT*SQR4PI*SQRT(FLOAT(LAM2 +1)*FLOAT(IQN(2)+1))
     1 *VCC(LAM2 ,LVR2,LTR2 ,0,0)
     2 *RACAH(IQN(2),IQN(2),LAM2 ,LTR2 ,LVR2,IQN(3))
      RME=0.0
      IF(IQN(2).EQ.0) RME=1.0
      IF(IQN(2).EQ.2) RME=SQRT(FLOAT(IS1*(IS1+2)))
      RME=RME*SQRT(FLOAT((IS1+1)*(IQN(2)+1)))
     1*SQRT(FLOAT((LLX+1)*(LAM2+1)))*VCC(LLX,LAM2,LPX,0,0)
     2 *SQRT(FLOAT(JLX+1)*FLOAT(JPX+1)*FLOAT(IQN(3)+1))
     3 *WINEJ(LPX,LLX,LAM2 ,IS2,IS1,IQN(2),JPX,JLX,IQN(3))
     4 *PHASEF((JI(1)+IQN(3)-JLX-JJ(2))/2)
     5 *RACAH(JPX,JJ(2),JLX,JJ(1),JI(1),IQN(3))
     6 *SQRT(FLOAT(JJ(1)+1)*FLOAT(JJ(2)+1))
      VOPT=VOPT*RME
     1 *(DRF/DRX)**3
      SL=0.0
      IF(VOPT.EQ.0.0) GO TO 2380
      IF(ICODE.EQ.6) GO TO 2365
C
      CALL SLATR (KT,KMAX,DRX,VB(  1),     LTR1,FMUV,ICODE)
      CALL SLATR (KT,KMAX,DRX,VB(2*kmax+1),LAM ,FMUV,ICODE)
C
C
      CALL RADIN
     1(KT,KMAX,DRX,FMUV,VB,UB(1,1),UB(kmax+1,1),SL,VOPT,SI(IFLAG),KM,SK)
C
      GO TO 2380
 2365 CONTINUE
      SCALE=VOPT/SQR4PI**2
      R=0.0
      MK=IFLAG
      DO 2370 M=1,KT
      R=R+DRF
      TEMP=UB(M    ,1)*UB(M+kmax,1)*SCALE
      SI(MK  )=SI(MK  )+TEMP
      SL      =SL      +TEMP*R**2
      IF(M.EQ.KM) SK=TEMP
      MK=MK+2
 2370 CONTINUE
      SL=SL*DRF
 2380 CONTINUE
      I=LAM -1
      if(icon19  .eq.0) then
      WRITE(ioutp,9102)SL,IQN(1),I,RME
      WRITE(ioutp,9002)KM,SK
      endif
 2390 CONTINUE
      GO TO 2990
c
c     TNT off-shell FR correction
c
 2400 continue
      itnt = .true.
      call fnrgcor(iqn,kmax,kt,opt,ub,si(iflag),qnum,km,sl)
      if(icon19  .eq.0) then
      write(ioutp,9002)km,sl
      endif
      iqflg=1
      go to 2990
C
C     HERE FOR TWO NUCLEON TRANSFER
C
 2500 CONTINUE
c     c1 = R1 scale
c     c2 = R2 scale
c     c3 = r1 scale
c     c4 = r2 scale
c     c5 = r0, integration scale length
c     c6 = Pauli flag
c     c7 = order  of integration
c     c8 = number of integration points
c
      itnt = .true.
      IF(FMUV.EQ.0.0) FMUV=1.7
      CON(1)= 1.0
      CON(2)= 1.0
      CON(3)= 0.5
      CON(4)=-0.5
      CON(5)= 2.0*fmuv
      con(6)= 0.0
      con(7)= 0.5
      con(8)= 0.0
      IQN(4) = 0
      OPT=OPT*16.0/PI
      CALL DSTRIP
     1(IQN,DRF,KT,UB(1,1),UB(kmax+1,1),SI(IFLAG),QNUM,OPT,KM,SL,CON)
      if(icon19  .eq.0) then
      WRITE(ioutp,9002)KM,SL
      endif
      IQFLG=1
      GO TO 2990
c
c     (e,e'D) Form factor
c
 2900 continue
      itnt = .true.
      ktemp=12
        CON(1)= 1.0
        CON(2)= 1.0
        CON(3)= 0.5
        CON(4)=-0.5
        CON(6)= 0.0
        con(8)= ktemp
c     icode =  9 and 29  Gauss - Laguerre
c           = 10 and 30  Gauss - Hermite
      if(icode.eq.9.or.icode.eq.29) then
        con(7)= 2.0
      else
        con(7)= 0.5
      endif
c
      DO 2903 M=1,KMAX
      MK=M+M-1
      VB(MK  )=0.0
      VB(MK+1)=0.0
 2903 CONTINUE
C
      itemp=-1
      if(c(10).ne.0.0) then
        drx=c(10)
      else
        drx=0.20
      endif
      drtemp=drf
      CALL KINEX (BKIN,VB(1),ECMT,drx,itemp,KMAX,IREL,icon19,CONST)
      kz=0
      CALL POTENT(BKIN,VB(1),VB(2*kmax-1),kz,KMAX
     1      ,iz,iz,iz,zero,IERR,icon19)
c
c     Set lrel for deuteron overlap
c
      IQN(4) = JI(1)
      lrel   = JI(1)     
      if(lrel .ne. 0) IQN(3) = JJ(1)
      dr0=drf
      drf=drtemp
      if(icon19  .eq.0) then
      WRITE(ioutp,9501)
      IF(IPRT.NE.2) then
      R=DR0
      K2=2*kz
      DO 2910 M=1,K2,20
      MK4=MIN0(M+19,K2)
      WRITE(ioutp,9052)R,(VB(N),N=M,MK4,2)
      R=R+DR0*10.0
 2910 CONTINUE
      endif
      endif
      kk=kz
      k0=max0(kk    /2,kk    -30)
      k1=max0(ik9(1)/2,ik9(1)-30)
      k2=max0(ik9(2)/2,ik9(2)-30)
      k4=max0(1,min0(k0/2,5))
      k5=k0+k4
      k6=k1+k4
      k7=k2+k4
c
      if(fmuv.eq.0.0) then
c     find fall-off length
      if    (icode.eq. 9.or.icode.eq.29) then
        flog0=    (float(k0+1)*vb(2*k0+1  )/(float(k5+ 1)*vb(2*k5+1  )))
        if(flog0.gt.0.0) then
          flam0=log(flog0)
     1  /(dr0*float(k4))
        else
          flam0=0.0
        endif
        flam1= log(float(k1+1)*ub(k1+     1,1)
     1  /(float(k6+ 1)*ub(k6+     1,1)))/(drf*float(k4))
c
        flam2= log(float(k2+1)*ub(k2+kmax+1,1)
     1  /(float(k7+ 1)*ub(k7+kmax+1,1)))/(drf*float(k4))
c
        r0  =1.0/(flam0+flam1/2.0+flam2/2.0)
      elseif(icode.eq.10.or.icode.eq.30) then
        flam0= log(vb(2*k0+  1)/vb(2*k5+        1))
     1       /(float(k5+ 1)**2-float(k0+ 1)**2)/dr0**2
c
        flam1= log(ub(  k1+  1,1)/ub(k6+        1,1))
     1       /(float(k6+ 1)**2-float(k1+ 1)**2)/drf**2
c
        flam2= log(ub(  k2+kmax+1,1)/ub(k7+kmax+1,1))
     1       /(float(k7+ 1)**2-float(k2+ 1)**2)/drf**2
c
	r0  =1.0/sqrt(flam0+flam1/4.0+flam2/4.0)
      else
        r0 =1.0
      endif
      else
        r0  =1.0/fmuv
      endif
      con(5)= r0
      jerr=0
      call gaussr(ktemp,ks,con(7),con(9),con(9+ktemp),jerr,cut)
      if(jerr.ne.0)
     1   write(ioutp,'(a,i3)')' Gaussr does not converge for root ',jerr
c	write(20,'(1h ,1p10e12.4)')con
c	write(20,*)flam0,flam1,flam2
c	write(20,*)' r0 = ',r0,'   ks =', ks,'   fk0 =',fk0
c     interpolate to get weights
c
      sum1=0.0
      sum2=0.0
      do 2920 i=1,ks
      t1=con(i+8)*r0/dr0
      j= max(int(t1),2)
      p= t1-float(j)
      if(j.lt.kk-2) then
c     4 point formula
        ff = (-        p*(1.-p)*(2.-p)  *vb(2*j-3)*float(j-1)/6.0
     1        + (1.+p)  *(1.-p)*(2.-p)  *vb(2*j-1)*float(j  )/2.0
     2        + (1.+p)*p        *(2.-p) *vb(2*j+1)*float(j+1)/2.0
     3        - (1.+p)*p*(1.-p)         *vb(2*j+3)*float(j+2)/6.0)/t1
        con(i+ks+8)=con(i+ktemp+8)*ff
        if    (icode.eq. 9.or.icode.eq.29) then
          con(i+ks+8)=con(i+ks+8)*exp(con(i+8)   )*r0**3
        elseif(icode.eq.10.or.icode.eq.30) then
          con(i+ks+8)=con(i+ks+8)*exp(con(i+8)**2)*r0**3
        else
          write(*,'(a,i3)')' Icode error in Cathen, icode =',icode
        endif
          arg = fk0*r0*con(i+8)/2.
          t1=cos(arg)/arg
          t2=sin(arg)/arg
          do 2915 j=1,lrel
          t0 = t1
          t1 = t2
          t2 = float(2*j-1)*t1/arg-t0
 2915     continue
          t3 = t2*con(i+ks+8)*float(2*lrel+1)
        if(icode.eq. 9.or.icode.eq.10) then
          con(i+ks+8)=t3
        endif
      else
        con(i+ks+8)=0.0 
      endif
      sum1=sum1+con(i+8+ks)
      sum2=sum2+t3*ff
c	write(20,'(2i3,1p10e12.4)')i,int(con(8)),con(i+8),arg
c     1,      t2,ff,con(i+8+ks)
c
c     N_strength   \int \phi_d**2 exp(iq*r/2) dr_vec 
c     D_strength   \int \phi_d exp(iq*r/2) d r_vec   icode =  9 and 10
c     D_strength   \int \phi_d  d r_vec              icode = 29 and 30
c
 2920 continue
      if(icon19  .eq.0) then
      write(ioutp,'(3(a,1pe12.4))')' N_strength =',sum2
     1     ,'    D_strength =',sum1*sqr4pi,  '    Scale =',r0
      endif
      CALL DSTRIP
     1(IQN,DRF,KT,UB(1,1),UB(kmax+1,1),SI(IFLAG),QNUM,OPT,KM,SL,CON)
      if(icon19  .eq.0) then
      WRITE(ioutp,9002)KM,SL
      endif
      IQFLG=1
      go to 2990
 2990 CONTINUE
 3000 CONTINUE
      IF(CNTROL.GT.0.0) GO TO 100
      IF(IQFLG.EQ.1) KT=KT-5
      K=KT
 3010 CONTINUE
      DRT=BKIN(16)
      if(itnt) then
        sum1=0.0
        r=0.0
        do 3020 m=1,kt
        r=r+drf
        mk = 2*m -2 +iflag
        sum1=sum1+(si(mk)*r)**2
 3020   continue
        if(icon19  .eq.0) then
        write(ioutp,'(a,1pe12.4)')' FF_norm    =',sum1*drf
        endif
      endif
      RETURN
 9001 FORMAT(10F8.4)
 9002 FORMAT(15H FORM FACTOR,M=,I3,1PE18.6)
 9051 FORMAT(18H0PARAMETERS-2 nucl,9H   QCODE=,F9.4,9H   RANGE=,F9.4
     1,9h   VZERO=,F9.4)
 9052 FORMAT(1H ,F6.2,1p10E12.4)
 9100 FORMAT(27H0SINGLE PARTICLE FUNCTIONS ,8HCOUPLING ,6H 2*J1=,F3.0
     1 ,9H    2*J2=,F3.0,9H    2*JI=,F3.0)
 9102 FORMAT(19H VOLUME INTEGRAL = ,F9.4,8H    LTR=,I3,8H    LAM=,I3
     1,8H    RME=,F8.4)
 9501 FORMAT(1H ,18HFORM FACTOR DATA  )
 9998 FORMAT(1H1)
      END

      subroutine fnrgcor(iq,kmax,kt,opt,f,fr,qnum,km,sl)
c
      dimension f(kmax,2),fr(4*kmax),qnum(4,2),iq(3)
c
      data sqr4pi /3.5449078/
c
      l =iq(1)
      is=iq(2)
      jj=iq(3)
      l1 =qnum(2,1)
      l2 =qnum(2,2)
      is1=qnum(4,1)
      is2=qnum(4,2)
      j1 =qnum(3,1)
      j2 =qnum(3,2)
      ll =l +l
      ll1=l1+l1
      ll2=l2+l2
c
      fnorm=opt
     1     *sqrt(float((jj+1)*(is+1)*(j1+1)*(j2+1)*(ll1+1)*(ll2+1)))
     2     *winej(ll2,is2,j2,ll1,is1,j2,ll,is,jj)
     3     *vcc(ll2,ll1,ll,0,0)/sqr4pi
     4     *phasef((2+j1+j2-jj)/2)
c
      do 100 m=1,kt
      mk=2*m-1
      t1 = fnorm*f(m,1)*f(m,2)
      fr(mk)=fr(mk) - t1
      if(m.eq.km) sl = t1
  100 continue
      return
      end

