Commit 1265a7e8 authored by laura.gassner's avatar laura.gassner

Merge branch 'develop' of git.scc.kit.edu:GPIAG-Software/IFOS2D into develop

parents 5ae8a37e e2632b7a
...@@ -169,7 +169,6 @@ ...@@ -169,7 +169,6 @@
"F_LOW_PASS_END" : "75.0", "F_LOW_PASS_END" : "75.0",
"F_LOW_PASS_INCR" : "10.0", "F_LOW_PASS_INCR" : "10.0",
"ORDER" : "2", "ORDER" : "2",
"ZERO_PHASE" : "0",
"FREQ_FILE" : "frequencies.dat", "FREQ_FILE" : "frequencies.dat",
"WRITE_FILTERED_DATA" : "0", "WRITE_FILTERED_DATA" : "0",
......
...@@ -53,7 +53,7 @@ int main(int argc, char **argv){ ...@@ -53,7 +53,7 @@ int main(int argc, char **argv){
int sum_killed_traces=0, sum_killed_traces_testshots=0, killed_traces=0, killed_traces_testshots=0; 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; 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 = 0.0;
float energy_SH, energy_sum_SH, energy_all_shots_SH, energy_sum_all_shots_SH; float energy_SH, energy_sum_SH, energy_all_shots_SH, energy_sum_all_shots_SH;
float L2_SH, L2sum_SH, L2_all_shots_SH, L2sum_all_shots_SH; float L2_SH, L2sum_SH, L2_all_shots_SH, L2sum_all_shots_SH;
...@@ -157,9 +157,15 @@ int main(int argc, char **argv){ ...@@ -157,9 +157,15 @@ int main(int argc, char **argv){
int nfrq=0; int nfrq=0;
int FREQ_NR=1; int FREQ_NR=1;
float JOINT_EQUAL_PSV=0.0, JOINT_EQUAL_SH=0.0;
float JOINT_EQUAL_PSV_all=0.0, JOINT_EQUAL_SH_all=0.0;
int JOINT_EQUAL_new_max=1;
FILE *fprec, *FPL2; FILE *fprec, *FPL2;
FILE *FPL2_JOINT;
char L2_joint_log[STRING_SIZE];
/* General parameters */ /* General parameters */
int nt_out; int nt_out;
...@@ -1072,11 +1078,19 @@ int main(int argc, char **argv){ ...@@ -1072,11 +1078,19 @@ int main(int argc, char **argv){
C_rho = rho_avg*rho_avg; C_rho = rho_avg*rho_avg;
} }
/* Seperate PSV and SH logging in case of a joint inversion */
if(WAVETYPE==3){
sprintf(L2_joint_log,"%s_JOINT",MISFIT_LOG_FILE);
}
/* Open Log File for L2 norm */ /* Open Log File for L2 norm */
if(FORWARD_ONLY!=1){ if(!FORWARD_ONLY && MYID==0){
if(MYID==0){
if(iter==1){ if(iter==1){
FPL2=fopen(MISFIT_LOG_FILE,"w"); FPL2=fopen(MISFIT_LOG_FILE,"w");
/* Write header for misfit log file */ /* Write header for misfit log file */
if(GRAD_METHOD==1&&VERBOSE) { if(GRAD_METHOD==1&&VERBOSE) {
if (TIME_FILT==0){ if (TIME_FILT==0){
...@@ -1085,14 +1099,20 @@ int main(int argc, char **argv){ ...@@ -1085,14 +1099,20 @@ int main(int argc, char **argv){
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"); 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");
} }
} }
}
if(iter>1){FPL2=fopen(MISFIT_LOG_FILE,"a");} if(WAVETYPE==3) FPL2_JOINT=fopen(L2_joint_log,"w");
} else {
FPL2=fopen(MISFIT_LOG_FILE,"a");
if(WAVETYPE==3) FPL2_JOINT=fopen(L2_joint_log,"a");
} }
} }
/* initialization of L2 calculation */ /* initialization of L2 calculation */
L2=0.0; L2=0.0;
Lcount=0;
energy=0.0; energy=0.0;
L2_all_shots=0.0; L2_all_shots=0.0;
energy_all_shots=0.0; energy_all_shots=0.0;
...@@ -2820,6 +2840,7 @@ int main(int argc, char **argv){ ...@@ -2820,6 +2840,7 @@ int main(int argc, char **argv){
if(MYID==0&&(WAVETYPE==3)) printf("\n SH: L2=%f",L2sum_all_shots_SH/energy_sum_all_shots_SH); if(MYID==0&&(WAVETYPE==3)) printf("\n SH: L2=%f",L2sum_all_shots_SH/energy_sum_all_shots_SH);
} }
sum_killed_traces=0; sum_killed_traces=0;
MPI_Allreduce(&killed_traces,&sum_killed_traces,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD); MPI_Allreduce(&killed_traces,&sum_killed_traces,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
sum_killed_traces_testshots=0; sum_killed_traces_testshots=0;
...@@ -2830,6 +2851,27 @@ int main(int argc, char **argv){ ...@@ -2830,6 +2851,27 @@ int main(int argc, char **argv){
case 2: case 2:
L2t[1]=0.0; L2t[4]=0.0; L2t[1]=0.0; L2t[4]=0.0;
if(JOINT_EQUAL_WEIGHTING){
if(JOINT_EQUAL_new_max){
JOINT_EQUAL_PSV=L2sum/energy_sum;
JOINT_EQUAL_SH=L2sum_SH/energy_sum_SH;
JOINT_EQUAL_PSV_all=L2sum_all_shots/energy_sum_all_shots;
JOINT_EQUAL_SH_all=L2sum_all_shots_SH/energy_sum_all_shots_SH;
JOINT_EQUAL_new_max=0;
}
L2t[1]+=(L2sum/energy_sum)/JOINT_EQUAL_PSV;
L2t[4]+=(L2sum_all_shots/energy_sum_all_shots)/JOINT_EQUAL_PSV_all;
L2t[1]+=(L2sum_SH/energy_sum_SH)/JOINT_EQUAL_SH;
L2t[4]+=(L2sum_all_shots_SH/energy_sum_all_shots_SH)/JOINT_EQUAL_SH_all;
break;
}
if(WAVETYPE==1||WAVETYPE==3){ if(WAVETYPE==1||WAVETYPE==3){
L2t[1]+=L2sum/energy_sum; L2t[1]+=L2sum/energy_sum;
L2t[4]+=L2sum_all_shots/energy_sum_all_shots; L2t[4]+=L2sum_all_shots/energy_sum_all_shots;
...@@ -2839,6 +2881,7 @@ int main(int argc, char **argv){ ...@@ -2839,6 +2881,7 @@ int main(int argc, char **argv){
L2t[1]+=L2sum_SH/energy_sum_SH; L2t[1]+=L2sum_SH/energy_sum_SH;
L2t[4]+=L2sum_all_shots_SH/energy_sum_all_shots_SH; L2t[4]+=L2sum_all_shots_SH/energy_sum_all_shots_SH;
} }
if(MYID==0&&(WAVETYPE==3)) printf("\n Sum: L2=%f",L2t[4]); if(MYID==0&&(WAVETYPE==3)) printf("\n Sum: L2=%f",L2t[4]);
break; break;
...@@ -3117,6 +3160,10 @@ int main(int argc, char **argv){ ...@@ -3117,6 +3160,10 @@ int main(int argc, char **argv){
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); 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);
} }
if(WAVETYPE==3 && MYID==0){
fprintf(FPL2_JOINT,"%d \t %f \t %f\n",iter,L2sum_all_shots/energy_sum_all_shots,L2sum_all_shots_SH/energy_sum_all_shots_SH);
}
/* No update is done here, however model fils are written to disk for easy post processing */ /* No update is done here, however model fils are written to disk for easy post processing */
alpha_SL=0.0; alpha_SL=0.0;
calc_mat_change_test(waveconv_up,waveconv_rho_up,waveconv_u_up,prhonp1,prho,ppinp1,ppi,punp1,pu,iter,1,FORWARD_ONLY,alpha_SL,0,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start); calc_mat_change_test(waveconv_up,waveconv_rho_up,waveconv_u_up,prhonp1,prho,ppinp1,ppi,punp1,pu,iter,1,FORWARD_ONLY,alpha_SL,0,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
...@@ -3158,6 +3205,10 @@ int main(int argc, char **argv){ ...@@ -3158,6 +3205,10 @@ int main(int argc, char **argv){
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); 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);
} }
if(WAVETYPE==3 && MYID==0){
fprintf(FPL2_JOINT,"%d \t %f \t %f\n",iter,L2sum_all_shots/energy_sum_all_shots,L2sum_all_shots_SH/energy_sum_all_shots_SH);
}
/* initiate variables for next iteration */ /* initiate variables for next iteration */
if(use_wolfe_failsafe==1) { if(use_wolfe_failsafe==1) {
L2_hist[iter]=L2_SL_new; L2_hist[iter]=L2_SL_new;
...@@ -3597,6 +3648,15 @@ int main(int argc, char **argv){ ...@@ -3597,6 +3648,15 @@ int main(int argc, char **argv){
case 2: case 2:
L2t[itest]=0.0; L2t[itest]=0.0;
if(JOINT_EQUAL_WEIGHTING){
L2t[itest]+=(L2sum/energy_sum)/JOINT_EQUAL_PSV;
L2t[itest]+=(L2sum_SH/energy_sum_SH)/JOINT_EQUAL_SH;
break;
}
if(WAVETYPE==1||WAVETYPE==3){ if(WAVETYPE==1||WAVETYPE==3){
L2t[itest]+=L2sum/energy_sum; L2t[itest]+=L2sum/energy_sum;
} }
...@@ -3796,6 +3856,9 @@ int main(int argc, char **argv){ ...@@ -3796,6 +3856,9 @@ int main(int argc, char **argv){
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]);} 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{ 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],F_LOW_PASS);} 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);}
if(WAVETYPE==3 && MYID==0){
fprintf(FPL2_JOINT,"%d \t %f \t %f\n",iter,L2sum_all_shots/energy_sum_all_shots,L2sum_all_shots_SH/energy_sum_all_shots_SH);
}
} }
/* saving history of final L2*/ /* saving history of final L2*/
...@@ -3830,6 +3893,9 @@ int main(int argc, char **argv){ ...@@ -3830,6 +3893,9 @@ int main(int argc, char **argv){
if(MYID==0){ if(MYID==0){
fclose(FPL2); fclose(FPL2);
} }
if(WAVETYPE==3 && MYID==0) {
fclose(FPL2_JOINT);
}
} }
if(iter==nfstart){ if(iter==nfstart){
...@@ -3885,7 +3951,7 @@ int main(int argc, char **argv){ ...@@ -3885,7 +3951,7 @@ int main(int argc, char **argv){
} }
/* abort criterion: did not found a step length which decreases the misfit*/ /* abort criterion: did not found a step length which decreases the misfit*/
if((step3==1)&&(TIME_FILT==0&&USE_WORKFLOW==0)){ if((step3==1||wolfe_SLS_failed)&&(TIME_FILT==0&&USE_WORKFLOW==0)){
if(MYID==0){ if(MYID==0){
printf("\n Did not find a step length which decreases the misfit.\n"); printf("\n Did not find a step length which decreases the misfit.\n");
} }
...@@ -3930,6 +3996,7 @@ int main(int argc, char **argv){ ...@@ -3930,6 +3996,7 @@ int main(int argc, char **argv){
wolfe_SLS_failed=0; wolfe_SLS_failed=0;
step3=0; step3=0;
JOINT_EQUAL_new_max=1;
} }
/* ------------------------------------------------- */ /* ------------------------------------------------- */
...@@ -3967,6 +4034,7 @@ int main(int argc, char **argv){ ...@@ -3967,6 +4034,7 @@ int main(int argc, char **argv){
alpha_SL_old=1; alpha_SL_old=1;
step3=0; step3=0;
JOINT_EQUAL_new_max=1;
} }
/* ------------------------------------------------- */ /* ------------------------------------------------- */
...@@ -4004,6 +4072,7 @@ int main(int argc, char **argv){ ...@@ -4004,6 +4072,7 @@ int main(int argc, char **argv){
alpha_SL_old=1; alpha_SL_old=1;
step3=0; step3=0;
JOINT_EQUAL_new_max=1;
} }
} }
......
...@@ -59,7 +59,7 @@ void exchange_par(void){ ...@@ -59,7 +59,7 @@ void exchange_par(void){
extern float npower, k_max_PML; extern float npower, k_max_PML;
extern int INV_STF, N_STF, N_STF_START; extern int INV_STF, N_STF, N_STF_START;
extern char PARA[STRING_SIZE]; extern char PARA[STRING_SIZE];
extern int TIME_FILT, ORDER, ZERO_PHASE,WRITE_FILTERED_DATA; extern int TIME_FILT, ORDER,WRITE_FILTERED_DATA;
extern float F_LOW_PASS_START, F_LOW_PASS_END, F_LOW_PASS_INCR, F_HIGH_PASS; extern float F_LOW_PASS_START, F_LOW_PASS_END, F_LOW_PASS_INCR, F_HIGH_PASS;
extern int LNORM, DTINV; extern int LNORM, DTINV;
extern int STEPMAX; extern int STEPMAX;
...@@ -95,6 +95,7 @@ void exchange_par(void){ ...@@ -95,6 +95,7 @@ void exchange_par(void){
extern int WAVETYPE; extern int WAVETYPE;
extern int SOURCE_SHAPE_SH; extern int SOURCE_SHAPE_SH;
extern int JOINT_INVERSION_PSV_SH_TYPE; extern int JOINT_INVERSION_PSV_SH_TYPE;
extern int JOINT_EQUAL_WEIGHTING;
/* Workflow */ /* Workflow */
extern char FILE_WORKFLOW[STRING_SIZE]; extern char FILE_WORKFLOW[STRING_SIZE];
extern int USE_WORKFLOW; extern int USE_WORKFLOW;
...@@ -325,7 +326,7 @@ void exchange_par(void){ ...@@ -325,7 +326,7 @@ void exchange_par(void){
idum[85] = NO_OF_TESTSHOTS; idum[85] = NO_OF_TESTSHOTS;
idum[86] = ZERO_PHASE; // idum[86] = EMPTY;
idum[87] = VELOCITY; idum[87] = VELOCITY;
...@@ -375,6 +376,8 @@ void exchange_par(void){ ...@@ -375,6 +376,8 @@ void exchange_par(void){
idum[115]=TRKILL_STF_OFFSET; idum[115]=TRKILL_STF_OFFSET;
idum[116]=TRKILL_STF_OFFSET_INVERT; idum[116]=TRKILL_STF_OFFSET_INVERT;
idum[117]=JOINT_EQUAL_WEIGHTING;
} /** if (MYID == 0) **/ } /** if (MYID == 0) **/
MPI_Barrier(MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD);
...@@ -609,7 +612,7 @@ void exchange_par(void){ ...@@ -609,7 +612,7 @@ void exchange_par(void){
NO_OF_TESTSHOTS = idum[85]; NO_OF_TESTSHOTS = idum[85];
ZERO_PHASE = idum[86]; // EMPTY = idum[86];
VELOCITY = idum[87]; VELOCITY = idum[87];
...@@ -660,6 +663,8 @@ void exchange_par(void){ ...@@ -660,6 +663,8 @@ void exchange_par(void){
TRKILL_STF_OFFSET=idum[115]; TRKILL_STF_OFFSET=idum[115];
TRKILL_STF_OFFSET_INVERT=idum[116]; TRKILL_STF_OFFSET_INVERT=idum[116];
JOINT_EQUAL_WEIGHTING=idum[117];
if ( MYID!=0 && L>0 ) { if ( MYID!=0 && L>0 ) {
FL=vector(1,L); FL=vector(1,L);
} }
......
...@@ -72,7 +72,7 @@ int TAPER_STF; ...@@ -72,7 +72,7 @@ int TAPER_STF;
int INV_STF, N_STF, N_STF_START,TW_IND; int INV_STF, N_STF, N_STF_START,TW_IND;
char PARA[STRING_SIZE]; char PARA[STRING_SIZE];
int TIME_FILT, ORDER, ZERO_PHASE; int TIME_FILT, ORDER;
int WRITE_FILTERED_DATA; int WRITE_FILTERED_DATA;
float F_LOW_PASS_START, F_LOW_PASS_END, F_LOW_PASS_INCR, F_HIGH_PASS; float F_LOW_PASS_START, F_LOW_PASS_END, F_LOW_PASS_INCR, F_HIGH_PASS;
...@@ -143,6 +143,8 @@ int VERBOSE; ...@@ -143,6 +143,8 @@ int VERBOSE;
int WAVETYPE; int WAVETYPE;
int JOINT_INVERSION_PSV_SH_TYPE; int JOINT_INVERSION_PSV_SH_TYPE;
int JOINT_EQUAL_WEIGHTING;
float JOINT_INVERSION_PSV_SH_ALPHA_VS; float JOINT_INVERSION_PSV_SH_ALPHA_VS;
float JOINT_INVERSION_PSV_SH_ALPHA_RHO; float JOINT_INVERSION_PSV_SH_ALPHA_RHO;
int SNAPSHOT_START,SNAPSHOT_END,SNAPSHOT_INCR; int SNAPSHOT_START,SNAPSHOT_END,SNAPSHOT_INCR;
\ No newline at end of file
...@@ -28,7 +28,7 @@ char ** varname_list,** value_list; ...@@ -28,7 +28,7 @@ char ** varname_list,** value_list;
void read_par_json(FILE *fp, char *fileinp){ void read_par_json(FILE *fp, char *fileinp){
/* declaration of extern variables */ /* declaration of extern variables */
extern int NX, NY, FDORDER, MAXRELERROR, SOURCE_SHAPE,SOURCE_SHAPE_SH, SOURCE_TYPE, SNAP, SNAP_FORMAT, ACOUSTIC, L, VERBOSE, WAVETYPE,JOINT_INVERSION_PSV_SH_TYPE; extern int NX, NY, FDORDER, MAXRELERROR, SOURCE_SHAPE,SOURCE_SHAPE_SH, SOURCE_TYPE, SNAP, SNAP_FORMAT, ACOUSTIC, L, VERBOSE, WAVETYPE,JOINT_INVERSION_PSV_SH_TYPE,JOINT_EQUAL_WEIGHTING;
extern float DH, TIME, DT, TS, *FL, TAU, VPPML, PLANE_WAVE_DEPTH, PHI, F_REF,JOINT_INVERSION_PSV_SH_ALPHA_VS,JOINT_INVERSION_PSV_SH_ALPHA_RHO; extern float DH, TIME, DT, TS, *FL, TAU, VPPML, PLANE_WAVE_DEPTH, PHI, F_REF,JOINT_INVERSION_PSV_SH_ALPHA_VS,JOINT_INVERSION_PSV_SH_ALPHA_RHO;
extern float XREC1, XREC2, YREC1, YREC2, FPML; extern float XREC1, XREC2, YREC1, YREC2, FPML;
extern float REC_ARRAY_DEPTH, REC_ARRAY_DIST; extern float REC_ARRAY_DEPTH, REC_ARRAY_DIST;
...@@ -62,7 +62,7 @@ void read_par_json(FILE *fp, char *fileinp){ ...@@ -62,7 +62,7 @@ void read_par_json(FILE *fp, char *fileinp){
extern int INV_STF, N_STF, N_STF_START; extern int INV_STF, N_STF, N_STF_START;
extern char PARA[STRING_SIZE]; extern char PARA[STRING_SIZE];
extern int TIME_FILT, ORDER, ZERO_PHASE,WRITE_FILTERED_DATA; extern int TIME_FILT, ORDER,WRITE_FILTERED_DATA;
extern float F_LOW_PASS_START, F_LOW_PASS_END, F_LOW_PASS_INCR, F_HIGH_PASS; extern float F_LOW_PASS_START, F_LOW_PASS_END, F_LOW_PASS_INCR, F_HIGH_PASS;
extern int LNORM, DTINV; extern int LNORM, DTINV;
...@@ -216,6 +216,12 @@ void read_par_json(FILE *fp, char *fileinp){ ...@@ -216,6 +216,12 @@ void read_par_json(FILE *fp, char *fileinp){
fprintf(fp,"For acoustic modelling WAVETYPE is set to %d.\n",WAVETYPE); fprintf(fp,"For acoustic modelling WAVETYPE is set to %d.\n",WAVETYPE);
} }
if(WAVETYPE==3) { if(WAVETYPE==3) {
if (get_int_from_objectlist("JOINT_EQUAL_WEIGHTING",number_readobjects,&JOINT_EQUAL_WEIGHTING,varname_list, value_list)){
JOINT_EQUAL_WEIGHTING=0;
fprintf(fp,"Variable JOINT_EQUAL_WEIGHTING is set to default value %d.\n",JOINT_EQUAL_WEIGHTING);
}
if (get_int_from_objectlist("JOINT_INVERSION_PSV_SH_TYPE",number_readobjects,&JOINT_INVERSION_PSV_SH_TYPE,varname_list, value_list)){ if (get_int_from_objectlist("JOINT_INVERSION_PSV_SH_TYPE",number_readobjects,&JOINT_INVERSION_PSV_SH_TYPE,varname_list, value_list)){
JOINT_INVERSION_PSV_SH_TYPE=1; JOINT_INVERSION_PSV_SH_TYPE=1;
fprintf(fp,"Variable JOINT_INVERSION_PSV_SH_TYPE is set to default value %d.\n",JOINT_INVERSION_PSV_SH_TYPE); fprintf(fp,"Variable JOINT_INVERSION_PSV_SH_TYPE is set to default value %d.\n",JOINT_INVERSION_PSV_SH_TYPE);
...@@ -847,9 +853,6 @@ void read_par_json(FILE *fp, char *fileinp){ ...@@ -847,9 +853,6 @@ void read_par_json(FILE *fp, char *fileinp){
declare_error("Variable ORDER could not be retrieved from the json input file!"); declare_error("Variable ORDER could not be retrieved from the json input file!");
} }
} }
if (get_int_from_objectlist("ZERO_PHASE",number_readobjects,&ZERO_PHASE,varname_list, value_list)){
ZERO_PHASE=0;
fprintf(fp,"Variable ZERO_PHASE is set to default value %i.\n",ZERO_PHASE);}
if (TIME_FILT==2) { if (TIME_FILT==2) {
if (get_float_from_objectlist("F_HIGH_PASS",number_readobjects,&F_HIGH_PASS,varname_list, value_list)){ if (get_float_from_objectlist("F_HIGH_PASS",number_readobjects,&F_HIGH_PASS,varname_list, value_list)){
F_HIGH_PASS=0.0; F_HIGH_PASS=0.0;
...@@ -858,9 +861,6 @@ void read_par_json(FILE *fp, char *fileinp){ ...@@ -858,9 +861,6 @@ void read_par_json(FILE *fp, char *fileinp){
declare_error("Variable FREQ_FILE could not be retrieved from the json input file!"); declare_error("Variable FREQ_FILE could not be retrieved from the json input file!");
if (get_int_from_objectlist("ORDER",number_readobjects,&ORDER,varname_list, value_list)) if (get_int_from_objectlist("ORDER",number_readobjects,&ORDER,varname_list, value_list))
declare_error("Variable ORDER could not be retrieved from the json input file!"); declare_error("Variable ORDER could not be retrieved from the json input file!");
if (get_int_from_objectlist("ZERO_PHASE",number_readobjects,&ZERO_PHASE,varname_list, value_list)){
ZERO_PHASE=0;
fprintf(fp,"Variable ZERO_PHASE is set to default value %i.\n",ZERO_PHASE);}
} }
} }
......
...@@ -39,7 +39,7 @@ void timedomain_filt(float ** data, float fc, int order, int ntr, int ns, int m ...@@ -39,7 +39,7 @@ void timedomain_filt(float ** data, float fc, int order, int ntr, int ns, int m
/* declaration of extern variables */ /* declaration of extern variables */
extern float DT, F_HIGH_PASS; extern float DT, F_HIGH_PASS;
extern int ZERO_PHASE, NT,MYID; extern int NT,MYID;
/* declaration of local variables */ /* declaration of local variables */
int itr, j, ns_reverse; int itr, j, ns_reverse;
...@@ -47,10 +47,8 @@ void timedomain_filt(float ** data, float fc, int order, int ntr, int ns, int m ...@@ -47,10 +47,8 @@ void timedomain_filt(float ** data, float fc, int order, int ntr, int ns, int m
double *seismogram_hp, *seismogram_reverse_hp, T0_hp; double *seismogram_hp, *seismogram_reverse_hp, T0_hp;
seismogram = dvector(1,ns); seismogram = dvector(1,ns);
if (ZERO_PHASE==1) seismogram_reverse = dvector(1,ns);
seismogram_hp = dvector(1,ns); seismogram_hp = dvector(1,ns);
if (ZERO_PHASE==1) seismogram_reverse_hp = dvector(1,ns);
T0=1.0/(double)fc; T0=1.0/(double)fc;
if(F_HIGH_PASS) if(F_HIGH_PASS)
...@@ -66,18 +64,6 @@ void timedomain_filt(float ** data, float fc, int order, int ntr, int ns, int m ...@@ -66,18 +64,6 @@ void timedomain_filt(float ** data, float fc, int order, int ntr, int ns, int m
seife_lpb(seismogram,ns+1,DT,T0,order); /* ns+1 because vector[0] is also allocated and otherwise seife_lpb do not filter the last sample */ seife_lpb(seismogram,ns+1,DT,T0,order); /* ns+1 because vector[0] is also allocated and otherwise seife_lpb do not filter the last sample */
if (ZERO_PHASE==1){
ns_reverse=ns;
for (j=1;j<=ns;j++) {
seismogram_reverse[ns_reverse]=seismogram[j];
ns_reverse--;}
seife_lpb(seismogram_reverse,ns+1,DT,T0,order);
ns_reverse=ns;
for (j=1;j<=ns;j++) {
seismogram[ns_reverse]=seismogram_reverse[j];
ns_reverse--;}
}
for (j=1;j<=ns;j++){ for (j=1;j<=ns;j++){
data[itr][j]=(float)seismogram[j]; data[itr][j]=(float)seismogram[j];
} }
...@@ -92,17 +78,6 @@ void timedomain_filt(float ** data, float fc, int order, int ntr, int ns, int m ...@@ -92,17 +78,6 @@ void timedomain_filt(float ** data, float fc, int order, int ntr, int ns, int m
seife_hpb(seismogram_hp,ns+1,DT,T0_hp,order); seife_hpb(seismogram_hp,ns+1,DT,T0_hp,order);
if (ZERO_PHASE==1){
ns_reverse=ns;
for (j=1;j<=ns;j++) {
seismogram_reverse_hp[ns_reverse]=seismogram_hp[j];
ns_reverse--;}
seife_hpb(seismogram_reverse_hp,ns+1,DT,T0_hp,order);
ns_reverse=ns;
for (j=1;j<=ns;j++) {
seismogram_hp[ns_reverse]=seismogram_reverse_hp[j];
ns_reverse--;}
}
for (j=1;j<=ns;j++){ for (j=1;j<=ns;j++){
data[itr][j]=(float)seismogram_hp[j]; data[itr][j]=(float)seismogram_hp[j];
} }
...@@ -110,9 +85,7 @@ void timedomain_filt(float ** data, float fc, int order, int ntr, int ns, int m ...@@ -110,9 +85,7 @@ void timedomain_filt(float ** data, float fc, int order, int ntr, int ns, int m
} /* end of itr<=ntr loop */ } /* end of itr<=ntr loop */
free_dvector(seismogram,1,ns); free_dvector(seismogram,1,ns);
if (ZERO_PHASE==1) free_dvector(seismogram_reverse,1,ns);
free_dvector(seismogram_hp,1,ns); free_dvector(seismogram_hp,1,ns);
if (ZERO_PHASE==1) free_dvector(seismogram_reverse_hp,1,ns);
} /* end of function */ } /* end of function */
...@@ -39,7 +39,7 @@ void timedomain_filt_vector(float * data, float fc, int order, int ns, int meth ...@@ -39,7 +39,7 @@ void timedomain_filt_vector(float * data, float fc, int order, int ns, int meth
/* declaration of external variables */ /* declaration of external variables */
extern float DT, F_HIGH_PASS; extern float DT, F_HIGH_PASS;
extern int ZERO_PHASE, NT,MYID; extern int NT,MYID;
/* declaration of local variables */ /* declaration of local variables */
int itr, j, ns_reverse; int itr, j, ns_reverse;
...@@ -50,10 +50,8 @@ void timedomain_filt_vector(float * data, float fc, int order, int ns, int meth ...@@ -50,10 +50,8 @@ void timedomain_filt_vector(float * data, float fc, int order, int ns, int meth
declare_error("Order of timedomain filter must be an even number!");} declare_error("Order of timedomain filter must be an even number!");}
seismogram = dvector(1,ns); seismogram = dvector(1,ns);
if (ZERO_PHASE==1) seismogram_reverse = dvector(1,ns);
seismogram_hp = dvector(1,ns); seismogram_hp = dvector(1,ns);
if (ZERO_PHASE==1) seismogram_reverse_hp = dvector(1,ns);
T0=1.0/(double)fc; T0=1.0/(double)fc;
if(F_HIGH_PASS) if(F_HIGH_PASS)
...@@ -68,18 +66,6 @@ void timedomain_filt_vector(float * data, float fc, int order, int ns, int meth ...@@ -68,18 +66,6 @@ void timedomain_filt_vector(float * data, float fc, int order, int ns, int meth
seife_lpb(seismogram,ns+1,DT,T0,order); /* ns+1 because vector[0] is also allocated and otherwise seife_lpb do not filter the last sample */ seife_lpb(seismogram,ns+1,DT,T0,order); /* ns+1 because vector[0] is also allocated and otherwise seife_lpb do not filter the last sample */
if (ZERO_PHASE==1){
ns_reverse=ns;
for (j=1;j<=ns;j++) {
seismogram_reverse[ns_reverse]=seismogram[j];
ns_reverse--;}
seife_lpb(seismogram_reverse,ns+1,DT,T0,order);
ns_reverse=ns;
for (j=1;j<=ns;j++) {
seismogram[ns_reverse]=seismogram_reverse[j];
ns_reverse--;}
}
for (j=1;j<=ns;j++){ for (j=1;j<=ns;j++){
data[j]=(float)seismogram[j]; data[j]=(float)seismogram[j];
} }
...@@ -92,26 +78,13 @@ void timedomain_filt_vector(float * data, float fc, int order, int ns, int meth ...@@ -92,26 +78,13 @@ void timedomain_filt_vector(float * data, float fc, int order, int ns, int meth
seife_hpb(seismogram_hp,ns+1,DT,T0_hp,order); seife_hpb(seismogram_hp,ns+1,DT,T0_hp,order);
if (ZERO_PHASE==1){
ns_reverse=ns;