C
C
      SUBROUTINE LSFAN (L,QKI,RKI,TABF,NTW,NEXP)
C     ******************************************
C
C     ------------------------------------------------------------------
C
C       SUBROUTINE PURPOSE:
C       USING LEAST SQUARE METHOD FOR FAN CURVE FITTING.
C
C     ------------------------------------------------------------------
C
      INCLUDE 'CMMNW.DAT'
      DIMENSION ZZ(IMX),S(IMX),ALPHA(IMX),BETA(IMX),SIGMA2(IMX),
     .PPP(IMX,IMX)
      DOUBLE PRECISION W(IMX)
C
      NFCW(L)=0
      NP=MPTS(L)
      IF (NTW.EQ.1) GO TO 145
      MARKXX=0
      DO 1 I=1,10
1     ALPHA(I)=0.0
      IF (NP.GT.2) GO TO 35
      ND(L)=1
      N1A=2
      NCOF(L)=N1A
      IF (NP.LT.1) THEN
         STOP 10
      ELSEIF (NP.EQ.1) THEN
         COF(L,1)=PF(L,1)
         COF(L,2)=0.
         GO TO 145
      ENDIF
      IF (QF(L,1).EQ.QF(L,2)) THEN
         STOP 25
      ELSE
         COF(L,1)=(PF(L,1)*QF(L,2)-PF(L,2)*QF(L,1))/(QF(L,2)-QF(L,1))
         COF(L,2)=(PF(L,2)-PF(L,1))/(QF(L,2)-QF(L,1))
         GO TO 145
      ENDIF
35    DSQ=0.0
      DO 40 J=1,NP
         DSQ=DSQ+PF(L,J)*PF(L,J)
         PPP(1,J)=0.
40       PPP(2,J)=1.
      W(1)=NP
      BETA(1)=0.
      N1A=4
      IF (NP.LT.5) N1A=NP-1
      IF (NEXP.GT.0) N1A=NEXP
      DO 85 I=1,N1A
         K=I+1
         ZZ(I)=0.
         DO 60 J=1,NP
60          ZZ(I)=ZZ(I)+PF(L,J)*PPP(K,J)
         S(I)=ZZ(I)/W(I)
         DSQ=DSQ-S(I)*S(I)*W(I)
         SIGMA2(I)=DSQ/FLOAT(NP-I)
         IF(SIGMA2(I).LE.1.D-6) GO TO 105
         IF (I.GE.N1A) GO TO 90
            ALPHA(I)=0.
            DO 75 J=1,NP
75          ALPHA(I)=ALPHA(I)+QF(L,J)*PPP(K,J)*PPP(K,J)
            ALPHA(I)=ALPHA(I)/W(I)
            W(K)=0.
            DO 80 J=1,NP
               PPP(K+1,J)=(QF(L,J)-ALPHA(I))*PPP(K,J)-BETA(I)*PPP(I,J)
80             W(K)=W(K)+PPP(K+1,J)*PPP(K+1,J)
85    BETA(K)=W(K)/W(I)
90    SMALL=SIGMA2(1)
      I=1
      DO 100 J=2,N1A
         IF (SMALL.LE.SIGMA2(J)) GO TO 100
         I=J
         SMALL=SIGMA2(J)
100   CONTINUE
105   WRITE (8,195) SIGMA2(I)
      N1A=I
      NCOF(L)=N1A
      N=N1A-1
      IF (N.LE.0) THEN
         ND(L)=1
         COF(L,1)=S(1)
         COF(L,2)=0.
         GO TO 145
      ENDIF
      ND(L)=N
      DO 125 I=1,N1A
         DO 120 J=1,I
120         PPP(I,J)=0.0
125      PPP(I,I+1)=1.
      DO 130 J=1,N
130      PPP(1,J+2)=-PPP(1,J+1)*ALPHA(J)-PPP(1,J)*BETA(J)
      DO 135 I=2,N
         DO 135 J=I,N
135         PPP(I,J+2)=PPP(I-1,J+1)-PPP(I,J+1)*ALPHA(J)-PPP(I,J)*BETA(J)
      DO 140 I=1,N1A
         COF(L,I)=0.
         DO 140 J=1,N1A
140   COF(L,I)=COF(L,I)+PPP(I,J+1)*S(J)
      NSKP(L)=1
      WRITE (8,*) 'ENTER FAN CURVE FITTING'
145   IF (MARKXX.EQ.0) WRITE (8,200) (COF(L,I),I=1,N1A)
      IF (ABS(QKI).LT.1.D-6) GO TO 160
      MARKXX=1
      IF (NPLOT.EQ.1) THEN
         IF (QKI.LT.QF(L,1).OR.QKI.GT.QF(L,NP)) NFCW(L)=1
         FANQ=QKI
         RGRAD(L)=0.0
         II=1
         GO TO 165
      ELSE IF (NPLOT.EQ.2) THEN
         IF (QKI.LT.QF(L,NP)) GO TO 146
         RKI=PF(L,NP)
         RGRAD(L)=0.0
         NFCW(L)=1
         RETURN
146      IF (QKI.GT.QF(L,1)) GO TO 147
         NFCW(L)=1
147      FANQ=QKI
         RGRAD(L)=0.0
         II=1
         GO TO 165
      ELSE IF (NPLOT.EQ.3) THEN
         IF (QKI.LT.QF(L,NP)) GO TO 150
         RKI=PF(L,NP)
         RGRAD(L)=0.0
         NFCW(L)=1
         RETURN
150      IF (QKI.GT.QF(L,1)) GO TO 155
         RKI=PF(L,1)
         RGRAD(L)=0.0
         NFCW(L)=1
         RETURN
155      FANQ=QKI
         RGRAD(L)=0.0
         II=1
         GO TO 165
      ENDIF
160   QCS=QF(L,NP)-QF(L,1)
      QCSA=QCS/120.
      RGRAD(L)=0.0
      Q0=QF(L,1)-5.*QCSA
      Q2(1)=Q0
      FANQ=Q0
      II=130
165   NE=ND(L)
      DO 180 I=1,II
         H2(I)=COF(L,NE+1)*FANQ
         IF (NE.LE.1) GO TO 175
         KT=NE
         DO 170 J=2,NE
            H2(I)=(H2(I)+COF(L,KT))*FANQ
            IF (MARKXX.EQ.1) RGRAD(L)=(RGRAD(L)+KT*COF(L,KT+1))*FANQ
170         KT=KT-1
175      H2(I)=H2(I)+COF(L,1)
         RGRAD(L)=RGRAD(L)+COF(L,2)
         IF (MARKXX.EQ.1) RKI=H2(1)
         Q2(I+1)=FANQ+QCSA
         FANQ=Q2(I+1)
180   CONTINUE
      IF (MARKXX.EQ.1) RETURN
      WRITE (6,190) (Q2(I),H2(I),I=1,130)
190   FORMAT (1X,5(F8.5,F6.2))
195   FORMAT (/18X,8HSIGMA2= ,1PE14.7)
200   FORMAT (1H0,17X,12HCOEFFICIENTS,//8X,1P6E18.7)
      RETURN
      END



