program transient_multidim_conduction_adi
    implicit none
    integer :: iostat_val, iter
    double precision :: Lx_mm, Ly_mm, k_cond, rho_dens, cp_spec, Ti_C, Tinf_C, h_conv, time_sec
    double precision :: Lx_m, Ly_m, alpha_diff, Bix, Biy, Fox, Foy
    double precision :: zeta_x, zeta_y, C_x, C_y, theta_x0, theta_y0, theta_x_surf, theta_y_surf
    double precision :: theta_center, theta_corner, T_center_C, T_corner_C, Q_ratio
    double precision :: z_low, z_high, z_mid, f_mid
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) Lx_mm      ! Half-width in X [mm] (e.g. 50.0)
    read(*,*,iostat=iostat_val) Ly_mm      ! Half-width in Y [mm] (e.g. 80.0)
    read(*,*,iostat=iostat_val) k_cond     ! Thermal Conductivity [W/(m.K)] (e.g. 45.0 for carbon steel)
    read(*,*,iostat=iostat_val) rho_dens   ! Density [kg/m3] (e.g. 7850.0)
    read(*,*,iostat=iostat_val) cp_spec    ! Specific Heat [J/(kg.K)] (e.g. 480.0)
    read(*,*,iostat=iostat_val) Ti_C       ! Initial Temperature [deg C] (e.g. 850.0)
    read(*,*,iostat=iostat_val) Tinf_C     ! Ambient Quench Temp [deg C] (e.g. 30.0)
    read(*,*,iostat=iostat_val) h_conv     ! Convection Coeff [W/(m2.K)] (e.g. 350.0)
    read(*,*,iostat=iostat_val) time_sec   ! Elapsed Time [s] (e.g. 300.0)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for 2D transient conduction calculation.'
        stop
    end if

    if (Lx_mm <= 0.0d0 .or. Ly_mm <= 0.0d0 .or. k_cond <= 0.0d0 .or. time_sec < 0.0d0) then
        write(*,*) 'ERROR: Dimensions, conductivity, and time must be positive.'
        stop
    end if

    Lx_m = Lx_mm * 1.0d-3
    Ly_m = Ly_mm * 1.0d-3
    alpha_diff = k_cond / (rho_dens * cp_spec) ! m2/s

    Bix = (h_conv * Lx_m) / k_cond
    Biy = (h_conv * Ly_m) / k_cond
    Fox = (alpha_diff * time_sec) / (Lx_m**2)
    Foy = (alpha_diff * time_sec) / (Ly_m**2)

    ! Solve transcendental root zeta_x * tan(zeta_x) = Bix
    z_low = 0.0001d0
    z_high = min(1.5707d0, max(0.1d0, atan(Bix)))
    if (Bix > 10.0d0) z_high = 1.5700d0
    do iter = 1, 60
        z_mid = 0.5d0 * (z_low + z_high)
        f_mid = z_mid * tan(z_mid) - Bix
        if (f_mid < 0.0d0) then
            z_low = z_mid
        else
            z_high = z_mid
        end if
    end do
    zeta_x = z_mid
    C_x = (4.0d0 * sin(zeta_x)) / (2.0d0 * zeta_x + sin(2.0d0 * zeta_x))

    ! Solve transcendental root zeta_y * tan(zeta_y) = Biy
    z_low = 0.0001d0
    z_high = min(1.5707d0, max(0.1d0, atan(Biy)))
    if (Biy > 10.0d0) z_high = 1.5700d0
    do iter = 1, 60
        z_mid = 0.5d0 * (z_low + z_high)
        f_mid = z_mid * tan(z_mid) - Biy
        if (f_mid < 0.0d0) then
            z_low = z_mid
        else
            z_high = z_mid
        end if
    end do
    zeta_y = z_mid
    C_y = (4.0d0 * sin(zeta_y)) / (2.0d0 * zeta_y + sin(2.0d0 * zeta_y))

    ! Product Solution for 2D Conduction
    theta_x0 = C_x * exp(-(zeta_x**2) * Fox)
    theta_y0 = C_y * exp(-(zeta_y**2) * Foy)
    if (theta_x0 > 1.0d0) theta_x0 = 1.0d0
    if (theta_y0 > 1.0d0) theta_y0 = 1.0d0

    theta_center = theta_x0 * theta_y0

    theta_x_surf = theta_x0 * cos(zeta_x)
    theta_y_surf = theta_y0 * cos(zeta_y)
    theta_corner = theta_x_surf * theta_y_surf

    T_center_C = Tinf_C + (Ti_C - Tinf_C) * theta_center
    T_corner_C = Tinf_C + (Ti_C - Tinf_C) * theta_corner

    ! Total energy removed ratio Q/Qmax
    Q_ratio = 1.0d0 - theta_center * (sin(zeta_x)*sin(zeta_y) / (zeta_x*zeta_y))

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — 2D TRANSIENT CONDUCTION ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.1,A)') 'Billet 2Lx x 2Ly          = ', 2.0d0*Lx_mm, ' mm x ', 2.0d0*Ly_mm, ' mm'
    write(*,'(A,F10.2,A,E12.4,A)')'Conductivity / Diffusivity= ', k_cond, ' W/m.K / ', alpha_diff, ' m2/s'
    write(*,'(A,F10.1,A,F10.1,A)') 'Initial / Ambient Temp    = ', Ti_C, ' deg C / ', Tinf_C, ' deg C'
    write(*,'(A,F10.1,A,F10.1,A)') 'Quench Time / Convection h= ', time_sec, ' s / ', h_conv, ' W/m2.K'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.4,A,F10.4)')   'Biot Numbers (Bix / Biy)  = ', Bix, ' / ', Biy
    write(*,'(A,F10.4,A,F10.4)')   'Fourier Numbers (Fox/Foy) = ', Fox, ' / ', Foy
    write(*,'(A,F10.2,A,F8.4)')    'CENTER TEMPERATURE (T0)   = ', T_center_C, ' deg C (theta = ', theta_center, ')'
    write(*,'(A,F10.2,A,F8.4)')    'CORNER TEMPERATURE (Tcorn)= ', T_corner_C, ' deg C (theta = ', theta_corner, ')'
    write(*,'(A,F10.2,A)')  'Energy Removed (Q/Qmax)   = ', Q_ratio * 100.0d0, ' %'
    write(*,'(A)') '============================================================'

end program transient_multidim_conduction_adi