SUBROUTINE GETADATA ( NCOH, TAVG, DRES, PHASE, VSQ, NAVG,
     $                      LASER0, PHIRMS, PHIAC, STATUS )
C
C      David Mozurkewich August 22, 1991
C
C      Get the 4 ms data for one coherent average.
C      average time and laser
C      Do astrometry and calculate the sideral time
C      Do coherent and incoherent averages
C      Return the phase in microns and the visibility amp squared
C
C      If end of star is reached, return no data.      
C       
C     ARGUEMENTS:
C       NCOH     number of points in coherent average
C       TAVG     averaged time
C       DRES     averaged residual one color delay ( laser position )
C       PHASE(4) averaged phase for each channel.
C       VSQ      averaged channel 1 visibility amplitude squared.
C       NAVG(4)  averaged number of photons in one coherent average
C       LASER0   averaged laser delay
C       PHIRMS   phase rms in radians
C       PHIAC    first lag of phase autocorrelation
C       NEWSTAR  True if a new star started before enough data
C                was returned.
C       MOREDATA False if there is no more data.
C
C       If NEWSTAR is true or MOREDATA is false, the data returned
C                  should not be used.
C
      IMPLICIT UNDEFINED (A-Z)
      SAVE
      REAL *8 PI, RAD, LON, LAT
      PARAMETER ( PI = 3.14159265358979D0 )
      PARAMETER ( RAD = 180.D0 / PI  )
      PARAMETER ( LON = 118.059167D0 )
      PARAMETER ( LAT =  34.2166944D0)
      INTEGER *4 CNT
      PARAMETER ( CNT = 1000 )
      INCLUDE 'AVER.INC'
      LOGICAL MOREDATA, NEWSTAR, LOSTLOCK
      INTEGER *4 NCOH, I, J, K, ICHAN, NJUMP, NPTS
      INTEGER *4 SPECT, LASTSTAR, LASTBASE, NTEMP
      INTEGER *4 ITEST
      REAL *4 PMRA4, PMDEC4, PAR4, RV4, MAG4, BV, ANGDIA
      REAL *4 PHASE(4), VSQ, DRES, XCOH(4), YCOH(4)
      REAL *8 PHIRMS, PHIAC, PHIAVG, CP, SP
      REAL *8 VNUM, NAVG(4), TAVG, LASER0, COEF(4)
      REAL *8 DELAY0(CNT), X0(4,CNT), Y0(4,CNT), N0(4,CNT), TIME(CNT)
      REAL *8 PHI(CNT), XTEMP, UT1, RAA, DECA, GAST, JD, TJD, DJUL
      CHARACTER *(15) STATUS
      EXTERNAL JD
      DATA NEWSTAR / .TRUE. /
      DATA LASTSTAR, LASTBASE, NTEMP /  0, 0, 0  /
      DATA ITEST / 0 /
C-----------------------------------------------------------------------
      STATUS = ' '
      TAVG = 0.
      DRES = 0.
      LASER0 = 0.
C-----------------------------------------------------------------------
C  Check input parameters
C
      IF ( IREC .GT. 220 ) ITEST = 1
      IF ( NCOH.LE.0 ) THEN
          WRITE(6,*) ' NCOH  < OR = TO ZERO '
          STATUS = 'ARG ERROR'
          GO TO 500
      ELSE IF ( NCOH .GT. CNT ) THEN
          WRITE(6,*) ' COHERENT INTEGRATION IS TOO LONG '
          WRITE(6,*) ' INCREASE THE SIZE OF CNT AND RECOMPILE '
          WRITE(6,*) ' GETADATA BOMBS *********************** '
          STATUS = 'ARG ERROR'
          GO TO 500
      END IF
C-----------------------------------------------------------------------
C   If this is a new star, make sure there is data for it.
C
      IF ( IREC .EQ. 0 ) NEWSTAR = .TRUE.
      DO WHILE ( NEWSTAR ) 
          CALL RDATA ( NEWSTAR, MOREDATA, LOSTLOCK )
          IF ( .NOT. MOREDATA ) THEN
              STATUS = 'NO MORE DATA'
              GO TO 500
          END IF
      END DO
C-----------------------------------------------------------------------
C  Get NCOH data samples.
C
  100 CONTINUE
          DO I = 1, NCOH
              DO ICHAN = 1, 4
                  X0(ICHAN,I) = DBLE( X(ICHAN) )
                  Y0(ICHAN,I) = DBLE( Y(ICHAN) )
                  N0(ICHAN,I) = DBLE( N(ICHAN) )
              END DO
              DELAY0(I) = LASER
              TIME (I) = MTIME
              CALL RDATA ( NEWSTAR, MOREDATA, LOSTLOCK )
              IF ( .NOT. MOREDATA ) THEN
                  STATUS = 'NO MORE DATA'
                  GO TO 500
              ELSE IF ( NEWSTAR ) THEN
                  STATUS = 'END OF STAR'
                  GO TO 500
              ELSE IF ( LOSTLOCK ) THEN
                  GO TO 100
              END IF
          END DO
C-----------------------------------------------------------------------
C   Update the current star/baseline information
C
      IF ( LASTSTAR .NE. ISTAR ) THEN
          CALL CATALOG( ISTAR, NAME, RAM, DECM, PMRA4, PMDEC4,
     $                    PAR4, RV4, MAG4, BV, SPECT, ANGDIA )
          write(6,'(a,i5,1x,a/)') ' Catalog called: ', ISTAR, NAME
          PMRA  = PMRA4
          PMDEC = PMDEC4
          PAR   = PAR4
          RV    = RV4
          LASTSTAR = ISTAR
      END IF
      IF ( LASTBASE .NE. IBASE ) THEN
          BX = OBX (IBASE)
          BY = OBY (IBASE)
          BZ = OBZ (IBASE)
          CON= OCON(IBASE)
          LASTBASE = IBASE
      END IF
C-----------------------------------------------------------------------
C     Make sure we know the filter wavelengths.  Determine the
C     conversions from phase to delay referenced to the tracking channel.
C     For channel 1, it is simply the conversion from phase(radians) to
C     delay(microns), using the wavelength (nm).
C
      IF ( ( FILTER(1) .EQ. -1 ) .OR.
     $     ( FILTER(2) .EQ. -1 ) .OR.
     $     ( FILTER(3) .EQ. -1 ) .OR.
     $     ( FILTER(4) .EQ. -1 ) ) THEN
          STATUS = 'NO FILTER ID'
          GO TO 500
      ELSE
          DO I = 1, 4
              WAVE(I) = WAVES(FILTER(I))
          END DO
          IF ( WAVE(1)*WAVE(2)*WAVE(3)*WAVE(4) .EQ. 0. ) THEN
              WRITE(6,*) ' BAD FILTER ID '
              DO I = 1, 4
                  WRITE(6,*) ' CHANNEL ', I, ' WAVE = ', WAVE(I)
              END DO
              STATUS = 'BAD FILTER ID'
              GO TO 500
          ELSE
              COEF(1) = 1.
              COEF(2) = WAVE(1) / WAVE(2)
              COEF(3) = WAVE(1) / WAVE(3)
              COEF(4) = WAVE(1) / WAVE(4)
          END IF
      END IF
C-----------------------------------------------------------------------
C     Determine the mean time and laser.
C
      TAVG = 0.
      LASER0 = 0.
      DO I = 2, NCOH
          TAVG = TAVG + TIME(I) - TIME(1)
          LASER0 = LASER0 + DELAY0(I) - DELAY0(1)
      END DO
      TAVG = TIME(1)     + TAVG / DBLE(NCOH)
      LASER0 = DELAY0(1) + LASER0 / DBLE(NCOH)
C-----------------------------------------------------------------------
C     Do the astrometry to determine the residual delay.
C
      DJUL = JD ( IDAY, IMON, IYEAR )
      UT1  = TAVG / 86400.D0
      TJD  = DJUL + UT1 + DBLE(TDTOFF)/86400.D0
      CALL APSTAR (TJD, 3, RAM, DECM, PMRA, PMDEC, PAR,
     $                         RV, RAA, DECA)
      CALL SIDTIM( DJUL, UT1, 1, GAST )
      CALL DIURN ( TJD, GAST, RAA, DECA, RA, DEC )
      HA = 15.D0 *( GAST - RA ) - LON
      IF ( HA .LT. -180.D0 ) HA = HA + 360.D0
      DRES = LASER0 - CON - BX*COS(DEC/RAD)*COS(HA/RAD)
     $                    + BY*COS(DEC/RAD)*SIN(HA/RAD)
     $                    - BZ*SIN(DEC/RAD)
C-----------------------------------------------------------------------
C    Unwrap the red phase.  Calculate the phase rms and the first lag
C    of the phase autocorrelation.  rotate X and Y to 0 red phase.
C    If there is no data, use the mean phase ( pi/4).
C
      NJUMP = 0
      IF ( (X0(1,1) .EQ. 0 ) .AND. ( Y0(1,1) .EQ. 0 )) THEN
          PHI(1) = PI/4.D0
          WRITE(6,'(A,/)') ' X=Y=0, CHANNEL 1 '
      ELSE
          PHI(1) = ATAN2 ( X0(1,1), Y0(1,1) ) + PI/4.D0
      END IF
      DO I = 2, NCOH
          IF ( (X0(1,I) .EQ. 0 ) .AND. ( Y0(1,I) .EQ. 0 )) THEN
              PHI(I) = PI/4.D0
              WRITE(6,'(A,/)') ' X=Y=0, CHANNEL 1 '
          ELSE
              PHI(I) = ATAN2( X0(1,I), Y0(1,I) ) 
     $               + PI/4.D0 + 2.*PI*DBLE(NJUMP)
          END IF
          IF ( (PHI(I)-PHI(I-1)) .GT. PI ) THEN
              NJUMP = NJUMP - 1
              PHI(I) = PHI(I) - 2.*PI
          ELSE IF ( (PHI(I)-PHI(I-1)) .LT. -PI ) THEN
              NJUMP = NJUMP + 1
              PHI(I) = PHI(I) + 2.*PI
          END IF
      END DO
C-----------------------------------------------------------------------
C   Calculate the phase rms and autocorrelation to estimate the seeing.
C
      PHIAVG  = 0.
      PHIRMS  = 0.
      PHIAC   = 0.
      DO I = 1, NCOH
          PHIAVG  = PHIAVG + PHI(I)
          PHIRMS  = PHIRMS + PHI(I)*PHI(I)
          PHIAC   = PHIAC  + PHI(I)*PHI(I-1)
      END DO
      PHIAVG = PHIAVG / FLOAT(NCOH)
      PHIRMS = PHIRMS/FLOAT(NCOH) - PHIAVG*PHIAVG
      IF ( PHIRMS .LE. 0. ) THEN
          WRITE(6,'(A,/)') ' *** PHASE RMS = 0 !!! '
          PHIRMS = 0.
      ELSE
          PHIRMS = SQRT ( PHIRMS )
      END IF
      PHIAC  = PHIAC / FLOAT(NCOH-1) - PHIAVG*PHIAVG
      IF ( PHIAC .GT. 0. ) THEN
          PHIAC = SQRT ( PHIAC )
      ELSE
          WRITE(6,'(A,/)') ' *** PHASE AC = 0 !!! '
          PHIAC = 0.0
      END IF
C-----------------------------------------------------------------------
C  Rotate X and Y by the channel 1 phase, in order to extend the 
C  coherent integration.  Channel 1 rotated phase is zero, but do
C  it anywhy to get the correct visibility.
C
      DO I = 1, NCOH
          DO ICHAN = 1, 4
              CP = COS( PHI(I) * COEF(ICHAN) )
              SP = SIN( PHI(I) * COEF(ICHAN) )
              XTEMP = X0(ICHAN,I)*CP - Y0(ICHAN,I)*SP
              Y0(ICHAN,I) = Y0(ICHAN,I)*CP + X0(ICHAN,I)*SP
              X0(ICHAN,I) = XTEMP
          END DO
      END DO
C-----------------------------------------------------------------------
C   Do the coherent integration. 
C
      DO ICHAN = 1, 4
          XCOH(ICHAN) = 0.
          YCOH(ICHAN) = 0.
          NAVG(ICHAN) = 0.
      END DO
      DO K = 1, NCOH
          DO ICHAN = 1, 4
              XCOH(ICHAN) = XCOH(ICHAN) + X0(ICHAN,K)
              YCOH(ICHAN) = YCOH(ICHAN) + Y0(ICHAN,K)
              NAVG(ICHAN) = NAVG(ICHAN) + N0(ICHAN,K)
          END DO
      END DO
C
C  Flip the sign on channels 2 and 4, since they come off the other
C  side of the beam splitter
C
      XCOH(2) = -XCOH(2)
      YCOH(2) = -YCOH(2)
      XCOH(4) = -XCOH(4)
      YCOH(4) = -YCOH(4)
      DO ICHAN = 1, 4
          PHASE(ICHAN) = ATAN2( XCOH(ICHAN), YCOH(ICHAN) )
      END DO
      VNUM = XCOH(1)**2 + YCOH(1)**2 - NAVG(1)
      IF ( NAVG(1) .GT. 0. ) THEN
          VSQ = 0.5*PI*PI*VNUM / ( NAVG(1)**2 )
      ELSE
          WRITE(6,'(A,/)') ' NO PHOTONS THIS COHERENT INTEGRATION '
          VSQ = 0.
      END IF
      PHASE(1) = PHIAVG
C-----------------------------------------------------------------------
C
C     There is a sign error in the definition of the phase above.
C     it should be atan2 ( y, x ).  To fix this, I added a - sign
C     to the phase in the next line
C-----------------------------------------------------------------------
C   Convert the phase, in radians, to microns.
C   The wavelength is in nanometers.
C
c     ichan=1
c     PHASE(ICHAN) = -(WAVE(ICHAN)/1000.) * (PHASE(ICHAN)+pi)/ (2.*PI)
      DO ICHAN = 1, 4
          PHASE(ICHAN) = -(WAVE(ICHAN)/1000.) * PHASE(ICHAN) / ( 2.*PI )
      END DO
C-----------------------------------------------------------------------
  500 CONTINUE
      RETURN
      END