program taylor_couette_flow
    implicit none
    integer :: iostat_val
    double precision :: R1_mm, R2_mm, H_len_mm, N1_rpm, N2_rpm, nu_cSt, rho_kgm3
    double precision :: R1_m, R2_m, H_m, d_m, eta_ratio, nu_m2s
    double precision :: Omega1_rads, Omega2_rads, U1_ms, Ta_num, Ta_crit, N1_crit_rpm
    double precision :: lambda_vortex_mm, num_vortex_pairs, torque_Nm, power_W
    character(len=32) :: flow_regime
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) R1_mm      ! Inner Radius [mm] (e.g. 50.0)
    read(*,*,iostat=iostat_val) R2_mm      ! Outer Radius [mm] (e.g. 60.0)
    read(*,*,iostat=iostat_val) H_len_mm   ! Cylinder Length/Height [mm] (e.g. 200.0)
    read(*,*,iostat=iostat_val) N1_rpm     ! Inner Speed [RPM] (e.g. 150.0)
    read(*,*,iostat=iostat_val) N2_rpm     ! Outer Speed [RPM] (e.g. 0.0)
    read(*,*,iostat=iostat_val) nu_cSt     ! Kinematic Viscosity [cSt] (e.g. 1.0)
    read(*,*,iostat=iostat_val) rho_kgm3   ! Fluid Density [kg/m3] (e.g. 998.0)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for Taylor-Couette calculation.'
        stop
    end if

    if (R1_mm <= 0.0d0 .or. R2_mm <= R1_mm .or. H_len_mm <= 0.0d0 .or. nu_cSt <= 0.0d0) then
        write(*,*) 'ERROR: Radii must satisfy R2 > R1 > 0, height and viscosity > 0.'
        stop
    end if

    R1_m = R1_mm * 1.0d-3
    R2_m = R2_mm * 1.0d-3
    H_m = H_len_mm * 1.0d-3
    d_m = R2_m - R1_m
    eta_ratio = R1_m / R2_m
    nu_m2s = nu_cSt * 1.0d-6

    Omega1_rads = (N1_rpm * 2.0d0 * PI) / 60.0d0
    Omega2_rads = (N2_rpm * 2.0d0 * PI) / 60.0d0
    U1_ms = Omega1_rads * R1_m

    ! Taylor Number
    Ta_num = (2.0d0 * (Omega1_rads**2) * R1_m * (d_m**3)) / ((nu_m2s**2) * (1.0d0 + eta_ratio))

    ! Critical Taylor Number
    Ta_crit = 1708.0d0 / (1.0d0 - 0.652d0 * (d_m / R1_m))
    if (Ta_crit < 1695.0d0) Ta_crit = 1708.0d0

    N1_crit_rpm = (sqrt((Ta_crit * (nu_m2s**2) * (1.0d0 + eta_ratio)) / &
                   (2.0d0 * R1_m * (d_m**3))) * 60.0d0) / (2.0d0 * PI)

    ! Flow Regime Classification
    if (Ta_num < Ta_crit) then
        flow_regime = 'Circular Couette Flow (Laminar)'
    else if (Ta_num < 1.2d0 * Ta_crit) then
        flow_regime = 'Taylor Vortex Flow (TVF Cells)'
    else if (Ta_num < 8.0d0 * Ta_crit) then
        flow_regime = 'Wavy Vortex Flow (WVF Mode)'
    else if (Ta_num < 40.0d0 * Ta_crit) then
        flow_regime = 'Modulated Wavy Vortex'
    else
        flow_regime = 'Turbulent Taylor-Couette Flow'
    end if

    lambda_vortex_mm = 2.0d0 * (d_m * 1000.0d0)
    num_vortex_pairs = H_len_mm / lambda_vortex_mm

    ! Viscous torque on inner cylinder [N.m] (Laminar baseline)
    torque_Nm = (4.0d0 * PI * (rho_kgm3 * nu_m2s) * (R1_m**2) * (R2_m**2) * (Omega1_rads - Omega2_rads) * H_m) / &
                ((R2_m**2) - (R1_m**2))
    if (Ta_num > Ta_crit) then
        ! Torque enhancement due to Taylor vortex circulation
        torque_Nm = torque_Nm * (1.0d0 + 0.15d0 * sqrt(max(0.0d0, Ta_num / Ta_crit - 1.0d0)))
    end if
    power_W = torque_Nm * Omega1_rads

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — TAYLOR-COUETTE STABILITY ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Radii: R1 / R2 (Gap d)    = ', R1_mm, ' mm / ', R2_mm, ' mm (d = ', d_m*1000.0d0, ' mm)'
    write(*,'(A,F10.2,A,F10.2,A)') 'Speed N1 / Critical N1_c  = ', N1_rpm, ' RPM / ', N1_crit_rpm, ' RPM'
    write(*,'(A,F10.3,A,F10.2,A)') 'Radius Ratio eta / Length = ', eta_ratio, ' / ', H_len_mm, ' mm'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F12.1)')   'TAYLOR NUMBER (Ta)        = ', Ta_num
    write(*,'(A,F12.1)')   'Critical Taylor (Ta_c)    = ', Ta_crit
    write(*,'(A,A)')       'Hydrodynamic Regime       = ', trim(flow_regime)
    write(*,'(A,F10.2,A,F6.1,A)') 'Vortex Pair Wavelength    = ', lambda_vortex_mm, ' mm (', num_vortex_pairs, ' pairs)'
    write(*,'(A,F10.4,A,F10.2,A)')'Viscous Torque / Power    = ', torque_Nm, ' N.m / ', power_W, ' Watts'
    write(*,'(A)') '============================================================'

end program taylor_couette_flow