Title                  :CLEDX
Keywords               :EDXS quantitative analysis, Cliff-Lorimer
Computer               :IBM PC family or clones
Operating System       :Microsoft DOS 3.x or higher
Programming Language   :Microsoft FORTRAN Optimizing compiler 4.10
                       :(large library, math-coprocessor emulator mode)
Hardware Requirements  :None
Author(s)              :Pierre TREBBIA
Correspondence Address :Laboratoire de Physique des Solides, Bat. 510
                       :F91405 ORSAY CEDEX FRANCE.
                       :Email : Bitnet address : TREBBIA@FRUPS51
 
Abstract:
 
From experimental conditions (angles, detection system, primary voltage,
line intensities & backgrounds), this program computes the weight &
atomic concentrations of up to 20 elements in a thin film X-Ray analysis
experiment using the Cliff-Lorimer procedure. An estimation of the
confidence intervals on the concentrations, for different first kind
risks, is also given. No special hardware is needed (i.e. math coprocessor
or graphic monitor). Results can be directed either to a lineprinter or to the
screen. From the original COLIJN & ROMIG program PCEDS, several improvements
are made concerning both the algorithms (parameterized maximum number
of elements to analyse, computation speed, analysis of the results
precision) and the accuracy of the formulae and data. Several misprints
were also corrected. Reference : "CLEDX : a stand-alone program for
quantitative X-Ray analysis of thin films using the Cliff-Lorimer
procedure", P. TREBBIA, Ultramicroscopy 27, 343-348 (1989).
Title                  :CLEDX
Keywords               :EDXS quantitative analysis, Cliff-Lorimer
Computer               :IBM PC family or clones
Operating System       :Microsoft DOS 3.x or higher
Programming Language   :Microsoft FORTRAN Optimizing compiler 4.10
                       :(large library, math-coprocessor emulator mode)
Hardware Requirements  :None
Author(s)              :Pierre TREBBIA
Correspondence Address :Laboratoire de Physique des Solides, Bat. 510
                       :F91405 ORSAY CEDEX FRANCE.
                       :Email : Bitnet address : TREBBIA@FRUPS51
 
Documentation:
 
CLEDX is a self-documented FORTRAN program. After compiling and
linking with the Microsoft Optimizing Compiler 4.10, you should
obtain an executable task CLEDX.EXE. Running this task will prompt
you with a lot of explicite questions containing always a suggested
(default) answer. No special user-guide is needed for using CLEDX.
Of course, I strongly suggest to test CLEDX using the following
output example : you must obtain the same results provided you enter
the same set-up and experimental data.
In case of problem, please refer to Ultramicroscopy 27, p.343-348
(1989), and/or print the listing of CLEDX : a lot of valuable
comments are included into the source code giving indications about
both the algorithms and the origin of the formulae and data.
 
Here is now a typical output from CLEDX :
 
                      DEMONSTRATION OF PROGRAM CLEDX
 
                       **** Angular Conditions ****
 
    Specimen tilt angles (deg.) =  25.00 18.00
    Detector elevation angle (deg.) = 33.00
    Detector azimuthal angle (deg.) =  90.00
    Detector inclination angle (deg.) =  30.00
    X-ray take-off angle (deg.) =  58.00
 
Global tilt (deg.) :  25.00
Real foil thickness = Apparent foil thickness *   0.91
 
                *****    AEM Thin Film EDX Analysis    *****
 
                      done on : 10-DEC-88 at 14:04:10
 
                  6 elements.      1 independant analysis.
 
                              *** CAUTION ***
     Remember that this is a compound in which ALL the elements
     are perhaps NOT detectable...: Carbide, Nitride, Oxide ...
 
Elements : As Mo Si Pr Bi Cd
 
 
Operating kV =  100.0
Effective carbon contamination (micr.) =  0.0000
Effective Be window (micr.) = 10.0
Effective Au layer (Angs.) =   300.
Effective dead layer (Angs.) =  3004.
Effective active layer (cm) =  0.300
Calculated C-L factors
NO absorption correction.
No user supplied bulk density
 
 
                   ***** Experiment number :   1 *****
 
 
                   Title : Demonstration output
 
 
Element  Atomic no.  Line type  Integrated Intensity  Integrated Background
 
  As        33        K alpha                 100.                  10.
  Mo        42        K beta                  200.                  15.
  Si        14        K total                 300.                  30.
  Pr        59        L alpha                 250.                  20.
  Bi        83        L beta                  300.                  28.
  Cd        48        L total                 150.                  12.
 
              *** Results of Cliff-Lorimer algorithm ***
 
                        beta = v/c = 0.5482
 
 
                            Element : As
Atomic nb : 33      Atomic weight :  74.922      X-Ray line type : K alpha
Mean wavelength (Ang.) :  1.1774        Mean X-Ray energy (keV) : 10.530
Ionization energy (keV) :  11.8667      Overvoltage :   6.4709
Cross section (cm2) :   0.1786E-21      Fluo yield : 0.5683
Weight of the line : 0.8640             Detection efficiency : 0.9914
                      Overall C-L factor :   1.0000
 
 
                            Element : Mo
Atomic nb : 42      Atomic weight :  95.940      X-Ray line type : K beta
Mean wavelength (Ang.) :  0.6308        Mean X-Ray energy (keV) : 19.653
Ionization energy (keV) :  19.9995      Overvoltage :   3.8395
Cross section (cm2) :   0.8685E-22      Fluo yield : 0.7638
Weight of the line : 0.1649             Detection efficiency : 0.9468
                      Overall C-L factor :  10.7459
 
 
                            Element : Si
Atomic nb : 14      Atomic weight :  28.086      X-Ray line type : K total
Mean wavelength (Ang.) :  7.1170        Mean X-Ray energy (keV) :  1.742
Ionization energy (keV) :   1.8389      Overvoltage :  41.7579
Cross section (cm2) :   0.2581E-20      Fluo yield : 0.0430
Weight of the line : 1.0000             Detection efficiency : 0.7335
                      Overall C-L factor :   0.4002
 
 
                            Element : Pr
Atomic nb : 59      Atomic weight : 140.908      X-Ray line type : L alpha
Mean wavelength (Ang.) :  2.4643        Mean X-Ray energy (keV) :  5.031
Ionization energy (keV) :   5.9643      Overvoltage :  12.8747
Cross section (cm2) :   0.1395E-20      Fluo yield : 0.1180
Weight of the line : 0.4823             Detection efficiency : 0.9385
                      Overall C-L factor :   2.1941
 
 
                            Element : Bi
Atomic nb : 83      Atomic weight : 208.980      X-Ray line type : L beta
Mean wavelength (Ang.) :  0.9519        Mean X-Ray energy (keV) : 13.025
Ionization energy (keV) :  15.7111      Overvoltage :   4.8875
Cross section (cm2) :   0.4765E-21      Fluo yield : 0.3730
Weight of the line : 0.4960             Detection efficiency : 0.9866
                      Overall C-L factor :   2.7877
 
 
                            Element : Cd
Atomic nb : 48      Atomic weight : 112.411      X-Ray line type : L total
Mean wavelength (Ang.) :  3.8565        Mean X-Ray energy (keV) :  3.215
Ionization energy (keV) :   3.5375      Overvoltage :  21.7070
Cross section (cm2) :   0.3006E-20      Fluo yield : 0.0560
Weight of the line : 1.0000             Detection efficiency : 0.8389
                      Overall C-L factor :   0.9234
 
 
              ***** Solution UNCORRECTED for absorption *****
 
 
                   **** Concentration Results ****
 
 
     *** Warning : elements Mo and Cd have been interchanged ***
 
 
              Estimation of Weight Concentrations
              -----------------------------------
 
Elem      Peak    Background    h       Var(P)      Wght %      Sigma(Wght %)
 
 As        100.         10.    2.0        120.       2.5690       0.2809
 Cd        150.         12.    2.0        174.       3.5585       0.3122
 Si        300.         30.    2.0        360.       3.0844       0.1949
 Pr        250.         20.    2.0        290.      14.0913       0.9380
 Bi        300.         28.    2.0        356.      21.4845       1.2896
 Mo        200.         15.    2.0        230.      55.2123       1.5919
 
  Estimation of covariances
  -------------------------
Covariance(As,Cd) =   -0.37849E-07
Covariance(As,Si) =   -0.14708E-07
Covariance(As,Pr) =   -0.35610E-06
Covariance(As,Bi) =   -0.70568E-06
Covariance(Cd,Si) =   -0.18186E-07
Covariance(Cd,Pr) =   -0.44032E-06
Covariance(Cd,Bi) =   -0.87257E-06
Covariance(Si,Pr) =   -0.17110E-06
Covariance(Si,Bi) =   -0.33907E-06
Covariance(Pr,Bi) =   -0.82096E-05
 
  Estimation of confidence intervals assuming NO ERROR on C-L coeff.
  ------------------------------------------------------------------
 
    At the confidence level alpha = 0.10 :
 
Element  As Interval :  2.569  + -  0.462
Element  Cd Interval :  3.558  + -  0.514
Element  Si Interval :  3.084  + -  0.321
Element  Pr Interval : 14.091  + -  1.543
Element  Bi Interval : 21.485  + -  2.121
Element  Mo Interval : 55.212  + -  2.619
 
    At the confidence level alpha = 0.05 :
 
Element  As Interval :  2.569  + -  0.551
Element  Cd Interval :  3.558  + -  0.612
Element  Si Interval :  3.084  + -  0.382
Element  Pr Interval : 14.091  + -  1.839
Element  Bi Interval : 21.485  + -  2.528
Element  Mo Interval : 55.212  + -  3.120
 
    At the confidence level alpha = 0.01 :
 
Element  As Interval :  2.569  + -  0.723
Element  Cd Interval :  3.558  + -  0.804
Element  Si Interval :  3.084  + -  0.502
Element  Pr Interval : 14.091  + -  2.415
Element  Bi Interval : 21.485  + -  3.321
Element  Mo Interval : 55.212  + -  4.099
 
Mu/ro(lamda(As),As) coeff :      37.8399
Mu/ro(lamda(As),Mo) coeff :      76.0130
Mu/ro(lamda(As),Si) coeff :      28.9956
Mu/ro(lamda(As),Pr) coeff :     194.6644
Mu/ro(lamda(As),Bi) coeff :     121.4439
Mu/ro(lamda(As),Cd) coeff :     110.5914
Mu/ro(lamda(Mo),As) coeff :      50.9196
Mu/ro(lamda(Mo),Mo) coeff :      13.9235
Mu/ro(lamda(Mo),Si) coeff :       5.1481
Mu/ro(lamda(Mo),Pr) coeff :      36.5585
Mu/ro(lamda(Mo),Bi) coeff :      94.3361
Mu/ro(lamda(Mo),Cd) coeff :      20.3841
Mu/ro(lamda(Si),As) coeff :    5141.5825
Mu/ro(lamda(Si),Mo) coeff :    1481.8705
Mu/ro(lamda(Si),Si) coeff :     322.7175
Mu/ro(lamda(Si),Pr) coeff :    4334.4004
Mu/ro(lamda(Si),Bi) coeff :    1918.2653
Mu/ro(lamda(Si),Cd) coeff :    2228.8499
Mu/ro(lamda(Pr),As) coeff :     284.2106
Mu/ro(lamda(Pr),Mo) coeff :     566.7220
Mu/ro(lamda(Pr),Si) coeff :     224.3123
Mu/ro(lamda(Pr),Pr) coeff :     275.0118
Mu/ro(lamda(Pr),Bi) coeff :     828.6407
Mu/ro(lamda(Pr),Cd) coeff :     818.4572
Mu/ro(lamda(Bi),As) coeff :     153.9934
Mu/ro(lamda(Bi),Mo) coeff :      42.6310
Mu/ro(lamda(Bi),Si) coeff :      16.0899
Mu/ro(lamda(Bi),Pr) coeff :     110.1079
Mu/ro(lamda(Bi),Bi) coeff :      69.8706
Mu/ro(lamda(Bi),Cd) coeff :      62.1560
Mu/ro(lamda(Cd),As) coeff :     965.2724
Mu/ro(lamda(Cd),Mo) coeff :    1916.1725
Mu/ro(lamda(Cd),Si) coeff :     775.6093
Mu/ro(lamda(Cd),Pr) coeff :     881.2006
Mu/ro(lamda(Cd),Bi) coeff :    1648.9332
Mu/ro(lamda(Cd),Cd) coeff :     453.1339
 
Mu/ro(lamda(As),specimen) coeff :     101.2927
Mu/ro(lamda(Mo),specimen) coeff :      35.2990
Mu/ro(lamda(Si),specimen) coeff :    2062.4331
Mu/ro(lamda(Pr),specimen) coeff :     573.0274
Mu/ro(lamda(Bi),specimen) coeff :      60.7288
Mu/ro(lamda(Cd),specimen) coeff :    1601.2467
 
         The calculated bulk density is  8.42 g/cc
 
 
  If the foil, at this analysis site, is greater than     497.Angstroms,
   an absorption correction may be required.
 
 
                    Absorption correction requested.
                    Apparent foil thickness (Angs.) :   800.
 
 
               ***** Solution CORRECTED for absorption *****
 
 
 
Iteration #   1
Mu/ro(lamda(As),specimen) coeff :     101.2927
Mu/ro(lamda(Mo),specimen) coeff :      35.2990
Mu/ro(lamda(Si),specimen) coeff :    2062.4331
Mu/ro(lamda(Pr),specimen) coeff :     573.0274
Mu/ro(lamda(Bi),specimen) coeff :      60.7288
Mu/ro(lamda(Cd),specimen) coeff :    1601.2467
    Absorption Correction (As,Mo):   1.002
    Absorption Correction (As,Si):   0.932
    Absorption Correction (As,Pr):   0.983
    Absorption Correction (As,Bi):   1.001
    Absorption Correction (As,Cd):   0.947
 
                   **** Concentration Results ****
 
 
     *** Warning : elements Mo and Cd have been interchanged ***
 
 
              Estimation of Weight Concentrations
              -----------------------------------
 
Elem      Peak    Background    h       Var(P)      Wght %      Sigma(Wght %)
 
 As        100.         10.    2.0        120.       2.5562       0.2795
 Cd        150.         12.    2.0        174.       3.7372       0.3278
 Si        300.         30.    2.0        360.       3.2935       0.2081
 Pr        250.         20.    2.0        290.      14.2611       0.9485
 Bi        300.         28.    2.0        356.      21.3459       1.2814
 Mo        200.         15.    2.0        230.      54.8060       1.5933
 
  Estimation of covariances
  -------------------------
Covariance(As,Cd) =   -0.41827E-07
Covariance(As,Si) =   -0.16803E-07
Covariance(As,Pr) =   -0.36544E-06
Covariance(As,Bi) =   -0.69796E-06
Covariance(Cd,Si) =   -0.23147E-07
Covariance(Cd,Pr) =   -0.50342E-06
Covariance(Cd,Bi) =   -0.96149E-06
Covariance(Si,Pr) =   -0.20223E-06
Covariance(Si,Bi) =   -0.38624E-06
Covariance(Pr,Bi) =   -0.84004E-05
 
  Estimation of confidence intervals assuming NO ERROR on C-L coeff.
  ------------------------------------------------------------------
 
    At the confidence level alpha = 0.10 :
 
Element  As Interval :  2.556  + -  0.460
Element  Cd Interval :  3.737  + -  0.539
Element  Si Interval :  3.294  + -  0.342
Element  Pr Interval : 14.261  + -  1.560
Element  Bi Interval : 21.346  + -  2.108
Element  Mo Interval : 54.806  + -  2.621
 
    At the confidence level alpha = 0.05 :
 
Element  As Interval :  2.556  + -  0.548
Element  Cd Interval :  3.737  + -  0.642
Element  Si Interval :  3.294  + -  0.408
Element  Pr Interval : 14.261  + -  1.859
Element  Bi Interval : 21.346  + -  2.511
Element  Mo Interval : 54.806  + -  3.123
 
    At the confidence level alpha = 0.01 :
 
Element  As Interval :  2.556  + -  0.720
Element  Cd Interval :  3.737  + -  0.844
Element  Si Interval :  3.294  + -  0.536
Element  Pr Interval : 14.261  + -  2.442
Element  Bi Interval : 21.346  + -  3.300
Element  Mo Interval : 54.806  + -  4.103
 
 
              Estimation of Atomic Concentrations
              -----------------------------------
 
Element : As   At. Wght. :  74.922   Density :   5.727   Atomic % :  3.5568
Element : Mo   At. Wght. :  95.940   Density :  10.200   Atomic % : 59.5528
Element : Si   At. Wght. :  28.086   Density :   2.330   Atomic % : 12.2252
Element : Pr   At. Wght. : 140.908   Density :   6.773   Atomic % : 10.5510
Element : Bi   At. Wght. : 208.980   Density :   9.800   Atomic % : 10.6484
Element : Cd   At. Wght. : 112.411   Density :   8.642   Atomic % :  3.4659
 
 
    At this analysis point, the apparent foil thickness is  800. Angstroms
             The mean X-Ray path length is  427. Angstroms.
                  The bulk density is  8.425 g/cc.
           **********************************************
Title                  :CLEDX
Keywords               :EDXS quantitative analysis, Cliff-Lorimer
Computer               :IBM PC family or clones
Operating System       :Microsoft DOS 3.x or higher
Programming Language   :Microsoft FORTRAN Optimizing compiler 4.10
                       :(large library, math-coprocessor emulator mode)
Hardware Requirements  :None
Author(s)              :Pierre TREBBIA
Correspondence Address :Laboratoire de Physique des Solides, Bat. 510
                       :F91405 ORSAY CEDEX FRANCE.
                       :Email : Bitnet address : TREBBIA@FRUPS51
SOURCE CODE:
--------------

      PROGRAM CLEDX
c
c     LAST VERSION 9-DEC-1988 14:56:24
c
c     purpose
c     -------
c       From experimental conditions (angles, detection system, primary voltage,
c     line intensities & backgrounds), this program computes the weight & atomic
c     concentrations of up to NMAX elements in a thin film X-Ray analysis
c     experiment using the Cliff-Lorimer procedure. An estimation of the
c     confidence intervals on the concentrations, for different first kind
c     risks, is also given.
c
c     subroutines and function subprograms required
c     ---------------------------------------------
c     CONFIG,TYPLIN,CLIFF,QUANT,THICK,ABCOR,DENSE,ATOMIC
c
c     comments
c     --------
c       This program is the heart of a set of 14 FORTRAN/77 program, functions
c     and subroutines. They have to be compiled and linked together.
c       There are 2 versions of this program : one for DEC computers (VAX or
c     LSI 11/73 for example), and one for IBM PC or compatible computers.
c       Depending on the compiler used to get an executive task, the maximum
c     number NMAX of elements the program can analyze at once can vary. This is
c     mainly due to the internal representation of the double precision real
c     variable DET in the subroutine INVMAT : for NMAX greater than 10, DET
c     can overflow 10**38. Such a problem does not occur with the PC version
c     and the Microsoft Compiler FORTRAN 4.1 (at least 20 elements can be
c     analyzed).
c       No special hardware is needed (i.e. math coprocessor or graphic
c     monitor). Results can be directed either to a lineprinter or to the
c     screen.
c       From the original COLIJN & ROMIG program PCEDS, several improvements
c     are made concerning both the algorithms (parameterized maximum number
c     of elements to analyse, computation speed, analysis of the results
c     precision) and the accuracy of the formulae and data. Several misprints
c     were also corrected.
c       For a detailed explanation of these improvements, please refer to the
c     following paper : "CLEDX : a stand-alone program for quantitative X-Ray
c     analysis of thin films using the Cliff-Lorimer procedure", P. TREBBIA,
c     submitted for publication in Ultramicroscopy (december 1988).
c
c       IN ORDER TO MODIFY THE MAXIMUM NUMBER OF ELEMENTS TO ANALYZE, MODIFY
c     THE FIRST STATEMENT "PARAMETER ...." IN EACH OF THE 14 PROGRAMS.
c
c     authors :
c     ---------
c       It's a super-coproduction involving :
c
c     1) A. D. ROMIG, JR., DIVISION 1832, SANDIA NATIONAL LABORATORIES,
c        ALBUQUERQUE, NEW MEXICO 87185 for the general scheme of the program
c        (written for VAX computers),
c     2) H. O. COLIJN, OHIO STATE UNIVERSITY, CENTRAL ELECTRON OPTICS FACILITY,
c        COLUMBUS, OH 43210 for the translation of Romig's original program to
c        PC computers and the inclusion of the CONFIG file,
c     3) N. ZALUZEC, ARGONNE NATIONAL LABORATORY, MATERIALS SCIENCE &
c        TECHNOLOGY DIVISION, ARGONNE, IL 60439 for accurate formulae,
c     4) Pierre TREBBIA, US 41 : "Microscopie Electronique Analytique
c        Quantitative", CNRS, Laboratoire de Physique des Solides, Bat. 510
c        Universite Paris-Sud, F91405 ORSAY Cedex, Phone : (33-1) 69 41 53 68
c        for a complete revision of the program PCEDS and algorithms.
c
      parameter (nmax=20)
      real l1,l2,l3,m1,m2,m3,m4,m5
      dimension ek(82),l1(82),L2(82),L3(82),m1(82),m2(82),m3(82),
     &m4(82),m5(82),ro(82)
      common /energ/ek,l1,l2,l3,m1,m2,m3,m4,m5
      common /dens/ro
      DIMENSION wconc(nmax),aconc(nmax),wnew(nmax),XI(nmax),BACK(nmax)
      common /conc/wconc,aconc
      DIMENSION na(nmax),li(nmax),alfw(nmax),CLKR(nmax),hegert(nmax)
      common /para1/ne,na,li
      common /para2/alfw,clkr
      DIMENSION ABCOF(nmax),ELABCF(nmax,nmax),ATWT(82),COR(nmax)
      COMMON /COEF/ ABCOF,elabcf
      COMMON /WT/ ATWT
      character lab(2,82)
      character*164 symb
      equivalence(symb,lab(1,1))
      common /label/lab
      common /ilp/ilp
      logical full,last,clflag,abflag,dnflag,muflag,limit,muprnt,denchk
      common /flags/full,last,clflag,abflag,dnflag,muflag,limit,muprnt,
     &              denchk
      REAL KV,LAMDA,lam
      DIMENSION LAMDA(10,82),lam(6,nmax),reli(5,nmax),factor(nmax)
      COMMON /PARAM/ KV
      COMMON /WAVE/ lam,reli
      common /lamda/lamda
      dimension relint(10,82)
      common /relint/relint
      COMMON /CONST/ BULKRO,PSIR
      COMMON /WIND/ WTH,CTH,AUTH,SITH,siacth
      COMMON /HOC/ ALPHA, PHI,detang
      INTEGER F12,F13,F23,W3
      DIMENSION NK3(33),w3(65),f12(65),f13(65),f23(65)
      common /nestor/nk3,w3,f12,f13,f23
      dimension wk1(7)
      common /wk/wk1
      CHARACTER CHAR,char2
      character btitle(40)
      character*7 lintyp
      common /lintyp/lintyp
      logical ready
c
c     table of the element symbols
c     ----------------------------
c
      symb(1:28)='NaMgAlSiP S ClArK CaScTiV Cr'
      symb(29:72)='MnFeCoNiCuZnGaGeAsSeBrKrRbSrY ZrNbMoTcRuRhPd'
      symb(73:116)='AgCdInSnSbTeI XeCsBaLaCePrNdPmSmEuGdTbDyHoEr'
      symb(117:164)='TmYbLuHfTaW ReOsIrPtAuHgTlPbBiPoAtRnFrRaAcThPaU '
c
c     set-up the experimental conditions
c     ----------------------------------
c
      CALL CONFIG
      ready=.true.
      char='Y'
      write(*,10) char
10    format(//' Are these parameters O.K. ? (Y/N) <',a1,'> ',$)
      read(*,21) char
      if((char.eq.'N').or.(char.eq.'n')) ready=.false.
c
c     set-up the analysis conditions
c     ------------------------------
c
      char='Y'
      ILP = 0
      WRITE (*,20) char
20    FORMAT (/,' Results on the printer ? (Y/N)<',a1,'>: ',$)
      READ (*,21) CHAR
21    FORMAT (A1)
      IF (.not.(CHAR .EQ. 'N'.OR. CHAR .EQ. 'n')) then
       ILP = 6
       OPEN (6, file='prn',STATUS='old')
      endif
      char='Y'
      full=.true.
      muprnt=.false.
      write (*,25) char
25    format(' Do you want a full report ? (Y/N)<',a1,'>: ',$)
      read(*,21) char
      if (char.eq.'N'.or.char.eq.'n') full=.false.
      if(full) then
       char='N'
       write(*,26) char
26     format(5x,'Including the mu/ro values ? (Y/N)<',a1,'>: ',$)
       read(*,21) char
       if (char.eq.'Y'.or.char.eq.'y') muprnt=.true.
      endif
30    nd=1
      WRITE (*,35) nd
35    FORMAT (/,' Number of independant experiments <',i1,'>: ',$)
      n = 0
      READ (*,36) n
36    format(i2)
      IF (n.NE. 0) nd=n
38    WRITE (*,40) nmax,NE
40    FORMAT(' Number of elements in the analysis (2-',i2,')<',i2,
     &'>: ',$)
      n = 0
      READ (*,36) n
c
c     check consistency
c
      if(n.gt.nmax) go to 38
c
c     resume
c
      if(n.ne.0) ne=n
c
c     check the experimental conditions & adjust the units
c     ----------------------------------------------------
c
      if(.not.ready) then
       WRITE (*,45) KV
45     FORMAT(//,' Operating kV <',F4.0,'>: ',$)
       var = 0.
       READ (*,46) var
46     format(f12.0)
       IF(var.NE. 0.) KV = var
      endif
      TILTX = 0.
      TILTY=0.
      WRITE (*,47) tiltx,tilty
47    FORMAT(/,' Enter now the 2 tilt angles (X,Y).',/,
     &' If you have only 1 tilt axis, answer X,0.',/,
     &' Specimen tilts (deg.)<',f3.0,',',f3.0,'>: ',$)
      vAR1 = 0.
      var2 = 0.
      READ (*,48) var1,var2
48    format(2f12.0)
      IF(var1 .NE. 0.) TILTX = var1
      if(var2.ne.0.) tilty=var2
      if(.not.ready) then
       WRITE (*,53) ALPHA
53     FORMAT(/,' You have defined the detector elevation angle as ',
     & f4.1,' deg.',/,15x,'I need confirmation :',/,
     & ' Detector elevation angle (deg.) : ',$)
       READ (*,46) alpha
       WRITE (*,55) PHI
55     FORMAT(/,' You have defined the detector azimuthal angle as ',
     & f4.1,' deg.',/,15x,'I need confirmation :',/,
     & ' Detector azimuthal angle (deg.) : ',$)
       READ (*,46) phi
       WRITE (*,57) detang
57     FORMAT(/,' You have defined the detector inclination angle as ',
     & f4.1,' deg.',/,15x,'I need confirmation :',/,
     & ' Detector inclination angle (deg.) : ',$)
       READ (*,46) detang
      endif
c
c     convert angles & calculate take-off angle
c     -----------------------------------------
c
      write(ilp,59)
59    format(//,24x,'**** Angular Conditions ****',/)
      if (full) then
       WRITE (ILP,60) TILTx,tilty
60     FORMAT(5X,'Specimen tilt angles (deg.) = ',2(F6.2))
       WRITE (ILP,61) ALPHA
61     FORMAT(5X,'Detector elevation angle (deg.) = ',F5.2)
       WRITE (ILP,62) PHI
62     FORMAT(5X,'Detector azimuthal angle (deg.) = ',F6.2)
       WRITE (ILP,63) detang
63     FORMAT(5X,'Detector inclination angle (deg.) = ',F6.2)
      endif
      PI = 3.1415926
      CONV = PI/180.
      PHIR = PHI*CONV
      THETXR = TILTX*CONV
      thetyr = tilty*conv
      ALPHAR = ALPHA*CONV
      var=((cos(alphar))**2)*((sin(thetyr)*cos(phir)+sin(thetxr)*
     &    cos(thetyr)*sin(phir))**2+(cos(thetxr)*cos(thetyr))**2)
      var=(cos(alphar))*(sin(thetyr)*cos(phir)+sin(thetxr)*
     &    cos(thetyr)*sin(phir))/(sqrt(var))
      tiltr=atan(var/sqrt(1.-var*var))
      tilt=tiltr/conv
      itilt=nint(tilt)
      psir=alphar+tiltr
      PSI    = PSIR/CONV
      WRITE (ILP,64) PSI
64    FORMAT(5X,'X-ray take-off angle (deg.) = ',F6.2,/)
      if((itilt.ne.0).and.full) then
       corth=cos(tiltr)
       write(ilp,65) tilt,corth
65     format(' Global tilt (deg.) : ',f6.2,/,' Real foil ',
     & 'thickness = Apparent foil thickness * ',f6.2)
      endif
c
c     check the detection conditions & adjust the units
c     -------------------------------------------------
c
      if(.not.ready) then
       WRITE (*,70) CTH
70     FORMAT (/,' You have defined the contamination layer thickness ',
     & 'as ',f4.1,' micr.',/,15x,'I need confirmation :',/,
     & ' Contamination layer thickness (micr.) : ',$)
       READ (*,71) cth
71     format(f12.4)
       WRITE (*,72) WTH
72     FORMAT (/,' You have defined the Be window thickness as ',
     & f4.1,' micr.',/,15x,'I need confirmation :',/,
     & ' Be window thickness (micr.) : ',$)
       READ (*,71) wth
       WRITE(*,73) AUTH
73     FORMAT(/' Gold contact layer thickness (Angs.)<',F5.0,'>: ',$)
       var = 0.
       READ (*,46) var
       IF (var .NE. 0.) AUTH = var
       WRITE (*,75) SITH
75     FORMAT(' Si dead layer thickness (Angs.)<',F6.0,'>: ',$)
       vAR = 0.
       READ (*,46) var
       IF (var .NE. 0.) SITH = var
       WRITE (*,77) siacth
77     FORMAT(' Si active layer thickness (cm)<',F7.4,'>: ',$)
       vAR = 0.
       READ (*,71) var
       IF (var .NE. 0.) siacth = var
      endif
      CTH=CTH*1.0E-4
      WTH=WTH*1.0E-4
      AUTH = AUTH * 1.0E-8
      SITH = SITH * 1.0E-8
c
c     compute the effective thicknesses of the detector
c     -------------------------------------------------
c
      if(detang.ne.alpha) then
       diff=abs(detang-alpha)
       diffr=diff*conv
       xfac=1./cos(diffr)
       cth=cth*xfac
       wth=wth*xfac
       auth=auth*xfac
       sith=sith*xfac
       siacth=siacth*xfac
       write(*,79) xfac
79     format(/,' *** All the detector thicknesses are multiplied by ',
     & f6.2,' ***',//)
      endif
c
c     input the elements to be analyzed
c     ---------------------------------
c
80    write(*,81)
81    format(/,15x,'*** YOUR ATTENTION PLEASE ***',//,15x,
     & 'For each input element,',/,' *) you MAY have to define ',
     & 'the line identifier (K or L)',/,' *) you MUST define the ',
     & 'line type :',/,5X,'A for alpha',/,5x,'B for beta',/,5X,
     & 'T for total if A and B are unresolved',/)
      DO 102 I=1,NE
82     WRITE (*,83) I
83     FORMAT(/,' Atomic no. of element #',i2,': ',$)
       na(i)=0
       READ (*,36) NA(I)
c
c      check consistency
c
       if((na(i).lt.11).or.(na(i).gt.92)) go to 82
c
c      offset to first element Na
c
       na(i)=na(i)-10
c
c      ask for the line type
c
       write(*,92) (lab(j,na(i)),j=1,2)
92     format(25x,'*** Element : ',2a1,' ***')
94     if(na(i).le.17) then
        li(i)=3
       else
        if(na(i).gt.61) then
         li(i)=6
        else
95       char='K'
         li(i)=3
         WRITE (*,96) char
96       FORMAT (' Line Identifier (K/L)<',a1,'>: ',$)
         READ (*,21) CHAR
         IF (CHAR.EQ.'L'.OR.CHAR.EQ.'l') then
          LI(I)=6
         endif
c
c        check consistency
c
         if(((li(i).eq.3).and.(na(i).gt.61)).or.
     &      ((li(i).eq.6).and.(na(i).lt.18))) go to 95
        endif
       endif
c
c      ask for the weight of the selected line
c
       if(li(i).eq.3) then
        char2='K'
       else
        char2='L'
       endif
       char='T'
       alfw(i) = 1.
       write(*,98) char2,char
98     format(' Which part of the ',a1,' line (A,B,T)<',a1,
     & '>: ',$)
       read (*,21) char
       if(char.eq.'A'.or.char.eq.'a') then
        alfw(i)=0.
        li(i)=li(i)-2
       endif
       if(char.eq.'B'.or.char.eq.'b') then
        alfw(i)=0.
        li(i)=li(i)-1
       endif
c
c      put the wavelengths in LAM
c
       do 99 j=1,5
99      lam(j,i)=0.
       if((li(i).eq.1).or.(li(i).eq.3)) then
        lam(1,i)=lamda(1,na(i))
         reli(1,i)=relint(1,na(i))
        lam(2,i)=lamda(2,na(i))
         reli(2,i)=relint(2,na(i))
       endif
       if((li(i).eq.2).or.(li(i).eq.3)) then
        lam(3,i)=lamda(3,na(i))
         reli(3,i)=relint(3,na(i))
        lam(4,i)=lamda(4,na(i))
         reli(4,i)=relint(4,na(i))
        lam(5,i)=lamda(5,na(i))
         reli(5,i)=relint(5,na(i))
       endif
       if((li(i).eq.4).or.(li(i).eq.6)) then
        lam(1,i)=lamda(6,na(i))
         reli(1,i)=relint(6,na(i))
        lam(2,i)=lamda(7,na(i))
         reli(2,i)=relint(7,na(i))
       endif
       if((li(i).eq.5).or.(li(i).eq.6)) then
        lam(3,i)=lamda(8,na(i))
         reli(3,i)=relint(8,na(i))
        lam(4,i)=lamda(9,na(i))
         reli(4,i)=relint(9,na(i))
        lam(5,i)=lamda(10,na(i))
         reli(5,i)=relint(10,na(i))
       endif
c
c      compute mean X-ray wavelength and put it in LAM(6,i)
c
       var1=0.
       var2=0.
       do 100 j=1,5
        if((lam(j,i).ne.0.).or.(reli(j,i).ne.0.)) then
         var1=var1+lam(j,i)*reli(j,i)
         var2=var2+reli(j,i)
        endif
100    continue
102    lam(6,i)=var1/var2
c
c     identifiers OK ?
c     ----------------
c
      char='Y'
      WRITE (*,103) char
103   FORMAT(5X,/,' Identifiers OK ? (Y/N)<',a1,'>',$)
      READ (*,21) char
      IF(char.EQ. 'N'.OR. char .EQ. 'n') THEN
105    WRITE (*,108)
108    FORMAT(//,' Re-Enter the data.',/)
       GO TO 80
      ENDIF
c
c     set-up the algorithms
c     ---------------------
c
      char='N'
      clflag=.true.
      WRITE (*,109) char
109   FORMAT (/,' User supplied Cliff-Lorimer factors ? (Y/N)<',a1,
     &'>: ',$)
      READ (*,21) CHAR
      IF (CHAR .EQ. 'Y'.OR. CHAR .EQ. 'y') clflag=.false.
      IF(.not.clflag) then
       clkr(1)=1.
       DO 111 I=2,NE
        WRITE (*,110) (lab(j,na(1)),j=1,2),(lab(j,na(i)),j=1,2)
110     FORMAT (' C-L Factor (',2a1,',',2a1,'): ',$)
111    READ (*,71) CLKR(i)
      endif
      char='N'
      abflag=.false.
      WRITE (*,112) char
112   FORMAT (' Correct for absorption ? (Y/N)<',a1,'>: ',$)
      READ (*,21) CHAR
      IF (CHAR .EQ. 'Y'.OR. CHAR .EQ. 'y') abflag=.true.
      if(abflag) then
       if(itilt.eq.0) then
        WRITE (*,113)
113     FORMAT (' Foil thickness (Angs.):  ',$)
       else
        WRITE (*,114)
114     FORMAT (' Apparent foil thickness (Angs.):  ',$)
       endif
       READ (*,46) TH
       THC=TH * (1.0E-08)
      endif
      char='N'
      dnflag=.true.
      WRITE (*,115) char
115   FORMAT (' User supplied density ? (Y/N)<',a1,'>: ',$)
      READ (*,21) CHAR
      IF (CHAR .EQ. 'Y'.OR. CHAR .EQ. 'y') THEN
        dnflag=.false.
        WRITE (*,117)
117     FORMAT(' Enter density (g/cc): ',$)
        READ (*,71) BULKRO
      ENDIF
c
c     everything OK ?
c     ---------------
c
      char='Y'
      WRITE (*,130) char
130   FORMAT(5X,/,' All data OK ? (Y/N)<',a1,'>',$)
      READ (*,21) char
      IF(char.EQ. 'N'.OR. char .EQ. 'n') THEN
       WRITE (*,108)
       GO TO 30
      ENDIF
c
c     abstract
c     --------
c
      WRITE (ILP,190)
190   format(/,17x,'*****    AEM Thin Film EDX Analysis    *****')
      WRITE (ILP,200) NE,ND
200   FORMAT(/,18X,I2,' elements.',5X,I2,' independant analysis.')
      WRITE(ILP,259)
259   FORMAT(/,31x,'*** CAUTION ***',/,6x,'Remember that this is a ',
     &'compound in which ALL the elements',/,6x,
     &'are perhaps NOT detectable...: Carbide, Nitride, Oxide ...',/)
      write(ilp,262) ((lab(j,NA(I)),j=1,2),i=1,ne)
262   format(' Elements :',23(' ',2a1))
      write(ilp,263)
263   format(/)
      if(full) then
       WRITE (ILP,265) KV
265    FORMAT(' Operating kV = ',F6.1)
       WRITE (ILP,270)CTH*1.0e4
270    FORMAT(' Effective carbon contamination (micr.) = ',F7.4)
       WRITE (ILP,275) WTH*1.0e4
275    FORMAT(' Effective Be window (micr.) = ',F4.1)
       write (ilp,276) auth*1.0e8
276    format(' Effective Au layer (Angs.) = ',f6.0)
       write(ilp,277) sith*1.0e8,siacth
277    format(' Effective dead layer (Angs.) = ',f6.0,/,
     & ' Effective active layer (cm) = ',f6.3)
      endif
      IF (clflag) THEN
        WRITE (ILP,*) 'Calculated C-L factors '
      ELSE
        WRITE (ILP,*) 'User supplied C-L factors :'
        DO 280 I=1,NE
         WRITE (ilp,278) (lab(j,na(1)),j=1,2),(lab(j,na(i)),j=1,2),
     &   clkr(i)
278      FORMAT (5x,'C-L factor (',2a1,',',2a1,'): ',f7.3)
280     CONTINUE
      ENDIF
      IF (abflag) THEN
       last=.false.
       WRITE (ILP,281)
281    format(//,20x,' Absorption correction requested.')
       if(itilt.eq.0) then
        write(ilp,282) th
282     format(20x,' Foil thickness (Angs.) : ',f6.0,/)
       else
        write(ilp,283) th
283     format(20x,' Apparent foil thickness (Angs.) : ',f6.0,/)
       endif
      ELSE
       last=.true.
       th=0.
       WRITE (ILP,*) 'NO absorption correction.'
      ENDIF
      IF (dnflag) THEN
       WRITE (ILP,*) 'No user supplied bulk density'
      ELSE
       WRITE (ILP,284) bulkro
284    format(' User supplied bulk density = ',f7.3,' g/cc')
      ENDIF
c
c     allow the computation of mu/ro coeff & the check of spurious densities
c     ----------------------------------------------------------------------
c
      muflag=.true.
      denchk=.true.
c
c     **************************************
c     Start the main loop for ND experiments
c     **************************************
c
      DO 2000 loop=1,ND
      write(*,400) loop
      if(ilp.eq.6) write(ilp,400) loop
400   format(//,20x,'***** Experiment number : ',i3,' *****')
c
c     ask for a title
c     ---------------
c
      WRITE (*,455)
455   FORMAT (/,' Enter a Title (1-40 char.): ',$)
      READ (*,456) (BTITLE(MM),MM=1,40)
456   FORMAT(40A1)
      WRITE(ILP,457) (BTITLE(MM),MM=1,40)
457   FORMAT(//,20X,'Title : ',40A1,/)
c
c     enable the detection limit test by QUANT & ask intensities & backgrounds
c     ------------------------------------------------------------------------
c
      limit=.true.
500   DO 510 I=1,NE
        call typlin(i)
502     WRITE (*,505) lintyp,(lab(j,NA(I)),j=1,2)
505     FORMAT (/,' Integrated ',a7,
     &  ' X-Ray peak intensity for element ',2a1,' : ',$)
        READ (*,46) XI(I)
c
c       check consistency
c
        if(xi(i).lt.1.) go to 502
c
c       resume
c
506     WRITE(*,507)
507     FORMAT(6x,' Integrated background intensity : ',$)
        READ(*,46) BACK(i)
c
c       check consistency
c
        if(back(i).lt.1.) go to 506
c
c       assume Poisson statistics
c
        hegert(i)=2.
510   continue
c
c     check X-ray intensity input
c     ---------------------------
c
      char='Y'
      WRITE (*,515) char
515   FORMAT(5X,/,' Are all the intensities OK ? (Y/N)<',a1,'> ',$)
      READ (*,21) char
      IF(char.EQ. 'N'.OR. char .EQ. 'n') THEN
       WRITE (*,108)
       GO TO 500
      ENDIF
c
c     output the intensities
c     ----------------------
c
      write (*,516)
516   format(/, ' Let me print your stuff...',/)
      WRITE(ILP,520)
520   FORMAT (/,' Element',2x,'Atomic no.',2X,'Line type',2x,
     &'Integrated Intensity',2x,'Integrated Background',/)
      do 530 i=1,ne
       call typlin(i)
530    WRITE(ILP,531) (lab(j,na(i)),j=1,2),NA(I)+10,lintyp,XI(i),
     & BACK(i)
531    FORMAT(3X,2a1,8x,i2,8x,a7,2(9x,F12.0))
c
c     compute the Cliff-Lorimer factor including the detector efficiency
c     ------------------------------------------------------------------
c
      if(clflag) CALL CLIFF
      if(.not.full) then
       DO 540 I=1,NE
540     WRITE (ilp,278) (lab(j,na(1)),j=1,2),(lab(j,na(i)),j=1,2),
     &  clkr(i)
      endif
c
c     compute the FACTOR f to be used by QUANT without absorption correction
c     ----------------------------------------------------------------------
c
      do 545 i=1,ne
545    factor(i)=1./clkr(i)
c
c     define if confidence levels have to be computed (last=.true.)
c     -------------------------------------------------------------
c
      if(abflag) then
       last=.false.
       WRITE(ILP,546)
546    FORMAT(/,16X,'***** Solution CORRECTED for absorption *****',/)
c
c      if new experiment, check the foil thickness
c
       if(loop.gt.1) then
        var=0.
        if(itilt.eq.0) then
         WRITE (*,547) th
547      FORMAT (' Foil thickness (Angs.) ? <',f6.0,'> :  ',$)
        else
         WRITE (*,548) th
548      FORMAT (' Apparent foil thickness (Angs.) ? <',f6.0,'> : ',$)
        endif
        read(*,46) var
        if(var.ne.0.) then
         th=var
         thc=th*1.0e-8
        endif
        if(itilt.eq.0) then
         write(ilp,282) th
        else
         write(ilp,283) th
        endif
       endif
      else
       last=.true.
       WRITE (ILP,549)
549    FORMAT(/,15x,'***** Solution UNCORRECTED for absorption *****',/)
      endif
c
c     first estimation of the concentrations
c     --------------------------------------
c
      call quant(na,ne,xi,back,hegert,factor,wconc)
c
c     if no absorption correction required, test the thin film criterion
c     ------------------------------------------------------------------
c
      IF (.not.abflag) then
       write(*,*) ' Checking the thin film criterion ...'
       CALL THICK
       char='Y'
       write(*,550) char
550    format(/,' Change your mind and need absorption correction ?',
     & ' (Y/N)<',a1,'> : ',$)
       read(*,21) char
       if((char.eq.'N').or.(char.eq.'n')) then
        abflag=.false.
        go to 812
       else
        abflag=.true.
        if(itilt.eq.0) then
         WRITE (*,113)
        else
         WRITE (*,114)
        endif
        READ (*,46) TH
        THC=TH * (1.0E-08)
        WRITE (ILP,281)
        if(itilt.eq.0) then
         write(ilp,282) th
        else
         write(ilp,283) th
        endif
        last=.false.
        WRITE(ILP,546)
       endif
      endif
c
c     Make the Absorption Correction.
c     -------------------------------
c
c     start of the iteration ITER for absorption correction
c
      DO 800 iter=1,10
       if(ilp.eq.6) then
        write(ilp,263)
        write(ilp,650) iter
       endif
       WRITE (*,650) iter
650    FORMAT(' Iteration #',I4)
c
c      compute mu/ro coefficients
c
       CALL ABCOR
c
c      compute specimen density with actual wght concentrations ?
c
       if(dnflag) CALL DENSE
c
c      compute the mean X-Ray path in the specimen (Cliff method)
c
       PATH = (TH*cos(tiltr))/2./SIN(PSIR)
       PATHC = PATH * 1.0E-08
c
c      calculate the correction coefficient COR & modify FACTOR for QUANT
c
       temp=EXP(-ABCOF(1)*BULKRO*PATHC)
       DO 660 I=2,NE
        cor(i)=(EXP(-ABCOF(i)*BULKRO*PATHC))/temp
        WRITE (ilp,655) (lab(j,na(1)),j=1,2),(lab(j,na(i)),j=1,2),
     &   cor(i)
655     FORMAT (5x,'Absorption Correction (',2a1,',',2a1,'): ',f7.3)
660    factor(i)=cor(i)/clkr(i)
c
c      compute new estimations WNEW of wght concentrations
c      ---------------------------------------------------
c
       if(abflag) last=.false.
       call quant(na,ne,xi,back,hegert,factor,wnew)
c
c      check convergence DELT of the iteration process
c      -----------------------------------------------
c
       delt=0.
       do 670 i=1,ne
        delt=delt+abs(wnew(i)-wconc(i))
670    wconc(i)=wnew(i)
       if(delt.lt.0.05) then
c
c       abort the iteration and output the results
c
        last=.true.
        call quant(na,ne,xi,back,hegert,factor,wconc)
        call atomic
        go to 813
       endif
c
c       otherwise, do another iteration
c
800   continue
c
c     Conclusion
c     ----------
c
810   if(iter.ge.10) then
       write(ilp,811)
811    format(' = Absorption correction aborted after 10 iterations =')
      endif
      last=.true.
      call quant(na,ne,xi,back,hegert,factor,wconc)
812   call atomic
813   if(th.gt.0.) then
       if(itilt.eq.0) then
        WRITE(ILP,815) TH
815     FORMAT(/,5X,'At this analysis point, the foil ',
     &  'thickness is',F6.0,' Angstroms')
       else
        write(ilp,816) th
816     FORMAT(/,5X,'At this analysis point, the apparent foil ',
     &  'thickness is',F6.0,' Angstroms')
       endif
       WRITE(ILP,817) PATH
817    FORMAT(14X,'The mean X-Ray path length is',f6.0,' Angstroms.')
      endif
      WRITE(ILP,820) BULKRO
820   FORMAT(19X,'The bulk density is',F7.3,' g/cc.')
      write(ilp,910)
910   format(12x,46('*'),//)
2000  CONTINUE
c
c     ********************
c     END OF THE MAIN LOOP
c     ********************
      if(ilp.eq.6) close(6)
      write(*,2010)
2010  format(32x,'That''s all folks !',//)
      STOP
      END
      BLOCK DATA INIT
c
c     LAST VERSION 9-DEC-1988 14:56:02
c
c     purpose
c     -------
c       To set up most of the physical constants used by the main program
c
c     comments
c     --------
c       All these parameters are taken from references quoted before each
c     DATA statement :
c       - Bearden & Burr : Atomic energy levels, US Atomic Energy Commission,
c                          NYO-2543-1, 1965
c       - N. Zaluzec : Program EDQNT, available from the author
c       - Gatan Inc : Inner Shell Loss Table, from Gatan Inc. 6678 Owens Drive,
c                     Pleasanton, CA 94566 USA.
c
c     authors :
c     ---------
c       It's a super-coproduction involving :
c
c     1) N. ZALUZEC, ARGONNE NATIONAL LABORATORY, MATERIALS SCIENCE &
c        TECHNOLOGY DIVISION, ARGONNE, IL 60439
c     2) Pierre TREBBIA, US 41 : "Microscopie Electronique Analytique
c        Quantitative", CNRS, Laboratoire de Physique des Solides, Bat. 510
c        Universite Paris-Sud, F91405 ORSAY Cedex, Phone : (33-1) 69 41 53 68
c
      parameter (nmax=20)
      REAL LAMDA,l1,l2,l3,m1,m2,m3,m4,m5
      DIMENSION LAMDA(10,82),xlin1(121),xlin2(90),xlin3(86),
     &xlin4(86),xlin5(86),xlin6(86),xlin7(86),xlin8(86),xlin9(93)
      equivalence(xlin1(1),lamda(1,1))
      equivalence(xlin2(1),lamda(2,13))
      equivalence(xlin3(1),lamda(2,22))
      equivalence(xlin4(1),lamda(8,30))
      equivalence(xlin5(1),lamda(4,39))
      equivalence(xlin6(1),lamda(10,47))
      equivalence(xlin7(1),lamda(6,56))
      equivalence(xlin8(1),lamda(2,65))
      equivalence(xlin9(1),lamda(8,73))
      common /lamda/lamda
      dimension relint(10,82),ri1(162),ri2(136),ri3(144),ri4(129),
     &ri5(144),ri6(105)
      equivalence(ri1(1),relint(1,1))
      equivalence(ri2(1),relint(3,17))
      equivalence(ri3(1),relint(9,30))
      equivalence(ri4(1),relint(3,45))
      equivalence(ri5(1),relint(2,58))
      equivalence(ri6(1),relint(6,72))
      COMMON /RELINT/RELINT
      DIMENSION ek(82),L1(82),L2(82),l3(82),m1(82),m2(82),m3(82),
     &m4(82),m5(82)
      common /energ/ek,l1,l2,l3,m1,m2,m3,m4,m5
      DIMENSION ATWT(82),ro(82)
      COMMON /WT/ATWT
      common /dens/ro
      INTEGER F12,F13,F23,W3
      DIMENSION NK3(33),w3(65),f12(65),f13(65),f23(65)
      common /nestor/nk3,w3,f12,f13,f23
      dimension wk1(7)
      common /wk/wk1
c
c DATA FOR ionization energies EK for K edges in keV for Z = 11,92
c                     (taken from Bearden & Burr & ZALUZEC)
c
      DATA EK /1.0721,1.3050,1.5596,1.8389,2.1490,2.4720,2.8224,3.2029,
     & 3.6074,4.0381,4.4928,4.9653,5.4651,5.9892,6.5390,7.1113,
     & 7.7089,8.3328,8.9787,9.6586,10.3671,11.1031,11.8667,12.6578,
     & 13.4737,14.3256,15.1997,16.1046,17.0384,17.9976,18.9856,19.9995,
     & 21.0440,22.1172,23.2199,24.3503,25.5140,26.7112,27.9399,29.2001,
     & 30.4912,31.8138,33.1694,34.5614,35.9846,37.4406,38.9246,40.4430,
     & 41.9906,43.5689,45.1840,46.8342,48.5190,50.2391,51.9957,53.7885,
     & 55.6177,57.4855,59.3896,61.3323,63.3138,65.3508,67.4164,69.5250,
     & 71.6764,73.8708,76.1110,78.3948,80.7249,83.1023,85.5304,88.0045,
     & 90.5259,93.1050,95.7299,98.404,101.137,103.9219,106.7553,
     & 109.6509,112.6014,115.6061 /
c
c DATA FOR ionization energies for L1 edges in keV for Z = 11,92
c                     (taken from Bearden & Burr)
c
      DATA L1 / 0.0641,0.0892,0.1174,0.152,0.190,0.232,0.277,0.326,
     & 0.379,0.434,0.493,0.556,0.6273,0.6964,0.7700,0.8454,0.9259,
     & 1.0090,1.0951,1.1936,1.2977,1.4143,1.5265,1.6539,1.7820,1.9210,
     & 2.0651,2.2163,2.3725,2.5316,2.6977,2.8655,3.0425,3.2240,3.4119,
     & 3.6043,3.8058,4.0180,4.2375,4.4647,4.6983,4.9392,5.1881,5.4528,
     & 5.7143,5.9888,6.2663,6.5488,6.8348,7.1260,7.4279,7.7368,8.0520,
     & 8.3756,8.7080,9.0458,9.3942,9.7513,10.1157,10.4864,10.8704,
     & 11.2707,11.6815,12.0998,12.5267,12.9680,13.4185,13.8799,14.3528,
     & 14.8393,15.3467,15.8608,16.3875,16.9393,17.493,18.049,18.639,
     & 19.2367,19.840,20.4721,21.1046,21.7574 /
c
c DATA FOR ionization energies for L2 edges in keV for Z = 11,92
c                     (taken from Bearden & Burr)
c
      DATA L2 / 0.0311,0.0514,0.0733,0.0995,0.1362,0.1654,
     & 0.2016,0.2473,0.2963,0.350,0.4067,0.4604,0.5205,0.5837,
     & 0.6514,0.7205,0.7936,0.8719,0.9509,1.0428,1.1423,1.2478,1.3586,
     & 1.4762,1.5960,1.7272,1.8639,2.0068,2.1555,2.3067,2.4647,2.6251,
     & 2.7932,2.9669,3.1461,3.3303,3.5237,3.7270,3.9380,4.1561,4.3804,
     & 4.6120,4.8521,5.1037,5.3594,5.6236,5.8906,6.1642,6.4404,6.7215,
     & 7.0128,7.3118,7.6171,7.9303,8.2516,8.5806,8.9178,9.2643,9.6169,
     & 9.9782,10.3486,10.7394,11.1361,11.5440,11.9587,12.3850,12.8241,
     & 13.2726,13.7336,14.2087,14.6979,15.2000,15.7111,16.2443,16.7847,
     & 17.3371,17.9065,18.4843,19.0832,19.6932,20.3137,20.9476 /
c
c DATA FOR ionization energies for L3 edges in keV for Z = 11,92
c                     (taken from Bearden & Burr)
c
      DATA L3 / 0.0311,0.0514,0.0729,0.0989,0.1353,0.1642,
     & 0.200,0.2452,0.2936,0.3464,0.4022,0.4545,0.5129,0.5745,
     & 0.6403,0.7075,0.7786,0.8547,0.9309,1.0197,1.1154,1.2167,1.3231,
     & 1.4358,1.5499,1.6749,1.8044,1.9396,2.0800,2.2223,2.3705,2.5202,
     & 2.6769,2.8379,3.0038,3.1733,3.3511,3.5375,3.7301,3.9288,4.1322,
     & 4.3414,4.5571,4.7822,5.0119,5.2470,5.4827,5.7234,5.9643,6.2079,
     & 6.4593,6.7162,6.9769,7.2428,7.5140,7.7901,8.0711,8.3579,8.6480,
     & 8.9436,9.2441,9.5607,9.8811,10.2068,10.5353,10.8709,11.2152,
     & 11.5637,11.9187,12.2839,12.6575,13.0352,13.4186,13.8138,14.2135,
     & 14.6194,15.0312,15.4444,15.8710,16.3003,16.7331,17.1663 /
c
c DATA FOR ionization energies for M1 edges in keV for Z = 11,92
c            (taken from Bearden & Burr & Gatan Inc. Data)
c
      DATA M1 / 0.0007,0.0021,0.0007,0.0076,0.0162,0.0158,0.0175,
     & 0.0253,0.0339,0.0437,0.0538,0.0592,0.0665,0.0741,0.0839,0.0923,
     & 0.1007,0.1118,0.1197,0.1359,0.1581,0.1800,0.2035,0.2315,0.2565,
     & 0.289,0.3221,0.3575,0.3946,0.4303,0.4684,0.5046,0.544,0.5850,
     & 0.6271,0.6699,0.7175,0.7702,0.8256,0.8838,0.9437,1.0060,1.0721,
     & 1.145,1.2171,1.2928,1.3613,1.4346,1.5110,1.5753,1.649,1.7228,
     & 1.8000,1.8808,1.9675,2.0468,2.1283,2.2065,2.3068,2.3981,2.4912,
     & 2.6009,2.7080,2.8196,2.9317,3.0485,3.1737,3.2960,3.4249,3.5616,
     & 3.7041,3.8507,3.9991,4.1494,4.317,4.482,4.652,4.8220,5.002,
     & 5.1823,5.3669,5.5480 /
c
c DATA FOR ionization energies for M2 edges in keV for Z = 11,92
c            (taken from Bearden & Burr)
c
      DATA M2 /0.,0.,0.,0.003,0.0099,0.008,0.0068,0.0124,0.0178,0.0254,
     & 0.0323,0.0335,0.0378,0.0425,0.0486,0.0533,0.0595,0.0681,0.0734,
     & 0.0866,0.1068,0.1279,0.1464,0.1682,0.1893,0.2227,0.2474,0.2798,
     & 0.3124,0.3442,0.3784,0.4097,0.4449,0.4828,0.521,0.5591,0.6024,
     & 0.6507,0.7022,0.7564,0.8119,0.8697,0.9305,0.999,1.065,1.1367,
     & 1.2044,1.2728,1.3374,1.4028,1.4714,1.5407,1.6139,1.6883,1.7677,
     & 1.8418,1.9228,2.0058,2.0898,2.173,2.2635,2.3654,2.4687,2.5749,
     & 2.6816,2.7922,2.9087,3.0265,3.1478,3.2785,3.4157,3.5542,3.6963,
     & 3.8541,4.008,4.159,4.327,4.4895,4.656,4.8304,5.0009,5.1822/
c
c DATA FOR ionization energies for M3 edges in keV for Z = 11,92
c            (taken from Bearden & Burr)
c
      DATA M3 /0.,0.,0.,0.003,0.0099,0.008,0.0068,0.0124,0.0178,0.0254,
     & 0.0323,0.0335,0.0378,0.0425,0.0486,0.0533,0.0595,0.0681,0.0734,
     & 0.0866,0.1029,0.1208,0.1405,0.1619,0.1815,0.2138,0.2385,0.2691,
     & 0.3003,0.3305,0.363,0.3923,0.425,0.4606,0.4962,0.5315,0.5714,
     & 0.6165,0.6643,0.7144,0.7656,0.8187,0.8746,0.937,0.9976,1.0622,
     & 1.1234,1.1854,1.2422,1.2974,1.3569,1.4198,1.4806,1.544,1.6113,
     & 1.6756,1.7412,1.8118,1.8845,1.9498,2.0236,2.1076,2.194,2.281,
     & 2.3673,2.4572,2.5507,2.6454,2.743,2.8471,2.9566,3.0664,3.1769,
     & 3.3019,3.426,3.538,3.663,3.7918,3.909,4.0461,4.1738,4.3034/
c
c DATA FOR ionization energies for M4 edges in keV for Z = 11,92
c            (taken from Bearden & Burr)
c
      DATA M4 /0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.0066,0.0026,0.0022,
     & 0.0023,0.0033,0.0029,0.0029,0.0036,0.0015,0.0081,0.0174,0.0287,
     & 0.0412,0.0567,0.0701,0.0889,0.1118,0.135,0.1596,0.1824,0.2074,
     & 0.2303,0.2564,0.2836,0.3117,0.34,0.3728,0.4105,0.4508,0.4933,
     & 0.5369,0.5825,0.6313,0.6723,0.7395,0.7961,0.8485,0.9013,0.9511,
     & 0.9999,1.0515,1.106,1.1606,1.2172,1.275,1.3325,1.3915,1.4533,
     & 1.5146,1.5763,1.6394,1.7164,1.7932,1.8716,1.9489,2.0308,2.1161,
     & 2.2019,2.2911,2.3849,2.4851,2.5856,2.6876,2.798,2.9087,3.0215,
     & 3.1362,3.2484,3.3702,3.4908,3.6112,3.7276/
c
c DATA FOR ionization energies for M5 edges in keV for Z = 11,92
c            (taken from Bearden & Burr)
c
      DATA M5 /0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.0066,0.0026,0.0022,
     & 0.0023,0.0033,0.0029,0.0029,0.0036,0.0015,0.0081,0.0174,0.0287,
     & 0.0412,0.0567,0.069,0.0889,0.1103,0.1331,0.1574,0.18,0.2046,
     & 0.227,0.2529,0.2794,0.307,0.3347,0.3667,0.4037,0.4431,0.4848,
     & 0.5275,0.5721,0.6194,0.6723,0.7255,0.7807,0.8317,0.8833,0.931,
     & 0.9777,1.0269,1.0802,1.1309,1.1852,1.2412,1.2949,1.3514,1.4093,
     & 1.4677,1.5278,1.5885,1.6617,1.7351,1.8092,1.8829,1.9601,2.0404,
     & 2.1216,2.2057,2.2949,2.3893,2.484,2.5796,2.683,2.7867,2.8924,
     & 2.9999,3.1049,3.219,3.332,3.4418,3.5517/
c
c DATA FOR average number of vacancies NK3 to be used for L cross sections
c      including Coster-Kronig transitions for even Z = 28 to 92
c           (taken from ZALUZEC)
c
      DATA NK3/6600,6560,6480,6360,6270,6160,6050,5960,5860,5760,
     &5700,5630,5554,5480,5423,5361,5309,5266,5217,5181,5133,5104,5063,
     &5030,5001,4975,4943,4918,4893,4868,4844,4810,4792/
c
c DATA FOR Coster-Kronig transitions probabilities F12,F13,F23 for Z = 28 to 92
c      F12 & F13 must be divided by 100, F23 by 1000
c                 (taken from ZALUZEC)
c
      DATA F12/30,30,29,29,28,28,28,28,27,27,27,26,26,10,10,10,10,10,
     &10,10,10,10,17,17,18,18,19,19,19,19,19,19,19,19,19,19,19,19,19,19,
     &19,19,19,19,18,18,17,16,16,15,14,14,13,13,12,11,11,10,10,10,9,9,9,
     &8,8/
      DATA F13/55,54,54,53,53,53,52,52,52,52,52,52,52,61,61,61,61,60,
     &60,59,59,59,27,28,28,28,28,28,28,29,29,29,30,30,30,30,30,30,30,30,
     &30,29,29,28,28,28,28,33,39,45,50,53,56,57,58,58,58,59,58,58,58,58,
     &57,58,57/
      DATA F23/28,28,26,32,50,63,76,88,100,109,117,126,132,137,141,
     &144,148,150,151,153,155,157,157,156,155,154,154,154,153,153,153,
     &153,152,151,150,149,147,145,143,142,140,139,138,136,135,134,133,
     &130,128,126,124,122,120,118,116,113,111,111,110,109,108,108,108,
     &139,167/
c
c DATA FOR fluo yields WK1 for K lines and Z = 12 to 18
c              (taken from ZALUZEC)
c
      data wk1/0.02475,0.0380,0.043,0.058,0.0795,0.1,0.1225/
c
c DATA FOR L3 fluorescent yield W3 for Z = 28 to 92
c    W3 MUST BE DIVIDED BY 1000. EXCEPT WHEN Z<30 (DIVIDE BY 10,000)
c            (taken from ZALUZEC)
c
      DATA W3/93,110,12,13,15,16,18,20,22,24,26,28,31,34,37,40,43,
     &46,49,52,56,60,64,69,74,79,85,91,97,104,111,118,125,132,139,147,
     &155,164,174,182,192,201,210,220,231,243,255,268,281,294,306,320,
     &333,347,360,373,386,399,411,424,437,450,463,476,489/
c
c DATA FOR atomic density RO for Z = 11 to 92
c             taken from Elsevier's periodic table 1987
c                 (Refs : Weast R.C. Handbook 1984)
c
      DATA RO /0.97,1.74,2.702,2.33,1.82,2.07,1.56,1.40,0.86,1.54,2.989,
     & 4.54,5.96,7.20,7.43,7.86,8.90,8.90,8.92,7.14,5.904,5.35,5.727,
     & 4.81,3.119,2.6,1.532,2.6,4.469,6.49,8.57,10.2,11.5,12.3,12.4,
     & 12.02,10.5,8.642,7.30,7.28,6.684,6.25,4.93,3.06,1.879,3.51,6.145,
     & 6.657,6.773,7.004,7.22,7.52,5.243,7.90,8.229,8.55,8.795,9.006,
     & 9.321,6.965,9.84,13.31,16.6,19.35,20.53,22.48,22.421,21.45,18.88,
     & 13.594,11.85,11.344,9.8,9.4,0.0,9.73,2.41,5.0,10.07,11.7,15.37,
     & 19.05 /
c
c DATA FOR atomic weight ATWT for Z = 11 to 92
c             taken from Elsevier's periodic table 1987
c    (refs : Holden & Martin 1984, IUPAC 1985 & IUPAC 1986.)
c
      DATA ATWT /22.9898,24.305,26.9815,28.0855,30.9738,
     & 32.066,35.4527,39.948,39.0983,40.078,44.956,47.88,50.9415,51.996,
     & 54.938,55.847,58.9332,58.69,63.546,65.39,69.723,72.61,74.9216,
     & 78.96,79.904,83.80,85.4678,87.62,88.906,91.224,92.906,95.94,98.,
     & 101.07,102.9055,106.42,107.868,112.411,114.82,118.710,121.75,
     & 127.60,126.9045,131.29,132.9054,137.327,138.9055,140.115,140.908,
     & 144.24,145.,150.36,151.965,157.25,158.925,162.50,164.93,167.26,
     & 168.934,173.04,174.967,178.49,180.948,183.85,186.21,190.2,192.22,
     & 195.08,196.967,200.59,204.3833,207.2,208.98,209.,210.,222.,223.,
     & 226.,227.,232.038,231.036,238.03 /
c
c DATA FOR X-ray wavelengths LAMDA(LINE,Z) (Angs.) for Z=11 to 92
c      and lines Kalpha1, Kalpha2, Kbeta1, Kbeta2, Kbeta3, Lalpha1, Lalpha2,
c       Lbeta1, Lbeta2, Lbeta3.
c                 (source : N. ZALUZEC)
c
      DATA XLIN1/11.90843,11.90843,11.61686,0.,0.,0.,0.,0.,0.,0.,
     &9.88846,9.88846,9.56934,0.,0.,0.,0.,0.,0.,0.,8.33780,8.34021,
     &7.98291,0.,7.98291,0.,0.,0.,0.,0.,7.12505,7.12905,6.77827,0.,
     &6.77827,0.,0.,0.,0.,0.,6.15870,6.16180,5.80412,0.,5.80412,0.,0.,
     &0.,0.,0.,5.37384,5.37624,5.03146,0.,5.03146,0.,0.,0.,0.,0.,
     &4.72827,4.73187,4.40409,0.,4.40409,0.,0.,0.,0.,0.,4.19260,4.19540,
     &3.88641,0.,3.88641,0.,0.,0.,0.,0.,3.74212,3.74552,3.45433,0.,
     &3.45433,0.,0.,0.,0.,0.,3.35884,3.36254,3.09015,0.,3.09015,
     &36.39161,36.39161,36.02087,0.,0.,3.03115,3.03485,2.77967,0.,
     &2.77967,31.38608,31.38608,30.99380,0.,27.06879,2.74887,2.75257,
     &2.51418,0.,2.51418,27.42807,27.42807,27.06879,0.,23.88724,2.50408/
      DATA XLIN2/2.50758,2.28479,0.,2.28479,24.26122,24.26122,23.88724,
     &0.,21.19228,2.28989,2.29369,2.08500,0.,2.08500,21.63615,21.63615,
     &21.26497,0.,18.95638,2.10200,2.10590,1.91051,0.,1.91051,19.46236,
     &19.46236,19.10247,0.,17.19487,1.93621,1.94011,1.75681,0.,1.75681,
     &17.58515,17.58515,17.26666,0.,15.65344,1.78921,1.79311,1.62102,0.,
     &1.62102,15.97613,15.97613,15.67324,0.,14.31581,1.65812,1.66182,
     &1.50033,1.48853,1.50033,14.56819,14.56819,14.26641,0.,13.17486,
     &1.54083,1.54463,1.39233,1.38092,1.39233,13.33065,13.33065,
     &13.04997,0.,12.11881,1.43543,1.43923,1.29544,1.28363,1.29544,
     &12.25051,12.25051,11.98982,0.,11.19916,1.34024,1.34423,1.20794,
     &1.19614,1.20844,11.29095,11.29095,11.01997,0.,10.35720,1.25414/
      DATA XLIN3/1.25814,1.12915,1.11695,1.12955,10.43559,10.43559,
     &10.17861,0.,9.58074,1.17604,1.18004,1.05745,1.04515,1.05795,
     &9.67043,9.67043,9.41344,0.,8.93197,1.10495,1.10905,0.99225,
     &0.98005,0.99285,8.99026,8.99026,8.73678,0.,8.32050,1.03985,
     &1.04395,0.93295,0.92056,0.93345,8.37669,8.37669,8.12421,0.,
     &7.76782,0.98015,0.98425,0.87866,0.86626,0.87916,7.81682,7.81682,
     &7.57793,0.,7.26705,0.92566,0.92985,0.82876,0.81641,0.82936,
     &7.31845,7.32714,7.07616,0.,6.78947,0.87536,0.87956,0.78284,
     &0.77072,0.78334,6.86467,6.87227,6.62618,0.,6.36749,0.82896,
     &0.83316,0.74066,0.72873,0.74116,6.45029,6.45699,6.21430,5.97471,
     &5.98331,0.78585,0.79005,0.70168,0.68986,0.70218,6.07131,6.07721/
      DATA XLIN4/5.83692,5.58693,5.63263,0.74607,0.75038,0.66571,
     &0.65409,0.66631,5.72372,5.73162,5.49293,5.23765,5.31174,0.70920,
     &0.71350,0.63224,0.62092,0.63274,5.40664,5.41614,5.17855,4.92356,
     &5.01316,0.67497,0.67927,0.60127,0.59017,0.60187,5.11445,5.12505,
     &4.88856,4.63628,4.74087,0.64297,0.64737,0.57243,0.56160,0.57303,
     &4.84657,4.85417,4.62078,4.37299,4.48698,0.61321,0.61762,0.54557,
     &0.53505,0.54617,4.59848,4.60528,4.37459,4.13110,4.25299,0.58535,
     &0.58976,0.52052,0.51020,0.52113,4.36839,4.37609,4.14630,3.90961,
     &4.03560,0.55930,0.56371,0.49698,0.48696,0.49758,4.15470,4.16300,
     &3.93571,3.70402,3.83351,0.53495,0.53936,0.47503,0.46532,0.47564,
     &3.95711,3.96591,3.73872,3.51403,3.64522,0.51211,0.51652,0.45449/
      DATA XLIN5/0.44498,0.45510,3.77282,3.78092,3.55533,3.33894,
     &3.47073,0.49057,0.49497,0.43516,0.42584,0.43586,3.60083,3.60913,
     &3.38544,3.17555,3.30604,0.47033,0.47473,0.41702,0.40790,0.41772,
     &3.43993,3.44853,3.22604,3.02375,3.15295,0.45129,0.45570,0.39999,
     &0.39107,0.40059,3.28934,3.29894,3.07705,2.88246,3.00905,0.43325,
     &0.43776,0.38386,0.37544,0.38456,3.14895,3.15865,2.93776,2.75077,
     &2.87446,0.41592,0.42043,0.36842,0.35991,0.36938,3.01715,3.02675,
     &2.80676,2.62767,2.74767,0.40029,0.40480,0.35430,0.34608,0.35500,
     &2.89256,2.90206,2.68407,2.51218,2.62877,0.38506,0.38967,0.34077,
     &0.33275,0.34147,2.77657,2.78596,2.56838,2.40448,2.51678,0.37073,
     &0.37524,0.32795,0.32013,0.32865,2.66617,2.67587,2.45938,2.30309/
      DATA XLIN6/2.41058,0.35710,0.36161,0.31582,0.30821,0.31652,
     &2.56198,2.57108,2.35649,2.20909,2.31129,0.34408,0.34869,0.30420,
     &0.29688,0.30490,2.46328,2.47308,2.25899,2.11960,2.21739,0.33185,
     &0.35650,0.29328,0.28626,0.29398,2.37089,2.38088,2.16699,2.03640,
     &2.12720,0.32073,0.32484,0.28289,0.27599,0.28359,2.28229,2.29289,
     &2.08010,1.95611,2.04240,0.30891,0.31362,0.27304,0.26632,0.27374,
     &2.20009,2.21109,1.99830,1.88241,1.96260,0.29849,0.30320,0.26362,
     &0.25691,0.26442,2.12100,2.13160,1.92061,1.81201,1.88701,0.28837,
     &0.29318,0.25440,0.24809,0.25520,2.04710,2.05800,1.84711,1.74562,
     &1.81511,0.27875,0.28346,0.24598,0.23957,0.24679,1.97660,1.98770,
     &1.77691,1.68332,1.74732,0.26953,0.27434,0.23757,0.23176,0.23837/
      DATA XLIN7/1.90911,1.92001,1.71092,1.62402,1.68242,0.26081,
     &0.26552,0.23009,0.22429,0.23079,1.84511,1.85621,1.64772,1.56732,
     &1.62062,0.25250,0.25711,0.22264,0.21713,0.22344,1.78461,1.79571,
     &1.58762,1.51413,1.56182,0.24438,0.24909,0.21532,0.21009,0.21603,
     &1.72692,1.73832,1.53053,1.46423,1.50653,0.23677,0.24148,0.20871,
     &0.20360,0.20961,1.67222,1.68312,1.47593,1.41573,1.45253,0.22925,
     &0.23406,0.20210,0.19689,0.20290,1.61972,1.63042,1.42373,1.37043,
     &1.40163,0.22214,0.22695,0.19558,0.19078,0.19619,1.56972,1.58072,
     &1.37433,1.32664,1.35313,0.21542,0.22023,0.19007,0.18506,0.19088,
     &1.52213,1.53323,1.32704,1.28484,1.30694,0.20901,0.21382,0.18436,
     &0.17945,0.18516,1.47663,1.48763,1.28194,1.24494,1.26284,0.20280/
      DATA XLIN8/0.20761,0.17885,0.17414,0.17965,1.43303,1.44413,
     &1.23874,1.20684,1.22044,0.19679,0.20160,0.17364,0.16903,0.17444,
     &1.39143,1.40243,1.19734,1.17014,1.17974,0.19098,0.19589,0.16853,
     &0.16402,0.16933,1.35133,1.36263,1.15804,1.13555,1.14094,0.18546,
     &0.19037,0.16362,0.15931,0.16452,1.31314,1.32454,1.12005,1.10225,
     &1.10365,0.18015,0.18506,0.15901,0.15470,0.15981,1.27654,1.28794,
     &1.08365,1.07045,1.06805,0.17509,0.17999,0.15449,0.15029,0.15489,
     &1.24134,1.25274,1.04875,1.04005,1.03375,0.17013,0.17504,0.15010,
     &0.14599,0.15100,1.20864,1.21894,1.01525,1.01055,1.00075,0.16533,
     &0.17023,0.14599,0.14198,0.14679,1.17514,1.18654,0.98295,0.98245,
     &0.96925,0.16082,0.16573,0.14198,0.13807,0.14278,1.14404,1.15554/
      DATA XLIN9/0.95187,0.95545,0.93865,0.15639,0.16129,0.13809,
     &0.13429,0.13889,1.11395,1.12565,0.92206,0.92966,0.90856,0.15209,
     &0.15699,0.13429,0.13059,0.13519,1.08515,1.09685,0.89356,0.90486,
     &0.88146,0.14799,0.15289,0.13069,0.12709,0.13159,1.05735,1.06915,
     &0.86616,0.88106,0.85446,0.14399,0.14899,0.12719,0.12369,0.12809,
     &1.03065,1.04245,0.83946,0.85826,0.82796,0.14009,0.14509,0.12379,
     &0.12039,0.12469,1.00485,1.01665,0.81370,0.83566,0.80268,0.13639,
     &0.14139,0.12059,0.11719,0.12139,0.98005,0.99185,0.78896,0.81396,
     &0.77816,0.13276,0.13777,0.11743,0.11412,0.11823,0.95605,0.96805,
     &0.76511,0.79346,0.75479,0.12929,0.13429,0.11429,0.11119,0.11519,
     &0.93295,0.94495,0.74216,0.77362,0.73214,0.12595,0.13096,0.11142,
     &0.10829,0.11232,0.91076,0.92266,0.71992,0.75458,0.71020/
c
c DATA FOR X-ray wavelength relative intensity RELINT(LINE,Z) for Z=11 to 92
c       and lines Kalpha1, Kalpha2, Kbeta1, Kbeta2, Kbeta3, Lalpha1, Lalpha2,
c       Lbeta1, Lbeta2, Lbeta3.
c                (source : N. ZALUZEC)
c
      DATA RI1/100.,50.,1.,0.,0.,0.,0.,0.,0.,0.,100.,50.,1.,0.,0.,0.,0.,
     &0.,0.,0.,100.,50.,1.32,0.,0.66,0.,0.,0.,0.,0.,100.,50.33,2.78,0.,
     &1.39,0.,0.,0.,0.,0.,100.,50.72,4.45,0.,2.25,0.,0.,0.,0.,0.,100.,
     &50.53,6.18,0.,3.11,0.,0.,0.,0.,0.,100.,50.54,8.06,0.,4.09,0.,0.,
     &0.,0.,0.,100.,50.42,10.17,0.,5.17,0.,0.,0.,0.,0.,100.,50.58,11.61,
     &0.,5.80,0.,0.,0.,0.,0.,100.,50.59,12.62,0.,6.31,100.,0.,50.,0.,0.,
     &100.,50.53,12.84,0.,6.42,100.,11.40,56.77,0.,24.12,100.,50.68,
     &12.99,0.,6.58,100.,11.37,56.70,0.,11.96,100.,50.99,13.24,0.,6.69,
     &100.,11.41,56.84,0.,7.56,100.,50.93,12.91,0.,6.51,100.,11.44,
     &57.27,0.,6.04,100.,50.97,13.40,0.,6.80,100.,11.38,56.49,0.,5.17,
     &100.,51.23,13.48,0.,6.84,100.,11.44,56.51,0.,4.79,100.,51.04/
      DATA RI2/13.54,0.,6.88,100.,11.46,56.68,0.,4.66,100.,51.18,13.58,
     &0.10,6.92,100.,11.45,51.86,0.,4.40,100.,51.40,13.41,0.10,6.84,
     &100.,11.47,51.28,0.,4.21,100.,51.43,13.71,0.10,7.01,100.,11.44,
     &51.40,0.,4.32,100.,51.34,14.06,0.16,7.17,100.,11.42,55.68,0.,4.61,
     &100.,51.50,14.36,0.47,7.33,100.,11.33,47.28,0.,4.43,100.,51.53,
     &14.59,0.82,7.50,100.,11.41,47.96,0.,4.77,100.,51.60,15.05,1.24,
     &7.69,100.,11.42,48.56,0.,4.70,100.,51.82,15.29,1.74,7.84,100.,
     &11.40,48.97,0.,4.63,100.,51.82,15.38,2.29,7.85,100.,11.40,48.79,
     &0.,4.65,100.,51.89,15.58,5.,7.99,100.,11.34,48.87,0.,4.71,100.,
     &52.09,15.,5.,15.,100.,11.34,47.47,0.,4.70,100.,52.31,15.,5.,15.,
     &100.,11.38,48.91,0.74,5.06,100.,50.,15.,5.,15.,100.,11.34,47.39/
      DATA RI3/1.77,5.22,100.,50.,15.,5.,15.,100.,11.36,40.65,3.30,6.43,
     &100.,50.,15.,5.,15.,100.,11.33,40.92,4.51,6.30,100.,50.,15.,5.,
     &15.,100.,11.35,41.19,5.84,6.44,100.,50.,15.,5.,15.,100.,11.34,
     &41.30,7.26,6.54,100.,50.,15.,5.,15.,100.,11.34,41.83,8.54,6.67,
     &100.,50.,15.,5.,15.,100.,11.34,42.97,10.35,6.82,100.,50.,15.,5.,
     &15.,100.,11.28,43.97,11.65,7.37,100.,50.,15.,5.,15.,100.,11.39,
     &44.63,12.88,7.72,100.,50.,15.,5.,15.,100.,11.32,45.20,13.71,8.78,
     &100.,50.,15.,5.,15.,100.,11.34,54.32,14.69,15.47,100.,50.,15.,5.,
     &15.,100.,11.35,53.45,15.56,15.10,100.,50.,15.,5.,15.,100.,11.28,
     &53.98,16.27,14.58,100.,50.,15.,5.,15.,100.,11.33,53.86,17.06,
     &14.64,100.,50.,15.,5.,15.,100.,11.31,52.81,17.70,14.12,100.,50./
      DATA RI4/15.,5.,15.,100.,11.34,53.73,19.59,13.99,100.,50.,15.,5.,
     &15.,100.,11.35,53.81,19.05,13.78,100.,50.,15.,5.,15.,100.,11.31,
     &53.29,19.58,13.41,100.,50.,15.,5.,15.,100.,11.35,53.75,19.40,
     &13.25,100.,50.,15.,5.,15.,100.,11.34,53.59,19.52,13.09,100.,50.,
     &15.,5.,15.,100.,11.32,53.34,19.57,12.87,100.,50.,15.,5.,15.,100.,
     &11.33,53.85,19.60,12.47,100.,50.,15.,5.,15.,100.,11.33,54.39,
     &19.77,12.67,100.,50.,15.,5.,15.,100.,11.35,54.88,19.85,12.65,100.,
     &50.,15.,5.,15.,100.,11.34,55.16,20.14,12.55,100.,50.,15.,5.,15.,
     &100.,11.34,55.13,19.93,12.40,100.,50.,15.,5.,15.,100.,11.34,55.50,
     &20.,12.53,100.,50.,15.,5.,15.,100.,11.34,56.32,19.95,12.55,100./
      DATA RI5/50.,15.,5.,15.,100.,11.34,56.58,20.05,12.58,100.,50.,15.,
     &5.,15.,100.,11.35,57.29,20.06,12.73,100.,50.,15.,5.,15.,100.,
     &11.34,57.78,20.17,12.79,100.,50.,15.,5.,15.,100.,11.35,58.72,
     &20.36,13.10,100.,50.,15.,5.,15.,100.,11.35,57.79,20.48,13.16,100.,
     &50.,15.,5.,15.,100.,11.36,57.81,20.76,13.33,100.,50.,15.,5.,15.,
     &100.,11.37,57.74,20.97,13.51,100.,50.,15.,5.,15.,100.,11.32,55.10,
     &21.22,12.22,100.,50.,15.,5.,15.,100.,11.37,54.01,21.46,10.24,100.,
     &50.,15.,5.,15.,100.,11.28,52.10,21.60,8.74,100.,50.,15.,5.,15.,
     &100.,11.39,51.10,21.83,7.83,100.,50.,15.,5.,15.,100.,11.39,50.19,
     &21.95,6.90,100.,50.,15.,5.,15.,100.,11.40,49.38,22.21,6.47,100.,
     &50.,15.,5.,15.,100.,11.38,48.89,22.29,6.07,100.,50.,15.,5.,15./
      DATA RI6/100.,11.31,47.95,22.44,6.05,100.,50.,15.,5.,15.,100.,
     &11.44,50.,22.78,6.07,100.,50.,15.,5.,15.,100.,11.36,48.17,22.89,
     &6.03,100.,50.,15.,5.,15.,100.,11.41,47.56,23.05,6.,100.,50.,15.,
     &5.,15.,100.,11.39,48.08,23.25,6.07,100.,50.,15.,5.,15.,100.,11.37,
     &47.91,23.37,5.97,100.,50.,15.,5.,15.,100.,11.35,50.,23.58,6.,100.,
     &50.,15.,5.,15.,100.,11.38,50.,23.60,6.,100.,50.,15.,5.,15.,100.,
     &11.41,50.,23.60,6.,100.,50.,15.,5.,15.,100.,11.43,50.,23.60,6.,
     &100.,50.,15.,5.,15.,100.,11.40,50.,23.60,6./
c
      end
      SUBROUTINE CONFIG
c
c     LAST VERSION 9-DEC-1988 15:59:30
c
c     purpose
c     -------
c       To keep in a special file CLEDX.CFG most of the experimental parameters
c    needed by the main program.
c
c     authors :
c     ---------
c       It's a super-coproduction involving :
c
c     1) H. O. COLIJN, OHIO STATE UNIVERSITY, CENTRAL ELECTRON OPTICS FACILITY,
c        COLUMBUS, OH 43210
c     2) Pierre TREBBIA, US 41 : "Microscopie Electronique Analytique
c        Quantitative", CNRS, Laboratoire de Physique des Solides, Bat. 510
c        Universite Paris-Sud, F91405 ORSAY Cedex, Phone : (33-1) 69 41 53 68
c
      PARAMETER (nmax=20)
      CHARACTER CHAR
      REAL KV
      COMMON /PARAM/ KV
      dimension na(nmax),li(nmax)
      common /para1/ne,na,li
      COMMON /WIND/ WTH,CTH,AUTH,SITH,siacth
      COMMON /HOC/ ALPHA, PHI,detang
c
c     set-up
c     ------
c
      do 7 i=1,5
7      write(*,5)
5      format(/////)
      WRITE(*,10)
10    format(//,9x,'******* AEM Thin Film EDX Quantitation *******',//)
      WRITE(*,15)
15    format(7x,'Default answers are given in angle brackets "<>"',//,
     &9x, 'If it is the first time you run this program,',/,
     &16x, 'answer "N" to the next question')
      CHAR = 'Y'
      WRITE (*,20) char
20    FORMAT(/,' Have you recorded your favorite parameters ? (Y/N)<',
     &a1,'>: ',$)
      READ(*,30) CHAR
30    FORMAT(A1)
      IF(CHAR.EQ.'N' .OR. CHAR.EQ.'n') THEN
c
c       need input
c       ----------
c
        WRITE (*,31)
31      FORMAT(//,10x,'Enter the parameters for your microscope.',/,
     &  10x,'under the same format as the default printed value',//)
c
        kv=100.
        var=0.
        WRITE (*,32) kv
32      FORMAT(' What primary beam voltage do you usually use ? (kV) ',
     &  '<',f4.0,'> : ',$)
        READ (*,132) var
132     format(f12.0)
        if(var.ne.0.) kv=var
c
        phi=0.
        var=0.
        WRITE (*,33) phi
33      FORMAT (/,' Your analytical microscope may have 2 distinct ',
     &  'tilt axes :',/,5X,'*) the X-tilt axis, normal to',/,5X,
     &  '*) the Y-tilt axis',/,' If you have a single tilt axis,',
     &  ' it is defined as the X-tilt axis.',//,
     &  ' The detector azimuthal angle is defined as the',
     &  ' angle between',/,' the projection, on a plane normal to ',
     &  'the incident beam,',/,' of the normal to the detector window',
     &  ' and the specimen X-tilt axis.',//,' What is your detector ',
     &  'azimuthal angle ? (deg.) <',f2.0,'> : ',$)
        READ (*,132) var
        if(var.ne.0.) phi=var
c
        alpha=0.
        var=0.
        WRITE (*,34) alpha
34      FORMAT(/,' The detector elevation angle is defined as the',
     &  ' angle between',/,' the mean detected X-Ray path and the',
     &  ' normal to the electron beam axis.',//,' What is your ',
     &  'detector elevation angle ? (deg.) <',f2.0,'> : ',$)
        READ (*,132) var
        if(var.ne.0.) alpha=var
c
        detang=0.
        var=0.
        WRITE (*,35) detang
35      FORMAT(/,' The detector inclination angle is defined as the',
     &  ' angle between',/,' the normal to the electron beam axis',
     &  ' and the normal to the detector window.',/,' Usually, it',
     &  ' should be equal to the elevation angle.',//,' What is your',
     &  ' detector inclination angle ? (deg.) <',f2.0,'> : ',$)
        READ (*,132) var
        if(var.ne.0.) detang=var
c
        cth=0.
        var=0.
        WRITE (*,36) cth
36      FORMAT(/,' What is your usual carbon contamination thickness ?',
     &  ' (micr.) <',f5.3,'> : ',$)
        READ (*,136) var
136     format(f12.3)
        if(var.ne.0.) cth=var
c
        wth=0.
        var=0.
        WRITE (*,37) wth
37      FORMAT(/,' If you don''t have a Be window, answer 0.',/,
     &  ' Otherwise, the Be window thickness is usually 8. micr.',//,
     &  ' What is your beryllium window thickness ? (micr.) ',
     &  '<',f5.3,'> : ',$)
        READ (*,136) var
        if(var.ne.0.) wth=var
c
        auth=0.
        var=0.
        WRITE(*,40) auth
40      FORMAT(/,' The detector gold contact layer thickness is ',
     &  'usually 200. Angs.',//,' What is your detector gold contact ',
     &  'layer thickness ? (Angs.) <',f2.0,'> : ',$)
        READ(*,132) var
        if(var.ne.0.) auth=var
c
        sith=0.
        var=0.
        WRITE(*,51) sith
51      FORMAT(/,' The silicon dead layer thickness is usually 1000.',
     &  ' Angs.',//,' What is your silicon dead layer thickness ? ',
     &  '(Angs.) <',f2.0,'> : ',$)
        READ(*,132) var
        if(var.ne.0.) sith=var
c
        siacth=0.
        var=0.
        WRITE(*,52) siacth
52      FORMAT(/,' The silicon active layer thickness is usually 0.3 ',
     &  'cm.',//,' What is your silicon active layer thickness ? (cm) ',
     &  '<',f5.3,'> : ',$)
        READ(*,136) var
        if(var.ne.0.) siacth=var
c
55      ne=nmax
        n=0
        WRITE (*,60) nmax,ne
60      FORMAT(/,' How many elements do you usually analyze for ',
     &  '(max : ',i2,') ? <',i2,'> : ',$)
        READ (*,61) N
61      format(i2)
        if(n.gt.nmax) go to 55
        if(n.ne.0) ne=n
c
c       save this input ?
c       -----------------
c
        CHAR = 'Y'
        WRITE(*,81) char
81      FORMAT(/,' Save these parameters to disk ? (Y/N)<',a1,'>',$)
        READ(*,30) CHAR
        IF(.NOT.(CHAR .EQ. 'N' .OR. CHAR .EQ.'n')) THEN
          OPEN (1,file='CLEDX.CFG',STATUS='NEW')
          WRITE (1,95) KV,ALPHA,PHI,WTH,AUTH,SITH,NE,CTH,siacth,detang
95        FORMAT (3(F5.0,/),F6.2,/,F5.0,/,F6.0,/,I2,/,3(F7.3,/))
          CLOSE (1)
        ENDIF
c
c     use previously stored data
c     --------------------------
c
      ELSE
        OPEN (1, file='CLEDX.CFG',STATUS='OLD')
        READ (1,*) KV,ALPHA,PHI,WTH,AUTH,SITH,NE,CTH,siacth,detang
        CLOSE (1)
      ENDIF
c
c     show set-up
c     -----------
c
      write(*,100) kv,phi,alpha,detang
100   format(//,' Here are your favorite parameters :',//,
     &  ' Primary beam voltage (kV) : ',f6.1,/,' Detector azimuthal ',
     &  'angle (deg.) : ',f4.0,/,' Detector elevation angle (deg.) : ',
     &  f4.0,/,' Detector inclination angle (deg.) : ',f4.0)
      write(*,110) cth,wth,auth,sith,siacth
110   format(' Contamination carbon layer (micr.) : ',f6.2,/
     &  ' Beryllium window (micr.) : ',f6.2,/,' Gold layer (Angs.) : ',
     &  f6.2,/,' Silicon dead layer (Angs.) : ',f8.2,/,' Silicon ',
     &  'active layer (cm) : ',f6.3)
      RETURN
      END
 
      SUBROUTINE CLIFF
c
c     LAST VERSION 9-DEC-1988 15:21:11
c
c     purpose
c     -------
c       To compute the so-called Cliff-Lorimer factors including the
c     coefficient of the detector absorption, for elements of atomic number 11
c     to 92.  K lines can be used for Z= 11 to 71, and L lines for Z= 28 to 92.
c
c     subroutines and function subprograms required
c     ---------------------------------------------
c     TYPLIN
c
c     comments
c     --------
c       The cross sections are calculated using the Zaluzec formulae (including
c     Coster-Kronig transition probabilities for the L shell). The fluorescence
c     yields are either tabulated or computed with the Bambinek formula,
c     the alpha weight of a line with the Schreiber formulae. The detection
c     efficiency takes into account the non-absorption probability of the
c     silicon active layer, the absorption by a possible carbon contamination
c     layer, all the specific thicknesses of the detector being corrected by
c     a possible misalignement. These absorption coefficients are calculated
c     for each original line in a series and weighted a-posteriori.
c
c     authors :
c     ---------
c       It's a super-coproduction involving :
c
c     1) A. D. ROMIG, JR., DIVISION 1832, SANDIA NATIONAL LABORATORIES,
c        ALBUQUERQUE, NEW MEXICO 87185
c     2) H. O. COLIJN, OHIO STATE UNIVERSITY, CENTRAL ELECTRON OPTICS FACILITY,
c        COLUMBUS, OH 43210
c     3) N. ZALUZEC, ARGONNE NATIONAL LABORATORY, MATERIALS SCIENCE &
c        TECHNOLOGY DIVISION, ARGONNE, IL 60439
c     4) Pierre TREBBIA, US 41 : "Microscopie Electronique Analytique
c        Quantitative", CNRS, Laboratoire de Physique des Solides, Bat. 510
c        Universite Paris-Sud, F91405 ORSAY Cedex, Phone : (33-1) 69 41 53 68
c
      parameter (nmax=20)
      DIMENSION na(nmax),li(nmax),alfw(nmax),CLKR(nmax),ATWT(82)
      common /para1/ne,na,li
      common /para2/alfw,clkr
      COMMON /WT/ ATWT
      character lab(2,82)
      common /label/lab
      common /ilp/ilp
      logical full,last,clflag,abflag,dnflag,muflag,limit,muprnt,denchk
      common /flags/full,last,clflag,abflag,dnflag,muflag,limit,muprnt,
     &              denchk
      REAL KV,lam,l1,l2,l3,m1,m2,m3,m4,m5
      DIMENSION lam(6,nmax),reli(5,nmax)
      COMMON /PARAM/ KV
      COMMON /WAVE/ lam,reli
      DIMENSION ek(82),L1(82),L2(82),l3(82),m1(82),m2(82),m3(82),m4(82),
     &m5(82)
      common /energ/ek,l1,l2,l3,m1,m2,m3,m4,m5
      COMMON /WIND/ WTH,CTH,AUTH,SITH,siacth
      INTEGER F12,F13,F23,W3
      DIMENSION NK3(33),w3(65),f12(65),f13(65),f23(65)
      common /nestor/nk3,w3,f12,f13,f23
      dimension wk1(7)
      common /wk/wk1
      dimension yfl(nmax),over(nmax),cross(nmax),detect(6,nmax),
     &enrj(nmax)
      dimension cmu(5,nmax),bemu(5,nmax),aumu(5,nmax),simu(5,nmax)
      character*7 lintyp
      common /lintyp/lintyp
c
c     set-up constants
c     ----------------
c
      clflag=.false.
      ROC=2.26
      ROBE=1.848
      ROAU = 19.32
      ROSI = 2.33
      const=6.5144
      betsq=1.-(1./(1.+kv/510.976))**2.
      beta=sqrt(betsq)
      relkv=255.488*betsq
c
c     beta = v/c
c     relkv = 0.5*Mo*v**2 (Mo = rest mass)
c
      do 50 i=1,ne
       z=float(na(i)+10)
c
c      the fluorescence yields : YFL
c      -----------------------------
c
       if(li(i).le.3) then
        xfluo=0.015+0.0327*z-0.64e-06*z*z*z
        temp=xfluo**4.
        yfl(i)= temp/(temp+1.)
        if((z.gt.11.).and.(z.lt.19.)) yfl(i)=wk1(na(i)-1)
       else
        if(na(i).lt.20) then
         div=10000.
        else
         div=1000.
        endif
        yfl(i)=float(w3(na(i)-17))/div
       endif
c
c      the weight of the line in the series : ALFW
c      -------------------------------------------
c
       if(alfw(i).eq.0.) then
        if(li(i).gt.3) then
c
c        L lines
c
         if(z.lt.51.) then
          alfw(i)=1.617-z*0.0398+z*z*3.766e-4
         else
          alfw(i)=0.609-z*1.619e-3-0.03248*sin(0.161*(z-51.))
         endif
c
c        beta lines
c
         if(li(i).eq.5) alfw(i)=1.-alfw(i)
        else
c
c       K lines
c
         if(z.lt.20.) then
          alfw(i)=1.052-z*z*4.39e-4
         else
          if(z.lt.30.) then
           alfw(i)=0.896-z*6.575e-4
          else
           if(z.lt.61.) then
            alfw(i)=1.0366-z*6.82e-3+z*z*4.815e-5
           endif
          if(z.gt.60.) alfw(i)=0.80
          endif
         endif
c
c        beta lines
c
         if(li(i).eq.2) alfw(i)=1.-alfw(i)
        endif
       endif
c
c      ionization energy ENERJ & overvoltage OVER
c      ------------------------------------------
c
       if(li(i).lt.4) then
c
c       take K shell
c
        enrj(i)=ek(na(i))
       else
c
c       take L3 for alpha and L2 for mean beta
c
        enrj(i)=l3(na(i))
        if(li(i).eq.5) enrj(i)=l2(na(i))
       endif
       over(i)=relkv/enrj(i)
c
c     the ionization cross-sections : CROSS
c     -------------------------------------
c
       if(li(i).le.3) then
c
c       K cross section
c
        BK=0.988-0.01883*z+3.066e-04*z*z-2.154e-06*z*z*z
        CK=0.2821+0.077*z-3.807e-03*z*z+8.262e-05*z*z*z
     &     -4.748e-07*z*z*z*z
        xfact=log(ck*relkv/enrj(i))-log(1.-betsq)-betsq
        CONK= const * 2. * BK
        cross(I)=CONK*xfact/enrj(i)/relkv
       else
c
c       L cross sections
c
c       1) compute vacancies VAC on L3 shell
c
        iz=na(i)+10
        irank=(iz-26)/2
        if(mod(iz,2).eq.0) then
c
c        the value is tabulated
c
         vac=float(nk3(irank))/10000.
        else
c
c        must interpolate
c
         ntemp=nk3(irank)+(nk3(irank+1)-nk3(irank))/2
         vac=float(ntemp)/10000.
        endif
c
c       2) compute Coster Kronig transitions probabilities CK1,CK2,CK3
c
        ck1=float(f12(na(i)-17))/100.
        ck2=float(f13(na(i)-17))/100.
        ck3=float(f23(na(i)-17))/1000.
c
c       3) compute individual cross sections for K, L1, L2, and L3 subshells
c
c           a) K shell SIGK
c
        ec=ek(na(i))
        BK=0.988-0.01883*z+3.066e-04*z*z-2.154e-06*z*z*z
        conk=const*2.*bk
        CK=0.2821+0.077*z-3.807e-03*z*z+8.262e-05*z*z*z
     &     -4.748e-07*z*z*z*z
        sigK=conk*(LOG(Ck*relkv/EC)-LOG(1-BETSQ)-BETSQ)/relkv/EC
c
c           b) L1 subshell SIGL1
c
        ec=l1(na(i))
        BL=0.05217 + 0.0037*Z -1.354E-4*Z*Z + 2.255E-6*Z*Z*Z
     &     - 1.742E-8*Z*Z*Z*Z + 5.072E-11*Z*Z*Z*Z*Z
        conl=const*8.*bl
        CL=15.9 - 0.1053*Z - 1.421E-3*Z*Z + 1.23E-5*Z*Z*Z
        sigL1=conl*(LOG(CL*relkv/EC)-LOG(1-BETSQ)-BETSQ)/relkv/EC
c
c           c) L2 subshell SIGL2
c
        ec=l2(na(i))
        BL=0.2776-6.609E-3*Z+8.711E-5*Z*Z - 3.599E-7*Z*Z*Z
        conl=const*8.*bl
        CL=-11.16+0.8708*Z-1.362E-2*Z*Z + 6.261E-5*Z*Z*Z
        sigl2=conl*(LOG(CL*relkv/EC)-LOG(1-BETSQ)-BETSQ)/relkv/EC
c
c           d) L3 subshell SIGL3
c
        ec=l3(na(i))
        BL=0.3962-9.284E-4*Z-1.555E-4*Z*Z+2.913E-6*Z*Z*Z
     &     -1.445E-8*Z*Z*Z*Z
        conl=const*8.*bl
        CL=5.167-0.4527*Z+0.02366*Z*Z-3.586E-4*Z*Z*Z+1.6700E-6*Z*Z*Z*Z
        sigl3=conl*(LOG(CL*relkv/EC)-LOG(1-BETSQ)-BETSQ)/relkv/EC
c
c       4) mix all for getting the L cross section CROSS
c
        cross(i)=sigl3+vac*sigK+(ck2+ck1*ck3)*sigl1+ck3*sigl2
       endif
c
c      mu/rho coefficients for the detector and for each line
c      ------------------------------------------------------
c
       var=0.
       detect(6,i)=0.
       do 40 k=1,5
        if((lam(k,i).ne.0.).and.(reli(k,i).ne.0.)) then
         xnrj=12.3981/lam(k,i)
c
c        use Tinh & Leroux parametrization if X-Ray energy > 1.486 keV
c        -------------------------------------------------------------
c
         if(xnrj.gt.1.486) then
c
c         carbon absorption
c
          IF(xnrj.GT.1.75) then
           Cmu(k,i)=3.8531*0.2838*(lam(k,i))**3.094
          else
           Cmu(k,i)=2.19*(lam(k,i))**2.7345
          endif
c
c         beryllium absorption
c
          BEmu(k,i)=2.4604*0.111*(lam(k,i))**3.03
c
c         gold absorption
c
          CAU=19.4943
          IF(xnrj.GT.14.3528) then
           aumu(k,i)=CAU*14.3528*(lam(k,i))**2.650
          else
           IF(xnrj.GT.13.7336) then
            aumu(k,i)=CAU*13.7336*(lam(k,i))**2.61439
           else
            IF(xnrj.GT.11.9187) then
             aumu(k,i)=CAU*11.9187*(lam(k,i))**2.3554
            else
             IF(xnrj.GT.3.4249) then
              aumu(k,i)=CAU*3.4249*(lam(k,i))**2.575
             else
              IF(xnrj.GT.3.1478) then
               aumu(k,i)=CAU*3.1478*(lam(k,i))**2.4471
              else
               IF(xnrj.GT.2.743) then
                aumu(k,i)=CAU*2.743*(lam(k,i))**2.4471
               else
                IF(xnrj.GT.2.2911) then
                 aumu(k,i)=CAU*2.2911*(lam(k,i))**2.4
                else
                 IF(xnrj.GT.2.2057) then
                  aumu(k,i)=CAU*2.2057*(lam(k,i))**2.2
                 else
                  aumu(k,i)=CAU*0.7588*(lam(k,i))**2.393
                 endif
                endif
               endif
              endif
             endif
            endif
           endif
          endif
c
c         silicon absorption
c
          CSI=9.1309
          IF(xnrj.GT.5.9) then
           SImu(k,i)=CSI*1.8389*(lam(k,i))**2.94
          else
           IF(xnrj.GT.1.8389) then
            SImu(k,i)=19.57*(lam(k,i))**2.7345
           else
            SImu(k,i)=CSI*0.1487*(lam(k,i))**2.79
           endif
          endif
c
         else
c
c         use Holton parametrization of Henke & Ebisu tables of mu/rho
c         ------------------------------------------------------------
c
c         carbon absorption
c
          IF(lam(k,i).GE.44.7) then
           cmu(k,i)=0.2572*(lam(k,i))**2.404
          else
           IF(lam(k,i).GE.17.6) then
            cmu(k,i)=4.129*(lam(k,i))**2.529
           else
            cmu(k,i)=1.966*(lam(k,i))**2.788
           endif
          endif
c
c         beryllium absorption
c
          IF(lam(k,i).GE.44.7) then
           bemu(k,i)=2.248*(lam(k,i))**2.419
          else
           IF(lam(k,i).GE.18.3) then
            bemu(k,i)=.506*(lam(k,i))**2.831
           else
            bemu(k,i)=.3102*(lam(k,i))**3.001
           endif
          endif
c
c         gold absorption
c
          IF(lam(k,i).GE.16.0) then
           aumu(k,i)=-4.411E-04*(lam(k,i))**5+9.878E-2*(lam(k,i))**4
     &               -8.218*(lam(k,i))**3+3.02E2*(lam(k,i))**2
     &               -4.466E3*lam(k,i)+2.966E4
          else
           aumu(k,i)=41.17*(lam(k,i))**1.906
          endif
c
c         silicon absorption
c
          IF(lam(k,i).GE.31.4) then
           simu(k,i)=26.2*(lam(k,i))**1.888
          else
           IF(lam(k,i).GE.18.3) then
            simu(k,i)=2.568*(lam(k,i))**2.574
           else
            simu(k,i)=1.759*(lam(k,i))**2.706
           endif
          endif
c
         endif
        endif
c
c       detection efficiency DETECT for each line
c       -----------------------------------------
c
        detect(k,i)=cmu(k,i)*roc*cth+bemu(k,i)*robe*wth+aumu(k,i)*roau
     &  *auth+simu(k,i)*rosi*sith
        detect(k,i)=exp(-detect(k,i))*(1.-exp(-simu(k,i)*rosi*siacth))
c
c       mean detection efficiency DETECT(6,i)
c       -------------------------------------
c
        detect(6,i)=detect(6,i)+detect(k,i)*reli(k,i)
40     var=var+reli(k,i)
       detect(6,i)=detect(6,i)/var
c
c      Cliff-Lorimer factors including detection efficiency CLKR
c      ---------------------------------------------------------
c
50    clkr(i)=atwt(na(i))/cross(i)/yfl(i)/alfw(i)/detect(6,i)
c
c     Normalize
c     ---------
c
      do 100 i=2,ne
       clkr(i)=clkr(i)/clkr(1)
c
c      Test if reasonnable values
c      --------------------------
c
       if((clkr(i).gt.100.).or.(clkr(i).lt.0.01)) then
        write(*,105) (lab(j,na(i)),j=1,2),clkr(i)
105     format(' Program aborted : ',2a1,' C-L factor : ',f12.4,' !!!')
        go to 1000
       endif
100   continue
      clkr(1)=1.
c
c     Print the results
c     -----------------
c
      if (full) then
       write(ilp,110) beta
110    format(/,15x,'*** Results of Cliff-Lorimer algorithm ***',//,
     & 25x,'beta = v/c = ',f6.4,/)
       do 150 i=1,ne
        xnrj=12.3981/lam(6,i)
        call typlin(i)
150     write(ilp,120) (lab(j,na(i)),j=1,2),na(i)+10,atwt(na(i)),lintyp,
     &  lam(6,i),xnrj,enrj(i),over(i),cross(i)*1.0e-20,yfl(i),alfw(i),
     &  detect(6,i),clkr(i)
120     format(/,29x,'Element : ',2a1,/,
     &  ' Atomic nb : ',i2,5x,' Atomic weight : ',f7.3,
     &  5X,' X-Ray line type : ',a7,/,
     &  ' Mean wavelength (Ang.) : ',f7.4,7x,' Mean X-Ray energy ',
     &  '(keV) : ',f6.3,/,' Ionization energy (keV) : ',
     &  f8.4,6x,'Overvoltage : ',f8.4,/,' Cross section (cm2) : ',
     &  e12.4,6x,'Fluo yield : ',f6.4,/,' Weight of the line : ',
     &  f6.4,12x,' Detection efficiency : ',f6.4,/,
     &  23x,'Overall C-L factor : ',f8.4,/)
      endif
      RETURN
1000  stop
      END
      SUBROUTINE TYPLIN(I)
c
c     VERSION 9-DEC-1988 15:52:18
c
c       purpose
c       -------
c       to code in a variable LINTYP the X-Ray line type
c
c       author :
c       --------
c       Pierre TREBBIA
c       US 41 : "Microscopie Electronique Analytique Quantitative"
c       Laboratoire de Physique des Solides, Bat. 510
c       Universite Paris-Sud, F91405 ORSAY Cedex
c       Phone : (33-1) 69 41 53 68
c
      parameter(nmax=20)
      dimension na(nmax),li(nmax)
      common/ para1/ne,na,li
      character*7 lintyp
      common /lintyp/lintyp
c
      lintyp='unknown'
      if(li(i).eq.1) lintyp='K alpha'
      if(li(i).eq.2) lintyp='K beta '
      if(li(i).eq.3) lintyp='K total'
      if(li(i).eq.4) lintyp='L alpha'
      if(li(i).eq.5) lintyp='L beta '
      if(li(i).eq.6) lintyp='L total'
      return
      end
      SUBROUTINE DENSE
c
c     LAST VERSION 9-DEC-1988 15:55:50
c
c     purpose
c     -------
c       This subroutine calculates the bulk density of an alloy.
c     The reciprocal bulk density is assumed to be equal to the sum
c     of the reciprocal densities of the components multiplied by the weight
c     concentrations
c
c     comments
c     --------
c       Some of elemental densities are not well known and their tabulated
c     values vary from author to author. These so-called "spurious" densities
c     are manually confirmed each time they are used.
c
c     author :
c     --------
c      Pierre TREBBIA, US 41 : "Microscopie Electronique Analytique
c      Quantitative", CNRS, Laboratoire de Physique des Solides, Bat. 510
c      Universite Paris-Sud, F91405 ORSAY Cedex, Phone : (33-1) 69 41 53 68
c
      parameter (nmax=20)
      DIMENSION wconc(nmax),aconc(nmax),ATWT(82),na(nmax),li(nmax)
      common /conc/wconc,aconc
      COMMON /WT/ATWT
      common /para1/ne,na,li
      dimension ro(82)
      common /dens/ro
      character char,lab(2,82)
      common /label/lab
      common /ilp/ilp
      COMMON /CONST/ BULKRO,PSIR
      logical full,last,clflag,abflag,dnflag,muflag,limit,muprnt,
     &        denchk,flag
      common /flags/full,last,clflag,abflag,dnflag,muflag,limit,muprnt,
     &              denchk
c
      TOT = 0
      DO 10 I=1,NE
       n=na(i)+10
c
c      if not already done, check spurious densities
c      ---------------------------------------------
c
       flag=.false.
       if(denchk) then
        if(n.eq.85) then
         flag=.true.
         write(*,1) (lab(j,na(i)),j=1,2)
1        format(/,'*** CAUTION : ',2a1,' DENSITY NOT DEFINED ***')
         write(*,2)
2        format(' Define the density value : ',$)
         read(*,*) ro(na(i))
        endif
        if((n.eq.17).or.(n.eq.18).or.(n.eq.36).or.(n.eq.43).or.(n.eq.54)
     &   .or.(n.eq.61).or.(n.eq.67).or.(n.eq.79).or.(n.eq.84).or.
     &   (n.eq.86).or.(n.eq.87).or.(n.eq.89)) then
         flag=.true.
         char='Y'
         write(*,3) (lab(j,na(i)),j=1,2),ro(na(i)),char
3        format(/,' Attention : ',2a1,' density : ',f7.3,
     &   ' O.K. ? (Y/N)<',a1,'> ',$)
         read(*,4) char
4        format(a1)
         if(char.eq.'N'.or.char.eq.'n') then
          write(*,5)
5         format(' Enter new value : ',$)
          read(*,*) ro(na(i))
         endif
        endif
        if(flag) write(ilp,9) (lab(j,na(i)),j=1,2),ro(na(i))
9       format(/,' Attention : ',2a1,' density redefined as : ',f7.3)
       endif
10     TOT = TOT + wconc(i)/ro(na(i))
      bulkro=1.0/tot
c
c     prevent another density check
c     -----------------------------
c
      denchk=.false.
      RETURN
      END
      SUBROUTINE QUANT(na,nd,p,b,h,f,c)
c
c     LAST VERSION 9-DEC-1988 15:34:55
c
c     purpose
c     -------
c      To estimate concentrations, variances and covariances from the following
c     input data : signal, background, h, f.
c
c               h = 1. + [variance(background)]/background
c               (Note : h = 2. if background obeys a Poisson distribution)
c               f = interaction-detection yield in arbitrary unit
c
c     subroutines and function subprograms required
c     ---------------------------------------------
c     invmat
c
c     comments
c     --------
c       This subroutine was built from the program QUANTI described in
c     "Unbiased Method...", P. TREBBIA, Ultramicroscopy, 24, 399-408 (1988).
c
c     author :
c     --------
c       Pierre TREBBIA
c       US 41 : "Microscopie Electronique Analytique Quantitative"
c       Laboratoire de Physique des Solides, Bat. 510
c       Universite Paris-Sud, F91405 ORSAY Cedex
c       Phone : (33-1) 69 41 53 68
c
      PARAMETER (NEL=20)
      logical full,last,clflag,abflag,dnflag,muflag,limit,muprnt,denchk
      common /flags/full,last,clflag,abflag,dnflag,muflag,limit,muprnt,
     &              denchk
      character lab(2,82)
      common /label/lab
      common /ilp/ilp
      double precision sigma,rtot,cftot,ptot,fact,test,try
      dimension sigma(nel-1,nel-1),nb(nel),p(1),b(1),h(1),na(1),
     &f(1),fn(nel),r(nel-1),vp(nel),c(nel),sc(nel),z(3),res(3,nel)
c
c     Put in z-array the Gaussian Zalpha values for alpha 0.10, 0.05, 0.01
c     --------------------------------------------------------------------
c
      z(1)=1.645
      z(2)=1.96
      z(3)=2.575
c
c     Compute peak variance
c     ---------------------
c
      if(last) write(ilp,10)
10    format(/,20x,'**** Concentration Results ****',/)
      norder=nel-1
      nfree=nd-1
      do 30 i=1,nd
       nb(i)=i
       vp(i)=p(i)+h(i)*b(i)
       if(limit.and.(p(i)*p(i)/vp(i).lt.13.)) then
        write(ilp,25) (lab(jj,na(nb(i))),jj=1,2)
25      format(' **** Caution : Check the detection limit of element ',
     &  2a1,' => Use program LIMIT !')
       endif
30    continue
c
c     disable the detection limit test
c
      limit=.false.
c
c     Interchange last element with element of lowest f*f/var(peak)
c     -------------------------------------------------------------
c
35    k=nd
      test=f(k)*f(k)/vp(k)
      do 40 i=1,nfree
       try=f(i)*f(i)/vp(i)
       if(try.lt.test) then
        k=i
        test=f(k)*f(k)/vp(k)
       endif
40    continue
      if(k.ne.nd) then
       if(last) write(ilp,45) (lab(jj,na(k)),jj=1,2),(lab(jj,na(nd)),
     &  jj=1,2)
45      format(/,6x,'*** Warning : elements ',2a1,' and ',2a1,
     &   ' have been interchanged ***',/)
       nb(k)=nd
       nb(nd)=k
       save=p(nd)
       p(nd)=p(k)
       p(k)=save
       save=b(nd)
       b(nd)=b(k)
       b(k)=save
       save=h(nd)
       h(nd)=h(k)
       h(k)=save
       save=f(nd)
       f(nd)=f(k)
       f(k)=save
       save=vp(nd)
       vp(nd)=vp(k)
       vp(k)=save
      endif
c
c     Compute the concentrations c(i)
c     -------------------------------
c
55    rtot=0.
      do 60 i=1,nfree
       r(i)=p(i)*f(nd)/(f(i)*p(nd))
60    rtot=rtot+r(i)
      do 65 i=1,nfree
65     c(i)=r(i)/(1.+rtot)
      c(nd)=1./(1.+rtot)
c
c     Compute the normalized FN(i)
c     ----------------------------
c
      if (last) then
       ptot=0.
       cftot=0.
       do 70 i=1,nd
        ptot=ptot+p(i)
70     cftot=cftot+c(i)*f(i)
       fact=ptot/cftot
       do 75 i=1,nd
75      fn(i)=f(i)*fact
c
c     Compute the sigma-1 matrix
c     --------------------------
c
       do 200 i=1,norder
        do 200 j=i,norder
         if((i.le.nfree).and.(j.le.nfree)) then
          if(i.eq.j) then
           sigma(i,i)=(fn(i)*fn(i)/vp(i))+(fn(nd)*fn(nd)/vp(nd))
          else
           sigma(i,j)=fn(nd)*fn(nd)/vp(nd)
          endif
         else
          if(i.eq.j) then
           sigma(i,i)=1.
          else
           sigma(i,j)=0.
          endif
         endif
         sigma(j,i)=sigma(i,j)
200    continue
c
c      Invert the sigma-1 matrix
c      -------------------------
c
       write(*,*) ' Coffee break ...'
       call invmat(sigma,norder)
       write(*,*) ' .... Wake up !!!'
c
c      Compute standard deviations
c      ---------------------------
c
       vcnd=0.
       do 320 i=1,nfree
        do 310 j=1,nfree
310      vcnd=vcnd+sigma(i,j)
320    sc(i)=sqrt(sigma(i,i))
       sc(nd)=sqrt(vcnd)
c
c      Compute confidence intervals at alpha levels 0.10, 0.05, 0.01
c      -------------------------------------------------------------
c
       do 350 i=1,nd
        do 350 j=1,3
350      res(j,i)=z(j)*sc(i)
c
c      Print the results
c      -----------------
c
       write(ilp,400)
400    format(/,15x,'Estimation of Weight Concentrations',
     & /,15x,'-----------------------------------',
     & //,' Elem',6x,'Peak',4x,'Background',4x,'h',7x,
     & 'Var(P)',6x,'Wght %',6x,'Sigma(Wght %)',/)
       do 410 i=1,nd
        write(ilp,405) (lab(jj,na(nb(i))),jj=1,2),p(i),b(i),h(i),
     &  vp(i),c(i)*100.,sc(i)*100.
405     format(2x,2a1,2f12.0,f7.1,f12.0,2(6x,f7.4))
410    continue
       if(full) then
        if(nfree.ge.2) then
         write(ilp,430)
430      format(/,3x,'Estimation of covariances',/,3x,
     &   '-------------------------')
         do 440 j=1,nfree-1
          do 440 i=j+1,nfree
440        write(ilp,441)(lab(jj,na(nb(j))),jj=1,2),
     &     (lab(jj,na(nb(i))),jj=1,2),sigma(j,i)
441       format(' Covariance(',2a1,',',2a1,') = ',e14.5)
        endif
       endif
       write(ilp,450)
450    format(/,3x,'Estimation of confidence intervals assuming NO',
     &  ' ERROR on C-L coeff.',
     &  /,3x,'----------------------------------------------',
     &  '--------------------',//,
     &  5x,'At the confidence level alpha = 0.10 :',/)
       do 465 i=1,nd
465     write(ilp,466) (lab(jj,na(nb(i))),jj=1,2),c(i)*100.,
     &  res(1,i)*100.
466     format(' Element  ',2a1,' Interval :',f7.3,'  + -',f7.3)
       if(full) then
        write(ilp,470)
470     format(/,5x,'At the confidence level alpha = 0.05 :',/)
        do 475 i=1,nd
475      write(ilp,476) (lab(jj,na(nb(i))),jj=1,2),c(i)*100.,
     &   res(2,i)*100.
476      format(' Element  ',2a1,' Interval :',f7.3,'  + -',f7.3)
        write(ilp,490)
490     format(/,5x,'At the confidence level alpha = 0.01 :',/)
        do 495 i=1,nd
495      write(ilp,496) (lab(jj,na(nb(i))),jj=1,2),c(i)*100.,
     &   res(3,i)*100.
496      format(' Element  ',2a1,' Interval :',f7.3,'  + -',f7.3)
       endif
      endif
c
c     Restore order if permutation occured
c     ------------------------------------
c
      if(k.ne.nd) then
       save=p(nd)
       p(nd)=p(k)
       p(k)=save
       save=b(nd)
       b(nd)=b(k)
       b(k)=save
       save=h(nd)
       h(nd)=h(k)
       h(k)=save
       save=f(nd)
       f(nd)=f(k)
       f(k)=save
       save=c(k)
       c(k)=c(nd)
       c(nd)=save
      endif
      return
      end
       SUBROUTINE INVMAT(ARRAY,NORDER)
c
c      LAST VERSION 9-DEC-1988 15:37:17
c
c       purpose
c       -------
c       To invert a symmetric matrix and calculates its determinant
c
c       description of parameters
c       -------------------------
c               input:
c               ------
c       array ---> input matrix which is replaced by its inverse
c       norder --> degree of matrix (order of determinant)
c
c               output:
c               -------
c       array ---> inverse of input matrix
c
c       comments
c       --------
c       1) Dimension statement valid for norder up to nmax = 9
c       2) Modify the statement PARAMETER to increase nmax.
c       3) Matrix ARRAY is in double precision floating point format.
c       4) This program is copied from " DATA REDUCTION AND ERROR ANALYSIS "
c          Author : P. R. BEVINGTON, McGraw-Hill Ed., 1969.
c
c       AN OVERFLOW CAN OCCUR ON VARIABLE DET IF NMAX IS SET TOO LARGE AND IF
c       THE COMPILER DOES NOT INCREASE THE ROOM TO STORE THE EXPONENT PART OF
c       THIS VARIABLE DET.
c
        PARAMETER (NMAX=19)
        double precision array,amax,save,det
        dimension array(nmax,nmax),ik(nmax),jk(nmax)
10      det=1.
11      do 100 k=1,norder
c
c        Find largest element array(i,j) in rest of matrix
c        -------------------------------------------------
c
         amax=0.
21       do 30 i=k,norder
          do 30 j=k,norder
23         if(dabs(amax)-dabs(array(i,j))) 24,24,30
24         amax=array(i,j)
           ik(k)=i
           jk(k)=j
30       continue
c
c        Interchange rows and columns to put amax in array(k,k)
c        ------------------------------------------------------
c
31       if(amax) 41,32,41
32       det=0.
         go to 140
41       i=ik(k)
         if(i-k) 21,51,43
43       do 50 j=1,norder
          save=array(k,j)
          array(k,j)=array(i,j)
50       array(i,j)=-save
51       j=jk(k)
         if(j-k) 21,61,53
53       do 60 i=1,norder
          save=array(i,k)
          array(i,k)=array(i,j)
60       array(i,j)=-save
c
c        Accumulate elements of inverse matrix
c        -------------------------------------
c
61       do 70 i=1,norder
          if(i-k) 63,70,63
63        array(i,k)=-array(i,k)/amax
70       continue
71       do 80 i=1,norder
          do 80 j=1,norder
           if(i-k) 74,80,74
74         if(j-k) 75,80,75
75         array(i,j)=array(i,j)+array(i,k)*array(k,j)
80       continue
81       do 90 j=1,norder
          if(j-k) 83,90,83
83        array(k,j)=array(k,j)/amax
90       continue
         array(k,k)=1./amax
100     det=det*amax
c
c       Restore ordering of matrix
c       --------------------------
c
101     do 130 l=1,norder
         k=norder-l+1
         j=ik(k)
         if(j-k) 111,111,105
105      do 110 i=1,norder
          save=array(i,k)
          array(i,k)=-array(i,j)
110       array(i,j)=save
111      i=jk(k)
         if(i-k) 130,130,113
113      do 120 j=1,norder
          save=array(k,j)
          array(k,j)=-array(i,j)
120       array(i,j)=save
130      continue
140     return
       end
      SUBROUTINE ATOMIC
c
c     LAST VERSION 9-DEC-1988 15:41:50
c
c     purpose
c     -------
c       To compute atomic concentrations from weight concentrations
c
c     author :
c     --------
c        Pierre TREBBIA, US 41 : "Microscopie Electronique Analytique
c        Quantitative", CNRS, Laboratoire de Physique des Solides, Bat. 510
c        Universite Paris-Sud, F91405 ORSAY Cedex, Phone : (33-1) 69 41 53 68
c
      parameter (nmax=20)
      dimension na(nmax),li(nmax),wconc(nmax),aconc(nmax)
      common /para1/ne,na,li
      common /conc/wconc,aconc
      character lab
      dimension lab(2,82),atwt(82),ro(82)
      common /label/lab
      common /ilp/ilp
      common /wt/atwt
      common /dens/ro
c
      tot=0.
      do 10 i=1,ne
       aconc(i)=wconc(i)/atwt(na(i))
10    tot=tot+aconc(i)
      do 20 i=1,ne
20     aconc(i)=aconc(i)/tot
      write(ilp,30)
30    format(//,15x,'Estimation of Atomic Concentrations',/,
     & 15x,'-----------------------------------',/)
      do 40 i=1,ne
40     write(ilp,50) (lab(j,na(i)),j=1,2),atwt(na(i)),ro(na(i)),
     & aconc(i)*100.
50     format(' Element : ',2a1,3x,'At. Wght. : ',f7.3,3x,
     & 'Density : ',f7.3,3x,'Atomic % : ',f7.4)
      write(ilp,*) ' '
      return
      end
 
      SUBROUTINE THICK
c
c     LAST VERSION 9-DEC-1988 14:09:31
c
c     purpose
c     -------
c       To estimate the maximum thickness beyond which an absorption correction
c     may be needed. This thickness depends on the actual sample composition.
c
c     subroutines and function subprograms required
c     ---------------------------------------------
c     DENSE, ABCOR
c
c     comments
c     --------
c       The criterion is : find Tmax such that for all element couples (i,j),
c           (chi(i) - chi(j)) * ro(specimen) * Tmax < 0.1
c
c     authors :
c     ---------
c       It's a super-coproduction involving :
c
c     1) A. D. ROMIG, JR., DIVISION 1832, SANDIA NATIONAL LABORATORIES,
c        ALBUQUERQUE, NEW MEXICO 87185
c     2) H. O. COLIJN, OHIO STATE UNIVERSITY, CENTRAL ELECTRON OPTICS FACILITY,
c        COLUMBUS, OH 43210
c     3) Pierre TREBBIA, US 41 : "Microscopie Electronique Analytique
c        Quantitative", CNRS, Laboratoire de Physique des Solides, Bat. 510
c        Universite Paris-Sud, F91405 ORSAY Cedex, Phone : (33-1) 69 41 53 68
c
      parameter(nmax=20)
      DIMENSION ABCOF(nmax),ELABCF(nmax,nmax),na(nmax),li(nmax)
      common /coef/abcof,elabcf
      common /para1/ne,na,li
      character lab(2,82)
      common /label/lab
      common /ilp/ilp
      logical full,last,clflag,abflag,dnflag,muflag,limit,muprnt,denchk
      common /flags/full,last,clflag,abflag,dnflag,muflag,limit,muprnt,
     &              denchk
      COMMON/CONST/ BULKRO,PSIR
c
c     do we know the specimen density ?
c     ---------------------------------
c
      if(dnflag) CALL DENSE
c
c     compute mu/ro coefficients ABCOF
c     --------------------------------
c
      CALL ABCOR
c
c     compute the critical thicknesses at level 0.1
c     ---------------------------------------------
c
      tmax=1000.
      con=0.1*sin(psir)/bulkro
      do 110 i=2,ne
       do 110 j=1,i-1
        temp=con/abs(abcof(i)-abcof(j))
        if(temp.lt.tmax) tmax=temp
110   continue
      thmax=tmax*1.0E+8
c
c     print the result
c     ----------------
c
      WRITE(ILP,121) BULKRO
121   FORMAT (/,10X,'The calculated bulk density is',F6.2,1X,
     &  'g/cc',/)
      write(*,122) thmax
      if (ilp.eq.6) WRITE(ILP,122) THmax
122   FORMAT (/,3X,'If the foil, at this analysis site, is greater ',
     & 'than',F9.0,'Angstroms,',/,3x,' an absorption correction may ',
     & 'be required.')
      RETURN
      END
      SUBROUTINE ABCOR
c
c     LAST VERSION 9-DEC-1988 15:56:30
c
c     purpose
c     -------
c       This subroutine calculates the mass absorption coefficients.
c     Elemental coefficients ELABCF are calculated by ABCORM.
c     Given the bulk composition WCONC, the absorption coeffiecents ABCOF in
c     the specimen are determined.
c
c     subroutines and function subprograms required
c     ---------------------------------------------
c     ABCORM
c
c     authors :
c     ---------
c       It's a super-coproduction involving :
c
c     1) A. D. ROMIG, JR., DIVISION 1832, SANDIA NATIONAL LABORATORIES,
c        ALBUQUERQUE, NEW MEXICO 87185
c     2) H. O. COLIJN, OHIO STATE UNIVERSITY, CENTRAL ELECTRON OPTICS FACILITY,
c        COLUMBUS, OH 43210
c     3) Pierre TREBBIA, US 41 : "Microscopie Electronique Analytique
c        Quantitative", CNRS, Laboratoire de Physique des Solides, Bat. 510
c        Universite Paris-Sud, F91405 ORSAY Cedex, Phone : (33-1) 69 41 53 68
c
      parameter (nmax=20)
      DIMENSION ABCOF(nmax),ELABCF(nmax,nmax),wconc(nmax),
     &aconc(nmax)
      common /coef/abcof,elabcf
      COMMON /COnc/wconc,aconc
      dimension na(nmax),li(nmax)
      common /para1/ne,na,li
      character lab(2,82)
      common /label/lab
      common /ilp/ilp
      logical full,last,clflag,abflag,dnflag,muflag,limit,muprnt,denchk
      common /flags/full,last,clflag,abflag,dnflag,muflag,limit,muprnt,
     &              denchk
c
c     if not already done, compute elemental mu/ro coefficients
c     ---------------------------------------------------------
c
      if(muflag) then
       write(ilp,*)' '
       CALL ABCORM
       write(ilp,*)' '
      endif
c
c     calculate the bulk specimen absorption coefficients
c     ---------------------------------------------------
c
      DO 100 I=1,NE
        ABCOF(i) = 0.
        DO 100 J=1,NE
100      ABCOF(i) = ABCOF(i) + wconc(j)*ELABCF(i,j)
c
c     print the result
c     ----------------
c
      if(muprnt) then
       do 200 i=1,ne
200      write(ilp,250) (lab(m,na(i)),m=1,2),abcof(i)
250     format(' Mu/ro(lamda(',2a1,'),specimen) coeff : ',f12.4)
      endif
      RETURN
      END
      SUBROUTINE ABCORM
c
c     LAST VERSION 9-DEC-1988 15:57:30
c
c     purpose
c     -------
c       To compute individual mu/ro coefficient ELABCF
c
c     subroutines and function subprograms required
c     ---------------------------------------------
c     COEFF,EXPON
c
c     author :
c     --------
c      Pierre TREBBIA, US 41 : "Microscopie Electronique Analytique
c      Quantitative", CNRS, Laboratoire de Physique des Solides, Bat. 510
c      Universite Paris-Sud, F91405 ORSAY Cedex, Phone : (33-1) 69 41 53 68
c
      parameter (nmax=20)
      DIMENSION ABCOF(nmax),ELABCF(nmax,nmax),na(nmax),li(nmax)
      common /coef/abcof,elabcf
      common /para1/ne,na,li
      character lab(2,82)
      common /label/lab
      common /ilp/ilp
      logical full,last,clflag,abflag,dnflag,muflag,limit,muprnt,denchk
      common /flags/full,last,clflag,abflag,dnflag,muflag,limit,muprnt,
     &              denchk
      REAL lam,l1,l2,l3,m1,m2,m3,m4,m5
      dimension lam(6,nmax),reli(5,nmax),en(9,nmax)
      common /wave/lam,reli
      dimension ek(82),l1(82),l2(82),l3(82),m1(82),m2(82),m3(82),m4(82),
     &m5(82)
      common /energ/ek,l1,l2,l3,m1,m2,m3,m4,m5
c
c     prevent redundant calculations
c     ------------------------------
c
      muflag=.false.
c
c     get ionization energies
c     -----------------------
c
      DO 20 I=1,NE
       en(1,i)=ek(na(i))
       en(2,i)=l1(na(i))
       en(3,i)=l2(na(i))
       en(4,i)=l3(na(i))
       en(5,i)=m1(na(i))
       en(6,i)=m2(na(i))
       en(7,i)=m3(na(i))
       en(8,i)=m4(na(i))
20    en(9,i)=m5(na(i))
c
c     get the most important line in the series & compute absorption coeff
c     --------------------------------------------------------------------
c
      do 110 i=1,ne
       if((li(i).eq.2).or.(li(i).eq.5)) then
        elam=12.3981/lam(3,i)
       else
        elam=12.3981/lam(1,i)
       endif
       DO 110 J=1,NE
        z=float(na(j)+10)
        do 100 n=1,9
         if(elam.gt.en(n,j)) then
          con=coeff(z,n)
          ex=expon(z,n)
          go to 101
         endif
100     continue
        con=coeff(z,10)
        ex=2.2
101     elabcf(i,j)=con*lam(6,i)**ex
c
c       print the result
c       ----------------
c
        if(muprnt) then
         write(ilp,105) (lab(m,na(i)),m=1,2),(lab(m,na(j)),m=1,2),
     &   elabcf(i,j)
105      format(' Mu/ro(lamda(',2a1,'),',2a1,') coeff : ',f12.4)
        endif
110   CONTINUE
      RETURN
      END
       FUNCTION EXPON(Z,N)
c
c       VERSION 9-DEC-1988 15:43:02
c
c       purpose
c       -------
c        To compute the exponent term to be used in the determination of
c       individual mu/rho coefficients.
c
c       comments
c       --------
c        Z is the atomic number of the absorbing medium, N (= 1 to 9) is
c       the order of the absorbing shell (K,L1,L2,L3,M1,M2,M3,M4,M5),
c       the case N = 10 is used if the X-ray energy is smaller than all
c       the shells of all the present species.
c
c       author :
c       --------
c       J.W. COLBY, in its program MAGIC IV
c
      GO TO(1,2,2,2,4,4,4,4,5,6),N
1     EXPON=2.9308262+Z*(-.020468343+Z*(.0010586311+Z*(-3.7419260E-05+
     &      Z*(7.2044980E-07-5.7357248E-09*Z))))
      GO TO 3
2     EXPON=2.7376804+Z*(-.0020643454+Z*(1.4651705E-04+
     &      Z*(-3.5898183E-06+Z*(2.8732263E-08-8.1864238E-11*Z))))
3     EXPON=FLOAT(IFIX(100.*EXPON+.5))/100.
      RETURN
4     EXPON=2.60
      RETURN
5     EXPON=2.33
      RETURN
6     EXPON=2.22
      RETURN
      END
       FUNCTION COEFF(Z,N)
c
c       VERSION 9-DEC-1988 15:51:01
c
c       purpose
c       -------
c        To compute the multiplicative term to be used in the determination of
c       individual mu/rho coefficients.
c
c       comments
c       --------
c        Z is the atomic number of the absorbing medium, N (= 1 to 9) is
c       the order of the absorbing shell (K,L1,L2,L3,M1,M2,M3,M4,M5),
c       the case N = 10 is used if the X-ray energy is smaller than all
c       the shells of all the present species.
c
c       author :
c       --------
c       J.W. COLBY, in its program MAGIC IV
c
      GO TO(1,2,3,4,5,6,7,8,9,10),N
1     COEFF=-.6396+Z*(.4649-Z*(.1357-Z*(2.404E-2-Z*(9.382E-4-Z*(1.937E-5
     &     -Z*1.552E-7)))))
      GO TO 11
2     COEFF=-.2887+Z*(.1264-Z*(2.218E-2-Z*(1.985E-3-Z*(3.189E-5
     &     -Z*(3.272E-7-Z*1.41E-9)))))
      GO TO 11
3     COEFF=-9.885+Z*(1.749-Z*(.1142-Z*(4.301E-3-Z*(6.378E-5-Z*(5.325E-7
     &      -Z*1.884E-9)))))
      GO TO 11
4     COEFF=-10.531+Z*(1.833-Z*(.119-Z*(4.277E-3-Z*(6.621E-5-Z*(5.606E-7
     &      -Z*1.974E-9)))))
      GO TO 11
5     COEFF=56.0595-Z*(8.3133-Z*(.48688-Z*(1.436E-2-Z*(2.3126E-4
     &      -Z*(1.8852E-6-Z*6.1809E-9)))))
      GO TO 11
6     COEFF=93.4532-Z*(12.0904-Z*(.63557-Z*(1.7339E-2-Z*(2.626E-4
     &      -Z*(2.0483E-6-Z*6.4866E-9)))))
      GO TO 11
7     COEFF=79.245-Z*(10.2893-Z*(.54254-Z*(1.4837E-2-Z*(2.2494E-4
     &      -Z*(1.7553E-6-Z*5.557E-9)))))
      GO TO 11
8     COEFF=68.993-Z*(9.0531-Z*(.48183-Z*(1.3277E-2-Z*(2.0225E-4
     &      -Z*(1.5824E-6-Z*5.0123E-9)))))
      GO TO 11
9     COEFF=2.676E-2-Z*(1.955E-3-Z*(3.973E-6+Z*(5.973E-5+Z*(1.19E-7
     &      -Z*(6.135E-10-Z*1.532E-12)))))
      GO TO 11
10    COEFF=-385.7058+Z*(37.148-Z*(1.4501-Z*(2.9444E-2-Z*(3.2486E-4
     &      -Z*(1.8472E-6-Z*4.2113E-9)))))
11    IF(COEFF.LT.0.) coeff=0.
      RETURN
      END