*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1996/04/01 15:02:28  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
#if defined(CERNLIB_CDC)
      SUBROUTINE PARLSQ
     1    (PARX2P,PARY2P,NPAR2P,COEF2P,SDEV2P)
C
C     TO FIT A PARABOLA TO NPAR2P POINTS
C
C   NPAR2P   NO. OF POINTS
C   PARX2P(I)   X VALUE OF POINT I
C   PARY2P(I)   Y VALUE OF POINT I
C
C   COEF2P(1...3)  COEFFICIENTS OF THE FITTED PARABOLA
C   Y=COEF2P(1) + COEF2P(2)*X + COEF2P(3)*X**2
C   SDEV2P= VARIANCE
C   METHOD : CHI**2 = MIN EQUATION SOLVED EXPLICITLY
      DIMENSION PARX2P(NPAR2P),PARY2P(NPAR2P),COEF2P(NPAR2P)
C
      DO 3  I=1,3
    3 COEF2P(I)=0.
      SDEV2P=0.
      IF(NPAR2P.LT.3) GO TO 10
      F=NPAR2P
C--- CENTER X VALUES FOR REASONS OF MACHINE PRECISION
      XM=0.
      DO 2  I=1,NPAR2P
    2 XM=XM+PARX2P(I)
      XM=XM/F
      X2=0.
      X3=0.
      X4=0.
      Y=0.
      Y2=0.
      XY=0.
      X2Y=0.
      DO 1  I=1,NPAR2P
      S=PARX2P(I)-XM
      T=PARY2P(I)
      S2=S*S
      X2=X2+S2
      X3=X3+S*S2
      X4=X4+S2*S2
      Y=Y+T
      Y2=Y2+T*T
      XY=XY+S*T
      X2Y=X2Y+S2*T
    1 CONTINUE
      A=(F*X4-X2**2)*X2-F*X3**2
      IF(A.EQ.0.)  GOTO 10
      COEF2P(3)=(X2*(F*X2Y-X2*Y)-F*X3*XY)/A
      COEF2P(2)=(XY-X3*COEF2P(3))/X2
      COEF2P(1)=(Y-X2*COEF2P(3))/F
      IF(NPAR2P.EQ.3)  GOTO 6
      SDEV2P=Y2-(COEF2P(1)*Y+COEF2P(2)*XY+COEF2P(3)*X2Y)
      IF(SDEV2P.LT.0.)  SDEV2P=0.
      SDEV2P=SDEV2P/(F-3)
    6 COEF2P(1)=COEF2P(1)+XM*(XM*COEF2P(3)-COEF2P(2))
      COEF2P(2)=COEF2P(2)-2.*XM*COEF2P(3)
   10 CONTINUE
      RETURN
      END
#endif
#if !defined(CERNLIB_CDC)
      SUBROUTINE PARLSQ
     1    (PARX2P,PARY2P,NPAR2P,COEF2P,SDEV2P)
C
C     TO FIT A PARABOLA TO NPAR2P POINTS
C
C   NPAR2P   NO. OF POINTS
C   PARX2P(I)   X VALUE OF POINT I
C   PARY2P(I)   Y VALUE OF POINT I
C
C   COEF2P(1...3)  COEFFICIENTS OF THE FITTED PARABOLA
C   Y=COEF2P(1) + COEF2P(2)*X + COEF2P(3)*X**2
C   SDEV2P= VARIANCE
C   METHOD : CHI**2 = MIN EQUATION SOLVED EXPLICITLY
      DIMENSION PARX2P(NPAR2P),PARY2P(NPAR2P),COEF2P(NPAR2P)
      DOUBLE PRECISION A,CZ(3),F,S,S2,T,XM,XY,X2,X2Y,X3,X4,Y,Y2
C
      DO 3  I=1,3
    3 CZ(I)=0.D0
      SDEV2P=0.
      IF(NPAR2P.LT.3) GO TO 10
      F=NPAR2P
C--- CENTER X VALUES FOR REASONS OF MACHINE PRECISION
      XM=0.D0
      DO 2  I=1,NPAR2P
    2 XM=XM+PARX2P(I)
      XM=XM/F
      X2=0.D0
      X3=0.D0
      X4=0.D0
      Y=0.D0
      Y2=0.D0
      XY=0.D0
      X2Y=0.D0
      DO 1  I=1,NPAR2P
      S=PARX2P(I)-XM
      T=PARY2P(I)
      S2=S*S
      X2=X2+S2
      X3=X3+S*S2
      X4=X4+S2*S2
      Y=Y+T
      Y2=Y2+T*T
      XY=XY+S*T
      X2Y=X2Y+S2*T
    1 CONTINUE
      A=(F*X4-X2**2)*X2-F*X3**2
      IF(A.EQ.0.D0)  GOTO 10
      CZ(3)=(X2*(F*X2Y-X2*Y)-F*X3*XY)/A
      CZ(2)=(XY-X3*CZ(3))/X2
      CZ(1)=(Y-X2*CZ(3))/F
      IF(NPAR2P.EQ.3)  GOTO 6
      SDEV2P=Y2-(CZ(1)*Y+CZ(2)*XY+CZ(3)*X2Y)
      IF(SDEV2P.LT.0.)  SDEV2P=0.
      SDEV2P=SDEV2P/(F-3.D0)
    6 CZ(1)=CZ(1)+XM*(XM*CZ(3)-CZ(2))
      CZ(2)=CZ(2)-2.*XM*CZ(3)
   10 CONTINUE
      DO 11  I=1,3
   11 COEF2P(I)=CZ(I)
      RETURN
      END
#endif