%Plotting routine to visialize Griggs Rig modelling results 
% -> time segment extracion routine
%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;

EvCyS = [2 20 ; 11 20 ; 16 20 ; 20 20];%cycles to evaluate if only end is monitored for MAX extraction 

savefigs = 1;
plotsingle = 0; %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)

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

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

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

ModName = '495BR';%part of output filename
ifiles = dir('*v19bHR_*1mm*551x*.mat');%file to be read

ifile = ifiles';
file = ifile.name%print filename!
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 %

%%
% if savefigs ==1 %create and enter subfolder for each plot file
%     mkdir(file); cd(file)
% 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

PSscale     = 1e9;
xyscale     = 1000;


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

for i = 1:length(EvCyS)
    EvCy = EvCyS (i,1):EvCyS (i,2)

%% max value extraction and calc to nodes    

RHO         = RhoPL(:,:,EvCy(end)); RHO=round(RHO*100)/100;


%load single plot variables 

Sig1mTS     = PrmTS + 2*SigiimTS;


PrmMAX      = max(PrmTS(:,EvCy),[],2); %
Sig1mMAX    = max(Sig1mTS(:,EvCy),[],2);
SigiimMAX   = max(SigiimTS(:,EvCy),[],2);
if EvCy(1) == 1%first and last step plotting needs special treatment as its a sum
    EiimSUM     = EiisummTS(:,EvCy(end));
else
    EiimSUM     = EiisummTS(:,EvCy(end)) - EiisummTS(:,(EvCy(1)-1)) ;
end

xm      = xmTS(:,EvCy(end));%define here which shape is used
ym      = ymTS(:,EvCy(end));


% interp markers to nodes
PrMxSUM     =zeros(Ny,Nx);
Sig1MxSUM   =zeros(Ny,Nx);
WTSUM       =zeros(Ny,Nx);
SigiiMxSUM  =zeros(Ny,Nx);
EiiSmSUM    =zeros(Ny,Nx);

PrMx    =zeros(Ny,Nx);
Sig1Mx  =zeros(Ny,Nx);
SigiiMx =zeros(Ny,Nx);
EiiSm   =zeros(Ny,Nx);
 
% Interpolationg Nodal values from markers
m=1;
    
for jm=1:1:marknum

    if(xm(m)>=0 && xm(m) <= xsize && ym(m) >=0 && ym(m) <=ysize)
    %calculate nodal value (from BI-lineraly wieghted markers) AFTER first time step
    % --- i,j ----- j+1 ---
    %      |    m    |
    % -- i+1,j--- i+1,j+1 -
    j=fix((xm(m)-x(1))/dx)+1 ;%find the left node, the right node is j+1
    i=fix((ym(m)-y(1))/dy)+1;
    %to check for exact marker/node overlap - OUT OF BORDERS:
        if(j<1)
            j=1;
        elseif(j>Nx-1)
            j=Nx-1;
        end 
        if(i<1)
            i=1;
        elseif(i>Ny-1)
            i=Ny-1;
        end
    %%% calculate weight factors
    wtmij=(1-(abs(xm(m)-x(j))/dx))*(1-(abs(ym(m)-y(i))/dy));%calculate weight factor for marker top left
    wtmi1j=(1-(abs(xm(m)-x(j))/dx))*(abs(ym(m)-y(i))/dy);% -"- bottom left
    wtmij1=(abs(xm(m)-x(j))/dx)*(1-(abs(ym(m)-y(i))/dy));% -"- top right
    wtmi1j1=(abs(xm(m)-x(j))/dx)*(abs(ym(m)-y(i))/dy);% -"- bottom right
    %%% calculate sums
    %P
    PrMxSUM(i,j)=PrMxSUM(i,j)+PrmMAX(m)*wtmij;%contribution to value of top left
    PrMxSUM(i+1,j)=PrMxSUM(i+1,j)+PrmMAX(m)*wtmi1j;% -"- bottom left
    PrMxSUM(i,j+1)=PrMxSUM(i,j+1)+PrmMAX(m)*wtmij1;% -"- top right
    PrMxSUM(i+1,j+1)=PrMxSUM(i+1,j+1)+PrmMAX(m)*wtmi1j1;% -"- bottom right
    %Sig1
    Sig1MxSUM(i,j)=Sig1MxSUM(i,j)+Sig1mMAX(m)*wtmij;%contribution to value of top left
    Sig1MxSUM(i+1,j)=Sig1MxSUM(i+1,j)+Sig1mMAX(m)*wtmi1j;% -"- bottom left
    Sig1MxSUM(i,j+1)=Sig1MxSUM(i,j+1)+Sig1mMAX(m)*wtmij1;% -"- top right
    Sig1MxSUM(i+1,j+1)=Sig1MxSUM(i+1,j+1)+Sig1mMAX(m)*wtmi1j1;% -"- bottom right
    %Sigii
    SigiiMxSUM(i,j)=SigiiMxSUM(i,j)+SigiimMAX(m)*wtmij;%contribution to value of top left
    SigiiMxSUM(i+1,j)=SigiiMxSUM(i+1,j)+SigiimMAX(m)*wtmi1j;% -"- bottom left
    SigiiMxSUM(i,j+1)=SigiiMxSUM(i,j+1)+SigiimMAX(m)*wtmij1;% -"- top right
    SigiiMxSUM(i+1,j+1)=SigiiMxSUM(i+1,j+1)+SigiimMAX(m)*wtmi1j1;% -"- bottom right
    %Eiii
    EiiSmSUM(i,j)=EiiSmSUM(i,j)+EiimSUM(m)*wtmij;%contribution to value of top left
    EiiSmSUM(i+1,j)=EiiSmSUM(i+1,j)+EiimSUM(m)*wtmi1j;% -"- bottom left
    EiiSmSUM(i,j+1)=EiiSmSUM(i,j+1)+EiimSUM(m)*wtmij1;% -"- top right
    EiiSmSUM(i+1,j+1)=EiiSmSUM(i+1,j+1)+EiimSUM(m)*wtmi1j1;% -"- bottom right
    %weight
    WTSUM(i,j)=WTSUM(i,j)+wtmij;%contribution to weight of top left
    WTSUM(i+1,j)=WTSUM(i+1,j)+wtmi1j;% -"- bottom left
    WTSUM(i,j+1)=WTSUM(i,j+1)+wtmij1;% -"- top right
    WTSUM(i+1,j+1)=WTSUM(i+1,j+1)+wtmi1j1;% -"- bottom right

    %%%counter
    m=m+1;
    end
end

%calculate final RHO on nodes from Sum and Wheight factor
%tests if computuation is possible, if not, keeps old value
%https://blogs.mathworks.com/loren/2008/05/14/acting-on-specific-elements/
indWTz = find(WTSUM ~= 0);%indices of WTSUM ~= 0
PrMx(indWTz) = PrMxSUM(indWTz) ./ WTSUM(indWTz);
Sig1Mx(indWTz) = Sig1MxSUM(indWTz) ./ WTSUM(indWTz);
SigiiMx(indWTz) = SigiiMxSUM(indWTz) ./ WTSUM(indWTz);
EiiSm(indWTz) = EiiSmSUM(indWTz) ./ WTSUM(indWTz);
    

%% Plot cycle --------------------------
sc = EvCy(end);

xp = x*xyscale;%axis scaling
yp = y*xyscale;

%-------%%% plotting single Figs
FSS = 12;
FSSs = 8;
LWW = 2;
LWW2 = 1.5;%define line-width for graphs

% axis and rotation definition 
tangle = SZang*pi/180;%matrix rotation angle

%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];
SZcenterNd  = round(SZcenter ./ [xyscale*dy xyscale*dx]');
RotM = [cosd(SZang) sind(SZang);-sind(SZang) cosd(SZang)];
SZrotCenter = round(RotM*SZcenter,2);%xyscale*dx*[probe_Cy(3); probe_Cx(3)],2);%center

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.55*[-SZmaxExtent SZmaxExtent] SZrotCenter(1)+xyscale*0.55*[-SZthick SZthick]];%zoom axis definitionis
end

% plot single - full size

%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

% SZ zoom stack distribution in 2D zoom

%Ptrans = Ptrans-0.13e9;

%settings for contour plots
nEx = 50; %color intervals
contIV = 100*1e6/PSscale; %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 + [-1250e6 1250e6])/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 = -0.5;%1 - 3*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

%
P_Cax= [1.5e9 3.5e9]/PSscale;
s1_Cax = P_Cax;
ds_Cax = [0 2e9]/PSscale;
Eii_Cax = [0 (ceil(EiiSm(SZcenterNd(1),SZcenterNd(2))*2*1.2*2))/2];%scaling: 2 is Eii to gamma, 1.3 is for plotting

P_bwrC = 2.75e9/PSscale;
s1_bwrC= P_bwrC;



%% -- %%% Stresses Sig1 , P, Sig3
gcf=figure(4); clf
ax1 = subplot(2,1,1);
t=hgtransform;
contourf(xp,yp, SigiiMx*2/PSscale,(lEx*ds_Cax(1) :contIV: hEx*ds_Cax(2)),'edgecolor','none','Parent',t); 
caxis(ds_Cax); 
colorbar; colormap(viridis)%(ax1,bluewhitered(ds_bwrC));
title('$\Delta\sigma$ $[GPa]$','FontSize',FSS)
axis equal;% shading interp;
%outline of phases
hold on; contour(xp,yp,RHO,[1.1 1.1],'edgecolor','k','linewidth',LWW,'Parent',t); contour(xp,yp,RHO,[1.1 1.1],'--','edgecolor','w','linewidth',LWW,'Parent',t);
hold on; contour(xp,yp,RHO,[6 6],'edgecolor','k','linewidth',LWW,'Parent',t); contour(xp,yp,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

ax2 = subplot(2,1,2);
t=hgtransform;
contourf(xp,yp,EiiSm*2,linspace(0,Eii_Cax(2),100),'edgecolor','none','Parent',t); 
caxis(Eii_Cax);
colorbar; colormap(viridis)%(ax2, bluewhitered(P_bwrC));
title('$\gamma$','FontSize',FSS)
axis equal; %shading interp;
%outline of phases
hold on; contour(xp,yp,RHO,[1.1 1.1],'edgecolor','k','linewidth',LWW,'Parent',t); contour(xp,yp,RHO,[1.1 1.1],'--','edgecolor','w','linewidth',LWW,'Parent',t);
hold on; contour(xp,yp,RHO,[6 6],'edgecolor','k','linewidth',LWW,'Parent',t); contour(xp,yp,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

PaperSizeZoom = [0 0 30 7]

set(gcf,'Units','centimeters');
set(gcf, 'Position', PaperSizeZoom + [0 13 0 0]); %size of paper
pos = get(gcf,'Position');
set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[pos(3), pos(4)])
if savefigs==1
    print(gcf,sprintf('TS-dS_E_%s_cy%d-%d.png', ModName, EvCy(1),EvCy(end)),'-dpng', '-r600')
end

end
