      program bsfit
c
c     Reads delay data from .AVG files and fits baseline and star positions
c     Uses HYCHOL package; not all modes allowed.
c
c     C.A. Hummel, 4 September 1993
c
      integer*4 maxfits,parlines,maxscans,maxstars,maxpars
      integer*4 smaxscans,mxrun,nbadmax
      parameter (maxfits=100)    ! number of nights or blocks
      parameter (parlines=22)
      parameter (maxscans=60000)! total number of scans 
      parameter (maxstars=100)
      parameter (maxpars=maxfits*4+maxstars*2)
      parameter (mxrun=50)
      parameter (nbadmax=maxscans/10)
      parameter (smaxscans=500)
      real*8 pi,rad,lon,lat
      PARAMETER ( PI = 3.141592653589793D0 )
      PARAMETER ( RAD = 180.D0 / PI  )
      PARAMETER ( LON = 118.059167D0 )
      PARAMETER ( LAT =  34.2166944D0)

c for HYCHOL
      integer*4 idim,ir
      parameter (idim=4) ! 4 baseline coordinates
      parameter (ir=idim)! number of right side columns
      real*8 a(idim,idim),b(idim,idim),c(idim,idim) ! sub-matrices
      real*8 bcoeff(idim),scoeff(idim)
      integer*4 ihbcol,ihzei,ihba1,ihbdi
      parameter (ihbcol=(maxstars-1)*2/idim+1)
      parameter (ihzei=maxfits+ihbcol)
      parameter (ihba1=ihzei+1)
      parameter (ihbdi=1)
      integer*4 iadm(ihzei,ihba1) ! address matrix
      real*8 ab(idim,idim,ihbcol),ar(idim,ir)

      integer*4 iou,igl,igr ! file unit specifiers
      integer*4 kkkk,nofil,lpru,ifile(7,2),lfile(3,2)
      common/out/iou,kkkk
      common/ibmf/lpru
      common/defie/nofil,ifile
      character*8 lfnam(2)

c from .AVG file
      integer*4 iyear(maxscans),imon(maxscans),iday(maxscans)
      integer*4 istar(maxscans),ibase,npts(maxscans) 
      integer*4 ifit(maxscans),nnscans(maxfits)
      real*8 hours(maxscans),delay(maxscans),sigma(maxscans)
c from parameter file
      integer*4 years(maxfits),months(maxfits),days(maxfits)
      integer*4 bstars(maxstars),fstars(maxstars)
      integer*4 badscans(nbadmax)
      integer*4 basenos(maxfits),startut(maxfits),stoput(maxfits)
      real*8 bx(maxfits),by(maxfits),bz(maxfits),bc(maxfits)
      real*4 ebx(maxfits),eby(maxfits),ebz(maxfits),ebc(maxfits)
      real*8 ut1utc(maxfits),tdtut1(maxfits),wtmode
c from catalog
      integer*4 spectrum
      real*8 ram,decm
      real*4 par4,rv4,pmra4,pmdec4,mag,angdia,bv
      real*8 par,rv,pmra,pmdec
c
      integer*4 stars(maxstars)
      integer*4 nfits,bfits,nfstars,nbstars
      real*8 stardra(maxstars),starddec(maxstars) ! [sec.] [arc sec.]
      real*8 ra,ha(maxscans),dec(maxscans)
      real*8 delayc(maxscans),dres(maxscans)
      real*8 chisq,weight,sweight
      real*8 jd,djul,tjd,ut1,utsec,raa,deca,gast,sinalt,cosalt
      real*8 cd,ch,sd,sh
      real*8 delra(maxstars),deldec(maxstars)     ! [sec.] [arc sec.]
      real*8 edelra(maxstars),edeldec(maxstars)
      real*4 zd4(maxscans),ha4(maxscans),ma4(maxscans)
      character infile*80,data*80,dev*1,plotfile*13
      character*1024 parline(parlines)
      character*8 name,fname(maxstars)
c
c for running average
      real*4 trun(mxrun),yrun(mxrun)
c     
c for PGPLOT
      integer*4 fknum(maxscans)
      real*4 hours4(maxscans),dres4(maxscans),sigma4(maxscans)
      real*4 shours4(smaxscans),sdres4(smaxscans)
      character header*40
c
c-----Zero arrays-------------------------------------------------------
      do j=1,maxstars
       stars(j)=0
       stardra(j)=0.D0
       starddec(j)=0.D0
       delra(j)=0.D0
       deldec(j)=0.D0
      enddo
      do j=1,maxfits
       nnscans(j)=0
      enddo
c
      open(3,file='bsfit.res',status='new')  ! Open output file
c
c-----Get user input----------------------------------------------------
   10 write(6,'(a,$)')' Please enter parameter file name: '
      read(5,1)infile
    1 format(a80)
      open(1,file=infile,status='old',err=10)
c.... Now read the parameter file in FIXED FORMAT!
c12345678901   --> all input values starting in column 11
cINFILE   =920808.avg        [.AVG file]
cNFITS    =2                 [number of baseline fits]
cUT1-UTC  =0.3,0.3           [seconds]
cTDT-UT1  =55,55             [seconds]
cDAYS     =8,15               
cMONTHS   =8,8             
cYEARS    =1992,1992
cSTARTUT  =0,0               [seconds]
cSTOPUT   =86400,86400       [seconds]
cBASENOS  =1,2               [1=NS,2=ES]
cBX       =                  [initial guess for Bx in microns]
cBY       =                          -"-        By
cBZ       =                          -"-        Bz
cBC       =                          -"-        B. constant
cNBSTARS  =12                [# of stars used for baseline fits]
cBSTARS   =1619,862,79,52,33,19,1568,1534,835,848,66,724
cNFSTARS  =12                [# of stars for which positions will be calc.]
cFSTARS   =1619,862,79,52,33,19,1568,1534,835,848,66,724
cWTMODE   =0                 [0=eq. weights,1=NPTS,2=rms,3=rms/sqrt(NPTS)]
cNBAD     =2                 [number of bad scans]
cBADSCANS =33,54             [scans to be dropped from fit]
cMODE     =1                 [1=RA centering (default)]
c
      do i=1,parlines       
       read(1,2,end=11)parline(i)
    2  format(a1024)
       write(3,1)parline(i)(1:1024)    
      enddo
      goto 12
   11 write(6,*)'EOF or read error on par file or data file!'
      stop
   12 continue
      read(parline(1) (11:1024),1,err=11)infile
      read(parline(2) (11:1024),*,err=11)nfits
      if(nfits.gt.maxfits)then
       write(6,*)'***ERROR: NFITS > ',maxfits,' (=MAXFITS)!'
       stop
      endif
      read(parline(3) (11:1024),*,err=11,end=11)(ut1utc(j),j=1,nfits)
      read(parline(4) (11:1024),*,err=11,end=11)(tdtut1(j),j=1,nfits)
      read(parline(5) (11:1024),*,err=11,end=11)(days(j),j=1,nfits)
      read(parline(6) (11:1024),*,err=11,end=11)(months(j),j=1,nfits)
      read(parline(7) (11:1024),*,err=11,end=11)(years(j),j=1,nfits)
      read(parline(8) (11:1024),*,err=11,end=11)(startut(j),j=1,nfits)
      read(parline(9) (11:1024),*,err=11,end=11)(stoput(j),j=1,nfits)
      read(parline(10)(11:1024),*,err=11,end=11)(basenos(j),j=1,nfits)
      read(parline(11)(11:1024),*,err=11,end=11)(bx(j),j=1,nfits)
      read(parline(12)(11:1024),*,err=11,end=11)(by(j),j=1,nfits)
      read(parline(13)(11:1024),*,err=11,end=11)(bz(j),j=1,nfits)
      read(parline(14)(11:1024),*,err=11,end=11)(bc(j),j=1,nfits)
      read(parline(15)(11:1024),*,err=11)nbstars
      if(nbstars.gt.maxstars)then
       write(6,*)'***ERROR: NBSTARS > ',maxstars,' (=MAXSTARS)!'
       stop
      endif
      read(parline(16)(11:1024),*,err=11,end=11)(bstars(j),j=1,nbstars)
      read(parline(17)(11:1024),*,err=11)nfstars
      if(nfstars.gt.maxstars)then
       write(6,*)'***ERROR: NFSTARS > ',maxstars,' (=MAXSTARS)!'
       stop
      endif
      read(parline(18)(11:1024),*,err=11,end=11)(fstars(j),j=1,nfstars)
      read(parline(19)(11:1024),*,err=11)wtmode
      read(parline(20)(11:1024),*,err=11)nbad
      if(nbad.gt.nbadmax)then
       write(6,*)'***ERROR: NBAD > ',NBADMAX,' (=NBADMAX)!'
       stop
      endif      
      read(parline(21)(11:1024),*,err=11,end=11)(badscans(j),j=1,nbad)
      read(parline(22)(11:1024),*,err=11)mode
      close(1)
      write(6,'(1x,a8,a40)')'INFILE= ',infile
      write(6,*)'Number of fits: ',nfits
      write(6,*)'UT1-UTC= ',(ut1utc(j),j=1,nfits)
      write(6,*)'TDT-UT1= ',(tdtut1(j),j=1,nfits)
      write(6,*)'DAYS= ',(days(j),j=1,nfits)
      write(6,*)'MONTHS= ',(months(j),j=1,nfits)
      write(6,*)'YEARS= ',(years(j),j=1,nfits)
      write(6,*)'STARTUT= ',(startut(j),j=1,nfits)
      write(6,*)'STOPUT= ',(stoput(j),j=1,nfits)
      write(6,*)'BASENOS= ',(basenos(j),j=1,nfits)
      write(6,*)'BX= ',(bx(j),j=1,nfits)
      write(6,*)'BY= ',(by(j),j=1,nfits)
      write(6,*)'BZ= ',(bz(j),j=1,nfits)
      write(6,*)'BC= ',(bc(j),j=1,nfits)
      write(6,*)'Number of stars for baseline fit: ',nbstars
      write(6,*)'BSTARS= ',(bstars(j),j=1,nbstars)
      write(6,*)'Number of stars for position fits: ',nfstars
      write(6,*)'FSTARS= ',(fstars(j),j=1,nfstars)
      write(6,*)'WTMODE= ',wtmode
      write(6,*)'Number of bad scans: ',nbad
      write(6,*)'BADSCANS= ',(badscans(j),j=1,nbad)
      write(6,*)'MODE=',mode
      write(6,*)
c
c     write(6,*)
      write(6,'(a,$)')
     &   ' Read MFIT.POS for initial position offsets? (1/0) '
      read(5,*)icom
      if(icom.eq.1)then
       open(16,file='mfit.pos',status='old',err=13)
       goto 14
   13  write(6,*)' *** ERROR *** File mfit.pos missing!'
       stop
   14  do j=1,maxscans
        read(16,1,end=15)data
        if(index(data,'CHANGE IN STAR POSTION').ne.0)then
         read(16,*)
         read(16,*)
         inum=0
         do i=1,maxstars
          read(16,1,end=15)data
          if(data.eq.'                                    ')goto 14
          if(data(1:80).ne.data(1:40))then
           inum=inum+1
           read(data,4,end=14)STARS(inum),STARDRA(inum),STARDDEC(inum)
    4      format(i5,25x,f8.4,10x,f7.3)
          endif
         enddo
         write(6,*)' *** WARNING *** More stars in list'
         goto 14
        endif
       enddo
c
   15  write(6,*)'Read ',inum,' valid positions from mfit.pos'
       close(16)
      endif
c
c-----Open data file and read data--------------------------------------
      open(1,file=infile,status='old',err=11)
      n=1
      write(6,*)'Reading data...'
      do 20 i=1,maxscans
       READ (1,21,END=23,err=11)IMON(n),IDAY(n),IYEAR(n),HOURS(n),
     $                          ISTAR(n),DELAY(n),ISCAN,IBASE,
     $                          SIGMA(n),NPTS(n)
   21  FORMAT (I2, 1X, I2.2, 1X, I2.2, F12.9, I5, F18.4, 1X, I5,
     $         I2, 1X, F7.2, I5 )
       IYEAR(n) = IYEAR(n) + 1900
       UTSEC = 3600.D0 * HOURS(n)
       sigma4(n)=real(sigma(n))
c......Check wether star is in any list BSTARS or FSTARS
       do j=1,nbstars
        if(istar(n).eq.bstars(j))goto 24
       enddo
       do j=1,nfstars
        if(istar(n).eq.fstars(j))goto 24
       enddo
       WRITE(6,*) ' STAR ',istar(n),' NOT IN FIT; SCAN ',I 
       goto 20
   24  IFIT(n) = 0
c......Determine id of fit and consult selection criteria
       DO 22 J = 1, NFITS
           DO L=1,NBAD
            IF ( I        .EQ. BADSCANS(L) ) GO TO 22
           ENDDO
           IF ( IMON(N)   .NE. MONTHS(J)   ) GO TO 22
           IF ( IDAY(N)   .NE. DAYS(J)     ) GO TO 22
           IF ( IYEAR(N)  .NE. YEARS(J)    ) GO TO 22
           IF ( IBASE     .NE. BASENOS(J)  ) GO TO 22
           IF ( int(UTSEC).LT. STARTUT(J)  ) GO TO 22
           IF ( int(UTSEC).GE. STOPUT(J)   ) GO TO 22
           IFIT(N) = J
           NNSCANS(J)=NNSCANS(J)+1
           N=N+1
           GO TO 20
   22  CONTINUE
       WRITE(6,*) ' SCAN ',I,' NOT IN A BASELINE FIT '
   20 continue
      write(6,*)'***ERROR: Number of scans > ',maxscans, 
     &          ' (=MAXSCANS)!'
      stop
c
   23 nscans=n-1
      write(6,*)'Total number of scans: ',i-1
      write(6,*)'Number of valid scans: ',nscans
      write(3,*)
      write(3,*)'Number of valid scans: ',nscans
      close(1)
c
c.....Check sequence
      n=0
      do j=1,nfits
       do i=1,nnscans(j)
        n=n+1
        if(ifit(n).ne.j)then
         write(6,*)'***ERROR: data sequence error!'
         stop
        endif
       enddo
      enddo
c
c.... Calculate some useful information on weights
      call mdian1(sigma4,nscans,smed)
      write(3,*)'Median delay uncertainty: ',smed
      sweight=0.D0
      do i=1,nscans
       sweight=sweight+1.D0/sigma(i)**2
      enddo
      sweight=sweight/dble(nscans)
      smean=real(1.D0/sqrt(sweight))
      write(3,*)'Mean   delay uncertainty: ',smean
c
c-----Get star positions from catalog and calculate expected delays-----
      write(6,*)'Calculating expected delays...'
      chisq=0.D0
      do i=1,nscans
       CALL CATALOG( ISTAR(i), NAME, RAM, DECM, PMRA4, PMDEC4,
     $               PAR4, RV4, MAG, BV, SPECTRUM, ANGDIA )
       pmra=dble(pmra4)
       pmdec=dble(pmdec4)
       par=dble(par4)
       rv=dble(rv4)
c......Fill in names of stars of input list
       do j=1,nfstars    
        if(istar(i).eq.fstars(j))fname(j)=name
       enddo
c......optionally include initial positional offsets
       do j=1,maxstars
        if(istar(i).eq.stars(j))then
         ram=ram+stardra(j)/3600.D0
         decm=decm+starddec(j)/3600.D0
        endif
       enddo
c......Do astrometry
       DJUL = JD ( IDAY(i), IMON(i), IYEAR(i) )
       UT1  = hours(i)/24.D0 
       TJD  = DJUL + UT1 + tdtut1(ifit(i))/86400.D0
       NBODY = 3
       CALL APSTAR (TJD, NBODY, RAM, DECM, PMRA, PMDEC, PAR,
     $                          RV, RAA, DECA)
       CALL SIDTIM ( DJUL, UT1, 1, GAST )
       CALL DIURN ( TJD, GAST, RAA, DECA, RA, DEC(i) )
       HA(i) = 15.D0 *( GAST - RA ) - LON
       IF ( HA(i) .LT. -180 ) HA(i) = HA(i) + 360.D0
       DELAYC(i) =  BC(IFIT(i))                 
     $            + BX(IFIT(i)) * COS(DEC(i)/RAD) * COS(HA(i)/RAD)
     $            - BY(IFIT(i)) * COS(DEC(i)/RAD) * SIN(HA(i)/RAD)
     $            + BZ(IFIT(i)) * SIN(DEC(i)/RAD)
       DRES(I)=DELAY(I)-DELAYC(I)

       ha4(i)=real(ha(i))

c......Calculate mirror angle and ...
       COSMA = real(SIN(HA(i)/RAD) * COS(DEC(i)/RAD))
       ma4(i) = real(RAD) * ACOS( COSMA )
c......zenith distance
       SINALT  = SIN(LAT/RAD) * SIN(DEC(i)/RAD) +
     $           COS(LAT/RAD) * COS(DEC(i)/RAD) * COS(HA(i)/RAD)
       COSALT  = COS( ASIN ( SINALT ) )
       ZD4(i)  = real(90.D0 - ASIN( SINALT )*RAD)
c
       IF ( WTMODE .EQ. 0 ) THEN
           sigma(i) = dble(smean)
       ELSE IF ( WTMODE .EQ. 1 ) THEN
           sigma(i) = 20.D0*sqrt(1.D0/dble(NPTS(i)))
       ELSE IF ( WTMODE .EQ. 2 ) THEN
           sigma(i)=max(sigma(i),0.5D0) 
       ELSE IF ( WTMODE .EQ. 3 ) THEN
           sigma(i)=max(sigma(i),0.5D0)*sqrt(1.D0/dble(NPTS(i)))
       ELSE
           write(6,*)'*** ERROR: WTMODE',wtmode,' not allowed'
           stop
       END IF
       chisq=chisq+(dres(i)/sigma(i))**2
      enddo
      write(6,*)'Initial Reduced Chisq = ',chisq/dble(nscans)
c
c-----Plot the residuals from the initial guess-------------------------
      ii=0
      xmin=real(hours(1))
      xmax=xmin
      ymin=real(dres(1))
      ymax=ymin
      ifito=ifit(1)
      write(6,'(a,$)')' Screen (s) or hardcopy (h,p,v)? '
      read(5,'(a)')dev
      plotfile='/ega'
      nx=1
      ny=1
      if(dev.eq.'h')plotfile='bsfit1.lpm/lpm'
      if(dev.eq.'p')plotfile='bsfit1.ps/ps'
      if(dev.eq.'v')plotfile='bsfit1.ps/vps'
      if(dev.eq.'p'.or.dev.eq.'v')then
       write(6,'(a,$)')' Number of plot in x and y direction: '
       read(5,*)nx,ny
      endif
      call pgbegin(0,plotfile,nx,ny)
      do i=1,nscans
       ii=ii+1                      
       dres4(ii)=real(dres(i))            
       hours4(ii)=real(hours(i))
       xmin=min(xmin,hours4(ii))
       xmax=max(xmax,hours4(ii))
       ymin=min(ymin,dres4(ii))
       ymax=max(ymax,dres4(ii))
       if(ifit(i+1).ne.ifito.or.i.eq.nscans)then
        call pgsci(5)
        call pgsch(1.5)
        call pgenv(xmin,xmax,ymin,ymax,0,1)
        header='                                        '
        write(header,35)basenos(ifito),years(ifito)-1900,months(ifito),
     &                  days(ifito)
   35   format('Initial guess: baseline ',i2,', on ',i2.2,i2.2,i2.2)
        call pglabel('UT','Residual delay O-C [microns]',header)
        call pgsci(3)
        call pgsch(1.0)
        call pgpoint(ii,hours4,dres4,2)
        ifito=ifit(i+1)
        ii=0
        xmin=real(hours(i+1))
        xmax=xmin
        ymin=real(dres(i+1))
        ymax=ymin
       endif
      enddo
      call pgend
C
c-----Initialize HYCHOL-------------------------------------------------
      write(6,*)'Initialize HYCHOL...'
      iou=6
      igl=17
      igr=18
      ntrsz=32258
      lpru=ntrsz/8
      OPEN(17,ACCESS='DIRECT',RECL=19040)
      OPEN(18,ACCESS='DIRECT',RECL=19040)
      call dinit(2)
      lfnam(1)='LE.SCR'
      lfile(1,1)=igl
      lfile(2,1)=idim*idim
      lfile(3,1)=ihzei*(ihba1-1)
      lfnam(2)='RI.SCR'
      lfile(1,2)=igr
      lfile(2,2)=idim*ir
      lfile(3,2)=ihzei
      call vorhy(iadm,ihzei,ihba1,lfile,lfnam)
c
c-----Now prepare normal equations NE for HYCHOL------------------------
      if(nbstars.eq.0)then
       bfits=0  ! Do not fit baseline
      else
       bfits=nfits
      endif
c     Number of parameters = BFITS* 4 + NFSTARS * 2
      npars=bfits*4+nfstars*2
      write(6,*)'Prepare and store sub-matrices...'
c
c.....Left side: diagonal and border, right side along border
c     There are NFITS symmetric sub-matrices in the diagonal of NE, and
c     NFITS*IHBRA sub-matrices in border band. Each one has to be stored
c     using PUTAR. All these have physical and logical dimension idim,
c     except the last hyper column of border which might have a smaller
c     logical dimension (ia). The border band has NFSTARS*2 columns.
c  
      ihbra=((nfstars-1)*2)/idim+1 ! # of sub-matrices in border row
      ia=idim-mod(nfstars*2,idim)  ! valid number of columns in last one
      maxzei=nfits+ihbra

c     write(6,*)'nfits=',nfits
c     do n=1,nfits
c      write(6,*)'fit ',n,' # scans = ',nnscans(n)
c     enddo
c     write(6,*)'ihzei=',ihzei,', idim=',idim,', ihba1=',ihba1
c     write(6,*)'ihbra=',ihbra,', ia=',ia,', maxzei=',maxzei
c     write(6,*)'ihbdi=',ihbdi,', ihbcol=',ihbcol

      DO N=1,NFITS

       iad=n
       jad=1

       do i=1,idim
        do iir=1,ir
         ar(i,iir)=0.d0  ! initialize right hand sub-matrix
        enddo
        if(n.eq.1)ar(i,i)=1.d0     ! initialize unit vectors in ar
        do j=1,idim
         a(i,j)=0.d0     ! initialize sub-matrix left side
         do k=1,ihbra
          ab(i,j,k)=0.d0 ! initialize ihbra sub-matrices in border
         enddo
        enddo
       enddo
       ar(1,1)=0.d0
c------
       istart=1
       do kk=1,n-1
        istart=istart+nnscans(kk)
       enddo
       do k=1,nnscans(n)  ! nnscans data-records starting at istart
        weight=1.d0/sigma(istart+k-1)**2
        CD = COS(DEC(istart+k-1)/RAD)
        SD = SIN(DEC(istart+k-1)/RAD)
        CH = COS(HA(istart+k-1) /RAD)
        SH = SIN(HA(istart+k-1) /RAD)

        bcoeff(1)= 1.D0    ! baseline coefficients
        bcoeff(2)= cd*ch
        bcoeff(3)=-cd*sh
        bcoeff(4)= sd

        do js=1,nfstars  ! determine no. of current star in array fstars 
         if(istar(istart+k-1).eq.fstars(js))goto 30
        enddo
   30   is=idim-(mod(js,2)*2)-1
        ib=(js-1)/2+1
        do l=1,idim
         scoeff(l)=0.d0
        enddo
        scoeff(is) = 1.D-7 *               
     $    (BZ(IFIT(ISTART+K-1))*CD
     $    +BY(IFIT(ISTART+K-1))*SD*SH
     $    -BX(IFIT(ISTART+K-1))*SD*CH)     ! declination
        scoeff(is+1) = 1.D-7 *             
     $    (BY(IFIT(ISTART+K-1))*CD*CH
     $    +BX(IFIT(ISTART+K-1))*CD*SH)     ! right ascension

        do i=1,idim
         ar(i,1)=ar(i,1)+bcoeff(i)*dres(istart+k-1)*weight ! r.side
         do j=1,idim
          a(i,j)=a(i,j)+bcoeff(i)*bcoeff(j)*weight
          ab(i,j,ib)=ab(i,j,ib)+bcoeff(i)*scoeff(j)*weight
         enddo
        enddo
       enddo
c------
       call putar(iadm,ihzei,iad,jad,a,idim,idim,idim,igl) ! diagonal
       do k=1,ihbra-1
        do i=1,idim
         do j=1,idim
          a(i,j)=ab(i,j,k)
         enddo
        enddo                                              ! border
        call putar(iadm,ihzei,iad,ihba1-1-ihbra+k,a,idim,idim,idim,igl)
       enddo
       do i=1,idim
        do j=1,idim
         a(i,j)=ab(i,j,ihbra)
        enddo
       enddo                                               ! last border
       call putar(iadm,ihzei,iad,ihba1-1,a,idim,idim,ia,igl)
       call putar(iadm,ihzei,iad,ihba1,ar,idim,ir,ir,igr)  ! right hand
      ENDDO
c
c.....Left side: lower right, right side lower right
c     There are IHBRA non-zero sub-matrices in the diagonal.
c     We add one equation, setting the right ascension offset of
c     the last star to zero, i.e. a line with only one non-zero
c     element at the end.

      DO N=1,IHBRA

       iad=ihzei-ihbra+n
       jad=ihba1-1-ihbra+n

       do i=1,idim
        do iir=1,ir
         ar(i,iir)=0.d0 ! initialize right hand sub-matrix
        enddo
        do j=1,idim
         a(i,j)=0.d0    ! initialize sub-matrix left side
        enddo
       enddo
c      do nn=1,ihbra-n  ! store zero matrices to the right of diagonal
c       call putar(iadm,ihzei,iad,jad+nn,a,idim,idim,idim,igl)
c      enddo

       do k=1,nscans
        weight=1.d0/sigma(k)**2
        CD = COS(DEC(k)/RAD)
        SD = SIN(DEC(k)/RAD)
        CH = COS(HA(k) /RAD)
        SH = SIN(HA(k) /RAD)

        do js=1,nfstars  ! determine no. of current star in array fstars 
         if(istar(k).eq.fstars(js))goto 31
        enddo
   31   is=idim-(mod(js,2)*2)-1
        ib=(js-1)/2+1
        if(ib.eq.n)then
         do l=1,idim
          scoeff(l)=0.d0
         enddo
         scoeff(is) = 1.D-7 *
     $     (BZ(IFIT(K))*CD
     $     +BY(IFIT(K))*SD*SH
     $     -BX(IFIT(K))*SD*CH)
         scoeff(is+1) = 1.D-7 *
     $     (BY(IFIT(K))*CD*CH
     $     +BX(IFIT(K))*CD*SH)
        
         do i=1,idim
          ar(i,1)=ar(i,1)+scoeff(i)*dres(k)*weight ! right hand
          do j=1,idim
           a(i,j)=a(i,j)+scoeff(i)*scoeff(j)*weight
          enddo
         enddo
        endif
       enddo

       if(n.eq.ihbra)then
        a(ia,ia)=a(ia,ia)+1.d0 ! right ascension zero-ing
        call putar(iadm,ihzei,iad,jad,a,idim,idim,ia,igl)
        call putar(iadm,ihzei,iad,ihba1,ar,idim,ir,ir,igr)  ! right hand
       else
        call putar(iadm,ihzei,iad,jad,a,idim,idim,idim,igl)
        call putar(iadm,ihzei,iad,ihba1,ar,idim,ir,ir,igr)  ! right hand
       endif

      ENDDO
c
c.....Now call HCHOL
      write(6,*)'Call HCHOL now...'
      call hchol(iadm,ihzei,ihba1,maxzei,ihbdi,ihbra,a,b,c,idim,
     &           ir,igl,igr)
      write(6,*)'HCHOL finished, error code = ',kkkk
c
c-----Get results and update baselines and star position offsets--------
      do n=1,nfits
       call getar(iadm,ihzei,n,ihba1,ar,idim,ir,igr)
       bc(n)=ar(1,1)+bc(n)
       bx(n)=ar(2,1)+bx(n)
       by(n)=ar(3,1)+by(n)
       bz(n)=ar(4,1)+bz(n)
       ebc(n)=0.0 ! HCHOL does not calculate errors, unfortunately
       ebx(n)=0.0
       eby(n)=0.0
       ebz(n)=0.0
      enddo
c
      do n=1,ihbra
       call getar(iadm,ihzei,ihzei-ihbra+n,ihba1,ar,idim,ir,igr)
       j=(n-1)*(idim/2)+1
       deldec(j) = ar(1,1) * RAD * 3600.D0        * 1.D-7
       delra(j)  = ar(2,1) * RAD * 3600.D0 /15.D0 * 1.D-7
       edeldec(j)= 0.001d0
       edelra(j) = 0.001d0
       if(n.eq.ihbra.and.ia.ne.4)goto 32
       deldec(j+1) = ar(3,1) * RAD * 3600.D0        * 1.D-7
       delra(j+1)  = ar(4,1) * RAD * 3600.D0 /15.D0 * 1.D-7
       edeldec(j+1)= 0.001d0
       edelra(j+1) = 0.001d0
      enddo
   32 write(6,*)'Solution arrays updated.'
c
c-----Calculate residuals and plot them---------------------------------
c.... Plot residuals versus time
      xmin=real(hours(1))
      xmax=xmin
      ymax=0.
      ymin=0.
      chisq=0.D0
      do i=1,nscans
       jj=1
       do j=1,nfstars
        if(istar(i).eq.fstars(j))jj=j
       enddo
       CD = COS((DEC(i)+deldec(jj)/3600.D0)/RAD)
       SD = SIN((DEC(i)+deldec(jj)/3600.D0)/RAD)
       CH = COS((HA(i)-delra(jj)/3600.D0*15.D0)/RAD)
       SH = SIN((HA(i)-delra(jj)/3600.D0*15.D0)/RAD)
       DELAYC(i) = BC(IFIT(i))                 
     $           + BX(IFIT(i)) * CD * CH
     $           - BY(IFIT(i)) * CD * SH
     $           + BZ(IFIT(i)) * SD
       chisq=chisq+((delay(i)-delayc(i))/sigma(i))**2
       dres4(i)=real(delay(i)-delayc(i))
       hours4(i)=real(hours(i))
c      fknum(ii)=istar(i)
c      dres4(ii)=real(delay(i)-delayc(i))             
       xmin=min(xmin,hours4(i))
       xmax=max(xmax,hours4(i))
       ymin=min(ymin,dres4(i))
       ymax=max(ymax,dres4(i))
      enddo
      ymin=-5.
      ymax=5.
      write(6,'(a,$)')' Screen (s) or hardcopy (h,p,v)? '
      read(5,'(a)')dev
      plotfile='/ega'
      ny=1
      if(dev.eq.'h')plotfile='bsfit2.lpm/lpm'
      if(dev.eq.'p')plotfile='bsfit2.ps/ps'
      if(dev.eq.'v')then
       plotfile='bsfit2.ps/vps'
       ny=2
      endif
      call pgbegin(0,plotfile,1,ny)
      call pgsci(5)
      call pgsch(1.5)
      call pgenv(xmin,xmax,ymin,ymax,0,1)
      header='                                        '
      write(header,36)real(chisq)/real(nscans)
   36 format('Res. delays (all baselines), Chisq=',f4.1)
      call pglabel('UT','Residual delay O-C [microns]',header)
      call pgsci(3)
      call pgsch(1.0)
      call pgpoint(nscans,hours4,dres4,2)
      call pgend
c
c.... Plot residuals versus zenith angle
      write(6,'(a,$)')' Screen (s) or hardcopy (h,p,v)? '
      read(5,'(a)')dev
      plotfile='/ega'
      ny=1
      if(dev.eq.'h')plotfile='bsfit3.lpm/lpm'
      if(dev.eq.'p')plotfile='bsfit3.ps/ps'
      if(dev.eq.'v')then
       plotfile='bsfit3.ps/vps'
       ny=2
      endif
      call pgbegin(0,plotfile,1,ny)
      call pgsci(5)
      call pgsch(1.5)
      xmin=0.
      xmax=60.
      call pgenv(xmin,xmax,ymin,ymax,0,1)
      header='                                        '
      call pglabel('Zenith angle (degrees)',
     &             'Residual delay O-C [microns]',' ')
      call pgsci(3)
      call pgsch(1.0)
      call pgpoint(nscans,zd4,dres4,2)
      call pgend
c
c.... Plot residuals versus hour angle
      write(6,'(a,$)')' Screen (s) or hardcopy (h,p,v)? '
      read(5,'(a)')dev
      plotfile='/ega'
      ny=1
      if(dev.eq.'h')plotfile='bsfit4.lpm/lpm'
      if(dev.eq.'p')plotfile='bsfit4.ps/ps'
      if(dev.eq.'v')then
       plotfile='bsfit4.ps/vps'
       ny=2
      endif
      call pgbegin(0,plotfile,1,ny)
      call pgsci(5)
      call pgsch(1.5)
      xmin=0.
      xmax=0.
      do i=1,nscans
       xmin=min(xmin,ha4(i))
       xmax=max(xmax,ha4(i))
      enddo
      call pgenv(xmin,xmax,ymin,ymax,0,1)
      header='                                        '
      call pglabel('Hour angle (degrees)',
     &             'Residual delay O-C [microns]',' ')
      call pgsci(3)
      call pgsch(1.0)
      call pgpoint(nscans,ha4,dres4,2)
      call pgend
c
c.... Plot residuals versus mirror angle
      write(6,'(a,$)')' Screen (s) or hardcopy (h,p,v)? '
      read(5,'(a)')dev
      plotfile='/ega'
      ny=1
      if(dev.eq.'h')plotfile='bsfit5.lpm/lpm'
      if(dev.eq.'p')plotfile='bsfit5.ps/ps'
      if(dev.eq.'v')then
       plotfile='bsfit5.ps/vps'
       ny=2
      endif
      call pgbegin(0,plotfile,1,ny)
      call pgsci(5)
      call pgsch(1.5)
      xmin=90.
      xmax=90.
      do i=1,nscans
       xmin=min(xmin,ma4(i))
       xmax=max(xmax,ma4(i))
      enddo
      call pgenv(xmin,xmax,ymin,ymax,0,1)
      header='                                        '
      call pglabel('Mirror angle (degrees)',
     &             'Residual delay O-C [microns]',' ')
      call pgsci(3)
      call pgsch(1.0)
      call pgpoint(nscans,ma4,dres4,2)
      call pgend
c
c-----Write output (and calculate mean and RMS of residuals)------------
c
      WRITE(3,2010)
 2010 FORMAT(/1X,'FIT BASELINE in microns ', /
     $       12X,'CON',17X,' BX', 17X, ' BY', 17X, ' BZ' )
c
      DO N=1,NFITS
         WRITE(3,2020) 
     $         BC(N), EBC(N), BX(N), EBX(N),
     $         BY(N), EBY(N), BZ(N), EBZ(N)
      ENDDO
 2020 FORMAT( 4(F12.2, '(',F5.2, ') ' ) )
      open(4,file='base.dat')
      write(4,2030)'BX       =',(bx(j),j=1,nfits)
      write(4,2030)'BY       =',(by(j),j=1,nfits)
      write(4,2030)'BZ       =',(bz(j),j=1,nfits)
      write(4,2030)'BC       =',(bc(j),j=1,nfits)
 2030 format(a10,1000(f12.2,1x))
      close(4)

      write(3,*)
      write(3,*)'Final Reduced Chisq = ',chisq/dble(nscans)

      write(3,1000)
 1000 FORMAT ( /3X, ' Summary of residuals  equal weight per point ' /
     $  / 13X, '       RESID DELAY      CHANGE IN STAR POSTION     '
     $  / 13X, ' NO.    (micron)   RA(sec of time)    DEC(arc sec) '
     $  / 13X, ' SCAN  mean   rms   delta   sigma     delta  sigma' )
c
c.....Calculate RMS and running averages
      chirun=0.
      timmin=hours4(1)
      timmax=timmin
      timstep=0.25
      do k=1,nfstars
       davg=0.
       drms=0.
       snpts=0.
       do i=1,nscans
        if(istar(i).eq.fstars(k))then
         davg=davg+dres4(i)
         drms=drms+dres4(i)**2
         snpts=snpts+1.
         timmin=min(timmin,hours4(i))
         timmax=max(timmax,hours4(i))
        endif
       enddo
       davg=davg/snpts
       drms=sqrt(drms/snpts-davg**2)
c
c      do nfit=1,nfits
c       do n=1,mxrun
c        trun(n)=timmin+timstep*(n-1)
c        yrun(n)=0.
c        swgn=0.
c        do i=1,nscans
c         if(abs(hours4(i)-trun(n)).lt.3)then
c          if(istar(i).eq.fstars(k).and.ifit(i).eq.nfit)then
c           wg=exp(-(hours4(i)-trun(n))**2)            ! Gaussian
c           wn=1./real(sigma(i)*sigma(i))              ! Noise
c           yrun(n)=yrun(n)+dres4(i)*wg*wn
c           swgn=swgn+wg*wn
c          endif
c         endif
c        enddo
c        if(swgn.ne.0)then
c         yrun(n)=yrun(n)/swgn
c        else
c         yrun(n)=0.
c        endif
c        if(trun(n).gt.timmax)goto 2040
c       enddo
c2040   nrun=min(n,mxrun)
c       do i=1,nscans
c        if(istar(i).eq.fstars(k).and.ifit(i).eq.nfit)then
c         j=int((hours4(i)-timmin)/timstep)+1
c         chirun=chirun+((dres4(i)-yrun(j))/real(sigma(i)))**2
c        endif
c       enddo
c      enddo
c
       hdelra=real(delra(k))
       hdeldec=real(deldec(k))
       do j=1,maxstars               ! account for initial offsets
        if(stars(j).eq.fstars(k))then
         hdelra =real(delra(k)  + stardra(j))
         hdeldec=real(deldec(k) + starddec(j))
        endif
       enddo
       WRITE(3,1100) FSTARS(K), FNAME(K), int(SNPTS), DAVG, DRMS,
     $               HDELRA, EDELRA(k),  HDELDEC, EDELDEC(k)
 1100  FORMAT ( I5, 1X, A8, I4, 2F6.1, 2F8.4, 2X,  2F7.3 )
      enddo 
c     chirun=chirun/real(nscans)
c     write(6,*)'Reduced Chisq relative to running average: ',chirun
c
c-----Recalculate residuals with new baseline and plot------------------
      write(6,'(a,$)')' Screen (s) or hardcopy (h,p,v) or exit (q) ? '
      read(5,'(a)')dev
      plotfile='/ega'
      nx=1
      ny=1
      if(dev.eq.'h')plotfile='bsfit6.lpm/lpm'
      if(dev.eq.'p')plotfile='bsfit6.ps/ps'
      if(dev.eq.'v')plotfile='bsfit6.ps/vps'
      if(dev.eq.'p'.or.dev.eq.'v')then
       write(6,'(a,$)')' Number of plot in x and y direction: '
       read(5,*)nx,ny
      endif
      if(dev.eq.'q'.or.dev.eq.'n')goto 900
   40 write(6,'(a,$)')' Which star should be plotted ? '
      read(5,*)ifknum
      if(ifknum.eq.0)goto 900
      call pgbegin(0,plotfile,nx,ny)
      ii=0
      jj=0
      xmin=real(hours(1))
      xmax=xmin
      ymin=0.
      ymax=0.
      ifito=ifit(1)
      do i=1,nscans
       CD = COS(DEC(i)/RAD)
       SD = SIN(DEC(i)/RAD)
       CH = COS(HA(i) /RAD)
       SH = SIN(HA(i) /RAD)
       DELAYC(i) = BC(IFIT(i))                 
     $           + BX(IFIT(i)) * CD * CH
     $           - BY(IFIT(i)) * CD * SH
     $           + BZ(IFIT(i)) * SD

       ii=ii+1                      
       fknum(ii)=istar(i)
       dres4(ii)=real(delay(i)-delayc(i))             
       if(istar(i).eq.ifknum)then
        jj=jj+1
        if(jj.gt.smaxscans)then
         write(6,*)'***ERROR: maximum number of scans/night exceeded!'
         stop
        endif
        do j=1,nfstars                           
         if(istar(i).eq.fstars(j))then
          sdres4(jj) = real(
     $      delra(j)/rad/3600.D0*15.D0 
     $            * (BY(IFIT(I))*CD*CH+BX(IFIT(I))*CD*SH)
     $     +deldec(j)/rad/3600.D0 
     $            * (BZ(IFIT(I))*CD+BY(IFIT(I))*SD*SH-BX(IFIT(I))*SD*CH)
     $                     )
          shours4(jj)=real(hours(i))
         endif
        enddo
       endif
       hours4(ii)=real(hours(i))
       sigma4(ii)=real(sigma(i))
       xmin=min(xmin,hours4(ii))
       xmax=max(xmax,hours4(ii))
       ymin=min(ymin,dres4(ii))
       ymax=max(ymax,dres4(ii))
       if(ifit(i+1).ne.ifito.or.i.eq.nscans)then
        call pgsci(5)
        call pgsch(1.5)
        call pgenv(xmin,xmax,ymin,ymax,0,1)
        header='                                        '
        write(header,45)basenos(ifito),years(ifito)-1900,months(ifito),
     &                  days(ifito)
   45   format('Fitted baseline ',i2,', on ',i2.2,i2.2,i2.2)
        call pglabel('UT','Residual delay O-C [microns]',header)
        call pgsci(3)
        call pgpoint(ii,hours4,dres4,2)
c...... Now plot the expected residuals for the positional offsets
        call pgsci(1)
        do j=1,ii                               
         if(fknum(j).eq.ifknum)then
          call pgpoint(1,hours4(j),dres4(j),2)
          ylo=dres4(j)-sigma4(j)
          yhi=dres4(j)+sigma4(j)
          call pgerry(1,hours4(j),ylo,yhi,1.0)
         endif
        enddo
        call pgsci(2)
        call pgline(jj,shours4,sdres4)
c        
c ..... Do a running average
c       timstep=0.25
c       do k=1,mxrun
c        trun(k)=shours4(1)+timstep*(k-1)
c        yrun(k)=0.
c        swgn=0.
c        do n=1,ii
c         if(fknum(n).eq.ifknum)then
c          wg=exp(-(hours4(n)-trun(k))**2)       ! Gaussian
c          wn=1./(sigma4(n)*sigma4(n))           ! Noise
c          yrun(k)=yrun(k)+dres4(n)*wg*wn
c          swgn=swgn+wg*wn
c         endif
c        enddo
c        if(swgn.ne.0)then
c         yrun(k)=yrun(k)/swgn
c        else
c         yrun(k)=0.
c        endif
c        if(trun(k).gt.shours4(jj))goto 50
c       enddo
c  50   nrun=min(k,mxrun)
c       call pgsls(2)
c       call pgsci(8)
c       call pgline(nrun,trun,yrun)
c       call pgsls(1)
c
c.......Prepare for next plot
        ifito=ifit(i+1)
        ii=0
        jj=0
        xmin=real(hours(i+1))
        xmax=xmin
        ymin=real(delay(i+1)-delayc(i+1))
        ymax=ymin
       endif
      enddo
c
      call pgend
      goto 40
c
  900 close(3)
      end


