mtga.c

Go to the documentation of this file.

/*****************************************************************************/
#include "libtpcmodel.h"
/*****************************************************************************/
 
/*****************************************************************************/
int patlak_data(
  int data_nr,
  double *i,
  double *ii,
  double *c,
  double *x,
  double *y
) {
  int fi, plot_nr=0;
  double divider_limit=1.0E-12;
 
  if(data_nr<0 || i==NULL || ii==NULL || c==NULL || x==NULL || y==NULL)
    return -1;
  for(fi=0; fi<data_nr; fi++) {
    // check that all measured data is available
    if(isnan(i[fi]) || isnan(ii[fi]) || isnan(c[fi])) continue;
    if(!(i[fi]>-1.0E+20 && i[fi]<+1.0E+20)) continue; 
    if(!(ii[fi]>-1.0E+20 && ii[fi]<+1.0E+20)) continue; 
    if(!(c[fi]>-1.0E+20 && c[fi]<+1.0E+20)) continue; 
    // check that integral has been >=0 all the time
    if(ii[fi]<0.0) {plot_nr=0; continue;}
    // check that dividers are not too close to zero
    if(fabs(i[fi])<divider_limit) continue;
    // calculate plot x axis value
    x[plot_nr]=ii[fi]/i[fi];
    if(!(x[plot_nr]>-1.0E+20 && x[plot_nr]<+1.0E+20)) continue;
    if(x[plot_nr]<0.0) continue;
    // calculate plot y axis value
    y[plot_nr]=c[fi]/i[fi];
    if(!(y[plot_nr]>-1.0E+20 && y[plot_nr]<+1.0E+20)) continue;
    // so this plot data point is fine
    plot_nr++;
  }
  return(plot_nr);
}
/*****************************************************************************/
 
/*****************************************************************************/
int logan_data(
  int data_nr,
  double *i,
  double *ii,
  double *c,
  double *ci,
  double k2,
  double *x,
  double *y
) {
  int fi, plot_nr=0;
  double divider_limit=1.0E-18;
 
  if(data_nr<0 || i==NULL || ii==NULL || c==NULL || ci==NULL) return -1;
  if(x==NULL || y==NULL) return -1;
 
  for(fi=0; fi<data_nr; fi++) {
    // check that all measured data is available
    if(isnan(i[fi]) || isnan(ii[fi]) || isnan(c[fi]) || isnan(ci[fi])) continue;
    if(!(i[fi]>-1.0E+30 && i[fi]<+1.0E+30)) continue; 
    if(!(ii[fi]>-1.0E+30 && ii[fi]<+1.0E+30)) continue; 
    if(!(c[fi]>-1.0E+30 && c[fi]<+1.0E+30)) continue; 
    if(!(ci[fi]>-1.0E+30 && ci[fi]<+1.0E+30)) continue; 
    // check that integrals have been >=0 all the time
    if(ii[fi]<0.0) {plot_nr=0; continue;}
    if(ci[fi]<0.0) {plot_nr=0; continue;}
    // check that dividers are not too close to zero
    if(fabs(c[fi])<divider_limit) continue;
    // calculate plot x axis value
    if(k2>0.0) x[plot_nr]=(ii[fi]+i[fi]/k2)/c[fi];
    else x[plot_nr]=ii[fi]/c[fi];
    if(!(x[plot_nr]>-1.0E+30 && x[plot_nr]<+1.0E+30)) continue;
    // calculate plot y axis value
    y[plot_nr]=ci[fi]/c[fi];
    if(!(y[plot_nr]>-1.0E+30 && y[plot_nr]<+1.0E+30)) continue;
    // so this plot data point is fine
    plot_nr++;
  }
  return(plot_nr);
}
/*****************************************************************************/
 
/*****************************************************************************/
int mtga_best_perp(
  double *x,
  double *y,
  int nr,
  double *slope,
  double *ic,
  double *ssd,
  int *fnr
) {
  int from, to, ret, from_min, to_min;
  double lic, lslope, lssd, ssd_min;
 
  /* Search the plot range that gives the lowest ssd */
  ssd_min=9.99E+99; from_min=to_min=-1;
  for(from=0, to=nr-1; ((to-from)+1)>=MTGA_BEST_MIN_NR; from++) {
    ret=llsqperp(x+from, y+from, (to-from)+1, &lslope, &lic, &lssd);
    if(ret==0 && lssd<ssd_min) {
      ssd_min=lssd; from_min=from; to_min=to;
      *slope=lslope; *ic=lic; if(ssd!=NULL) *ssd=lssd;
    }
  }
  if(from_min<0) {
    if(fnr!=NULL) *fnr=0;
    if(ssd!=NULL) *ssd=0.0;
    return(5);
  }
  if(fnr!=NULL) *fnr=(to_min-from_min)+1;
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/