Advection fields for the Monte Carlo Propagators

[ehlertdo]

Hey
Background
I am looking at scenarios where the cosmic rays propagate close to the source in environments where advective transport dominates over diffusion for $E\lesssim1~\mathrm{EeV}$. As far as I can tell, the advection fields currently implemented in CRPropa are only meant to be used with the DiffusionSDE solver. However, it appears to me that it's not possible to include interactions (photopion, photodisintegration, etc.) for that solver. Please correct me if I'm wrong on either of these two points.
My current approach
I have modified the PropagationBP code to accept and apply advection fields in a very approximate manner (simply move the particle by $t_\mathrm{step} \cdot v_\mathrm{adv}$ each step). See PropagationBP.txt ( I can also create a pull request if necessary). It's clear that this solution is not exact as the final propagation distance of particles (column 'D' in the output files) can become shorter than the system size. I assume that this occurs since the extra movement due to the advection is in addition to the regular step size. This also leads me to believe that the interactions are likely underestimated since the propagation distance assumed by CRPropa is smaller than the true path length. In addition, simply adding the normal/diffusive motion and advection does not respect $v_\mathrm{CR}\leq c$ by default.
My question / request
Is there a way to include the advection into either the Boris-Push or Cash-Carp solvers directly instead of adding it afterward as a separate step? (Alternatively my problem can also be solved if there is a way to include interactions with the DiffusionSDE module).

Best,
Domenik

[ehlertdo]

Short update: I have now included the advection in a way that is more realistic by taking it into account in terms of the "effective" propagation direction of the cosmic ray instead of simply adding the advection on top of the usual propagation as an additional step. See also the attached slide for more explanation.
If anyone has some comments on this that would be appreciated.

PropagationBP::Y PropagationBP::dY(Vector3d pos, Vector3d dir, double step,
			double z, double q, double m) const {
		// velocity of advection field [m/s]
		Vector3d vWind = getAdvFieldAtPosition(pos);
		// add advection vector [(vx vy vz) / c] to propagation direction
		Vector3d dir_tot = dir + vWind.getUnitVector() * (vWind.getR() / c_light);
		// renormalise to unit vector
		dir_tot = dir_tot.getUnitVector();

		// BORIS-PUSH ALGORITHM ========================
		// half leap frog step in the position
		pos += dir_tot * step / 2.;

		// get B field at particle position
		Vector3d B = getFieldAtPosition(pos, z);

		// Boris help vectors
		Vector3d t = B * q / 2 / m * step / c_light;
		Vector3d s = t * 2 / (1 + t.dot(t));
		Vector3d v_help;

		// Boris push
		v_help = dir + dir.cross(t);
		dir = dir + v_help.cross(s);

		// include advection for the second half step
		dir_tot = dir + vWind.getUnitVector() * (vWind.getR() / c_light);
		dir_tot = dir_tot.getUnitVector();

		// the other half leap frog step in the position
		pos += dir_tot * step / 2.;

		return Y(pos, dir);
	}```