SUBROUTINE VPPABN ( * * inputs * : CNT, CNTSD, APPA, TEMP, EPOCH, APERT, : NPTS, * * inputs/outputs * : TARGET, * * outputs * : FLTOUT, POL, PAOFF, INTEN, TMPAVE, : EPAVE, NROWS, STATUS) * * Module Number: 15.12.3.1 * * Module Name: polpav * * Keyphrase: * ---------- * Calculate the position angle offset for each filter * * Description: * ------------ * Sort input data to different filter, and calculate the position angle * offset for each filter. * 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: VPPABN.FOR * * Keywords of Accessed Files and Tables: * -------------------------------------- * Name I/O Description / Comments * * Keywords input from PAR file: * * Subroutines Called: * ------------------- * CDBS: * CDPTRG, CDPEF, CDP4AP * 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 *------------------------------------------------------------------------------- * *== input: * --count rate REAL CNT(1), * --standard deviation of count rate : CNTSD(1), * --position angle of aperture used : APPA(1), * --temperature : TEMP(1) * --epoch DOUBLE PRECISION EPOCH(1) * --aperture name CHARACTER*(*) APERT(1) * --number of input data points INTEGER NPTS * *== input/output: * --target name CHARACTER*(*) TARGET(1) * *== output: * --polarization REAL POL(1), * --position angle offset : PAOFF(1), * --intensity of target : INTEN(1), * --mean temperature : TMPAVE(1) * --mean epoch DOUBLE PRECISION EPAVE(1) * --number of output rows INTEGER NROWS, * --error status : STATUS * *== local: * --loop indices INTEGER I, J, K, M, N, * --number of entries in polarization * --target table : NPOL, * --number of entries in polarization * --efficiency table : NPEFF, * --number of different filters : NFLT, * --number of aperutres : NAPER(4), * --number of observation at a specific * --filter : OBFLT(4), * --number of observation at a * --specific aperture : OBAPER(4,4), * --filter index : INDEX1, * --aperture index : INDEX2, * --number of observation index : INDEX3, : INDEX(4, 4, 2000), : STAT(10), STATOK * --Stokes parameters REAL Q, U, * --polarization efficiency : PEFF(16), DPEFF(16), * --reference polarization and position angle : REFP(1000), REFPA(1000), * --errors of the above : DREFP(1000), DREFPA(1000), * --reference position angle : PA0, * --sum of count rates : SUM(4), * --sum of the weights : SUMWT(4), * --weight : WT, * --reference roll angles : APPAX(4) * --reference filter name CHARACTER*4 FILTEX(4), * --filter name of each individual observation : FILTER, : FLTOUT(4), * --filter name of standard stars : REFFLT(1000) CHARACTER*5 CHAR5 * --aperture name CHARACTER*10 APERX(16) * --reference target name CHARACTER*20 OBJ(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 NFLT = 0 DO 20 J = 1, 4 NAPER(J) = 0 OBFLT(J) = 0 TMPAVE(J) = 0. EPAVE(J) = 0.D0 * DO 10 K = 1, 4 OBAPER(J, K) = 0 10 CONTINUE 20 CONTINUE * DO 90 I = 1, NPTS WRITE (CHAR5, '(I5)') I * * indentify filter name (assume filter name has ten characters) corresponding * to the aperture * DO 30 M = 1, 7 IF (APERT(I)(M:M) .EQ. 'F') GO TO 40 30 CONTINUE * CONTXT = 'illegal aperture name at entry #' // CHAR5 CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) * 40 FILTER = APERT(I)(M:M+3) * * initialize for the first filter * IF (NFLT .EQ. 0) THEN NFLT = 1 FILTEX(NFLT) = FILTER END IF * * sort input data points based on filter: * index1 (and j) is the filter index, index2 (and k) is the aperture index * DO 50 J = 1, NFLT IF (FILTER .EQ. FILTEX(J)) THEN INDEX1 = J GO TO 60 END IF 50 CONTINUE * * if no matching filter, increment the filter number counter * NFLT = NFLT + 1 FILTEX(NFLT) = FILTER INDEX1 = NFLT * * increment number of observations for the filter * 60 OBFLT(INDEX1) = OBFLT(INDEX1) + 1 * * assign the target name as the reference object name for the new filter * also assign the roll angle as the reference roll angle for the new filter * IF (OBFLT(INDEX1) .EQ. 1) THEN OBJ(INDEX1) = TARGET(I) APPAX(INDEX1) = APPA(I) END IF * * check if the target name and roll angle is the same for the same filter, * if not, skip the data point and write warning message * IF (TARGET(I) .NE. OBJ(INDEX1) .OR. : APPA(I) .NE. APPAX(INDEX1)) THEN CONTXT = 'different targets/roll angles for the ' : // 'same filter at entry #' // CHAR5 CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) OBFLT(J) = OBFLT(J) - 1 GO TO 90 END IF * * sort input data points based on the aperture name * DO 70 M = 1, 10 IF (APERT(I)(M:M) .EQ. 'P') GO TO 80 70 CONTINUE * * flag illegal aperture name * CONTXT = 'illegal aperture name at entry #' // CHAR5 CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 90 * * identify polarizer name * 80 IF (APERT(I)(M+1:M+1) .EQ. '0') THEN INDEX2 = 1 ELSE IF (APERT(I)(M+1:M+2) .EQ. '90') THEN INDEX2 = 2 ELSE IF (APERT(I)(M+1:M+2) .EQ. '45') THEN INDEX2 = 3 ELSE IF (APERT(I)(M+1:M+3) .EQ. '135') THEN INDEX2 = 4 ELSE CONTXT = 'illegal aperture name at entry #' // CHAR5 CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 90 END IF * * if aperture name matches, increment the aperture number counter * OBAPER(INDEX1, INDEX2) = OBAPER(INDEX1, INDEX2) + 1 IF (OBAPER(INDEX1, INDEX2) .EQ. 1) THEN NAPER(INDEX1) = NAPER(INDEX1) + 1 END IF INDEX3 = OBAPER(INDEX1, INDEX2) INDEX(INDEX1, INDEX2, INDEX3) = I 90 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 * * read polarization efficiency table * CALL CDPEF ( APERX, PEFF, DPEFF, NPEFF, STATUS) IF (STATUS .NE. OK .OR. NPEFF .GT. 16) THEN CONTXT = 'error reading polarization efficiency table ' : // 'or it has more than 16 entries' STATUS = ERRNUM(1) GO TO 999 END IF * * process each filter * DO 200 J = 1, NFLT * * match with standard polarization target reference data * DO 110 I = 1, NPOL IF (FILTEX(J) .EQ. REFFLT(I)) THEN DO 100 N = 1, 5 IF (OBJ(J) .EQ. REFOBJ(N, I)) THEN PA0 = REFPA(I) GO TO 120 END IF 100 CONTINUE END IF 110 CONTINUE * * if no matching standard star data, skip, and write warning message * CONTXT = 'no reference polarization for ' // OBJ(J) CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 200 * * if data only obtained in less than 4 apertures, skip, and write * warning message * 120 IF (NAPER(J) .LT. 4) THEN CONTXT = FILTEX(J) // ' does not has data in all four ' : // 'apertures - cannot calculate polarization' CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 200 END IF * * process each aperture * DO 140 K = 1, 4 SUM(K) = 0 SUMWT(K) = 0 * DO 130 M = 1, OBAPER(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 SUM(K) = SUM(K) / SUMWT(K) 140 CONTINUE * * calculate Stokes parameters from measurements made at the 4 apertures of * the same filter * CALL CDP4AP (SUM, FILTEX(J), APERX, PEFF, NPEFF, : INTEN(J), Q, U, STAT(1)) IF (STAT(1) .NE. OK) THEN CONTXT = 'error calculating Stokes parameters ' : // 'for filter ' // FILTEX(J) CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 200 END IF * * calculate mean temperature and epoch * NROWS = NROWS + 1 TMPAVE(NROWS) = TMPAVE(J) / REAL(OBFLT(J)) EPAVE(NROWS) = EPAVE(J) / DBLE(OBFLT(J)) * INTEN(NROWS) = INTEN(J) FLTOUT(NROWS) = FILTEX(J) TARGET(NROWS) = OBJ(J) * * calculate polarization and position angle * POL(NROWS) = SQRT (Q**2 + U**2) PAOFF(NROWS) = 0.5 * ATAN2(U, Q) * 57.29577951 * * change p's unit to per cent * POL(NROWS) = POL(NROWS) * 100. PAOFF(NROWS) = PA0 - PAOFF(NROWS) - APPAX(J) * * restrict position angle offset to be between +90 and -90 degrees * 150 IF (PAOFF(NROWS) .GT. 90.) THEN PAOFF(NROWS) = PAOFF(NROWS) - 180. GO TO 150 END IF * 160 IF (PAOFF(NROWS) .LT. -90.) THEN PAOFF(NROWS) = PAOFF(NROWS) + 180. GO TO 160 END IF * 200 CONTINUE * STATUS = OK GO TO 1000 * 999 MESS = 'VPOLPA: ' // CONTXT CALL UMSPUT (MESS, DEST, PRIO, STATOK) * 1000 RETURN END