FML
Functions
FML::SURVEY Namespace Reference

This namespace contains methods for dealing with surveys. More...

Functions

template<class T , class U >
void EquitorialToCartesianCoordinates (T *galaxies_ra_dec_z, size_t ngalaxies, std::vector< U > &particles_xyz, std::function< double(double)> &hubble_over_c_of_z, std::pair< double, double > &min_max_x, std::pair< double, double > &min_max_y, std::pair< double, double > &min_max_z)
 Take a set of galaxies galaxies_ra_dec_z with (RA,DEC,z) and convert them to cartesian coordinates (x,y,z) stored in particles_xyz. More...
 
template<class T , class U >
void GalaxiesToBox (const T *galaxies_ra_dec_z, size_t ngalaxies, std::vector< U > &particles_xyz, std::function< double(double)> &hubble_over_c_of_z, double &boxsize, std::vector< double > &observer_position, bool shiftPositions, bool scalePositions, bool verbose)
 Transform galaxies with positions defined by RA,DEC,Z and transform these to cartesian positions in [0,1). More...
 
template<class T , class U >
void GalaxiesRandomsToBox (const T *galaxies_ra_dec_z, size_t ngalaxies, const T *randoms_ra_dec_z, size_t nrandoms, std::vector< U > &galaxies_xyz, std::vector< U > &randoms_xyz, std::function< double(double)> &hubble_over_c_of_z, double &boxsize, bool shiftPositions, bool scalePositions, std::vector< double > &observer_position, bool verbose)
 Transform galaxies and randos with positions defined by RA,DEC,Z and transform these to cartesian positions in [0,1. More...
 

Detailed Description

This namespace contains methods for dealing with surveys.

So far its not much here just some conversion from equitorial to cartesian coordinates.

Function Documentation

◆ EquitorialToCartesianCoordinates()

template<class T , class U >
void FML::SURVEY::EquitorialToCartesianCoordinates ( T *  galaxies_ra_dec_z,
size_t  ngalaxies,
std::vector< U > &  particles_xyz,
std::function< double(double)> &  hubble_over_c_of_z,
std::pair< double, double > &  min_max_x,
std::pair< double, double > &  min_max_y,
std::pair< double, double > &  min_max_z 
)

Take a set of galaxies galaxies_ra_dec_z with (RA,DEC,z) and convert them to cartesian coordinates (x,y,z) stored in particles_xyz.

The positions will be in the same units as the 1.0/hubble_over_c_of_z Gives back the min/max of the positions (useful for boxing the catalog)

Template Parameters
TParticle class for the galaxies
UParticle class for the particles we make from the galaxies
Parameters
[in]galaxies_ra_dec_zParticles with RA, DEC and Z.
[in]ngalaxiesNumber of galaxies
[out]particles_xyzVector with galaxies as particles with cartesian coordinates.
[in]hubble_over_c_of_zThis is the function \( H(z)/c \) used to compute the redshift-comobing distance relationship. Postions units is the same as those of \( c/H(z) \) (so e.g. if you want Mpc/h then we need H0 = 100, if you want kpc/s then use H0 = 10^5 and so on.
[out]min_max_xThe min/max values of x-postions
[out]min_max_yThe min/max values of x-postions
[out]min_max_zThe min/max values of x-postions

Definition at line 83 of file GalaxiesToBox.h.

89 {
90
91 // Sanity check on the particle class, must have get_RA, etc. methods
92 static_assert(FML::PARTICLE::has_get_RA<T>());
93 static_assert(FML::PARTICLE::has_get_DEC<T>());
94 static_assert(FML::PARTICLE::has_get_z<T>());
95
96 // Find maximum redshift
97 double z_max = 0.0;
98#ifdef USE_OMP
99#pragma omp parallel for
100#endif
101 for (size_t i = 0; i < ngalaxies; i++) {
102 const double z = FML::PARTICLE::GetRedshift(galaxies_ra_dec_z[i]);
103 z_max = std::max(z, z_max);
104 }
105
106 // Make a z-array
107 const int nz = 10000;
108 std::vector<double> z_arr(nz);
109 std::vector<double> r_arr(nz);
110 for (int i = 0; i < nz; i++) {
111 z_arr[i] = 1.1 * z_max * i / double(nz - 1);
112 }
113
114 // Solve the ODE for the co-moving distance
115 using ODESolver = FML::SOLVERS::ODESOLVER::ODESolver;
116 using ODEFunction = FML::SOLVERS::ODESOLVER::ODEFunction;
117 using DVector = FML::SOLVERS::ODESOLVER::DVector;
118 ODEFunction deriv = [&](double z, [[maybe_unused]] const double * y, double * dydx) {
119 dydx[0] = 1.0 / hubble_over_c_of_z(z);
120 return GSL_SUCCESS;
121 };
122
123 DVector r_ini{0.0};
124 ODESolver r_ode(1e-3, 1e-10, 1e-10);
125 r_ode.solve(deriv, z_arr, r_ini);
126 r_arr = r_ode.get_data_by_component(0);
127
128 // Make a spline of r(z)
129 FML::INTERPOLATION::SPLINE::Spline r_of_z_spline(z_arr, r_arr, "r_of_z_spline");
130
131 // Fetch ndim from particles and check that we have the right dimensions
132 assert_mpi(FML::PARTICLE::GetNDIM(U()) == 3,
133 "[EquitorialToCartesianCoordinates] Particles must have ndim = 3");
134
135 // Make particles and convert to cartesian coordinates
136 particles_xyz = std::vector<U>(ngalaxies);
137 double max_x = -1e100;
138 double max_y = -1e100;
139 double max_z = -1e100;
140 double min_x = +1e100;
141 double min_y = +1e100;
142 double min_z = +1e100;
143
144 const double degrees_to_radial = 2.0 * M_PI / 360.0;
145#ifdef USE_OMP
146#pragma omp parallel for reduction(max : max_x, max_y, max_z) reduction(min : min_x, min_y, min_z)
147#endif
148 for (size_t i = 0; i < ngalaxies; i++) {
149 auto * Pos = FML::PARTICLE::GetPos(particles_xyz[i]);
150 const double RA = FML::PARTICLE::GetRA(galaxies_ra_dec_z[i]);
151 const double DEC = FML::PARTICLE::GetDEC(galaxies_ra_dec_z[i]);
152 const double redshift = FML::PARTICLE::GetRedshift(galaxies_ra_dec_z[i]);
153 const double r = r_of_z_spline(redshift);
154
155 const double cosTheta = std::cos((90.0 - DEC) * degrees_to_radial);
156 const double sinTheta = std::sin((90.0 - DEC) * degrees_to_radial);
157 const double phi = RA * degrees_to_radial;
158
159 const double x = r * sinTheta * std::cos(phi);
160 const double y = r * sinTheta * std::sin(phi);
161 const double z = r * cosTheta;
162
163 // Compute max/min
164 max_x = std::max(x, max_x);
165 max_y = std::max(y, max_y);
166 max_z = std::max(z, max_z);
167 min_x = std::min(x, min_x);
168 min_y = std::min(y, min_y);
169 min_z = std::min(z, min_z);
170
171 // Assign positions
172 Pos[0] = x;
173 Pos[1] = y;
174 Pos[2] = z;
175 }
176
177 min_max_x = {min_x, max_x};
178 min_max_y = {min_y, max_y};
179 min_max_z = {min_z, max_z};
180 }
DVector
FML::INTERPOLATION::SPLINE::DVector DVector
Definition: Simulation.h:39
ODESolver
FML::SOLVERS::ODESOLVER::ODESolver ODESolver
Definition: Simulation.h:41
ODEFunction
FML::SOLVERS::ODESOLVER::ODEFunction ODEFunction
Definition: Simulation.h:42
FML::INTERPOLATION::SPLINE::GSLSpline
This is a wrapper class for easy use of GSL splines It is (OpenMP) thread safe as long as splines are...
Definition: Spline.h:80
FML::SOLVERS::ODESOLVER::ODESolver
This is a wrapper around the GSL library to easily solve ODEs and return the data in whatever format ...
Definition: ODESolver.h:87
FML::PARTICLE::GetPos
get_vel constexpr double * GetPos(...)
Definition: ReflectOnParticleMethods.h:128
FML::PARTICLE::GetDEC
constexpr double GetDEC(...)
Definition: ReflectOnParticleMethods.h:321
FML::PARTICLE::GetRA
set_RA set_DEC set_z set_weight constexpr double GetRA(...)
Definition: ReflectOnParticleMethods.h:313
FML::PARTICLE::GetRedshift
constexpr double GetRedshift(...)
Definition: ReflectOnParticleMethods.h:329
FML::SOLVERS::ODESOLVER::ODEFunction
std::function< int(double, const double *, double *)> ODEFunction
Definition: ODESolver.h:33
FML::SOLVERS::ODESOLVER::DVector
std::vector< double > DVector
Definition: ODESolver.h:29
FML::y
return y
Definition: Global.h:265
example.z
float z
Definition: example.py:56
test.x
x
Definition: test.py:18

◆ GalaxiesRandomsToBox()

template<class T , class U >
void FML::SURVEY::GalaxiesRandomsToBox ( const T *  galaxies_ra_dec_z,
size_t  ngalaxies,
const T *  randoms_ra_dec_z,
size_t  nrandoms,
std::vector< U > &  galaxies_xyz,
std::vector< U > &  randoms_xyz,
std::function< double(double)> &  hubble_over_c_of_z,
double &  boxsize,
bool  shiftPositions,
bool  scalePositions,
std::vector< double > &  observer_position,
bool  verbose 
)

Transform galaxies and randos with positions defined by RA,DEC,Z and transform these to cartesian positions in [0,1.

Useful if we have more catalogues like galaxy + randoms we should use the same shift and same boxsize so do it together.

Template Parameters
TParticle class for the galaxies/randoms
UParticle class for the particles we make from the galaxies/randoms
Parameters
[in]galaxies_ra_dec_zGalaxies with RA, DEC and Z.
[in]ngalaxiesNumber of galaxies
[in]randoms_ra_dec_zRandom galaxies with RA, DEC and Z.
[in]nrandomsNumber of random galaxies
[out]galaxies_xyzVector with galaxies as particles with cartesian coordinates.
[out]randoms_xyzVector with random galaxies as particles with cartesian coordinates.
[in]hubble_over_c_of_zThis is the function \( H(z)/c \) used to compute the redshift-comobing distance relationship. Postions units is the same as those of \( c/H(z) \) (so e.g. if you want Mpc/h then we need H0 = 100, if you want kpc/s then use H0 = 10^5 and so on.
[out]boxsizeThe boxsize we need to place the galaxies and randoms in a cubic box.
[in]shiftPositionsShift the positions such that all are >=0
[in]scalePositionsScale positions so that all are in [0,1). This requires also shifting positions.
[out]observer_positionIf we don't shift the positions the observer is at (0,0,0). Otherwise we return the position of the observer after shifting the box.
[in]verbosePrint info while doing this.

Definition at line 310 of file GalaxiesToBox.h.

321 {
322
323 verbose = verbose and FML::ThisTask == 0;
324 double boxsize1, boxsize2;
325
326 GalaxiesToBox(
327 galaxies_ra_dec_z, ngalaxies, galaxies_xyz, hubble_over_c_of_z, boxsize1, observer_position, false, false, verbose);
328 GalaxiesToBox(randoms_ra_dec_z, nrandoms, randoms_xyz, hubble_over_c_of_z, boxsize2, observer_position, false, false, verbose);
329
330 // Set observer position to be at origin
331 observer_position = std::vector<double>(3,0.0);
332
333 double max_x = -1e100, max_y = -1e100, max_z = -1e100;
334 double min_x = +1e100, min_y = +1e100, min_z = +1e100;
335#ifdef USE_OMP
336#pragma omp parallel for reduction(max : max_x, max_y, max_z) reduction(min : min_x, min_y, min_z)
337#endif
338 for (size_t i = 0; i < galaxies_xyz.size(); i++) {
339 auto * Pos = FML::PARTICLE::GetPos(galaxies_xyz[i]);
340 max_x = std::max(Pos[0], max_x);
341 max_y = std::max(Pos[1], max_y);
342 max_z = std::max(Pos[2], max_z);
343 min_x = std::min(Pos[0], min_x);
344 min_y = std::min(Pos[1], min_y);
345 min_z = std::min(Pos[2], min_z);
346 }
347#ifdef USE_OMP
348#pragma omp parallel for reduction(max : max_x, max_y, max_z) reduction(min : min_x, min_y, min_z)
349#endif
350 for (size_t i = 0; i < randoms_xyz.size(); i++) {
351 auto * Pos = FML::PARTICLE::GetPos(randoms_xyz[i]);
352 max_x = std::max(Pos[0], max_x);
353 max_y = std::max(Pos[1], max_y);
354 max_z = std::max(Pos[2], max_z);
355 min_x = std::min(Pos[0], min_x);
356 min_y = std::min(Pos[1], min_y);
357 min_z = std::min(Pos[2], min_z);
358 }
359
360 max_x -= min_x;
361 max_y -= min_y;
362 max_z -= min_z;
363
364 // The boxsize we use is the maximum of the two
365 boxsize = 1.05 * std::max(std::max(max_x, max_y), max_z);
366
367 if (shiftPositions) {
368 for (auto & p : galaxies_xyz) {
369 auto * Pos = FML::PARTICLE::GetPos(p);
370 Pos[0] += min_x;
371 Pos[1] += min_y;
372 Pos[2] += min_z;
373 }
374 for (auto & p : randoms_xyz) {
375 auto * Pos = FML::PARTICLE::GetPos(p);
376 Pos[0] += min_x;
377 Pos[1] += min_y;
378 Pos[2] += min_z;
379 }
380 observer_position[0] += min_x;
381 observer_position[1] += min_y;
382 observer_position[2] += min_z;
383 }
384
385 if (scalePositions) {
386 for (auto & p : galaxies_xyz) {
387 auto * Pos = FML::PARTICLE::GetPos(p);
388 Pos[0] /= boxsize;
389 Pos[1] /= boxsize;
390 Pos[2] /= boxsize;
391 }
392 for (auto & p : randoms_xyz) {
393 auto * Pos = FML::PARTICLE::GetPos(p);
394 Pos[0] /= boxsize;
395 Pos[1] /= boxsize;
396 Pos[2] /= boxsize;
397 }
398 observer_position[0] /= boxsize;
399 observer_position[1] /= boxsize;
400 observer_position[2] /= boxsize;
401 }
402
403 if (verbose) {
404 std::cout << "Boxsize for boxing galaxies and randoms: " << boxsize << "\n";
405 std::cout << "Observer at ( " << observer_position[0] << " , " << observer_position[1] << " , " << observer_position[2] << ") in code units\n";
406 }
407 }
FML::SURVEY::GalaxiesToBox
void GalaxiesToBox(const T *galaxies_ra_dec_z, size_t ngalaxies, std::vector< U > &particles_xyz, std::function< double(double)> &hubble_over_c_of_z, double &boxsize, std::vector< double > &observer_position, bool shiftPositions, bool scalePositions, bool verbose)
Transform galaxies with positions defined by RA,DEC,Z and transform these to cartesian positions in [...
Definition: GalaxiesToBox.h:202
FML::ThisTask
int ThisTask
The MPI task number.
Definition: Global.cpp:33
example.verbose
bool verbose
Definition: example.py:19

◆ GalaxiesToBox()

template<class T , class U >
void FML::SURVEY::GalaxiesToBox ( const T *  galaxies_ra_dec_z,
size_t  ngalaxies,
std::vector< U > &  particles_xyz,
std::function< double(double)> &  hubble_over_c_of_z,
double &  boxsize,
std::vector< double > &  observer_position,
bool  shiftPositions,
bool  scalePositions,
bool  verbose 
)

Transform galaxies with positions defined by RA,DEC,Z and transform these to cartesian positions in [0,1).

Template Parameters
TParticle class for the galaxies
UParticle class for the particles we make from the galaxies
Parameters
[in]galaxies_ra_dec_zParticles with RA, DEC and Z.
[in]ngalaxiesNumber of galaxies
[out]particles_xyzVector with galaxies as particles with cartesian coordinates.
[in]hubble_over_c_of_zThis is the function \( H(z)/c \) used to compute the redshift-comobing distance relationship. Postions units is the same as those of \( c/H(z) \) (so e.g. if you want Mpc/h then we need H0 = 100, if you want kpc/s then use H0 = 10^5 and so on.
[out]boxsizeThe boxsize we need to place the particles in a cubic box.
[in]shiftPositionsShift the positions such that all are >=0
[in]scalePositionsScale positions so that all are in [0,1). This requires also shifting positions.
[in]verbosePrint info while doing this.

Definition at line 202 of file GalaxiesToBox.h.

210 {
211
212 observer_position = {0.,0.,0.};
213
214 // If we are to scale the positions then they must be shifted so they lie in [0,1)
215 if (scalePositions)
216 assert(shiftPositions);
217
218 // Fetch ndim from particles and check that we have the right dimensions
219 assert_mpi(FML::PARTICLE::GetNDIM(U()) == 3,
220 "[GalaxiesToBox] Particles must have ndim = 3 ");
221
222#ifdef USE_MPI
223 // If we run with MPI only print once
224 if (verbose) {
225 int mpi_rank = 0, mpi_size = 1;
226 MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
227 MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
228 verbose = (mpi_rank == 0);
229 }
230#endif
231
232 // To cartesian coordinates
233 std::pair<double, double> min_max_x, min_max_y, min_max_z;
234 EquitorialToCartesianCoordinates(
235 galaxies_ra_dec_z, ngalaxies, particles_xyz, hubble_over_c_of_z, min_max_x, min_max_y, min_max_z);
236
237 double min_x = min_max_x.first;
238 double min_y = min_max_y.first;
239 double min_z = min_max_z.first;
240 double max_x = min_max_x.second;
241 double max_y = min_max_y.second;
242 double max_z = min_max_z.second;
243
244 max_x -= min_x;
245 max_y -= min_y;
246 max_z -= min_z;
247
248 // Compute the minimum boxsize we can fit these particle in
249 boxsize = (1. + 1e-10) * std::max(std::max(max_x, max_y), max_z);
250
251 if (verbose) {
252 std::cout << "Minimum Pos: " << min_x << " " << min_y << " " << min_z << "\n";
253 std::cout << "Maximum Pos: " << max_x << " " << max_y << " " << max_z << "\n";
254 std::cout << "Boxsize: " << boxsize << "\n";
255 }
256
257 // Shift positions and possily scale them so they are in [0,1)
258#ifdef USE_OMP
259#pragma omp parallel for
260#endif
261 for (size_t i = 0; i < ngalaxies; i++) {
262 auto * Pos = FML::PARTICLE::GetPos(particles_xyz[i]);
263
264 if (shiftPositions) {
265 Pos[0] -= min_x;
266 Pos[1] -= min_y;
267 Pos[2] -= min_z;
268 observer_position[0] -= min_x;
269 observer_position[1] -= min_y;
270 observer_position[2] -= min_z;
271 }
272
273 if (scalePositions) {
274 Pos[0] /= boxsize;
275 Pos[1] /= boxsize;
276 Pos[2] /= boxsize;
277 observer_position[0] /= boxsize;
278 observer_position[1] /= boxsize;
279 observer_position[2] /= boxsize;
280 }
281 }
282 }
FML::SURVEY::EquitorialToCartesianCoordinates
void EquitorialToCartesianCoordinates(T *galaxies_ra_dec_z, size_t ngalaxies, std::vector< U > &particles_xyz, std::function< double(double)> &hubble_over_c_of_z, std::pair< double, double > &min_max_x, std::pair< double, double > &min_max_y, std::pair< double, double > &min_max_z)
Take a set of galaxies galaxies_ra_dec_z with (RA,DEC,z) and convert them to cartesian coordinates (x...
Definition: GalaxiesToBox.h:83
Generated by doxygen 1.9.3