int find_peak(int n, float bw, float approx_wl, float **data)

{
	int i, peak_index;
	float y_max;
	
	y_max = 0;
	for(i = 1; i <= n; i++){
		if(data[i][1] > (approx_wl - bw)){
			/*printf("i = %d wl = %f sig = %f\n", i, data[i][1], data[i][2]);*/
			if(data[i][2] > y_max){ 
				y_max = data[i][2];
				peak_index = i;
			}
		}
		if(data[i][1] > (approx_wl + bw)) break;
	}
	
	if(y_max == 0) peak_index = 0;
	if(i == 1) peak_index = 0;
	
	return peak_index;
}

/* ----------------------------------------------------------------------- */

int find_limits_bw(int peak_index, int n, float bw, float **data, int *i_min, int *i_max)

{
	int i;
	
	i = peak_index;
	while(data[i][1] >= (data[peak_index][1] - bw)){
		*i_min = i;
		i--;
		if(i < 1) break;
	}
	
	i = peak_index;
	while(data[i][1] <= (data[peak_index][1] + bw)){
		*i_max = i;
		i++;
		if(i > n) break;
	}
	
	return 1;
}

/* --------------------------------------------------------------- */

int find_limits_fr(int peak_index, int n, float fraction, float **data, 
		int *i_min, int *i_max)
{
	int i;
	
	i = peak_index;
	while(data[i][2] >= (data[peak_index][2]*fraction)){
		*i_min = i;
		i--;
		if(i < 1) break;
	}
	
	i = peak_index;
	while(data[i][2] >= (data[peak_index][2]*fraction)){
		*i_max = i;
		i++;
		if(i > n) break;
	}
	
	return 1;
}

/* ----------------------------------------------------------------------- */

int bksub_calc(int n_wl, float bw, float line, float mult,
				float **data, float **bk_data)
				
{
	int i, j, k, u, i1_min, i1_max, i2_min, i2_max, i_pk;
	float y_min, y_max;
	
	i_pk = find_peak(n_wl, bw, line, data);
	if(i_pk != 0){
		find_limits_bw(i_pk, n_wl, bw, data, &i1_min, &i1_max);
		find_limits_bw(i_pk, n_wl, bw*mult, data, &i2_min, &i2_max);
		y_min = 1e10;
		for(i = i2_min; i <= i1_min; i++){
			if(i > 0){
				if(data[i][2] < y_min) y_min = data[i][2];
			}
		}
		for(i = i1_max; i <= i2_max; i++){
			if(i <= n_wl){
				if(data[i][2] < y_min) y_min = data[i][2];
			}
		}
		for(i = i2_min; i <= i2_max; i++){	
			bk_data[i][1] = data[i][1];
			bk_data[i][2] = data[i][2] - y_min;
			bk_data[i][3] = data[i][3];
		}
	}
	/*printf("Done with background subtraction\n");*/
	
	return 1;
}

/* ----------------------------------------------------------------------- */

int peak_intercept_calc(int i_pk, float bw, float **data, float *peak)

{
	int j, k, n1, n2;
	float *x1, *x2, *y1, *y2, *dy;
	float inter1, inter2, slope1, slope2, dum;
	
	x1 = vector(1, 20);
	x2 = vector(1, 20);
	y1 = vector(1, 20);
	y2 = vector(1, 20);
	
	j = i_pk - 1;
	k = 1;
	while(data[j][1] > (data[i_pk][1] - bw)){
		x1[k] = data[j][1];
		y1[k] = data[j][2];
		j--;
		k++;
	}
	n1 = k - 1;
	
	j = i_pk + 1;
	k = 1;
	while(data[j][1] < (data[i_pk][1] + bw)){
		x2[k] = data[j][1];
		y2[k] = data[j][2];
		j++;
		k++;
	}
	n2 = k - 1;
	
	fit(x1, y1, n1, dy, 0, &inter1, &slope1, &dum, &dum, &dum, &dum);
	fit(x2, y2, n2, dy, 0, &inter2, &slope2, &dum, &dum, &dum, &dum);
	*peak = (inter2 - inter1)/(slope1 - slope2);
	
	free_vector(x1, 1, 20);
	free_vector(x2, 1, 20);
	free_vector(y1, 1, 20);
	free_vector(y2, 1, 20);
	
	return 1;
}
/* ----------------------------------------------------------------------- */

int peak_fit_calc(int i_pk, float **data, float *peak)

{
	int j, k, n1, n2;
	float *x1, *y1, *y2;
	float min, wl_fit, fit_drv;
	
	x1 = vector(1, 20);
	y1 = vector(1, 20);
	y2 = vector(1, 20);
	
	k = 1;
	for(j = i_pk - 2; j <= i_pk + 2; j++){
		x1[k] = data[j][1];
		y1[k] = data[j][2];
		k++;
	}
	n1 = k - 1;
	min = 1e10;
	spline(x1, y1, n1, 1e40, 1e40, y2);
	for(wl_fit = x1[2]; wl_fit <= x1[4]; wl_fit += 0.001){
		spldrv(x1, y1, y2, n1, wl_fit, &fit_drv);
		if(fabs(fit_drv) < min){ 
			min = fabs(fit_drv);
			*peak = wl_fit;
		}
	}
	
	free_vector(x1, 1, 20);
	free_vector(y1, 1, 20);
	free_vector(y2, 1, 20);
	
	return 1;
}

/* ----------------------------------------------------------------------- */

int centroid_calc(int min_index, int max_index, float **data,
				float *centroid, float *centroid_unc)

{
	int i;
	double sum_y, sum_xy;
	double u;
	
	sum_y = 0, sum_xy = 0;
	for(i = min_index; i <= max_index; i++){	
		sum_y += data[i][2];
		sum_xy += data[i][1]*data[i][2];
	}
	*centroid = sum_xy/sum_y;
	
	*centroid_unc = 0;
	for(i = min_index; i <= max_index; i++){
		*centroid_unc += SQR((data[i][2]*sum_y - sum_xy)/SQR(sum_y)*data[i][3]);
	}
	*centroid_unc = sqrt(*centroid_unc);
	
	return 1;
}

/* ----------------------------------------------------------------- */

int fwhm_calc(int peak_index, float **data,
				float *min_halfmax, float *max_halfmax,		
				float *min_halfmax_unc, float *max_halfmax_unc)
				
{
	int i;
	int i_ss, i_sl, i_ls, i_ll;
	float y_max;
	
	y_max = data[peak_index][2];
	
	i_ss = 0, i = peak_index;
	while(i_ss == 0){
		if(data[i][2] < y_max/2.0){
			i_ss = i;
			i_sl = i+1;
		}
		else{
			i = i-1;
		}
	}
	
	i_ll = 0, i = peak_index;
	while(i_ll == 0){
		if(data[i][2] < y_max/2.0){
			i_ll = i;
			i_ls = i-1;
		}
		else{
			i = i+1;
		}
	}
	
	*min_halfmax = data[i_sl][1] + (data[i_ss][1]-data[i_sl][1])*
					((data[i_sl][2]-y_max/2.0)/(data[i_sl][2]-data[i_ss][2]));
	*max_halfmax = data[i_ll][1] + (data[i_ls][1]-data[i_ll][1])*
					((data[i_ll][2]-y_max/2.0)/(data[i_ll][2]-data[i_ls][2]));
					
	*min_halfmax_unc = 0;
	*max_halfmax_unc = 0;
	/*
	*min_halfmax_unc += SQR((x[i_ss]-x[i_sl])/(y[i_sl]-y[i_ss])/2.0*dy[peak_index]);
	*min_halfmax_unc += SQR((x[i_ss]-x[i_sl])*(y[i_sl]-y_max/2.0)/
							SQR(y[i_sl]-y[i_ss])*dy[i_ss]);
	*min_halfmax_unc += SQR((x[i_ss]-x[i_sl])*(y_max/2.0-y[i_ss])/
							SQR(y[i_sl]-y[i_ss])*dy[i_sl]);
	*min_halfmax_unc = sqrt(*min_halfmax_unc);
	
	*max_halfmax_unc += SQR((x[i_ls]-x[i_ll])/(y[i_ll]-y[i_ls])/2.0*dy[peak_index]);
	*max_halfmax_unc += SQR((x[i_ls]-x[i_ll])*(y[i_ll]-y_max/2.0)/
							SQR(y[i_ll]-y[i_ls])*dy[i_ls]);
	*max_halfmax_unc += SQR((x[i_ls]-x[i_ll])*(y_max/2.0-y[i_ls])/
							SQR(y[i_ll]-y[i_ls])*dy[i_ll]);
	*max_halfmax_unc = sqrt(*max_halfmax_unc);
	*/
	return 1;
}


/* ---------------------------------------------------------- */

int calc_filter_cents(FILE *fp_out, int i_pk, int n_wl, float bw, 
		float **trans)
				
{
	int i;
	int i_min, i_max;
	float *trans_unc;
	float fraction;
	float centroid, centroid_unc;
	float min, max, min_unc, max_unc;
	float fwhm;
	char st[50];
	
	trans_unc = vector(1, n_wl);
	
	for(i = 1; i <= n_wl; i++){
		trans_unc[i] = 0;
	}
	
	fraction = 0.5;
	find_limits_fr(i_pk, n_wl, fraction, trans, &i_min, &i_max);
	centroid_calc(i_min, i_max, trans, &centroid, &centroid_unc);
	printf("centroid\t%f\t+-%f\n", centroid, centroid_unc);
	fprintf(fp_out, "%.2f\t", centroid);
	
	fraction = 0.25;
	find_limits_fr(i_pk, n_wl, fraction, trans, &i_min, &i_max);
	centroid_calc(i_min, i_max, trans, &centroid, &centroid_unc);
	printf("centroid\t%f\t+-%f\n", centroid, centroid_unc);
	fprintf(fp_out, "%.2f\t", centroid);
	
	fraction = 0.1;
	find_limits_fr(i_pk, n_wl, fraction, trans, &i_min, &i_max);
	centroid_calc(i_min, i_max, trans, &centroid, &centroid_unc);
	printf("centroid\t%f\t+-%f\n", centroid, centroid_unc);
	fprintf(fp_out, "%.2f\t", centroid);
	
	find_limits_bw(i_pk, n_wl, bw*1.5, trans, &i_min, &i_max);
	centroid_calc(i_min, i_max, trans, &centroid, &centroid_unc);
	printf("centroid\t%f\t+-%f\n", centroid, centroid_unc);
	fprintf(fp_out, "%.2f\t", centroid);
	
	centroid_calc(1, n_wl, trans, &centroid, &centroid_unc);
	fprintf(fp_out, "%.2f\t", centroid);
	
	fwhm_calc(i_pk, trans, &min, &max, &min_unc, &max_unc);
	centroid = (max + min)/2.0;
	printf("centroid\t%f\t+-%f\n", centroid, centroid_unc);
	fprintf(fp_out, "%.2f\t", centroid);
	
	fwhm = max - min;
	printf("bandwidth\t%f\t+-%f\n", fwhm);
	fprintf(fp_out, "%.2f\t", fwhm);
	
	fprintf(fp_out, "\n");
	
	free_vector(trans_unc, 1, n_wl);
	
	return 1;
}