%-----------------------%
%   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


TcasSwitch = false;

legStr      = {' stressed' ' unstressed'};
expID       = 'Ref.ws1';
colMapCol   = 'rbCo';
sameColBar  = false;
printSwitch = true; 
funcLegend  = true;
obsRegression = true;
whiteDots   = false;    perc = 2;



inputFolder   = './Output/';
inputFile   = 'trop.CEBa.Ref.rawTemp.mat';





load([inputFolder inputFile]);
TcanRef = Tcan;
if(TcasSwitch)
    TairRef = Tcas;
else
    TairRef = Tair;
end

inputFolder   = './Output/';
inputFile   = 'trop.CEBa.WatStr1.rawTemp.mat';


load([inputFolder inputFile]);

if(TcasSwitch)
    Tair = Tcas;
end

minAngle = 10;
cosLimit = cos( (90 - minAngle)/180 * pi );
isDay   = cosz > cosLimit & ismember(month(Time),[4 5 6 7 8 9]); 
isNight = cosz <= cosLimit;


% x and y labels
xLabels = {['T_{air} [' char(176) 'C]'] ...
           ['T_{air} [' char(176) 'C]'] ...
           ['T_{air} [' char(176) 'C]'] ...
            'g_{s} [mm/s]'              ...
            'g_{s} [mm/s]'              ...  
            ['T_{leaf} [' char(176) 'C]']};
yLabels = {['T_{leaf} ' char(8211) ' T_{air} [K]'] ...
           ['T_{leaf} ' char(8211) ' T_{air} [K]'] ...
           ['T_{leaf} ' char(8211) ' T_{air} [K]'] ...
           ['T_{leaf} ' char(8211) ' T_{air} [K]'] ...
           ['T_{leaf} [' char(176) 'C]']           ... 
           'g_{s} [mm/s]'                 };

% Experiment ID
if(TcasSwitch)
    expStr  = {['Trop.' expID '.Tcas.Sep2001-Oct2003']};
else
    expStr  = {['Trop.' expID '.Sep2001-Oct2003']};
end

% x and y Data; red (summer, day, ...) + blue (winter, night, ...)
yDataS  = { TcanRef(isDay)   - TairRef(isDay)};
yDataW  = { Tcan(isDay) - Tair(isDay) };
xDataS  = { TairRef(isDay)-T0};
xDataW  = { Tair(isDay)-T0}; 
% x and y limits
if(TcasSwitch)
   xLimits = {[20 34]};
else
    xLimits = {[20 34]};
end
   
yLimits = {[-5 5]};       
       
% number of bins in x direction (2d)       
numBinX =  [91];
% number of bins in y direction (2d)       
numBinY =  [90];                    

% x and ytick length
xtLen   = [2];
ytLen   = [1];

% limits of secondary y axis
y2Limits = {[0 9]};                
y4Limits = {[0 3]}; 

if(TcasSwitch)
histMaxS = [0.090];
histMaxW = [0.090];
else
histMaxS = [0.090];
histMaxW = [0.090];    
end

if(TcasSwitch)
   skipXTick = [1];
else
   skipXTick = [2];
end


skipYTick = [1 ];


idx = 0;
for iVar = 1%2:length(yDataS)
    
    if(ismember(iVar,[2 3 4]))
        nanIdx = draf.gc(isDay)<3e-9;
        xDataS{iVar}(nanIdx) = NaN;
        xDataW{iVar}(nanIdx) = NaN;
        yDataS{iVar}(nanIdx) = NaN;
        yDataW{iVar}(nanIdx) = NaN;
    end
    
    switch colMapCol
        case 'rebl'
            n       = 9;
            colMap1 = cmocean('amp');
            colMap2 = flipud(cmocean('ice'));
            colMap1 = [[1 1 1];colMap1];
            colMap2 = [[1 1 1];colMap2];
            colMap1 = colMap1(1:end-40,:);
            colMap2 = colMap2(1:end-40,:);
            colStr  = 'rebl';
            contSwi = false;
        case 'rebl2'
            n       = 9;
            colMap1 = cmocean('amp');
            colMap2 = flipud(cmocean('ice'));
            colMap1 = [[1 1 1];colMap1];
            colMap2 = [[1 1 1];colMap2];
            colMap1 = colMap1(1:end-40,:);
            colMap2 = colMap2(1:end-40,:);
            colMap1 = discColMap(colMap1,n);
            colMap2 = discColMap(colMap2,n);
            colStr  = 'rebl2';
            contSwi = false;
        case 'rbCo'
            n       = 9;
            colMap1 = cmocean('amp');
            colMap2 = flipud(cmocean('ice'));
            colMap1 = [[1 1 1];colMap1];
            colMap2 = [[1 1 1];colMap2];
            colMap1 = colMap1(1:end-40,:);
            colMap2 = colMap2(1:end-40,:);
            colMap1 = discColMap(colMap1,n);
            colMap2 = discColMap(colMap2,n);
            colStr  = 'rbCo';
            contSwi = true;
        otherwise
            error('Choose valid color map option!')
    end


    dataS = yDataS{iVar};  
    dataW = yDataW{iVar}; 
    WspdS = xDataS{iVar};   
    WspdW = xDataW{iVar};  

  
    xW    = WspdW;
    yW    = dataW;
    xS    = WspdS;
    yS    = dataS;
    
    if(sum(yLimits{iVar}==[-7.5 7.5])==2)
        nanIdxS  = yS>7.5 | yS<-7.5;
        nanIdxW  = yW>7.5 | yW<-7.5;
        yS = yS(~nanIdxS);
        yW = yW(~nanIdxW);
        xS = xS(~nanIdxS);
        xW = xW(~nanIdxW);
    end
    
    if(sum(xLimits{iVar}==[-1 7])==2)
        nanIdxS  = xS>7;
        nanIdxW  = xW>7;
        yS = yS(~nanIdxS);
        yW = yW(~nanIdxW);
        xS = xS(~nanIdxS);
        xW = xW(~nanIdxW);
    end
    
    if(sum(yLimits{iVar}==[-1 7])==2)
        nanIdxS  = yS>7;
        nanIdxW  = yW>7;
        yS = yS(~nanIdxS);
        yW = yW(~nanIdxW);
        xS = xS(~nanIdxS);
        xW = xW(~nanIdxW);
    end
    
    dataS = dataS(~isnan(dataS));
    dataW = dataW(~isnan(dataW));
    WspdS = WspdS(~isnan(WspdS));
    WspdW = WspdW(~isnan(WspdW));

    pdfS      = pdfFit(dataS);
    pdfW      = pdfFit(dataW);
    pdfWspdS  = pdfFit(WspdS);
    pdfWspdW  = pdfFit(WspdW);


    X1 = dataS;
    X2 = dataW;
    X3 = WspdS;
    X4 = WspdW;

    Y1 = pdfS;
    Y2 = pdfW;
    Y3 = pdfWspdS;
    Y4 = pdfWspdW;

    
    plot2dTrop
    
end

        


