cpixe_fif/cpixe_fif.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/>.

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


#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include "../libcpixe/libcpixe.h"
#include "../libcpixe/compilopt.h"





//Function prototypes:
int translatesample(int NFoil, int Nels, const int *pFOILZ,
                    const float *pFOILM, const float *pAMNT, const float *pFOILTHICK,
          int *MaxZinsample, foil **MatArray);

int translatesample2(int NFoil, int MaxNels, const int *pNels, const int *pFOILZ,
                    const float *pFOILM, const float *pAMNT, const float *pFOILTHICK,
          int *MaxZinsample, foil **MatArray);

int translateresults(int dimSums, int dimResults, const int *PresentElems,
                     const CalibYld *XYldSums, CPIXERESULTS *pRESULTS);

int translateresultsNDF(int Nels, const int *pFOILZ, const CalibYld *XYldSums, CPIXERESULTS *pRESULTS);

//Functions:

/*Translation of sample structures. */
int translatesample(int NFoil, int Nels, const int *pFOILZ,
                    const float *pFOILM, const float *pAMNT, const float *pFOILTHICK,
          int *MaxZinsample, foil **MatArray)
{
  int i,j,nelem,iels;
  double sumM;
  foil *pfoil;
  COMPOUND *c;

  /*Initialize  Sample*/
  *MatArray=(foil*)malloc((NFoil)*sizeof(foil));
  if (*MatArray==NULL){printf("\n Error allocating memory (MatArray)\n");exit(1);}
  for (*MaxZinsample=0,i = 0; i < NFoil ; i++) {

    //obtain the number of elements for this layer
    for(nelem=0, iels=0; iels<Nels ; iels++) if(pAMNT[iels+Nels*i]>0.) nelem++;  //CHECK THIS!!!
    if(nelem<1){fprintf(stderr,"\n ERROR: At least 1 element in each foil needed\n");exit(1);}

    pfoil=&((*MatArray)[i]);
    pfoil->nfoilelm=nelem;
    pfoil->thick=(double)pFOILTHICK[i];
    c=&pfoil->comp;

    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(j=0, c->sumX=0., sumM=0., iels=0 ; iels<Nels ; iels++){
      if(pAMNT[iels+Nels*i]>0.){
        c->X[j]= (double) (pAMNT[iels+Nels*i]);  //NOTE!!: The indexes are switched and "offseted" to deal with the the Fortran convention!!
      c->elem[j].Z=pFOILZ[iels];
        c->elem[j].A=Z2mass(c->elem[j].Z, &c->elem[j].IM, 'm');
        c->elem[j].M=(double)pFOILM[iels];  //The mass passed by NDF is put on the "natural mass"
        c->sumX+=c->X[j];
        sumM+=c->elem[j].M*c->X[j];
        if(c->elem[j].Z>*MaxZinsample) *MaxZinsample=c->elem[j].Z;
        j++;
      }
    }
    /*Normalize concentrations and calculate the MASS concentration*/
    for(j=0;j<nelem;j++){
      c->xn[j]=c->X[j]/c->sumX;
      c->w[j]=c->X[j]*c->elem[j].M/sumM;
    }
  }
//  for(i=0;i<NFoil;i++){
//    printf("\nFoil %d: th=%le nelem=%d",i,(*MatArray)[i].thick,(*MatArray)[i].nfoilelm);
//    for(j=0;j<(*MatArray)[i].nfoilelm;j++){
//      printf("\nTrans: (i=%d , j=%d) (*MatArray)[i].comp.elem[j].Z=%d ",i,j,(*MatArray)[i].comp.elem[j].Z);
//    printf("\nTrans:             (*MatArray)[i].comp.elem[j].M=%lf ",(*MatArray)[i].comp.elem[j].M);
//    printf("\nTrans:             (*MatArray)[i].comp.elem[j].IM=%lf ",(*MatArray)[i].comp.elem[j].IM);
//    printf("\nTrans:             (*MatArray)[i].comp.elem[j].A=%d ",(*MatArray)[i].comp.elem[j].A);
//    printf("\nTrans:      X=%lf   xn=%lf   w=%lf",(*MatArray)[i].comp.X[j],(*MatArray)[i].comp.xn[j],(*MatArray)[i].comp.w[j]);
//    }
//  }
  return(0);

}

/*Alternative translation of sample structures. */
int translatesample2(int NFoil, int MaxNels, const int *pNels, const int *pFOILZ,
                    const float *pFOILM, const float *pAMNT, const float *pFOILTHICK,
          int *MaxZinsample, foil **MatArray)
{
  int i,j,nelem;
  double sumM;
  foil *pfoil;
  COMPOUND *c;

  #if CPIXEVERBOSITY >1
  printf("\n\nDEBUG:Translatesample2(): NFoil=%d    MaxNels=%d",NFoil,MaxNels);
  for(i=0;i<NFoil;i++){
    printf("\nDEBUG:Translatesample2(): Foil %d  Nels=%d",i,pNels[i]);
    for(j=0;j<pNels[i];j++)printf("\nDEBUG:Translatesample2():Elem %d --> Z=%d",j,pFOILZ[j+MaxNels*i]);
  }
  #endif
  /*Initialize  Sample*/
  *MatArray=(foil*)malloc((NFoil)*sizeof(foil));
  if (*MatArray==NULL){printf("\n Error allocating memory (MatArray)\n");exit(1);}
  for (*MaxZinsample=0,i = 0; i < NFoil ; i++) {

    //obtain the number of elements for this layer
    nelem=pNels[i];
    if(nelem<1){fprintf(stderr,"\n ERROR: At least 1 element in each foil needed\n");exit(1);}

    pfoil=&((*MatArray)[i]);
    pfoil->nfoilelm=nelem;
    pfoil->thick=(double)pFOILTHICK[i];
    c=&pfoil->comp;

    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(j=0, c->sumX=0., sumM=0. ; j<nelem ; j++){
        c->X[j]= (double) (pAMNT[j+MaxNels*i]);  //NOTE!!: The indexes are switched and "offseted" to deal with the the Fortran convention!!
      c->elem[j].Z=pFOILZ[j+MaxNels*i];
        c->elem[j].A=Z2mass(c->elem[j].Z, &c->elem[j].IM, 'm');
        c->elem[j].M=(double)pFOILM[j+MaxNels*i];  //The mass passed by NDF is put on the "natural mass"
        c->sumX+=c->X[j];
        sumM+=c->elem[j].M*c->X[j];
        if(c->elem[j].Z>*MaxZinsample) *MaxZinsample=c->elem[j].Z;
    }
    /*Normalize concentrations and calculate the MASS concentration*/
    for(j=0;j<nelem;j++){
      c->xn[j]=c->X[j]/c->sumX;
      c->w[j]=c->X[j]*c->elem[j].M/sumM;
    }
  }
//  for(i=0;i<NFoil;i++){
//    printf("\nFoil %d: th=%le nelem=%d",i,(*MatArray)[i].thick,(*MatArray)[i].nfoilelm);
//    for(j=0;j<(*MatArray)[i].nfoilelm;j++){
//      printf("\nTrans: (i=%d , j=%d) (*MatArray)[i].comp.elem[j].Z=%d ",i,j,(*MatArray)[i].comp.elem[j].Z);
//    printf("\nTrans:             (*MatArray)[i].comp.elem[j].M=%lf ",(*MatArray)[i].comp.elem[j].M);
//    printf("\nTrans:             (*MatArray)[i].comp.elem[j].IM=%lf ",(*MatArray)[i].comp.elem[j].IM);
//    printf("\nTrans:             (*MatArray)[i].comp.elem[j].A=%d ",(*MatArray)[i].comp.elem[j].A);
//    printf("\nTrans:      X=%lf   xn=%lf   w=%lf",(*MatArray)[i].comp.X[j],(*MatArray)[i].comp.xn[j],(*MatArray)[i].comp.w[j]);
//    }
//  }
  return(0);


}


/*Copies the results of XYldSums to the array pRESULTS */
int translateresults(int dimSums, int dimResults, const int *PresentElems, const CalibYld *XYldSums, CPIXERESULTS *pRESULTS)
{
  int Z,ielem;
  for(ielem=0, Z=1 ; Z < dimSums ; Z++){
    if(PresentElems[Z]){
      ielem++;
      if(ielem>dimResults){fprintf(stderr,"\n Error: cpixefif: 'dimResults' too small\n");exit(1);}
      pRESULTS[ielem].simareas=XYldSums[Z];
      pRESULTS[ielem].atnum=Z;
    }
  }
  return(0);
}


int translateresultsNDF(int Nels, const int *pFOILZ, const CalibYld *XYldSums, CPIXERESULTS *pRESULTS)
{
  int i,Z;
  for(i=0;i<Nels;i++){
    Z=pFOILZ[i];
    pRESULTS[i].atnum=Z;
    pRESULTS[i].simareas=XYldSums[Z];  
  }
  return(0);
}

void TESTFIF()
{
  printf("\n testfif() is working...\n");
}

void CPIXEMAIN(int NCALL)
{

  STATFILENAME calibfilename="";
  STATFILENAME fckfilename="";
  STATFILENAME abscoeffilename="";
  STATFILENAME dbpath="./";


  EXP_PARAM *pExpPar;
  SIM_PARAM *pSimPar;

  foil *MatArray;
  int NFoil, NFoilUsed, dimSums;
  LYR *lyr;
  CalibYld *XYldSums;
  int i,j,l;

  /*Static variables (Those that should be maintained between calls of CPIXEMAIN()*/
  static int *PresentElems;
  static XrayYield *XYldArray;
  static FluorCKCoef *FCKCoefArray;
  static AbsCoef *TotAbsCoefArray;
  static SPT *SPTArray;
  static int lastNcall=-1;
  static FILTER Filter;

  /*Extern Variables (defined in the CPIXEIF Fortran module)*/
  extern EXP_PARAM CPIXEIF_mp_EXPPAR;
  extern char CPIXEIF_mp_CALFILENM[256];
  extern char CPIXEIF_mp_FLTFILENM[256];
  extern char CPIXEIF_mp_FCKFILENM[256];
  extern char CPIXEIF_mp_ABSFILENM[256];
  extern char CPIXEIF_mp_DBPATH[256];
  extern int CPIXEIF_mp_CPNLAY;
  extern int CPIXEIF_mp_CPNELS;
  extern int CPIXEIF_mp_DIMRESULTS;
  extern int CPIXEIF_mp_NLAYFILTER;
  extern int CPIXEIF_mp_MAXNELSFILTER;
  extern int CPIXEIF_mp_FILTERCHANGES;

  extern int CPIXEIF_mp_POINT_FOILZ, CPIXEIF_mp_POINT_FILTERZ,
             CPIXEIF_mp_POINT_FILTERNELS;     //these are DEC pointers to ints
  extern int CPIXEIF_mp_POINT_FOILTHICK , CPIXEIF_mp_POINT_FOILM,
             CPIXEIF_mp_POINT_FOILAMNT , CPIXEIF_mp_POINT_FILTERM,
             CPIXEIF_mp_POINT_FILTERTHICK,CPIXEIF_mp_POINT_FILTERAMNT ;  //these are DEC pointers to floats
  extern int CPIXEIF_mp_POINT_PIXERESULTS; //This is a DEC pointer to a CPIXERESULT
  extern int CPIXEIF_mp_POINT_PIXECALC;

  int *pFOILZ, *pFILTERZ, *pFILTERNELS; //Those are the "true" C pointers to ints
  float *pFOILTHICK ,*pFOILM ,*pAMNT,*pFILTERTHICK ,*pFILTERM ,*pFILTERAMNT;
  CPIXERESULTS *pRESULTS;
  CPIXERESULTS *pCALCFLAGS;

  #if CPIXEVERBOSITY > 0
  printf("\n\n\n************ BEGIN OF CPIXE OUTPUT ****************\n\n");
  #endif

  //Transforming fortran (DEC) pointers (which are integers) into true C pointers
  pFOILZ=(int*)CPIXEIF_mp_POINT_FOILZ;
  pFOILTHICK=(float*)CPIXEIF_mp_POINT_FOILTHICK;
  pFOILM=(float*)CPIXEIF_mp_POINT_FOILM;
  pAMNT=(float*)CPIXEIF_mp_POINT_FOILAMNT;
  pFILTERNELS=(int*)CPIXEIF_mp_POINT_FILTERNELS;
  pFILTERZ=(int*)CPIXEIF_mp_POINT_FILTERZ;
  pFILTERTHICK=(float*)CPIXEIF_mp_POINT_FILTERTHICK;
  pFILTERM=(float*)CPIXEIF_mp_POINT_FILTERM;
  pFILTERAMNT=(float*)CPIXEIF_mp_POINT_FILTERAMNT;
  pRESULTS=(CPIXERESULTS*) CPIXEIF_mp_POINT_PIXERESULTS;
  pCALCFLAGS=(CPIXERESULTS*)CPIXEIF_mp_POINT_PIXECALC;


// printf("\nFOILTHICK->%d    NELEM-->%d\n",pFOILTHICK,pNELEM);
//  printf("\nNCALL=%d\t Z=%d\n",NCALL,CPIXEIF_mp_EXPPAR.ion.Z);
//  for(i=0;i<=10;i++)for(j=0;j<3;j++)
//    printf("\nF.coef[%d][%d]=%lf",i,j,CPIXEIF_mp_SIMPAR.Filter.coef[i][j]);
//  for(j=0;j<3;j++)printf("\nF.h[%d]=%lf\n",j,CPIXEIF_mp_SIMPAR.Filter.h[j]);
//
//  printf("\nCPIXE: '%s'",CPIXEIF_mp_CALFILENM);
//  printf("\nCPIXE: 'MAXELEM=%d'",CPIXEIF_mp_NELS);
//
//  for(i=0 ; i<=3; i++) printf("\nCPIXE: '%s'",CPIXEIF_mp_CALFILENM);
//
//      /*K lines:*/
//      for (j = 1; j <= 3; j++)  printf("\n%10s %1.1lf ", LineNm[1][j], pRESULTS->K_[j] );
//      /*L-Lines*/
//      for (j = 1; j <= 3; j++)
//        for (l = 1; l <= 3; l++)
//          printf("\n%10s %1.1lf ", LineNm[j+1][l], pRESULTS->L_[j][l] );
//      /*M lines:*/
//      for (j = 1; j <= 3; j++)  printf("\n%10s %1.1lf ", LineNm[5][j], pRESULTS->M_[j] );
// printf("\n!!!!!!!!!!!\n");
// getchar();



  /*Note: The following initialization needs to be done each time the function is called*/
  MatArray=NULL;
  lyr=NULL;

  //make a convenience pointer to the EXPPAR and SIMPAR structures defined in the fortran module
  pExpPar=&CPIXEIF_mp_EXPPAR;
  pSimPar=&pExpPar->simpar;

  if(NCALL>0){
    if(lastNcall<0){fprintf(stderr,"\n ERROR: CPIXEMAIN: CPIXEMAIN was not initialized.\n");exit(1);}
    /******* BEGIN CALCULATION*********/
    /*
    These calls are all what is needed to do for each iteration
    (But don't forget also the call to freeReusable() before starting a new iteration)
    */

    //translate the sample definition made in NDF to that used in CPIXE
    NFoil=CPIXEIF_mp_CPNLAY;
    translatesample(NFoil,CPIXEIF_mp_CPNELS,pFOILZ,pFOILM,pAMNT,pFOILTHICK,&pSimPar->MaxZinsample,&MatArray);

  //translate Filter definition
    if(Filter.changes){
      Filter.nlyr=CPIXEIF_mp_NLAYFILTER;
      Filter.geomcorr=1; //filter must be perpendicular
      translatesample2(Filter.nlyr, CPIXEIF_mp_MAXNELSFILTER , pFILTERNELS, pFILTERZ,
                        pFILTERM, pFILTERAMNT, pFILTERTHICK,
                        &Filter.MaxZ, &Filter.foil );
      createPresentElems(Filter.MaxZ, Filter.nlyr, Filter.foil, &Filter.FilterElems);
//       FilterTrans(pExpPar, XYldArray, TotAbsCoefArray, PresentElems, &Filter);
//      FilterTrans(pExpPar->simpar.MaxZinsample, LinesEnerArray, TotAbsCoefArray, PresentElems, &Filter);
      //the Filter.changes variable is set externally
      Filter.changes=CPIXEIF_mp_FILTERCHANGES;
    }

//     initlyrarray(pExpPar, MatArray, XYldArray, TotAbsCoefArray, SPTArray, PresentElems, NFoil, &NFoilUsed, &lyr);
    dimSums=pSimPar->MaxZinsample+1;
    XYldSums=(CalibYld*)calloc(dimSums,sizeof(CalibYld));  //Note that XYldSums is 0-initialized
    integrate_Simpson(pExpPar, TotAbsCoefArray, FCKCoefArray, XYldArray ,NFoilUsed, &Filter, lyr, XYldSums);
    translateresultsNDF(CPIXEIF_mp_CPNELS, pFOILZ, XYldSums, pRESULTS);
      /******END CALCULATION**********/


    #if CPIXEVERBOSITY > 0
    /**********BEGIN DEBUG LINES*******/
    /* Just write stuff for debugging*/
    printf("\n\nDEBUG:  Ebeam=%lf IncAng=%lf  DetAng=%lf \n",
          pExpPar->BeamEner, pExpPar->IncAng, pExpPar->DetAng);
    printf("\nDEBUG:  Charge=%lf   DTCC=%lf DetColFac=%lf\n",
          pExpPar->simpar.ColCharge, pExpPar->simpar.DTCC, pExpPar->DetColFac);

    printf("\n\nDEBUG: Iteration %d\n",NCALL);

    for(i=0;i<NFoil;i++){
      printf("\nMatArray[%d]  (%d elem)\n",i,MatArray[i].nfoilelm);
      printf("\n---Eff. thickness=%.3le  (%d sublyrs)\n",lyr[i+1].ThickIn,lyr[i+1].FESxlen);

      for(j=0;j<MatArray[i].comp.nelem;j++){
      printf("'%s'(Z=%d) X[%d]=%lf (%lfat%% / %lfM%%)\n",
          XYldArray[MatArray[i].comp.elem[j].Z].symb,MatArray[i].comp.elem[j].Z,j,MatArray[i].comp.X[j],100*MatArray[i].comp.xn[j],100*MatArray[i].comp.w[j]);
      }
      printf("\nFinal Energy in this layer=%lf keV\n\n",lyr[i+1].FoilOutEner);
    }

    for(i=0;i<Filter.nlyr;i++){
      printf("\nFilter[%d] (%lf 1e15at/cm2)  (%d elem)  \n",i,Filter.foil[i].thick,Filter.foil[i].nfoilelm);
      for(j=0;j<Filter.foil[i].comp.nelem;j++){
      printf("(Z=%d) X[%d]=%lf (%lfat%% / %lfM%%)\n",
          Filter.foil[i].comp.elem[j].Z,j,Filter.foil[i].comp.X[j],100*Filter.foil[i].comp.xn[j],100*Filter.foil[i].comp.w[j]);
      }
    }

    for(i=0;i<dimSums;i++){
      if(i<=pSimPar->MaxZinsample && PresentElems[i]){
        printf("\n\n ************RESULTS for %s (Z=%d)\n",XYldArray[i].symb,i);
        /*K lines:*/
        for (j = 1; j <= 3; j++)
          if(XYldSums[i].K_[j] > 0.)
            printf("\n%10s (%2.2lfkeV)\t%le counts   (cal=%.3le Flt=%.2le)", LineNm[1][j], XYldArray[i].ener.K_[j], XYldSums[i].K_[j],XYldArray[i].XYld.K_[j],Filter.Trans[i].K_[j]);
        /*L-Lines*/
        for (j = 1; j <= 3; j++)
          for (l = 1; l <= 3; l++)
            if(XYldSums[i].L_[j][l] > 0.)
              printf("\n%10s (%2.2lfkeV)\t%le counts   (cal=%.3le Flt=%.2le)", LineNm[j+1][l], XYldArray[i].ener.L_[j][l], XYldSums[i].L_[j][l],XYldArray[i].XYld.L_[j][l],Filter.Trans[i].L_[j][l]);
        /*M lines:*/
        for (j = 1; j <= 3; j++)
          if(XYldSums[i].M_[j] > 0.)
            printf("\n%10s (%2.2lfkeV)\t%le counts   (cal=%.3le Flt=%.2le)",LineNm[5][j], XYldArray[i].ener.M_[j], XYldSums[i].M_[j],XYldArray[i].XYld.M_[j],Filter.Trans[i].M_[j]);
      }
    }
//    getchar();
    /**********END DEBUG LINES*******/
    #endif

    /*It is VERY important to call freeReusable() before starting a new iteration*/
    freeReusable(NFoil,&lyr,&MatArray,&XYldSums);
//    printf("\n!!!!!!! %d  %d",NCALL,Filter.changes);
    if(Filter.changes)freeFilter(&Filter);

  }

  else if(NCALL==0){
    /**********BEGIN INITIALIZATION*********/
    /*
    This Block is run only when initialization is required.
    It reads the databases (Fluorescence, Stoppings,...) and performs some operations on data
    other data that never change within a fit.
    */


    strcpy(dbpath,CPIXEIF_mp_DBPATH);
    strcpy(calibfilename,CPIXEIF_mp_CALFILENM);
    strcpy(fckfilename  ,CPIXEIF_mp_FCKFILENM);
    strcpy(abscoeffilename,CPIXEIF_mp_ABSFILENM);

    #if CPIXEVERBOSITY > 0
    printf("\nDEBUG: Initializing. Reading databases...\n\n");
    printf("\nDEBUG: Stopping Coefs DB: '%sSCOEF.95[A/B]'",dbpath);
    printf("\nDEBUG: Absorption Coefs DB: '%s'",abscoeffilename);
    printf("\nDEBUG: Fluorescence Coefs DB: '%s'",fckfilename);
    printf("\nDEBUG: Detector Calibration DB: '%s'\n",calibfilename);
    #endif

    //translate the sample definition made in NDF to that used in CPIXE
    NFoil=CPIXEIF_mp_CPNLAY;
    translatesample(NFoil,CPIXEIF_mp_CPNELS,pFOILZ,pFOILM,pAMNT,pFOILTHICK,&pSimPar->MaxZinsample,&MatArray);

    //translate Filter definition
    Filter.nlyr=CPIXEIF_mp_NLAYFILTER;
    translatesample2(Filter.nlyr, CPIXEIF_mp_MAXNELSFILTER , pFILTERNELS, pFILTERZ,
                      pFILTERM, pFILTERAMNT, pFILTERTHICK,
                      &Filter.MaxZ, &Filter.foil );
    createPresentElems(Filter.MaxZ, Filter.nlyr, Filter.foil, &Filter.FilterElems);
    Filter.changes=1;    //First time at least Filter should be read:

    createPresentElems(pSimPar->MaxZinsample, NFoil, MatArray, &PresentElems);
    readXYld(calibfilename, PresentElems, pCALCFLAGS, pSimPar->MaxZinsample, &pSimPar->CalEner, &XYldArray);
    readFCK(fckfilename, pSimPar->MaxZinsample,&FCKCoefArray);
    //readAbsCoef_OLD(abscoeffilename, pSimPar->MaxZinsample, PresentElems, &Filter, &TotAbsCoefArray);
    readAbsCoef(abscoeffilename, pSimPar->MaxZinsample, &TotAbsCoefArray);
    createSPTs(dbpath,pExpPar, pSimPar->MaxZinsample, PresentElems, (double)(2*pExpPar->ion.A), &SPTArray);

    safefree((void*)&MatArray); //Clear the initialized Matrix Array.
    safefree((void*)&Filter.foil); //Clear the initialized Filter Foil.
    /********** END INITIALIZATION *********/
  }

  else if(NCALL<0){
    /******** BEGIN TIDING UP BEFORE FINAL EXIT*********/
    /*Free everything for a clean exit*/
    #if CPIXEVERBOSITY > 0
    printf("\nDEBUG: Cleaning CPIXE Memory usage...\n\n");
    #endif

    //Note: Filter.changes=0 implies that freeFilter was not called in last iteration:
    if(!Filter.changes)freeFilter(&Filter);

    for(i=1; i<=pSimPar->MaxZinsample ;i++){
      if(PresentElems[i]){
        free(SPTArray[i].S);
        free(SPTArray[i].E);
        free(SPTArray[i].dSdE);
      }
    }
    free(SPTArray);
    free(PresentElems);
    free(XYldArray);
    free(FCKCoefArray);
    free(TotAbsCoefArray);
    free(XYldSums);
    /******** END TIDING UP BEFORE FINAL EXIT*********/
  }

//  printf("\n\nDEBUG: End of program reached\n\n");getchar();

  /*Store Number of Ncall*/
  lastNcall=NCALL;

  #if CPIXEVERBOSITY > 0
  printf("\n\n\n************ END OF CPIXE OUTPUT ****************\n\n");
  #endif

  return;
}

---