Skip to content
Snippets Groups Projects
Commit 6225c9d4 authored by joeri.brackenhoff's avatar joeri.brackenhoff
Browse files

3D

parent 617c2b82
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
fdelmodc3D/grid3D.c 0 → 100644
+ 243
− 0
View file @ 6225c9d4
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
/**
* fills the gridded model below the interface zp used in makemod.
* Vertical and horizontal gradients are taken into account
*
* AUTHOR:
* Jan Thorbecke (janth@xs4all.nl)
* The Netherlands
**/
#ifndef ABS
#define ABS(x) ((x) < 0 ? -(x) : (x))
#endif
#ifndef MAX
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#endif
#ifndef MIN
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#endif
#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
void grid3D(float **gridcp, float **gridcs, float **gridro, long *zp,
float **cp, float **cs, float **ro, long minx, long maxx, long miny, long maxy,
long optgrad, float gradlen, float gradcp, float gradcs, float gradro, float dx,
float dy, float dz, long nz)
{
long g, ngrad, gradend, i, k, j, l;
float deltcp, deltcs, deltro, co;
if (gridcs == NULL && gridro == NULL) {
if (optgrad == 1) {
g = NINT(gradlen/dz)-1;
for (l = miny; l < maxy; l++) {
for (i = minx; i < maxx; i++) {
deltcp = (cp[1][l*maxx+i] - cp[0][l*maxx+i])/g;
if (zp[l*maxx+i] == 0) k = 1;
else k = 0;
ngrad = zp[l*maxx+i] + g;
gradend = MIN(ngrad, nz);
for (j = zp[l*maxx+i]; j <= gradend; j++) {
gridcp[l*maxx+i][j] = cp[0][l*maxx+i]+deltcp*k;
k += 1;
}
for (j = gradend+1, k = 0; j < nz; j++, k++) {
gridcp[l*maxx+i][j] = cp[1][l*maxx+i] + (float)k*gradcp;
}
}
}
}
else if (optgrad == 2) {
g = NINT(gradlen/dz)-1;
for (l = miny; l < maxy; l++) {
for (i = minx; i < maxx; i++) {
deltcp = (cp[1][l*maxx+i] - cp[0][l*maxx+i])/g;
if (zp[l*maxx+i] == 0) k = 1;
else k = 0;
ngrad = zp[l*maxx+i] + g;
gradend = MIN(ngrad, nz);
for (j = zp[l*maxx+i]; j <= gradend; j++) {
co = -g*(cos(M_PI*((float)k/(float)g))-1)/2.0;
gridcp[l*maxx+i][j] = cp[0][l*maxx+i]+deltcp*co;
k += 1;
}
for (j = gradend+1, k = 0; j < nz; j++, k++) {
gridcp[l*maxx+i][j] = cp[1][l*maxx+i] + (float)k*gradcp;
}
}
}
}
else if (optgrad == 3) {
for (l = miny; l < maxy; l++) {
for (i = minx; i < maxx; i++) {
for (j = zp[l*maxx+i], k = 0; j < nz; j++, k++) {
gridcp[l*maxx+i][j] = cp[1][l*maxx+i] + (float)k*gradcp;
}
}
}
}
else if (optgrad == 4) {
for (l = miny; l < maxy; l++) {
for (i = minx; i < maxx; i++) {
for (j = zp[l*maxx+i], k = 0; j < nz; j++, k++) {
gridcp[l*maxx+i][j] = cp[1][l*maxx+i] + (float)(drand48()-0.5)*2.0*gradcp;
}
}
}
}
return;
}
else if (gridcs == NULL) {
if (optgrad == 1) {
g = NINT(gradlen/dz)-1;
for (l = miny; l < maxy; l++) {
for (i = minx; i < maxx; i++) {
deltcp = (cp[1][l*maxx+i] - cp[0][l*maxx+i])/g;
deltro = (ro[1][l*maxx+i] - ro[0][l*maxx+i])/g;
if (zp[l*maxx+i] == 0) k = 1;
else k = 0;
ngrad = zp[l*maxx+i] + g;
gradend = MIN(ngrad, nz);
for (j = zp[l*maxx+i]; j <= gradend; j++) {
gridcp[l*maxx+i][j] = cp[0][l*maxx+i]+deltcp*k;
gridro[l*maxx+i][j] = ro[0][l*maxx+i]+deltro*k;
k += 1;
}
for (j = gradend+1, k = 0; j < nz; j++, k++) {
gridcp[l*maxx+i][j] = cp[1][l*maxx+i] + (float)k*gradcp;
gridro[l*maxx+i][j] = ro[1][l*maxx+i] + (float)k*gradro;
}
}
}
}
else if (optgrad == 2) {
g = NINT(gradlen/dz)-1;
for (l = miny; l < maxy; l++) {
for (i = minx; i < maxx; i++) {
deltcp = (cp[1][l*maxx+i] - cp[0][l*maxx+i])/g;
deltro = (ro[1][l*maxx+i] - ro[0][l*maxx+i])/g;
if (zp[l*maxx+i] == 0) k = 1;
else k = 0;
ngrad = zp[l*maxx+i] + g;
gradend = MIN(ngrad, nz);
for (j = zp[l*maxx+i]; j <= gradend; j++) {
co = -g*(cos(M_PI*((float)k/(float)g))-1)/2.0;
gridcp[l*maxx+i][j] = cp[0][l*maxx+i]+deltcp*co;
gridro[l*maxx+i][j] = ro[0][l*maxx+i]+deltro*co;
k += 1;
}
for (j = gradend+1, k = 0; j < nz; j++, k++) {
gridcp[l*maxx+i][j] = cp[1][l*maxx+i] + (float)k*gradcp;
gridro[l*maxx+i][j] = ro[1][l*maxx+i] + (float)k*gradro;
}
}
}
}
else if (optgrad == 3) {
for (l = miny; l < maxy; l++) {
for (i = minx; i < maxx; i++) {
for (j = zp[l*maxx+i], k = 0; j < nz; j++, k++) {
gridcp[l*maxx+i][j] = cp[1][l*maxx+i] + (float)k*gradcp;
gridro[l*maxx+i][j] = ro[1][l*maxx+i] + (float)k*gradro;
}
}
}
}
else if (optgrad == 4) {
for (l = miny; l < maxy; l++) {
for (i = minx; i < maxx; i++) {
for (j = zp[l*maxx+i], k = 0; j < nz; j++, k++) {
gridcp[l*maxx+i][j] = cp[1][l*maxx+i] + (float)(drand48()-0.5)*2.0*gradcp;
gridro[l*maxx+i][j] = ro[1][l*maxx+i] + (float)(drand48()-0.5)*2.0*gradro;
}
}
}
}
return;
}
else {
if (optgrad == 1) {
g = NINT(gradlen/dz)-1;
for (l = miny; l < maxy; l++) {
for (i = minx; i < maxx; i++) {
deltcp = (cp[1][l*maxx+i] - cp[0][l*maxx+i])/g;
deltcs = (cs[1][l*maxx+i] - cs[0][l*maxx+i])/g;
deltro = (ro[1][l*maxx+i] - ro[0][l*maxx+i])/g;
if (zp[l*maxx+i] == 0) k = 1;
else k = 0;
ngrad = zp[l*maxx+i] + g;
gradend = MIN(ngrad, nz);
for (j = zp[l*maxx+i]; j <= gradend; j++) {
gridcp[l*maxx+i][j] = cp[0][l*maxx+i]+deltcp*k;
gridcs[l*maxx+i][j] = cs[0][l*maxx+i]+deltcs*k;
gridro[l*maxx+i][j] = ro[0][l*maxx+i]+deltro*k;
k += 1;
}
for (j = gradend+1, k = 0; j < nz; j++, k++) {
gridcp[l*maxx+i][j] = cp[1][i] + (float)k*gradcp;
gridcs[l*maxx+i][j] = cs[1][i] + (float)k*gradcs;
gridro[l*maxx+i][j] = ro[1][i] + (float)k*gradro;
}
}
}
}
else if (optgrad == 2) {
g = NINT(gradlen/dz)-1;
for (l = miny; l < maxy; l++) {
for (i = minx; i < maxx; i++) {
deltcp = (cp[1][l*maxx+i] - cp[0][l*maxx+i])/g;
deltcs = (cs[1][l*maxx+i] - cs[0][l*maxx+i])/g;
deltro = (ro[1][l*maxx+i] - ro[0][l*maxx+i])/g;
if (zp[l*maxx+i] == 0) k = 1;
else k = 0;
ngrad = zp[l*maxx+i] + g;
gradend = MIN(ngrad, nz);
for (j = zp[l*maxx+i]; j <= gradend; j++) {
co = -g*(cos(M_PI*((float)k/(float)g))-1)/2.0;
gridcp[l*maxx+i][j] = cp[0][l*maxx+i]+deltcp*co;
gridcs[l*maxx+i][j] = cs[0][l*maxx+i]+deltcs*co;
gridro[l*maxx+i][j] = ro[0][l*maxx+i]+deltro*co;
k += 1;
}
for (j = gradend+1, k = 0; j < nz; j++, k++) {
gridcp[l*maxx+i][j] = cp[1][l*maxx+i] + (float)k*gradcp;
gridcs[l*maxx+i][j] = cs[1][l*maxx+i] + (float)k*gradcs;
gridro[l*maxx+i][j] = ro[1][l*maxx+i] + (float)k*gradro;
}
}
}
}
else if (optgrad == 3) {
for (l = miny; l < maxy; l++) {
for (i = minx; i < maxx; i++) {
for (j = zp[l*maxx+i], k = 0; j < nz; j++, k++) {
gridcp[l*maxx+i][j] = cp[1][l*maxx+i] + (float)k*gradcp;
gridcs[l*maxx+i][j] = cs[1][l*maxx+i] + (float)k*gradcs;
gridro[l*maxx+i][j] = ro[1][l*maxx+i] + (float)k*gradro;
}
}
}
}
else if (optgrad == 4) {
for (l = miny; l < maxy; l++) {
for (i = minx; i < maxx; i++) {
for (j = zp[l*maxx+i], k = 0; j < nz; j++, k++) {
gridcp[l*maxx+i][j] = cp[1][l*maxx+i] + (float)(drand48()-0.5)*2.0*gradcp;
gridcs[l*maxx+i][j] = cs[1][l*maxx+i] + (float)(drand48()-0.5)*2.0*gradcs;
gridro[l*maxx+i][j] = ro[1][l*maxx+i] + (float)(drand48()-0.5)*2.0*gradro;
}
}
}
}
return;
}
}
fdelmodc3D/interpolation3D.c 0 → 100644
+ 231
− 0
View file @ 6225c9d4
#include <math.h>
#include <stdlib.h>
/**
* Interpolates the interface defined by the input parameters to all grid points
* 4 different interpolation schemes can be chosen
* - linear
* - polynomal
* - cubic spline
* Used in makemod
*
* AUTHOR:
* Jan Thorbecke (janth@xs4all.nl)
* The Netherlands
**/
#ifndef ABS
#define ABS(x) ((x) < 0 ? -(x) : (x))
#endif
#ifndef MAX
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#endif
#ifndef MIN
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#endif
#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
void polint(float xa[], float ya[], long n, float x, float *y, float *dy);
void spline(float x[], float y[], long n, float yp1, float ypn, float y2[]);
void splint(float xa[], float ya[], float y2a[], long n, float x, float *y);
void interpolation(float *x, float *y, float *z, long nxp, long nx, long ny,
long poly, long *pminx, long *pmaxx, float dx, float dy, float **cp, float **cs, float **ro,
long nvel, float *interface)
{
long i, j, l, ndeltx, np;
float deltx, delty, deltz, xprev, yprev, zprev, minx, maxx, miny, maxy;
float deltcp, deltcs, deltro, cpprev, csprev, roprev;
float *xa, *ya, *za, dyp, xp, yp1, ypn, *y2;
if (poly == 0) {
np = 0;
xprev = yprev = zprev = 0.0;
minx = nx*dx;
miny = ny*dy;
maxx = 0;
maxy = 0;
for (i = 0; i < nxp; i++) {
if (x[i] < minx) {
xprev = x[i];
minx = x[i];
*pminx = NINT(minx/dx);
np = *pminx;
}
if (x[i] > maxx) {
maxx = x[i];
*pmaxx = MIN(NINT((maxx+dx)/dx),nx);
}
if (y[i] < miny) {
yprev = y[i];
miny = y[i];
*pminy = NINT(miny/dy);
np = *pminy;
}
if (x[i] > maxx) {
maxx = x[i];
*pmaxx = MIN(NINT((maxx+dx)/dx),nx);
}
deltx = x[i] - xprev;
deltz = z[i] - zprev;
if (i == 0) ndeltx = -1;
else ndeltx = NINT(ABS(deltx/dx));
for (j = 0; j < ndeltx && np < nx; j++) {
interface[np++] = zprev + (j*dx*deltz)/deltx;
}
xprev = x[i];
zprev = z[i];
}
for (j = np; j < *pmaxx; j++) interface[j] = z[nxp-1];
}
else if (poly == 1) {
xa = (float *)malloc((nxp+1)*sizeof(float));
za = (float *)malloc((nxp+1)*sizeof(float));
for (i = 1; i <= nxp; i++) xa[i] = x[i-1];
for (i = 1; i <= nxp; i++) za[i] = z[i-1];
np = 0;
minx = nx*dx;
maxx = 0;
for (i = 0; i < nxp; i++) {
if (x[i] < minx) {
xprev = x[i];
minx = x[i];
*pminx = NINT(minx/dx);
np = *pminx;
}
if (x[i] > maxx) {
maxx = x[i];
*pmaxx = MIN(NINT((maxx+dx)/dx),nx);
}
}
for (j = *pminx; j < *pmaxx; j++) {
xp = j*dx;
polint(xa, za, nxp, xp, &interface[j], &dyp);
}
free(xa);
free(za);
}
else if (poly == 2) {
xa = (float *)malloc((nxp+1)*sizeof(float));
za = (float *)malloc((nxp+1)*sizeof(float));
for (i = 1; i <= nxp; i++) xa[i] = x[i-1];
for (i = 1; i <= nxp; i++) za[i] = z[i-1];
np = 0;
minx = nx*dx;
maxx = 0;
for (i = 0; i < nxp; i++) {
if (x[i] < minx) {
xprev = x[i];
minx = x[i];
*pminx = NINT(minx/dx);
np = *pminx;
}
if (x[i] > maxx) {
maxx = x[i];
*pmaxx = MIN(NINT((maxx+dx)/dx),nx);
}
}
y2 = (float *)malloc((nxp+1)*sizeof(float));
yp1 = ypn = 1e30;
spline(xa, za, nxp, yp1, ypn, y2);
for (j = *pminx; j < *pmaxx; j++) {
xp = j*dx;
splint(xa, za, y2, nxp, xp, &interface[j]);
}
free(y2);
free(xa);
free(za);
}
if (nvel != 1) {
if (nvel != 2) {
for (j = 0; j < nx; j++) {
cp[0][j] = cp[1][j];
cs[0][j] = cs[1][j];
ro[0][j] = ro[1][j];
}
np = 0;
xprev = cpprev = csprev = roprev = 0.0;
minx = nx*dx;
maxx = 0;
for (i = 0; i < nxp; i++) {
if (x[i] < minx) {
xprev = x[i];
minx = x[i];
*pminx = NINT(minx/dx);
np = *pminx;
}
if (x[i] > maxx) {
maxx = x[i];
*pmaxx = MIN(NINT((maxx+dx)/dx),nx);
}
deltx = x[i] - xprev;
deltcp = cp[2][i] - cpprev;
deltcs = cs[2][i] - csprev;
deltro = ro[2][i] - roprev;
if (i == 0) ndeltx = -1;
else ndeltx = NINT(ABS(deltx/dx))-1;
if (ndeltx==0) ndeltx = 1;
for (j = 0; j <= ndeltx; j++) {
cp[1][np] = cpprev + (j*dx*deltcp)/deltx;
cs[1][np] = csprev + (j*dx*deltcs)/deltx;
ro[1][np] = roprev + (j*dx*deltro)/deltx;
np += 1;
}
xprev = x[i];
cpprev = cp[2][i];
csprev = cs[2][i];
roprev = ro[2][i];
}
cp[1][np] = cpprev;
cs[1][np] = csprev;
ro[1][np] = roprev;
}
else {
for (j = 0; j < nx; j++) {
cp[0][j] = cp[1][j];
cs[0][j] = cs[1][j];
ro[0][j] = ro[1][j];
}
np = 0;
xprev = 0.0;
minx = nx*dx;
maxx = 0;
cpprev = cp[2][0];
csprev = cs[2][0];
roprev = ro[2][0];
deltcp = cp[2][1] - cpprev;
deltcs = cs[2][1] - csprev;
deltro = ro[2][1] - roprev;
minx = x[0];
maxx = x[nxp-1];
deltx = maxx - minx;
np = NINT(minx/dx);
ndeltx = NINT(ABS(deltx/dx));
for (j = 0; j <= ndeltx; j++) {
cp[1][np] = cpprev + (j*dx*deltcp)/deltx;
cs[1][np] = csprev + (j*dx*deltcs)/deltx;
ro[1][np] = roprev + (j*dx*deltro)/deltx;
np += 1;
}
}
}
else {
for (j = 0; j < nx; j++) {
cp[0][j] = cp[1][j];
cs[0][j] = cs[1][j];
ro[0][j] = ro[1][j];
cp[1][j] = cp[2][0];
cs[1][j] = cs[2][0];
ro[1][j] = ro[2][0];
}
}
return;
}
fdelmodc3D/makemod3D.c 0 → 100644
+ 1004
− 0
View file @ 6225c9d4
This diff is collapsed.
marchenko3D/marchenko3D.c
+ 29
− 7
View file @ 6225c9d4
...@@ -742,13 +742,24 @@ int main (int argc, char **argv) ...@@ -742,13 +742,24 @@ int main (int argc, char **argv)
for (l=0; l<nyim; l++){ for (l=0; l<nyim; l++){
for (j=0; j<nxim; j++){ for (j=0; j<nxim; j++){
hdrs_Nfoc[l*nxim+j].ns = nzim; hdrs_Nfoc[l*nxim+j].ns = nzim;
hdrs_Nfoc[l*nxim+j].sx = xsyn[j]; hdrs_Nfoc[l*nxim+j].fldr = 1;
hdrs_Nfoc[l*nxim+j].sy = ysyn[l]; hdrs_Nfoc[l*nxim+j].tracl = 1;
hdrs_Nfoc[l*nxim+j].sdepth = zsyn[l]; hdrs_Nfoc[l*nxim+j].tracf = l*nxim+j+1;
hdrs_Nfoc[l*nxim+j].trid = 1;
hdrs_Nfoc[l*nxim+j].scalco = -1000;
hdrs_Nfoc[l*nxim+j].scalel = -1000;
hdrs_Nfoc[l*nxim+j].sx = xsyn[j]*(1e3);
hdrs_Nfoc[l*nxim+j].sy = ysyn[l]*(1e3);
hdrs_Nfoc[l*nxim+j].gx = xsyn[j]*(1e3);
hdrs_Nfoc[l*nxim+j].gy = ysyn[l]*(1e3);
hdrs_Nfoc[l*nxim+j].sdepth = zsyn[l]*(1e3);
hdrs_Nfoc[l*nxim+j].f1 = zsyn[0]; hdrs_Nfoc[l*nxim+j].f1 = zsyn[0];
hdrs_Nfoc[l*nxim+j].f2 = xsyn[0];
hdrs_Nfoc[l*nxim+j].d1 = zsyn[1]-zsyn[0]; hdrs_Nfoc[l*nxim+j].d1 = zsyn[1]-zsyn[0];
hdrs_Nfoc[l*nxim+j].d2 = xsyn[1]-xsyn[0];
hdrs_Nfoc[l*nxim+j].dt = (int)(hdrs_Nfoc[l*nxim+j].d1*(1E6)); hdrs_Nfoc[l*nxim+j].dt = (int)(hdrs_Nfoc[l*nxim+j].d1*(1E6));
hdrs_Nfoc[l*nxim+j].trwf = nxim*nyim; hdrs_Nfoc[l*nxim+j].trwf = nxim*nyim;
hdrs_Nfoc[l*nxim+j].ntr = nxim*nyim;
} }
} }
// Write the data // Write the data
...@@ -775,13 +786,24 @@ int main (int argc, char **argv) ...@@ -775,13 +786,24 @@ int main (int argc, char **argv)
for (l=0; l<nyim; l++){ for (l=0; l<nyim; l++){
for (j=0; j<nxim; j++){ for (j=0; j<nxim; j++){
hdrs_Nfoc[l*nxim+j].ns = nzim; hdrs_Nfoc[l*nxim+j].ns = nzim;
hdrs_Nfoc[l*nxim+j].sx = xsyn[j]; hdrs_Nfoc[l*nxim+j].fldr = 1;
hdrs_Nfoc[l*nxim+j].sy = ysyn[l]; hdrs_Nfoc[l*nxim+j].tracl = 1;
hdrs_Nfoc[l*nxim+j].sdepth = zsyn[l]; hdrs_Nfoc[l*nxim+j].tracf = l*nxim+j+1;
hdrs_Nfoc[l*nxim+j].trid = 1;
hdrs_Nfoc[l*nxim+j].scalco = -1000;
hdrs_Nfoc[l*nxim+j].scalel = -1000;
hdrs_Nfoc[l*nxim+j].sx = xsyn[j]*(1e3);
hdrs_Nfoc[l*nxim+j].sy = ysyn[l]*(1e3);
hdrs_Nfoc[l*nxim+j].gx = xsyn[j]*(1e3);
hdrs_Nfoc[l*nxim+j].gy = ysyn[l]*(1e3);
hdrs_Nfoc[l*nxim+j].sdepth = zsyn[l]*(1e3);
hdrs_Nfoc[l*nxim+j].f1 = zsyn[0]; hdrs_Nfoc[l*nxim+j].f1 = zsyn[0];
hdrs_Nfoc[l*nxim+j].f2 = xsyn[0];
hdrs_Nfoc[l*nxim+j].d1 = zsyn[1]-zsyn[0]; hdrs_Nfoc[l*nxim+j].d1 = zsyn[1]-zsyn[0];
hdrs_Nfoc[l*nxim+j].d2 = xsyn[1]-xsyn[0];
hdrs_Nfoc[l*nxim+j].dt = (int)(hdrs_Nfoc[l*nxim+j].d1*(1E6)); hdrs_Nfoc[l*nxim+j].dt = (int)(hdrs_Nfoc[l*nxim+j].d1*(1E6));
hdrs_Nfoc[l*nxim+j].trwf = nxim*nyim; hdrs_Nfoc[l*nxim+j].trwf = nxim*nyim;
hdrs_Nfoc[l*nxim+j].ntr = nxim*nyim;
} }
} }
// Write out image // Write out image
...@@ -946,4 +968,4 @@ long unique_elements(float *arr, long len) ...@@ -946,4 +968,4 @@ long unique_elements(float *arr, long len)
if (is_unique) ++unique; if (is_unique) ++unique;
} }
return unique; return unique;
} }
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment