Skip to content

Commit 2bf2625

Browse files
sfomelsfomel
committed
convection
1 parent 60a5a5a commit 2bf2625

File tree

4 files changed

+679
-17
lines changed

4 files changed

+679
-17
lines changed

‎user/fomels/Mcnvpaint.c‎

Lines changed: 116 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,116 @@
1+
/* Conection painting. */
2+
/*
3+
Copyright (C) 2004 University of Texas at Austin
4+
5+
This program is free software; you can redistribute it and/or modify
6+
it under the terms of the GNU General Public License as published by
7+
the Free Software Foundation; either version 2 of the License, or
8+
(at your option) any later version.
9+
10+
This program is distributed in the hope that it will be useful,
11+
but WITHOUT ANY WARRANTY; without even the implied warranty of
12+
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13+
GNU General Public License for more details.
14+
15+
You should have received a copy of the GNU General Public License
16+
along with this program; if not, write to the Free Software
17+
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18+
*/
19+
#include <math.h>
20+
21+
#include <rsf.h>
22+
#include "predict.h"
23+
24+
int main (int argc, char* argv[])
25+
{
26+
int order, n1, n2, n12, i1, i2, i0, nw;
27+
char *label, *unit;
28+
float **u, ***conv, *trace;
29+
float o1, d1, o2, d2, eps, *time;
30+
sf_file cnv, out, seed;
31+
32+
sf_init(argc, argv);
33+
cnv = sf_input("in");
34+
out = sf_output("out");
35+
36+
if (SF_FLOAT != sf_gettype(cnv)) sf_error("Need float input");
37+
if (!sf_histint(cnv,"n1",&n1)) sf_error("Need n1= in input");
38+
39+
if (!sf_histint(cnv,"n3",&n2)) sf_error("Need n3= in input");
40+
if (!sf_histfloat(cnv,"o3",&o2)) o2=0.;
41+
if (!sf_histfloat(cnv,"d3",&d2)) d2=1.;
42+
43+
n12 = n1*n2;
44+
45+
if (!sf_histint(cnv,"n2",&nw)) sf_error("Need n2= in input");
46+
47+
sf_putint(out,"n2",n2);
48+
sf_putfloat(out,"d2",d2);
49+
sf_putfloat(out,"o2",o2);
50+
if (NULL != (label= sf_histstring(cnv,"label3"))) sf_putstring(out,"label2",label);
51+
if (NULL != (unit= sf_histstring(cnv,"unit3"))) sf_putstring(out,"unit2",unit);
52+
53+
sf_putint(out,"n3",1);
54+
55+
order = nw/2;
56+
57+
if (NULL != sf_getstring("seed")) {
58+
seed = sf_input("seed");
59+
time = NULL;
60+
} else {
61+
seed = NULL;
62+
time = sf_floatalloc(n1);
63+
if (!sf_histfloat(cnv,"o1",&o1)) o1=0.;
64+
if (!sf_histfloat(cnv,"d1",&d1)) d1=1.;
65+
for (i1=0; i1 < n1; i1++) {
66+
time[i1] = o1+i1*d1;
67+
}
68+
}
69+
70+
if (!sf_getfloat("eps",&eps)) eps=0.01;
71+
/* regularization */
72+
73+
if (!sf_getint("i0",&i0)) i0=0;
74+
/* reference trace */
75+
76+
predict_init (n1, n2, eps*eps, order, 1, false);
77+
78+
u = sf_floatalloc2(n1,n2);
79+
conv = sf_floatalloc3(n1,nw,n2);
80+
trace = sf_floatalloc(n1);
81+
82+
sf_floatread(conv[0][0],n12*nw,cnv);
83+
84+
if (NULL != seed) {
85+
sf_floatread(trace,n1,seed);
86+
} else {
87+
for (i1=0; i1 < n1; i1++) {
88+
trace[i1] = time[i1];
89+
}
90+
}
91+
92+
for (i1=0; i1 < n1; i1++) {
93+
u[i0][i1] = trace[i1];
94+
}
95+
for (i2=i0-1; i2 >= 0; i2--) {
96+
predict_step(false,false,trace,conv[i2]);
97+
for (i1=0; i1 < n1; i1++) {
98+
u[i2][i1] = trace[i1];
99+
}
100+
}
101+
for (i1=0; i1 < n1; i1++) {
102+
trace[i1] = u[i0][i1];
103+
}
104+
for (i2=i0+1; i2 < n2; i2++) {
105+
predict_step(false,true,trace,conv[i2-1]);
106+
for (i1=0; i1 < n1; i1++) {
107+
u[i2][i1] = trace[i1];
108+
}
109+
}
110+
sf_floatwrite(u[0],n1*n2,out);
111+
112+
exit(0);
113+
}
114+
115+
116+

‎user/fomels/SConstruct‎

Lines changed: 17 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -6,23 +6,23 @@ progs = '''
66
abalance analytical angle angle2 approx arrival bdix beamform1 bil1
77
bilat2 blur boxcascade cbeamform1 causinv ccausint cchebyshevp
88
cfftexpa-dev cfftwave1 cdivn cflow chain2dfft chebvc chebyshev
9-
chebyshevp clpf cltft cltftfft cnvd constperm constpermh constpermh1
10-
convection cosftwave1 cpef cr ctf2dprec deblur diphase distance divn
11-
dix donut dpeiko eikonal eikonalvti eno2 erf exgr fedchain fedchain1
12-
fedchain2 fedchain21 fft2 fftexp0 fftexp1 fftexp3 fftexpa fftexp0a
13-
fftone ffttest fftwave1 fftwave2 fftwave3 findmin2 focus fpow freqest
14-
gaussmooth gbeamform imray interf interp2 interpt iphase kdsort kdtree
15-
kron label legacy lfftexp0 llpf localskew locov lpf lrmig0 lsfit max2
16-
median mffit mig3 miss3 morph nconv nnint nnshape nnshapet nsmooth
17-
nsmooth1 ocparcel octentwt ofilp ofsemb ortho patch pchain pchain1
18-
phaserot pick pick2 pick3 plane poly polyfit psmig2 regr rdiv rect1
19-
reshape riesz rsin sc seislet1 semblance semblancew shape shapeagc
20-
shapefill shearer shift1 sh1 sh2 similarity similarity2 simenv
21-
slschain2d smoothderw smoothreg smspray stfchain stfchain2 stcontrib
22-
stltft stpf stphase stream streamiss sttimefreq taupfit tfchain
23-
tf2dprec thin timecont timefreq tomo tree trace2 tricascade tristack
24-
tristack2 twofreq2 upgrad var2 velcon velinv vipmig0 vofzperm warp1
25-
warpadd warpscan zero zmarch ztrace
9+
chebyshevp clpf cltft cltftfft cnvd cnvpaint constperm constpermh
10+
constpermh1 convection cosftwave1 cpef cr ctf2dprec deblur diphase
11+
distance divn dix donut dpeiko eikonal eikonalvti eno2 erf exgr
12+
fedchain fedchain1 fedchain2 fedchain21 fft2 fftexp0 fftexp1 fftexp3
13+
fftexpa fftexp0a fftone ffttest fftwave1 fftwave2 fftwave3 findmin2
14+
focus fpow freqest gaussmooth gbeamform imray interf interp2 interpt
15+
iphase kdsort kdtree kron label legacy lfftexp0 llpf localskew locov
16+
lpf lrmig0 lsfit max2 median mffit mig3 miss3 morph nconv nnint
17+
nnshape nnshapet nsmooth nsmooth1 ocparcel octentwt ofilp ofsemb ortho
18+
patch pchain pchain1 phaserot pick pick2 pick3 plane poly polyfit
19+
psmig2 regr rdiv rect1 reshape riesz rsin sc seislet1 semblance
20+
semblancew shape shapeagc shapefill shearer shift1 sh1 sh2 similarity
21+
similarity2 simenv slschain2d smoothderw smoothreg smspray stfchain
22+
stfchain2 stcontrib stltft stpf stphase stream streamiss sttimefreq
23+
taupfit tfchain tf2dprec thin timecont timefreq tomo tree trace2
24+
tricascade tristack tristack2 twofreq2 upgrad var2 velcon velinv
25+
vipmig0 vofzperm warp1 warpadd warpscan zero zmarch ztrace
2626
'''
2727

2828
ccprogs = 'anisolr2 isolr2 isolr3 lrvc0 permlr1 permlr2 permlr3'

‎user/fomels/cnv.c‎

Lines changed: 164 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,164 @@
1+
/* Matrix definition for convection */
2+
/*
3+
Copyright (C) 2025 University of Texas at Austin
4+
5+
This program is free software; you can redistribute it and/or modify
6+
it under the terms of the GNU General Public License as published by
7+
the Free Software Foundation; either version 2 of the License, or
8+
(at your option) any later version.
9+
10+
This program is distributed in the hope that it will be useful,
11+
but WITHOUT ANY WARRANTY; without even the implied warranty of
12+
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13+
GNU General Public License for more details.
14+
15+
You should have received a copy of the GNU General Public License
16+
along with this program; if not, write to the Free Software
17+
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18+
*/
19+
20+
#include <rsf.h>
21+
/*^*/
22+
23+
#include "cnv.h"
24+
25+
#ifndef _cnv_h
26+
27+
typedef struct Cnv *cnv; /* abstract data type */
28+
/*^*/
29+
30+
#endif
31+
32+
struct Cnv {
33+
int n, na;
34+
float **a, *b;
35+
};
36+
37+
cnv cnv_init(int n1 /* trace length */,
38+
int nw /* filter order */)
39+
/*< initialize >*/
40+
{
41+
cnv w;
42+
43+
w = (cnv) sf_alloc(1,sizeof(*w));
44+
w->n = n1;
45+
w->na = 2*nw+1;
46+
47+
w->a = sf_floatalloc2 (n1,w->na);
48+
w->b = sf_floatalloc (w->na);
49+
50+
return w;
51+
}
52+
53+
void cnv_close (cnv w)
54+
/*< free allocated storage >*/
55+
{
56+
free (w->a[0]);
57+
free (w->a);
58+
free (w->b);
59+
free (w);
60+
}
61+
62+
void cnv_define (bool adj /* adjoint flag */,
63+
cnv w /* cnv object */,
64+
float** cc /* convection */,
65+
float* diag /* defined diagonal */,
66+
float** offd /* defined off-diagonal */)
67+
/*< fill the matrix >*/
68+
{
69+
int i, j, k, m, n, iw, nw;
70+
float am, aj;
71+
72+
nw = (w->na-1)/2;
73+
n = w->n;
74+
75+
for (i=0; i < n; i++) {
76+
w->b[nw] = 1.0f;
77+
for (iw=1; iw <= nw; iw++) {
78+
w->b[nw-iw] = cc[2*iw-1][i];
79+
w->b[nw+iw] = cc[2*iw-2][i];
80+
w->b[nw] -= w->b[nw-iw] + w->b[nw+iw];
81+
}
82+
83+
for (j=0; j < w->na; j++) {
84+
if (adj) {
85+
w->a[j][i] = w->b[w->na-1-j];
86+
} else {
87+
w->a[j][i] = w->b[j];
88+
}
89+
}
90+
}
91+
92+
for (i=0; i < n; i++) {
93+
for (j=0; j < w->na; j++) {
94+
k = i+j-nw;
95+
if (k >= nw && k < n-nw) {
96+
aj = w->a[j][k];
97+
diag[i] += aj*aj;
98+
}
99+
}
100+
for (m=0; m < 2*nw; m++) {
101+
for (j=m+1; j < w->na; j++) {
102+
k = i+j-nw;
103+
if (k >= nw && k < n-nw) {
104+
aj = w->a[j][k];
105+
am = w->a[j-m-1][k];
106+
offd[m][i] += am*aj;
107+
}
108+
}
109+
}
110+
}
111+
}
112+
113+
void cnv_set (bool adj /* adjoint flag */,
114+
cnv w /* cnv object */,
115+
float* inp /* input */,
116+
float* out /* output */,
117+
float* tmp /* temporary storage */)
118+
/*< matrix multiplication >*/
119+
{
120+
int i, j, k, n, nw;
121+
122+
nw = (w->na-1)/2;
123+
n = w->n;
124+
125+
if (adj) {
126+
for (i=0; i < n; i++) {
127+
tmp[i]=0.;
128+
}
129+
for (i=0; i < n; i++) {
130+
for (j=0; j < w->na; j++) {
131+
k = i+j-nw;
132+
if (k >= nw && k < n-nw)
133+
tmp[k] += w->a[j][k]*out[i];
134+
}
135+
}
136+
for (i=0; i < n; i++) {
137+
inp[i]=0.;
138+
}
139+
for (i=nw; i < n-nw; i++) {
140+
for (j=0; j < w->na; j++) {
141+
k = i+j-nw;
142+
inp[k] += w->a[j][i]*tmp[i];
143+
}
144+
}
145+
} else {
146+
for (i=0; i < n; i++) {
147+
tmp[i] = 0.;
148+
}
149+
for (i=nw; i < n-nw; i++) {
150+
for (j=0; j < w->na; j++) {
151+
k = i+j-nw;
152+
tmp[i] += w->a[j][i]*inp[k];
153+
}
154+
}
155+
for (i=0; i < n; i++) {
156+
out[i] = 0.;
157+
for (j=0; j < w->na; j++) {
158+
k = i+j-nw;
159+
if (k >= nw && k < n-nw)
160+
out[i] += w->a[j][k]*tmp[k];
161+
}
162+
}
163+
}
164+
}

0 commit comments

Comments
 (0)