|
|
|
|
@ -1,8 +1,7 @@
|
|
|
|
|
#include <stdio.h>
|
|
|
|
|
#include <stdarg.h>
|
|
|
|
|
#include <stddef.h>
|
|
|
|
|
#include <stdio.h>
|
|
|
|
|
#include <string.h>
|
|
|
|
|
#include <stdarg.h>
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
#ifndef _GEOSTAT_H
|
|
|
|
|
#define _GEOSTAT_H
|
|
|
|
|
@ -27,7 +26,6 @@
|
|
|
|
|
/* cdf_mod */
|
|
|
|
|
/* realization_mod */
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/* List of functions: */
|
|
|
|
|
/* ------------------ */
|
|
|
|
|
|
|
|
|
|
@ -44,9 +42,6 @@
|
|
|
|
|
/* test_fract, trun1, trungasdev,vec_vec, */
|
|
|
|
|
/* vf2gthres,polint */
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*STRUCTURES*/
|
|
|
|
|
/*----------*/
|
|
|
|
|
/*variogram */
|
|
|
|
|
@ -75,8 +70,6 @@ struct vario_mod {
|
|
|
|
|
double* var;
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*variogram table */
|
|
|
|
|
/*Nvario: number of combined variogram models */
|
|
|
|
|
/*vario: model of variogram per variogram model */
|
|
|
|
|
@ -104,11 +97,6 @@ struct variotable_mod {
|
|
|
|
|
float* var;
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*grid */
|
|
|
|
|
/*NX: number of gridblocks along the X axis*/
|
|
|
|
|
/*NY: number of gridblocks along the Y axis*/
|
|
|
|
|
@ -125,7 +113,6 @@ struct grid_mod {
|
|
|
|
|
double Xo, Yo, Zo;
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*well data */
|
|
|
|
|
/*nwell: number of wells */
|
|
|
|
|
/*n: number of measurement points per well */
|
|
|
|
|
@ -164,8 +151,6 @@ struct welldata_mod {
|
|
|
|
|
float* measure;
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*volume fractions for facies */
|
|
|
|
|
/*ncat: number of facies */
|
|
|
|
|
/*nblock: number of gridblocks with different */
|
|
|
|
|
@ -180,7 +165,6 @@ struct statfacies_mod {
|
|
|
|
|
float* vf;
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*inequalities for truncated plurigaussian realizations*/
|
|
|
|
|
/*only two basic realizations Y1 and Y2 are considered */
|
|
|
|
|
/*Y1 and Y2 are independent */
|
|
|
|
|
@ -210,7 +194,6 @@ struct inequalities_mod {
|
|
|
|
|
int* address_sY2;
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*statistical data */
|
|
|
|
|
/*type --> 0 : normal */
|
|
|
|
|
/* --> 1 : natural log */
|
|
|
|
|
@ -235,7 +218,6 @@ struct statistic_mod {
|
|
|
|
|
double* variance;
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*gradual deformation parameters */
|
|
|
|
|
/*Nadded: number of complementary realizations */
|
|
|
|
|
/*NZONES: number of subregions */
|
|
|
|
|
@ -251,7 +233,6 @@ struct grad_mod {
|
|
|
|
|
int *cellini, *cellfin;
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*gradient structures */
|
|
|
|
|
/*Nparam : number of parameters for which gradients are */
|
|
|
|
|
/* required */
|
|
|
|
|
@ -268,8 +249,6 @@ struct gradients_mod {
|
|
|
|
|
float* grad;
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*description of discretized cumulative distributions */
|
|
|
|
|
/*n: number of points */
|
|
|
|
|
/*x: values along the x axis i = [0...n-1] */
|
|
|
|
|
@ -281,7 +260,6 @@ struct cdf_mod {
|
|
|
|
|
float* fx;
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*realization */
|
|
|
|
|
/*n: number of components */
|
|
|
|
|
/*code: status of the realization */
|
|
|
|
|
@ -312,18 +290,15 @@ struct realization_mod {
|
|
|
|
|
/*FUNCTIONS*/
|
|
|
|
|
/*---------*/
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*normalization of the anostropy axes */
|
|
|
|
|
/*ap: anisotropy axes */
|
|
|
|
|
/*scf: correlation lengths */
|
|
|
|
|
/* The returned normalized axes are in ap */
|
|
|
|
|
void axes(double* ap, double* scf, int N);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*cardsin covariance value for lag h*/
|
|
|
|
|
double cardsin(double h);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*Cholesky decomposition of matrix C */
|
|
|
|
|
/* C : symetric positive-definite matrix recorded */
|
|
|
|
|
/* (per raws) as a vector with only components */
|
|
|
|
|
@ -333,7 +308,6 @@ double cardsin(double h);
|
|
|
|
|
/* C is turned into the lower triangular cholesky matrix*/
|
|
|
|
|
void choldc(double* C, int n);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*computes the coordinates of a given cell */
|
|
|
|
|
/*as numbers of cells along the X,Y and Z axes*/
|
|
|
|
|
/*maille = i[0]+1+i[1]*NX+i[2]*NX*NY */
|
|
|
|
|
@ -344,7 +318,6 @@ void choldc(double *C, int n);
|
|
|
|
|
/*i: vector with the coordinates */
|
|
|
|
|
void coordinates(int maille, int i[3], struct grid_mod grid);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*builds the sampled covariance function */
|
|
|
|
|
/*dimensions are even */
|
|
|
|
|
/*covar: covariance array, vector of size*/
|
|
|
|
|
@ -365,7 +338,6 @@ void covariance(double *covar,struct vario_mod variogram, struct grid_mod grid,
|
|
|
|
|
/*grid: structure defined above */
|
|
|
|
|
void cov_matrix(double* C, struct vario_mod variogram, struct welldata_mod well, struct grid_mod grid);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*calculation of the covariance value for a distance h */
|
|
|
|
|
/*defined by i,j,k */
|
|
|
|
|
/*available variogram model: */
|
|
|
|
|
@ -388,7 +360,6 @@ double cov_value(struct vario_mod variogram,double di,double dj,double dk);
|
|
|
|
|
/*cubic covariance value for lag h*/
|
|
|
|
|
double cubic(double h);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*truncation of the power spectrum to remove */
|
|
|
|
|
/*high frequencies - isotropic case */
|
|
|
|
|
/*covar: power spectrum */
|
|
|
|
|
@ -398,7 +369,6 @@ double cubic(double h);
|
|
|
|
|
/*n[3]: number of cells along the X, Y and Z axes*/
|
|
|
|
|
void cutspectr(float* covar, int kx, int ky, int kz, int n[3]);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*defines the threshold interval for a facies x*/
|
|
|
|
|
/*lim_inf: lower bound */
|
|
|
|
|
/*lim_sup: upper bound */
|
|
|
|
|
@ -408,7 +378,6 @@ void cutspectr(float *covar, int kx, int ky, int kz, int n[3]);
|
|
|
|
|
/*nblock: gridcell number of point x */
|
|
|
|
|
void deflimit(double* plim_inf, double* plim_sup, float x, float* thresholds, struct statfacies_mod facies, int nblock);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*kriges the realization considering weights */
|
|
|
|
|
/*realin: input realization */
|
|
|
|
|
/*variogram: structure defining the variogram */
|
|
|
|
|
@ -419,12 +388,9 @@ void deflimit(double *plim_inf, double *plim_sup, float x, float *thresholds, st
|
|
|
|
|
/*The kriged realization is stored in realout */
|
|
|
|
|
void dual_kri(struct realization_mod* realin, struct vario_mod variogram, struct welldata_mod well, struct grid_mod grid, double* D, struct realization_mod* realout);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*exponential covariance value for lag h*/
|
|
|
|
|
double exponential(double h);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*Fast Fourier Transform - Cooley-Tukey algorithm */
|
|
|
|
|
/*datar: real part vector - to be transformed */
|
|
|
|
|
/*datai: imaginary part vector - to be transformed */
|
|
|
|
|
@ -437,27 +403,21 @@ double exponential(double h);
|
|
|
|
|
/*The transformed data are returned in datar and datai*/
|
|
|
|
|
void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*calculates F(x) = (1/a)*exp(-x*x/2)*/
|
|
|
|
|
double funtrun1(double x);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*cumulative standard normal value*/
|
|
|
|
|
float G(float x);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*gamma covariance value for lag h and exponent alpha*/
|
|
|
|
|
double gammf(double h, double alpha);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*returns the value ln(G(x))*/
|
|
|
|
|
float gammln(float xx);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*incomplete gamma fnction*/
|
|
|
|
|
float gammp(float a, float x);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*returns a normally distributed deviate with 0 mean*/
|
|
|
|
|
/*and unit variance, using ran1(idum) as the source */
|
|
|
|
|
/*of uniform deviates */
|
|
|
|
|
@ -467,13 +427,11 @@ double gasdev(long *idum, long *idum2, long *iy, long *iv, int *iset);
|
|
|
|
|
/*gaussian covariance value for lag h*/
|
|
|
|
|
double gaussian(double h);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*incomplete gamma function evaluated by its continued */
|
|
|
|
|
/*fraction represented as gammcf, also returns ln(G(a))*/
|
|
|
|
|
/*as gln */
|
|
|
|
|
void gcf(float* gammcf, float a, float x, float* gln);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*computation of the covariance matrix for the well data*/
|
|
|
|
|
/*well coordinates have no specific unit */
|
|
|
|
|
/*The dimension of C is given by n */
|
|
|
|
|
@ -485,7 +443,6 @@ void gcf(float *gammcf, float a, float x, float *gln);
|
|
|
|
|
/*well: structure defined above */
|
|
|
|
|
void gen_cov_matrix(double* C, struct vario_mod variogram, struct welldata_mod well, int n);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*Ginv */
|
|
|
|
|
/*Computes the inverse of the standard normal cumulative*/
|
|
|
|
|
/*distribution function with a numerical approximation */
|
|
|
|
|
@ -495,7 +452,6 @@ void gen_cov_matrix(double *C, struct vario_mod variogram, struct welldata_mod w
|
|
|
|
|
/*p: cumulative probability value */
|
|
|
|
|
float Ginv(float p);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*gradual combination of 1 realization + Nadded */
|
|
|
|
|
/*complementary realizations */
|
|
|
|
|
/*rho: gradual deformation parameters */
|
|
|
|
|
@ -513,7 +469,6 @@ float Ginv(float p);
|
|
|
|
|
/*grid: grid definition */
|
|
|
|
|
void gradual(struct grad_mod grad, float* Zo, float* Z, float* Zfinal, int n, struct grid_mod grid);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*computes the size of the underlying grid for FFTs*/
|
|
|
|
|
/*input: */
|
|
|
|
|
/*variogram: structure defining the variogram model*/
|
|
|
|
|
@ -529,7 +484,6 @@ void cgrid(struct vario_mod variogram, struct grid_mod grid, int n[3]);
|
|
|
|
|
/*as gln */
|
|
|
|
|
void gser(float* gamser, float a, float x, float* gln);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*calculates x so that x = invF(u) */
|
|
|
|
|
/*F is the cumulative density function for the */
|
|
|
|
|
/*function approximating the Gaussian function */
|
|
|
|
|
@ -539,12 +493,10 @@ void gser(float *gamser, float a, float x, float *gln);
|
|
|
|
|
/*C: normalizing constant */
|
|
|
|
|
double invtrun1(double u, double lim_inf, double lim_sup, double C);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*computes the kriging mean and variance*/
|
|
|
|
|
/*for the vector bi, i = [0...n-1] */
|
|
|
|
|
void krig_stat(float* b, int n, struct vario_mod variogram, struct welldata_mod well, float* mean, float* var);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/* computes the number of gridblocks for one dimension*/
|
|
|
|
|
/*N: initial number of gridblocks */
|
|
|
|
|
/*i: considered direction */
|
|
|
|
|
@ -573,39 +525,31 @@ int maxfactor(int n);
|
|
|
|
|
/*metrop is only true with probability "ratio" */
|
|
|
|
|
int metrop(double ratio, long* idum, long* idum2, long* iy, long* iv);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*2-norm of vector b */
|
|
|
|
|
/* b : vector */
|
|
|
|
|
/* n : length of b, bi, i = [0...n-1]*/
|
|
|
|
|
/*returns the norm of b */
|
|
|
|
|
double norm(double* b, int n);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*value of g(x) where g is the normal function*/
|
|
|
|
|
double normal(double x);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*nugget covariance value for lag h*/
|
|
|
|
|
double nugget(double h);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*power covariance value for lag h and exponent alpha*/
|
|
|
|
|
double power(double h, double alpha);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*generates uniform deviates between 0 and 1*/
|
|
|
|
|
/*idum: seed */
|
|
|
|
|
double ran2(long* idum, long* idum2, long* iy, long* iv);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*calculates bt.b */
|
|
|
|
|
/* b : vector, bi, i = [0...n-1] */
|
|
|
|
|
/* n : length of b */
|
|
|
|
|
/*returns the scalar product of b*/
|
|
|
|
|
double scal_vec(double* b, int n);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*solves the set of n linear equations Cx = D */
|
|
|
|
|
/* C : symmetric positive-definite matrix recorded */
|
|
|
|
|
/* (per raws) as a vector with only components */
|
|
|
|
|
@ -617,31 +561,25 @@ double scal_vec(double *b,int n);
|
|
|
|
|
/* CONJUGATE GRADIENT method */
|
|
|
|
|
void solve3(double* C, double* D, int n);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*sorts an array [0...n-1] into ascending order using */
|
|
|
|
|
/*shell */
|
|
|
|
|
void sort(float n, float* arr);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*spherical covariance value for lag h*/
|
|
|
|
|
double spherical(double h);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*stable covariance value for lag h and exponent alpha*/
|
|
|
|
|
double stable(double h, double alpha);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*conversion of log mean and variance to nor*/
|
|
|
|
|
void statlog2nor(struct statistic_mod* pstat);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*tries factor */
|
|
|
|
|
/*pnum: number to be decomposed */
|
|
|
|
|
/*fact: suggested factor */
|
|
|
|
|
/*pmaxfac: memory to keep the greatest factor*/
|
|
|
|
|
int test_fact(int* pnum, int fact, int* pmaxfac);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*calculates the integrale of an approximate function*/
|
|
|
|
|
/*for the Gaussian function over an interval defined */
|
|
|
|
|
/*by lim_inf and lim_sup */
|
|
|
|
|
@ -649,7 +587,6 @@ int test_fact(int *pnum, int fact, int *pmaxfac);
|
|
|
|
|
/*lim_sup: upper bound of the considered interval */
|
|
|
|
|
double trun1(double lim_inf, double lim_sup);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*draws a truncated gaussian variable between lim_inf*/
|
|
|
|
|
/*and lim_sup */
|
|
|
|
|
/*idum: seed */
|
|
|
|
|
@ -662,7 +599,6 @@ double trungasdev(long *idum,double lim_inf,double lim_sup,long *idum2, long *iy
|
|
|
|
|
/*returns the norm the product tb1.b2 */
|
|
|
|
|
double vec_vec(double* b1, double* b2, int n);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*turns the volume fractions into Gaussian thresholds */
|
|
|
|
|
/*facies: structure defined above */
|
|
|
|
|
/*thresholds: output threshold vector fot i = [0...n-1]*/
|
|
|
|
|
@ -672,6 +608,4 @@ void vf2gthres(struct statfacies_mod facies,float *thresholds);
|
|
|
|
|
|
|
|
|
|
void polint(float xa[], float ya[], int n, float x, float* y, float* dy);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
#endif
|
|
|
|
|
|
chortas
3 years ago