S H I N E

Sunday, November 10, 2024

Friday, August 2, 2024

GSL integration

//https://stackoverflow.com/questions/47038457/gsl-integration-within-a-function-c

#include 
#include 
#include 
#include 

namespace details
{
  extern "C" {
  double gsl_function_cpp_wrapper_helper(double z, void *cpp_f)
  {
    const auto p_cpp_f =
        reinterpret_cast *>(cpp_f);
    return (*p_cpp_f)(z);
  }
  }
}

class gsl_function_cpp_wrapper
{
  using Cpp_F = std::function;

 public:
  operator const gsl_function_struct *() const { return &gsl_f_; };

  gsl_function_cpp_wrapper(Cpp_F &&cpp_f) : cpp_f_(std::move(cpp_f))
  {
    gsl_f_.function = details::gsl_function_cpp_wrapper_helper;
    gsl_f_.params = &cpp_f_;
  }

  gsl_function_cpp_wrapper(double(f)(double))
      : gsl_function_cpp_wrapper(Cpp_F(f))
  {
  }

  template 
  gsl_function_cpp_wrapper(double (OBJ::*method)(double) const, const OBJ *obj)
      : gsl_function_cpp_wrapper(
            Cpp_F(std::bind(method, obj, std::placeholders::_1)))
  {
  }
  template 
  gsl_function_cpp_wrapper(double (OBJ::*method)(double), OBJ *obj)
      : gsl_function_cpp_wrapper(
            Cpp_F(std::bind(method, obj, std::placeholders::_1)))
  {
  }

 protected:
  Cpp_F cpp_f_;
  gsl_function_struct gsl_f_;
};

//----------------

class Class_Example
{
 public:
  double f(double x) const { return std::log(alpha_ * x) / std::sqrt(x); }

 protected:
  double alpha_ = 1;
};

double free_function_example(double x)
{
  const double alpha = 1;
  return std::log(alpha * x) / std::sqrt(x);
}

//----------------

int main()
{
  double result, error;
  gsl_integration_workspace *const w = gsl_integration_workspace_alloc(1000);
  assert(w != nullptr);

  //----------------

  Class_Example class_example;

  gsl_function_cpp_wrapper wrapper_1(&Class_Example::f, &class_example);

  gsl_integration_qags(wrapper_1, 0, 1, 0, 1e-7, 1000, w, &result, &error);

  std::printf("result          = % .18f\n", result);
  std::printf("estimated error = % .18f\n", error);

  //----------------

  gsl_function_cpp_wrapper wrapper_2(free_function_example);

  gsl_integration_qags(wrapper_2, 0, 1, 0, 1e-7, 1000, w, &result, &error);

  std::printf("result          = % .18f\n", result);
  std::printf("estimated error = % .18f\n", error);

  //----------------

  gsl_integration_workspace_free(w);

  return EXIT_SUCCESS;
}

GSL interpolation

#include 
#include 
#include 


 double INTERPOLATE1D( double logx,const vector& logx0, const vector& logy0)
 {
     const gsl_interp_type *T = gsl_interp_cspline;//gsl_interp_linear;
     gsl_spline *s = gsl_spline_alloc(T,logx0.size());
     gsl_interp_accel *acc = gsl_interp_accel_alloc();
     gsl_spline_init(s,&logx0[0],&logy0[0],logx0.size());
 
     double logy = gsl_spline_eval(s,logx,acc);
     gsl_spline_free(s);
     gsl_interp_accel_free(acc);
     return pow(10,logy);
 }

Thursday, August 1, 2024

pass parameters in GSL integration

struct my_f_params {int a; int b;}; 
gsl_function F1; 
struct my_f_params alpha = {2,2};               
F1.function = &f1; 
F1.params = & alpha;

struct my_f_params * params = (struct my_f_params *)p;
int n = (params->a);
int m = (params->b);

----------

double params[] = { y, z };
F.params = params;

double y = ((double *)params)[0];
double z =  *((double *)p+1);