program ellfit
      parameter (maxp=100,nterms=5)
      real x(maxp),y(maxp),r(maxp),p(maxp)
      real sigmax(maxp),sigmay(maxp),sigmar(maxp)
      real rfit(maxp)
      real a(5),deltaa(5),sigmaa(5)
      character*80 text,answer
      common/data/x,y,sigmax,sigmay
      pi=3.141592654
      degrad=pi/180.
      mode=0
      chisqr=0.
c--open and read data file with npts points
   10 call questn('file with ellipse points? ',answer)
      open(1,file=answer,status='OLD',err=11)      
      do i=1,maxp
       read(1,*,err=11,end=12)x(i),y(i)
      enddo
   11 call putout('*** File does not exist or is corrupt! ')
      goto 10
   12 npts=i-1
      write(text,2)npts
    2 format(i2,' points read')
      call putout(text)
c--get start values of parameters
      call questn('guess for x-coord. of center and increment? ',
     &                  answer)
      read(answer,*)a(1),deltaa(1)
      call questn('guess for y-coord. of center and increment? ',
     &                  answer)
      read(answer,*)a(2),deltaa(2)
      call questn('guess for pos. angle of ellipse and increment? ',
     &                  answer)
      read(answer,*)a(3),deltaa(3)
      a(3)=a(3)*degrad
      deltaa(3)=deltaa(3)*degrad
      call questn('guess for maj. axis of ellipse and increment? ',
     &                  answer)
      read(answer,*)a(4),deltaa(4)
      call questn('guess for min. axis of ellipse and increment? ',
     &                  answer)
      read(answer,*)a(5),deltaa(5)
      call questn('fractional change of chisqr for convergence? ',
     &                  answer)
      read(answer,*)chifr
      call fixell(p,r,npts,a,rfit)
      nfree=npts-nterms
      chiold=fchisq(r,sigmar,npts,nfree,mode,rfit)
      write(text,4)chiold
    4 format('Value of chisqr: ',f5.1)
      call putout(text)
      do iter=1,50
       call gridls(p,r,sigmar,npts,nterms,mode,
     &  a,deltaa,sigmaa,rfit,chisqr)
       write(text,3)a(1),a(2),a(3)/degrad,a(4),a(5),chisqr
       call putout(text)
       if(1.-chisqr/chiold.lt.chifr)goto 111
       chiold=chisqr
      enddo
  111 write(text,3)a(1),a(2),a(3)/degrad,a(4),a(5),chisqr
    3 format('x= ',f5.1,' y= ',f5.1,' pa= ',f5.1,' a= ',f4.1,' b= ',
     & f4.1,' chisqr= ',f6.3)
      call putout(text)
      open(2,file='ellfit.out')
      do i=1,npts
      dx=rfit(i)*cos(p(i)+a(3))+a(1)
      dy=rfit(i)*sin(p(i)+a(3))+a(2)
      write(2,*)x(i),dx,y(i),dy
      enddo
      end