% JaeSub Hong, 2003-2005, version 1.5
% Please report any problem or suggestion at jaesub@head.cfa.harvard.edu

define imcounts (x,y) {

	variable i,n_data;
	variable ans;

	n_data = length(y);
	ans = Double_Type [n_data];

	ans[0] = (x.hi[0]-x.lo[0])*y[0];
	for (i=1;i<n_data;i++){
		ans[i] = (x.hi[i]-x.lo[i])*y[i]+ans[i-1];
	}

	for (i=0;i<n_data;i++){
		ans[i] = ans[i]/ans[n_data-1];
	}

	return ans;
}

define imcounts_raw (x,y) {

	variable i,n_data;
	variable ans;

	n_data = length(y);
	ans = Double_Type [n_data];

	ans[0] = 1.0*y[0];
	for (i=1;i<n_data;i++){
		ans[i] = 1.0*y[i]+ans[i-1];
	}

	for (i=0;i<n_data;i++){
		ans[i] = ans[i]/ans[n_data-1];
	}

	return ans;
}

define value_locate (x, nx, n_data) {

	variable middle;
	variable ilo, ihi;

	ihi = n_data;
	ilo = 0;

	if (nx >=x[ihi]) return ihi;
	if (nx <=x[ilo]) return ilo;

	forever {
	    	middle = typecast ((ilo + ihi)/2., Integer_Type);
	    	if (middle == ilo)  return ilo; 
	    	if (nx < x[middle]) ihi = middle; else ilo = middle; 
	}

}

define interpol (y, x, nx) {

	variable last, n_req;
	variable i,j,k;
	variable dy, dx;
	variable ans;

	last = length(x)-1;
	n_req = length(nx);
	ans = Double_Type [n_req];

	for (i=0;i<n_req;i++){
		j = value_locate(x, nx[i], last);
		if (j >= last) j = last-1;
		k = j + 1;
		dx = 1.0*(x[k]-x[j]);
		dy = (y[k]-y[j])/dx;
		ans[i]  = dy*(nx[i]-x[j]) + y[j];
	}

	return ans;

}

define mod_quantile(data, frac) {

	variable x,y,z,ans;

	x=get_energy_axes(data);
	y=get_mcounts(data);
	z=imcounts(x,y);
	ans = interpol(x.hi,z,frac);
	return ans;

}

define mod_quantile_raw(data, frac) {

	variable x,y,z,ans;

	x=get_raw_axes(data);
	y=get_data(data);
	z=imcounts_raw(x,y);
	ans = interpol(x.mid,z,frac);
	return ans;

}

define pri_quantile(quan) {

	variable i,n_frac,ans;

	n_frac = length(quan);
	ans = string(quan[0]);

	for (i=1;i<n_frac;i++){
		ans=ans+"\t"+string(quan[i]);
	}
%	print(ans);

	return ans;

}