subroutine rtfit
c
c     Consists of two parts. The first one fits radius and position angle 
c     of binary to data of individual nights. Gets initial values from the 
c     binary model. The second part plots fit to the observed visibilities.
c
      include 'constnts.inc'
      include 'data.inc'
      include 'model.inc'
c
      integer*4 nangles
      parameter (nangles=36)

      real*8 xdata(mxscn2),ydata(mxscn2),sigma(mxscn2),yfit(mxscn2)
      real*8 p(2),dp(2),sigmap(2),chisq,chisqo,chif,chifr,dangle
      real*8 er(nangles),ep(nangles),ea(5)
      integer*4 ndate(mxdates)

      real*8 modvis2(mxscn,mxfilt),modphas(mxscn,mxfilt)
      real*8 u(mxscn,mxfilt),v(mxscn,mxfilt)
      real*8 mhours(mxscn),lambda(mxscn,mxfilt),dlambda(mxscn,mxfilt)
      real*8 cchisq(mxfilt)
      integer*4 mdate(mxscn),mbaseline(mxscn),mstars(mxscn)
      integer*4 nchisq(mxfilt)

c     These are defined for plotting purposes
      real*4 y(mxpl),t(mxpl),vis2max(mxfilt),vis2min(mxfilt)
      real*4 timmax,timmin,timstep
      integer*4 firstid,lastid
      logical docalc
      character header*80,device*4,plotfile*15
      character orbsense*8,starn*8
      character fknam1*1,fknam2*2,fknam3*3,fknam4*4

c     These are defined for plotting the error ellipses
      real*4 ex4(nangles),ey4(nangles),dx(mxepl),dy(mxepl)
      real*8 radius,phi,phistep,esqr

c     These are defined for subroutine true2app
      real*8 jd,djul,a(7),dyda(7,2),ra,dec
      real*8 seprad(mxdates),parad(mxdates)
      integer*4 pr
c
      pr=6
      a(1)=a_true*mas2rad
      a(2)=e_true
      a(3)=period
      a(4)=epoch
      a(5)=inclin*degrad
      a(6)=asc_node*degrad
      a(7)=arg_peri*degrad
c
      do j=1,mxdates
       ndate(j)=0
      enddo
c
      write(outc,*)'Please enter your choice: '
      write(outc,*)'-1- Fit (r,theta) values'
      write(outc,*)'-2- Read from file and plot'
      write(outc,'(a,$)')' >'
      read(inc,2)device
      if(device.eq.'2')then
       docalc=.false.
       goto 1000
      endif
      docalc=.true.
c
c     The fitted r-theta models will be plotted, prepare plot here
c
      write(outc,*)'Screen (s) or hardcopy (h,p,v) or exit (n) ?'
      read(inc,2)device
    2 format(a4)
      if(device.eq.'h'.or.device.eq.'H')then
       plotfile='plotrte.lpm'//'/lpm'
      else if(device.eq.'p'.or.device.eq.'P')then
       plotfile='plotrte.ps'//'/ps'
      else if(device.eq.'v'.or.device.eq.'V')then
       plotfile='plotrte.ps'//'/vps'
      else if(device.eq.'s'.or.device.eq.'S')then
       plotfile='/xs'
      else
       return
      endif
c
c     Open .DAT file for output of r-theta values and error ellipses
      write(fknam4,'(i4)')starno
      if(log10(real(starno)).lt.1)then
       write(fknam1,'(i1)')starno
       open(4,file=fknam1//'.dat')
      else if(log10(real(starno)).lt.2)then
       write(fknam2,'(i2)')starno
       open(4,file=fknam2//'.dat')
      else if(log10(real(starno)).lt.3)then
       write(fknam3,'(i3)')starno
       open(4,file=fknam3//'.dat')
      else
       open(4,file=fknam4//'.dat')
      endif
      if(inclin.gt.90.D0)then
       orbsense='RETROGRADE'
      else
       orbsense='PROGRADE'
      endif
      starn=star
      do i=1,8
       if(starn(i:i).eq.' ')starn(i:i)='_'
      enddo
      write(4,4)starn,starno,period,orbsense
    4 format(a8,1x,i4,2x,f8.2,2x,a8) 
      call pgbegin(0,plotfile,1,1)   ! Initialize plot of error ellipses
      call pgask(.false.)            ! Advance plots automatically
c 
c Loop over all days and tweak r and theta using a grid search
c
      firstid=1
      nd=0
      DO 100 ID=1,NDATA      ! To do this, go through all data
       IF(DATE(ID).NE.DATE(ID+1).or.baseline(id+1).ne.baseline(id)
     &                          .OR.ID+1.GT.NDATA)THEN
        lastid=id
        ndat=lastid-firstid+1
        if(ndat.lt.3)then
         write(outc,*)
         write(outc,*)'Too few data points! Skip this night!'
         goto 101
        endif
        nd=nd+1
        if(ndat*3.gt.mxscn2)then
         write(outc,*)' RTFIT: Array too small!'
         write(outc,*)' Increase mxpl and recompile!'
         stop
        endif
        do ig=1,nfilt
         do iid=1,ndat
          xdata((ig-1)*ndat+iid)=dble(iid+firstid-1) !pass this index for
c                                                    identification purposes,
c                                                    xdata is not used in GRIDLS
          ydata((ig-1)*ndat+iid)=obsvis2(iid+firstid-1,ig)
          sigma((ig-1)*ndat+iid)=viserr(iid+firstid-1,ig)
         enddo
        enddo                           
        ndat3=nfilt*ndat
c
        IDAY  = MOD(DATE(ID),       100)
        IMONTH= MOD(DATE(ID)/100,   100)
        IYEAR = MOD(DATE(ID)/10000, 100) + 1900
        ihour=8 ! Midnight on Mt Wilson
        djul=(jd(iday,imonth,iyear)-2440000.D0)+dble(ihour)/24.D0
        call true2app(djul,a,seprad(nd),parad(nd),ra,dec,dyda,pr)
        p(1)=seprad(nd)
        p(2)=parad(nd)
c
c       Prepare for GRIDLS...
        nterms=2
        mode=0      ! tells GRIDLS to call rtfuncs instead of fixell
        dp(1)=p(1)/20.D0
        dp(1)=max(dp(1),0.1D0*mas2rad)
        dp(1)=min(dp(1),1.0D0*mas2rad)
        dp(2)=1.0D0*degrad
        chisqo=1.D9
        chifr=0.0001D0
c
c       ...and tweak binary position
        do iter=1,maxitg
         call gridls(xdata,ydata,sigma,ndat3,nterms,mode,p,dp,
     &               sigmap,yfit,chisq)
         p(1)=abs(p(1))         ! sometimes, radius gets negative
         p(2)=mod(p(2),twopi)   ! index 2: these are all angles
         dp(2)=mod(dp(2),twopi)
         sigmap(2)=mod(sigmap(2),twopi)
         if(1.D0-chisq/chisqo.lt.chifr.and.iter.gt.4)then
          write(outc,*)'RTFIT: number of iterations made:',iter
          goto 111
         endif
         chisqo=chisq
        enddo
        write(outc,*)'RTFIT: did not converge!'
        write(outc,*)'Continuing...'
  111   seprad(nd)=p(1)
        parad(nd)=p(2)
        ndate(nd)=date(id)
c
c       Now that you have found best-fit r and theta, determine error ellipse..
c       ...but first, we want to scale the visibility errors so that X^2=1
c
        do i=1,ndat3
         sigma(i)=sigma(i)*sqrt(chisq)
        enddo
        sigmap(1)=sigmap(1)*sqrt(chisq)
        sigmap(2)=sigmap(2)*sqrt(chisq)
        chif=min(1.d0,sqrt(chisq))
c       chif=sqrt(chisq)
c
        dp(1)=sqrt(sigmap(1) * sigmap(2)*p(1)) ! geometric average 2 axes
        dp(1)=max(dp(1),mas2rad/100.D0)        ! prevent too small a value
        dp(2)=1.D0
        dangle=2.D0*pi/dble(nangles)
c
        write(outc,*)'Determining Chi^2 boundary...'
        do iangle=1,nangles
         ep(iangle)=dp(2)+pi/2.D0  ! convention in errell is the standard one
         if(ep(iangle).gt.2.D0*pi)ep(iangle)=ep(iangle)-2.D0*pi
         call rterr(xdata,ydata,sigma,ndat3,p,dp,yfit)
         er(iangle)=dp(1)/chif
         ex4(iangle)=real(er(iangle)*cos(ep(iangle))/mas2rad)
         ey4(iangle)=real(er(iangle)*sin(ep(iangle))/mas2rad)
         dp(2)=dp(2)+dangle
        enddo
        write(outc,*)'Fitting error ellipse...'
        call errell(nangles,er,ep,ea)
c       We want to calculate the number of seconds since UT=0 hours
        utsec=0.
        do iid=firstid,lastid
         utsec=utsec+real(hours(iid))
        enddo
        utsec=utsec/(lastid-firstid+1)
        utsec=8. ! Midnight on Mt. Wilson (this is the correct choice)
        utsec=utsec*3600.

c Now write .DAT file in ELLIPSE-format; this will also be read by PLOTORB
        write(4,5)seprad(nd)/mas2rad,parad(nd)/degrad,
     &   ea(4)/mas2rad,ea(5)/mas2rad,ea(3)/degrad-90.,
     &   date(id),utsec,ea(1)/mas2rad,ea(2)/mas2rad
    5   format(f6.2,1x,f7.2,2x,f6.3,1x,f6.3,1x,f6.1,2x,i6,1x,f6.0,
     &         3x,f5.2,1x,f5.2)

c Now plot the ellipse; advance plot automatically
c--centx=a1,centy=a2,position angle=a3,major axis=a4,minor axis=a5
        esqr=1.D0-ea(5)*ea(5)/(ea(4)*ea(4))
        dxmax=0
        dymax=0
        dxmin=0
        dymin=0
        phistep=360.D0*degrad/dble(mxepl-1)
        do k=1,mxepl
         phi=dble(k)*phistep
         radius=sqrt(ea(5)*ea(5)/(1.D0-esqr*cos(phi)*cos(phi)))
         dx(k)=real((radius*cos(phi+ea(3))+ea(1))/mas2rad)
         dy(k)=real((radius*sin(phi+ea(3))+ea(2))/mas2rad)
         dxmax=max(dxmax,dx(k))
         dxmin=min(dxmin,dx(k))
         dymax=max(dymax,dy(k))
         dymin=min(dymin,dy(k))
        enddo
        dxmax=dxmax+(dxmax-dxmin)/10.
        dxmin=dxmin-(dxmax-dxmin)/10.
        dymax=dymax+(dymax-dymin)/10.
        dymin=dymin-(dymax-dymin)/10.
c      
c       Plot the ellipse
c      
        call pgsci(4)
        call pgsch(1.5)
        call pgenv(dxmin,dxmax,dymin,dymax,1,2)
        call pgsci(6)
        write(header,6)date(id),p(1)/mas2rad,p(2)/degrad,chisq
    6   format('Error ellipse on ',i6,', fitted r=',f5.1,
     &         ', theta=',f6.1,', \gx\u2\d=',f5.2)
        call pglabel('Right Ascension (mas)','Declination (mas)',
     &               header)
        call pgsci(7)
        call pgpoint(nangles,ex4,ey4,2)
        call pgsci(9) 
        call pgline(mxepl,dx,dy)
c
c       Prepare for next night
  101   firstid=lastid+1
c
       ENDIF
  100 CONTINUE

      close(4)
      call pgend

 1000 if(.not.docalc)then
       write(fknam4,'(i4)')starno
       if(log10(real(starno)).lt.1)then
        write(fknam1,'(i1)')starno
        open(4,file=fknam1//'.dat',status='old',err=52)
       else if(log10(real(starno)).lt.2)then
        write(fknam2,'(i2)')starno
        open(4,file=fknam2//'.dat',status='old',err=52)
       else if(log10(real(starno)).lt.3)then
        write(fknam3,'(i3)')starno
        open(4,file=fknam3//'.dat',status='old',err=52)
       else
        open(4,file=fknam4//'.dat',status='old',err=52)
       endif
       read(4,*,end=52)     ! skip header
       do j=1,mxdates 
        read(4,*,end=53)seprad(j),parad(j),emaj,emin,epa,ndate(j)
        seprad(j)=seprad(j)*mas2rad
        parad(j)=parad(j)*degrad
       enddo
       goto 53
   52  write(outc,'(a,a,a)')' File: ',fknam4,'.dat  missing or empty.'
       close(4)
       return
   53  close(4)
      endif
c
c     Initialize plotting of data --------------------------------------
c
      write(outc,*)'Screen (s) or hardcopy (h,p,v) or exit (n) ?'
      read(inc,2)device
      if(device.eq.'h'.or.device.eq.'H')then
       plotfile='plotrtf.lpm'//'/lpm'
       iy=4
       csize=2.5
      else if(device.eq.'p'.or.device.eq.'P')then
       plotfile='plotrtf.ps'//'/ps'
       iy=4
       csize=2.5
      else if(device.eq.'v'.or.device.eq.'V')then
       plotfile='plotrtf.ps'//'/vps'
       iy=4
       csize=2.5
      else if(device.eq.'s'.or.device.eq.'S')then
       plotfile='/xs'
       iy=2
       csize=2.5
      else
       return
      endif
      call pgbegin(0,plotfile,3,iy)
      call pgask(.true.)
      call pgsch(csize)
c
      firstid=1
      lastid=1
      do ig=1,nfilt
       vis2max(ig)=0
       vis2min(ig)=0
      enddo
      timmin=real(hours(1))
      timmax=timmin
c
      DO 200 ID=1,NDATA
       do ig=1,nfilt
        if(viserr(id,ig).ne.2)then
         vis2max(ig)=max(vis2max(ig),real(obsvis2(id,ig)))
         vis2min(ig)=min(vis2min(ig),real(obsvis2(id,ig)))
        endif
       enddo
       timmax=max(timmax,real(hours(id)))
       timmin=min(timmin,real(hours(id)))
c
       IF(DATE(ID).NE.DATE(ID+1).or.baseline(id+1).ne.baseline(id)
     &                          .OR.ID+1.GT.NDATA)THEN
        lastid=id
        ndat=lastid-firstid+1

        timr=timmax-timmin
        timmin=timmin-timr/10.
        timmax=timmax+timr/10.

        nd=0
        do j=1,mxdates
         if(ndate(j).eq.date(id))nd=j
        enddo
        if(nd.ne.0)then
         do iid=1,ndat
          mstars(iid)=stars(firstid+iid-1)
          mdate(iid)=date(firstid+iid-1)
          mbaseline(iid)=baseline(firstid+iid-1)
          mhours(iid)=hours(firstid+iid-1)
          do ig=1,nfilt
           lambda(iid,ig)=lambda0(firstid+iid-1,ig)
           dlambda(iid,ig)=dlambda0(firstid+iid-1,ig)
          enddo
         enddo
         call calcuv(ndat,nfilt,mxscn,mxfilt,mstars,
     &               mdate,mbaseline,mhours,lambda,u,v)
         mode=1
         call binavis(ndat,nfilt,mxscn,mxfilt,mdate,
     &                mbaseline,mhours,lambda,dlambda,
     &                modvis2,modphas,u,v,seprad(nd),parad(nd),
     &                mode,ftcurve)
         do ig=1,nfilt
          cchisq(ig)=0.d0
          nchisq(ig)=0
          do iid=1,ndat
           if(viserr(firstid+iid-1,ig).ne.2)then
            cchisq(ig)=cchisq(ig)
     &                +((obsvis2(firstid+iid-1,ig)-modvis2(iid,ig))
     &                /viserr(firstid+iid-1,ig))**2
            nchisq(ig)=nchisq(ig)+1
           endif
          enddo
         enddo
c
         timstep=0.2
         do k=1,mxpl
          mstars(k)=stars(id)
          mdate(k)=date(id)
          mbaseline(k)=baseline(id)
          mhours(k)=dble(timmin+timstep*(k-1))
          do ig=1,nfilt
           lambda(k,ig)=lambda0(id,ig)
           dlambda(k,ig)=dlambda0(id,ig)
          enddo
          if(mhours(k).gt.dble(timmax))goto 40
         enddo
   40    mdat=min(k,mxpl)
         call calcuv(mdat,nfilt,mxscn,mxfilt,mstars,mdate,mbaseline,
     &               mhours,lambda,u,v)
         mode=1
         call binavis(mdat,nfilt,mxscn,mxfilt,mdate,mbaseline,mhours,
     &                lambda,dlambda,modvis2,modphas,u,v,
     &                seprad(nd),parad(nd),mode,ftcurve)
        endif

        DO IG=1,NFILT

         if(ig.eq.1)then
          write(header,43)starno,int(lambda0(id,ig))
   43     format('FK5 ',i4,3x,i3,'nm')
          if(calib)write(header(1:8),'(a)')'calibr. '
         elseif(ig.eq.2)then
          write(header,44)date(id),int(lambda0(id,ig))
   44     format(2x,i6,3x,i3,'nm')
         elseif(ig.eq.3)then
          write(header,45)baseline(id),int(lambda0(id,ig))
   45     format(3x,'Bl ',i2,3x,i3,'nm')
         endif
         call pgsci(5)
         vr=vis2max(ig)-vis2min(ig)
         vis2max(ig)=vis2max(ig)+vr/8.
         call pgenv(timmin,timmax,
     &              vis2min(ig),vis2max(ig),0,1)
         csize=csize+1
         call pgsch(csize)
         call pglabel('','',header)
         csize=csize-1
         call pgsch(csize)
         if(ig.eq.1)call pgsci(2)
         if(ig.eq.2)call pgsci(3)
         if(ig.eq.3)call pgsci(11)
c
         do n=1,ndat
          iid=n+firstid-1
          y(n)=real(obsvis2(iid,ig))
          t(n)=real(hours(iid))
          if(viserr(iid,ig).ne.2)then
           call pgpoint(1,t(n),y(n),1)
           obsvis2lo=y(n)-real(viserr(iid,ig))
           obsvis2hi=y(n)+real(viserr(iid,ig))
           call pgerry(1,t(n),obsvis2lo,obsvis2hi,1.0)
          endif
          y(n)=0
         enddo

         header='                                        '//
     &          '                                        '
         if(nchisq(ig).ne.0)then
          call pgsci(6)
          write(header,46)real(cchisq(ig))/real(nchisq(ig))
   46     format('\gx\u2\d/\gn = ',f6.2)
          call pgmtext('T',1.0,0.5,0.5,header)
         endif
c
c        Plot model visibilities here and now
         if(nd.ne.0)then
          call pgsci(6)
          do k=1,mdat
           y(k)=real(modvis2(k,ig))
           t(k)=real(mhours(k))
          enddo
          call pgline(mdat,t,y)
         endif
c       
        enddo      ! End of ig-loop
c
c       Now prepare for the next date
        do ig=1,nfilt
         vis2max(ig)=0
         vis2min(ig)=0
        enddo
        timmin=real(hours(id+1))
        timmax=timmin
        firstid=lastid+1
       ENDIF
  200 CONTINUE
      call pgend
      call pgask(.true.)
      return
      end