program pipe_support_span
    implicit none
    integer :: iostat_val
    double precision :: Do_mm, wt_mm, Di_mm, E_GPa, Sh_MPa, S_allow_MPa
    double precision :: rho_pipe_kgm3, rho_fluid_kgm3, tins_mm, rho_ins_kgm3, W_conc_N
    double precision :: delta_allow_mm, Do_m, Di_m, Dins_m, wt_m
    double precision :: As_m2, Af_m2, A_ins_m2, I_m4, Z_m3
    double precision :: w_pipe_Nm, w_fluid_Nm, w_ins_Nm, w_tot_Nm
    double precision :: L_defl_m, L_stress_m, L_max_m
    double precision, parameter :: PI = 3.141592653589793d0, g = 9.80665d0

    ! Read inputs
    read(*,*,iostat=iostat_val) Do_mm           ! Outside Diameter [mm]
    read(*,*,iostat=iostat_val) wt_mm           ! Wall Thickness [mm]
    read(*,*,iostat=iostat_val) E_GPa           ! Elastic Modulus E [GPa] (e.g. 200 for steel)
    read(*,*,iostat=iostat_val) Sh_MPa          ! Material Hot Allowable Stress Sh [MPa]
    read(*,*,iostat=iostat_val) rho_fluid_kgm3  ! Fluid Density [kg/m3] (1000=water, 0=gas)
    read(*,*,iostat=iostat_val) tins_mm         ! Insulation Thickness [mm]
    read(*,*,iostat=iostat_val) rho_ins_kgm3    ! Insulation Density [kg/m3] (e.g. 150)
    read(*,*,iostat=iostat_val) delta_allow_mm  ! Allowable Sag Deflection [mm] (std 2.5 mm)

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

    rho_pipe_kgm3 = 7850.0d0 ! Carbon steel
    S_allow_MPa = 0.23d0 * Sh_MPa ! Allowable weight bending stress per ASME B31.1/B31.3

    Do_m = Do_mm * 1.0d-3
    wt_m = wt_mm * 1.0d-3
    Di_m = max(1.0d-3, Do_m - 2.0d0 * wt_m)
    Dins_m = Do_m + 2.0d0 * (tins_mm * 1.0d-3)

    ! Cross-sectional properties
    As_m2 = (PI / 4.0d0) * (Do_m**2 - Di_m**2)
    Af_m2 = (PI / 4.0d0) * (Di_m**2)
    A_ins_m2 = (PI / 4.0d0) * (Dins_m**2 - Do_m**2)

    I_m4 = (PI / 64.0d0) * (Do_m**4 - Di_m**4)
    Z_m3 = (2.0d0 * I_m4) / Do_m

    ! Weight per unit meter [N/m]
    w_pipe_Nm = As_m2 * rho_pipe_kgm3 * g
    w_fluid_Nm = Af_m2 * rho_fluid_kgm3 * g
    w_ins_Nm = A_ins_m2 * rho_ins_kgm3 * g
    w_tot_Nm = w_pipe_Nm + w_fluid_Nm + w_ins_Nm

    ! Deflection-limited span L_defl (delta = 5 w L^4 / 384 E I)
    L_defl_m = ((384.0d0 * (E_GPa * 1.0d9) * I_m4 * (delta_allow_mm * 1.0d-3)) / (5.0d0 * w_tot_Nm))**0.25d0

    ! Bending stress-limited span L_stress (sigma = w L^2 / (8 Z) <= S_allow)
    L_stress_m = sqrt((8.0d0 * Z_m3 * (S_allow_MPa * 1.0d6)) / w_tot_Nm)

    ! Maximum allowable span
    L_max_m = min(L_defl_m, L_stress_m)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — ASME B31.1 / B31.3 PIPE SPAN ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A)')  'Pipe Outside Diameter     = ', Do_mm, ' mm'
    write(*,'(A,F10.2,A)')  'Pipe Wall Thickness       = ', wt_mm, ' mm'
    write(*,'(A,F10.2,A)')  'Total Distributed Weight  = ', w_tot_Nm, ' N/m'
    write(*,'(A,F10.2,A)')  'Allowable Deflection Sag  = ', delta_allow_mm, ' mm'
    write(*,'(A,F10.2,A)')  'Allowable Bending Stress  = ', S_allow_MPa, ' MPa'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'Deflection-Limited Span   = ', L_defl_m, ' m'
    write(*,'(A,F10.2,A)')  'Stress-Limited Span       = ', L_stress_m, ' m'
    write(*,'(A,F10.2,A)')  'GOVERNING MAXIMUM SPAN    = ', L_max_m, ' m'
    write(*,'(A,F10.2,A)')  'Maximum Span in Feet      = ', L_max_m * 3.28084d0, ' ft'
    write(*,'(A)') '============================================================'

end program pipe_support_span