      module size
!============================================
! Choose input!
!============================================
      integer, parameter :: nprecision = 100
      integer, parameter :: nintegprec = 10
      integer, parameter :: ntotintprec = 1
      real (8), parameter :: emax = 12.d0!maximum energy for gauss integration (smearing)
!============================================
      end module size

      module storage
      use size
      real(8), dimension (6) :: coef
      real(8), dimension (5) :: choice
      character(5), dimension (1) :: strings
      save coef
      save choice
      save strings
      end module storage

!     ------------------------------------------------------------------------
      include 'useful_constants.f'
      include 'models.f'
      include 'intx.f'
!     ------------------------------------------------------------------------

      program prepare_plot
      use size
      use storage
      use hadron_properties
      implicit none
!     input
      integer        :: unitinput
      integer	     :: unitout_dxs,unitout_dxsang,unitout_totxs
      integer	     :: unitout_kll,unitout_kllang
      integer	     :: unitout_rho00,unitout_rho00ang
      integer	     :: unitout_rhop0,unitout_rhop0ang
      integer	     :: unitout_rhopm,unitout_rhopmang
      character (60) :: inputfilename
      character (60) :: outname_dxs,
     c  outname_dxsang,outname_totxs,
     c  outname_kll,outname_kllang,
     c  outname_rho00,outname_rho00ang,
     c  outname_rhop0,outname_rhop0ang,
     c  outname_rhopm,outname_rhopmang
      character (5) :: spin
!     internal variables
      integer :: iloop1,iloop2,iloop3,nzoo,observable
      real (8) :: coef1,coef2,coef3,coef4,coef5,
     c     coef6,coef7,coef8,a12flag,a32flag,smearing,a12in,enermax
      real (8) :: dxs_temp,ang_temp,totxs_temp
      real (8) :: kll_temp,kllang_temp
      complex (8) :: rho00_temp,rho00ang_temp
      complex (8) :: rhop0_temp,rhop0ang_temp
      complex (8) :: rhopm_temp,rhopmang_temp
      real (8) :: xgamr,xmassr,ener,tman,xcos,enerpeak
      real (8), external :: dxs_spin32,dxs_spin52,totxs
     c		,dxs_spin32_smear,dxs_spin52_smear,kll32
     c		,kll32_smear,kll52,kll52_smear
      complex (8) :: rho0032
     c		,rho0032_smear,rho0052,rho0052_smear
      complex (8) :: rhop032
     c		,rhop032_smear,rhop052,rhop052_smear
      complex (8) :: rhopm32
     c		,rhopm32_smear,rhopm52,rhopm52_smear
      real(8) :: x0,x1,res,xtol,res2,dran
      logical :: t
!     code
      unitinput      = 1
      unitout_dxs     = 2
      unitout_totxs     = 3
      unitout_dxsang     = 4
      unitout_kll     = 5
      unitout_kllang     = 6
      unitout_rho00     = 7
      unitout_rho00ang     = 8
      unitout_rhop0     = 9
      unitout_rhop0ang     = 10
      unitout_rhopm     = 11
      unitout_rhopmang     = 12
      inputfilename  = './Input/test.txt'
      outname_dxs = '../gnuplot/Data/dxs.txt'
      outname_dxsang = '../gnuplot/Data/ang.txt'
      outname_totxs = '../gnuplot/Data/totxs.txt'
      outname_kll = '../gnuplot/Data/kll.txt'
      outname_kllang = '../gnuplot/Data/kllang.txt'
      outname_rho00 = '../gnuplot/Data/rho00.txt'
      outname_rho00ang = '../gnuplot/Data/rho00ang.txt'
      outname_rhop0 = '../gnuplot/Data/rhop0.txt'
      outname_rhop0ang = '../gnuplot/Data/rhop0ang.txt'
      outname_rhopm = '../gnuplot/Data/rhopm.txt'
      outname_rhopmang = '../gnuplot/Data/rhopmang.txt'
      open (unitinput, file = inputfilename )
      read (unitinput,*)
     c        coef1,coef2,coef3,coef4,coef5,coef6,coef7,coef8,
     c        spin,smearing,a12in,observable,enermax
      close (unitinput)	
!======================================================================
! Test
!======================================================================
      a12flag=a12in
      a32flag=dsqrt(1.d0-a12in**2)
      coef(1)=coef1
      coef(2)=coef2
      coef(3)=coef3
      coef(4)=coef4
      coef(5)=coef5
      coef(6)=coef6
      xmassr = coef7
      xgamr = coef8
      strings(1) = spin
      choice(1) = a12flag
      choice(2) = a32flag
      choice(3) = observable
      choice(4) = smearing
      choice(5) = enermax
      enerpeak =(xmassr**2-xmpro**2)/(2.d0*xmpro)
!======================================================================
! Gives the DXS for each parameter set, as well as mean, and sigmas, as
! a function of energies close to threshold.
!======================================================================
      if (observable.eq.1) then
         open (unitout_dxs, file = outname_dxs)
         tman = 1.d0
         nzoo = 2
         dran = 3.d0
         if (smearing.eq.0.d0) then
            do iloop1 = 1, nzoo*nprecision
               ener=(enerpeak-1.5d0)+dran
     c		*iloop1/dble(nzoo*nprecision)
               if (spin.eq.'32') then
                  dxs_temp = dxs_spin32(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52') then
                  dxs_temp = dxs_spin52(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_dxs,1002)ener,dxs_temp
            enddo
            do iloop1 = 1, nprecision
               ener=(enerpeak+1.5d0)+(enermax-(enerpeak+1.5d0))
     c		*iloop1/dble(nprecision)
               if (spin.eq.'32') then
                  dxs_temp = dxs_spin32(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52') then
                  dxs_temp = dxs_spin52(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_dxs,1002)ener,dxs_temp
            enddo
	 elseif (smearing.gt.0.d0) then
            do iloop1 = 1, nzoo*nprecision
               ener=(enerpeak-1.5d0)+dran*iloop1/dble(nzoo*nprecision)
               if (spin.eq.'32') then
                  dxs_temp = dxs_spin32_smear(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52') then
                  dxs_temp = dxs_spin52_smear(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_dxs,1002)ener,dxs_temp
            enddo
         endif
         close (unitout_dxs)
         endif

!======================================================================
! Gives the total XS for each parameter set, as well as mean,
! and sigmas, as a function of energies close to threshold.
!======================================================================
      if (observable.eq.2) then
         open (unitout_totxs, file = outname_totxs)
         do iloop2 = 1, 2*nprecision
               ener=(enerpeak-1.5d0)+3.d0*iloop2/dble(2*nprecision)
               totxs_temp = totxs(ener,xgamr,xmassr)
            write (unitout_totxs,1002) ener
     c		,totxs_temp
         enddo
         close (unitout_totxs)
      endif

!======================================================================
! Gives the angular behavior of DXS for each parameter set, as well as
! mean, and sigmas, at the resonance-peak energy.
!======================================================================
      if (observable.eq.3) then
         open (unitout_dxsang, file = outname_dxsang)
         do iloop3 = 1, nprecision-1
            xcos = -1.d0 + 2.d0*iloop3/dble(nprecision)
               ener = (xmassr**2-xmpro**2)/(2.d0*xmpro)
               tman =  (-2.d0*ener**2*xmpro + xmjpsi**2*xmpro + 

     -    	ener*(xmjpsi**2 + xcos*dsqrt(xmpro*(2.d0*ener + xmpro))
     -        	*dsqrt((xmjpsi**4 + 4.d0*ener**2*xmpro**2 - 
     -          4.d0*xmjpsi**2*xmpro*(ener + xmpro))/
     -		(xmpro*(2.d0*ener + xmpro)))))/
     -  	(2.d0*ener + xmpro)
               if (spin.eq.'32'.and.smearing.eq.0.d0) then
                  ang_temp = dxs_spin32(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'32'.and.smearing.gt.0.d0) then
                  ang_temp = dxs_spin32_smear(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52'.and.smearing.eq.0.d0) then
                  ang_temp = dxs_spin52(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52'.and.smearing.gt.0.d0) then
                  ang_temp = dxs_spin52_smear(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_dxsang,1003) xcos,tman
     c			,ener,ang_temp
         enddo
         close (unitout_dxsang)
      endif

!======================================================================
! Gives the KLL for each parameter set, as well as mean, and sigmas, as
! a function of energies close to threshold.
!======================================================================
      if (observable.eq.4) then
         open (unitout_kll, file = outname_kll)
         tman = 1.d0
         nzoo = 2
         dran = 3.d0
         if (smearing.eq.0.d0) then
            do iloop1 = 1, nzoo*nprecision
               ener=(enerpeak-1.5d0)+dran
     c		*iloop1/dble(nzoo*nprecision)
               if (spin.eq.'32') then
                  kll_temp = kll32(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52') then
                  kll_temp = kll52(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_kll,1002)ener,kll_temp
            enddo
            do iloop1 = 1, nprecision
               ener=(enerpeak+1.5d0)+(enermax-(enerpeak+1.5d0))
     c		*iloop1/dble(nprecision)
               if (spin.eq.'32') then
                  kll_temp = kll32(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52') then
                  kll_temp = kll52(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_kll,1002)ener,kll_temp
            enddo
	 elseif (smearing.gt.0.d0) then
            do iloop1 = 1, nzoo*nprecision
               ener=(enerpeak-1.5d0)+dran*iloop1/dble(nzoo*nprecision)
               if (spin.eq.'32') then
                  kll_temp = kll32_smear(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52') then
                  kll_temp = kll52_smear(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_kll,1002)ener,kll_temp
            enddo
         endif
         close (unitout_kll)
         endif

!======================================================================
! Gives the angular behavior of KLL for each parameter set, as well as
! mean, and sigmas, at the resonance-peak energy.
!======================================================================
      if (observable.eq.5) then
         open (unitout_kllang, file = outname_kllang)
         do iloop3 = 1, nprecision-1
            xcos = -1.d0 + 2.d0*iloop3/dble(nprecision)
               ener = (xmassr**2-xmpro**2)/(2.d0*xmpro)
               tman =  (-2.d0*ener**2*xmpro + xmjpsi**2*xmpro + 

     -    	ener*(xmjpsi**2 + xcos*dsqrt(xmpro*(2.d0*ener + xmpro))
     -        	*dsqrt((xmjpsi**4 + 4.d0*ener**2*xmpro**2 - 
     -          4.d0*xmjpsi**2*xmpro*(ener + xmpro))/
     -		(xmpro*(2.d0*ener + xmpro)))))/
     -  	(2.d0*ener + xmpro)
               if (spin.eq.'32'.and.smearing.eq.0.d0) then
                  kllang_temp = kll32(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'32'.and.smearing.gt.0.d0) then
                  kllang_temp = kll32_smear(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52'.and.smearing.eq.0.d0) then
                  kllang_temp = kll52(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52'.and.smearing.gt.0.d0) then
                  kllang_temp = kll52_smear(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_kllang,1003) xcos,tman
     c			,ener,kllang_temp
         enddo
         close (unitout_kllang)
      endif

!======================================================================
! Gives rho00 for each parameter set, as well as mean, and sigmas, as
! a function of energies close to threshold.
!======================================================================
      if (observable.eq.6) then
         open (unitout_rho00, file = outname_rho00)
         tman = 1.d0
         nzoo = 2
         dran = 3.d0
         if (smearing.eq.0.d0) then
            do iloop1 = 1, nzoo*nprecision
               ener=(enerpeak-1.5d0)+dran
     c		*iloop1/dble(nzoo*nprecision)
               if (spin.eq.'32') then
                  rho00_temp = rho0032(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52') then
                  rho00_temp = rho0052(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_rho00,1002)ener,rho00_temp
            enddo
            do iloop1 = 1, nprecision
               ener=(enerpeak+1.5d0)+(enermax-(enerpeak+1.5d0))
     c		*iloop1/dble(nprecision)
               if (spin.eq.'32') then
                  rho00_temp = rho0032(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52') then
                  rho00_temp = rho0052(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_rho00,1002)ener,rho00_temp
            enddo
	 elseif (smearing.gt.0.d0) then
            do iloop1 = 1, nzoo*nprecision
               ener=(enerpeak-1.5d0)+dran*iloop1/dble(nzoo*nprecision)
               if (spin.eq.'32') then
                  rho00_temp = rho0032_smear(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52') then
                  rho00_temp = rho0052_smear(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_rho00,1002)ener,rho00_temp
            enddo
         endif
         close (unitout_rho00)
         endif

!======================================================================
! Gives the angular behavior of rho00 for each parameter set, as well as
! mean, and sigmas, rho00at the resonance-peak energy.
!======================================================================
      if (observable.eq.7) then
         open (unitout_rho00ang, file = outname_rho00ang)
         do iloop3 = 1, nprecision-1
            xcos = -1.d0 + 2.d0*iloop3/dble(nprecision)
               ener = (xmassr**2-xmpro**2)/(2.d0*xmpro)
               tman =  (-2.d0*ener**2*xmpro + xmjpsi**2*xmpro + 
     -    	ener*(xmjpsi**2 + xcos*dsqrt(xmpro*(2.d0*ener + xmpro))
     -        	*dsqrt((xmjpsi**4 + 4.d0*ener**2*xmpro**2 - 
     -          4.d0*xmjpsi**2*xmpro*(ener + xmpro))/
     -		(xmpro*(2.d0*ener + xmpro)))))/
     -  	(2.d0*ener + xmpro)
               if (spin.eq.'32'.and.smearing.eq.0.d0) then
                  rho00ang_temp = rho0032(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'32'.and.smearing.gt.0.d0) then
                  rho00ang_temp = rho0032_smear(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52'.and.smearing.eq.0.d0) then
                  rho00ang_temp = rho0052(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52'.and.smearing.gt.0.d0) then
                  rho00ang_temp = rho0052_smear(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_rho00ang,1003) xcos,tman
     c			,ener,rho00ang_temp
         enddo
         close (unitout_rho00ang)
      endif



!======================================================================
! Gives rhop0 for each parameter set, as well as mean, and sigmas, as
! a function of energies close to threshold.
!======================================================================
      if (observable.eq.8) then
         open (unitout_rhop0, file = outname_rhop0)
         tman = 1.d0
         nzoo = 2
         dran = 3.d0
         if (smearing.eq.0.d0) then
            do iloop1 = 1, nzoo*nprecision
               ener=(enerpeak-1.5d0)+dran
     c		*iloop1/dble(nzoo*nprecision)
               if (spin.eq.'32') then
                  rhop0_temp = rhop032(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52') then
                  rhop0_temp = rhop052(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_rhop0,1002)ener,rhop0_temp
            enddo
            do iloop1 = 1, nprecision
               ener=(enerpeak+1.5d0)+(enermax-(enerpeak+1.5d0))
     c		*iloop1/dble(nprecision)
               if (spin.eq.'32') then
                  rhop0_temp = rhop032(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52') then
                  rhop0_temp = rhop052(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_rhop0,1002)ener,rhop0_temp
            enddo
	 elseif (smearing.gt.0.d0) then
            do iloop1 = 1, nzoo*nprecision
               ener=(enerpeak-1.5d0)+dran*iloop1/dble(nzoo*nprecision)
               if (spin.eq.'32') then
                  rhop0_temp = rhop032_smear(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52') then
                  rhop0_temp = rhop052_smear(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_rhop0,1002)ener,rhop0_temp
            enddo
         endif
         close (unitout_rhop0)
         endif

!======================================================================
! Gives the angular behavior of rhop0 for each parameter set, as well as
! mean, and sigmas, rhop0at the resonance-peak energy.
!======================================================================
      if (observable.eq.9) then
         open (unitout_rhop0ang, file = outname_rhop0ang)
         do iloop3 = 1, nprecision-1
            xcos = -1.d0 + 2.d0*iloop3/dble(nprecision)
               ener = (xmassr**2-xmpro**2)/(2.d0*xmpro)
               tman =  (-2.d0*ener**2*xmpro + xmjpsi**2*xmpro + 
     -    	ener*(xmjpsi**2 + xcos*dsqrt(xmpro*(2.d0*ener + xmpro))
     -        	*dsqrt((xmjpsi**4 + 4.d0*ener**2*xmpro**2 - 
     -          4.d0*xmjpsi**2*xmpro*(ener + xmpro))/
     -		(xmpro*(2.d0*ener + xmpro)))))/
     -  	(2.d0*ener + xmpro)
               if (spin.eq.'32'.and.smearing.eq.0.d0) then
                  rhop0ang_temp = rhop032(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'32'.and.smearing.gt.0.d0) then
                  rhop0ang_temp = rhop032_smear(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52'.and.smearing.eq.0.d0) then
                  rhop0ang_temp = rhop052(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52'.and.smearing.gt.0.d0) then
                  rhop0ang_temp = rhop052_smear(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_rhop0ang,1003) xcos,tman
     c			,ener,rhop0ang_temp
         enddo
         close (unitout_rhop0ang)
      endif

!======================================================================
! Gives rhopm for each parameter set, as well as mean, and sigmas, as
! a function of energies close to threshold.
!======================================================================
      if (observable.eq.10) then
         open (unitout_rhopm, file = outname_rhopm)
         tman = 1.d0
         nzoo = 2
         dran = 3.d0
         if (smearing.eq.0.d0) then
            do iloop1 = 1, nzoo*nprecision
               ener=(enerpeak-1.5d0)+dran
     c		*iloop1/dble(nzoo*nprecision)
               if (spin.eq.'32') then
                  rhopm_temp = rhopm32(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52') then
                  rhopm_temp = rhopm52(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_rhopm,1002)ener,rhopm_temp
            enddo
            do iloop1 = 1, nprecision
               ener=(enerpeak+1.5d0)+(enermax-(enerpeak+1.5d0))
     c		*iloop1/dble(nprecision)
               if (spin.eq.'32') then
                  rhopm_temp = rhopm32(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52') then
                  rhopm_temp = rhopm52(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_rhopm,1002)ener,rhopm_temp
            enddo
	 elseif (smearing.gt.0.d0) then
            do iloop1 = 1, nzoo*nprecision
               ener=(enerpeak-1.5d0)+dran*iloop1/dble(nzoo*nprecision)
               if (spin.eq.'32') then
                  rhopm_temp = rhopm32_smear(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52') then
                  rhopm_temp = rhopm52_smear(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_rhopm,1002)ener,rhopm_temp
            enddo
         endif
         close (unitout_rhopm)
         endif

!======================================================================
! Gives the angular behavior of rhopm for each parameter set, as well as
! mean, and sigmas, rhopmat the resonance-peak energy.
!======================================================================
      if (observable.eq.11) then
         open (unitout_rhopmang, file = outname_rhopmang)
         do iloop3 = 1, nprecision-1
            xcos = -1.d0 + 2.d0*iloop3/dble(nprecision)
               ener = (xmassr**2-xmpro**2)/(2.d0*xmpro)
               tman =  (-2.d0*ener**2*xmpro + xmjpsi**2*xmpro + 
     -    	ener*(xmjpsi**2 + xcos*dsqrt(xmpro*(2.d0*ener + xmpro))
     -        	*dsqrt((xmjpsi**4 + 4.d0*ener**2*xmpro**2 - 
     -          4.d0*xmjpsi**2*xmpro*(ener + xmpro))/
     -		(xmpro*(2.d0*ener + xmpro)))))/
     -  	(2.d0*ener + xmpro)
               if (spin.eq.'32'.and.smearing.eq.0.d0) then
                  rhopmang_temp = rhopm32(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'32'.and.smearing.gt.0.d0) then
                  rhopmang_temp = rhopm32_smear(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52'.and.smearing.eq.0.d0) then
                  rhopmang_temp = rhopm52(ener
     c		,tman,xgamr,xmassr)
               elseif (spin.eq.'52'.and.smearing.gt.0.d0) then
                  rhopmang_temp = rhopm52_smear(ener
     c		,tman,xgamr,xmassr)
               endif
               write (unitout_rhopmang,1003) xcos,tman
     c			,ener,rhopmang_temp
         enddo
         close (unitout_rhopmang)
      endif
!======================================================================
! Formats
!======================================================================
      stop
 1002 format (9F15.8)
 1003 format (10F15.8)
      end program prepare_plot
