%Plotting routine to visialize Griggs Rig modelling results 
%working for Simple shear, pure shear, inclusions in PS
%based on: Gerya, 2009. Introduction to numerical geodynamic modelling
%this script is available at Zenodo: https://doi.org/10.5281/zenodo.6951760
%input file is generated with: calc_Griggs_v1_0.m
%by Sebastian Cionoiu 
% --- readme ---
% lines 12 - 32 are used to define general plotting sttings
% later in script material-dependent plotting is defined


clear all;

scPlotBigFig =  [12];%set to 0 for plotting all figs, some value/vector (including 1) for big fig : [1 17]
savefigs = 1;
savecycleOV = 0;%12;%11;%15%overwirte for partial plots - print if overwriting
plotsingle = 2; %set to 0 = all parameters + final stress-strain curve, 1 = stress-time + stress-strain curve,  
                %2 = only final stress-strain curve, 3 = only width curve;
                %3 = only width plot; 4 = only histograms
                %5 = EXPORT dSIG for grain-size (only last in list scPlotBigFig)

calcintv = 1; %for 1 inbetween step: 8 on new runs (97/160 MB), 30 on old runs (300 MB) 
plotSteps = 1; %inbetween plot steps (not first and last) - usually 1 or 3

globalPscale = 0;   % scale for zoom plots (P and sig) 1: all plots equal, 0 individual scales

Ptrans_scale = 0;  % scale  bar for single plots: adapt around phase transition
% Ptrans defined after load!
P_shift = 180e6; %confining pressure modification (does not affect stress and strain) 511: 180 | 528: 60

ModName = '528SC';%part of output filename
ifiles = dir('*220729_184818.mat');%%file to be read

for ifile = ifiles'
    file = ifile.name%print filename!
    clearvars -except file i ifile plotloop savefigs savecycleOV plotsingle makevid...
        plotintv calcintv scPlotBigFig pCscf globalPscale Ptrans Ptrans_scale ModName plotSteps P_shift% clear vars to be sure to not inherit
    close all
    
    if (savefigs==1 && scPlotBigFig(1)==0 && plotsingle==0)
        h1 = figure(1); clf;
        %set(h1,'visible','off');
        set(gcf, 'Position', get(0,'Screensize')) 
        drawnow
        h2 = figure(2); clf;
        set(gcf, 'Position', get(0,'Screensize'))
        %set(h2,'visible','off');
        drawnow
        figure(3); clf;
        set(gcf, 'Position', get(0,'Screensize'))
        drawnow
        if plotsingle == 1
            figure(4),clf
            set(gcf, 'Position', get(0,'Screensize'))
            drawnow
        end
        figure(5); clf;
        set(gcf, 'Position', get(0,'Screensize'))
        drawnow
        figure(7); clf;
        set(gcf, 'Position', get(0,'Screensize'))
        drawnow
        figure(8); clf;
        set(gcf, 'Position', get(0,'Screensize'))
        drawnow
    end
    
%colormap for material plot ; https://de.mathworks.com/matlabcentral/answers/265914-how-to-make-my-own-custom-colormap
MatCmap =  [216 233 233;...%salt
            034 91 170;...%sample
            092 203 034;...%platinum
            211 188 095] / 255;%piston

%file = sprintf('v8_PS_qz_n10_2e7.mat');
load(file) %_1025_elli_500
file=file(1:end-4);%erase(file,'.mat');%for plot title %=file;%
filestr=strrep(file(1:end-16),'_','\_');%for plot title %

%correct confinign P if needed
Pconf       = Pconf + P_shift;%for consistency!!!
PrPL        = PrPL + P_shift;

if savefigs ==1 %create and enter subfolder for each plot file
    mkdir(file); cd(file)
end

FIGps=[100, 100, 1800, 1400];

if savecycleOV ~= 0%overwrite savecyclcle for partial time plots
    savecycle = savecycleOV;
end

plotintv = (savecycle-1)/(plotSteps+1); %for 1 inbetween step: 8 on new runs (97/160 MB), 30 on old runs (300 MB)
if ~mod(plotintv,1) == 0; warning('plot steps are not integer. savacycle and intv and will be adapted'); end %test if plotintv is integer

%disregard plotinv for only final plots
if scPlotBigFig(1) == 0 && plotsingle == 0 || plotsingle == 4
    plotintv = floor(plotintv)
    savecycle = plotintv * (plotSteps+1) +1
else
    plotintv = savecycle - 1
    savecycle
end


%scaling factor for contour plots - for calcite (n>5 usually lower streses)
if nvv > 5 %calcite
    pCscf = 0.2;%Cal: 0.1 - 0.2
else % quartz
    pCscf = 0.6;%Qz: 0.5-1
end

%scaling for phase-transition plots
if nvv > 5          %Qz750°C: 2.75 GPa, Cal600°C:1.52GPa
    Ptrans = 1.52e9; %Cal
    raPTpm = 0.2; %extended range of contours around colorbar 
else
    Ptrans = 2.75e9;% Qz
    raPTpm = 0.2;
end

%for EULER use - in order to have colorbar in working directory
%copyfile ./../plotreq/bluewhitered.m bluewhitered.m

%Bugfix for Version 13-14 (time-vectors were allocated as Ny*Nx*sc matrix)
exptimePL   = exptimePL(1:savecycle);
elapsedTimePL = elapsedTimePL(1:savecycle);

%scaling factors for invariant to differential stress or strain
if SZang == 0
    gamSF = 2;
    sigSF = 2;
else
    gamSF = 2;
    sigSF = 2;
end

PSscale     = 1e6;
xyscale     = 1000;
EXPsig1     = zeros(1,savecycle);
EXPsig1_acf = zeros(1,savecycle);
EXPdevsig1  = zeros(1,savecycle);

uPtYV       = zeros(1,savecycle);
gSZth       = zeros(1,savecycle);
gSZlength   = zeros(1,savecycle);
gSZds_inc   = zeros(1,savecycle);
ginc        = zeros(1,savecycle);

SIG1mag     = zeros(Ny,Nx);
SIG1mag_A   = zeros(Ny,Nx,savecycle);
SIG3mag     = zeros(Ny,Nx);
SIG3mag_A   = zeros(Ny,Nx,savecycle);

ModeMeanEiiSZ   = zeros(savecycle,2);
ModeMeanETASZ   = zeros(savecycle,2);

DSIG            = zeros(Ny,Nx,2,2);
SIG1deg         = zeros(Ny,Nx);
SIG3deg         = zeros(Ny,Nx);
SIGnTAU         = zeros(Ny,Nx,2,2);
Pminmax         = zeros(savecycle,4);
SIG1magminmax   = zeros(savecycle,2);
SIG3magminmax   = zeros(savecycle,2);
SIG1degNminmax  = zeros(savecycle,2);
SIGiiminmax    = zeros(savecycle,4);
EiiSUMminmax     = zeros(savecycle,2);
Eiiminmax       = zeros(savecycle,2);
ETAminmax       = zeros(savecycle,2);
PdSIG1minmax    = zeros(savecycle,2);

if SZang ~= 0 %defining size of readout point array
    probelength = 5;
else 
    probelength = 6;
end
probe_CxA       = zeros(savecycle,probelength);
probe_CyA       = zeros(savecycle,probelength);
probe_TxA       = zeros(savecycle,4);
probe_TyA       = zeros(savecycle,4);
P_probe_Cp      = zeros(savecycle,probelength);
P_probe_conf    = zeros(1,savecycle);
Sig1_probe_Cp   = zeros(savecycle,probelength);
Sig3_probe_Cp   = zeros(savecycle,probelength);

scplots         = [1:1:savecycle]';%counter for number of calc-cycles
%EDIT1b
%add checking of max number of CAL+Elli nodes to define storage vector length
nRHO_rd = round(RhoPL(:,:,savecycle)*100)/100;
nSZ     =   round(sum(sum(nRHO_rd == RHO_sz))*1.2);%number of points in deforming zone (RHO=2) plus safety (e.g. ellips)
if SZang == 0%fraction of Max-Min values to evaluate for plotting points
    nMM     =   max(10,round(nSZ/1.2*0.01));%in pure shear
else
    nMM     =   max(10,round(nSZ/1.2*0.05));%in simple shear
end
%fraction of Max-Min values to evaluate for range plots mode
nMM_range   = 0.03;
EiiSZ   =   zeros(savecycle,nSZ);
ETASZ   =   zeros(savecycle,nSZ);

dx      =   xsize/(Nx-1); %grid step, always one step less than points
dy      =   ysize/(Ny-1);
x_m     =   0:dx:xsize; %horizontal coordinates of nodes in meters
y_m     =   0:dy:ysize; %vertical coordinates of nodes in meters
x       =   x_m*xyscale;
y       =   y_m*xyscale;
qs      =   fix(Nx/30)*4;%quiver stepsize 2 at low res; 6 at middle (100x150) 10 at higher;
% prepare grids for interpolation of Vx/Vy for Quiver
[msX,msY] = meshgrid(x,y);
[msXq,msYq] = meshgrid(x(1:qs:end), y(1:qs:end));

%scaling for zoom-plots - after re-definition of x    
scX=max(Nx)/max(x);
scY=max(Ny)/max(y);

set(0,'defaultTextInterpreter','latex')%for typesetting all titles properly
fontname = 'Arial'; 
set(0,'DefaultAxesFontName',fontname,'DefaultTextFontName',fontname);


%% readout, calc and plot cycle  
maxminEND = 1;%set switch to enter loop to read out final max-min values
 
%differentiate  beteween loop for plotting all figs in regular steps and single fig for defined cycle
if scPlotBigFig == 0
    scVect= [0, 1:calcintv:savecycle];%iteration 0 is replaced by last run to calculate max-min values for plots
else
    scVect= [0,scPlotBigFig];%iteration 0 needed 
end
    
for sc = scVect
    
    %EDIT7a
    %add readout of max - min Values for plots
    %before first full calculation and plot run, compute values of final
    %plots and store the max and min values
    if maxminEND == 1
        sc = max(savecycle);  
    end
    if sc > savecycle
        cd ..%back to root folder when looping through files
        error('manually set variable "sc" grater than number of savecycle - correct!')
    end
   
    %load single plot variables 
    RHO         = RhoPL(:,:,sc); RHO=round(RHO*100)/100;
    P           = PrPL(:,:,sc);
    ETA         = EtaPL(:,:,sc);
    EiiSUM      = EiiSUMPL(:,:,sc);           
    Vx          = VxPL(:,:,sc);
    Vy          = VyPL(:,:,sc);
    Wdot        = WdotPL(:,:,sc);
    WSUM        = WSUMPL(:,:,sc);
    %load and transform corner markers
    MCx         = MCxPL(:,sc);
    MCy         = MCyPL(:,sc);
    MCNx        = round(MCx*Nx/xsize)+[1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 0 0]';%move markers closer to interface
    MCNy        = round(MCy*Ny/ysize)+[1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 0 0 -1 1]';

    %% Calculate Strain and stress
    %all calculated to pressure pints (center of grid=(Nodal points = Np)-crosses )
    Exx = zeros(Ny,Nx);
    Exx(:,2:end) = diff(Vx(:,:),1,2)/dx;%Exx in P points
    ExxNp = (Exx(2:end-1,2:end-1)+Exx(2:end-1,3:end)+Exx(3:end,2:end-1)+Exx(3:end,3:end))/4;%average ETAxx to nodal points
    
    Eyy = zeros(Ny,Nx);
    Eyy(2:end,:) = diff(Vy(:,:),1,1)/dy;%Eyy in P points
    EyyNp = (Eyy(2:end-1,2:end-1)+Eyy(2:end-1,3:end)+Eyy(3:end,2:end-1)+Eyy(3:end,3:end))/4;%average ETAxx to nodal points
    
    ExyNp = 0.5 * (diff(Vx(2:end,2:end-1),1,1)/dy + diff(Vy(2:end-1,2:end),1,2)/dx);%Exy in internal nodal points, boundaries missing
    Exy = zeros(Ny,Nx);
    Exy(3:end-1,3:end-1) = (ExyNp(1:end-1,1:end-1)+ExyNp(1:end-1,2:end)+ExyNp(2:end,1:end-1)+ExyNp(2:end,2:end))/4;%average SIGxx to P points
    Exy(end,:)= 2 * Exy(end-1,:) -  Exy(end-2,:);%extrapolate to border rows
    Exy(:,end)= 2 * Exy(:,end-1) -  Exy(:,end-2);
    Exy(2,:)= 2 * Exy(3,:) -  Exy(4,:);
    Exy(:,2)= 2 * Exy(:,3) -  Exy(:,4);
    
    %second strain invariant
    Eii = sqrt(0.5 * (Exx.^2 + Eyy.^2) + Exy.^2);%strain rate inv. on P points

    %stress should be calculated on nodal pointe(where ETA is defined) to reduce artefacts
    SIGxxpNp = 2 * ETA(2:end-1,2:end-1) .* ExxNp;
    SIGyypNp = 2 * ETA(2:end-1,2:end-1) .* EyyNp;
    SIGxyNp = 2 * ETA(2:end-1,2:end-1) .* ExyNp;
    
    %then interpolate stress to P-points
    SIGxxp = zeros(Ny,Nx);
    SIGxxp(3:end-1,3:end-1) = (SIGxxpNp(1:end-1,1:end-1)+SIGxxpNp(1:end-1,2:end)+SIGxxpNp(2:end,1:end-1)+SIGxxpNp(2:end,2:end))/4;%average SIGxx to pressure points
    SIGxxp(end,:)= 2 * SIGxxp(end-1,:) -  SIGxxp(end-2,:);%extrapolate to border rows
    SIGxxp(:,end)= 2 * SIGxxp(:,end-1) -  SIGxxp(:,end-2);
    SIGxxp(2,:)= 2 * SIGxxp(3,:) -  SIGxxp(4,:);
    SIGxxp(:,2)= 2 * SIGxxp(:,3) -  SIGxxp(:,4);
    
    SIGyyp = zeros(Ny,Nx);
    SIGyyp(3:end-1,3:end-1) = (SIGyypNp(1:end-1,1:end-1)+SIGyypNp(1:end-1,2:end)+SIGyypNp(2:end,1:end-1)+SIGyypNp(2:end,2:end))/4;%average SIGyy to pressure points
    SIGyyp(end,:)= 2 * SIGyyp(end-1,:) -  SIGyyp(end-2,:);%extrapolate to border rows
    SIGyyp(:,end)= 2 * SIGyyp(:,end-1) -  SIGyyp(:,end-2);
    SIGyyp(2,:)= 2 * SIGyyp(3,:) -  SIGyyp(4,:);
    SIGyyp(:,2)= 2 * SIGyyp(:,3) -  SIGyyp(:,4);
    
    SIGxy = zeros(Ny,Nx);
    SIGxy(3:end-1,3:end-1) = (SIGxyNp(1:end-1,1:end-1)+SIGxyNp(1:end-1,2:end)+SIGxyNp(2:end,1:end-1)+SIGxyNp(2:end,2:end))/4;%average SIGxy to pressure points
    SIGxy(end,:)= 2 * SIGxy(end-1,:) -  SIGxy(end-2,:);%extrapolate to border rows
    SIGxy(:,end)= 2 * SIGxy(:,end-1) -  SIGxy(:,end-2);
    SIGxy(2,:)= 2 * SIGxy(3,:) -  SIGxy(4,:);
    SIGxy(:,2)= 2 * SIGxy(:,3) -  SIGxy(:,4);
      
    %second stress invariant
    SIGii = sqrt(0.5 * (SIGxxp.^2 + SIGyyp.^2) + SIGxy.^2);
      
    %full stress tensor
    SIGxx   = SIGxxp - P;
    SIGyy   = SIGyyp - P;

    
    %Calculate principal stress
    %
    %Calculate shear stress and normal stress in shear zone
    %http://www.continuummechanics.org/stressxforms.html for matrix notation
    %of EQ 12.18-12.20 from Gerya (2010)
    %redefine SZ-angle to have same reference frame (see example)
    TRang = 90-SZang;
    %
    for i=1:1:Ny
      for j=1:1:Nx
         ASIG = [SIGxxp(i,j) SIGxy(i,j) ; SIGxy(i,j) SIGyyp(i,j)];
         %principal stresses by eigen vector
         [VSIG, DSIG(i,j,:,:)] = eig(ASIG);
         SIG1deg(i,j) = atan2d(VSIG(2,1),VSIG(1,1));
         SIG3deg(i,j) = atan2d(VSIG(2,2),VSIG(1,2));
         %normal and shear stress
         SIGnTAU(i,j,:,:) = [cosd(TRang) sind(TRang); -sind(TRang) cosd(TRang)] *...
         ASIG * [cosd(TRang) -sind(TRang); sind(TRang) cosd(TRang)];
      end
    end

    
    SIGnSZ = SIGnTAU(:,:,1,1);
    TAUSZ = SIGnTAU(:,:,1,2);
    
    %Extract magnitude and direction of sig1 and sig3
    SIG1mag(:,:) = -DSIG(:,:,1,1)+P;
    SIG1mag_A(:,:,sc) = -DSIG(:,:,1,1)+P;
    SIG3mag(:,:) = -DSIG(:,:,2,2)+P;
    SIG3mag_A(:,:,sc) = -DSIG(:,:,2,2)+P;
    %normalize Sigma1 angle to obtain better plots
    SIG1degN =  SIG1deg;
    SIG1degN (SIG1deg<0) = SIG1deg(SIG1deg<0)+180;
   
        %EDIT7b
    %add readout of max - min Values - in calcite! -> for plots
    %MaxMinsig1=...
    
    %% Max-Min values extraction 
    %for matlab 2017 and higher maybe: [maxPval, maxPind] = maxk(P(RHO==3),round(nSZ*0.02))
    %define polygon strictly between pistons
    [NyG, NxG] = ndgrid(1:Ny, 1:Nx);
    Ny_overlap = abs(MCNx(6)-MCNx(7))/abs(MCNx(8)-MCNx(7))*abs(MCNy(8)-MCNy(7));
    inSZ = inpolygon(NyG, NxG, MCNy([7 7 6 6])'+[0 Ny_overlap 0 -Ny_overlap] ,MCNx([7 6 6 7])');
    
    %%% P
    Pdummy=zeros(size(P));  Pdummy(:,:)=NaN;    Pdummy(RHO==RHO_sz)=P(RHO==RHO_sz);%create a matrix containing only shear zone values
    Pdummylin=reshape(Pdummy,numel(Pdummy),1); %linearize this matrix, for sorting purpose
    num_PSZ=numel(Pdummylin)-sum(isnan(Pdummylin));%number of non NaN values in the matrix
    %sort and extract the top nMM fraction of max and min values, and according positions
    [P_FULLsort,P_FULLidx] = sort(Pdummylin);%https://ch.mathworks.com/matlabcentral/answers/48359-how-to-find-the-top-10-value-of-the-matrix
    [minPy, minPx]=ind2sub(size(Pdummy),P_FULLidx(1:nMM));%multiple min values locations - extract first element if needed
    [maxPy, maxPx]=ind2sub(size(Pdummy),P_FULLidx(num_PSZ-nMM:num_PSZ));%multiple max values locations
    %clf,pcolor(P), caxis([minPSZval maxPSZval]), shading flat, colorbar, hold on, plot(minPx,minPy,'oc'), plot(maxPx,maxPy,'om')
    P_SZsort = P_FULLsort(1:num_PSZ);%sample is full without NaNs (RHO def.)
    maxPSZval = median(P_SZsort(round(length(P_SZsort)*(1-nMM_range)):end));
    minPSZval = median(P_SZsort(1:round(length(P_SZsort)*nMM_range)));
    % min-max P in inner SZ region   
    Pdummy_inSZ=Pdummy;  Pdummy_inSZ(inSZ==0)=NaN;  
    Pdummy_inSZlin = sort( Pdummy_inSZ(~isnan(Pdummy_inSZ)));%inner SZ without NaNs (inSZ def.)
    maxP_inSZval = median(Pdummy_inSZlin( round(length(Pdummy_inSZlin)*(1-nMM_range)):end));
    minP_inSZval = median(Pdummy_inSZlin( 1:round(length(Pdummy_inSZlin)*nMM_range)));
    Pminmax(sc,:) = [minPSZval maxPSZval minP_inSZval maxP_inSZval];%full sample / in SZ
    %%% SIG1mag
    SIG1magdummy=zeros(size(SIG1mag));  SIG1magdummy(:,:)=NaN;    SIG1magdummy(RHO==RHO_sz)=SIG1mag(RHO==RHO_sz);%create a matrix containing only shear zone values
    SIG1magdummylin=reshape(SIG1magdummy,numel(SIG1magdummy),1); %linearize this matrix, for sorting purpose
    num_SIG1magSZ=numel(SIG1magdummylin)-sum(isnan(SIG1magdummylin));%number of non NaN values in the matrix
    %sort and extract the top nMM fraction of max and min values, and according positions
    [SIG1mag_FULLsort,SIG1mag_FULLidx] = sort(SIG1magdummylin);%https://ch.mathworks.com/matlabcentral/answers/48359-how-to-find-the-top-10-value-of-the-matrix
    maxSIG1magSZval = mean(SIG1mag_FULLsort(num_SIG1magSZ-nMM:num_SIG1magSZ));
    minSIG1magSZval = mean(SIG1mag_FULLsort(1:nMM));
    [minSIG1magy, minSIG1magx]=ind2sub(size(SIG1magdummy),SIG1mag_FULLidx(1:nMM));%multiple min values locations - extract first element if needed
    [maxSIG1magy, maxSIG1magx]=ind2sub(size(SIG1magdummy),SIG1mag_FULLidx(num_SIG1magSZ-nMM:num_SIG1magSZ));%multiple max values locations
    %clf,pcolor(SIG1mag), caxis([minSIG1magSZval maxSIG1magSZval]), shading flat, colorbar, hold on, plot(minSIG1magx,minSIG1magy,'oc'), plot(maxSIG1magx,maxSIG1magy,'om')
    SIG1magminmax(sc,:) = [minSIG1magSZval maxSIG1magSZval];
    %%% SIG3mag
    SIG3magdummy=zeros(size(SIG3mag));  SIG3magdummy(:,:)=NaN;    SIG3magdummy(RHO==RHO_sz)=SIG3mag(RHO==RHO_sz);%create a matrix containing only shear zone values
    SIG3magdummylin=reshape(SIG3magdummy,numel(SIG3magdummy),1); %linearize this matrix, for sorting purpose
    num_SIG3magSZ=numel(SIG3magdummylin)-sum(isnan(SIG3magdummylin));%number of non NaN values in the matrix
    %sort and extract the top nMM fraction of max and min values, and according positions
    [SIG3mag_FULLsort,SIG3mag_FULLidx] = sort(SIG3magdummylin);%https://ch.mathworks.com/matlabcentral/answers/48359-how-to-find-the-top-10-value-of-the-matrix
    maxSIG3magSZval = mean(SIG3mag_FULLsort(num_SIG3magSZ-nMM:num_SIG3magSZ));
    minSIG3magSZval = mean(SIG3mag_FULLsort(1:nMM));
    [minSIG3magy, minSIG3magx]=ind2sub(size(SIG3magdummy),SIG3mag_FULLidx(1:nMM));%multiple min values locations - extract first element if needed
    [maxSIG3magy, maxSIG3magx]=ind2sub(size(SIG3magdummy),SIG3mag_FULLidx(num_SIG3magSZ-nMM:num_SIG3magSZ));%multiple max values locations
    %clf,pcolor(SIG3mag), caxis([minSIG3magSZval maxSIG3magSZval]), shading flat, colorbar, hold on, plot(minSIG3magx,minSIG3magy,'oc'), plot(maxSIG3magx,maxSIG3magy,'om')
    SIG3magminmax(sc,:) = [minSIG3magSZval maxSIG3magSZval];
    %%% SIG1degN
    SIG1degNdummy=zeros(size(SIG1degN));  SIG1degNdummy(:,:)=NaN;    SIG1degNdummy(RHO==RHO_sz)=SIG1degN(RHO==RHO_sz);%create a matrix containing only shear zone values
    SIG1degNdummylin=reshape(SIG1degNdummy,numel(SIG1degNdummy),1); %linearize this matrix, for sorting purpose
    num_SIG1degNSZ=numel(SIG1degNdummylin)-sum(isnan(SIG1degNdummylin));%number of non NaN values in the matrix
    %sort and extract the top nMM fraction of max and min values, and according positions
    [SIG1degN_FULLsort,SIG1degN_FULLidx] = sort(SIG1degNdummylin);%https://ch.mathworks.com/matlabcentral/answers/48359-how-to-find-the-top-10-value-of-the-matrix
    maxSIG1degNSZval = mean(SIG1degN_FULLsort(num_SIG1degNSZ-nMM:num_SIG1degNSZ));
    minSIG1degNSZval = mean(SIG1degN_FULLsort(1:nMM));
    [minSIG1degNy, minSIG1degNx]=ind2sub(size(SIG1degNdummy),SIG1degN_FULLidx(1:nMM));%multiple min values locations - extract first element if needed
    [maxSIG1degNy, maxSIG1degNx]=ind2sub(size(SIG1degNdummy),SIG1degN_FULLidx(num_SIG1degNSZ-nMM:num_SIG1degNSZ));%multiple max values locations
    %clf,pcolor(SIG1degN), caxis([minSIG1degNSZval maxSIG1degNSZval]), shading flat, colorbar, hold on, plot(minSIG1degNx,minSIG1degNy,'oc'), plot(maxSIG1degNx,maxSIG1degNy,'om')
    SIG1degNminmax(sc,:) = [minSIG1degNSZval maxSIG1degNSZval];
  
    %%% SIGiid
    SIGiidummy=zeros(size(SIGii));  SIGiidummy(:,:)=NaN;    SIGiidummy(RHO==RHO_sz)=SIGii(RHO==RHO_sz);%create a matrix containing only shear zone values
    SIGiidummylin=reshape(SIGiidummy,numel(SIGiidummy),1); %linearize this matrix, for sorting purpose
    num_SIGiiSZ=numel(SIGiidummylin)-sum(isnan(SIGiidummylin));%number of non NaN values in the matrix
    %sort and extract the top nMM fraction of max and min values, and according positions
    [SIGii_FULLsort,SIGii_FULLidx] = sort(SIGiidummylin);%https://ch.mathworks.com/matlabcentral/answers/48359-how-to-find-the-top-10-value-of-the-matrix
    [minSIGiidy, minSIGiidx]=ind2sub(size(SIGiidummy),SIGii_FULLidx(1:nMM));%multiple min values locations - extract first element if needed
    [maxSIGiidy, maxSIGiidx]=ind2sub(size(SIGiidummy),SIGii_FULLidx(num_SIGiiSZ-nMM:num_SIGiiSZ));%multiple max values locations
    %clf,pcolor(SIGii), caxis([minSIGiiSZval maxSIGiiSZval]), shading flat, colorbar, hold on, plot(minSIGiix,minSIGiiy,'oc'), plot(maxSIGiix,maxSIGiiy,'om')
    SIGii_SZsort = SIGii_FULLsort(1:num_SIGiiSZ);%sample is full without NaNs (RHO def.)
    maxSIGiiSZval = median(SIGii_SZsort(round(length(SIGii_SZsort)*(1-nMM_range)):end));
    minSIGiiSZval = median(SIGii_SZsort(1:round(length(SIGii_SZsort)*nMM_range)));
    % min-max SIGiid in inner SZ region   
    SIGiidummy_inSZ=SIGiidummy;  SIGiidummy_inSZ(inSZ==0)=NaN;  
    SIGiidummy_inSZlin = sort( SIGiidummy_inSZ(~isnan(SIGiidummy_inSZ)));%inner SZ without NaNs (inSZ def.)
    maxSIGii_inSZval = median(SIGiidummy_inSZlin( round(length(SIGiidummy_inSZlin)*(1-nMM_range)):end));
    minSIGii_inSZval = median(SIGiidummy_inSZlin( 1:round(length(SIGiidummy_inSZlin)*nMM_range)));
    SIGiiminmax(sc,:) = [minSIGiiSZval maxSIGiiSZval minSIGii_inSZval maxSIGii_inSZval];%full sample / in SZ
    
    %%% EiiSUM
    EiiSUMdummy=zeros(size(EiiSUM));  EiiSUMdummy(:,:)=NaN;    EiiSUMdummy(RHO==RHO_sz)=EiiSUM(RHO==RHO_sz);%create a matrix containing only shear zone values
    EiiSUMdummylin=reshape(EiiSUMdummy,numel(EiiSUMdummy),1); %linearize this matrix, for sorting purpose
    num_EiiSUMSZ=numel(EiiSUMdummylin)-sum(isnan(EiiSUMdummylin));%number of non NaN values in the matrix
    %sort and extract the top nMM fraction of max and min values, and according positions
    [EiiSUM_FULLsort,EiiSUM_FULLidx] = sort(EiiSUMdummylin);%https://ch.mathworks.com/matlabcentral/answers/48359-how-to-find-the-top-10-value-of-the-matrix
    maxEiiSUMSZval = mean(EiiSUM_FULLsort(num_EiiSUMSZ-nMM:num_EiiSUMSZ));
    minEiiSUMSZval = mean(EiiSUM_FULLsort(1:nMM));
    [minEiiSUMy, minEiiSUMx]=ind2sub(size(EiiSUMdummy),EiiSUM_FULLidx(1:nMM));%multiple min values locations - extract first element if needed
    [maxEiiSUMy, maxEiiSUMx]=ind2sub(size(EiiSUMdummy),EiiSUM_FULLidx(num_EiiSUMSZ-nMM:num_EiiSUMSZ));%multiple max values locations
    %clf,pcolor(EiiSUM), caxis([minEiiSUMSZval maxEiiSUMSZval]), shading flat, colorbar, hold on, plot(minEiiSUMx,minEiiSUMy,'oc'), plot(maxEiiSUMx,maxEiiSUMy,'om')
    EiiSUMminmax(sc,:) = [minEiiSUMSZval maxEiiSUMSZval];
    %%% Eii
    Eiidummy=zeros(size(Eii));  Eiidummy(:,:)=NaN;    Eiidummy(RHO==RHO_sz)=Eii(RHO==RHO_sz);%create a matrix containing only shear zone values
    Eiidummylin=reshape(Eiidummy,numel(Eiidummy),1); %linearize this matrix, for sorting purpose
    num_EiiSZ=numel(Eiidummylin)-sum(isnan(Eiidummylin));%number of non NaN values in the matrix
    %sort and extract the top nMM fraction of max and min values, and according positions
    [Eii_FULLsort,Eii_FULLidx] = sort(Eiidummylin);%https://ch.mathworks.com/matlabcentral/answers/48359-how-to-find-the-top-10-value-of-the-matrix
    maxEiiSZval = mean(Eii_FULLsort(num_EiiSZ-nMM:num_EiiSZ));
    minEiiSZval = mean(Eii_FULLsort(1:nMM));
    [minEiiy, minEiix]=ind2sub(size(Eiidummy),Eii_FULLidx(1:nMM));%multiple min values locations - extract first element if needed
    [maxEiiy, maxEiix]=ind2sub(size(Eiidummy),Eii_FULLidx(num_EiiSZ-nMM:num_EiiSZ));%multiple max values locations
    %clf,pcolor(Eii), caxis([minEiiSZval maxEiiSZval]), shading flat, colorbar, hold on, plot(minEiix,minEiiy,'oc'), plot(maxEiix,maxEiiy,'om')
    Eiiminmax(sc,:) = [minEiiSZval maxEiiSZval];
    %%% ETA
    ETAdummy=zeros(size(ETA));  ETAdummy(:,:)=NaN;    ETAdummy(RHO==RHO_sz)=ETA(RHO==RHO_sz);%create a matrix containing only shear zone values
    ETAdummylin=reshape(ETAdummy,numel(ETAdummy),1); %linearize this matrix, for sorting purpose
    num_ETASZ=numel(ETAdummylin)-sum(isnan(ETAdummylin));%number of non NaN values in the matrix
    %sort and extract the top nMM fraction of max and min values, and according positions
    [ETA_FULLsort,ETA_FULLidx] = sort(ETAdummylin);%https://ch.mathworks.com/matlabcentral/answers/48359-how-to-find-the-top-10-value-of-the-matrix
    maxETASZval = mean(ETA_FULLsort(num_ETASZ-nMM:num_ETASZ));
    minETASZval = mean(ETA_FULLsort(1:nMM));
    [minETAy, minETAx]=ind2sub(size(ETAdummy),ETA_FULLidx(1:nMM));%multiple min values locations - extract first element if needed
    [maxETAy, maxETAx]=ind2sub(size(ETAdummy),ETA_FULLidx(num_ETASZ-nMM:num_ETASZ));%multiple max values locations
    %clf,pcolor(ETA), caxis([minETASZval maxETASZval]), shading flat, colorbar, hold on, plot(minETAx,minETAy,'oc'), plot(maxETAx,maxETAy,'om')
    ETAminmax(sc,:) = [minETASZval maxETASZval];
    %EDIT1a
    PdSIG1minmax(sc,1)=min(min(Pdummy./SIG1magdummy));
    PdSIG1minmax(sc,2)=max(max(Pdummy./SIG1magdummy));
    

             
    %% save Viscosity and strain rate within SZ for mean and mode
        
    %analysys (and development) -> histplots.m script
    EiiSZdummy = zeros(nSZ,1); EiiSZdummy(:)=NaN; EiiSZdummy(1:length(Eii(RHO==RHO_sz))) = Eii(RHO==RHO_sz); %RHO_sz-0.01 < RHO & RHO < RHO_sz+0.01
    EiiSZ(sc,:) = EiiSZdummy;
    ETASZdummy = zeros(nSZ,1); ETASZdummy(:)=NaN; ETASZdummy(1:length(ETA(RHO==RHO_sz))) = ETA(RHO==RHO_sz);
    ETASZ(sc,:) = ETASZdummy;
    
    EiiSZdummy = zeros(nSZ,1); EiiSZdummy(:)=NaN; EiiSZdummy(1:length(Eii(inSZ==1 & RHO==RHO_sz))) = Eii(inSZ==1 & RHO==RHO_sz); %RHO_sz-0.01 < RHO & RHO < RHO_sz+0.01
    Eii_inSZ(sc,:) = EiiSZdummy;
    ETASZdummy = zeros(nSZ,1); ETASZdummy(:)=NaN; ETASZdummy(1:length(ETA(inSZ==1 & RHO==RHO_sz))) = ETA(inSZ==1 & RHO==RHO_sz);
    ETA_inSZ(sc,:) = ETASZdummy;

    for schi=1:size(ETASZ,1)
    %find mode (better than mean) of Viscosity and strain rate
    [countsEii, binsEii]=histcounts(log10(EiiSZ(schi,:)),200);
    [ixEii, idxEii] = max(countsEii);
    ModeMeanEiiSZ(schi,1) = 10^(binsEii(idxEii)/2+binsEii(idxEii+1)/2);
    
    [countsETA, binsETA]=histcounts(log10(ETASZ(schi,:)),200);
    [ixETA, idxETA] = max(countsETA);
    ModeMeanETASZ(schi,1) = 10^(binsETA(idxETA)/2+binsETA(idxETA+1)/2);
    %mean in SZ
    ModeMeanEiiSZ(schi,2) = mean(EiiSZ(schi,:),'omitnan');
    ModeMeanETASZ(schi,2) = mean(ETASZ(schi,:),'omitnan');
    end

    
    %%%calculate and check force balance
%     sum(SIGyy,2);%sum lines of force (yy)
%     sum(SIGxx,1);%sum columns of force (xx)
        
        
    %% extract coordinates and plot size for Pressure/Stress plots
    
    %find internal BC - as defined in model - if changed in model adapt here
    %otherwise try to search for sigmaXX value jump
    uPty = round((MCNy(1)+MCNy(2))/2);
    uPtYV(sc) = uPty; 
    cPtBC=fix((iBC+uPty)/2);%y coord: middle between internal BC and Pt capsule - readout for bulk force
    lAlx = min(find((RHO(iBC,:)>6.5) & (RHO(iBC,:)<7.1)));% left x coord for bulk force readout (in piston)
    rAlx = max(find((RHO(iBC,:)>6.5) & (RHO(iBC,:)<7.1)));% right x coord -"-

    %% bulk experiment values
    EXPsig1(sc) = abs(mean(mean(SIGyy(cPtBC-2:cPtBC+2,lAlx:rAlx))));
    EXPdevsig1(sc) = abs(mean(mean(SIGyyp(cPtBC-2:cPtBC+2,lAlx:rAlx))));
    
    aSZang      = abs(SZang);%absolute of shear zone angel (for geometrical calculations)
    Apiston = MCx(4) - MCx(3);%piston width
    Ainter = Apiston/cosd(aSZang);%shear zone interface area
    
    
    %% extract stress values at readout points and calculate strains 
    %different loops for PS and SS geometry
    %EDIT8
    %readout points for pure shear experiments (array in shearzone)
    % ----------------
    % |1  2  3       | 50 95% vertical 
    % |4  5  6       | 0 25 50 % horizontal
    % |              |
    % ----------------
    if SZang == 0 %&&  elliONOFF ==0 %enter pure shear zone mode
        
                %%create array of readout points to linerarize positions - for simple shear mode
        %array of 5 points at vertical center of SZ, center to lower piston right end
        PElSZ = MCNx(5); PErSZ = MCNx(6);
        PEbSZ = MCNy(7); PEtSZ = MCNy(5);
                
        probe_Cy    =   round([ (PEtSZ-PEbSZ)*0.95+PEbSZ, ...
                                (PEtSZ-PEbSZ)*0.95+PEbSZ, ...
                                (PEtSZ-PEbSZ)*0.95+PEbSZ, ...
                                (PEtSZ-PEbSZ)*0.5+PEbSZ, ...
                                (PEtSZ-PEbSZ)*0.5+PEbSZ, ...
                                (PEtSZ-PEbSZ)*0.5+PEbSZ]);%y-coords
        probe_Cx    =   round([ (PErSZ-PElSZ)*0+PElSZ, ...
                                (PErSZ-PElSZ)*0.25+PElSZ, ...
                                (PErSZ-PElSZ)*0.5+PElSZ, ...
                                (PErSZ-PElSZ)*0.+PElSZ, ...
                                (PErSZ-PElSZ)*0.25+PElSZ, ...
                                (PErSZ-PElSZ)*0.5+PElSZ]);%x-coords
  
        %%%calculate strain (sum of incremantal strains - like rigB program (Renee Heilbronner)
        gSZlength(sc) = MCy(15)-MCy(16);%thickness of SZ
        if sc>1%calculate difference of SZ length = displacement in Shear direction
            gSZds_inc(sc) = gSZlength(sc-calcintv) - gSZlength(sc);
            ginc(sc) = ginc(sc-calcintv) + gSZds_inc(sc)/gSZlength(sc);%add 
        else% 0 displacemt in 1st step
            gSZds_inc(sc) = 0;
            ginc(sc) = 0;
        end
        ginc_val = ginc(1:calcintv:end);%extract only values for calculated steps, later interpolate

        if calcintv == 1
            ginc_smooth = smooth(ginc,10);%moving average smoothing of calculated strain, number is amount used for moving avg.
        else
            ginc_smooth = interp1(1:calcintv:max(savecycle),ginc_val,1:1:max(savecycle));%interpolate strain if bigger steps are used
        end
    
        %prepare ratio calculation for supported SZ width
        dSZwidthVc(sc)  =   MCNx(14)-MCNx(13);%SZ width at vertical center
        dSZsupSIG1(sc)  =   sum(SIG1mag(MCNy(14),MCNx(13):MCNx(14))>(Pconf+EXPsig1(sc))/2)-1;%sig1 supporting area
        dSZsupP(sc)     =   sum(P(MCNy(14),MCNx(13):MCNx(14))>Pconf)-1;%P supporting area
        dNOMsup(sc)     =   dSZwidthVc(1)/(gSZlength(sc)/gSZlength(1));
        
        %starin from markers
        %diff(MCyPL(5,:)-MCyPL(7,:))./(MCyPL(5,1:end-1)-MCyPL(7,1:end-1))./elapsedTimePL(2:end) %strain rate
        

    end%will later be linearized and read out

        
    %EDIT5
    %readout points for simple shear experiments (array in shearzone)
    if SZang ~= 0 %enter simple shear zone mode  
        
        SZthickNmax    =   round(1.5*SZthick*Ny/ysize);%maximum thickness of search for shear zone
        
        % Defining lagrangian marker points on corners to track motion of pistons
        % and shear zone
        % simple-shear  % pure-shear (coaxial)
        %   |  |        %   |  |
        %  1----2       %  1----2 
        %  3----4       %  3----4   
        %   |  |6       %  5|--|6  15 at center (5-15-6)
        %   | /|14      %   |  | 
        %  5|//|8       % 13|  |14 
        % 13|/ |        %   |  |
        %  7|  |        %  7|--|8  16 at center (7-16-8)
        %  9----10      %  9----10 
        % 11----12      % 11----12   
        %   |  |        %   |  | 
        
        SZthickNl   =   sum(RHO(MCNy(7):MCNy(7)+SZthickNmax,MCNx(7)+2)<RHO_sz+2);%SZ thinckness ax 2 px right of piston overlap (more stable)
        SZthickNr   =   sum(RHO(MCNy(6)-SZthickNmax:MCNy(6),MCNx(6)-2)<RHO_sz+2);%SZ thinckness ax 2 px left of piston overlap (more stable)
             
        MC5ay       =   MCNy(7)+ SZthickNl;
        MC8ay       =   MCNy(6)- SZthickNr;
        
        %array of 5 points at vertical center of SZ, left to right end
        probe_Cx    =   round(linspace(MCNx(7),MCNx(6),5));%x-coords
        probe_Cy    =   fix(linspace(round((MCNy(7)+MC5ay)/2),round((MCNy(6)+MC8ay)/2),5))-1;%y-coords
        
        %%create array of readout points to linerarize positions - for simple shear mode
        %array of 5 points at vertical center of SZ, center to lower piston right end
        probe_CxA(sc,:)    =   (linspace(MCNx(7),MCNx(6),5));%x-coords
        probe_CyA(sc,:)    =   (linspace(round((MCNy(7)+MC5ay)/2),round((MCNy(6)+MC8ay)/2),5))-1;%y-coords
        
        %piston tips readout - structure as above for center diagonal
        probe_Tx    =   [MCNx(7)-1 MCNx(7)+1 MCNx(6)-1 MCNx(6)+1];%x-coords
        probe_Ty    =   [ceil(MCNy(7)+0.9*SZthickNl) floor(MCNy(7)+0.1*SZthickNl) ...
                         ceil(MCNy(6)-0.1*SZthickNr) floor(MCNy(6)-0.9*SZthickNr)];
        
        probe_TxA(sc,:)    =   [MCNx(7)-1 MCNx(7)+1 MCNx(6)-1 MCNx(6)+1];%x-coords
        probe_TyA(sc,:)    =   [(MCNy(7)+0.8*SZthickNl) (MCNy(7)+0.1*SZthickNl)...
                                (MCNy(6)-0.1*SZthickNr) (MCNy(6)-0.8*SZthickNr)];


        if(sc == max(savecycle) && calcintv ==1 && maxminEND==0 && scPlotBigFig(1) == 0)
            %create matrix indices in 3D to collect moving readout points at each timestep
            probe_CpA = sub2ind(size(PrPL),round([smooth(probe_CyA(:,1)) smooth(probe_CyA(:,2)) smooth(probe_CyA(:,3))...
                smooth(probe_CyA(:,4)) smooth(probe_CyA(:,5))]),...
                round([smooth(probe_CxA(:,1)) smooth(probe_CxA(:,2)) smooth(probe_CxA(:,3)) smooth(probe_CxA(:,4)) smooth(probe_CxA(:,5))]),...
                [([1:max(savecycle)]') ([1:max(savecycle)]') ([1:max(savecycle)]') ([1:max(savecycle)]') ([1:max(savecycle)]')]);
            
            probe_TpA = sub2ind(size(PrPL),round([smooth(probe_TyA(:,1)) smooth(probe_TyA(:,2)) smooth(probe_TyA(:,3))...
                smooth(probe_TyA(:,4))]),...
                [ceil(smooth(probe_TxA(:,1))) floor(smooth(probe_TxA(:,2))) ceil(smooth(probe_TxA(:,3))) floor(smooth(probe_TxA(:,4)))],...
                [([1:max(savecycle)]') ([1:max(savecycle)]') ([1:max(savecycle)]') ([1:max(savecycle)]')]);

            %extract values Pressure dig center
            P_probe_Cp(:,:)=PrPL(probe_CpA);%SZ center array
            %extract values Sigma1
            Sig1_probe_Cp(:,:)=SIG1mag_A(probe_CpA);%SZ center array
            %extract values Sigma3
            Sig3_probe_Cp(:,:)=SIG3mag_A(probe_CpA);%SZ center array
            %extract values Pressure tips
            P_probe_Tp(:,:)=PrPL(probe_TpA);%SZ center array
            %extract values Sigma1
            Sig1_probe_Tp(:,:)=SIG1mag_A(probe_TpA);%SZ center array
            %extract values Sigma3
            Sig3_probe_Tp(:,:)=SIG3mag_A(probe_TpA);%SZ center array
            %extract values Sigma3
            Sigii_probe_Tp(:,:)=SIGiiPL(probe_TpA);%SZ center array
        end


        %%%calculate shear strain (sum of incremantal shear strains - like rigB program (Renee Heilbronner)
        gSZth(sc) = ((SZthickNl + SZthickNr )/2) * dy * cosd(SZang);%thickness of SZ - averaged on both ends
        gSZlength(sc) = (MCx(6)-MCx(7)) / cosd(SZang);%length of SZ
        if sc>1%calculate difference of SZ length = displacement in Shear direction
            gSZds_inc(sc) = gSZlength(sc) - gSZlength(sc-calcintv);
            ginc(sc) = ginc(sc-calcintv) + gSZds_inc(sc)/gSZth(sc);%add 
        else% 0 displacemt in 1st step
            gSZds_inc(sc) = 0;
            ginc(sc) = 0;
        end
        ginc_val = ginc(1:calcintv:end);%extract only values for calculated steps, later interpolate

        if calcintv == 1
            ginc_smooth = - smooth(ginc,10);%moving average smoothing of calculated strain, number is amount used for moving avg.
        else
            ginc_smooth = - interp1(1:calcintv:max(savecycle),ginc_val,1:1:max(savecycle));%interpolate strain if bigger steps are used
        end
    end% simple shear zone calculations

    %if elliONOFF == 0  %extract value for pure and simple shear arrays 
        %calculate corestponding matrix positions: equals sub2ind function
        probe_Cp=(probe_Cx-1)*Ny+probe_Cy;%center array
        %extract values
        P_probe_conf(sc)=P(5,5);%point in salt near corner
        
%         P_probe_Cp(sc,:)=P(probe_Cp);%SZ center array
%         Sig1_probe_Cp(sc,:)=SIG1mag(probe_Cp);%SZ center array
%         Sig3_probe_Cp(sc,:)=SIG3mag(probe_Cp);%SZ center array
    %end
    
   
    %% calculate ACF correction (gemoetrical correction for piston overlap - asdone for rig data)
    acfSZlength=abs(0.0066/cosd(SZang)+SZthick);%interface length + SZthick; 0.0066 is piston diamater [m]
    acfSZdisp=abs(linspace(0,exptimePL(end)*Vyset/cosd(90-SZang),savecycle));%displacemant parallel to shear direction
    acfMAXfactor=1/((acfSZlength-acfSZdisp(end))/acfSZlength);
    acfFactor_OLD=ones(1,savecycle);
    acfFactor_OLD(acfSZdisp>SZthick)=linspace(1,acfMAXfactor,length(acfFactor_OLD(acfSZdisp>SZthick)));
    
    acfcontact = acfSZdisp ./ acfSZlength ;
    acfFactor = cosd(acfcontact*90).^2;
    
    %define exp-force acf for c-axis in plots
    if SZang == 0
        EXPsig1_acf = EXPsig1;%./(2-gSZlength/gSZlength(1));%standard cylinder
        EXPsig1_acfp7 = EXPsig1./(2-((gSZlength/gSZlength(1)).^7/4+0.75));%best fit power 7 
        EXPsig3_acf = ones(1,length(EXPsig1_acf))*Pconf;
    elseif SZang ~= 0 %RigB procedure follows here
        EXPsig1_acf_OLD = P_probe_conf+(EXPsig1-P_probe_conf).*acfFactor_OLD;
                
        FNA = (EXPsig1 * Apiston * sind(aSZang))  +  (P_probe_conf * Apiston * cosd(aSZang));%FNsig1 + FNsig3
        FSA = (EXPsig1 * Apiston * cosd(aSZang))  -  (P_probe_conf * Apiston * sind(aSZang));% FSsig1 + FSsig3

        FNpc = P_probe_conf .* Ainter .* (1-acfFactor);%confining pressure supporting overlapping piston tip
        FNAcorr = FNA - FNpc;
        
        EXPsign = FNAcorr./(Ainter*acfFactor);
        EXPtau  = FSA./(Ainter*acfFactor);
        
        EXPTaumax = EXPtau ./ cosd(90 - 2*aSZang); 
        
        EXPsig1_acf = EXPsign + EXPTaumax;
        EXPdsig_acf = 2 * EXPTaumax;
        EXPsig3_acf = EXPsign - EXPTaumax;
    end
    
    %% interpolate Vx/Vy to readable arrow size (plot preparation)
    
    Vxq = interp2(msX,msY,Vx, msXq,msYq);
    Vyq = interp2(msX,msY,Vy, msXq,msYq);
    
    %% % maxminEND readout breaking execution and returning to start
    if maxminEND == 1
        maxminEND = 0;%dont enter maxminEND readout again
        continue%break for loop here and return to start
    end 
        
    
%% Plot cycle --------------------------
%%remap colors for Material
RHO((RHO_al)<RHO&RHO<(RHO_elli+1)) = RHO_al-1;%ellipse
RHO(RHO>(RHO_elli)) = RHO_al-3;%platinum
%-------%%% plotting single Figs
FSS = 12;
FSSs = 8;
LWW = 2;
LWW2 = 1.5;%define line-width for graphs
gy = [0.7 0.7 0.7];%grey line color

% axis and rotation definition 
tangle = SZang*pi/180;%matrix rotation angle
P_nom(sc) = (EXPsig1(sc) + P_probe_conf(sc)) /2;%define P_nom without ACF
P_acf(sc) = (EXPsig1_acf(sc) + EXPsig3_acf(sc)) /2;%define P_acf with ACF
Dsig_nom(sc) = (EXPsig1(sc) - P_probe_conf(sc));%define P_nom witout ACF
Dsig_acf(sc) = (EXPsig1_acf(sc) - EXPsig3_acf(sc));%define P_nom witout ACF
CzoomIntv = [-round(1*Dsig_acf(sc),-7)*0.75 round(1*Dsig_acf(sc),-7)*0.75];%c-axis width is full ACF corrected differential stress

%define center of SZ in rotated coordinate frame [mm]
SZcenter    = xyscale*[(MCyPL(13,sc) + MCyPL(14,sc))/2; (MCxPL(13,sc) + MCxPL(14,sc))/2];
RotM = [cosd(SZang) sind(SZang);-sind(SZang) cosd(SZang)];
SZrotCenter = round(RotM*SZcenter,2);%xyscale*dx*[probe_Cy(3); probe_Cx(3)],2);%center
%sample points in rot frame
SZrotPC =  round(RotM*(xyscale*[dy*probe_CyA(sc,:); dx*probe_CxA(sc,:)]),2);
SZrotPT =  round(RotM*(xyscale*[dy*probe_TyA(sc,:); dx*probe_TxA(sc,:)]),2);
%definition of rotation: https://ch.mathworks.com/help/images/matrix-representation-of-geometric-transformations.html
if SZang ~= 0
    SZmaxExtent = round(xyscale*diff(MCxPL([5 8],savecycle))/cosd(SZang),2);%lateral width beween piston corners
    SZzoomaxis  = [SZrotCenter(2)+0.5*[-SZmaxExtent SZmaxExtent] SZrotCenter(1)+xyscale*0.9*[-SZthick SZthick]];%zoom axis definitionis
else
    SZmaxExtent = round(xyscale*diff(MCxPL([13 14],savecycle))/cosd(SZang),2);%lateral width beween piston corners
    SZzoomaxis  = [SZrotCenter(2)+0.65*[-SZmaxExtent SZmaxExtent] SZrotCenter(1)+xyscale*0.55*[-SZthick SZthick]];%zoom axis definitionis
end

if scPlotBigFig(1) ~= 0 && plotsingle ~= 5%only plot single if requested
%% plot single - full siuze

%paper size for single figs
if SZang == 0
    PaperSizeZoom = ([0 0 20 20]); %size of paper %nature half-col: [0 0 8.9 12.3] % anture full-col: [0 0 18.3 12.3] height: 247
else
    PaperSizeZoom = ([0 0 30 10]); %size of paper %typically: 0 0 30 10
end

%P-distribution in 2D
gcf=figure(100+sc); clf
siPax   =   round([P_acf(sc)/PSscale+CzoomIntv/PSscale],1);

contourf(x(2:end),y,(P(:,2:end))/PSscale,(0.8*siPax(1) :25: 1.2*siPax(2)),'edgecolor','none')
hold on
plot([lAlx:rAlx]*dx*xyscale, cPtBC*dy*xyscale, '.r');%highlight force piston measurement zone
% plot(minPx/scX,minPy/scY,'.','color',[0.1 0.8 0.9],'markersize',5)%readout points max-min 
% plot(maxPx/scX,maxPy/scY,'.','color',[0.8 0.1 0.1],'markersize',5)%readout points max-min                                                                                   -
shading flat; axis image; axis equal; caxis(siPax)


xlabel('diameter [mm]','FontSize',FSS); ylabel('height [mm]','FontSize',FSS)
%title('Pressure $[GPa]$','FontSize',FSS); 
%colorbar; 
colormap(viridis); %colormap(jet); %for different colorbar, add in octave
axis off, axis tight
drawnow;

set(gcf,'Units','centimeters');
set(gcf, 'Position', [0 0 8.9 12.3]);
pos = get(gcf,'Position');
set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[pos(3), pos(4)])
if savefigs==1
    print(gcf,sprintf('P-OV_%s_%d_vir_minmax_br.png', ModName, sc),'-dpng', '-r600')
end

%% Material

gcf=figure(400+sc); clf
pcolor(x,y,RHO)
shading flat; hold on
%plot(MCx*xyscale,MCy*xyscale,'.m')%zone definitioin points
% if SZang~=0
% plot([MCNx(7) MCNx(6)]*dx*xyscale,[MC5ay MC8ay]*dy*xyscale,'.r')
% end
plot([lAlx:rAlx]*dx*xyscale, iBC*dy*xyscale, '.r');%highlight internal bc zone
%plot([lAlx:rAlx]*dx*xyscale, cPtBC*dy*xyscale, '.k');%highlight force piston measurement zone
quiver(x(1:qs:end),y(1:qs:end),Vxq,Vyq,0.7,'k')

hold on; contour(x,y,inSZ)

axis image; axis equal;
xlabel('diameter [mm]','FontSize',FSS); ylabel('height [mm]','FontSize',FSS)
%title('Materials of model','FontSize',FSS)
colormap(MatCmap);
axis off, axis tight
drawnow;

set(gcf,'Units','centimeters');
set(gcf, 'Position', [0 0 8.9 12.3]);
pos = get(gcf,'Position');
set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[pos(3), pos(4)])
if savefigs==1
    print(gcf,sprintf('Material_%s_%d.png', ModName, sc),'-dpng', '-r600')%
end
%%
if sc == max(scPlotBigFig)%zoom figs only in  last stage
%% SZ zoom stack distribution in 2D zoom

%settings for contour plots
nEx = 50; %color intervals
contIV = 25; % contour interval
% Ptrans = 1.54e9; % defined at start
if Ptrans_scale == 1%manual definition of c-axis scale and white center of bwr colorbar
    P_Cax   =   (Ptrans + [-120e6 130e6])/PSscale;%round([Ptrans/PSscale*0.7 Ptrans/PSscale*1.3],1);%[1.48 1.79];%
    s1_Cax  =   P_Cax;
    s3_Cax  =   P_Cax;
    P_bwrC  =   Ptrans/PSscale;%colorbar center mean(P_Cax);%
    s1_bwrC =   P_bwrC;
    s3_bwrC =   P_bwrC;
    lEx = 1 - raPTpm; % lower extension of contour levels rel to c-axis
    hEx = 1 + raPTpm; % upper extension
else
    P_Cax   =   round([P_acf(sc)/PSscale+CzoomIntv/PSscale],1);
    s1_Cax  =   round([EXPsig1_acf(sc)/PSscale+CzoomIntv/PSscale],1);
    s3_Cax  =   round([EXPsig3_acf(sc)/PSscale+CzoomIntv/PSscale],1);  
    P_bwrC  =   P_acf(sc)/PSscale;%colorbar center
    s1_bwrC =   EXPsig1_acf(sc)/PSscale;
    s3_bwrC =   EXPsig3_acf(sc)/PSscale;
    lEx = 0.8; % lower extendion rel to c-axis
    hEx = 1.2; % upper extension
end

% -- %%% Stresses Sig1 , P, Sig3
gcf=figure(200+sc); clf
ax1 = subplot(3,1,1);
t=hgtransform;
contourf(x,y,(SIG1mag)/PSscale,(lEx*s1_Cax(1) :contIV: hEx*s1_Cax(2)),'edgecolor','none','Parent',t); caxis(s1_Cax); 
colorbar; colormap(ax1,bluewhitered(s1_bwrC));
title('$\sigma_1$ $[GPa]$','FontSize',FSS)
axis equal; shading interp;
%outline of phases
hold on; contour(x,y,RHO,[1.1 1.1],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[1.1 1.1],'--','edgecolor','w','linewidth',LWW,'Parent',t);
hold on; contour(x,y,RHO,[6 6],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[6 6],'--','edgecolor','w','linewidth',LWW,'Parent',t);
%xlabel('width [mm]','FontSize',FSS); ylabel('height [mm]','FontSize',FSS)
 plot(SZrotPC(2,3),SZrotPC(1,3),'k.','markersize',15);%sample points
 plot(SZrotPC(2,3),SZrotPC(1,3),'w.','markersize',5);%sample points
Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
axis(SZzoomaxis); axis off

ax2 = subplot(3,1,2);
t=hgtransform;
contourf(x,y,(P)/PSscale,(lEx*P_Cax(1) :contIV: hEx*P_Cax(2)),'edgecolor','none','Parent',t); caxis(P_Cax);
colorbar; colormap(ax2, bluewhitered(P_bwrC));
title('$P$ $[GPa]$','FontSize',FSS)
axis equal; shading interp;
%outline of phases
hold on; contour(x,y,RHO,[1.1 1.1],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[1.1 1.1],'--','edgecolor','w','linewidth',LWW,'Parent',t);
hold on; contour(x,y,RHO,[6 6],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[6 6],'--','edgecolor','w','linewidth',LWW,'Parent',t);
% plot(SZrotPC(2,1),SZrotPC(1,1),'k<');%sample points
% plot(SZrotPC(2,1),SZrotPC(1,1),'w<','markersize',3);%sample points
% plot(SZrotPC(2,2),SZrotPC(1,2),'kx');%sample points
% plot(SZrotPC(2,2),SZrotPC(1,2),'wx','markersize',3);%sample points
 plot(SZrotPC(2,3),SZrotPC(1,3),'k.','markersize',15);%sample points
 plot(SZrotPC(2,3),SZrotPC(1,3),'w.','markersize',5);%sample points
% % plot(SZrotPC(2,4),SZrotPC(1,4),'kd');%sample points
% % plot(SZrotPC(2,4),SZrotPC(1,4),'wd','markersize',3);%sample points
% % plot(SZrotPC(2,5),SZrotPC(1,5),'k>');%sample points
% % plot(SZrotPC(2,5),SZrotPC(1,5),'w>','markersize',3);%sample points
% %xlabel('width [mm]','FontSize',FSS); ylabel('height [mm]','FontSize',FSS)


Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
axis(SZzoomaxis); axis off

ax3 = subplot(3,1,3);
t=hgtransform;
contourf(x,y,(SIG3mag)/PSscale,(lEx*s3_Cax(1) :contIV: hEx*s3_Cax(2))-3,'edgecolor','none','Parent',t); caxis(s3_Cax); 
colorbar; colormap(ax3, bluewhitered(s3_bwrC)); 
title('$\sigma_3$ $[GPa]$','FontSize',FSS)
axis equal; shading interp;
%outline of phases
hold on; contour(x,y,RHO,[1.1 1.1],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[1.1 1.1],'--','edgecolor','w','linewidth',LWW,'Parent',t);
hold on; contour(x,y,RHO,[6 6],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[6 6],'--','edgecolor','w','linewidth',LWW,'Parent',t);
%xlabel('width [mm]','FontSize',FSS); ylabel('height [mm]','FontSize',FSS)
 plot(SZrotPC(2,3),SZrotPC(1,3),'k.','markersize',15);%sample points
 plot(SZrotPC(2,3),SZrotPC(1,3),'w.','markersize',5);%sample points
Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
axis(SZzoomaxis); axis off

set(gcf,'Units','centimeters');
set(gcf, 'Position', PaperSizeZoom); %size of paper
pos = get(gcf,'Position');
set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[pos(3), pos(4)])
if savefigs==1
    print(gcf,sprintf('P-zoom_%s_%d.png', ModName, sc),'-dpng', '-r600')
end

%% -- %%% P and Grain-size
gcf=figure(260+sc); clf

ax1 = subplot(3,1,1);
t=hgtransform;
contourf(x,y,(P)/PSscale,(0.8*siPax(1) :contIV: 1.2*siPax(2)),'edgecolor','none','Parent',t); caxis(P_acf(sc)/PSscale+CzoomIntv/PSscale);
colorbar; colormap(ax1,viridis);
title('$P$ $[GPa]$','FontSize',FSS)
axis equal; shading interp;
%outline of phases
hold on; contour(x,y,RHO,[1.1 1.1],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[1.1 1.1],'--','edgecolor','w','linewidth',LWW,'Parent',t);
hold on; contour(x,y,RHO,[6 6],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[6 6],'--','edgecolor','w','linewidth',LWW,'Parent',t);
%xlabel('width [mm]','FontSize',FSS); ylabel('height [mm]','FontSize',FSS)
Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
axis(SZzoomaxis); axis off

GSaR_m =        2239;
GSbR_m =     -0.8432;

GSaR_r =       848.7;
GSbR_r =     -0.8872;

ax2 = subplot(3,1,2);%Rutter average flow-law
t=hgtransform;
contourf(x,y,(((SIGii*2/1e6)/GSaR_m).^(1/GSbR_m)+((SIGii*2/1e6)/GSaR_r).^(1/GSbR_r))/2,[0:1:200],'edgecolor','none','Parent',t); caxis([4 35]); 
colorbar; colormap(ax2, flipud(viridis)); %for different colorbar, add in octave
title('Grain-size $[\mu m]$','FontSize',FSS)%- reference (white): $\sigma_{yy}$ with overlap corr.
axis equal; shading interp;
%outline of phases
hold on; contour(x,y,RHO,[1.1 1.1],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[1.1 1.1],'--','edgecolor','w','linewidth',LWW,'Parent',t);
hold on; contour(x,y,RHO,[6 6],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[6 6],'--','edgecolor','w','linewidth',LWW,'Parent',t);
Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
axis(SZzoomaxis); axis off

ax3 = subplot(3,1,3);
t=hgtransform;
contourf(x,y,(2*SIGii)/PSscale,[0:5:400],'edgecolor','none','Parent',t); caxis(Dsig_acf(sc)/1e6+0.5*CzoomIntv/PSscale); 
colorbar; colormap(ax3, viridis); %for different colorbar, add in octave
title('$\Delta\sigma$ $[MPa]$','FontSize',FSS)%- reference (white): $\sigma_{yy}$ with overlap corr.
axis equal; shading interp;
%outline of phases
hold on; contour(x,y,RHO,[1.1 1.1],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[1.1 1.1],'--','edgecolor','w','linewidth',LWW,'Parent',t);
hold on; contour(x,y,RHO,[6 6],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[6 6],'--','edgecolor','w','linewidth',LWW,'Parent',t);
Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
axis(SZzoomaxis); axis off
drawnow

set(gcf,'Units','centimeters');
set(gcf, 'Position', PaperSizeZoom); %size of paper
pos = get(gcf,'Position');
set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[pos(3), pos(4)])
if savefigs==1
    print(gcf,sprintf('P_vir-GS_%s_%d.png', ModName, sc),'-dpng', '-r600')
end

%% -- %%% D-Stress and Strain
gcf=figure(220+sc); clf

% gcf.Renderer ='Painters'; to get vector output: https://de.mathworks.com/matlabcentral/answers/2755-printing-figure-to-pdf-produces-bitmap-instead-of-vector
ax1 = subplot(3,1,1);
t=hgtransform;
contourf(x,y,(2*SIGii)/PSscale,[0:25e6:Dsig_acf(sc)*2]/PSscale,'edgecolor','none','Parent',t); caxis((Dsig_acf(sc)+CzoomIntv)/PSscale); 
colorbar; colormap(ax1, bluewhitered(Dsig_acf(sc)/1e6)); %for different colorbar, add in octave
title('$\Delta\sigma$ $[MPa]$','FontSize',FSS)%- reference (white): $\sigma_{yy}$ with overlap corr.
axis equal; shading interp;
%outline of phases
hold on; contour(x,y,RHO,[1.1 1.1],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[1.1 1.1],'--','edgecolor','w','linewidth',LWW,'Parent',t);
hold on; contour(x,y,RHO,[6 6],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[6 6],'--','edgecolor','w','linewidth',LWW,'Parent',t);
% plot(SZrotPC(2,1),SZrotPC(1,1),'k<');%sample points
% plot(SZrotPC(2,1),SZrotPC(1,1),'w<','markersize',3);%sample points
% plot(SZrotPC(2,2),SZrotPC(1,2),'kx');%sample points
% plot(SZrotPC(2,2),SZrotPC(1,2),'wx','markersize',3);%sample points
plot(SZrotPC(2,3),SZrotPC(1,3),'k.','markersize',15);%sample points
plot(SZrotPC(2,3),SZrotPC(1,3),'w.','markersize',5);%sample points
% plot(SZrotPC(2,4),SZrotPC(1,4),'kd');%sample points
% plot(SZrotPC(2,4),SZrotPC(1,4),'wd','markersize',3);%sample points
% plot(SZrotPC(2,5),SZrotPC(1,5),'k>');%sample points
% plot(SZrotPC(2,5),SZrotPC(1,5),'w>','markersize',3);%sample points
% plot(SZrotPT(2,1),SZrotPT(1,1),'c<');%tip points
% plot(SZrotPT(2,2),SZrotPT(1,2),'m<');%tip points
% plot(SZrotPT(2,3),SZrotPT(1,3),'m>');%tip points
% plot(SZrotPT(2,4),SZrotPT(1,4),'c>');%tip points


Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
axis(SZzoomaxis); axis off

ax3 = subplot(3,1,2);
t=hgtransform;
contourf(x,y,log10(2*Eii),[-20:0.2:-2],'edgecolor','none','Parent',t); caxis([-8 -3]); 
colorbar; colormap(ax3, bluewhitered(log10(1.7e-5))); %for different colorbar, add in octave for 528: 2.1e-5 for 511: 1.67e-5
title('log strain rate $\dot{\epsilon}$ $[s^{-1}]$','FontSize',FSS)
axis equal; shading interp;
%outline of phases
hold on; contour(x,y,RHO,[1.1 1.1],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[1.1 1.1],'--','edgecolor','w','linewidth',LWW,'Parent',t);
hold on; contour(x,y,RHO,[6 6],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[6 6],'--','edgecolor','w','linewidth',LWW,'Parent',t);
Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
axis(SZzoomaxis); axis off

ax4 = subplot(3,1,3);
t=hgtransform;
contourf(x,y,2*EiiSUM,[0:0.5:15],'edgecolor','none','Parent',t); caxis([0 10]); 
colorbar; colormap(ax4, bluewhitered(4.4));
title('SUM of inc. strain $\epsilon$','FontSize',FSS)
axis equal; shading interp;
%outline of phases
hold on; contour(x,y,RHO,[1.1 1.1],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[1.1 1.1],'--','edgecolor','w','linewidth',LWW,'Parent',t);
hold on; contour(x,y,RHO,[6 6],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[6 6],'--','edgecolor','w','linewidth',LWW,'Parent',t);
Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
axis(SZzoomaxis); axis off
drawnow

set(gcf,'Units','centimeters');
set(gcf, 'Position', PaperSizeZoom); %size of paper
pos = get(gcf,'Position');
set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[pos(3), pos(4)])
if savefigs==1
    print(gcf,sprintf('sig-E-zoom_%s_%d.png', ModName, sc),'-dpng', '-r600')%, 
end


%% -- %%% P-Sig rate, Visc, mat
gcf=figure(280+sc); clf

tangle = SZang*pi/180;

ax1 = subplot(3,1,1);
t=hgtransform;
contourf(x,y,(2*SIGii)./P,[0:0.002:0.2],'edgecolor','none','Parent',t); caxis([0 0.2])
colorbar; colormap(ax1, viridis)
title('$\Delta\sigma$ / $P$ ','FontSize',FSS)
axis equal; shading interp;
%outline of phases
hold on; contour(x,y,RHO,[1.1 1.1],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[1.1 1.1],'--','edgecolor','w','linewidth',LWW,'Parent',t);
hold on; contour(x,y,RHO,[6 6],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[6 6],'--','edgecolor','w','linewidth',LWW,'Parent',t);
Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
axis(SZzoomaxis); axis off

ax2 = subplot(3,1,2);
t=hgtransform;
contourf(x,y,log10(ETA),[8:0.1:20],'edgecolor','none','Parent',t); caxis([12 16]);%[(floor(log10(ETAminmax(sc,1))) ceil(log10(ETAminmax(sc,2)))]
colorbar; colormap(ax2, viridis);
title('$log_{10}$ Viscosity $\eta$ $[Pa s]$','FontSize',FSS)
axis equal; shading interp;
%outline of phases
hold on; contour(x,y,RHO,[1.1 1.1],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[1.1 1.1],'--','edgecolor','w','linewidth',LWW,'Parent',t);
hold on; contour(x,y,RHO,[6 6],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[6 6],'--','edgecolor','w','linewidth',LWW,'Parent',t);
Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
axis(SZzoomaxis); axis off

ax3 = subplot(3,1,3);
t=hgtransform;
contourf(x,y,(RHO),'edgecolor','none','Parent',t);
colorbar; colormap(ax3,'default')
title('Materials','FontSize',FSS)
axis equal; shading interp;
%outline of phases
hold on; contour(x,y,RHO,[1.1 1.1],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[1.1 1.1],'--','edgecolor','w','linewidth',LWW,'Parent',t);
hold on; contour(x,y,RHO,[6 6],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[6 6],'--','edgecolor','w','linewidth',LWW,'Parent',t);
Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
axis(SZzoomaxis); axis off


set(gcf,'Units','centimeters');
set(gcf, 'Position', PaperSizeZoom); %size of paper
pos = get(gcf,'Position');
set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[pos(3), pos(4)])
if savefigs==1
    print(gcf,sprintf('ETA-zoom_%s_%d.png', ModName, sc),'-dpng', '-r600')
end


%% -- %%% Work / work-rate
gcf=figure(240+sc); clf

tangle = SZang*pi/180;

ax1 = subplot(3,1,1);
t=hgtransform;
contourf(x,y,log10(Wdot),[-9:0.1:4],'edgecolor','none','Parent',t); caxis([-3 4])
colorbar; colormap(ax1, viridis)
title('$log_{10}$ work rate $\dot{W}$ [$W/m^3$]','FontSize',FSS)
axis equal; shading interp;
%outline of phases
hold on; contour(x,y,RHO,[1.1 1.1],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[1.1 1.1],'--','edgecolor','w','linewidth',LWW,'Parent',t);
hold on; contour(x,y,RHO,[6 6],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[6 6],'--','edgecolor','w','linewidth',LWW,'Parent',t);
Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
axis(SZzoomaxis); axis off

ax2 = subplot(3,1,2);
t=hgtransform;
contourf(x,y,log10(WSUM),[5:0.1:9],'edgecolor','none','Parent',t); caxis([6 9])
colorbar; colormap(ax2, viridis)
title('$log_{10}$ integral of work ${W}$ [$J/m^3$]','FontSize',FSS)
axis equal; shading interp;
%outline of phases
hold on; contour(x,y,RHO,[1.1 1.1],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[1.1 1.1],'--','edgecolor','w','linewidth',LWW,'Parent',t);
hold on; contour(x,y,RHO,[6 6],'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,[6 6],'--','edgecolor','w','linewidth',LWW,'Parent',t);
Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
axis(SZzoomaxis); axis off


set(gcf,'Units','centimeters');
set(gcf, 'Position', PaperSizeZoom); %size of paper
pos = get(gcf,'Position');
set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[pos(3), pos(4)])
if savefigs==1
    print(gcf,sprintf('W-zoom_%s_%d.png', ModName, sc),'-dpng', '-r600')
end


end
%%   ***   ***  ***
%P_distribution histogram
gcf=figure(300+sc); clf
histogram(Pdummylin/PSscale,'Normalization','probability')
%caxis([1.45e9, 1.6e9])%keep constant scaling on P colorbar
hold on
text(P_probe_conf(sc)/PSscale,0.04,'\begin{tabular}{c} $\sigma_3$ \\ $\downarrow$ \end{tabular}',...
    'Color','r','HorizontalAlignment','center','FontSize',FSS)
text((P_probe_conf(sc)*2/3+EXPsig1_acf(sc)/3)/PSscale,0.045,'\begin{tabular}{c} $P_{mean}$ \\ $\downarrow$ \end{tabular}',...
    'Color','r','HorizontalAlignment','center','FontSize',FSS)
text(EXPsig1_acf(sc)/PSscale,0.04,'\begin{tabular}{c} $\sigma_1$ \\ $\downarrow$ \end{tabular}',...
    'Color','r','HorizontalAlignment','center','FontSize',FSS)
text(1.55,0.07,'\begin{tabular}{c} $Cal \rightarrow Arg$ \\ \vline \end{tabular}',...
    'Color','b','HorizontalAlignment','center','FontSize',FSS)

%plot(rSZxB*dx*xyscale,bSZyRb*dy*xyscale,'oy');corner bottom piston
axis tight; 
xlabel('Pressure[GPa]','FontSize',FSS); 
ylabel('relative proprtion of occurences','FontSize',FSS)
title(sprintf('timestep: %d, Pressure variation in deforming Calcite $[GPa]$',sc),'FontSize',FSS)

set(gcf,'Units','centimeters');
set(gcf, 'Position', [0 0 18.3 12.3]);
pos = get(gcf,'Position');
set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[pos(3), pos(4)])
if savefigs==1
    print(gcf,sprintf('P-hist_%s_%d.pdf', ModName,sc),'-dpdf', '-r300')
end
       
% 
% %sig3_distribution histogram
% gcf=figure(330+sc); clf
% histogram(SIG3magdummylin/PSscale,'Normalization','probability')
% %caxis([1.45e9, 1.6e9])%keep constant scaling on P colorbar
% hold on
% text(P_probe_conf(sc)/PSscale,0.04,'\begin{tabular}{c} $\sigma_3$ \\ $\downarrow$ \end{tabular}',...
%     'Color','r','HorizontalAlignment','center','FontSize',FSS)
% text((P_probe_conf(sc)*2/3-EXPsig1(sc)/3)/PSscale,0.03,'\begin{tabular}{c} $P_{mean}$ \\ $\downarrow$ \end{tabular}',...
%     'Color','r','HorizontalAlignment','center','FontSize',FSS)
% text(-EXPsig1(sc)/PSscale,0.04,'\begin{tabular}{c} $\sigma_1$ \\ $\downarrow$ \end{tabular}',...
%     'Color','r','HorizontalAlignment','center','FontSize',FSS)
% text(1.55,0.06,'\begin{tabular}{c} $Cal \rightarrow Arg$ \\ \vline \end{tabular}',...
%     'Color','b','HorizontalAlignment','center','FontSize',FSS)
% 
% %plot(rSZxB*dx*xyscale,bSZyRb*dy*xyscale,'oy');corner bottom piston
% shading flat; axis tight; 
% xlabel('Pressure[GPa]','FontSize',FSS); 
% ylabel('relative proprtion of occurences','FontSize',FSS)
% title(sprintf('timestep: %d, $\\sigma_{3}$ variation in deforming Calcite $[GPa]$',sc),'FontSize',FSS)
% 
% set(gcf,'Units','centimeters');
% set(gcf, 'Position', [0 0 29.7 21]);
% pos = get(gcf,'Position');
% set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[pos(3), pos(4)])
% print(gcf,sprintf('sig3-hist_%s_%d.pdf', ModName, sc),'-dpdf', '-r300')
% 
% %sig1_distribution histogram
% gcf=figure(360+sc); clf
% histogram(SIG1magdummylin/PSscale,'Normalization','probability')
% %caxis([1.45e9, 1.6e9])%keep constant scaling on P colorbar
% hold on
% text(P_probe_conf(sc)/PSscale,0.04,'\begin{tabular}{c} $\sigma_3$ \\ $\downarrow$ \end{tabular}',...
%     'Color','r','HorizontalAlignment','center','FontSize',FSS)
% text((P_probe_conf(sc)*2/3-EXPsig1(sc)/3)/PSscale,0.03,'\begin{tabular}{c} $P_{mean}$ \\ $\downarrow$ \end{tabular}',...
%     'Color','r','HorizontalAlignment','center','FontSize',FSS)
% text(-EXPsig1(sc)/PSscale,0.04,'\begin{tabular}{c} $\sigma_1$ \\ $\downarrow$ \end{tabular}',...
%     'Color','r','HorizontalAlignment','center','FontSize',FSS)
% text(1.55,0.06,'\begin{tabular}{c} $Cal \rightarrow Arg$ \\ \vline \end{tabular}',...
%     'Color','b','HorizontalAlignment','center','FontSize',FSS)
% 
% %plot(rSZxB*dx*xyscale,bSZyRb*dy*xyscale,'oy');corner bottom piston
% shading flat; axis tight; 
% xlabel('Pressure[GPa]','FontSize',FSS); 
% ylabel('relative proprtion of occurences','FontSize',FSS)
% title(sprintf('timestep: %d, $\\sigma_1$ variation in deforming Calcite $[GPa]$',sc),'FontSize',FSS)
% 
% set(gcf,'Units','centimeters');
% set(gcf, 'Position', [0 0 29.7 21]);
% pos = get(gcf,'Position');
% set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[pos(3), pos(4)])
% print(gcf,sprintf('sig1-hist_%s_%d.pdf', ModName, sc),'-dpdf', '-r300')

    %% if Single Figures are plotted loop breaking and returning to start
    continue%break for loop here and return to start
end
        
       
%-------%%% plotting multi Figs
        %clf%not using figure(), to get rig of taking over every time the screen
        %% overview Figures
    
    if((mod(sc-1,plotintv)==0)||sc==1)%plot every n'th plt (plotintv defined at beginning)
        if plotsingle==0
        h1 = figure(1); clf;
        %set(h1,'visible','on');
        FIGps=get(gcf,'Position'); %needed for keeping resolution exactly constant fo AVI making
        set(gcf, 'Position', FIGps);%or 
        %set(gcf, 'Position', get(0,'Screensize'));%or as above 
        %colormap(viridis); %
        colormap('jet'); %for different colorbar, add in octave
                
        subplot(3,5,1)
        pcolor(x,y,RHO)
        shading flat; hold on
        plot(MCNx*dx*xyscale,MCNy*dy*xyscale,'.m')
        if SZang~=0
        plot([MCNx(7) MCNx(6)]*dx*xyscale,[MC5ay MC8ay]*dy*xyscale,'.r')
        plot(probe_Tx*dx*xyscale,probe_Ty*dy*xyscale,'*g');
        end
        plot([lAlx:rAlx]*dx*xyscale, iBC*dy*xyscale, '.r');%highlight internal bc zone
        quiver(x(1:qs:end),y(1:qs:end),Vxq,Vyq)
        axis image; axis equal; colorbar
        xlabel('diameter [mm]'); ylabel('height [mm]')
        title('Density $\rho$ $[kg/m^3]$')

        subplot(3,5,2)
        pcolor(x,y,Vx)
        shading flat; axis image; axis equal; colorbar
        xlabel('diameter [mm]'); ylabel('height [mm]')
        title('$v_x$ $[m/s]$')

        subplot (3,5,3)
        pcolor(x,y,Vy)
        shading flat; axis image; axis equal; colorbar
        xlabel('diameter [mm]'); ylabel('height [mm]')
        title('$v_y$ $[m/s]$')
        
        subplot(3,5,4)
        pcolor(x,y,SIGxxp)
        shading flat; axis image; axis equal; colorbar
        xlabel('diameter [mm]'); ylabel('height [mm]')
        title('$\sigma^\prime_{xx}$ $[Pa]$')

        subplot(3,5,5)
        pcolor(x,y,SIGyyp)
        shading flat; axis image; axis equal; colorbar
        xlabel('diameter [mm]'); ylabel('height [mm]')
        title('$\sigma^\prime_{yy}$ $[Pa]$')

        
        subplot (3,5,6)
        pcolor(x(2:end),y,P(:, 2:end))
        %caxis([1.25e9, 2.1e9])%keep constant scaling on P colorbar
        hold on
        %EDIT2
        plot([lAlx:rAlx]*dx*xyscale, cPtBC*dy*xyscale, '.r');%highlight force piston measurement zone
        plot(minPx/scX,minPy/scY,'oc'), plot(maxPx/scX,maxPy/scY,'om')
        if elliONOFF==0
        plot(probe_Cx(1)*dx*xyscale,probe_Cy(1)*dy*xyscale,'*b');
        plot(probe_Cx(2)*dx*xyscale,probe_Cy(2)*dy*xyscale,'*r');
        plot(probe_Cx(3)*dx*xyscale,probe_Cy(3)*dy*xyscale,'*b');
        plot(probe_Cx(4)*dx*xyscale,probe_Cy(4)*dy*xyscale,'*r');
        plot(probe_Cx(5)*dx*xyscale,probe_Cy(5)*dy*xyscale,'*g');
        if SZang == 0
        plot(probe_Cx(6)*dx*xyscale,probe_Cy(6)*dy*xyscale,'*r');
        else
        plot(probe_Tx*dx*xyscale,probe_Ty*dy*xyscale,'*r');
        end
        end
        %plot(rSZxB*dx*xyscale,bSZyRb*dy*xyscale,'oy');corner bottom piston
        shading flat; axis image; axis equal; colorbar
        xlabel('diameter [mm]'); ylabel('height [mm]')
        title('Pressure $P$ $[Pa]$')
        %colormap('viridis')
        %%%clf; plot(P(200,:)), axis([1 251 1.6e9 2e9]); grid on ; grid minor

        ETA(ETA<1)=1;%%%ignoring all ETA<1, not perfect but got some strange results at contact in first cycle only
        subplot(3,5,7)
        pcolor(x,y,log10(ETA))
        shading flat; hold on
        %caxis([12.3 13])
        %quiver(x(1:qs:end),y(1:qs:end),Vxq,Vyq)
        axis image; axis equal; colorbar
        xlabel('diameter [mm]'); ylabel('height [mm]')
        title('log10 Viscosity $\eta$ $[Pa \cdot s]$')

        subplot(3,5,8)
        pcolor(x,y,log10(Eii))
        shading flat;
        axis image; axis equal; colorbar;
        %caxis([-6 -4.5])
        xlabel('diameter [mm]'); ylabel('height [mm]')
        title('(tot.) $2^{nd}$strain rate invariant $\dot{\epsilon}_{II}$')

        subplot(3,5,9)
        pcolor(x,y,(SIGii))
        %caxis([6 8])%[20e6 70e6])
        hold on
        plot(probe_Tx*dx*xyscale,probe_Ty*dy*xyscale,'*k')
        shading flat; axis image; axis equal; colorbar
        xlabel('diameter [mm]'); ylabel('height [mm]')
        title('(dev.) $2^{nd}$stress invariant ${\sigma}^\prime_{II}$')  
        
        subplot (3,5,10)
        pcolor(x,y,SIGxy)
        shading flat; axis image; axis equal; colorbar
        xlabel('diameter [mm]'); ylabel('height [mm]')
        title('$\sigma_{xy}$ $[Pa]$')
             
        subplot(3,5,11)
        %extract only values for calcite (RHO=2)
        SIG1degSZ=SIG1degN;
        %SIG1degSZ(RHO~=RHO_sz)=NaN;
        pcolor(x,y,SIG1degSZ)
        %caxis([80 100])
        shading flat; hold on;
        plot([lAlx:rAlx]*dx*xyscale, iBC*dy*xyscale, '.r');%highlight force piston measurement zone
        axis image; axis equal; colorbar
        xlabel('diameter [mm]'); ylabel('height [mm]')
        title('$\sigma_1$ angle [$^\circ$]')
        colormap('jet')
        
        subplot(3,5,12)
        pcolor(x,y,SIG1mag)
        shading flat; 
        axis image; axis equal; colorbar
        %caxis([1.4e9 1.6e9])
        xlabel('diameter [mm]'); ylabel('height [mm]')
        title('$\sigma_1$ magnitude $[Pa]$')
        
        subplot(3,5,13)
        pcolor(x,y,EiiSUM)
        shading flat;
        axis image; axis equal; colorbar;
        %caxis([0 2])%sz: 1-3
        xlabel('diameter [mm]'); ylabel('height [mm]')
        title('SUM of inc. strain $e$')

        subplot(3,5,14)
        pcolor(x,y,-SIGnTAU(:,:,1,1))
        shading flat; axis image; axis equal; colorbar
        %caxis([1.45e9 1.75e9])
        xlabel('diameter [mm]'); ylabel('height [mm]')
        title('$\sigma_{n}$ in SZ orient. $[Pa]$')
        
        subplot(3,5,15)
        pcolor(x,y,-SIGnTAU(:,:,1,2))
        shading flat; axis image; axis equal; colorbar
        %caxis([20e6 70e6])
        xlabel('diameter [mm]'); ylabel('height [mm]')
        title('$\tau$ in SZ orient. $[Pa]$')
        
        drawnow
        if savefigs==1
            print(strcat(file,'_full_',num2str(sc),'.jpg'),'-djpeg')
        end
%%
        %% zoom figures
        h2 = figure(2); clf;
        %set(h2,'visible','off');
        %colormap(viridis); %
        colormap('jet'); %for different colorbar, add in octave
        pltsc = savecycle; %switch between  sc (new scale each cycle) or savecycle (scale from last cycle for all) 
        %sc gives adapted scale per step, savecycle gives same scale for all plots (based on last one)
        
        % for connected colorbar: https://stackoverflow.com/questions/8073113/matlab-how-to-assign-different-colormaps-colorbars-to-different-surfaces-in-the
        %Ptrans = 1.47e9;
        %P_Cax_Ptrans=(Ptrans-P_Cax(1))/(P_Cax(2)-P_Cax(1));
        %colormap(axp,[copper(round(P_Cax_Ptrans*128));flipud(viridis(128-round(P_Cax_Ptrans*128)))])
        scR = 5e7; %scaling factor for rounding
        if globalPscale == 1%manual definition of c-axis scale
            %P_Cax   =   [scR*round(0.92*Pconf/scR) scR*round(1.08*Pconf/scR)];
            P_Cax   =   [-1e9 7e9];
            s1_Cax  =   [-1e9 7e9];%P_Cax;
            s3_Cax  =   [-1e9 7e9];%P_Cax;
        else
            P_Cax   =   [scR*round(Pminmax(pltsc,1)/scR) scR*round(Pminmax(pltsc,2)/scR)];
            s1_Cax  =   [scR*round(SIG1magminmax(pltsc,1)/scR) scR*round(SIG1magminmax(pltsc,2)/scR)];
            s3_Cax  =   [scR*round(SIG3magminmax(pltsc,1)/scR) scR*round(SIG3magminmax(pltsc,2)/scR)];          
        end
        %1e8*round(max(Pminmax(pltsc,1),0.8*Pconf)/1e8) 1e8*round(Pminmax(pltsc,2)/1e8)
            
       
        ax1 = subplot(3,3,1);
        t=hgtransform;
        contourf(x,y,(P)/PSscale,0.5+[0:pCscf*0.05:pCscf*7]*1e9/PSscale,'edgecolor','none','Parent',t); caxis(P_acf(sc)/PSscale+CzoomIntv/PSscale);
        %c1 = colorbar; c1.Ticks = [c1.Ticks(1) round(P_nom(sc)/PSscale,2)  c1.Ticks(end)];
        colorbar; colormap(ax1, bluewhitered(P_acf(sc)/PSscale));
        title('$P$ $[GPa]$','FontSize',FSS)
        axis equal; shading interp;
        %outline of phases
        hold on; contour(x,y,RHO,2,'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,2,'--','edgecolor','w','linewidth',LWW,'Parent',t);
        plot(SZrotPC(2,:),SZrotPC(1,:),'k*');%sample points
        %xlabel('width [mm]','FontSize',FSS); ylabel('height [mm]','FontSize',FSS)
        Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
        axis(SZzoomaxis); 

        ax2 = subplot(3,3,2);
        t=hgtransform;
        contourf(x,y,(SIG1mag)/PSscale,0.5+[0:pCscf*0.05:pCscf*7]*1e9/PSscale,'edgecolor','none','Parent',t); caxis(EXPsig1_acf(sc)/PSscale+CzoomIntv/PSscale); 
        %c2 = colorbar; c2.Ticks = [c2.Ticks(1) round(EXPsig1_acf(sc)/PSscale,2)  c2.Ticks(end)];
        colorbar; colormap(ax2,bluewhitered(EXPsig1_acf(sc)/PSscale));
        title('$\sigma_1$ $[GPa]$','FontSize',FSS)
        axis equal; shading interp;
        %outline of phases
        hold on; contour(x,y,RHO,2,'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,2,'--','edgecolor','w','linewidth',LWW,'Parent',t);
        %xlabel('width [mm]','FontSize',FSS); ylabel('height [mm]','FontSize',FSS)
        Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
        axis(SZzoomaxis);

        %         subplot(3,3,3)
        %         pcolor(x,y,SIG1degN)

        ax3 = subplot(3,3,3);
        t=hgtransform;
        contourf(x,y,log10(ETA),[11:0.01:17],'edgecolor','none','Parent',t); caxis([floor(log10(ETAminmax(sc,1))) ceil(log10(ETAminmax(sc,2)))]);
        colorbar; colormap(ax3, jet);
        title('$log_{10}$ Viscosity $\eta$ $[Pa s]$','FontSize',FSS)
        axis equal; shading interp;
        %outline of phases
        hold on; contour(x,y,RHO,2,'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,2,'--','edgecolor','w','linewidth',LWW,'Parent',t);
        Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
        axis(SZzoomaxis);

        ax4 = subplot(3,3,4);
        t=hgtransform;
        contourf(x,y,(SIG3mag)/PSscale,0.5+[0:pCscf*0.05:pCscf*7]*1e9/PSscale,'edgecolor','none','Parent',t); caxis(P_probe_conf(sc)/PSscale+CzoomIntv/PSscale); 
        %c4 = colorbar; c4.Ticks = [c4.Ticks(1) round(P_probe_conf/PSscale,2)  c4.Ticks(end)];
        colorbar; colormap(ax4, bluewhitered(P_probe_conf(sc)/PSscale)); 
        title('$\sigma_3$ $[GPa]$','FontSize',FSS)
        axis equal; shading interp;
        %outline of phases
        hold on; contour(x,y,RHO,2,'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,2,'--','edgecolor','w','linewidth',LWW,'Parent',t);
        %xlabel('width [mm]','FontSize',FSS); ylabel('height [mm]','FontSize',FSS)
        Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
        axis(SZzoomaxis);

        ax5 = subplot(3,3,5);
        t=hgtransform;
        contourf(x,y,(sigSF*SIGii)/1e6,[0:pCscf*50:pCscf*4000],'edgecolor','none','Parent',t);caxis((EXPsig1_acf(sc) - P_probe_conf(sc))/1e6 +CzoomIntv/1e6); 
        colorbar; colormap(ax5, bluewhitered((EXPsig1_acf(sc) - P_probe_conf(sc))/1e6)); %for different colorbar, add in octave
        title('$\Delta\sigma$ $[MPa]$- reference (white): $\sigma_{yy}$ no overlap corr','FontSize',FSS)
        axis equal; shading interp;
        %outline of phases
        hold on; contour(x,y,RHO,2,'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,2,'--','edgecolor','w','linewidth',LWW,'Parent',t);
        Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
        axis(SZzoomaxis);

        ax6 = subplot(3,3,6);
        t=hgtransform;
        contourf(x,y,log10(gamSF*Eii),[-10:0.05:-2],'edgecolor','none','Parent',t); caxis([-8 -4]); 
        colorbar; colormap(ax6, bluewhitered(-6)); %for different colorbar, add in octave
        title('log strain rate $\dot{\epsilon}$ $[s^{-1}]$','FontSize',FSS)
        axis equal; shading interp;
        %outline of phases
        hold on; contour(x,y,RHO,2,'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,2,'--','edgecolor','w','linewidth',LWW,'Parent',t);
        Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
        axis(SZzoomaxis);

       
        subplot(3,3,7)
        plot(scplots(Pminmax(:,1)~=0,1),Pminmax(Pminmax(:,1)~=0,1)/PSscale,'c',...
            scplots(Pminmax(:,1)~=0,1),Pminmax(Pminmax(:,2)~=0,2)/PSscale,'m')
        hold on
        plot(scplots(Pminmax(:,1)~=0,1),SIG1magminmax(Pminmax(:,1)~=0,1)/PSscale,'c.-.',...
            scplots(Pminmax(:,1)~=0,1),SIG1magminmax(Pminmax(:,2)~=0,2)/PSscale,'m.-.')
        plot(scplots(Pminmax(:,1)~=0,1),SIG3magminmax(Pminmax(:,1)~=0,1)/PSscale,'c:o',...
            scplots(Pminmax(:,1)~=0,1),SIG3magminmax(Pminmax(:,2)~=0,2)/PSscale,'m:o')
        legend('P_{min}','P_{max}','{\sigma}_{1-min}','{\sigma}_{1-max}','{\sigma}_{3-min}',...
            '{\sigma}_{3-max}','Location','bestoutside')
        if elliONOFF ==0
            if SZang == 0
                 plot(scplots(Pminmax(:,1)~=0,1),P_probe_Cp(Pminmax(:,1)~=0,6)/PSscale,'g',...
                    scplots(Pminmax(:,1)~=0,1),P_probe_Cp(Pminmax(:,2)~=0,1)/PSscale,'b')
            else  
                plot(scplots(Pminmax(:,1)~=0,1),P_probe_Cp(Pminmax(:,1)~=0,3)/PSscale,'g',...
                    scplots(Pminmax(:,1)~=0,1),P_probe_Cp(Pminmax(:,2)~=0,1)/PSscale,'b')
            end
            legend('P_{min}','P_{max}','{\sigma}_{1-min}','{\sigma}_{1-max}','{\sigma}_{3-min}',...
                '{\sigma}_{3-max}','P_{center}','P_{diagNW}','Location','bestoutside')
        end
        plot(scplots(Pminmax(:,1)~=0,1),P_probe_conf(Pminmax(:,1)~=0)/PSscale,':g', 'linewidth',2)
        plot(scplots(Pminmax(:,1)~=0,1),EXPsig1(Pminmax(:,1)~=0)/PSscale,':r', 'linewidth',2)
        plot(scplots(Pminmax(:,1)~=0,1),EXPsig1_acf(Pminmax(:,1)~=0)/PSscale,'--r', 'linewidth',2)

        %text(0.2,1.8,'test')
        xlabel('calculation cycles'); ylabel('Force [GPa]')
        title(sprintf('$%s$; cycle: %d; time: %.1f [h]; piston disp: %.2f [mm]',filestr, sc, exptimePL(sc)/3600, -((uPtYV(sc)-uPtYV(1))/scY)))%text(file,'interpreter','none')
        xlim([0 savecycle+2])
        grid ON

        ax8 = subplot(3,3,8);
        t=hgtransform;
        contourf(x,y,(2*SIGii)./P,[0:pCscf*0.005:pCscf*1],'edgecolor','none','Parent',t); caxis([0 pCscf*0.5])
        colorbar; colormap(ax8, jet)
        title('$\Delta\sigma$ / $P$ ','FontSize',FSS)
        axis equal; shading interp;
        %outline of phases
        hold on; contour(x,y,RHO,2,'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,2,'--','edgecolor','w','linewidth',LWW,'Parent',t);
        Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
        axis(SZzoomaxis);
        
        ax9 = subplot(3,3,9);
        t=hgtransform;
        contourf(x,y,2*EiiSUM,[0:0.1:20],'edgecolor','none','Parent',t); caxis([0 round(2.5*abs(ginc(end)))]); 
        colorbar; colormap(ax9, bluewhitered(round(abs(ginc(end)),2)));
        title('SUM of inc. strain $\epsilon$','FontSize',FSS)
        axis equal; shading interp;
        %outline of phases
        hold on; contour(x,y,RHO,2,'edgecolor','k','linewidth',LWW,'Parent',t); contour(x,y,RHO,2,'--','edgecolor','w','linewidth',LWW,'Parent',t);
        Rz = makehgtform('zrotate',tangle); t.Matrix = Rz;
        axis(SZzoomaxis);

        
        drawnow
        if savefigs==1
            print(strcat(file,'_zoom_',num2str(sc),'.jpg'),'-djpeg')
        end
            
% save figure
% print(figure(1), '-dtiff', '0419-61.tiff',['-r',num2str(600)])
       elseif plotsingle == 1%for plotsinge 
       %% stress-time graph plotting
            

            %if sc < max(savecycle)%plot stress-time plot until last calculation is done
                figure(4); clf
                clf
                %plot(exptimePL/3600,EXPdevYY,'*')
                hold on
                plot(exptimePL/3600,EXPsig1,'*r')
                plot(exptimePL/3600,P_probe_conf,'+g')
                plot(exptimePL/3600+1,P_probe_Cp,'o')
                legend('{\sigma}^‎\prime_{yy}','{\sigma}_{yy} = {\sigma}^‎\prime_{yy}-P')
                xlabel('exp time [h]'); ylabel('Force [Pa]')
                %ylim([1 2])
                %xlim([0 max(exptimePL/3600)])
                title('Stress-time curve')
                if savefigs==1
                    print(strcat(file,'_','stress-time-curve.jpg'),'-djpeg')
                end
            %else 
            %   close()
            %end
        end
        
        % text output min-max-mean values
        fprintf('exp. time[min]: %.2f --- cycle: %d of %d \n'...
            , exptimePL(sc)/60, sc, savecycle); 
        fprintf('mode strain rate: %g --- mode viscosity: %g  --- mean strain rate: %g --- mean viscosity: %g \n'...
            , ModeMeanEiiSZ(sc,1), ModeMeanETASZ(sc,1) , ModeMeanEiiSZ(sc,2), ModeMeanETASZ(sc,2) );
        fprintf('min strain rate: %g --- max strain rate: %g --- min viscosity: %g --- max viscosity: %g \n'...
            , Eiiminmax(sc,1), Eiiminmax(sc,2), ETAminmax(sc,1), ETAminmax(sc,2) ); 
        
    end %plotintv loop not plotting every iteration
        
        %% stress-strain graph plotting
        if sc == max(savecycle)

            %acf calculationt moved to front

            ixE=(EXPsig1==0);  % get the logical vector of zero locations
            P_probe_conf(ixE)=[];         % remove those from that array
            EXPsig1_cr=EXPsig1;
            EXPsig1_cr(ixE)=[];  
            EXPsig1_acf_cr = EXPsig1_acf;
            EXPsig1_acf_cr(ixE)=[]; 
            if SZang == 0
                EXPsig1_acfp7_cr = EXPsig1_acfp7;
                EXPsig1_acfp7_cr(ixE)=[];
            end
            Pminmax_cr=Pminmax;
            Pminmax_cr(ixE,:)=[];
            SIG1magminmax_cr=SIG1magminmax;
            SIG1magminmax_cr(ixE,:)=[];
            SIG3magminmax_cr=SIG3magminmax;
            SIG3magminmax_cr(ixE,:)=[];
            ginc_smooth(ixE)=[];           % and from the time vector, too.
            exptimePL_cr=exptimePL;
            exptimePL_cr(ixE)=[];
            Pist_disp=-((uPtYV-uPtYV(1))/scY);
            Pist_disp(ixE)=[];

            %% %plot stress-strain plot in last cycle
            FSS = 12;
            if plotsingle ~= 3 && plotsingle ~= 5
                
             gcf=figure(3); clf;
                
            subplot(3,1,[1 2])

            plot(ginc_smooth,EXPsig1_acf_cr/PSscale,'-.r')
            hold on

            plot(ginc_smooth,EXPsig3_acf/PSscale,'-.','color',[0 0.5 0])
            %core
            plot(ginc_smooth,P_probe_Cp(:,3)/PSscale,'-b', 'linewidth',LWW2)
            plot(ginc_smooth,Sig1_probe_Cp(:,3)/PSscale,'-r')
            plot(ginc_smooth,Sig3_probe_Cp(:,3)/PSscale,'-','color',[0 0.5 0])   
            
            plot(ginc_smooth,Pminmax(:,3)/PSscale,'color', gy)
            plot(ginc_smooth,Pminmax(:,4)/PSscale,'color', gy)
            
            % Legend SZ symm
            legend('^{sa}\sigma_{1}','^{sa}\sigma_{3}', 'P core',...
                '{\sigma}_1 core', '{\sigma}_3 core','Location','bestoutside')
            
            xlim([0 max(ginc_smooth)])
            ylabel('$\sigma$ (MPa)','FontSize',FSSs)%Stress; Pressure

            subplot(3,1,[3 ])
            
            plot(ginc_smooth,(EXPsig1_acf-EXPsig3_acf)/PSscale,'-.k')
            hold on
            plot(ginc_smooth,(Sig1_probe_Cp(:,3)-Sig3_probe_Cp(:,3))/PSscale,'-k', 'linewidth',LWW2)
            plot(ginc_smooth,2*SIGiiminmax(:,3)/PSscale,'color', gy)
            plot(ginc_smooth,2*SIGiiminmax(:,4)/PSscale,'color', gy)

            
            %%% Legend SZ symm
            legend('^{sa}\Delta\sigma', '\Delta\sigma core', 'Location','bestoutside')

            
%             %%% addon asym
%             plot(ginc_smooth,(Sig1_probe_Cp(:,4)-Sig3_probe_Cp(:,4))/PSscale,'--dk')
%             plot(ginc_smooth,(Sig1_probe_Cp(:,5)-Sig3_probe_Cp(:,5))/PSscale,'-->k')
% 
%             plot(ginc_smooth,(Sig1_probe_Tp(:,4)-Sig3_probe_Tp(:,4))/PSscale,'-->c')
%             plot(ginc_smooth,(Sig1_probe_Tp(:,3)-Sig3_probe_Tp(:,3))/PSscale,'-->m')
% 
%             %%% Legend asymm
%             legend('^{sa}\Delta\sigma', '\Delta\sigma core',...
%                 '\Delta\sigma mid', '\Delta\sigma rim','\Delta\sigma rim low','\Delta\sigma rim high',... 
%                 '\Delta\sigma mid', '\Delta\sigma rim','\Delta\sigma rim low','\Delta\sigma rim high',...
%                 'Location','bestoutside')

            
            
            ylim([0 300])            
            xlim([0 max(ginc_smooth)])
            ylabel('$\Delta\sigma$ (MPa)','FontSize',FSSs)%Diff. stress
            xlabel('$\gamma$','FontSize',FSSs)    

            drawnow, pause(0.001)

            set(gcf,'Units','centimeters');
            set(gcf, 'Position', [0 0 16.5 9]);%[0 0 18.3 12.3]
            pos = get(gcf,'Position');
            set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[pos(3), pos(4)])
            print(gcf,sprintf('%s_stress-strain-curve_RANGE.pdf',file),'-dpdf', '-r300')
            
                
         %%       
            gcf=figure(32); clf;   %all readout points evolution                  
            subplot(3,1,[1 2])

            %plot(ginc_smooth,EXPsig1_cr/PSscale,'.:r', 'linewidth',LWW2,'markersize',12)
            plot(ginc_smooth,EXPsig1_acf_cr/PSscale,'-.r')
            hold on

%             plot(ginc_smooth,EXPsig1_acfp7_cr/PSscale,'o-r', 'linewidth',LWW2)
%             plot(ginc_smooth,P_probe_conf/PSscale,'.:', 'linewidth',LWW2,'color',[0 0.5 0],'markersize',12)
            plot(ginc_smooth,EXPsig3_acf/PSscale,'-.','color',[0 0.5 0])

            plot(ginc_smooth,P_probe_Cp(:,3)/PSscale,'-b', 'linewidth',LWW2)
            plot(ginc_smooth,P_probe_Cp(:,2)/PSscale,':xb')
            plot(ginc_smooth,P_probe_Cp(:,1)/PSscale,':<b')
            
            %%% symmetric
            plot(ginc_smooth,Sig1_probe_Cp(:,3)/PSscale,'-r')
            plot(ginc_smooth,Sig3_probe_Cp(:,3)/PSscale,'-','color',[0 0.5 0])    
            
            % Legend SZ symm
            legend('^{sa}\sigma_{1}','^{sa}\sigma_{3}', 'P core     .', 'P mid', 'P rim',...
                '{\sigma}_1 core', '{\sigma}_3 core','Location','bestoutside')
            
%             %%% A-symmetric
%             plot(ginc_smooth,P_probe_Cp(:,4)/PSscale,'--db')
%             plot(ginc_smooth,P_probe_Cp(:,5)/PSscale,'-->b')
%             
%             %%% Legend asymm
%             legend('^{sa}\sigma_{1}','^{sa}\sigma_{3}', 'P core        .', 'P mid', 'P rim',...
%                   'P mid', 'P rim','Location','bestoutside')
            
%             legend('bulk \sigma_1: \sigma_{yy} in piston','homog.SZ \sigma_1: \sigma_{yy} ACF corr.',...
%                 'bulk \sigma_3:  P(salt)', 'homog.SZ \sigma_3:  P(salt) ACF corr.', ...%'SZ rim(y), center(x)','SZ rim(y), 2(x)','SZ rim(y), 3(x)','SZ rim(y), right(x)', ... 
%                 'local P: rim', 'local P: mid', 'local P: center',...
%                 'local {\sigma}_1: center', 'local {\sigma}_3: center',...
%                 'Location','bestoutside')
            
            xlim([0 max(ginc_smooth)])
            ylabel('$\sigma$ (MPa)','FontSize',FSSs)%Stress; Pressure
            %title(sprintf('$%s$; Stress-strain curve', filestr),'FontSize',FSS);

            subplot(3,1,[3 ])
            %plot(ginc_smooth,(EXPsig1-Pconf)/PSscale,':m')
            
            plot(ginc_smooth,(EXPsig1_acf-EXPsig3_acf)/PSscale,'-.k')
            hold on
            plot(ginc_smooth,(Sig1_probe_Cp(:,3)-Sig3_probe_Cp(:,3))/PSscale,'-k', 'linewidth',LWW2)
            plot(ginc_smooth,(Sig1_probe_Cp(:,2)-Sig3_probe_Cp(:,2))/PSscale,':xk')
            plot(ginc_smooth,(Sig1_probe_Cp(:,1)-Sig3_probe_Cp(:,1))/PSscale,':<k')
            

            plot(ginc_smooth,(Sig1_probe_Tp(:,1)-Sig3_probe_Tp(:,1))/PSscale,':<c')
            plot(ginc_smooth,(Sig1_probe_Tp(:,2)-Sig3_probe_Tp(:,2))/PSscale,':<m')
            
            %%% Legend SZ symm
            legend('^{sa}\Delta\sigma', '\Delta\sigma core', '\Delta\sigma mid', '\Delta\sigma rim',...
               '\Delta\sigma rim low','\Delta\sigma rim high', 'Location','bestoutside')

%             %%% addon asym
%             plot(ginc_smooth,(Sig1_probe_Cp(:,4)-Sig3_probe_Cp(:,4))/PSscale,'--dk')
%             plot(ginc_smooth,(Sig1_probe_Cp(:,5)-Sig3_probe_Cp(:,5))/PSscale,'-->k')
% 
%             plot(ginc_smooth,(Sig1_probe_Tp(:,4)-Sig3_probe_Tp(:,4))/PSscale,'-->c')
%             plot(ginc_smooth,(Sig1_probe_Tp(:,3)-Sig3_probe_Tp(:,3))/PSscale,'-->m')
% 
%             %%% Legend asymm
%             legend('^{sa}\Delta\sigma', '\Delta\sigma core',...
%                 '\Delta\sigma mid', '\Delta\sigma rim','\Delta\sigma rim low','\Delta\sigma rim high',... 
%                 '\Delta\sigma mid', '\Delta\sigma rim','\Delta\sigma rim low','\Delta\sigma rim high',...
%                 'Location','bestoutside')
        
            ylim([0 300])            
            xlim([0 max(ginc_smooth)])
            ylabel('$\Delta\sigma$ (MPa)','FontSize',FSSs)%Diff. stress
            xlabel('$\gamma$','FontSize',FSSs)    

            drawnow, pause(0.001)

            set(gcf,'Units','centimeters');
            set(gcf, 'Position', [0 0 16.5 9]);%[0 0 18.3 12.3]
            pos = get(gcf,'Position');
            set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[pos(3), pos(4)])
            print(gcf,sprintf('%s_stress-strain-curve.pdf',file),'-dpdf', '-r300')
            
            
           %%  
            gcf=figure(31); clf;
                
            subplot(2,1,[1])

            plot(ginc_smooth,EXPsig1_acf_cr/PSscale,'-r', 'linewidth',LWW2)
            hold on
            plot(ginc_smooth,EXPsig1_cr/PSscale,'.:r', 'linewidth',LWW2, 'markersize',12)
%             plot(ginc_smooth,EXPsig1_acfp7_cr/PSscale,'o-r', 'linewidth',LWW2)
            %plot(ginc_smooth,P_probe_conf/PSscale,'.:', 'linewidth',LWW2, 'color',[0 0.5 0], 'markersize',12)
            plot(ginc_smooth,EXPsig3_acf/PSscale,'-', 'linewidth',LWW2, 'color',[0 0.5 0])

            legend('^{sa}\sigma_{1}','^{piston}\sigma','^{conf}P = ^{sa}\sigma_{3}', 'Location','bestoutside')
            
            ylim([1400 1800])
            xlim([0 max(ginc_smooth)])
            ylabel('$\sigma$ (MPa)','FontSize',FSSs)%Stress; Pressure
%             title(sprintf('$%s$; Stress-strain curve', filestr),'FontSize',FSSs);
%             legend('bulk \sigma_1: \sigma_{yy} in piston','homog.SZ \sigma_1: \sigma_{yy} ACF corr.',...
%                 'bulk \sigma_3:  P(salt)', 'homog.SZ \sigma_3:  P(salt) ACF corr.', ...
%                 'Location','W')


            subplot(2,1,[2])
            plot(ginc_smooth,(EXPsig1_acf-EXPsig3_acf)/PSscale,'-k', 'linewidth',LWW2)
            hold on
            plot(ginc_smooth,(EXPsig1-Pconf)/PSscale,'.:k', 'linewidth',LWW2, 'markersize',12)
            
            
            legend('^{sa}\Delta\sigma', '^{piston}\Delta\sigma     .',  'Location','bestoutside')
            
            ylim([150 250])
            xlim([0 max(ginc_smooth)])
            ylabel('$\Delta\sigma$ (MPa)','FontSize',FSSs)%Diff. stress
            xlabel('$\gamma$','FontSize',FSSs)%strain
%             legend('bulk \sigma_d: \sigma_{yy} - P(salt)', 'h.SZ \sigma_d: [\sigma_{yy} - P(salt)] ACF corr.',...
%                 'Location','NW')
            drawnow

            set(gcf,'Units','centimeters');
            set(gcf, 'Position', [0 0 16.5 8]);%[0 0 18.3 12.3]
            pos = get(gcf,'Position');
            set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[pos(3), pos(4)])
            print(gcf,sprintf('%s_BULK_stress-strain-curve.pdf',file),'-dpdf', '-r300')
            end
           %% plot for pure shear cylinder width evolution
            if plotsingle == 3 || plotsingle == 0 && SZang==0
                dSZwidthVc_cr=dSZwidthVc;
                dSZwidthVc_cr(ixE)=[];  
                dSZsupSIG1_cr=dSZsupSIG1;
                dSZsupSIG1_cr(ixE)=[];  
                dSZsupP_cr=dSZsupP;
                dSZsupP_cr(ixE)=[]; 
                dNOMsup_cr=dNOMsup;
                dNOMsup_cr(ixE)=[]; 
                %%
                figure(5); clf
                plot(ginc_smooth,dSZwidthVc_cr,'-ob')
                hold on
                plot(ginc_smooth,dSZsupSIG1_cr,'-.r')
                plot(ginc_smooth,dSZsupP_cr,'-.g')
                plot(ginc_smooth,dNOMsup_cr,'-','Color',[255 125 0]/255)
                dASYMP     =   dSZwidthVc(1)./((gSZlength/gSZlength(1)).^4/2+0.5-(gSZlength/gSZlength(1))*0.5+0.5);
                plot(ginc_smooth,dASYMP,'-*m')
                dASYMPs     =   dSZwidthVc(1)./((gSZlength/gSZlength(1)).^7/4+0.75);
                plot(ginc_smooth,dASYMPs,'--om','linewidth',LWW2)
                legend('max-d','sig1','P','cylinder','complex fit','(\Sigma\gamma)^7 fit','location','nw')
                xlabel('sum of incremental strain $\Sigma\gamma_{inc}$'); ylabel('width of cylinder [px]')
                title(sprintf('$%s$; width evolution', filestr));
                drawnow, pause(0.001)
                if savefigs==1
                    print(strcat(file,'_','SZ_width.jpg'),'-djpeg')
                end
            end

        %% HISTPLOTS strain rate and viscosity
        LcolSe = num2cell(lines((savecycle-1)/plotintv+1),2);%create (number of lines) colors for plot from 'lines' colormap

        if plotsingle == 4 || plotsingle == 0
            ylimmax = 0;
            figure(7); clf;
            figure(8); clf;
            for i=1:plotintv:savecycle %old versions calcintv wat plotintv
                figure(7);
                hold on
                [a, b] = histcounts(log10(EiiSZ(i,:)),400);
                h1(i) = plot((diff(b)+b(1:end-1)),smooth(a)/round(sum(EiiSZ(i,:),'omitnan')+1),'-','linewidth',LWW2-0.5);
                [aisz, bisz] = histcounts(log10(Eii_inSZ(i,:)),400);
                scf= max(a)/max(aisz);
                h2(i) = plot((diff(bisz)+bisz(1:end-1)),smooth(aisz)*scf/round(sum(EiiSZ(i,:),'omitnan')+1),':','linewidth',LWW2+1);
                if i == savecycle %plot legend at end only to not get warning
                    legendCell = cellstr(num2str([1:plotintv:savecycle]', 'cycle %-d'));
                    legendCellinSZ = strcat(legendCell,' in SZ');
                    legend([h1(1:plotintv:savecycle) h2(1:plotintv:savecycle)],[legendCell; legendCellinSZ])
                    set([h1(1:plotintv:savecycle) h2(1:plotintv:savecycle)], {'color'}, [LcolSe; LcolSe]);
                end
                title(sprintf('$%s$, EiiSZ cycle: %d', filestr, i))
                %axis([0 2e14 0 5000])
                drawnow


                figure(8);
                hold on
                [a, b] = histcounts(log10(ETASZ(i,:)),400);
                h3(i) = plot((diff(b)+b(1:end-1)),smooth(a)/round(sum(ETASZ(i,:),'omitnan')/1e15+1),'-','linewidth',LWW2-0.5);
                [aisz, bisz] = histcounts(log10(ETA_inSZ(i,:)),400);
                scf= max(a)/max(aisz);
                h4(i) = plot((diff(bisz)+bisz(1:end-1)),smooth(aisz)*scf/round(sum(ETASZ(i,:),'omitnan')/1e15+1),':','linewidth',LWW2+1);

                if i == savecycle %plot legend at end only to not get warning
                    legendCell = cellstr(num2str([1:plotintv:savecycle]', 'cycle %-d'));
                    legendCellinSZ = strcat(legendCell,' in SZ');
                    legend([h3(1:plotintv:savecycle) h4(1:plotintv:savecycle)],[legendCell; legendCellinSZ])
                    set([h3(1:plotintv:savecycle) h4(1:plotintv:savecycle)], {'color'}, [LcolSe; LcolSe]);
                end
                title(sprintf('$%s$, ETA cycle: %d',filestr, i))
                if i>1 %adapt scaling after first cycle to max of all following
                    ylimmax=max([a/round(sum(ETASZ(i,:),'omitnan')/1e15+1) ylimmax]);
                    ylim([0 ylimmax+0.05])
                end
            end %iteration histplots
            drawnow
            
            if savefigs==1
                figure(7)
                print(strcat(file,'_','Eii-hist-Evo.jpg'),'-djpeg')
                figure(8)
                print(strcat(file,'_','ETA-hist-Evo.jpg'),'-djpeg')
            end
        end%end plot only histplots     
        %output message on timesteps:
%         fprintf('imposed time step was: %d maximum allowed time step was: %g\n', dtmax, max(dtSAVEmax));
            
        %save plot variables for re-plotting
        savename = sprintf('P-S-e_data_%s', file);
        save(savename,'file','filestr','ginc_smooth','P_probe_Cp','P_probe_conf','EXPsig1','Sig1_probe_Cp','acfFactor','calcintv','EiiSZ','ETASZ')
        end %max savecycle
        
        if plotsingle == 5 && sc == max(scPlotBigFig)  % save stress data
            dSIGname= sprintf('dSIG_%d_%s',sc,file);
            save(dSIGname,'RHO','SIGii','Eii','Temp','x','y','SZzoomaxis', 'SZang')
        end
end


%% save run-data for re-plotting
if scPlotBigFig == 0
    if length(exptimePL_cr)==savecycle
        rundata=[exptimePL; Pist_disp; ginc_smooth'; EXPsig1; P_probe_conf; Pminmax'; SIG1magminmax'; SIG3magminmax'];
        rundataname= sprintf('rundata_%s.txt',file);
        fid = fopen(rundataname,'w');
        fprintf(fid, 'exp_time Pist_disp ginc_smooth EXPsig1 P_probe_conf Pmin Pmax S1min S1max S3min S3max\n');
        fprintf(fid, '%g %g %g %g %g %g %g %g  %g %g %g\n', rundata);
        fclose(fid);
    end
end

%reset matlab to normal
%set(0,'defaultTextInterpreter','tex')

cd ..%back to root folder when looping through files


end%ifiles
