cpixe/cpixe.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"


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

STATFILENAME inputfilename="input.in";

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

int cpixemain(int NCALL)
{
  extern STATFILENAME inputfilename;
  STATFILENAME fckfilename="";
  STATFILENAME abscoeffilename="";
  STATFILENAME xenerfilename="";


  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 CalibYld *LinesEnerArray;
  static CalibYld *EffArray;
  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);
      FilterTrans(pExpPar->simpar.MaxZinsample, LinesEnerArray, TotAbsCoefArray, PresentElems, &Filter);
      Filter.changes=0;
    }

//     initlyrarray(pExpPar, Sample, XYldArray, TotAbsCoefArray, SPTArray, PresentElems, NFoil, &NFoilUsed, &lyr);
    initlyrarray(pExpPar, Sample, LinesEnerArray, 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);
    integrate_Simpson2(pExpPar, TotAbsCoefArray, FCKCoefArray, EffArray ,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",
          ChemicalSymbol[Sample[i].comp.elem[j].Z],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",ChemicalSymbol[i],i);
        /*K lines:*/
        for (j = 1; j <= 3; j++)
          if(XYldSums[i].K_[j] > 0.) {
              fprintf(stdout,"\n%10s (%2.2lfkeV)\t%le counts   (Eff=%.3le Flt=%.2le)", LineNm[1][j], LinesEnerArray[i].K_[j], XYldSums[i].K_[j],EffArray[i].K_[j],Filter.Trans[i].K_[j]);
              fprintf(stdout,"\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(stdout,"\n%10s (%2.2lfkeV)\t%le counts   (cal=%.3le Flt=%.2le)", LineNm[j+1][l], LinesEnerArray[i].L_[j][l], XYldSums[i].L_[j][l],EffArray[i].L_[j][l],Filter.Trans[i].L_[j][l]);
                fprintf(stdout,"\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(stdout,"\n%10s (%2.2lfkeV)\t%le counts   (cal=%.3le Flt=%.2le)",LineNm[5][j], LinesEnerArray[i].M_[j], XYldSums[i].M_[j],EffArray[i].M_[j],Filter.Trans[i].M_[j]);
              fprintf(stdout,"\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",
            ChemicalSymbol[Sample[i].comp.elem[j].Z],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",ChemicalSymbol[i],i);
          /*K lines:*/
          for (j = 1; j <= 3; j++)
            if(XYldSums[i].K_[j] > 0.) {
              fprintf(outputfile,"\n%10s (%2.2lfkeV)\t%le counts   (Eff=%.3le Flt=%.2le)", LineNm[1][j], LinesEnerArray[i].K_[j], XYldSums[i].K_[j],EffArray[i].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], LinesEnerArray[i].L_[j][l], XYldSums[i].L_[j][l],EffArray[i].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], LinesEnerArray[i].M_[j], XYldSums[i].M_[j],EffArray[i].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 that never change within a fit.
    */

    readINPUT(inputfilename, pExpPar, &ExtraInfo);
    /*Disallow the DEPRECATED CalibYld Method*/
    if(!ExtraInfo.useefficiencies){
      fprintf(stderr,"\nERROR: The DATTPIXE calibration scheme is now deprecated. Use Efficiencies instead.\n");exit(1);}

    readsample(ExtraInfo.SampleFileNm, &pSimPar->MaxZinsample, &NFoil, &Sample);
    readFilter(ExtraInfo.FilterFileNm, &Filter);
    //translate Filter definition
    Filter.changes=1;    //First time at least Filter should be read:

    snprintf(fckfilename, FILENMLENGTH,"%sSC_PRA74.fck",ExtraInfo.DBpath);
    snprintf(abscoeffilename, FILENMLENGTH,"%sxabs.abs",ExtraInfo.DBpath);
    snprintf(xenerfilename, FILENMLENGTH,"%sxener.dat",ExtraInfo.DBpath);

    #if CPIXEVERBOSITY > 0
    printf("\nInitializing. Reading databases...\n\n");
    printf("\nStopping Coefs DB: '%sSCOEF.95[A/B]'",ExtraInfo.DBpath);
    printf("\nAbsorption Coefs DB: '%s'",abscoeffilename);
    printf("\nFluorescence Coefs DB: '%s'",fckfilename);
    printf("\nX-ray Energies DB: '%s'",xenerfilename);
    printf("\nDetector Calibration DB: '%s'",ExtraInfo.CalibFileNm);
    printf("\nLines to calculate defined in: '%s'\n",ExtraInfo.AreasFileNm);
    #endif

    createPresentElems(pSimPar->MaxZinsample, NFoil, Sample, &PresentElems);
    readCalcFlags(ExtraInfo.AreasFileNm, PresentElems, pSimPar->MaxZinsample, ExtraInfo.AreasFormat, &CalcFlags);

    //Read Efficiencies for the lines used
    //2006/04/05: changed this to get the efficiencies from a list (see below)
    //readEff(ExtraInfo.CalibFileNm, PresentElems, pSimPar->MaxZinsample, &EffArray);

    //Read the energies for the lines used
    readLinesEner(xenerfilename, PresentElems, pSimPar->MaxZinsample, &LinesEnerArray);

    //2006/04/05:Read the efficiencies for the lines used (from a two-column "Energy|Effi" formated file)
    //readEff2(ExtraInfo.CalibFileNm, PresentElems, pSimPar->MaxZinsample, LinesEnerArray, &EffArray);



    readFCK(fckfilename, pSimPar->MaxZinsample,&FCKCoefArray);
//     readAbsCoef_OLD(abscoeffilename, pSimPar->MaxZinsample, PresentElems, &Filter, &TotAbsCoefArray);
    readAbsCoef(abscoeffilename, pSimPar->MaxZinsample, &TotAbsCoefArray);
    createSPTs(ExtraInfo.DBpath,pExpPar, pSimPar->MaxZinsample, PresentElems, (double)(2*pExpPar->ion.A), &SPTArray);

    //2006/04/07:Read the efficiencies with readeff3
    readEff3(ExtraInfo.CalibFileNm, PresentElems, pSimPar->MaxZinsample, LinesEnerArray, &EffArray);

    //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(ExtraInfo.OutputFileNm, "w");
      if(outputfile== NULL)
        {fprintf(stderr,"\nERROR: Could not open file '%s' for write.\n",ExtraInfo.OutputFileNm);exit(1);}
      fprintf(outputfile,"\nCPIXE Output File\n\n");
      fprintf(outputfile,"\nStopping Coefs DB: '%sSCOEF.95[A/B]'",ExtraInfo.DBpath);
      fprintf(outputfile,"\nAbsorption Coefs DB: '%s'",abscoeffilename);
      fprintf(outputfile,"\nFluorescence Coefs DB: '%s'",fckfilename);
      fprintf(outputfile,"\nX-ray Energies DB:'%s'",xenerfilename);
      fprintf(outputfile,"\nDetector Calibration DB: '%s'",ExtraInfo.CalibFileNm);
      fprintf(outputfile,"\nLines to calculate defined in: '%s'\n",ExtraInfo.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(EffArray);
    free(LinesEnerArray);
    free(FCKCoefArray);
    free(TotAbsCoefArray);
    free(CalcFlags);


    if(ExtraInfo.WantOutputfile)fclose(outputfile);
    freeExtraInfo(&ExtraInfo);

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

---