SUBROUTINE IDENT_3(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 d2cp(nmax,4),pspacing(10,3),davg
      REAL *4 VAVG, CHISQ(5),NEWCHI, DELTA1, DELTA2, P(6),chisq3(3)
      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
      IF ( NPTS .GT. NMAX ) THEN
          WRITE(6,*) ' TOO MUCH DATA ', NPTS, ' IDENT_3 DIES '
          STOP
      ELSE
          WRITE(12,*) ' IDENT_3 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,*)' Fitted dispersion constants will be used'
      else
       write(12,*)' Fitted 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
C
C-----------------------------------------------------------------------
C     Calculate position of central fringe in channel 4
      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
      p(3)=dummy2
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C     Calculate the phase differences and the two color delays,
c     but instead of projecting onto delay axis, proj. on phase axis
C
      DAVG = 0.
      DO I = 1, NPTS
          DAVG = DAVG + D1C(I)
      END DO
      DAVG = DAVG / FLOAT( NPTS )
      DO I = 1, NPTS
          DO J= 1, 2
              d2cp(i,j) = phi(i,j) - 1./disp(j)*(d1c(i)-davg)
          END DO
          d2cp(i,3) = d2c(i,4)     ! shift channel 4 to 3 !
      END DO
c     dummy2=0
c     dummy1=0.07
c     do j=1,3
c      do i=1,npts
c       data(i)=d2cp(i,j)
c      enddo
c      CALL HIST_2(NPTS,DATA,DUMMY2,DUMMY1,'CHAN #  C PEAK',plotyes)
c     enddo
C     Reset spacing values
      do i=1,2
       pspacing(1,i)=0
       pspacing(2,i)=0
       do j=1,9
        pspacing(j+1,i)=spacing(j,i)/disp(i)
       enddo
      enddo
C     The following pspacing values assume 800 nm reference and
C     730 nm broad band tracking channel !!!!
      pspacing(1,3)=0
      pspacing(2,3)=0
      pspacing(3,3)=-0.660
      pspacing(4,3)= 0.070
      pspacing(5,3)= 0.660
      pspacing(6,3)=-0.070
      pspacing(7,3)=-0.590
      pspacing(8,3)= 0.140
      pspacing(9,3)= 0.590
      pspacing(10,3)=-.140
C     Calculate peak values in phase
      p(4)=p(1)          ! save p(1)
      p(1)=(davg-p(1))/disp(1)
      p(2)=(davg-p(2))/disp(2)
      WRITE(12,*) ' PSPACING = '
      WRITE(12,'(10F8.2)') (PSPACING(k,1),k=1,10)
      WRITE(12,'(10F8.2)') (PSPACING(k,2),k=1,10)
      WRITE(12,'(10F8.2)') (PSPACING(k,3),k=1,10) 
      write(12,*)' peaks are: ',p(1),p(2),p(3)
C    
C     Now identify points by calculating the chi^2 for all possible
C     fringes and choosing the identification with the smallest chi^2
C
      do i=1,npts
       do j=1,5       ! each fringe id (0,-1,+1,-2,+2)
        chisq(j)=0
        do k=1,3       ! three channels
         offset1=abs(d2cp(i,k)-p(k)-pspacing(2*j-1,k))
         offset2=abs(d2cp(i,k)-p(k)-pspacing(2*j  ,k))
         if(k.eq.3)then
          offset1=offset1*3.       ! Enhance weight because in this
          offset2=offset2*3.       ! channel, spacings are small
         endif
         chisq(j)=chisq(j)+min(offset1,offset2)**2
        enddo
       enddo
       do j=1,5
        if(chisq(j).eq.
     $     min(chisq(1),chisq(2),chisq(3),chisq(4),chisq(5)))then
           IMARK(I) = IFRINGE(j)
        endif
       enddo
      enddo
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-----------------------------------------------------------------------
chummel  plot histogram of rms delta-2-color delays
      p(1)=p(4)
      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
c     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
      return
      end