cpixe_fif.c

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2005 by Carlos Pascual-Izarra                           *
 *   carlos.pascual@itn.pt                                                 *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/


 
/**\file cpixe_fif.c
This is a fortran interface, not properly part of the LibCPIXE lybrary.

I wrote it to facilitate the integration with the NDF program by N.P. Barradas but it might be useful as well in the case that other people wants to link libCPIXE to a Fortran program.

I try to maintain cpixe_fif.c with at least the same funcionalities available in the cpixe.c code.


*/


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

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

#include "libcpixe.h"
#include "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];  ///IMPORTANT!!!!!!: Each isotope gets the FULL number of counts for the element!!!
  }
  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);
      //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(abscoeffilename, pSimPar->MaxZinsample, PresentElems, &Filter, &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;
}

---