function NodPlotResults17(pathexcel,figpth)


%% Locate data etc
%if ~exist('mpd')
%    mpd=uigetdir(pwd,['Select the main project folder']);
%end
% if ~exist([mpd '\Excel files']);
%     disp('Excel files folder does not exist')
% end; %cd([mpd '\Excel files'])
% 
% if ~exist([mpd '\Figures2']);
%     mkdir([mpd '\Figures2']);   
% end; %cd([mpd '\Figures'])


%% SET FIGURE PARAMETERS
w=3; h=3;
fmts={'fig','eps','png'};
%figpth=[mpd '\Figures2'];
xlab='Lateral root axis (cm)';
ylab='Primary root axis (cm)';
fsz=10;
%figpth=pathfigs;

%% NODULE AND ROOT DENSITY MAPS
N='NoduleCoordinates';
clearvars -except mpd N w h fsz fmts figpth xlab ylab  ax pathexcel
import_as_column_vectors2([pathexcel '\' N '.xls'])
%ax=[0  floor(max(max(DistAlongLateral))) 0 floor(max(max(DistAlongPri)))];
ax=[0 20 0 23];
nodRtID=RootID(RootID>0);
xpos=DistAlongLateral(RootID>0);
ypos=DistAlongPri(RootID>0);
nodPltID=PlantID(RootID>0);
w=3; h=3;

N='RootData';
import_as_column_vectors2([pathexcel '\' N '.xls'])
RootDistAlongPri=RootDistAlongPri+.00001;

inc=1; disp('Densities are computed per grid square per plant, each grid square corresponding to 1 cm2');
%NOTE: values are computed for each sqaure grid of 1 cm2, so densities are 'per cm2', the program is not currently designed to handle 
%square grids of other sizes than 1cm2 - if user would like to have another
%grid size, e.g. 3x3 cm, the formulas would have to be adjusted accordingly

xg=0:inc:ax(2)+inc;
yg=0:inc:ax(4)+inc;
NodIn=[]; RtIn=[]; NodRtRatio=[];
Plts=unique(PlantID);
for plti=1:length(Plts)
    for xi=1:length(xg)-1
        for yi=1:length(yg)-1
            NodIn(yi,xi,plti)=sum(inpolygon(xpos(nodPltID==Plts(plti)),ypos(nodPltID==Plts(plti)),[xg(xi) xg(xi+1) xg(xi+1) xg(xi)],[yg(yi) yg(yi) yg(yi+1) yg(yi+1)]));
            RtIn(yi,xi,plti)=sum(RootLength>xg(xi) & RootDistAlongPri>yg(yi) & RootDistAlongPri<=yg(yi+1) & PlantID==Plts(plti));
            %AG adding info on root length in square
            indx=find(RootLength>xg(xi) & RootDistAlongPri>yg(yi) & RootDistAlongPri<=yg(yi+1) & PlantID==Plts(plti));
            RtL(yi,xi,plti)=0;
            for cp=1:length(indx)
                if RootLength(indx(cp))>=xg(xi+1)
                    RtL(yi,xi,plti)=RtL(yi,xi,plti)+1;
                else
                    RtL(yi,xi,plti)=RtL(yi,xi,plti)+RootLength(indx(cp))-xg(xi);
                end
            end
            NodRtnbRatio(yi,xi,plti)=NodIn(yi,xi,plti)./RtIn(yi,xi,plti);
            NodpermmRoot(yi,xi,plti)=NodIn(yi,xi,plti)./RtL(yi,xi,plti);   
        end
    end
end
%[xx yy]=meshgrid(0:inc:ax(2),0:inc:ax(4));

% uchoosex=input(['Choose length (cm) of lateral root (x) axis for nodulation zone maps (min ' num2str(ax(2)+inc) '): ']);
% uchoosey=input(['Choose length (cm) of primary root (y) axis for nodulation zone maps (min ' num2str(ax(4)+inc) '): ']);
% ax1=[0 uchoosex 0 uchoosey]; 



ax1=[0 ax(2)+3 0 ax(4)+3];
T='Map of mean nodule density ie mean Number of nodules per cm2(=per grid square)';
M=mean(NodIn,3);
cmax=max(max(M));

draw_gridmap2(M,cmax,size(M,1),size(M,2),0);
SaveFigs2(w,h,fsz,xlab,ylab,ax1,fmts,figpth,T)

%xlswrite([excelpth '\' T],xx,'Sheet1','A1');
%xlswrite([excelpth '\' T],yy','Sheet1','A2');
%xlswrite([excelpth '\' T],M,'Sheet1','B2');
xlswrite([pathexcel '\' T],M,'Sheet1','A1');

T='Map of standard error of nodule density';
M=std(NodIn,1,3)/sqrt(size(NodIn,3));
draw_gridmap2(M,cmax,size(M,1),size(M,2),0);
SaveFigs2(w,h,fsz,xlab,ylab,ax1,fmts,figpth,T)

% xlswrite([excelpth '\' T],xg,'Sheet1','A1');
% xlswrite([excelpth '\' T],yg','Sheet1','A2');
% xlswrite([excelpth '\' T],M,'Sheet1','B2');
xlswrite([pathexcel '\' T],M,'Sheet1','A1');



T='Map of root density ie mean Root length in cm per cm2(=per grid square)';
M=mean(RtL,3);
cmax=max(max(M));
draw_gridmap2(M,cmax,size(M,1),size(M,2),0);
SaveFigs2(w,h,fsz,xlab,ylab,ax1,fmts,figpth,T)

% xlswrite([excelpth '\' T],xg,'Sheet1','A1');
% xlswrite([excelpth '\' T],yg','Sheet1','A2');
% xlswrite([excelpth '\' T],M,'Sheet1','B2');
xlswrite([pathexcel '\' T],M,'Sheet1','A1');

T='Map of standard error of root density';
M=std(RtL,1,3)/sqrt(size(NodIn,3));
draw_gridmap2(M,cmax,size(M,1),size(M,2),0);
SaveFigs2(w,h,fsz,xlab,ylab,ax1,fmts,figpth,T)

% xlswrite([excelpth '\' T],xg,'Sheet1','A1');
% xlswrite([excelpth '\' T],yg','Sheet1','A2');
% xlswrite([excelpth '\' T],M,'Sheet1','B2');
xlswrite([pathexcel '\' T],M,'Sheet1','A1');

T='Map of nodule-to-root ratio ie mean nb nodule per cm of Root per cm2(=per grid square)';
NodpermmRoot(isnan(NodpermmRoot))=0;
M=mean(NodpermmRoot,3);
cmax=max(max(M));
draw_gridmap2(M,cmax,size(M,1),size(M,2),0);
SaveFigs2(w,h,fsz,xlab,ylab,ax1,fmts,figpth,T)

% xlswrite([excelpth '\' T],xg,'Sheet1','A1');
% xlswrite([excelpth '\' T],yg','Sheet1','A2');
% xlswrite([excelpth '\' T],M,'Sheet1','B2');
xlswrite([pathexcel '\' T],M,'Sheet1','A1');

T='Map of standard error of nodule-to-root ratio';
M=std(NodpermmRoot,1,3)/sqrt(size(NodIn,3));
draw_gridmap2(M,cmax,size(M,1),size(M,2),0);
SaveFigs2(w,h,fsz,xlab,ylab,ax1,fmts,figpth,T)

% xlswrite([excelpth '\' T],xg,'Sheet1','A1');
% xlswrite([excelpth '\' T],yg','Sheet1','A2');
% xlswrite([excelpth '\' T],M,'Sheet1','B2');
xlswrite([pathexcel '\' T],M,'Sheet1','A1');


%% LATERAL ROOT PLOT
N='RootData';
clearvars -except mpd N w h fsz fmts figpth xlab ylab  ax pathexcel
import_as_column_vectors2([pathexcel '\' N '.xls'])
RootDistAlongPri=RootDistAlongPri+.001;
w=6; h=6; ax=[0  floor(max(max(RootLength))) 0 floor(max(max(RootDistAlongPri)))];
figure; hold on
%colmatch=1:max(RootLength)/64:max(RootLength);
colmatch2=min(RootLength):(max(RootLength)-min(RootLength))/(64-1):max(RootLength);

matrgb=jet;
for k=1:length(RootLength)
    [v colind]=min(abs(RootLength(k)-colmatch2));
    plot([0 RootLength(k)],[RootDistAlongPri(k) RootDistAlongPri(k)],'color', matrgb(colind,:))
end

caxis([min(RootLength) max(RootLength)]);
colorbar('location','SouthOutside');%,'XTick',cat(ceil(min(RootLength)):2:floor(max(RootLength))));
%colorbar('YTick',cat(2,min(RootLength), round(min(RootLength))+2:2:round(max(RootLength))-2, max(RootLength)));

disp(['Total number of lateral roots from all plants all planes: ' num2str(length(RootDistAlongPri(RootID>0)))])
%plotnod=input('Plot nodules on lateral root plot? y/n','s');
plotnod='y';

if plotnod=='y'
    N2='Nodulecoordinates';
    clearvars -except mpd N2 N w h fsz fmts figpth xlab ylab  ax xlab ylab pathexcel
    import_as_column_vectors2([pathexcel '\' N2 '.xls'])
    %ax=[0  floor(max(max(DistAlongLateral))) 0 floor(max(max(DistAlongPri)))];
    %figure
    
    %h=circle6([DistAlongLateral DistAlongPri],.15,16,'k',1);
    plot(DistAlongLateral,DistAlongPri,'.k');    % h=plot(DistAlongLateral,DistAlongPri,'ok')
    %set(h,'MarkerFaceColor',[0 0 0],'MarkerSize',[3]);
    %circle6([DistAlongLateral DistAlongPri],.15,16,'k',1/length(DistAlongPri)*200);
    %SaveFigs(w,h,fsz,xlab,ylab,ax,fmts,figpth,N2)
    disp(['Total number of nodules from all plants all planes: ' num2str(length(DistAlongPri))])
    T='Lateral root plot with nodules';
else
    T='Lateral root plot';
end

SaveFigs3(w,h,fsz,xlab,ylab,ax,fmts,figpth,T)

close all
