mnpfit.f

Go to the documentation of this file.

*
* $Id: mnpfit.F,v 1.1.1.1 1996/03/07 14:31:31 mclareni Exp $
*
* $Log: mnpfit.F,v $
* Revision 1.1.1.1  1996/03/07 14:31:31  mclareni
* Minuit
*
*
      SUBROUTINE mnpfit(PARX2P,PARY2P,NPAR2P,COEF2P,SDEV2P)
*
* $Id: d506dp.inc,v 1.1.1.1 1996/03/07 14:31:32 mclareni Exp $
*
* $Log: d506dp.inc,v $
* Revision 1.1.1.1  1996/03/07 14:31:32  mclareni
* Minuit
*
*
*
*
* d506dp.inc
*
C ************ DOUBLE PRECISION VERSION *************
      IMPLICIT DOUBLE PRECISION (a-h,o-z)
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)
      dimension cz(3)
C
      DO 3  i=1,3
    3 cz(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
      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.)
    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