I first calculate a perfect PSF for the telescope in question (by
default using the analytic expression for an obscured circular pupil,
though there's a switch in my code to make it use the FFT of a
pixellated image of the Keck pupil to generate Keck PSFs). I register
the aberrated PSF and the perfect PSF to subpixel resolution by
cross-correlation.  I then shift the aberrated PSF to center the star
the same as for my perfect PSF. I perform photometry on both images
using a 2" aperture and the Goddard aper.pro procedure, and then scale
the perfect PSF to have the same flux inside that aperture as the
aberrated PSF. Finally, I compute the Strehl ratio as the ratio of peak
pixel intensities in both images.

In retrospect, I might get better performance by, after registering,
calculating a perfect PSF at the same location as the aberrated PSF,
rather than numerically shifting the aberrated PSF. I wrote the code the
way I did for performance reasons, as it allowed me to use one perfect
PSF to compute Strehl for many aberrated PSFs more rapidly. I'll try
doing it the other way and seeing if the results are any different.
I also might improve things by generating the perfect PSF oversampled
and then binning down to the proper pixel scale.

"(Results) with original version of my fitstrehl.pro code, which just
  evaluated the analytical expression for the Airy function..."

and his new version of the code...

"...which evaluates the same (analytic) expression on a lambda/8D grid
  and then rebins down to the lambda/2D or lambda/4D grid."