/*----------------------------------------------------------------------------------------- * Copyright (C) 2016 For the list of authors, see file AUTHORS. * * This file is part of IFOS. * * IFOS is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, version 2.0 of the License only. * * IFOS is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with IFOS. See file COPYING and/or . -----------------------------------------------------------------------------------------*/ /* *------------------------------------------------------------- * * perform wavefield separation in FK domain and save it to file * * *------------------------------------------------------------- */ #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 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; 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); /* 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 p-data to file*/ if (ishot==1 && iter==1) { 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 vz-data to file*/ if (ishot==1 && iter==1) { 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++) for (j=1;j<=ns;j++) fwrite(&data_vy[i][j],sizeof(float),1,file); fclose(file); } /* temporary array */ filtwk=matrix(1,ny_pot/2,1,nx_pot/2); filtwk_full=matrix(1,ny_pot,1,nx_pot); data_pim=matrix(1,ny_pot,1,nx_pot); data_vyim=matrix(1,ny_pot,1,nx_pot); data_pup=matrix(1,ny_pot,1,nx_pot); data_pupim=matrix(1,ny_pot,1,nx_pot); data_pt=matrix(1,ny_pot,1,nx_pot); data_vyt=matrix(1,ny_pot,1,nx_pot); w=vector(1,ny_pot/2); kx=vector(1,nx_pot/2); /* calculate variable for wave field separation */ 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/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)*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 */ /* 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++) { if(kx[i] < ka*sin(angi)) { filtwk[k][i] = (DENS*fabs(w[k]))/sqrt((ka*ka)-(kx[i]*kx[i])); }else{ if(kx[i] < ka*sin(ango) && kx[i] > ka*sin(angi)){ filtwk[k][i] = 0.5*(1.0+cos(PI/((ka*sin(ango))-(ka*sin(angi))) * (kx[i]-(ka*sin(angi))))) * (DENS*fabs(w[k]))/sqrt((ka*ka)-(kx[i]*kx[i])); }else { filtwk[k][i] = 0.0; } } } } /* 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) { filtwk_full[k][i]=filtwk[(ny_pot/2)-(k-1)][(nx_pot/2)-(i-1)]; }else if(k<=ny_pot/2 && i>nx_pot/2) { filtwk_full[k][i]=filtwk[(ny_pot/2)-(k-1)][i-(nx_pot/2)]; }else if(k>ny_pot/2 && i<=nx_pot/2) { filtwk_full[k][i]=filtwk[k-(ny_pot/2)][(nx_pot/2)-(i-1)]; }else if(k>ny_pot/2 && i>nx_pot/2) { filtwk_full[k][i]=filtwk[k-(ny_pot/2)][i-(nx_pot/2)]; } } } /* write filter to file*/ if (iter==1 && ishot==1) { sprintf(pf1,"%s_spectrum.filter.bin",SEIS_FILE); file = fopen(pf1,"wb"); for (i=1;i<=nx_pot;i++) for (j=1;j<=ny_pot;j++) fwrite(&filtwk_full[j][i],sizeof(float),1,file); fclose(file); } /* copy data to temp array */ for(k=1;k<=ntr_glob;k++) { for (i=1;i<=ns;i++) { data_pt[i][k]=data_p[k][i]; data_vyt[i][k]=data_vy[k][i]; } } 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); /* write spectrum to file */ if (ishot==1 && iter==1) { sprintf(pf2,"%s_spectrum.p.real.shot%i.it%i.bin",SEIS_FILE,ishot, iter); file = fopen(pf2,"wb"); for (i=1;i<=nx_pot;i++) for (j=1;j<=ny_pot;j++) fwrite(&data_pt[j][i],sizeof(float),1,file); fclose(file); } fft2(data_vyt, data_vyim, ny_pot, nx_pot,1); /* write spectrum to file */ if (ishot==1 && iter==1) { sprintf(pf2,"%s_spectrum.vz.real.shot%i.it%i.bin",SEIS_FILE,ishot, iter); file = fopen(pf2,"wb"); for (i=1;i<=nx_pot;i++) for (j=1;j<=ny_pot;j++) fwrite(&data_vyt[j][i],sizeof(float),1,file); fclose(file); } fprintf(FP,"Perform wavefield separation... \n"); /* perform wavefield separation */ for (k=1;k<=nx_pot;k++) { for (i=1;i<=ny_pot;i++) { data_pup[i][k]=0.5*(data_pt[i][k]-(filtwk_full[i][k]*data_vyt[i][k])); data_pupim[i][k]=0.5*(data_pim[i][k]-(filtwk_full[i][k]*data_vyim[i][k])); } } /* write spectrum to file */ if (ishot==1 && iter==1) { sprintf(pf2,"%s_spectrum.pup.real.shot%i.it%i.bin",SEIS_FILE,ishot, iter); file = fopen(pf2,"wb"); for (i=1;i<=nx_pot;i++) for (j=1;j<=ny_pot;j++) fwrite(&data_pup[j][i],sizeof(float),1,file); fclose(file); } fprintf(FP,"Reverse 2D FFT...\n "); /* reverse 2d fft */ fft2(data_pup, data_pupim, ny_pot, nx_pot,-1); fprintf(FP,"Write PUP data to file...\n "); /* write spectrum to file */ if (ishot==1 && iter==1) { sprintf(pf2,"%s_data_pup.shot%i.it%i.bin",SEIS_FILE,ishot, iter); file = fopen(pf2,"wb"); for (i=1;i<=nx_pot;i++) for (j=1;j<=ny_pot;j++) fwrite(&data_pup[j][i],sizeof(float),1,file); fclose(file); } fprintf(FP,"Write PUP-data into p-data array...\n "); /* overwrite p data */ for (k=1;k<=ntr_glob;k++) { for (i=1;i<=ns;i++) { data_p[k][i]=data_pup[i][k]; } } /* save pup*/ if (sw == 1){ sprintf(pf,"%s_pup.su.shot%i",SEIS_FILE,ishot); } if (sw == 2){ /* for step length calculation*/ sprintf(pf,"%s_pup.su.step.shot%i",SEIS_FILE,ishot); } outseis_glob(fp,fopen(pf,"w"), 0, data_p,recpos,recpos_loc,ntr,srcpos,1,ns,SEIS_FORMAT,ishot,1); free_matrix(filtwk, 1,ny_pot/2,1,nx_pot/2); free_matrix(filtwk_full, 1,ny_pot,1,nx_pot); free_matrix(data_pim, 1,ny_pot,1,nx_pot); free_matrix(data_vyim, 1,ny_pot,1,nx_pot); free_matrix(data_pup, 1,ny_pot,1,nx_pot); free_matrix(data_pupim, 1,ny_pot,1,nx_pot); free_matrix(data_pt, 1,ny_pot,1,nx_pot); free_matrix(data_vyt, 1,ny_pot,1,nx_pot); free_vector(w,1,ny_pot/2); free_vector(kx,1,ny_pot/2); }