C C Subroutines for the ST-ECF impol package. C Written by J. R. Walsh, ST-ECF, ESO (jwalsh@eso.org) December 1999 C Version 1.00 C SUBROUTINE INSTHEAD(N,IMDSCR,INSTNAME,FILTNAME,THETA,REFANG, : ORIENT,CORANG,IN,MTHETA,IMTHETA,LSTAT) C C This subroutine searches the keywords in the header of an C image whose descriptor is IMDSCR for the instrument name. C The expected instrument names are FOC, WFPC2, NICMOS, ACS C and SPECIAL. Depending on the instrument name the image C header is then searched for the values of the following: C the name of the colour filter (FILTNAME), C the polarizer angle (THETA) given either by the C polarizing filter name or for POLQ, the detector number C the correction to THETA (REFANG) for angles relative to C the s-axis of the pick-off mirror C the orientation angle ORIENT (from PA_V3) of the C instrument C the correction (CORANG) to ORIENT for angles relative C to the s-axis of the pick-off mirror C the polarizer angle (MTHETA) - only required for FOC and C ACS (IMTHETA 2 and 4 respectively) C integer indicator (IMTHETA) to indicate use of THETA by C instrument C and the IN'th value of the output arrays are set. If an error C occurs either in instrument name or keyword reading LSTAT is C set false C IMPLICIT NONE INTEGER N ! No. of elements in arrays INTEGER IMDSCR ! Image descriptor INTEGER IN ! Index of arrays to operate on INTEGER IMTHETA ! Flag to indicate cases for THETA setting: 1=WFPC2; 2=FOC; 3=NICMOS; 4=ACS REAL THETA(N) ! The position angle of the polarizer REAL REFANG(N) ! Reference angle to subtract from THETA REAL ORIENT(N) ! The position angle of the telescope/instrument REAL CORANG(N) ! Additive angle correction to ORIENT REAL MTHETA ! The position angle of the polarizer CHARACTER*16 INSTNAME ! Name of instrument CHARACTER*16 FILTNAME ! Name of colour filter LOGICAL LSTAT ! Logical status return (=FALSE if error) C C Depending on the instrument read the polarization parameters C from the header C IF (INSTNAME(:3).EQ.'ACS') THEN CALL SET_ACS(IMDSCR,FILTNAME,THETA(IN),REFANG(IN), : ORIENT(IN),CORANG(IN),LSTAT) MTHETA=THETA(IN) IMTHETA=4 ENDIF IF (INSTNAME(:3).EQ.'FOC') THEN CALL SET_FOC(IMDSCR,FILTNAME,THETA(IN),REFANG(IN), : ORIENT(IN),CORANG(IN),LSTAT) MTHETA=THETA(IN) IMTHETA=2 ENDIF IF (INSTNAME(:6).EQ.'NICMOS') THEN CALL SET_NICMOS(IMDSCR,FILTNAME,THETA(IN),REFANG(IN), : ORIENT(IN),CORANG(IN),LSTAT) IMTHETA=3 ENDIF IF (INSTNAME(:5).EQ.'WFPC2') THEN CALL SET_WFPC2(IMDSCR,FILTNAME,THETA(IN),REFANG(IN), : ORIENT(IN),CORANG(IN),LSTAT) IMTHETA=1 ENDIF IF (INSTNAME(:7).EQ.'SPECIAL') THEN CALL SET_SPECIAL(IMDSCR,FILTNAME,THETA(IN),REFANG(IN), : ORIENT(IN),CORANG(IN),LSTAT) MTHETA=THETA(IN) IMTHETA=5 ENDIF 99 END SUBROUTINE SET_ACS(IMDSCR,FILTNAME,THETA,REFANG, : ORIENT,CORANG,LSTAT) C C This subroutine searches the keywords in the header of a C ACS image for: C the name of the colour filter (FILTNAME), C the polarizer angle (THETA) given by the C polarizing filter name C the correction to THETA (REFANG) for angles relative to C the s-axis of the pick-off mirror C the orientation angle ORIENT (from PA_V3) of the C instrument C the correction (CORANG) to ORIENT for angles relative C to the s-axis of the pick-off mirror C and sets the output parameters. If an error occurs, C then LSTAT is set false C IMPLICIT NONE INTEGER IMDSCR ! Image descriptor REAL THETA ! The position angle of the polarizer REAL REFANG ! Reference angle to subtract from THETA REAL ORIENT ! The position angle of the telescope/instrument REAL CORANG ! Additive angle correction to ORIENT CHARACTER*16 FILTNAME ! Name of colour filter LOGICAL LSTAT ! Logical status return (=FALSE if error) C C Local variables C INTEGER STAT INTEGER STATO INTEGER ILEN1 INTEGER ILEN2 CHARACTER*16 CBUFF1 CHARACTER*16 CBUFF2 CHARACTER*16 CBUFF3 CHARACTER*80 TEXT LSTAT=.FALSE. C C Set correction angles for ACS C REFANG=135.0 ! Need to be checked ?? CORANG=90.0 ! Need to be checked ?? C C Get first filtername C CALL UHDGST(IMDSCR,'FILTNAM1',CBUFF1,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,21) 21 FORMAT(' ERROR READING FILTER NAME') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. GO TO 99 ENDIF IF (CBUFF1(:1).EQ.'F') THEN ILEN1=INDEX(CBUFF1,' ') - 1 FILTNAME=CBUFF1(:ILEN1) LSTAT=.TRUE. ENDIF IF (CBUFF1(:3).EQ.'POL') THEN ILEN2=INDEX(CBUFF1,' ') -1 CBUFF3=CBUFF1(:ILEN2) LSTAT=.TRUE. ENDIF C C Get second filtername C CALL UHDGST(IMDSCR,'FILTNAM2',CBUFF2,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,21) CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. GO TO 99 ENDIF IF (CBUFF2(:1).EQ.'F') THEN ILEN1=INDEX(CBUFF2,' ') - 1 FILTNAME=CBUFF2(:ILEN1) LSTAT=.TRUE. ENDIF IF (CBUFF2(:3).EQ.'POL') THEN ILEN2=INDEX(CBUFF2,' ') -1 CBUFF3=CBUFF2(:ILEN2) LSTAT=.TRUE. ENDIF C C Directly convert polarizer filter name to polarizer angle C IF (LSTAT) THEN IF (CBUFF3(:4).EQ.'POL0') THEN THETA=0.0 ENDIF IF (CBUFF3(:5).EQ.'POL60') THEN THETA=60.0 ENDIF IF (CBUFF3(:6).EQ.'POL120') THEN THETA=120.0 ENDIF ENDIF C C Get orientation angle (ORIENTAT) C CALL UHDGSR(IMDSCR,'ORIENTAT',ORIENT,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,41) 41 FORMAT(' ERROR READING PA_V3 ANGLE') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. GO TO 99 ENDIF 99 END SUBROUTINE SET_FOC(IMDSCR,FILTNAME,THETA,REFANG, : ORIENT,CORANG,LSTAT) C C This subroutine searches the keywords in the header of a C FOC image for: C the name of the colour filter (FILTNAME), C the polarizer angle (THETA) given by the C polarizing filter name C the orientation angle ORIENT (from PA_V3) of the C instrument C and sets the output parameters. If an error occurs, C then LSTAT is set false C IMPLICIT NONE INTEGER IMDSCR ! Image descriptor REAL THETA ! The position angle of the polarizer REAL REFANG ! Reference angle to subtract from THETA REAL ORIENT ! The position angle of the telescope/instrument REAL CORANG ! Additive angle correction to ORIENT CHARACTER*16 FILTNAME ! Name of colour filter LOGICAL LSTAT ! Logical status return (=FALSE if error) C C Local variables C INTEGER STAT INTEGER STATO INTEGER ILEN1 INTEGER ILEN2 CHARACTER*16 CBUFF1 CHARACTER*16 CBUFF2 CHARACTER*16 CBUFF3 CHARACTER*80 TEXT LSTAT=.FALSE. C C Set correction angles for FOC C REFANG=0.0 CORANG=90.0 ! Required to produce the correct polarization position angles C C Get first filtername C CALL UHDGST(IMDSCR,'FILTNAM1',CBUFF1,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,21) 21 FORMAT(' ERROR READING FILTER NAME') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. GO TO 99 ENDIF IF (CBUFF1(:1).EQ.'F') THEN ILEN1=INDEX(CBUFF1,' ') - 1 FILTNAME=CBUFF1(:ILEN1) LSTAT=.TRUE. ENDIF IF (CBUFF1(:3).EQ.'POL') THEN ILEN2=INDEX(CBUFF1,' ') -1 CBUFF3=CBUFF1(:ILEN2) LSTAT=.TRUE. ENDIF C C Get second filtername C CALL UHDGST(IMDSCR,'FILTNAM2',CBUFF2,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,21) CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. GO TO 99 ENDIF IF (CBUFF2(:1).EQ.'F') THEN ILEN1=INDEX(CBUFF2,' ') - 1 FILTNAME=CBUFF2(:ILEN1) LSTAT=.TRUE. ENDIF IF (CBUFF2(:3).EQ.'POL') THEN ILEN2=INDEX(CBUFF2,' ') -1 CBUFF3=CBUFF2(:ILEN2) LSTAT=.TRUE. ENDIF C C Directly convert polarizer filter name to polarizer angle C IF (LSTAT) THEN IF (CBUFF3(:ILEN2).EQ.'POL0') THEN THETA=0.0 ENDIF IF (CBUFF3(:ILEN2).EQ.'POL60') THEN THETA=60.0 ENDIF IF (CBUFF3(:ILEN2).EQ.'POL120') THEN THETA=120.0 ENDIF ENDIF C C Get orientation angle as position angle of the C image +Y axis (ORIENTAT) C CALL UHDGSR(IMDSCR,'ORIENTAT',ORIENT,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,41) 41 FORMAT(' ERROR READING ORIENTATION ANGLE') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. GO TO 99 ENDIF 99 END SUBROUTINE SET_NICMOS(IMDSCR,FILTNAME,THETA,REFANG, : ORIENT,CORANG,LSTAT) C C This subroutine searches the keywords in the header of a C NICMOS image for: C the name of the polarizer filter (FILTNAME), C the polarizer angle (THETA) given by the C polarizing filter name (not set in this routine) C the polarizer angle (THETA) given by the C polarizing filter name (not set in this routine) C the orientation angle ORIENT (from PA_V3) of the C instrument C and sets the output parameters. If an error occurs, C then LSTAT is set false C IMPLICIT NONE INTEGER IMDSCR ! Image descriptor REAL THETA ! The position angle of the polarizer REAL REFANG ! Reference angle to subtract from THETA REAL ORIENT ! The position angle of the telescope/instrument REAL CORANG ! Additive angle correction to ORIENT CHARACTER*16 FILTNAME ! Name of polarizer filter LOGICAL LSTAT ! Logical status return (=FALSE if error) C C Local variables C INTEGER STAT INTEGER STATO INTEGER ILEN CHARACTER*16 CBUFF CHARACTER*80 TEXT LSTAT=.FALSE. C C Set correction angles for NICMOS C REFANG=0.0 CORANG=0.0 C C Get filtername C CALL UHDGST(IMDSCR,'FILTER',CBUFF,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,21) 21 FORMAT(' ERROR READING FILTER NAME') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. GO TO 99 ENDIF IF (CBUFF(:3).EQ.'POL') THEN ILEN=INDEX(CBUFF,' ') -1 FILTNAME=CBUFF(:ILEN) LSTAT=.TRUE. ELSE LSTAT=.FALSE. ENDIF C C Get orientation angle (call ORIENTAT) C CALL UHDGSR(IMDSCR,'ORIENTAT',ORIENT,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,41) 41 FORMAT(' ERROR READING ORIENTATION ANGLE') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. GO TO 99 ENDIF 99 END SUBROUTINE SET_WFPC2(IMDSCR,FILTNAME,THETA,REFANG, : ORIENT,CORANG,LSTAT) C C This subroutine searches the keywords in the header of a C WFPC2 image for: C the name of the colour filter (FILTNAME), C the polarizer angle (THETA) given either by the C polarizing filter name or for POLQ, the detector number C the correction to THETA (REFANG) for angles relative to C the s-axis of the pick-off mirror C the orientation angle ORIENT (from PA_V3) of the C instrument C the correction (CORANG) to ORIENT for angles relative C to the s-axis of the pick-off mirror C and sets the output parameters. If an error occurs, or C a filter name POLQ is not found, LSTAT is set false C IMPLICIT NONE INTEGER IMDSCR ! Image descriptor REAL THETA ! The position angle of the polarizer REAL REFANG ! Reference angle to subtract from THETA REAL ORIENT ! The position angle of the telescope/instrument REAL CORANG ! Additive angle correction to ORIENT CHARACTER*16 FILTNAME ! Name of colour filter LOGICAL LSTAT ! Logical status return (=FALSE if error) C C Local variables C INTEGER STAT INTEGER STATO INTEGER ILEN INTEGER IBUFF CHARACTER*16 CBUFF1 CHARACTER*16 CBUFF2 CHARACTER*16 CBUFF3 CHARACTER*80 TEXT LSTAT=.FALSE. C C Set correction angles for WFPC2 C REFANG=135.0 CORANG=90.0 C REFANG=0.0 C CORANG=0.0 C call umsput(' WARNING - REFANG and CORANG set to 0',1,0,stat) C C Get first filtername C CALL UHDGST(IMDSCR,'FILTNAM1',CBUFF1,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,21) 21 FORMAT(' ERROR READING FILTER NAME') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. GO TO 99 ENDIF IF (CBUFF1(:1).EQ.'F') THEN ILEN=INDEX(CBUFF1,' ') - 1 FILTNAME=CBUFF1(:ILEN) ENDIF C C Get second filtername C CALL UHDGST(IMDSCR,'FILTNAM2',CBUFF2,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,21) CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. GO TO 99 ENDIF IF (CBUFF2(:1).EQ.'F') THEN ILEN=INDEX(CBUFF2,' ') - 1 FILTNAME=CBUFF2(:ILEN) ENDIF C C Check that POL filter is at least one of the C filters C IF (CBUFF1(:3).EQ.'POL'.OR.CBUFF2(:3).EQ.'POL') THEN LSTAT=.TRUE. ENDIF IF (CBUFF1(:3).EQ.'POL') THEN ILEN=INDEX(CBUFF1,' ') -1 CBUFF3=CBUFF1(:ILEN) ENDIF IF (CBUFF2(:3).EQ.'POL') THEN ILEN=INDEX(CBUFF2,' ') -1 CBUFF3=CBUFF2(:ILEN) ENDIF C C Deduce polarizer angle from polarizing filter name C IF (CBUFF3(:ILEN).EQ.'POLQ') THEN C C Need to get detector number to determine THETA C CALL UHDGSI(IMDSCR,'DETECTOR',IBUFF,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,31) 31 FORMAT(' ERROR READING DETECTOR NUMBER') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. GO TO 99 ENDIF C C Convert detector number to polarizer position angle C IF (IBUFF.EQ.1) THEN THETA=135.0 ! PC ENDIF IF (IBUFF.EQ.2) THEN THETA=0.0 ! WF2 ENDIF IF (IBUFF.EQ.3) THEN THETA=45.0 ! WF3 ENDIF IF (IBUFF.EQ.4) THEN THETA=90.0 ! WF4 ENDIF ELSE C C Directly convert polarizer filter name to polarizer angle C IF (CBUFF3(:ILEN).EQ.'POLQN33') THEN THETA=102.0 ENDIF IF (CBUFF3(:ILEN).EQ.'POLQN18') THEN THETA=117.0 ENDIF IF (CBUFF3(:ILEN).EQ.'POLQP15') THEN THETA=15.0 ENDIF IF (CBUFF3(:ILEN).EQ.'POLQP15P') THEN THETA=15.0 ENDIF IF (CBUFF3(:ILEN).EQ.'POLQP15W') THEN THETA=15.0 ENDIF ENDIF C C Get orientation angle (PA_V3) C CALL UHDGSR(IMDSCR,'PA_V3',ORIENT,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,41) 41 FORMAT(' ERROR READING PA_V3 ANGLE') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. GO TO 99 ENDIF 99 END SUBROUTINE SET_SPECIAL(IMDSCR,FILTNAME,THETA,REFANG, : ORIENT,CORANG,LSTAT) C C This subroutine searches the keywords in the header of a C SPECIAL image for: C the name of the polarizer filter (FILTNAME), C the polarizer angle (THETA) given by the parameter C POLANG C the orientation angle ORIENT (from PA_V3) of the C instrument C and sets the output parameters. If an error occurs, C then LSTAT is set false C IMPLICIT NONE INTEGER IMDSCR ! Image descriptor REAL THETA ! The position angle of the polarizer REAL REFANG ! Reference angle to subtract from THETA REAL ORIENT ! The position angle of the telescope/instrument REAL CORANG ! Additive angle correction to ORIENT CHARACTER*16 FILTNAME ! Name of polarizer filter LOGICAL LSTAT ! Logical status return (=FALSE if error) C C Local variables C INTEGER STAT INTEGER STATO INTEGER ILEN CHARACTER*16 CBUFF CHARACTER*80 TEXT LSTAT=.FALSE. C C Set correction angles for SPECIAL instrument C REFANG=0.0 CORANG=0.0 C C Get filtername C CALL UHDGST(IMDSCR,'FILTER',CBUFF,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,21) 21 FORMAT(' ERROR READING FILTER NAME') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. GO TO 99 ENDIF ILEN=INDEX(CBUFF,' ') -1 FILTNAME=CBUFF(:ILEN) LSTAT=.TRUE. C C Get polarizer angle C CALL UHDGSR(IMDSCR,'POLANG',THETA,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,31) 31 FORMAT(' ERROR READING POLARIZER ANGLE VALUE') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. GO TO 99 ENDIF C C Get orientation angle (ORIENTAT) C CALL UHDGSR(IMDSCR,'ORIENTAT',ORIENT,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,41) 41 FORMAT(' ERROR READING ORIENTATION ANGLE') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. GO TO 99 ENDIF 99 END SUBROUTINE INTOCUB(N1,N2,ARRY,M1,M2,M3,CUB,KS) C C This subroutine places the array ARRY into the C KS'th plane of the 3-D array CUB C IMPLICIT NONE INTEGER N1 INTEGER N2 INTEGER M1 INTEGER M2 INTEGER M3 INTEGER KS REAL ARRY(N1,N2) REAL CUB(M1,M2,M3) C C Local variables C INTEGER I,J C C Copy the array ARRY to the KS'th section of the 3-D C array CUB C DO J=1,N2,1 DO I=1,N1,1 CUB(I,J,KS)=ARRY(I,J) ENDDO ENDDO 99 END SUBROUTINE READ_POLTAB(TBDSCR,FILTNAME,MTHETA,IMTHETA,NR, : AFNAME,APANG,ATPAR,ATPER,ARS,ARP, : ADELPH,AINSPO,AINSPA,LFLAG,M,ON, : THETA,TPAR,TPER,RS,RP,DELPH,INSPO, : INSPA,LSTAT) C C This subroutine reads columns of the table file, whose descriptor C is TBDSCR, into the following arrays: C column FILTER to array AFNAME C column POLANG to array APANG C column TRANSPAR to array ATPAR C column TRANSPER to array ATPER C column REFLECTRS to array ARS C column REFLECTRP to array ARP C column RETARDPHI to array ADELPH C column INSTPOL to array AINSPO C column INSTPA to array AINSPA C A match is then sought between FILTNAME and the list of C filter names AFNAME. If IMTHETA is 2, then the polarizer C angle is also matched with MTHETA (FOC). If a match C occurs the corresponding values are read into the ON'th C values of the output arrays THETA, TPAR, TPER, RS, RP, C DELPH, INSPO and INSPA respectively. C If no match occurs then LSTAT is set FALSE C IMPLICIT NONE INTEGER TBDSCR ! Descriptor of table to be read INTEGER IMTHETA ! Distinguishes cases for matching: 1=WFPC2; 2=FOC; 3=NICMOS; 4=ACS; 5=SPECIAL INTEGER NR ! Number of rows in the table to be read INTEGER M ! First dimension of output arrays INTEGER ON ! Index number of output arrays for writing values REAL MTHETA ! Polarizer angle to match with APANG REAL APANG(NR) ! Array of table values for polarizer angle REAL ATPAR(NR) ! Array of table values for polarizer parallel transmission REAL ATPER(NR) ! Array of table values for polarizer perpendicular transmission REAL ARS(NR) ! Array of table values of mirror reflectance for s-wave (E vector) REAL ARP(NR) ! Array of table values of mirror reflectance for p-wave (M vector) REAL ADELPH(NR) ! Array of table values of mirror retardance REAL AINSPO(NR) ! Array of table values of instrumental polarization REAL AINSPA(NR) ! Array of table values of position angle of instrumental polarization REAL THETA(M) ! The position angle of the polarizer filter REAL TPAR(M) ! The parallel component of polarizer transmission REAL TPER(M) ! The perpendicular component of polarizer transmission REAL RS(M) ! Reflectance of s wave REAL RP(M) ! Reflectance of p wave REAL DELPH(M) ! Retardance of s wave relative to p wave REAL INSPO(M) ! Instrumental polarization (%) REAL INSPA(M) ! Position angle of instrumental polarization (deg.) CHARACTER*16 FILTNAME ! Name of colour/polarizer filter CHARACTER*16 AFNAME(NR) ! Array of table values of name of colour/polarizer filter LOGICAL LFLAG(NR) ! Array of flags for status in reading table column values LOGICAL LSTAT ! True if FILTER matched with AFNAME (and MTHETA with APANG) C C Local variables C INTEGER I INTEGER ILEN INTEGER STAT INTEGER STATO INTEGER COLID1 INTEGER COLID2 INTEGER COLID3 INTEGER COLID4 INTEGER COLID5 INTEGER COLID6 INTEGER COLID7 INTEGER COLID8 INTEGER COLID9 CHARACTER*12 BNAME CHARACTER*80 TEXT LSTAT=.TRUE. C C Match the column names one by one and read the values into C appropriate arrays C C Get column identifier for column FILTER C CALL UTCFND(TBDSCR,'FILTER',1,COLID1,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,21) 21 FORMAT(' Failed to find column named FILTER ') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. ENDIF C C Read all values for this column into array AFNAME C CALL UTCGTT(TBDSCR,COLID1,1,NR,AFNAME,LFLAG,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,22) 22 FORMAT(' Failed to read values from column named FILTER ') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. ENDIF C C Get column identifier for column POLANG C CALL UTCFND(TBDSCR,'POLANG',1,COLID2,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,31) 31 FORMAT(' Failed to find column named POLANG ') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. ENDIF C C Read all values for this column into array APANG C CALL UTCGTR(TBDSCR,COLID2,1,NR,APANG,LFLAG,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,32) 32 FORMAT(' Failed to read values from column named POLANG ') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. ENDIF C C Get column identifier for column TRANSPAR C CALL UTCFND(TBDSCR,'TRANSPAR',1,COLID3,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,41) 41 FORMAT(' Failed to find column named TRANSPAR ') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. ENDIF C C Read all values for this column into array ATPAR C CALL UTCGTR(TBDSCR,COLID3,1,NR,ATPAR,LFLAG,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,42) 42 FORMAT(' Failed to read values from column named TRANSPAR ') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. ENDIF C C Get column identifier for column TRANSPER C CALL UTCFND(TBDSCR,'TRANSPER',1,COLID4,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,51) 51 FORMAT(' Failed to find column named TRANSPER ') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. ENDIF C C Read all values for this column into array ATPER C CALL UTCGTR(TBDSCR,COLID4,1,NR,ATPER,LFLAG,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,52) 52 FORMAT(' Failed to read values from column named TRANSPER ') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. ENDIF C C Get column identifier for column REFLECTRS C CALL UTCFND(TBDSCR,'REFLECTRS',1,COLID5,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,61) 61 FORMAT(' Failed to find column named REFLECTRS ') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. ENDIF C C Read all values for this column into array ATRS C CALL UTCGTR(TBDSCR,COLID5,1,NR,ARS,LFLAG,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,72) 72 FORMAT(' Failed to read values from column named REFLECTRS ') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. ENDIF C C Get column identifier for column REFLECTRP C CALL UTCFND(TBDSCR,'REFLECTRP',1,COLID6,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,81) 81 FORMAT(' Failed to find column named REFLECTRP ') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. ENDIF C C Read all values for this column into array ARP C CALL UTCGTR(TBDSCR,COLID6,1,NR,ARP,LFLAG,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,92) 92 FORMAT(' Failed to read values from column named REFLECTRP ') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. ENDIF C C Get column identifier for column RETARDPHI C CALL UTCFND(TBDSCR,'RETARDPHI',1,COLID7,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,101) 101 FORMAT(' Failed to find column named RETARDPHI ') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. ENDIF C C Read all values for this column into array ADELPH C CALL UTCGTR(TBDSCR,COLID7,1,NR,ADELPH,LFLAG,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,112) 112 FORMAT(' Failed to read values from column named RETARDPHI ') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. ENDIF C C Get column identifier for column INSTPOL C CALL UTCFND(TBDSCR,'INSTPOL',1,COLID8,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,121) 121 FORMAT(' Failed to find column named INSTPOL ') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. ENDIF C C Read all values for this column into array AINSPO C CALL UTCGTR(TBDSCR,COLID8,1,NR,AINSPO,LFLAG,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,132) 132 FORMAT(' Failed to read values from column named INSTPOL ') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. ENDIF C C Get column identifier for column INSTPA C CALL UTCFND(TBDSCR,'INSTPA',1,COLID9,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,141) 141 FORMAT(' Failed to find column named INSTPA ') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. ENDIF C C Read all values for this column into array AINSPO C CALL UTCGTR(TBDSCR,COLID9,1,NR,AINSPA,LFLAG,STAT) IF (STAT.NE.0) THEN WRITE(TEXT,152) 152 FORMAT(' Failed to read values from column named INSTPA ') CALL UMSPUT(TEXT,3,0,STATO) LSTAT=.FALSE. ENDIF C C Attempt to match the filter name FILTNAM with the values C in AFNAME. If match found copy corresponding values in C all arrays to output values. C 200 LSTAT=.FALSE. ILEN=INDEX(FILTNAME,' ') -1 IF (IMTHETA.EQ.1) THEN ! WFPC2 case - dont set THETA DO I=1,NR,1 BNAME=AFNAME(I) IF (FILTNAME(:ILEN).EQ.BNAME(:ILEN)) THEN TPAR(ON)=ATPAR(I) TPER(ON)=ATPER(I) RS(ON)=ARS(I) RP(ON)=ARP(I) DELPH(ON)=ADELPH(I) INSPO(ON)=AINSPO(I) INSPA(ON)=AINSPA(I) LSTAT=.TRUE. ENDIF ENDDO ENDIF IF (IMTHETA.EQ.2) THEN ! FOC case - set THETA to APANG DO I=1,NR,1 BNAME=AFNAME(I) IF (FILTNAME(:ILEN).EQ.BNAME(:ILEN).AND. : MTHETA.EQ.APANG(I)) THEN THETA(ON)=APANG(I) TPAR(ON)=ATPAR(I) TPER(ON)=ATPER(I) RS(ON)=ARS(I) RP(ON)=ARP(I) DELPH(ON)=ADELPH(I) INSPO(ON)=AINSPO(I) INSPA(ON)=AINSPA(I) LSTAT=.TRUE. ENDIF ENDDO ENDIF IF (IMTHETA.EQ.3) THEN ! NICMOS case - set THETA to APANG DO I=1,NR,1 BNAME=AFNAME(I) IF (FILTNAME(:ILEN).EQ.BNAME(:ILEN)) THEN THETA(ON)=APANG(I) TPAR(ON)=ATPAR(I) TPER(ON)=ATPER(I) RS(ON)=ARS(I) RP(ON)=ARP(I) DELPH(ON)=ADELPH(I) INSPO(ON)=AINSPO(I) INSPA(ON)=AINSPA(I) LSTAT=.TRUE. ENDIF ENDDO ENDIF IF (IMTHETA.EQ.4) THEN ! ACS case - set THETA to APANG DO I=1,NR,1 BNAME=AFNAME(I) IF (FILTNAME(:ILEN).EQ.BNAME(:ILEN).AND. : MTHETA.EQ.APANG(I)) THEN THETA(ON)=APANG(I) TPAR(ON)=ATPAR(I) TPER(ON)=ATPER(I) RS(ON)=ARS(I) RP(ON)=ARP(I) DELPH(ON)=ADELPH(I) INSPO(ON)=AINSPO(I) INSPA(ON)=AINSPA(I) LSTAT=.TRUE. ENDIF ENDDO ENDIF IF (IMTHETA.EQ.5) THEN ! SPECIAL case - set THETA TO APANG DO I=1,NR,1 BNAME=AFNAME(I) IF (FILTNAME(:ILEN).EQ.BNAME(:ILEN).AND. : MTHETA.EQ.APANG(I)) THEN THETA(ON)=APANG(I) TPAR(ON)=ATPAR(I) TPER(ON)=ATPER(I) RS(ON)=ARS(I) RP(ON)=ARP(I) DELPH(ON)=ADELPH(I) INSPO(ON)=AINSPO(I) INSPA(ON)=AINSPA(I) LSTAT=.TRUE. ENDIF ENDDO ENDIF 999 END SUBROUTINE POLEFFIC(NV,THETA,REFANG,TPAR,TPER,RS,RP,DELPH, : ORIENT,CORANG,M,N,MATPOL,MATROP,MATMIR, : MATROT,WORK,MATSYS,EFFI,EFFQ,EFFU) C C This subroutine determines the system efficiency for Stokes C I (EFFI), for Stokes Q (EFFQU) and for Stokes U (EFFUQ) by C determining the system matrix of effective instrumental C (de)polarization (see e.g. Biretta & McMaster, 1997, WFPC2 C ISR 97-11). The polarizer rotation matrix is not included C since the cos2*theta and sin2*theta terms are explicitly C solved for in the least sqares procedure C IMPLICIT NONE INTEGER NV ! Number of polarizer elements to calculate INTEGER M ! No. of rows of matrixes MAT... INTEGER N ! No. of columns of matrices MAT... REAL THETA(NV) REAL REFANG(NV) REAL TPAR(NV) ! The parallel component of polarizer transmission REAL TPER(NV) ! The perpendicular component of polarizer transmission REAL RS(NV) ! Reflectance of s wave REAL RP(NV) ! Reflectance of s wave REAL DELPH(NV) ! Retardance of s wave relative to p wave REAL ORIENT(NV) ! The position angle of the telescope/instrument REAL CORANG(NV) ! Additive angle correction REAL MATROP(M,N) REAL MATPOL(M,N) ! Mueller matrix of polarizer REAL MATMIR(M,N) ! Mueller matrix of mirror polarization REAL MATROT(M,N) ! Mueller matrix of PA rotation REAL WORK(M,N) ! Work array for matrix multiplications REAL MATSYS(M,N) ! Mueller matrix of system polarization REAL EFFI(NV) ! Output efficiency for Stokes I detection for NV polarizer positions REAL EFFQ(NV) ! Output efficiency for Stokes Q detection for NV polarizer positions REAL EFFU(NV) ! Output efficiency for Stokes U detection for NV polarizer positions C C Local variables C INTEGER I REAL DELPHI C C Calculate the system efficiency for Stokes I (EFFI), for C Stokes Q (EFFQU) and for Stokes U (EFFUQ) for all number of C polarizer elements C DO I=1,NV,1 IF (DELPH(I).NE.0.0) THEN DELPHI=180.0-DELPH(I) ! Correct Delta Phi for WFPC2 ELSE DELPHI=0.0 ENDIF CALL SYSMATRIX(THETA(I),REFANG(I),TPAR(I),TPER(I),RS(I), : RP(I),DELPHI,ORIENT(I),CORANG(I),M,N,MATPOL, : MATROP,MATMIR,MATROT,WORK,MATSYS) C C Copy system efficiencies for I, Q and U to output arrays C EFFI(I)=MATSYS(1,1) EFFQ(I)=MATSYS(1,2) EFFU(I)=MATSYS(1,3) ENDDO 99 END SUBROUTINE SYSMATRIX(THETA,REFANG,TPAR,TPER,RS,RP,DELPH, : ORIENT,CORANG,M,N,MATPOL,MATROP, : MATMIR,MATROT,WORK,MATSYS) C C This program sets the elements of the Mueller polarization C matrices C TPAR and TPER for MATPOL - polarizer matrix C RS, RP and DELPH for MATMIR - mirror retardance matrix C ORIENT and CORANG for MATROT - PA rotation matrix C and multiples them together to form the output matreix of C the poalrization response of the whole instrument MATSYS C IMPLICIT NONE INTEGER M ! No. of rows of matrix MATPOL INTEGER N ! No. of columns of matrix MATPOL REAL THETA ! Position angle of polarizer REAL REFANG ! Reference angle to subtract from THETA REAL TPAR ! The parallel component of polarizer transmission REAL TPER ! The perpendicular component of polarizer transmission REAL RS ! Reflectance of s wave REAL RP ! Reflectance of s wave REAL DELPH ! Retardance of s wave relative to p wave REAL ORIENT ! The position angle of the telescope/instrument REAL CORANG ! Additive angle correction REAL MATROP(M,N) ! Mueller matrix of polarizer rotation REAL MATPOL(M,N) ! Mueller matrix of polarizer REAL MATMIR(M,N) ! Mueller matrix of mirror polarization REAL MATROT(M,N) ! Mueller matrix of PA rotation REAL WORK(M,N) ! Work array for matrix multiplications REAL MATSYS(M,N) ! Output Mueller matrix of system polarization C C Form the polarizer matrix MATPOL C CALL MATRIX_POL(M,N,TPAR,TPER,MATPOL) CALL MATRIX_POLROT(M,N,THETA,REFANG,MATROP) C C Form the pick-off mirror retardance matrix C CALL MATRIX_MIRPOL(M,N,RS,RP,DELPH,MATMIR) C C Form the instrument rotation matrix C CALL MATRIX_PAROT(M,N,ORIENT,CORANG,MATROT) C C Multiply the three matrices together C CALL MATRIX_PROD4(M,N,MATPOL,MATROP,MATMIR,MATROT, : WORK,MATSYS) 99 END SUBROUTINE MATRIX_POL(M,N,TPAR,TPER,MATPOL) C C This subroutine sets the terms of the polarizer Mueller C matrix MATPOL using the values of the polarizer transmission C parallel (TPAR) and perpendicular (TPER) to the axis of C the polarizer C IMPLICIT NONE INTEGER M ! No. of rows of matrix MATPOL INTEGER N ! No. of columns of matrix MATPOL REAL TPAR ! The parallel component of polarizer transmission REAL TPER ! The perpendicular component of polarizer transmission REAL MATPOL(M,N) ! Output Mueller matrix of polarizer C C Local variables C INTEGER I,J C C Initialize output array C DO I=1,M,1 DO J=1,N,1 MATPOL(I,J)=0.0 ENDDO ENDDO C C Set the terms C MATPOL(1,1)=0.50*(TPAR+TPER) MATPOL(1,2)=0.50*(TPAR-TPER) MATPOL(2,1)=0.50*(TPAR-TPER) MATPOL(2,2)=0.50*(TPAR+TPER) MATPOL(3,3)=SQRT(TPAR*TPER) MATPOL(4,4)=SQRT(TPAR*TPER) 99 END SUBROUTINE MATRIX_POLROT(M,N,THETA,REFANG,MATROP) C C This subroutine sets the terms of the polarizer rotation C Mueller matrix MATROP using the values of the polarizer C rotation angle THETA and the reference direction (to C correct for instrument rotation) C IMPLICIT NONE INTEGER M ! No. of rows of matrix MATRPO INTEGER N ! No. of columns of matrix MATROP REAL THETA ! The position angle of the polarizer REAL REFANG ! Reference angle to subtract from THETA REAL MATROP(M,N) ! Output Mueller matrix of polarizer C C Local variables C INTEGER I,J REAL SINRD REAL COSRD C C Initialize output array C DO I=1,M,1 DO J=1,N,1 MATROP(I,J)=0.0 ENDDO ENDDO C C Set the terms C MATROP(1,1)=1.00 MATROP(2,2)=COSRD(2.*(THETA-REFANG)) MATROP(2,3)=SINRD(2.*(THETA-REFANG)) MATROP(3,2)=-SINRD(2.*(THETA-REFANG)) MATROP(3,3)=COSRD(2.*(THETA-REFANG)) MATROP(4,4)=1.00 99 END SUBROUTINE MATRIX_MIRPOL(M,N,RS,RP,DELPH,MATMIR) C C This subroutine sets the terms of the mirror Mueller C matrix MATMIR which describes the polarization induced C by mirro reflection (WFPC2 pick-off mirror in particular). C RS is the transverse electric (s) wave reflectance and C RP the transverse magnetic (p) wave reflectance C and DELPH (degrees) is the circular retardance of the s C wave relative to the p wave C IMPLICIT NONE INTEGER M ! No. of rows of matrix MATPOL INTEGER N ! No. of columns of matrix MATPOL REAL RS ! Reflectance of s wave REAL RP ! Reflectance of s wave REAL DELPH ! Retardance of s wave relative to p wave REAL MATMIR(M,N) ! Output Mueller matrix of mirror polarization C C Local variables C INTEGER I,J REAL A,B,C,D REAL SINRD REAL COSRD C C Initialize output array C DO I=1,M,1 DO J=1,N,1 MATMIR(I,J)=0.0 ENDDO ENDDO A=0.50*(RS+RP) B=0.50*(RS-RP) C=SQRT(RS*RP)*COSRD(DELPH) D=SQRT(RS*RP)*SINRD(DELPH) C C Set the terms C MATMIR(1,1)=A MATMIR(1,2)=B MATMIR(2,1)=B MATMIR(2,2)=A MATMIR(3,3)=C MATMIR(3,4)=-D MATMIR(4,3)=D MATMIR(4,4)=C 99 END SUBROUTINE MATRIX_PAROT(M,N,PA_V3,CORANG,MATROT) C C This subroutine sets the terms of the rotation Mueller C matrix MATROT to describe the rotation of the C telescope/instrument (PA_V3). PA_V3 is the telescope C position angle and CORANG is an additive correction C to give the orientation in the correct frame C IMPLICIT NONE INTEGER M ! No. of rows of matrix MATROT INTEGER N ! No. of columns of matrix MATROT REAL PA_V3 ! The position angle of the telescope/instrument REAL CORANG ! Additive angle correction REAL MATROT(M,N) ! Output Mueller matrix of polarizer C C Local variables C INTEGER I,J REAL SINRD REAL COSRD C C Initialize output array C DO I=1,M,1 DO J=1,N,1 MATROT(I,J)=0.0 ENDDO ENDDO C C Set the terms C MATROT(1,1)=1.00 MATROT(2,2)=COSRD(2.*(PA_V3+CORANG)) MATROT(2,3)=SINRD(2.*(PA_V3+CORANG)) MATROT(3,2)=-SINRD(2.*(PA_V3+CORANG)) MATROT(3,3)=COSRD(2.*(PA_V3+CORANG)) MATROT(4,4)=1.00 99 END SUBROUTINE MATRIX_PROD4(M,N,MATIN1,MATIN2,MATIN3,MATIN4, : WORK,MATOUT) C C This subroutine succesively multiples the four input C MxN matrices MATIN1, MATIN2, MATIN3 and MATIN4 together C and forms the output matrix MATOUT as the product C IMPLICIT NONE INTEGER M ! No. of rows of matrices INTEGER N ! No. of columns of matrices REAL MATIN1(M,N) ! First input matrix for product REAL MATIN2(M,N) ! Second input matrix for product REAL MATIN3(M,N) ! Third input matrix for product REAL MATIN4(M,N) ! Forth input matrix for product REAL WORK(M,N) ! Work array for intermediate results REAL MATOUT(M,N) ! Output product matrix C C Local variables C INTEGER I,J C C Initialize output and work arrays C DO I=1,M,1 DO J=1,N,1 WORK(I,J)=0.0 MATOUT(I,J)=0.0 ENDDO ENDDO C C Multiply the matrices MATIN1 AND MATIN2 together C CALL MATRIX_MULT(M,N,MATIN1,MATIN2,MATOUT) C C Multiply the matrix MATIN3 by the running product C CALL MATRIX_MULT(M,N,MATOUT,MATIN3,WORK) C C Multiply the matrix MATIN4 by the running product C CALL MATRIX_MULT(M,N,WORK,MATIN4,MATOUT) 99 END SUBROUTINE MATRIX_MULT(M,N,MATIN1,MATIN2,MATOUT) C C This subroutine multiples the square matrix MATIN1 C by MATIN2 writing the result to MATOUT C IMPLICIT NONE INTEGER M ! No. of rows of matrices INTEGER N ! No. of columns of matrices REAL MATIN1(M,N) ! First input matrix for product REAL MATIN2(M,N) ! Second input matrix for product REAL MATOUT(M,N) ! Output product matrix C C Local variables C INTEGER I,J,K REAL SUMP C C Initialize output array C DO I=1,M,1 DO J=1,N,1 MATOUT(I,J)=0.0 ENDDO ENDDO C C Compute the terms of the product matrix MATIN1*MATIN2 C DO I=1,M,1 DO J=1,N,1 SUMP=0.0 DO K=1,M,1 SUMP=SUMP+MATIN1(I,K)*MATIN2(K,J) ENDDO MATOUT(I,J)=SUMP ENDDO ENDDO 99 END SUBROUTINE EFFCHK(N,EFF1,EFF2,EFF3,LCHK) C C This subroutine returns LCHK true if and only if the C N values of EFF1 are equal AND the N values of EFF2 C are equal AND the N values of EFF3 are equal C IMPLICIT NONE INTEGER N REAL EFF1(N) REAL EFF2(N) REAL EFF3(N) LOGICAL LCHK C C Local variables C INTEGER I LOGICAL LCHK1 LOGICAL LCHK2 LOGICAL LCHK3 C C Initialize logical variables C LCHK=.FALSE. LCHK1=.TRUE. LCHK2=.TRUE. LCHK3=.TRUE. DO I=1,N-1,1 IF (EFF1(I).NE.EFF1(I+1)) THEN LCHK1=.FALSE. ENDIF IF (EFF2(I).NE.EFF2(I+1)) THEN LCHK2=.FALSE. ENDIF IF (EFF3(I).NE.EFF3(I+1)) THEN LCHK3=.FALSE. ENDIF ENDDO IF (LCHK1.AND.LCHK2.AND.LCHK3) THEN LCHK=.TRUE. ENDIF 99 END SUBROUTINE SORT2(N,ARR,BRR,ISTACK) C C This is the Numerical Recipes Quicksort algorithm to C sort two input arrays in the ascending numerical order C of the first (ARR) C IMPLICIT NONE INTEGER N ! Number of points in arrays to be sorted INTEGER ISTACK(N) ! Work array for sorting index REAL ARR(N) ! First input array to be sorted REAL BRR(N) ! Second input array to be sorted C C Local variables C INTEGER M,NSTACK,STAT PARAMETER (M=7) ! M is size of sorted subarrays INTEGER I,IR,J,JSTACK,K,L REAL A,B,TEMP NSTACK=N JSTACK=0 L=1 IR=N C C Perform insertion sort when subarray small enough C 1 IF (IR-L.LT.M) THEN DO J=L+1,IR,1 A=ARR(J) B=BRR(J) DO I=J-1,L,-1 IF (ARR(I).LE.A) THEN GO TO 2 ENDIF ARR(I+1)=ARR(I) BRR(I+1)=BRR(I) ENDDO I=L-1 2 ARR(I+1)=A BRR(I+1)=B ENDDO C C Pop stack and begin a new round of partitioning C IF (JSTACK.EQ.0) THEN GO TO 99 ENDIF IR=ISTACK(JSTACK) L=ISTACK(JSTACK-1) JSTACK=JSTACK-2 ELSE C C Choose median of left, centre and right elements as C partitioning element A. Also rearrange so that A(L+1) C A C 3 I=I+1 IF (ARR(I).LT.A) THEN GO TO 3 ENDIF 4 J=J-1 IF (ARR(J).GT.A) THEN GO TO 4 ENDIF IF (J.LT.I) THEN GO TO 5 ENDIF C C Exchange elements of both arrays C TEMP=ARR(I) ARR(I)=ARR(J) ARR(J)=TEMP TEMP=BRR(I) BRR(I)=BRR(J) BRR(J)=TEMP GO TO 3 C C End of innermost loop C Insert partitioning element in both arrays C 5 ARR(L)=ARR(J) ARR(J)=A BRR(L)=BRR(J) BRR(J)=B JSTACK=JSTACK+2 C C Push pointers to larher subarray on stack, process C smaller subarray C IF (JSTACK.GT.NSTACK) THEN CALL UMSPUT(' Error in sorting',3,0,STAT) GO TO 99 ENDIF IF (IR-I+1.GE.J-l) THEN ISTACK(JSTACK)=IR ISTACK(JSTACK-1)=I IR=J-1 ELSE ISTACK(JSTACK)=J-1 ISTACK(JSTACK-1)=L L=I ENDIF ENDIF GO TO 1 99 END SUBROUTINE POLCALF(N,INTV,ERRV,ANANG,EFFI,EFFQ,EFFU, : INSTPO,INSTPA,LPLIM,ARRI,COVI,INDX, : HEFFI,HEFFQ,HEFFU,DAQ,X,Y,SIGY,YFIT, : WT,STI,STIER,STQ,STQER,STU,STUER, : LPOL,LPOLER,LPA,LPAER,IFAIL) C C This subroutine calculates the linear polarization (%) C and error and the position angle (degrees) and error and C the Stokes I parameter, the normalised Q and U Stokes C parameters for the N images taken at polarizer position C angles given by ANANG from a fit to: C I=INTI(EFFI + EFFQ*Q + EFFU*U) C If N only equals 3, then the set of linear equations are C solved as a matrix equation by LU decomposition (SOLVIQU3). C If N is greater than 3 then the set of equations is C solved by least squares using CURFIT (called by SOLVIQU4). C C If the polarization error is greater than LPLIM then the C polarization and position angle are not copied to the C output images. The errors are calculated from the input C errors ERRV. The instrumental polarization is corrected. C A pixel with a polarization error of 0 signifies either C no data or not sufficient signal-to-noise to compute C polarization C IMPLICIT NONE INTEGER N ! Number of input values for polarization determination INTEGER INDX(3) ! Vector of row permution of array ARRI in LU decomposition INTEGER DAQ(N) ! Array of data quality of input values INTEGER IFAIL ! Error flag for least squares estimation; non zero value implies error REAL INTV(N) ! Input array of observed signal values with ANANG REAL ERRV(N) ! Input array of error values with ANANG REAL ANANG(N) ! Input array of polarizer angles of signal values REAL EFFI(N) ! Efficiency for Stokes I detection REAL EFFQ(N) ! Efficiency for Stokes Q detection REAL EFFU(N) ! Efficiency for Stokes U detection REAL HEFFI(N) ! Holding array for efficiency for Stokes I detection REAL HEFFQ(N) ! Holding array for efficiency for Stokes Q detection REAL HEFFU(N) ! Holding array for efficiency for Stokes U detection REAL INSTPO ! Instrumental (fractional) polarization REAL INSTPA ! Position angle of instrumental polarization REAL LPLIM ! Limiting polarization error to allow ouput poalrization REAL STI ! Output total intensity (Stokes I) REAL STIER ! Output total intensity error REAL STQ ! Output normalised Stokes Q parameter REAL STQER ! Output normalised Stokes Q parameter error REAL STU ! Output normalised Stokes U parameter REAL STUER ! Output normalised Stokes U parameter error REAL LPOL ! Ouput fractional linear polarization REAL LPOLER ! Ouput fractional linear polarization error REAL LPA ! Outut polarization position angle REAL LPAER ! Outut polarization position angle error DOUBLE PRECISION ARRI(3,3) ! Array of LU decomposition of system efficiences array DOUBLE PRECISION COVI(3,3) ! Array of covariance matrix (inverse of system efficiences array) DOUBLE PRECISION X(N) ! Array if X values (ANG) for fitting DOUBLE PRECISION Y(N) ! Array of Y values (INTV) for fitting DOUBLE PRECISION SIGY(N) ! Array of error values (ERRV) for fitting DOUBLE PRECISION YFIT(N) ! Array of fitted Y values in CURFIT DOUBLE PRECISION WT(N) ! Array of weights for Y values in CURFIT C C Local variables C INTEGER NP INTEGER MODE REAL RAD2 REAL MEANI REAL MEANIER REAL QI REAL QIER REAL UI REAL UIER REAL Q REAL QER REAL U REAL UER REAL POL REAL POLB REAL ERPOL REAL TH REAL ERTH LOGICAL LNOZERO RAD2=180.0/3.1415926535898 C C Normalise the non-zero signal-sky values by the mean C and copy the normalised values, propagated errors, C data quality and weights to the arrays Y, SIGY, C DQ and WT for fitting. Copy the polarizer angle C values to array X. The system efficiencies are C copied to holding arrays to cover the cases when C the number of returned points is less than the C number of input values. Return an estimate of the C mean and its error C MODE=1 ! Weight points by 1/err**2 50 CALL POLNORMF(N,INTV,ERRV,ANANG,MODE,EFFI,EFFQ, : EFFU,NP,X,Y,SIGY,DAQ,WT,HEFFI,HEFFQ, : HEFFU,MEANI,MEANIER,LNOZERO) C C If number of non-zero points is 3 but input number was C greater, need to form the matrix of system efficiences, C the LU decomposition of this matrix and the inverse of C the matrix of system efficiences for the covariance matrix C IF (NP.EQ.3.AND.N.GT.3) THEN CALL LUSETUP(NP,HEFFI,HEFFQ,HEFFU,ARRI,COVI,INDX) ENDIF C C If less than three points then not possible to perform C fit. Zero the output values and exit C IF (NP.LT.3) THEN MEANI=0.0 MEANIER=0.0 Q=0.0 QER=0.0 U=0.0 UER=0.0 POL=0.0 ERPOL=0.0 TH=0.0 ERTH=0.0 GO TO 80 ENDIF C C If three points only then need to solve for I, Q and U C as a set of linear equations using the LU decomposition C (ARRI) and the covariance matrix (COVI) C IF (NP.EQ.3) THEN CALL SOLVIQU3(NP,Y,SIGY,ARRI,COVI,INDX,MEANI, : MEANIER,QI,QIER,UI,UIER) ENDIF C C If more than three points then solve the system using C least squares with the subroutine CURFIT C IF (NP.GT.3) THEN CALL SOLVIQU4(NP,X,Y,SIGY,WT,DAQ,HEFFI,HEFFQ, : HEFFU,YFIT,MEANI,MEANIER,QI, : QIER,UI,UIER,IFAIL) IF (IFAIL.NE.0) THEN MEANI=0.0 MEANIER=0.0 Q=0.0 QER=0.0 U=0.0 UER=0.0 POL=0.0 ERPOL=0.0 TH=0.0 ERTH=0.0 GO TO 80 ENDIF ENDIF C C Set to zero the ouput Q and U errors if no input errors C available C IF (.NOT.LNOZERO) THEN QIER=0.0 UIER=0.0 ENDIF C C Check that MEANI is not zero C IF (MEANI.LE.0.0) THEN MEANIER=0.0 Q=0.0 QER=0.0 U=0.0 UER=0.0 POL=0.0 ERPOL=0.0 TH=0.0 ERTH=0.0 GO TO 80 ENDIF C C Correct Stokes Q and U and errors for the instrumental C polarization. If instrumental polarization zero then C don't attempt correction. Deal with case when errors zero C IF (INSTPO.NE.0.0) THEN CALL INSPOL(QI,UI,QIER,UIER,INSTPO,INSTPA,Q,U,QER,UER) ELSE Q=QI QER=QIER U=UI UER=UIER ENDIF IF (.NOT.LNOZERO) THEN QER=0.0 UER=0.0 ENDIF C C Compute polarization and error corrected for instrumental C polarization C ERPOL=SQRT(( ((Q*QER)**2.) + ((U*UER)**2.) )/ : (Q**2. + U**2.)) POLB=SQRT(Q*Q + U*U) C C Correct polarization for bias C IF (ERPOL.NE.0.0) THEN CALL DEBIAS(POLB,ERPOL,POL) ELSE POL=POLB ENDIF C C If polarization error greater than limit, set polarization C and position angle and errors to zero C IF (ERPOL.GT.LPLIM.AND.ERPOL.NE.0.0) THEN MEANI=0.0 MEANIER=0.0 Q=0.0 QER=0.0 U=0.0 UER=0.0 POL=0.0 ERPOL=0.0 TH=0.0 ERTH=0.0 GO TO 80 ENDIF C C Compute position angle corrected for instrumental C polarization C TH=0.5*RAD2*ATAN2(U,Q) C C Distinguish the possible cases: C P zero and SIGMA(P) zero C P non-zero and SIGMA(P) zero C P zero and SIGMA(P) non-zero C P non-zero and SIGMA(P) non-zero C If P/SIGMA(P) < 8 the distribution function of C sigma(theta) is used to calculate the 1sigma_theta(v) C value using Naghizadeh-Khouei & Clarke, 1993, equ. 5, C otherwise the analystic expression is used C IF (POL.EQ.0.0.AND.ERPOL.EQ.0.0) THEN ERTH=0.0 GO TO 80 ENDIF IF (POL.NE.0.0.AND.ERPOL.EQ.0.0) THEN ERTH=0.0 GO TO 80 ENDIF IF (POL.EQ.0.0.AND.ERPOL.NE.0.0) THEN ERTH=180.0/SQRT(12.0) ENDIF IF (POL.NE.0.0.AND.ERPOL.NE.0.0) THEN IF ((POL/ERPOL).LT.8.0) THEN CALL PPASIG(POL,ERPOL,TH,ERTH) ELSE ERTH=0.50*RAD2*SQRT((Q*UER)**2. + (U*QER)**2.)/ : (Q**2. + U**2.) ENDIF GO TO 80 ENDIF 80 STI=MEANI STIER=ABS(MEANIER) STQ=Q STQER=ABS(QER) STU=U STUER=ABS(UER) LPOL=POL*100.0 LPOLER=ABS(ERPOL)*100.0 LPA=TH LPAER=ABS(ERTH) 99 END SUBROUTINE POLNORMF(N,OBJV,ERRV,ANG,MODE,EFFI,EFFQ,EFFU, : NP,X,Y,DY,DQ,WT,HEFFI,HEFFQ,HEFFU, : SMEAN,SMEANER,LNOZERO) C C This subroutine takes the input object and error values and C outputs the values and their errors. If a value has both C object and error values zero, then the output array X and C Y values are not filled. The data quality (DQ) and weight C arrays (WT) are also set and the number of valid output X,Y C pairs. The values of the instrumental efficiencies EFFI, C EFFQ and EFFU corresponding to the non-zero values are C written to the output arrays HEFFI, HEFFQ and HEFFU. C If both HEFFQ and HEFFU are zero then the flag LEFF is C set true. The mean Y value is also returned as SMEAN C and its error as SMEANER. If the error values are all C zero, LNOZERO is returned false C IMPLICIT NONE INTEGER N ! No. of points in input arrays INTEGER MODE ! Method of weighting data points: C C Weights: MODE=-1 -> weight=1/Y; MODE=0 -> weight=0; MODE=1 -> weight=1/error**2 C INTEGER NP ! No. of points in output arrays INTEGER DQ(N) ! Output array of data quality of points REAL OBJV(N) ! Input array of signal values REAL ERRV(N) ! Input array signal error values REAL ANG(N) ! Input polarizer angles correspomdomng to signal values REAL EFFI(N) ! Efficiency for Stokes I detection REAL EFFQ(N) ! Efficiency for Stokes Q detection REAL EFFU(N) ! Efficiency for Stokes U detection REAL HEFFI(N) ! Output efficiency for Stokes I detection REAL HEFFQ(N) ! Output efficiency for Stokes Q detection REAL HEFFU(N) ! Output efficiency for Stokes U detection REAL SMEAN ! Output mean signal value REAL SMEANER ! Output rms on the mean signal value DOUBLE PRECISION X(N) ! Output array of X (polarizer angle) DOUBLE PRECISION Y(N) ! Output array of Y (normalized signal) DOUBLE PRECISION DY(N) ! Output array of Y error DOUBLE PRECISION WT(N) ! Array of weights for Y values in CURFIT LOGICAL LNOZERO ! Flag to indicate output error array all zeroes (=FALSE) C C Local variables C INTEGER I INTEGER J INTEGER ISUM REAL SUM C C Form the mean of the non-zero values. Only include C values for which EFFQ and EFFU are non-zero (have C to be non-zero to solve for polarization) C ISUM=0 SUM=0.0 DO I=1,N,1 IF ((OBJV(I).GT.1.0E-6). : AND.(.NOT.(EFFQ(I).EQ.0.0.AND.EFFU(I).EQ.0.0))) THEN SUM=SUM + OBJV(I) ISUM=ISUM+1 ENDIF ENDDO C C Form the mean Y value C IF (ISUM.NE.0) THEN SMEAN=SUM/REAL(ISUM) ELSE NP=0 GO TO 99 ENDIF IF (SMEAN.LE.0.0) THEN NP=0 GO TO 99 ENDIF C C Compute the rms on the mean C ISUM=0 SUM=0.0 DO I=1,N,1 IF ((OBJV(I).GT.1.0E-6). : AND.(.NOT.(EFFQ(I).EQ.0.0.AND.EFFU(I).EQ.0.0))) THEN SUM=SUM + (OBJV(I)-SMEAN)**2 ISUM=ISUM+1 ENDIF ENDDO C C Form the RMS on the mean Y value C IF (ISUM.NE.0) THEN IF (SUM.GT.0.0) THEN SMEANER=SQRT(SUM/REAL(ISUM)) ELSE SMEANER=0.0 ENDIF ELSE SMEANER=0.0 ENDIF C C For the non-zero input signal values set the data qaulity C to good, form the weight and normalise the object-sky C values by the mean and propagate the errors C J=1 DO I=1,N,1 IF ((OBJV(I).GT.0.0). : AND.(.NOT.(EFFQ(I).EQ.0.0.AND.EFFU(I).EQ.0.0))) THEN C C Set the data quality to good (=0) C DQ(J)=0 C C Evaluate weights neglecting data quality C IF (MODE.EQ.-1) THEN C C Weight by 1/value C IF (OBJV(I).EQ.0.0) THEN WT(J)=1.0D0 ELSE WT(J)=1.0D0/DBLE(OBJV(I)) ENDIF ENDIF IF (MODE.EQ.0) THEN C C No weighting C WT(J)=1.0D0 ENDIF IF (MODE.EQ.1) THEN C C Weight by 1/sigmay**2 C IF (ERRV(I).EQ.0.0) THEN WT(J)=1.0D0 ELSE WT(J)=1.0D0/DBLE(ERRV(I)*ERRV(I)) ENDIF ENDIF C C Copy the angle, value and error to output arrays C X(J)=DBLE(ANG(I)) Y(J)=DBLE(OBJV(I)) DY(J)=DBLE(ERRV(I)) C C Set the system efficiencies to those appropriate C to the non-zero ouput points C HEFFI(J)=EFFI(I) HEFFQ(J)=EFFQ(I) HEFFU(J)=EFFU(I) c IF (HEFFQ(J).EQ.0.0.AND.HEFFU(J).EQ.0.0) THEN c LEFF=.TRUE. c ENDIF J=J+1 c ELSE C C Set the data quality to not good (not =0) C c DQ(J)=1 ENDIF ENDDO NP=J-1 C C Check that all errors not zero. C LNOZERO=.FALSE. DO I=1,NP,1 IF (DY(I).NE.0.0D0) THEN LNOZERO=.TRUE. ENDIF ENDDO 99 END SUBROUTINE LUSETUP(N,EFFI,EFFQ,EFFU,ARRI,COVI,INDX) C C This subroutine sets up the LU decomposition of the C array of system efficiences from EFFI, EFFQ and EFFU C and forms the inverse matrix COVI of the matrix of C system efficiences for the covariance matrix C IMPLICIT NONE INTEGER N ! No of points in input arrays INTEGER INDX(N) ! Output vector of row permution in LU decomposition REAL EFFI(N) ! Array of efficiences for Stokes I REAL EFFQ(N) ! Array of efficiences for Stokes Q REAL EFFU(N) ! Array of efficiences for Stokes U DOUBLE PRECISION ARRI(N,N) ! Output array of LU decomposition of system efficiences array DOUBLE PRECISION COVI(N,N) ! Output array of covariance matrix (inverse of system efficiences array) C C Local variables C INTEGER I,J DOUBLE PRECISION DS C C Copy the efficiency arrays to a single 3x3 array ARR C DO I=1,N,1 ARRI(I,1)=DBLE(EFFI(I)) ENDDO DO I=1,N,1 ARRI(I,2)=DBLE(EFFQ(I)) ENDDO DO I=1,N,1 ARRI(I,3)=DBLE(EFFU(I)) ENDDO C C Form the LU decomposition of ARR C CALL LUDCMP(N,N,ARRI,INDX,DS) C C Initialize the covariance matrix COVI, which is the C inverse of ARRI C DO I=1,N,1 DO J=1,N,1 COVI(I,J)=0.0D0 ENDDO COVI(I,I)=1.0D0 ENDDO C C Invert the matrix ARRI to COVI C DO J=1,N,1 CALL LUBKSB(N,N,ARRI,INDX,COVI(1,J)) ENDDO 99 END SUBROUTINE SOLVIQU3(N,Y,SIGY,ARRI,COVI,INDX,MEANI,MEANIER, : QI,QIER,UI,UIER) C C This subroutine calculates the Stokes I, Q and U C (MEANI, QI and UI) and their errors from the C LU decomposition of the system of three linear C equations: C Si = effi(i)*I + effq(i)*Q + effu(i)*U C given by the array ARRI. The errors on Stokes I, C Q and U (MEANIER, QIER and UIER) are calculated from C the covariance matrix COVI which is the inverse of the C matrix of the array of system efficiences C IMPLICIT NONE INTEGER N ! No of points in input arrays INTEGER INDX(N) ! ! Vector of row permution applied for LU decomposition REAL MEANI ! Output Stokes total intensity REAL MEANIER ! Output error on Stokes I REAL QI ! Output normalised Stokes Q REAL QIER ! Output error on Stokes Q REAL UI ! Output normalised Stokes U REAL UIER ! Output error on Stokes U DOUBLE PRECISION Y(N) ! Array of input observed signal values DOUBLE PRECISION SIGY(N) ! Array of input error on observed signal DOUBLE PRECISION ARRI(N,N) ! Array of LU decomposition of system efficiences array DOUBLE PRECISION COVI(N,N) ! Array of covariance matrix (inverse of system efficiences array) C C Local variables C INTEGER I DOUBLE PRECISION AA(3) DOUBLE PRECISION DD(3) LOGICAL LNZERO C C Copy the signal values to AA in preparation for solution of C the system of linear equations. Trap the case where all values C are zero C LNZERO=.FALSE. DO I=1,N,1 AA(I)=Y(I) IF (AA(I).NE.0.0D0) THEN LNZERO=.TRUE. ENDIF ENDDO IF (.NOT.LNZERO) THEN MEANI=0.0 MEANIER=0.0 QI=0.0 QIER=0.0 UI=0.0 UIER=0.0 GO TO 99 ENDIF C C Solve the system of linear equations from the LU decomposition C CALL LUBKSB(N,N,ARRI,INDX,AA) C C Copy the solution to the output variables for Stokes C I and normalised Q and U C IF (AA(1).NE.0.0D0) THEN MEANI=REAL(AA(1)) QI=REAL(AA(2)/AA(1)) UI=REAL(AA(3)/AA(1)) ELSE MEANI=0.0 MEANIER=0.0 QI=0.0 QIER=0.0 UI=0.0 UIER=0.0 GO TO 99 ENDIF C C Call the subroutine COVAR3 to estimate the errors on Stokes C I, Q and U from the array of the covariance matrix C CALL COVAR3(N,SIGY,COVI,DD) C C Form the output Stokes I and normalised Stokes Q and U C errors C MEANIER=REAL(DD(1)) IF (AA(2).NE.0.0D0.AND.AA(1).NE.0.0D0) THEN QIER=QI*REAL(DSQRT( (DD(2)/AA(2))**2. + (DD(1)/AA(1))**2. )) ELSE IF (AA(1).NE.0.0D0) THEN QIER=QI*REAL(DSQRT((DD(1)/AA(1))**2.)) ELSE QIER=0.0D0 ENDIF ENDIF IF (AA(3).NE.0.0D0) THEN UIER=UI*REAL(DSQRT( (DD(3)/AA(3))**2. + (DD(1)/AA(1))**2. )) ELSE IF (AA(1).NE.0.0D0) THEN UIER=UI*REAL(DSQRT((DD(1)/AA(1))**2.)) ELSE UIER=0.0D0 ENDIF ENDIF 99 END SUBROUTINE LUDCMP(NP,N,A,INDX,D) C C This is the Numerical Recipes subroutine to preform LU C of a matrix. C Given a matrix A(1:N,1:N), with physical dimension NP by NP, this C routine replaces it by the LU decomposition of a rowwise permutation C of itself. A and N are input. A is output, arranged by columns from C left to right and within each column from top to bottom (Crout's C method). INDX(1:N) is an output vector that records the row C permutation affected by the partial pivoting; D is output as + C or - depending on whether the number of row interchnages was even or C odd respectively. This routine is used in combination with C LUBKSB to solve a suystem of linear equations or invert a matrix C IMPLICIT NONE INTEGER N ! Square region of matrix to apply decomposition (1:N) INTEGER NP ! Input dimension of square matrix INTEGER INDX(N) ! Output vector of row permution applied DOUBLE PRECISION A(NP,NP) ! Input array to be decomposed DOUBLE PRECISION D ! Indicator of sign for oven/odd row interchanges C C Local variables C INTEGER I,J,K INTEGER IMAX INTEGER NMAX INTEGER STAT PARAMETER (NMAX=3) DOUBLE PRECISION AAMAX DOUBLE PRECISION DUM DOUBLE PRECISION SUM DOUBLE PRECISION VV(NMAX) DOUBLE PRECISION TINY PARAMETER (TINY=1.0D-25) CHARACTER*80 OUTEXT C C Loop over the rows to get the implicit scaling information C D=1.0D0 DO I=1,N,1 AAMAX=0.0D0 DO J=1,N,1 IF (ABS(A(I,J)).GT.AAMAX) THEN AAMAX=ABS(A(I,J)) ENDIF ENDDO C C Exit if no non-zero largest element C IF (AAMAX.EQ.0.0D0) THEN WRITE(OUTEXT,15) 15 FORMAT(' Error in Matrix decomposition - singlular matrix') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 99 ENDIF VV(I)=1.0D0/AAMAX ENDDO C C Loop over the columsn for Crout's method C DO J=1,N,1 DO I=1,J-1,1 SUM=A(I,J) DO K=1,I-1,1 SUM=SUM-A(I,K)*A(K,J) ENDDO A(I,J)=SUM ENDDO AAMAX=0.0D0 DO I=J,N,1 SUM=A(I,J) DO K=1,J-1,1 SUM=SUM - A(I,K)*A(K,J) ENDDO A(I,J)=SUM DUM=VV(I)*ABS(SUM) IF (DUM.GE.AAMAX) THEN IMAX=I AAMAX=DUM ENDIF ENDDO C C Check if rows need to be interchanged C IF (J.NE.IMAX) THEN DO K=1,N,1 DUM=A(IMAX,K) A(IMAX,K)=A(J,K) A(J,K)=DUM ENDDO D=-D VV(IMAX)=VV(J) ENDIF INDX(J)=IMAX C C If pivot element is zero matrix is singlular to accuracy of C TINT. TINY could be made smaller C IF (A(J,J).EQ.0.0D0) THEN A(J,J)=TINY ENDIF C C Divide by the pivot element C IF (J.NE.N) THEN DUM=1.0D0/A(J,J) DO I=J+1,N,1 A(I,J)=A(I,J)*DUM ENDDO ENDIF ENDDO 99 END SUBROUTINE LUBKSB(NP,N,A,INDX,B) C C This is the Numerical Recipes subroutine to solve the set of C linear equation A * X = B, where A is the LU decomposition of C a matrix determined by LUDCMP. INDX is the permutation vector C returned by LUDCMP. B(1:N) is input as the right hand side C vector B and returns the solution vector X. A is not modified C by this subroutine. C IMPLICIT NONE INTEGER N ! Square region of matrix to apply decomposition (1:N) INTEGER NP ! Input dimension of square matrix INTEGER INDX(N) ! Input vector of row permution applied DOUBLE PRECISION A(NP,NP) ! Input array to be decomposed DOUBLE PRECISION B(N) ! Input right hand side vector; output solution vector C C Local variables C INTEGER I,J INTEGER II,LL DOUBLE PRECISION SUM C C Perform forward substitution, unscrambling the permutation C II=0 DO I=1,N,1 LL=INDX(I) SUM=B(LL) B(LL)=B(I) IF (II.NE.0) THEN DO J=II,I-1,1 SUM=SUM-A(I,J)*B(J) ENDDO ELSE IF (SUM.NE.0.0D0) THEN II=I ENDIF B(I)=SUM ENDDO C C Perform the backsubstitution C DO I=N,1,-1 SUM=B(I) DO J=I+1,N,1 SUM=SUM-A(I,J)*B(J) ENDDO B(I)=SUM/A(I,I) ENDDO 99 END SUBROUTINE COVAR3(N,SIGY,ARR,DD) C C This subroutine computes the errors estimates on C the Stokes parameters I, Q and U as the array DD C from the errors on the N signal values SIGY and C the array of polarizer system efficiences (Mueller C matrix) ARR C INTEGER N ! Input dimension of arrays DOUBLE PRECISION SIGY(N) ! Input array of sigmas on observed signal DOUBLE PRECISION ARR(N,N) ! Input array of Mueller matrix of system efficiencies DOUBLE PRECISION DD(N) ! Output array of Stokes I,Q,U propagated errors C C Local variables C INTEGER I C C Compute the errors on Stokes parameters from the covariance C matrix following Sparks and Axon, 1999, equation 5 C DO I=1,N,1 DD(I)=SQRT( ARR(I,1)*ARR(I,1)*SIGY(1)*SIGY(1) + : ARR(I,2)*ARR(I,2)*SIGY(2)*SIGY(2) + : ARR(I,3)*ARR(I,3)*SIGY(3)*SIGY(3) ) ENDDO 99 END SUBROUTINE SOLVIQU4(N,X,Y,SIGY,WT,DAQ,EFFI,EFFQ, : EFFU,YFIT,MEANI,MEANIER,QI, : QIER,UI,UIER,IFAIL) C C This subroutine passes the input values (X,Y), Y C error values, weight (WT) and data quality (DQ) C for the least squares determination of I, Q and U. C An estimate of the initial solution values are found C and the least squares fitted values of Stokes I (MEANI) C and error (MEANIER), normalised Stokes Q (QI) and C error (QIER) and Stokes U (UI) and error (UIER) are C returned C IMPLICIT NONE INTEGER N ! No. of inputs points to fit INTEGER DAQ(N) ! Array of data quality of input values INTEGER IFAIL ! Error flag; value not zero implies error REAL EFFI(N) ! Array of efficiences for Stokes I REAL EFFQ(N) ! Array of efficiences for Stokes Q REAL EFFU(N) ! Array of efficiences for Stokes U REAL MEANI ! Output Stokes total intensity REAL MEANIER ! Output error on Stokes I REAL QI ! Output normalised Stokes Q REAL QIER ! Output error on Stokes Q REAL UI ! Output normalised Stokes U REAL UIER ! Output error on Stokes U DOUBLE PRECISION X(N) ! Array if X values (ANG) for fitting DOUBLE PRECISION Y(N) ! Array of Y values (INTV) for fitting DOUBLE PRECISION SIGY(N) ! Array of error values (ERRV) for fitting DOUBLE PRECISION WT(N) ! Array of weights for Y values in CURFIT DOUBLE PRECISION YFIT(N) ! Array of fitted Y values in CURFIT C C Local variables C INTEGER I,J INTEGER NTERMS INTEGER BAD INTEGER NPC ! 1st dimension of common arrays CEFFI, CEFFQ and CEFFU PARAMETER (NPC=64) INTEGER NITER DOUBLE PRECISION CEFFI(NPC) DOUBLE PRECISION CEFFQ(NPC) DOUBLE PRECISION CEFFU(NPC) DOUBLE PRECISION A(3) DOUBLE PRECISION SIGA(3) DOUBLE PRECISION DELTAA(3) DOUBLE PRECISION FLAMDA DOUBLE PRECISION CHISQR DOUBLE PRECISION OCHISQ DOUBLE PRECISION DELCHI COMMON/LSQ/CEFFI,CEFFQ,CEFFU C C Fix parameters for fitting routine CURFIT C NTERMS=3 BAD=1 FLAMDA=1.0D-3 ! Initial value of proportion of gradient search C C Find the value of X (angle ANANG) corresponding to the C maximum absolute value of the Stokes parameters as a C first guess on the value of the value of the polarization C cosine 2*theta curve. From the maximum polarization estimate C initial values of Q (A(1)) and U (A(2)) C CALL MAXCURVE(N,X,Y,MEANI,A) C C Copy the arrays of system efficiences for Stokes I, C Q and U to holding arrays for COMMON C J=1 DO I=1,N,1 c IF (DAQ(I).EQ.0) THEN CEFFI(J)=DBLE(EFFI(I)) CEFFQ(J)=DBLE(EFFQ(I)) CEFFU(J)=DBLE(EFFU(I)) J=J+1 c ENDIF ENDDO C C Set up for multiple call of CURFIT until chi-squared C (CHISQ) is less than 1 or 100 iterations have been reached C NITER=100 OCHISQ=0.0D0 C C Call the non-linear least squares fitting function C CURFIT to determine the polarization and position C angle which fit the data points by a cosine curve C DO J=1,NITER,1 OCHISQ=CHISQR CALL CURFIT(N,X,Y,WT,DAQ,BAD,NTERMS,A,DELTAA, : SIGA,FLAMDA,YFIT,CHISQR,IFAIL) DELCHI=ABS(OCHISQ-CHISQR) IF (DELCHI.LT.1.0D-6) THEN c WRITE(OUTEXT,151) J,CHISQR,DELCHI c151 FORMAT(' NO. ITER, CHI-SQ, DELTA CHI-SQ ', c :I4,2X,F11.5,2X,E10.4) c CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 50 ENDIF c WRITE(OUTEXT,151) J,CHISQR,DELCHI c151 FORMAT(' NO. ITER, CHI-SQ, DELTA CHI-SQ ',I4,2X,F10.5,2X,E10.4) c CALL UMSPUT(OUTEXT,1,0,STAT) IF (IFAIL.NE.0) THEN MEANI=0.0 MEANIER=0.0 QI=0.0 QIER=0.0 UI=0.0 UIER=0.0 GO TO 99 ENDIF IF (CHISQR.LE.1.0D0) THEN GO TO 50 ENDIF ENDDO J=NITER C C Report number of iterations, final chi-squared and C last change in chi-squared per iteration C c50 IF (J.GT.1) THEN c WRITE(OUTEXT,51) J,CHISQR,DELCHI c51 FORMAT(' NO. ITER, CHI-SQ, DELTA CHI-SQ ',I4,2X,F10.5,2X,E10.4) c CALL UMSPUT(OUTEXT,1,0,STAT) c ENDIF C C Call CURFIT with FLAMDA=0 to get matrix of covariances C and errors on A C 50 FLAMDA=0.0D0 CALL CURFIT(N,X,Y,WT,DAQ,BAD,NTERMS,A,DELTAA, : SIGA,FLAMDA,YFIT,CHISQR,IFAIL) C C Form mean intensity (Stokes I) and error and normalised C Stokes Q and U and errors C IF (A(1).NE.0.0D0) THEN MEANI=REAL(A(1)) MEANIER=REAL(SIGA(1)) QI=REAL(A(2)/A(1)) IF (A(2).NE.0.0D0) THEN QIER=QI*REAL(DSQRT( (SIGA(2)/A(2))**2. + (SIGA(1)/A(1))**2. )) ELSE QIER=QI*REAL(DSQRT((SIGA(1)/A(1))**2.)) ENDIF UI=REAL(A(3)/A(1)) IF (A(3).NE.0.0D0) THEN UIER=UI*REAL(DSQRT( (SIGA(3)/A(3))**2. + (SIGA(1)/A(1))**2. )) ELSE UIER=UI*REAL(DSQRT((SIGA(1)/A(1))**2.)) ENDIF ELSE MEANI=0.0 MEANIER=0.0 QI=0.0 QIER=0.0 UI=0.0 UIER=0.0 ENDIF 99 END SUBROUTINE MAXCURVE(N,X,Y,MEANY,A) C C This subroutine determines a starting estimate of C the position angle X corresponding to the maximum C absolute value of the polarization cosine 2*theta C curve of the input array Y. C The polarization is estimated from the interpolated C maximum amplitude and used to provide initial values C of Stokes Q (A(1)) and Stokes U (A(2)) C IMPLICIT NONE INTEGER N ! No. of points in input arrays X and Y REAL MEANY ! Mean of Y values DOUBLE PRECISION X(N) ! Input array of position angles DOUBLE PRECISION Y(N) ! Input array of signal values DOUBLE PRECISION A(3) ! Output array of estimated Stokes Q and U parameters C C Local variables C INTEGER I INTEGER IPMIN1 INTEGER IPMIN2 INTEGER IPMAX1 INTEGER IPMAX2 DOUBLE PRECISION PMIN1 DOUBLE PRECISION PMIN2 DOUBLE PRECISION PMAX1 DOUBLE PRECISION PMAX2 DOUBLE PRECISION XMAX DOUBLE PRECISION XMIN DOUBLE PRECISION SINDD DOUBLE PRECISION COSDD LOGICAL LPMAX C C Determine maximum value of Y C PMIN1=1.0D32 PMAX1=-1.0D32 DO I=1,N,1 IF (Y(I).GE.PMAX1) THEN PMAX1=Y(I) IPMAX1=I ENDIF IF (Y(I).LE.PMIN1) THEN PMIN1=Y(I) IPMIN1=I ENDIF ENDDO C C Determine if maximum or minimum is C maximum absolute value C IF (ABS(PMAX1).GT.ABS(PMIN1)) THEN LPMAX=.TRUE. ELSE LPMAX=.FALSE. ENDIF C C Depending on whether maximum or minimum has C greatest absolute value, find the 2nd highest C absolute value C PMIN2=1.0D32 PMAX2=-1.0D32 IF (LPMAX) THEN DO I=1,N,1 IF (Y(I).GE.PMAX2.AND.I.NE.IPMAX1) THEN PMAX2=Y(I) IPMAX2=I ENDIF ENDDO ENDIF IF (.NOT.LPMAX) THEN DO I=1,N,1 IF (Y(I).LE.PMIN2.AND.I.NE.IPMIN1) THEN PMIN2=Y(I) IPMIN2=I ENDIF ENDDO ENDIF C C Linearly interpolate between two maximum absolute C values for estimated position angle of peak of C polarization cosine 2*theta curve C IF (LPMAX) THEN IF ((PMAX1 + PMAX2).NE.0.0D0) THEN XMAX=(PMAX1*X(IPMAX1) + PMAX2*X(IPMAX2))/(PMAX1 + PMAX2) A(2)=ABS(PMAX1-DBLE(MEANY))*COSDD(2.*XMAX) A(3)=ABS(PMAX1-DBLE(MEANY))*SINDD(2.*XMAX) ELSE A(2)=0.0D0 A(3)=0.0D0 ENDIF ENDIF IF (.NOT.LPMAX) THEN IF ((PMIN1 + PMIN2).NE.0.0D0) THEN XMIN=(PMIN1*X(IPMIN1) + PMIN2*X(IPMIN2))/(PMIN1 + PMIN2) A(2)=-ABS(PMIN1-DBLE(MEANY))*COSDD(2.*XMIN) A(3)=-ABS(PMIN1-DBLE(MEANY))*SINDD(2.*XMIN) ELSE A(2)=0.0D0 A(3)=0.0D0 ENDIF ENDIF 99 END SUBROUTINE CURFIT(NPTS,X,Y,WEIGHT,DAQ,BAD,NTERMS, : A,DELTAA,SIGMAA,FLAMDA,YFIT,CHISQR, : IFAIL) C C Performs a least squares fit to a non-linear function with C a linearization of the fitting function. Points with 'bad' C data quality are excluded from the fit. C The program is taken from Bevington, 1968, p.237-9. C C Subprograms called: C FUNCTN(NPTS,X,I,N,A) - evaluates the fitting function for the C Ith term C FCHISQ(NPTS,Y,SIGMAY,DAQ,WEIGHT,NFREE,MODE,YFIT) - evaluates C reduced chi-squared for fit C FDERIV(NPTS,X,I,A,DELTAA,NTERMS,DERIV) - evaluates the derivatives C of the fitting function for the Ith term with C respect to each parameter C MATINV(ARRAY,NTERMS,DET) - inverts a symmetric two dimensional C matrix of degree NTERMS and calculates its C determinant C IMPLICIT NONE INTEGER NPTS ! Number of input points to fit INTEGER DAQ(NPTS) ! Input array of data quality of points INTEGER BAD ! Lowest quality value indicating rejected data INTEGER NTERMS ! No. of variables for least squares fit INTEGER IFAIL ! Status flag for CURFIT call (1= <4 points; 2= ALPHA array has zero value diagonal elements) DOUBLE PRECISION X(NPTS) ! Array of input points - independent variable DOUBLE PRECISION Y(NPTS) ! Array of input points - dependent variable DOUBLE PRECISION WEIGHT(NPTS) ! Weight of input points - dependent variable DOUBLE PRECISION A(NTERMS) ! Output array of terms of least squares fit DOUBLE PRECISION DELTAA(NTERMS) ! Input array of deltas for numerical derivative determination DOUBLE PRECISION SIGMAA(NTERMS) ! Output array of standard deviations of terms of fit DOUBLE PRECISION FLAMDA ! Proportion of gradient search included DOUBLE PRECISION YFIT(NPTS) ! Ouput fitted values for dependent variable C C Local variables C INTEGER I INTEGER J INTEGER K INTEGER NJ INTEGER NFREE INTEGER STAT DOUBLE PRECISION ALPHA(3,3) DOUBLE PRECISION BETA(3) DOUBLE PRECISION DERIV(3) DOUBLE PRECISION ARRAY(3,3) DOUBLE PRECISION B(3) DOUBLE PRECISION DET DOUBLE PRECISION CHISQ1 DOUBLE PRECISION CHISQR DOUBLE PRECISION FUNCTN DOUBLE PRECISION FCHISQ EXTERNAL FUNCTN,FCHISQ C C Get number of degrees of freedom taking account of the data C quality and test C J=1 DO I=1,NPTS,1 IF (DAQ(I).GT.BAD) THEN J=J+1 ENDIF ENDDO NJ=J-1 NFREE=NPTS-NJ-NTERMS IF (NFREE.LE.0) THEN CALL UMSPUT(' Too few points for fit',3,0,STAT) CHISQR=0.0D0 IFAIL=1 GO TO 999 ENDIF IFAIL=0 C C Evaluate ALPHA and BETA matrices C DO J=1,NTERMS,1 BETA(J)=0.0D0 DO K=1,J,1 ALPHA(J,K)=0.0D0 ENDDO ENDDO DO I=1,NPTS,1 IF (DAQ(I).LE.BAD) THEN CALL FDERIV(NPTS,X,I,NTERMS,A,DELTAA,DERIV) DO J=1,NTERMS,1 BETA(J)=BETA(J) + : ( WEIGHT(I)*( Y(I)-FUNCTN(NPTS,X,I,NTERMS,A) )* : DERIV(J) ) DO K=1,J,1 ALPHA(J,K)=ALPHA(J,K) + (WEIGHT(I)*DERIV(J)*DERIV(K)) ENDDO ENDDO ENDIF ENDDO DO J=1,NTERMS,1 DO K=1,J,1 ALPHA(K,J)=ALPHA(J,K) ENDDO ENDDO C C Evaluate chi-squared at starting point C DO I=1,NPTS,1 IF (DAQ(I).LE.BAD) THEN YFIT(I)=FUNCTN(NPTS,X,I,NTERMS,A) ENDIF ENDDO CHISQ1=FCHISQ(NPTS,Y,DAQ,BAD,WEIGHT,NFREE,YFIT) IF (CHISQ1.EQ.0.0D0) THEN CALL UMSPUT(' Warning. Chi-squared = 0',3,0,STAT) ENDIF C C Invert modified curvature matrix to find new parameters C 71 DO J=1,NTERMS,1 DO K=1,NTERMS,1 IF (ALPHA(J,J).EQ.0.0D0.OR.ALPHA(K,K).EQ.0.0D0) THEN IFAIL=2 GO TO 999 ENDIF ARRAY(J,K)=ALPHA(J,K)/DSQRT(ALPHA(J,J)*ALPHA(K,K)) ENDDO ARRAY(J,J)=1.0D0 + FLAMDA ENDDO CALL MATINV(ARRAY,NTERMS,DET) DO J=1,NTERMS,1 B(J)=A(J) DO K=1,NTERMS,1 B(J)=B(J) + BETA(K)*ARRAY(J,K)/DSQRT(ALPHA(J,J)*ALPHA(K,K)) ENDDO ENDDO C C If chi-squared increased, increase Flamda and retry C DO I=1,NPTS,1 IF (DAQ(I).NE.BAD) THEN YFIT(I)=FUNCTN(NPTS,X,I,NTERMS,B) ENDIF ENDDO CHISQR=FCHISQ(NPTS,Y,DAQ,BAD,WEIGHT,NFREE,YFIT) IF (CHISQR.EQ.0.0D0) THEN CALL UMSPUT(' Warning. Chi-squared = 0',3,0,STAT) ENDIF IF (FLAMDA.EQ.0.0D0) THEN GO TO 101 ENDIF IF ((CHISQ1-CHISQR).LT.0.0D0) THEN FLAMDA=10.0D0*FLAMDA GO TO 71 ENDIF C C Evaluate parameters and errors C 101 DO J=1,NTERMS,1 A(J)=B(J) SIGMAA(J)=DSQRT(DABS(ARRAY(J,J)/ALPHA(J,J))) ENDDO FLAMDA=FLAMDA/10.0D0 999 END DOUBLE PRECISION FUNCTION FUNCTN(NP,X,IN,NTERMS,A) C C Function subprogram to evaluate terms of a function of the C form Y = coeffi*I + coeffQ.Q + coeffU.U C Taken from Bevington, 1969, p.214 C IMPLICIT NONE INTEGER NP ! No. of points in input array X INTEGER IN ! Index number of point being evalauted in array X INTEGER NTERMS ! No. of terms in function DOUBLE PRECISION X(NP) ! Array of independent variable at which to evaluate derivative DOUBLE PRECISION A(NTERMS) ! Array of terms of function C C Variables in common C INTEGER NPC ! 1st dimension of common arrays COEFQU and COEFUQ PARAMETER (NPC=64) DOUBLE PRECISION EFFI(NPC) ! Array of coefficients for Stokes I parameter term DOUBLE PRECISION EFFQ(NPC) ! Array of coefficients for Stokes Q parameter term (cos2*theta) DOUBLE PRECISION EFFU(NPC) ! Array of oefficients for Stokes U parameter term (sin2*theta) COMMON/LSQ/EFFI,EFFQ,EFFU FUNCTN=A(1)*EFFI(IN) + A(2)*EFFQ(IN) + A(3)*EFFU(IN) 99 END SUBROUTINE FDERIV(NP,X,IN,NTERMS,A,DELTAA,DERIV) C C Subroutine subprogram to evaluate the derivatives of C a function which is of the form C Y = effi*I + effq*U + effu*U C Adapted from Bevington, 1969, p.241 C IMPLICIT NONE INTEGER NP ! No. of points in input array X INTEGER IN ! Index number of point being evalauted in array X INTEGER NTERMS ! No. of terms in function DOUBLE PRECISION X(NP) ! Array of independent variable at which to evaluate derivative DOUBLE PRECISION A(NTERMS) ! Array of terms of function DOUBLE PRECISION DELTAA(NTERMS) ! Input array of deltas for numerical derivative determination DOUBLE PRECISION DERIV(NTERMS) ! Output array of derivatives of terms of function C C Variables in common C INTEGER NPC ! 1st dimension of common arrays COEFQU and COEFUQ PARAMETER (NPC=64) DOUBLE PRECISION EFFI(NPC) ! Array of coefficients for Stokes I parameter term DOUBLE PRECISION EFFQ(NPC) ! Array of coefficients for Stokes Q parameter term (cos2*theta) DOUBLE PRECISION EFFU(NPC) ! Array of oefficients for Stokes U parameter term (sin2*theta) COMMON/LSQ/EFFI,EFFQ,EFFU DERIV(1)=EFFI(IN) DERIV(2)=EFFQ(IN) DERIV(3)=EFFU(IN) 99 END DOUBLE PRECISION FUNCTION FCHISQ(NPTS,Y,DAQ,BAD,WEIGHT, : NFREE,YFIT) C C Function subprogram to evaluate reduced chi-squared for fit to C data with rejection of bad points C Adapted from Bevington, 1969, p. 194. C IMPLICIT NONE INTEGER NPTS ! Number of input points to fit INTEGER DAQ(NPTS) ! Input array of data quality of points INTEGER BAD ! Lowest quality value indicating rejected data INTEGER NFREE ! Value of no. of varaibles - no. of terms (=degrees of freedom) DOUBLE PRECISION Y(NPTS) ! Array of input points - dependent variable DOUBLE PRECISION WEIGHT(NPTS) ! Weight of input points - dependent variable DOUBLE PRECISION YFIT(NPTS) ! Ouput fitted values for dependent variable C C Local variables C INTEGER I INTEGER NS DOUBLE PRECISION CHISQ CHISQ=0.0D0 IF (NFREE.LE.0) THEN FCHISQ=0.0D0 GO TO 99 ENDIF C C Accumulate chi-squared C NS=0 DO I=1,NPTS,1 IF (DAQ(I).NE.BAD) THEN CHISQ=CHISQ + ( WEIGHT(I)*((Y(I)-YFIT(I))**2) ) ELSE NS=NS+1 ENDIF ENDDO C C Divide by number of degrees of freedom C IF ((NFREE-NS).GT.0) THEN FCHISQ=CHISQ/DBLE(NFREE-NS) ELSE FCHISQ=0.0D0 ENDIF 99 END SUBROUTINE MATINV(ARRAY,NORDER,DET) C C Subroutine subprogram to invert a symmetric matrix and calculate C its determinant C The input matrix ARRAY is replaced by its inverse C Taken from Bevington, 1969, p.302-3 C IMPLICIT NONE INTEGER NORDER ! Order of input and output arrays DOUBLE PRECISION ARRAY(NORDER,NORDER) ! Input array to be inverted=Inverted output array DOUBLE PRECISION DET ! Determinant of input matrix C C Local variables C INTEGER I INTEGER J INTEGER K INTEGER L INTEGER JK(6) INTEGER IK(6) DOUBLE PRECISION AMAX DOUBLE PRECISION SAVE DET=1.0D0 DO K=1,NORDER,1 C C Find largest element in rest of matrix C AMAX=0.0D0 21 DO I=K,NORDER,1 DO J=K,NORDER,1 IF (DABS(AMAX)-DABS(ARRAY(I,J)).LE.0.0D0) THEN AMAX=ARRAY(I,J) IK(K)=I JK(K)=J ENDIF ENDDO ENDDO C C Interchange rows and columns to put AMAX in ARRAY(K,K) C IF (AMAX.EQ.0.0D0) THEN DET=0.0D0 GO TO 999 ENDIF I=IK(K) IF ((I-K).LT.0) THEN GO TO 21 ENDIF IF ((I-K).EQ.0) THEN GO TO 51 ENDIF IF ((I-K).GT.0) THEN GO TO 43 ENDIF 43 DO J=1,NORDER,1 SAVE=ARRAY(K,J) ARRAY(K,J)=ARRAY(I,J) ARRAY(I,J)=-SAVE ENDDO 51 J=JK(K) IF ((J-K).LT.0) THEN GO TO 21 ENDIF IF ((J-K).EQ.0) THEN GO TO 61 ENDIF IF ((J-K).GT.0) THEN GO TO 53 ENDIF 53 DO I=1,NORDER,1 SAVE=ARRAY(I,K) ARRAY(I,K)=ARRAY(I,J) ARRAY(I,J)=-SAVE ENDDO C C Accumulate elements of inverse matrix C 61 DO I=1,NORDER,1 IF ((I-K).NE.0) THEN ARRAY(I,K)=-ARRAY(I,K)/AMAX ENDIF ENDDO 71 DO I=1,NORDER,1 DO J=1,NORDER,1 IF ((I-K).EQ.0) THEN GO TO 80 ENDIF IF ((J-K).EQ.0) THEN GO TO 80 ENDIF ARRAY(I,J)=ARRAY(I,J) + ( ARRAY(I,K)*ARRAY(K,J) ) 80 CONTINUE ENDDO ENDDO 81 DO J=1,NORDER,1 IF ((J-K).NE.0) THEN ARRAY(K,J)=ARRAY(K,J)/AMAX ENDIF 90 CONTINUE ENDDO ARRAY(K,K)=1.0D0/AMAX DET=DET*AMAX ENDDO C C Restore ordering of matrix C DO L=1,NORDER,1 K=NORDER - L + 1 J=IK(K) IF ((J-K).LE.0) THEN GO TO 111 ENDIF DO I=1,NORDER,1 SAVE=ARRAY(I,K) ARRAY(I,K)=-ARRAY(I,J) ARRAY(I,J)=SAVE ENDDO 111 I=JK(K) IF ((I-K).LE.0) THEN GO TO 130 ENDIF DO J=1,NORDER,1 SAVE=ARRAY(K,J) ARRAY(K,J)=-ARRAY(I,J) ARRAY(I,J)=SAVE ENDDO 130 CONTINUE ENDDO 999 END SUBROUTINE DEBIAS(POL,POLERR,POLDB) C C This subroutine returns the debiased value of the C fractional polarization POL as POLDB using the C estimate of the most probable value of POL from C the fit to the peak of the Rice distribution given by C solving Equ. A2 of Wardle & Kronberg (ApJ, 194, 249, C 1974). C An initial estimate of POLDB is given by determining C the Bessel functions of (POLobs*POLobs/POLERR*POLERR). C Then the initial estimate of POLDB allows improved C values for the Bessel functions of (POLobs*POLreal/ C POLERR*POLERR) to be obtained. The procedure is C iterated until convergence in successive values of C POLDB is achieved to 1% of POLERR C IMPLICIT NONE REAL POL ! Initial value of polarization REAL POLERR ! Polarization error REAL POLDB ! Output debiased polarization C C Local variables C DOUBLE PRECISION RR DOUBLE PRECISION POLO DOUBLE PRECISION BESSI0 DOUBLE PRECISION BESSI1 RR=DBLE(POL*POL/(POLERR*POLERR)) IF (RR.GT.81.0D0) THEN C C Values of Modified Bessell functions > 6E33, C so use asymptotic expression for POLDB C GO TO 50 ENDIF C C Calculate the debiassed value of polarization C POLDB=REAL(( (BESSI0(RR)/BESSI1(RR))* : DBLE(POL - (POLERR*POLERR/POL)) )) IF (POLDB.GT.0.0) THEN POLO=POLDB ELSE POLDB=0.0 GO TO 99 ENDIF C C Recalculate RR with the new value of POL C 30 RR=DBLE(POL*POLDB/(POLERR*POLERR)) IF (RR.GT.81.0D0) THEN C C Values of Modified Bessell functions > 6E33, C so use asymptotic expression for POLDB C GO TO 50 ENDIF C C Calculate the new debiassed value of polarization C POLDB=REAL(( (BESSI0(RR)/BESSI1(RR))* : DBLE(POL - (POLERR*POLERR/POL)) )) IF (POLDB.LE.0.0) THEN POLDB=0.0 GO TO 99 ENDIF IF (ABS(POLDB-POLO).LT.(0.01*POLERR)) THEN GO TO 99 ! Converged ELSE POLO=POLDB GO TO 30 ! Retry ENDIF C C If RR>81 then Modified Bessel function values > ^E33. C Use approximate (asymptotic) value for POLDB C 50 IF (((POL*POL)-(POLERR*POLERR)).GT.0.0) THEN POLDB=SQRT(( (POL*POL) - (POLERR*POLERR) )) ELSE POLDB=0.0 ENDIF 99 END DOUBLE PRECISION FUNCTION BESSI0(X) C C This is the Numerical Recipes function for calulating C the modified Bessell function of Io(x) for any X C IMPLICIT NONE DOUBLE PRECISION X ! Input value of independent variable C C Local variables C DOUBLE PRECISION AX DOUBLE PRECISION P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,Y SAVE P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9 DATA P1,P2,P3,P4,P5,P6,P7/1.0D0,3.5156229D0,3.0899424D0, :1.2067492D0,0.2659732D0,0.360768D-1,0.45813D-2/ DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,0.1328592D-1, :0.225319D-2,-0.157565D-2,0.916281D-2,-0.2057706D-1,0.2635537D-1, :-0.1647633D-1,0.392377D-2/ IF (ABS(X).LT.3.75) then Y=(X/3.75)**2 BESSI0=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))) ELSE AX=ABS(X) Y=3.75/AX BESSI0=(EXP(AX)/SQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y* : (Q7+Y*(Q8+Y*Q9)))))))) ENDIF 99 END DOUBLE PRECISION FUNCTION BESSI1(X) C C This is the Numerical Recipes function to compute the C modified Bessell function I1(x) for any x C IMPLICIT NONE DOUBLE PRECISION X ! Input value of independent variable C C Local variables C DOUBLE PRECISION AX DOUBLE PRECISION P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,Y SAVE P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9 DATA P1,P2,P3,P4,P5,P6,P7/0.5D0,0.87890594D0,0.51498869D0, :0.15084934D0,0.2658733D-1,0.301532D-2,0.32411D-3/ DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,-0.3988024D-1, :-0.362018D-2,0.163801D-2,-0.1031555D-1,0.2282967D-1,-0.2895312D-1, :0.1787654D-1,-0.420059D-2/ IF (ABS(X).LT.3.75) THEN Y=(X/3.75)**2 BESSI1=X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))) ELSE AX=ABS(X) Y=3.75/AX BESSI1=(EXP(AX)/SQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y* :(Q7+Y*(Q8+Y*Q9)))))))) IF (X.LT.0.0D0) THEN BESSI1=-BESSI1 ENDIF ENDIF 99 END SUBROUTINE INSPOL(OBQ,OBU,OBSIQ,OBSIU,INPO,INPA, : COQ,COU,COSIQ,COSIU) C C This subroutine corrects the observed Stokes parameters C OBQ and OBU, and errors OBSIQ and OBSIU, for the C instrumental polarization INPO at position angle INPA. C The corrected Stokes parameters COQ and COU, and errors C COSIQ and COSIU, are returned. C The treatment of the instrumental polarization follows C the Appendix of Goodrich, ApJ, 311, 882, 1986 C IMPLICIT NONE REAL OBQ ! Input observed Stokes Q parameter REAL OBU ! Input observed Stokes U parameter REAL OBSIQ ! Input observed Stokes Q parameter standard deviation REAL OBSIU ! Input observed Stokes U parameter standard deviation REAL INPO ! Fractional instrumental polarization REAL INPA ! position angle of instrumental polarization REAL COQ ! Output instrument corrected Stokes Q parameter REAL COU ! Output instrument corrected Stokes U parameter REAL COSIQ ! Output instrument corrected Stokes Q parameter standard deviation REAL COSIU ! Output instrument corrected Stokes U parameter standard deviation C C Local variables C REAL A REAL B REAL C REAL D REAL E REAL F REAL CC REAL SS REAL SINRD REAL COSRD CC=COSRD(2.*INPA) SS=SINRD(2.*INPA) A=(1.0-INPO) B=-INPO*CC*(1.0-INPO) C=-INPO*SS*(1.0-INPO) D=(1.0 - INPO*CC*CC - INPO*INPO*SS*SS) E=-INPO*CC*SS*(1.0-INPO) F=(1.0 - INPO*SS*SS - INPO*INPO*CC*CC) COQ=(OBQ*(D/A) - INPO*CC*SS*OBU - INPO*CC)/ : (1.0 - INPO*CC*OBQ - INPO*SS*OBU) COU=(OBU*(F/A) - INPO*CC*SS*OBQ - INPO*SS)/ : (1.0 - INPO*CC*OBQ - INPO*SS*OBU) COSIQ=SQRT( ( (((A*D-B*B)+(C*D-B*E)*OBU)**2.*OBSIQ*OBSIQ) : + (((A*E-B*C)+(B*E-C*D)*OBQ)**2.*OBSIU*OBSIU) )/ :((A+B*OBQ+C*OBU)**4.) ) COSIU=SQRT( ( (((A*E-B*C)+(C*E-B*F)*OBU)**2.*OBSIQ*OBSIQ) : + (((A*F-C*C)+(B*F-C*E)*OBQ)**2.*OBSIU*OBSIU) )/ :((A+B*OBQ+C*OBU)**4.) ) 99 END SUBROUTINE PPASIG(P0,SIGP0,TH0,SIGTH) C C This subroutine returns the value of sigma(theta) (v) C for the error on the polarization position angle by C numerical integration of the variance of the distribution C of theta. The form of the distribution function is given C by Naghizadeh-Khouei and Clarke (A&A, 274, 968, 1993), C Equations 5 and 3. C Romberg integration is employed to integrate the function. C IMPLICIT NONE REAL P0 ! Input linear polarization REAL SIGP0 ! Input linear polarization error REAL TH0 ! Input value of polarization position angle REAL SIGTH ! Output value of polarization position angle error C C Local variables C DOUBLE PRECISION P00 DOUBLE PRECISION TH00 DOUBLE PRECISION PIE DOUBLE PRECISION PI2 DOUBLE PRECISION STH PARAMETER (PIE=3.1415926535898D0) PI2=PIE/2.0D0 P00=DBLE(P0/SIGP0) TH00=0.0D0 ! In radians CALL THROMB(TH00,P00,-PI2,PI2,STH) SIGTH=REAL(DSQRT(ABS(STH))*180.0D0/PIE) ! In degrees 99 END SUBROUTINE THROMB(TH0,P0,A,B,SS) C C This is the Numerical Recipes subroutine to perform C Romberg integration of the variance of the distribution C function for the polarization position angle TH0 over the C range A to B C IMPLICIT NONE DOUBLE PRECISION TH0 ! Input polarization position angle DOUBLE PRECISION P0 ! Input linear polarization DOUBLE PRECISION A ! Lower limit of numerical integration DOUBLE PRECISION B ! Upper limit of numerical integration DOUBLE PRECISION SS ! Output value of integrand C C Local variables C INTEGER J INTEGER JMAX INTEGER JMAXP INTEGER K INTEGER KM DOUBLE PRECISION EPS DOUBLE PRECISION S(21) DOUBLE PRECISION H(21) DOUBLE PRECISION DSS PARAMETER (EPS=1.D-6, :JMAX=12, ! Maximum number of interior points ( 2**(JMAX-2) ) to add :JMAXP=JMAX+1, :K=5, :KM=K-1) H(1)=1. DO J=1,JMAX,1 CALL TRAPZD(TH0,P0,A,B,S(J),J) IF (J.GE.K) THEN CALL POLINT(H(J-KM),S(J-KM),K,0.0D0,SS,DSS) IF (ABS(DSS).LT.EPS*ABS(SS)) THEN GO TO 99 ENDIF ENDIF S(J+1)=S(J) H(J+1)=0.25D0*H(J) ENDDO 99 END SUBROUTINE TRAPZD(TH0,P0,A,B,S,N) C C This is the Numerical Recipes extended Trapezoidal C integration routine to calculate the variance of the C distribution function for the polarization position angle C over the range A to B C IMPLICIT NONE INTEGER N ! No. of steps for numerical integration DOUBLE PRECISION TH0 ! Input polarization position angle DOUBLE PRECISION P0 ! Input linear polarization DOUBLE PRECISION A ! Lower limit of numerical integration DOUBLE PRECISION B ! Upper limit of numerical integration DOUBLE PRECISION S ! Output value of integrand C C Local variables C INTEGER J INTEGER IT DOUBLE PRECISION TNM DOUBLE PRECISION DEL DOUBLE PRECISION SUM DOUBLE PRECISION X DOUBLE PRECISION FUNC EXTERNAL FUNC IF (N.EQ.1) THEN S=0.50D0*(B - A)*(FUNC(P0,TH0,A) + FUNC(P0,TH0,B)) ELSE IT=2**(N-2) TNM=DBLE(IT) DEL=(B - A)/TNM X=A + 0.50D0*DEL SUM=0. DO J=1,IT,1 SUM=SUM + FUNC(P0,TH0,X) X=X+DEL ENDDO S=0.50D0*(S + ((B - A)*SUM/TNM)) ENDIF 99 END SUBROUTINE POLINT(XA,YA,N,X,Y,DY) C C This is the Numerical Recipes polynomial interpolation routine C IMPLICIT NONE INTEGER N ! No. of input values DOUBLE PRECISION XA(N) ! Input array - independent variable DOUBLE PRECISION YA(N) ! Input array - dependent variable DOUBLE PRECISION X ! Input value of independent variable to interpolate DOUBLE PRECISION Y ! Output interpolated value of dependent variable DOUBLE PRECISION DY ! Error estimate on interpolated value C C Local variables C INTEGER I INTEGER NMAX INTEGER NS INTEGER M INTEGER STAT DOUBLE PRECISION DIF DOUBLE PRECISION DIFT DOUBLE PRECISION C(10) DOUBLE PRECISION D(10) DOUBLE PRECISION HO DOUBLE PRECISION HP DOUBLE PRECISION W DOUBLE PRECISION DEN PARAMETER (NMAX=10) NS=1 DIF=ABS(X-XA(1)) DO I=1,N,1 DIFT=ABS(X-XA(I)) IF (DIFT.LT.DIF) THEN NS=I DIF=DIFT ENDIF C(I)=YA(I) D(I)=YA(I) ENDDO Y=YA(NS) NS=NS-1 DO M=1,N-1,1 DO I=1,N-M,1 HO=XA(I)-X HP=XA(I+M)-X W=C(I+1)-D(I) DEN=HO-HP IF (DEN.EQ.0.0D0) THEN CALL UMSPUT(' Problems calculating sigam(theta)', : 3,0,STAT) GO TO 99 ENDIF DEN=W/DEN D(I)=HP*DEN C(I)=HO*DEN ENDDO IF (2*NS.LT.N-M) THEN DY=C(NS+1) ELSE DY=D(NS) NS=NS-1 ENDIF Y=Y+DY ENDDO 99 END DOUBLE PRECISION FUNCTION FUNC(P0,TH0,TH) C C This function returns the value of the function for C the variance of the polarization position angle C distribution for 1sigma. It is taken from C Naghizadeh-Khouei and Clarke (1993), Equations 5 C and 3. C IMPLICIT NONE DOUBLE PRECISION P0 ! Input linear polarization DOUBLE PRECISION TH0 ! Input reference position angle DOUBLE PRECISION TH ! Input position angle for variance determination C C Local varaibles C DOUBLE PRECISION PI DOUBLE PRECISION ETA DOUBLE PRECISION ERF EXTERNAL ERF PARAMETER (PI=3.1415926535898D0) ETA=(P0/DSQRT(2.0D0))*DCOS(2.0D0*(TH-TH0)) FUNC=(1.0D0/DSQRT(PI))*( (1.0D0/DSQRT(PI)) + : (ETA*DEXP(ETA**2)*(1.0D0+ERF(ETA))) )*DEXP(-((P0**2)/2.0D0))* :(TH-TH0)**2 99 END DOUBLE PRECISION FUNCTION ERF(X) C C This is the series expansion approximation to evaluate the C Error Function (Erf) C IMPLICIT NONE DOUBLE PRECISION X ! Value of independent variable C INTEGER I,N DOUBLE PRECISION PI DOUBLE PRECISION TERM DOUBLE PRECISION SUM DOUBLE PRECISION NFAC PI=3.1415926535898D0 N=100 NFAC=1.0D0 SUM=X DO I=1,N,1 NFAC=DBLE(I)*NFAC TERM=(-1.0D0)**(I) * X**(2*I+1) / (NFAC*DBLE(2*I+1)) IF (ABS(TERM/SUM).LT.1.0D-30) THEN GO TO 50 ENDIF SUM=SUM+TERM ENDDO 50 ERF=(2.0D0/DSQRT(PI))*SUM 99 END REAL FUNCTION NORDEV(IDUM,MEAN,SIGMA) C C This is the Numerical Recipes function to return a normally C distributed deviate with mean value MEAN and sigma SIGMA, C using RAN1(IDUM) as a source of uniform deviates C IMPLICIT NONE INTEGER IDUM REAL MEAN,SIGMA INTEGER ISET REAL FAC,GSET,RSQ,V1,V2,RAN1,GASDEV SAVE ISET,GSET DATA ISET/0/ IF (ISET.EQ.0) THEN 10 V1=2.*RAN1(IDUM) - 1.0 V2=2.*RAN1(IDUM) - 1.0 RSQ=V1**2 + V2**2 IF (RSQ.GE.1.0.OR.RSQ.EQ.0.0) THEN GO TO 10 ENDIF FAC=SQRT(-2.*LOG(RSQ)/RSQ) GSET=V1*FAC GASDEV=V2*FAC ISET=1 ELSE GASDEV=GSET ISET=0 ENDIF C C GASDEV has zero mean and unit variance - convert C to a value with MEAN and SIGMA C NORDEV=MEAN + SIGMA*GASDEV 99 END REAL FUNCTION RAN1(IDUM) C C This is the Numerical Recipes random number generator C of Park amd Miller with Bays-Durham shuffle safeguards. C IMPLICIT NONE INTEGER IDUM INTEGER IA,IM,IQ,IR,NTAB,NDIV REAL AM,EPS,RNMX PARAMETER (IA=16807,IM=2147483647,AM=1./IM, :IQ=127773,IR=2836,NTAB=32,NDIV=1+(IM-1)/NTAB, :EPS=1.2E-7,RNMX=1.-EPS) INTEGER J,K,IV(NTAB),IY SAVE IV,IY DATA IV /NTAB*0/, IY /0/ IF (IDUM.LE.0.OR.IY.EQ.0) then IDUM=MAX(-IDUM,1) DO 11 J=NTAB+8,1,-1 K=IDUM/IQ IDUM=IA*(IDUM-K*IQ)-IR*K IF (IDUM.LT.0) THEN IDUM=IDUM+IM ENDIF IF (J.LE.NTAB) THEN IV(J)=IDUM ENDIF 11 CONTINUE IY=IV(1) ENDIF K=IDUM/IQ IDUM=IA*(IDUM-K*IQ)-IR*K IF (IDUM.LT.0) THEN IDUM=IDUM+IM ENDIF J=1+IY/NDIV IY=IV(J) IV(J)=IDUM RAN1=MIN(AM*IY,RNMX) 99 END DOUBLE PRECISION FUNCTION SINDD(X) C C This function returns the sine of X, where X is C in degrees. C It is required since LINUX does not support SIND C IMPLICIT NONE DOUBLE PRECISION X C C Local variables C DOUBLE PRECISION PIE PARAMETER (PIE=3.1415926535898D0) SINDD=SIN(X*PIE/180.0D0) 99 END DOUBLE PRECISION FUNCTION COSDD(X) C C This function returns the cosine of X, where X is C in degrees. C It is required since LINUX does not support COSD C IMPLICIT NONE DOUBLE PRECISION X C C Local variables C DOUBLE PRECISION PIE PARAMETER (PIE=3.1415926535898D0) COSDD=COS(X*PIE/180.0D0) 99 END REAL FUNCTION SINRD(X) C C This function returns the sine of X, where X is C in degrees. C It is required since LINUX does not support SIND C IMPLICIT NONE REAL X C C Local variables C REAL PIE PARAMETER (PIE=3.1415926535898) SINRD=SIN(X*PIE/180.0) 99 END REAL FUNCTION COSRD(X) C C This function returns the cosine of X, where X is C in degrees. C It is required since LINUX does not support COSD C IMPLICIT NONE REAL X C C Local variables C REAL PIE PARAMETER (PIE=3.1415926535898) COSRD=COS(X*PIE/180.0) 99 END