libcpixe.h

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2005 by Carlos Pascual-Izarra & Miguel A. Reis          *
 *   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.             *
 ***************************************************************************/

/**\file libcpixe.h
libcpixe.h is the header file for the LibCPIXE library.
All constants, structs, and functions are declared here.
This must be included by any program using libcpixe.c
*/

/** @todo Modify fscanf functions to check for length when reading filenames of fixed maximum length
*/

#include "compilopt.h"

#define LONGSTRINGLENGTH 256
#define FILENMLENGTH LONGSTRINGLENGTH

//conversion factor from degrees to rads
#define DEG2RAD (M_PI/180.0)
//Energy (in MeV) corresponding to 1 amu through E=mc2
#define MEVAMU 931.52181 
#define UATMASS (1./1822.887)

//Barn to m2 conversion
#define BARNTOM2  1.0e-28


typedef char STATFILENAME[FILENMLENGTH] ; ///<Fixed length string for use with filenames

typedef char ChemSymb[3];///<Fixed length string for use with chem. symbols

typedef int atomicnumber;///<Convenient definition of atomic numbers as ints

//typedef double tenarray[11];

//typedef double ninearray[10];

//typedef double threearray[4];

//typedef double LArray[4][4];

//typedef char LineNmType[6][4][11];

typedef double SecXL[4];

typedef enum { K, L, M } FwTipo;

/*Structures*/

/**Chemical element*/
typedef struct ELEMENT
  {
  char symbol[3]; ///<Chemical symbol
  int Z;          ///<Atomic number
  int A;          ///<Number of nucleons
  double M;       ///<Natural Atomic weight, in amu 
  double IM;      ///<Isotopic mass, in amu
  }ELEMENT;

/**Compound formed by various elements*/
typedef struct COMPOUND
  {
  char name[30];  /**<Name of the compound*/
  double sumX;    /**<sum of all ammounts in the compound (useful for normalizing)*/
  int nelem;      /**<number of elements present in the compound*/
  ELEMENT *elem;  /**<Pointer to an array of elements*/
  double *X;      /**<Pointer to an array containing the atomic ammount of each element*/
  double *xn;     /**<Pointer to an array containing the NORMALIZED atomic ammount of each element*/
  double *w;      /**<Pointer to an array containing the NORMALIZED MASS ammount of each element**/
  }COMPOUND;

/**Stopping Power Table.
  This is a modified version from that in hotstop.h*/
typedef struct SPT
  {
  char ID[30];
  int nrows;    ///<number of rows of the table (dimension)
  int logmode;  ///<flag. 0 if constant Energy step (linear) and 1 if Energy is logarithmic 
  double Emax;  ///<Maximum energy (minimum is always 0)
  double Estep; ///<Step between energies, in keV (log step in case of logmode=1)
  double Smax;  ///<Maximum valu of the stopping power in the region E=(0,Emax)
  double *E;    ///<CAUTION: E will be assumed to have following units: keV
  double *S;    ///<CAUTION: S will be assumed to have following units: keV/1e15atcm-2
  double *dSdE; ///< derivative for linear interpolations: dSdE[i]=(S[i+1]-S[i])/(E[i+1]-E[i]) 
  }SPT;

/**Filter description
@todo document this structure*/
typedef struct absorber
  {
  int number;
  char name[6];
  int nreg;
  double ereg[11];
  double coef[11][3];   
  double h[3];
  } absorber;


/**Simulation parameters*/
typedef struct SIM_PARAM
  {
  int MaxZinsample;  ///<larger Atomic number present in the whole sample
  int AllowSXFCorr;   ///<Flag for secondary fluorescence calculation.1=true, 0=false
  int AllowXEqCalc;    ///<Flag for X-ray equivalent calculation. 1=true, 0=false
  double DTCC;        ///<Dead Time Correction coefficient
  //double DTCCerr;     ///<Error in Dead Time Correction coefficient (Not really used)
  double ColCharge; ///<Collected charge (fluence*charge_state), in uC.
  double CalEner;   ///<Energy at which the calibration was done.
  int useFilter; ///<Flag: whether a filter on the detector is used (1) or not (0).
  }SIM_PARAM;

  
      
/**Experimental parameters*/
typedef struct EXP_PARAM
  {
  ELEMENT ion;      ///<Beam type @todo For the moment only H ions are supported.
  double BeamEner;  ///<Beam Energy, in keV
  double BeamCol;   ///<colimator aperture (diameter of the beam), in mm
  double DetColFac; ///<Relative diameter of det colimator to that used during calib
  double IncAng;    ///<Incident angle (beam to surface normal), in degrees 
  double DetAng;    ///<Detector Angle (detector to surface normal), in degrees
  double cosDet;    ///<For optimization. Cos(IncDet)
  double cosInc;    ///<For optimization. Cos(IncAng)
  double cosFac;    ///<For optimization. Cos(IncAng)/Cos(IncDet)
  double BeamCross; ///<Beam cross section, in cm2
  double FinalEner; ///<Minimum energy considered, in keV
  SIM_PARAM simpar; ///<Collection of parameters regarding to simulation
  }EXP_PARAM;

/**Structure for passing info read during from input files*/
typedef struct EXTRAINFO  
  {
   char *SampleFileNm;  ///<pointer to name of sample definition file
   char *FilterFileNm;  ///<pointer to name of filter definition file
   char *XYldFileNm;    ///<pointer to name of detector efficiency calibration file
   char *AreasFileNm;   ///<pointer to name of areas file (also used for calculation flags)
   char *DBpath;        ///<pointer to Path name for databases (for stoppings, FCK coefs, etc)
   char *OutputFileNm;   ///<pointer to name of output file
   int AreasFormat;     ///<Defines the format of the Areas File. 1 for raw areas format, 2 for DATTPIXE Areas Format
   int WantOutputfile;  ///<Flag indicating if an Output file is required (1) or not (0)
   int DoCalibration;   ///<Flag indicating that the input data is for performing a calibration (1) or not (0)
  }EXTRAINFO;

/**Peak areas*/
typedef struct PeakArea {
  ChemSymb symb;   ///< Symbol
  atomicnumber atnum; ///<Atomic number
  int area[5][3]; ///< Peak areas: [0][j]-->K  [1][j]-->LIII , [2][j]-->LII , [3][j]-->LI , [4][j]-->M
  int erro[5][3]; ///< Error in peak areas (same convention as in ".area")
} PeakArea;

/**L,K & M lines intensities*/
typedef struct CalibYld {
  double K_[4];   ///< K Lines 
  double L_[4][4];   ///< 1-LIII , 2-LII , 3-LI 
  double M_[4];   ///< M  Lines
} CalibYld;

/** Emission energies and other ... 
@todo document this structure*/
typedef struct XrayYield {
  ChemSymb symb;
  atomicnumber atnum;
  CalibYld ener;
  CalibYld XYld;
} XrayYield;

/**Fluorescence and Coster-Kronig coefficients*/
typedef struct FluorCKCoef {
  atomicnumber atnum;  
  double w[10];         
  double ck[3];        
  CalibYld k;
} FluorCKCoef;

/**@todo document this structure*/
typedef struct SRType {
  CalibYld AbsFact;
  CalibYld Xeq;
  CalibYld SFCr;
  CalibYld DetLim;
  CalibYld ConcErr;
  CalibYld MDV;
} SRType;


/**@todo document this structure*/
typedef struct SecResType {
  PeakArea Pk;
  SRType SR;
} SecResType;



/**@todo document this structure*/
typedef struct ESxType {
  double ep;
  double stpp;
  double x;
} ESxType;
 
/**Basic sample foil description*/
typedef struct foil {
  double thick; ///<Thickness, in 1e15at/cm2
  int nfoilelm; ///<Number of elements (redundant with foil.comp.nelem)
  COMPOUND comp; ///<compound definition
} foil;

 
/**Variables that depend on the layer*/  
typedef struct 
  {
  int NumOfTrc;                ///<Number of "traces" (i.e. elements)
  int *TrcAtnum;               ///<Atomic numbers for each trace.
  int dimTrans;                ///<Dimension of Transmission vector
  CalibYld *Trans;             ///<Vector of calculated transmission for the lines of each element in this layer.
//  PeakArea *AreasArray;        ///<Experimental Peak Areas. Not used for simulation. @todo check if it makes sense to maintain it here
//  XrayYield *ResArray;        ///<@todo description?
  XrayYield *ResYldArray;     ///<Results (areas and energies of the calculated lines)
  SecResType *SecResArray;    ///<Same as ResYldArray but for secondary fluorescence results only
  int *TrcUse;          ///<Pointer to an array of flags controlling if calcs should be done for each given "trace" of this layer
  int *NeedSFC;         ///<Pointer to an array of flags controlling if sec. fluor. calcs should be done for each given element present in the SAMPLE
  ESxType *ESxArray;    ///<Array containing Energy, Stopping and position of each sublayer. 
  int FESxlen;          ///<Number of sublayers created for this layer (dimension of ESxArray)
  CalibYld *SSTrsArray; ///<Array containing Transmission coefs for each element in each sublayer
  CalibYld *SSYldArray; ///<Array containing X-ray production cross section for each element in each sublayer
  double ThickIn;       ///<In-going path length for the layer (and the curren geometry)
  double FoilInEner;    ///<Energy (keV) when entering the layer
  double FoilOutEner;    ///<Energy (keV) when exiting the layer
  double absolutepos;   ///<absolute in-going path length, in 1e15at/cm2, from sample surface. i.e. sumation of ThickIn for all layers below this one.
  const foil *pFoil;          ///<pointer to the foil associated to this layer (see foil structure)
  } LYR;
  

/**X-ray Absortion coefficients*/
typedef struct AbsCoef {
  atomicnumber atnum;     ///<Atomic number
  double coefstd[23];     ///<23 absorption coefficients, each one for a given standard energy tabulated in stdtable[] in TotAbsor() 
  double enr[10];          ///<Array transition energies (energies at which a sharp edge in absorption occurs): from 1 to 9: K LI LII LIII MI MII MIII MIV MV  (absorption E is decreasing)
  double coefenr[9][2];   ///<abs coefs for E just below (coefenr[i][0]) and just above (coefenr[i][1]) the transition energy enr[i].
} AbsCoef;

/**Definition of absorber filter. The filter can be defined in 2 alternative ways: 1.- Just like a you define a sample,  layer by layer (the transmission being calculated for each relevant line). 2.- Defining its transmission function by giving a series of coefficients that fit T(E) (this was the default way in DATPIXE).*/   
typedef struct FILTER {
  int nlyr;         ///<number of layers defining the sample
  int MaxZ;         ///<maximum atomic number of elements present in sample
  int *FilterElems;  ///<Vector of flags indicating if a certain element is present in the filter (1) or not (0)
  double geomcorr;  ///<Geometrical correction to foil thickness (e.g. if filter is tilted: geomcorr=1/cos(tilt))
  foil *foil;       ///<definition of each layer
  int dimTrans;     ///<Dimension of Transmission vector
  CalibYld *Trans;  ///<Vector of calculated transmission for the lines of each element in the whole filter
  int changes;      ///<Flag indicating that Filter has changed (1) or not (0)
  int FilterType;   ///<Flag indicating type of filter used: 0->No filter 1->Foil-defined  2->Calibrated 
  absorber Calibflt;///<Structure to hold coeffs if alternative #2 is chosen. @todo (not supported yet)
}FILTER;

/**Calculated areas for a given element.*/
typedef struct CPIXERESULTS{
  int atnum;          ///<atomic number of the element
  CalibYld simareas;  ///<K ,L & M yields 
  CalibYld err;       ///<K ,L & M yield uncertainties 
}CPIXERESULTS;



/**Each element which fluoresces knows which elements (and which lines) are the origin*/
typedef struct SFCListElem {
  atomicnumber Za;  ///<Z of emissor of primary X rays
  atomicnumber Zb;  ///<Z of the element that fluoresces (secondary emissor)
  int epri;         ///<primary number identifying the emission group: 0=K, 1=LIII, 2=LII, 3=LI, 4=M
  int esec;         ///<secondary number identifying the emission line inside the group 
  int abs;          ///<Index of absorption energies. From 1 to 9: K LI LII LIII MI MII MIII MIV MV
  double Ea;        ///<Energy of the emission line (in keV)
  double sigmaphoto;///<photoemission cross section (@todo units?)
} SFCListElem;

/*Constants and Initializators*/  
/**Limits for atomic numbers to be used
 @todo These values should be read instead of hardcoded. Maybe in SIMPARAM?
 */
//static atomicnumber minK = 11, maxK = 54, minL = 30, maxL = 92, minM = 60;

static const atomicnumber minK = 11, maxK = 54, minL = 30, minM = 60;

static const CalibYld CYldNul = { { 0.0, 0.0, 0.0, 0.0 },
                            { { 0.0, 0.0, 0.0, 0.0 },
                              { 0.0, 0.0, 0.0, 0.0 },
                              { 0.0, 0.0, 0.0, 0.0 },
                              { 0.0, 0.0, 0.0, 0.0 }},
                            { 0.0, 0.0, 0.0, 0.0 }};
  
/**Names for X-ray emmission lines*/

//static LineNmType LineNm = {
static char LineNm[6][4][11] = {
  { "--", "  --      ", "  --      ", "  --      " },
  { "--", "K_alpha1,2", "K_beta_1  ", "K_beta_2  " },
  { "--", "L_alpha1,2", "L_beta_2  ", "L_l       " },
  { "--", "L_beta_1  ", "L_gamma_1 ", "L_eta     " },
  { "--", "L_beta_3  ", "L_beta_4  ", "L_gamma_3 " },
  { "--", "M_alpha1,2", "M_beta    ", "M_gamma   " }
};

  
  
/*Function declarations*/
  
void testlib(int **testarray,int *n);

int symbol2Z(char *symbol);

int Z2mass(int Z, double *mass, char option);

int readCOMPOUND(FILE *f, int nelem, COMPOUND *c);

int readINPUT(const char *InputFileNm, EXP_PARAM *pexppar, EXTRAINFO *pExtraInfo);

int readEXP_PARAM(char *ExpParFileNm, EXP_PARAM *pexppar);

int readSIM_PARAM(char *SampleParamFileNm, SIM_PARAM *simpar,char *IniMatFileNm, char *FilterFileNm, char *XYldFileNm);

int readsample(char *SampleDefFileNm, int *MaxZ, int *NFoil, foil **Sample);

int readFilter(const char *FilterDefFileNm, FILTER *Filter);

int createPresentElems(int MaxZ, int NFoil,const foil *MatArray, int **PresentElems);

int readCalcFlags(const char *CalcFlagsFileNm, const int *PresentElems, int MaxZinsample, int AreasFormat, CPIXERESULTS **CalcFlags);

int readXYld(const char *XYldFileNm, const int *PresentElems, const CPIXERESULTS *CalcFlags, int MaxZinsample, double *CalEner, XrayYield **XYldArray);

int readFCK(char *FCKCoefFileNm, int MaxZinsample, FluorCKCoef **FCKCoefArray);

int readAbsCoef(char *TotAbsCoefFileNm, int MaxZinsample, const int *PresentElems, const FILTER *Filter, AbsCoef **TotAbsCoefArray);

int createSPTs(const char *path, const EXP_PARAM *pexp, int MaxZinsample, const int *PresentElems, double step, SPT **SPTArray);

double getSP(const COMPOUND *pcmp, const SPT *SPTArray, double E);

int initlyrarray(const EXP_PARAM *pexp, const foil *MatArray, const XrayYield *XYldArray, const AbsCoef *TotAbsCoefArray, const SPT *SPTArray, const int *PresentElems, int NFoil, int *NFoilUsed, LYR **plyrarray);

int initlyr(const EXP_PARAM *pexp, const XrayYield *XYldArray, const AbsCoef *TotAbsCoefArray, const int *PresentElems, LYR *plyr, double *MACoef);

int FilterTrans(const EXP_PARAM *pexp, const XrayYield *XYldArray, const AbsCoef *TotAbsCoefArray, const int *PresentElems, FILTER *Filter);

int Transmission(const XrayYield *XYldArray, const AbsCoef *TotAbsCoefArray, int Z, const foil *pFoil, double geomcorr, const CalibYld *pTransOld, CalibYld *pTrans);

double TotAbsor(const AbsCoef *TotAbsCoefArray, const COMPOUND *cmp, double Xray);

double TotAbsor_elemental(int iener, const AbsCoef *Absco, int Z, double Xray);

int createsublyrs(const EXP_PARAM *pexp, const foil *pMat, const SPT *SPTArray, double MajAbsCoef, LYR *plyr);

int SSThick(const EXP_PARAM *pexp, const foil *pMat, const SPT *SPTArray, double MajAbsCoef, double LayerThickness, ESxType ESxin,
        ESxType *ESxfin, int *pFpos, double *thick, ESxType **ESxA);

int integrate_Simpson(const EXP_PARAM *pexp, const AbsCoef *TotAbsCoefArray, const FluorCKCoef *FCKCoefArray, const XrayYield *XYldArray, int NFoilUsed, const FILTER *Filter, LYR *plyrarray, CalibYld *XYldSums);

int Xprod(const AbsCoef *AbsC, const FluorCKCoef *pFCK, atomicnumber Z1, atomicnumber Z2, double M2 , double ener, CalibYld *XYld);

double PaulX(double ener, atomicnumber z);

double PaulX_y(double MeV, atomicnumber z);

double PaulX_lp(double x, long p);

int ReisX(const AbsCoef *AbsC, const FluorCKCoef *pFCK, double ener, atomicnumber z, double M2, double *sigmaXL);

double ReisX_gs(FwTipo T, int SS, double ksis);

double ReisX_hs(FwTipo T, int SS, double ksih, double thet);

double ReisX_En(double z, int niu);

double ReisX_polisec(int ssind, double kz, double tz);

double ReisX_g(int ss, atomicnumber Zg, double xi);

void PenInteg(atomicnumber atnumb, const CalibYld *absc, const ESxType *ESA,
        const CalibYld *YldA, const CalibYld *TrsA, const CalibYld *pTrs0, 
        int FExlen, int NeedSFC, int AllowXEqCalc,
        double x0, double CosInc,
        CalibYld *XYld, CalibYld *XSFCr, CalibYld *XYldxmed);

double Simps(double a, double b, double fa, double fi, double fb);

void deNormalize(const EXP_PARAM *pexp, const AbsCoef *AbsC, 
               const FluorCKCoef *pFCK, atomicnumber Z2, double M2,  double attfraction, 
               const CalibYld *pXYld, XrayYield *ResY, CalibYld *XYldSum);

void freeFilter(FILTER *Filter);    
               
void freeReusable(int NFoils, LYR **plyrarray, foil **MatArray, CalibYld **XYldSums );

void safefree(void **ptr);

void fprintCALIBYLD(FILE *f,const CalibYld *pCYld);

int createSFCList(const EXP_PARAM *pexp, const int *PresentElems, const XrayYield *XYldArray, const AbsCoef *TotAbsCoefArray, int *dimList, SFCListElem **SFCList);

int needSFC(atomicnumber Za, atomicnumber Zb, const CalibYld *enerA, const AbsCoef *AbsC, int *dimList, SFCListElem **SFCList );

double sigmaphoto(int iabs, const AbsCoef *AbsC, atomicnumber Zb, double Xray);

int inrange(double value,double min, double max);

---