From f4ee2b84663149327fed244a61bb621d8f62a528 Mon Sep 17 00:00:00 2001
From: JanThorbecke <janth@xs4all.nl>
Date: Thu, 4 Jun 2015 07:36:50 +0200
Subject: [PATCH] added receievers on elipse and bug fixes
---
fdelmodc/decomposition.c | 2 +-
fdelmodc/fdelmodc.h | 1 +
fdelmodc/getParameters.c | 13 +++++++------
fdelmodc/getRecTimes.c | 1 +
fdelmodc/recvPar.c | 24 ++++++++++++++----------
fdelmodc/writeSrcRecPos.c | 7 +++++++
utils/marchenko.c | 36 ++++++++++++++++++------------------
7 files changed, 49 insertions(+), 35 deletions(-)
diff --git a/fdelmodc/decomposition.c b/fdelmodc/decomposition.c
index e9c8a0a..fccee65 100644
--- a/fdelmodc/decomposition.c
+++ b/fdelmodc/decomposition.c
@@ -187,7 +187,7 @@ void decud(float om, float rho, float cp, float dx, int nkx, float kangle, float
kxnyq = M_PI/dx;
if (kpos > kxnyq) kpos = kxnyq;
band = kpos;
- filterpoints = (int)fabs((int)(perc*band/dkx));
+ filterpoints = (int)abs((int)(perc*band/dkx));
kfilt = fabs(dkx*filterpoints);
if (kpos+kfilt < kxnyq) {
kxfmax = kpos+kfilt;
diff --git a/fdelmodc/fdelmodc.h b/fdelmodc/fdelmodc.h
index 7be9c50..6ed46c5 100644
--- a/fdelmodc/fdelmodc.h
+++ b/fdelmodc/fdelmodc.h
@@ -21,6 +21,7 @@ typedef struct _receiverPar { /* Receiver Parameters */
int nt;
int delay;
int skipdt;
+ int max_nrec;
int *z;
int *x;
float *zr;
diff --git a/fdelmodc/getParameters.c b/fdelmodc/getParameters.c
index 41dd704..79fabb6 100644
--- a/fdelmodc/getParameters.c
+++ b/fdelmodc/getParameters.c
@@ -53,7 +53,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
float *xsrca, *zsrca, rrcv;
float rsrc, oxsrc, ozsrc, dphisrc, ncsrc;
size_t nsamp;
- int i, j, max_nrec;
+ int i, j;
int boundary, ibnd, cfree;
int ntaper,tapleft,tapright,taptop,tapbottom;
int nxsrc, nzsrc;
@@ -1017,17 +1017,18 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
if (!getparint("rec_int_p",&rec->int_p)) rec->int_p=0;
if (!getparint("rec_int_vx",&rec->int_vx)) rec->int_vx=0;
if (!getparint("rec_int_vz",&rec->int_vz)) rec->int_vz=0;
- if (!getparint("max_nrec",&max_nrec)) max_nrec=15000;
+ if (!getparint("max_nrec",&rec->max_nrec)) rec->max_nrec=15000;
if (!getparint("scale",&rec->scale)) rec->scale=0;
if (!getparfloat("dxspread",&dxspread)) dxspread=0;
if (!getparfloat("dzspread",&dzspread)) dzspread=0;
rec->nt=MIN(rec->nt, NINT((mod->tmod-rdelay)/dtrcv)+1);
rec->delay=NINT(rdelay/mod->dt);
- rec->x = (int *)malloc(max_nrec*sizeof(int));
- rec->z = (int *)malloc(max_nrec*sizeof(int));
- rec->xr = (float *)malloc(max_nrec*sizeof(float));
- rec->zr = (float *)malloc(max_nrec*sizeof(float));
+ rec->max_nrec += rec->max_nrec+1;
+ rec->x = (int *)malloc(rec->max_nrec*sizeof(int));
+ rec->z = (int *)malloc(rec->max_nrec*sizeof(int));
+ rec->xr = (float *)malloc(rec->max_nrec*sizeof(float));
+ rec->zr = (float *)malloc(rec->max_nrec*sizeof(float));
/* calculates the receiver coordinates */
diff --git a/fdelmodc/getRecTimes.c b/fdelmodc/getRecTimes.c
index 4c784d5..b1da156 100644
--- a/fdelmodc/getRecTimes.c
+++ b/fdelmodc/getRecTimes.c
@@ -3,6 +3,7 @@
#include<math.h>
#include<assert.h>
#include"fdelmodc.h"
+#include"par.h"
/**
* Stores the wavefield at the receiver positions.
diff --git a/fdelmodc/recvPar.c b/fdelmodc/recvPar.c
index 0cfb7df..1bf638a 100644
--- a/fdelmodc/recvPar.c
+++ b/fdelmodc/recvPar.c
@@ -29,22 +29,23 @@ int recvPar(recPar *rec, float sub_x0, float sub_z0, float dx, float dz, int nx,
double circ, h, a, b, e, s, xr, zr, dr, srun, phase;
float xrange, zrange;
int Nx1, Nx2, Nz1, Nz2, Ndx, Ndz, iarray, nskip, nrec, nh;
- int nxrcv, nzrcv, ncrcv, nrcv, max_nrec;
+ int nxrcv, nzrcv, ncrcv, nrcv;
float *xrcva, *zrcva;
if(!getparint("verbose", &verbose)) verbose = 0;
- if (!getparint("max_nrec",&max_nrec)) max_nrec=15000;
-
nrec=0;
/* check if receiver positions on a circle are defined */
if (getparfloat("rrcv", &rrcv)) {
if (!getparfloat("dphi",&dphi)) dphi=2.0;
+ ncrcv = NINT(360.0/dphi);
+ if (ncrcv >= rec->max_nrec) {
+ verr("Number of receivers is larger than max: increase parameter max_nrec= %d to %d",rec->max_nrec, ncrcv);
+ }
if (!getparfloat("oxrcv",&oxrcv)) oxrcv=0.0;
if (!getparfloat("ozrcv",&ozrcv)) ozrcv=0.0;
if (!getparfloat("arcv",&arcv)) {
arcv=rrcv;
- ncrcv = NINT(360.0/dphi);
for (ix=0; ix<ncrcv; ix++) {
rec->xr[ix] = oxrcv-sub_x0+rrcv*cos(((ix*dphi)/360.0)*(2.0*M_PI));
rec->zr[ix] = ozrcv-sub_z0+arcv*sin(((ix*dphi)/360.0)*(2.0*M_PI));
@@ -57,7 +58,6 @@ int recvPar(recPar *rec, float sub_x0, float sub_z0, float dx, float dz, int nx,
}
else { /* an ellipse */
/* simple numerical solution to find equidistant points on an ellipse */
- ncrcv = NINT(360.0/dphi);
nh = ncrcv*1000; /* should be fine enough for most configurations */
h = 2.0*M_PI/nh;
a = MAX(rrcv, arcv);
@@ -108,6 +108,9 @@ int recvPar(recPar *rec, float sub_x0, float sub_z0, float dx, float dz, int nx,
verr("Number of receivers in array xrcva (%d), zrcva(%d) are not equal",nxrcv, nzrcv);
}
+ if ((nrec+nxrcv) >= rec->max_nrec) {
+ verr("Number of receivers is larger than max: increase parameter max_nrec= %d to %d",rec->max_nrec, nrec+nxrcv);
+ }
if (nxrcv != 0) {
/* receiver array is defined */
xrcva = (float *)malloc(nxrcv*sizeof(float));
@@ -140,6 +143,7 @@ int recvPar(recPar *rec, float sub_x0, float sub_z0, float dx, float dz, int nx,
Nx1=1;
}
+
if (Nx1!=0) {
xrcv1 = (float *)malloc(Nx1*sizeof(float));
xrcv2 = (float *)malloc(Nx1*sizeof(float));
@@ -206,7 +210,7 @@ int recvPar(recPar *rec, float sub_x0, float sub_z0, float dx, float dz, int nx,
xrange = (xrcv2[iarray]-xrcv1[iarray]);
zrange = (zrcv2[iarray]-zrcv1[iarray]);
if (dxr[iarray] != 0.0) {
- nrcv = NINT(abs(xrange/dxr[iarray]))+1;
+ nrcv = NINT(fabs(xrange/dxr[iarray]))+1;
dxrcv=dxr[iarray];
dzrcv = zrange/(nrcv-1);
if (dzrcv != dzr[iarray]) {
@@ -218,7 +222,7 @@ int recvPar(recPar *rec, float sub_x0, float sub_z0, float dx, float dz, int nx,
if (dzr[iarray] == 0) {
verr("For receiver array %d: receiver distance dzrcv is not given", iarray);
}
- nrcv=NINT(abs(zrange/dzr[iarray]))+1;
+ nrcv=NINT(fabs(zrange/dzr[iarray]))+1;
dxrcv = xrange/(nrcv-1);
dzrcv = dzr[iarray];
if (dxrcv != dxr[iarray]) {
@@ -235,9 +239,9 @@ int recvPar(recPar *rec, float sub_x0, float sub_z0, float dx, float dz, int nx,
rec->x[nrec]=NINT((rec->xr[nrec])/dx);
rec->z[nrec]=NINT((rec->zr[nrec])/dz);
nrec++;
- }
- if (nrec >= max_nrec) {
- verr("Number of receivers in arrays xrcv1,xrcv2,zrcv1,zrcv2 are larger than max: increas max_nrec %d",max_nrec);
+ if (nrec > rec->max_nrec) {
+ verr("Number of receivers is larger than max: increase parameter max_nrec=%d to %d",rec->max_nrec, nrec+nxrcv);
+ }
}
}
diff --git a/fdelmodc/writeSrcRecPos.c b/fdelmodc/writeSrcRecPos.c
index 0d4886e..c9cec66 100644
--- a/fdelmodc/writeSrcRecPos.c
+++ b/fdelmodc/writeSrcRecPos.c
@@ -36,6 +36,11 @@ int writeSrcRecPos(modPar *mod, recPar *rec, srcPar *src, shotPar *shot)
sub_x0 = mod->x0;
sub_z0 = mod->z0;
+// ibndx = mod.ioPx;
+// ibndz = mod.ioPz;
+// if (bnd.lef==4 || bnd.lef==2) ibndx += bnd.ntap;
+// if (bnd.top==4 || bnd.top==2) ibndz += bnd.ntap;
+
/* write velocity field with positions of the sources */
dum = (float *)calloc(nx*nz, sizeof(float));
vmess("Positions: shot=%d src=%d rec=%d", shot->n, src->n, rec->n);
@@ -119,6 +124,8 @@ int writeSrcRecPos(modPar *mod, recPar *rec, srcPar *src, shotPar *shot)
dum[(MIN(nx-1,rec->x[is]+1))*nz+rec->z[is]] = -1.0;
dum[rec->x[is]*nz+MAX(0,rec->z[is]-1)] = -1.0;
dum[rec->x[is]*nz+MIN(nz-1,rec->z[is]+1)] = -1.0;
+
+// vmess("receiver position %d at grid[ix=%d, iz=%d] = (x=%f z=%f)", ir, ix+ioPx, rec.z[ir]+ioPz, rec.xr[ir]+mod.x0, rec.zr[ir]+mod.z0);
if (rec->int_vx==3) {
fprintf(fp, "(%f, %f)\n", rec->xr[is]*dx+sub_x0, rec->zr[is]*dz+sub_z0);
}
diff --git a/utils/marchenko.c b/utils/marchenko.c
index 19da0aa..5b53e49 100644
--- a/utils/marchenko.c
+++ b/utils/marchenko.c
@@ -131,7 +131,7 @@ int main (int argc, char **argv)
if (verbose) vwarn("parameter file_green not found, assume pipe");
file_green = NULL;
}
- if (!getparint("nshots", &nshots)) verr("nshots should be given");
+ //if (!getparint("nshots", &nshots)) verr("nshots should be given");
if (!getparfloat("fmin", &fmin)) fmin = 0.0;
if (!getparfloat("fmax", &fmax)) fmax = 70.0;
if (!getparint("ixa", &ixa)) ixa = 0;
@@ -308,7 +308,7 @@ int main (int argc, char **argv)
vmess("direction of increasing traces = %d", di);
vmess("number of time samples(fft) s = %d %d", nt, optn);
vmess("time sampling = %e ", dt);
- if (file_greeb != NULL) vmess("Green output file = %s ", file_gmin);
+ if (file_green != NULL) vmess("Green output file = %s ", file_green);
if (file_gmin != NULL) vmess("Gmin output file = %s ", file_gmin);
if (file_gplus != NULL) vmess("Gplus output file = %s ", file_gplus);
if (file_pmin != NULL) vmess("Pmin output file = %s ", file_pmin);
@@ -656,7 +656,7 @@ int main (int argc, char **argv)
fp_gplus = fopen(file_gplus, "w+");
if (fp_gplus==NULL) verr("error on creating output file %s", file_gplus);
}
- if (file_pplus!=NULL) {
+ if (file_pplus != NULL) {
fp_pplus = fopen(file_pplus, "w+");
if (fp_pplus==NULL) verr("error on creating output file %s", file_pplus);
}
@@ -664,11 +664,11 @@ int main (int argc, char **argv)
fp_pmin = fopen(file_pmin, "w+");
if (fp_pmin==NULL) verr("error on creating output file %s", file_pmin);
}
- if (file_f1plus!=NULL) {
+ if (file_f1plus != NULL) {
fp_f1plus = fopen(file_f1plus, "w+");
if (fp_f1plus==NULL) verr("error on creating output file %s", file_f1plus);
}
- if (file_f1min!=NULL) {
+ if (file_f1min != NULL) {
fp_f1min = fopen(file_f1min, "w+");
if (fp_f1min==NULL) verr("error on creating output file %s", file_f1min);
}
@@ -707,23 +707,23 @@ int main (int argc, char **argv)
ret = writeData(fp_out, (float *)&green[l*size], hdrs_out, n1, n2);
if (ret < 0 ) verr("error on writing output file.");
- if (fp_gmin != NULL) {
+ if (file_gmin != NULL) {
ret = writeData(fp_gmin, (float *)&Gmin[l*size], hdrs_out, n1, n2);
if (ret < 0 ) verr("error on writing output file.");
}
- if (fp_gplus != NULL) {
+ if (file_gplus != NULL) {
ret = writeData(fp_gplus, (float *)&Gplus[l*size], hdrs_out, n1, n2);
if (ret < 0 ) verr("error on writing output file.");
}
- if (fp_pplus != NULL) {
+ if (file_pplus != NULL) {
ret = writeData(fp_pplus, (float *)&pplus[l*size], hdrs_out, n1, n2);
if (ret < 0 ) verr("error on writing output file.");
}
- if (fp_pmin != NULL) {
+ if (file_pmin != NULL) {
ret = writeData(fp_pmin, (float *)&pmin[l*size], hdrs_out, n1, n2);
if (ret < 0 ) verr("error on writing output file.");
}
- if (fp_f1plus != NULL) {
+ if (file_f1plus != NULL) {
/* rotate to get t=0 in the middle */
for (i = 0; i < n2; i++) {
hdrs_out[i].f1 = -n1*0.5*dt;
@@ -738,7 +738,7 @@ int main (int argc, char **argv)
ret = writeData(fp_f1plus, (float *)&f1plus[l*size], hdrs_out, n1, n2);
if (ret < 0 ) verr("error on writing output file.");
}
- if (fp_f1min != NULL) {
+ if (file_f1min != NULL) {
/* rotate to get t=0 in the middle */
for (i = 0; i < n2; i++) {
hdrs_out[i].f1 = -n1*0.5*dt;
@@ -755,12 +755,12 @@ int main (int argc, char **argv)
}
}
ret = fclose(fp_out);
- if (fp_gplus != NULL) ret += fclose(fp_gplus);
- if (fp_gmin != NULL) ret += fclose(fp_gmin);
- if (fp_pplus != NULL) ret += fclose(fp_pplus);
- if (fp_pmin != NULL) ret += fclose(fp_pmin);
- if (fp_f1plus != NULL) ret += fclose(fp_f1plus);
- if (fp_f1min != NULL) ret += fclose(fp_f1min);
+ if (file_gplus != NULL) {ret += fclose(fp_gplus);}
+ if (file_gmin != NULL) {ret += fclose(fp_gmin);}
+ if (file_pplus != NULL) {ret += fclose(fp_pplus);}
+ if (file_pmin != NULL) {ret += fclose(fp_pmin);}
+ if (file_f1plus != NULL) {ret += fclose(fp_f1plus);}
+ if (file_f1min != NULL) {ret += fclose(fp_f1min);}
if (ret < 0) verr("err %d on closing output file",ret);
if (verbose) {
@@ -814,7 +814,7 @@ void synthesis(complex *shots, complex *syncdata, float *syndata, int nx, int nt
ixsrc = NINT((xsrc[k] - fxs)/dxs);
- if (abs(xrcv[k*nx+0]-xsrc[k]) > 0.5*nx*dx) { iox = 0; inx = nx-off; }
+ if (fabs(xrcv[k*nx+0]-xsrc[k]) > 0.5*nx*dx) { iox = 0; inx = nx-off; }
else { iox = off; inx = nx; }
/*================ SYNTHESIS ================*/
--
GitLab