/* Jacobi polynomials */

float Jacobi(alpha,beta,q,x)
float alpha,beta,x;
int q;
{
double ab,ab1,ab2,a2b2,y,ym1,di,a1,a2,a3,a4,yp1;
int i;
float yy;

  ab=alpha+beta;
  ab1=ab+1.0;
  ab2=ab+2.0;
  a2b2=alpha*alpha-beta*beta;
  y=1.0;
  ym1=0.0;
  for(i=0;i<q;i++)
  { di=(double)i;
    a1=2.0*(di+1.0)*(di+ab1)*(2.0*di+ab);
    a2=(2.*di+ab1)*a2b2;
    a3=(2.*di+ab)*(2.*di+ab1)*(2.*di+ab2);
    a4=2.0*(di+alpha)*(di+beta)*(2.*di+ab2);
    yp1=((a2+a3*x)*y-a4*ym1)/a1;
    ym1=y;
    y=yp1;
  }
  yy=(float)y;
  return yy;
  
}