program inclined_plate
    implicit none
    double precision :: Lp, Wp, th, Ts, Ta, rho, mu, kf, Pr, cp, beta
    integer :: config_type, fl_type
    
    double precision :: Tf, dT, nu_v, alpha, Lc, As, Ra, Nu, h, Q, g, pi, g_eff, temp_val
    integer :: iostat_val, i
    double precision :: dTs, Ras, Nus, hs, Qs
    
    g = 9.80665d0
    pi = 3.141592653589793d0
    
    ! Read inputs sequentially
    read(*,*,iostat=iostat_val) Lp
    read(*,*,iostat=iostat_val) Wp
    read(*,*,iostat=iostat_val) th
    read(*,*,iostat=iostat_val) Ts
    read(*,*,iostat=iostat_val) Ta
    read(*,*,iostat=iostat_val) config_type
    read(*,*,iostat=iostat_val) fl_type
    
    ! Read custom parameters (hardcoded as 0 from php except beta)
    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 ! beta from php
    
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all input values.'
        stop
    end if
    
    Tf = (Ts + Ta) / 2.0d0
    dT = abs(Ts - Ta)
    
    ! Select defaults based on fluid type
    if (fl_type == 1) then
        ! Air
        rho = 1.177d0
        mu = 1.85d-5
        kf = 0.0263d0
        Pr = 0.71d0
        cp = 1007d0
        beta = 1.0d0 / (Tf + 273.15d0)
    else if (fl_type == 2) then
        ! Water
        rho = 997d0
        mu = 8.9d-4
        kf = 0.613d0
        Pr = 6.13d0
        cp = 4180d0
        beta = 2.1d-4
    else
        ! Engine Oil
        rho = 870d0
        mu = 0.05d0
        kf = 0.14d0
        Pr = 500d0
        cp = 2000d0
        beta = 7.0d-4
    end if
    
    ! Override if user provided custom beta
    if (temp_val > 0.0d0) beta = temp_val
    
    nu_v = mu / rho
    alpha = kf / (rho * cp)
    As = Lp * Wp
    
    if (th < 89.9d0) then
        ! Vertical or Inclined Plate (theta < 90)
        ! Use vertical plate correlation but with effective gravity g*cos(theta)
        g_eff = g * cos(th * pi / 180.0d0)
        Lc = Lp
        Ra = g_eff * beta * dT * Lc**3 / (nu_v * alpha)
        
        ! Churchill and Chu (1975) Vertical Plate correlation
        Nu = (0.825d0 + (0.387d0 * Ra**(1.0d0/6.0d0)) / &
             ((1.0d0 + (0.492d0/Pr)**(9.0d0/16.0d0))**(8.0d0/27.0d0)))**2
    else
        ! Horizontal Plate (theta = 90)
        ! Characteristic length Lc = Area / Perimeter
        Lc = As / (2.0d0 * (Lp + Wp))
        Ra = g * beta * dT * Lc**3 / (nu_v * alpha)
        
        if (config_type == 1) then
            ! Upper Hot or Lower Cold (Upper surface of hot plate / lower surface of cold plate)
            if (Ra < 1.0d7) then
                Nu = 0.54d0 * Ra**0.25d0
            else
                Nu = 0.15d0 * Ra**(1.0d0/3.0d0)
            end if
        else
            ! Lower Hot or Upper Cold (Lower surface of hot plate / upper surface of cold plate)
            Nu = 0.27d0 * Ra**0.25d0
        end if
    end if
    
    h = Nu * kf / Lc
    Q = h * As * dT
    
    write(*,'(A)') '============================================'
    write(*,'(A)') '  INCLINED & HORIZONTAL PLATE CONVECTION'
    write(*,'(A)') '============================================'
    write(*,*)
    write(*,'(A,ES14.4)') '  Rayleigh Number Ra      = ', Ra
    write(*,'(A,F12.4)')    '  Nusselt Number Nu       = ', Nu
    write(*,'(A,F12.4,A)')  '  Coeff h                 = ', h, ' W/m2K'
    write(*,'(A,F12.4,A)')  '  Transfer Q              = ', Q, ' W'
    write(*,*)
    
    write(*,'(A)') '--- DELTA-T SWEEP ---'
    write(*,'(A)') '  dT[C]      Ra           Nu        h[W/m2K]   Q[W]'
    write(*,'(A)') '  ------------------------------------------------------------'
    do i=1,25
        dTs = 1.0d0 + (max(dT, 10.0d0)*3.0d0 - 1.0d0)*dble(i-1)/24.0d0
        if (th < 89.9d0) then
            Ras = g_eff * beta * dTs * Lc**3 / (nu_v * alpha)
            Nus = (0.825d0 + (0.387d0 * Ras**(1.0d0/6.0d0)) / &
                  ((1.0d0 + (0.492d0/Pr)**(9.0d0/16.0d0))**(8.0d0/27.0d0)))**2
        else
            Ras = g * beta * dTs * Lc**3 / (nu_v * alpha)
            if (config_type == 1) then
                if (Ras < 1.0d7) then
                    Nus = 0.54d0 * Ras**0.25d0
                else
                    Nus = 0.15d0 * Ras**(1.0d0/3.0d0)
                end if
            else
                Nus = 0.27d0 * Ras**0.25d0
            end if
        end if
        hs = Nus * kf / Lc
        Qs = hs * As * dTs
        write(*,'(2X,F8.2,2X,ES10.3,2X,F9.3,2X,F10.4,2X,F12.4)') dTs, Ras, Nus, hs, Qs
    end do
    write(*,*)
end program inclined_plate