SUBROUTINE VRSEN * * Module Number: 15.8 * * Module Name: relsen * * Keyphrase: * ---------- * Calculate the comparative sensitivities relative to the reference aperture. * * Description: * ------------ * Calculate the sensitivities of apertures with the same filter * relative to the "reference aperture", i. e. count rate relative to that * obtained at the reference aperture. Input data must have the * same target, use the same filter, and roughly obtained at the same * temperature/epoch. * Maximum number of input data points is 2000. * * FORTRAN Name: VRSEN.FOR * * Keywords of Accessed Files : * -------------------------- * Name I/O Description / Comments * * Modules Called: * --------------- * CDBS: * VRSNGT, VRSNPT, VREFAP * SDAS: * UMSPUT * OTHERS: * none * * History: * -------- * Version Date Author Description * 1 05-30-85 J.-C. Hsu Design and coding * 2 10-05-87 J.-C. Hsu F77 standard * 3 10-18-89 J.-C. Hsu error propagation *------------------------------------------------------------------------------- * *== local: * --size of input data arrays INTEGER SIZE PARAMETER (SIZE = 2000) * --observed quantity (e. g. count rate) * --and its standard deviation REAL COUNT(SIZE), DCOUNT(SIZE), * --relative sensitivity and its standard * --deviation : RSEN(SIZE), DRSEN(SIZE), * --temperature : TEMP(SIZE), * --summation of data taken at reference * --aperture : SUM, SUMWT, SUMSQ, WT, REFCNT, ERR * --epoch of the observation DOUBLE PRECISION EPOCH(SIZE) * --number of input data points INTEGER NPTS, * --number of output points : NP, NREF, * --indices corresponding to rsen : INDEX(SIZE), * --error status : STATUS, STATOK, * --loop index : I, J * --point source flag CHARACTER*1 PSFLAG * --reference filter name CHARACTER*4 REFFLT(100), * --filter name of this data set : FILTER CHARACTER*5 CHAR5 * --data type CHARACTER*7 TYPE * --reference aperture names CHARACTER*10 REFAP(100), APERX, * --aperture names of input data : APERT(SIZE) * --column name of the temperature CHARACTER*16 TNAME * --target names of input data CHARACTER*20 TRGT * --output table name CHARACTER*128 OFILE * --error message context CHARACTER*130 CONTXT, MESS *=========================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 CL parameters and input data * CALL VRSNGT ( COUNT, DCOUNT, TEMP, EPOCH, APERT, TRGT, PSFLAG, : NPTS, OFILE, TYPE, TNAME, STATUS) IF (STATUS .NE. OK .OR. NPTS .GT. SIZE) THEN WRITE (CHAR5, '(I5)') SIZE CONTXT = 'cannot get input parameters/data or more than ' : // CHAR5 // ' input points' GO TO 999 END IF * * read reference aperture table * CALL VREFAP ( REFFLT, REFAP, NREF, STATUS) IF (STATUS .NE. OK) THEN CONTXT = 'error getting reference apertures' GO TO 999 END IF * * identify filter name from first aperture name and get its corresponding * reference aperture name * DO 10 I = 1, 10 IF (APERT(1)(I:I) .EQ. 'F') GO TO 20 10 CONTINUE * STATUS = ERRNUM(1) CONTXT = 'illegal aperture name ' // APERT(1) // ' at the first' : // ' row' GO TO 999 * 20 FILTER = APERT(1)(I:I+3) * DO 30 I = 1, NREF IF (FILTER .EQ. REFFLT(I)) GO TO 40 30 CONTINUE * * average count rates taken at reference aperture and calculate relative * sensitivities * 40 APERX = REFAP(I) * * initialization * SUM = 0. SUMWT = 0. SUMSQ = 0. NREF = 0 * * pick out the measurements made at the reference aperture, calculate * the weighted average count rate and its error * DO 50 I = 1, NPTS IF (APERT(I) .EQ. APERX) THEN WT = 1. / DCOUNT(I)**2 SUMWT = SUMWT + WT SUM = SUM + WT * COUNT(I) SUMSQ = SUMSQ + WT * COUNT(I)**2 NREF = NREF + 1 J = I END IF 50 CONTINUE * * check if there is any data taken at the reference aperture * IF (NREF .LE. 0) THEN STATUS = ERRNUM(2) CONTXT = 'no data at the reference aperture' GO TO 999 END IF * REFCNT = SUM / SUMWT IF (NREF .EQ. 1) THEN ERR = DCOUNT(J) ELSE ERR = SQRT((SUMSQ/SUMWT - REFCNT**2) / : FLOAT(NREF-1)) END IF * * calculate the relative sensitivity of each data point * NP = 0 DO 80 I = 1, NPTS WRITE (CHAR5, '(I5)') I * DO 60 J = 1, 10 IF (APERT(I)(J:J) .EQ. 'F') GO TO 70 60 CONTINUE CONTXT = 'illegal aperture name ' // APERT(I) // ' at row ' : // CHAR5 CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 80 * * check if the observation was made with the correct filter * 70 IF (APERT(I)(J:J+3) .NE. FILTER) THEN CONTXT = 'different filter ' // APERT(I)(J:J+3) // : ' is used at row ' // CHAR5 CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 80 END IF * * skip observations made with the reference aperture * IF (APERT(I) .NE. APERX) THEN NP = NP + 1 INDEX(NP) = I RSEN(NP) = COUNT(I) / REFCNT DRSEN(NP) = RSEN(NP) * SQRT((ERR / REFCNT)**2 + : (DCOUNT(I) / COUNT(I))**2) END IF 80 CONTINUE * * check if there is any data points besides reference aperture * IF (NP .EQ. 0) THEN STATUS = ERRNUM(2) CONTXT = 'no data point besides reference aperture' GO TO 999 END IF * * write the result to the output table * CALL VRSNPT (OFILE, APERT, TRGT, PSFLAG, RSEN, DRSEN, TEMP, : EPOCH, INDEX, NP, TYPE, TNAME, STATUS) IF (STATUS .NE. OK) THEN GO TO 999 END IF * STATUS = OK GO TO 1000 * * write error message * 999 MESS = 'VRSEN: ' // CONTXT CALL UMSPUT (MESS, DEST, PRIO, STATOK) * 1000 RETURN END