Skip to content
Snippets Groups Projects
Commit 8c79c788 authored by Jan Thorbecke's avatar Jan Thorbecke
Browse files

before pull

parent 474f5c1d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
fdelmodc/fdelmodc
+ 0
− 0
View file @ 8c79c788
No preview for this file type
marchenko/marchenko.c
+ 20
− 14
View file @ 8c79c788
......@@ -56,7 +56,7 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
*reci_xsrc, int *reci_xrcv, float *ixmask, int verbose);
void synthesisPosistions(int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nfoc, float *xrcv, float *xsrc, int *xnx,
float fxse, float fxsb, float dxs, float dxsrc, float dx, int nshots, int *ixpos, int *npos, int *isxcount, int reci, int verbose);
float fxse, float fxsb, float dxs, float dxsrc, float dx, int nshots, int *ixpos, int *npos, int *isxcount, int countmin, int reci, int verbose);
int linearsearch(int *array, size_t N, int value);
......@@ -91,6 +91,7 @@ char *sdoc[] = {
" scale=2 .................. scale factor of R for summation of Ni with G_d",
" pad=0 .................... amount of samples to pad the reflection series",
" reci=0 ................... 1; add receivers as shots 2; only use receivers as shot positions",
" countmin=0 ............... 0.3*nxrcv; minumum number of reciprocal traces for a contribution",
" OUTPUT DEFINITION ",
" file_green= .............. output file with full Green function(s)",
" file_gplus= .............. output file with G+ ",
......@@ -116,7 +117,7 @@ int main (int argc, char **argv)
int i, j, l, ret, nshots, Nfoc, nt, nx, nts, nxs, ngath;
int size, n1, n2, ntap, tap, di, ntraces, pad;
int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
int reci, mode, n2out, verbose, ntfft;
int reci, countmin, mode, n2out, verbose, ntfft;
int iter, niter, tracf, *muteW;
int hw, smooth, above, shift, *ixpos, npos, ix;
int nshots_r, *isxcount, *reci_xsrc, *reci_xrcv;
......@@ -196,6 +197,7 @@ int main (int argc, char **argv)
nw_high = MIN((int)(fmax*ntfft*dt), nfreq-1);
nw = nw_high - nw_low + 1;
scl = 1.0/((float)ntfft);
if (!getparint("countmin", &countmin)) countmin = 0.3*nx;
/*================ Allocating all data arrays ================*/
......@@ -370,7 +372,7 @@ int main (int argc, char **argv)
/* dry-run of synthesis to get all x-positions calcalated by the integration */
synthesisPosistions(nx, nt, nxs, nts, dt, xsyn, Nfoc, xrcv, xsrc, xnx, fxse, fxsb,
dxs, dxsrc, dx, nshots, ixpos, &npos, isxcount, reci, verbose);
dxs, dxsrc, dx, nshots, ixpos, &npos, isxcount, countmin, reci, verbose);
if (verbose) {
vmess("synthesisPosistions: nshots=%d npos=%d", nshots, npos);
}
......@@ -442,12 +444,13 @@ int main (int argc, char **argv)
/* N_k(x,t) = -N_(k-1)(x,-t) */
/* p0^-(x,t) += iRN = (R * T_d^inv)(t) */
for (l = 0; l < Nfoc; l++) {
energyNi = 0.0;
for (i = 0; i < npos; i++) {
j = 0;
ix = ixpos[i];
Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+ix*nts+j];
pmin[l*nxs*nts+i*nts+j] += iRN[l*nxs*nts+ix*nts+j];
energyNi = iRN[l*nxs*nts+ix*nts+j]*iRN[l*nxs*nts+ix*nts+j];
energyNi += iRN[l*nxs*nts+ix*nts+j]*iRN[l*nxs*nts+ix*nts+j];
for (j = 1; j < nts; j++) {
Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+ix*nts+nts-j];
pmin[l*nxs*nts+i*nts+j] += iRN[l*nxs*nts+ix*nts+j];
......@@ -828,7 +831,7 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
#endif
} /* end of parallel region */
if (verbose>3) vmess("*** Shot gather %d processed ***", k);
if (verbose>4) vmess("*** Shot gather %d processed ***", k);
} /* end of nshots (k) loop */
} /* end of if reci
......@@ -862,14 +865,14 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
for (i = 0; i < inx; i++) {
il = reci_xrcv[ixsrc*nxs+i];
ik = reci_xsrc[ixsrc*nxs+i];
//ix = ixrcv[ik*nx+ik];
ix = NINT((xsrc[il] - fxsb)/dxs);
//if (j==nw_low) {
//fprintf(stderr,"ixsrc=%d with il=%d and ik=%d and ix=%d mask=%f\n", ixsrc, il, ik, ix, ixmask[ixsrc]);
// fprintf(stderr,"ixsrc=%d with il=%d and ik=%d and ix=%d mask=%f\n", ixsrc, il, ik, ix, ixmask[ixsrc]);
//}
sum[j].r += Refl[ik*nw*nx+m*nx+il].r*Fop[l*nw*nxs+m*nxs+ik].r -
Refl[ik*nw*nx+m*nx+il].i*Fop[l*nw*nxs+m*nxs+ik].i;
sum[j].i += Refl[ik*nw*nx+m*nx+il].i*Fop[l*nw*nxs+m*nxs+ik].r +
Refl[ik*nw*nx+m*nx+il].r*Fop[l*nw*nxs+m*nxs+ik].i;
sum[j].r += Refl[il*nw*nx+m*nx+ik].r*Fop[l*nw*nxs+m*nxs+ix].r -
Refl[il*nw*nx+m*nx+ik].i*Fop[l*nw*nxs+m*nxs+ix].i;
sum[j].i += Refl[il*nw*nx+m*nx+ik].i*Fop[l*nw*nxs+m*nxs+ix].r +
Refl[il*nw*nx+m*nx+ik].r*Fop[l*nw*nxs+m*nxs+ix].i;
}
}
......@@ -899,7 +902,7 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
}
void synthesisPosistions(int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nfoc, float *xrcv, float *xsrc, int *xnx,
float fxse, float fxsb, float dxs, float dxsrc, float dx, int nshots, int *ixpos, int *npos, int *isxcount, int reci, int verbose)
float fxse, float fxsb, float dxs, float dxsrc, float dx, int nshots, int *ixpos, int *npos, int *isxcount, int countmin, int reci, int verbose)
{
int i, j, l, ixsrc, ixrcv, dosrc, k, *count;
float x0, x1;
......@@ -951,7 +954,7 @@ float fxse, float fxsb, float dxs, float dxsrc, float dx, int nshots, int *ixpos
/* if reci=1 or reci=2 source-receive reciprocity is used and new (reciprocal-)sources are added */
if (reci != 0) {
for (k=0; k<nxs; k++) { /* check count in total number of shots added by reciprocity */
if (isxcount[k] != 0) {
if (isxcount[k] >= countmin) {
j = linearsearch(ixpos, *npos, k);
if (j < *npos) { /* the position (at j) is already included */
count[j] += isxcount[k];
......@@ -962,6 +965,9 @@ float fxse, float fxsb, float dxs, float dxsrc, float dx, int nshots, int *ixpos
*npos += 1;
}
}
else {
isxcount[k] = 0;
}
}
} /* end of reci branch */
} /* end of Nfoc loop */
......@@ -974,7 +980,7 @@ float fxse, float fxsb, float dxs, float dxsrc, float dx, int nshots, int *ixpos
free(count);
/* sort ixpos into increasing values */
qsort(ixpos, *npos, sizeof(int), compareInt);
// qsort(ixpos, *npos, sizeof(int), compareInt);
return;
......
marchenko/readShotData.c
+ 11
− 14
View file @ 8c79c788
......@@ -26,12 +26,10 @@ int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx
int *isx, *igx, k, l, m, j, nreci;
int samercv, samesrc, nxrk, nxrm, maxtraces, ixsrc;
float scl, scel, *trace, dt;
//float fxsb,dxs;
complex *ctrace;
/* Reading first header */
//fxsb = -2250; dxs=5;
if (filename == NULL) fp = stdin;
else fp = fopen( filename, "r" );
if ( fp == NULL ) {
......@@ -126,8 +124,7 @@ int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx
if ((sx_shot != hdr.sx) || (fldr_shot != hdr.fldr) ) break;
}
if (verbose>2) {
fprintf(stderr,"finished reading shot %d (%d) with %d traces\n",sx_shot,igath,itrace);
//disp_fileinfo(filename, nt, xnx[igath], hdr.f1, xrcv[igath*nxs], d1, d2, &hdr);
vmess("finished reading shot %d (%d) with %d traces",sx_shot,igath,itrace);
}
if (itrace != 0) { /* end of shot record */
......@@ -147,12 +144,13 @@ int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx
if (reci != 0) {
for (k=0; k<nshots; k++) {
ixsrc = NINT((xsrc[k] - fxsb)/dxs);
ixmask[ixsrc] = 1.0;
nxrk = xnx[k];
for (l=0; l<nxrk; l++) {
samercv = 0;
samesrc = 0;
for (m=0; m<nshots; m++) {
if (igx[k*nx+l] == isx[m] && reci == 1) { // receiver position already known as source
if (igx[k*nx+l] == isx[m] && reci == 1) { // receiver position already present as source position m
nxrm = xnx[m];
for (j=0; j<nxrm; j++) { // check if receiver l with source k is also present in shot m
if (isx[k] == igx[m*nx+j]) { // shot k with receiver l already known as receiver j in shot m: same data
......@@ -160,13 +158,13 @@ int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx
break;
}
}
if (samercv == 0) { // source k of receiver l -> accept trace as new receiver position for source m
ixsrc = NINT((xsrc[m] - fxsb)/dxs);
if (samercv == 0) { // source k of receiver l -> accept trace as new receiver position for source l
ixsrc = NINT((xrcv[k*nx+l] - fxsb)/dxs);
if ((ixsrc >= 0) && (ixsrc < nxs)) {
reci_xrcv[ixsrc*nxs+isxcount[ixsrc]] = l;
reci_xsrc[ixsrc*nxs+isxcount[ixsrc]] = k;
reci_xrcv[ixsrc*nxs+isxcount[ixsrc]] = k;
reci_xsrc[ixsrc*nxs+isxcount[ixsrc]] = l;
isxcount[ixsrc] += 1;
ixmask[ixsrc] = 0.5; // traces are added to already existing traces and must be scaled
if (reci==1) ixmask[ixsrc] = 0.5; // traces are added to already existing traces and must be scaled
}
}
samesrc = 1;
......@@ -177,10 +175,9 @@ int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx
//fprintf(stderr,"not a samesrc for receiver l=%d for source k=%d\n", l,k);
ixsrc = NINT((xrcv[k*nx+l] - fxsb)/dxs);
if ((ixsrc >= 0) && (ixsrc < nxs)) { // is this correct or should k and l be reversed: rcv=l src=k
reci_xrcv[ixsrc*nxs+isxcount[ixsrc]] = l;
reci_xsrc[ixsrc*nxs+isxcount[ixsrc]] = k;
reci_xrcv[ixsrc*nxs+isxcount[ixsrc]] = k;
reci_xsrc[ixsrc*nxs+isxcount[ixsrc]] = l;
isxcount[ixsrc] += 1;
ixmask[ixsrc] = 1.0;
}
}
}
......@@ -190,7 +187,7 @@ int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx
if (isxcount[k] != 0) {
maxtraces = MAX(maxtraces,isxcount[k]);
nreci++;
vmess("reciprocal receiver at %f (%d) has %d sources contributing", k, k*dxs+fxsb, isxcount[k]);
if (verbose>1) vmess("reciprocal receiver at %f (%d) has %d sources contributing", k, k*dxs+fxsb, isxcount[k]);
}
}
*nshots_r = nreci;
......
utils/getpars.c
+ 1
− 1
View file @ 8c79c788
......@@ -354,7 +354,7 @@ int getnpar (int n, char *name, char *type, void *ptr)
break;
case 'a':
{ char *tmpstr="";
tmpstr = (char *)calloc(strlen(aval),1);
tmpstr = (char *)calloc(strlen(aval)+1,1);
strchop(aval,tmpstr);
*(char**)ptr = tmpstr;
......
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