Commit dd778c01 authored by Florian Wittkamp's avatar Florian Wittkamp

Switched naming of frequency filter variables

I switched the naming of the variables of the frequency filtering to a more clear and easy one.
The change affects the JSON input file and the code itself.
The following renaming is used:
FC_START -> F_LOW_PASS_START
FC_END   -> F_LOW_PASS_END
FC_INCR  -> F_LOW_PASS_INCR
FC_EXT   -> F_LOW_PASS_EXT
FC       -> F_LOW_PASS
F_HP     -> F_HIGH_PASS

However, old JSON files will also work for now. If one of the old variables is not found, IFOS will search for the corresponding old one. This option will be removed in a future release.
parent 879a9627
......@@ -686,10 +686,10 @@ For more information see chapter \ref{cha:STF-Inversion}.
{\color{blue}{\begin{verbatim}
"Frequency filtering during inversion" : "comment",
"TIME_FILT" : "0",
"F_HP" : "1",
"FC_START" : "10.0",
"FC_END" : "75.0",
"FC_INCR" : "10.0",
"F_HIGH_PASS" : "1",
"F_LOW_PASS_START" : "10.0",
"F_LOW_PASS_END" : "75.0",
"F_LOW_PASS_INCR" : "10.0",
"ORDER" : "2",
"ZERO_PHASE" : "0",
"FREQ_FILE" : "frequencies.dat",
......@@ -706,13 +706,13 @@ Default values are:
MIN_ITER=0
\end{verbatim}}}
TIME\_FILT = 1 can be set to use frequency filtering. The parameter FC\_START defines the corner frequency of the Butterworth low pass filter at the beginning of the inversion. The parameter FC\_END defines the maximum corner frequency used in the inversion. The parameter FC\_INCR controls in which steps the bandwidth is increased during the inversion.
TIME\_FILT = 1 can be set to use frequency filtering. The parameter F\_LOW\_PASS\_START defines the corner frequency of the Butterworth low pass filter at the beginning of the inversion. The parameter F\_LOW\_PASS\_END defines the maximum corner frequency used in the inversion. The parameter F\_LOW\_PASS\_INCR controls in which steps the bandwidth is increased during the inversion.
If TIME\_FILT = 2 individual frequencies for each step can be read from FREQ\_FILE. In this file the first entry must be the number of frequencies used for filtering. Each frequency in Hz has to be specified in a row. The example file frequencies.dat can be found in \texttt{trunk/par}.
The parameter ORDER defines the order of the Butterworth low pass filter. If the variable ZERO\_PHASE is set to one a zero phase filter is applied. It is realized by filtering the traces in both forward and reverse direction with the defined Butterworth low pass filter. Therefore, the effective order of the low pass filter is doubled.
With F\_HP an additional high pass filter can be applied, where F\_HP is the corner frequency in Hz.
With F\_HIGH\_PASS an additional high pass filter can be applied, where F\_HIGH\_PASS is the corner frequency in Hz.
With the parameter PRO (see~\ref{json:abort_criterion}) one has to adjust the criterion that defines at which points the bandwidth of the signals are increased.
......
......@@ -164,10 +164,10 @@
"Frequency filtering during inversion" : "comment",
"TIME_FILT" : "0",
"F_HP" : "1",
"FC_START" : "10.0",
"FC_END" : "75.0",
"FC_INCR" : "10.0",
"F_HIGH_PASS" : "1",
"F_LOW_PASS_START" : "10.0",
"F_LOW_PASS_END" : "75.0",
"F_LOW_PASS_INCR" : "10.0",
"ORDER" : "2",
"ZERO_PHASE" : "0",
"FREQ_FILE" : "frequencies.dat",
......
......@@ -127,8 +127,8 @@
"Frequency filtering during inversion" : "comment",
"TIME_FILT" : "1",
"FC_START" : "10.0",
"FC_END" : "70.0",
"FC_INCR" : "10.0",
"F_LOW_PASS_END" : "70.0",
"F_LOW_PASS_INCR" : "10.0",
"ORDER" : "4",
"Minimum number of iteration per frequency" : "comment",
......
......@@ -139,9 +139,9 @@
"Frequency filtering during inversion" : "comment",
"TIME_FILT" : "0",
"FC_START" : "5.0",
"FC_END" : "5.0",
"FC_INCR" : "10.0",
"F_LOW_PASS_START" : "5.0",
"F_LOW_PASS_END" : "5.0",
"F_LOW_PASS_INCR" : "10.0",
"ORDER" : "4",
"Minimum number of iteration per frequency" : "comment",
......
......@@ -136,10 +136,10 @@
"Frequency filtering during inversion" : "comment",
"TIME_FILT" : "0",
"F_HP" : "1",
"FC_START" : "10.0",
"FC_END" : "70.0",
"FC_INCR" : "10.0",
"F_HIGH_PASS" : "1",
"F_LOW_PASS_START" : "10.0",
"F_LOW_PASS_END" : "70.0",
"F_LOW_PASS_INCR" : "10.0",
"ORDER" : "4",
"ZERO_PHASE" : "1",
......
......@@ -149,8 +149,8 @@ int main(int argc, char **argv){
WORKFLOW_STAGE=1;
/* variable for time domain filtering */
float FC;
float *FC_EXT=NULL;
float F_LOW_PASS;
float *F_LOW_PASS_EXT=NULL;
int nfrq=0;
int FREQ_NR=1;
......@@ -907,9 +907,9 @@ int main(int argc, char **argv){
RECINC=1;
switch(TIME_FILT){
case 1: FC=FC_START; break;
case 1: F_LOW_PASS=F_LOW_PASS_START; break;
/*read frequencies from file*/
case 2: FC_EXT=filter_frequencies(&nfrq); FC=FC_EXT[FREQ_NR]; break;
case 2: F_LOW_PASS_EXT=filter_frequencies(&nfrq); F_LOW_PASS=F_LOW_PASS_EXT[FREQ_NR]; break;
}
/* Save old SOURCE_SHAPE, which is needed for STF */
......@@ -927,7 +927,7 @@ int main(int argc, char **argv){
// At each iteration the workflow is applied
if(USE_WORKFLOW&&(FORWARD_ONLY==0)){
apply_workflow(workflow,workflow_lines,workflow_header,&iter,&FC,wavetype_start,&change_wavetype_iter,&LBFGS_iter_start);
apply_workflow(workflow,workflow_lines,workflow_header,&iter,&F_LOW_PASS,wavetype_start,&change_wavetype_iter,&LBFGS_iter_start);
}
......@@ -1084,7 +1084,7 @@ int main(int argc, char **argv){
if (TIME_FILT==0){
fprintf(FPL2,"opteps_vp \t epst1[1] \t epst1[2] \t epst1[3] \t L2t[1] \t L2t[2] \t L2t[3] \t L2t[4] \n");}
else{
fprintf(FPL2,"opteps_vp \t epst1[1] \t epst1[2] \t epst1[3] \t L2t[1] \t L2t[2] \t L2t[3] \t L2t[4] \t FC \n");
fprintf(FPL2,"opteps_vp \t epst1[1] \t epst1[2] \t epst1[3] \t L2t[1] \t L2t[2] \t L2t[3] \t L2t[4] \t F_LOW_PASS \n");
}
}
}
......@@ -1447,21 +1447,21 @@ int main(int argc, char **argv){
if(WAVETYPE==1 || WAVETYPE==3){
if ((ADJOINT_TYPE==1)|| (ADJOINT_TYPE==2)){
inseis(fprec,ishot,sectionvy_obs,ntr_glob,ns,2,iter);
timedomain_filt(sectionvy_obs,FC,ORDER,ntr_glob,ns,1);
timedomain_filt(sectionvy_obs,F_LOW_PASS,ORDER,ntr_glob,ns,1);
}
if (ADJOINT_TYPE==4){
inseis(fprec,ishot,sectionp_obs,ntr_glob,ns,9,iter);
timedomain_filt(sectionp_obs,FC,ORDER,ntr_glob,ns,1);
timedomain_filt(sectionp_obs,F_LOW_PASS,ORDER,ntr_glob,ns,1);
}
}
if(WAVETYPE==2 || WAVETYPE==3){
inseis(fprec,ishot,sectionvz_obs,ntr_glob,ns,10,iter);
timedomain_filt(sectionvz_obs,FC,ORDER,ntr_glob,ns,1);
timedomain_filt(sectionvz_obs,F_LOW_PASS,ORDER,ntr_glob,ns,1);
}
printf("\n ====================================================================================================== \n");
printf("\n Time Domain Filter is used for the inversion: lowpass filter, corner frequency of %.2f Hz, order %d\n",FC,ORDER);
printf("\n Time Domain Filter is used for the inversion: lowpass filter, corner frequency of %.2f Hz, order %d\n",F_LOW_PASS,ORDER);
printf("\n ====================================================================================================== \n");
if(iter==1){
......@@ -1475,15 +1475,15 @@ int main(int argc, char **argv){
if(WAVETYPE==1 || WAVETYPE==3){
if ((ADJOINT_TYPE==1)|| (ADJOINT_TYPE==2)){
stf(FP,fulldata_vy,sectionvy_obs,sectionvy_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,FC,0,nsrc_glob);
stf(FP,fulldata_vy,sectionvy_obs,sectionvy_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,F_LOW_PASS,0,nsrc_glob);
}
if (ADJOINT_TYPE==4){
stf(FP,fulldata_p,sectionp_obs,sectionp_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,FC,0,nsrc_glob);
stf(FP,fulldata_p,sectionp_obs,sectionp_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,F_LOW_PASS,0,nsrc_glob);
}
}
if(WAVETYPE==2 || WAVETYPE==3){
stf(FP,fulldata_vz,sectionvz_obs,sectionvz_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,FC,1,nsrc_glob);
stf(FP,fulldata_vz,sectionvz_obs,sectionvz_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,F_LOW_PASS,1,nsrc_glob);
}
......@@ -1513,15 +1513,15 @@ int main(int argc, char **argv){
}
if ((ADJOINT_TYPE==1)|| (ADJOINT_TYPE==2)){
stf(FP,fulldata_vy,sectionvy_obs,sectionvy_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,FC,0,nsrc_glob);
stf(FP,fulldata_vy,sectionvy_obs,sectionvy_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,F_LOW_PASS,0,nsrc_glob);
}
if (ADJOINT_TYPE==4){
stf(FP,fulldata_p,sectionp_obs,sectionp_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,FC,0,nsrc_glob);
stf(FP,fulldata_p,sectionp_obs,sectionp_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,F_LOW_PASS,0,nsrc_glob);
}
}
if(WAVETYPE==2 || WAVETYPE==3){
inseis(fprec,ishot,sectionvz_obs,ntr_glob,ns,10,iter);
stf(FP,fulldata_vz,sectionvz_obs,sectionvz_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,FC,1,nsrc_glob);
stf(FP,fulldata_vz,sectionvz_obs,sectionvz_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,F_LOW_PASS,1,nsrc_glob);
}
}
}
......@@ -1592,11 +1592,11 @@ int main(int argc, char **argv){
/*------------------------------------------------------------------------------*/
if (((TIME_FILT==1) || (TIME_FILT==2)) && (SOURCE_SHAPE!=6) && (INV_STF==0)){
fprintf(FP,"\n Time Domain Filter applied: Lowpass with corner frequency of %.2f Hz, order %d\n",FC,ORDER);
fprintf(FP,"\n Time Domain Filter applied: Lowpass with corner frequency of %.2f Hz, order %d\n",F_LOW_PASS,ORDER);
/*time domain filtering of the source signal */
if(WAVETYPE==1||WAVETYPE==3) timedomain_filt(signals,FC,ORDER,nsrc_loc,ns,1);
if(WAVETYPE==2||WAVETYPE==3) timedomain_filt(signals_SH,FC,ORDER,nsrc_loc,ns,1);
if(WAVETYPE==1||WAVETYPE==3) timedomain_filt(signals,F_LOW_PASS,ORDER,nsrc_loc,ns,1);
if(WAVETYPE==2||WAVETYPE==3) timedomain_filt(signals_SH,F_LOW_PASS,ORDER,nsrc_loc,ns,1);
}
/*------------------------------------------------------------------------------*/
......@@ -1989,7 +1989,7 @@ int main(int argc, char **argv){
if((ADJOINT_TYPE==1)||(ADJOINT_TYPE==3)){ /* if ADJOINT_TYPE */
inseis(fprec,ishot,sectionread,ntr_glob,ns,1,iter);
if ((TIME_FILT==1 )|| (TIME_FILT==2)){
timedomain_filt(sectionread,FC,ORDER,ntr_glob,ns,1);
timedomain_filt(sectionread,F_LOW_PASS,ORDER,ntr_glob,ns,1);
}
h=1;
for(i=1;i<=ntr;i++){
......@@ -2012,7 +2012,7 @@ int main(int argc, char **argv){
if((ADJOINT_TYPE==1)||(ADJOINT_TYPE==2)){ /* if ADJOINT_TYPE */
inseis(fprec,ishot,sectionread,ntr_glob,ns,2,iter);
if ((TIME_FILT==1 )|| (TIME_FILT==2)){
timedomain_filt(sectionread,FC,ORDER,ntr_glob,ns,1);
timedomain_filt(sectionread,F_LOW_PASS,ORDER,ntr_glob,ns,1);
}
h=1;
for(i=1;i<=ntr;i++){
......@@ -2034,7 +2034,7 @@ int main(int argc, char **argv){
if(ADJOINT_TYPE==4){ /* if ADJOINT_TYPE */
inseis(fprec,ishot,sectionread,ntr_glob,ns,9,iter);
if ((TIME_FILT==1 )|| (TIME_FILT==2)){
timedomain_filt(sectionread,FC,ORDER,ntr_glob,ns,1);
timedomain_filt(sectionread,F_LOW_PASS,ORDER,ntr_glob,ns,1);
}
h=1;
for(i=1;i<=ntr;i++){
......@@ -2056,7 +2056,7 @@ int main(int argc, char **argv){
if(WAVETYPE==2 || WAVETYPE==3){
inseis(fprec,ishot,sectionread,ntr_glob,ns,10,iter);
if ((TIME_FILT==1 )|| (TIME_FILT==2)){
timedomain_filt(sectionread,FC,ORDER,ntr_glob,ns,1);
timedomain_filt(sectionread,F_LOW_PASS,ORDER,ntr_glob,ns,1);
}
h=1;
for(i=1;i<=ntr;i++){
......@@ -2902,8 +2902,8 @@ int main(int argc, char **argv){
/* ----------------------------------------------------------------------*/
if((GRAD_METHOD==1)){
if (WAVETYPE==1 || WAVETYPE==3) PCG(waveconv, taper_coeff, nsrc, srcpos, recpos, ntr_glob, iter, C_vp, gradp, nfstart_jac, waveconv_u, C_vs, gradp_u, waveconv_rho, C_rho, gradp_rho,Vs_avg,FC);
if (WAVETYPE==2 || WAVETYPE==3) PCG_SH(taper_coeff, nsrc, srcpos, recpos, ntr_glob, iter, nfstart_jac, waveconv_u_z, C_vs, gradp_u_z, waveconv_rho_z, C_rho, gradp_rho_z,Vs_avg,FC);
if (WAVETYPE==1 || WAVETYPE==3) PCG(waveconv, taper_coeff, nsrc, srcpos, recpos, ntr_glob, iter, C_vp, gradp, nfstart_jac, waveconv_u, C_vs, gradp_u, waveconv_rho, C_rho, gradp_rho,Vs_avg,F_LOW_PASS);
if (WAVETYPE==2 || WAVETYPE==3) PCG_SH(taper_coeff, nsrc, srcpos, recpos, ntr_glob, iter, nfstart_jac, waveconv_u_z, C_vs, gradp_u_z, waveconv_rho_z, C_rho, gradp_rho_z,Vs_avg,F_LOW_PASS);
}
......@@ -2959,7 +2959,7 @@ int main(int argc, char **argv){
if(WAVETYPE==1 || WAVETYPE==3){
if (SWS_TAPER_FILE){ /* read taper from BIN-File*/
taper_grad(waveconv,taper_coeff,srcpos,nsrc,recpos,ntr_glob,4);}
if(GRAD_FILTER==1){smooth(waveconv,1,1,Vs_avg,FC);}
if(GRAD_FILTER==1){smooth(waveconv,1,1,Vs_avg,F_LOW_PASS);}
}
if (SWS_TAPER_FILE && !ACOUSTIC){ /* read taper from BIN-File*/
......@@ -2968,8 +2968,8 @@ int main(int argc, char **argv){
if (SWS_TAPER_FILE){ /* read taper from BIN-File*/
taper_grad(waveconv_rho,taper_coeff,srcpos,nsrc,recpos,ntr_glob,6);}
if(GRAD_FILTER==1&& !ACOUSTIC){smooth(waveconv_u,2,1,Vs_avg,FC);}
if(GRAD_FILTER==1){smooth(waveconv_rho,3,1,Vs_avg,FC);}
if(GRAD_FILTER==1&& !ACOUSTIC){smooth(waveconv_u,2,1,Vs_avg,F_LOW_PASS);}
if(GRAD_FILTER==1){smooth(waveconv_rho,3,1,Vs_avg,F_LOW_PASS);}
if(WOLFE_CONDITION) {
for (j=1;j<=NY;j=j+IDY){
......@@ -3093,7 +3093,7 @@ int main(int argc, char **argv){
if (TIME_FILT==0){
if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \n",0.0,iter,wolfe_sum_FWI,0.0,countstep-1,L2_SL_old,L2_SL_old);}
else{
if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \t %f\n",0.0,iter,wolfe_sum_FWI,0.0,countstep-1,L2_SL_old,L2_SL_old,FC);
if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \t %f\n",0.0,iter,wolfe_sum_FWI,0.0,countstep-1,L2_SL_old,L2_SL_old,F_LOW_PASS);
}
/* No update is done here, however model fils are written to disk for easy post processing */
......@@ -3134,7 +3134,7 @@ int main(int argc, char **argv){
if (TIME_FILT==0){
if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \n",alpha_SL,iter,wolfe_sum_FWI,diff,countstep-1,L2_SL_old,L2_SL_new);}
else{
if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \t %f\n",alpha_SL,iter,wolfe_sum_FWI,diff,countstep-1,L2_SL_old,L2_SL_new,FC);
if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \t %f\n",alpha_SL,iter,wolfe_sum_FWI,diff,countstep-1,L2_SL_old,L2_SL_new,F_LOW_PASS);
}
/* initiate variables for next iteration */
......@@ -3289,7 +3289,7 @@ int main(int argc, char **argv){
/*time domain filtering of the source signal */
if(WAVETYPE==1||WAVETYPE==3){
if(((TIME_FILT==1) || (TIME_FILT==2)) && (INV_STF==0)){
timedomain_filt(signals,FC,ORDER,nsrc_loc,ns,1);
timedomain_filt(signals,F_LOW_PASS,ORDER,nsrc_loc,ns,1);
}
}
......@@ -3297,7 +3297,7 @@ int main(int argc, char **argv){
/*time domain filtering of the source signal */
if(WAVETYPE==2||WAVETYPE==3){
if(((TIME_FILT==1) || (TIME_FILT==2)) && (INV_STF==0)){
timedomain_filt(signals_SH,FC,ORDER,nsrc_loc,ns,1);
timedomain_filt(signals_SH,F_LOW_PASS,ORDER,nsrc_loc,ns,1);
}
}
......@@ -3486,7 +3486,7 @@ int main(int argc, char **argv){
if((ADJOINT_TYPE==1)||(ADJOINT_TYPE==3)){ /* if ADJOINT_TYPE */
inseis(fprec,ishot,sectionread,ntr_glob,ns,1,iter);
if ((TIME_FILT==1 )|| (TIME_FILT==2)){
timedomain_filt(sectionread,FC,ORDER,ntr_glob,ns,1);
timedomain_filt(sectionread,F_LOW_PASS,ORDER,ntr_glob,ns,1);
}
h=1;
for(i=1;i<=ntr;i++){
......@@ -3505,7 +3505,7 @@ int main(int argc, char **argv){
if((ADJOINT_TYPE==1)||(ADJOINT_TYPE==2)){ /* if ADJOINT_TYPE */
inseis(fprec,ishot,sectionread,ntr_glob,ns,2,iter);
if ((TIME_FILT==1 )|| (TIME_FILT==2)){
timedomain_filt(sectionread,FC,ORDER,ntr_glob,ns,1);
timedomain_filt(sectionread,F_LOW_PASS,ORDER,ntr_glob,ns,1);
}
h=1;
for(i=1;i<=ntr;i++){
......@@ -3523,7 +3523,7 @@ int main(int argc, char **argv){
if(ADJOINT_TYPE==4){ /* if ADJOINT_TYPE */
inseis(fprec,ishot,sectionread,ntr_glob,ns,9,iter);
if ((TIME_FILT==1 )|| (TIME_FILT==2)){
timedomain_filt(sectionread,FC,ORDER,ntr_glob,ns,1);
timedomain_filt(sectionread,F_LOW_PASS,ORDER,ntr_glob,ns,1);
}
h=1;
for(i=1;i<=ntr;i++){
......@@ -3542,7 +3542,7 @@ int main(int argc, char **argv){
if(WAVETYPE==2 || WAVETYPE==3){
inseis(fprec,ishot,sectionread,ntr_glob,ns,10,iter);
if ((TIME_FILT==1 )|| (TIME_FILT==2)){
timedomain_filt(sectionread,FC,ORDER,ntr_glob,ns,1);
timedomain_filt(sectionread,F_LOW_PASS,ORDER,ntr_glob,ns,1);
}
h=1;
for(i=1;i<=ntr;i++){
......@@ -3774,7 +3774,7 @@ int main(int argc, char **argv){
if (TIME_FILT==0){
fprintf(FPL2,"%e \t %e \t %e \t %e \t %e \t %e \t %e \t %e \n",opteps_vp,epst1[1],epst1[2],epst1[3],L2t[1],L2t[2],L2t[3],L2t[4]);}
else{
fprintf(FPL2,"%e \t %e \t %e \t %e \t %e \t %e \t %e \t %e \t %f\n",opteps_vp,epst1[1],epst1[2],epst1[3],L2t[1],L2t[2],L2t[3],L2t[4],FC);}
fprintf(FPL2,"%e \t %e \t %e \t %e \t %e \t %e \t %e \t %e \t %f\n",opteps_vp,epst1[1],epst1[2],epst1[3],L2t[1],L2t[2],L2t[3],L2t[4],F_LOW_PASS);}
}
/* saving history of final L2*/
......@@ -3795,12 +3795,12 @@ int main(int argc, char **argv){
/* ------------------------------------*/
if (FORWARD_ONLY==0 && (opteps_vp>0.0 || WOLFE_CONDITION)){
if(!ACOUSTIC){
if(WAVETYPE==1||WAVETYPE==3) if(MODEL_FILTER)smooth(ppi,4,2,Vs_avg,FC);
if(MODEL_FILTER)smooth(pu,5,2,Vs_avg,FC);
if(MODEL_FILTER)smooth(prho,6,2,Vs_avg,FC);
if(WAVETYPE==1||WAVETYPE==3) if(MODEL_FILTER)smooth(ppi,4,2,Vs_avg,F_LOW_PASS);
if(MODEL_FILTER)smooth(pu,5,2,Vs_avg,F_LOW_PASS);
if(MODEL_FILTER)smooth(prho,6,2,Vs_avg,F_LOW_PASS);
}else{
if(WAVETYPE==1||WAVETYPE==3) if(MODEL_FILTER)smooth(ppi,4,2,Vp_avg,FC);
if(MODEL_FILTER)smooth(prho,6,2,Vp_avg,FC);
if(WAVETYPE==1||WAVETYPE==3) if(MODEL_FILTER)smooth(ppi,4,2,Vp_avg,F_LOW_PASS);
if(MODEL_FILTER)smooth(prho,6,2,Vp_avg,F_LOW_PASS);
}
}
......@@ -3917,19 +3917,19 @@ int main(int argc, char **argv){
printf("\n Did not find a step length which decreases the misfit.\n");}
}
FC=FC+FC_INCR;
F_LOW_PASS=F_LOW_PASS+F_LOW_PASS_INCR;
do_stf=1;
min_iter_help=0;
min_iter_help=iter+MIN_ITER;
if(FC>FC_END){
if(F_LOW_PASS>F_LOW_PASS_END){
if(MYID==0){
printf("\n Reached the maximum frequency of %4.2f Hz \n",FC);
printf("\n Reached the maximum frequency of %4.2f Hz \n",F_LOW_PASS);
}
break;
}
if(MYID==0) printf("\n Changing to corner frequency of %4.2f Hz \n",FC);
if(MYID==0) printf("\n Changing to corner frequency of %4.2f Hz \n",F_LOW_PASS);
/* Restart L-BFGS at next iteration */
LBFGS_iter_start=iter+1;
......@@ -3950,17 +3950,17 @@ int main(int argc, char **argv){
if(FREQ_NR==nfrq){
if(MYID==0){
printf("\n Finished at the maximum frequency of %4.2f Hz \n",FC);
printf("\n Finished at the maximum frequency of %4.2f Hz \n",F_LOW_PASS);
}
break;
}
FREQ_NR=FREQ_NR+1;
FC=FC_EXT[FREQ_NR];
F_LOW_PASS=F_LOW_PASS_EXT[FREQ_NR];
do_stf=1;
min_iter_help=0;
min_iter_help=iter+MIN_ITER;
if(MYID==0) printf("\n Changing to corner frequency of %4.2f Hz \n",FC);
if(MYID==0) printf("\n Changing to corner frequency of %4.2f Hz \n",F_LOW_PASS);
/* Restart L-BFGS at next iteration */
LBFGS_iter_start=iter+1;
......@@ -4348,7 +4348,7 @@ int main(int argc, char **argv){
if(TIME_FILT==2){
free_vector(FC_EXT,1,nfrq);
free_vector(F_LOW_PASS_EXT,1,nfrq);
}
if (nsrc_loc>0){
......
......@@ -24,7 +24,7 @@
#include "fd.h"
void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, float C_vp, float ** gradp, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float FC){
void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, float C_vp, float ** gradp, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS){
extern int NX, NY, IDX, IDY, SPATFILTER, GRAD_FILTER;
extern int FORWARD_ONLY, SWS_TAPER_GRAD_VERT, SWS_TAPER_GRAD_HOR, SWS_TAPER_GRAD_SOURCES, SWS_TAPER_FILE;
......@@ -99,7 +99,7 @@ void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int
spat_filt(waveconv,iter,1);}
/* apply 2D-Gaussian filter*/
if(GRAD_FILTER==1){smooth(waveconv,1,1,Vs_avg,FC);}
if(GRAD_FILTER==1){smooth(waveconv,1,1,Vs_avg,F_LOW_PASS);}
/* output of the preconditioned gradient */
for (i=1;i<=NX;i=i+IDX){
......@@ -331,7 +331,7 @@ void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int
spat_filt(waveconv_u,iter,2);}
/* apply 2D-Gaussian filter*/
if(GRAD_FILTER==1){smooth(waveconv_u,2,1,Vs_avg,FC);}
if(GRAD_FILTER==1){smooth(waveconv_u,2,1,Vs_avg,F_LOW_PASS);}
/* output of the preconditioned gradient */
for (i=1;i<=NX;i=i+IDX){
......@@ -558,7 +558,7 @@ void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int
spat_filt(waveconv_rho,iter,3);}
/* apply 2D-Gaussian filter*/
if(GRAD_FILTER==1){smooth(waveconv_rho,3,1,Vs_avg,FC);}
if(GRAD_FILTER==1){smooth(waveconv_rho,3,1,Vs_avg,F_LOW_PASS);}
/* output of the preconditioned gradient */
for (i=1;i<=NX;i=i+IDX){
......
......@@ -24,7 +24,7 @@
#include "fd.h"
void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float FC){
void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS){
extern int NX, NY, IDX, IDY, SPATFILTER, GRAD_FILTER;
extern int FORWARD_ONLY, SWS_TAPER_GRAD_VERT, SWS_TAPER_GRAD_HOR, SWS_TAPER_GRAD_SOURCES, SWS_TAPER_FILE;
......@@ -91,7 +91,7 @@ void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int
spat_filt(waveconv_u,iter,2);}
/* apply 2D-Gaussian filter*/
if(GRAD_FILTER==1){smooth(waveconv_u,2,1,Vs_avg,FC);}
if(GRAD_FILTER==1){smooth(waveconv_u,2,1,Vs_avg,F_LOW_PASS);}
/* output of the preconditioned gradient */
for (i=1;i<=NX;i=i+IDX){
......@@ -320,7 +320,7 @@ void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int
spat_filt(waveconv_rho,iter,3);}
/* apply 2D-Gaussian filter*/
if(GRAD_FILTER==1){smooth(waveconv_rho,3,1,Vs_avg,FC);}
if(GRAD_FILTER==1){smooth(waveconv_rho,3,1,Vs_avg,F_LOW_PASS);}
/* output of the preconditioned gradient */
for (i=1;i<=NX;i=i+IDX){
......
......@@ -23,7 +23,7 @@
#include "fd.h"
void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[STRING_SIZE],int *iter,float *FC,int wavetype_start, int * change_wavetype_iter, int * LBFGS_iter_start){
void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[STRING_SIZE],int *iter,float *F_LOW_PASS,int wavetype_start, int * change_wavetype_iter, int * LBFGS_iter_start){
/* local variables */
int x;
......@@ -90,12 +90,12 @@ void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[ST
/* Frequency filtering */
if(TIME_FILT==1) {
TIME_FILT=workflow[WORKFLOW_STAGE][6];
if(*FC>workflow[WORKFLOW_STAGE][7]&&(workflow[WORKFLOW_STAGE][6]>0)) {
if(MYID==0)printf("\n Due to the abort criteriom FC is already higher than specified in workflow\n");
if(MYID==0)printf(" therefore instead of %.2f HZ FC=%.2f HZ is used\n",workflow[WORKFLOW_STAGE][7],*FC);
if(*F_LOW_PASS>workflow[WORKFLOW_STAGE][7]&&(workflow[WORKFLOW_STAGE][6]>0)) {
if(MYID==0)printf("\n Due to the abort criteriom F_LOW_PASS is already higher than specified in workflow\n");
if(MYID==0)printf(" therefore instead of %.2f HZ F_LOW_PASS=%.2f HZ is used\n",workflow[WORKFLOW_STAGE][7],*F_LOW_PASS);
} else {
if(*FC!=workflow[WORKFLOW_STAGE][7]) *LBFGS_iter_start=*iter;
*FC=workflow[WORKFLOW_STAGE][7];
if(*F_LOW_PASS!=workflow[WORKFLOW_STAGE][7]) *LBFGS_iter_start=*iter;
*F_LOW_PASS=workflow[WORKFLOW_STAGE][7];
}
} else {
if(MYID==0&&(workflow[WORKFLOW_STAGE][6]>0))printf("\n TIME_FILT cannot be activated due to it is not activated in the JSON File \n");
......
......@@ -26,7 +26,7 @@ double calc_energy(float **sectiondata, int ntr, int ns, float energy, int ntr_g
/* declaration of variables */
extern float DT;
extern int MYID, USE_WORKFLOW, WORKFLOW_STAGE;
extern int TRKILL, NORMALIZE, FC, TIMEWIN;
extern int TRKILL, NORMALIZE, F_LOW_PASS, TIMEWIN;
extern char TRKILL_FILE[STRING_SIZE];
extern int VELOCITY;
int i, j;
......
......@@ -26,7 +26,7 @@ double calc_misfit(float **sectiondata, float **section, int ntr, int ns, int LN
/* declaration of variables */
extern float DT;
extern int MYID, USE_WORKFLOW, WORKFLOW_STAGE;
extern int TRKILL, NORMALIZE, FC, TIMEWIN;
extern int TRKILL, NORMALIZE, F_LOW_PASS, TIMEWIN;
extern char TRKILL_FILE[STRING_SIZE];
extern int VELOCITY;
......
......@@ -26,7 +26,7 @@ double calc_res(float **sectiondata, float **section, float **sectiondiff, float
/* declaration of variables */
extern float DT, WATERLEVEL_LNORM8;
extern int REC1, REC2, MYID, ACOUSTIC;
extern int TRKILL, NORMALIZE, FC, TIMEWIN;
extern int TRKILL, NORMALIZE, F_LOW_PASS, TIMEWIN;
extern char TRKILL_FILE[STRING_SIZE];
extern int VELOCITY, USE_WORKFLOW, WORKFLOW_STAGE;
float RMS, signL1;
......
......@@ -60,7 +60,7 @@ void exchange_par(void){
extern int INV_STF, N_STF, N_STF_START;
extern char PARA[STRING_SIZE];
extern int TIME_FILT, ORDER, ZERO_PHASE,WRITE_FILTERED_DATA;
extern float FC_START, FC_END, FC_INCR, F_HP;
extern float F_LOW_PASS_START, F_LOW_PASS_END, F_LOW_PASS_INCR, F_HIGH_PASS;
extern int LNORM, DTINV;
extern int STEPMAX;
extern float EPS_SCALE, SCALEFAC;
......@@ -175,9 +175,9 @@ void exchange_par(void){
fdum[39] = npower;
fdum[40] = k_max_PML;
fdum[42] = FC_START;
fdum[43] = FC_END;
fdum[44] = FC_INCR;
fdum[42] = F_LOW_PASS_START;
fdum[43] = F_LOW_PASS_END;
fdum[44] = F_LOW_PASS_INCR;
fdum[45] = EPS_SCALE;
fdum[46] = SCALEFAC;
......@@ -203,7 +203,7 @@ void exchange_par(void){
fdum[58] = A;
fdum[59] = F_HP;
fdum[59] = F_HIGH_PASS;
fdum[60] = JOINT_INVERSION_PSV_SH_ALPHA_VS;
fdum[61] = JOINT_INVERSION_PSV_SH_ALPHA_RHO;
......@@ -457,9 +457,9 @@ void exchange_par(void){
k_max_PML = fdum[40];
FC_START = fdum[42];
FC_END = fdum[43];
FC_INCR = fdum[44];
F_LOW_PASS_START = fdum[42];
F_LOW_PASS_END = fdum[43];
F_LOW_PASS_INCR = fdum[44];
EPS_SCALE = fdum[45];
SCALEFAC = fdum[46];
......@@ -486,7 +486,7 @@ void exchange_par(void){
A = fdum[58];
F_HP = fdum[59];
F_HIGH_PASS = fdum[59];
JOINT_INVERSION_PSV_SH_ALPHA_VS = fdum[60];
JOINT_INVERSION_PSV_SH_ALPHA_RHO = fdum[61];
......
......@@ -33,7 +33,7 @@ void window_cos(float **win, int npad, int nsrc, float it1, float it2, float it3
void catseis(float **data, float **fulldata, int *recswitch, int ntr_glob, MPI_Comm newcomm_nodentr);
void stf(FILE *fp, float **sectionvy, float ** sectionvy_obs, float ** sectionvy_conv, float * source_time_function, int **recpos, int **recpos_loc,
int ntr_glob,int ntr, float ** srcpos, int ishot, int ns, int iter, int nshots, float FC, int SH,int nsrc_glob);
int ntr_glob,int ntr, float ** srcpos, int ishot, int ns, int iter, int nshots, float F_LOW_PASS, int SH,int nsrc_glob);
int **splitrec(int **recpos,int *ntr_loc, int ntr, int *recswitch);
......@@ -186,9 +186,9 @@ void taper(float *section, int ns, float fc);
void output_source_signal(FILE *fp, float **signals, int ns, int seis_form);
void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, float C_vp, float ** gradp, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float FC);
void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, float C_vp, float ** gradp, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS);
void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float FC);
void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS);
void PML_pro(float * d_x, float * K_x, float * alpha_prime_x, float * a_x, float * b_x,
float * d_x_half, float * K_x_half, float * alpha_prime_x_half, float * a_x_half, float * b_x_half,
......@@ -426,7 +426,7 @@ void zero_fdveps_visc(int ny1, int ny2, int nx1, int nx2, float ** vx, float **
void FLnode(float ** rho, float ** pi, float ** u, float ** taus, float ** taup, float * eta);
void smooth(float ** mat, int sws, int filter, float Vs_avg, float FC);
void smooth(float ** mat, int sws, int filter, float Vs_avg, float F_LOW_PASS);
/* declaration of functions for parser*/
......@@ -515,7 +515,7 @@ void read_workflow(char file_in[STRING_SIZE],float *** workflow, int *workflow_l
float ** joint_inversion_grad ( float ** gradiant_1,float ** gradiant_2, float alpha, int joint_type);
void snap_SH(FILE *fp,int nt, int nsnap, float ** vz, float **u, float **pi, float *hc,int ishot);
void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[STRING_SIZE],int *iter,float *FC,int wavetype_start, int * change_wavetype_iter, int * LBFGS_iter_start);
void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[STRING_SIZE],int *iter,float *F_LOW_PASS,int wavetype_start, int * change_wavetype_iter, int * LBFGS_iter_start);
void eprecond(float ** W, float ** vx, float ** vy);
void eprecond_SH(float ** W, float ** vz);
......
......@@ -31,7 +31,7 @@ extern FILE *FP;
FILE *freqf;
char cline[256];
int IT=1;
float *FC;
float *F_LOW_PASS;
/*----------------------------------- open FREQ_FILE ------------------------------*/
......@@ -48,8 +48,8 @@ fprintf(FP,"\n Reading frequencies from %s \n",FREQ_FILE);
fscanf(freqf,"%i",nfrq);
/*----------------------------alocate frequency array------------------------------*/
FC=vector(1,*nfrq);
//FC = malloc((*nfrq+1) * sizeof(float));
F_LOW_PASS=vector(1,*nfrq);
//F_LOW_PASS = malloc((*nfrq+1) * sizeof(float));
rewind (freqf);
fprintf(FP," Number of freqencies specified: %d \n",*nfrq);
......@@ -57,16 +57,16 @@ fprintf(FP," Number of freqencies specified: %d \n",*nfrq);
/*----------------------Read frequencies from FREQ_FILE----------------------------*/
for (IT=1;IT<=(*nfrq);IT++){
fgets(cline,255,freqf);
fscanf(freqf,"%f",&FC[IT]);
fscanf(freqf,"%f",&F_LOW_PASS[IT]);
//Read only numbers from FREQ_FILE
if(FC[IT]==0){
if(F_LOW_PASS[IT]==0){
IT=IT-1;}else{
fprintf(FP,"\n %d. %.2f Hz",IT,FC[IT]);
fprintf(FP,"\n %d. %.2f Hz",IT,F_LOW_PASS[IT]);
}
}
fclose(freqf);
return(FC);
return(F_LOW_PASS);
}
......@@ -74,7 +74,7 @@ char PARA[STRING_SIZE];
int TIME_FILT, ORDER, ZERO_PHASE;
int WRITE_FILTERED_DATA;
float FC_START, FC_END, FC_INCR, F_HP;
float F_LOW_PASS_START, F_LOW_PASS_END, F_LOW_PASS_INCR, F_HIGH_PASS;
int LNORM, DTINV;
......
......@@ -63,7 +63,7 @@ void read_par_json(FILE *fp, char *fileinp){
extern char PARA[STRING_SIZE];