Adding support for beam precession #137
Conversation
|
A little difficult to structure a response, as there is a lot of very good stuff here. So I've just worked through it point by point, I think that will get us going in the right direction.
The first three points look good, I'll have a look at the changes for diffraction library generation once I've finished this comment
I'm not really working on this at the moment, but once we had your gridding in diffsims, it seemed worth making it clear in pyxem that the workflow you describe is an option (and probably a target)
So I think you can create those simulations with the code as it currently is, but it's clear that we want to improve the pipeline in this regard
I have thoughts about this, but I think I'll leave them to somewhere else
This one is a mixture of pyxem and orix (CrystalMap) and should be fairly simple to deal with. pyxem/pyxem#654 will be of interest.
In my experience this tends to happens when your relrods are too short.
Think this is probably a case of the more the merrier
This is a bug,
Philosophical. The
Does the solution proposed above cover these needs?
Have you seen https://journals-iucr-org.ezp.lib.cam.ac.uk/b/issues/2019/04/00/je5015/ ? I think that would save you doing the integrations (and could be written in as a "shape factor model"). |
|
Thanks for the quick thorough response, sorry for always hitting you with these long posts. They help me keep the overview myself.
Yes, exactly, which got me looking into how it was working in the background. It then seems that the default
Ok, that clarifies a lot for me, thanks. But if the sample dependent things are to be split of, would the following also not fit best in
Do you have a DOI link or so to this? I don't seem to have access. Indeed having an approximation for the precession wouldn't be a bad idea, especially for the library generation. But you might as well have the option for a full calculation if you only need a few DP's for some reason.
No worries, as I said it's still a preliminary mess, but in this way you can check the progress if you want |
I think for a while we had written in the documentation described as often 1/sample_thickness and 1e-1 corresponded to a very thin sample (10nm maybe, I would have to check)
I'm not sure. They are currently quiet underdeveloped parts of the code and I would have to investigate before coming to a conclusion.
I think adding a **kwargs options to diff_lib = lib_gen.get_diffraction_library(library_austenite,
calibration=diffraction_calibration,
reciprocal_radius=reciprocal_radius,
half_shape=(half_size, half_size),
with_direct_beam=False,
shape_factor_model="linear")to work.
https://doi.org/10.1107/S2052520619007534 Sorry, forgot that my browser extension makes a mess of copy+pasting links.
Yes, I agree. |
|
@pc494 Thanks for the link to the paper, I've implemented both the lorentzian and the closed form approximation for precession using the lorentzian. I modified the function slightly from the paper so that it also takes the parameter
I've decided that shape_factor_model is rather something related to the model we choose, rather than something dependent on the material. max_excitation_error and debye_waller_factors do depend on the sample. To me this is a logical split, and I think it is now more or less ready for my purposes. All the docstring documentation should be ok. Before I start messing around in the tests, I wanted to hear your thoughts on this way of organizing the parameters. If ok, I will update the broken tests and write new tests for the stuff I added. |
|
@din14970, at a first glance this look like a good set of decisions to me. I'll try to get a review of the code back to you later today or tomorrow if that's okay with you. |
| return sinc(excitation_error, max_excitation_error, minima_number)**2 | ||
|
|
||
|
|
||
| def atanc(excitation_error, max_excitation_error): |
There was a problem hiding this comment.
this should have a minima_number I think?
There was a problem hiding this comment.
to be consistent with the sinc functions yes, but its a bit awkward since it's a smooth function. And just like max_excitation error, the parameter just rescales in the x-direction. I added it, it can be removed later if necessary. In any case, I would expect the Lorentzian to become the default, since there is a reference for it.
|
Just to say - this looks great, thanks! I'd quite like to review it though at some point, can someone ping me when appropriate. |
|
@dnjohnstone @pc494 I would tentatively say it's ready for a more thorough review. The functionality that I wanted is there. The test coverage is again 100%, but correctness of the physics is of course best verified with additional pairs of eyes. @pc494 I didn't implement 1 of your comments above, hope you understand my reasoning. |
|
I will aim to review this tomorrow, and then ping dnjohnstone, it might take a week or so to get this all done, I hope that's okay with you @din14970. |
|
@pc494 I'm in no rush, it's your package take your time. I'm working on the next step for the indexation. |
Just to say @din14970 - I think we're getting to the point where if you're planning to continue contributing around the level you have been recently, then we can have a conversation about making "your package" more like "our package" - ping us an e-mail at pyxem.team@gmail.com and we can set that chat up... It would certainly be good for us to add another person to the core developer team for diffsims before too long if something like that might fit with your personal roadmap. |
|
I made an algebraic error, which means I thought the function below was incorrect, it's not, but for clarity I reproduce my working and thinking, the conclusion of which is that while this could be better documented, it is essentially correct. def _z_sphere_precession(phi, r_spot, wavelength, theta):
"""
Returns the z-coordinate in reciprocal space of the Ewald sphere at
distance r from the origin along the x-axis
Parameters
----------
phi : float
The angle the beam is currently precessed to in degrees.
If the beam were projected on the x-y plane, the angle
is the angle between this projection with the x-axis.
r_spot : float
The distance to the point on the x-axis where we calculate z in A^-1
wavelength : float
The electron wavelength in A^-1
theta : float
The precession angle in degrees (angle between beam and optical axis)
Returns
-------
z : float
The height of the ewald sphere at the point r in A^-1
"""
phi = np.deg2rad(phi)
r = 1 / wavelength
theta = np.deg2rad(theta)
return -np.sqrt(
r ** 2 * (1 - np.sin(theta) ** 2 * np.sin(phi) ** 2)
- (r_spot - r * np.sin(theta) * np.cos(phi)) ** 2
) + r * np.cos(theta)So that we're clear on our terms. So we arrive at perhaps my first misunderstanding, So, I'm now going to try and reproduce the final formula. Firstly, our origin is the origin of the reciprocal lattice (0,0,0). For our given geometry the center of the Ewald sphere is a function of (R sin theta cos phi , R sin theta sin phi R cos theta) = {for readability later} (A,B,C) Meanwhile our spot lies at ( X , Y, Z ) where Z will be our output. Assuming what I wrote about (r_spot, 0 , Z) We want to solve for Z, which we do by stating that the distance between these two things will be R So we get to R**2 = (r_spot - A) ** 2 + (Z-C) ** 2 + B ** 2 Compared to: -np.sqrt(
r ** 2 * (1 - np.sin(theta) ** 2 * np.sin(phi) ** 2)
- (r_spot - r * np.sin(theta) * np.cos(phi)) ** 2
) + r * np.cos(theta)So the last term in the python code is C (great), and the Excellent, and expensive set of typing for an algebra mistake on my part, but we can now be fairly sure that this function is correct. Going to go do something else for a bit and then will see what's next for review. |
| scattering_params : str | ||
| "lobato" or "xtables" | ||
| minimum_intensity : float |
There was a problem hiding this comment.
Would it be better if this were relative to the most intense spot? I'm not actually sure personally, but it's tricky to know how to set this without units... Probably it's fine as is.
There was a problem hiding this comment.
before it was just hard-coded as 1e-20 inside the function, I figured it might as well be a parameter. I guess some absolute value cut-off is most fair across different diffraction patterns in the library, resembling a kind of detection limit of a detector. On the other hand I guess different phases would require different cut-offs. Is this logic correct?
| scattering_params="lobato", | ||
| precession_angle=0, | ||
| shape_factor_model="linear", |
There was a problem hiding this comment.
Do we want to change the default here now that we have better options?
There was a problem hiding this comment.
My vote is on Lorenzian since we can point to a reference and should approximate a 2-beam rel-rod. Linear seems a bit arbitrary.
There was a problem hiding this comment.
Yes. Let's change this to Lorenzian.
|
Also, don't forget a changelog bump |
|
@din14970 @pc494 - this isn't really the right place, but thought you might both like a look at this paper if you haven't seen it yet: https://scripts.iucr.org/cgi-bin/paper?iv5009 |
Which step have you moved on to? |
While I don't mind helping out, I don't know whether I can commit to this for certain at this point. I'm now dealing with a specific problem which I would like solved within reasonable time, and pyxem/diffsims provides the perfect scaffold for this. It is mostly my frustrations with the rigidity of ASTAR that is driving me at this point. Still I'd like to hear from you where you think the package is headed. I'll send you an e-mail.
The indexation which simultaneously solves for in plane angle. I'm close to a kind of prototype which I'll share maybe in a Pyxem issue. Then we can discuss how to best integrate it there. |
|
email has arrived in the correct box and we can continue that conversation there.
Yep sounds good to me. Fingers crossed, it should be fairly easy to slot into the framework we currently have. |
What I'm working on here is very premature, opening this PR to get some opinions/feedback and to give you a preview of the code. Apologies for wall of text.
What does this PR do? Please describe and/or link to an open issue.
beam_directions_grid_to_eulerfunction.More details
In the context of my issue #124 and the other discussions we've had elsewhere, I've started looking into the entire chain of operations I would need for template matching with Pyxem. It seems @pc494 you are also kind of working on something in this direction judging by pyxem/pyxem#656?
The things I need for this are:
The simulated diffraction patterns generated with diffsims looked very odd for my material, in some cases only one diffraction spot was in the pattern. So I looked into:
shape_factor_modelsDiffractionGenerator.calculate_ed_dataWith regards to the shape factor function, I saw the default shape is linear. I looked into W&C and at least for a 2-beam case, Howie-Whelan predicts a
(sin(x)/x)^2- like function, but the exact shape depends also on the extinction distance and thickness. I messed around a bit with the shape factor functions, which you can see here:sincwas a function that was already there, and seems to follow Howie-Whelan quite well. I also triedsinc^2but this tapers off the tails too strongly. To get a smooth function that approximately follows the shape of the Howie-Whelan rel-rod I also implemented anarctan(x)/xfunction. I changed all the functions so that the maximum at s=0 hits 1. I figured it doesn't hurt to have a couple more shape factor functions.To the second point, I felt like the options available for generating a diffraction library "the official way" are a bit lacking:
Many of the options in
calculate_ed_datadon't seem to be available this way. I saw that one used to be able to setmax_excitation_errorinDiffractionGenerator.__init__but no longer, what is the reason for this?Personally, I would like to have full control over the way the patterns are simulated. So I'm inclined to change everything either in the init of
DiffractionGeneratoror in theget_diffraction_libraryfunction.In the process of looking into these diffraction patterns, I added the option to include the effect of beam precession. To calculate this I only use geometry; I don't make any additional assumptions. The Ewald sphere is precessed through the rel rods and the shape factor of each spot is calculated by averaging the shape factor function via integration. The difference is minor but there is some:
Especially further away reflections become a bit brighter with precession. The influence depends on the chosen shape factor model - above I use the
atancfunction. Below is withlinearI think the math on the precession is sound but a second pair of eyes might verify this. Of course with this implementation, precession is much more expensive, since I do an integration for each relevant diffraction spot.
Are there any known issues? Do you need help?
The main thing I'm unsure about is which arguments should go where, and why certain arguments were removed.
Work in progress?
For reviewers
later.
__init__.py.the unreleased section in
CHANGELOG.md.