
      SUBROUTINE DSTRIP(IQ,DR,K,F1,F2,FR,QNUM,OPT,KM,SL,C)
c
c     Calculates two nucleon transfer form factor (Bayman and Kallio)
c
c     IQ    Input quantum numbers for form factor
c     DR    Step size
c     K     Number of steps
c     F1    First  orbital
c     F2    Second orbital
c     FR    Output form factor
c     QNUM  Quantum numbers for orbitals
c     OPT
c
      parameter (maxg=10, maxr=12)
      DIMENSION F1(400),F2(400),FR(800),QNUM(4,2),IQ(4),TVCC(-9:9)
     1         ,D1(-9:9),D2(-9:9),C(32),AG(maxg),WG(maxg),BG(maxg)
     2         ,AR(maxr),WR(maxr)
C
      DATA KR,KX,krel/0,0,0/
      data ag, wg, bg, ar, wr/maxg*0.,maxg*0.,maxg*0.,maxr*0.,maxr*0./
C
c     c1 = R1 scale
c     c2 = R2 scale
c     c3 = r1 scale
c     c4 = r2 scale
c     c5 = r0, integration scale
c     c6 = Pauli flag
c     c7 = order of integration
c     c8 = number of integration points
c
      IPAULI=C(6)
      R0    =C(5)
      L    = IQ(1)
      IS   = IQ(2)
      JJ   = IQ(3)
      lrel = iq(4)
      N1=QNUM(1,1)
      N2=QNUM(1,2)
      L1=QNUM(2,1)
      L2=QNUM(2,2)
      J1=QNUM(3,1)
      J2=QNUM(3,2)
      IS1=QNUM(4,1)
      IS2=QNUM(4,2)
      c3=c(3)
      c4=c(4)
      dr1=dr*c(1)
      dr2=dr*c(2)
      ITEMP=N1+N1+N2+N2
      NX=(((ITEMP+L1+L2+L+2)/4+2)+3+4*lrel)/2+2*IPAULI
      if(c(8).eq.0.0) then
        NR=(((ITEMP+L1+L2-L+3)/4+8)+1+2)/2+2*IPAULI
        NR=MIN0(NR,maxr)
        IF(NR.NE.KR) then
          CUT=0.0
          IERR=0
          ALFA=C(7)
          CALL GAUSSR(NR,KS,ALFA,AR,WR,IERR,CUT)
          NR=KS
          kr=nr
        endif
      else
        nr = c(8)
        do 45 i=1,nr
        ar(i)=c(i+8)
        wr(i)=c(i+nr+8)
        if(wr(i).ne.0.0) kr=i
   45   continue
      nr=kr
      endif
      NX=MIN0(NX,maxg)
      IF(NX .NE. KX .or. lrel .ne. krel) then
        CALL LEGAUS(2*NX,AG,WG)
        KX=NX
        krel=lrel
        DL1=sqrt(YXFCT(0,2*lrel))*YXFCT(lrel,2*lrel)/2.0**lrel
     1     *sqrt(float(2*lrel+1))*phasef(lrel)
        do 47 i=1,nx
        BG(i) = sqrt(1.0-ag(i)**2)
        DM1=0.0
        DO 46 LI=1,lrel
        DL1=DL1*BG(i)
   46 CONTINUE
        wg(i) = wg(i)*DL1
   47   continue
      endif
      WRITE(6,9000)NR,NX
 9000 FORMAT(15H0 NO. R STEPS =,I3,18H     NO. X STEPS =,I3)
c
      LL=L+L
      LL1=L1+L1
      LL2=L2+L2
      FACT1=YXFCT(L1,LL1)/2.0**L1
      FACT2=YXFCT(L2,LL2)/2.0**L2
      FNORM=OPT*SQRT(FLOAT((LL1+1)*(LL2+1))/FLOAT(LL+1))
     1     /vcc(LL,2*lrel,LL,0,2*lrel)
      IF(IPAULI. eq. 0) then
        TEMP=2.0
        DO 50 I=1,4
        IF(QNUM(I,1).NE.QNUM(I,2)) GO TO 55
   50   CONTINUE
        TEMP=4.0
   55   CONTINUE
        FNORM=FNORM/SQRT(TEMP)
      endif
c
      FNORM=FNORM*SQRT(FLOAT((LL+1)*(IS+1)*(J1+1)*(J2+1)))
     1     *WINEJ(LL1,IS1,J1,LL2,IS2,J2,LL,IS,JJ)
      if(l1 .le. l2) then
        mmin= -min0(l1,l2-lrel)
        mmax= l1
      else
        mmin=-min0(l2,l1-lrel)
        mmax= l2+lrel
      endif
      DO 80 M=mmin,mmax
      TVCC( M)=VCC(LL1,LL2,LL, 2*M,-2*M+2*lrel)
     1   *SQRT(yxfct(l1+abs( m     ),l1-abs( m     ))
     2        *yxfct(l2+abs(-m+lrel),l2-abs(-m+lrel)))
c	write(20,'(a,2i4,2f10.5)')' m1, m2, Tvcc(m1,m2)'
c     1	     ,m,-m+lrel,Tvcc( m)
   80 CONTINUE
C
C     RS=r
C
C     R1=| C1*R+C3*r |
C     R2=| C2*R+C4*r |
C
C
C     CENTER OF MASS R LOOP
C
      R=0.0
      S=0.0
      DO 500 M =1,K
      R=R+DR1
      S=S+DR2
      RSQ=R**2
      SSQ=S**2
      SUMR=0.0
C
C     RELATIVE R LOOP
C
c     RS = r/2
      DO 400 MR=1,KR
      RS=AR(MR)*R0
      SS=RS*C3
      RS=RS*C4
      RSSQ=RS**2+RSQ
      RPROD=2.0*R*RS
      SSSQ=SS**2+SSQ
      SPROD=2.0*S*SS
C
C     RELATIVE R ANGLE LOOP
C
      SUMX=0.0
      DO 300 MX=1,KX
      X =AG(MX)
      y =BG(mx)
      IX=0
  110 CONTINUE
      R1=SQRT(RSSQ+RPROD*X)
      R2=SQRT(SSSQ+SPROD*X)
      FK1=R1/DR
      K1=FK1
      K1=MAX0(K1,2)
      FK1=FK1-FLOAT(K1)
      IF(K1.GT.K) GO TO 300
      FK2=R2/DR
      K2=FK2
      K2=MAX0(K2,2)
      FK2=FK2-FLOAT(K2)
      IF(K2.GT.K) GO TO 300
      COS1=(R+RS*X)/R1
      COS2=(S+SS*X)/R2
      SIN1=abs(rs*y/r1)
      SIN2=abs(ss*y/r2)
c
  120 CONTINUE
      DM1=0.0
      DL1=FACT1
      DO 140 LI=1,L1
      DL1=DL1*SIN1
  140 CONTINUE
c     time reversal phase
      D1( L1)=DL1*phasef(L1)
      D1(-L1)=DL1
      DM2=0.0
      DL2=FACT2
      DO 150 LI=1,L2
      DL2=DL2*SIN2
  150 CONTINUE
c     time reversal phase * exp(im*pi)
      D2( L2)=DL2
      D2(-L2)=DL2*phasef(L2)
      FJ1=1.0
      FL1=LL1
      FM1=LL1
      DO 170 m1=L1-1,0,-1
      DK1=(FM1*COS1*DL1/SIN1-DM1)/(FJ1*FL1)
      FJ1=FJ1+1.0
      FL1=FL1-1.0
      FM1=FM1-2.0
      DM1=DL1
      DL1=DK1
c     time reversal phase
      D1( m1)=DL1*phasef(m1)
      D1(-m1)=DL1
  170 CONTINUE
c	if(m .eq. 1 .and. mr .eq. 1) then
c	write(20,'(a,10(i4,f10.5))')' m1, D1 ='
c     1	     ,(m1,d1(m1),m1=l1,-l1,-1)
c	endif
      FJ2=1.0
      FL2=LL2
      FM2=LL2
      DO 180 m2=L2-1,0,-1
      DK2=(FM2*COS2*DL2/SIN2-DM2)/(FJ2*FL2)
      FJ2=FJ2+1.0
      FL2=FL2-1.0
      FM2=FM2-2.0
      DM2=DL2
      DL2=DK2
c     time reversal phase * exp(im*pi)
      D2( m2)=DL2
      D2(-m2)=DL2*phasef(m2)
  180 CONTINUE
c	if(m .eq. 1 .and. mr .eq. 1) then
c	write(20,'(a,10(i4,f10.5))')' m2, D2 ='
c     1	     ,(m2,d2(m2),m2=l2,-l2,-1)
c	endif
      PROD=0.0
      DO 185 LI=mmin,mmax
      PROD=PROD+D1( LI)*D2(-LI+lrel)*TVCC(LI)
c	if(m .eq. 1 .and. mr .eq. 1) then
c	write(20,'(a,2i4,4f10.5)')' m1, m2, D1, D2, Tvcc, Prod ='
c     1	     ,li,-li+lrel,d1(li),d2(-li+lrel),Tvcc(li),prod
c	endif
  185 CONTINUE
  280 CONTINUE
C     INTERPOLATE - 4 POINT FORMULA
      FT1=-FK1*(FK1-1.)*(FK1-2.)*F1(K1-1)/6.
     1  +   (FK1**2-1.)*(FK1-2.)*F1(K1  )/2.
     2  -  FK1*(FK1+1.)*(FK1-2.)*F1(K1+1)/2.
     3  +        FK1*(FK1**2-1.)*F1(K1+2)/6.
      FT2=-FK2*(FK2-1.)*(FK2-2.)*F2(K2-1)/6.
     1  +   (FK2**2-1.)*(FK2-2.)*F2(K2  )/2.
     2  -  FK2*(FK2+1.)*(FK2-2.)*F2(K2+1)/2.
     3  +        FK2*(FK2**2-1.)*F2(K2+2)/6.
      SUMX=SUMX+WG(MX)*PROD*FT1*FT2
c	if(m.eq.34) write(20,'(2i4,1p10e12.4)')mr,mx,x,r1,r2,r1**2+r2**2
c     1    ,ft1,ft2,ft1*ft2,prod,wg(mx),sumx
      IF(IX.eq.0) then
        IX=1
        ITEMP=K1
        K1=K2
        K2=ITEMP
        ATEMP=FK1
        FK1=FK2
        FK2=ATEMP
        IF(L1.EQ.L2) GO TO 280
        ATEMP=COS1
        COS1=COS2
        COS2=ATEMP
        ATEMP=SIN1
        SIN1=SIN2
        SIN2=ATEMP
        GO TO 120
      endif
  300 CONTINUE
      SUMR=SUMR+WR(MR)*SUMX
	if(m.eq.34) write(20,'(2i4,1p10e12.4)')m,mr,rs,sumx,wr(mr)
     1	   ,sumr,fnorm,fnorm*sumr
  400 CONTINUE
      SUMR=SUMR*FNORM
      FR(2*M-1)=FR(2*M-1)+SUMR
      IF(M.EQ.KM) SL=SUMR
  500 CONTINUE
 1000 CONTINUE
      RETURN
      END

