SUBROUTINE IDENT_2(D2C,VSQ,N1,N2,SPACING,PEAK,FRINGE_OFFSET,
     $ IMARK,NPTS,PLOT,scanid,other,disp,phi,d1c)
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)
      REAL *4 D2C(N1,N2),VSQ(*),SPACING(9,2),PEAK(*),FRINGE_OFFSET(3)
      real *4 d2ch(nmax,4)
      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),dispo(4)
      real*4 dx(nmax),dy(nmax),sig(nmax),da,db,siga,sigb,chi2,q
      integer*4 nd2ct,mwt,dn,iother,ifit
      character scanid*11
      LOGICAL PLOT, plotyes, other
      DATA IFRINGE / 0, 1, -1, 2, -2 /
      plotyes=.true.
      iother=0
      ifit=0
      plotyes=plot
      plot=.false.
c---------------------------
c     copy array d2c to d2ch
      do i=1,nmax
       do j=1,4
        d2ch(i,j)=d2c(i,j)
       enddo
      enddo
      do j=1,4
       dispo(j)=disp(j)
      enddo
C-----------------------------------------------------------------------
C     Find principal peak for channel 1
C
 9999 IF ( NPTS .GT. NMAX ) THEN
          WRITE(6,*) ' TOO MUCH DATA ', NPTS, ' IDENT_2 DIES '
          STOP
      ELSE
          WRITE(12,*) ' IDENT_2 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, D2CH(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( D2CH(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) = D2CH(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
      if(ifit.eq.0)then
      mwt=0   ! No uncertainties sig available
      call fit(dx,dy,dn,sig,mwt,da,db,siga,sigb,chi2,q)
      disp(1)=1./db
      endif
C-----------------------------------------------------------------------
C     Find principal peak for channel 2
C
      K = 0
      DO I = 1, NPTS
          IF ( ABS( D2CH(I,1)-P(1) ) .LT. DELTA1 ) THEN
              K = K + 1
              DATA(K) = D2CH(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
      if(ifit.eq.0)then
      dn=0
      do i=1,npts
          if(abs(d2ch(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,db,siga,sigb,chi2,q)
      disp(2)=1./db
      write(12,*)' Fitted new dispersion values (ch 1 & 2): ',disp
      endif
      ifit=ifit+1
C-----------------------------------------------------------------------
C.......................................................................
C     Calculate two color delays for identification purposes only
C
      if(ifit.eq.1)then
      DO I = 1, NMAX
          DO J= 1, 2
              D2CH(I,J) = D1C(I) - DISP(J) * PHI(I,J)
          END DO
      ENDDO
      plot=plotyes
      goto 9999
      endif
C........................................................................
C     Find secondary peak for channel 2
C
      K = 0
      DO I = 1, NPTS
          IF ( ABS( D2CH(I,1)-P(2) ) .LT. DELTA1 ) THEN
              K = K + 1
              DATA(K) = D2CH(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( D2CH(I,1)-P(1) ) .LT. DELTA1 ) THEN
              K = K + 1
              DATA(K) = D2CH(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( D2CH(I,1)-P(2) ) .LT. DELTA1 ) THEN
              K = K + 1
              DATA(K) = D2CH(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, D2CH(1,1), IMARK, 'FIRST  FILTER PAIR')
       CALL PLOT_FRINGE( NPTS, D2CH(1,2), IMARK, 'SECOND FILTER PAIR')
       CALL PLOT_FRINGE( NPTS, D2CH(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(1)-P(3)) .GT. 
     $      ABS(P(2)-FRINGE_OFFSET(1)-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(1)-P(3)) .GT. 
     $      ABS(P(2)-FRINGE_OFFSET(1)-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(1))
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( D2CH(I,1)-P(3) ) .LT. DELTA1 ) .AND.
     $         ( ABS( D2CH(I,2)-P(4) ) .LT. DELTA2 ) .AND.
     $          ( VSQ(I) .GT. VAVG ) ) THEN
              P(1) = P(1) + D2CH(I,1)
              P(2) = P(2) + D2CH(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(1)-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( D2CH(I,1)-P(1) ) .LT. DELTA1 ) .AND.
     $         ( ABS( D2CH(I,2)-P(2) ) .LT. DELTA2 ) ) THEN
              IMARK(I) = 1
          ELSE
              IMARK(I) = 0
          END IF
      END DO
c     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...
      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 ( D2CH(I,1) - P(1) ) .LT. DELTA1 )   THEN
             if(d2ch(i,4).gt.dummy1)then
              imark(i)=5
             else if(d2ch(i,4).lt.dummy3)then
              imark(i)=4
             else
              imark(i)=1
              k=k+1
              dummy2=dummy2+d2ch(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 ( D2CH(I,1) - P(1) ) .LT. DELTA1 )   THEN
             if(d2ch(i,4).gt.dummy1)then
              imark(i)=5
             else if(d2ch(i,4).lt.dummy3)then
              imark(i)=4
             else ! if( ABS( D2CH(I,2)-P(2) ) .LT. DELTA2 ) THEN
              imark(i)=1
              k=k+1
              dummy2=dummy2+d2ch(i,4)
c            else if(d2ch(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(d2ch(i,1)-p(1)).ge.delta1)then
        if(d2ch(i,4).gt.dummy2)then
         imark(i)=3
        else
         imark(i)=2
        endif
       endif
       IMARK(I) = IFRINGE( IMARK(I) )
      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 = D2CH(I,1) - P(1)
          OFFSET2 = D2CH(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, D2CH(1,1), IMARK, 'FIRST  FILTER PAIR')
       CALL PLOT_FRINGE( NPTS, D2CH(1,2), IMARK, 'SECOND FILTER PAIR')
       CALL PLOT_FRINGE( NPTS, D2CH(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)=(d2ch(i,1)-fringe_offset(1)-d2ch(i,2))**2
c    $        +(d2ch(i,1)-fringe_offset(2)-d2ch(i,3))**2
c    $        +(d2ch(i,2)-fringe_offset(3)-d2ch(i,3))**2
c      data(i)=data(i)/3.
       data(i)=sqrt(data(i))
       if(data(i).lt.rmsdd2cm)then
        d2ct=d2ct+d2ch(i,1)
        d2ctw=d2ctw+d2ch(i,1)*vsq(i)
        wsum=wsum+vsq(i)
        nd2ct=nd2ct+1
        d2cta(nd2ct)=d2ch(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(1) 
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!!!'
       iother=iother+1
       if(scanid(8:11).eq.'1000')iother=iother+1
       scanid(9:9)='1'
       other=.true.
       if(iother.eq.1)then
        write(12,*)' *** Try other peak!'
        goto 9999
       endif
      endif
      endif
      do j=1,4
       disp(j)=dispo(j)
      enddo
C-----------------------------------------------------------------------
      RETURN
      END
      SUBROUTINE IDENTIFY ( D2C, N1, N2, SPACING, PEAK, IMARK, NPTS )
C
C     Decide if this data is on fringe -2, -1, 0, +1, or +2.
C     by chosing the fringe that gives the smallest deviation.
C     
      IMPLICIT UNDEFINED (A-Z)
      INTEGER *4 IMARK(*), I, J, K, NPTS, N1, N2, IFRINGE(5)
      REAL *4 D2C(N1,N2), PEAK(*), CHISQ
      REAL *4 NEWCHI, DELTA1, DELTA2, SPACING(9,2)
      DATA IFRINGE / 0, 1, -1, 2, -2 /
      DO I = 1, NPTS
          DO K = 1, 2
              DELTA1 = ABS( D2C(I,K) - PEAK(K) - SPACING(1,1) )
              DELTA2 = ABS( D2C(I,K) - PEAK(K) - SPACING(1,2) )
          END DO
          CHISQ = MAX ( DELTA1, DELTA2 )
          IMARK(I) = 1
          DO J = 2, 5
              NEWCHI = 0.
              DELTA1 = MIN ( 
     $               ABS( D2C(I,1) - PEAK(1) - SPACING(2*J-1,1) ),
     $               ABS( D2C(I,1) - PEAK(1) - SPACING(2*J-2,1) ) )
              DELTA2 = MIN ( 
     $               ABS( D2C(I,2) - PEAK(2) - SPACING(2*J-1,2) ),
     $               ABS( D2C(I,2) - PEAK(2) - SPACING(2*J-2,2) ) )
              NEWCHI = MAX( DELTA1, DELTA2 )
              IF ( NEWCHI .LT. CHISQ ) THEN
                  IMARK(I) = J
                  CHISQ = NEWCHI
              END IF
          END DO
          IMARK(I) = IFRINGE(IMARK(I))
      END DO
      RETURN
      END
      SUBROUTINE AIR ( FSET, DISP, SPACING )
C
C   Calculate the appropriate dispersions for the two color method.
C
      IMPLICIT UNDEFINED (A-Z)
      INTEGER *4 I, J, K, FSET
      REAL *8 WAVE(2), WREF(2), T, P, E, INDEX(2), S1(2), S2(2)
      REAL *8 LAMBDA(4), FR1, FR2
      REAL *4 SPACING(9,2), DISP(4)
      DATA WAVE / .550, .500 /
      DATA WREF / .700, .800 /
      DATA T, P, E / 273., 1013., 0. /
      write(6,'(a,5f8.3)') ' Wavelengths are: ', wave, wref(fset)
      CALL AINDEX ( P, T, E, WREF(FSET), INDEX(2) )
      DO I = 1, 2
          CALL AINDEX ( P, T, E, WAVE(I), INDEX(1) )
          DISP(I) =  INDEX(1) / ( INDEX(1) - INDEX(2) )
          SPACING(1,I) = 0
          DO J = 1, 2
              FR1 = FLOAT(J) * WREF(FSET)
              DO WHILE ( FR1 .GT. 0. )
                  FR1 = FR1 - WAVE(I)
              END DO
              SPACING(4*J-2,I) =  FR1
              SPACING(4*J-1,I) =  FR1 + WAVE(I)
              SPACING(4*J  ,I) = -FR1
              SPACING(4*J+1,I) = -FR1 - WAVE(I)
          END DO
          DO J = 1, 9
              SPACING(J,I) = SPACING(J,I) * DISP(I)
          END DO
      END DO
      if(fset.eq.2)then
       call aindex ( p, t, e, wref(1), index(1) )
       disp(4) = index(1) / ( index(1) - index(2) )
      endif
      CALL AINDEX ( P, T, E, WAVE(1), INDEX(2) )
      CALL AINDEX ( P, T, E, WAVE(2), INDEX(1) )
      DISP(3) =  INDEX(1) / ( INDEX(1) - INDEX(2) )
      RETURN
      END
      SUBROUTINE AINDEX ( P, T, E, WAVE, INDEX )
C
C   Calculates the index of refraction (n-1) of air for pressure P (mb)
C   temperature T (kelvins) partial pressure of water E (mb) and
C   wavelength WAVE (microns).
C
      IMPLICIT UNDEFINED (A-Z)
      REAL *8 P, T, E, INDEX, A, B, K, WAVE
C
      K = 1./(WAVE*WAVE)
      A = 1.E-8*( 2372.434 + 684255.24/(130.-K) + 4549.40/(38.9-K) )
      B = A - 1.E-8*(6487.31 +K*( 58.058 +K*(-0.71150 +K*0.08851 ) ) )
      INDEX = ( A*P - B*E ) / T
      RETURN
      END
      SUBROUTINE PLOT_FRINGE ( N, DATA, IMARK, TITLE )
      IMPLICIT UNDEFINED ( A-Z )
      INTEGER *4 NMAX
      PARAMETER (NMAX=4000)
      INTEGER *4 N, I, IMARK(*)
      REAL *4 DATA(*), XPLOT(NMAX), XMIN, XMAX, YMIN, YMAX
      CHARACTER *(*) TITLE
      IF ( N .GT. NMAX ) THEN
          WRITE(6,*) ' PLOT_FRINGE FAILS '
          WRITE(6,*) ' CHANGE THE SIZE OF NMAX AND RECOMPILE '
      END IF
      XMIN = 0.
      XMAX = FLOAT(N)
      YMIN = DATA(1)
      YMAX = DATA(1)
      DO I = 1, N
          YMIN = MIN ( YMIN, DATA(I) )
          YMAX = MAX ( YMAX, DATA(I) )
          XPLOT(I) = I
      END DO
      CALL PGSCI ( 5 )
      CALL PGSLS ( 1 )
      CALL PGSCH ( 2.0 )
      CALL PGENV ( XMIN, XMAX, YMIN, YMAX, 0, 1 )
      CALL PGMTEXT ( 'T', 1.5, 0.5, 0.5, TITLE )
      DO I = 1, N
          CALL PGSCI ( IMARK(I)+3 )
          CALL PGPOINT ( 1, XPLOT(I), DATA(I) )
      END DO
      RETURN
      END
c*********************************************************************
      SUBROUTINE FIT(X,Y,NDATA,SIG,MWT,A,B,SIGA,SIGB,CHI2,Q)
      real*4 X(NDATA),Y(NDATA),SIG(NDATA)
      SX=0.
      SY=0.
      ST2=0.
      B=0.
      IF(MWT.NE.0) THEN
        SS=0.
        DO 11 I=1,NDATA
          WT=1./(SIG(I)**2)
          SS=SS+WT
          SX=SX+X(I)*WT
          SY=SY+Y(I)*WT
11      CONTINUE
      ELSE
        DO 12 I=1,NDATA
          SX=SX+X(I)
          SY=SY+Y(I)
12      CONTINUE
        SS=FLOAT(NDATA)
      ENDIF
      SXOSS=SX/SS
      IF(MWT.NE.0) THEN
        DO 13 I=1,NDATA
          T=(X(I)-SXOSS)/SIG(I)
          ST2=ST2+T*T
          B=B+T*Y(I)/SIG(I)
13      CONTINUE
      ELSE
        DO 14 I=1,NDATA
          T=X(I)-SXOSS
          ST2=ST2+T*T
          B=B+T*Y(I)
14      CONTINUE
      ENDIF
      B=B/ST2
      A=(SY-SX*B)/SS
      SIGA=SQRT((1.+SX*SX/(SS*ST2))/SS)
      SIGB=SQRT(1./ST2)
      CHI2=0.
      IF(MWT.EQ.0) THEN
        DO 15 I=1,NDATA
          CHI2=CHI2+(Y(I)-A-B*X(I))**2
15      CONTINUE
        Q=1.
        SIGDAT=SQRT(CHI2/(NDATA-2))
        SIGA=SIGA*SIGDAT
        SIGB=SIGB*SIGDAT
      ELSE
        DO 16 I=1,NDATA
          CHI2=CHI2+((Y(I)-A-B*X(I))/SIG(I))**2
16      CONTINUE
c        Q=GAMMQ(0.5*(NDATA-2),0.5*CHI2)
      ENDIF
      RETURN
      END