diff --git a/bin/Murat_checks.m b/bin/Murat_checks.m
index f98b0d5..ab217e9 100644
--- a/bin/Murat_checks.m
+++ b/bin/Murat_checks.m
@@ -90,8 +90,8 @@
 
 if availableVelocity ==  0
     
-    qLat                =   mean(origin(1),ending(1));
-    qLon                =   mean(origin(2),ending(2));
+    qLat                =   mean([origin(1),ending(1)]);
+    qLon                =   mean([origin(2),ending(2)]);
 
 
     gridPropagation.x   =   xM';
diff --git a/bin/Murat_inversionQ.m b/bin/Murat_inversionQ.m
index 0a36dff..d1d8b6a 100644
--- a/bin/Murat_inversionQ.m
+++ b/bin/Murat_inversionQ.m
@@ -76,6 +76,8 @@
         exitflag        =   'Tikhonov';
         output          =   [];
         saveFigureAsImage(pathFolder);
+        savefig(gcf, [pathFolder '.fig']);
+        close(gcf);
         fval            =   fval*obj0;
     
     case 'Particle'
diff --git a/bin/Murat_inversionQc.m b/bin/Murat_inversionQc.m
index d2a3d52..0797059 100644
--- a/bin/Murat_inversionQc.m
+++ b/bin/Murat_inversionQc.m
@@ -72,7 +72,9 @@
             Murat_tikhonovQc(PlotI,W*Ac_k,W*Qm_k,dampValue,x0.Qc);
         eflag           =   'Tikhonov';
         output          =   [];
-        saveFigureAsImage(pathFolder);        
+        saveFigureAsImage(pathFolder);
+        savefig(gcf, [pathFolder '.fig']);
+        close(gcf);
         fval            =   fval*obj0;
 
     case 'Particle'
diff --git a/bin/Murat_plot.m b/bin/Murat_plot.m
index cbbad2b..b7d3df4 100644
--- a/bin/Murat_plot.m
+++ b/bin/Murat_plot.m
@@ -68,6 +68,20 @@
     rtQk        =   retainQ(:,k);
     rtQck       =   retainQc(:,k);
     rcQck       =   ray_crosses_Qc(:,k);
+    modv_pd_k   =   modv_pd(:,:,k);
+    modv_Qc_k   =   modv_Qc(:,:,k);
+    modv_Q_k    =   modv_Q(:,:,k);
+    [X,Y,Z1,mPD]=   Murat_fold(x,y,z,modv_pd_k(:,4));
+    [~,~,~,mQc] =   Murat_fold(x,y,z,modv_Qc_k(:,4));
+    [~,~,~,mQ]  =   Murat_fold(x,y,z,modv_Q_k(:,4));
+    Z           =   Z1/1000;
+    evst_Qc     =   evst(rtQck,:);
+    evst_pd     =   evst(rtpdk,:);
+    evst_Q      =   evst(rtQk,:);
+    Qm_k        =   Qm(rtQck,k);
+    RZZ_k       =   RZZ(rtQck,k);
+    avQcFreq(1,k)   =   sum(RZZ_k.*Qm_k)/sum(RZZ_k);
+    avQcFreq(2,k)   =   std(Qm_k);
     
     if PlotRays == 1
     % Murat_plot starts plotting the ray distribution if asked by the user.
@@ -77,7 +91,7 @@
     % Peak Delay rays
     FName_peakDelay = ['Rays_PeakDelay_' fcName '_Hz'];
     rma_pd      = rma(:,2:4,rtpdk)/1000;
-    evst_pd     =   evst(rtpdk,:);
+    
     rays_peakDelay  =   Murat_imageRays(rma_pd,origin,ending,evst_pd,...
         x,y,z,FName_peakDelay);
     saveFigureAsImage(makePath(storeFolder, FName_peakDelay));   
@@ -87,7 +101,7 @@
     % The next figure shows the rays for the total attenuation (Q)
     FName_Q     =   ['Rays_Q_' fcName '_Hz'];
     rma_Q       =   rma(:,2:4,rtQk)/1000;
-    evst_Q      =   evst(rtQk,:);
+    
     rays_Q      =   Murat_imageRays(rma_Q,origin,ending,evst_Q,x,y,z,...
         FName_Q);
     saveFigureAsImage(makePath(storeFolder, FName_Q));
@@ -107,14 +121,12 @@
     
     % Qc test
     storeFolder =   fullfile('Tests','Qc');
-    Qm_k        =   Qm(rtQck,k);
-    RZZ_k       =   RZZ(rtQck,k);
+    
     residualQc_k=   residualQc(k);
     luntot_Qc   =   luntot(rtQck)/1000;
     Ac          =   Ac_i(rtQck,rcQck);
 
-    avQcFreq(1,k)   =   sum(RZZ_k.*Qm_k)/sum(RZZ_k);
-    avQcFreq(2,k)   =   std(Qm_k);
+    
 
     Qc_title    =   ['Qc check ' fcName ' Hz'];
     Qc_analysis =   Murat_imageCheckQc(Qm_k,RZZ_k,residualQc_k,...
@@ -168,14 +180,7 @@
     % work with the function "slice". All stored in the sub-folder.
     storeFolder = fullfile('Results','PeakDelay');
 
-    modv_pd_k   =   modv_pd(:,:,k);
-    modv_Qc_k   =   modv_Qc(:,:,k);
-    modv_Q_k    =   modv_Q(:,:,k);
-    [X,Y,Z1,mPD]=   Murat_fold(x,y,z,modv_pd_k(:,4));
-    [~,~,~,mQc] =   Murat_fold(x,y,z,modv_Qc_k(:,4));
-    [~,~,~,mQ]  =   Murat_fold(x,y,z,modv_Q_k(:,4));
-    Z           =   Z1/1000;
-    evst_Qc     =   evst(rtQck,:);
+    
     
     % Peak delays results, using interpolation defined by 'divi'.
     divi            =   5;
diff --git a/bin/saveFigureAsImage.m b/bin/saveFigureAsImage.m
index a58f580..56a878e 100644
--- a/bin/saveFigureAsImage.m
+++ b/bin/saveFigureAsImage.m
@@ -12,6 +12,55 @@
     fig = gcf
 end
 
+% Force white figure background
+fig.Color = [1 1 1];
+
+% Find all axes (including tiled, polar, etc.)
+ax = findall(fig, 'Type', 'axes');
+for k = 1:numel(ax)
+    % Axes background white
+    ax(k).Color = [1 1 1];
+    % Axes lines, tick labels, and title/label text to black
+    ax(k).XColor = [0 0 0];
+    ax(k).YColor = [0 0 0];
+    if isprop(ax(k),'ZColor')       % cartesian 3D axes
+        ax(k).ZColor = [0 0 0];
+    end
+    % Set label and title colors (covers xlabel/ylabel/title objects)
+    try
+        ax(k).Title.Color = [0 0 0];
+    catch
+    end
+    try
+        ax(k).XLabel.Color = [0 0 0];
+        ax(k).YLabel.Color = [0 0 0];
+        if isprop(ax(k),'ZLabel')
+            ax(k).ZLabel.Color = [0 0 0];
+        end
+    catch
+    end
+    % Tick label interpreter may be preserved; ensure color for text objects
+    tickText = findall(ax(k), 'Type', 'text');
+    for t = 1:numel(tickText)
+        tickText(t).Color = [0 0 0];
+    end
+end
+
+% Configure legends
+lg = findall(fig, 'Type', 'legend');
+for L = 1:numel(lg)
+    lg(L).TextColor   = [0 0 0];   % legend text black
+    lg(L).Color       = [1 1 1];   % legend background white
+    lg(L).EdgeColor   = [0 0 0];   % legend border black
+    % If legend contains title
+    try
+        lg(L).Title.Color = [0 0 0];
+    catch
+    end
+end
+
+drawnow; % ensure updates take effect
+
 % Resolve figure handle
 if isempty(fig) || ~ishandle(fig)
     fig = gcf;