analytical.cpp

Go to the documentation of this file.

#include "methods/analytical.hpp"
 
#include <cmath>
#include <numbers> // C++20 for std::numbers::pi
 
double Analytical::solve(double t, double x, const PhysParams& phys)
{
    // If t=0, return initial condition (handling discontinuity at boundaries if needed)
    // But strictly speaking, at t=0, T=Tin everywhere inside.
    if (t <= 0.0)
        return phys.Tin;
 
    // Constants
    const double pi = std::numbers::pi;
    const double L = phys.L_cm;
    const double D = phys.D_cm2h;
 
    double sum = 0.0;
    const int max_k = 200; // Number of terms for convergence
 
    for (int k = 0; k < max_k; ++k)
    {
        double term_k = 2.0 * k + 1.0;
        double lambda = (term_k * pi) / L;
 
        double sine_term = std::sin(lambda * x);
        double exp_term = std::exp(-D * lambda * lambda * t);
 
        sum += (1.0 / term_k) * sine_term * exp_term;
    }
 
    // T(x,t) = Tsur + (Tin - Tsur) * (4/pi) * sum
    return phys.Tsur + (phys.Tin - phys.Tsur) * (4.0 / pi) * sum;
}
 
std::vector<double>
Analytical::get_profile(double t, const std::vector<double>& x_grid, const PhysParams& phys)
{
    std::vector<double> T(x_grid.size());
    for (size_t i = 0; i < x_grid.size(); ++i)
    {
        // Handle boundaries explicitly to avoid series artifacts (Gibbs phenomenon) at x=0, x=L
        if (i == 0 || i == x_grid.size() - 1)
        {
            T[i] = phys.Tsur;
        }
        else
        {
            T[i] = solve(t, x_grid[i], phys);
        }
    }
    return T;
}