sfc.c

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2005 by Carlos Pascual-Izarra                           *
 *   carlos.pascual@itn.pt                                                 *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include "libcpixe.h"
#include "compilopt.h"  //here we store variables for the C preprocessor for conditional compiling


/**\file sfc.c
sfc.c contains the portion of libcpixe related to secondary fluorescence corrections.

This is still (in v1.0) under development so nothing is really used yet.
*/


int createSFCList(const EXP_PARAM *pexp, const int *PresentElems, const XrayYield *XYldArray, const AbsCoef *TotAbsCoefArray, int *nSFCList, SFCListElem **SFCList)
{
  int Za,Zb;
  int dim;
  
  *nSFCList=0;
  *SFCList=NULL;
  
  dim=pexp->simpar.MaxZinsample+1;
  
  for(Zb=minK;Zb<dim;Zb++){
    if(PresentElems[Zb]){ 
      for(Za=minK;Za<dim;Za++){
        if(PresentElems[Za]){
          needSFC(Za, Zb, &XYldArray[Za].ener, &TotAbsCoefArray[Zb], nSFCList, SFCList);
        }
      }  
    }
  }
  return(0);
}

int needSFC(atomicnumber Za, atomicnumber Zb, const CalibYld *enerA, const AbsCoef *AbsC, int *nSFCList, SFCListElem **SFCList )
{
  int j,l,i,ib,nabs;
  int n;
  double enerB,enerB_K,enerB_LIII;
  SFCListElem tempList[135]; ///<@todo 135=9x15 (absorption lines x emission_lines). If more lines or absorptions are added, this should be raised. Note that a conditional compilation for the case in which M absorption lines are not used, could reduce it to 4x15 or even to 4x12 if M emission lines are neither considered.
  
  if(Za<minK || Zb <minK ) {
    return(0);
  }
  
  
  nabs=4; ///<@todo Only K and L absorption lines are considered (since M lines are hard to quantify anyway, we don't waste time). If SFC is desired for M lines too, change nabs to 9 here. Also, a conditional compilation could be implemented depending on a preproc variable called ,i.e., QUANTIFYMLINES
  
  if(Zb<minL)nabs=1;
  
  enerB_K=AbsC->enr[0];
  if(enerB_K<0.001)enerB_K=HUGE;
  enerB_LIII=AbsC->enr[4];
  if(enerB_LIII<0.001)enerB_LIII=HUGE;
  
  //loop through all the absorption energies of the B element
  for (n=0, ib=1;ib<=nabs;ib++){
    enerB=AbsC->enr[ib]; //energy of absorption edge
    if(enerB>0){
      //K emission lines:
      for (j = 1; j <= 3; j++) {
        if(enerA->K_[j]>enerB){
          tempList[n].Za=Za;
          tempList[n].Zb=Zb;
          tempList[n].epri=0;  //this means K lines
          tempList[n].esec=j;   //this is the index inside the K group 
          tempList[n].abs=ib;  //Absorptions for: ib=1-->K ; ib=2,..4--> L ; ib=5,..9-->M 
          tempList[n].Ea=enerA->K_[j];
          tempList[n].sigmaphoto=sigmaphoto(ib, AbsC, Zb, enerA->K_[j]);
          n++;
        }
      }
      //L emission lines
      if(Za>minL) for (j = 1; j <= 3 ; j++) for (l = 1; l <= 3; l++) {
        if(inrange(enerA->L_[j][l],enerB,enerB_K)) {
          tempList[n].Za=Za;
          tempList[n].Zb=Zb;
          tempList[n].epri=j;  //j=1-->LIII, j=2-->LII, j=3-->LI 
          tempList[n].esec=l;   //this is the index inside the Lj group
          tempList[n].abs=ib;
          tempList[n].Ea=enerA->L_[j][l];
          tempList[n].sigmaphoto=sigmaphoto(ib, AbsC, Zb, enerA->L_[j][l]);
          n++;
        }
      }
      //M emission lines
      ///@todo In case M emission lines should not be considered for SFC, this for loop can be eliminated by a conditional compilation
      if(Za>minM) for (j = 1; j <= 3; j++) {
        if(inrange(enerA->M_[j],enerB,enerB_LIII) ) {
          tempList[n].Za=Za;
          tempList[n].Zb=Zb;
          tempList[n].epri=4;  //this means M lines 
          tempList[n].esec=j;   //this is the index inside the M group 
          tempList[n].abs=ib;
          tempList[n].Ea=enerA->M_[j];
          tempList[n].sigmaphoto=sigmaphoto(ib, AbsC, Zb, enerA->M_[j]);
          n++;
        }
      }
    }
  }
  
  //If any pair of emission of A and absorption by B matched the criteria, reallocate the SFClist and write the new entries.
  if(n>0){
    *SFCList =(SFCListElem*)realloc(*SFCList, (*nSFCList + n)*sizeof(SFCListElem));
    if (*SFCList==NULL){fprintf(stderr,"\n Error allocating memory (SFCList)\n");exit(1);}
  }
  for(i=0;i<n;i++)(*SFCList)[*nSFCList+i]=tempList[i];
  
  *nSFCList+=n;
  
  return(n); 
}

/*(see) eq. 2.18 of ref [1]*/
double sigmaphoto(int iabs, const AbsCoef *AbsC, atomicnumber Zb, double Xray)
{
  double result;
  double abs;
  int iener;
  int i;
  double stdener[23] = {      1.00, 1.25, 1.50, 1.75, 2.00, 
                        2.50, 3.00, 3.50, 4.00, 4.50, 5.00, 
                        5.50, 6.00, 6.50, 7.00, 8.00, 10.0,
                        12.5, 15.0, 17.5, 20.0, 25.0, 30.0};
  
  i=iabs-1;
  
  /*  (S-1)/S = (1-1/S) = ( 1 - 1/(coefener[1]/coefener[0]) ) = ( 1 - coefener[0]/coefener[1] ) */
  
  switch(iabs){
    case 1:  //K absorption
      result=( 1.-AbsC->coefenr[0][0]/AbsC->coefenr[0][1] );
    break;
    case 2:  //LI absorption
      result=( 1. - AbsC->coefenr[1][0]/AbsC->coefenr[1][1] );
    break;
    case 3:  //LII absorption
      result=( 1.-AbsC->coefenr[2][0]/AbsC->coefenr[2][1] );
      if(Xray > AbsC->enr[2]) result*= AbsC->coefenr[1][0]/AbsC->coefenr[1][1]; //Xray > E(LI)
    break;
    case 4:  //LIII absorption
      result=( 1.-AbsC->coefenr[3][0]/AbsC->coefenr[3][1] );
      if(Xray > AbsC->enr[3]) result*= AbsC->coefenr[2][0]/AbsC->coefenr[2][1]; //Xray > E(LII)
      if(Xray > AbsC->enr[2]) result*= AbsC->coefenr[1][0]/AbsC->coefenr[1][1]; //Xray > E(LI)
    break;
    case 5:  //MI absorption
    case 6:  //MII absorption
    case 7:  //MIII absorption
    case 8:  //MIV absorption
    case 9:  //MV absorption
      result=( 1.-AbsC->coefenr[4][0]/AbsC->coefenr[i][1] ); 
    break;
  }
  
  //locate the "std energy index"
  for(iener=1;Xray > stdener[iener];iener++);
  
  //calculate the absorption cross section by the element Zb 
  abs=TotAbsor_elemental(iener, AbsC, Zb, Xray);  
  
  //return the sigmaphoto value
  return(result*abs);  
}


/* *Calculates the probability for an X-ray (from A) being absorbed (by B) and produce a  from a given line of the element B af
 * 
 * See eq. (2.5) of ref [1]
 *
 *Ref. [1]: "Correcoes de Fluorescencia Secundaria em PIXE", M.A. Reis, Master Thesis (1993). 
 *
 * @param TotAbsCoefArray 
 * @param pSFC 
 * @param plyr 
 * @return Rba*mu(A) in cm2/1e15at
 */
/*
double SFC_RbaMu(const AbsCoef *TotAbsCoefArray, const SFCListElem *pSFC, const LYR *plyr)
{
}
*/

/*
void Rbeta(Char *FAbNm, FluorCKCoef FCK, atomicnumber AtNTra,
     double XEner, CalibYld *Rba)
{
  int i, j;
  double SigBA[10];

  memcpy(&SigBA[1], &NineArrayNul[1], sizeof(double) * 9);
      //   1-K  2-LI  3-LII  4-LIII   5...-M... 
  FotSec(FAbNm, AtNTra, FCK.ck[0], FCK.ck[1], FCK.ck[2], XEner, &SigBA[1]);
  *Rba = CYldNul;
  
  if (SigBA[1] > 0) {
    for (i = 1; i <= 3; i++)
      Rba->K_[i] = SigBA[1] * FCK.w[1] * FCK.k.K_[i] ;
  }
  if (SigBA[5] > 0) {
    for (i = 1; i <= 3; i++)
      Rba->M_[i] = SigBA[5] * FCK.w[5] * FCK.k.M_[i] ;
  }
  for (i = 1; i <= 3; i++) {
    if (SigBA[5 - i] > 0) {
      for (j = 1; j <= 3; j++)
        Rba->L_[i][j] = SigBA[5 - i] * FCK.w[5 - i] * FCK.k.L_[i][j] ;
    }
  }
}
*/


/**Checks if a value is within limits
 * 
 * @param value (I) Value to check 
 * @param min (I) minimum limit 
 * @param max (I) maximum limit
 * @return 1 if value is inside range, 0 otherwise
 */
int inrange(double value,double min, double max)
{
  if(value<min || value>max) return(0);
  else return(1);
}

---