sim1cm.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 simMBF(
  double *t,
  double *ci,
  const int nr,
  const double k1,
  const double k2,
  const double Vfit,
  double *ct
) {
  int i;
  double dt2;
  double cii, ci_last, t_last;
  double ct_last, cti, cti_last;
 
 
  /* Check for data */
  if(nr<2) return(1);
  if(t==NULL || ci==NULL || ct==NULL) return(2);
 
  /* Calculate curves */
  t_last=0.0; cii=ci_last=0.0;
  cti=ct_last=cti_last=0.0;
  for(i=0; i<nr; i++) {
    /* delta time / 2 */
    dt2=0.5*(t[i]-t_last);
    /* calculate values */
    if(dt2<0.0) {
      return(5);
    } else if(dt2>0.0) {
      /* input integral */
      cii+=(ci[i]+ci_last)*dt2;
      /* Tissue compartment and its integral */
      ct[i] = (Vfit*ci[i] + k1*cii - k2*(cti_last+dt2*ct_last)) / (1.0 + dt2*k2);
      cti = cti_last + dt2*(ct_last+ct[i]);
    }
    /* set very small values to zero
    if(fabs(ct[i])<1.0e-12) ct[i]=0.0; */
    /* prepare to the next loop */
    t_last=t[i]; ci_last=ci[i];
    ct_last=ct[i]; cti_last=cti;
  }
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int simC1(
  double *t,
  double *ca,
  const int nr,
  const double k1,
  const double k2,
  double *ct
) {
  int i;
  double dt2;
  double cai, ca_last, t_last;
  double ct1, ct1_last;
  double ct1i, ct1i_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(t==NULL || ca==NULL || ct==NULL) return 2;
 
  /* Check actual parameter number */
  if(!(k1>=0.0)) return 3;
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0]; 
  cai=ca_last=0.0;
  ct1_last=ct1i_last=0.0;
  ct1=ct1i=0.0;
  for(i=0; i<nr; i++) {
    /* delta time / 2 */
    dt2=0.5*(t[i]-t_last);
    /* calculate values */
    if(dt2<0.0) {
      return 5;
    } else if(dt2>0.0) {
      /* arterial integral */
      cai+=(ca[i]+ca_last)*dt2;
      /* tissue compartment and its integral */
      ct1 = (k1*cai - k2*(ct1i_last+dt2*ct1_last)) / (1.0 + dt2*k2);
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
    }
    /* copy values to argument arrays */
    ct[i]=ct1;
    /* prepare to the next loop */
    t_last=t[i]; ca_last=ca[i];
    ct1_last=ct1; ct1i_last=ct1i;
  }
 
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int simC1_i(
  double *t,
  double *cai,
  const int nr,
  const double k1,
  const double k2,
  double *ct
) {
  int i;
  double dt2;
  double t_last;
  double ct1, ct1_last;
  double ct1i, ct1i_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(t==NULL || cai==NULL || ct==NULL) return 2;
 
  /* Check actual parameter number */
  if(!(k1>=0.0)) return 3;
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0]; 
  ct1_last=ct1i_last=0.0;
  ct1=ct1i=0.0;
  for(i=0; i<nr; i++) {
    /* delta time / 2 */
    dt2=0.5*(t[i]-t_last);
    /* calculate values */
    if(dt2<0.0) {
      return 5;
    } else if(dt2>0.0) {
      /* tissue compartment and its integral */
      ct1 = (k1*cai[i] - k2*(ct1i_last+dt2*ct1_last)) / (1.0 + dt2*k2);
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
    }
    /* copy values to argument arrays */
    ct[i]=ct1;
    /* prepare to the next loop */
    t_last=t[i]; ct1_last=ct1; ct1i_last=ct1i;
  }
 
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int simC1_d(
  double *t,
  double *ca,
  const int nr,
  const double k1,
  const double k2,
  double *ct
) {
  int i;
  double dt2;
  double cap, ctp, tp; // 'p' marks previous
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(t==NULL || ca==NULL || ct==NULL) return 2;
 
  /* Check actual parameter number */
  if(!(k1>=0.0)) return 3;
 
  /* Set 'previous' sample */
  tp=0.0; if(t[0]<tp) tp=t[0]; 
  cap=ctp=0.0;
  /* Calculate curves */
  for(i=0; i<nr; i++) {
    /* delta time / 2 */
    dt2=0.5*(t[i]-tp);
    /* calculate values */
    if(dt2<0.0) {
      return 5;
    } else if(dt2>=0.0) { // must be >= because ct[0] is not set elsewhere
      ct[i]=(k1*dt2*(cap+ca[i]) + (1.0 - k2*dt2)*ctp) / (1.0 + k2*dt2);
    }
    /* prepare to the next sample */
    tp=t[i]; cap=ca[i];
    ctp=ct[i];
  }
 
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/