libCPIXE: libcpixe/libcpixe.h Source File

00001 
00002 /***************************************************************************
00003     Copyright (C) 2007 by Carlos Pascual-Izarra
00004     carlos.pascual_AT_users.sourceforge.net                                                  *
00005 
00006     This program is free software: you can redistribute it and/or modify
00007     it under the terms of the GNU General Public License as published by
00008     the Free Software Foundation, either version 3 of the License, or
00009     (at your option) any later version.
00010 
00011     This program is distributed in the hope that it will be useful,
00012     but WITHOUT ANY WARRANTY; without even the implied warranty of
00013     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014     GNU General Public License for more details.
00015 
00016     You should have received a copy of the GNU General Public License
00017     along with this program.  If not, see <http://www.gnu.org/licenses/>.
00018 
00019  ***************************************************************************/
00020 
00030 #include "compilopt.h"
00031 
00033 #define LONGSTRINGLENGTH 256
00034 #define FILENMLENGTH LONGSTRINGLENGTH
00035 
00037 #define DEG2RAD (M_PI/180.0)
00038 
00040 #define MEVAMU 931.52181
00041 #define UATMASS (1./1822.887)
00042 
00044 #define BARNTOM2  1.0e-28
00045 
00047 #define ELECTRONCHARGE_IN_UC 1.60217646e-13
00048 
00050 #define EFFIFILEVERSION1 1
00052 #define EFFIFILEVERSION2 2
00054 #define EFFIFILEVERSION3 3
00055 
00056 typedef char STATFILENAME[FILENMLENGTH] ; 
00057 
00058 typedef char ChemSymb[3];
00059 
00060 typedef int atomicnumber;
00061 
00062 //typedef double tenarray[11];
00063 
00064 //typedef double ninearray[10];
00065 
00066 //typedef double threearray[4];
00067 
00068 //typedef double LArray[4][4];
00069 
00070 //typedef char LineNmType[6][4][11];
00071 
00072 typedef double SecXL[4];
00073 
00074 typedef enum { K, L, M } FwTipo;
00075 
00076 /*Structures*/
00077 
00079 typedef struct ELEMENT
00080   {
00081   char symbol[3]; 
00082   int Z;          
00083   int A;          
00084   double M;       
00085   double IM;      
00086   }ELEMENT;
00087 
00089 typedef struct COMPOUND
00090   {
00091   char name[30];  
00092   double sumX;    
00093   int nelem;      
00094   ELEMENT *elem;  
00095   double *X;      
00096   double *xn;     
00097   double *w;      
00098   }COMPOUND;
00099 
00102 typedef struct SPT
00103   {
00104   char ID[30];
00105   int nrows;    
00106   int logmode;  
00107   double Emax;  
00108   double Estep; 
00109   double Smax;  
00110   double *E;    
00111   double *S;    
00112   double *dSdE; 
00113   }SPT;
00114 
00118 typedef struct absorber
00119   {
00120   int number;
00121   char name[6];
00122   int nreg;
00123   double ereg[11];
00124   double coef[11][3];
00125   double h[3];
00126   } absorber;
00127 
00128 
00130 typedef struct SIM_PARAM
00131   {
00132   int MaxZinsample;   
00133   int AllowSXFCorr;   
00134   int AllowXEqCalc;   
00135   double DTCC;        
00136   //double DTCCerr;     ///<Error in Dead Time Correction coefficient (Not really used)
00137   double ColCharge;   
00138   double CalEner;     
00139   int useFilter;      
00140   }SIM_PARAM;
00141 
00142 
00143 
00145 typedef struct EXP_PARAM
00146   {
00147   ELEMENT ion;      
00148   double ioncharge; 
00149   double BeamEner;  
00150   double BeamCol;   
00151   double DetColFac; 
00152   double IncAng;    
00153   double DetAng;    
00154   double cosDet;    
00155   double cosInc;    
00156   double cosFac;    
00157   double BeamCross; 
00158   double FinalEner; 
00159   double SAngFract; 
00160   double Fluence;   
00161   SIM_PARAM simpar; 
00162   }EXP_PARAM;
00163 
00165 typedef struct EXTRAINFO
00166   {
00167    char *SampleFileNm;  
00168    char *FilterFileNm;  
00169    char *DetectorFileNm;
00170    char *CalibFileNm;    
00171    char *AreasFileNm;   
00172    char *DBpath;        
00173    char *OutputFileNm;  
00174    int AreasFormat;     
00175    int WantOutputfile;  
00176    int givendetectorfile;
00177    int useefficiencies; 
00178   }EXTRAINFO;
00179 
00181 typedef struct PeakArea {
00182   ChemSymb symb;   
00183   atomicnumber atnum; 
00184   int area[5][3]; 
00185   int erro[5][3]; 
00186 } PeakArea;
00187 
00189 typedef struct CalibYld {
00190   double K_[4];   
00191   double L_[4][4];   
00192   double M_[4];   
00193 } CalibYld;
00194 
00197 typedef struct XrayYield {
00198   ChemSymb symb;
00199   atomicnumber atnum;
00200   CalibYld ener;
00201   CalibYld XYld;
00202 } XrayYield;
00203 
00205 typedef struct FluorCKCoef {
00206   atomicnumber atnum;
00207   double w[10];
00208   double ck[3];
00209   CalibYld k;
00210 } FluorCKCoef;
00211 
00213 typedef struct SRType {
00214   CalibYld AbsFact;
00215   CalibYld Xeq;
00216   CalibYld SFCr;
00217   CalibYld DetLim;
00218   CalibYld ConcErr;
00219   CalibYld MDV;
00220 } SRType;
00221 
00222 
00224 typedef struct SecResType {
00225   PeakArea Pk;
00226   SRType SR;
00227 } SecResType;
00228 
00229 
00230 
00232 typedef struct ESxType {
00233   double ep;
00234   double stpp;
00235   double x;
00236 } ESxType;
00237 
00239 typedef struct foil {
00240   double thick;   
00241   int nfoilelm;   
00242   COMPOUND comp;  
00243 } foil;
00244 
00245 
00247 typedef struct
00248   {
00249   int NumOfTrc;            
00250   int *TrcAtnum;           
00251   int dimTrans;            
00252   CalibYld *Trans;         
00253 //  PeakArea *AreasArray;  ///<Experimental Peak Areas. Not used for simulation. @todo check if it makes sense to maintain it here
00254 //  XrayYield *ResArray;   ///<@todo description?
00255   XrayYield *ResYldArray;  
00256   SecResType *SecResArray; 
00257   int *TrcUse;             
00258   int *NeedSFC;            
00259   ESxType *ESxArray;       
00260   int FESxlen;             
00261   CalibYld *SSTrsArray;    
00262   CalibYld *SSYldArray;    
00263   double ThickIn;          
00264   double FoilInEner;       
00265   double FoilOutEner;      
00266   double absolutepos;      
00267   const foil *pFoil;       
00268   } LYR;
00269 
00270 
00272 typedef struct AbsCoef {
00273   atomicnumber atnum;     
00274   double coefstd[23];     
00275   double enr[10];         
00276   double coefenr[9][2];   
00277 } AbsCoef;
00278 
00280 typedef struct FILTER {
00281   int nlyr;         
00282   int MaxZ;         
00283   int *FilterElems; 
00284   double geomcorr;  
00285   foil *foil;       
00286   int dimTrans;     
00287   CalibYld *Trans;  
00288   int changes;      
00289   int FilterType;   
00290   absorber Calibflt;
00291 }FILTER;
00292 
00294 typedef struct CPIXERESULTS{
00295   int atnum;          
00296   CalibYld simareas;  
00297   CalibYld err;       
00298 }CPIXERESULTS;
00299 
00301 typedef struct SFCListElem {
00302   atomicnumber Za;  
00303   atomicnumber Zb;  
00304   int epri;         
00305   int esec;         
00306   int abs;          
00307   double Ea;        
00308   double sigmaphoto;
00309 } SFCListElem;
00310 
00311 
00313 typedef struct IntCalibYld {
00314   int K_[4];   
00315   int L_[4][4];   
00316   int M_[4];   
00317 } IntCalibYld;
00318 
00320 typedef struct TwoCol {
00321   double x;   
00322   double y;   
00323 } TwoCol;
00324 
00325 
00326 // /** Structure for use in when using efficiency-based calibration.*/
00327 // typedef struct  {
00328 //   atomicnumber atnum;   ///<atomic number for the element
00329 //   CalibYld ener;        ///<Energy for each line
00330 //   CalibYld raweffi;     ///<raw efficiency (=num_of_detected_rays/num_of_rays_arriving_the_detector)
00331 //   CalibYld correffi;    ///<corrected efficiency (=raweffi/detector_transmission)
00332 //   IntCalibYld exp_flag; ///<Flag indicating wether effi was determined experimentally for this line (1) or not (0)
00333 // } CPIXECALIB;
00334 
00335 
00336 /*Constants and Initializators*/
00340 //static atomicnumber minK = 11, maxK = 54, minL = 30, maxL = 92, minM = 60;
00341 static const atomicnumber minK = 11, maxK = 54, minL = 30, minM = 60;
00342 
00343 static const CalibYld CYldNul = { { 0.0, 0.0, 0.0, 0.0 },
00344                             { { 0.0, 0.0, 0.0, 0.0 },
00345                               { 0.0, 0.0, 0.0, 0.0 },
00346                               { 0.0, 0.0, 0.0, 0.0 },
00347                               { 0.0, 0.0, 0.0, 0.0 }},
00348                             { 0.0, 0.0, 0.0, 0.0 }};
00349 
00352 //static LineNmType LineNm = {
00353 static char LineNm[6][4][11] = {
00354   { "--", "  --      ", "  --      ", "  --      " },
00355   { "--", "K_alpha1,2", "K_beta_1  ", "K_beta_2  " },
00356   { "--", "L_alpha1,2", "L_beta_2  ", "L_l       " },
00357   { "--", "L_beta_1  ", "L_gamma_1 ", "L_eta     " },
00358   { "--", "L_beta_3  ", "L_beta_4  ", "L_gamma_3 " },
00359   { "--", "M_alpha1,2", "M_beta    ", "M_gamma   " }
00360 };
00361 
00362 /*Chemical symbols. Z-ordered (e.g. ChemicalSymbol[14]="Si")*/
00363 static ChemSymb ChemicalSymbol[110] = {
00364     "--", "H" , "He", "Li", "Be", "B" , "C" , "N" , "O" , "F" , "Ne", "Na",
00365     "Mg", "Al", "Si", "P" , "S" , "Cl", "Ar", "K" , "Ca", "Sc", "Ti", "V" ,
00366     "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se", "Br",
00367     "Kr", "Rb", "Sr", "Y" , "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag",
00368     "Cd", "In", "Sn", "Sb", "Te", "I" , "Xe", "Cs", "Ba", "La", "Ce", "Pr",
00369     "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu",
00370     "Hf", "Ta", "W" , "Re", "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi",
00371     "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th", "Pa", "U" , "Np", "Pu", "Am",
00372     "Cm", "Bk", "Cf", "Es", "Fm", "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh",
00373     "Hs", "Mt"
00374   };
00375 
00376 /*Function declarations*/
00377 
00378 //void testlib(int **testarray,int *n);
00379 
00380 int symbol2Z(char *symbol);
00381 
00382 int Z2mass(int Z, double *mass, char option);
00383 
00384 int readCOMPOUND(FILE *f, int nelem, COMPOUND *c);
00385 
00386 int readINPUT(const char *InputFileNm, EXP_PARAM *pexppar, EXTRAINFO *pExtraInfo);
00387 
00388 // int readEXP_PARAM(char *ExpParFileNm, EXP_PARAM *pexppar);
00389 
00390 // int readSIM_PARAM(char *SampleParamFileNm, SIM_PARAM *simpar,char *IniMatFileNm, char *FilterFileNm, char *XYldFileNm);
00391 
00392 int readsample(char *SampleDefFileNm, int *MaxZ, int *NFoil, foil **Sample);
00393 
00394 int readFilter(const char *FilterDefFileNm, FILTER *Filter);
00395 
00396 int createPresentElems(int MaxZ, int NFoil,const foil *MatArray, int **PresentElems);
00397 
00398 int readCalcFlags(const char *CalcFlagsFileNm, const int *PresentElems, int MaxZinsample, int AreasFormat, CPIXERESULTS **CalcFlags);
00399 
00400 int readXYld(const char *XYldFileNm, const int *PresentElems, const CPIXERESULTS *CalcFlags, int MaxZinsample, double *CalEner, XrayYield **XYldArray);
00401 
00402 int readEff(const char *CalibFileNm, const int *PresentElems, int MaxZinsample, CalibYld **EffArray);
00403 
00404 int readLinesEner(const char *LinesEnerFileNm, const int *PresentElems, int MaxZinsample, CalibYld **LinesEnerArray);
00405 
00406 int readEff2(const char *CalibFileNm, const int *PresentElems, int MaxZinsample, const CalibYld *LinesEnerArray, CalibYld **EffArray);
00407 
00408 int readEff3(const char *CalibFileNm, const int *PresentElems, int MaxZinsample, const CalibYld *LinesEnerArray, CalibYld **EffArray);
00409 
00410 double interpolate_efficiency(int ndata, const double *Xarray, const double *Yarray, double X);
00411 
00412 int CreateEff(const char *TotAbsCoefFileNm, const  char *LinesEnerFileNm, const char *WindowDefFileNm, const char *CrystalDefFileNm, double factor ,CalibYld **EffBlockArray, int *dimEffTable, TwoCol **EffTable);
00413 
00414 int readFCK(char *FCKCoefFileNm, int MaxZinsample, FluorCKCoef **FCKCoefArray);
00415 
00416 //int readAbsCoef_OLD(char *TotAbsCoefFileNm, int MaxZinsample, const int *PresentElems, const FILTER *Filter, AbsCoef **TotAbsCoefArray);
00417 
00418 int readAbsCoef(const char *TotAbsCoefFileNm, int MaxZinsample, AbsCoef **TotAbsCoefArray);
00419 
00420 int createSPTs(const char *path, const EXP_PARAM *pexp, int MaxZinsample, const int *PresentElems, double step, SPT **SPTArray);
00421 
00422 double getSP(const COMPOUND *pcmp, const SPT *SPTArray, double E);
00423 
00424 //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);
00425 
00426 int initlyrarray(const EXP_PARAM *pexp, const foil *MatArray, const CalibYld *LinesEnerArray, const AbsCoef *TotAbsCoefArray, const SPT *SPTArray, const int *PresentElems, int NFoil, int *NFoilUsed, LYR **plyrarray);
00427 
00428 //int initlyr(const EXP_PARAM *pexp, const XrayYield *XYldArray, const AbsCoef *TotAbsCoefArray, const int *PresentElems, LYR *plyr, double *MACoef);
00429 
00430 int initlyr(const EXP_PARAM *pexp, const CalibYld *LinesEnerArray, const AbsCoef *TotAbsCoefArray, const int *PresentElems, LYR *plyr, double *MACoef);
00431 
00432 int AllFilterTrans(const char *TotAbsCoefFileNm, const  char *LinesEnerFileNm, const char *FilterDefFileNm, CalibYld *AFTArray);
00433 
00434 //int FilterTrans(const EXP_PARAM *pexp, const XrayYield *XYldArray, const AbsCoef *TotAbsCoefArray, const int *PresentElems, FILTER *Filter);
00435 int FilterTrans(int MaxZinsample, const CalibYld *LinesEner, const AbsCoef *TotAbsCoefArray, const int *PresentElems, FILTER *Filter);
00436 
00437 //int Transmission(const XrayYield *XYldArray, const AbsCoef *TotAbsCoefArray, int Z, const foil *pFoil, double geomcorr, const CalibYld *pTransOld, CalibYld *pTrans);
00438 int Transmission(const CalibYld *LinesEner,  const AbsCoef *TotAbsCoefArray, int Z, const foil *pFoil, double geomcorr, const CalibYld *pTransOld, CalibYld *pTrans);
00439 
00440 double TotAbsor(const AbsCoef *TotAbsCoefArray, const COMPOUND *cmp, double Xray);
00441 
00442 double TotAbsor_elemental(int iener, const AbsCoef *Absco, int Z, double Xray);
00443 
00444 int createsublyrs(const EXP_PARAM *pexp, const foil *pMat, const SPT *SPTArray, double MajAbsCoef, LYR *plyr);
00445 
00446 int SSThick(const EXP_PARAM *pexp, const foil *pMat, const SPT *SPTArray, double MajAbsCoef, double LayerThickness, ESxType ESxin,
00447         ESxType *ESxfin, int *pFpos, double *thick, ESxType **ESxA);
00448 
00449 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);
00450 
00451 int integrate_Simpson2(const EXP_PARAM *pexp, const AbsCoef *TotAbsCoefArray,
00452                        const FluorCKCoef *FCKCoefArray, const CalibYld *RawEffiArray,
00453                        int NFoilUsed, const FILTER *Filter, LYR *plyrarray, CalibYld *XYldSums);
00454 
00455 int Xprod(const AbsCoef *AbsC, const FluorCKCoef *pFCK, atomicnumber Z1, atomicnumber Z2, double M2 , double ener, CalibYld *XYld);
00456 
00457 double PaulX(double ener, atomicnumber z);
00458 
00459 double PaulX_y(double MeV, atomicnumber z);
00460 
00461 double PaulX_lp(double x, long p);
00462 
00463 int ReisX(const AbsCoef *AbsC, const FluorCKCoef *pFCK, double ener, atomicnumber z, double M2, double *sigmaXL);
00464 
00465 double ReisX_gs(FwTipo T, int SS, double ksis);
00466 
00467 double ReisX_hs(FwTipo T, int SS, double ksih, double thet);
00468 
00469 double ReisX_En(double z, int niu);
00470 
00471 double ReisX_polisec(int ssind, double kz, double tz);
00472 
00473 double ReisX_g(int ss, atomicnumber Zg, double xi);
00474 
00475 void PenInteg(atomicnumber atnumb, const CalibYld *absc, const ESxType *ESA,
00476         const CalibYld *YldA, const CalibYld *TrsA, const CalibYld *pTrs0,
00477         int FExlen, int NeedSFC, int AllowXEqCalc,
00478         double x0, double CosInc,
00479         CalibYld *XYld, CalibYld *XSFCr, CalibYld *XYldxmed);
00480 
00481 double Simps(double a, double b, double fa, double fi, double fb);
00482 
00483 void deNormalize(const EXP_PARAM *pexp, const AbsCoef *AbsC,
00484                const FluorCKCoef *pFCK, atomicnumber Z2, double M2,  double attfraction,
00485                const CalibYld *pXYld, XrayYield *ResY, CalibYld *XYldSum);
00486 
00487 void deNormalize2(const EXP_PARAM *pexp, atomicnumber Z2, double attfraction,
00488                const CalibYld *pEff, const CalibYld *PenInteg, CalibYld *XYld, CalibYld *XYldSum);
00489 
00490 void freeFilter(FILTER *Filter);
00491 
00492 void freeExtraInfo(EXTRAINFO *extrainfo);
00493 
00494 void freeReusable(int NFoils, LYR **plyrarray, foil **MatArray, CalibYld **XYldSums );
00495 
00496 void safefree(void **ptr);
00497 
00498 void fprintCALIBYLD(FILE *f,const CalibYld *pCYld);
00499 
00500 void fscanCALIBYLD(FILE *f, CalibYld *pCYld);
00501 
00502 int createSFCList(const EXP_PARAM *pexp, const int *PresentElems, const XrayYield *XYldArray, const AbsCoef *TotAbsCoefArray, int *dimList, SFCListElem **SFCList);
00503 
00504 int needSFC(atomicnumber Za, atomicnumber Zb, const CalibYld *enerA, const AbsCoef *AbsC, int *dimList, SFCListElem **SFCList );
00505 
00506 double sigmaphoto(int iabs, const AbsCoef *AbsC, atomicnumber Zb, double Xray);
00507 
00508 int inrange(double value,double min, double max);
00509 
00510 int compare_Twocol_x (const TwoCol *a, const TwoCol *b);
00511