function eventHandler(win,x,y,ibut)
	global("glb");
	if glb.locked; return; end;
	glb.locked = %t; 
	if ibut==-1000; return; end;
	sca(glb.axs.da);
	[x,y] = xchange(x,y,'i2f');
	
	select ibut;
		case 0; //Left pressed, set zoom start
			glb.zoomList = [glb.zoomList; x, %nan];
		case -5; // Left released, set zoom end
			glb.zoomList($,2) = x;
			doZoom();
		case 5; //Right clicked - zoom out
			if size(glb.zoomList, 1) > 1;
				glb.zoomList = glb.zoomList(1:($-1), :);
				doZoom();
			end;
		case 2; //Right press, set DA threshold start
			glb.lThreshold = y;
		case -3; //Right release, set DA threshold end
			if y < glb.lThreshold;
				glb.uThreshold = glb.lThreshold
				glb.lThreshold = y;
			else
				glb.uThreshold = y;
			end;
			calcDAlphaMarkers();
			plotDAlphaMarkers();
		case 10; //Left double-click			
			glb.win(1) = x;
			calcDAlphaMarkers();
			plotDAlphaMarkers();
		case 12; //Right double-click
			glb.win(2) = x;
			calcDAlphaMarkers();
			plotDAlphaMarkers();
		case 52; // 3 AD34
			switchDAlphaChannel(34);
		case 53; // 3 AD35
			switchDAlphaChannel(35);
		case 54; // 3 AD36
			switchDAlphaChannel(36);
//		case 114; // 'r': set slice radius at mouse pos
// ------ doesn't work ------
//			if glb.axs.neAllAxIsSet;
//				sca(glb.axs.neAll);
//				[x,y] = xchange(x,y,'i2f');
//				xinfo(msprintf("Future HRTS slice R = %f", x));
//				glb.sliceR = x;
//			end;
// --------------------------
		case -1; //mouse move, ignore
		else //anything else
			printf("%i ", ibut);

	end;

	glb.locked = %f;
endfunction

function plotDAlphaMarkers()
	global("glb");
	drawlater;
	ax = glb.axs.da;
	sca(ax);
	if size(ax.children, 1) > 1;
		delete(ax.children(1:($-1)));
	end;

	plot([glb.win(1), glb.win(1)], [min(glb.da), max(glb.da)], "r");
	plot([glb.win(2), glb.win(2)], [min(glb.da), max(glb.da)], "r");
	plot([min(glb.daTime), max(glb.daTime)], [glb.lThreshold, glb.lThreshold], "--k");
	plot([min(glb.daTime), max(glb.daTime)], [glb.uThreshold, glb.uThreshold], "--k");

	if length(glb.elmTimes) > 0;
		plot([glb.elmTimes, glb.elmTimes], [min(glb.da), max(glb.da)], "g");
	end;

	doZoom();

endfunction

// Sets the zoom to the current end of the zoom list
function doZoom()
	global("glb");
	drawlater;
	subplot(3,1,1);

	if length(glb.da) > 0;
		ax = glb.axs.da;
		i = glb.daTime >= glb.zoomList($,1) & glb.daTime < glb.zoomList($,2);
		minY = min(glb.da(i)); maxY = max(glb.da(i));
		ax.data_bounds = [glb.zoomList($, 1), minY; glb.zoomList($,2), maxY];
		ax.tight_limits = "on";
	end;

	if length(glb.wdia) > 0;
		ax = glb.axs.wdia;
		i = glb.wdiaTime >= glb.zoomList($,1) & glb.wdiaTime < glb.zoomList($,2);
		minY = min(glb.wdia(i)); maxY = max(glb.wdia(i));
		ax.data_bounds = [glb.zoomList($, 1), minY; glb.zoomList($,2), maxY];
		ax.tight_limits = "on";
	end;

	if length(glb.lid4) > 0;
		ax = glb.axs.lid4;
		i = glb.lid4Time >= glb.zoomList($,1) & glb.lid4Time < glb.zoomList($,2);
		minY = min(glb.lid4(i)); maxY = max(glb.lid4(i));
		ax.data_bounds = [glb.zoomList($, 1), minY; glb.zoomList($,2), maxY];
		ax.tight_limits = "on";
	end;

	if type(glb.axs.markers) == 9;
		ax = glb.axs.markers;
		ax.data_bounds = [glb.zoomList($, 1), 0; glb.zoomList($,2), 1];
		ax.tight_limits = "on";
	end;

	drawnow;
endfunction

function calcDAlphaMarkers()
	global("glb");

	glb.elmTimes = [];

	if glb.lThreshold == 0 | glb.uThreshold == 0; return; end;

	inRange = glb.daTime(1:($-1)) >= glb.win(1) & glb.daTime(2:$) < glb.win(2);

	//find all points where thresholds are crossed
	lCrossedUp = find(glb.da(1:($-1)) < glb.lThreshold & glb.da(2:$) >= glb.lThreshold & inRange);
	uCrossedUp = find(glb.da(1:($-1)) < glb.uThreshold & glb.da(2:$) >= glb.uThreshold & inRange);
	lCrossedDown = find(glb.da(1:($-1)) >= glb.lThreshold & glb.da(2:$) < glb.lThreshold & inRange);
	uCrossedDown = find(glb.da(1:($-1)) >= glb.uThreshold & glb.da(2:$) < glb.uThreshold & inRange);

	if length(lCrossedUp) <= 0 | length(uCrossedUp) <= 0; return; end;

	//put all the crossings together and sort, make upper crossings very slightly higher, so order it kept
	crossIdx = [lCrossedUp, uCrossedUp + 1e-5, lCrossedDown + 2e-5, uCrossedDown + 3e-5];
	crossType = [ones(lCrossedUp)*1, ones(uCrossedUp)*2, ones(lCrossedDown)*3, ones(uCrossedDown)*4];
	[a,i] = sort(crossIdx);
	i=i($:-1:1);
	crossIdx = crossIdx(i);
	crossType = crossType(i);

	//look for pairs of crossing both lower and upper upwards, in a row
	crossBoth = find(crossType(1:($-1)) == 1 & crossType(2:$) == 2);

	if length(crossBoth) <= 0; return; end;

	//turn back into an index into DA (well, DA differences anyway)
	crossBothAvIdx = int((crossIdx(crossBoth) + crossIdx(crossBoth+1))/2);

	glb.elmTimes = (glb.daTime(crossBothAvIdx) + glb.daTime(crossBothAvIdx+1)) / 2;

	xinfo(msprintf("%i ELMs marked in time window %6.3f-%6.3f", length(glb.elmTimes), glb.win(1), glb.win(2))); 

endfunction

function selectAD34(a, b); if glb.locked; return; end; glb.locked = %t; switchDAlphaChannel(34); glb.locked = %f; endfunction
function selectAD35(a, b); if glb.locked; return; end; glb.locked = %t; switchDAlphaChannel(35); glb.locked = %f; endfunction
function selectAD36(a, b); if glb.locked; return; end; glb.locked = %t; switchDAlphaChannel(36); glb.locked = %f; endfunction

function loadOthers(ch)
	global("glb");

	if length(glb.wdia) <= 0;
		[glb.wdia, glb.wdiaTime, jnk, glb.wdiaStat] = getSig("jet", glb.pulse, "ppf-chain1/btli/wdia");
		if length(glb.wdia) > 0; 
			ax=newaxes();
			plot(glb.wdiaTime, glb.wdia, "k");
			ax.axes_bounds = glb.axs.da.axes_bounds;
			ax.background = -99;
			ax.y_location = "right";
			glb.axs.wdia = ax;

		else
			printf("No data in ''"+string(pulse)+"/ppf-chain1/btli/wdia''.\n");
		end;
	end;

	if length(glb.lid4) <= 0;
		[glb.lid4, glb.lid4Time, jnk, glb.lid4Stat] = getSig("jet", glb.pulse, "ppf/kg1v/lid4");
		if length(glb.lid4) > 0; 
			ax=newaxes();
			plot(glb.lid4Time, glb.lid4, "r");
			ax.axes_bounds = glb.axs.da.axes_bounds;
			ax.background = -99;
			ax.y_location = "right";
			glb.axs.lid4 = ax;

		else
			printf("No data in ''"+string(pulse)+"/ppf-/kg1v/lid4''.\n");
		end;
	end;
	
endfunction

function switchDAlphaChannel(ch)
	global("glb");
	drawlater;
	[glb.da glb.daTime] = getSig("jet", pulse, "ppf/s3ad/ad" + string(ch));


	clf;
	if type(glb.axs.da) == 9;
		ax = glb.axs.da;
//		delete(ax.children);
	else
		subplot(3,1,1);
		ax=gca();
		glb.axs.da = ax;		
	end;
	ax.background = -99;

	plot(glb.daTime, glb.da, "b");

	glb.zoomList = [min(glb.daTime), max(glb.daTime)];
	if length(glb.win) < 2; glb.win = glb.zoomList; end;

	loadOthers();
	doZoom();

endfunction


function sync=elmSync(times)
	global("glb");

	n = length(times);

	sync.prevElmIdx = -ones(n, 1);
	sync.nextElmIdx = -ones(n, 1);
	sync.tFromElm = -ones(n, 1);
	sync.tToElm = -ones(n, 1);
	sync.elmFrac = -ones(n, 1);

	if length(glb.elmTimes) <= 0; return; end;

	for i=1:n;
		t = times(i);

		j = find(glb.elmTimes < t);
		if length(j) <= 0; continue; end;
		sync.prevElmIdx(i) = j($);
		sync.tFromElm(i) = t - glb.elmTimes(sync.prevElmIdx(i));

		j = find(glb.elmTimes >= t);
		if length(j) <= 0; continue; end;
		sync.nextElmIdx(i) = j(1);
		sync.tToElm(i) = glb.elmTimes(sync.nextElmIdx(i)) - t;

		sync.elmFrac(i) = sync.tFromElm(i) / (sync.tFromElm(i) + sync.tToElm(i));
		
		if modulo(i, int(n/100)) == 0; printf(" %i%%", i*100.0/n); end;
	end;

endfunction

function clearSubAxes();
	global("glb");
	//get rid of HRTS sync plots, if present
	if type(glb.axs.neAll) == 9; delete(glb.axs.neAll); glb.axs.neAll = []; end;
	if type(glb.axs.neAvg) == 9; delete(glb.axs.neAvg); glb.axs.neAvg = []; end;
	if type(glb.axs.TeAll) == 9; delete(glb.axs.TeAll); glb.axs.TeAll = []; end;
	if type(glb.axs.TeAvg) == 9; delete(glb.axs.TeAvg); glb.axs.TeAvg = []; end;
	if type(glb.axs.neSliceSync) == 9; delete(glb.axs.neSliceSync); glb.axs.neSliceSync = []; end;
	if type(glb.axs.TeSliceSync) == 9; delete(glb.axs.TeSliceSync); glb.axs.TeSliceSync = []; end;

	if type(glb.axs.wdiaSync) ==9; delete(glb.axs.wdiaSync); glb.axs.wdiaSync = []; end;
	if type(glb.axs.lid4Sync) ==9; delete(glb.axs.lid4Sync); glb.axs.lid4Sync = []; end;

	if type(glb.axs.markers) == 9; delete(glb.axs.markers); glb.axs.markers = []; end;
endfunction

function doLWSort(a,b);
	global("glb");

	if a <> 99 & glb.locked; return; end;
	glb.locked = %t;

	drawlater;

	clearSubAxes();

	if length(glb.wdia) > 0;
		printf("WDIA: ");
		glb.wdiaSync = elmSync(glb.wdiaTime);
		printf("\n");

		subplot(2,2,3); ax=gca();
		i = glb.wdiaSync.tFromElm > 0;
		plot2d(glb.wdiaSync.tFromElm(i), glb.wdia(i), -1);
		ax.data_bounds(:,1) = [0; 0.1];
		glb.axs.wdiaSync = ax;
	end;

	if length(glb.lid4) > 0;
		printf("LID4: ");
		glb.lid4Sync = elmSync(glb.lid4Time);
		printf("\n");

		subplot(2,2,4); ax=gca();
		i = glb.lid4Sync.tFromElm > 0;
		plot2d(glb.lid4Sync.tFromElm(i), glb.lid4(i), -1);
		ax.data_bounds(:,1) = [0; 0.1];
		glb.axs.lid4Sync = ax;
	end;

	plotDAlphaMarkers();
	doZoom();

	glb.locked = %f;
endfunction

function doHRTSSort(a,b)
	global("glb");
	drawlater;

	if a <> 99 & glb.locked; return; end;
	glb.locked = %t;

	if length(glb.elmTimes) <= 0; xinfo("No ELMs marked"); drawnow; return; end;

	glb.hrtsSync = elmSync(glb.hrtsTimes);

	nHRTS = length(glb.hrtsTimes);

	clearSubAxes();
	
	glb.axs.markers = newaxes();
	glb.axs.markers.background = -99;
	glb.axs.markers.axes_bounds = glb.axs.da.axes_bounds;

	subplot(3,3,4); glb.axs.neAll=gca();
	subplot(3,3,5);	glb.axs.neAvg = gca();
	subplot(3,3,7);	glb.axs.TeAll = gca();
	subplot(3,3,8);	glb.axs.TeAvg = gca();

	nPlots = 6;
	df = 1 / nPlots;
	col = [color(128,0,128), color(255,0,0), color(128,128,0), color(0,255,0), color(0,128,128), color(0,0,255)];
	for i=nPlots:-1:1;
		f0 = (i-1)*df;
		f1 = i*df;

		j = find(glb.hrtsSync.elmFrac >= f0 & glb.hrtsSync.elmFrac < f1);

		sca(glb.axs.markers);
		plot2d([glb.hrtsTimes(j), glb.hrtsTimes(j)]', [zeros(j); ones(j)], col(i)+0*glb.hrtsTimes(j));
		glb.axs.markers.children(1).children(:).line_style = 3;
		glb.axs.markers.children(1).children(:).thickness = 2;

		sca(glb.axs.neAll);
		plot2d(glb.hrtsR' * ones(1, length(j)), glb.hrtsNe(j, :)', col(i) * ones(1, length(j)));

		sca(glb.axs.neAvg);
		plot2d(glb.hrtsR, mean(glb.hrtsNe(j, :), "r"), col(i));

		sca(glb.axs.TeAll);
		plot2d(glb.hrtsR' * ones(1, length(j)), glb.hrtsTe(j, :)', col(i) * ones(1, length(j)));

		sca(glb.axs.TeAvg);
		plot2d(glb.hrtsR, mean(glb.hrtsTe(j, :), "r"), col(i));

	end;
	glb.axs.neAll.data_bounds = [3.0, 0; 4.0, 1.5e20];
	glb.axs.neAvg.data_bounds = [3.0, 0; 4.0, 1.5e20];
	glb.axs.TeAll.data_bounds = [3.0, 0; 4.0, 8000];
	glb.axs.TeAvg.data_bounds = [3.0, 0; 4.0, 8000];
	//glb.axs.neAll.data_bounds = [3.7, 0; 4.0, 1.5e20];
	//glb.axs.neAvg.data_bounds = [3.7, 0; 4.0, 1.5e20];
	//glb.axs.TeAll.data_bounds = [3.7, 0; 4.0, 2000];
	//glb.axs.TeAvg.data_bounds = [3.7, 0; 4.0, 2000];

	i = glb.hrtsSync.elmFrac > 0;

	//index of slice
	sliceIdx0 = find(glb.hrtsR(1:($-1)) < glb.sliceR & glb.hrtsR(2:$) > glb.sliceR);
	sliceFrac = (glb.sliceR - glb.hrtsR(sliceIdx0)) / (glb.hrtsR(sliceIdx0+1) - glb.hrtsR(sliceIdx0));
	sliceNe = (1-sliceFrac)*glb.hrtsNe(i, sliceIdx0) + sliceFrac*glb.hrtsNe(i, sliceIdx0+1);
	sliceTe = (1-sliceFrac)*glb.hrtsTe(i, sliceIdx0) + sliceFrac*glb.hrtsTe(i, sliceIdx0+1);

	subplot(3,3,6);
	ax=gca();
	if size(ax.children, 1) > 0; delete(ax.children); end;
	plot2d(glb.hrtsSync.elmFrac(i), sliceNe, -1);
//	plot2d(glb.hrtsSync.tFromElm(i), glb.hrtsNe(i, 53) - glb.hrtsNe(i, 45), -1);
	//ax=gca();
	//ax.data_bounds = [0, 0; 1, 1.5e20];
	glb.axs.neSliceSync = ax;

	subplot(3,3,9);
	ax=gca();
	if size(ax.children, 1) > 0; delete(ax.children); end;
	plot2d(glb.hrtsSync.elmFrac(i), sliceTe, -1);
//	plot2d(glb.hrtsSync.tFromElm(i), glb.hrtsTe(i, 53) - glb.hrtsTe(i, 45), -1);
	//ax=gca();
	//ax.data_bounds = [0, 0; 1, 1000];
	glb.axs.TeSliceSync = ax;

	doZoom();

	glb.locked = %f;
endfunction

function doKE9Plots(a,b)
	global("glb");

	if a <> 99 & glb.locked; return; end;
	glb.locked = %t;

	drawlater;

	if glb.ke9UID <> ""; uidStr = "-" + glb.ke9UID; else uidStr = ""; end;

	[glb.ke9Ne, glb.ke9Times, glb.ke9R, glb.ke9NeStat] = getSig("jet", pulse, "ppf"+uidStr+"/ke9d/ne");
	[glb.ke9Te, glb.ke9Times, glb.ke9R, glb.ke9TeStat] = getSig("jet", pulse, "ppf"+uidStr+"/ke9d/te");
	[glb.ke9Z, jnk, glb.ke9R] = getSig("jet", pulse, "ppf"+uidStr+"/ke9d/zvec");

	// for every ke9 time point, we need both the ke9 and hrts psiN for remapping
	printf("Remapping KE9...");
	glb.ke9PsiN = zeros(length(glb.ke9Times), length(glb.ke9R));
	glb.ke9R_HRTS = glb.ke9PsiN;
	for i=1:length(glb.ke9Times);
		ret = flushReq(glb.pulse, glb.ke9Times(i), [glb.ke9R, glb.hrtsR], [glb.ke9Z, glb.hrtsZ], "EFIT", "");
		glb.ke9PsiN(i, :) = ret(1:length(glb.ke9R),3)';

		hrtsPsiN = ret((length(glb.ke9R)+1):$,3)';


		glb.ke9R_HRTS(i, :) = interpln([hrtsPsiN; glb.hrtsR], glb.ke9PsiN(i, :));
		printf(".");
	end;
	printf("Done.");

	glb.ke9Sync = elmSync(glb.ke9Times);

	if type(glb.axs.markers) <> 9;
		glb.axs.markers = newaxes();
		glb.axs.markers.background = -99;
		glb.axs.markers.axes_bounds = glb.axs.da.axes_bounds;
	end;
	if type(glb.axs.neAll) <> 9; subplot(3,3,4); glb.axs.neAll=gca(); end;
	if type(glb.axs.neAvg) <> 9; subplot(3,3,5); glb.axs.neAvg = gca(); end;
	if type(glb.axs.TeAll) <> 9; subplot(3,3,7); glb.axs.TeAll = gca(); end;
	if type(glb.axs.TeAvg) <> 9; subplot(3,3,8); glb.axs.TeAvg = gca(); end;

	nPlots = 6;
	df = 1 / nPlots;
	col = [color(128,0,128), color(255,0,0), color(128,128,0), color(0,255,0), color(0,128,128), color(0,0,255)];
	for i=nPlots:-1:1;
		f0 = (i-1)*df;
		f1 = i*df;

		j = find(glb.ke9Sync.elmFrac >= f0 & glb.ke9Sync.elmFrac < f1);
		if length(j) <= 0; continue; end
printf("A");
		sca(glb.axs.markers);
		plot2d([glb.ke9Times(j), glb.ke9Times(j)]', [zeros(j); ones(j)], col(i)+0*glb.ke9Times(j));
		glb.axs.markers.children(1).children(:).line_style = 3;
		glb.axs.markers.children(1).children(:).thickness = 2;
printf("B");
		sca(glb.axs.neAvg);
		plot2d(glb.ke9R_HRTS(j, :)', glb.ke9Ne(j, :)', col(i) * ones(1, length(j)));
printf("C");
		//sca(glb.axs.neAvg);
		//plot2d(glb.ke9R_HRTS(j, :), mean(glb.ke9Ne(j, :), "r"), col(i));
printf("D");
		sca(glb.axs.TeAvg);
		plot2d(glb.ke9R_HRTS(j, :)', glb.ke9Te(j, :)', col(i) * ones(1, length(j)));
printf("E");
		//sca(glb.axs.TeAvg);
		//plot2d(glb.ke9R_HRTS(j, :), mean(glb.ke9Te(j, :), "r"), col(i));
printf("F");
	end;
	
	doZoom();

	glb.locked = %f;
endfunction




printf("pulse = %i\n", pulse);



if ~isdef("glb");
	glb.pulse = pulse;
	glb.da = [];
	glb.daTime = [];
	glb.zoomList = [];
	glb.win = [];
	glb.lThreshold = 0;
	glb.uThreshold = 0;
	glb.elmTimes = [];
	glb.sliceR = 3.81;
	glb.lid4 = [];
	glb.lid4Time = [];
	glb.wdia = [];
	glb.wdiaTime = [];
	glb.hrtsPsiN = [];

	if isdef("ke9UID"); glb.ke9UID = ke9UID; else glb.ke9UID = ""; end;
end;

global("glb");
clf;
glb.da = [];
glb.lid4 = [];
glb.wdia = [];

//clear all axes
glb.axs.da = 0;
glb.axs.wdia = 0;
glb.axs.lid4 = 0;
glb.axs.neAll = 0;
glb.axs.TeAll = 0;
glb.axs.neAvg = 0;
glb.axs.TeAvg = 0;
glb.axs.neSliceSync = 0;
glb.axs.TeSliceSync = 0;
glb.axs.wdiaSync = 0;
glb.axs.lid4Sync = 0;
glb.axs.markers = 0;

glb.locked = %f;

[glb.hrtsNe, glb.hrtsTimes, glb.hrtsR, glb.hrtsNeStat] = getSig("jet", pulse, "ppf/hrts/ne");
[glb.hrtsTe, glb.hrtsTimes, glb.hrtsR, glb.hrtsTeStat] = getSig("jet", pulse, "ppf/hrts/te");
[glb.hrtsZ, jnk, glb.hrtsR] = getSig("jet", pulse, "ppf/hrts/z");
[glb.ke9Ne, glb.ke9Times, glb.ke9R, glb.ke9NeStat] = getSig("jet", pulse, "ppf/ke9/ne");
[glb.ke9Te, glb.ke9Times, glb.ke9R, glb.ke9TeStat] = getSig("jet", pulse, "ppf/ke9/te");
//glb.ke9RMag

// too slow
//printf("Remapping HRTS...");
//glb.hrtsPsiN = zeros(length(glb.hrtsTimes), length(glb.hrtsR));
//for i=1:length(glb.hrtsTimes);
//	ret = flushReq(glb.pulse, glb.hrtsTimes(i), glb.hrtsR, glb.hrtsZ, "EFIT", "");
//	glb.hrtsPsiN(i, :) = ret(:,3)';
//	
//	if modulo(i, 10) == 0; printf("(%i/%i)", i, length(glb.hrtsTimes)); end;
//	printf(".");
//end;
//printf("Done.");


clf;

addmenu(0, "4) DAlpha=AD34", list(2,"selectAD34"));
addmenu(0, "5) DAlpha=AD34", list(2,"selectAD35"));
addmenu(0, "6) DAlpha=AD34", list(2,"selectAD36"));
addmenu(0, "H) Sort HRTS", list(2,"doHRTSSort"));
addmenu(0, "L) Sort LID4/WDIA", list(2,"doLWSort"));
addmenu(0, "9) Plot KE9", list(2,"doKE9Plots"));


seteventhandler('eventHandler')


switchDAlphaChannel(35);
loadOthers();
calcDAlphaMarkers();
plotDAlphaMarkers();

