C
C     ..................................................................
C
C        SUBROUTINE DAPCH
C
C        PURPOSE
C           SET UP NORMAL EQUATIONS OF LEAST SQUARES FIT IN TERMS OF
C           CHEBYSHEV POLYNOMIALS FOR A GIVEN DISCRETE FUNCTION
C
C        USAGE
C           CALL DAPCH(DATI,N,IP,XD,X0,WORK,IER)
C
C        DESCRIPTION OF PARAMETERS
C           DATI  - VECTOR OF DIMENSION 3*N (OR DIMENSION 2*N+1)
C                   CONTAINING THE GIVEN ARGUMENTS, FOLLOWED BY THE
C                   FUNCTION VALUES AND N (RESPECTIVELY 1) WEIGHT
C                   VALUES. THE CONTENT OF VECTOR DATI REMAINS
C                   UNCHANGED.
C                   DATI MUST BE OF DOUBLE PRECISION
C           N     - NUMBER OF GIVEN POINTS
C           IP    - DIMENSION OF LEAST SQUARES FIT, I.E. NUMBER OF
C                   CHEBYSHEV POLYNOMIALS USED AS FUNDAMENTAL FUNCTIONS
C                   IP SHOULD NOT EXCEED N
C           XD    - RESULTANT MULTIPLICATIVE CONSTANT FOR LINEAR
C                   TRANSFORMATION OF ARGUMENT RANGE
C                   XD MUST BE DOUBLE PRECISION
C           X0    - RESULTANT ADDITIVE CONSTANT FOR LINEAR
C                   TRANSFORMATION OF ARGUMENT RANGE
C                   X0 MUST BE DOUBLE PRECISION
C           WORK  - WORKING STORAGE OF DIMENSION (IP+1)*(IP+2)/2
C                   ON RETURN WORK CONTAINS THE SYMMETRIC COEFFICIENT
C                   MATRIX OF THE NORMAL EQUATIONS IN COMPRESSED FORM
C                   FOLLOWED IMMEDIATELY BY RIGHT HAND SIDE
C                   AND SQUARE SUM OF FUNCTION VALUES
C                   WORK MUST BE OF DOUBLE PRECISION
C           IER   - RESULTING ERROR PARAMETER
C                   IER =-1 MEANS FORMAL ERRORS IN DIMENSION
C                   IER = 0 MEANS NO ERRORS
C                   IER = 1 MEANS COINCIDING ARGUMENTS
C
C        REMARKS
C           NO WEIGHTS ARE USED IF THE VALUE OF DATI(2*N+1) IS
C           NOT POSITIVE.
C           EXECUTION OF SUBROUTINE DAPCH IS A PREPARATORY STEP FOR
C           CALCULATION OF LEAST SQUARES FITS IN CHEBYSHEV POLYNOMIALS
C           IT SHOULD BE FOLLOWED BY EXECUTION OF SUBROUTINE DAPFS
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           THE LEAST SQUARE FIT IS DETERMINED USING CHEBYSHEV
C           POLYNOMIALS AS FUNDAMENTAL FUNCTION SYSTEM.
C           THE METHOD IS DISCUSSED IN THE ARTICLE
C           A.T.BERZTISS, LEAST SQUARES FITTING TO IRREGULARLY SPACED
C           DATA, SIAM REVIEW, VOL.6, ISS.3, 1964, PP. 203-227.
C
C     ..................................................................
C
      SUBROUTINE DAPCH(DATI,N,IP,XD,X0,WORK,IER)
C
C
C       DIMENSIONED DUMMY VARIABLES
      DIMENSION DATI(1),WORK(1)
      DOUBLE PRECISION DATI,WORK,XD,X0,XA,XE,XM,DF,T,SUM
C
C        CHECK FOR FORMAL ERRORS IN SPECIFIED DIMENSIONS
      IF(N-1)19,20,1
    1 IF(IP)19,19,2
C
C        SEARCH SMALLEST AND LARGEST ARGUMENT
    2 IF(IP-N)3,3,19
    3 XA=DATI(1)
      X0=XA
      XE=0.D0
      DO 7 I=1,N
      XM=DATI(I)
      IF(XA-XM)5,5,4
    4 XA=XM
    5 IF(X0-XM)6,7,7
    6 X0=XM
    7 CONTINUE
C
C        INITIALIZE CALCULATION OF NORMAL EQUATIONS
      XD=X0-XA
      M=(IP*(IP+1))/2
      IEND=M+IP+1
      MT2=IP+IP
      MT2M=MT2-1
C
C        SET WORKING STORAGE AND RIGHT HAND SIDE TO ZERO
      DO 8 I=1,IP
      J=MT2-I
      WORK(J)=0.D0
      WORK(I)=0.D0
      K=M+I
    8 WORK(K)=0.D0
C
C        CHECK FOR DEGENERATE ARGUMENT RANGE
      IF(XD)20,20,9
C
C        CALCULATE CONSTANTS FOR REDUCTION OF ARGUMENTS
    9 X0=-(X0+XA)/XD
      XD=2.D0/XD
      SUM=0.D0
C
C        START GREAT LOOP OVER ALL GIVEN POINTS
      DO 15 I=1,N
      T=DATI(I)*XD+X0
      J=I+N
      DF=DATI(J)
C
C        CALCULATE AND STORE VALUES OF CHEBYSHEV POLYNOMIALS
C        FOR ARGUMENT T
      XA=1.D0
      XM=T
      IF(DATI(2*N+1))11,11,10
   10 J=J+N
      XA=DATI(J)
      XM=T*XA
   11 T=T+T
      SUM=SUM+DF*DF*XA
      DF=DF+DF
      J=1
   12 K=M+J
      WORK(K)=WORK(K)+DF*XA
   13 WORK(J)=WORK(J)+XA
      IF(J-MT2M)14,15,15
   14 J=J+1
      XE=T*XM-XA
      XA=XM
      XM=XE
      IF(J-IP)12,12,13
   15 CONTINUE
      WORK(IEND)=SUM+SUM
C
C        CALCULATE MATRIX OF NORMAL EQUATIONS
      LL=M
      KK=MT2M
      JJ=1
      K=KK
      DO 18 J=1,M
      WORK(LL)=WORK(K)+WORK(JJ)
      LL=LL-1
      IF(K-JJ)16,16,17
   16 KK=KK-2
      K=KK
      JJ=1
      GOTO 18
   17 JJ=JJ+1
      K=K-1
   18 CONTINUE
      IER=0
      RETURN
C
C        ERROR RETURN IN CASE OF FORMAL ERRORS
   19 IER=-1
      RETURN
C
C        ERROR RETURN IN CASE OF COINCIDING ARGUMENTS
   20 IER=1
      RETURN
      END