sim3cmp.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 simC3p(
  double *t,
  double *ca,
  const int nr,
  const double k1,
  const double k2,
  const double k3,
  const double k4,
  const double k5,
  const double k6,
  double *ct,
  double *cta,
  double *ctb,
  double *ctc
) {
  int i;
  double dt2, r, s, u, v, w;
  double cai, ca_last, t_last;
  double ct1, ct1_last, ct2, ct2_last, ct3, ct3_last;
  double ct1i, ct1i_last, ct2i, ct2i_last, ct3i, ct3i_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(ct==NULL) return 2;
 
  /* Check parameters */
  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=ct2_last=ct3_last=ct1i_last=ct2i_last=ct3i_last=0.0;
  ct1=ct2=ct3=ct1i=ct2i=ct3i=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;
      /* Calculate partial results */
      r=1.0+k4*dt2;
      s=1.0+k6*dt2;
      u=ct1i_last+dt2*ct1_last;
      v=ct2i_last+dt2*ct2_last;
      w=ct3i_last+dt2*ct3_last;
      /* 1st tissue compartment and its integral */
      ct1 = ( k1*cai - (k2 + (k3/r) + (k5/s))*u + (k4/r)*v + (k6/s)*w )
            / ( 1.0 + dt2*(k2 + (k3/r) + (k5/s)) );
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
      /* 2nd tissue compartment and its integral */
      ct2 = (k3*ct1i - k4*v) / r;
      ct2i = ct2i_last + dt2*(ct2_last+ct2);
      /* 3rd tissue compartment and its integral */
      ct3 = (k5*ct1i - k6*w) / s;
      ct3i = ct3i_last + dt2*(ct3_last+ct3);
    }
    /* copy values to argument arrays; set very small values to zero */
    ct[i]=ct1+ct2+ct3;
    if(cta!=NULL) cta[i]=ct1;
    if(ctb!=NULL) ctb[i]=ct2;
    if(ctc!=NULL) ctc[i]=ct3;
    /* prepare to the next loop */
    t_last=t[i]; ca_last=ca[i];
    ct1_last=ct1; ct1i_last=ct1i;
    ct2_last=ct2; ct2i_last=ct2i;
    ct3_last=ct3; ct3i_last=ct3i;
  }
 
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int simC3vp(
  double *t,
  double *ca,
  double *cb,
  const int nr,
  const double k1,
  const double k2,
  const double k3,
  const double k4,
  const double k5,
  const double k6,
  const double f,
  const double vb,
  const double fa,
  const int vvm,
  double *cpet,
  double *cta,
  double *ctb,
  double *ctc,
  double *ctab,
  double *ctvb
) {
  int i;
  double dt2, r, s, u, v, w, va, vv;
  double cai, ca_last, t_last, dct, cvb;
  double ct1, ct1_last, ct2, ct2_last, ct3, ct3_last;
  double ct1i, ct1i_last, ct2i, ct2i_last, ct3i, ct3i_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(cpet==NULL) return 2;
 
  /* Check parameters */
  if(k1<0.0) return 3;
  if(vb<0.0 || vb>=1.0) return 4;
  if(fa<=0.0 || fa>1.0) return 5;
  va=fa*vb; vv=(1.0-fa)*vb;
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0];
  cai=ca_last=0.0;
  ct1_last=ct2_last=ct3_last=ct1i_last=ct2i_last=ct3i_last=0.0;
  ct1=ct2=ct3=ct1i=ct2i=ct3i=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;
      /* Calculate partial results */
      r=1.0+k4*dt2;
      s=1.0+k6*dt2;
      u=ct1i_last+dt2*ct1_last;
      v=ct2i_last+dt2*ct2_last;
      w=ct3i_last+dt2*ct3_last;
      /* 1st tissue compartment and its integral */
      ct1 = ( k1*cai - (k2 + (k3/r) + (k5/s))*u + (k4/r)*v + (k6/s)*w )
            / ( 1.0 + dt2*(k2 + (k3/r) + (k5/s)) );
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
      /* 2nd tissue compartment and its integral */
      ct2 = (k3*ct1i - k4*v) / r;
      ct2i = ct2i_last + dt2*(ct2_last+ct2);
      /* 3rd tissue compartment and its integral */
      ct3 = (k5*ct1i - k6*w) / s;
      ct3i = ct3i_last + dt2*(ct3_last+ct3);
    }
    /* Venous curve */
    if(f>0.) {dct = k1*ca[i] - k2*ct1; cvb = cb[i] - dct/f;} else cvb=cb[i];
    /* copy values to argument arrays */
    if(vvm==0) cpet[i]= va*cb[i] + vv*cvb + (1.0-vb)*(ct1+ct2+ct3);
    else cpet[i]= va*cb[i] + vv*cvb + (ct1+ct2+ct3);
    if(cta!=NULL) {cta[i]=ct1; if(vvm==0) cta[i]*=(1.0-vb);}
    if(ctb!=NULL) {ctb[i]=ct2; if(vvm==0) ctb[i]*=(1.0-vb);}
    if(ctc!=NULL) {ctc[i]=ct3; if(vvm==0) ctc[i]*=(1.0-vb);}
    if(ctab!=NULL) {ctab[i]=va*cb[i];}
    if(ctvb!=NULL) {ctvb[i]=vv*cvb;}
    /* prepare to the next loop */
    t_last=t[i]; ca_last=ca[i];
    ct1_last=ct1; ct1i_last=ct1i;
    ct2_last=ct2; ct2i_last=ct2i;
    ct3_last=ct3; ct3i_last=ct3i;
  }
 
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/