program two_phase_flow
    implicit none
    integer :: i, iostat_val, n_points
    double precision :: D_mm, L_m, eps_mm, theta_deg, mdot_total, x_quality
    double precision :: rho_L, mu_L, sigma, rho_G, mu_G
    double precision :: D_m, eps_m, theta_rad, Area, G_flux
    double precision :: j_L, j_G, j_tot, Re_L0, Re_G0, f_L0, f_G0
    double precision :: rho_H, Fr_H, We_L, E_fr, F_fr, H_fr, phi2_friedel
    double precision :: X_tt, phi2_lockhart, alpha_void, HL_holdup, rho_2phi
    double precision :: dp_dz_L0, dp_fric_Pa, dp_grav_Pa, dp_total_Pa
    double precision :: x_var, x_tt_var, p2_var, void_var
    double precision, parameter :: PI = 3.141592653589793d0
    double precision, parameter :: g = 9.80665d0
    character(len=64) :: regime_name

    ! Read inputs
    read(*,*,iostat=iostat_val) D_mm
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid pipe diameter.'
        stop
    end if
    read(*,*,iostat=iostat_val) L_m
    read(*,*,iostat=iostat_val) eps_mm
    read(*,*,iostat=iostat_val) theta_deg
    read(*,*,iostat=iostat_val) mdot_total
    read(*,*,iostat=iostat_val) x_quality
    read(*,*,iostat=iostat_val) rho_L
    read(*,*,iostat=iostat_val) mu_L
    read(*,*,iostat=iostat_val) sigma
    read(*,*,iostat=iostat_val) rho_G
    read(*,*,iostat=iostat_val) mu_G

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all two-phase flow inputs.'
        stop
    end if

    ! Validation
    if (D_mm <= 0.0d0 .or. L_m <= 0.0d0 .or. mdot_total <= 0.0d0) then
        write(*,*) 'ERROR: Diameter, length, and mass flow rate must be positive.'
        stop
    end if
    if (x_quality < 0.0d0 .or. x_quality > 1.0d0) then
        write(*,*) 'ERROR: Vapor quality x must be between 0.0 and 1.0.'
        stop
    end if
    if (rho_L <= 0.0d0 .or. mu_L <= 0.0d0 .or. rho_G <= 0.0d0 .or. mu_G <= 0.0d0) then
        write(*,*) 'ERROR: Liquid and gas properties must be positive.'
        stop
    end if

    ! Conversions
    D_m = D_mm * 1.0d-3
    eps_m = eps_mm * 1.0d-3
    theta_rad = theta_deg * PI / 180.0d0
    Area = (PI / 4.0d0) * (D_m**2)
    G_flux = mdot_total / Area ! kg/(m2.s)

    ! Superficial Velocities
    j_L = (G_flux * (1.0d0 - x_quality)) / rho_L
    j_G = (G_flux * x_quality) / rho_G
    j_tot = j_L + j_G

    ! Single-phase All-Liquid and All-Gas Reynolds & Friction factors
    Re_L0 = (G_flux * D_m) / mu_L
    Re_G0 = (G_flux * D_m) / mu_G

    ! Churchill friction factor for all-liquid
    f_L0 = churchill_friction(Re_L0, eps_m / D_m)
    f_G0 = churchill_friction(Re_G0, eps_m / D_m)

    ! Single-Phase Liquid Pressure Gradient (dP/dz)_L0
    dp_dz_L0 = f_L0 * (G_flux**2) / (2.0d0 * rho_L * D_m)

    ! Homogeneous Density & Non-dimensional groups for Friedel
    if (x_quality > 0.0d0 .and. x_quality < 1.0d0) then
        rho_H = 1.0d0 / ((x_quality / rho_G) + ((1.0d0 - x_quality) / rho_L))
    else if (x_quality <= 0.0d0) then
        rho_H = rho_L
    else
        rho_H = rho_G
    end if

    Fr_H = (G_flux**2) / (g * D_m * (rho_H**2))
    if (sigma > 0.0d0) then
        We_L = (G_flux**2 * D_m) / (sigma * rho_H)
    else
        We_L = 100.0d0
    end if

    ! Friedel Two-Phase Multiplier Phi_L0^2
    E_fr = ((1.0d0 - x_quality)**2) + (x_quality**2) * ((rho_L * f_G0) / (rho_G * f_L0))
    F_fr = (x_quality**0.78d0) * ((1.0d0 - x_quality)**0.224d0)
    H_fr = ((rho_L / rho_G)**0.91d0) * ((mu_G / mu_L)**0.19d0) * ((1.0d0 - (mu_G / mu_L))**0.7d0)
    
    if (x_quality > 0.0d0 .and. x_quality < 1.0d0) then
        phi2_friedel = E_fr + (3.24d0 * F_fr * H_fr) / ((max(1.0d-4, Fr_H)**0.045d0) * (max(1.0d-4, We_L)**0.035d0))
    else
        phi2_friedel = 1.0d0
    end if

    ! Lockhart-Martinelli Parameter X_tt
    if (x_quality > 0.001d0 .and. x_quality < 0.999d0) then
        X_tt = (((1.0d0 - x_quality) / x_quality)**0.9d0) * &
               ((rho_G / rho_L)**0.5d0) * ((mu_L / mu_G)**0.1d0)
        phi2_lockhart = 1.0d0 + (20.0d0 / max(1.0d-3, X_tt)) + (1.0d0 / max(1.0d-6, X_tt**2))
    else if (x_quality <= 0.001d0) then
        X_tt = 1000.0d0
        phi2_lockhart = 1.0d0
    else
        X_tt = 0.001d0
        phi2_lockhart = (rho_L * f_G0) / (rho_G * f_L0)
    end if

    ! Void Fraction (Rouhani-Axelsson Drift Flux Model)
    if (x_quality > 0.0d0 .and. x_quality < 1.0d0) then
        alpha_void = (x_quality / rho_G) / &
            ((1.0d0 + 0.12d0 * (1.0d0 - x_quality)) * ((x_quality / rho_G) + ((1.0d0 - x_quality) / rho_L)) + &
             (1.18d0 * (1.0d0 - x_quality) * ((g * sigma * (rho_L - rho_G))**0.25d0)) / (G_flux * sqrt(rho_L)))
        alpha_void = min(0.999d0, max(0.001d0, alpha_void))
    else if (x_quality <= 0.0d0) then
        alpha_void = 0.0d0
    else
        alpha_void = 1.0d0
    end if

    HL_holdup = 1.0d0 - alpha_void
    rho_2phi = (alpha_void * rho_G) + (HL_holdup * rho_L)

    ! Pressure Drops over Length L
    dp_fric_Pa = phi2_friedel * dp_dz_L0 * L_m
    dp_grav_Pa = rho_2phi * g * L_m * sin(theta_rad)
    dp_total_Pa = dp_fric_Pa + dp_grav_Pa

    ! Flow Regime Identification
    if (alpha_void < 0.20d0 .and. j_L > 0.5d0) then
        regime_name = 'Bubbly Flow'
    else if (alpha_void < 0.70d0 .and. j_tot < 4.0d0) then
        if (abs(theta_deg) < 15.0d0) then
            regime_name = 'Stratified / Plug Flow'
        else
            regime_name = 'Intermittent / Slug Flow'
        end if
    else if (alpha_void >= 0.70d0 .and. x_quality < 0.85d0) then
        regime_name = 'Annular-Mist Flow'
    else
        regime_name = 'Dispersed Mist Flow'
    end if

    ! Output Formatted Report
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — TWO-PHASE PIPE FLOW SOLVER'
    write(*,'(A)') '============================================================'
    write(*,'(A,ES12.4,A)') 'Pipe Diameter D           = ', D_mm, ' mm'
    write(*,'(A,ES12.4,A)') 'Pipe Length L             = ', L_m, ' m'
    write(*,'(A,ES12.4,A)') 'Inclination Angle         = ', theta_deg, ' deg'
    write(*,'(A,ES12.4,A)') 'Total Mass Flow mdot      = ', mdot_total, ' kg/s'
    write(*,'(A,ES12.4,A)') 'Mass Flux G               = ', G_flux, ' kg/(m2.s)'
    write(*,'(A,ES12.4)')   'Vapor/Gas Quality x       = ', x_quality
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,ES12.4,A)') 'Superficial Liquid Vel jL = ', j_L, ' m/s'
    write(*,'(A,ES12.4,A)') 'Superficial Gas Vel jG    = ', j_G, ' m/s'
    write(*,'(A,ES12.4)')   'All-Liquid Reynolds ReL0  = ', Re_L0
    write(*,'(A,ES12.4)')   'All-Gas Reynolds ReG0     = ', Re_G0
    write(*,'(A,ES12.4)')   'Lockhart-Martinelli Xtt   = ', X_tt
    write(*,'(A,ES12.4)')   'Friedel Multiplier Phi2L0 = ', phi2_friedel
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,ES12.4)')   'Void Fraction alpha       = ', alpha_void
    write(*,'(A,ES12.4)')   'Liquid Holdup HL          = ', HL_holdup
    write(*,'(A,ES12.4,A)') 'In-situ Mixture Density   = ', rho_2phi, ' kg/m3'
    write(*,'(A,A)')        'Predicted Flow Regime     = ', trim(regime_name)
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,ES12.4,A)') 'Frictional Pressure Drop  = ', dp_fric_Pa * 1.0d-3, ' kPa'
    write(*,'(A,ES12.4,A)') 'Hydrostatic Pressure Drop = ', dp_grav_Pa * 1.0d-3, ' kPa'
    write(*,'(A,ES12.4,A)') 'Total Two-Phase Delta-P   = ', dp_total_Pa * 1.0d-3, ' kPa'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- TWO-PHASE MULTIPLIER AND VOID FRACTION VS QUALITY ---'
    write(*,'(A)') '  Quality_x    Phi2_Friedel  Void_Fraction'
    n_points = 20
    do i = 1, n_points
        x_var = dble(i) / dble(n_points + 1)
        E_fr = ((1.0d0 - x_var)**2) + (x_var**2) * ((rho_L * f_G0) / (rho_G * f_L0))
        F_fr = (x_var**0.78d0) * ((1.0d0 - x_var)**0.224d0)
        p2_var = E_fr + (3.24d0 * F_fr * H_fr) / ((max(1.0d-4, Fr_H)**0.045d0) * (max(1.0d-4, We_L)**0.035d0))
        void_var = (x_var / rho_G) / &
            ((1.0d0 + 0.12d0 * (1.0d0 - x_var)) * ((x_var / rho_G) + ((1.0d0 - x_var) / rho_L)) + &
             (1.18d0 * (1.0d0 - x_var) * ((g * sigma * (rho_L - rho_G))**0.25d0)) / (G_flux * sqrt(rho_L)))
        write(*,'(2X,F10.3,2X,F12.4,2X,F12.4)') x_var, max(1.0d0, p2_var), min(1.0d0, max(0.0d0, void_var))
    end do

contains

    function churchill_friction(Re, rel_rough) result(f)
        double precision, intent(in) :: Re, rel_rough
        double precision :: f, A_term, B_term
        if (Re < 1.0d-4) then
            f = 0.02d0
            return
        end if
        A_term = (2.457d0 * log(1.0d0 / ((7.0d0 / Re)**0.9d0 + 0.27d0 * rel_rough)))**16
        B_term = (37530.0d0 / Re)**16
        f = 8.0d0 * (( (8.0d0 / Re)**12 + 1.0d0 / ((A_term + B_term)**1.5d0) )**(1.0d0 / 12.0d0))
    end function churchill_friction

end program two_phase_flow
