'----------------------------------------------------
' Cycle Script
' 6/17/05
'----------------------------------------------------
'
'C H A N N E L S :
'------------------

'Cycle requires 5 channels of input: 1 unit and 4 analog signal channels.
'You are prompted to arrange the channel order

'10/25/06 - Jim wants the control channel zero crossings to be detected 
'from a velocity channel so horizontally-offset oscillation can be analyzed.




const
	kBOBmode								:= 0,					'0=booth mode
'	kBOBmode								:= 1,					'1=development mode
	kNaN 									:= 9999999,
	kpi 									:= 3.14159265359,
	k2pi									:= 3.14159265359 * 2,
	kRadians								:= 360.0 / (3.14159265359 * 2),
	kTrue									:= 1,
	kFalse								:= 0,
	kChannelsPlusOne					:= 5+1,

	'change these so they are x-Single-x
	kCtlTrial							:= 1,
	kEyeTrial							:= 2,

	kCtlVelTrial						:= 1,
	kEyeVelTrial						:= 2,
 
	kAux1Trial							:= 1,
	kAux2Trial							:= 1,

	kUnitTrial							:= 1,
	kISItrial							:= 2;

var
	'buttons
	gQuitButton							:= 1,
	gOpenButton							:= 2,
	gFinalButton  						:= 3,
	gBinButton							:= 4,
	gAverageButton						:= 6,
	gCurrentButton						:= 7,
	gResetButton						:= 9,
	gRejectButton						:=	10,
	gAcceptButton						:= 11,
	gDesaccadeButton					:=	12,

	'file stuff
	gDatafilepath$,
	gDatafileName$,

	'views
	gViewData%,
	gTrialView%,
	gAux1View%,
	gAux2View%,
	gVelocityView%,
	gUnitView%,
	gTrialAverageView%,
	gVelocityAverageView%,
	gAux1AverageView%,
	gAux2AverageView%,
	gUnitHistoAverageView%,
	gUnitISIAverageView%,
	gHandleAnalysis%,
	gHandleAnalysisAve%,
	gHandleAnalysisFinal%,

	gDurMS,
	gMidCycle,
	gBegCycle							:= 0,
	gEndCycle							:= 0,
	gStopAnalysis  					:= 0,
	gBegDsac[20],
	gEndDsac[20],
	gDurCycle,
	gDurCycleLast,
	gCyclesShown%,
	gCyclesAccepted%,

	gCycleIndFitHeading$					[37+1], 'resize aMean[], aSD[] similarly
	gCycleIndFitData			[100]		[37+1], 'resize aMean[], aSD[] similarly
	gCycleAvgFitData						[53+1], 
'*begin*  all must have the same size
	gMaxCycleDur%						 := 6000,	'change this and you should change *ALL* the values below
	gCycleKeepTime							[6000],	'arrays sized to deal with the largest possible cycles
	gAveCtlContributors					[6000],
	gAveCtlAccum							[6000],
	gAveCtlMean								[6000],
	gAveEyeContributors					[6000],
	gAveEyeAccum							[6000],
	gAveEyeMean								[6000],
	gAveEyeVelContributors				[6000],
	gAveEyeVelAccum						[6000],
	gAveEyeVelMean							[6000],
	gAveCtlVelContributors				[6000],
	gAveCtlVelAccum						[6000],
	gAveCtlVelMean							[6000],
	gAveAux1Contributors					[6000],
	gAveAux1Accum							[6000],
	gAveAux1Mean							[6000],
	gAveAux2Contributors					[6000],
	gAveAux2Accum							[6000],
	gAveAux2Mean							[6000],
	gAveUnitContributors					[6000],
	gAccumulatedAcceptancePulses		[6000],
	gHistogram								[6000],
	'gSingleTraceSpikeTimes				[6000],	'the spike times in one cycle
	gSingleTraceDesaccadedSpikeTimes	[6000],	'the spike times in one cycle after desaccading
	gxAveArray								[6000],							
'*end*  all must have the same size

	gCtlPosFreq,
	gCtlPosAmp,
	gCtlPosPhase,
	gCtlPosOffset,

	
	gNumSpikes%,										'how many spikes in the cycle available to be accepted (desaccaded)
	gMaxAcceptedPts%,
	gLastNumDataPts%,
	gAveXYnum%,


	gSampleInterval,
	gSampleFreq,
	gHistoNumBin%						:= 256,							
	pDirection							:= 0,
	pPhase								:= 0,

	'below are the user definable variables
	gDataX1								:= 0,			'scrolling data window
	gDataY1								:= 0,
	gDataX2								:= 60,
	gDataY2								:= 100,

	gTrialX1								:= 60,		'control/eye data window
	gTrialY1								:= 0,
	gTrialX2								:= 100,
	gTrialY2								:= 25,

	gVelX1								:= 60,		'velocity window
	gVelY1								:= 25,
	gVelX2								:= 100,
	gVelY2								:= 50,

	gUnitX1								:= 60,		'unit data window
	gUnitY1								:= 50,
	gUnitX2								:= 100,
	gUnitY2								:= 60,

	gAux1X1								:= 60,		'aux1 data window
	gAux1Y1								:= 60,
	gAux1X2								:= 100,
	gAux1Y2								:= 70,

	gAux2X1								:= 60,		'aux1 data window
	gAux2Y1								:= 70,
	gAux2X2								:= 100,
	gAux2Y2								:= 80,

	gLogX1								:= 60,		'Log window
	gLogY1								:= 80,
	gLogX2								:= 100,
	gLogY2								:= 100,

	gAux1AveX1							:= 60,		'Ave aux1 data window
	gAux1AveY1							:= 50,
	gAux1AveX2							:= 100,
	gAux1AveY2							:= 60,

	gAux2AveX1							:= 60,		'Ave aux1 data window
	gAux2AveY1							:= 60,
	gAux2AveX2							:= 100,
	gAux2AVeY2							:= 70,

	gUnitBINX1							:= 5,			'Ave BIN data window
	gUnitBINY1							:= 5,
	gUnitBINX2							:= 55,
	gUnitBINY2							:= 50,

	gUnitISIX1							:= 5,			'Ave ISI data window
	gUnitISIY1							:= 50,
	gUnitISIX2							:= 55,
	gUnitISIY2							:= 95,

	gLogAveX1							:= 60,		'Ave Log window
	gLogAveY1							:= 70,
	gLogAveX2							:= 100,
	gLogAveY2							:= 100,

	
	gMemChEyeVel%,
	gMemChCtlVel%,
	gMemChSlipVel%,
	gMemChRectVel%,
	gMemChRectGaussSlipVel%,
	gMemChCtlVelSmooth%,
	gMemChSmoothAccel%,
	chUnitDS%,
	chSpike%,
	gLastDesaccadeIndex%,
	gChFile[10],	'more than large enough;
	gChUnit%,
	gChOne%,
	gChTwo%,
	gChThree%,
	gChFour%;


'----------------------------------------------------
'							Main
'6/17/05
'----------------------------------------------------

	gChFile[1]								:= 29;	'initial channel assignments
	gChFile[2]								:= 10;
	gChFile[3]								:= 2;
	gChFile[4]								:= 3;
	gChFile[5]								:= 4;
	gChUnit%									:= gChFile[1];	'suggested channel assignments
	gChOne%									:= gChFile[2];
	gChTwo%									:= gChFile[3];
	gChThree%								:= gChFile[4];
	gChFour%									:= gChFile[5];

	FileClose(-1,-1);

	SetToolbarOpenFile(0);							'Setup the toolbar
	DoToolbar();										'Setup the toolbar

	ToolbarVisible(1);								'Make toolbar visible always

	




'----------------------------------------------------
'						DoToolbar
'6/20/05
'----------------------------------------------------
proc DoToolbar()
	Toolbar("", 1023);	'Wait here until quit is pressed
end	'DoToolbar




'----------------------------------------------------
'						SetToolbarAnalysis
'6/20/05
'----------------------------------------------------
proc SetToolbarAnalysis(aMode%)
	var
		a%;
		
	ToolbarClear();												'eliminate all the buttons

	'set up the basic buttons
	ToolbarSet(0,"",Idle%);										'Call Idle%() whenever there is free time
	ToolbarSet(gQuitButton,"Quit",Quit%);
	ToolbarSet(gOpenButton,"Open",DialogOpen%);
	ToolbarSet(gFinalButton,"Perform Final Analysis",FinalAnalysis%);
	ToolbarSet(gBinButton,"Bins",DialogSetBins%);
	ToolbarSet(gAverageButton,"View Average Fit",ShowAverageViews%);
	ToolbarSet(gCurrentButton,"View Current Fit",ShowCurrentViews%);
	ToolbarSet(gResetButton,"Reset",ResetCycle%);
	ToolbarSet(gRejectButton,"&Reject",RejectCycle%);
	ToolbarSet(gAcceptButton,"&Accept",AcceptTrial%);
	ToolbarSet(gDesaccadeButton,"&Desaccade",Desaccade%);

	ToolbarEnable(gQuitButton,1);
	ToolbarEnable(gOpenButton,1);
	ToolbarEnable(gFinalButton,0);
	ToolbarEnable(gBinButton,0);
	ToolbarEnable(gAverageButton,0);
	ToolbarEnable(gCurrentButton,0);
	ToolbarEnable(gResetButton,0);
	ToolbarEnable(gRejectButton,0);
	ToolbarEnable(gAcceptButton,0);
	ToolbarEnable(gDesaccadeButton,0);
end	'SetToolbarAnalysis



'----------------------------------------------------
'					SetToolbarOpenFile
'1/6/06
'----------------------------------------------------
proc SetToolbarOpenFile(aMode%)
	var
		a%;
		
	ToolbarClear();												'eliminate all the buttons

	ToolbarSet(0,"",Idle%);										'Call Idle%() whenever there is free time
	ToolbarSet(1,"Quit",Quit%);
	ToolbarSet(2,"Open a File",DialogOpen%);
end	'SetToolbarOpenFile



'----------------------------------------------------
'					SetToolbarBeginAnalysis
'1/6/06
'----------------------------------------------------
proc SetToolbarBeginAnalysis(aMode%)
	var
		a%;
		
	ToolbarClear();												'eliminate all the buttons

	ToolbarSet(0,"",Idle%);										'Call Idle%() whenever there is free time
	ToolbarSet(1,"Quit",Quit%);
	ToolbarSet(2,"You are looking at the whole file. Position the cursors where you want to BEGIN and END the analysis. Change the X-axis scale if you need better resolution.  Then click here.",DialogAnalyze%);
end	'SetToolbarBeginAnalysis



'----------------------------------------------------
'							DialogSetBins
'8/31/05
'----------------------------------------------------
func DialogSetBins%()
	var
		ok%;

	dlgcreate("Select Number of Unit Bins");

	dlginteger(1,"Bins", 1, 5000);						'5000 bin max, want more?
	ok% := dlgShow(gHistoNumBin%);
	if ok% > 0 then
		BinAndFitUnit(gCyclesAccepted%);																			'writes to Ave Analysis window
		WriteAveValues();										'write the now values
	endif

	return kTrue;												'don't leave toolbar
end	'DialogSetBins



'----------------------------------------------------
'						DialogSetChannels
'6/20/05
'----------------------------------------------------
func DialogSetChannels%()
	var
		i%,
		j%,
		aDuplicateCh, aMissingCh, aProblem,
		aCh[kChannelsPlusOne],
		ok%,
		aRow, 
		aTab1		:= 30,
		list%[100];

	repeat
		for i% := 1 to kChannelsPlusOne-1 do
			aCh[i%] := gChFile[i%];
		next
	
		dlgcreate("Display Options");
		'dlgbutton(0,"");					'eliminate the Cancel button

		aRow := 1;
		dlginteger(1,"Unit  ( <1 means no unit)", -1, 99, aTab1, aRow);
		aRow += 1;
		dlginteger(2,"Control", 1, 99, aTab1, aRow);
		aRow += 1;
		dlginteger(3,"Eye", 1, 99, aTab1, aRow);
		aRow += 1;
		dlginteger(4,"Aux1", 1, 99, aTab1, aRow);
		aRow += 1;
		dlginteger(5,"Aux2", 1, 99, aTab1, aRow);
		aRow += 2;

		dlglist(6,"Display cycles","positive-going|negative-going", 99, aTab1, aRow);
		aRow += 1;
		dlglist(7,"Display cycles starting from","0 degrees|90 degrees", 99, aTab1, aRow);
	
		ok% := dlgShow(aCh[1],aCh[2],aCh[3],aCh[4],aCh[5],pDirection,pPhase);

		if ok% > 0 then
		
			aDuplicateCh := kFalse;
			for i% := 1 to kChannelsPlusOne-1 do
				for j% := 1 to kChannelsPlusOne-1 do
					if (i% <> j%) and (aCh[i%] = aCh[j%]) then
						aDuplicateCh := kTrue;
					endif
				next
			next
	
			aMissingCh := kFalse;
			for i% := 1 to kChannelsPlusOne-1 do
				if aCh[1] < 1 then
					'this is how a user specifies no unit trace
				else
					if (chankind(aCh[i%]) <= 0) then
						aMissingCh := kTrue;
					endif
				endif
			next
	
	
			aProblem := aMissingCh + aDuplicateCh;
	
			if (aProblem = kFalse) then
				for i% := 1 to kChannelsPlusOne-1 do
					gChFile[i%] := aCh[i%];
				next
				gChUnit%				:= gChFile[1];
				gChOne%				:= gChFile[2];
				gChTwo%				:= gChFile[3];
				gChThree%			:= gChFile[4];
				gChFour%				:= gChFile[5];
			else
				if aDuplicateCh = kTrue then
					Message("Some of the channels were duplicated so your changes were not implemented.");
				endif
				if aMissingCh = kTrue then
					Message("Some of the channels were missing so your changes were not implemented.");
				endif
			endif
		else
			halt;				'user clicked CANCEL
		endif
	until (aProblem = kFalse) or (ok% = 0);
	return kTrue;												'don't leave toolbar
end	'DialogSetChannels



'----------------------------------------------------
'								Quit
'6/20/05
'----------------------------------------------------
func Quit%()													'If "Quit" is pressed
	var 
		list%[100],
		i%;

	ViewList(list%[]);										'get a list of the views that match 
	for i% := 1 to list%[0] do
		if ViewKind(list%[i%]) >= 0 then
			View(list%[i%]);									'pull up those views
			FileClose(0, -1);									'close current, don't query
		endif
	next

	View(LogHandle());										'Make log view the current view
	EditSelectAll();											'Select all text in log view
	EditClear();												'Delete text

	gDatafilepath$ := "";
	return kFalse;												'leave toolbar
end	'Quit



'----------------------------------------------------
'						FinalAnalysis
'Ends the analysis, computes some averages and 
'presents an output document
'10/31/06
'----------------------------------------------------
func FinalAnalysis%()												'If "Final Analysis" is pressed
	var
		i%, j%, aStr$, aResult%, aResult$, aMean[37+1], aSD[37+1];


	arrconst(aMean,0);
	arrconst(aSD,0);

	'real frequency


	for j% := 0 to 37 do
		for i% := 1 to gCyclesAccepted% do
			aMean[j%] := aMean[j%] + gCycleIndFitData[i%][j%];
		next  'j%
		aMean[j%] := aMean[j%]/gCyclesAccepted%;
	next  'i%

	for j% := 0 to 37 do
		for i% := 1 to gCyclesAccepted% do
			aSD[j%] := aSD[j%] + pow(gCycleIndFitData[i%][j%] - aMean[j%],2) ;
		next  'j%
		aSD[j%] := sqrt(aSD[j%]/(gCyclesAccepted%-1));		'NOTE: use the (n-1) method like Excel's STDEV function
	next  'i%


	gHandleAnalysisFinal% := FileNew(1);						'create a new "LOG" type window for the average text information
	Window(gUnitBINX1,gUnitBINY1,gUnitISIX2,gUnitISIY2);	'Display it at the bottom of the screen
	WindowVisible(1);
	WindowTitle$("Final Analysis");


	aStr$ := "";
	aStr$ := aStr$ + gDatafileName$ + "\n";					'write file name
	aStr$ := aStr$ + str$(gCyclesAccepted%) + " of " + str$(gCyclesShown%-1) + " cycles accepted \n\n\n";	'don't count last cycle shown because it is neither accepted or rejected
	aStr$ := aStr$ + "LEAST SQUARES FIT OF INDIVIDUAL CYCLES (velocity trace):\n";
	aStr$ := aStr$ + "\t";
	aStr$ := aStr$ + print$("%8s","gain") 			+ "\t";		'gain
	aStr$ := aStr$ + print$("%8s","phase") 		+ "\t";		'phase
	aStr$ := aStr$ + print$("%8s","offset") 		+ "\t";		'offset
	aStr$ := aStr$ + print$("%8s","%pts left") 	+ "\t";		'pts remaining
	aStr$ := aStr$ + print$("%8s","num spks") 	+ "\t";		'spikes accepted
	aStr$ := aStr$ + "\n";

	for i% := 1 to gCyclesAccepted% do
		aStr$ := aStr$ + "\t";
		aStr$ := aStr$ + print$("%8.3f",gCycleIndFitData[i%][18]) 			+ "\t";	'gain
		aStr$ := aStr$ + print$("%8.3f",gCycleIndFitData[i%][14]/kRadians)+ "\t";	'phase
		aStr$ := aStr$ + print$("%8.3f",gCycleIndFitData[i%][16]) 			+ "\t";	'offset
		aStr$ := aStr$ + print$("%8.3f",gCycleIndFitData[i%][19]) 			+ "\t";	'pts remaining
		aStr$ := aStr$ + print$("%8.3f",gCycleIndFitData[i%][20]) 			+ "\t";	'spikes accepted
		aStr$ := aStr$ + "\n";
	next
	aStr$ := aStr$ + "\n";

	aStr$ := aStr$ + "mean" + "\t"; 
	aStr$ := aStr$ + print$("%8.3f",aMean[18]) 			+ "\t";	'mean gain
	aStr$ := aStr$ + print$("%8.3f",aMean[14]/kRadians)+ "\t";	'mean phase
	aStr$ := aStr$ + print$("%8.3f",aMean[16]) 			+ "\t";	'mean offset
	aStr$ := aStr$ + print$("%8.3f",aMean[19]) 			+ "\t";	'mean phase
	aStr$ := aStr$ + print$("%8.3f",aMean[20]) 			+ "\t";	'mean offset
	aStr$ := aStr$ + "\n";

	aStr$ := aStr$ + "SD" + "\t"; 
	aStr$ := aStr$ + print$("%8.3f",aSD[18]) + "\t";		'gain SD
	aStr$ := aStr$ + print$("%8.3f",aSD[14]/kRadians) + "\t";	'phase SD
	aStr$ := aStr$ + print$("%8.3f",aSD[16]) + "\t";		'offset SD
	aStr$ := aStr$ + "\n\n\n";

	aStr$ := aStr$ + "LEAST SQUARES FIT OF ALL ACCEPTED CYCLES:\n";
	aStr$ := aStr$ + "\t" + "     amp" + "\t" + "   phase" + "\t" + "  offset" + "\n";

	aStr$ := aStr$ + "Eye Pos" + "\t" 	+ print$("%8.3f",gCycleAvgFitData[30]) + "\t" + print$("%8.3f",gCycleAvgFitData[34]/kRadians)+ "\t" + print$("%8.3f",gCycleAvgFitData[36]) + "\n";
	aStr$ := aStr$ + "SD" 		+ "\t" 	+ print$("%8.3f",gCycleAvgFitData[31]) + "\t" + print$("%8.3f",gCycleAvgFitData[35]) 			+ "\t" + print$("%8.3f",gCycleAvgFitData[37]) + "\n\n";
	aStr$ := aStr$ + "Ctl Pos" + "\t"  	+ print$("%8.3f",gCycleAvgFitData[22]) + "\t" + print$("%8.3f",gCycleAvgFitData[26]/kRadians)+ "\t" + print$("%8.3f",gCycleAvgFitData[28]) + "\n";
	aStr$ := aStr$ + "SD" 		+ "\t" 	+ print$("%8.3f",gCycleAvgFitData[23]) + "\t" + print$("%8.3f",gCycleAvgFitData[27]) 			+ "\t" + print$("%8.3f",gCycleAvgFitData[29]) + "\n\n";
	aStr$ := aStr$ + "Eye Vel" + "\t"  	+ print$("%8.3f",gCycleAvgFitData[10]) + "\t" + print$("%8.3f",gCycleAvgFitData[14]/kRadians)+ "\t" + print$("%8.3f",gCycleAvgFitData[16]) + "\n";
	aStr$ := aStr$ + "SD" 		+ "\t" 	+ print$("%8.3f",gCycleAvgFitData[11]) + "\t" + print$("%8.3f",gCycleAvgFitData[15]) 			+ "\t" + print$("%8.3f",gCycleAvgFitData[17]) + "\n\n";
	aStr$ := aStr$ + "Ctl Vel"	+ "\t"  	+ print$("%8.3f",gCycleAvgFitData[02]) + "\t" + print$("%8.3f",gCycleAvgFitData[06]/kRadians)+ "\t" + print$("%8.3f",gCycleAvgFitData[08]) + "\n";
	aStr$ := aStr$ + "SD" 		+ "\t" 	+ print$("%8.3f",gCycleAvgFitData[03]) + "\t" + print$("%8.3f",gCycleAvgFitData[07]) 			+ "\t" + print$("%8.3f",gCycleAvgFitData[09]) + "\n\n";

	if gChUnit% > 0 then
		aStr$ := aStr$ + "Bin Mod" 	+ "\t"	+ print$("%8.3f",gCycleAvgFitData[38]) + "\t" + print$("%8.3f",gCycleAvgFitData[42]/kRadians)+ "\t" + print$("%8.3f",gCycleAvgFitData[44]) + "\n";
		aStr$ := aStr$ + "SD" 			+ "\t" 	+ print$("%8.3f",gCycleAvgFitData[39]) + "\t" + print$("%8.3f",gCycleAvgFitData[43]) 			+ "\t" + print$("%8.3f",gCycleAvgFitData[45]) + "\n\n";
		aStr$ := aStr$ + "ISI Mod" 	+ "\t" 	+ print$("%8.3f",gCycleAvgFitData[46]) + "\t" + print$("%8.3f",gCycleAvgFitData[50]/kRadians)+ "\t" + print$("%8.3f",gCycleAvgFitData[52]) + "\n";
		aStr$ := aStr$ + "SD" 			+ "\t" 	+ print$("%8.3f",gCycleAvgFitData[47]) + "\t" + print$("%8.3f",gCycleAvgFitData[51]) 			+ "\t" + print$("%8.3f",gCycleAvgFitData[53]) + "\n\n";
	endif
	aStr$ := aStr$ + "\n";

	aStr$ := aStr$ + "LEAST SQUARES FIT COMPARISONS:\n";
	aStr$ := aStr$ + "\t" + " amp/amp" + "\t" + " phs-phs" + "\n";
	aStr$ := aStr$ + "Eye Pos vs Ctl Pos:\n";
	aStr$ := aStr$ + "Fit" 			+ "\t" 	+ print$("%8.3f",gCycleAvgFitData[30]/gCycleAvgFitData[22])	+ "\t" + print$("%8.3f",(gCycleAvgFitData[34]-gCycleAvgFitData[26])/kRadians)	+ "\n\n";
	aStr$ := aStr$ + "Eye Vel vs Ctl Vel:\n";
	aStr$ := aStr$ + "Bin Fit" 	+ "\t" 	+ print$("%8.3f",gCycleAvgFitData[10]/gCycleAvgFitData[2]) 	+ "\t" + print$("%8.3f",(gCycleAvgFitData[14]-gCycleAvgFitData[6])/kRadians) 	+ "\n\n";

	if gChUnit% > 0 then
		aStr$ := aStr$ + "Unit Mod vs Ctl Pos:\n";
		aStr$ := aStr$ + "Bin Fit" 	+ "\t" 	+ print$("%8.3f",gCycleAvgFitData[38]/gCycleAvgFitData[22]) + "\t" + print$("%8.3f",(gCycleAvgFitData[42]-gCycleAvgFitData[26])/kRadians) 	+ "\n";
		aStr$ := aStr$ + "ISI Fit" 	+ "\t" 	+ print$("%8.3f",gCycleAvgFitData[46]/gCycleAvgFitData[22]) + "\t" + print$("%8.3f",(gCycleAvgFitData[50]-gCycleAvgFitData[26])/kRadians) 	+ "\n\n";
		aStr$ := aStr$ + "Unit Mod vs Ctl Vel:\n";
		aStr$ := aStr$ + "Bin Fit" 	+ "\t" 	+ print$("%8.3f",gCycleAvgFitData[38]/gCycleAvgFitData[2]) 	+ "\t" + print$("%8.3f",(gCycleAvgFitData[42]-gCycleAvgFitData[6])/kRadians) 		+ "\n";
		aStr$ := aStr$ + "ISI Fit" 	+ "\t" 	+ print$("%8.3f",gCycleAvgFitData[46]/gCycleAvgFitData[2]) 	+ "\t" + print$("%8.3f",(gCycleAvgFitData[50]-gCycleAvgFitData[6])/kRadians) 		+ "\n\n";
		aStr$ := aStr$ + "Unit Mod vs Eye Pos:\n";
		aStr$ := aStr$ + "Bin Fit" 	+ "\t" 	+ print$("%8.3f",gCycleAvgFitData[38]/gCycleAvgFitData[30]) + "\t" + print$("%8.3f",(gCycleAvgFitData[42]-gCycleAvgFitData[34])/kRadians) 	+ "\n";
		aStr$ := aStr$ + "ISI Fit" 	+ "\t" 	+ print$("%8.3f",gCycleAvgFitData[46]/gCycleAvgFitData[30]) + "\t" + print$("%8.3f",(gCycleAvgFitData[50]-gCycleAvgFitData[34])/kRadians) 	+ "\n\n";
		aStr$ := aStr$ + "Unit Mod vs Eye Vel:\n";
		aStr$ := aStr$ + "Bin Fit" 	+ "\t" 	+ print$("%8.3f",gCycleAvgFitData[38]/gCycleAvgFitData[10]) + "\t" + print$("%8.3f",(gCycleAvgFitData[42]-gCycleAvgFitData[14])/kRadians) 	+ "\n";
		aStr$ := aStr$ + "ISI Fit" 	+ "\t" 	+ print$("%8.3f",gCycleAvgFitData[46]/gCycleAvgFitData[10]) + "\t" + print$("%8.3f",(gCycleAvgFitData[50]-gCycleAvgFitData[14])/kRadians) 	+ "\n\n";
	endif

	print("%s", aStr$);	'fill it with a string
	aStr$ := "final_" + gDatafileName$ + "_" + print$("%.3f",gCycleIndFitData[1][0]) + "-" + print$("%.3f",gCycleIndFitData[gCyclesAccepted%][1]) + ".txt";	'beg,end
	aResult% 	:= filesaveas(aStr$,1,0);			'save it, ask to overwrite
	return kTrue;
end	'FinalAnalysis  



'----------------------------------------------------
'								Idle
'
'Control comes here when Spike2 has nothing to do.  
'7/7/05
'----------------------------------------------------
func Idle%()
	'nothing to do
	return kTrue;												'Stay in toolbar
end 'Idle



'----------------------------------------------------
'							DialogOpen
'Open the file and prepare traces		
'6/17/05
'----------------------------------------------------
func DialogOpen%();
	var
		i%,
		aResult%,
		aStr$, 
		bStr$,
		aResponse%,	
		aLittleBit,
		aCenter,
		hi, lo, 
		aType;


	if ViewKind(gViewData%) >= 0 then
		aResponse% := query("\nOpen a new file now and your current analysis will not be saved.\n", "Open a new file", "Cancel");
		if aResponse% = 0 then														'be sure the user really wants to quit
			return kTrue;
		else
			View(gViewData%);
			FileClose(-1,-1);
		endif
	endif

	gDurCycle			:= 0;
	gMaxAcceptedPts% 	:= 0;
	gCyclesShown% 		:= 0;				'no trials shown yet
	gCyclesAccepted%	:= 0;				'don't use [0]
	
	arrconst(gAveCtlContributors, 0);
	arrconst(gAveCtlAccum, 0);
	arrconst(gAveCtlMean,0);

	arrconst(gAveEyeContributors, 0);
	arrconst(gAveEyeAccum, 0);
	arrconst(gAveEyeMean,0);

	arrconst(gAveCtlVelContributors, 0);
	arrconst(gAveCtlVelAccum, 0);
	arrconst(gAveCtlVelMean,0);

	arrconst(gAveEyeVelContributors, 0);
	arrconst(gAveEyeVelAccum, 0);
	arrconst(gAveEyeVelMean,0);

	arrconst(gAveCtlContributors, 0);
	arrconst(gAveCtlAccum, 0);
	arrconst(gAveCtlMean,0);

	arrconst(gAveAux1Contributors, 0);
	arrconst(gAveAux1Accum, 0);
	arrconst(gAveAux1Mean,0);

	arrconst(gAveAux2Contributors, 0);
	arrconst(gAveAux2Accum, 0);
	arrconst(gAveAux2Mean,0);

	arrconst(gAveUnitContributors, 0);
	arrconst(gAccumulatedAcceptancePulses, 0);
	arrconst(gHistogram,0);


	arrconst(gCycleIndFitData,0);

	gViewData% := FileOpen("",0,1,"Select a file to analyze");

	if gViewData% < 0 then
		return gViewData%;
	else	
		SetToolbarBeginAnalysis(0);

		'create the scrolling data view
		'-----------------------------------------------------------------
		View(gViewData%);
		FrontView(gViewData%);
		gDatafilepath$ := FileName$(); 														'full file path
		gDatafileName$ := FileName$(3); 														'just the file name, no extension
		WindowVisible(1);
		Window(0,0,100,100);																		'fill the screen
		XRange(0,maxtime());
		HCursorDelete(-1);																		'delete all horizontal cursors
		CursorDelete(-1);																			'delete all vertical cursors

		aLittleBit := maxtime()/100;
		aCenter := maxtime()/2;
		aResult% := CursorNew(aCenter - aLittleBit, 1);									'position the cursor in the middle of the screen
		CursorLabel(4,1,"Move cursor to where you want to BEGIN the analysis");
		aResult% := CursorNew(aCenter + aLittleBit, 2);									'position the cursor in the middle of the screen
		CursorLabel(4,2,"Move cursor to where you want to END the analysis");
		CursorLabelPos(1,10);
		CursorLabelPos(2,12);


		for i% := 1 to 40 do
			aType := chankind(i%);
			if (aType=1) or (aType=6) or (aType=7) or (aType=9) then
				optimise(i%);				'optimize the trace
				hi := yHigh(i%);
				lo := yLow(i%);
				
				if hi < 10 then
					hi := 10;				'the HI scale should never be less than 10
				endif

				if lo > -10 then
					lo := -10				'the LO scale should never be more than -10
				endif
				yrange(i%, lo, hi);		'adjust the scales

			endif
		next
	endif
	return gViewData%;
end 'DialogOpen



'----------------------------------------------------
'						DialogAnalyze 
'Open the file and prepare traces		
'1/16/06
'----------------------------------------------------
func DialogAnalyze%();
	var
		i%,
		aResult%,
		aStr$, 
		bStr$;		


	gBegCycle 		:= cursor(1);			'user wants to start analysis here
	gEndCycle 		:= cursor(1);
	gStopAnalysis 	:= cursor(2);			'user wants to end analysis here
	gDurCycleLast	:= kNaN;

	gMaxAcceptedPts% := 0;
	arrconst(gAveCtlContributors, 0);
	arrconst(gAveCtlAccum, 0);
	arrconst(gAveCtlMean,0);

	arrconst(gAveEyeContributors, 0);
	arrconst(gAveEyeAccum, 0);
	arrconst(gAveEyeMean,0);

	arrconst(gAveEyeVelContributors, 0);
	arrconst(gAveEyeVelAccum, 0);
	arrconst(gAveEyeVelMean,0);

	arrconst(gAveCtlVelContributors, 0);
	arrconst(gAveCtlVelAccum, 0);
	arrconst(gAveCtlVelMean,0);

	arrconst(gAveAux1Contributors, 0);
	arrconst(gAveAux1Accum, 0);
	arrconst(gAveAux1Mean,0);

	arrconst(gAveAux2Contributors, 0);
	arrconst(gAveAux2Accum, 0);
	arrconst(gAveAux2Mean,0);

	arrconst(gAveUnitContributors, 0);
	arrconst(gAccumulatedAcceptancePulses, 0);
	arrconst(gHistogram,0);


	

	DialogSetChannels%();
	SetToolbarAnalysis(0);

	ToolbarEnable(gQuitButton,1);
	ToolbarEnable(gOpenButton,1);
	'ToolbarEnable(gBinButton,1);
	ToolbarEnable(gResetButton,1);
	ToolbarEnable(gRejectButton,1);
	ToolbarEnable(gAcceptButton,1);
	ToolbarEnable(gDesaccadeButton,1);

	'create the scrolling data view
	'-----------------------------------------------------------------
	View(gViewData%);
	FrontView(gViewData%);
	gDatafilepath$ := FileName$(); 											'full file path
	gDatafileName$ := FileName$(3); 											'just the file name, no extension
	WindowVisible(1);
	Window(gDataX1,gDataY1,gDataX2,gDataY2);
	XRange(0,1);
	HCursorDelete(-1);															'delete all horizontal cursors
	CursorDelete(-1);																'delete all vertical cursors

	ChanColour(gChOne%,1,16);
	ChanTitle$(gChOne%,"Control");
	ChanUnits$(gChOne%,"deg");
	YRange(gChOne%,-21,21);
	HcursorNew(gChOne%, 0);

	ChanColour(gChTwo%,1,1);
	ChanTitle$(gChTwo%,"Eye");
	ChanUnits$(gChTwo%,"deg");
	YRange(gChTwo%,-21,21);
	HcursorNew(gChTwo%, 0);

	ChanColour(gChThree%,1,38);
	ChanTitle$(gChThree%,"Aux1");
	ChanUnits$(gChThree%,"deg");
	YRange(gChThree%,-21,21);
	'HcursorNew(gChThree%, 0);

	ChanColour(gChFour%,1,32);
	ChanTitle$(gChFour%,"Aux2");
	ChanUnits$(gChFour%,"deg");
	YRange(gChFour%,-21,21);

	if chankind(gChUnit%) > 0 then
		ChanColour(gChUnit%,1,2);
		ChanTitle$(gChUnit%,"Unit");
		ChanUnits$(gChUnit%,"volts");
		YRange(gChUnit%,0,1);
		ChanWeight(gChUnit%, 0.25);
	endif

	'Compute eye velocity for desaccading
	gMemChEyeVel% := MemChan(9,0,BinSize(gChOne%)); 									'401
	ChanDiffRobi(gMemChEyeVel%, gChTwo%, 0, MaxTime());									
	ChanTitle$(gMemChEyeVel%,"EyeVel");
	YRange(gMemChEyeVel%,-250,250);
	ChanShow(gMemChEyeVel%);
	HcursorNew(gMemChEyeVel%, 0);

	'Compute control velocity for desaccading
	gMemChCtlVel% := MemChan(9,0,BinSize(gChOne%)); 									'402
	ChanDiffRobi(gMemChCtlVel%, gChOne%, 0, MaxTime());									
	ChanTitle$(gMemChCtlVel%,"ControlVel");
	YRange(gMemChCtlVel%,-250,250);

	'Compute rectified eye velocity for desaccading
	gMemChRectVel% := MemChan(9,0,BinSize(gChOne%)); 									'403
	ChanAbsRobi(gMemChRectVel%, gMemChEyeVel%, 0, MaxTime());					
	ChanTitle$(gMemChRectVel%,"RectEyeVel");
	YRange(gMemChRectVel%,-5,250);
	aResult%:=HcursorNew(gMemChRectVel%, 75);
	ChanShow(gMemChRectVel%);
	HcursorNew(gMemChRectVel%, 0);

	'Compute control-eye slip velocity for desaccading
	gMemChSlipVel% := MemChan(9,0,BinSize(gChOne%)); 									'404
	ChanSubRobi(gMemChSlipVel%, gMemChCtlVel%, gMemChEyeVel%, 0, MaxTime());	
	ChanTitle$(gMemChSlipVel%,"SlipVel");
	YRange(gMemChSlipVel%,-250,250);

	'smooth rectified control velocity for desaccading
	gMemChCtlVelSmooth% := MemChan(9,0,BinSize(gChOne%)); 							'405
	ChanSmoothRobi(gMemChCtlVelSmooth%, gMemChCtlVel%, 0, MaxTime());			
	ChanTitle$(gMemChCtlVelSmooth%,"SmoothCtlVel");
	YRange(gMemChCtlVelSmooth%,-250,250);

'	'smooth control-eye rectified acceleration for desaccading
'	gMemChSmoothAccel% := MemChan(9,0,BinSize(gChOne%)); 								'406
'	ChanDiffRobi(gMemChSmoothAccel%, gMemChCtlVelSmooth%, 0, MaxTime());			
'	ChanTitle$(gMemChSmoothAccel%,"SmoothSlipAccel");
'	YRange(gMemChSmoothAccel%,-250,250);
'
'		'Compute control velocity for desaccading
'		gMemChCtlVel% := MemChan(9,0,BinSize(gChOne%)); 								'401
'		ChanDiffRobi(gMemChCtlVel%, gChOne%, 0, MaxTime());									
'		ChanTitle$(gMemChCtlVel%,"EyeVel");
'		YRange(gMemChCtlVel%,-250,250);
'		ChanShow(gMemChCtlVel%);
'		HcursorNew(gMemChCtlVel%, 0);

	'-1=drop before  0=drop on  1=drop after
	if chankind(gChUnit%) > 0 then
		chSpike%	:= MemChan(5);																		'407  show acceptance pulses
		MemImport(chSpike%, gChUnit%, 0, MaxTime(), 2, 0, 0.5);							'positive threshold

		ChanOrder(gChOne%, 0, gChTwo%);
		ChanOrder(gChTwo%, 1, gChUnit%);
		ChanOrder(gChUnit%, 1, gMemChEyeVel%);
	endif
	ChanOrder(gMemChEyeVel%, 1, gMemChRectVel%);
	ChanOrder(gChThree%, 1, gChFour%);

	'create the to-be ACCEPTed views
	'-----------------------------------------------------------------
	gTrialView% := FileNew(12);													'XY view where cycles are accepted or rejected
	XYsetChan(0);
	WindowTitle$(gDatafileName$ + " - position");
	ChanTitle$(0, "Degrees");														'equivalent to Ytitle$ which doesn't exist
	Xtitle$("seconds");
	XYcolour(1,16);
	XYcolour(2,1);
	XYDrawMode(1,2,2);																'drawn as size 2 dots
	XYDrawMode(2,2,2);																'drawn as size 2 dots
	XYkey(kTrue,kTrue);																'show border
	XYkey(3,0);																			'no key border
	XYkey(2,1);																			'key background transparent
	ChanTitle$(1,"ctl");
	ChanTitle$(2,"eye");

	gVelocityView% := FileNew(12);												'XY view where cycles are accepted or rejected
	XYsetChan(0);
	WindowTitle$(gDatafileName$ + " - velocity");
	ChanTitle$(0, "Degrees");														'equivalent to Ytitle$ which doesn't exist
	Xtitle$("seconds");
	XYcolour(1,16);
	XYcolour(2,1);
	XYDrawMode(1,2,2);																'drawn as size 2 dots
	XYDrawMode(2,2,2);																'drawn as size 2 dots
	XYkey(kTrue,kTrue);																'show border
	XYkey(3,0);																			'no key border
	XYkey(2,1);																			'key background transparent
	ChanTitle$(1,"ctl vel");
	ChanTitle$(2,"eye vel");

	gAux1View% := FileNew(12);														'XY view where cycles are accepted or rejected
	XYsetChan(0);
	WindowTitle$(gDatafileName$ + " - aux1");
	ChanTitle$(0, "Degrees");														'equivalent to Ytitle$ which doesn't exist
	Xtitle$("seconds");
	XYcolour(1,38);
	XYDrawMode(1,2,2);																'drawn as size 2 dots

	gAux2View% := FileNew(12);														'XY view where cycles are accepted or rejected
	XYsetChan(0);
	WindowTitle$(gDatafileName$ + " - aux2");
	ChanTitle$(0, "Degrees");														'equivalent to Ytitle$ which doesn't exist
	Xtitle$("seconds");
	XYcolour(1,32);
	XYDrawMode(1,2,2);																'drawn as size 2 dots

	gUnitView% := FileNew(12);
	XYsetChan(0);																		'ch to hold the non-desaccaded ISI
	ChanHide(2);																		'users don't want to see it
	Window(gUnitX1,gUnitY1,gUnitX2,gUnitY2);
	WindowTitle$(gDatafileName$ + " - unit");
	WindowVisible(1);																	'show window
	XYcolour(1,1);
	XYDrawMode(1,1,8);   															'draw as vertical lines
	XYDrawMode(1,2,15);   															'make lines long

	gHandleAnalysis% := FileNew(1);												'create a new "LOG" type window for the average text information
	Window(gLogX1,gLogY1,gLogX2,gLogY2);										'Display it at the bottom of the screen
	WindowVisible(1);
	WindowTitle$("Analysis");

	
	'create the AVERAGE ACCEPTed views
	'-----------------------------------------------------------------
	gTrialAverageView% := FileNew(12);											'XY view where cycles are accepted or rejected
	XYsetChan(0);
	WindowTitle$(gDatafileName$ + " - AVERAGE position");
	ChanTitle$(0, "Degrees");														'equivalent to Ytitle$ which doesn't exist
	Xtitle$("seconds");
	XYcolour(1,16);
	XYcolour(2,1);
	XYDrawMode(1,2,2);																'drawn as size 2 dots
	XYDrawMode(2,2,2);																'drawn as size 2 dots
	XYkey(kTrue,kTrue);																'show border
	XYkey(3,0);																			'no key border
	XYkey(2,1);																			'key background transparent
	ChanTitle$(1,"ctl");
	ChanTitle$(2,"eye");

	gVelocityAverageView% := FileNew(12);										'XY view where cycles are accepted or rejected
	XYsetChan(0);
	WindowTitle$(gDatafileName$ + " - AVERAGE velocity");
	ChanTitle$(0, "Degrees");														'equivalent to Ytitle$ which doesn't exist
	Xtitle$("seconds");
	XYcolour(1,16);
	XYcolour(2,1);
	XYDrawMode(1,2,2);																'drawn as size 2 dots
	XYDrawMode(2,2,2);   															'Point size
	XYkey(kTrue,kTrue);																'show border
	XYkey(3,0);																			'no key border
	XYkey(2,1);																			'key background transparent
	ChanTitle$(1,"ctl vel");
	ChanTitle$(2,"eye vel");

	gAux1AverageView% := FileNew(12);											'XY view where cycles are accepted or rejected
	WindowTitle$(gDatafileName$ + " - AVERAGE aux1");
	ChanTitle$(0, "Degrees");														'equivalent to Ytitle$ which doesn't exist
	Xtitle$("seconds");
	XYcolour(1,38);
	XYDrawMode(1,2,2);																'drawn as size 2 dots

	gAux2AverageView% := FileNew(12);											'XY view where cycles are accepted or rejected
	XYsetChan(0);
	WindowTitle$(gDatafileName$ + " - AVERAGE aux2");
	ChanTitle$(0, "Degrees");														'equivalent to Ytitle$ which doesn't exist
	Xtitle$("seconds");
	XYcolour(1,32);
	XYDrawMode(1,2,2);																'drawn as size 2 dots

	if gChUnit% > 0 then
		gUnitHistoAverageView% := FileNew(12);
		WindowTitle$(gDatafileName$ + " - AVERAGE unit histo");
		ChanTitle$(0, "Spikes/Sec");													'equivalent to Ytitle$ which doesn't exist
		WindowVisible(1);																	'show window
		XYcolour(1,2);
		XYcolour(2,2);
		XYDrawMode(1,2,1);   															'dots
		XYJoin(1,1);
		XYDrawMode(1,2,2);																'drawn as size 2 dots
		XYDrawMode(2,2,1);   															'dots
	
	
		gUnitISIAverageView% := FileNew(12);
		WindowTitle$(gDatafileName$ + " - AVERAGE unit ISI");
		ChanTitle$(0, "Spikes/Sec");													'equivalent to Ytitle$ which doesn't exist
		WindowVisible(1);																	'show window
		XYcolour(1,2);
		XYcolour(2,2);
		XYDrawMode(1,2,1);   															'dots
		XYJoin(1,0);																		'don't join
		XYDrawMode(1,2,2);																'drawn as size 2 dots
		XYDrawMode(2,2,1);   															'dots
	endif

	gHandleAnalysisAve% := FileNew(1);											'create a new "LOG" type window for the average text information
	Window(gLogX1,gLogY1,gLogX2,gLogY2);										'Display it at the bottom of the screen
	WindowVisible(1);
	WindowTitle$("Average Analysis");




	'load up the ACCEPT views
	'-----------------------------------------------------------------
	ShowAverageViews%();
	ShowCurrentViews%();
	view(gViewData%);
	gSampleInterval := binsize(gChOne%);
	gSampleFreq := 1/gSampleInterval;

	for i% := 0 to gMaxCycleDur% - 1 do
		gxAveArray[i%] := i% * gSampleInterval;								'array of X values for desaccaded and non-desaccaded data
	next

	GetNextCycle();
	return gViewData%;
end 'DialogAnalyze



'----------------------------------------------------
'							ShowCurrent
'Shows the windows with the current plots		
'8/3/05
'----------------------------------------------------
func ShowCurrentViews%();

	'show the ACCEPTed views
	'-----------------------
	FrontView(gTrialView%);
	Window(gTrialX1,gTrialY1,gTrialX2,gTrialY2);
	FrontView(gVelocityView%);
	Window(gVelX1,gVelY1,gVelX2,gVelY2);
	FrontView(gAux1View%);
	Window(gAux1X1,gAux1Y1,gAux1X2,gAux1Y2);
	FrontView(gAux2View%);
	Window(gAux2X1,gAux2Y1,gAux2X2,gAux2Y2);
	FrontView(gUnitView%);
	Window(gUnitX1,gUnitY1,gUnitX2,gUnitY2);

	FrontView(gViewData%);
	FrontView(gHandleAnalysis%);
	Window(gLogX1,gLogY1,gLogX2,gLogY2);					'Display it at the bottom of the screen

	ToolbarEnable(gAverageButton,1);
	ToolbarEnable(gCurrentButton,0);

	ToolbarEnable(gBinButton,0);								'always turn off

	return kTrue;
end 'ShowCurrentViews



'----------------------------------------------------
'						ShowAverage
'Shows the windows with the average plots		
'8/3/05
'----------------------------------------------------
func ShowAverageViews%();

	'show the AVERAGE ACCEPTed views
	'--------------------------------
	FrontView(gTrialAverageView%);
	Window(gTrialX1,gTrialY1,gTrialX2,gTrialY2);
	FrontView(gVelocityAverageView%);
	Window(gVelX1,gVelY1,gVelX2,gVelY2);
	FrontView(gAux1AverageView%);
	Window(gAux1AveX1,gAux1AveY1,gAux1AveX2,gAux1AveY2);
	FrontView(gAux2AverageView%);
	Window(gAux2AveX1,gAux2AveY1,gAux2AveX2,gAux2AveY2);

	if gChUnit% > 0 then
		FrontView(gUnitISIAverageView%);
		Window(gUnitISIX1,gUnitISIY1,gUnitISIX2,gUnitISIY2);
		FrontView(gUnitHistoAverageView%);
		Window(gUnitBINX1,gUnitBINY1,gUnitBINX2,gUnitBINY2);
	endif

	FrontView(gHandleAnalysisAve%);
	Window(gLogAveX1,gLogAveY1,gLogAveX2,gLogAveY2);	'Display it at the bottom of the screen

	ToolbarEnable(gAverageButton,0);
	ToolbarEnable(gCurrentButton,1);

	if gCyclesAccepted% > 0 then								'can't bin until there are units to bin
		if gChUnit% > 0 then
			ToolbarEnable(gBinButton,1);
		endif
	endif

	return kTrue;
end 'ShowAverageViews



'----------------------------------------------------
'						DemarcateCycle	
'Moves vertical cursors to show the beginning and end
'of a cycle.  Used when a cycle is found at various
'times during analysis to move the cursor back in case
'the user moved it.  Someday CED will provide a cursor
'feature that will "lock" a cursor in place.
'
' 11/29/06
'----------------------------------------------------
proc DemarcateCycle();
	var
		aResult%;

	'display the ends of a cycle with vertical cursors #1 & #2
	View(gViewData%);
	aResult% := Cursor(1,gBegCycle);												'move cursor if one exists
	if aResult% <= 0 then
		CursorNew(gBegCycle, 1);													'if one does not exist then create a new cursor
		CursorLabel(4,1,"Beg Control");
	endif

	aResult% := Cursor(2, gEndCycle);											'move cursor if one exists
	if aResult% <= 0 then
		CursorNew(gEndCycle, 2);													'if one does not exist then create a new cursor
		CursorLabel(4,2,"End Control");
	endif
end	'DemarcateCycle



'----------------------------------------------------
'						GetNextCycle	
'Searches for the next cycle of data.  Records the 
'beginning, middle and and of the cycle then positions 
'the vertical cursors.  Also fills the time array used 
'for plotting the data in the XY plots.
'
' 6/27/05
'----------------------------------------------------
proc GetNextCycle()
	const
		kNoisyZeroCrossing	:= 0.100;	'skip 10ms to avoid any zero-crossings caused by noise

	var
		aResult%,
		aMargin,
		i%,
		n,
		aHalfCycle,
		aBeg,aEnd,
		aBegEndDir%,
		aMidCycle,
		aMidDir%,
		aPhase%,
		lt,rt,
		aDifference;


	ArrConst(gBegDsac[],kNAN);														'reset the desaccaded zone list
	ArrConst(gEndDsac[],kNAN);														'reset the desaccaded zone list
	gLastDesaccadeIndex% := 0;

	'find the next cycle


	if pPhase = 0 then
		aPhase% 				:= 0;		'in-phase

		if pDirection = 0 then
			aBegEndDir% 	:= 7;		'positive-going
			aMidDir%			:= 8;		'mid
		else
			aBegEndDir%		:= 8;		'negative-going
			aMidDir%			:= 7;		'mid
		endif
	else
		aPhase% 				:= 1;		'out-of-phase
		if pDirection = 0 then
			aBegEndDir% 	:= 8;		'negative-going
			aMidDir%			:= 7;		'mid
		else
			aBegEndDir%		:= 7;		'positive-going
			aMidDir%			:= 8;		'mid
		endif
	endif

	View(gViewData%);
	gBegCycle := ((gEndCycle-gBegCycle)*0.875)+gBegCycle;		'start searching from just before the end of the last cycle 

	gBegCycle := ChanSearch(gChOne%,aBegEndDir%,gBegCycle+kNoisyZeroCrossing,MaxTime(),0);	'find the next control position zero-crossing
	aMidCycle := ChanSearch(gChOne%,aMidDir%,gBegCycle+kNoisyZeroCrossing,MaxTime(),0);		'find the next control position zero-crossing
	gEndCycle := ChanSearch(gChOne%,aBegEndDir%,aMidCycle+kNoisyZeroCrossing,MaxTime(),0);	'find the subsequent control signal zero-crossing

	if gDurCycle = 0 then
		gDurCycle := gEndCycle - gBegCycle;		'all cycles must match the length of the first cycle
	else
		gEndCycle := gBegCycle + gDurCycle;		'compute uniform end of subsequent cycles
	endif

	if aPhase% > 0 then
		aHalfCycle := gDurCycle/2;
		aBeg := gBegCycle+aHalfCycle;
		aEnd := gEndCycle+aHalfCycle;
		'9.9.05  insert here the capability to find the next PEAKS or TROUGHS for out of phase presentations
		gBegCycle := ChanSearch(gChOne%,3,gBegCycle,aBeg);			'find the next control signal zero-crossing
		gEndCycle := ChanSearch(gChOne%,3,gEndCycle,aEnd);			'find the subsequent control signal zero-crossing
	endif


	if (gBegCycle >= 0) and (gEndCycle >= 0) then
	
		'test that we don't exceed storage arrays
		if gDurCycle/gSampleInterval > gMaxCycleDur% then
			Message("Cycles in file exceed array size allocated to handle them.  Have Bob increase gMaxCycleDur% if you want to proceed.");
		endif

		if gEndCycle > gStopAnalysis then
				Message("This cycle is beyond the cursor set as the end of analysis.  Consider clicking on Final Analysis button.");
		endif

		if gDurCycleLast <> kNaN then
			aDifference := ((gDurCycleLast-gDurCycle)/gDurCycleLast)*100;
			if (abs(aDifference) > 10) then
				Message("This cycle is ", aDifference, "% longer than the last cycle.");
			endif
		endif
		gDurCycleLast	:= gDurCycle;
		aMargin 			:= gDurCycle/10;
		Xrange(gBegCycle-aMargin, gEndCycle+aMargin);
	
	
		'load the time axis array
		for i% := 0 to gMaxCycleDur%-1 do
			n := gBegCycle + (i%*gSampleInterval);'0.001);						'fill the array with X values, it harmlessly goes beyond the cycle duration to the end of the array
	
			if n > gEndCycle then
				gCycleKeepTime[i%] := kNAN;											'NAN indicates unused data point, desaccade or end of array
			else
				gCycleKeepTime[i%] := n;												'save ms time of this data point
			endif
		next
	
	
		DemarcateCycle();
'		'display the ends of a cycle with vertical cursors #1 & #2
'		View(gViewData%);
'		aResult% := Cursor(1,gBegCycle);												'move cursor if one exists
'		if aResult% <= 0 then
'			CursorNew(gBegCycle, 1);													'otherwise make a new cursor
'			CursorLabel(4,1,"Beg Control");
'		endif
'	
'		aResult% := Cursor(2, gEndCycle);											'move cursor if one exists
'		if aResult% <= 0 then
'			CursorNew(gEndCycle, 2);													'otherwise make a new cursor
'			CursorLabel(4,2,"End Control");
'		endif
	
		
		'display the desaccading vertical cursors, #3 & #4
		View(gViewData%);
		gDurMS := gDurCycle*1000;
		gMidCycle := gBegCycle + gDurCycle/2;
		lt := gMidCycle - 0.005;
		rt := gMidCycle + 0.005;
	
		View(gViewData%);
		aResult% := Cursor(3,lt);														'move cursor if one exists
		if aResult% <= 0 then
			CursorNew(lt, 3);																'otherwise make a new cursor
			CursorLabel(4,3,"Desaccade");
			CursorLabelPos(3,10);
		endif
	
		aResult% := Cursor(4, rt);														'move cursor if one exists
		if aResult% <= 0 then
			CursorNew(rt, 4);																'otherwise make a new cursor
			CursorLabel(4,4,"Desaccade");
			CursorLabelPos(4,12);
		endif
	
		gCyclesShown% := gCyclesShown% + 1;											'another trial has been found and shown
		DrawSingleCycle();
	else
			Message("S T O P !   Only a partial cycle remains.");
	endif
end		'GetNextCycle



'----------------------------------------------------
'							ResetCycle	
'Lets the user start over again on a cycle
'
' 8/2/05
'----------------------------------------------------
func ResetCycle%()
	var
		aResult%,
		aMargin,
		i%,
		n,
		lt,rt;


	ArrConst(gBegDsac[],kNAN);														'reset the desaccaded zone list
	ArrConst(gEndDsac[],kNAN);														'reset the desaccaded zone list
	gLastDesaccadeIndex% := 0;


	'display the ends of a cycle with vertical cursors #1 & #2
	View(gViewData%);
	aResult% := Cursor(1,gBegCycle);												'move cursor if one exists
	if aResult% <= 0 then
		CursorNew(gBegCycle, 1);													'otherwise make a new cursor
	endif

	aResult% := Cursor(2, gEndCycle);											'move cursor if one exists
	if aResult% <= 0 then
		CursorNew(gEndCycle, 2);													'otherwise make a new cursor
	endif

	
	'display the desaccading vertical cursors, #3 & #4
	View(gViewData%);
	gDurMS := gDurCycle*1000;
	gMidCycle := gBegCycle + gDurCycle/2;
	lt := gMidCycle - 0.005;
	rt := gMidCycle + 0.005;

	View(gViewData%);
	aResult% := Cursor(3,lt);														'move cursor if one exists
	if aResult% <= 0 then
		CursorNew(lt, 3);																'otherwise make a new cursor
	endif

	aResult% := Cursor(4, rt);														'move cursor if one exists
	if aResult% <= 0 then
		CursorNew(rt, 4);																'otherwise make a new cursor
	endif

	DrawSingleCycle();
	return kTrue;																		'Stay in toolbar

end		'ResetCycle



'----------------------------------------------------
'						DrawSingleCycle	
'
'7/8/05
'----------------------------------------------------
proc DrawSingleCycle()
	var
		a%,b%,c%,i%,j%,
		aNext%,
		aNumDataPts%,
		aNumDataPts2use%,
		aEyeMean, aEyeSD,
		aVelMean, aVelSD,
		aCtlMean, aCtlSD,
		result,
		remaining%,
		err;



	aNumDataPts%  		:= gDurCycle*gSampleFreq;								'for making arrays
	gLastNumDataPts%	:= aNumDataPts%;
	aNumDataPts2use% 	:= 0;
	

	'only show one complete control cycle
	View(gTrialView%);
	XYdelete(1);																		'clear the old trace
	XYdelete(2);																		'clear the old trace	
	
	View(gVelocityView%);
	XYdelete(1);																		'clear the old trace
	XYdelete(2);																		'clear the old trace

	View(gAux1View%);
	XYdelete(1);																		'clear the old trace

	View(gAux2View%);
	XYdelete(1);																		'clear the old trace

	View(gUnitView%);
	XYdelete(1);																		'clear the old trace


	var
		aYdata[aNumDataPts%],														'just enough for one cycle, even the desaccaded areas
		aYsin[aNumDataPts%],
		aMask[aNumDataPts%];

	arrconst(aMask[],1);																'intialize the mask to 1 (use that point, indicating no desaccading)


'determine how much to desaccade
'------------------------------------
'two things come from here:
'   aMask[], selectively removes data
'   a%, num points saved, neded for creating result views
'
	'desaccade mask by zeroing each region that has been desaccaded
	for i% := 1 to gLastDesaccadeIndex% do
		for j% := gBegDsac[i%] to gEndDsac[i%] do
			aMask[j%] := 0;  'zero = don't use
		next
	next
	'count how many data points haven't been desaccaded
	a% := 0;
	for i% := 0 to aNumDataPts%-1 do 
		if aMask[i%] > 0 then
			a% += 1;
		endif
	next
	remaining% := a%;

	'create properly sized data arrays.
	var
		aYtemp[remaining%],
		aXtime[remaining%],
		aAmp, 
		aFreq, 
		aPhase, 
		aOff,
		aLSfit,
		aCovar[4][4],

		aCtlPosFreq,
		aCtlPosAmp,
		aCtlPosPhase,
		aCtlPosOffset,
		aCtlPosCovar[4][4],

		aCtlVelFreq,
		aCtlVelAmp,
		aCtlVelPhase,
		aCtlVelOffset,
		aCtlVelCovar[4][4],

		aEyePosFreq,
		aEyePosAmp,
		aEyePosPhase,
		aEyePosOffset,
		aEyePosCovar[4][4],

		aEyeVelFreq,
		aEyeVelAmp,
		aEyeVelPhase,
		aEyeVelOffset,
		aEyeVelCovar[4][4],

		aUnitFreq,
		aUnitAmp,
		aUnitPhase,
		aUnitOffset,
		aUnitCovar[4][4];


'display ctl cycle to be accepted
'------------------------------------
	View(gViewData%);
	arrconst(aYdata,0);
	chanData(gChOne%,aYdata,gBegCycle,gEndCycle);

	a% := 0;
	for i% := 0 to aNumDataPts%-1 do 
		if aMask[i%] > 0 then
			aYtemp[a%] := aYdata[i%];
			aXtime[a%] := gCycleKeepTime[i%];									'this really only needs to be done once per cycle
			a% += 1;
		endif
	next

	View(gTrialView%);
	XYaddData(kCtlTrial,aXtime[],aYtemp[]);									'add the data and spread it with a 1ms spread
	Xrange(gBegCycle, gEndCycle);													'make the plot as wide as the data		

	aFreq := kNaN;	'kNaN indicates this is the control position signal
	ComputeAfit(gTrialView%, kCtlTrial, gBegCycle, gEndCycle, aAmp, aFreq, aPhase, aOff, aLSfit, aCtlPosCovar[][]);
	aCtlPosFreq 	:= aFreq;
	aCtlPosAmp 		:= aAmp;
	aCtlPosPhase 	:= aPhase;
	aCtlPosOffset 	:= aOff;

	'save the control fit values to seed other fits
	gCtlPosFreq 	:= aFreq;
	gCtlPosAmp 		:= aAmp;
	gCtlPosPhase 	:= aPhase;
	gCtlPosOffset 	:= aOff;


'display eye cycle to be accepted
'-------------------------------------
	View(gViewData%);
	arrconst(aYdata,0);
	chanData(gChTwo%,aYdata,gBegCycle,gEndCycle);

	a% := 0;
	for i% := 0 to aNumDataPts%-1 do 
		if aMask[i%] > 0 then
			aYtemp[a%] := aYdata[i%];
			aXtime[a%] := gCycleKeepTime[i%];									'this really only needs to be done once per cycle
			a% += 1;
		endif
	next
	
	View(gTrialView%);
	XYaddData(kEyeTrial, aXtime[], aYtemp[]);									'add the data and spread it with a 1ms spread
	Xrange(gBegCycle, gEndCycle);													'make the plot as wide as the data		

	aFreq := aCtlPosFreq;
	ComputeAfit(gTrialView%, kEyeTrial, gBegCycle, gEndCycle, aAmp, aFreq, aPhase, aOff, aLSfit, aEyePosCovar[][]);
	aEyePosFreq		:= aFreq;
	aEyePosAmp 		:= aAmp;
	aEyePosPhase 	:= aPhase;
	aEyePosOffset 	:= aOff;

	optimise(-1);		'set the Y scale, include all traces


'display ctl velocity cycle to be accepted (11/20/06 use smoothed velocity)
'------------------------------------
	View(gViewData%);
	arrconst(aYdata,0);
	chanData(gMemChCtlVelSmooth%,aYdata,gBegCycle,gEndCycle);
	a% := 0;
	for i% := 0 to aNumDataPts%-1 do 
		if aMask[i%] > 0 then
			aYtemp[a%] := aYdata[i%];
			aXtime[a%] := gCycleKeepTime[i%];									'this really only needs to be done once per cycle
			a% += 1;
		endif
	next

	View(gVelocityView%);
	XYaddData(kCtlTrial,aXtime[],aYtemp[]);									'add the data and spread it with a 1ms spread
	Xrange(gBegCycle, gEndCycle);													'make the plot as wide as the data		

	aFreq := aCtlPosFreq;
	ComputeAfit(gVelocityView%, kCtlTrial, gBegCycle, gEndCycle, aAmp, aFreq, aPhase, aOff, aLSfit, aCtlVelCovar[][]);
	aCtlVelFreq 	:= aFreq;
	aCtlVelAmp 		:= aAmp;
	aCtlVelPhase 	:= aPhase;
	aCtlVelOffset 	:= aOff;

	optimise(-1);		'set the Y scale, include all traces


'display eye velocity cycle to be accepted
'------------------------------------
	View(gViewData%);
	arrconst(aYdata,0);
	chanData(gMemChEyeVel%,aYdata,gBegCycle,gEndCycle);

	a% := 0;
	for i% := 0 to aNumDataPts%-1 do 
		if aMask[i%] > 0 then
			aYtemp[a%] := aYdata[i%];
			aXtime[a%] := gCycleKeepTime[i%];									'this really only needs to be done once per cycle
			a% += 1;
		endif
	next

	View(gVelocityView%);
	XYaddData(kEyeTrial,aXtime[],aYtemp[]);									'add the data and spread it with a 1ms spread
	Xrange(gBegCycle, gEndCycle);													'make the plot as wide as the data		

	aFreq := aCtlPosFreq;
	ComputeAfit(gVelocityView%, kEyeTrial, gBegCycle, gEndCycle, aAmp, aFreq, aPhase, aOff, aLSfit, aEyeVelCovar[][]);
	aEyeVelFreq 	:= aFreq;
	aEyeVelAmp 		:= aAmp;
	aEyeVelPhase 	:= aPhase;
	aEyeVelOffset 	:= aOff;

	optimise(-1);		'set the Y scale, include all traces


'display aux1 cycle to be accepted
'-------------------------------------
	View(gViewData%);
	arrconst(aYdata,0);
	chanData(gChThree%, aYdata, gBegCycle, gEndCycle);

	a% := 0;
	for i% := 0 to aNumDataPts%-1 do 
		if aMask[i%] > 0 then
			aYtemp[a%] := aYdata[i%];
			aXtime[a%] := gCycleKeepTime[i%];									'this really only needs to be done once per cycle
			a% += 1;
		endif
	next

	View(gAux1View%);
	XYaddData(kAux1Trial,aXtime[],aYtemp[]);									'add the data and spread it with a 1ms spread
	Xrange(gBegCycle, gEndCycle);													'make the plot as wide as the data		

	optimise(-1);		'set the Y scale, include all traces


'display aux2 cycle to be accepted
'-------------------------------------
	View(gViewData%);
	arrconst(aYdata,0);
	chanData(gChFour%,aYdata,gBegCycle,gEndCycle);

	a% := 0;
	for i% := 0 to aNumDataPts%-1 do 
		if aMask[i%] > 0 then
			aYtemp[a%] := aYdata[i%];
			aXtime[a%] := gCycleKeepTime[i%];									'this really only needs to be done once per cycle
			a% += 1;
		endif
	next

	View(gAux2View%);
	XYaddData(kAux2Trial,aXtime[],aYtemp[]);									'add the data and spread it with a 1ms spread
	Xrange(gBegCycle, gEndCycle);													'make the plot as wide as the data		

	optimise(-1);		'set the Y scale, include all traces
	
	'if chankind(gChUnit%) > 0 then
	if gChUnit% > 0 then
		'display unit data to be accepted
		'-------------------------------------
		View(gViewData%);
		var
			aSpikeTimes		[aNumDataPts%],	'the spike times in one cycle
			aXbinData		[aNumDataPts%],	'the spike times in one cycle
			aYbinData		[aNumDataPts%],	'the spike times in one cycle
			aXisiData		[aNumDataPts%],	'the spike times in one cycle
			aYisiData		[aNumDataPts%],	'the spike times in one cycle
			aISIflag,
			aLastSpikeTime;
	
		
		gNumSpikes% := chanData(chSpike%,aSpikeTimes,gBegCycle,gEndCycle);		'get ALL the spike times from the main data window
		aISIflag := kTrue;																		'end spikes not needed for ISI computation
		a%	:=0;
		b% := 0;
		for j% := 0 to gNumSpikes%-1 do														'use each spike
			for i% := 0 to aNumDataPts%-1 do 												'determine which analog sample point it is associated with
				if (aSpikeTimes[j%] >= gCycleKeepTime[i%])  and  (aSpikeTimes[j%] < gCycleKeepTime[i%+1]) then
					if aMask[i%] > 0 then														'continue if not desaccaded
						aXbinData[b%]	:= aSpikeTimes[j%];
						aYbinData[b%]	:= 1;
						b% += 1;
						if aISIflag = kFalse then
							aXisiData[a%]	:= aSpikeTimes[j%];
							aYisiData[a%]	:= 1/(aSpikeTimes[j%]-aLastSpikeTime);		'compute instantaneous freq	'compute ISI
							a% += 1;
						else
							aISIflag := kFalse
						endif
						aLastSpikeTime	:= aSpikeTimes[j%];
		 			else
						aISIflag := kTrue;
					endif
				endif
			next
		next
		gNumSpikes% := b%;																		'hold number of spikes left after desaccading

		a% := a%-1;
		b% := b%-1;
		View(gUnitView%);
		xydelete(1);
		xydelete(2);
		XYaddData(kISItrial,aXisiData[0:a%],aYisiData[0:a%]);							'add the data with resolution of sampling rate
		XYaddData(kUnitTrial,aXbinData[0:b%],aYbinData[0:b%]);						'add the data with resolution of sampling rate
	endif

	Xrange(gBegCycle, gEndCycle);															'make the plot as wide as the data		
	YRange(1,-2,2);


	'these values will be written and re-written in the same space until the cycle is accepted
	gCycleIndFitHeading$[00]					:= "ctl begin time";					'names assigned here so list is as current as the var assignments
	gCycleIndFitHeading$[01]					:= "ctl end time";

	gCycleIndFitHeading$[02]					:= "ctl vel amp";
	gCycleIndFitHeading$[03]					:= "ctl vel amp SD";
	gCycleIndFitHeading$[04]					:= "ctl vel freq";
	gCycleIndFitHeading$[05]					:= "ctl vel freq SD";
	gCycleIndFitHeading$[06]					:= "ctl vel phase";
	gCycleIndFitHeading$[07]					:= "ctl vel phase SD";
	gCycleIndFitHeading$[08]					:= "ctl vel offset";
	gCycleIndFitHeading$[09]					:= "ctl vel offset SD";

	gCycleIndFitHeading$[10]					:= "eye vel amp";
	gCycleIndFitHeading$[11]					:= "eye vel amp SD";
	gCycleIndFitHeading$[12]					:= "eye vel freq";
	gCycleIndFitHeading$[13]					:= "eye vel freq SD";
	gCycleIndFitHeading$[14]					:= "eye vel phase";
	gCycleIndFitHeading$[15]					:= "eye vel phase SD";
	gCycleIndFitHeading$[16]					:= "eye vel offset";
	gCycleIndFitHeading$[17]					:= "eye vel offset SD";

	gCycleIndFitHeading$[18]					:= "gain";
	gCycleIndFitHeading$[19]					:= "% pts remaining";
	gCycleIndFitHeading$[20]					:= "num spikes";
	gCycleIndFitHeading$[21]					:= "num spikes SD";

	gCycleIndFitHeading$[22]					:= "ctl pos amp";
	gCycleIndFitHeading$[23]					:= "ctl pos amp SD";
	gCycleIndFitHeading$[24]					:= "ctl pos freq";
	gCycleIndFitHeading$[25]					:= "ctl pos freq SD";
	gCycleIndFitHeading$[26]					:= "ctl pos phase";
	gCycleIndFitHeading$[27]					:= "ctl pos phase SD";
	gCycleIndFitHeading$[28]					:= "ctl pos offset";
	gCycleIndFitHeading$[29]					:= "ctl pos offset SD";

	gCycleIndFitHeading$[30]					:= "eye pos amp";
	gCycleIndFitHeading$[31]					:= "eye pos amp SD";
	gCycleIndFitHeading$[32]					:= "eye pos freq";
	gCycleIndFitHeading$[33]					:= "eye pos freq SD";
	gCycleIndFitHeading$[34]					:= "eye pos phase";
	gCycleIndFitHeading$[35]					:= "eye pos phase SD";
	gCycleIndFitHeading$[36]					:= "eye pos offset";
	gCycleIndFitHeading$[37]					:= "eye pos offset SD";


	'write-over the same elements until one is accepted
	aNext% := gCyclesAccepted% + 1;

	'save eye analysis values final screen averaging
	gCycleIndFitData[aNext%][00] 	:= gBegCycle;						'time where the cycle started
	gCycleIndFitData[aNext%][01] 	:= gEndCycle;						'time where the cycle ended

	gCycleIndFitData[aNext%][02] 	:= aCtlVelAmp;						'Ctl amplitude
	gCycleIndFitData[aNext%][03] 	:= sqrt(aCtlVelCovar[0][0]);	'Ctl amplitude SD
	gCycleIndFitData[aNext%][04] 	:= aCtlVelFreq;					'Ctl frequency
	gCycleIndFitData[aNext%][05] 	:= sqrt(aCtlVelCovar[1][1]);	'Ctl frequency SD
	gCycleIndFitData[aNext%][06] 	:= aCtlVelPhase;					'Ctl phase
	gCycleIndFitData[aNext%][07] 	:= sqrt(aCtlVelCovar[2][2]);	'Ctl phase SD
	gCycleIndFitData[aNext%][08] 	:= aCtlVelOffset;					'Ctl offset
	gCycleIndFitData[aNext%][09] 	:= sqrt(aCtlVelCovar[3][3]);	'Ctl offset SD
	
	gCycleIndFitData[aNext%][10] 	:= aEyeVelAmp;						'eye amplitude
	gCycleIndFitData[aNext%][11] 	:= sqrt(aEyeVelCovar[0][0]);	'eye amplitude SD
	gCycleIndFitData[aNext%][12] 	:= aEyeVelFreq;					'eye frequency
	gCycleIndFitData[aNext%][13] 	:= sqrt(aEyeVelCovar[1][1]);	'eye frequency SD
	gCycleIndFitData[aNext%][14] 	:= aEyeVelPhase;					'eye phase
	gCycleIndFitData[aNext%][15] 	:= sqrt(aEyeVelCovar[2][2]);	'eye phase SD
	gCycleIndFitData[aNext%][16] 	:= aEyeVelOffset;					'eye offset
	gCycleIndFitData[aNext%][17] 	:= sqrt(aEyeVelCovar[3][3]);	'eye offset SD

	gCycleIndFitData[aNext%][18] 	:= aEyeVelAmp/aCtlVelAmp;		'gain
	gCycleIndFitData[aNext%][19] 	:= remaining% * 100 / aNumDataPts%;	'percent points left after desaccading

	gCycleIndFitData[aNext%][20] 	:= gNumSpikes%;					'number of spikes accepted
	gCycleIndFitData[aNext%][21] 	:= 0;									'number of spikes accepted SD

	gCycleIndFitData[aNext%][22] 	:= aCtlPosAmp;						'Ctl amplitude
	gCycleIndFitData[aNext%][23] 	:= sqrt(aCtlPosCovar[0][0]);	'Ctl amplitude SD
	gCycleIndFitData[aNext%][24] 	:= aCtlPosFreq;					'Ctl frequency
	gCycleIndFitData[aNext%][25] 	:= sqrt(aCtlPosCovar[1][1]);	'Ctl frequency SD
	gCycleIndFitData[aNext%][26] 	:= aCtlPosPhase;					'Ctl phase
	gCycleIndFitData[aNext%][27] 	:= sqrt(aCtlPosCovar[2][2]);	'Ctl phase SD
	gCycleIndFitData[aNext%][28] 	:= aCtlPosOffset;					'Ctl offset
	gCycleIndFitData[aNext%][29] 	:= sqrt(aCtlPosCovar[3][3]);	'Ctl offset SD
	
	gCycleIndFitData[aNext%][30] 	:= aEyePosAmp;						'eye amplitude
	gCycleIndFitData[aNext%][31] 	:= sqrt(aEyePosCovar[0][0]);	'eye amplitude SD
	gCycleIndFitData[aNext%][32] 	:= aEyePosFreq;					'eye frequency
	gCycleIndFitData[aNext%][33] 	:= sqrt(aEyePosCovar[1][1]);	'eye frequency SD
	gCycleIndFitData[aNext%][34] 	:= aEyePosPhase;					'eye phase
	gCycleIndFitData[aNext%][35] 	:= sqrt(aEyePosCovar[2][2]);	'eye phase SD
	gCycleIndFitData[aNext%][36] 	:= aEyePosOffset;					'eye offset
	gCycleIndFitData[aNext%][37] 	:= sqrt(aEyePosCovar[3][3]);	'eye offset SD

	'gPercentLeft := gCycleIndFitData[aNext%][19]

	'update cycle info in Log window
	View(gHandleAnalysis%);																	'Make log view the current view
	EditSelectAll();																			'Select all text in log view
	EditClear();																				'Delete it
	'FontSet(1, 9, 0);	'system font, 7pt, normal text
	Print("cycle #"+ str$(gCyclesShown%) + "= " + print$("%5.3f", 1/gDurCycle) + " Hz, " + print$("%5.3f", gDurCycle) + " sec long, " + print$("%5.0f", remaining% * 100 / aNumDataPts%) + "%% data remaining after desaccading \n");

	Print("Fit Data               Amp       Freq      Phase     Offset \n");
	Print("EyePos Values:  %10.3f %10.3f %10.3f %10.3f \n", 	aEyePosAmp, aEyePosFreq/k2pi, aEyePosPhase*kRadians, aEyePosOffset);
	Print("EyePos StdDev:  %10.3f %10.3f %10.3f %10.3f \n", 	sqrt(aEyePosCovar[0][0]), sqrt(aEyePosCovar[1][1]), sqrt(aEyePosCovar[2][2]), sqrt(aEyePosCovar[3][3]));
	Print("EyeVel Values:  %10.3f %10.3f %10.3f %10.3f \n", 	aEyeVelAmp, aEyeVelFreq/k2pi, aEyeVelPhase*kRadians, aEyeVelOffset);
	Print("EyeVel StdDev:  %10.3f %10.3f %10.3f %10.3f \n", 	sqrt(aEyeVelCovar[0][0]), sqrt(aEyeVelCovar[1][1]), sqrt(aEyeVelCovar[2][2]), sqrt(aEyeVelCovar[3][3]));
	Print("CtlPos Values:  %10.3f %10.3f %10.3f %10.3f \n", 	aCtlPosAmp, aCtlPosFreq/k2pi, aCtlPosPhase*kRadians, aCtlPosOffset);
	Print("CtlPos StdDev:  %10.3f %10.3f %10.3f %10.3f \n", 	sqrt(aCtlPosCovar[0][0]), sqrt(aCtlPosCovar[1][1]), sqrt(aCtlPosCovar[2][2]), sqrt(aCtlPosCovar[3][3]));
	Print("CtlVel Values:  %10.3f %10.3f %10.3f %10.3f \n", 	aCtlVelAmp, aCtlVelFreq/k2pi, aCtlVelPhase*kRadians, aCtlVelOffset);
	Print("CtlVel StdDev:  %10.3f %10.3f %10.3f %10.3f", 		sqrt(aCtlVelCovar[0][0]), sqrt(aCtlVelCovar[1][1]), sqrt(aCtlVelCovar[2][2]), sqrt(aCtlVelCovar[3][3]));

	'standard deviation is defined as the square root of the variance. This means it is the root mean square (RMS) 
	'deviation from the average. It is defined this way in order to give us a measure of dispersion that is (1) a 
	'non-negative number, and (2) has the same units as the data.
	FrontView(gViewData%);																'make the main window active
end		'DrawSingleCycle
	

	
'----------------------------------------------------
'						ComputeAfit
' Calls Spike2's ChanFit() routine to get a sin fit.  
' Frequency = NaN signals the fit is being made to a 
' control position signal to which all other signals 
' must use the same frequency value.
'		
'10/24/06
'----------------------------------------------------
func ComputeAfit(aView%, aCh, aBeg, aEnd, &aAmp, &aFreq, &aPhase, &aOff, aLSfit, &aCovar[][])

	'NOTE: aFreq & aPhase are in RADIANS
	
	var
		aHold%[4],
		aIter,
		aResult,
		aFreqMax;

	arrconst(aHold%, 0);											'set to no holds

	View(aView%);
	ChanFit		(aCh,	4,	1);   								'Initialize for sine fit

	'coefficients:
	'  0 = amplitude	in my units -- degrees
	'  1 = frequency	in Radians
	'  2 = phase		in Radians
	'  3 = offset		in my units -- degrees

	'12/4/06  Greg Smith at CED advised me not to limit fits unless absolutely necessary
	'as it may prevent the process from passing through the limited region and prevent
	'finding the best fit.
	'ChanFitCoef	(aCh,	0,	10,	-100,	100); 			'coeff #0:  +/-100, amplitude
	'ChanFitCoef	(aCh,	2,	0,		-kPi,	kPi);   			'coeff #2, phase:   +/-3.14  or +/-180deg
	'ChanFitCoef	(aCh,	3,	0,		-100,	100);				'coeff #3, offset:  +/-100deg, offset
	if aFreq = kNaN then											'signals this fit for control sig, use seed freq, but let it find best freq
		aFreq 	:= (1/gDurCycle) * k2Pi;					'seed the fitting process with the zero-crossing computed frequency
		aFreqMax := 10 * k2Pi;									'10 Hz in radians
		ChanFitCoef	(aCh,	1,	aFreq,	0,		aFreqMax);	'coeff #1, frequency:  0Hz->aFreq->10Hz, frequency in Radians -> div by 2Pi to get Hz
		aHold%[1] := 0;											'set freq to float
	else
		ChanFitCoef	(aCh,	1,	aFreq,	aFreq,	aFreq); 	'coeff #1:  lock it to the control frequency passed to this routine (seed,lo,hi)
		aHold%[1] := 1;											'set freq to HOLD (perhaps redundant, prob better than above)
	endif
	ChanFitCoef (aCh, aHold%[]);
	ChanFitShow	(aCh,	1);   									'Set up fit display

	aResult := ChanFit		(aCh,	3,	aBeg,	aEnd, aBeg, aLSfit, 100, aIter, aCovar[][]);			'Fit to data, reference beginning of data
	colour(35,8);

	aAmp 		:= ChanFitCoef	(aCh, 0);
	aFreq 	:= ChanFitCoef	(aCh, 1);						'frequency in Radians -> div by 2Pi to get Hz
	aPhase	:= ChanFitCoef	(aCh, 2);
	aOff		:= ChanFitCoef	(aCh, 3);
	
	if (aPhase > kPi) or (aPhase < -kPi) or (aResult <> 0)then
		aIter := aIter;
	endif

	'page 8-3 in v5.12 of "The Script Language"
'	aSD := sqrt(covar[0][0]);'[0]=amplitude  [1]=freq  [2]=phase

	return kTrue;
end	'ComputeAfit



'----------------------------------------------------
'							Desaccade
'Responds when user clicks the Desaccade button.  The 
'desaccaded regions (to-from) are stored in an array.
'		
'6/20/05
'----------------------------------------------------
func Desaccade%()
	var
		aBegSac, aEndSac, aTmpSac,
		aBegCycle, aEndCycle,
		i;

	i := gSampleInterval;

	DemarcateCycle();						'draw the vertical cursors showing the beginning and end of the cycle.
	View(gViewData%);
	'aBegCycle	:= cursor(1);
	'aEndCycle	:= cursor(2);
	aBegCycle	:= gBegCycle;			'11/29/06  can't use cursor(1) and cursor(2) because user can move them
	aEndCycle	:= gEndCycle;
	aBegSac		:= cursor(3);																	'get the cursor locations
	aEndSac		:= cursor(4);

	if aBegSac > aEndSac then																	'make sure cursor(3) points to the earliest time
		aTmpSac := aEndSac;																		'if not, reverse them
		aEndSac := aBegSac;
		aBegSac := aTmpSac;
	endif

	'8/18/05  I'm sure this will come up again . . . if the last point is desaccaded then there 
	'is an index too large error in DrawSingleCycle().  If the first point is desaccaded then
	'the fit traces are shifted too far to the left.  Until this is understood, just prevent these
	'numbers from coming into play.
	'
	'11/29/06  This means an area desaccaded from the beginning of a trace will have the first original point
	'even if it was desaccaded out.  This also means that an area desaccaded from the end of a trace will have 
	'the last original point even if it was desaccade out.  Waiting for good alternative.
	if aBegSac <= aBegCycle then																'no desaccading before the cycle
		aBegSac := aBegCycle+i;
	endif
	if aEndSac < aBegCycle then
		aEndSac := aBegCycle+i;
	endif

	if aBegSac >= aEndCycle then																'no desaccading after the cycle
		aBegSac := aEndCycle;
	endif
	if aEndSac >= aEndCycle then
		aEndSac := aEndCycle-i;
	endif

	gLastDesaccadeIndex% += 1;																	'keep track of how many areas have been desaccaded
	gBegDsac[gLastDesaccadeIndex%] := trunc((aBegSac-gBegCycle)*gSampleFreq);	'store the beginning of the desaccaded region
	gEndDsac[gLastDesaccadeIndex%] := trunc((aEndSac-gBegCycle)*gSampleFreq);	'store the end of the desaccaded region

	DrawSingleCycle();																			'show the desaccaded cycle
	return(kTrue);
end		'desaccade



'----------------------------------------------------
' 							RejectCycle
'		
' 8/17/05
'----------------------------------------------------
func RejectCycle%()

	GetNextCycle();	'get another cycle
	return(kTrue);
end	'RejectCycle



'----------------------------------------------------
' 							Accept
'		
' 6/27/05
'----------------------------------------------------
func AcceptTrial%()

	var
		aXYnum%,
		i%,
		aX, aY,
		aCtlMean, aCtlSD,
		aEyeMean, aEyeSD,
		aVelMean, aVelSD,
		aCtlCoefs[4],
		aEyeCoefs[4],
		aVelCoefs[4],

		aAmp, 
		aFreq, 
		aPhase, 
		aOff,		
		aLSfit,

		aCtlPosFreq,
		aCtlPosAmp,
		aCtlPosPhase,
		aCtlPosOffset,
		aCtlPosCovar[4][4],

		aEyePosFreq,
		aEyePosAmp,
		aEyePosPhase,
		aEyePosOffset,
		aEyePosCovar[4][4],

		aCtlVelFreq,
		aCtlVelAmp,
		aCtlVelPhase,
		aCtlVelOffset,
		aCtlVelCovar[4][4],

		aEyeVelFreq,
		aEyeVelAmp,
		aEyeVelPhase,
		aEyeVelOffset,
		aEyeVelCovar[4][4],

		aISIFreq,
		aISIAmp,
		aISIPhase,
		aISIOffset,
		aISICovar[4][4],

		aUnitFreq,	'is this being used?
		aUnitAmp,
		aUnitPhase,
		aUnitOffset,
		aUnitCovar[4][4];

	gCyclesAccepted% := gCyclesAccepted% + 1;											'a cycle has been accepted, remember that[zero] is unused
	if gLastNumDataPts% > gMaxAcceptedPts% then										'must accomodate the largest cycle, shouldn't vary by more than a few data points
		gMaxAcceptedPts% := gLastNumDataPts%;
	endif

	var
		xArray		[gMaxAcceptedPts%],													'array size is the width of the cycle and is not affected by desaccading
		yArray		[gMaxAcceptedPts%],													'array size based on cycle duration and sample frequency
		aYsin			[gMaxAcceptedPts%];


	View(gTrialAverageView%);
	XYdelete(1);																				'clear the old trace
	XYdelete(2);																				'clear the old trace	

	View(gVelocityAverageView%);
	XYdelete(1);																				'clear the old trace
	XYdelete(2);																				'clear the old trace

	View(gAux1AverageView%);
	XYdelete(1);																				'clear the old trace

	View(gAux2AverageView%);
	XYdelete(1);																				'clear the old trace

	if gChUnit% > 0 then
		View(gUnitHistoAverageView%);
		XYdelete(1);																			'clear the old trace
		XYdelete(2);																			'clear the old trace
	
		View(gUnitISIAverageView%);
	endif

'get the acccepted control trace and add it to the previously accepted control traces

'display the average accepted CONTROL cycles
'-------------------------------------------
	View(gTrialView%);
	aXYnum% := XYGetData(kCtlTrial, xArray[], yArray[]);							'get the X & Y control data
 	AccumulateAve(aXYnum%, xArray[], yArray[], gAveCtlMean[], gAveCtlAccum[], gAveCtlContributors[],xArray[0]);

	gAveXYnum% := SizeAveArrays(gAveCtlContributors[]);							'all traces are desaccaded the same way so they all have the same number of elements
	var
		xAveArrayDS	[gAveXYnum%], 															'needs to be dimensioned to the size of the array - desaccades
		yAveArrayDS	[gAveXYnum%];
	LoadAveArrays(xAveArrayDS[], yAveArrayDS[], gAveCtlMean[], gAveCtlContributors[]);

	View(gTrialAverageView%);
	XYaddData(kCtlTrial,xAveArrayDS[],yAveArrayDS[]);								'display the ave control data

	aFreq := kNaN;	'kNaN indicates this is the control position signal
	ComputeAfit(gTrialAverageView%, kCtlTrial, 0, gEndCycle-gBegCycle, aAmp, aFreq, aPhase, aOff, aLSfit, aCtlPosCovar[][]);
	aCtlPosFreq 	:= aFreq;
	aCtlPosAmp 		:= aAmp;
	aCtlPosPhase 	:= aPhase;
	aCtlPosOffset 	:= aOff;
	optimise(-1);		'set the Y scale, include all traces


'display the average accepted EYE cycles
'-------------------------------------------
	View(gTrialView%);
	aXYnum% := XYGetData(kEyeTrial, xArray[], yArray[]);							'get the X & Y eye data
 	AccumulateAve(aXYnum%, xArray[], yArray[], gAveEyeMean[], gAveEyeAccum[], gAveEyeContributors[], xArray[0]);
	LoadAveArrays(xAveArrayDS[], yAveArrayDS[], gAveEyeMean[], gAveEyeContributors[]);

	View(gTrialAverageView%);
	XYaddData(kEyeTrial,xAveArrayDS[],yAveArrayDS[]);								'display the ave control data

	aFreq := aCtlPosFreq;
	ComputeAfit(gTrialAverageView%, kEyeTrial, 0, gEndCycle-gBegCycle, aAmp, aFreq, aPhase, aOff, aLSfit, aEyePosCovar[][]);
	aEyePosFreq 	:= aFreq;
	aEyePosAmp 		:= aAmp;
	aEyePosPhase 	:= aPhase;
	aEyePosOffset 	:= aOff;
	optimise(-1);		'set the Y scale, include all traces


'display the average accepted CTL VEL cycles
'-------------------------------------------
	View(gVelocityView%);
	aXYnum% := XYGetData(kCtlVelTrial, xArray[], yArray[]);					'get the X & Y vel data
 	AccumulateAve(aXYnum%, xArray[], yArray[], gAveCtlVelMean[], gAveCtlVelAccum[], gAveCtlVelContributors[], xArray[0]);
	LoadAveArrays(xAveArrayDS[], yAveArrayDS[], gAveCtlVelMean[], gAveCtlVelContributors[]);

	View(gVelocityAverageView%);
	XYaddData(kCtlVelTrial,xAveArrayDS[],yAveArrayDS[]);						'display the ave vel data

	aFreq := aCtlPosFreq;
	ComputeAfit(gVelocityAverageView%, kCtlVelTrial, 0, gEndCycle-gBegCycle, aAmp, aFreq, aPhase, aOff, aLSfit, aCtlVelCovar[][]);
	aCtlVelFreq 	:= aFreq;
	aCtlVelAmp 		:= aAmp;
	aCtlVelPhase 	:= aPhase;
	aCtlVelOffset 	:= aOff;
	optimise(-1);		'set the Y scale, include all traces
	
	
'display the average accepted EYE VEL cycles
'-------------------------------------------
	View(gVelocityView%);
	aXYnum% := XYGetData(kEyeVelTrial, xArray[], yArray[]);					'get the X & Y vel data
 	AccumulateAve(aXYnum%, xArray[], yArray[], gAveEyeVelMean[], gAveEyeVelAccum[], gAveEyeVelContributors[], xArray[0]);
	LoadAveArrays(xAveArrayDS[], yAveArrayDS[], gAveEyeVelMean[], gAveEyeVelContributors[]);

	View(gVelocityAverageView%);
	XYaddData(kEyeVelTrial,xAveArrayDS[],yAveArrayDS[]);						'display the ave vel data

	aFreq := aCtlPosFreq;
	ComputeAfit(gVelocityAverageView%, kEyeVelTrial, 0, gEndCycle-gBegCycle, aAmp, aFreq, aPhase, aOff, aLSfit, aEyeVelCovar[][]);
	aEyeVelFreq 	:= aFreq;
	aEyeVelAmp 		:= aAmp;
	aEyeVelPhase 	:= aPhase;
	aEyeVelOffset 	:= aOff;
	optimise(-1);		'set the Y scale, include all traces
	
	
'display the average accepted AUX1 cycles
'-------------------------------------------
	View(gAux1View%);
	aXYnum% := XYGetData(kAux1Trial, xArray[], yArray[]);							'get the X & Y Aux1 data
 	AccumulateAve(aXYnum%, xArray[], yArray[], gAveAux1Mean[], gAveAux1Accum[], gAveAux1Contributors[], xArray[0]);
	LoadAveArrays(xAveArrayDS[], yAveArrayDS[], gAveAux1Mean[], gAveAux1Contributors[]);

	View(gAux1AverageView%);
	XYaddData(kAux1Trial,xAveArrayDS[],yAveArrayDS[]);								'display the ave Aux1 data
	optimise(-1);		'set the Y scale, include all traces


'display the average accepted AUX2 cycles
'-------------------------------------------
	View(gAux2View%);
	aXYnum% := XYGetData(kAux2Trial, xArray[], yArray[]);							'get the X & Y Aux2 data
 	AccumulateAve(aXYnum%, xArray[], yArray[], gAveAux2Mean[], gAveAux2Accum[], gAveAux2Contributors[], xArray[0]);
	LoadAveArrays(xAveArrayDS[], yAveArrayDS[], gAveAux2Mean[], gAveAux2Contributors[]);

	View(gAux2AverageView%);
	XYaddData(kAux2Trial,xAveArrayDS[],yAveArrayDS[]);								'display the ave Aux2 data
	optimise(-1);		'set the Y scale, include all traces


	if gChUnit% > 0 then
		'display the average accepted unit activity
		'-------------------------------------------
		'	'ISI analysis
		ArrConst(xArray[],0);
		ArrConst(yArray[],0);
		ArrConst(gHistogram[],0); 
		ArrConst(gAccumulatedAcceptancePulses[],0);
		ArrConst(gAveUnitContributors[],0);
		View(gUnitView%);
		aXYnum% := XYGetData(kISItrial, xArray[], yArray[]);							'get the ISI Unit data
	 	AccumulateAve(aXYnum%, xArray[], yArray[], gHistogram[], gAccumulatedAcceptancePulses[], gAveUnitContributors[],xLow());
	
		View(gUnitISIAverageView%);
		XYaddData(1,xArray[], yArray[]);														'save it in the average ISI Unit XY-plot
		Xrange(0,gxAveArray[gMaxAcceptedPts%]);											'limit the x axis
	
		aFreq := aCtlPosFreq;
		ComputeAfit(gUnitISIAverageView%, 1, 0, gEndCycle-gBegCycle, aAmp, aFreq, aPhase, aOff, aLSfit, aISICovar[][]);
		aISIFreq 	:= aFreq;
		aISIAmp 		:= aAmp;
		aISIPhase 	:= aPhase;
		aISIOffset 	:= aOff;
		optimise(-1);		'set the Y scale, include all traces
	
		'BIN analysis
		ArrConst(xArray[],0);
		ArrConst(yArray[],0);
		ArrConst(gHistogram[],0); 
		ArrConst(gAccumulatedAcceptancePulses[],0);
		ArrConst(gAveUnitContributors[],0);
		View(gUnitView%);
		aXYnum% := XYGetData(kUnitTrial, xArray[], yArray[]);							'get the Bin Unit data
	 	AccumulateAve(aXYnum%, xArray[], yArray[], gHistogram[], gAccumulatedAcceptancePulses[], gAveUnitContributors[], xLow());
		BinAndFitUnit(gCyclesAccepted%);														'writes to Ave Analysis window
		'must be subroutine so Bins can be changed by a button click
	endif

	'save eye analysis values final screen averaging
	'gCycleAvgFitData[00] 	:= gBegCycle;						'time where the cycle started
	'gCycleAvgFitData[01] 	:= gEndCycle;						'time where the cycle ended

	gCycleAvgFitData[02] 	:= aCtlVelAmp;						'Ctl amplitude
	gCycleAvgFitData[03] 	:= sqrt(aCtlVelCovar[0][0]);	'Ctl amplitude SD
	gCycleAvgFitData[04] 	:= aCtlVelFreq;					'Ctl frequency
	gCycleAvgFitData[05] 	:= sqrt(aCtlVelCovar[1][1]);	'Ctl frequency SD
	gCycleAvgFitData[06] 	:= aCtlVelPhase;					'Ctl phase
	gCycleAvgFitData[07] 	:= sqrt(aCtlVelCovar[2][2]);	'Ctl phase SD
	gCycleAvgFitData[08] 	:= aCtlVelOffset;					'Ctl offset
	gCycleAvgFitData[09] 	:= sqrt(aCtlVelCovar[3][3]);	'Ctl offset SD
	
	gCycleAvgFitData[10] 	:= aEyeVelAmp;						'eye amplitude
	gCycleAvgFitData[11] 	:= sqrt(aEyeVelCovar[0][0]);	'eye amplitude SD
	gCycleAvgFitData[12] 	:= aEyeVelFreq;					'eye frequency
	gCycleAvgFitData[13] 	:= sqrt(aEyeVelCovar[1][1]);	'eye frequency SD
	gCycleAvgFitData[14] 	:= aEyeVelPhase;					'eye phase
	gCycleAvgFitData[15] 	:= sqrt(aEyeVelCovar[2][2]);	'eye phase SD
	gCycleAvgFitData[16] 	:= aEyeVelOffset;					'eye offset
	gCycleAvgFitData[17] 	:= sqrt(aEyeVelCovar[3][3]);	'eye offset SD

	gCycleAvgFitData[18] 	:= aEyeVelAmp/aCtlVelAmp;		'gain
	gCycleAvgFitData[19] 	:= gCycleAvgFitData[19] + gCycleIndFitData[gCyclesAccepted%][19];	'gAveXYnum% * 100 / gMaxAcceptedPts%;	'percent points left after desaccading

	gCycleAvgFitData[20] 	:= 0;									'number of spikes accepted
	gCycleAvgFitData[21] 	:= 0;									'number of spikes accepted SD

	gCycleAvgFitData[22] 	:= aCtlPosAmp;						'Ctl amplitude
	gCycleAvgFitData[23] 	:= sqrt(aCtlPosCovar[0][0]);	'Ctl amplitude SD
	gCycleAvgFitData[24] 	:= aCtlPosFreq;					'Ctl frequency
	gCycleAvgFitData[25] 	:= sqrt(aCtlPosCovar[1][1]);	'Ctl frequency SD
	gCycleAvgFitData[26] 	:= aCtlPosPhase;					'Ctl phase
	gCycleAvgFitData[27] 	:= sqrt(aCtlPosCovar[2][2]);	'Ctl phase SD
	gCycleAvgFitData[28] 	:= aCtlPosOffset;					'Ctl offset
	gCycleAvgFitData[29] 	:= sqrt(aCtlPosCovar[3][3]);	'Ctl offset SD
	
	gCycleAvgFitData[30] 	:= aEyePosAmp;						'eye amplitude
	gCycleAvgFitData[31] 	:= sqrt(aEyePosCovar[0][0]);	'eye amplitude SD
	gCycleAvgFitData[32] 	:= aEyePosFreq;					'eye frequency
	gCycleAvgFitData[33] 	:= sqrt(aEyePosCovar[1][1]);	'eye frequency SD
	gCycleAvgFitData[34] 	:= aEyePosPhase;					'eye phase
	gCycleAvgFitData[35] 	:= sqrt(aEyePosCovar[2][2]);	'eye phase SD
	gCycleAvgFitData[36] 	:= aEyePosOffset;					'eye offset
	gCycleAvgFitData[37] 	:= sqrt(aEyePosCovar[3][3]);	'eye offset SD

	gCycleAvgFitData[46] 	:= aISIAmp;							'isi amplitude
	gCycleAvgFitData[47] 	:= sqrt(aISICovar[0][0]);		'isi amplitude SD
	gCycleAvgFitData[48] 	:= aISIFreq;						'isi frequency
	gCycleAvgFitData[49] 	:= sqrt(aISICovar[1][1]);		'isi frequency SD
	gCycleAvgFitData[50] 	:= aISIPhase;						'isi phase
	gCycleAvgFitData[51] 	:= sqrt(aISICovar[2][2]);		'isi phase SD
	gCycleAvgFitData[52] 	:= aISIOffset;						'isi offset
	gCycleAvgFitData[53] 	:= sqrt(aISICovar[3][3]);		'isi offset SD 

	WriteAveValues();


	if gCyclesAccepted% > 2 then									'can't do STD DEV computation until there are 2 saved trials (n-1)
		ToolbarEnable(gFinalButton, 1);							'now there is data to present in the final screen
	endif
	GetNextCycle();													'get another cycle

	return(kTrue);
end	'Accept



'----------------------------------------------------
'						BinAndFitUnit
' Created so the average unit binning can be called
' after a trial has been accepted.
' 8/31/05
'----------------------------------------------------
proc BinAndFitUnit(aNumTrials)
	var
		aAmp, 
		aFreq, 
		aPhase, 
		aOff,
		aLSfit,
		aCovar[4][4];


	ComputeUnitHistogram(gHistogram[], gAccumulatedAcceptancePulses[], aNumTrials);
	View(gUnitHistoAverageView%);
	XYdelete(1);															'clear the old trace
	XYaddData(kUnitTrial,gxAveArray[],gHistogram[]);			'display the ave control data
	Xrange(0,gEndCycle-gBegCycle);									'limit the x axis
	Optimise();

	'11/16/06  must seed the fit with reasonable values to get the process started.  Control values the best choice.
	aFreq			:= gCtlPosFreq;
	aAmp			:= gCtlPosAmp;
	aPhase		:= gCtlPosPhase;
	aOff			:= gCtlPosOffset;
	ComputeAfit(gUnitHistoAverageView%, kUnitTrial, 0, gEndCycle-gBegCycle, aAmp, aFreq, aPhase, aOff, aLSfit, aCovar[][]);

	'update cycle info for the Log window
	gCycleAvgFitData[38] 	:= aAmp;									'bin amplitude
	gCycleAvgFitData[39] 	:= sqrt(aCovar[0][0]);				'bin amplitude SD
	gCycleAvgFitData[40] 	:= aFreq;								'bin frequency
	gCycleAvgFitData[41] 	:= sqrt(aCovar[1][1]);				'bin frequency SD
	gCycleAvgFitData[42] 	:= aPhase;								'bin phase
	gCycleAvgFitData[43] 	:= sqrt(aCovar[2][2]);				'bin phase SD
	gCycleAvgFitData[44] 	:= aOff;									'bin offset
	gCycleAvgFitData[45] 	:= sqrt(aCovar[3][3]);				'bin offset SD
end	'BinAndFitUnit



'----------------------------------------------------
'						WriteAveValues
' 11/17/06
'----------------------------------------------------
proc WriteAveValues()

	'update cycle info in Log window
	view(gHandleAnalysisAve%);
	EditSelectAll();													'Select all text in log view
	EditClear();														'Delete it
	Print(str$(gCyclesAccepted%) + " cycle avg= " + print$("%5.3f", 1/gDurCycle) + " Hz, " + print$("%5.3f", gDurCycle) + " sec long, " + print$("%5.0f", gCycleAvgFitData[19]/gCyclesAccepted%) + "%% data remaining after desaccading \n");

	Print("Ave Fit Data               Amp       Freq      Phase     Offset \n");
	Print("Ave EyePos  Values: %10.3f %10.3f %10.3f %10.3f \n", gCycleAvgFitData[30], gCycleAvgFitData[32]/k2Pi,gCycleAvgFitData[34]*kRadians,gCycleAvgFitData[36]);
	Print("Ave EyePos  StdDev: %10.3f %10.3f %10.3f %10.3f \n", gCycleAvgFitData[31], gCycleAvgFitData[33], 		gCycleAvgFitData[35], 			gCycleAvgFitData[37]);
	Print("Ave EyeVel  Values: %10.3f %10.3f %10.3f %10.3f \n", gCycleAvgFitData[10], gCycleAvgFitData[12]/k2Pi,gCycleAvgFitData[14]*kRadians,gCycleAvgFitData[16]);
	Print("Ave EyeVel  StdDev: %10.3f %10.3f %10.3f %10.3f \n", gCycleAvgFitData[11], gCycleAvgFitData[13], 		gCycleAvgFitData[15], 			gCycleAvgFitData[17]);
	Print("Ave CtlPos  Values: %10.3f %10.3f %10.3f %10.3f \n", gCycleAvgFitData[22], gCycleAvgFitData[24]/k2Pi,gCycleAvgFitData[26]*kRadians,gCycleAvgFitData[28]);
	Print("Ave CtlPos  StdDev: %10.3f %10.3f %10.3f %10.3f \n", gCycleAvgFitData[23], gCycleAvgFitData[25], 		gCycleAvgFitData[27], 			gCycleAvgFitData[29]);
	Print("Ave CtlVel  Values: %10.3f %10.3f %10.3f %10.3f \n", gCycleAvgFitData[02], gCycleAvgFitData[04]/k2Pi,gCycleAvgFitData[06]*kRadians,gCycleAvgFitData[08]);
	Print("Ave CtlVel  StdDev: %10.3f %10.3f %10.3f %10.3f \n", gCycleAvgFitData[03], gCycleAvgFitData[05], 		gCycleAvgFitData[07], 			gCycleAvgFitData[09]);
	Print("Ave UnitBin Values: %10.3f %10.3f %10.3f %10.3f \n", gCycleAvgFitData[38], gCycleAvgFitData[40]/k2Pi,gCycleAvgFitData[42]*kRadians,gCycleAvgFitData[44]);
	Print("Ave UnitBin StdDev: %10.3f %10.3f %10.3f %10.3f \n", gCycleAvgFitData[39], gCycleAvgFitData[41], 		gCycleAvgFitData[43], 			gCycleAvgFitData[45]);
	Print("Ave UnitISI Values: %10.3f %10.3f %10.3f %10.3f \n", gCycleAvgFitData[46], gCycleAvgFitData[48]/k2Pi,gCycleAvgFitData[50]*kRadians,gCycleAvgFitData[52]);
	Print("Ave UnitISI StdDev: %10.3f %10.3f %10.3f %10.3f"   , gCycleAvgFitData[47], gCycleAvgFitData[49], 		gCycleAvgFitData[51], 			gCycleAvgFitData[53]);

	'standard deviation is defined as the square root of the variance. This means it is the root mean square (RMS) 
	'deviation from the average. It is defined this way in order to give us a measure of dispersion that is (1) a 
	'non-negative number, and (2) has the same units as the data.
end		'WriteAveValues



'----------------------------------------------------
'						ComputeUnitHistogram
' Generate the histogram
' 8/16/05
'----------------------------------------------------
proc ComputeUnitHistogram(aHistogram[], aAccumulatedAcceptancePulses[], aNumTrials)
	var
		indexArray%, k%,
		aSamplesPerBin%, aSecPerBin,
		binAccum,
		aStart,
		aTotalMS, aNumBins,
		aEnd,
		aLength;

	aTotalMS 			:= gLastNumDataPts%;
	aNumBins 			:= gHistoNumBin%;
	aSamplesPerBin%	:= round(aTotalMS / aNumBins);
	aSecPerBin			:= (aTotalMS / aNumBins) / 1000;
	aStart 				:= 0;
	binAccum				:= 0;
	
	aLength := len(aHistogram);
	for indexArray% := 0 to gMaxCycleDur% - 1 do
		
		binAccum += aAccumulatedAcceptancePulses[indexArray%];

		if (indexArray%+1) mod aSamplesPerBin% = 0 then
			aEnd := indexArray%+1;
			if (aStart<aLength) and (aEnd<aLength) then
				for k% := aStart to aEnd do
					aHistogram[k%] := (binAccum / aNumTrials) / aSecPerBin;		'result = spikes/sec   
'					if aNumTrials <> 1 then
'						aNumTrials := aNumTrials;
'					endif
'					if (binAccum > 0)  and  (binAccum < 1) then
'						aNumTrials := aNumTrials;
'					endif
				next
			endif
			aStart 	:= indexArray%+1+1;
			binAccum := 0;
		endif
	next
	'must process a partial bin that may be left over
end	'ComputeUnitHistogram



'----------------------------------------------------
'						AccumulateAve
' 
' 8/16/05
'----------------------------------------------------
proc AccumulateAve(aXYnum%, xArray[], yArray[], aAveMean[], aAveAccum[], aAveContributors[], aZeroTime)
	var
		i%,j;

	'the X array is in milliseconds.  it has to be shifted to 0.0 and 
	'converted to whole numbers so it is compatible with array indexing
	'ArrSub(xArray, xArray[0]);																'shift the data so it starts at 0.0 seconds
	ArrSub(xArray, aZeroTime);																'shift the data so it starts at 0.0 seconds

	'Summing traces like this aligns the beginning of all the accepted traces.  
	'If there is any variability in length the end of the traces
	'may not be perfectly aligned.  Leo says this is fine.  8/4/05
	for i% := 0 to aXYnum%-1 do
		j := round(xArray[i%] / gSampleInterval);										'convert 0.001 values to whole numbers
		aAveAccum[j] += yArray[i%];														'add the current y position to the accumulator
		aAveContributors[j] += 1;															'keep track of how many values there are at this location
		aAveMean[j] := aAveAccum[j] / aAveContributors[j];							'compute the current mean value at the location
	next
end	'AccumulateAve



'----------------------------------------------------
'						SizeAveArrays
' build arrays that ignore the X & Y data from the 
' desaccaded regions
' 8/16/05
'----------------------------------------------------
func SizeAveArrays(aAveContributors[])
	var
		i%,j;

	j := 0;	'initialize to none
	for i% := 0 to gMaxCycleDur% - 1 do
		if aAveContributors[i%] > 0 then
			j += 1;
		endif
	next
	return j;		'num of ave points filled with data
end	'SizeAveArrays



'----------------------------------------------------
'						LoadAveArrays
' build arrays that ignore the X & Y data from the 
' desaccaded regions
' 8/16/05
'----------------------------------------------------
proc LoadAveArrays(xAveArrayDS[], yAveArrayDS[], aAveMean[], aAveContributors[])
	var
		i%,j;

	j := -1;
	for i% := 0 to gMaxCycleDur% - 1 do
		if aAveContributors[i%] > 0 then
			j += 1;
			xAveArrayDS[j] := i% * gSampleInterval;									'array of X values associated with non-desaccaded data
			yAveArrayDS[j] := aAveMean[i%];												'array of Y values from the non-desaccaded data
		endif
	next
end	'LoadAveArrays
 


'----------------------------------------------------
'							ChanDiffRobi
' 5 coefficient FIR derivative (Finite Impulse Response)
'
'< 3/10/04
'----------------------------------------------------
func ChanSmoothRobi(destCh%, sourceCh%, sTime, eTime)
	var 
		dataLength%;
	dataLength% := (eTime-sTime)/BinSize(sourceCh%) + 1;
	
	var 
		dataArray[dataLength%];
	ChanData(sourceCh%, dataArray[], sTime, eTime);							'fill the data array from the source channel

	var 
		coeffFIR[5];
	coeffFIR[0] := 0.4;																'load the filter coefficients
	coeffFIR[1] := 0.1;
	coeffFIR[2] := 0.0;
	coeffFIR[3] := 0.1;
	coeffFIR[4] := 0.4;

	ArrFilt(dataArray[], coeffFIR[]);											'filter the array
	
	ChanWriteWave(destCh%, dataArray[], sTime);								'write the data array to the destination channel
	return dataLength%;
end	'ChanSmoothRobi



'----------------------------------------------------
'							ChanAbsRobi
' rectify a channel and put the results in another channel
'< 3/10/04
'----------------------------------------------------
func ChanAbsRobi(dest%, source%, sTime, eTime)
	var 
		dataLength%;
	
	dataLength% := (eTime-sTime)/BinSize(source%) + 1;
	
	var 
		data[dataLength%];

	ChanData(source%, data[], sTime, eTime);
	Abs(data[]);
	ChanWriteWave(dest%, data[], sTime);
	return dataLength%;
end	'ChanAbsRobi



'----------------------------------------------------
'							ChanSubRobi
' subtract two arrays
'< 3/10/04
'----------------------------------------------------
func ChanSubRobi(destCh%, source1Ch%, source2Ch%, sTime, eTime)
	var 
		dataLength%;

	dataLength% := (eTime-sTime)/BinSize(source1Ch%) + 1;
	
	var
		datadest[dataLength%],
		datasource[dataLength%];

	ChanData(source1Ch%, datadest[], sTime, eTime);
	ChanData(source2Ch%, datasource[], sTime, eTime);
	ArrSub(datadest[], datasource[]);
	ChanWriteWave(destCh%, datadest[], sTime);
	return dataLength%;
end



'----------------------------------------------------
'							ChanDiffRobi
' 5 coefficient FIR derivative (Finite Impulse Response)
'
'< 3/10/04
'----------------------------------------------------
func ChanDiffRobi(destCh%, sourceCh%, sTime, eTime)
	var 
		dataLength%;

	dataLength% := (eTime-sTime)/BinSize(sourceCh%) + 1;
	
	var 
		dataArray[dataLength%];

	ChanData(sourceCh%, dataArray[], sTime, eTime);		'fill the data array from the source channel

	var 
		coeffFIR[5];

	coeffFIR[0] := -0.2;											'load the filter coefficients
	coeffFIR[1] := -0.1;
	coeffFIR[2] :=  0.0;
	coeffFIR[3] :=  0.1;
	coeffFIR[4] :=  0.2;

	ArrFilt(dataArray[], coeffFIR[]);						'filter the array
	ArrDiv(dataArray[], BinSize(sourceCh%));
	
	
	ChanWriteWave(destCh%, dataArray[], sTime);			'write the data array to the destination channel

	return dataLength%;
end	'ChanDiffRobi