forked from sbird/GenPK
-
Notifications
You must be signed in to change notification settings - Fork 0
/
gen-pk.h
128 lines (112 loc) · 6.2 KB
/
gen-pk.h
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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
/* Copyright (c) 2009, Simeon Bird <[email protected]>
*
* Permission to use, copy, modify, and/or distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */
/** \file
* GenPK Header*/
#ifndef GEN_PK_HEADER
#define GEN_PK_HEADER
#include <fftw3.h>
#include "gadgetreader.hpp"
#include <string>
//This is the type for reading particles
#ifdef DOUBLE_PRECISION_SNAP
#define GENPK_FLOAT_TYPE double
#else
#define GENPK_FLOAT_TYPE float
#endif
//This is the type for internal computations
#define GENFLOAT double
/** Returns the next power of two greater than its argument*/
int nexttwo(int);
/** Function to print the calculated power spectrum to a file
* Will print three columns, k_eff P(k) and N_modes.
* @param filename File to print to
* @param nrbins Number of bins to print
* @param keffs Effective k; this is the center of the bins, weighted by the fact that there are more
* modes at higher k
* @param power Power spectrum
* @param count Number of modes in each bin*/
int print_pk(std::string filename, int nrbins, double * keffs, double * power, int* count);
/** Prints usage info */
void help(void);
/** Function to convert particle type to a string
* @param type Particle type, such as BARYON_TYPE*/
std::string type_str(int type);
/** Wrapper function to read the data and pass it to fieldize in chunks.
* Returns 1 is some error occured allocating memory, and zero if successful.
* @param field Pointer to memory for output field
* @param snap GadgetReader object for snapshot
* @param type Particle type to read; since POS always has all types, complications are minimal
* @param box Boxsize, in kpc
* @param field_dims Length of side of the field above. */
int read_fieldize(GENFLOAT * field, GadgetReader::GSnap* snap, int type, double box, int field_dims, double * total_mass);
/** Wrapper function to read the data from an HDF5 snapshot and pass it to fieldize.
* Returns 1 is some error occured allocating memory, and zero if successful.
* @param field Pointer to memory for output field
* @param ffname filename of hdf snapshot
* @param type Particle type to read
* @param box Boxsize, in kpc
* @param field_dims Length of side of the field above.
* @param fileno number of file to use.
* @param total_mass Variable to store total mass of particles */
int read_fieldize_hdf5(GENFLOAT * field, const char *ffname, int type, double box, int field_dims, double * total_mass,int fileno);
/* this routine loads header data from the first file of an HDF5 snapshot.*/
int load_hdf5_header(const char *ffname, double *atime, double *redshift, double *box100, double *h100, int64_t *npart_all, double * mass);
/** Finds a snapshot set of hdf files from the given initial filename
*/
std::vector<std::string> find_hdf_set(const std::string& infname);
/** Takes an array of particle positions and places them into a grid, ready for use by FFTW.
* OpenMP-parallel. Returns 0. Uses a Cloud-In-Cell algorithm.
* @param boxsize The size of the box, in kpc.
* @param dims The size of the grid to use
* @param out Pointer to output array, of size dims**3. Should be already initialised.
* @param total_particles Total number of particles that will be placed on the grid. Used for density calculation.
* @param segment_particles Number of particles to place on the grid in this call.
* @param positions Array of particle positions, of size 3*segment_particles (like the output of GetBlocks)
* @param masses Array containing particle masses, if variable. Null if particle masses are constant.
* @param mass Particle mass if constant for all of this type.
* @param extra If this is 1, assume that the output is about to be handed to an FFTW in-place routine,
* and make it skip the last 2 places of the each row in the last dimension */
int fieldize(double boxsize, int dims, GENFLOAT *out, int64_t segment_particles, GENPK_FLOAT_TYPE *positions, GENPK_FLOAT_TYPE * masses, double mass, int extra);
#ifdef __cplusplus
extern "C" {
#endif
/** The inverse window function associated with the Cloud-in-cell algorithm.
* The power spectrum needs to deconvolve this.
* @param kx k_x
* @param ky k_y
* @param kz k_z
* @param n Total size of grid */
GENFLOAT invwindow(int kx, int ky, int kz, int n);
/** Wrapper around FFTW to calculate the 3D power spectrum from a 3D field.
* Returns 0.
* @param dims Size of grid to FFT.
* @param outfield Pointer to memory where the FFT stores its output, as was specified to the FFTW plan.
* This is likely to be an in-place transform, in which case
* we will need some contiguous memory space after the actual data in field.
* The real input data has size dims**3
* The output has size dims*dims*(dims/2+1) *complex* values
* So we must allocate 2*dims*dims*(dims/2+1)*sizeof(GENFLOAT).
* @param outfield2 Pointer to second array, just like outfield. What is computed is outfield*outfield2. May be the same array.
* @param nrbins Number of bins in the output.
* @param power Pointer to memory to output powerspectrum to. Needs to store nrbins values.
* @param count Ditto for modes per bin
* @param keffs Ditto for effective k.*/
int powerspectrum(int dims, fftw_complex* outfield, fftw_complex* outfield2, int nrbins, double *power, int *count, double *keffs, double total_mass, double total_mass2);
int is_bigfile(const char * infile);
int load_bigfile_header(const char *fname, double *atime, double *redshift, double *box100, double *h100, int64_t *npart_all, double * mass, double *Omega0);
int read_fieldize_bigfile(GENFLOAT * field, const char *fname, int type, double box, int field_dims, double *total_mass, int64_t* npart_all, double * mass, double Omega0);
#ifdef __cplusplus
}
#endif
#endif