-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathrndField.cpp
112 lines (83 loc) · 2.47 KB
/
rndField.cpp
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
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
#include <iostream>
#include "math.h"
#include "healpix_base.h"
#include "healpix_map.h"
#include "planck_rng.h"
#include "pointing.h"
#include "error_handling.h"
#include "TFile.h"
#include "TApplication.h"
#include "TTree.h"
#include "TH2F.h"
using namespace std;
class rndField {
public:
rndField(int);
double *getRaDec();
void genField(int, int);
double* getMap();
private:
Healpix_Base _base;
int _nside;
double* _map;
};
rndField::rndField(int nside) : _nside(nside) {
_base = Healpix_Base(_nside,RING);
}
void rndField::genField(int cell, int values) {
_map = new double[2*values];
pointing p = _base.pix2ang(cell);
double pi = acos(-1.0);
double th_c = p.theta;
double ph_c = p.phi;
double cth_min = cos(th_c-7.0*pi/180.0);
double cth_max = cos(th_c+7.0*pi/180.0);
if (cth_min > cth_max) {
int tmp = cth_max;
cth_max = cth_min;
cth_min = tmp;
}
double ph_min = ph_c-7.0*pi/180.0;
double ph_max = ph_c+7.0*pi/180.00;
planck_rng rng;
int count = 0;
while (count < values) {
double phi = (ph_max-ph_min)*rng.rand_uni()+ph_min;
double theta = acos( (cth_max-cth_min)*rng.rand_uni()+cth_min );
pointing pt(theta,phi);
int pix = _base.ang2pix(pt);
if (pix != cell) {
continue;
}
_map[count] = phi;
_map[values+count+1] = pi/2.0 - theta;
count++;
cout << count << " " << theta*180.0/pi << " - " << phi*180.0/pi << " " << pix << endl;
}
}
double* rndField::getMap() {
return _map;
}
int main(int argc, char* argv[]) {
double pi = acos(-1.0);
planck_rng rng;
// Healpix_Map<double> map = Healpix_Map<double>();
int nsides = 3;
rndField* field = new rndField(nsides);
Healpix_Base base = Healpix_Base(3,RING);
double cth_min = cos(53.0*pi/180.0);
double cth_max = cos(66.0*pi/180.0);
double ph_min = 156.0*pi/180.0;
double ph_max = 170.0*pi/180.00;
int count = 0;
TApplication* rootapp = new TApplication("Example",&argc,argv);
TH2* hdd = new TH2F("hdd", "Random galaxy phi / theta distribution", 100, 53, 66, 100, 156, 170);
int cell = 190;
int nvalues = 10;
field->genField(cell,nvalues);
double* map = field->getMap();
for (int i=0; i<nvalues; i++) {
cout << i << " " << map[i]*180.0/pi << " " << map[nvalues+i+1]*180.0/pi << " " << endl;
}
return 0;
}