%-----------------------%
%   Read JSBACH Ouput   %
%-----------------------%

clc;                % Refresh Command Window
clear;              % clear Workspace
close all;          % close all Figures
format long g       % 14 decimal places

addpath(genpath('../Functions')) % Load functions

physicalConstants


printSwitch = true; 
plotStyle   = 'diff'; % choose between 'lin', 'diff' and 'rel'

mask    = 'maskCEBa';
redID   = 'ref';
bluID   = 'ws1';
outFold = 'Plots';







%%% Read Red Data %%%

redExp  = 'CEBa3h';
redFold = 'outdata_20191104_finalThesis';
redName = 'ref_LTfalse.ref_LTtrue';

redPath = './Output/';
redFile = [redName '.GPP.1979-2008.' mask '.Apr-Sep.mat'];
redData = load([redPath '/' redFile]);

%%% LTfalse %%%
redF.gpp = double(redData.lgpp);
redF.npp = double(redData.lnpp);
redF.Cos = double(redData.lCo);

%%% LTtrue %%%
redT.gpp = double(redData.cgpp);
redT.npp = double(redData.cnpp);
redT.Cos = double(redData.cCo);



%%% Read Blue Data %%%

bluExp  = 'rs0CEBa3h';
bluFold = 'outdata_20200820_CEBa_ws1';
bluName = 'ws1_LTfalse.ws1_LTtrue';

bluPath = './Output/';
bluFile = [bluName '.GPP.1979-2008.' mask '.Apr-Sep.mat'];
bluData = load([bluPath '/' bluFile]);

%%% LTfalse %%%
bluF.gpp = double(bluData.lgpp);
bluF.npp = double(bluData.lnpp);
bluF.Cos = double(bluData.lCo);

%%% LTtrue %%%
bluT.gpp = double(bluData.cgpp);
bluT.npp = double(bluData.cnpp);
bluT.Cos = double(bluData.cCo);


minAngle = 10;
cosLimit = cos( (90 - minAngle)/180 * pi );

redIsDay = redF.Cos > cosLimit & redT.Cos > cosLimit;
bluIsDay = bluF.Cos > cosLimit & bluT.Cos > cosLimit;


c = 1e06;

redGppT = redT.gpp(redIsDay)*c;
redNppT = redT.npp(redIsDay)*c;
redGppF = redF.gpp(redIsDay)*c;
redNppF = redF.npp(redIsDay)*c;

bluGppT = bluT.gpp(bluIsDay)*c;
bluNppT = bluT.npp(bluIsDay)*c;
bluGppF = bluF.gpp(bluIsDay)*c;
bluNppF = bluF.npp(bluIsDay)*c;


x1 = redGppF;
x2 = bluGppF;

switch plotStyle
    case 'lin'
        y1 = redGppT;
        y2 = bluGppT;
    case 'diff'
        y1 = redGppT - redGppF;
        y2 = bluGppT - bluGppF;
    case 'rel'
        y1 = (redGppT-redGppF)./redGppF*100;
        y2 = (bluGppT-bluGppF)./bluGppF*100;
    otherwise
        error('Wrong plot style')
end


% Figure size
fig.xPos     = 200;
fig.yPos     = 100;
fig.width    = 1600;
fig.height   = 800;
fig.Position = [fig.xPos fig.yPos fig.width fig.height];

figure;
set(gcf,'Position',fig.Position)


cZero  = 'k';
wZero  = 2.0;
aZero  = 0.4;

cScat1 = newColor(2);
cScat2 = newColor(1);
aScat1 = 0.8;
aScat2 = 0.8;
sScat  = 100;
wScat  = 1.5;

cFit1  = 'k';
cFit2  = 'k';
aFit1  = 1.0;
aFit2  = 1.0;
sFit   = '--';
wFit   = 2.5;
nDig   = 3;

xLim   = [0 13];
xTick  = 0:2:12;
xSkip  = 2;
unit   = [ char(956) 'mol/(m^2s)'];

xLab   = ['GPP(T_{air})   [' unit ']'];




alLeg  = 'SouthEast';
asLeg1 = ' stressed';
asLeg2 = ' unstressed';


switch plotStyle
    case 'lin'
        lm = 0.055;     
        yLim  = [0 13];
        ySkip = 2;
        xZero = [-1 13];
        yZero = [-1 13];
        pFit  = [1 0];
        fFit  = 'a*x+b';
        tFit  = 'lin';
        yLab  = ['GPP(T_{leaf})   [' unit ']'];
        lLeg  = 'NorthWest';
        xtext = 11.0;
        ytext = 12.6;
    case 'diff'
        lm = 0.07; 
        yLim  = [-1.2 1.2];
        ySkip = 2;
        xZero = [-1 13];
        yZero = [ 0  0];
        pFit  = [1 1 1];
        fFit  = 'a*x^2+b*x+c';
        tFit  = 'lin';
        yLab  = ['GPP(T_{leaf}) ' char(8211) ' GPP(T_{air})   [' unit ']'];
        lLeg   = 'NorthWest';
        xtext = 12;
        ytext = 1.13;
    case 'rel'
        lm = 0.05; 
        yLim  = [-40 40];
        ySkip = 2;
        xZero = [-1 13];
        yZero = [ 0  0];
        pFit  = [1 0];
        fFit  = 'a/x+b';
        tFit  = 'hyp';
        yLab  = 'GPP [relative change in %]';
        lLeg   = 'North';
        xtext  = 12;
        ytext = 37;
    otherwise
        error('Wrong plot style')
end






% Axes Position
tm = 0.02;      % top margin
bm = 0.11;
rm = 0.01;
mm = 0.01;
yh = 1-bm-tm;
xw = (1-lm-mm-rm)/2;
xm = 1-rm-xw;

ax.Pos1  = [lm bm xw yh];    % left
ax.Pos2  = [xm bm xw yh];    % right
fontSize = 26;
fontName = 'Verdana';




%%% Red Plot (left) %%%
ax1 = subplot('Position',ax.Pos1);
s1 = scatter(x1,y1);
hold on
s1.MarkerEdgeColor = cScat1;
s1.MarkerEdgeAlpha = aScat1;
s1.LineWidth       = wScat;
s1.SizeData        = sScat;


z1           = plot(xZero,yZero);
z1.Color     = cZero;
z1.Color     = [z1.Color,aZero];
z1.LineWidth = wZero;

fit1            = fit(x1,y1,fFit,'StartPoint',pFit);
line1           = plot(fit1);
switch plotStyle
    case 'rel'
    line1.YData     = line1.YData(line1.XData>0);
    line1.XData     = line1.XData(line1.XData>0);
end
line1.Color     = cFit1;
line1.Color     = [line1.Color,aFit1];
line1.LineStyle = sFit;
line1.LineWidth = wFit;


ax1.FontSize = fontSize;
ax1.FontName = fontName;
ax1.XLim     = xLim;
ax1.YLim     = yLim;
switch plotStyle
    case 'diff'
        ax1.YTick    = -1.2:0.2:1.2;
end
ax1.XTick    = xTick;
skipXTickLabel(gca,xSkip)
skipYTickLabel(gca,ySkip)
xl1          = xlabel(xLab);
yl1          = ylabel(yLab);
xl1.FontSize = fontSize-2;
yl1.FontSize = fontSize-2;

text1 = text(xtext,ytext,'(a)');
text1.FontName = 'Verdana';
text1.FontSize = fontSize+4;

leg1 = legend(line1,'boxoff');
legend boxoff
leg1.FontSize  = fontSize;
leg1.FontName  = fontName;
switch plotStyle
    case {'lin','diff'}
        leg1.String    = linearFitLabel(fit1,nDig);
    case 'rel'
        leg1.String    = functionLabel(fit1.a ,fit1.b,nDig,tFit);
end
leg1.Location  = lLeg;

ax1a            = axes('Position',ax.Pos1);
ax1a.Visible    = 'off';         
leg1a           = legend(ax1a,s1);
legend boxoff
leg1a.FontSize  = fontSize+4;
leg1a.FontName  = fontName;
leg1a.String    = asLeg1;
leg1a.Location  = alLeg;

offsetAxesOrig(ax1)




%%% Blue Plot (right) %%%
ax2 = subplot('Position',ax.Pos2);
s2 = scatter(x2,y2);
hold on
s2.MarkerEdgeColor = cScat2;
s2.MarkerEdgeAlpha = aScat2;
s2.LineWidth       = wScat;
s2.SizeData        = sScat;

z2           = plot(xZero,yZero);
z2.Color     = cZero;
z2.Color     = [z2.Color,aZero];
z2.LineWidth = wZero;

fit2            = fit(x2,y2,fFit,'StartPoint',pFit);
line2           = plot(fit2);
switch plotStyle
    case 'rel'
    line2.YData     = line2.YData(line2.XData>0);
    line2.XData     = line2.XData(line2.XData>0);
end
line2.Color     = cFit2;
line2.Color     = [line2.Color,aFit2];
line2.LineStyle = sFit;
line2.LineWidth = wFit;


ax2.FontSize   = fontSize;
ax2.FontName   = fontName;
ax2.XLim       = xLim;
ax2.YLim       = yLim;
switch plotStyle
    case 'diff'
        ax2.YTick    = -1.2:0.2:1.2;
end
ax2.XTick      = xTick;
ax2.YTickLabel = {''}; 
skipXTickLabel(gca,xSkip)
skipYTickLabel(gca,ySkip)
xl2            = xlabel(xLab);
yl2            = ylabel('');
xl2.FontSize   = fontSize-2;

text2 = text(xtext,ytext,'(b)');
text2.FontName = 'Verdana';
text2.FontSize = fontSize+4;

leg2 = legend(ax2,line2);
legend boxoff
leg2.FontSize  = fontSize;
leg2.FontName  = fontName;
switch plotStyle
    case {'lin','diff'}
        leg2.String    = linearFitLabel(fit2,nDig);
    case 'rel'
        leg2.String    = functionLabel(fit2.a ,fit2.b,nDig,tFit);
end
leg2.Location  = lLeg;

ax2a            = axes('Position',ax.Pos2);
legend boxoff
ax2a.Visible    = 'off';         
leg2a           = legend(ax2a,s2);
leg2a.FontSize  = fontSize+4;
leg2a.FontName  = fontName;
leg2a.String    = asLeg2;
leg2a.Location  = alLeg;

offsetAxesOrig(ax2)


if(printSwitch)
   outFile = ['gpp.' plotStyle '.' redID '.' bluID '.' mask '.png'];
   drucken([outFold '/' outFile])
end












