Option Explicit 'Para que sea necesario declarar todas las variables y no quede a interpretación del VBA. Acelera el programa
Option Base 1 'Para que la primera posición de la matriz sea 1
'No considera el equilibrio
'Los gases a la salida son ideales
'El intervalo del cp del gas ideal es menor al valor calculado de Tf
'Mezcla de CH4 y C2H6 como combustible

Public Function fu(Te, T0, T, A, B, Ns, Ne, Hr, R) As Double
    'Balance de energía
    
    
    Dim tau, Mcpe, Mcps, tem  As Double
    
    
    
    tau = Te / T0
    
    tem = A(1) _
    + A(2) / 2 * T0 * (tau + 1) _
    + A(3) / 3 * T0 * T0 * (tau * tau + tau + 1) _
    + A(4) / 4 * T0 * T0 * T0 * (tau * tau * tau + tau * tau + tau + 1) _
    + A(5) / 5 * T0 * T0 * T0 * T0 * (tau * tau * tau * tau + tau * tau * tau + tau * tau + tau + 1)
    
    Mcpe = R * tem
    
    tau = T / T0
    tem = B(1) _
    + B(2) / 2 * T0 * (tau + 1) _
    + B(3) / 3 * T0 * T0 * (tau * tau + tau + 1) _
    + B(4) / 4 * T0 * T0 * T0 * (tau * tau * tau + tau * tau + tau + 1) _
    + B(5) / 5 * T0 * T0 * T0 * T0 * (tau * tau * tau * tau + tau * tau * tau + tau * tau + tau + 1)
    
    Mcps = R * tem
    
    fu = Hr + (T - T0) * (Ns * Mcps) - Ne * Mcpe * (Te - T0)

End Function


Function Quemador(ID As Integer, Te As Double, Ne As Double, Conversion As Range, _
                  Tabye As Range, Tabcp As Range, Dhf As Range, Tabalpha As Range) As Double

    Dim n, m, nR, i, j As Integer
    Dim cp(), ye(), A(), Xi(), alpha(), Ns, Reac, ys(), B() As Double
    Dim Hr, He, Hs, Dhr, T, T0, tau, tem, Mcpe, Mcps, Error, Qmax As Double
    Dim Indice As Long
    Dim R As Double

    nR = Tabalpha.Columns.Count 'número de reacciones
    n = Tabcp.Rows.Count 'numero reglones = numero de componentes
    m = Tabcp.Columns.Count 'numero de columnas = número de constantes

    ReDim ye(n), ys(n), Ccp(n, m), A(m), B(m) 'son m=5 consntante del cp
    ReDim Xi(nR), alpha(n, nR), E(nR)


    T0 = 298.15 'K temperatura de referencia

    'Asignación de valores

    For i = 1 To n
        ye(i) = Tabye.Cells(i)
        For j = 1 To m
            Ccp(i, j) = Tabcp.Cells(i, j)
        Next j
    Next i



    For j = 1 To nR
        For i = 1 To n
            alpha(i, j) = Tabalpha.Cells(i, j)
        Next i
    Next j

    Xi(1) = -Conversion.Cells(1) / 100 * ye(1) * Ne / alpha(1, 1)
    Xi(2) = -Conversion.Cells(2) / 100 * ye(2) * Ne / alpha(2, 2)



    'Valores reales de los coeficientes

    For j = 1 To m    'sobre las constantes del cp
        A(j) = 0
        For i = 1 To n    'sobre los componentes
            A(j) = A(j) + ye(i) * Ccp(i, j)
        Next i
    Next j

    A(1) = A(1)
    A(2) = A(2) * 10 ^ (-3)
    A(3) = A(3) * 10 ^ (-5)
    A(4) = A(4) * 10 ^ (-8)
    A(5) = A(5) * 10 ^ (-11)

    'Inicio Balance de materia


    'Flujo total a la salida

    Ns = 0
    For i = 1 To n
        Reac = 0
        For j = 1 To nR
            Reac = Reac + alpha(i, j) * Xi(j)
        Next j
        Ns = Ns + Ne * ye(i) + Reac
    Next i

    'Fracción molar a la salida

    For i = 1 To n
        Reac = 0
        For j = 1 To nR
            Reac = Reac + alpha(i, j) * Xi(j)
        Next j
        ys(i) = (Ne * ye(i) + Reac) / Ns
    Next i

    'Fin Balance de materia


    For j = 1 To m 'sobre las constantes del cp
        B(j) = 0
        For i = 1 To n 'sobre los componentes
            B(j) = B(j) + ys(i) * Ccp(i, j)
        Next i
    Next j

    'Valores reales de los coeficientes
    B(2) = B(2) * 10 ^ (-3)
    B(3) = B(3) * 10 ^ (-5)
    B(4) = B(4) * 10 ^ (-8)
    B(5) = B(5) * 10 ^ (-11)


    ' Inicio Balance de Energía

    Hr = 0
    For j = 1 To nR
        Dhr = 0
        For i = 1 To n
            Dhr = Dhr + alpha(i, j) * Dhf.Cells(i)
        Next i
        Hr = Hr + Dhr * Xi(j) 'entalpia de reaccion
    Next j


    T = T0 + 0.5 * T0
    R = 0.008314 'kJ/mol/K

    Error = 100
    Indice = 0

    'Inicio Secante

    Dim fa, fb, fc, Tc, Ti, Tf, xtol, ytol As Double
    'a=Ti
    'b=Tf
    'c=Tc
    Ti = T0
    Tf = 4000 'K
    xtol = 0.0001
    ytol = 0.0001
    'fu(Te, T0, T, A, B, Ns, Ne, Hr, R)
    fa = fu(Te, T0, Ti, A, B, Ns, Ne, Hr, R)
    fb = fu(Te, T0, Tf, A, B, Ns, Ne, Hr, R)
    Tc = Tf - fb * (Tf - Ti) / (fb - fa)
    fc = fu(Te, T0, Tc, A, B, Ns, Ne, Hr, R)
    Indice = 0
    Do While Abs(Tf - Tc) >= xtol And Abs(fc) >= ytol And Indice <= 1000000
        Ti = Tf
        fa = fb
        Tf = Tc
        fb = fc
        Tc = Tf - fb * (Tf - Ti) / (fb - fa)
        fc = fu(Te, T0, Tc, A, B, Ns, Ne, Hr, R)
        Indice = Indice + 1
    Loop
    'Fin Secante

    ' Fin Balance de Energía


    Dim Propiedad() As Double
    Dim nn As Integer

    nn = n + 3 ' n = numero de y(i)
    ' 3: Tflama, Qflam, Errflama

    ReDim Propiedad(nn)



    tau = Tc / T0
    tem = B(1) _
        + B(2) / 2 * T0 * (tau + 1) _
        + B(3) / 3 * T0 * T0 * (tau * tau + tau + 1) _
        + B(4) / 4 * T0 * T0 * T0 * (tau * tau * tau + tau * tau + tau + 1) _
        + B(5) / 5 * T0 * T0 * T0 * T0 * (tau * tau * tau * tau + tau * tau * tau + tau * tau + tau + 1)
    Mcps = R * tem
    Qmax = Ns * (Tc - T0) * Mcps

    Propiedad(1) = Tc
    Propiedad(2) = fc
    Propiedad(3) = Qmax

    For i = 1 To n
     Propiedad(3 + i) = ys(i)
    Next i

    
    
    Quemador = Propiedad(ID)

End Function