program semi_infinite_solid
    implicit none
    
    ! Input variables
    integer :: bc_type
    real(8) :: T_i, bc_param1, bc_param2, k, alpha, x, t
    
    ! Calculation variables
    real(8) :: T_xt, T_0t, q_0t, Q_t, delta_p
    real(8) :: pi, arg, arg1, arg2, arg3
    real(8) :: T_s, q_s, h, T_inf
    real(8) :: x_probe, T_probe, arg_p, arg3_p
    integer :: i
    
    pi = 3.141592653589793d0
    
    ! Read inputs
    read(*,*) bc_type
    read(*,*) T_i
    read(*,*) bc_param1
    read(*,*) bc_param2
    read(*,*) k
    read(*,*) alpha
    read(*,*) x
    read(*,*) t
    
    ! Input safety checks
    if (alpha <= 0.0d0) alpha = 1.0d-7
    if (k <= 0.0d0) k = 1.0d0
    if (t <= 0.0d0) t = 1.0d-5
    if (x < 0.0d0) x = 0.0d0
    
    delta_p = 4.0d0 * sqrt(alpha * t)
    
    if (bc_type == 1) then
        ! Constant Surface Temperature
        T_s = bc_param1
        arg = x / (2.0d0 * sqrt(alpha * t))
        T_xt = T_s + (T_i - T_s) * erf(arg)
        T_0t = T_s
        q_0t = k * (T_s - T_i) / sqrt(pi * alpha * t)
        Q_t = 2.0d0 * k * (T_s - T_i) * sqrt(t / (pi * alpha))
        
    else if (bc_type == 2) then
        ! Constant Surface Heat Flux
        q_s = bc_param1
        arg = x / (2.0d0 * sqrt(alpha * t))
        T_xt = T_i + (2.0d0 * q_s * sqrt(alpha * t / pi) / k) * exp(-arg**2) - (q_s * x / k) * erfc(arg)
        T_0t = T_i + (2.0d0 * q_s / k) * sqrt(alpha * t / pi)
        q_0t = q_s
        Q_t = q_s * t
        
    else if (bc_type == 3) then
        ! Surface Convection
        h = bc_param1
        T_inf = bc_param2
        arg1 = x / (2.0d0 * sqrt(alpha * t))
        arg2 = h * sqrt(alpha * t) / k
        arg3 = h * x / k + (h**2 * alpha * t) / (k**2)
        T_xt = T_i + (T_inf - T_i) * (erfc(arg1) - exp(arg3) * erfc(arg1 + arg2))
        T_0t = T_i + (T_inf - T_i) * (1.0d0 - exp(arg2**2) * erfc(arg2))
        q_0t = h * (T_inf - T_0t)
        Q_t = (T_inf - T_i) * (k**2 / (h * alpha)) * (exp(arg2**2) * erfc(arg2) - 1.0d0 + 2.0d0 * arg2 / sqrt(pi))
        
    else
        ! Default fallback
        bc_type = 1
        T_s = T_i
        T_xt = T_i
        T_0t = T_i
        q_0t = 0.0d0
        Q_t = 0.0d0
    end if
    
    ! Print output report
    write(*,*) "=================================================="
    write(*,*) "    SEMI-INFINITE SOLID CONDUCTION REPORT"
    write(*,*) "=================================================="
    write(*,*) ""
    write(*,*) "THERMO-PHYSICAL PROPERTIES:"
    write(*,*) "--------------------------------------------------"
    write(*, '(A, F10.4, A)') "  Conductivity (k):        ", k, " W/m.K"
    write(*, '(A, E12.4, A)') "  Diffusivity (alpha):     ", alpha, " m2/s"
    write(*, '(A, F10.2, A)') "  Initial Temp (T_i):      ", T_i, " C"
    write(*,*) ""
    
    write(*,*) "BOUNDARY CONDITIONS:"
    write(*,*) "--------------------------------------------------"
    if (bc_type == 1) then
        write(*,*) "  Type: Constant Surface Temperature (T_s)"
        write(*, '(A, F10.2, A)') "  Surface Temp (T_s):      ", T_s, " C"
    else if (bc_type == 2) then
        write(*,*) "  Type: Constant Surface Heat Flux (q''s)"
        write(*, '(A, F10.2, A)') "  Heat Flux (q''s):        ", q_s, " W/m2"
    else if (bc_type == 3) then
        write(*,*) "  Type: Surface Convection (h, T_inf)"
        write(*, '(A, F10.2, A)') "  Convection Coeff (h):    ", h, " W/m2.K"
        write(*, '(A, F10.2, A)') "  Ambient Temp (T_inf):    ", T_inf, " C"
    end if
    write(*,*) ""
    
    write(*,*) "TRANSIENT SIMULATION INPUTS:"
    write(*,*) "--------------------------------------------------"
    write(*, '(A, F10.4, A)') "  Target Depth (x):        ", x, " m"
    write(*, '(A, F10.2, A)') "  Target Time (t):         ", t, " s"
    write(*,*) ""
    
    write(*,*) "CALCULATED TRANSIENT METRICS:"
    write(*,*) "--------------------------------------------------"
    write(*, '(A, F10.4, A)') "  Penetration Depth (delta):", delta_p, " m"
    write(*, '(A, F10.2, A)') "  Surface Temp T(0, t):     ", T_0t, " C"
    write(*, '(A, F14.2, A)') "  Surface Heat Flux q''(0,t):", q_0t, " W/m2"
    write(*, '(A, ES14.4, A)') "  Cumulative Heat Trans (Q):", Q_t, " J/m2"
    write(*, '(A, F10.2, A)') "  Temp at x,t T(x, t):      ", T_xt, " C"
    write(*,*) ""
    
    write(*,*) "TEMPERATURE PROFILE TABLE T(x,t) at t =", t, "s"
    write(*,*) "--------------------------------------------------"
    write(*,*) "   Depth (mm)      Depth (m)       Temp (C)"
    write(*,*) "--------------------------------------------------"
    
    do i = 0, 10
        x_probe = (dble(i) / 10.0d0) * max(x, delta_p)
        if (bc_type == 1) then
            arg_p = x_probe / (2.0d0 * sqrt(alpha * t))
            T_probe = T_s + (T_i - T_s) * erf(arg_p)
        else if (bc_type == 2) then
            arg_p = x_probe / (2.0d0 * sqrt(alpha * t))
            T_probe = T_i + (2.0d0 * q_s * sqrt(alpha * t / pi) / k) * exp(-arg_p**2) - (q_s * x_probe / k) * erfc(arg_p)
        else if (bc_type == 3) then
            arg_p = x_probe / (2.0d0 * sqrt(alpha * t))
            arg3_p = h * x_probe / k + (h**2 * alpha * t) / (k**2)
            T_probe = T_i + (T_inf - T_i) * (erfc(arg_p) - exp(arg3_p) * erfc(arg_p + arg2))
        end if
        write(*, '(f12.3, f15.4, f15.2)') x_probe * 1000.0d0, x_probe, T_probe
    end do
    write(*,*) "=================================================="
    
end program semi_infinite_solid