module yaeos__models_ar_cubic_mixing_sddlc
   use yaeos__constants, only: pr, R
   use yaeos__models_ar_cubic_quadratic_mixing, only: QMRTD

   type, extends(QMRTD) :: sDDLC
      !! Segmented Density-Dependent Local-Composition Mixing Rule.
      real(pr), allocatable :: q(:) !! Segment size
   contains
      procedure :: Dmix => ddlc_Dmix
   end type sDDLC

   interface sDDLC
      module procedure :: init_sddlc
   end interface sDDLC

contains

   type(sDDLC) function init_sddlc(k, k0, tref, l, q) result(mixrule)
      real(pr), intent(in) :: k(:, :)
      real(pr), intent(in) :: k0(:, :)
      real(pr), intent(in) :: tref(:, :)
      real(pr), intent(in) :: l(:, :)
      real(pr), intent(in) :: q(:)

      mixrule%k = k
      mixrule%k0 = k0
      mixrule%tref =  tref
      mixrule%l =  l
      mixrule%q =  q
      mixrule%is_D_ddlc = .true.
   end function init_sddlc

   subroutine ddlc_Dmix(&
      self, n, V, T, &
      ai, daidt, daidt2, &
      D, &
      dDdV, dDdT, dDdV2, dDdT2, dDi, dDdTV, dDidV, dDidT, dDij &
      )
      !! s-DDLC D mixing rule including V, n, and T derivatives.
      !!
      !! The local-composition D is:
      !! \[
      !!   D = \frac{\sum_{ij} n_i n_j E_{ij}(V,T)\,a_{ij}(T)}{\Lambda(V,T,n)}
      !! \]
      !! where \(E_{ij} = \exp\!\left(\frac{a_{ij}\sum_k n_k q_k}{q_i q_j RTV}\right)\)
      !! and \(\Lambda = \sum_{ij} s_i s_j E_{ij}\), with segment fractions
      !! \(s_i = n_i q_i / \sum_k n_k q_k\).
      class(sDDLC), intent(in)  :: self
      real(pr),            intent(in)  :: n(:), V, T
      real(pr),            intent(in)  :: ai(:), daidt(:), daidt2(:)
      real(pr),            intent(out) :: D, dDdT, dDdT2, dDdV, dDdV2, dDdTV
      real(pr),            intent(out) :: dDi(:), dDidT(:), dDidV(:), dDij(:,:)

      integer  :: i, j, k, m, nc
      real(pr) :: RTV, sumnq

      ! -- aij and its T-derivatives --
      real(pr) :: aij(size(n),size(n))
      real(pr) :: daijdT(size(n),size(n)), daijdT2(size(n),size(n))

      ! -- Segment fractions and composition derivatives --
      real(pr) :: sfrac(size(n))
      real(pr) :: dsfdn(size(n),size(n))
      real(pr) :: d2sfdn2(size(n),size(n),size(n))

      ! -- E matrix and all its derivatives --
      real(pr) :: rata(size(n),size(n))
      real(pr) :: E(size(n),size(n))
      real(pr) :: dEdn(size(n),size(n),size(n))
      real(pr) :: d2Edn2(size(n),size(n),size(n),size(n))
      real(pr) :: dEdT(size(n),size(n)), d2EdT2(size(n),size(n))
      real(pr) :: dEdV(size(n),size(n)), d2EdV2(size(n),size(n))
      real(pr) :: d2EdnT(size(n),size(n),size(n))
      real(pr) :: d2EdnV(size(n),size(n),size(n))
      real(pr) :: d2EdVT(size(n),size(n))

      ! -- Denominator Den and its derivatives --
      real(pr) :: Den
      real(pr) :: dDendn(size(n)), d2Dendn2(size(n),size(n))
      real(pr) :: d2DendnT(size(n)), d2DendnV(size(n))
      real(pr) :: dDendV, d2DendV2, d2DendVT, dDendT, d2DendT2

      ! -- Loop temporaries --
      real(pr) :: aux, auxV, auxVV, auxT, auxTT, auxVT
      real(pr) :: auxjmi, auxjmiT, auxjmiV
      real(pr) :: auxdfE, auxd2fE, auxdfdE, auxfd2E, auxd2Ek
      real(pr) :: auxTij, auxdE, auxd2E

      nc  = size(n)
      RTV = R*T*V

      ! -----------------------------------------------------------------------
      ! Build aij matrix (reusing the shared helper)
      ! -----------------------------------------------------------------------
      call self%aij(T, ai, daidt, daidt2, aij, daijdT, daijdT2)

      ! -----------------------------------------------------------------------
      ! Segment fractions and their n-derivatives
      ! -----------------------------------------------------------------------
      sumnq = sum(n*self%q)

      do i = 1, nc
         sfrac(i) = n(i)*self%q(i)/sumnq
         do j = 1, nc
            dsfdn(i,j) = -sfrac(i)*self%q(j)/sumnq
            if (i == j) dsfdn(i,i) = dsfdn(i,i) + self%q(i)/sumnq
            do k = j, nc
               d2sfdn2(i,j,k) = 2.0_pr*sfrac(i)*self%q(j)*self%q(k)/sumnq**2
               if (i == j) d2sfdn2(i,i,k) = d2sfdn2(i,i,k) &
                  - self%q(i)*self%q(k)/sumnq**2
               if (i == k) d2sfdn2(i,j,i) = d2sfdn2(i,j,i) &
                  - self%q(i)*self%q(j)/sumnq**2
               d2sfdn2(i,k,j) = d2sfdn2(i,j,k)
            end do
         end do
      end do

      ! -----------------------------------------------------------------------
      ! E matrix, V/n/T-derivatives, and accumulation of Den
      ! -----------------------------------------------------------------------
      Den        = 0.0_pr
      dDendn     = 0.0_pr
      d2DendnT   = 0.0_pr
      d2DendnV   = 0.0_pr
      dDendV     = 0.0_pr
      d2DendV2   = 0.0_pr
      d2DendVT   = 0.0_pr
      dDendT     = 0.0_pr
      d2DendT2   = 0.0_pr
      d2Edn2     = 0.0_pr
      dEdT       = 0.0_pr
      d2EdT2     = 0.0_pr
      d2EdVT     = 0.0_pr
      d2EdnT     = 0.0_pr

      do k = 1, nc
         do m = k, nc
            rata(k,m) = aij(k,m)/(self%q(k)*self%q(m)*RTV)
            E(k,m)    = exp(rata(k,m)*sumnq)
            E(m,k)    = E(k,m)

            ! -- n-derivatives of E --
            do i = 1, nc
               dEdn(k,m,i)     = E(k,m)*rata(k,m)*self%q(i)
               dEdn(m,k,i)     = dEdn(k,m,i)
               d2Edn2(k,m,i,i) = dEdn(k,m,i)**2/E(k,m)
               d2Edn2(m,k,i,i) = d2Edn2(k,m,i,i)
               do j = 1, i-1
                  d2Edn2(k,m,i,j) = dEdn(k,m,i)*dEdn(k,m,j)/E(k,m)
                  d2Edn2(k,m,j,i) = d2Edn2(k,m,i,j)
                  d2Edn2(m,k,i,j) = d2Edn2(k,m,i,j)
                  d2Edn2(m,k,j,i) = d2Edn2(k,m,i,j)
               end do
            end do

            ! -- V-derivatives of E --
            dEdV(k,m)     = -E(k,m)*rata(k,m)*sumnq/V
            d2EdV2(k,m)   =  dEdV(k,m)**2/E(k,m) - 2.0_pr*dEdV(k,m)/V
            d2EdnV(k,m,:) =  self%q(:)*dEdV(k,m)/sumnq &
               + dEdV(k,m)*dEdn(k,m,:)/E(k,m)
            dEdV(m,k)     = dEdV(k,m)
            d2EdV2(m,k)   = d2EdV2(k,m)
            d2EdnV(m,k,:) = d2EdnV(k,m,:)

            ! -- T-derivatives of E (computed only when requested) --
            dEdT(k,m)    = E(k,m)*rata(k,m)*sumnq &
               *(daijdT(k,m)/aij(k,m) - 1.0_pr/T)
            d2EdT2(k,m)  = dEdT(k,m)**2/E(k,m) - 2.0_pr*dEdT(k,m)/T &
               + E(k,m)*rata(k,m)*sumnq*daijdT2(k,m)/aij(k,m)
            d2EdVT(k,m)  = dEdV(k,m)* &
               (dEdT(k,m)/E(k,m) + daijdT(k,m)/aij(k,m) - 1.0_pr/T)
            d2EdnT(k,m,:)= self%q(:)*dEdT(k,m)/sumnq &
               + dEdT(k,m)*dEdn(k,m,:)/E(k,m)
            dEdT(m,k)    = dEdT(k,m)
            d2EdT2(m,k)  = d2EdT2(k,m)
            d2EdVT(m,k)  = d2EdVT(k,m)
            d2EdnT(m,k,:)= d2EdnT(k,m,:)

            ! -- Accumulate Den and its derivatives --
            if (k == m) then
               Den      = Den      + sfrac(k)**2*E(k,k)
               dDendV   = dDendV   + sfrac(k)**2*dEdV(k,k)
               d2DendV2 = d2DendV2 + sfrac(k)**2*d2EdV2(k,k)
               dDendT   = dDendT   + sfrac(k)**2*dEdT(k,k)
               d2DendT2 = d2DendT2 + sfrac(k)**2*d2EdT2(k,k)
               d2DendVT = d2DendVT + sfrac(k)**2*d2EdVT(k,k)
               do i = 1, nc
                  dDendn(i)   = dDendn(i) &
                     + 2.0_pr*sfrac(k)*dsfdn(k,i)*E(k,k) &
                     + sfrac(k)**2*dEdn(k,k,i)
                  d2DendnV(i) = d2DendnV(i) &
                     + 2.0_pr*sfrac(k)*dsfdn(k,i)*dEdV(k,k) &
                     + sfrac(k)**2*d2EdnV(k,k,i)
                  d2DendnT(i) = d2DendnT(i) &
                     + 2.0_pr*sfrac(k)*dsfdn(k,i)*dEdT(k,k) &
                     + sfrac(k)**2*d2EdnT(k,k,i)
               end do
            else
               Den      = Den      + 2.0_pr*sfrac(k)*sfrac(m)*E(k,m)
               dDendV   = dDendV   + 2.0_pr*sfrac(k)*sfrac(m)*dEdV(k,m)
               d2DendV2 = d2DendV2 + 2.0_pr*sfrac(k)*sfrac(m)*d2EdV2(k,m)
               dDendT   = dDendT   + 2.0_pr*sfrac(k)*sfrac(m)*dEdT(k,m)
               d2DendT2 = d2DendT2 + 2.0_pr*sfrac(k)*sfrac(m)*d2EdT2(k,m)
               d2DendVT = d2DendVT + 2.0_pr*sfrac(k)*sfrac(m)*d2EdVT(k,m)
               do i = 1, nc
                  dDendn(i)   = dDendn(i) &
                     + 2.0_pr*sfrac(k)*sfrac(m)*dEdn(k,m,i) &
                     + 2.0_pr*(sfrac(k)*dsfdn(m,i) + sfrac(m)*dsfdn(k,i))*E(k,m)
                  d2DendnV(i) = d2DendnV(i) &
                     + 2.0_pr*sfrac(k)*sfrac(m)*d2EdnV(k,m,i) &
                     + 2.0_pr*(sfrac(k)*dsfdn(m,i)+sfrac(m)*dsfdn(k,i))*dEdV(k,m)
                  d2DendnT(i) = d2DendnT(i) &
                     + 2.0_pr*sfrac(k)*sfrac(m)*d2EdnT(k,m,i) &
                     + 2.0_pr*(sfrac(k)*dsfdn(m,i)+sfrac(m)*dsfdn(k,i))*dEdT(k,m)
               end do
            end if
         end do
      end do

      ! -----------------------------------------------------------------------
      ! Second composition derivative of Den
      ! -----------------------------------------------------------------------
      d2Dendn2 = 0.0_pr
      do i = 1, nc
         do j = i, nc
            do k = 1, nc
               auxdfE  = 0.0_pr
               auxd2fE = 0.0_pr
               auxdfdE = 0.0_pr
               auxfd2E = 0.0_pr
               do m = 1, nc
                  auxdfE  = auxdfE  + dsfdn(m,i)*E(k,m)
                  auxd2fE = auxd2fE + d2sfdn2(m,i,j)*E(k,m)
                  auxdfdE = auxdfdE &
                     + dsfdn(m,i)*dEdn(k,m,j) + dsfdn(m,j)*dEdn(k,m,i)
                  auxfd2E = auxfd2E + sfrac(m)*d2Edn2(k,m,i,j)
               end do
               d2Dendn2(i,j) = d2Dendn2(i,j) &
                  + 2.0_pr*dsfdn(k,j)*auxdfE &
                  + sfrac(k)*(2.0_pr*auxd2fE + 2.0_pr*auxdfdE + auxfd2E)
            end do
            d2Dendn2(j,i) = d2Dendn2(i,j)
         end do
      end do

      ! -----------------------------------------------------------------------
      ! Accumulate D and its n/V/T derivatives (numerator terms, before /Den)
      ! -----------------------------------------------------------------------
      D      = 0.0_pr
      dDdT   = 0.0_pr
      dDdT2  = 0.0_pr
      dDdV     = 0.0_pr
      dDdV2    = 0.0_pr
      dDdTV    = 0.0_pr
      dDi    = 0.0_pr
      dDidV    = 0.0_pr
      dDidT  = 0.0_pr

      do i = 1, nc
         aux   = 0.0_pr; auxV  = 0.0_pr; auxVV = 0.0_pr
         auxT  = 0.0_pr; auxTT = 0.0_pr; auxVT = 0.0_pr
         do j = 1, nc
            auxjmi  = 0.0_pr; auxjmiT = 0.0_pr; auxjmiV = 0.0_pr
            do m = 1, nc
               auxjmi  = auxjmi  + n(m)*dEdn(m,j,i)*aij(m,j)
               auxjmiV = auxjmiV + n(m)*d2EdnV(m,j,i)*aij(m,j)
               auxjmiT = auxjmiT + n(m)*( &
                  dEdn(m,j,i)*daijdT(m,j) + d2EdnT(m,j,i)*aij(m,j))
            end do

            dDi(i) = dDi(i) + 2.0_pr*n(j)*E(i,j)*aij(i,j) + n(j)*auxjmi
            dDidV(i) = dDidV(i) + 2.0_pr*n(j)*dEdV(i,j)*aij(i,j) + n(j)*auxjmiV

            auxTij  = E(i,j)*daijdT(i,j) + dEdT(i,j)*aij(i,j)
            auxT    = auxT + n(j)*auxTij
            dDidT(i)= dDidT(i) + 2.0_pr*n(j)*auxTij + n(j)*auxjmiT
            auxTT   = auxTT + n(j)*( &
               E(i,j)*daijdT2(i,j) + 2.0_pr*dEdT(i,j)*daijdT(i,j) &
               + d2EdT2(i,j)*aij(i,j))
            auxVT   = auxVT + n(j)*( &
               d2EdVT(i,j)*aij(i,j) + dEdV(i,j)*daijdT(i,j))

            aux   = aux   + n(j)*E(i,j)*aij(i,j)
            auxV  = auxV  + n(j)*dEdV(i,j)*aij(i,j)
            auxVV = auxVV + n(j)*d2EdV2(i,j)*aij(i,j)
         end do

         D   = D   + n(i)*aux
         dDdV  = dDdV  + n(i)*auxV
         dDdV2 = dDdV2 + n(i)*auxVV
         dDdT  = dDdT  + n(i)*auxT
         dDdT2 = dDdT2 + n(i)*auxTT
         dDdTV   = dDdTV   + n(i)*auxVT
      end do

      ! -----------------------------------------------------------------------
      ! Complete D and its derivatives (divide by Den, apply quotient rule)
      ! -----------------------------------------------------------------------
      D   =  D/Den
      dDdV  = (dDdV  - D*dDendV)/Den
      dDdV2 = (dDdV2 - 2.0_pr*dDdV*dDendV - D*d2DendV2)/Den

      dDi(1:nc) = (dDi(1:nc)  - D*dDendn(1:nc))/Den
      dDidV(1:nc) = (dDidV(1:nc)  - dDdV*dDendn(1:nc) - D*d2DendnV(1:nc) &
         - dDi(1:nc)*dDendV)/Den

      dDdT  = (dDdT  - D*dDendT)/Den
      dDdT2 = (dDdT2 - 2.0_pr*dDdT*dDendT - D*d2DendT2)/Den
      dDdTV   = (dDdTV   - dDdT*dDendV - dDdV*dDendT - D*d2DendVT)/Den
      dDidT(1:nc) = (dDidT(1:nc) - dDdT*dDendn(1:nc) - D*d2DendnT(1:nc) &
         - dDi(1:nc)*dDendT)/Den

      ! -----------------------------------------------------------------------
      ! Second composition derivative of D (= Dij)
      ! -----------------------------------------------------------------------
      do i = 1, nc
         do j = i, nc
            auxdE  = 0.0_pr
            auxd2E = 0.0_pr
            do k = 1, nc
               auxd2Ek = 0.0_pr
               do m = 1, nc
                  auxd2Ek = auxd2Ek + n(m)*d2Edn2(k,m,i,j)*aij(k,m)
               end do
               auxdE  = auxdE + n(k)*( &
                  dEdn(i,k,j)*aij(i,k) + dEdn(k,j,i)*aij(k,j))
               auxd2E = auxd2E + n(k)*auxd2Ek
            end do
            dDij(i,j) = (2.0_pr*E(i,j)*aij(i,j) + 2.0_pr*auxdE + auxd2E &
               - dDi(j)*dDendn(i) - dDi(i)*dDendn(j) &
               - D*d2Dendn2(i,j)) / Den
            dDij(j,i) = dDij(i,j)
         end do
      end do
   end subroutine ddlc_Dmix
end module yaeos__models_ar_cubic_mixing_sddlc

