program boiling_condensation
    implicit none

    ! Inputs
    integer :: mode         ! 1 = Pool Boiling, 2 = Film Condensation
    integer :: geom_type    ! For Cond: 1=Vertical Plate, 2=Cylinder, 3=Sphere. For Boiling: ignored.
    double precision :: dim1, dim2
    double precision :: T_s, T_sat
    double precision :: rho_l, mu_l, k_l, Pr_l, Cp_l
    double precision :: rho_v, mu_v, k_v, Pr_v, Cp_v
    double precision :: h_fg, sigma, C_sf, emiss

    ! Constants
    double precision, parameter :: g = 9.80665d0
    double precision, parameter :: sigma_sb = 5.67d-8
    double precision, parameter :: PI = 3.141592653589793d0

    ! Solved variables
    double precision :: dT_e, dT, h_fg_mod, h_avg, Q, A_s, flux, q_max, q_min, L_c
    double precision :: h_film, h_rad, T_s_K, T_sat_K
    double precision :: mass_flow

    ! Local profile
    integer :: i, n_points
    double precision :: x, dx, local_dT_e, local_flux, log_dT, local_delta
    character(len=50) :: regime_name
    integer :: iostat_val

    ! Read inputs
    read(*,*,iostat=iostat_val) mode
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid mode input.'
        stop
    end if

    read(*,*,iostat=iostat_val) geom_type
    read(*,*,iostat=iostat_val) dim1
    read(*,*,iostat=iostat_val) dim2
    read(*,*,iostat=iostat_val) T_s
    read(*,*,iostat=iostat_val) T_sat
    read(*,*,iostat=iostat_val) rho_l
    read(*,*,iostat=iostat_val) mu_l
    read(*,*,iostat=iostat_val) k_l
    read(*,*,iostat=iostat_val) Pr_l
    read(*,*,iostat=iostat_val) Cp_l
    read(*,*,iostat=iostat_val) rho_v
    read(*,*,iostat=iostat_val) mu_v
    read(*,*,iostat=iostat_val) k_v
    read(*,*,iostat=iostat_val) Pr_v
    read(*,*,iostat=iostat_val) Cp_v
    read(*,*,iostat=iostat_val) h_fg
    read(*,*,iostat=iostat_val) sigma
    read(*,*,iostat=iostat_val) C_sf
    read(*,*,iostat=iostat_val) emiss

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all simulation parameters.'
        stop
    end if

    ! Basic validations
    if (dim1 <= 0.0d0) dim1 = 0.01d0 ! Safe cylinder diameter default
    if (rho_l <= 0.0d0 .or. h_fg <= 0.0d0 .or. sigma <= 0.0d0) then
        write(*,*) 'ERROR: Density, latent heat, and surface tension must be positive.'
        stop
    end if

    if (mode == 1) then
        ! ============================================================
        ! POOL BOILING ANALYSIS
        ! ============================================================
        dT_e = T_s - T_sat
        if (dT_e < 0.0d0) then
            dT_e = 0.0d0
        end if

        ! Zuber CHF calculation
        q_max = 0.149d0 * h_fg * rho_v * ((g * sigma * (rho_l - rho_v)) / (rho_v**2))**0.25d0

        ! Leidenfrost point minimum heat flux
        q_min = 0.09d0 * rho_v * h_fg * ((g * sigma * (rho_l - rho_v)) / ((rho_l + rho_v)**2))**0.25d0

        ! Saturated pool boiling curve piecewise modeling
        T_s_K = T_s + 273.15d0
        T_sat_K = T_sat + 273.15d0

        ! Nucleate boiling Rohsenow correlation
        if (dT_e > 0.0d0) then
            flux = mu_l * h_fg * sqrt((g * (rho_l - rho_v)) / sigma) * &
                   ((Cp_l * dT_e) / (C_sf * h_fg * Pr_l))**3
        else
            flux = 0.0d0
        end if

        ! Check Bromley film boiling (at Leidenfrost or higher)
        h_fg_mod = h_fg + 0.8d0 * Cp_v * dT_e
        h_film = 0.62d0 * ((g * rho_v * (rho_l - rho_v) * k_v**3 * h_fg_mod) / &
                           (mu_v * dim1 * max(dT_e, 1.0d-5)))**0.25d0
        h_rad = (emiss * sigma_sb * (T_s_K**4 - T_sat_K**4)) / max(dT_e, 1.0d-5)
        h_avg = h_film + 0.75d0 * h_rad

        ! Piecewise heat flux selection matching Nukiyama regimes
        if (dT_e < 5.0d0) then
            regime_name = "Natural Convection Boiling"
            flux = 500.0d0 * dT_e**1.25d0
        else if (dT_e < 30.0d0) then
            regime_name = "Nucleate Boiling"
            ! Rohsenow flux matches standard nucleate regime
            if (flux > q_max) flux = q_max
        else if (dT_e < 120.0d0) then
            regime_name = "Transition Boiling"
            ! Log-log interpolation between CHF and Leidenfrost points
            log_dT = (log10(dT_e) - log10(30.0d0)) / (log10(120.0d0) - log10(30.0d0))
            flux = 10.0d0**(log10(q_max) + log_dT * (log10(q_min) - log10(q_max)))
        else
            regime_name = "Film Boiling"
            flux = h_avg * dT_e
            if (flux < q_min) flux = q_min
        end if

        ! Output Boiling results
        write(*,'(A)') '============================================================'
        write(*,'(A)') '   POOL BOILING HEAT TRANSFER CURVE ANALYSIS'
        write(*,'(A)') '============================================================'
        write(*,*)
        write(*,'(A,A)')        '  Boiling Regime          = ', trim(regime_name)
        write(*,'(A,F12.2,A)')  '  Excess Temperature (dTe)= ', dT_e, ' deg-C'
        write(*,'(A,F12.2,A)')  '  Surface Temperature (Ts)= ', T_s, ' deg-C'
        write(*,'(A,F12.2,A)')  '  Saturation Temp (Tsat)  = ', T_sat, ' deg-C'
        write(*,*)
        write(*,'(A)') '--- TRANSITION BOUNDARIES -----------------------------------'
        write(*,'(A,ES12.4,A)') '  Critical Heat Flux (CHF)= ', q_max, ' W/m2'
        write(*,'(A,ES12.4,A)') '  Leidenfrost Flux (min)  = ', q_min, ' W/m2'
        write(*,*)
        write(*,'(A)') '--- HEAT TRANSFER RESULTS -----------------------------------'
        write(*,'(A,ES12.4,A)') '  Calculated Heat Flux (q")= ', flux, ' W/m2'
        write(*,'(A,F12.4,A)')  '  Film HTC (Convective)   = ', h_film, ' W/m2.K'
        write(*,'(A,F12.4,A)')  '  Film HTC (Radiation)    = ', h_rad, ' W/m2.K'
        write(*,'(A,F12.4,A)')  '  Total HTC (h_total)     = ', h_avg, ' W/m2.K'
        write(*,*)

        ! Generate log boiling curve dataset
        write(*,'(A)') '--- LOCAL PROFILE ALONG SURFACE ----------------------------'
        write(*,'(A)') '  dTe [C]     Heat_Flux [W/m2]'
        write(*,'(A)') '  ----------------------------------------------------------'
        n_points = 40
        do i = 0, n_points
            log_dT = 3.0d0 * dble(i) / dble(n_points) ! dT_e from 10^0 (1) to 10^3 (1000)
            local_dT_e = 10.0d0**log_dT
            
            ! Piecewise calculation for profile
            if (local_dT_e < 5.0d0) then
                local_flux = 500.0d0 * local_dT_e**1.25d0
            else if (local_dT_e < 30.0d0) then
                local_flux = mu_l * h_fg * sqrt((g * (rho_l - rho_v)) / sigma) * &
                             ((Cp_l * local_dT_e) / (C_sf * h_fg * Pr_l))**3
                if (local_flux > q_max) local_flux = q_max
            else if (local_dT_e < 120.0d0) then
                log_dT = (log10(local_dT_e) - log10(30.0d0)) / (log10(120.0d0) - log10(30.0d0))
                local_flux = 10.0d0**(log10(q_max) + log_dT * (log10(q_min) - log10(q_max)))
            else
                h_fg_mod = h_fg + 0.8d0 * Cp_v * local_dT_e
                h_film = 0.62d0 * ((g * rho_v * (rho_l - rho_v) * k_v**3 * h_fg_mod) / &
                                   (mu_v * dim1 * local_dT_e))**0.25d0
                h_rad = (emiss * sigma_sb * (((T_sat + local_dT_e + 273.15d0)**4) - &
                        (T_sat + 273.15d0)**4)) / local_dT_e
                local_flux = (h_film + 0.75d0 * h_rad) * local_dT_e
                if (local_flux < q_min) local_flux = q_min
            end if

            write(*,'(F10.4,4X,ES14.6)') local_dT_e, local_flux
        end do

    else
        ! ============================================================
        ! FILM CONDENSATION ANALYSIS
        ! ============================================================
        dT = T_sat - T_s
        if (dT <= 0.0d0) then
            write(*,*) 'ERROR: Surface must be colder than saturation temp for condensation.'
            stop
        end if

        h_fg_mod = h_fg + 0.68d0 * Cp_l * dT

        select case (geom_type)
        case (1)
            regime_name = "Vertical Plate Condensation"
            L_c = dim1
            if (dim2 <= 0.0d0) dim2 = 1.0d0 ! width default
            A_s = dim1 * dim2
            h_avg = 0.943d0 * ((g * rho_l * (rho_l - rho_v) * k_l**3 * h_fg_mod) / &
                               (mu_l * L_c * dT))**0.25d0
        case (2)
            regime_name = "Horizontal Cylinder Condensation"
            L_c = dim1 ! Diameter
            if (dim2 <= 0.0d0) dim2 = 1.0d0 ! cylinder length
            A_s = PI * dim1 * dim2
            h_avg = 0.729d0 * ((g * rho_l * (rho_l - rho_v) * k_l**3 * h_fg_mod) / &
                               (mu_l * L_c * dT))**0.25d0
        case (3)
            regime_name = "Sphere Condensation"
            L_c = dim1
            A_s = PI * dim1**2
            h_avg = 0.815d0 * ((g * rho_l * (rho_l - rho_v) * k_l**3 * h_fg_mod) / &
                               (mu_l * L_c * dT))**0.25d0
        case default
            write(*,*) 'ERROR: Invalid geometry type for condensation.'
            stop
        end select

        Q = h_avg * A_s * dT
        mass_flow = Q / h_fg_mod

        ! Output Condensation results
        write(*,'(A)') '============================================================'
        write(*,'(A)') '   FILM CONDENSATION HEAT TRANSFER SIZING'
        write(*,'(A)') '============================================================'
        write(*,*)
        write(*,'(A,A)')        '  Condensation Regime     = ', trim(regime_name)
        write(*,'(A,F12.4,A)')  '  Characteristic Length   = ', L_c, ' m'
        write(*,'(A,F12.4,A)')  '  Heat Transfer Area (As) = ', A_s, ' m2'
        write(*,'(A,F12.2,A)')  '  Saturation Temp (Tsat)  = ', T_sat, ' deg-C'
        write(*,'(A,F12.2,A)')  '  Surface Temperature (Ts)= ', T_s, ' deg-C'
        write(*,'(A,F12.2,A)')  '  Temp Difference (dT)    = ', dT, ' deg-C'
        write(*,*)
        write(*,'(A)') '--- HEAT TRANSFER RESULTS -----------------------------------'
        write(*,'(A,F12.4,A)')  '  Modified Latent Heat    = ', h_fg_mod / 1000.0d0, ' kJ/kg'
        write(*,'(A,F12.4,A)')  '  Average Condensation h  = ', h_avg, ' W/m2.K'
        write(*,'(A,F12.4,A)')  '  Total Heat Transfer (Q) = ', Q, ' W'
        write(*,'(A,ES12.4,A)') '  Condensation Mass Rate  = ', mass_flow, ' kg/s'
        write(*,*)

        ! Generate local boundary layer profile data (film thickness delta along x)
        write(*,'(A)') '--- LOCAL PROFILE ALONG SURFACE ----------------------------'
        write(*,'(A)') '  x [m]       film_thickness [mm]'
        write(*,'(A)') '  ----------------------------------------------------------'
        n_points = 40
        dx = L_c / dble(n_points)
        do i = 1, n_points
            x = dble(i) * dx
            ! Nusselt film thickness relation
            local_delta = ((4.0d0 * mu_l * k_l * dT * x) / (g * rho_l * (rho_l - rho_v) * h_fg_mod))**0.25d0
            write(*,'(F8.4,4X,F12.6)') x, local_delta * 1000.0d0
        end do
    end if

end program boiling_condensation