Installation

From sources

The sources for ibs can be downloaded from the Github repo.

You can either clone the public repository:

$ git clone git://github.com/tomerten/ibs

Detailed info about Installation

The generation of this package is based on cmake_cpp_pybind11_tutorial.

This IBS Library can be used as a independent C++ Library and as Python Library. Both are generated using CMake.

The C++ Library

We start by setting up the directory structure for the C++ Library.

$ mkdir cpp
$ cd cpp
$ mkdir include
$ mkdir src
$ cd include
$ mkdir ibs_bits

As a next step we add some code, the content of the first file will be explained later.

$ touch ibs
$ cd ibs_bits
$ touch twiss.hpp
$ cd ../../src
$ touch twiss.cpp

We now add a CMakeLists.txt file.

$ cd ..
$ touch CMakeLists.txt

The final directory structure should look like:

cpp/CMakeLists.txt
cpp/include/ibs
cpp/include/ibs_bits/twiss.hpp
cpp/src/twiss.cpp

The CMakeLists.txt file in the cpp folder will allow to build and install the independent C++ Library. Below we show the content of this file for convenience and later reference.

cmake_minimum_required(VERSION 3.10.2)

set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)

project(IBSLib VERSION 0.2.0)

# include dir
include_directories(/usr/local/include)

# Src
AUX_SOURCE_DIRECTORY(src SRC_FILES)

# Headers
set(PROJECT_SOURCE_DIR "src")
set(PROJECT_INCLUDE_DIR "include/ibs_bits")

# Source files
set(SOURCE_FILES
    ${PROJECT_INCLUDE_DIR}/twiss.hpp
    ${PROJECT_INCLUDE_DIR}/RadiationDamping.hpp
    ${PROJECT_INCLUDE_DIR}/NumericFunctions.hpp
    ${PROJECT_INCLUDE_DIR}/CoulombLogFunctions.hpp
    ${PROJECT_INCLUDE_DIR}/Integrators.hpp
    ${PROJECT_INCLUDE_DIR}/Models.hpp
    ${PROJECT_INCLUDE_DIR}/OrdDiffEq.hpp
    ${PROJECT_SOURCE_DIR}/twiss.cpp
    ${PROJECT_SOURCE_DIR}/RadiationDamping.cpp
    ${PROJECT_SOURCE_DIR}/NumericFunctions.cpp
    ${PROJECT_SOURCE_DIR}/CoulombLogFunctions.cpp
    ${PROJECT_SOURCE_DIR}/Integrators.cpp
    ${PROJECT_SOURCE_DIR}/Models.cpp
    ${PROJECT_SOURCE_DIR}/OrdDiffEq.cpp
)

#file (GLOB SOURCE_FILES "${PROJECT_INCLUDE_DIR}/*.hpp" "${PROJECT_SOURCE_DIR}/*.cpp")

# Set up such that XCode organizes the files correctly
source_group(TREE ${CMAKE_CURRENT_SOURCE_DIR} FILES ${SOURCE_FILES})

# Add Library
add_library(${PROJECT_NAME}  SHARED ${SOURCE_FILES})

# Include directories
target_include_directories(${PROJECT_NAME} PRIVATE include/)

# Install
install(TARGETS ${PROJECT_NAME} DESTINATION lib)

# Install the headers
install(FILES include/ibs DESTINATION include)

# create base directory
install(DIRECTORY include/ DESTINATION include)

We also need to have the main header file find all the actual libraries (we use the cpp/include/ibs file).

#ifndef IBS_LIBRARY_H
#define IBS_LIBRARY_H

#include "ibs_bits/twiss.hpp"
#include "ibs_bits/NumericFunctions.hpp"
#include "ibs_bits/RadiationDamping.hpp"
#include "ibs_bits/CoulombLogFunctions.hpp"
#include "ibs_bits/Integrators.hpp"
#include "ibs_bits/Models.hpp"
#include "ibs_bits/OrdDiffEq.hpp"

#endif

Building the C++ Library

$ cd cpp
$ mkdir build
$ cd build
$ cmake ..
$ make
$ make install

For generating an XCode project IBSLib.xcodeproject in the build director:

$ cd cpp/build
$ cmake .. -GXcode

Testing the C++ Library

Note

The tests are not actual proper unit tests, but rather a place where some test code can be put to test the library functions.

We again start by creating the directory and adding files to them:

$ cd cpp
$ mkdir tests
$ cd tests
$ touch CMakeLists.txt
$ mkdir src
$ cd src
$ touch test_cpp.cpp

The specific content for the a quick test of the Twiss library is shown below.

#include <ibs>
#include <iostream>
#include <map>
#include <stdio.h>
#include <string>
#include <vector>

int main() {
  string twissfilename = "../src/b2_design_lattice_1996.twiss";
  map<string, double> twissheadermap;
  twissheadermap = GetTwissHeader(twissfilename);

  printf("Test %12.6e\n", twissheadermap["GAMMA"]);
  return 0;
}

The content of the corresponding CMakeLists.txt file:

cmake_minimum_required(VERSION 3.10.2)

set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)


project(test_cpp)

include_directories(/usr/local/include)

set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/../bin)
set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_RELEASE ${CMAKE_BINARY_DIR}/../bin)
set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_DEBUG ${CMAKE_BINARY_DIR}/../bin)

find_library(IBSLIB_LIB IBSLib HINTS /usr/local/lib/)

add_executable(test_cpp src/Demo_twiss.cpp)
add_executable(test_radiation_damping_cpp src/DemoRadDamping.cpp)
add_executable(test_coulomblog_functions_cpp src/DemoCoulombLog.cpp)
add_executable(test_integrators_cpp src/DemoIntegrators.cpp)
add_executable(test_ibs_models_cpp src/DemoIBS.cpp)
add_executable(test_ibs_ode_cpp src/DemoODE.cpp)


target_link_libraries(test_cpp PUBLIC ${IBSLIB_LIB})
target_link_libraries(test_radiation_damping_cpp PUBLIC ${IBSLIB_LIB})
target_link_libraries(test_coulomblog_functions_cpp PUBLIC ${IBSLIB_LIB})
target_link_libraries(test_integrators_cpp PUBLIC ${IBSLIB_LIB})
target_link_libraries(test_ibs_models_cpp PUBLIC ${IBSLIB_LIB})
target_link_libraries(test_ibs_ode_cpp PUBLIC ${IBSLIB_LIB})

Now we are prepared to build and run the tests.

$ cd cpp/tests
$ mkdir build
$ cd build
$ cmake ..
$ cd ../bin
4 ./test_cpp

Python wrapper

PYbind11 and git

For the Python wrapper we rely on pybind11. There are several methods to this package (e.g. using conda install) but here we decided to add it in third party folder using git submodules. Therefor to get it one needs to run:

$ git init
$ git add .
$ git commit -am "initial commit"
$ git submodule add https://github.com/pybind/pybind11.git third_party/pybind11-2.6.2
$ git add .

Note

Note that the order of commits and immediately after adding the submodule doing a new commit is important to not get errors.

As git submodules are not updated automatically we need to add some extra code to the CMakeLists.txt file:

PROJECT(IBSLib VERSION 0.2.0)

# DOWNLOAD THE GIT SUBMODULES (PYBIND11 HERE)
find_package(Git QUIET)

if (GIT_FOUND AND EXISTS "${PROJECT_SOURCE_DIR}/.git")
    # update as needed
    option(GIT_SUBMODULE "Check submodules during build" ON)

    if(GIT_SUBMODULE)
        message(STATUS "Updating submodules...")
        execute_process(COMMAND ${GIT_EXECUTABLE} submodule update --init --recursive
            WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
            RESULT_VARIABLE GIT_SUBMOD_RESULT)
        if(NOT GIT_SUBMOD_RESULT EQUAL "0")
            message(FATAL_ERROR "git submodule update --init failed with ${GIT_SUBMOD_RESULT}, please check...")
        endif()
    endif()
endif()

# CHECK ALL THE SUBMODULES

Wrapping the code

In order to wrap the C++ code we need to add the following extra code into the Python directory.

$ touch python/pytwiss.cpp

The content of this file is given by:

//#include "../cpp/include/ibs_bits/CoulombLogFunctions.hpp"
//#include "../cpp/include/ibs_bits/NumericFunctions.hpp"
//#include "../cpp/include/ibs_bits/RadiationDamping.hpp"
//#include "../cpp/include/ibs_bits/twiss.hpp"
#include <ibs>
#include <pybind11/functional.h>
#include <pybind11/numpy.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>

// PYBIND11_MAKE_OPAQUE(std::vector<double>);
namespace py = pybind11;

PYBIND11_MODULE(IBSLib, m) {
  m.doc() = "Python wrapper around C++ IBSLib.";
  /*
================================================================================
                                    TWISS
================================================================================
  */
  m.def("GetTwissHeader", &GetTwissHeader,
        "Get the twiss header as a dictionary.", py::arg("filename"));

  m.def("GetTwissTable", &GetTwissTableAsMap, "Get the twiss data table.",
        py::arg("filename"));

  m.def("updateTwiss",
        [](map<string, vector<double>> &table) {
          updateTwiss(table);
          return table;
        },
        "Extend Twiss Table with rad int, CS gamma, curly H and rho.",
        py::arg("table"));

  m.def("printTwissColumn", &printTwissMap, "Print Twiss column",
        py::arg("columnName"), py::arg("twissTableMap"));
  /*
================================================================================
                                    CONSTANTS
================================================================================
  */
  m.attr("clight") = py::float_(clight);
  m.attr("hbarGeV") = py::float_(hbar);
  m.attr("electron_mass") = py::float_(emass);
  m.attr("proton_mass") = py::float_(pmass);
  m.attr("neutron_mass") = py::float_(nmass);
  m.attr("mu_mass") = py::float_(mumass);
  m.attr("atomic_mass_unit") = py::float_(atomicmassunit);
  m.attr("pi") = py::float_(pi);
  m.attr("electric_charge") = py::float_(ec);
  m.attr("euler") = py::float_(euler);
  m.attr("electron_radius") = py::float_(erad);
  m.attr("proton_radius") = py::float_(prad);
  /*
================================================================================
                             RADIATION DAMPING
================================================================================
 */
  m.def("RadiationDampingApprox",
        [](double latticeLength, double gamma, double gammaTransition,
           double dipoleBendingRadius, double betax, double betay,
           py::array_t<double> out) {
          double *radint;
          radint = RadiationDampingApprox(latticeLength, gamma, gammaTransition,
                                          dipoleBendingRadius, betax, betay);
          // convert double pointer to numpy array
          auto buf_out = out.request();
          double *ptr_out = static_cast<double *>(buf_out.ptr);

          ptr_out[0] = radint[0];
          ptr_out[1] = radint[1];
          ptr_out[2] = radint[2];
          ptr_out[3] = radint[3];
          ptr_out[4] = radint[4];
          ptr_out[5] = radint[5];
          ptr_out[6] = radint[6];
        },
        "Radiation Damping using Ring averages.", py::arg("latticelength"),
        py::arg("gamma"), py::arg("gammaTransition"), py::arg("bendingradius"),
        py::arg("betax"), py::arg("betay"), py::arg("outputArray"));

  m.def("RadiationDampingLattice",
        [](map<string, vector<double>> &table, py::array_t<double> out) {
          double *radint;
          radint = RadiationDampingLattice(table);

          auto buf_out = out.request();
          double *ptr_out = static_cast<double *>(buf_out.ptr);
          ptr_out[0] = radint[0];
          ptr_out[1] = radint[1];
          ptr_out[2] = radint[2];
          ptr_out[3] = radint[3];
          ptr_out[4] = radint[4];
          ptr_out[5] = radint[5];
          ptr_out[6] = radint[6];
        },
        "Radiation damping weighted per element.", py::arg("twissTableMap"),
        py::arg("outputArray"));

  m.def(
      "RadiationDampingEquilibria",
      [](map<string, double> &twissheadermap, py::array_t<double> radint,
         double aatom, double qs, py::array_t<double> out) {
        double *equi;
        double r[7];
        auto buf = radint.request();
        double *ptr = static_cast<double *>(buf.ptr);
        for (int i = 0; i < 7; i++) {
          r[i] = ptr[i];
        }
        equi =
            RadiationDampingLifeTimesAndEquilibriumEmittancesWithPartitionNumbers(
                twissheadermap, r, aatom, qs);
        auto buf_out = out.request();
        double *ptr_out = static_cast<double *>(buf_out.ptr);
        ptr_out[0] = equi[0];
        ptr_out[1] = equi[1];
        ptr_out[2] = equi[2];
        ptr_out[3] = equi[3];
        ptr_out[4] = equi[4];
        ptr_out[5] = equi[5];
        ptr_out[6] = equi[6];
        ptr_out[7] = equi[7];
        ptr_out[8] = equi[8];
      },
      "Get radiation damping equilibs.", py::arg("twissHeaderMap"),
      py::arg("RadIntArray"), py::arg("AtomicMassNumber"),
      py::arg("synchrotronTune"), py::arg("outputArray"));

  m.def("get_energy_loss_per_turn", &RadiationLossesPerTurn,
        "Get the energy loss per turn.", py::arg("twissTableMap"),
        py::arg("I2"), py::arg("AtomicMassNumber"));

  /*
 ================================================================================
                              NUMERIC FUNCTIONS
 ================================================================================
  */
  m.def("beta_relativistic_from_gamma", &BetaRelativisticFromGamma,
        "Beta relativistic.", py::arg("Gamma"));

  m.def("eta", &eta, "Slip factor", py::arg("Gamma"),
        py::arg("GammaTransition"));

  m.def("particle_radius", &ParticleRadius,
        "Particle radius from charge and atomic mass.", py::arg("charge"),
        py::arg("AtomicNumber"));

  m.def("dee_to_dpp", &dee_to_dpp, "DE/E to DP/P", py::arg("dee"),
        py::arg("beta"));

  m.def("dpp_to_dee", &dpp_to_dee, "DP/P to DE/E", py::arg("dpp"),
        py::arg("beta"));

  /*
 ================================================================================
                              LONGITUDINAL / RF
 ================================================================================
  */
  m.def("sige_from_sigs", &sigefromsigs, "Energy spread from bunch length.",
        py::arg("omega0"), py::arg("sigs"), py::arg("qs"), py::arg("gamma"),
        py::arg("gammatr"));

  m.def("sigs_from_sige", &sigsfromsige, "Bunch length from energy spread.",
        py::arg("sige"), py::arg("gamma"), py::arg("gammatr"),
        py::arg("omegas"));

  m.def("rf_voltage_in_ev",
        [](double phi, double c, std::vector<double> h, std::vector<double> v) {
          return EffectiveRFVoltageInElectronVolt(phi, c, h.size(), h.data(),
                                                  v.data());
        },
        "RF energy gain per turn [eV]", py::arg("phi"), py::arg("charge"),
        py::arg("harmonicNumbers"), py::arg("voltages"));

  m.def("rf_voltage_in_ev_prime",
        [](double phi, double c, std::vector<double> h, std::vector<double> v) {
          return EffectiveRFVoltageInElectronVoltPrime(phi, c, h.size(),
                                                       h.data(), v.data());
        },
        "Derivative of RF energy gain per turn [eV]", py::arg("phi"),
        py::arg("charge"), py::arg("harmonicNumbers"), py::arg("voltages"));

  m.def("rf_voltage_in_ev_with_rad_losses",
        [](double phi, double U0, double c, std::vector<double> h,
           std::vector<double> v) {
          return VeffRFeVRadlosses(phi, c, U0, h.size(), h.data(), v.data());
        },
        "RF energy gain per turn [eV] with Radiation Losses", py::arg("phi"),
        py::arg("U0"), py::arg("charge"), py::arg("harmonicNumbers"),
        py::arg("voltages"));

  m.def("get_synchronuous_phase",
        [](double target, double init_phi, double U0, double c,
           std::vector<double> h, std::vector<double> v, double e) {
          return SynchronuousPhase(target, init_phi, U0, c, h.size(), h.data(),
                                   v.data(), e);
        },
        "Synchronuous Phase", py::arg("targetEnergyGain"), py::arg("init_phi"),
        py::arg("U0"), py::arg("charge"), py::arg("harmonicNumbers"),
        py::arg("voltages"), py::arg("epsilon"));

  m.def("rf_voltage_with_potential_well_distortion",
        [](double target, double U0, double c, std::vector<double> h,
           std::vector<double> v, double L, double N, double sigs, double pc) {
          return VeffRFeVPotentialWellDistortion(
              target, U0, c, h.size(), h.data(), v.data(), L, N, sigs, pc);
        },
        "RF energy gain [eV] with Potential Well Distortion", py::arg("phi"),
        py::arg("U0"), py::arg("charge"), py::arg("harmonicNumbers"),
        py::arg("voltages"), py::arg("L"), py::arg("N"), py::arg("sigs"),
        py::arg("pc"));

  m.def("rf_voltage_with_potential_well_distortion_prime",
        [](double target, double U0, double c, std::vector<double> h,
           std::vector<double> v, double L, double N, double sigs, double pc) {
          return VeffRFeVPotentialWellDistortionPrime(
              target, U0, c, h.size(), h.data(), v.data(), L, N, sigs, pc);
        },
        "Derivative of RF energy gain [eV] with Potential Well Distortion",
        py::arg("phi"), py::arg("U0"), py::arg("charge"),
        py::arg("harmonicNumbers"), py::arg("voltages"), py::arg("L"),
        py::arg("N"), py::arg("sigs"), py::arg("pc"));

  m.def("get_synchronuous_phase_with_potential_well_distortion",
        [](double target, double init_phi, double U0, double c,
           std::vector<double> h, std::vector<double> v, double L, double N,
           double sigs, double pc, double e) {
          return SynchronuousPhaseWithPWD(target, init_phi, U0, c, h.size(),
                                          h.data(), v.data(), L, N, sigs, pc,
                                          e);
        },
        "Synchronuous Phase with Potential Well Distortion",
        py::arg("targetEnergyGain"), py::arg("init_phi"), py::arg("U0"),
        py::arg("charge"), py::arg("harmonicNumbers"), py::arg("voltages"),
        py::arg("L"), py::arg("N"), py::arg("sigs"), py::arg("pc"),
        py::arg("epsilon"));

  m.def("get_synchrotron_tune",
        [](double omega0, double U0, double c, std::vector<double> h,
           std::vector<double> v, double phis, double eta, double pc) {
          return SynchrotronTune(omega0, U0, c, h.size(), h.data(), v.data(),
                                 phis, eta, pc);
        },
        "Synchronuous Tune", py::arg("omega0"), py::arg("U0"),
        py::arg("charge"), py::arg("harmonicNumbers"), py::arg("voltages"),
        py::arg("phis"), py::arg("eta"), py::arg("pc"));

  m.def("get_synchrotron_tune_with_potential_well",
        [](double omega0, double U0, double c, std::vector<double> h,
           std::vector<double> v, double L, double N, double sigs, double phis,
           double eta, double pc) {
          return SynchrotronTunePWD(omega0, U0, c, h.size(), h.data(), v.data(),
                                    L, N, sigs, phis, eta, pc);
        },
        "Synchronuous Tune ith Potential Well Distortion", py::arg("omega0"),
        py::arg("U0"), py::arg("charge"), py::arg("harmonicNumbers"),
        py::arg("voltages"), py::arg("L"), py::arg("N"), py::arg("sigs"),
        py::arg("phis"), py::arg("eta"), py::arg("pc"));

  m.def("sige_form_sigs_using_rf",
        [](double sigs, double U0, double c, std::vector<double> h,
           std::vector<double> v, double gamma, double gammatr, double pc,
           double circ, double phis, bool printout) {
          return SigeFromRFAndSigs(sigs, U0, c, h.size(), h.data(), v.data(),
                                   gamma, gammatr, pc, circ, phis, printout);
        },
        "Sigma E from Sigma s using RF settings.", py::arg(""), py::arg("U0"),
        py::arg("charge"), py::arg("harmonicNumbers"), py::arg("voltages"),
        py::arg("gamma"), py::arg("gammaTransition"), py::arg("pc"),
        py::arg("circ"), py::arg("phis"), py::arg("printout"));

  /*
================================================================================
                       IBS NUMERIC FUNCTIONS
================================================================================
*/
  m.def("fmohl", &fmohl, "Fmohl function");

  m.def("rds", &rds, "Nagaitsev paper rds function");

  /*
================================================================================
                       COULOMB LOG
================================================================================
*/
  m.def("twclog",
        [](double pnumber, double bx, double by, double dx, double dy,
           double ex, double ey, double r0, double gamma, double charge,
           double en0, double amass, double sige, double sigt,
           py::array_t<double> npclog) {
          double clog[2];
          twclog(pnumber, bx, by, dx, dy, ex, ey, r0, gamma, charge, en0, amass,
                 sige, sigt, clog);
          auto buf_clog = npclog.request();
          double *ptr_clog = static_cast<double *>(buf_clog.ptr);
          ptr_clog[0] = clog[0];
          ptr_clog[1] = clog[1];
        },
        "Calculate the Coulomb Log.", py::arg("pnumber"), py::arg("betax"),
        py::arg("betay"), py::arg("dispersionx"), py::arg("disperisony"),
        py::arg("emitx"), py::arg("emity"), py::arg("classicalRadius"),
        py::arg("gamma"), py::arg("charge"), py::arg("energy"),
        py::arg("restEnergy_GeV"), py::arg("sige"), py::arg("sigt"),
        py::arg("outputArray"));

  m.def("twclogtail",
        [](double pnumber, double l, double bx, double by, double dx,
           double dpx, double dy, double dpy, double ax, double ay,
           double angle, double k1l, double k1sl, double ex, double ey,
           double r0, double aatom, double gamma, double en0, double len,
           double amass, double charge, double sige, double sigt,
           py::array_t<double> npclog) {
          double clog[2];
          twclogtail(pnumber, l, bx, by, dx, dpx, dy, dpy, ax, ay, angle, k1l,
                     k1sl, ex, ey, r0, aatom, gamma, en0, len, amass, charge,
                     sige, sigt, clog);
          auto buf_clog = npclog.request();
          double *ptr_clog = static_cast<double *>(buf_clog.ptr);
          ptr_clog[0] = clog[0];
          ptr_clog[1] = clog[1];
        },
        "Calculate the Coulomb Log with Tailcut.", py::arg("pnumber"),
        py::arg("elementLength"), py::arg("betax"), py::arg("betay"),
        py::arg("dispersionx"), py::arg("dispersionx_der"),
        py::arg("dispersiony"), py::arg("dispersiony_der"), py::arg("alphax"),
        py::arg("alphay"), py::arg("angle"), py::arg("k1l"), py::arg("k1sl"),
        py::arg("emitx"), py::arg("emity"), py::arg("classicalRadius"),
        py::arg("AtomicMassNumber"), py::arg("gamma"), py::arg("energy"),
        py::arg("acceleratorLength"), py::arg("restEnergy_GeV"),
        py::arg("charge"), py::arg("sige"), py::arg("sigt"),
        py::arg("outputArray"));

  m.def("coulomblog",
        [](double pnumber, double ex, double ey,
           map<string, double> &twissheader, double sige, double sigt,
           double r0, bool printout, py::array_t<double> npclog) {
          double clog[2];
          CoulombLog(pnumber, ex, ey, twissheader, sige, sigt, r0, printout,
                     clog);
          auto buf_clog = npclog.request();
          double *ptr_clog = static_cast<double *>(buf_clog.ptr);
          ptr_clog[0] = clog[0];
          ptr_clog[1] = clog[1];
        },
        "Calculate the Coulomb Log using ring average.", py::arg("pnumber"),
        py::arg("emitx"), py::arg("emity"), py::arg("twissheader"),
        py::arg("sige"), py::arg("sigt"), py::arg("classialRadius"),
        py::arg("printout"), py::arg("outputArray"));

  m.def("coulomblogtail",
        [](double pnumber, double ex, double ey,
           map<string, double> &twissheader, double sige, double sigt,
           double tauradx, double taurady, double taurads, double r0,
           bool printout, py::array_t<double> npclog) {
          double clog[2];
          TailCutCoulombLog(pnumber, ex, ey, twissheader, sige, sigt, tauradx,
                            taurady, taurads, r0, printout, clog);
          auto buf_clog = npclog.request();
          double *ptr_clog = static_cast<double *>(buf_clog.ptr);
          ptr_clog[0] = clog[0];
          ptr_clog[1] = clog[1];
        },
        "Calculate the Coulomb Log with tailcut using ring average.",
        py::arg("pnumber"), py::arg("emitx"), py::arg("emity"),
        py::arg("twissheader"), py::arg("sige"), py::arg("sigt"),
        py::arg("tauradx"), py::arg("taurady"), py::arg("taurads"),
        py::arg("classialRadius"), py::arg("printout"), py::arg("outputArray"));
  /*
================================================================================
                    INTEGRATORS
================================================================================
*/
  m.def("integrator_simpson_decade",
        [](double a, double b, double c, double cl, double cx, double cy,
           double cprime, double cyy, double tl1, double tl2, double tx1,
           double tx2, double ty1, double ty2, py::array_t<double> tau) {
          double alpha[3];
          SimpsonDecade(a, b, c, cl, cx, cy, cprime, cyy, tl1, tl2, tx1, tx2,
                        ty1, ty2, alpha);
          auto buf_tau = tau.request();
          double *ptr_tau = static_cast<double *>(buf_tau.ptr);
          ptr_tau[0] = alpha[0];
          ptr_tau[1] = alpha[1];
          ptr_tau[2] = alpha[2];
        },
        "Simpson integrator with decade splitting.", py::arg("a"), py::arg("b"),
        py::arg("c"), py::arg("cl"), py::arg("cx"), py::arg("cy"),
        py::arg("cprime"), py::arg("cyy"), py::arg("tl1"), py::arg("tl2"),
        py::arg("tx1"), py::arg("tx2"), py::arg("ty1"), py::arg("ty2"),
        py::arg("tau"));
  m.def("integrator_twsint",
        [](double pnumber, double ex, double ey, double sigs, double sige,
           double gammas, double betax, double betay, double alx, double aly,
           double dx, double dpx, double dy, double dpy,
           py::array_t<double> tau) {
          double alpha[3];
          twsint(pnumber, ex, ey, sigs, sige, gammas, betax, betay, alx, aly,
                 dx, dpx, dy, dpy, alpha);
          auto buf_tau = tau.request();
          double *ptr_tau = static_cast<double *>(buf_tau.ptr);
          ptr_tau[0] = alpha[0];
          ptr_tau[1] = alpha[1];
          ptr_tau[2] = alpha[2];
        },
        "MADX Simpson Decade with scaling integrator", py::arg("pnumber"),
        py::arg("emitx"), py::arg("emity"), py::arg("sigs"), py::arg("sige"),
        py::arg("gamma"), py::arg("betax"), py::arg("betay"), py::arg("alphax"),
        py::arg("alphay"), py::arg("dispersionx"), py::arg("dispersionx_der"),
        py::arg("dispersiony"), py::arg("dispersiony_der"), py::arg("tau"));
  m.def("integrand", &IBSIntegralIntegrand, "IBS integral integrand ",
        py::arg("int_var"), py::arg("ax"), py::arg("bx"), py::arg("a"),
        py::arg("b"), py::arg("c"));
  m.def("integrator_simpson", &simpson,
        "standard simpson integration on variable integrand (function)"); /*
[](const std::function<double(double, double, double, double, double,
        double)> &ibsintegrand,
double ax, double bx, double a, double b, double c, double al,
double bl,
int n) { return simpson(ibsintegrand, ax, bx, a, b, c, al, bl, n);
}, "Standard Simpson integrator");*/
  m.def("integrator_simpson_bjorken_mtingwa",
        [](const std::function<double(double, double, double, double, double,
                                      double)> &ibsintegrand,
           double ax, double bx, double ay, double by, double as, double bs,
           double a, double b, double ci, py::array_t<double> tau) {
          double alpha[3];
          intSimpson(ibsintegrand, ax, bx, ay, by, as, bs, a, b, ci, alpha);
          auto buf_tau = tau.request();
          double *ptr_tau = static_cast<double *>(buf_tau.ptr);
          ptr_tau[0] = alpha[0];
          ptr_tau[1] = alpha[1];
          ptr_tau[2] = alpha[2];
        },
        "Standard Simpson integral on Bjorken-Mtingwa integrand.");
  m.def("integral_bjorken_mtingwa",
        [](double pnumber, double ex, double ey, double sigs, double sige,
           double gammas, double betx, double bety, double alx, double aly,
           double dx, double dpx, double dy, double dpy,
           py::array_t<double> tau) {
          double alpha[3];
          BjorkenMtingwaInt(pnumber, ex, ey, sigs, sige, gammas, betx, bety,
                            alx, aly, dx, dpx, dy, dpy, alpha);
          auto buf_tau = tau.request();
          double *ptr_tau = static_cast<double *>(buf_tau.ptr);
          ptr_tau[0] = alpha[0];
          ptr_tau[1] = alpha[1];
          ptr_tau[2] = alpha[2];
        },
        "Bjorken-Mtingwa integral calculated using Simpson Decade",
        py::arg("pnumber"), py::arg("emitx"), py::arg("emity"), py::arg("sigs"),
        py::arg("sige"), py::arg("gamma"), py::arg("betax"), py::arg("betay"),
        py::arg("alphax"), py::arg("alphay"), py::arg("dispersionx"),
        py::arg("dispersionx_der"), py::arg("dispersiony"),
        py::arg("dispersiony_der"), py::arg("tau"));
  m.def("integral_conte_martini",
        [](double pnumber, double ex, double ey, double sigs, double sige,
           double gammas, double betx, double bety, double alx, double aly,
           double dx, double dpx, double dy, double dpy,
           py::array_t<double> tau) {
          double alpha[3];
          ConteMartiniInt(pnumber, ex, ey, sigs, sige, gammas, betx, bety, alx,
                          aly, dx, dpx, dy, dpy, alpha);
          auto buf_tau = tau.request();
          double *ptr_tau = static_cast<double *>(buf_tau.ptr);
          ptr_tau[0] = alpha[0];
          ptr_tau[1] = alpha[1];
          ptr_tau[2] = alpha[2];
        },
        "Conte-Martini integral calculated using Simpson Decade",
        py::arg("pnumber"), py::arg("emitx"), py::arg("emity"), py::arg("sigs"),
        py::arg("sige"), py::arg("gamma"), py::arg("betax"), py::arg("betay"),
        py::arg("alphax"), py::arg("alphay"), py::arg("dispersionx"),
        py::arg("dispersionx_der"), py::arg("dispersiony"),
        py::arg("dispersiony_der"), py::arg("tau"));
  m.def("integral_zimmerman",
        [](double pnumber, double ex, double ey, double sigs, double sige,
           double gammas, double betx, double bety, double alx, double aly,
           double dx, double dpx, double dy, double dpy,
           py::array_t<double> tau) {
          double alpha[3];
          MadxInt(pnumber, ex, ey, sigs, sige, gammas, betx, bety, alx, aly, dx,
                  dpx, dy, dpy, alpha);
          auto buf_tau = tau.request();
          double *ptr_tau = static_cast<double *>(buf_tau.ptr);
          ptr_tau[0] = alpha[0];
          ptr_tau[1] = alpha[1];
          ptr_tau[2] = alpha[2];
        },
        "Zimmerman integral calculated using Simpson Decade",
        py::arg("pnumber"), py::arg("emitx"), py::arg("emity"), py::arg("sigs"),
        py::arg("sige"), py::arg("gamma"), py::arg("betax"), py::arg("betay"),
        py::arg("alphax"), py::arg("alphay"), py::arg("dispersionx"),
        py::arg("dispersionx_der"), py::arg("dispersiony"),
        py::arg("dispersiony_der"), py::arg("tau"));
  /*
================================================================================
                 IBS MODELS
================================================================================
*/
  m.def("PiwinskiSmooth",
        [](double pnumber, double ex, double ey, double sigs, double dponp,
           map<string, double> &header, double r0, py::array_t<double> out) {
          double *ibs;
          ibs = PiwinskiSmooth(pnumber, ex, ey, sigs, dponp, header, r0);

          auto buf_out = out.request();
          double *ptr_out = static_cast<double *>(buf_out.ptr);
          ptr_out[0] = ibs[0];
          ptr_out[1] = ibs[1];
          ptr_out[2] = ibs[2];
        },
        "Piwinski smooth", py::arg("pnumber"), py::arg("emitx"),
        py::arg("emity"), py::arg("bunchLength"), py::arg("dpop"),
        py::arg("twissHeaderMap"), py::arg("classicalRadius"),
        py::arg("outputArray"));

  m.def("PiwinskiLattice",
        [](double pnumber, double ex, double ey, double sigs, double dponp,
           map<string, double> &header, map<string, vector<double>> &table,
           double r0, py::array_t<double> out) {
          double *ibs;
          ibs =
              PiwinskiLattice(pnumber, ex, ey, sigs, dponp, header, table, r0);

          auto buf_out = out.request();
          double *ptr_out = static_cast<double *>(buf_out.ptr);
          ptr_out[0] = ibs[0];
          ptr_out[1] = ibs[1];
          ptr_out[2] = ibs[2];
        },
        "Piwinski Lattice", py::arg("pnumber"), py::arg("emitx"),
        py::arg("emity"), py::arg("bunchLength"), py::arg("dpop"),
        py::arg("twissHeaderMap"), py::arg("twissTableMap"),
        py::arg("classicalRadius"), py::arg("outputArray"));

  m.def("PiwinskiLatticeModified",
        [](double pnumber, double ex, double ey, double sigs, double dponp,
           map<string, double> &header, map<string, vector<double>> &table,
           double r0, py::array_t<double> out) {
          double *ibs;
          ibs = PiwinskiLatticeModified(pnumber, ex, ey, sigs, dponp, header,
                                        table, r0);

          auto buf_out = out.request();
          double *ptr_out = static_cast<double *>(buf_out.ptr);
          ptr_out[0] = ibs[0];
          ptr_out[1] = ibs[1];
          ptr_out[2] = ibs[2];
        },
        "Piwinski Lattice Modified", py::arg("pnumber"), py::arg("emitx"),
        py::arg("emity"), py::arg("bunchLength"), py::arg("dpop"),
        py::arg("twissHeaderMap"), py::arg("twissTableMap"),
        py::arg("classicalRadius"), py::arg("outputArray"));

  m.def("Nagaitsev",
        [](double pnumber, double ex, double ey, double sigs, double dponp,
           map<string, double> &header, map<string, vector<double>> &table,
           double r0, py::array_t<double> out) {
          double *ibs;
          ibs = Nagaitsev(pnumber, ex, ey, sigs, dponp, header, table, r0);

          auto buf_out = out.request();
          double *ptr_out = static_cast<double *>(buf_out.ptr);
          ptr_out[0] = ibs[0];
          ptr_out[1] = ibs[1];
          ptr_out[2] = ibs[2];
        },
        "Nagaitsev", py::arg("pnumber"), py::arg("emitx"), py::arg("emity"),
        py::arg("bunchLength"), py::arg("dpop"), py::arg("twissHeaderMap"),
        py::arg("twissTableMap"), py::arg("classicalRadius"),
        py::arg("outputArray"));

  m.def("NagaitsevTailcut",
        [](double pnumber, double ex, double ey, double sigs, double dponp,
           map<string, double> &header, map<string, vector<double>> &table,
           double r0, double aatom, py::array_t<double> out) {
          double *ibs;
          ibs = Nagaitsevtailcut(pnumber, ex, ey, sigs, dponp, header, table,
                                 r0, aatom);

          auto buf_out = out.request();
          double *ptr_out = static_cast<double *>(buf_out.ptr);
          ptr_out[0] = ibs[0];
          ptr_out[1] = ibs[1];
          ptr_out[2] = ibs[2];
        },
        "Nagaitsev Tailcut", py::arg("pnumber"), py::arg("emitx"),
        py::arg("emity"), py::arg("bunchLength"), py::arg("dpop"),
        py::arg("twissHeaderMap"), py::arg("twissTableMap"),
        py::arg("classicalRadius"), py::arg("AtomicMassNumber"),
        py::arg("outputArray"));

  m.def("Zimmerman",
        [](double pnumber, double ex, double ey, double sigs, double sige,
           map<string, double> &header, map<string, vector<double>> &table,
           double r0, bool printout, py::array_t<double> out) {
          double *ibs;
          ibs =
              ibsmadx(pnumber, ex, ey, sigs, sige, header, table, r0, printout);

          auto buf_out = out.request();
          double *ptr_out = static_cast<double *>(buf_out.ptr);
          ptr_out[0] = ibs[0];
          ptr_out[1] = ibs[1];
          ptr_out[2] = ibs[2];
        },
        "Zimmerman", py::arg("pnumber"), py::arg("emitx"), py::arg("emity"),
        py::arg("bunchLength"), py::arg("sige"), py::arg("twissHeaderMap"),
        py::arg("twissTableMap"), py::arg("classicalRadius"),
        py::arg("printoutflag"), py::arg("outputArray"));

  m.def("ZimmermanTailcut",
        [](double pnumber, double ex, double ey, double sigs, double sige,
           map<string, double> &header, map<string, vector<double>> &table,
           double r0, double aatom, py::array_t<double> out) {
          double *ibs;
          ibs = ibsmadxtailcut(pnumber, ex, ey, sigs, sige, header, table, r0,
                               aatom);

          auto buf_out = out.request();
          double *ptr_out = static_cast<double *>(buf_out.ptr);
          ptr_out[0] = ibs[0];
          ptr_out[1] = ibs[1];
          ptr_out[2] = ibs[2];
        },
        "Zimmerman Tailcut", py::arg("pnumber"), py::arg("emitx"),
        py::arg("emity"), py::arg("bunchLength"), py::arg("sige"),
        py::arg("twissHeaderMap"), py::arg("twissTableMap"),
        py::arg("classicalRadius"), py::arg("AtomicMassNumber"),
        py::arg("outputArray"));

  m.def(
      "BjorkenMtingwaSimpson",
      [](double pnumber, double ex, double ey, double sigs, double dponp,
         map<string, double> &header, map<string, vector<double>> &table,
         double r0, py::array_t<double> out) {
        double *ibs;
        ibs = BjorkenMtingwa2(pnumber, ex, ey, sigs, dponp, header, table, r0);

        auto buf_out = out.request();
        double *ptr_out = static_cast<double *>(buf_out.ptr);
        ptr_out[0] = ibs[0];
        ptr_out[1] = ibs[1];
        ptr_out[2] = ibs[2];
      },
      "Bjorken-Mtingwa using Standard Simpson integration.", py::arg("pnumber"),
      py::arg("emitx"), py::arg("emity"), py::arg("bunchLength"),
      py::arg("dpop"), py::arg("twissHeaderMap"), py::arg("twissTableMap"),
      py::arg("classicalRadius"), py::arg("outputArray"));

  m.def("BjorkenMtingwaSimpsonDecade",
        [](double pnumber, double ex, double ey, double sigs, double sige,
           map<string, double> &header, map<string, vector<double>> &table,
           double r0, py::array_t<double> out) {
          double *ibs;
          ibs = BjorkenMtingwa(pnumber, ex, ey, sigs, sige, header, table, r0);

          auto buf_out = out.request();
          double *ptr_out = static_cast<double *>(buf_out.ptr);
          ptr_out[0] = ibs[0];
          ptr_out[1] = ibs[1];
          ptr_out[2] = ibs[2];
        },
        "Bjorken-Mtingwa using Simpson Decade integration.", py::arg("pnumber"),
        py::arg("emitx"), py::arg("emity"), py::arg("bunchLength"),
        py::arg("dpop"), py::arg("twissHeaderMap"), py::arg("twissTableMap"),
        py::arg("classicalRadius"), py::arg("outputArray"));

  m.def("BjorkenMtingwaTailcutSimpsonDecade",
        [](double pnumber, double ex, double ey, double sigs, double dponp,
           map<string, double> &header, map<string, vector<double>> &table,
           double r0, double aatom, py::array_t<double> out) {
          double *ibs;
          ibs = BjorkenMtingwatailcut(pnumber, ex, ey, sigs, dponp, header,
                                      table, r0, aatom);

          auto buf_out = out.request();
          double *ptr_out = static_cast<double *>(buf_out.ptr);
          ptr_out[0] = ibs[0];
          ptr_out[1] = ibs[1];
          ptr_out[2] = ibs[2];
        },
        "Bjorken-Mtingwa using Simpson Decade integration with Tailcut",
        py::arg("pnumber"), py::arg("emitx"), py::arg("emity"),
        py::arg("bunchLength"), py::arg("dpop"), py::arg("twissHeaderMap"),
        py::arg("twissTableMap"), py::arg("classicalRadius"),
        py::arg("AtomicMassNumber"), py::arg("outputArray"));

  m.def("ConteMartiniSimpsonDecade",
        [](double pnumber, double ex, double ey, double sigs, double sige,
           map<string, double> &header, map<string, vector<double>> &table,
           double r0, py::array_t<double> out) {
          double *ibs;
          ibs = ConteMartini(pnumber, ex, ey, sigs, sige, header, table, r0);

          auto buf_out = out.request();
          double *ptr_out = static_cast<double *>(buf_out.ptr);
          ptr_out[0] = ibs[0];
          ptr_out[1] = ibs[1];
          ptr_out[2] = ibs[2];
        },
        "Conte-Martini using Simpson Decade integration.", py::arg("pnumber"),
        py::arg("emitx"), py::arg("emity"), py::arg("bunchLength"),
        py::arg("dpop"), py::arg("twissHeaderMap"), py::arg("twissTableMap"),
        py::arg("classicalRadius"), py::arg("outputArray"));

  m.def("ConteMartiniTailcutSimpsonDecade",
        [](double pnumber, double ex, double ey, double sigs, double sige,
           map<string, double> &header, map<string, vector<double>> &table,
           double r0, double aatom, py::array_t<double> out) {
          double *ibs;
          ibs = ConteMartinitailcut(pnumber, ex, ey, sigs, sige, header, table,
                                    r0, aatom);

          auto buf_out = out.request();
          double *ptr_out = static_cast<double *>(buf_out.ptr);
          ptr_out[0] = ibs[0];
          ptr_out[1] = ibs[1];
          ptr_out[2] = ibs[2];
        },
        "Conte-Martini using Simpson Decade integration with Tailcut",
        py::arg("pnumber"), py::arg("emitx"), py::arg("emity"),
        py::arg("bunchLength"), py::arg("dpop"), py::arg("twissHeaderMap"),
        py::arg("twissTableMap"), py::arg("classicalRadius"),
        py::arg("AtomicMassNumber"), py::arg("outputArray"));

  m.def("ZimmermanSimpsonDecade",
        [](double pnumber, double ex, double ey, double sigs, double dponp,
           map<string, double> &header, map<string, vector<double>> &table,
           double r0, py::array_t<double> out) {
          double *ibs;
          ibs = MadxIBS(pnumber, ex, ey, sigs, dponp, header, table, r0);

          auto buf_out = out.request();
          double *ptr_out = static_cast<double *>(buf_out.ptr);
          ptr_out[0] = ibs[0];
          ptr_out[1] = ibs[1];
          ptr_out[2] = ibs[2];
        },
        "Zimmerman using Simpson Decade", py::arg("pnumber"), py::arg("emitx"),
        py::arg("emity"), py::arg("bunchLength"), py::arg("dpop"),
        py::arg("twissHeaderMap"), py::arg("twissTableMap"),
        py::arg("classicalRadius"), py::arg("outputArray"));

  m.def("runODE",
        [](map<string, double> &twiss, map<string, vector<double>> &twissdata,
           vector<double> h, vector<double> v, vector<double> &t,
           vector<double> &ex, vector<double> &ey, vector<double> &sigs,
           vector<double> sige, int model, double pnumber,
           int couplingpercentage, double threshold, string method) {
          ODE(twiss, twissdata, h.size(), h.data(), v.data(), t, ex, ey, sigs,
              sige, model, pnumber, couplingpercentage, threshold, method);
          map<string, vector<double>> res;
          res["t"] = t;
          res["ex"] = ex;
          res["ey"] = ey;
          res["sigs"] = sigs;
          return res;
        },
        "Run ODE simulation using auto time step.", py::arg("twissheader"),
        py::arg("twisstable"), py::arg("harmonic_rf"), py::arg("voltages_rf"),
        py::arg("t"), py::arg("ex"), py::arg("ey"), py::arg("sigs"),
        py::arg("sige"), py::arg("model"), py::arg("pnumber"),
        py::arg("couplingPercentage"), py::arg("threshold"),
        py::arg("simulationMethod"));
  m.def("runODE",
        [](map<string, double> &twiss, map<string, vector<double>> &twissdata,
           vector<double> h, vector<double> v, vector<double> &t,
           vector<double> &ex, vector<double> &ey, vector<double> &sigs,
           vector<double> sige, int model, double pnumber, int nsteps,
           double stepsize, int couplingpercentage, string method) {
          ODE(twiss, twissdata, h.size(), h.data(), v.data(), t, ex, ey, sigs,
              sige, model, pnumber, nsteps, stepsize, couplingpercentage,
              method);
          map<string, vector<double>> res;
          res["t"] = t;
          res["ex"] = ex;
          res["ey"] = ey;
          res["sigs"] = sigs;
          return res;
        },
        "Run ODE simulation with fixed number of steps and stepsize.",
        py::arg("twissheader"), py::arg("twisstable"), py::arg("harmonic_rf"),
        py::arg("voltages_rf"), py::arg("t"), py::arg("ex"), py::arg("ey"),
        py::arg("sigs"), py::arg("sige"), py::arg("model"), py::arg("pnumber"),
        py::arg("nsteps"), py::arg("stepsize"), py::arg("couplingPercentage"),
        py::arg("simulationMethod"));
}

creating the Python function wrappers around the C++ functions. As before we use CMake to build and install the code. Due to the fact we use a third party library and we need also Python the CMakeLists.txt is more complicated (for more info and a tutorial on how to use this with CMake see Cmake tutorial c++.).

Using the method describe above to include pybind11 one needs to make sure pybind11 can be found by CMake which can be accomplished by adding the following:

SET(MY_PYBIND ${CMAKE_CURRENT_SOURCE_DIR}/third_party/pybind11-2.5.0)
add_subdirectory(${MY_PYBIND})

As we rely on Python their will also be extra code to set the appropriate directories for installing the wrapper function and making them available from within python.

# find python
find_program(
    PYTHON_EXECUTABLE
    NAMES python
)
message("Python Executable ${PYTHON_EXECUTABLE}")
# if python found
# get the PYTHON_LIBRARY_DIR to allow to install files
if(PYTHON_EXECUTABLE)
    message("Python used: ${PYTHON_EXECUTABLE}")
    execute_process(
        COMMAND "${PYTHON_EXECUTABLE}" -c
        "from distutils import sysconfig;print(sysconfig.get_config_var('LIBDIR'))"
            OUTPUT_VARIABLE PYTHON_LIBRARY_DIR
            RESULT_VARIABLE PYTHON_LIBRARY_DIR_RESULT
            ERROR_QUIET)
    string(STRIP ${PYTHON_LIBRARY_DIR} PYTHON_LIBRARY_DIR)
    if(NOT PYTHON_LIBRARY_DIR_RESULT MATCHES 0)
        message(SEND_ERROR "Failed to determine PYTHON_LIBRARY")
    endif()
    execute_process(
        COMMAND "${PYTHON_EXECUTABLE}" -c "import site; print(site.getsitepackages()[0])"
        OUTPUT_VARIABLE _site_packages
        OUTPUT_STRIP_TRAILING_WHITESPACE
        ERROR_QUIET
    )
else()
    message(FATAL_ERROR "python executable not found.")
endif()

The full CMakeLists.txt is given below:

cmake_minimum_required(VERSION 3.10.2)

SET(CMAKE_CXX_STANDARD 17)
SET(CMAKE_CXX_STANDARD_REQUIRED ON)
SET(CMAKE_CXX_EXTENSIONS OFF)
SET(GIT_SUBMODULE ON)

# setting build type if none is given
if(NOT CMAKE_BUILD_TYPE)
    set(CMAKE_BUILD_TYPE Release)
endif()

# setting compiler flags
SET(GCC_OPENMP_COMPILE_FLAGS "-fopenmp -D use_openmp -Ofast -march=native -ffast-math")
SET(GCC_OPENMP_LINK_FLAGS "-fopenmp ")

SET(CMAKE_CXX_FLAGS  "${CMAKE_CXX_FLAGS} ${GCC_OPENMP_COMPILE_FLAGS}")
SET(CMAKE_CXX_FLAGS_RELEASE  "${CMAKE_CXX_FLAGS} ${GCC_OPENMP_COMPILE_FLAGS}")
SET(CMAKE_EXE_LINKER_FLAGS  "${CMAKE_EXE_LINKER_FLAGS} ${GCC_OPENMP_LINK_FLAGS}")


PROJECT(IBSLib VERSION 0.2.0)

# DOWNLOAD THE GIT SUBMODULES (PYBIND11 HERE)
find_package(Git QUIET)

if (GIT_FOUND AND EXISTS "${PROJECT_SOURCE_DIR}/.git")
    # update as needed
    option(GIT_SUBMODULE "Check submodules during build" ON)

    if(GIT_SUBMODULE)
        message(STATUS "Updating submodules...")
        execute_process(COMMAND ${GIT_EXECUTABLE} submodule update --init --recursive
            WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
            RESULT_VARIABLE GIT_SUBMOD_RESULT)
        if(NOT GIT_SUBMOD_RESULT EQUAL "0")
            message(FATAL_ERROR "git submodule update --init failed with ${GIT_SUBMOD_RESULT}, please check...")
        endif()
    endif()
endif()

# CHECK ALL THE SUBMODULES
if(NOT EXISTS "${PROJECT_SOURCE_DIR}/third_party/pybind11-2.6.2/CMakeLists.txt")
    message(FATAL_ERROR "The PYBIND11 submodule was not downloaded! GIT_SUBMODULE was turned off or failed.")
endif()

# find python
find_program(
    PYTHON_EXECUTABLE
    NAMES python
)

# if python found 
# get the PYTHON_LIBRARY_DIR to allow to install files
if(PYTHON_EXECUTABLE)
    message(STATUS "PYTHON_EXECUTABLE = ${PYTHON_EXECUTABLE}")

    execute_process(
        COMMAND "${PYTHON_EXECUTABLE}" -c
        "from distutils import sysconfig;print(sysconfig.get_config_var('LIBDIR'))"
            OUTPUT_VARIABLE PYTHON_LIBRARY_DIR
            RESULT_VARIABLE PYTHON_LIBRARY_DIR_RESULT
            ERROR_QUIET)    
    string(STRIP ${PYTHON_LIBRARY_DIR} PYTHON_LIBRARY_DIR)
    if(NOT PYTHON_LIBRARY_DIR_RESULT MATCHES 0)
        message(SEND_ERROR "Failed to determine PYTHON_LIBRARY")
    endif()
    execute_process(
        COMMAND "${PYTHON_EXECUTABLE}" -c "import site; print(site.getsitepackages()[0])"
        OUTPUT_VARIABLE _site_packages
        OUTPUT_STRIP_TRAILING_WHITESPACE
        ERROR_QUIET
    )
else()
    message(FATAL_ERROR "python executable not found.")
endif()


# adding pybind11
# git submodule add https://github.com/pybind/pybind11.git third_party/pybind11-2.5.0
SET(MY_PYBIND ${CMAKE_CURRENT_SOURCE_DIR}/third_party/pybind11-2.6.2)
add_subdirectory(${MY_PYBIND})


include_directories("${CMAKE_SOURCE_DIR}/cpp/include/ibs_bits")
include_directories("${CMAKE_SOURCE_DIR}/python")

file (GLOB SOURCE_FILES "cpp/src/*.cpp")
file (GLOB HEADER_FILES "cpp/include/ibs_bits/*.hpp")
file (GLOB PYTHON_FILES "python/*.cpp" "python/*.hpp")

message(STATUS "SOURCE_FILES = ${SOURCE_FILES}")
message(STATUS "HEADER_FILES = ${HEADER_FILES}")
message(STATUS "PYTHON_FILES = ${PYTHON_FILES}")

# Set up such that XCode organizes the files
source_group(TREE ${CMAKE_CURRENT_SOURCE_DIR} FILES ${SOURCE_FILES} ${HEADER_FILES} ${PYTHON_FILES} )

# include(pybind11.cmake)
pybind11_add_module(${PROJECT_NAME} 
${SOURCE_FILES}
${HEADER_FILES}
${PYTHON_FILES}
)


target_link_libraries(${PROJECT_NAME}  PUBLIC)

if(NOT _site_packages )
    set(_site_packages  ${PYTHON_LIBRARY_DIR})
endif()

# show the PYTHON_LIBRARY_DIR
message(STATUS "PYTHON_LIBARAY_DIR = ${PYTHON_LIBRARY_DIR}")
message(STATUS "_site_package      = ${_site_packages}")

install(TARGETS ${PROJECT_NAME} 
COMPONENT python
#LIBRARY DESTINATION "${PYTHON_LIBRARY_DIR}"
LIBRARY DESTINATION "${_site_packages}"
)