SUBROUTINE VPEFBN ( * * inputs * : CNT, CNTSD, APPA, TEMP, EPOCH, NPTS, * * inputs/outputs * : APERT, TARGET, * * outputs * : PEFF, PADIF, INTEN, DPEFF, DPADIF, : DINTEN, TMPAVE, EPAVE, NROWS, STATUS) * * Module Number: 15.12.2.1 * * Module Name: poleffv * * Keyphrase: * ---------- * Calculate the polarization efficiency for each aperture. * * Description: * ------------ * Sort input data to different polarizer apertures, and calculate * polarization efficiency for each aperture. * Maximum number of standard targets is 250 (i. e. total number of entries in * the standard target table can not be more than 1000, since each target has * four entries corresponding to each filter). * * FORTRAN Name: VPEFBN.FOR * * Keywords of Accessed Files and Tables: * -------------------------------------- * Name I/O Description / Comments * * Subroutines Called: * ------------------- * CDBS: * CDPTRG, CDROLL, CDP3AN * SDAS: * UMSPUT * OTHERS: * None * * History: * -------- * Version Date Author Description * 1 07-15-86 J.-C. Hsu Design and coding * 2 10-05-87 J.-C. Hsu F77 standard * 3 11-08-89 J.-C. Hsu Error propagation * 4 12-18-89 J.-C. Hsu change V3 position angle to aperture * position angle *------------------------------------------------------------------------------- * *== input: * --count rate and its standard deviation REAL CNT(1), CNTSD(1), * --position angle of the aperture used : APPA(1), * --temperature : TEMP(1) * --epoch DOUBLE PRECISION EPOCH(1) * --number of input data points INTEGER NPTS * *== inputs/outputs * --aperture name CHARACTER*(*) APERT(1), * --target name : TARGET(1) * *== output: * --intensity, polarization and its * --position angle REAL INTEN(1), PEFF(1), PADIF(1), * --error of the above : DINTEN(1), DPEFF(1), DPADIF(1), * --average temperature : TMPAVE(1) * --average epoch DOUBLE PRECISION EPAVE(1) * --number of rows to be put in output table INTEGER NROWS, * --error status : STATUS * *== local: * --loop indices INTEGER I, J, K, M, N, * --number of entries in standard * --polarization target table : NPOL, * --number of different apertures : NAPER, * --number of roll angles : NANGLE(20), * --number of observation at a specific * --aperture : OBAPER(20), * --number of observation at a * --specific roll angle of one aperture : OBANG(20,3), * --aperture index : INDEX1, * --roll angle index : INDEX2, * --number of observations index : INDEX3, : INDEX(20, 3, 2000), * --roll angle indices : ROLLID(3) * --Stokes parameters REAL Q, U, * --polarization and position angle read * --from standard target table : REFP(1000), REFPA(1000), * --errors of the above : DREFP(1000), DREFPA(1000), * --standard polarization and position * --angle for a specific target : P0, PA0, * --sum of count rate, weights : SUM(3), SUMWT(3), * --weights : WT, * --roll angles : ROLL(3), APPAX(20, 3), THETA1, THETA2 INTEGER STAT(10), STATOK * --filter name CHARACTER*4 FILTER, * --filter name of the standard stars : REFFLT(1000) CHARACTER*5 CHAR5 * --aperture name CHARACTER*10 APERX(20) * --target name CHARACTER*20 OBJ(20) * --reference target name CHARACTER*20 REFOBJ(5,1000) * --error message 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=========================================== *------------------------------------------------------------------------------- * * initialize * NROWS = 0 DO 20 J = 1, 20 NANGLE(J) = 0 OBAPER(J) = 0 TMPAVE(J) = 0. EPAVE(J) = 0.D0 DO 10 K = 1, 3 OBANG(J, K) = 0 10 CONTINUE 20 CONTINUE * * sort input data points according to their aperture: * index1 (and j) is the aperture index, index2 (and k) is the roll angle index * DO 70 I = 1, NPTS * * initialize for the first data point * IF (I .EQ. 1) THEN NAPER = 1 APERX(NAPER) = APERT(I) END IF * * match with known apertures * DO 30 J = 1, NAPER IF (APERT(I) .EQ. APERX(J)) THEN INDEX1 = J GO TO 40 END IF 30 CONTINUE * * if no matching aperture, increment the aperture number counter * NAPER = NAPER + 1 APERX(NAPER) = APERT(I) INDEX1 = NAPER * * increment number of observations for the aperture * 40 OBAPER(INDEX1) = OBAPER(INDEX1) + 1 * * assign the target name as the reference object name for the new aperture * IF (OBAPER(INDEX1) .EQ. 1) THEN OBJ(INDEX1) = TARGET(I) END IF * * check if the target name is the same for the same aperture, if not, * skip the data point and write warning message * IF (TARGET(I) .NE. OBJ(INDEX1)) THEN WRITE (CHAR5, '(I5)') I CONTXT = 'VPEFBN: different target for the same ' : // 'aperture at entry #' // CHAR5 CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) OBAPER(J) = OBAPER(J) - 1 ELSE * * sort input data points according to their roll angle * * initialize for the first entry of the specified roll angle * IF (NANGLE(INDEX1) .EQ. 0) THEN NANGLE(INDEX1) = 1 APPAX(INDEX1, NANGLE(INDEX1)) = APPA(I) END IF * * match with know roll angles * DO 50 K = 1, NANGLE(INDEX1) IF (APPA(I) .EQ. APPAX(INDEX1, K)) THEN INDEX2 = K GO TO 60 END IF 50 CONTINUE * * if no matching roll angle, increment the roll angle number counter * NANGLE(INDEX1) = NANGLE(INDEX1) + 1 APPAX(INDEX1, NANGLE(INDEX1)) = APPA(I) INDEX2 = NANGLE(INDEX1) * * check if within maximum number (=3) of roll angles: * 60 IF (INDEX2 .GT. 3) THEN CONTXT = 'VPEFBN: ' // APERT(I) // ' has more than ' : // '3 roll angles' CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) ELSE OBANG(INDEX1, INDEX2) = OBANG(INDEX1, INDEX2) + 1 INDEX3 = OBANG(INDEX1, INDEX2) INDEX(INDEX1, INDEX2, INDEX3) = I END IF END IF 70 CONTINUE * * obtain reference data from standard polarization target table * CALL CDPTRG ( REFOBJ, REFFLT, REFP, DREFP, REFPA, DREFPA, : NPOL, STATUS) IF (STATUS .NE. OK .OR. NPOL .GT. 1000) THEN CONTXT = 'error reading standard polarization target table ' : // 'or too many standard targets (>1000)' GO TO 999 END IF * * process each aperture * DO 200 J = 1, NAPER * * identify filter name (assume filter name has ten characters) corresponding * to the aperture * DO 80 I = 1, 8 IF (APERX(J)(I:I) .EQ. 'F') GO TO 90 80 CONTINUE * STATUS = ERRNUM(1) CONTXT = 'illegal aperture name ' // APERX(J) GO TO 999 * 90 FILTER = APERX(J)(I:I+3) * * match with standard polarization target data * DO 110 I = 1, NPOL IF (FILTER .EQ. REFFLT(I)) THEN DO 100 N = 1, 5 IF (OBJ(J) .EQ. REFOBJ(N, I)) THEN P0 = REFP(I) PA0 = REFPA(I) GO TO 120 END IF 100 CONTINUE END IF 110 CONTINUE * * if no matching standard star data, skip, and write nasty message * CONTXT = 'VPEFBN: no reference polarization for ' // OBJ(J) : // ' at aperture ' // APERX(J) CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 200 * * if data only obtained in less than 3 roll angles, skip, and write * error message * 120 IF (NANGLE(J) .LT. 3) THEN CONTXT = 'VPEFBN: Aperture ' // APERX(J) // : ' has less than three roll angles' CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 200 END IF * * process each roll angle, if there are more than 3 roll angles, only * the first 3 will be processed * DO 140 K = 1, 3 SUM(K) = 0 SUMWT(K) = 0 DO 130 M = 1, OBANG(J, K) WT = 1. / CNTSD(INDEX(J, K, M))**2 SUM(K) = CNT(INDEX(J, K, M)) * WT + SUM(K) SUMWT(K) = WT + SUMWT(K) TMPAVE(J) = TMPAVE(J) + TEMP(INDEX(J, K, M)) EPAVE(J) = EPAVE(J) + EPOCH(INDEX(J, K, M)) 130 CONTINUE * * if no observation in the roll angle bin, put error message * SUM(K) = SUM(K) / SUMWT(K) ROLL(K) = APPAX(J, K) 140 CONTINUE * * determine separations among the three roll angles * CALL CDROLL (ROLL, : ROLLID, THETA1, THETA2, STAT(1)) IF (STAT(1) .NE. OK) THEN CONTXT = 'VPEFBN: ' // APERX(J) // 'has inadequate ' : // 'roll angles' CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 200 END IF * * calculate Stokes parameters from measurements made at the three roll angles * CALL CDP3AN (SUM(ROLLID(1)), SUM(ROLLID(2)), SUM(ROLLID(3)), : THETA1, THETA2, : INTEN(J), Q, U, STAT(2)) IF (STAT(2) .NE. OK) THEN CONTXT = 'VPEFBN: error calculating Stokes parameters ' : // 'for aperture ' // APERX(J) CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 200 END IF * * mean temperature and epoch * NROWS = NROWS + 1 TMPAVE(NROWS) = TMPAVE(J) / REAL(OBAPER(J)) EPAVE(NROWS) = EPAVE(J) / DBLE(OBAPER(J)) * INTEN(NROWS) = INTEN(J) APERT(NROWS) = APERX(J) TARGET(NROWS) = OBJ(J) DINTEN(NROWS) = 0. * * calculate polarization efficiency * PEFF(NROWS) = SQRT (Q**2 + U**2) PADIF(NROWS) = 0.5 * ATAN2(U, Q) * 57.29577951 * * reference polarization is assumed in PER CENTS * PEFF(NROWS) = PEFF(NROWS) / (P0 * 0.01) DPEFF(NROWS) = 0. * * angle diff = observed angle - reference angle + aperture position angle of * zero roll * PADIF(NROWS) = PADIF(NROWS) - PA0 + ROLL(ROLLID(1)) DPADIF(NROWS) = 0. * * restrict angle difference to be between +90 and -90 degrees * 150 IF (PADIF(NROWS) .GT. 90.) THEN PADIF(NROWS) = PADIF(NROWS) - 180. GO TO 150 END IF * 160 IF (PADIF(NROWS) .LT. -90.) THEN PADIF(NROWS) = PADIF(NROWS) + 180. GO TO 160 END IF * 200 CONTINUE * STATUS = OK GO TO 1000 * 999 MESS = 'VPEFBN: ' // CONTXT CALL UMSPUT (MESS, DEST, PRIO, STATOK) * 1000 RETURN END