Skip to content

GitLab

Commit 0ca42abe authored by niklas.thiel's avatar niklas.thiel
Browse files
Options

fixed wrong calculation of filter + minor changes

parent f1dc19ef
Hide whitespace changes
Inline Side-by-side
Showing with 67 additions and 47 deletions
src/IFOS2D.c
View file @ 0ca42abe
......@@ -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++;
}
}
/*-------------------------------------------------------------------------------*/
......
src/inseis.c
View file @ 0ca42abe
......@@ -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 */
......
src/pup.c
View file @ 0ca42abe
......@@ -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);
......
src/read_par_json.c
View file @ 0ca42abe
......@@ -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);}
}
}
......
Markdown is supported
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