Cleaning

bin/Murat_body.m

@@ -20,18 +20,22 @@
 %    energyRatioBodyCoda_i:     energy ratio between body and coda waves
 %    energyRatioCodaNoise_i:    energy ratio between coda waves and noise
 
+% Redefine windows with sampling for noise and direct waves
 int                     =   bodyWindow*srate_i;
 lc                      =   length(cursorCodaStart_i:cursorCodaEnd_i);
 cursor0                 =	floor(startNoise*srate_i);
 cursor0_1               =   floor(cursor0 + int);
 cursor2                 =   floor(cursorPick_i + int-1);
 
+% Direct wave normalised by window length
 spamp                   =...
     trapz(sp_i(cursorPick_i:cursor2,:))/int;
 
+% Noise pre-pick normalised by window length
 spampn                  =...
     trapz(sp_i(cursor0:cursor0_1,:))/int;
 
+% Coda normalised by window length
 spampc                  =...
     trapz(sp_i(cursorCodaStart_i:cursorCodaEnd_i,:))/lc;
 

bin/Murat_checks.m

@@ -1,9 +1,7 @@
 % ADDITIONAL input variables that are not set by the user.
 function Murat                  =   Murat_checks(Murat)
 
-% INPUTS
-%dataDirectory                   =   ['./' Murat.input.dataDirectory];
-dataDirectory                   =   [Murat.input.dataDirectory]; %if data is outside of current folder
+dataDirectory                   =   [Murat.input.dataDirectory];
 FPath                           =   './';
 FLabel                          =   Murat.input.label;
 PTime                           =   ['SAChdr.times.' Murat.input.PTime];

bin/Murat_codaCheck.m

@@ -25,31 +25,37 @@
 %    cursorCodaStart_i:     coda starting time after check on trace
 %    cursorCodaEnd_i:       coda end time after check on trace
 
+%Define envelope duration
 t00                                             =   tempis(1);
 lengthTempis                                    =   length(tempis);
 
+% Method peak is generally valid for active seismics
 if isequal(peakDelayMethod,'Peak')
 
     tCoda_i                                     =...
         (pktime_i-originTime_i)+peakDelay_i;
-
+
+% Method constant is the most used in small local tomography 
 elseif isequal(peakDelayMethod,'Constant')
 
     tCoda_i                                     =   tCm;
-
+
+% Method travel is the standard at the regional scale
 elseif isequal(peakDelayMethod,'Travel')
 
     tCoda_i                                     =...
         originTime_i+ tCm*(pktime_i-originTime_i);
 
 end
 
+% Define the indexes along the seismogram
 cursorCodaStart_i                               =...
     floor((originTime_i - t00 + tCoda_i) * srate_i - 1);
 
 cursorCodaEnd_i                                 =...
     floor(cursorCodaStart_i + tWm * srate_i - 1);
 
+% Some seismograms are not long enough so you consider a shorter window
 if cursorCodaEnd_i > lengthTempis
     cursorCodaEnd_i                             =   lengthTempis;
 end

bin/Murat_cornering.m

@@ -14,8 +14,8 @@
 % Output parameters:
 %    ip:        index of minimum x grid coord - West;
 %    jp:        index of minimum y grid coord - South;
-%    kp:     i  ndex of maximum z grid coord - shallowest.
-
+%    kp:        index of maximum z grid coord - shallowest.
+%    flag:      flags if ray outside of grid
 
 % Coordinates of the propagation grid are unwrapped
 xGrid                       =   gridD.x;

bin/Murat_data.m

@@ -21,6 +21,7 @@
 cf                                      =   Murat.input.centralFrequency;
 lcf                                     =   length(cf);
 
+workers                                 =   Murat.input.workers;
 origin                                  =   Murat.input.origin;
 originTime                              =   Murat.input.originTime;
 PTime                                   =   Murat.input.PTime;
@@ -68,7 +69,7 @@
 % on peak delays and coda attenuation.
 count_trash                             =   0;
 
-for i = 1:lengthData
+parfor (i = 1:lengthData,workers)
 
     if isequal(mod(i,100),0)
 
@@ -182,21 +183,21 @@
 end
 
 % Setting up the final data vectors and matrices with checks on values
-Murat.data.locationsDeg                 =   locationDeg;
-Murat.data.locationsM                   =   locationM;
-Murat.data.theoreticalTime              =   theoreticalTime;
-Murat.data.peakDelay                    =   peakDelay;
-Murat.data.inversionMatrixPeakDelay     =   inversionMatrixPeakDelay;
-Murat.data.inversionMatrixQ             =   inversionMatrixQ;
-Murat.data.totalLengthRay               =   totalLengthRay;
-Murat.data.raysPlot                     =   raysPlot;
-Murat.data.rayCrossing                  =   sum(rayCrossing);
-Murat.data.inverseQc                    =   inverseQc; 
-Murat.data.uncertaintyQc                =   uncertaintyQc; 
-Murat.data.inversionMatrixQc            =   inversionMatrixQc;
-Murat.data.energyRatioBodyCoda          =   energyRatioBodyCoda; 
-Murat.data.energyRatioCodaNoise         =   energyRatioCodaNoise;
-Murat.data.tCoda                        =   tCoda;
+Murat.rays.locationsDeg                 =   locationDeg;
+Murat.rays.locationsM                   =   locationM;
+Murat.rays.theoreticalTime              =   theoreticalTime;
+Murat.PD.peakDelay                      =   peakDelay;
+Murat.PD.inversionMatrixPeakDelay       =   inversionMatrixPeakDelay;
+Murat.Q.inversionMatrixQ                =   inversionMatrixQ;
+Murat.rays.totalLengthRay               =   totalLengthRay;
+Murat.rays.raysPlot                     =   raysPlot;
+Murat.rays.rayCrossing                  =   sum(rayCrossing);
+Murat.Qc.inverseQc                      =   inverseQc; 
+Murat.Qc.uncertaintyQc                  =   uncertaintyQc; 
+Murat.Qc.inversionMatrixQc              =   inversionMatrixQc;
+Murat.Q.energyRatioBodyCoda             =   energyRatioBodyCoda; 
+Murat.Q.energyRatioCodaNoise            =   energyRatioCodaNoise;
+Murat.Qc.tCoda                          =   tCoda;
 
 Murat                                   =   Murat_selection(Murat);
 

bin/Murat_dataWarning.m

@@ -19,6 +19,7 @@
 %    rapspcn:           coda to noise ratio
 %    comp:           	components
 %    flag:           	flag to change between one or more components
+%    QcM:           	chosen method to measure Qc between Lin e NonLin
 %
 % Output parameters:
 %    problempd:         stores the data not used because of peak delays

bin/Murat_declustering.m

@@ -1,21 +1,24 @@
 function Murat = Murat_declustering(Murat,factor)
 % function Murat = Murat_declustering(Murat,factor)
 %
-% PLOTS rays before and after declustering
+% SAVES declustered variables and PLOTS rays before and after declustering.
+% The declustering is based on dividing the original grid into finer node
+% spacing and selecting those with the best coda decay. 
 %
 % Input parameters:
 %    Murat:                Murat structure variable
-%    factor:               factor by which original grid is divided into
+%    factor:               input factor used to divide the original grid
+
 components                  =   Murat.input.components;
 origin                      =   Murat.input.origin;
 ending                      =   Murat.input.end;
 x                           =   Murat.input.x;
 y                           =   Murat.input.y;
 z                           =   Murat.input.z;
 
-locationsDeg                =   Murat.data.locationsDeg;
-uncertaintyQc               =   Murat.data.uncertaintyQc;
-Murat.data.locDegOriginal   =   locationsDeg;
+locationsDeg                =   Murat.rays.locationsDeg;
+uncertaintyQc               =   Murat.Qc.uncertaintyQc;
+Murat.rays.locDegOriginal   =   locationsDeg;
 
 disp(['used waveforms before declustering: ',...
     num2str(length(locationsDeg)*components)])
@@ -67,28 +70,26 @@
 
 %% remove deleted data from Murat structure variable
 ind_to_keep                 =   new_evestaz(:,8);
-Murat.data.energyRatioBodyCoda  =...
-    Murat.data.energyRatioBodyCoda(ind_to_keep,:);
-Murat.data.energyRatioCodaNoise =...
-    Murat.data.energyRatioCodaNoise(ind_to_keep,:);
-Murat.data.inverseQc        =   Murat.data.inverseQc(ind_to_keep,:);
-Murat.data.inversionMatrixPeakDelay	=...
-    Murat.data.inversionMatrixPeakDelay(ind_to_keep,:);
-Murat.data.inversionMatrixQ =...
-    Murat.data.inversionMatrixQ(ind_to_keep,:);
-Murat.data.inversionMatrixQc=...
-    Murat.data.inversionMatrixQc(ind_to_keep,:);
-Murat.data.locationsDeg     =   Murat.data.locationsDeg(ind_to_keep,:);
-Murat.data.peakDelay        =   Murat.data.peakDelay(ind_to_keep,:);
-Murat.data.retainPeakDelay  =   Murat.data.retainPeakDelay(ind_to_keep,:);
-Murat.data.retainQ          =   Murat.data.retainQ(ind_to_keep,:);
-Murat.data.retainQc         =   Murat.data.retainQc(ind_to_keep,:);
-Murat.data.tCoda            =   Murat.data.tCoda(ind_to_keep,:);
-Murat.data.totalLengthRay   =   Murat.data.totalLengthRay(ind_to_keep,:);
-Murat.data.travelTime       =   Murat.data.travelTime(ind_to_keep,:);
-Murat.data.uncertaintyQc    =   Murat.data.uncertaintyQc(ind_to_keep,:);
-Murat.data.variationPeakDelay =...
-    Murat.data.variationPeakDelay(ind_to_keep,:);
+Murat.Q.energyRatioBodyCoda =...
+                            Murat.Q.energyRatioBodyCoda(ind_to_keep,:);
+Murat.Q.energyRatioCodaNoise=...
+    Murat.Q.energyRatioCodaNoise(ind_to_keep,:);
+Murat.Qc.inverseQc          =   Murat.Qc.inverseQc(ind_to_keep,:);
+Murat.PD.inversionMatrixPeakDelay...
+                            =...
+    Murat.PD.inversionMatrixPeakDelay(ind_to_keep,:);
+Murat.Q.inversionMatrixQ    =   Murat.Q.inversionMatrixQ(ind_to_keep,:);
+Murat.Qc.inversionMatrixQc  =   Murat.Qc.inversionMatrixQc(ind_to_keep,:);
+Murat.rays.locationsDeg     =   Murat.rays.locationsDeg(ind_to_keep,:);
+Murat.PD.peakDelay          =   Murat.PD.peakDelay(ind_to_keep,:);
+Murat.PD.retainPeakDelay    =   Murat.PD.retainPeakDelay(ind_to_keep,:);
+Murat.Q.retainQ             =   Murat.Q.retainQ(ind_to_keep,:);
+Murat.Qc.retainQc           =   Murat.Qc.retainQc(ind_to_keep,:);
+Murat.Qc.tCoda              =   Murat.Qc.tCoda(ind_to_keep,:);
+Murat.rays.totalLengthRay   =   Murat.rays.totalLengthRay(ind_to_keep,:);
+Murat.rays.travelTime       =   Murat.rays.travelTime(ind_to_keep,:);
+Murat.Qc.uncertaintyQc      =   Murat.Qc.uncertaintyQc(ind_to_keep,:);
+Murat.PD.variationPeakDelay =   Murat.PD.variationPeakDelay(ind_to_keep,:);
 
 disp(['used waveforms after declustering: ',...
     num2str(length(ind_to_keep)*components)])

bin/Murat_image3D.m

@@ -1,7 +1,7 @@
 function [image,mVp,xp,yp,zp,Xp,Yp,Zp]                    = ...
     Murat_image3D(X,Y,Z,V,color,sections,evestaz,x,y,z,divi,name)
 % function image          = ...
-%     Murat_image3D(X,Y,Z,V,color,sections,evestaz,x,y,z,name)
+%     Murat_image3D(X,Y,Z,V,color,sections,evestaz,x,y,z,divi,name)
 %
 % PLOTS a 3D image of a field on slices.
 %
@@ -50,13 +50,13 @@
 if max_scale            <   1
     max_scale           =   round(max_scale,2);
     if max_scale        ==  0
-        caxis([-1 1])
+        clim([-1 1])
     else
-        caxis([-max_scale max_scale])
+        clim([-max_scale max_scale])
     end
 elseif max_scale        >   5
     max_scale           =   ceil(max_scale);
-    caxis([-max_scale max_scale])
+    clim([-max_scale max_scale])
 end
 
 colormap(color);

bin/Murat_image3D_2panels.m

@@ -1,7 +1,7 @@
 function image          = ...
     Murat_image3D_2panels(X,Y,Z,V,color,sections,evestaz,x,y,z)
 % function image          = ...
-%     Murat_image3D(X,Y,Z,V,color,sections,evestaz,x,y,z,name)
+%     Murat_image3D(X,Y,Z,V,color,sections,evestaz,x,y,z)
 %
 % PLOTS a 3D image of a field on slices on two subpanels.
 %
@@ -16,7 +16,6 @@
 %    x:         x vector
 %    y:         y vector
 %    z:         z vector
-%    name:      name of the title of the figure
 %
 % Output parameters:
 %    image:     image produced in one panel

bin/Murat_image3Dparameters.m

@@ -46,10 +46,10 @@
 max_scale               =   max(scale_mVp);
 if max_scale            <   1
     max_scale           =   round(max_scale,2);
-    caxis([-max_scale max_scale])
+    clim([-max_scale max_scale])
 elseif max_scale        >   5
     max_scale           =   ceil(max_scale);
-    caxis([-max_scale max_scale])
+    clim([-max_scale max_scale])
 end
 
 colormap(color);

bin/Murat_imageCheckQc.m

@@ -12,7 +12,7 @@
 %    Ac:            inverse coda attenuation matrix
 %    sizeTitle:     size of title font
 %    Qc_title:      title of Qc figure
-%    QcM:           Linearized or Non Linear measurement
+%    QcM:           Linearized or Non Linear Qc measurement
 %
 % Output parameters:
 %    Qc_analysis:   figure for Qc check

bin/Murat_imageDeclustering.m

@@ -1,7 +1,7 @@
 function clustering       =   Murat_imageDeclustering(locDegOriginal,...
     locDegrees,origin,ending,name)
 % function clustering       =...
-%     Murat_imageDeclustering(locDegOriginal,locDegrees,origin,ending)
+%     Murat_imageDeclustering(locDegOriginal,locDegrees,origin,ending,name)
 %
 % PLOTS 2D rays before and after declustering.
 %

bin/Murat_imageParameters.m

@@ -10,6 +10,7 @@
 %    z:                 z vector
 %    modv_pd_k:         results of peak delay imaging
 %    modv_Qc_k:         results of Qc imaging
+%    sizeTitle:         size of the title font
 %
 % Output parameters:
 %    image:             image produced

bin/Murat_imageParametersMaps.m

@@ -5,18 +5,24 @@
 %     Murat_imageParametersMaps(par,para_map,x,y,z,X,Y,Z,evestaz_Qc,...
 %     sections,sizeTitle,FName_Parameters)
 %
-% PLOTS the parameter models in their space and on a 3D map
+% PLOTS the parameter models in a 3D map
 %
 % Input parameters:
+%    par:               condition on the parameters to image
+%    para_map:          parameter map
 %    x:                 x vector
 %    y:                 y vector
 %    z:                 z vector
-%    modv_pd_k:         results of peak delay imaging
-%    modv_Qc_k:         results of Qc imaging
+%    X:                 X matrix
+%    y:                 Y matrix
+%    z:                 Z matrix
+%    evestaz_Qc:        events and stations for Qc
+%    sections:          input sections
+%    sizeTitle:         font size for title
+%    FName_Parameters:  title of the image
 %
 % Output parameters:
-%    image:             image produced
-%    para_condition:    condition on the parameters to image
+%    ParaMap:           parameter map 3D image
 %    para_map:          parameter map
 
 for n = 1:length(par(:,1))