libcpixe/libcpixe.c

Go to the documentation of this file.

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

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.

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






#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include "libcpixe.h"
#include "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*/


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;
  pexppar->ioncharge=-1; //initialize the ion charge to an impossible value
  //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
    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. */
        arg=strtok( NULL,"\n\t \r");
        strncpy(pexppar->ion.symbol,arg,3);
        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. */
        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,"IONCHARGE")){
        /*Syntax:
        IONCHARGE charge   <---"charge" is the beam ion charge in electron charge units.*/
        arg=strtok( NULL,"\n\t \r");
        pexppar->ioncharge=abs(atof(arg));
        /*Note that this is an optional command, no flag is changed. Default value=1 */
      }
      else if(!strcmp(command,"KEV")){
        /*Syntax:
        KEV Ebeam   <---"Ebeam" is Beam energy in 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*/
        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. */
        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 */
        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
        (Use this if using DATTPIXE-style calibration files... [DEPRECATED] )*/
        arg=strtok( NULL,"\n\t \r");
        pexppar->DetColFac=atof(arg);
        flag[7]=1;  //note that this flag can also be raised true (but 2) with SOLIDANGLE
        fprintf(stderr,"\n WARNING: Deprecated command '%s'. It may trigger bugs.\n Press <INTRO> to continue ...***\n",command);
        getchar();
      }
      else if(!strcmp(command,"SOLIDANGLE")){
        /*Syntax:
        SOLIDANGLE factor   <---Detector Solid Angle in msr [DEPRECATED]
        (Use this if using Efficiency-style calibration files [Recomended])*/
        arg=strtok( NULL,"\n\t \r");
        pexppar->SAngFract=atof(arg)/(4000.*M_PI);
        flag[7]=2;  //note that this flag be raised true with DETCOLFAC, SOLIDANGLE, SOLIDANGLEMSR or SOLIDANGLEFRAC
        fprintf(stderr,"\n WARNING: Deprecated command '%s'. Use SOLIDANGLEFRAC or SOLIDANGLEMSR instead.\n Press <INTRO> to continue ...***\n",command);
        getchar();
      }
      else if(!strcmp(command,"SOLIDANGLEMSR")){
        /*Syntax:
        SOLIDANGLEMSR msr   <---Detector Solid Angle in msr
        (Use this if using Efficiency-style calibration files [Recomended])*/
        arg=strtok( NULL,"\n\t \r");
        pexppar->SAngFract=atof(arg)/(4000.*M_PI);
        flag[7]=3;  
      }
      else if(!strcmp(command,"SOLIDANGLEFRAC")){
        /*Syntax:
        SOLIDANGLEFRAC factor   <---Detector Solid Angle fraction (solid angle/4PI)
        (Use this if using Efficiency-style calibration files [Recomended])*/
        arg=strtok( NULL,"\n\t \r");
        pexppar->SAngFract=atof(arg);
        flag[7]=4;  //note that this flag be raised true with DETCOLFAC, SOLIDANGLE, SOLIDANGLEMSR or SOLIDANGLEFRAC
      }
      else if(!strcmp(command,"CHARGE")){
        /*Syntax:
        CHARGE ColCharge   <---Collected charge, in uC */
        arg=strtok( NULL,"\n\t \r");
        pexppar->simpar.ColCharge=atof(arg);
        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);
        flag[9]=1;
      }
      else if(!strcmp(command,"ALLOWSXFCORR")){
        /*Syntax:
        ALLOWSXFCORR flag   <---Flag allowing for calculation of Sec. Fluorescence  */
        arg=strtok( NULL,"\n\t \r"); 
        pexppar->simpar.AllowSXFCorr=atoi(arg);
        flag[10]=1;
      }
      else if(!strcmp(command,"ALLOWXEQCALC")){
        /*Syntax:
        ALLOWXEQCALC flag   <---Deat time correction (Integrated"Charge"/IntegratedCounts) */
        arg=strtok( NULL,"\n\t \r"); 
        pexppar->simpar.AllowXEqCalc=atoi(arg);
        flag[11]=1;
      }
      else if(!strcmp(command,"CALIBFILE") || !strcmp(command,"DTCALIBFILE")){
        /*Syntax:
        DTCALIBFILE Calfilename     <---Name of the calibration  (XrayYld) file
        (Use this if using DATTPIXE-style calibration files [DEPRECATED])
        Note CALIBFILE is, at this moment an alias for DTCALIBFILE, but this will change to EFFCALIBFILE in a future
        */
        arg=strtok( NULL,"\n\t \r");
        if(!(pExtraInfo->CalibFileNm=(char*)calloc(LONGSTRINGLENGTH,sizeof(char)))){fprintf(stderr,"\n Error allocating memory\n");exit(1);}
        strncpy(pExtraInfo->CalibFileNm,arg,LONGSTRINGLENGTH);
        flag[12]=1;
        fprintf(stderr,"\n WARNING: Deprecated command '%s'. It may trigger bugs.\n Press <INTRO> to continue ...***\n",command);
      }
      else if(!strcmp(command,"EFFCALIBFILE")){
        /*Syntax:
        EFFCALIBFILE Calfilename     <---Name of the Efficiency calibration file
        (Use this if using Efficiency-style calibration files [Recomended])
        */
        arg=strtok( NULL,"\n\t \r");
        if(!(pExtraInfo->CalibFileNm=(char*)calloc(LONGSTRINGLENGTH,sizeof(char)))){fprintf(stderr,"\n Error allocating memory\n");exit(1);}
        strncpy(pExtraInfo->CalibFileNm,arg,LONGSTRINGLENGTH);
        flag[12]=2;
      }
      else if(!strcmp(command,"SAMPLEFILE")){
        /*Syntax:
        SAMPLEFILE SampleFilename     <---Name of the Sample definition file
        */
        arg=strtok( NULL,"\n\t \r");
        if(!(pExtraInfo->SampleFileNm=(char*)calloc(LONGSTRINGLENGTH,sizeof(char)))){fprintf(stderr,"\n Error allocating memory\n");exit(1);}
        strncpy(pExtraInfo->SampleFileNm,arg,LONGSTRINGLENGTH);
        flag[13]=1;
      }
      else if(!strcmp(command,"FILTERFILE")){
        /*Syntax:
        FILTERFILE Filterfilename     <---Name of the Filter definition file
        */
        arg=strtok( NULL,"\n\t \r");
        if(!(pExtraInfo->FilterFileNm=(char*)calloc(LONGSTRINGLENGTH,sizeof(char)))){fprintf(stderr,"\n Error allocating memory\n");exit(1);}
        strncpy(pExtraInfo->FilterFileNm,arg,LONGSTRINGLENGTH);
        flag[14]=1;
      }
      else if(!strcmp(command,"CFLAGSFILE")){
        /*Syntax:
        CFLAGSFILE CFlagsfilename     <---Name of the Calculation flags definition file
        */
        arg=strtok( NULL,"\n\t \r");
        if(!(pExtraInfo->AreasFileNm=(char*)calloc(LONGSTRINGLENGTH,sizeof(char)))){fprintf(stderr,"\n Error allocating memory\n");exit(1);}
        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)
        */
        arg=strtok( NULL,"\n\t \r");
        if(!(pExtraInfo->AreasFileNm=(char*)calloc(LONGSTRINGLENGTH,sizeof(char)))){fprintf(stderr,"\n Error allocating memory\n");exit(1);}
        strncpy(pExtraInfo->AreasFileNm,arg,LONGSTRINGLENGTH);
        pExtraInfo->AreasFormat=2;
        flag[15]=1;
      }
      else if(!strcmp(command,"DBPATH")){
        /*Syntax:
        DBPATH Path   <---Path to DataBase files (optional)
        */
        arg=strtok( NULL,"\n\t \r");
        if(!(pExtraInfo->DBpath=(char*)calloc(LONGSTRINGLENGTH,sizeof(char)))){fprintf(stderr,"\n Error allocating memory\n");exit(1);}
        strncpy(pExtraInfo->DBpath,arg,LONGSTRINGLENGTH);
      }
      else if(!strcmp(command,"OUTPUTFILE")){
        /*Syntax:
        OUTPUTFILE Outputfilename   <---Name for the output file
        */
        arg=strtok( NULL,"\n\t \r");
        if(!(pExtraInfo->OutputFileNm=(char*)calloc(LONGSTRINGLENGTH,sizeof(char)))){fprintf(stderr,"\n Error allocating memory\n");exit(1);}
        strncpy(pExtraInfo->OutputFileNm,arg,LONGSTRINGLENGTH);
        pExtraInfo->WantOutputfile=1;
        /*note: this command is optional, but a flag indicating if it was used is passed*/
      }
      else if(!strcmp(command,"DETECTORFILE")){
        /*Syntax:
        DETECTORFILE Detectorfilename   <---Name for the detector definition file (Optional)
        */
        arg=strtok( NULL,"\n\t \r");
        if(!(pExtraInfo->DetectorFileNm=(char*)calloc(LONGSTRINGLENGTH,sizeof(char)))){fprintf(stderr,"\n Error allocating memory\n");exit(1);}
        strncpy(pExtraInfo->DetectorFileNm,arg,LONGSTRINGLENGTH);
        pExtraInfo->givendetectorfile=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.;

  /*Calculate Fluence*/
  if (pexppar->ioncharge <0.){
    pexppar->ioncharge=1.;
    if (pexppar->ion.Z>1) printf("\n Warning: Ion charge not provided, assuming =1\n");
    }
  pexppar->Fluence=pexppar->simpar.ColCharge*pexppar->simpar.DTCC/(pexppar->ioncharge*ELECTRONCHARGE_IN_UC);

  /*Decide which scheme for calibration will be used*/
  if (flag[7]==1 && flag[12]==1) pExtraInfo->useefficiencies=0;
  else if (flag[7]>=2 && flag[12]==2) pExtraInfo->useefficiencies=1;
  else {fprintf(stderr,"\n ERROR:use either [DETCOLFAC+CALIBFILE] or [SOLIDANGLE+EFFCALIBFILE] \n");exit(1);}

  return(0);

}




int symbol2Z
(char *symbol){
  int Z;
/*  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,ChemicalSymbol[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);
}


int Z2mass(int Z, double *mass, char option)
{
  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]);
}


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);
}




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);
}

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);
}


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 for Z=%2d (%2s) read", pCFlags->atnum,tempchem);
      #endif
      #if CPIXEVERBOSITY > 0
      printf("\n");
      fprintCALIBYLD(stdout,&pCFlags->simareas);
      #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);
}


int readEff(const char *CalibFileNm, const int *PresentElems, int MaxZinsample, CalibYld **EffArray)
{
  FILE *CalibFile;
  char dummy[LONGSTRINGLENGTH];
  int maxZ;
  int i,tempZ,result,version;
  ChemSymb tempchem;
  double trash;

  CalibFile = fopen(CalibFileNm, "r");
  if (CalibFile == NULL)
    {fprintf(stderr,"\nERROR: Could not open file '%s' for read.\n",CalibFileNm);exit(1);}

  /*Check file version*/
  do {
    fscanf(CalibFile,"%255s",dummy);
    if(feof(CalibFile)) {fprintf(stderr,"\nERROR: 'EFFIFILEVERSION' not found in '%s'\n",CalibFileNm);exit(1);}
  } while (strcasecmp(dummy, "EFFIFILEVERSION"));

  fscanf(CalibFile, "%d", &version);
  if(version!=EFFIFILEVERSION1){fprintf(stderr,"\nERROR: Calib. file '%s' version=%d!=%d\n",CalibFileNm,version,EFFIFILEVERSION1);exit(1);}

  /*Read Max Z*/
  do {
    fscanf(CalibFile,"%255s",dummy);
    if(feof(CalibFile)) {fprintf(stderr,"\nERROR: 'MAXZ' not found in '%s'\n",CalibFileNm);exit(1);}
  } while (strcasecmp(dummy, "MAXZ"));

  fscanf(CalibFile, "%d", &maxZ);
  if(MaxZinsample>maxZ){fprintf(stderr,"\nERROR: Calib. file '%s' maxZ=%d<%d\n",CalibFileNm,maxZ,MaxZinsample);exit(1);}

  /*Initialize EffArray */
  *EffArray=(CalibYld*)calloc(MaxZinsample+1,sizeof(CalibYld));
  if (*EffArray==NULL){fprintf(stderr,"\n Error allocating memory (EffArray)\n");exit(1);}

  do {
    fscanf(CalibFile,"%255s",dummy);
    if(feof(CalibFile)) {fprintf(stderr,"\nERROR: 'DATA' not found in '%s'\n",CalibFileNm);exit(1);}
  } while (strcasecmp(dummy, "DATA"));

  /*Read Data blocks*/
  for(result=0,tempZ=0;!feof(CalibFile) && tempZ!=MaxZinsample;){
    fscanf(CalibFile, "%d", &tempZ);
    if(tempZ<1 || tempZ>MaxZinsample){
      fprintf(stderr,"\nERROR: Bad format in Calib. file '%s' (Z=%d)\n",CalibFileNm,tempZ);
      exit(1);
    }
    fscanf(CalibFile,"%3s",tempchem);
    if(tempZ != symbol2Z(tempchem)){
      fprintf(stderr,"\nERROR: Bad format in Calib. file '%s' (Z=%d)\n",CalibFileNm,tempZ);
      exit(1);
    }
    if(PresentElems[tempZ]){
      result++;

      fscanCALIBYLD(CalibFile,&(*EffArray)[tempZ]); //Detector efficiency

      #if CPIXEVERBOSITY > 0
      printf("\n Efficiency for Z=%2d (%2s) read", tempZ,tempchem);
      #endif
      #if CPIXEVERBOSITY > 1
      printf("\n");
      fprintCALIBYLD(stdout,&(*EffArray)[tempZ]);
      #endif

    }
    else for(i=0;i<15;i++)fscanf(CalibFile, "%lg",&trash);//just discard the 15 following numbers
  }
  fclose(CalibFile);

  return(result);
}



int readEff2(const char *CalibFileNm, const int *PresentElems, int MaxZinsample, const CalibYld *LinesEnerArray, CalibYld **EffArray)
{
  FILE *CalibFile;
  char dummy[LONGSTRINGLENGTH];
  int i,jj,ll,tempZ,result,version,ndata;
  double tempEner,tempEff;
  double *tempEnerArray,*tempEffArray;
  const CalibYld *pLineEner;
  CalibYld *pEff;

  CalibFile = fopen(CalibFileNm, "r");
  if (CalibFile == NULL)
    {fprintf(stderr,"\nERROR: Could not open file '%s' for read.\n",CalibFileNm);exit(1);}

  /*Check file version*/
  do {
    fscanf(CalibFile,"%255s",dummy);
    if(feof(CalibFile)) {fprintf(stderr,"\nERROR: 'EFFIFILEVERSION' not found in '%s'\n",CalibFileNm);exit(1);}
  } while (strcasecmp(dummy, "EFFIFILEVERSION"));

  fscanf(CalibFile, "%d", &version);
  if(version!=EFFIFILEVERSION2){fprintf(stderr,"\nERROR: Calib. file '%s' version=%d!=%d\n",CalibFileNm,version,EFFIFILEVERSION2);exit(1);}

  /*Initialize EffArray */
  *EffArray=(CalibYld*)calloc(MaxZinsample+1,sizeof(CalibYld));
  if (*EffArray==NULL){fprintf(stderr,"\n Error allocating memory (EffArray)\n");exit(1);}

  do {
    fscanf(CalibFile,"%255s",dummy);
    if(feof(CalibFile)) {fprintf(stderr,"\nERROR: 'DATA' not found in '%s'\n",CalibFileNm);exit(1);}
  } while (strcasecmp(dummy, "DATA"));

  /*Check the length of the data (number of data points) and leave the file
   cursor at the same point for further reading*/
  for (ndata=-1;!feof(CalibFile);ndata++){
    fscanf(CalibFile,"%lf%lf",&tempEner,&tempEff);
  }
  fclose(CalibFile);
  CalibFile = fopen(CalibFileNm, "r");
  if (CalibFile == NULL)
     {fprintf(stderr,"\nERROR: Could not open file '%s' for read.\n",CalibFileNm);exit(1);}
  do {
    fscanf(CalibFile,"%255s",dummy);
    if(feof(CalibFile)) {fprintf(stderr,"\nERROR: 'DATA' not found in '%s'\n",CalibFileNm);exit(1);}
  } while (strcasecmp(dummy, "DATA"));

  /*Initialize tempEnerArray and tempEeffArray*/
  tempEnerArray=(double*)calloc(ndata,sizeof(double));
  if (tempEnerArray==NULL){fprintf(stderr,"\n Error allocating memory (tempEnerArray)\n");exit(1);}
  tempEffArray=(double*)calloc(ndata,sizeof(double));
  if (tempEffArray==NULL){fprintf(stderr,"\n Error allocating memory (tempEffArray)\n");exit(1);}

  /*Load the efficiencies versus energy in memory*/
  for(i=0,tempEner=-1.;i<ndata;i++){
    fscanf(CalibFile,"%lf%lf",&tempEnerArray[i],&tempEffArray[i]);
    //printf("\nDEBUG: %d)  %lf\t%lf",i,tempEnerArray[i],tempEffArray[i]);
    if (tempEner>tempEnerArray[i] || tempEnerArray[i]<0 )
       {fprintf(stderr,"\nERROR: format not valid in '%s' (energies must be sorted and not negative).\n",CalibFileNm);exit(1);}
  }

  /*Fill the Effarray blocks by interpolating each line energy from the tempEffArray values*/
  for(result=0,tempZ=0; tempZ<=MaxZinsample;tempZ++){
    if(PresentElems[tempZ]){
      result++;
      pLineEner=&LinesEnerArray[tempZ];
      pEff=&(*EffArray)[tempZ];

      /*Fill K lines*/
      for (jj = 1; jj <= 3; jj++) pEff->K_[jj]=interpolate_efficiency(ndata,tempEnerArray,tempEffArray,pLineEner->K_[jj]);

      /*Fill M lines*/
      for (jj = 1; jj <= 3; jj++) pEff->M_[jj]=interpolate_efficiency(ndata,tempEnerArray,tempEffArray,pLineEner->M_[jj]);

      /*Fill L lines*/
      for (jj = 1; jj <= 3; jj++) for (ll = 1; ll <= 3; ll++)
        pEff->L_[jj][ll]=interpolate_efficiency(ndata,tempEnerArray,tempEffArray,pLineEner->L_[jj][ll]);

      #if CPIXEVERBOSITY > 0
      fprintf(stdout,"\n Efficiencies for Z=%2d read", tempZ);
      #endif
      #if CPIXEVERBOSITY > 1
      fprintf(stdout,"\n");
      fprintCALIBYLD(stdout,&(*EffArray)[tempZ]);
      #endif

    }
  }
  fclose(CalibFile);
  free(tempEnerArray);
  free(tempEffArray);

  return(result);
}

double interpolate_efficiency(int ndata, const double *Xarray, const double *Yarray, double X)
{
  int i;
  double result;

  if(X<=0.) return 0.;
  if(X>Xarray[ndata-1]) {
    #if CPIXEVERBOSITY > 0
    printf("\nWARNING:Line energy (%g keV) exceded maximum energy in efficiency database (%g keV)", X,Xarray[ndata]);
    #endif
    return 0.;
  }
  /*finds the index of the energy inmediately superior*/
  for(i=0;Xarray[i]<=X && i<ndata;i++);
  /*returns the interpolated value*/
  result=(X-Xarray[i-1])*(Yarray[i]-Yarray[i-1])/(Xarray[i]-Xarray[i-1])+Yarray[i-1];
  //printf("\nDEBUG:--> EL=%lg =>i=%d   E=%lg - %lf   eff=%lg - %lg",X,i,Xarray[i-1],Xarray[i],Yarray[i-1],Yarray[i]);
  //printf("\nDEBUG:--> Interpol:%lg ",result);
  return result;
}



int readEff3(const char *CalibFileNm, const int *PresentElems, int MaxZinsample, const CalibYld *LinesEnerArray, CalibYld **EffArray)
{
  FILE *CalibFile;
  char dummy[LONGSTRINGLENGTH],dummy2[LONGSTRINGLENGTH];
  int i,jj,ll,tempZ,result,version,ndata,singlelineflag,exceptionflag;
  double tempEner,tempEff;
  double *tempEnerArray,*tempEffArray;
  ChemSymb tempchem="";
  const CalibYld *pLineEner;
  CalibYld *pEff;
  double trash;

  CalibFile = fopen(CalibFileNm, "r");
  if (CalibFile == NULL)
    {fprintf(stderr,"\nERROR: Could not open file '%s' for read.\n",CalibFileNm);exit(1);}

  /*Check file version*/
  do {
    fscanf(CalibFile,"%255s",dummy);
    if(feof(CalibFile)) {fprintf(stderr,"\nERROR: 'EFFIFILEVERSION' not found in '%s'\n",CalibFileNm);exit(1);}
  } while (strcasecmp(dummy, "EFFIFILEVERSION"));

  fscanf(CalibFile, "%d", &version);
  if(version!=EFFIFILEVERSION3){fprintf(stderr,"\nERROR: Calib. file '%s' version=%d!=%d\n",CalibFileNm,version,EFFIFILEVERSION3);exit(1);}

  /*Initialize EffArray */
  *EffArray=(CalibYld*)calloc(MaxZinsample+1,sizeof(CalibYld));
  if (*EffArray==NULL){fprintf(stderr,"\n Error allocating memory (EffArray)\n");exit(1);}

  do {
    fscanf(CalibFile,"%255s",dummy);
    if(feof(CalibFile)) {fprintf(stderr,"\nERROR: 'DETECTORCURVE' not found in '%s'\n",CalibFileNm);exit(1);}
  } while (strcasecmp(dummy, "DETECTORCURVE"));

  /*Check the length of the data (number of data points) and leave the file
   cursor at the same point for further reading*/
  for (ndata=-1;!feof(CalibFile) && strcasecmp(dummy, "LINEEFFICIENCIES") ;ndata++){
    fscanf(CalibFile,"%255s%255s",dummy,dummy2);
    if(strcasecmp(dummy, "LINEEFFICIENCIES"))exceptionflag=1;
  }
  fclose(CalibFile);
  CalibFile = fopen(CalibFileNm, "r");
  if (CalibFile == NULL)
     {fprintf(stderr,"\nERROR: Could not open file '%s' for read.\n",CalibFileNm);exit(1);}
  do {
    fscanf(CalibFile,"%255s",dummy);
    if(feof(CalibFile)) {fprintf(stderr,"\nERROR: 'DETECTORCURVE' not found in '%s'\n",CalibFileNm);exit(1);}
  } while (strcasecmp(dummy, "DETECTORCURVE"));

  //fprintf(stderr,"\n\nDEBUG--->%d\n",ndata);getchar();

  /*Initialize tempEnerArray and tempEeffArray*/
  tempEnerArray=(double*)calloc(ndata,sizeof(double));
  if (tempEnerArray==NULL){fprintf(stderr,"\n Error allocating memory (tempEnerArray)\n");exit(1);}
  tempEffArray=(double*)calloc(ndata,sizeof(double));
  if (tempEffArray==NULL){fprintf(stderr,"\n Error allocating memory (tempEffArray)\n");exit(1);}

  /*Load the efficiencies versus energy in memory*/
  for(i=0,tempEner=-1.;i<ndata;i++){
    fscanf(CalibFile,"%lf%lf",&tempEnerArray[i],&tempEffArray[i]);
    //printf("\nDEBUG: %d)  %lf\t%lf",i,tempEnerArray[i],tempEffArray[i]);
    if (tempEner>tempEnerArray[i] || tempEnerArray[i]<0 )
       {fprintf(stderr,"\nERROR: format not valid in '%s' (energies must be sorted and not negative).\n",CalibFileNm);exit(1);}
  }

  /*Fill the Effarray blocks by interpolating each line energy from the tempEffArray values*/
  for(result=0,tempZ=0; tempZ<=MaxZinsample;tempZ++){
    if(PresentElems[tempZ]){
      result++;
      pLineEner=&LinesEnerArray[tempZ];
      pEff=&(*EffArray)[tempZ];

      /*Fill K lines*/
      for (jj = 1; jj <= 3; jj++) pEff->K_[jj]=interpolate_efficiency(ndata,tempEnerArray,tempEffArray,pLineEner->K_[jj]);

      /*Fill M lines*/
      for (jj = 1; jj <= 3; jj++) pEff->M_[jj]=interpolate_efficiency(ndata,tempEnerArray,tempEffArray,pLineEner->M_[jj]);

      /*Fill L lines*/
      for (jj = 1; jj <= 3; jj++) for (ll = 1; ll <= 3; ll++)
        pEff->L_[jj][ll]=interpolate_efficiency(ndata,tempEnerArray,tempEffArray,pLineEner->L_[jj][ll]);

      #if CPIXEVERBOSITY > 0
      fprintf(stdout,"\n Efficiencies for %2s extracted from detector curve", ChemicalSymbol[tempZ]);
      #endif
      #if CPIXEVERBOSITY > 1
      fprintf(stdout,"\n");
      fprintCALIBYLD(stdout,&(*EffArray)[tempZ]);
      #endif

    }
  }

  /*Now overwrite with exceptions if present*/
  if(exceptionflag){
    fscanf(CalibFile,"%255s",dummy); //this just reads and discards the LINEFFICIENCIES keyword
    for(result=0,tempZ=0;!feof(CalibFile); ){
      fscanf(CalibFile, "%d", &tempZ);
      fscanf(CalibFile,"%3s",tempchem);
      if(tempZ==0){ //This means that a single line will be modified (instead of a block)
        tempZ=symbol2Z(tempchem);
        singlelineflag=1;
        //printf("\nDEBUG: :::::::::: '%s'----'%d'\n",tempchem,tempZ);
      }

      if(tempZ != symbol2Z(tempchem)){
        fprintf(stderr,"\nERROR: Bad format in Calib. file '%s' (Z=%d)\n",CalibFileNm,tempZ);exit(1);
        exit(1);
      }
      if(singlelineflag){  /*Read Single line*/
        fscanf(CalibFile,"%255s%lf",dummy,&tempEff);
        if( tempZ<=MaxZinsample && PresentElems[tempZ]){  //This relies in the fact that the left check is done BEFORE the right one
          pEff=&(*EffArray)[tempZ];
          if      (!strcasecmp(dummy, "Ka12" ))pEff->K_[1]=tempEff;
          else if (!strcasecmp(dummy, "Kb1"))pEff->K_[2]=tempEff;
          else if (!strcasecmp(dummy, "Kb2"))pEff->K_[3]=tempEff;

          else if (!strcasecmp(dummy, "La12" ))pEff->L_[1][1]=tempEff;
          else if (!strcasecmp(dummy, "Lb2"))pEff->L_[1][2]=tempEff;
          else if (!strcasecmp(dummy, "Ll" ))pEff->L_[1][3]=tempEff;

          else if (!strcasecmp(dummy, "Lb1"))pEff->L_[2][1]=tempEff;
          else if (!strcasecmp(dummy, "Lg1"))pEff->L_[2][2]=tempEff;
          else if (!strcasecmp(dummy, "Le" ))pEff->L_[2][3]=tempEff;

          else if (!strcasecmp(dummy, "Lb3"))pEff->L_[3][1]=tempEff;
          else if (!strcasecmp(dummy, "Lb4"))pEff->L_[3][2]=tempEff;
          else if (!strcasecmp(dummy, "Lg3"))pEff->L_[3][3]=tempEff;

          else if (!strcasecmp(dummy, "Ma" ))pEff->M_[1]=tempEff;
          else if (!strcasecmp(dummy, "Mb" ))pEff->M_[2]=tempEff;
          else if (!strcasecmp(dummy, "Mg" ))pEff->M_[3]=tempEff;
          else {
            fprintf(stderr,"\nERROR: Bad format in Calib. file '%s'. Unrecognized line name '%s'\n",CalibFileNm,dummy);
            exit(1);
          }
          #if CPIXEVERBOSITY > 0
          printf("\n Exception Efficiency for %2s (Line %s) read", tempchem,dummy);
          #endif
          #if CPIXEVERBOSITY > 1
          printf("\n");
          fprintCALIBYLD(stdout,&(*EffArray)[tempZ]);
          #endif
        }
        else{
          #if CPIXEVERBOSITY > 0
          printf("\n Warning: Exception for %2s ignored (not present in sample).",tempchem);
          #endif
        }
      }
      else{  /*Read Data blocks*/
        if(tempZ<=MaxZinsample && PresentElems[tempZ]){
          result++;
          fscanCALIBYLD(CalibFile,&(*EffArray)[tempZ]); //Detector efficiency

          #if CPIXEVERBOSITY > 0
          printf("\n Exception Efficiency Block for Z=%2d (%2s) read", tempZ,tempchem);
          #endif
          #if CPIXEVERBOSITY > 1
          printf("\n");
          fprintCALIBYLD(stdout,&(*EffArray)[tempZ]);
          #endif
        }
        else for(i=0;i<15;i++)fscanf(CalibFile, "%lg",&trash);//just discard the 15 following numbers
      }
    }
  }
  fclose(CalibFile);
  free(tempEnerArray);
  free(tempEffArray);

  return(result);
}

int readLinesEner(const char *LinesEnerFileNm, const int *PresentElems, int MaxZinsample, CalibYld **LinesEnerArray)
{
  FILE *LinesEnerFile;
  char dummy[LONGSTRINGLENGTH];
  int maxZ;
  int i,tempZ,result;
  ChemSymb tempchem;
  double trash;

  LinesEnerFile = fopen(LinesEnerFileNm, "r");
  if (LinesEnerFile == NULL)
    {fprintf(stderr,"\nERROR: Could not open file '%s' for read.\n",LinesEnerFileNm);exit(1);}


  /*Read Max Z*/
  do {
    fscanf(LinesEnerFile,"%255s",dummy);
    if(feof(LinesEnerFile)) {fprintf(stderr,"\nERROR: 'MAXZ' not found in '%s'\n",LinesEnerFileNm);exit(1);}
  } while (strcasecmp(dummy, "MAXZ"));

  fscanf(LinesEnerFile, "%d", &maxZ);
  if(MaxZinsample>maxZ){fprintf(stderr,"\nERROR: while reading file '%s' maxZ=%d<%d\n",LinesEnerFileNm,maxZ,MaxZinsample);exit(1);}

  /*Initialize LinesEnerArray */
  *LinesEnerArray=(CalibYld*)calloc(MaxZinsample+1,sizeof(CalibYld));
  if (*LinesEnerArray==NULL){fprintf(stderr,"\n Error allocating memory (LinesEnerArray)\n");exit(1);}

  do {
    fscanf(LinesEnerFile,"%255s",dummy);
    if(feof(LinesEnerFile)) {fprintf(stderr,"\nERROR: 'DATA' not found in '%s'\n",LinesEnerFileNm);exit(1);}
  } while (strcasecmp(dummy, "DATA"));

  /*Read Data blocks*/
  for(result=0,tempZ=0;!feof(LinesEnerFile) && tempZ!=MaxZinsample;){
    fscanf(LinesEnerFile, "%d", &tempZ);
    if(tempZ<1 || tempZ>MaxZinsample){
      fprintf(stderr,"\nERROR: Bad format in LinesEner file '%s' (Z=%d)\n",LinesEnerFileNm,tempZ);
      exit(1);
    }
    fscanf(LinesEnerFile,"%3s",tempchem);
    if(tempZ != symbol2Z(tempchem)){
      fprintf(stderr,"\nERROR: Bad format in LinesEner file '%s' (Z=%d)\n",LinesEnerFileNm,tempZ);
      exit(1);
    }
    if(PresentElems[tempZ]){
      result++;

      fscanCALIBYLD(LinesEnerFile,&(*LinesEnerArray)[tempZ]);    //emission energies
      //fscanCALIBYLD(LinesEnerFile,&pEff->raweffi); //Detector efficiency

      #if CPIXEVERBOSITY > 0
      printf("\n Lines Energy for Z=%2d (%2s) read", tempZ,tempchem);
      #endif
      #if CPIXEVERBOSITY > 1
      printf("\n");
      fprintCALIBYLD(stdout,&(*LinesEnerArray)[tempZ]);
      #endif

    }
    else for(i=0;i<15;i++)fscanf(LinesEnerFile, "%lg",&trash);//just discard the 15 following numbers
  }
  fclose(LinesEnerFile);

  return(result);
}


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);
}




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);
}


int readAbsCoef(const char *TotAbsCoefFileNm, int MaxZinsample, AbsCoef **TotAbsCoefArray)
{
  FILE *TotAbsCoefFile;
  char dummy[LONGSTRINGLENGTH];
  int maxZ,Nread;
  AbsCoef *pTotAbsCoef;
  int i,tempZ;
  ChemSymb tempchem;
  double 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 to hold the WHOLE table*/
  *TotAbsCoefArray=(AbsCoef*)calloc(maxZ+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(Nread=0,tempZ=0;!feof(TotAbsCoefFile) && tempZ!=maxZ;Nread++){
    fscanf(TotAbsCoefFile, "%d", &tempZ);
    if(tempZ<1 || tempZ>maxZ || (*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);
    }
    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 > 1
    printf("\n TotAbsCoef[%2d] (%2s) read", pTotAbsCoef->atnum,tempchem);
    #endif
  }
  fclose(TotAbsCoefFile);
  #if CPIXEVERBOSITY > 0
  printf("\n Absorption Coefs. for %d elements read", Nread);
  #endif

  return(0);
}


int readAbsCoef_OLD(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);
}



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


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);
}



int initlyrarray(const EXP_PARAM *pexp, const foil *MatArray, const CalibYld *LinesEnerArray, const AbsCoef *TotAbsCoefArray, const SPT *SPTArray, const int *PresentElems, int NFoil, int *NFoilUsed, LYR **plyrarray)
{
  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 index 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, LinesEnerArray, 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);
}


int initlyr(const EXP_PARAM *pexp, const CalibYld *LinesEnerArray, 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);}

  plyr->TrcUse=(int*)calloc(plyr->NumOfTrc,sizeof(int));
  if (plyr->TrcUse==NULL){fprintf(stderr,"\n Error allocating memory (plyr->TrcUse)\n");exit(1);}

   plyr->NeedSFC=(int*)calloc(plyr->NumOfTrc,sizeof(int));
   if (plyr->NeedSFC==NULL){fprintf(stderr,"\n Error allocating memory (plyr->NeedSFC)\n");exit(1);}
//   plyr->NeedSFC=NULL;

  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(LinesEnerArray, 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*/
    Z=plyr->pFoil->comp.elem[i].Z;
    plyr->TrcUse[i]= (Z > minK);
    if (plyr->TrcUse[i]) {
      /*Initialization for plyr->ResXYldArray */
      ResYld=&(plyr->ResYldArray[i]);
      ResYld->ener=LinesEnerArray[Z];
      ResYld->atnum=Z;
      /*Some (stupid?) initialization for SecResArray*/
      plyr->SecResArray[i].Pk.atnum=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 ONLY FOR THOSE ELEMENTS THAT WILL BE NEEDED*/
      /*Note: */

      plyr->TrcAtnum[i] = Z;   //<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 && Z > minM) {
        for (jj = 1; jj <= 3; jj++) {
          AbsFac->M_[jj] = TotAbsor(TotAbsCoefArray, &plyr->pFoil->comp, ResYld->ener.M_[jj]);
          /* if *MACoef<AbsC[ii].M[jj] then *MACoef:=AbsC[ii].M[jj] ;  */
        }
      }
      if (ResYld->ener.L_[1][1] > 0 && Z > 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.;  

  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);
}

int CreateEff(const char *TotAbsCoefFileNm, const  char *LinesEnerFileNm, const char *WindowDefFileNm, const char *CrystalDefFileNm, double factor ,CalibYld **EffBlockArray, int *dimEffTable, TwoCol **EffTable)
{
  int maxZ=92, nlines=15;
  int ielem,jj,ll,*AllPresentElems;
  CalibYld *pEff, *pwin, *pcry, *pLinesEner, *LinesEner,*WindowAFTArray,*CrystalAFTArray;
  TwoCol *pEffTable;

  /*The maximum dimension of Energy and efficiency vectors*/
  *dimEffTable=maxZ*nlines;

  /*initialize The arrays that should be returned*/
  *EffTable=(TwoCol*)calloc(*dimEffTable,sizeof(double));
  if (*EffTable==NULL){fprintf(stderr,"\n Error allocating memory (EffTable)\n");exit(1);}
  *EffBlockArray=(CalibYld*)calloc(maxZ+1,sizeof(CalibYld));
  if (*EffBlockArray==NULL){fprintf(stderr,"\n Error allocating memory (EffBlockArray)\n");exit(1);}

  /*initialize the arrays that won't be returned (local use only)*/
  WindowAFTArray=(CalibYld*)calloc(maxZ+1,sizeof(CalibYld));
  if (WindowAFTArray==NULL){fprintf(stderr,"\n Error allocating memory (WindowAFTArray)\n");exit(1);}
  CrystalAFTArray=(CalibYld*)calloc(maxZ+1,sizeof(CalibYld));
  if (CrystalAFTArray==NULL){fprintf(stderr,"\n Error allocating memory (CrystalAFTArray)\n");exit(1);}
  AllPresentElems=(int*)calloc(maxZ+1,sizeof(int));
  if (AllPresentElems==NULL){fprintf(stderr,"\n Error allocating memory (AllPresentElems)\n");exit(1);}
  for(AllPresentElems[0]=0,jj=1 ; jj<=maxZ ; jj++)AllPresentElems[jj]=1; //create fake PresentElems (all present)
  readLinesEner(LinesEnerFileNm, AllPresentElems, maxZ, &LinesEner);

  AllFilterTrans(TotAbsCoefFileNm, LinesEnerFileNm, WindowDefFileNm, WindowAFTArray);
  AllFilterTrans(TotAbsCoefFileNm, LinesEnerFileNm, CrystalDefFileNm, CrystalAFTArray);

  /*Reset the dimension for Energy and Efficiency*/
  *dimEffTable=0;

  for (ielem=1,pEffTable=*EffTable;ielem<maxZ;ielem++){
    pEff=&(*EffBlockArray)[ielem];
    pwin=&WindowAFTArray[ielem];
    pcry=&CrystalAFTArray[ielem];
    pLinesEner=&LinesEner[ielem];
    /*Fill K lines*/
    for (jj = 1; jj <= 3; jj++) {
      pEff->K_[jj]=pwin->K_[jj] * (1.-pcry->K_[jj]) * factor;
      if( pLinesEner->K_[jj]>0.){
        pEffTable->x=pLinesEner->K_[jj];
        pEffTable->y=pEff->K_[jj];
        pEffTable++;
      }
    }
    /*Fill M lines*/
    for (jj = 1; jj <= 3; jj++){
      pEff->M_[jj]=pwin->M_[jj] * (1.-pcry->M_[jj]) * factor;
      if( pLinesEner->M_[jj]>0.){
        pEffTable->x=pLinesEner->M_[jj];
        pEffTable->y=pEff->M_[jj];
        pEffTable++;
      }
    }
    /*Fill L lines*/
    for (jj = 1; jj <= 3; jj++) for (ll = 1; ll <= 3; ll++) {
      pEff->L_[jj][ll]=pwin->L_[jj][ll] * (1.-pcry->L_[jj][ll]) * factor;
      if( pLinesEner->L_[jj][ll]>0.){
        pEffTable->x=pLinesEner->L_[jj][ll];
        pEffTable->y=pEff->L_[jj][ll];
        pEffTable++;
      }
    }

    #if CPIXEVERBOSITY > 0
    fprintf(stdout,"\n Efficiency for %2s:", ChemicalSymbol[ielem]);
    #endif
    #if CPIXEVERBOSITY > 1
    fprintf(stdout,"\n");
    fprintCALIBYLD(stdout,&(*EffBlockArray)[ielem]);
    #endif

  }

  /*Reallocate the Energy and Efficiency vectors to free unused space*/
  *dimEffTable=pEffTable-(*EffTable);
  *EffTable=(TwoCol*)realloc(*EffTable,(*dimEffTable)*sizeof(TwoCol));
  if (*EffTable==NULL){fprintf(stderr,"\n Error reallocating memory (Efftable)\n");exit(1);}

  /*Sort the arrays by ascending energy*/
  qsort(*EffTable,*dimEffTable,sizeof(TwoCol),(void*)compare_Twocol_x);

  #if CPIXEVERBOSITY > 0
  fprintf(stdout,"\n Efficiency for %d lines generated", *dimEffTable);
  #endif
  #if CPIXEVERBOSITY > 1
  for(jj=0; jj<*dimEffTable; jj++)fprintf(stdout,"\n%le\t\%le",(*EffTable)[jj].x,(*EffTable)[jj].y);
  #endif

  free(WindowAFTArray);
  free(CrystalAFTArray);
  free(LinesEner);
  free(AllPresentElems);
  return(maxZ);
}

int AllFilterTrans(const char *TotAbsCoefFileNm, const  char *LinesEnerFileNm, const char *FilterDefFileNm, CalibYld *AFTArray)
{
  int i;
  int maxZ=92;
  int *AllPresentElems;
  CalibYld *LinesEnerArray;
  AbsCoef *TotAbsCoefArray;
  FILTER Filter;

  AllPresentElems=(int*)calloc(maxZ+1,sizeof(int));
  if (AllPresentElems==NULL){fprintf(stderr,"\n Error allocating memory (AllPresentElems)\n");exit(1);}
  for(AllPresentElems[0]=0,i=1 ; i<=maxZ ; i++)AllPresentElems[i]=1; //create fake PresentElems

  /*Read the energies for All the lines*/
  readLinesEner(LinesEnerFileNm, AllPresentElems, maxZ, &LinesEnerArray);
  /*Read the Absorption Coefs*/
  readAbsCoef(TotAbsCoefFileNm, maxZ, &TotAbsCoefArray);
  /*Read the Filter*/
  readFilter(FilterDefFileNm, &Filter);

  /*Calculate transmission for all the elements*/
  FilterTrans(maxZ, LinesEnerArray, TotAbsCoefArray, AllPresentElems, &Filter);
  for(i=1;i<maxZ;i++)AFTArray[i]=Filter.Trans[i];

  free(AllPresentElems);
  free(TotAbsCoefArray);
  free(LinesEnerArray);
  freeFilter(&Filter);
  return(maxZ);
}


int FilterTrans(int MaxZinsample, const CalibYld *LinesEnerArray, 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 }};

  //Initialize the Filter->Trans array
  Filter->dimTrans=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(LinesEnerArray, 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);
}


int Transmission(const CalibYld *LinesEnerArray, 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=&LinesEnerArray[Z];  //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);
}

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;
  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);
}

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

  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);
}



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;

  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);
}


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)*/
   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)++;


//      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);
}



int integrate_Simpson2(const EXP_PARAM *pexp, const AbsCoef *TotAbsCoefArray,
                       const FluorCKCoef *FCKCoefArray, const CalibYld *RawEffiArray,
                       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 *pEff;
  SecResType *pSecRes;
  CalibYld tempPenInt;
  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);
        pEff=&RawEffiArray[tempZ];
        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);

        tempPenInt=pResYld->XYld; //make a copy of the penetration integral results for this element
        deNormalize2(pexp, tempZ, attfraction, pEff, &tempPenInt, &pResYld->XYld, pXYldSum);
      }
    }
  }


  return(0);
}



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);
}

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: 
  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) { 
      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) {
      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);
}

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)));
}



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));
}

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);
}


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;            
  cau = 137.03604;         
  a02 = 2.800283608e-21;   

  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;   
  a2s = (double)(Ncam * Ncam) / Z2b;
  TEMP = 1 + EnrSm1;
  TEMP = 1 / (TEMP * TEMP);
  v1 = cau * sqrt(1 - TEMP);   
  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);
      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;
      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);
}

double ReisX_gs(FwTipo T, int SS, double ksis)
{
  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);
}


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));
}


double ReisX_En(double z, int niu)
{
  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+= 1+niu/aux2+niu*(niu-2*z)/aux4; //patch contributed by AT (was: Eaux++)
  Eaux = exp(-z) * Eaux / aux;
  return(Eaux);
}

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);
}

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);
  }
}



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


}


double Simps(double a, double b, double fa, double fi, double fb)
{
  return ((b - a) / 6. * (fa + 4. * fi + fb));
}



void deNormalize(const EXP_PARAM *pexp, const AbsCoef *AbsC, const FluorCKCoef *pFCK, atomicnumber Z2, double M2, double attfraction, const CalibYld *pXYld, XrayYield *ResY, CalibYld *XYldSum)
{
  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){
    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); //Calculate cross section at the calibration energy

  switch (ResY->atnum) {

  case 10:
  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) {
      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) {
      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;
  }
}




void deNormalize2(const EXP_PARAM *pexp, atomicnumber Z2, double attfraction, const CalibYld *pEff, const CalibYld *PenInteg, CalibYld *XYld, CalibYld *XYldSum)
{
  int jj, ll;
  const SIM_PARAM *psim;

  psim=&pexp->simpar;
  /*K lines*/
  if (Z2 >= minK && Z2 <maxK) {
    for (jj = 1; jj <= 3; jj++) {
      XYld->K_[jj]=pexp->SAngFract * pexp->Fluence * attfraction * pEff->K_[jj] * PenInteg->K_[jj];
      XYldSum->K_[jj] +=XYld->K_[jj];
    }
  }
  /*L lines*/
  if (Z2 >= minL ) {
    for (jj = 1; jj <= 3; jj++) {
      for (ll = 1; ll <= 3; ll++) {
        XYld->L_[jj][ll]=pexp->SAngFract * pexp->Fluence * attfraction *
                          pEff->L_[jj][ll] * PenInteg->L_[jj][ll];
        XYldSum->L_[jj][ll] += XYld->L_[jj][ll];
      }
    }
  }
  /*M lines*/
  if (Z2 >= minM) { 
    for (jj = 1; jj <= 3; jj++) {
      XYld->M_[jj]=pexp->SAngFract * pexp->Fluence * attfraction * pEff->M_[jj] * PenInteg->M_[jj];
      XYldSum->M_[jj] +=XYld->M_[jj];
    }
  }
}








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

}

void freeExtraInfo(EXTRAINFO *extrainfo)
{
  safefree((void*)(&extrainfo->SampleFileNm));
  safefree((void*)(&extrainfo->FilterFileNm));
  safefree((void*)(&extrainfo->DetectorFileNm));
  safefree((void*)(&extrainfo->CalibFileNm));
  safefree((void*)(&extrainfo->AreasFileNm));
  safefree((void*)(&extrainfo->DBpath));
  safefree((void*)(&extrainfo->OutputFileNm));
}

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);
}


void safefree(void **ptr){
  if(*ptr!=NULL){
    //printf("\nDEBUG: freing %x ...",*ptr);fflush(stdout);
    free(*ptr);
    //printf(" done (%x)",*ptr);fflush(stdout);
    *ptr=NULL;
  }
}


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");
}

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

}

int compare_Twocol_x (const TwoCol *a, const TwoCol *b)
{
  double temp=a->x - b->x;
  if (temp > 0)
    return 1;
  else if (temp < 0)
    return -1;
  else
    return 0;
}

---