simdispersion.c File Reference

Simulation of dispersion. More...

#include "tpcclibConfig.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include "tpccm.h"

Go to the source code of this file.

Functions
int	simDispersion (double *x, double *y, const int n, const double tau1, const double tau2, double *tmp)
int	corDispersion (double *x, double *y, const int n, const double tau, double *tmp)
int	simTTM (double *t, double *c0, const int n, const double k, const int cn, double *cout)
int	simTTM_i (double *t, double *c0i, const int n, const double k, const int cn, double *cout)

Detailed Description

Simulation of dispersion.

Definition in file simdispersion.c.

Function Documentation

corDispersion()

int corDispersion	(	double *	x,
		double *	y,
		const int	n,
		const double	tau,
		double *	tmp )

Correct time-activity curve for the effect of dispersion.

Data must be noise-free and have very short sampling intervals. The units of rate constants must be related to the TAC time units; 1/min and min, or 1/sec and sec.

Returns

Function returns 0 when successful, else a value >= 1.

Author

Vesa Oikonen

See also

simDispersion, simDelay

Parameters

x	Array of sample times; must be in increasing order and >=0.
y	Array of sample values, which will be replaced here by dispersion corrected values.
n	Nr of samples.
tau	Dispersion time constant (zero if no dispersion); in same time unit as sample times.
tmp	Array for temporary data, for at least n samples; enter NULL to let function to allocate and free the temporary space.

Definition at line 88 of file simdispersion.c.

  {
  /* 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;
}

simDispersion()

int simDispersion	(	double *	x,
		double *	y,
		const int	n,
		const double	tau1,
		const double	tau2,
		double *	tmp )

Simulate the effect of dispersion on a time-activity curve.

The units of rate constants must be related to the TAC time units; 1/min and min, or 1/sec and sec.

Returns

Function returns 0 when successful, else a value >= 1.

Author

Vesa Oikonen

See also

simC1, simC1_i, corDispersion, simDelay, simTTM, simBTAC

Parameters

x	Array of sample times; must be in increasing order.
y	Array of sample values, which will be replaced here by dispersion added values.
n	Nr of samples.
tau1	First dispersion time constant (zero if no dispersion); in same time unit as sample times.
tau2	2nd dispersion time constant (zero if no dispersion); in same time unit as sample times.
tmp	Array for temporary data, for at least n samples; enter NULL to let function to allocate and free the temporary space.

Definition at line 26 of file simdispersion.c.

  {
  /* 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;
}

simTTM()

int simTTM	(	double *	t,
		double *	c0,
		const int	n,
		const double	k,
		const int	cn,
		double *	cout )

Simulate output of n-compartmental transit-time model, at input TAC sample times.

The units of rate constant must be related to the time unit; 1/min and min, or 1/sec and sec.

See also

simTTM_i, simDispersion, simC1, simSamples, simDelay

Returns

Function returns 0 when successful, else a value >= 1.

Author

Vesa Oikonen

Parameters

t	Array of time values.
c0	Array of input values.
n	Number of values in TACs.
k	Rate constant of the model.
cn	Number of compartments (cn>0).
cout	Pointer for TAC array to be simulated; must be allocated.

Definition at line 144 of file simdispersion.c.

  {
  /* 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);
}

simTTM_i()

int simTTM_i	(	double *	t,
		double *	c0i,
		const int	n,
		const double	k,
		const int	cn,
		double *	cout )

Simulate output of n-compartmental transit-time model, at input TAC sample times.

This version uses integral of input function, enabling user to fully control the calculation of integral, which is more precise when input is based on an integrable mathematical function.

The units of rate constant must be related to the time unit; 1/min and min, or 1/sec and sec.

See also

simTTM, simDispersion, simC1_i, simSamples

Returns

Function returns 0 when successful, else a value >= 1.

Author

Vesa Oikonen

Parameters

t	Array of time values.
c0i	Array of AUC 0-t of input function.
n	Number of values in TACs.
k	Rate constant of the model.
cn	Number of compartments (cn>0).
cout	Pointer for TAC array to be simulated; must be allocated.

Definition at line 227 of file simdispersion.c.

  {
  /* 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);
}