function [target_window,circle_r] = PlotTuningCurve(ev,ts,cnd,varargin)
%PlotTuningCurve(ev, ts , cnd,varargin)
%   Functions for ploting tuning curve, the ev is the event time, ts is the
%   spike's timestamp
%
iod = @utils.inputordefault;
pre = iod('pre',0,varargin);
post = iod('post',1,varargin);
binsz = iod('binsz',0.01,varargin);
krn = iod('krn',0.1,varargin);
tar_pre = iod('tar_pre',0.25,varargin);
tar_post = iod('tar_post',0.25,varargin);
target_time = iod('target_time',[],varargin);% the time point as the center of the tuning curve time window
clr = iod('color',[0 0 0],varargin);
titlestr = iod('titlestr',[],varargin);
sub_cnd = iod('sub_cnd',[],varargin);
r_lim = iod('rlim',[],varargin);
isplot = iod('isplot',1,varargin);
savedir = iod('savedir',[],varargin);
cellid = iod('cellid',0,varargin);
ax = iod('ax',[],varargin);
circle_r = iod('circle_r',[],varargin);
[headlength,varargin]=iod('headlength',3,varargin);
[headstyle,varargin]=iod('headstyle','cback3',varargin);
[linewidth,varargin] = iod('linewidth',2,varargin);
if iscell(cnd)
    %cnd_nan = cellfun(@(x)isnan(x), cnd);
    cnd_nan = cellfun(@(x)strcmp(x,'null'), cnd);
    %cnd(cnd_nan) = 'NaN';
    cnd = categorical(cnd);
    n_cnd = categories(cnd);
else
    cnd = categorical(cnd);
    n_cnd = categories(cnd);
end

if ~isempty(sub_cnd) && iscell(cnd)
    %cnd_nan = cellfun(@(x)isnan(x), cnd);
    sub_cnd_nan = cellfun(@(x)strcmp(x,'null'), sub_cnd);
    %cnd(cnd_nan) = 'NaN';
    sub_cnd = categorical(sub_cnd);
    n_sub_cnd = categories(sub_cnd);
else
    sub_cnd = categorical(sub_cnd);
    n_sub_cnd = categories(sub_cnd);
end

if isscalar(krn)
    dx=ceil(5*krn);
    kx=-dx:binsz:dx;
    krn=normpdf(kx,0, krn);
	if isempty(find(kx==0, 1))
		error('Your binsz needs to divide 1 second into interger # of bins');
	end
    krn(kx<0)=0;
    krn=krn/sum(krn);
end

[Y,x]=stats.spike_filter(ev,ts,krn,'pre',pre,'post',post,'kernel_bin_size',binsz);
% [refY,x] = stats.spike_filter(ref_ev,ts,krn,'pre',pre,'post',post,'kernel_bin_size',binsz);
ymn = nan(numel(n_cnd),numel(x));
for ci=1:numel(n_cnd)
    ref=cnd==n_cnd(ci); 
    y=Y(ref,:); 
    ymn(ci,:) = nanmean(y,1);
end

if isempty(target_time) %is the user is not prividing target_time for doing tuning curve, we find a time window with maxinum variance
    % take nanvar?
    [~,I] = sort(nanvar(ymn));
    target_time = x(I(end));
end
target_window = [max([x(1),target_time-tar_pre]),min([x(end),target_time+tar_post])];

% get firing rate in the target_window
tar_fr = nan(numel(n_cnd),1);
tar_se = nan(numel(n_cnd),1);

for ci = 1:numel(n_cnd)
    ref = cnd==n_cnd(ci);
    tar_fr(ci) = nanmean(stats.qcount(ts,ev(ref)+target_window(1),ev(ref)+target_window(2))./(target_window(2)-target_window(1)));
    tar_se(ci) = nanstd(stats.qcount(ts,ev(ref)+target_window(1),ev(ref)+target_window(2))./(target_window(2)-target_window(1)))/(nansum(ref)^0.5);
end
if isempty(sub_cnd)
    angles = cellfun(@(x)utils.ExtractAngleInfo('MidC',cellstr(x),'Rat'), n_cnd);
else
    angles = cellfun(@(x,y)utils.ExtractAngleInfo(cellstr(x),cellstr(y),'Rat'), n_cnd,{n_sub_cnd{1};n_sub_cnd{2};n_sub_cnd{1};n_sub_cnd{2}});
end
[~ ,sort_index]=sort(angles);% angles from MidR, counterclockwise
sort_index = [sort_index;sort_index(1)];
theta = angles(sort_index);                                                           
[x1, y1] = pol2cart(theta, tar_fr(sort_index));         % Convert To Cartesian
vecsum = nanmean([x1,y1],1);
[x2, y2] = pol2cart(theta, tar_fr(sort_index)-tar_se(sort_index));
[x3,y3] = pol2cart(theta,tar_fr(sort_index)+tar_se(sort_index));
[~,temp] = max(tar_fr(sort_index));
max_theta = theta(temp);

if isplot
    if isempty(ax)
        ax = draw.jaxes();
    else
        axes(ax);
    end
hold on;

if isempty(circle_r)
circle_r = ceil(nanmax(abs(tar_fr))*1.1);
end
lim = circle_r;% axis limits

set(gca,'XTickLabel',[],'YTickLabel',[],'XTick',[],'YTick',[],'XColor','none','YColor','none','XLim',[-lim lim],'YLim',[-lim lim]);
% start plotting

draw_circle(ax,0,0,circle_r);
draw_lines(ax,0,0,angles,circle_r,headstyle,linewidth,headlength);
p_out = patch(x3,y3,clr + [0.8 0.8 0.8], 'EdgeColor','none') ;    % Outer std  % Fill Area Between Radius Limits
p_in = patch(x2,y2,[1 1 1],'EdgeColor','none');% Inner std
p_mean = plot(x1, y1, '-o','MarkerSize',2,'MarkerFaceColor',clr,'MarkerEdgeColor',clr,'Color',clr); %mean line


text(0,-lim,[titlestr,' ',num2str(circle_r),' Hz '],'HorizontalAlignment','center','VerticalAlignment','bottom');

else
    h = nan;
end

if ~isempty(savedir)
   print('-dpdf','-r300',figh, fullfile(savedir,sprintf('tune.cell.%d.pdf',cellid)));
   saveas(figh, fullfile(savedir,sprintf('tune.cell.%d.png',cellid)));

end
end
function draw_circle(ax,x,y,r)
%x and y are the coordinates of the center of the circle
%r is the radius of the circle
%0.01 is the angle step, bigger values will draw the circle faster but
%you might notice imperfections (not very smooth)
ang=0:0.01:2*pi; 
xp=r*cos(ang);
yp=r*sin(ang);
plot(ax,x+xp,y+yp,'Color','k');
plot(ax,x,y,'Color','k'); % Center of the Circle
end
function draw_lines(ax,x,y,angles,r,headstyle,linewidth,headlength)
ax_left = ax.Position(1);
ax_bot = ax.Position(2);
ax_wid = ax.Position(3);
ax_hei = ax.Position(4);
xp = r.* cos(angles);
yp = r.* sin(angles);
for ii = 1:6
     clr = draw.get_dirpos_clrs(utils.ang2name(angles(ii)));
    x_start = (x+xp(ii)*0.9-ax.XLim(1))/(ax.XLim(2)-ax.XLim(1))*ax_wid+ax_left;
    x_end = (x+xp(ii)*1.2-ax.XLim(1))/(ax.XLim(2)-ax.XLim(1))*ax_wid+ax_left;
    y_start = (y+yp(ii)*0.9-ax.YLim(1))/(ax.YLim(2)-ax.YLim(1))*ax_hei+ax_bot;
    y_end = (y+yp(ii)*1.2-ax.YLim(1))/(ax.YLim(2)-ax.YLim(1))*ax_hei+ax_bot;
    annotation('arrow',[x_start,x_end],[y_start,y_end],'Color',clr{1},'HeadStyle',headstyle,'LineWidth',linewidth,'HeadLength',headlength,'HeadWidth',5);
end
end