program fluid_statics
    implicit none

    ! Gravitational acceleration [m/s2]
    double precision, parameter :: g = 9.80665d0

    integer :: mode
    
    ! Mode A variables (Layered Fluids)
    double precision :: P0, z_interest, P_final
    integer :: N_layers
    double precision, allocatable :: rho_layer(:), h_layer(:), p_interface(:)
    double precision :: current_depth, accum_pressure
    integer :: i, k

    ! Mode B variables (Submerged Surfaces)
    integer :: shape
    double precision :: dim1, dim2, dim3, dim4, dim5
    double precision :: theta, hc, rho_fluid, P_atm
    double precision :: theta_rad, area, Ixx_c, yc, ycp, hcp, FR_gauge, FR_abs
    double precision :: shape_centroid_dist

    ! Mode C variables (Manometers)
    integer :: N_left, N_right
    double precision, allocatable :: rho_L(:), h_L(:)
    double precision, allocatable :: rho_R(:), h_R(:)
    double precision :: rho_m, dh, sum_L, sum_R, delta_P

    ! Read mode
    read(*,*,iostat=i) mode
    if (i /= 0) then
        print *, "ERROR: Invalid calculation mode."
        stop
    end if

    select case(mode)
    case(1)
        ! ==========================================
        ! MODE 1: HYDROSTATIC PRESSURE (LAYERED)
        ! ==========================================
        read(*,*,iostat=i) P0
        if (i /= 0) then
            print *, "ERROR: Failed to read surface pressure P0."
            stop
        end if

        read(*,*,iostat=i) N_layers
        if (i /= 0 .or. N_layers <= 0 .or. N_layers > 20) then
            print *, "ERROR: Invalid number of layers."
            stop
        end if

        allocate(rho_layer(N_layers))
        allocate(h_layer(N_layers))
        allocate(p_interface(0:N_layers))

        do i = 1, N_layers
            read(*,*,iostat=k) rho_layer(i), h_layer(i)
            if (k /= 0 .or. rho_layer(i) < 0.0d0 .or. h_layer(i) <= 0.0d0) then
                print *, "ERROR: Invalid layer data for layer", i
                stop
            end if
        end do

        read(*,*,iostat=i) z_interest
        if (i /= 0 .or. z_interest < 0.0d0) then
            print *, "ERROR: Invalid depth of interest."
            stop
        end if

        ! Compute interface pressures
        p_interface(0) = P0
        accum_pressure = P0
        do i = 1, N_layers
            accum_pressure = accum_pressure + (rho_layer(i) * g * h_layer(i)) / 1000.0d0  ! in kPa
            p_interface(i) = accum_pressure
        end do

        ! Compute pressure at z_interest
        P_final = P0
        current_depth = 0.0d0
        do i = 1, N_layers
            if (z_interest <= current_depth + h_layer(i)) then
                ! Target lies in current layer
                P_final = p_interface(i-1) + (rho_layer(i) * g * (z_interest - current_depth)) / 1000.0d0
                exit
            else
                current_depth = current_depth + h_layer(i)
                if (i == N_layers) then
                    ! Depth is beyond all layers, extrapolate with bottom layer density
                    P_final = p_interface(N_layers) + (rho_layer(N_layers) * g * (z_interest - current_depth)) / 1000.0d0
                end if
            end if
        end do

        ! Print results
        print *, "--- RESULTS MODE A ---"
        print '(A,F12.4,A)', "Surface Pressure (P0)      = ", P0, " kPa"
        print '(A,F12.4,A)', "Target Depth               = ", z_interest, " m"
        print '(A,F12.4,A)', "Pressure at Target Depth   = ", P_final, " kPa"
        print *
        print *, "--- INTERFACE PRESSURE DISTRIBUTION ---"
        print '(A12, A15, A18)', "Boundary", "Depth [m]", "Pressure [kPa]"
        current_depth = 0.0d0
        print '(I12, F15.3, F18.4)', 0, current_depth, p_interface(0)
        do i = 1, N_layers
            current_depth = current_depth + h_layer(i)
            print '(I12, F15.3, F18.4)', i, current_depth, p_interface(i)
        end do

    case(2)
        ! ==========================================
        ! MODE 2: SUBMERGED PLANE SURFACES
        ! ==========================================
        read(*,*,iostat=i) shape
        if (i /= 0 .or. shape < 1 .or. shape > 4) then
            print *, "ERROR: Invalid shape code (1=Rect, 2=Circ, 3=Tri, 4=Trap)."
            stop
        end if

        read(*,*,iostat=i) dim1, dim2, dim3, dim4, dim5
        if (i /= 0) then
            print *, "ERROR: Failed to read dimensions."
            stop
        end if

        read(*,*,iostat=i) theta
        if (i /= 0 .or. theta < 0.0d0 .or. theta > 90.0d0) then
            print *, "ERROR: Inclination angle must be in [0, 90] degrees."
            stop
        end if

        read(*,*,iostat=i) hc
        if (i /= 0 .or. hc < 0.0d0) then
            print *, "ERROR: Centroid depth must be non-negative."
            stop
        end if

        read(*,*,iostat=i) rho_fluid
        if (i /= 0 .or. rho_fluid <= 0.0d0) then
            print *, "ERROR: Fluid density must be positive."
            stop
        end if

        read(*,*,iostat=i) P_atm
        if (i /= 0 .or. P_atm < 0.0d0) then
            P_atm = 101.325d0 ! Default standard atm in kPa
        end if

        ! Calculate shape properties
        select case(shape)
        case(1)
            ! Rectangle: dim1 = width b, dim2 = height h
            if (dim1 <= 0.0d0 .or. dim2 <= 0.0d0) then
                print *, "ERROR: Rectangle dimensions must be positive."
                stop
            end if
            area = dim1 * dim2
            Ixx_c = (dim1 * dim2**3) / 12.0d0
        case(2)
            ! Circle: dim1 = Radius R
            if (dim1 <= 0.0d0) then
                print *, "ERROR: Circle radius must be positive."
                stop
            end if
            area = 3.141592653589793d0 * dim1**2
            Ixx_c = (3.141592653589793d0 * dim1**4) / 4.0d0
        case(3)
            ! Triangle: dim1 = base b, dim2 = height h
            if (dim1 <= 0.0d0 .or. dim2 <= 0.0d0) then
                print *, "ERROR: Triangle dimensions must be positive."
                stop
            end if
            area = 0.5d0 * dim1 * dim2
            Ixx_c = (dim1 * dim2**3) / 36.0d0
        case(4)
            ! Trapezoid: dim1 = top base a, dim2 = bottom base b, dim3 = height h
            if (dim1 <= 0.0d0 .or. dim2 <= 0.0d0 .or. dim3 <= 0.0d0) then
                print *, "ERROR: Trapezoid dimensions must be positive."
                stop
            end if
            area = 0.5d0 * (dim1 + dim2) * dim3
            Ixx_c = (dim3**3 * (dim1**2 + 4.0d0 * dim1 * dim2 + dim2**2)) / (36.0d0 * (dim1 + dim2))
        end select

        theta_rad = theta * 3.141592653589793d0 / 180.0d0

        ! Resultant Forces
        ! Force = pressure at centroid * Area
        FR_gauge = (rho_fluid * g * hc * area) / 1000.0d0  ! in kN
        FR_abs = ((rho_fluid * g * hc + P_atm * 1000.0d0) * area) / 1000.0d0 ! in kN

        ! Center of Pressure
        if (theta < 0.001d0) then
            ! Horizontal surface
            yc = hc
            ycp = yc
            hcp = hc
        else
            yc = hc / sin(theta_rad)
            ycp = yc + Ixx_c / (yc * area)
            hcp = ycp * sin(theta_rad)
        end if

        ! Print results
        print *, "--- RESULTS MODE B ---"
        print '(A,I2)',    "Shape Code                 = ", shape
        print '(A,F12.4,A)', "Surface Area (A)           = ", area, " m2"
        print '(A,E14.6,A)', "Moment of Inertia (Ixx,c)  = ", Ixx_c, " m4"
        print '(A,F12.4,A)', "Inclination Angle (theta)  = ", theta, " deg"
        print '(A,F12.4,A)', "Centroid Depth (hc)        = ", hc, " m"
        print '(A,F12.4,A)', "Centroid distance (yc)     = ", yc, " m"
        print '(A,F12.4,A)', "Resultant Force (Gauge)    = ", FR_gauge, " kN"
        print '(A,F12.4,A)', "Resultant Force (Absolute) = ", FR_abs, " kN"
        print '(A,F12.4,A)', "Center of Pressure (ycp)   = ", ycp, " m"
        print '(A,F12.4,A)', "Center of Pressure depth(h)= ", hcp, " m"

    case(3)
        ! ==========================================
        ! MODE 3: MANOMETER DIFFERENTIAL PRESSURE
        ! ==========================================
        read(*,*,iostat=i) N_left
        if (i /= 0 .or. N_left < 0 .or. N_left > 10) then
            print *, "ERROR: Invalid number of fluids in left leg."
            stop
        end if

        allocate(rho_L(N_left))
        allocate(h_L(N_left))

        do i = 1, N_left
            read(*,*,iostat=k) rho_L(i), h_L(i)
            if (k /= 0 .or. rho_L(i) < 0.0d0 .or. h_L(i) < 0.0d0) then
                print *, "ERROR: Invalid left leg fluid data at index", i
                stop
            end if
        end do

        read(*,*,iostat=i) N_right
        if (i /= 0 .or. N_right < 0 .or. N_right > 10) then
            print *, "ERROR: Invalid number of fluids in right leg."
            stop
        end if

        allocate(rho_R(N_right))
        allocate(h_R(N_right))

        do i = 1, N_right
            read(*,*,iostat=k) rho_R(i), h_R(i)
            if (k /= 0 .or. rho_R(i) < 0.0d0 .or. h_R(i) < 0.0d0) then
                print *, "ERROR: Invalid right leg fluid data at index", i
                stop
            end if
        end do

        read(*,*,iostat=i) rho_m, dh
        if (i /= 0 .or. rho_m <= 0.0d0) then
            print *, "ERROR: Invalid manometer fluid density."
            stop
        end if

        ! Calculate stages
        sum_L = 0.0d0
        do i = 1, N_left
            sum_L = sum_L + (rho_L(i) * g * h_L(i)) / 1000.0d0  ! in kPa
        end do

        sum_R = 0.0d0
        do i = 1, N_right
            sum_R = sum_R + (rho_R(i) * g * h_R(i)) / 1000.0d0  ! in kPa
        end do

        ! Delta P = P_left - P_right
        ! P_left + sum_L = P_right + sum_R + rho_m*g*dh
        ! Delta P = sum_R + (rho_m*g*dh)/1000 - sum_L
        delta_P = sum_R + (rho_m * g * dh) / 1000.0d0 - sum_L

        ! Print results
        print *, "--- RESULTS MODE C ---"
        print '(A,F12.4,A)', "Left Leg Fluids head       = ", sum_L, " kPa"
        print '(A,F12.4,A)', "Right Leg Fluids head      = ", sum_R, " kPa"
        print '(A,F12.4,A)', "Manometer deflection head  = ", (rho_m * g * dh) / 1000.0d0, " kPa"
        print '(A,F12.4,A)', "Pressure Difference (DP)   = ", delta_P, " kPa"
        print '(A,F12.4,A)', "Pressure Difference (mH2O) = ", delta_P / (9.80665d0), " mH2O"

    case default
        print *, "ERROR: Calculator mode not recognized."
    end select

end program fluid_statics