libcpixe/sfc.c

Go to the documentation of this file.

/***************************************************************************
    Copyright (C) 2007 by Carlos Pascual-Izarra
    carlos.pascual_AT_users.sourceforge.net                                                  *

    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 3 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, see <http://www.gnu.org/licenses/>.

 ***************************************************************************/

#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


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]; 

  if(Za<minK || Zb <minK ) {
    return(0);
  }

  nabs=4; 

  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
      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.18a-h 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] ;
    }
  }
}
*/


int inrange(double value,double min, double max)
{
  if(value<min || value>max) return(0);
  else return(1);
}

---