from functools import partial
from pathlib import Path
import numpy as np
from dautil.IO import makedirs
c.f. 21cmFAST 2018 paper, Table 1
Followings are first cleaned up using gcc -fpreprocessed -dD -E -P
ANAL_PARAMS_H = '''#ifndef _ANAL_PARAM_H_
#define _ANAL_PARAM_H_
#define N_GAMMA_UV (float) (5000)
#define STELLAR_BARYON_FRAC (double) ({STELLAR_BARYON_FRAC})
#define STELLAR_BARYON_PL (double) ({STELLAR_BARYON_PL})
#define ESC_FRAC (double) ({ESC_FRAC})
#define ESC_PL (double) ({ESC_PL})
#define M_TURNOVER (double) ({M_TURNOVER})
#define INHOMO_RECO (int) (1)
#define t_STAR (float) ({t_STAR})
#define ALPHA_UVB (float) (5)
#define SHARP_CUTOFF (int) (0)
#define R_BUBBLE_MAX (float) (50)
#define EVOLVE_DENSITY_LINEARLY (int) (0)
#define SMOOTH_EVOLVED_DENSITY_FIELD (int) (1)
#define R_smooth_density (float) (0.2)
#define SECOND_ORDER_LPT_CORRECTIONS (int) (1)
#define HII_ROUND_ERR (float) (1e-3)
#define FIND_BUBBLE_ALGORITHM (int) (2)
#define R_BUBBLE_MIN (float) (L_FACTOR*1)
#define USE_HALO_FIELD (int) (0)
#define N_POISSON (int) (5)
#define T_USE_VELOCITIES (int) (1)
#define MAX_DVDR (float) (0.2)
#define VELOCITY_COMPONENT (int) (3)
#define SUBCELL_RSD (int) (1)
#define N_RSD_STEPS (int) (20)
#define DIMENSIONAL_T_POWER_SPEC (int) (1)
#define DELTA_R_FACTOR (float) (1.1)
#define DELTA_R_HII_FACTOR (float) (1.1)
#define R_OVERLAP_FACTOR (float) (1.0)
#define DELTA_CRIT_MODE (int) (1)
#define HALO_FILTER (int) (0)
#define HII_FILTER (int) (1)
#define OPTIMIZE (int) (0)
#define OPTIMIZE_MIN_MASS (float) (1e11)
#define SIZE_RANDOM_SEED (-23456789)
#define LOS_RANDOM_SEED (-123456789)
#define INITIAL_REDSHIFT (float) 300
#define HII_EFF_FACTOR (float) ( N_GAMMA_UV * STELLAR_BARYON_FRAC * ESC_FRAC)
#define L_FACTOR (float) (0.620350491)
#define HII_D (unsigned long long) (HII_DIM)
#define HII_MIDDLE (HII_DIM/2)
#define HII_MID ((unsigned long long)HII_MIDDLE)
#define HII_TOT_NUM_PIXELS (unsigned long long)(HII_D*HII_D*HII_D)
#define HII_TOT_FFT_NUM_PIXELS ((unsigned long long)(HII_D*HII_D*2llu*(HII_MID+1llu)))
#define HII_KSPACE_NUM_PIXELS ((unsigned long long)(HII_D*HII_D*(HII_MID+1llu)))
#define HII_C_INDEX(x,y,z) ((unsigned long long)((z)+(HII_MID+1llu)*((y)+HII_D*(x))))
#define HII_R_FFT_INDEX(x,y,z) ((unsigned long long)((z)+2llu*(HII_MID+1llu)*((y)+HII_D*(x))))
#define HII_R_INDEX(x,y,z) ((unsigned long long)((z)+HII_D*((y)+HII_D*(x))))
#endif
'''
HEAT_PARAMS_H = '''#ifndef _HEAT_PARAMS_H_
#define _HEAT_PARAMS_H_
#define USE_TS_IN_21CM (int) (1)
#define L_X (double) ({L_X})
#define NU_X_BAND_MAX (double) (2000*NU_over_EV)
#define NU_X_THRESH (double) (500*NU_over_EV)
#define X_RAY_SPEC_INDEX (double) (1.0)
#define HEAT_FILTER (int) 0
#define CLUMPING_FACTOR (double) (2)
#define Z_HEAT_MAX (float) (30)
#define R_XLy_MAX (float) (500)
#define NUM_FILTER_STEPS_FOR_Ts (int) (40)
#define ZPRIME_STEP_FACTOR (float) (1.02)
#define TK_at_Z_HEAT_MAX (double) (-1)
#define XION_at_Z_HEAT_MAX (double) (-1)
#define Ts_verbose (int) (1)
#define RECFAST_FILENAME (const char *) "../External_tables/recfast_LCDM.dat"
#define STELLAR_SPECTRA_FILENAME (const char *) "../External_tables/stellar_spectra.dat"
#define KAPPA_EH_FILENAME (const char *) "../External_tables/kappa_eH_table.dat"
#define KAPPA_PH_FILENAME (const char *) "../External_tables/kappa_pH_table.dat"
#define Pop (int) (2)
#define Pop2_ion (float) (N_GAMMA_UV)
#define Pop3_ion (float) (50000)
#define DEBUG_ON (int) (0)
#endif
'''
# first value is default
# rest are alternative value for each realization
PARAMS = {
'STELLAR_BARYON_FRAC': (0.05, 0.001, 1.),
'STELLAR_BARYON_PL': (0.5, -0.5, 1.),
'ESC_FRAC': (0.01, 0.001, 1.),
'ESC_PL': (-0.5, -1., 0.5),
'M_TURNOVER': (5e8, 1e8, 1e10),
't_STAR': (0.5, 0., 1.),
'L_X': (40.5, 38., 42.)
}
def get_param(params, idxs):
return {key: value[idx] for idx, (key, value) in zip(idxs, params.items())}
def gen_idxes(n=7, idxs=(1, 2)):
'''`n`: no. of levels
`idxs`: possible index values for alternative realization
'''
def gen_idxs(n, level, idx):
idxs = np.zeros(n, dtype=np.int32)
idxs[level] = idx
return idxs
# first one is all defaults
return [np.zeros(n, dtype=np.int32)] + [gen_idxs(n, level, idx) for level in range(n) for idx in idxs]
def gen_h_files(params):
# ANAL_PARAMS_H
anal_dict = {
key: params[key]
for key in ('STELLAR_BARYON_FRAC', 'STELLAR_BARYON_PL', 'ESC_FRAC', 'ESC_PL', 'M_TURNOVER', 't_STAR')
}
heat_dict = {'L_X': params['L_X']}
return ANAL_PARAMS_H.format(**anal_dict), HEAT_PARAMS_H.format(**heat_dict)
def gen_dirname(basedir, idxs, subdir='Parameter_files'):
return basedir / '_'.join(map(str, idxs)) / subdir
def main(basedir):
'''write header files for all realization parameters to `basedir`.
'''
idxes = gen_idxes()
for idxs in idxes:
dirname = gen_dirname(basedir, idxs)
params = get_param(PARAMS, idxs)
file_contents = gen_h_files(params)
makedirs(dirname)
for name, content in zip(('ANAL_PARAMS.H', 'HEAT_PARAMS.H'), file_contents):
with open(dirname / name, 'w') as f:
print(content, file=f)
basedir = Path('~/21cmfast/21cmFAST/Parameter_spaces').expanduser()
main(basedir)
idxes = gen_idxes()
for idxs in idxes:
print('_'.join(map(str, idxs)), end=' ')