SUBROUTINE WORTHO 1 (X,Y,W,NPTS,COEFF,N,ERROR,WORK) C C This subroutine calculates the coefficients of the Nth order C polynomial which provides the best least-squares fit to the C weighted set of points, (X,Y):W. C C X Array of independent variables. C Y Array of dependent variables. C W Array of weights. C NPTS In absolute value, the number of data points. If negative C the error will be expressed as a percentage. C COEFF Coefficients of the fitted polynomial, stored such that C Y(X) = COEFF(1) + ... + COEFF(N+1)*X**N. COEFF must C be dimensioned at least N+1. C N Order of the polynomial, N < 6. C ERROR Array containing the differences between the actual C values and those calculated by the polynomial. C Percentage error are stored if NPTS < 0. Also C stored are the maximum error in ERROR( IABS(NPTS) + 1 ) C and the sigma error in ERROR( IABS(NPTS) + 2 ). ERROR C must be dimensioned at least IABS(NPTS)+2. C WORK Working array, dimensioned at least 2*IABS(NPTS). C REAL A(6,6),B(6),C(6),D(6),S(6) REAL X(1),Y(1),W(1),COEFF(1),ERROR(1),WORK(1) MPTS = IABS(NPTS) NORDER = N+1 C C - ZERO WORKING ARRAYS. DO 1000 J=1,NORDER B(J) = 0.0 COEFF(J) = 0.0 D(J) = 0.0 S(J) = 0.0 DO 1000 I=1,NORDER 1000 A(I,J) = 0.0 C C - CALCULATE SUM OF X AND Y VALUES FOR AVERAGE. DO 2000 I=1,MPTS B(1) = B(1)+X(I)*W(I) D(1) = D(1)+Y(I)*W(I) 2000 S(1) = S(1) + W(I) C(1) = 0.0 D(1) = D(1)/S(1) C C - CALCULATE DIFFERENCE (ERROR) BETWEEN Y VALUES AND Y AVERAGE. DO 3000 I=1,MPTS 3000 ERROR(I) = Y(I)-D(1) IF( NORDER .EQ. 1 ) GO TO 7100 B(1) = B(1)/S(1) C C - CALCALUTE FIRST ORDER MOMENT ARMS. DO 4000 I=1,MPTS WORK( MPTS+I ) = 1. 4000 WORK(I) = X(I)-B(1) J=1 4100 J = J+1 C C - CALCULATE JTH ORDER MOMENTS. DO 5000 I=1,MPTS D(J) = D(J)+ERROR(I)*WORK(I)*W(I) B(J) = B(J)+X(I)*WORK(I)**2 5000 S(J) = S(J)+WORK(I)**2 D(J) = D(J)/S(J) C C - CALCULATE ERROR FROM JTH ORDER CURVE. DO 6000 I=1,MPTS 6000 ERROR(I) = ERROR(I)-D(J)*WORK(I) IF( J .EQ. NORDER ) GO TO 7100 B(J) = B(J)/S(J) C(J) = S(J)/S(J-1) C C - CALCULATE NEXT ORDER MOMENT ARMS. DO 7000 I=1,MPTS WW = WORK(I) WORK(I) = (X(I)-B(J))*WORK(I)-C(J)*WORK( MPTS+I ) 7000 WORK( MPTS+I ) = WW GO TO 4100 C C - SET UP COEFFICIENT ARRAY. 7100 A(1,1) = 1.0 IF( NORDER .EQ. 1 ) GO TO 8100 A(2,2) = A(1,1) A(2,1) = -B(1)*A(1,1) IF( NORDER .EQ. 2 ) GO TO 8100 DO 8000 J=3,NORDER A(J,1) = -B(J-1)*A(J-1,1)-C(J-1)*A(J-2,1) DO 8000 M=2,NORDER 8000 A(J,M) = A(J-1,M-1)-B(J-1)*A(J-1,M)-C(J-1)*A(J-2,M) C C - REDUCE COEFFICIENT MATRIX TO OBTAIN COEFFICIENTS. 8100 DO 9000 I=1,NORDER DO 9000 J=I,NORDER 9000 COEFF(I) = COEFF(I)+D(J)*A(J,I) C C - CALCULATE ERROR. ERROR(MPTS+1) = -1. ERROR(MPTS+2) = 0.0 DO 10000 I=1,MPTS ERROR(MPTS+2) = ERROR(MPTS+2)+ERROR(I)**2 IF( NPTS .GT. 0 ) GO TO 9200 C C - CALCULATE PERCENTAGE ERROR. IF( ABS(Y(I)) .GT. 1.E-6 ) GO TO 9100 ERROR(I) = 1.E10 GO TO 10000 9100 ERROR(I) = 1.0-ERROR(I)/Y(I) 9200 IF( ABS(ERROR(I)) .LE. ERROR(MPTS+1) ) GO TO 10000 C C - FIND MAXIMUM ERROR. ERROR(MPTS+1) = ABS(ERROR(I)) 10000 CONTINUE C C - CACULATE SIGMA. ERROR(MPTS+2) = SQRT(ERROR(MPTS+2)/MPTS) C RETURN END