# This Python module retrieves spectrum data from a spectral processor GBT # FITS file and frequency information from the corresponding LO1A file and # computes and displays the spectrum of the system temperature as a function # of sky frequency. It uses the FITS retrieval routine in the module # 'get_fits'. The FITS files are assumed to reside in the directory # /home/gbtdata under the project name directory. A Python interface to # the PGPLOT library is used for plotting. from get_fits import * #from pgplot import cpgplot import ppgplot import Numeric # The arguments are the project name, e.g., 'AGBT02A_069_01', the FITS file # name of the data to be plotted, e.g., '2003_06_17_09:22:34.fits', and the # calibration noise source values for each receiver channel, e.g., [1.7, 1.7]. # The noise values are defaulted to 1.0. We assume that the only phase # switching in the data is cal no, cal off. def plot_sp_tsys ( project, file, tcal=[1.0, 1.0] ) : # assemble the directory paths and retrieve the spectral processor and # LO1A data dictionaries. See the 'get_fits' module for the details of # the directory structure and a routine for listing a FITS file's # contents base_path = '/home/gbtdata/' + project + '/' spft = get_fits(base_path + 'SpectralProcessor/' + file) loft = get_fits(base_path + 'LO1A/' + file) # get the first mixer's sideband string rf_sideband = loft['main_header']['SIDEBAND'] # compute the observed center frequency from the LO frequency at the # beginning of the scan, the first IF center frequency, and the sideband # and convert from Hz to MHz if rf_sideband == 'LOWER' : ctr_freq = (loft['LO1TBL']['Column']['LO1FREQ'][0] - loft['main_header']['IFFREQ']) * 1e-6 else : ctr_freq = (loft['LO1TBL']['Column']['LO1FREQ'][0] + loft['main_header']['IFFREQ']) * 1e-6 # get the spacing of the spectral processor spectrum in MHz ch_spacing = spft['RECEIVER']['Column']['FREQRES'][0] * 1e-6 # get the spectral processor input sideband: 0 = upper so that increasing # output spectrum channel number corresponds to frequency at the spectral # processor input. if_sideband = 1 is the reverse. if_sideband = spft['RECEIVER']['Column']['IFSIDE'][0] # get the dimensions of the spectrum data array and give them physical # meaning data_shape = spft['DATA']['Column']['DATA'].shape num_records = data_shape[0] # number of integrations in the scan num_rcvrs = data_shape[1] # number of IF channels into spectral processor num_phases = data_shape[2] # switching phases, i.e., cal on, cal off num_chan = data_shape[3] # number of frequency bins in each spectrum # compute the array of frquency offsets of the spectral channels with # respect to the reference channel, num_chan / 2, where the channel # numbers run from 0 to num_chan - 1. freq_offset = (Numeric.arange(0,num_chan) - (num_chan / 2)) * ch_spacing # the GBT has three mixer stages before the spectral processor input, # but the second and third mixer conversions are both always lower sideband # so the following frequency offset to RF frequency rules apply: if ((rf_sideband == 'LOWER') and (if_sideband == 0)) or \ ((rf_sideband == 'UPPER') and (if_sideband == 1)) : freq = ctr_freq - freq_offset else : freq = ctr_freq + freq_offset # if the frequency array runs from higher to lower, flip it and set a # flag for flipping the spectrum arrays as well reverse_flag = 0 if freq[0] > freq[-1] : freq = freq[-1::-1] reverse_flag = 1 # set spectrum array phases for the two cal states (1 = on, 0 = off) if spft['STATE']['Column']['CAL'][0,0] == 0 : on_state = 1 off_state = 0 else : on_state = 0 off_state = 1 # open the PGPLOT window ppgplot.pgopen("/xwindow") # draw the fequency and temperature axes and label them ppgplot.pgenv(freq[0], freq[-1], 0.0, 30.0, 0, 1) ppgplot.pglab('Frequency in MHz', 'System Temperature in Kelvins', 'Spectral Processor Tsys Plot ' + file) # set the color of the first spectrum plot and the Y position of the # text legend color = 2 # red txty = 1.5 # process the spectra for each receiver for rcvr_num in range(0,num_rcvrs) : # if not enough cal values have been specified, default it to 1.0 if len(tcal) < rcvr_num : tcal.append(1.0) # create empty data arrays for the cal on and cal of spectrum # accumulations cal_on = Numeric.zeros(num_chan, Numeric.Float32) cal_off = Numeric.zeros(num_chan, Numeric.Float32) # accumulate the spectra from all integrations in the scan for rec_num in range(0,num_records) : cal_on += spft['DATA']['Column']['DATA'][rec_num, rcvr_num, on_state, :] cal_off += spft['DATA']['Column']['DATA'][rec_num, rcvr_num, off_state, :] # renormalize the spectra cal_on = cal_on / num_records cal_off = cal_off / num_records # compute the channel by channel system temperature and clip # out-of-range values tsys = Numeric.clip(tcal[rcvr_num] * cal_off / (cal_on - cal_off), 0.0, 30.0) # reverse the spectrum, if necessary if reverse_flag: tsys = tsys[-1::-1] # set the ploting color ppgplot.pgsci(color) # draw the spectrum ppgplot.pgline(freq, tsys) # set the legend text X position txtx = freq[0] + 30 * ch_spacing # write the legend on the plot ppgplot.pgtext(txtx, txty, 'Rcvr ' + str(rcvr_num + 1)) # increment color and legend position color += 1 # 3 = green, 4 = blue, 5 = blue-green txty += 1.5 # close the plotting window after requesting keyboard acknowledgement ppgplot.pgask(1) ppgplot.pgclos() return # To test uncomment the plot call below and run : # import plot_tsys # or # reload(plot_tsys) #plot_sp_tsys('AGBT02A_069_01', '2003_06_17_09:22:34.fits', [1.7, 1.7])