/
Mgetcregather.c
147 lines (118 loc) · 4.3 KB
/
Mgetcregather.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
/* Build CRE gather given CMP X Offset CRE trajectory coordinates and interpolated data cube
This program searches for the closest trace to the CRE trajectory to build the CRE Gather for each (m, h) pair given in the interpolated data cube.
Programmer: Rodolfo A. C. Neves (Dirack) 04/09/2019
Email: rodolfo_profissional@hotmail.com
License: GPL-3.0 <https://www.gnu.org/licenses/gpl-3.0.txt>.
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <rsf.h>
int main(int argc, char* argv[])
{
float*** t; // Data cube A(m,h,t)
float*** creGather; // CRE gathers
float** m; // CMPs
bool verb; // Verbose parameters
float dm; // CMP sampling
float om; // CMP axis origin
int nm; // Number of CMP samples
float dh; // Offset sampling
float oh; // Offset axis origin
int nh; // Offset number of samples
float dt; // Time sampling
float ot; // Time axis origin
int nt; // Number of time samples
int i,j,k,l; // loop counter
int cn1; // Number of traces in CRE vector
float cd1; // CRE Gather sampling
float co1; // CRE Gather axis origin
int cn2; // Number of m0s x t0s pairs
int trac_m; // CMP sample index
int aperture; // Number of traces in CRE Gather
float mMax; // maximum CMP coordinate of the model
int nm0; // Number of m0s
int nt0; // Number of t0s
/* RSF files I/O */
sf_file in, out, out_m, cremh;
/* RSF files axis */
sf_axis ax,ay,az,am1,am2;
sf_init(argc,argv);
in = sf_input("in"); // Data cube A(m,h,t)
cremh = sf_input("cremh"); /* CRE Gather m(h) coordinates */
out = sf_output("out"); // CRE Gather
out_m = sf_output("m"); // CRE Gather CMP coordinate
/* seismic data cube A(m,h,t) */
if (!sf_histint(in,"n1",&nt)) sf_error("No n1= in input file");
if (!sf_histfloat(in,"d1",&dt)) sf_error("No d1= in input file");
if (!sf_histfloat(in,"o1",&ot)) sf_error("No o1= in input file");
if (!sf_histint(in,"n2",&nh)) sf_error("No n2= in input file");
if (!sf_histfloat(in,"d2",&dh)) sf_error("No d2= in input file");
if (!sf_histfloat(in,"o2",&oh)) sf_error("No o2= in input file");
if (!sf_histint(in,"n3",&nm)) sf_error("No n3= in input file");
if (!sf_histfloat(in,"d3",&dm)) sf_error("No d3= in input file");
if (!sf_histfloat(in,"o3",&om)) sf_error("No o3= in input file");
/* cre trajectories m(h) */
if(!sf_histint(cremh,"n1",&cn1)) sf_error("No n1= in cremh file");
if(!sf_histfloat(cremh,"d1",&cd1)) sf_error("No d1= in cremh file");
if(!sf_histfloat(cremh,"o1",&co1)) sf_error("No o1= in cremh file");
if(!sf_histint(cremh,"n2",&cn2)) sf_error("No n2= in cremh file");
if(!sf_getint("nm0",&nm0)) sf_error("Need nm0");
/* Number of central CMPs in cremh file */
if(!sf_getint("nt0",&nt0)) sf_error("Need nt0");
/* Number of t0s in cremh file */
if(!sf_getint("aperture",&aperture)) aperture=1;
/* Number of traces to put in a CRE Gather*/
if(aperture > cn1){
sf_error("The aperture can't be > n1 in cremh file\naperture=%i n2=%i",aperture,cn1);
}
if(! sf_getbool("verb",&verb)) verb=0;
/* 1: active mode; 0: quiet mode */
if (verb) {
sf_warning("Active mode on!!!");
sf_warning("CRE gather coordinates m(h) (cremh file): ");
sf_warning("n1=%d d1=%f o1=%f",cn1,cd1,co1);
sf_warning("n2=%d",cn2);
sf_warning("Input file dimensions: ");
sf_warning("n1=%d d1=%f o1=%f",nt,dt,ot);
sf_warning("n2=%d d2=%f o2=%f",nh,dh,oh);
sf_warning("n3=%d d3=%f o3=%f",nm,dm,om);
}
/* Read data cube */
t=sf_floatalloc3(nt,nh,nm);
sf_floatread(t[0][0],nh*nm*nt,in);
/* Read cre trajectories */
m = sf_floatalloc2(cn1,cn2);
sf_floatread(m[0],cn1*cn2,cremh);
creGather = sf_floatalloc3(nt,aperture,cn2);
mMax = om+dm*nm;
for(l=0;l<nm0;l++){
for(k=0;k<nt0;k++){
for(i=0;i<aperture;i++){
trac_m = (int)((double)m[(l*nt0)+k][i]/dm);
for(j=0;j<nt;j++){
creGather[(l*nt0)+k][i][j] =
(m[(l*nt0)+k][i] <= mMax)?
t[trac_m][i][j] : 0.;
}
}
}/* loop over t0s */
}/* loop over m0s */
/* eixo = sf_maxa(n,o,d)*/
ax = sf_maxa(nt, ot, dt);
ay = sf_maxa(aperture, oh, dh);
az = sf_maxa(cn2, 0, 1);
/* sf_oaxa(arquivo, eixo, índice do eixo) */
sf_oaxa(out,ax,1);
sf_oaxa(out,ay,2);
sf_oaxa(out,az,3);
sf_floatwrite(creGather[0][0],cn2*aperture*nt,out);
/* eixo do vetor m */
am1 = sf_maxa(cn1,oh,dh);
am2 = sf_maxa(1,0,1);
sf_oaxa(out_m,am1,1);
sf_oaxa(out_m,az,2);
sf_oaxa(out_m,am2,3);
sf_floatwrite(m[0],cn1*cn2,out_m);
exit(0);
}