program binafit
c
c Version for HP Fortran on FORNAX (4-Oct-1995).
c
c This program fits the system and orbit parameters of a binary star 
c directly to the visibilities measured during several nights - preferably 
c providing a good coverage of the orbit. It takes a data file containing 
c all .CAL appended to each other. It also reads in a model file with the 
c parameters (or initial guesses). The Levenberg-Marquardt algorithm 
c (Press et al.) is used for the non-linear least squares fitting routine. 
c Note: currently, visibilities larger than 1.2 or less than 0 are flagged:
c a flagged visibility has value=2.D0 and error=2.D0; the same happens to 
c visibilities with zero error (see READVIS).
c
c     Christian A. Hummel      
c
c These are the sign convention for the BINAFIT and BINFIT software:
c (0 deg up, 90 deg left, 180 deg down, 270 deg right for ref. dirs.)
c    +DEC       +RA          -DEC             -RA
c
c In the u-v plots, the position angle increases from the 270 deg dir.
c Retrograd means: the position angle decreases with time.
c In retrograd motion, the star orbits clockwise (indicated by arrow).
c The following statements are preliminary: 
c With the period always positive, the inclination is > 90 for retrograd
c orbits, and < 90 for prograde orbits. 
c
c Main program and 'non-basic'-subroutines use the only two main
c common blocks /data/ and /model/ (via data.inc and model.inc).
c 'Basic'-subroutines get their parameters from the subroutine call
c
      save modvis2
c
      include 'constnts.inc'
      include 'data.inc'
      include 'model.inc'

c ig=filter number(1-3),id=data number, nrec=total number of visibilities
      integer*4 id,ig,nrec

      real*8 a(mxparm),mseprad,mparad
      real*8 modvis2(mxrec,mxfilt),modphas(mxrec,mxfilt)
      real*8 chisq,chisqo,cchisq(mxfilt),mcjd,jd,tjd

c Miscellaneous for RTFIT
      integer*4 nangles
      parameter (nangles=36)  ! number of points on error ellipse
      real*8 ex(nangles),ey(nangles)
      common/edata/ex,ey      ! used only in errell and fixell

c Miscellaneous for restarting MTCARLO
      common/random/iran

      character data*80,prognm*8,versn*3,hist*80,command*1
      character device*4
      character*3 months(12)
      logical modplot,modfit,query,adjust,sorry
      logical mcinit,mcrun
      data months/'Jan','Feb','Mar','Apr','May','Jun',
     &            'Jul','Aug','Sep','Oct','Nov','Dec'/
      data iran /0/

      mcinit=.true.
      modfit=.false.
      modplot=.false.
      starno=0
      nrec=0
      mode=0     ! 0 = Gaussian model in FCHISQ
c                ! 1 = double exponential
c                ! 2 = Lorentzian

      prognm=' BINAFIT'
      versn='4.8'      !  2 November 1994
      write(outc,*)
      write(outc,'(a,a,a,a)')
     &     'This is program ',prognm,', Version ',versn
      write(outc,*)
c Open output file 3:
      open(outdat,file='binafit.out',status='new',err=40)
      goto 60
   40 if(query('Output file exists! Overwrite ? (y/n) '))goto 50
      write(outc,*)'Enter file name! '
      read(inc,'(a)')data
      open(outdat,file=data)
      goto 60
   50 open(outdat,file='binafit.out')

c Write history record
   60 call mkhist(prognm,versn,hist)
      write(outdat,'(a,a70)')'!',hist
c-----------------------------------------------------------------------
c     Now print a selection of options:
c
   30 write(outc,*)
      write(outc,'(a)')'Enter your selection:'
      write(outc,'(a)')
      write(outc,'(a)')'-1- Read data'
      write(outc,'(a)')'-2- Read model'
      write(outc,'(a)')'-3- Plot orbit'
      write(outc,'(a)')'-4- Plot data'
      write(outc,'(a)')'-5- Fit model to data'
      write(outc,'(a)')'-6- Plot statistics of data'
      write(outc,'(a)')'-7- Fit and/or plot r-t values'
      write(outc,'(a)')'-8- Plot predicted visibilities'
      write(outc,'(a)')'-9- Initialize MtC / create data'
      write(outc,'(a)')'-m- Do Monte-Carlo simulations'
      write(outc,'(a)')'-p- Plot Monte-Carlo results'
      write(outc,'(a)')'-e- Calculate formal errors'
      write(outc,'(a)')'-u- Vis vs baseline length'
      write(outc,'(a)')'-r- RESET and read data'
      write(outc,'(a)')'-w- Write data'
      write(outc,'(a)')'-q- Quit'
      write(outc,*)
      write(outc,'(a,$)')'>'
      read(inc,9)command
    9 format(a)
c   1 write(outc,*)''
    1 continue
c
      if(command.eq.'r')then
       modfit=.false.
       modplot=.false.
       starno=0
       nrec=0
       command='1'
       goto 1
      endif
      if(command.eq.'q')goto 300
      if(command.eq.'1')then
       mcrun=.false.
       call readvis
       call calcuv(ndata,nfilt,mxrec,mxfilt,stars,date,baseline,hours,
     &             lambda0,udata,vdata)
       nrec=nfilt*ndata
      elseif(command.eq.'d')then
       call dspot
      elseif(command.eq.'w')then
       call writed(modvis2,starno)
      elseif(command.eq.'2')then
       call readmod
       modplot=.true.
       a(1)=diam(1,1)          ! save model parameters
       a(2)=diam(2,1)
       a(3)=dratio(1,2)
       a(4)=dratio(1,3)
       a(5)=dratio(2,2)
       a(6)=dratio(2,3)
       a(7)=delmag(2,1)
       a(8)=delmag(2,2)
       a(9)=delmag(2,3)
       a(10)=a_true
       a(11)=e_true
       a(12)=period
       a(13)=epoch
       a(14)=inclin
       a(15)=asc_node
       a(16)=arg_peri
      elseif(command.eq.'3')then
       if(.not.modplot)then
        write(outc,'(a)')'*** Sorry: you have to read the model first!'
       else
        if(.not.calib)call plotorb(chisq)
       endif
      elseif(command.eq.'4')then
       if(nrec.eq.0)then
        write(outc,'(a)')'*** Sorry: you have to read the data first!'
       else
        imode=0
        call plotvis(modplot,modvis2,modphas,imode,device)
       endif
      elseif(command.eq.'u')then
       if(nrec.eq.0)then
        write(outc,'(a)')'*** Sorry: you have to read the data first!'
       else
        call uvdist(modplot)
       endif
      elseif(command.eq.'5')then
       if(.not.sorry(nrec,modplot,outc))then
        call mrqfit(modvis2,modphas,adjust,mcrun)
        a(1)=diam(1,1)          ! save model parameters
        a(2)=diam(2,1)
        a(3)=dratio(1,2)
        a(4)=dratio(1,3)
        a(5)=dratio(2,2)
        a(6)=dratio(2,3)
        a(7)=delmag(2,1)
        a(8)=delmag(2,2)
        a(9)=delmag(2,3)
        a(10)=a_true
        a(11)=e_true
        a(12)=period
        a(13)=epoch
        a(14)=inclin
        a(15)=asc_node
        a(16)=arg_peri
       endif
      elseif(command.eq.'6')then
       if(.not.sorry(nrec,modplot,outc))then
        call errstat(modvis2,chisq)
       endif
      elseif(command.eq.'7')then
       if(.not.sorry(nrec,modplot,outc))then
        if(.not.calib)call rtfit
       endif
      elseif(command.eq.'8')then
       if(.not.modplot)then
        write(outc,'(a)')'*** Sorry: you have to read the model first!'
       else
        if(starno.eq.0)then
         write(outc,'(a,$)')'Please enter FK5 number of star: '
         read(inc,'(a)')starno
        endif
        if(.not.calib)then
         call plotmod
        endif
       endif
      elseif(command.eq.'9')then
       if(.not.sorry(nrec,modplot,outc))then
        if(query('Initialize Monte Carlo Method ? (y/n) '))then
         mcinit=.true.
        else
         if(mcinit)then
          write(outc,'(a)')'*** Sorry, no initialization done so far'
          goto 30
         else
          mcinit=.false.
         endif
        endif
        call mtcarlo(modvis2,modphas,mcinit,mcrun)
        chisq=0.D0
        do ig=1,nfilt
         cchisq(ig)=0.D0
         do id=1,ndata
          call fchisq(obsvis2(id,ig),viserr(id,ig),1,1,
     &                modvis2(id,ig),chisqo,mode)
          chisq=chisq+chisqo
          cchisq(ig)=cchisq(ig)+chisqo
         enddo
        enddo
        write(outc,*)
        write(outc,'(a,f7.1)')' Total chisq = ',chisq
        write(outc,'(a,3(f7.1))')' Channel chisq = ',cchisq
       endif
      elseif(command.eq.'m')then
       if(mcinit)then
        write(outc,'(a,a)')'*** Sorry: you have to initialize the',
     &                  ' Monte-Carlo method first!'
       else
        mcrun=.true.
        do i=1,maxitc            ! Clear MONTE CARLO results array
         do j=1,mxparm
          mca(i,j)=0.D0
         enddo
        enddo
   10   call mkhist(prognm,versn,hist)
        write(outc,'(a,a70)')' Time check: ',hist
        write(outc,'(a,a)')'Enter date and time by which program ',
     &               'should end '
        write(outc,'(a,a)')'(month day year hour) e.g. 12 16 1960 19 ! '
        read(inc,*,end=10,err=10)mcmonth,mcday,mcyear,mchour
        if(mcmonth.gt.12)goto 10
        if(mcday.gt.31)goto 10
        if(mchour.gt.24)goto 10
        if(abs(mcyear-1992).gt.20)goto 10
        mcjd=jd(mcday,mcmonth,mcyear)+dble(mchour)/24.D0
        do mci=1,maxitc
         diam(1,1)=a(1)          ! restore old model parameters
         diam(2,1)=a(2)
         dratio(1,2)=a(3)
         dratio(1,3)=a(4)
         dratio(2,2)=a(5)
         dratio(2,3)=a(6)
         delmag(2,1)=a(7)
         delmag(2,2)=a(8)
         delmag(2,3)=a(9)
         a_true=a(10)
         e_true=a(11)
         period=a(12)
         epoch=a(13)
         inclin=a(14)
         asc_node=a(15)
         arg_peri=a(16)
         call cadjmod
         imode=0
         call binavis(ndata,nfilt,mxrec,mxfilt,date,baseline,hours,
     &                lambda0,dlambda0,modvis2,modphas,udata,vdata,
     &                mseprad,mparad,imode,ftcurve)
 
         call mtcarlo(modvis2,modphas,mcinit,mcrun) 
         call mrqfit(modvis2,modphas,adjust,mcrun) 
 
         mca(mci,1)=diam(1,1)        ! save results from iteration mci
         mca(mci,2)=diam(2,1)
         mca(mci,3)=dratio(1,2)
         mca(mci,4)=dratio(1,3)
         mca(mci,5)=dratio(2,2)
         mca(mci,6)=dratio(2,3)
         mca(mci,7)=delmag(2,1)
         mca(mci,8)=delmag(2,2)
         mca(mci,9)=delmag(2,3)
         mca(mci,10)=a_true
         mca(mci,11)=e_true
         mca(mci,12)=period
         mca(mci,13)=epoch
         mca(mci,14)=inclin
         mca(mci,15)=asc_node
         mca(mci,16)=arg_peri
c ...... Check whether current date and time passed the deadline
         call mkhist(prognm,versn,hist)
         do i=1,12
          j=index(hist,months(i))
          if(j.ne.0)then
           read(hist(j-3:j-2),*)iday
           read(hist(j+4:j+7),*)iyear
           read(hist(j+9:j+10),*)ihour
           read(hist(j+12:j+13),*)imin
           imonth=i
          goto 11
          endif
         enddo
   11    write(outc,13)mci,imonth,iday,iyear,ihour,imin,
     &                    mcmonth,mcday,mcyear,mchour     
   13    format(' MONTE CARLO iteration ',i3,' completed',/,
     &   ' Time: ',i2,'/',i2,'/',i4,1x,i2,':',i2.2,', will end: ',
     &             i2,'/',i2,'/',i4,1x,i2,':00')
         tjd=jd(iday,imonth,iyear)+dble(ihour)/24.D0+dble(imin)/3600.D0
         if(tjd.gt.mcjd)goto 14
        enddo
   14   write(outdat,12)hist(38:57)
   12   format('! MONTE CARLO iterations stop ',a20)
        diam(1,1)=a(1)          ! restore old model parameters
        diam(2,1)=a(2)
        dratio(1,2)=a(3)
        dratio(1,3)=a(4)
        dratio(2,2)=a(5)
        dratio(2,3)=a(6)
        delmag(2,1)=a(7)
        delmag(2,2)=a(8)
        delmag(2,3)=a(9)
        a_true=a(10)
        e_true=a(11)
        period=a(12)
        epoch=a(13)
        inclin=a(14)
        asc_node=a(15)
        arg_peri=a(16)
        call cadjmod
       endif
      elseif(command.eq.'p')then
       call mcplot
      elseif(command.eq.'e')then
       if(.not.sorry(nrec,modplot,outc))then
        write(outc,'(a,a)')'This takes a while...'
        call calcerr(modvis2)
       endif
      endif                ! End command interpreter
c-----------------------------------------------------------------------
c     Here we calculate the model visibilities,
c     once data and model have been read and whenever needed.
c
      if(command.eq.'1'.or.command.eq.'2'.or.command.eq.'5')then
      if(nrec.ne.0.and.modplot)then

c      Adjust epoch to nearest cycle (only binaries)
       if(.not.calib)then
        IDAY  = MOD(DATE(1),       100)
        IMON  = MOD(DATE(1)/100,   100)
        IYEAR = MOD(DATE(1)/10000, 100) + 1900
        TJD   = JD ( IDAY, IMON, IYEAR ) - 2440000.D0
        icycles=int((tjd-epoch)/period)
        epoch=epoch+dble(icycles)*period
        a(13)=epoch
       endif

       if(command.ne.'5')then
        adjust=.false.
        do i=1,mxdates
         do ig=1,mxfilt
          vis2rat(i,ig)=0.d0
         enddo
        enddo
        if(query('Adjust calibration constants ? (y/n) '))
     &            adjust=.true.
       endif
       if(.not.calib)then
        imode=0
        call binavis(ndata,nfilt,mxrec,mxfilt,date,baseline,hours,
     &               lambda0,dlambda0,modvis2,modphas,udata,vdata,
     &               mseprad,mparad,imode,ftcurve)
       else
        call calvis(ndata,nfilt,mxrec,mxfilt,date,baseline,hours,
     &              lambda0,modvis2,udata,vdata,
     &              stars,calstars,caldiams)
       endif
       if(adjust)then
        call vadjust(modvis2)
       endif
       chisq=0.D0
       do ig=1,nfilt
        cchisq(ig)=0.D0
        do id=1,ndata
         call fchisq(obsvis2(id,ig),viserr(id,ig),1,1,
     &               modvis2(id,ig),chisqo,mode)
         chisq=chisq+chisqo
         cchisq(ig)=cchisq(ig)+chisqo
        enddo
       enddo
       write(outdat,'(a,f7.1)')'! Total chisq = ',chisq
       write(outdat,'(a,3(f7.1))')'! Channel chisq = ',cchisq
       write(outc,*)
       write(outc,'(a,f7.1)')' Total chisq = ',chisq
       write(outc,'(a,3(f7.1))')' Channel chisq = ',cchisq
      endif
      endif
c-----------------------------------------------------------------------
c     Return to selection menu
      goto 30
c
c     Finish up
c
  300 write(outc,'(a,a)')'Your output is in file: binafit.out.'
      if(iran.ne.0)write(outdat,303)iran
  303 format('! Simulations: ',i7,' random numbers created by RAN1.')
      write(outdat,301)
  301 format('! BINAFIT ends.')
  302 format(a80)
      data='                                        '//
     &     '                                        '
c     Overwrite stuff that is left from the old file
      do i=1,2
       write(outdat,302)data 
      enddo
      close(outdat)
      end