Modular computation of diffusion matrices

[lfarv]

While the tracking in AT is modular (pass methods may be modified or added easily), the computation of the diffusion matrices used in ohmi_envelope is centralised and difficult to maintain.

The computation is done in findmpoleraddiffmatrix.c for all magnets, straight or curved. Wigglers are ignored. This C function combines the tracking and the computation of diffusion matrices and is based on Bnd/StrMPoleSymplectic4RadPass.

This PR is a proof of concept where the computation of diffusion matrices is included as optional in each *RadPass passmethod. This allows the computation to be consistent for each tracking variant, and has the following main advantages:

• it allows the computation for wigglers (to be done),
• it gives a better accuracy for large off-axis orbits with pass methods different from Bnd/StrMPoleSymplectic4RadPass (on-axis, all methods should be equivalent).

As an example, in this branch, the following passmethods include the computation of the diffusion matrix: BndMPoleSymplectic4RadPass, StrMPoleSymplectic4RadPass (with strictly identical results as before) and ExactMultipoleRadPass (with a slightly different result for large off-momentum). A new C function diffusion_matrix.c is added to replace findmpoleraddiffmatrix.c using the new pass methods while keeping a similar interface. ohmi_envelope is modified to call this new function, with a test mode allowing to switch easily between findmpoleraddiffmatrix.c and diffusion_matrix.c.

This approach will require all *RadPass methods to be modified, which is a significant task. As it is now, a non-modified method generates a zero contribution to the diffusion (but still propagates correctly the cumulated matrix): as if the element does not contribute to the emittance.

This PR is open to discuss the potential of this method and identify the difficulties before going further on. It follows the discussion in #715 and integrates the proposal of @carmignani (#715 (comment)) which makes the transition easier.

[lfarv]

The processing of the diffusion matrix requires two things:

1. For each integration step, one needs to propagate the current diffusion matrix. This requires the 6x6 transfer matrix of each of the elementary steps, computed around the closed orbit. Up to now, the transfer matrix is computed analytically for the following steps:

   • “classical” drift
   • “exact” drift
   • kick in straight geometry (same for "classical" and "exact")
   • kick in curved geometry (same for "classical" and "exact")
   • “classical” pole face angle (magnetic wedge + fringe field)
   • Y rotation

   Still missing for implementing the “exact” bends:

   • Curved drift
   • “exact” magnetic wedge
   • “exact” dipole fringe field

   Building the transfer matrix analytically becomes difficult for those: taking the derivatives of E. Forest’s formulae is tedious. It’s possible to build the matrix numerically in the C integrator, but this might be too slow at run time.

2. For the “radiating” steps, one needs to compute the increment to be added to the propagated diffusion matrix. This is already done for:

   • “classical” kick in straight geometry
   • “classical” kick in curved geometry
   • “exact” kick in straight geometry

   Still missing is:

   • “exact” kick in curved geometry

So there is a significant work to provide the “exact” bends. ON the other hand I think that all ingredients are there for the wiggler.

[lfarv]

As it is now, this is already better than the present computation: one could redirect the missing (non-upgraded) passmethods to BndMPoleSymplectic4RadPass, with a result identical to what we have now. But this implies that we keep a hard-coded list of upgraded methods. Otherwise, one can easily check that a passmethod does not modify the diffusion matrix, and in this case revert to BndMPoleSymplectic4RadPass, but this doubles the processing time…

The hard-coded list avoids that, but it has to be maintained… If you think of another way, let us know.

If this is solved, one could think of merging this as it is, or with the GWigSymplecticPass, and have a smooth transition while other passmethods are upgraded. Please advise!

[lfarv]

This study revealed three problems with the present implementation:

1. The rotation around the z-axis at entrance or exit is not taken into account. So a vertical bend is not correctly taken into account. This is easy to correct,
2. The entrance and exit face angles correctly propagate the diffusion matrix, but they do not generate diffusion. For off-axis closed orbit, I think they should: it adds or remove some magnetic field, so it should radiate, it's similar to quadrupole radiation. Implementing this is not obvious,
3. findmpoleraddiffmatrix.c still uses old physical constants. For comparison with the new method, I had to modify it to use the new ones. Results are then identical to machine precision.

I propose to open another pull request to solve points 1 and 3: it will allow vertical bends and make comparisons straightforward. What's your opinion on that?