Title                  :QUANTI
Keywords               :EELS,EDX,AEM,CTEM,STEM,CONCENTRATIONS,STATISTICS
Computer               :DEC LSI-11/73
Operating System       :RSX-11M-PLUS
Programming Language   :DEC-FORTRAN IV
Hardware Requirements  :None
Author(s)              :Dr. Pierre TREBBIA
Correspondence Address :Laboratoire de Physique des Solides
                       :Bat. 510, F91405 ORSAY CEDEX FRANCE
Abstract: The purpose of this program is quite general. It can be used not only 
for the determination of atomic concentrations in EELS experiments, but also, 
for example, for the determination of weight concentrations in EDX 
spectroscopy. For given SIGNAL, BACKGROUND, H = 1 + VARIANCE(BACK)/BACKGROUND 
and F = Interaction-detection yield (in arbitrary units) values, it returns 
concentrations values, concentration covariances and confidence intervals at 
three different confidence levels : 0.1, 0.05 and 0.01.
    A complete description of the algorithm is given in "Unbiased 
Method for Signal Estimation in Electron Energy Loss Spectroscopy...", 
P. TREBBIA, Ultramicroscopy (1988).
Title                  :QUANTI
Keywords               :EELS,EDX,AEM,CTEM,STEM,CONCENTRATIONS,STATISTICS
Computer               :DEC LSI-11/73
Operating System       :RSX-11M-PLUS
Programming Language   :DEC-FORTRAN IV
Hardware Requirements  :None
Author(s)              :Dr. Pierre TREBBIA
Correspondence Address :Laboratoire de Physique des Solides
                       :Bat. 510, F91405 ORSAY CEDEX FRANCE
Description : see comments line in the source code text
References : see comments line in the source code text
Compilation procedure : compile the main program QUANTI and the subroutine 
                        INVMAT
Linking procedure : Link the 2 resulting object files.
Test data : QUANTI.DAT is the output that QUANTI should give 
            if the appropriate parameters are entered.
General comments : If the compiler is not a DEC-FORTRAN IV compiler, do check 
                   the source code : for example, the ! sign in a line means 
                   that the remaining text is comments, or the LOG function is 
                   implicitely changed to DLOG if the variable is double 
                   precision.
Title                  :QUANTI
Keywords               :EELS,EDX,AEM,CTEM,STEM,CONCENTRATIONS,STATISTICS
Computer               :DEC LSI-11/73
Operating System       :RSX-11M-PLUS
Programming Language   :DEC-FORTRAN IV
Hardware Requirements  :None
Author(s)              :Dr. Pierre TREBBIA
Correspondence Address :Laboratoire de Physique des Solides
                       :Bat. 510, F91405 ORSAY CEDEX FRANCE

 
        PROGRAM QUANTI
c
c       LAST VERSION : 3-SEP-1987 10:51:14
c
c       purpose
c       -------
c       estimates 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       1) This program cannot handle more than 21 detected elements.
c       2) Modify the PARAMETER statement to increase this limit in BOTH
c          QUANTI and INVMAT programs.
c       3) All the f values MUST be expressed in the SAME arbitrary unit.
c          The program computes the normalization factor FACT to be multiplied
c          to the f(i) for getting FN(i) = Nelec.Natoms.Probint(i).Probdet(i)
c               with :  Nelec   = Number of incident electrons.
c                       Natoms  = Number of probed atoms in the specimen.
c                       Probint = Probability of electron-atom interaction.
c                       Probdet = Probability of event detection.
c       4) In order to minimize variances, input elements are reorganized
c          so that the last element, whose concentration is calculated by
c          difference to 1, is the one with lowest f.f/var(peak) value.
c       5) The program performs a test on detection limit and prompts the
c          advice to check the spurious cases with the programm LIMIT.
c       6) The program allows to test different f values for the same
c          (peak,background,h) set.
c       7) Type CTRL Z to exit.
c       8) In order to get an executive task, one must :
c               a) compile the following FORTRAN programs :
c                       QUANTI,INVMAT
c               b) link the resulting object files.
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=21)
        double precision sigma,rtot,cftot,ptot,fact,test,try
        dimension sigma(nel-1,nel-1),nb(nel),p(nel),b(nel),h(nel),
     &  f(nel),fn(nel),r(nel-1),vp(nel),c(nel),sc(nel),z(3),res(3,nel)
        logical pr,again
        norder=nel-1                    ! order of matrix sigma
        pr=.true.                       ! enable record of data
        again=.false.                   ! disable the test of others f(i)
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       Results storage ?
c       -----------------
c
        type 2
2       format(' Do you want a report ? [default : Y] ? ',$)
        read(5,1010,end=800) rep
        if ((rep.eq.'n').or.(rep.eq.'N')) pr=.false.
        if(pr) then
         open(4,name='QUANTI.DAT',status='new')
         type *,' ** Your report is in the file QUANTI.DAT **'
        endif
c
c       Ask for input values and compute peak variance
c       ----------------------------------------------
c
5       write(5,10) nel
10      format(3x,'How many detected elements (max:',i3,') ? ',$)
        read(5,*,end=800) nd
        if((nd.gt.nel).or.(nd.lt.2)) go to 5
        nfree=nd-1                      ! degree of freedom
        do 30 i=1,nd
         nb(i)=i
15       write(5,20) nb(i)
20       format(' Element # ',i2,' : Input P,B,h,f ',$)
         read(5,*,end=800) p(i),b(i),h(i),f(i)
         if(h(i).lt.1.) then
          type *,' *** Incorrect value for h ***'
          go to 15
         endif
         vp(i)=p(i)+h(i)*b(i)
         if(p(i)*p(i)/vp(i).lt.13.) then
          write(5,25) nb(i)
          if(pr) write(4,25) nb(i)
25        format(' **** Caution : Check the detection limit of element',
     &    ' # ',i2,' Use program LIMIT !')
         endif
30      continue
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
         write(5,45) k,nd
         if(pr) write(4,45) k,nd
45       format(/,' *** Warning : elements ',i2,' and ',i2,' 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
        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
        type *,' Coffee break ...'
        call invmat(sigma,norder,det)
        type *,' .... 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       Output the results
c       ------------------
c
        write(5,400)
        if(pr) write (4,400)
400     format(/,3x,'Estimation of concentrations',
     &  /,3x,'----------------------------',
     &  //,' Elem',6x,'Peak',4x,'Background',4x,'h',7x,
     &  'Var(P)',6x,'Conc',8x,'Sigma(Conc)',10x,'f')
        do 410 i=1,nd
         write(5,405) nb(i),p(i),b(i),h(i),vp(i),c(i),sc(i),f(i)
         if(pr) write (4,405) nb(i),p(i),b(i),h(i),vp(i),c(i),sc(i),f(i)
405      format(i4,2f12.0,f8.3,f12.0,3e14.5)
410     continue
        write(5,430)
        if(pr) write(4,430)
430     format(/,3x,'Estimation of covariances',/,3x,
     &  '-------------------------')
        do 440 j=1,nfree-1
         do 440 i=j+1,nfree
          write(5,435) nb(j),nb(i),sigma(j,i)
          if(pr) write(4,435) nb(j),nb(i),sigma(j,i)
435       format(' Covariance(',i2,',',i2,') = ',e14.5)
440     continue
        write(5,450)
        if(pr) write(4,450)
450     format(/,3x,'Estimation of confidence intervals',
     &  /,3x,'----------------------------------',//,
     &  5x,'At the confidence level alpha = 0.10 :',/)
        do 465 i=1,nd
         write(5,460) nb(i),c(i),res(1,i)
         if(pr) write(4,460) nb(i),c(i),res(1,i)
460      format(' Element',i3,' Interval :',e14.5,'  + -',e14.5)
465     continue
        write(5,470)
        if(pr) write(4,470)
470     format(/,5x,'At the confidence level alpha = 0.05 :',/)
        do 475 i=1,nd
         write(5,472) nb(i),c(i),res(2,i)
         if(pr) write(4,472) nb(i),c(i),res(2,i)
472      format(' Element',i3,' Interval :',e14.5,'  + -',e14.5)
475     continue
        write(5,490)
        if(pr) write(4,490)
490     format(/,5x,'At the confidence level alpha = 0.01 :',/)
        do 495 i=1,nd
         write(5,492) nb(i),c(i),res(3,i)
         if(pr) write(4,492) nb(i),c(i),res(3,i)
492      format(' Element',i3,' Interval :',e14.5,'  + -',e14.5)
495     continue
        write(5,500)
        if(pr) write(4,500)
500     format(/,'**************************************************',//)
c
c       Test others f(i) ?
c       ------------------
c
        again=.false.
        write(5,600)
600     format(/,' Do you want to change f values [default : N] ? ',$)
        read(5,1010,end=800) new
        if ((new.eq.'y').or.(new.eq.'Y')) again=.true.
        if(again) then
         do 640 i=1,nd
          write(5,620) nb(i),f(i)
620       format(' Element # ',i2,'  Old f :',f9.2,' Input new f : ',$)
          read(5,*,end=800) f(i)
640      continue
         go to 35
        endif
c
c       Run again
c       ---------
        type *,' Start again. Input new data. Enter CTRL Z for exit'
        go to 5
c
c       Program exit
c       ------------
c
800     if(pr) then
         close(5)
         open(5,status='new')
         write(5,805)
805      format(10x,'Do you want your report to be printed ? ',
     &   '[default : N] ',$)
         read(5,1010,end=1000) pri
         if((pri.eq.'y').or.(pri.eq.'Y')) then
          close(4,dispose='print')
          type *,'Don''t forget your listing !'
         endif
        endif
1000    stop
1010    format(a1)
        end


 
        SUBROUTINE INVMAT(ARRAY,NORDER,DET)
c
c       LAST VERSION : 15-JUL-1987 09:24:21
c
c       purpose
c       -------
c       inverts 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       det -----> determinant of input matrix
c
c       subroutines and function subprograms required
c       ---------------------------------------------
c       none
c
c       comments
c       --------
c       1) Dimension statement valid for norder up to nmax = 20.
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
        PARAMETER (NMAX=20)
        double precision array,amax,save
        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


                     DEMONSTRATION OF PROGRAM QUANTI

 Caution : Check the detection limit of element #  2 Use program LIMIT !

  Estimation of concentrations
  ----------------------------

Elem  Peak    Background    h       Var(P)     
 1  299507.     663668.   4.320    3166553.  
 2    6000.     663668.   4.320    2873046.  
 3  299507.     663668.   4.320    3166553.  
 4  299507.     663668.  12.000    8263523.  

   Conc        Sigma(Conc)         f
 0.42492E+00   0.22353E-02   0.10000E+01
 0.85125E-02   0.21562E-02   0.10000E+01
 0.14164E+00   0.83138E-03   0.30000E+01
 0.42492E+00   0.26934E-02   0.10000E+01

  Estimation of covariances
  -------------------------
Covariance( 1, 2) =   -0.12495E-05
Covariance( 1, 3) =   -0.15302E-06
Covariance( 2, 3) =   -0.13884E-06

  Estimation of confidence intervals
  ----------------------------------

    At the confidence level alpha = 0.10 :

Element  1 Interval :   0.42492E+00  + -   0.36771E-02
Element  2 Interval :   0.85125E-02  + -   0.35470E-02
Element  3 Interval :   0.14164E+00  + -   0.13676E-02
Element  4 Interval :   0.42492E+00  + -   0.44306E-02

    At the confidence level alpha = 0.05 :

Element  1 Interval :   0.42492E+00  + -   0.43812E-02
Element  2 Interval :   0.85125E-02  + -   0.42262E-02
Element  3 Interval :   0.14164E+00  + -   0.16295E-02
Element  4 Interval :   0.42492E+00  + -   0.52790E-02

    At the confidence level alpha = 0.01 :

Element  1 Interval :   0.42492E+00  + -   0.57559E-02
Element  2 Interval :   0.85125E-02  + -   0.55522E-02
Element  3 Interval :   0.14164E+00  + -   0.21408E-02
Element  4 Interval :   0.42492E+00  + -   0.69354E-02

*************************************************