text:
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_expint.h>
#include "ccl.h"

//maths helper function for NFW profile
static double helper_fx(double x){
    double f;
    f = log(1.+x)-x/(1.+x);
    return f;
}

//cosmo: ccl cosmology object containing cosmological parameters
//c: halo concentration, needs to be consistent with halo size definition
//halomass: halo mass
//massdef_delta: mass definition, overdensity relative to mean matter density
//a: scale factor
//r: radii at which to calculate output
//nr: number of radii for calculation
//rho_r: stores densities at r
//returns void
void ccl_halo_profile_nfw(ccl_cosmology *cosmo, double c, double halomass,
                          double massdef_delta_m, double a, double *r, int nr,
                          double *rho_r, int *status) {

    //haloradius: halo radius for mass definition
    //rs: scale radius
    double haloradius, rs;

    haloradius = r_delta(cosmo, halomass, a, massdef_delta_m, status);
    rs = haloradius/c;

    //rhos: NFW density parameter
    double rhos;

    rhos = halomass/(4.*M_PI*(rs*rs*rs)*helper_fx(c));

    int i;
    double x;
    for(i=0; i < nr; i++) {
        x = r[i]/rs;
        rho_r[i] = rhos/(x*(1.+x)*(1.+x));
    }
    return;
}

//cosmo: ccl cosmology object containing cosmological parameters
//c: halo concentration, needs to be consistent with halo size definition
//halomass: halo mass
//massdef_delta: mass definition, overdensity relative to mean matter density
//a: scale factor
//rp: radius at which to calculate output
//nr: number of radii for calculation
//sigma_r: stores surface mass density (integrated along line of sight) at given projected rp
//returns void
void ccl_projected_halo_profile_nfw(ccl_cosmology *cosmo, double c,
                                    double halomass, double massdef_delta_m,
                                    double a, double *rp, int nr, double *sigma_r,
                                    int *status) {

    //haloradius: halo radius for mass definition
    //rs: scale radius
    double haloradius, rs;

    haloradius = r_delta(cosmo, halomass, a, massdef_delta_m, status);
    rs = haloradius/c;

    //rhos: NFW density parameter
    double rhos;

    rhos = halomass/(4.*M_PI*(rs*rs*rs)*helper_fx(c));

    double x;

    int i;
    for(i=0; i < nr; i++){

        x = rp[i]/rs;
        if (x==1.){
            sigma_r[i] = 2.*rs*rhos/3.;
        }
        else if (x<1.){
            sigma_r[i] = 2.*rs*rhos*(1.-2.*atanh(sqrt(fabs((1.-x)/(1.+x))))/sqrt(fabs(1.-x*x)))/(x*x-1.);
        }
        else {
            sigma_r[i] = 2.*rs*rhos*(1.-2.*atan(sqrt(fabs((1.-x)/(1.+x))))/sqrt(fabs(1.-x*x)))/(x*x-1.);
        }
    }
    return;

}

//maths helper function assuming NFW approximation
static double helper_solve_cvir(double c, void *rhs_pointer){
    double rhs = *(double*)rhs_pointer;
    return (log(1.+c)-c/(1.+c))/(c*c*c) - rhs;
}

//solve for cvir from different mass definition iteratively, assuming NFW profile.
static double solve_cvir(ccl_cosmology *cosmo, double rhs, double initial_guess, int *status){

    int rootstatus;
    int iter = 0, max_iter = 100;
    const gsl_root_fsolver_type *T;
    gsl_root_fsolver *s;
    double cvir;
    double c_lo = 0.1, c_hi = initial_guess*100.;
    gsl_function F;

    F.function = &helper_solve_cvir;
    F.params = &rhs;

    T = gsl_root_fsolver_brent;
    s = gsl_root_fsolver_alloc (T);
    if (s != NULL) {
      gsl_root_fsolver_set (s, &F, c_lo, c_hi);

      do
        {
          iter++;
          gsl_root_fsolver_iterate (s);
          cvir = gsl_root_fsolver_root (s);
          c_lo = gsl_root_fsolver_x_lower (s);
          c_hi = gsl_root_fsolver_x_upper (s);
          rootstatus = gsl_root_test_interval (c_lo, c_hi, 0, 0.0001);
        }
      while (rootstatus == GSL_CONTINUE && iter < max_iter);

      gsl_root_fsolver_free(s);

      // Check for errors
      if (rootstatus != GSL_SUCCESS) {
          ccl_raise_gsl_warning(rootstatus, "ccl_haloprofile.c: solve_cvir():");
          *status = CCL_ERROR_PROFILE_ROOT;
          ccl_cosmology_set_status_message(cosmo, "ccl_haloprofile.c: solve_cvir(): Root finding failure\n");
          return NAN;
      } else {
          return cvir;
      }
    }
    else {
      *status = CCL_ERROR_MEMORY;
      return NAN;
    }
}

// Structure to hold parameters of integrand_einasto
typedef struct{
  ccl_cosmology *cosmo;
  double alpha, rs;
} Int_Einasto_Par;

static double integrand_einasto(double r, void *params){

    Int_Einasto_Par *p = (Int_Einasto_Par *)params;
    double alpha = p->alpha;
    double rs = p->rs;

    return 4.*M_PI*r*r*exp(-2.*(pow(r/rs,alpha)-1)/alpha);
}

//integrate einasto profile to get normalization of rhos
static double integrate_einasto(ccl_cosmology *cosmo, double R, double alpha, double rs, int *status){
    int qagstatus;
    double result = 0, eresult;
    Int_Einasto_Par ipar;
    gsl_function F;
    gsl_integration_workspace *w = NULL;

    w = gsl_integration_workspace_alloc(1000);

    if (w == NULL) {
      *status = CCL_ERROR_MEMORY;
    }

    if (*status == 0) {
      // Structure required for the gsl integration
      ipar.cosmo = cosmo;
      ipar.alpha = alpha;
      ipar.rs = rs;
      F.function = &integrand_einasto;
      F.params = &ipar;

      // Actually does the integration
      qagstatus = gsl_integration_qag(
        &F, 0, R, 0, 0.0001, 1000, 3, w,
        &result, &eresult);

      // Check for errors
      if (qagstatus != GSL_SUCCESS) {
        ccl_raise_gsl_warning(qagstatus, "ccl_haloprofile.c: integrate_einasto():");
        *status = CCL_ERROR_PROFILE_INT;
        ccl_cosmology_set_status_message(cosmo, "ccl_haloprofile.c: integrate_einasto(): Integration failure\n");
        result = NAN;
      }
    }

  // Clean up
  gsl_integration_workspace_free(w);

  return result;
}

//cosmo: ccl cosmology object containing cosmological parameters
//c: halo concentration, needs to be consistent with halo size definition
//halomass: halo mass
//massdef_delta: mass definition, overdensity relative to mean matter density
//a: scale factor
//r: radii at which to calculate output
//nr: number of radii for calculation
//rho_r: stores densities at r
//returns void
void ccl_halo_profile_einasto(ccl_cosmology *cosmo, double c, double halomass,
                              double massdef_delta_m, double a, double *r, int nr,
                              double *rho_r, int *status) {

    //haloradius: halo radius for mass definition
    //rs: scale radius
    double haloradius, rs;

    haloradius = r_delta(cosmo, halomass, a, massdef_delta_m, status);
    rs = haloradius/c;

    //nu: peak height, https://arxiv.org/pdf/1401.1216.pdf eqn1
    //alpha: Einasto parameter, https://arxiv.org/pdf/1401.1216.pdf eqn5
    double nu;
    double alpha;  //calibrated relation with nu, with virial mass
    double Mvir;
    double Delta_v;

    Delta_v = Dv_BryanNorman(cosmo, a, status); //virial definition, if odelta is this, the definition is virial.

    if (massdef_delta_m<Delta_v+1 && massdef_delta_m>Delta_v-1){   //allow rounding of virial definition
        Mvir = halomass;
    }
    else{
        double rhs; //NFW equation for cvir: f(Rvir/Rs)/f(c)=(Rvir^3*Delta_v)/(R^3*massdef) -> f(cvir)/cvir^3 = rhs
        rhs = (helper_fx(c)*(rs*rs*rs)*Delta_v)/((haloradius*haloradius*haloradius)*massdef_delta_m);
        Mvir = halomass*helper_fx(solve_cvir(cosmo, rhs, c, status))/helper_fx(c);
    }

    nu = 1.686/ccl_sigmaM(cosmo, log10(Mvir), a, status); //delta_c_Tinker
    alpha = 0.155 + 0.0095*nu*nu;

    //rhos: scale density
    double rhos;

    rhos = halomass/integrate_einasto(cosmo, haloradius, alpha, rs, status); //normalize

    int i;
    for(i=0; i < nr; i++) {
        rho_r[i] = rhos*exp(-2.*(pow(r[i]/rs,alpha)-1.)/alpha);
    }

    return;

}

// Structure to hold parameters of integrand_hernquist
typedef struct{
  ccl_cosmology *cosmo;
  double rs;
} Int_Hernquist_Par;

static double integrand_hernquist(double r, void *params){

    Int_Hernquist_Par *p = (Int_Hernquist_Par *)params;
    double rs = p->rs;

    return 4.*M_PI*r*r/((r/rs)*(1.+r/rs)*(1.+r/rs)*(1.+r/rs));
}

//integrate hernquist profile to get normalization of rhos
static double integrate_hernquist(ccl_cosmology *cosmo, double R, double rs, int *status){
    int qagstatus;
    double result = 0, eresult;
    Int_Hernquist_Par ipar;
    gsl_function F;
    gsl_integration_workspace *w = NULL;

    w = gsl_integration_workspace_alloc(1000);

    if (w == NULL) {
      *status = CCL_ERROR_MEMORY;
    }

    if (*status == 0) {
      // Structure required for the gsl integration
      ipar.cosmo = cosmo;
      ipar.rs = rs;
      F.function = &integrand_hernquist;
      F.params = &ipar;

      // Actually does the integration
      qagstatus = gsl_integration_qag(
        &F, 0, R, 0, 0.0001, 1000, 3, w,
        &result, &eresult);

      // Check for errors
      if (qagstatus != GSL_SUCCESS) {
        ccl_raise_gsl_warning(qagstatus, "ccl_haloprofile.c: integrate_hernquist():");
        *status = CCL_ERROR_PROFILE_INT;
        ccl_cosmology_set_status_message(cosmo, "ccl_haloprofile.c: integrate_hernquist(): Integration failure\n");
        result = NAN;
      }
    }

    // Clean up
    gsl_integration_workspace_free(w);

    return result;
}

//cosmo: ccl cosmology object containing cosmological parameters
//c: halo concentration, needs to be consistent with halo size definition
//halomass: halo mass
//massdef_delta: mass definition, overdensity relative to mean matter density
//a: scale factor
//r: radii at which to calculate output
//nr: number of radii for calculation
//rho_r: stores densities at r
//returns void
void ccl_halo_profile_hernquist(ccl_cosmology *cosmo, double c, double halomass,
                                double massdef_delta_m, double a, double *r, int nr,
                                double *rho_r, int *status) {

    //haloradius: halo radius for mass definition
    //rs: scale radius
    double haloradius, rs;

    haloradius = r_delta(cosmo, halomass, a, massdef_delta_m, status);
    rs = haloradius/c;

    //rhos: scale density
    double rhos;

    rhos = halomass/integrate_hernquist(cosmo, haloradius, rs, status); //normalize

    int i;
    double x;
    for(i=0; i < nr; i++) {
        x = r[i]/rs;
        rho_r[i] = rhos/(x*pow((1.+x),3));
    }

    return;

}

max_stars_repo_name: benediktdiemer/CCL

max_issues_repo_name: benediktdiemer/CCL

max_forks_repo_name: benediktdiemer/CCL


text:
#ifndef CWANNIER_DOS_UTIL
#define CWANNIER_DOS_UTIL

#include <stdlib.h>
#include <stdio.h>
#include "bstrlib/bstrlib.h"
#include <gsl/gsl_matrix.h>

gsl_matrix* parse_R_from_bs(char *b1, char *b2, char *b3);
int set_R_row(gsl_matrix *R, char *b, int row);
double* linspace(double a, double b, int num);

#endif //CWANNIER_DOS_UTIL

ext: h

max_stars_repo_name: tflovorn/cwannier

max_issues_repo_name: tflovorn/cwannier

max_forks_repo_name: tflovorn/cwannier


text:
#pragma once

#ifndef UPCA_ARCH_H
#error "Don't include this file directly. Use #include <upca/upca.h>."
#endif

#include <iostream>
#include <memory>
#include <sstream>
#include <string>

#include <gsl/gsl>

#include <fcntl.h>
#include <stdio.h>
#include <sys/stat.h>
#include <sys/types.h>

#ifdef JEVENTS_FOUND

#include <inttypes.h>
#include <unistd.h>

#ifdef __cplusplus
extern "C" {
#endif
#include <jevents.h>
#include <rdpmc.h>
#ifdef __cplusplus
}
#endif

#include <linux/perf_event.h>
#include <sys/ioctl.h>

#endif

namespace upca {

std::ostream &operator<<(std::ostream &os, const struct perf_event_attr &c);

namespace arch {
namespace x86_64 {

class x86_64_base_pmu {
public:
  uint64_t timestamp_begin() {
    unsigned high, low;
    __asm__ volatile("CPUID\n\t"
                     "RDTSC\n\t"
                     "mov %%edx, %0\n\t"
                     "mov %%eax, %1\n\t"
                     : "=r"(high), "=r"(low)::"%rax", "%rbx", "%rcx", "%rdx");

    return static_cast<uint64_t>(high) << 32 | low;
  }

  uint64_t timestamp_end() {
    unsigned high, low;
    __asm__ volatile("RDTSCP\n\t"
                     "mov %%edx,%0\n\t"
                     "mov %%eax,%1\n\t"
                     "CPUID\n\t"
                     : "=r"(high), "=r"(low)::"%rax", "%rbx", "%rcx", "%rdx");
    return static_cast<uint64_t>(high) << 32 | low;
  }
};

#ifdef JEVENTS_FOUND

#ifndef __APPLE__

/*
 * SYSCALL_HANDLED(601, pmc_init) int counter, int type, int mode
 * SYSCALL_HANDLED(602, pmc_start)
 * SYSCALL_HANDLED(603, pmc_stop)
 * SYSCALL_HANDLED(604, pmc_reset)
 *
 * Yes, the interface is that fucked up. counter is either int or unsigned long,
 * depending on the syscall …
 */

static inline long mck_pmc_init(int counter, int type, unsigned mode) {
  return syscall(601, counter, type, mode);
}

static inline long mck_pmc_start(unsigned long counter) {
  return syscall(602, counter);
}
static inline long mck_pmc_stop(unsigned long counter) { return syscall(603, counter); }
static inline long mck_pmc_reset(int counter) { return syscall(604, counter); }

static inline int mck_is_mckernel() { return syscall(732) == 0; }

#else

static inline int mck_is_mckernel() { return 0; }

#endif

class resolver {
public:
  resolver();

  struct perf_event_attr resolve(const std::string &name) const;
};

class fd {
  int fd_;

public:
  explicit fd(int fd) : fd_(fd) {}
  fd(std::nullptr_t = nullptr) : fd_(-1) {}
  ~fd() {
    if (fd_ != -1) {
      close(fd_);
    }
  }
  fd(fd &&l) {
    fd_ = l.fd_;
    l.fd_ = -1;
  }
  explicit operator bool() { return fd_ != -1; }
  explicit operator int() { return fd_; }
  friend bool operator==(const fd &l, const fd &r) { return l.fd_ == r.fd_; }
  friend bool operator!=(const fd &l, const fd &r) { return !(l == r); }
};

static inline std::unique_ptr<fd> make_perf_fd(struct perf_event_attr &attr) {
  const int fd = perf_event_open(&attr, 0, -1, -1, 0);

  if (fd < 0) {
    std::ostringstream os;
    os << "perf_event_open() failed: " << errno << strerror(errno) << std::endl;
    throw std::runtime_error(os.str());
  }
  return std::make_unique<class fd>(fd);
}

static uint64_t rdpmc(const uint32_t counter) {
  uint32_t low, high;
  __asm__ volatile("rdpmc" : "=a"(low), "=d"(high) : "c"(counter));
  return static_cast<uint64_t>(high) << 32 | low;
}

enum msrs {
  IA32_PMC_BASE = 0x0c1,
  IA32_PERFEVTSEL_BASE = 0x186,
  IA32_FIXED_CTR_CTRL = 0x38d,
  IA32_PERF_GLOBAL_STATUS = 0x38e,
  IA32_PERF_GLOBAL_CTRL = 0x38f,
  IA32_PERF_GLOBAL_OVF_CTRL = 0x390,
  IA32_PERF_CAPABILITIES = 0x345,
  IA32_DEBUGCTL = 0x1d9,
};

static inline void cpuid(const int code, uint32_t *a, uint32_t *b, uint32_t *c,
                         uint32_t *d) {
  asm volatile("cpuid" : "=a"(*a), "=b"(*b), "=c"(*c), "=d"(*d) : "a"(code));
}

#define MASK(v, high, low) ((v >> low) & ((1 << (high - low + 1)) - 1))

template <typename POLICY> static void pmu_info(const POLICY &msr) {
  uint32_t eax, ebx, ecx, edx;

  cpuid(0x1, &eax, &ebx, &ecx, &edx);

  std::cout << std::hex;

  std::cout << "EAX: " << eax << std::endl;
  std::cout << "EBX: " << ebx << std::endl;
  std::cout << "ECX: " << ecx << std::endl;
  std::cout << "EDX: " << edx << std::endl;

  if (ecx & (1 << 15)) {
    const uint64_t caps = msr.rdmsr(IA32_PERF_CAPABILITIES);
    const int vmm_freeze = MASK(caps, 12, 12);
    std::cout << "Caps: " << caps << std::endl;
    std::cout << "VMM Freeze: " << std::boolalpha
              << static_cast<bool>(vmm_freeze) << std::noboolalpha << std::endl;
    if (vmm_freeze) {
      const uint64_t debugctl = msr.rdmsr(IA32_DEBUGCTL);
      msr.wrmsr(IA32_DEBUGCTL, debugctl & ~(1 << 14));
    }
  }

  cpuid(0xa, &eax, &ebx, &ecx, &edx);

  const unsigned version = MASK(eax, 7, 0);
  const unsigned counters = MASK(eax, 15, 8);
  const unsigned width = MASK(eax, 23, 16);

  const unsigned ffpc = MASK(edx, 4, 0);
  const unsigned ff_width = MASK(edx, 12, 5);

  std::cout << "EAX: " << eax << std::endl;
  std::cout << "EBX: " << ebx << std::endl;
  std::cout << "ECX: " << ecx << std::endl;
  std::cout << "EDX: " << edx << std::endl;

  std::cout << std::dec;

  std::cout << "PMC version " << version << std::endl;
  std::cout << "PMC counters: " << counters << std::endl;
  std::cout << "PMC width: " << width << std::endl;

  std::cout << "FFPCs: " << ffpc << ", width: " << ff_width << std::endl;
}

struct x86_pmc_base {
  virtual ~x86_pmc_base();
  virtual gsl::span<uint64_t>::index_type start(gsl::span<uint64_t>) = 0;
  virtual gsl::span<uint64_t>::index_type stop(gsl::span<uint64_t>) = 0;
};

struct mck_mck_policy {
  static long init(const int active, const uint64_t config) {
    return mck_pmc_init(active, gsl::narrow<int>(config), 0x4);
  }

  /* writes MSR_IA32_PMC0 + i to 0 */
  static void reset(const int i) { mck_pmc_reset(i); }
  /* sets and clears bits in MSR_PERF_GLOBAL_CTRL */
  static void start(const unsigned long mask) {
    const auto err = mck_pmc_start(mask);
    if (err) {
      std::cerr << "Error starting PMCs" << std::endl;
    }
  }
  static void stop(const unsigned long mask) {
    const auto err = mck_pmc_stop(mask);
    if (err) {
      std::cerr << "Error stopping PMCs" << std::endl;
    }
  }
  static uint64_t read(const unsigned i) { return rdpmc(i); }
};

class msr_mck {
public:
  msr_mck() {}
  uint64_t rdmsr(const uint32_t reg) const {
    return static_cast<unsigned long>(syscall(850, reg));
  }
  auto wrmsr(const uint32_t reg, const uint64_t val) const {
    /* always returns 0 */
    return syscall(851, reg, val);
  }
};

class msr_linux {
  std::unique_ptr<fd> fd_;

  static unsigned cpu() {
    int cpu = sched_getcpu();
    if (cpu < 0) {
      using namespace std::string_literals;
      throw std::runtime_error("Error getting CPU number: "s + strerror(errno) +
                               " (" + std::to_string(errno) + ")\n");
    }
    /* TODO: check thread affinity mask
     * - has a single CPU
     * - this single CPU is this CPU
     */
    return static_cast<unsigned>(cpu);
  }

public:
  msr_linux()
      : fd_(std::make_unique<fd>(open(
            ("/dev/cpu/" + std::to_string(cpu()) + "/msr").c_str(), O_RDWR))) {}

  uint64_t rdmsr(const uint32_t reg) const {
    uint64_t v;
    const ssize_t ret = pread(static_cast<int>(*fd_), &v, sizeof(v), reg);
    if (ret != sizeof(v)) {
      std::cerr << "Reading MSR " << std::hex << reg << std::dec << " failed."
                << std::endl;
    }
    return v;
  }

  auto wrmsr(const uint32_t reg, const uint64_t val) const {
    const ssize_t ret = pwrite(static_cast<int>(*fd_), &val, sizeof(val), reg);
    if (ret != sizeof(val)) {
      std::cerr << "Writing MSR " << std::hex << reg << std::dec << " failed."
                << std::endl;
    }
    return ret;
  }
};

template <typename MSR> class rawmsr_policy : MSR {
  uint64_t global_ctrl_;

public:
  rawmsr_policy() : global_ctrl_(this->rdmsr(IA32_PERF_GLOBAL_CTRL)) {
    // pmu_info(static_cast<MSR&>(*this));
  }
  ~rawmsr_policy() { this->wrmsr(IA32_PERF_GLOBAL_CTRL, global_ctrl_); }
  int init(const int i, const uint64_t config) {
    const uint64_t v = (1 << 22) | (1 << 16) | config;
    this->wrmsr(IA32_PERFEVTSEL_BASE + static_cast<uint32_t>(i), v);
    return 0;
  }

  void reset(const int i) {
    this->wrmsr(IA32_PMC_BASE + static_cast<uint32_t>(i), 0);
  }
  void start(const unsigned long mask) {
    this->wrmsr(IA32_PERF_GLOBAL_OVF_CTRL, 0);
    this->wrmsr(IA32_PERF_GLOBAL_CTRL, mask);
  }
  void stop(const unsigned long /* mask */) {
    this->wrmsr(IA32_PERF_GLOBAL_CTRL, 0);
    const uint64_t ovf = this->rdmsr(IA32_PERF_GLOBAL_STATUS);
    if (ovf) {
      std::cout << "Overflow: " << std::hex << ovf << std::dec << std::endl;
      this->wrmsr(IA32_PERF_GLOBAL_OVF_CTRL, 0);
    }
  }
  /* Should be equivalent to rdpmc instruction */
  uint64_t read(const unsigned i) { return this->rdmsr(IA32_PMC_BASE + i); }
};

template <typename ACCESS> class msr_pmc final : ACCESS, public x86_pmc_base {
  unsigned long active_mask_ = 0;
  int active = 0;

public:
  using resolver_type = resolver;

  template <typename T> msr_pmc(const T &pmcs) {
    for (const auto &pmc : pmcs) {
      const auto err = ACCESS::init(active, pmc.data().config);
      if (err) {
        std::cerr << "Error configuring PMU " << pmc.name() << " with "
                  << std::hex << pmc.data().config << std::dec << std::endl;
        continue;
      }
      ++active;
    }

    active_mask_ = static_cast<unsigned long>((1 << active) - 1);
  }

  gsl::span<uint64_t>::index_type start(gsl::span<uint64_t>) override {
    for (int i = 0; i < active; ++i) {
      ACCESS::reset(i);
    }
    ACCESS::start(active_mask_);

    return active;
  }

  gsl::span<uint64_t>::index_type stop(gsl::span<uint64_t> buf) override {
    ACCESS::stop(active_mask_);
    for (int i = 0; i < active; ++i) {
      buf[i] = ACCESS::read(static_cast<uint32_t>(i));
    }
    return active;
  }
};

using mck_rawmsr = msr_pmc<rawmsr_policy<msr_mck>>;
using mck_mckmsr = msr_pmc<mck_mck_policy>;
using linux_rawmsr = msr_pmc<rawmsr_policy<msr_linux>>;

class linux_jevents final : public x86_pmc_base {
  struct rdpmc_t {
    struct rdpmc_ctx ctx;
    int close = 0;

    rdpmc_t(struct perf_event_attr attr) {
      const int err = rdpmc_open_attr(&attr, &ctx, nullptr);
      if (err) {
        throw std::runtime_error(std::to_string(errno) + ' ' + strerror(errno));
      }

      close = 1;
    }
    ~rdpmc_t() {
      if (close) {
        rdpmc_close(&ctx);
      }
    }
    rdpmc_t(const rdpmc_t &) = delete;
    rdpmc_t(rdpmc_t &&rhs) {
      this->ctx = rhs.ctx;
      this->close = rhs.close;
      rhs.close = 0;
    }

    uint64_t read() { return rdpmc_read(&ctx); }
  };

  std::vector<rdpmc_t> rdpmc_ctxs;

public:
  using resolver_type = resolver;

  template <typename T> linux_jevents(const T &pmcs) {
    for (const auto &pmc : pmcs) {
      rdpmc_ctxs.emplace_back(pmc.data());
    }
  }

  gsl::span<uint64_t>::index_type start(gsl::span<uint64_t> buf) override {
    int idx = 0;
    for (auto &&pmc : rdpmc_ctxs) {
      buf[idx] = pmc.read();
      ++idx;
    }
    return idx;
  }

  gsl::span<uint64_t>::index_type stop(gsl::span<uint64_t> buf) override {
    int idx = 0;
    for (auto &&pmc : rdpmc_ctxs) {
      buf[idx] = pmc.read() - buf[idx];
      ++idx;
    }
    return idx;
  }
};

class linux_perf final : public x86_pmc_base {
  std::vector<std::unique_ptr<fd>> perf_fds;

public:
  using resolver_type = resolver;

  template <typename T> linux_perf(const T &pmcs) {
    for (const auto &pmc : pmcs) {
      struct perf_event_attr pe;
      memset(&pe, 0, sizeof(pe));
      const auto &attr = pmc.data();
      pe.config = attr.config;
      pe.config1 = attr.config1;
      pe.config2 = attr.config2;
      pe.type = attr.type;
      pe.size = attr.size;
      // attr.disabled = 1;
      pe.exclude_kernel = 1;
      pe.exclude_hv = 1;
      try {
        perf_fds.push_back(make_perf_fd(pe));
      } catch (const std::exception &e) {
        std::cerr << pmc.name() << ": " << e.what() << std::endl;
      }
    }
  }

  gsl::span<uint64_t>::index_type start(gsl::span<uint64_t>) override {
    for (const auto &perf_fd : perf_fds) {
      if (ioctl(static_cast<int>(*perf_fd), PERF_EVENT_IOC_RESET, 0)) {
        std::cerr << "Error in ioctl\n";
      }
    }

    return static_cast<gsl::span<uint64_t>::index_type>(perf_fds.size());
  }

  gsl::span<uint64_t>::index_type stop(gsl::span<uint64_t> buf) override {
    uint64_t v;
    int idx = 0;
    for (const auto &perf_fd : perf_fds) {
      const auto ret = read(static_cast<int>(*perf_fd), &v, sizeof(v));
      if (ret < 0) {
        std::cerr << "perf: Error reading performance counter: " << errno
                  << ": " << strerror(errno) << std::endl;
      } else if (ret != sizeof(v)) {
        std::cerr << "perf: Error reading " << sizeof(v) << " bytes."
                  << std::endl;
      }
      buf[idx] = v;
      ++idx;
    }
    return idx;
  }
};

template <typename MCK, typename LINUX>
class x86_linux_mckernel : public x86_64_base_pmu {
  std::unique_ptr<x86_pmc_base> backend_;

public:
  using resolver_type = resolver;

  template <typename T>
  x86_linux_mckernel(const T &pmcs)
      : backend_(mck_is_mckernel() ? static_cast<std::unique_ptr<x86_pmc_base>>(
                                         std::make_unique<MCK>(pmcs))
                                   : static_cast<std::unique_ptr<x86_pmc_base>>(
                                         std::make_unique<LINUX>(pmcs))) {}

  auto start(gsl::span<uint64_t> o) { return backend_->start(o); }
  auto stop(gsl::span<uint64_t> o) { return backend_->stop(o); }
};

#else

struct Reason {
  static constexpr const char reason[] =
      "PMC support not compiled in; libjevents missing";
};

class x86_64_pmu : public x86_64_base_pmu {
public:
  using resolver_type = upca::arch::detail::null_resolver<Reason>;
  template <typename T> x86_64_pmu(const T &) {}

  gsl::span<uint64_t>::index_type start(gsl::span<uint64_t>) { return 0; }
  gsl::span<uint64_t>::index_type stop(gsl::span<uint64_t>) { return 0; }
};

#endif /* JEVENTS_FOUND */

} // namespace x86_64
} // namespace arch
} // namespace upca

max_stars_repo_name: hannesweisbach/ucpa

max_issues_repo_name: hannesweisbach/ucpa

max_forks_repo_name: hannesweisbach/ucpa


text:
/******************************************************************************
CosmoLike Configuration Space Covariances for Projected Galaxy 2-Point Statistics
https://github.com/CosmoLike/CosmoCov
by CosmoLike developers
******************************************************************************/

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>

#include <fftw3.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_sf_trig.h>
#include <gsl/gsl_sf_erf.h>
#define NR_END 1
#define FREE_ARG char*

#define EXIT_MISSING_FILE(ein, purpose, filename) if (!ein) {fprintf(stderr, "Could not find %s file %s\n",purpose,filename);exit(1);}

static double darg __attribute__((unused)),maxarg1 __attribute__((unused)), maxarg2 __attribute__((unused));

#define FMAX(a,b) (maxarg1=(a), maxarg2=(b), (maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2))
#define FMIN(a,b) (maxarg1=(a), maxarg2=(b), (maxarg1) < (maxarg2) ? (maxarg1) : (maxarg2))

//Note: SQR*SQR....gives undefined warning
static double sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
static double cubearg;
#define CUBE(a) ((cubearg=(a)) == 0.0 ? 0.0 : cubearg*cubearg*cubearg)
static double pow4arg;
#define POW4(a) ((pow4arg=(a)) == 0.0 ? 0.0 : pow4arg*pow4arg*pow4arg*pow4arg)

void SVD_inversion(gsl_matrix *cov, gsl_matrix *inverseSVD,int Nmatrix);
double interpol2d(double **f, int nx, double ax, double bx, double dx, double x,   int ny, double ay, double by, double dy, double y, double lower, double upper);
double interpol2d_fitslope(double **f, int nx, double ax, double bx, double dx, double x, int ny, double ay, double by, double dy, double y, double lower);

double interpol(double *f, int n, double a, double b, double dx, double x, double lower, double upper);
double interpol_fitslope(double *f, int n, double a, double b, double dx, double x, double lower);
void free_double_vector(double *v, long nl, long nh);
long *long_vector(long nl, long nh);
int *int_vector(long nl, long nh);
double *create_double_vector(long nl, long nh);
void free_double_matrix(double **m, long nrl, long nrh, long ncl, long nch);
double **create_double_matrix(long nrl, long nrh, long ncl, long nch);
int line_count(char *filename);

void error(char *s);
void cdgamma(fftw_complex x, fftw_complex *res);

double int_gsl_integrate_insane_precision(double (*func)(double, void*),void *arg,double a, double b, double *error, int niter);
double int_gsl_integrate_high_precision(double (*func)(double, void*),void *arg,double a, double b, double *error, int niter);
double int_gsl_integrate_medium_precision(double (*func)(double, void*),void *arg,double a, double b, double *error, int niter);
double int_gsl_integrate_low_precision(double (*func)(double, void*),void *arg,double a, double b, double *error, int niter);
double int_gsl_integrate_cov_precision(double (*func)(double, void*),void *arg,double a, double b, double *error, int niter);

typedef struct {
     double pi;
     double pi_sqr;
     double twopi;
     double ln2;
     double arcmin;
     double lightspeed;
}con;

con constants = {
     3.14159265358979323846, //pi
     9.86960440108935861883, //pisqr
     6.28318530717958647693, //twopi
     0.69314718,
     2.90888208665721580e-4, //arcmin
     299792.458 //speed of light km/s
};

typedef struct {
     double low;
     double medium;
     double high;
     double insane;
}pre;

pre precision= {
  1e-2, //low
  1e-3, //medium
  1e-5, //high
  1e-7 //insane
};

typedef struct {
     double a_min;
     double k_min_mpc;
     double k_max_mpc;
     double k_max_mpc_class;
     double k_min_cH0;
     double k_max_cH0;
     double P_2_s_min;
     double P_2_s_max;
     double xi_via_hankel_theta_min;
     double xi_via_hankel_theta_max;
     double xi_3d_rmin;
     double xi_3d_rmax;
     double M_min;
     double M_max;
}lim;

lim limits = {
//	0.19, //a_min (in order to compute z=4 WFIRST)
    1./(1.+10.), //a_min (z=10, needed for CMB lensing)
	6.667e-6, //k_min_mpc
	1.e3, //k_max_mpc
	50., //k_max_mpc_class
	2.e-2, //k_min_cH0
	3.e+6,  //k_max_cH0
	0.1,//P_2_s_min
	1.0e5,//P_2_s_max
	3.0e-7,//xi_via_hankel_theta_min
	0.12, //xi_via_hankel_theta_max
    1.e-5,//xi_3d_rmin
    1.e+0,//xi_3d_rmax
	1.e+6, //M_min
	1.e+17,//M_max
};

typedef struct {
     int N_a ;
     int N_k_lin;
     int N_k_nlin;
     int N_ell;
     int N_theta;
     int N_thetaH;
     int N_S2;
     int N_DS;
     int N_norm;
     int N_r_3d;
     int N_k_3d;
     int N_a_halo;
}Ntab;

Ntab Ntable = {
100, //N_a
500, //N_k_lin
500, //N_k_nlin
200, //N_ell
200, //N_theta
2048, //N_theta for Hankel
1000, //N_S2
1000, //N_DS
50, //N_norm
50, //N_r_3d
25, //N_k_3d
20, //N_a_halo
};

struct cos{
     int ORDER ;
     double vt_max;
     double vt_min;
     double vt_bin_max;
     double vt_bin_min;
     int ni ;
     int nj ;
};

void SVD_inversion(gsl_matrix *cov, gsl_matrix *inverseSVD,int Nmatrix)
{
    int i,j;
    gsl_matrix *V=gsl_matrix_calloc(Nmatrix,Nmatrix);
    gsl_matrix *U=gsl_matrix_calloc(Nmatrix,Nmatrix);
    gsl_vector *S=gsl_vector_calloc(Nmatrix);
    gsl_vector *work=gsl_vector_calloc(Nmatrix);

    gsl_matrix_memcpy(U,cov);

    gsl_linalg_SV_decomp(U,V,S,work);

    for (i=0;i<Nmatrix;i++){
        gsl_vector *b=gsl_vector_calloc(Nmatrix);
        gsl_vector *x=gsl_vector_calloc(Nmatrix);
        gsl_vector_set(b,i,1.0);
        gsl_linalg_SV_solve(U,V,S,b,x);
        for (j=0;j<Nmatrix;j++){
            gsl_matrix_set(inverseSVD,i,j,gsl_vector_get(x,j));
        }
        gsl_vector_free(b);
        gsl_vector_free(x);
    }
    gsl_matrix_free(V);
    gsl_matrix_free(U);
    gsl_vector_free(S);
    gsl_vector_free(work);
}

void invert_matrix_colesky(gsl_matrix *A)
{
  gsl_linalg_cholesky_decomp (A);  // Adummy will be overwritten  */
  gsl_linalg_cholesky_invert (A);
}

double int_gsl_integrate_insane_precision(double (*func)(double, void*),void *arg,double a, double b, double *error, int niter)
{
  double res, err;
  gsl_integration_cquad_workspace *w = gsl_integration_cquad_workspace_alloc(niter);
  gsl_function F;
  F.function = func;
  F.params  = arg;
 gsl_integration_cquad(&F,a,b,0,precision.insane,w,&res,&err,0);
  if(NULL!=error)
    *error=err;
  gsl_integration_cquad_workspace_free(w);
  return res;
}

double int_gsl_integrate_high_precision(double (*func)(double, void*),void *arg,double a, double b, double *error, int niter)
{
  double res, err;
  gsl_integration_cquad_workspace *w = gsl_integration_cquad_workspace_alloc(niter);
  gsl_function F;
  F.function = func;
  F.params  = arg;
 gsl_integration_cquad(&F,a,b,0,precision.high,w,&res,&err,0);
  if(NULL!=error)
    *error=err;
  gsl_integration_cquad_workspace_free(w);
  return res;
}

double int_gsl_integrate_medium_precision(double (*func)(double, void*),void *arg,double a, double b, double *error, int niter)
{
  double res, err;
  gsl_integration_cquad_workspace *w = gsl_integration_cquad_workspace_alloc(niter);
  gsl_function F;
  F.function = func;
  F.params  = arg;
 gsl_integration_cquad(&F,a,b,0,precision.medium,w,&res,&err,0);
  if(NULL!=error)
    *error=err;
  gsl_integration_cquad_workspace_free(w);
  return res;
}

double int_gsl_integrate_low_precision(double (*func)(double, void*),void *arg,double a, double b, double *error, int niter)
{
  double res, err;
  gsl_integration_cquad_workspace *wcrude = gsl_integration_cquad_workspace_alloc(niter);
  gsl_function F;
  F.function = func;
  F.params  = arg;
 gsl_integration_cquad(&F,a,b,0,precision.low,wcrude,&res,&err,0);
  if(NULL!=error)
    *error=err;
  gsl_integration_cquad_workspace_free(wcrude);
  return res;
}

double int_gsl_integrate_cov_precision(double (*func)(double, void*),void *arg,double a, double b, double *error, int niter)
{
  double res, err;
  gsl_integration_cquad_workspace *wcrude = gsl_integration_cquad_workspace_alloc(niter);
  gsl_function F;
  F.function = func;
  F.params  = arg;
 gsl_integration_cquad(&F,a,b,1e-16,precision.high,wcrude,&res,&err,0);
  if(NULL!=error)
    *error=err;
  gsl_integration_cquad_workspace_free(wcrude);
  return res;
}

int line_count(char *filename)
{
  FILE *n;
  n = fopen(filename,"r");
  if (!n){printf("line_count: %s not found!\nEXIT!\n",filename);exit(1);}
  int ch =0, prev =0,number_of_lines = 0;

  do
    {
    prev = ch;
    ch = fgetc(n);
    if(ch == '\n')
    	number_of_lines++;
    } while (ch != EOF);
  fclose(n);
  // last line might not end with \n, but if previous character does, last line is empty
  if(ch != '\n' && prev !='\n' && number_of_lines != 0)
    number_of_lines++;
  return number_of_lines;
}

double **create_double_matrix(long nrl, long nrh, long ncl, long nch)
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
	double **m;

	/* allocate pointers to rows */
	m=(double **) calloc(nrow+NR_END,sizeof(double*));
	if (!m) error("allocation failure 1 in create_double_matrix()");
	m += NR_END;
	m -= nrl;

	/* allocate rows and set pointers to them */
	m[nrl]=(double *) calloc(nrow*ncol+NR_END,sizeof(double));
	if (!m[nrl]) error("allocation failure 2 in create_double_matrix()");
	m[nrl] += NR_END;
	m[nrl] -= ncl;

	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

	/* return pointer to array of pointers to rows */
	return m;
}

void free_double_matrix(double **m, long nrl, long nrh, long ncl, long nch)
		/* free a double matrix allocated by create_double_matrix() */
{
	free((FREE_ARG) (m[nrl]+ncl-NR_END));
	free((FREE_ARG) (m+nrl-NR_END));
}

double *create_double_vector(long nl, long nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
	double *v;

	v=(double *) calloc(nh-nl+1+NR_END,sizeof(double));
	if (!v) error("allocation failure in double vector()");
	return v-nl+NR_END;
}

int *int_vector(long nl, long nh)
			/* allocate a int vector with subscript range v[nl..nh] */
{
	int *v;

	v=(int *)calloc(nh-nl+1+NR_END, sizeof(int));
	if (!v) error("allocation failure in int vector()");
	return v-nl+NR_END;
}

long *long_vector(long nl, long nh)
			/* allocate a int vector with subscript range v[nl..nh] */
{
	long *v;

	v=(long *)calloc(nh-nl+1+NR_END,sizeof(long));
	if (!v) error("allocation failure in int vector()");
	return v-nl+NR_END;
}

void free_double_vector(double *v, long nl, long nh)
		/* free a double vector allocated with vector() */
{
	free((FREE_ARG) (v+nl-NR_END));
}

void error(char *s)
{
	printf("error:%s\n ",s);
	exit(1);
}

/* ============================================================ *
 * Interpolates f at the value x, where f is a double[n] array,	*
 * representing a function between a and b, stepwidth dx.	*
 * 'lower' and 'upper' are powers of a logarithmic power law	*
 * extrapolation. If no	extrapolation desired, set these to 0	*
 * ============================================================ */
double interpol(double *f, int n, double a, double b, double dx, double x, double lower, double upper)
{
	double r;
	int  i;
	if (x < a) {
		if (lower==0.) {
			//error("value too small in interpol");
			return 0.0;
		}
		return f[0] + lower*(x - a);
	}
	r = (x - a)/dx;
	i = (int)(floor(r));
	if (i+1 >= n) {
		if (upper==0.0) {
			if (i+1==n) {
				return f[i];  /* constant extrapolation */
			} else {
				//error("value too big in interpol");
				return 0.0;
			}
		} else {
			return f[n-1] + upper*(x-b); /* linear extrapolation */
		}
	} else {
		return (r - i)*(f[i+1] - f[i]) + f[i]; /* interpolation */
	}
}
double interpol_fitslope(double *f, int n, double a, double b, double dx, double x, double lower)
{
	double r;
	int  i,fitrange;
	if (x < a) {
		if (lower==0.) {
			//error("value too small in interpol");
			return 0.0;
		}
		return f[0] + lower*(x - a);
	}
	r = (x - a)/dx;
	i = (int)(floor(r));
	if (i+1 >= n) {
		if (n > 50){fitrange =5;}
		else{fitrange = (int)floor(n/10);}
		double upper = (f[n-1] - f[n-1-fitrange])/(dx*fitrange);
		return f[n-1] + upper*(x-b); /* linear extrapolation */

	} else {
		return (r - i)*(f[i+1] - f[i]) + f[i]; /* interpolation */
	}
}

/* ============================================================ *
 * like interpol, but f beeing a 2d-function			*
 * 'lower' and 'upper' are the powers of a power law extra-	*
 * polation in the second argument				*
 * ============================================================ */
double interpol2d(double **f,
		      int nx, double ax, double bx, double dx, double x,
		      int ny, double ay, double by, double dy, double y,
		      double lower, double upper)
{
	double t, dt, s, ds;
	int i, j;
	if (x < ax) {
		return 0.;
//		error("value too small in interpol2d");
	}

	if (x > bx) {
		return 0.;//
//		printf("%le %le\n",x,bx);
	//	error("value too big in interpol2d");
	}
	t = (x - ax)/dx;
	i = (int)(floor(t));
	dt = t - i;
	if (y < ay) {
		return ((1.-dt)*f[i][0] + dt*f[i+1][0]) + (y-ay)*lower;
	} else if (y > by) {
		return ((1.-dt)*f[i][ny-1] + dt*f[i+1][ny-1]) + (y-by)*upper;
	}
	s = (y - ay)/dy;
	j = (int)(floor(s));
	ds = s - j;
	if ((i+1==nx)&&(j+1==ny)) {
	  //printf("%d %d\n",i+1,j+1);
	  return (1.-dt)*(1.-ds)*f[i][j];

	}
	if (i+1==nx){
	  //printf("%d %d\n",i+1,j+1);
	  return (1.-dt)*(1.-ds)*f[i][j]+ (1.-dt)*ds*f[i][j+1];
	}
	if (j+1==ny){
	  //printf("%d %d\n",i+1,j+1);
	  return (1.-dt)*(1.-ds)*f[i][j]+ dt*(1.-ds)*f[i+1][j];
	}
	return (1.-dt)*(1.-ds)*f[i][j] +(1.-dt)*ds*f[i][j+1] + dt*(1.-ds)*f[i+1][j] + dt*ds*f[i+1][j+1];
}
double interpol2d_fitslope(double **f,
		      int nx, double ax, double bx, double dx, double x,
		      int ny, double ay, double by, double dy, double y,
		      double lower)
{
	double t, dt, s, ds, upper;
	int i, j, fitrange;
	if (x < ax) {
		return 0.;
//		error("value too small in interpol2d");
	}

	if (x > bx) {
		return 0.;
//	printf("%le %le\n",x,bx);
//	  error("value too big in interpol2d");
	}
	t = (x - ax)/dx;
	i = (int)(floor(t));
	dt = t - i;
	if (y < ay) {
		return ((1.-dt)*f[i][0] + dt*f[i+1][0]) + (y-ay)*lower;
	} else if (y > by) {
		if (ny > 25){fitrange =5;}
		else{fitrange = (int)floor(ny/5);}
		upper = ((1.-dt)*(f[i][ny-1] - f[i][ny-1-fitrange])+dt*(f[i+1][ny-1] - f[i+1][ny-1-fitrange]))/(dy*fitrange);
		return ((1.-dt)*f[i][ny-1] + dt*f[i+1][ny-1]) + (y-by)*upper;
	}
	s = (y - ay)/dy;
	j = (int)(floor(s));
	ds = s - j;
	if ((i+1==nx)&&(j+1==ny)) {
	  //printf("%d %d\n",i+1,j+1);
	  return (1.-dt)*(1.-ds)*f[i][j];

	}
	if (i+1==nx){
	  //printf("%d %d\n",i+1,j+1);
	  return (1.-dt)*(1.-ds)*f[i][j]+ (1.-dt)*ds*f[i][j+1];
	}
	if (j+1==ny){
	  //printf("%d %d\n",i+1,j+1);
	  return (1.-dt)*(1.-ds)*f[i][j]+ dt*(1.-ds)*f[i+1][j];
	}
	return (1.-dt)*(1.-ds)*f[i][j] +(1.-dt)*ds*f[i][j+1] + dt*(1.-ds)*f[i+1][j] + dt*ds*f[i+1][j+1];
}

void cdgamma(fftw_complex x, fftw_complex *res)
{
	double  xr, xi, wr, wi, ur, ui, vr, vi, yr, yi, t;

	xr = (double) x[0];
	xi = (double) x[1];

	if (xr<0) {
		wr = 1 - xr;
		wi = -xi;
	} else {
		wr = xr;
		wi = xi;
	}

	ur = wr + 6.00009857740312429;
	vr = ur * (wr + 4.99999857982434025) - wi * wi;
	vi = wi * (wr + 4.99999857982434025) + ur * wi;
	yr = ur * 13.2280130755055088 + vr * 66.2756400966213521 +
			0.293729529320536228;
	yi = wi * 13.2280130755055088 + vi * 66.2756400966213521;
	ur = vr * (wr + 4.00000003016801681) - vi * wi;
	ui = vi * (wr + 4.00000003016801681) + vr * wi;
	vr = ur * (wr + 2.99999999944915534) - ui * wi;
	vi = ui * (wr + 2.99999999944915534) + ur * wi;
	yr += ur * 91.1395751189899762 + vr * 47.3821439163096063;
	yi += ui * 91.1395751189899762 + vi * 47.3821439163096063;
	ur = vr * (wr + 2.00000000000603851) - vi * wi;
	ui = vi * (wr + 2.00000000000603851) + vr * wi;
	vr = ur * (wr + 0.999999999999975753) - ui * wi;
	vi = ui * (wr + 0.999999999999975753) + ur * wi;
	yr += ur * 10.5400280458730808 + vr;
	yi += ui * 10.5400280458730808 + vi;
	ur = vr * wr - vi * wi;
	ui = vi * wr + vr * wi;
	t = ur * ur + ui * ui;
	vr = yr * ur + yi * ui + t * 0.0327673720261526849;
	vi = yi * ur - yr * ui;
	yr = wr + 7.31790632447016203;
	ur = log(yr * yr + wi * wi) * 0.5 - 1;
	ui = atan2(wi, yr);
	yr = exp(ur * (wr - 0.5) - ui * wi - 3.48064577727581257) / t;
	yi = ui * (wr - 0.5) + ur * wi;
	ur = yr * cos(yi);
	ui = yr * sin(yi);
	yr = ur * vr - ui * vi;
	yi = ui * vr + ur * vi;
	if (xr<0) {
		wr = xr * 3.14159265358979324;
		wi = exp(xi * 3.14159265358979324);
		vi = 1 / wi;
		ur = (vi + wi) * sin(wr);
		ui = (vi - wi) * cos(wr);
		vr = ur * yr + ui * yi;
		vi = ui * yr - ur * yi;
		ur = 6.2831853071795862 / (vr * vr + vi * vi);
		yr = ur * vr;
		yi = ur * vi;
	}

	(*res)[0]=yr; (*res)[1]=yi;
}

ext: c

max_stars_repo_name: joezuntz/CosmoCov

max_issues_repo_name: joezuntz/CosmoCov

max_forks_repo_name: joezuntz/CosmoCov


text:
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* Set of C routines to use FFTW1.3 library
* 2D- FFT's for any size
*
*
* JLP
* Version 12/12/2006 (from my BORLAND version)
*--------------------------------------------------------*/
#include <stdio.h>
#include <math.h>
#include "jlp_fftw.h"
#include "jlp_num_rename.h"

/* JLP2001: To be compatible with "matlab"
* I should not normalize by sqrt(nx*ny): */
/*
#define NORMALIZE_SQRT
*/

/* Content of jlp_fftw.h:
#include <fftw.h>

int FFTW_1D_Y_FLT(float *re, float *im, INT4 *nx, INT4 *ny, INT4 *idim,
                   INT4 *direct);
int fftw_1D_Y_float(float *re, float *im, int nx, int ny, int direct);
int FFTW_1D_Y_DBLE(double *re, double *im, INT4 *nx, INT4 *ny, INT4 *idim,
             INT4 *direct);
int fftw_1D_Y_double(double *re, double *im, int nx, int ny, int direct);
int FFTW_2D_DBLE(double *re, double *im, int *nx, int *ny, int *direct);
int fftw_2D_double(double *re, double *im, int nx, int ny, int direct);
int FFTW_2D_FLT(float *re, float *im, INT4 *nx, INT4 *ny, INT4 *idim,
                  INT4 *kod);
int fftw_2D_float(float *re, float *im, int nx, int ny, int direct);
int fftw_setup(int nx, int ny);
int FFTW_SETUP(int *nx, int *ny);
int fftw_fast(FFTW_COMPLEX *image, int nx, int ny, int direct);
int fftw_shutdown();
*/

/* Static variables: */
static fftwnd_plan plan_fwd, plan_bkwd;

/*
#define MAIN_TEST
*/

#ifdef MAIN_TEST
main()
{
register int i;
int nx = 9, ny = 8;
double re[128], im[128];
FFTW_COMPLEX* image;
char s[80];

for(i = 0; i < nx * ny; i++)
  {
  re[i] = i;
  im[i] = -i;
  }

fftw_2D_double(re, im, nx, ny, 1);
fftw_2D_double(re, im, nx, ny, -1);

for(i = 0; i < nx * ny; i++)
  {
  printf(" i=%d re=%f im=%f \n", i, re[i], im[i]);
  }

printf(" Now going to fast fft...");
gets(s);

nx = 388; ny = 128;
image = (FFTW_COMPLEX *) malloc(nx * ny * sizeof(FFTW_COMPLEX));
for(i = 0; i < nx * ny; i++)
  {
  c_re(image[i]) = i;
  c_im(image[i]) = -i;
  }
  printf("OK, I go on, with nx=%d ny=%d",nx,ny);

fftw_setup(nx, ny);
fftw_fast(image, nx, ny, 1);
fftw_fast(image, nx, ny, -1);
fftw_shutdown();
for(i = 0; i < 200; i++)
  {
  printf(" i=%d re=%f im=%f \n", i, c_re(image[i]), c_im(image[i]));
  }

gets(s);
free(image);
}
#endif
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* Interface with Fortran routines (with arrays of first dimension = idim)
* kod=1, or 2 direct
* kod=-1, or -2 inverse
* kod=-2, or 2: power spectrum, and phase
*--------------------------------------------------------------*/
int FFTW_2D_FLT(float *re, float *im, INT4 *nx, INT4 *ny, INT4 *idim,
                  INT4 *kod)
{
int nx1, ny1, direct, status;
if(*idim != *nx)
 {printf("FFTW_2D/Fatal error idim=%d while nx=%d \n", *idim, *nx);
  exit(-1);
 }
if(*kod != 1 && *kod != -1)
 {printf("FFTW_2D/Fatal error kod=%d (option not yet implemented)\n", *kod);
  exit(-1);
 }
nx1 = *nx; ny1 = *ny; direct = *kod;
status = fftw_2D_float(re, im, nx1, ny1, direct);
return(status);
}
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* Interface with Fortran routines
* kod=1, or 2 direct
* kod=-1, or -2 inverse
* kod=-2, or 2: power spectrum, and phase
* --------------------------------------------------------------*/
int FFTW_1D_Y_DBLE(double *re, double *im, INT4 *nx, INT4 *ny,
                   INT4 *idim, INT4 *kod)
{
int nx1, ny1, direct, status;
if(*idim != *nx)
 {printf("FFTW_1D_Y/Fatal error idim=%d while nx=%d \n", *idim, *nx);
  exit(-1);
 }
if(*kod != 1 && *kod != -1)
 {printf("FFTW_1D_Y/Fatal error kod=%d (option not yet implemented)\n", *kod);
  exit(-1);
 }
nx1 = *nx; ny1 = *ny; direct = *kod;
printf("FFTW_1D_Y/ nx1=%d ny1=%d direct=%d\n",nx1,ny1,direct);
status = fftw_1D_Y_double(re, im, nx1, ny1, direct);
return(status);
}
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* Interface with Fortran routines
* kod=1, or 2 direct
* kod=-1, or -2 inverse
* kod=-2, or 2: power spectrum, and phase
* --------------------------------------------------------------*/
int FFTW_1D_Y_FLT(float *re, float *im, INT4 *nx, INT4 *ny,
                   INT4 *idim, INT4 *kod)
{
int nx1, ny1, direct, status;

#ifndef TOTO_
double sum;
register int i;
#endif

if(*idim != *nx)
 {printf("FFTW_1D_Y_FLT/Fatal error idim=%d while nx=%d \n", *idim, *nx);
  exit(-1);
 }
if(*kod != 1 && *kod != -1)
 {printf("FFTW_1D_Y_FLT/Fatal error kod=%d (option not yet implemented)\n", *kod);
  exit(-1);
 }
nx1 = *nx; ny1 = *ny; direct = *kod;
#ifdef DEBUG
printf("FFTW_1D_Y_FLT/ nx1=%d ny1=%d direct=%d (fftw)\n",nx1,ny1,direct);
#endif

#ifndef TOTO_
sum =0;
for(i=0; i < ny1; i++) sum += re[nx1/2 + i * nx1];
printf(" Sum of central line = %e and sum/sqrt(ny)=%e \n",
        sum,sum/sqrt((double)ny1));
#endif

status = fftw_1D_Y_float(re, im, nx1, ny1, direct);

/* JLP99: to check that FFT is correct: */
/* Not recentred yet !! */
#ifndef TOTO_
for(i=0; i <= 2; i++)
printf("re,im [ixc,%d]: %e %e \n",i,
        re[nx1/2 + i*nx1],im[nx1/2 + i*nx1]);
#endif

return(status);
}
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* 2D FFT routine for double precision arrays
* direct = 1  : direct or forward
* direct = -1 : reverse or backward
* (FORTRAN version)
* --------------------------------------------------------------*/
int FFTW_2D_DBLE(double *re, double *im, int *nx, int *ny, int *direct)
{
int status;
status = fftw_2D_double(re, im, *nx, *ny, *direct);
return(status);
}
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* 2D FFT routine for double precision arrays
* direct = 1  : direct or forward
* direct = -1 : reverse or backward
*
* --------------------------------------------------------------*/
int fftw_2D_double(double *re, double *im, int nx, int ny, int direct)
{
register int i;
double norm;
int isize;
fftwnd_plan plan;
fftw_direction dir;
FFTW_COMPLEX* in_out;

isize = nx * ny;
in_out = (FFTW_COMPLEX *) malloc(isize * sizeof(FFTW_COMPLEX));

/* Transfer to complex array structure: */
for(i = 0; i < nx * ny; i++)
  {
   c_re(in_out[i]) = re[i];
   c_im(in_out[i]) = im[i];
  }

/* direct = 1  : forward
   direct = -1 : backward
*/
dir = (direct == 1) ? FFTW_FORWARD : FFTW_BACKWARD;
/* Warning: inversion nx-ny here: */
plan = fftw2d_create_plan(ny, nx, dir,
                              FFTW_ESTIMATE | FFTW_IN_PLACE);
if(plan == NULL)
{printf("FFTW_2D_DBLE: fatal error creating plan\n"); exit(-1);}

/* Compute the FFT: */
   fftwnd(plan, 1, in_out, 1, 0, 0, 0, 0);

/* Transfer back to original arrays (and normalize by sqrt(nx * ny): */
#ifdef NORMALIZE_SQRT
norm = nx * ny; norm = sqrt(norm);
for(i = 0; i < nx * ny; i++)
  {
  re[i] = c_re(in_out[i]) /norm;
  im[i] = c_im(in_out[i]) /norm;
  }
#else
if(direct == 1)
  {
  for(i = 0; i < nx * ny; i++)
    {
    re[i] = c_re(in_out[i]);
    im[i] = c_im(in_out[i]);
    }
  }
else
  {
  norm = nx * ny;
  for(i = 0; i < nx * ny; i++)
    {
    re[i] = c_re(in_out[i]) /norm;
    im[i] = c_im(in_out[i]) /norm;
    }
  }
#endif

free(in_out);

/* Delete plan: */
  fftwnd_destroy_plan(plan);

return(0);
}
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* 2D FFT routine for single precision arrays (for which idim=nx)
* (well suited to C arrays since nx=idim...)
* direct = 1  : direct or forward
* direct = -1 : reverse or backward
* --------------------------------------------------------------*/
int fftw_2D_float(float *re, float *im, int nx, int ny, int direct)
{
register int i;
double norm;
int isize;
fftwnd_plan plan;
fftw_direction dir;
FFTW_COMPLEX* in_out;

isize = nx * ny;
in_out = (FFTW_COMPLEX *) malloc(isize * sizeof(FFTW_COMPLEX));

/* Transfer to complex array structure: */
for(i = 0; i < nx * ny; i++)
  {
   c_re(in_out[i]) = re[i];
   c_im(in_out[i]) = im[i];
  }

/* direct = 1  : forward
   direct = -1 : backward
*/
dir = (direct == 1) ? FFTW_FORWARD : FFTW_BACKWARD;
/* Warning: inversion nx-ny here: */
plan = fftw2d_create_plan(ny, nx, dir,
                              FFTW_ESTIMATE | FFTW_IN_PLACE);
if(plan == NULL)
{printf("fftw_2D_float: fatal error creating plan\n"); exit(-1);}

/* Compute the FFT: */
   fftwnd(plan, 1, in_out, 1, 0, 0, 0, 0);

/* Transfer back to original arrays (and normalize by sqrt(nx * ny): */
#ifdef NORMALIZE_SQRT
norm = nx * ny; norm = sqrt(norm);
for(i = 0; i < nx * ny; i++)
  {
  re[i] = c_re(in_out[i]) /norm;
  im[i] = c_im(in_out[i]) /norm;
  }
#else
if(direct == 1)
  {
  for(i = 0; i < nx * ny; i++)
    {
    re[i] = c_re(in_out[i]);
    im[i] = c_im(in_out[i]);
    }
  }
else
  {
  norm = nx * ny;
  for(i = 0; i < nx * ny; i++)
    {
    re[i] = c_re(in_out[i]) /norm;
    im[i] = c_im(in_out[i]) /norm;
    }
  }
#endif

free(in_out);

/* Delete plan: */
  fftwnd_destroy_plan(plan);

return(0);
}
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* 1D FFT routine for single precision arrays
* along the columns (for spectroscopic mode)
* direct = 1  : direct or forward
* direct = -1 : reverse or backward
*
* --------------------------------------------------------------*/
int fftw_1D_Y_float(float *re, float *im, int nx, int ny, int direct)
{
register int i, j;
double norm;
fftw_plan plan;
fftw_direction dir;
FFTW_COMPLEX* in_out;

in_out = (FFTW_COMPLEX *) malloc(ny * sizeof(FFTW_COMPLEX));
if(in_out == NULL)
{printf("fftw_1D_Y_float: fatal error allocating memory\n"); exit(-1);}

/* direct = 1  : forward
   direct = -1 : backward
*/
dir = (direct == 1) ? FFTW_FORWARD : FFTW_BACKWARD;
plan = fftw_create_plan(ny, dir, FFTW_ESTIMATE | FFTW_IN_PLACE);
if(plan == NULL)
{printf("fftw_1D_Y_float: fatal error creating plan\n"); exit(-1);}

for(i = 0; i < nx; i++)
{
/* Transfer to complex array structure: */
  for(j = 0; j < ny; j++)
  {
   c_re(in_out[j]) = re[i + j * nx];
   c_im(in_out[j]) = im[i + j * nx];
  }

/* Compute the FFT: */
   fftw(plan, 1, in_out, 1, 0, 0, 0, 0);

/* Transfer back to original arrays (and normalize by sqrt(ny): */
#ifdef NORMALIZE_SQRT
norm = sqrt((double)ny);
for(j = 0; j < ny; j++)
  {
  re[i + j * nx] = c_re(in_out[j]) /norm;
  im[i + j * nx] = c_im(in_out[j]) /norm;
  }
#else
if(direct == 1)
  {
  for(j = 0; j < ny; j++)
    {
    re[i + j * nx] = c_re(in_out[j]);
    im[i + j * nx] = c_im(in_out[j]);
    }
  }
else
  {
  norm = (double)ny;
  for(j = 0; j < ny; j++)
    {
    re[i + j * nx] = c_re(in_out[j]) /norm;
    im[i + j * nx] = c_im(in_out[j]) /norm;
    }
  }
#endif
/* End of loop on i (from 0 to nx-1) */
}

free(in_out);

/* Delete plan: */
  fftw_destroy_plan(plan);

return(0);
}
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* 1D FFT routine for double precision arrays
* along the columns (for spectroscopic mode)
* direct = 1  : direct or forward
* direct = -1 : reverse or backward
*
* --------------------------------------------------------------*/
int fftw_1D_Y_double(double *re, double *im, int nx, int ny, int direct)
{
register int i, j;
double norm;
fftw_plan plan;
fftw_direction dir;
FFTW_COMPLEX* in_out;

in_out = (FFTW_COMPLEX *) malloc(ny * sizeof(FFTW_COMPLEX));

/* direct = 1  : forward
   direct = -1 : backward
*/
dir = (direct == 1) ? FFTW_FORWARD : FFTW_BACKWARD;
plan = fftw_create_plan(ny, dir, FFTW_ESTIMATE | FFTW_IN_PLACE);
if(plan == NULL)
{printf("fftw_1D_Y: fatal error creating plan\n"); exit(-1);}

for(i = 0; i < nx; i++)
{
/* Transfer to complex array structure: */
  for(j = 0; j < ny; j++)
  {
   c_re(in_out[j]) = re[i + j * nx];
   c_im(in_out[j]) = im[i + j * nx];
  }

/* Compute the FFT: */
   fftw(plan, 1, in_out, 1, 0, 0, 0, 0);

/* Transfer back to original arrays (and normalize by sqrt(ny): */
#ifdef NORMALIZE_SQRT
norm = sqrt((double)ny);
for(j = 0; j < ny; j++)
  {
  re[i + j * nx] = c_re(in_out[j]) /norm;
  im[i + j * nx] = c_im(in_out[j]) /norm;
  }
#else
if(direct == 1)
  {
  for(j = 0; j < ny; j++)
    {
    re[i + j * nx] = c_re(in_out[j]);
    im[i + j * nx] = c_im(in_out[j]);
    }
  }
else
  {
  norm = (double)ny;
  for(j = 0; j < ny; j++)
    {
    re[i + j * nx] = c_re(in_out[j]) /norm;
    im[i + j * nx] = c_im(in_out[j]) /norm;
    }
  }
#endif
/* End of loop on i (from 0 to nx-1) */
}

free(in_out);

/* Delete plan: */
  fftw_destroy_plan(plan);

return(0);
}
/****************************************************************
* Fortran interface to fftw_setup
****************************************************************/
int FFTW_FSETUP(int *nx, int *ny)
{
fftw_setup(*nx, *ny);
return(0);
}
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// fftw_setup, to initialize FFTW routines
// Create forward and backward plans
//
// --------------------------------------------------------------*/
int fftw_setup(int nx, int ny)
{
char fwd_name[60], bkwd_name[60];
FILE *fp_wisd;
int ToOutput;

/********************** With "forward"  ***************************/
#ifdef BORLAND
sprintf(fwd_name,"c:\\tc\\fftw13\\jlp\\fwd_%d.wis",nx);
#else
sprintf(fwd_name,"/d/fftw/fwd_%d%d.wis",nx,ny);
#endif

/* Check if wisdom file already exits: */
printf("fftw_setup/Read wisdom file: %s\n",fwd_name);
if((fp_wisd = fopen(fwd_name,"r")) == NULL)
  {
  printf("fftw_setup/Failure to open wisdom file: %s\n",fwd_name);
  ToOutput = 1;
  }
else
  {
  ToOutput = 0;
  if(fftw_import_wisdom_from_file(fp_wisd) == FFTW_FAILURE) ToOutput = 1;
  fclose(fp_wisd);
  }

/* Create  "plan" for fftw (speed is wanted for subsequent FFT's) */
/* Warning: inversion nx-ny here: */
plan_fwd = fftw2d_create_plan(ny, nx, FFTW_FORWARD,
                              FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);

/* Output wisdom file if needed: */
if(ToOutput)
  {
  printf("fftw_setup/Write wisdom file: %s\n",fwd_name);
  if((fp_wisd = fopen(fwd_name,"w")) != NULL)
                                       fftw_export_wisdom_to_file(fp_wisd);
  fclose(fp_wisd);
  }

/********************** With "backward"  ***************************/
#ifdef BORLAND
sprintf(bkwd_name,"c:\\tc\\fftw13\\jlp\\bk_%d%d.wis",nx,ny);
#else
sprintf(bkwd_name,"/d/fftw/bk_%d%d.wis",nx,ny);
#endif

/* Check if wisdom file already exits: */
printf("fftw_setup/Read wisdom file: %s\n",bkwd_name);
if((fp_wisd = fopen(bkwd_name,"r")) == NULL)
  {
  printf("fftw_setup/Failure to open wisdom file: %s\n",bkwd_name);
  ToOutput = 1;
  }
else
  {
  ToOutput = 0;
  if(fftw_import_wisdom_from_file(fp_wisd) == FFTW_FAILURE) ToOutput = 1;
  fclose(fp_wisd);
  }

/* Create  "plan" for fftw (speed is wanted for subsequent FFT's) */
/* Warning: inversion nx-ny here: */
plan_bkwd = fftw2d_create_plan(ny, nx, FFTW_BACKWARD,
                              FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);

/* Output wisdom file if needed: */
if(ToOutput)
  {
  printf("fftw_setup/Write wisdom file: %s\n",bkwd_name);
  if((fp_wisd = fopen(bkwd_name,"w")) != NULL)
                                        fftw_export_wisdom_to_file(fp_wisd);
  fclose(fp_wisd);
  }

return(0);
}
/*************************************************************
*
**************************************************************/
int FFTW_SHUTDOWN()
{
/* Delete plans: */
  fftwnd_destroy_plan(plan_fwd);
  fftwnd_destroy_plan(plan_bkwd);
return(0);
}
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* 2D FFT routine for double precision fftw_complex arrays
* direction = 1  : forward
* direction = -1 : backward
*  Should be called after fftw_setup and before fftw_shutdown
* --------------------------------------------------------------*/
int fftw_fast(FFTW_COMPLEX *image, int nx, int ny, int dir)
{
register int i;
double norm;

/* Compute the FFT: */
if(dir == 1)
   fftwnd(plan_fwd, 1, image, 1, 0, 0, 0, 0);
else
   fftwnd(plan_bkwd, 1, image, 1, 0, 0, 0, 0);

/* Normalize by sqrt(nx * ny): */
#ifdef NORMALIZE_SQRT
norm = nx * ny; norm = sqrt(norm);
for(i = 0; i < nx * ny; i++)
  {
  c_re(image[i]) /= norm;
  c_im(image[i]) /= norm;
  }
#else
if(dir == -1)
  {
  norm = (double)(nx * ny);
  for(i = 0; i < nx * ny; i++)
    {
    c_re(image[i]) /= norm;
    c_im(image[i]) /= norm;
    }
  }
#endif

return(0);
}

max_stars_repo_name: jlprieur/jlplib

max_issues_repo_name: jlprieur/jlplib

max_forks_repo_name: jlprieur/jlplib


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sort.h>

#include "cosmocalc.h"

gsl_spline *cosmocalc_aexpn2comvdist_spline = NULL;
gsl_interp_accel *cosmocalc_aexpn2comvdist_acc = NULL;
gsl_spline *cosmocalc_comvdist2aexpn_spline = NULL;
gsl_interp_accel *cosmocalc_comvdist2aexpn_acc = NULL;
double DH_sqrtok;

/* function for integration using gsl integration */
static double comvdist_integ_funct(double a, void *p)
{
  return 1.0/a/a/hubble_noscale(a);
}

double comvdist_exact(double a)
{
#define WORKSPACE_NUM 100000
#define ABSERR 0.0
#define RELERR 1e-8

  gsl_integration_workspace *workspace;
  gsl_function F;
  double result,abserr;

  workspace = gsl_integration_workspace_alloc((size_t) WORKSPACE_NUM);

  F.function = &comvdist_integ_funct;
  gsl_integration_qag(&F,a,1.0,ABSERR,RELERR,(size_t) WORKSPACE_NUM,GSL_INTEG_GAUSS51,workspace,&result,&abserr);

  gsl_integration_workspace_free(workspace);

#undef ABSERR
#undef RELERR
#undef WORKSPACE_NUM

  return result*DH;
}

/* init function  - some help from Gadget-2 applied here */
void init_cosmocalc_comvdist_table(void)
{
#define WORKSPACE_NUM 100000
#define ABSERR 0.0
#define RELERR 1e-8

  static int initFlag = 1;
  static int currCosmoNum;

  gsl_integration_workspace *workspace;
  gsl_function F;
  long i;
  double result,abserr,afact;
  double comvdist_table[COSMOCALC_COMVDIST_TABLE_LENGTH];
  double aexpn_table[COSMOCALC_COMVDIST_TABLE_LENGTH];
  size_t sindex[COSMOCALC_COMVDIST_TABLE_LENGTH];
  double tmpDouble[COSMOCALC_COMVDIST_TABLE_LENGTH];

  if(initFlag == 1 || currCosmoNum != cosmoData.cosmoNum)
    {
      initFlag = 0;
      currCosmoNum = cosmoData.cosmoNum;

      if(cosmoData.OmegaK != 0.0)
	DH_sqrtok = DH/sqrt(fabs(cosmoData.OmegaK));
      else
	DH_sqrtok = 1.0;

      workspace = gsl_integration_workspace_alloc((size_t) WORKSPACE_NUM);

      for(i=0;i<COSMOCALC_COMVDIST_TABLE_LENGTH-1;++i)
	{
	  afact = (AEXPN_MAX - AEXPN_MIN_DIST)/(COSMOCALC_COMVDIST_TABLE_LENGTH-1.0)*((double) i) + AEXPN_MIN_DIST;
	  F.function = &comvdist_integ_funct;
	  F.params = &(cosmoData.OmegaM);
	  gsl_integration_qag(&F,afact,1.0,ABSERR,RELERR,(size_t) WORKSPACE_NUM,GSL_INTEG_GAUSS51,workspace,&result,&abserr);
	  aexpn_table[i] = afact;
	  comvdist_table[i] = result*DH;
	}
      aexpn_table[i] = 1.0;
      comvdist_table[i] = 0.0;

      gsl_integration_workspace_free(workspace);
#undef ABSERR
#undef RELERR
#undef WORKSPACE_NUM

      //init the spline and accelerators
      if(cosmocalc_aexpn2comvdist_spline != NULL)
	gsl_spline_free(cosmocalc_aexpn2comvdist_spline);
      cosmocalc_aexpn2comvdist_spline = gsl_spline_alloc(GSL_SPLINE_TYPE,(size_t) (COSMOCALC_COMVDIST_TABLE_LENGTH));
      gsl_spline_init(cosmocalc_aexpn2comvdist_spline,aexpn_table,comvdist_table,(size_t) (COSMOCALC_COMVDIST_TABLE_LENGTH));
      if(cosmocalc_aexpn2comvdist_acc != NULL)
	gsl_interp_accel_reset(cosmocalc_aexpn2comvdist_acc);
      else
	cosmocalc_aexpn2comvdist_acc = gsl_interp_accel_alloc();

      //sort the distances properly
      gsl_sort_index(sindex,comvdist_table,(size_t) 1,(size_t) (COSMOCALC_COMVDIST_TABLE_LENGTH));
      for(i=0;i<COSMOCALC_COMVDIST_TABLE_LENGTH;++i)
	tmpDouble[i] = comvdist_table[sindex[i]];
      for(i=0;i<COSMOCALC_COMVDIST_TABLE_LENGTH;++i)
	comvdist_table[i] = tmpDouble[i];
      for(i=0;i<COSMOCALC_COMVDIST_TABLE_LENGTH;++i)
	tmpDouble[i] = aexpn_table[sindex[i]];
      for(i=0;i<COSMOCALC_COMVDIST_TABLE_LENGTH;++i)
	aexpn_table[i] = tmpDouble[i];

      if(cosmocalc_comvdist2aexpn_spline != NULL)
	gsl_spline_free(cosmocalc_comvdist2aexpn_spline);
      cosmocalc_comvdist2aexpn_spline = gsl_spline_alloc(GSL_SPLINE_TYPE,(size_t) (COSMOCALC_COMVDIST_TABLE_LENGTH));
      gsl_spline_init(cosmocalc_comvdist2aexpn_spline,comvdist_table,aexpn_table,(size_t) (COSMOCALC_COMVDIST_TABLE_LENGTH));
      if(cosmocalc_comvdist2aexpn_acc != NULL)
	gsl_interp_accel_reset(cosmocalc_comvdist2aexpn_acc);
      else
	cosmocalc_comvdist2aexpn_acc = gsl_interp_accel_alloc();
    }
}

double acomvdist(double dist)
{
  static int initFlag = 1;
  static int currCosmoNum;

  if(initFlag == 1 || currCosmoNum != cosmoData.cosmoNum)
    {
      initFlag = 0;
      currCosmoNum = cosmoData.cosmoNum;
      init_cosmocalc_comvdist_table();
    }

  return gsl_spline_eval(cosmocalc_comvdist2aexpn_spline,dist,cosmocalc_comvdist2aexpn_acc);
}

double comvdist(double a)
{
  static int initFlag = 1;
  static int currCosmoNum;

  if(initFlag == 1 || currCosmoNum != cosmoData.cosmoNum)
    {
      initFlag = 0;
      currCosmoNum = cosmoData.cosmoNum;
      init_cosmocalc_comvdist_table();
    }

  return gsl_spline_eval(cosmocalc_aexpn2comvdist_spline,a,cosmocalc_aexpn2comvdist_acc);
}

//NOTE: you must call comvdist *FIRST* to init DH_sqrtok!
double angdist(double a)
{
  double cd;
  cd = comvdist(a);

  if(cosmoData.OmegaK > 0.0)
    return DH_sqrtok*sinh(cd/DH_sqrtok)*a;
  else if(cosmoData.OmegaK < 0)
    return DH_sqrtok*sin(cd/DH_sqrtok)*a;
  else
    return cd*a;
}

double lumdist(double a)
{
  double cd;
  cd = comvdist(a);

  if(cosmoData.OmegaK > 0.0)
    return DH_sqrtok*sinh(cd/DH_sqrtok)/a;
  else if(cosmoData.OmegaK < 0)
    return DH_sqrtok*sin(cd/DH_sqrtok)/a;
  else
    return cd/a;
}

double angdistdiff(double amin, double amax)
{
  double cdmin,cdmax;
  assert(amin <= amax);
  cdmin = comvdist(amin);
  cdmax = comvdist(amax);

  if(cosmoData.OmegaK > 0.0)
    return DH_sqrtok*sinh((cdmin-cdmax)/DH_sqrtok)*amin;
  else if(cosmoData.OmegaK < 0)
    return DH_sqrtok*sin((cdmin-cdmax)/DH_sqrtok)*amin;
  else
    return (cdmin-cdmax)*amin;
}

max_stars_repo_name: beckermr/cosmocalc

max_issues_repo_name: beckermr/cosmocalc

max_forks_repo_name: beckermr/cosmocalc


text:
#include <math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>

#define max(A, B) ((A) > (B) ? (A) : (B))
#define min(A, B) ((A) < (B) ? (A) : (B))

double R0(double flock_size, double housing_period, double prob_dam_to_lamb, int maxage, double *params)
{
    int n_age_classes = (maxage+1)*12;
    int c, a, i, y, t;
    int infectious_period, years_infectious, period_infectious_to_lambs;

    int latent_period = lrint(params[5]);
    double beta_field = params[0];
    double beta_housed = params[1];

    double p_mort, p_cull, sum, lambs_per_dam;
    double removal_rate[n_age_classes], infection_rate;
    double cohort_size_age_0;
    double cohort_size[maxage+1][n_age_classes];
    double infecteds[maxage+1][n_age_classes];
    double R0[n_age_classes][n_age_classes];

    // age-specific mortality and cull probabilities from Andrew's data
    double prob_mortality[] = {
        0.01,
        0.012,
        0.027,
        0.041,
        0.019,
        0.021,
        0.034,
    };

    double prob_culled[] = {
        0,
        0.046,
        0.065,
        0.093,
        0.134,
        0.196,
        0.295,
    };

    // create vector of removal rates of a cohort from ages 0 to n_age_classes
    for (a = 0; a < n_age_classes; a++) {
        p_mort = prob_mortality[a/12];
        p_cull = prob_culled[a/12];
        removal_rate[a] = -(log(1.-p_mort)+log(1.-p_cull))/12.;
    }

    // these two ways are equivalent
    // create vector of infection rates within a year starting in April
    infection_rate = ((12-housing_period)*beta_field+housing_period*beta_housed)/12.;
    // for (a = 0; a < 12; a++)
        // if (a < 12-housing_period)
        //     infection_rate[a] = beta_field;
        // else
        //     infection_rate[a] = beta_housed;

    // find intial cohort size in a non-infected flock
    // take initial cohort size as 1
    // find size of cohort as it ages
    cohort_size[0][0] = 1;
    for (a = 1; a < n_age_classes; a++)
        cohort_size[0][a] = cohort_size[0][a-1]*exp(-removal_rate[a]);
    // sum size of all existent cohorts in April
    sum = 0;
    for (y = 0; y <= maxage; y++)
        sum += cohort_size[0][12*y];
    cohort_size_age_0 = flock_size/sum;

    // for cohorts of y years old at first exposure
    // assuming that the infectious ewe is starts to be infectious in April - but
    // this is not necessary, just convenient
    for (y = 0; y <= maxage; y++) {
        // cohort size at age 0
        cohort_size[y][0] = cohort_size_age_0;

        // just removal for first y years before exposure
        for (a = 1; a < 12*y; a++)
            cohort_size[y][a] = cohort_size[y][a-1]*exp(-removal_rate[a]);
        // then removal and infection for rest of lifespan
        for (; a < n_age_classes; a++)
            cohort_size[y][a] = cohort_size[y][a-1]*exp(-removal_rate[a]-infection_rate);

        // calculate number of infecteds
        infecteds[y][0] = 0;
        for (a = 1; a < n_age_classes; a++)
            infecteds[y][a] = infection_rate/(removal_rate[a]+infection_rate)*(cohort_size[y][a-1]-cohort_size[y][a]);
    }

    // set matrix elements to zero
    for (a = 0; a < n_age_classes; a++)
        for (i = 0; i < n_age_classes; i++)
            R0[a][i] = 0;

    // for a ewe infected at age a, sum the number of ewes it infects of age i
    for (a = 0; a < n_age_classes; a++) {
        infectious_period = max(0, n_age_classes-a-latent_period);
        for (t = 0; t < infectious_period; t++) {
            // for all cohorts from past to future
            for (y = maxage; y >= -maxage; y--) {
                // get cohort
                // +ve y are cohorts in the past
                // -ve y are future cohorts (identical to cohort 0)
                c = max(0, y);
                // age of cohort y at time t
                i = 12*y+t;
                if (0 <= i && i < n_age_classes)
                    R0[a][i] += infecteds[c][i];
            }
        }
    }

    // contribution of vertical transmission to an infectious ewe's lambs
    lambs_per_dam = cohort_size_age_0/(flock_size-cohort_size_age_0);

    for (a = 0; a < n_age_classes; a++) {
        // a ewe infected at age a is infectious for max_age-a-1-latent_period months
        infectious_period = max(0, n_age_classes-a-latent_period);
        if (a < 24)
            // only 2 year olds and older give birth
            // so ewes infectious before 24 months cannot infect their lambs
            // until they are 2 years old
            period_infectious_to_lambs = max(0, n_age_classes-24-latent_period);
        else
            period_infectious_to_lambs = infectious_period;
        years_infectious = period_infectious_to_lambs/12;
        R0[a][0] += lambs_per_dam*prob_dam_to_lamb*years_infectious;
    }

    gsl_matrix *A = gsl_matrix_alloc(n_age_classes, n_age_classes);
    gsl_vector_complex *v = gsl_vector_complex_alloc(n_age_classes);
    gsl_eigen_nonsymm_workspace *w = gsl_eigen_nonsymm_alloc(n_age_classes);
    for (a = 0; a < n_age_classes; a++)
        for (i = 0; i < n_age_classes; i++)
            gsl_matrix_set(A, a, i, R0[a][i]);
    gsl_eigen_nonsymm(A, v, w);
    double max_eigenvalue = 0, e;
    for (a = 0; a < n_age_classes; a++) {
        e = GSL_REAL(gsl_vector_complex_get(v, a));
        if (e > max_eigenvalue)
            max_eigenvalue = e;
    }
    gsl_eigen_nonsymm_free(w);
    gsl_vector_complex_free(v);
    gsl_matrix_free(A);

    return max_eigenvalue;
}

max_stars_repo_name: nicksavill/maedi-visna-epidemiology

max_issues_repo_name: nicksavill/maedi-visna-epidemiology

max_forks_repo_name: nicksavill/maedi-visna-epidemiology