program errell
c
c     This program reads the output from errbin (subroutine chibndry)
c     and prepares a plot of the error ellipse in del Dec vs. del RA. 
c     It also fits an ellipse to the data and prints out the estimates 
c     of major and minor axis, as well as the position angle.
c
      parameter(maxfiles=3)
      character*64 string
      character*12 moddsn(maxfiles),datdsn(maxfiles)
      character*80 xtext,ytext,text
      character*80 data
      character*8 starnam(maxfiles)
      parameter (pi=3.141592654)
      parameter (degrad=pi/180.)
      parameter (maxp=8,nterms=5)
      real x(maxp),y(maxp),r(maxp),p(maxp)
      real xf(18),yf(18)
      real dx(360),dy(360)
      real sigmax(maxp),sigmay(maxp),sigmar(maxp)
      real rfit(maxp)
      real a(5),deltaa(5),sigmaa(5)
      common/data/x,y,sigmax,sigmay
c
      iun=0
      do iunit=1,maxfiles
   10  write(6,1)iunit
    1  format(' Enter name of file',i2,' ! ',$)
       read(5,2,end=11)string
    2  format(a64)
       open(iunit,file=string,status='old',err=10)
       iun=iun+1
      enddo
   11 continue
c     write(6,*)' Enter variable parameter id to select!'
c     read(5,*)id
c
c     Now read the files to get the error ellipse points
c
      num=0
      numf=0
      xmax=0
      ymax=0
      do 100 iunit=1,iun
       string='Input data file'
       call findnstr(iunit,string,data,ier,ipos)
       if(ier.ne.0)then
       write(6,3)iunit
    3  format(' Did not find header of file',i2,' ! Stop now!')
       stop
       endif
       ipos1=ipos+19
       ipos2=64
       read(data(ipos1:ipos2),'(a)')datdsn(iunit)
       string='Final model'
       call findnstr(iunit,string,data,ier,ipos)
       ipos1=ipos+19
       ipos2=ipos1+11
       read(data(ipos1:ipos2),'(a)')moddsn(iunit)
       string='Source'
       call findnstr(iunit,string,data,ier,ipos)
       ipos1=ipos+8
       ipos2=ipos1+7
       read(data(ipos1:ipos2),'(a)')starnam(iunit)
       if(iunit.gt.1)then
        if(datdsn(iunit).ne.datdsn(iunit-1))goto 13
        if(moddsn(iunit).ne.moddsn(iunit-1))goto 13
        if(starnam(iunit).ne.starnam(iunit-1))goto 13
        goto 12
   13   write(6,*)' Not the same data or model or star!! Stop now'
        stop
       endif
   12  continue
       string='Final model'
       call findnstr(iunit,string,data,ier,ipos)
       call readnrt(iunit,rbest,tbest,ier)
       xbest=rbest*sin(tbest*degrad)
       ybest=rbest*cos(tbest*degrad)
       string='PARM1D'
       do k=1,30
       call findnstr(iunit,string,data,ier,ipos)
       if(ier.ne.0)goto 100
       ipos1=ipos+7
       ipos2=ipos1+2
       read(data(ipos1:ipos2),*)iparm 
       if(iparm.eq.2.or.iparm.eq.8)then
        call readnrt(iunit,radius,theta)
        num=num+1
        x(num)=radius*sin(theta*degrad)-xbest
        y(num)=radius*cos(theta*degrad)-ybest
        xmax=max(xmax,x(num))
        ymax=max(ymax,y(num))
       else
        call readnrt(iunit,radius,theta)
        numf=numf+1
        xf(numf)=radius*sin(theta*degrad)-xbest
        yf(numf)=radius*cos(theta*degrad)-ybest
       endif
       enddo
  100 enddo
c
c     Now you've read the ivar ellipse points, determine major and minor
c     axis and position angle of the major axis.
c
      mode=0
      chisqr=0.
c--open and read data file with npts points
      npts=num
      write(text,4)npts
    4 format(i2,' points read')
      call putout(text)
c--set start values of parameters
c     guess for maj. axis of ellipse and increment?
      a(4)=xmax
      deltaa(4)=a(4)/20.
c     guess for min. axis of ellipse and increment?
      a(5)=ymax
      deltaa(5)=a(5)/20.
c     guess for x-coord. of center and increment?
      a(1)=0
      deltaa(1)=deltaa(4)
c     guess for y-coord. of center and increment?
      a(2)=0
      deltaa(2)=deltaa(5)
c     guess for pos. angle of ellipse and increment?
      a(3)=0
      deltaa(3)=5.
      a(3)=a(3)*degrad
      deltaa(3)=deltaa(3)*degrad
c     fractional change of chisqr for convergence? 
      chifr=0.001
      call fixell(p,r,npts,a,rfit)
      nfree=npts-nterms
      chiold=fchisq(r,sigmar,npts,nfree,mode,rfit)
      write(text,5)chiold
    5 format('Value of chisqr: ',f5.1)
c     call putout(text)
      do iter=1,50
       call gridls(p,r,sigmar,npts,nterms,mode,
     &  a,deltaa,sigmaa,rfit,chisqr)
       write(text,6)a(1),a(2),a(3)/degrad,a(4),a(5),chisqr
c      call putout(text)
       if(1.-chisqr/chiold.lt.chifr)goto 111
       chiold=chisqr
      enddo
  111 write(text,6)a(1),a(2),a(3)/degrad,a(4),a(5),chisqr
    6 format('x= ',f5.3,' y= ',f5.3,' pa= ',f5.1,' a= ',f4.2,' b= ',
     & f4.2,' chisqr= ',f7.5)
      call putout(text)
      write(xtext,7)xbest
    7 format('RA - ',f6.2,' (mas)')
      write(ytext,8)ybest
    8 format('DEC - ',f6.2,' (mas)')
      write(text,9)starnam(1),moddsn(1),a(4),a(5),90.-a(3)/degrad
    9 format('Error ellipse for ',a8,', Model: ',a12,', a=',f4.2,', b=',
     & f4.2,', pa=',f5.1)
c--centx=a1,centy=a2,position angle=a3,major axis=a4,minor axis=a5
      esqr=1.-a(5)*a(5)/(a(4)*a(4))
      dxmax=0
      dymax=0
      do k=1,360
       phi=float(k)*degrad
       radius=sqrt(a(5)*a(5)/(1.-esqr*cos(phi)*cos(phi)))
       dx(k)=radius*cos(phi+a(3))+a(1)
       dy(k)=radius*sin(phi+a(3))+a(2)
       dxmax=max(dxmax,abs(dx(k)))
       dymax=max(dymax,abs(dy(k)))
      enddo
c
c     Plot the ellipse, also make hardcopy
c
      do iplot=1,2
      call pgbegin(0,'?',1,1)
      xplmin=-dxmax-dxmax/5.
      xplmax=-xplmin
      yplmin=-dymax-dymax/5.
      yplmax=-yplmin
      call pgsci(4)
      call pgenv(xplmin,xplmax,yplmin,yplmax,0,2)
      call pgsci(5)
      call pglabel(xtext,ytext,text)
      call pgsci(7)
      call pgpoint(num,x,y,17)
      call pgsci(3)
      call pgpoint(numf,xf,yf,21)
      call pgsci(9) 
      call pgline(360,dx,dy)
      call pgend
      write(6,*)' Now select /lpm to get hardcopy !'
      enddo
      end
c
      subroutine findnstr(iunit,string,data,ier,ipos)
      character*64 string
      character*80 data
      ier=0
      do i=1,100
       read(iunit,1,end=111)data
    1  format(a80)
       ipos=index(data,string(1:len1(string)))
       if(ipos.ne.0)return
      enddo
  111 ier=1
      return
      end
c
      subroutine readnrt(iunit,r,t)
      character*64 string
      character*80 data
      string='Lambda'
      call findnstr(iunit,string,data,ier,ipos)
      read(iunit,*)
      read(iunit,*)ncomp1,wavel1
      do j=1,3
       read(iunit,*,end=111,err=111)ncomp2,wavel2
       if(wavel2.eq.wavel1.and.ncomp2.eq.2)then
       backspace iunit
       read(iunit,*)ncomp2,wavel2,r,t
       goto 112
       endif
      enddo
  111 write(6,8)iunit
    8 format(' Something is wrong with the data file ',i2)
      stop
  112 return
      end