program sphere_convection
    implicit none
    double precision :: D, V, Ti, Ts, rho, mu, kf, Pr, cp, beta
    integer :: fl_type
    
    double precision :: Tf, dT, nu, Re, NuW, NuR, h, Q, As, g, pi, mu_s, temp_val
    integer :: iostat_val, i
    double precision :: Vs, Res, NuWs, NuRs, hs, Qs
    
    g = 9.80665d0
    pi = 3.141592653589793d0
    
    ! Read inputs sequentially
    read(*,*,iostat=iostat_val) D
    read(*,*,iostat=iostat_val) V
    read(*,*,iostat=iostat_val) Ti
    read(*,*,iostat=iostat_val) Ts
    read(*,*,iostat=iostat_val) fl_type
    
    ! Read custom parameters
    read(*,*,iostat=iostat_val) rho
    read(*,*,iostat=iostat_val) mu
    read(*,*,iostat=iostat_val) kf
    read(*,*,iostat=iostat_val) Pr
    read(*,*,iostat=iostat_val) cp
    read(*,*,iostat=iostat_val) temp_val ! mu_s from php
    
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all input values.'
        stop
    end if
    
    ! Evaluate bulk properties at free-stream temperature T_inf = Ti
    if (fl_type == 1) then
        ! Air
        rho = 1.177d0
        mu = 1.85d-5
        kf = 0.0263d0
        Pr = 0.71d0
        cp = 1007d0
        mu_s = 1.85d-5 * ((Ts + 273.15d0) / (Ti + 273.15d0))**0.7d0
    else if (fl_type == 2) then
        ! Water
        rho = 997d0
        mu = 8.9d-4
        kf = 0.613d0
        Pr = 6.13d0
        cp = 4180d0
        mu_s = 2.414d-5 * 10.0d0**(247.8d0 / (Ts + 273.15d0 - 140.0d0))
    else
        ! Engine Oil
        rho = 870d0
        mu = 0.05d0
        kf = 0.14d0
        Pr = 500d0
        cp = 2000d0
        mu_s = 0.05d0 * exp(-0.04d0 * (Ts - Ti))
    end if
    
    ! Override custom mu_s if provided
    if (temp_val > 0.0d0) mu_s = temp_val
    
    nu = mu / rho
    Re = rho * V * D / mu
    
    ! Nusselt Whitaker correlation (Eq. 7.56)
    NuW = 2.0d0 + (0.4d0 * Re**0.5d0 + 0.06d0 * Re**(2.0d0/3.0d0)) * Pr**0.4d0 * (mu / mu_s)**0.25d0
    
    ! Nusselt Ranz-Marshall correlation
    NuR = 2.0d0 + 0.6d0 * Re**0.5d0 * Pr**(1.0d0/3.0d0)
    
    h = NuW * kf / D
    As = pi * D**2
    dT = abs(Ts - Ti)
    Q = h * As * dT
    
    write(*,'(A)') '============================================'
    write(*,'(A)') '  CONVECTION OVER SPHERE ENGINE'
    write(*,'(A)') '============================================'
    write(*,*)
    write(*,'(A,ES14.4)') '  Reynolds Re             = ', Re
    write(*,'(A,F12.4)')    '  Whitaker Nu             = ', NuW
    write(*,'(A,F12.4)')    '  Ranz Nu                 = ', NuR
    write(*,'(A,F12.4,A)')  '  Coeff h                 = ', h, ' W/m2K'
    write(*,'(A,F12.4,A)')  '  Transfer Q              = ', Q, ' W'
    write(*,*)
    
    write(*,'(A)') '--- VELOCITY SWEEP ---'
    write(*,'(A)') '  V[m/s]     Re           Nu_Whitaker  Nu_Ranz      h[W/m2K]   Q[W]'
    write(*,'(A)') '  -------------------------------------------------------------------------'
    do i=1,25
        Vs = 0.05d0 + (max(V,0.2d0)*4.0d0 - 0.05d0)*dble(i-1)/24.0d0
        Res = rho * Vs * D / mu
        NuWs = 2.0d0 + (0.4d0 * Res**0.5d0 + 0.06d0 * Res**(2.0d0/3.0d0)) * Pr**0.4d0 * (mu / mu_s)**0.25d0
        NuRs = 2.0d0 + 0.6d0 * Res**0.5d0 * Pr**(1.0d0/3.0d0)
        hs = NuWs * kf / D
        Qs = hs * As * dT
        write(*,'(2X,F8.3,2X,ES10.3,2X,F11.3,2X,F11.3,2X,F10.4,2X,F12.4)') Vs, Res, NuWs, NuRs, hs, Qs
    end do
    write(*,*)
end program sphere_convection