!     -----------------------------------------------------------------------
      module Constants
      real    (8), parameter :: mpi=0.1380386666666667d0
      real    (8), parameter :: momega0=0.78257d0, mphi0=1.01946d0
      end module Constants
!     =======================================================================
      program ktamplitude
      use Constants
      implicit none
      real    (8) :: s, t, u, x, y, param, mv0, norm0
      real    (8) :: s_temp, t_temp, u_temp, mv_temp, norm_temp         !temp s0,t0,u0 for output
      complex (8) :: result
	  integer     :: mvcase, threebody
	  

      ! for Dalitz plot
      real (8)           :: smin, smax, ymin, ymax
      real (8)           :: step_s, step_t, step_x, step_y
      real (8), external :: tmin, tmax, xmin, xmax, norm
      integer            :: i, k, ndalitz

      integer :: infile=100, outfile=101
      integer :: DalitzfileST=11, DalitzfileXY=12
      open(infile, file='param.txt',form='formatted',status ='old')
      open(outfile,file='output.txt',form='formatted',status ='unknown')
	  open(DalitzfileST, file='DalitzST.dat',form='formatted')
	  open(DalitzfileXY, file='DalitzXY.dat',form='formatted')


      read(infile,*) mvcase, threebody, param, s_temp,t_temp,u_temp ! with/ without 3b rescattering effects, param - parameter that is responsible for the shape of Dalitz plot
      call upload_omegaphi(threebody) ! Upload amplitudes for different MV=[0.6,1.11] and energy and then interpolate

!     -----------------------------------------------------------------------
      
      call amp(result,s_temp,t_temp,u_temp,param) ! result - is F(s,t,u) from Eq.(37) (arXiv:1409.7708)
                                      ! to get the intensity one has to multiply it by the p-wave phase space P(s,t,u) from Eq.(37)
      mv_temp=dsqrt(s_temp+t_temp+u_temp-3.d0*mpi**2)                                
 	  norm_temp=norm(mv_temp,param)
      write(*,*) 's=',s_temp,'t=',t_temp,'u=',u_temp
      write(*,*) 'mv=',mv_temp, 'par=',param      
      write(*,*) 'F(s,t,u)=',result/dsqrt(norm_temp)
      write(outfile,*) result/dsqrt(norm_temp)


      If (mvcase == 0) then 
       mv0=momega0
      else if (mvcase == 1) then
       mv0=mphi0
      else 
       write(*,*) 'wrong mvcase'
      end if
      norm0=norm(mv0,param)
!     -----------------------------------------------------------------------
!	  Dalitz plot (s,t) for fixed mv
!     -----------------------------------------------------------------------      
      ndalitz=100
      smin=4*mpi**2
      smax=(mv0-mpi)**2
      step_s = (smax - smin)/dble(ndalitz)
      s      = smin
      do i=1,ndalitz
         step_t=(tmax(s,mv0)-tmin(s,mv0))/dble(ndalitz)
         t=tmin(s,mv0)
         do k=1,ndalitz
          u=3.d0*mpi**2+mv0**2-s-t
		  call amp(result,s,t,u,param)
          write(DalitzfileST,*) s,t,abs(result)**2/norm0
          t = t + step_t
         end do
         write(DalitzfileST,*) ''
         s = s + step_s
      end do
      close(DalitzfileST)

      
!     -----------------------------------------------------------------------
!	  Dalitz plot (x,y) for fixed mv
!     -----------------------------------------------------------------------      
      ndalitz=100
      ymin=-1.d0
      ymax=(mv0 + 3.d0*mpi)/(2.d0*mv0)
      step_y = (ymax - ymin)/dble(ndalitz)
      y      = ymin
      do i=1,ndalitz
         step_x=(xmax(y,mv0)-xmin(y,mv0))/dble(ndalitz)
         x=xmin(y,mv0)
         do k=1,ndalitz  
          s= 1/3.d0*(mv0**2*(1.d0 - 2.d0*y)
     $       + 6.d0*mv0*mpi*y + 3.d0*mpi**2)
          t= 1/3.d0*(mv0**2*(dsqrt(3.d0)*x + y + 1.d0)
     $       - 3.d0*mv0*mpi*(dsqrt(3.d0)*x + y) + 3*mpi**2)
          u=3.d0*mpi**2+mv0**2-s-t
		  call amp(result,s,t,u,param)
          write(DalitzfileXY,*) x,y,abs(result)**2/norm0
          x = x + step_x
         end do
         write(DalitzfileXY,*) ''
         y = y + step_y
      end do
      close(DalitzfileXY)

!     Make the plots
!      call SYSTEM('/usr/local/bin/gnuplot scriptplotxy.txt')


      
      stop
      end program ktamplitude

!     -----------------------------------------------------------------------
!     Setting up the normalization at x=y=0 for any mv, par
!     -----------------------------------------------------------------------
	  real (8) function norm(mv,par)
	  use Constants
	  implicit none 
      real (8) :: s,t,u,x,y,mv,par
      complex (8) :: res
	  y=0.d0
	  x=0.d0
          s= 1/3.d0*(mv**2*(1.d0 - 2.d0*y)
     $       + 6.d0*mv*mpi*y + 3.d0*mpi**2)
          t= 1/3.d0*(mv**2*(dsqrt(3.d0)*x + y + 1.d0)
     $       - 3.d0*mv*mpi*(dsqrt(3.d0)*x + y) + 3*mpi**2)
          u=3.d0*mpi**2+mv**2-s-t
          
	  call amp(res,s,t,u,par)
	  norm=abs(res)**2
      end function norm
!     -----------------------------------------------------------------------



!     -----------------------------------------------------------------------
!	  Boundaries of Dalitz plot
!     -----------------------------------------------------------------------
      real (8) function tmin(s,mv0)
      use Constants
      implicit none
      real (8) :: s, mv0
      tmin=-(dsqrt(-mpi**2+(mv0**2-mpi**2-s)**2/(4.d0*s))
     $	   + dsqrt(-mpi**2+s/4.d0))**2
     $     + ((mv0**2-mpi**2-s)/(2.d0*dsqrt(s))+dsqrt(s)/2.d0)**2
      end function tmin
!     -----------------------------------------------------------------------
      real (8) function tmax(s,mv0)
      use Constants
      implicit none
      real (8) :: s, mv0
      tmax=-(-dsqrt(-mpi**2+(mv0**2-mpi**2-s)**2/(4.d0*s))
     $     + dsqrt(-mpi**2+s/4.d0))**2
     $     + ((mv0**2-mpi**2-s)/(2.d0*dsqrt(s))+dsqrt(s)/2.d0)**2
      end function tmax


!     -----------------------------------------------------------------------
      real (8) function xmin(y,mv0)
      use Constants
      implicit none
      real (8) :: y, mv0
      xmin=-((dsqrt(-9*mpi**2 + mv0*(mv0 - 2*mv0*y 
     $       + 6*mpi*y))*
     $      dsqrt((mv0**2*(mv0 - 3*mpi)*(1 + y)*
     $          (mv0 + 3*mpi + mv0*y - 3*mpi*y))/
     $        (3*mpi**2 + mv0*(mv0 - 2*mv0*y 
     $        + 6*mpi*y))))/
     $    (dsqrt(3.d0)*mv0*(mv0 - 3*mpi)))
      end function xmin
      
!     -----------------------------------------------------------------------
      real (8) function xmax(y,mv0)
      use Constants
      implicit none
      real (8) :: y, mv0
      xmax=(dsqrt(-9*mpi**2 + mv0*(mv0 - 2*mv0*y
     $       + 6*mpi*y))*
     $      dsqrt((mv0**2*(mv0 - 3*mpi)*(1 + y)*
     $          (mv0 + 3*mpi + mv0*y - 3*mpi*y))/
     $        (3*mpi**2 + mv0*(mv0 - 2*mv0*y 
     $       + 6*mpi*y))))/
     $    (dsqrt(3.d0)*mv0*(mv0 - 3*mpi))
      end function xmax




!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine amp(ret,s,t,u,param)
      use Constants
      implicit none
!     input
      real (8) :: s,t,u,param
!     output
      complex (8) :: ret
!     internal variables
      real (8) :: mv
      complex (8) :: f0s,f1s,f0u,f1u,f0t,f1t

      mv = dsqrt(s+t+u-3.d0*mpi**2)

      call f0f1(s,mv,f0s,f1s)
      call f0f1(u,mv,f0u,f1u)
      call f0f1(t,mv,f0t,f1t)
 
      ret  = f0s+f0t+f0u+ param*(f1s + f1t + f1u) 
      return
      end subroutine amp

      module KTmodel
      complex (8), parameter :: xr=cmplx(1.d0,0.d0), xi=cmplx(0.d0,1.d0)
      integer, parameter     :: nintpl1 = 201, nv = 52 !related to step and mvmin
      real (8), parameter    :: step = 0.010d0, mvmin = 0.6d0

      real (8), dimension (nintpl1) :: sn
      real (8), dimension (nintpl1,nv) ::  storeref0, storebref0
      real (8), dimension (nintpl1,nv) :: storecref0, storedref0
      real (8), dimension (nintpl1,nv) ::  storeimf0, storebimf0
      real (8), dimension (nintpl1,nv) :: storecimf0, storedimf0
      real (8), dimension (nintpl1,nv) ::  storeref1, storebref1
      real (8), dimension (nintpl1,nv) :: storecref1, storedref1
      real (8), dimension (nintpl1,nv) ::  storeimf1, storebimf1
      real (8), dimension (nintpl1,nv) :: storecimf1, storedimf1

      save sn
      save storeref0, storebref0, storecref0, storedref0
      save storeimf0, storebimf0, storecimf0, storedimf0
      save storeref1, storebref1, storecref1, storedref1
      save storeimf1, storebimf1, storecimf1, storedimf1
      end module KTmodel
      
!     -----------------------------------------------------------------------
!     Upload files
!     -----------------------------------------------------------------------
      subroutine upload_omegaphi(threebody)
      use KTmodel
      implicit none
      integer :: threebody
      character (10) :: genericname
      character (1)  :: namethreebody,numeric1,numeric2
      character (4)  :: extension
      character (17), dimension(100) :: filename ! 17=10+1+1+1+4; 100 (file names) related to i,j
      integer :: i,j,mvindex,ind
      real (8), dimension (nintpl1) :: ref0, bref0, cref0, dref0
      real (8), dimension (nintpl1) :: imf0, bimf0, cimf0, dimf0
      real (8), dimension (nintpl1) :: ref1, bref1, cref1, dref1
      real (8), dimension (nintpl1) :: imf1, bimf1, cimf1, dimf1

!     generate the name of the import files (see below open(...) )
      genericname = 'ListF_iter'
      if (threebody == 1) then
         namethreebody = char(51) !51 -> 3
      else if (threebody == 0) then
         namethreebody = char(48) !48 -> 0
      else
         write(*,*) 'wrong threebody case'   
      endif
      extension = '.dat'
      ind = 0
      do i=0,9
         do j=0,9
            ind = ind + 1
            numeric1 = char(48+i)
            numeric2 = char(48+j)
            filename(ind) = genericname // namethreebody // numeric1 
     c           // numeric2 // extension
         enddo
      enddo


!     read import and store the interpolated results
      do mvindex = 1,nv
         open(1,file=filename(mvindex),status='old')
         do i=1, nintpl1
            read(1,*) sn(i),ref0(i),imf0(i),ref1(i),imf1(i)
         enddo
         close(1)
!     interpolation (determination of the coefficients)
         call spline(nintpl1,sn,ref0,bref0,cref0,dref0) 
         call spline(nintpl1,sn,imf0,bimf0,cimf0,dimf0) 
         call spline(nintpl1,sn,ref1,bref1,cref1,dref1) 
         call spline(nintpl1,sn,imf1,bimf1,cimf1,dimf1) 
         
         do i=1,nintpl1
            storeref0(i,mvindex)  =  ref0(i)
            storebref0(i,mvindex) = bref0(i)
            storecref0(i,mvindex) = cref0(i)
            storedref0(i,mvindex) = dref0(i)
            
            storeimf0(i,mvindex)  =  imf0(i)
            storebimf0(i,mvindex) = bimf0(i)
            storecimf0(i,mvindex) = cimf0(i)
            storedimf0(i,mvindex) = dimf0(i)
            
            storeref1(i,mvindex)  =  ref1(i)
            storebref1(i,mvindex) = bref1(i)
            storecref1(i,mvindex) = cref1(i)
            storedref1(i,mvindex) = dref1(i)
            
            storeimf1(i,mvindex)  =  imf1(i)
            storebimf1(i,mvindex) = bimf1(i)
            storecimf1(i,mvindex) = cimf1(i)
            storedimf1(i,mvindex) = dimf1(i)
         enddo
      enddo

      return
      end subroutine upload_omegaphi


!     -----------------------------------------------------------------------
      subroutine f0f1(s,mv,f0,f1)
      use KTmodel
      implicit none
!     input
      real (8) :: s,mv
!     output
      complex (8) :: f0,f1
!     internal variables
      integer :: i,mvindex
      real (8), dimension (nintpl1) :: ref0, bref0, cref0, dref0
      real (8), dimension (nintpl1) :: imf0, bimf0, cimf0, dimf0
      real (8), dimension (nintpl1) :: ref1, bref1, cref1, dref1
      real (8), dimension (nintpl1) :: imf1, bimf1, cimf1, dimf1
      real (8) :: ref0sd,imf0sd,ref1sd,imf1sd
      real (8) :: ref0su,imf0su,ref1su,imf1su
      real (8) :: mvd,mvu
      complex (8) :: f0su,f0sd,f1su,f1sd
!     external functions
      real(8), external :: interpolate

      mvindex = 1+int((mv-mvmin)/step)
      mvd = mvmin + dble(mvindex-1)*step
      mvu = mvmin + dble(mvindex)*step

      do i=1,nintpl1
          ref0(i) = storeref0(i,mvindex)  
         bref0(i) = storebref0(i,mvindex) 
         cref0(i) = storecref0(i,mvindex) 
         dref0(i) = storedref0(i,mvindex) 
         
          imf0(i) = storeimf0(i,mvindex) 
         bimf0(i) = storebimf0(i,mvindex) 
         cimf0(i) = storecimf0(i,mvindex) 
         dimf0(i) = storedimf0(i,mvindex)
         
          ref1(i) = storeref1(i,mvindex)  
         bref1(i) = storebref1(i,mvindex) 
         cref1(i) = storecref1(i,mvindex)
         dref1(i) = storedref1(i,mvindex)
          
          imf1(i) = storeimf1(i,mvindex) 
         bimf1(i) = storebimf1(i,mvindex)
         cimf1(i) = storecimf1(i,mvindex) 
         dimf1(i) = storedimf1(i,mvindex)         
      enddo

      ref0sd = interpolate(nintpl1,s,sn,ref0,bref0,cref0,dref0)
      imf0sd = interpolate(nintpl1,s,sn,imf0,bimf0,cimf0,dimf0)
      ref1sd = interpolate(nintpl1,s,sn,ref1,bref1,cref1,dref1)
      imf1sd = interpolate(nintpl1,s,sn,imf1,bimf1,cimf1,dimf1)

      f0sd = ref0sd+xi*imf0sd
      f1sd = ref1sd+xi*imf1sd

      mvindex = mvindex + 1
      do i=1,nintpl1
          ref0(i) =  storeref0(i,mvindex)  
         bref0(i) = storebref0(i,mvindex) 
         cref0(i) = storecref0(i,mvindex) 
         dref0(i) = storedref0(i,mvindex) 
         
          imf0(i) =  storeimf0(i,mvindex) 
         bimf0(i) = storebimf0(i,mvindex) 
         cimf0(i) = storecimf0(i,mvindex) 
         dimf0(i) = storedimf0(i,mvindex)
         
          ref1(i) =  storeref1(i,mvindex)  
         bref1(i) = storebref1(i,mvindex) 
         cref1(i) = storecref1(i,mvindex)
         dref1(i) = storedref1(i,mvindex)
          
          imf1(i) =  storeimf1(i,mvindex) 
         bimf1(i) = storebimf1(i,mvindex)
         cimf1(i) = storecimf1(i,mvindex) 
         dimf1(i) = storedimf1(i,mvindex)         
      enddo

      ref0su = interpolate(nintpl1,s,sn,ref0,bref0,cref0,dref0)
      imf0su = interpolate(nintpl1,s,sn,imf0,bimf0,cimf0,dimf0)
      ref1su = interpolate(nintpl1,s,sn,ref1,bref1,cref1,dref1)
      imf1su = interpolate(nintpl1,s,sn,imf1,bimf1,cimf1,dimf1)

      f0su = ref0su+xi*imf0su
      f1su = ref1su+xi*imf1su

      f0 = (f0su-f0sd)*mv/step + (f0sd*mvu-f0su*mvd)/step
      f1 = (f1su-f1sd)*mv/step + (f1sd*mvu-f1su*mvd)/step

      return
      end subroutine f0f1
!     -----------------------------------------------------------------------
!     interpolation subroutine
!     -----------------------------------------------------------------------
      subroutine spline(n,x,y,b,c,d)
      implicit double precision (a-h,o-z)

      real (8) x(n),y(n),b(n),c(n),d(n)

      nm1 = n-1
      if(n.lt.2) return
      if(n.lt.3) goto 50

      d(1)=x(2)-x(1)
      c(2)=(y(2)-y(1))/d(1)
      do 10 i=2,nm1
         d(i)=x(i+1)-x(i)
         b(i)=2.0*(d(i-1)+d(i))
         c(i+1)=(y(i+1)-y(i))/d(i)
         c(i)=c(i+1)-c(i)
 10   continue
      
      b(1)=-d(1)
      b(n)=-d(n-1)
      c(1)=0.0
      c(n)=0.0
      if(n.eq.3) goto 15
      c(1)=c(3)/(x(4)-x(2))-c(2)/(x(3)-x(1))
      c(n)=c(n-1)/(x(n)-x(n-2))-c(n-2)/(x(n-1)-x(n-3))
      c(1)=c(1)*d(1)**2/(x(4)-x(1))
      c(n)=-c(n)*d(n-1)**2/(x(n)-x(n-3))
      
 15   do 20 i=2,n
         t=d(i-1)/b(i-1)
         b(i)=b(i)-t*d(i-1)
         c(i)=c(i)-t*c(i-1)
 20   continue
      
      c(n)=c(n)/b(n)
      do 30 ib=1,nm1
         i=n-ib
         c(i)=(c(i)-d(i)*c(i+1))/b(i)
 30   continue

      b(n)=(y(n)-y(nm1))/d(nm1)+d(nm1)*(c(nm1)+2.0*c(n))
      do 40 i=1,nm1
         b(i)=(y(i+1)-y(i))/d(i)-d(i)*(c(i+1)+2.0*c(i))
         d(i)=(c(i+1)-c(i))/d(i)
         c(i)=3.0*c(i)
 40   continue
      c(n)=3.0*c(n)
      d(n)=d(n-1)
      return
        
 50   b(1)=(y(2)-y(1))/(x(2)-x(1))
      c(1)=0.0
      d(1)=0.0
      b(2)=b(1)
      c(2)=0.0
      d(2)=0.0
      return
      end
      
      real (8) function interpolate(n,u,x,y,b,c,d)
     *     result (output)
      implicit double precision (a-h,o-z)
      real (8) x(n),y(n),b(n),c(n),d(n)
      data i /1/
      
      if(i.ge.n) i=1
      if(u.lt.x(i)) goto 10
      if(u.le.x(i+1)) goto 30
      
 10   i=1
      j=n+1
 20   k=(i+j)/2
      if(u.lt.x(k)) j=k
      if(u.ge.x(k)) i=k
      if(j.gt.i+1) goto 20
      
 30   dx=u-x(i)
      output=y(i)+dx*(b(i)+dx*(c(i)+dx*d(i)))
      return
      end
      
