From 6acb6d353a7ee110bda440a0be0ae0245eddf946 Mon Sep 17 00:00:00 2001
From: Jan Thorbecke <janth@xs4all.nl>
Date: Thu, 9 Aug 2018 16:35:28 +0200
Subject: [PATCH] added option to output frequency slices
---
extrap3d/lib/dataFileIO.c | 91 ++++++++++++++++++++++++-----------
extrap3d/main/cfpmod3d.c | 8 +--
extrap3d/main/extrap3d.c | 8 +--
extrap3d/main/zomodel3d.c | 11 +++--
fdelmodc/demo/modelOilGas.scr | 8 +--
5 files changed, 85 insertions(+), 41 deletions(-)
diff --git a/extrap3d/lib/dataFileIO.c b/extrap3d/lib/dataFileIO.c
index 84b910c..4612557 100644
--- a/extrap3d/lib/dataFileIO.c
+++ b/extrap3d/lib/dataFileIO.c
@@ -107,7 +107,7 @@ int read_FFT_DataFile(FILE *fp, complex *data, Area vel_area, int nfft, int nw,
}
-int write_FFT_DataFile(FILE *fp, complex *data, Area data_area, int fldr, int nt, int nfft, int nw, int nw_low, float dt, int out_su, int conjg, int verbose)
+int write_FFT_DataFile(FILE *fp, complex *data, Area data_area, int fldr, int nt, int nfft, int nw, int nw_low, float dt, int out_su, int conjg, int freq, int verbose)
{
int ix, iy, iw, i, nx, ny, sxy, nfreq, pos, sign;
int xvmin, yvmin;
@@ -142,37 +142,72 @@ int write_FFT_DataFile(FILE *fp, complex *data, Area data_area, int fldr, int nt
hdr[0].d2 = dx;
hdr[0].scalco = -1000;
- for (iy=0; iy<ny; iy++) {
- for (ix=0; ix<nx; ix++) {
- pos = iy*nx+ix;
- for (i=0; i<nfreq; i++)
- ctrace[i].r = ctrace[i].i = 0.0;
- for (iw=0; iw<nw; iw++) {
- ctrace[iw+nw_low].r = data[iw*sxy+pos].r;
- ctrace[iw+nw_low].i = conjg*data[iw*sxy+pos].i;
- }
- /* transform back to time */
- cr1fft(ctrace, trace, nfft, sign);
-
- /* write to output file */
- if (out_su) {
- hdr[0].gx = (int)(xvmin+ix*dx)*1000;
- hdr[0].gy = (int)(yvmin+iy*dy)*1000;
- hdr[0].f2 = xvmin+ix*dx;
- hdr[0].tracf = ix+1;
- hdr[0].tracl = iy*nx+ix+1;
- nwrite = fwrite(hdr, 1, TRCBYTES, fp);
- assert( nwrite == TRCBYTES );
- nwrite = fwrite(trace, sizeof(float), nt, fp);
- assert( nwrite == nt );
+ if (freq==0) {
+ for (iy=0; iy<ny; iy++) {
+ for (ix=0; ix<nx; ix++) {
+ pos = iy*nx+ix;
+ for (i=0; i<nfreq; i++)
+ ctrace[i].r = ctrace[i].i = 0.0;
+ for (iw=0; iw<nw; iw++) {
+ ctrace[iw+nw_low].r = data[iw*sxy+pos].r;
+ ctrace[iw+nw_low].i = conjg*data[iw*sxy+pos].i;
+ }
+
+ /* transform back to time */
+ cr1fft(ctrace, trace, nfft, sign);
+
+ /* write to output file */
+ if (out_su) {
+ hdr[0].gx = (int)(xvmin+ix*dx)*1000;
+ hdr[0].gy = (int)(yvmin+iy*dy)*1000;
+ hdr[0].f2 = xvmin+ix*dx;
+ hdr[0].tracf = ix+1;
+ hdr[0].tracl = iy*nx+ix+1;
+ nwrite = fwrite(hdr, 1, TRCBYTES, fp);
+ assert( nwrite == TRCBYTES );
+ nwrite = fwrite(trace, sizeof(float), nt, fp);
+ assert( nwrite == nt );
+ }
+ else {
+ nwrite = fwrite(trace, sizeof(float), nt, fp);
+ assert( nwrite == nt );
+ }
}
- else {
- nwrite = fwrite(trace, sizeof(float), nt, fp);
- assert( nwrite == nt );
+ }
+ }
+ else { /* Frequency output slices */
+ hdr[0].trid = 122;
+ hdr[0].fldr = fldr;
+ hdr[0].trwf = ny;
+ hdr[0].ntr = nw*ny;
+ hdr[0].ns = nx;
+ hdr[0].dt = 1000000*dt;
+ hdr[0].f1 = 0.0;
+ hdr[0].f2 = yvmin;
+
+ for (iw=0; iw<nw; iw++) {
+ for (iy=0; iy<ny; iy++) {
+ /* write to output file */
+ if (out_su) {
+ hdr[0].gx = (int)(xvmin)*1000;
+ hdr[0].gy = (int)(yvmin+iy*dy)*1000;
+ hdr[0].f1 = xvmin+ix*dx;
+ hdr[0].f2 = yvmin+iy*dy;
+ hdr[0].tracf = iy+1;
+ hdr[0].tracl = iw*ny+iy+1;
+ nwrite = fwrite(hdr, 1, TRCBYTES, fp);
+ assert( nwrite == TRCBYTES );
+ nwrite = fwrite(&data[iw*sxy+iy*nx], sizeof(complex), nx, fp);
+ assert( nwrite == nx );
+ }
+ else {
+ nwrite = fwrite(&data[iw*sxy+iy*nx], sizeof(complex), nx, fp);
+ assert( nwrite == nx );
+ }
}
}
- }
+ }
fflush(fp);
free(trace);
diff --git a/extrap3d/main/cfpmod3d.c b/extrap3d/main/cfpmod3d.c
index 6882e55..2eb38af 100644
--- a/extrap3d/main/cfpmod3d.c
+++ b/extrap3d/main/cfpmod3d.c
@@ -24,7 +24,7 @@ void cr1fft(complex *cdata, float *rdata, int n, int sign);
/****** IO routines *******/
int openVelocityFile(char *file_vel, FILE **fp, Area *vel_area, int verbose);
void readVelocitySlice(FILE *fp, float *velocity, int iz, int nyv, int nxv);
-int write_FFT_DataFile(FILE *fp, complex *data, Area data_area, int fldr, int nt, int nfft, int nw, int nw_low, float dt, int out_su, int conjg, int verbose);
+int write_FFT_DataFile(FILE *fp, complex *data, Area data_area, int fldr, int nt, int nfft, int nw, int nw_low, float dt, int out_su, int conjg, int freq, int verbose);
int write_ImageFile(FILE *fp, float *data, Area data_area, int fldr, int d, int out_su, int verbose);
@@ -129,6 +129,7 @@ char *sdoc[] = {
" filter_inc=1 ............. the increment in dz for the Li filter (0=off)",
" nterms=1 ................. the number of terms in the paraxial expansion",
" OUTPUT DEFINITION ",
+" freq=0 ................... output in frequency slices",
" snap=0 ................... snapshots (not yet implemented)",
" beam=0 ................... beams ",
" verbose=0 ................ >1: shows various parameters and results",
@@ -152,7 +153,7 @@ int main(int argc, char *argv[])
{
FILE *out_file, *vel_file, *src_file, *beam_file;
size_t nread, size, size_out;
- int verbose, method, ntraces, oper_opt, MB, out_su, verb_root;
+ int verbose, method, ntraces, oper_opt, MB, out_su, out_freq, verb_root;
int nxv, nyv, nzv, dstep, fldr;
int nt, err;
int ntap, tap_opt, order, McC, oplx, oply, fine;
@@ -239,6 +240,7 @@ int main(int argc, char *argv[])
if(!getparint("fine", &fine)) fine = 2;
if(!getparint("beam", &beam)) beam = 0;
if(!getparint("add", &add)) add = 0;
+ if(!getparint("freq", &out_freq)) out_freq = 0;
if(!getparint("verbose", &verbose)) verbose = 0;
if(!ISODD(oplx)) oplx += 1;
@@ -651,7 +653,7 @@ int main(int argc, char *argv[])
fldr = is+1;
write_FFT_DataFile(out_file, rec_all, shot_area, fldr,
- nt, nfft, nw, nw_low, dt, out_su, conjg, verbose);
+ nt, nfft, nw, nw_low, dt, out_su, conjg, out_freq, verbose);
fclose(out_file);
if (verbose)
diff --git a/extrap3d/main/extrap3d.c b/extrap3d/main/extrap3d.c
index 17f6dab..6a4d6d3 100644
--- a/extrap3d/main/extrap3d.c
+++ b/extrap3d/main/extrap3d.c
@@ -22,7 +22,7 @@ int openVelocityFile(char *file_vel, FILE **fp, Area *vel_area, int verbose);
void readVelocitySlice(FILE *fp, float *velocity, int iz, int nyv, int nxv);
-int write_FFT_DataFile(FILE *fp, complex *data, Area data_area, int fldr, int nt, int nfft, int nw, int nw_low, float dt, int out_su, int conjg, int verbose);
+int write_FFT_DataFile(FILE *fp, complex *data, Area data_area, int fldr, int nt, int nfft, int nw, int nw_low, float dt, int out_su, int conjg, int freq, int verbose);
int read_FFT_DataFile(FILE *fp, complex *data, Area vel_area, int nfft, int nw, int nw_low, int *tr_read_in, int *tr_shot, int *ixmin, int *ixmax, int *iymin, int *iymax, int *sx, int *sy, int conjg, int verbose);
@@ -107,6 +107,7 @@ char *sdoc[] = {
" filter_inc=1 ............. the increment in dz for the Li filter (0=off)",
" nterms=1 ................. the number of terms in the paraxial expansion",
" OUTPUT DEFINITION ",
+" freq=0 ................... output in frequency slices",
" snap=0 ................... snapshots (not yet implemented)",
" beam=0 ................... beams ",
" verbose=0 ................ >1: shows various parameters and results",
@@ -136,7 +137,7 @@ int main(int argc, char *argv[])
size_t nread, bytes, size, trace_sz, size_out;
int verbose, method, ntraces, verb_root;
int nxv, nyv, nzv, binary_file, dstep, id, id1, id2;
- int d, nt, ndepth, i, j, conjg, conjgs, mode, out_su;
+ int d, nt, ndepth, i, j, conjg, conjgs, mode, out_su, out_freq;
int ntap, tap_opt, order, McC, oplx, oply, fine, MB;
int stackmigr, imc, area, ixmin, ixmax, iymin, iymax, ns;
int nfft, nfreq, nw_high, nw_low, nw, sx, sy, ix, iy;
@@ -206,6 +207,7 @@ int main(int argc, char *argv[])
if(!getparfloat("weights", &weights)) weights = 1e-2;
if(!getparint("fine", &fine)) fine = 2;
if(!getparint("beam", &beam)) beam = 0;
+ if(!getparint("freq", &out_freq)) out_freq = 0;
if(!getparint("verbose", &verbose)) verbose = 0;
if(!ISODD(oplx)) oplx += 1;
@@ -602,7 +604,7 @@ int main(int argc, char *argv[])
assert( out_file );
write_FFT_DataFile(out_file, rec_all, shot_area, (is+1),
- nt, nfft, nw, nw_low, dt, out_su, conjg, verb_root);
+ nt, nfft, nw, nw_low, dt, out_su, conjg, out_freq, verb_root);
fclose(out_file);
}
diff --git a/extrap3d/main/zomodel3d.c b/extrap3d/main/zomodel3d.c
index 39f66a2..1bb9a53 100644
--- a/extrap3d/main/zomodel3d.c
+++ b/extrap3d/main/zomodel3d.c
@@ -24,7 +24,7 @@ void cr1fft(complex *cdata, float *rdata, int n, int sign);
/****** IO routines *******/
int openVelocityFile(char *file_vel, FILE **fp, Area *vel_area, int verbose);
void readVelocitySlice(FILE *fp, float *velocity, int iz, int nyv, int nxv);
-int write_FFT_DataFile(FILE *fp, complex *data, Area data_area, int fldr, int nt, int nfft, int nw, int nw_low, float dt, int out_su, int conjg, int verbose);
+int write_FFT_DataFile(FILE *fp, complex *data, Area data_area, int fldr, int nt, int nfft, int nw, int nw_low, float dt, int out_su, int conjg, int freq, int verbose);
/***** Operator tables *****/
void tablecalc_2D(int oplx, int oply, int nx, float dx, int ny, float dy,
@@ -64,7 +64,6 @@ char *sdoc[] = {
" Optional parameters:",
" ",
" conjg=0 .................. 1: take complex conjugate of input data",
-" verbose=0 ................ >1: shows various parameters and results",
" VELOCITY MODEL ",
" dxv=dxv .................. stepsize in x-direction of velocity model ",
" dyv=dyv .................. stepsize in y-direction of velocity model ",
@@ -113,6 +112,9 @@ char *sdoc[] = {
" fine=2 ................... fine sampling in operator table",
" filter_inc=1 ............. the increment in dz for the Li filter (0=off)",
" nterms=1 ................. the number of terms in the paraxial expansion",
+" OUTPUT DEFINITION ",
+" freq=0 ................... output in frequency slices",
+" verbose=0 ................ >1: shows various parameters and results",
" ",
" Options for method:",
" - 1 = Direct 2D convolution (Walter's scheme)",
@@ -148,7 +150,7 @@ int main(int argc, char *argv[])
int nfft, nfreq, nw_high, nw_low, nw, ix, iy;
int sxy, iw, sign, ixv, iyv;
int nxy, pos, id, idp, oper_opt;
- int out_su, nterms, filter_inc, interpol;
+ int out_freq, out_su, nterms, filter_inc, interpol;
Area shot_area;
float alpha, weight;
float depth1,depth2;
@@ -217,6 +219,7 @@ int main(int argc, char *argv[])
if(!getparint("fine", &fine)) fine = 2;
if(!getparfloat("dt", &dt)) dt = 0.004;
if(!getparint("nt", &nt)) nt = 512;
+ if(!getparint("freq", &out_freq)) out_freq = 0;
if(!getparint("verbose", &verbose)) verbose = 0;
if(!ISODD(oplx)) oplx += 1;
@@ -699,7 +702,7 @@ int main(int argc, char *argv[])
out_file = fopen( file_out, "w+" ); assert( out_file );
write_FFT_DataFile(out_file, rec_field, shot_area, 1,
- nt, nfft, nw, nw_low, dt, out_su, conjg, verbose);
+ nt, nfft, nw, nw_low, dt, out_su, conjg, out_freq, verbose);
fclose(out_file);
}
diff --git a/fdelmodc/demo/modelOilGas.scr b/fdelmodc/demo/modelOilGas.scr
index 5005543..5fa2690 100755
--- a/fdelmodc/demo/modelOilGas.scr
+++ b/fdelmodc/demo/modelOilGas.scr
@@ -1,6 +1,6 @@
#!/bin/bash
-#export PATH=$HOME/bin:$HOME/src/OpenSource/utils:$PATH:
+export PATH=$HOME/bin:$HOME/src/OpenSource/utils:$PATH:
cp=2000
rho=2500
@@ -16,12 +16,14 @@ makemod sizex=5000 sizez=2500 dx=$dx dz=$dx cp0=$cp ro0=$rho \
makewave w=g1 fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
-export OMP_NUM_THREADS=8
+export OMP_NUM_THREADS=1
zsrc=1100
zsrc=0
-~/bin/fdelmodc \
+which fdelmodc
+
+fdelmodc \
file_cp=syncl_cp.su ischeme=1 iorder=4 \
file_den=syncl_ro.su \
file_src=wave.su \
--
GitLab