SUBROUTINE VEPOCH * * Module Number: 15.1.0.2 * * Module Name: polyepoch * * Keyphrase: * ---------- * Analyze a set of physical quantity as polynomial functions of epoch. * * Description: * ------------ * Fit a set of quantities (e. g. temperature dependence coefficients of * dark count) as polynomial functions of epoch. This routine is usually used * to obtain the 16 polynomial coefficients RSDP needs to calibrate the dark * signal, the pre-amp noise, the CVC offset, and other quantities. We * will hereafter refer to the coefficients to be fitted as temperature * dependence coefficients and they are to be expressed as polynomials of * epoch. However, this routine is NOT limited to such application, as its * flexibility is reflected in the user-defined column names of input and * output tables and degrees of the input and output polynomial functions. * * Minimum order of the fitting polynomial is 1, and maximum is 9, default is 3. * Maximum number of input data points (number of rows in the input table) * is 2000. * Input temperature dependence coefficients are in single precision but are * double precision during the polynomial regression. Epoch is in double * precision throughout. * If the weighting flag, WFALG = 0, equal weights are applied to each * data point, if = -2, the inverses of input data's uncertainty is used as * weights. * * FORTRAN Name: VEPOCH.FOR * * Keywords of Accessed Files and Tables: * -------------------------------------- * Name I/O Description / Comments * * Subroutines Called: * ------------------- * CDBS: * CDPFIT, CDPBAS, VEPGT, VEPPT * 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(4*NMAX), * --standard deviation of y : YERR(4*NMAX), * --coefficients of the fitting polynomial : COEFF(100), * --lower and upper limits of temperature : TMIN(NMAX), TMAX(NMAX), TTMIN, TTMAX, * --base temperature : TEMP0, * --weighting flag : WFLAG, * --attributes : VOLT, GAIN, THRESH * --epoch of the observation DOUBLE PRECISION EPOCH(NMAX), * --lower and upper limits of epoch : EPMIN(NMAX), EPMAX(NMAX), EMIN, EMAX, * --base time : TIME0, * --x and y array in double precision : XDP(NMAX), YDP(NMAX), * --y standard deviation array in double precision : YERRDP(NMAX), * --polynomial and their error matrix : COEF(0:NPMAX), : ERMTRX(0:NPMAX, 0:NPMAX), * --orthogonal polynomial construct : ORTHO(0:NPMAX, 0:NPMAX), * --coefficients of the orthogonal * --polynomials and their errors : A(0:NPMAX), SIG(0:NPMAX), * --chi-squared of the fits : CHISQ(10) * --number of input data points INTEGER NROWS, * --number of successfully read points : NPTS, * --maximum order of the fitting polynomial : ORDER, * --number of input coefficients : NCOEFF, * --detector ID : DETID, * --number of attributes to be copied : NATTRB, * --loop indices : I, J, K, * --return status : STATUS, STATOK * --data type CHARACTER*7 TYPE CHARACTER*10 APERT * --attributes to be copied from input * --table to output table CHARACTER*16 ATTRB(5), * --column name of the base value : BASEKY * --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 VEPGT (Y, YERR, EPOCH, TMIN, TMAX, EPMIN, EPMAX, : TEMP0, TIME0, ORDER, WFLAG, TYPE, NPTS, NROWS, : NCOEFF, NATTRB, ATTRB, BASEKY, : 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 * * fit each set of input coefficients * DO 40 J = 1, NCOEFF * * convert Y to double precision, and fit the data with a polynomial * offset epoch with base time * DO 20 I = 1, NPTS K = I + (J-1) * NROWS XDP(I) = EPOCH(I) - TIME0 YDP(I) = DBLE(Y(K)) YERRDP(I) = DBLE(YERR(K)) 20 CONTINUE * CALL CDPFIT (XDP, YDP, YERRDP, NPTS, WFLAG, ORDER, NPMAX, : ORTHO, A, SIG, CHISQ(J)) * * convert from orthogonal basis to regular basis * CALL CDPBAS (A, SIG, ORDER, NPMAX, ORTHO, COEF, ERMTRX) * DO 30 I = 0, ORDER K = I + (J-1) * (ORDER+1) + 1 COEFF(K) = SNGL(COEF(I)) 30 CONTINUE 40 CONTINUE * * find the range of temperature and epoch * TTMIN = TMIN(1) TTMAX = TMAX(1) EMIN = EPMIN(1) EMAX = EPMAX(1) * DO 50 I = 1, NPTS IF (TMIN(I) .LT. TTMIN) THEN TTMIN = TMIN(I) END IF IF (TMAX(I) .GT. TTMAX) THEN TTMAX = TMAX(I) END IF * IF (EPMIN(I) .LT. EMIN) THEN EMIN = EPMIN(I) END IF IF (EPMAX(I) .GT. EMAX) THEN EMAX = EPMAX(I) END IF 50 CONTINUE * * write result to output table * CALL VEPPT (OFILE, ORDER, NPTS, WFLAG, TYPE, : COEFF, CHISQ, BASEKY, : TEMP0, TIME0, TTMIN, TTMAX, EMIN, EMAX, : DETID, APERT, VOLT, GAIN, THRESH, : NATTRB, ATTRB, 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 ('VEPOCH: ' // CONTXT, DEST, PRIO, STATOK) * 1000 RETURN END