SUBROUTINE VPNFT2 * * Module number: * * Module name: twodpolyfit * * Keyphrase: * ---------- * Perform 2D polynomial fitting by using singular value decomposition * * Description: * ------------ * * FORTRAN name: VPNFT2.FOR * * Keywords of accessed files and tables: * -------------------------------------- * Name I/O Description / Comments * * Subroutines Called: * ------------------- * CDBS: * VPNGT2, CDSVFT, CDSVDV, VRSCAL, CDGAMD, VPNPT2 * SDAS: * UDMGET, UDMFRE, UMSPUT * Others: * None * * History: * -------- * Version Date Author Description * 1 12-15-88 J.-C. HSU design and coding *------------------------------------------------------------------------------- * *== local: * --pointers of the dependent variable * --array Y and its standard deviation INTEGER Y, YERR, * --pointers of the independent variable * --array X : X, * --pointers of the covariance matrix : COVAR, * --pointers of the decomposition * --component matrices : U, V, W, * --pointers of the working space : WORKB, * --dimension of X : NX, * --total number of valid input data points : NPTS, * --maximum number of X's dimension : NXMAX, * --maximum number of coefficients : MAMAX, * --detector ID : DETID, * --number of attributes : NATTRB, * --error status : STATUS PARAMETER (NXMAX = 7) PARAMETER (MAMAX = 50) * --order of fitting polynomial INTEGER IORDER(NXMAX), IPOWER(NXMAX), : MORDER(NXMAX), * --number of coefficients : MA, * --loop indices, number of input columns : I, J, K, L, M, * --status : STAT(30), STATOK * --high voltage, gain setting, and * --discriminator threshold setting REAL VOLT, GAIN, THRESH, * --reduced chi-squared : RCHISQ(0:9, 0:9), * --minimum probability : PROB0, MAXPRB * --X offsets and extremes DOUBLE PRECISION X0(NXMAX), XMIN(NXMAX), XMAX(NXMAX), * --X scale factors : XS(NXMAX), * --coefficients : COEFF(MAMAX), * --work area : WORK(MAMAX), * --probability : PROB(0:9,0:9), * --chi-squared : CHISQ CHARACTER*3 CHAR3 * --data type CHARACTER*7 TYPE * --aperture name CHARACTER*10 APERT * --attribute names CHARACTER*16 ATTRB(5), * --scheme of executing this task : SCHEME * --output file name CHARACTER*128 OFILE * --error message context CHARACTER*130 CONTXT, MESS *==========================begin iraf77.inc (without INTEGER*2)================= * Include file for the iraf77 FORTRAN interface to the IRAF VOS * Get IRAF common into main program * LOGICAL MEMB(1) INTEGER MEMI(1) INTEGER MEML(1) REAL MEMR(1) DOUBLE PRECISION MEMD(1) COMPLEX MEMX(1) EQUIVALENCE (MEMB, MEMI, MEML, MEMR, MEMD, MEMX) COMMON /MEM/ MEMD * * File I/O access modes * INTEGER RDONLY PARAMETER (RDONLY = 1) INTEGER RDWRIT PARAMETER (RDWRIT = 2) INTEGER WRONLY PARAMETER (WRONLY = 3) INTEGER APPEND PARAMETER (APPEND = 4) INTEGER NEWFIL PARAMETER (NEWFIL = 5) INTEGER TMPFIL PARAMETER (TMPFIL = 6) INTEGER NEWCPY PARAMETER (NEWCPY = 7) INTEGER NEWIMG PARAMETER (NEWIMG = 5) * * codes for data types * INTEGER TYBOOL PARAMETER (TYBOOL = 1) INTEGER TYCHAR PARAMETER (TYCHAR = 2) INTEGER TYSHOR PARAMETER (TYSHOR = 3) INTEGER TYINT PARAMETER (TYINT = 4) INTEGER TYLONG PARAMETER (TYLONG = 5) INTEGER TYREAL PARAMETER (TYREAL = 6) INTEGER TYDOUB PARAMETER (TYDOUB = 7) INTEGER TYCPLX PARAMETER (TYCPLX = 8) INTEGER TYUSHT PARAMETER (TYUSHT = 11) INTEGER TYUBYT PARAMETER (TYUBYT = 12) * * TYTEXT is a special code for the iraf77 interface; it is not in the VOS * INTEGER TYTEXT PARAMETER (TYTEXT = 13) *========================end iraf77.inc========================================= *=========================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=========================================== *------------------------------------------------------------------------------ * * input parameters and data * CALL VPNGT2 ( X, Y, YERR, NX, X0, XMIN, XMAX, IORDER, TYPE, : NPTS, SCHEME, PROB0, NATTRB, ATTRB, : DETID, APERT, VOLT, GAIN, THRESH, OFILE, STATUS) * IF (STATUS .NE. OK) THEN CONTXT = 'cannot input parameters or data' GO TO 999 END IF * * determine the factor used to scale X vector * MA = 1 DO 10 I = 1, NX MA = MA * (IORDER(I) + 1) XS(I) = (XMAX(I) - XMIN(I)) / 2.D0 IF (XS(I) .EQ. 0.D0) XS(I) = 1.D0 10 CONTINUE * IF (MA .GT. MAMAX) THEN STATUS = ERRNUM(1) WRITE (CHAR3, '(I3)') MAMAX CONTXT = 'number of coefficients is greater than ' // CHAR3 GO TO 999 END IF * * allocate dynamic memory for (decomposition) arrays * CALL UDMGET (NPTS*MA, TYDOUB, U, STAT(1)) CALL UDMGET (MA**2, TYDOUB, V, STAT(2)) CALL UDMGET (MA, TYDOUB, W, STAT(3)) CALL UDMGET (NPTS, TYDOUB, WORKB, STAT(4)) CALL UDMGET (MA**2, TYDOUB, COVAR, STAT(5)) * DO 20 I = 1, 5 IF (STAT(I) .NE. OK) THEN STATUS = ERRNUM(2) CONTXT = 'cannot allocate dynamic memory' GO TO 999 END IF 20 CONTINUE * * reset dynamic memory to zero * DO 50 I = 1, MA DO 30 J = 1, MA L = MA * (I-1) + J MEMD(V+L-1) = 0.D0 MEMD(COVAR+L-1) = 0.D0 30 CONTINUE * MEMD(W+I-1) = 0.D0 * DO 40 K = 1, NPTS M = NPTS * (I-1) + K MEMD(U+M-1) = 0.D0 40 CONTINUE 50 CONTINUE * DO 60 I = 1, NPTS MEMD(WORKB+I-1) = 0.D0 60 CONTINUE * * if the scheme is 'auto' or 'exam', do all polynomial fits with orders no * larger than the specified orders * IF (SCHEME(1:4) .EQ. 'auto' .OR. SCHEME(1:4) .EQ. 'exam') THEN * DO 80 I = 0, IORDER(1) IPOWER(1) = I DO 70 J = 0, IORDER(2) IPOWER(2) = J * MA = (IPOWER(1)+1) * (IPOWER(2)+1) * * perform the SVD polynomial fit * CALL CDSVFT (MEMD(X), MEMD(Y), MEMD(YERR), X0, XS, : IPOWER, WORK, NPTS, NX, NX, MA, NPTS, : MA, MEMD(WORKB), COEFF, MEMD(U), : MEMD(V), MEMD(W), CHISQ, STATUS) IF (STATUS .NE. OK) THEN CONTXT = 'SVD fitting fails' GO TO 999 END IF * * reduced chi-squared * RCHISQ(I, J) = REAL (CHISQ / (NPTS - MA)) * * calculate chi-squared probability function * IF (NPTS .GT. MA) THEN CALL CDGAMD (0.5D0 * DBLE (NPTS - MA), : 0.5D0 * CHISQ, PROB(I,J), STATUS) PROB(I,J) = 1.D0 - PROB(I,J) IF (STATUS .NE. OK) GO TO 999 ELSE PROB(I,J) = 1.D0 END IF * 70 CONTINUE 80 CONTINUE END IF * MAXPRB = 0. MORDER(1) = 0 MORDER(2) = 0 * * find the largest probability and use the corresponding orders * to redo the fit as final result * DO 100 J = 0, IORDER(2) DO 90 I = 0, IORDER(1) IF (PROB(I,J) .GE. PROB0 .OR. PROB(I,J) .GT. MAXPRB) : THEN MAXPRB = PROB(I,J) MORDER(1) = I MORDER(2) = J * * if probability larger than threshold, use the corresponding orders * note that this loop ensures that the smallest time-dependence order is used * IF (PROB(I,J) .GE. PROB0) GO TO 110 END IF 90 CONTINUE 100 CONTINUE 110 CONTINUE * * write probabilities and reduced chi-square to the screen and log file * IF (SCHEME(1:4) .EQ. 'exam') THEN CALL UMSPUT ('chi-square probabilities: x1 down, x2 across', : DEST, PRIO, STATOK) WRITE (CONTXT, '(3X, 10I11)') (J, J = 0, IORDER(2)) CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) * DO 120 I = 0, IORDER(1) WRITE (CONTXT, '(I3, 10E11.3)') I, : (PROB(I, J), J = 0, IORDER(2)) CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) 120 CONTINUE * CALL UMSPUT ('reduced chi-squares: x1 down, x2 across', : DEST, PRIO, STATOK) WRITE (CONTXT, '(3X, 10I11)') (J, J = 0, IORDER(2)) CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) * DO 130 I = 0, IORDER(1) WRITE (CONTXT, '(I3, 10E11.3)') I, : (RCHISQ(I, J), J = 0, IORDER(2)) CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) 130 CONTINUE * ELSE IF (SCHEME(1:4) .EQ. 'auto') THEN IORDER(1) = MORDER(1) IORDER(2) = MORDER(2) END IF * MA = (IORDER(1)+1) * (IORDER(2)+1) * * perform another SVD polynomial fit * CALL CDSVFT (MEMD(X), MEMD(Y), MEMD(YERR), X0, XS, : IORDER, WORK, NPTS, NX, NX, MA, NPTS, : MA, MEMD(WORKB), COEFF, MEMD(U), : MEMD(V), MEMD(W), CHISQ, STATUS) IF (STATUS .NE. OK) THEN CONTXT = 'SVD fitting fails' GO TO 999 END IF * * calculate covariance matrix * CALL CDSVDV (MEMD(V), MA, MA, MEMD(W), WORK, : MA, MEMD(COVAR)) * * unscale the coefficients and covariance matrix * CALL VRSCAL (NX, XS, IORDER, MA, COEFF, MEMD(COVAR)) * * calculate chi-squared probability function * IF (NPTS .GT. MA) THEN CALL CDGAMD (0.5D0 * DBLE (NPTS - MA), : 0.5D0 * CHISQ, WORK(1), STATUS) MAXPRB = 1. - WORK(1) ELSE MAXPRB = 1. END IF * CHISQ = CHISQ / DBLE(NPTS - MA) * * put the result in a table * CALL VPNPT2 (OFILE, IORDER, NPTS, NX, MA, TYPE, COEFF, : MEMD(COVAR), CHISQ, MAXPRB, X0, XMIN, XMAX, : DETID, APERT, VOLT, GAIN, THRESH, NATTRB, : ATTRB, STATUS) * IF (STATUS .NE. OK) THEN CONTXT = 'cannot write result to output table' GO TO 999 END IF END IF * * deallocate dynamic memory * CALL UDMFRE (X, TYDOUB, STAT(1)) CALL UDMFRE (Y, TYDOUB, STAT(2)) CALL UDMFRE (YERR, TYDOUB, STAT(3)) CALL UDMFRE (U, TYDOUB, STAT(4)) CALL UDMFRE (V, TYDOUB, STAT(5)) CALL UDMFRE (W, TYDOUB, STAT(6)) CALL UDMFRE (WORKB, TYDOUB, STAT(7)) CALL UDMFRE (COVAR, TYDOUB, STAT(8)) * DO 190 I = 1, 8 IF (STAT(I) .NE. OK) THEN STATUS = ERRNUM(1) CONTXT = 'cannot free dynamic memory' GO TO 999 END IF 190 CONTINUE * STATUS = OK GO TO 1000 * * write error message * 999 MESS = 'VPNFT2: ' // CONTXT CALL UMSPUT (MESS, DEST, PRIO, STATOK) * 1000 RETURN END