Commit 95a33157 authored by Florian Wittkamp's avatar Florian Wittkamp

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.
parent 0f8fe79c
...@@ -95,9 +95,9 @@ int main(int argc, char **argv){ ...@@ -95,9 +95,9 @@ int main(int argc, char **argv){
int step1, step2, step3=0, itests, iteste, stepmax, countstep; int step1, step2, step3=0, itests, iteste, stepmax, countstep;
float scalefac; float scalefac;
/* Variables for Pseudo-Hessian calculation */
int RECINC, ntr1; int RECINC, ntr1;
int SOURCE_SHAPE_OLD; int SOURCE_SHAPE_OLD=0;
int SOURCE_SHAPE_OLD_SH=0;
/* Variables for L-BFGS */ /* Variables for L-BFGS */
int LBFGS_NPAR=3; int LBFGS_NPAR=3;
...@@ -903,11 +903,6 @@ int main(int argc, char **argv){ ...@@ -903,11 +903,6 @@ int main(int argc, char **argv){
MPI_Barrier(MPI_COMM_WORLD); 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; SHOTINC=1;
RECINC=1; RECINC=1;
...@@ -917,7 +912,9 @@ int main(int argc, char **argv){ ...@@ -917,7 +912,9 @@ int main(int argc, char **argv){
case 2: FC_EXT=filter_frequencies(&nfrq); FC=FC_EXT[FREQ_NR]; break; 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; SOURCE_SHAPE_OLD = SOURCE_SHAPE;
if(WAVETYPE==2 || WAVETYPE==3) SOURCE_SHAPE_OLD_SH=SOURCE_SHAPE_SH;
nt_out=10000; nt_out=10000;
if(!VERBOSE) nt_out=1e5; if(!VERBOSE) nt_out=1e5;
...@@ -1163,7 +1160,8 @@ int main(int argc, char **argv){ ...@@ -1163,7 +1160,8 @@ int main(int argc, char **argv){
for (ishot=1;ishot<=nshots;ishot+=SHOTINC){ for (ishot=1;ishot<=nshots;ishot+=SHOTINC){
SOURCE_SHAPE = SOURCE_SHAPE_OLD; SOURCE_SHAPE = SOURCE_SHAPE_OLD;
if(WAVETYPE==2 || WAVETYPE==3) SOURCE_SHAPE_SH=SOURCE_SHAPE_OLD_SH;
/*------------------------------------------------------------------------------*/ /*------------------------------------------------------------------------------*/
/*----------- Start of inversion of source time function -----------------------*/ /*----------- Start of inversion of source time function -----------------------*/
/*------------------------------------------------------------------------------*/ /*------------------------------------------------------------------------------*/
...@@ -1558,10 +1556,19 @@ int main(int argc, char **argv){ ...@@ -1558,10 +1556,19 @@ int main(int argc, char **argv){
srcpos_loc = splitsrc(srcpos,&nsrc_loc, nsrc); srcpos_loc = splitsrc(srcpos,&nsrc_loc, nsrc);
} }
/*-------------------*/
/* Use STF wavelet */
/*-------------------*/
if(INV_STF){ if(INV_STF){
SOURCE_SHAPE=7; 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==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==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); MPI_Barrier(MPI_COMM_WORLD);
...@@ -2043,6 +2050,9 @@ int main(int argc, char **argv){ ...@@ -2043,6 +2050,9 @@ int main(int argc, char **argv){
} /* end ADJOINT_TYPE */ } /* end ADJOINT_TYPE */
} }
/* --------------------------------- */
/* read seismic data from SU file vz */
/* --------------------------------- */
if(WAVETYPE==2 || WAVETYPE==3){ if(WAVETYPE==2 || WAVETYPE==3){
inseis(fprec,ishot,sectionread,ntr_glob,ns,10,iter); inseis(fprec,ishot,sectionread,ntr_glob,ns,10,iter);
if ((TIME_FILT==1 )|| (TIME_FILT==2)){ if ((TIME_FILT==1 )|| (TIME_FILT==2)){
...@@ -3526,6 +3536,9 @@ int main(int argc, char **argv){ ...@@ -3526,6 +3536,9 @@ int main(int argc, char **argv){
} }
} }
/* --------------------------------- */
/* read seismic data from SU file vz */
/* --------------------------------- */
if(WAVETYPE==2 || WAVETYPE==3){ if(WAVETYPE==2 || WAVETYPE==3){
inseis(fprec,ishot,sectionread,ntr_glob,ns,10,iter); inseis(fprec,ishot,sectionread,ntr_glob,ns,10,iter);
if ((TIME_FILT==1 )|| (TIME_FILT==2)){ if ((TIME_FILT==1 )|| (TIME_FILT==2)){
......
...@@ -62,14 +62,9 @@ float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH){ ...@@ -62,14 +62,9 @@ float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH){
fc=srcpos_loc[5][k]; fc=srcpos_loc[5][k];
a=srcpos_loc[6][k]; a=srcpos_loc[6][k];
ts=1.0/fc; ts=1.0/fc;
ag = PI*PI*fc*fc;
switch (SOURCE_SHAPE){ switch (SOURCE_SHAPE){
case 1 : 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 */ /* New Ricker Wavelet, equal to SOFI2D */
tau=PI*(t-1.5*ts-tshift)/(ts); tau=PI*(t-1.5*ts-tshift)/(ts);
amp=(((1.0-2.0*tau*tau)*exp(-tau*tau))); 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){ ...@@ -78,22 +73,17 @@ float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH){
if ((t<tshift) || (t>(tshift+ts))) amp=0.0; if ((t<tshift) || (t>(tshift+ts))) amp=0.0;
else amp=((sin(2.0*PI*(t-tshift)*fc) else amp=((sin(2.0*PI*(t-tshift)*fc)
-0.5*sin(4.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; break;
case 3 : case 3 :
/* source wavelet from file SOURCE_FILE */
if (nt<=nts) amp=psource[nt]; if (nt<=nts) amp=psource[nt];
else amp=0.0; else amp=0.0;
// amp=psource[nt]; break;
break; /* source wavelet from file SOURCE_FILE */
case 4 : case 4 :
/*tau=PI*(t-ts-tshift)/(1.5*ts);*/ /* Ricker */ /* sinus raised to the power of three */
/*amp=((t-ts-tshift)*exp(-2.0*tau*tau));*/
if ((t<tshift) || (t>(tshift+ts))) amp=0.0; if ((t<tshift) || (t>(tshift+ts))) amp=0.0;
else amp=pow(sin(PI*(t+tshift)/ts),3.0); else amp=pow(sin(PI*(t+tshift)/ts),3.0);
break; /* sinus raised to the power of three */ break;
break; break;
case 5 : case 5 :
...@@ -147,10 +137,15 @@ float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH){ ...@@ -147,10 +137,15 @@ float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH){
/* ---------------------------- */ /* ---------------------------- */
if (SH==1) { 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){ 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);
}
signals=fmatrix(1,nsrc,1,NT); signals=fmatrix(1,nsrc,1,NT);
...@@ -162,14 +157,9 @@ float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH){ ...@@ -162,14 +157,9 @@ float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH){
fc=srcpos_loc[5][k]; fc=srcpos_loc[5][k];
a=srcpos_loc[6][k]; a=srcpos_loc[6][k];
ts=1.0/fc; ts=1.0/fc;
ag = PI*PI*fc*fc;
switch (SOURCE_SHAPE_SH){ switch (SOURCE_SHAPE_SH){
case 1 : 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 */ /* New Ricker Wavelet, equal to SOFI2D */
tau=PI*(t-1.5*ts-tshift)/(ts); tau=PI*(t-1.5*ts-tshift)/(ts);
amp=(((1.0-2.0*tau*tau)*exp(-tau*tau))); 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){ ...@@ -178,22 +168,17 @@ float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH){
if ((t<tshift) || (t>(tshift+ts))) amp=0.0; if ((t<tshift) || (t>(tshift+ts))) amp=0.0;
else amp=((sin(2.0*PI*(t-tshift)*fc) else amp=((sin(2.0*PI*(t-tshift)*fc)
-0.5*sin(4.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; break;
case 3 : case 3 :
/* source wavelet from file SOURCE_FILE */
if (nt<=nts) amp=psource[nt]; if (nt<=nts) amp=psource[nt];
else amp=0.0; else amp=0.0;
// amp=psource[nt]; break;
break; /* source wavelet from file SOURCE_FILE */
case 4 : case 4 :
/*tau=PI*(t-ts-tshift)/(1.5*ts);*/ /* Ricker */ /* sinus raised to the power of three */
/*amp=((t-ts-tshift)*exp(-2.0*tau*tau));*/
if ((t<tshift) || (t>(tshift+ts))) amp=0.0; if ((t<tshift) || (t>(tshift+ts))) amp=0.0;
else amp=pow(sin(PI*(t+tshift)/ts),3.0); else amp=pow(sin(PI*(t+tshift)/ts),3.0);
break; /* sinus raised to the power of three */ break;
break; break;
case 5 : case 5 :
......
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