c
      subroutine VBag(npts, P, ze, Zk, Vsuper)
c
c***********************************************************************
c     VBag calculates super potential matrix. 
c        Vsuper_ba(k', k, l), where l stands for eigen-channel (up to 6) 
c     Cplv, Cplct and Cplcs are SU(N) coupling constants in isospin base.
c***********************************************************************
c     LOCAL variables
      integer first, nch, npts, npts1, clmin, clmax, cl
      integer ich, jch, ipt, jpt, ichpt, jchpt
      double precision zero, two
c     SHARED variables
      integer Nifty, Fidx
      double precision P(64), Wm(64), Wmst(64), rcbm 
      double precision f0, mgb, alpha, c1, c2, Mbare, mpion, Mass
      double precision Uborn(64), Uct(64,64,0:1), Ucs(64,64,0:1)
      double precision Urop(64)
      double complex Wm0(2), Wm0st(2), Vsuper(128,128,6)
      double complex ze, Zk(2), zUborn(2),zUrop(2)
      double complex zUct1(2,64,0:1),zUcs1(2,64,0:1)
      double complex zUct2(64,2,0:1),zUcs2(64,2,0:1)
      double complex zUct0(2,2,0:1),zUcs0(2,2,0:1)

      parameter (zero=0.0d0, two=2.0d0)
      common /control/ Nifty(15)
      common /bag/ rcbm,f0,mgb,alpha,c1,c2,Mbare(3),Fidx(3)
      common /mass/ mpion, Mass(2)
      common /cpl/Uborn
      data first /1/
c********************

      npts1 = npts+1

c     do following only at first time or after model parameters change 
      if (first .eq. 1) then
         nch = Nifty(2)
         clmax = Nifty(6)
	 clmin = 1
	 if (Nifty(4) .eq. 4) clmin = clmax

c        calc meson energy for each grid point using relat kinematics
         do ipt = 1, npts
              Wm(ipt)=sqrt(mpion**2+P(ipt)*P(ipt))
              Wmst(ipt)=sqrt(two*Wm(ipt))
         enddo
	 first = 0
	endif
c     end if (first)

      call FfInt(npts,P,Zk,Uborn,Urop,Uct,Ucs,
     & zUborn,zUrop,zUct1,zUct2,zUct0,zUcs1,zUcs2,zUcs0)

c     Compute meson onshell point energy and its square root
c     Zk(nch) is on-shell channel momentum, complex, used in LS equ.
      do ich = 1, nch
            Wm0(ich) = sqrt(Zk(ich)**2+mpion**2)
            Wm0st(ich) = sqrt(two*Wm0(ich))
      enddo

c********************************************************************
c       The pot matrix Vsuper is of nch*nch block. Each block (i,j)  
c     is a npts1*npts1 complex matrix corresp'g to the trans'n pot 
c     from j- to i-chan'l. The 1st npts points are gaussian grids 
c     and the last one is on the shell needed for calc of T-matrix. 
c     For calc pot only, the on-shell momentum is set to P(npts)
c       By sym, only lower half (inc diagnal) of Vsuper are computed.
c     Upper half are set to lower half after all diag h'n included.
c
c     initialize the whole potential matrix
      do cl = 1, clmax
         do jchpt = 1, nch*npts1
            do ichpt = 1, nch*npts1
               Vsuper(ichpt, jchpt, cl) = zero
            enddo
         enddo
      enddo

c     The order of 6 partial waves: l=S11,S31,P11,P31,P13,P33
c     start to call Diag [s, u, t]

      call DiagSU(npts,P,ze,Zk,Wm,Wmst,Wm0,Wm0st
     &             ,Mbare,f0,Uborn,zUborn,Urop,zUrop,Vsuper)
       if (Nifty(10) .ge. 3) then
           call DiagT(npts,P,ze,Zk,Wm,Wmst,Wm0,Wm0st,f0,Uct,Ucs
     &             ,zUct1,zUct2,zUct0,zUcs1,zUcs2,zUcs0,Vsuper)
       endif


c     now, we get another half of potential matrix
      do cl = clmin, clmax
c         do j = 1, nch*npts1
c            do i = 1, j-1
c	donot use symmetry of potntial
c               Vsuper(i, j, cl) = Vsuper(j, i, cl)
c           enddo
c         enddo

        do 103 jch = 1, nch
        do 104 ich = 1, nch
            do jpt = 1, npts
               jchpt = jpt + (jch-1)*npts1
               do ipt = 1, jpt
                  ichpt = ipt + (ich-1)*npts1
                  Vsuper(jchpt,ichpt,cl) = Vsuper(ichpt,jchpt,cl)
               enddo
               ichpt = ich*npts1
               Vsuper(jchpt,ichpt,cl) = Vsuper(ichpt,jchpt,cl)
            enddo

            jchpt = jch*npts1
            do ipt = 1, npts
               ichpt = ipt + (ich-1)*npts1
               Vsuper(jchpt,ichpt,cl) = Vsuper(ichpt,jchpt,cl)
            enddo
            ichpt = ich*npts1
            jchpt = jch*npts1
            Vsuper(jchpt,ichpt,cl) = Vsuper(ichpt,jchpt,cl)
  104    continue
  103    continue
	enddo
      return
      end