SUBROUTINE NLFIT (NDATA, X, Y, DY, NMODEL, A, DA, CHISQ, FLAMDA) C C************************************************************************ C * C Performs a fit of the non-linear model F[A(1),A(2),...,A(N);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 F is a user * C defined function. * C * C NOTE: This routine does not check for the goodness of the fit. * C The calling routine should call NLFIT repetitively, * C checking the value of CHISQ on successive calls until * C it does not change significantly. The variable FLAMDA * C determines the amount of gradient search included. * C The starting value of FLAMDA should be 0.001; NLFIT * C will adjust FLAMDA as required between each call. * C It is also necessary to supply an initial guess at * C the values of the model parameters; the better this * C guess, the quicker will be the fit. * 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 NMODEL - The number of model parameters to be fit. * C A - An initial guess at the values of the parameters. * C FLAMDA - Proportion of gradient search included. * C A recommended starting value of FLAMDA is 0.001 * C * C OUTPUTS * C A - The fitted parameters. (dimension >= NMODEL) * C DA - Parameter uncertainties. (dimension >= NMODEL) * C CHISQ - Chi-squared per degree of freedom for the fit. * C * C OTHER ROUTINES REQUIRED * C FUNCTN (X, A, NMODEL) * C Evaluates the model function at a given X for * C a given set of model parameters. * C FCHISQ (X, Y, DY, A, NMODEL) * C Evaluates reduced chi-squared for a given set * C of model parameters. * C FDERIV (DERIV, X, A, NMODEL) * C Evaluates the derivatives of the fitting function * C with respect to the model parameters evaluated at * C a particular value of X and for a particular set * C of model parameters. * C MATINV (ARRAY, NORDER, DET) * C Inverts a symmetric double precision matrix of * C specified order. * C * C RESTRICTIONS * C NMODEL should not exceed 10. * 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(100), ARRAY(100), G(10), DERIV(10), 1 B(10), C(10), WI, FUNCTN EXTERNAL FUNCTN C C.......10 model parameters is the maximum allowed C CHISQ = 0. NORDER = NMODEL IF (NORDER .GT. 10) RETURN C C.......Be sure there is enough data to fit C NFREE = NDATA - NMODEL IF (NDATA .LT. NMODEL) RETURN C C.......Normalize the uncertainties C DYNORM = DY(1) DO 10 I = 1,NDATA DY(I) = DY(I) / DYNORM 10 CONTINUE C C.......Initialize the G vector and the H matrix C 20 CONTINUE DO 40 K = 1,NMODEL B(K) = A(K) G(K) = 0D0 KK = (K-1)*NMODEL DO 30 J=1,K H(J+KK) = 0D0 30 CONTINUE 40 CONTINUE C C.......Compute G vector and half of H matrix C DO 70 I=1, NDATA CALL FDERIV (DERIV, X(I), B, NMODEL) YMODEL = FUNCTN (X(I), B, NMODEL) WI = 1D0 / DBLE(DY(I))**2 DO 60 K = 1,NMODEL G(K) = G(K) + DERIV(K) * DBLE(Y(I) - YMODEL) * WI KK = (K-1)*NMODEL DO 50 J = 1,K H(J+KK) = H(J+KK) + DERIV(J) * DERIV(K) * WI 50 CONTINUE 60 CONTINUE 70 CONTINUE C C.......Fill in rest of H matrix C DO 90 J = 2,NMODEL JJ = (J-1)*NMODEL DO 80 K = 1,J-1 KK = (K-1)*NMODEL H(J+KK) = H(K+JJ) 80 CONTINUE 90 CONTINUE C C.......Evaluate initial chi-squared C CHISQI = FCHISQ (NDATA, X, Y, DY, NMODEL, B) C C.......Normalize H and add FLAMDA to diagonal C 100 CONTINUE DO 110 J=1, NMODEL JJ = (J-1)*NMODEL C(J) = DSQRT( DABS( H(J+JJ) ) ) 110 CONTINUE C 120 CONTINUE DO 140 J=1, NMODEL-1 JJ = (J-1)*NMODEL ARRAY(J+JJ) = DBLE(1.+FLAMDA) DO 130 K=J+1, NMODEL KK = (K-1)*NMODEL ARRAY(J+KK) = H(J+KK) / (C(J)*C(K)) ARRAY(K+JJ) = ARRAY(J+KK) 130 CONTINUE 140 CONTINUE ARRAY(NMODEL*NMODEL) = DBLE(1.+FLAMDA) C C.......Invert the curvature matrix C CALL MATINV (ARRAY, NMODEL, DET) C C.......Compute the new parameters in B C DO 160 J=1, NMODEL B(J) = A(J) DO 150 K=1, NMODEL KK = (K-1)*NMODEL B(J) = B(J) + G(K)*ARRAY(J+KK)/(C(J)*C(K)) 150 CONTINUE 160 CONTINUE C C.......If chi-squared went up, multiply FLAMDA by 10 and try again C CHISQ = FCHISQ (NDATA, X, Y, DY, NMODEL, B) IF (CHISQI .GE. CHISQ) GO TO 170 FLAMDA = 10.*FLAMDA GO TO 120 170 CONTINUE C C.......Restore parameters and calculate their uncertainties C DO 180 J=1, NMODEL JJ = (J-1)*NMODEL A(J) = B(J) DA(J) = SQRT( ARRAY(J+JJ)/H(J+JJ) ) * DYNORM 180 CONTINUE C DO 190 I=1,NDATA DY(I) = DY(I)*DYNORM 190 CONTINUE C CHISQ = CHISQ / DYNORM**2 FLAMDA = FLAMDA/10. C RETURN END