!     --------------------------------------------------------------------------------------
!     --------------------------------------------------------------------------------------
!
!     Total Cross Section calculator for KN scattering
!     
!     by the Joint Physics Analysis Center (JPAC)
!     
!     Version: September 1, 2015
!
!     --------------------------------------------------------------------------------------
!     --------------------------------------------------------------------------------------
!
!     Authors:
!
!     Cesar Fernandez-Ramirez (Jefferson Lab)
!     Igor V. Danilkin (Jefferson Lab)
!     Vincent Mathieu (Indiana University)
!     Adam P. Szczepaniak (Indiana University and Jefferson Lab)
!
!     Citation: Fernandez-Ramirez et al., arxiv:
!
!     First version: Cesar Fernandez-Ramirez (Jefferson Lab)
!     This version: Cesar Fernandez-Ramirez (Jefferson Lab)
!
!     Contact: cefera@gmail.com (Cesar Fernandez-Ramirez)
!
!     Disclaimers:
!
!     1 - This code follows the 'garbage in, garbage out' philosphy. If your
!     parameters do not make sense, the output will not make sense either.
!     2 - You can use, share and modify this code under your own responsability.
!     3 - This code is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
!     4 - No PhD students or postdocs were severely damaged during the development
!     of this project. 
!
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      module constantes
      integer, parameter :: max_coefs = 558
      integer, parameter :: ncoefs = max_coefs
      end module constantes
!     -----------------------------------------------------------------------
      module coeficientes
      use constantes
      real (8), dimension (ncoefs) :: coeficiente
      save coeficiente
      end module coeficientes
!     -----------------------------------------------------------------------
      module matrixsizes
      integer, parameter :: matrixsize = 16
      end module matrixsizes
!     -----------------------------------------------------------------------
      module resonancesizes
      integer, parameter :: nresonances = 59
      integer, parameter :: nstables = 14
      integer, parameter :: ndecays  = 16
      integer, parameter :: ndecaysm = 6*ndecays, ndecaysw = 2*ndecays
      integer, parameter :: max_data_points = 1000
      integer, parameter :: max_archivos = 12
      end module resonancesizes
!     -----------------------------------------------------------------------
      module canalesop
      integer, parameter :: nres = 6
      integer, dimension(nres):: cs01,cp01,cp03,cd03,cd05,cf05,cf07,cg07
      integer, dimension(nres):: cs11,cp11,cp13,cd13,cd15,cf15,cf17,cg17
      save cs01,cp01,cp03,cd03,cd05,cf05,cf07,cg07
      save cs11,cp11,cp13,cd13,cd15,cf15,cf17,cg17
      end module canalesop
!     -----------------------------------------------------------------------
      module matrixstorage
      use resonancesizes
      use matrixsizes
      integer,  dimension (nresonances) :: momangresonancias
      real (8), dimension (nstables) :: masas,widths
      real (8), dimension (nresonances) :: masasresonancias,isresonancia
      real (8), dimension (6*ndecays) :: channelmass
      real (8), dimension (2*ndecays) :: channelwidths
      real (8), dimension(nresonances,ndecays) :: branchings,phimres
      real (8), dimension(ndecays,max_data_points) :: phiphasestorage
      complex (8), dimension(ndecays,max_data_points) :: phistorage
      complex (8), dimension (matrixsize,matrixsize) :: id,scatmat
      save momangresonancias
      save masas,widths
      save masasresonancias,attraction,isresonancia
      save channelmass
      save channelwidths
      save branchings,phimres
      save phistorage
      save phiphasestorage
      save id,scatmat
      end module matrixstorage
!      -----------------------------------------------------------------------
      module superdata
      use resonancesizes
      integer :: indice_maximo
      character (3) :: pwchoice
      real (8), dimension (max_data_points) :: svector
      save indice_maximo
      save pwchoice
      save svector
      end module superdata
!     -----------------------------------------------------------------------
      module channelassignment
      integer, parameter :: jprot       = 1
      integer, parameter :: jkaon       = 2
      integer, parameter :: jpion       = 3
      integer, parameter :: jeta        = 4
      integer, parameter :: jkstar      = 5
      integer, parameter :: jdelta      = 6
      integer, parameter :: jlambda     = 7
      integer, parameter :: jlambda1520 = 8
      integer, parameter :: jsigma      = 9
      integer, parameter :: jsigma1385  = 10
      integer, parameter :: jdummy      = 11
      integer, parameter :: jpion0      = 12
      integer, parameter :: jkav        = 13
      integer, parameter :: jnucl       = 14
!     active channels
      integer, parameter :: kbarn         = 1
      integer, parameter :: pisigma       = 2
      integer, parameter :: pilambda      = 3
      integer, parameter :: etalambda     = 4
      integer, parameter :: etasigma      = 5
      integer, parameter :: kbar1starn    = 6
      integer, parameter :: kbar3starnm   = 7
      integer, parameter :: kbar3starnp   = 8
      integer, parameter :: pisigma1385p  = 9
      integer, parameter :: pisigma1385m  = 10
      integer, parameter :: kbardeltam    = 11
      integer, parameter :: kbardeltap    = 12
      integer, parameter :: pilambda1520m = 13
      integer, parameter :: pilambda1520p = 14
      integer, parameter :: dummylambda   = 15
      integer, parameter :: dummysigma    = 16
      end module channelassignment
!     -----------------------------------------------------------------------
      module usefulnumbers
      real (8), parameter :: pi = 4.d0*datan(1.d0)
      real (8), parameter :: milibarn = 10.d0/25.68189994d0
      real (8), parameter :: sq2 = dsqrt(2.d0), epsilon = 1.d-20
      real (8), parameter :: s0 = 1.d0, w00 = 2.d0, cutoff = 10.d0
      real (8), parameter :: fermi2 = (1.d0/0.19733d0)**2
      complex (8), parameter :: xzero = dcmplx(0.d0,0.d0)
      complex (8), parameter :: xr = dcmplx(1.d0,0.d0)
      complex (8), parameter :: xi = dcmplx(0.d0,1.d0)
      end module usefulnumbers
!     -----------------------------------------------------------------------
      module masses             ! gev
      real (8), parameter :: mphoton = 0.d0, gphoton  = 0.d0 ! ok
      real (8), parameter :: mg      = mphoton
      real (8), parameter :: mpi     = 0.139d0,   gpi     = 0.d0 ! ok
      real (8), parameter :: mpi0    = 0.13498d0, gpi0    = 0.d0 ! ok
      real (8), parameter :: mkm     = 0.49368d0, gkm     = 0.d0 ! ok
      real (8), parameter :: mk      = 0.49368d0, gk      = 0.d0 ! ok
      real (8), parameter :: meta    = 0.54751d0, geta    = 0.d0 ! ok
      real (8), parameter :: mkstar  = 0.8939d0,  gkstar  = 0.0508d0 ! ok 
      real (8), parameter :: metap   = 0.958d0,   getap   = 0.d0
      real (8), parameter :: mrho    = 0.771d0,   grho    = 0.1492d0 ! ok
      real (8), parameter :: momega  = 0.7826d0,  gomega  = 0.d0 ! ok
      real (8), parameter :: mdummy  = 2.d0*mpi,  gdummy  = 0.800d0 ! ok
      real (8), parameter :: mp      = 0.93827d0, gprot   = 0.d0 ! ok
      real (8), parameter :: mnuc    = 0.939d0,   gnuc    = 0.d0 ! ok
      real (8), parameter :: mlambda = 1.11568d0, glambda = 0.d0 ! ok 
      real (8), parameter :: mdelta  = 1.232d0,   gdelta  = 0.120d0 ! ok
      real (8), parameter :: mroper  = 1.440d0,   groper  = 0.200d0
      real (8), parameter :: msigma  = 1.19264d0, gsigma  = 0.d0 ! ok
      real (8), parameter :: msigma1385 = 1.3837d0, gsigma1385 = 0.036d0 ! ok
      real (8), parameter :: mlambda1520= 1.5195d0, glambda1520=0.0156d0 ! ok
      end module masses
!     -----------------------------------------------------------------------
      module matrixoperations
      implicit none
      contains
!     -----------------------------------------------------------------------
      function smatrix(id,rhophase,tmat)
      use usefulnumbers
      use matrixsizes
      implicit none
      complex (8), dimension(matrixsize,matrixsize) :: id,rhophase,tmat
      complex (8), dimension(matrixsize,matrixsize) :: smatrix
      smatrix = id + 2.d0*xi*matmul(rhophase,matmul(tmat,rhophase))
      return
      end function smatrix
!     -----------------------------------------------------------------------
      function rhomatrix(indice)
      use usefulnumbers
      use matrixstorage
      implicit none
      integer :: indice,j,k
      complex (8), dimension (matrixsize,matrixsize) :: rhomatrix
      do j=1,matrixsize
         do k=j,matrixsize
            if (j.eq.k) then
               rhomatrix(j,k) = phistorage(k,indice)
            else
               rhomatrix(k,j) = xzero
               rhomatrix(j,k) = xzero
            endif
         enddo
      enddo
      return
      end function rhomatrix
!     -----------------------------------------------------------------------
      function rhomatrixphase(indice)
      use usefulnumbers
      use matrixstorage
      implicit none
      integer :: indice,j,k
      complex (8), dimension (matrixsize,matrixsize) :: rhomatrixphase
      do j=1,matrixsize
         do k=j,matrixsize
            if (j.eq.k) then
               rhomatrixphase(j,k) = xr*dsqrt(phiphasestorage(j,indice))
            else
               rhomatrixphase(k,j) = xzero
               rhomatrixphase(j,k) = xzero
            endif
         enddo
      enddo
      return
      end function rhomatrixphase
!     -----------------------------------------------------------------------
      function tr6fast(s,indice,a1,a2,a3,a4,a5,a6)
      use matrixsizes
      implicit none
      real (8), parameter :: units1 = 1000.d0
      integer :: j,k,a1,a2,a3,a4,a5,a6,indice
      real (8) :: s
      complex (8) :: x11,x22,x33,x44,x55,x66,x12,x13,x14,x15,x16,
     c           x23,x24,x25,x26,x34,x35,x36,x45,x46,x56
      complex (8) :: a11,a22,a33,a44,a55,a66,a12,a13,a14,a15,a16,
     c     a23,a24,a25,a26,a34,a35,a36,a45,a46,a56,d6
      complex (8), dimension (matrixsize,matrixsize) :: tr6fast
      call amatrix(s,indice,a1,a2,a3,a4,a5,a6,
     c     a11,a22,a33,a44,a55,a66,a12,a13,a14,a15,a16,
     c     a23,a24,a25,a26,a34,a35,a36,a45,a46,a56,d6)
      do j=1,matrixsize
         do k=j,matrixsize
            call xxmatrix(j,k,a1,a2,a3,a4,a5,a6,
     c           x11,x22,x33,x44,x55,x66,x12,x13,x14,x15,x16,
     c           x23,x24,x25,x26,x34,x35,x36,x45,x46,x56)
            tr6fast(j,k) = ( ( a11*x11 
     c           + a22*x22 + a33*x33 + a44*x44 + a55*x55 + a66*x66
     c           + a12*x12 + a13*x13 + a14*x14 + a15*x15 + a16*x16
     c           + a23*x23 + a24*x24 + a25*x25 + a26*x26
     c           + a34*x34 + a35*x35 + a36*x36
     c           + a45*x45 + a46*x46
     c           + a56*x56 )/d6 )/units1
            tr6fast(k,j) = tr6fast(j,k)
         enddo
      enddo
      return
      end function tr6fast
      end module matrixoperations
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      program xsecef
      use usefulnumbers
      use channelassignment
      use matrixstorage
      use superdata
      implicit none
      character (15) :: name
      character (11) :: process
      integer :: i,j,i1,i2,index,kinchoice,npoints
      real (8) :: sini,sfin,elabini,elabfin,plabini,plabfin,step
      real (8) :: s,m1,m2,elab,plab,q2
      real (8) :: kiniinput,kfininput,c0,c1,output
      call initall
      open (2,file='file.inp',status='old')
      read (2,*)process
      read (2,*)kinchoice
      read (2,*)kiniinput
      read (2,*)kfininput
      read (2,*)npoints
      close (2)
      indice_maximo = 1
      index = indice_maximo
      i1 =  kbarn
      m1 = masas(jkaon)
      m2 = masas(jprot)
      open (5,file='title.txt')
      select case (process)
      case ('prk-toprok-')
         c0 =  0.5d0
         c1 =  0.5d0
         i2 =  kbarn
         write(5,*)'''K^- p -> K^- p'''
      case ('prk-toneuk0')
         c0 = -0.5d0
         c1 =  0.5d0
         i2 =  kbarn
         write(5,*)'''K^- p -> barK^0 n'''
      case ('prk-tolapi0')
         c0 =  0.d0
         c1 =  1.d0/dsqrt(2.d0)
         i2 =  pilambda
         write(5,*)'''K^- p -> pi^0 Lambda'''
      case ('prk-tos+pi-')
         c0 = -1.d0/dsqrt(6.d0) 
         c1 = -0.5d0
         i2 =  pisigma
         write(5,*)'''K^- p -> pi^- Sigma^+'''
      case ('prk-tos0pi0')
         c0 =  1.d0/dsqrt(6.d0)
         c1 =  0.d0
         i2 =  pisigma
         write(5,*)'''K^- p -> pi^0 Sigma^0'''
      case ('prk-tos-pi+')
         c0 = -1.d0/dsqrt(6.d0)
         c1 =  0.5d0
         i2 =  pisigma
         write(5,*)'''K^- p -> pi^+ Sigma^-'''
      end select
      close (5)
      select case(kinchoice)
      case (1)                  ! input is s
         sini    = kiniinput
         sfin    = kfininput
      case (2)                  ! input is plab
         plabini = kiniinput
         plabfin = kfininput
         elabini = dsqrt(plabini**2 + m1**2)
         elabfin = dsqrt(plabfin**2 + m1**2)
         sini    = 2.d0*elabini*m2 + m1**2 + m2**2
         sfin    = 2.d0*elabfin*m2 + m1**2 + m2**2
      case (3)                  ! input is elab
         elabini = kiniinput
         elabfin = kfininput
         sini    = 2.d0*elabini*m2 + m1**2 + m2**2
         sfin    = 2.d0*elabfin*m2 + m1**2 + m2**2
      end select
      name = 'output.txt'
      open (1,file=name)
      write (1,'(A1,A14,4A15)')' #','s (GeV^2)','elab (GeV)',
     c     'plab (GeV)','q^2 (GeV^2)','xsec (mb)'
      step = (sfin-sini)/dble(npoints-1)
      s = sini
      do i=1,npoints
         svector(index) = s
         elab = (s - m1**2 - m2**2 )/2.d0/m2
         plab = dsqrt( elab**2 - m1**2 )
         q2   = ( s-( m1+m2 )**2 )*( s - ( m1-m2 )**2 )/4.d0/s
         call calculation_xsec(i1,i2,s,index,c0,c1,output)
         write (1,'(5F15.8)')s,elab,plab,q2,output
         s = s + step
      enddo
      close (1)
      stop
      end program xsecef
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine calculation_xsec(i1,i2,s,index,c0,c1,secef)
      use usefulnumbers
      use channelassignment
      use matrixstorage
      implicit none
      character (5) :: casechannel
      integer :: index,j,i1,i2
      real (8) :: s,m1,m2,c0,c1,amp2,q2,f2,g2,secef
      complex (8) :: f,g
      complex (8), dimension (5) :: f0,f1,g0,g1
      m1 = masas(jkaon)
      m2 = masas(jprot)
      q2 = ( s - ( m1+m2 )**2 )*( s - ( m1-m2 )**2 )/4.d0/s
      call pwcalculator(i1,i2,s,index,f0,f1,g0,g1)
      f2 = 0.d0
      g2 = 0.d0
      do j=1,5
         f = c0*f0(j) + c1*f1(j)
         g = c0*g0(j) + c1*g1(j)
         f2 = f2 +   dble(j)*dble(f*dconjg(f))
         g2 = g2 + dble(j-1)*dble(g*dconjg(g))
      enddo
      amp2 = f2 + g2
      secef = 4.d0*pi*milibarn*amp2/q2
      return
      end subroutine calculation_xsec
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine initall
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: i,j
      call identityinit
      call assignmasses
      call channelsassignment
      do i=1,nresonances
         do j=1,ndecays
            branchings(i,j) = 0.d0
         enddo
      enddo
      do i=1,ncoefs
         coeficiente(i) = 0.d0
      enddo
      call parametersuploader
      return
      end subroutine initall
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine pwcalculator(i1,i2,s,i,f0,f1,g0,g1)
      use usefulnumbers
      use channelassignment
      use matrixstorage
      use canalesop
      implicit none
      character (3) :: pwchoice
      integer :: i,i1,i2,a1,a2,a3,a4,a5,a6,l
      real (8) :: s
      complex (8) :: pws01,pwp01,pwp03,pwd03,pwd05,pwf05,pwf07,pwg07
      complex (8) :: pws11,pwp11,pwp13,pwd13,pwd15,pwf15,pwf17,pwg17
      complex (8), dimension (5) :: f0,f1,g0,g1
      complex (8), external :: partialwave
!     ---------------- s01
      pwchoice = 's01'
      l  = 0
      a1 = cs01(1)
      a2 = cs01(2)
      a3 = cs01(3)
      a4 = cs01(4)
      a5 = cs01(5)
      a6 = cs01(6)
      call phicalculator(l)
      pws01 = partialwave(i1,i2,s,i,pwchoice,a1,a2,a3,a4,a5,a6)
!     ---------------- p01
      pwchoice = 'p01'
      l  = 1
      a1 = cp01(1)
      a2 = cp01(2)
      a3 = cp01(3)
      a4 = cp01(4)
      a5 = cp01(5)
      a6 = cp01(6)
      call phicalculator(l)
      pwp01 = partialwave(i1,i2,s,i,pwchoice,a1,a2,a3,a4,a5,a6)
!     ---------------- p03
      pwchoice = 'p03'
      l  = 1
      a1 = cp03(1)
      a2 = cp03(2)
      a3 = cp03(3)
      a4 = cp03(4)
      a5 = cp03(5)
      a6 = cp03(6)
      call phicalculator(l)
      pwp03 = partialwave(i1,i2,s,i,pwchoice,a1,a2,a3,a4,a5,a6)
!     ---------------- d03
      pwchoice = 'd03'
      l  = 2
      a1 = cd03(1)
      a2 = cd03(2)
      a3 = cd03(3)
      a4 = cd03(4)
      a5 = cd03(5)
      a6 = cd03(6)
      call phicalculator(l)
      pwd03 = partialwave(i1,i2,s,i,pwchoice,a1,a2,a3,a4,a5,a6)
!     ---------------- d05
      pwchoice = 'd05'
      l  = 2
      a1 = cd05(1)
      a2 = cd05(2)
      a3 = cd05(3)
      a4 = cd05(4)
      a5 = cd05(5)
      a6 = cd05(6)
      call phicalculator(l)
      pwd05 = partialwave(i1,i2,s,i,pwchoice,a1,a2,a3,a4,a5,a6)
!     ---------------- f05
      pwchoice = 'f05'
      l  = 3
      a1 = cf05(1)
      a2 = cf05(2)
      a3 = cf05(3)
      a4 = cf05(4)
      a5 = cf05(5)
      a6 = cf05(6)
      call phicalculator(l)
      pwf05 = partialwave(i1,i2,s,i,pwchoice,a1,a2,a3,a4,a5,a6)
!     ---------------- f07
      pwchoice = 'f07'
      l  = 3
      a1 = cf07(1)
      a2 = cf07(2)
      a3 = cf07(3)
      a4 = cf07(4)
      a5 = cf07(5)
      a6 = cf07(6)
      call phicalculator(l)
      pwf07 = partialwave(i1,i2,s,i,pwchoice,a1,a2,a3,a4,a5,a6)
!     ---------------- g07
      pwchoice = 'g07'
      l  = 4
      a1 = cg07(1)
      a2 = cg07(2)
      a3 = cg07(3)
      a4 = cg07(4)
      a5 = cg07(5)
      a6 = cg07(6)
      call phicalculator(l)
      pwg07 = partialwave(i1,i2,s,i,pwchoice,a1,a2,a3,a4,a5,a6)
!     ---------------- s11
      pwchoice = 's11'
      l  = 0
      a1 = cs11(1)
      a2 = cs11(2)
      a3 = cs11(3)
      a4 = cs11(4)
      a5 = cs11(5)
      a6 = cs11(6)
      call phicalculator(l)      
      pws11 = partialwave(i1,i2,s,i,pwchoice,a1,a2,a3,a4,a5,a6)
!     ---------------- p11
      pwchoice = 'p11'
      l  = 1
      a1 = cp11(1)
      a2 = cp11(2)
      a3 = cp11(3)
      a4 = cp11(4)
      a5 = cp11(5)
      a6 = cp11(6)
      call phicalculator(l)
      pwp11 = partialwave(i1,i2,s,i,pwchoice,a1,a2,a3,a4,a5,a6)
!     ---------------- p13
      pwchoice = 'p13'
      l  = 1
      a1 = cp13(1)
      a2 = cp13(2)
      a3 = cp13(3)
      a4 = cp13(4)
      a5 = cp13(5)
      a6 = cp13(6)
      call phicalculator(l)
      pwp13 = partialwave(i1,i2,s,i,pwchoice,a1,a2,a3,a4,a5,a6)
!     ---------------- d13
      pwchoice = 'd13'
      l  = 2
      a1 = cd13(1)
      a2 = cd13(2)
      a3 = cd13(3)
      a4 = cd13(4)
      a5 = cd13(5)
      a6 = cd13(6)
      call phicalculator(l)
      pwd13 = partialwave(i1,i2,s,i,pwchoice,a1,a2,a3,a4,a5,a6)
!     ---------------- d15
      pwchoice = 'd15'
      l  = 2
      a1 = cd15(1)
      a2 = cd15(2)
      a3 = cd15(3)
      a4 = cd15(4)
      a5 = cd15(5)
      a6 = cd15(6)
      call phicalculator(l)
      pwd15 = partialwave(i1,i2,s,i,pwchoice,a1,a2,a3,a4,a5,a6)
!     ---------------- f15
      pwchoice = 'f15'
      l  = 3
      a1 = cf15(1)
      a2 = cf15(2)
      a3 = cf15(3)
      a4 = cf15(4)
      a5 = cf15(5)
      a6 = cf15(6)
      call phicalculator(l)
      pwf15 = partialwave(i1,i2,s,i,pwchoice,a1,a2,a3,a4,a5,a6)
!     ---------------- f17
      pwchoice = 'f17'
      l  = 3
      a1 = cf17(1)
      a2 = cf17(2)
      a3 = cf17(3)
      a4 = cf17(4)
      a5 = cf17(5)
      a6 = cf17(6)
      call phicalculator(l)
      pwf17 = partialwave(i1,i2,s,i,pwchoice,a1,a2,a3,a4,a5,a6)
!     ---------------- g17
      pwchoice = 'g17'
      l  = 4
      a1 = cg17(1)
      a2 = cg17(2)
      a3 = cg17(3)
      a4 = cg17(4)
      a5 = cg17(5)
      a6 = cg17(6)
      call phicalculator(l)
      pwg17 = partialwave(i1,i2,s,i,pwchoice,a1,a2,a3,a4,a5,a6)
!     t_l+ isospin 0
      f0(1) = pws01
      f0(2) = pwp03 
      f0(3) = pwd05
      f0(4) = pwf07
      f0(5) = xzero
!     t_l+ isospin 1
      f1(1) = pws11
      f1(2) = pwp13 
      f1(3) = pwd15
      f1(4) = pwf17
      f1(5) = xzero
!     t_l- isospin 0
      g0(1) = xzero
      g0(2) = pwp01
      g0(3) = pwd03
      g0(4) = pwf05
      g0(5) = pwg07
!     t_l- isospin 1
      g0(1) = xzero
      g1(2) = pwp11
      g1(3) = pwd13
      g1(4) = pwf15
      g1(5) = pwg17
      return
      end subroutine pwcalculator
!     -----------------------------------------------------------------------
      complex (8) function partialwave(i1,i2,s,index,pwchoice,
     c     a1,a2,a3,a4,a5,a6) result (output)
      use usefulnumbers
      use channelassignment
      use matrixstorage
      implicit none
      character (3) :: pwchoice
      integer :: i1,i2,index,a1,a2,a3,a4,a5,a6
      real (8) :: s
      complex (8) :: teoria
      call tpwcalculation(s,index,pwchoice,a1,a2,a3,a4,a5,a6)
      if (i1.eq.i2) then
         teoria = (scatmat(i1,i2) - xr)/2.d0/xi
      else
         teoria = scatmat(i1,i2)/2.d0/xi
      endif
      output = teoria
      return
      end function partialwave
!     -----------------------------------------------------------------------
      subroutine parametersuploader
      implicit none
      character (18) :: name
      logical :: upload
      logical, external :: params
      name = 'parameters.s01.inp'
      upload =  params(name,1,43)
!      if (upload) print*,'uploaded ',name
      name = 'parameters.p01.inp'
      upload =  params(name,44,85)   
!      if (upload) print*,'uploaded ',name
      name = 'parameters.p03.inp'
      upload = params(name,86,121)
!      if (upload) print*,'uploaded ',name
      name = 'parameters.d03.inp'
      upload = params(name,122,149)
!      if (upload) print*,'uploaded ',name
      name = 'parameters.d05.inp'
      upload = params(name,150,167)
!      if (upload) print*,'uploaded ',name
      name = 'parameters.f05.inp'
      upload = params(name,168,194)
!      if (upload) print*,'uploaded ',name
      name = 'parameters.f07.inp'
      upload = params(name,195,204)
!      if (upload) print*,'uploaded ',name
      name = 'parameters.g07.inp'
      upload = params(name,205,218)
!      if (upload) print*,'uploaded ',name
      name = 'parameters.s11.inp'
      upload = params(name,219,284)
!      if (upload) print*,'uploaded ',name
      name = 'parameters.p11.inp'
      upload = params(name,285,334)
!      if (upload) print*,'uploaded ',name
      name = 'parameters.p13.inp'
      upload = params(name,335,406)
!      if (upload) print*,'uploaded ',name
      name = 'parameters.d13.inp'
      upload = params(name,407,448)
!      if (upload) print*,'uploaded ',name
      name = 'parameters.d15.inp'
      upload = params(name,449,484)
!      if (upload) print*,'uploaded ',name
      name = 'parameters.f15.inp'
      upload = params(name,485,523)
!      if (upload) print*,'uploaded ',name
      name = 'parameters.f17.inp'
      upload = params(name,524,547)
!      if (upload) print*,'uploaded ',name
      name = 'parameters.g17.inp'
      upload = params(name,548,558)   
!      if (upload) print*,'uploaded ',name
      return
      end subroutine parametersuploader
!     -----------------------------------------------------------------------
      logical function params(name,i1,i2) result (output)
      use coeficientes
      implicit none
      logical :: uptemp
      character (18) :: name
      integer :: i1,i2,ivoid,i,j
      real (8) :: rvoid,rvoid2
      uptemp = .false.
      open(1,file=name,status='old')
      read (1,*)ivoid
      read (1,*)
      do j=i1,i2
         read (1,*)ivoid,rvoid
         coeficiente(j) = rvoid
      enddo
      uptemp = .true.
      output = uptemp
      close (1)
      return
      end function params
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      complex (8) function a16function(t1,t2,t3,t4,t5,t6,
     c     e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
     c     result (output)
      implicit none
      complex (8), dimension (6) :: tamp
      complex (8), dimension (15) :: epsilons
      complex (8) :: t1,t2,t3,t4,t5,t6
      complex (8) :: e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,
     c     e45,e46,e56
      complex (8), external :: bfunction
      tamp(1) = t1
      tamp(2) = t2
      tamp(3) = t3
      tamp(4) = t4
      tamp(5) = t5
      tamp(6) = t6
      epsilons(1)  = e12
      epsilons(2)  = e13
      epsilons(3)  = e14
      epsilons(4)  = e15
      epsilons(5)  = e16
      epsilons(6)  = e23
      epsilons(7)  = e24
      epsilons(8)  = e25
      epsilons(9)  = e26
      epsilons(10) = e34
      epsilons(11) = e35
      epsilons(12) = e36
      epsilons(13) = e45
      epsilons(14) = e46
      epsilons(15) = e56
      output  = bfunction(tamp,epsilons)
      return
      end function a16function
!     -----------------------------------------------------------------------
      complex (8) function a26function(t1,t2,t3,t4,t5,t6,
     c     e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
     c     result (output)
      implicit none
      complex (8), dimension (6) :: tamp
      complex (8), dimension (15) :: epsilons
      complex (8) :: t1,t2,t3,t4,t5,t6
      complex (8) :: e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,
     c     e45,e46,e56
      complex (8), external :: bfunction
      tamp(1) = t2 !t1 
      tamp(2) = t1 !t2
      tamp(3) = t3
      tamp(4) = t4
      tamp(5) = t5
      tamp(6) = t6
      epsilons(1)  = e12 !e12
      epsilons(2)  = e23 !e13
      epsilons(3)  = e24 !e14
      epsilons(4)  = e25 !e15
      epsilons(5)  = e26 !e16
      epsilons(6)  = e13 !e23
      epsilons(7)  = e14 !e24
      epsilons(8)  = e15 !e25
      epsilons(9)  = e16 !e26
      epsilons(10) = e34
      epsilons(11) = e35
      epsilons(12) = e36
      epsilons(13) = e45
      epsilons(14) = e46
      epsilons(15) = e56
      output  = bfunction(tamp,epsilons)
      return
      end function a26function
!     -----------------------------------------------------------------------
      complex (8) function a36function(t1,t2,t3,t4,t5,t6,
     c     e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
     c     result (output)
      implicit none
      complex (8), dimension (6) :: tamp
      complex (8), dimension (15) :: epsilons
      complex (8) :: t1,t2,t3,t4,t5,t6
      complex (8) :: e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,
     c     e45,e46,e56
      complex (8), external :: bfunction
      tamp(1) = t3 !t1
      tamp(2) = t2
      tamp(3) = t1 !t3
      tamp(4) = t4
      tamp(5) = t5
      tamp(6) = t6
      epsilons(1)  = e23 !e12
      epsilons(2)  = e13 !e13
      epsilons(3)  = e34 !e14
      epsilons(4)  = e35 !e15
      epsilons(5)  = e36 !e16
      epsilons(6)  = e12 !e23
      epsilons(7)  = e24
      epsilons(8)  = e25
      epsilons(9)  = e26
      epsilons(10) = e14 !e34
      epsilons(11) = e15 !e35
      epsilons(12) = e16 !e36
      epsilons(13) = e45
      epsilons(14) = e46
      epsilons(15) = e56
      output  = bfunction(tamp,epsilons)
      return
      end function a36function
!     -----------------------------------------------------------------------
      complex (8) function a46function(t1,t2,t3,t4,t5,t6,
     c     e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
     c     result (output)
      implicit none
      complex (8), dimension (6) :: tamp
      complex (8), dimension (15) :: epsilons
      complex (8) :: t1,t2,t3,t4,t5,t6
      complex (8) :: e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,
     c     e45,e46,e56
      complex (8), external :: bfunction
      tamp(1) = t4 !t1
      tamp(2) = t2
      tamp(3) = t3
      tamp(4) = t1 !t4
      tamp(5) = t5
      tamp(6) = t6
      epsilons(1)  = e24 !e12
      epsilons(2)  = e34 !e13
      epsilons(3)  = e14 !e14
      epsilons(4)  = e45 !e15
      epsilons(5)  = e46 !e16
      epsilons(6)  = e23
      epsilons(7)  = e12 !e24
      epsilons(8)  = e25
      epsilons(9)  = e26
      epsilons(10) = e13 !e34
      epsilons(11) = e35
      epsilons(12) = e36
      epsilons(13) = e15 !e45
      epsilons(14) = e16 !e46
      epsilons(15) = e56
      output  = bfunction(tamp,epsilons)
      return
      end function a46function
!     -----------------------------------------------------------------------
      complex (8) function a56function(t1,t2,t3,t4,t5,t6,
     c     e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
     c     result (output)
      implicit none
      complex (8), dimension (6) :: tamp
      complex (8), dimension (15) :: epsilons
      complex (8) :: t1,t2,t3,t4,t5,t6
      complex (8) :: e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,
     c     e45,e46,e56
      complex (8), external :: bfunction
      tamp(1) = t5 !t1
      tamp(2) = t2
      tamp(3) = t3
      tamp(4) = t4
      tamp(5) = t1 !t5
      tamp(6) = t6
      epsilons(1)  = e25 !e12
      epsilons(2)  = e35 !e13
      epsilons(3)  = e45 !e14
      epsilons(4)  = e15 !e15
      epsilons(5)  = e56 !e16
      epsilons(6)  = e23
      epsilons(7)  = e24
      epsilons(8)  = e12 !e25
      epsilons(9)  = e26
      epsilons(10) = e34
      epsilons(11) = e13 !e35
      epsilons(12) = e36
      epsilons(13) = e14 !e45
      epsilons(14) = e46
      epsilons(15) = e16 !e56
      output  = bfunction(tamp,epsilons)
      return
      end function a56function
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      complex (8) function a12function(t1,t2,t3,t4,t5,t6,
     c     e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
     c     result (output)
      implicit none
      complex (8), dimension (6) :: tamp
      complex (8), dimension (15) :: epsilons
      complex (8) :: t1,t2,t3,t4,t5,t6
      complex (8) :: e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,
     c     e45,e46,e56
      complex (8), external :: bfunction
      tamp(1) = t1
      tamp(2) = t6 !t2
      tamp(3) = t3
      tamp(4) = t4
      tamp(5) = t5
      tamp(6) = t2 !t6
      epsilons(1)  = e16 !e12
      epsilons(2)  = e13
      epsilons(3)  = e14
      epsilons(4)  = e15
      epsilons(5)  = e12 !e16
      epsilons(6)  = e36 !e23
      epsilons(7)  = e46 !e24
      epsilons(8)  = e56 !e25
      epsilons(9)  = e26 !e26
      epsilons(10) = e34
      epsilons(11) = e35
      epsilons(12) = e23 !e36
      epsilons(13) = e45
      epsilons(14) = e24 !e46
      epsilons(15) = e25 !e56
      output  = bfunction(tamp,epsilons)
      return
      end function a12function
!     -----------------------------------------------------------------------
      complex (8) function a13function(t1,t2,t3,t4,t5,t6,
     c     e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
     c     result (output)
      implicit none
      complex (8), dimension (6) :: tamp
      complex (8), dimension (15) :: epsilons
      complex (8) :: t1,t2,t3,t4,t5,t6
      complex (8) :: e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,
     c     e45,e46,e56
      complex (8), external :: bfunction
      tamp(1) = t1
      tamp(2) = t2
      tamp(3) = t6 !t3
      tamp(4) = t4
      tamp(5) = t5
      tamp(6) = t3 !t6
      epsilons(1)  = e12
      epsilons(2)  = e16 !e13
      epsilons(3)  = e14
      epsilons(4)  = e15
      epsilons(5)  = e13 !e16
      epsilons(6)  = e26 !e23
      epsilons(7)  = e24
      epsilons(8)  = e25
      epsilons(9)  = e23 ! e26
      epsilons(10) = e46 !e34
      epsilons(11) = e56 !e35
      epsilons(12) = e36 !e36
      epsilons(13) = e45
      epsilons(14) = e34 !e46
      epsilons(15) = e35 !e56
      output  = bfunction(tamp,epsilons)
      return
      end function a13function
!     -----------------------------------------------------------------------
      complex (8) function a14function(t1,t2,t3,t4,t5,t6,
     c     e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
     c     result (output)
      implicit none
      complex (8), dimension (6) :: tamp
      complex (8), dimension (15) :: epsilons
      complex (8) :: t1,t2,t3,t4,t5,t6
      complex (8) :: e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,
     c     e45,e46,e56
      complex (8), external :: bfunction
      tamp(1) = t1
      tamp(2) = t2
      tamp(3) = t3
      tamp(4) = t6 !t4
      tamp(5) = t5
      tamp(6) = t4 !t6
      epsilons(1)  = e12
      epsilons(2)  = e13
      epsilons(3)  = e16 !e14
      epsilons(4)  = e15
      epsilons(5)  = e14 !e16
      epsilons(6)  = e23
      epsilons(7)  = e26 !e24
      epsilons(8)  = e25
      epsilons(9)  = e24 !e26
      epsilons(10) = e36 !e34
      epsilons(11) = e35
      epsilons(12) = e34 !e36 
      epsilons(13) = e56 !e45
      epsilons(14) = e46 !e46
      epsilons(15) = e45 !e56
      output  = bfunction(tamp,epsilons)
      return
      end function a14function
!     -----------------------------------------------------------------------
      complex (8) function a15function(t1,t2,t3,t4,t5,t6,
     c     e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
     c     result (output)
      implicit none
      complex (8), dimension (6) :: tamp
      complex (8), dimension (15) :: epsilons
      complex (8) :: t1,t2,t3,t4,t5,t6
      complex (8) :: e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,
     c     e45,e46,e56
      complex (8), external :: bfunction
      tamp(1) = t1
      tamp(2) = t2
      tamp(3) = t3
      tamp(4) = t4
      tamp(5) = t6 !t5
      tamp(6) = t5 !t6
      epsilons(1)  = e12
      epsilons(2)  = e13
      epsilons(3)  = e14
      epsilons(4)  = e16 !e15
      epsilons(5)  = e15 !e16
      epsilons(6)  = e23
      epsilons(7)  = e24
      epsilons(8)  = e26 !e25
      epsilons(9)  = e25 !e26
      epsilons(10) = e34
      epsilons(11) = e36 !e35
      epsilons(12) = e35 !e36
      epsilons(13) = e46 !e45
      epsilons(14) = e45 !e46
      epsilons(15) = e56 !e56
      output  = bfunction(tamp,epsilons)
      return
      end function a15function
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      complex (8) function a23function(t1,t2,t3,t4,t5,t6,
     c     e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
     c     result (output)
      implicit none
      complex (8), dimension (6) :: tamp
      complex (8), dimension (15) :: epsilons
      complex (8) :: t1,t2,t3,t4,t5,t6
      complex (8) :: e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,
     c     e45,e46,e56
      complex (8), external :: bfunction
      tamp(1) = t2 !t1
      tamp(2) = t1 !t2
      tamp(3) = t6 !t3
      tamp(4) = t4
      tamp(5) = t5
      tamp(6) = t3 !t6
      epsilons(1)  = e12 !e12
      epsilons(2)  = e26 !e13
      epsilons(3)  = e24 !e14
      epsilons(4)  = e25 !e15
      epsilons(5)  = e23 !e16
      epsilons(6)  = e16 !e23
      epsilons(7)  = e14 !e24
      epsilons(8)  = e15 !e25
      epsilons(9)  = e13 !e26
      epsilons(10) = e46 !e34
      epsilons(11) = e56 !e35
      epsilons(12) = e36 !e36
      epsilons(13) = e45
      epsilons(14) = e34 !e46
      epsilons(15) = e35 !e56
      output  = bfunction(tamp,epsilons)
      return
      end function a23function
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      complex (8) function a34function(t1,t2,t3,t4,t5,t6,
     c     e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
     c     result (output)
      implicit none
      complex (8), dimension (6) :: tamp
      complex (8), dimension (15) :: epsilons
      complex (8) :: t1,t2,t3,t4,t5,t6
      complex (8) :: e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,
     c     e45,e46,e56
      complex (8), external :: bfunction
      tamp(1) = t3 !t1
      tamp(2) = t2
      tamp(3) = t1 !t3
      tamp(4) = t6 !t4
      tamp(5) = t5
      tamp(6) = t4 !t6
      epsilons(1)  = e23 !e12
      epsilons(2)  = e13 !e13
      epsilons(3)  = e36 !e14
      epsilons(4)  = e35 !e15
      epsilons(5)  = e34 !e16
      epsilons(6)  = e12 !e23
      epsilons(7)  = e26 !e24
      epsilons(8)  = e25
      epsilons(9)  = e24 !e26
      epsilons(10) = e16 !e34
      epsilons(11) = e15 !e35
      epsilons(12) = e14 !e36
      epsilons(13) = e56 !e45
      epsilons(14) = e46 !e46
      epsilons(15) = e45 !e56
      output  = bfunction(tamp,epsilons)
      return
      end function a34function
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      complex (8) function a24function(t1,t2,t3,t4,t5,t6,
     c     e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
     c     result (output)
      implicit none
      complex (8), dimension (6) :: tamp
      complex (8), dimension (15) :: epsilons
      complex (8) :: t1,t2,t3,t4,t5,t6
      complex (8) :: e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,
     c     e45,e46,e56
      complex (8), external :: bfunction
      tamp(1) = t2 !t1
      tamp(2) = t1 !t2
      tamp(3) = t3
      tamp(4) = t6 !t4
      tamp(5) = t5
      tamp(6) = t4 !t6
      epsilons(1)  = e12 !e12
      epsilons(2)  = e23 !e13
      epsilons(3)  = e26 !e14
      epsilons(4)  = e25 !e15
      epsilons(5)  = e24 !e16
      epsilons(6)  = e13 !e23
      epsilons(7)  = e16 !e24
      epsilons(8)  = e15 !e25
      epsilons(9)  = e14 !e26
      epsilons(10) = e36 !e34
      epsilons(11) = e35
      epsilons(12) = e34 !e36
      epsilons(13) = e56 !e45
      epsilons(14) = e46 !e46
      epsilons(15) = e45 !e56
      output  = bfunction(tamp,epsilons)
      return
      end function a24function
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      complex (8) function a25function(t1,t2,t3,t4,t5,t6,
     c     e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
     c     result (output)
      implicit none
      complex (8), dimension (6) :: tamp
      complex (8), dimension (15) :: epsilons
      complex (8) :: t1,t2,t3,t4,t5,t6
      complex (8) :: e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,
     c     e45,e46,e56
      complex (8), external :: bfunction
      tamp(1) = t2 !t1 
      tamp(2) = t1 !t2
      tamp(3) = t3
      tamp(4) = t4
      tamp(5) = t6 !t5
      tamp(6) = t5 !t6
      epsilons(1)  = e12 !e12
      epsilons(2)  = e23 !e13
      epsilons(3)  = e24 !e14
      epsilons(4)  = e26 !e15
      epsilons(5)  = e25 !e16
      epsilons(6)  = e13 !e23
      epsilons(7)  = e14 !e24
      epsilons(8)  = e16 !e25
      epsilons(9)  = e15 !e26
      epsilons(10) = e34
      epsilons(11) = e36 !e35
      epsilons(12) = e35 !e36
      epsilons(13) = e46 !e45
      epsilons(14) = e45 !e46
      epsilons(15) = e56 !e56
      output  = bfunction(tamp,epsilons)
      return
      end function a25function
!     -----------------------------------------------------------------------
      complex (8) function a35function(t1,t2,t3,t4,t5,t6,
     c     e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
     c     result (output)
      implicit none
      complex (8), dimension (6) :: tamp
      complex (8), dimension (15) :: epsilons
      complex (8) :: t1,t2,t3,t4,t5,t6
      complex (8) :: e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,
     c     e45,e46,e56
      complex (8), external :: bfunction
      tamp(1) = t3 !t1
      tamp(2) = t2
      tamp(3) = t1 !t3
      tamp(4) = t4
      tamp(5) = t6 !t5
      tamp(6) = t5 !t6
      epsilons(1)  = e23 !e12
      epsilons(2)  = e13 !e13
      epsilons(3)  = e34 !e14
      epsilons(4)  = e36 !e15
      epsilons(5)  = e35 !e16
      epsilons(6)  = e12 !e23
      epsilons(7)  = e24
      epsilons(8)  = e26 !e25
      epsilons(9)  = e25 !e26
      epsilons(10) = e14 !e34
      epsilons(11) = e16 !e35
      epsilons(12) = e15 !e36
      epsilons(13) = e46 !e45
      epsilons(14) = e45 !e46
      epsilons(15) = e56 !e56
      output  = bfunction(tamp,epsilons)
      return
      end function a35function
!     -----------------------------------------------------------------------
      complex (8) function a45function(t1,t2,t3,t4,t5,t6,
     c     e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
     c     result (output)
      implicit none
      complex (8), dimension (6) :: tamp
      complex (8), dimension (15) :: epsilons
      complex (8) :: t1,t2,t3,t4,t5,t6
      complex (8) :: e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,
     c     e45,e46,e56
      complex (8), external :: bfunction
      tamp(1) = t4 !t1
      tamp(2) = t2
      tamp(3) = t3
      tamp(4) = t1 !t4
      tamp(5) = t6 !t5
      tamp(6) = t5 !t6
      epsilons(1)  = e24 !e12
      epsilons(2)  = e34 !e13
      epsilons(3)  = e14 !e14
      epsilons(4)  = e46 !e15
      epsilons(5)  = e45 !e16
      epsilons(6)  = e23
      epsilons(7)  = e12 !e24
      epsilons(8)  = e26 !e25
      epsilons(9)  = e25 !e26
      epsilons(10) = e13 !e34
      epsilons(11) = e36 !e35
      epsilons(12) = e35 !e36
      epsilons(13) = e16 !e45
      epsilons(14) = e15 !e46
      epsilons(15) = e56 !e56
      output  = bfunction(tamp,epsilons)
      return
      end function a45function
!     -----------------------------------------------------------------------
      complex (8) function afunction(t1,t2,t3,t4,t5,t6,
     c     e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
      implicit none
      complex (8), parameter :: xr = dcmplx(1.d0,0.d0)
      complex (8), parameter :: xi = dcmplx(0.d0,1.d0)
      real (8), parameter :: units1 = 1000.d0
      real (8), parameter :: units2 = units1*units1
      real (8), parameter :: units3 = units1*units2
      real (8), parameter :: units4 = units1*units3
      real (8), parameter :: units5 = units1*units4
      real (8), parameter :: units6 = units1*units5
      complex (8) :: t1,t2,t3,t4,t5,t6,salida
      complex (8) :: e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,
     c     e45,e46,e56
      complex (8) :: t23456,t2456,t3456,t2356,t2345,t2346,t234,t235,
     c     t345,t245,t346,t246,t256,t356,t236,t456
      complex (8) :: f23456,f2456,f3456,f2356,f2345,f2346
      t23456 = t2*t3*t4*t5*t6
      t2456 = t2*t4*t5*t6
      t3456 = t3*t4*t5*t6
      t2356 = t2*t3*t5*t6
      t2345 = t2*t3*t4*t5
      t2346 = t2*t3*t4*t6
      t234  = t2*t3*t4
      t235  = t2*t3*t5
      t345  = t3*t4*t5
      t245  = t2*t4*t5 
      t346  = t3*t4*t6 
      t246  = t2*t4*t6
      t256  = t2*t5*t6
      t356  = t3*t5*t6
      t236  = t2*t3*t6
      t456  = t4*t5*t6
      f2345 = e25**2*e34**2 - 2.d0*e24*e25*e34*e35 + 
     -     e24**2*e35**2 - 2.d0*e23*e25*e34*e45 - 
     -     2.d0*e23*e24*e35*e45 + e23**2*e45**2
      f23456 = e26**2*e34*e35*e45 - e25*e26*e34*e36*e45 - 
     -     e24*e26*e35*e36*e45 + e24*e25*e36**2*e45 + 
     -     e23*e26*e36*e45**2 - e25*e26*e34*e35*e46 + 
     -     e24*e26*e35**2*e46 + e25**2*e34*e36*e46 - 
     -     e24*e25*e35*e36*e46 - e23*e26*e35*e45*e46 - 
     -     e23*e25*e36*e45*e46 + e23*e25*e35*e46**2 + 
     -     e25*e26*e34**2*e56 - e24*e26*e34*e35*e56 - 
     -     e24*e25*e34*e36*e56 + e24**2*e35*e36*e56 - 
     -     e23*e26*e34*e45*e56 -  e23*e24*e36*e45*e56 - 
     -     e23*e25*e34*e46*e56 -  e23*e24*e35*e46*e56 + 
     -     e23**2*e45*e46*e56 +  e23*e24*e34*e56**2
      f2346 = e26**2*e34**2- 2.d0*e24*e26*e34*e36 + 
     -     e24**2*e36**2 - 2.d0*e23*e26*e34*e46 - 
     -     2.d0*e23*e24*e36*e46 + e23**2*e46**2
      f2356 = e26**2*e35**2 - 2.d0*e25*e26*e35*e36 + e25**2*e36**2 - 
     -     2.d0*e23*e26*e35*e56 - 2.d0*e23*e25*e36*e56 + 
     -     e23**2*e56**2
      f2456 = e26**2*e45**2 - 2.d0*e25*e26*e45*e46 + e25**2*e46**2 - 
     -     2.d0*e24*e26*e45*e56 - 2.d0*e24*e25*e46*e56 + 
     -     e24**2*e56**2
      f3456 = e36**2*e45**2 - 2.d0*e35*e36*e45*e46 +  e35**2*e46**2 - 
     -     2.d0*e34*e36*e45*e56 - 2.d0*e34*e35*e46*e56 + 
     -     e34**2*e56**2
      salida =     xr + ( e23**2*t2*t3 + e24**2*t2*t4 + 
     -     e34**2*t3*t4 + e25**2*t2*t5 + e35**2*t3*t5 + 
     -     e45**2*t4*t5 + e26**2*t2*t6 + e36**2*t3*t6 + 
     -     e46**2*t4*t6 + e56**2*t5*t6 )/units2 +
     -     ( f2345*t2345 + t2346*f2346 + f2356*t2356 +
     -     t2456*f2456 + f3456*t3456 )/units4 +
     -     2.d0*xi*( e23*e24*e34*t234 + e34*e36*e46*t346 + 
     -     e35*e36*e56*t356 + e45*e46*e56*t456 + e23*e25*e35*t235 +
     -     e23*e26*e36*t236 + e24*e26*e46*t246 + e25*e26*e56*t256 + 
     -     e34*e35*e45*t345 + e24*e25*e45*t245 )/units3 +
     -     2.d0*xi*f23456*t23456/units5
      afunction = t1*salida
      return
      end function afunction
!     -----------------------------------------------------------------------
      complex (8) function bfunction(tamp,epsilons)
      implicit none
      complex (8), parameter :: xr = dcmplx(1.d0,0.d0)
      complex (8), parameter :: xi = dcmplx(0.d0,1.d0)
      real (8), parameter :: units1 = 1000.d0
      real (8), parameter :: units2 = units1*units1
      real (8), parameter :: units3 = units1*units2
      real (8), parameter :: units4 = units1*units3
      real (8), parameter :: units5 = units1*units4
      real (8), parameter :: units6 = units1*units5
      complex (8), dimension (6) :: tamp
      complex (8), dimension (15) :: epsilons
      complex (8) :: t1,t2,t3,t4,t5,t6,salida
      complex (8) :: e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,
     c     e45,e46,e56
      complex (8) :: t2345,t345,t245,t235,t234
      complex (8) :: f2345,f345,f245,f235,f234
      t1 = tamp(1)
      t2 = tamp(2)
      t3 = tamp(3)
      t4 = tamp(4)
      t5 = tamp(5)
      t6 = tamp(6)
      e12 = epsilons(1) 
      e13 = epsilons(2) 
      e14 = epsilons(3)
      e15 = epsilons(4)
      e16 = epsilons(5) 
      e23 = epsilons(6) 
      e24 = epsilons(7) 
      e25 = epsilons(8) 
      e26 = epsilons(9) 
      e34 = epsilons(10) 
      e35 = epsilons(11)
      e36 = epsilons(12) 
      e45 = epsilons(13)
      e46 = epsilons(14) 
      e56 = epsilons(15) 
      t2345 = t2*t3*t4*t5
      t345 = t3*t4*t5
      t245 = t2*t4*t5
      t235 = t2*t3*t5
      t234 = t2*t3*t4
      f2345 = e16*e25**2*e34**2 - e15*e25*e26*e34**2 - 
     -       2.d0*e16*e24*e25*e34*e35 + e15*e24*e26*e34*e35 + 
     -       e14*e25*e26*e34*e35 + e16*e24**2*e35**2 - 
     -       e14*e24*e26*e35**2 + e15*e24*e25*e34*e36 - 
     -       e14*e25**2*e34*e36 - e15*e24**2*e35*e36 + 
     -       e14*e24*e25*e35*e36 - 2.d0*e16*e23*e25*e34*e45 + 
     -       e15*e23*e26*e34*e45 + e13*e25*e26*e34*e45 - 
     -       2.d0*e16*e23*e24*e35*e45 + e14*e23*e26*e35*e45 + 
     -       e13*e24*e26*e35*e45 - 2.d0*e12*e26*e34*e35*e45 + 
     -       e15*e23*e24*e36*e45 + e14*e23*e25*e36*e45 - 
     -       2.d0*e13*e24*e25*e36*e45 + e12*e25*e34*e36*e45 + 
     -       e12*e24*e35*e36*e45 + e16*e23**2*e45**2 - 
     -       e13*e23*e26*e45**2 - e12*e23*e36*e45**2 + 
     -       e15*e23*e25*e34*e46 - e13*e25**2*e34*e46 + 
     -       e15*e23*e24*e35*e46 - 2.d0*e14*e23*e25*e35*e46 + 
     -       e13*e24*e25*e35*e46 + e12*e25*e34*e35*e46 - 
     -       e12*e24*e35**2*e46 - e15*e23**2*e45*e46 + 
     -       e13*e23*e25*e45*e46 + e12*e23*e35*e45*e46 - 
     -       2.d0*e15*e23*e24*e34*e56 + e14*e23*e25*e34*e56 + 
     -       e13*e24*e25*e34*e56 - e12*e25*e34**2*e56 + 
     -       e14*e23*e24*e35*e56 - e13*e24**2*e35*e56 + 
     -       e12*e24*e34*e35*e56 - e14*e23**2*e45*e56 + 
     -       e13*e23*e24*e45*e56 + e12*e23*e34*e45*e56
       salida = xi*e16 -
     -     ( e12*e26*t2 + e14*e46*t4 + e15*e56*t5 + e13*e36*t3 )/units1
     -     + (- 2*e16*e23*e24*e34*t234 + e14*e23*e26*e34*t234 + 
     -     e13*e24*e26*e34*t234 - e12*e26*e34**2*t234 + 
     -     e14*e23*e24*e36*t234 - e13*e24**2*e36*t234 + 
     -     e12*e24*e34*e36*t234 - e14*e23**2*e46*t234 + 
     -     e13*e23*e24*e46*t234 + e12*e23*e34*e46*t234 -
     -     2*e16*e23*e25*e35*t235 + e15*e23*e26*e35*t235 + 
     -     e13*e25*e26*e35*t235 - e12*e26*e35**2*t235 + 
     -     e15*e23*e25*e36*t235 - e13*e25**2*e36*t235 + 
     -     e12*e25*e35*e36*t235 - e15*e23**2*e56*t235 + 
     -     e13*e23*e25*e56*t235 + e12*e23*e35*e56*t235 - 
     -     2*e16*e24*e25*e45*t245 + e15*e24*e26*e45*t245 + 
     -     e14*e25*e26*e45*t245 - e12*e26*e45**2*t245 + 
     -     e15*e24*e25*e46*t245 - e14*e25**2*e46*t245 + 
     -     e12*e25*e45*e46*t245 - e15*e24**2*e56*t245 + 
     -     e14*e24*e25*e56*t245 + e12*e24*e45*e56*t245  + 
     -     e15*e34*e36*e45*t345 + e14*e35*e36*e45*t345 - 
     -     e13*e36*e45**2*t345 + e15*e34*e35*e46*t345 - 
     -     e14*e35**2*e46*t345 + e13*e35*e45*e46*t345 - 
     -     e15*e34**2*e56*t345 + e14*e34*e35*e56*t345 + 
     -     e13*e34*e45*e56*t345 - 2*e16*e34*e35*e45*t345 )/units3
     -     + xi*( e16*e23**2*t2*t3 - e13*e23*e26*t2*t3 - 
     -     e12*e23*e36*t2*t3 + e16*e24**2*t2*t4 - e14*e24*e26*t2*t4 - 
     -     e12*e24*e46*t2*t4 + e16*e34**2*t3*t4 - e14*e34*e36*t3*t4 - 
     -     e13*e34*e46*t3*t4 + e16*e25**2*t2*t5 - e15*e25*e26*t2*t5 - 
     -     e12*e25*e56*t2*t5 + e16*e35**2*t3*t5 - e15*e35*e36*t3*t5 - 
     -     e13*e35*e56*t3*t5 + e16*e45**2*t4*t5 - e15*e45*e46*t4*t5 - 
     -     e14*e45*e56*t4*t5 )/units2 + xi*t2345*f2345/units4
      bfunction = t1*t6*salida
      return
      end function bfunction
!     -----------------------------------------------------------------------
      complex (8) function determinante(t1,t2,t3,t4,t5,t6,
     c     e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,
     c     e36,e45,e46,e56) result (output)
      implicit none
      complex (8), parameter :: xr = dcmplx(1.d0,0.d0)
      complex (8), parameter :: xi = dcmplx(0.d0,1.d0)
      real (8), parameter :: units1 = 1000.d0
      real (8), parameter :: units2 = units1*units1
      real (8), parameter :: units3 = units1*units2
      real (8), parameter :: units4 = units1*units3
      real (8), parameter :: units5 = units1*units4
      real (8), parameter :: units6 = units1*units5
      complex (8) :: t1,t2,t3,t4,t5,t6,salida
      complex (8) :: e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,
     c     e45,e46,e56
      complex (8) :: t123456,t23456,t13456,t12456,t12356,t12345,
     c     t12346,t3456,t2456,t1456,t1356,t1346,t1256,t1246,t1236,t1345,
     c     t1245,t2345,t2356,t2346,t1235,t1234
      complex (8) :: f123456,f23456,f13456,f12456,f12356,f12345,f12346,
     c     f3456,f2456,f1456,f1356,f1346,f1256,f1246,f1236,f1345,f1245,
     c     f2345,f2356,f2346,f1235,f1234
      complex (8) :: f01,f04,f021,f022,f023,f024,f025,f026,f027,f028
      t123456 = t1*t2*t3*t4*t5*t6
      t23456 = t2*t3*t4*t5*t6
      t13456 = t1*t3*t4*t5*t6
      t12456 = t1*t2*t4*t5*t6
      t12356 = t1*t2*t3*t5*t6
      t12345 = t1*t2*t3*t4*t5
      t12346 = t1*t2*t3*t4*t6 
      t3456 = t3*t4*t5*t6 
      t2456 = t2*t4*t5*t6 
      t1456 = t1*t4*t5*t6 
      t1356 = t1*t3*t5*t6 
      t1346 = t1*t3*t4*t6 
      t1256 = t1*t2*t5*t6 
      t1246 = t1*t2*t4*t6 
      t1236 = t1*t2*t3*t6 
      t1345 = t1*t3*t4*t5
      t1245 = t1*t2*t4*t5
      t2345 = t2*t3*t4*t5
      t2356 = t2*t3*t5*t6 
      t2346 = t2*t3*t4*t6 
      t1235 = t1*t2*t3*t5
      t1234 = t1*t2*t3*t4
      f1234 = e14**2*e23**2 - 2.d0*e13*e14*e23*e24 + 
     -     e13**2*e24**2 - 2.d0*e12*e14*e23*e34 - 
     -     2.d0*e12*e13*e24*e34 + e12**2*e34**2
      f1235 = e15**2*e23**2 - 2.d0*e13*e15*e23*e25 + 
     -     e13**2*e25**2 - 2.d0*e12*e15*e23*e35 - 
     -     2.d0*e12*e13*e25*e35 + e12**2*e35**2
      f1236 = e16**2*e23**2 - 2.d0*e13*e16*e23*e26 + e13**2*e26**2 -
     -     2.d0*e12*e16*e23*e36 - 2.d0*e12*e13*e26*e36 + e12**2*e36**2
      f1245 = e15**2*e24**2 - 
     -     2.d0*e14*e15*e24*e25 + e14**2*e25**2 - 2.d0*e12*e15*e24*e45 - 
     -     2.d0*e12*e14*e25*e45 + e12**2*e45**2
      f1246 = e16**2*e24**2 - 2.d0*e14*e16*e24*e26 + 
     -     e14**2*e26**2 - 2.d0*e12*e16*e24*e46 - 
     -     2.d0*e12*e14*e26*e46 + e12**2*e46**2
      f1256 = e16**2*e25**2 - 2.d0*e15*e16*e25*e26 + 
     -     e15**2*e26**2 - 2.d0*e12*e16*e25*e56 - 2.d0*e12*e15*e26*e56 + 
     -     e12**2*e56**2
      f1345 = e15**2*e34**2 - 2.d0*e14*e15*e34*e35 + 
     -     e14**2*e35**2 - 2.d0*e13*e15*e34*e45 - 
     -     2.d0*e13*e14*e35*e45 + e13**2*e45**2
      f1346 = e16**2*e34**2 - 2.d0*e14*e16*e34*e36 + 
     -     e14**2*e36**2 - 2.d0*e13*e16*e34*e46 - 
     -     2.d0*e13*e14*e36*e46 + e13**2*e46**2
      f1356 = e16**2*e35**2 - 2.d0*e15*e16*e35*e36 + 
     -     e15**2*e36**2 - 2.d0*e13*e16*e35*e56 - 
     -     2.d0*e13*e15*e36*e56 + e13**2*e56**2
      f1456 = e16**2*e45**2 - 2.d0*e15*e16*e45*e46 + 
     -     e15**2*e46**2 - 2.d0*e14*e16*e45*e56 - 
     -     2.d0*e14*e15*e46*e56 + e14**2*e56**2
      f2345 = e25**2*e34**2 - 2.d0*e24*e25*e34*e35 + 
     -     e24**2*e35**2 - 2.d0*e23*e25*e34*e45 - 
     -     2.d0*e23*e24*e35*e45 + e23**2*e45**2
      f3456 = e36**2*e45**2 - 2.d0*e35*e36*e45*e46 + 
     -     e35**2*e46**2 - 2.d0*e34*e36*e45*e56 - 
     -     2.d0*e34*e35*e46*e56 + e34**2*e56**2 
      f2346 = e26**2*e34**2 - 2.d0*e24*e26*e34*e36 + 
     -     e24**2*e36**2 - 2.d0*e23*e26*e34*e46 - 
     -     2.d0*e23*e24*e36*e46 + e23**2*e46**2 
      f2356 = e26**2*e35**2 - 2.d0*e25*e26*e35*e36 + 
     -     e25**2*e36**2 - 2.d0*e23*e26*e35*e56 - 
     -     2.d0*e23*e25*e36*e56 + e23**2*e56**2 
      f2456 = e26**2*e45**2 - 2.d0*e25*e26*e45*e46 + 
     -     e25**2*e46**2 - 2.d0*e24*e26*e45*e56 - 
     -     2.d0*e24*e25*e46*e56 + e24**2*e56**2
      f23456 = (e26**2*e34*e35*e45 - e25*e26*e34*e36*e45- 
     -     e24*e26*e35*e36*e45 + e24*e25*e36**2*e45 + 
     -     e23*e26*e36*e45**2 - e25*e26*e34*e35*e46 + 
     -     e24*e26*e35**2*e46 + e25**2*e34*e36*e46 - 
     -     e24*e25*e35*e36*e46 - e23*e26*e35*e45*e46 - 
     -     e23*e25*e36*e45*e46 + e23*e25*e35*e46**2 + 
     -     e25*e26*e34**2*e56 - e24*e26*e34*e35*e56 - 
     -     e24*e25*e34*e36*e56 + e24**2*e35*e36*e56 - 
     -     e23*e26*e34*e45*e56 - e23*e24*e36*e45*e56 - 
     -     e23*e25*e34*e46*e56 - e23*e24*e35*e46*e56 + 
     -     e23**2*e45*e46*e56 + e23*e24*e34*e56**2 )/units5
      f13456 = (e16**2*e34*e35*e45 - 
     -     e15*e16*e34*e36*e45 - e14*e16*e35*e36*e45 + 
     -     e14*e15*e36**2*e45 + e13*e16*e36*e45**2 - 
     -     e15*e16*e34*e35*e46 + e14*e16*e35**2*e46 + 
     -     e15**2*e34*e36*e46 - e14*e15*e35*e36*e46 - 
     -     e13*e16*e35*e45*e46 - e13*e15*e36*e45*e46 + 
     -     e13*e15*e35*e46**2 + e15*e16*e34**2*e56 - 
     -     e14*e16*e34*e35*e56 - e14*e15*e34*e36*e56 + 
     -     e14**2*e35*e36*e56 - e13*e16*e34*e45*e56 - 
     -     e13*e14*e36*e45*e56 - e13*e15*e34*e46*e56 - 
     -     e13*e14*e35*e46*e56 + e13**2*e45*e46*e56 + 
     -     e13*e14*e34*e56**2)/units5
      f12456 = (e16**2*e24*e25*e45 - e15*e16*e24*e26*e45 - 
     -     e14*e16*e25*e26*e45 + e14*e15*e26**2*e45 + 
     -     e12*e16*e26*e45**2 - e15*e16*e24*e25*e46 + 
     -     e14*e16*e25**2*e46 + e15**2*e24*e26*e46 - 
     -     e14*e15*e25*e26*e46 - e12*e16*e25*e45*e46 - 
     -     e12*e15*e26*e45*e46 + e12*e15*e25*e46**2 + 
     -     e15*e16*e24**2*e56 - e14*e16*e24*e25*e56 - 
     -     e14*e15*e24*e26*e56 + e14**2*e25*e26*e56 - 
     -     e12*e16*e24*e45*e56 - e12*e14*e26*e45*e56 - 
     -     e12*e15*e24*e46*e56 - e12*e14*e25*e46*e56 + 
     -     e12**2*e45*e46*e56 +  e12*e14*e24*e56**2)/units5
      f12345 = (e15**2*e23*e24*e34 - 
     -     e14*e15*e23*e25*e34 - e13*e15*e24*e25*e34 + 
     -     e13*e14*e25**2*e34 + e12*e15*e25*e34**2 - 
     -     e14*e15*e23*e24*e35 + e13*e15*e24**2*e35 + 
     -     e14**2*e23*e25*e35 - e13*e14*e24*e25*e35 - 
     -     e12*e15*e24*e34*e35 - e12*e14*e25*e34*e35 + 
     -     e12*e14*e24*e35**2 + e14*e15*e23**2*e45 - 
     -     e13*e15*e23*e24*e45 - e13*e14*e23*e25*e45 + 
     -     e13**2*e24*e25*e45 - e12*e15*e23*e34*e45 - 
     -     e12*e13*e25*e34*e45 - e12*e14*e23*e35*e45 - 
     -     e12*e13*e24*e35*e45 + e12**2*e34*e35*e45 + 
     -     e12*e13*e23*e45**2)/units5
      f12346 = (e16**2*e23*e24*e34 - 
     -     e14*e16*e23*e26*e34 - e13*e16*e24*e26*e34 + 
     -     e13*e14*e26**2*e34 + e12*e16*e26*e34**2 - 
     -     e14*e16*e23*e24*e36 + e13*e16*e24**2*e36 + 
     -     e14**2*e23*e26*e36 - e13*e14*e24*e26*e36 - 
     -     e12*e16*e24*e34*e36 - e12*e14*e26*e34*e36 + 
     -     e12*e14*e24*e36**2 + e14*e16*e23**2*e46 - 
     -     e13*e16*e23*e24*e46 - e13*e14*e23*e26*e46 + 
     -     e13**2*e24*e26*e46 - e12*e16*e23*e34*e46 - 
     -     e12*e13*e26*e34*e46 - e12*e14*e23*e36*e46 - 
     -     e12*e13*e24*e36*e46 + e12**2*e34*e36*e46 + 
     -     e12*e13*e23*e46**2)/units5
      f12356 = (e16**2*e23*e25*e35 - e15*e16*e23*e26*e35 - 
     -     e13*e16*e25*e26*e35 + e13*e15*e26**2*e35 + 
     -     e12*e16*e26*e35**2 - e15*e16*e23*e25*e36 + 
     -     e13*e16*e25**2*e36 + e15**2*e23*e26*e36 - 
     -     e13*e15*e25*e26*e36 - e12*e16*e25*e35*e36 - 
     -     e12*e15*e26*e35*e36 + e12*e15*e25*e36**2 + 
     -     e15*e16*e23**2*e56 - e13*e16*e23*e25*e56 -
     -     e13*e15*e23*e26*e56 + e13**2*e25*e26*e56 - 
     -     e12*e16*e23*e35*e56 - e12*e13*e26*e35*e56 -
     -     e12*e15*e23*e36*e56 - e12*e13*e25*e36*e56 +
     -     e12**2*e35*e36*e56 + e12*e13*e23*e56**2)/units5
      f01 = e16**2*e25**2*e34**2 + e15**2*e26**2*e34**2 + 
     -     e16**2*e24**2*e35**2 + e14**2*e26**2*e35**2 + 
     -     e15**2*e24**2*e36**2 + e14**2*e25**2*e36**2 +
     -     e16**2*e23**2*e45**2 + e12**2*e36**2*e45**2 + 
     -     e13**2*e26**2*e45**2 + e14**2*e23**2*e56**2 +
     -     e13**2*e25**2*e46**2 + e15**2*e23**2*e46**2 +
     -     e13**2*e24**2*e56**2 + e12**2*e34**2*e56**2 +
     -     e12**2*e35**2*e46**2
      f04 = e12*e13*e23*e45*e46*e56 + e12*e14*e24*e35*e36*e56 + 
     -     e12*e16*e26*e34*e35*e45 + e13*e16*e24*e25*e36*e45 + 
     -     e14*e15*e23*e26*e36*e45 + e14*e16*e23*e25*e35*e46 + 
     -     e13*e15*e24*e26*e35*e46 + e12*e15*e25*e34*e36*e46 + 
     -     e15*e16*e23*e24*e34*e56 + e13*e14*e25*e26*e34*e56
      f021 = -e15*e16*e25*e26*e34**2 - e16**2*e24*e25*e34*e35 + 
     -     e15*e16*e24*e26*e34*e35 + e14*e16*e25*e26*e34*e35 - 
     -     e14*e15*e26**2*e34*e35 -  e14*e16*e24*e26*e35**2 + 
     -     e15*e16*e24*e25*e34*e36 - e14*e16*e25**2*e34*e36 - 
     -     e15**2*e24*e26*e34*e36 +  e14*e15*e25*e26*e34*e36 - 
     -     e15*e16*e24**2*e35*e36
      f022 = e14*e16*e24*e25*e35*e36 + e14*e15*e24*e26*e35*e36 - 
     -     e14**2*e25*e26*e35*e36 - e14*e15*e24*e25*e36**2 - 
     -     e16**2*e23*e25*e34*e45 + e15*e16*e23*e26*e34*e45 + 
     -     e13*e16*e25*e26*e34*e45 - e13*e15*e26**2*e34*e45 - 
     -     e16**2*e23*e24*e35*e45 + e14*e16*e23*e26*e35*e45 + 
     -     e13*e16*e24*e26*e35*e45 - e13*e14*e26**2*e35*e45 
      f023 = e15*e16*e23*e24*e36*e45 + e14*e16*e23*e25*e36*e45 + 
     -     e13*e15*e24*e26*e36*e45 + e13*e14*e25*e26*e36*e45 + 
     -     e12*e16*e25*e34*e36*e45 + e12*e15*e26*e34*e36*e45 + 
     -     e12*e16*e24*e35*e36*e45 + e12*e14*e26*e35*e36*e45 - 
     -     e12*e15*e24*e36**2*e45 - e12*e14*e25*e36**2*e45 - 
     -     e13*e16*e23*e26*e45**2 - e12*e16*e23*e36*e45**2 - 
     -     e12*e13*e26*e36*e45**2 + e15*e16*e23*e25*e34*e46 - 
     -     e13*e16*e25**2*e34*e46 - e15**2*e23*e26*e34*e46 
      f024 = e13*e15*e25*e26*e34*e46 + e15*e16*e23*e24*e35*e46 + 
     -     e13*e16*e24*e25*e35*e46 + e14*e15*e23*e26*e35*e46 + 
     -     e13*e14*e25*e26*e35*e46 + e12*e16*e25*e34*e35*e46 + 
     -     e12*e15*e26*e34*e35*e46 - e12*e16*e24*e35**2*e46 - 
     -     e12*e14*e26*e35**2*e46 - e15**2*e23*e24*e36*e46 + 
     -     e14*e15*e23*e25*e36*e46 + e13*e15*e24*e25*e36*e46 - 
     -     e13*e14*e25**2*e36*e46 
      f025 = e12*e15*e24*e35*e36*e46 + e12*e14*e25*e35*e36*e46 - 
     -     e15*e16*e23**2*e45*e46 + e13*e16*e23*e25*e45*e46 + 
     -     e13*e15*e23*e26*e45*e46 - e13**2*e25*e26*e45*e46 + 
     -     e12*e16*e23*e35*e45*e46 + e12*e13*e26*e35*e45*e46 + 
     -     e12*e15*e23*e36*e45*e46 + e12*e13*e25*e36*e45*e46 - 
     -     e12**2*e35*e36*e45*e46 - e13*e15*e23*e25*e46**2 - 
     -     e12*e15*e23*e35*e46**2 - e12*e13*e25*e35*e46**2
      f026 = e14*e16*e23*e25*e34*e56 + e13*e16*e24*e25*e34*e56 + 
     -     e14*e15*e23*e26*e34*e56 + e13*e15*e24*e26*e34*e56 - 
     -     e12*e16*e25*e34**2*e56 - e12*e15*e26*e34**2*e56 + 
     -     e14*e16*e23*e24*e35*e56 - e13*e16*e24**2*e35*e56 - 
     -     e14**2*e23*e26*e35*e56 + e13*e14*e24*e26*e35*e56 + 
     -     e12*e16*e24*e34*e35*e56 + e12*e14*e26*e34*e35*e56 
      f027 = e14*e15*e23*e24*e36*e56 - e13*e15*e24**2*e36*e56 - 
     -     e14**2*e23*e25*e36*e56 + e13*e14*e24*e25*e36*e56 + 
     -     e12*e15*e24*e34*e36*e56 + e12*e14*e25*e34*e36*e56 - 
     -     e14*e16*e23**2*e45*e56 + e13*e16*e23*e24*e45*e56 + 
     -     e13*e14*e23*e26*e45*e56 - e13**2*e24*e26*e45*e56 + 
     -     e12*e16*e23*e34*e45*e56 + e12*e13*e26*e34*e45*e56 
      f028 = e12*e14*e23*e36*e45*e56 + e12*e13*e24*e36*e45*e56 - 
     -     e12**2*e34*e36*e45*e56 - e14*e15*e23**2*e46*e56 + 
     -     e13*e15*e23*e24*e46*e56 + e13*e14*e23*e25*e46*e56 - 
     -     e13**2*e24*e25*e46*e56 + e12*e15*e23*e34*e46*e56 + 
     -     e12*e13*e25*e34*e46*e56 + e12*e14*e23*e35*e46*e56 + 
     -     e12*e13*e24*e35*e46*e56 - e12**2*e34*e35*e46*e56 - 
     -     e13*e14*e23*e24*e56**2 - e12*e14*e23*e34*e56**2 - 
     -     e12*e13*e24*e34*e56**2 
      f123456 = ( f01 - 4.0d0*f04
     -     + 2.d0*( f021+f022+f023+f024+f025+f026+f027+f028 ) )/units6
      salida = xr + ( e12**2*t1*t2 + e13**2*t1*t3 + e23**2*t2*t3 +
     -     e24**2*t2*t4 + e14**2*t1*t4 + e56**2*t5*t6 + 
     -     e34**2*t3*t4 + e25**2*t2*t5 + e15**2*t1*t5 +
     -     e35**2*t3*t5 + e45**2*t4*t5 + e36**2*t3*t6 +
     -     e16**2*t1*t6 + e26**2*t2*t6 + e46**2*t4*t6 )/units2 +
     -     2.d0*xi*( e12*e13*e23*t1*t2*t3 + e12*e14*e24*t1*t2*t4 + 
     -     e13*e14*e34*t1*t3*t4 + e23*e24*e34*t2*t3*t4 + 
     -     e12*e15*e25*t1*t2*t5 + e13*e15*e35*t1*t3*t5 + 
     -     e23*e25*e35*t2*t3*t5 + e14*e15*e45*t1*t4*t5 + 
     -     e24*e25*e45*t2*t4*t5 + e34*e35*e45*t3*t4*t5 + 
     -     e12*e16*e26*t1*t2*t6 + e13*e16*e36*t1*t3*t6 + 
     -     e23*e26*e36*t2*t3*t6 + e14*e16*e46*t1*t4*t6 + 
     -     e24*e26*e46*t2*t4*t6 + e34*e36*e46*t3*t4*t6 + 
     -     e15*e16*e56*t1*t5*t6 + e25*e26*e56*t2*t5*t6 + 
     -     e35*e36*e56*t3*t5*t6 + e45*e46*e56*t4*t5*t6 )/units3 +
     -     2.d0*xi*( t12345*f12345 + t12346*f12346 + t12356*f12356 +
     -     t12456*f12456 + t13456*f13456 + t23456*f23456 ) +
     -     ( f1234*t1234 + t1235*f1235 + t1236*f1236 + t1245*f1245 + 
     -     t1246*f1246 + t1256*f1256 + t1345*f1345 + t1346*f1346 +
     -     t1356*f1356 + t1456*f1456 + t2345*f2345 + t3456*f3456 +
     -     t2346*f2346 + t2356*f2356 + t2456*f2456 )/units4 +
     -     t123456*f123456
      output = salida
      return
      end function determinante
!     -----------------------------------------------------------------------
      subroutine assignmasses
      use masses
      use channelassignment
      use matrixstorage
      implicit none
!     masses and widths
!     mesons
      masas(jprot)   = mp       ! proton
      widths(jprot)  = gprot
      masas(jkaon)   = mkm      ! negative kaon
      widths(jkaon)  = gkm
      masas(jpion)   = mpi      ! pion
      widths(jpion)  = gpi
      masas(jeta)    = meta     ! eta
      widths(jeta)   = geta
      masas(jkstar)  = mkstar   ! kbar*
      widths(jkstar) = gkstar
      masas(jpion0)  = mpi0     ! neutral pion
      widths(jpion0) = gpi0
      masas(jkav)    = mk       ! averaged kaon
      widths(jkav)   = gk
      masas(jdummy)  = mdummy   ! dummy meson
      widths(jdummy) = gdummy
!     baryons
      masas(jnucl)        = mnuc     ! nucleon
      widths(jnucl)       = gnuc
      masas(jdelta)       = mdelta   ! delta
      widths(jdelta)      = gdelta
      masas(jlambda)      = mlambda  ! lambda
      widths(jlambda)     = glambda
      masas(jlambda1520)  = mlambda1520 ! lambda (1520)
      widths(jlambda1520) = glambda1520
      masas(jsigma)       = msigma      ! sigma
      widths(jsigma)      = gsigma
      masas(jsigma1385)   = msigma1385  ! sigma (1385)
      widths(jsigma1385)  = gsigma1385
      return
      end subroutine assignmasses
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     resonance -> j1 + j2
!     j1 -> j3 + j4
!     j2 -> j5 + j6
      subroutine assignchan(k,j1,j2,j3,j4,j5,j6)
      use matrixstorage
      implicit none
      integer :: k,j1,j2,j3,j4,j5,j6
      if (j1.ne.0) channelmass(k)           = masas(j1)
      if (j2.ne.0) channelmass(ndecays+k)   = masas(j2)
      if (j3.ne.0) channelmass(2*ndecays+k) = masas(j3)
      if (j4.ne.0) channelmass(3*ndecays+k) = masas(j4)
      if (j5.ne.0) channelmass(4*ndecays+k) = masas(j5)
      if (j6.ne.0) channelmass(5*ndecays+k) = masas(j6)
      if (j1.ne.0) channelwidths(k)         = widths(j1)
      if (j2.ne.0) channelwidths(ndecays+k) = widths(j2)
      return
      end subroutine assignchan
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine channelsassignment
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k,a0
      a0 = 0
!     channels
      do k=1,ndecaysm
         channelmass(k) = 0.d0
      enddo
      do k=1,ndecaysw
         channelwidths(k) = 0.d0
      enddo
      call assignchan
     c     (kbarn,jkaon,jprot,a0,a0,a0,a0) ! kbarn
      call assignchan
     c     (pilambda,jpion0,jlambda,a0,a0,a0,a0) ! pilambda
      call assignchan
     c     (pisigma,jpion0,jsigma,a0,a0,a0,a0) ! pisigma
      call assignchan
     c     (etalambda,jeta,jlambda,a0,a0,a0,a0) ! etalambda
      call assignchan
     c     (etasigma,jeta,jsigma,a0,a0,a0,a0) ! etasigma
      call assignchan
     c     (kbar1starn,jkstar,jnucl,jkaon,jpion0,a0,a0) ! kbar1starn
      call assignchan
     c     (kbar3starnm,jkstar,jnucl,jkaon,jpion0,a0,a0) ! kbar3starnm
      call assignchan
     c     (kbar3starnp,jkstar,jnucl,jkaon,jpion0,a0,a0) ! kbar3starnm
      call assignchan
     c     (kbardeltam,jdelta,jkaon,jnucl,jpion0,a0,a0) ! (kbardelta)-
      call assignchan
     c     (kbardeltap,jdelta,jkaon,jnucl,jpion0,a0,a0) ! (kbardelta)-
      call assignchan
     c     (pisigma1385m,jsigma1385,jpion0,jlambda,jpion0,a0,a0) ! (pisigma1385)-
      call assignchan
     c     (pisigma1385p,jsigma1385,jpion0,jlambda,jpion0,a0,a0) ! (pisigma1385)+
      call assignchan
     c     (pilambda1520m,jlambda1520,jpion0,jkaon,jnucl,a0,a0) ! (pilambda1520)-
      call assignchan
     c     (pilambda1520p,jlambda1520,jpion0,jkaon,jnucl,a0,a0) ! (pilambda1520)+
      call assignchan
     c     (dummylambda,jdummy,jlambda,jpion0,jpion0,a0,a0) ! dummylambda
      call assignchan
     c     (dummysigma,jdummy,jsigma,jpion0,jpion0,a0,a0) ! dummysigma
      return
      end subroutine channelsassignment
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     resonances parameters upload
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine resonancesparametersinput
      use canalesop
      use channelassignment
      use matrixstorage
      implicit none
      integer :: indres,iloop1,iloop2
      real (8), external :: phi
!     initialization of resonances
      do iloop1=1,nresonances
         masasresonancias(iloop1)  = 0.d0
         momangresonancias(iloop1) = 0
         isresonancia(iloop1) = 0.d0
         do iloop2=1,ndecays
            branchings (iloop1,iloop2) =  0.d0
         enddo
      enddo
!     -----------------------------------------------------------------------
!     upload of the parameters of the resonances and backgrounds
      indres = 0                   ! resonance storage index
!     -----------------------------------------------------------------------
!     isospin 0
!     -----------------------------------------------------------------------
!     s01 res
      indres = indres + 1
      cs01(1) = indres 
      call resonancelambda_s01_1(indres)
      indres = indres + 1
      cs01(2) = indres
      call resonancelambda_s01_2(indres)
      indres = indres + 1
      cs01(3) = indres
      call resonancelambda_s01_3(indres)
      indres = indres + 1
      cs01(4) = indres
      call resonancelambda_s01_4(indres)
      indres = indres + 1
      cs01(5) = indres
      call resonancelambda_s01_5(indres)
      indres = indres + 1
      cs01(6) = indres
      call resonancelambda_s01_6(indres)
!     -----------------------------------------------------------------------
!     p01 res
      indres = indres + 1
      cp01(1) = indres
      call resonancelambda_p01_1(indres)
      indres = indres + 1
      cp01(2) = indres
      call resonancelambda_p01_2(indres)
      indres = indres + 1
      cp01(3) = indres
      call resonancelambda_p01_3(indres)
      indres = indres + 1
      cp01(4) = indres
      call resonancelambda_p01_4(indres)
      indres = indres + 1
      cp01(5) = indres
      call resonancelambda_p01_5(indres)
      indres = indres + 1
      cp01(6) = indres
      call resonancelambda_p01_6(indres)
!     -----------------------------------------------------------------------
!     p03 res
      indres = indres + 1
      cp03(1) = indres
      call resonancelambda_p03_1(indres)
      indres = indres + 1
      cp03(2) = indres
      call resonancelambda_p03_2(indres)
      indres = indres + 1
      cp03(3) = indres
      call resonancelambda_p03_3(indres)
      indres = indres + 1
      cp03(4) = indres
      call resonancelambda_p03_4(indres)
      cp03(5) = 0
      cp03(6) = 0
!     -----------------------------------------------------------------------
!     d03 res
      indres = indres + 1
      cd03(1) = indres
      call resonancelambda_d03_1(indres)
      indres = indres + 1
      cd03(2) = indres
      call resonancelambda_d03_2(indres)
      indres = indres + 1
      cd03(3) = indres
      call resonancelambda_d03_3(indres)
      indres = indres + 1
      cd03(4) = indres
      call resonancelambda_d03_4(indres)
      cd03(5) = 0
      cd03(6) = 0
!     -----------------------------------------------------------------------
!     d05 res
      indres = indres + 1
      cd05(1) = indres
      call resonancelambda_d05_1(indres)
      indres = indres + 1
      cd05(2) = indres
      call resonancelambda_d05_2(indres)
      indres = indres + 1
      cd05(3) = indres
      call resonancelambda_d05_3(indres)
      cd05(4) = 0
      cd05(5) = 0
      cd05(6) = 0
!     -----------------------------------------------------------------------
!     f05 res
      indres = indres + 1
      cf05(1) = indres
      call resonancelambda_f05_1(indres)
      indres = indres + 1
      cf05(2) = indres
      call resonancelambda_f05_2(indres)
      indres = indres + 1
      cf05(3) = indres
      call resonancelambda_f05_3(indres)
      cf05(4) = 0
      cf05(5) = 0
      cf05(6) = 0
!     -----------------------------------------------------------------------
!     f07 res
      indres = indres + 1
      cf07(1) = indres
      call resonancelambda_f07_1(indres)
      indres = indres + 1
      cf07(2) = indres
      call resonancelambda_f07_2(indres)
      cf07(3) = 0
      cf07(4) = 0
      cf07(5) = 0
      cf07(6) = 0
!     -----------------------------------------------------------------------
!     g07 res
      indres = indres + 1
      cg07(1) = indres
      call resonancelambda_g07_1(indres)
      indres = indres + 1
      cg07(2) = indres
      call resonancelambda_g07_2(indres)
      cg07(3) = 0
      cg07(4) = 0
      cg07(5) = 0
      cg07(6) = 0
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     isospin 1, sigma resonances
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     s11 res
      indres = indres + 1
      cs11(1) = indres
      call resonancesigma_s11_1(indres)
      indres = indres + 1
      cs11(2) = indres
      call resonancesigma_s11_2(indres)
      indres = indres + 1
      cs11(3) = indres
      call resonancesigma_s11_3(indres)
      indres = indres + 1
      cs11(4) = indres
      call resonancesigma_s11_4(indres)
      indres = indres + 1
      cs11(5) = indres
      call resonancesigma_s11_5(indres)
      indres = indres + 1
      cs11(6) = indres
      call resonancesigma_s11_6(indres)
!     -----------------------------------------------------------------------
!     p11 res
      indres = indres + 1
      cp11(1) = indres
      call resonancesigma_p11_1(indres)
      indres = indres + 1
      cp11(2) = indres
      call resonancesigma_p11_2(indres)
      indres = indres + 1
      cp11(3) = indres
      call resonancesigma_p11_3(indres)
      indres = indres + 1
      cp11(4) = indres
      call resonancesigma_p11_4(indres)
      indres = indres + 1
      cp11(5) = indres
      call resonancesigma_p11_5(indres)
      cp11(6) = 0
!     -----------------------------------------------------------------------
!     p13 res
      indres = indres + 1
      cp13(1) = indres
      call resonancesigma_p13_1(indres)
      indres = indres + 1
      cp13(2) = indres
      call resonancesigma_p13_2(indres)
      indres = indres + 1
      cp13(3) = indres
      call resonancesigma_p13_3(indres)
      indres = indres + 1
      cp13(4) = indres
      call resonancesigma_p13_4(indres)
      indres = indres + 1
      cp13(5) = indres
      call resonancesigma_p13_5(indres)
      indres = indres + 1
      cp13(6) = indres
      call resonancesigma_p13_6(indres)
!     -----------------------------------------------------------------------
!     d13 res
      indres = indres + 1
      cd13(1) = indres
      call resonancesigma_d13_1(indres)
      indres = indres + 1
      cd13(2) = indres 
      call resonancesigma_d13_2(indres)
      indres = indres + 1 
      cd13(3) = indres 
      call resonancesigma_d13_3(indres)
      cd13(4) = 0
      cd13(5) = 0
      cd13(6) = 0
!     -----------------------------------------------------------------------
!     d15 res
      indres = indres + 1
      cd15(1) = indres
      call resonancesigma_d15_1(indres)
      indres = indres + 1
      cd15(2) = indres
      call resonancesigma_d15_2(indres)
      indres = indres + 1
      cd15(3) = indres
      call resonancesigma_d15_3(indres)
      cd15(4) = 0
      cd15(5) = 0
      cd15(6) = 0
!     -----------------------------------------------------------------------
!     f15 res
      indres = indres + 1
      cf15(1) = indres
      call resonancesigma_f15_1(indres)
      indres = indres + 1
      cf15(2) = indres
      call resonancesigma_f15_2(indres)
      indres = indres + 1
      cf15(3) = indres
      call resonancesigma_f15_3(indres)
      cf15(4) = 0
      cf15(5) = 0
      cf15(6) = 0
!     -----------------------------------------------------------------------
!     f17 res
      indres = indres + 1
      cf17(1) = indres
      call resonancesigma_f17_1(indres)
      indres = indres + 1
      cf17(2) = indres
      call resonancesigma_f17_2(indres)
      cf17(3) = 0
      cf17(4) = 0
      cf17(5) = 0
      cf17(6) = 0
!     -----------------------------------------------------------------------
!     g17 res
      indres = indres + 1
      cg17(1) = indres
      call resonancesigma_g17_1(indres)
      cg17(2) = 0
      cg17(3) = 0
      cg17(4) = 0
      cg17(5) = 0
      cg17(6) = 0
      return
      end subroutine resonancesparametersinput
!     -----------------------------------------------------------------------
!     resonance widths normalization (not wasting time)
!     -----------------------------------------------------------------------
      subroutine phiresnorm(a1,a2,a3,a4,a5,a6)
      use canalesop
      use channelassignment
      use matrixstorage
      implicit none
      integer :: a1,a2,a3,a4,a5,a6,iloop1,iloop2,l
      real (8) :: m0,m1,m2,res
      real (8), external :: phinormfunc
!     1st resonance
      iloop1 = a1
      l  = momangresonancias(iloop1) 
      m0 = masasresonancias(iloop1)
      res = isresonancia(iloop1)
      if (res.gt.0.d0) then
         do iloop2=1,matrixsize
            m1 = channelmass(iloop2)
            m2 = channelmass(ndecays+iloop2)
            phimres(iloop1,iloop2) = phinormfunc(l,iloop2,m0*m0)
         enddo
      endif
!     2nd resonance
      if (a2.ne.0) then
         iloop1 = a2
         l  = momangresonancias(iloop1) 
         m0 = masasresonancias(iloop1)
         res = isresonancia(iloop1)
         if (res.gt.0.d0) then
            do iloop2=1,matrixsize
               m1 = channelmass(iloop2)
               m2 = channelmass(ndecays+iloop2)
               phimres(iloop1,iloop2) = phinormfunc(l,iloop2,m0*m0)
            enddo
         endif
      endif
!     3rd resonance
      if (a3.ne.0) then
         iloop1 = a3
         l  = momangresonancias(iloop1) 
         m0 = masasresonancias(iloop1)
         res = isresonancia(iloop1)
         if (res.gt.0.d0) then
            do iloop2=1,matrixsize
               m1 = channelmass(iloop2)
               m2 = channelmass(ndecays+iloop2)
               phimres(iloop1,iloop2) = phinormfunc(l,iloop2,m0*m0)
            enddo
         endif
      endif
!     4th resonance
      if (a4.ne.0) then
         iloop1 = a4
         l  = momangresonancias(iloop1) 
         m0 = masasresonancias(iloop1)
         res = isresonancia(iloop1)
         if (res.gt.0.d0) then
            do iloop2=1,matrixsize
               m1 = channelmass(iloop2)
               m2 = channelmass(ndecays+iloop2)
               phimres(iloop1,iloop2) = phinormfunc(l,iloop2,m0*m0)
            enddo
         endif
      endif
!     5th resonance
      if (a5.ne.0) then
         iloop1 = a5
         l  = momangresonancias(iloop1) 
         m0 = masasresonancias(iloop1)
         res = isresonancia(iloop1)
         if (res.gt.0.d0) then
            do iloop2=1,matrixsize
               m1 = channelmass(iloop2)
               m2 = channelmass(ndecays+iloop2)
               phimres(iloop1,iloop2) = phinormfunc(l,iloop2,m0*m0)
            enddo
         endif
      endif
!     6th resonance
      if (a6.ne.0) then
         iloop1 = a6
         l  = momangresonancias(iloop1) 
         m0 = masasresonancias(iloop1)
         res = isresonancia(iloop1)
         if (res.gt.0.d0) then
            do iloop2=1,matrixsize
               m1 = channelmass(iloop2)
               m2 = channelmass(ndecays+iloop2)
               phimres(iloop1,iloop2) = phinormfunc(l,iloop2,m0*m0)
            enddo
         endif
      endif
      return
      end subroutine phiresnorm
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!
!     lambda resonances
!
!     -----------------------------------------------------------------------
!     s01 resonances 
!     -----------------------------------------------------------------------
      subroutine resonancelambda_s01_1(k) ! lambda(1405)s01
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)  = 0
      isresonancia(k)       = 1.d0
      masasresonancias(k)   = coeficiente(1) 
      branchings(k,kbarn)   = coeficiente(2) ! kbarn
      branchings(k,pisigma) = coeficiente(3) ! pisigma
      return
      end subroutine resonancelambda_s01_1
!     -----------------------------------------------------------------------
      subroutine resonancelambda_s01_2(k) ! lambda(1670)s01
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       = 0
      isresonancia(k)            = 1.d0
      masasresonancias(k)        = coeficiente(4) 
      branchings(k,kbarn)        = coeficiente(5) ! kbarn
      branchings(k,pisigma)      = coeficiente(6) ! pisigma
      branchings(k,etalambda)    = coeficiente(7) ! etalambda
      branchings(k,kbar1starn)   = coeficiente(8) ! kbar1starn
      branchings(k,kbar3starnp)  = coeficiente(9) ! (kbar3starn)d
      branchings(k,pisigma1385p) = coeficiente(10) ! (pisigma1385)d
      branchings(k,dummylambda)  = coeficiente(11) ! dummylambda
      return
      end subroutine resonancelambda_s01_2
!     -----------------------------------------------------------------------
      subroutine resonancelambda_s01_3(k) ! lambda(1780)s01
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       = 0
      isresonancia(k)            = 1.d0
      masasresonancias(k)        = coeficiente(12) 
      branchings(k,kbarn)        = coeficiente(13) ! kbarn
      branchings(k,pisigma)      = coeficiente(14) ! pisigma
      branchings(k,etalambda)    = coeficiente(15) ! etalambda
      branchings(k,kbar1starn)   = coeficiente(16) ! kbar1starn
      branchings(k,kbar3starnp)  = coeficiente(17) ! (kbar3starn)d
      branchings(k,pisigma1385p) = coeficiente(18) ! (pisigma1385)d
      branchings(k,dummylambda)  = coeficiente(19) ! dummylambda
      return
      end subroutine resonancelambda_s01_3
!     -----------------------------------------------------------------------
      subroutine resonancelambda_s01_4(k) ! lambda(2020)s01
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       = 0
      isresonancia(k)            = 1.d0
      masasresonancias(k)        = coeficiente(20)
      branchings(k,kbarn)        = coeficiente(21) ! kbarn
      branchings(k,pisigma)      = coeficiente(22) ! pisigma
      branchings(k,etalambda)    = coeficiente(23) ! etalambda
      branchings(k,kbar1starn)   = coeficiente(24) ! kbar1starn
      branchings(k,kbar3starnp)  = coeficiente(25) ! (kbar3starn)d
      branchings(k,pisigma1385p) = coeficiente(26) ! (pisigma1385)d
      branchings(k,dummylambda)  = coeficiente(27) ! dummylambda
      return
      end subroutine resonancelambda_s01_4
!     -----------------------------------------------------------------------
      subroutine resonancelambda_s01_5(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       = 0
      isresonancia(k)            = -1.d0
      masasresonancias(k)        = coeficiente(28)
      branchings(k,kbarn)        = coeficiente(29) ! kbarn
      branchings(k,pisigma)      = coeficiente(30) ! pisigma
      branchings(k,etalambda)    = coeficiente(31) ! etalambda
      branchings(k,kbar1starn)   = coeficiente(32) ! kbar1starn
      branchings(k,kbar3starnp)  = coeficiente(33) ! (kbar3starn)d
      branchings(k,pisigma1385p) = coeficiente(34) ! (pisigma1385)d
      branchings(k,dummylambda)  = coeficiente(35) ! dummylambda
      return
      end subroutine resonancelambda_s01_5
!     -----------------------------------------------------------------------
      subroutine resonancelambda_s01_6(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       = 0
      isresonancia(k)            = -1.d0
      masasresonancias(k)        = coeficiente(36)
      branchings(k,kbarn)        = coeficiente(37) ! kbarn
      branchings(k,pisigma)      = coeficiente(38) ! pisigma
      branchings(k,etalambda)    = coeficiente(39) ! etalambda
      branchings(k,kbar1starn)   = coeficiente(40) ! kbar1starn
      branchings(k,kbar3starnp)  = coeficiente(41) ! (kbar3starn)d
      branchings(k,pisigma1385p) = coeficiente(42) ! (pisigma1385)d
      branchings(k,dummylambda)  = coeficiente(43) ! dummylambda
      return
      end subroutine resonancelambda_s01_6
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     p01 resonances
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine resonancelambda_p01_1(k) ! lambda(1600)p01
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       = 1
      isresonancia(k)            = 1.d0
      masasresonancias(k)        = coeficiente(44) 
      branchings(k,kbarn)        = coeficiente(45) ! kbarn
      branchings(k,pisigma)      = coeficiente(46) ! pisigma
      branchings(k,kbar1starn)   = coeficiente(47) ! kbar1starn
      branchings(k,kbar3starnp)  = coeficiente(48) ! (kbar3starn)d
      branchings(k,pisigma1385p) = coeficiente(49) ! (pisigma1385)d
      branchings(k,dummylambda)  = coeficiente(50) ! dummylambda
      return
      end subroutine resonancelambda_p01_1
!     -----------------------------------------------------------------------
      subroutine resonancelambda_p01_2(k) ! lambda(1710)p01
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       = 1
      isresonancia(k)            = 1.d0
      masasresonancias(k)        = coeficiente(51) 
      branchings(k,kbarn)        = coeficiente(52) ! kbarn
      branchings(k,pisigma)      = coeficiente(53) ! pisigma
      branchings(k,kbar1starn)   = coeficiente(54) ! kbar1starn
      branchings(k,kbar3starnp)  = coeficiente(55) ! (kbar3starn)d
      branchings(k,pisigma1385p) = coeficiente(56) ! (pisigma1385)d
      branchings(k,dummylambda)  = coeficiente(57) ! dummylambda
      return
      end subroutine resonancelambda_p01_2
!     -----------------------------------------------------------------------
      subroutine resonancelambda_p01_3(k) ! lambda(1810)p01
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       = 1
      isresonancia(k)            = 1.d0
      masasresonancias(k)        = coeficiente(58) 
      branchings(k,kbarn)        = coeficiente(59) ! kbarn
      branchings(k,pisigma)      = coeficiente(60) ! pisigma
      branchings(k,kbar1starn)   = coeficiente(61) ! kbar1starn
      branchings(k,kbar3starnp)  = coeficiente(62) ! (kbar3starn)d
      branchings(k,pisigma1385p) = coeficiente(63) ! (pisigma1385)d
      branchings(k,dummylambda)  = coeficiente(64) ! dummylambda
      return
      end subroutine resonancelambda_p01_3
!     -----------------------------------------------------------------------
      subroutine resonancelambda_p01_4(k) ! lambda(2200)p01
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       = 1
      isresonancia(k)            = 1.d0
      masasresonancias(k)        = coeficiente(65) 
      branchings(k,kbarn)        = coeficiente(66) ! kbarn
      branchings(k,pisigma)      = coeficiente(67) ! pisigma
      branchings(k,kbar1starn)   = coeficiente(68) ! kbar1starn
      branchings(k,kbar3starnp)  = coeficiente(69) ! (kbar3starn)d
      branchings(k,pisigma1385p) = coeficiente(70) ! (pisigma1385)d
      branchings(k,dummylambda)  = coeficiente(71) ! dummylambda
      return
      end subroutine resonancelambda_p01_4
!     -----------------------------------------------------------------------
      subroutine resonancelambda_p01_5(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       =  1
      isresonancia(k)            = -1.d0
      masasresonancias(k)        = coeficiente(72)
      branchings(k,kbarn)        = coeficiente(73) ! kbarn
      branchings(k,pisigma)      = coeficiente(74) ! pisigma
      branchings(k,kbar1starn)   = coeficiente(75) ! kbar1starn
      branchings(k,kbar3starnp)  = coeficiente(76) ! (kbar3starn)d
      branchings(k,pisigma1385p) = coeficiente(77) ! (pisigma1385)d
      branchings(k,dummylambda)  = coeficiente(78) ! dummylambda
      return
      end subroutine resonancelambda_p01_5
!     -----------------------------------------------------------------------
      subroutine resonancelambda_p01_6(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       =  1
      isresonancia(k)            = -1.d0
      masasresonancias(k)        = coeficiente(79)
      branchings(k,kbarn)        = coeficiente(80) ! kbarn
      branchings(k,pisigma)      = coeficiente(81) ! pisigma
      branchings(k,kbar1starn)   = coeficiente(82) ! kbar1starn
      branchings(k,kbar3starnp)  = coeficiente(83) ! (kbar3starn)d
      branchings(k,pisigma1385p) = coeficiente(84) ! (pisigma1385)d
      branchings(k,dummylambda)  = coeficiente(85) ! dummylambda
      return
      end subroutine resonancelambda_p01_6
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     p03 resonances
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine resonancelambda_p03_1(k) ! lambda(1890)p03
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       = 1
      isresonancia(k)            = 1.d0
      masasresonancias(k)        = coeficiente(86)
      branchings(k,kbarn)        = coeficiente(87) ! kbarn
      branchings(k,pisigma)      = coeficiente(88) ! pisigma
      branchings(k,kbar1starn)   = coeficiente(89) ! kbar1starn
      branchings(k,kbar3starnm)  = coeficiente(90) ! (kbar3starn)p
      branchings(k,kbar3starnp)  = coeficiente(91) ! (kbar3starn)f
      branchings(k,pisigma1385m) = coeficiente(92) ! (pisigma1385)p
      branchings(k,pisigma1385p) = coeficiente(93) ! (pisigma1385)f
      branchings(k,dummylambda)  = coeficiente(94) ! dummylambda
      return
      end subroutine resonancelambda_p03_1
!     -----------------------------------------------------------------------
      subroutine resonancelambda_p03_2(k) ! lambda(2300)p03
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       = 1
      isresonancia(k)            = 1.d0
      masasresonancias(k)        = coeficiente(95) 
      branchings(k,kbarn)        = coeficiente(96) ! kbarn
      branchings(k,pisigma)      = coeficiente(97) ! pisigma
      branchings(k,kbar1starn)   = coeficiente(98) ! kbar1starn
      branchings(k,kbar3starnm)  = coeficiente(99) ! (kbar3starn)p
      branchings(k,kbar3starnp)  = coeficiente(100) ! (kbar3starn)f
      branchings(k,pisigma1385m) = coeficiente(101) ! (pisigma1385)p
      branchings(k,pisigma1385p) = coeficiente(102) ! (pisigma1385)f
      branchings(k,dummylambda)  = coeficiente(103) ! dummylambda
      return
      end subroutine resonancelambda_p03_2
!     -----------------------------------------------------------------------
      subroutine resonancelambda_p03_3(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       =  1
      isresonancia(k)            = -1.d0
      masasresonancias(k)        = coeficiente(104) 
      branchings(k,kbarn)        = coeficiente(105) ! kbarn
      branchings(k,pisigma)      = coeficiente(106) ! pisigma
      branchings(k,kbar1starn)   = coeficiente(107) ! kbar1starn
      branchings(k,kbar3starnm)  = coeficiente(108) ! (kbar3starn)p
      branchings(k,kbar3starnp)  = coeficiente(109) ! (kbar3starn)f
      branchings(k,pisigma1385m) = coeficiente(110) ! (pisigma1385)p
      branchings(k,pisigma1385p) = coeficiente(111) ! (pisigma1385)f
      branchings(k,dummylambda)  = coeficiente(112) ! dummylambda
      return
      end subroutine resonancelambda_p03_3
!     -----------------------------------------------------------------------
      subroutine resonancelambda_p03_4(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       =  1
      isresonancia(k)            = -1.d0
      masasresonancias(k)        = coeficiente(113) 
      branchings(k,kbarn)        = coeficiente(114) ! kbarn
      branchings(k,pisigma)      = coeficiente(115) ! pisigma
      branchings(k,kbar1starn)   = coeficiente(116) ! kbar1starn
      branchings(k,kbar3starnm)  = coeficiente(117) ! (kbar3starn)p
      branchings(k,kbar3starnp)  = coeficiente(118) ! (kbar3starn)f
      branchings(k,pisigma1385m) = coeficiente(119) ! (pisigma1385)p
      branchings(k,pisigma1385p) = coeficiente(120) ! (pisigma1385)f
      branchings(k,dummylambda)  = coeficiente(121) ! dummylambda
      return
      end subroutine resonancelambda_p03_4
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     d03 resonances
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine resonancelambda_d03_1(k) ! lambda(1520)d03
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      real (8) :: sum
      momangresonancias(k)       =  2
      isresonancia(k)            =  1.d0
      masasresonancias(k)        = coeficiente(122) 
      branchings(k,kbarn)        = coeficiente(123) ! kbarn
      branchings(k,pisigma)      = coeficiente(124) ! pisigma
      branchings(k,kbar1starn)   = coeficiente(125) ! kbar1starn ! subthreshold channel
      branchings(k,pisigma1385m) = coeficiente(126) ! (pisigma1385)s
      branchings(k,pisigma1385p) = coeficiente(127) ! (pisigma1385)d
      branchings(k,dummylambda)  = coeficiente(128) ! dummylambda
      return
      end subroutine resonancelambda_d03_1
!     -----------------------------------------------------------------------
      subroutine resonancelambda_d03_2(k) ! lambda(1690)d03
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       =  2
      isresonancia(k)            =  1.d0
      masasresonancias(k)        = coeficiente(129)  
      branchings(k,kbarn)        = coeficiente(130) ! kbarn
      branchings(k,pisigma)      = coeficiente(131) ! pisigma
      branchings(k,kbar1starn)   = coeficiente(132) ! kbar1starn
      branchings(k,pisigma1385m) = coeficiente(133) ! (pisigma1385)s
      branchings(k,pisigma1385p) = coeficiente(134) ! (pisigma1385)d
      branchings(k,dummylambda)  = coeficiente(135) ! dummylambda
      return
      end subroutine resonancelambda_d03_2
!     -----------------------------------------------------------------------
      subroutine resonancelambda_d03_3(k) ! lambda(2050)d03
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       =  2
      isresonancia(k)            =  1.d0
      masasresonancias(k)        = coeficiente(136)  
      branchings(k,kbarn)        = coeficiente(137) ! kbarn
      branchings(k,pisigma)      = coeficiente(138) ! pisigma
      branchings(k,kbar1starn)   = coeficiente(139) ! kbar1starn
      branchings(k,pisigma1385m) = coeficiente(140) ! (pisigma1385)s
      branchings(k,pisigma1385p) = coeficiente(141) ! (pisigma1385)d
      branchings(k,dummylambda)  = coeficiente(142) ! dummylambda
      return
      end subroutine resonancelambda_d03_3
!     -----------------------------------------------------------------------
      subroutine resonancelambda_d03_4(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       =  2
      isresonancia(k)            = -1.d0
      masasresonancias(k)        = coeficiente(143) 
      branchings(k,kbarn)        = coeficiente(144) ! kbarn
      branchings(k,pisigma)      = coeficiente(145) ! pisigma
      branchings(k,kbar1starn)   = coeficiente(146) ! kbar1starn
      branchings(k,pisigma1385m) = coeficiente(147) ! (pisigma1385)s
      branchings(k,pisigma1385p) = coeficiente(148) ! (pisigma1385)d
      branchings(k,dummylambda)  = coeficiente(149) ! dummylambda
      return
      end subroutine resonancelambda_d03_4
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     d05 resonances
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine resonancelambda_d05_1(k) ! lambda(1830)d05
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       =  2
      isresonancia(k)            =  1.d0
      masasresonancias(k)        = coeficiente(150)  
      branchings(k,kbarn)        = coeficiente(151) ! kbarn
      branchings(k,pisigma)      = coeficiente(152) ! pisigma
      branchings(k,pisigma1385m) = coeficiente(153) ! (pisigma1385)d
      branchings(k,pisigma1385p) = coeficiente(154) ! (pisigma1385)g
      branchings(k,dummylambda)  = coeficiente(155) ! dummylambda
      return
      end subroutine resonancelambda_d05_1
!     -----------------------------------------------------------------------
      subroutine resonancelambda_d05_2(k) ! lambda(2300)d05
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       =  2
      isresonancia(k)            =  1.d0
      masasresonancias(k)        = coeficiente(156)
      branchings(k,kbarn)        = coeficiente(157) ! kbarn
      branchings(k,pisigma)      = coeficiente(158) ! pisigma
      branchings(k,pisigma1385m) = coeficiente(159) ! (pisigma1385)d
      branchings(k,pisigma1385p) = coeficiente(160) ! (pisigma1385)g
      branchings(k,dummylambda)  = coeficiente(161) ! dummylambda
      return
      end subroutine resonancelambda_d05_2
!     -----------------------------------------------------------------------
      subroutine resonancelambda_d05_3(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       =  2
      isresonancia(k)            = -1.d0
      masasresonancias(k)        = coeficiente(162)  
      branchings(k,kbarn)        = coeficiente(163) ! kbarn
      branchings(k,pisigma)      = coeficiente(164) ! pisigma
      branchings(k,pisigma1385m) = coeficiente(165) ! (pisigma1385)d
      branchings(k,pisigma1385p) = coeficiente(166) ! (pisigma1385)g
      branchings(k,dummylambda)  = coeficiente(167) ! dummylambda
      return
      end subroutine resonancelambda_d05_3
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     f05 resonances
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine resonancelambda_f05_1(k) ! lambda(1820)f05
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       =  3
      isresonancia(k)            =  1.d0
      masasresonancias(k)        = coeficiente(168) 
      branchings(k,kbarn)        = coeficiente(169) ! kbarn
      branchings(k,pisigma)      = coeficiente(170) ! pisigma
      branchings(k,kbar1starn)   = coeficiente(171) ! kbar1starn
      branchings(k,kbar3starnm)  = coeficiente(172) ! (kbar3starn)p
      branchings(k,kbar3starnp)  = coeficiente(173) ! (kbar3starn)f
      branchings(k,pisigma1385m) = coeficiente(174) ! (pisigma1385)p
      branchings(k,pisigma1385p) = coeficiente(175) ! (pisigma1385)f
      branchings(k,dummylambda)  = coeficiente(176) ! dummylambda
      return
      end subroutine resonancelambda_f05_1
!     -----------------------------------------------------------------------
      subroutine resonancelambda_f05_2(k) ! lambda(2110)f05
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       = 3
      isresonancia(k)            = 1.d0
      masasresonancias(k)        = coeficiente(177) 
      branchings(k,kbarn)        = coeficiente(178) ! kbarn
      branchings(k,pisigma)      = coeficiente(179) ! pisigma
      branchings(k,kbar1starn)   = coeficiente(180) ! kbar1starn
      branchings(k,kbar3starnm)  = coeficiente(181) ! (kbar3starn)p
      branchings(k,kbar3starnp)  = coeficiente(182) ! (kbar3starn)f
      branchings(k,pisigma1385m) = coeficiente(183) ! (pisigma1385)p
      branchings(k,pisigma1385p) = coeficiente(184) ! (pisigma1385)f
      branchings(k,dummylambda)  = coeficiente(185) ! dummylambda
      return
      end subroutine resonancelambda_f05_2
!     -----------------------------------------------------------------------
      subroutine resonancelambda_f05_3(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)       =  3
      isresonancia(k)            = -1.d0
      masasresonancias(k)        = coeficiente(186)
      branchings(k,kbarn)        = coeficiente(187) ! kbarn
      branchings(k,pisigma)      = coeficiente(188) ! pisigma
      branchings(k,kbar1starn)   = coeficiente(189) ! kbar1starn
      branchings(k,kbar3starnm)  = coeficiente(190) ! (kbar3starn)p
      branchings(k,kbar3starnp)  = coeficiente(191) ! (kbar3starn)f
      branchings(k,pisigma1385p) = coeficiente(192) ! (pisigma1385)p
      branchings(k,pisigma1385m) = coeficiente(193) ! (pisigma1385)f
      branchings(k,dummylambda)  = coeficiente(194) ! dummylambda
      return
      end subroutine resonancelambda_f05_3
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     f07 resonances
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine resonancelambda_f07_1(k) ! lambda(2020)f07
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)      = 3
      isresonancia(k)           = 1.d0
      masasresonancias(k)       = coeficiente(195)
      branchings(k,kbarn)       = coeficiente(196) ! kbarn
      branchings(k,pisigma)     = coeficiente(197) ! pisigma
      branchings(k,kbar1starn)  = coeficiente(198) ! kbar1starn
      branchings(k,dummylambda) = coeficiente(199) ! dummylambda
      return
      end subroutine resonancelambda_f07_1
!     -----------------------------------------------------------------------
      subroutine resonancelambda_f07_2(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)      =  3
      isresonancia(k)           = -1.d0
      masasresonancias(k)       = coeficiente(200)
      branchings(k,kbarn)       = coeficiente(201) ! kbarn
      branchings(k,pisigma)     = coeficiente(202) ! pisigma
      branchings(k,kbar1starn)  = coeficiente(203) ! kbar1starn
      branchings(k,dummylambda) = coeficiente(204) ! dummylambda
      return
      end subroutine resonancelambda_f07_2
!     -----------------------------------------------------------------------
!     g07 resonances
!     -----------------------------------------------------------------------
      subroutine resonancelambda_g07_1(k) ! lambda(2100)g07
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)      =  4
      isresonancia(k)           =  1.d0
      masasresonancias(k)       = coeficiente(205) 
      branchings(k,kbarn)       = coeficiente(206) ! kbarn
      branchings(k,pisigma)     = coeficiente(207) ! pisigma
      branchings(k,kbar1starn)  = coeficiente(208) ! kbar1starn
      branchings(k,kbar3starnm) = coeficiente(209) ! (kbar3starn)d
      branchings(k,kbar3starnp) = coeficiente(210) ! (kbar3starn)g
      branchings(k,dummylambda) = coeficiente(211) ! dummylambda
      return
      end subroutine resonancelambda_g07_1
!     -----------------------------------------------------------------------
      subroutine resonancelambda_g07_2(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)      =  4
      isresonancia(k)           = -1.d0
      masasresonancias(k)       = coeficiente(212) 
      branchings(k,kbarn)       = coeficiente(213) ! kbarn
      branchings(k,pisigma)     = coeficiente(214) ! pisigma
      branchings(k,kbar1starn)  = coeficiente(215) ! kbar1starn
      branchings(k,kbar3starnm) = coeficiente(216) ! (kbar3starn)d
      branchings(k,kbar3starnp) = coeficiente(217) ! (kbar3starn)g
      branchings(k,dummylambda) = coeficiente(218) ! dummylambda
      return
      end subroutine resonancelambda_g07_2
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!
!     sigma resonances
!
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     s11 resonances 
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine resonancesigma_s11_1(k) ! sigma(1620)s11
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  0
      isresonancia(k)             =  1.d0
      masasresonancias(k)         = coeficiente(219) 
      branchings(k,kbarn)         = coeficiente(220) ! kbarn
      branchings(k,pilambda)      = coeficiente(221) ! pilambda
      branchings(k,pisigma)       = coeficiente(222) ! pisigma
      branchings(k,pisigma1385p)  = coeficiente(223) ! (pisigma1385)d
      branchings(k,pilambda1520p) = coeficiente(224) ! (pilambda1520)p
      branchings(k,kbardeltap)    = coeficiente(225) ! (kbardelta)d
      branchings(k,kbar1starn)    = coeficiente(226) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(227) ! (kbar3starn)d
      branchings(k,etasigma)      = coeficiente(228) ! etasigma ! subthreshold channel
      branchings(k,dummysigma)    = coeficiente(229) ! dummysigma
      return
      end subroutine resonancesigma_s11_1
!     -----------------------------------------------------------------------
      subroutine resonancesigma_s11_2(k) ! sigma(1750)s11
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  0
      isresonancia(k)             =  1.d0
      masasresonancias(k)         = coeficiente(230) 
      branchings(k,kbarn)         = coeficiente(231) ! kbarn
      branchings(k,pilambda)      = coeficiente(232) ! pilambda
      branchings(k,pisigma)       = coeficiente(233) ! pisigma
      branchings(k,pisigma1385p)  = coeficiente(234) ! (pisigma1385)d
      branchings(k,pilambda1520p) = coeficiente(235) ! (pilambda1520)p
      branchings(k,kbardeltap)    = coeficiente(236) ! (kbardelta)d
      branchings(k,kbar1starn)    = coeficiente(237) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(238) ! (kbar3starn)d
      branchings(k,etasigma)      = coeficiente(239) ! etasigma ! subthreshold channel
      branchings(k,dummysigma)    = coeficiente(240) ! dummysigma
      return
      end subroutine resonancesigma_s11_2
!     -----------------------------------------------------------------------
      subroutine resonancesigma_s11_3(k) ! sigma(1900)s11
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  0
      isresonancia(k)             =  1.d0
      masasresonancias(k)         = coeficiente(241) 
      branchings(k,kbarn)         = coeficiente(242) ! kbarn
      branchings(k,pilambda)      = coeficiente(243) ! pilambda
      branchings(k,pisigma)       = coeficiente(244) ! pisigma
      branchings(k,pisigma1385p)  = coeficiente(245) ! (pisigma1385)d
      branchings(k,pilambda1520p) = coeficiente(246) ! (pilambda1520)p
      branchings(k,kbardeltap)    = coeficiente(247) ! (kbardelta)d
      branchings(k,kbar1starn)    = coeficiente(248) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(249) ! (kbar3starn)d
      branchings(k,etasigma)      = coeficiente(250) ! etasigma
      branchings(k,dummysigma)    = coeficiente(251) ! dummysigma
      return
      end subroutine resonancesigma_s11_3
!     -----------------------------------------------------------------------
      subroutine resonancesigma_s11_4(k) ! sigma(2000)s11
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  0
      isresonancia(k)             =  1.d0
      masasresonancias(k)         = coeficiente(252) 
      branchings(k,kbarn)         = coeficiente(253) ! kbarn
      branchings(k,pilambda)      = coeficiente(254) ! pilambda
      branchings(k,pisigma)       = coeficiente(255) ! pisigma
      branchings(k,pisigma1385p)  = coeficiente(256) ! (pisigma1385)d
      branchings(k,pilambda1520p) = coeficiente(257) ! (pilambda1520)p
      branchings(k,kbardeltap)    = coeficiente(258) ! (kbardelta)d
      branchings(k,kbar1starn)    = coeficiente(259) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(260) ! (kbar3starn)d
      branchings(k,etasigma)      = coeficiente(261) ! etasigma
      branchings(k,dummysigma)    = coeficiente(262) ! dummysigma
      return
      end subroutine resonancesigma_s11_4
!     -----------------------------------------------------------------------
      subroutine resonancesigma_s11_5(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  0
      isresonancia(k)             = -1.d0
      masasresonancias(k)         = coeficiente(263) 
      branchings(k,kbarn)         = coeficiente(264) ! kbarn
      branchings(k,pilambda)      = coeficiente(265) ! pilambda
      branchings(k,pisigma)       = coeficiente(266) ! pisigma
      branchings(k,pisigma1385p)  = coeficiente(267) ! (pisigma1385)d
      branchings(k,pilambda1520p) = coeficiente(268) ! (pilambda1520)p
      branchings(k,kbardeltap)    = coeficiente(269) ! (kbardelta)d
      branchings(k,kbar1starn)    = coeficiente(270) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(271) ! (kbar3starn)d
      branchings(k,etasigma)      = coeficiente(272) ! etasigma
      branchings(k,dummysigma)    = coeficiente(273) ! dummysigma
      return
      end subroutine resonancesigma_s11_5
!     -----------------------------------------------------------------------
      subroutine resonancesigma_s11_6(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  0
      isresonancia(k)             = -1.d0
      masasresonancias(k)         = coeficiente(274) 
      branchings(k,kbarn)         = coeficiente(275) ! kbarn
      branchings(k,pilambda)      = coeficiente(276) ! pilambda
      branchings(k,pisigma)       = coeficiente(277) ! pisigma
      branchings(k,pisigma1385p)  = coeficiente(278) ! (pisigma1385)d
      branchings(k,pilambda1520p) = coeficiente(279) ! (pilambda1520)p
      branchings(k,kbardeltap)    = coeficiente(280) ! (kbardelta)d
      branchings(k,kbar1starn)    = coeficiente(281) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(282) ! (kbar3starn)d
      branchings(k,etasigma)      = coeficiente(283) ! etasigma
      branchings(k,dummysigma)    = coeficiente(284) ! dummysigma
      return
      end subroutine resonancesigma_s11_6
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     p11 resonances 
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine resonancesigma_p11_1(k) ! sigma(1770)p11
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      isresonancia(k)             =  1.d0
      momangresonancias(k)        =  1
      masasresonancias(k)         = coeficiente(285) 
      branchings(k,kbarn)         = coeficiente(286) ! kbarn
      branchings(k,pilambda)      = coeficiente(287) ! pilambda
      branchings(k,pisigma)       = coeficiente(288) ! pisigma
      branchings(k,pisigma1385p)  = coeficiente(289) ! (pisigma1385)p
      branchings(k,pilambda1520p) = coeficiente(290) ! (pilambda1520)d
      branchings(k,kbardeltap)    = coeficiente(291) ! (kbardelta)p
      branchings(k,kbar1starn)    = coeficiente(292) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(293) ! (kbar3starn)p
      branchings(k,dummysigma)    = coeficiente(294) ! dummysigma
      return
      end subroutine resonancesigma_p11_1
!     -----------------------------------------------------------------------
      subroutine resonancesigma_p11_2(k) ! sigma(1880)p11
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  1
      isresonancia(k)             =  1.d0
      masasresonancias(k)         = coeficiente(295)
      branchings(k,kbarn)         = coeficiente(296) ! kbarn
      branchings(k,pilambda)      = coeficiente(297) ! pilambda
      branchings(k,pisigma)       = coeficiente(298) ! pisigma
      branchings(k,pisigma1385p)  = coeficiente(299) ! (pisigma1385)p
      branchings(k,pilambda1520p) = coeficiente(300) ! (pilambda1520)d
      branchings(k,kbardeltap)    = coeficiente(301) ! (kbardelta)p
      branchings(k,kbar1starn)    = coeficiente(302) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(303) ! (kbar3starn)p
      branchings(k,dummysigma)    = coeficiente(304) ! dummysigma
      return
      end subroutine resonancesigma_p11_2
!     -----------------------------------------------------------------------
      subroutine resonancesigma_p11_3(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  1
      isresonancia(k)             = -1.d0
      masasresonancias(k)         = coeficiente(305)
      branchings(k,kbarn)         = coeficiente(306) ! kbarn
      branchings(k,pilambda)      = coeficiente(307) ! pilambda
      branchings(k,pisigma)       = coeficiente(308) ! pisigma
      branchings(k,pisigma1385p)  = coeficiente(309) ! (pisigma1385)p
      branchings(k,pilambda1520p) = coeficiente(310) ! (pilambda1520)d
      branchings(k,kbardeltap)    = coeficiente(311) ! (kbardelta)p
      branchings(k,kbar1starn)    = coeficiente(312) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(313) ! (kbar3starn)p
      branchings(k,dummysigma)    = coeficiente(314) ! dummysigma
      return
      end subroutine resonancesigma_p11_3

!     -----------------------------------------------------------------------
      subroutine resonancesigma_p11_4(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  1
      isresonancia(k)             = -1.d0
      masasresonancias(k)         = coeficiente(315)
      branchings(k,kbarn)         = coeficiente(316) ! kbarn
      branchings(k,pilambda)      = coeficiente(317) ! pilambda
      branchings(k,pisigma)       = coeficiente(318) ! pisigma
      branchings(k,pisigma1385p)  = coeficiente(319) ! (pisigma1385)p
      branchings(k,pilambda1520p) = coeficiente(320) ! (pilambda1520)d
      branchings(k,kbardeltap)    = coeficiente(321) ! (kbardelta)p
      branchings(k,kbar1starn)    = coeficiente(322) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(323) ! (kbar3starn)p
      branchings(k,dummysigma)    = coeficiente(324) ! dummysigma
      return
      end subroutine resonancesigma_p11_4
!     -----------------------------------------------------------------------
      subroutine resonancesigma_p11_5(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  1
      isresonancia(k)             = -1.d0
      masasresonancias(k)         = coeficiente(325)
      branchings(k,kbarn)         = coeficiente(326) ! kbarn
      branchings(k,pilambda)      = coeficiente(327) ! pilambda
      branchings(k,pisigma)       = coeficiente(328) ! pisigma
      branchings(k,pisigma1385p)  = coeficiente(329) ! (pisigma1385)p
      branchings(k,pilambda1520p) = coeficiente(330) ! (pilambda1520)d
      branchings(k,kbardeltap)    = coeficiente(331) ! (kbardelta)p
      branchings(k,kbar1starn)    = coeficiente(332) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(333) ! (kbar3starn)p
      branchings(k,dummysigma)    = coeficiente(334) ! dummysigma
      return
      end subroutine resonancesigma_p11_5
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     p13 resonances 
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine resonancesigma_p13_1(k) ! sigma(1730)p13
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  1
      isresonancia(k)             =  1.d0
      masasresonancias(k)         = coeficiente(335) 
      branchings(k,kbarn)         = coeficiente(336) ! kbarn
      branchings(k,pilambda)      = coeficiente(337) ! pilambda
      branchings(k,pisigma)       = coeficiente(338) ! pisigma
      branchings(k,pisigma1385m)  = coeficiente(339) ! (pisigma1385)p
      branchings(k,pisigma1385p)  = coeficiente(340) ! (pisigma1385)f
      branchings(k,pilambda1520m) = coeficiente(341) ! (pilambda1520)s
      branchings(k,pilambda1520p) = coeficiente(342) ! (pilambda1520)d
      branchings(k,kbardeltap)    = coeficiente(343) ! (kbardelta)p
      branchings(k,kbar1starn)    = coeficiente(344) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(345) ! (kbar3starn)p
      branchings(k,dummysigma)    = coeficiente(346) ! dummysigma
      return
      end subroutine resonancesigma_p13_1
!     -----------------------------------------------------------------------
      subroutine resonancesigma_p13_2(k) ! sigma(1940)p13
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  1
      isresonancia(k)             =  1.d0
      masasresonancias(k)         = coeficiente(347) 
      branchings(k,kbarn)         = coeficiente(348) ! kbarn
      branchings(k,pilambda)      = coeficiente(349) ! pilambda
      branchings(k,pisigma)       = coeficiente(350) ! pisigma
      branchings(k,pisigma1385m)  = coeficiente(351) ! (pisigma1385)p
      branchings(k,pisigma1385p)  = coeficiente(352) ! (pisigma1385)f
      branchings(k,pilambda1520m) = coeficiente(353) ! (pilambda1520)s
      branchings(k,pilambda1520p) = coeficiente(354) ! (pilambda1520)d
      branchings(k,kbardeltap)    = coeficiente(355) ! (kbardelta)p
      branchings(k,kbar1starn)    = coeficiente(356) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(357) ! (kbar3starn)p
      branchings(k,dummysigma)    = coeficiente(358) ! dummysigma
      return
      end subroutine resonancesigma_p13_2
!     -----------------------------------------------------------------------
      subroutine resonancesigma_p13_3(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  1
      isresonancia(k)             = -1.d0
      masasresonancias(k)         = coeficiente(359)
      branchings(k,kbarn)         = coeficiente(360) ! kbarn
      branchings(k,pilambda)      = coeficiente(361) ! pilambda
      branchings(k,pisigma)       = coeficiente(362) ! pisigma
      branchings(k,pisigma1385m)  = coeficiente(363) ! (pisigma1385)p
      branchings(k,pisigma1385p)  = coeficiente(364) ! (pisigma1385)f
      branchings(k,pilambda1520m) = coeficiente(365) ! (pilambda1520)s
      branchings(k,pilambda1520p) = coeficiente(366) ! (pilambda1520)d
      branchings(k,kbardeltap)    = coeficiente(367) ! (kbardelta)p
      branchings(k,kbar1starn)    = coeficiente(368) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(369) ! (kbar3starn)p
      branchings(k,dummysigma)    = coeficiente(370) ! dummysigma
      return
      end subroutine resonancesigma_p13_3
!     -----------------------------------------------------------------------
      subroutine resonancesigma_p13_4(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  1
      isresonancia(k)             = -1.d0
      masasresonancias(k)         = coeficiente(371)
      branchings(k,kbarn)         = coeficiente(372) ! kbarn
      branchings(k,pilambda)      = coeficiente(373) ! pilambda
      branchings(k,pisigma)       = coeficiente(374) ! pisigma
      branchings(k,pisigma1385m)  = coeficiente(375) ! (pisigma1385)p
      branchings(k,pisigma1385p)  = coeficiente(376) ! (pisigma1385)f
      branchings(k,pilambda1520m) = coeficiente(377) ! (pilambda1520)s
      branchings(k,pilambda1520p) = coeficiente(378) ! (pilambda1520)d
      branchings(k,kbardeltap)    = coeficiente(379) ! (kbardelta)p
      branchings(k,kbar1starn)    = coeficiente(380) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(381) ! (kbar3starn)p
      branchings(k,dummysigma)    = coeficiente(382) ! dummysigma
      return
      end subroutine resonancesigma_p13_4
!     -----------------------------------------------------------------------
      subroutine resonancesigma_p13_5(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  1
      isresonancia(k)             = -1.d0
      masasresonancias(k)         = coeficiente(383)
      branchings(k,kbarn)         = coeficiente(384) ! kbarn
      branchings(k,pilambda)      = coeficiente(385) ! pilambda
      branchings(k,pisigma)       = coeficiente(386) ! pisigma
      branchings(k,pisigma1385m)  = coeficiente(387) ! (pisigma1385)p
      branchings(k,pisigma1385p)  = coeficiente(388) ! (pisigma1385)f
      branchings(k,pilambda1520m) = coeficiente(389) ! (pilambda1520)s
      branchings(k,pilambda1520p) = coeficiente(390) ! (pilambda1520)d
      branchings(k,kbardeltap)    = coeficiente(391) ! (kbardelta)p
      branchings(k,kbar1starn)    = coeficiente(392) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(393) ! (kbar3starn)p
      branchings(k,dummysigma)    = coeficiente(394) ! dummysigma
      return
      end subroutine resonancesigma_p13_5
!     -----------------------------------------------------------------------
      subroutine resonancesigma_p13_6(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  1
      isresonancia(k)             = -1.d0
      masasresonancias(k)         = coeficiente(395)
      branchings(k,kbarn)         = coeficiente(396) ! kbarn
      branchings(k,pilambda)      = coeficiente(397) ! pilambda
      branchings(k,pisigma)       = coeficiente(398) ! pisigma
      branchings(k,pisigma1385m)  = coeficiente(399) ! (pisigma1385)p
      branchings(k,pisigma1385p)  = coeficiente(400) ! (pisigma1385)f
      branchings(k,pilambda1520m) = coeficiente(401) ! (pilambda1520)s
      branchings(k,pilambda1520p) = coeficiente(402) ! (pilambda1520)d
      branchings(k,kbardeltap)    = coeficiente(403) ! (kbardelta)p
      branchings(k,kbar1starn)    = coeficiente(404) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(405) ! (kbar3starn)p
      branchings(k,dummysigma)    = coeficiente(406) ! dummysigma
      return
      end subroutine resonancesigma_p13_6
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     d13 resonances 
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine resonancesigma_d13_1(k) ! sigma(1670)d13
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  2
      isresonancia(k)             =  1.d0
      masasresonancias(k)         = coeficiente(407) 
      branchings(k,kbarn)         = coeficiente(408) ! kbarn
      branchings(k,pilambda)      = coeficiente(409) ! pilambda
      branchings(k,pisigma)       = coeficiente(410) ! pisigma
      branchings(k,pisigma1385m)  = coeficiente(411) ! (pisigma1385)s
      branchings(k,pisigma1385p)  = coeficiente(412) ! (pisigma1385)d
      branchings(k,pilambda1520m) = coeficiente(413) ! (pilambda1520)p
      branchings(k,pilambda1520p) = coeficiente(414) ! (pilambda1520)f
      branchings(k,kbardeltam)    = coeficiente(415) ! (kbardelta)s
      branchings(k,kbardeltap)    = coeficiente(416) ! (kbardelta)d
      branchings(k,kbar1starn)    = coeficiente(417) ! kbar1starn
      branchings(k,kbar3starnm)   = coeficiente(418) ! (kbar3starn)s
      branchings(k,kbar3starnp)   = coeficiente(419) ! (kbar3starn)d
      branchings(k,dummysigma)    = coeficiente(420) ! dummysigma
      return
      end subroutine resonancesigma_d13_1
!     -----------------------------------------------------------------------
      subroutine resonancesigma_d13_2(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  2
      isresonancia(k)             = -1.d0
      masasresonancias(k)         = coeficiente(421)
      branchings(k,kbarn)         = coeficiente(422) ! kbarn
      branchings(k,pilambda)      = coeficiente(423) ! pilambda
      branchings(k,pisigma)       = coeficiente(424) ! pisigma
      branchings(k,pisigma1385m)  = coeficiente(425) ! (pisigma1385)s
      branchings(k,pisigma1385p)  = coeficiente(426) ! (pisigma1385)d
      branchings(k,pilambda1520m) = coeficiente(427) ! (pilambda1520)p
      branchings(k,pilambda1520p) = coeficiente(428) ! (pilambda1520)f
      branchings(k,kbardeltam)    = coeficiente(429) ! (kbardelta)s
      branchings(k,kbardeltap)    = coeficiente(430) ! (kbardelta)d
      branchings(k,kbar1starn)    = coeficiente(431) ! kbar1starn
      branchings(k,kbar3starnm)   = coeficiente(432) ! (kbar3starn)s
      branchings(k,kbar3starnp)   = coeficiente(433) ! (kbar3starn)d
      branchings(k,dummysigma)    = coeficiente(434) ! dummysigma
      return
      end subroutine resonancesigma_d13_2
!     -----------------------------------------------------------------------
      subroutine resonancesigma_d13_3(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  2
      isresonancia(k)             = -1.d0
      masasresonancias(k)         = coeficiente(435)
      branchings(k,kbarn)         = coeficiente(436) ! kbarn
      branchings(k,pilambda)      = coeficiente(437) ! pilambda
      branchings(k,pisigma)       = coeficiente(438) ! pisigma
      branchings(k,pisigma1385m)  = coeficiente(439) ! (pisigma1385)s
      branchings(k,pisigma1385p)  = coeficiente(440) ! (pisigma1385)d
      branchings(k,pilambda1520m) = coeficiente(441) ! (pilambda1520)p
      branchings(k,pilambda1520p) = coeficiente(442) ! (pilambda1520)f
      branchings(k,kbardeltam)    = coeficiente(443) ! (kbardelta)s
      branchings(k,kbardeltap)    = coeficiente(444) ! (kbardelta)d
      branchings(k,kbar1starn)    = coeficiente(445) ! kbar1starn
      branchings(k,kbar3starnm)   = coeficiente(446) ! (kbar3starn)s
      branchings(k,kbar3starnp)   = coeficiente(447) ! (kbar3starn)d
      branchings(k,dummysigma)    = coeficiente(448) ! dummysigma
      return
      end subroutine resonancesigma_d13_3
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     d15 resonances 
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine resonancesigma_d15_1(k) ! sigma(1775)d15
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  2
      isresonancia(k)             =  1.d0
      masasresonancias(k)         = coeficiente(449) 
      branchings(k,kbarn)         = coeficiente(450) ! kbarn
      branchings(k,pilambda)      = coeficiente(451) ! pilambda
      branchings(k,pisigma)       = coeficiente(452) ! pisigma
      branchings(k,pisigma1385m)  = coeficiente(453) ! (pisigma1385)d
      branchings(k,pisigma1385p)  = coeficiente(454) ! (pisigma1385)g
      branchings(k,pilambda1520m) = coeficiente(455) ! (pilambda1520)p
      branchings(k,pilambda1520p) = coeficiente(456) ! (pilambda1520)f
      branchings(k,kbardeltap)    = coeficiente(457) ! (kbardelta)d
      branchings(k,kbar1starn)    = coeficiente(458) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(459) ! (kbar3starn)d
      branchings(k,dummysigma)    = coeficiente(460) ! dummysigma
      return
      end subroutine resonancesigma_d15_1
!     -----------------------------------------------------------------------
      subroutine resonancesigma_d15_2(k) ! sigma(2200)d15
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  2
      isresonancia(k)             =  1.d0
      masasresonancias(k)         = coeficiente(461) 
      branchings(k,kbarn)         = coeficiente(462) ! kbarn
      branchings(k,pilambda)      = coeficiente(463) ! pilambda
      branchings(k,pisigma)       = coeficiente(464) ! pisigma
      branchings(k,pisigma1385m)  = coeficiente(465) ! (pisigma1385)d
      branchings(k,pisigma1385p)  = coeficiente(466) ! (pisigma1385)g
      branchings(k,pilambda1520m) = coeficiente(467) ! (pilambda1520)p
      branchings(k,pilambda1520p) = coeficiente(468) ! (pilambda1520)f
      branchings(k,kbardeltap)    = coeficiente(469) ! (kbardelta)d
      branchings(k,kbar1starn)    = coeficiente(470) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(471) ! (kbar3starn)d
      branchings(k,dummysigma)    = coeficiente(472) ! dummysigma
      return
      end subroutine resonancesigma_d15_2
!     -----------------------------------------------------------------------
      subroutine resonancesigma_d15_3(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  2
      isresonancia(k)             = -1.d0
      masasresonancias(k)         = coeficiente(473)
      branchings(k,kbarn)         = coeficiente(474) ! kbarn
      branchings(k,pilambda)      = coeficiente(475) ! pilambda
      branchings(k,pisigma)       = coeficiente(476) ! pisigma
      branchings(k,pisigma1385m)  = coeficiente(477) ! (pisigma1385)d
      branchings(k,pisigma1385p)  = coeficiente(478) ! (pisigma1385)g
      branchings(k,pilambda1520m) = coeficiente(479) ! (pilambda1520)p
      branchings(k,pilambda1520p) = coeficiente(480) ! (pilambda1520)f
      branchings(k,kbardeltap)    = coeficiente(481) ! (kbardelta)d
      branchings(k,kbar1starn)    = coeficiente(482) ! kbar1starn
      branchings(k,kbar3starnp)   = coeficiente(483) ! (kbar3starn)d
      branchings(k,dummysigma)    = coeficiente(484) ! dummysigma
      return
      end subroutine resonancesigma_d15_3
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     f15 resonances 
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine resonancesigma_f15_1(k) ! sigma(1915)f15
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  3
      isresonancia(k)             =  1.d0
      masasresonancias(k)         = coeficiente(485) 
      branchings(k,kbarn)         = coeficiente(486) ! kbarn
      branchings(k,pilambda)      = coeficiente(487) ! pilambda
      branchings(k,pisigma)       = coeficiente(488) ! pisigma
      branchings(k,pisigma1385m)  = coeficiente(489) ! (pisigma1385)p
      branchings(k,pisigma1385p)  = coeficiente(490) ! (pisigma1385)f
      branchings(k,pilambda1520m) = coeficiente(491) ! (pilambda1520)d
      branchings(k,pilambda1520p) = coeficiente(492) ! (pilambda1520)g
      branchings(k,kbardeltap)    = coeficiente(493) ! (kbardelta)p
      branchings(k,kbar1starn)    = coeficiente(494) ! kbar1starn
      branchings(k,kbar3starnm)   = coeficiente(495) ! (kbar3starn)p
      branchings(k,kbar3starnp)   = coeficiente(496) ! (kbar3starn)f
      branchings(k,dummysigma)    = coeficiente(497) ! dummysigma
      return
      end subroutine resonancesigma_f15_1
!     -----------------------------------------------------------------------
      subroutine resonancesigma_f15_2(k) ! sigma(2100)f15
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  3
      isresonancia(k)             =  1.d0
      masasresonancias(k)         = coeficiente(498) 
      branchings(k,kbarn)         = coeficiente(499) ! kbarn
      branchings(k,pilambda)      = coeficiente(500) ! pilambda
      branchings(k,pisigma)       = coeficiente(501) ! pisigma
      branchings(k,pisigma1385m)  = coeficiente(502) ! (pisigma1385)p
      branchings(k,pisigma1385p)  = coeficiente(503) ! (pisigma1385)f
      branchings(k,pilambda1520m) = coeficiente(504) ! (pilambda1520)d
      branchings(k,pilambda1520p) = coeficiente(505) ! (pilambda1520)g
      branchings(k,kbardeltap)    = coeficiente(506) ! (kbardelta)p
      branchings(k,kbar1starn)    = coeficiente(507) ! kbar1starn
      branchings(k,kbar3starnm)   = coeficiente(508) ! (kbar3starn)p
      branchings(k,kbar3starnp)   = coeficiente(509) ! (kbar3starn)f
      branchings(k,dummysigma)    = coeficiente(510) ! dummysigma
      return
      end subroutine resonancesigma_f15_2
!     -----------------------------------------------------------------------
      subroutine resonancesigma_f15_3(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  3
      isresonancia(k)             = -1.d0
      masasresonancias(k)         = coeficiente(511)
      branchings(k,kbarn)         = coeficiente(512) ! kbarn
      branchings(k,pilambda)      = coeficiente(513) ! pilambda
      branchings(k,pisigma)       = coeficiente(514) ! pisigma
      branchings(k,pisigma1385m)  = coeficiente(515) ! (pisigma1385)p
      branchings(k,pisigma1385p)  = coeficiente(516) ! (pisigma1385)f
      branchings(k,pilambda1520m) = coeficiente(517) ! (pilambda1520)d
      branchings(k,pilambda1520p) = coeficiente(518) ! (pilambda1520)g
      branchings(k,kbardeltap)    = coeficiente(519) ! (kbardelta)p
      branchings(k,kbar1starn)    = coeficiente(520) ! kbar1starn
      branchings(k,kbar3starnm)   = coeficiente(521) ! (kbar3starn)p
      branchings(k,kbar3starnp)   = coeficiente(522) ! (kbar3starn)f
      branchings(k,dummysigma)    = coeficiente(523) ! dummysigma
      return
      end subroutine resonancesigma_f15_3
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     f17 resonances 
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine resonancesigma_f17_1(k) ! sigma(2030)f17
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  3
      isresonancia(k)             =  1.d0
      masasresonancias(k)         = coeficiente(524) 
      branchings(k,kbarn)         = coeficiente(525) ! kbarn
      branchings(k,pilambda)      = coeficiente(526) ! pilambda
      branchings(k,pisigma)       = coeficiente(527) ! pisigma
      branchings(k,pisigma1385p)  = coeficiente(528) ! (pisigma1385)f
      branchings(k,pilambda1520m) = coeficiente(529) ! (pilambda1520)d
      branchings(k,pilambda1520p) = coeficiente(530) ! (pilambda1520)g
      branchings(k,kbardeltap)    = coeficiente(531) ! (kbardelta)f
      branchings(k,kbar1starn)    = coeficiente(532) ! kbar1starn
      branchings(k,kbar3starnm)   = coeficiente(533) ! (kbar3starn)f
      branchings(k,kbar3starnp)   = coeficiente(534) ! (kbar3starn)h
      branchings(k,dummysigma)    = coeficiente(535) ! dummysigma
      return
      end subroutine resonancesigma_f17_1
!     -----------------------------------------------------------------------
      subroutine resonancesigma_f17_2(k)
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)         =  3
      isresonancia(k)              = -1.d0
      masasresonancias(k)          = coeficiente(536)
      branchings(k,kbarn)          = coeficiente(537) ! kbarn
      branchings(k,pilambda)       = coeficiente(538) ! pilambda
      branchings(k,pisigma)        = coeficiente(539) ! pisigma
      branchings(k,pisigma1385p)   = coeficiente(540) ! (pisigma1385)g
      branchings(k,pilambda1520m)  = coeficiente(541) ! (pilambda1520)d
      branchings(k,pilambda1520p)  = coeficiente(542) ! (pilambda1520)g
      branchings(k,kbardeltap)     = coeficiente(543) ! (kbardelta)f
      branchings(k,kbar1starn)     = coeficiente(544) ! kbar1starn
      branchings(k,kbar3starnm)    = coeficiente(545) ! (kbar3starn)f
      branchings(k,kbar3starnp)    = coeficiente(546) ! (kbar3starn)h
      branchings(k,dummysigma)     = coeficiente(547) ! dummysigma
      return
      end subroutine resonancesigma_f17_2
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     g17 resonances 
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine resonancesigma_g17_1(k) ! sigma(2250)g17
      use coeficientes
      use channelassignment
      use matrixstorage
      implicit none
      integer :: k
      momangresonancias(k)        =  4
      isresonancia(k)             =  1.d0
      masasresonancias(k)         = coeficiente(548) 
      branchings(k,kbarn)         = coeficiente(549) ! kbarn
      branchings(k,pilambda)      = coeficiente(550) ! pilambda
      branchings(k,pisigma)       = coeficiente(551) ! pisigma
      branchings(k,pilambda1520m) = coeficiente(552) ! (pilambda1520)f
      branchings(k,pilambda1520p) = coeficiente(553) ! (pilambda1520)h
      branchings(k,kbardeltam)    = coeficiente(554) ! (kbardelta)g
      branchings(k,kbardeltap)    = coeficiente(555) ! (kbardelta)f
      branchings(k,kbar3starnm)   = coeficiente(556) ! (kbar3starn)f
      branchings(k,kbar3starnp)   = coeficiente(557) ! (kbar3starn)h
      branchings(k,dummysigma)    = coeficiente(558) ! dummysigma
      return
      end subroutine resonancesigma_g17_1
!     -----------------------------------------------------------------------
!     Phi calculator
!     -----------------------------------------------------------------------
      subroutine phicalculator(l)
      use usefulnumbers
      use superdata
      use matrixstorage
      implicit none
      integer :: l,iloop1,iloop2,iloop3,iloop4
      real (8) :: s,s00
      real (8), external :: phiphase,phinormfunc
      complex (8), external :: phi
      do iloop2 = 1,matrixsize
         do iloop3 = 1,indice_maximo
            s = svector(iloop3)
            phiphasestorage(iloop2,iloop3) = phiphase(l,iloop2,s)
            phistorage(iloop2,iloop3) = phi(l,iloop2,s)
         enddo
         do iloop4 = indice_maximo+1,max_data_points
            phistorage(iloop2,iloop4) = dcmplx(0.d0,0.d0)
            phiphasestorage(iloop2,iloop4) = dcmplx(0.d0,0.d0)
         enddo
      enddo
!     -----------------------------------------------------------------------
!     background widths normalization and initialization
!     -----------------------------------------------------------------------
      s00 = 4.d0
      do iloop1=1,nresonances
         do iloop2=1,matrixsize
            phimres(iloop1,iloop2) = phinormfunc(l,iloop2,s00)
         enddo
         do iloop2=matrixsize+1,ndecays
            phimres(iloop1,iloop2) = 0.d0
         enddo
      enddo
      return
      end subroutine phicalculator
!     -----------------------------------------------------------------------
!     s matrix
!     -----------------------------------------------------------------------
      subroutine tpwcalculation(s,indice,pwchoice,a1,a2,a3,a4,a5,a6)
      use canalesop
      use matrixoperations
      use matrixstorage
      implicit none
      character (3) :: pwchoice
      integer :: i,j
      integer :: indice,l,a1,a2,a3,a4,a5,a6
      real (8) :: s,numero1,numero2
      integer, dimension (nres) :: cres
      complex (8), dimension (matrixsize,matrixsize) :: phasemat,rhomat,
     c     tresmat
!     -----------------------------------------------------------------------
      call resonancesparametersinput
!     -----------------------------------------------------------------------
      select case (pwchoice)
      case ('s01')
         cres(1) = cs01(1)
         cres(2) = cs01(2)
         cres(3) = cs01(3)
         cres(4) = cs01(4)
         cres(5) = cs01(5)
         cres(6) = cs01(6)
      case ('p01')
         cres(1) = cp01(1)
         cres(2) = cp01(2)
         cres(3) = cp01(3)
         cres(4) = cp01(4)
         cres(5) = cp01(5)
         cres(6) = cp01(6)
      case ('p03')
         cres(1) = cp03(1)
         cres(2) = cp03(2)
         cres(3) = cp03(3)
         cres(4) = cp03(4)
         cres(5) = cp03(5)
         cres(6) = cp03(6)
      case ('d03')
         cres(1) = cd03(1)
         cres(2) = cd03(2)
         cres(3) = cd03(3)
         cres(4) = cd03(4)
         cres(5) = cd03(5)
         cres(6) = cd03(6)
      case ('d05')
         cres(1) = cd05(1)
         cres(2) = cd05(2)
         cres(3) = cd05(3)
         cres(4) = cd05(4)
         cres(5) = cd05(5)
         cres(6) = cd05(6)
      case ('f05')
         cres(1) = cf05(1)
         cres(2) = cf05(2)
         cres(3) = cf05(3)
         cres(4) = cf05(4)
         cres(5) = cf05(5)
         cres(6) = cf05(6)
      case ('f07')
         cres(1) = cf07(1)
         cres(2) = cf07(2)
         cres(3) = cf07(3)
         cres(4) = cf07(4)
         cres(5) = cf07(5)
         cres(6) = cf07(6)
      case ('g07')
         cres(1) = cg07(1)
         cres(2) = cg07(2)
         cres(3) = cg07(3)
         cres(4) = cg07(4)
         cres(5) = cg07(5)
         cres(6) = cg07(6)
!     -----------------------------------------------------------------------
      case ('s11')
         cres(1) = cs11(1)
         cres(2) = cs11(2)
         cres(3) = cs11(3)
         cres(4) = cs11(4)
         cres(5) = cs11(5)
         cres(6) = cs11(6)
      case ('p11')
         cres(1) = cp11(1)
         cres(2) = cp11(2)
         cres(3) = cp11(3)
         cres(4) = cp11(4)
         cres(5) = cp11(5)
         cres(6) = cp11(6)
      case ('p13')
         cres(1) = cp13(1)
         cres(2) = cp13(2)
         cres(3) = cp13(3)
         cres(4) = cp13(4)
         cres(5) = cp13(5)
         cres(6) = cp13(6)
      case ('d13')
         cres(1) = cd13(1)
         cres(2) = cd13(2)
         cres(3) = cd13(3)
         cres(4) = cd13(4)
         cres(5) = cd13(5)
         cres(6) = cd13(6)
      case ('d15')
         cres(1) = cd15(1)
         cres(2) = cd15(2)
         cres(3) = cd15(3)
         cres(4) = cd15(4)
         cres(5) = cd15(5)
         cres(6) = cd15(6)
      case ('f15')
         cres(1) = cf15(1)
         cres(2) = cf15(2)
         cres(3) = cf15(3)
         cres(4) = cf15(4)
         cres(5) = cf15(5)
         cres(6) = cf15(6)
      case ('f17')
         cres(1) = cf17(1)
         cres(2) = cf17(2)
         cres(3) = cf17(3)
         cres(4) = cf17(4)
         cres(5) = cf17(5)
         cres(6) = cf17(6)
      case ('g17')
         cres(1) = cg17(1)
         cres(2) = cg17(2)
         cres(3) = cg17(3)
         cres(4) = cg17(4)
         cres(5) = cg17(5)
         cres(6) = cg17(6)
      end select
!     -----------------------------------------------------------------------
      a1 = cres(1)
      a2 = cres(2)
      a3 = cres(3)
      a4 = cres(4)
      a5 = cres(5)
      a6 = cres(6)
      call phiresnorm(a1,a2,a3,a4,a5,a6)
      l = momangresonancias(a1)
      phasemat = rhomatrixphase(indice)
      rhomat   = rhomatrix(indice)
      tresmat  = tr6fast(s,indice,a1,a2,a3,a4,a5,a6)
      scatmat  = smatrix(id,phasemat,tresmat)
!      print*,pwchoice,l,a1,a2,a3,a4,a5,a6
!      print*,indice,s,tracec(matmul(scatmat,dagger(scatmat)))-matrixsize
      return
      end subroutine tpwcalculation
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      subroutine amatrix(s,indice,i1,i2,i3,i4,i5,i6,
     c     a11,a22,a33,a44,a55,a66,a12,a13,a14,a15,a16,
     c     a23,a24,a25,a26,a34,a35,a36,a45,a46,a56,d6)
      use usefulnumbers
      use matrixstorage
      implicit none
      integer  :: indice,j,k,i1,i2,i3,i4,i5,i6
      real (8), parameter :: units1 = 1000.d0, units2 = dsqrt(1000.d0)
      real (8) :: s, mr1,mr2,mr3,mr4,mr5,mr6,isres
      complex (8) :: gr1,gr2,gr3,gr4,gr5,gr6,anch1,anch2,anch3,anch4,
     c     anch5,anch6
      complex (8) :: e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,
     c     e45,e46,e56
      complex (8) :: a11,a22,a33,a44,a55,a66,a12,a13,a14,a15,a16,a23,
     c     a24,a25,a26,a34,a35,a36,a45,a46,a56,t1,t2,t3,t4,t5,t6,d6
      complex (8), external :: breitwigner,epsilon2dres
      complex (8), external :: a12function,a13function,a14function,
     c     a15function,a16function,a23function,a24function,a25function,
     c     a26function,a34function,a35function,a36function,a45function,
     c     a46function,a56function
      complex (8), external :: afunction,determinante
!     resonance 1
      gr1     = epsilon2dres(i1,i1,indice)
      anch1   = gr1/units1
      mr1     = masasresonancias(i1)
      isres   = isresonancia(i1)
      t1      = breitwigner(s,mr1,anch1,isres)
!     resonance 2
      if (i2.ne.0) then
         mr2   = masasresonancias(i2)
         gr2   = epsilon2dres(i2,i2,indice)
         anch2 = gr2/units1
         e12   = epsilon2dres(i1,i2,indice)
         isres = isresonancia(i2)
         t2    = breitwigner(s,mr2,anch2,isres)
      else
         mr2   = 0.d0
         gr2   = xzero
         anch2 = xzero
         e12   = xzero
         t2    = xzero
      endif
!     resonance 3
      if (i3.ne.0) then
         mr3   = masasresonancias(i3)
         gr3   = epsilon2dres(i3,i3,indice)
         anch3 = gr3/units1
         e13   = epsilon2dres(i1,i3,indice)
         e23   = epsilon2dres(i2,i3,indice)
         isres = isresonancia(i3)
         t3    = breitwigner(s,mr3,anch3,isres)
      else
         mr3   = 0.d0
         gr3   = xzero
         anch3 = xzero
         e13   = xzero
         e23   = xzero
         t3    = xzero
      endif
!     resonance 4
      if (i4.ne.0) then
         mr4   = masasresonancias(i4)
         gr4   = epsilon2dres(i4,i4,indice)
         anch4 = gr4/units1
         e14   = epsilon2dres(i1,i4,indice)
         e24   = epsilon2dres(i2,i4,indice)
         e34   = epsilon2dres(i3,i4,indice)
         isres = isresonancia(i4)
         t4    = breitwigner(s,mr4,anch4,isres)
      else
         mr4   = 0.d0
         gr4   = xzero
         anch4 = xzero
         e14   = xzero
         e24   = xzero
         e34   = xzero
         t4    = xzero
      endif
!     resonance 5
      if (i5.ne.0) then
         mr5   = masasresonancias(i5)
         gr5   = epsilon2dres(i5,i5,indice)
         anch5 = gr5/units1
         e15   = epsilon2dres(i1,i5,indice)
         e25   = epsilon2dres(i2,i5,indice)
         e35   = epsilon2dres(i3,i5,indice)
         e45   = epsilon2dres(i4,i5,indice)
         isres = isresonancia(i5)
         t5    = breitwigner(s,mr5,anch5,isres)
      else
         mr5   = 0.d0
         gr5   = xzero
         anch5 = xzero
         e15   = xzero
         e25   = xzero
         e35   = xzero
         e45   = xzero
         t5    = xzero
      endif
!     resonance 6
      if (i6.ne.0) then
         mr6   = masasresonancias(i6)
         gr6   = epsilon2dres(i6,i6,indice)
         anch6 = gr6/units1
         e16   = epsilon2dres(i1,i6,indice)
         e26   = epsilon2dres(i2,i6,indice)
         e36   = epsilon2dres(i3,i6,indice)
         e46   = epsilon2dres(i4,i6,indice)
         e56   = epsilon2dres(i5,i6,indice)
         isres = isresonancia(i6)
         t6    = breitwigner(s,mr6,anch6,isres)
      else
         mr6   = 0.d0
         gr6   = xzero
         anch6 = xzero
         e16   = xzero
         e26   = xzero
         e36   = xzero
         e46   = xzero
         e56   = xzero
         t6    = xzero
      endif
      d6 = determinante(t1,t2,t3,t4,t5,t6,
     c     e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
      a11 = afunction(t1,t2,t3,t4,t5,t6,
     c   e12,e13,e14,e15,e16,e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
      a22 = xzero
      a12 = xzero
      if (i2.ne.0) then
         a12 = a12function(t1,t2,t3,t4,t5,t6,e12,e13,e14,e15,e16,e23,
     c        e24,e25,e26,e34,e35,e36,e45,e46,e56)
         a22 = afunction(t2,t1,t3,t4,t5,t6,e12,e23,e24,e25,e26,e13,
     c        e14,e15,e16,e34,e35,e36,e45,e46,e56)
      endif
      a33 = xzero
      a13 = xzero
      a23 = xzero
      if (i3.ne.0) then
         a13 = a13function(t1,t2,t3,t4,t5,t6,e12,e13,e14,e15,e16,e23,
     c        e24,e25,e26,e34,e35,e36,e45,e46,e56)
         a23 = a23function(t1,t2,t3,t4,t5,t6,e12,e13,e14,e15,e16,e23,
     c        e24,e25,e26,e34,e35,e36,e45,e46,e56)
         a33 = afunction(t3,t2,t1,t4,t5,t6,e23,e13,e34,e35,e36,e12,
     c        e24,e25,e26,e14,e15,e16,e45,e46,e56)
      endif
      a44 = xzero
      a14 = xzero
      a24 = xzero
      a34 = xzero
      if (i4.ne.0) then
         a14 = a14function(t1,t2,t3,t4,t5,t6,e12,e13,e14,e15,e16,e23,
     c        e24,e25,e26,e34,e35,e36,e45,e46,e56)
         a24 = a24function(t1,t2,t3,t4,t5,t6,e12,e13,e14,e15,e16,e23,
     c        e24,e25,e26,e34,e35,e36,e45,e46,e56)
         a34 = a34function(t1,t2,t3,t4,t5,t6,e12,e13,e14,e15,e16,e23,
     c        e24,e25,e26,e34,e35,e36,e45,e46,e56)
         a44 = afunction(t4,t2,t3,t1,t5,t6,e24,e34,e14,e45,e46,e23,
     c        e12,e25,e26,e13,e35,e36,e15,e16,e56)
      endif
      a55 = xzero
      a15 = xzero
      a25 = xzero
      a35 = xzero
      a45 = xzero
      if (i5.ne.0) then
         a15 = a15function(t1,t2,t3,t4,t5,t6,e12,e13,e14,e15,e16,
     c        e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
         a25 = a25function(t1,t2,t3,t4,t5,t6,e12,e13,e14,e15,e16,
     c        e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
         a35 = a35function(t1,t2,t3,t4,t5,t6,e12,e13,e14,e15,e16,
     c        e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
         a45 = a45function(t1,t2,t3,t4,t5,t6,e12,e13,e14,e15,e16,
     c        e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
         a55 = afunction(t5,t2,t3,t4,t1,t6,e25,e35,e45,e15,e56,
     c        e23,e24,e12,e26,e34,e13,e36,e14,e46,e16)
      endif
      a66 = xzero
      a16 = xzero
      a26 = xzero
      a36 = xzero
      a46 = xzero
      a56 = xzero
      if (i6.ne.0) then
         a16 = a16function(t1,t2,t3,t4,t5,t6,e12,e13,e14,e15,e16,
     c        e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
         a26 = a26function(t1,t2,t3,t4,t5,t6,e12,e13,e14,e15,e16,
     c        e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
         a36 = a36function(t1,t2,t3,t4,t5,t6,e12,e13,e14,e15,e16,
     c        e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
         a46 = a46function(t1,t2,t3,t4,t5,t6,e12,e13,e14,e15,e16,
     c        e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
         a56 = a56function(t1,t2,t3,t4,t5,t6,e12,e13,e14,e15,e16,
     c        e23,e24,e25,e26,e34,e35,e36,e45,e46,e56)
         a66 = afunction(t6,t3,t2,t4,t5,t1,e36,e26,e46,e56,e16,
     c        e23,e34,e35,e13,e24,e25,e12,e45,e14,e15)
      endif
      return
      end subroutine amatrix
!     -----------------------------------------------------------------------
!     Resonant T matrix
!     -----------------------------------------------------------------------
      subroutine xxmatrix(j,k,i1,i2,i3,i4,i5,i6,
     c           x11,x22,x33,x44,x55,x66,x12,x13,x14,x15,x16,
     c           x23,x24,x25,x26,x34,x35,x36,x45,x46,x56)
      use usefulnumbers
      use matrixstorage
      implicit none
      integer  :: j,k,i1,i2,i3,i4,i5,i6
      real (8), parameter :: units1 = 1000.d0
      real (8), parameter :: units2 = units1*units1
      real (8) :: x1j,x1k,x2j,x2k,x3j,x3k,x4j,x4k,x5j,x5k,x6j,x6k
      complex (8) :: x11,x22,x33,x44,x55,x66,x12,x13,x14,x15,x16,x23,
     c     x24,x25,x26,x34,x35,x36,x45,x46,x56
!     resonance 1
      x1j     = branchings(i1,j)/dsqrt(phimres(i1,j))
      x1k     = branchings(i1,k)/dsqrt(phimres(i1,k))
!     resonance 2
      if (i2.ne.0) then
         x2j   = branchings(i2,j)/dsqrt(phimres(i2,j))
         x2k   = branchings(i2,k)/dsqrt(phimres(i2,k))
      else
         x2j   = 0.d0
         x2k   = 0.d0
      endif
!     resonance 3
      if (i3.ne.0) then
         x3j   = branchings(i3,j)/dsqrt(phimres(i3,j))
         x3k   = branchings(i3,k)/dsqrt(phimres(i3,k))
      else
         x3j   = 0.d0
         x3k   = 0.d0
      endif
!     resonance 4
      if (i4.ne.0) then
         x4j   = branchings(i4,j)/dsqrt(phimres(i4,j))
         x4k   = branchings(i4,k)/dsqrt(phimres(i4,k))
      else
         x4j   = 0.d0
         x4k   = 0.d0
      endif
!     resonance 5
      if (i5.ne.0) then
         x5j   = branchings(i5,j)/dsqrt(phimres(i5,j))
         x5k   = branchings(i5,k)/dsqrt(phimres(i5,k))
      else
         x5j   = 0.d0
         x5k   = 0.d0
      endif
!     resonance 6
      if (i6.ne.0) then
         x6j   = branchings(i6,j)/dsqrt(phimres(i6,j))
         x6k   = branchings(i6,k)/dsqrt(phimres(i6,k))
      else
         x6j   = 0.d0
         x6k   = 0.d0
      endif
!     output
      x11 = xr*x1j*x1k
      x22 = xr*x2j*x2k
      x33 = xr*x3j*x3k
      x44 = xr*x4j*x4k
      x55 = xr*x5j*x5k
      x66 = xr*x6j*x6k
      x12 = xr*( x1j*x2k + x2j*x1k )/units1
      x13 = xr*( x1j*x3k + x3j*x1k )/units1
      x14 = xr*( x1j*x4k + x4j*x1k )/units1
      x15 = xr*( x1j*x5k + x5j*x1k )/units1
      x16 = xr*( x1j*x6k + x6j*x1k )/units1
      x23 = xr*( x2j*x3k + x3j*x2k )/units1
      x24 = xr*( x2j*x4k + x4j*x2k )/units1
      x25 = xr*( x2j*x5k + x5j*x2k )/units1
      x26 = xr*( x2j*x6k + x6j*x2k )/units1
      x34 = xr*( x3j*x4k + x4j*x3k )/units1
      x35 = xr*( x3j*x5k + x5j*x3k )/units1
      x36 = xr*( x3j*x6k + x6j*x3k )/units1
      x45 = xr*( x4j*x5k + x5j*x4k )/units1
      x46 = xr*( x4j*x6k + x6j*x4k )/units1
      x56 = xr*( x5j*x6k + x6j*x5k )/units1
      return
      end subroutine xxmatrix
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      complex (8) function breitwigner(s,m1,g1,isres) 
     c     result (output)
      use usefulnumbers
      implicit none
      real (8) :: s,m1,isres
      complex (8) :: w,g1,salida
      salida = xr*m1/(m1*m1 - isres*s - xi*m1*g1)
      output = salida
      return
      end function breitwigner
!     -----------------------------------------------------------------------
!     epsilon_(ab)
!     -----------------------------------------------------------------------
!     nmatrix is the size of the t matrix (number of decay channels)
!     i and j fix the resonances
      complex (8) function epsilon2dres(i,j,indice) result (output)
      use usefulnumbers
      use matrixstorage
      implicit none
      integer :: i,j,indice,n
      real (8) :: s,coupling,phinorm
      complex (8) :: xsum
      xsum = dcmplx(0.d0,0.d0)
      do n=1,matrixsize
         coupling = branchings(i,n)*branchings(j,n)
         phinorm  = dsqrt(phimres(i,n)*phimres(j,n))
         xsum     = xsum + phistorage(n,indice)*coupling/phinorm
      enddo
      output = xsum
      return
      end function epsilon2dres
!     ----------------------------------------------------------------------
!     ----------------------------------------------------------------------
!     Rho matrix calculation (complex quantity)
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     k: resonance
!     j: channel
!     s: s mandelstam variable
      complex (8) function phi(l,j,s) result (output)
      use usefulnumbers
      use channelassignment
      use matrixstorage
      implicit none
      integer :: l,j
      real (8) :: s,m1,m2
      complex (8) :: salida
      complex (8), external :: stable_channel
      m1 = channelmass(j)
      m2 = channelmass(ndecays+j)
      output = stable_channel(l,m1,m2,s)
      return
      end function phi
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     l: angular momentum of the resonance
!     m1: mass of decay product 1
!     m2: mass of decay product 2
!     w: sqrt(s mandelstam variable)
      complex (8) function stable_channel(l,m1,m2,sr) result (output)
      use usefulnumbers
      implicit none
      integer :: l
      real (8) :: sr,m1,m2,st,a,a0
      complex (8) :: s,salida
      st = (m1+m2)**2
      a  = fermi2*m1*m2/st
      a0 = dsqrt(m1*m2/st)/w00
      s  = sr + xi*epsilon
      select case(l)
      case (0)
        salida = a0*(s-st)/zsqrt(-s + st)
      case (1)
        salida = a0*a*(s-st)/(xr*dsqrt(a) + a*zsqrt(-s + st))
      case (2)
        salida = a0*(a**2)*(s-st)*((-s + st)**1.5*(2.d0 + (1.d0 
     c    + 3.d0*a*(s - st))/(a*(-s + st))**1.5))
     c    /(2.d0*(1.d0 + a*(s - st))**2)
      case (3)
        salida = a0*(a**3)*(s-st)*((-s + st)**2.5*(-8.d0 
     c    + (3.d0*(1.d0 + (10.d0*a*(s - st))/3.d0 + 5.d0*a**2
     c    *(s - st)**2))/(a*(-s + st))**2.5))
     c    /(8.d0*(1.d0 + a*(s - st))**3)
      case (4)
        salida = a0*(a**4)*(s-st)*((-s + st)**3.5*(16.d0 
     c    + (5.d0*(1.d0 + (21.d0*a*(s - st))/5.d0 + 7.d0*a**2
     c    *(s - st)**2 + 7.d0*a**3*(s - st)**3))/
     c    (a*(-s + st))**3.5))/(16.d0*(1.d0 + a*(s - st))**4)
      case default
        stop 'wrong l, stable_channel'
      end select
      output = -xi*salida
      return
      end function stable_channel
!     ----------------------------------------------------------------------
!     ----------------------------------------------------------------------
!     reonance/bkg normalization zabs(complex quantity)
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
!     k: resonance
!     j: channel
!     s: s mandelstam variable
      real (8) function phinormfunc(l,j,s) result (output)
      use usefulnumbers
      use channelassignment
      use matrixstorage
      implicit none
      integer :: l,j
      real (8) :: s,m1,m2,q2fermi,st
      complex (8) :: q2
      m1 = channelmass(j)
      m2 = channelmass(ndecays+j)
      st = (m1+m2)**2
      q2      = dcmplx(m1*m2*( s - st )/st,epsilon)
      q2fermi = dble(q2*fermi2)
      output  = zabs(xi*zsqrt(-xr*q2)*((q2fermi/(1.d0+q2fermi))**l)/w00)
      return
      end function phinormfunc
!     ----------------------------------------------------------------------
!     ----------------------------------------------------------------------
!     Phase space calculation  (real quantity and above threshold)
!     -----------------------------------------------------------------------
!     -----------------------------------------------------------------------
      real (8) function phiphase(l,j,s) result (output)
      use usefulnumbers
      use channelassignment
      use matrixstorage
      implicit none
      integer :: l,j
      real (8) :: m1,m2,s,st,q2,q2fermi
      m1 = channelmass(j)
      m2 = channelmass(ndecays+j)
      st = (m1+m2)**2
      q2 = m1*m2*( s - st )/st
      if (q2.lt.0.d0) q2 = 0.d0      
      q2fermi = q2*fermi2
      output  = dsqrt(q2)/w00*((q2fermi/(1.d0+q2fermi))**l)
      return
      end function phiphase
!     -----------------------------------------------------------------------
!     identity matrix
!     -----------------------------------------------------------------------
      subroutine identityinit
      use usefulnumbers
      use matrixstorage
      implicit none
      integer :: iloop1,iloop2
      do iloop1=1,matrixsize
         do iloop2=iloop1,matrixsize
            id(iloop1,iloop2) = xzero
            id(iloop2,iloop1) = xzero
         enddo
         id(iloop1,iloop1) = xr
      enddo
      return
      end subroutine identityinit
!     -----------------------------------------------------------------------
