#include "track.h"
float search(fm,qsum,tphi,tpsi,tgamph,tgamps,g,h,freq,n)
float *fm,qsum;
float tphi[50][5][18],tpsi[50][6][18],tgamph[50][5];
float tgamps[50][6],g[],h[],freq[50];
int n;
{
	float funmin = 1.e10,filfun(); 
	float fun[50];
	float f1,f2,f3,x0,x1,x2,x3,a,b,c,ftemp;
	int istar,i;
	filfun(tphi,tgamph,tpsi,tgamps,fun,freq,g,h,qsum,n);
	istar =  -1;
	for(i=0 ; i<50 ; ++i)
	{
		ftemp=fun[i];
		if(ftemp < funmin)
		{
			funmin = ftemp;
			istar = i;
		}
	}
	if( istar < 0 )
	{
		fprintf(stderr,"ptrack: search: error in search\n");
		exit();
	}
	if( istar == 0 || istar == 49 )
	{
		*fm = fun[istar];
		return (freq[istar]);
	} else {
		x1 = freq[istar-1];
		f1 = fun[istar-1];
		x2 = freq[istar];
		f2 = fun[istar];
		x3 = freq[istar+1];
		f3 = fun[istar+1];
		a = f3/((x3-x1)*(x3-x2));
		b = f2/((x2-x1)*(x2-x3));
		c = f1/((x1-x2)*(x1-x3));
		x0 = .5*(a*(x1+x2)+b*(x1+x3)+c*(x2+x3))/(a+b+c);
		*fm = a*(x0-x1)*(x0-x2)+b*(x0-x1)*(x0-x3)+c*(x0-x2)*(x0-x3);
		return(x0);
	}
}
