sim3cms.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 simC3s(
  double *t,
  double *ca,
  const int nr,
  double k1,
  double k2,
  double k3,
  double k4,
  double k5,
  double k6,
  double *ct,
  double *cta,
  double *ctb,
  double *ctc
) {
  int i;
  double b, c, d, w, z, dt2;
  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 actual parameter number */
  if(k1<0.0) return 3;
  if(k3<=0.0) {k3=0.0; if(k2<=0.0) k2=0.0;}
  else if(k5<=0.0) {k5=0.0; if(k4<=0.0) k4=0.0;}
  else {if(k6<=0.0) k6=0.0;}
 
  /* 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;
      /* partial results */
      b=ct1i_last+dt2*ct1_last;
      c=ct2i_last+dt2*ct2_last;
      d=ct3i_last+dt2*ct3_last;
      w=k4 + k5 - (k5*k6*dt2)/(1.0+k6*dt2);
      z=1.0+w*dt2;
      /* 1st tissue compartment and its integral */
      ct1 = (
          + k1*z*cai + (k3*k4*dt2 - (k2+k3)*z)*b
          + k4*c + k4*k6*dt2*d/(1.0+k6*dt2)
        ) / ( z*(1.0 + dt2*(k2+k3)) - k3*k4*dt2*dt2 );
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
      /* 2nd tissue compartment and its integral */
      ct2 = (k3*ct1i - w*c + k6*d/(1.0+k6*dt2)) / z;
      ct2i = ct2i_last + dt2*(ct2_last+ct2);
      /* 3rd tissue compartment and its integral */
      ct3 = (k5*ct2i - k6*d) / (1.0 + k6*dt2);
      ct3i = ct3i_last + dt2*(ct3_last+ct3);
    }
    /* copy values to argument arrays */
    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 simC3vs(
  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 b, c, d, w, z, dt2, 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;
      /* partial results */
      b=ct1i_last+dt2*ct1_last;
      c=ct2i_last+dt2*ct2_last;
      d=ct3i_last+dt2*ct3_last;
      w=k4 + k5 - (k5*k6*dt2)/(1.0+k6*dt2);
      z=1.0+w*dt2;
      /* 1st tissue compartment and its integral */
      ct1 = (
          + k1*z*cai + (k3*k4*dt2 - (k2+k3)*z)*b
          + k4*c + k4*k6*dt2*d/(1.0+k6*dt2)
        ) / ( z*(1.0 + dt2*(k2+k3)) - k3*k4*dt2*dt2 );
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
      /* 2nd tissue compartment and its integral */
      ct2 = (k3*ct1i - w*c + k6*d/(1.0+k6*dt2)) / z;
      ct2i = ct2i_last + dt2*(ct2_last+ct2);
      /* 3rd tissue compartment and its integral */
      ct3 = (k5*ct2i - k6*d) / (1.0 + k6*dt2);
      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;
}
/*****************************************************************************/
 
/*****************************************************************************/