C----------------------------------------------------------------------- C; Copyright (C) 1995,1996,1997,1998,2001 C; Associated Universities, Inc. Washington DC, USA. C; C; This program is free software; you can redistribute it and/or C; modify it under the terms of the GNU General Public License as C; published by the Free Software Foundation; either version 2 of C; the License, or (at your option) any later version. C; C; This program is distributed in the hope that it will be useful, C; but WITHOUT ANY WARRANTY; without even the implied warranty of C; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C; GNU General Public License for more details. C; C; You should have received a copy of the GNU General Public C; License along with this program; if not, write to the Free C; Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, C; MA 02139, USA. C C For further information about this software contact: C Internet email: bcotton@nrao.edu. C Postal address: Bill Cotton C National Radio Astronomy Observatory C 520 Edgemont Road C Charlottesville, VA 22903-2475 USA C----------------------------------------------------------------------- SUBROUTINE CORERR (RA, DEC, PEAK, MAJOR, MINOR, POSANG, * QCENT, UCENT, PFLUX, IRMS, PRMS, BEAM, FITTED, DORAW, FBLANK, * FLUX, EFLUX, EPFLUX, CHPANG, CHEPAN, * ERRRA, ERRDEC, CMAJOR, CMINOR, CPA, EMAJOR, EMINOR, EPA) C----------------------------------------------------------------------- C Routine to apply any corrections to NVSS fitted source parameters C and determine errors. Unless FITTED, the model parameters returned C are deconvolved using the restoring BEAM. C Now allows for blanking of PFLUX, QCENT and UCENT C Now computes effective NVSS point source response rather than using C BEAM. C Inputs: C RA D RA (deg) C DEC D Dec (deg) C PEAK R Peak Ipol (Jy/beam) C MAJOR R Fitted major axis size C MINOR R Fitted minor axis size C POSANG R Fitted PA C QCENT R Center Q flux density C UCENT R Center U flux density C PFLUX R Integrated polarized flux density C IRMS R RMS (sigma) in Ipol. C PRMS R RMS (sigma) in Qpol and Upol. C BEAM R(3) Restoring beam major, minor axes and position angle C (all deg) (Now ignored) C FITTED L If true return fitted values else deconvolved C DORAW L If true return raw values else bias corrected. C Generally, FITTED should be true if DORAW is. C FBLANK R Magic value blanking value C Outputs: C FLUX R Model peak/integrated Flux density (Jy) C EFLUX R Error in FLUX C PFLUX R Polarized flux density (Jy) C Now derived from QCENT, UCENT C EPFLUX R Error in PFLUX C CHPANG C*6 Polarization angle or blank if probabliity of a C false detection exceeds 2% C CHEPAN C*6 Error in CHPANG or blank C ERRRA R Error (sec of time) of Right Ascension C ERRDEC R Error (asec) of Declination C CMAJOR C*6 Major axis size or limit as string (asec) C CMINOR C*6 Minor axis size or limit as string (asec) C CPA C*6 Position angle or blank as string (asec) C EMAJOR C*6 Error of major axis size or limit as string (asec) C EMINOR C*6 Error of Minor axis size or limit as string (asec) C EPA C*6 Error of position angle or blank as string (asec) C----------------------------------------------------------------------- DOUBLE PRECISION RA, DEC REAL PEAK, MAJOR, MINOR, POSANG, QCENT, UCENT, PFLUX, IRMS, * PRMS, BEAM(3), FBLANK, FLUX, EFLUX, EPFLUX LOGICAL FITTED, DORAW CHARACTER CHPANG*6, CHEPAN*6, CMAJOR*6, CMINOR*6, CPA*6, * EMAJOR*6, EMINOR*6, EPA*6 C REAL BEMRAT, SNR, SNRAMP, SNRMAJ, SNRMIN, ERRPEK, ERRPK2, * SINC, COSC, TMAJ, TMIN, ERRMAJ, ERRMIN, ERRPA, ERRX2, ERRY2, * ERRRA, ERRDEC, TFLUX, ETFLUX, PAMP, * ERRPAN, POLANG, PEAKP, PEAKX, PERR, * DECRAD, SPMJY, PSF(3) INTEGER IER LOGICAL RESOLV, HAFRES(3) C Calibration component of C position error (squared) REAL CALERA, CALEDE PARAMETER (CALERA = (0.45 / 3600.0) ** 2) PARAMETER (CALEDE = (0.56 / 3600.0) ** 2) C Mean and uncertainty of bias REAL BIASAV, BIASER PARAMETER (BIASAV = 0.00030) PARAMETER (BIASER = 0.00030) C Amplitude calibration C uncertainty REAL CALAER PARAMETER (CALAER = 0.03) C Polarization calibration error C (Wild guess) REAL CALPER PARAMETER (CALPER = 0.003) C Position bias DOUBLE PRECISION BIASRA, BIASDE PARAMETER (BIASRA = -0.025D0/3600.0D0) PARAMETER (BIASDE = 0.113D0/3600.0D0) C Degrees per radian DOUBLE PRECISION DG2RAD PARAMETER (DG2RAD = 1.745329252E-2) C----------------------------------------------------------------------- C Compute effective PSF DECRAD = DEC * DG2RAD SPMJY = PEAK * 1000.0 CALL NEWPS (DECRAD, SPMJY, PSF(1), PSF(2), PSF(3)) C Convert sizes to degrees PSF(1) = PSF(1) / 3600.0 PSF(2) = PSF(2) / 3600.0 C Force fitted value to be at C least as large as the PSF MAJOR = MAX (MAJOR, PSF(1)) MINOR = MAX (MINOR, PSF(2)) C Save peak Flux density PEAKX = PEAK C Correct position bias IF (.NOT.DORAW) THEN RA = RA + BIASRA / COS (DEC * DG2RAD) DEC = DEC + BIASDE END IF C Beam ratio BEMRAT = (MAJOR/PSF(1)) * (MINOR/PSF(2)) C Trap under size beams BEMRAT = MAX (BEMRAT, 1.0) C Effective SNRs^2 to account for C correlated noise. SNR = PEAKX / IRMS C SNR**2 for amplitude errors SNRAMP = (MAJOR*MINOR/(4.0*PSF(1)*PSF(1))) * * ((1.0+(PSF(1)/MAJOR)**2)**1.5) * * ((1.0+(PSF(1)/MINOR)**2)**1.5) * * SNR * SNR C SNR**2 for major axis error SNRMAJ = (MAJOR*MINOR/(4.0*PSF(1)*PSF(1))) * * ((1.0+(PSF(1)/MAJOR)**2)**2.5) * * ((1.0+(PSF(1)/MINOR)**2)**0.5) * * SNR * SNR C SNR**2 for minor axis/PA errors SNRMIN = (MAJOR*MINOR/(4.0*PSF(1)*PSF(1))) * * ((1.0+(PSF(1)/MAJOR)**2)**0.5) * * ((1.0+(PSF(1)/MINOR)**2)**2.5) * * SNR * SNR C Bias correct peak flux Cold IF (.NOT.DORAW) PEAKX = PEAKX + BIASAV - 1.0E-3*(BIASFA / PEAK IF (.NOT.DORAW) PEAKX = PEAKX + BIASAV - (IRMS * IRMS / PEAKX) C Flux errors C Add bias uncertainty ERRPEK = SQRT ((2.0 * PEAKX*PEAKX / SNRAMP) + BIASER**2) C Add Flux calibration (3%) error ERRPEK = SQRT (ERRPEK*ERRPEK + (CALAER * PEAKX)**2) C C Polarization flux, angle and C error IF ((QCENT.NE.FBLANK) .AND. (UCENT.NE.FBLANK)) THEN POLANG = 28.6479 * ATAN2 (UCENT, QCENT+1.0E-20) PAMP = SQRT (QCENT*QCENT + UCENT*UCENT) C Bias correction IF (.NOT.DORAW) CALL PDBIAS (PAMP, PRMS) C PAMP = MAX (1.0E-15, PAMP) ERRPAN = 28.6479 * PRMS / MAX (1.0E-15, PAMP) ERRPAN = MIN (ERRPAN, 90.0) ELSE PAMP = 0.0 POLANG = 0.0 ERRPAN = -1.0 END IF C Use central polarization rather C than integral PFLUX = PAMP IF (.NOT.(FITTED.OR.DORAW)) THEN PFLUX = PFLUX * BEMRAT PRMS = PRMS * BEMRAT END IF C Bad hair? (no polarization) IF ((QCENT.EQ.FBLANK) .OR. (UCENT.EQ.FBLANK)) PFLUX = FBLANK C Only quote Polarization angle if C Poln SNR>2.6 C Poln uncertainty adding C calibration term PERR = SQRT (PRMS*PRMS * (CALPER * PEAK)**2) C Probability of false detection < IF ((PFLUX.NE.FBLANK) .AND. ((PFLUX/PERR).GT.2.6)) THEN WRITE(CHPANG,1050) POLANG WRITE(CHEPAN,1050) ERRPAN ELSE IF (PFLUX.EQ.FBLANK) THEN CHPANG = 'Blank ' CHEPAN = ' ' ELSE CHPANG = ' ' CHEPAN = ' ' END IF C Put position angle in range (+/- C 90 deg) IF (POSANG.LT.-90.0) POSANG = POSANG + 180.0 IF (POSANG.GT.90.0) POSANG = POSANG - 180.0 C Error in polarized flux EPFLUX = PRMS * 1.41421 * BEMRAT C Errors SINC = SIN (DG2RAD*POSANG) COSC = COS (DG2RAD*POSANG) C Trap under size beams TMAJ = MAX (MAJOR, PSF(2)) TMIN = MAX (MINOR, PSF(2)) C Axis sizes include 2% C calibration error. ERRMAJ = SQRT (((2.0 * TMAJ**2) / SNRMAJ) + (0.02*PSF(1))**2) ERRMIN = SQRT (((2.0 * TMIN**2) / SNRMIN) + (0.02*PSF(1))**2) ERRPA = 57.296 * SQRT ((TMAJ*TMIN/ * (TMAJ*TMAJ - TMIN*TMIN + 1.0E-20))**2 * 4.0 / SNRMIN) ERRPA = MIN (ERRPA, 90.0) C Position errors ERRX2 = 2.0 * (2.0 * TMAJ**2) / (8.0 * LOG (2.0) * SNRMAJ) ERRY2 = 2.0 * (2.0 * TMIN**2) / (8.0 * LOG (2.0) * SNRMIN) C Include calibration ERRRA = SQRT (CALERA + ERRX2*SINC*SINC + * ERRY2*COSC*COSC) ERRDEC = SQRT (CALEDE + ERRY2*SINC*SINC + * ERRX2*COSC*COSC) C RA error in seconds of time ERRRA = ERRRA / (15.0 * COS (DG2RAD * DEC)) C Deconvolve CALL BMVAL (MAJOR, MINOR, POSANG, ERRMAJ, ERRMIN, * ERRPA, PSF(1), PSF(2), PSF(3), * CMAJOR, CMINOR, CPA, EMAJOR, EMINOR, EPA, RESOLV, HAFRES, * IER) C Convert to strings IF (FITTED) THEN C Fitted values WRITE(CMAJOR,1050) MAJOR * 3600.0 WRITE(CMINOR,1050) MINOR * 3600.0 WRITE(CPA,1050) POSANG WRITE(EMAJOR,1050) ERRMAJ * 3600.0 WRITE(EMINOR,1050) ERRMIN * 3600.0 WRITE(EPA,1050) ERRPA IF (RESOLV.OR.DORAW) THEN FLUX = PEAKX EFLUX = ERRPEK ELSE C Correct if unresolved FLUX = PEAKX * SQRT (BEMRAT) EFLUX = SQRT (((CALAER * FLUX)**2) + * (FLUX*FLUX/SNRAMP) + BIASER**2) END IF ELSE C Total flux density C Trap resolved on one axis IF (HAFRES(1)) THEN C Pseudo peak flux IF (HAFRES(2)) THEN C Major axis only PEAKP = PEAKX * SQRT (MAX (1.0, MINOR/PSF(2))) TFLUX = PEAKP * MAX (1.0, MAJOR/PSF(1)) ELSE C Minor axis only PEAKP = PEAKX * SQRT (MAX (1.0, MAJOR/PSF(1))) TFLUX = PEAKP * MAX (1.0, MINOR/PSF(2)) END IF C Error in pseudo peak flux ERRPK2 = SQRT (((CALAER * PEAKP)**2) + * (3.0 * PEAKP*PEAKP / (2.0 * SNRAMP)) + * BIASER**2) C Integrated flux error ETFLUX = SQRT (TFLUX*TFLUX * * ((ERRPK2*ERRPK2 / (PEAKP*PEAKP)) + * (PSF(1)/MAJOR) * ((ERRMAJ*ERRMAJ/(MAJOR*MAJOR))))) ELSE C Fully resolved values TFLUX = PEAKX * BEMRAT ETFLUX = SQRT (TFLUX*TFLUX * * ((ERRPEK*ERRPEK/(PEAKX*PEAKX)) + * (1.0 / BEMRAT) * ((ERRMAJ*ERRMAJ/(MAJOR*MAJOR)) + * (ERRMIN*ERRMIN/(MINOR*MINOR))))) END IF IF (RESOLV) THEN FLUX = TFLUX EFLUX = ETFLUX ELSE C Correct if unresolved FLUX = PEAKX * SQRT (BEMRAT) EFLUX = SQRT (((CALAER * FLUX)**2) + * (FLUX*FLUX/SNRAMP) + BIASER**2) END IF END IF C Convert units for output ERRRA = ERRRA * 3600.0 ERRDEC = ERRDEC * 3600.0 C 999 RETURN C----------------------------------------------------------------------- 1050 FORMAT (F6.1) 1051 FORMAT ('<', F5.1) END SUBROUTINE PDBIAS (P, RMS) C----------------------------------------------------------------------- C Estimates the polarization bias in a polarization amplitude, P, C measured in the presence of Q and U RMS noise, RMS. Returns the C corrected value. C The bias correction is such that the average bias is removed; C thus the average in the absence of a signal is zero. Does table C look of values calculated by J. Condon in the range of P/RMS of C 1.253 (the no signal limit) and 4 (the high SNR regime). Does C second order Lagrange interpolation. At lower values of P/RMS the C bias is a constant 1.253*RMS. Above a signal-to-noise ratio of 4, C use the formula: C normalized bias = 1 / (2 * s) + 1 / (8 * s**3), C where s is the true normalized flux density, iterating once to C estimate s from the normalized map flux density. "Normalized" means C divided by the rms noise in the q and u maps. C Inputs: C P R On input, P is the observed total polarized intensity C RMS R The standard deviation of the (assumed equal) C Gaussian distributions of the Stokes Q or U maps. C Output: C P R On output, P is the estimated intrinsic total C polarized intensity. C----------------------------------------------------------------------- REAL P, RMS C INTEGER I, INDEX, I1, I2, I3 REAL TABLE(2,40), PNORM, BIAS, D1, D2, D3, WT1, WT2, WT3, * SUM, SUMWT C (map_flux,map_bias) pairs DATA TABLE / * 1.253,1.253, 1.256,1.156, 1.266,1.066, 1.281,0.9814, * 1.303,0.9030, 1.330,0.8304, 1.364,0.7636, 1.402,0.7023, * 1.446,0.6462, 1.495,0.5951, 1.549,0.5486, 1.606,0.5064, * 1.668,0.4683, 1.734,0.4339, 1.803,0.4028, 1.875,0.3749, * 1.950,0.3498, 2.027,0.3273, 2.107,0.3070, 2.189,0.2888, * 2.272,0.2724, 2.358,0.2576, 2.444,0.2442, 2.532,0.2321, * 2.621,0.2212, 2.711,0.2112, 2.802,0.2021, 2.894,0.1938, * 2.986,0.1861, 3.079,0.1791, 3.173,0.1726, 3.267,0.1666, * 3.361,0.1610, 3.456,0.1557, 3.551,0.1509, 3.646,0.1463, * 3.742,0.1420, 3.838,0.1380, 3.934,0.1342, 4.031,0.1306/ C----------------------------------------------------------------------- C Check RMS IF (RMS.LE.0.0) GO TO 999 PNORM = P / RMS C Which regime? IF (PNORM.LE.TABLE(1,1)) THEN C Low (no) SNR case P = P - TABLE(2,1) * RMS ELSE IF (PNORM.GE.TABLE(1,40)) THEN C High SNR BIAS = 1.0 / (2.0 * PNORM) + 1.0 / (8.0 * PNORM**3) PNORM = PNORM - BIAS BIAS = 1.0 / (2.0 * PNORM) + 1.0 / (8.0 * PNORM**3) C Correct for bias P = P - BIAS * RMS ELSE C Middle, interpolate in table INDEX = 2 DO 20 I = 3,39 IF (PNORM.LT.TABLE(1,I)) GO TO 30 INDEX = I 20 CONTINUE C Lagrange interpolation 30 I1 = INDEX - 1 I2 = INDEX I3 = INDEX + 1 D1 = (TABLE(1,I1) - TABLE(1,I2)) * (TABLE(1,I1) - TABLE(1,I3)) D2 = (TABLE(1,I2) - TABLE(1,I1)) * (TABLE(1,I2) - TABLE(1,I3)) D3 = (TABLE(1,I3) - TABLE(1,I1)) * (TABLE(1,I3) - TABLE(1,I2)) WT1 = (PNORM - TABLE(1,I2)) * (PNORM - TABLE(1,I3)) / D1 WT2 = (PNORM - TABLE(1,I1)) * (PNORM - TABLE(1,I3)) / D2 WT3 = (PNORM - TABLE(1,I1)) * (PNORM - TABLE(1,I2)) / D3 SUM = TABLE(2,I1) * WT1 + TABLE(2,I2) * WT2 + TABLE(2,I3) * WT3 SUMWT = WT1 + WT2 + WT3 IF (SUMWT.GT.0.0) THEN BIAS = SUM / SUMWT ELSE C Shouldn't ever get here but do C something reasonable. BIAS = TABLE(2,I2) END IF C Correct for bias P = P - BIAS * RMS END IF C 999 RETURN END SUBROUTINE BMVAL (BMAJ, BMIN, BPA, BMAJE, BMINE, BPAE, CBMAJ, * CBMIN, CBPA, CMAJOR, CMINOR, CPA, EMAJOR, EMINOR, EPA, RESOLV, * HAFRES, IER) C----------------------------------------------------------------------- C Subroutine BMVAL deconvolves the fitted beam from the clean beam and C also generates appropriate errors. C Note: this only really works properly for round beams. C Inputs: C BMAJ R Fitted major axis C BMIN R Fitted minor axis C BPA R Fitted pos. angle (deg) C BMAJE R Fitted major axis error C BMINE R Fitted minor axis error C BPAE R Fitted pos. angle error (deg) C CBMAJ R Clean beam major axis C CBMIN R Clean beam minor axis C CBPA R Clean beam pos. angle (deg) C Outputs: C CMAJOR C*6 string with Major axis in asec C CMINOR C*6 string with Minor axis in asec C CPA C*6 string with position angle in deg C EMAJOR C*6 string with error of Major axis in asec C EMINOR C*6 string with error of Minor axis in asec C EPA C*6 string with error of position angle in deg C RESOLV L If true source was resolved. C HAFRES L(3) If true source was resolved on: C 1) one axis C 2) Major axis only C 3) Minor axis only C IER I Error return 0-> Can completely deconvolve C 1-> Cannot deconv some limits C 2-> Cannot deconv fitted source C----------------------------------------------------------------------- REAL BMAJ, BMIN, BPA, BMAJE, BMINE, BPAE, CBMAJ, CBMIN, CBPA CHARACTER CMAJOR*6, CMINOR*6, CPA*6, EMAJOR*6, EMINOR*6, EPA*6 LOGICAL RESOLV, HAFRES(3) INTEGER IER C REAL R(3,3), ERRPA, SRCMAJ, USRCMJ, SRCMIN, USRCMN, SRCPA C----------------------------------------------------------------------- C Get deconvolved sizes and C errors. CALL DECONV (CBMAJ, CBMIN, CBPA, * BMAJ, BMAJE, BMIN, BMINE, BPA, BPAE, * SRCMAJ, USRCMJ, SRCMIN, USRCMN, SRCPA) R(1,1) = SRCMAJ R(1,2) = USRCMJ R(1,3) = USRCMJ R(2,1) = SRCMIN R(2,2) = USRCMN R(2,2) = USRCMN R(3,1) = SRCPA R(3,2) = BPAE R(3,3) = BPAE RESOLV = .FALSE. HAFRES(1) = .FALSE. HAFRES(2) = .FALSE. HAFRES(3) = .FALSE. C Convert to strings 75 IF ((R(1,2).GT.0.0) .AND. (R(2,2).GT.0.0)) THEN C Deconvolved both axes WRITE(CMAJOR,1050) R(1,1) * 3600.0 WRITE(CMINOR,1050) R(2,1) * 3600.0 WRITE(CPA,1050) R(3,1) WRITE(EMAJOR,1050) (R(1,2) + R(1,3)) * 0.5 * 3600.0 WRITE(EMINOR,1050) (R(2,2) + R(2,3)) * 0.5 * 3600.0 ERRPA = MIN (90.0, 0.5 * ABS (R(3,2) + R(3,3))) WRITE(EPA,1050) ERRPA RESOLV = .TRUE. ELSE IF ((R(1,2).LT.0.0) .AND. (R(2,2).LT.0.0)) THEN C Upper limits only IF ((R(1,1)*3600.0).LE.99.0) THEN WRITE(CMAJOR,1051) R(1,1) * 3600.0 ELSE WRITE(CMAJOR,1052) R(1,1) * 3600.0 END IF IF ((R(2,1)*3600.0).LE.99.0) THEN WRITE(CMINOR,1051) R(2,1) * 3600.0 ELSE WRITE(CMINOR,1052) R(2,1) * 3600.0 END IF CPA = ' ' EMAJOR = ' ' EMINOR = ' ' EPA = ' ' ELSE C Some upper limits C Major axis IF (R(1,2).GT.0.0) THEN WRITE(CMAJOR,1050) R(1,1) * 3600.0 WRITE(EMAJOR,1050) (R(1,2) + R(1,3)) * 0.5 * 3600.0 RESOLV = .TRUE. HAFRES(2) = .TRUE. ELSE IF ((R(1,1)*3600.0).LE.99.0) THEN WRITE(CMAJOR,1051) R(1,1) * 3600.0 ELSE WRITE(CMAJOR,1052) R(1,1) * 3600.0 END IF EMAJOR = ' ' END IF C Minor axis IF (R(2,2).GT.0.0) THEN WRITE(CMINOR,1050) R(2,1) * 3600.0 WRITE(EMINOR,1050) (R(2,2) + R(2,3)) * 0.5 * 3600.0 RESOLV = .TRUE. HAFRES(3) = .TRUE. ELSE IF ((R(2,1)*3600.0).LE.99.0) THEN WRITE(CMINOR,1051) R(2,1) * 3600.0 ELSE WRITE(CMINOR,1052) R(2,1) * 3600.0 END IF EMINOR = ' ' END IF C Position angle IF (RESOLV) THEN WRITE(CPA,1050) R(3,1) ERRPA = MIN (90.0, 0.5 * ABS (R(3,2) + R(3,3))) WRITE(EPA,1050) ERRPA ELSE CPA = ' ' EPA = ' ' END IF END IF C One axis resolved? HAFRES(1) = HAFRES(2) .OR. HAFRES(3) C RETURN C----------------------------------------------------------------------- 1050 FORMAT (F6.1) 1051 FORMAT (' <', F4.1) 1052 FORMAT (' <', F4.0) END SUBROUTINE NEWPS (DECRAD, SPMJY, PSMAJ, PSMIN, PSPA) C----------------------------------------------------------------------- C J. J. Condon 2001 Feb 9 C To correct the NVSS "Wall problem", the FORTRAN subroutine newps C calculates the NVSS point-source response in terms of an C "effective" elliptical Gaussian with C psmaj = FWHM major axis (arcsec) C psmin = FWHM minor axis (arcsec) C pspa = major axis position angle (deg) C east of north C from the source parameters C decrad = J2000 declination (radians) C spmjy = raw peak flux density (mJy) C C----------------------------------------------------------------------- REAL DECRAD, SPMJY, PSMAJ, PSMIN, PSPA C REAL VLALAT, ZARAD, OMEGA C VLA latitude, radians DATA VLALAT /0.5948/ C----------------------------------------------------------------------- C first calculate dirty-beam parameters for the C 1 mJy/beam uncleaned "pedestal". C DnC configuration: IF (DECRAD .LE. -0.176719 .OR. DECRAD .GT. 1.361255) THEN C minor axis = 51 arcsec, in units of 45 arcsec C restoring beam PSMIN = 1.133333 C zenith angle in radians ZARAD = DECRAD - VLALAT C dirty beam solid angle from NVSS raw fits C (/aten/ATEN_1/isot.dir/beamsize.sm, BEAMSIZE.OUT) OMEGA = 1.405 + 0.065 / COS (ZARAD) PSMAJ = OMEGA / PSMIN PSPA = 0. ELSE C D configuration: C minor axis = 54 arcsec, in units of 45 arcsec C restoring beam PSMIN = 1.2000 ZARAD = DECRAD - VLALAT C beam is elongated north-south near transit PSPA = 0. C "zenith" band 6, observed off transit: IF (DECRAD .GT. 0.437220 .AND. DECRAD .LE. 0.768818) THEN C average zenith angle is 27 deg ZARAD = 0.471 C beam is elongated east-west instead of north-south PSPA = 90. END IF OMEGA = 0.91 + 0.61 / COS (ZARAD) PSMAJ = OMEGA / PSMIN END IF C next calculate point-source response C for whole source, cleaned down to 1 mJy/beam C equation derived in notes 010201 R = SPMJY - 1. PSMIN = (R + PSMIN*PSMIN*PSMIN) / (R + PSMIN) PSMIN = SQRT (PSMIN) * 45.0 PSMAJ = (R + PSMAJ*PSMAJ*PSMAJ) / (R + PSMAJ) PSMAJ = SQRT (PSMAJ) * 45.0 RETURN END SUBROUTINE DECONV (PTSMAJ, PTSMIN, PTSPA, * FITMAJ, UFITMJ, FITMIN, UFITMN, FITPA, * SRCMAJ, USRCMJ, SRCMIN, USRCMN, SRCPA) C----------------------------------------------------------------------- C Thus subroutine deconvolves a gaussian point-source response from a C fitted elliptical gaussian to yield the gaussian source parameters. C This calculation is from AJ, 109, 2318 and my notes of 010209. C DECONV also determines whether each source axis is significantly C (98% confidence = 2.33 sigma) resolved. If so, the rms error C in that axis is calculated. If not, the 98% confidence upper limit C to that axis is calculated and its rms error is set to -1. C These error calculations are based on the old DECONV subroutine C and the NVSS paper (AJ, 115, 1693). C----------------------------------------------------------------------- C INPUT ARGUMENTS: C ptsmaj = major-axis of pt src response (arcsec) C ptsmin = minor-axis of pt src response (arcsec) C ptspa = position angle of pt src response (DEG) C MUST BE IN RANGE -180 DEG TO +180 DEG C fitmaj = major-axis of raw fit (arcsec) C ufitmaj = rms error in fitted major axis (arcsec) C fitmin = minor-axis of raw fit (arcsec) C ufitmin = rms error in fitted minor axis (arcsec) C fitpa = position angle of raw fit (DEG) C MUST BE IN RANGE -180 DEG TO +180 DEG C OUTPUT ARGUMENTS: C srcmaj = major-axis of deconv source (arcsec) C or 98% confidence upper limit C usrcmaj = rms error in srcmaj if resolved, C or -1 if srcmaj is an upper limit C srcmin = minor-axis of deconv source (arcsec) C or 98% confidence upper limit C usrcmin = rms error in srcmin if resolved, C or -1 if srcmaj in an upper limit C srcpa = position angle of deconv source (deg) C RETURNED IN RANGE -90 TO +90 deg C----------------------------------------------------------------------- REAL PTSMAJ, PTSMIN, PTSPA, FITMAJ, UFITMJ, FITMIN, UFITMN, * FITPA, SRCMAJ, USRCMJ, SRCMIN, USRCMN, SRCPA C REAL PI, DEGRAD, TEMP, SIGMAX, PTPRAD, FITPRD, PHIRAD, * PTMJSQ, PTMNSQ, FTMJSQ, FTMNSQ, DELTAR, SCMJSQ, * SCMNSQ, PHIRD2, CS2PHI, SN2PHI, DENOM, XNUM, * COSPHI, SINPHI, SRCPAR, PTSPAR, RADARG, SCMJMX, * UMSMIN, SCMNMX, TEST C confidence-level parameters for "significant" C resolution; 2.33 => 98% confidence PARAMETER (SIGMAX = 2.33) DATA PI /3.14159265/ C----------------------------------------------------------------------- DEGRAD = 180. / PI C fix any reversed input major and minor axes IF (PTSMAJ .LT. PTSMIN) THEN TEMP = PTSMAJ PTSMAJ = PTSMIN PTSMIN = TEMP PTSPA = PTSPA + 90. IF (PTSPA .GT. 180.) PTSPA = PTSPA - 360. END IF IF (FITMAJ .LT. FITMIN) THEN TEMP = FITMAJ FITMAJ = FITMIN FITMIN = TEMP FITPA = FITPA + 90. IF (FITPA .GT. 180.) FITPA = FITPA - 360. END IF C convert PA'S to radians, in range 0 to pi PTPRAD = PTSPA / DEGRAD IF (PTPRAD .LT. 0.) PTPRAD = PTPRAD + PI FITPRD = FITPA / DEGRAD IF (FITPRD .LT. 0.) FITPRD = FITPRD + PI C calculate PA difference phirad (radians) C between raw fit and point-source response PHIRAD = FITPRD - PTPRAD C make range of phirad from -pi/2 to +pi/2 IF (PHIRAD .GT. PI/2.) PHIRAD = PHIRAD - PI IF (PHIRAD .LT. -PI/2.) PHIRAD = PHIRAD + PI C calculate squares of fwhm sizes PTMJSQ = PTSMAJ * PTSMAJ PTMNSQ = PTSMIN * PTSMIN FTMJSQ = FITMAJ * FITMAJ FTMNSQ = FITMIN * FITMIN C C do various special cases for which the C general formula is not needed or is C not well defined (division by zero) C C ptsmajsq = ptsminsq? (circular pt-src response) IF (ABS (PTMJSQ - PTMNSQ) .LT. 1.E-3) THEN C delta is the position angle difference C between the deconvolved src major axis and C the pt src response major axis DELTAR = PHIRARAD SCMJSQ = FTMJSQ - PTMJSQ SCMNSQ = FTMNSQ - PTMNSQ GO TO 90 END IF C fitmajsq = fitminsq? (circular raw fit) IF (ABS (FTMJSQ - FTMNSQ) .LT. 1.E-3) THEN DELTAR = PI / 2. SCMJSQ = FTMJSQ - PTMNSQ SCMNSQ = FTMNSQ - PTMJSQ GO TO 90 END IF C phirad = 0? (fit parallel to pt-src response) IF (ABS (PHIRAD) .LT. 1.E-5) THEN SCMJSQ = FTMJSQ - PTMJSQ SCMNSQ = FTMNSQ - PTMNSQ IF (SCMJSQ .GE. SCMNSQ) THEN DELTAR = 0. ELSE TEMP = SCMJSQ SCMJSQ = SCMNSQ SCMNSQ = TEMP DELTAR = PI / 2. END IF GO TO 90 END IF C phirad = +- pi/2? (fit perpendicular to beam) IF (ABS (PHIRAD - PI/2.) .LT. 1.E-5 .OR. * ABS (PHIRAD + PI/2.) .LT. 1.E-5) THEN SCMJSQ = FTMJSQ - PTMNSQ SCMNSQ = FTMNSQ - PTMJSQ DELTAR = PHIRAD GO TO 90 END IF C C end of special cases C C calculate deltarad PHIRD2 = 2. * PHIRAD CS2PHI = COS (PHIRD2) SN2PHI = SIN (PHIRD2) DENOM = (FTMJSQ - FTMNSQ) * CS2PHI * - (PTMJSQ - PTMNSQ) XNUM = (FTMJSQ - FTMNSQ) * SN2PHI C calculate deltarad IF (ABS(DENOM) .GE. 1.E-5) THEN DELTAR = 0.5 * ATAN (XNUM / DENOM) ELSE IF (XNUM .GT. 0.) DELTAR = PI / 4. IF (XNUM .LE. 0.) DELTAR = -PI / 4. END IF C range is now +- pi/4. C resolve ambiguities to make range +- pi/2. IF (DENOM .LT. 0.) THEN IF (XNUM .GT. 0.) DELTAR = DELTAR + PI / 2. IF (XNUM .LE. 0.) DELTAR = DELTAR - PI / 2. END IF C calculate srcmajsq COSPHI = COS (PHIRAD) SINPHI = SIN (PHIRAD) SCMJSQ = (FTMJSQ - FTMNSQ) * COSPHI * SINPHI SCMJSQ = SCMJSQ / (COS (DELTAR) * SIN (DELTAR)) SCMJSQ = 0.5 * (SCMJSQ + FTMJSQ + FTMNSQ * - (PTMJSQ + PTMNSQ)) C calculate srcminsq SCMNSQ = FTMJSQ + FTMNSQ - (PTMJSQ + PTMNSQ) * - SCMJSQ 90 CONTINUE C srcmajsq < srcminsq? IF (SCMJSQ .LT. SCMNSQ) THEN TEMP = SCMJSQ SCMJSQ = SCMNSQ SCMNSQ = TEMP DELTAR = DELTAR + PI / 2. END IF C if deconvolution fails (that is, the C square of the source size is negative, C set source size negative IF (SCMJSQ .GT. 0.) THEN SRCMAJ = SQRT (SCMJSQ) ELSE SRCMAJ = -SQRT(-SCMJSQ) END IF IF (SCMNSQ .GT. 0.) THEN SRCMIN = SQRT (SCMNSQ) ELSE SRCMIN = -SQRT(-SCMNSQ) END IF C calculate source position angle SRCPAR = DELTAR + PTPRAD SRCPA = SRCPAR * DEGRAD C make sure srcpa in range -90 to +90 deg IF (SRCPA .LT. -90.) SRCPA = SRCPA + 180. IF (SRCPA .GE. 90.) SRCPA = SRCPA - 180. C C end of fit calculation C next calculate fit uncertainties C C test for significant major-axis resolution C (see section 5.2.4 in AJ, 115, 1693) PTSPAR = SQRT (PTSMAJ*PTSMAJ*COSPHI*COSPHI + * PTSMIN*PTSMIN*SINPHI*SINPHI) IF (FITMAJ .LT. PTSPAR) THEN TEMP = PTSPAR ELSE TEMP = FITMAJ END IF C is fit major axis > 2.33*sigma + beam C projected onto fit major axis? TEST = TEMP - UFITMJ * SIGMAX - PTSPAR IF (TEST .GT. 0.) THEN C src major axis is significantly resolved C calculate rms error in src major axis C NVSS paper, eq. 32a RADARG = (TEMP - UFITMJ) * (TEMP - UFITMJ) - * PTSPAR * PTSPAR UMSMAJ = SRCMAJ - SQRT(RADARG) RADARG = (TEMP + UFITMJ) * (TEMP + UFITMJ) - * PTSPAR * PTSPAR UPSMAJ = SQRT(RADARG) - SRCMAJ USRCMJ = 0.5 * (UMSMAJ + UPSMAJ) ELSE C src major axis is not significantly resolved C calculate source major axis C 98% confidence upper limit SCMJMX = TEMP + SIGMAX * UFITMJ RADARG = SCMJMX * SCMJMX - PTSPAR * PTSPAR SCMJMX = SQRT (RADARG) C Set source size limit to srcmajmax SRCMAJ = SCMJMX USRCMJ = -1. END IF C test for significant minor-axis resolution PTSPAR = SQRT (PTSMAJ*PTSMAJ*SINPHI*SINPHI + * PTSMIN*PTSMIN*COSPHI*COSPHI) IF (FITMIN .LT. PTSPAR) THEN TEMP = PTSPAR ELSE TEMP = FITMIN END IF TEST = TEMP - UFITMN * SIGMAX - PTSPAR IF (TEST .GT. 0.) THEN C src minor axis is significantly resolved C calculate rms error in src minor axis RADARG = (TEMP - UFITMN) * (TEMP - UFITMN) - * PTSPAR * PTSPAR UMSMIN = SRCMIN - SQRT(RADARG) RADARG = (TEMP + UFITMN) * (TEMP + UFITMN) - * PTSPAR * PTSPAR UPSMIN = SQRT(RADARG) - SRCMIN USRCMN = 0.5 *(UPSMIN + UMSMIN) ELSE C src minor axis is not significantly resolved C calculate source minor axis upper limt SCMNMX = TEMP + SIGMAX * UFITMN RADARG = SCMNMX * SCMNMX - PTSPAR * PTSPAR SCMNMX = SQRT (RADARG) C set source size limit to srcminmax SRCMIN = SCMNMX USRCMN = -1. END IF C make sure larger upper limit is called srcmaj IF (USRCMN .EQ. -1. .AND. USRCMJ .EQ. -1.) THEN IF (SRCMIN .GT. SRCMAJ) THEN TEMP = SRCMIN SRCMIN = SRCMAJ SRCMAJ = TEMP END IF END IF RETURN END