After spending a little more time with the data, I'm beginnig to think the scan was too fast to be able to do a reasonable particle-size calculation. The error bars on the data may be too large to do a meaningful analysis of the various contributions to the peak broadening (see post#4, above).
In any case, here's how I would do it from the given data:
1. For each of the 5 peaks, do a Gaussian fit to the data (limit it to about 3-sigma from the mean) and determine the width(\beta) and the peak value (\theta) from the Gaussian (not from the raw data), as well as the errors in these values \delta\beta~,~~ \delta \theta. For the baseline subtraction in half-maximum determination, I would use two sets of baseline data about 3 degrees on either side of each peak and use a mean of both sets (this is to minimize any angular dependence in the background).
2. Calculate the 5 values of \beta ^2 cos ^2\theta and plot these points (along with their error bars) against sin^2 \theta, for each of the two data sets (as received, and 10-hrs milled).
3. If you are lucky, and the error bars are small enough, you will be able to fit a straight line to both sets and extract a slope and intercept for each one.
4. What these 4 numbers represent:
(from post #4)
\beta^2 = \left( \frac{A}{cos \theta} \right) ^2 + (Btan \theta)^2 + C^2
\implies \beta^2 cos^2 \theta = A^2 + (Bsin \theta)^2 + (Ccos \theta )^2
= A^2 + (Bsin \theta)^2 + C^2(1-sin^2 \theta )
= A^2 + C^2 + sin^2 \theta(B^2 - C^2)
From the linear fit, you have slope=B^2 - C^2~, ~~intercept=A^2 + C^2.
(i) For the as received powder, B (the stress-related broadening coefficient) should be small and hence it might be reasonable to neglect it. Note that, this would mean the slope for this line should be negative. If you do not have a negative slope, you can not neglect the stress in the as received sample; this makes the calculation more tricky.
(ii) The broadening coefficient from instrumental factors should be sample independent, so C is the same number for both samples. If the intercepts for both samples are close to each other, the broadening is probably dominated by intrument-related factors.
That leaves you with 4 unknowns : A(as rec'd), A'(milled), B'(milled) and C(instr). From the two data sets you have 4 equations for these 4 unknowns (2 slopes and 2 intercepts).
5. Finally, using A=0.9 \lambda / d, you can determine the approximate particle size, d.