CylindricalBesselGenerator.cpp

#include <cassert>

#include "compat.h"

#ifndef NO_BESSEL_TABLE
#include <greens_functions/CylindricalBesselTable.hpp>
#endif

#include "CylindricalBesselGenerator.hpp"

namespace greens_functions{

#ifndef NO_BESSEL_TABLE
static inline double hermite_interp(double x, 
                                    double x0, double dx, 
                                    double const* y_array)
#else
static inline double hermite_interp(double x, 
                                    double x0, double dx, 
                                    std::vector<double> const& y_array)
#endif
{
    const double hinv = 1.0 / dx;

    const size_t i = static_cast<size_t>((x - x0 ) * hinv);
    const size_t index = i * 2;

    const double x_lo = (x - x0) * hinv - i;
    const double x_hi =  1.0 - x_lo;

    const double y_lo = y_array[index];
    const double ydot_lo = y_array[index + 1] * dx;
    const double y_hi = y_array[index + 2];
    const double ydot_hi = y_array[index + 3] * dx;
    
    return x_hi * x_hi * (y_lo + x_lo * (2 * y_lo + ydot_lo)) 
        + x_lo * x_lo * (y_hi + x_hi * (2 * y_hi - ydot_hi));
}

#ifndef NO_BESSEL_TABLE
inline static Real interp(Real x_start, Real delta_x,
                          Real const* yTable, Real x)
#else
inline static Real interp(Real x_start, Real delta_x,
                          std::vector<double> const& yTable, Real x)
#endif
{
    return hermite_interp(x, x_start, delta_x, yTable);
}

static Real _J(UnsignedInteger n, Real z)
{
    return gsl_sf_bessel_Jn(n, z);
}

static Real _Y(UnsignedInteger n, Real z)
{
    return gsl_sf_bessel_Yn(n, z);
}

CylindricalBesselGenerator const& CylindricalBesselGenerator::instance()
{
    static const CylindricalBesselGenerator cylindricalBesselGenerator;
    return cylindricalBesselGenerator;
}

UnsignedInteger CylindricalBesselGenerator::getMinNJ()
{
    return cb_table::cj_table_min;
}

UnsignedInteger CylindricalBesselGenerator::getMinNY()
{
    return cb_table::cy_table_min;
}

UnsignedInteger CylindricalBesselGenerator::getMaxNJ()
{
    return cb_table::cj_table_max;
}

UnsignedInteger CylindricalBesselGenerator::getMaxNY()
{
    return cb_table::cy_table_max;
}

#ifndef NO_BESSEL_TABLE
static cb_table::Table const* getCJTable(UnsignedInteger n)
{
    return cb_table::cj_table[n];
}


static cb_table::Table const* getCYTable(UnsignedInteger n)
{
    return cb_table::cy_table[n];
}

static inline Real _J_table(UnsignedInteger n, Real z)
{
    cb_table::Table const* tablen(getCJTable(n));

    return interp(tablen->x_start, tablen->delta_x, tablen->y, z);
}

static inline Real _Y_table(UnsignedInteger n, Real z)
{
    cb_table::Table const* tablen(getCYTable(n));

    return interp(tablen->x_start, tablen->delta_x, tablen->y, z);
}
#else
cb_table::Table const* CylindricalBesselGenerator::getCJTable(UnsignedInteger n) const
{
    return &cj_table_[n];
}


cb_table::Table const* CylindricalBesselGenerator::getCYTable(UnsignedInteger n) const
{
    return &cy_table_[n];
}

Real CylindricalBesselGenerator::_J_table(UnsignedInteger n, Real z) const
{
    cb_table::Table const* tablen(getCJTable(n));

    return interp(tablen->x_start, tablen->delta_x, tablen->y, z);
}

Real CylindricalBesselGenerator::_Y_table(UnsignedInteger n, Real z) const
{
    cb_table::Table const* tablen(getCYTable(n));

    return interp(tablen->x_start, tablen->delta_x, tablen->y, z);
}
#endif

Real CylindricalBesselGenerator::J(UnsignedInteger n, Real z) const
{
    if(n > getMaxNJ())
    {
        return _J(n, z);
    }
    
    const cb_table::Table* table(getCJTable(n));
    assert(table != 0);

    const Real minz(table->x_start + table->delta_x * 3);
    const Real maxz(table->x_start + table->delta_x * (table->N-3));
    
    if(z >= minz && z < maxz)
    {
        return _J_table(n, z);
    }
    else
    {
        return _J(n, z);
    }
}

Real CylindricalBesselGenerator::Y(const UnsignedInteger n, const Real z) const
{
    if(n > getMaxNY())
    {
        return _Y(n, z);
    }
    
    const cb_table::Table* table(getCYTable(n));
    assert(table != 0);
    
    const Real minz(table->x_start + table->delta_x * 3);
    const Real maxz(table->x_start + table->delta_x * (table->N-3));
    
    if(z >= minz && z < maxz)
    {
        return _Y_table(n, z);
    }
    else
    {
        return _Y(n, z);
    }
}

}