header.h

Go to the documentation of this file.

// Includes 
#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>


// Constants

#define MAXNUM_LINES 10000

#define MAX_LINELENGTH 100

#define MAXNUM_COMMENTLINES 1000

#define MAX_COMMAND_LENGTH 1000

#define ADDITIONAL_SOURCES_STRING_LENGTH 500

#define FALSE 0

#define TRUE 1

#ifdef M_PI
# define PI M_PI
#else
# define PI (4*atan(1.0))
#endif

#define FARADAY 96485.3399

#ifdef DBL_EPSILON
# define SMALLNUM DBL_EPSILON
#elif defined GSL_DBL_EPSILON
# define SMALLNUM GSL_DBL_EPSILON
#else
# define SMALLNUM 2.2204460492503131e-16
#endif


// Macros

#define IS_ZERO(x) (fabs(x) < (SMALLNUM) ? (TRUE) : (FALSE))

#define STREQ(s1,s2) (strcmp(s1,s2) == 0)

#define SQR(x) ((x) * (x))

#define MAX(x,y) ((x) > (y) ? (x) : (y))

#define MIN(x,y) ((x) < (y) ? (x) : (y))


/* Define a macro for calculating the 1D index of the arrays
 * from the 2D indices i and j, in order to make the code simpler.
 * The 2D indices are: i corresponds to z and j corresponds to r.
 * The number of rows in most of these matrices is (nr+1)
 */
#define INDEX(i,j) ((i)*(nr+1)+(j))


/* Define macros for calculating the 1D indices of the arrays
 * that have -rmax < r < rmax, which I output as images in model.c.
 * The parts of the array corresponding to r >= 0 are taken from the
 * concentration arrays in cylindrical coordinates (indexed with
 * INDEX(i,j)) and the parts of the array corresponding to r < 0
 * are obtained from symmetry. The first macro is for r >=0 and
 * the second is for r < 0.
 */
#define INDEX_FULL_P(i,j) ((i)*(2*nr-1)+nr+(j)-2)

#define INDEX_FULL_N(i,j) ((i)*(2*nr-1)+nr-(j))


// Struct typedefs

typedef struct {
    int nt;                    
    double spdist;             
    double samplitude;         
    double sdelay;             
    double sduration;          
    double kappa;              
    double dfree;              
    double alpha;              
    double theta;              
    double *t;                 
    double *p_model;           
    double *p_theory;          
} mse_rti_params_struct_type;

typedef struct {
    double sz;         
    double sr;         
    double crnt;       
} source_struct_type;

typedef struct {
    int n;                       
    source_struct_type *source;  
} more_sources_struct_type;




// Function prototypes

// convo.c
void convolve3(int M, int N, double *a, double scale1, double scale2, double *invr, double *out);

// extras.c
void error(char *errorstring, ...);

void print_usage(char *program);

double *create_array(int N, char *string);

//io.c
void get_filename(char *in, char *out);

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

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

double read_source_parameter(char *string, int nsource);

// model.c
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);

// rti-theory.c
double calc_mse_rti(const gsl_vector *x, void *params);

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);