'this script uses the phase plot of the CRab Analyzer
' it divides each cycle of the rhythm into the same NUMBER of Bins, that is, the BINwidth may change from cycle to cycle
' it then calculates a histogram of
' 1. the total number of spikes in one bin (counts all spike of all cycles)
' 2. the average number of spikes per bin (divides the total number by the number of cycles)
' 3. the normalized number of spikes per bin (normalized to the average number of spikes / cycle). The average numberof spikes/cycle = 1, all other Bins contain a fraction of this
' 4. the average firing frequencies per bin 

'HERE, all spikes are counted for each Bin, summed up and averaged later

var a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,period,bin,spikefreq,allspikefrequ,allspikefrequs[1000],together;
var spikenum,allspike,allspikeno[1001],allspikebin[1001];
var kanaeleanzahl,spikekanaele[10],kanaeletext$,mainwin%,log%,yes$,logfilename$;

var bins[1001];





mainwin%:=View();
log%:=LogHandle(); 'get the logfile handle

'Clear log file
View(log%); 'go to logfile
editselectall();
editclear();
View(mainwin%); 'and back to main window
logfilename$:=filename$(3);


message("Use Crab Analyzer to do classical phase analysis first!");
'a:=input("number of channel with spike events?",401);
'b:=a+1;
b:=input("number of reference channel with bursts for phase analysis?",b);

c:=input("give % BINsize for histogram (min. 0.1% !)",2);
c:=c/100; 'calculate in real nunmber no per cent 


d:=0;
f:=1; 'just to make it bigger than 0 

'aks for number of channel that contains the bursts of the reference neuron

kanaeleanzahl:=input("How many neurons do you want to analyze?",1);

for u:=1 to kanaeleanzahl do
	kanaeletext$:="Give channel number for spike events of neuron " + STr$(u);
	spikekanaele[u]:=input(kanaeletext$,spikekanaele[u-1]);
next;

for u:=1 to kanaeleanzahl do 'Loop for analyzing each neuron



 '	FIND Bursts on reference channel (we repeat this for each neuron, although not really necessary)

	While f>0 do
		'search for burst begin and end on reference channel
		d:=nexttime(b,d); 'set d to begin of burst
		e:=nexttime(b,d); 'set e to end of burst
		f:=nexttime(b,e); 'set f to next begin of burst

		
		if f>0 then 'if last burst was found do not calculate this

	
			period:=f-d; 'calculate period
	
			'divide period in % bins
			bin:=period*c;
	
			' ok , we have 1/c numbers of bins, each bin has the size of bin
			' now, lets see if we can find spikes in these bins and remember them later.
		
			g:=bin; 'g holds the time that will be incremented during a cycle (the bin width)
		
			REPEAT 'step through all BINs and count spikes and calculate firing frequencies
				h:=h+1; 'h holds the number of the bin
				'Cursor(3,d+g-bin);
				'Cursor(4,d+g);
				'Spikecount *******
				spikenum:=count(spikekanaele[u],d+g-bin,d+g); 'count the spikes in this bin
				allspike:=spikenum+allspike; 'count all spikes in one cycle
				allspikebin[h]:=allspikebin[h]+spikenum; 'fills allspikebin array with spikecount and remembers it
	
				'spikefrequencies ********
				spikefreq:=spikenum; 'copy the number of spikes to spikefrequ
				spikefreq:=spikefreq/bin; 'calculate firing frequencies
				bins[h]:=bins[h]+spikefreq; 'fill bins array with spikefrequency and remember it

				g:=g+bin; 'increment g by bin size -> jump to next bin
			UNTIl g>period-bin; 'until g is bigger than the period minus one bin size. AFter that, we need to start a new cycle.
		
			h:=0; 'reset bin number
			g:=bin; 'reset incremented g
			d:=e; 'set d to end of burst, so that nexttime in next loop will find correct burst begin and end
			n:=n+1; 'n counts the number of bursts / periods that are calculated.
			allspikeno[n]:=allspike; 'remember number of spikes for this cycle
			allspike:=0; 'reset counter for spikes for the next cycle
		endif;

	Wend;

	getout();

Next; 



'****************** go to logfile and save the damn thing


View(log%); ' show it
editselectall();
editcopy();

yes$:=Input$("results were copied to clipboard, save (Excel) file now? (y/n)","n",1);

if yes$="y" then
FileSaveAs(logfilename$,1); ' save it
endif;

View(mainwin%); 'and back to main window

halt;
'******************************




'get data out of spike2 
proc getout();
	Printlog("Analyzing neuron	",STR$(u));
	Printlog("Neuron name:	", Chantitle$(spikekanaele[u]));
	printlog("the_n-number_is	",n,"	and_the_rel_bin_size_in_%_is	",c*100);
	Printlog("number of spikes / cycle:");
	Printlog("\n");
	
	for x:=1 to n do
	
		printlog(x,"	",allspikeno[x]);
		together:=together+allspikeno[x];
	next;
	
	Printlog("total_number	",together);
	Printlog("average	",together/n);
	Printlog("\n");	
	
	Printlog("bin#	","total_number_of_spikes_for_this_BIN	","average_number_of_spikes_for_this_BIN	","normalized_to_average_spikecount_per_cycle	","average_spikefrequency/BIN");
	
	for x:=1 to 1/c+1 do '1/c = number of bins
 	
		printlog(x,"	",allspikebin[x],"	",allspikebin[x]/n,"	",(allspikebin[x]/n)/(together/n),"	",bins[x]/n); 
	
	next;
	
	Printlog("\n");

	Printlog("XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX");
	Printlog("\n");

	'RESET VARIABLES
		d:=0;
		
		for x:=0 to 1/c+1 do
			bins[x]:=0;
		next;

		h:=0;
		spikefreq:=0;
		allspike:=0;
		spikenum:=0;
		g:=0;
		
		for x:=0 to 1/c+1 do
			allspikebin[x]:=0;
		next;
		
		n:=0;
		f:=1; 'just to make it bigger than 0;
		together:=0;


return;

end;