From 95a331579e7fca33c3fb6a3f88036713639d7be9 Mon Sep 17 00:00:00 2001 From: Florian Wittkamp Date: Wed, 9 Mar 2016 15:38:16 +0100 Subject: [PATCH] BUGFIX: STF and SH waves In case of the SH waves, the STF wavelet was NOT used. To calculate the synthetic data, a new synthetic wavelet was calculated. --- src/IFOS2D.c | 33 +++++++++++++++++++++++---------- src/wavelet.c | 45 +++++++++++++++------------------------------ 2 files changed, 38 insertions(+), 40 deletions(-) diff --git a/src/IFOS2D.c b/src/IFOS2D.c index d779c84..9407459 100644 --- a/src/IFOS2D.c +++ b/src/IFOS2D.c @@ -95,9 +95,9 @@ int main(int argc, char **argv){ int step1, step2, step3=0, itests, iteste, stepmax, countstep; float scalefac; - /* Variables for Pseudo-Hessian calculation */ int RECINC, ntr1; - int SOURCE_SHAPE_OLD; + int SOURCE_SHAPE_OLD=0; + int SOURCE_SHAPE_OLD_SH=0; /* Variables for L-BFGS */ int LBFGS_NPAR=3; @@ -903,11 +903,6 @@ int main(int argc, char **argv){ MPI_Barrier(MPI_COMM_WORLD); - /* comunication initialisation for persistent communication */ - /*comm_ini(bufferlef_to_rig, bufferrig_to_lef, buffertop_to_bot, bufferbot_to_top, req_send, req_rec);*/ - - snapseis=1; - snapseis1=1; SHOTINC=1; RECINC=1; @@ -917,7 +912,9 @@ int main(int argc, char **argv){ case 2: FC_EXT=filter_frequencies(&nfrq); FC=FC_EXT[FREQ_NR]; break; } + /* Save old SOURCE_SHAPE, which is needed for STF */ SOURCE_SHAPE_OLD = SOURCE_SHAPE; + if(WAVETYPE==2 || WAVETYPE==3) SOURCE_SHAPE_OLD_SH=SOURCE_SHAPE_SH; nt_out=10000; if(!VERBOSE) nt_out=1e5; @@ -1163,7 +1160,8 @@ int main(int argc, char **argv){ for (ishot=1;ishot<=nshots;ishot+=SHOTINC){ SOURCE_SHAPE = SOURCE_SHAPE_OLD; - + if(WAVETYPE==2 || WAVETYPE==3) SOURCE_SHAPE_SH=SOURCE_SHAPE_OLD_SH; + /*------------------------------------------------------------------------------*/ /*----------- Start of inversion of source time function -----------------------*/ /*------------------------------------------------------------------------------*/ @@ -1558,10 +1556,19 @@ int main(int argc, char **argv){ srcpos_loc = splitsrc(srcpos,&nsrc_loc, nsrc); } + /*-------------------*/ + /* Use STF wavelet */ + /*-------------------*/ if(INV_STF){ + SOURCE_SHAPE=7; - if(WAVETYPE==1||WAVETYPE==3) fprintf(FP,"\n Using optimized source time function located in %s.shot%d \n",SIGNAL_FILE,ishot); - if(WAVETYPE==2||WAVETYPE==3) fprintf(FP,"\n Using optimized source time function located in %s.shot%d \n",SIGNAL_FILE_SH,ishot); + if( WAVETYPE==1 || WAVETYPE==3 ) fprintf(FP,"\n Using optimized source time function located in %s.shot%d \n",SIGNAL_FILE,ishot); + + if( WAVETYPE==2 || WAVETYPE==3 ) { + SOURCE_SHAPE_SH=7; + fprintf(FP,"\n Using optimized source time function located in %s.shot%d \n",SIGNAL_FILE_SH,ishot); + } + } MPI_Barrier(MPI_COMM_WORLD); @@ -2043,6 +2050,9 @@ int main(int argc, char **argv){ } /* end ADJOINT_TYPE */ } + /* --------------------------------- */ + /* read seismic data from SU file vz */ + /* --------------------------------- */ if(WAVETYPE==2 || WAVETYPE==3){ inseis(fprec,ishot,sectionread,ntr_glob,ns,10,iter); if ((TIME_FILT==1 )|| (TIME_FILT==2)){ @@ -3526,6 +3536,9 @@ int main(int argc, char **argv){ } } + /* --------------------------------- */ + /* read seismic data from SU file vz */ + /* --------------------------------- */ if(WAVETYPE==2 || WAVETYPE==3){ inseis(fprec,ishot,sectionread,ntr_glob,ns,10,iter); if ((TIME_FILT==1 )|| (TIME_FILT==2)){ diff --git a/src/wavelet.c b/src/wavelet.c index f1dcf8d..f69bee1 100644 --- a/src/wavelet.c +++ b/src/wavelet.c @@ -62,14 +62,9 @@ float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH){ fc=srcpos_loc[5][k]; a=srcpos_loc[6][k]; ts=1.0/fc; - ag = PI*PI*fc*fc; 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))); @@ -78,22 +73,17 @@ float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH){ 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 : + /* source wavelet from file SOURCE_FILE */ if (nt<=nts) amp=psource[nt]; else amp=0.0; - // amp=psource[nt]; - break; /* source wavelet from file SOURCE_FILE */ + break; case 4 : - /*tau=PI*(t-ts-tshift)/(1.5*ts);*/ /* Ricker */ - /*amp=((t-ts-tshift)*exp(-2.0*tau*tau));*/ - + /* 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; /* sinus raised to the power of three */ + break; break; case 5 : @@ -147,10 +137,15 @@ float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH){ /* ---------------------------- */ if (SH==1) { - if (SOURCE_SHAPE_SH==3) psource=rd_sour(&nts,fopen(SIGNAL_FILE_SH,"r")); + + if (SOURCE_SHAPE_SH==3) { + psource=rd_sour(&nts,fopen(SIGNAL_FILE_SH,"r")); + } + if (SOURCE_SHAPE_SH==7){ psource=vector(1,NT); - inseis_source_wavelet(psource,NT,ishot,SH);} + inseis_source_wavelet(psource,NT,ishot,SH); + } signals=fmatrix(1,nsrc,1,NT); @@ -162,14 +157,9 @@ float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH){ fc=srcpos_loc[5][k]; a=srcpos_loc[6][k]; ts=1.0/fc; - ag = PI*PI*fc*fc; switch (SOURCE_SHAPE_SH){ 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))); @@ -178,22 +168,17 @@ float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH){ 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 : + /* source wavelet from file SOURCE_FILE */ if (nt<=nts) amp=psource[nt]; else amp=0.0; - // amp=psource[nt]; - break; /* source wavelet from file SOURCE_FILE */ + break; case 4 : - /*tau=PI*(t-ts-tshift)/(1.5*ts);*/ /* Ricker */ - /*amp=((t-ts-tshift)*exp(-2.0*tau*tau));*/ - + /* 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; /* sinus raised to the power of three */ + break; break; case 5 : -- GitLab