By careful analysis of the halo around stars, Jupiter and the Moon, we have arrived at an understanding that the PSF is related to a power-function 1/r^alfa, where alfa is some number below 3.
It is also possible to estimnate the PSF using deconvolution techniques. In particular, the 'multiframe blind deconvolution' method produces estimates of both the deconvolved image and the PSF, without any other input than the requirement that the output image and the PSF be positive everywhere.
There is matlab software provided here: http://pixel.kyb.tuebingen.mpg.de/obd/
and a paper here: http://www.aanda.org/articles/aa/pdf/2011/07/aa13955-09.pdf
Fredrik Boberg and I - well, Fredrik - has tested the software on images pulled from a 100-stack on (nonaligned, but almost aligned) images. The output comes in the form of estimates of the PSF and estimates of the deconvolved image. The deconvolved image has an unnaturally sharp limb - but more on that somewhere else - for now, let us see the PSF!
The estimated PSF is highly non-rotationally symmetric, but does have a central peak. The 'skirt' wavers up and down over many orders of magnitude. A radial plot of the PSF and a surface plut looks like this:
The red line shows the 1/r³ profile. The elliptical shape of the PSF causes the broad tunnel seen above - the divergence at radii from 30 pixels and out is the 'wavering skirt'.
The PSF does appear to generally be a 1/r³ power law, which is encouraging. It appears to be quite flat out to 2 or 3 pixels - which also generally matches what we have seen in the PSFs empirically generated from imags of stars.
Test on a synthetic image convolved by a known PSF.
pre-align the images in the stack - perhaps a slow drift causes the non-round psf. Repeat the above for several stack from different nights.
Repeat the above for stacks taken through different filters - any indication that PSFs are different?