Skip to content
Snippets Groups Projects
Commit bb7b0ad9 authored by Joeri Brackenhoff's avatar Joeri Brackenhoff
Browse files

update

parent 2bae37b2
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 380 additions and 391 deletions
marchenko_applications/JespersRayTracer.c
+ 2
− 2
View file @ bb7b0ad9
......@@ -123,7 +123,7 @@ int getWaveParameter(float *slowness, icoord size, float dgrid, fcoord s, fcoord
if (ray.useT2 != 0)
T2 += getdT2(x, z, so, angle, ds, nRayTmp, s, r, rayReference3D, slowness, size);
if (ray.geomspread != 0) {
/*if (ray.geomspread != 0) {
if (so <= 0) {
dQdPhi = 0;
}
......@@ -131,7 +131,7 @@ int getWaveParameter(float *slowness, icoord size, float dgrid, fcoord s, fcoord
greentmp = greenIntP(lengthRefRay, so, lengthRefRay, slowness, size, nRayTmp, r, s);
dQdPhi += greentmp*secondDerivativeU1(slowness, size, x, z, angle, r, s)*ds/so;
}
}
}*/
}
if (ray.useT2)
......
marchenko_applications/Makefile
+ 26
− 8
View file @ bb7b0ad9
......@@ -7,7 +7,7 @@ LIBS += -L$L -lgenfft -lm $(LIBSM)
#ALL: fmute marchenko marchenko2
ALL: fmute marchenko
ALL: fmute marchenko_app combine
SRCJ = fmute.c \
getFileInfo.c \
......@@ -51,7 +51,19 @@ SRCH = marchenko.c \
gaussGen.c \
iterations.c \
imaging.c \
threadAffinity.c
threadAffinity.c \
makeWindow.c \
homogeneousg.c
SRCC = combine.c \
getFileInfo.c \
writeData.c \
wallclock_time.c \
getpars.c \
verbosepkg.c \
atopkge.c \
docpkge.c \
readSnapData.c
OBJJ = $(SRCJ:%.c=%.o)
......@@ -60,20 +72,26 @@ fmute: $(OBJJ)
OBJH = $(SRCH:%.c=%.o)
marchenko: $(OBJH) raytime.h
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko $(OBJH) $(LIBS)
marchenko_app: $(OBJH) raytime.h
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko_app $(OBJH) $(LIBS)
OBJC = $(SRCC:%.c=%.o)
combine: $(OBJC)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o combine $(OBJC) $(LIBS)
install: fmute marchenko
install: fmute marchenko_app combine
cp fmute $B
cp marchenko $B
cp marchenko_app $B
cp combine $B
# cp marchenko2 $B
clean:
rm -f core fmute $(OBJJ) marchenko $(OBJH)
rm -f core fmute $(OBJJ) marchenko_app $(OBJH) combine $(OBJC)
realclean: clean
rm -f $B/fmute $B/marchenko
rm -f $B/fmute $B/marchenko_app $B/combine
......
marchenko_applications/marchenko
+ 0
− 0
View file @ bb7b0ad9
No preview for this file type
marchenko_applications/marchenko.c
+ 193
− 119
View file @ bb7b0ad9
This diff is collapsed.
marchenko_applications/marchenko.h
+ 1
− 1
View file @ bb7b0ad9
......@@ -10,7 +10,7 @@
#define WAVEPAR
typedef struct WaveParameters {
int nt, shift, inv, scfft, cm, cn, wav;
float dt, fp, fmin, flef, frig, fmax, t0, db, scale, eps;
float dt, fp, fmin, flef, frig, fmax, t0, db, scale, eps, xloc, zloc;
char w[10], *file_wav;
} WavePar;
#endif
......
marchenko_applications/raytime.c
+ 10
− 1
View file @ bb7b0ad9
......@@ -42,7 +42,7 @@ int defineSource(wavPar wav, srcPar src, modPar mod, float **src_nwav, int rever
int writeSrcRecPos(modPar *mod, recPar *rec, srcPar *src, shotPar *shot);
int raytime(float *time, float *ampl, int *xnx, float *xrcv, float *xsrc, float *zsrc)
int raytime(float *time, float *ampl, int *xnx, float *xrcv, float *xsrc, float *zsrc, float xloc, float zloc)
{
modPar mod;
recPar rec;
......@@ -70,6 +70,15 @@ int raytime(float *time, float *ampl, int *xnx, float *xrcv, float *xsrc, float
if (!getparint("verbose",&verbose)) verbose=0;
getParameters(&mod, &rec, &src, &shot, &ray, verbose);
if (xloc!=-123456.0 && zloc!=-123456.0) {
if (verbose > 3) vmess("Setting source ray to x = %.3f, z = %.3f",xloc,zloc);
shot.nx = 1;
shot.nz = 1;
shot.n = 1;
shot.x[0] = NINT((xloc-mod.x0)/mod.dx);
shot.z[0] = NINT((zloc-mod.z0)/mod.dz);
}
/* allocate arrays for model parameters: the different schemes use different arrays */
n1 = mod.nz;
......
marchenko_applications/readTinvData.c
+ 135
− 247
View file @ bb7b0ad9
......@@ -25,276 +25,164 @@ void findShotInMute(float *xrcvMute, float xrcvShot, int nxs, int *imute);
int readSnapData(char *filename, float *data, segy *hdrs, int nsnaps, int nx, int nz, int sx, int ex, int sz, int ez);
int raytime(float *amp, float *time, int *xnx, float *xrcv, float *xsrc, float *zsrc);
int readTinvData(char *filename, WavePar WP, char *file_ray, char *file_amp, float dt, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nsyn, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose)
int readTinvData(char *filename, float dt, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nsyn, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose)
{
FILE *fp;
segy hdr, *hdrs_mute, *hdrs_amp, *hdrs_wav;
segy hdr, *hdrs_mute;
size_t nread;
char *file_cp;
int fldr_shot, sx_shot, itrace, one_shot, ig, isyn, i, j;
int end_of_file, nt, gx0, gx1, nfreq, geosp;
int end_of_file, nt, gx0, gx1;
int nx1, jmax, imax, tstart, tend, nwav;
float xmax, tmax, lmax, *wavelet, *wavelet2;
float scl, scel, *trace, dxrcv, *timeval, dw, *amp;
complex *cmute, *cwav;
float xmax, tmax, lmax;
float scl, scel, *trace, dxrcv;
if (!getparstring("file_cp", &file_cp)) file_cp=NULL;
if (!getparint("geomspread",&geosp)) geosp=1;
if (file_cp==NULL) geosp=0;
/*Check wheter the raytime is used or not*/
if (file_ray!=NULL || file_cp!=NULL) {
/*Define parameters*/
nfreq = ntfft/2+1;
wavelet = (float *)calloc(ntfft,sizeof(float));
cwav = (complex *)malloc(nfreq*sizeof(complex));
cmute = (complex *)malloc(nfreq*sizeof(complex));
dw = 2*M_PI/(ntfft*dt);
/* Reading first header */
/*Create wavelet using parameters or read in wavelet*/
if (WP.wav) {
if (WP.file_wav == NULL) {
if (verbose>0) vmess("Modeling wavelet");
freqwave(wavelet, WP.nt, WP.dt, WP.fp, WP.fmin, WP.flef, WP.frig, WP.fmax,
WP.t0, WP.db, WP.shift, WP.cm, WP.cn, WP.w, WP.scale, WP.scfft, WP.inv, WP.eps, verbose);
}
else {
if (verbose>0) vmess("Reading in wavelet");
fp = fopen( WP.file_wav, "r" );
if ( fp == NULL ) {
perror("Error opening file containing wavelet");
}
fclose(fp);
wavelet2= (float *)calloc(ntfft,sizeof(float));
hdrs_wav = (segy *)calloc(1, sizeof(segy));
readSnapData(WP.file_wav, wavelet2, hdrs_wav, Nsyn, 1, ntfft, 0, 1, 0, ntfft);
nwav = hdrs_wav[0].ns/2;
for (i=0; i<nwav; i++) {
wavelet[i] = wavelet2[i];
wavelet[ntfft-1-i] = wavelet2[hdrs_wav[0].ns-1-i];
}
}
rc1fft(wavelet,cwav,ntfft,-1);
free(wavelet);
}
timeval = (float *)calloc(Nsyn*nx,sizeof(float));
amp = (float *)calloc(Nsyn*nx,sizeof(float));
if (file_ray!=NULL) {
/* Defining mute window using raytimes */
vmess("Using raytime for mutewindow");
hdrs_mute = (segy *) calloc(Nsyn,sizeof(segy));
fp = fopen( file_ray, "r" );
if ( fp == NULL ) {
perror("Error opening file containing wavelet");
}
fclose(fp);
readSnapData(file_ray, timeval, hdrs_mute, Nsyn, 1, nx, 0, 1, 0, nx);
/*Check whether the amplitude is also used*/
if (file_amp != NULL) {
hdrs_amp = (segy *) calloc(Nsyn,sizeof(segy));
fp = fopen( file_amp, "r" );
if ( fp == NULL ) {
perror("Error opening file containing wavelet");
}
fclose(fp);
readSnapData(file_amp, amp, hdrs_amp, Nsyn, 1, nx, 0, 1, 0, nx);
}
/*Define source and receiver locations from the raytime*/
for (isyn=0; isyn<Nsyn; isyn++) {
for (itrace=0; itrace<nx; itrace++) {
xrcv[isyn*nx+itrace] = (hdrs_mute[isyn].f1 + hdrs_mute[isyn].d1*((float)itrace));
}
xnx[isyn]=nx;
xsrc[isyn] = hdrs_mute[isyn].sx;
zsrc[isyn] = hdrs_mute[isyn].sdepth;
}
}
else {
raytime(timeval,amp,xnx,xrcv,xsrc,zsrc);
}
if (filename == NULL) fp = stdin;
else fp = fopen( filename, "r" );
if ( fp == NULL ) {
fprintf(stderr,"input file %s has an error\n", filename);
perror("error in opening file: ");
fflush(stderr);
return -1;
}
/*Determine the mutewindow*/
for (j=0; j<Nsyn; j++) {
for (i=0; i<nx; i++) {
maxval[j*nx+i] = (int)roundf(timeval[j*nx+i]/dt);
if (maxval[j*nx+i] > ntfft-1) maxval[j*nx+i] = ntfft-1;
if (WP.wav) { /*Apply the wavelet to create a first arrival*/
if (file_amp != NULL || geosp==1) {
for (ig=0; ig<nfreq; ig++) {
cmute[ig].r = (cwav[ig].r*cos(ig*dw*timeval[j*nx+i]-M_PI/4.0)-cwav[ig].i*sin(ig*dw*timeval[j*nx+i]-M_PI/4.0))/amp[j*nx+i];
cmute[ig].i = (cwav[ig].i*cos(ig*dw*timeval[j*nx+i]-M_PI/4.0)+cwav[ig].r*sin(ig*dw*timeval[j*nx+i]-M_PI/4.0))/amp[j*nx+i];
}
}
else { /*Use the raytime only to determine the mutewindow*/
for (ig=0; ig<nfreq; ig++) {
cmute[ig].r = cwav[ig].r*cos(ig*dw*timeval[j*nx+i]-M_PI/4.0)-cwav[ig].i*sin(ig*dw*timeval[j*nx+i]-M_PI/4.0);
cmute[ig].i = cwav[ig].i*cos(ig*dw*timeval[j*nx+i]-M_PI/4.0)+cwav[ig].r*sin(ig*dw*timeval[j*nx+i]-M_PI/4.0);
}
}
cr1fft(cmute,&tinv[j*nx*ntfft+i*ntfft],ntfft,1);
tinv[j*nx*ntfft+i*ntfft] /= ntfft;
}
}
}
}
fseek(fp, 0, SEEK_SET);
nread = fread( &hdr, 1, TRCBYTES, fp );
assert(nread == TRCBYTES);
if (hdr.scalco < 0) scl = 1.0/fabs(hdr.scalco);
else if (hdr.scalco == 0) scl = 1.0;
else scl = hdr.scalco;
if (hdr.scalel < 0) scel = 1.0/fabs(hdr.scalel);
else if (hdr.scalel == 0) scel = 1.0;
else scel = hdr.scalel;
fseek(fp, 0, SEEK_SET);
if (WP.wav == 0 && file_ray==NULL && filename!=NULL) {
nt = hdr.ns;
trace = (float *)calloc(ntfft,sizeof(float));
/* Reading first header */
end_of_file = 0;
one_shot = 1;
isyn = 0;
if (filename == NULL) fp = stdin;
else fp = fopen( filename, "r" );
if ( fp == NULL ) {
fprintf(stderr,"input file %s has an error\n", filename);
perror("error in opening file: ");
fflush(stderr);
return -1;
}
/* Read shots in file */
fseek(fp, 0, SEEK_SET);
while (!end_of_file) {
/* start reading data (shot records) */
itrace = 0;
nread = fread( &hdr, 1, TRCBYTES, fp );
assert(nread == TRCBYTES);
if (hdr.scalco < 0) scl = 1.0/fabs(hdr.scalco);
else if (hdr.scalco == 0) scl = 1.0;
else scl = hdr.scalco;
if (hdr.scalel < 0) scel = 1.0/fabs(hdr.scalel);
else if (hdr.scalel == 0) scel = 1.0;
else scel = hdr.scalel;
fseek(fp, 0, SEEK_SET);
nt = hdr.ns;
trace = (float *)calloc(ntfft,sizeof(float));
end_of_file = 0;
one_shot = 1;
isyn = 0;
/* Read shots in file */
if (nread != TRCBYTES) { /* no more data in file */
break;
}
while (!end_of_file) {
/* start reading data (shot records) */
itrace = 0;
sx_shot = hdr.sx;
fldr_shot = hdr.fldr;
gx0 = hdr.gx;
xsrc[isyn] = sx_shot*scl;
zsrc[isyn] = hdr.selev*scel;
xnx[isyn] = 0;
ig = isyn*nx*ntfft;
while (one_shot) {
xrcv[isyn*nx+itrace] = hdr.gx*scl;
nread = fread( trace, sizeof(float), nt, fp );
assert (nread == hdr.ns);
/* copy trace to data array */
memcpy( &tinv[ig+itrace*ntfft], trace, nt*sizeof(float));
gx1 = hdr.gx;
itrace++;
/* read next hdr of next trace */
nread = fread( &hdr, 1, TRCBYTES, fp );
if (nread != TRCBYTES) { /* no more data in file */
if (nread != TRCBYTES) {
one_shot = 0;
end_of_file = 1;
break;
}
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,isyn,itrace);
//disp_fileinfo(filename, nt, xnx[isyn], hdr.f1, xrcv[isyn*nxm], d1, d2, &hdr);
}
sx_shot = hdr.sx;
fldr_shot = hdr.fldr;
gx0 = hdr.gx;
xsrc[isyn] = sx_shot*scl;
zsrc[isyn] = hdr.selev*scel;
xnx[isyn] = 0;
ig = isyn*nx*ntfft;
while (one_shot) {
xrcv[isyn*nx+itrace] = hdr.gx*scl;
nread = fread( trace, sizeof(float), nt, fp );
assert (nread == hdr.ns);
/* copy trace to data array */
memcpy( &tinv[ig+itrace*ntfft], trace, nt*sizeof(float));
gx1 = hdr.gx;
itrace++;
/* read next hdr of next trace */
nread = fread( &hdr, 1, TRCBYTES, fp );
if (nread != TRCBYTES) {
one_shot = 0;
end_of_file = 1;
break;
}
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,isyn,itrace);
//disp_fileinfo(filename, nt, xnx[isyn], hdr.f1, xrcv[isyn*nxm], d1, d2, &hdr);
}
/* look for maximum in shot record to define mute window */
/* find consistent (one event) maximum related to maximum value */
nx1 = itrace;
xnx[isyn]=nx1;
if (file_ray==NULL) {/*Use the raytime to determine the mutewindow instead of searching*/
/* alternative find maximum at source position */
dxrcv = (gx1 - gx0)*scl/(float)(nx1-1);
imax = NINT(((sx_shot-gx0)*scl)/dxrcv);
tmax=0.0;
jmax = 0;
for (j = 0; j < nt; j++) {
lmax = fabs(tinv[ig+imax*ntfft+j]);
if (lmax > tmax) {
jmax = j;
tmax = lmax;
if (lmax > xmax) {
xmax=lmax;
}
}
}
maxval[isyn*nx+imax] = jmax;
if (verbose >= 3) vmess("Mute max at src-trace %d is sample %d", imax, maxval[imax]);
/* search forward in trace direction from maximum in file */
for (i = imax+1; i < nx1; i++) {
tstart = MAX(0, (maxval[isyn*nx+i-1]-hw));
tend = MIN(nt-1, (maxval[isyn*nx+i-1]+hw));
jmax=tstart;
tmax=0.0;
for(j = tstart; j <= tend; j++) {
lmax = fabs(tinv[ig+i*ntfft+j]);
if (lmax > tmax) {
jmax = j;
tmax = lmax;
}
}
maxval[isyn*nx+i] = jmax;
}
/* search backward in trace direction from maximum in file */
for (i = imax-1; i >=0; i--) {
tstart = MAX(0, (maxval[isyn*nx+i+1]-hw));
tend = MIN(nt-1, (maxval[isyn*nx+i+1]+hw));
jmax=tstart;
tmax=0.0;
for (j = tstart; j <= tend; j++) {
lmax = fabs(tinv[ig+i*ntfft+j]);
if (lmax > tmax) {
jmax = j;
tmax = lmax;
}
}
maxval[isyn*nx+i] = jmax;
}
}
if (itrace != 0) { /* end of shot record, but not end-of-file */
fseek( fp, -TRCBYTES, SEEK_CUR );
isyn++;
}
else {
end_of_file = 1;
}
/* look for maximum in shot record to define mute window */
/* find consistent (one event) maximum related to maximum value */
nx1 = itrace;
xnx[isyn]=nx1;
/* alternative find maximum at source position */
dxrcv = (gx1 - gx0)*scl/(float)(nx1-1);
imax = NINT(((sx_shot-gx0)*scl)/dxrcv);
tmax=0.0;
jmax = 0;
for (j = 0; j < nt; j++) {
lmax = fabs(tinv[ig+imax*ntfft+j]);
if (lmax > tmax) {
jmax = j;
tmax = lmax;
if (lmax > xmax) {
xmax=lmax;
}
}
}
maxval[isyn*nx+imax] = jmax;
if (verbose >= 3) vmess("Mute max at src-trace %d is sample %d", imax, maxval[imax]);
/* search forward in trace direction from maximum in file */
for (i = imax+1; i < nx1; i++) {
tstart = MAX(0, (maxval[isyn*nx+i-1]-hw));
tend = MIN(nt-1, (maxval[isyn*nx+i-1]+hw));
jmax=tstart;
tmax=0.0;
for(j = tstart; j <= tend; j++) {
lmax = fabs(tinv[ig+i*ntfft+j]);
if (lmax > tmax) {
jmax = j;
tmax = lmax;
}
}
maxval[isyn*nx+i] = jmax;
}
/* search backward in trace direction from maximum in file */
for (i = imax-1; i >=0; i--) {
tstart = MAX(0, (maxval[isyn*nx+i+1]-hw));
tend = MIN(nt-1, (maxval[isyn*nx+i+1]+hw));
jmax=tstart;
tmax=0.0;
for (j = tstart; j <= tend; j++) {
lmax = fabs(tinv[ig+i*ntfft+j]);
if (lmax > tmax) {
jmax = j;
tmax = lmax;
}
}
maxval[isyn*nx+i] = jmax;
}
if (itrace != 0) { /* end of shot record, but not end-of-file */
fseek( fp, -TRCBYTES, SEEK_CUR );
isyn++;
}
else {
end_of_file = 1;
}
/* copy trace to data array for mode=-1 */
/* time reverse trace */
if (mode==-1) {
for (i = 0; i < nx1; i++) {
memcpy( trace, &tinv[ig+i*ntfft], ntfft*sizeof(float));
j=0;
tinv[ig+i*ntfft+j] = trace[j];
for (j=1; j<ntfft; j++) tinv[ig+i*ntfft+ntfft-j] = trace[j];
}
/* copy trace to data array for mode=-1 */
/* time reverse trace */
if (mode==-1) {
for (i = 0; i < nx1; i++) {
memcpy( trace, &tinv[ig+i*ntfft], ntfft*sizeof(float));
j=0;
tinv[ig+i*ntfft+j] = trace[j];
for (j=1; j<ntfft; j++) tinv[ig+i*ntfft+ntfft-j] = trace[j];
}
}
free(trace);
}
free(trace);
return 0;
}
......
marchenko_full/Makefile
+ 7
− 7
View file @ bb7b0ad9
......@@ -7,7 +7,7 @@ OPTC += -g -O0 -Wall
#ALL: fmute marchenko marchenko2
ALL: fmute marchenko
ALL: fmute marchenko_full
SRCJ = fmute.c \
getFileInfo.c \
......@@ -58,20 +58,20 @@ fmute: $(OBJJ)
OBJH = $(SRCH:%.c=%.o)
marchenko: $(OBJH) raytime.h
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko $(OBJH) $(LIBS)
marchenko_full: $(OBJH) raytime.h
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko_full $(OBJH) $(LIBS)
install: fmute marchenko
install: fmute marchenko_full
cp fmute $B
cp marchenko $B
cp marchenko_full $B
# cp marchenko2 $B
clean:
rm -f core fmute $(OBJJ) marchenko $(OBJH)
rm -f core fmute $(OBJJ) marchenko_full $(OBJH)
realclean: clean
rm -f $B/fmute $B/marchenko
rm -f $B/fmute $B/marchenko_full
......
marchenko_full/fmute
+ 0
− 0
View file @ bb7b0ad9
No preview for this file type
marchenko_full/marchenko deleted 100755 → 0
+ 0
− 0
View file @ 2bae37b2
File deleted
marchenko_full/marchenko.c
+ 6
− 6
View file @ bb7b0ad9
......@@ -49,6 +49,7 @@ int readData(FILE *fp, float *data, segy *hdrs, int n1);
int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
double wallclock_time(void);
int getParameters(modPar *mod, recPar *rec, srcPar *src, shotPar *shot, rayPar *ray, int verbose);
void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, double *tfft, int verbose);
......@@ -214,7 +215,7 @@ int main (int argc, char **argv)
double t0, t1, t2, t3, tsyn, tread, tfft, tcopy, energyNi, *J;
float d1, d2, f1, f2, fxs, ft, fx, *xsyn, dxsrc, Q, f0, *Costdet;
float *green, *f2p, *pmin, *G_d, dt, dx, dxs, scl, mem;
float *f1plus, *f1min, *iRN, *Ni, *trace, *Gmin, *Gplus, *Gm0;
float *f1plus, *f1min, *iRN, *Ni, *Gmin, *Gplus, *Gm0;
float xmin, xmax, weight, tsq, *Gd, *amp, bstart, bend, db, *bdet, bp, b, bmin;
complex *Refl, *Fop;
char *file_tinv, *file_shot, *file_green, *file_iter, *file_wav, *file_ray, *file_amp;
......@@ -367,7 +368,6 @@ int main (int argc, char **argv)
Ni = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
G_d = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
muteW = (int *)calloc(Nsyn*nxs,sizeof(int));
trace = (float *)malloc(ntfft*sizeof(float));
ixpossyn = (int *)malloc(nxs*sizeof(int));
xrcvsyn = (float *)calloc(Nsyn*nxs,sizeof(float));
xsyn = (float *)malloc(Nsyn*sizeof(float));
......@@ -978,7 +978,7 @@ for (ib=0; ib<=nb; ib++) {
}
if (file_f1plus != NULL) {
/* rotate to get t=0 in the middle */
for (i = 0; i < n2; i++) {
/*for (i = 0; i < n2; i++) {
hdrs_out[i].f1 = -n1*0.5*dt;
memcpy(&trace[0],&f1plus[l*size+i*nts],nts*sizeof(float));
for (j = 0; j < n1/2; j++) {
......@@ -987,13 +987,13 @@ for (ib=0; ib<=nb; ib++) {
for (j = n1/2; j < n1; j++) {
f1plus[l*size+i*nts+j-n1/2] = trace[j];
}
}
}*/
ret = writeData(fp_f1plus, (float *)&f1plus[l*size], hdrs_out, n1, n2);
if (ret < 0 ) verr("error on writing output file.");
}
if (file_f1min != NULL) {
/* rotate to get t=0 in the middle */
for (i = 0; i < n2; i++) {
/*for (i = 0; i < n2; i++) {
hdrs_out[i].f1 = -n1*0.5*dt;
memcpy(&trace[0],&f1min[l*size+i*nts],nts*sizeof(float));
for (j = 0; j < n1/2; j++) {
......@@ -1002,7 +1002,7 @@ for (ib=0; ib<=nb; ib++) {
for (j = n1/2; j < n1; j++) {
f1min[l*size+i*nts+j-n1/2] = trace[j];
}
}
}*/
ret = writeData(fp_f1min, (float *)&f1min[l*size], hdrs_out, n1, n2);
if (ret < 0 ) verr("error on writing output 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