>
> |
%META:TOPICINFO{author="FrankGhigo" date="1214934540" format="1.0" version="1.1"}%
C 2004 Aug 19 version J. J. Condon
c calfind2.f finds suitable GBT calibrators from file PCALS3.3
C 2005 Aug 19 version: changed from PCALS3.3 to PCALS4.0
C and prints them on the screen.
C 2007 Jul 25 version: changed from PCALS4.0 to PCALS4.1
C and recognizes calcode 'P' for PTCS "gold" calibrators
C usable for Ka-band pointing calibration.
C 2008 Jun 3 version: limits slew time and zenith angle
C as well as angular separation
C and relaxes criteria if needed to find a calibrator.
C
character*1 decsgn, gold
write (*, *) 'CALFIND: I find suitable GBT pointing calibrators no
*rth of -40 deg dec.'
write (*, *)
C observing frequency in GHz, needed to calculate GBT FWHM beamwidth
write (*, *) 'Enter the observing frequency (GHz)'
read (*, *) fghz
C enter the calibrator search radius (deg)
write (*, *)
write (*, *) 'Enter the calibration-source search radius (deg)'
read (*, *) searchdeg
C enter the calibrator minimum flux density
sjymin = 0.
write (*, *)
write (*, *) 'Enter the minimum calibrator flux density (Jy) or ze
*ro for no flux limit'
read (*, *) sjymin
C return only PTCS "gold" calibrators? (Y/N)
C (These are for PTCS; most users don't need these and should answer 'N')
write (*, *)
write (*, *) 'Return only PTCS "gold" pointing calibrators? (Y/N)'
write (*, *) '(Most observers should answer "N")'
read (*, *) gold
if (gold .eq. 'Y' .or. gold .eq. 'y') gold = 'Y'
C enter search center position
write (*, *)
write (*, *) 'Enter the J2000 search center in the columns below:'
write (*, *) 'HH MM SS.SS -DD MM SS.S'
read (*, 1030) ihra, imra, sra, decsgn, idd, imd, sd
C enter Local Sidereal Time (LST) in decimal hours
C (needed for determining the azimuth and elevation of the calibrator)
C If you just want a list of all suitable calibrators within the search
C area, enter -1)
write (*, *)
write (*, *) 'Enter LST (decimal hours), 0. < LST < 24.',
* ' to limit slew time and tracking error.'
write (*, *) 'Note: Enter -1. to get all nearby calibrators',
* ' regardless of azimuth or elevation'
read (*, *) hlst
write (*, *)
C call the CALFIND subroutine
call CALFIND (fghz, searchdeg, sjymin, ihra, imra, sra,
* decsgn, idd, imd, sd, gold, hlst)
stop
C format statements
1030 format (1x, i2.2, 1x, i2.2, 1x, f5.2, 1x,
* a1, i2.2, 1x, i2.2, 1x, f4.1)
1040 format ('LST = ', f6.1, ' hours?')
end
C
C ----------------------------------------------------------------------
subroutine CALFIND (fghz, searchdeg, sjymin, ihra, imra, sra,
* decsgn, idd, imd, sd, gold, hlst)
C ----------------------------------------------------------------------
C
C 2004 Aug 19 version J. J. Condon
C (Bug found by Dana Balser fixed)
C (Estimated flux densities at observing frequency printed)
C
C 2008 Jun 3 version J. J. Condon
C (An option to select only those calibrators reachable with short slew
C times was added, to avoid long azimuth slews near the zenith.
C This requires a new calling argument, hlst, the LST in decimal hours.
C If hlst .lt. 0. then slew times and calibrator elevations are ignored.
C If hlst .ge. 0. then calfind rejects calibrators with slew times
C exceeding the larger of
C 1 minute or the 0.1 minutes * search radius (in degrees).
C This avoids long slews near the zenith.
C Also, if hlst .gt. 0. and the calibrator elevation angle is too
C high (> 88 deg) for tracking, the calibrator is rejected.
C If no suitable calibrator is found, the search area and flux
C criteria are relaxed until a suggested calibrator is found
C and a warning message is printed.)
C
C subroutine CALFIND finds suitable GBT pointing calibration sources
C in the disk file PCALS4.1 that satisfy user-specified constraints on
C observing frequency, flux density, search radius, and center position.
C
C calling arguments:
C fghz = observing frequency (GHz)
C searchdeg = search radius (degrees)
C sjymin = minimum calibrator flux density (Jy) at frequency fghz
C ihra = search center hours of right ascension
C imra = search center minutes of right ascension
C sra = search center seconds of right ascension
C decsgn = search center sign of declination ('+', '-', or ' ')
C idd = search center degrees of declination
C imd = search center arcminutes of declination
C sd = search center arcseconds of declination
C gold = 'Y' for PTCS "gold" calibrators only
C hlst = LST in decimal hours, needed to calculate slew times
C if hlst .lt. 0., slew times are ignored
C
C calling argument and local variable types:
character*1 decsgn, col(90), decsgnc, gold, goldmark
character*9 iaunam
integer*4 ihra, imra, idd, imd
real*4 fghz, searchdeg, sjyest, sjymin, sra, sd
C
real*4 pi/3.1415926535/, sjy(8), freq(7)
C GBT telescope latitude (radians), needed for slew time calculation
tellat = 0.6681
C list of frequencies (GHz) at which we have flux densities
C this may have to be expanded when the 30 GHz receiver goes up
data freq /1.4, 4.9, 8.6, 15., 22., 43., 86./
C goldmark ="P" in PCALS4.1 for PTCS "Gold" pointing calibrators
goldmark = ' '
C
C input data file of calibration sources
open (unit = 21, file = 'PCALS4.1', status = 'old')
twopi = 2. * pi
piby2 = pi / 2.
write (*, *) '----------------------------------------------------
*----'
C calculate the true GBT FWHM beamwidth in arcsec
fwhmtrue = 740. / fghz
fwhm = fwhmtrue
C calculate wavelength in cm
cmlambda = 30. / fghz
if (fwhm .gt. 540.) then
write (* ,*) 'CALFIND: Warning! Large GBT beam; use only strong
* '
write (*, *) ' calibrators to minimize confusion'
C estimate rms confusion on the GBT (see Eq. 10 of
C "Single Dish Radio Astronomy", ASP Conference Series, vol. 278,
C page 155)
rmsjy = 0.030 * fghz**(-2.7)
write (*, 1024) rmsjy
C reset larger true beam to 539. arcsec to avoid excluding
C all calibrators at low frequencies
fwhm = 539.
end if
C
C choose the nearest band (20, 6, 3.5, 2, 1.3, 0.7, or 0.3 cm) for
C estimating flux at fghz
iband = 8
if (fghz .lt. 86.) iband = 7
if (fghz .lt. 43.) iband = 6
if (fghz .lt. 22.) iband = 5
if (fghz .lt. 15.) iband = 4
if (fghz .lt. 8.5) iband = 3
if (fghz .lt. 4.9) iband = 2
if (fghz .lt. 1.4) iband = 1
C
C reset the input search radius to 180 deg
C if it was larger than 180 deg. This will search the whole sky.
if (searchdeg .gt. 180.) searchdeg = 180.
C convert the search radius from degrees to radians
searchrad = searchdeg * pi / 180.
C convert the search radius from radians to
C cartesian offset on unit sphere
offsetmax = 2. * sin (searchrad / 2.)
C
C check for a bad input search-center position
if (ihra .gt. 23 .or. imra .gt. 59 .or. sra .gt. 60.0 .or.
* ihra .lt. 0 .or. imra .lt. 0 .or. sra .lt. 0.0 .or.
* idd .gt. 90 .or. imd .gt. 59 .or. sd .gt. 60.0 .or.
* idd .lt. -90 .or. imd .lt. 0 .or. sd .lt. 0.0) then
write (*, *) 'Bad input position:'
write (*, 1030) ihra, imra, sra, decsgn, idd, imd, sd
return
end if
if (decsgn .ne. ' ' .and. decsgn .ne. '+' .and.
* decsgn .ne. '-') then
write (*, *) 'Bad input position:'
write (*, 1030) ihra, imra, sra, decsgn, idd, imd, sd
return
end if
C
C convert the search center to radians of ra, dec
racenter = 3600 * ihra + 60 * imra
racenter = (racenter + sra) * pi / 43200.
deccenter = 3600 * idd + 60 * imd + sd
deccenter = deccenter * pi / 648000.
if (decsgn .eq. '-') deccenter = -deccenter
C
if (hlst .ge. 0.) then
C calculate LST in radians
radlst = hlst * pi / 12.
C calculate the hour angle (radians) of the search center
hacenter = radlst - racenter
C calculate the current elevation (radians) of the search center
elcenter = asin (sin (tellat) * sin (deccenter) +
* cos (tellat) * cos (deccenter) * cos (hacenter) )
C calculate the current azimuth (radians) of the search center
denom = cos (tellat) * sin (deccenter) -
* sin (tellat) * cos (deccenter) * cos (hacenter)
xnum = -sin(hacenter) * cos (deccenter)
azcenter = atan2 (xnum,denom)
C convert azimuth range from (-pi to +pi) to (0 to 2pi)
if (azcenter .lt. 0.) azcenter = azcenter + 2. * pi
C (see http://aa.usno.navy.mil/faq/docs/Alt_Az.php)
end if
C
C convert the search center to cartesian coordinates on a unit sphere
C to avoid polar coordinate poles and ra wrap ambiguities
x = cos (racenter) * cos (deccenter)
y = sin (racenter) * cos (deccenter)
z = sin (deccenter)
C
C search for suitable calibrators near the search center
C
iwarning = 0
5 continue
C Write warning that search failed and criteria had to be relaxed
if (iwarning .eq. 1) then
write (*, *) 'WARNING: No calibrator met all of your search',
*' criteria.'
write (*, *) ' Use these suggested calibrator(s) or'
write (*, *) ' retry using weaker criteria.'
end if
C set the maximum slew times in minutes
tslewmax = searchdeg / 10.
if (tslewmax .lt. 1.) tslewmax = 1.
C how many calibrators have been found?
ncals = 0
C skip the two header lines of the input file
read (21, 1160) col
read (21, 1160) col
C
C
C loop over the calibrator list in the input data file
do 100 inline = 1, 100000
C
C read the data for one calibrator source
read (21, 1110, end = 200) iaunam, ihrac, imrac, srac,
* decsgnc, iddc, imdc, sdc, min, max, (sjy(i), i = 1, 6),
* goldmark
sjy(7) = 0.
sjy(8) = 0.
C
C convert the calibrator position to radians of arc
racal = 3600 * ihrac + 60 * imrac
racal = (racal + srac) * pi / 43200.
deccal = iddc * 3600 + imdc * 60
deccal = (deccal + sdc) * pi / 648000.
if (decsgnc .eq. '-') deccal = -deccal
C convert the calibrator position to cartesian coordinates
xcal = cos (racal) * cos (deccal)
ycal = sin (racal) * cos (deccal)
zcal = sin (deccal)
C calculate the calibrator offset from the input position
C in cartesian coordinates
offset = (xcal - x)*(xcal - x) + (ycal - y)*(ycal - y) +
* (zcal - z)*(zcal -z)
offset = sqrt(offset)
C
C Is this calibrator too far away?
if (offset .gt. offsetmax) go to 100
C Calibrator is inside search circle. Convert its offset to degrees.
offsetrad = 2. * asin (offset / 2.)
offsetdeg = offsetrad * 180. / pi
C
C Is the GBT beam too large for this calibrator?
xmax = max
C
C If no maximum beam size is specified (max = 0) in the calibrator
C list, then set max = 45 arcsec
if (max .eq. 0) xmax = 45.0
if (fwhm .gt. xmax) go to 100
C
C Is the GBT beam too small for this calibrator?
xmin = min
if (fwhm .lt. xmin) go to 100
C
C Estimate the calibrator flux (Jy) at frequency fghz
C by interpolating or extrapolating measured flux densities
C listed in the calibrator data file.
C
sjyest = 0.
C
C Is the frequency fghz below 1.4 GHz?
if (iband .eq. 1) then
C Is there a measured flux at 1.4 GHz?
if (sjy(iband) .gt. 0.) then
C Extrapolate the measured 1.4 GHz flux assuming alpha = 0.7
C (that is, flux declines as the 0.7 power of frequency, a
C normal "steep-spectrum" source).
sjyest = (fghz / 1.4)**(-0.7) * sjy(1)
C Is there a measured flux at a higher frequency?
do 10 i = 2, 6
if (sjy(i) .gt. 0.) then
C Extrapolate using measured alpha instead of assuming 0.7
alpha = -alog(sjy(i)/sjy(1))/alog(freq(i)/freq(1))
sjyest = (fghz / 1.4)**(-alpha) * sjy(1)
go to 90
end if
10 continue
end if
go to 90
end if
C
C Is the frequency fghz between 1.4 and 86 GHz?
if (iband .le. 7) then
C Get nonzero flux at highest frequency below fghz
do 20 i = 1, iband-1
if (sjy(i) .gt. 0.) lowband = i
20 continue
C Extrapolate sjy(lowband) using alpha = 0.7
sjyest = sjy(lowband) * (fghz / freq(lowband))**(-0.7)
C Get measured flux (if any) at lowest frequency above fghz
ihiband = 0
do 30 i = 7, iband, -1
if (sjy(i) .gt. 0.) ihiband = i
30 continue
C Interpolate between lowband and ihiband assuming a
C power-law spectrum
if (ihiband .ne. 0) then
alpha = -alog(sjy(lowband)/sjy(ihiband)) /
* alog(freq(lowband)/freq(ihiband))
sjyest = (fghz / freq(lowband))**(-alpha) *
* sjy(lowband)
end if
go to 90
end if
C
C Is the frequency fghz above 86 GHz?
if (iband .eq. 8) then
C Extrapolate highest-frequency nonzero flux assuming alpha = 0.7
do 40 i = 1, 7
if (sjy(i) .gt. 0.) lowband = i
40 continue
sjyest = sjy(lowband) * (fghz / freq(lowband))**(-0.7)
end if
90 continue
C
C Is the calibrator too weak at frequency fghz?
if (sjyest .lt. sjymin) go to 100
C
C Is the calibrator a "gold" calibrator?
if (gold .eq. 'Y' .and. goldmark .ne. 'P') go to 100
C
C ********** start new section added 2008 Jun 3
C Is the slew time larger than the limit? Is the elevation too high?
C (skip this if hlst < 0.)
if (hlst .ge. 0.) then
C
C calculate the hour angle (radians) of the calibrator
hacal = radlst - racal
C calculate the current elevation (radians) of the calibrator
elcal = asin (sin (tellat) * sin (deccal) +
* cos (tellat) * cos (deccal) * cos (hacal) )
C reject calibrators with elevation angle > 88 deg = 1.5359 rad
if (elcal .gt. 1.5359) go to 100
C calculate the current azimuth (radians) of the calibrator
denom = cos (tellat) * sin (deccal) -
* sin (tellat) * cos (deccal) * cos (hacal)
xnum = -sin(hacal) * cos (deccal)
azcal = atan2(xnum,denom)
C convert azimuth range from (-pi to +pi) to (0 to 2pi)
if (azcal .lt. 0.) azcal = azcal + 2. * pi
C (see http://aa.usno.navy.mil/faq/docs/Alt_Az.php)
C Estimate the current GBT slew time (minutes) assuming slew rates of
C 40 deg/min = 0.6981 rad/min in azimuth and
C 20 deg/min = 0.3491 rad/min in elevation
zslew = abs ( (elcal - elcenter) / 0.3491 )
azdiff = abs (azcal - azcenter)
if (azdiff .gt. pi) azdiff = 2. * pi - azdiff
azslew = azdiff / 0.6981
tslew = zslew
if (azslew .gt. tslew) tslew = azslew
C Limit slew time (minutes) from search center to calibrator to the larger of
C 1 minute or
C the search radius in degrees divided by 10 (half the elevation slew rate).
tslewmax = searchdeg / 10.
if (tslewmax .lt. 1.) tslewmax = 1.
if (tslew .gt. tslewmax) go to 100
end if
C ********** end new section added 2008 Jun 3
C
C The calibrator is suitable. Reread and print out the calibrator data.
ncals = ncals + 1
backspace 21
C write the output header
if (ncals .eq. 1) then
if (fghz .lt. 99.5) write (*, 1070) fghz
if (fghz .ge. 99.5) then
ifghz = fghz
write (*, 1071) ifghz
end if
write (*, *) 'HHMM+DDMM HH MM SS.SSSS DDD MM SS.SSS Jy
*deg'
write (*, *) '----------------------------------------------------
*----'
end if
read (21, 1160) col
write (*, 1170) (col(i), i = 1, 9), (col(i), i = 10, 39),
* sjyest, offsetdeg
100 continue
C
C end of calibrator file reached
200 continue
C
C ********** start new section added 2008 Jun 3
C ensure that at least one calibrator is found
C by relaxing the search criteria on flux and search radius
if (ncals .eq. 0) then
iwarning = iwarning + 1
C Double search area
searchrad = searchrad * sqrt(2.)
offsetmax = 2. * sin (searchrad / 2.)
C Lower minimum flux
sjymin = sjymin / sqrt(2.)
rewind 21
go to 5
end if
C ********** end new section added 2008 Jun 3
C
write (*, *) '----------------------------------------------------
*----'
close (unit = 21)
C normal end
return
C
C format statements
1024 format (' (rms = ', f5.3,
* ' Jy/beam at this frequency).')
1030 format (1x, i2.2, 1x, i2.2, 1x, f5.2, 1x,
* a1, i2.2, 1x, i2.2, 1x, f4.1)
1070 format (' IAU NAME RA (J2000) DEC S(',
*f4.1,' GHz) Offset')
1071 format (' IAU NAME RA (J2000) DEC S(',
*i3,' GHz) Offset')
1110 format (a9, 1x,
* i2, 1x, i2, 1x, f7.4, 1x, a1, i2, 1x, i2, 1x, f6.3,
* 1x, i3, 1x, i3, 6(1x, f6.2), 1x, a1)
1160 format (90a1)
1170 format (1x, 9a1, 2x, 30a1, 2x, f5.2, 2x, f5.1)
end
-- FrankGhigo - 01 Jul 2008 |