Layers
Diffusion in heterogeneous environments
Data Structures | Macros | Functions
header.h File Reference

Header file for program 3layer. More...

#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <getopt.h>
#include <string.h>
#include <strings.h>
#include <time.h>
#include <math.h>
#include <gsl/gsl_multimin.h>
Include dependency graph for header.h:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Data Structures

struct  mse_rti_params_struct_type
 
struct  source_struct_type
 
struct  more_sources_struct_type
 

Macros

#define MAXNUM_LINES   10000
 Maximum number of lines in input file.
 
#define MAX_LINELENGTH   100
 Maximum length of any line in input file.
 
#define MAXNUM_COMMENTLINES   1000
 Maximum number of comment lines of input file to copy to output file.
 
#define MAX_COMMAND_LENGTH   1000
 Maximum number of characters of command to copy to output file.
 
#define ADDITIONAL_SOURCES_STRING_LENGTH   500
 Maximum length of string argument to additional_sources option.
 
#define FALSE   0
 FALSE assigned to 0.
 
#define TRUE   1
 TRUE assigned to 1.
 
#define PI   (4*atan(1.0))
 Sets constant PI to M_PI if that is defined; otherwise, computes PI as $ 4 \: tan^{-1}(1) $.
 
#define FARADAY   96485.3399
 Faraday constant in C/mol.
 
#define SMALLNUM   2.2204460492503131e-16
 For comparing doubles with 0: Sets SMALLNUM to DBL_EPSILON or to GSL_DBL_EPSILON if they're defined; otherwise set to another small number.
 
#define IS_ZERO(x)   (fabs(x) < (SMALLNUM) ? (TRUE) : (FALSE))
 Evaluates to TRUE if |x| < SMALLNUM, FALSE otherwise.
 
#define STREQ(s1, s2)   (strcmp(s1,s2) == 0)
 Compares strings s1 and s2; evaluates to TRUE if strings are equal, FALSE if not.
 
#define SQR(x)   ((x) * (x))
 Computes the square of x.
 
#define MAX(x, y)   ((x) > (y) ? (x) : (y))
 Computes the maximum of x and y.
 
#define MIN(x, y)   ((x) < (y) ? (x) : (y))
 Computes the minimum of x and y.
 
#define INDEX(i, j)   ((i)*(nr+1)+(j))
 Computes the 1D index for concentration or alpha arrays, given pseudo-2D indices i (z index) and j (r index); uses column-major ordering.
 
#define INDEX_FULL_P(i, j)   ((i)*(2*nr-1)+nr+(j)-2)
 Computes the 1D index for r >= 0 part of a concentration image, given pseudo-2D indices i (z index) and j (r index); uses column-major ordering.
 
#define INDEX_FULL_N(i, j)   ((i)*(2*nr-1)+nr-(j))
 Computes the 1D index for r >= 0 part of a concentration image, given pseudo-2D indices i (z index) and j (r index); uses column-major ordering.
 

Functions

void convolve3 (int M, int N, double *a, double scale1, double scale2, double *invr, double *out)
 Calculates updates to the concentration in a layer by applying the Laplacian in cylindrical coordinates; also scales by $ \theta D_{free} $. More...
 
void error (char *errorstring,...)
 Print an error message to stderr and exit with status EXIT_FAILURE. More...
 
void print_usage (char *program)
 Print usage message and exit with status EXIT_FAILURE. More...
 
double * create_array (int N, char *string)
 Create an array of doubles of the specified size. More...
 
void get_filename (char *in, char *out)
 Make sure that the filename is not too long. More...
 
void get_io_filenames (char *argstring, const char *inf_extension, const char *outf_extension, char *infilename, char *outfilename)
 Determine the input filename and the default output file name. More...
 
int assemble_command (int argc, char *argv[], char *command)
 Create a string containing the command used to run the program. More...
 
double read_source_parameter (char *string, int nsource)
 Read a parameter from the string argument to the additional_sources option. More...
 
void calc_diffusion_curve_layer (int nt, int nz, int nr, int iprobe, int jprobe, int iz1, int iz2, int nolayer, double dt, double dr, double sdelay, double sduration, double alpha_so, double theta_so, double kappa_so, double alpha_sp, double theta_sp, double kappa_sp, double alpha_sr, double theta_sr, double kappa_sr, double dfree, double *t, double *s, double *invr, char *imagebasename, double image_spacing, double *p)
 Calculates the concentration as a function of space and time and returns the probe concentration as a function of time; also outputs the concentration as images if that option was chosen. More...
 
double calc_mse_rti (const gsl_vector *x, void *params)
 Mean squared error function for simplex fitting. More...
 
void rti_theory (int nt, double spdist, double samplitude, double sdelay, double sduration, double kappa, double dfree, double alpha, double theta, double *t, double *p_theory)
 Calculates RTI data for diffusion in an isotropic, homogeneous environment (direct calculation from an equation) More...
 

Detailed Description

Header file for program 3layer.

Author
David Lewis, CABI, NKI
Copyright
GNU Public License
Date
2012-2013

Definition in file header.h.

Function Documentation

int assemble_command ( int  argc,
char *  argv[],
char *  command 
)

Create a string containing the command used to run the program.

It can be very useful for the program's output file to have a comment line that contains the command used to run the program. This function assembles that command into a string from the argc and argv parameters so it can be added to the output file.

Parameters
[in]argcArgument count (parameter passed to main)
[in]argvArgument vector (parameter passed to main)
[out]commandString containing the command line
Returns
Number of elements on the command line

Definition at line 91 of file io.c.

Here is the caller graph for this function:

void calc_diffusion_curve_layer ( int  nt,
int  nz,
int  nr,
int  iprobe,
int  jprobe,
int  iz1,
int  iz2,
int  nolayer,
double  dt,
double  dr,
double  sdelay,
double  sduration,
double  alpha_so,
double  theta_so,
double  kappa_so,
double  alpha_sp,
double  theta_sp,
double  kappa_sp,
double  alpha_sr,
double  theta_sr,
double  kappa_sr,
double  dfree,
double *  t,
double *  s,
double *  invr,
char *  imagebasename,
double  image_spacing,
double *  p 
)

Calculates the concentration as a function of space and time and returns the probe concentration as a function of time; also outputs the concentration as images if that option was chosen.

This function solves the diffusion equation in each layer $ k $ ,

\[ \frac{\partial c_k(\vec r, t)}{\partial t} = D_{free} \theta_k \nabla^2{c_k(\vec r, t)} + s/\alpha_k - \kappa_k c_k(\vec r, t) \quad , \]

subject to the continuity conditions at the interfaces between layers

\[ c(\vec r_k, t) = c(\vec r_{k+1}, t) \quad \]

\[ \alpha_k \theta_k \nabla c(\vec r_k, t) = \alpha_{k+1} \theta_{k+1} \nabla c(\vec r_{k+1}, t) \quad \]

and the boundary condition $ c(\vec r, t) = 0 $ (total absorption) at the top, the bottom, and the side of the cylinder, where

$ c_k(\vec r, t) $ = concentration in layer $ k $

$ D_{free} $ = free diffusion coefficient

$ \theta_k $ = permeability in layer $ k $

$ s_k(\vec r, t) $ = source in layer $ k $

$ \alpha_k $ = extracellular volume fraction in layer $ k $

$ \kappa_k $ = nonspecific clearance factor in layer $ k $

This function calls convolve3() to compute the Laplacian in cylindrical coordinates.

Parameters
[in]ntNumber of support points in time
[in]nzNumber of rows of concentration matrix
[in]nrNumber of columns of concentration matrix
[in]iprobez-index of probe location
[in]jprober-index of probe location
[in]iz1z-index of SR-SP boundary
[in]iz2z-index of SP-SO boundary
[in]nolayerFlag for whether to model a single homogeneous environment (true) or to model a 3-layer environment
[in]dtSpacing in time ( $ t_{i+1} = t_i + \Delta t $)
[in]drSpacing in r (in this program, $ \Delta z = \Delta r $)
[in]sdelaySource delay (time before source starts)
[in]sdurationDuration of source
[in]alpha_soExtracellular volume fraction in SO layer
[in]theta_soPermeability in SO layer
[in]kappa_soNonspecific clearance factor in SO layer
[in]alpha_spExtracellular volume fraction in SP layer
[in]theta_spPermeability in SP layer
[in]kappa_spNonspecific clearance factor in SP layer
[in]alpha_srExtracellular volume fraction in SR layer
[in]theta_srPermeability in SR layer
[in]kappa_srNonspecific clearance factor in SR layer
[in]dfreeFree diffusion coefficient
[in]tTime array
[in]sSource array
[in]invrArray for $ 1/r $ values
[out]pProbe array (concentration as a function of time)

Definition at line 101 of file model.c.

Here is the call graph for this function:

Here is the caller graph for this function:

double calc_mse_rti ( const gsl_vector *  x,
void *  params 
)

Mean squared error function for simplex fitting.

The 3layer program uses a minimization function in GSL to minimize the mean squared error between the diffusion curve calculated from the multilayer model and a diffusion curve calculated from a formula that assumes an isotropic homogeneous environment. It does so to calculate the apparent parameters and characteristic diffusion curve (for comparison purposes).

This function calculates the mean squared error between the two curves. The GSL minimization routine needs the first parameter to be a vector of doubles which represent the parameters to be fit and the second parameter to be a pointer to parameters of the function to fit (rti_theory).

The second parameter actually points to a struct of parameters. The characteristic curve calculated from rti_theory is stored in the p_theory array of this struct.

Parameters
[in]xVector of doubles representing parameters to fit (alpha and theta)
[in,out]paramsVector of parameters needed by rti_theory
Returns
Mean squared error between multilayer curve and characteristic curve

Definition at line 42 of file rti-theory.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void convolve3 ( int  M,
int  N,
double *  a,
double  scale1,
double  scale2,
double *  invr,
double *  out 
)

Calculates updates to the concentration in a layer by applying the Laplacian in cylindrical coordinates; also scales by $ \theta D_{free} $.

Calculates the update terms

\[ \mathbf{out} = s_1 \mathbf{L} * \mathbf{a} + \frac{s_2}{r} \mathbf{D} * \mathbf{a} \qquad (r \neq 0) \]

or

\[ \mathbf{out} = s_1 \mathbf{L0} * \mathbf{a} \qquad (r = 0) \]

where

$ \mathbf{out} $ = output array

$ s_1, s_2 $ = scaling factors

$ \mathbf{a} $ = input array (concentration)

$ * $ means convolution

$ \mathbf{L} $ = 3x3 Laplace operator (see above)

$ \mathbf{L0} $ = 3x3 Laplace operator modified for $ r = 0 $

$ \mathbf{D} $ = 3x3 derivative operator for $ r $

Note that the index $ j = 1 $ corresponds to $ r = 0 $.

Parameters
[in]MNumber of columns of input matrix (z)
[in]NNumber of rows of input matrix (r)
[in]aInput matrix (concentration, c)
[in]scale1Scaling factor #1
[in]scale2Scaling factor #2
[in]invrVector of 1/r values
[out]outOutput matrix

Definition at line 133 of file convo.c.

double* create_array ( int  N,
char *  string 
)

Create an array of doubles of the specified size.

Also initialize the array to all 0.

Parameters
[in]NNumber of elements in array
[in]stringString containing a description of the array (for debugging purposes)
Returns
Pointer to the new array

Definition at line 105 of file extras.c.

void error ( char *  errorstring,
  ... 
)

Print an error message to stderr and exit with status EXIT_FAILURE.

Parameters
[in]errorstringString containing error message (variable-length argument list, like printf)

Definition at line 20 of file extras.c.

void get_filename ( char *  in,
char *  out 
)

Make sure that the filename is not too long.

If the input string is not to big, copy it to the output string. Otherwise, exit with an error message. This function is used to check that the user input filename is not too long.

Parameters
[in]inInput string = the user input filename (or base filename)
[out]outOutput string = the input string, if not too long

Definition at line 25 of file io.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void get_io_filenames ( char *  argstring,
const char *  inf_extension,
const char *  outf_extension,
char *  infilename,
char *  outfilename 
)

Determine the input filename and the default output file name.

From the last argument on the command line, determine the input file name and the default output file name. (If the output file name is specified on the command line, it will be used instead of the default output file name determined here.)

If the argument has an extension, it is the input filename; remove the extension to get the basename and add the default output filename extension to get the default output filename.

If the argument has no extension, it is the basename; then add the default extensions to get the input filename and the default output filename.

Parameters
[in]argstringFinal command-line argument
[in]inf_extensionDefault extension for input file
[in]outf_extensionDefault extension for output file
[out]infilenameName of input file
[out]outfilenameName of output file

Definition at line 60 of file io.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void print_usage ( char *  program)

Print usage message and exit with status EXIT_FAILURE.

Parameters
[in]programString containing the name of the program

Definition at line 40 of file extras.c.

Here is the caller graph for this function:

double read_source_parameter ( char *  string,
int  nsource 
)

Read a parameter from the string argument to the additional_sources option.

The additional_sources option to 3layer takes several arguments. (The total number of arguments depends on the number of extra sources.) The getopt_long function (used in the main() function to parse 3layer's options) takes only one argument per option, so for additional_sources's argument use a string that contains the several arguments that additional_sources needs.

This function uses strtok to parse the string to obtain the next parameter. The two arguments for this function are used only for diagnostic purposes – if there is an error, the function aborts with an error message that lists these two arguments.

The first element of the string (the number of extra sources) was previously obtained with the first call to strtok, and subsequent calls to strtok (via this function) use NULL rather than the name of the string.

Parameters
[in]stringString containing the type of parameter to be read
[out]nsourceThe number of the source whose parameter is being read
Returns
Value of the parameter that is read

Definition at line 153 of file io.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void rti_theory ( int  nt,
double  spdist,
double  samplitude,
double  sdelay,
double  sduration,
double  kappa,
double  dfree,
double  alpha,
double  theta,
double *  t,
double *  p_theory 
)

Calculates RTI data for diffusion in an isotropic, homogeneous environment (direct calculation from an equation)

Todo:
Change formula to take kappa into account
Parameters
[in]ntNumber of time points of calculation
[in]spdistDistance between source and probe
[in]samplitudeAmplitude of source
[in]sdelaySource delay (time before source starts)
[in]sdurationDuration of source
[in]kappaNonspecific clearance factor
[in]dfreeFree diffusion coefficient
[in]alphaExtracellular volume fraction
[in]thetaPermeability
[in]tTime array
[out]p_theoryProbe array (concentration as a function of time)

Definition at line 88 of file rti-theory.c.

Here is the caller graph for this function:

Generated on Sun Aug 6 2017 18:21:12 for Layers by   doxygen 1.8.11