SUBROUTINE VFLEX * * Module Number: * * Module Name: flex * * Keyphrase: * ---------- * Analyze the deflection offset as a polynomial function of temperatures * due to temperature differences caused optical bench flexure ("banana" mode * correction). * * Description: * ------------ * (refere to 4/5/1988 Memo by Rick White) * The deflection corrections are calculated by the application processor (AP) * using the following equations: * * T = c1 * T1 + c2 * T2 * dx = a0 + a1 * T + a2 * T^2 + a3 * T^3 * dy = b0 + b1 * T + b2 * T^2 + b3 * T^3 * * where T1 abd T2 are two temperature from HSP telemetry and are specified in * the parameter file, dx and dy are measured deflection offsets in X and Y * directions respectively. c1 and c2 are free parameters which can be * specified in the parameter file, they represent the relative importance * between T1 and T2; a's and b's are the coefficients this task is trying to * calculate. * The order of the fitting polynomial can be any number 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: VFLEX.FOR * * Keywords of Accessed Files and Tables: * -------------------------------------- * Name I/O Description / Comments * * Subroutines Called: * ------------------- * CDBS: * CDPFIT, CDPBAS, VFLXGT, VFLXPT * SDAS: * UMSPUT * OTHERS: * None * * History: * -------- * Version Date Author Description * 1 10-01-88 J.-C. Hsu Design and coding * *------------------------------------------------------------------------------- * * === local: * --maximum number of order and points INTEGER NPMAX, NMAX PARAMETER (NPMAX = 9) PARAMETER (NMAX = 2000) * --deflextions in x and y REAL DX(NMAX), DY(NMAX), * --standard deviation of DX, and DY : DXERR(NMAX), DYERR(NMAX), * --temperatures : T1(NMAX), T2(NMAX), * --temperature parameters : C1, C2, * --lower limit of T1 an T2 : T1MIN, T2MIN, * --upper limit of T1 and T2 : T1MAX, T2MAX, * --nominal X and Y deflections : X0, Y0, * --weighting flag : WFLAG * --epoch of the observation DOUBLE PRECISION EPOCH(NMAX), * --lower and upper limit of epoch : EPMIN, EPMAX, * --temperature in double precision : TDP(NMAX), * --DX and DY arrays in double precision : DXDP(NMAX), DYDP(NMAX), * --DXERR, and DYERR arrays in double precision : DXERDP(NMAX), DYERDP(NMAX), * --orthogonal polynomial construct : ORTHOA(0:NPMAX, 0:NPMAX), : ORTHOB(0:NPMAX, 0:NPMAX), * --coefficients of the orthogonal * --polynomials and their errors : A(0:NPMAX), SIGA(0:NPMAX), : B(0:NPMAX), SIGB(0:NPMAX), * --coefficients of the "regular" * --polynomial and their error matrix : COEFFA(0:NPMAX), MTRXA(0:NPMAX,0:NPMAX), : COEFFB(0:NPMAX), MTRXB(0:NPMAX,0:NPMAX), * --reduced chi-squared of the fit : CHISQA, CHISQB * --number of successfully read points INTEGER NPTS, * --maximum order of the fitting polynomial : ORDER, * --loop indices : I, * --return status : STATUS, STATOK CHARACTER*5 CHAR5 * --aperture name CHARACTER*10 APERT * --column name of reference x value CHARACTER*16 X0NAME, Y0NAME, * --column names of T1 and T2 : T1NAME, T2NAME * --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 VFLXGT (T1, T2, DX, DY, DXERR, DYERR, EPOCH, NPTS, : X0, Y0, C1, C2, APERT, ORDER, WFLAG, : T1NAME, T2NAME, X0NAME, Y0NAME, 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) WRITE (CHAR5, '(I5)') ORDER CONTXT = 'improper input polynomial order ' // CHAR5 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 * * convert deflection data to double precision * DO 10 I = 1, NPTS DXDP(I) = DBLE(DX(I)) DYDP(I) = DBLE(DY(I)) DXERDP(I) = DBLE(DXERR(I)) DYERDP(I) = DBLE(DYERR(I)) 10 CONTINUE * * find the range of T1, T2, and epoch, and calculate the average epoch * T1MIN = T1(1) T1MAX = T1(1) T2MIN = T2(1) T2MAX = T2(1) EPMIN = EPOCH(1) EPMAX = EPOCH(1) * DO 20 I = 1, NPTS IF (T1(I) .LT. T1MIN) THEN T1MIN = T1(I) ELSE IF (T1(I) .GT. T1MAX) THEN T1MAX = T1(I) END IF * IF (T2(I) .LT. T2MIN) THEN T2MIN = T2(I) ELSE IF (T2(I) .GT. T2MAX) THEN T2MAX = T2(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 * * calculate T = c1 * T1 + c2 * T2 and convert to double precision * DO 30 I = 1, NPTS TDP(I) = DBLE(C1 * T1(I) + C2 * T2(I)) 30 CONTINUE * * polynomial fit * CALL CDPFIT (TDP, DXDP, DXERDP, NPTS, WFLAG, ORDER, NPMAX, : ORTHOA, A, SIGA, CHISQA) CALL CDPFIT (TDP, DYDP, DYERDP, NPTS, WFLAG, ORDER, NPMAX, : ORTHOB, B, SIGB, CHISQB) * * convert from orthogonal basis to regular basis * CALL CDPBAS (A, SIGA, ORDER, NPMAX, ORTHOA, COEFFA, MTRXA) CALL CDPBAS (B, SIGB, ORDER, NPMAX, ORTHOB, COEFFB, MTRXB) * * write result to output table * CALL VFLXPT (OFILE, ORDER, NPTS, WFLAG, : COEFFA, CHISQA, COEFFB, CHISQB, : C1, C2, X0, Y0, X0NAME, Y0NAME, APERT, : T1MIN, T1MAX, T2MIN, T2MAX, T1NAME, T2NAME, : EPMIN, EPMAX, 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 ('VFLEX: ' // CONTXT, DEST, PRIO, STATOK) * 1000 RETURN END