libcpixe.c

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.c
LibCPIXE is a collection of C functions which are useful for Particle
Induced X-ray Emmission (PIXE) data analysis. The main objective is to produce a
publicly available set of functions which can be incorporated in IBA data
analysis codes. 

News, documentation and latest releases can be found at: http://cpixe.sourceforge.net 

A number of papers will be published and should be cited if libcpixe is used for producing scientific publications.

The reference paper for LibCPIXE and CPIXE is:
"LibCPIXE: an open-source library for PIXE simulation and Analysis", Carlos Pascual-Izarra, Nuno P. Barradas and Miguel A. Reis, To be published in Proceedings of the IBA2005 (NIMB)

A paper showing the integration of LibCPIXE into the Data Furnace code is:
"Simultaneous PIXE and RBS analysis using Bayesian Inference" Carlos Pascual-Izarra, Miguel A. Reis and Nuno P. Barradas, To be published in Proceedings of the IBA2005 (NIMB)

The seed for this code is on a translation from Pascal to C of
parts of the code writen by M.A. Reis for the DatPixe-v5.3 code. The Pascal code
was initially translated using p2c (an automatic Pascal-to-C free translator)
and afterwards heavily modified.
*/

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




/*Global Constants*/

/*Physical constants of interest*/
// double m0c2 = 1.503186433e-0010, keV = 1.60621e-0016, uAtmass = 1 / 1822.887,
//        M1sme = 1836.076012, m0pamu = 1.007276470, m0pMeV = 938.3,
//        MeVpamu = 931.52181, v0 = 2.187578544e6, v02 = 4.785499886e12,
//        c2 = 8.98740441e16, cau = 137.03604, pi = 3.141592654,
//        a02 = 2.800283608e-21, Ry = 13.6, Hr = 27.2116, barn = 1.0e-28;



/*Functions*/


/**Reads an input file consisting on a series of commands.
 * 
 * @param InputFileNm (I) Name of the input file (no more than LONGSTRINGLENGTH characters)
 * @param pexppar (O) pointer to the experimental parameter structure
 * @param pExtraInfo (O) pointer to structure containing misc data read from input file 
 * @return Always 0
 */
int readINPUT(const char *InputFileNm, EXP_PARAM *pexppar, EXTRAINFO *pExtraInfo)
{
  FILE *f;
  int i;
  char *command,*arg;
  char line[LONGSTRINGLENGTH]="";
  int flag[16]; /*We use 11 control flags to check that the reading was fine (they must all end with a value different than 0)*/
  
  
  /*Reset control flags*/
  for(i=0;i<16;i++)flag[i]=0;
  pExtraInfo->WantOutputfile=0;
  pExtraInfo->DoCalibration=0;
  //set to -1 (which is true) some flags regarding to non-mandatory commands (and set default value for those commands); 
  pexppar->simpar.AllowSXFCorr=0;flag[10]=-1; 
  pexppar->simpar.AllowXEqCalc=0; flag[11]=-1;
  
  f = fopen(InputFileNm, "r");
  if (f == NULL)
    {fprintf(stderr,"\nERROR: Could not open file '%s' for read.\n",InputFileNm);exit(1);}
  
  /*Skip until "BEGIN_INPUT"*/
  for(;!feof(f) && strcasecmp(line, "BEGIN_INPUT");){
   fscanf(f,"%255s",line);
   flag[0]=1;
  }  
  
  for(;!feof(f);){
    //reads an entire line of at much LONGSTRINGLENGTH characters
    fgets(line,sizeof(line),f);   
    
    //Now we parse the line looking for a known command and discarding comments
    //fscanf(f,"%255s",command);
    command = strtok( line, "\n\t \r" ) ;
    if( command != NULL && command[0] != '#' ){
      
      
      
      /*"END_INPUT" command just exits the for loop*/
      if(!strcmp(command,"END_INPUT")){
        flag[1]=1;
        break;
      }
      else if(!strcmp(command,"ION")){
        /*Syntax: 
        ION symbol   <---"symbol" is the beam ion chem. symbol. The Most abundant isotope will be used by default. */
//         fscanf(f,"%3s",pexppar->ion.symbol);
        arg=strtok( NULL,"\n\t \r");
        strncpy(pexppar->ion.symbol,arg,3);
        ///@todo For the moment, only H ions are supported
        if (strncasecmp(pexppar->ion.symbol,"H",3)) 
          {fprintf(stderr,"\nERROR: Currently, only H ions are supported\n");exit(1);} 
        pexppar->ion.Z=symbol2Z(pexppar->ion.symbol);
        pexppar->ion.A=Z2mass(pexppar->ion.Z, &pexppar->ion.IM,'m');
        Z2mass(pexppar->ion.Z, &pexppar->ion.M,'n');
        flag[2]=1;
        //printf("\n\n!!!!!!!!!!!!!!!!!!!!!!!! ION='%s'\n\n",pexppar->ion.symbol);getchar();
      }
      else if(!strcmp(command,"IONMASS")){
        /*Syntax: 
        IONMASS mass   <---"mass" is the beam ion weight in amu. This can be used for defining isotopes. */
//         fscanf(f,"%lf", &pexppar->ion.M);
        arg=strtok( NULL,"\n\t \r");
        pexppar->ion.M=atof(arg);
        /*Note that this is an optional command, no flag is changed*/
      }
      else if(!strcmp(command,"KEV")){
        /*Syntax: 
        KEV Ebeam   <---"Ebeam" is Beam energy in keV.  */
//         fscanf(f, "%lg", &pexppar->BeamEner); //keV
        arg=strtok( NULL,"\n\t \r");
        pexppar->BeamEner=atof(arg);
        flag[3]=1;
      }
      else if(!strcmp(command,"INCANG")){
        /*Syntax: 
        INCANG thetavalue   <---Incience angle: Sample normal to beam*/
//         fscanf(f, "%lg", &pexppar->IncAng);
        arg=strtok( NULL,"\n\t \r");
        pexppar->IncAng=atof(arg);
        pexppar->IncAng *= DEG2RAD;
        pexppar->cosInc = cos(pexppar->IncAng);
        flag[4]=1;
        }
      else if(!strcmp(command,"EXITANG")){ 
        /*Syntax: 
        EXITANG phivalue   <---Exit angle. Sample normal to detector. */
//         fscanf(f, "%lg", &pexppar->DetAng); 
        arg=strtok( NULL,"\n\t \r");
        pexppar->DetAng=atof(arg);
        pexppar->DetAng *= DEG2RAD;
        pexppar->cosDet = cos(pexppar->DetAng);
        flag[5]=1;
      }
      else if(!strcmp(command,"BEAMCOL")){ 
        /*Syntax: 
        BEAMCOL diameter   <--Beam diameter, in mm */
//         fscanf(f, "%lg", &pexppar->BeamCol);
        arg=strtok( NULL,"\n\t \r");
        pexppar->BeamCol=atof(arg);
        pexppar->BeamCross = M_PI * (pexppar->BeamCol*pexppar->BeamCol) / 400;   /* beam cross-section in cm2 */
        flag[6]=1;
      }
      else if(!strcmp(command,"DETCOLFAC")){ 
        /*Syntax: 
        DETCOLFAC factor   <---Detector colimator factor... really this should be solid angle */
        arg=strtok( NULL,"\n\t \r");
        pexppar->DetColFac=atof(arg);
//         fscanf(f, "%lg", &pexppar->DetColFac);  
        flag[7]=1;
      }
      else if(!strcmp(command,"CHARGE")){ 
        /*Syntax: 
        CHARGE ColCharge   <---Collected charge, in uC */
        arg=strtok( NULL,"\n\t \r");
        pexppar->simpar.ColCharge=atof(arg);
//         fscanf(f, "%lg", &pexppar->simpar.ColCharge);  
        flag[8]=1;
      }
      else if(!strcmp(command,"LIVETIME")){ 
        /*Syntax: 
        LIVETIME factor   <---Dead time correction (Integrated"Charge"/IntegratedCounts) */
        arg=strtok( NULL,"\n\t \r");
        pexppar->simpar.DTCC=atof(arg);
//         fscanf(f, "%lg", &pexppar->simpar.DTCC);  
        flag[9]=1;
      }
      else if(!strcmp(command,"ALLOWSXFCORR")){ 
        /*Syntax: 
        ALLOWSXFCORR flag   <---Flag allowing for calculation of Sec. Fluorescence  */  
        arg=strtok( NULL,"\n\t \r"); ///@todo: SXF not supported yet
        pexppar->simpar.AllowSXFCorr=atoi(arg);
//        fscanf(f, "%d", &pexppar->simpar.AllowSXFCorr);   
        flag[10]=1;
      }
      else if(!strcmp(command,"ALLOWXEQCALC")){ 
        /*Syntax: 
        ALLOWXEQCALC flag   <---Deat time correction (Integrated"Charge"/IntegratedCounts) */
        arg=strtok( NULL,"\n\t \r"); ///@todo: Xeq calculation not
        pexppar->simpar.AllowXEqCalc=atoi(arg);
//        fscanf(f, "%d", &pexppar->simpar.AllowXEqCalc);   
        flag[11]=1;
      }
      else if(!strcmp(command,"CALIBFILE")){
        /*Syntax:
        CALIBFILE Calfilename     <---Name of the calibration  (XrayYld) file
        */
        arg=strtok( NULL,"\n\t \r");
        strncpy(pExtraInfo->XYldFileNm,arg,LONGSTRINGLENGTH);
//         fscanf(f,"%255s",pExtraInfo->XYldFileNm);
        flag[12]=1;
        pExtraInfo->DoCalibration=0;
      }  
      else if(!strcmp(command,"DOCALIBRATION")){
        /*Syntax:
        DOCALIBRATION flag    <---Flag indicating whether the intention is to simulate (0) or to do a calibration (1)
        */
        arg=strtok( NULL,"\n\t \r");
        pExtraInfo->DoCalibration=atoi(arg);
//         fscanf(f, "%d", &pExtraInfo->DoCalibration);
        flag[12]=1;
      } 
      else if(!strcmp(command,"SAMPLEFILE")){
        /*Syntax:
        SAMPLEFILE SampleFilename     <---Name of the Sample definition file
        */
//         fscanf(f,"%255s",pExtraInfo->SampleFileNm);
        arg=strtok( NULL,"\n\t \r");
        strncpy(pExtraInfo->SampleFileNm,arg,LONGSTRINGLENGTH);
        flag[13]=1;
      } 
      else if(!strcmp(command,"FILTERFILE")){
        /*Syntax:
        FILTERFILE Filterfilename     <---Name of the Filter definition file
        */
//         fscanf(f,"%255s",pExtraInfo->FilterFileNm);
        arg=strtok( NULL,"\n\t \r");
        strncpy(pExtraInfo->FilterFileNm,arg,LONGSTRINGLENGTH);
        flag[14]=1;
      }  
      else if(!strcmp(command,"CFLAGSFILE")){
        /*Syntax:
        CFLAGSFILE CFlagsfilename     <---Name of the Calculation flags definition file
        */
//         fscanf(f,"%255s",pExtraInfo->AreasFileNm);
        arg=strtok( NULL,"\n\t \r");
        strncpy(pExtraInfo->AreasFileNm,arg,LONGSTRINGLENGTH);
        pExtraInfo->AreasFormat=1;
        flag[15]=1;
      }
      else if(!strcmp(command,"DTAREASFILE")){
        /*Syntax:
        DTAREASFILE DTAreasfilename     <---Name of the Areas file (in DATTPIXE format)
        */
//         fscanf(f,"%255s",pExtraInfo->AreasFileNm);
        arg=strtok( NULL,"\n\t \r");
        strncpy(pExtraInfo->AreasFileNm,arg,LONGSTRINGLENGTH);
        pExtraInfo->AreasFormat=2;
        flag[15]=1;
      } 
      else if(!strcmp(command,"DBPATH")){
        /*Syntax:
        DBPATH Path   <---Path to DataBase files (optional)
        */
//         fscanf(f,"%255s",pExtraInfo->DBpath);
        arg=strtok( NULL,"\n\t \r");
        strncpy(pExtraInfo->DBpath,arg,LONGSTRINGLENGTH);
      }  
      else if(!strcmp(command,"OUTPUTFILE")){
        /*Syntax:
        OUTPUTFILE Outputfilename   <---Name for the output file
        */
//         fscanf(f,"%255s",pExtraInfo->OutputFileNm);
        arg=strtok( NULL,"\n\t \r");
        strncpy(pExtraInfo->OutputFileNm,arg,LONGSTRINGLENGTH);
        pExtraInfo->WantOutputfile=1;
        /*note: this command is optional, but a flag indicating if it was used is passed*/
      }  
      /*If "command" is not valid...*/
      else {
        fprintf(stderr,"\n ERROR:ReadInput: command '%s' not known. Press <INTRO> to continue ...***\n",command);
        getchar();
      }
    }
  }
  
  //close input file
  fclose(f);
  
  /*check control flags*/  
  for(i=0;i<16;i++)if(flag[i]==0){fprintf(stderr,"\n ERROR:ReadInput: Missing command in inputfile (#%i)\n",i);exit(1);}
  
  pexppar->cosFac = pexppar->cosInc / pexppar->cosDet;
  
  pexppar->FinalEner = 0.1 * pexppar->BeamEner;  //below this energy no calculations will be done
  if(pexppar->FinalEner>100.)pexppar->FinalEner = 100.; 
  
  return(0); 
  
}




/**Fills an EXP_PARAM struct with data from a file.
  The file must have the following format:
  ------------------
  BeamEner
  BeamCol
  DetColFac
  IncAng
  DetAng
  ------------------

  * 
  * @param ExpParFileNm (I) Name of the file to be read.
  * @param pexppar (O) Pointer to struct containing the experimental parameters
  * @return Always 0
  */
 int readEXP_PARAM(char *ExpParFileNm, EXP_PARAM *pexppar)
{
  //EXP_PARAM pexppar;
  FILE *ExpParFile;
    
  ExpParFile = fopen(ExpParFileNm, "r");
  if (ExpParFile == NULL)
    {fprintf(stderr,"\nERROR: Could not open file '%s' for read.\n",ExpParFileNm);exit(1);}
  
  
  strncpy(pexppar->ion.symbol,"H",2);  ///@todo For the moment, only H ions are supported
  pexppar->ion.Z=symbol2Z(pexppar->ion.symbol);
  pexppar->ion.A=Z2mass(pexppar->ion.Z, &pexppar->ion.IM,'m');
  Z2mass(pexppar->ion.Z, &pexppar->ion.M,'n');
 
  fscanf(ExpParFile, "%lg", &pexppar->BeamEner); //keV
  fscanf(ExpParFile, "%lg", &pexppar->BeamCol);  
  fscanf(ExpParFile, "%lg", &pexppar->DetColFac); 
  fscanf(ExpParFile, "%lg", &pexppar->IncAng);  
  fscanf(ExpParFile, "%lg", &pexppar->DetAng);   
  
  fclose(ExpParFile);
  
  pexppar->IncAng *= DEG2RAD;
  pexppar->cosInc = cos(pexppar->IncAng);
  pexppar->DetAng *= DEG2RAD;
  pexppar->cosDet = cos(pexppar->DetAng);
  pexppar->BeamCross = M_PI * (pexppar->BeamCol*pexppar->BeamCol) / 400;   /* beam cross-section in cm2 */
  

  pexppar->cosFac = pexppar->cosInc / pexppar->cosDet;
  
  pexppar->FinalEner = 0.1 * pexppar->BeamEner;  //below this energy no calculations will be done
  if(pexppar->FinalEner>100.)pexppar->FinalEner = 100.; 

      
  return(0);
}

/**Fills a SIM_PARAM struct with data from a file.
  The file must have the following format:
  ------------------
  Samplefilename
  absorbertype
  AllowSXFCorr
  AllowXEqCalc
  GenPulNumb
  GenInteg
  ColCharge
  ------------------
 * 
 * @param SampleParamFileNm (I) Name of the file to be read.
 * @param simpar (O) Pointer to struct containing the simulation parameters
 * @param IniMatFileNm (O) Name of the file defining sample
 * @param FilterFileNm (O) Name of the filter definition file
 * @param XYldFileNm (O) Name of the calibration yield file
 * @return Always 0
 */
int readSIM_PARAM(char *SampleParamFileNm, SIM_PARAM *simpar,char *IniMatFileNm, char *FilterFileNm, char *XYldFileNm)
{
  
  FILE *SampleParamFile;
  char AbsType[6];
//   double GenPulNumb,GenInteg;
//   absorber FilterNul ={ 0, "00000", 1,
//                       { 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, 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 }  };
  
    
    
  SampleParamFile = fopen(SampleParamFileNm, "r");
  if (SampleParamFile == NULL)
    {fprintf(stderr,"\nERROR: Could not open file '%s' for read.\n",SampleParamFileNm);exit(1);}

 
  fscanf(SampleParamFile,"%255s",IniMatFileNm);  //name of sample file
  fscanf(SampleParamFile,"%255s",FilterFileNm);  //name of Filter file
  fscanf(SampleParamFile,"%255s",XYldFileNm);  //name of calibration file
  fscanf(SampleParamFile, "%5s", AbsType);   // is Char[6]
  fscanf(SampleParamFile, "%d",  &simpar->AllowSXFCorr);   //1=true, 0=false
  fscanf(SampleParamFile, "%d",  &simpar->AllowXEqCalc);   //1=true, 0=false
  fscanf(SampleParamFile, "%lg", &simpar->DTCC); //deadtime correction (Integrated"Charge"/IntegratedCounts)
//   fscanf(SampleParamFile, "%lg", &GenPulNumb);
//   fscanf(SampleParamFile, "%lg", &GenInteg);
//   simpar->DTCC = GenInteg / GenPulNumb;
  //simpar->DTCCerr = 2 * sqrt(GenPulNumb - GenInteg) / GenPulNumb;
  fscanf(SampleParamFile, "%lg", &simpar->ColCharge); //Colected charge, in ?? uC ??
    
  fclose(SampleParamFile);
  
  ///@todo complete here the initialization for the case when a filter is used. UPDATE: this should be managed elsewhere, probably in FilterTrans() or readFilter()
  if (strcasecmp(AbsType, "NoAbs")) {
    simpar->useFilter=1;
    /*
    strcpy(FiltersCoefFile.name, FiltersCoefFileNm);
    if (*FiltersCoefFile.name != '\0') {
      if (FiltersCoefFile.f != NULL)
  FiltersCoefFile.f = freopen(FiltersCoefFile.name, "rb",
            FiltersCoefFile.f);
      else
  FiltersCoefFile.f = fopen(FiltersCoefFile.name, "rb");
    } else
      rewind(FiltersCoefFile.f);
    if (FiltersCoefFile.f == NULL)
      _EscIO2(FileNotFound, FiltersCoefFile.name);
    RESETBUF(FiltersCoefFile.f, absorber);
    fseek(FiltersCoefFile.f, (AbsN - 1) * sizeof(absorber), SEEK_SET);
    SETUPBUF(FiltersCoefFile.f, absorber);
    fread(&Filter, sizeof(absorber), 1, FiltersCoefFile.f);
    fclose(FiltersCoefFile.f);
    FiltersCoefFile.f = NULL;
    strcpy(AbsType, Filter.name);
    */
  }else {
//     simpar->Filter = FilterNul;
    simpar->useFilter=0;
  }
  return(0);  
}




/**Returns the atomic number for a certain element (up to Z=109)
 * This is a faster version than that of hstoplib.c because it does not access any file.
 * 
 * @param symbol (I) Chemical symbol of the element (case insensitive)
 * @return Atomic Number
 */
int symbol2Z(char *symbol){
  int Z;
  /**Array containing the chemical symbols from H(Z=1) to Mt(Z=109) */
  ChemSymb chsymbols[110] = {
    "--", "H" , "He", "Li", "Be", "B" , "C" , "N" , "O" , "F" , "Ne", "Na",
    "Mg", "Al", "Si", "P" , "S" , "Cl", "Ar", "K" , "Ca", "Sc", "Ti", "V" , 
    "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se", "Br", 
    "Kr", "Rb", "Sr", "Y" , "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", 
    "Cd", "In", "Sn", "Sb", "Te", "I" , "Xe", "Cs", "Ba", "La", "Ce", "Pr", 
    "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", 
    "Hf", "Ta", "W" , "Re", "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi", 
    "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th", "Pa", "U" , "Np", "Pu", "Am", 
    "Cm", "Bk", "Cf", "Es", "Fm", "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", 
    "Hs", "Mt"
  };
  for(Z=1;strncasecmp(symbol,chsymbols[Z],3) && Z<=99;Z++);
  if(Z<1 || Z> 109) {fprintf(stderr,"\nERROR: Chemical Symbol '%s' not known.\n",symbol);exit(1);}
  else return(Z);
}


/**Returns mass information for a given element (specified by its atomic number).
 * Currently the range is limited from Z=1 (H) to Z=92 (U).
 * It provides the Natural (average) weight, Most Abundant Isotope (MAI) weight and
 * Most Abundant Isotope number of nucleons
 * 
 * @param Z (I) Atomic number
 * @param mass (O) mass in amu. See also the "option" argument.
 * @param option (I) Either 'n' or 'm' to select "Natural" or "MAI" mass, respectively
 * @return Number of nucleons of the Most Abundant Isotope
 */
int Z2mass(int Z, double *mass, char option)
{
  /**Array containing the chemical symbols from H(Z=1) to Es(Z=99) */
  int maiAArray[93]={
      0 ,       1 ,       4 ,       7 ,       9 ,      11 ,      12 ,      14 ,      
     16 ,      19 ,      20 ,      23 ,      24 ,      27 ,      28 ,      31 ,      
     32 ,      35 ,      40 ,      39 ,      40 ,      45 ,      48 ,      51 ,      
     52 ,      55 ,      56 ,      59 ,      58 ,      63 ,      64 ,      69 ,      
     74 ,      75 ,      80 ,      79 ,      84 ,      85 ,      88 ,      89 ,      
     90 ,      93 ,      98 ,      97 ,     102 ,     103 ,     106 ,     107 ,     
    114 ,     115 ,     120 ,     121 ,     130 ,     127 ,     132 ,     133 ,     
    138 ,     139 ,     140 ,     141 ,     142 ,     148 ,     152 ,     153 ,     
    158 ,     159 ,     164 ,     165 ,     166 ,     169 ,     174 ,     175 ,     
    180 ,     181 ,     184 ,     187 ,     192 ,     193 ,     195 ,     197 ,     
    202 ,     205 ,     208 ,     209 ,     209 ,     210 ,     222 ,     223 ,     
    226 ,     227 ,     232 ,     231 ,     238};
  
  double maiMArray[93]={
     0.000  ,   1.008 ,   4.003 ,   7.016 ,   9.012 ,  11.009 ,  12.000 ,  14.003 ,  
     15.995 ,  18.998 ,  19.992 ,  22.990 ,  23.985 ,  26.982 ,  27.977 ,  30.974 ,  
     31.972 ,  34.969 ,  39.962 ,  38.964 ,  39.963 ,  44.956 ,  47.950 ,  50.940 ,  
     51.940 ,  54.940 ,  55.935 ,  58.930 ,  57.940 ,  62.930 ,  63.930 ,  68.930 ,  
     73.920 ,  74.920 ,  79.920 ,  78.920 ,  83.912 ,  84.910 ,  87.910 ,  88.906 ,  
     89.900 ,  92.910 ,  97.905 ,  97.000 , 101.900 , 102.900 , 105.900 , 106.900 , 
    113.900 , 114.900 , 119.900 , 120.900 , 129.906 , 126.900 , 131.904 , 132.905 , 
    137.905 , 139.000 , 140.000 , 141.000 , 142.000 , 148.000 , 152.000 , 153.000 , 
    157.900 , 158.925 , 164.000 , 165.000 , 166.000 , 169.000 , 174.000 , 175.000 , 
    180.000 , 181.000 , 184.000 , 187.000 , 192.000 , 193.000 , 195.000 , 197.000 , 
    202.000 , 205.000 , 208.000 , 209.000 , 208.982 , 210.000 , 222.000 , 223.000 , 
    226.000 , 227.000 , 232.000 , 231.000 , 238.040 };
  
  double natMArray[93]={
      0.000 ,   1.008 ,   4.003 ,   6.941 ,   9.012 ,  10.811 ,  12.011 ,  14.007 ,
     15.999 ,  18.998 ,  20.180 ,  22.990 ,  24.305 ,  26.982 ,  28.086 ,  30.974 ,  
     32.066 ,  35.453 ,  39.948 ,  39.098 ,  40.080 ,  44.956 ,  47.900 ,  50.942 ,  
     51.996 ,  54.938 ,  55.847 ,  58.933 ,  58.690 ,  63.546 ,  65.390 ,  69.720 ,  
     72.610 ,  74.922 ,  78.960 ,  79.904 ,  83.800 ,  85.470 ,  87.620 ,  88.905 ,  
     91.220 ,  92.906 ,  95.940 ,  97.000 , 101.070 , 102.910 , 106.400 , 107.870 , 
    112.400 , 114.820 , 118.710 , 121.750 , 127.600 , 126.900 , 131.300 , 132.910 , 
    137.327 , 138.910 , 140.120 , 140.910 , 144.240 , 148.000 , 150.360 , 151.970 , 
    157.250 , 158.930 , 162.500 , 164.930 , 167.260 , 168.930 , 173.040 , 174.970 , 
    178.490 , 180.950 , 183.850 , 186.200 , 190.200 , 192.200 , 195.080 , 196.970 , 
    200.590 , 204.380 , 207.190 , 208.980 , 210.000 , 210.000 , 222.000 , 223.000 , 
    226.000 , 227.000 , 232.000 , 231.000 , 238.030 };
  
  if(Z<1 || Z> 92) {fprintf(stderr,"\nERROR: Data for Atomic number '%d' not available.\n",Z);exit(1);}
  
  switch(option){
    case 'n':  // Natural Mass (i.e., isotopic average) , in amu
      *mass=natMArray[Z];
      break;
    case 'm':  // Most Abundant Isotope mass, in amu
      *mass=maiMArray[Z];
      break;
    default:
      {fprintf(stderr,"\nERROR: Z2mass() Option argument can only be 'n' (natural) or 'm' (most abundant) \n");exit(1);}
      break;
  }
  
  return(maiAArray[Z]);
}


/**Reads a compound of 'nelem' elements from the stream f and fills with this info  the the compound c
 * This is a slightly modified version from that of hstoplib.c 
 * It fills the normalized atomic and mass concentration.
 * Mass of the elements is obtained by calling Z2mass()
 *
 * The data from the stream must have the following format (repeated nelem times):
 * [element] [atomic_concentration]
 * ...
 * 
 * @param f (I) Pointer to a stream from which the compound definition is read.
 * @param nelem (I) Number of elements in the compound
 * @param c (O) Pointer to structure defining the compound
 * @return Maximum Atomic number read 
 */
int readCOMPOUND(FILE *f, int nelem, COMPOUND *c)
{
  int i,nchar1,nchar2,maxZ;
  double sumM;
  char temp[30],temp2[30];
  
  //if(nelem==0)readCOMPOUND2(f, c);  //Modified from the original
  if(nelem<1){fprintf(stderr,"\n ERROR: readCOMPOUND must read at least 1 element\n");exit(1);}
  else{
    c->nelem=nelem;
    if(!(c->elem=(ELEMENT*)calloc(nelem,sizeof(ELEMENT)))){fprintf(stderr,"\n Error allocating memory\n");exit(1);}
    if(!(c->X=(double*)calloc(nelem,sizeof(double)))){fprintf(stderr,"\n Error allocating memory\n");exit(1);}
    if(!(c->xn=(double*)calloc(nelem,sizeof(double)))){fprintf(stderr,"\n Error allocating memory\n");exit(1);}
    if(!(c->w=(double*)calloc(nelem,sizeof(double)))){fprintf(stderr,"\n Error allocating memory\n");exit(1);}
    for(maxZ=0.,c->sumX=0.,sumM=0.,i=0,nchar1=0,nchar2=0 ; i<nelem ; i++){
      fscanf(f,"%2s%lf",(c->elem[i].symbol),&(c->X[i]));
      c->elem[i].Z=symbol2Z(c->elem[i].symbol);
      c->elem[i].A=Z2mass(c->elem[i].Z, &c->elem[i].IM, 'm');
      Z2mass(c->elem[i].Z, &c->elem[i].M, 'n');
      
      
      if(c->elem[i].Z==0){fprintf(stderr,"\n ERROR: Chemical symbol '%s' not valid\n",c->elem[i].symbol);exit(1);}
      c->sumX+=c->X[i];
      sumM+=c->elem[i].M*c->X[i];
      if(c->elem[i].Z>maxZ)maxZ=c->elem[i].Z;
      nchar1+=strlen(c->elem[i].symbol);
      nchar2+=snprintf(temp,30,"%4lf",c->X[i]);
    }
    /*Normalize concentrations and calculate the MASS concentration*/
    for(i=0;i<nelem;i++){
      c->xn[i]=c->X[i]/c->sumX;
      c->w[i]=c->X[i]*c->elem[i].M/sumM;
    }
    /*Build the name of the compound*/
    strcpy(temp2,"");
    strcpy(temp,"");
    if(nchar1+nchar2 < 30 && nelem<4){
      for(i=0;i<nelem;i++){
        if(c->X[i]!=1.)snprintf(temp,30,"%s%1.lf",c->elem[i].symbol,c->X[i]);
        else snprintf(temp,30,"%s",c->elem[i].symbol);
        strcat(temp2,temp);
      }
    }
    else if(nchar1<30){
      for(i=0;i<nelem;i++){
        snprintf(temp,30,"%s",c->elem[i].symbol);
        strcat(temp2,temp);
      }
    }
    else{
      for(i=0;i<nelem;i++){
        snprintf(temp2,25,"%s%s",temp,c->elem[i].symbol);
        strcpy(temp,temp2);
      }
      snprintf(temp2,30,"%s+",temp); 
    } 
    strcpy(c->name,temp2);
  }
  return(maxZ);
}




/**Allocates and initializes the sample definition from a file. 
   It also returns Number of foils and Maximum atomic number present in sample.
  
  The File must have the following format:
  -----------------
  NUMBER_OF_FOILS [nfoils]
  FOIL [thickness_in_ug/cm2] 
  [numberofelems]
  [symbol] [concentration]
  [symbol] [concentration]
  ...
  FOIL CalEner
  [numberofelems]
  [symbol] [concentration]
  ...
  -----------------
  
  Note:Any number of foils and any number of elements are allowed.

 * 
 * @param SampleDefFileNm (I) Name of the file to be read.
 * @param MaxZ (O) Largest atomic number present in sample (Dimension of PresentElems).
 * @param NFoil (O) Number of target foils
 * @param Sample (O) Pointer to Array of files (Dim=NFoil)
 * @return Always 0
 */
int readsample(char *SampleDefFileNm, int *MaxZ, int *NFoil, foil **Sample)
{
  
  FILE *SampleDefFile;
  char dummy[LONGSTRINGLENGTH];
  int i,tempZ,j;
  foil *pfoil;
  double ug2at,natmass;
  SampleDefFile = fopen(SampleDefFileNm, "r");
  if (SampleDefFile == NULL)
    {fprintf(stderr,"\nERROR: Could not open file '%s' for read.\n",SampleDefFileNm);exit(1);}
      
  /*Skip until "NUMBER OF FOILS"*/
  do {
    fscanf(SampleDefFile,"%255s",dummy);   
  } while (strcasecmp(dummy, "NUMBER_OF_FOILS"));  
  
  fscanf(SampleDefFile, "%d", NFoil);
  
  /*Initialize  Sample*/
  //if(*Sample!=NULL)free(*Sample);
  *Sample=(foil*)malloc((*NFoil)*sizeof(foil));
  if (*Sample==NULL){fprintf(stderr,"\n Error allocating memory (Sample)\n");exit(1);}
  
  for (*MaxZ=0,i = 0; i < *NFoil ; i++) {
    pfoil=&((*Sample)[i]);
    do {
      fscanf(SampleDefFile,"%255s",dummy);   
      } while (strcasecmp(dummy, "FOIL"));  
    
    fscanf(SampleDefFile, "%lg %10s", &pfoil->thick, dummy);
    if(strcasecmp(dummy,"cm2") && strcasecmp(dummy,"mg")){
      fprintf(stderr,"\n Error: Bad syntax in sample definition file: '%s' is not a known thickness unit\n",dummy);exit(1);
    }
    fscanf(SampleDefFile, "%d",  &pfoil->nfoilelm);
    tempZ=readCOMPOUND(SampleDefFile, pfoil->nfoilelm, &pfoil->comp);
    if(tempZ>*MaxZ) *MaxZ=tempZ;
    
    if(!strcasecmp(dummy,"mg")){
      /*Convert the value in pfoil->thick from mg/cm2 to 1e15at/cm2 */
      for(natmass=0, j=0 ; j<pfoil->nfoilelm ; j++) natmass += pfoil->comp.xn[j] * pfoil->comp.elem[j].M;
      ug2at=natmass/602.2141; // (Pm [g/mol] *1e6 [ug/g] / (Navo [at/mol])) * 1e15  
                              // ==> ug2at is in [ug/1e15at]
      pfoil->thick *= (1e3/ug2at);  //Before this, thick was in mg/cm2, hence multiply by 1e3 to get ug/cm2 
                                     //and then divide by ug2at to get 1e15at/cm2
    }
    
  }
  fclose(SampleDefFile);
 
  return(0);
}

/**Creates an array of integer flags indicating whether an element
 * is present in the sample (1) or not (0). 
 * The dimension of the array is that of the maximum atomic number present in the sample. 
 * For example, if the sample contains only C and Si, the array will be:
 * PresentElems[6]= PresentElems[14]=1 (and all the other=0)
 * 
 * @param MaxZ (I) Largest atomic number present in sample (Dimension of PresentElems).
 * @param NFoil (I) Number of target foils
 * @param MatArray (I) Array of foils defining the current sample
 * @param PresentElems (O) Pointer to array of flags indicating presence for each element
 * @return Always 0
 */
int createPresentElems(int MaxZ, int NFoil, const foil *MatArray, int **PresentElems)
{  
  const COMPOUND *pcmp;
  int tempZ;
  int i,j;
  
  /*Initialize and fill the PresentElems Array (it is calloc'ed so by default contains 0's*/
  *PresentElems=(int*)calloc((MaxZ+1),sizeof(int));
  if (*PresentElems==NULL){fprintf(stderr,"\n Error allocating memory (PresentElems)\n");exit(1);}
  for(i=0;i<NFoil; i++) {
    pcmp=&MatArray[i].comp;
    for(j=0 ; j< pcmp->nelem ; j++)  {
      tempZ=pcmp->elem[j].Z;
      if(tempZ<=MaxZ) (*PresentElems)[tempZ]=1;
      else {fprintf(stderr,"\n ERROR: in createPresentElems(): %d > MaxZ (=%d)\n",tempZ,MaxZ);exit(1);}
    }
  }
  return(0);
}


/**Reads files containing the Calculation Flags (i.e. which lines should be calculated) or,
 *alternativelly reads areas files from a file in the DATTPIXE Areas format. The behaviour
 * depends on the DTAreasFlag (see below)
 * 
 * The Format of the file is as follows:
      -----------------
      (any header)
      DATA 
      (datablocks)
      -----------------
      
  if AreasFormat=1, The datablock for each element has the following structure:
    
    [ChemSymb] 
    [Kalpha1,2]  [Kbeta1]   [Kbeta2]  
    [Lalpha1,2]  [Lbeta2]   [Ll]  
    [Lbeta1]     [Lgamma1]  [Leta] 
    [Lbeta3]     [Lbeta4]   [Lgamma3]
    [Malpha1,2]  [Mbeta]    [Mgamma]
 
 Note: the flags are: 0 if not to be calculated, 1 (or >1) if should be calculated
 
 If AreasFormat=2 the datablock format is:
    [ChemSymb] 
    [Kalpha1,2]  [err]  [Kbeta1]  [err]  [Kbeta2]  [err]
    [Lalpha1,2]  [err]  [Lbeta2]  [err]  [Ll]      [err]  
    [Lbeta1]     [err]  [Lgamma1] [err]  [Leta]    [err] 
    [Lbeta3]     [err]  [Lbeta4]  [err]  [Lgamma3] [err]
    [Malpha1,2]  [err]  [Mbeta]   [err]  [Mgamma]  [err]
 *
 * Note: The information for elements not present in the sample is just ignored.
 *
 * @param CalcFlagsFileNm (I) Name of the input file
 * @param PresentElems (I) Array of flags indicating presence for each element. See createPresentElems()
 * @param MaxZinsample (I) Largest atomic number present in sample (dimension of PresentElems).
 * @param AreasFormat (I) Defines the format of the Areas File. 1 for raw areas format, 2 for DATTPIXE Areas Format
 * @param CalcFlags (O) Array (dimension is allocated to MaxZinsample+1) of structures containing the flags.
 * @return Always 0
 */
int readCalcFlags(const char *CalcFlagsFileNm, const int *PresentElems, int MaxZinsample, int AreasFormat,  CPIXERESULTS **CalcFlags){
  
  int i,j,tempZ,nread=0;
  char dummy[LONGSTRINGLENGTH];
  double trash;
  ChemSymb tempchem;
  CPIXERESULTS *pCFlags;
  FILE *CalcFlagsFile;
  /*Open File*/
  CalcFlagsFile = fopen(CalcFlagsFileNm, "r");
  if (CalcFlagsFile== NULL)
    {fprintf(stderr,"\nERROR: Could not open file '%s' for read.\n",CalcFlagsFileNm);exit(1);}
  
  /*Initialize CalcFlags */
  *CalcFlags=(CPIXERESULTS*)calloc(MaxZinsample+1,sizeof(CPIXERESULTS));
  if (*CalcFlags==NULL){fprintf(stderr,"\n Error allocating memory (CalcFlags)\n");exit(1);}
  
  /*Initialize the atnums of the present elems (to prevent crashes afterwards) */
  for(i=0;i<=MaxZinsample;i++)if(PresentElems[i])(*CalcFlags)[i].atnum=i;
    
  do {
    fscanf(CalcFlagsFile,"%255s",dummy);   
    if(feof(CalcFlagsFile)) {fprintf(stderr,"\nERROR: 'DATA' not found in '%s'\n",CalcFlagsFileNm);exit(1);}
  } while (strcasecmp(dummy, "DATA"));
 
  /*Read Data blocks*/
  for(;!feof(CalcFlagsFile);){
    fscanf(CalcFlagsFile,"%255s",dummy);
    
    if( (!strcmp(dummy,"END_DATA") && AreasFormat==1) || (AreasFormat==2 && !strcmp(dummy,"CALIB" )) ) {
      break;
    }
    /*
    else if(!strcmp(dummy,"INCLUDE")){
      fscanf(CalcFlagsFileNm,"%3s",tempchem);
    }
    */
    
    
    //If the reading continues, "dummy" contains either a chem. symbol or a Z value  
    if(AreasFormat==1){
      strncpy(tempchem,dummy,3);
      tempZ = symbol2Z(tempchem);
    }
    else if (AreasFormat==2){
      strncpy(tempchem,"-",3);
      tempZ = atoi(dummy);
    }
    else tempZ=0;
    
    if(tempZ<=MaxZinsample && PresentElems[tempZ]){
      nread++;

      pCFlags=&(*CalcFlags)[tempZ];
      pCFlags->atnum=tempZ;
      
      /*Read K-lines*/
      for(i=1;i<=3;i++) {
        fscanf(CalcFlagsFile, "%lg", &pCFlags->simareas.K_[i]);
        if(AreasFormat==2)fscanf(CalcFlagsFile, "%lg", &pCFlags->err.K_[i]);
      }
      /*Read L-lines*/
      for(i=1;i<=3;i++){
        for(j=1;j<=3;j++){
          fscanf(CalcFlagsFile, "%lg", &pCFlags->simareas.L_[i][j]);
          if(AreasFormat==2)fscanf(CalcFlagsFile, "%lg", &pCFlags->err.L_[i][j]);
        }
      }
      /*read M-lines*/
      for(i=1;i<=3;i++) {
        fscanf(CalcFlagsFile, "%lg", &pCFlags->simareas.M_[i]);
        if(AreasFormat==2)fscanf(CalcFlagsFile, "%lg", &pCFlags->err.M_[i]);
      }
      #if CPIXEVERBOSITY > 0  
      printf("\n CalcFlags[%2d] (%2s) read", pCFlags->atnum,tempchem);
      #endif
    }
    else {
      for(i=0;i<15;i++)fscanf(CalcFlagsFile, "%lg",&trash);//just discard the 3+9+3 following numbers
      if(AreasFormat==2)for(i=0;i<15;i++)fscanf(CalcFlagsFile, "%lg",&trash); //plus 15 more if err are present
      #if CPIXEVERBOSITY > 0  
      printf("\n Warning: Element '%s' ignored (not present in sample).",tempchem);
      #endif
    }
  }
  fclose(CalcFlagsFile);
  //Do some checks on the read data
  if(AreasFormat==2 && strcmp(dummy,"CALIB")) 
    {fprintf(stderr,"\nERROR: Bad DT format in '%s'\n",CalcFlagsFileNm);exit(1);}
  if(AreasFormat==1 && strcmp(dummy,"END_DATA"))
    {fprintf(stderr,"\nERROR: 'END_DATA' not found in '%s'\n",CalcFlagsFileNm);exit(1);}
  if(nread==0){fprintf(stderr,"\nERROR: '%s' contained no valid data\n",CalcFlagsFileNm);exit(1);}
  
  return(0);
}


/**Reads from an ascii file the calibration yields and also returns CalEner.
    @todo Modify to read only the present elems.
    
    Note: the index of XYldArray starts in 1.
    The File must have the following format:
      -----------------
      (any header)
      CALENER [calener]
      MAXZ [maxZ]
      
      DATA 
      (datablocks)
      -----------------
      
   The datablock for each element has the following structure:
    
    [AtNum]     [ChemSymb] 
    
    [Kalpha1,2]  [Kbeta1]   [Kbeta2]   <----Energies
    [Kalpha1,2]  [Kbeta1]   [Kbeta2]   <----Yields 
    
    [Lalpha1,2]  [Lbeta2]   [Ll]  <---------Energies      
    [Lalpha1,2]  [Lbeta2]   [Ll]  <---------Yields      
    [Lbeta1]     [Lgamma1]  [Leta] <--------Energies     
    [Lbeta1]     [Lgamma1]  [Leta] <--------Yields     
    [Lbeta3]     [Lbeta4]   [Lgamma3]  <----Energies
    [Lbeta3]     [Lbeta4]   [Lgamma3]  <----Yields
    
    [Malpha1,2]  [Mbeta]    [Mgamma]  <-----Energies  
    [Malpha1,2]  [Mbeta]    [Mgamma]  <-----Yields  

  
  -----------------

 * 
 * @param XYldFileNm (I) Name of the file to be read.
 * @param PresentElems (I) Array of flags indicating presence for each element. See createPresentElems()
 * @param CalcFlags (I) array of structures indicating which lines are to be calculated (simareas values>0) and which not (simareas values=-1). See readCalcFlags()
 * @param MaxZinsample (I) Largest atomic number present in sample (dimension of PresentElems).
 * @param CalEner (O) Energy at which the calibration was performed
 * @param XYldArray (O) Pointer to array of structs containing the efficiencies for each element
 * @return Always 0
 *
 * NOTE that the calib yield files from DT are in (cm2/ug)/uC  and should be passed to (cm2/1e15at)/uC
 */
int readXYld(const char *XYldFileNm, const int *PresentElems, const CPIXERESULTS *CalcFlags, int MaxZinsample, double *CalEner, XrayYield **XYldArray)
{
  FILE *XYldFile;
  char dummy[LONGSTRINGLENGTH];
  int maxZ;
  XrayYield *pXYld;
  int i,j,tempZ,result,iel;
  ChemSymb tempchem;
  double natmass,ug2at, trash;
  const CalibYld *pFlag;
  int atomicunits=0;
  
  XYldFile = fopen(XYldFileNm, "r");
  if (XYldFile == NULL)
    {fprintf(stderr,"\nERROR: Could not open file '%s' for read.\n",XYldFileNm);exit(1);}
    
  do {
    fscanf(XYldFile,"%255s",dummy);   
    if(feof(XYldFile)) {fprintf(stderr,"\nERROR: 'CALENER' not found in '%s'\n",XYldFileNm);exit(1);}
  } while (strcasecmp(dummy, "CALENER"));
  
  fscanf(XYldFile, "%lf", CalEner);
  
  do {
    fscanf(XYldFile,"%255s",dummy);   
    if(feof(XYldFile)) {fprintf(stderr,"\nERROR: 'MAXZ' not found in '%s'\n",XYldFileNm);exit(1);}
  } while (strcasecmp(dummy, "MAXZ"));
  
  fscanf(XYldFile, "%d", &maxZ);
  if(MaxZinsample>maxZ){fprintf(stderr,"\nERROR: Calib. file '%s' maxZ=%d<%d\n",XYldFileNm,maxZ,MaxZinsample);exit(1);}
  
  /*Initialize XYldArray */
  *XYldArray=(XrayYield*)calloc(MaxZinsample+1,sizeof(XrayYield));
  if (*XYldArray==NULL){fprintf(stderr,"\n Error allocating memory (XYldArray)\n");exit(1);}
  
  do {
    fscanf(XYldFile,"%255s",dummy);   
    if(feof(XYldFile)) {fprintf(stderr,"\nERROR: 'UNITS' not found in '%s'\n",XYldFileNm);exit(1);}
  } while (strcasecmp(dummy, "UNITS"));
  
  fscanf(XYldFile, "%255s",dummy);
  if(!strcasecmp(dummy, "atom"))atomicunits=1;
  else if(!strcasecmp(dummy, "mass"))atomicunits=0;
  else {fprintf(stderr,"\nERROR: unknown UNITS ('%s') in '%s'. Valid values are 'atom' or mass'\n",dummy,XYldFileNm);exit(1);}
  
  do {
    fscanf(XYldFile,"%255s",dummy);
    if(feof(XYldFile)) {fprintf(stderr,"\nERROR: 'DATA' not found in '%s'\n",XYldFileNm);exit(1);}
  } while (strcasecmp(dummy, "DATA"));
  
  /*Read Data blocks*/
  for(result=0,tempZ=0;!feof(XYldFile) && tempZ!=MaxZinsample;){
    fscanf(XYldFile, "%d", &tempZ);
    if(tempZ<1 || tempZ>MaxZinsample || (*XYldArray)[tempZ].atnum != 0){
      fprintf(stderr,"\nERROR: Bad format in Calib. file '%s' (Z=%d)\n",XYldFileNm,tempZ);
      exit(1);
    }
    fscanf(XYldFile,"%3s",tempchem);
    if(tempZ != symbol2Z(tempchem)){
      fprintf(stderr,"\nERROR: Bad format in Calib. file '%s' (Z=%d)\n",XYldFileNm,tempZ);
      exit(1);
    }
    if(PresentElems[tempZ]){
      result++;
      pXYld=&(*XYldArray)[tempZ];
      pXYld->atnum=tempZ;
      strcpy(pXYld->symb, tempchem);
      for(iel=0;CalcFlags[iel].atnum!=tempZ;iel++){} //findout the iel corresponding to tempZ
      
      pFlag=&CalcFlags[iel].simareas;
      
      if(!atomicunits){
        Z2mass(tempZ, &natmass,'n');
        ug2at=natmass/602.2141; // (Pm [g/mol] *1e6 [ug/g] / (Navo [at/mol])) * 1e15  
                                // ==> ug2at is in [ug/1e15at]
      }
      else ug2at=1; //If Units are already atomic, the conversion shouldn't be done
      
      /*Read K-lines*/
      for(i=1;i<=3;i++) fscanf(XYldFile, "%lg", &pXYld->ener.K_[i]);
      for(i=1;i<=3;i++) {
        fscanf(XYldFile, "%lg", &pXYld->XYld.K_[i]);
        pXYld->XYld.K_[i]*=ug2at;  // ug2at is in [ug/1e15at]
        //If the line is not to be used, put it at 0
        if(!(int)pFlag->K_[i]){ pXYld->XYld.K_[i]=0. ; pXYld->ener.K_[i]=0.;} 
      }
      /*Read L-lines*/
      for(i=1;i<=3;i++){
        for(j=1;j<=3;j++){
          fscanf(XYldFile, "%lg", &pXYld->ener.L_[i][j]);
        }
        for(j=1;j<=3;j++){
          fscanf(XYldFile, "%lg", &pXYld->XYld.L_[i][j]);
          pXYld->XYld.L_[i][j]*=ug2at;
          if(!(int)pFlag->L_[i][j]){ pXYld->XYld.L_[i][j]=0. ; pXYld->ener.L_[i][j]=0.;}
        }
      }
    
      /*read M-lines*/
      for(i=1;i<=3;i++) fscanf(XYldFile, "%lg", &pXYld->ener.M_[i]);
      for(i=1;i<=3;i++) {
        fscanf(XYldFile, "%lg", &pXYld->XYld.M_[i]);
        pXYld->XYld.M_[i]*=ug2at;
        if(!(int)pFlag->M_[i]){ pXYld->XYld.M_[i]=0. ; pXYld->ener.M_[i]=0.;}
      }
      #if CPIXEVERBOSITY > 0  
      printf("\n XyldArray[%2d] (%2s) read", pXYld->atnum,tempchem);
      #endif
    }
    else for(i=0;i<30;i++)fscanf(XYldFile, "%lg",&trash);//just discard the 3*2+9*2+3*2 following numbers
  }
  fclose(XYldFile);
  
  return(result);
}


/**Reads from an ascii file the Fluorescence and Coster-Kronig coefficients.
 Note: the index of FCKCoefArray starts in 1.
 @todo Modify to read only the present elems.
 
 The File must have the following format:
 
  -----------------
  (any header)
  
  MAXZ [maxZ]
  
  DATA 
  (datablocks)
  -----------------
  
  The datablock for each element has the following structure:
  
  
  [atnum]  [chemsymb]
  
  (Fluorescence Coeficients)
  [K] [LI] [LII] [LIII] [MI]... @todo Substitute the "..." by explicit labels
  
  [f12] [f13] [f23]
  
  (Line fractions for:) 
  [Kalpha1,2]  [Kbeta1]   [Kbeta2]   
  [Lalpha1,2]  [Lbeta2]   [Ll]        
  [Lbeta1]     [Lgamma1]  [Leta]      
  [Lbeta3]     [Lbeta4]   [Lgamma3]  
  [Malpha1,2]  [Mbeta]    [Mgamma]    


  -----------------
 
 * 
 * @param FCKCoefFileNm (I) Name of the file to be read.
 * @param MaxZinsample (I) Largest atomic number present in sample.
 * @param FCKCoefArray (O) Pointer to array of structs containing Fluorescence and C-K Coefs
 * @return Always 0
 */
int readFCK(char *FCKCoefFileNm, int MaxZinsample, FluorCKCoef **FCKCoefArray)
{
  FILE *FCKCoefFile;
  char dummy[LONGSTRINGLENGTH];
  int maxZ;
  FluorCKCoef *pFCKCoef;
  int i,tempZ;
  ChemSymb tempchem;
 
  maxZ=0;
  FCKCoefFile= fopen(FCKCoefFileNm, "r");
  if (FCKCoefFile == NULL)
    {fprintf(stderr,"\nERROR: Could not open file '%s' for read.\n",FCKCoefFileNm);exit(1);}
  
  do {
  fscanf(FCKCoefFile,"%255s",dummy);
  if(feof(FCKCoefFile)) {fprintf(stderr,"\nERROR: 'MAXZ' not found in '%s'\n",FCKCoefFileNm);exit(1);}
  } while (strcasecmp(dummy, "MAXZ"));
  
  fscanf(FCKCoefFile, "%d", &maxZ);
  if(MaxZinsample>maxZ){fprintf(stderr,"\nERROR: FCK file '%s' maxZ=%d<%d\n",FCKCoefFileNm,maxZ,MaxZinsample);exit(1);}
  
  /*Initialize FCKCoefArray*/
  *FCKCoefArray=(FluorCKCoef*)calloc(MaxZinsample+1,sizeof(FluorCKCoef));
  if (*FCKCoefArray==NULL){fprintf(stderr,"\n Error allocating memory (FCKCoefArray)\n");exit(1);}
  
  do {
  fscanf(FCKCoefFile,"%255s",dummy);   
  if(feof(FCKCoefFile)) {fprintf(stderr,"\nERROR: 'DATA' not found in '%s'\n",FCKCoefFileNm);exit(1);}
  } while (strcasecmp(dummy, "DATA"));
  
  /*Read Data blocks*/
  for(tempZ=0;!feof(FCKCoefFile) && tempZ!=MaxZinsample;){
    fscanf(FCKCoefFile, "%d", &tempZ);
    if(tempZ<1 || tempZ>MaxZinsample || (*FCKCoefArray)[tempZ].atnum != 0){
      fprintf(stderr,"\nERROR: Bad format in FCKCoef file '%s' (Z=%d)\n",FCKCoefFileNm,tempZ);
      exit(1);
    }
    fscanf(FCKCoefFile,"%3s",tempchem);
    if(tempZ != symbol2Z(tempchem)){
      fprintf(stderr,"\nERROR: Bad format in FCKCoef file '%s' (Z=%d)\n",FCKCoefFileNm,tempZ);
      exit(1);
    }
    pFCKCoef=&(*FCKCoefArray)[tempZ];
    pFCKCoef->atnum=tempZ;
    
    /*Read Fluorescence coeffs*/
    for (i = 1; i <= 9; i++) fscanf(FCKCoefFile, "%lg", &pFCKCoef->w[i]);
    
    /*Read f12,f13,f23 coeffs (Coster-Kronig)*/
    fscanf(FCKCoefFile, "%lg%lg%lg", &pFCKCoef->ck[0], &pFCKCoef->ck[1], &pFCKCoef->ck[2]);
    
    /*Read Line fractions*/
    fscanf(FCKCoefFile, "%lg%lg%lg", &pFCKCoef->k.K_[1], &pFCKCoef->k.K_[2], &pFCKCoef->k.K_[3]);
    for (i = 1; i <= 3; i++) fscanf(FCKCoefFile, "%lg%lg%lg", &pFCKCoef->k.L_[i][1], &pFCKCoef->k.L_[i][2], &pFCKCoef->k.L_[i][3]);
    fscanf(FCKCoefFile, "%lg%lg%lg", &pFCKCoef->k.M_[1], &pFCKCoef->k.M_[2], &pFCKCoef->k.M_[3]);
   
//     printf("\n DEBUG: FCKCoef[%2d] (%2s) read", pFCKCoef->atnum,tempchem);
  }
  fclose(FCKCoefFile);
  return(0);
}



/**Reads from an ascii file the X-ray Absorption coefficients.
 The File must have the following format:
 Note: the index of TotAbsCoefArray starts in 1.
 
  -----------------
  (any header)
  
  MAXZ [maxZ]
  
  DATA 
  (datablocks)
  -----------------
  
  The datablock for each element has the following structure:
  
  [atnum]  [chemsymb]
  
  (23 coefs each one for a given standard energy (see stdtable in TotAbsor() ) 
  [stdcoef[0]]....[stdcoef[22]]
  
  (9 triplets --1 for K, 3 for L and 5 for M-- describing absorp edges for K, L & M (1+3+5 triplets). First member of triplet is the transition ener, 2nd is abs coef for E below E_edge and 3rd is  abs coef for E above E_edge )
  [enr[1]] [coefenr[0][0]] [coefenr[0][1]]
  ...
  [enr[9]] [coefenr[8][0]] [coefenr[8][1]]
  
  -----------------
    
  @todo Note that enr indexing starts in 1 while the stdcoeffs & coefenr start in 0. A conversion from 1-->0 could be done with little work.
  
  IMPORTANT: The Absorption coeffs are tabulated using cm2/mg as unit but they need to be converted to 
             (cm2/1e15at) 
 
 * 
 * @param TotAbsCoefFileNm (I) Name of the file to be read. Coefs are expected in (cm2/mg)
 * @param MaxZinsample (I) Largest atomic number present in sample.
 * @param PresentElems (I) Array of flags indicating presence for each element. See createPresentElems()
 * @param Filter (I) Pointer to filter structure.
 * @param TotAbsCoefArray (O) Pointer to array of structs containing the Absorption Coefs (in cm2/1e15at ) (see readAbsCoef() ) 
 * @return  Always 0
 */
int readAbsCoef(char *TotAbsCoefFileNm, int MaxZinsample, const int *PresentElems, const FILTER *Filter, AbsCoef **TotAbsCoefArray)
{
  FILE *TotAbsCoefFile;
  char dummy[LONGSTRINGLENGTH];
  int maxZ;
  AbsCoef *pTotAbsCoef;
  int i,tempZ;
  ChemSymb tempchem;
  double trash, natmass;
  double mg2at; //conversion factor from cm2/mg to cm2/(1e15at) , which depends on the mass. 
   
  TotAbsCoefFile= fopen(TotAbsCoefFileNm, "r");
  if (TotAbsCoefFile == NULL)
    {fprintf(stderr,"\nERROR: Could not open file '%s' for read.\n",TotAbsCoefFileNm);exit(1);}
  
  do {
  fscanf(TotAbsCoefFile,"%255s",dummy);
  if(feof(TotAbsCoefFile)) {fprintf(stderr,"\nERROR: 'MAXZ' not found in '%s'\n",TotAbsCoefFileNm);exit(1);}
  } while (strcasecmp(dummy, "MAXZ"));
  
  fscanf(TotAbsCoefFile, "%d", &maxZ);
  if(MaxZinsample>maxZ){fprintf(stderr,"\nERROR: X-Absorption file '%s' maxZ=%d<%d\n",TotAbsCoefFileNm,maxZ,MaxZinsample);exit(1);}
  
  /*Initialize TotAbsCoefArray*/
  *TotAbsCoefArray=(AbsCoef*)calloc(MaxZinsample+1,sizeof(AbsCoef));
  if (*TotAbsCoefArray==NULL){fprintf(stderr,"\n Error allocating memory (TotAbsCoefArray)\n");exit(1);}
  
  do {
  fscanf(TotAbsCoefFile,"%255s",dummy);   
  if(feof(TotAbsCoefFile)) {fprintf(stderr,"\nERROR: 'DATA' not found in '%s'\n",TotAbsCoefFileNm);exit(1);}
  } while (strcasecmp(dummy, "DATA"));
  
  /*Read Data blocks*/
  for(tempZ=0;!feof(TotAbsCoefFile) && tempZ!=MaxZinsample;){
    fscanf(TotAbsCoefFile, "%d", &tempZ);
    if(tempZ<1 || tempZ>MaxZinsample || (*TotAbsCoefArray)[tempZ].atnum != 0){
      fprintf(stderr,"\nERROR: Bad format in X-Absorption file '%s' (Z=%d)\n",TotAbsCoefFileNm,tempZ);
      exit(1);
    }
    fscanf(TotAbsCoefFile,"%3s",tempchem);
    if(tempZ != symbol2Z(tempchem)){
      fprintf(stderr,"\nERROR: Bad format in X-Absorption file '%s' (Z=%d)\n",TotAbsCoefFileNm,tempZ);
      exit(1);
    }
    if(PresentElems[tempZ] || ( tempZ<=Filter->MaxZ && Filter->FilterElems[tempZ] ) ){
      pTotAbsCoef=&(*TotAbsCoefArray)[tempZ];
      pTotAbsCoef->atnum=tempZ;
      
      Z2mass(tempZ, &natmass,'n');
      mg2at=natmass/602214.1; // (Pm [g/mol] *1e3 [mg/g] / (Navo [at/mol])) * 1e15  
                              // ==> mg2at is in [mg/1e15at]
      
      for (i = 0; i < 23; i++){
        fscanf(TotAbsCoefFile, "%lg", &pTotAbsCoef->coefstd[i]);
        /*Transform: from cm2/mg to cm2/(1e15at)*/
        pTotAbsCoef->coefstd[i]*=mg2at;
      }
      for (i = 1; i <= 9; i++){
        fscanf(TotAbsCoefFile, "%lg%lg%lg", &pTotAbsCoef->enr[i], &pTotAbsCoef->coefenr[i-1][0],&pTotAbsCoef->coefenr[i-1][1]);
        /*Transform: from cm2/mg to cm2/(1e15at)*/
         pTotAbsCoef->coefenr[i-1][0]*=mg2at;
         pTotAbsCoef->coefenr[i-1][1]*=mg2at;
      }
      #if CPIXEVERBOSITY > 0
      printf("\n TotAbsCoef[%2d] (%2s) read", pTotAbsCoef->atnum,tempchem);
      #endif
    }
    else for (i = 0; i < 50; i++) fscanf(TotAbsCoefFile, "%lg", &trash); //just discard the 23+9*3 following numbers
  }
  fclose(TotAbsCoefFile);
  return(0);
}



/**Generates the stopping force tables 
 * Returns an array of pointers referencing to one table for each element present
 * in sample)
 * 
 * IMPORTANT: The tables are generated in the following units: keV/(1e15at/cm2) !!!
 *
 * @param path (I) String containing the path to look for the data base files.
 *              (e.g. for SCOEF.95A and SCOEF95B in the case of ZBL96 stoppings)
 *               Note that the path MUST contain the last "directory separator" 
 * @param pexp (I) Pointer to structure containing experimental parameters.
 * @param MaxZinsample (I) Largest atomic number present in sample.
 * @param PresentElems (I) Array of flags indicating presence for each element. See createPresentElems()
 * @param step (I) Energy step to be used (in keV) in the case of linear energy scale
 *             or the multiplicative factor in the case logarithmic energy scale.
 * @param SPTArray (O) Array of Stopping power table structs for each present element. 
 *                 IMPORTANT: The tables are generated in the following units: keV/(1e15at/cm2) !!!
 * @return Always 0
 */
int createSPTs(const char *path, const EXP_PARAM *pexp, int MaxZinsample, const int *PresentElems, double step, SPT **SPTArray)
{
  int Z,j;
  SPT *pSPT;
  double E;
  double *nucSP;
  char afilename[FILENMLENGTH], bfilename[FILENMLENGTH];  
  /**@todo readstopcoef assumes that the files SCOEF.95A and SCOEF.95B are located where the program is called. This must be generalized (In fact, the coefficients could be hardcoded using an include file)*/
  snprintf(afilename,FILENMLENGTH,"%sSCOEF.95A",path);
  snprintf(bfilename,FILENMLENGTH,"%sSCOEF.95B",path);
  readstopcoef(afilename, bfilename);
  
  /*Initialize SPTables*/
  /*Calloc is used, so all the non-explicitely initialized, remain 0/NULL */
  *SPTArray=(SPT*)calloc(MaxZinsample+1,sizeof(SPT));
  if (*SPTArray==NULL){fprintf(stderr,"\n Error allocating memory (SPTArray)\n");exit(1);}
  
  for(Z=1; Z<=MaxZinsample ;Z++){
    if(PresentElems[Z]){
      pSPT=&(*SPTArray)[Z];
      
      pSPT->Emax=1.1*pexp->BeamEner;
      pSPT->Estep=step;  
      pSPT->logmode=0;  ///<For the moment, only linear Energy scale is implemented
      pSPT->nrows=1+(int)ceil(pSPT->Emax/pSPT->Estep);
      
      /*Initialize the arrays*/
      pSPT->E=(double*)malloc(pSPT->nrows*sizeof(double));
      if (pSPT->E==NULL){fprintf(stderr,"\n Error allocating memory (SPTArray[%d].E)\n",Z);exit(1);}
      pSPT->S=(double*)malloc(pSPT->nrows*sizeof(double));
      if (pSPT->S==NULL){fprintf(stderr,"\n Error allocating memory (SPTArray[%d].S)\n",Z);exit(1);}
      pSPT->dSdE=(double*)malloc(pSPT->nrows*sizeof(double));
      if (pSPT->dSdE==NULL){fprintf(stderr,"\n Error allocating memory (SPTArray[%d].dSdE)\n",Z);exit(1);}
      
      nucSP=(double*)malloc(pSPT->nrows*sizeof(double));
      if (nucSP==NULL){fprintf(stderr,"\n Error allocating memory (nucSP)\n");exit(1);}
      
      for(j=0 ,E=0; j< pSPT->nrows ; j++, E+=step)pSPT->E[j]=E;
      
      stop96d(pexp->ion.Z, Z, pexp->ion.IM,  pSPT->E, pSPT->nrows, pSPT->S);
      nuclearstopping_ZBL(pexp->ion.Z, Z, pexp->ion.IM, (double)(-1.), pSPT->E, pSPT->nrows, nucSP);
      
      for(pSPT->Smax=0., j=0 ; j< pSPT->nrows ; j++){
        pSPT->S[j] += nucSP[j];
        pSPT->S[j] *=0.001; //Return stopping in keV/(1e15at/cm2)
        if(pSPT->Smax < pSPT->S[j]) pSPT->Smax = pSPT->S[j];  //obtain the maximum of S
      }
      for(j=0 ; j< pSPT->nrows-1 ; j++){
      pSPT->dSdE[j]=(pSPT->S[j+1] - pSPT->S[j])/(pSPT->E[j+1] - pSPT->E[j]);
      }
      pSPT->dSdE[pSPT->nrows-1]=pSPT->dSdE[pSPT->nrows -2];
      
      free(nucSP);
      #if CPIXEVERBOSITY > 0
      printf("\n Created SP table for Z=%2d (%d rows)",Z,pSPT->nrows);
      #endif
    }
  }
  
  return(0);
}


/**Returns stopping Using previously calculated tables.
 * Performs a linear interpolation and a Bragg adition.
 * 
 * @param pcmp (I) Pointer to target definition
 * @param SPTArray (I) Pointer to precalculated tables of stoppings (see createSPT() )
 * @param E (I) Energy of the ion, in keV
 * @return Stopping Power, in keV/(1e15at/cm2)  IMPORTANT!!! note that it's keV and ot eV
 */
double getSP(const COMPOUND *pcmp, const SPT *SPTArray, double E)
{
  int i,j;
  double SP;
  const SPT *pspt;
  
  /*For each element in the compound, obtain the stopping power from its table and apply bragg rule*/
  for(SP=0, i=0 ; i< pcmp->nelem ; i++){
    pspt= &SPTArray[pcmp->elem[i].Z];
    j=(int)floor(E / pspt->Estep);  // j is the index in the table (index of nearest lower tabulated energy)
    SP+= (pspt->S[j] + pspt->dSdE[j]*(E - pspt->E[j]) )* pcmp->xn[i];    //linear interpolation + Bragg rule
  }
  return (SP);
}



/**Initializes the lyr array 
 * 
 * @param pexp (I) Pointer to structure containing experimental parameters.
 * @param MatArray (I) Array of foils defining the current sample
 * @param XYldArray (I) Array of Xray yield coefs. See readXYld()
 * @param TotAbsCoefArray (I) Array containing a tables of Absorption Coefs (see readAbsCoef() )
 * @param SPTArray (I) Pointer to precalculated tables of stoppings (see createSPT() )
 * @param NFoil (I) Number of target foils
 * @param PresentElems (O) Pointer to array of flags indicating presence for each element
 * @param NFoilUsed (O) Number of foils finally used in the calcs 
 * @param plyrarray (O) Pointer to the layer array.
 * @return Always 0
 */
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 i;
  int Eovermin;
  double MajAbsCoef;
  LYR *plyr, *plyrold;
    
  *plyrarray=(LYR*)calloc(NFoil+1,sizeof(LYR));
  if (*plyrarray==NULL){fprintf(stderr,"\n Error allocating memory (lyrarray)\n");exit(1);}
  
  /*The following definitions for layer 0 are done for convenience even when first "real" layer has label 1.*/
  plyr=&(*plyrarray)[0];
  
//   (*plyrarray)[0].FoilOutEner=pexp->BeamEner; 
//   (*plyrarray)[0].absolutepos=0.;
//   (*plyrarray)[0].ThickIn=0.;
  plyr->FoilOutEner=pexp->BeamEner;
  plyr->absolutepos=0.;
  plyr->ThickIn=0.;

  for(Eovermin=1, i=1; i<=NFoil && Eovermin; i++){
//     printf("\n\nDEBUG:  Layer %d: \n",i);   
//     (*plyrarray)[i].absolutepos=(*plyrarray)[i-1].absolutepos + (*plyrarray)[i-1].ThickIn;
//     (*plyrarray)[i].FoilInEner=(*plyrarray)[i-1].FoilOutEner;
//     initlyr(pexp, &MatArray[i-1], XYldArray, TotAbsCoefArray, &(*plyrarray)[i], &MajAbsCoef);
//     createsublyrs(pexp, &MatArray[i-1], SPTArray, MajAbsCoef, &(*plyrarray)[i] );
//     Eovermin=((*plyrarray)[i].FoilOutEner > pexp->FinalEner);
    plyrold=plyr;
    plyr=&(*plyrarray)[i];
    plyr->absolutepos=plyrold->absolutepos + plyrold->ThickIn;
    plyr->FoilInEner=plyrold->FoilOutEner;
    plyr->pFoil=&MatArray[i-1];
    initlyr(pexp, XYldArray, TotAbsCoefArray, PresentElems, plyr, &MajAbsCoef);
    createsublyrs(pexp, &MatArray[i-1], SPTArray, MajAbsCoef, plyr);
    Eovermin=(plyr->FoilOutEner > pexp->FinalEner);
  }
  
  /*if the lyr initialization loop was exited without initializing some layers, 
    reallocate the lyr structure and change the (used) number of foils.
    */
  *NFoilUsed=i-1;
  if(*NFoilUsed<NFoil){
    #if CPIXEVERBOSITY >0
    printf("\n\nWarning: Layer %d (and deeper) Have been ignored during calculation\n",*NFoilUsed+1);
    #endif
  }
  
  return(0);
}


/** Initializes ONE given layer.
 *  It allocates space for various arrays and sets initial values for  TrcUse, ResYldArray,...
 *  It also calculates the absorption and returns the Maximum absorption coefficient (MACoef).
 * 
 * Note: it uses the global vars minK, maxK, minM, minL,...
 * 
 * @param pexp (I) Pointer to structure containing experimental parameters.
 * @param XYldArray (I) Array of Xray yield coefs. See readXYld()
 * @param TotAbsCoefArray (I) Array containing a tables of Absorption Coefs (see readAbsCoef() )
 * @param plyr (I+O) Pointer to current layer struct.
 * @param PresentElems (O) Pointer to array of flags indicating presence for each element
 * @param MACoef (I+O) Maximum absorption coefficient to be used in calcs (presently only cares about K&Lcoefs) 
 * @return Always 0
 */
int initlyr(const EXP_PARAM *pexp, const XrayYield *XYldArray, const AbsCoef *TotAbsCoefArray, const int *PresentElems, LYR *plyr, double *MACoef)
{
  int i,j,jj,ll,Z;
  /*Auxiliary vars (or pointers):*/
  CalibYld *AbsFac;  
  XrayYield *ResYld;
  double geomcorr;
  
  
  plyr->NumOfTrc=plyr->pFoil->nfoilelm;
  
//   plyr->ResArray=(XrayYield*)calloc(plyr->NumOfTrc,sizeof(XrayYield));
//   if (plyr->ResArray==NULL){fprintf(stderr,"\n Error allocating memory (plyr->ResArray)\n");exit(1);}
  
  plyr->ResYldArray=(XrayYield*)calloc(plyr->NumOfTrc,sizeof(XrayYield));
  if (plyr->ResYldArray==NULL){fprintf(stderr,"\n Error allocating memory (plyr->ResYldArray)\n");exit(1);}
  
  plyr->SecResArray=(SecResType*)calloc(plyr->NumOfTrc,sizeof(SecResType));
  if (plyr->SecResArray==NULL){fprintf(stderr,"\n Error allocating memory (plyr->SecResArray)\n");exit(1);}
  
  /**CAUTION: Change with respect to CDT: TrcUse now starts with index 0*/
  plyr->TrcUse=(int*)calloc(plyr->NumOfTrc,sizeof(int));
  if (plyr->TrcUse==NULL){fprintf(stderr,"\n Error allocating memory (plyr->TrcUse)\n");exit(1);}
  
  /**CAUTION: Change with respect to CDT: NeedSFC now starts with index 0*/
   plyr->NeedSFC=(int*)calloc(plyr->NumOfTrc,sizeof(int));
   if (plyr->NeedSFC==NULL){fprintf(stderr,"\n Error allocating memory (plyr->NeedSFC)\n");exit(1);}
  ///Note: @todo plyr->NeedSFC is just a convenience pointer to the Sample-general NeedSFC. This should be changed  and not used, but for the moment it works.
//   plyr->NeedSFC=NULL;   

   /**CAUTION: Change with respect to CDT: TrcAtnum now starts with index 0*/
  plyr->TrcAtnum=(int*)calloc(plyr->NumOfTrc,sizeof(int));
  if (plyr->TrcAtnum==NULL){fprintf(stderr,"\n Error allocating memory (plyr->TrcAtnum)\n");exit(1);}
    
  
  plyr->dimTrans=pexp->simpar.MaxZinsample + 1;  //The dimenssion of the Transmission vector must be enough to hold the element with Z=MaxZinsample
  plyr->Trans=(CalibYld*)calloc(plyr->dimTrans,sizeof(CalibYld));
  if (plyr->Trans==NULL){fprintf(stderr,"\n Error allocating memory (plyr->Trans)\n");exit(1);}
  
  /*Calculation of Transmission coefficient in this layer for each element present in the sample*/
  geomcorr= 1. / pexp->cosDet; //the thickness will be corrected considering the detector direction
  for(Z=1;Z<plyr->dimTrans;Z++){
    if(PresentElems[Z]){
      Transmission(XYldArray, TotAbsCoefArray, Z, plyr->pFoil, geomcorr, (CalibYld*)NULL, &plyr->Trans[Z]);
    }
  } 
  
  
  /*Initializations for each element in the layer*/
  for (*MACoef=0, i = 0; i < plyr->NumOfTrc; i++){
    /*Select which elements should be simulated*/
    plyr->TrcUse[i]= (plyr->pFoil->comp.elem[i].Z > minK);
    if (plyr->TrcUse[i]) {
      /**Note: At this point in the CDT code, CalYld is defined, but it makes no sense since 
      it is just a copy of XYldArray[].XYld (which is already there and never changing...
      So whenever CalYld[i] is needed, use XYldArray[i].XYld instead*/
      /*Initialization for XYldArray*/
      plyr->ResYldArray[i]=XYldArray[plyr->pFoil->comp.elem[i].Z];
      /*Some (stupid?) initialization for SecResArray*/
      plyr->SecResArray[i].Pk.atnum=plyr->pFoil->comp.elem[i].Z;
      for(j=0;j<3;j++){plyr->SecResArray[i].Pk.area[j][0]=1; plyr->SecResArray[i].Pk.area[j][1]=1;}
      
      /*Calculate absorption coefficients*/
      ResYld=&(plyr->ResYldArray[i]);
      plyr->TrcAtnum[i] = ResYld->atnum;   //<Fills the TrcAtnum Array
      AbsFac= &(plyr->SecResArray[i].SR.AbsFact);
      
      if (ResYld->ener.K_[1] > 0 && ResYld->atnum < maxK) {
        for (jj = 1; jj <= 3; jj++) {
          AbsFac->K_[jj] = TotAbsor(TotAbsCoefArray, &plyr->pFoil->comp, ResYld->ener.K_[jj]);
          if (*MACoef < AbsFac->K_[jj]) *MACoef = AbsFac->K_[jj];
        }
      }
      if (ResYld->ener.M_[1] > 0 && ResYld->atnum > minM) {
        for (jj = 1; jj <= 3; jj++) {
          AbsFac->M_[jj] = TotAbsor(TotAbsCoefArray, &plyr->pFoil->comp, ResYld->ener.M_[jj]);
          ///@todo TODO:Check why was the following line commented out (in the Pascal Code):
          /* if *MACoef<AbsC[ii].M[jj] then *MACoef:=AbsC[ii].M[jj] ;  */
        }
      }
      if (ResYld->ener.L_[1][1] > 0 && ResYld->atnum > minL) {
        for (jj = 1; jj <= 3; jj++) {
          for (ll = 1; ll <= 3; ll++) {
            AbsFac->L_[jj][ll] = TotAbsor(TotAbsCoefArray, &plyr->pFoil->comp, ResYld->ener.L_[jj][ll]);
            if (*MACoef < AbsFac->L_[jj][ll])  *MACoef = AbsFac->L_[jj][ll];
          }
        }
      }
    }
  }
  return(0);
}


int readFilter(const char *FilterDefFileNm, FILTER *Filter)
{
  
  FILE *FilterDefFile;
  char dummy[LONGSTRINGLENGTH];
  int i,tempZ,j;
  foil *pfoil;
  double mg2at,natmass;
  
  
  Filter->geomcorr=1.;  ///<@todo For the moment, the geometric correction is just initialized to 1
  
  FilterDefFile = fopen(FilterDefFileNm, "r");
  if (FilterDefFile == NULL)
    {fprintf(stderr,"\nERROR: Could not open file '%s' for read.\n",FilterDefFileNm);exit(1);}
      
  /*Skip until "NUMBER OF FOILS"*/
  do {
    fscanf(FilterDefFile,"%255s",dummy);   
  } while (strcasecmp(dummy, "NUMBER_OF_FOILS"));  
  
  fscanf(FilterDefFile, "%d", &Filter->nlyr);
  
  /*Initialize  Sample*/
  //if(*Sample!=NULL)free(*Sample);
  Filter->foil=(foil*)malloc((Filter->nlyr)*sizeof(foil));
  if (Filter->foil==NULL){fprintf(stderr,"\n Error allocating memory (Filter->foil)\n");exit(1);}
  
  for (Filter->MaxZ=0, i = 0; i < Filter->nlyr ; i++) {
    pfoil=&(Filter->foil[i]);
    do {
      fscanf(FilterDefFile,"%255s",dummy);   
      } while (strcasecmp(dummy, "FOIL"));  
    
    fscanf(FilterDefFile, "%lg %10s", &pfoil->thick, dummy);
    if(strcasecmp(dummy,"cm2") && strcasecmp(dummy,"mg")){
      fprintf(stderr,"\n Error: Bad syntax in filter definition file: '%s' is not a known thickness unit\n",dummy);exit(1);
    }
    fscanf(FilterDefFile, "%d",  &pfoil->nfoilelm);
    tempZ=readCOMPOUND(FilterDefFile, pfoil->nfoilelm, &pfoil->comp);
    if(Filter->MaxZ < tempZ) Filter->MaxZ = tempZ;
    
    if(!strcasecmp(dummy,"mg")){
      /*Convert the value in pfoil->thick from mg/cm2 to 1e15at/cm2 */
      for(natmass=0, j=0 ; j<pfoil->nfoilelm ; j++) natmass += pfoil->comp.xn[j] * pfoil->comp.elem[j].M;
      mg2at=natmass/602214.1; // (Pm [g/mol] *1e3 [mg/g] / (Navo [at/mol])) * 1e15  
                              // ==> mg2at is in [mg/1e15at]
      pfoil->thick /= mg2at;  //Before this, thick was in mg/cm2, hence divide by mg2at to get 1e15at/cm2
    }
    
  }
  fclose(FilterDefFile);
  
  createPresentElems(Filter->MaxZ, Filter->nlyr, Filter->foil, &Filter->FilterElems);
  
  return(0);
}

/**@todo Document this function*/
int FilterTrans(const EXP_PARAM *pexp, const XrayYield *XYldArray, const AbsCoef *TotAbsCoefArray, const int *PresentElems, FILTER *Filter)
{
  int Z,i;
  CalibYld auxTrans, *pTrans, *pauxTrans;
  CalibYld unitTrans={ { 0.0, 1.0, 1.0, 1.0 },
                       { { 0.0, 0.0, 0.0, 0.0 },
                         { 0.0, 1.0, 1.0, 1.0 },
                         { 0.0, 1.0, 1.0, 1.0 },
                         { 0.0, 1.0, 1.0, 1.0 }},
                       { 0.0, 1.0, 1.0, 1.0 }};
  
  Filter->dimTrans=pexp->simpar.MaxZinsample + 1;  //The dimenssion of the Transmission vector must be enough to hold the element with Z=MaxZinsample
  Filter->Trans=(CalibYld*)calloc(Filter->dimTrans,sizeof(CalibYld));
  if (Filter->Trans==NULL){fprintf(stderr,"\n Error allocating memory (Filter->Trans)\n");exit(1);}
   
  for(Z=1;Z<Filter->dimTrans;Z++){
    for(Z=1;Z<Filter->dimTrans;Z++)
    if(PresentElems[Z]){
      pauxTrans=NULL; //Initially, we have nothing to accumulate 
      pTrans=&Filter->Trans[Z];
      *pTrans=unitTrans;
   
      for (i=0 ; i < Filter->nlyr ; i++){
        Transmission(XYldArray, TotAbsCoefArray, Z, &Filter->foil[i], Filter->geomcorr, pauxTrans, pTrans);
        auxTrans= *pTrans;  //backup the last pTrans values (not the pointer but the values!)
        pauxTrans=&auxTrans; //pauxTrans is set to the last foil transmission for accumulating
      }
    }
  }
  
  return(0);
}


/** Calculates the Transmission coefficients for a given foil (whith effective thickness calculated
 * using a geometric correction).
 * 
 * @param XYldArray (I) Array of Xray yield coefs. See readXYld()
 * @param TotAbsCoefArray (I) Array containing a tables of Absorption Coefs (see readAbsCoef() )
 * @param Z (I) Atomic number of the element for whose lines the transmission is calculated
 * @param pFoil (I) pointer to the current foil definition 
 * @param geomcorr (I) geometric correction to the foil thickness (allows tiltings)
 * @param pTransOld (I) If this pointer is not null, pTrans will be multiplied by it line-by-line. 
 * @param pTrans (O) Structure where the transmission coefs for each line are returned.
 * @return Always 0
 */
int Transmission(const XrayYield *XYldArray, const AbsCoef *TotAbsCoefArray, int Z, const foil *pFoil, double geomcorr, const CalibYld *pTransOld, CalibYld *pTrans)
{
  int jj,ll;
  double thickout;
  const CalibYld *pXYldener;  
  
  thickout=pFoil->thick * geomcorr;  //Effective thickness 
  pXYldener=&XYldArray[Z].ener;  //energies of lines for element with at.number Z
  
  /*Transmission coefs for K lines*/
  if (pXYldener->K_[1] > 0 && Z < maxK) {
    for (jj = 1; jj <= 3; jj++) {
      pTrans->K_[jj] = exp(-thickout * TotAbsor(TotAbsCoefArray, &pFoil->comp, pXYldener->K_[jj]) );
      if(pTransOld!=NULL) 
        pTrans->K_[jj] *= pTransOld->K_[jj];
    }
  }
  /*Transmission coefs for M lines*/
  if (pXYldener->M_[1] > 0 && Z > minM) {
    for (jj = 1; jj <= 3; jj++) {
      pTrans->M_[jj] = exp(-thickout * TotAbsor(TotAbsCoefArray, &pFoil->comp, pXYldener->M_[jj]) );
      if(pTransOld!=NULL) pTrans->M_[jj] *= pTransOld->M_[jj];
    }
  }
  /*Transmission coefs for L lines*/
  if (pXYldener->L_[1][1] > 0 && Z > minL) {
    for (jj = 1; jj <= 3; jj++) {
      for (ll = 1; ll <= 3; ll++) {
        pTrans->L_[jj][ll] = exp(-thickout * TotAbsor(TotAbsCoefArray, &pFoil->comp, pXYldener->L_[jj][ll]) );
        if(pTransOld!=NULL) pTrans->L_[jj][ll] *= pTransOld->L_[jj][ll];
      }
    }
  }
  return(0);
}

/**Calculates the total absoption (it returns it as a double) for a given X-ray line in a given compound
 * 
 * @param TotAbsCoefArray (I) Array containing a tables of Absorption Coefs (see readAbsCoef() )
 * @param cmp (I) Definition of the absorber material.
 * @param Xray (I) Energy of the X-ray, in keV  
 * @return Total absorption coefficient
 */
double TotAbsor(const AbsCoef *TotAbsCoefArray, const COMPOUND *cmp, double Xray)
{
  int i,iener;
  double absaux;
  int dimstdener;
  const AbsCoef *Absco;
  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};
  
  dimstdener=23-1;                    
  /**<@todo This function returns 0 if the Xray energy is not in the range covered by stdener.
    Maybe an error (and exit) should be triggered instead. Moreover: Is this
    absolutely general? I think this limits the generality of the program since
    elements below Na cannot be treated (and this is not desirable from a
    generalistic point of view). Therefore, these limits should be coded as
    variables just as minK, maxK,...
    UPDATE 1: The limits can be expanded, but the X-absorption data table should
    be altered accordingly (and dimstdener too).
    UPDATE 2: Better to show a warning than to exit in case that the limits are
    violated (TODO). 
    */
  if (Xray < stdener[0] || Xray >= stdener[dimstdener])  return (0.);
  
  for(iener=1;Xray > stdener[iener];iener++);
  
  
  for (absaux=0, i = 0; i < cmp->nelem; i++) {
    Absco=&(TotAbsCoefArray[cmp->elem[i].Z]); 
    absaux += cmp->xn[i] * TotAbsor_elemental(iener, Absco, cmp->elem[i].Z, Xray);
  }
  
  return(absaux);
}

/**Calculates the total absoption coefficient (it returns it as a double) for a given X-ray line in a given element.
 *
 * @param iener (I) Index of the absorption line to be considered
 * @param Absco (I) Pointer to structure of Absorption Coefs for the absorber element
 * @param Z (I) Atomic number of the absorber element
 * @param Xray (I) Energy of the X-ray, in keV  
 * @return Total absorption coefficient
 */
double TotAbsor_elemental(int iener, const AbsCoef *Absco, int Z, double Xray)
{
  int j,limj;
  double absaux = 0.0;
  double enermaj, enermin, a, b, majcoef, mincoef;
  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};
                        
  /**@todo This function returns 0 if the Xray energy is not in the 1-30 range.
    Maybe an error (and exit) should be triggered instead. Moreover: Is this
    absolutely general? I think this limits the generality of the program since
    elements below Na cannot be treated (and this is not desirable from a
    generalistic point of view). Therefore, these limits should be coded as
    variables just as minK, maxK,...
    UPDATE 1: The limits can be expanded, but the X-absorption data table should
    be altered accordingly.
    UPDATE 2: Better to show a warning than to exit in case that the limits are
    violated (TODO). 
    */

  enermaj = stdener[iener];
  enermin = stdener[iener-1];
  majcoef = Absco->coefstd[iener];
  mincoef = Absco->coefstd[iener-1];
  
  if(Z<11) limj=0;     //No absorp edges in the 1-30keV range for elements of Z<11
  else{
  if(Z<28) limj=1;     //Only K absorp edges in the 1-30keV range for elements of Z<28
    else{
      if(Z<52) limj=4; //Only K & L absorp edges in the 1-30keV range for elements of Z<52
      else limj=9;                  // K, L & M  absorp edges in the 1-30keV range for elements of Z>52
    }  
  }
  /*Check wether the energy is affected by an absorption edge*/
  for (j = 1; j <= limj; j++) {
    if (Absco->enr[j] > 0 && Absco->enr[j] > enermin && Absco->enr[j] < Xray) {
      enermin = Absco->enr[j];
      mincoef = Absco->coefenr[j-1][1];
    }
    if (Absco->enr[j] > 0 && Absco->enr[j] < enermaj && Absco->enr[j] > Xray) {
      enermaj = Absco->enr[j];
      majcoef = Absco->coefenr[j-1][0];
    }
  }
  /*Logarithmic average*/
  a = log(majcoef / mincoef) / log(enermaj / enermin);
  b = log(majcoef) - a * log(enermaj);
  absaux = exp(a * log(Xray) + b);

  return(absaux);
}



/**Manages the creation of internal sublayers in each layer
 * 
 * @param pexp (I) Pointer to structure containing experimental parameters.
 * @param pMat (I) Pointer to current foil structure.
 * @param SPTArray (I) Pointer to precalculated tables of stoppings (see createSPT() )
 * @param MajAbsCoef (I) Maximum absorption coef for this layer.
 * @param plyr (I+O) Pointer to current layer struct.
 * @return Always 0
 */
int createsublyrs(const EXP_PARAM *pexp, const foil *pMat, const SPT *SPTArray, double MajAbsCoef, LYR *plyr)
{
  ESxType ESx1,ESx2,ESx3;
  int i,Fpos;
  double Range,range1,range2, Smax;
  
  //Before starting, find out the maximum stopping force for this layer (Smax)
  for(Smax=0., i=0 ; i < pMat->nfoilelm ; i++) Smax += SPTArray[pMat->comp.elem[i].Z].Smax * pMat->comp.xn[i];
   
 
  plyr->ThickIn=pMat->thick/pexp->cosInc;

  Fpos=1;
  plyr->ESxArray=(ESxType*)malloc((Fpos)*sizeof(ESxType));
  if (plyr->ESxArray==NULL){fprintf(stderr,"\n Error allocating memory (plyr->ESxArray)\n");exit(1);}
  
  ESx1.ep = plyr->FoilInEner;     //note that FoilInEner has been already initialized in initlyrarray()
  ESx1.stpp = getSP(&pMat->comp, SPTArray, ESx1.ep);
  ESx1.x = 0.0;
  plyr->ESxArray[Fpos-1]=ESx1;
  
  ///@todo Discuss this initialization with M.R.
  ESx2.ep = plyr->FoilInEner * 0.90;  
  ESx2.stpp = getSP(&pMat->comp, SPTArray, ESx2.ep);
  ESx2.x = 0.0; 
  
  ESx3.ep=pexp->FinalEner;
  ESx3.stpp = getSP(&pMat->comp, SPTArray, ESx3.ep);
  ESx3.x = 0.0;

//   ESx2.ep = plyr->FoilInEner - plyr->ThickIn * Smax; 
//   if(ESx2.ep < pexp->FinalEner)ESx2.ep= pexp->FinalEner;
//   ESx2.stpp = getSP(&pMat->comp, SPTArray, ESx2.ep);
//   ESx2.x = plyr->ThickIn; 
//
//   ESx3.ep=pexp->FinalEner;
//   ESx3.stpp = getSP(&pMat->comp, SPTArray, ESx3.ep);
//   ESx3.x = (plyr->FoilInEner - pexp->FinalEner)/Smax;

  
  SSThick(pexp, pMat, SPTArray, MajAbsCoef, plyr->ThickIn, ESx1, &ESx2, &Fpos, &range1, &plyr->ESxArray);

  ESx2=  plyr->ESxArray[Fpos - 1];

  if (ESx2.x == plyr->ThickIn) {
    plyr->FoilOutEner = ESx2.ep;
    Range = ESx2.x;
    plyr->FESxlen= Fpos;
  } 
  else {
    SSThick(pexp, pMat, SPTArray, MajAbsCoef, plyr->ThickIn, ESx2, &ESx3, &Fpos, &range2, &plyr->ESxArray);
    ESx3=plyr->ESxArray[Fpos-1];
    
    plyr->FoilOutEner = ESx3.ep;
    Range = range1 + range2;
    if (Range < plyr->ThickIn) plyr->ThickIn = Range;
    plyr->FESxlen = Fpos;
  }
  
  return(0);
}


/**Recursive function which divides a sublayer in two if either energy loss or absorption are too large
 * 
 * @param pexp (I) Pointer to structure containing experimental parameters.
 * @param pMat (I) Pointer to current foil structure.
 * @param SPTArray (I) Pointer to precalculated tables of stoppings (see createSPT() )
 * @param MajAbsCoef (I) Maximum absorption coef for this layer.
 * @param LayerThickness (I) Effective Thickness (in-going path) for the layer containing this sublayer.
 * @param ESxin (I) Struct of energy, stopping and depth (ESx) for first half of this sublayer
 * @param ESxfin (I+O) Pointer to struct of energy, stopping and depth (ESx) for 2nd half of this sublayer 
 * @param pFpos (I+O) Pointer to index of sublayer position into the layer.
 * @param thick (O) Pointer to thickness of current sublayer.
 * @param ESxA (O) Pointer to Array of sublayers in the current layer.
 * @return always 0.
 */
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)
{
  double tck1, tck2, auxdecis, auxstpchange,auxaux;
  ESxType ESxm1, ESxaux;

  auxaux=exp(-MajAbsCoef * ESxin.x * pexp->cosFac);
  if (auxaux > 1e-7){
    /*calculate the relative transmission 
      Note: auxdecis is the relative contribution to the transmission due to the present layer
            compared to the transmission from the present layer to the surface.
            The cosFac transforms distances from in-going to out-going paths.
    */ 
    auxdecis = exp(MajAbsCoef * (ESxin.x - ESxfin->x) * pexp->cosFac) /auxaux;
         
  }
  else auxdecis = 1.0;
    
  /*Calculate the relative change in stopping forces (with an ad-hoc weight based on the relative energy)
  */ 
  auxstpchange= (ESxin.ep / pexp->BeamEner)* fabs(ESxin.stpp - ESxfin->stpp) / ((ESxin.stpp + ESxfin->stpp)/2.);

  /* Exit condition (no more recursive calls)*/
  /**@todo: check if the thressholds of these conditions are reasonable or should be more conservative.
   */
   if ( auxstpchange < 0.005  && auxdecis > 0.997 ){    
//    if ( auxstpchange < 0.05   ){    //DEBUG!!!!!!!
    ESxaux.ep = (ESxin.ep + ESxfin->ep) / 2.; //mean energy
    ESxaux.stpp = getSP(&pMat->comp, SPTArray, ESxaux.ep);//stopping at mean energy 
    *thick = (ESxin.ep - ESxfin->ep) / ((ESxin.stpp + 4. *ESxaux.stpp  + ESxfin->stpp) / 6.); //
    ESxfin->x = ESxin.x + *thick;
    
    if (ESxfin->x > LayerThickness) {        
      *thick = LayerThickness - ESxin.x;
      ESxfin->x = LayerThickness;
      ESxfin->ep = ESxin.ep - 0.5*(ESxin.stpp + ESxfin->stpp) * (*thick) ;
      ESxaux.ep = (ESxin.ep + ESxfin->ep) / 2.;
      ESxaux.stpp = getSP(&pMat->comp, SPTArray, ESxaux.ep);
    }

    ESxaux.x = ESxin.x + *thick / 2.;
    
    //Reallocating ESxA
    *ESxA=(ESxType*)realloc(*ESxA,((*pFpos)+2)*sizeof(ESxType));
    if (*ESxA==NULL){fprintf(stderr,"\n Error allocating memory (ESxA)\n");exit(1);}
       
    (*ESxA)[*pFpos]=ESxaux;
    (*pFpos)++;
    
    (*ESxA)[*pFpos]= *ESxfin; 
    (*pFpos)++;

    
    /**@todo The reallocation of ESxArray could be done inteligently in increments larger than 2, thus sacrificing some memory to speed up the process. Some extra variable would also be needed to remember the actual dimension of the matrix and comparing it with pFpos.*/
//      printf("\nDEBUG :  x[%d]=%.3le (E=%.2lf)\tx[%d]=%.3le (E=%.2lf) ",*pFpos-2,ESxaux.x,ESxaux.ep,*pFpos-1, ESxfin->x,ESxfin->ep);
//     printf("\nDEBUG : %d) xin=%.3le (E=%.2lf)\txfin=%.3le (E=%.2lf) ",*pFpos,ESxin.x,ESxin.ep,ESxfin->x,ESxfin->ep); 
  }
  else{ //split in two and call SSThick() for each part
//      if(auxdecis<0.97) printf("\nDEBUG :  R[%d]=A (%.4le)  x=%.4le",*pFpos,auxdecis,ESxfin->x);
//      else printf("\nDEBUG :  R[%d]=S",*pFpos);
    
    ESxm1.ep = (ESxin.ep + ESxfin->ep) / 2;
    ESxm1.stpp = getSP(&pMat->comp, SPTArray, ESxm1.ep);
    ESxm1.x = ESxin.x + (ESxin.ep - ESxm1.ep) / ESxin.stpp;
    //if (ESxm1.x > LayerThickness)ESxm1.x = LayerThickness; //note: this value of ESxm1.x is only used for the calc of auxdecis, but does not affect to the stopping and E calcs, so we don't need to update ESxm1.ep and ESxm1.stpp
    
    SSThick(pexp, pMat, SPTArray, MajAbsCoef, LayerThickness, ESxin, &ESxm1, pFpos, &tck1, ESxA);
    
    if (ESxm1.x < LayerThickness) {
      SSThick(pexp, pMat, SPTArray, MajAbsCoef, LayerThickness, ESxm1, ESxfin, pFpos, &tck2, ESxA);
      *thick = tck1 + tck2;
    }
    else *thick = tck1;
  }  
  return(0);
}



/**Performs Calculation of Transmission and Yields and then calculates the Penetration Integral.
 Uses Simpson integration method.

TotAbsCoefArray
XYldSums (O)is a matrix which contains, for each element, the total sum of XRayYlds for the whole sample.


@todo document this function

NOTE: This function uses the global vars maxK, minM & minL
*/
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 ilyr, itrans, ii, jj, ll,mm;
  LYR *plyr;
  CalibYld *AbsFac;   
  CalibYld *pSSYld, *pSSTrs;  
  const AbsCoef *AbsC;
  const FluorCKCoef *pFCKTrc;
  XrayYield *pResYld;
  const CalibYld *pXYld;
  SecResType *pSecRes;
  CalibYld *pXYldSum;
  double tckout,tempE,tempM,attfraction;
  int tempZ, FirstReg;
  const SIM_PARAM *psim;
  CalibYld tempTr0;   
  psim=&pexp->simpar;
  
  for(ilyr=1; ilyr<=NFoilUsed; ilyr++){
    plyr=&plyrarray[ilyr]; 
  
    /*Initialize (with 0 values) the SSTrsArray. (Simpson Transmission)*/
    plyr->SSTrsArray=(CalibYld*)calloc((plyr->NumOfTrc * plyr->FESxlen),sizeof(CalibYld));
    if (plyr->SSTrsArray==NULL){fprintf(stderr,"\n Error allocating memory (lyr[layer].SSTrsArray)\n");exit(1);}
    /*Initialize (with 0 values) the SSYldArray. (Simpson Yield) */
    plyr->SSYldArray=(CalibYld*)calloc((plyr->NumOfTrc * plyr->FESxlen),sizeof(CalibYld));
    if (plyr->SSYldArray==NULL){fprintf(stderr,"\n Error allocating memory (lyr[layer].SSYldArray)\n");exit(1);}
    
          
    /*Calculate the transmission array (all sublayers x all elements) FOR THE GIVEN LAYER ONLY*/
  
    for (jj = 0; jj < plyr->NumOfTrc ; jj++) {
      if (plyr->TrcUse[jj]) {
        
        AbsFac= &(plyr->SecResArray[jj].SR.AbsFact);
        tempZ=plyr->TrcAtnum[jj];
        
        for (ii = 0; ii < plyr->FESxlen; ii++) { 
          tckout=plyr->ESxArray[ii].x * pexp->cosFac;  //effective thicknes to the surface of this layer from the sublayer ii  
          pSSTrs=&(plyr->SSTrsArray[ii + plyr->FESxlen * jj]);  //points to the "transmission matrix element"
          
          if (tempZ < maxK) {
            for (ll = 1; ll <= 3; ll++) {
              if (AbsFac->K_[ll] > 0)
                pSSTrs->K_[ll] = exp(-AbsFac->K_[ll] * tckout);
            }
          }
          if (tempZ > minM) {
            for (ll = 1; ll <= 3; ll++) {
              if (AbsFac->M_[ll] > 0)
                pSSTrs->M_[ll] = exp(-AbsFac->M_[ll] * tckout);
            }
          }
          if (tempZ > minL) {
            for (ll = 1; ll <= 3; ll++) {
              for (mm = 1; mm <= 3; mm++) {
                if (AbsFac->L_[ll][mm] > 0)
                  pSSTrs->L_[ll][mm] = exp(-AbsFac->L_[ll][mm] * tckout);
              }
            }
          }
        }
      }
    }
    
    /*Perform "Simpson" Calculation of yield. 
    (Originally this was done in a separated function called SimpSSYld()*/
    for (jj = 0; jj < plyr->NumOfTrc; jj++) {
      if (plyr->TrcUse[jj]) {
        
        tempZ=plyr->TrcAtnum[jj];
        Z2mass(tempZ, &tempM, 'n');
        pFCKTrc=&FCKCoefArray[tempZ];
        AbsC=&TotAbsCoefArray[tempZ];
        
        for (ii = 0; ii < plyr->FESxlen; ii++) {
          tempE=plyr->ESxArray[ii].ep;
          pSSYld=&( plyr->SSYldArray[ (ii + plyr->FESxlen * jj) ] );
          
          Xprod(AbsC, pFCKTrc, pexp->ion.Z, tempZ, tempM, tempE, pSSYld);
          
        }
      }
    }     
  
  }
  
  /*clears any previously used Sum*/
  
  /*Calculates the Penetration Integral*/
  for(ilyr=1; ilyr<=NFoilUsed; ilyr++){
    plyr=&plyrarray[ilyr]; 
    
    for (ii = 0; ii < plyr->NumOfTrc ; ii++) {
    
      if (plyr->TrcUse[ii]) {
        tempZ=plyr->TrcAtnum[ii];
        Z2mass(tempZ, &tempM, 'n');
        pFCKTrc=&FCKCoefArray[tempZ];
        AbsC=&TotAbsCoefArray[tempZ];
        pResYld=&plyr->ResYldArray[ii];
        pSecRes=&plyr->SecResArray[ii];
        AbsFac= &(plyr->SecResArray[ii].SR.AbsFact);
        pXYld=&XYldArray[tempZ].XYld;
        pXYldSum=&XYldSums[tempZ];
        attfraction=plyr->pFoil->comp.xn[ii];

//         printf("\n DEBUG: Calculating Penetration Integral for : %s in layer %d\n",pResYld->symb,ilyr) ;

        //pSSTrs and pSSYld are reused to hold now one-dimensional arrays of CalibYld structures containing the Transmission and yield values for the ii element in all sublayers. More precisely, they point to subsections (rows) of the 2-dim arrays SSTrsArray and SSYldArray that contain the mentioned info.
        FirstReg = (ii) * plyr->FESxlen;  
        pSSTrs=&plyr->SSTrsArray[FirstReg]; 
        pSSYld=&plyr->SSYldArray[FirstReg];
        
        // keep in mind that YFirstReg and TFirstReg (now combined in just FirstReg) have been shifted down by one to accomodate to the redefiniton of the SSTrsArray and Yield matrices (which start in 0 for both indexes) 
                
        //correction due to transmission in layers over the one we are considering
        tempTr0=plyr->SSTrsArray[FirstReg];  //the relevant lines become initialized at value=1
        //For K lines...
        if (tempZ < maxK) {
          for (ll = 1; ll <= 3; ll++) {
            tempTr0.K_[ll] *= Filter->Trans[tempZ].K_[ll]; //Apply transmission of filter
            for(itrans=ilyr-1 ; itrans>0 ; itrans--){ //apply transmission of every layer on top of current one
              tempTr0.K_[ll] *= plyrarray[itrans].Trans[tempZ].K_[ll];
            }
          }
        }
        //Same for M lines...
        if (tempZ > minM) {
          for (ll = 1; ll <= 3; ll++) {
            tempTr0.M_[ll] *= Filter->Trans[tempZ].M_[ll]; 
            for(itrans=ilyr-1 ; itrans>0 ; itrans--){
              tempTr0.M_[ll] *= plyrarray[itrans].Trans[tempZ].M_[ll];
            }
          }
        }
        //Same for L lines...
        if (tempZ > minL) {
          for (ll = 1; ll <= 3; ll++) {
            for (mm = 1; mm <= 3; mm++) {
              tempTr0.L_[ll][mm] *= Filter->Trans[tempZ].L_[ll][mm]; 
              for(itrans=ilyr-1 ; itrans>0 ; itrans--){
                tempTr0.L_[ll][mm] *= plyrarray[itrans].Trans[tempZ].L_[ll][mm];
              }
            }
          }
        }  
          
        PenInteg(tempZ, AbsFac, plyr->ESxArray, pSSYld,
                pSSTrs, &tempTr0, plyr->FESxlen, 
                plyr->NeedSFC[ii], psim->AllowXEqCalc, plyr->absolutepos, pexp->cosInc,
                &pResYld->XYld, &pSecRes->SR.SFCr, &pSecRes->SR.Xeq);
          
        deNormalize(pexp, AbsC ,pFCKTrc, tempZ, tempM, attfraction, pXYld, pResYld, pXYldSum);
        
      }
    }
  }

  
  return(0);
}

/**Calculates X-Ray production efficiency for a given beam and target, in cm2/(1e15at).
 * 
 * @param AbsC (I) pointer to Absorption coef structure for Z2
 * @param pFCK (I) pointer to Fluorescence and Coster-Kronig coefficients for Z2
 * @param Z1 (I) atomic number of the ion (note: currently only Z1=1 is supported)
 * @param Z2 (I) atomic number of the target
 * @param M2 (I) mass of the target atom
 * @param ener (I) energy of the incident ion, in keV
 * @param XYld (O) Pointer to return the X-ray production efficiency, in cm2/(1e15at)
 * @return Always 0
 *
 *@todo This function only works for proton beams!!
 */
int Xprod(const AbsCoef *AbsC, const FluorCKCoef *pFCK, atomicnumber Z1, atomicnumber Z2, double M2 , double ener, CalibYld *XYld)
{
  int i, j;
  SecXL XLaux;
  double auxsec;
  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 }};
  double barn_to_cm21e15;
  
  barn_to_cm21e15=1e-9;  // 1 barn/at = 1e-24 cm2/at= 1e-9 cm2/(1e15at)

  
  if(Z1!=1){fprintf(stderr,"\n Error: Xprod() currently supports only H ions\n");exit(1);}
    
  *XYld = CYldNul;
  switch (Z2) {

  case 10: ///@todo Check hardcoded Z limits 
  case 11:
  case 12:
  case 13:
  case 14:
    auxsec = PaulX(ener, Z2);  // barn
    XYld->K_[1] = auxsec * pFCK->k.K_[1] * barn_to_cm21e15;  // cm2/1e15at
    break;

  case 54:
  case 55:
  case 56:
  case 57:
  case 58:
  case 59:
    ReisX(AbsC,pFCK, ener, Z2, M2, XLaux);
    for (i = 1; i <= 3; i++) {
      for (j = 1; j <= 3; j++)
        XYld->L_[i][j] = XLaux[4 - i] * pFCK->k.L_[i][j] * barn_to_cm21e15;  // cm2/1e15at
    }
    break;

  default:
    if (Z2 >= 15 && Z2 <minL) { ///@todo Check hardcoded Z limits 
      auxsec = PaulX(ener, Z2);
      XYld->K_[1] = auxsec * pFCK->k.K_[1] * barn_to_cm21e15;  // cm2/1e15at
      XYld->K_[2] = auxsec * pFCK->k.K_[2] * barn_to_cm21e15;  // cm2/1e15at
    } 
    else if (Z2 >= minL && Z2 < maxK) {
      auxsec = PaulX(ener, Z2);
      ReisX(AbsC,pFCK, ener, Z2, M2, XLaux);
      for (i = 1; i <= 3; i++) {
        XYld->K_[i] = auxsec * pFCK->k.K_[i] * barn_to_cm21e15;  // cm2/1e15at
        for (j = 1; j <= 3; j++)
          XYld->L_[i][j] = XLaux[4 - i] * pFCK->k.L_[i][j] * barn_to_cm21e15;  // cm2/1e15at
      }
    }
    else if (Z2 >= minM && Z2 <= 99) {///@todo Check hardcoded Z limits 
      ReisX(AbsC,pFCK, ener, Z2, M2, XLaux);
      for (i = 1; i <= 3; i++) {
        for (j = 1; j <= 3; j++)
          XYld->L_[i][j] = XLaux[4 - i] * pFCK->k.L_[i][j] * barn_to_cm21e15;  // cm2/1e15at
      }
    }
    break;
  }
  return(0);
}

/**
 * 
 * @param ener 
 * @param z Atomic number of Target
 * @return X-ray cross-section, in barn
 *
 *@todo Document this function (description, units for ener & return and references)
 *
 *@todo This function only works for proton beams!!
 */
double PaulX(double ener, atomicnumber z)
{
  static double sc1 = -2.1717, sc2 = 10.8883, sc3 = 9.45875, sc4 = 0.975316,
    sc5 = 0.0165458, sc6 = 1.00859, sc7 = 0.0474606;

  static double ylim[4] = { 0., -0.86, -0.582, -0.037  };  //note, index 0 not used

  static double C[8][7] = {//note, indexes 0 not used
    { 0.,  0.       ,  0.        , 0.         ,  0.        ,  0.        ,  0.         },
    { 0.,  4.97159  , -0.0334597 ,  4.56721e-3, -4.17618e-6, -0.0152467 ,  9.11304e-4 },
    { 0.,  -3.89274 , -0.883283  , -0.0177272 ,  1.03168e-4,  0.824937  ,  0.0110194  },
    { 0.,  597612.0 ,  597919.0  ,  25220.1   , -37.2554   , -362099.0  , -13942.3    },
    { 0.,  0.107444 , -4.47727e-3,  1.30581e-4, -1.9793e-6 ,  1.00964e-8,  0.0        },
    { 0.,  0.113657 , -8.4197e-3 ,  2.40606e-4, -2.95528e-6,  1.26726e-8,  0.0        },
    { 0.,  0.0105593,  7.4928e-4 , -1.30255e-5,  9.19824e-8,  0.0       ,  0.0        },
    { 0., -0.0306492,  1.00377e-3, -1.68953e-5,  8.74937e-8,  0.0       ,  0.0        }};

  long i;
  double sc = 0.0;
  double MeV,f, pauly, paule, xx, auxpp1, auxpp2, auxpp3;
  double auxl[7];
  double b[8];

  
  MeV = ener*0.001;   // converting keV-->MeV
  paule = log10(MeV / (z * z));
  xx = paule / 1.15 + 2.22;
  pauly = PaulX_y(MeV, z);
  if (pauly < ylim[1])   // This correction is not in the PAul's paper
    sc = 0.0;
  if (pauly >= ylim[1] && pauly <= ylim[2])
    sc = sc1 - sc2 * pauly - sc3 * pauly * pauly;
  if (pauly > ylim[2] && pauly <= ylim[3])
    sc = sc4 - sc5 * cos(13.6 * (pauly + 0.393));
  if (pauly > ylim[3])
    sc = sc6 + sc7 * cos(6.23 * (pauly - 0.33));
  for (i = 1; i <= 7; i++) {
    b[i] = C[i][1] + C[i][2] * z + C[i][3] * z * z + C[i][4] * z * z * z;
    switch (i) {

    case 1:
    case 2:
    case 3:
      b[i] /= 1 + C[i][5] * z + C[i][6] * z * z;
      break;

    case 4:
    case 5:
      b[i] += C[i][5] * z * z * z * z;
      break;
    }
  }
  for (i = 3; i <= 6; i++)
    auxl[i] = PaulX_lp(xx, i);
  auxpp1 = paule - b[3];
  auxpp1 *= auxpp1;
  auxpp2 = b[4] * PaulX_lp(xx, 3);
  auxpp3 = b[5] * PaulX_lp(xx, 4) + b[6] * PaulX_lp(xx, 5) + b[7] * PaulX_lp(xx, 6);
  f = b[1] + b[2] * auxpp1 + auxpp2 + auxpp3;
  return (sc * exp(f * log(10.0) - 2.2 * log(z)));
}



/**
 * 
 * @param MeV Energy of the ion, in MeV
 * @param z Atomic number of Target
 * @return 
 *
 *@todo Document this function (description, units for ener & return and references)
 *
 *@todo This function only works for proton beams!!
 */
double PaulX_y(double MeV, atomicnumber z)
{
  static double teta1 = 0.313076, teta2 = 0.149231, teta3 = 5.54982e-5,
    teta4 = 3.7093e-6, teta5 = 0.166608;
  double auxz = z;
  double eta, theta, TEMP;

  TEMP = auxz - 0.3;
  eta = 40.0283 * MeV / (TEMP * TEMP);
  theta = teta1 + teta2 * auxz - teta3 * auxz * auxz + teta4 * auxz * auxz * auxz;
  theta /= 1 + teta5 * auxz;
  return (log10(2 * sqrt(eta) / theta));
}

/**
 * 
 * @param x 
 * @param p
 * @return 
 *
 *@todo Document this function (description, units for ener & return and references)
 */
double PaulX_lp(double x, long p)
{
  double prr = 0.0;
  double TEMP, TEMP1;

  switch (p) {

  case 3:
    prr = 0.5 * (5 * x * x * x - 3 * x);
    break;

  case 4:
    TEMP = x * x;
    prr = 0.125 * (35 * (TEMP * TEMP) - 30 * x * x + 3);
    break;

  case 5:
    TEMP = x * x;
    prr = 0.125 * (63 * x * (TEMP * TEMP) - 70 * x * x * x + 15 * x);
    break;

  case 6:
    TEMP = x * x;
    TEMP1 = x * x;
    prr = 0.0625 * (231 * x * x * (TEMP * TEMP) - 315 * (TEMP1 * TEMP1) +
        105 * x * x - 5);
    break;
  }
  return (prr);
}


/**Returns X-ray fluorescence cross-sections for L lines. See paper ??? @todo(cite M.Reis paper named ???)
 * 
 * @param AbsC (I) Pointer to Absorption coef structure for element with Z2
 * @param pFCK (I) Pointer to structure containing the Fluorescence and Coster-Kronig coefs
 * @param ener (I) Proton energy, in keV
 * @param z (I) Atomic number of Target
 * @param M2 (I) Atomic mass of Target (in amu)
 * @param sigmaXL (O) X-ray cross sections, in barn
 * @return Always 0
 *
 *@todo This function is only for proton beams!! 
 */
int ReisX(const AbsCoef *AbsC, const FluorCKCoef *pFCK, double ener, atomicnumber z, double M2, double *sigmaXL)
{  

  SecXL ion, ionuni;
  int i;
  double xi[4];
  int seted = 0;
  double TEMP;
  double enr;
  FwTipo tipo;
  
  int Ncam, Zpr, nium;
  double MeV, Hr, cau, a02, Z2b, z2, z4, miu, EnrSm1, Mpr, MprMeV, MprUAt, a2s, v1,Iabs, q0sb,
         theta, ksi, zeta, ys, ys2, d, mRs, ksir, xcb, zcbs, cbe, etam, ratnorm;
  
  MeV = ener * 1e-3;   //converting from keV to MeV
  Zpr = 1;   //Only protons!
  Mpr = 1.007;   //Only protons!
  Ncam = 2;   //L atomic layer
  
  Hr = 27.2116;            ///@todo What is Hr ???
  cau = 137.03604;         ///@todo What is cau ???
  a02 = 2.800283608e-21;   ///@todo What is a02 ???
  
  tipo = L;
  Z2b = (double)(z) - 4.15;
  z2 = Z2b * Z2b;
  z4 = z2 * z2;
  
  MprMeV = Mpr * MEVAMU;   //MEVAMU is defined at the beginning of the file
  MprUAt = Mpr / UATMASS ; //UATMASS is defined at the beginning of the file
  EnrSm1 = MeV / MprMeV;
  miu = Mpr / (1 + Mpr / M2) / UATMASS;   /** U.Atomicas de massa  **/
  a2s = (double)(Ncam * Ncam) / Z2b;
  TEMP = 1 + EnrSm1;
  TEMP = 1 / (TEMP * TEMP);
  v1 = cau * sqrt(1 - TEMP);   /** U.Atomicas  **/
  d = (double)Zpr * z / (miu * (v1 * v1));
  for (i = 1; i <= 3; i++) {
    enr=AbsC->enr[i+1];
    if (enr > 0) {
      seted = 1;
      Iabs = enr * 1e3;   // in eV 
      q0sb = Iabs / (Hr * v1);
      theta = (double)(Ncam * Ncam * 2) * Iabs / (z2 * Hr);
      ksi = 2 * v1 / (theta * (Z2b / (double)Ncam));
      TEMP = ReisX_gs(tipo, i, ksi);
      TEMP -= ReisX_hs(tipo, i, ksi, theta);
      zeta = 1. + (double)(Zpr * 2) * TEMP / (Z2b * theta);
    /** parametro de correccao de polarizacao-ligacao **/
      TEMP=z2/(cau*cau*ksi/zeta);
      if (i == 1)  ys = TEMP * 0.40  / (double)Ncam ;
      else  ys =  TEMP * 0.15 ;
      ys2 = ys * ys;
      mRs = sqrt(1. + 1.1 * ys2) + ys;
      ksir = sqrt(mRs) * ksi;
      xi[i] = 1. / ksir;
      /** correccao de defleccao Coulombiana  **/
      xcb = 2. * M_PI * d * q0sb * zeta;
      zcbs = sqrt(1. - 4. * zeta / (theta * miu * ksi * ksi * mRs));
      switch (i) {
        case 1:
          nium = 9;
          break;
        case 2:
        case 3:
          nium = 11;
          break;
        case 4:
        case 5:
          nium = 13;
          break;
      }
      TEMP = xcb / (zcbs * (1 + zcbs));
      cbe = (double)nium * ReisX_En(TEMP, nium + 1);
      TEMP = theta * ksir / ((double)(Ncam * 2));
      etam = TEMP * TEMP;
      TEMP = (double)Zpr / Z2b;
      ratnorm = 8. * M_PI * (a02 / (etam * theta)) * (TEMP * TEMP) * cbe / BARNTOM2;
      ionuni[i] = ReisX_polisec(i, ksir, theta);
      ion[i] = ionuni[i] * ratnorm;
    }
  }
  sigmaXL[1] = ReisX_g(1, z, xi[1]) * pFCK->w[2] * ion[1];
  sigmaXL[2] = ReisX_g(2, z, xi[2]) * pFCK->w[3] * (pFCK->ck[0] * ion[1] + ion[2]);
  sigmaXL[3] = ReisX_g(3, z, xi[3]) * pFCK->w[4] * 
               ((pFCK->ck[1] + pFCK->ck[0] * pFCK->ck[2]) * ion[1] + pFCK->ck[2] * ion[2] + ion[3]);
  
  return(0);
}

/**Auxiliar function to ReisX()
 * 
 * @param T 
 * @param SS 
 * @param ksis 
 * @return 
 * 
 *@todo check that the units of the hardcoded constants are appropriate 
 *
 *
 *@todo document this function
 */
double ReisX_gs(FwTipo T, int SS, double ksis)
{
  /** ref. Brandt and Lapicki,Phys.Rev.,A20(1979)465 **/
  static double gK[9] = {
    1.0, 1.0, 9.0, 31.0, 98.0, 12.0, 25.0, 4.2, 0.515
  };

  static double gL1[9] = {
    1.0, 1.0, 9.0, 31.0, 49.0, 162.0, 63.0, 18.0, 1.97
  };

  static double gL23[10] = {
    1.0, 1.0, 10.0, 45.0, 102.0, 331.0, 6.7, 58.0, 7.8, 0.888
  };

  double ksis2, ksis3, ksis4, ksis5, ksis6, ksis7, ksis8;
  double gsaux = 0.0;
  double gsaux2;

  ksis2 = ksis * ksis;
  ksis3 = ksis * ksis2;
  ksis4 = ksis2 * ksis2;
  ksis5 = ksis2 * ksis3;
  ksis6 = ksis3 * ksis3;
  ksis7 = ksis3 * ksis4;
  ksis8 = ksis4 * ksis4;
  switch (T) {

  case K:
    gsaux = gK[1] + gK[2] * ksis + gK[3] * ksis2 + gK[4] * ksis3;
    gsaux += gK[5] * ksis4 + gK[6] * ksis5 + gK[7] * ksis6 + gK[8] * ksis7;
    gsaux2 = 1 + ksis;
    gsaux2 *= gsaux2 * gsaux2;
    gsaux2 *= gsaux2 * gsaux2;
    gsaux /= gsaux2;
    break;

  case L:
    if (SS == 1) {
      gsaux = gL1[1] + gL1[2] * ksis + gL1[3] * ksis2 + gL1[4] * ksis3;
      gsaux += gL1[5] * ksis4 + gL1[6] * ksis5 + gL1[7] * ksis6 + gL1[8] * ksis7;
      gsaux2 = 1 + ksis;
      gsaux2 *= gsaux2 * gsaux2;
      gsaux2 *= gsaux2 * gsaux2;
      gsaux /= gsaux2;
    } else {
      gsaux = gL23[1] + gL23[2] * ksis + gL23[3] * ksis2 + gL23[4] * ksis3;
      gsaux += gL23[5] * ksis4 + gL23[6] * ksis5 + gL23[7] * ksis6 +
         gL23[8] * ksis7;
      gsaux += gL23[9] * ksis8;
      gsaux2 = 1 + ksis;
      gsaux2 *= gsaux2 * gsaux2;
      gsaux2 *= gsaux2 * gsaux2;
      gsaux2 *= 1 + ksis;
      gsaux /= gsaux2;
    }
    break;
  default:break;
  }
  return (gsaux);
} 


/** Auxiliar function to ReisX()
 * 
 * @param T 
 * @param SS 
 * @param ksih 
 * @param thet 
 * @return 
 * 
 *@todo check that the units of the hardcoded constants are appropriate 
 *
 *
 *@todo document this function
 */
double ReisX_hs(FwTipo T, int SS, double ksih, double thet)
{
  static double IhsC[6] = { 0.031, 0.031, 0.210, 0.005, -0.069, 0.324  };
  double x, x12, ksih3;
  double Ihs = 0.0, cs = 0.0;
  int n2hs = 0;

  switch (T) {

  case K:
    n2hs = 1;
    cs = 3.0 / 2.;
    break;

  case L:
    n2hs = 2;
    if (SS == 1) cs = 3.0 / 2.;
    else  cs = 5.0 / 4.;
    break;
    
  default:break;  
  }
  ksih3 = ksih * ksih * ksih;
  x = cs * n2hs / ksih;
  x12 = sqrt(x);
  if (x > 0 && x <= 0.035)
    Ihs = 3 * M_PI / 4 * log(1 / (x * x) - 1);
  if (x > 0.035 && x <= 3.1) {
    Ihs = IhsC[1] + IhsC[2] * x12 + IhsC[3] * x + IhsC[4] * x * x12 + IhsC[5] * x * x;
    Ihs = exp(-2 * x) / Ihs;
  }
  if (x > 3.1 && x < 11)
    Ihs = 2 * exp(-2 * x) / exp(1.6 * log(x));
    
  return(n2hs * 2 * Ihs / (thet * ksih3));
} 


/**Auxiliar function to ReisX()
 * 
 * @param z 
 * @param niu 
 * @return 
 * 
 *@todo check that the units of the hardcoded constants are appropriate 
 *
 *
 *@todo document this function
 */
double ReisX_En(double z, int niu)
{
  /** aproximation for niu>> **/
  /** vide Abramowitz & Stegun , Dover 1Ed. 1965 **/
  double Eaux, z2, z3, niu2, aux, aux2, aux4, aux6;

  z2 = z * z;
  z3 = z2 * z; 
  niu2 = niu * niu;
  aux = z + niu;
  aux2 = aux * aux;
  aux4 = aux2 * aux2;
  aux6 = aux4 * aux2;
  Eaux = niu * (6 * z2 - 8 * niu * z + niu2) / aux6;
  Eaux++;
  Eaux = exp(-z) * Eaux / aux;
  return(Eaux);
}

/**Auxiliar function to ReisX()
 * 
 * @param ssind 
 * @param kz 
 * @param tz 
 * @return 
 * 
 * 
 *@todo check that the units of the hardcoded constants are appropriate 
 *
 *
 *@todo document this function
 */
double ReisX_polisec(int ssind, double kz, double tz)
{
  double Result = 0.0;

  static double fc1[8] = {
     737.658767, -5422.24598, 17070.7835, -29572.7846, 
     30490.2736, -18719.3387, 6343.25887, -916.108807
  };

  static double fc2[8] = {
    -1513.63296, 5152.69169, -7347.41315, 5748.20489,
    -2665.54100, 733.610039, -111.064039, 7.14045527
  };

  static double fc3[8] = {
     13.1119073, -41.7263393, 97.1163461, -100.159645, 
     60.9662790, -21.6315896, 4.12048128, -0.325361134
  };

  static double fc4[8] = {
     18.9172599, -59.1559460, 126.801210, -127.067875, 
     74.6217068, -25.5046721, 4.68627281, -0.357639954
  };

  double p1;
  double x = 0.0;
  double coef[8];
  int i;
  
  switch (ssind) {

  case 1:
    x = 1 / sqrt(kz * exp(0.3 * log(tz)));
    if (x < 1.35)
      memcpy(coef, fc1, sizeof(double) * 8);
    else
      memcpy(coef, fc2, sizeof(double) * 8);
    break;

  case 2:
    x = 1 / sqrt(kz * exp(0.2 * log(tz)));
    memcpy(coef, fc3, sizeof(double) * 8);
    break;

  case 3:
    x = 1 / sqrt(kz * exp(0.32 * log(tz)));
    memcpy(coef, fc4, sizeof(double) * 8);
    break;
  }
  
  for (p1 = coef[7], i = 6; i >=0; i--)  p1 = x * p1 + coef[i]; //polynomial calculation
  
  switch (ssind) {

  case 1:
    Result = exp(-p1) / exp(3.8 * log(tz));
    break;

  case 2:
    Result = exp(-p1) / exp(2.5 * log(tz));
    break;

  case 3:
    Result = exp(-p1) / exp(6.5 * log(tz));
    break;
  }
  return(Result);
}

/**Auxiliar function to ReisX()
 * 
 * @param ss 
 * @param Zg 
 * @param xi 
 * @return 
 * 
 *@todo check that the units of the hardcoded constants are appropriate 
 *
 *
 *@todo document this function
 */
double ReisX_g(int ss, atomicnumber Zg, double xi)
{
  static double g3pol[2][8] = {
    { 4.980374404, -18.84795179, 37.53135546, -39.83029409, 23.94060788,
      -8.048533664, 1.397129084, -0.09706728869 },
    { -1.859600742, 12.30009972, -21.67920610, 20.38505540, -11.17032345,
      3.604521250, -0.6347055945, 0.04676377726 }
  };

  static double g12pol[2][8] = {
    { 0.7800031227, 0.5779760376, 0.01209015664, 1.382072250, -2.555731314,
      1.534786106, -0.3844825563, 0.03457310793 },
    { -4.910259919, 27.43184219, -50.44519949, 48.46725025, -26.40954176,
      8.230146858, -1.365807669, 0.09332547710 }
  };

  int i;
  int polsel = 0;
  double gaux = 0.0;

  if (xi < 0.6 || xi > 2.6)
    return(1.0);
  else {
    if (ss == 3) {
      if (Zg < 47)  return(1.0);
      
      if (Zg < 64) polsel = 0;
      else polsel = 1;        
      
      for (gaux = g3pol[polsel][7],i = 6; i >= 0; i--)  gaux = gaux * xi + g3pol[polsel][i];
      return(gaux);
    }
    if (ss == 1) {
      if (Zg > 61 && Zg < 72)  polsel = 0;
      else return(1.0);
    }
    if (ss == 2) {
      if (Zg > 54 && Zg < 75)   polsel = 1;
      else return(1.0);
    }
    
    for (gaux = g12pol[polsel][7], i = 6; i >= 0; i--)  gaux = gaux * xi + g12pol[polsel][i];
    return(gaux);
  }
}


/** Calculates a normalized "Penetration Integral" for a given X-ray emission line of a given element. 
See "Reis et al. NIM B 109/110 (1996) 134-138"
See section 1.4 of "Correccoes de Fluorescencia Secundaria em PIXE" Master Thesis by M.A. Reis

It calculates the integral of eq. (4) of the cited paper except that it does not include the normalization constant "sigmaX_ij(Ep)"

  * 
  * @param atnumb (I) Atomic number of the target element
  * @param AbsFac (I) Absorption coefficients structure for the target element
  * @param ESA  (I) Pointer to Array of sublayers in the current layer.
  * @param YldA (I) Array of X-ray production cross section for the target element in each sublayer
  * @param TrsA (I) Array of Transmission coefs for the target element in each sublayer
  * @param pTrs0 (I) Array of Transmission factors due tolayers over the one we are considering.
  * @param FExlen (I) Number of sublayers in this layer
  * @param NeedSFC (I) Flag indicating whether secondary fluorescence should be calculated (1) or not (0)
  * @param AllowXEqCalc (I) Flag indicating whether "equivalent thickness" should be calculated (1) or not (0)
  * @param x0 (I) absolute in-going path length, in 1e15at/cm2, from sample surface (Only relevant if AllowXEqCalc=1)
  * @param CosInc Cosine of the incident angle (Only relevant if AllowXEqCalc=1)
  * @param XYld (O) Penetration integral results.
  * @param XSFCr (O) Correction factor applied due to secondary fluorescence  (Only relevant if NeedSFC=1)
  * @param XYldxmed (O) Equivalent thickness results (Only relevant if AllowXEqCalc=1)

*/
void PenInteg(atomicnumber atnumb, const CalibYld *AbsFac, 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)
{
  int i, j, l, limit;
  double auxfa, auxfb, auxfi, diffl;
  CalibYld auxxmed, auxy2;
  const CalibYld *pTrs1, *pYld1, *pSFC1, *pTrs2, *pYld2, *pSFC2, *pTrs3, *pYld3, *pSFC3;
  
  const ESxType *pEF1, *pEF2, *pEF3;
  
/*  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 }};*/
        
  /*These variables need to be 0-initialized because they will be used for sumation*/

  *XYld=CYldNul;
  
  if(AllowXEqCalc) auxxmed = CYldNul;
  
  
  if (NeedSFC) {
    fprintf(stderr,"\n Error: Secondary Fluorescence not implemented yet\n"); exit(1);
    /*Note: when implementing Secondary fluorescence correction, 
    SFC1, SFC2 and SFC3 should be read here*/
    auxy2 = CYldNul;
    }
  else {
    pSFC1=&CYldNul;
    pSFC2=&CYldNul;
    pSFC3=&CYldNul;
  }
  
  pEF1=&ESA[0];
  pTrs1=&TrsA[0];
  pYld1=&YldA[0];
  
  /*Note: The following loop starts in 1 and ends 1 before the max to make possible
    the definition of Trs1, Trs2,Trs3 and so on. Also note that its step is 2*/
  limit=FExlen-1;
  for(i=1 ; i<limit ; i+=2){    
    
    pEF2=&ESA[i];
    pTrs2=&TrsA[i];
    pYld2=&YldA[i];

    pEF3=&ESA[i+1];
    pTrs3=&TrsA[i +1];
    pYld3=&YldA[i +1];
    
    if (atnumb < maxK) {
      for (j = 1; j <= 3; j++) {
        if (AbsFac->K_[j] > 0. && pTrs2->K_[j] > 1e-5) {
          
          auxfa = (pTrs1->K_[j] * pYld1->K_[j] + pSFC1->K_[j]) / pEF1->stpp;
          auxfi = (pTrs2->K_[j] * pYld2->K_[j] + pSFC2->K_[j]) / pEF2->stpp;
          auxfb = (pTrs3->K_[j] * pYld3->K_[j] + pSFC3->K_[j]) / pEF3->stpp;
          XYld->K_[j] -= pTrs0->K_[j] * Simps(pEF1->ep, pEF3->ep, auxfa, auxfi, auxfb);
          #if CPIXEVERBOSITY > 1
          if(j==1) printf("\n Z=%d   XYld=%le  (XYld/Dx=%le)   x1=%le    ",
                          atnumb, pTrs0->K_[j] * 
                          Simps(pEF1->ep, pEF3->ep, auxfa, auxfi, auxfb),
                          Simps(pEF1->ep, pEF3->ep, auxfa, auxfi, auxfb)/(pEF1->x - pEF3->x),
                          pEF1->x);
          #endif
          ///@todo IMPORTANT: check here the use of cosInc in relation with the definition of x0
          if (AllowXEqCalc) {
            auxfa = (pTrs1->K_[j] * pYld1->K_[j] + pSFC1->K_[j]) * (pEF1->x + x0) * CosInc / pEF1->stpp;
            auxfi = (pTrs2->K_[j] * pYld2->K_[j] + pSFC2->K_[j]) * (pEF2->x + x0) * CosInc / pEF2->stpp;
            auxfb = (pTrs3->K_[j] * pYld3->K_[j] + pSFC3->K_[j]) * (pEF3->x + x0) * CosInc / pEF3->stpp;
            auxxmed.K_[j] -= pTrs0->K_[j] * Simps(pEF1->ep, pEF3->ep, auxfa, auxfi, auxfb);
          }
          if (NeedSFC){
            auxfa = pSFC1->K_[j] / pEF1->stpp;
            auxfi = pSFC2->K_[j] / pEF2->stpp;
            auxfb = pSFC3->K_[j] / pEF3->stpp;
            auxy2.K_[j] -= pTrs0->K_[j] * Simps(pEF1->ep, pEF3->ep, auxfa, auxfi, auxfb);
          }
        }
      }
    }
    if (atnumb > minM) {
      for (j = 1; j <= 3; j++) {
        if (AbsFac->M_[j] > 0. && pTrs2->M_[j] > 1e-5) {
          auxfa = (pTrs1->M_[j] * pYld1->M_[j] + pSFC1->M_[j]) / pEF1->stpp;
          auxfi = (pTrs2->M_[j] * pYld2->M_[j] + pSFC2->M_[j]) / pEF2->stpp;
          auxfb = (pTrs3->M_[j] * pYld3->M_[j] + pSFC3->M_[j]) / pEF3->stpp;
          XYld->M_[j] -= pTrs0->M_[j] * Simps(pEF1->ep, pEF3->ep, auxfa, auxfi, auxfb);
          if (AllowXEqCalc) {
            auxfa = (pTrs1->M_[j] * pYld1->M_[j] + pSFC1->M_[j]) * (pEF1->x + x0) * CosInc / pEF1->stpp;
            auxfi = (pTrs2->M_[j] * pYld2->M_[j] + pSFC2->M_[j]) * (pEF2->x + x0) * CosInc / pEF2->stpp;
            auxfb = (pTrs3->M_[j] * pYld3->M_[j] + pSFC3->M_[j]) * (pEF3->x + x0) * CosInc / pEF3->stpp;
            auxxmed.M_[j] -= pTrs0->M_[j] * Simps(pEF1->ep, pEF3->ep, auxfa, auxfi, auxfb);
          }
          if (NeedSFC){
            auxfa = pSFC1->M_[j] / pEF1->stpp;
            auxfb = pSFC2->M_[j] / pEF2->stpp;
            auxy2.M_[j] -= pTrs0->M_[j] * Simps(pEF1->ep, pEF3->ep, auxfa, auxfi, auxfb);
          }
        }
      }
    }
    if (atnumb > minL) {
      for (j = 1; j <= 3; j++) {
        for (l = 1; l <= 3; l++) {
          if (AbsFac->L_[j][l] > 0. && pTrs2->L_[j][l] > 1e-5) {
            auxfa = (pTrs1->L_[j][l] * pYld1->L_[j][l] + pSFC1->L_[j][l]) / pEF1->stpp;
            auxfi = (pTrs2->L_[j][l] * pYld2->L_[j][l] + pSFC2->L_[j][l]) / pEF2->stpp;
            auxfb = (pTrs3->L_[j][l] * pYld3->L_[j][l] + pSFC3->L_[j][l]) / pEF3->stpp;
            XYld->L_[j][l] -= pTrs0->L_[j][l] * Simps(pEF1->ep, pEF3->ep, auxfa, auxfi, auxfb);
            if (AllowXEqCalc) {
              auxfa = (pTrs1->L_[j][l] * pYld1->L_[j][l] + pSFC1->L_[j][l]) * (pEF1->x + x0) * CosInc / pEF1->stpp;
              auxfi = (pTrs2->L_[j][l] * pYld2->L_[j][l] + pSFC2->L_[j][l]) * (pEF2->x + x0) * CosInc / pEF2->stpp;
              auxfb = (pTrs3->L_[j][l] * pYld3->L_[j][l] + pSFC3->L_[j][l]) * (pEF3->x + x0) * CosInc / pEF3->stpp;
              auxxmed.L_[j][l] -= pTrs0->L_[j][l] * Simps(pEF1->ep, pEF3->ep, auxfa, auxfi, auxfb);
            }
            if (NeedSFC){
              auxfa = pSFC1->L_[j][l] / pEF1->stpp;
              auxfi = pSFC2->L_[j][l] / pEF2->stpp;
              auxfb = pSFC3->L_[j][l] / pEF3->stpp;
              diffl = 0.0;
              auxy2.L_[j][l] -= pTrs0->L_[j][l] * Simps(pEF1->ep, pEF3->ep, auxfa, auxfi, auxfb);
            }
          }
        }
      }
    }
    pYld1 = pYld3;
    pTrs1 = pTrs3;
    pEF1 = pEF3;
    pSFC1 = pSFC3;
  }
  
  if (NeedSFC) {
    for (i = 1; i <= 3; i++) {
      if (XYld->K_[i] > 0.)  XSFCr->K_[i] = auxy2.K_[i] / XYld->K_[i];
      if (XYld->M_[i] > 0.)  XSFCr->M_[i] = auxy2.M_[i] / XYld->M_[i];
      for (j = 1; j <= 3; j++) {
        if (XYld->L_[i][j] > 0.)  XSFCr->L_[i][j] = auxy2.L_[i][j] / XYld->L_[i][j];
      }
    }
  } 
  else  *XSFCr = CYldNul;
  
  if (AllowXEqCalc) {
    for (i = 1; i <= 3; i++) {
      if (XYld->K_[i] > 0.)  XYldxmed->K_[i] = auxxmed.K_[i] / XYld->K_[i];
      if (XYld->M_[i] > 0.)  XYldxmed->M_[i] = auxxmed.M_[i] / XYld->M_[i];
      for (j = 1; j <= 3; j++) {
        if (XYld->L_[i][j] > 0.)  XYldxmed->L_[i][j] = auxxmed.L_[i][j] / XYld->L_[i][j];
      }
    }
  }
  
  
}


/** Returns area of trapezoid, useful for Simpson Integration.
 *
 * See "Corrected Trapezoidal Rule" in 
 *    Conte & de Boor "Elementar Numerical Analysis" 3rd Ed. McGraw-Hill,1980,,p309
 *    ExpCTR:=(b-a)/2*(fa+fb)  
 *
 * @param a First point of the trapezoid base
 * @param b Second point of the Trapezoid base
 * @param fa 
 * @param fi 
 * @param fb 
 * @return 
 *
 *@todo IMPORTANT: Why the 6 and the 4??????? Has something to do with exponential function integration? If not, it should be rather: return(  (b-a)/4 * (fa+2*fi+fb)  ) 
 UPDATE: This has to do with the line in SSthick where "thick" is calculated.
 *
 *@todo document this function
 */
double Simps(double a, double b, double fa, double fi, double fb)
{
  
  return ((b - a) / 6. * (fa + 4. * fi + fb));
}



/*
useFilter is a flag
Filter is the filter struct
*/
/**Calculates the yield values using the penetration integral results, i.e. Applies all the factors such as solid angle, atomic fraction, charge,... etc

See eq. (1) of "Reis et al. NIM B 109/110 (1996) 134-138"

@todo This function does not support M-lines
@todo Note: This function is only valid for 10  Z2 in range [10,99]
 * 
 * @param pexp (I) Pointer to structure containing experimental parameters.
 * @param AbsC (I) Pointer to Absorption coef structure for element with Z2
 * @param pFCK (I) pointer to Fluorescence and Coster-Kronig coefficients for Z2
 * @param Z2 (I) atomic number of the target element
 * @param M2 (I) mass of the target element
 * @param attfraction (I) atomic fraction of the target element in the current layer
 * @param pXYld (I) Penetration integral values.
 * @param ResY (I+O) The calculated yields (in counts) for each layer will be returned in ResY->XYld
 * @param XYldSum (O) The total yield (sum for all layers processed up to this moment of ResY->XYld) is acumulated in XYldSum. 
 
 Note: deNormalize() returns both the partial and the total spectra through ResY->XYld and XYldSum, respectively.
 */
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)
{
  double aux1;
  long jj, ll;
  CalibYld TransFact, Xaux;
  const SIM_PARAM *psim;
  
  psim=&pexp->simpar;
  
  /*Set the Transfererence factor of the filter (=1 in case no filter)*/
  if (psim->useFilter){
    ///@todo do this part when the filters are supported
    fprintf(stderr,"\n Error: Filters not implemented yet\n"); exit(1);
//     for (jj = 1; jj <= 3; jj++) {
//       if (ResY->ener.K_[jj] > 0) TransFact.K_[jj] = FiltTrs(psim->Filter, ResY->ener.K_[jj]);
//       if (ResY->ener.M_[jj] > 0) TransFact.M_[jj] = FiltTrs(psim->Filter, ResY->ener.M_[jj]);
//     }
//     for (jj = 1; jj <= 3; jj++) {
//       for (ll = 1; ll <= 3; ll++) {
//         if (ResY->ener.L_[jj][ll] > 0) TransFact.L_[jj][ll] = FiltTrs(psim->Filter, ResY->ener.L_[jj][ll]);
//       }
//     }
  }
  else{
    for (jj = 1; jj <= 3; jj++) {
      TransFact.K_[jj] = 1.0;
      TransFact.M_[jj] = 1.0;
    }
    for (jj = 1; jj <= 3; jj++) {
      for (ll = 1; ll <= 3; ll++) TransFact.L_[jj][ll] = 1.0;
    }
  }
    
  Xprod(AbsC, pFCK, pexp->ion.Z, Z2, M2, psim->CalEner, &Xaux);
  
  switch (ResY->atnum) {

  case 10:///@todo Check hardcoded Z limits 
  case 11:
  case 12:
  case 13:
  case 14:
    if (Xaux.K_[1] > 0) {
      aux1 = ResY->XYld.K_[1] / Xaux.K_[1];
      ResY->XYld.K_[1] = pXYld->K_[1] * aux1 * psim->ColCharge * pexp->DetColFac * 
                         TransFact.K_[1] * psim->DTCC * attfraction;
      XYldSum->K_[1] +=ResY->XYld.K_[1];
    }
    break;

  case 54:
  case 55:
  case 56:
  case 57:
  case 58:
  case 59:
    for (jj = 1; jj <= 3; jj++) {
      for (ll = 1; ll <= 3; ll++) {
        if (Xaux.L_[jj][ll] > 0){
          ResY->XYld.L_[jj][ll] = pXYld->L_[jj][ll] * (ResY->XYld.L_[jj][ll] / Xaux.L_[jj][ll]) *
                                  psim->ColCharge * pexp->DetColFac * TransFact.L_[jj][ll] *
                                  psim->DTCC * attfraction;
          XYldSum->L_[jj][ll] += ResY->XYld.L_[jj][ll];
        }
      }
    }
    break;

  default:
    if (ResY->atnum >= 15 && ResY->atnum <minL) {///@todo Check hardcoded Z limits 
      for (jj = 1; jj <= 2; jj++) {
        if (Xaux.K_[jj] > 0) {
          aux1 = ResY->XYld.K_[jj] / Xaux.K_[jj];
          ResY->XYld.K_[jj] = pXYld->K_[jj] * aux1 * psim->ColCharge *
                              pexp->DetColFac * TransFact.K_[jj] * psim->DTCC * attfraction;
          XYldSum->K_[jj] +=ResY->XYld.K_[jj];
        }
      }
    } 
    else if (ResY->atnum >= minL && ResY->atnum < maxK) {
      for (jj = 1; jj <= 3; jj++) {
        if (Xaux.K_[jj] > 0) {
          aux1 = ResY->XYld.K_[jj] / Xaux.K_[jj];
          ResY->XYld.K_[jj] = pXYld->K_[jj] * aux1 * psim->ColCharge * 
                              pexp->DetColFac * TransFact.K_[jj] * psim->DTCC * attfraction;
          XYldSum->K_[jj] +=ResY->XYld.K_[jj];
        }
      }
      for (jj = 1; jj <= 3; jj++) {
        for (ll = 1; ll <= 3; ll++) {
          if (Xaux.L_[jj][ll] > 0){
            ResY->XYld.L_[jj][ll] = pXYld->L_[jj][ll] * (ResY->XYld.L_[jj][ll] / Xaux.L_[jj][ll]) *
                                    psim->ColCharge * pexp->DetColFac * TransFact.L_[jj][ll] * 
                                    psim->DTCC * attfraction;
            XYldSum->L_[jj][ll] +=ResY->XYld.L_[jj][ll];
          }
        }
      }
    } 
    else if (ResY->atnum >= minM && ResY->atnum <= 99) {///@todo Check hardcoded Z limits 
      for (jj = 1; jj <= 3; jj++) {
        for (ll = 1; ll <= 3; ll++) {
          if (Xaux.L_[jj][ll] > 0){
            ResY->XYld.L_[jj][ll] = pXYld->L_[jj][ll] * (ResY->XYld.L_[jj][ll] / Xaux.L_[jj][ll]) *
                                    psim->ColCharge * pexp->DetColFac * TransFact.L_[jj][ll] * 
                                     psim->DTCC * attfraction;
            XYldSum->L_[jj][ll] +=ResY->XYld.L_[jj][ll];
          }
        }
      }
      /**@todo Support also M lines here*/
    }
    break;
  }
}








/** Frees all memory allocated for the filter structure.
 * 
 * @param Filter (I+O) Pointer to filter structure.
 */
void freeFilter(FILTER *Filter)
{
  int i;
  foil *pFoil;
  
  for(i=0 ; i < Filter->nlyr ; i++){
    pFoil=&Filter->foil[i];
    safefree((void*)(&pFoil->comp.elem));
    safefree((void*)(&pFoil->comp.X));
    safefree((void*)(&pFoil->comp.xn));
    safefree((void*)(&pFoil->comp.w));
    
  }
  safefree((void*)(&Filter->FilterElems));
  safefree((void*)(&Filter->foil));
  safefree((void*)(&Filter->Trans));
    
}

/**Frees the space  allocated for arrays that are to be reused in subsequent iterations
 * 
 * @param NFoils (I)  Number of target foils
 * @param plyrarray (I+O)  Pointer to array of layers. Dim=NFoils+1
 * @param MatArray (I+O)  Pointer to array of foils. Dim=NFoils
 * @param XYldSums (I+O)  Pointer to the Output Array of summed X-Ray Yields
 */
void freeReusable(int NFoils, LYR **plyrarray, foil **MatArray, CalibYld **XYldSums)
{
  int i;
  LYR *plyr;
  foil *pMat;
  
  /*Clean the layer array and the Sample Array as well as all its members which have been allocated*/
  for(i=1;i<=NFoils;i++){
    plyr=&(*plyrarray)[i];
    safefree((void*)(&plyr->TrcAtnum));
    safefree((void*)(&plyr->Trans));
//     safefree((void*)(&plyr->AreasArray));
    safefree((void*)(&plyr->ResYldArray));
    safefree((void*)(&plyr->SecResArray));
    safefree((void*)(&plyr->TrcUse));
    safefree((void*)(&plyr->NeedSFC));
    safefree((void*)(&plyr->ESxArray));
    safefree((void*)(&plyr->SSTrsArray));
    safefree((void*)(&plyr->SSYldArray));

    
    pMat=&(*MatArray)[i-1];
    safefree((void*)(&pMat->comp.elem));
    safefree((void*)(&pMat->comp.X));
    safefree((void*)(&pMat->comp.xn));
    safefree((void*)(&pMat->comp.w));
    
  }
   safefree((void*)plyrarray);

   safefree((void*)MatArray);
   
   safefree((void*)XYldSums);
}


/**Frees memory only if *ptr is not NULL
 * 
 * @param ptr (I+O) Pointer to array to be freed
 */
void safefree(void **ptr){
  if(*ptr!=NULL){
    //printf("\nDEBUG: freing %x ...",*ptr);fflush(stdout);
    free(*ptr);
    //printf(" done (%x)",*ptr);fflush(stdout);
    *ptr=NULL;
  }
}


/**Prints a CalibYld structure on a given stream
 * 
 * @param f (I+O) Stream to which the structure is printed 
 * @param pCYld (I) Structure to print
 */
void fprintCALIBYLD(FILE *f,const CalibYld *pCYld)
{ 
  int j,l;
  for (j = 1; j <= 3; j++) fprintf(f,"\t%le",pCYld->K_[j]);
  for (j = 1; j <= 3; j++) { 
    fprintf(f,"\n");
    for (l = 1; l <= 3; l++)fprintf(f,"\t%le",pCYld->L_[j][l]);
  }
  fprintf(f,"\n");
  for (j = 1; j <= 3; j++) fprintf(f,"\t%le",pCYld->M_[j]);
  
  fprintf(f,"\n");
}

/**Reads a CalibYld structure from a given stream. The expected format is that as printed by fprintCALIBYLD
 * 
 * @param f (I) Stream from which the structure is read 
 * @param pCYld (O) Structure to store the read values
 */
void fscanCALIBYLD(FILE *f, CalibYld *pCYld)
{ 
  int j,l;
  for (j = 1; j <= 3; j++) fscanf(f,"%lf",&pCYld->K_[j]);
  for (j = 1; j <= 3; j++) { 
    for (l = 1; l <= 3; l++)fscanf(f,"%lf",&pCYld->L_[j][l]);
  }
  for (j = 1; j <= 3; j++) fscanf(f,"%lf",&pCYld->M_[j]);
  
}

---