Fixing bug in redshift at infall calculation

Fixing a bug in the recalculation of satellite subhalo properties which
was assuming infall_t was a snapshot rather than a redshift.

I also put back in the code the definition of subhalo properties in the
merger_tree_reader because the redefinition is optional.

include/subhalo.h

@@ -234,7 +234,7 @@ class Subhalo : public Identifiable<subhalo_id_t>, public Spatial<float> {
 	float lambda = 0;
 	/// redshift at which the subhalo became a type > 0.
 	float infall_t = 0;
-	/// redshift at which the subhalo became a satellite of its z=0 host.
+	/// snapshot at which the subhalo became a satellite of its z=0 host.
 	float infall_t_z0host = 0;
 	/// halo mass and stellar mass of central galaxy at infall_t.
 	float Mvir_infall = 0;

src/galaxy_writer.cpp

@@ -364,7 +364,7 @@ void HDF5GalaxyWriter::write_galaxies(hdf5::Writer &file, int snapshot, const st
 
 			descendant_id.push_back(subhalo->descendant_id);
 			infall_time_subhalo.push_back(subhalo->infall_t);
-			infall_time_subhalo_z0host.push_back(subhalo->infall_t_z0host);
+			infall_time_subhalo_z0host.push_back(sim_params.redshifts[subhalo->infall_t_z0host]);
 
 			int m = 0;
 			if(subhalo->main_progenitor){

src/merger_tree_reader.cpp

@@ -236,9 +236,23 @@ const std::vector<SubhaloPtr> SURFSReader::read_subhalos(unsigned int batch)
 		//Assign circular velocity
 		subhalo->Vcirc = Vcirc[i];
 
+		auto z = simulation_params.redshifts[subhalo->snapshot];
+		subhalo->concentration = dark_matter_halos->nfw_concentration(subhalo->Mvir, z);
+
+		if (subhalo->concentration < 1) {
+			throw invalid_argument("concentration is <1, cannot continue. Please check input catalogue");
+		}
+
+		double npart = Mvir[i]/simulation_params.particle_mass;
+
+		subhalo->lambda = dark_matter_halos->halo_lambda(*subhalo, Mvir[i], z, npart);
+
+		// Calculate virial velocity from the virial mass and redshift.
+		subhalo->Vvir = dark_matter_halos->halo_virial_velocity(subhalo->Mvir, z);
+
 		if (!transients_prefix.empty()){
-		        // create flag to indicate this subhalo is transient
-		        subhalo->transient = std::find(std::begin(transientsIndex), std::end(transientsIndex), nodeIndex[i]) != std::end(transientsIndex);
+			// create flag to indicate this subhalo is transient
+			subhalo->transient = std::find(std::begin(transientsIndex), std::end(transientsIndex), nodeIndex[i]) != std::end(transientsIndex);
 		}
 
 		// Done, save it now

src/tree_builder.cpp

@@ -569,10 +569,16 @@ void TreeBuilder::define_properties_satellite_subhalos(const std::vector<MergerT
 				for(auto &halo: tree->halos_at(snapshot)){
 					for(auto &subhalo: halo->all_subhalos()){
 						if(subhalo->subhalo_type == Subhalo::SATELLITE){
-							// in the case of satellite subhalos, we use the virial mass of its host at the time of infall, and the time of infall
+
+							double mvir = subhalo->Mvir;
+							double z = sim_params.redshifts[subhalo->snapshot];
+
+							// in the case of satellite subhalos with a well defined infall mass, we use the virial mass of its host at the time of infall, and the time of infall
 							// to define other properties.
-							double mvir = subhalo->Mvir_infall;
-							double z = sim_params.redshifts[subhalo->infall_t];
+							if(subhalo->Mvir_infall > 0){
+								mvir = subhalo->Mvir_infall;
+								z = subhalo->infall_t;
+							}
 
 							double npart = mvir/sim_params.particle_mass;