Update spin/jet AGN model

This implements lots of small updates that have accumulated in the last few months. I will try to group these into separate categories for clarity. I know that this is quite a lot, sorry. I probably shouldn't have waited this long with all the changes, I didn't quite realize how many there were until I started making this MR.

Improvements to kicking scheme for jets:

1. There was an unnecessary loop for sorting particles in black_holes_iact.h (old lines 370-376). This loop effectively meant that the same particle occupied the entire ray. The end result was that whenever the code wanted to kick more than one pair of particles in a time-step, it would just give the extra energy to the same two particles. The loop has been removed.
2. One-sided kicks (where only one kick happens, from one side of the BH) are now prevented from happening. This is accomplished by lines 1224-1249 in the new version of black_holes.h, where the condition (bp->rays_jet[i].id_min_length != -1) && (bp->rays_jet_pos[i].id_min_length != -1) makes sure that both hemispheres have to have two particles to kick.
3. Kicks are now perfectly antisymmetrical, as suggested by Joop recently. They are still random, but whenever a pair is kicked, their added momentum is along a line separated by 180 degrees. This, however, requires that the direction for the kicks is randomly drawn at some point earlier than in runner_iact_nonsym_bh_gas_feedback() within black_holes_iact.h. The way I implemented it now is that I define a new variable (bp->jet_direction[3]) for each BH that is randomly drawn every time the code decides a BH will hand out some number of kicks (this is done in black_holes.h, new lines 1269-1274). Since this is now done at the level of the BH and not individual particles, I had to also change the function random_direction_in_cone() to accept one fewer argument, since the particle ID (of the gas particle) is now unnecessary. This is why random.h and tests/testRandomCone.c change.
4. The ray_jet_correction term, which is added to quantities when sorting particles in the jet rays to make sure that already kicked particles end up at the latter end of the rays, are now included less liberally. Previously, this term was added (it was non-zero) when the relative velocity between the BH and the gas particle was a relatively small fraction of the target jet velocity, i.e. we used it only if relative_velocity > 0.3 * bi->v_jet. However, in realistic cases, bi->v_jet can be only a few thousand km/s, meaning that particles with velocities of order 500 km/s relative to the BH were being excluded from being considered for kicks. Now we use relative_velocity > 0.8 * bi->v_jet, which should still make sure that already kicked particles don't get kicked again, but it should not be as restrictive.

Changes to the super-Eddington regime:

1. We now assume that the slim disc modelling of super-Eddington accretion is active once the disc becomes super-Eddington, but using a spin-dependant Eddington ratio, rather than the usual one (that uses 0.1 as the radiative efficiency in the definition of the Eddington accretion rate): new lines 515-520 in black_holes_spin.h. Previously this was being done in a more complicated way by using the parameter TD_SD_eps_r_threshold, which is now removed.
2. The aspect ratio H/R of the slim disc is now assumed to be the same as that of the thick disc at very low accretion rates (new lines 561-563 in black_holes_spin.h). This means that jet efficiencies and spindown are the same in the two regimes. Previously there was a small difference in the two aspect ratios, but for simplicity and because of theoretical uncertainties they are now the same.
3. The slim disc modelling is now on by default in the model; include_slim_disk is set to 1 in the example parameter file and example parameter block in doc.

Changes to the variable velocity schemes:

1. A maximum jet velocity has been added for the variable velocity cases. This is currently set to 1e5 km/s, i.e. c/3. This has been added to the parameter example file and code block in doc.
2. The 'Local' and 'SoundSpeed' velocity schemes have been merged into a single one, called just 'Local'. This implements both of the conditions that we discussed when I was in Leiden: the one based on the replenishment time-scale of the BH kernel, and the one based on the crossing time-scale. Both are implemented as minima to the jet velocity. The replenishment part of this scheme uses a hot gas sound speed, which is calculated using only particles in the BH kernel whose temperature is above some minimum temperature_hot_gas_min_K (1e4 at the moment, matching the ISM) that is set in the parameter files. The velocity dispersion of gas is also included in this calculation of the replenishment time-scale.
3. Added the last jet velocity last_jet_kick_velocity as one of the tracer variables, for these variable schemes.

Additional changes to parameters:

1. The mdot_crit_ADAF parameter, which is the critical Eddington ratio at which modelling of the accretion disc changes from thick to thin (and feedback from jets to thermal) has now been correctly set in the parameter files (it was missing previously, but it was in the code).
2. The fiducial velocity scheme is now 'Constant', and the value is now set to 3160 km/s (the fiducial value for low-res COLIBRE at the moment).
3. The subgrid accretion disc viscosity parameter alpha_acc has been set to 0.2 instead of 0.1.

With regards to tracers, it seems that the jet-related ones weren't properly included at all for the EAGLE scheme, so I have now done that. This is why the tracer files change (other than the new tracer last_jet_kick_velocity).

Finally, I added an accretion efficiency parameter. This is applied only if the BH is in the thick disc regime, i.e. if the Eddington ratio is below the already-mentioned mdot_crit_ADAF. I have removed all the 'ADIOS'-related parameters; these were just a more complicated (and realistic, but we don't need it) way of accounting for an accretion efficiency.

Edit: I added some further changes:

1. In the previous version, before spin was updated in black_holes.h, some other quantities (like the jet/radiative efficiencies and the increase to the jet/thermal reservoirs) were already computed before it was updated. The code now updates spin first, and only later computes the other quantities.
2. I added another updating of the jet and radiative efficiencies right after the line decide_mode(bp, props) near the end of black_holes_prepare_feedback in black_holes.h. This is a very important change. The decide_mode(bp, props) function can make the BH change feedback modes, but that is not reflected properly until the jet and radiative efficiency fields of the BHs are updated, so it must be done straight after decide_mode(bp, props). If this is not done (like in the previous version), the BH may have, for example, a high jet efficiency and a zero radiative efficiency, when it should have a zero jet efficiency and a non-zero radiative one. The effect is that the BH timesteps were being computed incorrectly, since black_holes_compute_timestep uses the fields bp->jet_efficiency and bp->radiative_efficiency to estimate when the next feedback event will be.
3. I have added an additional spin-related time-step. This one makes sure that the spin vector (and therefore jet direction) don't flip by a huge amount over a single time step.

I have tested the above three changes and things seem stable. It is possible that the last one slows the simulations down somewhat (it doesn't seem to at low-res and in small boxes, but that may change).

Edit 2: I have modified the calculation of the orbital angular momentum of two merging black holes. Previously this was being calculated in the frame of one of the BHs, but I think that's incorrect. Now, the centre of mass is calculated and the coordinates and velocities of each BHs calculated to be in the centre-of-mass frame. The total orbital angular momentum is then simply the sum of each of those for the 2 BHs.

Also, when calculating the total angular momentum (including the orbital contribution and the two spins), I wasn't correctly calculating the angular momentum of each BH (the one coming from spin), which has now been corrected.

Edit 3: The parameter fix_jet_efficiency used to not only do what it says, but it also fixed the jet directions to be along the z-axis. This has now been changed, a new parameter, fix_jet_direction, is used to do that.

Edit 4: I updated the idealized cluster example to do the following: 1) Have the necessary parameters in the parameter files to also run with jets, 2) modified the README to let the user know how to configure the subgrid model, either EAGLE or SPIN_JET_EAGLE.