close all;
% load('N:\kropff\analyisis\figs\fig2a\all open data_q95.db','-mat');
% load('N:\kropff\analyisis\figs\fig2a\all car data_q95.db','-mat');
% % % load('all open data.db','-mat');
% % % load('all car data.db','-mat');
% db_all.car=db_car;
% db_all.open=db_open;
% clear db_car; clear db_open;
%     addpath('/gpfs/work/kropff/common_functions/');

dats={'ec', 'hc'};
% Haccum={};
w1=0; w2=3;
jet2=interp1(1:64,jet,(1+w1:(64-w1-w2-1)/(64-1):64-w2));

% load(['all ' dat ' data.db'],'-mat')
dang=pi/50;
dacc=20;%15;
accbins=0:dacc:200; %350
angbins=-pi:dang:pi; angbins(end)=[]; angbins=angbins+dang/2;  
cols=[ [1 0 0]; [255 175 3]/255];


for d=1:2%1:2
% dat='open';
% dat='car';
dat=dats{d};

% close all;



 
% close all;
histlim=5*dacc;



angbins=-pi:dang:pi; angbins(end)=[]; angbins=angbins+dang/2;  
Ha=zeros(numel(angbins),numel(accbins)); Oa=Ha;
Hva=Ha;

myfl=['data_fig_3.db'];

nr=numel(unique(who('-file',myfl)));
L={};
L.Fr=cell(1,-1+2*numel(accbins));
L.BFr=cell(-1+2*numel(accbins),7);
% L.Occ=zeros(1,21);
% L.Fr=L.Occ; L.BFr=L.Occ;
for i=[1 3:nr] %2 hippocampal drive not in hippocampus
    rat=['rat' num2str(i)];
    db_all.(dat)=load(myfl,rat,'-mat');
    db_all.(dat)=db_all.(dat).(rat);

    sel=true(size(db_all.(dat).accel));
    if isempty(sel)
        warning([num2str(i) 'empty']);
        continue;
    end
    

    sel=sel&db_all.(dat).t_dist>5;

    sel1=conv(double(sel), ones(25,1)/25,'same');
    sel=sel1>.99;
    selx=abs(cos(db_all.(dat).anga-db_all.(dat).angv))>.9;
    xxx=db_all.(dat).accel.*sign(cos(db_all.(dat).anga-db_all.(dat).angv));

    frw=50*5;
%     [Hs,Occs,Hv]=freq_hist(db_all.(dat).accel(sel),mod(db_all.(dat).anga(sel)-db_all.(dat).angv(sel)+pi,2*pi)-pi,db_all.(dat).freq(sel)-fr_baseline(sel),accbins,angbins);
    [Hs,Occs,Hv]=freq_hist(db_all.(dat).accel(sel),mod(db_all.(dat).anga(sel)-db_all.(dat).angv(sel)+pi,2*pi)-pi,db_all.(dat).freq(sel),accbins,angbins);

    % Linear results
    
    rxxx=round(xxx/dacc);
%     L.Occ=arrayfun(@(x) sum(rxxx(sel & selx)==x),-10:10);
    ws=50*[1 2 5 10 20 50 200];
    aend=accbins(end)/dacc;
    for w=1:numel(ws)
        fr_baseline=nanconv(db_all.(dat).freq,ones(ws(w),1)/ws(w)); %10 s window
        
        for x=-aend:aend
            xj=aend+1+x;
            if w==1
                L.Fr{xj}=[L.Fr{xj}; db_all.(dat).freq(rxxx==x & sel & selx)];
            end
            L.BFr{xj,w}=[L.BFr{xj,w}; fr_baseline(rxxx==x & sel & selx)];        
        end    
    end
    
    Ha=Ha+Hs.*Occs; Oa=Oa+Occs; Hva=Hva+Hv.*Occs;
    Haccum{d}=Ha;
    Oaccum{d}=Oa;
    Laccum{d}=L;

end

end
%%
close all
for d=1:2
angbins(101)=angbins(1)+2*pi;

f=figure;
hold off;
H=Haccum{d}; O=Oaccum{d};
H(O>0)=H(O>0)./O(O>0);
histlim=max(O(:))/10000;

H(O<histlim)=nan;
Hs=gauss_smoothing_2D(repmat(H,3,1),[1.5 3]);
h1=size(H,1);
Hs=Hs(h1+(1:h1),:);

polar3d(flipud(Hs'),angbins(1),angbins(end),accbins(1),accbins(size(Hs,2)),1,'surf'); shading flat;view(2); 
axis equal; view([-90,90]); axis tight; 
hold on;

colormap(jet2);

cc=colorbar;
set(cc,'YTick',7.5:.5:10);
name=['fig2a_' dat];
plot3([0 100],200*[1 1],[15 15],'k')
caxis([8.1 9.7])



name=['colormap_' dats{d}];
end
%%
f=figure;
accbins2=[fliplr(-accbins(2:end)) accbins];
for d=2
    L=Laccum{d};
    sel=cellfun(@numel,L.Fr);
    sel=sel>max(sel)/10000;
    mn=cellfun(@mean,L.Fr);
    eom=cellfun(@std,L.Fr)./sqrt(cellfun(@numel,L.Fr));
    errorbar(accbins2(sel),mn(sel),eom(sel),'Color',cols(d,:));
    hold on;
    plot(accbins2(sel),mn(sel),'Color',cols(d,:),'LineWidth',3);
    dat_acc=arrayfun(@(x) accbins2(x)*ones(size(L.Fr{x})),1:numel(L.Fr),'UniformOutput',0);
 
end
axis([-220 240 7.95 10]);
set(gca,'YTick',6:.5:12)

%%
f=figure;
accbins2=[fliplr(-accbins(2:end)) accbins];
for d=2
    subplot(1,2,d)
    for j=[3 ]
        
    L=Laccum{d};
    L.Fr2=arrayfun(@(x) L.Fr{x}-L.BFr{x,j},1:numel(L.Fr),'UniformOutput',0);
    sel=cellfun(@numel,L.Fr2);
    sel=sel>max(sel)/10000 & abs(accbins2)<210;
    mn=cellfun(@mean,L.Fr2);
    eom=cellfun(@std,L.Fr2)./sqrt(cellfun(@numel,L.Fr2));
    errorbar(accbins2(sel),mn(sel),eom(sel),'Color',cols(d,:)*(j-1)/6);
    hold on;
    plot(accbins2(sel),mn(sel),'Color',cols(d,:)*(j-1)/6,'LineWidth',3);
    end

    
    axis([-220 240 -.25 2]);
    golden(2.5)
set(gca,'YTick',-2:.5:12)

    
    dat_acc=arrayfun(@(x) accbins2(x)*ones(size(L.Fr2{x})),1:numel(L.Fr2),'UniformOutput',0);    
%     [cr_neg,p_neg]=corr(cell2mat(dat_acc(sel & accbins2<0)'),cell2mat(L.Fr2(sel & accbins2<0)'))
%     [cr_neg,p_neg]=corr(accbins2(sel & accbins2<0)',mn(sel & accbins2<0)');
    lneg_bas = fitlm(accbins2(sel & accbins2<0)',mn(sel & accbins2<0)','linear')
    lpos_bas = fitlm(accbins2(sel & accbins2>0)',mn(sel & accbins2>0)','linear')
    
end