program dependent_k
    implicit none

    ! Constants
    real, parameter :: PI = 3.141592653589793

    ! Inputs
    integer :: geom         ! 1 = Plane, 2 = Cylinder, 3 = Sphere
    real :: L_thick, Area   ! Plane wall parameters (thickness, area)
    real :: r1, r2, H_len   ! Cylinder/Sphere parameters (radii, length)
    integer :: model        ! 1 = Linear, 2 = Polynomial, 3 = Tabulated
    real :: k0, beta        ! Linear parameters
    real :: poly_a, poly_b, poly_c ! Polynomial parameters: k = a + b*T + c*T^2
    integer :: num_pts      ! Number of tabulated points
    real, allocatable :: tab_T(:), tab_k(:), tab_theta(:) ! Tabulated arrays
    real :: T1, T2          ! Boundary temperatures (C)

    ! Outputs & Solved variables
    real :: theta1, theta2
    real :: Q, k_eff
    real :: T_node(11), x_node(11), theta_node(11)
    real :: Q_const, k_const, T_const(11)

    ! Temporary/Iterative variables
    integer :: i, j, temp_idx
    real :: temp_val, k_val, m_val, dT
    real :: r_val, pos_ratio
    real :: f_val, df_val, T_guess

    ! 1. Read Inputs from stdin
    read(*,*) geom
    if (geom == 1) then
        read(*,*) L_thick
        read(*,*) Area
    else if (geom == 2) then
        read(*,*) r1
        read(*,*) r2
        read(*,*) H_len
    else if (geom == 3) then
        read(*,*) r1
        read(*,*) r2
    end if

    read(*,*) model
    if (model == 1) then
        read(*,*) k0
        read(*,*) beta
    else if (model == 2) then
        read(*,*) poly_a
        read(*,*) poly_b
        read(*,*) poly_c
    else if (model == 3) then
        read(*,*) num_pts
        allocate(tab_T(num_pts))
        allocate(tab_k(num_pts))
        allocate(tab_theta(num_pts))
        do i = 1, num_pts
            read(*,*) tab_T(i)
            read(*,*) tab_k(i)
        end do
        
        ! Sort tabulated points in increasing order of T (simple bubble sort)
        do i = 1, num_pts - 1
            do j = i + 1, num_pts
                if (tab_T(i) > tab_T(j)) then
                    temp_val = tab_T(i)
                    tab_T(i) = tab_T(j)
                    tab_T(j) = temp_val

                    temp_val = tab_k(i)
                    tab_k(i) = tab_k(j)
                    tab_k(j) = temp_val
                end if
            end do
        end do
        
        ! Compute cumulative integrals (theta) at tabulated temperatures
        tab_theta(1) = 0.0
        do i = 2, num_pts
            tab_theta(i) = tab_theta(i-1) + 0.5 * (tab_k(i-1) + tab_k(i)) * (tab_T(i) - tab_T(i-1))
        end do
    end if

    read(*,*) T1
    read(*,*) T2

    ! 2. Compute Surface Thetas (theta1 and theta2)
    theta1 = get_theta(T1)
    theta2 = get_theta(T2)

    ! 3. Compute Heat Transfer Rate (Q) & Mean Effective Conductivity (k_eff)
    if (geom == 1) then
        Q = (Area / L_thick) * (theta1 - theta2)
        k_eff = (theta1 - theta2) / (T1 - T2)
    else if (geom == 2) then
        Q = (2.0 * PI * H_len / log(r2 / r1)) * (theta1 - theta2)
        k_eff = (theta1 - theta2) / (T1 - T2)
    else if (geom == 3) then
        Q = (4.0 * PI * r1 * r2 / (r2 - r1)) * (theta1 - theta2)
        k_eff = (theta1 - theta2) / (T1 - T2)
    end if

    ! 4. Generate Node Profiles (11 nodes from inner to outer boundary)
    do i = 1, 11
        pos_ratio = real(i - 1) / 10.0
        
        if (geom == 1) then
            x_node(i) = pos_ratio * L_thick
            theta_node(i) = theta1 - pos_ratio * (theta1 - theta2)
        else
            r_val = r1 + pos_ratio * (r2 - r1)
            x_node(i) = r_val
            if (geom == 2) then
                theta_node(i) = theta1 - (theta1 - theta2) * (log(r_val / r1) / log(r2 / r1))
            else
                theta_node(i) = theta1 - (theta1 - theta2) * ((1.0/r1 - 1.0/r_val) / (1.0/r1 - 1.0/r2))
            end if
        end if

        T_node(i) = get_temp_from_theta(theta_node(i))
    end do

    ! 5. Comparison with Constant-k Solver (evaluating at average temperature)
    k_const = get_k_at_temp(0.5 * (T1 + T2))
    if (geom == 1) then
        Q_const = (k_const * Area / L_thick) * (T1 - T2)
        do i = 1, 11
            T_const(i) = T1 - (real(i - 1) / 10.0) * (T1 - T2)
        end do
    else if (geom == 2) then
        Q_const = (2.0 * PI * k_const * H_len / log(r2 / r1)) * (T1 - T2)
        do i = 1, 11
            r_val = r1 + (real(i - 1) / 10.0) * (r2 - r1)
            T_const(i) = T1 - (T1 - T2) * (log(r_val / r1) / log(r2 / r1))
        end do
    else if (geom == 3) then
        Q_const = (4.0 * PI * k_const * r1 * r2 / (r2 - r1)) * (T1 - T2)
        do i = 1, 11
            r_val = r1 + (real(i - 1) / 10.0) * (r2 - r1)
            T_const(i) = T1 - (T1 - T2) * ((1.0/r1 - 1.0/r_val) / (1.0/r1 - 1.0/r2))
        end do
    end if

    ! 6. Output Engineering Report
    write(*,*) "=========================================================================="
    write(*,*) "         STEADY CONDUCTION WITH VARIABLE PROPERTIES k = f(T)              "
    write(*,*) "=========================================================================="
    
    if (geom == 1) then
        write(*, '(A)') " Geometry:               PLANE WALL"
        write(*, '(A, F10.4, A, F10.4, A)') " Dimensions:             Thickness = ", L_thick, " m, Area = ", Area, " m2"
    else if (geom == 2) then
        write(*, '(A)') " Geometry:               CYLINDER"
        write(*, '(A, F10.4, A, F10.4, A, F10.4, A)') &
            " Dimensions:             r1 = ", r1, " m, r2 = ", r2, " m, Length = ", H_len, " m"
    else if (geom == 3) then
        write(*, '(A)') " Geometry:               SPHERE"
        write(*, '(A, F10.4, A, F10.4, A)') " Dimensions:             r1 = ", r1, " m, r2 = ", r2, " m"
    end if

    if (model == 1) then
        write(*, '(A)') " Property Model:         Linear k(T) = k0*(1 + beta*T)"
        write(*, '(A, F10.4, A, E13.4, A)') " Parameters:             k0 = ", k0, " W/(m-K), beta = ", beta, " 1/K"
    else if (model == 2) then
        write(*, '(A)') " Property Model:         Polynomial k(T) = a + b*T + c*T^2"
        write(*, '(A, F10.4, A, E13.4, A, E13.4, A)') &
            " Parameters:             a = ", poly_a, ", b = ", poly_b, ", c = ", poly_c
    else if (model == 3) then
        write(*, '(A)') " Property Model:         Tabulated / Interpolated points"
        write(*, '(A, I3)') " Tabulated points:       ", num_pts
        do i = 1, num_pts
            write(*, '(A, I2, A, F8.2, A, F8.3)') "   Point ", i, ": T = ", tab_T(i), " C, k = ", tab_k(i)
        end do
    end if

    write(*, '(A, F8.2, A, F8.2, A)') " Boundary Temperatures:  T1 = ", T1, " C, T2 = ", T2, " C"
    write(*,*) "--------------------------------------------------------------------------"
    write(*,*) " THERMAL PERFORMANCE ANALYSIS"
    write(*,*) "--------------------------------------------------------------------------"
    write(*, '(A, F14.2, A)') " Heat Transfer Rate (Kirchhoff Q):     ", Q, " W"
    write(*, '(A, F14.2, A)') " Heat Transfer Rate (Constant k_avg Q):", Q_const, " W"
    write(*, '(A, F14.2, A)') " Absolute Difference in Q:             ", abs(Q - Q_const), " W"
    write(*, '(A, F14.4, A)') " Mean Effective Conductivity (k_eff):  ", k_eff, " W/(m-K)"
    write(*, '(A, F14.4, A)') " Constant Conductivity (k_avg):        ", k_const, " W/(m-K)"

    write(*,*) "--------------------------------------------------------------------------"
    write(*,*) " SPATIAL PROFILE DATA"
    write(*,*) "--------------------------------------------------------------------------"
    if (geom == 1) then
        write(*,*) " Node   Position x(m)   Transformed theta   Temp T(C)   Temp Constant-k(C)"
    else
        write(*,*) " Node   Radius r(m)     Transformed theta   Temp T(C)   Temp Constant-k(C)"
    end if
    
    do i = 1, 11
        write(*, '(I4.2, F16.4, F20.4, F12.2, F21.2)') &
            i, x_node(i), theta_node(i), T_node(i), T_const(i)
    end do
    write(*,*) "=========================================================================="

contains

    ! Helper: evaluates theta(T) = \int_0^T k(T') dT'
    real function get_theta(temp)
        real, intent(in) :: temp
        integer :: k_idx
        real :: m_local, dT_local

        if (model == 1) then
            get_theta = k0 * (temp + 0.5 * beta * temp**2)
        else if (model == 2) then
            get_theta = poly_a * temp + 0.5 * poly_b * temp**2 + (poly_c / 3.0) * temp**3
        else if (model == 3) then
            ! For tabulated model, T_ref = tab_T(1)
            if (temp <= tab_T(1)) then
                ! Extrapolate using first interval gradient
                m_local = (tab_k(2) - tab_k(1)) / (tab_T(2) - tab_T(1))
                dT_local = temp - tab_T(1)
                get_theta = tab_k(1) * dT_local + 0.5 * m_local * dT_local**2
            else if (temp >= tab_T(num_pts)) then
                ! Extrapolate using last interval gradient
                m_local = (tab_k(num_pts) - tab_k(num_pts-1)) / (tab_T(num_pts) - tab_T(num_pts-1))
                dT_local = temp - tab_T(num_pts)
                get_theta = tab_theta(num_pts) + tab_k(num_pts) * dT_local + 0.5 * m_local * dT_local**2
            else
                ! Find interval
                k_idx = 1
                do while (tab_T(k_idx+1) < temp .and. k_idx < num_pts - 1)
                    k_idx = k_idx + 1
                end do
                m_local = (tab_k(k_idx+1) - tab_k(k_idx)) / (tab_T(k_idx+1) - tab_T(k_idx))
                dT_local = temp - tab_T(k_idx)
                get_theta = tab_theta(k_idx) + tab_k(k_idx) * dT_local + 0.5 * m_local * dT_local**2
            end if
        end if
    end function get_theta

    ! Helper: evaluates k(T) at a given temperature T
    real function get_k_at_temp(temp)
        real, intent(in) :: temp
        integer :: k_idx
        real :: m_local

        if (model == 1) then
            get_k_at_temp = k0 * (1.0 + beta * temp)
        else if (model == 2) then
            get_k_at_temp = poly_a + poly_b * temp + poly_c * temp**2
        else if (model == 3) then
            if (temp <= tab_T(1)) then
                m_local = (tab_k(2) - tab_k(1)) / (tab_T(2) - tab_T(1))
                get_k_at_temp = tab_k(1) + m_local * (temp - tab_T(1))
            else if (temp >= tab_T(num_pts)) then
                m_local = (tab_k(num_pts) - tab_k(num_pts-1)) / (tab_T(num_pts) - tab_T(num_pts-1))
                get_k_at_temp = tab_k(num_pts) + m_local * (temp - tab_T(num_pts))
            else
                k_idx = 1
                do while (tab_T(k_idx+1) < temp .and. k_idx < num_pts - 1)
                    k_idx = k_idx + 1
                end do
                m_local = (tab_k(k_idx+1) - tab_k(k_idx)) / (tab_T(k_idx+1) - tab_T(k_idx))
                get_k_at_temp = tab_k(k_idx) + m_local * (temp - tab_T(k_idx))
            end if
        end if
    end function get_k_at_temp

    ! Helper: inverts theta to retrieve T
    real function get_temp_from_theta(theta_val)
        real, intent(in) :: theta_val
        real :: disc, root
        integer :: iter, k_idx
        real :: m_local, T_val, f_val, df_local

        if (model == 1) then
            if (beta == 0.0) then
                get_temp_from_theta = theta_val / k0
            else
                disc = 1.0 + 2.0 * beta * theta_val / k0
                if (disc >= 0.0) then
                    ! Standard physical root
                    root = sqrt(disc)
                    if (1.0 + beta * 0.5 * (T1 + T2) >= 0.0) then
                        get_temp_from_theta = (-1.0 + root) / beta
                    else
                        get_temp_from_theta = (-1.0 - root) / beta
                    end if
                else
                    get_temp_from_theta = theta_val / k0 ! fallback
                end if
            end if
        else if (model == 2) then
            ! Solve a*T + (b/2)*T^2 + (c/3)*T^3 - theta = 0 using Newton's method
            T_val = 0.5 * (T1 + T2) ! initial guess
            iter = 0
            do while (iter < 100)
                iter = iter + 1
                f_val = poly_a * T_val + 0.5 * poly_b * T_val**2 + (poly_c / 3.0) * T_val**3 - theta_val
                df_local = poly_a + poly_b * T_val + poly_c * T_val**2
                if (abs(df_local) < 1e-6) exit
                dT = f_val / df_local
                T_val = T_val - dT
                if (abs(dT) < 1e-5) exit
            end do
            get_temp_from_theta = T_val
        else if (model == 3) then
            ! Tabulated inversion
            if (theta_val <= 0.0) then
                ! Extrapolate below tab_T(1)
                m_local = (tab_k(2) - tab_k(1)) / (tab_T(2) - tab_T(1))
                if (m_local == 0.0) then
                    get_temp_from_theta = tab_T(1) + theta_val / tab_k(1)
                else
                    disc = tab_k(1)**2 + 2.0 * m_local * theta_val
                    if (disc >= 0.0) then
                        get_temp_from_theta = tab_T(1) + (-tab_k(1) + sqrt(disc)) / m_local
                    else
                        get_temp_from_theta = tab_T(1) + theta_val / tab_k(1)
                    end if
                end if
            else if (theta_val >= tab_theta(num_pts)) then
                ! Extrapolate above tab_T(num_pts)
                m_local = (tab_k(num_pts) - tab_k(num_pts-1)) / (tab_T(num_pts) - tab_T(num_pts-1))
                if (m_local == 0.0) then
                    get_temp_from_theta = tab_T(num_pts) + (theta_val - tab_theta(num_pts)) / tab_k(num_pts)
                else
                    disc = tab_k(num_pts)**2 + 2.0 * m_local * (theta_val - tab_theta(num_pts))
                    if (disc >= 0.0) then
                        get_temp_from_theta = tab_T(num_pts) + (-tab_k(num_pts) + sqrt(disc)) / m_local
                    else
                        get_temp_from_theta = tab_T(num_pts) + (theta_val - tab_theta(num_pts)) / tab_k(num_pts)
                    end if
                end if
            else
                ! Find interval
                k_idx = 1
                do while (tab_theta(k_idx+1) < theta_val .and. k_idx < num_pts - 1)
                    k_idx = k_idx + 1
                end do
                m_local = (tab_k(k_idx+1) - tab_k(k_idx)) / (tab_T(k_idx+1) - tab_T(k_idx))
                if (m_local == 0.0) then
                    get_temp_from_theta = tab_T(k_idx) + (theta_val - tab_theta(k_idx)) / tab_k(k_idx)
                else
                    disc = tab_k(k_idx)**2 + 2.0 * m_local * (theta_val - tab_theta(k_idx))
                    if (disc >= 0.0) then
                        get_temp_from_theta = tab_T(k_idx) + (-tab_k(k_idx) + sqrt(disc)) / m_local
                    else
                        get_temp_from_theta = tab_T(k_idx) + (theta_val - tab_theta(k_idx)) / tab_k(k_idx)
                    end if
                end if
            end if
        end if
    end function get_temp_from_theta

end program dependent_k
