simrtcm.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 simRTCM(
  double *t,
  double *cr,
  const int nr,
  const double R1,
  const double k2,
  const double k3,
  const double k4,
  double *ct,
  double *cta,
  double *ctb
) {
  int i;
  double f, b, w, dt2;
  double cri, cr_last, t_last;
  double cf, cf_last, cb, cb_last;
  double cfi, cfi_last, cbi, cbi_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(ct==NULL) return 2;
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0];
  cri=cr_last=0.0; cf_last=cb_last=cfi_last=cbi_last=cf=cb=cfi=cbi=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) {
      /* reference integral */
      cri+=(cr[i]+cr_last)*dt2;
      /* partial results */
      f=cfi_last+dt2*cf_last;
      b=cbi_last+dt2*cb_last;
      w=k2 + k3 + k2*k4*dt2;
      /* 1st tissue compartment and its integral */
      cf = ( (1.0 + k4*dt2)*(R1*cr[i] + k2*cri) + k4*b - w*f ) / ( 1.0 + dt2*(w+k4) );
      cfi = cfi_last + dt2*(cf_last+cf);
      /* 2nd tissue compartment and its integral */
      cb = (k3*cfi - k4*b) / (1.0 + k4*dt2);
      cbi = cbi_last + dt2*(cb_last+cb);
    }
    /* copy values to argument arrays */
    ct[i]=cf+cb;
    if(cta!=NULL) cta[i]=cf;
    if(ctb!=NULL) ctb[i]=cb;
    /* prepare to the next loop */
    t_last=t[i]; cr_last=cr[i];
    cf_last=cf; cfi_last=cfi;
    cb_last=cb; cbi_last=cbi;
  }
 
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int simSRTM(
  double *t,
  double *cr,
  const int nr,
  const double R1,
  const double k2,
  const double BP,
  double *ct
) {
  int i;
  double dt2;
  double cri, cr_last, t_last;
  double ct_last, cti, cti_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(ct==NULL) return 2;
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0];
  cri=cr_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) {
      /* reference integral */
      cri+=(cr[i]+cr_last)*dt2;
      /* Tissue compartment and its integral */
      ct[i] = ( R1*cr[i] + k2*cri - (k2/(1.0+BP))*(cti_last+dt2*ct_last) ) /
              ( 1.0 + dt2*(k2/(1.0+BP)) );
      cti = cti_last + dt2*(ct_last+ct[i]);
    } else { // dt==0
      ct[i]=ct_last;
    }
    /* prepare to the next loop */
    t_last=t[i]; cr_last=cr[i];
    ct_last=ct[i]; cti_last=cti;
  }
 
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int simTRTM(
  double *t,
  double *cr,
  const int nr,
  const double R1,
  const double k2,
  const double k3,
  double *ct
) {
  int i;
  double dt2;
  double cri, cr_last, t_last;
  double ct_last, cti, cti_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(ct==NULL) return 2;
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0];
  cri=cr_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) {
      /* reference integral */
      cri+=(cr[i]+cr_last)*dt2;
      /* Tissue compartment and its integral */
      ct[i] = ( R1*cr[i] + R1*k3*cri - (k2+k3)*(cti_last+dt2*ct_last) ) / ( 1.0 + dt2*(k2+k3) );
      cti = cti_last + dt2*(ct_last+ct[i]);
    } else { // dt==0
      ct[i]=ct_last;
    }
    /* prepare to the next loop */
    t_last=t[i]; cr_last=cr[i];
    ct_last=ct[i]; cti_last=cti;
  }
 
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/