Compare commits

..

36 Commits

Author SHA1 Message Date
Gervasio Daniel Pérez 4397259bf1 Usar Process para correr realization 1 year ago
Gervasio Daniel Pérez 74aae9f91d Actualizar 'README-TUPAC.md' 2 years ago
Gervasio Daniel Pérez cbb0768707 Update README TUPAC 2 years ago
Gervasio Daniel Pérez 7f0f57b58f Actualizar 'README-TUPAC.md' 2 years ago
Gervasio Daniel Pérez 949b7853ab Readme + plantillas SLURM
+ Makefile con posibilidad de correr en slurm
2 years ago
Gervasio Daniel Pérez fc50b6b600 Actualizar 'README.md' 2 years ago
Gervasio Daniel Pérez debbcaca50 README y Makefile actualizados 2 years ago
Gervasio Daniel Pérez f775ea5faf Slurm example 2 years ago
Gervasio Daniel Pérez 0aac969f23 Extras para NIX 2 years ago
Sebastián Santisi f2e15adaba Deleting Conda env 2 years ago
Sebastián Santisi 0dcc4054bb Adding benchmarker as a requirement 2 years ago
Sebastián Santisi bd43d201f9 Adding venv to run_simulation.sh 2 years ago
Sebastián Santisi ca1a550bd3 Prerequisites 2 years ago
Sebastián Santisi e6a68668f5 chmod to script_install.sh 2 years ago
Sebastián Santisi f323721159 Upgrading README 2 years ago
Sebastián Santisi c641704839 Install only via venv 2 years ago
Oli d1af191f70 explain branches 3 years ago
Oli a541475bc0 explain branches 3 years ago
chortas 84c6356613 Fix numpy 3 years ago
chortas 621f1509b9 Change install rule 3 years ago
Cecilia Hortas 96db89d10c
Merge pull request #3 from medios-porosos-fiuba/milestone_4
Milestone 4 - Part 1
3 years ago
chortas fffbb8f24d Delete travis file (it requires card now) 3 years ago
chortas 7c40bde80c Add conditional decorator 3 years ago
chortas ab87cc8099 Remove extra prints and change generation 3 years ago
chortas 17cb088c42 Fixing conflicts 3 years ago
chortas 03c9901e88 Fixing conflicts 3 years ago
Oli 96dd67da96 Merge branch 'milestone_4' of github.com:medios-porosos-fiuba/simulacion-permeabilidad into milestone_4 3 years ago
Oli df7e1f631e add memory profiling in generation 3 years ago
chortas bea64b75da Add connectivity perf test 3 years ago
chortas 6524fae7a3 Add travis file 3 years ago
Oli 8f1d1f3df3 Merge branch 'milestone_4' of github.com:medios-porosos-fiuba/simulacion-permeabilidad into milestone_4 3 years ago
Oli d3716caa23 fix bug. not nice 3 years ago
chortas 928c56ce87 Update README 3 years ago
Oli 1f613c87a8 change fttma invocation 3 years ago
Oli 10cb61d6a7 generation performance tests. wip 3 years ago
Cecilia Hortas 916decaeca
Milestone 2 and 3 #2 3 years ago

15
.gitignore vendored

@ -2,5 +2,18 @@ petsc/
petsc4py/
openmpi-2.1.1/
openmpi-2.1.1.tar.bz2
tools/connec/__pycache__/
tools/postprocessK/__pycache__/
tools/postprocessK/kperm/__pycache__/
tools/solver/__pycache__/
tools/__pycache__/
tools/generation/__pycache__/
output/
tests/integration/tmp_output/
output/
tests/performance/tmp_gen_output/
output/
tests/integration/__pycache__/
fftma_module/gen/build/
tools/connec/conec2d
tools/connec/conec3d
utilities/__pycache__/

@ -1,15 +1,23 @@
SHELL:= zsh
env:
source ./env.sh
install:
source ./script_install.sh $(SHELL)
bash ./script_install.sh
run:
bash ./run_simulation.sh
test:
cd tests/integration && python3 -m unittest test.py
perf:
cd tests/performance && python3 generation.py
cd tests/performance && python3 connectivity.py
run: env
./run_simulation.sh
run-slurm:
sbatch run-simulation-slurm
binaries: ./script_fortran.sh
test-slurm:
sbatch run-test-slurm
test: env binaries
cd tests/integration && python test.py
perf-slurm:
sbatch run-perf-slurm

@ -0,0 +1,73 @@
# Paralelización de simulación de permeabilidad
## Clonar este repositorio
```
# Por HTTPS
git clone https://git.csc.gob.ar/ssantisi/simulacion-permeabilidad.git
# Por SSH
git clone git@git.csc.gob.ar:ssantisi/simulacion-permeabilidad.git
```
Entrar a la carpeta del proyecto
```
cd simulacion-permeabilidad
```
## Instalación de las librerías (Python 3)
Como requisitos previos a la instalación hace falta tener:
* `python3`
* `python3-pip`
* `gfortran`
* `gcc`
* `build-essentials`
```
make install
```
El instalador genera un entorno virtual en `venv`.
Para activar el entorno virtual:
```
source venv/bin/activate
```
Para desactivarlo `deactivate`.
## Correr la simulación
```
make run
```
## Correr los casos de prueba
```
make test
```
## Correr las pruebas de performance
```
make perf
```
## Branches de github
Explicacion sobre las branches que quedaron en github para entender cuales usar:
- `main` -> solo la migracion a python3 sin ninguna cambio sobre el modulo de FFTMA (sin mejoras)
- `improvement_NOMBRE` -> branch con el `NOMBRE` de la mejora que esta explicado que hace cada una en el analisis
donde la primera es `remove_generate_array` y la ultima es `multiple_buffers`
- `improvement_NOMBRE-logs` -> misma branch de mejora pero agregando informacion de logs que se utilizo para realizar los analisis y no deberian utilizarse porque degradan la performance considerablemente.
- `migrate_fortran` -> una branch donde se trato de migrar de fortran a C pero quedo inconcluso el trabajo y no se continuo.

@ -0,0 +1,126 @@
# Paralelización de simulación de permeabilidad en Cluster HPC (SLURM)
## ANTES DEL PRIMER USO
Seguir los siguientes pasos del tutorial:
1. Clonar este repositorio
2. Entrar al entorno NIX
3. Generar entorno Python para la simulación
## USO HABITUAL:
Habiendo completado las instrucciones previas al primer uso:
1. Entrar al entorno NIX
2. Activar el entorno python
3. Ejecutar las siimulaciones de interés
## Clonar este repositorio
```
# Por HTTPS
git clone https://git.csc.gob.ar/ssantisi/simulacion-permeabilidad.git
# Por SSH
git clone git@git.csc.gob.ar:ssantisi/simulacion-permeabilidad.git
```
Entrar a la carpeta del proyecto
```
cd simulacion-permeabilidad
```
## Entrar al entorno NIX
```
nix-shell
```
Para más info ver el archivo `shell.nix`.
## Generar entorno Python para la simulación
Los siguentes requisitos ya vienen provistos por `nix-shell`:
* `python3`
* `python3-pip`
* `gfortran`
* `gcc`
* `build-essentials`
```
make install
```
El instalador genera un entorno virtual en `venv`.
## Activar el entorno python
Para activar el entorno virtual (LUEGO DE HABER ACTIVADO EL ENTORNO NIX):
```
source venv/bin/activate
```
Para desactivarlo `deactivate`.
## Correr simulaciones de interés
Los siguientes comandos encolan trabajos en el cluster, que correrán tan pronto como les toque su lugar en la cola de trabajos.
Se puede monitorear el estado de el/los trabajos encolados con `squeue`.
Se puede cancelar un trabajo encolado con `scancel numero_de_trabajo`
### Correr la simulación (en SLURM)
Por defecto: 2 nodos, 32 tasks por nodo.
Para cambiarlo editar archivo `run-simulation-slurm`.
```
make run-slurm
```
Una vez en ejecución, la salida se escribe a los siguientes archivos:
* `simulation-nro-de-trabajo.out` la salida del trabajo.
* `simulation-nro-de-trabajo.err` los errores del proceso
### Correr los casos de prueba (en SLURM)
Por defecto: 1 nodo, 64 cpus por nodo.
Para cambiarlo editar archivo `run-tests-slurm`.
```
make test-slurm
```
Una vez en ejecución, la salida se escribe a los siguientes archivos:
* `test-nro-de-trabajo.out` la salida del trabajo.
* `test-nro-de-trabajo.err` los errores del proceso
### Correr las pruebas de performance (en SLURM)
Por defecto: 1 nodo, 64 cpus por nodo.
Para cambiarlo editar archivo `run-perf-slurm`.
```
make perf-slurm
```
Una vez en ejecución, la salida se escribe a los siguientes archivos:
* `perf-nro-de-trabajo.out` la salida del trabajo.
* `perf-nro-de-trabajo.err` los errores del proceso
## Branches de github
Explicacion sobre las branches que quedaron en github para entender cuales usar:
- `main` -> solo la migracion a python3 sin ninguna cambio sobre el modulo de FFTMA (sin mejoras)
- `improvement_NOMBRE` -> branch con el `NOMBRE` de la mejora que esta explicado que hace cada una en el analisis
donde la primera es `remove_generate_array` y la ultima es `multiple_buffers`
- `improvement_NOMBRE-logs` -> misma branch de mejora pero agregando informacion de logs que se utilizo para realizar los analisis y no deberian utilizarse porque degradan la performance considerablemente.
- `migrate_fortran` -> una branch donde se trato de migrar de fortran a C pero quedo inconcluso el trabajo y no se continuo.

@ -1 +1,9 @@
# simulacion-permeabilidad
# Paralelización de simulación de permeabilidad
## Corrida en PC de escritorio/laptop
Ver archivo **README-LOCAL.md**.
## Corrida en cluster TUPAC
Ver archivo **README-TUPAC.md**.

@ -1,5 +0,0 @@
#!/bin/bash
source ~/miniconda3/etc/profile.d/conda.sh
conda activate py2_7
python --version

@ -0,0 +1,47 @@
#include <stdlib.h>
#define SWITCH(a,b,temp) temp=b; b=a; a=temp;
int med5(double* A, long k, int* pj)
{
return 0;
}
int mom(double* A, long N, long i)
{
long k=0, q, rest;
int j;
double* Anext;
rest=N%5;
Anext=(double*)malloc(sizeof(double)*(q+1?0));///???
q=N-rest;
while(k<q)
{
med5(A,k,&j)
k+=5;
}
return 0;
}

@ -0,0 +1,39 @@
import numpy as np
A = [1,2,3,4,5,4.000000001,1000,8,18,9,991473,74168,6387,3678,8796,1343]
def median_of_medians(A, i):
#divide A into sublists of len 5
sublists = [A[j:j+5] for j in range(0, len(A), 5)]
medians = [sorted(sublist)[len(sublist)/2] for sublist in sublists]
if len(medians) <= 5:
pivot = sorted(medians)[len(medians)/2]
else:
#the pivot is the median of the medians
pivot = median_of_medians(medians, len(medians)/2)
#partitioning step
low = [j for j in A if j < pivot]
high = [j for j in A if j > pivot]
k = len(low)
if i < k:
return median_of_medians(low,i)
elif i > k:
return median_of_medians(high,i-k-1)
else: #pivot = k
return pivot
#Here are some example lists you can use to see how the algorithm works
#A = [1,2,3,4,5,1000,8,9,99]
#B = [1,2,3,4,5,6]
#print median_of_medians(A, 0) #should be 1
#print median_of_medians(A,7) #should be 99
#print median_of_medians(B,4) #should be 5
b=np.sort(A)
c=np.array([median_of_medians(A, i) for i in range(len(A))])
print b
print c
print bool(np.sum(b==c))

@ -0,0 +1,31 @@
import numpy as np
import random as rd
MAX=100000
def med3(l):
if l[0]<l[1]:
if l[1]<l[2]:
med=l[1]
else:
if l[0]<l[2]:
med=l[2]
else:
med=l[0]
else:
if l[0]<l[2]:
med=l[0]
else:
if l[1]<l[2]:
med=l[2]
else:
med=l[1]
return med
S=0
for i in range(1000000):
l=[float(rd.randint(0,MAX))/MAX,float(rd.randint(0,MAX))/MAX,float(rd.randint(0,MAX))/MAX]
#print l
#L=np.sort(np.array(l))
#print L
S+=np.sort(np.array(l))[1]==med3(l)
print S

@ -0,0 +1,36 @@
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))

@ -0,0 +1,442 @@
#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_

@ -0,0 +1,163 @@
/*
* 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

@ -0,0 +1,690 @@
#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

@ -0,0 +1,21 @@
#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

@ -0,0 +1,75 @@
/*
* 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

@ -0,0 +1,80 @@
#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

@ -0,0 +1,91 @@
#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

@ -0,0 +1,13 @@
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

@ -0,0 +1,21 @@
# 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

@ -0,0 +1,63 @@
/*
* 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);
}

@ -0,0 +1,163 @@
/*
* 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

@ -0,0 +1,77 @@
/*
* 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);
}

@ -0,0 +1,73 @@
/*
* 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

@ -0,0 +1,162 @@
/*
* 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);
}

@ -0,0 +1,159 @@
/*
* 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

@ -0,0 +1,157 @@
/*
* 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);
}

@ -0,0 +1,75 @@
/*
* 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

@ -0,0 +1,122 @@
/*
* 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;
}

@ -0,0 +1,115 @@
/*
* 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

@ -0,0 +1,246 @@
/*
* 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));
}

@ -0,0 +1,243 @@
/*
* 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

@ -0,0 +1,180 @@
/*
* 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);
}

@ -0,0 +1,85 @@
/*
* 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

@ -0,0 +1,30 @@
########################################################################
#
# 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

@ -0,0 +1,106 @@
#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;
}

@ -0,0 +1,41 @@
#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;
}

@ -0,0 +1,18 @@
#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);
}

@ -0,0 +1,43 @@
#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;
}

@ -0,0 +1,55 @@
#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);
}

@ -0,0 +1,90 @@
#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;
}

@ -0,0 +1,17 @@
#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);
}

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

@ -0,0 +1,185 @@
#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;
}

@ -0,0 +1,591 @@
#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;
}

@ -0,0 +1,16 @@
#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);
}

@ -0,0 +1,31 @@
#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);
}
}

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

@ -0,0 +1,43 @@
#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;
}

@ -0,0 +1,689 @@
#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

@ -0,0 +1,39 @@
#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);
}

@ -0,0 +1,36 @@
#
########################################################################
#
# 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

@ -0,0 +1,35 @@
#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);
}

@ -0,0 +1,78 @@
#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;
}

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

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

@ -0,0 +1,50 @@
#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);
}

@ -0,0 +1,17 @@
#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);
}

@ -0,0 +1,10 @@
#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)));
}

@ -0,0 +1,26 @@
#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);
}

@ -0,0 +1,110 @@
#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;
}

@ -0,0 +1,30 @@
########################################################################
#
# 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

@ -0,0 +1,119 @@
#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;
}

@ -0,0 +1,45 @@
#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;
}

@ -0,0 +1,43 @@
#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;
}

@ -0,0 +1,105 @@
#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;
}

@ -0,0 +1,79 @@
#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;
}

@ -0,0 +1,76 @@
#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;
}

@ -0,0 +1,54 @@
#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;
}

@ -0,0 +1,82 @@
#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

@ -0,0 +1,162 @@
#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;
}

@ -0,0 +1,157 @@
#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;
}

@ -0,0 +1,68 @@
#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;
}

@ -0,0 +1,34 @@
########################################################################
#
# 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

@ -0,0 +1,66 @@
#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;
}

@ -0,0 +1,117 @@
#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;
}

@ -0,0 +1,50 @@
#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;
}

@ -0,0 +1,216 @@
#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;
}
}

@ -0,0 +1,33 @@
/*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;
}

@ -0,0 +1,48 @@
#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;
}

@ -0,0 +1,99 @@
#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;
}

@ -0,0 +1,59 @@
#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

@ -0,0 +1,35 @@
#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;
}

@ -0,0 +1,50 @@
#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;
}

@ -0,0 +1,156 @@
//#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);
}*/

@ -0,0 +1,31 @@
#
########################################################################
#
# 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

@ -0,0 +1,16 @@
#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;
}

@ -0,0 +1,73 @@
#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;
}

@ -0,0 +1,122 @@
#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;
}

@ -0,0 +1,71 @@
#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;
}

@ -0,0 +1,21 @@
#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

@ -0,0 +1,204 @@
#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;
}

@ -0,0 +1,205 @@
#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;
}

@ -0,0 +1,206 @@
#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;
}

@ -0,0 +1,35 @@
#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;
}

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

Loading…
Cancel
Save