### Merge branch 'feature/STF_external_wavelet' into develop

parents 9be0a542 002d49f3
 ... ... @@ -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 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} 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} \end{equation} Fuchs-M\"uller wavelet (SOURCE\_SHAPE=2): SOURCE\_SHAPE=2, Fuchs-M\"uller wavelet: \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 \label{eq_fm} \end{equation} $sin^3$ wavelet (SOURCE\_SHAPE=4): SOURCE\_SHAPE=4, $sin^3$ wavelet: \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 \label{eq_s3} \end{equation} First derivative of a Gaussian function (SOURCE\_SHAPE=5): SOURCE\_SHAPE=5, First derivative of a Gaussian function: \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 \label{eq_deriv_of_gaussian} \end{equation} Delta pulse (SOURCE\_SHAPE=6): 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\_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.\\ 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{center} ... ... @@ -162,8 +163,9 @@ spectrum, c) phase spectrum. } \label{fig_source_wavelets_json} \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 your source wavelet in ASCII-format in SIGNAL\_FILE. SIGNAL\_FILE should contain one sample per line. It should thus look like: \newpage 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} 0.0 ... ... @@ -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. \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. 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.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\_start.su and the inverted wavelets will be stored in SIGNAL\_FILE.shot.su (or SIGNAL\_FILE...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 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 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: {\color{blue}{\begin{verbatim} ... ... @@ -620,6 +627,7 @@ smaller than PRO the inversion aborts or in case of using frequency filtering (T \newpage \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).\\ {\color{blue}{\begin{verbatim} ... ...
 ... ... @@ -1187,7 +1187,7 @@ int main(int argc, char **argv){ 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); ... ... @@ -1197,12 +1197,12 @@ int main(int argc, char **argv){ /* calculate wavelet for each source point P SV */ if(WAVETYPE==1||WAVETYPE==3){ 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 */ if(WAVETYPE==2||WAVETYPE==3){ 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){ /* calculate wavelet for each source point P SV */ if(WAVETYPE==1||WAVETYPE==3){ 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 */ if(WAVETYPE==2||WAVETYPE==3){ 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){ /* calculate wavelet for each source point P SV */ if(WAVETYPE==1||WAVETYPE==3){ 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 */ if(WAVETYPE==2||WAVETYPE==3){ 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, 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); ... ... @@ -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 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_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, float ** psxy, float ** prho, 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); void writebufs(float ** sxx, float ** syy, ... ...
 ... ... @@ -22,7 +22,7 @@ #include "fd.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 */ extern int MYID; ... ... @@ -38,17 +38,25 @@ void inseis_source_wavelet(float *section, int ns, int ishot, int SH){ char data[STRING_SIZE]; FILE *fpdata; if(SH==0) { if(USE_WORKFLOW){ sprintf(data,"%s.stage%d.shot%d.su",SIGNAL_FILE,WORKFLOW_STAGE,ishot); if (STF==0){ /* reading inverted signals */ if(SH==0) { 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 { 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 { if(USE_WORKFLOW){ sprintf(data,"%s.stage%d.shot%d.su",SIGNAL_FILE_SH,WORKFLOW_STAGE,ishot); } else { /* reading signals for STF inversion */ if(SH==0) { sprintf(data,"%s.shot%d_start.su",SIGNAL_FILE,ishot); } 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 /* declaration of global variables */ 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 int TRKILL_STF, NORMALIZE, USE_WORKFLOW, WORKFLOW_STAGE; extern char TRKILL_FILE_STF[STRING_SIZE]; extern char SIGNAL_FILE[STRING_SIZE]; extern char SIGNAL_FILE_SH[STRING_SIZE]; extern FILE *FP; extern int TRKILL_STF_OFFSET; extern float TRKILL_STF_OFFSET_LOWER; ... ... @@ -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==7){ inseis_source_wavelet(psource,ns,ishot,SH); inseis_source_wavelet(psource,ns,ishot,SH,1); } /* calculating wavelet SIN**3 for convoling with STF */ tshift=srcpos[ishot]; fc=srcpos[ishot]; ts=1.0/fc; for (nt=1;nt<=ns;nt++){ t=(float)nt*DT; switch (SOURCE_SHAPE){ case 1 : /* Old Ricker Wavelet */ /* tau=PI*(t-ts-tshift)/(1.5*ts); amp=(((1.0-4.0*tau*tau)*exp(-2.0*tau*tau))); */ /* New Ricker Wavelet, equal to SOFI2D */ tau=PI*(t-1.5*ts-tshift)/(ts); amp=(((1.0-2.0*tau*tau)*exp(-tau*tau))); break; case 2 : if ((t(tshift+ts))) amp=0.0; else amp=((sin(2.0*PI*(t-tshift)*fc) -0.5*sin(4.0*PI*(t-tshift)*fc))); /*amp=((sin(2.0*PI*(t+tshift)*fc) -0.5*sin(4.0*PI*(t+tshift)*fc)));*/ break; case 3 : if (nt<=nts) amp=psource[nt]; else amp=0.0; // amp=psource[nt]; break; /* source wavelet from file SOURCE_FILE */ case 4 : /*tau=PI*(t-ts-tshift)/(1.5*ts);*/ /* Ricker */ /*amp=((t-ts-tshift)*exp(-2.0*tau*tau));*/ if ((t(tshift+ts))) amp=0.0; else amp=pow(sin(PI*(t+tshift)/ts),3.0); break; /* sinus raised to the power of three */ break; case 5 : /* first derivative of a Gaussian */ ts=1.2/fc; ag = PI*PI*fc*fc; amp = - 2.0 * ag * (t-ts) * exp(-ag*(t-ts)*(t-ts)); break; case 6 : /* Bandlimited Spike */ amp=0.0; if(nt==1+iround(tshift/DT)){ amp = 1.0;} break; case 7 : /* source wavelet from file SOURCE_FILE */ amp=psource[nt]; break; case 8 : /* integral of sinus raised to the power of three */ if (t=tshift) && (t<=(tshift+ts))){ amp=(ts/(0.75*PI))*(0.5-0.75*cos(PI*(t-tshift)/ts)+0.25*pow(cos(PI*(t-tshift)/ts),3.0));} if (t>(tshift+ts)) {amp=ts/(0.75*PI);} break; default : declare_error("Which source-wavelet ? "); }/* end of switch (SOURCE_SHAPE) */ wavelet[nt]=amp; }/* end of for (nt=1;nt<=ns;nt++) */ for (nt=1;nt<=ns;nt++){ t=(float)nt*DT; switch (SOURCE_SHAPE){ case 1 : /* New Ricker Wavelet, equal to SOFI2D */ tau=PI*(t-1.5*ts-tshift)/(ts); amp=(((1.0-2.0*tau*tau)*exp(-tau*tau))); break; case 2 : if ((t(tshift+ts))) amp=0.0; else amp=((sin(2.0*PI*(t-tshift)*fc) -0.5*sin(4.0*PI*(t-tshift)*fc))); break; case 3 : /* source wavelet from file SOURCE_FILE */ if (nt<=nts) amp=psource[nt]; else amp=0.0; break; case 4 : /* sinus raised to the power of three */ if ((t(tshift+ts))) amp=0.0; else amp=pow(sin(PI*(t+tshift)/ts),3.0); break; break; case 5 : /* first derivative of a Gaussian */ ts=1.2/fc; ag = PI*PI*fc*fc; amp = - 2.0 * ag * (t-ts) * exp(-ag*(t-ts)*(t-ts)); break; case 6 : /* Bandlimited Spike */ amp=0.0; if(nt==1+iround(tshift/DT)){ amp = 1.0;} break; case 7 : /* source wavelet from file SOURCE_FILE */ amp=psource[nt]; break; case 8 : /* integral of sinus raised to the power of three */ if (t=tshift) && (t<=(tshift+ts))){ amp=(ts/(0.75*PI))*(0.5-0.75*cos(PI*(t-tshift)/ts)+0.25*pow(cos(PI*(t-tshift)/ts),3.0));} if (t>(tshift+ts)) {amp=ts/(0.75*PI);} 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 */ 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 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); /*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(); free(data.triples);*/ ... ...
 ... ... @@ -24,10 +24,10 @@ #include "segy.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 order : order of filter ntr : number of traces ... ... @@ -36,40 +36,82 @@ void timedomain_filt_vector(float * data, float fc, int order, int ntr, int ns, 1: lowpass 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 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 */ extern float DT; /* declaration of local variables */ int itr, j; double *seismogram, T0; seismogram=dvector(1,ns); T0=1.0/fc; if ((order%2)!=0){ declare_error("Order of timedomain filter must be an even number!");} for (j=1;j<=ns;j++){ seismogram[j]=(double)data[j];} if (method==1){ /*lowpass filter*/ seife_lpb(seismogram,ns,DT,T0,order);} if (method==2){ /*highpass filter*/ seife_hpb(seismogram,ns,DT,T0,order);} for (j=1;j<=ns;j++){ data[j]=(float)seismogram[j];} free_dvector(seismogram,1,ns); } for (j=1;j<=ns;j++){ data[j]=(float)seismogram[j]; } } if ((method==2)||(F_HP)){ /*highpass filter*/ for (j=1;j<=ns;j++){ seismogram_hp[j]=(double)data[j]; } seife_hpb(seismogram_hp,ns+1,DT,T0_hp,order); if (ZERO_PHASE==1){ ns_reverse=ns; 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); ns_reverse=ns; for (j=1;j<=ns;j++) { seismogram_hp[ns_reverse]=seismogram_reverse_hp[j]; ns_reverse--;} } for (j=1;j<=ns;j++){ data[j]=(float)seismogram_hp[j]; } } free_dvector(seismogram,1,ns); if (ZERO_PHASE==1) free_dvector(seismogram_reverse,1,ns); free_dvector(seismogram_hp,1,ns); if (ZERO_PHASE==1) free_dvector(seismogram_reverse_hp,1,ns); } \ No newline at end of file
 ... ... @@ -26,7 +26,7 @@ #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 */ ... ... @@ -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==7){ 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); ... ... @@ -144,7 +144,7 @@ float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH){ if (SOURCE_SHAPE_SH==7){ 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); ... ...
 ... ... @@ -153,7 +153,7 @@ void write_par(FILE *fp){ fprintf(fp," Bandlimited Spike \n"); break; 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; case 8 : 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!