Currently a function is in development that solves for FRD of an arbitrarily centered PSF generated from Neven's code, using the scipy shgo optimizer (https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.shgo.html). The current version of this code solves for FRD with plus/minus 0.2 mrad error with about 30 iterations. I am currently investigating how to reduce iteration time using the insight from my previous FRD analysis

The current version of my code is available at https://github.com/bbelland/FRDtesting. This code can run and solve for an FRD, but while comparing this code to my previous work's precision, I found that the uncertainties are larger than expected from my previous FRD analysis. I am working to address this issue to demonstrate FRDs can be solved for properly.

My initial thoughts after testing code today:

1. (minor issue) in FRDClass_TestNotebook  variables are called ``knownFRDlist'', ``imagelist'', even thought the code demands them te be arrays and not lists.

2. The uncertainty that is reported from the algorithm depends on the points that have been provided to the algorithm. I think it would be better to report where the interpolation crosses the limit. I.e., the reported uncertainty should not directly depend on the provided points, but one should estimate at which point the critical chi**2 values is reached.

3. ``uncertaintyrange'' variable returned by the algorithm is a single number, even though the distribution is nonsymmetric and reported uncertainty has different lower and upper bounds. I would assume that the ``uncertaintyrange'' variable would therefore also be consisting of two numbers.

4.  (further thoughts, not critical to this code) - at this point the code return basically what image that is provided has the best fit - i.e., the result depends critically on the density of sampling. It might be beneficial to try to extrapolate the trends as one is approaching to the correct solution from the left and the right side to extrapolate what the correct answer might be, or at least provide the user with a hint for the second round of sampling, i.e., in which range to simulate more images. 

I'm addressing the first three points, and they're almost all done.

1. Variable names have been renamed to reflect that they are arrays. This should be done now.
2. I've added an interpolation function as suggested. I need to test this to make sure it works generally though.
3. Uncertaintyrange has been changed from the difference between the lower and upper FRD ranges to a tuple with the (lowerFRD, upperFRD) at the chi squared limit, to reflect the fact that the there tends to be a larger uncertainty toward larger FRD values than lower FRD values as you indicated.

I'll be thinking about how to address point 4 but focusing on PIPE2D-623 as discussed in the telecom.

Some further thoughts:

1. I believe that the reason why you code creating images is quite slow is due to the fact you are creating 50x50 images when creating focused images. As we discussed few times, one can use smaller images which will be faster and also more realistic as in real life we will be using the code on relatively small post-stamp images (probably 10x10 or even less).
2. As a first idea to understand if the sampling has been sufficient, one can measure the gradient of the points around the deduced best solution to measure if the result would likely improve by denser sampling. If the nearby sampled points have similar quality the improved sampling would likely not improve solution much, but improved sampling would be beneficial if this gradient analysis shows that the algorithm has not reached stable solution.