C----------------------------------------------------------------------- C; Copyright (C) 2002 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, * IRMS, BEAM, FITTED, DORAW, FBLANK, FLUX, EFLUX, ERRRA, ERRDEC, * CMAJOR, CMINOR, CPA, EMAJOR, EMINOR, EPA) C----------------------------------------------------------------------- C Routine to apply any corrections to 4MASS fitted source parameters C and determine errors. Unless FITTED, the model parameters returned C are deconvolved using the restoring BEAM. 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 IRMS R RMS (sigma) in Ipol. 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 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, IRMS, BEAM(3), FBLANK, FLUX, * EFLUX LOGICAL FITTED, DORAW CHARACTER 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, PEAKX, PERR, * DECRAD, SPMJY, PSF(3), PEAKP INTEGER IER LOGICAL RESOLV, HAFRES(3) C Calibration component of C position error (squared) REAL CALERA, CALEDE PARAMETER (CALERA = (0.0) ** 2) PARAMETER (CALEDE = (0.0) ** 2) C Amplitude calibration C uncertainty REAL CALAER PARAMETER (CALAER = 0.03) C Position bias DOUBLE PRECISION BIASRA, BIASDE PARAMETER (BIASRA = 0.0D0) PARAMETER (BIASDE = 0.0D0) C Degrees per radian DOUBLE PRECISION DG2RAD PARAMETER (DG2RAD = 1.745329252E-2) C----------------------------------------------------------------------- C Convert sizes to degrees PSF(1) = BEAM(1) PSF(2) = BEAM(2) 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(2)/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(2)/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(2)/MINOR)**2)**2.5) * * SNR * SNR C Flux errors ERRPEK = SQRT ((2.0 * PEAKX*PEAKX / SNRAMP)) C Add Flux calibration (3%) error ERRPEK = SQRT (ERRPEK*ERRPEK + (CALAER * PEAKX)**2) C 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 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(2))**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*3600.0, MINOR*3600.0, POSANG, * ERRMAJ*3600.0, ERRMIN*3600.0, ERRPA, * PSF(1)*3600.0, PSF(2)*3600.0, 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)) 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))) 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)) END IF END IF C Convert units for output ERRRA = ERRRA * 3600.0 ERRDEC = ERRDEC * 3600.0 C debug C WRITE(0,1234) BEMRAT, FLUX, PEAK, TFLUX C 1234 FORMAT('BEMRAT,FLUX,PEAK,TFLUX =',1PE12.5,3E12.5) C 999 RETURN C----------------------------------------------------------------------- 1050 FORMAT (F6.1) 1051 FORMAT ('<', F5.1) 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 Inputs: C BMAJ R Fitted major axis (asec) C BMIN R Fitted minor axis (asec) C BPA R Fitted pos. angle (deg) C BMAJE R Fitted major axis error (asec) C BMINE R Fitted minor axis error (asec) C BPAE R Fitted pos. angle error (deg) C CBMAJ R Clean beam major axis (asec) C CBMIN R Clean beam minor axis (asec) 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, * 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,3) = 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) WRITE(CMINOR,1050) R(2,1) WRITE(CPA,1050) R(3,1) WRITE(EMAJOR,1050) (R(1,2) + R(1,3)) * 0.5 WRITE(EMINOR,1050) (R(2,2) + R(2,3)) * 0.5 ERRPA = MIN (90.0, 0.5 * ABS (R(3,2) + R(3,3))) WRITE(EPA,1050) ERRPA RESOLV = .TRUE. ELSE IF ((R(1,2).LE.0.0) .AND. (R(2,2).LE.0.0)) THEN C Upper limits only IF (R(1,1).LE.99.0) THEN WRITE(CMAJOR,1051) R(1,1) ELSE WRITE(CMAJOR,1052) R(1,1) END IF IF (R(2,1).LE.99.0) THEN WRITE(CMINOR,1051) R(2,1) ELSE WRITE(CMINOR,1052) R(2,1) 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) WRITE(EMAJOR,1050) (R(1,2) + R(1,3)) * 0.5 RESOLV = .TRUE. HAFRES(2) = .TRUE. ELSE IF (R(1,1).LE.99.0) THEN WRITE(CMAJOR,1051) R(1,1) ELSE WRITE(CMAJOR,1052) R(1,1) END IF EMAJOR = ' ' END IF C Minor axis IF (R(2,2).GT.0.0) THEN WRITE(CMINOR,1050) R(2,1) WRITE(EMINOR,1050) (R(2,2) + R(2,3)) * 0.5 RESOLV = .TRUE. HAFRES(3) = .TRUE. ELSE IF (R(2,1).LE.99.0) THEN WRITE(CMINOR,1051) R(2,1) ELSE WRITE(CMINOR,1052) R(2,1) 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 (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, UMSMAJ, UPSMAJ, UPSMIN 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 COSPHI = COS (PHIRAD) SINPHI = SIN (PHIRAD) 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 = PHIRAD 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 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