program heat_generation
    implicit none

    integer :: geom_type, n_points, i
    real(8) :: dim_char, k_mat, q_gen, h_conv, temp_inf
    real(8) :: cyl_length
    real(8) :: volume, area_surf, temp_surface, temp_max, delta_temp
    real(8) :: q_total, x_pos, temp_x, ratio_x
    real(8) :: pi
    character(len=20) :: geom_name

    pi = 3.14159265358979d0

    ! Read inputs
    read *, geom_type    ! 1=Plane Wall, 2=Long Cylinder, 3=Sphere
    read *, dim_char     ! Half-thickness L (wall) or radius r0 [m]

    select case (geom_type)
        case (1)
            geom_name = "PLANE WALL"
        case (2)
            geom_name = "LONG CYLINDER"
            read *, cyl_length   ! Cylinder length [m]
        case (3)
            geom_name = "SPHERE"
        case default
            print *, 'ERROR: Invalid geometry type'
            stop
    end select

    read *, k_mat        ! Thermal conductivity [W/m.K]
    read *, q_gen        ! Volumetric heat generation [W/m3]
    read *, h_conv       ! Convection coefficient [W/m2.K]
    read *, temp_inf     ! Ambient temperature [deg-C]

    ! =============================================
    ! GEOMETRY CALCULATIONS
    ! =============================================
    select case (geom_type)
        case (1) ! Plane Wall (per unit area, symmetric)
            volume = 2.0d0 * dim_char * 1.0d0  ! per m2 of wall face
            area_surf = 2.0d0 * 1.0d0           ! both faces per m2
        case (2) ! Long Cylinder
            volume = pi * dim_char**2 * cyl_length
            area_surf = 2.0d0 * pi * dim_char * cyl_length
        case (3) ! Sphere
            volume = (4.0d0 / 3.0d0) * pi * dim_char**3
            area_surf = 4.0d0 * pi * dim_char**2
    end select

    ! =============================================
    ! TEMPERATURE CALCULATIONS
    ! =============================================
    select case (geom_type)
        case (1) ! Plane Wall
            temp_surface = temp_inf + q_gen * dim_char / h_conv
            temp_max = temp_surface + q_gen * dim_char**2 / (2.0d0 * k_mat)
        case (2) ! Long Cylinder
            temp_surface = temp_inf + q_gen * dim_char / (2.0d0 * h_conv)
            temp_max = temp_surface + q_gen * dim_char**2 / (4.0d0 * k_mat)
        case (3) ! Sphere
            temp_surface = temp_inf + q_gen * dim_char / (3.0d0 * h_conv)
            temp_max = temp_surface + q_gen * dim_char**2 / (6.0d0 * k_mat)
    end select

    delta_temp = temp_max - temp_surface
    q_total = q_gen * volume

    ! =============================================
    ! DISPLAY RESULTS
    ! =============================================
    print *, '================================================='
    print *, '   CONDUCTION WITH INTERNAL HEAT GENERATION'
    print *, '================================================='
    print *, ''
    print *, '   Geometry: ', trim(geom_name)
    print *, ''

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

    select case (geom_type)
        case (1)
            print '(A,F12.6,A)', '   Half-thickness (L):        ', dim_char, ' m'
            print '(A,F12.4,A)', '                              ', dim_char*1000.0d0, ' mm'
        case (2)
            print '(A,F12.6,A)', '   Radius (r0):               ', dim_char, ' m'
            print '(A,F12.4,A)', '                              ', dim_char*1000.0d0, ' mm'
            print '(A,F12.4,A)', '   Cylinder Length:            ', cyl_length, ' m'
        case (3)
            print '(A,F12.6,A)', '   Radius (r0):               ', dim_char, ' m'
            print '(A,F12.4,A)', '                              ', dim_char*1000.0d0, ' mm'
    end select

    print '(A,F12.4,A)',  '   Thermal Conductivity (k):  ', k_mat, ' W/m.K'
    print '(A,ES14.4,A)', '   Heat Generation (q_gen):   ', q_gen, ' W/m3'
    print '(A,F12.2,A)',  '   Convection Coeff. (h):     ', h_conv, ' W/m2.K'
    print '(A,F12.2,A)',  '   Ambient Temp. (T_inf):     ', temp_inf, ' deg-C'
    print *, ''

    print *, '================================================='
    print *, '         GEOMETRIC PROPERTIES'
    print *, '================================================='
    print *, ''
    print '(A,ES14.6,A)', '   Volume:                    ', volume, ' m3'
    print '(A,ES14.6,A)', '   Surface Area:              ', area_surf, ' m2'
    print *, ''

    ! =============================================
    ! TEMPERATURE RESULTS
    ! =============================================
    print *, '================================================='
    print *, '       TEMPERATURE DISTRIBUTION'
    print *, '================================================='
    print *, ''

    select case (geom_type)
        case (1)
            print *, '   Formula (Plane Wall, symmetric):'
            print *, '   T(x) = Ts + q_gen/(2k) * (L^2 - x^2)'
            print *, '   Ts = T_inf + q_gen*L/h'
        case (2)
            print *, '   Formula (Long Cylinder):'
            print *, '   T(r) = Ts + q_gen/(4k) * (r0^2 - r^2)'
            print *, '   Ts = T_inf + q_gen*r0/(2h)'
        case (3)
            print *, '   Formula (Sphere):'
            print *, '   T(r) = Ts + q_gen/(6k) * (r0^2 - r^2)'
            print *, '   Ts = T_inf + q_gen*r0/(3h)'
    end select

    print *, ''
    print '(A,F12.2,A)', '   Surface Temperature (Ts):  ', temp_surface, ' deg-C'
    print '(A,F12.2,A)', '   Maximum Temperature (Tmax):', temp_max, ' deg-C'
    print '(A,F12.2,A)', '   Temperature Drop (dT):     ', delta_temp, ' deg-C'
    print '(A,F12.2,A)', '   (Tmax - Ts)'
    print *, ''

    ! =============================================
    ! HEAT TRANSFER
    ! =============================================
    print *, '================================================='
    print *, '         HEAT TRANSFER ANALYSIS'
    print *, '================================================='
    print *, ''
    print '(A,ES14.6,A)', '   Total Heat Generated:      ', q_total, ' W'

    if (q_total >= 1000.0d0) then
        print '(A,F14.4,A)', '                              ', q_total / 1000.0d0, ' kW'
    end if

    print '(A,ES14.6,A)', '   Heat Flux at Surface:      ', q_total / area_surf, ' W/m2'
    print *, ''

    ! Energy balance verification
    print *, '   Energy Balance Verification:'
    print '(A,ES14.6,A)', '   Q_generated = q_gen * V =  ', q_total, ' W'
    print '(A,ES14.6,A)', '   Q_convection = h*As*dTs =  ', &
        h_conv * area_surf * (temp_surface - temp_inf), ' W'
    print *, ''

    ! =============================================
    ! TEMPERATURE PROFILE
    ! =============================================
    print *, '================================================='
    print *, '       TEMPERATURE PROFILE'
    print *, '================================================='
    print *, ''

    select case (geom_type)
        case (1)
            print *, '   x/L       |  x (mm)      |  T (deg-C)   |  (T-Ts)/dT'
        case default
            print *, '   r/r0      |  r (mm)      |  T (deg-C)   |  (T-Ts)/dT'
    end select
    print *, '   --------------------------------------------------------'

    n_points = 25
    do i = 0, n_points
        ratio_x = dble(i) / dble(n_points)
        x_pos = ratio_x * dim_char

        select case (geom_type)
            case (1)
                temp_x = temp_surface + q_gen / (2.0d0 * k_mat) * &
                         (dim_char**2 - x_pos**2)
            case (2)
                temp_x = temp_surface + q_gen / (4.0d0 * k_mat) * &
                         (dim_char**2 - x_pos**2)
            case (3)
                temp_x = temp_surface + q_gen / (6.0d0 * k_mat) * &
                         (dim_char**2 - x_pos**2)
        end select

        write(*, '(3X,F8.4,5X,A,2X,F10.4,5X,A,2X,F10.2,5X,A,2X,F8.5)') &
            ratio_x, '|', x_pos*1000.0d0, '|', temp_x, '|', &
            (temp_x - temp_surface) / max(delta_temp, 1.0d-10)
    end do
    print *, ''

    ! =============================================
    ! THERMAL RESISTANCE ANALYSIS
    ! =============================================
    print *, '================================================='
    print *, '     THERMAL RESISTANCE BREAKDOWN'
    print *, '================================================='
    print *, ''

    select case (geom_type)
        case (1)
            print '(A,ES14.6,A)', '   R_conduction (L/k):        ', &
                dim_char / k_mat, ' m2.K/W'
            print '(A,ES14.6,A)', '   R_convection (1/h):        ', &
                1.0d0 / h_conv, ' m2.K/W'
            print '(A,F12.4)',    '   Ratio R_cond/R_conv:       ', &
                (dim_char / k_mat) / (1.0d0 / h_conv)
        case (2)
            print '(A,ES14.6,A)', '   R_convection:              ', &
                1.0d0 / (h_conv * area_surf), ' K/W'
        case (3)
            print '(A,ES14.6,A)', '   R_convection:              ', &
                1.0d0 / (h_conv * area_surf), ' K/W'
    end select
    print *, ''

    ! =============================================
    ! SAFETY ANALYSIS
    ! =============================================
    print *, '================================================='
    print *, '         SAFETY ANALYSIS'
    print *, '================================================='
    print *, ''

    print '(A,F12.2,A)', '   Maximum Temperature:       ', temp_max, ' deg-C'
    print *, ''

    if (temp_max > 1000.0d0) then
        print *, '   CRITICAL: Tmax > 1000 deg-C'
        print *, '   -------------------------------------------'
        print *, '   Temperature exceeds most material limits!'
        print *, '   Risk of melting or structural failure.'
    else if (temp_max > 500.0d0) then
        print *, '   WARNING: Tmax > 500 deg-C'
        print *, '   -------------------------------------------'
        print *, '   High temperature may cause material'
        print *, '   degradation. Verify material compatibility.'
    else if (temp_max > 200.0d0) then
        print *, '   CAUTION: Tmax > 200 deg-C'
        print *, '   -------------------------------------------'
        print *, '   Moderate temperature. Check thermal limits'
        print *, '   for polymers and organic materials.'
    else
        print *, '   STATUS: Tmax < 200 deg-C'
        print *, '   -------------------------------------------'
        print *, '   Temperature is within safe limits for'
        print *, '   most engineering materials.'
    end if
    print *, ''

    ! =============================================
    ! PARAMETRIC STUDY
    ! =============================================
    print *, '================================================='
    print *, '       PARAMETRIC SENSITIVITY'
    print *, '================================================='
    print *, ''
    print *, '   Effect of doubling key parameters on Tmax:'
    print *, ''

    ! Double q_gen
    select case (geom_type)
        case (1)
            print '(A,F12.2,A)', '   2x q_gen -> Tmax =         ', &
                temp_inf + 2.0d0*q_gen*dim_char/h_conv + &
                2.0d0*q_gen*dim_char**2/(2.0d0*k_mat), ' deg-C'
        case (2)
            print '(A,F12.2,A)', '   2x q_gen -> Tmax =         ', &
                temp_inf + 2.0d0*q_gen*dim_char/(2.0d0*h_conv) + &
                2.0d0*q_gen*dim_char**2/(4.0d0*k_mat), ' deg-C'
        case (3)
            print '(A,F12.2,A)', '   2x q_gen -> Tmax =         ', &
                temp_inf + 2.0d0*q_gen*dim_char/(3.0d0*h_conv) + &
                2.0d0*q_gen*dim_char**2/(6.0d0*k_mat), ' deg-C'
    end select

    ! Double h
    select case (geom_type)
        case (1)
            print '(A,F12.2,A)', '   2x h     -> Tmax =         ', &
                temp_inf + q_gen*dim_char/(2.0d0*h_conv) + &
                q_gen*dim_char**2/(2.0d0*k_mat), ' deg-C'
        case (2)
            print '(A,F12.2,A)', '   2x h     -> Tmax =         ', &
                temp_inf + q_gen*dim_char/(4.0d0*h_conv) + &
                q_gen*dim_char**2/(4.0d0*k_mat), ' deg-C'
        case (3)
            print '(A,F12.2,A)', '   2x h     -> Tmax =         ', &
                temp_inf + q_gen*dim_char/(6.0d0*h_conv) + &
                q_gen*dim_char**2/(6.0d0*k_mat), ' deg-C'
    end select

    ! Double k
    select case (geom_type)
        case (1)
            print '(A,F12.2,A)', '   2x k     -> Tmax =         ', &
                temp_inf + q_gen*dim_char/h_conv + &
                q_gen*dim_char**2/(4.0d0*k_mat), ' deg-C'
        case (2)
            print '(A,F12.2,A)', '   2x k     -> Tmax =         ', &
                temp_inf + q_gen*dim_char/(2.0d0*h_conv) + &
                q_gen*dim_char**2/(8.0d0*k_mat), ' deg-C'
        case (3)
            print '(A,F12.2,A)', '   2x k     -> Tmax =         ', &
                temp_inf + q_gen*dim_char/(3.0d0*h_conv) + &
                q_gen*dim_char**2/(12.0d0*k_mat), ' deg-C'
    end select

    print *, ''
    print '(A,F12.2,A)', '   Current Tmax =             ', temp_max, ' deg-C'
    print *, ''

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

    print *, '   To REDUCE maximum temperature:'
    print *, '   - Decrease object size (reduce L or r0)'
    print *, '   - Increase thermal conductivity (k)'
    print *, '   - Increase convection coefficient (h)'
    print *, '   - Reduce heat generation rate (q_gen)'
    print *, ''

    if (delta_temp > (temp_surface - temp_inf)) then
        print *, '   NOTE: Internal resistance dominates'
        print *, '   -> Focus on increasing k or reducing size'
    else
        print *, '   NOTE: External resistance dominates'
        print *, '   -> Focus on increasing h (better cooling)'
    end if
    print *, ''

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

end program heat_generation