! Programa elaborado por Nalleli Iridian Avila García (1941603)
! Materia Física Computacional
! Tarea #09

program main
    implicit none
    ! Definición de variables
    double precision :: T1, T2, T3, cp1, cp2, cp3
    double precision :: T_usuario, capacidad_calorifica_interpolada
    double precision :: T_inicial, T_final, delta_H, Q, n
    double precision :: max_third_derivative, error_absoluto

    ! Datos de temperatura y capacidad calorífica molar a presión constante
    T1 = 300.0d0; T2 = 500.0d0; T3 = 2100.0d0
    cp1 = 6.81495d0; cp2 = 8.25375d0; cp3 = 8.6055d0

    ! Permitir al usuario ingresar un nuevo valor de temperatura (T)
    print*, "Ingrese una nueva temperatura (T) para calcular la capacidad calorífica interpolada: "
    read*, T_usuario

    ! Calcular la capacidad calorífica interpolada para la temperatura ingresada por el usuario
    capacidad_calorifica_interpolada = lagrange_polynomial(T_usuario, T1, T2, T3, cp1, cp2, cp3)
    print*, "La capacidad calorífica interpolada a ", T_usuario, " K es: ", capacidad_calorifica_interpolada

    ! Datos del proceso
    T_inicial = 350.0d0  ! Temperatura inicial en Kelvin
    T_final = 800.0d0  ! Temperatura final en Kelvin
    n = 23.0d0  ! Cantidad de sustancia en moles

    ! Calcular el cambio en la entalpía mediante la integral
    delta_H = quad(lagrange_polynomial(T_inicial), lagrange_polynomial(T_final))

    ! Calor absorbido es igual al cambio en la entalpía
    Q = delta_H * n
    print*, "El calor absorbido por los ", n, " moles de nitrógeno en el proceso isobarico cuasiestático es: ", Q

    ! Estimación máxima del valor absoluto de la tercera derivada
    max_third_derivative = 1e-13

    ! Calcular el error absoluto
    error_absoluto = n * (T_final - T_inicial)**4 * max_third_derivative / 24.0d0
    print*, "Estimación del error absoluto: ", error_absoluto
end program main

! Función para calcular el polinomio de segundo grado de Lagrange
double precision function lagrange_polynomial(T)
   implicit none
   double precision :: term1, term2, term3
   
   term1 = ((T - T2) * (T - T3)) / ((T1 - T2) * (T1 - T3)) * cp1
   term2 = ((T - T1) * (T - T3)) / ((T2 - T1) * (T2 - T3)) * cp2
   term3 = ((T - T1) * (T - T2)) / ((T3 - T1) * (T3 - T2)) * cp3
   
   lagrange_polynomial = term1 + term2 + term3
end function lagrange_polynomial

! Función para calcular la integral definida usando el método del trapecio.
double precision function quad(fa, fb)
   implicit none
   double precision :: fa, fb
   
   quad = 0.5d0 * (fb - fa) * (fa + fb)
end function quad