Adjusted L2 weighting of SH and PSV waves

Both L2-misfits are now normalized at each iteration. Therefore, both objective functions should be weighted equally.

src/IFOS2D.c

@@ -157,6 +157,9 @@ int main(int argc, char **argv){
     int nfrq=0;
     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;
     
@@ -2838,6 +2841,27 @@ int main(int argc, char **argv){
                             case 2:
                             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){
                                 L2t[1]+=L2sum/energy_sum;
                                 L2t[4]+=L2sum_all_shots/energy_sum_all_shots;
@@ -2847,6 +2871,7 @@ int main(int argc, char **argv){
                                 L2t[1]+=L2sum_SH/energy_sum_SH;
                                 L2t[4]+=L2sum_all_shots_SH/energy_sum_all_shots_SH;
                             }
+                            
                         if(MYID==0&&(WAVETYPE==3)) printf("\n Sum: L2=%f",L2t[4]);
                             
                             break;
@@ -3613,6 +3638,15 @@ int main(int argc, char **argv){
                         case 2:
                             L2t[itest]=0.0;
                             
+                            if(JOINT_EQUAL_WEIGHTING){
+                                
+                                L2t[1]+=(L2sum/energy_sum)/JOINT_EQUAL_PSV;
+                                
+                                L2t[1]+=(L2sum_SH/energy_sum_SH)/JOINT_EQUAL_SH;
+                                
+                                break;
+                            }
+                            
                             if(WAVETYPE==1||WAVETYPE==3){
                                 L2t[itest]+=L2sum/energy_sum;
                             }
@@ -3950,6 +3984,8 @@ int main(int argc, char **argv){
                 PCG_iter_start=iter+1;
                 
                 wolfe_SLS_failed=0;
+                
+                JOINT_EQUAL_new_max=1;
             }
             
             /* ------------------------------------------------- */
@@ -3985,6 +4021,8 @@ int main(int argc, char **argv){
                 
                 wolfe_SLS_failed=0;
                 alpha_SL_old=1;
+                
+                JOINT_EQUAL_new_max=1;
             }
             
             /* ------------------------------------------------- */
@@ -4020,6 +4058,8 @@ int main(int argc, char **argv){
                 
                 wolfe_SLS_failed=0;
                 alpha_SL_old=1;
+                
+                JOINT_EQUAL_new_max=1;
             }
             
         }

src/exchange_par.c

@@ -95,6 +95,7 @@ void exchange_par(void){
     extern int WAVETYPE;
     extern int SOURCE_SHAPE_SH;
     extern int JOINT_INVERSION_PSV_SH_TYPE;
+    extern int JOINT_EQUAL_WEIGHTING;
     /* Workflow  */
     extern char FILE_WORKFLOW[STRING_SIZE];
     extern int USE_WORKFLOW;
@@ -375,6 +376,8 @@ void exchange_par(void){
         idum[115]=TRKILL_STF_OFFSET;
         idum[116]=TRKILL_STF_OFFSET_INVERT;
         
+        idum[117]=JOINT_EQUAL_WEIGHTING;
+        
     } /** if (MYID == 0) **/
     
     MPI_Barrier(MPI_COMM_WORLD);
@@ -660,6 +663,8 @@ void exchange_par(void){
     TRKILL_STF_OFFSET=idum[115];
     TRKILL_STF_OFFSET_INVERT=idum[116];
     
+    JOINT_EQUAL_WEIGHTING=idum[117];
+    
     if ( MYID!=0 && L>0 ) {
         FL=vector(1,L);
     }

src/globvar.h

@@ -143,6 +143,8 @@ int VERBOSE;
 int WAVETYPE;
 int JOINT_INVERSION_PSV_SH_TYPE;
 
+int JOINT_EQUAL_WEIGHTING;
+
 float JOINT_INVERSION_PSV_SH_ALPHA_VS;
 float JOINT_INVERSION_PSV_SH_ALPHA_RHO;
 int SNAPSHOT_START,SNAPSHOT_END,SNAPSHOT_INCR;
\ No newline at end of file

src/read_par_json.c

@@ -28,7 +28,7 @@ char ** varname_list,** value_list;
 void read_par_json(FILE *fp, char *fileinp){
     
     /* 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 XREC1, XREC2, YREC1, YREC2, FPML;
     extern float REC_ARRAY_DEPTH, REC_ARRAY_DIST;
@@ -216,6 +216,12 @@ void read_par_json(FILE *fp, char *fileinp){
                 fprintf(fp,"For acoustic modelling WAVETYPE is set to %d.\n",WAVETYPE);
             }
             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)){
                     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);