subroutine Obs(er1, er2, ner)

c***********************************************************************
c     Obs calculates the observables, including phase shifts, cross section
c     and scattering length etc.
c***********************************************************************
c     LOCAL variables
      integer i, nch, clmax, cl
      double precision hbarc, pi, zero, two, three
      double precision x1, x2, stmp, etmp, rho, psdrop, psold(6)
      double precision SL(6), Phase(6), Eta(6)
      double precision dm, p11, p13, p31, p33
c     SHARED variables
      integer Nifty, npts, ner, nfit, first
      double precision er1, er2, P, W, mpion, Mass, tlab
      double precision Efit(20), exp_xs, th_xs
      double complex enew, eold, Knew(2), Knew2(2), Dk2de(2), Wgf(128)
      double complex Vsuper(128,128,6), Tonshl(2,2,6)

      common /control/ Nifty(15)
      common /mass/ mpion, Mass(2)
      common /grid/ P(64), W(64), npts
      common /data/ exp_xs(20, 6), th_xs(20, 6), nfit

      parameter (hbarc=1.9732858d2, pi=3.1415926d0)
      parameter (zero=0.0d0, two=2.0d0, three=3.0d0)
      data first/1/
c***************************************************
         
      nch = Nifty(2)
      clmax = Nifty(6)
      eold = dcmplx(er1, zero)

      if (Nifty(4) .eq. 2) then
c	 calc phase shift from energy er1 to er2
c        pi-N threshold 1078.9 MeV
         if (er1 .lt. 1078.0) then
            write(*,*)'starting energy (er1) too small, try threshold'
         endif
	 
         do cl = 1, 6
	   psold(cl) = 0.
	 enddo

	 open(11, file='o.on-T', status='unknown')
	 open(2,  file='o.ps-s', status='unknown')
         open(12, file='o.ps-p', status='unknown')
	 open(3,  file='o.sc-len', status='unknown')

	 write(11,100) 
100 	 format('# Partial wave amplitude (onshl Tmatrix)'/,
     &   '# 1st col is energy, followed by Re/Im tmatrix for eigenchannels
     &  S11,S31,P11,P31,P13,P33')

         etmp = er1 - (er2-er1)/(ner-1)
         do 10010 i = 1, ner 
	    etmp = etmp + (er2-er1)/(ner-1)
            enew = dcmplx(etmp, zero)

            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           Devide by relativistic density of state (real), PiN->PiN
	    do cl = 1, clmax
	      rho = Knew(1)*Dk2de(1)
	      Tonshl(1,1,cl) = -Tonshl(1,1,cl) * rho
	    enddo

c           print out for ploting, S11, S31, P11, P31, P13, P33 Re/Im
c           with phase shift anlysis: Arndt,Ford and Roper, PRD32(85)1085

  	    write(11,101) etmp, (Tonshl(1,1,cl),cl=1,clmax) 
101         format(f6.1,12e11.3)

	    do cl = 1, clmax
c             Compute phase shift and inelastic coef.
	      x1 = 2.d0*real(Tonshl(1,1,cl))
              x2 = 1.d0-2.d0*imag(Tonshl(1,1,cl))
	      Phase(cl) = 0.5d0 * atan2(x1,x2) * 180.d0/pi
	      Eta(cl) = sqrt(x1*x1 + x2*x2) 
	      psdrop = Phase(cl) - psold(cl)
	      if (psdrop .lt. -20.) then
		 Phase(cl)=Phase(cl)+180.d0
	      endif
	      psold(cl) = Phase(cl)

c	      Compute scattering length
	      stmp = tan(Phase(cl)*pi/180.0)
	      if (cl .gt. 2) then
	        if (stmp.lt.0.0d0) then
		  stmp = -(-stmp)**(0.333333d0)
	          else
		    stmp = stmp**(0.333333d0)
		    endif
	      endif
	      SL(cl) = stmp*hbarc/real(Knew(1))
	    enddo
	    write(2,201) etmp, (Phase(cl),Eta(cl),cl=1,2)
            write(12,201)etmp, (Phase(cl),Eta(cl),cl=3,6),
     &			dsin(Phase(6)*pi/180.0)**2
	    tlab = 0.5*(etmp+mpion+Mass(1))*(etmp-mpion-Mass(1))/Mass(1)
	    write(3,301)tlab, real(Knew(1)), (SL(cl),cl=1,clmax)
201	    format(f6.1,6(e12.3,f6.2),e10.3)
301	    format(f6.1,e10.3,6e12.4)
     
c            if (Nifty(9) .eq. 1) then
c            calc total cross sections sigma(nch)?
c               call Xsctn(enew, Knew, Dk2de, Tonshl, tlm, sigma)
c	    
c            else if (Nifty(9) .eq. 2) then
c            calc diff cross sections and polarizations
c               write(8, *) 'X and pola'
c            endif
	
10010    continue
c	 end loop for energy axis
	close(11)
	close(2)
	close(3)
	
	else

c	  calc phase shifts for selected energy points for fitting
c         p11,p13,p31,p33 are phase shifts in units of degree
c         dm is a useless var only for convenience (data format)
c         parameters er1, er2, and ner are irrelevant in this case

	  if (first .eq. 1) then
	    nfit = 5

	    open(4,file='datafit',status='unknown')

	    do i = 1, nfit
	     read(4,*) dm, p11,dm, p13,dm, p31,dm, p33,dm, Efit(i)
	     exp_xs(i,3) = dsin(p11*pi/180.)**2
             exp_xs(i,6) = dsin(p33*pi/180.)**2
	     write(8,900) Efit(i), p11,exp_xs(i,3), p33,exp_xs(i,6)
	    enddo
900	    format('Fit XS:', 5e13.5)
	    close(4)
	    first = 0
	  endif

	  do 10020 i = 1, nfit
            etmp = Efit(i)
            enew = dcmplx(etmp, zero)

            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           Devide by relativistic density of state (real), PiN->PiN
            do cl = 3, clmax, 3
              rho = Knew(1)*Dk2de(1)
              Tonshl(1,1,cl) = -Tonshl(1,1,cl) * rho

c             Compute phase shift and inelastic coef.
              x1 = 2.d0*real(Tonshl(1,1,cl))
              x2 = 1.d0-2.d0*imag(Tonshl(1,1,cl))
              th_xs(i,cl) = dsin(0.5d0 * atan2(x1,x2))**2
	    enddo
10020	  continue  
	endif

c       print out the potential parameters
c        call PrtPot()

        return
        end