From 2a3163806087ae6e17bd6d1d6f5d3eb7db9b4112 Mon Sep 17 00:00:00 2001 From: Niklas Thiel Date: Mon, 29 Feb 2016 11:13:39 +0100 Subject: [PATCH] wavefield seperation added --- par/in_and_out/IFOS2D_FW_all_parameters.json | 7 + par/in_and_out/IFOS2D_INV_all_parameters.json | 8 + src/IFOS2D.c | 44 +- src/Makefile | 5 +- src/exchange_par.c | 21 + src/fd.h | 32 +- src/fft2_filt.c | 143 ++++ src/fft2d.c | 624 ++++++++++++++++++ src/globvar.h | 5 +- src/inseis.c | 8 + src/pup.c | 170 +++++ src/read_par_json.c | 28 + src/seismo_ssg.c | 4 +- src/util.c | 20 + 14 files changed, 1104 insertions(+), 15 deletions(-) create mode 100644 src/fft2_filt.c create mode 100644 src/fft2d.c create mode 100644 src/pup.c diff --git a/par/in_and_out/IFOS2D_FW_all_parameters.json b/par/in_and_out/IFOS2D_FW_all_parameters.json index 065ad2b..e203ced 100644 --- a/par/in_and_out/IFOS2D_FW_all_parameters.json +++ b/par/in_and_out/IFOS2D_FW_all_parameters.json @@ -65,6 +65,13 @@ "XREC2, YREC2" : "93.0 , 0.2", "NGEOPH" : "80", +"Wavefield Separation" : "comment", + "WAVESEP" : "0", + "VEL" : "1500.0", + "DENS" : "1000.0", + "ANGTAPI" : "0.0", + "ANGTAPO" : "70.0", + "Seismograms" : "comment", "NDT" : "1", "SEIS_FORMAT" : "1", diff --git a/par/in_and_out/IFOS2D_INV_all_parameters.json b/par/in_and_out/IFOS2D_INV_all_parameters.json index bc45d83..3ac26fa 100644 --- a/par/in_and_out/IFOS2D_INV_all_parameters.json +++ b/par/in_and_out/IFOS2D_INV_all_parameters.json @@ -65,6 +65,14 @@ "XREC2, YREC2" : "93.0 , 0.2", "NGEOPH" : "80", + +"Wavefield Separation" : "comment", + "WAVESEP" : "0", + "VEL" : "1500", + "DENS" : "1000", + "ANGTAPI" : "0", + "ANGTAPO" : "70", + "Seismograms" : "comment", "NDT" : "1", "SEIS_FORMAT" : "1", diff --git a/src/IFOS2D.c b/src/IFOS2D.c index 919076b..98eb4a6 100644 --- a/src/IFOS2D.c +++ b/src/IFOS2D.c @@ -699,21 +699,17 @@ int main(int argc, char **argv){ fulldata_curl = matrix(1,ntr_glob,1,NT); break; case 5 : /* everything except curl and div*/ - switch (WAVETYPE) { - case 1: + if (WAVETYPE==1) { fulldata_vx = matrix(1,ntr_glob,1,NT); fulldata_vy = matrix(1,ntr_glob,1,NT); - break; - - case 2: + } + if (WAVETYPE==2) { fulldata_vz = matrix(1,ntr_glob,1,NT); - break; - - case 3: + } + if (WAVETYPE==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; @@ -1442,7 +1438,18 @@ int main(int argc, char **argv){ } catseis(sectionp, fulldata_p, recswitch, ntr_glob, MPI_COMM_WORLD); if (MYID==0) saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,fulldata_curl,fulldata_div,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1); + + /* 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); + } + /* overwrite sectionp with pup data */ + inseis(fprec,ishot,sectionp,ntr_glob,ns,3,iter); + } break; + } /* end of switch (SEISMO) */ @@ -1985,7 +1992,15 @@ int main(int argc, char **argv){ catseis(sectionvz, fulldata_vz, recswitch, ntr_glob, MPI_COMM_WORLD); } catseis(sectionp, fulldata_p, recswitch, ntr_glob, MPI_COMM_WORLD); - if (MYID==0) saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,fulldata_curl,fulldata_div,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1); + if (MYID==0) saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,fulldata_curl,fulldata_div,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1); + + /* wavefield separation */ + if (WAVESEP==1) { + if (MYID == 0) fprintf(FP,"\n Testposition out pup"); + 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); + } break; } /* end of switch (SEISMO) */ @@ -3485,6 +3500,15 @@ int main(int argc, char **argv){ } + /* wavefield separation */ + if (WAVESEP==1) { + catseis(sectionvz, fulldata_vz, recswitch, ntr_glob, MPI_COMM_WORLD); + catseis(sectionp, fulldata_p, recswitch, ntr_glob, MPI_COMM_WORLD); + if (MYID == 0) pup(fulldata_p, fulldata_vz, FP, ntr_glob,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,2); + /* overwrite sectionp with pup data */ + inseis(fprec,ishot,sectionp,ntr_glob,ns,3,iter); + } + /*-------------------------------------------------------------------------------*/ /*------------------ end loop over timesteps (forward model) step length -------*/ /*-------------------------------------------------------------------------------*/ diff --git a/src/Makefile b/src/Makefile index ff2e6ea..5acb442 100644 --- a/src/Makefile +++ b/src/Makefile @@ -163,7 +163,10 @@ IFOS2D= \ readmod_viscac.c \ time_window_glob.c \ create_trkill_table.c \ - filter_frequencies.c + filter_frequencies.c \ + pup.c \ + fft2_filt.c \ + fft2d.c SNAPMERGE_OBJ = $(SNAPMERGE_SCR:%.c=%.o) diff --git a/src/exchange_par.c b/src/exchange_par.c index 6e1ed48..1f94944 100644 --- a/src/exchange_par.c +++ b/src/exchange_par.c @@ -117,6 +117,9 @@ void exchange_par(void){ extern int WOLFE_TRY_OLD_STEPLENGTH; extern float WOLFE_C1_SL; extern float WOLFE_C2_SL; + + extern int WAVESEP; + extern float VEL, DENS, ANGTAPI, ANGTAPO; /* definition of local variables */ @@ -218,6 +221,12 @@ void exchange_par(void){ fdum[68]=TRKILL_OFFSET_LOWER; fdum[69]=TRKILL_OFFSET_UPPER; + fdum[70]=VEL; + fdum[71]=DENS; + fdum[72]=ANGTAPI; + fdum[73]=ANGTAPO; + + /***********/ /* Integer */ /***********/ @@ -371,6 +380,10 @@ void exchange_par(void){ idum[114]=TRKILL_OFFSET; idum[115]=TRKILL_STF_OFFSET; + + idum[116]=WAVESEP; + + } /** if (MYID == 0) **/ @@ -498,6 +511,12 @@ void exchange_par(void){ TRKILL_STF_OFFSET_UPPER=fdum[67]; TRKILL_OFFSET_LOWER=fdum[68]; TRKILL_OFFSET_UPPER=fdum[69]; + + VEL=fdum[70]; + DENS=fdum[71]; + ANGTAPI=fdum[72]; + ANGTAPO=fdum[73]; + /***********/ /* Integer */ @@ -655,6 +674,8 @@ void exchange_par(void){ TRKILL_OFFSET=idum[114]; TRKILL_STF_OFFSET=idum[115]; + + WAVESEP=idum[116]; MPI_Bcast(&FL[1],L,MPI_FLOAT,0,MPI_COMM_WORLD); diff --git a/src/fd.h b/src/fd.h index c9c2035..83efd66 100644 --- a/src/fd.h +++ b/src/fd.h @@ -10,6 +10,7 @@ #include #include #include +#include #include #define iround(x) ((int)(floor)(x+0.5)) @@ -24,8 +25,32 @@ #define REQUEST_COUNT 4 #define WORKFLOW_MAX_VAR 12 -/* declaration of functions */ +/***************/ +/* FFT2D stuff */ +/***************/ +typedef struct {float re; float im;} COMPLEX; +typedef struct {double re; double im;} DCOMPLEX; + +#ifndef ERROR +#define ERROR -1 +#define NO_ERROR 0 +#endif +#define FFT_FORWARD 0 +#define FFT_INVERSE 1 + +#define TWOPI 6.2831853071795865 /* 2.0 * PI */ +#define HALFPI 1.5707963267948966 /* PI / 2.0 */ +#define PI8 0.392699081698724 /* PI / 8.0 */ +#define RT2 1.4142135623731 /* sqrt(2.0) */ +#define IRT2 0.707106781186548 /* 1.0/sqrt(2.0) */ + +extern int forward_fft2f(COMPLEX *array, int rows, int cols); /* Perform forward 2D transform on a COMPLEX array. */ +extern int inverse_fft2f(COMPLEX *array, int rows, int cols); /* Perform inverse 2D transform on a COMPLEX array. */ + +/****************************/ +/* declaration of functions */ +/****************************/ void window_cos(float **win, int npad, int nsrc, float it1, float it2, float it3, float it4); @@ -125,6 +150,8 @@ void FFT_filt(float ** data, float freqshift, int ntr, int ns,int method); void FFT(int isign, unsigned long nlog, float *re, float *im); /* NR version*/ +void fft2(float **array, float **arrayim, int NYG, int NXG, int dir); /* 2d fft filtering*/ + float *filter_frequencies(int *nfrq); float *holbergcoeff(void); @@ -201,6 +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 ); float *rd_sour(int *nts,FILE* fp_source); @@ -464,6 +492,7 @@ void err(char err_text[]); void warning(char warn_text[]); double maximum(float **a, int nx, int ny); float *vector(int nl, int nh); +COMPLEX *cplxvector(int nstart, int nend); int *ivector(int nl, int nh); double *dvector(int nl, int nh); float **fmatrix(int nrl, int nrh, int ncl, int nch); @@ -473,6 +502,7 @@ float **matrix(int nrl, int nrh, int ncl, int nch); int **imatrix(int nrl, int nrh, int ncl, int nch); float ***f3tensor(int nrl, int nrh, int ncl, int nch,int ndl, int ndh); void free_vector(float *v, int nl, int nh); +void free_cplxvector(COMPLEX *v, int nstart, int nend); void free_dvector(double *v, int nl, int nh); void free_ivector(int *v, int nl, int nh); void free_matrix(float **m, int nrl, int nrh, int ncl, int nch); diff --git a/src/fft2_filt.c b/src/fft2_filt.c new file mode 100644 index 0000000..b88fa5c --- /dev/null +++ b/src/fft2_filt.c @@ -0,0 +1,143 @@ +/*----------------------------------------------------------------------------------------- + * 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 . + -----------------------------------------------------------------------------------------*/ + +/*------------------------------------------------------------------------ + * 2D-FFT preparation + * ----------------------------------------------------------------------*/ + +#include "fd.h" + +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; + + 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); + + + + /* 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]; + C[k].im = 0.0; + k++; + } + forward_fft2f(C, nxy, nxy); + 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; + 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]; + k++; + } + inverse_fft2f(C, nxy, nxy); + 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; + 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 diff --git a/src/fft2d.c b/src/fft2d.c new file mode 100644 index 0000000..928702c --- /dev/null +++ b/src/fft2d.c @@ -0,0 +1,624 @@ +/* + ------------------------------------------------------------------------------ + * + * Copyright (C) 2013, see file AUTHORS for list of authors/contributors + * + * This file is part of PROTEUS. + * + * PROTEUS 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. + * + * PROTEUS 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 PROTEUS. See file COPYING and/or + * . + * + ------------------------------------------------------------------------------ +*/ + + +/************************************************************** + * kube-gustavson-fft.c + * + * Forward and inverse discrete 2D Fourier transforms. + * + * Originally a FORTRAN implementation published in: + * Programming for Digital Signal Processing, IEEE Press 1979, + * Section 1, by G. D. Bergland and M. T. Dolan. + * Ported long ago from Fortran to old-fashioned Standard C by + * someone who failed to put his or her name in the source. + * Ported from Standard C to ANSI C by Stefan Gustavson + * (stegu@itn.liu.se) 2003-10-20. + * + * This is a pretty fast FFT algorithm. Not super-optimized + * lightning-fast, but very good considering its age and + * relative simplicity. There are several places in the code + * where clock cycles could be saved, for example by getting rid + * of some copying of data back and forth, but the savings + * would not be that great. If you want optimally fast FFTs, + * consider the excellent FFTW package. (http://www.fftw.org) + * + * This software is in the public domain. + * + ************************************************************** + * + * int forward_fft2f(COMPLEX *array, int rows, int cols) + * int inverse_fft2f(COMPLEX *array, int rows, int cols) + * int forward_fft2d(DCOMPLEX *array, int rows, int cols) + * int inverse_fft2d(DCOMPLEX *array, int rows, int cols) + * + * These functions compute the forward and inverse DFT's, respectively, + * of a single-precision COMPLEX or DCOMPLEX array by means of an + * FFT algorithm. + * + * The result is a COMPLEX/DCOMPLEX array of the same size, returned + * in the same space as the input array. That is, the original array + * is overwritten and destroyed. + * + * Rows and columns must each be an integral power of 2. + * + * These routines return integer value ERROR if an error was detected, + * NO_ERROR otherwise. + * + * This particular implementation of the DFT uses the transform pair + * defined as follows: + * Let there be two complex arrays each with n rows and m columns. + * Index them as + * f(x,y): 0 <= x <= m - 1, 0 <= y <= n - 1 + * F(u,v): -m/2 <= u <= m/2 - 1, -n/2 <= v <= n/2 - 1 + * + * Then the forward and inverse transforms are related as follows. + * Forward: + * F(u,v) = \sum_{x=0}^{m-1} \sum_{y=0}^{n-1} + * f(x,y) \exp{-2\pi i (ux/m + vy/n)} + * + * Inverse: + * f(x,y) = 1/(mn) \sum_{u=-m/2}^{m/2-1} \sum_{v=-n/2}^{n/2-1} + * F(u,v) \exp{2\pi i (ux/m + vy/n)} + * + * Therefore, the transforms have these properties: + * 1. \sum_x \sum_y f(x,y) = F(0,0) + * 2. m n \sum_x \sum_y |f(x,y)|^2 = \sum_u \sum_v |F(u,v)|^2 + * + */ +#include "fd.h" + +/* + * These somewhat awkward macros move data between the input/output array + * and stageBuff, while also performing any conversions float<->double. + */ + +#define LoadRow(buffer,row,cols) if (1){\ + register int j,k;\ + for (j = row*cols, k = 0 ; k < cols ; j++, k++){\ + stageBuff[k].re = buffer[j].re;\ + stageBuff[k].im = buffer[j].im;\ + }\ + } else {} + +#define StoreRow(buffer,row,cols) if (1){\ + register int j,k;\ + for (j = row*cols, k = 0 ; k < cols ; j++, k++){\ + buffer[j].re = stageBuff[k].re;\ + buffer[j].im = stageBuff[k].im;\ + }\ + } else {} + +#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;\ + }\ + } 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;\ + }\ + } else {} + + +/* Do something useful with an error message */ +#define handle_error(msg) fprintf(stderr,msg) + +DCOMPLEX *stageBuff; /* 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 */ + + +/* Allocate space for stageBuff */ +int allocateBuffer(int size) +{ + stageBuff = (DCOMPLEX*) malloc(size*sizeof(DCOMPLEX)); + if(stageBuff==(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); +} + + +/* See if exactly one bit is set in integer argument */ +int power_of_2(int n) +{ + int bits=0; + while(n) { + bits += n & 1; + n >>= 1; + } + return(bits==1); +} + + +/* Get binary log of integer argument - exact if n is a power of 2 */ +int fastlog2(int n) +{ + int log = -1; + while(n) { + log++; + n >>= 1; + } + return(log); +} + + +/* Radix-2 iteration subroutine */ +void R2TX(int nthpo, DCOMPLEX *c0, DCOMPLEX *c1) +{ + int k,kk; + double *cr0, *ci0, *cr1, *ci1, r1, fi1; + + cr0 = &(c0[0].re); + ci0 = &(c0[0].im); + cr1 = &(c1[0].re); + ci1 = &(c1[0].im); + + for(k=0; k1) + { + cr1[kk] = c4*(br0-br1) - s4*(bi0-bi1); + cr2[kk] = c2*(br2-bi3) - s2*(bi2+br3); + cr3[kk] = c6*(br2+bi3) - s6*(bi2-br3); + ci1[kk] = c4*(bi0-bi1) + s4*(br0-br1); + ci2[kk] = c2*(bi2+br3) + s2*(br2-bi3); + ci3[kk] = c6*(bi2-br3) + s6*(br2+bi3); + tr = IRT2*(br5-bi5); + ti = IRT2*(br5+bi5); + cr4[kk] = c1*(br4+tr) - s1*(bi4+ti); + ci4[kk] = c1*(bi4+ti) + s1*(br4+tr); + cr5[kk] = c5*(br4-tr) - s5*(bi4-ti); + ci5[kk] = c5*(bi4-ti) + s5*(br4-tr); + tr = -IRT2*(br7+bi7); + ti = IRT2*(br7-bi7); + cr6[kk] = c3*(br6+tr) - s3*(bi6+ti); + ci6[kk] = c3*(bi6+ti) + s3*(br6+tr); + cr7[kk] = c7*(br6-tr) - s7*(bi6-ti); + ci7[kk] = c7*(bi6-ti) + s7*(br6-tr); + } + else + { + cr1[kk] = br0 - br1; + cr2[kk] = br2 - bi3; + cr3[kk] = br2 + bi3; + ci1[kk] = bi0 - bi1; + ci2[kk] = bi2 + br3; + ci3[kk] = bi2 - br3; + tr = IRT2*(br5-bi5); + ti = IRT2*(br5+bi5); + cr4[kk] = br4 + tr; + ci4[kk] = bi4 + ti; + cr5[kk] = br4 - tr; + ci5[kk] = bi4 - ti; + tr = -IRT2*(br7+bi7); + ti = IRT2*(br7-bi7); + cr6[kk] = br6 + tr; + ci6[kk] = bi6 + ti; + cr7[kk] = br6 - tr; + ci7[kk] = bi6 - ti; + } + } + } +} + + +/* + * FFT842 (Name kept from the original Fortran version) + * This routine replaces the input DCOMPLEX vector by its + * finite discrete complex fourier transform if in==FFT_FORWARD. + * It replaces the input DCOMPLEX vector by its finite discrete + * complex inverse fourier transform if in==FFT_INVERSE. + * + * The implementation is a radix-2 FFT, but with faster shortcuts for + * radix-4 and radix-8. It performs as many radix-8 iterations as + * possible, and then finishes with a radix-2 or -4 iteration if needed. + */ +void FFT842(int direction, int n, DCOMPLEX *b) +/* direction: FFT_FORWARD or FFT_INVERSE + * n: length of vector + * *b: input vector */ +{ + double fn, r, fi; + + int L[16],L1,L2,L3,L4,L5,L6,L7,L8,L9,L10,L11,L12,L13,L14,L15; +/* int j0,j1,j2,j3,j4,j5,j6,j7,j8,j9,j10,j11,j12,j13,j14;*/ + int j1,j2,j3,j4,j5,j6,j7,j8,j9,j10,j11,j12,j13,j14; + int i, j, ij, ji, ij1, ji1; +/* int nn, n2pow, n8pow, nthpo, ipass, nxtlt, length;*/ + int n2pow, n8pow, nthpo, ipass, nxtlt, length; + + n2pow = fastlog2(n); + nthpo = n; + fn = 1.0 / (double)nthpo; /* Scaling factor for inverse transform */ + + if(direction==FFT_FORWARD) + /* Conjugate the input */ + for(i=0;icols ? rows : cols; + errflag = allocateBuffer(maxsize); + if(errflag != NO_ERROR) return(errflag); + + /* Compute transform row by row */ + if(cols>1) + for(i=0;i1) + for(i=0;icols ? rows : cols; +// errflag = allocateBuffer(maxsize); +// if(errflag != NO_ERROR) return(errflag); +// +// /* Compute transform row by row */ +// if(cols>1) +// for(i=0;i1) +// for(i=0;i. +-----------------------------------------------------------------------------------------*/ + +/* + *------------------------------------------------------------- + * + * 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; + 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; + 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); + + w=vector(1,ns/2); + kx=vector(1,ntr_glob/2); + + 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]); + + angi = ANGTAPI * PI/180.0; + ango = ANGTAPO * PI/180.0; + +// 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); + + /* 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]); + 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++) { + 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; + } + } + } + } + + 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)]; + } + + } + } + + /* 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); + 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); + fclose(file); + + /* write spectrum to file */ + 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); + fclose(file); + } + + /* 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])); + } + } + + fprintf(FP,"\n Testposition 4\n "); + + /* reverse 2d fft */ + fft2(data_pup, data_pupim, ns, ntr_glob,-1); + + /* save pup*/ + if (sw == 1){ + sprintf(pf,"%s_pup.su.shot%i.it%i",SEIS_FILE,ishot,iter); + } + if (sw == 2){ /* for step length calculation*/ + sprintf(pf,"%s_pup.su.step.shot%i.it%i",SEIS_FILE,ishot,iter); + } + outseis_glob(fp,fopen(pf,"w"), 0, data_pup,recpos,recpos_loc,ntr,srcpos,1,ns,SEIS_FORMAT,ishot,1); + + /* overwrite p data */ + for (k=1;k<=ns;k++) { + for (i=1;i<=ntr_glob;i++) { + data_p[k][i]=data_pup[k][i]; + } + } + + fprintf(FP,"\n Testposition 5\n "); + + free_matrix(filtwk, 1,ns/2,1,ntr_glob/2); + free_matrix(filtwk_full, 1,ns,1,ntr_glob); + free_matrix(data_pim, 1,ns,1,ntr_glob); + free_matrix(data_vyim, 1,ns,1,ntr_glob); + free_matrix(data_pup, 1,ns,1,ntr_glob); + free_matrix(data_pupim, 1,ns,1,ntr_glob); + free_vector(w,1,ns/2); + free_vector(kx,1,ntr_glob/2); +} \ No newline at end of file diff --git a/src/read_par_json.c b/src/read_par_json.c index 3e31fc9..1c843f8 100644 --- a/src/read_par_json.c +++ b/src/read_par_json.c @@ -127,6 +127,9 @@ void read_par_json(FILE *fp, char *fileinp){ extern int WOLFE_TRY_OLD_STEPLENGTH; extern float WOLFE_C1_SL; extern float WOLFE_C2_SL; + + extern int WAVESEP; + extern float VEL, DENS, ANGTAPI, ANGTAPO; /* definition of local variables */ @@ -328,6 +331,31 @@ void read_par_json(FILE *fp, char *fileinp){ if (get_int_from_objectlist("IDY",number_readobjects,&IDY,varname_list, value_list)){ IDY=1; fprintf(fp,"Variable IDY is set to default value %d.\n",IDY);} + + + /*================================= + section wavefield separation + =================================*/ + if (get_int_from_objectlist("WAVESEP",number_readobjects,&WAVESEP,varname_list, value_list)){ + WAVESEP=0; + fprintf(fp,"Variable WAVESEP is set to default value %d.\n",WAVESEP);} + else { + if (WAVESEP>0) { + 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);} + 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);} + 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);} + 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);} + } + } + /*================================= section seismogramm parameters diff --git a/src/seismo_ssg.c b/src/seismo_ssg.c index b3fe083..8851884 100644 --- a/src/seismo_ssg.c +++ b/src/seismo_ssg.c @@ -27,7 +27,7 @@ void seismo_ssg(int lsamp, int ntr, int **recpos, float **sectionvx, float **sectionvy,float **sectionvz, float **sectionp, float **sectioncurl, float **sectiondiv, float **vx, float **vy, float **vz, float **sxx, float **syy, float **sp, float **pi, float **u, float *hc){ - extern int NDT, SEISMO, FDORDER, ACOUSTIC,WAVETYPE; + extern int NDT, SEISMO, FDORDER, ACOUSTIC,WAVETYPE,WAVESEP; extern float DH, DT; extern FILE *FP; int i,j, itr, ins, nxrec, nyrec, m, fdoh; @@ -129,7 +129,7 @@ void seismo_ssg(int lsamp, int ntr, int **recpos, float **sectionvx, else sectionp[itr][ins]=-sp[nyrec][nxrec]; } - if (WAVETYPE==2 || WAVETYPE==3) { + if (WAVETYPE==2 || WAVETYPE==3 || WAVESEP==1) { sectionvz[itr][ins]=vz[nyrec][nxrec]; } break; diff --git a/src/util.c b/src/util.c index cd04a69..1bb0b1e 100644 --- a/src/util.c +++ b/src/util.c @@ -102,6 +102,21 @@ float *vector(int ni, int nj){ return a-ni+SHIFT_IND; } +COMPLEX *cplxvector(int nstart, int nend){ + /* allocate a complex vector with subscript range v[nstart..nend] and initializing + this vector, eg. vector[nstart..nend]=0.0 */ + COMPLEX *v; + int i; + + v=(COMPLEX *)malloc((size_t) ((nend-nstart+1+SHIFT_IND)*sizeof(COMPLEX))); + if (!v) err("allocation failure in function cplxvector()"); + for (i=0;i<(nend-nstart+1+SHIFT_IND);i++){ + v[i].re=0.0; + v[i].im=0.0; + } + return v-nstart+SHIFT_IND; +} + int *ivector(int ni, int nj){ @@ -349,6 +364,11 @@ void free_vector(float *a, int ni, int nj){ free((FREE_ARGUMENT) (a+ni-SHIFT_IND)); } +void free_cplxvector(COMPLEX *v, int nstart, int nend){ + /* free a complex vector allocated with vector() */ + free((char*) (v+nstart-SHIFT_IND)); +} + void free_ivector(int *a, int ni, int nj){ free((FREE_ARGUMENT) (a+ni-SHIFT_IND)); } -- GitLab