absscal.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <ctype.h>
#include <time.h>
#include <math.h>
/*****************************************************************************/
#include "libtpcmisc.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Extraction and calibration of blood TACs measured using",
  "automatic blood sampling systems (ABSS, \"blood pump\") (1).",
  "Application name was previously blo2kbq, version 3.8.1.",
  " ",
  "Usage: @P [Options] ABSS_file [BTAC_file]",
  " ",
  "Options:",
  " -c=<filename>",
  "     Calibration coefficients for ABSS and well-counter should be specified",
  "     with calibration file (2), given with option -c.",
  " -i=<isotope>",
  "     If ABSS data file does not contain the isotope code or halflife, then",
  "     it can be given with this option, using codes",
  "     C-11, F-18, Ga-68, Ge-68, O-15, ...",
  " -decay=<on|off>",
  "     Blood data is corrected for radioactive decay (on, default), or not",
  "     corrected (off).",
//"     Possible background subtraction is done before decay correction.",
  " -start=<t>",
  "     If ABSS measurement was not started at the time of injection,",
  "     the delay t (sec) must be entered with option -start=<t>. If data",
  "     collection was started before injection, enter a negative value for t;",
  "     activities before that time are removed from the output data.",
  " -density[=<blood density>]",
  "     With this option blood radioactivity concentrations (kBq/mL) are",
  "     divided by the given blood density (by default 1.06)",
  "     to get concentrations in units kBq/g.",
  " -min",
  "     Sample times are written in minutes (in seconds by default).",
  " -mean, -left, -right, -both",
  "     The counts or radioactivity concentrations are by default calculated",
  "     as the mean of the channels (-mean). With options -left and -right ",
  "     only the left (ch2) or right (ch1) channel is used.",
  "     With option -both the results from both channels are saved.",
  "     Only applicable to Scanditronics and GEMS data.",
/*
  " -z",
  "     For reading of the measurement date, local time zone can be used with",
  "     option -z instead of Greenwich main time (default).",
*/
  " -abss=<device>",
  "     ABSS produces device-specific file formats, and detectors have",
  "     calibration factors. This option is needed if ABSS can not be correctly",
  "     determined from the filename (*.lis, *.bld, *.alg, or *.txt).",
  "     Accepted devices are Scanditronics, GEMS, Allogg1, and Allogg2, or",
  "     numbers 1-4, accordingly.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Filename for corrected blood TAC data can be entered as the last command-line",
  "argument. If not given, then file is written to the directory where ABSS file",
  "is; calibrated file is named *.kbq and non-calibrated *.dat by default.",
  " ",
  "Example:",
  "  @P -i=O-15 -c=S:/Lab/plasma/bsampler_calibration/pump_cal.dat sr452_blo.bld",
  " ",
//"In case of the GEMS online sampler (GE Advange/PET CT), the injection",
//"time is always taken from the first sample time in the *.bld file,",
//"not from the header of the file, like with Scanditronics data.",
//" ",
  "References:",
  "1. https://www.turkupetcentre.net/petanalysis/input_abss.html",
  "2. https://www.turkupetcentre.net/petanalysis/input_abss_calibration.html",
  " ",
  "See also: absszero, absstime, abssexam, abssfch, taccat, disp4dft, fitdelay",
  " ",
  "Keywords: blood, input, ABSS, calibration",
  0};
/*****************************************************************************/
 
/*****************************************************************************/
/* Types and defaults */
#define MAX_PUMPS 4
#define BLOOD_DENSITY 1.06
typedef enum {SCANDITRONICS,ALLOGG,ALLOGG2} Format;
static int both_cols=0;
static int gm_timez=1; /*time zone = Greenwich main time*/
static int pump=0;
char isotope[10]="UNKNOWN";
double halflife=-1.0;
char studynr[MAX_STUDYNR_LEN+1];
/*****************************************************************************/
/* Results data */
struct BDataResults {
  double t, ch1, ch2;
};
/*****************************************************************************/
/* Scanditronics data */
struct BDataS {
  struct BDataResults results;
  double time_of_day;
  double time_since_start;
  double meas_interval;
  int pair1_coincid;
  int pair1_det1_cnts;
  int pair1_det2_cnts;
  int pair2_coincid;
  int pair2_det1_cnts;
  int pair2_det2_cnts;
  int aux_cnts;
};
 
/* Allogg data */
struct BDataA {
  struct BDataResults results;
  double time_since_start;
  double meas_interval; /* This is calculated from the data */
  int data[2];
};
 
struct DataList {
  union {
    struct BDataS scanditronics;
    struct BDataA allogg;
    struct BDataResults results;
  };
  struct DataList *prev, *next;
};
 
struct Blood {
  struct tm date;
  time_t time;
  char bgdate[11];
  float background;    
  struct DataList *data;
};
/*****************************************************************************/
/* Local functions */
int read_scanditronics_data(
  const char *inputfile, struct Blood *blood, int verbose);
int read_allogg_data(
  const char *inputfile, struct Blood *blood, int verbose);
int read_allogg2_data(
  char *inputfile, struct Blood *blood, int verbose, char *msg);
void correct_scanditronics(struct Blood *blood, int verbose);
int correct_allogg(struct Blood *blood);
int correct_allogg2(struct Blood *blood, int verbose);
void subtract_background(struct Blood *blood);
int get_calfactors(const char *file,time_t bloodt, double *wcf,
  double *pumpf, int pump, char *isotope, char *caldate, int verbose);
int find_date(char *str, struct Blood *blood, int verbose);
int find_halflife(char *str, int verbose);
/*****************************************************************************/
 
/*****************************************************************************/
/* Turn on the globbing of the command line, since it is disabled by default in
   mingw-w64 (_dowildcard=0); in MinGW32 define _CRT_glob instead, if necessary;
   In Unix&Linux wildcard command line processing is enabled by default. */
/*
#undef _CRT_glob
#define _CRT_glob -1
*/
int _dowildcard = -1;
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char *argv[])
{
  int ai, help=0, version=0, verbose=1;
 
  double startTime=0.0, density=1.0, wc_factor=-1.0, pump_factor=-1.0;
  int decay_correct=1; // 0=no, 1=yes
  char pumpfile[FILENAME_MAX], calfile[FILENAME_MAX], outfile[FILENAME_MAX];
  char tmp[FILENAME_MAX];
  float sum1, sum2;
  Format format=SCANDITRONICS;
  struct DataList *ptr;
  int timeInMinutes=0;
  struct Blood blood;
  char caldate[25];
  int column=0; // use mean (0), left (1), or right (2) channel
  char *cptr;
  FILE *fp;
  int r;
  int zero1=1, zero2=1;
  time_t t_injection, t_samp;
  struct tm tm_injection, tm_samp;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  pumpfile[0]=calfile[0]=outfile[0]=(char)0;
  studynr[0]=(char)0;
 
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    //printf("argv[%d] := '%s'\n", ai, argv[ai]);
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strncasecmp(cptr, "C=", 2)==0 && strlen(cptr)>2) {
      strlcpy(calfile, cptr+2, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "O=", 2)==0 && strlen(cptr)>2) { // deprecated
      strlcpy(outfile, cptr+2, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "I=", 2)==0 && strlen(cptr)>2) {
      if(hlCorrectIsotopeCode(cptr+2)!=NULL) {
        strcpy(isotope, hlCorrectIsotopeCode(cptr+2));
        halflife=hlFromIsotope(isotope);
        continue;
      }
      printf("Error: Incorrect isotope code '%s'\n", cptr+2);
      return(1);
    } else if(strncasecmp(cptr, "DECAY=", 6)==0 && strlen(cptr)>6) {
      cptr+=6;
      if(strcasecmp(cptr, "OFF")==0) {decay_correct=0; continue;}
      if(strcasecmp(cptr, "ON")==0) {decay_correct=1; continue;}
    } else if(strcasecmp(cptr, "DENSITY")==0) {
      density=BLOOD_DENSITY; continue;
    } else if(strncasecmp(cptr, "DENSITY=", 8)==0) {
      cptr+=8; density=atof_dpi(cptr); if(density>0.0) continue;
    } else if(strncasecmp(cptr, "MIN", 1)==0) {
      timeInMinutes=1; continue;
    } else if(strcasecmp(cptr, "SEC")==0) {
      timeInMinutes=0; continue;
    } else if(strncasecmp(cptr, "Z", 1)==0) {
      gm_timez=0; continue;
    } else if(strncasecmp(cptr, "LEFT", 1)==0) {
      column=1; continue;
    } else if(strncasecmp(cptr, "RIGHT", 1)==0) {
      column=2; continue;
    } else if(strcasecmp(cptr, "BOTH")==0) {
      both_cols=1; continue;
    } else if(strncasecmp(cptr, "PUMP=", 5)==0 && strlen(cptr)>5) {
      cptr+=5;
      if(strncasecmp(cptr, "SCANDITRONICS", 1)==0) {format=SCANDITRONICS; pump=1; continue;}
      if(strncasecmp(cptr, "GEMS", 1)==0) {format=SCANDITRONICS; pump=2; continue;}
      if(strncasecmp(cptr, "ALLOGG1", 1)==0) {format=ALLOGG; pump=3; continue;}
      if(strcasecmp(cptr, "ALLOGG2")==0) {format=ALLOGG2; pump=4; continue;}
      pump=atoi(cptr); 
      if(pump==1 || pump==2) {format=SCANDITRONICS; continue;}
      if(pump==3) {format=ALLOGG; continue;}
      if(pump==4) {format=ALLOGG2; continue;}
    } else if(strncasecmp(cptr, "START=", 6)==0 && strlen(cptr)>6) {
      startTime=atof_dpi(cptr+6); continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
  
  /* Print help or version? */
  if(help==2) {tpcHtmlUsage(argv[0], info, ""); return(0);}
  if(help) {tpcPrintUsage(argv[0], info, stdout); return(0);}
  if(version) {tpcPrintBuild(argv[0], stdout); return(0);}
 
  /* Process other arguments, starting from the first non-option */
  if(ai<argc) strlcpy(pumpfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(outfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {fprintf(stderr, "Error: too many arguments: '%s'.\n", argv[ai]); return(1);}
 
  /* Is something missing? */
  if(!pumpfile[0]) {
    fprintf(stderr, "Error: missing ABSS file name.\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("pumpfile := %s\n", pumpfile);
    printf("isotope := %s\n", isotope);
    printf("halflife := %g\n", halflife);
    printf("decay_correct := %d\n", decay_correct);
    printf("startTime := %g\n", startTime);
    printf("density := %g\n", density);
    printf("timeInMinutes := %d\n", timeInMinutes);
    printf("column := %d\n", column);
    printf("both_cols := %d\n", both_cols);
    fflush(stdout);
  }
 
  /* Decide which pump we are dealing with, if not given by user */
  if(pump==0) {
    cptr=pumpfile+strlen(pumpfile)-4;
    if(strcasecmp(cptr, ".lis")==0) {pump=1; format=SCANDITRONICS;}
    else if(strcasecmp(cptr, ".bld")==0) {pump=2; format=SCANDITRONICS;}
    else if(strcmp(cptr, ".alg")==0) {pump=3; format=ALLOGG;}
    else if(strcmp(cptr, ".txt")==0) {pump=4; format=ALLOGG2;}
    else {
      fprintf(stderr, "Error: cannot identify the pump;\n");
      fprintf(stderr, "please rename the file to *.lis, *.bld, *.alg or *.txt\n");
      fprintf(stderr, "or use option -abss=<id>.\n");
      return(1);
    }
  }
  if(verbose>1) {
    printf("pump := %d\n", pump);
    printf("format := %d\n", format);
  }
 
 
  /* 
   *  Read in data
   */
  if(verbose>0) printf("reading %s\n", pumpfile);
  switch(format) {
    case SCANDITRONICS:
      if(read_scanditronics_data(pumpfile, &blood, verbose-2)) return 2;
      break;
    case ALLOGG:
      if(read_allogg_data(pumpfile, &blood, verbose-2)) return 2;
      break;
    case ALLOGG2:
      if(read_allogg2_data(pumpfile, &blood, verbose-2, tmp)) {
        fprintf(stderr, "Error in reading %s: %s.\n", pumpfile, tmp);
        return 2;
      }
      break;
  }
  /* Check that isotope is specified with decay correction */
  if(decay_correct!=0) {
    if(halflife<=0.0 || strcasecmp(isotope, "UNKNOWN")==0) {
      fprintf(stderr,"Error: isotope must be specified for decay correction.\n");
      return 1;
    }
  }
  /* Check that isotope is specified with calfile */
  if(calfile[0] && strcasecmp(isotope, "UNKNOWN")==0) {
    fprintf(stderr,"Error: isotope must be specified for calibration.\n");
    return 1;
  }
 
  /* Process the data */
  if(verbose>0) printf("processing the data\n");
  switch(format) {
    case SCANDITRONICS:
      correct_scanditronics(&blood, verbose-2);
      break;
    case ALLOGG:
      if(correct_allogg(&blood) == 1) return 3;
      break;
    case ALLOGG2:
      if(correct_allogg2(&blood, verbose-2) == 1) return 3;
      break;
  }
 
  /* Check that neither of the channels has all values zero*/
  if(verbose>1) printf("verifying that data is not all zeroes\n");
  zero1=zero2=1;
  if(format==SCANDITRONICS) { 
    ptr=blood.data; 
    while(ptr) {
      if(fabs(ptr->results.ch1) > 1.0e-10) zero1=0;
      if(fabs(ptr->results.ch2) > 1.0e-10) zero2=0;
      ptr=ptr->next;
    }
    if(zero1){
      fprintf(stderr, "Error: All values from channel 1 are zeros. Contact the physicist!\n");
      if(column!=2) return 4;
    }
    if(zero2){
      fprintf(stderr, "Error: All values from channel 2 are zeros. Contact the physicist!\n");
      if(column!=1) return 4;
    }
  }
 
  /* Subtract background radiation; must be done before decay correction */
  if(blood.background>0.0) {
    if(verbose>=0) printf("subtracting background\n"); 
    subtract_background(&blood);
  } else if(blood.background<-1.0) {
    if(format==ALLOGG || format==ALLOGG2) {
      if(calfile[0]) {
        fprintf(stderr, "Error: %s does not contain background.\n", pumpfile);
        return(5);
      } else
        fprintf(stderr, "Warning: %s does not contain background.\n", pumpfile);        
    }
  }
 
 
  /* 
   *  Calibration 
   */
  if(calfile[0]) {
  
    /* Check that we have the measurement date */
    if(blood.date.tm_year==0) {
      fprintf(stderr, "\nError: Blood file must contain the measurement date!\n");
      fprintf(stderr, "Open the blood file in a text editor and add the date to\n");
      fprintf(stderr, "the beginning as its own comment line, e.g. # 1999-11-31\n\n");
      return(6);
    }
  
    if(verbose>0) printf("calibrating\n");
    ptr=blood.data;
    if(get_calfactors(calfile, timegm(&blood.date), &wc_factor, &pump_factor,
                      pump, isotope, caldate, verbose-1))
//    if(get_calfactors(calfile,mktime(&blood.date),&wc_factor,&pump_factor,
//                      pump,isotope,caldate))
      return 7;
    if(verbose>=0)
      printf("  PC := %g\n  WC/BR := %g\n  calibration_date := %s\n",
        pump_factor, wc_factor, caldate);
    while(ptr) {
      ptr->results.ch1*=wc_factor*pump_factor;
      ptr->results.ch2*=wc_factor*pump_factor;
      ptr=ptr->next;
    }
  } else if(verbose>=0) {
    printf("Warning: no calibration file was entered; data is not calibrated.\n");
  }
 
 
  /* Correct the start time for decay correction */
  /* If negative, then correct them back again later */
  /* Set earlier activities to 0 */
  if(startTime!=0.0) {
    if(verbose>1) printf("correcting the start time\n");
    ptr=blood.data;
    while(ptr) {
      ptr->results.t+=startTime;
      if(format==SCANDITRONICS)
        ptr->scanditronics.time_since_start+=startTime;
      else
        ptr->allogg.time_since_start+=startTime;
      if(ptr->results.t<0.0) {
        ptr->results.ch1=0;
        ptr->results.ch2=0;
      }
      ptr=ptr->next;
    }
  }
 
  /* Print the time of the first sample */
  if(format==SCANDITRONICS) {
    ptr=blood.data; t_samp=ptr->scanditronics.time_of_day;
    if(gmtime_r(&t_samp, &tm_samp)==NULL) {
      fprintf(stderr, "Warning: invalid sample time.\n");
    } else {
      tm_samp.tm_year=blood.date.tm_year;
      tm_samp.tm_mon=blood.date.tm_mon;
      tm_samp.tm_mday=blood.date.tm_mday;
      t_samp=timegm(&tm_samp);
      char buf[32];
      if(!ctime_r_int(&t_samp, buf)) {
        fprintf(stderr, "Warning: invalid sample time.\n");
      } else {
        fprintf(stdout, "Time point of the first sample is: %s\n", buf);
      }
    }
  }
 
  /* Decay correction */
  if(decay_correct!=0) {
    double lambda=hl2lambda(halflife*60.0);
    double dc=0;
 
    ptr=blood.data;
    if(verbose>=0) printf("decay correction with half-life %.2f min.\n", halflife);
    while(ptr) {
      dc=1.0;
      if(format==SCANDITRONICS) {
        if(ptr->scanditronics.time_since_start>=0.0) {
          dc = hlLambda2factor(lambda,
                 ptr->scanditronics.time_since_start, ptr->scanditronics.meas_interval);
        }
      } else {
        if(ptr->allogg.time_since_start>=0.0) {
          dc = hlLambda2factor(lambda, ptr->allogg.time_since_start, ptr->allogg.meas_interval);
          if(verbose>5) printf("dc:%f %f\n",dc,ptr->allogg.time_since_start);
        }
      }
      ptr->results.ch1*=dc;
      ptr->results.ch2*=dc;
      ptr=ptr->next;
    }
  }
 
#if(0)  // REMOVED IN BLO2KBQ VERSION 3.4.3
  /* If start time is negative, then set orig times */
  if(startTime<0.0) {
    ptr=blood.data;
    while(ptr) {
      ptr->results.t-=startTime;
      ptr=ptr->next;
    }
  }
#endif
 
  /* Blood density correction */
  if(density!=1.0) {
    if(verbose>0)
      printf("kBq/mL to kBq/g correction with blood density of %.2f g/mL.\n", density);
    ptr=blood.data;
    while(ptr) {
      ptr->results.ch1/=density;
      ptr->results.ch2/=density;
      ptr=ptr->next;
    }
  }
 
  /* Check that channels are not too different */
  ptr=blood.data; sum1=sum2=0.0;
  while(ptr) {
    if(format==ALLOGG || format==ALLOGG2) {
      if(ptr->results.ch1 < -2.0)
        printf("Warning: channel below zero:%lf\n", ptr->results.ch1);
    } else {
      if(ptr->results.ch1 < 0)
        printf("Warning: channel one below zero:%lf\n", ptr->results.ch1);
      if(ptr->results.ch2 < 0)
        printf("Warning: channel two below zero:%lf\n", ptr->results.ch2);
    }
    sum1+=ptr->results.ch1;
    sum2+=ptr->results.ch2;
    ptr=ptr->next;
  }
  if(column!=1 && column!=2 && (fabs(sum1-sum2)/(1.0+sum1+sum2))>0.3)
    fprintf(stderr, "Warning: Large difference between channels of online blood sampler!\n");
 
 
  /*
   *  Write output 
   */
 
  /* Create output filename, if not given */
  if(!outfile[0]) {
    if(studynr_validity_check2(studynr, 1)) strcpy(outfile, studynr);
    else strcpy(outfile, pumpfile);
    cptr=strrchr(outfile,'.'); if(cptr!=NULL) *cptr='\0';
    if(calfile[0]) strcat(outfile,".kbq"); else strcat(outfile,".dat");
    if(verbose>1) printf("outputfile := %s\n", outfile);
  }
  /* Check if output file exists; backup, if necessary */
  if(backupExistingFile(outfile, NULL, tmp)) {
    fprintf(stderr, "Error: %s\n", tmp);
  }
  /* Open the file */
  if(verbose>=0) printf("writing output in file %s\n", outfile);
  fp=fopen(outfile, "w");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot write %s\n", outfile);
    perror("fopen");
    return 11;
  }
 
  /*
   *  Write header information
   */
  if(verbose>2) printf("writing 'header'\n");
 
  /* measurement start time */
  if(gmtime_r(&t_injection, &tm_samp)==NULL) {
    fprintf(stderr, "Warning: invalid measurement start time.\n");
  }
#if(0)
  } else {
    /* Greenwich main time->Hki time zone*/
    if(blood.date.tm_hour==tm_samp.tm_hour-2) 
      blood.date.tm_hour+=2; 
    else if(blood.date.tm_hour==tm_samp.tm_hour-3) 
      blood.date.tm_hour+=3; 
  }
#endif
  fprintf(fp, "# measurement_start_time := %.4d-%.2d-%.2d %.2d:%.2d:%.2d\n",
    blood.date.tm_year+1900,blood.date.tm_mon+1,blood.date.tm_mday,
    blood.date.tm_hour,blood.date.tm_min,blood.date.tm_sec);
  /* detector id */
  fprintf(fp, "# detector := %d\n", pump);
  /* sample time unit */
  if(timeInMinutes==1) fprintf(fp, "# time_unit := min\n");
  else fprintf(fp, "# time_unit := sec\n");
  /* calibration */
  if(calfile[0]) {
    fprintf(fp, "# calibrated := yes\n");
    fprintf(fp, "# unit := kBq/ml\n");
    fprintf(fp, "# calibration_date := %s\n", caldate);
    fprintf(fp, "# pump_calibration_factor := %g\n", pump_factor);
    fprintf(fp, "# wc_calibration_factor_per_br := %g\n", wc_factor);
    fprintf(fp, "# calibration_file := %s\n", calfile);
  } else {
    fprintf(fp, "# calibrated := no\n");
    fprintf(fp, "# unit := cps\n");
  }
  /* isotope */
  fprintf(fp, "# isotope := %s\n", isotope);
  /* decay correction */
  if(decay_correct!=0) {
    fprintf(fp, "# decay_correction := yes\n");
    fprintf(fp, "# halflife := %g [min]\n", halflife);
  } else {
    fprintf(fp, "# decay_correction := no\n");
  }
  /* injection time */
  r=3600*blood.date.tm_hour + 60*blood.date.tm_min + blood.date.tm_sec;
  if(startTime==0.0) { /* injection time = measurement start time */
    fprintf(fp, "# injection_time := %.4d-%.2d-%.2d %.2d:%.2d:%.2d\n",
      blood.date.tm_year+1900, blood.date.tm_mon+1, blood.date.tm_mday,
      blood.date.tm_hour, blood.date.tm_min, blood.date.tm_sec);
  } else if(r==0) {
    fprintf(fp, "# injection_time := %g [sec]\n", startTime);
  } else {
    t_injection=timegm(&blood.date);
    if(t_injection==-1) {
      fprintf(stderr, "Error: invalid injection time.\n");
      fclose(fp); return(12);
    } else {
      t_injection-=startTime;
      if(gmtime_r(&t_injection, &tm_injection)==NULL) {
        fprintf(stderr, "Error: invalid injection time.\n");
        fclose(fp); return(12);
      }
    }
    fprintf(fp, "# injection_time := %.4d-%.2d-%.2d %.2d:%.2d:%.2d\n",
      tm_injection.tm_year+1900, tm_injection.tm_mon+1, tm_injection.tm_mday,
      tm_injection.tm_hour, tm_injection.tm_min, tm_injection.tm_sec);
  }
  /* background activity */
  if(blood.background>0.0) {
    if(blood.bgdate[0])
      fprintf(fp,"# background_date := %s\n", blood.bgdate);
    fprintf(fp,"# background := %g [cps]\n", blood.background);
  }
  /* Study number */
  if(studynr_validity_check2(studynr, 1)) {
    fprintf(fp, "# studynr := %s\n", studynr);
  }
 
  /*
   *  Write the data lines
   */
  if(verbose>2) printf("Writing data lines\n");
  ptr=blood.data;
  while(ptr) {
    /* Don't print if sample time is negative */
    if(ptr->results.t<0.0) {ptr=ptr->next; continue;}
    /* Write the sample time */
    if(timeInMinutes) fprintf(fp, "%12.4f", ptr->results.t/60.0);
    else fprintf(fp, "%12.3f", ptr->results.t);
    /* Write the channel avg, or both separately, or only one */
    if(both_cols) {
      fprintf(fp, " %12.3f %12.3f\n", ptr->results.ch1, ptr->results.ch2);
    } else if(column==1) {
      fprintf(fp, " %12.3f\n", ptr->results.ch1);
    } else if(column==2) {
      fprintf(fp, " %12.3f\n", ptr->results.ch2);
    } else {
      fprintf(fp, " %12.3f\n", 0.5*(ptr->results.ch1+ptr->results.ch2));
    }
    ptr=ptr->next;
  }
  
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int find_date(
  char *str,
  struct Blood *blood,
  int verbose
) {
  int c,yyyy,MM,dd,hh=0,mm=0,ss=0;
  if(str[6]=='-' && (str[8]=='-'||str[9]=='-')) {
    c=sscanf(str,"# %d-%d-%d %d:%d:%d",&yyyy,&MM,&dd,&hh,&mm,&ss);
    /* Found all fields */
    if(c==6 || c==3) { 
      blood->date.tm_year=yyyy-1900;
      blood->date.tm_mday=dd;
      blood->date.tm_mon=MM-1;
      blood->date.tm_hour=hh;
      blood->date.tm_min=mm;
      blood->date.tm_sec=ss;
      blood->date.tm_isdst=-1;
      if(verbose>0) 
        printf("find_date(): %04d-%02d-%02d %02d:%02d:%02d\n", yyyy, MM, dd, hh, mm, ss);
      return 0;
    }
  }
  return 1;
}
/*****************************************************************************/
 
/*****************************************************************************/
int find_halflife(
  char *str,
  int verbose
) {
  char *cptr=strcasestr(str, "half-life: ");
  if(cptr!=NULL) {
    cptr+=11; 
  } else {
    cptr=strcasestr(str, "halflife: ");
    if(cptr!=NULL) {
      cptr+=10; 
    } else {
      return(1);
    }
  }
  if(verbose>0) printf("%s(): %s\n", __func__, cptr);
  double f=atof_dpi(cptr); if(!isnormal(f)) return(1);
  if(!isnan(halflife) && halflife>0.0 && !doubleMatch(halflife, f, 0.1)) {
    fprintf(stderr, "Error: different isotope in file and given by user.\n");
    exit(1);
  }
  halflife=f;
  int c=hlIsotopeFromHalflife(halflife); if(c<0) return(1);
  strcpy(isotope, hlIsotopeCode(c));
  return(1);
}
/*****************************************************************************/
 
/*****************************************************************************/
int find_background(
  const char *str,
  struct Blood *blood
) {
  char tmp[11];
 
  while(isspace((int)*str) || *str=='#') str++;
  if(*str) {
    if(strncasecmp(str, "Background:", 11)==0) {
      str+=11;
      if(sscanf(str, "%f", &blood->background) == 1) {
        return 0;
      }
    } else if(strncasecmp(str, "Background on", 13)==0) {
      str+=13;
      if(sscanf(str,"%8s: %f", tmp, &blood->background) == 2) {
        sprintf(blood->bgdate, "%4.4s-%2.2s-%2.2s", tmp, tmp+4, tmp+6);
        return 0;
      } else if(sscanf(str, "%10s: %f", blood->bgdate, &blood->background)==2) {
        return 0;
      }
    }
  }
  return 1;
}
/*****************************************************************************/
 
/*****************************************************************************/
int read_scanditronics_data(
  const char *inputfile,
  struct Blood *blood,
  int verbose
) {
  char line[512],*str;
  struct DataList *points;
  struct BDataS *bdata;
  int datenotfound=1;
  int hlnotfound=1;
  FILE *fp;
  int n;
  //time_t t;
 
  if(verbose>0) printf("%s(%s)\n", __func__, inputfile);
  blood->bgdate[0]='\0';
  blood->background = 0.0;
 
  fp=fopen(inputfile,"r");
  if(!fp) {
    fprintf(stderr,"Error: cannot open \"%s\".\n",inputfile);
    perror("fopen");
    return 1;
  }
  /*
  t=0; if(gm_timez) blood->date=*gmtime(&t); else blood->date=*localtime(&t);
  */
  points=malloc(sizeof(struct DataList));
  points->prev=NULL; points->next=NULL;
  bdata=&points->scanditronics;
  memset(bdata,0,sizeof(struct BDataS));
  blood->data=points;
  /* Read data */
  while(fgets(line,sizeof(line),fp)) {
    str=line;
    while(*str && isspace((unsigned char)*str)) str++;
    if(*str=='#') { /* Comment line; find date */
      if(datenotfound) datenotfound=find_date(str, blood, verbose);
      if(hlnotfound) hlnotfound=find_halflife(str, verbose);
      continue;
    }
    if(datenotfound) {
      fprintf(stderr, "Error: valid measurement date was not found in header.\n");
      fclose(fp); return 5;
    }
    if(*str=='*' || *str=='\0') continue;
    n=sscanf(str, "%lf %lf %lf %i %i %i %i %i %i %i",
          &bdata->time_of_day, &bdata->time_since_start, &bdata->meas_interval, 
          &bdata->pair1_coincid, &bdata->pair1_det1_cnts, &bdata->pair1_det2_cnts,
          &bdata->pair2_coincid, &bdata->pair2_det1_cnts, &bdata->pair2_det2_cnts, 
          &bdata->aux_cnts);
    if(n==10) {
      points->next=malloc(sizeof(struct DataList));
      points->next->prev=points;
      points=points->next;
      points->next=NULL;
      bdata=&points->scanditronics;
      memset(bdata,0,sizeof(struct BDataS));
    } else {
      if(n>0 && n!=7) fprintf(stderr,"Warning: malformed line!:{%s}\n",str);
    }
  }
  fclose(fp);
 
  if(bdata->time_of_day==0 && bdata->meas_interval==0) {
    /* Remove extra entry if any */
    struct DataList *prev=points->prev;
    free(points);
    points=prev;
    if(points) {
      points->next=NULL;
    } else {
      blood->data=NULL;
    }
  }
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int read_allogg_data(
  const char *inputfile,
  struct Blood *blood,
  int verbose
) {
  char line[512],*str;
  struct DataList *points;
  struct BDataA *bdata;
  int datenotfound=1, bgnotfound=1, timenotfound=1;
  FILE *fp;
  int n;
  //time_t t;
 
  if(verbose>0) printf("%s(%s)\n", __func__, inputfile);
 
  blood->bgdate[0] = '\0';
  blood->background = -1.0E20;
 
  fp=fopen(inputfile,"r");
  if(!fp) {
    fprintf(stderr,"Error: cannot open \"%s\".\n",inputfile);
    perror("fopen");
    return 1;
  }
  /*
  t=0; if(gm_timez) blood->date=*gmtime(&t); else blood->date=*localtime(&t);
  */
  memset(&blood->date,0,sizeof(struct tm));
  points=malloc(sizeof(struct DataList));
  points->prev=NULL; points->next=NULL;
  bdata=&points->allogg;
  memset(bdata,0,sizeof(struct BDataA));
  blood->data=points;
  /* Read data */
  while(fgets(line,sizeof(line),fp)) {
    str=line;
    while(isspace((unsigned char)*str)) str++;
    if(*str=='#') { /* Comment line; find date */
      if(datenotfound) datenotfound=find_date(str, blood, verbose-2);
      if(bgnotfound) bgnotfound=find_background(str, blood);
      continue;
    }
    if(*str=='\0') continue;
    n=sscanf(str,"%lf %d %d", &bdata->time_since_start,&bdata->data[0],&bdata->data[1]);
    if(timenotfound && strlen(str)>6&&strlen(str)<9) {
            blood->date.tm_hour=(str[0]-'0')*10+(str[1]-'0');
            blood->date.tm_min=(str[2]-'0')*10+(str[3]-'0');
            blood->date.tm_sec=(str[4]-'0')*10+(str[5]-'0');
            blood->date.tm_isdst=-1; timenotfound=0;
    } else if(n<2) {
      fprintf(stderr,"Error: malformed line %sCorrect syntax is hhmmss.\n",str);
      return 2;
    } else {
      points->next=malloc(sizeof(struct DataList));
      points->next->prev=points;
      points=points->next;
      points->next=NULL;
      bdata=&points->allogg;
      memset(bdata,0,sizeof(struct BDataA));
    }
  }
  fclose(fp);
  if(bdata->time_since_start==0) { /* Remove extra entry if any */
    struct DataList *prev=points->prev;
    free(points);
    points=prev;
    if(points) {
      points->next=NULL;
    } else {
      blood->data=NULL;
    }
  }
 
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int read_allogg2_data(
  char *inputfile,
  struct Blood *blood,
  int verbose,
  char *msg
) {
  IFT ift;
  int ii, n, ret;
  int singles, coincidents;
  char tmp[1024];
  double f, tas;
  //time_t t;
  struct DataList *points;
  struct BDataA *bdata;
 
 
  if(verbose>0) printf("%s(%s)\n", __func__, inputfile);
  /* Check the arguments */
  if(inputfile==NULL || blood==NULL || strlen(inputfile)<1) {
    if(msg!=NULL) strcpy(msg, "program error"); 
    return(1);
  }
  blood->bgdate[0] = '\0';
  blood->background = -1.0E20;
 
  /* Read Allogg2 file into IFT struct */
  iftInit(&ift); ret=iftRead(&ift, inputfile, 0, 0);
  if(ret) {if(msg!=NULL) strcpy(msg, ift.status); iftEmpty(&ift); return(2);}
  if(verbose>=12) {iftWrite(&ift, "stdout", 0); fflush(stdout);}
 
  /* Check that this is TPC Allogg2 */
  strcpy(tmp, "System ID"); ii=iftGet(&ift, tmp, 0);
  if(ii<0) {
    if(msg!=NULL) strcpy(msg, "wrong format"); 
    iftEmpty(&ift); return(3);
  }
  if(strcmp(ift.item[ii].value, "ABSS09282")!=0) {
    if(msg!=NULL) strcpy(msg, "wrong sampler ID");
    iftEmpty(&ift); return(4);
  }
 
  /* Read the background */
  strcpy(tmp, "//Average background counts"); ii=iftFindNthKey(&ift, tmp, 1, 0);
  if(ii>=0) {
    if(verbose>3) printf("key := '%s'\n", ift.item[ii].key);
    /* Get background date from the key */
    ret=isdate(ift.item[ii].key+28);
    if(ret!=0) ret=isdate2(ift.item[ii].key+28, tmp);
    if(ret!=0) ret=isdate3(ift.item[ii].key+28, tmp);
    if(ret==0) {
      if(verbose>1) printf("background_date := %s\n", tmp);
      strcpy(blood->bgdate, tmp);
    } else if(verbose>0) printf("background_date not valid\n");
    /* Get background counts/sec from value */
    if(verbose>1) printf("background_value := %s\n", ift.item[ii].value);
    blood->background=atof(ift.item[ii].value);
  } else if(verbose>0) printf("background not found\n");
 
  /* Read the isotope (half-life) */
  strcpy(tmp, "HalfTime"); ii=iftGet(&ift, tmp, 0);
  if(ii>=0) {
    f=atof(ift.item[ii].value);
    ii=hlIsotopeFromHalflife(f);
    /* Use it only if user did not specify the isotope */
    if(strcmp(isotope, "UNKNOWN")==0) {
      strcpy(isotope, hlIsotopeCode(ii));
      halflife=hlFromIsotope(isotope);
      if(verbose) printf("half-life := %g min\n", halflife);
    }
  }
 
  /* Read study number */
  strcpy(tmp, "Run number"); ii=iftGet(&ift, tmp, 0);
  if(ii>=0) {
#if(1)
    studynr_from_fname2(ift.item[ii].value, studynr, 0);
#else
    int len, i, j;
    strncpy(studynr, ift.item[ii].value, MAX_STUDYNR_LEN);
    /* Remove possible extra tail */
    len=strlen(studynr);
    for(i=0; i<len; i++) if(!isalnum((int)studynr[i])) {strcpy(studynr, ""); len=0;}
    for(i=0; i<len; i++) if(!isalpha((int)studynr[i])) break;
    if(i<1 || i>5) {strcpy(studynr, ""); len=0;}
    for(j=0; (i+j)<len; j++) if(!isdigit((int)studynr[i+j])) {studynr[i+j]=(char)0; break;}
    if(j<1 || j>5) {strcpy(studynr, ""); len=0;}
#endif
    if(verbose>0) printf("studynr := %s\n", studynr);
  }
 
 
  /* Find the first line which starts with date and time, that is, */
  /* the first count data line */
  for(ii=0; ii<ift.keyNr; ii++) if(!ift.item[ii].key[0]) {
    if(isdatetime(ift.item[ii].value, NULL)!=0) continue;
    if(verbose>7) printf("%s\n", ift.item[ii].value);
    break;
  } 
  /* Get start time from it */
  /*
  t=0; if(gm_timez) blood->date=*gmtime(&t); else blood->date=*localtime(&t);
  */
  //memset(&blood->date,0,sizeof(struct tm));
  //find_date(ift.item[ii].value, blood);
  ret=get_datetime(ift.item[ii].value, &(blood->date), verbose-4);
  if(ret!=0 && verbose>1) {
    fprintf(stderr, "Error %d in reading date and time in %s\n", ret, inputfile);
  }
 
  /* Read the count data lines */
  if(verbose>1) printf("reading count data\n");
  /* Initiate the list */
  points=malloc(sizeof(struct DataList));
  points->prev=NULL; points->next=NULL;
  bdata=&points->allogg;
  memset(bdata,0,sizeof(struct BDataA));
  blood->data=points;
  n=0;
  for(; ii<ift.keyNr; ii++) {
    /* Stop if item key is found */
    if(ift.item[ii].key[0]) break;
    /* Ignore comment lines */
    if(strncmp(ift.item[ii].value, "//", 2)==0) continue;
    if(strncmp(ift.item[ii].value, "#", 1)==0) continue;
    if(strncmp(ift.item[ii].value, ";", 1)==0) continue;
    /* Stop if line does not start with valid date and time */
    if(isdatetime(ift.item[ii].value, NULL)!=0) break;
    /* Read time after start, singles and coincidents */
    ret=sscanf(ift.item[ii].value+19, "%lf %d %d", &tas, &singles, &coincidents);
    if(ret!=3) continue;
    if(verbose>10 && n<4) printf("    %f  %d\n", tas, coincidents);
    /* Increase the list */
    if(n>0) {
      points->next=malloc(sizeof(struct DataList));
      points->next->prev=points;
      points=points->next;
      points->next=NULL;
      bdata=&points->allogg;
      memset(bdata,0,sizeof(struct BDataA));
    }
    /* Copy to list */
    bdata->time_since_start=tas;
    bdata->data[0]=coincidents;
    bdata->data[1]=0;
    n++;
  }
 
  iftEmpty(&ift);
  if(n<1) {if(msg!=NULL) strcpy(msg, "no valid count data"); return(9);}
  if(msg!=NULL) strcpy(msg, "ok");
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
void correct_scanditronics(
  struct Blood *blood,
  int verbose
) {
  struct DataList *data;
  struct tm date; //*date;
  time_t tim;
 
  if(verbose>0) printf("%s()\n", __func__);
 
  data=blood->data;
  tim=data->scanditronics.time_of_day;
  if(verbose>2) {
    char tmp[32];
    if(!ctime_r_int(&tim, tmp)) strcpy(tmp, "1900-01-01 00:00:00");
    printf("   %s\n", tmp);
  }
  if(verbose>3) {
    char buf[32];
    strftime(buf, 32, "%Y-%m-%d %H:%M:%S", &blood->date);
    printf("   blood->date: %s\n", buf);
  }
  /* Set time if it was not properly set in the comments */
  if(pump==2 || (blood->date.tm_hour==0 && blood->date.tm_min==0 && blood->date.tm_sec==0)) {
    //if(gm_timez) date=gmtime(&tim); else date=localtime(&tim);
    if(gmtime_r(&tim, &date)!=NULL) {
      blood->date.tm_hour=date.tm_hour;
      blood->date.tm_min=date.tm_min;
      blood->date.tm_sec=date.tm_sec;
    } else {
      fprintf(stderr, "Warning: invalid date in file.\n");
    }
  }
  while(data) {
    struct BDataS *d=&data->scanditronics;
    d->results.t=d->time_since_start+0.5*d->meas_interval;
    d->results.ch1=(double)d->pair1_coincid;
    d->results.ch2=(double)d->pair2_coincid;
    if(d->meas_interval>0) {
      d->results.ch1/=(double)d->meas_interval;
      d->results.ch2/=(double)d->meas_interval;
    }
    data=data->next;
  }
}
/*****************************************************************************/
 
/*****************************************************************************/
int correct_allogg(
  struct Blood *blood
) {
  struct DataList *data;
  double meas_interval=0.0;
 
  data=blood->data;
  if(data->next)
    meas_interval=data->next->allogg.time_since_start-data->allogg.time_since_start;
  else {
    printf("Error: Only one blood data line found, cannot calculate frame time.\n");
    return 1;
  }
  while(data) {
    struct BDataA *d=&data->allogg;
    d->meas_interval=meas_interval;
    d->results.t=d->time_since_start-0.5*meas_interval; /* '+'->'-' by VO */
#if 0
    d->results.ch1=(d->data[0]-d->data[1])/d->meas_interval;
    d->results.ch2=(d->data[0]-d->data[1])/d->meas_interval;
#else /* Use only second discriminator */
    if(meas_interval == 0) {
      printf("Error: zero measuring interval found, would lead to divide by zero\n");
      return 1;
    }
    d->results.ch1=d->data[1]/d->meas_interval;
    d->results.ch2=d->data[1]/d->meas_interval;
#endif
    data=data->next;
  }
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int correct_allogg2(
  struct Blood *blood,
  int verbose
) {
  if(verbose>0) printf("%s()\n", __func__);
 
  struct DataList *data;
  double meas_interval=0.0;
 
  data=blood->data;
  if(data->next) {
    meas_interval=data->next->allogg.time_since_start-data->allogg.time_since_start;
    if(verbose>1) {
      printf("time_since_start := %g\n", data->allogg.time_since_start);
      printf("time_since_start_next := %g\n", data->next->allogg.time_since_start);
      printf("data0 := %d\n", data->allogg.data[0]);
      printf("data0 next := %d\n", data->next->allogg.data[0]);
    }
  } else {
    printf("Error: Only one blood data line found, cannot calculate frame time.\n");
    return 1;
  }
  while(data) {
    struct BDataA *d=&data->allogg;
    d->meas_interval=meas_interval;
    if(meas_interval == 0) {
      printf("Error: zero measuring interval found, would lead to divide by zero\n");
      if(verbose>1) {
        printf("time_since_start := %g\n", data->allogg.time_since_start);
        printf("time_since_start_next := %g\n", data->next->allogg.time_since_start);
      }
      return 1;
    }
    d->results.t=d->time_since_start+0.5*meas_interval;
    d->results.ch1=d->data[0]/d->meas_interval;
    d->results.ch2=d->data[0]/d->meas_interval;
    data=data->next;
  }
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
void subtract_background(
  struct Blood *blood
) {
  struct DataList *data=blood->data;
  while(data) {
    struct BDataResults *res=&data->results;
 
    res->ch1 -= blood->background;
    res->ch2 -= blood->background;
    data=data->next;
  }
}
/*****************************************************************************/
 
/*****************************************************************************/
int get_calfactors(
  const char *file,
  time_t bloodt, 
  double *wcf, 
  double *pumpf,
  int pump, 
  char *isotope,
  char *caldate,
  int verbose
) {
  if(verbose>0) printf("%s(%s, ..., %d, %s, ...)\n", __func__, file, pump, isotope);
 
  const time_t too_old = 60*60*24*31*6;  /* 6 months */
  struct tm calibrdate; //*calibrdate;
  char line[512],*str;
  int i,yyyy,mm,dd;
  time_t calibrt, last_date=0;
  char sdate[25];
  FILE *fp;
 
  strcpy(sdate, "");
  if(pump<1) {
    fprintf(stderr,"Warning: pump number %d too low, set to 1\n",pump);
    pump=1;
  } else if(pump>MAX_PUMPS) {
    fprintf(stderr,"Warning: pump number %d too high, set to %d\n", pump, MAX_PUMPS);
    pump=MAX_PUMPS;
  }
 
  fp=fopen(file,"r");
  if(!fp) {
    fprintf(stderr,"Error: cannot open file %s\n",file);
    perror("fopen");
    return 1;
  }
  while(fgets(line,sizeof(line),fp)) {
    str=line;
    while(*str && isspace((unsigned char)*str)) str++;
    if(*str=='#' || *str=='\0') continue;
 
    double p[MAX_PUMPS]; /* Pumps */
    double w; /* Well-counter */
    i=sscanf(str,"%d-%d-%d %lf %lf %lf %lf %lf",
             &yyyy, &mm, &dd, &p[0], &p[1], &p[2], &p[3], &w);
    if(i<5) {
      printf("Error: Malformed line in calibration file, cannot continue\n");
      fclose(fp); return 1;
    }
    if(pump>i-4) {
      printf("Error: no calibration for this blood sampler in %s\n", file);
      fclose(fp); return 2;
    }
    if(i==5) {w=p[1];} else if(i==6) {w=p[2];} else if(i==7) {w=p[3];}
 
    calibrdate.tm_mday=dd; calibrdate.tm_mon=mm-1;
    calibrdate.tm_year=yyyy-1900;
    calibrdate.tm_hour=calibrdate.tm_min=calibrdate.tm_sec=0;
    calibrdate.tm_isdst=-1;
    calibrt=timegm(&calibrdate);
    if(calibrt==-1) calibrt=0;
    if(calibrt>=last_date) last_date=calibrt;
/*
    calibrt=time(NULL); calibrdate=localtime(&calibrt);
    calibrdate->tm_mday=dd; calibrdate->tm_mon=mm-1;
    calibrdate->tm_year=yyyy-1900;
    calibrdate->tm_isdst=-1;
    calibrt=mktime(calibrdate);
    if(calibrt>last_date) last_date=calibrt;
*/
 
    /* Check that calibration date is not later than measurement date */
    if(calibrt>bloodt) break;
    //if(calibrt!=(time_t)-1 && difftime(bloodt, calibrt)<0.0) break;
    strftime(sdate, 24, "%Y-%m-%d", &calibrdate);
 
    /* It was not, so set calibration factors and try the next line */
    *pumpf=p[pump-1];
    *wcf=w;
  }
 
  if(verbose>1) {
    printf("    pump_factor: %g\n", *pumpf);
    printf("    wc_factor: %g\n", *wcf);
  }
 
  /* Correct for branching ratio */
  if(strcasecmp(isotope, "C-11")==0) *wcf/=BRANCHING_C;
  else if(strcasecmp(isotope, "F-18")==0) *wcf/=BRANCHING_F; 
  else if(strcasecmp(isotope, "O-15")==0) *wcf/=BRANCHING_O; 
  else if(strcasecmp(isotope, "GE-68")==0) *wcf/=BRANCHING_Ge;
  else if(strcasecmp(isotope, "GA-68")==0) *wcf/=BRANCHING_Ga;
  else if(strcasecmp(isotope, "Cu-64")==0) *wcf/=BRANCHING_Cu64;
  else fprintf(stderr, "Warning: branching factor not included in calibration.\n");
 
  if(verbose>1) {
    printf("    wc_factor_including_branching_correction: %g\n", *wcf);
  }
 
  fclose(fp);
  /* Warn user if latest calibration date is older than 6 months */
  if(bloodt>last_date && bloodt-last_date > too_old) {
    fprintf(stderr, "Warning: check that %s contains the latest calibration coefficients.\n", file);
  }
 
  /* Set return value for calibration date */
  if(caldate!=NULL) strcpy(caldate, sdate);
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/