Skip to content

GitLab

Commit 1e085c34 authored by niklas.thiel's avatar niklas.thiel
Browse files

BUGFIX 

make it possible to use number of rows unequal number of columns
parent 74505328
Hide whitespace changes
Inline Side-by-side
src/IFOS2D.c
View file @ 1e085c34
......@@ -1430,10 +1430,7 @@ int main(int argc, char **argv){
/* wavefield separation*/
if (WAVESEP==1) {
if (MYID == 0) {
fprintf(FP,"\n Testposition 0 \n");
pup(fulldata_p, fulldata_vz, FP, ntr_glob,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);
}
if (MYID == 0) pup(fulldata_p, fulldata_vz, FP, ntr_glob,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);
/* overwrite sectionp with pup data */
inseis(fprec,ishot,sectionp,ntr_glob,ns,3,iter);
}
......@@ -1960,10 +1957,10 @@ int main(int argc, char **argv){
/* wavefield separation */
if (WAVESEP==1) {
if (MYID == 0) fprintf(FP,"\n Testposition out pup");
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);
/* overwrite sectionp with pup data */
inseis(fprec,ishot,sectionp,ntr_glob,ns,3,iter);
if (MYID == 0) fprintf(FP,"\n wavefield separation finished");
// inseis(fprec,ishot,sectionp,ntr_glob,ns,3,iter);
}
break;
......
src/fft2_filt.c
View file @ 1e085c34
......@@ -29,115 +29,56 @@ void fft2(float **array, float **arrayim, int NYG, int NXG, int dir) {
COMPLEX *C;
int j, i, k, nx, ny, nxy;
float **RE, **IM;
fprintf(FP,"\n Testposition fft2_filt 1\n ");
usleep(900000);
/****************************************************/
/* recompute the dimensions to become values of 2^m */
nx = (int)(pow(2,(ceil(log(NXG)/log(2)))));
ny = (int)(pow(2,(ceil(log(NYG)/log(2)))));
nxy = max(nx,ny);
nxy *= 2;
fprintf(FP,"\n Testposition fft2_filt 1, nxy: %i\n ",nxy);
usleep(900000);
RE = matrix(1, nxy, 1, nxy);
IM = matrix(1, nxy, 1, nxy);
C = cplxvector(1, nxy*nxy);
char pf1[STRING_SIZE];
FILE *file;
extern char SEIS_FILE[STRING_SIZE];
nx = NXG;
ny = NYG;
C = cplxvector(1, ny*nx);
if (dir==1) {
/* copy array data to the temp-array */
for (i=1;i<=NXG;i++)
for (j=1;j<=NYG;j++){
RE[j][i] = array[j][i];
IM[j][i] = arrayim[j][i];
}
fprintf(FP,"\n Testposition fft2_filt 3\n ");
usleep(900000);
if (dir==1) {
k = 1;
for (i=1;i<=nxy;i++)
for (j=1;j<=nxy;j++) {
C[k].re = RE[j][i];
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, nxy, nxy);
forward_fft2f(C, ny, nx);
k = 1;
for (i=1;i<=nxy;i++)
for (j=1;j<=nxy;j++) {
RE[j][i] = C[k].re;
IM[j][i] = C[k].im;
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++;
}
} else {
k = 1;
for (i=1;i<=nxy;i++)
for (j=1;j<=nxy;j++) {
C[k].re = RE[j][i];
C[k].im = IM[j][i];
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, nxy, nxy);
inverse_fft2f(C, ny, nx);
k = 1;
for (i=1;i<=nxy;i++)
for (j=1;j<=nxy;j++) {
RE[j][i] = C[k].re;
IM[j][i] = 0.0;
for (i=1;i<=nx;i++)
for (j=1;j<=ny;j++) {
array[j][i] = C[k].re;
arrayim[j][i] = 0.0;
k++;
}
}
/* write result back into array*/
for (i=1;i<=NXG;i++)
for (j=1;j<=NYG;j++){
array[j][i]=RE[j][i];
arrayim[j][i]=IM[j][i];
}
free_cplxvector(C, 1, nxy*nxy);
free_matrix(RE, 1,nxy,1,nxy);
free_matrix(IM, 1,nxy,1,nxy);
}
/********************************************************/
/* 2D-FFTSHIFT */
/********************************************************/
// void fft2shift(float **RE, float **IM, int ny, int nx) {
//
// int j, i;
// float **A;
//
//
// A = matrix(1, ny, 1, nx);
//
// /* shift the real part */
// for (i=1;i<=nx;i++)
// for (j=1;j<=ny;j++) A[j][i] = RE[j][i];
// for (i=1;i<=nx/2;i++)
// for (j=1;j<=ny/2;j++) {
// RE[j][i] = A[j+ny/2][i+nx/2];
// RE[j+ny/2][i+nx/2] = A[j][i];
// RE[j+ny/2][i] = A[j][i+nx/2];
// RE[j][i+nx/2] = A[j+ny/2][i];
// }
//
// /* shift the imaginary part */
// for (i=1;i<=nx;i++)
// for (j=1;j<=ny;j++) A[j][i] = IM[j][i];
// for (i=1;i<=nx/2;i++)
// for (j=1;j<=ny/2;j++) {
// IM[j][i] = A[j+ny/2][i+nx/2];
// IM[j+ny/2][i+nx/2] = A[j][i];
// IM[j+ny/2][i] = A[j][i+nx/2];
// IM[j][i+nx/2] = A[j+ny/2][i];
// }
//
// free_matrix(A, 1, ny, 1, nx);
// }
\ No newline at end of file
free_cplxvector(C, 1, ny*nx);
}
\ No newline at end of file
src/fft2d.c
View file @ 1e085c34
......@@ -93,17 +93,17 @@
* and stageBuff, while also performing any conversions float<->double.
*/
#define LoadRow(buffer,row,cols) if (1){\
#define LoadRow(buffer,row,rows,cols) if (1){\
register int j,k;\
for (j = row*cols, k = 0 ; k < cols ; j++, k++){\
for (j = row, k = 0 ; k < cols ; j+=rows, k++){\
stageBuff[k].re = buffer[j].re;\
stageBuff[k].im = buffer[j].im;\
}\
} else {}
#define StoreRow(buffer,row,cols) if (1){\
#define StoreRow(buffer,row,rows,cols) if (1){\
register int j,k;\
for (j = row*cols, k = 0 ; k < cols ; j++, k++){\
for (j = row, k = 0 ; k < cols ; j+=rows, k++){\
buffer[j].re = stageBuff[k].re;\
buffer[j].im = stageBuff[k].im;\
}\
......@@ -111,27 +111,29 @@
#define LoadCol(buffer,col,rows,cols) if (1){\
register int j,k;\
for (j = i,k = 0 ; k < rows ; j+=cols, k++){\
stageBuff[k].re = buffer[j].re;\
stageBuff[k].im = buffer[j].im;\
for (j = col*rows,k = 0 ; k < rows ; j++, k++){\
stageBuff1[k].re = buffer[j].re;\
stageBuff1[k].im = buffer[j].im;\
}\
} else {}
#define StoreCol(buffer,col,rows,cols) if (1){\
register int j,k;\
for (j = i,k = 0 ; k < rows ; j+=cols, k++){\
buffer[j].re = stageBuff[k].re;\
buffer[j].im = stageBuff[k].im;\
for (j = col*rows,k = 0 ; k < rows ; j++, k++){\
buffer[j].re = stageBuff1[k].re;\
buffer[j].im = stageBuff1[k].im;\
}\
} else {}
/* Do something useful with an error message */
#define handle_error(msg) fprintf(stderr,msg)
#define handle_error(msg) fprintf(FP,msg)
DCOMPLEX *stageBuff; /* buffer to hold a row or column at a time */
COMPLEX *bigBuff; /* a pointer to a float input array */
DCOMPLEX *stageBuff1; /* buffer to hold a row or column at a time */
COMPLEX *bigBuff; /* a pointer to a float input array */
DCOMPLEX *bigBuffd; /* a pointer to a double input array */
extern FILE *FP;
/* Allocate space for stageBuff */
......@@ -146,11 +148,25 @@ int allocateBuffer(int size)
else return(NO_ERROR);
}
/* Allocate space for stageBuff */
int allocateBuffer1(int size)
{
stageBuff1 = (DCOMPLEX*) malloc(size*sizeof(DCOMPLEX));
if(stageBuff1==(DCOMPLEX*)NULL)
{
handle_error("Insufficient storage for fft buffers");
return(ERROR);
}
else return(NO_ERROR);
}
/* Free space for stageBuff */
void freeBuffer(void)
{
free((char*)stageBuff);
free((char*)stageBuff1);
}
......@@ -520,26 +536,29 @@ void FFT842(int direction, int n, DCOMPLEX *b)
*/
int fft2f(COMPLEX *array, int rows, int cols, int direction)
{
int i, maxsize, errflag;
int i;
if(!power_of_2(rows) || !power_of_2(cols)) {
handle_error("fft: input array must have dimensions a power of 2");
handle_error("\n ****************************************************");
handle_error("\n fft: input array must have dimensions a power of 2! \n");
handle_error(" ****************************************************\n");
return(ERROR);
}
/* Allocate 1D buffer */
bigBuff = array;
maxsize = rows>cols ? rows : cols;
errflag = allocateBuffer(maxsize);
if(errflag != NO_ERROR) return(errflag);
// maxsize = rows>cols ? rows : cols;
allocateBuffer(cols);
allocateBuffer1(rows);
// if(errflag != NO_ERROR) return(errflag);
/* Compute transform row by row */
if(cols>1)
for(i=0;i<rows;i++)
{
LoadRow(bigBuff,i,cols);
LoadRow(bigBuff,i,rows,cols);
FFT842(direction,cols,stageBuff);
StoreRow(bigBuff,i,cols);
StoreRow(bigBuff,i,rows,cols);
}
/* Compute transform column by column */
......@@ -547,7 +566,7 @@ int fft2f(COMPLEX *array, int rows, int cols, int direction)
for(i=0;i<cols;i++)
{
LoadCol(bigBuff,i,rows,cols);
FFT842(direction,rows,stageBuff);
FFT842(direction,rows,stageBuff1);
StoreCol(bigBuff,i,rows,cols);
}
......@@ -555,48 +574,7 @@ int fft2f(COMPLEX *array, int rows, int cols, int direction)
return(NO_ERROR);
}
/*
* Compute 2D DFT on double data:
* forward if direction==FFT_FORWARD,
* inverse if direction==FFT_INVERSE.
*/
// int fft2d(DCOMPLEX *array, int rows, int cols, int direction)
// {
// int i, maxsize, errflag;
//
// if(!power_of_2(rows) || !power_of_2(cols)) {
// handle_error("fft: input array must have dimensions a power of 2");
// return(ERROR);
// }
//
// /* Allocate 1D buffer */
// bigBuffd = array;
// maxsize = rows>cols ? rows : cols;
// errflag = allocateBuffer(maxsize);
// if(errflag != NO_ERROR) return(errflag);
//
// /* Compute transform row by row */
// if(cols>1)
// for(i=0;i<rows;i++)
// {
// LoadRow(bigBuffd,i,cols);
// FFT842(direction,cols,stageBuff);
// StoreRow(bigBuffd,i,cols);
// }
//
// /* Compute transform column by column */
// if(rows>1)
// for(i=0;i<cols;i++)
// {
// LoadCol(bigBuffd,i,rows,cols);
// FFT842(direction,rows,stageBuff);
// StoreCol(bigBuffd,i,rows,cols);
// }
//
// freeBuffer();
// return(NO_ERROR);
// }
//
// /* Finally, the entry points we announce in kube_fft.h */
/* Perform forward 2D transform on a COMPLEX array. */
......@@ -610,15 +588,3 @@ int inverse_fft2f(COMPLEX *array, int rows, int cols)
{
return(fft2f(array, rows, cols, FFT_INVERSE));
}
// /* Perform forward 2D transform on a DCOMPLEX array. */
// int forward_fft2d(DCOMPLEX *array, int rows, int cols)
// {
// return(fft2d(array, rows, cols, FFT_FORWARD));
// }
//
// /* Perform inverse 2D transform on a DCOMPLEX array. */
// int inverse_fft2d(DCOMPLEX *array, int rows, int cols)
// {
// return(fft2d(array, rows, cols, FFT_INVERSE));
//}
src/inseis.c
View file @ 1e085c34
......@@ -40,14 +40,14 @@ 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%d",DATA_DIR,comp);
sprintf(data,"%s_pup.su.shot%i.it%i",DATA_DIR,comp,iter);
}
if(sws==4){ /* open seismic data pup used for step length calculation */
sprintf(data,"%s_pup.su.step.shot%d",DATA_DIR,comp);
sprintf(data,"%s_pup.su.step.shot%i.it%i",DATA_DIR,comp,iter);
}
/* sws 3 -- 6 not used */
/* sws 5 -- 6 not used */
if(sws==7){ /* open convolved seismic data vx */
sprintf(data,"%s_vx.su.conv.shot%d",DATA_DIR,comp);
......
src/pup.c
View file @ 1e085c34
......@@ -34,48 +34,71 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
extern int SEIS_FORMAT;
extern FILE *FP;
int i, j, k;
int i, j, k, nx_pot, ny_pot,nxy;
float nyq_k, nyq_f, dk, df, angi, ango, ka;
float *w, *kx;
float **filtwk, **filtwk_full;
float **data_pim,**data_vyim,**data_pup,**data_pupim;
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;
/* temporary array */
filtwk=matrix(1,ns/2,1,ntr_glob/2);
filtwk_full=matrix(1,ns,1,ntr_glob);
data_pup=matrix(1,ns,1,ntr_glob);
data_pim=matrix(1,ns,1,ntr_glob);
data_vyim=matrix(1,ns,1,ntr_glob);
data_pupim=matrix(1,ns,1,ntr_glob);
/* 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*/
w=vector(1,ns/2);
kx=vector(1,ntr_glob/2);
/* test: write filter to file*/
sprintf(pf1,"%s_data_p_output.bin",SEIS_FILE);
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*/
sprintf(pf1,"%s_data_vz_output.bin",SEIS_FILE);
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*DH); /* Nyquist of data in first dimension */
nyq_f = 1/(2*DT); /* Nyquist of data in time dimension */
dk = 1/(ntr_glob*DH); /* Wavenumber increment */
df = 1/(ns*DT); /* Frequency increment */
for (j=1;j<=ns/2;j++) w[j]=df*(j-1); /* vector of angular frequencies*/
for (j=1;j<=ntr_glob/2;j++) kx[j]=dk*(j-1); /*vector of angular wavenumbers*/
fprintf(FP,"\n kx(8): %7.4f\n",kx[8]);
dk = 1/(nx_pot*DH); /* Wavenumber increment */
df = 1/(ny_pot*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*/
angi = ANGTAPI * PI/180.0;
ango = ANGTAPO * PI/180.0;
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: 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) */
for (k=1;k<=ns/2;k++) {
if (k==1) fprintf(FP,"\n Testposition 2.5: w(%i): %f \n",k,w[k]);
for (k=1;k<=ny_pot/2;k++) {
ka = fabs(w[k])/VEL;
if (k==3000) fprintf(FP,"\n Testposition 2.5: w(%i): %f \n",k,w[k]);
for (i=1;i<=ntr_glob/2;i++) {
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{
......@@ -87,59 +110,98 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
}
}
}
fprintf(FP,"\n Testposition 3: angi: %f; ango: %f.\n ",angi, ango);
/* expand filter */
for (k=1;k<=ns;k++) {
fprintf(FP,"k: %i \t",k);
for (i=1;i<=ntr_glob;i++) {
if (k==3000) fprintf(FP,"i: %i \n",i);
if(k<=ns/2 && i<=ntr_glob/2) {
filtwk_full[k][i]=filtwk[(ns/2)-(k-1)][(ntr_glob/2)-(i-1)];
}else if(k<=ns/2 && i>ntr_glob/2) {
filtwk_full[k][i]=filtwk[(ns/2)-(k-1)][i-(ntr_glob/2)];
}else if(k>ns/2 && i<=ntr_glob/2) {
filtwk_full[k][i]=filtwk[k-(ns/2)][(ntr_glob/2)-(i-1)];
}else if(k>ns/2 && i>ntr_glob/2) {
filtwk_full[k][i]=filtwk[k-(ns/2)][i-(ntr_glob/2)];
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)];
}
}
}
/* transformation into fk domain */
fft2(data_p, data_pim, ns, ntr_glob,1);
fft2(data_vy, data_vyim, ns, ntr_glob,1);
/* write filter to file*/
if (iter==1){
sprintf(pf1,"$s_spectrum.filt.shot%i.it%i.bin",SEIS_FILE,ishot, iter);
/* test: write filter to file*/
sprintf(pf1,"%s_spectrum.filt.shot%i.it%i.bin",SEIS_FILE,ishot,iter);
file = fopen(pf1,"wb");
for (i=1;i<=ns;i++)
for (j=1;j<=ntr_glob;j++) fwrite(&filtwk_full[j][i],sizeof(float),1,FP);
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 */
sprintf(pf2,"$s_spectrum.p.real.shot%i.it%i.bin",SEIS_FILE,ishot, iter);
sprintf(pf2,"%s_spectrum.p.real.shot%i.it%i.bin",SEIS_FILE,ishot, iter);
file = fopen(pf2,"wb");
for (i=1;i<=ns;i++)
for (j=1;j<=ntr_glob;j++) fwrite(&data_p[j][i],sizeof(float),1,FP);
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 */
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<=ns;k++) {
for (i=1;i<=ntr_glob;i++) {
data_pup[k][i]=0.5*(data_p[k][i]-(filtwk_full[k][i]*data_vy[k][i]));
data_pupim[k][i]=0.5*(data_pim[k][i]-(filtwk_full[k][i]*data_vyim[k][i]));
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 */
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,"\n Testposition 4\n ");
fprintf(FP,"\n Reverse 2D FFT...\n ");
/* reverse 2d fft */
fft2(data_pup, data_pupim, ns, ntr_glob,-1);
fft2(data_pup, data_pupim, ny_pot, nx_pot,-1);
fprintf(FP,"\n Write PUP data to file...\n ");
/* write spectrum to file */
sprintf(pf2,"%s_pup.hand.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);