FEATURE: Added option to write filtered data

Added the option WRITE_FILTERED_DATA to write time filtered measured data to disk.
If WRITE_FILTERED_DATA is not set in the JSON file, it will be set to 0.

src/denise.c

@@ -2075,6 +2075,7 @@ int main(int argc, char **argv){
                                     itestshot+=TESTSHOT_INCR;
                                 }
                                 
+                                /* Write differences between measured and synthetic seismogramms to disk */
                                 if (SEISMO){
                                     if(WAVETYPE==1 || WAVETYPE==3){
                                         catseis(sectionvxdiff, fulldata_vx, recswitch, ntr_glob, MPI_COMM_NTR);
@@ -2089,6 +2090,22 @@ int main(int argc, char **argv){
                                         saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,sectionpdiff,sectionpdiff,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,-1);
                                     }
                                 }
+                                
+                                /* Write measured filtered seismogramms to disk */
+                                if (SEISMO && TIME_FILT && WRITE_FILTERED_DATA){
+                                    if(WAVETYPE==1 || WAVETYPE==3){
+                                        catseis(sectionvxdata, fulldata_vx, recswitch, ntr_glob, MPI_COMM_NTR);
+                                        catseis(sectionvydata, fulldata_vy, recswitch, ntr_glob, MPI_COMM_NTR);
+                                        if(ACOUSTIC)
+                                            catseis(sectionpdata, fulldata_p, recswitch, ntr_glob, MPI_COMM_NTR);
+                                    }
+                                    if(WAVETYPE==2 || WAVETYPE==3){
+                                        catseis(sectionvzdata, fulldata_vz, recswitch, ntr_glob, MPI_COMM_NTR);
+                                    }
+                                    if(myid_ntr==0){
+                                        saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,sectionpdiff,sectionpdiff,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,-2);
+                                    }
+                                }
                             }
                             
                             

src/exchange_par.c

@@ -64,7 +64,7 @@ void exchange_par(void){
     extern float npower, k_max_PML;
     extern int INV_STF, N_STF, N_STF_START;
     extern char PARA[STRING_SIZE];
-    extern int TIME_FILT, ORDER, ZERO_PHASE;
+    extern int TIME_FILT, ORDER, ZERO_PHASE,WRITE_FILTERED_DATA;
     extern float FC_START, FC_END, FC_INCR, F_HP;
     extern int LNORM, DTINV;
     extern int STEPMAX;
@@ -342,6 +342,8 @@ void exchange_par(void){
         idum[108]=WOLFE_NUM_TEST;
         idum[109]=WOLFE_TRY_OLD_STEPLENGTH;
         
+        idum[110]=WRITE_FILTERED_DATA;
+        
     } /** if (MYID == 0) **/
     
     if (MYID != 0) FL=vector(1,L);
@@ -607,6 +609,8 @@ void exchange_par(void){
     
     WOLFE_TRY_OLD_STEPLENGTH=idum[109];
     
+    WRITE_FILTERED_DATA=idum[110];
+    
     MPI_Bcast(&FL[1],L,MPI_FLOAT,0,MPI_COMM_WORLD);
     
 }

src/fd.h

@@ -21,7 +21,7 @@
 #define fsign(x) ((x<0.0)?(-1):1)    
 
 #define PI (3.141592653589793)
-#define NPAR 118
+#define NPAR 119
 #define STRING_SIZE 74
 #define STRING_SIZE2 256
 #define REQUEST_COUNT 4

src/globvar.h

@@ -74,6 +74,7 @@ int INV_STF, N_STF, N_STF_START,TW_IND;
 char PARA[STRING_SIZE];
 
 int TIME_FILT, ORDER, ZERO_PHASE;
+int WRITE_FILTERED_DATA;
 float FC_START, FC_END, FC_INCR, F_HP;
 
 int LNORM, DTINV;

src/read_par_json.c

@@ -64,7 +64,7 @@ void read_par_json(FILE *fp, char *fileinp){
     extern int INV_STF, N_STF, N_STF_START;
     extern char PARA[STRING_SIZE];
     
-    extern int TIME_FILT, ORDER, ZERO_PHASE;
+    extern int TIME_FILT, ORDER, ZERO_PHASE,WRITE_FILTERED_DATA;
     extern float FC_START, FC_END, FC_INCR, F_HP;
     
     extern int LNORM, DTINV;
@@ -764,6 +764,9 @@ void read_par_json(FILE *fp, char *fileinp){
                     TIME_FILT=0;
                     fprintf(fp,"Variable TIME_FILT is set to default value %d.\n",TIME_FILT);}
                 else {
+                    if (get_int_from_objectlist("WRITE_FILTERED_DATA",number_readobjects,&WRITE_FILTERED_DATA,varname_list, value_list)){
+                        WRITE_FILTERED_DATA=0;
+                    }
                     if (get_float_from_objectlist("F_HP",number_readobjects,&F_HP,varname_list, value_list)){
                         F_HP=0.0;
                         fprintf(fp,"Variable F_HP is set to default value %f.\n",F_HP);}

src/saveseis_glob.c

@@ -2,37 +2,35 @@
  * 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>.
------------------------------------------------------------------------------------------*/
+ -----------------------------------------------------------------------------------------*/
 
 /*------------------------------------------------------------------------
- *   write seismograms to files 
+ *   write seismograms to files
  *   last update Oct 2015 F.Wittkamp
  *  ----------------------------------------------------------------------*/
 
 #include "fd.h"
 
-void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectionvz,float **sectionp,
-float **sectioncurl, float **sectiondiv, int  **recpos, int  **recpos_loc, 
-int ntr, float ** srcpos, int ishot, int ns, int iter){ 
-		
-	extern int SEISMO, SEIS_FORMAT, RUN_MULTIPLE_SHOTS, WAVETYPE, VERBOSE;
-	extern char  SEIS_FILE_VX[STRING_SIZE], SEIS_FILE_VY[STRING_SIZE], SEIS_FILE_VZ[STRING_SIZE];
-	extern char  SEIS_FILE_CURL[STRING_SIZE], SEIS_FILE_DIV[STRING_SIZE], SEIS_FILE_P[STRING_SIZE];
-
-        char vxf[STRING_SIZE], vyf[STRING_SIZE],vzf[STRING_SIZE], curlf[STRING_SIZE], divf[STRING_SIZE], pf[STRING_SIZE];
-        int nsrc=1;
+void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectionvz,float **sectionp,float **sectioncurl, float **sectiondiv, int  **recpos, int  **recpos_loc,int ntr, float ** srcpos, int ishot, int ns, int iter){
+    
+    extern int SEISMO, SEIS_FORMAT, RUN_MULTIPLE_SHOTS, WAVETYPE, VERBOSE;
+    extern char  SEIS_FILE_VX[STRING_SIZE], SEIS_FILE_VY[STRING_SIZE], SEIS_FILE_VZ[STRING_SIZE];
+    extern char  SEIS_FILE_CURL[STRING_SIZE], SEIS_FILE_DIV[STRING_SIZE], SEIS_FILE_P[STRING_SIZE];
+    
+    char vxf[STRING_SIZE], vyf[STRING_SIZE],vzf[STRING_SIZE], curlf[STRING_SIZE], divf[STRING_SIZE], pf[STRING_SIZE];
+    int nsrc=1;
     
     if(iter>=0){
         sprintf(vxf,"%s.shot%d.it%d",SEIS_FILE_VX,ishot,iter);
@@ -56,6 +54,16 @@ int ntr, float ** srcpos, int ishot, int ns, int iter){
         sprintf(curlf,"%s.shot%d_adjoint_src",SEIS_FILE_CURL,ishot);
     }
     
+    if(iter==-2){
+        sprintf(vxf,"%s_measured.shot%d.it%d",SEIS_FILE_VX,ishot,iter);
+        sprintf(vyf,"%s_measured.shot%d.it%d",SEIS_FILE_VY,ishot,iter);
+        if(WAVETYPE==2 || WAVETYPE==3) {
+            sprintf(vzf,"%s_measured.shot%d.it%d",SEIS_FILE_VZ,ishot,iter);
+        }
+        sprintf(pf,"%s_measured.shot%d.it%d",SEIS_FILE_P,ishot,iter);
+        sprintf(divf,"%s_measured.shot%d.it%d",SEIS_FILE_DIV,ishot,iter);
+        sprintf(curlf,"%s_measured.shot%d.it%d",SEIS_FILE_CURL,ishot,iter);
+    }
     
     switch (SEISMO){
         case 1 : /* particle velocities only */
@@ -118,5 +126,5 @@ int ntr, float ** srcpos, int ishot, int ns, int iter){
             }
             break;
             
-    }     
+    }
 }