Improving/bugfixing how we simulate precession

[pc494]

I think our current precession simulations could be improved a little. Currently we use "max_excitation_error" to eliminate spots from our precessed simulations, but this isn't quiet accurate as it doesn't account for the size of the g vector sensibly.

Case 1 (no approximation):
This finds the average excitation error, but a spot can have an average excitation error greater than max_excitation_error and we would still expect to see it in a pattern as long as somewhere on the circuit it had a sufficiently small excitation error. This is a bug (I think).

Case 2 (approximation):
These spots are eliminated before we even consider that there is precession going on.

Proposed solution:

Working from Equation 10 of [1] which provides the intensity of a spot as a function of:

g-vector (we have this)
precession angle (we have this)
the unprecessed excitation error, which is exactly what we have
the extinction distance (we can use max_excitation_error as a proxy)

(and is currently used in "lorenztian precession)

What Figure 3 (same manuscript) shows is that changing max_excitation_error doesn't change the width of the camel plot significantly. This is likely because we are in the regime (below equation 11) where we have full integration of our spots. We could create two cases (basically by hand) and use them as the eliminations. They would be:

~ numerical_factor * g * phi (precesison regime)
~ The 98th percentile (or similar) of the phi = 0 lorenztian

I would take this opportunity to remove the "manual integration" as I'm not sure this is a particularly good option.

[1] - Specifics of the data processing of precession electron diffraction tomography data and their implementation in the program PETS2.0 - Lukáš Palatinus et al.

Code

As it currently is.

# Identify points intersecting the Ewald sphere within maximum
        # excitation error and store the magnitude of their excitation error.
        r_sphere = 1 / wavelength
        r_spot = np.sqrt(np.sum(np.square(cartesian_coordinates[:, :2]), axis=1))
        z_spot = cartesian_coordinates[:, 2]

        if self.precession_angle > 0 and not self.approximate_precession:
            # We find the average excitation error - this step can be
            # quite expensive
            excitation_error = _average_excitation_error_precession(
                z_spot,
                r_spot,
                wavelength,
                self.precession_angle,
            )
        else:
            z_sphere = -np.sqrt(r_sphere ** 2 - r_spot ** 2) + r_sphere
            excitation_error = z_sphere - z_spot
        # Mask parameters corresponding to excited reflections.
        intersection = np.abs(excitation_error) < max_excitation_error
        intersection_coordinates = cartesian_coordinates[intersection]
        excitation_error = excitation_error[intersection]

[pc494]

@din14970 this is your chance for objections, as I am planning to implement this over the weekend

[din14970]

@pc494 thanks for the ping. Maybe I'm a bit dense, not sure I fully understand the issue, some plots or specific lines of code may help. If you find there is a bug feel free to modify or remove. The integration approach is useless from the point of view of creating a template library, but maybe someone cares to check the difference between this "full calculation" and the approximation?

As far as max_excitation_error is concerned, some time ago I thought about the fact that we use this both for defining the width/shape of the rel-rod and also as a kind of cut-off. I thought then it might be more appropriate to split this further into two parameters but didn't look at it since. So definitely there is still definitely a bug in there, but I forgot to report it. Don't know if we are talking about the same things here.

[pc494]

Okay, this is actually two bugs as far as I'm concerned, but the most obvious illustration of the one I'm going to fix is the following:

When we precess we should see more spots, here we're slightly titled away from zone axis but both precessed and unprecessed have exactly the same spots.

[HelvoortTon]

Hi, I find s_max a difficult concept. Have not a direct input on raised issues, sound logic. However, Endre Jacobsen has in his master tested the effect of chosen s-max on phase mapping of ppts in 6xxx alloys (contact me or Endre for a copy). Seems there is in the best case an overlap for s-max ranges for different phases (worse case: no overlap, choice of s-max determines what phases "found"). Could test s_max to find suitable ranges be interesting (if it could be done practically, Endre had a reduced library) beside a calculation/estimation of s_max?
(For accuracy orientation mapping the approach as in PETS might likely be the best).

[pc494]

I think for unprecessed data s_max is likely to be a real issue, however I have found that for sufficient precession angles some much of this problem goes away (the associated rocking curves are not very dependent on s_max). This is likely to be written up for my own thesis.

In general, following PETS seems to be the right idea.

[din14970]

Sorry for the delay, I looked into the code for dealing with #163 and re-read this thread. Indeed there was a bug where the spots were filtered based on unprecessed excitation error, so I think you @pc494 fixed this in #161 . Note however that the else branch in calculate_ed_data will likely be more computationally expensive now since you calculate this on the entire "sphere" of points instead of just a small slice. Perhaps it is still worth thinking about whether a fast pre-filtering is possible. If we assume the ewald sphere is flat, then we could separate out those points that intersect with and are between the two cones swept out by the ewald surface. That means if abs(z_reflection) - max_excitation_error < r_spot * tan(precession_angle). If the ewald sphere is very curved for very small acceleration voltages then it may break down.

To resolve the second issue I brought up I figured it is easiest to introduce a parameter shape_factor_width which determines how wide the shape functions should be, whereas max_excitation_error only determines the cut-off. If no shape_factor_width is provided then it is set by default to max_excitation_error giving the same behavior as before, but now with the flexibility that you can independently control shape factor and cut-off. I delved through the entire simulation chain and couldn't find other glaring bugs. In the process I also implemented simulation whereby the atomic scattering factor is ignored (you can pass scattering_param = None), strongly mimicking ASTAR behavior:

For reference, the simulations by ASTAR for the same crystal and similar conditions:

I will start a PR with these notes and update the tests later.

[pc494]

this is considered fixed (#161 and #166), closing