https://github.com/cran/gstat
Tip revision: 17af27a3a93374490b213395eea35272fbb3bfec authored by Edzer Pebesma on 04 April 2011, 07:38:27 UTC
version 0.9-80
version 0.9-80
Tip revision: 17af27a
vario.h
#ifndef VARIO_H
# define VARIO_H /* avoid multiple inclusion */
typedef enum {
NOTSPECIFIED = 0,
SEMIVARIOGRAM,
CROSSVARIOGRAM,
COVARIOGRAM,
CROSSCOVARIOGRAM
} SAMPLE_VGM_TYPE;
extern const char *vgm_type_str[];
#define LENGTH_OF_MODEL 100 /* max string length for one variogram model */
typedef enum {
ZERO_DEFAULT = 0,
ZERO_INCLUDE,
ZERO_AVOID,
ZERO_SPECIAL
} DO_AT_ZERO;
typedef struct {
double tm[3][3]; /* 3D-transformation matrix */
double angle[3]; /* angle in <x,y>, ccl from pos x; angle up; rot. angle */
double ratio[2]; /* ratio axis2:axis1, ratio axis axis3:axis1 */
} ANIS_TM;
typedef enum {
NO_FIT = 0,
WLS_FIT,
WLS_FIT_MOD,
WLS_GNUFIT,
WLS_GNUFIT_MOD,
MIVQUE_FIT,
OLS_FIT,
WLS_NHH
} FIT_TYPE;
typedef struct {
int n_est, n_max, cloud, plot_numbers, is_asym;
int recalc, refit, pseudo, is_directional;
double *gamma, *dist;
unsigned long *nh;
double cutoff, iwidth;
SAMPLE_VGM_TYPE evt;
FIT_TYPE fit;
DO_AT_ZERO zero;
void *map, /* variogram map structure, i/o using files */
*S_grid /* variogram map structure, passed from S interface*/ ;
struct { double x, y, z; } direction;
DPOINT ***pairs;
/* optionally, the point pair list -- for j in [ 0, nh[i] >
((DPOINT ***)pairs)[i][j*2] and its successor are two
pointers to a pair of data points that were used to calculate gamma[i].
The length of pairs is (at least) nh[i] * 2. See register_pairs() */
} SAMPLE_VGM; /* a sample variogram */
typedef enum {
NOT_SP = 0,
NUGGET,
EXPONENTIAL,
SPHERICAL,
GAUSSIAN,
EXCLASS,
#ifdef USING_R
MATERN,
STEIN,
#endif
CIRCULAR,
LINEAR,
BESSEL,
PENTASPHERICAL,
PERIODIC,
HOLE,
LOGARITHMIC,
POWER,
SPLINE,
LEGENDRE,
MERROR,
INTERCEPT
} VGM_MODEL_TYPE;
typedef struct {
VGM_MODEL_TYPE model;
const char *name, *name_long, *v_gnuplot, *c_gnuplot;
double (*fn)(double h, double *r), /* variogram value at h of basic model */
(*da_fn)(double h, double *r); /* it's derivative to the range parm. */
} V_MODEL;
extern const V_MODEL v_models[];
#define NRANGEPARS 2 /* number of range parameters in variogram models */
typedef struct {
VGM_MODEL_TYPE model;
int fit_sill, fit_range, id;
double range[NRANGEPARS], sill,
(*fnct)(double h, double *r), /* (partial) unit variogram function */
(*da_fnct)(double h, double *r); /* (partial) derivative to range of unit variogram */
ANIS_TM *tm_range;
} VGM_MODEL;
typedef struct {
long n; /* length */
double maxdist, *values;
ANIS_TM *tm_range;
} COV_TABLE;
#define COV_TABLE_VALUE(tablep, dist) \
(dist >= tablep->maxdist ? tablep->values[tablep->n - 1] : \
tablep->values[(int) floor(tablep->n * (dist / tablep->maxdist))])
#define SEM_TABLE_VALUE(tablep, dist) \
(tablep->values[0] - COV_TABLE_VALUE(tablep, dist))
typedef struct {
char *descr, *fname, *fname2; /* descript. and sample variogram (maps) */
int n_models, max_n_models, n_fit, id, id1, id2,
block_semivariance_set, block_covariance_set, isotropic,
is_valid_covariance, fit_is_singular;
VGM_MODEL *part; /* the basic models */
COV_TABLE *table; /* covariance value table */
double block_semivariance, /* average within-block semivariance */
block_covariance, /* average within-block covariance */
max_range, /* maximum range: where sill is reached */
sum_sills, /* sum of partial sill's */
measurement_error, /* measurement error value--def. zero */
max_val, /* maximum value that is ever reached */
min_val, /* minimum value that is ever reached */
SSErr; /* fit result */
SAMPLE_VGM *ev;
} VARIOGRAM;
#define dist2(x,y,z) (x*x+y*y+z*z)
#define TM_IS3D(tm) \
(tm->angle[1] != 0.0 || tm->angle[2] != 0.0 || tm->ratio[1] < 1.0)
#define relative_norm(v,x,y,z) ((v == NULL || dist2(x,y,z) == 0.0) ? 1.0 : \
(transform_norm(NULL,x,y,z)/transform_norm(v,x,y,z)))
#define UnitCovariance(part,x,y,z) \
(part.model == INTERCEPT ? (1.0) :\
(1.0 - part.fnct(transform_norm(part.tm_range,x,y,z), part.range)))
#define UnitSemivariance(part,x,y,z) \
(part.fnct(transform_norm(part.tm_range,x,y,z), part.range))
#define Covariance(part,x,y,z) (part.sill * UnitCovariance(part,x,y,z))
#define Semivariance(part,x,y,z) (part.sill * UnitSemivariance(part,x,y,z))
#define DA_DELTA 0.001
#define da_Semivariance(part,x,y,z) \
(part.da_fnct != NULL ? \
(part.sill * part.da_fnct(transform_norm(part.tm_range,x,y,z), part.range)) : \
da_general(&(part), transform_norm(part.tm_range, x,y,z)))
/*
(part.sill * \
(part.fnct(transform_norm(part.tm_range,x,y,z), part.range * (1.0 + DA_DELTA)) - \
part.fnct(transform_norm(part.tm_range,x,y,z), part.range * (1.0 - DA_DELTA))) \
/ (2 * part.range * DA_DELTA)))
*/
#define EPSILON 1.0e-30 /* for small, non-zero anisotropy ratios */
#define is_covariogram(v) \
((v->ev->evt==COVARIOGRAM||v->ev->evt==CROSSCOVARIOGRAM))
#define is_variogram(v) \
((v->ev->evt==SEMIVARIOGRAM||v->ev->evt==CROSSVARIOGRAM))
#define is_direct(v) \
((v->ev->evt==COVARIOGRAM||v->ev->evt==SEMIVARIOGRAM))
#define is_cross(v) \
((v->ev->evt==CROSSCOVARIOGRAM||v->ev->evt==CROSSVARIOGRAM))
#define MODELHASNORANGE(m) (m == NUGGET || m == INTERCEPT)
#define PARTHASNORANGE(m) (m->model == NUGGET || m->model == INTERCEPT || \
(m->model == LINEAR && m->range == 0.0))
#if defined(__cplusplus)
extern "C" {
#endif
extern const char *v_model_gnuplot[];
extern const char *v_model_names[];
extern const char *c_model_gnuplot[];
VARIOGRAM *init_variogram(VARIOGRAM *v);
SAMPLE_VGM *init_ev(void);
void vgm_init_block_values(VARIOGRAM *v);
void free_variogram(VARIOGRAM *v);
void logprint_variogram(const VARIOGRAM *v, int verbose);
void fprint_variogram(FILE *f, const VARIOGRAM *v, int verbose);
const char *sprint_variogram(const VARIOGRAM *v, int verbose);
double get_semivariance(const VARIOGRAM *v, double dx, double dy, double dz);
double get_covariance (const VARIOGRAM *v, double dx, double dy, double dz);
double transform_norm(const ANIS_TM *tm, double dx, double dy, double dz);
void check_variography(const VARIOGRAM **v, int n);
void update_variogram(VARIOGRAM *vp);
double get_max_sills(int n);
double da_general(VGM_MODEL *p, double h);
double effective_range(const VARIOGRAM *v);
int push_variogram_model(VARIOGRAM *v, VGM_MODEL part);
VGM_MODEL_TYPE which_variogram_model(const char *m);
double relative_nugget(VARIOGRAM *v);
int vario(int argc, char **argv);
FIT_TYPE fit_int2enum(int fit);
DO_AT_ZERO zero_int2enum(int zero);
DO_AT_ZERO zero_shift(DO_AT_ZERO now, int next);
FIT_TYPE fit_shift(FIT_TYPE now, int next);
VGM_MODEL_TYPE model_shift(VGM_MODEL_TYPE now, int next);
int get_n_variogram_models(void);
void push_to_v(VARIOGRAM *v, const char *mod, double sill, double *range,
int nrangepars, double *d, int fit_sill, int fit_range);
void push_to_v_table(VARIOGRAM *v, double maxdist, int length, double *values,
double *anis);
#if defined(__cplusplus)
}
#endif
#endif /* VARIO_H */