#include #include #include #include #ifndef _GEOSTAT_H #define _GEOSTAT_H #define NTAB 32 /* Modified january, the 22th 2004 */ /* from float to double */ /* List of structures: */ /* ------------------- */ /* vario_mod */ /* variotable_mod*/ /* grid_mod */ /* welldata_mod */ /* statfacies_mod */ /* inequalities_mod */ /* statistic_mod */ /* grad_mod */ /* gradients_mod */ /* cdf_mod */ /* realization_mod */ /* List of functions: */ /* ------------------ */ /* axes, cardsin, choldc, coordinates */ /* covariance, cov_matrix, cov_value, */ /* cubic, cutspectr, deflimit, dual_kri */ /* exponential, fourt, funtrun1, G, gammf */ /* gammln, gammp, gasdev, gaussian, gcf, */ /* gen_cov_matrix, Ginv, gradual, cgrid, */ /* gser, invtrun1, krig_stat, length, */ /* LtimeZ, mat_vec, maxfactor, metrop, norm */ /* normal, nugget, power, ran2, scal_vect, */ /* solve3, sort, spherical, stable, statlog2nor */ /* test_fract, trun1, trungasdev,vec_vec, */ /* vf2gthres,polint */ /*STRUCTURES*/ /*----------*/ /*variogram */ /*Nvario: number of combined variogram models */ /*vario: model of variogram per variogram model */ /* 1 --> exponential */ /* 2 --> gaussian */ /* 3 --> spherical */ /* 4 --> cardsin */ /* 5 --> stable */ /* 6 --> gamma */ /* 7 --> cubic */ /* 8 --> nugget */ /* 9 --> power */ /*alpha: exponent for the stable and gamma variogram */ /* per variogram model */ /*ap: anisotropy axes per variogram model */ /*scf: correlation lengths per variogram model */ /*var: normalized variance per variogram model(sum = 1)*/ struct vario_mod { int Nvario; int* vario; double* alpha; double* ap; double* scf; double* var; }; /*variogram table */ /*Nvario: number of combined variogram models */ /*vario: model of variogram per variogram model */ /* 1 --> exponential */ /* 2 --> gaussian */ /* 3 --> spherical */ /* 4 --> cardsin */ /* 5 --> stable */ /* 6 --> gamma */ /* 7 --> cubic */ /* 8 --> nugget */ /* 9 --> power */ /*alpha: exponent for the stable and gamma variogram */ /* per variogram model */ /*ap: anisotropy axes per variogram model */ /*scf: correlation lengths per variogram model */ /*var: normalized variance per variogram model(sum = 1)*/ struct variotable_mod { int number_of_variograms; int* Nvario; int* vario; float* alpha; float* ap; float* scf; float* var; }; /*grid */ /*NX: number of gridblocks along the X axis*/ /*NY: number of gridblocks along the Y axis*/ /*NZ: number of gridblocks along the Z axis*/ /*DX: gridblock length along the X axis */ /*DY: gridblock length along the Y axis */ /*DZ: gridblock length along the Z axis */ /*Xo: X-cell number of the origin cell */ /*Yo: Y-cell number of the origin cell */ /*Zo: Z-cell number of the origin cell */ struct grid_mod { int NX, NY, NZ; double DX, DY, DZ; double Xo, Yo, Zo; }; /*well data */ /*nwell: number of wells */ /*n: number of measurement points per well */ /* i = [0...nwell-1] */ /*ntype: number of measurement types */ /*code: status of the measurements i=[0...ntype-1] */ /* --> 0 : Gaussian white noise */ /* --> 1: standard Normal */ /* --> 2: non standard Normal */ /* --> 3: lognormal (neperien) */ /* --> 4: lognormal (log10) */ /* --> 5: facies */ /* --> 6: uniform */ /* --> 7: any */ /*x: X-coordinates of the measurements */ /* i = [0 ... n[0]-1 n[0] ... n[0]+n[1]-1...sum(n[k])-1]*/ /*y: Y-coordinates of the measurements */ /* i = [0 ... n[0]-1 n[0] ... n[0]+n[1]-1...sum(n[k])-1]*/ /*z: Z-coordinates of the measurements */ /* i = [0 ... n[0]-1 n[0] ... n[0]+n[1]-1...sum(n[k])-1]*/ /*measure: values of the measurements */ /* same kind of indexation, but repeated per type of */ /* measurement */ /* type 1 : */ /* i = [0 ... n[0]-1 n[0] ... n[0]+n[1]-1...sum(n[k])-1]*/ /* type 2 : */ /* i=[sum(n[k])... sum(n[k])+n[0]-1 ... 2*(sum(n[k])-1)]*/ struct welldata_mod { int nwell; int* n; int ntype; int* code; float* x; float* y; float* z; float* measure; }; /*volume fractions for facies */ /*ncat: number of facies */ /*nblock: number of gridblocks with different */ /* volume fractions */ /* = 1 --> stationary case */ /* = NX*NY*NZ --> nonstationary case */ /*vf: volume fractions for the first ncat-1 facies*/ /* i = [0...ncat-2]*nblock */ struct statfacies_mod { int ncat; int nblock; float* vf; }; /*inequalities for truncated plurigaussian realizations*/ /*only two basic realizations Y1 and Y2 are considered */ /*Y1 and Y2 are independent */ /*nsY1: number of unknown thresholds for Y1 */ /*nsY2: number of unknown thresholds for Y2 */ /*thresholds: vector with the thresholds and -10 and 10*/ /* the output values are the Gaussian */ /* thresholds */ /* i = [0 ... n+1], n = nsY1+nsY2 */ /* thresholds[n] = -10,thresholds[n+1] = 10 */ /* given at the beginning */ /*address_sY1: successive upper and lower bounds for */ /* the different facies for Y1 */ /* i = [0 ... 2*ncat-1] */ /* the values in address_sY1 are integers */ /* ranging from 0 to n+1 with n = nsY1+nsY2*/ /*address_sY2: successive upper and lower bounds for */ /* the different facies for Y2 */ /* i = [0 ... 2*ncat-1] */ /* the values in address_sY2 are integers */ /* ranging from 0 to n+1 with n = nsY1+nsY2*/ struct inequalities_mod { int nsY1; int nsY2; float* thresholds; int* address_sY1; int* address_sY2; }; /*statistical data */ /*type --> 0 : normal */ /* --> 1 : natural log */ /* --> 2 : log 10 */ /*nblock_mean: number of gridblocks with different */ /* means */ /* = 1 --> stationary case */ /* = NX*NY*NZ --> nonstationary case */ /*mean: mean of the variable i = [0...nblock_mean] */ /* DHF CHANGE: FOR TYPE 1 AND TYPE 2, MEAN IS LOG MEAN ! */ /*nblock_var: number of gridblocks with different */ /* variances */ /* = 1 --> stationary case */ /* = NX*NY*NZ --> nonstationary case */ /*variance: variance of the variable i = [0...nblock_var]*/ /* DHF CHANGE: FOR TYPE 1 AND TYPE 2, VAR IS LOG VAR ! */ struct statistic_mod { int type; int nblock_mean; double* mean; int nblock_var; double* variance; }; /*gradual deformation parameters */ /*Nadded: number of complementary realizations */ /*NZONES: number of subregions */ /*rho: gradual deformation parameters */ /*rho[NZONES*(0...Nadded)] */ /*cellini[NZONES*(0...2)] lower cell bound for */ /*for subregions along axes X,Y,Z */ /*cellfin[NZONES*(0...2)] upper cell bound for */ /*for subregions along axes X,Y,Z */ struct grad_mod { int Nadded, NZONES; float* rho; int *cellini, *cellfin; }; /*gradient structures */ /*Nparam : number of parameters for which gradients are */ /* required */ /*Ncells : number of cells for which gradients are */ /* calculated */ /*grad : vector with the calculated gradients */ /* dimension = Nparam*Ncells */ /* this vector is organized as */ /* 0 1...Ncells-1 for the first parameter followed*/ /* Ncells....2*Ncells-1 for the second parameter */ /* and so on */ struct gradients_mod { int Nparam, Ncells; float* grad; }; /*description of discretized cumulative distributions */ /*n: number of points */ /*x: values along the x axis i = [0...n-1] */ /*fx: corresponding values for the cumulative */ /* distribution i = [0...n-1] */ struct cdf_mod { int n; float* x; float* fx; }; /*realization */ /*n: number of components */ /*code: status of the realization */ /* --> 0 : Gaussian white noise */ /* --> 1: standard Normal */ /* --> 2: non standard Normal */ /* --> 3: lognormal (neperien) */ /* --> 4: lognormal (log10) */ /* --> 5: facies */ /* --> 6: conditional standard Normal */ /* --> 7: conditional non standard Normal */ /* --> 8: conditional lognormal (neperien) */ /* --> 9: conditional lognormal (log10) */ /* --> 10: conditional facies */ /* --> 11: uniform numbers */ /* --> 12: conditional uniform numbers */ /* --> 13: unconditional any type */ /* --> 14: conditional any type */ /*vector: realization vector i = [0...n-1] */ struct realization_mod { int n; int code; double* vector; }; /*=====================================================*/ /*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, int cores); /*Cholesky decomposition of matrix C */ /* C : symetric positive-definite matrix recorded */ /* (per raws) as a vector with only components */ /* Cij so that j <= i, 0 <= i <= n */ /* n : dimension of matrix Cij */ /* */ /* 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 */ /*input: */ /*maille: number of the cell to identify */ /*grid: structure defining the grid */ /*output: */ /*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*/ /*1+NX*NY*NZ, covar[0] is a dead cell */ /*variogram: structure defined above */ /*grid: structure defined above */ /*n: number of gridblocks along X,Y and Z*/ void covariance(double* covar, struct vario_mod variogram, struct grid_mod grid, int n[3], int cores); /*computation of the covariance matrix for the well data*/ /*well coordinates are given as a number of cells */ /*The dimension of C is given by well.nwell */ /*C is recorded as a vector so that */ /*C[k] = Cij with i = [0...nwell-1], j = [0...nwell-1] */ /*and k = j+i(i+1)/2 */ /*variogram: structure defined above */ /*well: structure defined above */ /*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: */ /* 1 -> exponential */ /* 2 -> gaussian */ /* 3 -> spherical */ /* 4 -> cardsin */ /* 5 -> stable */ /* 6 -> gamma */ /* 7 -> cubic */ /* 8 -> nugget */ /* 9 -> power */ /*variogram: variogram with the structure defined above*/ /*di: distance along the X axis */ /*dj: distance along the Y axis */ /*dk: distance along the Z axis */ /* The returned value is the computed covariance value */ double cov_value(struct vario_mod variogram, double di, double dj, double dk, int cores); /*cubic covariance value for lag h*/ double cubic(double h, int cores); /*truncation of the power spectrum to remove */ /*high frequencies - isotropic case */ /*covar: power spectrum */ /*kx: number of cells to save along the x-axis */ /*ky: number of cells to save along the y-axis */ /*kz: number of cells to save along the z-axis */ /*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 */ /*x: facies */ /*thresholds: Gaussian threshold vector */ /*facies: structure defined above */ /*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 */ /*k: type of measure */ /*well: structure defining the well data */ /*grid: structure defined above */ /*D: weight vector of length Ndata, Di, i = 0...Ndata-1*/ /*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 */ /*nn: number of gridblocks along the X,Y and Z axes */ /*ndim: number of dimensions */ /*ifrwd: non-zero for forward transform, 0 for inverse*/ /*icplx: non-zero for complex data, 0 for real */ /*workr: utility real part vector for storage */ /*worki: utility imaginary part vector for storage */ /*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, int cores); /*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, int cores); /*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 */ /*idum: seed */ double gasdev(long* idum, long* idum2, long* iy, long* iv, int cores); /*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 */ /*C defines the correlation between the first n wells */ /*C is recorded as a vector so that */ /*C[k] = Cij with i = [0...nwell-1], j = [0...nwell-1] */ /*and k = j+i(i+1)/2 */ /*variogram: structure defined above */ /*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 */ /*from Statistical Computing,by W.J. Kennedy, Jr. and */ /*James E. Gentle, 1980, p. 95. */ /*input : */ /*p: cumulative probability value */ float Ginv(float p); /*gradual combination of 1 realization + Nadded */ /*complementary realizations */ /*rho: gradual deformation parameters */ /*rho[0...NZONES*Nadded-1] */ /*rangement des vecteurs colonnes les uns apres */ /*les autres */ /*Zo: starting realization */ /*Zo[0...n] */ /*Z: complementary realizations all stored */ /*Z[0...n-1 n...2n-1 ...(Nadded-1)n...Nadded.n-1]*/ /*sequentially in a single vector */ /*Nadded: number of complementary realizations */ /*n: number of components per realization */ /*NZONES: number of subregions */ /*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*/ /*grid: structure defining the actual grid */ /*output: */ /*n: vector with the number of cells along the */ /* X, Y and Z axes for the underlying grid */ /* i = [0 1 2] */ void cgrid(struct vario_mod variogram, struct grid_mod grid, int n[3], int cores); /*incomplete gamma function evaluated by its series*/ /*representation as gamser, also returns ln(G(a)) */ /*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 */ /*u: uniform deviate between 0 and 1 */ /*lim_inf: lower bound of the considered interval*/ /*lim_sup: upper bound of the considered interval*/ /*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 */ /*scf: correlation length */ /*ap: normalized anisotropy axes */ int length(int N, int i, double* scf, double* ap, double D, int Nvari, int cores); /*calculates L.Z/ /* L : lower triangular matrix recorded */ /* (per raws) as a vector with only components */ /* Lij so that j <= i, i = [0...n-1] */ /* Z : vector, Zi avec i = [0...n-1] */ /* b : vector, bi avec i = [0...n-1] */ /* n : dimension of matrix Lij */ /* */ /* The solution vector is returned in b */ void LtimeZ(double* L, float* Z, float* b, int n); /*determines the greatest prime factor of an integer*/ int maxfactor(int n, int cores); /*metrop returns a boolean varible that issues a */ /*verdict on whether to accept a reconfiguration */ /*defined by the probability ratio "ratio". */ /*If ratio >= 1, metrop = 1(true), while if ratio < 1, */ /*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, int cores); /*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 */ /* Cij so that j <= i, i = [0...n-1] */ /* D : right-hand side vector, Di avec i = [0...n-1]*/ /* n : dimension of matrix Cij */ /* */ /* The solution vector is returned in D */ /* 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 */ /*lim_inf: lower bound of the considered interval */ /*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 */ double trungasdev(long* idum, double lim_inf, double lim_sup, long* idum2, long* iy, long iv[NTAB]); /* tb1.b2 */ /* b1 : vector */ /* b2 : vector */ /* n : length of b1 and b2, bi, i = [0...n-1]*/ /*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]*/ /*where n is the number of facies-1 */ /*USES normal*/ void vf2gthres(struct statfacies_mod facies, float* thresholds); void polint(float xa[], float ya[], int n, float x, float* y, float* dy); #endif