Armstrong absorption correction

[Brian Joy]

Hi John,

Looking for some guidance...

I’ve written a matrix correction program (in Fortran 90/95) so that I can re-process output from the JEOL PC-EPMA software, as I don’t necessarily like the options offered by JEOL.  Where possible, I’ve been testing output from my program against output from CalcZAF.

For the PAP, XPP, and FRAME-ZAF models, I can duplicate f(chi) calculated by CalcZAF precisely for X-rays in both pure elements and mixtures.  However, I’m having problems doing the same using the Armstrong absorption model.  I’ve attempted to locate the source(s) of the discrepancy by combing through the code that you’ve posted on-line, primarily in the file ZAF.bas.  However, I have no substantive knowledge of Visual Basic syntax or program structure, and so I’ve had some difficulty following equations through the code; for instance, I can’t seem to find exactly where the integrals of phi(rho*z) and phi(rho*z)*(-chi)*rho*z are evaluated.  Could you point me in the right direction?

When I calculate f(chi) for K lines of pure elements using the Armstrong absorption model, I get exactly the same output as generated by CalcZAF.  However, for K lines of elements in mixtures, my f(chi) results differ from the CalcZAF output (example below).  In addition, when I calculate f(chi) for M lines of pure elements, my output always differs at least slightly from the CalcZAF output.  For instance, for W Ma in pure W, I get f(chi) = 0.6798, while CalcZAF gives 0.6727 (accelerating potential = 15 kV, takeoff angle = 40 degrees, MAC database = MAC30).  For Ce Ma in pure Ce, I get f(chi) = 0.5158, while CalcZAF gives 0.4969.  For L lines of pure elements, I get f(chi) that differs by no more than 0.0001 from CalcZAF output.

For the case of the K lines in mixtures, the discrepancy could not be due, for instance, to use of a different “mixing rule” for the mean ionization energy.  This would impact the calculation of stopping power as well, but the value that I calculate for this latter quantity in the mixture is identical to the value calculated by CalcZAF.

The expressions below (which I’ve proofread many times) are consistent with those presented in Armstrong’s paper in Electron Probe Quantitation, though his presentation isn’t really complete (for instance, he doesn’t specify the Berger-Seltzer expression for mean ionization energy).  I’ve replaced the ambiguous E in the expression for alpha with E0, and this appears to be consistent with the equation given in line 1130 in ZAF.bas.

Regarding my calculation of f(chi), I’ve checked results from the analytical expressions that I give below against results obtained using Simpson’s rule (via subroutine qsimp from Numerical Recipes).





Here I give an example of the differences in my calculated f(chi) versus those from CalcZAF for the case of Fe2SiO4 analyzed at 15 kV potential and 40-degree takeoff angle.  Some very minor discrepancies are present in the MACs (which I’ve calculated using Heinrich’s MAC30 equations), but these differences are too small to have any impact on calculated f(chi) at the illustrated four digits of precision.  For Si Ka in the mixture, I calculate f(chi) = 0.6303, while CalcZAF gives 0.6342.  For O Ka in the mixture, I calculate f(chi) = 0.3861, while CalcZAF gives 0.3904.

CalcZAF output:



My output:



Any thoughts?  Would you be willing to let me have a look at the older CITZAF code, which might be easier for me to pick through?  I’ve been tearing my hair out trying to figure out why I’m getting different results.

[Probeman]

Hi Brian,
Wow, you amaze me sir!     

This all reminds me very much of my early efforts when I took John Armstrong's QuickBasic code and translated it into Microsoft FORTRAN on one of the first IBM AT computers!   Ironically, 10 years later I translated all that FORTRAN back into Visual Basic to get it running as a native Windows application!

And here you are (with apparently more time than money just as I was!), translating the Visual Basic back into FORTRAN!  I really think you should just buy John Donovan's Probe for EPMA software for your JEOL 8230/8530 instrument and be done with it!   

Here are some suggestions:  First of all, you do know that you can run CalcZAF to reprocess bulk intensity data from a text file into an output file?  See the File | Open CalcZAF Input Data File menu...  The input file format is nicely documented in the Help file.

Second, the Phi code is in the ZAFPhiCal.  The complete CalcZAF physics code is on Github here:

https://github.com/openmicroanalysis/calczaf

Third, feel free to call me at my office or lab to chat. My phone numbers are at the bottom of this page:

http://epmalab.uoregon.edu/

I'm usually at UofO by 10 AM Pacific time.

[Probeman]

Hi Brian,
Two more thoughts:  if I were you I'd run a simple composition in CalcZAF in VerboseMode (it's already in DebugMode by default) and compare the intermediate calculations with your FORTRAN.  See the Output menu.

Second: you are aware that I offer a matrix correction API?  It's called Matrix (appropriately enough!), and is described here:

http://probesoftware.com/smf/index.php?board=8.0

Basically all the matrix correction options in calcZAF, including fast Monte-Carlo methods are available through this interface.  It can be called from any OLE (COM) compliant container such as Excel, MatLab, etc.

I'm here at UofO if you want to chat
john

[Brian Joy]

Hi John,

Thanks for the input.  I actually programmed it from scratch, as I wanted to familiarize myself with the different models.  This is why I’m running into the problem with consistency versus CalcZAF.
 
My other motivation for writing the program originally was to add in Claude Merlet’s X-Phi model (which seems to work very nicely for a wide range of compounds, but for which I only have an ancient standalone program written for Windows 3.1).  However, I haven’t been able to get the details of his atomic number correction (see, for instance, Mikrochimica Acta, v. 114/115, p. 366).

Yes, now I see the evaluations of the integrals of phi(rho*z) and phi(rho*z)*(-chi)*rho*z in lines 1250-1260 of ZAF.bas.  I overlooked these lines while getting frustrated trying to get used to the VB syntax.

I’ll try running CalcZAF in verbose mode (hadn’t thought about this) and see if I can identify what I’m doing wrong.  I wasn’t aware of the matrix correction API.  When I click on the forum link to the Matrix.msi installer, I end up with a “page not found” error.

I have in fact tried processing batch data with CalcZAF (as well as with CITZAF).  However, I wanted to be able to use the JEOL output with only the most minimal modifications and then produce my own output organized exactly the way I wanted it.

Thanks for the help.  Maybe I’ll give you a holler later this week or next after I have a chance to look through the CalcZAF code some more.

Brian

[Probeman]

I’ll try running CalcZAF in verbose mode (hadn’t thought about this) and see if I can identify what I’m doing wrong.  I wasn’t aware of the matrix correction API.  When I click on the forum link to the Matrix.msi installer, I end up with a “page not found” error.

I'll check the forum link to matrix.mdb...  but if you use this link it should work:

http://www.probesoftware.com/ftp/V11/Matrix.msi

But you will need a login and password. Contact me by email and I'll get that to you.
john

[John Donovan]

I’ll try running CalcZAF in verbose mode (hadn’t thought about this) and see if I can identify what I’m doing wrong.  I wasn’t aware of the matrix correction API.  When I click on the forum link to the Matrix.msi installer, I end up with a “page not found” error.

Hi Brian,
Thanks for catching that.  The update page got renamed from .htm to .html!   I fixed the link you mentioned and here is the corrected link:

http://www.probesoftware.com/Update.html

Darn computers... they do exactly what you tell them to do!   

[Brian Joy]

Hi John,

I figured out why my program produces different values of f(chi) for X-rays in mixtures than CalcZAF does when using the Armstrong absorption model:

When CalcZAF determines the mixture mean ionization energy for the purpose of calculating f(chi), it calls subroutine ZAFMip2 (in zaf.bas), and the following code is executed:

962      Sub ZAFMip2(ii As Integer, zz As Single, ww As Single)
970      ww! = (9.76 * zz! + 58.5 * (1# / zz! ^ 0.19)) / 1000#
1032    Sub ZAFPhiCal(zafinit As Integer)
1082    zz! = 0#
1088    For i1% = 1 To zaf.in0%
1089    zz! = zz! + zaf.conc!(i1%) * zaf.Z%(i1%) / zaf.atwts!(i1%)
1090    aa! = aa! + zaf.conc!(i1%)
1091    pp! = pp! + zaf.conc!(i1%) / zaf.atwts!(i1%)
1092    Next i1%
1094    zz! = zz! / pp!
1095    aa! = aa! / pp!
1100    Call ZAFMip2(i%, zz!, ww!)
1130    tx4! = Log(1.166 * v0! / ww!) / (v0! - v1!)
1132    x4! = 297000# * zz! ^ 1.05 / (aa! * v0! ^ 1.25) * (tx4!) ^ 0.5
1133    x3! = 850000# * zz! * zz! / (aa! * v0! * v0! * (X2! - 1#))

In contrast, when CalcZAF determines the mixture mean ionization energy for the purpose of calculating stopping power, it calls subroutine ZAFMip, and the following code is executed (for the Love-Scott model):

833     Sub ZAFMip(zafinit As Integer)
855     For i% = 1 To zaf.in0%
856     jm!(i%) = (9.76 * zaf.Z%(i%) + 58.5 * (1# / zaf.Z%(i%) ^ 0.19)) / 1000#
858     Next i%
929     m5! = 0#
930     m6! = 0#
931     For i% = 1 To zaf.in0%
932     m5! = m5! + zaf.conc!(i%) * zaf.Z%(i%) / zaf.atwts!(i%)
933     m6! = m6! + zaf.conc!(i%) * zaf.Z%(i%) * Log(jm!(i%)) / zaf.atwts!(i%)
934     Next i%
938     jbar! = Exp(m6! / m5!)
2774   Sub ZAFStp(zafinit As Integer)
2954   m6! = jbar!
2955   zaf.stp!(i%) = (1# + 16.05 * (m6! / zaf.eC!(i%)) ^ 0.5 * ((zaf.v!(i%) ^ 0.5 - 1) / (zaf.v!(i%) - 1)) ^ 1.07) / m5!
2956   zaf.stp!(i%) = 1# / zaf.stp!(i%)

The latter block of code is consistent with the work of Berger and Seltzer (1964, NASA SP-3012, p. 6), while the former is not.  When I use the former, i.e., when I substitute Sum(C(i)*Z(i)/A(i))/Sum(C(i)/A(i)) for Z in the equation of Berger and Seltzer when calculating f(chi) (and then use the latter expression for calculation of stopping power), then I get exactly the same results as produced by CalcZAF when using the Armstrong/Love-Scott model (at least for K lines).

So then I guess this begs the obvious question:  Why is the mixture mean ionization energy calculated differently in different parts of the program?

[Probeman]

Hi Brian,
This is very interesting.

As I mentioned above, besides adding some additional features, I pretty am pretty much utilizing John Armstrong's and Paul Carpenter's original code from their Cal Tech days.  Translating code from one language to another is tedious but it does allow for such "sanity checks".   In fact I found a few bugs in their code during such a translation process myself! . One of these bugs was  in the original FRAME fluorescence correction, and a few others from typos in the various published papers...

So this could very well be another bug or even a publication typo.  It appears to only make a difference in the 2nd or 3rd decimal place of the absorption correction, but we should try and get to the bottom of it. 

I've written John and Paul and hopefully they will chime in with their thoughts on this.

Thanks for noticing this.

[Brian Joy]

Hi John,

Thanks for checking on this.  I've been scratching my head over this discrepancy in my calculated Armstrong f(chi) versus CalcZAF for at least the past year.  I really appreciate the fact that you've posted the CalcZAF code, as now I can check my work.  Don't worry, I'll have more questions soon.

Like you note, for practical purposes, the discrepancy doesn't make much difference in the absorption correction, especially since the standard isn't likely to be a pure element (at least not in cases in which absorption is likely to be pronounced, e.g., X-rays of light elements).  I'm just trying to reproduce the model as precisely as possible simply for the sake of understanding the numbers.

Brian

[Probeman]

Hi Brian,
I spoke with John Armstrong this week and he told me he will be looking at your work this weekend and also checking the TryZAF/CalcZAF code.

He thinks it is very possible you are on to something, but that he has had to sort through so many typos in the original papers, that it's a non-trivial problem to reproduce some of this published work.

That said, I thought I would start by posting some "baseline" error distributions, so if we do end up making some changes to the code, we can document any degradation/improvement in these calculations.

So here are some calculations on the Bastin binary k-ratio database (the full database you sent me with 876 binaries and attached below) using the current release of CalcZAF.  Remember, this operation is easily performed using this menu in CalcZAF:



First the JTA correction using Henke MACs:



Not too bad, but with a 5.5% variance we can do better.  Let's try with FFAST MACs:



3.75% variance is an improvement.  Now let's try PAP, first with Henke:



Slightly better with a 3.72% variance. Now PAP with FFAST:



About the same but slightly worse (3.75% variance).  Interesting, but this is mainly to provide a visible "baseline" for any changes we make to the matrix correction codes...
john

[Probeman]

By the way, I find this little bump interesting:



Here are the largest outliers:

Problematic k-ratio errors (< 0.8 or > 1.2)
  Line             ConcA   ConcB     TOA      eO      Uo   K-Exp   K-Cal   K-Err
378   Al ka in Fe .100000 .900000 20.0000 35.0000 22.4359 .012600 .015186 1.20524
379   Al ka in Fe .100000 .900000 20.0000 40.0000 25.6410 .011400 .013768 1.20773
651   Al ka in Ni .049000 .951000 40.0000 26.6000 17.0513 .009800 .011836 1.20776
652   Al ka in Ni .049000 .951000 40.0000 31.9000 20.4487 .007500 .009277 1.23698

[Probeman]

Hi Brian,
In case you might be interested I've attached the GMRFILM FORTRAN source code (in a ZIP file) below (remember to log in to see it).

According to Paul Carpenter Waldo's source contains both an implementation of PAP and Proza, so maybe looking at this code will help you figure out the typos in the various publications.
john

[Brian Joy]

I figured out why I was getting different results for f(chi) of M lines using the Armstrong model compared to output from CalcZAF.  It is an issue of numerical precision.

Initially I had been using the following expressions to determine f(chi) for the Armstrong (modified Packwood-Brown) model:



However, for heavy elements, the beta parameter is large because it scales with the square of atomic number.  For instance, in the case of W Ma @ 15 kV, its value (units omitted) according to the Armstrong model is 48982 (and alpha = 2455.0).  Thus beta^2/(4*alpha^2) = 99.520, and exp(99.520) = 1.6634*10^43.  For this case, Fortran returns a value for ERF(beta/(2*alpha)) precisely equal to one, when in truth it should be very, very slightly less than one (beyond the limit of attainable precision).  Thus the expression for F(0) reduces to F(0) = (sqrt(Pi)*gamma0 )/(2*alpha), which is not correct.  Using ERFC instead of ERF produces a better result (since the returned value is not precisely zero), but evaluation of F(0) and F(chi) still requires working with multiplication of extremely large numbers by extremely small numbers.

The solution to the problem is to avoid exponentiation, as is done in CalcZAF.  This is achieved by using the approach of Bastin et al. (1984, X-ray Spectrometry 13:91-97), in which only the polynomial that multiplies the exponential factor in the ERFC approximation is evaluated, noting that the exponential factor in the ERFC approximation is the inverse of the exponential factor that multiplies each respective ERFC factor in the expressions for F(0) and F(chi).  (See the attached document if you want.)

Lesson learned.