Adding most bound BH ID and position to output

Added the position and ID of the BH that is closest to the most bound dark matter particle.
This was tested in COLIBRE and works well.

src/allvars.cxx

@@ -1811,6 +1811,9 @@ void PropDataHeader::declare_all_datasets(const Options &opt)
 #ifdef BHON
     declare_dataset<unsigned long>("n_bh");
     declare_dataset("Mass_bh", MASS);
+    //save ID most bound BH and its position
+    declare_dataset<long>("ID_mbp_bh");
+    declare_XYZ_datasets("cmbp_bh", LENGTH);
 #endif // BHON
 
 #ifdef HIGHRES

src/allvars.h

@@ -1621,6 +1621,10 @@ struct PropData
     ///mean accretion rate, metallicty
     Double_t acc_bh, acc_bh_mostmassive;
 
+    ///ID of most bound BH and its position
+    Int_t ibound_bh = 0;
+    Coordinate gposmbp_bh;
+
     ///\name blackhole aperture/radial profiles
     //@{
     vector<int> aperture_npart_bh;
@@ -2478,6 +2482,7 @@ struct PropData
 #endif
 #ifdef BHON
         M_bh*=opt.h;
+        gposmbp_bh*=opt.h/opt.a;
 #endif
 #ifdef HIGHRES
         M_interloper*=opt.h;

src/io.cxx

@@ -2177,6 +2177,15 @@ void WriteProperties(Options &opt, const Int_t ngroups, PropData *pdata){
 
         for (Int_t i=0;i<ngroups;i++) ((Double_t*)data)[i]=pdata[i+1].M_bh;
         Fhdf.write_dataset(opt, head.headerdatainfo[itemp],ng,data,head.hdfpredtypeinfo[itemp]);itemp++;
+
+        for (Int_t i=0;i<ngroups;i++) ((long long*)data)[i]=pdata[i+1].ibound_bh;
+        Fhdf.write_dataset(opt, head.headerdatainfo[itemp],ng,data,head.hdfpredtypeinfo[itemp]);itemp++;
+
+        for (int k=0;k<3;k++){
+            for (Int_t i=0;i<ngroups;i++) ((Double_t*)data)[i]=pdata[i+1].gposmbp_bh[k];
+            Fhdf.write_dataset(opt, head.headerdatainfo[itemp],ng,data,head.hdfpredtypeinfo[itemp]);itemp++;
+        }
+
 #endif
 #ifdef HIGHRES
         for (Int_t i=0;i<ngroups;i++) ((unsigned long*)data)[i]=pdata[i+1].n_interloper;

src/substructureproperties.cxx

@@ -1136,6 +1136,14 @@ private(EncMassSF,EncMassNSF,Krot_sf,Krot_nsf,Ekin_sf,Ekin_nsf)
                 mval = opt.MassValue;
 #endif
                 pdata[i].M_bh+=mval;
+
+		//save particle information if the most bound id has not been assigned.
+		//note that this works because particles has been previously ordered by distance to defined centre.
+		if(pdata[i].ibound_bh == 0){
+			pdata[i].ibound_bh = Pval->GetPID();
+			for (k=0;k<3;k++) pdata[i].gposmbp_bh[k] = Pval->GetPosition(k) + pdata[i].gposmbp[k];
+		}
+
             }
         }
 #endif
@@ -2073,6 +2081,14 @@ private(j,Pval,x,y,z,vx,vy,vz,jval,jzval,zdist,Rdist)
                 mval = opt.MassValue;
 #endif
                 pdata[i].M_bh+=mval;
+
+                //save particle information if the most bound id has not been assigned.
+                //note that this works because particles has been previously ordered by distance to defined centre.
+                if(pdata[i].ibound_bh == 0){
+                        pdata[i].ibound_bh = Pval->GetPID();
+                        for (k=0;k<3;k++) pdata[i].gposmbp_bh[k] = Pval->GetPosition(k) + pdata[i].gposmbp[k];
+                }
+
             }
         }
 #endif