taccross.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <math.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpcift.h"
#include "tpctac.h"
#include "tpcli.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Determines the time when TAC in file 2 crosses over the TAC in file 1.",
  "The first time point (x) of file 2 after which at least the given nr of",
  "samples (y values) become all lower or higher than the samples in file 1",
  "is written to stdout.",
  "Only the first TAC of data files is used.",
  " ",
  "Usage: @P [options] file1 file2 nr",
  " ",
  "Options:",
  " -ift",
  "     Cross-over point is written in stdout with explanations.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "     @P iea345_lvcav.bld iea345_vena_cava.bld 3",
  " ",
  "See also: taccut, tacsplit, tacadd, taccut, tacdel, tacsort",
  " ",
  "Keywords: TAC, input, tool",
  0};
/*****************************************************************************/
 
/*****************************************************************************/
/* 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;
  int   ret;
  char *cptr, tacfile1[FILENAME_MAX], tacfile2[FILENAME_MAX];
  TAC   tac1, tac2;
  int   crossLimit=0;
  int   outputType=0; 
 
  
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacInit(&tac1); tacInit(&tac2);
  tacfile1[0]=tacfile2[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(!*cptr) continue;
    if(strcasecmp(cptr, "IFT")==0) {
      outputType=1; continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break; // tac name argument may start with '-'
  
  /* 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);}
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("outputType := %d\n", outputType);
  }
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-5;
 
  /* Arguments */
  if(ai<argc) strlcpy(tacfile1, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(tacfile2, argv[ai++], FILENAME_MAX);
  if(ai<argc) crossLimit=atoi(argv[ai++]);
  if(ai<argc) {fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]); return(1);}
 
  /* Is something missing? */
  if(!tacfile2[0] || crossLimit<1) {tpcPrintUsage(argv[0], info, stdout); return(1);}
  if(crossLimit<1) {fprintf(stderr, "Error: missing or invalid cross number.\n"); return(1);}
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    for(ai=0; ai<argc; ai++) printf("%s ", argv[ai]); 
    printf("\n");
    printf("tacfile1 := %s\n", tacfile1);
    printf("tacfile2 := %s\n", tacfile2);
    printf("crossLimit := %d\n", crossLimit);
    fflush(stdout);
  }
 
 
  /*
   *  Read the files
   */
  if(verbose>1) printf("reading %s\n", tacfile1);
  ret=tacRead(&tac1, tacfile1, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
    tacFree(&tac1); return(2);
  }
  if(verbose>3) {
    printf("fileformat1 := %s\n", tacFormattxt(tac1.format));
    printf("tacNr1 := %d\n", tac1.tacNr);
    printf("sampleNr1 := %d\n", tac1.sampleNr);
    printf("xunit1 := %s\n", unitName(tac1.tunit));
    printf("yunit1 := %s\n", unitName(tac1.cunit));
  }
  /* Use only the first TAC */
  if(tac1.tacNr>1) tac1.tacNr=1;
 
  if(verbose>1) printf("reading %s\n", tacfile2);
  ret=tacRead(&tac2, tacfile2, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
    tacFree(&tac1); tacFree(&tac2); return(2);
  }
  if(verbose>3) {
    printf("fileformat2 := %s\n", tacFormattxt(tac2.format));
    printf("tacNr2 := %d\n", tac2.tacNr);
    printf("sampleNr2 := %d\n", tac2.sampleNr);
    printf("xunit2 := %s\n", unitName(tac2.tunit));
    printf("yunit2 := %s\n", unitName(tac2.cunit));
  }
  /* Use only the first TAC */
  if(tac2.tacNr>1) tac2.tacNr=1;
 
 
  /* Convert TAC 1 into units of TAC 2 */
  ret=tacXUnitConvert(&tac1, tac2.tunit, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Warning: could not convert times from '%s' to '%s'.\n",
            unitName(tac1.tunit), unitName(tac2.tunit));
  }
  ret=tacYUnitConvert(&tac1, tac2.cunit, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Warning: could not convert concentrations from '%s' to '%s'.\n",
            unitName(tac1.cunit), unitName(tac2.cunit));
  }
 
  /* Interpolate TAC 1 to sample times of TAC 2 */
  tacSortByTime(&tac1, &status);
  tacSortByTime(&tac2, &status);
  ret=tacAllocateMore(&tac2, 1);
  if(ret==TPCERROR_OK)
    ret=liInterpolate(tac1.x, tac1.c[0].y, tac1.sampleNr, 
                      tac2.x, tac2.c[1].y, NULL, NULL, tac2.sampleNr, 
                      4, 1, verbose-8);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot interpolate data.\n");
    tacFree(&tac1); tacFree(&tac2); return(3);
  }
  tacFree(&tac1);
 
  /*
   *  Search for the crossing
   */
  int i, consecNr=0, crossNr=0, above=0;
  i=0; if(tac2.c[0].y[i]>tac2.c[1].y[i]) above=1;
  if(verbose>2) printf("  %d: %d   %g %g\n", 1+i, above, tac2.c[0].y[i], tac2.c[1].y[i]);
  for(i=1; i<tac2.sampleNr; i++) {
    if(tac2.c[0].y[i]>tac2.c[1].y[i]) { // above
      if(above) {consecNr++;} else {consecNr=0; crossNr++;}
      above=1;
    } else { // not above
      if(!above) {consecNr++;} else {consecNr=0; crossNr++;}
      above=0;
    }
    if(verbose>2) printf("  %d: %d   %g %g\n", 1+i, above, tac2.c[0].y[i], tac2.c[1].y[i]);
    if(crossNr>0 && consecNr>=crossLimit) break;
  }
  if(crossNr==0 || consecNr<crossLimit) {
    fprintf(stderr, "Error: cannot determine cross-over time.\n");
    tacFree(&tac2); return(4);
  }
  i-=consecNr;
  if(verbose>2) printf("cross-over_point := %d\n", 1+i);
  if(outputType==0) {
    printf("%g\n", tac2.x[i]);
  } else {
    printf("x := %g\n", tac2.x[i]);
    printf("y := %g\n", tac2.c[0].y[i]);
    printf("index := %d\n", i);
  }
  tacFree(&tac2);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/