program moody_diagram
    implicit none
    integer :: i, iostat_val
    double precision :: Re, epsD, D, rho, mu, f_cw, f_sj, V, dPperL, nu
    double precision :: re_i, f_lam_i, f_cw_i, f_smooth_i
    character(len=80) :: regime

    read(*,*,iostat=iostat_val) Re
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid Reynolds number input.'
        stop
    end if
    read(*,*,iostat=iostat_val) epsD
    read(*,*,iostat=iostat_val) D
    read(*,*,iostat=iostat_val) rho
    read(*,*,iostat=iostat_val) mu
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all Moody diagram inputs.'
        stop
    end if
    if (Re <= 0.0d0) then
        write(*,*) 'ERROR: Reynolds number must be positive.'
        stop
    end if
    if (epsD < 0.0d0) then
        write(*,*) 'ERROR: Relative roughness cannot be negative.'
        stop
    end if

    if (Re < 2300.0d0) then
        f_cw = 64.0d0/Re
        f_sj = f_cw
        regime = 'Laminar flow (Re < 2300)'
    else if (Re < 4000.0d0) then
        f_cw = colebrook_solve(Re, epsD)
        f_sj = swamee_jain(Re, epsD)
        regime = 'Transition zone (2300 < Re < 4000)'
    else
        f_cw = colebrook_solve(Re, epsD)
        f_sj = swamee_jain(Re, epsD)
        if (epsD < 1.0d-8) then
            regime = 'Turbulent - hydraulically smooth'
        else
            regime = 'Turbulent - rough wall'
        end if
    end if

    V = 0.0d0
    dPperL = 0.0d0
    if (D > 0.0d0 .and. rho > 0.0d0 .and. mu > 0.0d0) then
        nu = mu/rho
        V = Re*nu/D
        dPperL = f_cw*rho*V*V/(2.0d0*D)
    end if

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   MOODY DIAGRAM EXPLORER ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- OPERATING POINT -----------------------------------------'
    write(*,'(A,ES12.4)')   '  Reynolds Number          = ', Re
    write(*,'(A,ES12.4)')   '  Relative Roughness       = ', epsD
    write(*,'(A,A)')        '  Flow Regime              = ', trim(regime)
    write(*,*)
    write(*,'(A)') '--- FRICTION FACTOR RESULTS ---------------------------------'
    write(*,'(A,ES12.6)')   '  Colebrook-White f        = ', f_cw
    write(*,'(A,ES12.6)')   '  Swamee-Jain f            = ', f_sj
    if (f_cw > 0.0d0) then
        write(*,'(A,ES12.4,A)') '  Swamee-Jain Error        = ', &
            abs(f_sj-f_cw)/f_cw*100.0d0, ' percent'
    end if
    write(*,*)
    write(*,'(A)') '--- PIPE FLOW ESTIMATES -------------------------------------'
    write(*,'(A,ES12.4,A)') '  Pipe Diameter            = ', D, ' m'
    write(*,'(A,ES12.4,A)') '  Fluid Density            = ', rho, ' kg/m3'
    write(*,'(A,ES12.4,A)') '  Dynamic Viscosity        = ', mu, ' Pa.s'
    write(*,'(A,ES12.4,A)') '  Bulk Velocity            = ', V, ' m/s'
    write(*,'(A,ES12.4,A)') '  Pressure Drop per Length = ', dPperL, ' Pa/m'
    write(*,*)

    write(*,'(A)') '--- MOODY CURVE DATA ----------------------------------------'
    write(*,'(A)') '  Re             f_laminar     f_colebrook   f_smooth'
    write(*,'(A)') '  ----------------------------------------------------------'
    do i=1,200
        re_i = 500.0d0 * (1.0d8/500.0d0)**(dble(i-1)/199.0d0)
        if (re_i < 2300.0d0) then
            f_lam_i = 64.0d0/re_i
            f_cw_i = 0.0d0
            f_smooth_i = 0.0d0
        else
            f_lam_i = 0.0d0
            f_cw_i = colebrook_solve(re_i, epsD)
            f_smooth_i = colebrook_solve(re_i, 0.0d0)
        end if
        write(*,'(ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') &
            re_i, f_lam_i, f_cw_i, f_smooth_i
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Colebrook-White: 1/sqrt(f) = -2 log10(epsD/3.7 + 2.51/(Re*sqrt(f))).'
    write(*,'(A)') '  Swamee-Jain: f = 0.25/[log10(epsD/3.7 + 5.74/Re^0.9)]^2.'
    write(*,'(A)') '  Laminar: f = 64/Re.'

contains
    double precision function colebrook_solve(Re_in, eD)
        implicit none
        double precision, intent(in) :: Re_in, eD
        integer :: it
        double precision :: f, rhs
        f = swamee_jain(Re_in, eD)
        do it=1,80
            rhs = -2.0d0*log10(eD/3.7d0 + 2.51d0/(Re_in*sqrt(f)))
            f = 1.0d0/(rhs*rhs)
        end do
        colebrook_solve = f
    end function colebrook_solve

    double precision function swamee_jain(Re_in, eD)
        implicit none
        double precision, intent(in) :: Re_in, eD
        double precision :: v
        v = log10(eD/3.7d0 + 5.74d0/(Re_in**0.9d0))
        swamee_jain = 0.25d0/(v*v)
    end function swamee_jain
end program moody_diagram