program xlk

c       -- based on the famous lk program - mk


      implicit real*8 (a-h,o-z)
      
      real*8 rct(6000),rphase(6000),resid1(6000),resid2(6000),
     +     xindex(6000),dayofyr(6000), error(6000)
      
      real*8 manbase, sub

      parameter (manbase=39126.0d0)

      character*1 an
      character infile*32, labelx*80,labely*80, term*10, arg*32
      
      data infile/'resid2.tmp'/

      logical dump, first
      
      data dump/.false./
      data first/.true./

      write(*,*) 'Use: xlk [-f<file>] [-ps] [-o]'
      
      term='/xw'
      do i=1, iargc()
         call getarg(i,arg)
         if (arg(1:3).eq.'-ps') term='/ps'
         if (arg(1:2).eq.'-o') dump=.true.
         if (arg(1:2).eq.'-f') infile=arg(3:)
      enddo
      
      open(32,file=infile,form='unformatted',status='old')

      if (dump) then
         open(36,file='res.out',form='formatted',status='unknown')
      endif
      

 1    continue
      rewind 32
      ctmin=1.d30
      ctmax=0.

      do  i=1,5999
         read(32,end=4) ct,dt2,dt2sec,phase,freq,w,terr,y
c         ct=ct+manbase ! commented out since tempo 11 uses mjd already
         rct(i)=ct
         if(ct.lt.ctmin) ctmin=ct
         if(ct.gt.ctmax) ctmax=ct
         p0=dt2sec/dt2
         rphase(i)=phase
         resid1(i)=y*1.0d6*p0            ! prefit
         resid2(i)=1.0d6*dt2sec          ! postfit
         error(i)=terr
         xindex(i)=i
         dayofyr(i)=mod(ct,365.25d0)
         if (dump) write(36,*) rct(i), resid2(i), error(i), 
     +        resid1(i),phase
      enddo
      

 4    continue


      if (dump) then
         close(36)         
         stop 'TOAs dumped to res.out'
      endif

      npts=i-1

      sub=int(ctmin)
      do i=1, npts
         rct(i)=rct(i)-sub
      enddo
      ctmin=ctmin-sub
      ctmax=ctmax-sub
      
 10   continue
      xnpts=npts

      r1=0.
      r2=0.
      do i=1,npts
         r1=dmax1(r1,dabs(resid1(i)))
         r2=dmax1(r2,dabs(resid2(i)))
      enddo

      call pgbegin(0,term,1,1)
      call pgpaper(7.0,0.666)
      call pgask(.false.)

 12   continue
      
      if (.not.first) then
         print*,'Enter command: '
         read(*,1014) an
 1014    format(a1)
      else
         an = '1'
         first = .false.
      endif

      if ((an.eq.'x').or.(an.eq.'0')) goto 100

      if(an.eq.'1') then
         labely='Pre-fit residuals [usec]'
         write(labelx,99) sub
 99      format('MJD [Epoch - ',f7.1,']')
         call plt(rct,ctmin,ctmax,resid1,-r1,r1,npts,rct,p0,
     +            error,labelx,labely)
      endif

      if(an.eq.'3') then
         labely='Pre-fit residuals [usec]'
         labelx='Orbital phase'
         call plt(rphase,0.d0,1.d0,resid1,-r1,r1,npts,rct,p0,
     +            error,labelx,labely)
      endif

      if(an.eq.'2') then
         labely='Post-fit residuals [usec]'
         write(labelx,99) sub
         call plt(rct,ctmin,ctmax,resid2,-r2,r2,npts,rct,p0,
     +            error,labelx,labely)
      endif

      if(an.eq.'4') then
         labely='Post-fit residuals [usec]'
         labelx='Orbital phase'
         call plt(rphase,0.d0,1.d0,resid2,-r2,r2,npts,rct,p0,
     +            error,labelx,labely)
      endif

      if(an.eq.'5') then
         labely='Pre-fit residuals [usec]'
         labelx='TOA number'
         call plt(xindex,1.0d0,xnpts,resid1,-r1,r1,npts,rct,p0,
     +            error,labelx,labely)
      endif

      if(an.eq.'6') then
         labely='Post-fit residuals [usec]'
         labelx='TOA number'
         call plt(xindex,1.0d0,xnpts,resid2,-r2,r2,npts,rct,p0,
     +            error,labelx,labely)
      endif

      if(an.eq.'7') then
         labely='Pre-fit residuals [usec]'
         labelx='Day of year'
         call plt(dayofyr,0.d0,365.25d0,resid1,-r1,r1,npts,rct,p0,
     +            error,labelx,labely)
      endif

      if(an.eq.'8') then
         labely='Post-fit residuals [usec]'
         labelx='Day of year'
         call plt(dayofyr,0.d0,365.25d0,resid2,-r2,r2,npts,rct,p0,
     +            error,labelx,labely)
      endif

      go to 12

 100  continue
      call pgend(.false.)
      
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine plt(x,xmin,xmax,y,ymin,ymax,npts,ct,p0,error,
     +               labelx,labely)
      
      implicit real*8 (a-h,o-z)
      real*8 x(npts),y(npts),ct(npts), error(npts), xgap, ygap
      character*26 letter
      character labelx*(*), labely*(*), label*120
      integer mark, markOld

      real xp(6000), yp(6000), dy1(6000), dy2(6000)

      xgap=(xmax-xmin)/30.0d0

      ygap=(ymax-ymin)/30.0d0

      call pgpage
      call pgenv(sngl(xmin-xgap),sngl(xmax+xgap),
     +           sngl(ymin-ygap),sngl(ymax+ygap),0,0)
      call pgsls(2)
      call pgmove(sngl(xmin-xgap), 0.0)
      call pgdraw(sngl(xmax+xgap), 0.0)
      call pgsls(1)

      markOld=0
      sum=0.
      sq=0.
      do i=1,npts
         sum=sum+y(i)
         sq=sq+y(i)**2
         mark=mod(ifix(sngl(ct(i))),26)+1
         xp(i)=sngl(x(i))
         yp(i)=sngl(y(i))
         dy1(i)=yp(i)+sngl(error(i))
         dy2(i)=yp(i)-sngl(error(i))
         if (markold.ne.mark) then
            call pgptext(xp(i),sngl(ymax+1.5*ygap),0.0,0.5,'|')
            markold=mark
         endif
      enddo
      call pgpoint(npts,xp,yp,17)
      call pgerry(npts,xp,dy1,dy2,1.0)

      rms=dsqrt(sq/npts - (sum/npts)**2)


      write(label,1050) ymax*1.0e-6/p0,ymax,rms*1.0e-6/p0,rms
 1050 format('Max:',f15.9,' P0,',f9.3,' us.',5x,'RMS:',f15.9,
     +     ' P0,',f9.3,' us.')
               
      call pglabel(labelx, labely,label)

      return
      end