decay.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
/*****************************************************************************/
#include "tpcisotope.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *decay_correction[] = {
"unknown",
"not corrected",
"corrected",
"corrected to scan start time",
"corrected to administration time",
0};
/*****************************************************************************/
 
/*****************************************************************************/
char *decayDescr(
  decaycorrection d
) {
  if(d>4) return(decay_correction[0]);
  return(decay_correction[d]);
}
/*****************************************************************************/
 
/*****************************************************************************/
double lambdaFromHalflife(
  double halflife
) {
  if(halflife>0.0) return(M_LN2/halflife);  // M_LN2=log(2.0)
  else return nan("");
}
/*****************************************************************************/
 
/*****************************************************************************/
double lambdaFromIsotope(
  int isotope
) {
  double hl=isotopeHalflife(isotope); if(isnan(hl) || hl<=0.0) return nan("");
  return(M_LN2/hl);
}
/*****************************************************************************/
 
/*****************************************************************************/
double decayCorrectionFactorFromLambda(
  double lambda,
  double starttime,
  double duration
) {
  if(starttime<0.0) return nan("");
  if(duration<0.0) duration=0.0;
  if(fabs(lambda)<1.0E-100) return nan("");
  double cf=exp(lambda*starttime);
  if(duration>1.0E-10) {
    double ff=fabs(lambda)*duration/(1.0-exp(-fabs(lambda)*duration));
    if(lambda<0.0) cf/=ff; else cf*=ff;
  }
  return(cf);
}
/*****************************************************************************/
 
/*****************************************************************************/
double decayCorrectionFactorFromIsotope(
  int isotope,
  double starttime,
  double duration
) {
  double lambda=lambdaFromIsotope(isotope); if(isnan(lambda)) return nan("");
  return decayCorrectionFactorFromLambda(lambda, starttime, duration);
}
/*****************************************************************************/
 
/*****************************************************************************/