As feared, this has turned out to be quite a bit of work. Up to now, the code was matching my model and the data by essentially doing the brute force approach, in which it tested all the possible delta_x and delta_y parameters in order to determine the best fit. 

It turns out that is not possible to extend this to the case where we are fitting two donuts at one - it simply takes too much time, because the problem scales as N**2.

I have also discovered that my previous interpolation which I used to interpolate at scales lower than my oversampling scales (x10 at the moment) was not entirely correct. I semi-knew that before as I implemented the code very crudely initially. 

I have been experimenting with the approaches in which I use a small optimizing loop in order to find the best centering (instead of previous brute force approach). Given the wide variety of possible combinations (defocus to focus, different values of separations between spots) it is not obvious to me that I can find a fast clean algorithm that would always converge to correct solution under different circumstances.

The plan is now to measure the effective ratio and separation between close spots that I would like to study. The fitting algorithm will then accept those values to put tight priors on those centering parameters during fitting, which should help with speed and convergence.  On some level, this information is contained in https://github.com/Subaru-PFS/obs_pfs/blob/master/pfs/lineLists/ArCdHgKrNeXe.txt 

I have implemented the possibility to fit 2 donuts at the same time in the algorithm (e.g., see Figure ``2imagesExample''). 

In the process, while changing on how the interpolation and centering is conducted, I have significantly slowed down the algorithm (even in the case of single psf available in the image). The slowdown is of the order of 50% to 100%. I will investigate and hopefully speed it up in a separate ticket.