Commit c7172e61 authored by Florian Wittkamp's avatar Florian Wittkamp

Automatic offset based tracekill

parent ba1116d6
...@@ -736,6 +736,10 @@ If you use the option "Workflow" (section \ref{sec:workflow}) it is possible to ...@@ -736,6 +736,10 @@ If you use the option "Workflow" (section \ref{sec:workflow}) it is possible to
"Trace killing" : "comment", "Trace killing" : "comment",
"TRKILL" : "0", "TRKILL" : "0",
"TRKILL_FILE" : "./trace_kill/trace_kill", "TRKILL_FILE" : "./trace_kill/trace_kill",
"TRKILL_OFFSET" : "0",
"TRKILL_OFFSET_LOWER" : "20",
"TRKILL_OFFSET_UPPER" : "100",
\end{verbatim}}} \end{verbatim}}}
{\color{red}{\begin{verbatim} {\color{red}{\begin{verbatim}
...@@ -747,6 +751,7 @@ If you don't want to use single traces or a subset of your data, the parameter T ...@@ -747,6 +751,7 @@ If you don't want to use single traces or a subset of your data, the parameter T
If you use the option "Workflow" (section \ref{sec:workflow}) it is possible to define separate files for each workflow stage. They have to be named TRKILL\_FILE\_<workflowstage>.dat. If you don't provide separate files for some or all stages the standard file (TRKILL\_FILE.dat) is used. This means a standard file should exist in any case. If you use the option "Workflow" (section \ref{sec:workflow}) it is possible to define separate files for each workflow stage. They have to be named TRKILL\_FILE\_<workflowstage>.dat. If you don't provide separate files for some or all stages the standard file (TRKILL\_FILE.dat) is used. This means a standard file should exist in any case.
Moreover, you can directly define a offset range that should be killed. Simply set TRKILL\_OFFSET and TRKILL to one and all traces with a offset between TRKILL\_OFFSET\_LOWER and TRKILL\_OFFSET\_UPPER will be killed.
\section{Gradient manipulation} \section{Gradient manipulation}
\subsection{Definition of a gradient taper} \subsection{Definition of a gradient taper}
......
...@@ -2032,10 +2032,10 @@ int main(int argc, char **argv){ ...@@ -2032,10 +2032,10 @@ int main(int argc, char **argv){
} }
h++; h++;
} }
L2=calc_res(sectionvxdata,sectionvx,sectionvxdiff,sectionvxdiffold,ntr,ns,LNORM,L2,0,1,swstestshot,ntr_glob,recpos_loc,nsrc_glob,ishot,iter); L2=calc_res(sectionvxdata,sectionvx,sectionvxdiff,sectionvxdiffold,ntr,ns,LNORM,L2,0,1,swstestshot,ntr_glob,recpos_loc,nsrc_glob,ishot,iter,srcpos,recpos);
if(swstestshot==1){energy=calc_energy(sectionvxdata,ntr,ns,energy, ntr_glob, recpos_loc, nsrc_glob, ishot,iter);} if(swstestshot==1){energy=calc_energy(sectionvxdata,ntr,ns,energy, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);}
L2_all_shots=calc_misfit(sectionvxdata,sectionvx,ntr,ns,LNORM,L2_all_shots,0,1,1, ntr_glob, recpos_loc, nsrc_glob, ishot,iter); L2_all_shots=calc_misfit(sectionvxdata,sectionvx,ntr,ns,LNORM,L2_all_shots,0,1,1, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
energy_all_shots=calc_energy(sectionvxdata,ntr,ns,energy_all_shots, ntr_glob, recpos_loc, nsrc_glob, ishot,iter); energy_all_shots=calc_energy(sectionvxdata,ntr,ns,energy_all_shots, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
/*fprintf(FP,"Energy vxdata for PE %d: %f\n\n", MYID,energy);*/ /*fprintf(FP,"Energy vxdata for PE %d: %f\n\n", MYID,energy);*/
} /* end ADJOINT_TYPE */ } /* end ADJOINT_TYPE */
...@@ -2055,10 +2055,10 @@ int main(int argc, char **argv){ ...@@ -2055,10 +2055,10 @@ int main(int argc, char **argv){
} }
h++; h++;
} }
L2=calc_res(sectionvydata,sectionvy,sectionvydiff,sectionvydiffold,ntr,ns,LNORM,L2,0,1,swstestshot,ntr_glob,recpos_loc,nsrc_glob,ishot,iter); L2=calc_res(sectionvydata,sectionvy,sectionvydiff,sectionvydiffold,ntr,ns,LNORM,L2,0,1,swstestshot,ntr_glob,recpos_loc,nsrc_glob,ishot,iter,srcpos,recpos);
if(swstestshot==1){energy=calc_energy(sectionvydata,ntr,ns,energy, ntr_glob, recpos_loc, nsrc_glob, ishot,iter);} if(swstestshot==1){energy=calc_energy(sectionvydata,ntr,ns,energy, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);}
L2_all_shots=calc_misfit(sectionvydata,sectionvy,ntr,ns,LNORM,L2_all_shots,0,1,1, ntr_glob, recpos_loc, nsrc_glob, ishot,iter); L2_all_shots=calc_misfit(sectionvydata,sectionvy,ntr,ns,LNORM,L2_all_shots,0,1,1, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
energy_all_shots=calc_energy(sectionvydata,ntr,ns,energy_all_shots, ntr_glob, recpos_loc, nsrc_glob, ishot,iter); energy_all_shots=calc_energy(sectionvydata,ntr,ns,energy_all_shots, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
/*fprintf(FP,"Energy vydata for PE %d: %f\n\n", MYID,energy); */ /*fprintf(FP,"Energy vydata for PE %d: %f\n\n", MYID,energy); */
} /* end ADJOINT_TYPE */ } /* end ADJOINT_TYPE */
...@@ -2077,10 +2077,10 @@ int main(int argc, char **argv){ ...@@ -2077,10 +2077,10 @@ int main(int argc, char **argv){
} }
h++; h++;
} }
L2=calc_res(sectionpdata,sectionp,sectionpdiff,sectionpdiffold,ntr,ns,LNORM,L2,0,1,swstestshot,ntr_glob,recpos_loc,nsrc_glob,ishot,iter); L2=calc_res(sectionpdata,sectionp,sectionpdiff,sectionpdiffold,ntr,ns,LNORM,L2,0,1,swstestshot,ntr_glob,recpos_loc,nsrc_glob,ishot,iter,srcpos,recpos);
if(swstestshot==1){energy=calc_energy(sectionpdata,ntr,ns,energy, ntr_glob, recpos_loc, nsrc_glob, ishot,iter);} if(swstestshot==1){energy=calc_energy(sectionpdata,ntr,ns,energy, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);}
L2_all_shots=calc_misfit(sectionpdata,sectionp,ntr,ns,LNORM,L2_all_shots,0,1,1, ntr_glob, recpos_loc, nsrc_glob, ishot,iter); L2_all_shots=calc_misfit(sectionpdata,sectionp,ntr,ns,LNORM,L2_all_shots,0,1,1, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
energy_all_shots=calc_energy(sectionpdata,ntr,ns,energy_all_shots, ntr_glob, recpos_loc, nsrc_glob, ishot,iter); energy_all_shots=calc_energy(sectionpdata,ntr,ns,energy_all_shots, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
} /* end ADJOINT_TYPE */ } /* end ADJOINT_TYPE */
} }
...@@ -2096,15 +2096,15 @@ int main(int argc, char **argv){ ...@@ -2096,15 +2096,15 @@ int main(int argc, char **argv){
} }
h++; h++;
} }
L2_SH=calc_res(sectionvzdata,sectionvz,sectionvzdiff,sectionvzdiffold,ntr,ns,LNORM,L2_SH,0,1,swstestshot,ntr_glob,recpos_loc,nsrc_glob,ishot,iter); L2_SH=calc_res(sectionvzdata,sectionvz,sectionvzdiff,sectionvzdiffold,ntr,ns,LNORM,L2_SH,0,1,swstestshot,ntr_glob,recpos_loc,nsrc_glob,ishot,iter,srcpos,recpos);
if(swstestshot==1){energy_SH=calc_energy(sectionvzdata,ntr,ns,energy_SH, ntr_glob, recpos_loc, nsrc_glob, ishot,iter);} if(swstestshot==1){energy_SH=calc_energy(sectionvzdata,ntr,ns,energy_SH, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);}
L2_all_shots_SH=calc_misfit(sectionvzdata,sectionvz,ntr,ns,LNORM,L2_all_shots_SH,0,1,1, ntr_glob, recpos_loc, nsrc_glob, ishot,iter); L2_all_shots_SH=calc_misfit(sectionvzdata,sectionvz,ntr,ns,LNORM,L2_all_shots_SH,0,1,1, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
energy_all_shots_SH=calc_energy(sectionvzdata,ntr,ns,energy_all_shots_SH, ntr_glob, recpos_loc, nsrc_glob, ishot,iter); energy_all_shots_SH=calc_energy(sectionvzdata,ntr,ns,energy_all_shots_SH, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
} }
// Tracekill // Tracekill
if (TRKILL){ if (TRKILL){
count_killed_traces(ntr,swstestshot,ntr_glob,recpos_loc,nsrc_glob,ishot,ptr_killed_traces,ptr_killed_traces_testshots); count_killed_traces(ntr,swstestshot,ntr_glob,recpos_loc,nsrc_glob,ishot,ptr_killed_traces,ptr_killed_traces_testshots,srcpos,recpos);
} }
if((ishot==itestshot)&&(ishot<=TESTSHOT_END)){ if((ishot==itestshot)&&(ishot<=TESTSHOT_END)){
...@@ -3523,7 +3523,7 @@ int main(int argc, char **argv){ ...@@ -3523,7 +3523,7 @@ int main(int argc, char **argv){
} }
h++; h++;
} }
L2=calc_res(sectionvxdata,sectionvx,sectionvxdiff,sectionvxdiffold,ntr,ns,LNORM,L2,0,1,1,ntr_glob,recpos_loc,nsrc_glob,ishot,iter); L2=calc_res(sectionvxdata,sectionvx,sectionvxdiff,sectionvxdiffold,ntr,ns,LNORM,L2,0,1,1,ntr_glob,recpos_loc,nsrc_glob,ishot,iter,srcpos,recpos);
} }
/* --------------------------------- */ /* --------------------------------- */
...@@ -3542,7 +3542,7 @@ int main(int argc, char **argv){ ...@@ -3542,7 +3542,7 @@ int main(int argc, char **argv){
} }
h++; h++;
} }
L2=calc_res(sectionvydata,sectionvy,sectionvydiff,sectionvydiffold,ntr,ns,LNORM,L2,0,1,1,ntr_glob,recpos_loc,nsrc_glob,ishot,iter); L2=calc_res(sectionvydata,sectionvy,sectionvydiff,sectionvydiffold,ntr,ns,LNORM,L2,0,1,1,ntr_glob,recpos_loc,nsrc_glob,ishot,iter,srcpos,recpos);
} }
/* --------------------------------- */ /* --------------------------------- */
...@@ -3560,7 +3560,7 @@ int main(int argc, char **argv){ ...@@ -3560,7 +3560,7 @@ int main(int argc, char **argv){
} }
h++; h++;
} }
L2=calc_res(sectionpdata,sectionp,sectionpdiff,sectionpdiffold,ntr,ns,LNORM,L2,0,1,1,ntr_glob,recpos_loc,nsrc_glob,ishot,iter); L2=calc_res(sectionpdata,sectionp,sectionpdiff,sectionpdiffold,ntr,ns,LNORM,L2,0,1,1,ntr_glob,recpos_loc,nsrc_glob,ishot,iter,srcpos,recpos);
} }
} }
...@@ -3576,7 +3576,7 @@ int main(int argc, char **argv){ ...@@ -3576,7 +3576,7 @@ int main(int argc, char **argv){
} }
h++; h++;
} }
L2_SH=calc_res(sectionvzdata,sectionvz,sectionvzdiff,sectionvzdiffold,ntr,ns,LNORM,L2_SH,0,1,1,ntr_glob,recpos_loc,nsrc_glob,ishot,iter); L2_SH=calc_res(sectionvzdata,sectionvz,sectionvzdiff,sectionvzdiffold,ntr,ns,LNORM,L2_SH,0,1,1,ntr_glob,recpos_loc,nsrc_glob,ishot,iter,srcpos,recpos);
} }
} }
......
...@@ -162,6 +162,7 @@ IFOS2D= \ ...@@ -162,6 +162,7 @@ IFOS2D= \
prepare_update_p.c \ prepare_update_p.c \
readmod_viscac.c \ readmod_viscac.c \
time_window_glob.c \ time_window_glob.c \
create_trkill_table.c \
filter_frequencies.c filter_frequencies.c
SNAPMERGE_OBJ = $(SNAPMERGE_SCR:%.c=%.o) SNAPMERGE_OBJ = $(SNAPMERGE_SCR:%.c=%.o)
......
...@@ -14,7 +14,7 @@ ...@@ -14,7 +14,7 @@
* *
* You should have received a copy of the GNU General Public License * You should have received a copy of the GNU General Public License
* along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>. * along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/ -----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------ /*------------------------------------------------------------------------
* Calculate energy of measured data * Calculate energy of measured data
...@@ -22,38 +22,56 @@ ...@@ -22,38 +22,56 @@
* ----------------------------------------------------------------------*/ * ----------------------------------------------------------------------*/
#include "fd.h" #include "fd.h"
double calc_energy(float **sectiondata, int ntr, int ns, float energy, int ntr_glob, int **recpos_loc, int nsrc_glob, int ishot, int iter){ double calc_energy(float **sectiondata, int ntr, int ns, float energy, int ntr_glob, int **recpos_loc, int nsrc_glob, int ishot, int iter, float ** srcpos, int ** recpos){
/* declaration of variables */ /* declaration of variables */
extern float DT; extern float DT;
extern int MYID, USE_WORKFLOW, WORKFLOW_STAGE; extern int MYID, USE_WORKFLOW, WORKFLOW_STAGE;
extern int TRKILL, NORMALIZE, FC, TIMEWIN; extern int TRKILL, NORMALIZE, FC, TIMEWIN;
extern char TRKILL_FILE[STRING_SIZE]; extern char TRKILL_FILE[STRING_SIZE];
extern int VELOCITY; extern int VELOCITY;
int i, j; int i, j;
float energy_dummy, intseis; float energy_dummy, intseis;
int umax=0, h; int umax=0, h;
float intseis_sd; float intseis_sd;
float **intseis_sectiondata=NULL; float **intseis_sectiondata=NULL;
float *picked_times=NULL; float *picked_times=NULL;
intseis_sectiondata = matrix(1,ntr,1,ns); /* declaration of variables for integration */ intseis_sectiondata = matrix(1,ntr,1,ns); /* declaration of variables for integration */
if(TIMEWIN) picked_times = vector(1,ntr); /* declaration of variables for TIMEWIN */ if(TIMEWIN) picked_times = vector(1,ntr); /* declaration of variables for TIMEWIN */
/*extern int MYID;*/
/*printf("MYID: %d; ns=%d\n",MYID,ns);
printf("MYID: %d; ntr=%d\n",MYID,ntr);*/
/* TRACE KILLING */
int ** kill_tmp, *kill_vector; /* declaration of variables for trace killing */
char trace_kill_file[STRING_SIZE];
FILE *ftracekill;
extern int TRKILL_OFFSET;
extern float TRKILL_OFFSET_LOWER;
extern float TRKILL_OFFSET_UPPER;
/*-----------------*/
/* Start Tracekill */
/*-----------------*/
if(TRKILL){
kill_tmp = imatrix(1,ntr_glob,1,nsrc_glob);
kill_vector = ivector(1,ntr);
/*extern int MYID;*/ if(TRKILL_OFFSET) {
/*printf("MYID: %d; ns=%d\n",MYID,ns); if(MYID==0) {
printf("MYID: %d; ntr=%d\n",MYID,ntr);*/ printf("Automatic offset based TraceKill: Experimental feature");
}
/* TRACE KILLING */ /* Generate TraceKill file on the fly */
int ** kill_tmp, *kill_vector; /* declaration of variables for trace killing */ create_trkill_table(kill_tmp,ntr_glob,recpos,nsrc_glob,srcpos,ishot,TRKILL_OFFSET_LOWER,TRKILL_OFFSET_UPPER);
char trace_kill_file[STRING_SIZE];
FILE *ftracekill;
if(TRKILL==1){ } else {
kill_tmp = imatrix(1,ntr_glob,1,nsrc_glob); /* READ TraceKill file from disk */
kill_vector = ivector(1,ntr);
if(USE_WORKFLOW){ if(USE_WORKFLOW){
sprintf(trace_kill_file,"%s_%i.dat",TRKILL_FILE,WORKFLOW_STAGE); sprintf(trace_kill_file,"%s_%i.dat",TRKILL_FILE,WORKFLOW_STAGE);
...@@ -81,15 +99,21 @@ if(TRKILL==1){ ...@@ -81,15 +99,21 @@ if(TRKILL==1){
fclose(ftracekill); fclose(ftracekill);
}
/* Generate local kill vector */
h=1; h=1;
for(i=1;i<=ntr;i++){ for(i=1;i<=ntr;i++){
kill_vector[h] = kill_tmp[recpos_loc[3][i]][ishot]; kill_vector[h] = kill_tmp[recpos_loc[3][i]][ishot];
h++; h++;
} }
} /* end if(TRKILL)*/ }
/*---------------*/
/* End Tracekill */
/*---------------*/
/* Integration of measured data */ /* Integration of measured data */
for(i=1;i<=ntr;i++){ for(i=1;i<=ntr;i++){
intseis_sd=0; intseis_sd=0;
...@@ -104,37 +128,37 @@ for(i=1;i<=ntr;i++){ ...@@ -104,37 +128,37 @@ for(i=1;i<=ntr;i++){
intseis_sectiondata[i][j]=sectiondata[i][j]; intseis_sectiondata[i][j]=sectiondata[i][j];
} }
} }
} }
/* TIME WINDOWING */ /* TIME WINDOWING */
if(TIMEWIN==1) time_window(intseis_sectiondata, iter, ntr_glob,recpos_loc, ntr, ns, ishot); if(TIMEWIN==1) time_window(intseis_sectiondata, iter, ntr_glob,recpos_loc, ntr, ns, ishot);
/* NORMALIZE TRACES */ /* NORMALIZE TRACES */
if(NORMALIZE==1) normalize_data(intseis_sectiondata,ntr,ns); if(NORMALIZE==1) normalize_data(intseis_sectiondata,ntr,ns);
/* calculate energy of measured data */ /* calculate energy of measured data */
for(i=1;i<=ntr;i++){ for(i=1;i<=ntr;i++){
if((TRKILL==1)&&(kill_vector[i]==1)) if((TRKILL==1)&&(kill_vector[i]==1))
continue; continue;
for(j=1;j<=ns;j++){ for(j=1;j<=ns;j++){
energy += intseis_sectiondata[i][j]*intseis_sectiondata[i][j]; energy += intseis_sectiondata[i][j]*intseis_sectiondata[i][j];
} }
} }
energy_dummy=energy; energy_dummy=energy;
/* printf("\n MYID = %i IN CALC_ENERGY: energy = %10.12f \n",MYID,energy); */ /* printf("\n MYID = %i IN CALC_ENERGY: energy = %10.12f \n",MYID,energy); */
/* free memory for integrated seismogram */ /* free memory for integrated seismogram */
free_matrix(intseis_sectiondata,1,ntr,1,ns); free_matrix(intseis_sectiondata,1,ntr,1,ns);
/* free memory for time windowing and trace killing */ /* free memory for time windowing and trace killing */
if(TIMEWIN==1) free_vector(picked_times,1,ntr); if(TIMEWIN==1) free_vector(picked_times,1,ntr);
if(TRKILL==1){ if(TRKILL==1){
free_imatrix(kill_tmp,1,ntr_glob,1,nsrc_glob); free_imatrix(kill_tmp,1,ntr_glob,1,nsrc_glob);
free_ivector(kill_vector,1,ntr);} free_ivector(kill_vector,1,ntr);}
return energy_dummy; return energy_dummy;
} /* end of function */ } /* end of function */
...@@ -14,7 +14,7 @@ ...@@ -14,7 +14,7 @@
* *
* You should have received a copy of the GNU General Public License * You should have received a copy of the GNU General Public License
* along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>. * along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/ -----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------ /*------------------------------------------------------------------------
* Calculate Misfit * Calculate Misfit
...@@ -22,39 +22,59 @@ ...@@ -22,39 +22,59 @@
* ----------------------------------------------------------------------*/ * ----------------------------------------------------------------------*/
#include "fd.h" #include "fd.h"
double calc_misfit(float **sectiondata, float **section, int ntr, int ns, int LNORM, float L2, int itest, int sws, int swstestshot,int ntr_glob, int **recpos_loc, int nsrc_glob, int ishot, int iter){ double calc_misfit(float **sectiondata, float **section, int ntr, int ns, int LNORM, float L2, int itest, int sws, int swstestshot,int ntr_glob, int **recpos_loc, int nsrc_glob, int ishot, int iter, float ** srcpos, int ** recpos){
/* declaration of variables */ /* declaration of variables */
extern float DT; extern float DT;
extern int MYID, USE_WORKFLOW, WORKFLOW_STAGE; extern int MYID, USE_WORKFLOW, WORKFLOW_STAGE;
extern int TRKILL, NORMALIZE, FC, TIMEWIN; extern int TRKILL, NORMALIZE, FC, TIMEWIN;
extern char TRKILL_FILE[STRING_SIZE]; extern char TRKILL_FILE[STRING_SIZE];
extern int VELOCITY; extern int VELOCITY;
int i,j;
float l2; int i,j;
int umax=0, h; float l2;
float abs_section, abs_sectiondata; int umax=0, h;
float intseis_s, intseis_sd; float abs_section, abs_sectiondata;
float *picked_times=NULL; float intseis_s, intseis_sd;
float **intseis_section=NULL, **intseis_sectiondata=NULL; float *picked_times=NULL;
float **intseis_sectiondata_envelope=NULL, **intseis_section_envelope=NULL; float **intseis_section=NULL, **intseis_sectiondata=NULL;
if(LNORM==8){ float **intseis_sectiondata_envelope=NULL, **intseis_section_envelope=NULL;
if(LNORM==8){
intseis_section_envelope = matrix(1,ntr,1,ns); intseis_section_envelope = matrix(1,ntr,1,ns);
intseis_sectiondata_envelope = matrix(1,ntr,1,ns); intseis_sectiondata_envelope = matrix(1,ntr,1,ns);
} }
intseis_section = matrix(1,ntr,1,ns); /* declaration of variables for integration */ intseis_section = matrix(1,ntr,1,ns); /* declaration of variables for integration */
intseis_sectiondata = matrix(1,ntr,1,ns); intseis_sectiondata = matrix(1,ntr,1,ns);
if(TIMEWIN) picked_times = vector(1,ntr); /* declaration of variables for TIMEWIN */ if(TIMEWIN) picked_times = vector(1,ntr); /* declaration of variables for TIMEWIN */
/* TRACE KILLING */ /* TRACE KILLING */
int ** kill_tmp, *kill_vector; /* declaration of variables for trace killing */ int ** kill_tmp, *kill_vector; /* declaration of variables for trace killing */
char trace_kill_file[STRING_SIZE]; char trace_kill_file[STRING_SIZE];
FILE *ftracekill; FILE *ftracekill;
extern int TRKILL_OFFSET;
if(TRKILL){ extern float TRKILL_OFFSET_LOWER;
extern float TRKILL_OFFSET_UPPER;
/*-----------------*/
/* Start Tracekill */
/*-----------------*/
if(TRKILL){
kill_tmp = imatrix(1,ntr_glob,1,nsrc_glob); kill_tmp = imatrix(1,ntr_glob,1,nsrc_glob);
kill_vector = ivector(1,ntr); kill_vector = ivector(1,ntr);
if(TRKILL_OFFSET) {
if(MYID==0) {
printf("Automatic offset based TraceKill: Experimental feature");
}
/* Generate TraceKill file on the fly */
create_trkill_table(kill_tmp,ntr_glob,recpos,nsrc_glob,srcpos,ishot,TRKILL_OFFSET_LOWER,TRKILL_OFFSET_UPPER);
} else {
/* READ TraceKill file from disk */
if(USE_WORKFLOW){ if(USE_WORKFLOW){
sprintf(trace_kill_file,"%s_%i.dat",TRKILL_FILE,WORKFLOW_STAGE); sprintf(trace_kill_file,"%s_%i.dat",TRKILL_FILE,WORKFLOW_STAGE);
ftracekill=fopen(trace_kill_file,"r"); ftracekill=fopen(trace_kill_file,"r");
...@@ -81,15 +101,21 @@ if(TRKILL){ ...@@ -81,15 +101,21 @@ if(TRKILL){
fclose(ftracekill); fclose(ftracekill);
}
/* Generate local kill vector */
h=1; h=1;
for(i=1;i<=ntr;i++){ for(i=1;i<=ntr;i++){
kill_vector[h] = kill_tmp[recpos_loc[3][i]][ishot]; kill_vector[h] = kill_tmp[recpos_loc[3][i]][ishot];
h++; h++;
} }
} /* end if(TRKILL)*/ }
/*---------------*/
/* End Tracekill */
/*---------------*/
/* Integration of measured and synthetic data */ /* Integration of measured and synthetic data */
for(i=1;i<=ntr;i++){ for(i=1;i<=ntr;i++){
intseis_s=0; intseis_s=0;
intseis_sd=0; intseis_sd=0;
...@@ -108,30 +134,30 @@ for(i=1;i<=ntr;i++){ ...@@ -108,30 +134,30 @@ for(i=1;i<=ntr;i++){
} }
} }
} }
/* TIME WINDOWING */ /* TIME WINDOWING */
if(TIMEWIN==1){ if(TIMEWIN==1){
time_window(intseis_section, iter, ntr_glob,recpos_loc, ntr, ns, ishot); time_window(intseis_section, iter, ntr_glob,recpos_loc, ntr, ns, ishot);
time_window(intseis_sectiondata, iter, ntr_glob,recpos_loc, ntr, ns, ishot); time_window(intseis_sectiondata, iter, ntr_glob,recpos_loc, ntr, ns, ishot);
} }
/* NORMALIZE TRACES */ /* NORMALIZE TRACES */
if(NORMALIZE==1){ if(NORMALIZE==1){
normalize_data(intseis_section,ntr,ns); normalize_data(intseis_section,ntr,ns);
normalize_data(intseis_sectiondata,ntr,ns); normalize_data(intseis_sectiondata,ntr,ns);
} }
/* calculate envelope for LNORM==8 */ /* calculate envelope for LNORM==8 */
if (LNORM==8){ if (LNORM==8){
calc_envelope(intseis_sectiondata,intseis_sectiondata_envelope,ns,ntr); calc_envelope(intseis_sectiondata,intseis_sectiondata_envelope,ns,ntr);
calc_envelope(intseis_section,intseis_section_envelope,ns,ntr); calc_envelope(intseis_section,intseis_section_envelope,ns,ntr);
} }
/* end of calculate envelope for LNORM==8 */ /* end of calculate envelope for LNORM==8 */
/* calculate kind of "energy" */ /* calculate kind of "energy" */
for(i=1;i<=ntr;i++){ for(i=1;i<=ntr;i++){
if((TRKILL==1)&&(kill_vector[i]==1)) if((TRKILL==1)&&(kill_vector[i]==1))
continue; continue;
...@@ -162,24 +188,24 @@ for(i=1;i<=ntr;i++){ ...@@ -162,24 +188,24 @@ for(i=1;i<=ntr;i++){
L2 += (intseis_section_envelope[i][j]-intseis_sectiondata_envelope[i][j]) * (intseis_section_envelope[i][j]-intseis_sectiondata_envelope[i][j]);} L2 += (intseis_section_envelope[i][j]-intseis_sectiondata_envelope[i][j]) * (intseis_section_envelope[i][j]-intseis_sectiondata_envelope[i][j]);}
} }
} }
l2=L2; l2=L2;
/* printf("\n MYID = %i IN CALC_MISFIT: L2 = %10.12f \n",MYID,l2); */ /* printf("\n MYID = %i IN CALC_MISFIT: L2 = %10.12f \n",MYID,l2); */
/* free memory for integrated seismograms */ /* free memory for integrated seismograms */
free_matrix(intseis_section,1,ntr,1,ns); free_matrix(intseis_section,1,ntr,1,ns);
free_matrix(intseis_sectiondata,1,ntr,1,ns); free_matrix(intseis_sectiondata,1,ntr,1,ns);
/* free memory for time windowing and trace killing */ /* free memory for time windowing and trace killing */
if(TIMEWIN==1) free_vector(picked_times,1,ntr); if(TIMEWIN==1) free_vector(picked_times,1,ntr);
if(TRKILL==1){ if(TRKILL==1){
free_imatrix(kill_tmp,1,ntr_glob,1,nsrc_glob); free_imatrix(kill_tmp,1,ntr_glob,1,nsrc_glob);
free_ivector(kill_vector,1,ntr);} free_ivector(kill_vector,1,ntr);}
if(LNORM==8){ if(LNORM==8){
free_matrix(intseis_sectiondata_envelope,1,ntr,1,ns); free_matrix(intseis_sectiondata_envelope,1,ntr,1,ns);
free_matrix(intseis_section_envelope,1,ntr,1,ns); free_matrix(intseis_section_envelope,1,ntr,1,ns);
} }
return l2; return l2;
} /* end of function */ } /* end of function */
...@@ -14,7 +14,7 @@ ...@@ -14,7 +14,7 @@
* *
* You should have received a copy of the GNU General Public License * You should have received a copy of the GNU General Public License
* along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>. * along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/ -----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------ /*------------------------------------------------------------------------
* Calculate Data Residuals * Calculate Data Residuals
...@@ -22,52 +22,67 @@ ...@@ -22,52 +22,67 @@
* ----------------------------------------------------------------------*/ * ----------------------------------------------------------------------*/
#include "fd.h" #include "fd.h"
double calc_res(float **sectiondata, float **section, float **sectiondiff, float **sectiondiffold, int ntr, int ns, int LNORM, float L2, int itest, int sws, int swstestshot, int ntr_glob, int **recpos_loc, int nsrc_glob, int ishot, int iter){ double calc_res(float **sectiondata, float **section, float **sectiondiff, float **sectiondiffold, int ntr, int ns, int LNORM, float L2, int itest, int sws, int swstestshot, int ntr_glob, int **recpos_loc, int nsrc_glob, int ishot, int iter, float ** srcpos, int ** recpos){
/* declaration of variables */ /* declaration of variables */
extern float DT, WATERLEVEL_LNORM8; extern float DT, WATERLEVEL_LNORM8;
extern int REC1, REC2, MYID, ACOUSTIC; extern int REC1, REC2, MYID, ACOUSTIC;
extern int TRKILL, NORMALIZE, FC, TIMEWIN; extern int TRKILL, NORMALIZE, FC, TIMEWIN;
extern char TRKILL_FILE[STRING_SIZE]; extern char TRKILL_FILE[STRING_SIZE];
extern int VELOCITY, USE_WORKFLOW, WORKFLOW_STAGE; extern int VELOCITY, USE_WORKFLOW, WORKFLOW_STAGE;
float RMS, signL1, intseis; float RMS, signL1, intseis;
int Lcount,i,j,invtime,k,h, umax=0; int Lcount,i,j,invtime,k,h, umax=0;
float l2; float l2;
float abs_section, abs_sectiondata, sectiondata_mult_section; float abs_section, abs_sectiondata, sectiondata_mult_section;
float intseis_s, intseis_sd; float intseis_s, intseis_sd;
float *picked_times=NULL; float *picked_times=NULL;
float **intseis_section=NULL, **intseis_sectiondata=NULL; float **intseis_section=NULL, **intseis_sectiondata=NULL;
float **intseis_sectiondata_envelope=NULL, **intseis_section_envelope=NULL, **intseis_section_hilbert=NULL, **dummy_1=NULL, **dummy_2=NULL; float **intseis_sectiondata_envelope=NULL, **intseis_section_envelope=NULL, **intseis_section_hilbert=NULL, **dummy_1=NULL, **dummy_2=NULL;
if(LNORM==8){ if(LNORM==8){
intseis_section_envelope = matrix(1,ntr,1,ns); intseis_section_envelope = matrix(1,ntr,1,ns);
intseis_sectiondata_envelope = matrix(1,ntr,1,ns); intseis_sectiondata_envelope = matrix(1,ntr,1,ns);
intseis_section_hilbert = matrix(1,ntr,1,ns); intseis_section_hilbert = matrix(1,ntr,1,ns);
dummy_1 = matrix(1,ntr,1,ns); dummy_1 = matrix(1,ntr,1,ns);
dummy_2 = matrix(1,ntr,1,ns); dummy_2 = matrix(1,ntr,1,ns);
} }
intseis_section = matrix(1,ntr,1,ns); /* declaration of variables for integration */
intseis_sectiondata = matrix(1,ntr,1,ns);
if(TIMEWIN) picked_times = vector(1,ntr); /* declaration of variables for TIMEWIN */
/* sectiondiff will be set to zero */
umax=ntr*ns;
zero(&sectiondiff[1][1],umax);
/* TRACE KILLING */ intseis_section = matrix(1,ntr,1,ns); /* declaration of variables for integration */
int ** kill_tmp, *kill_vector; /* declaration of variables for trace killing */ intseis_sectiondata = matrix(1,ntr,1,ns);
char trace_kill_file[STRING_SIZE]; if(TIMEWIN) picked_times = vector(1,ntr); /* declaration of variables for TIMEWIN */
FILE *ftracekill;
if(TRKILL==1){
/* sectiondiff will be set to zero */ /* sectiondiff will be set to zero */
umax=ntr*ns; umax=ntr*ns;
zero(&sectiondiff[1][1],umax); zero(&sectiondiff[1][1],umax);
/* TRACE KILLING */
int ** kill_tmp, *kill_vector; /* declaration of variables for trace killing */
char trace_kill_file[STRING_SIZE];
FILE *ftracekill;
extern int TRKILL_OFFSET;
extern float TRKILL_OFFSET_LOWER;