Compare commits

...

91 Commits

Author SHA1 Message Date
chortas 7b1dfecc1d Merge branch 'main' of github.com:medios-porosos-fiuba/simulacion-permeabilidad into improvement-double_with_dirty_flag 3 years ago
chortas ca9c01c962 Change chunk size 3 years ago
chortas 3246963979 Add tupac sh 3 years ago
Oli f5a43fe455 add mem profile from outside 3 years ago
Oli d464ac5cb9 change dirty flag to false on update 3 years ago
Oli 2fd7368bca add dirty flag 3 years ago
chortas 64a2c0f3c4 Delete core variable in generation file 3 years ago
chortas 18a2d52f2b Change tolerance test 3 years ago
Oli c0d4ec18ab change sample size 3 years ago
Oli 1b8c56489c update makefile 3 years ago
Oli 8399e1971d add perf test 3 years ago
Oli 213612e597 cleaning 3 years ago
Oli c3f2939bf8 chunk array clean 3 years ago
Oli 92fcc3bb5b fix 3 years ago
Oli bcbfaf9f75 chunk not working 3 years ago
Oli 90c8673110 double chunk not working 3 years ago
chortas 7d68fc075e Fix test 3 years ago
Cecilia Hortas 35570fb088
Merge pull request #7 from medios-porosos-fiuba/removing-logs 3 years ago
Cecilia Hortas f49258f4aa
Merge pull request #6 from medios-porosos-fiuba/debugging-fourt
Debugging fourt
3 years ago
Cecilia Hortas 36e64fa03e
Merge pull request #5 from medios-porosos-fiuba/removing-vector-2
Removing vector 2
3 years ago
chortas 83affab77b Working for 16 3 years ago
chortas 6b4bb901e2 Revert "wip"
This reverts commit 897690a1e4.
3 years ago
chortas 897690a1e4 wip 3 years ago
Oli e626ff6c93 covar test 3 years ago
chortas 44a16ce9b9 Add more prints.. 3 years ago
chortas 90f76bb61d Add more prints.. 3 years ago
chortas 0a7f2fe1ba Add prints for saving as well 3 years ago
chortas 9d92fa5be1 Add corrida 3 years ago
chortas 30a14002ae Add indexes to print 3 years ago
chortas 4a90597609 Add numeration to prints and remove extra prints 3 years ago
chortas 7ef3c2d2e5 Debugging fourt 3 years ago
chortas c64a61e071 Add logs to debug and fix seg fault 3 years ago
chortas dcbaafdb1a At least compiling. Segfault 3 years ago
chortas 0d40eae129 Commit to debug 3 years ago
chortas 8ab0002e44 Remove logs 3 years ago
chortas bfdc81ee4b Fix 3 years ago
chortas e54ac71486 First approach 3 years ago
chortas f35afe06f2 Add improvements 3 years ago
Oli f4ccbbff6f add more conclusions 3 years ago
Oli 466b480109 add return values analysis 3 years ago
chortas b5ae9abdc6 Add cpu analysis 3 years ago
chortas a2f794772e Fix santisi bug and reorganize analysis 3 years ago
Oli 9c5f09e2f2 add memory analisis and overall comparisons 3 years ago
chortas 077bef4395 Fix cpu bug 3 years ago
Oli e0f7d0000e add 128 plots 3 years ago
Oli 245e171e46 add params distribution 3 years ago
Oli 9fbf57b8c8 fix 128 3 years ago
Oli 3deb635c41 add treemap and make cleaner piecharts 3 years ago
chortas ca2a4f133d Fill analysis and remove extra newline 3 years ago
chortas 3be55fe2a1 Add CPU to all files 3 years ago
chortas ec130dddd5 Fix bug in gasdev 3 years ago
chortas f43eec7034 Fix seg fault 3 years ago
chortas 96bcbf2c95 Add CPU analysis 3 years ago
chortas 036d958eb8 Add merge of 256 3 years ago
chortas 757a8befe2 Delete extra data 3 years ago
chortas 5c5930376a Fix merge conflict 3 years ago
chortas e88a091b10 Add 256 case 3 years ago
Oli d7c2fb69c7 add gasdev pie chart 3 years ago
Oli f4efc59d85 add pie charts to compare proportions 3 years ago
chortas 0f8cbc686b Save logs with information and reorganize notebook 3 years ago
Oli 4e338ab748 plot memory 3 years ago
chortas a3d7a23da1 Add memory diff instead of total 3 years ago
chortas 6510ea7a88 Add memory profiler to a file 3 years ago
Oli e5c70006f4 Merge branch 'milestone_5' of github.com:medios-porosos-fiuba/simulacion-permeabilidad into milestone_5 3 years ago
Oli eaeb01d343 plot times 3 years ago
chortas 0c109fb4db Add memory profiler to all files 3 years ago
chortas fb20a57560 Add analysis 3 years ago
chortas 6567ec3afc Add elapsed_time to prebuild_gwn 3 years ago
chortas 5914f9dd07 Make logs optional 3 years ago
chortas 14c8458474 Add configurable number of cores 3 years ago
chortas 8970e8927c Add memory profiler 3 years ago
chortas 2ec7f40125 Add elapsed time in logs 3 years ago
chortas 1e9aefebf0 Add logs 3 years ago
Oli e52630f86d garbage collect array after generating 3 years ago
Oli 96e1dcd704 change function to create array 3 years ago
Oli fd71c3faab add small script to test fftma 3 years ago
chortas 9ebb8b1e24 Free some memory 3 years ago
chortas ab8d288b09 Add compiler flags 3 years ago
chortas 4fd1d5ddc6 Linter to c files 3 years ago
chortas 83d7d29273 Linter to c files 3 years ago
Oli 3210b268aa rm refine module 3 years ago
Oli f152649da5 rm files 3 years ago
Oli 6b97a9dfeb rm files 3 years ago
Oli 562b948508 rm files 3 years ago
Oli 6ac4de81a9 rm files 3 years ago
Oli 2593883e22 rm files 3 years ago
Oli e01ef454b0 rm files 3 years ago
Oli 15a19c5935 fix 3 years ago
Oli 028d0a8d9f intermediate commit 3 years ago
Oli de5f4c8fe6 delete extra files 3 years ago
chortas 6d62c0e798 Run linter 3 years ago

5
.gitignore vendored

@ -17,3 +17,8 @@ fftma_module/gen/build/
tools/connec/conec2d
tools/connec/conec3d
utilities/__pycache__/
fftma_module/gen/log_*
fftma_module/gen/out*.npy
.ipynb_checkpoints/analysis-checkpoint.ipynb
fftma_module/gen/.ipynb_checkpoints/
__pycache__/

@ -2,7 +2,7 @@
install:
. ./script_install.sh
run:
run: binaries
./run_simulation.sh
fftma:
@ -11,9 +11,15 @@ fftma:
binaries:
./script_fortran.sh
test: binaries
test: binaries fftma
cd tests/integration && python3 -m unittest test.py
perf: binaries
cd tests/performance && python3 generation.py
cd tests/performance && python3 connectivity.py
cd tests/performance && python3 connectivity.py
test-gen: fftma
cd tests/stages/generation && python3 -m unittest test.py
time-gen: fftma
cd tests/performance/generation && python3 time.py $(N)

@ -18,18 +18,32 @@ make binaries
make run
```
## Correr los casos de prueba
## Pruebas de integración
```
make test
```
## Pruebas unitarias
### Generacion
```
make test-gen
```
## Correr las pruebas de performance
```
make perf
```
### Obtener tiempos de generacion
```
make time-gen N=<tamaño>
```
## Instalar el binario FFTMA
```

@ -2,16 +2,16 @@
simDir=output
startJob=0
[Iterables]
p=[10,39,15]
seeds=[5462,2]
p=[10,39,1]
seeds=[5462,1]
lc=[4]
connectivity=[1,2,3]
connectivity=[1]
variances=[1]
[Generation]
Nx = 16
Ny = 16
Nz = 16
Nx = 8
Ny = 8
Nz = 8
variogram_type=1
binary = yes
kh = 100
@ -22,23 +22,23 @@ genera=yes
[Connectivity]
keep_aspect= yes
keep_aspect= no
block_size = 4
indicators_MinBlockSize =4
Max_sample_size = 12
compGconec= 1
conec=yes
conec=no
[Solver]
num_of_cores = 1
ref = 2
solve = yes
solve = no
rtol = 1e-4
[K-Postprocess]
MinBlockSize =1
Max_sample_size = 4
kperm=yes
postprocess=yes
SaveVfield=yes
kperm=no
postprocess=no
SaveVfield=no

@ -1,36 +0,0 @@
from FFTMA import gen
from refine import refine as ref
import numpy as np
from random import randint as rdi
N=128
for i in range(1):
nx, ny, nz = N,N,N
dx, dy, dz = 1.0, 1.0, 1.0
seed= 1548762 #rdi(10000,99999)
var=1
vario=2
alpha=1
lcx=2
lcy=4
lcz=16
ap1x=1
ap1y=0
ap1z=0
ap2x=0
ap2y=1
ap2z=0
v1 = (var, vario, alpha, lcx, lcy, lcz, ap1x, ap1y, ap1z, ap2x, ap2y, ap2z)
variograms = [v1]
mean=15.3245987
variance=3.5682389
typ=1
k=gen(nx, ny, nz, dx, dy, dz, seed, variograms, mean, variance, typ)
np.save("out"+str(i)+".npy",ref(k,4,4,4))

@ -1,442 +0,0 @@
#include "geostat.h"
#ifndef _CONDOR_H_
#define _CONDOR_H_
/*=======================================================*/
/*FUNCTIONS*/
/*---------*/
/*addstat */
/*adds mean and variance effects to the N(0,1) realization*/
/*input: */
/*realin: structure defining a realization - */
/* must have zero mean and unit variance */
/*stat: structure defining the mean and variance */
/*output: */
/*realout: structure defining a realization - */
void addstat(struct realization_mod *realin,struct statistic_mod stat ,struct realization_mod *realout);
/*condit */
/*conditioning of a Gaussian realization using */
/*the primal kriging approach */
/*inversion method = CONJUGATE GRADIENTS */
/*NB: The realization and the well data are */
/*reduced to the standard Gaussian N(0,1) */
/*INPUT */
/*k: measurement type */
/*well: structure with the well data */
/*variogram: structure describing the variogram model*/
/*realin: structure defining the input realization */
/*grid: structure defining the grid */
/*OUTPUT */
/*realout: structure defining the output realization */
void condit(int k,struct welldata_mod well,struct vario_mod variogram,struct grid_mod grid,struct realization_mod *realin,struct realization_mod *realout);
/*d_kriging */
/*interpolation from ponctual data with simple dual */
/*kriging. d_kriging is used to interpolate standard */
/*Gaussian data only */
/*inversion method = CONJUGATE GRADIENTS */
/*INPUT */
/*k: measurement type */
/*well: structure with the well data */
/*variogram: structure describing the variogram model*/
/*grid: structure defining the grid */
/*OUTPUT */
/*realout: structure defining the output realization */
void d_kriging(int k,struct welldata_mod well,struct vario_mod variogram,struct grid_mod grid,struct realization_mod *realout);
/*dual_condit */
/*conditioning of a Gaussian realization using */
/*the dual kriging approach */
/*inversion method = CONJUGATE GRADIENTS */
/*NB: The realization and the well data are */
/*reduced to the standard Gaussian N(0,1) */
/*INPUT */
/*k: measurement type */
/*well: structure with the well data */
/*variogram: structure describing the variogram model*/
/*realin: structure defining the input realization */
/*OUTPUT */
/*realout: structure defining the output realization */
void dual_condit(int k,struct welldata_mod well,struct vario_mod variogram,struct grid_mod grid,struct realization_mod *realin,struct realization_mod *realout);
/*dual_condit_gen */
/*conditioning of a Gaussian realization using */
/*the dual kriging approach */
/*considered cases : Gaussian data */
/* Uniform data */
/* Log */
/*inversion method = CONJUGATE GRADIENTS */
/*INPUT */
/*k: measurement type */
/*wellin: structure with the well data */
/*variogram: structure describing the variogram model*/
/*realin: structure defining the input realization */
/*stat : structure defining the statistical data */
/*OUTPUT */
/*realout: structure defining the output realization */
void dual_condit_gen(int k,struct welldata_mod wellin,struct vario_mod variogram,struct grid_mod grid,struct realization_mod *realin,struct statistic_mod stat, struct realization_mod *realout);
/*FAST FOURIER TRANSFORM MOVING AVERAGE METHOD */
/*Turns a Gaussian white noise vector into a */
/*spatially correlated vector */
/*input: */
/*variogram: structure defining the variogram */
/* model */
/*grid: structure defining the grid */
/*n: vector with the number of cells along the */
/* X, Y and Z axes for the underlying grid */
/* i = [0 1 2] */
/* --> 0 0 0 : n will be computed and */
/* updated as output */
/* --> nx ny nz: these dimensions are used */
/*realin: structure defining a realization - */
/* must be a Gaussian white noise */
/*output: */
/*realout: structure defining a realization - */
void FFTMA(struct vario_mod variogram,struct grid_mod grid,int n[3],struct realization_mod *realin,struct realization_mod *realout);
/*FAST FOURIER TRANSFORM MOVING AVERAGE METHOD */
/*PLUS GRADIENT CALCULATIONS - GRADIENTS ARE */
/*CALCULATED FOR THE GRADUAL DEFORMATION PARAMETERS */
/* */
/*INPUT: */
/*------- */
/*variogram: structure defining the variogram */
/* model */
/*grid: structure defining the grid */
/*n: vector with the number of cells along the */
/* X, Y and Z axes for the underlying grid */
/* i = [0 1 2] */
/* --> 0 0 0 : n will be computed and */
/* updated as output */
/* --> nx ny nz: these dimensions are used */
/*realin: structure defining a realization - */
/* must be a Gaussian white noise */
/*DZoverDRHO: gradients calculated durin the gradual */
/* deformation process */
/* */
/*OUTPUT: */
/*------- */
/*realout: structure defining a realization - */
/*DYoverDRHO : structure with the output gradients */
void FFTMA_gradient(struct vario_mod variogram,struct grid_mod grid,int n[3],struct realization_mod *realin,struct gradients_mod *DZoverDRHO,struct realization_mod *realout,struct gradients_mod *DYoverDRHO);
/* GENERATION OF A GAUSSIAN WHITE NOISE VECTOR */
/*input: */
/* seed: seed */
/* n: number of components in the vector */
/*output: */
/* realization: structure defining the realization*/
void generate(long *seed, int n, struct realization_mod *realization);
void ikrig(int option, struct statfacies_mod facies, int kk,struct welldata_mod data, struct variotable_mod variogram,struct grid_mod grid,float *krigout);
/*kriging */
/*interpolation from ponctual data with simple */
/*kriging. kriging is used to interpolate standard */
/*Gaussian data only */
/*inversion method = CONJUGATE GRADIENTS */
/*INPUT */
/*k: measurement type */
/*well: structure with the well data */
/*variogram: structure describing the variogram model*/
/*grid: structure defining the grid */
/*OUTPUT */
/*realout: structure defining the output realization */
void kriging(int k,struct welldata_mod data,struct vario_mod variogram,struct grid_mod grid,struct realization_mod *realout);
/*CHOLESKY SIMULATION TECHNIQUE - LUSIM */
/*appropriate for regular grids and small number */
/*of points */
/*input: */
/*variogram: structure defining the variogram model */
/*realin: structure defining a realization - */
/* must be a Gaussian white noise */
/*output: */
/*realout: structure defining a realization - */
void LUSIM(struct vario_mod variogram,struct grid_mod grid,struct realization_mod *realin,struct realization_mod *realout);
void nor2any(struct cdf_mod *pcumulf,struct realization_mod *realin,struct realization_mod *realout);
/*TURNS NORMAL NUMBERS INTO LOGNORMAL NUMBERS */
/*input: */
/*realin: structure defining a realization - */
/* normal numbers */
/*typelog: --> 3: lognormal (natural) */
/* --> 4: lognormal (log10) */
/*output: */
/*realout: structure defining a realization - */
/* lognormal numbers */
void nor2log(struct realization_mod *realin, int typelog, struct realization_mod *realout);
/*TURNS NORMAL NUMBERS INTO UNIFORM NUMBERS */
/*input: */
/*realin: structure defining a realization - */
/* must have zero mean and unit variance */
/*output: */
/*realout: structure defining a realization - */
/* uniform numbers */
void nor2unif(struct realization_mod *realin,struct realization_mod *realout);
/*o_kriging */
/*interpolation from ponctual data with ordinary kriging */
/*o_kriging is used to interpolate standard Gaussian data*/
/*only */
/*INPUT */
/*k: measurement type or rank */
/*variogram: structure describing the variogram model */
/*data: structure with the well data */
/*grid: structure defining the grid */
/*OUTPUT */
/*realout: structure defining the kriged realization */
void o_kriging(int k,struct vario_mod variogram,struct welldata_mod data,struct grid_mod grid,struct realization_mod *realout);
/*od_kriging */
/*interpolation from ponctual data with ordinary dual*/
/*kriging. od_kriging is used to interpolate */
/*Gaussian data only */
/*inversion method = CONJUGATE GRADIENTS */
/*INPUT */
/*k: measurement type */
/*well: structure with the well data */
/*variogram: structure describing the variogram model*/
/*grid: structure defining the grid */
/*OUTPUT */
/*realout: structure defining the output realization */
void od_kriging(int k,struct welldata_mod well,struct vario_mod variogram,struct grid_mod grid,struct realization_mod *realout);
/*PLURIGAUSSIAN METHOD */
/*Starting from two independent Gaussian realizations Y1*/
/*and Y2 with zero mean and unit variance, but possibly */
/*different variogram models, a categorial realization is */
/*built. */
/*Truncation is performed on the basis of a mask with the */
/*x-axis associated to Y1 and the y-axis associated to Y2 */
/*This truncation may be stationary or non stationary */
/*The mask consists of lines perpendicular to the Y1 and */
/*Y2 axes. They define rectangular subregions attributed */
/*to categories i = [1...ncat] */
/*INPUT: */
/*facies: structure defining the facies proportions in the*/
/* whole reservoir model */
/*ineq: structure defining the mask for truncation */
/*Y1: first starting Gaussian realization */
/*Y2: second starting realization - same length as Y1 */
/*OUTPUT: */
/*Y: output realization */
void plurigau(struct statfacies_mod facies,struct inequalities_mod ineq,struct realization_mod *Y1, struct realization_mod *Y2, struct realization_mod *Y);
/*tr_kriging */
/*kriging with an external drift */
/*The trend is limited to a 2-term function */
/*m(u) = a0+a1y(u) where y(u) is a secondary */
/*external variable (e.g., seismic) */
/*INPUT: */
/*k: measurement type or rank */
/*variogram: structure describing the variogram */
/* model */
/*data: structure with the well data */
/*grid: structure defining the grid */
/*y: secondary external variable */
/*OUTPUT: */
/*realout: kriged realization */
void tr_kriging(int k, struct vario_mod variogram,struct welldata_mod data,struct grid_mod grid,struct realization_mod *Y,struct realization_mod *realout);
/*truncat */
/*TRUNCATES A STANDARD GAUSSIAN REALIZATION */
/*WITH THE TRUNCATED GAUSSIAN METHOD */
/*truncation may be stationary or nonstationary*/
/*INPUT: */
/*facies: description of the facies proportions*/
/*thresholds: may be an input if they have been*/
/* previously computed */
/*realin: structure defining a realization - */
/* must be a standard normal */
/*OUTPUT: */
/*realout: structure defining a realization - */
/*thresholds: Gaussian thresholds corresponding*/
/* to the facies volume fractions if*/
void truncat(struct statfacies_mod facies,struct realization_mod *realin,struct realization_mod *realout, struct statfacies_mod *thresholds);
/*TURNS UNIFORM NUMBERS INTO NORMAL NUMBERS*/
/*INPUT: */
/*realin: structure defining a realization - */
/* uniform numbers */
/*OUTPUT: */
/*realout: structure defining a realization - */
/* must have zero mean and unit variance */
void unif2nor(struct realization_mod *realin,struct realization_mod *realout);
/*Wany2nor */
/*converts any kind of continuous data to Gaussian data*/
/*INPUT: */
/*k: type of measurement */
/*wellin: structure with the input well data */
/*seed: seed */
/*OUTPUT: */
/*wellout: structure with the output well data */
/*pcumulf: structure describing the experimental */
/* cumulative distribution */
void Wany2nor(int k, struct welldata_mod wellin, long *seed,struct welldata_mod *wellout, struct cdf_mod *pcumulf);
/*Wfac2nor0 */
/*turns the facies well measurements into Gaussian data */
/*intervals only are respected - spatial variability is */
/*disregarded */
/*USES vf2thres, deflimit,trungasdev */
/*This step is required before performing kriging */
/*input: */
/*k: type measurement, may be 1,2,3 ... or ntype */
/*wellin: structure defining well data */
/* code = 5 --> facies */
/*grid: structure defining the grid */
/*facies: structure defining the facies proportions */
/*seed: seed */
/*output: */
/*wellout: structure defining well data - the kth vector*/
/*of measures are converted to standard normal data */
void Wfac2nor0(int k,struct grid_mod grid,struct welldata_mod wellin,struct statfacies_mod facies,long *seed,struct welldata_mod *wellout);
/*Wfac2nor1 */
/*turns the facies well measurements into Gaussian data */
/*Proceeds in 2 steps */
/* 1- builds a realization Y = LZ (variability ensured) */
/* 2- proposes successively local changes for Z to */
/* respect the intervals-Metropoils-Hastings approach */
/*input: */
/*k: type measurement, may be 1,2,3 ... or ntype */
/*variogram: structure defining the variogram model */
/*wellin: structure defining well data */
/* code = 5 --> facies */
/*grid: structure defining the grid */
/*facies: structure defining the facies proportions */
/*seed: seed */
/*output: */
/*wellout: structure defining well data - the kth vector*/
/*of measures are converted to standard normal data */
void Wfac2nor1(int k,struct vario_mod variogram,struct grid_mod grid,struct welldata_mod wellin,struct statfacies_mod facies,long *seed,struct welldata_mod *wellout);
/*Wfac2norPG0 */
/*turns the facies well measurements into Gaussian data */
/*intervals only are respected - spatial variability is */
/*disregarded */
/*method developed for pluriGaussian realizations */
/*stationary and nonstationary truncation */
/*This step is required before performing kriging */
/*INPUT: */
/*k: type measurement, may be 1,2,3 ... or ntype */
/*wellin: structure defining well data */
/* code = 5 --> facies */
/*grid: structure defining the grid */
/*facies: structure defining the facies proportions */
/*seed: seed */
/*idreal: identification of the starting realization */
/* = 1 for realization 1 */
/* = 2 for realization 2 */
/*ineq: structure defining the mask for truncation */
/*OUTPUT: */
/*wellout: structure defining well data - the kth vector*/
/*of measures are converted to standard normal data */
void Wfac2norPG0(int k,struct inequalities_mod ineq,struct grid_mod grid,struct welldata_mod wellin,struct statfacies_mod facies,long *seed,int idreal,struct welldata_mod *wellout);
/*Wfac2norPG1 */
/*turns the facies well measurements into Gaussian data */
/*Proceeds in 2 steps */
/* 1- builds a realization Y = LZ (variability ensured) */
/* 2- proposes successively local changes for Z to */
/* respect the intervals-Metropoils-Hastings approach */
/*method developed for pluriGaussian realizations */
/*STATIONARY TRUNCATION */
/*INPUT: */
/*k: type measurement, may be 1,2,3 ... or ntype */
/*wellin: structure defining well data */
/* code = 5 --> facies */
/*grid: structure defining the grid */
/*facies: structure defining the facies proportions */
/*seed: seed */
/*ineq: structure describing the truncation mask */
/*variogram1: structure describing the variogram model */
/* for realization Y1 */
/*variogram2: structure describing the variogram model */
/* for realization Y2 */
/*ineq: structure defining the mask for truncation */
/*OUTPUT: */
/*wellout1: structure with converted well data for */
/*realization 1 - the kth vector */
/*wellout2: structure with converted well data for */
/*realization 2 - the kth vector */
void Wfac2norPG1(int k,struct inequalities_mod ineq,struct vario_mod variogram1,struct vario_mod variogram2,struct grid_mod grid,struct welldata_mod wellin,struct statfacies_mod facies,long *seed,struct welldata_mod *wellout1,struct welldata_mod *wellout2);
/*Wlog2nor */
/*turns the log well measurements into std Gaussian data*/
/*This step is required before performing kriging */
/*input: */
/*k: type measurement, may be 1,2,3 ... or ntype */
/*wellin: structure defining well data */
/* code = 3 --> natural log */
/* code = 4 --> logarithm 10 */
/*grid: structure defining the grid */
/*stat: structure defining the statistics */
/*output: */
/*wellout: structure defining well data - the kth vector*/
/*of measures are converted to standard normal data */
void Wlog2nor(int k, struct welldata_mod wellin,struct grid_mod grid,struct statistic_mod stat,struct welldata_mod *wellout);
/*Wnor2nor */
/*converts Well Normal numbers to the standard Gaussian */
/*input: */
/*k: type measurement, may be 1,2,3 ... or ntype */
/*wellin: structure defining well data */
/*grid: structure defining the grid */
/*stat: structure defining the statistics */
/*output: */
/*wellout: structure defining well data - the kth vector*/
/*of measures are converted to standard normal data */
void Wnor2nor(int k, struct welldata_mod wellin,struct grid_mod grid,struct statistic_mod stat,struct welldata_mod *wellout);
/*Wunif2nor */
/*converts Well uniform numbers to the standard Gaussian*/
/*input: */
/*k: type measurement, may be 1,2,3 ... or ntype */
/*wellin: structure defining well data */
/*output: */
/*wellout: structure defining well data - the kth vector*/
/*of measures are converted to standard normal data */
void Wunif2nor(int k, struct welldata_mod wellin,struct welldata_mod *wellout);
#endif // define _CONDOR_H_

@ -1,163 +0,0 @@
/*
* File: genlib.h
* Version: 1.0
* Last modified on Fri Jul 15 15:45:52 1994 by eroberts
* -----------------------------------------------------
* This file contains several definitions that form the
* core of a general-purpose ANSI C library developed by Eric
* Roberts. The goal of this library is to provide a basic
* set of tools and conventions that increase the readability
* of C programs, particularly as they are used in a teaching
* environment.
*
* The basic definitions provided by genlib.h are:
*
* 1. Declarations for several new "primitive" types
* (most importantly bool and string) that are
* used throughout the other libraries and
* applications as fundamental types.
*
* 2. A new set of functions for memory allocation.
*
* 3. A function for error handling.
*/
#ifndef _genlib_h
#define _genlib_h
#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
/* Section 1 -- Define new "primitive" types */
/*
* Type: bool
* ----------
* This type has two values, FALSE and TRUE, which are equal to 0
* and 1, respectively. Most of the advantage of defining this type
* comes from readability because it allows the programmer to
* provide documentation that a variable will take on only one of
* these two values. Designing a portable representation, however,
* is surprisingly hard, because many libraries and some compilers
* define these names. The definitions are usually compatible but
* may still be flagged as errors.
*/
#ifdef THINK_C
typedef int bool;
#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
#endif
/*
* Type: string
* ------------
* The type string is identical to the type char *, which is
* traditionally used in C programs. The main point of defining a
* new type is to improve program readability. At the abstraction
* levels at which the type string is used, it is usually not
* important to take the string apart into its component characters.
* Declaring it as a string emphasizes this atomicity.
*/
typedef char *string;
/*
* Constant: UNDEFINED
* -------------------
* Besides NULL, the only other constant of pointer type is
* UNDEFINED, which is used in certain packages as a special
* sentinel to indicate an undefined pointer value. In many
* such contexts, NULL is a legitimate data value and is
* therefore inappropriate as a sentinel.
*/
#define UNDEFINED ((void *) undefined_object)
extern char undefined_object[];
/* Section 2 -- Memory allocation */
/*
* General notes:
* --------------
* These functions provide a common interface for memory
* allocation. All functions in the library that allocate
* memory do so using GetBlock and FreeBlock. Even though
* the ANSI standard defines malloc and free for the same
* purpose, using GetBlock and FreeBlock provides greater
* compatibility with non-ANSI implementations, automatic
* out-of-memory error detection, and the possibility of
* substituting a garbage-collecting allocator.
*/
/*
* Function: GetBlock
* Usage: ptr = (type) GetBlock(nbytes);
* -------------------------------------
* GetBlock allocates a block of memory of the given size. If
* no memory is available, GetBlock generates an error.
*/
void *GetBlock(size_t nbytes);
/*
* Function: FreeBlock
* Usage: FreeBlock(ptr);
* ----------------------
* FreeBlock frees the memory associated with ptr, which must
* have been allocated using GetBlock, New, or NewArray.
*/
void FreeBlock(void *ptr);
/*
* Macro: New
* Usage: p = New(pointer-type);
* -----------------------------
* The New pseudofunction allocates enough space to hold an
* object of the type to which pointer-type points and returns
* a pointer to the newly allocated pointer. Note that
* "New" is different from the "new" operator used in C++;
* the former takes a pointer type and the latter takes the
* target type.
*/
#define New(type) ((type) GetBlock(sizeof *((type) NULL)))
/*
* Macro: NewArray
* Usage: p = NewArray(n, element-type);
* -------------------------------------
* NewArray allocates enough space to hold an array of n
* values of the specified element type.
*/
#define NewArray(n, type) ((type *) GetBlock((n) * sizeof(type)))
/* Section 3 -- Basic error handling */
/*
* Function: Error
* Usage: Error(msg, ...)
* ----------------------
* Error generates an error string, expanding % constructions
* appearing in the error message string just as printf does.
* After printing the error message, the program terminates.
*/
void Error(string msg, ...);
#endif

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

@ -1,21 +0,0 @@
#include <stdlib.h>
#include <string.h>
#include "genlib.h"
#include "geostat.h"
#ifndef _PRESSURE_H
#define _PRESSURE_H
/*STRUCTURES*/
/*----------*/
/* pressure_mod */
/* x: macroscopic gradient value in x direction */
/* y: macroscopic gradient value in x direction */
/* z: macroscopic gradient value in x direction */
struct pressure_mod {
double x,y,z;
};
#endif // define _PRESSURE_H

@ -1,75 +0,0 @@
/*
* File: simpio.h
* Version: 1.0
* Last modified on Wed Apr 27 07:29:13 1994 by eroberts
* -----------------------------------------------------
* This interface provides access to a simple package of
* functions that simplify the reading of input data.
*/
#ifndef _simpio_h
#define _simpio_h
#include "genlib.h"
/*
* Function: GetInteger
* Usage: i = GetInteger();
* ------------------------
* GetInteger reads a line of text from standard input and scans
* it as an integer. The integer value is returned. If an
* integer cannot be scanned or if more characters follow the
* number, the user is given a chance to retry.
*/
int GetInteger(void);
/*
* Function: GetLong
* Usage: l = GetLong();
* ---------------------
* GetLong reads a line of text from standard input and scans
* it as a long integer. The value is returned as a long.
* If an integer cannot be scanned or if more characters follow
* the number, the user is given a chance to retry.
*/
long GetLong(void);
/*
* Function: GetReal
* Usage: x = GetReal();
* ---------------------
* GetReal reads a line of text from standard input and scans
* it as a double. If the number cannot be scanned or if extra
* characters follow after the number ends, the user is given
* a chance to reenter the value.
*/
double GetReal(void);
/*
* Function: GetLine
* Usage: s = GetLine();
* ---------------------
* GetLine reads a line of text from standard input and returns
* the line as a string. The newline character that terminates
* the input is not stored as part of the string.
*/
string GetLine(void);
/*
* Function: ReadLine
* Usage: s = ReadLine(infile);
* ----------------------------
* ReadLine reads a line of text from the input file and
* returns the line as a string. The newline character
* that terminates the input is not stored as part of the
* string. The ReadLine function returns NULL if infile
* is at the end-of-file position.
*/
string ReadLine(FILE *infile);
#endif

@ -1,80 +0,0 @@
#include <stdlib.h>
#include <string.h>
#include "genlib.h"
#include "geostat.h"
#ifndef _TOOLSFFTMA_H
#define _TOOLSFFTMA_H
/* Create december, the 16th 2002 */
/* Modified december, the 9th 2002 */
/* List of functions: */
/* ------------------ */
/* kgeneration, FFTMA2, prebuild_gwn, build_real, addstat2, clean_real */
/*FUNCTIONS*/
/*----------*/
void kgeneration(long seed,struct grid_mod grid,struct statistic_mod stat,struct vario_mod variogram,string filename[8],struct realization_mod *Z,struct realization_mod *Y,struct realization_mod *Y2, int n[3], int *genere, int *gwnwrite, struct realization_mod *gwnoise);
/*FAST FOURIER TRANSFORM Pressure Simulation */
/*Turns a Gaussian white noise vector into a */
/*spatially correlated vector */
/*input: */
/*variogram: structure defining the variogram */
/* model */
/*grid: structure defining the grid */
/*n: vector with the number of cells along the */
/* X, Y and Z axes for the underlying grid */
/* i = [0 1 2] */
/* --> 0 0 0 : n will be computed and */
/* updated as output */
/* --> nx ny nz: these dimensions are used */
/*realin: structure defining a realization - */
/* must be a Gaussian white noise */
/*gradient: macroscopic gradient pression vector */
/*output: */
/*realout: structure defining a realization - */
/*realout2: structure defining a pressure field */
/*realout3: structure defining a xvelocity field */
/*realout4: structure defining a yvelocity field */
/*realout5: structure defining a zvelocity field */
void FFTMA2(struct vario_mod variogram,struct grid_mod grid,int n[3],struct realization_mod *realin,struct realization_mod *realout);
/* prebuild_gwn */
/* Produce a first construction in real space of the Gaussian white noise */
/* grid: structure defining the grid */
/*n: vector with the number of cells along the */
/* X, Y and Z axes for the underlying grid */
/* i = [0 1 2] */
/* --> 0 0 0 : n will be computed and */
/* updated as output */
/* --> nx ny nz: these dimensions are used */
/*realin: structure defining a realization - */
/* 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);
/* build_real */
/* build a realization in the spectral domain */
/*n: vector with the number of cells along the */
/* X, Y and Z axes for the underlying grid */
/* i = [0 1 2] */
/* --> 0 0 0 : n will be computed and */
/* updated as output */
/* --> nx ny nz: these dimensions are used */
/*covar: vector defining the covariance in spectral domain */
/*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 addstat2(struct realization_mod *realin,struct statistic_mod stat ,struct realization_mod *realout,struct realization_mod *realout2);
void clean_real(struct realization_mod *realin,int n[3],struct grid_mod grid,double *vectorresult,struct realization_mod *realout);
#endif // define _TOOLSFFTMA_H

@ -1,91 +0,0 @@
#include <stdlib.h>
#include <string.h>
#include "genlib.h"
#include "geostat.h"
#include "pressure.h"
#ifndef _TOOLSIO_H
#define _TOOLSIO_H
/* Create december, the 16th 2002 */
/* Modified december, the 9th 2002 */
/* List of functions: */
/* ------------------ */
/* writefile,writefile_bin,readfile_bin,inputdata,inputfiledata,readdata,debuginput */
/* testmemory, testopenfile */
/*FUNCTIONS*/
/*----------*/
/* Lecture dans un fichier */
/* Inputdata */
/* */
/* input data needed for realizations*/
/* seed: seed */
/* grid: structure defining the actual grid */
/* filename: name of files */
/* variogram: struture defining the variogram model */
/* stat: struture defining the mean and the variance of permeability realization */
/* pression: structure defining the gradient pressure */
void inputdata(long *seed,struct grid_mod *grid,string filename[7],struct vario_mod *variogram,struct statistic_mod *stat,struct pressure_mod *pression);
void inputfiledata(string inputfile,long *seed,struct grid_mod *grid,string filename[7],struct vario_mod *variogram,struct statistic_mod *stat,struct pressure_mod *pression);
/* void readdata(long *seed,struct grid_mod *grid,string filename[8],struct vario_mod *variogram,struct statistic_mod *stat,struct pressure_mod *pression,int *Ksolver,struct realization_mod *Kfield, char *argv[]); */
void readdata(long *seed,struct grid_mod *grid,string filename[8],struct vario_mod *variogram,struct statistic_mod *stat,struct pressure_mod *pression,int *Ksolver,int *genere, int *gwnwrite, struct realization_mod *Kfield,struct realization_mod *gwnoise, char *argv[]);
void readdata2(long *seed,struct grid_mod *grid,string filename[8],struct vario_mod *variogram,struct statistic_mod *stat,struct pressure_mod *pression,int *Ksolver,int *genere, int *gwnwrite, int *format_file, struct realization_mod *Kfield,struct realization_mod *gwnoise, char *argv[]);
/* readfile_bin */
/* */
/* read in the file "filename" the vector values of a */
/* realization_mod variable */
/* filename: explicit */
/* realin: structure defining a realization */
void readfile_bin(string filename, struct realization_mod *realin);
/* Writefile */
/* */
/* write in the file "filename" the vector values of a */
/* realization_mod variable */
/* filename: explicit */
/* realin: structure defining a realization */
void writefile( string filename, struct realization_mod *realin);
/* Writefile_bin */
/* */
/* write in the file "filename" the vector values of a */
/* realization_mod variable */
/* filename: explicit */
/* realin: structure defining a realization */
void writefile_bin( string filename, struct realization_mod *realin);
/* DebugInput */
/* */
/* Display the input data */
/* seed: seed */
/* grid: structure defining the actual grid */
/* filename: name of files */
/* variogram: struture defining the variogram model */
/* stat: struture defining the mean and the variance of permeability realization */
/* pression: structure defining the gradient pressure */
void debuginput(long *seed,struct grid_mod *grid,string filename[7],struct vario_mod *variogram,struct statistic_mod *stat,struct pressure_mod *pression,int *Ksolver,int *genere, int *gwnwrite);
/* Allocation test */
void testmemory(double *realint);
/* Test open file */
void testopenfile(FILE *fp);
#endif // define _TOOLSIO_H

@ -1,13 +0,0 @@
CC = cc
CCFLAG =
LIB = libCS106X_${ARCH}.a
%.o: %.c
$(CC) $(CCFLAG) -c $<
$(LIB) : genlib.o random.o simpio.o strlib.o symtab.o scanadt.o stack.o
ar r $@ $?
clean :
rm *.o

@ -1,21 +0,0 @@
# CC = gcc
# CCFLAG =
ifeq (${ARCH},sun)
CC = cc
CCFLAG = -g
endif
ifeq (${ARCH},linux)
CC = gcc
CCFLAG = -g
endif
LIB = libCS106X_debug_${ARCH}.a
%.o: %.c
$(CC) $(CCFLAG) -c $<
$(LIB) : genlib.o random.o simpio.o strlib.o symtab.o scanadt.o stack.o
ar r $@ $?
clean :
rm *.o

@ -1,63 +0,0 @@
/*
* File: genlib.c
* Version: 1.0
* Last modified on Fri Jul 15 15:45:52 1994 by eroberts
* -----------------------------------------------------
* This file implements the general C library package. See the
* interface description in genlib.h for details.
*/
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include "genlib.h"
/*
* Constants:
* ----------
* ErrorExitStatus -- Status value used in exit call
*/
#define ErrorExitStatus 1
/* Section 1 -- Define new "primitive" types */
/*
* Constant: UNDEFINED
* -------------------
* This entry defines the target of the UNDEFINED constant.
*/
char undefined_object[] = "UNDEFINED";
/* Section 2 -- Memory allocation */
void *GetBlock(size_t nbytes)
{
void *result;
result = malloc(nbytes);
if (result == NULL) Error("No memory available");
return (result);
}
void FreeBlock(void *ptr)
{
free(ptr);
}
/* Section 3 -- Basic error handling */
void Error(string msg, ...)
{
va_list args;
va_start(args, msg);
fprintf(stderr, "Error: ");
vfprintf(stderr, msg, args);
fprintf(stderr, "\n");
va_end(args);
exit(ErrorExitStatus);
}

@ -1,163 +0,0 @@
/*
* File: genlib.h
* Version: 1.0
* Last modified on Fri Jul 15 15:45:52 1994 by eroberts
* -----------------------------------------------------
* This file contains several definitions that form the
* core of a general-purpose ANSI C library developed by Eric
* Roberts. The goal of this library is to provide a basic
* set of tools and conventions that increase the readability
* of C programs, particularly as they are used in a teaching
* environment.
*
* The basic definitions provided by genlib.h are:
*
* 1. Declarations for several new "primitive" types
* (most importantly bool and string) that are
* used throughout the other libraries and
* applications as fundamental types.
*
* 2. A new set of functions for memory allocation.
*
* 3. A function for error handling.
*/
#ifndef _genlib_h
#define _genlib_h
#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
/* Section 1 -- Define new "primitive" types */
/*
* Type: bool
* ----------
* This type has two values, FALSE and TRUE, which are equal to 0
* and 1, respectively. Most of the advantage of defining this type
* comes from readability because it allows the programmer to
* provide documentation that a variable will take on only one of
* these two values. Designing a portable representation, however,
* is surprisingly hard, because many libraries and some compilers
* define these names. The definitions are usually compatible but
* may still be flagged as errors.
*/
#ifdef THINK_C
typedef int bool;
#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
#endif
/*
* Type: string
* ------------
* The type string is identical to the type char *, which is
* traditionally used in C programs. The main point of defining a
* new type is to improve program readability. At the abstraction
* levels at which the type string is used, it is usually not
* important to take the string apart into its component characters.
* Declaring it as a string emphasizes this atomicity.
*/
typedef char *string;
/*
* Constant: UNDEFINED
* -------------------
* Besides NULL, the only other constant of pointer type is
* UNDEFINED, which is used in certain packages as a special
* sentinel to indicate an undefined pointer value. In many
* such contexts, NULL is a legitimate data value and is
* therefore inappropriate as a sentinel.
*/
#define UNDEFINED ((void *) undefined_object)
extern char undefined_object[];
/* Section 2 -- Memory allocation */
/*
* General notes:
* --------------
* These functions provide a common interface for memory
* allocation. All functions in the library that allocate
* memory do so using GetBlock and FreeBlock. Even though
* the ANSI standard defines malloc and free for the same
* purpose, using GetBlock and FreeBlock provides greater
* compatibility with non-ANSI implementations, automatic
* out-of-memory error detection, and the possibility of
* substituting a garbage-collecting allocator.
*/
/*
* Function: GetBlock
* Usage: ptr = (type) GetBlock(nbytes);
* -------------------------------------
* GetBlock allocates a block of memory of the given size. If
* no memory is available, GetBlock generates an error.
*/
void *GetBlock(size_t nbytes);
/*
* Function: FreeBlock
* Usage: FreeBlock(ptr);
* ----------------------
* FreeBlock frees the memory associated with ptr, which must
* have been allocated using GetBlock, New, or NewArray.
*/
void FreeBlock(void *ptr);
/*
* Macro: New
* Usage: p = New(pointer-type);
* -----------------------------
* The New pseudofunction allocates enough space to hold an
* object of the type to which pointer-type points and returns
* a pointer to the newly allocated pointer. Note that
* "New" is different from the "new" operator used in C++;
* the former takes a pointer type and the latter takes the
* target type.
*/
#define New(type) ((type) GetBlock(sizeof *((type) NULL)))
/*
* Macro: NewArray
* Usage: p = NewArray(n, element-type);
* -------------------------------------
* NewArray allocates enough space to hold an array of n
* values of the specified element type.
*/
#define NewArray(n, type) ((type *) GetBlock((n) * sizeof(type)))
/* Section 3 -- Basic error handling */
/*
* Function: Error
* Usage: Error(msg, ...)
* ----------------------
* Error generates an error string, expanding % constructions
* appearing in the error message string just as printf does.
* After printing the error message, the program terminates.
*/
void Error(string msg, ...);
#endif

@ -1,77 +0,0 @@
/*
* File: random.c
* Version: 1.0
* Last modified on Mon Sep 13 10:42:45 1993 by eroberts
* -----------------------------------------------------
* This file implements the random.h interface.
*/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "genlib.h"
#include "random.h"
/*
* Function: Randomize
* -------------------
* This function operates by setting the random number
* seed to the current time. The srand function is
* provided by the <stdlib.h> library and requires an
* integer argument. The time function is provided
* by <time.h>.
*/
void Randomize(void)
{
srand((int) time(NULL));
}
/*
* Function: RandomInteger
* -----------------------
* This function first obtains a random integer in
* the range [0..RAND_MAX] by applying four steps:
* (1) Generate a real number between 0 and 1.
* (2) Scale it to the appropriate range size.
* (3) Truncate the value to an integer.
* (4) Translate it to the appropriate starting point.
*/
int RandomInteger(int low, int high)
{
int k;
double d;
d = (double) rand() / ((double) RAND_MAX + 1);
k = (int) (d * (high - low + 1));
return (low + k);
}
/*
* Function: RandomReal
* --------------------
* The implementation of RandomReal is similar to that
* of RandomInteger, without the truncation step.
*/
double RandomReal(double low, double high)
{
double d;
d = (double) rand() / ((double) RAND_MAX + 1);
return (low + d * (high - low));
}
/*
* Function: RandomChance
* ----------------------
* This function uses RandomReal to generate a number
* between 0 and 100, which it then compares to p.
*/
bool RandomChance(double p)
{
return (RandomReal(0, 1) < p);
}

@ -1,73 +0,0 @@
/*
* File: random.h
* Version: 1.0
* Last modified on Fri Jul 22 16:44:36 1994 by eroberts
* -----------------------------------------------------
* This interface provides several functions for generating
* pseudo-random numbers.
*/
#ifndef _random_h
#define _random_h
#include "genlib.h"
#include <stdlib.h>
/*
* Constant: RAND_MAX
* ------------------
* Unfortunately, several libraries that supposedly conform to
* the ANSI standard do not define RAND_MAX in <stdlib.h>. To
* reduce portability problems, this interface defines RAND_MAX
* to be the largest positive integer if it is undefined.
*/
#ifndef RAND_MAX
# define RAND_MAX ((int) ((unsigned) ~0 >> 1))
#endif
/*
* Function: Randomize
* Usage: Randomize();
* -------------------
* This function sets the random seed so that the random sequence
* is unpredictable. During the debugging phase, it is best not
* to call this function, so that program behavior is repeatable.
*/
void Randomize(void);
/*
* Function: RandomInteger
* Usage: n = RandomInteger(low, high);
* ------------------------------------
* This function returns a random integer in the range low to high,
* inclusive.
*/
int RandomInteger(int low, int high);
/*
* Function: RandomReal
* Usage: d = RandomReal(low, high);
* ---------------------------------
* This function returns a random real number in the half-open
* interval [low .. high), meaning that the result is always
* greater than or equal to low but strictly less than high.
*/
double RandomReal(double low, double high);
/*
* Function: RandomChance
* Usage: if (RandomChance(p)) . . .
* ---------------------------------
* The RandomChance function returns TRUE with the probability
* indicated by p, which should be a floating-point number between
* 0 (meaning never) and 1 (meaning always). For example, calling
* RandomChance(.30) returns TRUE 30 percent of the time.
*/
bool RandomChance(double p);
#endif

@ -1,162 +0,0 @@
/*
* File: scanadt.c
* ---------------
* This file implements the scanadt.h interface.
*/
#include <stdio.h>
#include <ctype.h>
#include "genlib.h"
#include "strlib.h"
#include "scanadt.h"
/*
* Type: scannerCDT
* ----------------
* This structure is the concrete representation of the type
* scannerADT, which is exported by this interface. Its purpose
* is to maintain the state of the scanner between calls. The
* details of the representation are invisible to the client,
* but consist of the following fields:
*
* str -- Copy of string passed to SetScannerString
* len -- Length of string, saved for efficiency
* cp -- Current character position in the string
* savedToken -- String saved by SaveToken (NULL indicates none)
* spaceOption -- Setting of the space option extension
*/
struct scannerCDT {
string str;
int len;
int cp;
string savedToken;
spaceOptionT spaceOption;
};
/* Private function prototypes */
static void SkipSpaces(scannerADT scanner);
static int ScanToEndOfIdentifier(scannerADT scanner);
/* Exported entries */
scannerADT NewScanner(void)
{
scannerADT scanner;
scanner = New(scannerADT);
scanner->str = NULL;
scanner->spaceOption = PreserveSpaces;
return (scanner);
}
void FreeScanner(scannerADT scanner)
{
if (scanner->str != NULL) FreeBlock(scanner->str);
FreeBlock(scanner);
}
void SetScannerString(scannerADT scanner, string str)
{
if (scanner->str != NULL) FreeBlock(scanner->str);
scanner->str = CopyString(str);
scanner->len = StringLength(str);
scanner->cp = 0;
scanner->savedToken = NULL;
}
string ReadToken(scannerADT scanner)
{
char ch;
string token;
int start, finish;
if (scanner->str == NULL) {
Error("SetScannerString has not been called");
}
if (scanner->savedToken != NULL) {
token = scanner->savedToken;
scanner->savedToken = NULL;
return (token);
}
if (scanner->spaceOption == IgnoreSpaces) SkipSpaces(scanner);
start = finish = scanner->cp;
if (start >= scanner->len) return (CopyString(""));
ch = scanner->str[scanner->cp];
if (isalnum(ch)) {
finish = ScanToEndOfIdentifier(scanner);
} else {
scanner->cp++;
}
return (SubString(scanner->str, start, finish));
}
bool MoreTokensExist(scannerADT scanner)
{
if (scanner->str == NULL) {
Error("SetScannerString has not been called");
}
if (scanner->savedToken != NULL) {
return (!StringEqual(scanner->savedToken, ""));
}
if (scanner->spaceOption == IgnoreSpaces) SkipSpaces(scanner);
return (scanner->cp < scanner->len);
}
void SaveToken(scannerADT scanner, string token)
{
if (scanner->str == NULL) {
Error("SetScannerString has not been called");
}
if (scanner->savedToken != NULL) {
Error("Token has already been saved");
}
scanner->savedToken = token;
}
void SetScannerSpaceOption(scannerADT scanner, spaceOptionT option)
{
scanner->spaceOption = option;
}
spaceOptionT GetScannerSpaceOption(scannerADT scanner)
{
return (scanner->spaceOption);
}
/* Private functions */
/*
* Function: SkipSpaces
* Usage: SkipSpaces(scanner);
* ---------------------------
* This function advances the position of the scanner until the
* current character is not a whitespace character.
*/
static void SkipSpaces(scannerADT scanner)
{
while (isspace(scanner->str[scanner->cp])) {
scanner->cp++;
}
}
/*
* Function: ScanToEndOfIdentifier
* Usage: finish = ScanToEndOfIdentifier(scanner);
* -----------------------------------------------
* This function advances the position of the scanner until it
* reaches the end of a sequence of letters or digits that make
* up an identifier. The return value is the index of the last
* character in the identifier; the value of the stored index
* cp is the first character after that.
*/
static int ScanToEndOfIdentifier(scannerADT scanner)
{
while (isalnum(scanner->str[scanner->cp])) {
scanner->cp++;
}
return (scanner->cp - 1);
}

@ -1,159 +0,0 @@
/*
* File: scanadt.h
* ---------------
* This file is the interface to a module that exports an abstract
* data type to facilitate dividing a string into logical units
* called "tokens," which are either
*
* 1. Strings of consecutive letters and digits representing words
* 2. One-character strings representing punctuation or separators
*
* To use this package, you must first create an instance of a
* scannerADT by calling
*
* scanner = NewScanner();
*
* All other calls to the scanner package take this variable as their
* first argument to identify a particular instance of the abstract
* scanner type.
*
* You initialize the scanner to hold a particular string by calling
*
* SetScannerString(scanner, str);
*
* where str is the string from which tokens should be read. To
* retrieve each individual token, you make the following call:
*
* token = ReadToken(scanner);
*
* To determine whether any tokens remain to be read, you can call
* the predicate function MoreTokensExist(scanner). The ReadToken
* function returns the empty string after the last token is read.
*
* The following code fragment serves as an idiom for processing
* each token in the string inputString:
*
* scanner = NewScanner();
* SetScannerString(scanner, inputString);
* while (MoreTokensExist(scanner)) {
* token = ReadToken(scanner);
* . . . process the token . . .
* }
*
* This version of scanadt.h also supports the following extensions,
* which are documented later in the interface:
*
* SaveToken
* SetScannerSpaceOption
*/
#ifndef _scanadt_h
#define _scanadt_h
#include "genlib.h"
/*
* Type: scannerADT
* ----------------
* This type is the abstract type used to represent a single instance
* of a scanner. As with any abstract type, the details of the
* internal representation are hidden from the client.
*/
typedef struct scannerCDT *scannerADT;
/*
* Function: NewScanner
* Usage: scanner = NewScanner();
* ------------------------------
* This function creates a new scanner instance. All other functions
* in this interface take this scanner value as their first argument
* so that they can identify what particular instance of the scanner
* is in use. This design makes it possible for clients to have more
* than one scanner process active at the same time.
*/
scannerADT NewScanner(void);
/*
* Function: FreeScanner
* Usage: FreeScanner(scanner);
* ----------------------------
* This function frees the storage associated with scanner.
*/
void FreeScanner(scannerADT scanner);
/*
* Function: SetScannerString
* Usage: SetScannerString(scanner, str);
* --------------------------------------
* This function initializes the scanner so that it will start
* extracting tokens from the string str.
*/
void SetScannerString(scannerADT scanner, string str);
/*
* Function: ReadToken
* Usage: token = ReadToken(scanner);
* ----------------------------------
* This function returns the next token from scanner. If
* ReadToken is called when no tokens are available, it returns
* the empty string. The token returned by ReadToken is always
* allocated in the heap, which means that clients can call
* FreeBlock when the token is no longer needed.
*/
string ReadToken(scannerADT scanner);
/*
* Function: MoreTokensExist
* Usage: if (MoreTokensExist(scanner)) . . .
* ------------------------------------------
* This function returns TRUE as long as there are additional
* tokens for the scanner to read.
*/
bool MoreTokensExist(scannerADT scanner);
/*
* Function: SaveToken
* Usage: SaveToken(scanner, token);
* ---------------------------------
* This function stores the token in the scanner data structure
* in such a way that the next time ReadToken is called, it will
* return that token without reading any additional characters
* from the input.
*/
void SaveToken(scannerADT scanner, string token);
/*
* Functions: SetScannerSpaceOption, GetScannerSpaceOption
* Usage: SetScannerSpaceOption(scanner, option);
* option = GetScannerSpaceOption(scanner);
* -----------------------------------------------
* The SetScannerSpaceOption function controls whether the scanner
* ignores whitespace characters or treats them as valid tokens.
* By default, the ReadToken function treats whitespace characters,
* such as spaces and tabs, just like any other punctuation mark.
* If, however, you call
*
* SetScannerSpaceOption(scanner, IgnoreSpaces);
*
* the scanner will skip over any white space before reading a token.
* You can restore the original behavior by calling
*
* SetScannerSpaceOption(scanner, PreserveSpaces);
*
* The GetScannerSpaceOption function returns the current setting
* of this option.
*/
typedef enum { PreserveSpaces, IgnoreSpaces } spaceOptionT;
void SetScannerSpaceOption(scannerADT scanner, spaceOptionT option);
spaceOptionT GetScannerSpaceOption(scannerADT scanner);
#endif

@ -1,157 +0,0 @@
/*
* File: simpio.c
* Version: 1.0
* Last modified on Fri Jul 15 14:10:41 1994 by eroberts
* -----------------------------------------------------
* This file implements the simpio.h interface.
*/
#include <stdio.h>
#include <string.h>
#include "genlib.h"
#include "strlib.h"
#include "simpio.h"
/*
* Constants:
* ----------
* InitialBufferSize -- Initial buffer size for ReadLine
*/
#define InitialBufferSize 120
/* Exported entries */
/*
* Functions: GetInteger, GetLong, GetReal
* ---------------------------------------
* These functions first read a line and then call sscanf to
* translate the number. Reading an entire line is essential to
* good error recovery, because the characters after the point of
* error would otherwise remain in the input buffer and confuse
* subsequent input operations. The sscanf line allows white space
* before and after the number but no other extraneous characters.
*/
int GetInteger(void)
{
string line;
int value;
char termch;
while (TRUE) {
line = GetLine();
switch (sscanf(line, " %d %c", &value, &termch)) {
case 1:
FreeBlock(line);
return (value);
case 2:
printf("Unexpected character: '%c'\n", termch);
break;
default:
printf("Please enter an integer\n");
break;
}
FreeBlock(line);
printf("Retry: ");
}
}
long GetLong(void)
{
string line;
long value;
char termch;
while (TRUE) {
line = GetLine();
switch (sscanf(line, " %ld %c", &value, &termch)) {
case 1:
FreeBlock(line);
return (value);
case 2:
printf("Unexpected character: '%c'\n", termch);
break;
default:
printf("Please enter an integer\n");
break;
}
FreeBlock(line);
printf("Retry: ");
}
}
double GetReal(void)
{
string line;
double value;
char termch;
while (TRUE) {
line = GetLine();
switch (sscanf(line, " %lf %c", &value, &termch)) {
case 1:
FreeBlock(line);
return (value);
case 2:
printf("Unexpected character: '%c'\n", termch);
break;
default:
printf("Please enter a real number\n");
break;
}
FreeBlock(line);
printf("Retry: ");
}
}
/*
* Function: GetLine
* -----------------
* This function is a simple wrapper; all the work is done by
* ReadLine.
*/
string GetLine(void)
{
return (ReadLine(stdin));
}
/*
* Function: ReadLine
* ------------------
* This function operates by reading characters from the file
* into a dynamically allocated buffer. If the buffer becomes
* full before the end of the line is reached, a new buffer
* twice the size of the previous one is allocated.
*/
string ReadLine(FILE *infile)
{
string line, nline;
int n, ch, size;
n = 0;
size = InitialBufferSize;
line = GetBlock(size + 1);
while ((ch = getc(infile)) != '\n' && ch != EOF) {
if (n == size) {
size *= 2;
nline = (string) GetBlock(size + 1);
strncpy(nline, line, n);
FreeBlock(line);
line = nline;
}
line[n++] = ch;
}
if (n == 0 && ch == EOF) {
FreeBlock(line);
return (NULL);
}
line[n] = '\0';
nline = (string) GetBlock(n + 1);
strcpy(nline, line);
FreeBlock(line);
return (nline);
}

@ -1,75 +0,0 @@
/*
* File: simpio.h
* Version: 1.0
* Last modified on Wed Apr 27 07:29:13 1994 by eroberts
* -----------------------------------------------------
* This interface provides access to a simple package of
* functions that simplify the reading of input data.
*/
#ifndef _simpio_h
#define _simpio_h
#include "genlib.h"
/*
* Function: GetInteger
* Usage: i = GetInteger();
* ------------------------
* GetInteger reads a line of text from standard input and scans
* it as an integer. The integer value is returned. If an
* integer cannot be scanned or if more characters follow the
* number, the user is given a chance to retry.
*/
int GetInteger(void);
/*
* Function: GetLong
* Usage: l = GetLong();
* ---------------------
* GetLong reads a line of text from standard input and scans
* it as a long integer. The value is returned as a long.
* If an integer cannot be scanned or if more characters follow
* the number, the user is given a chance to retry.
*/
long GetLong(void);
/*
* Function: GetReal
* Usage: x = GetReal();
* ---------------------
* GetReal reads a line of text from standard input and scans
* it as a double. If the number cannot be scanned or if extra
* characters follow after the number ends, the user is given
* a chance to reenter the value.
*/
double GetReal(void);
/*
* Function: GetLine
* Usage: s = GetLine();
* ---------------------
* GetLine reads a line of text from standard input and returns
* the line as a string. The newline character that terminates
* the input is not stored as part of the string.
*/
string GetLine(void);
/*
* Function: ReadLine
* Usage: s = ReadLine(infile);
* ----------------------------
* ReadLine reads a line of text from the input file and
* returns the line as a string. The newline character
* that terminates the input is not stored as part of the
* string. The ReadLine function returns NULL if infile
* is at the end-of-file position.
*/
string ReadLine(FILE *infile);
#endif

@ -1,122 +0,0 @@
/*
* File: stack.c
* -------------
* This file implements the stack.h interface. Note that the
* implementation is independent of the type of value stored
* in the stack. That type is defined by the interface as
* the type name stackElementT.
*/
#include <stdio.h>
#include "genlib.h"
#include "stack.h"
/*
* Constant: InitialStackSize
* --------------------------
* This constant defines the initial size of the stack array.
* Any positive value will work correctly, although changing
* this parameter can affect performance. Making this value
* larger postpones the first reallocation but causes stacks
* to consume more memory.
*/
#define InitialStackSize 100
/*
* Type: stackCDT
* --------------
* The type stackCDT is the concrete representation of the type
* stackADT defined by the interface. In this implementation,
* the elements are stored in a dynamic array that doubles in
* size if the old stack becomes full.
*/
struct stackCDT {
stackElementT *elements;
int count;
int size;
};
/* Private function prototypes */
static void ExpandStack(stackADT stack);
/* Exported entries */
stackADT NewStack(void)
{
stackADT stack;
stack = New(stackADT);
stack->elements = NewArray(InitialStackSize, stackElementT);
stack->count = 0;
stack->size = InitialStackSize;
return (stack);
}
void FreeStack(stackADT stack)
{
FreeBlock(stack->elements);
FreeBlock(stack);
}
void Push(stackADT stack, stackElementT element)
{
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");
return (stack->elements[--stack->count]);
}
bool StackIsEmpty(stackADT stack)
{
return (stack->count == 0);
}
bool StackIsFull(stackADT stack)
{
return (FALSE);
}
int StackDepth(stackADT stack)
{
return (stack->count);
}
stackElementT GetStackElement(stackADT stack, int index)
{
if (index < 0 || index >= stack->count) {
Error("Non-existent stack element");
}
return (stack->elements[stack->count - index - 1]);
}
/* Private functions */
/* Function: ExpandStack
* Usage: ExpandStack(stack);
* --------------------------
* This function expands a full stack by doubling the size of its
* dynamic array, copying the old elements to the new array, and
* then freeing the old array storage.
*/
static void ExpandStack(stackADT stack)
{
stackElementT *array;
int i, newSize;
newSize = stack->size * 2;
array = NewArray(newSize, stackElementT);
for (i = 0; i < stack->size; i++) {
array[i] = stack->elements[i];
}
FreeBlock(stack->elements);
stack->elements = array;
stack->size = newSize;
}

@ -1,115 +0,0 @@
/*
* File: stack.h
* -------------
* This interface defines an abstraction for stacks. In any
* single application that uses this interface, the values in
* the stack are constrained to a single type, although it
* is easy to change that type by changing the definition of
* stackElementT in this interface.
*/
#ifndef _stack_h
#define _stack_h
#include "genlib.h"
/*
* Type: stackElementT
* -------------------
* The type stackElementT is used in this interface to indicate
* the type of values that can be stored in the stack. Here the
* stack is used to store values of type double, but that can
* be changed by editing this definition line.
*/
typedef void *stackElementT;
/*
* Type: stackADT
* --------------
* The type stackADT represents the abstract type used to store
* the elements that have been pushed. Because stackADT is
* defined only as a pointer to a concrete structure that is not
* itself defined in the interface, clients have no access to
* the underlying fields.
*/
typedef struct stackCDT *stackADT;
/*
* Function: NewStack
* Usage: stack = NewStack();
* --------------------------
* This function allocates and returns a new stack, which is
* initially empty.
*/
stackADT NewStack(void);
/*
* Function: FreeStack
* Usage: FreeStack(stack);
* ------------------------
* This function frees the storage associated with the stack.
*/
void FreeStack(stackADT stack);
/*
* Function: Push
* Usage: Push(stack, element);
* ----------------------------
* This function pushes the specified element onto the stack.
*/
void Push(stackADT stack, stackElementT element);
/*
* Function: Pop
* Usage: element = Pop(stack);
* ----------------------------
* This function pops the top element from the stack and returns
* that value. The first value popped is always the last one
* that was pushed. If the stack is empty when Pop is called,
* the function calls Error with an appropriate message.
*/
stackElementT Pop(stackADT stack);
/*
* Functions: StackIsEmpty, StackIsFull
* Usage: if (StackIsEmpty(stack)) . . .
* if (StackIsFull(stack)) . . .
* -------------------------------------
* This functions test whether the stack is empty or full.
*/
bool StackIsEmpty(stackADT stack);
bool StackIsFull(stackADT stack);
/*
* Function: StackDepth
* Usage: depth = StackDepth(stack);
* ---------------------------------
* This function returns the number of elements currently pushed
* on the stack.
*/
int StackDepth(stackADT stack);
/*
* Function: GetStackElement
* Usage: element = GetStackElement(stack, index);
* -----------------------------------------------
* This function returns the element at the specified index in
* the stack, where the top of the stack is defined as index 0.
* For example, calling GetStackElement(stack, 0) returns the top
* element on the stack without removing it. If the caller tries
* to select an out-of-range element, GetStackElement calls Error.
* Note: This function is not a fundamental stack operation and
* is instead provided principally to facilitate debugging.
*/
stackElementT GetStackElement(stackADT stack, int index);
#endif

@ -1,246 +0,0 @@
/*
* File: strlib.c
* Version: 1.0
* Last modified on Fri Jul 15 14:10:41 1994 by eroberts
* -----------------------------------------------------
* This file implements the strlib.h interface.
*/
/*
* General implementation notes:
* -----------------------------
* This module implements the strlib library by mapping all
* functions into the appropriate calls to the ANSI <string.h>
* interface. The implementations of the individual functions
* are all quite simple and do not require individual comments.
* For descriptions of the behavior of each function, see the
* interface.
*/
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include "genlib.h"
#include "strlib.h"
/*
* Constant: MaxDigits
* -------------------
* This constant must be larger than the maximum
* number of digits that can appear in a number.
*/
#define MaxDigits 30
/* Private function prototypes */
static string CreateString(int len);
/* Section 1 -- Basic string operations */
string Concat(string s1, string s2)
{
string s;
int len1, len2;
if (s1 == NULL || s2 == NULL) {
Error("NULL string passed to Concat");
}
len1 = strlen(s1);
len2 = strlen(s2);
s = CreateString(len1 + len2);
strcpy(s, s1);
strcpy(s + len1, s2);
return (s);
}
char IthChar(string s, int i)
{
int len;
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");
}
return (s[i]);
}
string SubString(string s, int p1, int p2)
{
int len;
string result;
if (s == NULL) Error("NULL string passed to SubString");
len = strlen(s);
if (p1 < 0) p1 = 0;
if (p2 >= len) p2 = len - 1;
len = p2 - p1 + 1;
if (len < 0) len = 0;
result = CreateString(len);
strncpy(result, s + p1, len);
result[len] = '\0';
return (result);
}
string CharToString(char ch)
{
string result;
result = CreateString(1);
result[0] = ch;
result[1] = '\0';
return (result);
}
int StringLength(string s)
{
if (s == NULL) Error("NULL string passed to StringLength");
return (strlen(s));
}
string CopyString(string s)
{
string newstr;
if (s == NULL) Error("NULL string passed to CopyString");
newstr = CreateString(strlen(s));
strcpy(newstr, s);
return (newstr);
}
/* Section 2 -- String comparison functions */
bool StringEqual(string s1, string s2)
{
if (s1 == NULL || s2 == NULL) {
Error("NULL string passed to StringEqual");
}
return (strcmp(s1, s2) == 0);
}
int StringCompare(string s1, string s2)
{
if (s1 == NULL || s2 == NULL) {
Error("NULL string passed to StringCompare");
}
return (strcmp(s1, s2));
}
/* Section 3 -- Search functions */
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);
cptr = strchr(text + start, ch);
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);
cptr = strstr(text + start, str);
if (cptr == NULL) return (-1);
return ((int) (cptr - text));
}
/* Section 4 -- Case-conversion functions */
string ConvertToLowerCase(string s)
{
string result;
int i;
if (s == NULL) {
Error("NULL string passed to ConvertToLowerCase");
}
result = CreateString(strlen(s));
for (i = 0; s[i] != '\0'; i++) result[i] = tolower(s[i]);
result[i] = '\0';
return (result);
}
string ConvertToUpperCase(string s)
{
string result;
int i;
if (s == NULL) {
Error("NULL string passed to ConvertToUpperCase");
}
result = CreateString(strlen(s));
for (i = 0; s[i] != '\0'; i++) result[i] = toupper(s[i]);
result[i] = '\0';
return (result);
}
/* Section 5 -- Functions for converting numbers to strings */
string IntegerToString(int n)
{
char buffer[MaxDigits];
sprintf(buffer, "%d", n);
return (CopyString(buffer));
}
int StringToInteger(string s)
{
int result;
char dummy;
if (s == NULL) {
Error("NULL string passed to StringToInteger");
}
if (sscanf(s, " %d %c", &result, &dummy) != 1) {
Error("StringToInteger called on illegal number %s", s);
}
return (result);
}
string RealToString(double d)
{
char buffer[MaxDigits];
sprintf(buffer, "%G", d);
return (CopyString(buffer));
}
double StringToReal(string s)
{
double result;
char dummy;
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);
}
return (result);
}
/* Private functions */
/*
* Function: CreateString
* Usage: s = CreateString(len);
* -----------------------------
* This function dynamically allocates space for a string of
* len characters, leaving room for the null character at the
* end.
*/
static string CreateString(int len)
{
return ((string) GetBlock(len + 1));
}

@ -1,243 +0,0 @@
/*
* File: strlib.h
* Version: 1.0
* Last modified on Fri Jul 15 14:10:40 1994 by eroberts
* -----------------------------------------------------
* The strlib.h file defines the interface for a simple
* string library. In the context of this package, strings
* are considered to be an abstract data type, which means
* that the client relies only on the operations defined for
* the type and not on the underlying representation.
*/
/*
* Cautionary note:
* ----------------
* Although this interface provides an extremely convenient
* abstraction for working with strings, it is not appropriate
* for all applications. In this interface, the functions that
* return string values (such as Concat and SubString) do so
* by allocating new memory. Over time, a program that uses
* this package will consume increasing amounts of memory
* and eventually exhaust the available supply. If you are
* writing a program that runs for a short time and stops,
* the fact that the package consumes memory is not a problem.
* If, however, you are writing an application that must run
* for an extended period of time, using this package requires
* that you make some provision for freeing any allocated
* storage.
*/
#ifndef _strlib_h
#define _strlib_h
#include "genlib.h"
/* Section 1 -- Basic string operations */
/*
* Function: Concat
* Usage: s = Concat(s1, s2);
* --------------------------
* This function concatenates two strings by joining them end
* to end. For example, Concat("ABC", "DE") returns the string
* "ABCDE".
*/
string Concat(string s1, string s2);
/*
* Function: IthChar
* Usage: ch = IthChar(s, i);
* --------------------------
* This function returns the character at position i in the
* string s. It is included in the library to make the type
* string a true abstract type in the sense that all of the
* necessary operations can be invoked using functions. Calling
* IthChar(s, i) is like selecting s[i], except that IthChar
* checks to see if i is within the range of legal index
* positions, which extend from 0 to StringLength(s).
* IthChar(s, StringLength(s)) returns the null character
* at the end of the string.
*/
char IthChar(string s, int i);
/*
* Function: SubString
* Usage: t = SubString(s, p1, p2);
* --------------------------------
* SubString returns a copy of the substring of s consisting
* of the characters between index positions p1 and p2,
* inclusive. The following special cases apply:
*
* 1. If p1 is less than 0, it is assumed to be 0.
* 2. If p2 is greater than the index of the last string
* position, which is StringLength(s) - 1, then p2 is
* set equal to StringLength(s) - 1.
* 3. If p2 < p1, SubString returns the empty string.
*/
string SubString(string s, int p1, int p2);
/*
* Function: CharToString
* Usage: s = CharToString(ch);
* ----------------------------
* This function takes a single character and returns a
* one-character string consisting of that character. The
* CharToString function is useful, for example, if you
* need to concatenate a string and a character. Since
* Concat requires two strings, you must first convert
* the character into a string.
*/
string CharToString(char ch);
/*
* Function: StringLength
* Usage: len = StringLength(s);
* -----------------------------
* This function returns the length of s.
*/
int StringLength(string s);
/*
* Function: CopyString
* Usage: newstr = CopyString(s);
* ------------------------------
* CopyString copies the string s into dynamically allocated
* storage and returns the new string. This function is not
* ordinarily required if this package is used on its own,
* but is often necessary when you are working with more than
* one string package.
*/
string CopyString(string s);
/* Section 2 -- String comparison functions */
/*
* Function: StringEqual
* Usage: if (StringEqual(s1, s2)) ...
* -----------------------------------
* This function returns TRUE if the strings s1 and s2 are
* equal. For the strings to be considered equal, every
* character in one string must precisely match the
* corresponding character in the other. Uppercase and
* lowercase characters are considered to be different.
*/
bool StringEqual(string s1, string s2);
/*
* Function: StringCompare
* Usage: if (StringCompare(s1, s2) < 0) ...
* -----------------------------------------
* This function returns a number less than 0 if string s1
* comes before s2 in alphabetical order, 0 if they are equal,
* and a number greater than 0 if s1 comes after s2. The
* ordering is determined by the internal representation used
* for characters, which is usually ASCII.
*/
int StringCompare(string s1, string s2);
/* Section 3 -- Search functions */
/*
* Function: FindChar
* Usage: p = FindChar(ch, text, start);
* -------------------------------------
* Beginning at position start in the string text, this
* function searches for the character ch and returns the
* first index at which it appears or -1 if no match is
* found.
*/
int FindChar(char ch, string text, int start);
/*
* Function: FindString
* Usage: p = FindString(str, text, start);
* ----------------------------------------
* Beginning at position start in the string text, this
* function searches for the string str and returns the
* first index at which it appears or -1 if no match is
* found.
*/
int FindString(string str, string text, int start);
/* Section 4 -- Case-conversion functions */
/*
* Function: ConvertToLowerCase
* Usage: s = ConvertToLowerCase(s);
* ---------------------------------
* This function returns a new string with all
* alphabetic characters converted to lower case.
*/
string ConvertToLowerCase(string s);
/*
* Function: ConvertToUpperCase
* Usage: s = ConvertToUpperCase(s);
* ---------------------------------
* This function returns a new string with all
* alphabetic characters converted to upper case.
*/
string ConvertToUpperCase(string s);
/* Section 5 -- Functions for converting numbers to strings */
/*
* Function: IntegerToString
* Usage: s = IntegerToString(n);
* ------------------------------
* This function converts an integer into the corresponding
* string of digits. For example, IntegerToString(123)
* returns "123" as a string.
*/
string IntegerToString(int n);
/*
* Function: StringToInteger
* Usage: n = StringToInteger(s);
* ------------------------------
* This function converts a string of digits into an integer.
* If the string is not a legal integer or contains extraneous
* characters, StringToInteger signals an error condition.
*/
int StringToInteger(string s);
/*
* Function: RealToString
* Usage: s = RealToString(d);
* ---------------------------
* This function converts a floating-point number into the
* corresponding string form. For example, calling
* RealToString(23.45) returns "23.45". The conversion is
* the same as that used for "%G" format in printf.
*/
string RealToString(double d);
/*
* Function: StringToReal
* Usage: d = StringToReal(s);
* ---------------------------
* This function converts a string representing a real number
* into its corresponding value. If the string is not a
* legal floating-point number or if it contains extraneous
* characters, StringToReal signals an error condition.
*/
double StringToReal(string s);
#endif

@ -1,180 +0,0 @@
/*
* File: symtab.c
* --------------
* This file implements the symbol table abstraction.
*/
#include <stdio.h>
#include "genlib.h"
#include "strlib.h"
#include "symtab.h"
/*
* Constants
* ---------
* NBuckets -- Number of buckets in the hash table
*/
#define NBuckets 101
/*
* Type: cellT
* -----------
* This type defines a linked list cell for the symbol table.
*/
typedef struct cellT {
string key;
void *value;
struct cellT *link;
} cellT;
/*
* Type: symtabCDT
* ---------------
* This type defines the underlying concrete representation for a
* symtabADT. These details are not relevant to and therefore
* not exported to the client. In this implementation, the
* underlying structure is a hash table organized as an array of
* "buckets," in which each bucket is a linked list of elements
* that share the same hash code.
*/
struct symtabCDT {
cellT *buckets[NBuckets];
};
/* Private function declarations */
static void FreeBucketChain(cellT *cp);
static cellT *FindCell(cellT *cp, string s);
static int Hash(string s, int nBuckets);
/* Public entries */
symtabADT NewSymbolTable(void)
{
symtabADT table;
int i;
table = New(symtabADT);
for (i = 0; i < NBuckets; i++) {
table->buckets[i] = NULL;
}
return (table);
}
void FreeSymbolTable(symtabADT table)
{
int i;
for (i = 0; i < NBuckets; i++) {
FreeBucketChain(table->buckets[i]);
}
FreeBlock(table);
}
void Enter(symtabADT table, string key, void *value)
{
int bucket;
cellT *cp;
bucket = Hash(key, NBuckets);
cp = FindCell(table->buckets[bucket], key);
if (cp == NULL) {
cp = New(cellT *);
cp->key = CopyString(key);
cp->link = table->buckets[bucket];
table->buckets[bucket] = cp;
}
cp->value = value;
}
void *Lookup(symtabADT table, string key)
{
int bucket;
cellT *cp;
bucket = Hash(key, NBuckets);
cp = FindCell(table->buckets[bucket], key);
if (cp == NULL) return(UNDEFINED);
return (cp->value);
}
void MapSymbolTable(symtabFnT fn, symtabADT table,
void *clientData)
{
int i;
cellT *cp;
for (i = 0; i < NBuckets; i++) {
for (cp = table->buckets[i]; cp != NULL; cp = cp->link) {
fn(cp->key, cp->value, clientData);
}
}
}
/* Private functions */
/*
* Function: FreeBucketChain
* Usage: FreeBucketChain(cp);
* ---------------------------
* This function takes a chain pointer and frees all the cells
* in that chain. Because the package makes copies of the keys,
* this function must free the string storage as well.
*/
static void FreeBucketChain(cellT *cp)
{
cellT *next;
while (cp != NULL) {
next = cp->link;
FreeBlock(cp->key);
FreeBlock(cp);
cp = next;
}
}
/*
* Function: FindCell
* Usage: cp = FindCell(cp, key);
* ------------------------------
* This function finds a cell in the chain beginning at cp that
* matches key. If a match is found, a pointer to that cell is
* returned. If no match is found, the function returns NULL.
*/
static cellT *FindCell(cellT *cp, string key)
{
while (cp != NULL && !StringEqual(cp->key, key)) {
cp = cp->link;
}
return (cp);
}
/*
* Function: Hash
* Usage: bucket = Hash(key, nBuckets);
* ------------------------------------
* This function takes the key and uses it to derive a hash code,
* which is an integer in the range [0, nBuckets - 1]. The hash
* code is computed using a method called linear congruence. The
* choice of the value for Multiplier can have a significant effect
* on the performance of the algorithm, but not on its correctness.
*/
#define Multiplier -1664117991L
static int Hash(string s, int nBuckets)
{
int i;
unsigned long hashcode;
hashcode = 0;
for (i = 0; s[i] != '\0'; i++) {
hashcode = hashcode * Multiplier + s[i];
}
return (hashcode % nBuckets);
}

@ -1,85 +0,0 @@
/*
* File: symtab.h
* --------------
* This interface exports a simple symbol table abstraction.
*/
#ifndef _symtab_h
#define _symtab_h
#include "genlib.h"
/*
* Type: symtabADT
* ---------------
* This type is the ADT used to represent a symbol table.
*/
typedef struct symtabCDT *symtabADT;
/*
* Type: symtabFnT
* ---------------
* This type defines the class of functions that can be used to
* map over the entries in a symbol table.
*/
typedef void (*symtabFnT)(string key, void *value,
void *clientData);
/* Exported entries */
/*
* Function: NewSymbolTable
* Usage: table = NewSymbolTable();
* --------------------------------
* This function allocates a new symbol table with no entries.
*/
symtabADT NewSymbolTable(void);
/*
* Function: FreeSymbolTable
* Usage: FreeSymbolTable(table);
* ------------------------------
* This function frees the storage associated with the symbol table.
*/
void FreeSymbolTable(symtabADT table);
/*
* Function: Enter
* Usage: Enter(table, key, value);
* --------------------------------
* This function associates key with value in the symbol table.
* Each call to Enter supersedes any previous definition for key.
*/
void Enter(symtabADT table, string key, void *value);
/*
* Function: Lookup
* Usage: value = Lookup(table, key);
* ----------------------------------
* This function returns the value associated with key in the symbol
* table, or UNDEFINED, if no such value exists.
*/
void *Lookup(symtabADT table, string key);
/*
* Function: MapSymbolTable
* Usage: MapSymbolTable(fn, table, clientData);
* ---------------------------------------------
* This function goes through every entry in the symbol table
* and calls the function fn, passing it the following arguments:
* the current key, its associated value, and the clientData
* pointer. The clientData pointer allows the client to pass
* additional state information to the function fn, if necessary.
* If no clientData argument is required, this value should be NULL.
*/
void MapSymbolTable(symtabFnT fn, symtabADT table,
void *clientData);
#endif

@ -1,30 +0,0 @@
########################################################################
#
# Makefile for library
#
########################################################################
INTERFACE = ../include
INCLUDE = -I${INTERFACE}
LIBS = -lm -L../lib
CC = cc
CCFLAGS =
LIB = libFFTMA_${ARCH}.a
NOBJS= gammf.o fftma.o addstat.o axes.o cgrid.o covariance.o fourt.o length.o maxfactor.o test_fact.o cov_value.o generate.o gasdev.o ran2.o stable.o gaussian.o power.o cubic.o spherical.o nugget.o exponential.o cardsin.o nor2log.o
.c.o:
${CC} $(INCLUDE) $(CCFLAGS) -c $<
# LIBRARY
$(LIB) : $(NOBJS)
ar cr $(LIB) $(NOBJS)
ranlib $(LIB)
clean :
rm *.o

@ -1,106 +0,0 @@
#include <stdio.h>
#include <math.h>
#include "geostat.h"
#include "genlib.h"
/*addstat */
/*adds mean and variance effects to the N(0,1) realization*/
/*input: */
/*realin: structure defining a realization - */
/* must have zeio mean and unit variance */
/*stat: structure defining the mean and variance */
/*output: */
/*realout: structure defining a realization - */
void addstat(struct realization_mod *realin,struct statistic_mod stat ,struct realization_mod *realout)
{
int i,nblockm,nblockv;
double r,moy,var;
/*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)
Error("No memory available");
}
(*realout).n = (*realin).n;
/*test over the input realization*/
switch ((*realin).code) {
case 0:
case 1:
(*realout).code = 2;
break;
case 6:
(*realout).code = 7;
break;
default:
(*realout).code = (*realin).code;
break;
}
for (i = 0; i < (*realin).n; i++) {
/*mean*/
switch (stat.nblock_mean) {
case 1:
/*stationary case*/
nblockm = 1;
break;
default:
/*number of the cell - nonstationary case*/
nblockm = i+1;
break;
}
/*variance*/
switch (stat.nblock_var) {
case 1:
/*stationary case*/
nblockv = 1;
break;
default:
/*number of the cell - nonstationary case*/
nblockv = i+1;
break;
}
switch (stat.type) {
case 0:
/*normal case*/
moy = stat.mean[nblockm-1];
var = stat.variance[nblockv-1];
break;
case 1:
/*lognormal (natural) case*/
r = (double)sqrt(stat.variance[nblockv-1])/stat.mean[nblockm-1];
r *= r;
moy = (double)log(stat.mean[nblockm-1]/sqrt(1.0+r));
var = (double)log(1.0+r);
break;
case 2:
/*lognormal (log10) case*/
r = (double)sqrt(stat.variance[nblockv-1])/stat.mean[nblockm-1];
r *= r;
moy = (double)log10(stat.mean[nblockm-1]/sqrt(1.0+r));
var = (double)log10(1.0+r);
break;
default:
Error("Type not defined in addstat");
break;
}
(*realout).vector[i] = (double)sqrt(var)*(*realin).vector[i]+moy;
}
return;
}

@ -1,41 +0,0 @@
#include <math.h>
#include <stdlib.h>
/*normalizes anisotropy axes*/
void axes(double *ap,double *scf,int N)
{
double sclpdt, r, eps = 1.E-6;
int i,j,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];
}
}
}
}
return;
}

@ -1,18 +0,0 @@
#include <stdio.h>
#include <math.h>
#include "genlib.h"
/*cardsin covariance function*/
double cardsin(double h)
{
float delta = 20.371;
double z;
if (h != 0) {
z = (double)(h*delta);
z = sin(z)/z;
} else {
z = 1.;
}
return(z);
}

@ -1,43 +0,0 @@
#include <stdlib.h>
#include "geostat.h"
/*computes the size of the grid for FFTs*/
/*input: */
/*variogram: structure defining the variogram model*/
/*grid: structure defining the actual grid */
/*output: */
/*n: vector with the number of cells along the */
/* X, Y and Z axes for the underlying grid */
/* i = [0 1 2] */
void cgrid(struct vario_mod variogram, struct grid_mod grid, int n[3])
{
int 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]<grid.NX) || (n[1]<grid.NY) || (n[2]<grid.NZ)) {
printf("Indicated dimensions are inappropriate in cgrid");
exit;
}
}
return;
}

@ -1,55 +0,0 @@
#include <math.h>
#include "geostat.h"
#include "genlib.h"
/*selection of model covariance*/
double cov_value(struct vario_mod variogram,double di,double dj,double dk)
{
double hx,hy,hz,h;
double cov;
int k;
cov = 0.;
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;
}
}
return (cov);
}

@ -1,90 +0,0 @@
#include "geostat.h"
/*builds the sampled covariance function*/
/*dimensions are even*/
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 <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;
}

@ -1,17 +0,0 @@
#include <stdio.h>
#include <math.h>
#include "genlib.h"
/*cubic covariance function*/
double cubic(double h)
{
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);
}

@ -1,9 +0,0 @@
#include <stdio.h>
#include <math.h>
#include "genlib.h"
/*exponential covariance function*/
double exponential(double h)
{
return (exp(-3.*(double)h));
}

@ -1,185 +0,0 @@
#include <stdlib.h>
#include <math.h>
#include "geostat.h"
/*FAST FOURIER TRANSFORM MOVING AVERAGE METHOD */
/*Turns a Gaussian white noise vector into a */
/*spatially correlated vector */
/*input: */
/*variogram: structure defining the variogram */
/* model */
/*grid: structure defining the grid */
/*n: vector with the number of cells along the */
/* X, Y and Z axes for the underlying grid */
/* i = [0 1 2] */
/* --> 0 0 0 : n will be computed and */
/* updated as output */
/* --> nx ny nz: these dimensions are used */
/*realin: structure defining a realization - */
/* must be a Gaussian white noise */
/*output: */
/*realout: structure defining a realization - */
void FFTMA(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;
double *table,*covar,*workr,*worki,*realization,temp;
/*test over the input realization*/
/*if ((*realin).code != 0) {
printf("Input realizations in FFTMA must be Gaussian white noises");
exit;
}*/
/*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));
if (covar == NULL) {
printf("FFTMA.c: No memory available for covar");
exit;
}
table = (double *) malloc(ntot * sizeof(double));
if (table == NULL) {
printf("FFTMA.c: No memory available for table");
exit;
}
realization = (double *) malloc(ntot * sizeof(double));
if (realization == NULL) {
printf("FFTMA.c: No memory available for realization");
exit;
}
workr = (double *) malloc(nmax * sizeof(double));
if (workr == NULL) {
printf("FFTMA.c: No memory available for workr");
exit;
}
worki = (double *) malloc(nmax * sizeof(double));
if (worki == NULL) {
printf("FFTMA.c: No memory available for worki");
exit;
}
/*covariance function creation*/
covariance(covar,variogram,grid,n);
/*power spectrum*/
fourt(covar,table,n,NDIM,1,0,workr,worki);
/*organization of the input Gaussian white noise*/
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];
if (i <= grid.NX && j <= grid.NY && k <= grid.NZ) {
maille0 = i-1+(j-1+(k-1)*grid.NY)*grid.NX;
realization[maille1] = (*realin).vector[maille0];
} else {
realization[maille1] = 0.;
}
}
}
}
/*forward fourier transform of the GWN*/
fourt(realization,table,n,NDIM,1,0,workr,worki);
/*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;
table[maille1] *= temp;
}
}
}
free(covar);
/*backward fourier transform*/
fourt(realization,table,n,NDIM,0,1,workr,worki);
free(table);
free(workr);
free(worki);
/*output realization*/
/*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("FFTMA.c: No memory available");
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;
(*realout).vector[maille0] = realization[maille1];
}
}
}
free(realization);
return;
}

@ -1,591 +0,0 @@
#include <stdio.h>
#include <math.h>
/*fast fourier transform */
/* THE COOLEY-TUKEY FAST FOURIER TRANSFORM */
/* EVALUATES COMPLEX FOURIER SERIES FOR COMPLEX OR REAL FUNCTIONS. */
/* THAT IS, IT COMPUTES */
/* FTRAN(J1,J2,...)=SUM(DATA(I1,I2,...)*W1**(I1-1)*(J1-1) */
/* *W2**(I2-1)*(J2-1)*...), */
/* WHERE W1=EXP(-2*PI*SQRT(-1)/NN(1)), W2=EXP(-2*PI*SQRT(-1)/NN(2)), */
/* ETC. AND I1 AND J1 RUN FROM 1 TO NN(1), I2 AND J2 RUN FROM 1 TO */
/* NN(2), ETC. THERE IS NO LIMIT ON THE DIMENSIONALITY (NUMBER OF */
/* SUBSCRIPTS) OF THE ARRAY OF DATA. THE PROGRAM WILL PERFORM */
/* A THREE-DIMENSIONAL FOURIER TRANSFORM AS EASILY AS A ONE-DIMEN- */
/* SIONAL ONE, THO IN A PROPORTIONATELY GREATER TIME. AN INVERSE */
/* TRANSFORM CAN BE PERFORMED, IN WHICH THE SIGN IN THE EXPONENTIALS */
/* IS +, INSTEAD OF -. IF AN INVERSE TRANSFORM IS PERFORMED UPON */
/* AN ARRAY OF TRANSFORMED DATA, THE ORIGINAL DATA WILL REAPPEAR, */
/* MULTIPLIED BY NN(1)*NN(2)*... THE ARRAY OF INPUT DATA MAY BE */
/* REAL OR COMPLEX, AT THE PROGRAMMERS OPTION, WITH A SAVING OF */
/* ABOUT THIRTY PER CENT IN RUNNING TIME FOR REAL OVER COMPLEX. */
/* (FOR FASTEST TRANSFORM OF REAL DATA, NN(1) SHOULD BE EVEN.) */
/* THE TRANSFORM VALUES ARE ALWAYS COMPLEX, AND ARE RETURNED IN THE */
/* ORIGINAL ARRAY OF DATA, REPLACING THE INPUT DATA. THE LENGTH */
/* OF EACH DIMENSION OF THE DATA ARRAY MAY BE ANY INTEGER. THE */
/* PROGRAM RUNS FASTER ON COMPOSITE INTEGERS THAN ON PRIMES, AND IS */
/* PARTICULARLY FAST ON NUMBERS RICH IN FACTORS OF TWO. */
/* TIMING IS IN FACT GIVEN BY THE FOLLOWING FORMULA. LET NTOT BE THE */
/* TOTAL NUMBER OF POINTS (REAL OR COMPLEX) IN THE DATA ARRAY, THAT */
/* IS, NTOT=NN(1)*NN(2)*... DECOMPOSE NTOT INTO ITS PRIME FACTORS, */
/* SUCH AS 2**K2 * 3**K3 * 5**K5 * ... LET SUM2 BE THE SUM OF ALL */
/* THE FACTORS OF TWO IN NTOT, THAT IS, SUM2 = 2*K2. LET SUMF BE */
/* THE SUM OF ALL OTHER FACTORS OF NTOT, THAT IS, SUMF = 3*K3+5*K5+.. */
/* THE TIME TAKEN BY A MULTIDIMENSIONAL TRANSFORM ON THESE NTOT DATA */
/* IS T = T0 + T1*NTOT + T2*NTOT*SUM2 + T3*NTOT*SUMF. FOR THE PAR- */
/* TICULAR IMPLEMENTATION FORTRAN 32 ON THE CDC 3300 (FLOATING POINT */
/* ADD TIME = SIX MICROSECONDS), */
/* T = 3000 + 600*NTOT + 50*NTOT*SUM2 + 175*NTOT*SUMF MICROSECONDS */
/* ON COMPLEX DATA. */
/* IMPLEMENTATION OF THE DEFINITION BY SUMMATION WILL RUN IN A TIME */
/* PROPORTIONAL TO NTOT**2. FOR HIGHLY COMPOSITE NTOT, THE SAVINGS */
/* OFFERED BY COOLEY-TUKEY CAN BE DRAMATIC. A MATRIX 100 BY 100 WILL */
/* BE TRANSFORMED IN TIME PROPORTIONAL TO 10000*(2+2+2+2+5+5+5+5) = */
/* 280,000 (ASSUMING T2 AND T3 TO BE ROUGHLY COMPARABLE) VERSUS */
/* 10000**2 = 100,000,000 FOR THE STRAIGHTFORWARD TECHNIQUE. */
/* THE COOLEY-TUKEY ALGORITHM PLACES TWO RESTRICTIONS UPON THE */
/* NATURE OF THE DATA BEYOND THE USUAL RESTRICTION THAT */
/* THE DATA FROM ONE CYCLE OF A PERIODIC FUNCTION. THEY ARE-- */
/* 1. THE NUMBER OF INPUT DATA AND THE NUMBER OF TRANSFORM VALUES */
/* MUST BE THE SAME. */
/* 2. CONSIDERING THE DATA TO BE IN THE TIME DOMAIN, */
/* THEY MUST BE EQUI-SPACED AT INTERVALS OF DT. FURTHER, THE TRANS- */
/* FORM VALUES, CONSIDERED TO BE IN FREQUENCY SPACE, WILL BE EQUI- */
/* SPACED FROM 0 TO 2*PI*(NN(I)-1)/(NN(I)*DT) AT INTERVALS OF */
/* 2*PI/(NN(I)*DT) FOR EACH DIMENSION OF LENGTH NN(I). OF COURSE, */
/* DT NEED NOT BE THE SAME FOR EVERY DIMENSION. */
/* THE CALLING SEQUENCE IS-- */
/* CALL FOURT(DATAR,DATAI,NN,NDIM,IFRWD,ICPLX,WORKR,WORKI) */
/* DATAR AND DATAI ARE THE ARRAYS USED TO HOLD THE REAL AND IMAGINARY */
/* PARTS OF THE INPUT DATA ON INPUT AND THE TRANSFORM VALUES ON */
/* OUTPUT. THEY ARE FLOATING POINT ARRAYS, MULTIDIMENSIONAL WITH */
/* IDENTICAL DIMENSIONALITY AND EXTENT. THE EXTENT OF EACH DIMENSION */
/* IS GIVEN IN THE INTEGER ARRAY NN, OF LENGTH NDIM. THAT IS, */
/* NDIM IS THE DIMENSIONALITY OF THE ARRAYS DATAR AND DATAI. */
/* IFRWD IS AN INTEGER USED TO INDICATE THE DIRECTION OF THE FOURIER */
/* TRANSFORM. IT IS NON-ZERO TO INDICATE A FORWARD TRANSFORM */
/* (EXPONENTIAL SIGN IS -) AND ZERO TO INDICATE AN INVERSE TRANSFORM */
/* (SIGN IS +). ICPLX IS AN INTEGER TO INDICATE WHETHER THE DATA */
/* ARE REAL OR COMPLEX. IT IS NON-ZERO FOR COMPLEX, ZERO FOR REAL. */
/* IF IT IS ZERO (REAL) THE CONTENTS OF ARRAY DATAI WILL BE ASSUMED */
/* TO BE ZERO, AND NEED NOT BE EXPLICITLY SET TO ZERO. AS EXPLAINED */
/* ABOVE, THE TRANSFORM RESULTS ARE ALWAYS COMPLEX AND ARE STORED */
/* IN DATAR AND DATAI ON RETURN. WORKR AND WORKI ARE ARRAYS USED */
/* FOR WORKING STORAGE. THEY ARE NOT NECESSARY IF ALL THE DIMENSIONS */
/* OF THE DATA ARE POWERS OF TWO. IN THIS CASE, THE ARRAYS MAY BE */
/* REPLACED BY THE NUMBER 0 IN THE CALLING SEQUENCE. THUS, USE OF */
/* POWERS OF TWO CAN FREE A GOOD DEAL OF STORAGE. IF ANY DIMENSION */
/* IS NOT A POWER OF TWO, THESE ARRAYS MUST BE SUPPLIED. THEY ARE */
/* FLOATING POINT, ONE DIMENSIONAL OF LENGTH EQUAL TO THE LARGEST */
/* ARRAY DIMENSION, THAT IS, TO THE LARGEST VALUE OF NN(I). */
/* WORKR AND WORKI, IF SUPPLIED, MUST NOT BE THE SAME ARRAYS AS DATAR */
/* OR DATAI. ALL SUBSCRIPTS OF ALL ARRAYS BEGIN AT 1. */
/* THERE ARE NO ERROR MESSAGES OR ERROR HALTS IN THIS PROGRAM. THE */
/* PROGRAM RETURNS IMMEDIATELY IF NDIM OR ANY NN(I) IS LESS THAN ONE. */
/* 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)
{
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];
}
/*main loop for each dimension*/
np1 = 1;
for (idim = 1; idim <= ndim; idim++) {
n = nn[idim-1];
np2 = np1*n;
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;
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;
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;
L60:
ifact[iff] = m;
L70:
non2p = np2/ntwo;
/*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*/
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*/
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*/
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;
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;
L540:
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;
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*/
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;
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;
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];
}
L740:
np2 *= 2;
ntot *= 2;
j = ntot+1;
imax = ntot/2+1;
L745:
imin = imax-nhalf;
i = imin;
goto L755;
L750:
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;
}
L765:
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;
L780:
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*/
L900:
np0 = np1;
np1 = np2;
nprev = n;
}
L920:
return;
}

@ -1,16 +0,0 @@
#include <stdio.h>
#include <math.h>
#include "genlib.h"
/*gamma covariance function*/
double gammf(double h, double alpha)
{
float delta;
double z;
delta = pow(20.,1./alpha)-1.;
z = 1./(double)(pow(1.+h*delta,alpha));
return(z);
}

@ -1,31 +0,0 @@
#include <math.h>
#include "genlib.h"
#define NTAB 32
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;
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);
fac = sqrt(-2.0*log(rsq)/rsq);
gset = v1*fac;
iset = 1;
return (v2*fac);
} else {
iset = 0;
return (gset);
}
}

@ -1,10 +0,0 @@
#include <stdio.h>
#include <math.h>
#include "genlib.h"
/*gaussian covariance function*/
double gaussian(double h)
{
return (exp(-3.*(double)(h*h)));
}

@ -1,43 +0,0 @@
#include <stdlib.h>
#include <math.h>
#include "geostat.h"
/* GENERATION OF A GAUSSIAN WHITE NOISE VECTOR */
/*input: */
/* seed: seed */
/* n: number of components in the vector */
/*output: */
/* realization: structure defining the 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;
}

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

@ -1,39 +0,0 @@
#include <math.h>
/* compute the length for one dimension*/
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;
} 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);
}
}
return (n);
}

@ -1,36 +0,0 @@
#
########################################################################
#
# Makefile for library
#
########################################################################
#
INCLUDE = -I./
LIBS = -lm -L./
# CCFLAGS = -fast
ifeq (${ARCH},sun)
CC = cc
CCFLAGS = -g
endif
ifeq (${ARCH},linux)
CC = gcc
CCFLAGS = -g
endif
LIB = libFFTMA_debug_${ARCH}.a
NOBJS= gammf.o FFTMA.o addstat.o axes.o cgrid.o covariance.o fourt.o length.o maxfactor.o test_fact.o cov_value.o generate.o gasdev.o ran2.o stable.o gaussian.o power.o cubic.o spherical.o nugget.o exponential.o cardsin.o nor2log.o
.c.o:
${CC} $(INCLUDE) $(CCFLAGS) -c $<
# LIBRARY
$(LIB) : $(NOBJS)
ar cr $(LIB) $(NOBJS)
ranlib $(LIB)
clean :
rm *.o

@ -1,35 +0,0 @@
#include "genlib.h"
/*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;
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 != 1 ) {
if (lnum > maxfac)
maxfac = lnum;
}
}
return (maxfac);
}

@ -1,78 +0,0 @@
#include <stdlib.h>
#include <math.h>
#include "geostat.h"
/*TURNS NORMAL NUMBERS INTO LOGNORMAL NUMBERS */
/*input: */
/*realin: structure defining a realization - */
/* normal numbers */
/*typelog: --> 3: lognormal (natural) */
/* --> 4: lognormal (log10) */
/*output: */
/*realout: structure defining a realization - */
/* lognormal numbers */
void nor2log(struct realization_mod *realin, int typelog, struct realization_mod *realout)
{
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;
}
}
(*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;
}
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;
}
}
return;
}

@ -1,14 +0,0 @@
#include <stdio.h>
#include <math.h>
#include "genlib.h"
/*nugget covariance function*/
double nugget(double h)
{
if (h == 0) {
return (1.);
} else {
return(0.);
}
}

@ -1,10 +0,0 @@
#include <stdio.h>
#include <math.h>
#include "genlib.h"
/*power covariance function*/
double power(double h, double alpha)
{
return(pow(h,alpha));
}

@ -1,50 +0,0 @@
#include "genlib.h"
#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
double ran2(long *idum, long *idum2, long *iy, long iv[NTAB])
{
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;
}
*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);
}

@ -1,17 +0,0 @@
#include <stdio.h>
#include <math.h>
#include "genlib.h"
/*spherical covariance function*/
double spherical(double h)
{
double z;
if (h >= 1.) {
z = 0.;
} else {
z = 1.-1.5*(double)h+0.5*(double)(h*h*h);
}
return (z);
}

@ -1,10 +0,0 @@
#include <stdio.h>
#include <math.h>
#include "genlib.h"
/*stable covariance function*/
double stable(double h, double alpha)
{
return(exp(-3.*(double)pow(h,alpha)));
}

@ -1,26 +0,0 @@
#include "genlib.h"
/*tries factor*/
int test_fact(int *pnum, int fact, int *pmaxfac)
{
int power, t;
power = 0;
while ( ( t = *pnum / fact ) * fact == *pnum ) {
++power;
*pnum = t;
}
if ( power != 0 ) {
if (fact > *pmaxfac)
*pmaxfac = fact;
}
if ( t > fact ) {
return (1);
}
return (0);
}

@ -1,110 +0,0 @@
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "geostat.h"
/*FAST FOURIER TRANSFORM MOVING AVERAGE METHOD */
/*Turns a Gaussian white noise vector into a */
/*spatially correlated vector */
/*input: */
/*variogram: structure defining the variogram */
/* model */
/*grid: structure defining the grid */
/*n: vector with the number of cells along the */
/* X, Y and Z axes for the underlying grid */
/* i = [0 1 2] */
/* --> 0 0 0 : n will be computed and */
/* updated as output */
/* --> nx ny nz: these dimensions are used */
/*realin: structure defining a realization - */
/* must be a Gaussian white noise */
/*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)
{
int NTOT,i,j,k,NMAX,NDIM,ntot,nmax,NXYZ,nxyz,maille0,maille1;
int solver;
double temp;
double *ireal,*covar,*workr,*worki,*realization;
string variable;
/*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); */
allouememoire(covar,'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;
}

@ -1,30 +0,0 @@
########################################################################
#
# Makefile for library
#
########################################################################
INTERFACE = ../include
INCLUDE = -I${INTERFACE}
LIBS = -lm -L../lib
CC = cc
CCFLAGS =
LIB = libFFTMA2_${ARCH}.a
NOBJS= kgeneration.o kgeneration2.o fftma2.o prebuild_gwn.o build_real.o addstat2.o clean_real.o
.c.o:
${CC} $(INCLUDE) $(CCFLAGS) -c $<
# LIBRARY
$(LIB) : $(NOBJS)
ar cr $(LIB) $(NOBJS)
ranlib $(LIB)
clean :
rm *.o

@ -1,119 +0,0 @@
#include <stdio.h>
#include <math.h>
#include "genlib.h"
#include "geostat.h"
/*addstat */
/*adds mean and variance effects to the N(0,1) realization*/
/*input: */
/*realin: structure defining a realization - */
/* must have zeio mean and unit variance */
/*stat: structure defining the mean and variance */
/*output: */
/*realout: structure defining a realization - */
void addstat2(struct realization_mod *realin,struct statistic_mod stat ,struct realization_mod *realout,struct realization_mod *realout2)
{
int i,nblockm,nblockv;
double r,moy,var;
/*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)
Error("No memory available");
}
(*realout).n = (*realin).n;
if ((*realout2).vector == NULL || (*realout2).n != (*realin).n) {
(*realout2).vector = (double *) malloc((*realin).n * sizeof(double));
if ((*realout2).vector == NULL)
Error("No memory available");
}
(*realout2).n = (*realin).n;
/*test over the input realization*/
switch ((*realin).code) {
case 0:
case 1:
(*realout).code = 2;
(*realout2).code = 2;
break;
case 6:
(*realout).code = 7;
(*realout2).code = 7;
break;
default:
(*realout).code = (*realin).code;
(*realout2).code = (*realin).code;
break;
}
for (i = 0; i < (*realin).n; i++) {
/*mean*/
switch (stat.nblock_mean) {
case 1:
/*stationary case*/
nblockm = 1;
break;
default:
/*number of the cell - nonstationary case*/
nblockm = i+1;
break;
}
/*variance*/
switch (stat.nblock_var) {
case 1:
/*stationary case*/
nblockv = 1;
break;
default:
/*number of the cell - nonstationary case*/
nblockv = i+1;
break;
}
/* switch (stat.type) { */
/* case 0: */
/*normal case*/
moy = stat.mean[nblockm-1];
var = stat.variance[nblockv-1];
/* break; */
/* case 1: */
/*lognormal (natural) case*/
/* r = (double)sqrt(stat.variance[nblockv-1])/stat.mean[nblockm-1]; */
/* r *= r; */
/* moy = (double)log(stat.mean[nblockm-1]/sqrt(1.0+r)); */
/* var = (double)log(1.0+r); */
/* break; */
/* case 2: */
/*lognormal (log10) case*/
/* r = (double)sqrt(stat.variance[nblockv-1])/stat.mean[nblockm-1]; */
/* r *= r; */
/* moy = (double)log10(stat.mean[nblockm-1]/sqrt(1.0+r)); */
/* var = (double)log10(1.0+r); */
/* break; */
/* default: */
/* Error("Type not defined in addstat"); */
/* break; */
/* } */
(*realout).vector[i] = (double)sqrt(var)*(*realin).vector[i];
(*realout2).vector[i] = (double)sqrt(var)*(*realin).vector[i]+moy;
}
return;
}

@ -1,45 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include <math.h>
#include "geostat.h"
/* build_real */
/* build a realization in the spectral domain */
/*n: vector with the number of cells along the */
/* X, Y and Z axes for the underlying grid */
/* i = [0 1 2] */
/* --> 0 0 0 : n will be computed and */
/* updated as output */
/* --> nx ny nz: these dimensions are used */
/*covar: vector defining the covariance in spectral domain */
/*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)
{
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;
}
}
}
return;
}

@ -1,43 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include <math.h>
#include "geostat.h"
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*/
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];
}
}
}
return;
}

@ -1,105 +0,0 @@
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "geostat.h"
/*FAST FOURIER TRANSFORM MOVING AVERAGE METHOD */
/*Turns a Gaussian white noise vector into a */
/*spatially correlated vector */
/*input: */
/*variogram: structure defining the variogram */
/* model */
/*grid: structure defining the grid */
/*n: vector with the number of cells along the */
/* X, Y and Z axes for the underlying grid */
/* i = [0 1 2] */
/* --> 0 0 0 : n will be computed and */
/* updated as output */
/* --> nx ny nz: these dimensions are used */
/*realin: structure defining a realization - */
/* must be a Gaussian white noise */
/*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)
{
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;
}

@ -1,79 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include "genlib.h"
#include "simpio.h"
#include "geostat.h"
#include "toolsIO.h"
#include "toolsFFTMA.h"
/* 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 kgeneration(long seed,struct grid_mod grid,struct statistic_mod stat,struct vario_mod variogram,string filename[7],struct realization_mod *Z,struct realization_mod *Y,struct realization_mod *Y1, int n[3], int *genere, int *gwnwrite, struct realization_mod *gwnoise)
{
int i,N;
int typelog;
string file1,file2;
/*generate Gaussian white noise*/
N = grid.NX*grid.NY*grid.NZ;
n[0] = 0;
n[1] = 0;
n[2] = 0;
printf("\n\n\n");
switch (*genere)
{
case 0:
generate(&seed,N,Z);
/*save the Gaussian white noise file*/
if (*gwnwrite == 0)
{
writefile(filename[0],Z);
}
break;
case 1:
(*Z).vector=(double *) malloc(N*sizeof(double));
(*Z).n = N;
for (i=0;i<N;i++)
{
(*Z).vector[i]=(*gwnoise).vector[i];
}
break;
}
/*FFTMA*/
FFTMA2(variogram,grid,n,Z,Y);
/* file1="prout"; */
/* writefile(file1,Y); */
/*add the statistics*/
if (stat.mean[0] != 0 || stat.variance[0]!= 1)
addstat2(Y,stat,Y1,Y);
/* file2="prout2"; */
/* writefile(file2,Y); */
/* make a log normal realization */
if (stat.type==1 || stat.type==2){
typelog=stat.type+2;
nor2log(Y1,typelog,Y1);
nor2log(Y,typelog,Y);
}
/*save the realization*/
/* writefile(filename[0],Y1); */
/*save the permeability realization*/
writefile_bin(filename[1],Y);
readfile_bin(filename[1],Y);
return;
}

@ -1,76 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include "genlib.h"
#include "simpio.h"
#include "geostat.h"
#include "toolsIO.h"
#include "toolsFFTMA.h"
/* 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 kgeneration2(long seed,struct grid_mod grid,struct statistic_mod stat,struct vario_mod variogram,string filename[7],struct realization_mod *Z,struct realization_mod *Y,struct realization_mod *Y1, int n[3], int *genere, int *gwnwrite, struct realization_mod *gwnoise,int *format_file)
{
int i,N;
int typelog;
string file1,file2;
/*generate Gaussian white noise*/
N = grid.NX*grid.NY*grid.NZ;
n[0] = 0;
n[1] = 0;
n[2] = 0;
printf("\n\n\n");
switch (*genere)
{
case 0:
generate(&seed,N,Z);
/*save the Gaussian white noise file*/
if (*gwnwrite == 0)
{
writefile(filename[0],Z);
}
break;
case 1:
(*Z).vector=(double *) malloc(N*sizeof(double));
(*Z).n = N;
for (i=0;i<N;i++)
{
(*Z).vector[i]=(*gwnoise).vector[i];
}
break;
}
/*FFTMA*/
FFTMA2(variogram,grid,n,Z,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){
typelog=stat.type+2;
/* nor2log(Y1,typelog,Y1); */
nor2log(Y,typelog,Y);
}
switch (*format_file)
{
case 0:
writefile(filename[1],Y);
break;
case 1:
writefile_bin(filename[1],Y);
break;
}
return;
}

@ -1,54 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include <math.h>
#include "geostat.h"
/* prebuild_gwn */
/* Produce a first construction in real space of the Gaussian white noise */
/* grid: structure defining the grid */
/*n: vector with the number of cells along the */
/* X, Y and Z axes for the underlying grid */
/* i = [0 1 2] */
/* --> 0 0 0 : n will be computed and */
/* updated as output */
/* --> nx ny nz: these dimensions are used */
/*realin: structure defining a realization - */
/* 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)
{
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<ntot;i++)
{
realization[i+1]=(*realin).vector[i];
}
}
else
{
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];
if (i <= grid.NX && j <= grid.NY && k <= grid.NZ) {
maille0 = i-1+(j-1+(k-1)*grid.NY)*grid.NX;
realization[maille1] = (*realin).vector[maille0];
} else {
realization[maille1] = 0.;
}
}
}
}
}
return;
}

@ -1,82 +0,0 @@
#include <stdlib.h>
#include <string.h>
#include "genlib.h"
#include "geostat.h"
#ifndef _TOOLSFFTMA_H
#define _TOOLSFFTMA_H
/* Create december, the 16th 2002 */
/* Modified december, the 9th 2002 */
/* List of functions: */
/* ------------------ */
/* kgeneration, FFTMA2, prebuild_gwn, build_real, addstat2, clean_real */
/*FUNCTIONS*/
/*----------*/
void kgeneration(long seed,struct grid_mod grid,struct statistic_mod stat,struct vario_mod variogram,string filename[8],struct realization_mod *Z,struct realization_mod *Y,struct realization_mod *Y2, int n[3], int *genere, int *gwnwrite, struct realization_mod *gwnoise);
void kgeneration2(long seed,struct grid_mod grid,struct statistic_mod stat,struct vario_mod variogram,string filename[7],struct realization_mod *Z,struct realization_mod *Y,struct realization_mod *Y1, int n[3], int *genere, int *gwnwrite, struct realization_mod *gwnoise,int *format_file) ;
/*FAST FOURIER TRANSFORM Pressure Simulation */
/*Turns a Gaussian white noise vector into a */
/*spatially correlated vector */
/*input: */
/*variogram: structure defining the variogram */
/* model */
/*grid: structure defining the grid */
/*n: vector with the number of cells along the */
/* X, Y and Z axes for the underlying grid */
/* i = [0 1 2] */
/* --> 0 0 0 : n will be computed and */
/* updated as output */
/* --> nx ny nz: these dimensions are used */
/*realin: structure defining a realization - */
/* must be a Gaussian white noise */
/*gradient: macroscopic gradient pression vector */
/*output: */
/*realout: structure defining a realization - */
/*realout2: structure defining a pressure field */
/*realout3: structure defining a xvelocity field */
/*realout4: structure defining a yvelocity field */
/*realout5: structure defining a zvelocity field */
void FFTMA2(struct vario_mod variogram,struct grid_mod grid,int n[3],struct realization_mod *realin,struct realization_mod *realout);
/* prebuild_gwn */
/* Produce a first construction in real space of the Gaussian white noise */
/* grid: structure defining the grid */
/*n: vector with the number of cells along the */
/* X, Y and Z axes for the underlying grid */
/* i = [0 1 2] */
/* --> 0 0 0 : n will be computed and */
/* updated as output */
/* --> nx ny nz: these dimensions are used */
/*realin: structure defining a realization - */
/* 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);
/* build_real */
/* build a realization in the spectral domain */
/*n: vector with the number of cells along the */
/* X, Y and Z axes for the underlying grid */
/* i = [0 1 2] */
/* --> 0 0 0 : n will be computed and */
/* updated as output */
/* --> nx ny nz: these dimensions are used */
/*covar: vector defining the covariance in spectral domain */
/*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 addstat2(struct realization_mod *realin,struct statistic_mod stat ,struct realization_mod *realout,struct realization_mod *realout2);
void clean_real(struct realization_mod *realin,int n[3],struct grid_mod grid,double *vectorresult,struct realization_mod *realout);
#endif // define _TOOLSFFTMA_H

@ -1,162 +0,0 @@
#include <stdlib.h>
#include <math.h>
#include "geostat.h"
#include "toolsludo.h"
/*FAST FOURIER TRANSFORM Pressure Simulation */
/*Turns a Gaussian white noise vector into a */
/*spatially correlated vector */
/*input: */
/*variogram: structure defining the variogram */
/* model */
/*grid: structure defining the grid */
/*n: vector with the number of cells along the */
/* X, Y and Z axes for the underlying grid */
/* i = [0 1 2] */
/* --> 0 0 0 : n will be computed and */
/* updated as output */
/* --> nx ny nz: these dimensions are used */
/*realin: structure defining a realization - */
/* must be a Gaussian white noise */
/*gradient: macroscopic gradient pression vector */
/*output: */
/*realout: structure defining a realization - */
/*realout2: structure defining a pressure field */
/*realout3: structure defining a xvelocity field */
/*realout4: structure defining a yvelocity field */
/*realout5: structure defining a zvelocity field */
void FFTPSim(struct vario_mod variogram,struct statistic_mod stat,struct grid_mod grid,int n[3],struct realization_mod *realin,struct pressure_mod gradient,struct realization_mod *realout,struct realization_mod *realout2,struct realization_mod *realout3,struct realization_mod *realout4,struct realization_mod *realout5)
{
int NTOT,i,j,k,NMAX,NDIM,ntot,nmax,NXYZ,nxyz,maille0,maille1;
double *workr,*worki,temp,temp2,coeff;
double *covar,*realization,*pressure;
double *icovar,*ireal,*ipressure;
double *xvelocity,*ixvelocity,*yvelocity,*iyvelocity,*zvelocity,*izvelocity;
double ki,kj,kk;
/*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);
icovar = (double *) malloc(ntot * sizeof(double));
testmemory(icovar);
realization = (double *) malloc(ntot * sizeof(double));
testmemory(realization);
ireal = (double *) malloc(ntot * sizeof(double));
testmemory(ireal);
pressure = (double *) malloc(ntot * sizeof(double));
testmemory(pressure);
ipressure = (double *) malloc(ntot * sizeof(double));
testmemory(ipressure);
xvelocity = (double *) malloc(ntot * sizeof(double));
testmemory(xvelocity);
ixvelocity = (double *) malloc(ntot * sizeof(double));
testmemory(ixvelocity);
yvelocity = (double *) malloc(ntot * sizeof(double));
testmemory(yvelocity);
iyvelocity = (double *) malloc(ntot * sizeof(double));
testmemory(iyvelocity);
zvelocity = (double *) malloc(ntot * sizeof(double));
testmemory(zvelocity);
izvelocity = (double *) malloc(ntot * sizeof(double));
testmemory(izvelocity);
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,icovar,n,NDIM,1,0,workr,worki);
/*organization of the input Gaussian white noise*/
prebuild_gwn(grid,n,realin,realization);
/*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);
/* pressure calculation in the spectral domain*/
build_pressure(n,grid,gradient,realization,ireal,pressure,ipressure);
build_velocity(n,grid,stat,gradient,realization,ireal,xvelocity,ixvelocity,1);
build_velocity(n,grid,stat,gradient,realization,ireal,yvelocity,iyvelocity,2);
build_velocity(n,grid,stat,gradient,realization,ireal,zvelocity,izvelocity,3);
/*backward fourier transform*/
fourt(realization,ireal,n,NDIM,0,1,workr,worki);
fourt(pressure,ipressure,n,NDIM,0,1,workr,worki);
fourt(xvelocity,ixvelocity,n,NDIM,0,1,workr,worki);
fourt(yvelocity,iyvelocity,n,NDIM,0,1,workr,worki);
fourt(zvelocity,izvelocity,n,NDIM,0,1,workr,worki);
free(ireal);
free(ipressure);
free(ixvelocity);
free(iyvelocity);
free(izvelocity);
free(workr);
free(worki);
/*output realization*/
/*is the output realization already allocated?*/
/*if not, memory allocation*/
clean_real(realin,n,grid,realization,realout);
clean_real(realin,n,grid,pressure,realout2);
clean_real(realin,n,grid,xvelocity,realout3);
clean_real(realin,n,grid,yvelocity,realout4);
clean_real(realin,n,grid,zvelocity,realout5);
free(realization);
free(pressure);
free(xvelocity);
free(yvelocity);
free(zvelocity);
return;
}

@ -1,157 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include <math.h>
#include "geostat.h"
#include "pressure.h"
/*FAST FOURIER TRANSFORM Pressure Simulation */
/*Turns a Gaussian white noise vector into a */
/*spatially correlated vector */
/*input: */
/*variogram: structure defining the variogram */
/* model */
/*grid: structure defining the grid */
/*n: vector with the number of cells along the */
/* X, Y and Z axes for the underlying grid */
/* i = [0 1 2] */
/* --> 0 0 0 : n will be computed and */
/* updated as output */
/* --> nx ny nz: these dimensions are used */
/*realin: structure defining a realization - */
/* must be a Gaussian white noise */
/*gradient: macroscopic gradient pression vector */
/*output: */
/*realout: structure defining a realization - */
/*realout2: structure defining a pressure field */
/*realout3: structure defining a xvelocity field */
/*realout4: structure defining a yvelocity field */
/*realout5: structure defining a zvelocity field */
void FFTPressure(int n[3],struct grid_mod grid,struct realization_mod *realin,struct statistic_mod stat,struct pressure_mod gradient,struct realization_mod *realout,struct realization_mod *realout2,struct realization_mod *realout3,struct realization_mod *realout4,int solver)
{
int NTOT,i,j,k,NMAX,NDIM,ntot,nmax,NXYZ,nxyz,maille0,maille1;
double *workr,*worki,temp,temp2,coeff;
double *realization,*pressure;
double *ireal,*ipressure;
double *xvelocity,*ixvelocity,*yvelocity,*iyvelocity,*zvelocity,*izvelocity;
double ki,kj,kk;
FILE *fp;
/* string nomfichier; */
/*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;
realization = (double *) malloc(ntot * sizeof(double));
testmemory(realization);
ireal = (double *) malloc(ntot * sizeof(double));
testmemory(ireal);
pressure = (double *) malloc(ntot * sizeof(double));
testmemory(pressure);
ipressure = (double *) malloc(ntot * sizeof(double));
testmemory(ipressure);
xvelocity = (double *) malloc(ntot * sizeof(double));
testmemory(xvelocity);
ixvelocity = (double *) malloc(ntot * sizeof(double));
testmemory(ixvelocity);
yvelocity = (double *) malloc(ntot * sizeof(double));
testmemory(yvelocity);
iyvelocity = (double *) malloc(ntot * sizeof(double));
testmemory(iyvelocity);
zvelocity = (double *) malloc(ntot * sizeof(double));
testmemory(zvelocity);
izvelocity = (double *) malloc(ntot * sizeof(double));
testmemory(izvelocity);
workr = (double *) malloc(nmax * sizeof(double));
testmemory(workr);
worki = (double *) malloc(nmax * sizeof(double));
testmemory(worki);
/*organization of the realization*/
prebuild_gwn(grid,n,realin,realization,solver);
/*forward fourier transform of the GWN*/
fourt(realization,ireal,n,NDIM,1,0,workr,worki);
/* pressure calculation in the spectral domain*/
build_pressure(n,grid,gradient,realization,ireal,pressure,ipressure);
build_velocity(n,grid,stat,gradient,realization,ireal,xvelocity,ixvelocity,1);
build_velocity(n,grid,stat,gradient,realization,ireal,yvelocity,iyvelocity,2);
build_velocity(n,grid,stat,gradient,realization,ireal,zvelocity,izvelocity,3);
/*backward fourier transform*/
fourt(realization,ireal,n,NDIM,0,1,workr,worki);
fourt(pressure,ipressure,n,NDIM,0,1,workr,worki);
fourt(xvelocity,ixvelocity,n,NDIM,0,1,workr,worki);
fourt(yvelocity,iyvelocity,n,NDIM,0,1,workr,worki);
fourt(zvelocity,izvelocity,n,NDIM,0,1,workr,worki);
free(ireal);
free(ipressure);
free(ixvelocity);
free(iyvelocity);
free(izvelocity);
free(workr);
free(worki);
for (i=1;i<=NTOT;i++)
realization[i]=realization[i]/(double) NTOT;
fp = fopen("realization.test", "w");
for (i=1;i<=NTOT;i++)
fprintf(fp,"%f\n",realization[i]);
fclose(fp);
/* nomfichier="pression.test"; */
fp = fopen("pression.test", "w");
for (i=1;i<=NTOT;i++)
fprintf(fp,"%f\n",pressure[i]);
fclose(fp);
/*output realization*/
/*is the output realization already allocated?*/
/*if not, memory allocation*/
clean_real2(realin,n,grid,solver,pressure,realout);
clean_real2(realin,n,grid,solver,xvelocity,realout2);
clean_real2(realin,n,grid,solver,yvelocity,realout3);
clean_real2(realin,n,grid,solver,zvelocity,realout4);
/* free(realization); */
/* free(pressure); */
/* free(xvelocity); */
/* free(yvelocity); */
/* free(zvelocity); */
return;
}

@ -1,68 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include <math.h>
#include "geostat.h"
/*FAST FOURIER TRANSFORM Test */
void FFTtest(int n[3],struct grid_mod grid,struct realization_mod *realin,struct realization_mod *realout,int solver)
{
int NTOT,i,j,k,NMAX,NDIM,ntot,nmax,NXYZ,nxyz,maille0,maille1;
double *workr,*worki,temp,temp2,coeff;
double *realization;
double *ireal;
double ki,kj,kk;
FILE *fp;
/*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;
realization = (double *) malloc(ntot * sizeof(double));
testmemory(realization);
ireal = (double *) malloc(ntot * sizeof(double));
testmemory(ireal);
workr = (double *) malloc(nmax * sizeof(double));
testmemory(workr);
worki = (double *) malloc(nmax * sizeof(double));
testmemory(worki);
/*organization of the realization*/
prebuild_gwn(grid,n,realin,realization,solver);
fp = fopen("perm.test", "w");
for (i=1;i<=NTOT;i++)
fprintf(fp,"%f\n",realization[i]);
fclose(fp);
/*forward fourier transform of the GWN*/
fourt(realization,ireal,n,NDIM,1,0,workr,worki);
/*backward fourier transform*/
fourt(realization,ireal,n,NDIM,0,1,workr,worki);
fp = fopen("perm2.test", "w");
for (i=1;i<=NTOT;i++)
fprintf(fp,"%f\n",realization[i]);
fclose(fp);
return;
}

@ -1,34 +0,0 @@
########################################################################
#
# Makefile for library
#
########################################################################
#INTERFACE = ${BOSSE}/CODES/LIBS_C/Interface
#INTERFACE2= ${BOSSE}/CODES/Geostat/FFTPSim/Version_binaire2/src/Interface
INTERFACE = ../include
#INTERFACE2= ../Interface
INCLUDE = -I${INTERFACE}
LIBS = -lm -L../lib
CCFLAGS =
CC = cc
LIB = libFFTPSIM_${ARCH}.a
NOBJS= pgeneration.o pgeneration2.o FFTPressure.o FFTtest.o build_pressure.o build_velocity.o total_pressure.o total_velocity.o clean_real2.o waveVectorCompute3D.o mat_vec.o derivReal.o
.c.o:
${CC} $(INCLUDE) $(CCFLAGS) -c $<
# LIBRARY
$(LIB) : $(NOBJS)
ar cr $(LIB) $(NOBJS)
ranlib $(LIB)
clean :
rm *.o

@ -1,66 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include <math.h>
#include "geostat.h"
#include "pressure.h"
/* Build_pressure */
/* build pressure in spectral domain */
void build_pressure(int n[3],struct grid_mod grid,struct pressure_mod gradient,double *realization,double *ireal,double *pressure,double *ipressure)
{
int i,j,k,maille1;
double ki,kj,kk;
double coeff,temp,temp2;
/* pressure calculation 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];
if (n[0]==1)
{
ki=0.;
}
else
{
ki =(double)(i-1)/(double)(grid.NX*grid.DX);
}
if (n[1]==1)
{
kj=0.;
}
else
{
kj =(double)j/(double)(grid.NY*grid.DY);
}
if (n[2]==1)
{
kk=0.;
}
else
{
kk =(double)k/(double)(grid.NZ*grid.DZ);
}
coeff = (gradient.x*ki+gradient.y*kj+gradient.z*kk)/(ki*ki+kj*kj+kk*kk);
temp = realization[maille1];
temp2= ireal[maille1];
if (maille1==1)
{
pressure[maille1] =0.;
ipressure[maille1]=0.;
}
else
{
pressure[maille1] =-1./(2*3.14)*coeff*temp2;
ipressure[maille1]= 1./(2*3.14)*coeff*temp;
}
}
}
}
return;
}

@ -1,117 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <stdarg.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "geostat.h"
#include "pressure.h"
/* Build_velocity */
/* Build velocity in spectral domain */
void build_velocity(int n[3],struct grid_mod grid,struct statistic_mod stat,struct pressure_mod gradient,double *realization,double *ireal,double *xvelocity,double *ixvelocity,int direction)
{
int i,j,k,maille1;
int x,y,z;
double ki,kj,kk;
double coeff,temp,temp2;
double grad;
/* x-direction velocity calculation in the spectral domain*/
switch(direction)
{
case 1:
grad = gradient.x;
x=1;
y=0;
z=0;
break;
case 2:
grad = gradient.y;
x=0;
y=1;
z=0;
break;
case 3:
grad = gradient.z;
x=0;
y=0;
z=1;
break;
default :
printf("build_velocity.c: wrong velocity direction!!! direction: %d\n",direction);
break;
}
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];
/* if (i<n[0]/2) { */
/* ki =(double)i/(double)(grid.NX*grid.DX); */
/* } else { */
/* ki =-(double)(n[0]-i+1)/(double)(grid.NX*grid.DX); */
/* } */
/* if (j<n[1]/2) { */
/* kj =(double)j/(double)(grid.NY*grid.DY); */
/* } else { */
/* kj =-(double)(n[1]-j+1)/(double)(grid.NY*grid.DY); */
/* } */
/* if (k<n[2]/2) { */
/* kk =(double)k/(double)(grid.NZ*grid.DZ); */
/* } else { */
/* kk =-(double)(n[2]-k+1)/(double)(grid.NZ*grid.DZ); */
/* } */
if (n[0]==1)
{
ki=0;
}
else if (i<n[0]/2)
{
ki =(double)i/(double)(grid.NX*grid.DX);
}
else
{
ki =-(double)(n[0]-i+1)/(double)(grid.NX*grid.DX);
}
if (n[1]==1)
{
kj=0;
}
else if (j<n[1]/2)
{
kj =(double)j/(double)(grid.NY*grid.DY);
}
else
{
kj =-(double)(n[1]-j+1)/(double)(grid.NY*grid.DY);
}
if (n[2]==1)
{
kk=0;
}
else if (k<n[2]/2)
{
kk =(double)k/(double)(grid.NZ*grid.DZ);
}
else
{
kk =-(double)(n[2]-k+1)/(double)(grid.NZ*grid.DZ);
}
coeff = (gradient.x*ki+gradient.y*kj+gradient.z*kk)/(ki*ki+kj*kj+kk*kk);
temp = realization[maille1];
temp2= ireal[maille1];
xvelocity[maille1] = stat.mean[0]*(grad-(ki*x+kj*y+kk*z)*coeff)*temp;
ixvelocity[maille1]= stat.mean[0]*(grad-(ki*x+kj*y+kk*z)*coeff)*temp2;
}
}
}
return;
}

@ -1,50 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include <math.h>
#include "geostat.h"
void clean_real2(struct realization_mod *realin,int n[3],struct grid_mod grid,int solver,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*/
if ((*realout).vector == NULL || (*realout).n != (*realin).n)
{
(*realout).vector = (double *) malloc((*realin).n * sizeof(double));
if ((*realout).vector == NULL)
{
printf("Clean_real2.c: No memory available\n");
exit;
}
}
(*realout).n = (*realin).n;
(*realout).code = 1;
if (solver==1)
{
for ( k = 1; k <= NTOT;k++)
{
(*realout).vector[k-1] = vectorresult[k]/(double) NTOT;
}
}
else
{
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;
(*realout).vector[maille0] = vectorresult[maille1]/(double) NTOT;
}
}
}
}
return;
}

@ -1,216 +0,0 @@
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "condor.h"
#include "geostat.h"
/* Private functions */
void normAxes(double *vec, double *normed);
void derivReal(struct realization_mod *Z, struct realization_mod *dZ, double *dir, struct grid_mod *grid, struct vario_mod vario) {
int n[3],i,j,k,maille,extMaille;
int NTOT,ntot,NMAX,nmax,NDIM;
// int NXYZ,nxyz;
double nDir[3];
double *table,*workr,*worki,*realization,*waveVect;
int Ngrid;
double tmp;
double *ExtendedWaveVect;
/* Test the input real*/
/* if ((*Z).code != 1) { */
/* printf("Realization should be Standard Normal\n"); */
/* return; */
/* } */
printf("grid.Nx = %d\n",(*grid).NX);
printf("grid.Ny = %d\n",(*grid).NY);
printf("grid.Nz = %d\n",(*grid).NZ);
printf("vario.Nvario = %d\n",vario.Nvario);
for(i=0; i<vario.Nvario;i++) {
printf("Component %d:\n",i);
printf("Type: %d\n",vario.vario[i]);
printf("Exponent: %f\n",vario.alpha[i]);
printf("Anisotropy axes:");
for(j=0; j<6; j++) {
printf(" %f",vario.ap[6*i+j]);
}
printf("\nCorrelation length:\n");
for(j=0;j<3;j++) {
printf("axe%d :%f\n",j+1,vario.scf[3*i+j]);
}
printf("Variance: %f\n",vario.var[i]);
}
/* covariance axis normalisation */
axes(vario.ap,vario.scf,vario.Nvario);
n[0]=0;
n[1]=0;
n[2]=0;
/* pseudo-grid definition */
cgrid(vario,*grid,n);
printf("n[0] = %d\n", n[0]);
printf("n[1] = %d\n", n[1]);
printf("n[2] = %d\n", n[2]);
/* constants definition */
NTOT = n[0]*n[1]*n[2];
ntot = NTOT + 1;
NDIM = 3;
NMAX = n[0];
Ngrid = (*grid).NX*(*grid).NY*(*grid).NZ;
for (i=1;i<3;i++) {
if (n[i] > NMAX) NMAX = n[i];
if (n[i] == 1) NDIM = NDIM-1;
}
nmax = NMAX+1;
printf("NTOT = %d, ntot = %d, NMAX = %d\n",NTOT,ntot,NMAX);
/* wave vector computation */
normAxes(dir,nDir);
printf("Derivation direction (normed) %f %f %f\n",nDir[0],nDir[1],nDir[2]);
waveVect = (double *) malloc(Ngrid*sizeof(double));
if (waveVect == NULL) {
printf("waveVectCompute.cpp : No memory available for waveVect\n");
return;
}
waveVectorCompute3D((*grid).NX,(*grid).NY,(*grid).NZ,nDir,waveVect);
/* memory allocation */
table = (double *) malloc(ntot * sizeof(double));
if (table == NULL) {
printf("derivReal.cpp: No memory availble for table\n");
return;
}
realization = (double *) malloc(ntot * sizeof(double));
if (realization == NULL) {
printf("derivReal.cpp: No memory availble for realization\n");
return;
}
ExtendedWaveVect = (double *) malloc(ntot * sizeof(double));
if (ExtendedWaveVect == NULL) {
printf("derivReal.cpp: No memory availble for realization\n");
return;
}
workr = (double *) malloc(nmax * sizeof(double));
if (workr == NULL) {
printf("derivReal.cpp: No memory available for workr\n");
return;
}
worki = (double *) malloc(nmax * sizeof(double));
if (worki == NULL) {
printf("derivReal.cpp: No memory available for worki\n");
return;
}
extMaille =0;
/* organization of the extended realization */
for (k=1;k<=n[2];k++) {
for (j=1;j<=n[1];j++) {
for (i=1;i<=n[0];i++) {
extMaille = i + ((j-1) + (((k-1) * n[1]) * n[0]));
if (i <= (*grid).NX && j <= (*grid).NY && k <= (*grid).NZ) {
maille = i-1 + ((j-1) + ((k-1) * (*grid).NY) * (*grid).NX);
realization[extMaille] = (*Z).vector[maille];
ExtendedWaveVect[extMaille] = waveVect[maille];
} else {
realization[extMaille] = 0.0;
ExtendedWaveVect[extMaille] = 0.0;
}
}
}
}
/* forward fourier transform of the realization */
fourt(realization,table,n,NDIM,1,0,workr,worki);
FILE *wave;
wave = fopen("/home/irsrvshare1/R03/UPS_FLEX/waveVector.eas","w");
for (i=1;i<ntot;i++) {
tmp = realization[i];
realization[i] = - ExtendedWaveVect[i]*table[i];
table[i] = - ExtendedWaveVect[i]*tmp;
fprintf(wave,"%f\n", ExtendedWaveVect[i]);
}
fclose(wave);
free(waveVect);
/* backward fourier tranform */
fourt(realization,table,n,NDIM,0,1,workr,worki);
/* free the memory */
free(table);
free(workr);
free(worki);
printf("after second fourt \n");
(*dZ).n = (*Z).n;
(*dZ).code = (*Z).code;
(*dZ).vector = (double *) malloc(Ngrid*sizeof(double));
FILE *real;
FILE *derReal;
real = fopen("/home/irsrvshare1/R03/UPS_FLEX/realization.eas","w");
derReal = fopen("/home/irsrvshare1/R03/UPS_FLEX/derivative.eas","w");
printf("before saving file \n");
for(k=1;k<=(*grid).NZ;k++) {
for(j=1;j<=(*grid).NY;j++) {
for(i=1;i<=(*grid).NX;i++) {
extMaille = i+(j-1+(k-1)*n[1])*n[0];
maille = i-1+(j-1+(k-1)*(*grid).NY)*(*grid).NX;
(*dZ).vector[maille] = realization[extMaille]/(double)NTOT;
// printf("Z[%d] = %f, dZ[%d] = %f\n",maille,(*Z).vector[maille],maille,(*dZ).vector[maille]);
fprintf(real,"%f\n",(*Z).vector[maille]);
fprintf(derReal,"%f\n",(*dZ).vector[maille]);
}
}
}
printf("after saving file \n");
fclose(real);
fclose(derReal);
free(realization);
}
void normAxes(double *vec, double *normed) {
int i;
double r;
r = sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]);
if (normed == NULL) {
normed = (double *)malloc(3*sizeof(double));
if (normed == NULL) {
printf("derivReal.c: No memory available for normed\n");
}
}
for (i=0;i<3;i++) {
normed[i] = vec[i]/r;
}
}

@ -1,33 +0,0 @@
/*calculates C.x/
/* C : symetric positive-definite matrix recorded */
/* (per raws) as a vector with only components */
/* Cij so that j <= i, i = [0...n-1]*/
/* x : vector, xi avec i = [0...n-1]*/
/* b : vector, bi avec i = [0...n-1]*/
/* n : dimension of matrix Cij*/
/* */
/* The solution vector is returned in b*/
void mat_vec(double *C, double *x, double *b, int n)
{
int i,k,l;
double temp;
for (i = 0; i < n; i++) {
temp = 0.;
for (k = 0; k < n; k++) {
if ( k <= i ) {
l = k+i*(i+1)/2;
} else {
l = i+k*(k+1)/2;
}
temp += C[l]*x[k];
}
b[i] = temp;
}
return;
}

@ -1,48 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include "genlib.h"
#include "simpio.h"
#include "geostat.h"
#include "pressure.h"
#include "toolsIO.h"
void pgeneration(int n[3],struct grid_mod grid,struct statistic_mod stat,struct vario_mod variogram,string filename[7],struct pressure_mod pression,struct realization_mod *Y,struct realization_mod *P,struct realization_mod *VX,struct realization_mod *VY,struct realization_mod *VZ, int solver)
{
int i,ntot;
struct realization_mod GP;
FFTPressure(n,grid,Y,stat,pression,P,VX,VY,VZ,solver);
/*save the delta-pressure realization*/
writefile(filename[2],P);
total_pressure(grid,pression,P);
total_velocity(grid,stat.mean[0],1,pression,VX);
total_velocity(grid,stat.mean[0],2,pression,VY);
total_velocity(grid,stat.mean[0],3,pression,VZ);
/*save the total pressure realization*/
writefile(filename[3],P);
/*save the x-velocity realization*/
writefile(filename[4],VX);
/*save the y-velocity realization*/
writefile(filename[5],VY);
/*save the z-velocity realization*/
writefile(filename[6],VZ);
ntot=grid.NX*grid.NY*grid.NZ;
GP.vector = (double *) malloc(ntot * sizeof(double));
testmemory(GP.vector);
GP.n= ntot-1;
for (i = 0; i < ntot-1; i++)
{
GP.vector[i]=(*P).vector[i]-(*P).vector[i-1];
}
/*save the pressure gradient realization*/
writefile(filename[7],&GP);
return;
}

@ -1,99 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include "genlib.h"
#include "simpio.h"
#include "geostat.h"
#include "pressure.h"
#include "toolsIO.h"
void pgeneration2(int n[3],struct grid_mod grid,struct statistic_mod stat,struct vario_mod variogram,string filename[7],struct pressure_mod pression,struct realization_mod *Y,struct realization_mod *P,struct realization_mod *VX,struct realization_mod *VY,struct realization_mod *VZ, int solver, int format_file)
{
int i,ntot;
struct realization_mod GP;
FFTPressure(n,grid,Y,stat,pression,P,VX,VY,VZ,solver);
/*save the delta-pressure realization*/
switch (format_file)
{
case 0:
writefile(filename[2],P);
break;
case 1:
writefile_bin(filename[2],P);
break;
}
total_pressure(grid,pression,P);
total_velocity(grid,stat.mean[0],1,pression,VX);
total_velocity(grid,stat.mean[0],2,pression,VY);
total_velocity(grid,stat.mean[0],3,pression,VZ);
/*save the total pressure realization*/
switch (format_file)
{
case 0:
writefile(filename[3],P);
break;
case 1:
writefile_bin(filename[3],P);
break;
}
/*save the x-velocity realization*/
switch (format_file)
{
case 0:
writefile(filename[4],VX);
break;
case 1:
writefile_bin(filename[4],VX);
break;
}
/*save the y-velocity realization*/
switch (format_file)
{
case 0:
writefile(filename[5],VY);
break;
case 1:
writefile_bin(filename[5],VY);
break;
}
/*save the z-velocity realization*/
switch (format_file)
{
case 0:
writefile(filename[6],VZ);
break;
case 1:
writefile_bin(filename[6],VZ);
break;
}
ntot=grid.NX*grid.NY*grid.NZ;
GP.vector = (double *) malloc(ntot * sizeof(double));
testmemory(GP.vector);
GP.n= ntot-1;
for (i = 0; i < ntot-1; i++)
{
GP.vector[i]=(*P).vector[i]-(*P).vector[i-1];
}
/*save the pressure gradient realization*/
switch (format_file)
{
case 0:
writefile(filename[7],&GP);
break;
case 1:
writefile_bin(filename[7],&GP);
break;
}
return;
}

@ -1,59 +0,0 @@
#include <string.h>
#include "geostat.h"
#include "pressure.h"
#ifndef _TOOLSFFTPSIM_H
#define _TOOLSFFTPSIM_H
/* Create december, the 16th 2002 */
/* List of functions: */
/* ------------------ */
/* pgeneration, FFTPSim, FFTPressure, build_pressure, build_velocity,total_pressure,total_velocity, clean_real2 */
/*Functions */
/*----------*/
void pgeneration(int n[3],struct grid_mod grid,struct statistic_mod stat,struct vario_mod variogram,string filename[8],struct pressure_mod pression,struct realization_mod *Y,struct realization_mod *P,struct realization_mod *VX,struct realization_mod *VY,struct realization_mod *VZ, int solver);
void pgeneration2(int n[3],struct grid_mod grid,struct statistic_mod stat,struct vario_mod variogram,string filename[8],struct pressure_mod pression,struct realization_mod *Y,struct realization_mod *P,struct realization_mod *VX,struct realization_mod *VY,struct realization_mod *VZ, int solver, int format_file);
void FFTPSim(struct vario_mod variogram,struct statistic_mod stat,struct grid_mod grid,int n[3],struct realization_mod *realin,struct pressure_mod gradient,struct realization_mod *realout,struct realization_mod *realout2,struct realization_mod *realout3,struct realization_mod *realout4,struct realization_mod *realout5);
void FFTPressure(int n[3],struct grid_mod grid,struct realization_mod *realin,struct statistic_mod stat,struct pressure_mod gradient,struct realization_mod *realout,struct realization_mod *realout2,struct realization_mod *realout3,struct realization_mod *realout4, int solver);
void build_pressure(int n[3],struct grid_mod grid,struct pressure_mod gradient,double *realization,double *ireal,double *pressure,double *ipressure);
void build_velocity(int n[3],struct grid_mod grid,struct statistic_mod stat,struct pressure_mod gradient,double *realization,double *ireal,double *xvelocity,double *ixvelocity,int direction);
/* total_pressure */
/* Build total pressure */
/* grid: structure defining the grid */
/* pression:: structure defining the gradient pressure */
/* realout: structure defining a realization */
void total_pressure(struct grid_mod grid,struct pressure_mod gradient,struct realization_mod *realout);
/* total_velocity */
/* Build total velocity in one direction */
/* grid: structure defining the grid */
/* mean: permeability mean */
/* dgradient: macroscopic pressure gradient in wanted direction */
/* realout: structure defining a X velocity realization */
void total_velocity(struct grid_mod grid,double mean,int direction,struct pressure_mod gradient,struct realization_mod *realout);
void clean_real2(struct realization_mod *realin,int n[3],struct grid_mod grid,int solver,double *vectorresult,struct realization_mod *realout);
void normAxes(double *vec, double *normed);
void waveVectorCompute1D(int n,double *vec);
void waveVectorCompute3D(int nX,int nY, int nZ, /*float dX, float dY, float dZ,*/ double nDir[3], double *waveVect);
void derivReal(struct realization_mod *Z, struct realization_mod *dZ, double *dir, struct grid_mod *grid, struct vario_mod vario);
void mat_vec(double *C, double *x, double *b, int n);
#endif // define _TOOLSFFTPSIM_H

@ -1,35 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <stdarg.h>
#include <stdlib.h>
#include "geostat.h"
#include "pressure.h"
/* total_pressure */
/* Build total pressure */
/* grid: structure defining the grid */
/* pression: structure defining the gradient pressure */
/* realout: structure defining a realization */
void total_pressure(struct grid_mod grid,struct pressure_mod gradient,struct realization_mod *realout)
{
int i,j,k,maille0,maille1;
double temp,reference;
reference=abs(gradient.x)*grid.NX*grid.DX+abs(gradient.y)*grid.NY*grid.DY+abs(gradient.z)*grid.NZ*grid.DZ;
/* pressure reconstruction */
for ( k = 0; k < grid.NZ; k++) {
for (j = 0; j < grid.NY; j++) {
for (i = 0; i < grid.NX; i++) {
maille1 = i+(j+k*grid.NY)*grid.NX;
/* temp=reference+gradient.x*i*grid.DX+gradient.y*j*grid.DY+gradient.z*k*grid.DZ+(*realout).vector[maille1]; */
temp=-gradient.x*(i+1)/grid.NX-gradient.y*(j+1)/grid.NY-gradient.z*(k+1)/grid.NZ+(*realout).vector[maille1];
(*realout).vector[maille1]=temp;
}
}
}
return;
}

@ -1,50 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <stdarg.h>
#include <stdlib.h>
#include "geostat.h"
#include "pressure.h"
/* total_velocity */
/* Build total velocity */
/* grid: structure defining the grid */
/* mean: permeability mean */
/* dgradient: macroscopic pressure gradient in wanted direction */
/* realout: structure defining a X velocity realization */
void total_velocity(struct grid_mod grid,double mean,int direction,struct pressure_mod gradient,struct realization_mod *realout)
{
int i,j,k,maille0,maille1;
double temp,grad;
switch(direction)
{
case 1:
grad = gradient.x/(grid.NX*grid.DX);
break;
case 2:
grad = gradient.y/(grid.NY*grid.DY);/* *grid.NY*grid.DY; */
break;
case 3:
grad = gradient.z/(grid.NZ*grid.DZ); /* *grid.NZ*grid.DZ; */
break;
default :
printf("build_velocity.c: wrong velocity direction!!! direction: %d\n",direction);
break;
}
/* x velocity reconstruction */
for ( k = 0; k < grid.NZ; k++) {
for (j = 0; j < grid.NY; j++) {
for (i = 0; i < grid.NX; i++) {
maille1 = i+(j+k*grid.NY)*grid.NX;
temp=mean*grad+(*realout).vector[maille1];
/* temp=mean*grad; */
(*realout).vector[maille1]=temp;
}
}
}
return;
}

@ -1,156 +0,0 @@
//#include <fstream.h>
//#include <iostream>
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include <math.h>
//using namespace std;
/*Private functions */
void waveVectorCompute1D(int n,double *vec);
void waveVectorCompute3D(int nX,int nY, int nZ, /*double dX, double dY, double dZ,*/ double nDir[3], double *waveVect) {
int nTot;
double *vecX,*vecY,*vecZ;
double *waveVect3DX, *waveVect3DY, *waveVect3DZ;
int i, j, k, maille;
nTot = nX*nY*nZ;
// printf("nX = %d, nY = %d, nZ = %d, nTot = %d\n",nX,nY,nZ,nTot);
// printf("Entering waveVectorCompute3D\n");
/*Memory allocation */
waveVect3DX = (double *)malloc(nTot * sizeof(double));
if (waveVect3DX == NULL) {
printf("waveVectCompute.cpp : No memory available for waveVect3DX\n");
return;
}
waveVect3DY = (double *)malloc(nTot * sizeof(double));
if (waveVect3DY == NULL) {
printf("waveVectCompute.cpp : No memory available for waveVect3DY\n");
return;
}
waveVect3DZ = (double *)malloc(nTot * sizeof(double));
if (waveVect3DZ == NULL) {
printf("waveVectCompute.cpp : No memory available for waveVect3DZ\n");
return;
}
vecX = (double *)malloc(nX * sizeof(double));
if (vecX == NULL) {
printf("waveVectCompute.cpp : No memory available for vecX\n");
return;
}
vecY = (double *)malloc(nY * sizeof(double));
if (vecY == NULL) {
printf("waveVectCompute.cpp : No memory available for vecY\n");
return;
}
vecZ = (double *)malloc(nZ * sizeof(double));
if (vecZ == NULL) {
printf("waveVectCompute.cpp : No memory available for vecZ\n");
return;
}
// printf("memory allocation done.\n");
waveVectorCompute1D(nX,vecX);
waveVectorCompute1D(nY,vecY);
waveVectorCompute1D(nZ,vecZ);
// printf("waveVectorCompute1D: done!\n");
for(k = 0; k < nZ; k++) {
// printf("k=%d\n",k);
for(j = 0; j < nY; j++) {
// printf("j=%d: ",j);
for(i = 0; i < nX; i++) {
maille = i + (j + (k)*nY)*nX;
// waveVect3DX[maille] = vecX[i]/(nX*dX);
// waveVect3DY[maille] = vecY[j]/(nY*dY);
// waveVect3DZ[maille] = vecZ[k]/(nZ*dZ);
waveVect3DX[maille] = vecX[i];
waveVect3DY[maille] = vecY[j];
waveVect3DZ[maille] = vecZ[k];
// printf("i=%d: %f ",i,vecZ[k]);
}
// printf("\n");
}
}
// printf("waveVector3DX, Y and Z rearanged\n");
free(vecX);
free(vecY);
free(vecZ);
if (waveVect == NULL) {
printf("allocate memory for waveVect\n");
waveVect = (double *) malloc((nTot+1)*sizeof(double));
if (waveVect == NULL) {
printf("waveVectCompute.cpp : No memory available for waveVect\n");
return;
}
printf("after allocate memory for waveVect\n");
}
for (i=1;i<=nTot;i++) {
waveVect[i] = nDir[0]*waveVect3DX[i-1] + nDir[1]*waveVect3DY[i-1] + nDir[2]*waveVect3DZ[i-1];
// printf("waveVect[%d] = %f\n",i,waveVect[i]);
}
// printf("waveVector rearanged\n");
free(waveVect3DX);
free(waveVect3DY);
free(waveVect3DZ);
}
void waveVectorCompute1D(int n,double *vec){
int midPoint,k;
if(ceil(n/2) - floor(n/2) < 1.) {
// printf("coucou\n");
/*n is even */
midPoint = (int) floor(n/2);
vec[midPoint] = (double) midPoint;
} else {
/*n is odd */
midPoint = (int) floor(n/2)-1;
vec[midPoint] = (double) midPoint+1;
}
for (k=1;k<midPoint;k++) {
vec[k] = (double) k;
vec[n-k] = (double) -(k);
}
vec[0] = (double) 0;
}
/*int main() {
double *vecX;
int k;
vecX = (double *)malloc(10 * sizeof(double));
if (vecX == NULL) {
printf("waveVectCompute.cpp : No memory available for vecX\n");
return(0);
}
printf("avant\n");
waveVectorCompute1D(10,vecX);
printf("apres\n");
ofstream vec;
vec.open("waveVect1D.pol");
for (k=0;k<10;k++)
vec << vecX[k] << "\n";
vec.close();
return(1);
}*/

@ -1,31 +0,0 @@
#
########################################################################
#
# Makefile for library
#
########################################################################
#
#INTERFACE = ${BOSSE}/CODES/LIBS_C/Interface
INTERFACE = ../include
INCLUDE = -I${INTERFACE}
CC = cc
CCFLAGS =
LIB = libIO_${ARCH}.a
NOBJS= inputdata.o inputfiledata.o debuginput.o readdata.o readfile_bin.o writefile.o writefile_bin.o testmemory.o testopenfile.o readdata2.o readdata3.o
.c.o:
${CC} $(INCLUDE) $(CCFLAGS) -c $<
# LIBRARY
$(LIB) : $(NOBJS)
ar cr $(LIB) $(NOBJS)
ranlib $(LIB)
clean :
rm *.o

@ -1,16 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include <math.h>
void allouememoire(double *realint, string variable)
{
if (realint == NULL) {
printf("Testmemory.c: No memory available for %s\n",variable);
exit;
}
return;
}

@ -1,73 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include "genlib.h"
#include "geostat.h"
#include "toolsIO.h"
/* DebugInput */
/* */
/* Display the input data */
/* seed: seed */
/* grid: structure defining the actual grid */
/* filename: name of files */
/* variogram: struture defining the variogram model */
/* stat: struture defining the mean and the variance of permeability realization */
/* pression: structure defining the gradient pressure */
void debuginput(long *seed,struct grid_mod *grid,string filename[8],struct vario_mod *variogram,struct statistic_mod *stat,struct pressure_mod *pression,int *Ksolver,int *genere, int *gwnwrite)
{
int i;
/* debug du 8/7/2002 */
printf("\n\n");
printf("Starting seed (integer): %d\n",(*seed));
printf("Number of cells along the X axis: %d\n",(*grid).NX);
printf("Number of cells along the Y axis: %d\n",(*grid).NY);
printf("Number of cells along the Z axis: %d\n",(*grid).NZ);
printf("cell length along the X axis: %6.4f\n",(*grid).DX);
printf("cell length along the Y axis: %6.4f\n",(*grid).DY);
printf("cell length along the Z axis: %6.4f\n",(*grid).DZ);
printf("Number of structures for the variogram: %d\n",(*variogram).Nvario);
i=0;
printf("Weight: %6.4f\n",(*variogram).var[i]);
printf("Type of variogram: %d\n",(*variogram).vario[i]);
printf("Exponent: %6.4f\n",(*variogram).alpha[i]);
printf("Mean of the output realization: %15.8f\n",(*stat).mean[0]);
printf("Variance of the output realization: %6.4f\n",(*stat).variance[0]);
printf("Structure of the field (0-normal case 1-lognormal case 2-log10 case) : %d\n",(*stat).type);
/*output files*/
printf("output filename for permeability realization: %s\n", filename[1]);
printf("\n\n");
if (*gwnwrite == 0)
{
printf("K field generation with Gaussian white noise!\n");
printf("output filename for Gaussian white noise: %s\n",filename[0]);
}
printf("\n\n");
if ((*Ksolver == 1) | (*Ksolver == 2))
{
printf("P field generation: %d/n",*Ksolver);
printf("output filename for pressure realization: %s\n", filename[2]);
printf("output filename for pressure total realization: %s\n", filename[3]);
printf("output filename for x-velocity realization: %s\n", filename[4]);
printf("output filename for y-velocity realization: %s\n", filename[5]);
printf("output filename for z-velocity realization: %s\n", filename[6]);
printf("output filename for pressure gradient realization: %s\n", filename[7]);
/*Pressure data*/
printf("Pressure gradient in x direction: %6.4f\n",(*pression).x);
printf("Pressure gradient in y direction: %6.4f\n",(*pression).y);
printf("Pressure gradient in z direction: %6.4f\n",(*pression).z);
}
printf("\n\n");
return;
}

@ -1,122 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include "genlib.h"
#include "simpio.h"
#include "geostat.h"
#include "toolsIO.h"
/* Inputdata */
/* */
/* input data needed for realizations*/
/* seed: seed */
/* grid: structure defining the actual grid */
/* filename: name of files */
/* variogram: struture defining the variogram model */
/* stat: struture defining the mean and the variance of permeability realization */
/* pression: structure defining the gradient pressure */
void inputdata(long *seed,struct grid_mod *grid,string filename[7],struct vario_mod *variogram,struct statistic_mod *stat,struct pressure_mod *pression)
{
int i,j;
/*seed*/
printf("Starting seed (integer): \n");
*seed = GetInteger();
/*Grid description*/
printf("Number of cells along the X axis: \n");
(*grid).NX = GetInteger();
printf("Number of cells along the Y axis: \n");
(*grid).NY = GetInteger();
printf("Number of cells along the Z axis: \n");
(*grid).NZ = GetInteger();
printf("cell length along the X axis: \n");
(*grid).DX = GetReal();
printf("cell length along the Y axis: \n");
(*grid).DY = GetReal();
printf("cell length along the Z axis: \n");
(*grid).DZ = GetReal();
/*output file*/
printf("output filename for Gaussian white noise: \n");
filename[0] = GetLine();
/*variogram description*/
printf("Number of structures for the variogram: \n");
(*variogram).Nvario = GetInteger();
(*variogram).vario = (int *) malloc((*variogram).Nvario * sizeof(int));
(*variogram).alpha = (double *) malloc((*variogram).Nvario * sizeof(double));
(*variogram).ap = (double *) malloc(9*(*variogram).Nvario * sizeof(double));
(*variogram).scf = (double *) malloc(3*(*variogram).Nvario * sizeof(double));
(*variogram).var = (double *) malloc((*variogram).Nvario * sizeof(double));
for (i= 0; i < (*variogram).Nvario; i++) {
printf("i %d\n",i);
printf("Weight: \n");
(*variogram).var[i] = GetReal();
printf("Type of variogram: \n");
(*variogram).vario[i] = GetInteger();
printf("Exponent: \n");
(*variogram).alpha[i] = GetReal();
printf("Correlation lengths (3): \n");
for (j = 0; j < 3; j++)
(*variogram).scf[i*3+j]=GetReal();
printf("Coordinates of the first two main axes (first axis first, then second): \n");
for (j = 0; j < 6; j++)
(*variogram).ap[i*9+j] = GetReal();
}
/*statistical data*/
(*stat).nblock_mean = 1;
(*stat).nblock_var = 1;
(*stat).mean = (double *)malloc((*stat).nblock_mean * sizeof(double));
if ((*stat).mean == NULL)
Error("No memory available");
(*stat).variance = (double *)malloc((*stat).nblock_var * sizeof(double));
if ((*stat).variance == NULL)
Error("No memory available");
printf("Mean of the output realization: \n");
(*stat).mean[0] = GetReal();
printf("Variance of the output realization: \n");
(*stat).variance[0] = GetReal();
printf("Structure of the field (0-normal case 1-lognormal case 2-log10 case) : \n");
(*stat).type = GetInteger();
/*output file*/
printf("output filename for permeability realization: \n");
filename[1] = GetLine();
/*Pressure data*/
printf("Macroscopic pressure gradient in x direction: \n");
(*pression).x = GetReal();
printf("Macroscopic pressure gradient in y direction: \n");
(*pression).y = GetReal();
printf("Macroscopic pressure gradient in z direction: \n");
(*pression).z = GetReal();
/*output pressure file*/
printf("output filename for pressure realization: \n");
filename[2] = GetLine();
/*output pressure totale file*/
printf("output filename for total pressure realization: \n");
filename[3] = GetLine();
/*output x-velocity file*/
printf("output filename for x-velocity realization: \n");
filename[4] = GetLine();
/*output y-velocity file*/
printf("output filename for y-velocity realization: \n");
filename[5] = GetLine();
/*output z-velocity file*/
printf("output filename for z-velocity realization: \n");
filename[6] = GetLine();
return;
}

@ -1,71 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include "genlib.h"
#include "simpio.h"
#include "geostat.h"
#include "toolsIO.h"
void inputfiledata(string inputfile,long *seed,struct grid_mod *grid,string filename[7],struct vario_mod *variogram,struct statistic_mod *stat,struct pressure_mod *pression)
{
FILE *fp;
int i,j;
fp=fopen(inputfile,"r");
if(fp== NULL)
{
printf("Erreur d'ouverture du fichier\n");
exit(0);
}
*seed=atoi(ReadLine(fp));
(*grid).NX = atoi(ReadLine(fp));
(*grid).NY = atoi(ReadLine(fp));
(*grid).NZ = atoi(ReadLine(fp));
(*grid).DX = atof(ReadLine(fp));
(*grid).DY = atof(ReadLine(fp));
(*grid).DZ = atof(ReadLine(fp));
filename[0] = ReadLine(fp);
(*variogram).Nvario = atoi(ReadLine(fp));
(*variogram).vario = (int *) malloc((*variogram).Nvario * sizeof(int));
(*variogram).alpha = (double *) malloc((*variogram).Nvario * sizeof(double));
(*variogram).ap = (double *) malloc(9*(*variogram).Nvario * sizeof(double));
(*variogram).scf = (double *) malloc(3*(*variogram).Nvario * sizeof(double));
(*variogram).var = (double *) malloc((*variogram).Nvario * sizeof(double));
for (i= 0; i < (*variogram).Nvario; i++) {
(*variogram).var[i] = atof(ReadLine(fp));
(*variogram).vario[i] = atoi(ReadLine(fp));
(*variogram).alpha[i] = atof(ReadLine(fp));
for (j = 0; j < 3; j++)
(*variogram).scf[i*3+j]= atof(ReadLine(fp));
for (j = 0; j < 6; j++)
(*variogram).ap[i*9+j] = atof(ReadLine(fp));
}
/*statistical data*/
(*stat).nblock_mean = 1;
(*stat).nblock_var = 1;
(*stat).mean = (double *)malloc((*stat).nblock_mean * sizeof(double));
if ((*stat).mean == NULL)
Error("No memory available");
(*stat).variance = (double *)malloc((*stat).nblock_var * sizeof(double));
if ((*stat).variance == NULL)
Error("No memory available");
(*stat).mean[0] = atof(ReadLine(fp));
(*stat).variance[0] = atof(ReadLine(fp));
(*stat).type = atoi(ReadLine(fp));
filename[1] = ReadLine(fp);
(*pression).x = atof(ReadLine(fp));
(*pression).y = atof(ReadLine(fp));
(*pression).z = atof(ReadLine(fp));
filename[2] = ReadLine(fp);
filename[3] = ReadLine(fp);
filename[4] = ReadLine(fp);
filename[5] = ReadLine(fp);
filename[6] = ReadLine(fp);
fclose(fp);
return;
}

@ -1,21 +0,0 @@
#include <stdlib.h>
#include <string.h>
#include "genlib.h"
#include "geostat.h"
#ifndef _PRESSURE_H
#define _PRESSURE_H
/*STRUCTURES*/
/*----------*/
/* pressure_mod */
/* x: macroscopic gradient value in x direction */
/* y: macroscopic gradient value in x direction */
/* z: macroscopic gradient value in x direction */
struct pressure_mod {
double x,y,z;
};
#endif // define _PRESSURE_H

@ -1,204 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include "genlib.h"
#include "simpio.h"
#include "geostat.h"
#include "pressure.h"
#include "toolsIO.h"
void readdata(long *seed,struct grid_mod *grid,string filename[8],struct vario_mod *variogram,struct statistic_mod *stat,struct pressure_mod *pression,int *Ksolver,int *genere, int *gwnwrite, struct realization_mod *Kfield,struct realization_mod *gwnoise, char *argv[])
{
FILE *fp;
int i,j,n;
char *file1,*file2,*file3,*file4,*file5;
char *prog=argv[0];
double tmp;
file1=argv[1];
file2=argv[2];
/* Ouverture du fichier de données principal */
if ((fp=fopen(file1,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file1);
exit(1);
}
else
{
*Ksolver = atoi(ReadLine(fp));
*genere = atoi(ReadLine(fp));
*gwnwrite = atoi(ReadLine(fp));
if (*gwnwrite==0)
{
filename[0] = ReadLine(fp);
}
(*grid).NX = atoi(ReadLine(fp));
(*grid).NY = atoi(ReadLine(fp));
(*grid).NZ = atoi(ReadLine(fp));
(*grid).DX = atof(ReadLine(fp));
(*grid).DY = atof(ReadLine(fp));
(*grid).DZ = atof(ReadLine(fp));
fclose(fp);
}
n=(*grid).NX*(*grid).NY*(*grid).NZ;
/* Ouverture du fichier de données sur le champ de perméabilité */
if ((fp=fopen(file2,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file2);
exit(1);
}
else
{
*seed=atoi(ReadLine(fp));
(*variogram).Nvario = atoi(ReadLine(fp));
(*variogram).vario = (int *) malloc((*variogram).Nvario * sizeof(int));
(*variogram).alpha = (double *) malloc((*variogram).Nvario * sizeof(double));
(*variogram).ap = (double *) malloc(9*(*variogram).Nvario * sizeof(double));
(*variogram).scf = (double *) malloc(3*(*variogram).Nvario * sizeof(double));
(*variogram).var = (double *) malloc((*variogram).Nvario * sizeof(double));
for (i= 0; i < (*variogram).Nvario; i++) {
(*variogram).var[i] = atof(ReadLine(fp));
(*variogram).vario[i] = atoi(ReadLine(fp));
(*variogram).alpha[i] = atof(ReadLine(fp));
for (j = 0; j < 3; j++)
(*variogram).scf[i*3+j]= atof(ReadLine(fp));
for (j = 0; j < 6; j++)
(*variogram).ap[i*9+j] = atof(ReadLine(fp));
}
/*statistical data*/
(*stat).nblock_mean = 1;
(*stat).nblock_var = 1;
(*stat).mean = (double *)malloc((*stat).nblock_mean * sizeof(double));
if ((*stat).mean == NULL)
Error("No memory available");
(*stat).variance = (double *)malloc((*stat).nblock_var * sizeof(double));
if ((*stat).variance == NULL)
Error("No memory available");
(*stat).mean[0] = atof(ReadLine(fp));
(*stat).variance[0] = atof(ReadLine(fp));
(*stat).type = atoi(ReadLine(fp));
filename[1] = ReadLine(fp);
fclose(fp);
}
switch (*Ksolver)
{
case 0:
/* Ouverture du fichier contenant le gaussian white noise */
if (*genere == 1)
{
file3=argv[3];
if ((fp=fopen(file3,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file3);
exit(1);
}
else
{
/* Ouverture du champ de permeabilite K */
(*gwnoise).vector = (double *) malloc(n * sizeof(double));
for (i=0;i<n;i++)
{
/* (*gwnoise).vector[i]=atof(ReadLine(fp)); */
/* sscanf(ReadLine(fp)," %Lf %c",gwnoise.vector[i],&merde); */
fscanf(fp,"%lf",&tmp);
(*gwnoise).vector[i]= tmp;
}
fclose(fp);
}
}
break;
case 1:
file3=argv[3];
/* Ouverture du fichier de données de pression*/
if ((fp=fopen(file3,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file3);
exit(1);
}
else
{
/* testopenfile(fp); */
filename[2] = ReadLine(fp);
filename[3] = ReadLine(fp);
(*pression).x = atof(ReadLine(fp));
(*pression).y = atof(ReadLine(fp));
(*pression).z = atof(ReadLine(fp));
filename[4] = ReadLine(fp);
filename[5] = ReadLine(fp);
filename[6] = ReadLine(fp);
filename[7] = ReadLine(fp);
fclose(fp);
}
/* Ouverture du champ de permeabilite K */
(*Kfield).vector = (double *) malloc(n * sizeof(double));
file4=argv[4];
if ((fp=fopen(file4,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file4);
exit(1);
}
else
{
for (i=0;i<n;i++)
(*Kfield).vector[i]=atof(ReadLine(fp));
fclose(fp);
}
break;
case 2:
/* Ouverture du fichier de données de pression */
file3=argv[3];
if ((fp=fopen(file3,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file3);
exit(1);
}
else
{
/* testopenfile(fp); */
filename[2] = ReadLine(fp);
filename[3] = ReadLine(fp);
(*pression).x = atof(ReadLine(fp));
(*pression).y = atof(ReadLine(fp));
(*pression).z = atof(ReadLine(fp));
filename[4] = ReadLine(fp);
filename[5] = ReadLine(fp);
filename[6] = ReadLine(fp);
filename[7] = ReadLine(fp);
fclose(fp);
}
/* Ouverture du fichier contenant le gaussian white noise */
if (*genere == 1)
{
file4=argv[4];
if ((fp=fopen(file4,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file3);
exit(1);
}
else
{
/* Ouverture du champ de permeabilite K */
(*gwnoise).vector = (double *) malloc(n * sizeof(double));
for (i=0;i<n;i++)
(*gwnoise).vector[i]=atof(ReadLine(fp));
fclose(fp);
}
}
break;
}
return;
}

@ -1,205 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include "genlib.h"
#include "simpio.h"
#include "geostat.h"
#include "pressure.h"
#include "toolsIO.h"
void readdata2(long *seed,struct grid_mod *grid,string filename[8],struct vario_mod *variogram,struct statistic_mod *stat,struct pressure_mod *pression,int *Ksolver,int *genere, int *gwnwrite, int *format_file, struct realization_mod *Kfield,struct realization_mod *gwnoise, char *argv[])
{
FILE *fp;
int i,j,n;
char *file1,*file2,*file3,*file4,*file5;
char *prog=argv[0];
double tmp;
file1=argv[1];
file2=argv[2];
/* Ouverture du fichier de données principal */
if ((fp=fopen(file1,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file1);
exit(1);
}
else
{
*Ksolver = atoi(ReadLine(fp));
*genere = atoi(ReadLine(fp));
*gwnwrite = atoi(ReadLine(fp));
if (*gwnwrite==0)
{
filename[0] = ReadLine(fp);
}
*format_file=atoi(ReadLine(fp));
(*grid).NX = atoi(ReadLine(fp));
(*grid).NY = atoi(ReadLine(fp));
(*grid).NZ = atoi(ReadLine(fp));
(*grid).DX = atof(ReadLine(fp));
(*grid).DY = atof(ReadLine(fp));
(*grid).DZ = atof(ReadLine(fp));
fclose(fp);
}
n=(*grid).NX*(*grid).NY*(*grid).NZ;
/* Ouverture du fichier de données sur le champ de perméabilité */
if ((fp=fopen(file2,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file2);
exit(1);
}
else
{
*seed=atoi(ReadLine(fp));
(*variogram).Nvario = atoi(ReadLine(fp));
(*variogram).vario = (int *) malloc((*variogram).Nvario * sizeof(int));
(*variogram).alpha = (double *) malloc((*variogram).Nvario * sizeof(double));
(*variogram).ap = (double *) malloc(9*(*variogram).Nvario * sizeof(double));
(*variogram).scf = (double *) malloc(3*(*variogram).Nvario * sizeof(double));
(*variogram).var = (double *) malloc((*variogram).Nvario * sizeof(double));
for (i= 0; i < (*variogram).Nvario; i++) {
(*variogram).var[i] = atof(ReadLine(fp));
(*variogram).vario[i] = atoi(ReadLine(fp));
(*variogram).alpha[i] = atof(ReadLine(fp));
for (j = 0; j < 3; j++)
(*variogram).scf[i*3+j]= atof(ReadLine(fp));
for (j = 0; j < 6; j++)
(*variogram).ap[i*9+j] = atof(ReadLine(fp));
}
/*statistical data*/
(*stat).nblock_mean = 1;
(*stat).nblock_var = 1;
(*stat).mean = (double *)malloc((*stat).nblock_mean * sizeof(double));
if ((*stat).mean == NULL)
Error("No memory available");
(*stat).variance = (double *)malloc((*stat).nblock_var * sizeof(double));
if ((*stat).variance == NULL)
Error("No memory available");
(*stat).mean[0] = atof(ReadLine(fp));
(*stat).variance[0] = atof(ReadLine(fp));
(*stat).type = atoi(ReadLine(fp));
filename[1] = ReadLine(fp);
fclose(fp);
}
switch (*Ksolver)
{
case 0:
/* Ouverture du fichier contenant le gaussian white noise */
if (*genere == 1)
{
file3=argv[3];
if ((fp=fopen(file3,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file3);
exit(1);
}
else
{
/* Ouverture du champ de permeabilite K */
(*gwnoise).vector = (double *) malloc(n * sizeof(double));
for (i=0;i<n;i++)
{
/* (*gwnoise).vector[i]=atof(ReadLine(fp)); */
/* sscanf(ReadLine(fp)," %Lf %c",gwnoise.vector[i],&merde); */
fscanf(fp,"%lf",&tmp);
(*gwnoise).vector[i]= tmp;
}
fclose(fp);
}
}
break;
case 1:
file3=argv[3];
/* Ouverture du fichier de données de pression*/
if ((fp=fopen(file3,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file3);
exit(1);
}
else
{
/* testopenfile(fp); */
filename[2] = ReadLine(fp);
filename[3] = ReadLine(fp);
(*pression).x = atof(ReadLine(fp));
(*pression).y = atof(ReadLine(fp));
(*pression).z = atof(ReadLine(fp));
filename[4] = ReadLine(fp);
filename[5] = ReadLine(fp);
filename[6] = ReadLine(fp);
filename[7] = ReadLine(fp);
fclose(fp);
}
/* Ouverture du champ de permeabilite K */
(*Kfield).vector = (double *) malloc(n * sizeof(double));
file4=argv[4];
if ((fp=fopen(file4,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file4);
exit(1);
}
else
{
for (i=0;i<n;i++)
(*Kfield).vector[i]=atof(ReadLine(fp));
fclose(fp);
}
break;
case 2:
/* Ouverture du fichier de données de pression */
file3=argv[3];
if ((fp=fopen(file3,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file3);
exit(1);
}
else
{
/* testopenfile(fp); */
filename[2] = ReadLine(fp);
filename[3] = ReadLine(fp);
(*pression).x = atof(ReadLine(fp));
(*pression).y = atof(ReadLine(fp));
(*pression).z = atof(ReadLine(fp));
filename[4] = ReadLine(fp);
filename[5] = ReadLine(fp);
filename[6] = ReadLine(fp);
filename[7] = ReadLine(fp);
fclose(fp);
}
/* Ouverture du fichier contenant le gaussian white noise */
if (*genere == 1)
{
file4=argv[4];
if ((fp=fopen(file4,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file3);
exit(1);
}
else
{
/* Ouverture du champ de permeabilite K */
(*gwnoise).vector = (double *) malloc(n * sizeof(double));
for (i=0;i<n;i++)
(*gwnoise).vector[i]=atof(ReadLine(fp));
fclose(fp);
}
}
break;
}
return;
}

@ -1,206 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include "genlib.h"
#include "simpio.h"
#include "geostat.h"
#include "pressure.h"
#include "toolsIO.h"
void readdata3(long *seed,struct grid_mod *grid,string filename[8],struct vario_mod *variogram,struct statistic_mod *stat,struct pressure_mod *pression,int *Ksolver,int *genere, int *gwnwrite, int *format_file,int *Psolver, struct realization_mod *Kfield,struct realization_mod *gwnoise, char *argv[])
{
FILE *fp;
int i,j,n;
char *file1,*file2,*file3,*file4,*file5;
char *prog=argv[0];
double tmp;
file1=argv[1];
file2=argv[2];
/* Ouverture du fichier de données principal */
if ((fp=fopen(file1,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file1);
exit(1);
}
else
{
*Ksolver = atoi(ReadLine(fp));
*genere = atoi(ReadLine(fp));
*gwnwrite = atoi(ReadLine(fp));
if (*gwnwrite==0)
{
filename[0] = ReadLine(fp);
}
*format_file=atoi(ReadLine(fp));
*Psolver=atoi(ReadLine(fp));
(*grid).NX = atoi(ReadLine(fp));
(*grid).NY = atoi(ReadLine(fp));
(*grid).NZ = atoi(ReadLine(fp));
(*grid).DX = atof(ReadLine(fp));
(*grid).DY = atof(ReadLine(fp));
(*grid).DZ = atof(ReadLine(fp));
fclose(fp);
}
n=(*grid).NX*(*grid).NY*(*grid).NZ;
/* Ouverture du fichier de données sur le champ de perméabilité */
if ((fp=fopen(file2,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file2);
exit(1);
}
else
{
*seed=atoi(ReadLine(fp));
(*variogram).Nvario = atoi(ReadLine(fp));
(*variogram).vario = (int *) malloc((*variogram).Nvario * sizeof(int));
(*variogram).alpha = (double *) malloc((*variogram).Nvario * sizeof(double));
(*variogram).ap = (double *) malloc(9*(*variogram).Nvario * sizeof(double));
(*variogram).scf = (double *) malloc(3*(*variogram).Nvario * sizeof(double));
(*variogram).var = (double *) malloc((*variogram).Nvario * sizeof(double));
for (i= 0; i < (*variogram).Nvario; i++) {
(*variogram).var[i] = atof(ReadLine(fp));
(*variogram).vario[i] = atoi(ReadLine(fp));
(*variogram).alpha[i] = atof(ReadLine(fp));
for (j = 0; j < 3; j++)
(*variogram).scf[i*3+j]= atof(ReadLine(fp));
for (j = 0; j < 6; j++)
(*variogram).ap[i*9+j] = atof(ReadLine(fp));
}
/*statistical data*/
(*stat).nblock_mean = 1;
(*stat).nblock_var = 1;
(*stat).mean = (double *)malloc((*stat).nblock_mean * sizeof(double));
if ((*stat).mean == NULL)
Error("No memory available");
(*stat).variance = (double *)malloc((*stat).nblock_var * sizeof(double));
if ((*stat).variance == NULL)
Error("No memory available");
(*stat).mean[0] = atof(ReadLine(fp));
(*stat).variance[0] = atof(ReadLine(fp));
(*stat).type = atoi(ReadLine(fp));
filename[1] = ReadLine(fp);
fclose(fp);
}
switch (*Ksolver)
{
case 0:
/* Ouverture du fichier contenant le gaussian white noise */
if (*genere == 1)
{
file3=argv[3];
if ((fp=fopen(file3,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file3);
exit(1);
}
else
{
/* Ouverture du champ de permeabilite K */
(*gwnoise).vector = (double *) malloc(n * sizeof(double));
for (i=0;i<n;i++)
{
/* (*gwnoise).vector[i]=atof(ReadLine(fp)); */
/* sscanf(ReadLine(fp)," %Lf %c",gwnoise.vector[i],&merde); */
fscanf(fp,"%lf",&tmp);
(*gwnoise).vector[i]= tmp;
}
fclose(fp);
}
}
break;
case 1:
file3=argv[3];
/* Ouverture du fichier de données de pression*/
if ((fp=fopen(file3,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file3);
exit(1);
}
else
{
/* testopenfile(fp); */
filename[2] = ReadLine(fp);
filename[3] = ReadLine(fp);
(*pression).x = atof(ReadLine(fp));
(*pression).y = atof(ReadLine(fp));
(*pression).z = atof(ReadLine(fp));
filename[4] = ReadLine(fp);
filename[5] = ReadLine(fp);
filename[6] = ReadLine(fp);
filename[7] = ReadLine(fp);
fclose(fp);
}
/* Ouverture du champ de permeabilite K */
(*Kfield).vector = (double *) malloc(n * sizeof(double));
file4=argv[4];
if ((fp=fopen(file4,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file4);
exit(1);
}
else
{
for (i=0;i<n;i++)
(*Kfield).vector[i]=atof(ReadLine(fp));
fclose(fp);
}
break;
case 2:
/* Ouverture du fichier de données de pression */
file3=argv[3];
if ((fp=fopen(file3,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file3);
exit(1);
}
else
{
/* testopenfile(fp); */
filename[2] = ReadLine(fp);
filename[3] = ReadLine(fp);
(*pression).x = atof(ReadLine(fp));
(*pression).y = atof(ReadLine(fp));
(*pression).z = atof(ReadLine(fp));
filename[4] = ReadLine(fp);
filename[5] = ReadLine(fp);
filename[6] = ReadLine(fp);
filename[7] = ReadLine(fp);
fclose(fp);
}
/* Ouverture du fichier contenant le gaussian white noise */
if (*genere == 1)
{
file4=argv[4];
if ((fp=fopen(file4,"r")) == NULL)
{
fprintf(stderr,"%s :impossible d'ouvrir %s\n",prog, file3);
exit(1);
}
else
{
/* Ouverture du champ de permeabilite K */
(*gwnoise).vector = (double *) malloc(n * sizeof(double));
for (i=0;i<n;i++)
(*gwnoise).vector[i]=atof(ReadLine(fp));
fclose(fp);
}
}
break;
}
return;
}

@ -1,35 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include <math.h>
#include "genlib.h"
#include "geostat.h"
/*readfile */
/* read in the file "filename" the vector values of a */
/* realization_mod variable */
/* filename: explicit */
/* realin: structure defining a realization */
void readfile_bin(string filename, struct realization_mod *realin)
{
FILE *fp;
int i;
double prout;
/*read the permeability realization*/
fp = fopen(filename, "r");
for (i=0;i<(*realin).n;i++)
{
fread(&prout,sizeof(double),1,fp);
printf("Prout: %15.10f\n",prout);
}
fclose(fp);
return;
}

@ -1,15 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include <math.h>
void testmemory(double *realint)
{
if (realint == NULL) {
printf("Testmemory.c: No memory available \n");
exit;
}
return;
}

@ -1,16 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include <math.h>
void testopenfile(FILE *fp)
{
if (fp == NULL)
{
printf("Erreur d'ouverture de fichier\n");
exit(0);
}
return;
}

@ -1,93 +0,0 @@
#include <stdlib.h>
#include <string.h>
#include "genlib.h"
#include "geostat.h"
#include "pressure.h"
#ifndef _TOOLSIO_H
#define _TOOLSIO_H
/* Create december, the 16th 2002 */
/* Modified december, the 9th 2002 */
/* List of functions: */
/* ------------------ */
/* writefile,writefile_bin,readfile_bin,inputdata,inputfiledata,readdata,debuginput */
/* testmemory, testopenfile */
/*FUNCTIONS*/
/*----------*/
/* Lecture dans un fichier */
/* Inputdata */
/* */
/* input data needed for realizations*/
/* seed: seed */
/* grid: structure defining the actual grid */
/* filename: name of files */
/* variogram: struture defining the variogram model */
/* stat: struture defining the mean and the variance of permeability realization */
/* pression: structure defining the gradient pressure */
void inputdata(long *seed,struct grid_mod *grid,string filename[7],struct vario_mod *variogram,struct statistic_mod *stat,struct pressure_mod *pression);
void inputfiledata(string inputfile,long *seed,struct grid_mod *grid,string filename[7],struct vario_mod *variogram,struct statistic_mod *stat,struct pressure_mod *pression);
/* void readdata(long *seed,struct grid_mod *grid,string filename[8],struct vario_mod *variogram,struct statistic_mod *stat,struct pressure_mod *pression,int *Ksolver,struct realization_mod *Kfield, char *argv[]); */
void readdata(long *seed,struct grid_mod *grid,string filename[8],struct vario_mod *variogram,struct statistic_mod *stat,struct pressure_mod *pression,int *Ksolver,int *genere, int *gwnwrite, struct realization_mod *Kfield,struct realization_mod *gwnoise, char *argv[]);
void readdata2(long *seed,struct grid_mod *grid,string filename[8],struct vario_mod *variogram,struct statistic_mod *stat,struct pressure_mod *pression,int *Ksolver,int *genere, int *gwnwrite, int *format_file, struct realization_mod *Kfield,struct realization_mod *gwnoise, char *argv[]);
void readdata3(long *seed,struct grid_mod *grid,string filename[8],struct vario_mod *variogram,struct statistic_mod *stat,struct pressure_mod *pression,int *Ksolver,int *genere, int *gwnwrite, int *format_file,int *Psolver, struct realization_mod *Kfield,struct realization_mod *gwnoise, char *argv[]);
/* readfile_bin */
/* */
/* read in the file "filename" the vector values of a */
/* realization_mod variable */
/* filename: explicit */
/* realin: structure defining a realization */
void readfile_bin(string filename, struct realization_mod *realin);
/* Writefile */
/* */
/* write in the file "filename" the vector values of a */
/* realization_mod variable */
/* filename: explicit */
/* realin: structure defining a realization */
void writefile( string filename, struct realization_mod *realin);
/* Writefile_bin */
/* */
/* write in the file "filename" the vector values of a */
/* realization_mod variable */
/* filename: explicit */
/* realin: structure defining a realization */
void writefile_bin( string filename, struct realization_mod *realin);
/* DebugInput */
/* */
/* Display the input data */
/* seed: seed */
/* grid: structure defining the actual grid */
/* filename: name of files */
/* variogram: struture defining the variogram model */
/* stat: struture defining the mean and the variance of permeability realization */
/* pression: structure defining the gradient pressure */
void debuginput(long *seed,struct grid_mod *grid,string filename[7],struct vario_mod *variogram,struct statistic_mod *stat,struct pressure_mod *pression,int *Ksolver,int *genere, int *gwnwrite);
/* Allocation test */
void testmemory(double *realint);
/* Test open file */
void testopenfile(FILE *fp);
#endif // define _TOOLSIO_H

@ -1,32 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include <math.h>
#include "genlib.h"
#include "geostat.h"
/*writefile */
/* write in the file "filename" the vector values of a */
/* realization_mod variable */
/* filename: explicit */
/* realin: structure defining a realization */
void writefile(string filename, struct realization_mod *realin)
{
FILE *fp;
int i;
/* string filename;*/
/*struct realization_mod *realin;*/
/*save the permeability realization*/
fp = fopen(filename, "w");
/* printf("writefile.c : (*realin).vector[0] %f\n",(*realin).vector[0]); */
/* (*realin).vector[0]=1.1; */
for (i=0;i<(*realin).n;i++)
fprintf(fp,"%15.10f\n",(*realin).vector[i]);
fclose(fp);
return;
}

@ -1,33 +0,0 @@
#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include <math.h>
#include "genlib.h"
#include "geostat.h"
/*writefile */
/* write in the file "filename" the vector values of a */
/* realization_mod variable */
/* filename: explicit */
/* realin: structure defining a realization */
void writefile_bin(string filename, struct realization_mod *realin)
{
FILE *fp;
int i;
double *prout;
/*save the permeability realization*/
fp = fopen(filename, "w");
for (i=0;i<(*realin).n;i++)
{
prout=(*realin).vector;
fwrite(&prout[i],sizeof(double),1,fp);
}
fclose(fp);
return;
}

Some files were not shown because too many files have changed in this diff Show More

Loading…
Cancel
Save