import pyfits as pf import os import time import numpy as np class TwentyMeterRcvrPafFits : def __init__ ( self, main_directory, write_flag=False ) : if os.name == 'posix' : self.dir = main_directory + '/RcvrPAF/' else : self.dir = main_directory + '\\RcvrPAF\\' if not os.path.exists(self.dir) : os.makedirs(self.dir) self.software_start_time = self.make_time_string() self.write_flag = write_flag if not write_flag : print 'RcvrPafFits opened READ ONLY' self.file_name = '' self.file_open = False self.fits_struct_open = False self.fits_struct_unsaved = False self.prihdr = False def make_time_string ( self, file_name=False, year=False, month=False, day=False, hour=False, min=False, sec=False ) : if not year : td = time.gmtime() year = td[0] month = td[1] day = td[2] hour = td[3] min = td[4] sec = td[5] if file_name : tstr = '%04d_%02d_%02d_%02dh%02dm%02d' % (td[0], td[1], td[2], td[3], td[4], td[5]) else : tstr = '%04d-%02d-%02dT%02d:%02d:%02d' % (td[0], td[1], td[2], td[3], td[4], td[5]) return tstr def open_file ( self, file_name ) : if self.fits_struct_open : print 'A RcvrPAF FITS structure is already open' return if self.file_open : print 'A RcvrPAF FITS file is already open' % (self.file_name) return if not os.path.exists(self.dir + file_name) : print '%s does not exist' % (self.dir + file_name) return self.hdulist = pf.open(self.dir + file_name) self.file_name = file_name self.file_open = True def close_file ( self ) : if self.fits_struct_unsaved : self.save_fits_struct() self.hdulist.close() self.file_name = '' self.file_open = False def open_new_fits_struct ( self, date_time_str ) : if self.file_open : print 'A RcvrPAF FITS file is already open' % (self.file_name) return if self.fits_struct_open : print 'A RcvrPAF FITS structure is already open' return self.fits_struct_open = True self.file_name = date_time_str + '.fits' hdu0 = pf.PrimaryHDU() prihdr = hdu0.header self.prihdr = prihdr prihdr.add_comment('FITS (Flexible Image Transport System) format is defined in Astronomy') prihdr.add_comment('and Astrophysics, volume 376, page 359; bibcode: 2001A&A...376..359H') prihdr.update('ORIGIN', 'NRAO Green Bank') prihdr.update('INSTRUME', 'Measurements', 'device or program of origin') prihdr.update('GBTMCVER', '00.0 ', 'telescope control software release') prihdr.update('FITSVER', '0.0 ', 'FITS definition version for this device') prihdr.update('DATEBLD', '0000-00-00T00:00:00', 'time program was released') prihdr.update('SIMULATE', 0, 'Is the instrument in simulate mode') prihdr.update('DATE-OBS', self.software_start_time, 'Manager parameter startTime') prihdr.update('TIMESYS', 'UTC', 'time scale specification for DATE-OBS') prihdr.update('TELESCOP', 'NRAO_20M', 'Green Bank 20-meter') prihdr.update('OBJECT', 'none', 'Manager parameter source') prihdr.update('PROJID', 'A20Mxxxx', 'Manager parameter projectId') prihdr.update('OBSID', 'calibration', 'Manager parameter scanId') prihdr.update('SCAN', 0, 'integer scan identifier') prihdr.update('RECEIVER', 'RcvrPAF', 'receiver name') prihdr.update('MIN_FREQ', 1.3e+09, '') prihdr.update('MAX_FREQ', 1.75e+09, '') prihdr.update('NELEM', 19, 'number of elements') prihdr.update('LO1FREQ', 2.15e9, 'LO1 frequency in Hz') prihdr.update('IF1FREQ', 4.0e8, 'IF1 center frequency in Hz') prihdr.update('IF1BW', 5.0e6, 'IF1 bandwidth in Hz') prihdr.update('LO2FREQ', 3.96e8, 'LO2 frequency in Hz') prihdr.update('IF2FREQ', 2.81e6, 'IF2 center frequency in Hz') prihdr.update('IF2BW', 4.02e5, 'IF2 bandwidth in Hz') # double sample_clock_freq; // sample clock freq in MHz (2x sample rate) # double sample_clock_amplitude;// sample_clock amplitude in dBm # double marker_pulse_freq; // marker pulse frequency in MHz # double marker_pulse_max; // marker pulse maximum in volts # double marker_pulse_min; // marker pulse manimum in volts # int marker_pulse_cycles; // number of marker pulse cycles prihdr.update('SAMPCLKF', 2.5e6, 'sample clock freq Hz (2x samp rate') prihdr.update('SAMPCLKA', 0.0, 'sample_clock amplitude in dBm') prihdr.update('MKRPFREQ', 2e5, 'marker pulse frequency in Hz') prihdr.update('MKRAMAX', 5.0, 'marker pulse maximum in volts') prihdr.update('MKRAMIN', 0.0, 'marker pulse minimum in volts') prihdr.update('MKRCYCLS', 200, 'number of marker pulse cycles') prihdr.update('ETAL', 0.98, 'rear spillover and scattering efficiency') prihdr.update('APEREFF', 0.70, 'aperture efficiency') prihdr.update('BEAMEFF', 0.80, 'main beam efficiency') prihdr.update('TAUZENIT', 0.003, 'opacity per unit air mass') # tbhdu_no col -- np.array -- # calib coeff dimensions: center_freq, hot/cold, coef_index, freq_chan num_correl = 21 * 20 / 2 num_chan = 256 real_imag = 2 data = np.zeros(num_chan * num_correl * real_imag, dtype=np.float32) data.shape = (num_correl, num_chan, real_imag) for cor in range(num_correl) : data[cor,:,0] = np.arange(num_chan, dtype=np.float32) data[cor,:,1] = -np.arange(num_chan, dtype=np.float32) data.shape = (num_chan * num_correl * real_imag) x_adc_list = np.array([0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 1, 3, 5, 7, 9, 11, 13, 15, 17, 19]) + 1 y_adc_list = np.zeros(len(x_adc_list), dtype=np.int16) #y_adc_list = np.array([]) hdu1 = self.make_element_map_table(x_adc_list, y_adc_list) hdu2 = self.make_cal_table(test_date='2011-07-01', center_freq=1.45e9, bandwidth=1.0e7, num_correl=num_correl, num_chan=num_chan, engineer='B. Simon', technician='K. Ward', hot_matrix_x=data, cold_matrix_x=data, hot_matrix_y=data, cold_matrix_y=data) hdu3 = self.make_cal_table(test_date='2011-07-01', center_freq=1.55e9, bandwidth=1.0e7, num_correl=num_correl, num_chan=num_chan, engineer='B. Simon', technician='K. Ward', hot_matrix_x=data, cold_matrix_x=data, hot_matrix_y=data, cold_matrix_y=data) self.hdulist = pf.HDUList([hdu0, hdu1, hdu2, hdu3]) self.fits_struct_unsaved = True def make_element_map_table ( self, x_adc_list, y_adc_list ) : if len(x_adc_list) != len(y_adc_list) : print 'Element maps must be same length (%d, %d)' % \ (len(x_adc_list) != len(y_adc_list)) return col1 = pf.Column(name='XELEMMAP', format='%dI' % (len(x_adc_list)), unit='none', array=[x_adc_list]) col2 = pf.Column(name='YELEMMAP', format='%dI' % (len(y_adc_list)), unit='none', array=[y_adc_list]) cols = pf.ColDefs([col1, col2]) hdu = pf.new_table(cols) hdu.name = 'ELEM_MAP' hdr = hdu.header hdr.add_comment('maps element numbers to aggregated data ADC channels') hdr.add_comment('array index = element number,') hdr.add_comment(' value = ADC channel number, base 1') return hdu def make_cal_table ( self, test_date, center_freq, bandwidth, num_correl, num_chan, engineer, technician, hot_matrix_x, cold_matrix_x, hot_matrix_y, cold_matrix_y ) : col1 = pf.Column(name='HOT_X', format='%dE' % \ (num_correl * num_chan * 2), unit='corcoef', array=[hot_matrix_x]) col2 = pf.Column(name='COLD_X', format='%dE' % \ (num_correl * num_chan * 2), unit='corcoef', array=[cold_matrix_x]) col3 = pf.Column(name='HOT_Y', format='%dE' % \ (num_correl * num_chan * 2), unit='corcoef', array=[hot_matrix_y]) col4 = pf.Column(name='COLD_Y', format='%dE' % \ (num_correl * num_chan * 2), unit='corcoef', array=[cold_matrix_y]) cols = pf.ColDefs([col1, col2, col3, col4]) hdu = pf.new_table(cols) hdu.name = 'CAL_MTRX' hdr = hdu.header hdr.update('TESTDATE', test_date, 'date test data acquired') hdr.update('CNTRFREQ', center_freq, 'center frequency of sampled passband, Hz') hdr.update('BANDWDTH', bandwidth, 'bandwidth of correlated spectra, Hz') hdr.update('MTXSHAPE', '[%d,%d,%d]' % (num_correl, num_chan, 2), '[num_correl, num_chan, real_imag]') hdr.update('ENGINEER', engineer, 'receiver engineer') hdr.update('TECH', technician, 'receiver technician') hdr.update('EXTNAME', 'CAL_MTRX', 'name of this binary table extension') return hdu def save_fits_struct ( self ) : if not self.write_flag : print 'TwentyMeterRcvrPafFits class not open for writing' return file_path = self.dir + self.file_name if os.path.exists(file_path) : os.remove(file_path) self.hdulist.writeto(file_path) self.fits_struct_unsaved = False def show_hdu_info ( self ) : if (not self.file_open) and (not self.fits_struct_open): print 'RcvrPAF file not open' return print self.hdulist.info() def show_primary_card_list ( self ) : if (not self.file_open) and (not self.fits_struct_open): print 'RcvrPAF file not open' return print self.hdulist[0].header.ascardlist() def show_hdu_card_list ( self, hdu_num ) : if (not self.file_open) and (not self.fits_struct_open): print 'RcvrPAF file not open' return print self.hdulist[hdu_num].header.ascardlist() #rcvr_paf = TwentyMeterRcvrPafFits('.', False) rcvr_paf = TwentyMeterRcvrPafFits('.', True) #rcvr_paf.open_file('2010_11_07_00h39m08.fits') rcvr_paf.open_new_fits_struct('2010_11_07_00h39m08') rcvr_paf.show_hdu_info() #rcvr_paf.update_keyword('SUBMOTIN', 'somethng') rcvr_paf.show_primary_card_list() rcvr_paf.show_hdu_card_list(1) rcvr_paf.close_file()