-
Notifications
You must be signed in to change notification settings - Fork 11
/
Toa_smooth3d.c
executable file
·172 lines (135 loc) · 9.09 KB
/
Toa_smooth3d.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
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
#include<stdio.h>
#include<malloc.h>
#include<math.h>
#include<stdlib.h>
#include "/home/Toa/hc/cjbsegy.h"
#include "/home/Toa/hc/fft.c"
#include "/home/Toa/hc/alloc.c"
#include "/home/Toa/hc/complex.c"
void main()
{
float ***vel;
int i,j,k;
int nx,ny,nz,nn,mmz,n;
char FN1[250]={"vel707070.dat"};
char FN2[250]={"vel707070initial.dat"};
FILE *fpvel;
fpvel=fopen(FN1,"rb");
FILE *fpsmooth;
fpsmooth=fopen(FN2,"wb");
nx=70;
ny=70;
nz=70;
nn=20;
mmz=20;
vel=alloc3float(nz,ny,nx);
zero3float(vel,nz,ny,nx);
for(j=0;j<ny;j++)
{
for(i=0;i<nx;i++)
{
for(k=0;k<nz;k++)
{
fread(&vel[i][j][k],4L,1,fpvel);
}
}
}
/********************************************** smooth3 **********************************/
for(n=0;n<nn;n++)
{
//************* inner
for(i=1;i<nx-1;i++)
{
for(j=1;j<ny-1;j++)
{
for(k=mmz;k<nz-1;k++)
{
vel[i][j][k]=((vel[i-1][j-1][k-1]+vel[i-1][j][k-1]+vel[i-1][j+1][k-1]+vel[i][j-1][k-1]+vel[i][j][k-1]+vel[i][j+1][k-1]+vel[i+1][j-1][k-1]+vel[i+1][j][k-1]+vel[i+1][j+1][k-1])+(vel[i-1][j-1][k]+vel[i-1][j][k]+vel[i-1][j+1][k]+vel[i][j-1][k]+vel[i][j][k]+vel[i][j+1][k]+vel[i+1][j-1][k]+vel[i+1][j][k]+vel[i+1][j+1][k])+(vel[i-1][j-1][k+1]+vel[i-1][j][k+1]+vel[i-1][j+1][k+1]+vel[i][j-1][k+1]+vel[i][j][k+1]+vel[i][j+1][k+1]+vel[i+1][j-1][k+1]+vel[i+1][j][k+1]+vel[i+1][j+1][k+1]))/27;
}
}
}
//************* up and down
for(j=1;j<ny-1;j++)
{
for(i=1;i<nx-1;i++)
{
vel[i][j][mmz]=((vel[i-1][j-1][mmz]+vel[i-1][j][mmz]+vel[i-1][j+1][mmz]+vel[i][j-1][mmz]+vel[i][j][mmz]+vel[i][j+1][mmz]+vel[i+1][j-1][mmz]+vel[i+1][j][mmz]+vel[i+1][j+1][mmz])+(vel[i-1][j-1][mmz+1]+vel[i-1][j][mmz+1]+vel[i-1][j+1][mmz+1]+vel[i][j-1][mmz+1]+vel[i][j][mmz+1]+vel[i][j+1][mmz+1]+vel[i+1][j-1][mmz+1]+vel[i+1][j][mmz+1]+vel[i+1][j+1][mmz+1]))/18;
vel[i][j][nz-1]=((vel[i-1][j-1][nz-1]+vel[i-1][j][nz-1]+vel[i-1][j+1][nz-1]+vel[i][j-1][nz-1]+vel[i][j][nz-1]+vel[i][j+1][nz-1]+vel[i+1][j-1][nz-1]+vel[i+1][j][nz-1]+vel[i+1][j+1][nz-1])+(vel[i-1][j-1][nz-2]+vel[i-1][j][nz-2]+vel[i-1][j+1][nz-2]+vel[i][j-1][nz-2]+vel[i][j][nz-2]+vel[i][j+1][nz-2]+vel[i+1][j-1][nz-2]+vel[i+1][j][nz-2]+vel[i+1][j+1][nz-2]))/18;
}
}
//************* left and right
for(j=1;j<ny-1;j++)
{
for(k=mmz;k<nz-1;k++)
{
vel[0][j][k]=((vel[0][j-1][k-1]+vel[0][j-1][k]+vel[0][j-1][k+1]+vel[0][j][k-1]+vel[0][j][k]+vel[0][j][k+1]+vel[0][j+1][k-1]+vel[0][j+1][k]+vel[0][j+1][k+1])+(vel[1][j-1][k-1]+vel[1][j-1][k]+vel[1][j-1][k+1]+vel[1][j][k-1]+vel[1][j][k]+vel[1][j][k+1]+vel[1][j+1][k-1]+vel[1][j+1][k]+vel[1][j+1][k+1]))/18;
vel[nx-1][j][k]=((vel[nx-1][j-1][k-1]+vel[nx-1][j-1][k]+vel[nx-1][j-1][k+1]+vel[nx-1][j][k-1]+vel[nx-1][j][k]+vel[nx-1][j][k+1]+vel[nx-1][j+1][k-1]+vel[nx-1][j+1][k]+vel[nx-1][j+1][k+1])+(vel[nx-2][j-1][k-1]+vel[nx-2][j-1][k]+vel[nx-2][j-1][k+1]+vel[nx-2][j][k-1]+vel[nx-2][j][k]+vel[nx-2][j][k+1]+vel[nx-2][j+1][k-1]+vel[nx-2][j+1][k]+vel[nx-2][j+1][k+1]))/18;
}
}
//************** front and back
for(i=1;i<nx-1;i++)
{
for(k=mmz;k<nz;k++)
{
vel[i][0][k]=((vel[i-1][0][k-1]+vel[i-1][0][k]+vel[i-1][0][k+1]+vel[i][0][k-1]+vel[i][0][k]+vel[i][0][k+1]+vel[i+1][0][k-1]+vel[i+1][0][k]+vel[i+1][0][k+1])+(vel[i-1][1][k-1]+vel[i-1][1][k]+vel[i-1][1][k+1]+vel[i][1][k-1]+vel[i][1][k]+vel[i][1][k+1]+vel[i+1][1][k-1]+vel[i+1][1][k]+vel[i+1][1][k+1]))/18;
vel[i][ny-1][k]=((vel[i-1][ny-1][k-1]+vel[i-1][ny-1][k]+vel[i-1][ny-1][k+1]+vel[i][ny-1][k-1]+vel[i][ny-1][k]+vel[i][ny-1][k+1]+vel[i+1][ny-1][k-1]+vel[i+1][ny-1][k]+vel[i+1][ny-1][k+1])+(vel[i-1][ny-2][k-1]+vel[i-1][ny-2][k]+vel[i-1][ny-2][k+1]+vel[i][ny-2][k-1]+vel[i][ny-2][k]+vel[i][ny-2][k+1]+vel[i+1][ny-2][k-1]+vel[i+1][ny-2][k]+vel[i+1][ny-2][k+1]))/18;
}
}
//********* the eight point
vel[0][0][mmz]=(vel[1][0][mmz]+vel[0][1][mmz]+vel[1][1][mmz]+vel[0][0][mmz]+vel[1][0][mmz+1]+vel[0][1][mmz+1]+vel[1][1][mmz+1]+vel[0][0][mmz+1])/8;
vel[0][0][nz-1]=(vel[1][0][nz-1]+vel[0][1][nz-1]+vel[1][1][nz-1]+vel[0][0][nz-1]+vel[1][0][nz-2]+vel[0][1][nz-2]+vel[1][1][nz-2]+vel[0][0][nz-2])/8;
vel[nx-1][0][mmz]=(vel[nx-2][0][mmz]+vel[nx-1][1][mmz]+vel[nx-2][1][mmz]+vel[nx-1][0][mmz]+vel[nx-2][0][mmz+1]+vel[nx-1][1][mmz+1]+vel[nx-2][1][mmz+1]+vel[nx-1][0][mmz+1])/8;
vel[nx-1][0][nz-1]=(vel[nx-2][0][nz-1]+vel[nx-1][1][nz-1]+vel[nx-2][1][nz-1]+vel[nx-1][0][nz-1]+vel[nx-2][0][nz-2]+vel[nx-1][1][nz-2]+vel[nx-2][1][nz-2]+vel[nx-1][0][nz-2])/8;
vel[0][ny-1][mmz]=(vel[1][ny-1][mmz]+vel[0][ny-2][mmz]+vel[1][ny-2][mmz]+vel[0][ny-1][mmz]+vel[1][ny-1][mmz+1]+vel[0][ny-2][mmz+1]+vel[1][ny-2][mmz+1]+vel[0][ny-1][mmz+1])/8;
vel[0][ny-1][nz-1]=(vel[1][ny-1][nz-1]+vel[0][ny-2][nz-1]+vel[1][ny-2][nz-1]+vel[0][ny-1][nz-1]+vel[1][ny-1][nz-2]+vel[0][ny-2][nz-2]+vel[0][ny-1][nz-2]+vel[1][ny-2][nz-2])/8;
vel[nx-1][ny-1][mmz]=(vel[nx-2][ny-1][mmz]+vel[nx-1][ny-2][mmz]+vel[nx-2][ny-2][mmz]+vel[nx-1][ny-1][mmz]+vel[nx-2][ny-1][mmz+1]+vel[nx-1][ny-2][mmz+1]+vel[nx-2][ny-2][mmz+1]+vel[nx-1][ny-1][mmz+1])/8;
vel[nx-1][ny-1][nz-1]=(vel[nx-2][ny-1][nz-1]+vel[nx-1][ny-2][nz-1]+vel[nx-2][ny-2][nz-1]+vel[nx-1][ny-1][nz-1]+vel[nx-2][ny-1][nz-2]+vel[nx-1][ny-2][nz-2]+vel[nx-2][ny-2][nz-2]+vel[nx-1][ny-1][nz-2])/8;
for(k=mmz;k<nz-1;k++)
{
//************* x=0 & y=0
vel[0][0][k]=((vel[1][0][k]+vel[0][1][k]+vel[1][1][k]+vel[0][0][k])+(vel[1][0][k-1]+vel[0][1][k-1]+vel[1][1][k-1]+vel[0][0][k-1])+(vel[1][0][k+1]+vel[0][1][k+1]+vel[1][1][k+1]+vel[0][0][k+1]))/12;
//************* x=0 & y=ny-1
vel[0][ny-1][k]=((vel[1][ny-1][k]+vel[0][ny-2][k]+vel[1][ny-2][k]+vel[0][ny-1][k])+(vel[1][ny-1][k-1]+vel[0][ny-2][k-1]+vel[1][ny-2][k-1]+vel[0][ny-1][k-1])+(vel[1][ny-1][k+1]+vel[0][ny-2][k+1]+vel[1][ny-2][k+1]+vel[0][ny-1][k+1]))/12;
//************* x=nx-1 & y=0
vel[nx-1][0][k]=((vel[nx-2][0][k]+vel[nx-1][1][k]+vel[nx-2][1][k]+vel[nx-1][0][k])+(vel[nx-2][0][k-1]+vel[nx-1][1][k-1]+vel[nx-2][1][k-1]+vel[nx-1][0][k-1])+(vel[nx-2][0][k+1]+vel[nx-1][1][k+1]+vel[nx-2][1][k+1]+vel[nx-1][0][k+1]))/12;
//************* x=nx-1 & y=ny-1
vel[nx-1][ny-1][k]=((vel[nx-2][ny-1][k]+vel[nx-1][ny-2][k]+vel[nx-2][ny-2][k]+vel[nx-1][ny-1][k])+(vel[nx-2][ny-1][k-1]+vel[nx-1][ny-2][k-1]+vel[nx-2][ny-2][k-1]+vel[nx-1][ny-1][k-1])+(vel[nx-2][ny-1][k+1]+vel[nx-1][ny-2][k+1]+vel[nx-2][ny-2][k+1]+vel[nx-1][ny-1][k+1]))/12;
}
for(i=1;i<nx-1;i++)
{
//************** y=0 & z=mmz
vel[i][0][mmz]=((vel[i][1][mmz]+vel[i][0][mmz+1]+vel[i][1][mmz+1]+vel[i][0][mmz])+(vel[i-1][1][mmz]+vel[i-1][0][mmz+1]+vel[i-1][1][mmz+1]+vel[i-1][0][mmz])+(vel[i+1][1][mmz]+vel[i+1][0][mmz+1]+vel[i+1][1][mmz+1]+vel[i+1][0][mmz]))/12;
//************** y=0 & z=nz-1
vel[i][0][nz-1]=((vel[i][1][nz-1]+vel[i][0][nz-2]+vel[i][1][nz-2]+vel[i][0][nz-1])+(vel[i-1][1][nz-1]+vel[i-1][0][nz-2]+vel[i-1][1][nz-2]+vel[i-1][0][nz-1])+(vel[i+1][1][nz-1]+vel[i+1][0][nz-2]+vel[i+1][1][nz-2]+vel[i+1][0][nz-1]))/12;
//************** y=ny-1 & z=mmz
vel[i][ny-1][mmz]=((vel[i][ny-2][mmz]+vel[i][ny-1][mmz+1]+vel[i][ny-2][mmz+1]+vel[i][ny-1][mmz])+(vel[i-1][ny-2][mmz]+vel[i-1][ny-1][mmz+1]+vel[i-1][ny-2][mmz+1]+vel[i-1][ny-1][mmz])+(vel[i+1][ny-2][mmz]+vel[i+1][ny-1][mmz+1]+vel[i+1][ny-2][mmz+1]+vel[i+1][ny-1][mmz]))/12;
//************** y=ny-1 & z=nz-1
vel[i][ny-1][nz-1]=((vel[i][ny-2][nz-1]+vel[i][ny-1][nz-2]+vel[i][ny-2][nz-2]+vel[i][ny-1][nz-1])+(vel[i-1][ny-2][nz-1]+vel[i-1][ny-1][nz-2]+vel[i-1][ny-2][nz-2]+vel[i-1][ny-1][nz-1])+(vel[i+1][ny-2][nz-1]+vel[i+1][ny-1][nz-2]+vel[i+1][ny-2][nz-2]+vel[i+1][ny-1][nz-1]))/12;
}
for(j=1;j<ny-1;j++)
{
//************** x=0 & z=mmz
vel[0][j][mmz]=((vel[1][j][mmz]+vel[0][j][mmz+1]+vel[1][j][mmz+1]+vel[0][j][mmz])+(vel[1][j-1][mmz]+vel[0][j-1][mmz+1]+vel[1][j-1][mmz+1]+vel[0][j-1][mmz])+(vel[1][j+1][mmz]+vel[0][j+1][mmz+1]+vel[1][j+1][mmz+1]+vel[0][j+1][mmz]))/12;
//************** x=0 & z=nz-1
vel[0][j][nz-1]=((vel[1][j][nz-1]+vel[0][j][nz-2]+vel[1][j][nz-2]+vel[0][j][nz-1])+(vel[1][j-1][nz-1]+vel[0][j-1][nz-2]+vel[1][j-1][nz-2]+vel[0][j-1][nz-1])+(vel[1][j+1][nz-1]+vel[0][j+1][nz-2]+vel[1][j+1][nz-2]+vel[0][j+1][nz-1]))/12;
//************** x=nx-1 & z=mmz
vel[nx-1][j][mmz]=((vel[nx-2][j][mmz]+vel[nx-1][j][mmz+1]+vel[nx-2][j][mmz+1]+vel[nx-1][j][mmz])+(vel[nx-2][j-1][mmz]+vel[nx-1][j-1][mmz+1]+vel[nx-2][j-1][mmz+1]+vel[nx-1][j-1][mmz])+(vel[nx-2][j+1][mmz]+vel[nx-1][j+1][mmz+1]+vel[nx-2][j+1][mmz+1]+vel[nx-1][j+1][mmz]))/12;
//************** x=nx-1 &z=nz-1
vel[nx-1][j][nz-1]=((vel[nx-2][j][nz-1]+vel[nx-1][j][nz-2]+vel[nx-2][j][nz-2]+vel[nx-1][j][nz-1])+(vel[nx-2][j-1][nz-1]+vel[nx-1][j-1][nz-2]+vel[nx-2][j-1][nz-2]+vel[nx-1][j-1][nz-1])+(vel[nx-2][j+1][nz-1]+vel[nx-1][j+1][nz-2]+vel[nx-2][j+1][nz-2]+vel[nx-1][j+1][nz-1]))/12;
}
}
/*************************************************************************************************/
for(j=0;j<ny;j++)
{
for(i=0;i<nx;i++)
{
for(k=0;k<nz;k++)
{
vel[i][j][k]=vel[i][j][k];
fwrite(&vel[i][j][k],4L,1,fpsmooth);
}
}
}
fclose(fpvel);
free3float(vel);
}