tacnan.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include "tpcift.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
/*****************************************************************************/
#include "tpctac.h"
/*****************************************************************************/
 
/*****************************************************************************/
int tacXNaNs(
  TAC *tac
) {
  if(tac==NULL || tac->sampleNr<1) return(0);
  
  int fi, nan_nr=0;
  if(tac->isframe) {
    for(fi=0; fi<tac->sampleNr; fi++) {
      if(!isfinite(tac->x1[fi])) nan_nr++;
      if(!isfinite(tac->x2[fi])) nan_nr++;
    }
  } else {
    for(fi=0; fi<tac->sampleNr; fi++) {
      if(!isfinite(tac->x[fi])) nan_nr++;
    }
  }
  return(nan_nr);
}
/*****************************************************************************/
int tacYNaNs(
  TAC *tac,
  const int i
) {
  if(tac==NULL || tac->sampleNr<1 || tac->tacNr<1 || i>=tac->tacNr) return(0);
  
  int ri, fi, nan_nr=0;
  if(i>=0) {
    for(fi=0; fi<tac->sampleNr; fi++)
      if(!isfinite(tac->c[i].y[fi])) nan_nr++;
    return(nan_nr);
  }
  for(ri=0; ri<tac->tacNr; ri++)
    for(fi=0; fi<tac->sampleNr; fi++)
      if(!isfinite(tac->c[ri].y[fi])) nan_nr++;
  return(nan_nr);
}
/*****************************************************************************/
int tacNaNs(
  TAC *tac
) {
  return(tacXNaNs(tac)+tacYNaNs(tac, -1));
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacNotNaNs(
  TAC *tac,
  const int i
) {
  if(tac==NULL || tac->sampleNr<1 || tac->tacNr<1 || i>=tac->tacNr) return(0);
  
  int ri, fi, nan_nr=0;
  for(fi=0; fi<tac->sampleNr; fi++) {
    /* First, check x value(s) */
    if(tac->isframe) {
      if(!isfinite(tac->x1[fi])) {nan_nr++; continue;}
      if(!isfinite(tac->x2[fi])) {nan_nr++; continue;}
    } else {
      if(!isfinite(tac->x[fi])) {nan_nr++; continue;}
    }
    /* If ok, then check y value(s) */
    if(i>=0) {
      if(!isfinite(tac->c[i].y[fi])) nan_nr++; 
      continue;
    }
    for(ri=0; ri<tac->tacNr; ri++)
      if(!isfinite(tac->c[ri].y[fi])) {nan_nr++; break;}
  } // next sample
  return(tac->sampleNr-nan_nr);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacFixNaNs(
  TAC *tac
) {
  if(tac==NULL || tac->sampleNr<1 || tac->tacNr<1) return(1);
  if(tacXNaNs(tac)>0) return(2);
 
  int ri, fi, fj;
  double x1, x2, y1, y2, x, y;
  /* If NaNs with negative sample times, then set value to zero */
  for(fi=0; fi<tac->sampleNr; fi++) {
    if(tac->isframe) x=0.5*(tac->x1[fi]+tac->x2[fi]); else x=tac->x[fi];
    if(x>0.0) continue; 
    for(ri=0; ri<tac->tacNr; ri++) 
      if(!isfinite(tac->c[ri].y[fi])) tac->c[ri].y[fi]=0.0;
  }
  /* Go through the data, finding previous and next good value */
  for(ri=0; ri<tac->tacNr; ri++) if(tacYNaNs(tac, ri)>0) {
    x1=x2=y1=y2=nan("");
    for(fi=0; fi<tac->sampleNr; fi++) {
      if(tac->isframe) x=0.5*(tac->x1[fi]+tac->x2[fi]); else x=tac->x[fi];
      /* Save previous good point */
      if(isfinite(tac->c[ri].y[fi])) {y1=tac->c[ri].y[fi]; x1=x; continue;}
      /* We have missing value if we are here */
      /* Search for the next good point after this NaN */
      x2=y2=nan("");
      for(fj=fi+1; fj<tac->sampleNr; fj++) {
        if(!isfinite(tac->c[ri].y[fj])) continue;
        if(tac->isframe) x2=0.5*(tac->x1[fj]+tac->x2[fj]); else x2=tac->x[fj];
        y2=tac->c[ri].y[fj];
        break;
      }
      /* If no good points after current point, then that is bad */
      if(!isfinite(x2)) return(3);
      /* If the point before is not available, then use (0,0) */
      if(!isfinite(x1)) x1=y1=0.0;
      /* Now, interpolate value for current point */
      if(x2==x1) y=0.5*(y1+y2); else y=y2-(x2-x)*(y2-y1)/(x2-x1);
      tac->c[ri].y[fi]=y;
    } // next sample
  } // next TAC
  
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacNthSample(
  TAC *tac,
  const int sn,
  const int i
) {
  if(tac==NULL || tac->sampleNr<1 || tac->tacNr<1 || i>=tac->tacNr) return(-1);
  if(sn<1 || sn>tac->sampleNr) return(-1);
  
  int ri, fi, n=0;
  for(fi=0; fi<tac->sampleNr; fi++) {
    /* First, check x value(s) */
    if(tac->isframe) {
      if(!isfinite(tac->x1[fi]) || !isfinite(tac->x2[fi])) continue;
    } else {
      if(!isfinite(tac->x[fi])) continue;
    }
    /* If ok, then check y value(s) */
    if(i>=0) {
      if(!isfinite(tac->c[i].y[fi])) continue;
    } else {
      for(ri=0; ri<tac->tacNr; ri++) if(!isfinite(tac->c[ri].y[fi])) break;
      if(ri!=tac->tacNr) continue;
    }
    n++;
    if(n==sn) return(fi);
  } // next sample, if necessary
  return(-1);
}
/*****************************************************************************/
 
/*****************************************************************************/