cpixe.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.c
CPIXE is a small program for demonstrating (and testing) the capabilities of the LibCPIXE library.

See the libcpixe.c file for more information.

*/


#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"


int main(int argc, char *argv[]);
int cpixemain(int NCALL);

STATFILENAME inputfilename="input.in";

/**Main program:
 * The executable accepts 1 optional command line argument indicating the input
 * file name (which defaults to "input.in").
 *
 * @return EXIT_SUCCESS if normal finishing
 */
int main(int argc, char *argv[])
{
  extern STATFILENAME inputfilename;

  if(argc>1)strcpy(inputfilename,argv[1]);
  cpixemain(0);
  cpixemain(1);
  cpixemain(-1);
  
  return EXIT_SUCCESS;
}

/** Simulates PIXE yields. See NCALL for mode of opertaion.
 * 
 * @param NCALL (I) Flag controlling the mode of operation of cpixemain: \n
      If NCALL = 0, initialization is done \n
      If NCALL > 0, the calculations are done \n
      If NCALL < 0, cleaning of memory is performed 
 *       
 *
 *
 * @return EXIT_SUCCESS if no error ocurred and "-1" if some error happened. 
 */
int cpixemain(int NCALL)
{
  STATFILENAME dbpath="./";
  extern STATFILENAME inputfilename;
  STATFILENAME samplefilenm=""; 
  STATFILENAME filterfilenm=""; 
  STATFILENAME calibfilenm="";
  STATFILENAME areasfilenm="";
  STATFILENAME fckfilename="";
  STATFILENAME abscoeffilename=""; 
  STATFILENAME outputfilenm="";
  
  
  LYR *lyr;
  CalibYld *XYldSums;
  int i,j,l;
  EXP_PARAM *pExpPar;
  SIM_PARAM *pSimPar;
 
   /*Static variables (Those that should be maintained between calls of CPIXEMAIN()*/
  static int *PresentElems;
  static CPIXERESULTS *CalcFlags;
  static XrayYield *XYldArray;
  static FluorCKCoef *FCKCoefArray;
  static AbsCoef *TotAbsCoefArray;
  static SPT *SPTArray;
  static SFCListElem *SFCList;
  static int nSFCList;
  static int lastNcall=-1; 
  static FILTER Filter;
  
  static EXP_PARAM ExpPar;
  static foil *Sample;
  static int NFoil, NFoilUsed, dimSums;
  static FILE *outputfile;
  static EXTRAINFO ExtraInfo;
  
  #if CPIXEVERBOSITY > 0
  printf("\n\n\n************ BEGIN OF CPIXE OUTPUT ****************\n\n");
  #endif

  pExpPar=&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)
    */
    

    //Calculate Filter Transmission if the filter changed
    if(Filter.changes){
      FilterTrans(pExpPar, XYldArray, TotAbsCoefArray, PresentElems, &Filter);
      ///@todo This indicates that the filter is not going to be varied (Provisional)
      Filter.changes=0;
    }

    initlyrarray(pExpPar, Sample, 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);
//    SFCorr(pExpPar, TotAbsCoefArray, FCKCoefArray, XYldArray ,NFoilUsed, &Filter, lyr, XYldSums);
      /******END CALCULATION**********/
    
      
    #if CPIXEVERBOSITY > 0  
    /**********BEGIN DEBUG LINES*******/
    /* Just write stuff for debugging*/
    printf("\n\n Ebeam=%lf IncAng=%lf  DetAng=%lf \n",
          pExpPar->BeamEner, pExpPar->IncAng/DEG2RAD, pExpPar->DetAng/DEG2RAD);
    printf("\n Charge=%lf   DTCC=%lf DetColFac=%lf\n",
          pExpPar->simpar.ColCharge, pExpPar->simpar.DTCC, pExpPar->DetColFac);

    printf("\n\nIteration %d\n",NCALL);
    
    for(i=0;i<NFoil;i++){
      printf("\nSample[%d]  (%d elem)\n",i,Sample[i].nfoilelm);
      printf("\n---Eff. thickness=%.3le  (%d sublyrs)\n",lyr[i+1].ThickIn,lyr[i+1].FESxlen);
      
      for(j=0;j<Sample[i].comp.nelem;j++){
      printf("'%s'(Z=%d) X[%d]=%lf (%lfat%% / %lfM%%)\n",
          XYldArray[Sample[i].comp.elem[j].Z].symb,Sample[i].comp.elem[j].Z,j,Sample[i].comp.X[j],100*Sample[i].comp.xn[j],100*Sample[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]);
            printf("\n%10s \texp:\t%le counts   (err=%.1le) ", "", CalcFlags[i].simareas.K_[j],CalcFlags[i].err.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]); 
              printf("\n%10s \texp:\t%le counts   (err=%.1le) ", "", CalcFlags[i].simareas.L_[j][l],CalcFlags[i].err.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]);
            printf("\n%10s \texp:\t%le counts   (err=%.1le) ", "", CalcFlags[i].simareas.M_[j],CalcFlags[i].err.M_[j]);
          }
      }
    }
//    getchar();
    /**********END DEBUG LINES*******/
    #endif  
    
    /**********BEGIN FILE OUTPUT LINES*******/
 
    if(ExtraInfo.WantOutputfile){
      fprintf(outputfile,"\n\n Ebeam=%lf IncAng=%lf  DetAng=%lf \n",
            pExpPar->BeamEner, pExpPar->IncAng/DEG2RAD, pExpPar->DetAng/DEG2RAD);
      fprintf(outputfile,"\n Charge=%lf   DTCC=%lf DetColFac=%lf\n",
            pExpPar->simpar.ColCharge, pExpPar->simpar.DTCC, pExpPar->DetColFac);
  
      fprintf(outputfile,"\n\nIteration %d\n",NCALL);
      
      for(i=0;i<NFoil;i++){
        fprintf(outputfile,"\nSample[%d]  (%d elem)\n",i,Sample[i].nfoilelm);
        fprintf(outputfile,"\n---Eff. thickness=%.3le  (%d sublyrs)\n",lyr[i+1].ThickIn,lyr[i+1].FESxlen);
        
        for(j=0;j<Sample[i].comp.nelem;j++){
        fprintf(outputfile,"'%s'(Z=%d) X[%d]=%lf (%lfat%% / %lfM%%)\n",
            XYldArray[Sample[i].comp.elem[j].Z].symb,Sample[i].comp.elem[j].Z,j,Sample[i].comp.X[j],100*Sample[i].comp.xn[j],100*Sample[i].comp.w[j]);
        }
        fprintf(outputfile,"\nFinal Energy in this layer=%lf keV\n\n",lyr[i+1].FoilOutEner);
      }
    
      for(i=0;i<Filter.nlyr;i++){
        fprintf(outputfile,"\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++){
        fprintf(outputfile,"(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]){
          fprintf(outputfile,"\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.) {
              fprintf(outputfile,"\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]);
              fprintf(outputfile,"\n%10s \texp:\t%le counts   (err=%.1le) ", "", CalcFlags[i].simareas.K_[j],CalcFlags[i].err.K_[j]);
            }
          /*L-Lines*/
          for (j = 1; j <= 3; j++)
            for (l = 1; l <= 3; l++) 
              if(XYldSums[i].L_[j][l] > 0.){
                fprintf(outputfile,"\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]); 
                fprintf(outputfile,"\n%10s \texp:\t%le counts   (err=%.1le) ", "", CalcFlags[i].simareas.L_[j][l],CalcFlags[i].err.L_[j][l]);
              }
          /*M lines:*/
          for (j = 1; j <= 3; j++) 
            if(XYldSums[i].M_[j] > 0.) {
              fprintf(outputfile,"\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]);
              fprintf(outputfile,"\n%10s \texp:\t%le counts   (err=%.1le) ", "", CalcFlags[i].simareas.M_[j],CalcFlags[i].err.M_[j]);
            }
        }
      }
    }
    /**********END FILE OUTPUT LINES*******/
      
    /*It is VERY important to call freeReusable() before starting a new iteration*/
    freeReusable(NFoil,&lyr,&Sample,&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.    
    */
    
    //initialize the file name pointers of the ExtraInfo structure
    ExtraInfo.SampleFileNm=samplefilenm;
    ExtraInfo.FilterFileNm=filterfilenm;
    ExtraInfo.XYldFileNm=calibfilenm;
    ExtraInfo.AreasFileNm=areasfilenm;
    ExtraInfo.DBpath=dbpath;
    ExtraInfo.OutputFileNm=outputfilenm;

    readINPUT(inputfilename, pExpPar, &ExtraInfo);
    readsample(samplefilenm, &pSimPar->MaxZinsample, &NFoil, &Sample);
    readFilter(filterfilenm, &Filter);
    //translate Filter definition                       
    Filter.changes=1;    //First time at least Filter should be read:
    
    snprintf(fckfilename, FILENMLENGTH,"%sSC_PRA74.fck",dbpath);
    snprintf(abscoeffilename, FILENMLENGTH,"%sxabs.abs",dbpath); 
  
    #if CPIXEVERBOSITY > 0  
    printf("\nInitializing. Reading databases...\n\n");
    printf("\nStopping Coefs DB: '%sSCOEF.95[A/B]'",dbpath);
    printf("\nAbsorption Coefs DB: '%s'",abscoeffilename);
    printf("\nFluorescence Coefs DB: '%s'",fckfilename);
    printf("\nDetector Calibration DB: '%s'",calibfilenm);
    printf("\nLines to calculate defined in: '%s'\n",areasfilenm);
    #endif

    createPresentElems(pSimPar->MaxZinsample, NFoil, Sample, &PresentElems);
    readCalcFlags(areasfilenm, PresentElems, pSimPar->MaxZinsample, ExtraInfo.AreasFormat, &CalcFlags);
    
    #if CPIXEVERBOSITY > 0
    for(i=1;i<pSimPar->MaxZinsample+1;i++){
      if(PresentElems[i]){
        printf("\nZ=%d  (%d)\n ",i,CalcFlags[i].atnum);
        for (j = 1; j <= 3; j++) printf("\t%.0lf",CalcFlags[i].simareas.K_[j]);
        for (j = 1; j <= 3; j++) { 
          printf("\n");
          for (l = 1; l <= 3; l++)printf("\t%.0lf",CalcFlags[i].simareas.L_[j][l]);
        }
        printf("\n");
        for (j = 1; j <= 3; j++) printf("\t%.0lf",CalcFlags[i].simareas.M_[j]);
      }
    }
    #endif
    
//     printf("\n\n!!!!!!!!!!");getchar();
    readXYld(calibfilenm, PresentElems, CalcFlags, 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);
    //If Sec Fluorescence is required, determine which elements produce fluorescence on which
    if(pSimPar->AllowSXFCorr){
      createSFCList(pExpPar, PresentElems, XYldArray, TotAbsCoefArray, &nSFCList, &SFCList);
//      prepareSFC(pExpPar, PresentElems);
      };
    
    //safefree((void*)&Sample); //Clear the initialized Matrix Array.
    //safefree((void*)&Filter.foil); //Clear the initialized Filter Foil.
    /********** END INITIALIZATION *********/

     //Write initialization parameters to output file
    if(ExtraInfo.WantOutputfile){
      outputfile = fopen(outputfilenm, "w");
      if(outputfile== NULL)
        {fprintf(stderr,"\nERROR: Could not open file '%s' for write.\n",outputfilenm);exit(1);}
      fprintf(outputfile,"\nCPIXE Output File\n\n");
      fprintf(outputfile,"\nStopping Coefs DB: '%sSCOEF.95[A/B]'",dbpath);
      fprintf(outputfile,"\nAbsorption Coefs DB: '%s'",abscoeffilename);
      fprintf(outputfile,"\nFluorescence Coefs DB: '%s'",fckfilename);
      fprintf(outputfile,"\nDetector Calibration DB: '%s'",calibfilenm);
      fprintf(outputfile,"\nLines to calculate defined in: '%s'\n",areasfilenm);
 
    }
  }
  
  else if(NCALL<0){
    /******** BEGIN TIDING UP BEFORE FINAL EXIT*********/
    /*Free everything for a clean exit*/
    #if CPIXEVERBOSITY > 0  
    printf("\nCleaning 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(CalcFlags);
    
    
    if(ExtraInfo.WantOutputfile)fclose(outputfile);

    /******** END TIDING UP BEFORE FINAL EXIT*********/
  }
  
  /*Store Number of Ncall*/
  lastNcall=NCALL;
  
  #if CPIXEVERBOSITY > 0
  printf("\n\n\n************ END OF CPIXE OUTPUT ****************\n\n");
  #endif
  
  return EXIT_SUCCESS;
}

---