subroutine polfit
      common /fit/x,y,sigmay,npts,mode
      common /pfit/nterms,a,chisqr
      double precision x,y,sigmay,a,chisqr
      double precision sumx,sumy,xterm,yterm,array,chisq
      integer npts,mode
      dimension x(50),y(50),sigmay(50),a(10)
      dimension sumx(19),sumy(10),array(10,10)
11    nmax=2*nterms-1
      do 13 n=1,nmax
13    sumx(n)=0.
      do 15 j=1,nterms
15    sumy(j)=0.
      chisq=0.
21    do 50 i=1,npts
      xi=x(i)
      yi=y(i)
31    if (mode) 32,37,39
32    if (yi) 35,37,33
33    weight=1./yi
      goto 41
35    weight=1./(-yi)
      goto 41
37    weight=1.
      goto 41
39    weight=1./sigmay(i)**2
41    xterm=weight
      do 44 n=1,nmax
      sumx(n)=sumx(n)+xterm
44    xterm=xterm*xi
45    yterm=weight*yi
      do 48 n=1,nterms
      sumy(n)=sumy(n)+yterm
48    yterm=yterm*xi
49    chisq=chisq+weight*yi**2
50    continue
51    do 54 j=1,nterms
      do 54 k=1,nterms
      n=j+k-1
54    array(j,k)=sumx(n)
      detla=determ(array,nterms)
      if (delta) 61,57,61
57    chisqr=0.
      do 59 j=1,nterms
59    a(j)=0.
      goto 80
61    do 70 l=1,nterms
62    do 66 j=1,nterms
      do 65 k=1,nterms
      n=j+k-1
65    array(j,k)=sumx(n)
66    array(j,l)=sumy(j)
70    a(l)=determ(array,nterms)/delta
71    do 75 j=1,nterms
      chisq=chisq-2.*a(j)*sumy(j)
      do 75 k=1,nterms
      n=j+k-1
75    chisq=chisq+a(j)*a(k)*sumx(n)
76    free=npts-nterms
77    chisqr=chisq/free
80    return
      end