simdispersion.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 "tpccm.h"
/*****************************************************************************/
 
/*****************************************************************************/
int simDispersion(
  double *x,
  double *y,
  const int n, 
  const double tau1, 
  const double tau2,
  double *tmp
) {
  /* Check input */
  if(x==NULL || y==NULL || n<2) return 1;
  if(tau1<0.0 || tau2<0.0) return 2;
 
  /* Allocate memory if not allocated by user */
  double *buf;
  if(tmp!=NULL) buf=tmp; else buf=(double*)malloc(n*sizeof(double));
  if(buf==NULL) return 3;
 
  /* First dispersion */
  if(tau1>0.0) {
    double k=1.0/tau1;
    int ret=simC1(x, y, n, k, k, buf); 
    if(ret!=0) {
      if(tmp==NULL) free(buf);
      return 100+ret;
    }
    for(int i=0; i<n; i++) y[i]=buf[i];
  }
 
  /* Second dispersion */
  if(tau2>0.0) {
    double k=1.0/tau2;
    int ret=simC1(x, y, n, k, k, buf);
    if(ret!=0) {
      if(tmp==NULL) free(buf);
      return 200+ret;
    }
    for(int i=0; i<n; i++) y[i]=buf[i];
  }
 
  if(tmp==NULL) free(buf);
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int corDispersion(
  double *x,
  double *y,
  const int n, 
  const double tau,
  double *tmp
) {
  /* Check input */
  if(x==NULL || y==NULL || n<2) return 1;
  if(tau<0.0) return 2;
  if(x[0]<0.0) return 3;
 
  /* Allocate memory if not allocated by user */
  double *buf;
  if(tmp!=NULL) buf=tmp; else buf=(double*)malloc(n*sizeof(double));
  if(buf==NULL) return 3;
 
  /* Integrate measured data */
  buf[0]=0.5*y[0]*x[0];
  for(int i=1; i<n; i++) buf[i]=buf[i-1]+0.5*(y[i]+y[i-1])*(x[i]-x[i-1]);
 
  /* Calculate the integral of dispersion corrected data */
  for(int i=0; i<n; i++) buf[i]+=tau*y[i];
 
  /* Calculate dispersion corrected curve from its integral */
  {
    int i=0;
    double dt=x[i];
    if(dt>0.0) y[i]=2.0*buf[i]/dt; else y[i]=0.0;
    for(i=1; i<n-1; i++) {
      dt=x[i+1]-x[i-1]; if(dt>0.0) y[i]=(buf[i+1]-buf[i-1])/dt; else y[i]=y[i-1];
    }
    dt=x[i]-x[i-1]; if(dt>0.0) y[i]=(buf[i]-buf[i-1])/dt; else y[i]=y[i-1];
  }
 
  if(tmp==NULL) free(buf);
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int simTTM(
  double *t,
  double *c0,
  const int n,
  const double k,
  const int cn,
  double *cout
) {
  /* Check for data */
  if(n<2 || cn<1) return(1);
  if(t==NULL || c0==NULL || cout==NULL) return(2);
  /* Check the rate constant */
  if(!(k>=0.0)) return(3);
 
  double c[cn+1], ci[cn+1]; // Compartmental concentrations for current sample
  double c_last[cn+1], ci_last[cn+1]; // Compartmental concentrations for previous sample
  for(int j=0; j<=cn; j++) c[j]=c_last[j]=ci[j]=ci_last[j]=0.0;
 
  /* Calculate curves */
#if(0)
  printf("\nTime\tC0\tiC0");
  for(int j=1; j<=cn; j++) printf("\tC%d", j);
  printf("\n");
#endif
  double t_last=0.0; if(t[0]<t_last) t_last=t[0];
  double c0_last=0.0, c0i=0.0;
  for(int i=0; i<n; i++) { // loop through sample times
    /* delta time / 2 */
    double dt2=0.5*(t[i]-t_last); if(!(dt2>=0.0)) return(5);
    /* input integral */
    c0i+=(c0[i]+c0_last)*dt2;
    /* k/(1+k*(dt/2)) */
    double kdt=k/(1.0+dt2*k);
    /* calculate compartmental concentrations */
    ci[0]=c0i; 
    if(dt2>0.0) {
      for(int j=1; j<=cn; j++) {
        c[j] = kdt * (ci[j-1] - ci_last[j] - dt2*c_last[j]);
        ci[j] = ci_last[j] + dt2*(c_last[j]+c[j]);
      }
    } else { // sample time same as previously, thus concentration do not change either
      for(int j=0; j<=cn; j++) {
        c[j]=c_last[j];
        ci[j]=ci_last[j];
      }
    }
#if(0)
    printf("%g\t%g\t%g", t[i], c0[i], ci[0]);
    for(int j=1; j<=cn; j++) printf("\t%g", c[j]);
    printf("\n");
#endif
    cout[i]=c[cn];
    /* prepare to the next loop */
    t_last=t[i]; c0_last=c0[i];
    for(int j=0; j<=cn; j++) {
      c_last[j]=c[j];
      ci_last[j]=ci[j];
    }
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int simTTM_i(
  double *t,
  double *c0i,
  const int n,
  const double k,
  const int cn,
  double *cout
) {
  /* Check for data */
  if(n<2 || cn<1) return(1);
  if(t==NULL || c0i==NULL || cout==NULL) return(2);
  /* Check the rate constant */
  if(!(k>=0.0)) return(3);
 
  double c[cn+1], ci[cn+1]; // Compartmental concentrations for current sample
  double c_last[cn+1], ci_last[cn+1]; // Compartmental concentrations for previous sample
  for(int j=0; j<=cn; j++) c[j]=c_last[j]=ci[j]=ci_last[j]=0.0;
 
  /* Calculate curves */
#if(0)
  printf("\nTime\tiC0");
  for(int j=1; j<=cn; j++) printf("\tC%d", j);
  printf("\n");
#endif
  double t_last=0.0; if(t[0]<t_last) t_last=t[0];
  for(int i=0; i<n; i++) { // loop through sample times
    /* delta time / 2 */
    double dt2=0.5*(t[i]-t_last); if(!(dt2>=0.0)) return(5);
    /* k/(1+k*(dt/2)) */
    double kdt=k/(1.0+dt2*k);
    /* calculate compartmental concentrations */
    ci[0]=c0i[i]; 
    if(dt2>0.0) {
      for(int j=1; j<=cn; j++) {
        c[j] = kdt * (ci[j-1] - ci_last[j] - dt2*c_last[j]);
        ci[j] = ci_last[j] + dt2*(c_last[j]+c[j]);
      }
    } else { // sample time same as previously, thus concentration do not change either
      for(int j=0; j<=cn; j++) {
        c[j]=c_last[j];
        ci[j]=ci_last[j];
      }
    }
#if(0)
    printf("%g\t%g", t[i], ci[0]);
    for(int j=1; j<=cn; j++) printf("\t%g", c[j]);
    printf("\n");
#endif
    cout[i]=c[cn];
    /* prepare to the next loop */
    t_last=t[i];
    for(int j=0; j<=cn; j++) {
      c_last[j]=c[j];
      ci_last[j]=ci[j];
    }
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/