SUBROUTINE DEMOD(X,N,F) 
C
C     PURPOSE--THIS SUBROUTINE PERFORMS A COMPLEX DEMODULATION
C              ON THE DATA IN THE INPUT VECTOR X
C              AT THE INPUT DEMODULATION FREQUENCY = F.
C              THE COMPLEX DEMODULATION CONSISTS OF THE FOLLOWING--
C              1) AN AMPLITUDE VERSUS TIME PLOT;
C              2) A PHASE VERSUS TIME PLOT;
C              3) AN UPDATED DEMODULATION FREQUENCY ESTIMATE
C                 TO ASSIST THE ANALYST IN DETERMINING A
C                 MORE APPROPRIATE FREQUENCY AT WHICH
C                 TO DEMODULATE IN CASE THE SPECIFIED
C                 INPUT DEMODULATION FREQUENCY F
C                 DOES NOT FLATTEN SUFFICIENTLY THE
C                 PHASE PLOT. 
C
C              THE ALLOWABLE RANGE OF THE INPUT DEMODULATION
C              FREQUENCY F IS 0.0 TO 0.5 (EXCLUSIVELY).
C              THE INPUT DEMODULATION FREQUENCY F IS MEASURED  OF
C              IN UNITS OF CYCLES PER 'DATA POINT' OR,
C              MORE PRECISELY, IN CYCLES PER UNIT TIME WHERE
C              'UNIT TIME' IS DEFINED AS THE
C              ELAPSED TIME BETWEEN ADJACENT OBSERVATIONS.
C
C     INPUT ARGUMENTS--X      = THE SINGLE PRECISION VECTOR OF
C                               (UNSORTED) OBSERVATIONS.
C                      N      = THE INTEGER NUMBER OF OBSERVATIONS
C                               IN THE VECTOR X.
C                      F      = THE SINGLE PRECISION
C                               DEMODULATION FREQUENCY.
C                               F IS IN UNITS OF CYCLES PER DATA POINT.
C                               F IS BETWEEN 0.0 AND 0.5 (EXCLUSIVELY).
C     OUTPUT--2 PAGES OF AUTOMATIC PRINTOUT--
C             1) AN AMPLITUDE PLOT;
C             2) A PHASE PLOT; AND
C             3) AN UPDATED DEMODULATION FREQUENCY ESTIMATE.
C     PRINTING--YES.
C     RESTRICTIONS--THE MAXIMUM ALLOWABLE VALUE OF N
C                   FOR THIS SUBROUTINE IS 5000.
C                 --THE SAMPLE SIZE N MUST BE GREATER
C                   THAN OR EQUAL TO 3. 
C                 --THE INPUT FREQUENCY F MUST BE 
C                   GREATER THAN OR EQUAL TO 2/(N-2).
C                 --THE INPUT FREQUENCY F MUST BE 
C                   SMALLER THAN 0.5.
C     OTHER DATAPAC   SUBROUTINES NEEDED--PLOTX.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, SIN, COS, ATAN.
C     MODE OF INTERNAL OPERATIONS--SINGLE PRECISION.
C     LANGUAGE--ANSI FORTRAN. 
C     COMMENT--IN ORDER THAT THE RESULTS OF THE COMPLEX DEMODULATION
C              BE VALID AND PROPERLY INTERPRETED, THE INPUT DATA
C              IN X SHOULD BE EQUI-SPACED IN TIME 
C              (OR WHATEVER VARIABLE CORRESPONDS TO TIME).
C            --IF THE INPUT OBSERVATIONS IN X ARE CONSIDERED
C              TO HAVE BEEN COLLECTED 1 SECOND APART IN TIME,
C              THEN THE DEMODULATION FREQUENCY F
C              WOULD BE IN UNITS OF HERTZ
C              (= CYCLES PER SECOND).
C            --A FREQUENCY OF 0.0 CORRESPONDS TO A CYCLE
C              IN THE DATA OF INFINITE (= 1/(0.0))
C              LENGTH OR PERIOD.
C              A FREQUENCY OF 0.5 CORRESPONDS TO A CYCLE
C              IN THE DATA OF LENGTH = 1/(0.5) = 2 DATA POINTS.
C            --IN EXAMINING THE AMPLITUDE AND PHASE PLOTS,
C              ATTENTION SHOULD BE PAID NOT ONLY TO THE
C              STRUCTURE OF THE PHASE PLOT
C              (NEAR-ZERO SLOPE VERSUS NON-ZERO SLOPE)
C              BUT ALSO TO THE RANGE
C              OF VALUES ON THE VERTICAL AXIS.
C              A PLOT WITH MUCH STRUCTURE BUT
C              WITH A SMALL RANGE ON THE VERTICAL AXIS
C              IS USUALLY MORE INDICATIVE OF A
C              DEFINITE CYCLIC COMPONENT AT THE
C              SPECIFIED INPUT DEMODULATION FREQUENCY,
C              THAN IS A PLOT WITH LESS STRUCTURE BUT
C              A WIDER RANGE ON THE VERTICAL AXIS.
C            --INTERNAL TO THIS SUBROUTINE, 2 MOVING
C              AVERAGES ARE APPLIED, EACH OF LENGTH 1/F.
C              HENCE THE AMPLITUDE AND PHASE PLOTS
C              HAVE N - 2/F VALUES
C              (RATHER THAN N VALUES) ALONG THE
C              HORIZONTAL (TIME) AXIS.
C              IN ORDER THAT THE AMPLITUDE AND PHASE
C              PLOTS BE NON-EMPTY, AN INPUT
C              REQUIREMENT ON F FOR THIS SUBROUTINE
C              IS THAT THE SAMPLE SIZE N
C              AND THE DEMODULATION FREQUENCY F
C              MUST BE SUCH THAT
C              N - 2/F BE GREATER THAN ZERO.
C              FURTHER, SINCE A PLOT WITH BUT
C              1 POINT IS MEANINGLESS
C              AND OUGHT ALSO BE EXCLUDED,
C              THE REQUIREMENT IS EXTENDED
C              SO THAT N - 2/F MUST BE GREATER THAN 1.
C     REFERENCES--GRANGER AND HATANAKA, PAGES 170 TO 189,
C                 ESPECIALLY PAGES 173, 177, AND 182.
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING LABORATORY (205.03)
C                 NATIONAL BUREAU OF STANDARDS
C                 WASHINGTON, D. C. 20234
C                 PHONE--301-921-2315
C     ORIGINAL VERSION--NOVEMBER  1972. 
C     UPDATED         --NOVEMBER  1975. 
C     UPDATED         --FEBRUARY  1976. 
C
C---------------------------------------------------------------------
C
      DIMENSION X(1)
      DIMENSION Y1(5000),Y2(5000),Z(5000)
      COMMON /BLOCK2/ WS(15000)
      EQUIVALENCE (Y1(1),WS(1)),(Y2(1),WS(5001)),(Z(1),WS(10001))
      DATA PI/3.141592653/
C
      IPR=6
      ILOWER=3
      IUPPER=5000
      AN=N
      FMIN=2.0/(AN-2.0)
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(N.LT.ILOWER.OR.N.GT.IUPPER)GOTO50
      IF(F.LE.FMIN.OR.F.GE.0.5)GOTO60
      HOLD=X(1)
      DO65I=2,N
      IF(X(I).NE.HOLD)GOTO90
   65 CONTINUE
      WRITE(IPR, 9)HOLD
      RETURN
   50 WRITE(IPR,17)ILOWER,IUPPER
      WRITE(IPR,47)N
      RETURN
   60 WRITE(IPR,27)FMIN
      WRITE(IPR,46)F
      WRITE(IPR,28)FMIN,N
      RETURN
   90 CONTINUE
    9 FORMAT(1H ,109H***** NON-FATAL DIAGNOSTIC--THE FIRST  INPUT ARGUME
     1NT (A VECTOR) TO THE DEMOD  SUBROUTINE HAS ALL ELEMENTS = ,E15.8,6
     1H *****)
   17 FORMAT(1H , 96H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 DEMOD  SUBROUTINE IS OUTSIDE THE ALLOWABLE (,I6,1H,,I6,16H) INTER
     1VAL *****)
   27 FORMAT(1H , 96H***** FATAL ERROR--THE THIRD INPUT ARGUMENT TO THE
     1 DEMOD  SUBROUTINE IS OUTSIDE THE ALLOWABLE (,F10.8,6H,0.5) ,
     1 14HINTERVAL *****)
   28 FORMAT(1H ,42H                   THE ABOVE LOWER LIMIT (,F10.8, 
     1 46H) = 2/(N-2) WHERE N = THE INPUT SAMPLE SIZE = ,I8)
   46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****)
   47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8   ,6H *****)
C
C-----START POINT-----------------------------------------------------
C
C     FORM THE COSINE AND SINE SERIES
C
      DO100I=1,N
      AI=I
      Y1(I)=X(I)*COS(6.2831853*F*AI)
      Y2(I)=X(I)*SIN(6.2831853*F*AI)
  100 CONTINUE
C
C     DEFINE THE LENGTH OF THE 2 MOVING AVERAGES
C
      LENMA1=1.0/F
      LENMA2=1.0/F
      ALEN1=LENMA1
      ALEN2=LENMA2
      IMAX1=N-LENMA1
      IMAX2=IMAX1-LENMA2
C
C     FORM THE FIRST MOVING AVERAGE FOR THE COSINE SERIES
C
      DO200I=1,IMAX1
      ISTART=I+1
      IEND=I+LENMA1-1
      IENDP1=I+LENMA1
      SUM=0.0
      DO300J=ISTART,IEND
      SUM=SUM+Y1(J) 
  300 CONTINUE
      SUM=SUM+Y1(I)/2.0+Y1(IENDP1)/2.0
      Z(I)=SUM/ALEN1
  200 CONTINUE
C
C     FORM THE SECOND MOVING AVERAGE FOR THE COSINE SERIES
C
      DO400I=1,IMAX2
      ISTART=I+1
      IEND=I+LENMA2-1
      IENDP1=I+LENMA2
      SUM=0.0
      DO500J=ISTART,IEND
      SUM=SUM+Z(J)
  500 CONTINUE
      SUM=SUM+Z(I)/2.0+Z(IENDP1)/2.0
      Y1(I)=SUM/ALEN2
  400 CONTINUE
C
C     FORM THE FIRST MOVING AVERAGE FOR THE SINE SERIES
C
      DO800I=1,IMAX1
      ISTART=I+1
      IEND=I+LENMA1-1
      IENDP1=I+LENMA1
      SUM=0.0
      DO900J=ISTART,IEND
      SUM=SUM+Y2(J) 
  900 CONTINUE
      SUM=SUM+Y2(I)/2.0+Y2(IENDP1)/2.0
      Z(I)=SUM/ALEN1
  800 CONTINUE
C
C     FORM THE SECOND MOVING AVERAGE FOR THE SINE SERIES
C
      DO1000I=1,IMAX2
      ISTART=I+1
      IEND=I+LENMA1-1
      IENDP1=I+LENMA1
      SUM=0.0
      DO1100J=ISTART,IEND
      SUM=SUM+Z(J)
 1100 CONTINUE
      SUM=SUM+Z(I)/2.0+Z(IENDP1)/2.0
      Y2(I)=SUM/ALEN2
 1000 CONTINUE
C
C
C     FORM THE AMPLITUDES AND PLOT THEM 
C
      DO1500I=1,IMAX2
      Z(I)=2.0*SQRT(Y1(I)*Y1(I)+Y2(I)*Y2(I))
 1500 CONTINUE
      CALL PLOTX(Z,IMAX2)
      WRITE(IPR,2005)F
C
C     COMPUTE THE DIFFERENCE BETWEEN THE MAX AND MIN AMPLITUDES AND WRITE IT OUT
C
      ZMIN=Z(1)
      ZMAX=Z(1)
      DO1600I=1,IMAX2
      IF(Z(I).LT.ZMIN)ZMIN=Z(I)
      IF(Z(I).GT.ZMAX)ZMAX=Z(I)
 1600 CONTINUE
      RANGE=ZMAX-ZMIN
      WRITE(IPR,2010)ZMIN,ZMAX,RANGE
C
C     FORM THE PHASES AND PLOT THEM
C
      DO1700I=1,IMAX2
      Z(I)=ATAN(Y1(I)/Y2(I))
 1700 CONTINUE
      CALL PLOTX(Z,IMAX2)
      WRITE(IPR,2020)F
C
C     COMPUTE A NEW ESTIMATE FOR THE DEMODULATION FREQUENCY AND WRITE IT OUT
C
      AIMAX2=IMAX2
      IMAX2M=IMAX2-1
      IFLAG=0
      ZMIN=Z(1)
      ZMAX=Z(1)
      DO1800I=1,IMAX2M
      IP1=I+1
      DEL=Z(IP1)-Z(I)
      IF(DEL.GT.2.5)IFLAG=IFLAG-1
      IF(DEL.LT.-2.5)IFLAG=IFLAG+1
      AIFLAG=IFLAG
      ZNEW=Z(IP1)+AIFLAG*PI
      IF(ZNEW.LT.ZMIN)ZMIN=ZNEW
      IF(ZNEW.GT.ZMAX)ZMAX=ZNEW
 1800 CONTINUE
      RANGE=ZMAX-ZMIN
      SLOPER=RANGE/AIMAX2
      SLOPEH=SLOPER/(2.0*PI)
      FEST=F+SLOPEH 
      WRITE(IPR,2025)ZMIN,ZMAX,RANGE
      WRITE(IPR,2030)SLOPER,SLOPEH,FEST 
C
 2005 FORMAT(1H ,30X, 48HAMPLITUDE PLOT FOR THE DEMODULATION FREQUENCY =
     1 ,F8.6,21H CYCLES PER UNIT TIME)
 2010 FORMAT(1H ,9X,20HMINIMUM AMPLITUDE = ,E15.8,5X,20HMAXIMUM AMPLITUD
     1E = ,E15.8,5X,22HRANGE OF AMPLITUDES = ,E15.8)
 2020 FORMAT(1H ,32X, 44HPHASE PLOT FOR THE DEMODULATION FREQUENCY = ,F8
     1.6,21H CYCLES PER UNIT TIME)
 2025 FORMAT(1H ,3X,16HMINIMUM PHASE = ,E15.8,11H RADIANS   ,16HMAXIMUM
     1PHASE = ,E15.8,11H RADIANS   ,18HRANGE OF PHASES = ,E15.8,8H RADIA
     1NS) 
 2030 FORMAT(1H ,8HSLOPE = ,E14.8,11H RADIANS = ,E14.6,52H CYCLES PER UN
     1IT TIME    EST. OF NEW DEMOD. FREQ. = ,E15.8,15H CYC./UNIT TIME)
C
      RETURN
      END