You won't be able to get the best results from just theoretically analyzing one or more effects that cause blurring. The best way to proceed is by measuring the blurring and then do deconvolution based on the measured point spread function. You can print a some black and white pattern costing of lines in different directions, or circles on a piece of paper. A picture taken from relatively close up, say a meter away but such that you don't get much edge deformation can be used as a high resolution reference image. Then, if you take another picture from farther way (using exactly the same focal length) you'll have a lower resolution image. The mapping from the high resolution image to the low resolution image involves all the blurring effects that you want to correct for.
Now, as mentioned in the comments in Michael's answer, you must work with the original raw images without doing any demosaicing. For this you can use the dcraw program and use the -d or -D option to produce TIFF files. These Tiff files must be split into 3 parts by picking the pixels corresponding to R, G and B from the original image. You then need to make sure that the output you get is indeed linear as it should be, otherwise you must correct for that. So, you have to get this right before you can start with the measurement of the point spread function.
What you can do is take pictures of some objects in a room with two lights A and B using a tripod. You take pictures with only light A, only light B and both A and B. Then if the raw pictures are really exactly linear, you should be able to calculate exactly the A+B picture from the A and B picture. But you must take many pictures, align them and average over to make them almost totally noise free. You must also correct for the different exposure times, after stacking you end up with some unknown effective exposure time. But that doesn't matter, because A + B is supposed to be some linear combination of A and B, and if it isn't you can easily detect that.
Then assuming that you can produce non-demosaiced images that are almost exactly linear and noise free, you can produce such images for the patterns and then use these to calculate the point spread function for various camera settings. This should be done for different F-numbers and also for deliberate departures from perfect focus. The calculation of the point spread function in each case can be done using Fourier transform methods.
With the data of all the point spread functions, you can sharpen very accurately using deconvolution. If you take pictures of some object then you should to do that using image stacking methods to remove the noise as best as possible. Also, you should do this for slightly different focus settings within the margin what you would judge as perfectly sharp. Then use the non-demosaiced linear raw files (using whatever additional transform that you have found were necessary to get to an almost exactly linear output).
Then by studying the images you can extract the point spread functions to some approximation (e.g. by looking at high contrast edges). You then compare these point spread functions with the exact experimentally determined ones, to select which one in the best. It's best to have the experimental point spread functions processed in a mathematical form where the defocussing is described by a continuous parameter. As a function of this parameter, the point spread function will to a good approximation just scale.
You then select that image that corresponds to the narrowest point spread function and then deconvolve that image. The narrower the point spread function the less noise you get (the wider the point spread function, the narrower it is in k-space, which causes the high Fourier components to be amplified, which causes the noise to be amplified).