Commit 44acf6f6 authored by niklas.thiel's avatar niklas.thiel

BUGFIX with using NDT

+ test on input parameter
parent 4399e0af
......@@ -156,6 +156,8 @@ int main(int argc, char **argv){
float *F_LOW_PASS_EXT=NULL;
int nfrq=0;
int FREQ_NR=1;
char pf[STRING_SIZE];
FILE *fprec, *FPL2;
......@@ -764,6 +766,7 @@ int main(int argc, char **argv){
break;
}
sectionp=matrix(1,ntr,1,ns);
break;
}
}
......@@ -1427,14 +1430,14 @@ int main(int argc, char **argv){
/* wavefield separation*/
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 ((VERBOSE)&&(MYID == 0)) fprintf(FP,"\n wavefield separation finished");
if (MYID == 0) pup(fulldata_p, fulldata_vy, FP, ntr_glob,recpos,recpos_loc,srcpos,ishot,ns,iter,1);
if ((VERBOSE==1)&&(MYID == 0)) fprintf(FP,"\n wavefield separation finished");
MPI_Barrier(MPI_COMM_WORLD);
inseis(fprec,ishot,sectionread,ntr_glob,ns,3,iter);
MPI_Bcast( &fulldata_p[1][1], ntr_glob*ns, MPI_FLOAT, 0, MPI_COMM_WORLD);
h=1;
for(i=1;i<=ntr;i++){
for(j=1;j<=ns;j++){
sectionp[h][j]=sectionread[recpos_loc[3][i]][j];
sectionp[h][j]=fulldata_p[recpos_loc[3][i]][j];
}
h++;
}
......@@ -1973,24 +1976,27 @@ int main(int argc, char **argv){
/* wavefield separation */
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 ((VERBOSE) && (MYID == 0)) fprintf(FP,"wavefield separation finished.");
if (MYID == 0) pup(fulldata_p, fulldata_vy, FP, ntr_glob,recpos,recpos_loc,srcpos,ishot,ns,iter,1);
if ((VERBOSE==1) && (MYID == 0)) fprintf(FP,"wavefield separation finished.");
MPI_Barrier(MPI_COMM_WORLD);
MPI_Bcast( &fulldata_p[1][1], ntr_glob*ns, MPI_FLOAT, 0, MPI_COMM_WORLD);
if (FORWARD_ONLY==0) {
MPI_Barrier(MPI_COMM_WORLD);
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++;
}
}
if(ntr>0) {
h=1;
for(i=1;i<=ntr;i++){
for(j=1;j<=ns;j++){
sectionp[h][j]=fulldata_p[recpos_loc[3][i]][j];
}
h++;
}
}
}
}
break;
} /* end of switch (SEISMO) */
/*---------------------------------------------------------------*/
/*---------------------- start of inversion process ------------*/
/*---------------------------------------------------------------*/
......@@ -2067,6 +2073,8 @@ int main(int argc, char **argv){
if ((TIME_FILT==1 )|| (TIME_FILT==2)){
timedomain_filt(sectionread,F_LOW_PASS,ORDER,ntr_glob,ns,1);
}
sprintf(pf,"%s_pup_data.su.shot%i",SEIS_FILE,ishot);
outseis_glob(FP,fopen(pf,"w"), 0, sectionread,recpos,recpos_loc,ntr_glob,srcpos,1,ns,SEIS_FORMAT,ishot,1);
h=1;
for(i=1;i<=ntr;i++){
for(j=1;j<=ns;j++){
......@@ -3514,14 +3522,14 @@ int main(int argc, char **argv){
catseis(sectionvy, fulldata_vy, recswitch, ntr_glob, MPI_COMM_WORLD);
catseis(sectionp, fulldata_p, recswitch, ntr_glob, MPI_COMM_WORLD);
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 ((VERBOSE) && (MYID == 0)) fprintf(FP,"\n wavefield separation finished");
MPI_Barrier(MPI_COMM_WORLD);
inseis(fprec,ishot,sectionread,ntr_glob,ns,4,iter);
if (MYID == 0) pup(fulldata_p, fulldata_vy, FP, ntr_glob,recpos,recpos_loc,srcpos,ishot,ns,iter,2);
if ((VERBOSE==1) && (MYID == 0)) fprintf(FP,"\n wavefield separation finished");
MPI_Barrier(MPI_COMM_WORLD);
MPI_Bcast( &fulldata_p[1][1], ntr_glob*ns, MPI_FLOAT, 0, MPI_COMM_WORLD);
h=1;
for(i=1;i<=ntr;i++){
for(j=1;j<=ns;j++){
sectionp[h][j]=sectionread[recpos_loc[3][i]][j];
sectionp[h][j]=fulldata_p[recpos_loc[3][i]][j];
}
h++;
}
......@@ -4359,8 +4367,7 @@ int main(int argc, char **argv){
free_matrix(sectionvzdiff,1,ntr,1,ns);
free_matrix(sectionvzdiffold,1,ntr,1,ns);
}
free_matrix(sectionread,1,ntr_glob,1,ns);
free_matrix(sectionread,1,ntr_glob,1,ns);
if((INV_STF==1)||(TIME_FILT==1) || (TIME_FILT==2)) {
/* free memory for inversion of source time function */
......
......@@ -228,7 +228,7 @@ void psource(int nt, float ** sxx, float ** syy, float ** sp,
void psource_rsg(int nt, float ** sxx, float ** syy,
float ** srcpos_loc, float ** signals, int nsrc);
void pup(float **data_p, float **data_vz, FILE *fp, int ntr_glob, int **recpos, int **recpos_loc,int ntr, float ** srcpos, int ishot, int ns, int iter, int sw );
void pup(float **data_p, float **data_vz, FILE *fp, int ntr_glob, int **recpos, int **recpos_loc, float ** srcpos, int ishot, int ns, int iter, int sw );
float *rd_sour(int *nts,FILE* fp_source);
......
......@@ -24,37 +24,39 @@
void fft2(float **array, float **arrayim, int NYG, int NXG, int dir) {
extern FILE *FP;
COMPLEX *C;
int j, i, k, nx, ny, nxy;
float **RE, **IM;
char pf1[STRING_SIZE];
FILE *file;
extern char SEIS_FILE[STRING_SIZE];
nx = NXG;
ny = NYG;
C = cplxvector(1, ny*nx);
if (dir==1) {
extern FILE *FP;
COMPLEX *C;
int j, i, k, nx, ny, nxy;
float **RE, **IM;
char pf1[STRING_SIZE];
FILE *file;
extern char SEIS_FILE[STRING_SIZE];
nx = NXG;
ny = NYG;
C = cplxvector(1, ny*nx);
if (dir==1) {
k = 1;
for (i=1;i<=nx;i++)
for (j=1;j<=ny;j++) {
for (i=1; i<=nx; i++)
for (j=1; j<=ny; j++) {
C[k].re = array[j][i];
C[k].im = 0.0;
k++;
}
forward_fft2f(C, ny, nx);
k = 1;
for (i=1;i<=nx;i++)
for (j=1;j<=ny;j++) {
for (i=1; i<=nx; i++)
for (j=1; j<=ny; j++) {
array[j][i] = C[k].re;
arrayim[j][i] = C[k].im;
k++;
......@@ -63,22 +65,25 @@ void fft2(float **array, float **arrayim, int NYG, int NXG, int dir) {
} else {
k = 1;
for (i=1;i<=nx;i++)
for (j=1;j<=ny;j++) {
for (i=1; i<=nx; i++)
for (j=1; j<=ny; j++) {
C[k].re = array[j][i];
C[k].im = arrayim[j][i];
k++;
}
inverse_fft2f(C, ny, nx);
k = 1;
for (i=1;i<=nx;i++)
for (j=1;j<=ny;j++) {
for (i=1; i<=nx; i++)
for (j=1; j<=ny; j++) {
array[j][i] = C[k].re;
arrayim[j][i] = 0.0;
k++;
}
}
free_cplxvector(C, 1, ny*nx);
}
\ No newline at end of file
}
This diff is collapsed.
......@@ -41,15 +41,15 @@ void inseis(FILE *fp, int comp, float **section, int ntr, int ns, int sws, int
}
if(sws==3){ /* open forward modeled seismic data pup*/
sprintf(data,"%s_pup.su.shot%i",SEIS_FILE,comp);
sprintf(data,"%s_pup.su.shot%d",SEIS_FILE,comp);
}
if(sws==4){ /* open forward modeled seismic data pup used for step length calculation */
sprintf(data,"%s_pup.su.step.shot%i",SEIS_FILE,comp);
sprintf(data,"%s_pup.su.step.shot%d",SEIS_FILE,comp);
}
if(sws==5){ /* open seismic data pup */
sprintf(data,"%s_pup.su.shot%i",DATA_DIR,comp);
sprintf(data,"%s_pup.su.shot%d",DATA_DIR,comp);
}
/* SWS==6 not used */
......
......@@ -2,16 +2,16 @@
* 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 <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
......@@ -20,188 +20,205 @@
*-------------------------------------------------------------
*
* 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, VERBOSE;
extern FILE *FP;
void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos, int **recpos_loc, 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, VERBOSE;
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;
if (VERBOSE) 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)))));
}
/* 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) */
if (VERBOSE) 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;
}
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;
if (VERBOSE==1) { 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)))));
/* 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) */
if (VERBOSE==1) { 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)];
}
}
/* 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) {
}
/* 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];
}
}
if (VERBOSE) 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);
}
if (VERBOSE) 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);
}
if (VERBOSE) fprintf(FP,"Reverse 2D FFT...\n ");
/* reverse 2d fft */
fft2(data_pup, data_pupim, ny_pot, nx_pot,-1);
if (VERBOSE) 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];
}
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];
}
/* save pup*/
if (VERBOSE) fprintf(FP,"Write PUP data to file...\n ");
if (sw == 1){
sprintf(pf,"%s_pup.su.shot%i",SEIS_FILE,ishot);
}
if (VERBOSE==1) { 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);
}
if (VERBOSE==1) { 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]));
}
if (sw == 2){ /* for step length calculation*/
sprintf(pf,"%s_pup.su.step.shot%i",SEIS_FILE,ishot);
}
/* 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);
}
if (VERBOSE==1) { fprintf(FP,"Reverse 2D FFT...\n "); }
/* reverse 2d fft */
fft2(data_pup, data_pupim, ny_pot, nx_pot,-1);
if (VERBOSE==1) { 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];
}
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);
}
\ No newline at end of file
}
// /* save pup*/
// if (VERBOSE==1) { fprintf(FP,"Write PUP data to file...\n "); }
//
// if (sw == 1) {
// sprintf(pf,"%s_pup.su.shot%i",SEIS_FILE,ishot);