subroutine Amp(er1, er2, ner, ei1, ei2, nei)
c
c***********************************************************************
c     Amp calculates the scattering amplitudes (3D Tmatrix pole)
c     for given complex energy range er1, er2, ei1, ei2
c     if nif11=0, calc PiNN cpling const and ff, possibly for PiDD later
c	       1, 3D diagonal tmatrix vs E
c	       2, 3D inverse det(1-vg) vs E (incorrect)
c***********************************************************************
c     
      integer cl, i, istep, rstep,  c
      integer first, iflag, ii, ileft,iright
      double precision zero, hbarc,  de, de1, de2, tmp, rmd 
      double precision fpi1, ff2ren1(64),ff2ren2(64)
      double precision c0, s1, s2, Iprop, ff2bar1(64), ff2bar2(64)
c     
      integer Nifty, nch, ner, nei, npts
      double precision er1, er2, ei1, ei2, rsize, isize, er, ei, P, W
      double precision mpion, Mass, fact, pi
      double precision Evec(2), Detvec(2)
      double complex enew, eold, Knew(2), Knew2(2), Dk2de(2), Wgf(128)
      double complex Vsuper(128,128,6), Tonshl(2,2,6)
      double complex Thalf(64), Tdiag(64), Te(64,240), Ve(64,240)
      double complex epole, ebare
      common /half/ Thalf, Tdiag
      common /grid/ P(64), W(64), npts
      common /control/ Nifty(15)
      common /mass/ mpion, Mass(2)

      parameter (zero=0.0d0, hbarc=1.9732858d2, pi=3.1415926d0)
      data first/1/
c******************
c
      nch = Nifty(2)
      cl = Nifty(6)
      if (cl.eq.3) then
	rmd = 1.
	else if (cl.eq.6) then
	rmd = (2.*1.414/5.0)**2
	endif
      call PrtPot

      if (Nifty(11) .eq. 0) then
        open(9, file='PiNNee', status='unknown')
	write(9, 2000) er1, ei1
2000    format('# Energy dependnece of ff**2(upto a const):Resonance 
     &E =', 2e12.5/
     &	'#    E',4x,'~Tonshl(1,1,3)',2x,'~Tdiag(1)',2x,'~Vbare(1,1,3)')
c       plot diagonal t matrix, start 10 MeV below resonance and bare pole
c	use (er1,ei1) as physical pole, use (er2,ei2) as bare mass
c	epole and ebare is 0.001MeV away from exact
        epole = dcmplx(er1, ei1) + 0.001 
	ebare = dcmplx(er2, ei2) + 0.001
        ner = 40

c    renormalized cpling const	
	write(*,*) 'Calculate Renomalized Cpling and FF'
        er = er1 - 10.0
	eold = dcmplx(er, ei1)
        do    rstep = 1, ner
          er = er + 0.5
          enew = dcmplx(er, ei1)
          call OnShl(enew, eold, Knew, Knew2, Dk2de)
          call WtgFcn(npts, P, W, enew, Knew, Knew2, Dk2de, Wgf)
          call VBag(npts, P, enew, Knew, Vsuper)
          call TMtrx(npts, Dk2de, Wgf, Vsuper, Tonshl)
c          write(*,1005) er,Tdiag(1),Vsuper(1,1,3)
	  if (first.eq.1 .and. Nifty(10).eq.1) then
c           calc extra factor c0 in cpling const due to mass renormalization
            s1 = Iprop(er1-1.0, P, W, npts)
            s2 = Iprop(er1+1.0, P, W, npts)
            c0 = 1. + (s2 - s1)/2.0
	    write(*,*) '*** c0,s1,s2 ***', c0, s1, s2
	    first = 0
	  else 
	    c0 = 1.0
	  endif
c
	  fact = 3.*rmd*(P(1)/mpion)**2/(2.*sqrt(mpion**2+P(1)**2))
          fpi1 = Tdiag(1)*(enew - epole)/fact
          write(9,2002) er, fpi1/c0, fpi1
	  do i = 1, npts
            Te(i,rstep) = Tdiag(i)
          enddo
        enddo
	
c       subtract the ff(q)**2 at 7 MeV away from pole (both side)
	ii = 4
	ileft = 2
	de = er1 - 10.0 + ileft*0.5 - er1
	de1 = 1./de - 1./(de+0.5*ii)
	iright = 34
	de = er1 - 10.0 + iright*0.5 - er1
	de2 = 1./de - 1./(de+0.5*ii)

	do i = 1, npts
	  fact = 3.*rmd*(P(i)/mpion)**2/(2.*sqrt(mpion**2+P(i)**2))
	  ff2ren1(i) = (Te(i,ileft)-Te(i,ileft+ii))/de1/fact
 	  ff2ren2(i) = (Te(i,iright)-Te(i,iright+ii))/de2/fact
	enddo

c   bare cpling const
	write(*,*) 'Calculate Bare Cpling and FF'
	write(9,2003)
        er = er2 - 10.0
        eold = dcmplx(er, 0.0d0)
        do    rstep = 1, ner
          er = er + 0.500001
          enew = dcmplx(er, 0.0d0)
          call OnShl(enew, eold, Knew, Knew2, Dk2de)
          call WtgFcn(npts, P, W, enew, Knew, Knew2, Dk2de, Wgf)
          call VBag(npts, P, enew, Knew, Vsuper)
          fact = 3.*rmd*(P(1)/mpion)**2/(2.*sqrt(mpion**2+P(1)**2))
          fpi1 = Vsuper(1,1,cl)*(enew - ebare)/fact
	  write(*,1005) er, Vsuper(1,1,cl), fpi1
          write(9,2002) er, fpi1
          do i = 1, npts
            Ve(i,rstep) = Vsuper(i,i,cl)
          enddo
        enddo
	open(10, file='PiNNff', status='unknown')
        de = er2 - 10.0 + ileft*0.5 - er2
        de1 = 1./de - 1./(de+0.5*ii)
        de = er2 - 10.0 + iright*0.5 - er2
        de2 = 1./de - 1./(de+0.5*ii)

        do i = 1, npts
          fact = 3.*rmd*(P(i)/mpion)**2/(2.*sqrt(mpion**2+P(i)**2))
          ff2bar1(i) = (Ve(i,ileft)-Ve(i,ileft+ii))/de1/fact
          ff2bar2(i) = (Ve(i,iright)-Ve(i,iright+ii))/de2/fact
          write(10,2002) P(i),ff2ren1(i),ff2ren2(i),
     &ff2bar1(i),ff2bar2(i),ff2ren1(i)/ff2ren1(1),ff2bar1(i)/ff2bar1(1)
        enddo
2002      format(1x,e12.5,9e11.4)
2003	  format(/)
1005      format('Amp:er ',e12.4,6e11.3)
c********t*********2*********3*********4**********5**********6*********7
      else if (Nifty(11) .ge. 1) then
       	open(9, file='Amp.re', status='unknown')
      	open(10, file='Amp.im', status='unknown')
	c = Nifty(6)

c    	 calculating step size
        if (ner .ne. 1) then
          rsize = (er2-er1)/(ner-1)
      	  else
            rsize = zero
      	    endif
        if (nei .ne. 1) then
          isize = (ei2-ei1)/(nei-1)
          else
            isize = zero
            endif

        eold = dcmplx(er1, ei1)
        do istep = 1, nei
          ei = ei1 + isize*(istep-1)
          do    rstep = 1, ner
            er = er1 + rsize*(rstep-1)
	  if (Nifty(11) .eq. 1) then
            enew = dcmplx(er, ei)
            call OnShl(enew, eold, Knew, Knew2, Dk2de)
            call WtgFcn(npts, P, W, enew, Knew, Knew2, Dk2de, Wgf)
            call VBag(npts, P, enew, Knew, Vsuper)
            call TMtrx(npts, Dk2de, Wgf, Vsuper, Tonshl)

c           print out scattering amplitudes as TRIUMF's convention
c	    show pole using 1:2:4 3D plot
            call PrtKin(enew, Knew)
            write(9,1001) er, ei, -hbarc*dble(Dk2de(1)*Tonshl(1,1,c)),
     &	    -hbarc*dble(Dk2de(1)*Tdiag(1)),dble(Tdiag(1)),dble(Tdiag(2))
	    write(10,1001) er, ei, -hbarc*imag(Dk2de(1)*Tonshl(1,1,c)),
     &	    -hbarc*imag(Dk2de(1)*Tdiag(1)),imag(Tdiag(1)),imag(Tdiag(2))
          else if (Nifty(11) .eq. 2) then
c           Detvec(2) return det(1-vg)
            iflag = 1
            Evec(1) = er
            Evec(2) = ei
            call Det1Mvg(2, Evec, Detvec, iflag)
            tmp = Detvec(1)**2 + Detvec(2)**2
	    write(9,1001) er, ei, Detvec(1)/tmp, 1.0/Detvec(1)
	    write(10,1001) er, ei, -Detvec(2)/tmp, 1.0/Detvec(2)
	  endif
c
            eold = enew
	  enddo
	   write(9,1002)
	   write(10,1002)
        enddo
      endif
1001	format(8e12.4)
1002	format(/)

      return
      end
	
c     propogator integral for subtract the extra factor for cpling
	double precision function Iprop(ener, P, W, npts) 
        integer npts, ipt, Fidx
        double precision Uborn(64), Eg(64), ener
        double precision mpion, Mass, fact, pi, P(64), W(64)
        double precision f0, mgb, alpha, c1, c2, Mbare, rcbm, sum
        
        common /mass/ mpion, Mass(2)
        common /bag/ rcbm,f0,mgb,alpha,c1,c2,Mbare(3),Fidx(3)
	common /cpl/ Uborn
	parameter (pi=3.1415926d0)

c     relativistic kinematics
        do ipt = 1, npts
          Eg(ipt) = sqrt(P(ipt)**2 + mpion**2)
     &                     +sqrt(P(ipt)**2 + Mass(1)**2)
        enddo
c
	sum = 0.0
        do ipt = 1, npts
 	  fact = P(ipt)**4/sqrt(P(ipt)**2 + mpion**2)	  
          sum=sum + fact*Uborn(ipt)**2*W(ipt)/(ener-Eg(ipt))
        enddo
	Iprop = sum * 25.0/(12.*pi*f0**2)/4./pi
c	the extra 1/4pi factor is necessary to cancel the denominator of reson
	write(*,*) ener,Iprop,sum
	return
	end