C***********************************************************************
C    THIS SUBROUTINE CALCULATES AN ESTIMATE OF F(0) AND ITS VARIANCE
C    USING A FOURIER SERIES.  IT NEEDS AS INPUT THE UNGROUPED DATA
C    VECTOR - X,THE MAXIMUM VALUE OF THE DATA - WIDTH, THE NUMBER
C    OF DATA POINTS - CNT, AND WHETHER THE NUMBER OF PARAMETERS HAS
C    BEEN SPECIFIED (SET) AND HOW MANY IF IT IS (NPSET) - ALL ARE
C    IN COMMON.
C
C     SUBROUTINES CALLED: HEADER,NAME
C***********************************************************************
      SUBROUTINE FSEST (IND)
      INCLUDE 'PARMTR.INC'
C***********************************************************************
C     DECLARATIONS
C***********************************************************************
      INTEGER CNT
      LOGICAL YES, WARN, MSET, STRT, SET
C***********************************************************************
C     COMMON STATEMENTS
C***********************************************************************
      COMMON /PDOPT/ STRT, SET, PSTRT(10,MAXPAR), NPSET(10)
      COMMON /SOLN/ FZERO, VARF, FMAX, FMIN, FZ(100)
      COMMON /NUM/ XL(MAXLIN), WIDTH, N, CNT, CONV(3), VARN, IDF, WARN
      COMMON /MEASUR/ X(MAXOBJ), AA(MAXOBJ,4)
      DOUBLE PRECISION PAR, VCMAT, G, XLL
      COMMON /DPAR/ PAR(MAXPR2), VCMAT(MAXPAR,MAXPAR), G(MAXPAR),
     1    XLL, NPAR, INDEX
C***********************************************************************
C     INITIALIZE VARIABLES AND WRITE OUT HEADER STATEMENTS
C***********************************************************************
      NP=MAXPAR
      MSET=SET
      XN=FLOAT(CNT)
      YES=.FALSE.
      FZERO=1./WIDTH
      NPAR=0
      IF (NPSET(IND).EQ.0) MSET=.FALSE.
      IF (MSET) NPAR=NPSET(IND)
      CALL HEADER (1)
      CALL NAME (INDEX)
      WRITE (6,130)
      I=0
      STOP=SQRT(2./(XN+1.))/WIDTH
      DO 10 K=1,NP*2
   10 PAR(K)=0.D0
C***********************************************************************
C     CALCULATE PARAMETER ESTIMATES
C***********************************************************************
   20 I=I+1
      DO 30 K=1,CNT
   30 PAR(I)=PAR(I)+COS((FLOAT(I)*3.14159*X(K))/WIDTH)
      PAR(I)=2.*PAR(I)/(XN*WIDTH)
C***********************************************************************
C     THE LOGICAL VARIABLE YES IS SET TO TRUE IF THE NUMBER OF PARAMETER
C     IN THE MODEL HAS BEEN CALCULATED.
C     THEN PAR(NPAR+1),...,PAR(2*NPAR) ARE CALCULATED TO GET VARIANCES.
C***********************************************************************
      IF (YES.AND.(I.LT.2*NPAR)) GO TO 20
      IF (YES.AND.(I.GE.2*NPAR)) GO TO 60
      FZERO=FZERO+PAR(I)
         WRITE (6,100) I,PAR(I),FZERO,STOP
C***********************************************************************
C     IF THE NUMBER OF PARAMETERS HAS BEEN SET THEN SEE IF THAT NUMBER
C     HAS BEEN REACHED.
C***********************************************************************
      IF (.NOT.MSET) GO TO 40
      IF (I.EQ.NPAR) YES=.TRUE.
      WRITE (6,150)
      GO TO 20
C***********************************************************************
C     SEE IF THE ABSOLUTE VALUE OF THE PARAMETER ESTIMATE IS GREATER
C     THAN THE STOPPING RULE VALUE TO DETERMINE IF THE PARAMETER SHOULD
C     BE INCLUDED IN THE MODEL.
C***********************************************************************
   40 IF (DABS(PAR(I)).GE.STOP) GO TO 50
      YES=.TRUE.
      NPAR=I-1
      FZERO=FZERO-PAR(I)
      WRITE (6,150)
      WRITE (6,110) NPAR
      GO TO 20
C***********************************************************************
C     FOR THE STOPPING RULE APPROACH NEVER LET THE NUMBER OF PARAMETERS
C     EXCEED 6.
C***********************************************************************
   50 IF (I.LT.6) GO TO 20
      YES=.TRUE.
      NPAR=I
      WRITE (6,150)
      WRITE (6,120)
      GO TO 20
C***********************************************************************
C     NEXT CALCULATE THE COVARIANCE MATRIX OF THE COEFFICIENTS IN THE
C     SERIES AND PAR VARIANCE FOR F(0).
C***********************************************************************
   60 VARF=0.0
      CALL HEADER (1)
      CALL NAME (INDEX)
      IF ((NPAR.EQ.0).OR.(XN.EQ.1.)) RETURN
      DO 90 I=1,NPAR
      IF (I.EQ.NPAR) GO TO 80
      DO 70 J=I+1,NPAR
      VCMAT(I,J)=((PAR(I+J)+PAR(J-I))/DBLE(WIDTH)
     1   -PAR(I)*PAR(J))/DBLE(XN-1.)
      VCMAT(J,I)=VCMAT(I,J)
   70 VARF=VARF+2.*SNGL(VCMAT(I,J))
   80 I2=I*2
      VCMAT(I,I)=((PAR(I2)+2./DBLE(WIDTH))/DBLE(WIDTH)
     1   -PAR(I)*PAR(I))/DBLE(XN-1.)
   90 VARF=VARF+SNGL(VCMAT(I,I))
C***********************************************************************
C     FORMAT STATEMENTS
C***********************************************************************
      RETURN
C
  100 FORMAT (4X,'A(',I2.2,')',4X,G12.5,6X,G12.5,7X,G15.8)
  110 FORMAT ('0The absolute value of the parameter is less',
     1' than the stopping rule value.'/
     2' Thus ',I2,' term(s) will be used in the model.')
  120 FORMAT ('0The stopping rule may have allowed more than 6',
     1' terms in the model.'/
     2' This porgram will only use at the most 6 terms unless',
     3' the user specifies a'/
     4' set number of terms greater than 6.')
  130 FORMAT (//'0The ungrouped version of the Fourier series',
     1' estimator determines the'/
     2' number of terms by a stopping rule which minimizes',
     3' the mean integrated'/
     4' squared error.  Below the value for each parameter,',
     5' the cumulative value'/
     6' of F(0) and the stopping rule value are given.'/
     7'0Parameter     Estimate     Cumulative F(0)',
     5 4X,'Stopping Rule Value'/1X,66(1H-))
  150 FORMAT (1X,66(1H-))
      END
