# This Python module retrieves autocorellation data from a autocorrelation # spectrometer 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 autocorrelation correction functions from the # 'VanVleck' module, and the FFT routine from the Numeric FFT module. 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. import sys if not sys.path.__contains__('/users/rfisher/Applications/Python/AcfPowerCor') : sys.path.append('/users/rfisher/Applications/Python/AcfPowerCor') import acfpwrcor from get_fits import * #from pgplot import cpgplot import ppgplot import Numeric import FFT from gbt.spectrometer import cspectrometer # The arguments are the project name, e.g., 'AGBT02A_069', the FITS file # name of the spectrometer data to be plotted, e.g., # '2003_05_24_20:44:16A.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_ac_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 + '/' acft = get_fits(base_path + 'Spectrometer/' + file) # the spectrometer FITS file name has a bank designato suffix that needs # to be removed for use by other device file names base_file = get_base_file_name(file); loft = get_fits(base_path + 'LO1A/' + base_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 spectrum bandwidth in MHz. bandwidth = acft['PORT']['Column']['BANDWDTH'][0] * 1e-6 # get the dimensions of the spectrum data array and give them physical # meaning data_shape = acft['DATA']['Column']['DATA'].shape num_records = data_shape[0] # number of integrations in the scan num_phases = data_shape[1] # switching phases, i.e., cal on, cal off num_rcvrs = data_shape[2] # number of IF channels into spectral processor num_chan = data_shape[3] # number of lags in each ACF # compute the array of frquency offsets of the spectral channels with # respect to the reference channel, num_chan / 2 - 1, where the channel # numbers run from 0 to num_chan - 1. The number of frequency # channels after Fourier transforming the ACF is equal to the number of # ACF lags. ch_spacing = bandwidth / num_chan freq_offset = (Numeric.arange(0,num_chan) - (num_chan / 2) + 1) * 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. # The effective sideband of the autocorrelation spectrometer depends on its # bandwidth, since there is an extra lower-sideband conversion before the # 50- and 12.5-MHz samplers, and the samplers use one of the aliases # rather than baseband. The following rules apply: inversion_flag = 0 bwkhz = int(bandwidth * 1e3 + 0.5) # kHz if (bwkhz == 12500) or (bwkhz == 800000) : inversion_flag = 1 if ((rf_sideband == 'LOWER') and (not inversion_flag)) or \ ((rf_sideband == 'UPPER') and (inversion_flag)) : 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 acft['STATE']['Column']['CAL'][0] == 1 : on_state = 1 off_state = 0 else : on_state = 0 off_state = 1 # get the number of sampler levels used to use in thw correlation # quantization (Van Vleck) corrections. levels = acft['PORT']['Column']['LEVEL'][0] # open the PGPLOT window ppgplot.pgopen("/xwindow") # draw the fequency and temperature axes and label them ppgplot.pgenv(freq[0], freq[-1], 0, 30.0, 0, 1) ppgplot.pglab('Frequency in MHz', 'System Temperature in Kelvins', 'Autocorrelation Spectrometer 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 += get_spectrum(acft['DATA']['Column']['DATA'][rec_num, on_state, rcvr_num, :], levels) cal_off += get_spectrum(acft['DATA']['Column']['DATA'][rec_num, off_state, rcvr_num, :], levels) # renormalize the spectra cal_on = cal_on / num_records cal_off = cal_off / num_records # do a Hanning convolution cal_on[1:-1] = (0.5 * (cal_on[1:-1] + 0.5 * (cal_on[0:-2] + cal_on[2:]))).astype(Numeric.Float32) cal_off[1:-1] = (0.5 * (cal_off[1:-1] + 0.5 * (cal_off[0:-2] + cal_off[2:]))).astype(Numeric.Float32) # reverse the spectra, if necessary if reverse_flag: cal_on = cal_on[-1::-1] cal_off = cal_off[-1::-1] # 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) # set the plotting color ppgplot.pgsci(color) # draw the spectrum ppgplot.pgline(freq, tsys) # ppgplot.pgline(freq, 50.0 * cal_off) # ppgplot.pgline(freq, 50.0 * cal_on) # 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 # This function removes the GBT spectrometer bank suffix from the FITS file # name. def get_base_file_name ( file ) : if (string.find(file, 'A')) : return string.join(string.split(file, 'A'), '') if (string.find(file, 'B')) : return string.join(string.split(file, 'B'), '') if (string.find(file, 'C')) : return string.join(string.split(file, 'C'), '') if (string.find(file, 'D')) : return string.join(string.split(file, 'D'), '') else : return file # This function applies the quantization and power linearization correction # to the ACF and Fourier transforms it to return a power spectrum. # At the moment the corrections are a bit of a kluge because the # cspectrometer module does not return a normalized ACF as expected by # the power correction. def get_spectrum ( acf, levels ) : # 'optimum' is the maximum sensitivity threshold level for the ampler if (levels == 3) : optimum = 0.5405 # get the corrected ACF values. The returned ACF is normalized to a # zero-lag value of 1. pwr_cor = acfpwrcor.power_correction(acf[0], 3) cspectrometer.schwab3Level(acf) else : # the 9-level ACF received from the GBT spectrometer # is a factor of 16 smaller than the quantization mathematics assume acf = (16.0 * acf).astype(Numeric.Float32) pwr_cor = acfpwrcor.power_correction(acf[0], 9) optimum = 3.401 # get the corrected ACF values. The returned ACF is normalized to a # zero-lag value of 1. cspectrometer.schwab9Level(acf) # normalize the ACF to the zero-lag value as expected by the power # linearization correction acf = acf / acf[0] # apply the linearization correction acf = (acf * pwr_cor).astype(Numeric.Float32) # create a symmetric array out of the ACF by mirror image around its # zero-lag value dacf = Numeric.zeros(2 * len(acf), Numeric.Float32) dacf[len(acf):(2 * len(acf))] = acf dacf[1:len(acf)] = acf[-1:0:-1] dacf[0] = dacf[1]; # do the Fourier transform spec = FFT.real_fft(dacf) # return the amplitude of the complex vectors return (((spec[1:] * Numeric.conjugate(spec[1:])).real)**0.5).astype(Numeric.Float32) # To test uncomment the plot call below and run : # import plot_actsys # or # reload(plot_actsys) print 'The quantization corrections for 40 spectra takes about 18 seconds:' plot_ac_tsys('AGBT02A_069', '2003_05_24_20:44:16A.fits', [1.5, 1.5])