SUBROUTINE LIN FIT (NDATA, X, Y, DY, NMODEL, A, DA, CHISQ) C C************************************************************************ C * C Performs a fit of the linear model * C A(1)*F(1,X) + A(2)*F(2,X) + ... + A(NMODEL)*F(NMODEL,X) * C to the given data [X(i),Y(i),DY(i), i=1,NDATA) to determine * C the coefficients A(1),A(2),...A(NMODEL), where the F(J,X) * C are NMODEL user-defined functions. * C * C************************************************************************ C * C INPUTS * C NDATA - The number of data points to be fit * C X - Data x-values (dimension >= NDATA) * C Y - Data y-values (dimension >= NDATA) * C DY - Data uncertainties (dimension >= NDATA) * C NOTE: If the uncertainties are unknown, set * C them all equal to one. * C NMODEL - The number of model parameters to be fit * C * C OUTPUT * C A - The fitted parameters (dimension >= NMODEL) * C DA - Uncertainties in parameters(dimension >= NMODEL) * C CHISQ - Chi-squared per degree of freedom for the fit * C * C OTHER ROUTINES REQUIRED * C MATINV - Inverts symmetric matrix * C F - User-supplied function routine to return values * C of the NMODEL functions F(J,X) * C * C RESTRICTIONS * C NMODEL should not exceed 20. If larger values are * C desired, change the dimension of G to NMODEL and the * C dimension of A to NMODEL*NMODEL. * C * C AUTHOR * C Robert Walraven * C Department of Applied Science * C University of California * C Davis, CA 95616 * C (916) 752-0360 * C * C************************************************************************ C DIMENSION X(1), Y(1), DY(1), A(1), DA(1) DOUBLE PRECISION H(400), G(20), F, DUMMY EXTERNAL F C C.......Calculate the G vector C DO 20 K=1,NMODEL G(K) = 0. DO 10 I=1,NDATA DUMMY = DY(I) G(K) = G(K) + DBLE(Y(I)) * F(K,X(I)) / (DUMMY*DUMMY) 10 CONTINUE 20 CONTINUE C C.......Calculate the H matrix C DO 50 K=1,NMODEL KK = (K-1)*NMODEL DO 40 J=1,NMODEL JK = J+KK H(JK) = 0D0 DO 30 I=1,NDATA DUMMY = DY(I) H(JK) = H(JK) + F(J,X(I))*F(K,X(I)) / (DUMMY*DUMMY) 30 CONTINUE 40 CONTINUE 50 CONTINUE C C.......Invert H in place C CALL MATINV (H, NMODEL, DET) C C.......Compute parameter values C DO 70 K=1,NMODEL DUMMY = 0D0 KK = (K-1)*NMODEL DO 60 J=1,NMODEL DUMMY = DUMMY + G(J) * H(J+KK) 60 CONTINUE A(K) = DUMMY 70 CONTINUE C C.......Compute parameter uncertainties C S2 = 1. DO 80 I=1,NDATA IF (DY(I).NE.1.) GO TO 110 80 CONTINUE S2 = 0. DO 100 I=1,NDATA YM = 0. DO 90 K=1,NMODEL YM = YM + A(K)*F(K,X(I)) 90 CONTINUE S2 = S2 + (Y(I)-YM)**2 100 CONTINUE S2 = S2 / FLOAT(NDATA-NMODEL-1) C 110 CONTINUE DO 120 K=1,NMODEL DA(K) = SQRT(S2) * DSQRT( DABS(H(K+(K-1)*NMODEL)) ) 120 CONTINUE C C.......Compute chi-squared per degree of freedom C CHISQ = 0. DO 140 I=1,NDATA FIT = 0. DO 130 J=1,NMODEL FIT = FIT + A(J) * F(J,X(I)) 130 CONTINUE CHISQ = CHISQ + (Y(I)-FIT)**2/DY(I)**2 140 CONTINUE CHISQ = CHISQ / FLOAT(NDATA-NMODEL) / S2 C RETURN END