NRAO Home  >  Green Bank  |  Wiki Topic:    GB > Main > FrankGhigo_CalFind3Notes

Users | Groups | Offices | Changes | Index | Contents | Notify | Search | Jump to Topic:

Notes on Condon's calfind3.f and conversion to python

Herewith a description of Condon's calfind3.f program, referred to in the following MR.

• 4C408 : Modify Auto Peak Calibrator Selection

• The actual Fortran code of calfind3

Here are the steps it goes through.

• read in and check the parameters
• for each calibrator in the list, determine :
  • arc distance between calibrator and requested location
  • the angular size limits for the calibrator
  • flux density at requested frequency
  • The "gold" parameter for this calibrator.
  • if calibrator satisfies criteria (distance, angular size, flux, gold), then also check the move time
  • find the distance in az and el between requested location and calibrator
  • estimate move times
• if calibrator fails all the criteria, relax the criteria and repeat.
• Print out all calibrators that meet the search criteria. This is a difference with the requirements of AutoPeak, in which only the closest is selected.

Inputs to calfind3:

• fghz -- observing frequency in GHz.
• searchdeg -- search radius in degrees. If >180, set to 180 degrees, which includes the whole sky.
• sjymin -- minimum calibrator flux density in Jy.
• gold -- if "y" choose only "gold standard" calibrators.
• RA and DEC -- the search center in J2000.
• hlst -- the LST in hours.

Relevant Inputs to AutoPeak/Focus

• location -- The RA and DEC of the center of the search. At present this Location must be in J2000. Perhaps we should implement conversion from other coordinate modes.
• frequency -- This is given in MHz for AutoPeak/Focus, instead of GHz.
• flux -- Same meaning as for calfind3.
• radius -- Same meaning as for calfind3.
• gold -- "True" to choose only "gold standard" calibrators.
• LST -- The user can give an explicit LST (probably we should allow hh:mm:ss format); but if no, the default is for the program to user the current LST.

OK, Here are some details.

The Calibrator Catalog

• The calibrator catalog is located in directory /home/astro-util/astridcats/
• At present the name of the catalog is "pcals4_1.cat", which is also known as "pointing.cat" and "default.cat".
• The catalog is assumed to have a specific format. It conforms to the Astrid catalog conventions.
  • The header information is "FORMAT=SPHERICAL" and "COORDMODE=J2000"
  • and "HEAD = NAME RA DEC BMIN BMAX S20 S6 S3.5 S2 S1.3 S0.7 Gold"
  • The catalog may have columns past column 12, which are ignored.
  • It is expected that in a year or so there will be a new version of the catalog with an additional column "S0.3" for the 3mm flux densities, between S0.7 and Gold.
• A "P" in the "Gold" column means that it is a gold standard calibrator. If the user has set gold=True, then only these calibrators are considered.

The angular distance

• Convert the search radius to cartesian offset on the unit sphere.
  • offset_c = 2 sin(radius/2)
• Convert the search center location to cartesian to avoid polar coordinate poles and RA wrap ambibuities.
  • x = cos(racenter)cos(deccenter)
  • y = sin(racenter)cos(deccenter)
  • z = sin(deccenter)
• Do likewise for the calibrator coordinates, to obtain (xcal, ycal, zcal)
• calculate cartesian distance
  • cal_offset = sqrt( (x-xcal)² + (y-ycal)² + (z-zcal)² )
• if cal_offset > offset_c reject this calibrator.
  • otherwise, set angular distance ADIST = 2 * arcsin( cal_offset/2) and accept this calibrator.
  • ADIST will be printed out if the calibrator meets all criteria.

Source size and beam size criteria

• get the beam size, full width half max (fwhm) in arcseconds : fwhm = 754/fghz
  • if the beam size >= 540 arcsec, then issue a warning and print the confusion rms:
  • "Warning! Large GBT beam; use only strong calibrators to minimize confusion."
  • The confusion rms in Jy is 0.030 * fghz^-2.7
  • if fwhm >= 540 then set fwhm to 539.
• Test the calibrator size values BMIN and BMAX vs the fwhm.
  • If BMAX is zero then use BMAX = 45 arcsec.
  • If fwhm > BMAX or fwhm < BMIN, then reject this calibrator.

Interpolate catalog flux density to the required frequency (fghz)

The frequencies Fx = (1.4, 4.9, 8.6, 15.0, 22.0, 43.0) correspond to fluxes (S20, S6, S3.5, S2, S1.3, S0.7).

• If fghz below 1.4 GHz
  • If S20 is zero, reject this calibrator.
  • Use S20 and Sx = the next non-zero flux in (S6, S3.5, S2, S1.3, S0.7); Fx is the corresponding frequency.
  • if no non-zero Sx, then extrapolate: cal_flux = S20 * (fghz/1.4)^-0.7
  • Otherwise, set alpha = ln(Sx/S20)/ln(Fx/1.4)
  • and cal_flux = S20 * (fghz/1.4)^alpha

• if fghz ≥ 1.4 and < 86 GHz
  • Find the catalog flux for the highest frequency below fghz with a non-zero flux. This is flux Sa at frequency Fa.
  • Extrapolate to observing frequency: cal_flux = Sa * (fghz/Fa)^-0.7
  • Find the catalog flux for the lowest frequency above fghz with non-zero flux. This will be Sb at frequency Fb.
  • If no non-zero flux is found in the previous step, then use the extrapolated cal_flux from 2 steps ago.
  • If there is both non-zero Sa and Sb then interpolate as follows:
  • alpha = ln(Sb/Sa)/ln(Fb/Fa) and cal_flux = Sa * (fghz/Fa)^alpha

• if fghz ≥ 86GHz then
  • find the highest frequency non-zero flux in the catalog. This will be Shigh at frequency Fhigh.
  • Extrapolate assuming alpha=-0.7:
  • cal_flux = Shigh * (fghz/Fhigh)^-0.7

• Note that when we get a revised catalog with a column for S0.3, these steps will change accordingly.

• If the cal_flux < specified minimum flux, then reject this calibrator.

If calibrator passes all tests so far, then check move times.

• If the specified LST (hlst) is less than zero, skip this step.
• Calculate AZs and ELs
  • Let the calibrator position be (racal, deccal), and hacal = hlst-racal.
  • Let the center position, i.e., the Location input be (racenter, deccenter), and hacenter = hlst-racenter.
  • The GBT Latitude (gbtlat) is 38.4331°
  • Calculate az and el for the calibrator (azcal and elcal) and of the input Location (azcenter, elcenter).
  • EL calculation:
    • el = arcsin{ sin(gbtlat)*sin(dec) + cos(gbtlat)*cos(dec)*cos(ha) }
  • AZ calculation:
    • denom = cos(gbtlat)*sin(dec) - sin(gbtlat)*cos(dec)*cos(ha)
    • xnum = -sin(ha)*cos(dec)
    • az = atan2(xnum,denom)
    • if az <0 then add 2*π (put az range to 0 to 2π)

• If calibrator elevation (elcal) >88° then reject this calibrator. [note: I may change this criterion to depend on frequency]
• If calibrator elevation (elcal) <7° then reject this calibrator. New requirement!
• Estimate slew times
  • Use azrate = 36°/min, and elrate = 18°/min
  • But note that in cold weather conditions (Temp <16F) the azrate is lowered to 18°/min.
  • elslew = abs(elcal-elcenter)/elrate
  • azdiff = abs(azcal-azcenter); if azdiff >π then set azdiff = 2π-azdiff
  • azslew = azdiff/azrate
  • tslew = max(azslew, elslew)
• Limit slew time in minutes from search center to calibrator to the larger of one minute or the search radius in degrees divided by 10.
  • slewlim = radius(in degrees)/10.
  • if slewlim <1, then set slewlim = 1.
  • if tslew > slewlim then reject this calibrator.

If No Calibrator passes all tests, then relax the criteria.

• Increase search radius by a factor sqrt(2)
• decrease flux limit by a factor sqrt(2)
• Issue a warning message that search criteria have been changed, and print the new radius and flux limit.
• repeat the search through all calibrators in the list.
• repeat the relaxation of criteria until a calibrator is found.

Print info for all selected calibrators.

-- FrankGhigo - 01 Jul 2008