From 90a4bd2e8b020c074cd8fcf5ff7d1d58560958b7 Mon Sep 17 00:00:00 2001 From: Benjamin Moster Date: Wed, 8 Apr 2020 21:11:03 +0200 Subject: [PATCH] Output routines reorganized and bug for halo radius fixed (now in physical units rather than comoving) --- README.md | 2 +- src/allvars.h | 2 +- src/galaxies.c | 2 +- src/output.c | 750 ++++++++++++++++++++++++----------------------- src/read_trees.c | 2 +- src/setup.c | 25 +- user-guide.pdf | Bin 464563 -> 464558 bytes 7 files changed, 414 insertions(+), 369 deletions(-) diff --git a/README.md b/README.md index ab06cf5..a4b7eed 100644 --- a/README.md +++ b/README.md @@ -1040,7 +1040,7 @@ the redshift for all correlation functions must be the same (given at the start ## 6) Halo merger tree files In addition to observed galaxy data, Emerge relies on simulated dark matter halo merger trees that will be populated with galaxies. The -files that containt these merger trees and their locations can be specified in the parameter file with the parameter `TreefileName`. The +files that contain these merger trees and their locations can be specified in the parameter file with the parameter `TreefileName`. The suggested location for the merger tree files is the `trees` folder. If the trees are located in a single file, `TreefileName` corresponds to the name of this file including folder, e.g. `trees/P100`. If the trees are stored within multiple files, `TreefileName` gives the file base name, which will be appended with the file index, e.g. `trees/P100` corresponds to the files `P100.0`, `P100.1`, `P100.2`, `...`, diff --git a/src/allvars.h b/src/allvars.h index bede944..0c1bf3b 100644 --- a/src/allvars.h +++ b/src/allvars.h @@ -36,7 +36,7 @@ #define CODE "EMERGE" ///< Code name #define BRANCH "master" ///< Branch name -#define VERSION "1.0.0" ///< Code version +#define VERSION "1.0.1" ///< Code version #define SYMBOL '#' ///< Symbol used to start output lines #define NSTRING 200 ///< Number of elements used for a generic string diff --git a/src/galaxies.c b/src/galaxies.c index 7f2e83b..aeaf0fc 100644 --- a/src/galaxies.c +++ b/src/galaxies.c @@ -698,7 +698,7 @@ float get_tdf(int ihalo, int imain) dx = NEAREST(H[ihalo].pos[0]-H[imain].pos[0]); dy = NEAREST(H[ihalo].pos[1]-H[imain].pos[1]); dz = NEAREST(H[ihalo].pos[2]-H[imain].pos[2]); - rsat = sqrt(dx*dx+dy*dy+dz*dz); + rsat = sqrt(dx*dx+dy*dy+dz*dz)*H[imain].a; clog = log(1.0+(H[imain].mvir+H[imain].mstar)/(H[ihalo].mvir+H[ihalo].mstar)); mu = (H[ihalo].mvir+H[ihalo].mstar)/(H[imain].mvir+H[imain].mstar); tau = (H[imain].rvir*All.x_unit/vvir) * MPC_IN_KM / YR_IN_SEC / All.t_unit; diff --git a/src/output.c b/src/output.c index b310ce1..bb03993 100644 --- a/src/output.c +++ b/src/output.c @@ -35,7 +35,7 @@ */ void output_galaxies(void) { - int i, j, k, ihalo, iprog, filenr, masterTask, lastTask, ngroups, maxlen, task; + int i, j, ihalo, iprog, filenr, masterTask, lastTask, ngroups, maxlen, task; int left_to_send, p, send_this_turn, offset, nhaloout; float sigmaobs, mstarobs, sfrobs, scatter, amhalf; char buf[NSTRING]; @@ -44,12 +44,13 @@ void output_galaxies(void) //Structure that stores all halo/galaxy properties that will be saved struct haloout { - float mh,mdot,mhp,mhh,amp,amh,r,c,l,ms,sfr,icm,mso,sfro,x,y,z,u,v,w; + float mh,mdot,mhp,mpd,mhh,amp,amh,r,c,l,ms,sfr,icm,mso,sfro,x,y,z,u,v,w; unsigned short t; long long hid, did, uid; } *hp, hsend; #ifdef HDF5_SUPPORT + int k; hid_t gal_file, hdf5_dataspace_memory, hdf5_dataspace_in_file, hdf5_dtype, hdf5_ftype, hdf5_filespace; hid_t hdf5_properties, hdf5_dataset; hid_t hdf5_att_space, hdf5_att_redshift, hdf5_att_lbox, hdf5_att_hubble; @@ -78,17 +79,17 @@ void output_galaxies(void) distribute_file(All.NTaskPerUniverse, All.NumOutputFiles, 0, 0, All.NTaskPerUniverse - 1, &filenr, &masterTask, &lastTask); //Go through all output redshifts - for(i = 0; i < All.Noutputredshifts; i++) + for (i = 0; i < All.Noutputredshifts; i++) { //If the output format has been set to HDF5 if (All.OutputFormat == 2) { - //Check if HDF5 libraries are set + //If HDF5 libraries are set #ifdef HDF5_SUPPORT if (All.NumOutputFiles > 1) - sprintf(buf, "%s/galaxies.S%d.%d.h5", All.OutputDir, Output_iscale[i], filenr); + sprintf(buf, "%s/galaxies/galaxies.S%d.%d.h5", All.OutputDir, Output_iscale[i], filenr); else - sprintf(buf, "%s/galaxies.S%d.h5", All.OutputDir, Output_iscale[i]); + sprintf(buf, "%s/galaxies/galaxies.S%d.h5", All.OutputDir, Output_iscale[i]); //If not return #else if (ThisTask == 0) printf("%s Output format has been set to 2 but HDF5 support was not enabled.\n",All.startline); @@ -99,22 +100,22 @@ void output_galaxies(void) else { if (All.NumOutputFiles > 1) - sprintf(buf, "%s/galaxies.S%d.%d", All.OutputDir, Output_iscale[i], filenr); + sprintf(buf, "%s/galaxies/galaxies.S%d.%d", All.OutputDir, Output_iscale[i], filenr); else - sprintf(buf, "%s/galaxies.S%d.out", All.OutputDir, Output_iscale[i]); + sprintf(buf, "%s/galaxies/galaxies.S%d.out", All.OutputDir, Output_iscale[i]); } //Get number of groups ngroups = All.NumOutputFiles / All.NumFilesInParallel; - if((All.NumOutputFiles % All.NumFilesInParallel)) ngroups++; + if ((All.NumOutputFiles % All.NumFilesInParallel)) ngroups++; //For each group do... - for(j = 0; j < ngroups; j++) + for (j = 0; j < ngroups; j++) { //This task will be processed now - if((filenr / All.NumFilesInParallel) == j && MasterTask == 0) + if ((filenr / All.NumFilesInParallel) == j && MasterTask == 0) { //The masterTask opens the file - if(ThisTask == masterTask) + if (ThisTask == masterTask) { //If the hdf5 output format has been selected if (All.OutputFormat == 2) { //Check if the libraries have been included @@ -125,36 +126,40 @@ void output_galaxies(void) printf("%s Writing output file at z = %f: '%s' (file %d of %d)\n",All.startline, OutputRedshifts[i], buf, filenr+1, All.NumOutputFiles); //Specify data type hdf5_dtype = H5Tcreate(H5T_COMPOUND, sizeof(struct haloout)); - hdf5_status = H5Tinsert(hdf5_dtype, "Halo_mass", HOFFSET(struct haloout, mh), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Halo_growth_rate", HOFFSET(struct haloout, mdot), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Halo_mass_peak", HOFFSET(struct haloout, mhp), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Halo_mass_host", HOFFSET(struct haloout, mhh), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Scale_peak_mass", HOFFSET(struct haloout, amp), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Scale_half_mass", HOFFSET(struct haloout, amh), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Halo_radius", HOFFSET(struct haloout, r), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Concentration", HOFFSET(struct haloout, c), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Halo_spin", HOFFSET(struct haloout, l), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Stellar_mass", HOFFSET(struct haloout, ms), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "SFR", HOFFSET(struct haloout, sfr), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Intra_cluster_mass", HOFFSET(struct haloout, icm), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Stellar_mass_obs", HOFFSET(struct haloout, mso), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "SFR_obs", HOFFSET(struct haloout, sfro), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "X_pos", HOFFSET(struct haloout, x), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Y_pos", HOFFSET(struct haloout, y), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Z_pos", HOFFSET(struct haloout, z), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "X_vel", HOFFSET(struct haloout, u), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Y_vel", HOFFSET(struct haloout, v), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Z_vel", HOFFSET(struct haloout, w), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Type", HOFFSET(struct haloout, t), H5T_NATIVE_USHORT); - hdf5_status = H5Tinsert(hdf5_dtype, "Halo_ID", HOFFSET(struct haloout, hid), H5T_NATIVE_LLONG); - hdf5_status = H5Tinsert(hdf5_dtype, "Desc_ID", HOFFSET(struct haloout, did), H5T_NATIVE_LLONG); - hdf5_status = H5Tinsert(hdf5_dtype, "Up_ID", HOFFSET(struct haloout, uid), H5T_NATIVE_LLONG); + k = 0; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_mass", HOFFSET(struct haloout, mh), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_growth_rate", HOFFSET(struct haloout, mdot), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_mass_peak", HOFFSET(struct haloout, mhp), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_growth_peak", HOFFSET(struct haloout, mpd), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_mass_host", HOFFSET(struct haloout, mhh), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Scale_peak_mass", HOFFSET(struct haloout, amp), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Scale_half_mass", HOFFSET(struct haloout, amh), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_radius", HOFFSET(struct haloout, r), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Concentration", HOFFSET(struct haloout, c), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_spin", HOFFSET(struct haloout, l), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Stellar_mass", HOFFSET(struct haloout, ms), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "SFR", HOFFSET(struct haloout, sfr), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Intra_cluster_mass", HOFFSET(struct haloout, icm), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Stellar_mass_obs", HOFFSET(struct haloout, mso), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "SFR_obs", HOFFSET(struct haloout, sfro), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "X_pos", HOFFSET(struct haloout, x), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Y_pos", HOFFSET(struct haloout, y), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Z_pos", HOFFSET(struct haloout, z), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "X_vel", HOFFSET(struct haloout, u), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Y_vel", HOFFSET(struct haloout, v), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Z_vel", HOFFSET(struct haloout, w), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Type", HOFFSET(struct haloout, t), H5T_NATIVE_USHORT); k+=2; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_ID", HOFFSET(struct haloout, hid), H5T_NATIVE_LLONG); k+=8; + hdf5_status = H5Tinsert(hdf5_dtype, "Desc_ID", HOFFSET(struct haloout, did), H5T_NATIVE_LLONG); k+=8; + hdf5_status = H5Tinsert(hdf5_dtype, "Up_ID", HOFFSET(struct haloout, uid), H5T_NATIVE_LLONG); k+=8; + //Specify file type + hdf5_ftype = H5Tcreate(H5T_COMPOUND, k); k = 0; - hdf5_ftype = H5Tcreate(H5T_COMPOUND, 4 * 20 + 2 + 8 * 3); hdf5_status = H5Tinsert(hdf5_ftype, "Halo_mass", k, H5T_IEEE_F32LE); k+=4; hdf5_status = H5Tinsert(hdf5_ftype, "Halo_growth_rate", k, H5T_IEEE_F32LE); k+=4; hdf5_status = H5Tinsert(hdf5_ftype, "Halo_mass_peak", k, H5T_IEEE_F32LE); k+=4; + hdf5_status = H5Tinsert(hdf5_ftype, "Halo_growth_peak", k, H5T_IEEE_F32LE); k+=4; hdf5_status = H5Tinsert(hdf5_ftype, "Halo_mass_host", k, H5T_IEEE_F32LE); k+=4; hdf5_status = H5Tinsert(hdf5_ftype, "Scale_peak_mass", k, H5T_IEEE_F32LE); k+=4; hdf5_status = H5Tinsert(hdf5_ftype, "Scale_half_mass", k, H5T_IEEE_F32LE); k+=4; @@ -181,7 +186,7 @@ void output_galaxies(void) //Otherwise open standard ascii file else { //Open file - if(!(fd = fopen(buf, "w"))) + if (!(fd = fopen(buf, "w"))) { //If not possible print to screen and abort printf("%s Can't open file `%s' for writing output.\n",All.startline,buf); endrun("file open error"); @@ -193,6 +198,7 @@ void output_galaxies(void) fprintf(fd,"#Halo_mass(%d)",task); task++; fprintf(fd," Halo_growth_rate(%d)",task); task++; fprintf(fd," Halo_mass_peak(%d)",task); task++; + fprintf(fd," Halo_growth_peak(%d)",task); task++; fprintf(fd," Halo_mass_host(%d)",task); task++; fprintf(fd," Scale_peak_mass(%d)",task); task++; fprintf(fd," Scale_half_mass(%d)",task); task++; @@ -228,6 +234,7 @@ void output_galaxies(void) fprintf(fd,"#Halo_mass: Current virial mass of the halo (log Msun).\n"); fprintf(fd,"#Halo_growth_rate: Growth rate of the halo (Msun/yr)\n"); fprintf(fd,"#Halo_mass_peak: Peak virial mass of the halo through its history (log Msun).\n"); + fprintf(fd,"#Halo_growth_peak: Growth rate of the halo when peak mass is reached (Msun/yr).\n"); fprintf(fd,"#Halo_mass_host: Current virial mass of the host halo (log Msun).\n"); fprintf(fd,"#Scale_peak_mass: Scale Factor when halo had its peak mass.\n"); fprintf(fd,"#Scale_half_mass: Scale Factor when halo first had half of its peak mass.\n"); @@ -257,252 +264,251 @@ void output_galaxies(void) for (ihalo = 0; ihalo < Nhalos; ihalo++) if (H[ihalo].iscale == Output_iscale[i] && H[ihalo].gone == 0 && H[ihalo].mstar >= All.minmass) nhaloout++; - if (nhaloout > 0) +#ifdef HDF5_SUPPORT + if (ThisTask == masterTask && All.OutputFormat == 2) { + //Set the path name for the galaxy data set + sprintf(path,"/Galaxies"); + //Initial size is the number of haloes on the master task + dims[0] = nhaloout; + //Allow for variable total size + maxdims[0] = H5S_UNLIMITED; + //Create the data space + hdf5_dataspace_in_file = H5Screate_simple(1, dims, maxdims); + //Create the data properties + hdf5_properties = H5Pcreate(H5P_DATASET_CREATE); + hdf5_status = H5Pset_chunk(hdf5_properties, 1, dims); // set chunk size + hdf5_status = H5Pset_shuffle(hdf5_properties); // reshuffle bytes to get better compression ratio + hdf5_status = H5Pset_deflate(hdf5_properties, 9); // gzip compression level 9 + hdf5_status = H5Pset_fletcher32(hdf5_properties); // Fletcher32 checksum on dataset + //If filters could be set use them to create the data set + if (H5Pall_filters_avail(hdf5_properties)) + hdf5_dataset = H5Dcreate(gal_file, path, hdf5_ftype, hdf5_dataspace_in_file, H5P_DEFAULT, hdf5_properties, H5P_DEFAULT); + //Otherwise create the default data set + else + hdf5_dataset = H5Dcreate(gal_file, path, hdf5_ftype, hdf5_dataspace_in_file, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + //Initialise the number of haloes in the set to zero + sendsum = 0; + //Define Scalar for Attributes + hdf5_att_space = H5Screate(H5S_SCALAR); + //Write Redshift as attribute to dataset + hdf5_att_redshift = H5Acreate(hdf5_dataset, "Scale Factor", H5T_NATIVE_FLOAT, hdf5_att_space, H5P_DEFAULT, H5P_DEFAULT); + H5Awrite(hdf5_att_redshift, H5T_IEEE_F32LE, &ScaleFactor[Output_iscale[i]]); + H5Aclose(hdf5_att_redshift); + //Write Box Size as attribute to dataset + tmp = All.Lbox*All.x_unit; + hdf5_att_lbox = H5Acreate(hdf5_dataset, "Box Size", H5T_NATIVE_DOUBLE, hdf5_att_space, H5P_DEFAULT, H5P_DEFAULT); + H5Awrite(hdf5_att_lbox, H5T_IEEE_F64LE, &tmp); + H5Aclose(hdf5_att_lbox); + //Write Hubble parameter as attribute to dataset + hdf5_att_hubble = H5Acreate(hdf5_dataset, "Hubble Parameter", H5T_NATIVE_DOUBLE, hdf5_att_space, H5P_DEFAULT, H5P_DEFAULT); + H5Awrite(hdf5_att_hubble, H5T_IEEE_F64LE, &All.h_100); + H5Aclose(hdf5_att_hubble); + //Write Omega_0 as attribute to dataset + hdf5_att_omega0 = H5Acreate(hdf5_dataset, "Omega_0", H5T_NATIVE_DOUBLE, hdf5_att_space, H5P_DEFAULT, H5P_DEFAULT); + H5Awrite(hdf5_att_omega0, H5T_IEEE_F64LE, &All.Omega_0); + H5Aclose(hdf5_att_omega0); + //Write Omega_Lambda as attribute to dataset + hdf5_att_omegal = H5Acreate(hdf5_dataset, "Omega_Lambda", H5T_NATIVE_DOUBLE, hdf5_att_space, H5P_DEFAULT, H5P_DEFAULT); + H5Awrite(hdf5_att_omegal, H5T_IEEE_F64LE, &All.Omega_Lambda_0); + H5Aclose(hdf5_att_omegal); + //Write Omega_Baryon as attribute to dataset + hdf5_att_omegab = H5Acreate(hdf5_dataset, "Omega_Baryon", H5T_NATIVE_DOUBLE, hdf5_att_space, H5P_DEFAULT, H5P_DEFAULT); + H5Awrite(hdf5_att_omegab, H5T_IEEE_F64LE, &All.Omega_Baryon_0); + H5Aclose(hdf5_att_omegab); + //Write minimum stellar mass as attribute to dataset + tmp = log10(All.minmass*All.m_unit); + hdf5_att_mmin = H5Acreate(hdf5_dataset, "Minimum Stellar Mass", H5T_NATIVE_DOUBLE, hdf5_att_space, H5P_DEFAULT, H5P_DEFAULT); + H5Awrite(hdf5_att_mmin, H5T_IEEE_F64LE, &tmp); + H5Aclose(hdf5_att_mmin); + } +#endif -#ifdef HDF5_SUPPORT - if (ThisTask == masterTask && All.OutputFormat == 2) + //Now go through all tasks that are member of this group + for (task = masterTask, offset = 0; task <= lastTask; task++) + { + //if this task is processed + if (task == ThisTask) { - //Set the path name for the galaxy data set - sprintf(path,"/Galaxies"); - //Initial size is the number of haloes on the master task - dims[0] = nhaloout; - //Allow for variable total size - maxdims[0] = H5S_UNLIMITED; - //Create the data space - hdf5_dataspace_in_file = H5Screate_simple(1, dims, maxdims); - //Create the data properties - hdf5_properties = H5Pcreate(H5P_DATASET_CREATE); - hdf5_status = H5Pset_chunk(hdf5_properties, 1, dims); // set chunk size - hdf5_status = H5Pset_shuffle(hdf5_properties); // reshuffle bytes to get better compression ratio - hdf5_status = H5Pset_deflate(hdf5_properties, 9); // gzip compression level 9 - hdf5_status = H5Pset_fletcher32(hdf5_properties); // Fletcher32 checksum on dataset - //If filters could be set use them to create the data set - if(H5Pall_filters_avail(hdf5_properties)) - hdf5_dataset = H5Dcreate(gal_file, path, hdf5_ftype, hdf5_dataspace_in_file, H5P_DEFAULT, hdf5_properties, H5P_DEFAULT); - //Otherwise create the default data set - else - hdf5_dataset = H5Dcreate(gal_file, path, hdf5_ftype, hdf5_dataspace_in_file, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - //Initialise the number of haloes in the set to zero - sendsum = 0; - //Define Scalar for Attributes - hdf5_att_space = H5Screate(H5S_SCALAR); - //Write Redshift as attribute to dataset - hdf5_att_redshift = H5Acreate(hdf5_dataset, "Scale Factor", H5T_NATIVE_FLOAT, hdf5_att_space, H5P_DEFAULT, H5P_DEFAULT); - H5Awrite(hdf5_att_redshift, H5T_IEEE_F32LE, &ScaleFactor[Output_iscale[i]]); - H5Aclose(hdf5_att_redshift); - //Write Box Size as attribute to dataset - tmp = All.Lbox*All.x_unit; - hdf5_att_lbox = H5Acreate(hdf5_dataset, "Box Size", H5T_NATIVE_DOUBLE, hdf5_att_space, H5P_DEFAULT, H5P_DEFAULT); - H5Awrite(hdf5_att_lbox, H5T_IEEE_F64LE, &tmp); - H5Aclose(hdf5_att_lbox); - //Write Hubble parameter as attribute to dataset - hdf5_att_hubble = H5Acreate(hdf5_dataset, "Hubble Parameter", H5T_NATIVE_DOUBLE, hdf5_att_space, H5P_DEFAULT, H5P_DEFAULT); - H5Awrite(hdf5_att_hubble, H5T_IEEE_F64LE, &All.h_100); - H5Aclose(hdf5_att_hubble); - //Write Omega_0 as attribute to dataset - hdf5_att_omega0 = H5Acreate(hdf5_dataset, "Omega_0", H5T_NATIVE_DOUBLE, hdf5_att_space, H5P_DEFAULT, H5P_DEFAULT); - H5Awrite(hdf5_att_omega0, H5T_IEEE_F64LE, &All.Omega_0); - H5Aclose(hdf5_att_omega0); - //Write Omega_Lambda as attribute to dataset - hdf5_att_omegal = H5Acreate(hdf5_dataset, "Omega_Lambda", H5T_NATIVE_DOUBLE, hdf5_att_space, H5P_DEFAULT, H5P_DEFAULT); - H5Awrite(hdf5_att_omegal, H5T_IEEE_F64LE, &All.Omega_Lambda_0); - H5Aclose(hdf5_att_omegal); - //Write Omega_Baryon as attribute to dataset - hdf5_att_omegab = H5Acreate(hdf5_dataset, "Omega_Baryon", H5T_NATIVE_DOUBLE, hdf5_att_space, H5P_DEFAULT, H5P_DEFAULT); - H5Awrite(hdf5_att_omegab, H5T_IEEE_F64LE, &All.Omega_Baryon_0); - H5Aclose(hdf5_att_omegab); - //Write minimum stellar mass as attribute to dataset - tmp = log10(All.minmass*All.m_unit); - hdf5_att_mmin = H5Acreate(hdf5_dataset, "Minimum Stellar Mass", H5T_NATIVE_DOUBLE, hdf5_att_space, H5P_DEFAULT, H5P_DEFAULT); - H5Awrite(hdf5_att_mmin, H5T_IEEE_F64LE, &tmp); - H5Aclose(hdf5_att_mmin); + //We need to get those haloes + left_to_send = nhaloout; + + //Tell this to the other tasks + for (p = masterTask; p <= lastTask; p++) + if (p != ThisTask) + MPI_Send(&left_to_send, 1, MPI_INT, p, TAG_NHALOS, MPI_COMM_WORLD); } -#endif + //Reveive the number of haloes from active task + else MPI_Recv(&left_to_send, 1, MPI_INT, task, TAG_NHALOS, MPI_COMM_WORLD, &status); - //Now go through all tasks that are member of this group - for(task = masterTask, offset = 0; task <= lastTask; task++) + //Now send the haloes in increments + while (left_to_send > 0) { - //if this task is processed - if(task == ThisTask) - { - //We need to get those haloes - left_to_send = nhaloout; - - //Tell this to the other tasks - for(p = masterTask; p <= lastTask; p++) - if(p != ThisTask) - MPI_Send(&left_to_send, 1, MPI_INT, p, TAG_NHALOS, MPI_COMM_WORLD); - } - //Reveive the number of haloes from active task - else MPI_Recv(&left_to_send, 1, MPI_INT, task, TAG_NHALOS, MPI_COMM_WORLD, &status); + //Compute how many haloes we can send this turn + send_this_turn = left_to_send; + //If they do not all fit into the CommBuffer set to the maximum + if (send_this_turn > maxlen) send_this_turn = maxlen; - //Now send the haloes in increments - while(left_to_send > 0) + //Sending task + if (ThisTask == task) { - //Compute how many haloes we can send this turn - send_this_turn = left_to_send; - //If they do not all fit into the CommBuffer set to the maximum - if(send_this_turn > maxlen) send_this_turn = maxlen; - - //Sending task - if(ThisTask == task) - { - //Use the haloout struct for the CommBuffer - hp = (struct haloout *) CommBuffer; - //Initialise the number of systems sent to zero - ihalo = 0; - //Go through all systems until the maximum that can be sent - while (ihalo < send_this_turn) - { //If we are at the right scale factor and the halo is not gone send it to the master - if (H[offset].iscale == Output_iscale[i] && H[offset].gone == 0 && H[offset].mstar >= All.minmass) + //Use the haloout struct for the CommBuffer + hp = (struct haloout *) CommBuffer; + //Initialise the number of systems sent to zero + ihalo = 0; + //Go through all systems until the maximum that can be sent + while (ihalo < send_this_turn) + { //If we are at the right scale factor and the halo is not gone send it to the master + if (H[offset].iscale == Output_iscale[i] && H[offset].gone == 0 && H[offset].mstar >= All.minmass) + { + //Set the scatter in stellar mass + scatter = get_gaussian_random_number(offset); + //Compute the observational error and the observed stellar mass + sigmaobs = All.obssigma0+(1./H[offset].a-1.)*All.obssigmaz; + if (1./H[offset].a-1. > ZMAX_SMFERR) sigmaobs = All.obssigma0+(ZMAX_SMFERR)*All.obssigmaz; + mstarobs = H[offset].mstar * pow(10.,sigmaobs*scatter); + //Set the scatter in the star formation rate + scatter = get_gaussian_random_number(offset + RANDOM_NUMBER_TABLE/2); + //Compute the observed star formation rate + sfrobs = H[offset].sfr * pow(10.,sigmaobs*scatter); + + //Find scale factor when virial mass was half of maximum value + amhalf = ScaleFactor[H[H[offset].impeak].iscale]; + iprog = H[H[offset].impeak].iprog; + if (iprog >=0 ) { - //Set the scatter in stellar mass - scatter = get_gaussian_random_number(offset); - //Compute the observational error and the observed stellar mass - sigmaobs = All.obssigma0+(1./H[offset].a-1.)*All.obssigmaz; - if (1./H[offset].a-1. > ZMAX_SMFERR) sigmaobs = All.obssigma0+(ZMAX_SMFERR)*All.obssigmaz; - mstarobs = H[offset].mstar * pow(10.,sigmaobs*scatter); - //Set the scatter in the star formation rate - scatter = get_gaussian_random_number(offset + RANDOM_NUMBER_TABLE/2); - //Compute the observed star formation rate - sfrobs = H[offset].sfr * pow(10.,sigmaobs*scatter); - - //Find scale factor when virial mass was half of maximum value - amhalf = ScaleFactor[H[H[offset].impeak].iscale]; - iprog = H[H[offset].impeak].iprog; - if (iprog >=0 ) + while (H[iprog].mvir > 0.5 * H[H[offset].impeak].mvir) { - while (H[iprog].mvir > 0.5 * H[H[offset].impeak].mvir) - { - if (H[iprog].iprog < 0) break; - iprog = H[iprog].iprog; - } - amhalf = ScaleFactor[H[iprog].iscale]; + if (H[iprog].iprog < 0) break; + iprog = H[iprog].iprog; } + amhalf = ScaleFactor[H[iprog].iscale]; + } - //Set all halo/galaxy properties - hsend.mh = log10(H[offset].mvir*All.m_unit); - hsend.mdot = H[offset].mdotbary*All.m_unit/All.t_unit/All.f_baryon; - hsend.mhp = log10(H[H[offset].impeak].mvir*All.m_unit); + //Set all halo/galaxy properties + hsend.mh = log10(H[offset].mvir*All.m_unit); + hsend.mdot = H[offset].mdotbary*All.m_unit/All.t_unit/All.f_baryon; + hsend.mhp = log10(H[H[offset].impeak].mvir*All.m_unit); + hsend.mpd = H[H[offset].impeak].mdotbary*All.m_unit/All.t_unit/All.f_baryon; #ifdef COMPUTE_ICM - hsend.mhh = log10(H[H[offset].ihost].mvir*All.m_unit); + hsend.mhh = log10(H[H[offset].ihost].mvir*All.m_unit); #else - hsend.mhh = 0.0; + hsend.mhh = 0.0; #endif - hsend.amp = ScaleFactor[H[H[offset].impeak].iscale]; - hsend.amh = amhalf; - hsend.r = H[offset].rvir*All.x_unit*1.e3; - hsend.c = H[offset].c; - hsend.l = H[offset].lambda; - hsend.ms = log10(H[offset].mstar*All.m_unit); - hsend.sfr = H[offset].sfr*All.m_unit/All.t_unit; - hsend.icm = log10(H[offset].icm*All.m_unit); - hsend.mso = log10(mstarobs*All.m_unit); - hsend.sfro = sfrobs*All.m_unit/All.t_unit; - hsend.x = H[offset].pos[0]; - hsend.y = H[offset].pos[1]; - hsend.z = H[offset].pos[2]; - hsend.u = H[offset].vel[0]; - hsend.v = H[offset].vel[1]; - hsend.w = H[offset].vel[2]; - hsend.t = H[offset].type; - hsend.hid = (long long)(H[offset].haloid)-1; - hsend.did = (long long)(H[offset].descid)-1; - hsend.uid = (long long)(H[offset].upid)-1; - - //Write this system to the CommBuffer - *hp++ = hsend; - //Increment the number of systems that are being sent - ihalo++; - } - //Increment the halo index - offset++; + hsend.amp = ScaleFactor[H[H[offset].impeak].iscale]; + hsend.amh = amhalf; + hsend.r = H[offset].rvir*All.x_unit*1.e3; + hsend.c = H[offset].c; + hsend.l = H[offset].lambda; + hsend.ms = log10(H[offset].mstar*All.m_unit); + hsend.sfr = H[offset].sfr*All.m_unit/All.t_unit; + hsend.icm = log10(H[offset].icm*All.m_unit); + hsend.mso = log10(mstarobs*All.m_unit); + hsend.sfro = sfrobs*All.m_unit/All.t_unit; + hsend.x = H[offset].pos[0]; + hsend.y = H[offset].pos[1]; + hsend.z = H[offset].pos[2]; + hsend.u = H[offset].vel[0]; + hsend.v = H[offset].vel[1]; + hsend.w = H[offset].vel[2]; + hsend.t = H[offset].type; + hsend.hid = (long long)(H[offset].haloid)-1; + hsend.did = (long long)(H[offset].descid)-1; + hsend.uid = (long long)(H[offset].upid)-1; + + //Write this system to the CommBuffer + *hp++ = hsend; + //Increment the number of systems that are being sent + ihalo++; } + //Increment the halo index + offset++; } + } - //Receive haloes - if(ThisTask == masterTask && task != masterTask) MPI_Recv(CommBuffer, (sizeof(struct halo)) * send_this_turn, MPI_BYTE, task, TAG_HDATA, MPI_COMM_WORLD, &status); + //Receive haloes + if (ThisTask == masterTask && task != masterTask) MPI_Recv(CommBuffer, (sizeof(struct halo)) * send_this_turn, MPI_BYTE, task, TAG_HDATA, MPI_COMM_WORLD, &status); - //Send haloes - if(ThisTask != masterTask && task == ThisTask) MPI_Ssend(CommBuffer, (sizeof(struct halo)) * send_this_turn, MPI_BYTE, masterTask, TAG_HDATA, MPI_COMM_WORLD); + //Send haloes + if (ThisTask != masterTask && task == ThisTask) MPI_Ssend(CommBuffer, (sizeof(struct halo)) * send_this_turn, MPI_BYTE, masterTask, TAG_HDATA, MPI_COMM_WORLD); - //Collect haloes and write to file - if(ThisTask == masterTask) - { + //Collect haloes and write to file + if (ThisTask == masterTask) + { - //If the output format is hdf5 - if (All.OutputFormat == 2) - { + //If the output format is hdf5 + if (All.OutputFormat == 2) + { #ifdef HDF5_SUPPORT - //Set the starting index and the count - start[0] = sendsum; - count[0] = send_this_turn; - //Add the number of haloes that are sent this turn to the total number and set the size to the total - sendsum += send_this_turn; - dims[0] = sendsum; - //If the total size is larger than the initial set size for the data set extend it - if (sendsum > nhaloout) hdf5_status = H5Dset_extent(hdf5_dataset, dims); - //Get the data space of the data set - hdf5_filespace = H5Dget_space(hdf5_dataset); - //Select the hyperslab that corresponds to the current offset - hdf5_status = H5Sselect_hyperslab(hdf5_filespace, H5S_SELECT_SET, start, NULL, count, NULL); - //Set the size to the number of systems sent this turn - dims[0] = send_this_turn; - //Create a data space for the memory - hdf5_dataspace_memory = H5Screate_simple(1, dims, NULL); - //Write the data set - hdf5_status = H5Dwrite(hdf5_dataset, hdf5_dtype, hdf5_dataspace_memory, hdf5_filespace, H5P_DEFAULT, CommBuffer); - //Close the data space of the data set - H5Sclose(hdf5_filespace); - //Close the data space for the memory - H5Sclose(hdf5_dataspace_memory); + //Set the starting index and the count + start[0] = sendsum; + count[0] = send_this_turn; + //Add the number of haloes that are sent this turn to the total number and set the size to the total + sendsum += send_this_turn; + dims[0] = sendsum; + //If the total size is larger than the initial set size for the data set extend it + if (sendsum > nhaloout) hdf5_status = H5Dset_extent(hdf5_dataset, dims); + //Get the data space of the data set + hdf5_filespace = H5Dget_space(hdf5_dataset); + //Select the hyperslab that corresponds to the current offset + hdf5_status = H5Sselect_hyperslab(hdf5_filespace, H5S_SELECT_SET, start, NULL, count, NULL); + //Set the size to the number of systems sent this turn + dims[0] = send_this_turn; + //Create a data space for the memory + hdf5_dataspace_memory = H5Screate_simple(1, dims, NULL); + //Write the data set + hdf5_status = H5Dwrite(hdf5_dataset, hdf5_dtype, hdf5_dataspace_memory, hdf5_filespace, H5P_DEFAULT, CommBuffer); + //Close the data space of the data set + H5Sclose(hdf5_filespace); + //Close the data space for the memory + H5Sclose(hdf5_dataspace_memory); #endif - }//End HDF5 output - //Otherwise simply write the properties to an ascii file - else - { //Use the haloout struct for the CommBuffer - hp = (struct haloout *) CommBuffer; - //Go through all systems that were sent this turn - for (ihalo = 0; ihalo < send_this_turn; ihalo++) - { //Write a new line in the ascii file - fprintf(fd,"%f %e %f %f %f %f %f %f %f %f %e %f %f %e %f %f %f %f %f %f %hu %lld %lld %lld\n", - hp[ihalo].mh, - hp[ihalo].mdot, - hp[ihalo].mhp, - hp[ihalo].mhh, - hp[ihalo].amp, - hp[ihalo].amh, - hp[ihalo].r, - hp[ihalo].c, - hp[ihalo].l, - hp[ihalo].ms, - hp[ihalo].sfr, - hp[ihalo].icm, - hp[ihalo].mso, - hp[ihalo].sfro, - hp[ihalo].x, - hp[ihalo].y, - hp[ihalo].z, - hp[ihalo].u, - hp[ihalo].v, - hp[ihalo].w, - hp[ihalo].t, - hp[ihalo].hid, - hp[ihalo].did, - hp[ihalo].uid); - }//End loop through all haloes - }//End ascii output - }//End master task (writer task) - //Decrease the number of haloes left to send by the number sent this turn - left_to_send -= send_this_turn; - }//End halo increments - }//End loop through all tasks + } //End HDF5 output + //Otherwise simply write the properties to an ascii file + else + { //Use the haloout struct for the CommBuffer + hp = (struct haloout *) CommBuffer; + //Go through all systems that were sent this turn + for (ihalo = 0; ihalo < send_this_turn; ihalo++) + { //Write a new line in the ascii file + fprintf(fd,"%f %e %f %e %f %f %f %f %f %f %f %e %f %f %e %f %f %f %f %f %f %hu %lld %lld %lld", + hp[ihalo].mh, + hp[ihalo].mdot, + hp[ihalo].mhp, + hp[ihalo].mpd, + hp[ihalo].mhh, + hp[ihalo].amp, + hp[ihalo].amh, + hp[ihalo].r, + hp[ihalo].c, + hp[ihalo].l, + hp[ihalo].ms, + hp[ihalo].sfr, + hp[ihalo].icm, + hp[ihalo].mso, + hp[ihalo].sfro, + hp[ihalo].x, + hp[ihalo].y, + hp[ihalo].z, + hp[ihalo].u, + hp[ihalo].v, + hp[ihalo].w, + hp[ihalo].t, + hp[ihalo].hid, + hp[ihalo].did, + hp[ihalo].uid); + fprintf(fd,"\n"); + } //End loop through all haloes + } //End ascii output + } //End master task (writer task) + //Decrease the number of haloes left to send by the number sent this turn + left_to_send -= send_this_turn; + } //End halo increments + } //End loop through all tasks - }// nhaloout > 0 //If this is the master task (writer task) close the file (and other objects) - if(ThisTask == masterTask) + if (ThisTask == masterTask) { //For HDF5 output if (All.OutputFormat == 2) @@ -566,7 +572,7 @@ void output_halos(void) //Structure that stores all halo/galaxy properties that will be saved struct haloout { - float mvir,mdot,mpeak,ampeak,amhalf,rvir,c,l,x,y,z,u,v,w; + float mvir,mdot,mpeak,mpeakdot,ampeak,amhalf,rvir,c,l,x,y,z,u,v,w; unsigned short t; long long hid, did, uid; } *hp, hsend; @@ -600,17 +606,17 @@ void output_halos(void) distribute_file(All.NTaskPerUniverse, All.NumOutputFiles, 0, 0, All.NTaskPerUniverse - 1, &filenr, &masterTask, &lastTask); //Go through all output redshifts - for(i = 0; i < All.Noutputredshifts; i++) + for (i = 0; i < All.Noutputredshifts; i++) { //If the output format has been set to HDF5 if (All.OutputFormat == 2) { - //Check if HDF5 libraries are set + //If HDF5 libraries are set #ifdef HDF5_SUPPORT if (All.NumOutputFiles > 1) - sprintf(buf, "%s/haloes.S%d.%d.h5", All.OutputDir, Output_iscale[i], filenr); + sprintf(buf, "%s/haloes/haloes.S%d.%d.h5", All.OutputDir, Output_iscale[i], filenr); else - sprintf(buf, "%s/haloes.S%d.h5", All.OutputDir, Output_iscale[i]); + sprintf(buf, "%s/haloes/haloes.S%d.h5", All.OutputDir, Output_iscale[i]); //If not return #else if (ThisTask == 0) printf("%s Output format has been set to 2 but HDF5 support was not enabled.\n",All.startline); @@ -621,22 +627,22 @@ void output_halos(void) else { if (All.NumOutputFiles > 1) - sprintf(buf, "%s/haloes.S%d.%d", All.OutputDir, Output_iscale[i], filenr); + sprintf(buf, "%s/haloes/haloes.S%d.%d", All.OutputDir, Output_iscale[i], filenr); else - sprintf(buf, "%s/haloes.S%d.out", All.OutputDir, Output_iscale[i]); + sprintf(buf, "%s/haloes/haloes.S%d.out", All.OutputDir, Output_iscale[i]); } //Get number of groups ngroups = All.NumOutputFiles / All.NumFilesInParallel; - if((All.NumOutputFiles % All.NumFilesInParallel)) ngroups++; + if ((All.NumOutputFiles % All.NumFilesInParallel)) ngroups++; //For each group do... - for(j = 0; j < ngroups; j++) + for (j = 0; j < ngroups; j++) { //This task will be processed now - if((filenr / All.NumFilesInParallel) == j && MasterTask == 0) + if ((filenr / All.NumFilesInParallel) == j && MasterTask == 0) { //The masterTask opens the file - if(ThisTask == masterTask) + if (ThisTask == masterTask) { //If the hdf5 output format has been selected if (All.OutputFormat == 2) { //Check if the libraries have been included @@ -647,30 +653,33 @@ void output_halos(void) printf("%s Writing output file at z = %f: '%s' (file %d of %d)\n",All.startline, OutputRedshifts[i], buf, filenr+1, All.NumOutputFiles); //Specify data type hdf5_dtype = H5Tcreate(H5T_COMPOUND, sizeof(struct haloout)); - hdf5_status = H5Tinsert(hdf5_dtype, "Halo_mass", HOFFSET(struct haloout, mvir), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Halo_growth_rate", HOFFSET(struct haloout, mdot), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Halo_mass_peak", HOFFSET(struct haloout, mpeak), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Scale_peak_mass", HOFFSET(struct haloout, ampeak), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Scale_half_mass", HOFFSET(struct haloout, amhalf), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Halo_radius", HOFFSET(struct haloout, rvir), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Concentration", HOFFSET(struct haloout, c), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Halo_spin", HOFFSET(struct haloout, l), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "X_pos", HOFFSET(struct haloout, x), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Y_pos", HOFFSET(struct haloout, y), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Z_pos", HOFFSET(struct haloout, z), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "X_vel", HOFFSET(struct haloout, u), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Y_vel", HOFFSET(struct haloout, v), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Z_vel", HOFFSET(struct haloout, w), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Type", HOFFSET(struct haloout, t), H5T_NATIVE_USHORT); - hdf5_status = H5Tinsert(hdf5_dtype, "Halo_ID", HOFFSET(struct haloout, hid), H5T_NATIVE_LLONG); - hdf5_status = H5Tinsert(hdf5_dtype, "Desc_ID", HOFFSET(struct haloout, did), H5T_NATIVE_LLONG); - hdf5_status = H5Tinsert(hdf5_dtype, "Up_ID", HOFFSET(struct haloout, uid), H5T_NATIVE_LLONG); + k = 0; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_mass", HOFFSET(struct haloout, mvir), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_growth_rate", HOFFSET(struct haloout, mdot), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_mass_peak", HOFFSET(struct haloout, mpeak), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_growth_peak", HOFFSET(struct haloout, mpeakdot), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Scale_peak_mass", HOFFSET(struct haloout, ampeak), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Scale_half_mass", HOFFSET(struct haloout, amhalf), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_radius", HOFFSET(struct haloout, rvir), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Concentration", HOFFSET(struct haloout, c), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_spin", HOFFSET(struct haloout, l), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "X_pos", HOFFSET(struct haloout, x), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Y_pos", HOFFSET(struct haloout, y), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Z_pos", HOFFSET(struct haloout, z), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "X_vel", HOFFSET(struct haloout, u), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Y_vel", HOFFSET(struct haloout, v), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Z_vel", HOFFSET(struct haloout, w), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Type", HOFFSET(struct haloout, t), H5T_NATIVE_USHORT); k+=2; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_ID", HOFFSET(struct haloout, hid), H5T_NATIVE_LLONG); k+=8; + hdf5_status = H5Tinsert(hdf5_dtype, "Desc_ID", HOFFSET(struct haloout, did), H5T_NATIVE_LLONG); k+=8; + hdf5_status = H5Tinsert(hdf5_dtype, "Up_ID", HOFFSET(struct haloout, uid), H5T_NATIVE_LLONG); k+=8; //Specify file type + hdf5_ftype = H5Tcreate(H5T_COMPOUND, k); k = 0; - hdf5_ftype = H5Tcreate(H5T_COMPOUND, 4 * 14 + 2 + 8 * 3); hdf5_status = H5Tinsert(hdf5_ftype, "Halo_mass", k, H5T_IEEE_F32LE); k+=4; hdf5_status = H5Tinsert(hdf5_ftype, "Halo_growth_rate", k, H5T_IEEE_F32LE); k+=4; hdf5_status = H5Tinsert(hdf5_ftype, "Halo_mass_peak", k, H5T_IEEE_F32LE); k+=4; + hdf5_status = H5Tinsert(hdf5_ftype, "Halo_growth_peak", k, H5T_IEEE_F32LE); k+=4; hdf5_status = H5Tinsert(hdf5_ftype, "Scale_peak_mass", k, H5T_IEEE_F32LE); k+=4; hdf5_status = H5Tinsert(hdf5_ftype, "Scale_half_mass", k, H5T_IEEE_F32LE); k+=4; hdf5_status = H5Tinsert(hdf5_ftype, "Halo_radius", k, H5T_IEEE_F32LE); k+=4; @@ -691,7 +700,7 @@ void output_halos(void) //Otherwise open standard ascii file else { //Open file - if(!(fd = fopen(buf, "w"))) + if (!(fd = fopen(buf, "w"))) { //If not possible print to screen and abort printf("%s Can't open file `%s' for writing output.\n",All.startline,buf); endrun("file open error"); @@ -703,6 +712,7 @@ void output_halos(void) fprintf(fd,"#Halo_mass(%d)",task); task++; fprintf(fd," Halo_growth_rate(%d)",task); task++; fprintf(fd," Halo_mass_peak(%d)",task); task++; + fprintf(fd," Halo_growth_peak(%d)",task); task++; fprintf(fd," Scale_peak_mass(%d)",task); task++; fprintf(fd," Scale_half_mass(%d)",task); task++; fprintf(fd," Halo_radius(%d)",task); task++; @@ -731,6 +741,7 @@ void output_halos(void) fprintf(fd,"#Halo_mass: Current virial mass of the halo (log Msun).\n"); fprintf(fd,"#Halo_growth_rate: Growth rate of the halo (Msun/yr)\n"); fprintf(fd,"#Halo_mass_peak: Peak virial mass of the halo through its history (log Msun).\n"); + fprintf(fd,"#Halo_growth_peak: Growth rate of the halo when peak mass is reached (Msun/yr).\n"); fprintf(fd,"#Scale_peak_mass: Scale Factor when halo had its peak mass.\n"); fprintf(fd,"#Scale_half_mass: Scale Factor when halo first had half of its peak mass.\n"); fprintf(fd,"#Halo_radius: Current virial radius of the halo (kpc).\n"); @@ -752,7 +763,7 @@ void output_halos(void) //Determine number of systems that will be printed on each task (depending on redshift and mass) nhaloout = 0; for (ihalo = 0; ihalo < Nhalos; ihalo++) - if (H[ihalo].iscale == Output_iscale[i] && H[ihalo].gone == 0) nhaloout++; + if (H[ihalo].iscale == Output_iscale[i] && H[ihalo].gone == 0 && H[H[ihalo].impeak].mvir >= All.minmass) nhaloout++; if (nhaloout > 0) { @@ -775,7 +786,7 @@ void output_halos(void) hdf5_status = H5Pset_deflate(hdf5_properties, 9); // gzip compression level 9 hdf5_status = H5Pset_fletcher32(hdf5_properties); // Fletcher32 checksum on dataset //If filters could be set use them to create the data set - if(H5Pall_filters_avail(hdf5_properties)) + if (H5Pall_filters_avail(hdf5_properties)) hdf5_dataset = H5Dcreate(halo_file, path, hdf5_ftype, hdf5_dataspace_in_file, H5P_DEFAULT, hdf5_properties, H5P_DEFAULT); //Otherwise create the default data set else @@ -813,32 +824,32 @@ void output_halos(void) #endif //Now go through all tasks that are member of this group - for(task = masterTask, offset = 0; task <= lastTask; task++) + for (task = masterTask, offset = 0; task <= lastTask; task++) { //if this task is processed - if(task == ThisTask) + if (task == ThisTask) { //We need to get those haloes left_to_send = nhaloout; //Tell this to the other tasks - for(p = masterTask; p <= lastTask; p++) - if(p != ThisTask) + for (p = masterTask; p <= lastTask; p++) + if (p != ThisTask) MPI_Send(&left_to_send, 1, MPI_INT, p, TAG_NHALOS, MPI_COMM_WORLD); } //Reveive the number of haloes from active task else MPI_Recv(&left_to_send, 1, MPI_INT, task, TAG_NHALOS, MPI_COMM_WORLD, &status); //Now send the haloes in increments - while(left_to_send > 0) + while (left_to_send > 0) { //Compute how many haloes we can send this turn send_this_turn = left_to_send; //If they do not all fit into the CommBuffer set to the maximum - if(send_this_turn > maxlen) send_this_turn = maxlen; + if (send_this_turn > maxlen) send_this_turn = maxlen; //Sending task - if(ThisTask == task) + if (ThisTask == task) { //Use the haloout struct for the CommBuffer hp = (struct haloout *) CommBuffer; @@ -847,7 +858,7 @@ void output_halos(void) //Go through all systems until the maximum that can be sent while (ihalo < send_this_turn) { //If we are at the right scale factor and the halo is not gone send it to the master - if (H[offset].iscale == Output_iscale[i] && H[offset].gone == 0) + if (H[offset].iscale == Output_iscale[i] && H[offset].gone == 0 && H[H[ihalo].impeak].mvir >= All.minmass) { //Find scale factor when virial mass was half of current value amhalf = ScaleFactor[H[H[offset].impeak].iscale]; @@ -863,24 +874,25 @@ void output_halos(void) } //Set all halo properties - hsend.mvir = log10(H[offset].mvir*All.m_unit); - hsend.mdot = H[offset].mdotbary*All.m_unit/All.t_unit/All.f_baryon; - hsend.mpeak = log10(H[H[offset].impeak].mvir*All.m_unit); - hsend.ampeak = ScaleFactor[H[H[offset].impeak].iscale]; - hsend.amhalf = amhalf; - hsend.rvir = H[offset].rvir*All.x_unit*1.e3; - hsend.c = H[offset].c; - hsend.l = H[offset].lambda; - hsend.x = H[offset].pos[0]; - hsend.y = H[offset].pos[1]; - hsend.z = H[offset].pos[2]; - hsend.u = H[offset].vel[0]; - hsend.v = H[offset].vel[1]; - hsend.w = H[offset].vel[2]; - hsend.t = H[offset].type; - hsend.hid = (long long)(H[offset].haloid)-1; - hsend.did = (long long)(H[offset].descid)-1; - hsend.uid = (long long)(H[offset].upid)-1; + hsend.mvir = log10(H[offset].mvir*All.m_unit); + hsend.mdot = H[offset].mdotbary*All.m_unit/All.t_unit/All.f_baryon; + hsend.mpeak = log10(H[H[offset].impeak].mvir*All.m_unit); + hsend.mpeakdot = H[H[offset].impeak].mdotbary*All.m_unit/All.t_unit/All.f_baryon; + hsend.ampeak = ScaleFactor[H[H[offset].impeak].iscale]; + hsend.amhalf = amhalf; + hsend.rvir = H[offset].rvir*All.x_unit*1.e3; + hsend.c = H[offset].c; + hsend.l = H[offset].lambda; + hsend.x = H[offset].pos[0]; + hsend.y = H[offset].pos[1]; + hsend.z = H[offset].pos[2]; + hsend.u = H[offset].vel[0]; + hsend.v = H[offset].vel[1]; + hsend.w = H[offset].vel[2]; + hsend.t = H[offset].type; + hsend.hid = (long long)(H[offset].haloid)-1; + hsend.did = (long long)(H[offset].descid)-1; + hsend.uid = (long long)(H[offset].upid)-1; //Write this system to the CommBuffer *hp++ = hsend; @@ -893,13 +905,13 @@ void output_halos(void) } //Receive haloes - if(ThisTask == masterTask && task != masterTask) MPI_Recv(CommBuffer, (sizeof(struct halo)) * send_this_turn, MPI_BYTE, task, TAG_HDATA, MPI_COMM_WORLD, &status); + if (ThisTask == masterTask && task != masterTask) MPI_Recv(CommBuffer, (sizeof(struct halo)) * send_this_turn, MPI_BYTE, task, TAG_HDATA, MPI_COMM_WORLD, &status); //Send haloes - if(ThisTask != masterTask && task == ThisTask) MPI_Ssend(CommBuffer, (sizeof(struct halo)) * send_this_turn, MPI_BYTE, masterTask, TAG_HDATA, MPI_COMM_WORLD); + if (ThisTask != masterTask && task == ThisTask) MPI_Ssend(CommBuffer, (sizeof(struct halo)) * send_this_turn, MPI_BYTE, masterTask, TAG_HDATA, MPI_COMM_WORLD); //Collect haloes and write to file - if(ThisTask == masterTask) + if (ThisTask == masterTask) { //If the output format is hdf5 if (All.OutputFormat == 2) @@ -936,10 +948,12 @@ void output_halos(void) //Go through all systems that were sent this turn for (ihalo = 0; ihalo < send_this_turn; ihalo++) { //Write a new line in the ascii file - fprintf(fd,"%f %e %f %f %f %f %f %f %f %f %f %f %f %f %hu %lld %lld %lld\n", + fprintf(fd, + "%f %e %f %f %f %f %f %f %f %f %f %f %f %f %f %hu %lld %lld %lld", hp[ihalo].mvir, hp[ihalo].mdot, hp[ihalo].mpeak, + hp[ihalo].mpeakdot, hp[ihalo].ampeak, hp[ihalo].amhalf, hp[ihalo].rvir, @@ -955,6 +969,7 @@ void output_halos(void) hp[ihalo].hid, hp[ihalo].did, hp[ihalo].uid); + fprintf(fd,"\n"); }//End loop through all haloes }//End ascii output }//End master task (writer task) @@ -966,7 +981,7 @@ void output_halos(void) }// nhaloout > 0 //If this is the master task (writer task) close the file (and other objects) - if(ThisTask == masterTask) + if (ThisTask == masterTask) { //For HDF5 output if (All.OutputFormat == 2) @@ -1086,9 +1101,9 @@ void output_mainbranch(void) //Check if HDF5 libraries are set #ifdef HDF5_SUPPORT if (All.NumOutputFiles > 1) - sprintf(buf, "%s/mainbranches.S%d.%d.h5", All.OutputDir, All.MainBranch_iscale, filenr); + sprintf(buf, "%s/mainbranches/mainbranches.S%d.%d.h5", All.OutputDir, All.MainBranch_iscale, filenr); else - sprintf(buf, "%s/mainbranches.S%d.h5", All.OutputDir, All.MainBranch_iscale); + sprintf(buf, "%s/mainbranches/mainbranches.S%d.h5", All.OutputDir, All.MainBranch_iscale); //If not return #else if (ThisTask == 0) printf("%s Output format has been set to 2 but HDF5 support was not enabled.\n",All.startline); @@ -1099,22 +1114,22 @@ void output_mainbranch(void) else { if (All.NumOutputFiles > 1) - sprintf(buf, "%s/mainbranches.S%d.%d", All.OutputDir, All.MainBranch_iscale, filenr); + sprintf(buf, "%s/mainbranches/mainbranches.S%d.%d", All.OutputDir, All.MainBranch_iscale, filenr); else - sprintf(buf, "%s/mainbranches.S%d.out", All.OutputDir, All.MainBranch_iscale); + sprintf(buf, "%s/mainbranches/mainbranches.S%d.out", All.OutputDir, All.MainBranch_iscale); } //Get number of groups ngroups = All.NumOutputFiles / All.NumFilesInParallel; - if((All.NumOutputFiles % All.NumFilesInParallel)) ngroups++; + if ((All.NumOutputFiles % All.NumFilesInParallel)) ngroups++; //For each group do... - for(j = 0; j < ngroups; j++) + for (j = 0; j < ngroups; j++) { //This task will be processed now - if((filenr / All.NumFilesInParallel) == j && MasterTask == 0) + if ((filenr / All.NumFilesInParallel) == j && MasterTask == 0) { //The masterTask opens the file - if(ThisTask == masterTask) + if (ThisTask == masterTask) { //If the hdf5 output format has been selected if (All.OutputFormat == 2) { //Check if the libraries have been included @@ -1125,30 +1140,31 @@ void output_mainbranch(void) printf("%s Writing main branch for galaxies at z = %3.2f: '%s' (file %d of %d)\n",All.startline, All.mainBranchRedshift, buf, filenr+1, All.NumOutputFiles); //Specify data type hdf5_dtype = H5Tcreate(H5T_COMPOUND, sizeof(struct mbout)); - hdf5_status = H5Tinsert(hdf5_dtype, "Scale_factor", HOFFSET(struct mbout, a), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Halo_mass", HOFFSET(struct mbout, mh), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Halo_growth_rate", HOFFSET(struct mbout, mdot), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Halo_mass_peak", HOFFSET(struct mbout, mhp), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Scale_peak_mass", HOFFSET(struct mbout, amp), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Halo_radius", HOFFSET(struct mbout, r), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Concentration", HOFFSET(struct mbout, c), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Halo_spin", HOFFSET(struct mbout, l), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Stellar_mass", HOFFSET(struct mbout, ms), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "SFR", HOFFSET(struct mbout, sfr), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Intra_cluster_mass", HOFFSET(struct mbout, icm), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "X_pos", HOFFSET(struct mbout, x), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Y_pos", HOFFSET(struct mbout, y), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Z_pos", HOFFSET(struct mbout, z), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "X_vel", HOFFSET(struct mbout, u), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Y_vel", HOFFSET(struct mbout, v), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Z_vel", HOFFSET(struct mbout, w), H5T_NATIVE_FLOAT); - hdf5_status = H5Tinsert(hdf5_dtype, "Type", HOFFSET(struct mbout, t), H5T_NATIVE_USHORT); - hdf5_status = H5Tinsert(hdf5_dtype, "Halo_ID", HOFFSET(struct mbout, hid), H5T_NATIVE_LLONG); - hdf5_status = H5Tinsert(hdf5_dtype, "Desc_ID", HOFFSET(struct mbout, did), H5T_NATIVE_LLONG); - hdf5_status = H5Tinsert(hdf5_dtype, "Up_ID", HOFFSET(struct mbout, uid), H5T_NATIVE_LLONG); + k = 0; + hdf5_status = H5Tinsert(hdf5_dtype, "Scale_factor", HOFFSET(struct mbout, a), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_mass", HOFFSET(struct mbout, mh), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_growth_rate", HOFFSET(struct mbout, mdot), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_mass_peak", HOFFSET(struct mbout, mhp), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Scale_peak_mass", HOFFSET(struct mbout, amp), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_radius", HOFFSET(struct mbout, r), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Concentration", HOFFSET(struct mbout, c), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_spin", HOFFSET(struct mbout, l), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Stellar_mass", HOFFSET(struct mbout, ms), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "SFR", HOFFSET(struct mbout, sfr), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Intra_cluster_mass", HOFFSET(struct mbout, icm), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "X_pos", HOFFSET(struct mbout, x), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Y_pos", HOFFSET(struct mbout, y), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Z_pos", HOFFSET(struct mbout, z), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "X_vel", HOFFSET(struct mbout, u), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Y_vel", HOFFSET(struct mbout, v), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Z_vel", HOFFSET(struct mbout, w), H5T_NATIVE_FLOAT); k+=4; + hdf5_status = H5Tinsert(hdf5_dtype, "Type", HOFFSET(struct mbout, t), H5T_NATIVE_USHORT); k+=2; + hdf5_status = H5Tinsert(hdf5_dtype, "Halo_ID", HOFFSET(struct mbout, hid), H5T_NATIVE_LLONG); k+=8; + hdf5_status = H5Tinsert(hdf5_dtype, "Desc_ID", HOFFSET(struct mbout, did), H5T_NATIVE_LLONG); k+=8; + hdf5_status = H5Tinsert(hdf5_dtype, "Up_ID", HOFFSET(struct mbout, uid), H5T_NATIVE_LLONG); k+=8; //Specify file type + hdf5_ftype = H5Tcreate(H5T_COMPOUND, k); k = 0; - hdf5_ftype = H5Tcreate(H5T_COMPOUND, 4 * 17 + 2 + 8 * 3); hdf5_status = H5Tinsert(hdf5_ftype, "Scale_factor", k, H5T_IEEE_F32LE); k+=4; hdf5_status = H5Tinsert(hdf5_ftype, "Halo_mass", k, H5T_IEEE_F32LE); k+=4; hdf5_status = H5Tinsert(hdf5_ftype, "Halo_growth_rate", k, H5T_IEEE_F32LE); k+=4; @@ -1183,7 +1199,7 @@ void output_mainbranch(void) //Otherwise open standard ascii file else { //Open file - if(!(fd = fopen(buf, "w"))) + if (!(fd = fopen(buf, "w"))) { //If not possible print to screen and abort printf("%s Can't open file `%s' for writing output.\n",All.startline,buf); endrun("file open error"); @@ -1226,7 +1242,7 @@ void output_mainbranch(void) fprintf(fd,"#This file contains all main branches with a mass of %s at z = %4.2f\n",All.output_mass_mb,All.mainBranchRedshift); //Add which mass has been used for selecting the trees if (All.MainBranchMassType == 1) fprintf(fd,"#Stellar mass has been used for this selection.\n"); - else fprintf(fd,"# Halo mass has been used for this selection.\n"); + else fprintf(fd,"#Halo mass has been used for this selection.\n"); fprintf(fd,"#Halo_mass: Current virial mass of the halo (log Msun).\n"); fprintf(fd,"#Halo_growth_rate: Growth rate of the halo (Msun/yr)\n"); fprintf(fd,"#Halo_mass_peak: Peak virial mass of the halo through its history (log Msun).\n"); @@ -1315,9 +1331,9 @@ void output_mainbranch(void) imb = 0; //Now go through all tasks that are member of this group - for(task = masterTask; task <= lastTask; task++) + for (task = masterTask; task <= lastTask; task++) { //If this task is processed - if(task == ThisTask) + if (task == ThisTask) { //If we're beyond the master task, send number of root haloes to master task if (ThisTask != masterTask) MPI_Ssend(&Nroot, 1, MPI_INT, masterTask, TAG_NTREES, MPI_COMM_WORLD); //Go through all root haloes on this task @@ -1327,7 +1343,8 @@ void output_mainbranch(void) //Set flag for presence at selected redshift to 1 present_at_z = 1; //Set iprog to root halo - iprog = roothalos[iroot]; //Loop until main branch progenitor is at selected redshift + iprog = roothalos[iroot]; + //Loop until main branch progenitor is at selected redshift while (H[iprog].iscale > All.MainBranch_iscale) { //If progenitor does not exist if (H[iprog].iprog < 0) @@ -1385,7 +1402,7 @@ void output_mainbranch(void) Nmainbranch++; } //End loop through main branch //If we're beyond the master task, send number of haloes in main branch and then send full main branch - if(ThisTask != masterTask) + if (ThisTask != masterTask) { MPI_Ssend(&Nmainbranch, 1, MPI_INT, masterTask, TAG_NHALOS, MPI_COMM_WORLD); MPI_Ssend(mainbranch, (sizeof(struct mbout)) * Nmainbranch, MPI_BYTE, masterTask, TAG_HDATA, MPI_COMM_WORLD); @@ -1410,7 +1427,7 @@ void output_mainbranch(void) hdf5_status = H5Pset_deflate(hdf5_properties, 9); // gzip compression level 9 hdf5_status = H5Pset_fletcher32(hdf5_properties); // Fletcher32 checksum on dataset //If filters could be set use them to create the data set - if(H5Pall_filters_avail(hdf5_properties)) + if (H5Pall_filters_avail(hdf5_properties)) hdf5_dataset = H5Dcreate(mb_file, path, hdf5_ftype, hdf5_dataspace_in_file, H5P_DEFAULT, hdf5_properties, H5P_DEFAULT); //Otherwise create the default data set else @@ -1430,8 +1447,9 @@ void output_mainbranch(void) fprintf(fd,"#%03d_%07d\n",i,imb); //Then write all systems in the main branch for (ihalo = 0; ihalo < Nmainbranch; ihalo++) + { fprintf(fd, - "%f %f %e %f %f %f %f %f %f %e %f %f %f %f %f %f %f %hu %lld %lld %lld\n", + "%f %f %e %f %f %f %f %f %f %e %f %f %f %f %f %f %f %hu %lld %lld %lld", mainbranch[ihalo].a, mainbranch[ihalo].mh, mainbranch[ihalo].mdot, @@ -1453,6 +1471,8 @@ void output_mainbranch(void) mainbranch[ihalo].hid, mainbranch[ihalo].did, mainbranch[ihalo].uid); + fprintf(fd,"\n"); + } //And separate with an empty line fprintf(fd,"\n"); } //Done writing the main branch for master task @@ -1501,7 +1521,7 @@ void output_mainbranch(void) hdf5_status = H5Pset_deflate(hdf5_properties, 9); // gzip compression level 9 hdf5_status = H5Pset_fletcher32(hdf5_properties); // Fletcher32 checksum on dataset //If filters could be set use them to create the data set - if(H5Pall_filters_avail(hdf5_properties)) + if (H5Pall_filters_avail(hdf5_properties)) hdf5_dataset = H5Dcreate(mb_file, path, hdf5_ftype, hdf5_dataspace_in_file, H5P_DEFAULT, hdf5_properties, H5P_DEFAULT); //Otherwise create the default data set else @@ -1521,8 +1541,9 @@ void output_mainbranch(void) fprintf(fd,"#%03d_%07d\n",i,imb); //Then write all systems in the main branch for (ihalo = 0; ihalo < Nmainbranch; ihalo++) + { fprintf(fd, - "%f %f %e %f %f %f %f %f %f %e %f %f %f %f %f %f %f %hu %lld %lld %lld\n", + "%f %f %e %f %f %f %f %f %f %e %f %f %f %f %f %f %f %hu %lld %lld %lld", mainbranch[ihalo].a, mainbranch[ihalo].mh, mainbranch[ihalo].mdot, @@ -1544,6 +1565,8 @@ void output_mainbranch(void) mainbranch[ihalo].hid, mainbranch[ihalo].did, mainbranch[ihalo].uid); + fprintf(fd,"\n"); + } //And separate with an empty line fprintf(fd,"\n"); } //Done writing the main branch for master task @@ -1568,7 +1591,7 @@ void output_mainbranch(void) }//Done with this mass bins - go to next //If this is the master task (writer task) close the file (and other objects) - if(ThisTask == masterTask) + if (ThisTask == masterTask) { //For HDF5 output if (All.OutputFormat == 2) { //Check if libraries have been included @@ -1598,3 +1621,4 @@ void output_mainbranch(void) } #endif + diff --git a/src/read_trees.c b/src/read_trees.c index 4fa8eaf..3a549b9 100644 --- a/src/read_trees.c +++ b/src/read_trees.c @@ -732,7 +732,7 @@ void empty_read_buffer(int offset, int hc) H[offset + n].a = h[n].a; H[offset + n].mvir = h[n].mvir/All.h_100/All.m_unit;; - H[offset + n].rvir = h[n].rvir/All.h_100/All.x_unit/1000.0; + H[offset + n].rvir = h[n].rvir/All.h_100/All.x_unit/1000.0*h[n].a; H[offset + n].c = h[n].c; H[offset + n].lambda = h[n].lambda; diff --git a/src/setup.c b/src/setup.c index 05cca9b..114e584 100644 --- a/src/setup.c +++ b/src/setup.c @@ -45,7 +45,7 @@ void read_parameterfile(char *fname) #define MAXTAGS 300 ///< Maximum number of entries the parameter file can contain FILE *fd; - char buf[200], buf1[200], buf2[200], buf3[200]; + char buf[NSTRING], buf1[NSTRING], buf2[NSTRING], buf3[NSTRING]; int i, j, nt; int id[MAXTAGS]; void *addr[MAXTAGS]; @@ -534,7 +534,7 @@ void read_parameterfile(char *fname) void setup(void){ int i; - char buf[500], tmp1[NSTRING], *tmp2; + char buf[NSTRING], tmp1[NSTRING], *tmp2; //Print a line if (ThisTask == 0) printf("%s\n",All.fullline); @@ -653,6 +653,27 @@ void setup(void){ { if (All.verbose >= VERBOSE_MIN) printf("%s Creating output directory '%s'\n",All.startline,All.OutputDir); } +#ifdef WRITE_GALAXY_CATALOG + sprintf(buf, "%s/galaxies", All.OutputDir); + if (mkdir(buf, 02755) == 0) + { + if (All.verbose >= VERBOSE_MIN) printf("%s Creating galaxy catalogue directory '%s'\n",All.startline,buf); + } +#endif +#ifdef WRITE_HALO_CATALOG + sprintf(buf, "%s/haloes", All.OutputDir); + if (mkdir(buf, 02755) == 0) + { + if (All.verbose >= VERBOSE_MIN) printf("%s Creating halo catalogue directory '%s'\n",All.startline,buf); + } +#endif +#ifdef WRITE_MAINBRANCH + sprintf(buf, "%s/mainbranches", All.OutputDir); + if (mkdir(buf, 02755) == 0) + { + if (All.verbose >= VERBOSE_MIN) printf("%s Creating mainbranch directory '%s'\n",All.startline,buf); + } +#endif } //Open all log files diff --git a/user-guide.pdf b/user-guide.pdf index ff573aa53536d43ccd1d821954e9b983fdd5418c..eaa7164d28701eca6097957893df84e1968a5fc9 100644 GIT binary patch delta 5185 zcmV-H6u#@TtsJhc9DuX|vjdahy()hiS##sI`Q5)lGkqvM6Jv>k1elq=q}e3fHk)+f zw=~nz5+!p?i4IA{>;C$F2f#z3WbaxteUU&MAHMql5HB{1c=7$!H@8>szDwhUlp;$K zwYXg`l2oLbUO4PBc5%B}{2D29{oC!|7A8)G%nYog-fzV9jkZbj^L1t--~NBGZo8tv z*;nsxuYSLhFd4x2P-82MG>hf->et`m#R|s$wupt!+~SF5ZWqQ$;WB7ei(jt(c@mi*EdG1Xq9(4x9+O?`bODk zc?~}gEq5;SogbyKf11=2Yt5vb$W>@?Fdp}m}{7Ah%82MMMO;rFhZK9%& zbytM(yz6H+zH^J%?OQ*vtGfw@5~i+4R%;1=GTJ_Zq@wN`Fn#7;k`I42<1l+ovcnES zjFuwG*n=ofMi}L<_#+R~`r)QnwfE(^<8FY=jl8OOd|5Y5QRQ7(@BEbgZi!zWrGd(T zAhLKy7kEdNMs2>OmIc(3cvJp?v=%$=LO$Z8VRn}{9go$(o;5f7VWtiMNx{05WD9A9 zl8N!UWD*fuk3))9V6%Vdha5FXCt0*DnoWT;r@)h@D7crQwkNbxOJS5_h#yet03xVU z$viUFu|`2}eb*MvAIKrnHR##oRSsl)dgf2m>95chZ1kQ`;q}1r5@k*>{0^YW#A=WQ z19PytdEmiq-l2n`_B|1YA^3lIq>uWs-f966|G;g1;%)B(8z!zs_nc=l_WbAkwX@~4 zLOLDd#4V7wE~|gyBbvjEqJjT*aS~<&dlRFHy~qgAYMANVlppPWvng5P{syp9gyk!$p*013lKlbc33t;ngSz}HsH#Fn;}VSVbh)g>$+MM4G?Bz9DP_H z;N_){&AcfN2?`{Zo~|7KB-qGqQ+W<4V0WhVltaN8&e4CMgpX)@B23D_B1M{*u4<4sMHafN}8|-U42udq%qh*wc|Bk=C}b-sP=tp+JzQ;JBXdOG+p10|^uW zarP*rLCcC1qRE6W`J;C!yA|tY@Q*4nL7|f#1MDne<|o)#rP)bgN-KqR?v=tcke%7n z>WB!SI?iC!&TcN9NYfr5ZYyRt)~U#173R&#GX;P7Boj#TBqaSJBr}1ez7P_IN0Igl zo(yU73_RckxZ439WI72446RRjL?eY#@@0=`gw-=15u|g=nI!H(J_7!CHUEQfyQW-* zoJ>dc2pOJ97DsmHV|n9euKBFhgxU;w=u{}1{_)Y*Rg`3N{-n}Tkph3K8tH!8Rdr784*NW&-(bNRupB$9>m5e*Ny<)6ApTNF62ho8+SVD>O1^+o+Ld`J~Xutk4;BK z4R&&!35^1|?RblXN*6R2=L?A;F>2?km_1bv zpp%i$hQUCWJ2-aXQ7YqbjCSwH+gYokpLTwEyZ79hq0`)iWf<~QaP0+KoUcKE4nNxR zf5$xd=(IiPOC;m0Fv44m)*W$Yt+7)fPyW5xc1d7&e}Ic z075~0430^~7FQX{Dt^I79YG;$DCGh}m?E32qT_mo%AwEG;fVRw$-v@y)` zf*ByDXxDrK(dFn9IVFGORc-P69@NNBB+B9cn8`3TAe!@^M{~xRQXCtFa2uRF1pCq_ z5;V?y|N7)u~;7+%TF#4AOV3y6-sm?QfHXue>&!&%Y>211+jkzfT=Wx*t~ASSW)(hb-WOQ`wZPTO{Y1 z99c9&i&9$1i;4&#>SGx}_N(mcX;FB(b|j**W4VyyJEp}sVV%AtEdVi_Lc{?JML(3} zuO|kU8>6Pgysv*KDFU%{?J+Gn$%Jy3lnX*M?no}u0@Gek3aF{b%z<3~3OTALWyX(c zC-fId<;amD2=Li&a->#*?WX6kz>7YUc@V^egYzjZVl-EWM=AdLG(PS9zHG~B^_V48 z(@nPg*dN5q`=?8gDjgoe+UdqVvpwTYF0r*zU;>9UYjuCd?k}=yW*}HkHNDlKPjkg*@iG{$|O>T)$8nT()zG-;ghzu+h(^=GaETIcChWADt}b~+rkpE20ars~*W z0Ys|aR-GsS^>SbbJYO_dEg^5o{Ck_z=K3WY&Kv7^L%f{(@a<`1eLI=GsgAe(@NGv2 zNB&^(jLE(of!<)p2!O#-ElkIl@SmD92K;u`2lIa&qwu2nzMb^mc&EJgr2D5hqZ3?o zpXZJ5GjTdLJpf>0`6knoFyHFS7@q{2(Ej2;pEkZP+j8Fg##_=mrjy2Z2||1G8*j@; zlV_~%5mAvXTq z=Hu{W`{skG5%Fj0GW)r@uoF6X&r0Fen%bnp$ zcftE{bMYcw!{}ga}TUcEWnChCB8)VqS+P2%0s5}ewU7(B*I62R~7v*aY=(h9zK6F zTGqo|B|q5HNg9w&lBmgd8@jJawSVWEq@Fp+=t(E`wnfLigiaF9{ITOOuR3yvI;NLl zvszcxc&{;TC8sn|>P#fV<$k&L`}^ytIpzaTUC;Y}{Pgkr58wVW%X};J@FuB#q&p(M z@HZLwHhFMgPOmRa#IGgr!i6seixq!9AaNP%hcUgDn&i~09|-}_xc&iN2+v2M4Z>58 z-{Gw)?lgKOJ8SvYg4+!MhG0aT{|o>9vfuu@WC_p&AghIrIZ32O#HxR2R#g4fbNdy) z8YxBLS<#!Z&VbUy9KJu@qFOs&%lR8-Y#bDI&KGr#-%i zKQCMfk0JMa<9*c>tvsyqH>~|UZuA-L>@69+P%t|2SH$#2ou2IXz_9E5ZQ3Su4L51w z1UdXeO#UCE3;ocCixB~bixC2cixC5dixC8eixCC4ixCE}^(rwjFf=zaH!CnQConK4 zDGD!5Z)8MabY&nYL^?7sF*z|XGBGhQG&eIhw`B(h;s`=SGdMU!Ml(V=HA6&1MMf|- zIYC22H#s;oG%!L%G&DXSJVit^I5N}&aMp@kj_6ey)F(3V75y1{60VIt9uK{h5D(L^+3h`4b< zFoeWMqK2q(W#ZPzN|P~~xY0!K@4xxo@B7c3bLY(Wm~qaz|03riXD;_@;m3u^$VEvt znt=#y3att^1=48KXqC7bkVQL$R)L!XMQDEuXyv%Ypakst>DgpLuln_7jO@QBWT1Z zUdEjTd9)g|^SDRBF|=m1bGUP09<6^1?JVwba00C#Z3%Y)BxoaOXK+t~MYMfrr*Th# z(`ftAa5GwBu-cl@_^-b_y-Q)d#MiX-~BbZH{(L+o^q402CR8Y}1@8 zSAuNt*cdnnlqGw>IG6!cZG4Eyq)$+FlW;*bO^OCpGD#g&zm!2x)lxV?wMu`f1yv|@ z8B~{4V^B3xgMrQvsJu`EsJBoDsIt%isIAZhsHo5asH4yZx!s4lvM; zHEM39TAx4(g>FFMgdKo#34MSf2?K!A2!nuv2t$A}2%YB2Z|GH2Ko|w&eHaI1d)NcW z@j$`_GBJ#pt6I`y`+)L52{L~K_M>UJipOq1K?H?j8Y>S-)^H4vqG2A85#flr)_Z!r zGG`d*J5FHrQ&wllmkKQHGdUE=qGAZ-IV9!RNh&fMfP8YLZrvQD2x zUWGMt`Oh>;)&ge}&R0)_KqYn@{xZ{~mg{Zv*lNSnKM zTW$vAKq)9O_uO3##9$jJ2l_y@>oRjUKGc*-P|I4>52nCTv}RCm?&ed?Q%6?|Zv|@7 z?Le)$6R6R)n7dW4C)A(2dGiDo<($6&*1=_P6{r^vgI;rQzN}{z6g6aZV6|NJ&O)l( zc|}upgE64SIAQLCpaFmNQ#H_oV9MMl6&g@4nZ@E$^V8-YY-nHx&?575VAkB@Z#6(W z%xgD}gOgz1+}FbzrQzilz*(?p?#Guja0*-mh1~dgLE3q61+1ET`h&DJ7NnMIxP=t> zzdU-86e!g8q@+} z>Fdo#p4AkB>>B}*_MFb1Bh|N=%@k`2G55sW6LZfw9kWB8K=~fCWK{D=2Hy*auP2s% zC*WB0yUbSK*E~Y&iM}VGo{;%5vuoN(Ps}~(}J14wZ#O=&zXO{@s%{4U}^xGtOk&ny<4PFBJFih>0}~!e#z|qBaISnucK%UESo(X z(g5l03A|qg?+M2YscUxAEWxxMP)($WS=7Adg3Crd5Z3#D-gyxpi zfqKvYbmdBH0$gm87SIaX0LNm|0XhNqoJ1E=x>=%YDF2)E0%DYK?MViJE~JA%*HBKd vgcB#>)J-_N65U0)U?k&~%iedEEC0p)4=MY7jhBiR2OtSCI0_{tMNdWwNg3;j delta 5175 zcmV-76v*qYtsJwh9DuX|vjYJ&li|H7e;QqL_Ut4tw2PDMf0;>dzm?gT0l)$Cd>&~3P(EJY7_E|3(HGphtDEZ^ z<&xz!{Cw=Vb6IS|D2JCj38yv2q(mn)T7hzfx`L5`^ij zN|m#o{u(HK;$GSxe{7}^dQE$W z9Yh#SMUt@xNuG=_$#3|hh+xBTTdul?YTa`;VCF_q*F3&#+P18VzG}7sWWQbFS0L%Y zWh4+qJfjP|qRNu4*ig#>X*RvDen(!*Eq9?D5$FKj7j4gD4TxvW&0(0i14L4=F4=q` zolw%6pi9Pz)CC$cv;vnWf8Q6VLE7fYrfly^F{GtB60-QFo#BSmCB)>tI6gf3c&FW;d)%p#1~4 zB11qgn1wesxqS!|qZE&pSprA6g)pnALNkDrU6>lTg}M*4;!0bzkJT9st7t1#LJq=% zPl%eD1qXqR!4Kd^$Ly8M3du$OG%QpFP}VbOkX7RHEog7g8WTBW2o{NEKEk9TLX#p* z{j_#oNUvkqxCJiPe^p(6KTqaZ-Z0n+=aJi7h;ec92TDH~+s zf=sxJ%pDsQ($w_MOu0PW?0f%xNH!>x~W%X3yc|= zChylrbOj~AX3>_%3`G{Jz*m8P5=>@$Uk3>(5O=QiR6@}=f6nosgqQZkyW@p&11sb+ zgJLz?73!Zepl2Y7@Sq)32 z(7@_=AKv~X({zr5ObR>IFXwO@SiLe0!pr9zdYZ+((W&FT!M?7C7`5UyT1JKRZ}~g; zvA){X0XitxfBZ(sRH)Rycqbs0$_8kPE@0hVpp;2=2k}DZUD1UWZV1X$%H;8(qU_=! zvOw_=p+^x8T25pTTqbfUcEP1=SFD#d_y%Q7RA@U;fSnaU{s9ZCEI+%M(n;aGf8k~t zNY>nWbwmPC9cL)&W*3)UWZ8famld-M>rCXSig2_0f0@Dl*a;?i7L)!YCNsgLz7P{e zOOXub?^x2r48X>W3n8aR)M8HCW@ zpkBKSe-g_6q1oS~8Du8GP>lcV7FN~o)oNeV+!taFv>;+*+Dz>sCBY9~nlXfB_-lnx zT!W-yUV)u~ISMU?(-XqZIl?M@vWXLM`@VsMs#-!4{4|$ntRSW@QAEZE5t$ri>`!G3 z10-XBJdj~infwqFPsG8A=K+Q%mVfkZv83RJe`D*m*>J1izfnfWTMu(W7BU7}1>AK> zsnb0~cZMO-X4pfgLdo?HyHI9P&dvFg%0@K{{H+?~`*B-01-U!y3zUx2_$9~cRhaci zUiq%8dNP+QEQYbr{kVv&zn>JgDzhT9d9;p)zTbWI_U+^2qu6&?IEi|*7pwAZ$vF(8 zf8{dy#DtxK5o^b7P16e3S;^1ct?SnzZV5YlW!TY`97ZbML`L}0QCL@|h1GrL)Bc^Os_ zXNUY7vsosH_RSe~VgL~`6g#RNy_C>b;Yyq@Z)Kc^w?zAHX%H= zJ=Hka$z>-rN|d(eB`j4kXf9r|)gBU~Zmy`=Q}sY@3uy}<%Vx8o8VcI~*#}{8A)~>q zRJa{~j-VtNP8o(Gg&1vOu5^h+sy%;a3YA8pd1he|cbCvWGRN`Bmi>D@j^YlcsA6Q^OwQ_;1TY;mCm0Y1WL z$Nw?YG*(Z~M{tnNJE*blBbM@Df6;b0RkBnT%4EtDo`v0biG@~NE~-7YLE+Z>n%h+J z!$0FL2caCQq+@nL+qw$Mv+a^>M1?9K7M)IwiU4G@-7oun%hsqvH30=%@$ouTvk0|s z%ZP-6_!u3tiY=}(R8T@}F|a(C4t`MLxD04^b+L?#qgdw(aHJ{CQ9P2Qe^>1=n@%yo z9YQ>QpB&)ftl#!*NgwgNj-cWp7_|2t*FG{HU8M&--9_*_B zhS7tG(4@e~y6E~jMqw@TbVW<&1z9*Ri198P6&M-UKh3`94+!Xc$}0%$p+IWbo#Qtx z3{YWX4j41EYd(SK3-pOLe;l;ntgifR4{8(uiF5cr6*7QELJR&2WX?EKiesY??xK^2 zU|)qqg2q|yU!6VN*3ya`DzcX!ZeJP&9wng_f}|| zy_|Zg|D2hU|23Vdg4ej#k_-$Y_%>AOQ?8_=^P$Ueg$Q+)AG2i~l)<0F7RbS?e9G1> zvhzfaEScd&DJ|qhw+J!n(>;RhSGkw-qVRO>WQ(ev_k|qaDKFj&=k+Cd0gU++BMvwy z`k^F$IWsWb7&T?)e_c&U5o}A}o${h>E|kCIz92^9PxeJR5ZcRG0W}r5Iog-MLXH|( znen6A3;jjAa^gr41^Db&IZ`LVcC+WHzX@jK*8XKeJdxjD5_aE`fi zud-32Adb1q*l32n+MnIWXro`Ymep`dJ&j zgrK91P6U0%e;E#jdI@$1L!GRAY^cgRGF1HM+w5?MgB^6ohqqos8IQd)SKZm+u>Bc> z{cNsI4Hi&j>UGtL0#u(4>_F!y%~eatTk`Py=DfLn$%6C7I$jVj=iYyF-dJDH=3uJh zW#50()4@?VSbWB0Ur)ebuww$iV5t_SV@&wZ%rgf3e|phJ^Bt4$N%MU@>x1!5Wgl1% zFLOpGxadBgH@?rL>D2T9fl1}7Oiu#5)0Z(m2{xg_#eqI=d|$TYy!nlnq-RZMjqehS z4(2yrmX9W%vARpFJlNfMW#iMS;aQhGV|ka5d$heV8Jwf^P5+GfT|x$f0ge~{>G(oL z8BTMRe@Vwd623w#HJs-n96vb0o=;}BpUuFjwE_z1%xg{p|1T(Jb|T=y+JATo$nV}# z9j_Ep6P|JI4o^mowaMker^A&d#IdK{?5h?79#;)-9jf_98jyzg0oQ)pRrhqWhp%s| zZGYz3Ql^Gy9mACo@=&geeT{dMS(YW=k1qH4f655o9vKy0kRhx5ij!_O;j1qAEhAu> z-conm8Lvcc6j1HW@ynv7r<(%E^5$PI_z-X4f07^KH6>qD`|&k3d=NolYX1;H6i>OU z>KbT6W_?gkNIowt9gFS7tNVK&{TWs&?z9!S**{{ulaI75(2YS1*=UyVm#xsBHIL@n1 z+@VhCW!S9N)eYWjj9d1cCQ6+Pd$`=M)?xi{9W^I>kf|SJ|BoL(eE0sFpJ#>dgg(4U zYIbx-#25Z12j3=-?#t=*g-Q6e1YWrCf5l+A!UrTiXZF_^K=fm3Ua-=A(!t$nEF!VNPvj*2=TiaO8F&IuH=S;0Q$N%m+O zv`$_T(BwA2#Vg@4=3#BTubPro#985nb%5i>kkQWGlF8wneQ>YixPxE_EkV16dk74pokzQhI|4@0h(~+{cMKdxTSwc%9S29y zuApt=PJp9myJ#D@li*lh^A7DY?r|`MR)DsSI}J{tm7-n3Jqb>s5ut5xK-I=4m`wTvRW}J2RMVtrP$iSpLG?=+1XV4C6I82|e_Bw5QkOw>Ni_ylBQ+T4 z41vlERe*X6HGnD$b%5Fm4S_0(PO2G;&F_DR;*!j>y`QhN+@&! z3McFbluPIV6iMgy5$LNLlRJ$%Qck2U9DF@Z8Mg3qBoJMN|wdQU=(L8l@weV)3Cfy3un%jXIZIijz zYxRWsb0=>e!=jvXm%tXd3a$h7;z7`D?(WNaRzXoiRtHwgRqxEF$~&)W${{cU)ELLi zy&p87e}1Y4dK^re`?yR4>LoK+d}?mW+{0}ROaoeEZWhd#`|>*t&<=Cjjk90@%$fUU zP@^=w+&s7d7R~+gq6U`0WsuK}U+1M=1lPcZxhFqK+hjp%xrUohfu|2Ouw(Z1GgzJ! z@OwaHmU>%S8e{-5_r%v1nw8(vXb~Xro=E#*f3uD`4U_=F?JGewAeO$?Z0I>nA;`WS z5NXfp>^V|>i`jIErVw*a%snypoYOHoEns29|sda$vS0r^Rea;L9avZB+!vHV|IH$quOEuPB1k9O;!U)%~^5HE9DqfGb#>0%@R=ohw+J1zf!1LQn)^uoo1A z5}Ev--A#n@mD;OKLzZr~|rkN;Uv4 zHc1m`1}%VNF=+$sfO}4&3n|?!(KVF+O}YUwO1SnUeLxq|exPe8Cs@LXlW^)L9A1g; lqFgYNQOjmNb(Sst!~G9)kbPE{g%$@O2{JPZB_%~qMhZ=^vBm%Z