SUBROUTINE VTEMP * * Module Number: 15.1.0.1 * * Module Name: polyfit * * Keyphrase: * ---------- * Analyze a specified physical quantity as a polynomial function of another * specified quantity * * Description: * ------------ * Fit a specified physical quantity (e. g. dark count) as a polynomial of * another specified physical quantity (e. g. detector temperature). * The order of the fitting polynomial must be between 1 and 9. * Maximum number of input data points is 2000. * Input data are in single precision but are double precision during the * polynomial regression calculation. * * FORTRAN Name: VPOLYFIT.FOR * * Keywords of Accessed Files and Tables: * -------------------------------------- * Name I/O Description / Comments * * Subroutines Called: * ------------------- * CDBS: * CDPFIT, CDPBAS, VTMPGT, VTMPPT * SDAS: * UMSPUT * OTHERS: * None * * History: * -------- * Version Date Author Description * 1 10-15-85 J.-C. Hsu Design and coding * 2 07-15-87 J.-C. Hsu F77 standard * 3 12-01-87 J.-C. Hsu avoid IMSL routines * *------------------------------------------------------------------------------- * * === local: * --maximum number of order and points INTEGER NPMAX, NMAX PARAMETER (NPMAX = 9) PARAMETER (NMAX = 2000) * --the quantity of y-axis REAL Y(NMAX), * --standard deviation of y : YERR(NMAX), * --the quantity of x-axis : X(NMAX), * --lower limit of x : XMIN, * --upper limit of x : XMAX, * --reference value of x : XREF, * --weighting flag : WFLAG, * --attributes : VOLT, GAIN, THRESH * --epoch of the observation DOUBLE PRECISION EPOCH(NMAX), * --lower and upper limit of epoch : EPMIN, EPMAX, * --average value of epoch : EPAVE, * --x and y arrays in double precision : XDP(NMAX), YDP(NMAX), * --y standard deviation array in double precision : YERRDP(NMAX), * --orthogonal polynomial construct : ORTHO(0:NPMAX, 0:NPMAX), * --coefficients of the orthogonal * --polynomials and their errors : A(0:NPMAX), SIG(0:NPMAX), * --coefficients of the "regular" * --polynomial and their error matrix : COEFF(0:NPMAX), : ERMTRX(0:NPMAX, 0:NPMAX), * --chi-squared of the fit : CHISQ * --number of successfully read points INTEGER NPTS, * --maximum order of the fitting polynomial : ORDER, * --detector ID : DETID, * --number of attributes to be copied : NATTRB, * --loop indices : I, * --return status : STATUS, STATOK * --data type CHARACTER*7 TYPE CHARACTER*10 APERT * --column name of reference x value CHARACTER*16 XREFKY, * --column names of minimum and maximum x : XMINKY, XMAXKY, * --attributes to be copied from input * --table to output table : ATTRB(5) * --VMS file name of output table CHARACTER*128 OFILE * --error message context CHARACTER*130 CONTXT *=========================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 VTMPGT (X, Y, YERR, EPOCH, ORDER, WFLAG, TYPE, NPTS, : NATTRB, ATTRB, XREF, XMINKY, XMAXKY, XREFKY, : DETID, APERT, VOLT, GAIN, THRESH, : OFILE, STATUS) IF (STATUS .NE. OK .OR. NPTS .GT. NMAX) THEN CONTXT = 'cannot get parameters/data or too many ' : // 'input data points' GO TO 999 END IF * * check if the order is too high/low * IF (ORDER .LT. 0 .OR. ORDER .GT. NPMAX) THEN STATUS = ERRNUM(1) CONTXT = 'improper input polynomial order' GO TO 999 END IF * * check if there is enough data points * IF (NPTS .LT. ORDER+2) THEN STATUS = ERRNUM(2) CONTXT = 'not enough data points to fit the polynomial' GO TO 999 END IF * * offset x values, convert to double precision, and fit the data with * a polynomial * DO 10 I = 1, NPTS XDP(I) = DBLE(X(I) - XREF) YDP(I) = DBLE(Y(I)) YERRDP(I) = DBLE(YERR(I)) 10 CONTINUE * CALL CDPFIT (XDP, YDP, YERRDP, NPTS, WFLAG, ORDER, NPMAX, : ORTHO, A, SIG, CHISQ) * * convert from orthogonal basis to regular basis * CALL CDPBAS (A, SIG, ORDER, NPMAX, ORTHO, COEFF, ERMTRX) * * find the range of X and epoch, and calculate the average epoch * EPAVE = 0 XMIN = X(1) XMAX = X(1) EPMIN = EPOCH(1) EPMAX = EPOCH(1) * DO 20 I = 1, NPTS EPAVE = EPAVE + EPOCH(I) * IF (X(I) .LT. XMIN) THEN XMIN = X(I) ELSE IF (X(I) .GT. XMAX) THEN XMAX = X(I) END IF * IF (EPOCH(I) .LT. EPMIN) THEN EPMIN = EPOCH(I) ELSE IF (EPOCH(I) .GT. EPMAX) THEN EPMAX = EPOCH(I) END IF 20 CONTINUE * EPAVE = EPAVE / DBLE(NPTS) * * write result to output table * CALL VTMPPT (OFILE, ORDER, NPTS, WFLAG, TYPE, : COEFF, CHISQ, : XREF, XMIN, XMAX, EPMIN, EPMAX, EPAVE, : DETID, APERT, VOLT, GAIN, THRESH, : NATTRB, ATTRB, XREFKY, XMINKY, XMAXKY, : STATUS) IF (STATUS .NE. OK) THEN CONTXT = 'cannot write to output table' GO TO 999 END IF * STATUS = OK GO TO 1000 * * write error message * 999 CALL UMSPUT ('VTEMP: ' // CONTXT, DEST, PRIO, STATOK) * 1000 RETURN END