#include "track.h"

float
trigpo(omega,phi,psi,gampsi,gamphi,n)
float omega,psi[6][18],phi[5][18];
float gamphi[5],gampsi[6];
int n;
{
	int j=0;
	int k,np;
	double alpha,beta,gamma;
	double wcos[18],wsin[18];
	double p,z,a,b,yy;
	double sin(),cos();
	np=n+1;
	for (k=0 ; k<MM ; ++k)
	{
		yy = omega* (float) (k);
		wcos[k]=cos(yy);
		wsin[k]=sin(yy);
	}
	beta=gamma= ZERO;
	for (k=0 ; k<MM ; ++k)
	{
		p=wsin[k];
		z=p*p;
		beta+=z*wcos[k];
		gamma+=z;
		phi[0][k]=p;
	}
	gamphi[0]=gamma;
	a=TWO * beta/gamma;
	alpha=beta=gamma=ZERO;
	for (k=0 ; k<MM ; ++k)
	{
		p= (TWO * wcos[k]-a) * phi[0][k];
		alpha += wcos[k] * p * phi[0][k];
		beta += wcos[k]* ( p * p ) ;
		gamma += ( p * p );
		phi[1][k] = p ;
	}
	gamphi[1]=gamma;
	a = TWO * beta/gamma;
	b = TWO *alpha/gamphi[0];
	for (j=2 ; j<n ; ++j)
	{
		alpha=beta=gamma = ZERO;
		for (k=0 ; k< MM ; ++k)
		{
			p = (TWO * wcos[k] - a ) * phi[j-1][k] - b * phi[j-2][k];
			alpha += wcos[k] * p * phi[j-1][k];
			beta += wcos[k] * (p * p);
			gamma += (p * p) ;
			phi[j][k] = p ;
		}
		gamphi[j] = gamma;
		a = TWO * beta/gamma;
		b = TWO *alpha/gamphi[j-1];
	}
	beta = ZERO ;
	gamma = (double) MM;
	for ( k=0; k < MM ; ++k)
	{
		beta += wcos[k];
		psi[0][k]= ONE;
	}
	gampsi[0]=gamma;
	a=beta/gamma;
	alpha=beta=gamma= ZERO;
	for ( k=0 ; k < MM ; ++k)
	{
		p = wcos[k]-a;
		alpha += wcos[k] * p*psi[0][k];
		beta += wcos[k] * ( p * p );
		gamma += (p * p);
		psi[1][k] = p ;
	}
	gampsi[1] = gamma ;
	a = TWO * beta / gamma ;
	b = TWO * alpha / gampsi[0];
	for (j=2 ; j<np ;++j)
	{
		alpha=beta=gamma= ZERO;
		for ( k=0; k < MM ; ++k)
		{

			p=(TWO * wcos[k]-a)* psi[j-1][k]-b*psi[j-2][k];
			alpha += wcos[k]*p*psi[j-1][k];
			beta += wcos[k]* (p*p);
			gamma += (p*p);
			psi[j][k] = p;
		}
		gampsi[j]=gamma;
		a=TWO*beta/gamma;
		b=TWO*alpha/gampsi[j-1];
	}
	return;
}