SUBROUTINE CFI ( DRES, PHASE, VSQ, MAXSCAN, NPTS, FRINGE_OFFSET,
     $ ISCN,ISTAR,IBASE,LAMBDA,DELAY,D2RMS,N2C,UNIT,imark,noplot,scanid,
     $ dispf)
C
C     Central fringe identification.   Automated version.
C
C     This version should be used if channel 2 is the 700/40 filter.
C     In good seeing, the 500-700 phase should give good fringe
C     identification for all of the data.
C
C             VERSION 1.1    6 September 1991, David Mozurkewich
C
      IMPLICIT UNDEFINED (A-Z)
      SAVE
      REAL *4 PI
      PARAMETER ( PI = 3.1415926535897932384626 )
      INTEGER *4 NMAX
      PARAMETER ( NMAX = 4000 )
      INTEGER *4 NPTS, MAXSCAN, ISCN, ISTAR, IBASE, ITER, FSET
      INTEGER *4 SYMBOL(6), IMARK(NMAX), IPLOT, COUNT(5), N2C
      INTEGER *4 IC, N, I, J, K, IERR, MODE, UNIT
      REAL *4 DRES(*), PHASE(4,MAXSCAN), VSQ(MAXSCAN)
      REAL *4 D1C(NMAX), DATA(5), V2(5), FRINGE_OFFSET(3),dispf(4)
      REAL *4 D2C(NMAX,4), PEAK(2), SPACING(9,2)
      REAL *4 PHI(NMAX,4),FRINGE(4),RMS,XBAR(2),YBAR(2),DISP(4),dispo(4)
      REAL *4 XPLOT(NMAX), YPLOT(NMAX), DELAY(2), D2RMS(2), VIS(NMAX)
      REAL *4 XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX, CHISQ
      REAL *4 XLO, XHI, YLO, YHI, TEMP, DAVG, DRMS, V0, V1
      REAL *4 DELTA(2), RESID, RX, RY, LAMBDA(*), VLIMIT, csize
      CHARACTER *1 CH, ANS
      CHARACTER *40 TITLE(5), LABEL(6), FNAME, blank
      character scanid*11
      real*4 vrlim,ymean,d2rms1,d2rms2,delay1,delay2,ch2d2off,rrb,null
      integer*4 nymean,nimark6,iother
      parameter (vrlim=1.2)  ! limit for visibility ratio 0 / +/-2
      parameter (rrb= 1.10)  ! limit for visibility ratio -1 / +1
      LOGICAL FIRST, LPLOT, HCOPY, PLOT, stopplot, noplot, other, repeat
      logical firstdisp,fitdisp
      real*8 evapor
C-----------------------------------------------------------------------
      INCLUDE 'HARDCOPY.INC'
C-----------------------------------------------------------------------
      INTEGER *4 PGBEGIN, PGCURSE
      EXTERNAL PGBEGIN, PGCURSE
c     DATA SYMBOL /  2,    3,    6,   13,  16,  -1 /
      DATA SYMBOL /  2,    4,   17,    6 ,  5,  -1 /
      DATA MODE   /   2 /
      DATA FIRST, LPLOT, PLOT  / .TRUE., .false., .FALSE. /
      data null / 0./
      HCOPY = .FALSE.
      blank='                                        '
C----------------------------------------------------------------------
      IF ( NPTS .GT. NMAX ) THEN
          WRITE(6,*) ' TOO MANY POINTS IN THIS SCAN '
          WRITE(6,*) ' NUMBER OF POINTS = ', NPTS
          WRITE(6,*) ' CHANGE NMAX AND RECOMPILE '
          STOP
      END IF
C----------------------------------------------------------------------
      IF ( (LAMBDA(1) .EQ. 700.) .AND. 
     $     (LAMBDA(2) .EQ. 700.) .AND.
     $     (LAMBDA(3) .EQ. 500.) .AND.
     $     (LAMBDA(4) .EQ. 550.) ) THEN
          FSET = 1
      ELSE IF ( (LAMBDA(1) .EQ. 700.) .AND. 
     $          (LAMBDA(2) .EQ. 800.) .AND.
     $          (LAMBDA(3) .EQ. 500.) .AND.
     $          (LAMBDA(4) .EQ. 550.) ) THEN
          FSET = 2
      ELSE IF ( (LAMBDA(1) .EQ. 700.) .AND. 
     $          (LAMBDA(2) .EQ. 800.) .AND.
     $          (LAMBDA(3) .EQ. 550.) .AND.
     $          (LAMBDA(4) .EQ. 500.) ) THEN
       write(12,*)'*** Warning: channels 3 and 4 filter numbers',
     $            ' may be exchanged!!'
       write(12,*)'*** Stop now!'
       write(6,*)'*** Warning: channels 3 and 4 filter numbers',
     $            ' may be exchanged!!'
       write(6,*)'*** Stop now!'
       stop
      END IF
      FRINGE(1) = LAMBDA(4) / 2000.
      FRINGE(2) = LAMBDA(3) / 2000.
      fringe(4) = lambda(1) / 2000.
      CALL AIR ( FSET, DISP, SPACING, fitdisp, dispf )
      firstdisp=.true.
C
      WRITE(12,*) ' STARTING DATA FOR ', ISCN, ISTAR, IBASE
      WRITE(12,'(A,4F5.0)') ' FILTER WAVELENGTHS = ', (LAMBDA(k),k=1,4)
      WRITE(12,*) ' FILTER SET = ', FSET
      WRITE(12,*) ' DISP = ', DISP
      WRITE(12,*) ' SPACING = '
      WRITE(12,'(9F8.2)') (SPACING(k,1),k=1,9)
      WRITE(12,'(9F8.2)') (SPACING(k,2),k=1,9)
      WRITE(12,1200) LAMBDA(4), LAMBDA(2)
      WRITE(12,1200) LAMBDA(3), LAMBDA(2)
      WRITE(12,1200) LAMBDA(3), LAMBDA(4)
C
      WRITE ( TITLE(1), 1200) LAMBDA(4), LAMBDA(2)
      WRITE ( TITLE(2), 1200) LAMBDA(3), LAMBDA(2)
      WRITE ( TITLE(3), 1200) LAMBDA(3), LAMBDA(4)
      write ( title(4), 1200) lambda(1), lambda(2)
c1200 FORMAT ( F4.0, ' nm REFERENCED TO ', F4.0, ' nm' )
 1200 FORMAT ( 'd(',f4.0,'nm) - d(',f4.0,'nm) [microns]' )
C----------------------------------------------------------------------
C                      Set up graphics here
      IF ( FIRST ) THEN
          FIRST = .FALSE.
          stopplot=.false.
          if(noplot)stopplot=.true.
          NPLOTS = 0
          PLOT_DEV = '/vps'
C      ELSE
C          GO TO 100
      END IF
   80 CONTINUE
c     IF ( HCOPY .and..not.stopplot) THEN
c         IERR = 0
c         DO WHILE ( IERR .NE. 1 ) THEN
c             NPLOTS = NPLOTS + 1
c             WRITE(PLOTFILE(1:8),'(A5,I3.3)') 'PLOT.', NPLOTS
c             WRITE(PLOTFILE(1:8),'(A2,I3.3,a3)') 'PL', NPLOTS,'.ps'
c             WRITE(6,*) ' TRYING TO SEND HARDCOPY TO FILE ', PLOTFILE
c             IERR = PGBEGIN(0, PLOTFILE // PLOT_DEV, 3, 4 )
c         END DO
c         WRITE(6,*) ' HARDCOPY STARTED '
c     ELSE if(.not.stopplot)then
c         IERR = PGBEGIN(0, '?', 3, 2)
c     END IF
c     IF ( IERR .NE. 1 .and..not.stopplot) THEN
c         WRITE(6,*) ' PGBEGIN FAILED.  RETURN CODE = ', IERR
c         RETURN
c     END IF
  100 CONTINUE
C-----------------------------------------------------------------------
C     Calculate the mean one color delay.
C
      DAVG = 0.
      DRMS = 0.
      DO I = 1, NPTS
          D1C(I) = DRES(I) + PHASE(1,I) + phase(2,i) 
          DAVG = DAVG + D1C(I)
          DRMS = DRMS + D1C(I)*D1C(I)
      END DO
      DAVG = DAVG / FLOAT( NPTS )
      DRMS = SQRT ( DRMS/FLOAT(NPTS) - DAVG*DAVG)
      WRITE(12,'(A,F10.1,A,F6.1)') ' 1-C DELAY ', DAVG, ' +/- ', DRMS
C-----------------------------------------------------------------------
C     Calculate the phase differences and the two color delays.
C
      DO I = 1, NPTS
c     write(6,*)' tracking channel phase=',phase(1,i)
          DO J= 1, 2
              PHI(I,J) = PHASE(5-J,I) - PHASE(2,I)
              IF ( PHI(I,J) .GT. FRINGE(J) ) THEN
                  PHI(I,J) = PHI(I,J) - 2.*FRINGE(J)
              ELSE IF ( PHI(I,J) .LT. -FRINGE(J) ) THEN
                  PHI(I,J) = PHI(I,J) + 2.*FRINGE(J)
              END IF
              D2C(I,J) = D1C(I) - DISP(J) * PHI(I,J)
          END DO
          PHI(I,3) = PHI(I,2) - PHI(I,1)
          D2C(I,3) = D1C(I) - DISP(3) * PHI(I,3)
          phi(i,4) = - phase(2,i)
          if ( phi(i,4) .gt. fringe(4) ) then
              phi(i,4) = phi(i,4) - 2.*fringe(4)
          else if ( phi(i,4) .lt. -fringe(4) ) then
              phi(i,4) = phi(i,4) + 2.*fringe(4)
          endif
          if(fset.eq.2)then
c          d2c(i,4) = d1c(i) - disp(4) * phi(i,4)
c          Instead of projecting onto delay axis, proj. on phase axis
           d2c(i,4) = phi(i,4) - 1./disp(4)*(d1c(i)-davg)
          else
           d2c(i,4) = d1c(i)
          endif
      END DO
C-----------------------------------------------------------------------
C     Identify the fringes.
C
      other=.false.
      iother=0
      if(scanid(8:8).eq.'1')then
       other=.true.
      endif
c     scanid='GOOD ID0000'    ! passed to CFI now from TWOC
   90 PEAK(1) = DAVG
      PEAK(2) = DAVG - FRINGE_OFFSET(ibase)
C     write(12,*) ' identify mode = ', mode
      do i=1,4
       dispo(i)=disp(i)
      enddo
      if(fset.eq.2)then
      IF ( MODE .EQ. 1 ) THEN
          WRITE(12,'(A,2F9.2)') ' PEAK POSITIONS = ', PEAK(1), PEAK(2)
          CALL IDENT_1( D2C, NMAX, 4, SPACING, PEAK, IMARK, NPTS)
      ELSE IF ( MODE .EQ. 2 ) THEN
          WRITE(12,'(A,2F9.2)') ' PEAK POSITIONS = ', PEAK(1), PEAK(2)
          CALL IDENT_2(D2C,VSQ,NMAX,4,SPACING,PEAK,FRINGE_OFFSET,
     $     IMARK,NPTS,PLOT,scanid,other,disp,phi,d1c,ibase,fitdisp)
      ELSE IF ( MODE .EQ. 3 ) THEN
          WRITE(12,'(A,2F9.2)') ' PEAK POSITIONS = ', PEAK(1), PEAK(2)
          CALL IDENT_3(D2C,VSQ,NMAX,4,SPACING,PEAK,FRINGE_OFFSET,
     $     IMARK,NPTS,PLOT,scanid,other,disp,phi,d1c,ibase,fitdisp)
      ELSE
          WRITE(12,*) ' YOU CANNOT GET HERE BOZO, mode = ', mode
      END IF
      elseif(fset.eq.1)then
          WRITE(12,'(A,2F9.2)') ' PEAK POSITIONS = ', PEAK(1), PEAK(2)
          CALL IDENT_4(D2C,VSQ,NMAX,4,SPACING,PEAK,FRINGE_OFFSET,
     $     IMARK,NPTS,PLOT,scanid,other,disp,phi,d1c,ibase,fitdisp)
      else
          WRITE(12,*) ' YOU CANNOT GET HERE BOZO, mode = ', mode
      endif
      do i=1,4
       if(disp(i).ne.dispo(i).and.firstdisp)then
        firstdisp=.false.
        goto 100
       endif
      enddo
      DO I = 1, NPTS
          IMARK(I) = IMARK(I) + 3
      END DO
C-----------------------------------------------------------------------
C     Determine the mean visibility for each fringe.
C
      DO J = 1, 5
          COUNT(J) = 0
          V2(J) = 0.
          DO I = 1, NPTS
              IF ( IMARK(I) .EQ. J ) THEN
                  COUNT(J) = COUNT(J) + 1
                  V2(J) = V2(J) + VSQ(I)
              END IF
          END DO
          IF ( COUNT(J) .GT. 1 ) THEN
              V2(J) = V2(J) / COUNT(J)
          END IF
      END DO
      write(12,'(a,5i5)')  ' Number of points per fringe ', count
      write(12,'(a,5f5.2)')' Mean visibility of fringe   ', v2
c     if(v2(2).lt.v2(4)*rrb.and.count(2).ge.10.and.count(4).ge.10)then
c      write(12,*)' *** Bad ID! -1 fringe vis less than +1 fringe vis!'
c      if(.not.other.and.ch.ne.'o'.and.ch.ne.'O')then
c       write(12,*)' *** Try other peak.'
c       other=.true.
c       scanid(8:11)='1000'
c       goto 90
c      endif
c     endif
c     if(v2(2).lt.v2(4).and.count(2).ge.10.and.count(4).ge.10)
c    $                     then
c     scanid(11:11)='1'
c     endif
      if(count(1).ge.10.and.v2(3)/v2(1).le.vrlim)scanid(11:11)='1'
      if(count(5).ge.10.and.v2(3)/v2(5).le.vrlim)scanid(11:11)='1'
      if(count(2).ge.10.and.count(4).ge.10)then
       if(2.*v2(3)/(v2(2)+v2(4)).le.rrb)scanid(11:11)='1'
      endif
      if(scanid(11:11).eq.'1')then
       write(12,*)' *** Bad ID! Ratios of fringe vis. incorrect'
c      if(.not.other.and.ch.ne.'o'.and.ch.ne.'O')then
c       write(12,*)' *** Try other peak.'
c       other=.true.
c       scanid(8:11)='1000'
c       goto 90
c      endif
      endif
C-----------------------------------------------------------------------
C     Edit data based on visibility.  The visibility limit for the
C     central fringe is taken as the mean value of fringes +1 and -1.
C
      IF ( (COUNT(2).GT.5) .AND. (COUNT(4).GT.5) ) THEN
          VLIMIT = 0.5*( V2(2) + V2(4) )
      ELSE
          VLIMIT = 0.5 * V2(3)
      END IF
      DO I = 1, NPTS
          IF ( (IMARK(I).EQ.3) .AND. (VSQ(I).LT.VLIMIT) ) THEN
              IMARK(I) = 6
          END IF
      END DO
C----------------------------------------------------------------------
C
C     Recalculate the 2-color delays for points not on the central fringe
C
c     DO I = 1, NPTS
c         DO J= 1, 2
c             if(j.eq.1)then
c             PHI(I,J) = PHASE(5-J,I) - PHASE(2,I)
c    $        -float(imark(i)-3)*(lambda(2)-lambda(5-j))/1000.
c             else
c             PHI(I,J) = PHASE(5-J,I) - PHASE(2,I)
c    $        +float(imark(i)-3)*(lambda(2)-lambda(5-j))/1000.
c             endif
c             do k=1,2
c             IF ( PHI(I,J) .GT. FRINGE(J) ) THEN
c                 PHI(I,J) = PHI(I,J) - 2.*FRINGE(J)
c             ELSE IF ( PHI(I,J) .LT. -FRINGE(J) ) THEN
c                 PHI(I,J) = PHI(I,J) + 2.*FRINGE(J)
c             END IF
c             enddo
c             D2C(I,J) = D1C(I) - DISP(J) * PHI(I,J)
c         END DO
c         PHI(I,3) = PHI(I,2) - PHI(I,1)
c         D2C(I,3) = D1C(I) - DISP(3) * PHI(I,3)
c     ENDDO
C----------------------------------------------------------------------
c     Recalculate phases for central fringe points
c     (does not work)
c
c     DO I = 1, NPTS
c         DO J= 1, 2
c           if(imark(i).eq.3)then
c             if(phase(1,i).gt.fringe(j))then
c              D2C(I,J) = D1C(I) - DISP(J) * (PHI(I,J)+fringe(j))
c             else if(phase(1,i).lt.-fringe(j))then
c              D2C(I,J) = D1C(I) - DISP(J) * (PHI(I,J)-fringe(j))
c             endif
c           endif
c         END DO
c     ENDDO
C----------------------------------------------------------------------
C     Determine 2-color delays, which will be edited after this
C
      DO IPLOT = 1, 2
          DELAY(IPLOT) = 0.
          D2RMS(IPLOT) = 0.
          N2C = 0
          ch2d2off=fringe(2)*disp(2)
          DO I = 1, NPTS
            IF ( IMARK(I) .EQ. 3 ) THEN
               if(iplot.eq.1)then
                DELAY(IPLOT) = DELAY(IPLOT) + D2C(I,IPLOT)
                D2RMS(IPLOT) = D2RMS(IPLOT) + D2C(I,IPLOT)**2
                N2C = N2C + 1
               else
             if(d2c(i,iplot).lt.delay(1)-fringe_offset(ibase)-ch2d2off)
     $          d2c(i,iplot)=d2c(i,iplot)+2.0*ch2d2off
             if(d2c(i,iplot).gt.delay(1)-fringe_offset(ibase)+ch2d2off)
     $          d2c(i,iplot)=d2c(i,iplot)-2.0*ch2d2off
                DELAY(IPLOT) = DELAY(IPLOT) + D2C(I,IPLOT)
                D2RMS(IPLOT) = D2RMS(IPLOT) + D2C(I,IPLOT)**2
                N2C = N2C + 1
               endif
            END IF
          END DO
          IF ( N2C .GT. 1 ) THEN
              DELAY(IPLOT) = DELAY(IPLOT) / FLOAT(N2C)
              D2RMS(IPLOT) = D2RMS(IPLOT)/FLOAT(N2C) - DELAY(IPLOT)**2
              IF ( D2RMS(IPLOT) .GT. 0. ) THEN
                  D2RMS(IPLOT) = SQRT ( D2RMS(IPLOT) )
              END IF
          END IF
      ENDDO
      delay1=delay(1)
      d2rms1=d2rms(1)
      delay2=delay(2)
      d2rms2=d2rms(2)
C
C     Plot BOTH dispersion diagrams.
C     Determine the final two color delays.
C
      IF ( HCOPY .and..not.stopplot) THEN
          IERR = 0
          DO WHILE ( IERR .NE. 1 ) THEN
              NPLOTS = NPLOTS + 1
c             WRITE(PLOTFILE(1:8),'(A5,I3.3)') 'PLOT.', NPLOTS
              WRITE(PLOTFILE(1:8),'(A2,I3.3,a3)') 'PL', NPLOTS,'.ps'
              WRITE(6,*) ' TRYING TO SEND HARDCOPY TO FILE ', PLOTFILE
              IERR = PGBEGIN(0, PLOTFILE // PLOT_DEV, 2, 3 )
          END DO
          WRITE(6,*) ' HARDCOPY STARTED '
      ELSE if(.not.stopplot)then
          IERR = PGBEGIN(0, '/ega', 2, 2)
      END IF
      IF ( IERR .NE. 1 .and..not.stopplot) THEN
          WRITE(6,*) ' PGBEGIN FAILED.  RETURN CODE = ', IERR
          RETURN
      END IF
c
      DO IPLOT = 1, 2
C         write(12,*) ' starting plot number, iplot = ', iplot
C         write(12,*) ' number of points,      npts = ', npts
          XMAX = -1.E19
          XMIN =  1.E19
          YMAX = -1.E19
          YMIN =  1.E19
          DELAY(IPLOT) = 0.
          D2RMS(IPLOT) = 0.
          N2C = 0
          nimark6=0
          DO I = 1, NPTS
              XPLOT(I) = D1C(I) - PEAK(IPLOT)
              YPLOT(I) = PHI(I,IPLOT)
c             yplot(i) = d2c(i,iplot)
              IF ( IMARK(I) .EQ. 3 ) THEN
                if(iplot.eq.1)then
                  DELAY(IPLOT) = DELAY(IPLOT) + D2C(I,IPLOT)
                  D2RMS(IPLOT) = D2RMS(IPLOT) + D2C(I,IPLOT)**2
                  N2C = N2C + 1
                else if(iplot.eq.2
     $            .and.abs(d2c(i,2)-delay2).lt.d2rms1*3)then
                  DELAY(IPLOT) = DELAY(IPLOT) + D2C(I,IPLOT)
                  D2RMS(IPLOT) = D2RMS(IPLOT) + D2C(I,IPLOT)**2
                  N2C = N2C + 1
                else if(iplot.eq.2
     $            .and.abs(d2c(i,2)-delay2).ge.d2rms1*3)then
                  imark(i)=6
                  nimark6=nimark6+1
                endif
              END IF
              XMIN = MIN ( XMIN, XPLOT(I) )
              XMAX = MAX ( XMAX, XPLOT(I) )
              YMIN = MIN ( YMIN, YPLOT(I) )
              YMAX = MAX ( YMAX, YPLOT(I) )
          END DO
C         write(12,*) ' number of 2c points,   n2c  = ', n2c
          if(iplot.eq.2)
     $    write(12,'(a,10x,i5)')' Number of 0 fringe pts. ed. ',nimark6
          IF ( N2C .GT. 1 ) THEN
              DELAY(IPLOT) = DELAY(IPLOT) / FLOAT(N2C)
              D2RMS(IPLOT) = D2RMS(IPLOT)/FLOAT(N2C) - DELAY(IPLOT)**2
              IF ( D2RMS(IPLOT) .GT. 0. ) THEN
                  D2RMS(IPLOT) = SQRT ( D2RMS(IPLOT) )
c             write(12,*)' iplot=',iplot,' delay=',delay(iplot),
c    $                  ' d2rms=',d2rms(iplot),' n2c=',n2c
              END IF
          END IF
          LABEL(1) = 'residual one-color delay [microns]'
          WRITE(LABEL(2),1203) ISCN, ISTAR, IBASE
          LABEL(3) = TITLE(IPLOT)
         if(.not.stopplot)then         ! don't plot if(stopplot)
          CALL PGSCI ( 5 )
          CALL PGSLS ( 1 )
          CALL PGSCH  ( 2.0 )
          CALL PGRNGE ( XMIN, XMAX, XLO, XHI )
          CALL PGRNGE ( YMIN, YMAX, YLO, YHI )
          YLO = MIN ( YMIN, -FRINGE(IPLOT) )
          YHI = MAX ( YMAX,  FRINGE(IPLOT) )
          call pgqch(csize)
          call pgsch(1.4)
          CALL PGENV  ( XLO, XHI, YLO, YHI, 0, 1 )
          call pgsch(csize)
          call pgsch(1.5)
          CALL PGMTEXT ( 'B', 3.5, 0.5, 0.5, LABEL(1) )
          CALL PGMTEXT ( 'L', 3.0, 0.5, 0.5, LABEL(3) )
          if(iplot.eq.1)CALL PGMTEXT ( 'T', 1.5, 0.5, 0.5, 'a' )
          if(iplot.eq.2)CALL PGMTEXT ( 'T', 1.5, 0.5, 0.5, 'b' )
c         CALL PGMTEXT ( 'T', 1.5, 0.5, 0.5, blank )
          CALL PGSCH ( 1.0 )
          DO I = 1, NPTS
              CALL PGSCI ( IMARK(I) )
              CALL PGPOINT( 1, XPLOT(I), YPLOT(I), SYMBOL(imark(i)) )
          END DO
          IF ( LPLOT ) THEN
              CALL PGSLS ( 1 )
              CALL PGSCI ( 5 )
              XBAR(1) = XMIN
              XBAR(2) = XMAX
              DO I = 1, 9
                  YBAR(1) = (XBAR(1)-SPACING(I,IPLOT)) / DISP(IPLOT)
                  YBAR(2) = (XBAR(2)-SPACING(I,IPLOT)) / DISP(IPLOT)
                  CALL PGLINE ( 2, XBAR, YBAR )
                  CALL PGSLS ( 2 )
              END DO
          END IF
          CALL PGSLS ( 1 )
         endif               ! endif of if(.not.stopplot)
      END DO
c++++++++++++++++++++++++ This code now in twoc.f
c             If the delay difference is off by more than 2 micron from
c             fringe_offset, we change the id!
c
C             if(ch.ne.'o'.and.ch.ne.'O')then
C             if(abs(delay(1)-delay(2)-fringe_offset(ibase)).gt.2.0)then
C              other=.true.
C              iother=iother+1
C              if(iother.eq.1)then
C               write(12,'(a,a)')' *** Delay offset too large, ',
C    $                           'try other peak!'
C               if(.not.stopplot)then
C                repeat=.false.
C                scanid(8:8)='1'
C                scanid(10:11)='00'
C                goto 90
C               else
C                repeat=.true.
C               endif
C              endif
C             endif
C             endif
c+++++++++++++++++++++++
c
c     Plot 700 nm referenced to 800 nm
c
         if(.not.stopplot)then         ! don't plot if(stopplot)
          XMAX = -1.E19
          XMIN =  1.E19
          YMAX = -1.E19
          YMIN =  1.E19
          ymean=0
          nymean=0
          DO I = 1, NPTS
              XPLOT(I) = d1c(i)
              YPLOT(I) = PHI(I,4)
              yplot(i) = d2c(i,4)
c             yplot(i) = phase(2,i)
              if(imark(i).eq.3)then
               ymean=ymean+yplot(i)
               nymean=nymean+1
              endif
              XMIN = MIN ( XMIN, XPLOT(I) )
              XMAX = MAX ( XMAX, XPLOT(I) )
              YMIN = MIN ( YMIN, YPLOT(I) )
              YMAX = MAX ( YMAX, YPLOT(I) )
          END DO
          if(nymean.gt.0)then
           ymean=ymean/float(nymean)
           write(12,*)' Mean ch 2 phase in green fringe is: ',ymean
          else
           ymean=0
           write(12,*)' *** Warning: no central fringe points left!'
          endif
          CALL PGSCI ( 5 )
          CALL PGSLS ( 1 )
          CALL PGSCH  ( 2.0 )
          CALL PGRNGE ( XMIN, XMAX, XLO, XHI )
          CALL PGRNGE ( YMIN, YMAX, YLO, YHI )
c         YLO = MIN ( YMIN, -FRINGE(4) )
c         YHI = MAX ( YMAX,  FRINGE(4) )
          call pgqch(csize)
          call pgsch(1.4)
          CALL PGENV  ( XLO, XHI, YLO, YHI, 0, 1 )
          call pgsch(csize)
          call pgsch(1.5)
          CALL PGMTEXT ( 'B', 3.5, 0.5, 0.5, LABEL(1) )
          CALL PGMTEXT ( 'L', 3.0, 0.5, 0.5, title(4) )
          CALL PGMTEXT ( 'T', 1.5, 0.5, 0.5, 'c' )
          CALL PGSCH ( 1.0 )
          xbar(1)=xlo
          xbar(2)=xhi
          ybar(1)=ymean
          ybar(2)=ymean
          call pgline(2,xbar,ybar)
          DO I = 1, NPTS
              CALL PGSCI ( IMARK(I) )
              CALL PGPOINT( 1, XPLOT(I), YPLOT(I), SYMBOL(imark(i)) )
          END DO
         endif               ! else  of if(.not.stopplot)
         do iplot=1,5
          ymean=0
          nymean=0
          DO I = 1, NPTS
              yplot(i) = d2c(i,4)
              if(imark(i).eq.iplot)then
               ymean=ymean+yplot(i)
               nymean=nymean+1
              endif
          ENDDO
          if(nymean.ne.0)then
           ymean=ymean/float(nymean)
          else
           ymean=0
          endif                                             
          write(12,*)' Mean ch 2 phase in fringe ',iplot-3,' is: ',ymean
          if(iplot.eq.3)then
c          if(ymean.gt.0.07.or.ymean.lt.-0.07)scanid(10:10)='1'
          endif
         enddo
C-----------------------------------------------------------------------
C     Plot the visibility amplitude.  It should be highest for fringes
C     0 or 1 or -1.
C
C-----------------------------------------------------------------------
 1203     FORMAT('#INT =', I4, ' STAR =', I5, ' BL =', I3 )
C-----------------------------------------------------------------------
C     Plot visibility amplitude squared versus time (integraton).
C     Plot the one color delay versus time (integraton).
C     Plot the two color delays versus time (integraton).
C
      if(.not.stopplot)then
      DO IPLOT = 1, 4
          XMIN =  0.
          XMAX = -1.E19
          YMIN =  1.E19
          YMAX = -1.E19
          K = 0
          DO I = 1, NPTS
              IF ( ( IPLOT .LE. 2 ) .OR. ( IMARK(I) .EQ. 3 ) ) THEN
                  K = K + 1
                  XPLOT(K) = I
                  IF ( IPLOT .EQ. 1 ) THEN      ! squared visibility
                      YPLOT(K) = VSQ(I)
                  ELSE IF ( IPLOT .EQ. 2 ) THEN   ! one-color delay
                      YPLOT(K) = D1C(I) - PEAK(IPLOT)
                  ELSE                            ! two-color delays
                      YPLOT(K) = D2C(I,IPLOT-2) - DAVG
                  END IF
                  XMAX = MAX ( XMAX, XPLOT(K) )
                  YMIN = MIN ( YMIN, YPLOT(K) )
                  YMAX = MAX ( YMAX, YPLOT(K) )
              END IF
          END DO
          if(iplot.eq.2)then
           call pgend
           IF ( HCOPY .and..not.stopplot) THEN
               IERR = 0
               DO WHILE ( IERR .NE. 1 ) THEN
                  NPLOTS = NPLOTS + 1
c                 WRITE(PLOTFILE(1:8),'(A5,I3.3)') 'PLOT.', NPLOTS
                  WRITE(PLOTFILE(1:8),'(A2,I3.3,a3)') 'PL', NPLOTS,'.ps'
                  WRITE(6,*)
     &            ' TRYING TO SEND HARDCOPY TO FILE ', PLOTFILE
                  IERR = PGBEGIN(0, PLOTFILE // PLOT_DEV, 3, 4 )
               END DO
               WRITE(6,*) ' HARDCOPY STARTED '
           ELSE if(.not.stopplot)then
               IERR = PGBEGIN(0, '/ega', 3, 2)
           END IF
           IF ( IERR .NE. 1 .and..not.stopplot) THEN
               WRITE(6,*) ' PGBEGIN FAILED.  RETURN CODE = ', IERR
               RETURN
           END IF
          endif
          CALL PGSCI ( 5 )
          CALL PGSCH  ( 2.0 )
          CALL PGRNGE ( XMIN, XMAX, XLO, XHI )
          CALL PGRNGE ( YMIN, YMAX, YLO, YHI )
          if(iplot.ge.2)then   ! same scale for two-color plots
           ylo=-15.
           yhi=15.
          endif
c         if(iplot.ne.2)then        ! comment if/endif for 1-color plot
           call pgqch(csize)
           call pgsch(1.4)
           CALL PGENV  ( XLO, XHI, YLO, YHI, 0, 1 )
           call pgsch(csize)
           call pgsch(1.5)
           CALL PGMTEXT ( 'B', 3.5, 0.5, 0.5, 'coh. average no.' )
c         endif
          IF ( IPLOT .EQ. 1 ) THEN
           CALL PGMTEXT ('L', 3.0, 0.5, 0.5, 'squared vis., channel 1')
           CALL PGMTEXT ('T', 1.5, 0.5, 0.5, 'd' )
          ELSE IF ( IPLOT .EQ. 2 ) THEN    ! comment to not plot 1-color
           CALL PGMTEXT ('L',3.0,0.5,0.5, 'one-color delay [microns]')
           CALL PGMTEXT ('T',1.5,0.5,0.5, 'a' )
          ELSE IF ( IPLOT .GE. 3 ) THEN
           if(iplot.eq.3)then
            CALL PGMTEXT('L',3.0,0.5,0.5, 'two-color delay (550/800)')
            CALL PGMTEXT('T',1.5,0.5,0.5, 'b' )
           else
            CALL PGMTEXT('L',3.0,0.5,0.5, 'two-color delay (500/800)')
            CALL PGMTEXT('T',1.5,0.5,0.5, 'c' )
           endif
c          CALL PGMTEXT ( 'T', 1.5, 0.5, 0.5, TITLE(IPLOT-2) )
          END IF
          CALL PGSCH ( 1.0 )
          IF ( IPLOT .GT. 2 ) THEN
              CALL PGSCI ( 3 )
c             do i=1,k
c              call pgsci(imark(i))
               CALL PGPOINT( K, XPLOT, YPLOT, SYMBOL(3) )
c              call pgpoint(1, xplot(i),yplot(i),symbol(imark(i)))
c             enddo
              xbar(1)=xlo
              xbar(2)=xhi
              ybar(1)=delay(iplot-2)-davg
              ybar(2)=ybar(1)
              call pgsci(5)
              call pgline(2,xbar,ybar)
          ELSE  ! if(iplot.ne.2)then         ! comment 'if' to plot 1-color
              DO I = 1, K
                  CALL PGSCI ( IMARK(I) )
                  CALL PGPOINT( 1, XPLOT(I), YPLOT(I), SYMBOL(imark(i)))
              END DO
          END IF
      END DO
      endif
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      IF ( HCOPY ) THEN
          HCOPY = .NOT. HCOPY
          CALL PGEND
          WRITE(6,*) ' HARDCOPY FINISHED '
          GO TO 80
      END IF
  200 CONTINUE
      if(.not.stopplot)then
      CALL LOCATE(0,0)
      call ega_restore_default
      CALL WRITE_STRING ( 'Command:' )
      CALL PGCURSE ( RX, RY, CH )
      IF ( CH .EQ. '?' ) THEN
          CALL LOCATE (0, 1)
          CALL WRITE_STRING ( '?    Displays this Help       ' )
          CALL LOCATE (0, 2)
          CALL WRITE_STRING ( 'D    Change the Dispersion    ' )
          CALL LOCATE (0, 3)
          CALL WRITE_STRING ( 'L    Toggle Line option       ' )
          CALL LOCATE (0, 4)
          CALL WRITE_STRING ( 'F    Enter Fringe Offset      ' )
          CALL LOCATE (0, 5)
          CALL WRITE_STRING ( 'M    Toggle Fringe ID Mode    ' )
          CALL LOCATE (0, 6)
          CALL WRITE_STRING ( 'I    Show Fringe Information  ' )
          CALL LOCATE (0, 7)
          CALL WRITE_STRING ( 'H    Hardcopy                 ' )
          CALL LOCATE (0, 8)
          CALL WRITE_STRING ( 'P    Toggle plot option       ' )
          CALL LOCATE (0, 9)
          CALL WRITE_STRING ( 'N    Go to next scan          ' )
          CALL LOCATE (0,10)
          CALL WRITE_STRING ( 'Q    Quit program             ' )
          call locate (0,11)
          call write_string ( 'S    Stop plotting            ' )
          call locate (0,12)
          call write_string ( 'O    Try other peak           ' )
          call locate (0,13)
          call write_string ( 'E    Water vapor pressure     ' )
          call locate (0,14)
          call write_string ( 'W    Write 2-color data points' )
          GO TO 200
      ELSE IF ( CH .EQ. 'D' .or. ch .eq. 'd') THEN
          WRITE(6,*) ' INPUT CHANNEL, DISPERSION '
          WRITE(6,'(A,4F8.1)') ' values are ', DISP
          READ (5,*) I, TEMP
c         IF ( ( I .EQ. 1 ) .OR. ( I .EQ. 2 ) ) THEN
              DISP(I) = TEMP
c         END IF
      else if ( ch .eq. 'e' .or. ch .eq. 'E') then
          write(6,*) ' INPUT partial water vapor pressure'
          read(5,*)evapor
          CALL AIR2 ( FSET, DISP, SPACING, evapor )
      ELSE IF ( CH .EQ. 'L' .or. ch .eq. 'l') THEN
          LPLOT = .NOT. LPLOT
      else if ( ch .eq. 'w' .or. ch .eq. 'W') then
          do i=1,npts
           if(imark(i).eq.3)
     &     write(12,*)d2c(i,1)-davg,null,d2c(i,2)-davg,null
          enddo
          goto 200
      ELSE IF ( CH .EQ. 'F' .or. ch .eq. 'f') THEN
          WRITE(6,'(A,$)') ' INPUT FRINGE OFFSET '
          READ (5,*)  (FRINGE_OFFSET(j),j=1,3)
      ELSE IF ( CH .EQ. 'M' .or. ch .eq. 'm') THEN
          MODE = MOD(MODE, 2) + 1
      ELSE IF ( CH .EQ. 'I' .or. ch .eq. 'i') THEN
          WRITE(6,*) ' Fringe Information ', ISCN, ISTAR, IBASE
          WRITE(6,*) ' Fringe id mode = ', MODE
          WRITE(6,'(A,5F6.3)') ' fringe visibility = ', V2
          WRITE(6,'(A,5I6)'  ) ' points per fringe = ', COUNT
          WRITE(6,'(A,2I5)'  ) ' NUMBER OF POINTS  = ', NPTS, N2C
          WRITE(6,'(A,2F8.2)') ' ONE COLOR DELAY   = ', DAVG, DRMS
          WRITE(6,'(A,2F8.2)') ' Two color delays  = ', DELAY
          WRITE(6,'(A,2F8.2)') ' Two color D RMS   = ', D2RMS
          WRITE(12,*) ' Fringe Information '
          WRITE(12,*) ' Fringe id mode = ', MODE
          WRITE(12,'(A,5F6.3)') ' fringe visibility = ', V2
          WRITE(12,'(A,5I6)'  ) ' points per fringe = ', COUNT
          WRITE(12,'(A,2I5)'  ) ' NUMBER OF POINTS  = ', NPTS, N2C
          WRITE(12,'(A,2F8.2)') ' ONE COLOR DELAY   = ', DAVG, DRMS
          WRITE(12,'(A,2F8.2)') ' Two color delays  = ', DELAY
          WRITE(12,'(A,2F8.2)') ' Two color D RMS   = ', D2RMS
          GO TO 200
      ELSE IF ( CH .EQ. 'P' .or. ch .eq. 'p') THEN
          PLOT = .NOT. PLOT
      ELSE IF ( CH .EQ. 'H' .or. ch .eq. 'h') THEN
          CALL PGEND
          HCOPY = .TRUE.
          GO TO 80
      ELSE IF ( CH .EQ. 'N' .or. ch .eq. 'n') THEN
          GO TO 900
      ELSE IF ( CH .EQ. 'Q' .or. ch .eq. 'q') THEN
          GO TO 900
      else if ( ch .eq. 'S' .or. ch .eq. 's') then
          stopplot=.true.
          plot=.false.
          goto 900
      else if ( ch .eq. 'O' .or. ch .eq. 'o') then
          other=.not.other
          scanid(8:11)='1000'
          goto 90
      ELSE
          WRITE(6,*) ' BAD COMMAND '
          GO TO 200
      END IF
      GO TO 100
  900 CONTINUE
      CALL PGEND
      endif
      if(scanid(10:11).ne.'00')scanid(1:7)='BAD  ID'
      WRITE(UNIT,1330) ISCN, ISTAR, IBASE, DAVG, DELAY, DRMS, D2RMS,
     $                  V2, COUNT, VLIMIT, NPTS, N2C, MODE, scanid
      IF ( CH .EQ. 'Q' .or. ch .eq. 'q') THEN
          STOP
      END IF
      if(repeat)then  ! In case of a dubious id, we want both results
                 repeat=.false.
                 scanid(8:8)='1'
                 scanid(10:11)='00'
                 goto 90
      endif
      RETURN
C----------------------------------------------------------------------
 1300 FORMAT ( 5F10.3, I5, F6.3 )
 1330 FORMAT ( 2I5, I3, 3F8.2, 2X, 3F5.2, 2X, 5F5.3, 1X, 5I4, F6.3, 3I5,
     $ 1x,a11)
      END