mockgallib/src

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
#ifndef HALO_MASS_H
#define HALO_MASS_H 1

#include <gsl/gsl_spline.h>

class HaloMassFoF {
 public:
  HaloMassFoF(const char filename[], const double M0=1.0e14);
  ~HaloMassFoF();
  double mass(const int nfof) const {
    // statistics if poor for large rare haloes
    if(nfof > nfof0_) return m_ * nfof;
  
    return gsl_spline_eval(spline_, nfof, acc_);
  }
 private:
  gsl_spline* spline_;
  gsl_interp_accel* acc_;
  double nfof0_, M0_, m_;
};

class HaloMassFileError {

};

#endif

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
#include <vector>
#include <cstdio>
#include <cassert>
#include "msg.h"
#include "util.h"
#include "halo_mass.h"

using namespace std;

static void read_nfof_mass(const char filename[],
			   vector<double>& v_nfof, vector<double>& v_M);

HaloMassFoF::HaloMassFoF(const char filename[], const double M0) :
  M0_(M0)
{
  // Reads nfof -> M file
  
  vector<double> v_nfof, v_M;

  FILE* fp= fopen(filename, "r");
  if(fp == 0) {
    msg_printf(msg_fatal,
	       "Error: unable to open FoF -> M calibration file %s\n",
	       filename);
    throw HaloMassFileError();
  }

  char buf[128];
  while(fgets(buf, 127, fp)) {
    if(buf[0] == '#')
      continue;

    int nfof;
    double M;

    int ret= sscanf(buf, "%d %le\n", &nfof, &M); assert(ret == 2);
    v_nfof.push_back((double) nfof);
    v_M.push_back(M);

    if(M > M0_) {
      m_ = M/nfof;
      nfof0_= nfof;
      break;
    }
  }

  fclose(fp);

  spline_ = gsl_spline_alloc(gsl_interp_cspline, v_nfof.size());
  acc_ = gsl_interp_accel_alloc();
  gsl_spline_init(spline_, &v_nfof.front(), &v_M.front(), v_nfof.size());

}

HaloMassFoF::~HaloMassFoF()
{
  gsl_spline_free(spline_);
  gsl_interp_accel_free(acc_);
}