diff --git a/RTS/2.2-T_sys_T_nd/T_sys_T_nd_red.py b/RTS/2.2-T_sys_T_nd/T_sys_T_nd_red.py index e5c32f8f..1539990b 100644 --- a/RTS/2.2-T_sys_T_nd/T_sys_T_nd_red.py +++ b/RTS/2.2-T_sys_T_nd/T_sys_T_nd_red.py @@ -2,8 +2,8 @@ import matplotlib matplotlib.use('Agg') import argparse -#from katsdpscripts.RTS import diodelib -from katsdpscripts.reduction import diodelib +#from katsdpscripts.RTS import diodeib +import diodelib def parse_arguments(): parser = argparse.ArgumentParser(description=" This produces a pdf file with graphs verifying the ND model and Tsys for each antenna in the file") @@ -14,23 +14,22 @@ def parse_arguments(): parser.add_argument("--error_bars", action='store_true',default=False, help="Include error bars - Still in development") parser.add_argument("--off_target", default='off1', help="which of the two off targets to use") parser.add_argument("--write_nd", action='store_true', default=False, help="Write the Noise Diode temp to a file") + parser.add_argument("--rfi_mask",default="/var/kat/katsdpscripts/RTS/rfi_mask.pickle",help="Default is /var/kat/katsdpscripts/RTS/rfi_mask.pickle") + parser.add_argument("--receptor",default='',help="Receptor number to process e.g. 'm053', default is all") parser.add_argument("filename", nargs=1) - + args,unknown = parser.parse_known_args() if args.filename[0] == '': raise RuntimeError('Please specify an h5 file to load.') - + return args,unknown if __name__ == "__main__": args,unknown = parse_arguments() - print(unknown) + print unknown kwargs = dict(zip(unknown[0::2],unknown[1::2])) - diodelib.read_and_plot_data(args.filename[0],args.output_dir,args.pdf,args.Ku,args.verbose,args.error_bars,args.off_target,args.write_nd,**kwargs) - - - + diodelib.read_and_plot_data(args.filename[0],args.output_dir,args.pdf,args.Ku,args.verbose,args.error_bars,args.off_target,args.write_nd,args.rfi_mask,args.receptor,**kwargs) diff --git a/katsdpscripts/RTS/diodelib.py b/katsdpscripts/RTS/diodelib.py index 27888228..6783f270 100644 --- a/katsdpscripts/RTS/diodelib.py +++ b/katsdpscripts/RTS/diodelib.py @@ -1,6 +1,4 @@ #!/usr/bin/python -from __future__ import absolute_import -from __future__ import print_function import katdal as katfile import scape import numpy as np @@ -10,14 +8,132 @@ from katsdpscripts import git_info from scipy.signal import medfilt import logging -import scape -from six.moves import zip +import scikits.fitting as fit +import time, ephem + + +def save_ND(diode_filename,file_base,freq,Tdiode_pol ): + outfilename = diode_filename.split('/')[-1] + outfile = file(outfilename, 'w') + outfile.write('#Data from %s\n# Frequency [Hz], Temperature [K]\n'%file_base) + # Write CSV part of file + outfile.write(''.join(['%s, %s\n' % (entry[0], entry[1]) for entry in zip(freq,medfilt(Tdiode_pol))])) + outfile.close() + return outfilename + + +def get_nd_on_off(h5,buff = 2,log=None): + on = h5.sensor['Antennas/%s/nd_coupler'%(h5.ants[0].name)] + off = ~h5.sensor['Antennas/%s/nd_coupler'%(h5.ants[0].name)] + n_on = np.tile(False,on.shape[0]) + n_off = np.tile(False,on.shape[0]) + if not any(on): + if log is not None: + log.critical('No noise diode fired during track of %s'%target) + else : + print('No noise diode fired during track of %s'%target) + else: + #jumps = (np.diff(on).nonzero()[0] + 1).tolist() + n_on[on.nonzero()[0][buff:-buff]] = True + n_off[off.nonzero()[0][buff:-buff]] = True + return n_on,n_off + + +def plot_Tsys_eta_A(freq,Tsys,eta_A,TAc,Ku=False,Tsys_std=None,ant = '', file_base='.'): + fig = plt.figure(figsize=(20,5)) + pols = ['v','h'] + for p,pol in enumerate(pols) : + fig.add_subplot(1,2,p+1) + ax = plt.gca() + ax.text(0.95, 0.01,git_info(), horizontalalignment='right',fontsize=10,transform=ax.transAxes) + plt.title('%s $T_{sys}/eta_{A}$: %s pol: %s'%(ant,str(pol).upper(),file_base)) + plt.ylabel("$T_{sys}/eta_{A}$ [K]") + plt.xlabel('f [MHz]') + #if p == ant_num * 2 -1: plt.ylabel(ant) + if Tsys_std[pol] is not None : + plt.errorbar(freq/1e6,Tsys[pol],Tsys_std[pol],color = 'b',linestyle = '.',label='Measurement') + plt.plot(freq/1e6,Tsys[pol]/eta_A[pol],'b.',label='Measurement: Y-method') + if not(Ku): plt.plot(freq/1e6,TAc[pol]/eta_A[pol],'c.',label='Measurement: ND calibration') + plt.axhline(np.mean(Tsys[pol]/eta_A[pol]),linewidth=2,color='k',label='Mean: Y-method') + if freq.min() < 2090e6: + D = 13.5 + Ag = np.pi* (D/2)**2 # Antenna geometric area + spec_Tsys_eta = np.zeros_like(freq) + plt.ylim(15,50) + #plt.xlim(900,1670) + spec_Tsys_eta[freq<1420e6] = 42 # [R.T.P095] == 220 + spec_Tsys_eta[freq>=1420e6] = 46 # [R.T.P.096] == 200 + plt.plot(freq/1e6, spec_Tsys_eta,'r',linewidth=2,label='PDR Spec') + plt.plot(freq/1e6,np.interp(freq/1e6,[900,1670],[(64*Ag)/275.0, + (64*Ag)/410.0]),'g',linewidth=2,label="275-410 m^2/K at Receivers CDR") + plt.grid() + plt.legend(loc=2,fontsize=12) + + return fig -def read_and_plot_data(filename,output_dir='.',pdf=True,Ku = False,verbose = False,error_bars=False,target='off1',write_nd=False,**kwargs): - print('inside',kwargs) +def plot_nd(freq,Tdiode,nd_temp,ant = '', file_base=''): + fig = plt.figure(figsize=(20,5)) + pols = ['v','h'] + for p,pol in enumerate(pols) : + fig.add_subplot(1,2,p+1) # + ax = plt.gca() + ax.text(0.95, 0.01,git_info(), horizontalalignment='right',fontsize=10,transform=ax.transAxes) + plt.title('%s Coupler Diode: %s pol: %s'%(ant,str(pol).upper(),file_base)) + plt.ylim(0,50) + plt.ylabel('$T_{ND}$ [K]') + #plt.xlim(900,1670) + plt.xlabel('f [MHz]') + #plt.ylabel(ant) + plt.axhspan(14, 35, facecolor='g', alpha=0.5) + plt.plot(freq/1e6,Tdiode[pol],'b.',label='Measurement: Y-method') + plt.plot(freq/1e6,nd_temp[pol],'k.',label='Model: EMSS') + plt.grid() + plt.legend() + return fig + +def plot_ts(h5,on_ts=None): + import scape + fig = plt.figure(figsize=(20,5)) + a = h5.ants[0] + nd = scape.gaincal.NoiseDiodeModel(freq=[856,1712],temp=[20,20]) + d = scape.DataSet(h5,nd_h_model = nd, nd_v_model=nd ) + scape.plot_xyz(d,'time','amp',label='Average of the data') + if on_ts is not None : + on = on_ts + else : + on = h5.sensor['Antennas/'+a.name+'/nd_coupler'] + ts = h5.timestamps - h5.timestamps[0] + plt.plot(ts,np.array(on).astype(float)*4000,'g',label='katdal ND sensor') + plt.title("Timeseries for antenna %s - %s"%(a.name,git_info())) + plt.legend() + return fig + +def get_fit(x,y,order=0,true_equals=None): + """true_equals is a test that returns bool values to be fitted""" + if true_equals is None : + y = y.astype(float) + else : + y = y==true_equals + return fit.PiecewisePolynomial1DFit(order).fit(x,y)#(timestamps[:]) + +def Dmoon(observer): + """ The moon's apparante diameter changes by ~10% during over a year, and + an x% error in Dmoon leads to systemtic error in Thot scaling as 2*x/lambda! + @param observer: either an ephem.Observer or a suitable date (then at earth's centre). + @return [rad] diameter of the moon, from the earth on the specified date (ephem.Observer)""" + observer = observer if isinstance(observer,ephem.Observer) else ephem.Date(observer) + moon = ephem.Moon(observer) + D = moon.size/3600. # deg + return D*np.pi/180. # rad + + +def read_and_plot_data(filename,output_dir='.',pdf=True,Ku = False, + verbose = False,error_bars=False,target='off1', + write_nd=False,rfi_mask='/var/kat/katsdpscripts/RTS/rfi_mask.pickle',receptor='',**kwargs): + print 'inside',kwargs file_base = filename.split('/')[-1].split('.')[0] nice_filename = file_base + '_T_sys_T_nd' - + # Set up logging: logging everything (DEBUG & above), both to console and file logger = logging.root logger.setLevel(logging.DEBUG) @@ -26,59 +142,77 @@ def read_and_plot_data(filename,output_dir='.',pdf=True,Ku = False,verbose = Fal fh.setFormatter(logging.Formatter('%(levelname)s: %(message)s')) logger.addHandler(fh) logger.info('Beginning data processing with:\n%s'%git_info('standard')) - - h5 = katfile.open(filename,**kwargs) - if verbose: logger.debug(h5.__str__()) - ants = h5.ants - - pickle_file = open('/var/kat/katsdpscripts/RTS/rfi_mask.pickle') - rfi_static_flags = pickle.load(pickle_file) - pickle_file.close() - edge = np.tile(True,4096) - edge[slice(211,3896)] = False - static_flags = np.logical_or(edge,rfi_static_flags) + if Ku: logger.debug("Using Ku band ... unsetting L band RFI flags") - h5.spectral_windows[0].centre_freq = 12500.5e6 + h5 = katfile.open(filename,centre_freq = 12500.5e6 , **kwargs) + length = h5.shape[1] # Don't subtract half a channel width as channel 0 is centred on 0 Hz in baseband - h5.spectral_windows[0].channel_freqs = h5.spectral_windows[0].centre_freq + h5.spectral_windows[0].channel_width * (np.arange(h5.spectral_windows[0].num_chans) - h5.spectral_windows[0].num_chans / 2) - static_flags = edge - - n_ants = len(ants) - ant_ind = np.arange(n_ants) + rfi_static_flags = np.tile(True,length) + else : + h5 = katfile.open(filename,**kwargs) + length = h5.shape[1] + pickle_file = open(rfi_mask) + rfi_static_flags = pickle.load(pickle_file) + pickle_file.close() + # Now find the edges of the mask + rfi_freqs_width = (856000000.0/rfi_static_flags.shape[0]) + rfi_freqs_min = 856000000.0-rfi_freqs_width/2. # True Start of bin + rfi_freqs_max = rfi_freqs_min*2-rfi_freqs_width/2. # Middle of Max-1 bin + rfi_freqs = np.linspace(rfi_freqs_min,rfi_freqs_max,rfi_static_flags.shape[0]) + rfi_function = get_fit(np.r_[rfi_freqs,rfi_freqs+rfi_freqs_width*0.9999] , np.r_[rfi_static_flags,rfi_static_flags]) + rfi_static_flags = rfi_function(h5.channel_freqs) + #print ( (rfi_static_flags_new-rfi_static_flags)**2).sum() + if verbose: logger.debug(h5.__str__()) + edge = np.tile(True,length) + #edge[slice(211,3896)] = False #Old + edge[slice(int(round(length*0.0515)),int(round(0.9512*length)))] = False### + static_flags = np.logical_or(edge,rfi_static_flags) + + #ants = h5.ants[38:] + if receptor=='': + ants = h5.ants + print(ants) + raw_input("Press Enter to continue...") + else: + ant_names=[a.name for a in h5.ants] + #receptor_name='m053' + receptor_index=ant_names.index(receptor) + print(ant_names) + print('We are now serving customer number: ', receptor_index, ' Receptor: ', ant_names[receptor_index]) + ants = h5.ants[receptor_index:receptor_index+1] + print(ants) + raw_input("Press Enter to continue...") + colour = ['b', 'g', 'r', 'c', 'm', 'y', 'k'] pols = ['v','h'] - diode= 'coupler' - for a,col in zip(ants,colour): - if pdf: - pp = PdfPages(output_dir+'/'+nice_filename+'.'+a.name+'.pdf') - logger.debug("Created output PDF file: %s"%output_dir+'/'+nice_filename+'.'+a.name+'.pdf') - - if not(Ku): - fig1 = plt.figure(2,figsize=(20,5)) - fig2 = plt.figure(1,figsize=(20,5)) - - fig0 = plt.figure(0,figsize=(20,5)) + + for a in ants: + ant = a.name + try: + rx_sn = h5.receivers[ant] + except KeyError: + logger.error('Receiver serial number for antennna %s not found in the H5 file'%ant) + rx_sn = 'l.SN_NOT_FOUND' + band,SN = rx_sn.split('.') + if pdf: + pdf_filename = output_dir+'/'+nice_filename+'.'+rx_sn+'.'+a.name+'.pdf' + pp = PdfPages(pdf_filename) + logger.debug("Created output PDF file: %s"%pdf_filename) + + #fig0 = plt.figure(0,figsize=(20,5)) h5.select() h5.select(ants = a.name,channels=~static_flags) - d = scape.DataSet(h5) - scape.plot_xyz(d,'time','amp',label='Average of the data') - on = h5.sensor['Antennas/'+a.name+'/nd_coupler'] - ts = h5.timestamps - h5.timestamps[0] - plt.plot(ts,np.array(on).astype(float)*4000,'g',label='katdal ND sensor') - plt.title("Timeseries for antenna %s - %s"%(a.name,git_info())) - plt.legend() + observer = h5.ants[0].observer; observer.date = time.gmtime(h5.timestamps.mean())[:6] # katdal resets this date to now()! + fig0 = plot_ts(h5) + Tsys, TAc, Tsys_std = {}, {}, {} + eta_A = {} + Tdiode = {} + nd_temp = {} for pol in pols: logger.debug("Processing: %s%s"%(a.name,pol)) - ant = a.name - ant_num = int(ant[3]) - air_temp = h5.temperature.mean() + Tsys_std[pol] = None if not(Ku): - try: - rx_sn = h5.receivers[ant] - except KeyError: - logger.error('Receiver serial number for antennna %s not found in the H5 file'%ant) - rx_sn = 'SN_NOT_FOUND' diode_filename = '/var/kat/katconfig/user/noise-diode-models/mkat/rx.'+rx_sn+'.'+pol+'.csv' logger.info('Loading noise diode file %s from config'%diode_filename) try: @@ -88,148 +222,78 @@ def read_and_plot_data(filename,output_dir='.',pdf=True,Ku = False,verbose = Fal logger.error("Be sure to reprocess the data once the file is in the config") nd = scape.gaincal.NoiseDiodeModel(freq=[856,1712],temp=[20,20]) - s = h5.spectral_windows[0] - f_c = s.centre_freq #cold data logger.debug('Using off target %s'%target) h5.select(ants=a.name,pol=pol,channels=~static_flags, targets = target,scans='track') freq = h5.channel_freqs - if not(Ku): nd_temp = nd.temperature(freq / 1e6) + cold_data = h5.vis[:].real - on = h5.sensor['Antennas/'+ant+'/nd_coupler'] - n_on = np.tile(False,on.shape[0]) - n_off = np.tile(False,on.shape[0]) - buff = 5 - if not any(on): - logger.critical('No noise diode fired during track of %s'%target) - else: - jumps = (np.diff(on).nonzero()[0] + 1).tolist() - n_on[slice(jumps[0]+buff,jumps[1]-buff)] = True - n_off[slice(jumps[1]+buff,-buff)] = True - - cold_off = n_off - cold_on = n_on + cold_on,cold_off = get_nd_on_off(h5,log=logger) #hot data h5.select(ants=a.name,pol=pol,channels=~static_flags,targets = 'Moon',scans='track') + hot_on,hot_off = get_nd_on_off(h5,log=logger) hot_data = h5.vis[:].real - on = h5.sensor['Antennas/'+ant+'/nd_coupler'] - n_on = np.tile(False,on.shape[0]) - n_off = np.tile(False,on.shape[0]) - if not any(on): - logger.critical('No noise diode fired during track of %s'%target) - else: - jumps = (np.diff(on).nonzero()[0] + 1).tolist() - n_on[slice(jumps[0]+buff,jumps[1]-buff)] = True - n_off[slice(jumps[1]+buff,-buff)] = True - - hot_off = n_off - hot_on = n_on + cold_spec = np.mean(cold_data[cold_off,:,0],0) hot_spec = np.mean(hot_data[hot_off,:,0],0) cold_nd_spec = np.mean(cold_data[cold_on,:,0],0) hot_nd_spec = np.mean(hot_data[hot_on,:,0],0) - if error_bars: - cold_spec_std = np.std(cold_data[cold_off,:,0],0) - hot_spec_std = np.std(hot_data[hot_off,:,0],0) - cold_nd_spec_std = np.std(cold_data[cold_on,:,0],0) - hot_nd_spec_std = np.std(hot_data[hot_on,:,0],0) - if not(Ku): - TAh = hot_spec/(hot_nd_spec - hot_spec) * nd_temp # antenna temperature on the moon (from diode calibration) - TAc = cold_spec/(cold_nd_spec - cold_spec) * nd_temp # antenna temperature on cold sky (from diode calibration) (Tsys) + nd_temp[pol] = nd.temperature(freq / 1e6) + # antenna temperature on the moon (from diode calibration) + TAh = hot_spec/(hot_nd_spec - hot_spec) * nd_temp[pol] + # antenna temperature on cold sky (from diode calibration) (Tsys) + TAc[pol] = cold_spec/(cold_nd_spec - cold_spec) * nd_temp[pol] + print("Mean TAh = %f mean TAc = %f "%(TAh.mean(),TAc[pol].mean())) Y = hot_spec / cold_spec - if error_bars: Y_std = Y * np.sqrt((hot_spec_std/hot_spec)**2 + (cold_spec_std/cold_spec)**2) - D = 13.5 - lam = 3e8/freq + D = 13.5 # Efficiency tables are defined for 13.5 + lam = 299792458./freq HPBW = 1.18 *(lam/D) Om = 1.133 * HPBW**2 # main beam solid angle for a gaussian beam - R = np.radians(0.25) # radius of the moon - Os = np.pi * R**2 # disk source solid angle - _f_MHz, _eff_pct = np.loadtxt("/var/kat/katconfig/user/aperture-efficiency/mkat/ant_eff_L_%s_AsBuilt.csv"%pol.upper(), skiprows=2, delimiter="\t", unpack=True) - eta_A = np.interp(freq,_f_MHz,_eff_pct)/100 # EMSS aperture efficiency - if Ku: eta_A = 0.7 + R = 0.5*Dmoon(observer) # radius of the moon + Os = 2*np.pi*(1-np.cos(R)) # disk source solid angle + _f_MHz, _eff_pct = np.loadtxt("/var/kat/katconfig/user/aperture-efficiency/mkat/ant_eff_%s_%s_AsBuilt.csv"%(band.upper(),pol.upper()), skiprows=2, delimiter="\t", unpack=True) + eta_A[pol] = np.interp(freq,_f_MHz,_eff_pct)/100. # EMSS aperture efficiency + if Ku: eta_A[pol] = 0.7 Ag = np.pi* (D/2)**2 # Antenna geometric area - Ae = eta_A * Ag # Effective aperture + Ae = eta_A[pol] * Ag # Effective aperture x = 2*R/HPBW # ratio of source to beam K = ((x/1.2)**2) / (1-np.exp(-((x/1.2)**2))) # correction factor for disk source from Baars 1973 - TA_moon = 225 * (Os*Ae/(lam**2)) * (1/K) # contribution from the moon (disk of constant brightness temp) - if error_bars: Thot_std = 2.25 + TA_moon = 225 * (Os/Om) * (1/K) # contribution from the moon (disk of constant brightness temp) gamma = 1.0 - if error_bars: gamma_std = 0.01 - Tsys = gamma * (TA_moon)/(Y-gamma) # Tsys from y-method ... compare with diode TAc - if error_bars: Tsys_std = Tsys * np.sqrt((Thot_std/Thot)**2 + (Y_std/Y)**2 + (gamma_std/gamma)**2) + Thot = TA_moon + Tcold = 0 + Tsys[pol] = gamma * (Thot-Tcold)/(Y-gamma) # Tsys from y-method ... compare with diode TAc + if error_bars: + cold_spec_std = np.std(cold_data[cold_off,:,0],0) + hot_spec_std = np.std(hot_data[hot_off,:,0],0) + cold_nd_spec_std = np.std(cold_data[cold_on,:,0],0) + hot_nd_spec_std = np.std(hot_data[hot_on,:,0],0) + Y_std = Y * np.sqrt((hot_spec_std/hot_spec)**2 + (cold_spec_std/cold_spec)**2) + Thot_std = 2.25 + gamma_std = 0.01 + # This is not definded + raise NotImplementedError("The factor Thot has not been defined ") + Tsys_std[pol] = Tsys[pol] * np.sqrt((Thot_std/Thot)**2 + (Y_std/Y)**2 + (gamma_std/gamma)**2) + else : + Tsys_std[pol] = None if not(Ku): Ydiode = hot_nd_spec / hot_spec - Tdiode = (TA_moon + Tsys)*(Ydiode/gamma-1) - - p = 1 if pol == 'v' else 2 - if not(Ku): - plt.figure(2) - plt.subplot(1,2,p) - plt.ylim(0,50) - plt.ylabel('T_ND [K]') - plt.xlim(900,1670) - plt.xlabel('f [MHz]') - if p ==ant_num * 2-1: plt.ylabel(ant) - plt.axhspan(14, 35, facecolor='g', alpha=0.5) - plt.plot(freq/1e6,Tdiode,'b.',label='Measurement: Y-method') - if write_nd: - outfilename = diode_filename.split('/')[-1] - outfile = open(outfilename, 'w') - outfile.write('#Data from %s\n# Frequency [Hz], Temperature [K]\n'%file_base) - # Write CSV part of file - outfile.write(''.join(['%s, %s\n' % (entry[0], entry[1]) for entry in zip(freq,medfilt(Tdiode))])) - outfile.close() - logger.info('Noise temp data written to file %s'%outfilename) - plt.plot(freq/1e6,nd_temp,'k.',label='Model: EMSS') - plt.grid() - plt.legend() - - plt.figure(1) - plt.subplot(1,2,p) - if not(Ku): plt.ylim(15,50) - plt.ylabel('Tsys/eta_A [K]') - if not(Ku): plt.xlim(900,1670) - plt.xlabel('f [MHz]') - if p == ant_num * 2 -1: plt.ylabel(ant) - if error_bars: plt.errorbar(freq/1e6,Tsys,Tsys_std,color = 'b',linestyle = '.',label='Measurement') - plt.plot(freq/1e6,Tsys/eta_A,'b.',label='Measurement: Y-method') - if not(Ku): plt.plot(freq/1e6,TAc/eta_A,'c.',label='Measurement: ND calibration') - plt.axhline(np.mean(Tsys/eta_A),linewidth=2,color='k',label='Mean: Y-method') - spec_Tsys_eta = 0*freq - spec_Tsys_eta[freq<1420e6] = 42 # [R.T.P095] == 220 - spec_Tsys_eta[freq>=1420e6] = 46 # [R.T.P.096] == 200 - if not(Ku): plt.plot(freq/1e6, spec_Tsys_eta,'r',linewidth=2,label='PDR Spec') - if not(Ku): plt.plot(freq/1e6,np.interp(freq/1e6,[900,1670],[(64*Ag)/275.0,(64*Ag)/410.0]),'g',linewidth=2,label="275-410 m^2/K at Receivers CDR") - - plt.grid() - plt.legend(loc=2,fontsize=12) + Tdiode_h = (Tsys[pol]-Tcold+Thot)*(Ydiode/gamma-1) # Tsys as computed above includes Tcold + Ydiode = cold_nd_spec / cold_spec + Tdiode_c = Tsys[pol]*(Ydiode/gamma-1) + Tdiode[pol] = (Tdiode_h+Tdiode_c)/2. # Average two equally valid, independent results + if write_nd: + outfilename = save_ND(diode_filename,file_base,freq,Tdiode[pol] ) + logger.info('Noise temp data written to file %s'%outfilename) - if not(Ku): - plt.figure(2) - plt.subplot(1,2,1) - ax = plt.gca() - ax.text(0.95, 0.01,git_info(), horizontalalignment='right',fontsize=10,transform=ax.transAxes) - plt.title('%s Coupler Diode: V pol: %s'%(ant,file_base)) - plt.subplot(1,2,2) - ax = plt.gca() - ax.text(0.95, 0.01,git_info(), horizontalalignment='right',fontsize=10,transform=ax.transAxes) - plt.title('%s Coupler Diode: H pol: %s'%(ant,file_base)) - - plt.figure(1) - plt.subplot(1,2,1) - ax = plt.gca() - ax.text(0.95, 0.01,git_info(), horizontalalignment='right',fontsize=10,transform=ax.transAxes) - plt.title('%s Tsys/eta_A: V pol: %s'%(ant,file_base)) - plt.subplot(1,2,2) - ax = plt.gca() - ax.text(0.95, 0.01,git_info(), horizontalalignment='right',fontsize=10,transform=ax.transAxes) - plt.title('%s Tsys/eta_A: H pol: %s'%(ant,file_base)) + fig2 = plot_nd(freq,Tdiode,nd_temp,ant = ant, file_base=file_base) + fig1 = plot_Tsys_eta_A(freq,Tsys,eta_A,TAc,Tsys_std=Tsys_std,ant = ant, file_base=file_base,Ku=Ku) if pdf: if not(Ku): - fig1.savefig(pp,format='pdf') - fig2.savefig(pp,format='pdf') + fig2.savefig(pp,format='pdf') + fig1.savefig(pp,format='pdf') fig0.savefig(pp,format='pdf') pp.close() # close the pdf file plt.close("all") @@ -240,13 +304,21 @@ def read_and_plot_data(filename,output_dir='.',pdf=True,Ku = False,verbose = Fal # test main method for the library if __name__ == "__main__": #test the method with a know file - filename = '/var/kat/archive/data/RTS/telescope_products/2015/04/22/1429706275.h5' + filename = '/var/kat/archive/data/RTS/telescope_products/2016/01/19/1453216690.h5' + #filename = '/var/kat/archive2/data/MeerKATAR1/telescope_products/2016/07/21/1469134098.h5' pdf=True Ku=False verbose=False out = '.' error_bars = False target = 'off1' - write_nd = False - print('Performing test run with: ' + filename) + write_nd = True + print 'Performing test run with: ' + filename read_and_plot_data(filename,out,pdf,Ku,verbose,error_bars,target,write_nd) +# Output checsums of the noisediode files for 1453216690.h5s +#md5sum rx.* +#7545907bbe621f4e5937de7536be21f3 rx.l.4002.h.csv +#532332722c8bca9714f4e2d97978ccd6 rx.l.4002.v.csv +#66678871aeaadc4178e4420bef9861aa rx.l.4.h.csv +#af113bd4ab10bb61386521f61eb4d716 rx.l.4.v.csv +