program lumped_system_analysis
    implicit none

    ! Variable declarations
    integer :: geometry_type, n_points, i
    real(8) :: k, rho, cp, h, T_initial, T_ambient, T_target
    real(8) :: radius, length_dim, side_dim
    real(8) :: volume, A_surface, Lc, Bi, tau, alpha
    real(8) :: T_t, t, dt, t_max, time_target, Q_t, Q_max
    real(8) :: theta_ratio, pi
    character(len=20) :: geometry_name

    pi = 3.14159265358979d0

    ! =============================================
    ! READ INPUT PARAMETERS
    ! =============================================
    read *, geometry_type   ! 1=Sphere, 2=Long Cylinder, 3=Plane Wall, 4=Cube

    select case (geometry_type)
        case (1)
            geometry_name = "SPHERE"
            read *, radius          ! Radius [m]
        case (2)
            geometry_name = "LONG CYLINDER"
            read *, radius          ! Radius [m]
            read *, length_dim      ! Length [m]
        case (3)
            geometry_name = "PLANE WALL"
            read *, length_dim      ! Half-thickness L [m]
            read *, side_dim        ! Area of one face [m2]
        case (4)
            geometry_name = "CUBE"
            read *, side_dim        ! Side length [m]
        case default
            print *, 'ERROR: Invalid geometry type'
            stop
    end select

    read *, k               ! Thermal conductivity [W/m.K]
    read *, rho             ! Density [kg/m3]
    read *, cp              ! Specific heat capacity [J/kg.K]
    read *, h               ! Convection coefficient [W/m2.K]
    read *, T_initial       ! Initial temperature [deg-C]
    read *, T_ambient       ! Ambient temperature [deg-C]
    read *, T_target        ! Target temperature [deg-C] (0 = no target)

    ! =============================================
    ! CALCULATE GEOMETRY
    ! =============================================
    select case (geometry_type)
        case (1) ! Sphere
            volume = (4.0d0 / 3.0d0) * pi * radius**3
            A_surface = 4.0d0 * pi * radius**2
        case (2) ! Long Cylinder
            volume = pi * radius**2 * length_dim
            A_surface = 2.0d0 * pi * radius * length_dim + 2.0d0 * pi * radius**2
        case (3) ! Plane Wall (2L thickness, area A on each face)
            volume = 2.0d0 * length_dim * side_dim
            A_surface = 2.0d0 * side_dim
        case (4) ! Cube
            volume = side_dim**3
            A_surface = 6.0d0 * side_dim**2
    end select

    ! Characteristic length
    Lc = volume / A_surface

    ! Thermal diffusivity
    alpha = k / (rho * cp)

    ! Biot number
    Bi = h * Lc / k

    ! Time constant
    tau = rho * cp * volume / (h * A_surface)

    ! Maximum energy transfer
    Q_max = rho * cp * volume * abs(T_initial - T_ambient)

    ! =============================================
    ! DISPLAY RESULTS
    ! =============================================
    print *, '================================================='
    print *, '   LUMPED SYSTEM ANALYSIS'
    print *, '   (Transient Heat Conduction)'
    print *, '================================================='
    print *, ''
    print *, '   Geometry: ', trim(geometry_name)
    print *, ''

    print *, '================================================='
    print *, '          INPUT PARAMETERS'
    print *, '================================================='
    print *, ''

    select case (geometry_type)
        case (1)
            print '(A,F12.6,A)', '   Radius:                    ', radius, ' m'
            print '(A,F12.4,A)', '                              ', radius*1000.0d0, ' mm'
        case (2)
            print '(A,F12.6,A)', '   Radius:                    ', radius, ' m'
            print '(A,F12.4,A)', '   Length:                    ', length_dim, ' m'
        case (3)
            print '(A,F12.6,A)', '   Half-thickness (L):        ', length_dim, ' m'
            print '(A,F12.4,A)', '                              ', length_dim*1000.0d0, ' mm'
            print '(A,F12.4,A)', '   Face Area (A):             ', side_dim, ' m2'
        case (4)
            print '(A,F12.6,A)', '   Side Length:               ', side_dim, ' m'
            print '(A,F12.4,A)', '                              ', side_dim*1000.0d0, ' mm'
    end select

    print '(A,F12.4,A)', '   Thermal Conductivity (k):  ', k, ' W/m.K'
    print '(A,F12.2,A)', '   Density (rho):             ', rho, ' kg/m3'
    print '(A,F12.2,A)', '   Specific Heat (cp):        ', cp, ' J/kg.K'
    print '(A,F12.2,A)', '   Convection Coeff. (h):     ', h, ' W/m2.K'
    print '(A,F12.2,A)', '   Initial Temperature (Ti):  ', T_initial, ' deg-C'
    print '(A,F12.2,A)', '   Ambient Temperature (T_inf):', T_ambient, ' deg-C'
    if (T_target /= 0.0d0) then
        print '(A,F12.2,A)', '   Target Temperature:        ', T_target, ' deg-C'
    end if
    print *, ''

    print *, '================================================='
    print *, '         GEOMETRIC CALCULATIONS'
    print *, '================================================='
    print *, ''
    print '(A,ES14.6,A)', '   Volume (V):                ', volume, ' m3'
    print '(A,ES14.6,A)', '   Surface Area (As):         ', A_surface, ' m2'
    print '(A,ES14.6,A)', '   Characteristic Length (Lc):', Lc, ' m'
    print '(A,F12.4,A)',  '                              ', Lc*1000.0d0, ' mm'
    print '(A,ES14.6,A)', '   Thermal Diffusivity (a):   ', alpha, ' m2/s'
    print *, ''

    ! =============================================
    ! BIOT NUMBER ANALYSIS
    ! =============================================
    print *, '================================================='
    print *, '         BIOT NUMBER ANALYSIS'
    print *, '================================================='
    print *, ''
    print *, '   Bi = h * Lc / k'
    print '(A,ES14.6)',   '   Bi =                       ', Bi
    print *, ''

    if (Bi < 0.1d0) then
        print *, '   STATUS: Bi < 0.1'
        print *, '   -----------------------------------------------'
        print *, '   The lumped system analysis IS VALID.'
        print *, '   Internal temperature gradients are negligible.'
        print *, '   The object can be treated as spatially uniform.'
    else if (Bi < 1.0d0) then
        print *, '   WARNING: 0.1 < Bi < 1.0'
        print *, '   -----------------------------------------------'
        print *, '   Lumped analysis has LIMITED accuracy.'
        print *, '   Consider using the one-term approximation'
        print *, '   (Heisler charts) for better precision.'
        print *, '   Error may be up to 5%.'
    else
        print *, '   ERROR: Bi > 1.0'
        print *, '   -----------------------------------------------'
        print *, '   Lumped system analysis is NOT VALID!'
        print *, '   Internal temperature gradients are significant.'
        print *, '   Use the Heisler chart method or numerical'
        print *, '   methods for accurate results.'
    end if
    print *, ''

    ! =============================================
    ! TIME CONSTANT
    ! =============================================
    print *, '================================================='
    print *, '         TIME CONSTANT'
    print *, '================================================='
    print *, ''
    print *, '   tau = rho * cp * V / (h * As)'
    print '(A,F14.4,A)', '   tau =                      ', tau, ' s'

    if (tau >= 3600.0d0) then
        print '(A,F14.4,A)', '                              ', tau/3600.0d0, ' hours'
    else if (tau >= 60.0d0) then
        print '(A,F14.4,A)', '                              ', tau/60.0d0, ' min'
    end if
    print *, ''
    print *, '   Physical meaning:'
    print *, '   At t = tau, the object reaches 63.2% of the'
    print *, '   total temperature change.'
    print *, '   At t = 3*tau, it reaches 95.0%.'
    print *, '   At t = 5*tau, it reaches 99.3%.'
    print *, ''

    ! =============================================
    ! TARGET TEMPERATURE
    ! =============================================
    if (T_target /= 0.0d0) then
        theta_ratio = (T_target - T_ambient) / (T_initial - T_ambient)

        if (theta_ratio > 0.0d0 .and. theta_ratio < 1.0d0) then
            time_target = -tau * log(theta_ratio)

            print *, '================================================='
            print *, '      TIME TO REACH TARGET TEMPERATURE'
            print *, '================================================='
            print *, ''
            print '(A,F12.2,A)', '   Target Temperature:        ', T_target, ' deg-C'
            print '(A,F14.4,A)', '   Time Required:             ', time_target, ' s'

            if (time_target >= 3600.0d0) then
                print '(A,F14.4,A)', '                              ', time_target/3600.0d0, ' hours'
            else if (time_target >= 60.0d0) then
                print '(A,F14.4,A)', '                              ', time_target/60.0d0, ' min'
            end if

            ! Energy transferred at target time
            Q_t = Q_max * (1.0d0 - exp(-time_target / tau))
            print *, ''
            print '(A,F14.4,A)', '   Energy Transferred:        ', Q_t, ' J'
            if (Q_t >= 1000.0d0) then
                print '(A,F14.4,A)', '                              ', Q_t/1000.0d0, ' kJ'
            end if
            print *, ''
        else if (theta_ratio <= 0.0d0) then
            print *, '================================================='
            print *, '      TARGET TEMPERATURE ANALYSIS'
            print *, '================================================='
            print *, ''
            print *, '   WARNING: Target temperature is beyond the'
            print *, '   ambient temperature. The object will'
            print *, '   asymptotically approach T_ambient but'
            print *, '   never reach/cross it.'
            print *, ''
        else
            print *, '================================================='
            print *, '      TARGET TEMPERATURE ANALYSIS'
            print *, '================================================='
            print *, ''
            print *, '   NOTE: Target temperature is between Ti'
            print *, '   and T_ambient or equals Ti. Already reached.'
            print *, ''
        end if
    end if

    ! =============================================
    ! ENERGY ANALYSIS
    ! =============================================
    print *, '================================================='
    print *, '         ENERGY ANALYSIS'
    print *, '================================================='
    print *, ''
    print '(A,F14.4,A)', '   Max Energy Transfer (Qmax):', Q_max, ' J'
    if (Q_max >= 1000.0d0) then
        print '(A,F14.4,A)', '                              ', Q_max/1000.0d0, ' kJ'
    end if
    if (Q_max >= 1.0d6) then
        print '(A,F14.4,A)', '                              ', Q_max/1.0d6, ' MJ'
    end if
    print *, ''
    print *, '   Energy at key times:'
    print '(A,F10.2,A,F14.4,A)', '   At t = 1*tau (', tau, ' s): Q = ', Q_max * 0.6321d0, ' J'
    print '(A,F10.2,A,F14.4,A)', '   At t = 2*tau (', 2.0d0*tau, ' s): Q = ', Q_max * 0.8647d0, ' J'
    print '(A,F10.2,A,F14.4,A)', '   At t = 3*tau (', 3.0d0*tau, ' s): Q = ', Q_max * 0.9502d0, ' J'
    print '(A,F10.2,A,F14.4,A)', '   At t = 5*tau (', 5.0d0*tau, ' s): Q = ', Q_max * 0.9933d0, ' J'
    print *, ''

    ! =============================================
    ! TEMPERATURE DISTRIBUTION OVER TIME
    ! =============================================
    print *, '================================================='
    print *, '      TEMPERATURE vs TIME'
    print *, '================================================='
    print *, ''
    print *, '   Time (s)    |  T (deg-C)   |  theta/theta0  |  Q/Qmax'
    print *, '   ---------------------------------------------------------'

    n_points = 50
    t_max = 5.0d0 * tau
    dt = t_max / dble(n_points)

    do i = 0, n_points
        t = dble(i) * dt
        theta_ratio = exp(-t / tau)
        T_t = T_ambient + (T_initial - T_ambient) * theta_ratio
        Q_t = Q_max * (1.0d0 - theta_ratio)

        write(*, '(3X,F12.4,4X,A,2X,F10.3,6X,A,2X,F10.6,4X,A,2X,F8.5)') &
            t, '|', T_t, '|', theta_ratio, '|', Q_t / Q_max
    end do
    print *, ''

    ! =============================================
    ! HEAT TRANSFER RATE
    ! =============================================
    print *, '================================================='
    print *, '      INSTANTANEOUS HEAT TRANSFER RATE'
    print *, '================================================='
    print *, ''
    print *, '   dQ/dt = h * As * (T(t) - T_inf)'
    print *, ''

    ! At t=0
    print '(A,F14.4,A)', '   At t = 0:    dQ/dt =       ', &
        h * A_surface * abs(T_initial - T_ambient), ' W'
    ! At t=tau
    print '(A,F14.4,A)', '   At t = tau:  dQ/dt =       ', &
        h * A_surface * abs(T_initial - T_ambient) * exp(-1.0d0), ' W'
    ! At t=3tau
    print '(A,F14.4,A)', '   At t = 3tau: dQ/dt =       ', &
        h * A_surface * abs(T_initial - T_ambient) * exp(-3.0d0), ' W'
    print *, ''

    ! =============================================
    ! DESIGN RECOMMENDATIONS
    ! =============================================
    print *, '================================================='
    print *, '      DESIGN RECOMMENDATIONS'
    print *, '================================================='
    print *, ''

    if (Bi < 0.1d0) then
        print *, '   VALIDITY: Analysis is valid (Bi < 0.1)'
        print *, ''
    end if

    ! Suggestions to speed up or slow down
    if (T_initial > T_ambient) then
        print *, '   To SPEED UP cooling:'
    else
        print *, '   To SPEED UP heating:'
    end if
    print *, '   - Increase h (forced convection, fluid change)'
    print *, '   - Decrease object size (reduce V/As ratio)'
    print *, '   - Use material with lower rho*cp product'
    print *, ''

    if (T_initial > T_ambient) then
        print *, '   To SLOW DOWN cooling:'
    else
        print *, '   To SLOW DOWN heating:'
    end if
    print *, '   - Add insulation (decrease h)'
    print *, '   - Use larger objects (increase V/As ratio)'
    print *, '   - Use material with higher rho*cp product'
    print *, ''

    ! Convection regime recommendations
    print *, '   Current convection regime:'
    if (h < 10.0d0) then
        print *, '   - Natural convection in air (h < 10)'
        print *, '   - Consider forced air (h = 25-250) for faster results'
    else if (h < 250.0d0) then
        print *, '   - Forced convection in air (10 < h < 250)'
        print *, '   - Consider liquid cooling (h > 500) if needed'
    else if (h < 1000.0d0) then
        print *, '   - Moderate liquid convection (250 < h < 1000)'
    else
        print *, '   - Strong liquid convection or boiling (h > 1000)'
    end if
    print *, ''

    print *, '================================================='
    print *, '          CALCULATION COMPLETE'
    print *, '================================================='
    print *, ''

end program lumped_system_analysis