Commit 67c207a1 authored by laura.gassner's avatar laura.gassner

Added taper for inverted source time function (use switch TAPER_STF, see...

Added taper for inverted source time function (use switch TAPER_STF, see taper.c for definition) and individual time windowing (use switch TW_IND, then three columns are expected in PICK_FILE for each trace with picked time - twl plus - twl minus, in this case TWLENGTH_PLUS and TWLENGTH_MINUS are ignored). Time windowing is now also applied during the inversion of the source time function.
parent 30bd4c77
......@@ -200,6 +200,7 @@
"PARA" : "fdlsq:tshift=0.0",
"N_STF" : "10",
"N_STF_START" : "1",
"TAPER_STF" : "0",
"TRKILL_STF" : "0",
"TRKILL_FILE_STF" : "./trace_kill/trace_kill.dat",
......@@ -219,6 +220,7 @@
"Time windowing and damping" : "comment",
"TIMEWIN" : "0",
"TW_IND" : "0",
"PICKS_FILE" : "./picked_times/picks",
"TWLENGTH_PLUS" : "0.01",
"TWLENGTH_MINUS" : "0.01",
......
......@@ -124,9 +124,9 @@
"SRTRADIUS" : "3.0",
"FILTSIZE" : "1",
"SWS_TAPER_FILE" : "1",
"TAPER_FILE_NAME" : "taper/taper.bin",
"TAPER_FILE_NAME_U" : "taper/taper_u.bin",
"TAPER_FILE_NAME_RHO" : "taper/taper_rho.bin",
"TAPER_FILE_NAME" : "taper/taper_te_ac.bin",
"TAPER_FILE_NAME_U" : "taper/taper_te_ac_u.bin",
"TAPER_FILE_NAME_RHO" : "taper/taper_te_ac_rho.bin",
"Upper and lower limits for model parameters" : "comment",
"VPUPPERLIM" : "5000",
......@@ -145,6 +145,7 @@
"PARA" : "fdlsq:tshift=0.0",
"N_STF" : "10",
"N_STF_START" : "1",
"TAPER_STF" : "1",
"TRKILL_STF" : "0",
"TRKILL_FILE_STF" : "./trace_kill/trace_kill.dat",
......
......@@ -141,6 +141,7 @@ DENISE= \
timedomain_filt.c \
timedomain_filt_vector.c \
time_window.c \
time_window_glob.c \
util.c \
wavelet.c \
wavelet_stf.c \
......
......@@ -35,10 +35,8 @@ float energy_dummy, intseis;
int umax=0, h;
float intseis_sd;
float **intseis_sectiondata=NULL;
float *picked_times=NULL;
intseis_sectiondata = matrix(1,ntr,1,ns); /* declaration of variables for integration */
if(TIMEWIN) picked_times = vector(1,ntr); /* declaration of variables for TIMEWIN */
/*extern int MYID;*/
......@@ -93,7 +91,7 @@ for(i=1;i<=ntr;i++){
}
/* TIME WINDOWING */
if(TIMEWIN==1) time_window(intseis_sectiondata, picked_times, 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 */
if(NORMALIZE==1) normalize_data(intseis_sectiondata,ntr,ns);
......@@ -116,8 +114,7 @@ energy_dummy=energy;
/* free memory for integrated seismogram */
free_matrix(intseis_sectiondata,1,ntr,1,ns);
/* free memory for time windowing and trace killing */
if(TIMEWIN==1) free_vector(picked_times,1,ntr);
/* free memory for trace killing */
if(TRKILL==1){
free_imatrix(kill_tmp,1,ntr_glob,1,nsrc_glob);
free_ivector(kill_vector,1,ntr);}
......
......@@ -35,7 +35,6 @@ float l2;
int umax=0, h;
float abs_section, abs_sectiondata;
float intseis_s, intseis_sd;
float *picked_times=NULL;
float **intseis_section=NULL, **intseis_sectiondata=NULL;
float **intseis_sectiondata_envelope=NULL, **intseis_section_envelope=NULL;
if(LNORM==8){
......@@ -44,7 +43,6 @@ if(LNORM==8){
}
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 */
/* TRACE KILLING */
int ** kill_tmp, *kill_vector; /* declaration of variables for trace killing */
......@@ -98,8 +96,8 @@ for(i=1;i<=ntr;i++){
/* TIME WINDOWING */
if(TIMEWIN==1){
time_window(intseis_section, picked_times, iter, ntr_glob,recpos_loc, ntr, ns, ishot);
time_window(intseis_sectiondata, picked_times, 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);
}
/* NORMALIZE TRACES */
......@@ -158,8 +156,8 @@ l2=L2;
free_matrix(intseis_section,1,ntr,1,ns);
free_matrix(intseis_sectiondata,1,ntr,1,ns);
/* free memory for time windowing and trace killing */
if(TIMEWIN==1) free_vector(picked_times,1,ntr);
/* free memory for trace killing */
if(TRKILL==1){
free_imatrix(kill_tmp,1,ntr_glob,1,nsrc_glob);
free_ivector(kill_vector,1,ntr);}
......
......@@ -35,7 +35,6 @@ int Lcount,i,j,invtime,k,h, umax=0;
float l2;
float abs_section, abs_sectiondata, sectiondata_mult_section;
float intseis_s, intseis_sd;
float *picked_times=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;
......@@ -49,7 +48,6 @@ if(LNORM==8){
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;
......@@ -111,8 +109,8 @@ for(i=1;i<=ntr;i++){
/* TIME WINDOWING */
if(TIMEWIN==1){
time_window(intseis_section, picked_times, iter, ntr_glob,recpos_loc, ntr, ns, ishot);
time_window(intseis_sectiondata, picked_times, 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);
}
/* NORMALIZE TRACES */
......@@ -238,16 +236,10 @@ for(i=1;i<=ntr;i++){
l2=L2;
// taper(sectiondiff,ntr,ns);
/* printf("\n MYID = %i IN CALC_RES: L2 = %10.12f \n",MYID,l2); */
/* free memory for integrated seismograms */
free_matrix(intseis_section,1,ntr,1,ns);
free_matrix(intseis_sectiondata,1,ntr,1,ns);
/* free memory for time windowing and trace killing */
if(TIMEWIN==1) free_vector(picked_times,1,ntr);
if(TRKILL==1){
free_imatrix(kill_tmp,1,ntr_glob,1,nsrc_glob);
free_ivector(kill_vector,1,ntr);
......
......@@ -50,40 +50,38 @@ int main(int argc, char **argv){
char *buff_addr, ext[10], *fileinp;
char wave_forward[225], wave_recipro[225], wave_conv[225], jac[225], jac2[225], jacsum[225], dwavelet[225], vyf[STRING_SIZE];
double time1, time2, time3, time4, time5, time6, time7, time8,
time_av_v_update=0.0, time_av_s_update=0.0, time_av_v_exchange=0.0,
time_av_s_exchange=0.0, time_av_timestep=0.0;
double time1, time2, time3, time4, time5, time6, time7, time8, time_av_v_update=0.0, time_av_s_update=0.0, time_av_v_exchange=0.0, time_av_s_exchange=0.0, time_av_timestep=0.0;
float L2, L2sum, L2_all_shots, L2sum_all_shots, *L2t, alphanomsum, alphanom, alphadenomsum, alphadenom, scaleamp ,sdummy, lamr;
int sum_killed_traces=0, sum_killed_traces_testshots=0, killed_traces=0, killed_traces_testshots=0;
int *ptr_killed_traces=&killed_traces, *ptr_killed_traces_testshots=&killed_traces_testshots;
float energy, energy_sum, energy_all_shots, energy_sum_all_shots;
float energy, energy_sum, energy_all_shots, energy_sum_all_shots;
float ** psxx, ** psxy, ** psyy, **psp, ** ux, ** uy, ** uxy, ** uyx, ** Vp0, ** uttx, ** utty, ** Vs0, ** Rho0;
float ** pvx, ** pvy, **waveconv, **waveconv_lam, **waveconv_mu, **waveconv_rho, **waveconv_rho_s, **waveconv_u, **waveconvtmp, **wcpart, **wavejac;
float **waveconv_shot, **waveconv_u_shot, **waveconv_rho_shot;
float ** pvxp1, ** pvyp1, ** pvxm1, ** pvym1;
float ** gradg, ** gradp,** gradg_rho, ** gradp_rho, ** gradg_u, ** gradp_u;
float ** prho,** prhonp1, **prip=NULL, **prjp=NULL, **pripnp1=NULL, **prjpnp1=NULL, ** ppi, ** pu, ** punp1, ** puipjp, ** ppinp1;
float ** vpmat, ***forward_prop_x, ***forward_prop_y, ***forward_prop_rho_x, ***forward_prop_u, ***forward_prop_rho_y, ***forward_prop_p;
float ** psxx, ** psxy, ** psyy, **psp, ** ux, ** uy, ** uxy, ** uyx, ** Vp0, ** uttx, ** utty, ** Vs0, ** Rho0;
float ** pvx, ** pvy, ** waveconv, **waveconv_lam, **waveconv_mu, **waveconv_rho, **waveconv_rho_s, **waveconv_u, **waveconvtmp, **wcpart, **wavejac;
float ** waveconv_shot, ** waveconv_u_shot, **waveconv_rho_shot;
float ** pvxp1, ** pvyp1, ** pvxm1, ** pvym1;
float ** gradg, ** gradp, ** gradg_rho, ** gradp_rho, ** gradg_u, ** gradp_u;
float ** prho, ** prhonp1, ** prip=NULL, **prjp=NULL, **pripnp1=NULL, **prjpnp1=NULL, ** ppi, ** pu, ** punp1, ** puipjp, ** ppinp1;
float ** vpmat, ***forward_prop_x, ***forward_prop_y, ***forward_prop_rho_x, ***forward_prop_u, ***forward_prop_rho_y, ***forward_prop_p;
float ** sectionvx=NULL, ** sectionvy=NULL, ** sectionp=NULL, ** sectionpnp1=NULL,
** sectioncurl=NULL, ** sectiondiv=NULL, ** sectionvxdata=NULL, ** sectionvxdiff=NULL, ** sectionvxdiffold=NULL, ** sectionvydiffold=NULL, ** sectionpdata=NULL, ** sectionpdiff=NULL, ** sectionpdiffold=NULL,
** sectionvydiff=NULL, ** sectionvydata=NULL, ** sectionpn=NULL, ** sectionread=NULL, ** sectionvy_conv=NULL, ** sectionvy_obs=NULL, ** sectionvx_conv=NULL,** sectionvx_obs=NULL,
** sectionp_conv=NULL,** sectionp_obs=NULL, * source_time_function=NULL;
float ** sectionvx=NULL, ** sectionvy=NULL, ** sectionp=NULL, ** sectionpnp1=NULL, ** sectioncurl=NULL, ** sectiondiv=NULL, ** sectionvxdata=NULL, ** sectionvxdiff=NULL, ** sectionvxdiffold=NULL,
** sectionvydiffold=NULL, ** sectionpdata=NULL, ** sectionpdiff=NULL, ** sectionpdiffold=NULL, ** sectionvydiff=NULL, ** sectionvydata=NULL, ** sectionpn=NULL, ** sectionread=NULL,
** sectionvy_conv=NULL, ** sectionvy_obs=NULL, ** sectionvx_conv=NULL,** sectionvx_obs=NULL, ** sectionp_conv=NULL,** sectionp_obs=NULL, * source_time_function=NULL;
float ** absorb_coeff, ** taper_coeff, * epst1, * epst2, * epst3, * picked_times;
float ** srcpos=NULL, **srcpos_loc=NULL, ** srcpos1=NULL, **srcpos_loc_back=NULL, ** signals=NULL, ** signals_rec=NULL, *hc=NULL, ** dsignals=NULL;
int ** recpos=NULL, ** recpos_loc=NULL;
/*int ** tracekill=NULL, TRKILL, DTRKILL;*/
int * DTINV_help;
float ** bufferlef_to_rig, ** bufferrig_to_lef, ** buffertop_to_bot, ** bufferbot_to_top;
/* PML variables */
float * d_x, * K_x, * alpha_prime_x, * a_x, * b_x, * d_x_half, * K_x_half, * alpha_prime_x_half, * a_x_half, * b_x_half, * d_y, * K_y, * alpha_prime_y, * a_y, * b_y, * d_y_half, * K_y_half, * alpha_prime_y_half, * a_y_half, * b_y_half;
float ** psi_sxx_x, ** psi_syy_y, ** psi_sxy_y, ** psi_sxy_x, ** psi_vxx, ** psi_vyy, ** psi_vxy, ** psi_vyx, ** psi_vxxs;
float * d_x, * K_x, * alpha_prime_x, * a_x, * b_x, * d_x_half, * K_x_half, * alpha_prime_x_half, * a_x_half, * b_x_half,
* d_y, * K_y, * alpha_prime_y, * a_y, * b_y, * d_y_half, * K_y_half, * alpha_prime_y_half, * a_y_half, * b_y_half;
float ** psi_sxx_x, ** psi_syy_y, ** psi_sxy_y, ** psi_sxy_x, ** psi_vxx, ** psi_vyy, ** psi_vxy, ** psi_vyx, ** psi_vxxs;
/* Variables for viscoelastic modeling */
float **ptaus=NULL, **ptaup=NULL, *etaip=NULL, *etajm=NULL, *peta=NULL, **ptausipjp=NULL, **fipjp=NULL, ***dip=NULL, *bip=NULL, *bjm=NULL;
......
......@@ -71,7 +71,8 @@ void exchange_par(void){
extern float PRO;
extern int TRKILL, TRKILL_STF;
extern char TRKILL_FILE[STRING_SIZE], TRKILL_FILE_STF[STRING_SIZE];
extern int TIMEWIN, NORMALIZE;
extern int TAPER_STF;
extern int TIMEWIN, NORMALIZE, TW_IND;
extern float TWLENGTH_PLUS, TWLENGTH_MINUS, GAMMA;
extern char PICKS_FILE[STRING_SIZE];
extern char MISFIT_LOG_FILE[STRING_SIZE];
......@@ -289,6 +290,9 @@ void exchange_par(void){
idum[91] = ACOUSTIC;
idum[92] = TW_IND;
idum[93] = TAPER_STF;
} /** if (MYID == 0) **/
......@@ -524,6 +528,10 @@ void exchange_par(void){
GRAD_FILT_WAVELENGTH = idum[90];
ACOUSTIC = idum[91];
TW_IND = idum[92];
TAPER_STF = idum[93];
MPI_Bcast(&FL[1],L,MPI_FLOAT,0,MPI_COMM_WORLD);
......
......@@ -195,7 +195,7 @@ void inseis(FILE *fp, int comp, float **section, int ntr, int ns, int sws, int
void inseis_source_wavelet(float *section, int ns, int ishot);
void taper(float **sectionpdiff, int ntr, int ns);
void taper(float *section, int ns, float fc);
void output_source_signal(FILE *fp, float **signals, int ns, int seis_form);
......@@ -303,7 +303,8 @@ float *etajm, float *peta, float * hc, float * K_x, float * a_x, float * b_x, fl
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 time_window(float **sectiondata, float * picked_times, int iter, int ntr_glob, int **recpos_loc, int ntr, int ns, int ishot);
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);
void prepare_update_s(float *etajm, float *etaip, float *peta, float **fipjp, float **pu,
float **puipjp, float **ppi, float **prho, float **ptaus, float **ptaup,
......
......@@ -84,7 +84,9 @@ char TRKILL_FILE[STRING_SIZE];
int TRKILL_STF;
char TRKILL_FILE_STF[STRING_SIZE];
int TIMEWIN, NORMALIZE;
int TAPER_STF;
int TIMEWIN, NORMALIZE, TW_IND;
float TWLENGTH_PLUS, TWLENGTH_MINUS, GAMMA;
char PICKS_FILE[STRING_SIZE];
char MISFIT_LOG_FILE[STRING_SIZE];
......
......@@ -80,7 +80,9 @@ void read_par_json(FILE *fp, char *fileinp){
extern int TRKILL_STF;
extern char TRKILL_FILE_STF[STRING_SIZE];
extern int TIMEWIN, NORMALIZE;
extern int TAPER_STF;
extern int TIMEWIN, NORMALIZE, TW_IND;
extern float TWLENGTH_PLUS, TWLENGTH_MINUS, GAMMA;
extern char PICKS_FILE[STRING_SIZE];
......@@ -660,6 +662,10 @@ void read_par_json(FILE *fp, char *fileinp){
}
}
/* Taper STF */
if (get_int_from_objectlist("TAPER_STF",number_readobjects,&TAPER_STF,varname_list, value_list)){
TAPER_STF=0;
fprintf(fp,"Variable TAPER_STF is set to default value %d.\n",TAPER_STF);}
/* Frequency filtering during inversion */
if (get_int_from_objectlist("TIME_FILT",number_readobjects,&TIME_FILT,varname_list, value_list)){
......@@ -746,6 +752,9 @@ void read_par_json(FILE *fp, char *fileinp){
fprintf(fp,"Variable TIMEWIN is set to default value %d.\n",TIMEWIN);}
else {
if (TIMEWIN==1){
if (get_int_from_objectlist("TW_IND",number_readobjects,&TW_IND,varname_list, value_list)){
TW_IND=0;
fprintf(fp,"Variable TW_IND is set to default value %d.\n",TW_IND);}
if (get_string_from_objectlist("PICKS_FILE",number_readobjects,PICKS_FILE,varname_list, value_list))
err("Variable PICKS_FILE could not be retrieved from the json input file!");
if (get_float_from_objectlist("TWLENGTH_PLUS",number_readobjects,&TWLENGTH_PLUS,varname_list, value_list))
......
......@@ -30,7 +30,7 @@ int ntr_glob,int ntr, float ** srcpos, int ishot, int ns, int iter, int nshots,
/* declaration of global variables */
extern float DT, DH;
extern int SEIS_FORMAT, MYID, NT, QUELLART, TIME_FILT, TIMEWIN;
extern int SEIS_FORMAT, MYID, NT, QUELLART, TIME_FILT, TIMEWIN, TAPER_STF;
extern char SEIS_FILE_VY[STRING_SIZE], SEIS_FILE_P[STRING_SIZE], PARA[STRING_SIZE], DATA_DIR[STRING_SIZE];
extern int TRKILL_STF, NORMALIZE;
extern char TRKILL_FILE_STF[STRING_SIZE];
......@@ -38,11 +38,9 @@ int ntr_glob,int ntr, float ** srcpos, int ishot, int ns, int iter, int nshots,
/* declaration of variables for trace killing */
int ** kill_tmp, *kill_vector, h, j;
char trace_kill_file[STRING_SIZE];
char trace_kill_file[STRING_SIZE];
FILE *ftracekill;
float *picked_times=NULL;
/* --------------- declaration of variables --------------- */
unsigned int nrec, nsamp, i, npairs;
float dt;
......@@ -58,8 +56,6 @@ int ntr_glob,int ntr, float ** srcpos, int ishot, int ns, int iter, int nshots,
wavelet=vector(1,ns);
stf_conv_wavelet=vector(1,ns);
if(TIMEWIN) picked_times = vector(1,ntr); /* declaration of variables for TIMEWIN */
printf("\n================================================================================================\n\n");
printf("\n ***** Inversion of Source Time Function - shot: %d - it: %d ***** \n\n",ishot,iter);
......@@ -82,10 +78,10 @@ int ntr_glob,int ntr, float ** srcpos, int ishot, int ns, int iter, int nshots,
h=1;
for(i=1;i<=ntr_glob;i++){
kill_vector[h] = kill_tmp[i][ishot];
h++;
kill_vector[h] = kill_tmp[i][ishot];
h++;
}
} /* end if(TRKILL_STF)*/
} /* end if(TRKILL_STF)*/
if(TRKILL_STF){
for(i=1;i<=ntr_glob;i++){
......@@ -93,25 +89,25 @@ int ntr_glob,int ntr, float ** srcpos, int ishot, int ns, int iter, int nshots,
if(kill_vector[i]==1){
printf("%d \t",i);
for(j=1;j<=ns;j++){
sectionvy[i][j]=0.0;
sectionvy_obs[i][j]=0.0;
}
}
sectionvy[i][j]=0.0;
sectionvy_obs[i][j]=0.0;
}
}
if(i==ntr_glob)printf(" ***** \n\n");
}
}
/* trace killing ends here */
if(TIMEWIN==1){
time_window(sectionvy, picked_times, iter, ntr_glob,recpos_loc, ntr, ns, ishot);
time_window(sectionvy_obs, picked_times, iter, ntr_glob,recpos_loc, ntr, ns, ishot);
time_window_glob(sectionvy, iter, ntr_glob, ns, ishot);
time_window_glob(sectionvy_obs, iter, ntr_glob, ns, ishot);
}
/* NORMALIZE TRACES */
/*if(NORMALIZE==1){*/
normalize_data(sectionvy,ntr_glob,ns);
normalize_data(sectionvy_obs,ntr_glob,ns);
/*}*/
if(NORMALIZE==1){
normalize_data(sectionvy,ntr_glob,ns);
normalize_data(sectionvy_obs,ntr_glob,ns);
}
nrec=(unsigned int)ntr_glob;
nsamp=(unsigned int)ns;
......@@ -246,6 +242,8 @@ int ntr_glob,int ntr, float ** srcpos, int ishot, int ns, int iter, int nshots,
/* convolving wavelet with STF */
conv_FD(wavelet,source_time_function,stf_conv_wavelet,ns);
if(TAPER_STF)
taper(stf_conv_wavelet, ns, FC);
// /* --------------- writing out the observed seismograms --------------- */
// sprintf(obs_y_tmp,"%s.shot%d.obs",SEIS_FILE_P,ishot);
......@@ -280,9 +278,6 @@ int ntr_glob,int ntr, float ** srcpos, int ishot, int ns, int iter, int nshots,
free_vector(stf_conv_wavelet,1,ns);
free_vector(psource,1,ns);
/* free memory for time windowing and trace killing */
if(TIMEWIN==1) free_vector(picked_times,1,ntr);
/* free memory for trace killing */
if(TRKILL_STF){
free_imatrix(kill_tmp,1,ntr_glob,1,nshots);
......
......@@ -17,20 +17,21 @@
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* Write seismograms to disk
* last update 19/01/02, T. Bohlen
* Taper function now adapted for use in stf.c, only for single traces!
* last update 21/10/15, L. Gassner
* ----------------------------------------------------------------------*/
#include "fd.h"
#include "segy.h"
void taper(float **sectionpdiff, int ntr, int ns){
void taper(float *section, int ns, float fc){
/* declaration of extern variables */
extern int MYID, TAPERLENGTH;
extern int MYID;
extern float DT;
extern FILE *FP;
/* declaration of local variables */
int i,j, h;
int i,j, h, taperlength, taperduration;
int tracl1;
float a;
float damping, amp;
......@@ -39,31 +40,47 @@ void taper(float **sectionpdiff, int ntr, int ns){
window = vector(1,ns);
amp1 = vector(1,ns);
taperlength=(int)(ceil(2.0/fc/DT));
taperduration=2*taperlength;
/* "Cerjan"-Window */
damping=100.0;
damping=99.9;
amp=1.0-damping/100.0;
a=sqrt(-log(amp)/((TAPERLENGTH-1)*(TAPERLENGTH-1)));
for (i=1;i<=ns;i++){
a=sqrt(-log(amp)/((taperlength-1)*(taperlength-1)));
for (i=1;i<=ns;i++){
window[i]=1.0;
amp1[i]=0.0;
}
if (MYID==0){
fprintf(FP,"\n ntr: %d\n",ntr);
fprintf(FP,"\n ns: %d\n",ns);
fprintf(FP,"\n TAPERLENGTH: %d\n",TAPERLENGTH);
fprintf(FP,"\n fc: %f\n",fc);
fprintf(FP,"\n taperlength: %d\n",taperlength);
}
for (i=1;i<=taperlength;i++){
amp1[i]=exp(-(a*a*(taperlength-i)*(taperlength-i)));
}
// /* Taper at the beginning of the window*/
// for (i=1;i<=taperlength;i++){
// window[i]=amp1[i];
// }
h=1;
for (i=taperduration;i>=(taperduration-taperlength+3);i--){
window[i]=amp1[h];
h++;
}
for (i=1;i<=TAPERLENGTH;i++){amp1[i]=exp(-(a*a*(TAPERLENGTH-i)*(TAPERLENGTH-i)));}
for (i=1;i<=1;i++){window[i]=amp1[i];}
h=1;
for (i=ns;i>=(ns-TAPERLENGTH+3);i--){window[i]=amp1[h];h++;}
/* Apply taper window */
for(tracl1=1;tracl1<=ntr;tracl1++){
for(j=1;j<=ns;j++){
sectionpdiff[tracl1][j]*=window[j];
}
for (i=taperduration;i<=ns;i++){
window[i]=amp1[i];
}
for(j=1;j<=ns;j++){
section[j]*=window[j];
}
free_vector(window,1,ns);
free_vector(amp1,1,ns);
}
......@@ -23,76 +23,126 @@
* ----------------------------------------------------------------------*/
#include "fd.h"
void time_window(float **sectiondata, float *picked_times, int iter, int ntr_glob, int **recpos_loc, int ntr, int ns, int ishot){
void time_window(float **sectiondata, int iter, int ntr_glob, int **recpos_loc, int ntr, int ns, int ishot){
/* declaration of variables */
extern float DT;
extern int REC1, REC2, MYID;
extern float GAMMA, TWLENGTH_PLUS, TWLENGTH_MINUS;
extern int TW_IND;
extern char PICKS_FILE[STRING_SIZE];
int READ_PICKED_TIMES;
char pickfile_char[STRING_SIZE];
float time, dump, dump1, taper, taper1;
float *pick_tmp;
float time, dump, dump1, dump2, dump3, taper, taper1;
float *pick_tmp, **pick_tmp_m;
int i, j, h;
float *picked_times=NULL, **picked_times_m=NULL;
FILE *fptime;
READ_PICKED_TIMES=1; /*other options?*/
/* read picked first arrival times */
if(READ_PICKED_TIMES==1){
if(TW_IND){
picked_times_m = matrix(1,3,1,ntr_glob);
pick_tmp_m = matrix(1,3,1,ntr_glob);
}else{
picked_times = vector(1,ntr_glob);
pick_tmp = vector(1,ntr_glob);
sprintf(pickfile_char,"%s_%i.dat",PICKS_FILE,ishot);
fptime=fopen(pickfile_char,"r");
if (fptime == NULL) {
err(" picks_?.dat could not be opened !");
}
sprintf(pickfile_char,"%s_%i.dat",PICKS_FILE,ishot);
fptime=fopen(pickfile_char,"r");
if (fptime == NULL) {
err(" picks_?.dat could not be opened !");
}
if(TW_IND){
for(i=1;i<=ntr_glob;i++){
fscanf(fptime,"%f%f%f",&dump,&dump2,&dump3);
pick_tmp_m[1][i] = dump;
pick_tmp_m[2][i] = dump2;
pick_tmp_m[3][i] = dump3;
}
}else{
for(i=1;i<=ntr_glob;i++){
fscanf(fptime,"%f",&dump);
pick_tmp[i] = dump;
}
fclose(fptime);
/* distribute picks on CPUs */
h=1;
}
fclose(fptime);
/* distribute picks on CPUs */
h=1;
if(TW_IND){
for(i=1;i<=ntr;i++){
picked_times[h] = pick_tmp[recpos_loc[3][i]];
picked_times_m[1][h] = pick_tmp_m[1][recpos_loc[3][i]];
picked_times_m[2][h] = pick_tmp_m[2][recpos_loc[3][i]];
picked_times_m[3][h] = pick_tmp_m[3][recpos_loc[3][i]];
/*printf("MYID IS: %i REC1 is= %i REC2 is %i\n",MYID, REC1, REC2);
printf("MYID= %i pick_tmp[%i]= %f\n",MYID, i,pick_tmp[i]);*/
h++;
}
}else{
for(i=1;i<=ntr;i++){
picked_times[h] = pick_tmp[recpos_loc[3][i]];
h++;
}
}
if(TW_IND)
free_matrix(pick_tmp_m,1,3,1,ntr_glob);
else
free_vector(pick_tmp,1,ntr_glob);
} /* end of if(READTIMES==1) */
for(i=1;i<=ntr;i++){
for(j=2;j<=ns;j++){
time = (float)(j * DT);
dump = (time-picked_times[i]-TWLENGTH_PLUS);
taper = exp(-GAMMA*dump*dump);
dump1 = (time-picked_times[i]+TWLENGTH_MINUS);
taper1 = exp(-GAMMA*dump1*dump1);
if(time>=picked_times[i]+TWLENGTH_PLUS){
sectiondata[i][j] = sectiondata[i][j] * taper;}
if(time<=picked_times[i]-TWLENGTH_MINUS){
sectiondata[i][j] = sectiondata[i][j] * taper1;}
if(TW_IND){
for(i=1;i<=ntr;i++){
for(j=2;j<=ns;j++){
sectiondata[i][j] = sectiondata[i][j];
time = (float)(j * DT);
dump = (time-picked_times_m[1][i]-picked_times_m[2][i]);
taper = exp(-GAMMA*dump*dump);
dump1 = (time-picked_times_m[1][i]+picked_times_m[3][i]);
taper1 = exp(-GAMMA*dump1*dump1);
if(time>=picked_times_m[1][i]+picked_times_m[2][i]){
sectiondata[i][j] = sectiondata[i][j] * taper;}
if(time<=picked_times_m[1][i]-picked_times_m[3][i]){
sectiondata[i][j] = sectiondata[i][j] * taper1;}
sectiondata[i][j] = sectiondata[i][j];
}
}
}else{
for(i=1;i<=ntr;i++){
for(j=2;j<=ns;j++){
}
time = (float)(j * DT);
dump = (time-picked_times[i]-TWLENGTH_PLUS);
taper = exp(-GAMMA*dump*dump);
dump1 = (time-picked_times[i]+TWLENGTH_MINUS);
taper1 = exp(-GAMMA*dump1*dump1);
if(time>=picked_times[i]+TWLENGTH_PLUS){
sectiondata[i][j] = sectiondata[i][j] * taper;}
if(time<=picked_times[i]-TWLENGTH_MINUS){
sectiondata[i][j] = sectiondata[i][j] * taper1;}
sectiondata[i][j] = sectiondata[i][j];
}
}
}
if(TW_IND)
free_matrix(picked_times_m,1,3,1,ntr_glob);
else
free_vector(picked_times,1,ntr_glob);
} /* end of function time_window.c */
/*-----------------------------------------------------------------------------------------
* Copyright (C) 2013 For the list of authors, see file AUTHORS.
*
* This file is part of DENISE.
*
* DENISE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, version 2.0 of the License only.
*
* DENISE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* Apply time damping (after Brossier (2009))
* last update 31/08/11, D.Koehn
* modified 02/02/12, S.Heider
* ----------------------------------------------------------------------*/
#include "fd.h"
void time_window_glob(float **sectiondata, int iter, int ntr_glob, int ns, int ishot){