function [grand]=frate_simu;
% It plots all trajectories segments coincident with the spikes
% of the specified electrode and neuron class superimposed.
% The sample of the spike triggered average.
% Example:  frate_simu;

elext={'B' 'C' 'D' 'E' 'F' 'G' 'H' 'I'};
unitext={'A' 'B' 'C' 'D' 'E'};
electrodes=8;
nperel=4;

minx=-2000;
maxx=2000;
miny=-1500;
maxy=1500;
nbins=40;

fpath1='../analog/xy_today.';
fpath2='../final/ses.';

sessions=4;
from=[9 33 57 81];
to=[24 48 72 96];

xax=maxx-minx;
yax=maxy-miny;
unitx=xax/nbins;
unity=yax/nbins;

figure(3);clf;
figure(4);clf;
figure(5);clf;
figure(6);clf;
figure(7);clf;
figure(9);clf;
figure(10);clf;

grandsum=zeros(nbins,nbins);
grandfr=zeros(nbins,nbins);
shist=zeros(1,101);
msess=0;
% calculate the time density first ----------------------
t=zeros(nbins,nbins);
st=zeros(nbins,nbins);
for j=1:sessions
      for i=from(j):to(j)
	fnum=num2str(i);
	if (i<10) ext=strcat('00',fnum);end
	if (i<100 & i>9) ext=strcat('0',fnum);end
	if (i>99) ext=fnum;end
	fnamexy=strcat(fpath1,ext)
	[t]=map_timespent(fnamexy,0);
	st=st+t;
	clear t;
      end      
end

figure(4);
pcolor(st');
colormap(jet);
shading flat;
colorbar;
xlabel(strcat('x [',num2str(minx),',',num2str(maxx),']'));
ylabel(strcat('y [',num2str(miny),',',num2str(maxy),']'));
title(['Time density map  (File# ',num2str(from),'..',num2str(to),')']);
saveas(gcf, 'timedens.fig', 'fig');



% then the individual spike density maps for each cells---
for k=1:electrodes
  for l=1:nperel
%    tgv=zeros(1,10000);
    m=zeros(nbins,nbins);
    n=zeros(nbins,nbins);
    sm=zeros(nbins,nbins);
    fr=zeros(nbins,nbins);
    prevnspk=1;

    el=char(elext(k));
    nclas=char(unitext(l));
    nsess=0;
    for j=1:sessions
      for i=from(j):to(j)
	fnum=num2str(i);
	if (i<10) ext=strcat('00',fnum);end
	if (i<100 & i>9) ext=strcat('0',fnum);end
	if (i>99) ext=fnum;end
	fnamexy=strcat(fpath1,ext)
	fnamespk=strcat(fpath2,ext,'.ten.00',el,'.00',nclas)
	[m,n,tgv,nspk]=map_xyspk(fnamexy,fnamespk,'b-',15);
	if (nspk>0) nsess=nsess+1;
	sm=sm+m;
	[N,X]=hist(tgv,[0:0.01:1]);
	shist=shist+N;
%	endnspk=nspk+prevnspk-1;
%	tgvarray([prevnspk:endnspk])=tgv([1:nspk]);
%	prevnspk=endnspk+1;
        clear n
	end;
      end
    end
    a(k,l,[1:nbins],[1:nbins])=sm;
    grandsum=grandsum+sm;
    msess=msess+nsess;
    
    figure(1);
    H=shist/nsess;
    bar([1:101],H);
    title('tang vel');
    
    figure(3);
    subplot(electrodes,nperel,(k-1)*nperel+l);
    pcolor(sm');
    colormap(jet);
    shading flat;
    colorbar;
    if (k==1 & l==1) 
      xlabel(strcat('x [',num2str(minx),',',num2str(maxx),']'));
      ylabel(strcat('y [',num2str(miny),',',num2str(maxy),']'));
    else
      axis off
    end
    hold on;
    
% calculate the local firing rate on each pixel
% by deviding the spike density over pixel by the 
% time density (the time the trajectory spent inside 
% the pixel
    % firing rate = total spike count / time spent * 1000 [spike/sec]
    fr=(sm./st)*1000; 
    grandfr=grandfr+fr;
    b(k,l,[1:nbins],[1:nbins])=fr;

    
    figure(5);
    subplot(electrodes,nperel,(k-1)*nperel+l);
    pcolor(fr');
    colormap(jet);
    shading flat;
    me=mean(mean(fr(find(fr>0))));
    if (me>0) caxis([0 2*me]);end;
%    me=mean(mean(fr));
%    if (me~=0) caxis([0 7*me]); end;
%    caxis([0 0.01]);
    colorbar;
    if (k==1 & l==1) 
      xlabel(strcat('x [',num2str(minx),',',num2str(maxx),']'));
      ylabel(strcat('y [',num2str(miny),',',num2str(maxy),']'));
    else
      axis off;
    end
    hold on;

%    figure(6);
%    subplot(electrodes,nperel,(k-1)*nperel+l);
%    hist(tgvarray,40);
%    xlabel('tang. vel.');
%    xlabel('spikes');
%    hold on;
%    axis([0 15000 0 1000]);
  
  end
end

pp1=zeros(nbins,nbins);
sdenom=zeros(nbins,nbins);
pp2=zeros(nbins,nbins);
for k=1:electrodes
  for l=1:nperel
    pa=zeros(nbins,nbins);
    pb=zeros(nbins,nbins);
    sqa=zeros(nbins,nbins);
    sqb=zeros(nbins,nbins);
    for kk=1:electrodes
    for ll=1:nperel
      if ~(k==kk & l==ll)
	    pa=pa+(a(k,l).*a(kk,ll));
	    pb=pb+(b(k,l).*b(kk,ll));
	    sqa=sqa+(a(k,l).*a(k,l));
	    sqb=sqb+(b(kk,ll).*b(kk,ll));
      end
    end
    end
    pp1=pp1+pa;
    pp2=pp2+pb;
    sdenom=sdenom+sqrt(sqa)*sqrt(sqb);
end
end

     pred=pp1./sdenom;
     pred2=pp2./sdenom;

figure(3);
textglobal(0.1,0.95,['Spike density map  (File# ',num2str(from),'..',num2str(to),'   El: ',el,'   Ncl: ', nclas,')']);
saveas(gcf, 'spikedens.fig', 'fig');
figure(5);
textglobal(0.1,0.95,['Firing rate map (File# ',num2str(from),'..',num2str(to),'   El: ',el,'   Ncl: ', nclas,')']);
legend('[spike/sec]');
saveas(gcf, 'spikefreq.fig', 'fig');
%figure(6);
%textglobal(0.1,0.95,['Tang. vel. (File# ',num2str(from),'..',num2str(to),'   El: ',el,'   Ncl: ', nclas,')']);


figure(7);
clf;
%pcolor(pp1');
pcolor(pred');
colormap(jet);
shading flat;
colorbar;
      xlabel(strcat('x [',num2str(minx),',',num2str(maxx),']'));
      ylabel(strcat('y [',num2str(miny),',',num2str(maxy),']'));
      textglobal(0.1,0.95,['Prod. of spike densities superimposed (File# ',num2str(from),'..',num2str(to),')']);
saveas(gcf, 'allneursuperpos.fig', 'fig');

figure(9)
Mfr=grandfr/nsess;
pcolor(Mfr');
colormap(jet);
shading flat;
colorbar;
      xlabel(strcat('x [',num2str(minx),',',num2str(maxx),']'));
      ylabel(strcat('y [',num2str(miny),',',num2str(maxy),']'));
      title(['Superpos all firing rates / neuron (File# ',num2str(from),'..',num2str(to),')']);
saveas(gcf, 'fratesuperpos.fig', 'fig');

figure(10)
pcolor(grandsum');
colormap(jet);
shading flat;
colorbar;
      xlabel(strcat('x [',num2str(minx),',',num2str(maxx),']'));
      ylabel(strcat('y [',num2str(miny),',',num2str(maxy),']'));
      title(['Superpos all firing densities / neuron /pixels (File# ',num2str(from),'..',num2str(to),')']);
saveas(gcf, 'spkdenssuperpos.fig', 'fig');
save workspace.mat;