sites={'seitseminen','laegeren','dundo'};

font_size=6;

figure;

for isite = 1:length(sites)
    filename = char(strcat('D:\Data\LPJmL-FIT\Validation\trait_ind_',sites(isite),'_const.txt'));
    data = readtable(filename);

    data.SLA=data.SLA*455;
    data.Wooddens=data.Wooddens*(10^-6)/0.455;

    %normalization
    SLA=[0.005,0.07]*455;
    LL=[0.2,20];
    WD=[0.7e5,6.5e5]*(10^-6)/0.455;

    rows = data.Year==2013;
    data=data(rows,:);

    rows = data.Longevity<20;
    data=data(rows,:);

    rows = WD(1)<data.Wooddens;
    data=data(rows,:);
    rows = data.Wooddens<WD(2);
    data=data(rows,:);

    if isite==3
        rows = 5>data.Cell;
    else
        rows = 20>data.Cell;
    end
    data=data(rows,:);
    
    columns={'SLA','Wooddens','Longevity'};
    data_convhull=table2array(data(:,columns));
    
    colors=data.PFT;

    %[255, 255, 0]; %yellow, BL-SG
    col_BLSG=[255, 220, 0];
    col_BLEG=[255, 153, 51];
    col_BNLEG=[0, 153, 255];
    col_TNLEG=[255, 80, 80];
    
   %% distribution plot

    X = data_convhull;

    p_dist = subplot(3,4,(isite*4-3));
    hold on
    
    plot3(p_dist, X(colors==0,1), X(colors==0,3), X(colors==0,2), '.' ,'color',col_BLSG/255);
    plot3(p_dist, X(colors==1,1), X(colors==1,3), X(colors==1,2), '.' ,'color',col_BLEG/255);
    plot3(p_dist, X(colors==2,1), X(colors==2,3), X(colors==2,2), '.' ,'color',col_BNLEG/255);
    plot3(p_dist, X(colors==3,1), X(colors==3,3), X(colors==3,2), '.' ,'color',col_TNLEG/255);

    axis(p_dist,[ SLA LL WD ]);
    grid on
    view( -57,26 );

    title('');

    xlabel('SLA [mm²/mg]');
    ylabel('Leaf Longevity [a]');
    zlabel('Wood Density [g/cm³]');
    set( gca, 'FontName', 'Arial', 'FontSize', font_size ,'Xdir','reverse','Ydir','reverse');


    maxP = 0;
    for side = 1:3
        % Calculate the histogram.
        p = [0 hist(X(:,side), 15) 0];
        p = p / sum(p);
        maxP = max([maxP p]);
        dx = (max(X(:,side)) - min(X(:,side))) / numel(p) / 2.3;
        p2 = [zeros(1,numel(p)); p; p; zeros(1,numel(p))]; p2 = p2(:);
        x = linspace(min(X(:,side)), max(X(:,side)), numel(p));
        x2 = [x-dx; x-dx; x+dx; x+dx]; x2 = max(min(x2(:), max(X(:,side))), min(X(:,side)));

        % Calculate the curve.
        nPtsCurve = numel(p) * 10;
        xx = linspace(min(X(:,side)), max(X(:,side)), nPtsCurve);

        % Plot the curve and the histogram.
        if side == 1 %SLA
            plot3(p_dist,xx, ones(1, nPtsCurve)* 20 , spline(x,p,xx)*max(X(:,2))/maxP/2, 'r-', 'LineWidth', 2);
            %plot3(x2, ones(numel(p2), 1)* 20 , p2*max(X(:,2))/maxP/2, 'k-', 'LineWidth', 1);
        else
            if side ==2 %WD
                plot3(p_dist,ones(1, nPtsCurve)* max(SLA), spline(x,p,xx)*max(X(:,3))/maxP/2,xx, 'r-', 'LineWidth', 2);
                %plot3(ones(numel(p2), 1)* max(SLA), p2*max(X(:,3))/maxP/2, 'k-',x2, 'LineWidth', 1);
            else %LL
                   plot3(p_dist,ones(1, nPtsCurve)* max(SLA),xx,spline(x,p,xx)*max(X(:,2))/maxP/2, 'r-', 'LineWidth', 2);
                   %plot3(ones(numel(p2), 1)* max(SLA),x2, p2*max(X(:,2))/maxP/2, 'k-', 'LineWidth', 1);
            end
        end

    end

    hold off


    %% calculate Richness
    columns={'SLA','Wooddens','Longevity'};
    data_convhull=table2array(data(:,columns));
    [K, v] = convhull(data_convhull);


    p_rich = subplot(3,4,(isite*4-2));
    hold on
    
    plot3(p_rich, data_convhull(colors==0,1), data_convhull(colors==0,3), data_convhull(colors==0,2), '.' ,'color',col_BLSG/255);
    plot3(p_rich, data_convhull(colors==1,1), data_convhull(colors==1,3), data_convhull(colors==1,2), '.' ,'color',col_BLEG/255);
    plot3(p_rich, data_convhull(colors==2,1), data_convhull(colors==2,3), data_convhull(colors==2,2), '.' ,'color',col_BNLEG/255);
    plot3(p_rich, data_convhull(colors==3,1), data_convhull(colors==3,3), data_convhull(colors==3,2), '.' ,'color',col_TNLEG/255);
    tri = trisurf( K, data_convhull(:,1), data_convhull(:,3), data_convhull(:,2), 'FaceColor', 'r');

    colormap( gray );
    shading interp
    alpha( tri, 0.15 );
    camlight('headlight');
    material dull

    %axis image
    axis([ SLA LL WD ]);
    grid on
    view( -57,26 );

    title('');

    xlabel('SLA [mm²/mg]');
    ylabel('Leaf Longevity [a]');
    zlabel('Wood Density [g/cm³]');
    set( gca, 'FontName', 'Arial', 'FontSize', font_size ,'Xdir','reverse','Ydir','reverse');


    hold off
   

    %% plot functional divergence


    X = data_convhull;

    S = size( X, 1 );

    % FDiv Functional Divergence by Villéger, Mason, Mouillot 2008
    g = [ mean( X(:,1) ), mean( X(:,3) ), mean( X(:,2) ) ];

    dG = pdist( cat( 1, g, X ), 'euclidean' );
    dG = dG( 1:S );
    dG_mean = mean( dG );

    deltaD = sum( 1./S .* ( dG - dG_mean ) );
    deltaAbsD = sum( 1./S .* abs( dG - dG_mean ) );

    radius = ( deltaD + dG_mean ) / ( deltaAbsD + dG_mean ) / 2; % fdiv divided by 2 results in the radius of the sphere to be plotted

    subplot(3,4,(isite*4-1));
    hold on
    plot3( X(colors==0,1), X(colors==0,3), X(colors==0,2), '.' ,'color',col_BLSG/255);
    plot3( X(colors==1,1), X(colors==1,3), X(colors==1,2), '.' ,'color',col_BLEG/255);
    plot3( X(colors==2,1), X(colors==2,3), X(colors==2,2), '.' ,'color',col_BNLEG/255);
    plot3( X(colors==3,1), X(colors==3,3), X(colors==3,2), '.' ,'color',col_TNLEG/255);

    h = plot3( g(:,1), g(:,2), g(:,3), 'x', 'MarkerSize',20 );
    set( h, 'LineWidth', 5 );


    %ylim=get(gca,'YLim');
    %xlim=get(gca,'XLim');
    %zlim=get(gca,'ZLim');
    %rad_coef=max([diff(xlim),diff(ylim),diff(zlim)]);

    %radius=radius*rad_coef;

    [x,y,z] = sphere(100);
    sp = surf( x.*radius*diff(SLA) + g(:,1), y.*radius*diff(LL) + g(:,2), z.*radius*max(WD) + g(:,3) );

    colormap( gray );
    shading interp
    alpha( sp, 0.15 );
    camlight('headlight');
    material dull

    %axis image
    axis([ SLA LL WD ]);
    grid on
    view( -57,26 );

    title('');

    xlabel('SLA [mm²/mg]');
    ylabel('Leaf Longevity [a]');
    zlabel('Wood Density [g/cm³]');
    set( gca, 'FontName', 'Arial', 'FontSize', font_size ,'Xdir','reverse','Ydir','reverse');

    hold off

 
   
 %% plot functional evenness

    % FEeve Functional Evenness by Villéger, Mason, Mouillot 2008
    X=data_convhull;

    EW=cell2mat(struct2cell(load(char(strcat('D:\Data\LPJmL-FIT\Validation\',sites(isite),'_2013_20cells.mat')))));
   
    subplot(3,4,(isite*4));
    hold on
    plot3( data_convhull(colors==0,1), data_convhull(colors==0,3), data_convhull(colors==0,2), '.' ,'color',col_BLSG/255);
    plot3( data_convhull(colors==1,1), data_convhull(colors==1,3), data_convhull(colors==1,2), '.' ,'color',col_BLEG/255);
    plot3( data_convhull(colors==2,1), data_convhull(colors==2,3), data_convhull(colors==2,2), '.' ,'color',col_BNLEG/255);
    plot3( data_convhull(colors==3,1), data_convhull(colors==3,3), data_convhull(colors==3,2), '.' ,'color',col_TNLEG/255);

    for m = 1:size( EW, 1 )
        for n = 1:size( EW, 2 )
            if EW( n, m ) ~= 0
                plot3( cat(1,X(m,1),X(n,1)), cat(1,X(m,3),X(n,3)), cat(1,X(m,2),X(n,2)), 'k' );
            end
        end
    end

    %axis image
    axis([ SLA LL WD ]);
    grid on
    view( -57,26 );

    title('');

    xlabel('SLA [mm²/mg]');
    ylabel('Leaf Longevity [a]');
    zlabel('Wood Density [g/cm³]');
    set( gca, 'FontName', 'Arial', 'FontSize', font_size ,'Xdir','reverse','Ydir','reverse');

    hold off


    



end


    set(gcf, 'PaperSize', [30 14],...
        'PaperPosition',[0 0 30 14],...
        'PaperUnits','centimeters'); %

    %print (gcf, '-dpng' ,'-r600', [pathFig 'fric_' app{i} '.png']);
    print (gcf, '-dpng','-r1200', strcat('D:\Data\LPJmL-FIT\Validation\all_scatter_matlab.png'));
