SUBROUTINE IDENT_4(D2C,VSQ,N1,N2,SPACING,PEAK,FRINGE_OFFSET,
     $ IMARK,NPTS,PLOT,scanid,other,disp,phi,d1c,ibase,fitdisp)
C
C     Decide if this data is on fringe -2, -1, 0, +1, or +2.
C     The first fringe data is used to separate the even from the odd
C     fringes.  Then the second fringe is used to eliminate double
C     fringe jumps.
C
      IMPLICIT UNDEFINED (A-Z)
      INTEGER *4 NMAX
      PARAMETER ( NMAX = 4000 )
      INTEGER *4 N1, N2, IMARK(*), NPTS, I, J, K, IFRINGE(5),ibase
      REAL *4 D2C(N1,N2),VSQ(*),SPACING(9,2),PEAK(*),FRINGE_OFFSET(3)
      REAL *4 VAVG, CHISQ, NEWCHI, DELTA1, DELTA2, P(6)
      REAL *4 DATA(NMAX), OFFSET1, OFFSET2
      real*4 d2ct,rmsdd2cm,dummy1,dummy2,dummy3,mdian,d2ctw,wsum
      real*4 d2cta(nmax),disp(4),d1c(*),phi(nmax,4)
      real*4 dx(nmax),dy(nmax),sig(nmax),da,db1,db2,siga,sigb,chi2,q
      integer*4 nd2ct,mwt,dn,iother
      character scanid*11
      LOGICAL PLOT, plotyes, other,fitdisp
      DATA IFRINGE / 0, 1, -1, 2, -2 /
      plotyes=.true.
      iother=0
C-----------------------------------------------------------------------
C     Find principal peak for channel 1
C
 9999 IF ( NPTS .GT. NMAX ) THEN
          WRITE(6,*) ' TOO MUCH DATA ', NPTS, ' IDENT_4 DIES '
          STOP
      ELSE
          WRITE(12,*) ' IDENT_4 begins '
      END IF
      DELTA1 = 0.5 * MIN( ABS(SPACING(2,1)), ABS(SPACING(3,1)) )
      DELTA2 = 0.5 * MIN( ABS(SPACING(2,2)), ABS(SPACING(3,2)) )
      P(1) = PEAK(1)
      CALL HIST_2 (NPTS, D2C(1,1), P(1), 0.25*DELTA1, 'CH 1  ALL DATA',
     $             PLOT )
C-----------------------------------------------------------------------
C     Find secondary peak for channel 1
C
      K = 0
      dn=0
      WRITE(12,*) ' FIRST PEAK IS DATA BETWEEN ', P(1)-DELTA1, ' AND ',
     $            P(1)+DELTA1
      DO I = 1, NPTS
          IF  ( ABS( D2C(I,1)-P(1) ) .LT. DELTA1 )   THEN
              IMARK(I) = 1
c             Here prepare for a fit of new dispersion values
              dn=dn+1
              dx(dn)=d1c(i)-peak(1)
              dy(dn)=phi(i,1)
          ELSE
              IMARK(I) = 0
              K = K + 1
              DATA(K) = D2C(I,1)
          END IF
      END DO
      P(2) = PEAK(1)
      CALL HIST_2 ( K, DATA, P(2), 0.25*DELTA1, 'CHAN 1  PEAK 2', PLOT )
c
c     Now fit new dispersion value for channel 1 with largest fringe
c
      mwt=0   ! No uncertainties sig available
      call fit(dx,dy,dn,sig,mwt,da,db1,siga,sigb,chi2,q)
      db1=1./db1
      if(fitdisp)disp(1)=db1
C-----------------------------------------------------------------------
C     Find principal peak for channel 2
C
      K = 0
      DO I = 1, NPTS
          IF ( ABS( D2C(I,1)-P(1) ) .LT. DELTA1 ) THEN
              K = K + 1
              DATA(K) = D2C(I,2)
          END IF
      END DO
      P(3) = PEAK(2)
      CALL HIST_2 ( K, DATA, P(3), 0.25*DELTA2, 'CHAN 2  PEAK 1', PLOT )
c
c     Now calculate dispersion value for channel 2
c
      dn=0
      do i=1,npts
          if(abs(d2c(i,2)-p(3)).lt.delta2)then
               dn=dn+1
               dx(dn)=d1c(i)-peak(1)
               dy(dn)=phi(i,2) 
          endif
      enddo
      mwt=0
      call fit(dx,dy,dn,sig,mwt,da,db2,siga,sigb,chi2,q)
      db2=1./db2
      if(fitdisp)then
       disp(2)=db2
       write(12,*)' New dispersion constants will be used'
      else
       write(12,*)' New dispersion constants will not be used'
      endif
      write(12,*)' Fitted new dispersion values (ch 1 & 2): ',db1,db2
C-----------------------------------------------------------------------
C     Find secondary peak for channel 2
C
      K = 0
      DO I = 1, NPTS
          IF ( ABS( D2C(I,1)-P(2) ) .LT. DELTA1 ) THEN
              K = K + 1
              DATA(K) = D2C(I,2)
          END IF
      END DO
      P(4) = PEAK(2)
      CALL HIST_2 ( K, DATA, P(4), 0.25*DELTA2, 'CHAN 2  PEAK 2', PLOT )
C-----------------------------------------------------------------------
C     Find principal peak for channel 3
C
      K = 0
      DO I = 1, NPTS
          IF ( ABS( D2C(I,1)-P(1) ) .LT. DELTA1 ) THEN
              K = K + 1
              DATA(K) = D2C(I,3)
          END IF
      END DO
      P(5) = PEAK(2)
      CALL HIST_2 ( K, DATA, P(5), 0.25*DELTA2, 'CHAN 3  PEAK 1', PLOT )
C-----------------------------------------------------------------------
C     Find secondary peak for channel 3
C
      K = 0
      DO I = 1, NPTS
          IF ( ABS( D2C(I,1)-P(2) ) .LT. DELTA1 ) THEN
              K = K + 1
              DATA(K) = D2C(I,3)
          END IF
      END DO
      P(6) = PEAK(2)
      CALL HIST_2 ( K, DATA, P(6), 0.25*DELTA2, 'CHAN 3  PEAK 2', PLOT )
C-----------------------------------------------------------------------
      IF ( PLOT ) THEN
          CALL PLOT_FRINGE( NPTS, D2C(1,1), IMARK, 'FIRST  FILTER PAIR')
          CALL PLOT_FRINGE( NPTS, D2C(1,2), IMARK, 'SECOND FILTER PAIR')
          CALL PLOT_FRINGE( NPTS, D2C(1,3), IMARK, 'THIRD  FILTER PAIR')
      END IF
C-----------------------------------------------------------------------
C     Choose the best peak in channel 1.
C
      WRITE(12,*) ' Channel 1 peaks are at ', P(1), P(2)
      WRITE(12,*) ' Channel 2 peaks are at ', P(3), P(4)
      WRITE(12,*) ' Channel 3 peaks are at ', P(5), P(6)
      if(.not.other)then
       IF ( ABS(P(1)-FRINGE_OFFSET(ibase)-P(3)) .GT. 
     $      ABS(P(2)-FRINGE_OFFSET(ibase)-P(4)) ) THEN
           P(1) = P(2)
           P(2) = P(4)
       ELSE
           P(2) = P(3)
       END IF
      else
       IF ( ABS(P(1)-FRINGE_OFFSET(ibase)-P(3)) .GT. 
     $      ABS(P(2)-FRINGE_OFFSET(ibase)-P(4)) ) THEN
           P(2) = P(3)
       ELSE
           P(1) = P(2)
           P(2) = P(4)
       END IF
      endif
      WRITE(12,*) ' Adopted peaks are ', P(1), P(2)
c     dummy1=sign(1.0,p(1)-p(2))
c     dummy2=sign(1.0,fringe_offset(ibase))
c     if(dummy1.ne.dummy2)scanid(8:8)='1'
C-----------------------------------------------------------------------
C     Determine the mean visibility of the good data.
C
      K = 0
      VAVG = 0.
      DO I = 1, NPTS
          IF ( IMARK(I) .EQ. 1 ) THEN
              VAVG = VAVG + VSQ(I)
              K = K + 1
          END IF
      END DO
      IF ( K .GT. 1 ) THEN
          VAVG = VAVG / FLOAT(K)
      END IF
C-----------------------------------------------------------------------
C     Re-calculate the peak positions from the high visibility points.
C
      P(4) = P(2)
      P(3) = P(1)
      P(1) = 0.
      P(2) = 0.
      K = 0
      DO I = 1, NPTS
          IF ( ( ABS( D2C(I,1)-P(3) ) .LT. DELTA1 ) .AND.
     $         ( ABS( D2C(I,2)-P(4) ) .LT. DELTA2 ) .AND.
     $          ( VSQ(I) .GT. VAVG ) ) THEN
              P(1) = P(1) + D2C(I,1)
              P(2) = P(2) + D2C(I,2)
              K = K + 1
          END IF
      END DO
      DO J = 1, 2
          IF ( K .GT. 1 ) THEN
              P(J) = P(J) / FLOAT(K)
          ELSE
              P(J) = P(J+2)
          END IF
      END DO
      WRITE(12,*) ' Edited  peaks are ', P(1), P(2)
c     if(abs(p(1)-fringe_offset(ibase)-p(2)).gt.2)then
c      write(12,*)' *** Warning: peaks differ by more than 2 mu !'
c      scanid(9:9)='1'
c     endif
C-----------------------------------------------------------------------
C     Re-calculate IMARK.  
C
      DO I = 1, NPTS
          IF ( ( ABS( D2C(I,1)-P(1) ) .LT. DELTA1 ) .AND.
     $         ( ABS( D2C(I,2)-P(2) ) .LT. DELTA2 ) ) THEN
              IMARK(I) = 1
          ELSE
              IMARK(I) = 0
          END IF
      END DO
      goto 1235     ! goto 1235 if you want to use the old id. scheme
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Now, here try another identification scheme using the 700/800 phase
C Then skip the next section, which is the old identifier!
C
C First, select all the data in the central (0,+/-2) peak in channel 1...
C ...and use the median value of the ch4 value as the ch4 pos. of cent.fr.
C
      k=0
      DO I = 1, NPTS
       IF ( ABS ( D2C(I,1) - P(1) ) .LT. DELTA1 )   THEN
        k=k+1
        data(k)=d2c(i,4)
       endif
      enddo
      if(k.gt.1)then
       call mdian1(data,k,mdian)
      else if(k.eq.1)then
       mdian=data(1)
      else
       mdian=0
      endif
      dummy2=mdian
      dummy1=0.07
c     CALL HIST_2 ( K, DATA, DUMMY2, DUMMY1, 'CHAN 4  C PEAK', plotyes )
      write(12,*)' median dummy2=',dummy2
c     dummy2= 0.0        ! mean central fringe position in ch 4
      dummy1=dummy2+0.07 ! 2 times this value is the fringe spacing 0 to +/-2
      dummy3=dummy2-0.07
      dummy2=0
      k=0
      DO I = 1, NPTS
       IF ( ABS ( D2C(I,1) - P(1) ) .LT. DELTA1 )   THEN
             if(d2c(i,4).gt.dummy1)then
              imark(i)=5
             else if(d2c(i,4).lt.dummy3)then
              imark(i)=4
             else
              imark(i)=1
              k=k+1
              dummy2=dummy2+d2c(i,4)
             endif
       END IF
      END DO
      if(k.ne.0)then
       dummy2=dummy2/float(k)
      else
       dummy2=0
       write(12,*)' *** Warning: no central fringe points found!'
      endif
C
C     Go through the last step again, since the green fringe might be
C     offset from zero
      dummy1=dummy2+0.07
      dummy3=dummy2-0.07
      dummy2=0
      k=0
      DO I = 1, NPTS
       IF ( ABS ( D2C(I,1) - P(1) ) .LT. DELTA1 )   THEN
             if(d2c(i,4).gt.dummy1)then
              imark(i)=5
             else if(d2c(i,4).lt.dummy3)then
              imark(i)=4
             else ! if( ABS( D2C(I,2)-P(2) ) .LT. DELTA2 ) THEN
              imark(i)=1
              k=k+1
              dummy2=dummy2+d2c(i,4)
c            else if(d2c(i,4).gt.(dummy1+dummy3)/2)then
c             imark(i)=1
c            else
c             imark(i)=2
             endif
       END IF
      END DO
      if(k.ne.0)then
       dummy2=dummy2/float(k)
      else
       dummy2=0
       write(12,*)' *** Warning: no central fringe points found!'
      endif
C Now you have identified 0 and +/-2; now look at all other points
C to determine +/-1. Their half way distance in d2c should be at dummy2
      do i=1,npts
       if(abs(d2c(i,1)-p(1)).ge.delta1)then
        if(d2c(i,4).gt.dummy2)then
         imark(i)=3
        else
         imark(i)=2
        endif
       endif
       IMARK(I) = IFRINGE( IMARK(I) )
       if(imark(i).eq.0.and.abs(d2c(i,2)-p(2)).gt.delta2)imark(i)=4
      enddo
      goto 1234
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C-----------------------------------------------------------------------
C    P(1), P(2) are the positions of the central fringe in each channel
C       Fringes +1 and -1 are easilly separated from the central fringe,
C       but +2 and -2 are not easilly separated from either.
C    The fringe identification proceedure first separates fringe 0
C       from +1/-1.
C    Then, it separates out the groups 0, +2, -2 and +1, -1, +2, -2.
C
 1235 DO I = 1, NPTS
C
C     The position rel to central fringe.
C     Is this closer to 0 or +/-1?
C
          OFFSET1 = D2C(I,1) - P(1)
          OFFSET2 = D2C(I,2) - P(2)
          IF ( ( ABS(OFFSET1) .LT. DELTA1 ) .AND.
     $         ( ABS(OFFSET2) .LT. DELTA2 ) ) THEN
C
C     Closer to 0.
C
              IMARK(I) = 1
              CHISQ = OFFSET1**2 + OFFSET2**2
              DO J = 6, 7
                  DO K = 6, 7
                      NEWCHI = (OFFSET1-SPACING(J,1))**2
     $                       + (OFFSET2-SPACING(K,2))**2
                      IF ( NEWCHI .LT. CHISQ ) THEN
                          CHISQ = NEWCHI
                          IMARK(I) = 4
                      END IF
                  END DO
              END DO
              DO J = 8, 9
                  DO K = 8, 9
                      NEWCHI = (OFFSET1-SPACING(J,1))**2
     $                       + (OFFSET2-SPACING(K,2))**2
                      IF ( NEWCHI .LT. CHISQ ) THEN
                          CHISQ = NEWCHI
                          IMARK(I) = 5
                      END IF
                  END DO
              END DO
          ELSE
C
C     Closer to +1/-1
C
              CHISQ = 1.E19
              DO J = 2, 3
                  DO K = 2, 3
                      NEWCHI = (OFFSET1-SPACING(J,1))**2
     $                       + (OFFSET2-SPACING(K,2))**2
                      IF ( NEWCHI .LT. CHISQ ) THEN
                          CHISQ = NEWCHI
                          IMARK(I) = 2
                      END IF
                  END DO
              END DO
              DO J = 4, 5
                  DO K = 4, 5
                      NEWCHI = (OFFSET1-SPACING(J,1))**2
     $                       + (OFFSET2-SPACING(K,2))**2
                      IF ( NEWCHI .LT. CHISQ ) THEN
                          CHISQ = NEWCHI
                          IMARK(I) = 3
                      END IF
                  END DO
              END DO
              DO J = 6, 7
                  DO K = 6, 7
                      NEWCHI = (OFFSET1-SPACING(J,1))**2
     $                       + (OFFSET2-SPACING(K,2))**2
                      IF ( NEWCHI .LT. CHISQ ) THEN
                          CHISQ = NEWCHI
                          IMARK(I) = 4
                      END IF
                  END DO
              END DO
              DO J = 8, 9
                  DO K = 8, 9
                      NEWCHI = (OFFSET1-SPACING(J,1))**2
     $                       + (OFFSET2-SPACING(K,2))**2
                      IF ( NEWCHI .LT. CHISQ ) THEN
                          CHISQ = NEWCHI
                          IMARK(I) = 5
                      END IF
                  END DO
              END DO
          END IF
C         write(12,111) i, offset1, offset2, imark(i), ifringe(imark(i))
          IMARK(I) = IFRINGE( IMARK(I) )
          IF ( PLOT ) THEN
          END IF
      END DO
  111 FORMAT ( I5, 2F10.1, 2X, 2I5 )
C-----------------------------------------------------------------------
 1234 continue    ! Continue here in new identification scheme
      IF ( PLOT ) THEN
          CALL PLOT_FRINGE( NPTS, D2C(1,1), IMARK, 'FIRST  FILTER PAIR')
          CALL PLOT_FRINGE( NPTS, D2C(1,2), IMARK, 'SECOND FILTER PAIR')
          CALL PLOT_FRINGE( NPTS, D2C(1,3), IMARK, 'THIRD  FILTER PAIR')
      END IF
C-----------------------------------------------------------------------
chummel  plot histogram of rms delta-2-color delays
      rmsdd2cm=2
      nd2ct=0
      d2ct=0
      d2ctw=0
      wsum=0
      do i=1,npts
       data(i)=(d2c(i,1)-fringe_offset(ibase)-d2c(i,2))**2
c    $        +(d2c(i,1)-fringe_offset(2)-d2c(i,3))**2
c    $        +(d2c(i,2)-fringe_offset(3)-d2c(i,3))**2
c      data(i)=data(i)/3.
       data(i)=sqrt(data(i))
       if(data(i).lt.rmsdd2cm)then
        d2ct=d2ct+d2c(i,1)
        d2ctw=d2ctw+d2c(i,1)*vsq(i)
        wsum=wsum+vsq(i)
        nd2ct=nd2ct+1
        d2cta(nd2ct)=d2c(i,1)
       endif
      enddo
c     dummy1=rmsdd2cm
c     dummy2=1
c     call hist_2(npts,data,dummy1,dummy2,'RMS dD2C all d',plotyes)
      if(nd2ct.ne.0)then
       d2ct=d2ct/float(nd2ct)
      else
       d2ct=0
      endif
      if(wsum.ne.0)then
       d2ctw=d2ctw/wsum
      else
       d2ctw=0
      endif
      if(nd2ct.gt.1)then
       call mdian1(d2cta,nd2ct,mdian)
      else if(nd2ct.eq.1)then
       mdian=d2cta(1)
      else 
       mdian=0
       write(12,*)' *** Warning: nd2ct = 0 !'
      endif
      write(12,*)' d2ct=',d2ct,' nd2ct=',nd2ct,' rmsdd2cm=',rmsdd2cm
      write(12,*)' using fringe offsets 1-2=',fringe_offset(ibase) 
c    $                                ,'1-3=',fringe_offset(2),
c    $                                 '2-3=',fringe_offset(3)
      write(12,*)' Median value of d2ct is: ',mdian
      write(12,*)' Weighted mean value  is: ',d2ctw
      if(d2ct.ne.0)then
      if(abs(mdian-p(1)).gt.delta1/2.0)then
       write(12,*)' *** Warning: Hummels additional test failed!!!'
       scanid(9:9)='1'
C      iother=iother+1
C      other=.not.other
C      if(iother.eq.1)then
C       write(12,*)' *** Try other peak!'
C       goto 9999
C      endif
      endif
      endif
C-----------------------------------------------------------------------
      RETURN
      END