*======== SUBROUTINE VDTLS ( * * inputs * : X, Y, DX, DY, *==== : NPTS, SCHEME, FRAC, : NCOEFF, DIM, NVAR, * * input/output * : COEFF, TOLERN, ITER, * * outputs * : CHISQ, MATRIX, STATUS) * * Module number: 15.5.2.3 * * Module name: deadtimev * * Keyphrase: * ---------- *======== * Perform a least square fitting to dead time correction equation *==== * * Description: * ------------ * Least-square fit the function: * *======== * * F = A * X / (1 + A * X * tau) - Y * * COEFF(1) = tau, COEFF(2) = A, VAR(1) = X, VAR(2) = Y * * Maximum number of coefficients is 20. * Maximum number of variables is 10. * * FORTRAN name: VDTLS.FOR *==== * * Keywords of accessed files and tables: * -------------------------------------- * Name I/O Description / Comments * * Subroutines Called: * ------------------- * CDBS: *======== * VDTFN, CDINVS, CDWT *==== * SDAS: * None * Others: * None * * History: * -------- * Version Date Author Description * 1 02-16-87 J.-C. HSU design and coding * 2 05-16-87 J.-C. HSU F77 SDAS * *------------------------------------------------------------------------------- * *== input: *======== * --input x, and y array and their errors DOUBLE PRECISION X(1), Y(1), DX(1), DY(1) *==== * --number of input data points INTEGER NPTS * --data weighting scheme REAL SCHEME(1), * --specified fraction applied to the * --coefficients adjustment : FRAC * --number of fitting coefficients INTEGER NCOEFF, * --declared dimensions of MATRIX : DIM, * --number of variables : NVAR * *== input/output: * --coefficients of the fitting function DOUBLE PRECISION COEFF(1) * --tolerance of the degree of sigma-square * --convergence REAL TOLERN * --maximum number of iterations (input) * --number of iterations (output) INTEGER ITER * *== output: * --chi-squared DOUBLE PRECISION CHISQ, * --locally, normal equation matrix and * --its inverse, on output, the covariance * --matrix : MATRIX(DIM, DIM) * --error status INTEGER STATUS * *== local: * --partial derivatives of the fitting * --function relative to coefficients DOUBLE PRECISION FCOEFF(20), * --wieghted sum of functional value times * --FCOEFF : FCF0(20), * --coefficient adjustments : ACOEFF(20), LASTAC(20), * --variables, their errors, their * --weights, and their residuals : VAR(10), DVAR(10), WVAR(10), VVAR(10), * --partial derivatives of the fitting * --function relative to variables : FVAR(10), * --functional value of the fitting * --function and its weighted square : F0, F0F0, * --total weight : WEIGHT, * --sigma-squared of last iteration : LASTS, * --total sigma-squared : S, * --difference between current and * --previous sigma-square : DIFF, * --first order Taylor expansion of F0 : F0PRM * --loop indices INTEGER I, J, K, M * *=========================begin hsp.inc========================================= * --status return code INTEGER OK, ERRNUM(20) * --message destination and priority INTEGER DEST, PRIO DATA OK /0/ DATA ERRNUM /701, 702, 703, 704, 705, 706, 707, 708, 709, 710, : 711, 712, 713, 714, 715, 716, 717, 718, 719, 720/ DATA DEST, PRIO /1, 0/ *=========================end hsp.inc=========================================== *------------------------------------------------------------------------------ * * initialize * DO 10 J = 1, NCOEFF LASTAC(J) = 0.D0 10 CONTINUE * * iterate the least square procedure * DO 140 M = 1, ITER * * reset the normal equations * CHISQ = 0.D0 F0F0 = 0.D0 * DO 30 J = 1, NCOEFF FCF0(J) = 0.D0 ACOEFF(J) = 0.D0 * DO 20 K = 1, NCOEFF MATRIX(J, K) = 0.D0 20 CONTINUE 30 CONTINUE * * summation over each data point * DO 80 I = 1, NPTS * * get functional value and partial derivatives * *======== VAR(1) = X(I) VAR(2) = Y(I) DVAR(1) = DX(I) DVAR(2) = DY(I) * CALL VDTFN (VAR, COEFF, : F0, FVAR, FCOEFF) *==== * * determine the weights at each point * * CALL CDWT (SCHEME, FVAR, DVAR, NVAR, : WVAR, WEIGHT, STATUS) IF (STATUS .NE. OK) THEN GO TO 999 END IF * * calculate residuals (using coefficient adjustment of LAST iteration) and * chi-squared * DO 50 K = 1, NVAR IF (SCHEME(K) .GT. -1000.) THEN F0PRM = F0 DO 40 J = 1, NCOEFF F0PRM = F0PRM - FCOEFF(J) * LASTAC(J) 40 CONTINUE * VVAR(K) = WEIGHT * F0PRM * FVAR(K) / WVAR(K) CHISQ = (VVAR(K) / DVAR(K)) ** 2 + CHISQ END IF 50 CONTINUE * F0F0 = F0 * F0 * WEIGHT + F0F0 * * construct (upper right part of) the normal equations * DO 70 J = 1, NCOEFF FCF0(J) = FCOEFF(J) * F0 * WEIGHT + FCF0(J) * DO 60 K = J, NCOEFF MATRIX(J, K) = FCOEFF(J) * FCOEFF(K) * WEIGHT + : MATRIX(J, K) 60 CONTINUE 70 CONTINUE 80 CONTINUE * * equate the symmetric terms of the matrix (lower left half) * DO 100 J = 2, NCOEFF DO 90 K = 1, J-1 MATRIX(J, K) = MATRIX(K, J) 90 CONTINUE 100 CONTINUE * * invert the matrix * CALL CDINVS (MATRIX, NCOEFF, DIM, : MATRIX, STATUS) IF (STATUS .NE. OK) THEN GO TO 999 END IF * * calculate the adjustment of coefficients and the new coefficients * DO 120 J = 1, NCOEFF DO 110 K = 1, NCOEFF ACOEFF(J) = MATRIX(J, K) * FCF0(K) + ACOEFF(J) 110 CONTINUE * ACOEFF(J) = ACOEFF(J) * FRAC COEFF(J) = COEFF(J) - ACOEFF(J) LASTAC(J) = ACOEFF(J) 120 CONTINUE * * calculate sigma square * S = F0F0 DO 130 J = 1, NCOEFF S = S - FCF0(J) * ACOEFF(J) 130 CONTINUE * * decide whether to terminate the iteration * IF (M .EQ. 1) THEN LASTS = F0F0 END IF * DIFF = ABS(S - LASTS) / LASTS * IF (DIFF .LT. TOLERN .OR. S .LT. 1.D-30) THEN GO TO 150 END IF * LASTS = S 140 CONTINUE * 150 TOLERN = DIFF ITER = MIN (M, ITER) * * calculate the covariance matrix * S = S / DBLE(NPTS - NCOEFF) DO 170 J = 1, NCOEFF DO 160 K = 1, NCOEFF MATRIX(J, K) = S * MATRIX(J, K) 160 CONTINUE 170 CONTINUE * * calculate chi-squared per degree of freedom * CHISQ = CHISQ / DBLE(NPTS - NCOEFF) * STATUS = OK * 999 CONTINUE * RETURN END