      subroutine errell(npts,r,p,a)
c
c     Basic subroutine fits ellipse to positions (r,p)
c
      real*8 pi,degrad,mas2rad
      parameter (pi=3.141592653589793D0)
      parameter (degrad=pi/180.D0)
      parameter (mas2rad=pi/(180.D0*3600.D0*1000.D0))

      integer*4 npts,nangles
      real*8 r(*),p(*)

      parameter (nangles=36)
      real*8 ex(nangles),ey(nangles)
      common/edata/ex,ey

      integer*4 maxit
      parameter (maxit=50)   ! maximum number of iterations

      real*8 a(5),deltaa(5),sigmaa(5)
      real*8 chisqr,chiold,chifr,xmax,xmin,ymax,ymin,dist,disto
      real*8 sigmar(nangles),rfit(nangles)
c
      mode=1    ! tells GRIDLS to call fixell instead of rtfuncs
      modef=0   ! Gaussian model in FCHISQ
      chisqr=0.D0
      nterms=5
      nfree=npts-nterms
c--check for bad values
      do i=1,npts
       if(r(i).eq.4.848d-9)then   ! indicator for problems in rterr
        do j=1,npts
         if(r(j).ne.4.848d-9)then
          r(i)=r(j)
          p(i)=p(j)
         endif
        enddo
       endif
      enddo        
c--set start values of parameters
      xmax=r(1)*cos(p(1))
      xmin=xmax
      ymax=r(1)*sin(p(1))
      ymin=ymax
      do i=1,npts
       ex(i)=r(i)*cos(p(i))
       xmax=max(xmax,ex(i))
       xmin=min(xmin,ex(i))
       ey(i)=r(i)*sin(p(i))
       ymax=max(ymax,ey(i))
       ymin=min(ymin,ey(i))
      enddo
      a(1)=(xmin+xmax)/2.D0   ! x-coordinate of center
      a(2)=(ymin+ymax)/2.D0   ! y-coordinate of center
      deltaa(1)=sqrt((xmax-xmin)*(ymax-ymin))/10.D0
      deltaa(2)=deltaa(1)
c     guess for pos. angle of ellipse and increment?
      disto=0.D0
      do i=1,npts
       do j=1,npts
        dist=sqrt((ex(i)-ex(j))**2+(ey(i)-ey(j))**2)
        if(dist.gt.disto)then
         imax=i
         jmax=j
         disto=dist
        endif
       enddo
      enddo
      a(3)=atan2(ey(imax)-ey(jmax),ex(imax)-ex(jmax))
      deltaa(3)=5.D0*degrad
c     guess for axes and increments?
      xmax=ey(1)*sin(a(3))+ex(1)*cos(a(3))
      xmin=xmax
      ymax=ey(1)*cos(a(3))-ex(1)*sin(a(3))
      ymin=ymax
      do i=1,npts
       xmax=max(xmax,ey(i)*sin(a(3))+ex(i)*cos(a(3)))
       ymax=max(ymax,ey(i)*cos(a(3))-ex(i)*sin(a(3)))
       xmin=min(xmin,ey(i)*sin(a(3))+ex(i)*cos(a(3)))
       ymin=min(ymin,ey(i)*cos(a(3))-ex(i)*sin(a(3)))
      enddo
      a(4)=(xmax-xmin)/2.D0    ! major half axis
      deltaa(4)=a(4)/10.D0
      a(5)=(ymax-ymin)/2.D0    ! minor half axis
      deltaa(5)=a(5)/10.D0
c     set uncertainties to a mean value
      do i=1,npts
       sigmar(i)=deltaa(4)
      enddo
c--do the iterations
c     fractional change of chisqr for convergence? 
      chifr=0.001D0
      call fixell(p,r,npts,a,rfit)
      call fchisq(r,sigmar,npts,nfree,rfit,chiold,modef)
      do iter=1,maxit
       call gridls(p,r,sigmar,npts,nterms,mode,
     &             a,deltaa,sigmaa,rfit,chisqr)
       if(1.D0-chisqr/chiold.lt.chifr.and.iter.gt.2)goto 111
       chiold=chisqr
      enddo
  111 if(a(3).lt.0.D0)a(3)=a(3)+2.D0*pi
      return
      end

