!=======================================================================
! Module for angular momentum coupling coefficients
! Contains: 9j symbols, Racah coefficients, 6j symbols, 3j symbols,
!           and Clebsch-Gordan coefficients
!
! Converted from legacy Fortran 3j6j9jmod.f90
!=======================================================================

module j_coupling_symbols

  implicit none
  
  private
  
  ! Public subroutines and functions
  public :: coef9j, racah, sj2i, par, tj2i, tji, tj, clebr, clebi2, cleb
  public :: th_facinit
  
  ! Parameters
  integer, parameter :: lfactc = 2000
  
  ! Shared data for factorial table
  logical, save :: first = .true.
  integer, save :: lfact
  real(kind=8), allocatable, save :: faclog(:)

contains

  !=======================================================================
  ! Function coef9j: Calculate 9j symbols using Racah coefficients
  ! 
  ! The 9j symbol is arranged as:
  !     j1  j2  j3
  !     j4  j5  j6
  !     j7  j8  j9
  !
  ! Arguments use the convention that when numbered 1-9, they correspond
  ! to: j1, j2, j3, j4, j5, j6, j7, j8, j9
  !=======================================================================
  real(kind=8) function coef9j(j1, j2, j3, j4, j5, j6, j7, j8, j9)
    
    implicit none
    real(kind=8), intent(in) :: j1, j2, j3, j4, j5, j6, j7, j8, j9
    real(kind=8) :: lt(9)
    integer :: i, imin, m1, m2, m3, m4, kex, k1, k2, k3, nmin, nmax
    integer :: n, imin_adj, k
    real(kind=8) :: lmin, rac, ph, w1, kp, u9, l1, rt
    
    u9 = 0.0d0
    lt(1) = j1
    lt(2) = j2
    lt(3) = j4
    lt(4) = j5
    lt(5) = j3
    lt(6) = j6
    lt(7) = j7
    lt(8) = j8
    lt(9) = j9
    
    ! Find minimum element and its position
    lmin = lt(1)
    imin = 1
    do i = 2, 9
      if (lt(i) < lmin) then
        lmin = lt(i)
        imin = i
      end if
    end do
    
    kex = 0
    
    ! Perform rotations based on position of minimum element
    select case(imin)
    case(1:4)
      m1 = ((imin - 1) / 2) + 1
      m1 = m1 + m1 - 1
      m2 = m1 + 1
      m3 = ((imin - 1) / 2) + 1 + 4  ! m3 = mm + 4, where mm = ((imin-1)/2) + 1
      
      ! Swap elements
      l1 = lt(7)
      lt(7) = lt(m1)
      lt(m1) = l1
      l1 = lt(8)
      lt(8) = lt(m2)
      lt(m2) = l1
      l1 = lt(9)
      lt(9) = lt(m3)
      lt(m3) = l1
      
      imin = imin + (7 - m1)
      
      ! Continue with m1=5, m2=6
      m1 = 5
      m2 = 6
      m3 = imin - 6
      m4 = m3 + 2
      
    case(5:6)
      kex = 1
      m1 = 7
      m2 = 8
      m3 = imin + imin - 9
      m4 = m3 + 1
      
    case(7:8)
      kex = 1
      m1 = 5
      m2 = 6
      m3 = imin - 6
      m4 = m3 + 2
      
    case(9)
      ! imin == 9: No rearrangement needed, skip swaps
      m1 = 0  ! Sentinel value to skip swaps
      
    end select
    
    ! Perform final swaps (only if imin <= 8)
    ! Note: case(9) sets m1=0, which ensures this block is skipped correctly
    if (imin /= 9) then
      l1 = lt(m1)
      lt(m1) = lt(m3)
      lt(m3) = l1
      l1 = lt(m2)
      lt(m2) = lt(m4)
      lt(m4) = l1
      l1 = lt(9)
      lt(9) = lt(imin)
      lt(imin) = l1
    end if
    
    ! Special case: lt(9) <= 0 (lt(9) < 0 or lt(9) == 0)
    if (lt(9) <= 1.0d-10) then
      ! Check if lt(5) == lt(6) and lt(7) == lt(8)
      ! If not, return zero (go to 1000)
      if (abs(lt(5) - lt(6)) > 1.0d-10 .or. abs(lt(7) - lt(8)) > 1.0d-10) then
        u9 = 0.0d0
      else
        rt = (2.0d0 * lt(5) + 1.0d0) * (2.0d0 * lt(7) + 1.0d0)
        k = nint(lt(5) + lt(7) - lt(1) - lt(4))
        rac = racah(2.0d0*lt(1), 2.0d0*lt(2), 2.0d0*lt(3), 2.0d0*lt(4), 2.0d0*lt(5), 2.0d0*lt(7))
        ph = 1.0d0
        if (mod(k/2, 2) /= 0) ph = -1.0d0
        u9 = (rac / sqrt(rt)) * ph
      end if
    else
      ! General case: lt(9) > 0
      k1 = iabs(nint(lt(2) - lt(7)))
      k2 = iabs(nint(lt(3) - lt(5)))
      k3 = iabs(nint(lt(4) - lt(9)))
      nmin = max(k1, k2, k3)
      
      k1 = nint(lt(2) + lt(7))
      k2 = nint(lt(3) + lt(5))
      k3 = nint(lt(4) + lt(9))
      nmax = min(k1, k2, k3)
      
      if (nmin <= nmax) then
        do n = nmin, nmax, 2
          w1 = dble(2*n + 1)
          rac = racah(2.0d0*lt(2), 2.0d0*lt(5), 2.0d0*lt(7), 2.0d0*lt(3), 2.0d0*lt(1), 2.0d0*dble(n))
          if (abs(rac) > 1.0d-14) then
            w1 = w1 * rac
            rac = racah(2.0d0*lt(2), 2.0d0*lt(4), 2.0d0*lt(7), 2.0d0*lt(9), 2.0d0*lt(8), 2.0d0*dble(n))
            if (abs(rac) > 1.0d-14) then
              w1 = w1 * rac
              rac = racah(2.0d0*lt(3), 2.0d0*lt(4), 2.0d0*lt(5), 2.0d0*lt(9), 2.0d0*lt(6), 2.0d0*dble(n))
              if (abs(rac) > 1.0d-14) then
                u9 = u9 + w1 * rac
              end if
            end if
          end if
        end do
      end if
    end if
    
    ! Apply phase factor if rotation was done
    if (kex /= 0) then
      kp = 0
      do i = 1, 9
        kp = kp + nint(2.0d0 * lt(i))
      end do
      k = kp / 2  ! Sum of j_i
      ph = 1.0d0  ! Phase (-1)^(sum j_i)
      if (mod(k, 2) /= 0) ph = -1.0d0
      u9 = u9 * ph
    end if
    
    coef9j = u9
    
  end function coef9j

  !=======================================================================
  ! Function racah: Calculate Racah coefficients
  ! 
  ! Entry points:
  !   racah(jad, jbd, jcd, jdd, jed, jfd)   - arguments are 2*j (integers)
  !   racahi(ja1, jb1, jc1, jd1, je1, jf1)  - arguments are true integers
  !   racahr(a, b, c, d, e, f)              - arguments are real values
  !=======================================================================
  real(kind=8) function racah(jad, jbd, jcd, jdd, jed, jfd)
    
    implicit none
    real(kind=8), intent(in) :: jad, jbd, jcd, jdd, jed, jfd
    integer :: i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16
    integer :: k, n, il, j1, j2, j3, j4, j5, j6, j7, j5_adj, j6_adj, j7_adj
    integer :: j, kl, k_adj, iad, ibd, jj
    integer :: iv_vals(4)
    real(kind=8) :: ph, h, s, lmin_val
    
    ! Initialize factorial table if needed
    if (first) call th_facinit()
    
    ! Calculate intermediate 'i' values (these are related to sums/differences of j's).
s
    ! and cause many valid Racah coefficients to be zero.
    ! The parity of (j_a+j_b-j_e) etc. is implicitly handled by the triangle inequalities
    ! and the definition of the summation range.
    i1 = nint((jad + jbd - jed + 2.0d0) / 2.0d0)
    i4 = nint((jcd + jdd - jed + 2.0d0) / 2.0d0)
    i7 = nint((jad + jcd - jfd + 2.0d0) / 2.0d0)
    i10 = nint((jbd + jdd - jfd + 2.0d0) / 2.0d0)
    
    i13 = i1 + nint(jed)
    i14 = i4 + nint(jed)
    i15 = i7 + nint(jfd)
    i16 = i10 + nint(jfd)
    i2 = i13 - nint(jad)
    i3 = i13 - nint(jbd)
    i5 = i14 - nint(jcd)
    i6 = i14 - nint(jdd)
    i8 = i15 - nint(jad)
    i9 = i15 - nint(jcd)
    i11 = i16 - nint(jbd)
    i12 = i16 - nint(jdd)
    
    ! Check triangular inequalities
    n = min(i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12) - 1
    
    if (n < 0) then
      racah = 0.0d0
      return
    end if
    
    ! Find minimum value of summation index
    il = max(i13, i14, i15, i16)
    
    !   Check for minimum j value (lmin in original code)
    ! If any argument is zero, use special formula
    if (min(nint(jad), nint(jbd), nint(jcd), nint(jdd), nint(jed), nint(jfd)) == 0) then
      ! Zero argument case: use special formula
      iad = huge(0) ! Initialize with a very large number
      ibd = huge(0) ! Initialize with a very large number
      
      ! Find max and second-max values from i(13), i(14), i(15), i(16)
      do jj = 1, 4
        iv_vals(jj) = 0
      end do
      iv_vals(1) = i13
      iv_vals(2) = i14
      iv_vals(3) = i15
      iv_vals(4) = i16
      
      do jj = 1, 4
        if (iv_vals(jj) < iad) then
          ibd = iad  ! The old smallest becomes the new second smallest
          iad = iv_vals(jj) ! The current value is the new smallest
        else if (iv_vals(jj) < ibd) then
          ibd = iv_vals(jj) ! The current value is smaller than the second smallest
        end if
      end do
      j5 = i13 + i4 - il
      ph = 1.0d0
      if (mod(j5, 2) == 0) ph = -1.0d0
      racah = ph / sqrt(dble(iad * ibd))
      return
    end if
    
    if (il >= lfact) then
      stop 'racah: length factorial table insufficient'
    end if
    
    ! Calculate Racah coefficient
    j1 = il - i13 + 1
    j2 = il - i14 + 1
    j3 = il - i15 + 1
    j4 = il - i16 + 1
    j5 = i13 + i4 - il
    j6 = i15 + i5 - il
    j7 = i16 + i6 - il
    
    ph = 1.0d0
    if (mod(j5, 2) == 0) ph = -1.0d0
    
    h = ph * dexp((faclog(i1) + faclog(i2) + faclog(i3) - faclog(i13 + 1) + &
                   faclog(i4) + faclog(i5) + faclog(i6) - faclog(i14 + 1) + &
                   faclog(i7) + faclog(i8) + faclog(i9) - faclog(i15 + 1) + &
                   faclog(i10) + faclog(i11) + faclog(i12) - faclog(i16 + 1)) * 0.50d0 + &
                  faclog(il + 1) - faclog(j1) - faclog(j2) - faclog(j3) - faclog(j4) - &
                  faclog(j5) - faclog(j6) - faclog(j7))
    
    if (n <= 0) then
      racah = h
      return
    end if
    
    ! Sum over n terms
    s = 1.0d0
    k_adj = n - 1
    kl = il + 1
    j5_adj = j5 - 1
    j6_adj = j6 - 1
    j7_adj = j7 - 1
    
    do j = 1, n
      s = 1.0d0 - ((kl + k_adj) * (j5_adj - k_adj) * (j6_adj - k_adj) * &
                    (j7_adj - k_adj)) * s / ((j1 + k_adj) * (j2 + k_adj) * &
                    (j3 + k_adj) * (j4 + k_adj))
      k_adj = k_adj - 1
    end do
    
    racah = h * s
    
  end function racah

  !=======================================================================
  ! Subroutine th_facinit: Initialize factorial table
  ! Generates logarithm of factorials for numerical stability
  !=======================================================================
  subroutine th_facinit()
    
    implicit none
    integer :: i
    real(kind=8) :: fn
    
    if (.not. allocated(faclog)) then
      allocate(faclog(lfactc))
    end if
    
    first = .false.
    lfact = lfactc
    
    faclog(1) = 0.0d0
    faclog(2) = 0.0d0
    fn = 1.0d0
    
    do i = 3, lfactc
      fn = fn + 1.0d0
      faclog(i) = faclog(i - 1) + dlog(fn)
    end do
    
  end subroutine th_facinit

  !=======================================================================
  ! Function sj2i: Calculate 6j symbol from Racah coefficient
  ! Arguments are 2*j (integers)
  !=======================================================================
  real(kind=8) function sj2i(j1, j2, k2, k1, j3, k3)
    
    implicit none
    real(kind=8), intent(in) :: j1, j2, k2, k1, j3, k3
    
    sj2i = par((j1 + j2 + j3 + k1) / 2.0d0) * racah(j1, j2, j3, k1, k2, k3)
    
  end function sj2i

  !=======================================================================
  ! Function par: Calculate parity factor (-1)^a
  !=======================================================================
  real(kind=8) function par(a)
    
    implicit none
    real(kind=8), intent(in) :: a
    integer :: i
    real(kind=8) :: p, q
    
    i = nint(a)
    i = i / 2
    p = a / 2.0d0
    q = dble(i)
    
    if (abs(p - q) > 1.0d-10) then
      par = -1.0d0
    else
      par = 1.0d0
    end if
    
  end function par

  !=======================================================================
  ! Function tj2i / tji / tj: Calculate 3j symbol from Clebsch-Gordan
  ! 
  ! Entry points:
  !   tj2i(j21, j22, j23, m21, m22, m23) - arguments are 2*j (integers)
  !   tji(j1, j2, j3, m1, m2, m3)        - arguments are true integers
  !   tj(aa1, aa2, aa3, bb1, bb2, bb3)   - arguments are real values
  !=======================================================================
  real(kind=8) function tj2i(j21, j22, j23, m21, m22, m23)
    
    implicit none
    real(kind=8), intent(in) :: j21, j22, j23, m21, m22, m23
    real(kind=8) :: a1, b1, a2, b2, a3, b3
    
    a1 = j21 / 2.0d0
    b1 = m21 / 2.0d0
    a2 = j22 / 2.0d0
    b2 = m22 / 2.0d0
    a3 = j23 / 2.0d0
    b3 = m23 / 2.0d0
    
    tj2i = par(a1 - a2 - b3) * clebr(a1, b1, a2, b2, a3, -b3) / sqrt(2.0d0 * a3 + 1.0d0)
    
  end function tj2i

  !=======================================================================
  ! Function tji: Calculate 3j symbol (integer arguments)
  !=======================================================================
  real(kind=8) function tji(j1, j2, j3, m1, m2, m3)
    
    implicit none
    real(kind=8), intent(in) :: j1, j2, j3, m1, m2, m3
    real(kind=8) :: a1, b1, a2, b2, a3, b3
    
    a1 = j1
    b1 = m1
    a2 = j2
    b2 = m2
    a3 = j3
    b3 = m3
    
    tji = par(a1 - a2 - b3) * clebr(a1, b1, a2, b2, a3, -b3) / sqrt(2.0d0 * a3 + 1.0d0)
    
  end function tji

  !=======================================================================
  ! Function tj: Calculate 3j symbol (real arguments)
  !=======================================================================
  real(kind=8) function tj(aa1, aa2, aa3, bb1, bb2, bb3)
    
    implicit none
    real(kind=8), intent(in) :: aa1, aa2, aa3, bb1, bb2, bb3
    real(kind=8) :: a1, b1, a2, b2, a3, b3
    
    a1 = aa1
    b1 = bb1
    a2 = aa2
    b2 = bb2
    a3 = aa3
    b3 = bb3
    
    tj = par(a1 - a2 - b3) * clebr(a1, b1, a2, b2, a3, -b3) / sqrt(2.0d0 * a3 + 1.0d0)
    
  end function tj

  !=======================================================================
  ! Function clebr / clebi2 / cleb: Calculate Clebsch-Gordan coefficients
  ! 
  ! Entry points:
  !   clebr(a, b, c, d, e, f)       - arguments are real true values
  !   clebi2(ll1, lm1, ll2, lm2, ll3, lm3) - arguments are integers
  !   cleb(i2j1, i2m1, i2j2, i2m2, i2j3, i2m3) - arguments are 2*j
  !=======================================================================
  real(kind=8) function clebr(a, b, c, d, e, f)
    
    implicit none
    real(kind=8), intent(in) :: a, b, c, d, e, f
    integer :: ia, ib, ic, id, ie, ifv
    
    ia = nint(2.0d0 * a)
    ib = nint(2.0d0 * c)
    ic = nint(2.0d0 * e)
    id = nint(2.0d0 * b)
    ie = nint(2.0d0 * d)
    ifv = nint(2.0d0 * f)
    
    clebr = clebsch_gordan_calc(ia, ib, ic, id, ie, ifv)
    
  end function clebr

  !=======================================================================
  ! Function clebi2: Calculate Clebsch-Gordan (integer entry point)
  !=======================================================================
  real(kind=8) function clebi2(ll1, lm1, ll2, lm2, ll3, lm3)
    
    implicit none
    integer, intent(in) :: ll1, lm1, ll2, lm2, ll3, lm3
    integer :: ia, ib, ic, id, ie, ifv
    
    ia = 2 * ll1
    ib = 2 * ll2
    ic = 2 * ll3
    id = 2 * lm1
    ie = 2 * lm2
    ifv = 2 * lm3
    
    clebi2 = clebsch_gordan_calc(ia, ib, ic, id, ie, ifv)
    
  end function clebi2

  !=======================================================================
  ! Function cleb: Calculate Clebsch-Gordan (2*j integer entry point)
  !=======================================================================
  real(kind=8) function cleb(i2j1, i2m1, i2j2, i2m2, i2j3, i2m3)
    
    implicit none
    integer, intent(in) :: i2j1, i2m1, i2j2, i2m2, i2j3, i2m3
    
    cleb = clebsch_gordan_calc(i2j1, i2j2, i2j3, i2m1, i2m2, i2m3)
    
  end function cleb

  !=======================================================================
  ! Internal function: clebsch_gordan_calc
  ! Performs the actual Clebsch-Gordan coefficient calculation
  !=======================================================================
  real(kind=8) function clebsch_gordan_calc(ia_in, ib_in, ic_in, id_in, ie_in, ifv_in) result(rac)
    
    implicit none
    integer, intent(in) :: ia_in, ib_in, ic_in, id_in, ie_in, ifv_in
    integer :: ia, ib, ic, id, ie, ifv
    integer :: k1, k2, k3, k4, k5, k6, k7, nz, nzm1, nzt1, nzt2, nzt3
    integer :: nzt4, nzt5, iabcp, iabc, icab, ibca, iapd, iamd, ibpe, ibme, icpf, icmf
    integer :: nzmic2, nzmic3, nzmi, nzmx, iam, ibm, icm, idm, iem, ifm
    integer :: it
    real(kind=8) :: fibm, ficm, fc2
    real(kind=8) :: signfc, sqfclg, s, s1, ssterm, termlg, h, ph, p, q
    
    rac = 0.0d0
    
    ! Copy input arguments to local variables for modification
    ia = ia_in
    ib = ib_in
    ic = ic_in
    id = id_in
    ie = ie_in
    ifv = ifv_in
    
    ! Initialize factorial table if needed
    if (first) call th_facinit()
    
    ! Check M conservation
    if (id + ie - ifv /= 0) return
    
    ! Check parity of sum
    k1 = ia + ib + ic
    if (mod(k1, 2) /= 0) return
      
    k2 = ic + ia - ib
    k3 = ib + ic - ia
    k4 = ia - iabs(ib - ic)
    k5 = ib - iabs(ic - ia)
    k6 = ic - iabs(ia - ib)
    k7 = min(k1, k2, k3, k4, k5, k6)
    
    if (k7 < 0) return
    
    ! Check parity conditions
    if (mod(ia + id, 2) /= 0) return
    if (mod(ib + ie, 2) /= 0) return
    if (mod(ic + ifv, 2) /= 0) return
    
    ! Check M ranges
    if (ia < abs(id)) return
    if (ib < abs(ie)) return
    if (ic < abs(ifv)) return
    
    signfc = 1.0d0
    iam = ia
    ibm = ib
    icm = ic
    idm = id
    iem = ie
    ifm = ifv
    
    ! Order ia >= ib >= ic
    if (ia < ib) then
      if (ia < ic) then
        ! Swap ia and ib
        it = ia
        ia = ib
        ib = it
        it = id
        id = ie
        ie = it
        signfc = real((-1) ** ((ia + ib - ic) / 2), 8)
      else
        ! Swap ic and ib
        it = ic
        ic = ib
        ib = it
        it = ifv
        ifv = -ie
        ie = -it
        fibm = real(ibm + 1, 8)
        ficm = real(icm + 1, 8)
        signfc = real((-1) ** ((iam - idm) / 2), 8) * sqrt(ficm / fibm)
      end if
    else if (ic < ib) then
      ! Swap ic and ib
      it = ic
      ic = ib
      ib = it
      it = ifv
      ifv = -ie
      ie = -it
      fibm = real(ibm + 1, 8)
      ficm = real(icm + 1, 8)
      signfc = real((-1) ** ((iam - idm) / 2), 8) * sqrt(ficm / fibm)
    end if
    
    ! Special case: jb = 0
    if (ib == 0) then
      rac = signfc
      return
    end if
    
    ! Check if special case for me /= 0
    if (ie /= 0) then
      signfc = signfc * real((-1) ** ((ia + ib - ic) / 2), 8)
      id = -id
      ie = -ie
      ifv = -ifv
    end if
    
    ! Calculate factorials and prefactor
    fc2 = real(ic + 1, 8)
    iabcp = (ia + ib + ic) / 2 + 1
    iabc = iabcp - ic
    icab = iabcp - ib
    ibca = iabcp - ia
    iapd = (ia + id) / 2 + 1
    iamd = iapd - id
    ibpe = (ib + ie) / 2 + 1
    ibme = ibpe - ie
    icpf = (ic + ifv) / 2 + 1
    icmf = icpf - ifv
    
    sqfclg = 0.5d0 * (dlog(fc2) - faclog(iabcp + 1) + &
                      faclog(iabc) + faclog(icab) + faclog(ibca) + &
                      faclog(iapd) + faclog(iamd) + faclog(ibpe) + &
                      faclog(ibme) + faclog(icpf) + faclog(icmf))
    
    nzmic2 = (ib - ic - id) / 2
    nzmic3 = (ia - ic + ie) / 2
    nzmi = max(0, nzmic2, nzmic3) + 1
    nzmx = min(iabc, iamd, ibpe)
    
    if (nzmi > nzmx) return
    
    ! Sum over allowed nz values
    s = 0.0d0
    s1 = real((-1) ** (nzmi - 1), 8)
    
    do nz = nzmi, nzmx
      nzm1 = nz - 1
      nzt1 = iabc - nzm1
      nzt2 = iamd - nzm1
      nzt3 = ibpe - nzm1
      nzt4 = nz - nzmic2
      nzt5 = nz - nzmic3
      
      termlg = sqfclg - faclog(nz) - faclog(nzt1) - faclog(nzt2) - &
               faclog(nzt3) - faclog(nzt4) - faclog(nzt5)
      ssterm = s1 * dexp(termlg)
      s = s + ssterm
      s1 = -s1
    end do
    
    rac = signfc * s
    
  end function clebsch_gordan_calc

end module j_coupling_symbols
!!========================================
!subroutine test_3j6j9j_comparison()
!  !-------------------------------------------------------------------
!  ! Comprehensive test program for new 3j6j9jmod.f90
!  ! Tests with physically valid angular momentum coupling
!  ! (Triangle inequality satisfied, M conservation satisfied)
!  !-------------------------------------------------------------------
!  use j_coupling_symbols
!  implicit none
!  
!  real(kind=8) :: result
!  integer :: i, ntest_3j, ntest_6j, ntest_9j
!  real(kind=8), allocatable :: cases_3j(:,:), cases_6j(:,:), cases_9j(:,:)
!  
!  !===================================================================
!  ! 3j SYMBOL TESTS: (j1, m1; j2, m2; j3, m3)
!  ! Must satisfy: |j1-j2| <= j3 <= j1+j2, m1+m2+m3=0
!  !===================================================================
!  ntest_3j = 15
!  allocate(cases_3j(ntest_3j, 6))
!  
!  ! Row format: j1, m1, j2, m2, j3, m3 (for 3j symbol)
!  ! Case 1: (1/2, 1/2; 1/2, -1/2; 0, 0) -> Non-zero
!  cases_3j(1, :) = [0.5d0, 0.5d0, 0.5d0, -0.5d0, 0.0d0, 0.0d0]
!  
!  ! Case 2: (1/2, 1/2; 1/2, 1/2; 1, -1) -> Non-zero
!  cases_3j(2, :) = [0.5d0, 0.5d0, 0.5d0, 0.5d0, 1.0d0, -1.0d0]
!  
!  ! Case 3: (1, 1; 1, -1; 0, 0) -> Non-zero
!  cases_3j(3, :) = [1.0d0, 0.0d0, 1.0d0, 0.0d0, 0.0d0, 0.0d0]
!  
!  ! Case 4: (1, 1; 1, -1; 0, 0)
!  cases_3j(4, :) = [1.0d0, 1.0d0, 1.0d0, -1.0d0, 0.0d0, 0.0d0]
!  
!  ! Case 5: (1, 0; 1, 0; 2, 0)
!  cases_3j(5, :) = [1.0d0, 0.0d0, 1.0d0, 0.0d0, 2.0d0, 0.0d0]
!  
!  ! Case 6: (1, 1; 1, 0; 2, -1)
!  cases_3j(6, :) = [1.0d0, 1.0d0, 1.0d0, 0.0d0, 2.0d0, -1.0d0]
!  
!  ! Case 7: (1/2, 1/2; 1, -1/2; 1/2, 0)
!  cases_3j(7, :) = [0.5d0, 0.5d0, 1.0d0, -0.5d0, 0.5d0, 0.0d0]
!  
!  ! Case 8: (1, 1; 1/2, -1/2; 3/2, -0.5)
!  cases_3j(8, :) = [1.0d0, 1.0d0, 0.5d0, -0.5d0, 1.5d0, -0.5d0]
!  
!  ! Case 9: (3/2, 3/2; 1/2, -1/2; 1, -1.0)
!  cases_3j(9, :) = [1.5d0, 1.5d0, 0.5d0, -0.5d0, 1.0d0, -1.0d0]
!  
!  ! Case 10: (2, 2; 1, -1; 1, -1.0)
!  cases_3j(10, :) = [2.0d0, 2.0d0, 1.0d0, -1.0d0, 1.0d0, -1.0d0]
!  
!  ! Case 11: (3/2, 1/2; 1/2, 1/2; 1, -1)
!  cases_3j(11, :) = [1.5d0, 0.5d0, 0.5d0, 0.5d0, 1.0d0, -1.0d0]
!  
!  ! Case 12: (1, -1; 1, 1; 0, 0)
!  cases_3j(12, :) = [1.0d0, -1.0d0, 1.0d0, 1.0d0, 0.0d0, 0.0d0]
!  
!  ! Case 13: (2, 0; 1, 0; 1, 0)
!  cases_3j(13, :) = [2.0d0, 0.0d0, 1.0d0, 0.0d0, 1.0d0, 0.0d0]
!  
!  ! Case 14: (3/2, -1/2; 1, 1/2; 1/2, 0)
!  cases_3j(14, :) = [1.5d0, -0.5d0, 1.0d0, 0.5d0, 0.5d0, 0.0d0]
!  
!  ! Case 15: (2, 1; 1, -1; 1, 0)
!  cases_3j(15, :) = [2.0d0, 1.0d0, 1.0d0, -1.0d0, 1.0d0, 0.0d0]
!  
!  !===================================================================
!  ! 6j SYMBOL TESTS: {j1, j2, j3; j4, j5, j6}
!  ! Each triplet must satisfy triangle inequality
!  !===================================================================
!  ntest_6j = 12
!  allocate(cases_6j(ntest_6j, 6))
!  
!  ! Row format: j1, j2, j3, j4, j5, j6 (for 6j symbol)
!  ! Case 1: {1, 1, 1; 1, 1, 1} -> Standard non-zero case
!  cases_6j(1, :) = [1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0]
!  
!  ! Case 2: {1/2, 1/2, 1; 1/2, 1/2, 0} -> Non-zero
!  cases_6j(2, :) = [0.5d0, 0.5d0, 1.0d0, 0.5d0, 0.5d0, 0.0d0]
!  
!  ! Case 3: {1, 1, 2; 1, 1, 2}
!  cases_6j(3, :) = [1.0d0, 1.0d0, 2.0d0, 1.0d0, 1.0d0, 2.0d0]
!  
!  ! Case 4: {1/2, 1/2, 0; 1/2, 1/2, 0} -> Non-zero
!  cases_6j(4, :) = [0.5d0, 0.5d0, 0.0d0, 0.5d0, 0.5d0, 0.0d0]
!  
!  ! Case 5: {1, 1, 1; 1, 1, 1}
!  cases_6j(5, :) = [1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0]
!  
!  ! Case 6: {2, 2, 2; 2, 2, 2}
!  cases_6j(6, :) = [2.0d0, 2.0d0, 2.0d0, 2.0d0, 2.0d0, 2.0d0]
!  
!  ! Case 7: {1, 2, 1; 1, 2, 1}
!  cases_6j(7, :) = [1.0d0, 2.0d0, 1.0d0, 1.0d0, 2.0d0, 1.0d0]
!  
!  ! Case 8: {1.5, 1.5, 1; 1.5, 1.5, 1}
!  cases_6j(8, :) = [1.5d0, 1.5d0, 1.0d0, 1.5d0, 1.5d0, 1.0d0]
!  
!  ! Case 9: {3/2, 3/2, 1; 3/2, 3/2, 1}
!  cases_6j(9, :) = [1.5d0, 1.5d0, 1.0d0, 1.5d0, 1.5d0, 1.0d0]
!  
!  ! Case 10: {1, 2, 2; 2, 1, 2}
!  cases_6j(10, :) = [1.0d0, 2.0d0, 2.0d0, 2.0d0, 1.0d0, 2.0d0]
!  
!  ! Case 11: {1/2, 3/2, 1; 3/2, 1/2, 1}
!  cases_6j(11, :) = [0.5d0, 1.5d0, 1.0d0, 1.5d0, 0.5d0, 1.0d0]
!  
!  ! Case 12: {2, 2, 0; 2, 2, 0}
!  cases_6j(12, :) = [2.0d0, 2.0d0, 0.0d0, 2.0d0, 2.0d0, 0.0d0]
!  
!  !===================================================================
!  ! 9j SYMBOL TESTS (from coef9j function mapping)
!  ! j values arranged as:
!  ! j1  j2  j3     -> j1, j2, j4, j5, j3, j6, j7, j8, j9 in coef9j
!  ! j4  j5  j6
!  ! j7  j8  j9
!  !===================================================================
!  ntest_9j = 10
!  allocate(cases_9j(ntest_9j, 9))
!  
!  ! Physical valid 9j cases (constructed carefully)
!  ! Case 1: All j=1 (should be non-zero)
!  cases_9j(1, :) = [1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0]
!  
!  ! Case 2: All j=1/2, 1 (should be non-zero)
!  cases_9j(2, :) = [0.5d0, 0.5d0, 1.0d0, 0.5d0, 0.5d0, 1.0d0, 0.5d0, 0.5d0, 1.0d0]
!  
!  ! Case 3: j9=0 identity case (reduces to 6j, should be non-zero)
!  ! { j1 j2 j3 }
!  ! { j4 j5 j3 }
!  ! { j7 j7 0  }
!  cases_9j(3, :) = [1.0d0, 1.0d0, 0.0d0, 1.0d0, 1.0d0, 0.0d0, 1.0d0, 1.0d0, 0.0d0]
!  
!  ! Case 4: Another j9=0 identity (all 1/2s, 0s)
!  cases_9j(4, :) = [0.5d0, 0.5d0, 0.0d0, 0.5d0, 0.5d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0]
!  
!  ! Case 5: j3=0, j6=0, j9=0 (reduces to 6j, should be non-zero)
!  ! { j1 j1 0 }
!  ! { j4 j4 0 }
!  ! { j7 j7 0 }
!  cases_9j(5, :) = [1.0d0, 1.0d0, 0.0d0, 2.0d0, 2.0d0, 0.0d0, 1.0d0, 1.0d0, 0.0d0]
!  
!  ! Case 6: Mixed integers, half-integers (should be non-zero)
!  cases_9j(6, :) = [0.5d0, 0.5d0, 1.0d0, 1.0d0, 1.0d0, 0.0d0, 0.5d0, 0.5d0, 1.0d0]
!  
!  ! Case 7: Larger angular momenta (should be non-zero)
!  cases_9j(7, :) = [2.0d0, 1.0d0, 1.0d0, 1.0d0, 2.0d0, 1.0d0, 1.0d0, 2.0d0, 1.0d0]
!  
!  ! Case 8: From original, should be non-zero
!  cases_9j(8, :) = [0.5d0, 0.5d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0, 0.0d0]
!  
!  ! Case 9: From original, should be non-zero
!  cases_9j(9, :) = [1.0d0, 1.0d0, 2.0d0, 1.0d0, 1.0d0, 2.0d0, 1.0d0, 1.0d0, 2.0d0]
!  
!  ! Case 10: From original, should be non-zero
!  cases_9j(10, :) = [1.0d0, 0.5d0, 1.5d0, 0.5d0, 1.0d0, 1.5d0, 0.5d0, 1.5d0, 1.0d0]
!  
!  !===================================================================
!  ! TEST 3j SYMBOLS
!  !===================================================================
!  write(*, '(A)') ''
!  write(*, '(A)') '========================================='
!  write(*, '(A)') '3j SYMBOL TESTS (tj2i function)'
!  write(*, '(A)') 'Valid triangle inequality & M conservation'
!  write(*, '(A)') '========================================='
!  write(*, '(A)') ''
!  
!  do i = 1, ntest_3j
!    ! Fix: Passed in order (j1, j2, j3, m1, m2, m3) for tj2i
!    result = tj2i(2.0d0*cases_3j(i,1), 2.0d0*cases_3j(i,3), &
!                  2.0d0*cases_3j(i,5), 2.0d0*cases_3j(i,2), &
!                  2.0d0*cases_3j(i,4), 2.0d0*cases_3j(i,6))
!    
!    write(*, '(A,I2,A)') 'Test 3j-', i, ': (j1,m1; j2,m2; j3,m3) = '
!    write(*, '(A,3(F4.1,A,F5.1,A))') &
!      '  (', cases_3j(i,1), ',', cases_3j(i,2), '; ', &
!           cases_3j(i,3), ',', cases_3j(i,4), '; ', &
!           cases_3j(i,5), ',', cases_3j(i,6), ')'
!    write(*, '(A,F18.12)') '  Result: ', result
!    
!    if (abs(result) < 1.0d-14) then
!      write(*, '(A)') '  Status: ZERO'
!    else
!      write(*, '(A,E15.8)') '  Status: NON-ZERO', result
!    end if
!    write(*, '(A)') ''
!  end do
!  
!  !===================================================================
!  ! TEST 6j SYMBOLS
!  !===================================================================
!  write(*, '(A)') '========================================='
!  write(*, '(A)') '6j SYMBOL TESTS (sj2i function)'
!  write(*, '(A)') 'Format: {j1,j2,j3; j4,j5,j6}'
!  write(*, '(A)') '========================================='
!  write(*, '(A)') ''
!  
!  do i = 1, ntest_6j
!    ! For sj2i, arguments are 2*j values
!    result = sj2i(2.0d0*cases_6j(i,1), 2.0d0*cases_6j(i,2), &
!                  2.0d0*cases_6j(i,3), 2.0d0*cases_6j(i,4), &
!                  2.0d0*cases_6j(i,5), 2.0d0*cases_6j(i,6))
!    
!    write(*, '(A,I2,A)') 'Test 6j-', i, ':'
!    write(*, '(A,6F4.1,A)') '  {', cases_6j(i,:), '}'
!    write(*, '(A,F18.12)') '  Result: ', result
!    
!    if (abs(result) < 1.0d-14) then
!      write(*, '(A)') '  Status: ZERO'
!    else
!      write(*, '(A,E15.8)') '  Status: NON-ZERO', result
!    end if
!    write(*, '(A)') ''
!  end do
!  
!  !===================================================================
!  ! TEST 9j SYMBOLS
!  !===================================================================
!  write(*, '(A)') '========================================='
!  write(*, '(A)') '9j SYMBOL TESTS (coef9j function)'
!  write(*, '(A)') 'Arrangement: j1 j2 j3 / j4 j5 j6 / j7 j8 j9'
!  write(*, '(A)') '========================================='
!  write(*, '(A)') ''
!  
!  do i = 1, ntest_9j
!    result = coef9j(cases_9j(i,1), cases_9j(i,2), cases_9j(i,3), &
!                    cases_9j(i,4), cases_9j(i,5), cases_9j(i,6), &
!                    cases_9j(i,7), cases_9j(i,8), cases_9j(i,9))
!    
!    write(*, '(A,I2,A)') 'Test 9j-', i, ':'
!    write(*, '(A,3F4.1)') '  Row 1: ', cases_9j(i,1:3)
!    write(*, '(A,3F4.1)') '  Row 2: ', cases_9j(i,4:6)
!    write(*, '(A,3F4.1)') '  Row 3: ', cases_9j(i,7:9)
!    write(*, '(A,F18.12)') '  Result: ', result
!    
!    if (abs(result) < 1.0d-14) then
!      write(*, '(A)') '  Status: ZERO'
!    else
!      write(*, '(A,E15.8)') '  Status: NON-ZERO', result
!    end if
!    write(*, '(A)') ''
!  end do
!  
!  write(*, '(A)') '========================================='
!  write(*, '(A)') 'All comprehensive tests completed'
!  write(*, '(A)') '========================================='
!  write(*, '(A)') ''
!  
!  deallocate(cases_3j, cases_6j, cases_9j)
!    
!end subroutine test_3j6j9j_comparison
!!------------------------------------
!program main
!  implicit none
!
!  !-------------------------------------------------------------------
!  ! Test comparison between old and new 3j6j9j implementations
!  !-------------------------------------------------------------------
!  call test_3j6j9j_comparison()
!end program main