Commit 002d49f3 authored by laura.gassner's avatar laura.gassner Committed by Florian Wittkamp

Adaption of external wavelet with STF

# Conflicts:
#	src/IFOS2D.c
#	src/wavelet.c
parent 9be0a542
...@@ -120,37 +120,38 @@ Default values are: ...@@ -120,37 +120,38 @@ Default values are:
Five built-in wavelets of the seismic source are available. The corresponding source time functions are defined in \texttt{src/wavelet.c}. You may modify the time functions in this file and recompile to include your Five built-in wavelets of the seismic source are available. The corresponding source time functions are defined in \texttt{src/wavelet.c}. You may modify the time functions in this file and recompile to include your
own analytical wavelet or to modify the shape of the built-in wavelets. own analytical wavelet or to modify the shape of the built-in wavelets.
\newline
Ricker wavelet (SOURCE\_SHAPE=1): SOURCE\_SHAPE=1, Ricker wavelet:
\begin{equation} \begin{equation}
r(\tau)=\left(1-2\tau^2\right)\exp(-\tau^2) \quad \mbox{with} \quad \tau=\frac{\pi(t-1.5/f_c-t_d)}{1.0/f_c} r(\tau)=\left(1-2\tau^2\right)\exp(-\tau^2) \quad \mbox{with} \quad \tau=\frac{\pi(t-1.5/f_c-t_d)}{1.0/f_c}
\label{eq_ricker} \label{eq_ricker}
\end{equation} \end{equation}
Fuchs-M\"uller wavelet (SOURCE\_SHAPE=2): SOURCE\_SHAPE=2, Fuchs-M\"uller wavelet:
\begin{equation} \begin{equation}
f_m(t)=\sin(2\pi(t-t_d)f_c)-0.5\sin(4\pi(t-t_d)f_c) \quad \mbox{if} \quad t\in[t_d,t_d+1/fc] \quad \mbox{else} \quad fm(t)=0 f_m(t)=\sin(2\pi(t-t_d)f_c)-0.5\sin(4\pi(t-t_d)f_c) \quad \mbox{if} \quad t\in[t_d,t_d+1/fc] \quad \mbox{else} \quad fm(t)=0
\label{eq_fm} \label{eq_fm}
\end{equation} \end{equation}
$sin^3$ wavelet (SOURCE\_SHAPE=4): SOURCE\_SHAPE=4, $sin^3$ wavelet:
\begin{equation} \begin{equation}
s3(t)=0.75 \pi f_c \sin(\pi(t+t_d)f_c)^3\quad \mbox{if} \quad t \in[t_d,t_d+1/fc] \quad \mbox{else} \quad s3(t)=0 s3(t)=0.75 \pi f_c \sin(\pi(t+t_d)f_c)^3\quad \mbox{if} \quad t \in[t_d,t_d+1/fc] \quad \mbox{else} \quad s3(t)=0
\label{eq_s3} \label{eq_s3}
\end{equation} \end{equation}
First derivative of a Gaussian function (SOURCE\_SHAPE=5): SOURCE\_SHAPE=5, First derivative of a Gaussian function:
\begin{equation} \begin{equation}
f(t)= -2.0 a (t-t_s) \exp(-a (t-t_s)^2)\quad \mbox{with} \quad a=\pi^2 f_c^2 \quad \mbox{and} \quad t_s=1.2/f_c f(t)= -2.0 a (t-t_s) \exp(-a (t-t_s)^2)\quad \mbox{with} \quad a=\pi^2 f_c^2 \quad \mbox{and} \quad t_s=1.2/f_c
\label{eq_deriv_of_gaussian} \label{eq_deriv_of_gaussian}
\end{equation} \end{equation}
Delta pulse (SOURCE\_SHAPE=6): SOURCE\_SHAPE=6, delta pulse: Lowpass filtered delta pulse. Note, that it is not clear if the lowpass filter used in the current version works correctly for a delta pulse.\\
Lowpass filtered delta pulse. Note, that it is not clear if the lowpass filter used in the current version works correctly for a delta pulse.\\
Source time function from SIGNAL\_FILE in su format (SOURCE\_SHAPE=7).\\ % Source time function from SIGNAL\_FILE in su format (SOURCE\_SHAPE=7).\\
In these equations, t denotes time and $f_c=1/TS$ is the center frequency. $t_d$ is a time delay which can be defined for each source position in SOURCE\_FILE. Note that the symmetric (zero phase) Ricker signal is always delayed by $1.0/f_c$, which means that after one period the maximum amplitude is excited at the source location. Three of these 5 source wavelets and the corresponding amplitude spectra for a center frequency of $f_c=50$ Hz and a delay of $t_d=0$ are plotted in Figure \ref{fig_source_wavelets_json}. Note the delay of the Ricker signal described above. The Fuchs-M\"uller wavelet has a slightly higher center frequency and covers a broader frequency range. In these equations, t denotes time and $f_c$ is the center frequency. $t_d$ is a time delay which can be defined for each source position. Note that the symmetric (zero phase) Ricker signal is always delayed by $1.0/f_c$, which means that after one period the maximum amplitude is excited at the source location. Three of these 5 source wavelets and the corresponding amplitude spectra for a center frequency of $f_c=50$ Hz and $t_d=0$ are plotted in Figure \ref{fig_source_wavelets_json}. Note the delay of the Ricker signal described above. The Fuchs-M\"uller wavelet has a slightly higher center frequency and covers a broader frequency range.
\newline
\begin{figure} \begin{figure}
\begin{center} \begin{center}
...@@ -162,8 +163,9 @@ spectrum, c) phase spectrum. } ...@@ -162,8 +163,9 @@ spectrum, c) phase spectrum. }
\label{fig_source_wavelets_json} \label{fig_source_wavelets_json}
\end{figure} \end{figure}
You may also use your own time function as the source wavelet (for instance the signal of the first arrival recorded by a geophone at near offsets). Specify SOURCE\_SHAPE=3 and save the samples of \newpage
your source wavelet in ASCII-format in SIGNAL\_FILE. SIGNAL\_FILE should contain one sample per line. It should thus look like:
SOURCE\_SHAPE=3 allows you to use your own time function as the source wavelet stored in ASCII-format in SIGNAL\_FILE. SIGNAL\_FILE should then contain one sample per line. It should thus look like:
{\color{blue}{\begin{verbatim} {\color{blue}{\begin{verbatim}
0.0 0.0
...@@ -175,12 +177,17 @@ your source wavelet in ASCII-format in SIGNAL\_FILE. SIGNAL\_FILE should contain ...@@ -175,12 +177,17 @@ your source wavelet in ASCII-format in SIGNAL\_FILE. SIGNAL\_FILE should contain
The time interval between the samples must equal the time step interval (DT) of the FD simulation (see above)! Therefore it might be necessary to resample/interpolate a given source time function with a smaller sample rate. You may use the matlab script mfiles/resamp.m to resample your external source signal to the required sampling interval. The time interval between the samples must equal the time step interval (DT) of the FD simulation (see above)! Therefore it might be necessary to resample/interpolate a given source time function with a smaller sample rate. You may use the matlab script mfiles/resamp.m to resample your external source signal to the required sampling interval.
\newline \newline
It is also possible to read different external source wavelets for each shot. Specify SOURCE\_SHAPE=7 and save the wavelets in su format in SIGNAL\_FILE.shot<no\_of\_shot>. The wavelets in each su file must equal the time step intervel (DT) and the number of time steps of the FD simulation! SOURCE\_SHAPE=7 is used for reading different external source wavelets for each shot. The wavelets in SU-format need to be saved in SIGNAL\_FILE.shot<shotnumber>.su.
If you want to use the source time function inversion (INV\_STF==1, section \ref{sec:STF}) with an external wavelet, this wavelet needs to be provided at SIGNAL\_FILE.shot<shotnumber>\_start.su and the inverted wavelets will be stored in SIGNAL\_FILE.shot<shotnumber>.su (or SIGNAL\_FILE.<workflowstage>.<shotnumber>.su if you use the WORKFLOW option).
The wavelets in each su file must have the same number of samples as specified by TIME/DT!
\newline \newline
The following source types are availabe: explosive sources that excite compressional waves only (SOURCE\_TYPE=1), and point forces in the x- and y-direction (SOURCE\_TYPE=2,3). The following source types are availabe: explosive sources that excite compressional waves only (SOURCE\_TYPE=1), and point forces in the x- and y-direction (SOURCE\_TYPE=2,3).
The force sources excite both P- and S-waves. The explosive source is located at the same position as the diagonal elements of the stress tensor, i.e. at (i,j) (Figure \ref{fig_cell}). The force sources excite both P- and S-waves. The explosive source is located at the same position as the diagonal elements of the stress tensor, i.e. at (i,j) (Figure \ref{fig_cell}).
The forces are located at the same position as the corresponding components of particle velocity (Figure \ref{fig_cell}). If (x,y) denotes the position at which the source location is defined in source.dat, then the actual force in x-direction is located at (x+DX/2,y) and the actual force in y-direction is located at (x,y+DY/2). With SOURCE\_TYPE=4 a custom directive force can be defined by a force angle between y and x. The angle of the force must be specified in the SOURCE\_FILE after AMP. This force is not aligned along the main directions. The forces are located at the same position as the corresponding components of particle velocity (Figure \ref{fig_cell}). If (x,y) denotes the position at which the source location is defined in source.dat, then the actual force in x-direction is located at (x+DX/2,y) and the actual force in y-direction is located at (x,y+DY/2). With SOURCE\_TYPE=4 a custom directive force can be defined by a force angle between y and x. The angle of the force must be specified in the SOURCE\_FILE after AMP. This force is not aligned along the main directions.
\newline
The locations of multiple sources must be defined in an external ASCII file (SOURCE\_FILE) that has the following format: The locations of multiple sources must be defined in an external ASCII file (SOURCE\_FILE) that has the following format:
{\color{blue}{\begin{verbatim} {\color{blue}{\begin{verbatim}
...@@ -620,6 +627,7 @@ smaller than PRO the inversion aborts or in case of using frequency filtering (T ...@@ -620,6 +627,7 @@ smaller than PRO the inversion aborts or in case of using frequency filtering (T
\newpage \newpage
\section{Source wavelet inversion} \section{Source wavelet inversion}
\label{sec:STF}
To remove the contribution of the unknown source time function (STF) from the waveform residuals, it is necessary to design a filter which minimizes the misfit to the field recordings and raw synthetics. Therefore, a second forward simulation is applied. The first one is done with the wavelet specified in SOURCE\_SHAPE and the second one with the optimized source wavelet saved in SIGNAL\_FILE (see Section~\ref{sec:sources}). This optimized source wavelet is kept constant within N\_STF or within a frequency range (see below).\\ To remove the contribution of the unknown source time function (STF) from the waveform residuals, it is necessary to design a filter which minimizes the misfit to the field recordings and raw synthetics. Therefore, a second forward simulation is applied. The first one is done with the wavelet specified in SOURCE\_SHAPE and the second one with the optimized source wavelet saved in SIGNAL\_FILE (see Section~\ref{sec:sources}). This optimized source wavelet is kept constant within N\_STF or within a frequency range (see below).\\
{\color{blue}{\begin{verbatim} {\color{blue}{\begin{verbatim}
......
...@@ -1187,7 +1187,7 @@ int main(int argc, char **argv){ ...@@ -1187,7 +1187,7 @@ int main(int argc, char **argv){
srcpos_loc = splitsrc(srcpos,&nsrc_loc, nsrc); srcpos_loc = splitsrc(srcpos,&nsrc_loc, nsrc);
} }
if((SOURCE_SHAPE==7)||(SOURCE_SHAPE==3))declare_error("SOURCE_SHAPE==7 or SOURCE_SHAPE==3 isn't possible with INV_STF==1"); if(SOURCE_SHAPE==3) declare_error("SOURCE_SHAPE==3 isn't possible with INV_STF==1");
MPI_Barrier(MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD);
...@@ -1197,12 +1197,12 @@ int main(int argc, char **argv){ ...@@ -1197,12 +1197,12 @@ int main(int argc, char **argv){
/* calculate wavelet for each source point P SV */ /* calculate wavelet for each source point P SV */
if(WAVETYPE==1||WAVETYPE==3){ if(WAVETYPE==1||WAVETYPE==3){
signals=NULL; signals=NULL;
signals=wavelet(srcpos_loc,nsrc_loc,ishot,0); signals=wavelet(srcpos_loc,nsrc_loc,ishot,0,1);
} }
/* calculate wavelet for each source point SH */ /* calculate wavelet for each source point SH */
if(WAVETYPE==2||WAVETYPE==3){ if(WAVETYPE==2||WAVETYPE==3){
signals_SH=NULL; signals_SH=NULL;
signals_SH=wavelet(srcpos_loc,nsrc_loc,ishot,1); signals_SH=wavelet(srcpos_loc,nsrc_loc,ishot,1,1);
} }
...@@ -1579,12 +1579,12 @@ int main(int argc, char **argv){ ...@@ -1579,12 +1579,12 @@ int main(int argc, char **argv){
/* calculate wavelet for each source point P SV */ /* calculate wavelet for each source point P SV */
if(WAVETYPE==1||WAVETYPE==3){ if(WAVETYPE==1||WAVETYPE==3){
signals=NULL; signals=NULL;
signals=wavelet(srcpos_loc,nsrc_loc,ishot,0); signals=wavelet(srcpos_loc,nsrc_loc,ishot,0,0);
} }
/* calculate wavelet for each source point SH */ /* calculate wavelet for each source point SH */
if(WAVETYPE==2||WAVETYPE==3){ if(WAVETYPE==2||WAVETYPE==3){
signals_SH=NULL; signals_SH=NULL;
signals_SH=wavelet(srcpos_loc,nsrc_loc,ishot,1); signals_SH=wavelet(srcpos_loc,nsrc_loc,ishot,1,0);
} }
/*------------------------------------------------------------------------------*/ /*------------------------------------------------------------------------------*/
...@@ -3275,12 +3275,12 @@ int main(int argc, char **argv){ ...@@ -3275,12 +3275,12 @@ int main(int argc, char **argv){
/* calculate wavelet for each source point P SV */ /* calculate wavelet for each source point P SV */
if(WAVETYPE==1||WAVETYPE==3){ if(WAVETYPE==1||WAVETYPE==3){
signals=NULL; signals=NULL;
signals=wavelet(srcpos_loc,nsrc_loc,ishot,0); signals=wavelet(srcpos_loc,nsrc_loc,ishot,0,0);
} }
/* calculate wavelet for each source point SH */ /* calculate wavelet for each source point SH */
if(WAVETYPE==2||WAVETYPE==3){ if(WAVETYPE==2||WAVETYPE==3){
signals_SH=NULL; signals_SH=NULL;
signals_SH=wavelet(srcpos_loc,nsrc_loc,ishot,1); signals_SH=wavelet(srcpos_loc,nsrc_loc,ishot,1,0);
} }
/*------------------------------------------------------------------------------*/ /*------------------------------------------------------------------------------*/
......
...@@ -180,7 +180,7 @@ void outseis_vector(FILE *fp, FILE *fpdata, int comp, float *section, ...@@ -180,7 +180,7 @@ void outseis_vector(FILE *fp, FILE *fpdata, int comp, float *section,
void inseis(FILE *fp, int comp, float **section, int ntr, int ns, int sws, int iter); void inseis(FILE *fp, int comp, float **section, int ntr, int ns, int sws, int iter);
void inseis_source_wavelet(float *section, int ns, int ishot, int SH); void inseis_source_wavelet(float *section, int ns, int ishot, int SH, int STF);
void taper(float *section, int ns, float fc); void taper(float *section, int ns, float fc);
...@@ -280,7 +280,7 @@ void surface_elastic_PML(int ndepth, float ** vx, float ** vy, float ** sxx, flo ...@@ -280,7 +280,7 @@ void surface_elastic_PML(int ndepth, float ** vx, float ** vy, float ** sxx, flo
void surface_PML(int ndepth, float ** vx, float ** vy, float ** sxx, float ** syy, float ** sxy, float ** syz, float ***p, float ***q, float ** ppi, float ** pu, float **prho, float **ptaup, float **ptaus, float *etajm, float *peta, float * hc, float * K_x, float * a_x, float * b_x, float ** psi_vxx, float ** ux, float ** uy, float ** uxy, float ** uyz,float ** sxz,float **uxz); void surface_PML(int ndepth, float ** vx, float ** vy, float ** sxx, float ** syy, float ** sxy, float ** syz, float ***p, float ***q, float ** ppi, float ** pu, float **prho, float **ptaup, float **ptaus, float *etajm, float *peta, float * hc, float * K_x, float * a_x, float * b_x, float ** psi_vxx, float ** ux, float ** uy, float ** uxy, float ** uyz,float ** sxz,float **uxz);
void timedomain_filt(float ** data, float fc, int order, int ntr, int ns, int method); void timedomain_filt(float ** data, float fc, int order, int ntr, int ns, int method);
void timedomain_filt_vector(float * data, float fc, int order, int ntr, int ns, int method); void timedomain_filt_vector(float * data, float fc, int order, int ns, int method);
void time_window(float **sectiondata, int iter, int ntr_glob, int **recpos_loc, int ntr, int ns, int ishot); void time_window(float **sectiondata, int iter, int ntr_glob, int **recpos_loc, int ntr, int ns, int ishot);
void time_window_glob(float **sectiondata, int iter, int ntr_glob, int ns, int ishot); void time_window_glob(float **sectiondata, int iter, int ntr_glob, int ns, int ishot);
...@@ -400,7 +400,7 @@ void update_v_rsg_4th(int nx1, int nx2, int ny1, int ny2, int nt, ...@@ -400,7 +400,7 @@ void update_v_rsg_4th(int nx1, int nx2, int ny1, int ny2, int nt,
float ** psxy, float ** prho, float ** psxy, float ** prho,
float ** srcpos_loc, float ** signals, int nsrc, float ** absorb_coeff); float ** srcpos_loc, float ** signals, int nsrc, float ** absorb_coeff);
float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH); float ** wavelet(float ** srcpos_loc, int nsrc, int ishot, int SH, int STF);
float ** wavelet_stf(int nsrc, int ishot, float ** signals_stf); float ** wavelet_stf(int nsrc, int ishot, float ** signals_stf);
void writebufs(float ** sxx, float ** syy, void writebufs(float ** sxx, float ** syy,
......
...@@ -22,7 +22,7 @@ ...@@ -22,7 +22,7 @@
#include "fd.h" #include "fd.h"
#include "segy.h" #include "segy.h"
void inseis_source_wavelet(float *section, int ns, int ishot, int SH){ void inseis_source_wavelet(float *section, int ns, int ishot, int SH, int STF){
/* declaration of extern variables */ /* declaration of extern variables */
extern int MYID; extern int MYID;
...@@ -38,17 +38,25 @@ void inseis_source_wavelet(float *section, int ns, int ishot, int SH){ ...@@ -38,17 +38,25 @@ void inseis_source_wavelet(float *section, int ns, int ishot, int SH){
char data[STRING_SIZE]; char data[STRING_SIZE];
FILE *fpdata; FILE *fpdata;
if(SH==0) { if (STF==0){ /* reading inverted signals */
if(USE_WORKFLOW){ if(SH==0) {
sprintf(data,"%s.stage%d.shot%d.su",SIGNAL_FILE,WORKFLOW_STAGE,ishot); if(USE_WORKFLOW){
sprintf(data,"%s.stage%d.shot%d.su",SIGNAL_FILE,WORKFLOW_STAGE,ishot);
} else {
sprintf(data,"%s.shot%d.su",SIGNAL_FILE,ishot);
}
} else { } else {
sprintf(data,"%s.shot%d.su",SIGNAL_FILE,ishot); if(USE_WORKFLOW){
sprintf(data,"%s.stage%d.shot%d.su",SIGNAL_FILE_SH,WORKFLOW_STAGE,ishot);
} else {
sprintf(data,"%s.shot%d.su",SIGNAL_FILE_SH,ishot);
}
} }
} else { } else { /* reading signals for STF inversion */
if(USE_WORKFLOW){ if(SH==0) {
sprintf(data,"%s.stage%d.shot%d.su",SIGNAL_FILE_SH,WORKFLOW_STAGE,ishot); sprintf(data,"%s.shot%d_start.su",SIGNAL_FILE,ishot);
} else { } else {
sprintf(data,"%s.shot%d.su",SIGNAL_FILE_SH,ishot); sprintf(data,"%s.shot%d_start.su",SIGNAL_FILE_SH,ishot);
} }
} }
......
...@@ -28,12 +28,13 @@ void stf(FILE *fp, float **sectionvy, float ** sectionvy_obs, float ** sectionvy ...@@ -28,12 +28,13 @@ void stf(FILE *fp, float **sectionvy, float ** sectionvy_obs, float ** sectionvy
/* declaration of global variables */ /* declaration of global variables */
extern float DT, DH; extern float DT, DH;
extern int SEIS_FORMAT, MYID, NT, SOURCE_SHAPE, TIME_FILT, TIMEWIN, TAPER_STF; extern int SEIS_FORMAT, MYID, NT, SOURCE_SHAPE, TIME_FILT, TIMEWIN, TAPER_STF, ORDER;
extern char PARA[STRING_SIZE], DATA_DIR[STRING_SIZE]; extern char PARA[STRING_SIZE], DATA_DIR[STRING_SIZE];
extern int TRKILL_STF, NORMALIZE, USE_WORKFLOW, WORKFLOW_STAGE; extern int TRKILL_STF, NORMALIZE, USE_WORKFLOW, WORKFLOW_STAGE;
extern char TRKILL_FILE_STF[STRING_SIZE]; extern char TRKILL_FILE_STF[STRING_SIZE];
extern char SIGNAL_FILE[STRING_SIZE]; extern char SIGNAL_FILE[STRING_SIZE];
extern char SIGNAL_FILE_SH[STRING_SIZE]; extern char SIGNAL_FILE_SH[STRING_SIZE];
extern FILE *FP;
extern int TRKILL_STF_OFFSET; extern int TRKILL_STF_OFFSET;
extern float TRKILL_STF_OFFSET_LOWER; extern float TRKILL_STF_OFFSET_LOWER;
...@@ -211,81 +212,71 @@ void stf(FILE *fp, float **sectionvy, float ** sectionvy_obs, float ** sectionvy ...@@ -211,81 +212,71 @@ void stf(FILE *fp, float **sectionvy, float ** sectionvy_obs, float ** sectionvy
if (SOURCE_SHAPE==3) psource=rd_sour(&nts,fopen(SIGNAL_FILE,"r")); if (SOURCE_SHAPE==3) psource=rd_sour(&nts,fopen(SIGNAL_FILE,"r"));
if (SOURCE_SHAPE==7){ if (SOURCE_SHAPE==7){
inseis_source_wavelet(psource,ns,ishot,SH); inseis_source_wavelet(psource,ns,ishot,SH,1);
} }
/* calculating wavelet SIN**3 for convoling with STF */ /* calculating wavelet SIN**3 for convoling with STF */
tshift=srcpos[4][ishot]; tshift=srcpos[4][ishot];
fc=srcpos[5][ishot]; fc=srcpos[5][ishot];
ts=1.0/fc; ts=1.0/fc;
for (nt=1;nt<=ns;nt++){ for (nt=1;nt<=ns;nt++){
t=(float)nt*DT; t=(float)nt*DT;
switch (SOURCE_SHAPE){ switch (SOURCE_SHAPE){
case 1 : case 1 :
/* Old Ricker Wavelet */ /* New Ricker Wavelet, equal to SOFI2D */
/* tau=PI*(t-ts-tshift)/(1.5*ts); tau=PI*(t-1.5*ts-tshift)/(ts);
amp=(((1.0-4.0*tau*tau)*exp(-2.0*tau*tau))); */ amp=(((1.0-2.0*tau*tau)*exp(-tau*tau)));
break;
/* New Ricker Wavelet, equal to SOFI2D */ case 2 :
tau=PI*(t-1.5*ts-tshift)/(ts); if ((t<tshift) || (t>(tshift+ts))) amp=0.0;
amp=(((1.0-2.0*tau*tau)*exp(-tau*tau))); else amp=((sin(2.0*PI*(t-tshift)*fc)
break; -0.5*sin(4.0*PI*(t-tshift)*fc)));
case 2 : break;
if ((t<tshift) || (t>(tshift+ts))) amp=0.0; case 3 :
else amp=((sin(2.0*PI*(t-tshift)*fc) /* source wavelet from file SOURCE_FILE */
-0.5*sin(4.0*PI*(t-tshift)*fc))); if (nt<=nts) amp=psource[nt];
else amp=0.0;
/*amp=((sin(2.0*PI*(t+tshift)*fc) break;
-0.5*sin(4.0*PI*(t+tshift)*fc)));*/ case 4 :
break; /* sinus raised to the power of three */
case 3 : if ((t<tshift) || (t>(tshift+ts))) amp=0.0;
if (nt<=nts) amp=psource[nt]; else amp=pow(sin(PI*(t+tshift)/ts),3.0);
else amp=0.0; break;
// amp=psource[nt];
break; /* source wavelet from file SOURCE_FILE */ break;
case 4 : case 5 :
/*tau=PI*(t-ts-tshift)/(1.5*ts);*/ /* Ricker */ /* first derivative of a Gaussian */
/*amp=((t-ts-tshift)*exp(-2.0*tau*tau));*/ ts=1.2/fc;
ag = PI*PI*fc*fc;
if ((t<tshift) || (t>(tshift+ts))) amp=0.0; amp = - 2.0 * ag * (t-ts) * exp(-ag*(t-ts)*(t-ts));
else amp=pow(sin(PI*(t+tshift)/ts),3.0); break;
break; /* sinus raised to the power of three */ case 6 :
/* Bandlimited Spike */
break; amp=0.0;
case 5 : if(nt==1+iround(tshift/DT)){
/* first derivative of a Gaussian */ amp = 1.0;}
ts=1.2/fc; break;
ag = PI*PI*fc*fc; case 7 :
amp = - 2.0 * ag * (t-ts) * exp(-ag*(t-ts)*(t-ts)); /* source wavelet from file SOURCE_FILE */
break; amp=psource[nt];
case 6 : break;
/* Bandlimited Spike */ case 8 :
amp=0.0; /* integral of sinus raised to the power of three */
if(nt==1+iround(tshift/DT)){ if (t<tshift) {
amp = 1.0;} amp=0.0;}
break; if ((t>=tshift) && (t<=(tshift+ts))){
case 7 : amp=(ts/(0.75*PI))*(0.5-0.75*cos(PI*(t-tshift)/ts)+0.25*pow(cos(PI*(t-tshift)/ts),3.0));}
/* source wavelet from file SOURCE_FILE */ if (t>(tshift+ts))
amp=psource[nt]; {amp=ts/(0.75*PI);}
break; break;
case 8 : default :
/* integral of sinus raised to the power of three */ declare_error("Which source-wavelet ? ");
if (t<tshift) {
amp=0.0;}
if ((t>=tshift) && (t<=(tshift+ts))){ }/* end of switch (SOURCE_SHAPE) */
amp=(ts/(0.75*PI))*(0.5-0.75*cos(PI*(t-tshift)/ts)+0.25*pow(cos(PI*(t-tshift)/ts),3.0));} wavelet[nt]=amp;
if (t>(tshift+ts))
{amp=ts/(0.75*PI);} }/* end of for (nt=1;nt<=ns;nt++) */
break;
default :
declare_error("Which source-wavelet ? ");
}/* end of switch (SOURCE_SHAPE) */
wavelet[nt]=amp;
}/* end of for (nt=1;nt<=ns;nt++) */
/* convolving wavelet with STF */ /* convolving wavelet with STF */
conv_FD(wavelet,source_time_function,stf_conv_wavelet,ns); conv_FD(wavelet,source_time_function,stf_conv_wavelet,ns);
...@@ -341,9 +332,6 @@ void stf(FILE *fp, float **sectionvy, float ** sectionvy_obs, float ** sectionvy ...@@ -341,9 +332,6 @@ void stf(FILE *fp, float **sectionvy, float ** sectionvy_obs, float ** sectionvy
printf(" PE %d is writing source time function for shot = %d to\n\t %s \n",MYID,ishot,qw); printf(" PE %d is writing source time function for shot = %d to\n\t %s \n",MYID,ishot,qw);
outseis_vector(fp,fopen(qw,"w"),1,stf_conv_wavelet,recpos,recpos_loc,ntr,srcpos,0,ns,SEIS_FORMAT,ishot,0); outseis_vector(fp,fopen(qw,"w"),1,stf_conv_wavelet,recpos,recpos_loc,ntr,srcpos,0,ns,SEIS_FORMAT,ishot,0);
/*sprintf(conv_y_tmp,"%s.shot%d.forward",SEIS_FILE_VY,ishot);
printf(" PE %d is writing %d seismograms (vy) for shot = %d to\n\t %s \n",MYID,ntr_glob,ishot,conv_y_tmp);
outseis_glob(fp,fopen(conv_y_tmp,"w"),1,sectionvy,recpos,recpos_loc,ntr_glob,srcpos,0,ns,SEIS_FORMAT,ishot,0);*/
/*freestfinvengine(); /*freestfinvengine();
free(data.triples);*/ free(data.triples);*/
......
...@@ -24,10 +24,10 @@ ...@@ -24,10 +24,10 @@
#include "segy.h" #include "segy.h"
#include "cseife.h" #include "cseife.h"
void timedomain_filt_vector(float * data, float fc, int order, int ntr, int ns, int method){ void timedomain_filt_vector(float * data, float fc, int order, int ns, int method){
/* /*
data : 2-dimensional array containing seismograms ( data : 1-dimensional array containing seismograms (
fc : corner frequency in Hz fc : corner frequency in Hz
order : order of filter order : order of filter
ntr : number of traces ntr : number of traces
...@@ -36,40 +36,82 @@ void timedomain_filt_vector(float * data, float fc, int order, int ntr, int ns, ...@@ -36,40 +36,82 @@ void timedomain_filt_vector(float * data, float fc, int order, int ntr, int ns,
1: lowpass filter 1: lowpass filter
2: highpass filter 2: highpass filter
*/ */
/* declaration of external variables */
extern float DT, F_HP;
extern int ZERO_PHASE, NT,MYID;
/* declaration of local variables */
int itr, j, ns_reverse;
double *seismogram, *seismogram_reverse, T0;
double *seismogram_hp, *seismogram_reverse_hp, T0_hp;
if ((order%2)!=0){
declare_error("Order of timedomain filter must be an even number!");}
seismogram = dvector(1,ns);
if (ZERO_PHASE==1) seismogram_reverse = dvector(1,ns);
seismogram_hp = dvector(1,ns);
if (ZERO_PHASE==1) seismogram_reverse_hp = dvector(1,ns);
T0=1.0/(double)fc;
if(F_HP)
T0_hp=1.0/(double)F_HP;
if(method==2)
T0_hp=1.0/(double)fc;
if (method==1){ /* lowpass filter */
for (j=1;j<=ns;j++){
seismogram[j]=(double)data[j];
}
seife_lpb(seismogram,ns+1,DT,T0,order); /* ns+1 because vector[0] is also allocated and otherwise seife_lpb do not filter the last sample */
if (ZERO_PHASE==1){
ns_reverse=ns;
for (j=1;j<=ns;j++) {
seismogram_reverse[ns_reverse]=seismogram[j];
ns_reverse--;}
seife_lpb(seismogram_reverse,ns+1,DT,T0,order);
ns_reverse=ns;
for (j=1;j<=ns;j++) {
seismogram[ns_reverse]=seismogram_reverse[j];
ns_reverse--;}
}
/* declaration of extern variables */ for (j=1;j<=ns;j++){
extern float DT; data[j]=(float)seismogram[j];
}
/* declaration of local variables */ }
int itr, j;
double *seismogram, T0; if ((method==2)||(F_HP)){ /*highpass filter*/
for (j=1;j<=ns;j++){
seismogram_hp[j]=(double)data[j];
}
seismogram=dvector(1,ns);
seife_hpb(seismogram_hp,ns+1,DT,T0_hp,order);
T0=1.0/fc;
if (ZERO_PHASE==1){
if ((order%2)!=0){ ns_reverse=ns;
declare_error("Order of timedomain filter must be an even number!");} for (j=1;j<=ns;j++) {
seismogram_reverse_hp[ns_reverse]=seismogram_hp[j];
ns_reverse--;}
seife_hpb(seismogram_reverse_hp,ns+1,DT,T0_hp,order);
for (j=1;j<=ns;j++){ ns_reverse=ns;
seismogram[j]=(double)data[j];} for (j=1;j<=ns;j++) {
seismogram_hp[ns_reverse]=seismogram_reverse_hp[j];
if (method==1){ /*lowpass filter*/ ns_reverse--;}
seife_lpb(seismogram,ns,DT,T0,order);} }
for (j=1;j<=ns;j++){
data[j]=(float)seismogram_hp[j];
if (method==2){ /*highpass filter*/ }
seife_hpb(seismogram,ns,DT,T0,order);} }
for (j=1;j<=ns;j++){ free_dvector(seismogram,1,ns);
data[j]=(float)seismogram[j];} if (ZERO_PHASE==1) free_dvector(seismogram_reverse,1,ns);
free_dvector(seismogram_hp,1,ns);
free_dvector(seismogram,1,ns); if (ZERO_PHASE==1) free_dvector(seismogram_reverse_hp,1,ns);
}
}
\ No newline at end of file
...@@ -26,7 +26,7 @@ ...@@ -26,7 +26,7 @@
#include "fd.h" #include "fd.h"
float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH){ float ** wavelet(float ** srcpos_loc, int nsrc, int ishot, int SH, int STF){
/* extern variables */ /* extern variables */
...@@ -50,7 +50,7 @@ float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH){ ...@@ -50,7 +50,7 @@ float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH){
if (SOURCE_SHAPE==3) psource=rd_sour(&nts,fopen(SIGNAL_FILE,"r")); if (SOURCE_SHAPE==3) psource=rd_sour(&nts,fopen(SIGNAL_FILE,"r"));
if (SOURCE_SHAPE==7){ if (SOURCE_SHAPE==7){
psource=vector(1,NT); psource=vector(1,NT);
inseis_source_wavelet(psource,NT,ishot,SH);} inseis_source_wavelet(psource,NT,ishot,SH,STF);}
signals=fmatrix(1,nsrc,1,NT); signals=fmatrix(1,nsrc,1,NT);
...@@ -144,7 +144,7 @@ float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH){ ...@@ -144,7 +144,7 @@ float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH){
if (SOURCE_SHAPE_SH==7){ if (SOURCE_SHAPE_SH==7){
psource=vector(1,NT); psource=vector(1,NT);
inseis_source_wavelet(psource,NT,ishot,SH); inseis_source_wavelet(psource,NT,ishot,SH,STF);}
} }
signals=fmatrix(1,nsrc,1,NT); signals=fmatrix(1,nsrc,1,NT);
......
...@@ -153,7 +153,7 @@ void write_par(FILE *fp){ ...@@ -153,7 +153,7 @@ void write_par(FILE *fp){
fprintf(fp," Bandlimited Spike \n"); fprintf(fp," Bandlimited Spike \n");
break; break;
case 7 : case 7 :
fprintf(fp," reading from \n\t %s.shot* in su format (one file for each shot)\n",SIGNAL_FILE); fprintf(fp," reading from \n\t %s.shot*.su in su format (one file for each shot)\n",SIGNAL_FILE);
break; break;
case 8 : case 8 :
fprintf(fp," Integral of sin^3 function\n"); fprintf(fp," Integral of sin^3 function\n");
......
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