! Tarea#11 
! Fisica computacional (2:00-3:00 pm)
! Programa elaborado por Nalleli Iridian Avila Garcia (1941603)

program main
    implicit none
    double precision, dimension(3) :: T = [300,500,2100]
    double precision, dimension(3) :: Cp = [6.81495,8.25375,8.60055]
    double precision :: a, b, c, Ti, Tf, Q, error
    integer :: i

    print*, 'Este programa calcula los coeficientes de un polinomio de segundo grado que se ajusta a los datos de capacidad calorífica molar en función de la temperatura. Luego utiliza este polinomio para calcular el calor total absorbido por 23 moles de nitrógeno en un proceso isobárico cuasiestático. Finalmente calcula una cota del error absoluto cometido en el inciso anterior.'

    ! Aproximación de la función Cp(T) mediante una expresión polinomial de 2do grado
    call polyfit(T, Cp, a, b, c)

    ! Calor total absorbido por 23 moles de nitrógeno en un proceso isobárico cuasiestático
    Ti = 350
    Tf = 800
    Q = 23 * quad(a, b, c, Ti, Tf)

    ! Cota del error absoluto cometido en el inciso anterior
    error = abs((Tf - Ti) * 10**-13)

    print*, 'a) Coeficientes del polinomio: ', a, b, c
    print*, 'b) Calor total absorbido: ', Q
    print*, 'c) Cota del error absoluto: ', error

end program main

subroutine polyfit(x, y, a, b, c)
    implicit none
    double precision, dimension(:), intent(in) :: x, y
    double precision, intent(out) :: a, b, c
    double precision :: x1 = x(1), x2 = x(2), x3 = x(3)
    double precision :: y1 = y(1), y2 = y(2), y3 = y(3)

    a = ((y3 - y1)*(x1 - x2) - (y2 - y1)*(x1 - x3)) / ((x1**2 - x2**2)*(x1 - x3) - (x1**2 - x3**2)*(x1 - x2))
    b = ((y2 - y1) - a*(x1**2 - x2**2)) / (x1 - x2)
    c = y1 - a*x1**2 - b*x1

end subroutine polyfit

double precision function quad(a, b, c, Ti, Tf)
    implicit none
    double precision, intent(in) :: a, b, c, Ti, Tf

    quad = (a/3)*(Tf**3 - Ti**3) + 0.5*b*(Tf**2 - Ti**2) + c*(Tf - Ti)

end function quad