diff --git a/fftma_module/gen/lib_src/Py_getvalues.c b/fftma_module/gen/lib_src/Py_getvalues.c old mode 100755 new mode 100644 index d33f1dc..22560a1 --- a/fftma_module/gen/lib_src/Py_getvalues.c +++ b/fftma_module/gen/lib_src/Py_getvalues.c @@ -1,110 +1,110 @@ -#include -#include -#include -#include -#include -#include -#include #include "Py_py-api.h" #include "genlib.h" -#include "simpio.h" #include "geostat.h" #include "pressure.h" +#include "simpio.h" #include "toolsIO.h" +#include +#include +#include +#include +#include +#include +#include #if PY_MAJOR_VERSION >= 3 - #define PyIntObject PyLongObject - #define PyInt_Type PyLong_Type - #define PyInt_Check(op) PyLong_Check(op) - #define PyInt_CheckExact(op) PyLong_CheckExact(op) - #define PyInt_FromString PyLong_FromString - #define PyInt_FromUnicode PyLong_FromUnicode - #define PyInt_FromLong PyLong_FromLong - #define PyInt_FromSize_t PyLong_FromSize_t - #define PyInt_FromSsize_t PyLong_FromSsize_t - #define PyInt_AsLong PyLong_AsLong - #define PyInt_AS_LONG PyLong_AS_LONG - #define PyInt_AsSsize_t PyLong_AsSsize_t - #define PyInt_AsUnsignedLongMask PyLong_AsUnsignedLongMask - #define PyInt_AsUnsignedLongLongMask PyLong_AsUnsignedLongLongMask - #define PyNumber_Int PyNumber_Long +#define PyIntObject PyLongObject +#define PyInt_Type PyLong_Type +#define PyInt_Check(op) PyLong_Check(op) +#define PyInt_CheckExact(op) PyLong_CheckExact(op) +#define PyInt_FromString PyLong_FromString +#define PyInt_FromUnicode PyLong_FromUnicode +#define PyInt_FromLong PyLong_FromLong +#define PyInt_FromSize_t PyLong_FromSize_t +#define PyInt_FromSsize_t PyLong_FromSsize_t +#define PyInt_AsLong PyLong_AsLong +#define PyInt_AS_LONG PyLong_AS_LONG +#define PyInt_AsSsize_t PyLong_AsSsize_t +#define PyInt_AsUnsignedLongMask PyLong_AsUnsignedLongMask +#define PyInt_AsUnsignedLongLongMask PyLong_AsUnsignedLongLongMask +#define PyNumber_Int PyNumber_Long #endif #if PY_MAJOR_VERSION >= 3 - #define PyBoolObject PyLongObject +#define PyBoolObject PyLongObject #endif #if PY_MAJOR_VERSION >= 3 && CYTHON_COMPILING_IN_PYPY - #ifndef PyUnicode_InternFromString - #define PyUnicode_InternFromString(s) PyUnicode_FromString(s) - #endif +#ifndef PyUnicode_InternFromString +#define PyUnicode_InternFromString(s) PyUnicode_FromString(s) +#endif #endif - -int Py_getvalues(PyObject* args, long* seed,struct grid_mod* grid,struct vario_mod* variogram,struct statistic_mod* stat) +int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario_mod* variogram, struct statistic_mod* stat) { - int i, varioNargs=12, j=0; - PyObject* listvario; - PyObject* vgr; - //char* gwnfilename; - - stat->nblock_mean=1; - stat->nblock_var=1; - stat->mean=(double*)malloc(stat->nblock_mean * sizeof(double)); - if (stat->mean == NULL) return 0; - stat->variance=(double*)malloc(stat->nblock_var * sizeof(double)); - if (stat->variance == NULL) return 0; - - - if(!PyArg_ParseTuple(args, "iiidddlO!ddi", /*"iiidddslO!ddi",*/ - &(grid->NX), - &(grid->NY), - &(grid->NZ), - &(grid->DX), - &(grid->DY), - &(grid->DZ), - /*gwnfilename,*/ - seed, - &PyList_Type, &listvario, - stat->mean+0, - stat->variance+0, - &(stat->type))) return 0; - - - variogram->Nvario=PyList_Size(listvario); - - variogram->var=(double*)malloc(variogram->Nvario*sizeof(double)); - if(variogram->var==NULL) return 0; - variogram->vario=(int*)malloc(variogram->Nvario*sizeof(int)); - if(variogram->vario==NULL) return 0; - variogram->alpha=(double*)malloc(variogram->Nvario*sizeof(double)); - if(variogram->alpha==NULL) return 0; - variogram->scf=(double*)malloc(3*variogram->Nvario*sizeof(double)); - if(variogram->var==NULL) return 0; - variogram->ap=(double*)malloc(9*variogram->Nvario*sizeof(double)); - if(variogram->var==NULL) return 0; - for(i=0;iNvario;i++) - { - vgr=PyList_GetItem(listvario,i); - if(PyTuple_Size(vgr)!=12) return 0; - (variogram->var)[i]=PyFloat_AsDouble(PyTuple_GetItem(vgr,j++)); - (variogram->vario)[i]=(int)PyInt_AsLong(PyTuple_GetItem(vgr,j++)); - (variogram->alpha)[i]=PyFloat_AsDouble(PyTuple_GetItem(vgr,j++)); - (variogram->scf)[i*3+0]=PyFloat_AsDouble(PyTuple_GetItem(vgr,j++)); - (variogram->scf)[i*3+1]=PyFloat_AsDouble(PyTuple_GetItem(vgr,j++)); - (variogram->scf)[i*3+2]=PyFloat_AsDouble(PyTuple_GetItem(vgr,j++)); - (variogram->ap)[i*9+0]=PyFloat_AsDouble(PyTuple_GetItem(vgr,j++)); - (variogram->ap)[i*9+1]=PyFloat_AsDouble(PyTuple_GetItem(vgr,j++)); - (variogram->ap)[i*9+2]=PyFloat_AsDouble(PyTuple_GetItem(vgr,j++)); - (variogram->ap)[i*9+3]=PyFloat_AsDouble(PyTuple_GetItem(vgr,j++)); - (variogram->ap)[i*9+4]=PyFloat_AsDouble(PyTuple_GetItem(vgr,j++)); - (variogram->ap)[i*9+5]=PyFloat_AsDouble(PyTuple_GetItem(vgr,j++)); - } - - - - - - - return 1; + int i, varioNargs = 12, j = 0; + PyObject* listvario; + PyObject* vgr; + //char* gwnfilename; + + stat->nblock_mean = 1; + stat->nblock_var = 1; + stat->mean = (double*)malloc(stat->nblock_mean * sizeof(double)); + if (stat->mean == NULL) + return 0; + stat->variance = (double*)malloc(stat->nblock_var * sizeof(double)); + if (stat->variance == NULL) + return 0; + + if (!PyArg_ParseTuple(args, "iiidddlO!ddi", /*"iiidddslO!ddi",*/ + &(grid->NX), + &(grid->NY), + &(grid->NZ), + &(grid->DX), + &(grid->DY), + &(grid->DZ), + /*gwnfilename,*/ + seed, + &PyList_Type, &listvario, + stat->mean + 0, + stat->variance + 0, + &(stat->type))) + return 0; + + variogram->Nvario = PyList_Size(listvario); + + variogram->var = (double*)malloc(variogram->Nvario * sizeof(double)); + if (variogram->var == NULL) + return 0; + variogram->vario = (int*)malloc(variogram->Nvario * sizeof(int)); + if (variogram->vario == NULL) + return 0; + variogram->alpha = (double*)malloc(variogram->Nvario * sizeof(double)); + if (variogram->alpha == NULL) + return 0; + variogram->scf = (double*)malloc(3 * variogram->Nvario * sizeof(double)); + if (variogram->var == NULL) + return 0; + variogram->ap = (double*)malloc(9 * variogram->Nvario * sizeof(double)); + if (variogram->var == NULL) + return 0; + for (i = 0; i < variogram->Nvario; i++) { + vgr = PyList_GetItem(listvario, i); + if (PyTuple_Size(vgr) != 12) + return 0; + (variogram->var)[i] = PyFloat_AsDouble(PyTuple_GetItem(vgr, j++)); + (variogram->vario)[i] = (int)PyInt_AsLong(PyTuple_GetItem(vgr, j++)); + (variogram->alpha)[i] = PyFloat_AsDouble(PyTuple_GetItem(vgr, j++)); + (variogram->scf)[i * 3 + 0] = PyFloat_AsDouble(PyTuple_GetItem(vgr, j++)); + (variogram->scf)[i * 3 + 1] = PyFloat_AsDouble(PyTuple_GetItem(vgr, j++)); + (variogram->scf)[i * 3 + 2] = PyFloat_AsDouble(PyTuple_GetItem(vgr, j++)); + (variogram->ap)[i * 9 + 0] = PyFloat_AsDouble(PyTuple_GetItem(vgr, j++)); + (variogram->ap)[i * 9 + 1] = PyFloat_AsDouble(PyTuple_GetItem(vgr, j++)); + (variogram->ap)[i * 9 + 2] = PyFloat_AsDouble(PyTuple_GetItem(vgr, j++)); + (variogram->ap)[i * 9 + 3] = PyFloat_AsDouble(PyTuple_GetItem(vgr, j++)); + (variogram->ap)[i * 9 + 4] = PyFloat_AsDouble(PyTuple_GetItem(vgr, j++)); + (variogram->ap)[i * 9 + 5] = PyFloat_AsDouble(PyTuple_GetItem(vgr, j++)); + } + + return 1; } diff --git a/fftma_module/gen/lib_src/Py_kgeneration.c b/fftma_module/gen/lib_src/Py_kgeneration.c old mode 100755 new mode 100644 index fac79a9..6b0e4b3 --- a/fftma_module/gen/lib_src/Py_kgeneration.c +++ b/fftma_module/gen/lib_src/Py_kgeneration.c @@ -1,54 +1,49 @@ -#include -#include -#include -#include -#include -#include -#include #include "Py_py-api.h" #include "genlib.h" -#include "simpio.h" #include "geostat.h" -#include "toolsIO.h" +#include "simpio.h" #include "toolsFFTMA.h" +#include "toolsIO.h" +#include +#include +#include +#include +#include +#include +#include /* kgeneration */ /* Z is the GWN with 0-mean and 1-variance */ /* Y1 is the realization with 0-mean and variance wanted */ /* Y is the realization with mean and variance wanted */ -void Py_kgeneration(long seed,struct grid_mod grid,struct statistic_mod stat,struct vario_mod variogram,struct realization_mod *Z,struct realization_mod *Y,struct realization_mod *Y1, int n[3]) - { - int i,N; - int typelog; - - /*generate Gaussian white noise*/ - N = grid.NX*grid.NY*grid.NZ; - n[0] = 0; - n[1] = 0; - n[2] = 0; - - - - - generate(&seed,N,Z); +void Py_kgeneration(long seed, struct grid_mod grid, struct statistic_mod stat, struct vario_mod variogram, struct realization_mod* Z, struct realization_mod* Y, struct realization_mod* Y1, int n[3]) +{ + int i, N; + int typelog; + /*generate Gaussian white noise*/ + N = grid.NX * grid.NY * grid.NZ; + n[0] = 0; + n[1] = 0; + n[2] = 0; + generate(&seed, N, Z); - /*FFTMA*/ - FFTMA2(variogram,grid,n,Z,Y); + /*FFTMA*/ + FFTMA2(variogram, grid, n, Z, Y); - /*add the statistics*/ - //if (stat.mean[0] != 0 || stat.variance[0]!= 1) - //addstat2(Y,stat,Y1,Y); + /*add the statistics*/ + //if (stat.mean[0] != 0 || stat.variance[0]!= 1) + //addstat2(Y,stat,Y1,Y); - /* make a log normal realization */ - if (stat.type==1 || stat.type==2){ + /* make a log normal realization */ + if (stat.type == 1 || stat.type == 2) { - typelog=stat.type+2; - /* nor2log(Y1,typelog,Y1); */ - nor2log(Y,typelog,Y); - } + typelog = stat.type + 2; + /* nor2log(Y1,typelog,Y1); */ + nor2log(Y, typelog, Y); + } - return; - } + return; +} diff --git a/fftma_module/gen/lib_src/axes.c b/fftma_module/gen/lib_src/axes.c old mode 100755 new mode 100644 index b9e8351..e0f8216 --- a/fftma_module/gen/lib_src/axes.c +++ b/fftma_module/gen/lib_src/axes.c @@ -2,40 +2,37 @@ #include /*normalizes anisotropy axes*/ -void axes(double *ap,double *scf,int N) +void axes(double* ap, double* scf, int N) { - double sclpdt, r, eps = 1.E-6; - int i,j,k; + double sclpdt, r, eps = 1.E-6; + int i, j, k; + for (k = 0; k < N; k++) { - for (k = 0; k < N; k++) { - - r = sqrt(ap[9*k]*ap[9*k]+ap[9*k+1]*ap[9*k+1]+ap[9*k+2]*ap[9*k+2]); - ap[9*k] /= r; - ap[9*k+1] /= r; - ap[9*k+2] /= r; - sclpdt = ap[9*k]*ap[9*k+3]+ap[9*k+1]*ap[9*k+4]+ap[9*k+2]*ap[9*k+5]; - if (sclpdt > eps) { - printf("Non orthogonal axes"); - exit; - } else { - r = sqrt(ap[9*k+3]*ap[9*k+3]+ap[9*k+4]*ap[9*k+4]+ap[9*k+5]*ap[9*k+5]); - ap[9*k+3] /= r; - ap[9*k+4] /= r; - ap[9*k+5] /= r; - ap[9*k+6] = ap[9*k+1]*ap[9*k+5]-ap[9*k+2]*ap[9*k+4]; - ap[9*k+7] = ap[9*k+2]*ap[9*k+3]-ap[9*k]*ap[9*k+5]; - ap[9*k+8] = ap[9*k]*ap[9*k+4]-ap[9*k+1]*ap[9*k+3]; - for (i=0; i<3; i++) { - for (j=0; j<3; j++) { - if (scf[3*k+j] == 0.) - scf[3*k+j] = 1.; - ap[9*k+3*j+i] /= scf[3*k+j]; - } - } + r = sqrt(ap[9 * k] * ap[9 * k] + ap[9 * k + 1] * ap[9 * k + 1] + ap[9 * k + 2] * ap[9 * k + 2]); + ap[9 * k] /= r; + ap[9 * k + 1] /= r; + ap[9 * k + 2] /= r; + sclpdt = ap[9 * k] * ap[9 * k + 3] + ap[9 * k + 1] * ap[9 * k + 4] + ap[9 * k + 2] * ap[9 * k + 5]; + if (sclpdt > eps) { + printf("Non orthogonal axes"); + exit; + } else { + r = sqrt(ap[9 * k + 3] * ap[9 * k + 3] + ap[9 * k + 4] * ap[9 * k + 4] + ap[9 * k + 5] * ap[9 * k + 5]); + ap[9 * k + 3] /= r; + ap[9 * k + 4] /= r; + ap[9 * k + 5] /= r; + ap[9 * k + 6] = ap[9 * k + 1] * ap[9 * k + 5] - ap[9 * k + 2] * ap[9 * k + 4]; + ap[9 * k + 7] = ap[9 * k + 2] * ap[9 * k + 3] - ap[9 * k] * ap[9 * k + 5]; + ap[9 * k + 8] = ap[9 * k] * ap[9 * k + 4] - ap[9 * k + 1] * ap[9 * k + 3]; + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) { + if (scf[3 * k + j] == 0.) + scf[3 * k + j] = 1.; + ap[9 * k + 3 * j + i] /= scf[3 * k + j]; + } + } + } } - - } - return; + return; } - diff --git a/fftma_module/gen/lib_src/build_real.c b/fftma_module/gen/lib_src/build_real.c old mode 100755 new mode 100644 index 06a8651..f35b0c6 --- a/fftma_module/gen/lib_src/build_real.c +++ b/fftma_module/gen/lib_src/build_real.c @@ -1,10 +1,10 @@ -#include -#include -#include +#include "geostat.h" +#include #include +#include +#include #include -#include -#include "geostat.h" +#include /* build_real */ /* build a realization in the spectral domain */ @@ -18,27 +18,27 @@ /*realization: vector defining the real part */ /*ireal: vector defining the i-part */ -void build_real(int n[3],int NTOT,double *covar,double *realization,double *ireal) +void build_real(int n[3], int NTOT, double* covar, double* realization, double* ireal) { - int i,j,k,maille1; - double temp; + int i, j, k, maille1; + double temp; - /*decomposition and multiplication in the spectral domain*/ - for ( k = 1; k <= n[2]; k++) { - for (j = 1; j <= n[1]; j++) { - for (i = 1; i <= n[0]; i++) { - maille1 = i+(j-1+(k-1)*n[1])*n[0]; - temp = covar[maille1]; - if (temp > 0.) { - temp = sqrt(temp)/(double) NTOT; - } else if (temp < 0.) { - temp = sqrt(-temp)/(double) NTOT; - } - realization[maille1] *= temp; - ireal[maille1] *= temp; - } + /*decomposition and multiplication in the spectral domain*/ + for (k = 1; k <= n[2]; k++) { + for (j = 1; j <= n[1]; j++) { + for (i = 1; i <= n[0]; i++) { + maille1 = i + (j - 1 + (k - 1) * n[1]) * n[0]; + temp = covar[maille1]; + if (temp > 0.) { + temp = sqrt(temp) / (double)NTOT; + } else if (temp < 0.) { + temp = sqrt(-temp) / (double)NTOT; + } + realization[maille1] *= temp; + ireal[maille1] *= temp; + } + } } - } - - return; + + return; } diff --git a/fftma_module/gen/lib_src/cardsin.c b/fftma_module/gen/lib_src/cardsin.c old mode 100755 new mode 100644 index 504f286..8c4cd89 --- a/fftma_module/gen/lib_src/cardsin.c +++ b/fftma_module/gen/lib_src/cardsin.c @@ -1,18 +1,18 @@ -#include -#include #include "genlib.h" +#include +#include /*cardsin covariance function*/ double cardsin(double h) { - float delta = 20.371; - double z; + float delta = 20.371; + double z; - if (h != 0) { - z = (double)(h*delta); - z = sin(z)/z; - } else { - z = 1.; - } - return(z); + if (h != 0) { + z = (double)(h * delta); + z = sin(z) / z; + } else { + z = 1.; + } + return (z); } diff --git a/fftma_module/gen/lib_src/cgrid.c b/fftma_module/gen/lib_src/cgrid.c old mode 100755 new mode 100644 index 0d300aa..743b8cf --- a/fftma_module/gen/lib_src/cgrid.c +++ b/fftma_module/gen/lib_src/cgrid.c @@ -1,6 +1,5 @@ -#include #include "geostat.h" - +#include /*computes the size of the grid for FFTs*/ /*input: */ @@ -12,32 +11,32 @@ /* i = [0 1 2] */ void cgrid(struct vario_mod variogram, struct grid_mod grid, int n[3]) { - int i,N; - double D; + int i, N; + double D; - if (n == NULL || n[0] == 0 || n[1] == 0 || n[2] == 0) { - for (i = 0; i<3; i++) { - switch (i) { - case 0: - N = grid.NX; - D = grid.DX; - break; - case 1: - N = grid.NY; - D = grid.DY; - break; - case 2: - N = grid.NZ; - D = grid.DZ; - break; - } - n[i] = length(N,i,variogram.scf,variogram.ap,D,variogram.Nvario); - } - } else { - if ((n[0] -#include -#include +#include "geostat.h" +#include #include +#include +#include #include -#include -#include "geostat.h" +#include -void clean_real(struct realization_mod *realin,int n[3],struct grid_mod grid,double *vectorresult,struct realization_mod *realout) -{ - int i,j,k,maille0,maille1; - double NTOT; +void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, double* vectorresult, struct realization_mod* realout) +{ + int i, j, k, maille0, maille1; + double NTOT; - NTOT=n[0]*n[1]*n[2]; - /*is the output realization already allocated?*/ - /*if not, memory allocation*/ + NTOT = n[0] * n[1] * n[2]; + /*is the output realization already allocated?*/ + /*if not, memory allocation*/ - if ((*realout).vector == NULL || (*realout).n != (*realin).n) - { - (*realout).vector = (double *) malloc((*realin).n * sizeof(double)); - if ((*realout).vector == NULL) - { - printf("Clean_real.c: No memory available\n"); - exit; - } + if ((*realout).vector == NULL || (*realout).n != (*realin).n) { + (*realout).vector = (double*)malloc((*realin).n * sizeof(double)); + if ((*realout).vector == NULL) { + printf("Clean_real.c: No memory available\n"); + exit; + } } - - (*realout).n = (*realin).n; - (*realout).code = 1; - for ( k = 1; k <= grid.NZ; k++) { - for (j = 1; j <= grid.NY; j++) { - for (i = 1; i <= grid.NX; i++) { - maille1 = i+(j-1+(k-1)*n[1])*n[0]; - maille0 = i-1+(j-1+(k-1)*grid.NY)*grid.NX; - /* Modif du 18 juin 2003 */ - /*(*realout).vector[maille0] = vectorresult[maille1]/(double) NTOT;*/ - (*realout).vector[maille0] = vectorresult[maille1]; - } + + (*realout).n = (*realin).n; + (*realout).code = 1; + for (k = 1; k <= grid.NZ; k++) { + for (j = 1; j <= grid.NY; j++) { + for (i = 1; i <= grid.NX; i++) { + maille1 = i + (j - 1 + (k - 1) * n[1]) * n[0]; + maille0 = i - 1 + (j - 1 + (k - 1) * grid.NY) * grid.NX; + /* Modif du 18 juin 2003 */ + /*(*realout).vector[maille0] = vectorresult[maille1]/(double) NTOT;*/ + (*realout).vector[maille0] = vectorresult[maille1]; + } + } } - } - - return; + + return; } diff --git a/fftma_module/gen/lib_src/cov_value.c b/fftma_module/gen/lib_src/cov_value.c old mode 100755 new mode 100644 index 4a1398f..540c5b8 --- a/fftma_module/gen/lib_src/cov_value.c +++ b/fftma_module/gen/lib_src/cov_value.c @@ -1,54 +1,53 @@ -#include -#include "geostat.h" #include "genlib.h" +#include "geostat.h" +#include /*selection of model covariance*/ -double cov_value(struct vario_mod variogram,double di,double dj,double dk) +double cov_value(struct vario_mod variogram, double di, double dj, double dk) { - double hx,hy,hz,h; - double cov; - int k; - - cov = 0.; + double hx, hy, hz, h; + double cov; + int k; - for (k = 0; k < variogram.Nvario; k++) { + cov = 0.; - hx = di*variogram.ap[9*k]+dj*variogram.ap[9*k+1]+dk*variogram.ap[9*k+2]; - hy = di*variogram.ap[9*k+3]+dj*variogram.ap[9*k+4]+dk*variogram.ap[9*k+5]; - hz = di*variogram.ap[9*k+6]+dj*variogram.ap[9*k+7]+dk*variogram.ap[9*k+8]; - h = sqrt(hx*hx+hy*hy+hz*hz); + for (k = 0; k < variogram.Nvario; k++) { + hx = di * variogram.ap[9 * k] + dj * variogram.ap[9 * k + 1] + dk * variogram.ap[9 * k + 2]; + hy = di * variogram.ap[9 * k + 3] + dj * variogram.ap[9 * k + 4] + dk * variogram.ap[9 * k + 5]; + hz = di * variogram.ap[9 * k + 6] + dj * variogram.ap[9 * k + 7] + dk * variogram.ap[9 * k + 8]; + h = sqrt(hx * hx + hy * hy + hz * hz); - switch (variogram.vario[k]) { - case 1: - cov += variogram.var[k]*exponential(h); - break; - case 2: - cov += variogram.var[k]*gaussian(h); - break; - case 3: - cov += variogram.var[k]*spherical(h); - break; - case 4: - cov += variogram.var[k]*cardsin(h); - break; - case 5: - cov += variogram.var[k]*stable(h,variogram.alpha[k]); - break; - case 6: - cov += variogram.var[k]*gammf(h,variogram.alpha[k]); - break; - case 7: - cov += variogram.var[k]*cubic(h); - break; - case 8: - cov += variogram.var[k]*nugget(h); - break; - case 9: - cov += variogram.var[k]*power(h,variogram.alpha[k]); - break; + switch (variogram.vario[k]) { + case 1: + cov += variogram.var[k] * exponential(h); + break; + case 2: + cov += variogram.var[k] * gaussian(h); + break; + case 3: + cov += variogram.var[k] * spherical(h); + break; + case 4: + cov += variogram.var[k] * cardsin(h); + break; + case 5: + cov += variogram.var[k] * stable(h, variogram.alpha[k]); + break; + case 6: + cov += variogram.var[k] * gammf(h, variogram.alpha[k]); + break; + case 7: + cov += variogram.var[k] * cubic(h); + break; + case 8: + cov += variogram.var[k] * nugget(h); + break; + case 9: + cov += variogram.var[k] * power(h, variogram.alpha[k]); + break; + } } - } - return (cov); + return (cov); } diff --git a/fftma_module/gen/lib_src/covariance.c b/fftma_module/gen/lib_src/covariance.c old mode 100755 new mode 100644 index e8e739d..c3d13bd --- a/fftma_module/gen/lib_src/covariance.c +++ b/fftma_module/gen/lib_src/covariance.c @@ -2,88 +2,79 @@ /*builds the sampled covariance function*/ /*dimensions are even*/ -void covariance(double *covar, struct vario_mod variogram, struct grid_mod mesh, int n[3]) +void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh, int n[3]) { - int i,j,k,maille,n2[3],symmetric; - double di,dj,dk; - - for (i=0;i<3;i++) - n2[i] = n[i]/2; - - for (i=0; i<= n2[0]; i++) { - for (j=0; j<= n2[1]; j++) { - for (k=0; k<= n2[2]; k++) { - - - /*area 1*/ - maille = 1+i+n[0]*(j+n[1]*k); - di = (double)i*mesh.DX; - dj = (double)j*mesh.DY; - dk = (double)k*mesh.DZ; - covar[maille] = (double)cov_value(variogram,di,dj,dk); - - - if (k > 0 && k 0 && j 0 && i 0 && i 0 && k 0 && j 0 && i 0 && j < n2[1]) { - /*area 5*/ - di = -(double)i*mesh.DX; - dj = -(double)j*mesh.DY; - dk = (double)k*mesh.DZ; - maille = 1+(n[0]-i)+n[0]*(n[1]-j+n[1]*k); - covar[maille] = (double)cov_value(variogram,di,dj,dk); - } - - - if (k > 0 && k 0 && j < n2[1]) { - /*area 3*/ - di = (double)i*mesh.DX; - dj = -(double)j*mesh.DY; - dk = (double)k*mesh.DZ; - maille = 1+i+n[0]*(n[1]-j+n[1]*k); - covar[maille] = (double)cov_value(variogram,di,dj,dk); - } - - - if (k > 0 && k 0 && i 0 && k < n2[2] && j > 0 && j < n2[1] && i > 0 && i < n2[0]) { + /*area 2*/ + symmetric = 1 + n[0] - i + n[0] * (n[1] - j + n[1] * (n[2] - k)); + covar[symmetric] = covar[maille]; + } + + if (i > 0 && i < n2[0]) { + /*area 4*/ + di = -(double)i * mesh.DX; + dj = (double)j * mesh.DY; + dk = (double)k * mesh.DZ; + maille = 1 + (n[0] - i) + n[0] * (j + n[1] * k); + covar[maille] = (double)cov_value(variogram, di, dj, dk); + } + + if (k > 0 && k < n2[2] && j > 0 && j < n2[1]) { + /*area 8*/ + symmetric = 1 + i + n[0] * (n[1] - j + n[1] * (n[2] - k)); + covar[symmetric] = covar[maille]; + } + + if (i > 0 && i < n2[0] && j > 0 && j < n2[1]) { + /*area 5*/ + di = -(double)i * mesh.DX; + dj = -(double)j * mesh.DY; + dk = (double)k * mesh.DZ; + maille = 1 + (n[0] - i) + n[0] * (n[1] - j + n[1] * k); + covar[maille] = (double)cov_value(variogram, di, dj, dk); + } + + if (k > 0 && k < n2[2]) { + /*area 6*/ + symmetric = 1 + i + n[0] * (j + n[1] * (n[2] - k)); + covar[symmetric] = covar[maille]; + } + + if (j > 0 && j < n2[1]) { + /*area 3*/ + di = (double)i * mesh.DX; + dj = -(double)j * mesh.DY; + dk = (double)k * mesh.DZ; + maille = 1 + i + n[0] * (n[1] - j + n[1] * k); + covar[maille] = (double)cov_value(variogram, di, dj, dk); + } + + if (k > 0 && k < n2[2] && i > 0 && i < n2[0]) { + /*area 7*/ + symmetric = 1 + n[0] - i + n[0] * (j + n[1] * (n[2] - k)); + covar[symmetric] = covar[maille]; + } + } + } } - } - return; + return; } diff --git a/fftma_module/gen/lib_src/cubic.c b/fftma_module/gen/lib_src/cubic.c old mode 100755 new mode 100644 index f4247b2..9fbfd03 --- a/fftma_module/gen/lib_src/cubic.c +++ b/fftma_module/gen/lib_src/cubic.c @@ -1,17 +1,16 @@ -#include -#include #include "genlib.h" - +#include +#include /*cubic covariance function*/ double cubic(double h) { - double z; + double z; - if (h >= 1.) { - z = 0.; - } else { - z = 1.-7.*(double)(h*h)+(35./4.)*(double)(h*h*h)-3.5*(double)(h*h*h*h*h)+.75*(double)(h*h*h*h*h*h*h); - } - return (z); + if (h >= 1.) { + z = 0.; + } else { + z = 1. - 7. * (double)(h * h) + (35. / 4.) * (double)(h * h * h) - 3.5 * (double)(h * h * h * h * h) + .75 * (double)(h * h * h * h * h * h * h); + } + return (z); } diff --git a/fftma_module/gen/lib_src/exponential.c b/fftma_module/gen/lib_src/exponential.c old mode 100755 new mode 100644 index 003d1bb..c27ee4c --- a/fftma_module/gen/lib_src/exponential.c +++ b/fftma_module/gen/lib_src/exponential.c @@ -1,9 +1,9 @@ -#include -#include #include "genlib.h" +#include +#include /*exponential covariance function*/ double exponential(double h) { - return (exp(-3.*(double)h)); + return (exp(-3. * (double)h)); } diff --git a/fftma_module/gen/lib_src/fftma2.c b/fftma_module/gen/lib_src/fftma2.c old mode 100755 new mode 100644 index 7c069e8..2b4fd58 --- a/fftma_module/gen/lib_src/fftma2.c +++ b/fftma_module/gen/lib_src/fftma2.c @@ -1,8 +1,7 @@ +#include "geostat.h" +#include #include #include -#include -#include "geostat.h" - /*FAST FOURIER TRANSFORM MOVING AVERAGE METHOD */ /*Turns a Gaussian white noise vector into a */ @@ -22,84 +21,79 @@ /*output: */ /*realout: structure defining a realization - */ -void FFTMA2(struct vario_mod variogram,struct grid_mod grid,int n[3],struct realization_mod *realin,struct realization_mod *realout) +void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct realization_mod* realin, struct realization_mod* realout) { - int NTOT,i,j,k,NMAX,NDIM,ntot,nmax,NXYZ,nxyz,maille0,maille1; - int solver; - double temp; - double *ireal,*covar,*workr,*worki,*realization; - - /*covariance axis normalization*/ - axes(variogram.ap,variogram.scf,variogram.Nvario); - - /*pseudo-grid definition*/ - cgrid(variogram,grid,n); - - /*constant definition*/ - NTOT = n[0]*n[1]*n[2]; - ntot = NTOT+1; - NMAX = n[0]; - NDIM = 3; - for (i=1;i<3;i++) { - if (n[i] > NMAX) - NMAX = n[i]; - if (n[i] == 1) - NDIM--; - } - nmax = NMAX+1; - NXYZ = grid.NX*grid.NY*grid.NZ; - nxyz = NXYZ+1; - - /*array initialization*/ - covar = (double *) malloc(ntot * sizeof(double)); - testmemory(covar); - - ireal = (double *) malloc(ntot * sizeof(double)); - testmemory(ireal); - - realization = (double *) malloc(ntot * sizeof(double)); - testmemory(realization); - - workr = (double *) malloc(nmax * sizeof(double)); - testmemory(workr); - - worki = (double *) malloc(nmax * sizeof(double)); - testmemory(worki); - - /*covariance function creation*/ - covariance(covar,variogram,grid,n); - - /*power spectrum*/ - fourt(covar,ireal,n,NDIM,1,0,workr,worki); - - /*organization of the input Gaussian white noise*/ - solver=0; - prebuild_gwn(grid,n,realin,realization,solver); - - /*forward fourier transform of the GWN*/ - fourt(realization,ireal,n,NDIM,1,0,workr,worki); - - /* build realization in spectral domain */ - build_real(n,NTOT,covar,realization,ireal); - - free(covar); - - /*backward fourier transform*/ - fourt(realization,ireal,n,NDIM,0,1,workr,worki); - - - free(ireal); - free(workr); - free(worki); - - /*output realization*/ - clean_real(realin,n,grid,realization,realout); - - free(realization); - - return; -} + int NTOT, i, j, k, NMAX, NDIM, ntot, nmax, NXYZ, nxyz, maille0, maille1; + int solver; + double temp; + double *ireal, *covar, *workr, *worki, *realization; + + /*covariance axis normalization*/ + axes(variogram.ap, variogram.scf, variogram.Nvario); + + /*pseudo-grid definition*/ + cgrid(variogram, grid, n); + + /*constant definition*/ + NTOT = n[0] * n[1] * n[2]; + ntot = NTOT + 1; + NMAX = n[0]; + NDIM = 3; + for (i = 1; i < 3; i++) { + if (n[i] > NMAX) + NMAX = n[i]; + if (n[i] == 1) + NDIM--; + } + nmax = NMAX + 1; + NXYZ = grid.NX * grid.NY * grid.NZ; + nxyz = NXYZ + 1; + + /*array initialization*/ + covar = (double*)malloc(ntot * sizeof(double)); + testmemory(covar); + + ireal = (double*)malloc(ntot * sizeof(double)); + testmemory(ireal); + realization = (double*)malloc(ntot * sizeof(double)); + testmemory(realization); + workr = (double*)malloc(nmax * sizeof(double)); + testmemory(workr); + worki = (double*)malloc(nmax * sizeof(double)); + testmemory(worki); + /*covariance function creation*/ + covariance(covar, variogram, grid, n); + + /*power spectrum*/ + fourt(covar, ireal, n, NDIM, 1, 0, workr, worki); + + /*organization of the input Gaussian white noise*/ + solver = 0; + prebuild_gwn(grid, n, realin, realization, solver); + + /*forward fourier transform of the GWN*/ + fourt(realization, ireal, n, NDIM, 1, 0, workr, worki); + + /* build realization in spectral domain */ + build_real(n, NTOT, covar, realization, ireal); + + free(covar); + + /*backward fourier transform*/ + fourt(realization, ireal, n, NDIM, 0, 1, workr, worki); + + free(ireal); + free(workr); + free(worki); + + /*output realization*/ + clean_real(realin, n, grid, realization, realout); + + free(realization); + + return; +} diff --git a/fftma_module/gen/lib_src/fourt.c b/fftma_module/gen/lib_src/fourt.c old mode 100755 new mode 100644 index fe80035..8776660 --- a/fftma_module/gen/lib_src/fourt.c +++ b/fftma_module/gen/lib_src/fourt.c @@ -1,5 +1,5 @@ -#include #include +#include /*fast fourier transform */ /* THE COOLEY-TUKEY FAST FOURIER TRANSFORM */ @@ -89,503 +89,499 @@ /* PROGRAM MODIFIED FROM A SUBROUTINE OF BRENNER */ /* 10-06-2000, MLR */ - -void fourt(double *datar,double *datai, int nn[3], int ndim, int ifrwd, int icplx, double *workr,double *worki) +void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki) { - int ifact[21],ntot,idim,np1,n,np2,m,ntwo,iff,idiv,iquot,irem,inon2,non2p,np0,nprev,icase,ifmin,i,j,jmax,np2hf,i2,i1max,i3,j3,i1,ifp1,ifp2,i2max,i1rng,istep,imin,imax,mmax,mmin,mstep,j1,j2max,j2,jmin,j3max,nhalf; - double theta,wstpr,wstpi,wminr,wmini,wr,wi,wtemp,thetm,wmstr,wmsti,twowr,sr,si,oldsr,oldsi,stmpr,stmpi,tempr,tempi,difi,difr,sumr,sumi,TWOPI = 6.283185307179586476925286766559; + int ifact[21], ntot, idim, np1, n, np2, m, ntwo, iff, idiv, iquot, irem, inon2, non2p, np0, nprev, icase, ifmin, i, j, jmax, np2hf, i2, i1max, i3, j3, i1, ifp1, ifp2, i2max, i1rng, istep, imin, imax, mmax, mmin, mstep, j1, j2max, j2, jmin, j3max, nhalf; + double theta, wstpr, wstpi, wminr, wmini, wr, wi, wtemp, thetm, wmstr, wmsti, twowr, sr, si, oldsr, oldsi, stmpr, stmpi, tempr, tempi, difi, difr, sumr, sumi, TWOPI = 6.283185307179586476925286766559; ntot = 1; for (idim = 0; idim < ndim; idim++) { - ntot *= nn[idim]; + ntot *= nn[idim]; } /*main loop for each dimension*/ np1 = 1; for (idim = 1; idim <= ndim; idim++) { - n = nn[idim-1]; - np2 = np1*n; + n = nn[idim - 1]; + np2 = np1 * n; - if (n < 1) { - goto L920; - } else if (n == 1) { - goto L900; - } + if (n < 1) { + goto L920; + } else if (n == 1) { + goto L900; + } - /*is n a power of 2 and if not, what are its factors*/ - m = n; - ntwo = np1; - iff = 1; - idiv = 2; + /*is n a power of 2 and if not, what are its factors*/ + m = n; + ntwo = np1; + iff = 1; + idiv = 2; L10: - iquot = m/idiv; - irem = m-idiv*iquot; - if (iquot < idiv) - goto L50; - if (irem == 0) { - ntwo *= 2; - ifact[iff] = idiv; - iff++; - m= iquot; - goto L10; - } - idiv = 3; - inon2 = iff; + iquot = m / idiv; + irem = m - idiv * iquot; + if (iquot < idiv) + goto L50; + if (irem == 0) { + ntwo *= 2; + ifact[iff] = idiv; + iff++; + m = iquot; + goto L10; + } + idiv = 3; + inon2 = iff; L30: - iquot = m/idiv; - irem = m-idiv*iquot; - if (iquot < idiv) - goto L60; - if (irem == 0) { - ifact[iff] = idiv; - iff++; - m = iquot; - goto L30; - } - - idiv += 2; - goto L30; - + iquot = m / idiv; + irem = m - idiv * iquot; + if (iquot < idiv) + goto L60; + if (irem == 0) { + ifact[iff] = idiv; + iff++; + m = iquot; + goto L30; + } + + idiv += 2; + goto L30; + L50: - inon2 = iff; - if (irem != 0) - goto L60; - ntwo *= 2; - goto L70; - + inon2 = iff; + if (irem != 0) + goto L60; + ntwo *= 2; + goto L70; + L60: - ifact[iff] = m; - + ifact[iff] = m; + L70: - non2p = np2/ntwo; - + non2p = np2 / ntwo; - /*SEPARATE FOUR CASES-- + /*SEPARATE FOUR CASES-- 1. COMPLEX TRANSFORM 2. REAL TRANSFORM FOR THE 2ND, 3RD, ETC. DIMENSION. METHOD: TRANSFORM HALF THE DATA, SUPPLYING THE OTHER HALF BY CONJUGATE SYMMETRY. 3. REAL TRANSFORM FOR THE 1ST DIMENSION, N ODD. METHOD: SET THE IMAGINARY PARTS TO ZERO. 4. REAL TRANSFORM FOR THE 1ST DIMENSION, N EVEN. METHOD: TRANSFORM A COMPLEX ARRAY OF LENGTH N/2 WHOSE REAL PARTS ARE THE EVEN NUMBERED REAL VALUES AND WHOSE IMAGINARY PARTS ARE THE ODD-NUMBERED REAL VALUES. UNSCRAMBLE AND SUPPLY THE SECOND HALF BY CONJUGATE SYMMETRY. */ - - icase = 1; - ifmin = 1; - if (icplx != 0) - goto L100; - icase = 2; - if (idim > 1) - goto L100; - icase = 3; - if (ntwo <= np1) - goto L100; - icase = 4; - ifmin = 2; - ntwo /= 2; - n /= 2; - np2 /= 2; - ntot /= 2; - i = 1; - for (j = 1; j <= ntot; j++) { - datar[j] = datar[i]; - datai[j] = datar[i+1]; - i += 2; - } - - /*shuffle data by bit reversal, since n = 2^k. As the shuffling can be done by simple interchange, no working array is needed*/ + + icase = 1; + ifmin = 1; + if (icplx != 0) + goto L100; + icase = 2; + if (idim > 1) + goto L100; + icase = 3; + if (ntwo <= np1) + goto L100; + icase = 4; + ifmin = 2; + ntwo /= 2; + n /= 2; + np2 /= 2; + ntot /= 2; + i = 1; + for (j = 1; j <= ntot; j++) { + datar[j] = datar[i]; + datai[j] = datar[i + 1]; + i += 2; + } + + /*shuffle data by bit reversal, since n = 2^k. As the shuffling can be done by simple interchange, no working array is needed*/ L100: - if (non2p > 1) - goto L200; - np2hf = np2/2; - j = 1; - for (i2 = 1; i2 <= np2; i2 += np1) { - if (j >= i2) - goto L130; - i1max = i2+np1-1; - for (i1 = i2; i1 <= i1max; i1++) { - for (i3 = i1; i3 <= ntot; i3 += np2) { - j3 = j+i3-i2; - tempr = datar[i3]; - tempi = datai[i3]; - datar[i3] = datar[j3]; - datai[i3] = datai[j3]; - datar[j3] = tempr; - datai[j3] = tempi; - } - } - - L130: - m = np2hf; - - L140: - if (j <= m) { - j += m; - } else { - j -= m; - m /= 2; - if (m >= np1) - goto L140; - } - - } - goto L300; - - /*shuffle data by digit reversal for general n*/ + if (non2p > 1) + goto L200; + np2hf = np2 / 2; + j = 1; + for (i2 = 1; i2 <= np2; i2 += np1) { + if (j >= i2) + goto L130; + i1max = i2 + np1 - 1; + for (i1 = i2; i1 <= i1max; i1++) { + for (i3 = i1; i3 <= ntot; i3 += np2) { + j3 = j + i3 - i2; + tempr = datar[i3]; + tempi = datai[i3]; + datar[i3] = datar[j3]; + datai[i3] = datai[j3]; + datar[j3] = tempr; + datai[j3] = tempi; + } + } + + L130: + m = np2hf; + + L140: + if (j <= m) { + j += m; + } else { + j -= m; + m /= 2; + if (m >= np1) + goto L140; + } + } + goto L300; + + /*shuffle data by digit reversal for general n*/ L200: - for (i1 = 1; i1 <= np1; i1++) { - for (i3 = i1; i3 <= ntot; i3 += np2) { - j = i3; - for (i = 1; i <= n; i++) { - if (icase != 3) { - workr[i] = datar[j]; - worki[i] = datai[j]; - } else { - workr[i] = datar[j]; - worki[i] = 0.; - } - ifp2 = np2; - iff = ifmin; - L250: - ifp1 = ifp2/ifact[iff]; - j += ifp1; - if (j >= i3+ifp2) { - j -= ifp2; - ifp2 = ifp1; - iff += 1; - if (ifp2 > np1) - goto L250; - } - } - i2max = i3+np2-np1; - i = 1; - for (i2 = i3; i2 <= i2max; i2 += np1) { - datar[i2] = workr[i]; - datai[i2] = worki[i]; - i++; - } - - } - } - - /*special case-- W=1*/ + for (i1 = 1; i1 <= np1; i1++) { + for (i3 = i1; i3 <= ntot; i3 += np2) { + j = i3; + for (i = 1; i <= n; i++) { + if (icase != 3) { + workr[i] = datar[j]; + worki[i] = datai[j]; + } else { + workr[i] = datar[j]; + worki[i] = 0.; + } + ifp2 = np2; + iff = ifmin; + L250: + ifp1 = ifp2 / ifact[iff]; + j += ifp1; + if (j >= i3 + ifp2) { + j -= ifp2; + ifp2 = ifp1; + iff += 1; + if (ifp2 > np1) + goto L250; + } + } + i2max = i3 + np2 - np1; + i = 1; + for (i2 = i3; i2 <= i2max; i2 += np1) { + datar[i2] = workr[i]; + datai[i2] = worki[i]; + i++; + } + } + } + + /*special case-- W=1*/ L300: - i1rng = np1; - if (icase == 2) - i1rng = np0*(1+nprev/2); - if (ntwo <= np1) - goto L600; - for (i1 = 1; i1 <= i1rng; i1++) { - imin = np1+i1; - istep = 2*np1; - goto L330; - - L310: - j = i1; - for (i = imin; i <= ntot; i += istep) { - tempr = datar[i]; - tempi = datai[i]; - datar[i] = datar[j]-tempr; - datai[i] = datai[j]-tempi; - datar[j] = datar[j]+tempr; - datai[j] = datai[j]+tempi; - j += istep; - } - imin = 2*imin-i1; - istep *= 2; - - L330: - if (istep <= ntwo) - goto L310; - - /*special case-- W = -sqrt(-1)*/ - imin = 3*np1+i1; - istep = 4*np1; - goto L420; - - L400: - j = imin-istep/2; - for (i = imin; i <= ntot; i += istep) { - if (ifrwd != 0) { - tempr = datai[i]; - tempi = -datar[i]; - } else { - tempr = -datai[i]; - tempi = datar[i]; - } - datar[i] = datar[j]-tempr; - datai[i] = datai[j]-tempi; - datar[j] += tempr; - datai[j] += tempi; - j += istep; - } - - imin = 2*imin-i1; - istep *= 2; - - L420: - if (istep <= ntwo) - goto L400; - } - - /*main loop for factors of 2. W=EXP(-2*PI*SQRT(-1)*m/mmax) */ - theta = -TWOPI/8.; - wstpr = 0.; - wstpi = -1.; - if (ifrwd == 0) { - theta = -theta; - wstpi = 1.; - } - mmax = 8*np1; - goto L540; - + i1rng = np1; + if (icase == 2) + i1rng = np0 * (1 + nprev / 2); + if (ntwo <= np1) + goto L600; + for (i1 = 1; i1 <= i1rng; i1++) { + imin = np1 + i1; + istep = 2 * np1; + goto L330; + + L310: + j = i1; + for (i = imin; i <= ntot; i += istep) { + tempr = datar[i]; + tempi = datai[i]; + datar[i] = datar[j] - tempr; + datai[i] = datai[j] - tempi; + datar[j] = datar[j] + tempr; + datai[j] = datai[j] + tempi; + j += istep; + } + imin = 2 * imin - i1; + istep *= 2; + + L330: + if (istep <= ntwo) + goto L310; + + /*special case-- W = -sqrt(-1)*/ + imin = 3 * np1 + i1; + istep = 4 * np1; + goto L420; + + L400: + j = imin - istep / 2; + for (i = imin; i <= ntot; i += istep) { + if (ifrwd != 0) { + tempr = datai[i]; + tempi = -datar[i]; + } else { + tempr = -datai[i]; + tempi = datar[i]; + } + datar[i] = datar[j] - tempr; + datai[i] = datai[j] - tempi; + datar[j] += tempr; + datai[j] += tempi; + j += istep; + } + + imin = 2 * imin - i1; + istep *= 2; + + L420: + if (istep <= ntwo) + goto L400; + } + + /*main loop for factors of 2. W=EXP(-2*PI*SQRT(-1)*m/mmax) */ + theta = -TWOPI / 8.; + wstpr = 0.; + wstpi = -1.; + if (ifrwd == 0) { + theta = -theta; + wstpi = 1.; + } + mmax = 8 * np1; + goto L540; + L500: - wminr = cos(theta); - wmini = sin(theta); - wr = wminr; - wi = wmini; - mmin = mmax/2+np1; - mstep = np1*2; - for (m = mmin; m <= mmax; m += mstep) { - for (i1 = 1; i1 <= i1rng; i1++) { - istep = mmax; - imin = m+i1; - L510: - j = imin-istep/2; - for (i = imin; i <= ntot; i += istep) { - tempr = datar[i]*wr-datai[i]*wi; - tempi = datar[i]*wi+datai[i]*wr; - datar[i] = datar[j]-tempr; - datai[i] = datai[j]-tempi; - datar[j] += tempr; - datai[j] += tempi; - j += istep; - } - imin = 2*imin-i1; - istep *= 2; - if (istep <= ntwo) - goto L510; - } - wtemp = wr*wstpi; - wr = wr*wstpr-wi*wstpi; - wi = wi*wstpr+wtemp; - } - wstpr = wminr; - wstpi = wmini; - theta /= 2.; - mmax += mmax; + wminr = cos(theta); + wmini = sin(theta); + wr = wminr; + wi = wmini; + mmin = mmax / 2 + np1; + mstep = np1 * 2; + for (m = mmin; m <= mmax; m += mstep) { + for (i1 = 1; i1 <= i1rng; i1++) { + istep = mmax; + imin = m + i1; + L510: + j = imin - istep / 2; + for (i = imin; i <= ntot; i += istep) { + tempr = datar[i] * wr - datai[i] * wi; + tempi = datar[i] * wi + datai[i] * wr; + datar[i] = datar[j] - tempr; + datai[i] = datai[j] - tempi; + datar[j] += tempr; + datai[j] += tempi; + j += istep; + } + imin = 2 * imin - i1; + istep *= 2; + if (istep <= ntwo) + goto L510; + } + wtemp = wr * wstpi; + wr = wr * wstpr - wi * wstpi; + wi = wi * wstpr + wtemp; + } + wstpr = wminr; + wstpi = wmini; + theta /= 2.; + mmax += mmax; L540: - if (mmax <= ntwo) - goto L500; - - /*main loop for factors not equal to 2-- W=EXP(-2*PI*SQRT(-1)*(j2-i3)/ifp2)*/ + if (mmax <= ntwo) + goto L500; + + /*main loop for factors not equal to 2-- W=EXP(-2*PI*SQRT(-1)*(j2-i3)/ifp2)*/ L600: - if (non2p <= 1) - goto L700; - ifp1 = ntwo; - iff = inon2; + if (non2p <= 1) + goto L700; + ifp1 = ntwo; + iff = inon2; L610: - ifp2 = ifact[iff]*ifp1; - theta = -TWOPI/ (double)ifact[iff]; - if (ifrwd == 0) - theta = -theta; - thetm = theta/ (double)(ifp1/np1); - wstpr = cos(theta); - wstpi = sin(theta); - wmstr = cos(thetm); - wmsti = sin(thetm); - wminr = 1.; - wmini = 0.; - - for (j1 = 1; j1 <= ifp1; j1 += np1) { - i1max = j1+i1rng-1; - for (i1 = j1; i1 <= i1max; i1++) { - for (i3 = i1; i3 <= ntot; i3 += np2) { - i = 1; - wr = wminr; - wi = wmini; - j2max = i3+ifp2-ifp1; - for (j2 = i3; j2 <= j2max; j2 += ifp1) { - twowr = 2.*wr; - jmin = i3; - j3max = j2+np2-ifp2; - for (j3 = j2 ; j3 <= j3max; j3 += ifp2) { - j = jmin+ifp2-ifp1; - sr = datar[j]; - si = datai[j]; - oldsr = 0.; - oldsi = 0.; - j -= ifp1; - L620: - stmpr = sr; - stmpi = si; - sr = twowr*sr-oldsr+datar[j]; - si = twowr*si-oldsi+datai[j]; - oldsr = stmpr; - oldsi = stmpi; - j -= ifp1; - if (j > jmin) - goto L620; - workr[i] = wr*sr-wi*si-oldsr+datar[j]; - worki[i] = wi*sr+wr*si-oldsi+datai[j]; - jmin += ifp2; - i++; - } - wtemp = wr*wstpi; - wr = wr*wstpr-wi*wstpi; - wi = wi*wstpr+wtemp; - } - i = 1; - for (j2 = i3; j2 <= j2max; j2 += ifp1) { - j3max = j2+np2-ifp2; - for (j3 = j2; j3 <= j3max; j3 += ifp2) { - datar[j3] = workr[i]; - datai[j3] = worki[i]; - i++; - } - } - } - } - wtemp = wminr*wmsti; - wminr = wminr*wmstr-wmini*wmsti; - wmini = wmini*wmstr+wtemp; - } - iff++; - ifp1 = ifp2; - if (ifp1 < np2) - goto L610; - - /*complete a real transform in the 1st dimension, n even, by conjugate symmetries*/ + ifp2 = ifact[iff] * ifp1; + theta = -TWOPI / (double)ifact[iff]; + if (ifrwd == 0) + theta = -theta; + thetm = theta / (double)(ifp1 / np1); + wstpr = cos(theta); + wstpi = sin(theta); + wmstr = cos(thetm); + wmsti = sin(thetm); + wminr = 1.; + wmini = 0.; + + for (j1 = 1; j1 <= ifp1; j1 += np1) { + i1max = j1 + i1rng - 1; + for (i1 = j1; i1 <= i1max; i1++) { + for (i3 = i1; i3 <= ntot; i3 += np2) { + i = 1; + wr = wminr; + wi = wmini; + j2max = i3 + ifp2 - ifp1; + for (j2 = i3; j2 <= j2max; j2 += ifp1) { + twowr = 2. * wr; + jmin = i3; + j3max = j2 + np2 - ifp2; + for (j3 = j2; j3 <= j3max; j3 += ifp2) { + j = jmin + ifp2 - ifp1; + sr = datar[j]; + si = datai[j]; + oldsr = 0.; + oldsi = 0.; + j -= ifp1; + L620: + stmpr = sr; + stmpi = si; + sr = twowr * sr - oldsr + datar[j]; + si = twowr * si - oldsi + datai[j]; + oldsr = stmpr; + oldsi = stmpi; + j -= ifp1; + if (j > jmin) + goto L620; + workr[i] = wr * sr - wi * si - oldsr + datar[j]; + worki[i] = wi * sr + wr * si - oldsi + datai[j]; + jmin += ifp2; + i++; + } + wtemp = wr * wstpi; + wr = wr * wstpr - wi * wstpi; + wi = wi * wstpr + wtemp; + } + i = 1; + for (j2 = i3; j2 <= j2max; j2 += ifp1) { + j3max = j2 + np2 - ifp2; + for (j3 = j2; j3 <= j3max; j3 += ifp2) { + datar[j3] = workr[i]; + datai[j3] = worki[i]; + i++; + } + } + } + } + wtemp = wminr * wmsti; + wminr = wminr * wmstr - wmini * wmsti; + wmini = wmini * wmstr + wtemp; + } + iff++; + ifp1 = ifp2; + if (ifp1 < np2) + goto L610; + + /*complete a real transform in the 1st dimension, n even, by conjugate symmetries*/ L700: - switch (icase) { - case 1: - goto L900; - break; - case 2: - goto L800; - break; - case 3: - goto L900; - break; - } - - nhalf = n; - n += n; - theta = -TWOPI/ (double) n; - if (ifrwd == 0) - theta = -theta; - wstpr = cos(theta); - wstpi = sin(theta); - wr = wstpr; - wi = wstpi; - imin = 2; - jmin = nhalf; - goto L725; + switch (icase) { + case 1: + goto L900; + break; + case 2: + goto L800; + break; + case 3: + goto L900; + break; + } + + nhalf = n; + n += n; + theta = -TWOPI / (double)n; + if (ifrwd == 0) + theta = -theta; + wstpr = cos(theta); + wstpi = sin(theta); + wr = wstpr; + wi = wstpi; + imin = 2; + jmin = nhalf; + goto L725; L710: - j = jmin; - for (i = imin; i <= ntot; i += np2) { - sumr = (datar[i]+datar[j])/2.; - sumi = (datai[i]+datai[j])/2.; - difr = (datar[i]-datar[j])/2.; - difi = (datai[i]-datai[j])/2.; - tempr = wr*sumi+wi*difr; - tempi = wi*sumi-wr*difr; - datar[i] = sumr+tempr; - datai[i] = difi+tempi; - datar[j] = sumr-tempr; - datai[j] = tempi-difi; - j += np2; - } - imin++; - jmin--; - wtemp = wr*wstpi; - wr = wr*wstpr-wi*wstpi; - wi = wi*wstpr+wtemp; + j = jmin; + for (i = imin; i <= ntot; i += np2) { + sumr = (datar[i] + datar[j]) / 2.; + sumi = (datai[i] + datai[j]) / 2.; + difr = (datar[i] - datar[j]) / 2.; + difi = (datai[i] - datai[j]) / 2.; + tempr = wr * sumi + wi * difr; + tempi = wi * sumi - wr * difr; + datar[i] = sumr + tempr; + datai[i] = difi + tempi; + datar[j] = sumr - tempr; + datai[j] = tempi - difi; + j += np2; + } + imin++; + jmin--; + wtemp = wr * wstpi; + wr = wr * wstpr - wi * wstpi; + wi = wi * wstpr + wtemp; L725: - if (imin < jmin) { - goto L710; - } else if (imin > jmin) { - goto L740; - } - if (ifrwd == 0) - goto L740; - for (i = imin; i <= ntot; i += np2) { - datai[i] = -datai[i]; - } + if (imin < jmin) { + goto L710; + } else if (imin > jmin) { + goto L740; + } + if (ifrwd == 0) + goto L740; + for (i = imin; i <= ntot; i += np2) { + datai[i] = -datai[i]; + } L740: - np2 *= 2; - ntot *= 2; - j = ntot+1; - imax = ntot/2+1; + np2 *= 2; + ntot *= 2; + j = ntot + 1; + imax = ntot / 2 + 1; L745: - imin = imax-nhalf; - i = imin; - goto L755; + imin = imax - nhalf; + i = imin; + goto L755; L750: - datar[j] = datar[i]; - datai[j] = -datai[i]; + datar[j] = datar[i]; + datai[j] = -datai[i]; L755: - i++; - j--; - if (i < imax) - goto L750; - datar[j] = datar[imin]-datai[imin]; - datai[j] = 0.; - if (i >= j) { - goto L780; - } else { - goto L770; - } + i++; + j--; + if (i < imax) + goto L750; + datar[j] = datar[imin] - datai[imin]; + datai[j] = 0.; + if (i >= j) { + goto L780; + } else { + goto L770; + } L765: - datar[j] = datar[i]; - datai[j] = datai[i]; + datar[j] = datar[i]; + datai[j] = datai[i]; L770: - i--; - j--; - if (i > imin) - goto L765; - datar[j] = datar[imin]+datai[imin]; - datai[j] = 0.; - imax = imin; - goto L745; + i--; + j--; + if (i > imin) + goto L765; + datar[j] = datar[imin] + datai[imin]; + datai[j] = 0.; + imax = imin; + goto L745; L780: - datar[1] += datai[1]; - datai[1] = 0.; - goto L900; - - /*complete a real transform for the 2nd, 3rd, ... dimension by conjugate symmetries*/ + datar[1] += datai[1]; + datai[1] = 0.; + goto L900; + + /*complete a real transform for the 2nd, 3rd, ... dimension by conjugate symmetries*/ L800: - if (nprev <= 2) - goto L900; - for (i3 = 1; i3 <= ntot; i3 += np2) { - i2max = i3+np2-np1; - for (i2 = i3; i2 <= i2max; i2 += np1) { - imax = i2+np1-1; - imin = i2+i1rng; - jmax = 2*i3+np1-imin; - if (i2 > i3) - jmax += np2; - if (idim > 2) { - j = jmax+np0; - for (i = imin; i <= imax; i++) { - datar[i] = datar[j]; - datai[i] = -datai[j]; - j--; - } - } - j = jmax; - for (i = imin; i <= imax; i += np0) { - datar[i] = datar[j]; - datai[i] = -datai[j]; - j -= np0; - } - } - } - - /*end of loop on each dimension*/ + if (nprev <= 2) + goto L900; + for (i3 = 1; i3 <= ntot; i3 += np2) { + i2max = i3 + np2 - np1; + for (i2 = i3; i2 <= i2max; i2 += np1) { + imax = i2 + np1 - 1; + imin = i2 + i1rng; + jmax = 2 * i3 + np1 - imin; + if (i2 > i3) + jmax += np2; + if (idim > 2) { + j = jmax + np0; + for (i = imin; i <= imax; i++) { + datar[i] = datar[j]; + datai[i] = -datai[j]; + j--; + } + } + j = jmax; + for (i = imin; i <= imax; i += np0) { + datar[i] = datar[j]; + datai[i] = -datai[j]; + j -= np0; + } + } + } + + /*end of loop on each dimension*/ L900: - np0 = np1; - np1 = np2; - nprev = n; + np0 = np1; + np1 = np2; + nprev = n; } - L920: +L920: return; } diff --git a/fftma_module/gen/lib_src/gammf.c b/fftma_module/gen/lib_src/gammf.c old mode 100755 new mode 100644 index 97ff3f9..c2c31cf --- a/fftma_module/gen/lib_src/gammf.c +++ b/fftma_module/gen/lib_src/gammf.c @@ -1,16 +1,14 @@ -#include -#include #include "genlib.h" - +#include +#include /*gamma covariance function*/ double gammf(double h, double alpha) { - float delta; - double z; + float delta; + double z; - delta = pow(20.,1./alpha)-1.; - z = 1./(double)(pow(1.+h*delta,alpha)); - return(z); + delta = pow(20., 1. / alpha) - 1.; + z = 1. / (double)(pow(1. + h * delta, alpha)); + return (z); } - diff --git a/fftma_module/gen/lib_src/gasdev.c b/fftma_module/gen/lib_src/gasdev.c old mode 100755 new mode 100644 index a25dfa2..c66f9c9 --- a/fftma_module/gen/lib_src/gasdev.c +++ b/fftma_module/gen/lib_src/gasdev.c @@ -1,31 +1,31 @@ -#include #include "genlib.h" +#include #define NTAB 32 -double gasdev(long *idum,long *idum2, long *iy, long iv[NTAB]) +double gasdev(long* idum, long* idum2, long* iy, long iv[NTAB]) /*returns a normally distributed deviate with 0 mean*/ /*and unit variance, using ran2(idum) as the source */ /*of uniform deviates */ { - double ran2(long *idum,long *idum2, long *iy, long iv[NTAB]); - static int iset = 0; - static double gset; - double fac,rsq,v1,v2; + double ran2(long* idum, long* idum2, long* iy, long iv[NTAB]); + static int iset = 0; + static double gset; + double fac, rsq, v1, v2; if (iset == 0) { - do { - v1 = 2.0*ran2(idum,idum2,iy,iv)-1.0; - v2 = 2.0*ran2(idum,idum2,iy,iv)-1.0; - rsq = v1*v1+v2*v2; - } while (rsq >= 1.0 || rsq == 0.0); + do { + v1 = 2.0 * ran2(idum, idum2, iy, iv) - 1.0; + v2 = 2.0 * ran2(idum, idum2, iy, iv) - 1.0; + rsq = v1 * v1 + v2 * v2; + } while (rsq >= 1.0 || rsq == 0.0); - fac = sqrt(-2.0*log(rsq)/rsq); - gset = v1*fac; - iset = 1; - return (v2*fac); - } else { - iset = 0; - return (gset); - } + fac = sqrt(-2.0 * log(rsq) / rsq); + gset = v1 * fac; + iset = 1; + return (v2 * fac); + } else { + iset = 0; + return (gset); + } } diff --git a/fftma_module/gen/lib_src/gaussian.c b/fftma_module/gen/lib_src/gaussian.c old mode 100755 new mode 100644 index 5ef9911..a8f85f2 --- a/fftma_module/gen/lib_src/gaussian.c +++ b/fftma_module/gen/lib_src/gaussian.c @@ -1,10 +1,9 @@ -#include -#include #include "genlib.h" - +#include +#include /*gaussian covariance function*/ double gaussian(double h) { - return (exp(-3.*(double)(h*h))); + return (exp(-3. * (double)(h * h))); } diff --git a/fftma_module/gen/lib_src/generate.c b/fftma_module/gen/lib_src/generate.c old mode 100755 new mode 100644 index 46869d3..8f8f06a --- a/fftma_module/gen/lib_src/generate.c +++ b/fftma_module/gen/lib_src/generate.c @@ -1,43 +1,39 @@ -#include -#include #include "geostat.h" - +#include +#include /* GENERATION OF A GAUSSIAN WHITE NOISE VECTOR */ /*input: */ /* seed: seed */ /* n: number of components in the vector */ -/*output: */ +/*output: */ /* realization: structure defining the realization*/ -void generate(long *seed, int n, struct realization_mod *realization) +void generate(long* seed, int n, struct realization_mod* realization) { - int i; - long idum2 = 123456789,iy = 0; - long *iv; - int iset =0; - - iv = (long *) malloc(NTAB * sizeof(long)); - - /*negative seed*/ - if (*seed > 0.0) - *seed = -(*seed); - - /*realization definition*/ - (*realization).n = n; - (*realization).code = 0; - (*realization).vector = (double *) malloc(n * sizeof(double)); - if ((*realization).vector == NULL) { - printf("No memory available in generate"); - exit; - } - - /*Gaussian white noise generation*/ - for (i=0; i < n; i++) - (*realization).vector[i] = gasdev(seed,&idum2,&iy,iv,&iset); - - return; + int i; + long idum2 = 123456789, iy = 0; + long* iv; + int iset = 0; + + iv = (long*)malloc(NTAB * sizeof(long)); + + /*negative seed*/ + if (*seed > 0.0) + *seed = -(*seed); + + /*realization definition*/ + (*realization).n = n; + (*realization).code = 0; + (*realization).vector = (double*)malloc(n * sizeof(double)); + if ((*realization).vector == NULL) { + printf("No memory available in generate"); + exit; + } + + /*Gaussian white noise generation*/ + for (i = 0; i < n; i++) + (*realization).vector[i] = gasdev(seed, &idum2, &iy, iv, &iset); + + return; } - - - diff --git a/fftma_module/gen/lib_src/genlib.c b/fftma_module/gen/lib_src/genlib.c old mode 100755 new mode 100644 index 1a106f6..8575537 --- a/fftma_module/gen/lib_src/genlib.c +++ b/fftma_module/gen/lib_src/genlib.c @@ -7,10 +7,10 @@ * interface description in genlib.h for details. */ -#include +#include #include +#include #include -#include #include "genlib.h" @@ -34,16 +34,17 @@ char undefined_object[] = "UNDEFINED"; /* Section 2 -- Memory allocation */ -void *GetBlock(size_t nbytes) +void* GetBlock(size_t nbytes) { - void *result; + void* result; result = malloc(nbytes); - if (result == NULL) Error("No memory available"); + if (result == NULL) + Error("No memory available"); return (result); } -void FreeBlock(void *ptr) +void FreeBlock(void* ptr) { free(ptr); } diff --git a/fftma_module/gen/lib_src/genlib.h b/fftma_module/gen/lib_src/genlib.h old mode 100755 new mode 100644 index f5a490f..3f1fb2f --- a/fftma_module/gen/lib_src/genlib.h +++ b/fftma_module/gen/lib_src/genlib.h @@ -25,9 +25,9 @@ #ifndef _genlib_h #define _genlib_h +#include #include #include -#include /* Section 1 -- Define new "primitive" types */ @@ -45,20 +45,21 @@ */ #ifdef THINK_C - typedef int bool; +typedef int bool; +#else +#ifdef TRUE +#ifndef bool +#define bool int +#endif #else -# ifdef TRUE -# ifndef bool -# define bool int -# endif -# else -# ifdef bool -# define FALSE 0 -# define TRUE 1 -# else - typedef enum {FALSE, TRUE} bool; -# endif -# endif +#ifdef bool +#define FALSE 0 +#define TRUE 1 +#else +typedef enum { FALSE, + TRUE } bool; +#endif +#endif #endif /* @@ -72,7 +73,7 @@ * Declaring it as a string emphasizes this atomicity. */ -typedef char *string; +typedef char* string; /* * Constant: UNDEFINED @@ -84,7 +85,7 @@ typedef char *string; * therefore inappropriate as a sentinel. */ -#define UNDEFINED ((void *) undefined_object) +#define UNDEFINED ((void*)undefined_object) extern char undefined_object[]; @@ -111,7 +112,7 @@ extern char undefined_object[]; * no memory is available, GetBlock generates an error. */ -void *GetBlock(size_t nbytes); +void* GetBlock(size_t nbytes); /* * Function: FreeBlock @@ -121,7 +122,7 @@ void *GetBlock(size_t nbytes); * have been allocated using GetBlock, New, or NewArray. */ -void FreeBlock(void *ptr); +void FreeBlock(void* ptr); /* * Macro: New @@ -135,7 +136,7 @@ void FreeBlock(void *ptr); * target type. */ -#define New(type) ((type) GetBlock(sizeof *((type) NULL))) +#define New(type) ((type)GetBlock(sizeof *((type)NULL))) /* * Macro: NewArray @@ -145,7 +146,7 @@ void FreeBlock(void *ptr); * values of the specified element type. */ -#define NewArray(n, type) ((type *) GetBlock((n) * sizeof(type))) +#define NewArray(n, type) ((type*)GetBlock((n) * sizeof(type))) /* Section 3 -- Basic error handling */ diff --git a/fftma_module/gen/lib_src/geostat.h b/fftma_module/gen/lib_src/geostat.h old mode 100755 new mode 100644 index f5aeb1a..78b121a --- a/fftma_module/gen/lib_src/geostat.h +++ b/fftma_module/gen/lib_src/geostat.h @@ -1,8 +1,7 @@ -#include +#include #include +#include #include -#include - #ifndef _GEOSTAT_H #define _GEOSTAT_H @@ -27,7 +26,6 @@ /* cdf_mod */ /* realization_mod */ - /* List of functions: */ /* ------------------ */ @@ -36,7 +34,7 @@ /* cubic, cutspectr, deflimit, dual_kri */ /* exponential, fourt, funtrun1, G, gammf */ /* gammln, gammp, gasdev, gaussian, gcf, */ -/* gen_cov_matrix, Ginv, gradual, cgrid, */ +/* gen_cov_matrix, Ginv, gradual, cgrid, */ /* gser, invtrun1, krig_stat, length, */ /* LtimeZ, mat_vec, maxfactor, metrop, norm */ /* normal, nugget, power, ran2, scal_vect, */ @@ -44,9 +42,6 @@ /* test_fract, trun1, trungasdev,vec_vec, */ /* vf2gthres,polint */ - - - /*STRUCTURES*/ /*----------*/ /*variogram */ @@ -67,16 +62,14 @@ /*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; + 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 */ @@ -95,18 +88,15 @@ struct vario_mod { /*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; + 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*/ @@ -118,12 +108,11 @@ struct variotable_mod { /*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; + 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 */ @@ -153,18 +142,16 @@ struct grid_mod { /* 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; + 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 */ @@ -174,12 +161,11 @@ struct welldata_mod { /*vf: volume fractions for the first ncat-1 facies*/ /* i = [0...ncat-2]*nblock */ struct statfacies_mod { - int ncat; - int nblock; - float *vf; + 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 */ @@ -202,14 +188,13 @@ struct statfacies_mod { /* 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; + int nsY1; + int nsY2; + float* thresholds; + int* address_sY1; + int* address_sY2; }; - /*statistical data */ /*type --> 0 : normal */ /* --> 1 : natural log */ @@ -227,14 +212,13 @@ struct inequalities_mod { /*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; + int type; + int nblock_mean; + double* mean; + int nblock_var; + double* variance; }; - /*gradual deformation parameters */ /*Nadded: number of complementary realizations */ /*NZONES: number of subregions */ @@ -245,12 +229,11 @@ struct statistic_mod { /*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; + int Nadded, NZONES; + float* rho; + int *cellini, *cellfin; }; - /*gradient structures */ /*Nparam : number of parameters for which gradients are */ /* required */ @@ -263,24 +246,21 @@ struct grad_mod { /* Ncells....2*Ncells-1 for the second parameter */ /* and so on */ struct gradients_mod { - int Nparam,Ncells; - float *grad; + 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; + int n; + float* x; + float* fx; }; - /*realization */ /*n: number of components */ /*code: status of the realization */ @@ -301,9 +281,9 @@ struct cdf_mod { /* --> 14: conditional any type */ /*vector: realization vector i = [0...n-1] */ struct realization_mod { - int n; - int code; - double *vector; + int n; + int code; + double* vector; }; /*=====================================================*/ @@ -311,18 +291,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); - +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 */ @@ -330,8 +307,7 @@ double cardsin(double h); /* n : dimension of matrix Cij */ /* */ /* C is turned into the lower triangular cholesky matrix*/ -void choldc(double *C, int n); - +void choldc(double* C, int n); /*computes the coordinates of a given cell */ /*as numbers of cells along the X,Y and Z axes*/ @@ -341,8 +317,7 @@ void choldc(double *C, int n); /*grid: structure defining the grid */ /*output: */ /*i: vector with the coordinates */ -void coordinates(int maille, int i[3],struct grid_mod grid); - +void coordinates(int maille, int i[3], struct grid_mod grid); /*builds the sampled covariance function */ /*dimensions are even */ @@ -351,7 +326,7 @@ void coordinates(int maille, int i[3],struct grid_mod grid); /*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]); +void covariance(double* covar, struct vario_mod variogram, struct grid_mod grid, int n[3]); /*computation of the covariance matrix for the well data*/ /*well coordinates are given as a number of cells */ @@ -362,8 +337,7 @@ void covariance(double *covar,struct vario_mod variogram, struct grid_mod grid, /*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); - +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 */ @@ -382,12 +356,11 @@ void cov_matrix(double *C, struct vario_mod variogram, struct welldata_mod well, /*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); +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 */ @@ -395,8 +368,7 @@ double cubic(double h); /*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]); - +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 */ @@ -405,8 +377,7 @@ void cutspectr(float *covar, int kx, int ky, int kz, int n[3]); /*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); - +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 */ @@ -416,14 +387,11 @@ void deflimit(double *plim_inf, double *plim_sup, float x, float *thresholds, st /*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); - - +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 */ @@ -434,44 +402,36 @@ double exponential(double h); /*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); - +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 */ /*idum: seed */ -double gasdev(long *idum, long *idum2, long *iy, long *iv, int *iset); +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))*/ +/*fraction represented as gammcf, also returns ln(G(a))*/ /*as gln */ -void gcf(float *gammcf, float a, float x, float *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 */ @@ -482,19 +442,17 @@ void gcf(float *gammcf, float a, float x, float *gln); /*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); - +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 */ +/*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 */ @@ -510,8 +468,7 @@ float Ginv(float p); /*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); - +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: */ @@ -524,10 +481,9 @@ void gradual(struct grad_mod grad,float *Zo,float *Z,float *Zfinal,int n,struct void cgrid(struct vario_mod variogram, struct grid_mod grid, int n[3]); /*incomplete gamma function evaluated by its series*/ -/*representation as gamser, also returns ln(G(a)) */ +/*representation as gamser, also returns ln(G(a)) */ /*as gln */ -void gser(float *gamser, float a, float x, float *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 */ @@ -538,18 +494,16 @@ 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); - +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 length(int N, int i, double* scf, double* ap, double D, int Nvari); /*calculates L.Z/ /* L : lower triangular matrix recorded */ @@ -560,7 +514,7 @@ int length(int N, int i, double *scf, double *ap, double D, int Nvari); /* n : dimension of matrix Lij */ /* */ /* The solution vector is returned in b */ -void LtimeZ(double *L, float *Z, float *b, int n); +void LtimeZ(double* L, float* Z, float* b, int n); /*determines the greatest prime factor of an integer*/ int maxfactor(int n); @@ -570,40 +524,32 @@ int maxfactor(int n); /*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); - +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); - +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); - +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); - - +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); - +double scal_vec(double* b, int n); /*solves the set of n linear equations Cx = D */ /* C : symmetric positive-definite matrix recorded */ @@ -614,32 +560,26 @@ double scal_vec(double *b,int n); /* */ /* The solution vector is returned in D */ /* CONJUGATE GRADIENT method */ -void solve3(double *C, double *D, int n); - +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); - +/*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); - +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); - +int test_fact(int* pnum, int fact, int* pmaxfac); /*calculates the integrale of an approximate function*/ /*for the Gaussian function over an interval defined */ @@ -648,29 +588,25 @@ 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 */ -double trungasdev(long *idum,double lim_inf,double lim_sup,long *idum2, long *iy, long iv[NTAB]); +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); - +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); - +void vf2gthres(struct statfacies_mod facies, float* thresholds); +void polint(float xa[], float ya[], int n, float x, float* y, float* dy); #endif diff --git a/fftma_module/gen/lib_src/length.c b/fftma_module/gen/lib_src/length.c old mode 100755 new mode 100644 index c2dfc8b..7534a35 --- a/fftma_module/gen/lib_src/length.c +++ b/fftma_module/gen/lib_src/length.c @@ -1,39 +1,39 @@ #include /* compute the length for one dimension*/ -int length(int N, int i, double *scf, double *ap, double D, int Nvari) +int length(int N, int i, double* scf, double* ap, double D, int Nvari) { - int maxfactor(int n); - double temp1, temp2; - int n, j, k, nmax; - int nlimit = 13; - - if (N == 1) { - n = 1; + int maxfactor(int n); + double temp1, temp2; + int n, j, k, nmax; + int nlimit = 13; + + if (N == 1) { + n = 1; } else { - - for (k = 0; k < Nvari; k++) { - temp1 = fabs(ap[9*k+i])*scf[3*k]*scf[3*k]; - for (j = 1; j <3; j++) { - temp2 = fabs(ap[9*k+3*j+i])*scf[3*k+j]*scf[3*k+j]; - if (temp2 > temp1) - temp1 = temp2; - } - } - temp1 = sqrt(temp1); - temp1 /= (double)D; - if ((double)N/temp1 < 2.) { - n = N+(int)(2*temp1); - } else { - n = N+(int)temp1; - } - if ((n % 2) != 0) - n = n+1; - nmax = maxfactor(n); - while (nmax > nlimit) { - n += 2; - nmax = maxfactor(n); - } + + for (k = 0; k < Nvari; k++) { + temp1 = fabs(ap[9 * k + i]) * scf[3 * k] * scf[3 * k]; + for (j = 1; j < 3; j++) { + temp2 = fabs(ap[9 * k + 3 * j + i]) * scf[3 * k + j] * scf[3 * k + j]; + if (temp2 > temp1) + temp1 = temp2; + } + } + temp1 = sqrt(temp1); + temp1 /= (double)D; + if ((double)N / temp1 < 2.) { + n = N + (int)(2 * temp1); + } else { + n = N + (int)temp1; + } + if ((n % 2) != 0) + n = n + 1; + nmax = maxfactor(n); + while (nmax > nlimit) { + n += 2; + nmax = maxfactor(n); + } } - return (n); + return (n); } diff --git a/fftma_module/gen/lib_src/maxfactor.c b/fftma_module/gen/lib_src/maxfactor.c old mode 100755 new mode 100644 index 842cd8b..7658d93 --- a/fftma_module/gen/lib_src/maxfactor.c +++ b/fftma_module/gen/lib_src/maxfactor.c @@ -3,33 +3,33 @@ /*determines the greatest prime factor of an integer*/ int maxfactor(int n) { - int test_fact(int *pnum, int fact, int *pmaxfac); - int lnum, fact; - int maxfac; + int test_fact(int* pnum, int fact, int* pmaxfac); + int lnum, fact; + int maxfac; - maxfac = 1; - lnum = n; + maxfac = 1; + lnum = n; - if ( lnum != 0 && lnum != 1 ) { - fact = 2; - if ( test_fact( &lnum, fact,&maxfac)) { - fact = 3; - if ( test_fact( &lnum, fact,&maxfac)) { - fact = 5; - for (;;) { - if (!test_fact( &lnum, fact,&maxfac)) - break; - fact += 2; - if (!test_fact( &lnum, fact,&maxfac)) - break; - fact += 4; - } - } + if (lnum != 0 && lnum != 1) { + fact = 2; + if (test_fact(&lnum, fact, &maxfac)) { + fact = 3; + if (test_fact(&lnum, fact, &maxfac)) { + fact = 5; + for (;;) { + if (!test_fact(&lnum, fact, &maxfac)) + break; + fact += 2; + if (!test_fact(&lnum, fact, &maxfac)) + break; + fact += 4; + } + } + } + if (lnum != 1) { + if (lnum > maxfac) + maxfac = lnum; + } } - if ( lnum != 1 ) { - if (lnum > maxfac) - maxfac = lnum; - } - } - return (maxfac); + return (maxfac); } diff --git a/fftma_module/gen/lib_src/nor2log.c b/fftma_module/gen/lib_src/nor2log.c old mode 100755 new mode 100644 index 44b6568..72dd830 --- a/fftma_module/gen/lib_src/nor2log.c +++ b/fftma_module/gen/lib_src/nor2log.c @@ -1,7 +1,6 @@ -#include -#include #include "geostat.h" - +#include +#include /*TURNS NORMAL NUMBERS INTO LOGNORMAL NUMBERS */ /*input: */ @@ -13,66 +12,62 @@ /*realout: structure defining a realization - */ /* lognormal numbers */ -void nor2log(struct realization_mod *realin, int typelog, struct realization_mod *realout) +void nor2log(struct realization_mod* realin, int typelog, struct realization_mod* realout) { - int i; - double coeff; - - - coeff = log(10.0); + int i; + double coeff; + coeff = log(10.0); - /*Is the output realization allocated ?*/ - /*is its length equal to realin.n?*/ - if ((*realout).vector == NULL || (*realout).n != (*realin).n) { - (*realout).vector = (double *) malloc((*realin).n * sizeof(double)); - if ((*realout).vector == NULL) { - printf("No memory available"); - return; + /*Is the output realization allocated ?*/ + /*is its length equal to realin.n?*/ + if ((*realout).vector == NULL || (*realout).n != (*realin).n) { + (*realout).vector = (double*)malloc((*realin).n * sizeof(double)); + if ((*realout).vector == NULL) { + printf("No memory available"); + return; + } } - } - (*realout).n = (*realin).n; + (*realout).n = (*realin).n; - - switch ((*realin).code) { - case 0: - case 1: - case 2: - (*realout).code = typelog; - break; - case 6: - case 7: - if (typelog == 3) { - (*realout).code = 8; - } else if (typelog == 4) { - (*realout).code = 9; + switch ((*realin).code) { + case 0: + case 1: + case 2: + (*realout).code = typelog; + break; + case 6: + case 7: + if (typelog == 3) { + (*realout).code = 8; + } else if (typelog == 4) { + (*realout).code = 9; + } + break; + default: + printf("Input is not normal in nor2log"); + return; + break; } - break; - default: - printf("Input is not normal in nor2log"); - return; - break; - } - - /*anamorphose*/ - for (i = 0; i < (*realin).n; i++) { - switch (typelog) { - case 3: - /*natural logarithm*/ - (*realout).vector[i] = exp((*realin).vector[i]); - break; - case 4: - /*log10*/ - (*realout).vector[i] = exp((*realin).vector[i]*coeff); - break; - default: - printf("Unexpected case in nor2log"); - return; - break; + /*anamorphose*/ + for (i = 0; i < (*realin).n; i++) { + switch (typelog) { + case 3: + /*natural logarithm*/ + (*realout).vector[i] = exp((*realin).vector[i]); + break; + case 4: + /*log10*/ + (*realout).vector[i] = exp((*realin).vector[i] * coeff); + break; + default: + printf("Unexpected case in nor2log"); + return; + break; + } } - } - return; + return; } diff --git a/fftma_module/gen/lib_src/nugget.c b/fftma_module/gen/lib_src/nugget.c old mode 100755 new mode 100644 index cda4929..31e3e4c --- a/fftma_module/gen/lib_src/nugget.c +++ b/fftma_module/gen/lib_src/nugget.c @@ -1,14 +1,13 @@ -#include -#include #include "genlib.h" - +#include +#include /*nugget covariance function*/ double nugget(double h) { - if (h == 0) { - return (1.); - } else { - return(0.); - } + if (h == 0) { + return (1.); + } else { + return (0.); + } } diff --git a/fftma_module/gen/lib_src/power.c b/fftma_module/gen/lib_src/power.c old mode 100755 new mode 100644 index 88281f6..1026213 --- a/fftma_module/gen/lib_src/power.c +++ b/fftma_module/gen/lib_src/power.c @@ -1,10 +1,9 @@ -#include -#include #include "genlib.h" - +#include +#include /*power covariance function*/ double power(double h, double alpha) { - return(pow(h,alpha)); + return (pow(h, alpha)); } diff --git a/fftma_module/gen/lib_src/prebuild_gwn.c b/fftma_module/gen/lib_src/prebuild_gwn.c old mode 100755 new mode 100644 index 447e9e7..b7f03e3 --- a/fftma_module/gen/lib_src/prebuild_gwn.c +++ b/fftma_module/gen/lib_src/prebuild_gwn.c @@ -1,10 +1,10 @@ -#include -#include -#include +#include "geostat.h" +#include #include +#include +#include #include -#include -#include "geostat.h" +#include /* prebuild_gwn */ /* Produce a first construction in real space of the Gaussian white noise */ @@ -19,36 +19,32 @@ /* must be a Gaussian white noise */ /*realization: structure defining a realization*/ -void prebuild_gwn(struct grid_mod grid,int n[3],struct realization_mod *realin,double *realization,int solver) +void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, double* realization, int solver) - { - int i,j,k,maille0,maille1; +{ + int i, j, k, maille0, maille1; int ntot; - ntot=n[0]*n[1]*n[2]; - realization[0]=0.; - if (solver==1) - { - for (i=0;i -#include #include "genlib.h" #include "geostat.h" +#include +#include #ifndef _PRESSURE_H #define _PRESSURE_H @@ -14,8 +14,7 @@ /* y: macroscopic gradient value in x direction */ /* z: macroscopic gradient value in x direction */ struct pressure_mod { - double x,y,z; + double x, y, z; }; - #endif // define _PRESSURE_H diff --git a/fftma_module/gen/lib_src/ran2.c b/fftma_module/gen/lib_src/ran2.c old mode 100755 new mode 100644 index 8075f4f..db0750e --- a/fftma_module/gen/lib_src/ran2.c +++ b/fftma_module/gen/lib_src/ran2.c @@ -2,8 +2,8 @@ #define IM1 2147483563 #define IM2 2147483399 -#define AM (1.0/IM1) -#define IMM1 (IM1-1) +#define AM (1.0 / IM1) +#define IMM1 (IM1 - 1) #define IA1 40014 #define IA2 40692 #define IQ1 53668 @@ -11,40 +11,49 @@ #define IR1 12211 #define IR2 3791 #define NTAB 32 -#define NDIV (1+IMM1/NTAB) +#define NDIV (1 + IMM1 / NTAB) #define EPS 1.2e-7 -#define RNMX (1.0-EPS) +#define RNMX (1.0 - EPS) -double ran2(long *idum, long *idum2, long *iy, long iv[NTAB]) +double ran2(long* idum, long* idum2, long* iy, long iv[NTAB]) { - int j; - long k; - double temp; + int j; + long k; + double temp; - if (*idum <= 0) { - if (-(*idum) < 1) *idum = 1; - else *idum = -(*idum); - *idum2 = (*idum); - for (j = NTAB+7; j >= 0; j--) { - k = (*idum)/IQ1; - *idum = IA1*(*idum-k*IQ1)-k*IR1; - if (*idum < 0) *idum += IM1; - if (j < NTAB) iv[j] = *idum; + if (*idum <= 0) { + if (-(*idum) < 1) + *idum = 1; + else + *idum = -(*idum); + *idum2 = (*idum); + for (j = NTAB + 7; j >= 0; j--) { + k = (*idum) / IQ1; + *idum = IA1 * (*idum - k * IQ1) - k * IR1; + if (*idum < 0) + *idum += IM1; + if (j < NTAB) + iv[j] = *idum; + } + *iy = iv[0]; } - *iy = iv[0]; - } - k = (*idum)/IQ1; - *idum = IA1*(*idum-k*IQ1)-k*IR1; - if (*idum < 0) *idum += IM1; - k = *idum2/IQ2; - *idum2 = IA2*(*idum2-k*IQ2)-k*IR2; - if (*idum2 < 0) *idum2 += IM2; - j = (*iy)/NDIV; - *iy = iv[j]-(*idum2); - iv[j] = *idum; - if (*iy < 1) (*iy) += IMM1; - if ((temp = AM*(*iy)) > RNMX) return (RNMX); - else return (temp); + k = (*idum) / IQ1; + *idum = IA1 * (*idum - k * IQ1) - k * IR1; + if (*idum < 0) + *idum += IM1; + k = *idum2 / IQ2; + *idum2 = IA2 * (*idum2 - k * IQ2) - k * IR2; + if (*idum2 < 0) + *idum2 += IM2; + j = (*iy) / NDIV; + *iy = iv[j] - (*idum2); + iv[j] = *idum; + if (*iy < 1) + (*iy) += IMM1; + if ((temp = AM * (*iy)) > RNMX) + return (RNMX); + else + return (temp); } diff --git a/fftma_module/gen/lib_src/random.c b/fftma_module/gen/lib_src/random.c old mode 100755 new mode 100644 index 696312d..b315f2f --- a/fftma_module/gen/lib_src/random.c +++ b/fftma_module/gen/lib_src/random.c @@ -25,7 +25,7 @@ void Randomize(void) { - srand((int) time(NULL)); + srand((int)time(NULL)); } /* @@ -44,8 +44,8 @@ int RandomInteger(int low, int high) int k; double d; - d = (double) rand() / ((double) RAND_MAX + 1); - k = (int) (d * (high - low + 1)); + d = (double)rand() / ((double)RAND_MAX + 1); + k = (int)(d * (high - low + 1)); return (low + k); } @@ -60,7 +60,7 @@ double RandomReal(double low, double high) { double d; - d = (double) rand() / ((double) RAND_MAX + 1); + d = (double)rand() / ((double)RAND_MAX + 1); return (low + d * (high - low)); } diff --git a/fftma_module/gen/lib_src/random.h b/fftma_module/gen/lib_src/random.h old mode 100755 new mode 100644 index 881af5f..167f890 --- a/fftma_module/gen/lib_src/random.h +++ b/fftma_module/gen/lib_src/random.h @@ -23,7 +23,7 @@ */ #ifndef RAND_MAX -# define RAND_MAX ((int) ((unsigned) ~0 >> 1)) +#define RAND_MAX ((int)((unsigned)~0 >> 1)) #endif /* diff --git a/fftma_module/gen/lib_src/scanadt.c b/fftma_module/gen/lib_src/scanadt.c old mode 100755 new mode 100644 index bd887c5..032a70a --- a/fftma_module/gen/lib_src/scanadt.c +++ b/fftma_module/gen/lib_src/scanadt.c @@ -4,11 +4,11 @@ * This file implements the scanadt.h interface. */ -#include -#include +#include "scanadt.h" #include "genlib.h" #include "strlib.h" -#include "scanadt.h" +#include +#include /* * Type: scannerCDT @@ -53,13 +53,15 @@ scannerADT NewScanner(void) void FreeScanner(scannerADT scanner) { - if (scanner->str != NULL) FreeBlock(scanner->str); + if (scanner->str != NULL) + FreeBlock(scanner->str); FreeBlock(scanner); } void SetScannerString(scannerADT scanner, string str) { - if (scanner->str != NULL) FreeBlock(scanner->str); + if (scanner->str != NULL) + FreeBlock(scanner->str); scanner->str = CopyString(str); scanner->len = StringLength(str); scanner->cp = 0; @@ -80,9 +82,11 @@ string ReadToken(scannerADT scanner) scanner->savedToken = NULL; return (token); } - if (scanner->spaceOption == IgnoreSpaces) SkipSpaces(scanner); + if (scanner->spaceOption == IgnoreSpaces) + SkipSpaces(scanner); start = finish = scanner->cp; - if (start >= scanner->len) return (CopyString("")); + if (start >= scanner->len) + return (CopyString("")); ch = scanner->str[scanner->cp]; if (isalnum(ch)) { finish = ScanToEndOfIdentifier(scanner); @@ -100,7 +104,8 @@ bool MoreTokensExist(scannerADT scanner) if (scanner->savedToken != NULL) { return (!StringEqual(scanner->savedToken, "")); } - if (scanner->spaceOption == IgnoreSpaces) SkipSpaces(scanner); + if (scanner->spaceOption == IgnoreSpaces) + SkipSpaces(scanner); return (scanner->cp < scanner->len); } diff --git a/fftma_module/gen/lib_src/scanadt.h b/fftma_module/gen/lib_src/scanadt.h old mode 100755 new mode 100644 index 1da09c3..eaa59b5 --- a/fftma_module/gen/lib_src/scanadt.h +++ b/fftma_module/gen/lib_src/scanadt.h @@ -60,7 +60,7 @@ * internal representation are hidden from the client. */ -typedef struct scannerCDT *scannerADT; +typedef struct scannerCDT* scannerADT; /* * Function: NewScanner @@ -151,7 +151,8 @@ void SaveToken(scannerADT scanner, string token); * of this option. */ -typedef enum { PreserveSpaces, IgnoreSpaces } spaceOptionT; +typedef enum { PreserveSpaces, + IgnoreSpaces } spaceOptionT; void SetScannerSpaceOption(scannerADT scanner, spaceOptionT option); spaceOptionT GetScannerSpaceOption(scannerADT scanner); diff --git a/fftma_module/gen/lib_src/simpio.c b/fftma_module/gen/lib_src/simpio.c old mode 100755 new mode 100644 index 07cab6a..e6769a0 --- a/fftma_module/gen/lib_src/simpio.c +++ b/fftma_module/gen/lib_src/simpio.c @@ -10,8 +10,8 @@ #include #include "genlib.h" -#include "strlib.h" #include "simpio.h" +#include "strlib.h" /* * Constants: @@ -43,14 +43,14 @@ int GetInteger(void) while (TRUE) { line = GetLine(); switch (sscanf(line, " %d %c", &value, &termch)) { - case 1: + case 1: printf("GetInteger llama a FreeBlock\n"); FreeBlock(line); return (value); - case 2: + case 2: printf("Unexpected character: '%c'\n", termch); break; - default: + default: printf("Please enter an integer\n"); break; } @@ -69,13 +69,13 @@ long GetLong(void) while (TRUE) { line = GetLine(); switch (sscanf(line, " %ld %c", &value, &termch)) { - case 1: + case 1: FreeBlock(line); return (value); - case 2: + case 2: printf("Unexpected character: '%c'\n", termch); break; - default: + default: printf("Please enter an integer\n"); break; } @@ -94,13 +94,13 @@ double GetReal(void) while (TRUE) { line = GetLine(); switch (sscanf(line, " %lf %c", &value, &termch)) { - case 1: + case 1: FreeBlock(line); return (value); - case 2: + case 2: printf("Unexpected character: '%c'\n", termch); break; - default: + default: printf("Please enter a real number\n"); break; } @@ -131,7 +131,7 @@ string GetLine(void) * twice the size of the previous one is allocated. */ -string ReadLine(FILE *infile) +string ReadLine(FILE* infile) { string line, nline; int n, ch, size; @@ -142,7 +142,7 @@ string ReadLine(FILE *infile) while ((ch = getc(infile)) != '\n' && ch != EOF) { if (n == size) { size *= 2; - nline = (string) GetBlock(size + 1); + nline = (string)GetBlock(size + 1); strncpy(nline, line, n); FreeBlock(line); line = nline; @@ -154,7 +154,7 @@ string ReadLine(FILE *infile) return (NULL); } line[n] = '\0'; - nline = (string) GetBlock(n + 1); + nline = (string)GetBlock(n + 1); strcpy(nline, line); FreeBlock(line); return (nline); diff --git a/fftma_module/gen/lib_src/simpio.h b/fftma_module/gen/lib_src/simpio.h old mode 100755 new mode 100644 index ca1efc9..133f2a4 --- a/fftma_module/gen/lib_src/simpio.h +++ b/fftma_module/gen/lib_src/simpio.h @@ -70,6 +70,6 @@ string GetLine(void); * is at the end-of-file position. */ -string ReadLine(FILE *infile); +string ReadLine(FILE* infile); #endif diff --git a/fftma_module/gen/lib_src/spherical.c b/fftma_module/gen/lib_src/spherical.c old mode 100755 new mode 100644 index a36e492..18a3712 --- a/fftma_module/gen/lib_src/spherical.c +++ b/fftma_module/gen/lib_src/spherical.c @@ -1,17 +1,16 @@ -#include -#include #include "genlib.h" - +#include +#include /*spherical covariance function*/ double spherical(double h) { - double z; + double z; - if (h >= 1.) { - z = 0.; - } else { - z = 1.-1.5*(double)h+0.5*(double)(h*h*h); - } - return (z); + if (h >= 1.) { + z = 0.; + } else { + z = 1. - 1.5 * (double)h + 0.5 * (double)(h * h * h); + } + return (z); } diff --git a/fftma_module/gen/lib_src/stable.c b/fftma_module/gen/lib_src/stable.c old mode 100755 new mode 100644 index 45c5a29..36608a2 --- a/fftma_module/gen/lib_src/stable.c +++ b/fftma_module/gen/lib_src/stable.c @@ -1,10 +1,9 @@ -#include -#include #include "genlib.h" - +#include +#include /*stable covariance function*/ double stable(double h, double alpha) { - return(exp(-3.*(double)pow(h,alpha))); + return (exp(-3. * (double)pow(h, alpha))); } diff --git a/fftma_module/gen/lib_src/stack.c b/fftma_module/gen/lib_src/stack.c old mode 100755 new mode 100644 index 8bd5ee4..1a5bc45 --- a/fftma_module/gen/lib_src/stack.c +++ b/fftma_module/gen/lib_src/stack.c @@ -7,9 +7,9 @@ * the type name stackElementT. */ -#include -#include "genlib.h" #include "stack.h" +#include "genlib.h" +#include /* * Constant: InitialStackSize @@ -33,7 +33,7 @@ */ struct stackCDT { - stackElementT *elements; + stackElementT* elements; int count; int size; }; @@ -63,13 +63,15 @@ void FreeStack(stackADT stack) void Push(stackADT stack, stackElementT element) { - if (stack->count == stack->size) ExpandStack(stack); + if (stack->count == stack->size) + ExpandStack(stack); stack->elements[stack->count++] = element; } stackElementT Pop(stackADT stack) { - if (StackIsEmpty(stack)) Error("Pop of an empty stack"); + if (StackIsEmpty(stack)) + Error("Pop of an empty stack"); return (stack->elements[--stack->count]); } @@ -108,7 +110,7 @@ stackElementT GetStackElement(stackADT stack, int index) static void ExpandStack(stackADT stack) { - stackElementT *array; + stackElementT* array; int i, newSize; newSize = stack->size * 2; diff --git a/fftma_module/gen/lib_src/stack.h b/fftma_module/gen/lib_src/stack.h old mode 100755 new mode 100644 index 2f765d1..e281990 --- a/fftma_module/gen/lib_src/stack.h +++ b/fftma_module/gen/lib_src/stack.h @@ -22,7 +22,7 @@ * be changed by editing this definition line. */ -typedef void *stackElementT; +typedef void* stackElementT; /* * Type: stackADT @@ -34,7 +34,7 @@ typedef void *stackElementT; * the underlying fields. */ -typedef struct stackCDT *stackADT; +typedef struct stackCDT* stackADT; /* * Function: NewStack diff --git a/fftma_module/gen/lib_src/strlib.c b/fftma_module/gen/lib_src/strlib.c old mode 100755 new mode 100644 index 6f00938..a669790 --- a/fftma_module/gen/lib_src/strlib.c +++ b/fftma_module/gen/lib_src/strlib.c @@ -17,9 +17,9 @@ * interface. */ +#include #include #include -#include #include "genlib.h" #include "strlib.h" @@ -59,7 +59,8 @@ char IthChar(string s, int i) { int len; - if (s == NULL) Error("NULL string passed to IthChar"); + if (s == NULL) + Error("NULL string passed to IthChar"); len = strlen(s); if (i < 0 || i > len) { Error("Index outside of string range in IthChar"); @@ -72,12 +73,16 @@ string SubString(string s, int p1, int p2) int len; string result; - if (s == NULL) Error("NULL string passed to SubString"); + if (s == NULL) + Error("NULL string passed to SubString"); len = strlen(s); - if (p1 < 0) p1 = 0; - if (p2 >= len) p2 = len - 1; + if (p1 < 0) + p1 = 0; + if (p2 >= len) + p2 = len - 1; len = p2 - p1 + 1; - if (len < 0) len = 0; + if (len < 0) + len = 0; result = CreateString(len); strncpy(result, s + p1, len); result[len] = '\0'; @@ -96,7 +101,8 @@ string CharToString(char ch) int StringLength(string s) { - if (s == NULL) Error("NULL string passed to StringLength"); + if (s == NULL) + Error("NULL string passed to StringLength"); return (strlen(s)); } @@ -104,7 +110,8 @@ string CopyString(string s) { string newstr; - if (s == NULL) Error("NULL string passed to CopyString"); + if (s == NULL) + Error("NULL string passed to CopyString"); newstr = CreateString(strlen(s)); strcpy(newstr, s); return (newstr); @@ -132,27 +139,36 @@ int StringCompare(string s1, string s2) int FindChar(char ch, string text, int start) { - char *cptr; - - if (text == NULL) Error("NULL string passed to FindChar"); - if (start < 0) start = 0; - if (start > strlen(text)) return (-1); + char* cptr; + + if (text == NULL) + Error("NULL string passed to FindChar"); + if (start < 0) + start = 0; + if (start > strlen(text)) + return (-1); cptr = strchr(text + start, ch); - if (cptr == NULL) return (-1); - return ((int) (cptr - text)); + if (cptr == NULL) + return (-1); + return ((int)(cptr - text)); } int FindString(string str, string text, int start) { - char *cptr; - - if (str == NULL) Error("NULL pattern string in FindString"); - if (text == NULL) Error("NULL text string in FindString"); - if (start < 0) start = 0; - if (start > strlen(text)) return (-1); + char* cptr; + + if (str == NULL) + Error("NULL pattern string in FindString"); + if (text == NULL) + Error("NULL text string in FindString"); + if (start < 0) + start = 0; + if (start > strlen(text)) + return (-1); cptr = strstr(text + start, str); - if (cptr == NULL) return (-1); - return ((int) (cptr - text)); + if (cptr == NULL) + return (-1); + return ((int)(cptr - text)); } /* Section 4 -- Case-conversion functions */ @@ -166,7 +182,8 @@ string ConvertToLowerCase(string s) Error("NULL string passed to ConvertToLowerCase"); } result = CreateString(strlen(s)); - for (i = 0; s[i] != '\0'; i++) result[i] = tolower(s[i]); + for (i = 0; s[i] != '\0'; i++) + result[i] = tolower(s[i]); result[i] = '\0'; return (result); } @@ -180,7 +197,8 @@ string ConvertToUpperCase(string s) Error("NULL string passed to ConvertToUpperCase"); } result = CreateString(strlen(s)); - for (i = 0; s[i] != '\0'; i++) result[i] = toupper(s[i]); + for (i = 0; s[i] != '\0'; i++) + result[i] = toupper(s[i]); result[i] = '\0'; return (result); } @@ -222,7 +240,8 @@ double StringToReal(string s) double result; char dummy; - if (s == NULL) Error("NULL string passed to StringToReal"); + if (s == NULL) + Error("NULL string passed to StringToReal"); if (sscanf(s, " %lg %c", &result, &dummy) != 1) { Error("StringToReal called on illegal number %s", s); } @@ -242,5 +261,5 @@ double StringToReal(string s) static string CreateString(int len) { - return ((string) GetBlock(len + 1)); + return ((string)GetBlock(len + 1)); } diff --git a/fftma_module/gen/lib_src/symtab.c b/fftma_module/gen/lib_src/symtab.c old mode 100755 new mode 100644 index 6e27a10..f762d5c --- a/fftma_module/gen/lib_src/symtab.c +++ b/fftma_module/gen/lib_src/symtab.c @@ -4,10 +4,10 @@ * This file implements the symbol table abstraction. */ -#include +#include "symtab.h" #include "genlib.h" #include "strlib.h" -#include "symtab.h" +#include /* * Constants @@ -25,8 +25,8 @@ typedef struct cellT { string key; - void *value; - struct cellT *link; + void* value; + struct cellT* link; } cellT; /* @@ -41,13 +41,13 @@ typedef struct cellT { */ struct symtabCDT { - cellT *buckets[NBuckets]; + cellT* buckets[NBuckets]; }; /* Private function declarations */ -static void FreeBucketChain(cellT *cp); -static cellT *FindCell(cellT *cp, string s); +static void FreeBucketChain(cellT* cp); +static cellT* FindCell(cellT* cp, string s); static int Hash(string s, int nBuckets); /* Public entries */ @@ -74,15 +74,15 @@ void FreeSymbolTable(symtabADT table) FreeBlock(table); } -void Enter(symtabADT table, string key, void *value) +void Enter(symtabADT table, string key, void* value) { int bucket; - cellT *cp; + cellT* cp; bucket = Hash(key, NBuckets); cp = FindCell(table->buckets[bucket], key); if (cp == NULL) { - cp = New(cellT *); + cp = New(cellT*); cp->key = CopyString(key); cp->link = table->buckets[bucket]; table->buckets[bucket] = cp; @@ -90,22 +90,23 @@ void Enter(symtabADT table, string key, void *value) cp->value = value; } -void *Lookup(symtabADT table, string key) +void* Lookup(symtabADT table, string key) { int bucket; - cellT *cp; + cellT* cp; bucket = Hash(key, NBuckets); cp = FindCell(table->buckets[bucket], key); - if (cp == NULL) return(UNDEFINED); + if (cp == NULL) + return (UNDEFINED); return (cp->value); } void MapSymbolTable(symtabFnT fn, symtabADT table, - void *clientData) + void* clientData) { int i; - cellT *cp; + cellT* cp; for (i = 0; i < NBuckets; i++) { for (cp = table->buckets[i]; cp != NULL; cp = cp->link) { @@ -125,9 +126,9 @@ void MapSymbolTable(symtabFnT fn, symtabADT table, * this function must free the string storage as well. */ -static void FreeBucketChain(cellT *cp) +static void FreeBucketChain(cellT* cp) { - cellT *next; + cellT* next; while (cp != NULL) { next = cp->link; @@ -146,7 +147,7 @@ static void FreeBucketChain(cellT *cp) * returned. If no match is found, the function returns NULL. */ -static cellT *FindCell(cellT *cp, string key) +static cellT* FindCell(cellT* cp, string key) { while (cp != NULL && !StringEqual(cp->key, key)) { cp = cp->link; diff --git a/fftma_module/gen/lib_src/symtab.h b/fftma_module/gen/lib_src/symtab.h old mode 100755 new mode 100644 index c83c307..c03b93b --- a/fftma_module/gen/lib_src/symtab.h +++ b/fftma_module/gen/lib_src/symtab.h @@ -15,7 +15,7 @@ * This type is the ADT used to represent a symbol table. */ -typedef struct symtabCDT *symtabADT; +typedef struct symtabCDT* symtabADT; /* * Type: symtabFnT @@ -24,8 +24,8 @@ typedef struct symtabCDT *symtabADT; * map over the entries in a symbol table. */ -typedef void (*symtabFnT)(string key, void *value, - void *clientData); +typedef void (*symtabFnT)(string key, void* value, + void* clientData); /* Exported entries */ @@ -55,7 +55,7 @@ void FreeSymbolTable(symtabADT table); * Each call to Enter supersedes any previous definition for key. */ -void Enter(symtabADT table, string key, void *value); +void Enter(symtabADT table, string key, void* value); /* * Function: Lookup @@ -65,7 +65,7 @@ void Enter(symtabADT table, string key, void *value); * table, or UNDEFINED, if no such value exists. */ -void *Lookup(symtabADT table, string key); +void* Lookup(symtabADT table, string key); /* * Function: MapSymbolTable @@ -80,6 +80,6 @@ void *Lookup(symtabADT table, string key); */ void MapSymbolTable(symtabFnT fn, symtabADT table, - void *clientData); + void* clientData); #endif diff --git a/fftma_module/gen/lib_src/test_fact.c b/fftma_module/gen/lib_src/test_fact.c old mode 100755 new mode 100644 index d2456ee..21af512 --- a/fftma_module/gen/lib_src/test_fact.c +++ b/fftma_module/gen/lib_src/test_fact.c @@ -1,26 +1,24 @@ #include "genlib.h" /*tries factor*/ -int test_fact(int *pnum, int fact, int *pmaxfac) +int test_fact(int* pnum, int fact, int* pmaxfac) { - int power, t; - - power = 0; - while ( ( t = *pnum / fact ) * fact == *pnum ) { - ++power; - *pnum = t; + int power, t; + + power = 0; + while ((t = *pnum / fact) * fact == *pnum) { + ++power; + *pnum = t; } - if ( power != 0 ) { - if (fact > *pmaxfac) - *pmaxfac = fact; - } + if (power != 0) { + if (fact > *pmaxfac) + *pmaxfac = fact; + } - if ( t > fact ) { - return (1); - } + if (t > fact) { + return (1); + } - return (0); + return (0); } - - diff --git a/fftma_module/gen/lib_src/testmemory.c b/fftma_module/gen/lib_src/testmemory.c old mode 100755 new mode 100644 index 7cbd10e..3d7207f --- a/fftma_module/gen/lib_src/testmemory.c +++ b/fftma_module/gen/lib_src/testmemory.c @@ -1,15 +1,15 @@ -#include -#include -#include +#include #include +#include +#include #include -#include +#include -void testmemory(double *realint) +void testmemory(double* realint) { - if (realint == NULL) { - printf("Testmemory.c: No memory available \n"); - exit; - } - return; + if (realint == NULL) { + printf("Testmemory.c: No memory available \n"); + exit; + } + return; }