From 0ca42abeede95d2d9db2aba824e366d628baab99 Mon Sep 17 00:00:00 2001 From: Niklas Thiel Date: Wed, 9 Mar 2016 10:37:48 +0100 Subject: [PATCH] fixed wrong calculation of filter + minor changes --- src/IFOS2D.c | 43 ++++++++++++++++++++++++++++-------- src/inseis.c | 7 +++--- src/pup.c | 54 ++++++++++++++++++++------------------------- src/read_par_json.c | 10 ++++----- 4 files changed, 67 insertions(+), 47 deletions(-) diff --git a/src/IFOS2D.c b/src/IFOS2D.c index d188097..a4e5539 100644 --- a/src/IFOS2D.c +++ b/src/IFOS2D.c @@ -688,17 +688,21 @@ int main(int argc, char **argv){ fulldata_curl = matrix(1,ntr_glob,1,NT); break; case 5 : /* everything except curl and div*/ - if (WAVETYPE==1) { + switch (WAVETYPE) { + case 1: fulldata_vx = matrix(1,ntr_glob,1,NT); fulldata_vy = matrix(1,ntr_glob,1,NT); - } - if (WAVETYPE==2) { + break; + + case 2: fulldata_vz = matrix(1,ntr_glob,1,NT); - } - if (WAVETYPE==3) { + break; + + case 3: fulldata_vx = matrix(1,ntr_glob,1,NT); fulldata_vy = matrix(1,ntr_glob,1,NT); fulldata_vz = matrix(1,ntr_glob,1,NT); + break; } fulldata_p = matrix(1,ntr_glob,1,NT); break; @@ -1433,7 +1437,14 @@ int main(int argc, char **argv){ if (MYID == 0) fprintf(FP,"\n start wavefield separation..."); if (MYID == 0) pup(fulldata_p, fulldata_vy, FP, ntr_glob,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1); if (MYID == 0) fprintf(FP,"\n wavefield separation finished"); - inseis(fprec,ishot,sectionp,ntr_glob,ns,3,iter); + inseis(fprec,ishot,sectionread,ntr_glob,ns,3,iter); + h=1; + for(i=1;i<=ntr;i++){ + for(j=1;j<=ns;j++){ + sectionp[h][j]=sectionread[recpos_loc[3][i]][j]; + } + h++; + } } break; @@ -1960,8 +1971,15 @@ int main(int argc, char **argv){ if (WAVESEP==1) { if (MYID == 0) fprintf(FP,"\n start wavefield separation..."); if (MYID == 0) pup(fulldata_p, fulldata_vy, FP, ntr_glob,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1); - if (MYID == 0) fprintf(FP,"\n wavefield separation finished"); - inseis(fprec,ishot,sectionp,ntr_glob,ns,3,iter); + if (MYID == 0) fprintf(FP,"wavefield separation finished."); + inseis(fprec,ishot,sectionread,ntr_glob,ns,3,iter); + h=1; + for(i=1;i<=ntr;i++){ + for(j=1;j<=ns;j++){ + sectionp[h][j]=sectionread[recpos_loc[3][i]][j]; + } + h++; + } } break; @@ -3468,7 +3486,14 @@ int main(int argc, char **argv){ if (MYID == 0) fprintf(FP,"\n start wavefield separation..."); if (MYID == 0) pup(fulldata_p, fulldata_vy, FP, ntr_glob,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,2); if (MYID == 0) fprintf(FP,"\n wavefield separation finished"); - inseis(fprec,ishot,sectionp,ntr_glob,ns,3,iter); + inseis(fprec,ishot,sectionread,ntr_glob,ns,3,iter); + h=1; + for(i=1;i<=ntr;i++){ + for(j=1;j<=ns;j++){ + sectionp[h][j]=sectionread[recpos_loc[3][i]][j]; + } + h++; + } } /*-------------------------------------------------------------------------------*/ diff --git a/src/inseis.c b/src/inseis.c index bc1bc27..a6ddffc 100644 --- a/src/inseis.c +++ b/src/inseis.c @@ -17,7 +17,7 @@ -----------------------------------------------------------------------------------------*/ /*------------------------------------------------------------------------ - * Write seismograms to disk + * Read seismograms from disk * ----------------------------------------------------------------------*/ #include "fd.h" #include "segy.h" @@ -27,6 +27,7 @@ void inseis(FILE *fp, int comp, float **section, int ntr, int ns, int sws, int /* declaration of extern variables */ extern int NDT, MYID; extern char DATA_DIR[STRING_SIZE]; + extern char SEIS_FILE[STRING_SIZE]; extern float TIME, DH, DT, REFREC[4]; char data[STRING_SIZE]; FILE *fpdata; @@ -40,11 +41,11 @@ void inseis(FILE *fp, int comp, float **section, int ntr, int ns, int sws, int } if(sws==3){ /* open seismic data pup */ - sprintf(data,"%s_pup.su.shot%i.it%i",DATA_DIR,comp,iter); + sprintf(data,"%s_pup.su.shot%i.it%i",SEIS_FILE,comp,iter); } if(sws==4){ /* open seismic data pup used for step length calculation */ - sprintf(data,"%s_pup.su.step.shot%i.it%i",DATA_DIR,comp,iter); + sprintf(data,"%s_pup.su.step.shot%i.it%i",SEIS_FILE,comp,iter); } /* sws 5 -- 6 not used */ diff --git a/src/pup.c b/src/pup.c index 33418a3..86bf7d0 100644 --- a/src/pup.c +++ b/src/pup.c @@ -28,35 +28,36 @@ #include "fd.h" void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos, int **recpos_loc,int ntr, float ** srcpos, int ishot, int ns, int iter, int sw ) { - extern float DT, DH; + extern float DT,DH; extern float VEL, DENS, ANGTAPI, ANGTAPO; extern char SEIS_FILE[STRING_SIZE]; extern int SEIS_FORMAT; extern FILE *FP; - int i, j, k, nx_pot, ny_pot,nxy; - float nyq_k, nyq_f, dk, df, angi, ango, ka; + int i, j, k, nx_pot, ny_pot; + float nyq_k, nyq_f, dk, df, angi, ango, ka, dtr; float *w, *kx; float **filtwk, **filtwk_full; float **data_pim,**data_vyim,**data_pup,**data_pupim,**data_pt,**data_vyt; char pf[STRING_SIZE], pf1[STRING_SIZE], pf2[STRING_SIZE]; FILE *file; + /* calculate receiver sampling */ + dtr=(recpos[1][ntr_glob]-recpos[1][1])/(ntr_glob-1)*DH; + fprintf(FP,"\n Receiver sampling was calculatet to dtr: %f m \n ",dtr); - /* recompute the dimensions to become values of 2^m */ - nx_pot = (int)(pow(2,(ceil(log(ntr_glob)/log(2))))); - ny_pot = (int)(pow(2,(ceil(log(ns)/log(2))))); - nxy = max(nx_pot,ny_pot); - nx_pot *=2; /* double columns to avaoid filtering artefacts*/ + /* double the dimensions to avoid artefacts and recompute the dimensions to become values of 2^m */ + nx_pot = (int)(pow(2,(ceil(log(ntr_glob*2)/log(2))))); + ny_pot = (int)(pow(2,(ceil(log(ns*2)/log(2))))); - /* test: write filter to file*/ + /* test: write p-data to file*/ sprintf(pf1,"%s_data_p_.shot%i.it%i.bin",SEIS_FILE,ishot, iter); file = fopen(pf1,"wb"); for (i=1;i<=ntr_glob;i++) for (j=1;j<=ns;j++) fwrite(&data_p[i][j],sizeof(float),1,file); fclose(file); - /* test: write filter to file*/ + /* test: write vz-data to file*/ sprintf(pf1,"%s_data_vz_.shot%i.it%i.bin",SEIS_FILE,ishot, iter); file = fopen(pf1,"wb"); for (i=1;i<=ntr_glob;i++) @@ -77,25 +78,19 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos, kx=vector(1,nx_pot/2); /* calculate variable for wave field separation */ - nyq_k = 1/(2*DH); /* Nyquist of data in first dimension */ + nyq_k = 1/(2*dtr); /* Nyquist of data in first dimension */ nyq_f = 1/(2*DT); /* Nyquist of data in time dimension */ - dk = 1/(nx_pot*DH); /* Wavenumber increment */ - df = 1/(ny_pot*DT); /* Frequency increment */ + dk = 1/((nx_pot/2)*dtr); /* Wavenumber increment */ + df = 1/((ny_pot/2)*DT); /* Frequency increment */ - for (j=1;j<=ny_pot/2;j++) w[j]=df*(j-1); /* vector of angular frequencies*/ - for (j=1;j<=nx_pot/2;j++) kx[j]=dk*(j-1); /*vector of angular wavenumbers*/ + for (j=1;j<=ny_pot/2;j++) w[j]=df*(j-1)*2*PI; /* vector of angular frequencies*/ + for (j=1;j<=nx_pot/2;j++) kx[j]=dk*(j-1)*2*PI; /* vector of angular wavenumbers*/ angi = ANGTAPI * PI/180.0; /* angle at which taper begin */ ango = ANGTAPO * PI/180.0; /* angle at which taper ends */ - -// fprintf(FP,"\n Testposition 2: angi: %f; ango: %f.\n ",angi, ango); -// fprintf(FP,"\n Testposition 2: ns: %i; DT: %f.\n ",ns, DT); -// fprintf(FP,"\n Testposition 2: ntr_glob: %i; DH: %f.\n ",ntr_glob, DH); -// fprintf(FP,"\n Testposition 2: DENS: %f; VEL: %f.\n ",DENS, VEL); -// fprintf(FP,"\n Testposition 3: angi: %f; ango: %f.\n ",angi, ango); - fprintf(FP,"\n Calculate Filter for wavefield separation...\n ",angi, ango); - /* calculate Filter vor wavefield separation (after Kluever, 2008) */ + /* calculate Filter vor wavefield separation for first quadrant (after Kluever, 2008) */ + fprintf(FP,"Calculate Filter for wavefield separation... \n "); for (k=1;k<=ny_pot/2;k++) { ka = fabs(w[k])/VEL; for (i=1;i<=nx_pot/2;i++) { @@ -111,8 +106,7 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos, } } - - /* expand filter */ + /* expand filter to all 4 quadrants */ for (k=1;k<=ny_pot;k++) { for (i=1;i<=nx_pot;i++) { if(k<=ny_pot/2 && i<=nx_pot/2) { @@ -128,7 +122,7 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos, } } - /* test: write filter to file*/ + /* write filter to file*/ if (iter==1 && ishot==1) { sprintf(pf1,"%s_spectrum.filter.bin",SEIS_FILE); file = fopen(pf1,"wb"); @@ -145,7 +139,7 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos, } } - fprintf(FP,"2D FFT of p-data and vz-data... \n"); + fprintf(FP,"2D FFT of p-data and vz-data... \n"); /* transformation into fk domain */ fft2(data_pt, data_pim, ny_pot, nx_pot,1); @@ -184,12 +178,12 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos, fclose(file); - fprintf(FP,"\n Reverse 2D FFT...\n "); + fprintf(FP,"Reverse 2D FFT...\n "); /* reverse 2d fft */ fft2(data_pup, data_pupim, ny_pot, nx_pot,-1); - fprintf(FP,"\n Write PUP data to file...\n "); + fprintf(FP,"Write PUP data to file...\n "); /* write spectrum to file */ sprintf(pf2,"%s_data_pup.shot%i.it%i.bin",SEIS_FILE,ishot, iter); @@ -214,7 +208,7 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos, } outseis_glob(fp,fopen(pf,"w"), 0, data_p,recpos,recpos_loc,ntr,srcpos,1,ns,SEIS_FORMAT,ishot,1); - fprintf(FP,"\n Write PUP-data into p-data array...\n "); + fprintf(FP,"Write PUP-data into p-data array...\n "); free_matrix(filtwk, 1,ny_pot/2,1,nx_pot/2); diff --git a/src/read_par_json.c b/src/read_par_json.c index c9d3e9d..7f47dd0 100644 --- a/src/read_par_json.c +++ b/src/read_par_json.c @@ -404,19 +404,19 @@ void read_par_json(FILE *fp, char *fileinp){ fprintf(fp,"Variable WAVESEP is set to default value %i.\n",WAVESEP);} else { - if (WAVESEP=1) { + if (WAVESEP==1) { if (get_float_from_objectlist("VEL",number_readobjects,&VEL,varname_list, value_list)){ VEL=1500; - fprintf(fp,"Variable VEL is set to default value %d.\n",VEL);} + fprintf(fp,"Variable VEL is set to default value %f.\n",VEL);} if (get_float_from_objectlist("DENS",number_readobjects,&DENS,varname_list, value_list)){ DENS=1000; - fprintf(fp,"Variable DENS is set to default value %d.\n",DENS);} + fprintf(fp,"Variable DENS is set to default value %f.\n",DENS);} if (get_float_from_objectlist("ANGTAPI",number_readobjects,&ANGTAPI,varname_list, value_list)){ ANGTAPI=0; - fprintf(fp,"Variable ANGTAPI is set to default value %d.\n",ANGTAPI);} + fprintf(fp,"Variable ANGTAPI is set to default value %f.\n",ANGTAPI);} if (get_float_from_objectlist("ANGTAPO",number_readobjects,&ANGTAPO,varname_list, value_list)){ ANGTAPO=70; - fprintf(fp,"Variable ANGTAPO is set to default value %d.\n",ANGTAPO);} + fprintf(fp,"Variable ANGTAPO is set to default value %f.\n",ANGTAPO);} } } -- GitLab