! ==============================================================================
! Stefan Problem — Phase Change (Solidification/Melting) Calculator
! Neumann Analytical Solution for 1D Semi-Infinite Media
! References: Carslaw & Jaeger Ch.11, Alexiades & Solomon
! ==============================================================================
program stefan_problem
    implicit none
    
    ! Inputs
    integer :: mode ! 1=Solidification, 2=Melting
    real(8) :: Tm ! Phase change temp (C)
    real(8) :: L ! Latent heat of fusion (kJ/kg)
    real(8) :: Ti ! Initial temperature (C)
    real(8) :: Ts ! Wall/surface temperature (C)
    
    ! Solid properties
    real(8) :: ks ! conductivity (W/m-K)
    real(8) :: rhos ! density (kg/m3)
    real(8) :: Cps ! specific heat (J/kg-K)
    
    ! Liquid properties
    real(8) :: kl ! conductivity (W/m-K)
    real(8) :: rhol ! density (kg/m3)
    real(8) :: Cpl ! specific heat (J/kg-K)
    
    ! Dimensions & simulation time
    real(8) :: d_target ! target thickness for phase change (mm)
    real(8) :: x_eval ! evaluation depth (mm)
    real(8) :: t_sim ! elapsed time (s)
    
    ! Derived properties
    real(8) :: alphas ! solid diffusivity (m2/s)
    real(8) :: alphal ! liquid diffusivity (m2/s)
    real(8) :: Ste ! Stefan number
    real(8) :: lambda ! similarity growth constant
    
    ! Outputs
    real(8) :: front_pos ! front position s(t) (m)
    real(8) :: time_to_d ! solidification/melting time for d (s)
    real(8) :: temp_at_x ! temperature at x_eval (C)
    real(8) :: energy_transferred ! energy released/absorbed (kJ/m2)
    
    ! Constants & loop variables
    real(8) :: pi, eps
    integer :: iter
    real(8) :: a, b, c, fa, fc
    character(len=20) :: mode_str
    
    pi = 4.0d0 * datan(1.0d0)
    eps = 1.0d-12
    
    ! Read standard inputs
    read(*, *, iostat=iter) mode
    if (iter /= 0) then
        print *, "Error: Invalid input format (mode)."
        stop
    endif
    
    read(*, *) Tm
    read(*, *) L
    read(*, *) Ti
    read(*, *) Ts
    read(*, *) ks, rhos, Cps
    read(*, *) kl, rhol, Cpl
    read(*, *) d_target
    read(*, *) x_eval
    read(*, *) t_sim
    
    ! Sanity Checks & Input Validation
    if (L <= 0.0d0 .or. ks <= 0.0d0 .or. rhos <= 0.0d0 .or. Cps <= 0.0d0 &
        .or. kl <= 0.0d0 .or. rhol <= 0.0d0 .or. Cpl <= 0.0d0) then
        print *, "Error: Physical properties must be positive."
        stop
    endif
    
    if (d_target < 0.0d0 .or. x_eval < 0.0d0 .or. t_sim <= 0.0d0) then
        print *, "Error: Target thickness, depth, and time must be non-negative."
        stop
    endif
    
    if (mode == 1) then
        mode_str = "Solidification"
        if (Ts >= Tm) then
            print *, "Error: For Solidification, wall temperature (Ts) must be below melting temperature (Tm)."
            stop
        endif
        if (Ti < Tm) then
            print *, "Error: For Solidification, initial temperature (Ti) must be above or at melting temperature (Tm)."
            stop
        endif
    else if (mode == 2) then
        mode_str = "Melting"
        if (Ts <= Tm) then
            print *, "Error: For Melting, wall temperature (Ts) must be above melting temperature (Tm)."
            stop
        endif
        if (Ti > Tm) then
            print *, "Error: For Melting, initial temperature (Ti) must be below or at melting temperature (Tm)."
            stop
        endif
    else
        print *, "Error: Invalid Mode (must be 1 for Solidification or 2 for Melting)."
        stop
    endif
    
    ! Diffusivities
    alphas = ks / (rhos * Cps)
    alphal = kl / (rhol * Cpl)
    
    ! Stefan Number
    if (mode == 1) then
        Ste = Cps * (Tm - Ts) / (L * 1000.0d0)
    else
        Ste = Cpl * (Ts - Tm) / (L * 1000.0d0)
    endif
    
    ! Bisection Solver for Lambda
    a = 1.0d-8
    b = 10.0d0
    fa = evaluate_f(a)
    
    if (fa * evaluate_f(b) > 0.0d0) then
        print *, "Error: Root is not bracketed in [10^-8, 10.0]."
        stop
    endif
    
    do iter = 1, 100
        c = (a + b) / 2.0d0
        fc = evaluate_f(c)
        
        if (abs(fc) < eps .or. (b - a)/2.0d0 < eps) then
            lambda = c
            exit
        endif
        
        if (fa * fc < 0.0d0) then
            b = c
        else
            a = c
            fa = fc
        endif
    enddo
    
    ! Calculations
    if (mode == 1) then
        ! Growing phase is solid
        front_pos = 2.0d0 * lambda * sqrt(alphas * t_sim)
        time_to_d = (d_target / 1000.0d0)**2 / (4.0d0 * lambda**2 * alphas)
        
        ! Temperature at x_eval
        if (x_eval / 1000.0d0 < front_pos) then
            ! Solid region
            temp_at_x = Ts + (Tm - Ts) * erf( (x_eval / 1000.0d0) / (2.0d0 * sqrt(alphas * t_sim)) ) / erf(lambda)
        else
            ! Liquid region
            temp_at_x = Ti - (Ti - Tm) * erfc( (x_eval / 1000.0d0) / (2.0d0 * sqrt(alphal * t_sim)) ) / &
                        erfc(lambda * sqrt(alphas / alphal))
        endif
        
        ! Energy released (kJ/m2)
        energy_transferred = 2.0d0 * ks * (Tm - Ts) / erf(lambda) * sqrt(t_sim / (pi * alphas)) / 1000.0d0
    else
        ! Growing phase is liquid
        front_pos = 2.0d0 * lambda * sqrt(alphal * t_sim)
        time_to_d = (d_target / 1000.0d0)**2 / (4.0d0 * lambda**2 * alphal)
        
        ! Temperature at x_eval
        if (x_eval / 1000.0d0 < front_pos) then
            ! Liquid region
            temp_at_x = Ts - (Ts - Tm) * erf( (x_eval / 1000.0d0) / (2.0d0 * sqrt(alphal * t_sim)) ) / erf(lambda)
        else
            ! Solid region
            temp_at_x = Ti + (Tm - Ti) * erfc( (x_eval / 1000.0d0) / (2.0d0 * sqrt(alphas * t_sim)) ) / &
                        erfc(lambda * sqrt(alphal / alphas))
        endif
        
        ! Energy absorbed (kJ/m2)
        energy_transferred = 2.0d0 * kl * (Ts - Tm) / erf(lambda) * sqrt(t_sim / (pi * alphal)) / 1000.0d0
    endif
    
    ! ==========================================================================
    ! Print Results in a beautiful engineering report format
    ! ==========================================================================
    print *, "=========================================================================="
    print *, "             STEFAN PROBLEM PHASE CHANGE REPORT (NEUMANN SOLUTION)        "
    print *, "=========================================================================="
    print "(A, A20)", " Operation Mode:          ", trim(mode_str)
    print "(A, F10.2, A)", " Phase Change Temp (Tm):   ", Tm, " C"
    print "(A, F10.2, A)", " Latent Heat of Fusion (L):", L, " kJ/kg"
    print "(A, F10.2, A)", " Initial Temp (Ti):        ", Ti, " C"
    print "(A, F10.2, A)", " Wall/Surface Temp (Ts):   ", Ts, " C"
    print *
    print *, "--------------------------------------------------------------------------"
    print *, " MATERIAL PROPERTIES                                                      "
    print *, "--------------------------------------------------------------------------"
    print *, " Parameter               Solid Phase           Liquid Phase"
    print "(A, F18.3, F21.3)", " Conductivity (k): ", ks, kl
    print "(A, F18.1, F21.1)", " Density (rho):    ", rhos, rhol
    print "(A, F18.1, F21.1)", " Spec. Heat (Cp):  ", Cps, Cpl
    print "(A, ES18.4, ES21.4)", " Diffusivity (alpha):", alphas, alphal
    print *
    print *, "--------------------------------------------------------------------------"
    print *, " COMPUTATION & DIMENSIONLESS METRICS                                      "
    print *, "--------------------------------------------------------------------------"
    print "(A, F12.5)", " Stefan Number (Ste):           ", Ste
    print "(A, F12.6)", " Growth Eigenvalue (lambda):    ", lambda
    print *
    print *, "--------------------------------------------------------------------------"
    print *, " DETAILED PERFORMANCE RESULTS                                             "
    print *, "--------------------------------------------------------------------------"
    print "(A, F12.2, A)", " Elapsed Simulation Time:       ", t_sim, " s"
    print "(A, F12.2, A)", " Phase Interface Position s(t): ", front_pos * 1000.0d0, " mm"
    print "(A, F8.2, A, F12.2, A)", " Phase Change Time for d =", d_target, " mm: ", time_to_d, " s"
    print "(A, F8.2, A, F12.2, A)", " Temperature at depth x =", x_eval, " mm:   ", temp_at_x, " C"
    if (mode == 1) then
        print "(A, F12.2, A)", " Total Thermal Energy Released: ", energy_transferred, " kJ/m2"
    else
        print "(A, F12.2, A)", " Total Thermal Energy Absorbed: ", energy_transferred, " kJ/m2"
    endif
    print *, "=========================================================================="
    print *, " References:"
    print *, " 1. Carslaw, H. S. & Jaeger, J. C., Conduction of Heat in Solids, Ch. 11"
    print *, " 2. Alexiades, V. & Solomon, A. D., Mathematical Modeling of Phase Change"
    print *, "=========================================================================="
    
contains

    ! Helper function to evaluate the transcendental equation residual f(lambda)
    real(8) function evaluate_f(l)
        real(8), intent(in) :: l
        real(8) :: term1, term2, y
        
        ! Term 1: Ste / (exp(l^2) * erf(l))
        if (l > 20.0d0) then
            term1 = 0.0d0
        else
            term1 = Ste / (exp(l**2) * erf(l))
        endif
        
        ! Term 2: solid/liquid coupling term
        if (mode == 1) then
            ! Solidification
            y = l * sqrt(alphas / alphal)
            term2 = (kl * sqrt(alphas) * (Ti - Tm)) / &
                    (ks * sqrt(alphal) * (Tm - Ts)) * Ste * exp_erfc_reciprocal(y)
        else
            ! Melting
            y = l * sqrt(alphal / alphas)
            term2 = (ks * sqrt(alphal) * (Tm - Ti)) / &
                    (kl * sqrt(alphas) * (Ts - Tm)) * Ste * exp_erfc_reciprocal(y)
        endif
        
        evaluate_f = term1 - term2 - l * sqrt(pi)
    end function evaluate_f

    ! Helper to safely evaluate 1 / (exp(y^2) * erfc(y)) without overflow
    real(8) function exp_erfc_reciprocal(y)
        real(8), intent(in) :: y
        
        if (y < 10.0d0) then
            exp_erfc_reciprocal = 1.0d0 / (exp(y**2) * erfc(y))
        else
            ! Asymptotic expansion for y -> infinity
            exp_erfc_reciprocal = (y * sqrt(pi)) / &
                                  (1.0d0 - 0.5d0 / (y**2) + 0.75d0 / (y**4))
        endif
    end function exp_erfc_reciprocal

end program stefan_problem
