C----------------------------------------------------------------------- C; Copyright (C) 1995,1996,1997,1998 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 Inputs: C RA D RA (deg) C DEC D Dec (deg) C PEAK R Peak Ipol 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) 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 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 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/BEAM(1)) * (MINOR/BEAM(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*BEAM(1)*BEAM(1))) * * ((1.0+(BEAM(1)/MAJOR)**2)**1.5) * * ((1.0+(BEAM(1)/MINOR)**2)**1.5) * * SNR * SNR C SNR**2 for major axis error SNRMAJ = (MAJOR*MINOR/(4.0*BEAM(1)*BEAM(1))) * * ((1.0+(BEAM(1)/MAJOR)**2)**2.5) * * ((1.0+(BEAM(1)/MINOR)**2)**0.5) * * SNR * SNR C SNR**2 for minor axis/PA errors SNRMIN = (MAJOR*MINOR/(4.0*BEAM(1)*BEAM(1))) * * ((1.0+(BEAM(1)/MAJOR)**2)**0.5) * * ((1.0+(BEAM(1)/MINOR)**2)**2.5) * * SNR * SNR C Bias correct peak flux Cold IF (.NOT.DORAW) PEAKX = PEAKX + BIASAV - 1.0E-3*(BIASFA / PEAKX) 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 <2% 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, BEAM(2)) TMIN = MAX (MINOR, BEAM(2)) C Axis sizes include 2% C calibration error. ERRMAJ = SQRT (((2.0 * TMAJ**2) / SNRMAJ) + (0.02*BEAM(1))**2) ERRMIN = SQRT (((2.0 * TMIN**2) / SNRMIN) + (0.02*BEAM(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, BEAM(1), BEAM(2), BEAM(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/BEAM(2))) TFLUX = PEAKP * MAX (1.0, MAJOR/BEAM(1)) ELSE C Minor axis only PEAKP = PEAKX * SQRT (MAX (1.0, MAJOR/BEAM(1))) TFLUX = PEAKP * MAX (1.0, MINOR/BEAM(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)) + * (BEAM(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 C----------------------------------------------------------------------- C Get deconvolved sizes and C errors. CALL DECONV (CBMIN, BMAJ, BMAJE, BMIN, BMINE, R(1,1), R(1,2), * R(1,3), R(2,1), R(2,2), R(2,3)) R(3,1) = BPA 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 DECONV (BWID, FMAJ, UFMAJ, FMIN, UFMIN, SMAJ, UPSMAJ, * UMSMAJ, SMIN, UPSMIN, UMSMIN) C----------------------------------------------------------------------- C This subroutine deconvolves a fitted elliptical gaussian into an C elliptical gaussian source and a CIRCULAR gaussian beam, and it C calculates the errors in the deconvolved sizes. The arguments are C dimensionless angles, so any single unit (e.g., deg or arcsec) may C be used for all of them. The fitted major and minor axes are C compared with the beamwidth to see if the source is well resolved, C that is, the fitted axis is some number of sigma bigger than the C beam (e.g., SIGMIN = 2.33 => 98% confidence that the source is C resolved). If a source axis is well resolved, both its +rms error C and -rms are calculated. The +/-rms error is the difference between C the deconvolved size and the deconvolved size corresponding to the C fitted size +/- its rms error. Unless the source is just barely C resolved, the +rms and -rms errors are nearly equal and can be C averaged. If either source axis is not well resolved, then that axis C (either SMAJ or SMIN) is set to an upper limit whose significance is C determined by SIGMAX (e.g., SIGMAX = 2.33 => 98% confidence upper C limit). If the fitted size is significantly smaller than the beam C (an "impossble", or at least "improbable", source), the source axis C size is set to -1. The errors associated with unresolved source axes C (e.g., UMSMAJ, UMSMIN, UPSMAJ, UPSMIN) are set to -1. to indicate C that the quoted sizes are upper limits. C Version: 95/05/13 J. J. Condon C 97/12/10 Corrected error in calculation of C rms uncertainties in deconvolved sizes J. J. Condon C Inputs: C BWID R Beam circular gaussian FWHM C FMAJ R Fitted elliptical gaussian FWHM major axis C FMIN R Fitted elliptical gaussian FWHM minor axis C UFMAJ R RMS error in FMAJ C UFMIN R RMS error in FMIN C Outputs: C SMAJ R Source elliptical gaussian FWHM major axis C SMIN R Source elliptical gaussian FWHM minor axis C UMSMAJ R -RMS error in SMAJ C UMSMIN R -RMS error in SMIN C UPSMAJ R +RMS error in SMAJ C UPSMIN R +RMS error in SMIN C Control parameters: C SIGMAX R Number of UFMAJ, UFMIN to add to set source size limit C SIGMIN R Number of UFMAJ, UFMIN to decide if source well C resolved C (for both sigmax and sigmin, 2.33 => 98% confidence --- C see Bevington, P. R., Data Reduction and Error Analysis for the C Physical Sciences, Table C-2) C----------------------------------------------------------------------- REAL BWID, FMAJ, FMIN, UFMAJ, UFMIN, SMAJ, SMIN, UMSMAJ, * UMSMIN, UPSMAJ, UPSMIN C REAL SIGMAX, SIGMIN REAL TEST, RADARG C Confidence-level parameters C sigmax, sigmin: 1.28 => 80%, C 1.65 => 90%, 2.33 => 98% C Both must be > 1.00 for C equations to be valid. PARAMETER (SIGMAX = 2.33) PARAMETER (SIGMIN = 2.33) C----------------------------------------------------------------------- C Calculate source FWHM major axis SMAJ = FMAJ * FMAJ - BWID * BWID IF (SMAJ .GT. 0.) THEN SMAJ = SQRT (SMAJ) ELSE SMAJ = 0. END IF C Calculate source FWHM minor C axis SMIN = FMIN * FMIN - BWID * BWID IF (SMIN .GT. 0.) THEN SMIN = SQRT (SMIN) ELSE SMIN = 0. END IF C Test for significant C major-axis resolution TEST = (FMAJ - UFMAJ * SIGMIN) TEST = TEST * TEST - BWID * BWID IF (TEST .GT. 0.) THEN C Source major axis is well C resolved. Calculate source C major axis - rms error. RADARG = (FMAJ - UFMAJ) * (FMAJ - UFMAJ) - BWID * BWID UMSMAJ = SMAJ - SQRT (RADARG) C Calculate source major axis C +rms error. RADARG = (FMAJ + UFMAJ) * (FMAJ + UFMAJ) - BWID * BWID UPSMAJ = SQRT(RADARG) - SMAJ ELSE C Source major axis is not C significantly resolved. C Calculate source major axis C upper limit. SMAJ = FMAJ + SIGMAX * UFMAJ IF (SMAJ .GT. BWID) THEN RADARG = SMAJ * SMAJ - BWID * BWID SMAJ = SQRT (RADARG) ELSE SMAJ = -1. END IF C Set source uncertainties to C -1 UPSMAJ = -1. UMSMAJ = -1. END IF C Test for significant C minor-axis resolution TEST = (FMIN - UFMIN * SIGMIN) TEST = TEST * TEST - BWID * BWID IF (TEST .GT. 0.) THEN C Source minor axis is well C resolved. Calculate source C minor axis -rms error. RADARG = (FMIN - UFMIN) * (FMIN - UFMIN) - BWID * BWID UMSMIN = SMIN - SQRT (RADARG) C Calculate source minor axis C +rms error. RADARG = (FMIN + UFMIN) * (FMIN + UFMIN) - BWID * BWID UPSMIN = SQRT(RADARG) - SMIN ELSE C Source minor axis is not C significantly resolved. C Calculate source minor axis C upper limit SMIN = FMIN + SIGMAX * UFMIN IF (SMIN .GT. BWID) THEN RADARG = SMIN * SMIN - BWID * BWID SMIN = SQRT(RADARG) ELSE SMIN = -1. END IF C Set source uncertainties to C -1 UPSMIN = -1. UMSMIN = -1. END IF C 999 RETURN END