Index: tags/ohl/attic/omega-000.011beta/tools/testbed_old.f95
===================================================================
--- tags/ohl/attic/omega-000.011beta/tools/testbed_old.f95	(revision 8687)
+++ tags/ohl/attic/omega-000.011beta/tools/testbed_old.f95	(revision 8688)
@@ -1,1601 +0,0 @@
-! $Id: testbed_old.f95,v 1.1 2004/04/09 20:11:17 ohl Exp $
-
-module testbed_old
-   use omega_kinds
-   use omega95
-   use kinematics
-   use rambo
-   use omega_parameters
-   use omega_parameters_madgraph
-   implicit none
-   real(kind=omega_prec), parameter, private :: THRESHOLD_SIGMA = 1.0e-21_omega_prec
-   public :: read_parameters
-   public :: compare_sigma
-   public :: compare_squared
-   public :: compare
-   public :: compare_omega
-   public :: compare_omega_sum
-   public :: symmetry_omega
-   public :: check_omega
-   public :: ward_omega
-   public :: ward4
-   public :: compare4_madgraph
-   public :: compare5_madgraph
-   public :: compare6_madgraph
-   public :: compare7_madgraph
-   public :: compare8_madgraph
-   public :: compare_sum4_madgraph
-   public :: compare_sum5_madgraph
-   public :: compare_sum6_madgraph
-   public :: compare_sum7_madgraph
-   public :: compare_sum8_madgraph
-   public :: check4_madgraph
-   public :: check5_madgraph
-   public :: check6_madgraph
-   public :: check7_madgraph
-   public :: check8_madgraph
-contains
- 
-  subroutine read_parameters (roots, n, tolerance, mode)
-    real(kind=single), intent(out) :: roots
-    integer, intent(out) :: n, tolerance
-    character(len=8), intent(out) :: mode
-    real(kind=single) :: me, mw, gw, mz, gz, mt, gt, mh, gh
-    complex(kind=single) :: o1, o2, o3, o4, m1, m2, m3, m4
-    namelist /options/ mode, n, tolerance, roots, me, mw, gw, mz, gz, &
-         mt, gt, mh, gh, o1, o2, o3, o4, m1, m2, m3, m4
-    integer :: ios
-    character (len=128) :: cmd
-    call setup_parameters ()
-    mode = "compare"
-    n = 10
-    tolerance = 10000
-    roots = 800
-    me = mass(11)
-    mz = mass(23)
-    gz = width(23)
-    mw = mass(24)
-    gw = width(24)
-    mt = mass(6)
-    gt = width(6)
-    mh = mass(25)
-    gh = width(25)
-    o1 = fudge_o1
-    o2 = fudge_o2
-    o3 = fudge_o3
-    o4 = fudge_o4
-    m1 = fudge_m1
-    m2 = fudge_m2
-    m3 = fudge_m3
-    m4 = fudge_m4
-    write (unit = *, nml = options)
-    open (unit = 42, status = "scratch")
-    cmds: do
-       read (unit = *, fmt = "(A)", iostat = ios) cmd
-       if (ios /= 0 .or. trim(cmd) == "") exit cmds
-       rewind (unit = 42)
-       write (unit = 42, fmt = "('&options ',A,'/')") cmd
-       rewind (unit = 42)
-       read (unit = 42, nml = options)
-    end do cmds
-    close (unit = 42)
-    write (unit = *, nml = options)
-    mass(11) = me
-    mass(23) = mz
-    width(23) = gz
-    mass(24) = mw
-    width(24) = gw
-    mass(6) = mt
-    width(6) = gt
-    mass(25) = mh
-    width(25) = gh
-    fudge_o1 = o1
-    fudge_o2 = o2
-    fudge_o3 = o3
-    fudge_o4 = o4
-    fudge_m1 = m1
-    fudge_m2 = m2
-    fudge_m3 = m3
-    fudge_m4 = m4
-  end subroutine read_parameters
-
-  subroutine compare_sigma (n1, s1, n2, s2, threshold, tolerance)
-    real(kind=omega_prec), intent(in) :: s1, s2
-    character(len=*), intent(in) :: n1, n2
-    real(kind=omega_prec), intent(in) :: threshold
-    integer, intent(in), optional :: tolerance
-    real(kind=omega_prec) :: ds, tolerance_local
-    if (present (tolerance)) then
-       tolerance_local = tolerance * epsilon (tolerance_local)
-    else
-       tolerance_local = 1000 * epsilon (tolerance_local)
-    end if
-    if (max (s1, s2) >= threshold) then
-       ds = abs (s1 - s2) / (s2 + s2)
-       if (ds > tolerance_local) then
-          write (unit = *, fmt = &
-               "(2(2X,A1,': ',E17.10), 2X,A2,1X,E8.2, 1X,A1,1X,E8.2, A4)") &
-               n1, s1, n2, s2, "d=", ds, "=", ds / epsilon (ds), "*eps"
-       end if
-    end if
-  end subroutine compare_sigma
-
-  subroutine compare_squared (n1, s1, n2, s2, s, threshold, tolerance)
-    real(kind=omega_prec), intent(in) :: s1, s2
-    character(len=*), intent(in) :: n1, n2
-    integer, dimension(:), intent(in) :: s
-    real(kind=omega_prec), intent(in) :: threshold
-    integer, intent(in), optional :: tolerance
-    real(kind=omega_prec) :: ds, tolerance_local
-    character(len=2) :: num
-    if (present (tolerance)) then
-       tolerance_local = tolerance * epsilon (tolerance_local)
-    else
-       tolerance_local = 1000 * epsilon (tolerance_local)
-    end if
-    if (max (s1, s2) >= threshold) then
-       ds = abs (s1 - s2) / (s1 + s2)
-       if (ds > tolerance_local) then
-          write (unit = num, fmt = "(I2)") size (s)
-          write (unit = *, fmt = "(" // num // &
-               "I3, 2(2X,A1,': ',E17.10), 2X,A2,1X,E8.2, 1X,A1,1X,E8.2, A4)") &
-               s, n1, s1, n2, s2, "d=", ds, "=", ds / epsilon (ds), "*eps"
-       end if
-    end if
-  end subroutine compare_squared
-
-  subroutine compare (n1, a1, n2, a2, s, threshold, tolerance)
-    complex(kind=omega_prec), intent(in) :: a1, a2
-    character(len=*), intent(in) :: n1, n2
-    integer, dimension(:), intent(in) :: s
-    real(kind=omega_prec), intent(in) :: threshold
-    integer, intent(in), optional :: tolerance
-    real(kind=omega_prec) :: ds, tolerance_local
-    character(len=2) :: num
-    if (present (tolerance)) then
-       tolerance_local = tolerance * epsilon (tolerance_local)
-    else
-       tolerance_local = 1000 * epsilon (tolerance_local)
-    end if
-    if (max (abs (a1), abs (a2)) >= threshold) then
-       ds = abs (a1 - a2) / (abs (a1) + abs (a2))
-       if (ds > tolerance_local) then
-          write (unit = num, fmt = "(I2)") size (s)
-          write (unit = *, fmt = "(" // num // &
-               "I3, 2(2X,A1,': (',E10.3,',',E10.3,')'), 2X,A2,1X,E8.2, 1X,A1,1X,E8.2, A3)") &
-               s, n1, a1, n2, a2, "d=", ds, "=", ds / epsilon (ds), "eps"
-       end if
-    end if
-  end subroutine compare
-
-  subroutine compare_omega (n, omega1, omega2, roots, masses, states, tolerance)
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    interface
-       pure function omega1 (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega1
-       pure function omega2 (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega2
-    end interface
-    integer :: k, j
-    real(kind=omega_prec) :: threshold, threshold1, threshold2
-    complex(kind=omega_prec) :: a1, a2
-    real(kind=omega_prec), dimension(0:3,size(masses)) :: p
-    integer, dimension(size(masses)) :: s, nstates
-    if (present (states)) then
-       nstates = states
-    else
-       nstates = 2
-    end if
-    call beams (roots, masses(1), masses(2), p(:,1), p(:,2))
-    do k = 1, n
-       if (any (masses(3:) > 0)) then
-          call massive_decay (roots, masses(3:), p(:,3:))
-       else
-          call massless_isotropic_decay (roots, p(:,3:))
-       end if
-       threshold1 = omega_sum (omega1, p, states) &
-            / num_states (size(s) - 2, nstates(3:)) / 1000
-       threshold2 = omega_sum (omega2, p, states) &
-            / num_states (size(s) - 2, nstates(3:)) / 1000
-       threshold = max (threshold1, threshold2)
-       s = -1
-       loop_spins: do
-          a1 = omega1 (p, s)
-          a2 = omega2 (p, s)
-          call compare ("1", a1, "2", a2, s, threshold, tolerance)
-          do j = size (masses), 1, -1
-             select case (nstates (j))
-             case (3)
-                s(j) = modulo (s(j) + 2, 3) - 1
-             case (2)
-                s(j) = - s(j)
-             case (1)
-                s(j) = -1
-             case default
-                s(j) = -1
-             end select
-             if (s(j) /= -1) then
-                cycle loop_spins
-             end if
-          end do
-          exit loop_spins
-       end do loop_spins
-    end do
-  end subroutine compare_omega
-
-  subroutine compare_omega_sum (n, omega1, omega2, roots, masses, states, tolerance)
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    interface
-       pure function omega1 (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega1
-       pure function omega2 (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega2
-    end interface
-    integer :: k
-    real(kind=omega_prec) :: s1, s2
-    real(kind=omega_prec), dimension(0:3,size(masses)) :: p
-    integer, dimension(:), allocatable :: zero
-    allocate (zero(num_states(size(masses),states)))
-    call beams (roots, masses(1), masses(2), p(:,1), p(:,2))
-    zero = 0
-    do k = 1, n
-       if (any (masses(3:) > 0)) then
-          call massive_decay (roots, masses(3:), p(:,3:))
-       else
-          call massless_isotropic_decay (roots, p(:,3:))
-       end if
-       call omega_sum_nonzero (s1, omega1, p, zero, k, states)
-       call omega_sum_nonzero (s2, omega2, p, zero, k, states)
-       call compare_sigma ("1", s1, "2", s2, THRESHOLD_SIGMA, tolerance)
-    end do
-    deallocate (zero)
-  end subroutine compare_omega_sum
-
-  subroutine symmetry_omega (n, omega, roots, masses, sign, n1, n2, states, tolerance)
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, intent(in) :: sign, n1, n2
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    interface
-       pure function omega (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega
-    end interface
-    integer, dimension(size(masses)) :: s, sx, nstates
-    real(kind=omega_prec), dimension(0:3,size(masses)) :: p, px
-    integer :: k, j
-    complex(kind=omega_prec) :: a, ax
-    real(kind=omega_prec) :: threshold
-    if (present (states)) then
-       nstates = states
-    else
-       nstates = 2
-    end if
-    call beams (roots, masses(1), masses(2), p(:,1), p(:,2))
-    do k = 1, n
-       if (any (masses(3:) > 0)) then
-          call massive_decay (roots, masses(3:), p(:,3:))
-       else
-          call massless_isotropic_decay (roots, p(:,3:))
-       end if
-       threshold = omega_sum (omega, p, states) &
-            / num_states (size(s) - 2, nstates(3:)) / 1000
-       s = -1
-       loop_spins: do
-          px = p
-          px(:,n1) = p(:,n2)
-          px(:,n2) = p(:,n1)
-          sx = s
-          sx(n1) = s(n2)
-          sx(n2) = s(n1)
-          a = omega (p, s)
-          ax = sign * omega (px, sx)
-          call compare ("r", a, "i", ax, s, threshold, tolerance)
-          do j = size (masses), 1, -1
-             select case (nstates (j))
-             case (3)
-                s(j) = modulo (s(j) + 2, 3) - 1
-             case (2)
-                s(j) = - s(j)
-             case (1)
-                s(j) = -1
-             case default
-                s(j) = -1
-             end select
-             if (s(j) /= -1) then
-                cycle loop_spins
-             end if
-          end do
-          exit loop_spins
-       end do loop_spins
-    end do
-  end subroutine symmetry_omega
-  
-  subroutine check_omega (tag, n, omega1, omega2, &
-       roots, masses, symmetry, states, tolerance, mode)
-    character(len=*), intent(in) :: tag
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, dimension(0:,:), intent(in), optional :: symmetry
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    character(len=*), intent(in), optional :: mode
-    interface
-       pure function omega1 (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega1
-       pure function omega2 (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega2
-    end interface
-    integer :: i
-    character(len=8) :: mode_local
-    character(len=130) :: tags
-    print *, trim (tag) // ":"
-    call compare_omega_sum (n, omega1, omega2, roots, masses, states, tolerance)
-    return
-    print *, trim (tag) // " (polarized):"
-    call compare_omega (n, omega1, omega2, roots, masses, states, tolerance)
-    if (present (symmetry)) then
-       do i = 1, size (symmetry, dim=2)
-          write (unit = tags, fmt = "('(',I1,'<>',I1,')')") symmetry(1,i), symmetry(2,i)
-          if (symmetry(0,i) > 0) then
-             print *, trim (tag) // " - " // trim (tags) // ":"
-          else
-             print *, trim (tag) // " + " // trim (tags) // ":"
-          end if
-          call symmetry_omega (n, omega1, roots, masses, symmetry(0,i), &
-               symmetry(1,i), symmetry(2,i), states, tolerance)
-       end do
-    end if
-  end subroutine check_omega
-
-  subroutine ward_omega (n, omega, roots, masses, i, states, tolerance)
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, intent(in) :: i
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    interface
-       pure function omega (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega
-    end interface
-    integer, dimension(size(masses)) :: s, nstates
-    real(kind=omega_prec), dimension(0:3,size(masses)) :: p
-    integer :: k, j
-    complex(kind=omega_prec) :: a
-    real(kind=omega_prec) :: a2, threshold
-    if (present (states)) then
-       nstates = states
-    else
-       nstates = 2
-    end if
-    call beams (roots, masses(1), masses(2), p(:,1), p(:,2))
-    do k = 1, n
-       if (any (masses(3:) > 0)) then
-          call massive_decay (roots, masses(3:), p(:,3:))
-       else
-          call massless_isotropic_decay (roots, p(:,3:))
-       end if
-       threshold = omega_sum (omega, p, states) &
-            / num_states (size(s) - 2, nstates(3:)) / 1000
-       s = -1
-       loop_spins: do
-          s(i) = -1
-          a = omega (p, s)
-          a2 = a * conjg (a)
-          print *, s, a2
-          if (nstates(i) == 3) then
-             s(i) = 0
-             a = omega (p, s)
-             a2 = a * conjg (a)
-             print *, s, a2
-          end if
-          s(i) = 1
-          a = omega (p, s)
-          a2 = a * conjg (a)
-          print *, s, a2
-          s(i) = 4
-          a = omega (p, s)
-          a2 = a * conjg (a)
-          print *, s, a2
-          ! call compare ("r", a, "i", ax, s, threshold, tolerance)
-          do j = size (masses), 1, -1
-             if (j /= i) then
-                select case (nstates (j))
-                case (3)
-                   s(j) = modulo (s(j) + 2, 3) - 1
-                case (2)
-                   s(j) = - s(j)
-                case (1)
-                   s(j) = -1
-                case default
-                   s(j) = -1
-                end select
-                if (s(j) /= -1) then
-                   cycle loop_spins
-                end if
-             end if
-          end do
-          exit loop_spins
-       end do loop_spins
-    end do
-  end subroutine ward_omega
-  
-  subroutine ward4 (n, omega, madgraph, roots, masses, i, states, tolerance)
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, intent(in) :: i
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    interface
-       pure function omega (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega
-       function madgraph (p1, p2, p3, p4, hel) result (m)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:3), intent(in) :: p1, p2, p3, p4
-         integer, dimension(4), intent(in) :: hel
-       end function madgraph
-    end interface
-    integer, dimension(size(masses)) :: s, nstates
-    real(kind=omega_prec), dimension(0:3,size(masses)) :: p
-    integer :: k, j
-    complex(kind=omega_prec) :: a
-    real(kind=omega_prec) :: a2, m2, threshold
-    if (present (states)) then
-       nstates = states
-    else
-       nstates = 2
-    end if
-    call beams (roots, masses(1), masses(2), p(:,1), p(:,2))
-    do k = 1, n
-       if (any (masses(3:) > 0)) then
-          call massive_decay (roots, masses(3:), p(:,3:))
-       else
-          call massless_isotropic_decay (roots, p(:,3:))
-       end if
-       threshold = omega_sum (omega, p, states) &
-            / num_states (size(s) - 2, nstates(3:)) / 1000
-       s = -1
-       loop_spins: do
-          s(i) = -1
-          a = omega (p, s)
-          a2 = a * conjg (a)
-          m2 = madgraph (p(:,1), p(:,2), p(:,3), p(:,4), s)
-          print *, s, a2, m2
-          if (nstates(i) == 3) then
-             s(i) = 0
-             a = omega (p, s)
-             a2 = a * conjg (a)
-             m2 = madgraph (p(:,1), p(:,2), p(:,3), p(:,4), s)
-             print *, s, a2, m2
-          end if
-          s(i) = 1
-          a = omega (p, s)
-          a2 = a * conjg (a)
-          m2 = madgraph (p(:,1), p(:,2), p(:,3), p(:,4), s)
-          print *, s, a2, m2
-          s(i) = 4
-          a = omega (p, s)
-          a2 = a * conjg (a)
-          m2 = madgraph (p(:,1), p(:,2), p(:,3), p(:,4), s)
-          print *, s, a2, m2
-          ! call compare ("r", a, "i", ax, s, threshold, tolerance)
-          do j = size (masses), 1, -1
-             if (j /= i) then
-                select case (nstates (j))
-                case (3)
-                   s(j) = modulo (s(j) + 2, 3) - 1
-                case (2)
-                   s(j) = - s(j)
-                case (1)
-                   s(j) = -1
-                case default
-                   s(j) = -1
-                end select
-                if (s(j) /= -1) then
-                   cycle loop_spins
-                end if
-             end if
-          end do
-          exit loop_spins
-       end do loop_spins
-    end do
-  end subroutine ward4
-  
-  subroutine compare4_madgraph (n, omega, madgraph, roots, masses, states, tolerance)
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    interface
-       pure function omega (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega
-       function madgraph (p1, p2, p3, p4, hel) result (m)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:3), intent(in) :: p1, p2, p3, p4
-         integer, dimension(4), intent(in) :: hel
-       end function madgraph
-    end interface
-    integer :: k, j
-    real(kind=omega_prec) :: threshold, sm, so
-    complex(kind=omega_prec) :: ao
-    real(kind=omega_prec), dimension(0:3,4) :: p
-    integer, dimension(4) :: s, nstates
-    if (present (states)) then
-       nstates = states
-    else
-       nstates = 2
-    end if
-    call beams (roots, masses(1), masses(2), p(:,1), p(:,2))
-    do k = 1, n
-       if (any (masses(3:) > 0)) then
-          call massive_decay (roots, masses(3:), p(:,3:))
-       else
-          call massless_isotropic_decay (roots, p(:,3:))
-       end if
-       s = -1
-       threshold = omega_sum (omega, p, states) &
-            / num_states (size(s) - 2, nstates(3:)) / 1000
-       loop_spins: do
-          sm = madgraph (p(:,1), p(:,2), p(:,3), p(:,4), s)
-          ao = omega (p, s)
-          so = ao * conjg (ao)
-          call compare_squared ("o", so, "m", sm, s, threshold, tolerance)
-          do j = size (masses), 1, -1
-             select case (nstates (j))
-             case (3)
-                s(j) = modulo (s(j) + 2, 3) - 1
-             case (2)
-                s(j) = - s(j)
-             case (1)
-                s(j) = -1
-             case default
-                s(j) = -1
-             end select
-             if (s(j) /= -1) then
-                cycle loop_spins
-             end if
-          end do
-          exit loop_spins
-       end do loop_spins
-    end do
-  end subroutine compare4_madgraph
-
-  subroutine compare5_madgraph (n, omega, madgraph, roots, masses, states, tolerance)
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    interface
-       pure function omega (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega
-       function madgraph (p1, p2, p3, p4, p5, hel) result (m)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:3), intent(in) :: p1, p2, p3, p4, p5
-         integer, dimension(5), intent(in) :: hel
-       end function madgraph
-    end interface
-    integer :: k, j
-    real(kind=omega_prec) :: threshold, sm, so
-    complex(kind=omega_prec) :: ao
-    real(kind=omega_prec), dimension(0:3,5) :: p
-    integer, dimension(5) :: s, nstates
-    if (present (states)) then
-       nstates = states
-    else
-       nstates = 2
-    end if
-    call beams (roots, masses(1), masses(2), p(:,1), p(:,2))
-    do k = 1, n
-       if (any (masses(3:) > 0)) then
-          call massive_decay (roots, masses(3:), p(:,3:))
-       else
-          call massless_isotropic_decay (roots, p(:,3:))
-       end if
-       threshold = omega_sum (omega, p, states) &
-            / num_states (size(s) - 2, nstates(3:)) / 1000
-       s = -1
-       loop_spins: do
-          sm = madgraph (p(:,1), p(:,2), p(:,3), p(:,4), p(:,5), s)
-          ao = omega (p, s)
-          so = ao * conjg (ao)
-          call compare_squared ("o", so, "m", sm, s, threshold, tolerance)
-          do j = size (masses), 1, -1
-             select case (nstates (j))
-             case (3)
-                s(j) = modulo (s(j) + 2, 3) - 1
-             case (2)
-                s(j) = - s(j)
-             case (1)
-                s(j) = -1
-             case default
-                s(j) = -1
-             end select
-             if (s(j) /= -1) then
-                cycle loop_spins
-             end if
-          end do
-          exit loop_spins
-       end do loop_spins
-    end do
-  end subroutine compare5_madgraph
-
-  subroutine compare6_madgraph (n, omega, madgraph, roots, masses, states, tolerance)
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    interface
-       pure function omega (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega
-       function madgraph (p1, p2, p3, p4, p5, p6, hel) result (m)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:3), intent(in) :: &
-              p1, p2, p3, p4, p5, p6
-         integer, dimension(6), intent(in) :: hel
-       end function madgraph
-    end interface
-    integer :: k, j
-    real(kind=omega_prec) :: threshold, sm, so
-    complex(kind=omega_prec) :: ao
-    real(kind=omega_prec), dimension(0:3,6) :: p
-    integer, dimension(6) :: s, nstates
-    if (present (states)) then
-       nstates = states
-    else
-       nstates = 2
-    end if
-    call beams (roots, masses(1), masses(2), p(:,1), p(:,2))
-    do k = 1, n
-       if (any (masses(3:) > 0)) then
-          call massive_decay (roots, masses(3:), p(:,3:))
-       else
-          call massless_isotropic_decay (roots, p(:,3:))
-       end if
-       threshold = omega_sum (omega, p, states) &
-            / num_states (size(s) - 2, nstates(3:)) / 1000
-       s = -1
-       loop_spins: do
-          sm = madgraph (p(:,1), p(:,2), p(:,3), p(:,4), p(:,5), p(:,6), s)
-          ao = omega (p, s)
-          so = ao * conjg (ao)
-          call compare_squared ("o", so, "m", sm, s, threshold, tolerance)
-          do j = size (masses), 1, -1
-             select case (nstates (j))
-             case (3)
-                s(j) = modulo (s(j) + 2, 3) - 1
-             case (2)
-                s(j) = - s(j)
-             case (1)
-                s(j) = -1
-             case default
-                s(j) = -1
-             end select
-             if (s(j) /= -1) then
-                cycle loop_spins
-             end if
-          end do
-          exit loop_spins
-       end do loop_spins
-    end do
-  end subroutine compare6_madgraph
-
-  subroutine compare7_madgraph (n, omega, madgraph, roots, masses, states, tolerance)
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    interface
-       pure function omega (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega
-       function madgraph (p1, p2, p3, p4, p5, p6, p7, hel) result (m)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:3), intent(in) :: &
-              p1, p2, p3, p4, p5, p6, p7
-         integer, dimension(7), intent(in) :: hel
-       end function madgraph
-    end interface
-    integer :: k, j
-    real(kind=omega_prec) :: threshold, sm, so
-    complex(kind=omega_prec) :: ao
-    real(kind=omega_prec), dimension(0:3,7) :: p
-    integer, dimension(7) :: s, nstates
-    if (present (states)) then
-       nstates = states
-    else
-       nstates = 2
-    end if
-    call beams (roots, masses(1), masses(2), p(:,1), p(:,2))
-    do k = 1, n
-       if (any (masses(3:) > 0)) then
-          call massive_decay (roots, masses(3:), p(:,3:))
-       else
-          call massless_isotropic_decay (roots, p(:,3:))
-       end if
-       threshold = omega_sum (omega, p, states) &
-            / num_states (size(s) - 2, nstates(3:)) / 1000
-       s = -1
-       loop_spins: do
-          sm = madgraph (p(:,1), p(:,2), p(:,3), p(:,4), p(:,5), p(:,6), p(:,7), s)
-          ao = omega (p, s)
-          so = ao * conjg (ao)
-          call compare_squared ("o", so, "m", sm, s, threshold, tolerance)
-          do j = size (masses), 1, -1
-             select case (nstates (j))
-             case (3)
-                s(j) = modulo (s(j) + 2, 3) - 1
-             case (2)
-                s(j) = - s(j)
-             case (1)
-                s(j) = -1
-             case default
-                s(j) = -1
-             end select
-             if (s(j) /= -1) then
-                cycle loop_spins
-             end if
-          end do
-          exit loop_spins
-       end do loop_spins
-    end do
-  end subroutine compare7_madgraph
-
-  subroutine compare8_madgraph (n, omega, madgraph, roots, masses, states, tolerance)
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    interface
-       pure function omega (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega
-       function madgraph (p1, p2, p3, p4, p5, p6, p7, p8, hel) result (m)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:3), intent(in) :: &
-              p1, p2, p3, p4, p5, p6, p7, p8
-         integer, dimension(8), intent(in) :: hel
-       end function madgraph
-    end interface
-    integer :: k, j
-    real(kind=omega_prec) :: threshold, sm, so
-    complex(kind=omega_prec) :: ao
-    real(kind=omega_prec), dimension(0:3,8) :: p
-    integer, dimension(8) :: s, nstates
-    if (present (states)) then
-       nstates = states
-    else
-       nstates = 2
-    end if
-    call beams (roots, masses(1), masses(2), p(:,1), p(:,2))
-    do k = 1, n
-       if (any (masses(3:) > 0)) then
-          call massive_decay (roots, masses(3:), p(:,3:))
-       else
-          call massless_isotropic_decay (roots, p(:,3:))
-       end if
-       threshold = omega_sum (omega, p, states) &
-            / num_states (size(s) - 2, nstates(3:)) / 1000
-       s = -1
-       loop_spins: do
-          sm = madgraph (p(:,1), p(:,2), p(:,3), p(:,4), p(:,5), p(:,6), p(:,7), p(:,8), s)
-          ao = omega (p, s)
-          so = ao * conjg (ao)
-          call compare_squared ("o", so, "m", sm, s, threshold, tolerance)
-          do j = size (masses), 1, -1
-             select case (nstates (j))
-             case (3)
-                s(j) = modulo (s(j) + 2, 3) - 1
-             case (2)
-                s(j) = - s(j)
-             case (1)
-                s(j) = -1
-             case default
-                s(j) = -1
-             end select
-             if (s(j) /= -1) then
-                cycle loop_spins
-             end if
-          end do
-          exit loop_spins
-       end do loop_spins
-    end do
-  end subroutine compare8_madgraph
-
-  subroutine compare_sum4_madgraph (n, omega, madgraph, roots, masses, states, tolerance, mode)
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    character(len=*), intent(in), optional :: mode
-    interface
-       pure function omega (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega
-       function madgraph (p1, p2, p3, p4) result (s)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: s
-         real(kind=omega_prec), dimension(0:3), intent(in) :: &
-              p1, p2, p3, p4
-       end function madgraph
-    end interface
-    real(kind=omega_prec), save :: s = 0
-    integer :: k
-    character(len=8) :: mode_local
-    real(kind=omega_prec) :: sm, so
-    real(kind=omega_prec), dimension(0:3,4) :: p
-    integer, dimension(:), allocatable :: zero
-    if (present (mode)) then
-       mode_local = mode
-    else
-       mode_local = ""
-    end if
-    call beams (roots, masses(1), masses(2), p(:,1), p(:,2))
-    if (trim(mode_local) == "omega") then
-       allocate (zero(num_states(4,states)))
-       zero = 0
-       do k = 1, n
-          if (any (masses(3:) > 0)) then
-             call massive_decay (roots, masses(3:), p(:,3:))
-          else
-             call massless_isotropic_decay (roots, p(:,3:))
-          end if
-          call omega_sum_nonzero (so, omega, p, zero, k, states)
-          s = s + so
-       end do
-       deallocate (zero)
-    else if (trim(mode_local) == "madgraph") then
-       do k = 1, n
-          if (any (masses(3:) > 0)) then
-             call massive_decay (roots, masses(3:), p(:,3:))
-          else
-             call massless_isotropic_decay (roots, p(:,3:))
-          end if
-          s = s + madgraph (p(:,1), p(:,2), p(:,3), p(:,4))
-       end do
-    else
-       allocate (zero(num_states(4,states)))
-       zero = 0
-       do k = 1, n
-          if (any (masses(3:) > 0)) then
-             call massive_decay (roots, masses(3:), p(:,3:))
-          else
-             call massless_isotropic_decay (roots, p(:,3:))
-          end if
-          call omega_sum_nonzero (so, omega, p, zero, k, states)
-          sm = madgraph (p(:,1), p(:,2), p(:,3), p(:,4))
-          call compare_sigma ("o", so, "m", sm, THRESHOLD_SIGMA, tolerance)
-       end do
-       deallocate (zero)
-    end if
-  end subroutine compare_sum4_madgraph
-
-  subroutine compare_sum5_madgraph (n, omega, madgraph, roots, masses, states, tolerance, mode)
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    character(len=*), intent(in), optional :: mode
-    interface
-       pure function omega (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega
-       function madgraph (p1, p2, p3, p4, p5) result (s)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: s
-         real(kind=omega_prec), dimension(0:3), intent(in) :: &
-              p1, p2, p3, p4, p5
-       end function madgraph
-    end interface
-    real(kind=omega_prec), save :: s = 0
-    integer :: k
-    character(len=8) :: mode_local
-    real(kind=omega_prec) :: sm, so
-    real(kind=omega_prec), dimension(0:3,5) :: p
-    integer, dimension(:), allocatable :: zero
-    if (present (mode)) then
-       mode_local = mode
-    else
-       mode_local = ""
-    end if
-    call beams (roots, masses(1), masses(2), p(:,1), p(:,2))
-    if (trim(mode_local) == "omega") then
-       allocate (zero(num_states(5,states)))
-       zero = 0
-       do k = 1, n
-          if (any (masses(3:) > 0)) then
-             call massive_decay (roots, masses(3:), p(:,3:))
-          else
-             call massless_isotropic_decay (roots, p(:,3:))
-          end if
-          call omega_sum_nonzero (so, omega, p, zero, k, states)
-          s = s + so
-       end do
-       deallocate (zero)
-    else if (trim(mode_local) == "madgraph") then
-       do k = 1, n
-          if (any (masses(3:) > 0)) then
-             call massive_decay (roots, masses(3:), p(:,3:))
-          else
-             call massless_isotropic_decay (roots, p(:,3:))
-          end if
-          s = s + madgraph (p(:,1), p(:,2), p(:,3), p(:,4), p(:,5))
-       end do
-    else
-       allocate (zero(num_states(5,states)))
-       zero = 0
-       do k = 1, n
-          if (any (masses(3:) > 0)) then
-             call massive_decay (roots, masses(3:), p(:,3:))
-          else
-             call massless_isotropic_decay (roots, p(:,3:))
-          end if
-          call omega_sum_nonzero (so, omega, p, zero, k, states)
-          sm = madgraph (p(:,1), p(:,2), p(:,3), p(:,4), p(:,5))
-          call compare_sigma ("o", so, "m", sm, THRESHOLD_SIGMA, tolerance)
-       end do
-       deallocate (zero)
-    end if
-  end subroutine compare_sum5_madgraph
-
-  subroutine compare_sum6_madgraph (n, omega, madgraph, roots, masses, states, tolerance, mode)
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    character(len=*), intent(in), optional :: mode
-    interface
-       pure function omega (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega
-       function madgraph (p1, p2, p3, p4, p5, p6) result (s)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: s
-         real(kind=omega_prec), dimension(0:3), intent(in) :: &
-              p1, p2, p3, p4, p5, p6
-       end function madgraph
-    end interface
-    real(kind=omega_prec), save :: s = 0
-    integer :: k
-    character(len=8) :: mode_local
-    real(kind=omega_prec) :: sm, so
-    real(kind=omega_prec), dimension(0:3,6) :: p
-    integer, dimension(:), allocatable :: zero
-    if (present (mode)) then
-       mode_local = mode
-    else
-       mode_local = ""
-    end if
-    call beams (roots, masses(1), masses(2), p(:,1), p(:,2))
-    if (trim(mode_local) == "omega") then
-       allocate (zero(num_states(6,states)))
-       zero = 0
-       do k = 1, n
-          if (any (masses(3:) > 0)) then
-             call massive_decay (roots, masses(3:), p(:,3:))
-          else
-             call massless_isotropic_decay (roots, p(:,3:))
-          end if
-          call omega_sum_nonzero (so, omega, p, zero, k, states)
-          s = s + so
-       end do
-       deallocate (zero)
-    else if (trim(mode_local) == "madgraph") then
-       do k = 1, n
-          if (any (masses(3:) > 0)) then
-             call massive_decay (roots, masses(3:), p(:,3:))
-          else
-             call massless_isotropic_decay (roots, p(:,3:))
-          end if
-          s = s + madgraph (p(:,1), p(:,2), p(:,3), p(:,4), p(:,5), p(:,6))
-       end do
-    else
-       allocate (zero(num_states(6,states)))
-       zero = 0
-       do k = 1, n
-          if (any (masses(3:) > 0)) then
-             call massive_decay (roots, masses(3:), p(:,3:))
-          else
-             call massless_isotropic_decay (roots, p(:,3:))
-          end if
-          call omega_sum_nonzero (so, omega, p, zero, k, states)
-          sm = madgraph (p(:,1), p(:,2), p(:,3), p(:,4), p(:,5), p(:,6))
-          call compare_sigma ("o", so, "m", sm, THRESHOLD_SIGMA, tolerance)
-       end do
-       deallocate (zero)
-    end if
-  end subroutine compare_sum6_madgraph
-
-  subroutine compare_sum7_madgraph (n, omega, madgraph, roots, masses, states, tolerance, mode)
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    character(len=*), intent(in), optional :: mode
-    interface
-       pure function omega (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega
-       function madgraph (p1, p2, p3, p4, p5, p6, p7) result (s)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: s
-         real(kind=omega_prec), dimension(0:3), intent(in) :: &
-              p1, p2, p3, p4, p5, p6, p7
-       end function madgraph
-    end interface
-    real(kind=omega_prec), save :: s = 0
-    integer :: k
-    character(len=8) :: mode_local
-    real(kind=omega_prec) :: sm, so
-    real(kind=omega_prec), dimension(0:3,7) :: p
-    integer, dimension(:), allocatable :: zero
-    if (present (mode)) then
-       mode_local = mode
-    else
-       mode_local = ""
-    end if
-    call beams (roots, masses(1), masses(2), p(:,1), p(:,2))
-    if (trim(mode_local) == "omega") then
-       allocate (zero(num_states(7,states)))
-       zero = 0
-       do k = 1, n
-          if (any (masses(3:) > 0)) then
-             call massive_decay (roots, masses(3:), p(:,3:))
-          else
-             call massless_isotropic_decay (roots, p(:,3:))
-          end if
-          call omega_sum_nonzero (so, omega, p, zero, k, states)
-          s = s + so
-       end do
-       deallocate (zero)
-    else if (trim(mode_local) == "madgraph") then
-       do k = 1, n
-          if (any (masses(3:) > 0)) then
-             call massive_decay (roots, masses(3:), p(:,3:))
-          else
-             call massless_isotropic_decay (roots, p(:,3:))
-          end if
-          s = s + madgraph (p(:,1), p(:,2), p(:,3), p(:,4), p(:,5), p(:,6), p(:,7))
-       end do
-    else
-       allocate (zero(num_states(7,states)))
-       zero = 0
-       do k = 1, n
-          if (any (masses(3:) > 0)) then
-             call massive_decay (roots, masses(3:), p(:,3:))
-          else
-             call massless_isotropic_decay (roots, p(:,3:))
-          end if
-          call omega_sum_nonzero (so, omega, p, zero, k, states)
-          sm = madgraph (p(:,1), p(:,2), p(:,3), p(:,4), p(:,5), p(:,6), p(:,7))
-          call compare_sigma ("o", so, "m", sm, THRESHOLD_SIGMA, tolerance)
-       end do
-       deallocate (zero)
-    end if
-  end subroutine compare_sum7_madgraph
-
-  subroutine compare_sum8_madgraph (n, omega, madgraph, roots, masses, states, tolerance, mode)
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    character(len=*), intent(in), optional :: mode
-    interface
-       pure function omega (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega
-       function madgraph (p1, p2, p3, p4, p5, p6, p7, p8) result (s)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: s
-         real(kind=omega_prec), dimension(0:3), intent(in) :: &
-              p1, p2, p3, p4, p5, p6, p7, p8
-       end function madgraph
-    end interface
-    real(kind=omega_prec), save :: s = 0
-    integer :: k
-    character(len=8) :: mode_local
-    real(kind=omega_prec) :: sm, so
-    real(kind=omega_prec), dimension(0:3,8) :: p
-    integer, dimension(:), allocatable :: zero
-    if (present (mode)) then
-       mode_local = mode
-    else
-       mode_local = ""
-    end if
-    call beams (roots, masses(1), masses(2), p(:,1), p(:,2))
-    if (trim(mode_local) == "omega") then
-       allocate (zero(num_states(8,states)))
-       zero = 0
-       do k = 1, n
-          if (any (masses(3:) > 0)) then
-             call massive_decay (roots, masses(3:), p(:,3:))
-          else
-             call massless_isotropic_decay (roots, p(:,3:))
-          end if
-          call omega_sum_nonzero (so, omega, p, zero, k, states)
-          s = s + so
-       end do
-       deallocate (zero)
-    else if (trim(mode_local) == "madgraph") then
-       do k = 1, n
-          if (any (masses(3:) > 0)) then
-             call massive_decay (roots, masses(3:), p(:,3:))
-          else
-             call massless_isotropic_decay (roots, p(:,3:))
-          end if
-          s = s + madgraph (p(:,1), p(:,2), p(:,3), p(:,4), p(:,5), p(:,6), p(:,7), p(:,8))
-       end do
-    else
-       allocate (zero(num_states(8,states)))
-       zero = 0
-       do k = 1, n
-          if (any (masses(3:) > 0)) then
-             call massive_decay (roots, masses(3:), p(:,3:))
-          else
-             call massless_isotropic_decay (roots, p(:,3:))
-          end if
-          call omega_sum_nonzero (so, omega, p, zero, k, states)
-          sm = madgraph (p(:,1), p(:,2), p(:,3), p(:,4), p(:,5), p(:,6), p(:,7), p(:,8))
-          call compare_sigma ("o", so, "m", sm, THRESHOLD_SIGMA, tolerance)
-       end do
-       deallocate (zero)
-    end if
-  end subroutine compare_sum8_madgraph
-
-  subroutine check4_madgraph (tag, n, omega, smadgraph, madgraph, &
-       roots, masses, symmetry, states, tolerance, mode)
-    character(len=*), intent(in) :: tag
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, dimension(0:,:), intent(in), optional :: symmetry
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    character(len=*), intent(in), optional :: mode
-    interface
-       pure function omega (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega
-       function smadgraph (p1, p2, p3, p4) result (s)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: s
-         real(kind=omega_prec), dimension(0:3), intent(in) :: &
-              p1, p2, p3, p4
-       end function smadgraph
-       function madgraph (p1, p2, p3, p4, hel) result (m)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:3), intent(in) :: &
-              p1, p2, p3, p4
-         integer, dimension(4), intent(in) :: hel
-       end function madgraph
-    end interface
-    integer :: i
-    character(len=8) :: mode_local
-    character(len=130) :: tags
-    if (present (mode)) then
-       mode_local = mode
-    else
-       mode_local = "compare"
-    end if
-    if (trim (mode_local) == "compare") then
-       print *, trim (tag) // ":"
-       call compare_sum4_madgraph (n, omega, smadgraph, roots, masses, states, tolerance)
-       print *, trim (tag) // " (polarized):"
-       call compare4_madgraph (n, omega, madgraph, roots, masses, states, tolerance)
-       if (present (symmetry)) then
-          do i = 1, size (symmetry, dim=2)
-             write (unit = tags, fmt = "('(',I1,'<>',I1,')')") symmetry(1,i), symmetry(2,i)
-             if (symmetry(0,i) > 0) then
-                print *, trim (tag) // " - " // trim (tags) // ":"
-             else
-                print *, trim (tag) // " + " // trim (tags) // ":"
-             end if
-             call symmetry_omega (n, omega, roots, masses, symmetry(0,i), &
-                  symmetry(1,i), symmetry(2,i), states, tolerance)
-          end do
-       end if
-    else
-       print *, trim (tag) // " (" // trim (mode_local) // "):"
-       call compare_sum4_madgraph (n, omega, smadgraph, roots, masses, states, &
-            tolerance, mode = mode_local)
-    end if
-  end subroutine check4_madgraph
-
-  subroutine check5_madgraph (tag, n, omega, smadgraph, madgraph, &
-       roots, masses, symmetry, states, tolerance, mode)
-    character(len=*), intent(in) :: tag
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, dimension(0:,:), intent(in), optional :: symmetry
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    character(len=*), intent(in), optional :: mode
-    interface
-       pure function omega (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega
-       function smadgraph (p1, p2, p3, p4, p5) result (s)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: s
-         real(kind=omega_prec), dimension(0:3), intent(in) :: &
-              p1, p2, p3, p4, p5
-       end function smadgraph
-       function madgraph (p1, p2, p3, p4, p5, hel) result (m)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:3), intent(in) :: &
-              p1, p2, p3, p4, p5
-         integer, dimension(5), intent(in) :: hel
-       end function madgraph
-    end interface
-    integer :: i
-    character(len=8) :: mode_local
-    character(len=130) :: tags
-    if (present (mode)) then
-       mode_local = mode
-    else
-       mode_local = "compare"
-    end if
-    if (trim (mode_local) == "compare") then
-       print *, trim (tag) // ":"
-       call compare_sum5_madgraph (n, omega, smadgraph, roots, masses, states, tolerance)
-       print *, trim (tag) // " (polarized):"
-       call compare5_madgraph (n, omega, madgraph, roots, masses, states, tolerance)
-       if (present (symmetry)) then
-          do i = 1, size (symmetry, dim=2)
-             write (unit = tags, fmt = "('(',I1,'<>',I1,')')") symmetry(1,i), symmetry(2,i)
-             if (symmetry(0,i) > 0) then
-                print *, trim (tag) // " - " // trim (tags) // ":"
-             else
-                print *, trim (tag) // " + " // trim (tags) // ":"
-             end if
-             call symmetry_omega (n, omega, roots, masses, symmetry(0,i), &
-                  symmetry(1,i), symmetry(2,i), states, tolerance)
-          end do
-       end if
-    else
-       print *, trim (tag) // " (" // trim (mode_local) // "):"
-       call compare_sum5_madgraph (n, omega, smadgraph, roots, masses, states, &
-            tolerance, mode = mode_local)
-    end if
-  end subroutine check5_madgraph
-
-  subroutine check6_madgraph (tag, n, omega, smadgraph, madgraph, &
-       roots, masses, symmetry, states, tolerance, mode)
-    character(len=*), intent(in) :: tag
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, dimension(0:,:), intent(in), optional :: symmetry
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    character(len=*), intent(in), optional :: mode
-    interface
-       pure function omega (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega
-       function smadgraph (p1, p2, p3, p4, p5, p6) result (s)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: s
-         real(kind=omega_prec), dimension(0:3), intent(in) :: &
-              p1, p2, p3, p4, p5, p6
-       end function smadgraph
-       function madgraph (p1, p2, p3, p4, p5, p6, hel) result (m)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:3), intent(in) :: &
-              p1, p2, p3, p4, p5, p6
-         integer, dimension(6), intent(in) :: hel
-       end function madgraph
-    end interface
-    integer :: i
-    character(len=8) :: mode_local
-    character(len=130) :: tags
-    if (present (mode)) then
-       mode_local = mode
-    else
-       mode_local = "compare"
-    end if
-    if (trim (mode_local) == "compare") then
-       print *, trim (tag) // ":"
-       call compare_sum6_madgraph (n, omega, smadgraph, roots, masses, states, tolerance)
-       print *, trim (tag) // " (polarized):"
-       call compare6_madgraph (n, omega, madgraph, roots, masses, states, tolerance)
-       if (present (symmetry)) then
-          do i = 1, size (symmetry, dim=2)
-             write (unit = tags, fmt = "('(',I1,'<>',I1,')')") symmetry(1,i), symmetry(2,i)
-             if (symmetry(0,i) > 0) then
-                print *, trim (tag) // " - " // trim (tags) // ":"
-             else
-                print *, trim (tag) // " + " // trim (tags) // ":"
-             end if
-             call symmetry_omega (n, omega, roots, masses, symmetry(0,i), &
-                  symmetry(1,i), symmetry(2,i), states, tolerance)
-          end do
-       end if
-    else
-       print *, trim (tag) // " (" // trim (mode_local) // "):"
-       call compare_sum6_madgraph (n, omega, smadgraph, roots, masses, states, &
-            tolerance, mode = mode_local)
-    end if
-  end subroutine check6_madgraph
-
-  subroutine check7_madgraph (tag, n, omega, smadgraph, madgraph, &
-       roots, masses, symmetry, states, tolerance, mode)
-    character(len=*), intent(in) :: tag
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, dimension(0:,:), intent(in), optional :: symmetry
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    character(len=*), intent(in), optional :: mode
-    interface
-       pure function omega (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega
-       function smadgraph (p1, p2, p3, p4, p5, p6, p7) result (s)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: s
-         real(kind=omega_prec), dimension(0:3), intent(in) :: &
-              p1, p2, p3, p4, p5, p6, p7
-       end function smadgraph
-       function madgraph (p1, p2, p3, p4, p5, p6, p7, hel) result (m)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:3), intent(in) :: &
-              p1, p2, p3, p4, p5, p6, p7
-         integer, dimension(7), intent(in) :: hel
-       end function madgraph
-    end interface
-    integer :: i
-    character(len=8) :: mode_local
-    character(len=130) :: tags
-    if (present (mode)) then
-       mode_local = mode
-    else
-       mode_local = "compare"
-    end if
-    if (trim (mode_local) == "compare") then
-       print *, trim (tag) // ":"
-       call compare_sum7_madgraph (n, omega, smadgraph, roots, masses, states, tolerance)
-       print *, trim (tag) // " (polarized):"
-       call compare7_madgraph (n, omega, madgraph, roots, masses, states, tolerance)
-       if (present (symmetry)) then
-          do i = 1, size (symmetry, dim=2)
-             write (unit = tags, fmt = "('(',I1,'<>',I1,')')") symmetry(1,i), symmetry(2,i)
-             if (symmetry(0,i) > 0) then
-                print *, trim (tag) // " - " // trim (tags) // ":"
-             else
-                print *, trim (tag) // " + " // trim (tags) // ":"
-             end if
-             call symmetry_omega (n, omega, roots, masses, symmetry(0,i), &
-                  symmetry(1,i), symmetry(2,i), states, tolerance)
-          end do
-       end if
-    else
-       print *, trim (tag) // " (" // trim (mode_local) // "):"
-       call compare_sum7_madgraph (n, omega, smadgraph, roots, masses, states, &
-            tolerance, mode = mode_local)
-    end if
-  end subroutine check7_madgraph
-
-  subroutine check8_madgraph (tag, n, omega, smadgraph, madgraph, &
-       roots, masses, symmetry, states, tolerance, mode)
-    character(len=*), intent(in) :: tag
-    integer, intent(in) :: n
-    real(kind=omega_prec), intent(in) :: roots
-    real(kind=omega_prec), dimension(:), intent(in) :: masses
-    integer, dimension(0:,:), intent(in), optional :: symmetry
-    integer, dimension(:), intent(in), optional :: states
-    integer, intent(in), optional :: tolerance
-    character(len=*), intent(in), optional :: mode
-    interface
-       pure function omega (p, s) result (m)
-         use omega_kinds
-         implicit none
-         complex(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:,:), intent(in) :: p
-         integer, dimension(:), intent(in) :: s
-       end function omega
-       function smadgraph (p1, p2, p3, p4, p5, p6, p7, p8) result (s)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: s
-         real(kind=omega_prec), dimension(0:3), intent(in) :: &
-              p1, p2, p3, p4, p5, p6, p7, p8
-       end function smadgraph
-       function madgraph (p1, p2, p3, p4, p5, p6, p7, p8, hel) result (m)
-         use omega_kinds
-         implicit none
-         real(kind=omega_prec) :: m
-         real(kind=omega_prec), dimension(0:3), intent(in) :: &
-              p1, p2, p3, p4, p5, p6, p7, p8
-         integer, dimension(8), intent(in) :: hel
-       end function madgraph
-    end interface
-    integer :: i
-    character(len=8) :: mode_local
-    character(len=130) :: tags
-    if (present (mode)) then
-       mode_local = mode
-    else
-       mode_local = "compare"
-    end if
-    if (trim (mode_local) == "compare") then
-       print *, trim (tag) // ":"
-       call compare_sum8_madgraph (n, omega, smadgraph, roots, masses, states, tolerance)
-       print *, trim (tag) // " (polarized):"
-       call compare8_madgraph (n, omega, madgraph, roots, masses, states, tolerance)
-       if (present (symmetry)) then
-          do i = 1, size (symmetry, dim=2)
-             write (unit = tags, fmt = "('(',I1,'<>',I1,')')") symmetry(1,i), symmetry(2,i)
-             if (symmetry(0,i) > 0) then
-                print *, trim (tag) // " - " // trim (tags) // ":"
-             else
-                print *, trim (tag) // " + " // trim (tags) // ":"
-             end if
-             call symmetry_omega (n, omega, roots, masses, symmetry(0,i), &
-                  symmetry(1,i), symmetry(2,i), states, tolerance)
-          end do
-       end if
-    else
-       print *, trim (tag) // " (" // trim (mode_local) // "):"
-       call compare_sum8_madgraph (n, omega, smadgraph, roots, masses, states, &
-            tolerance, mode = mode_local)
-    end if
-  end subroutine check8_madgraph
-
-end module testbed_old