import numpy as np import matplotlib.pyplot as plt from matplotlib.pyplot import rcParams from matplotlib import cm import scipy.ndimage import os # BeamCrossPlots.py # # Written by Chat Hull, Jul 2015 # Final updates for data release: Apr 2020 ##### HOW TO RUN THE CODE IN CASA (VERSION <= 5.6.1) # (1) Place this script in the directory where your data subdirectory is located (see beamCrossRedux.py) # (2) Check that the *peak*.txt files were successfully copied into the directory where this script is by beamCrossRedux.py # Then: # execfile('../beamCrossPlots.py') # beamPlots(B3_3) # Above, replace B3_3 with whichever trial whose plots you want to create # Details for different trials B3_3 = { 'source' : '3c279', 'refant' : 'DV21', 'datasets' : ['uid___A002_X9b98ec_X94d.ms', 'uid___A002_X9b98ec_Xb98.ms', 'uid___A002_X9b98ec_Xde3.ms'], 'size' : '11x11', 'imsize' : [300], 'band' : '3', 'dir' : 'B3-3', 'freq_act' : 97.47895e9, 'freq_nom' : 86.24343e9, } B5_1 = { 'source' : '3c279', 'refant' : 'DV21', 'datasets' : ['uid___A002_Xba6edb_X31db.ms', 'uid___A002_Xba6edb_X3988.ms', 'uid___A002_Xba6edb_X3d57.ms'], 'date' : '10nov16', 'size' : '11x11', 'imsize' : [300], 'band' : '5', 'dir' : 'B5-1', 'freq_act' : 183.26111e9, 'freq_nom' : 177.26111e9, } B6_1 = { 'source' : '3c279', 'refant' : 'DA59', 'datasets' : ['uid___A002_Xa018c4_X743.ms', 'uid___A002_Xa018c4_Xa3d.ms', 'uid___A002_Xa018c4_Xdc7.ms'], 'size' : '11x11', 'imsize' : [300], 'band' : '6', 'dir' : 'B6-1', 'freq_act' : 233.00000e9, 'freq_nom' : 230.538e9, } B7_1 = { 'source' : '3c279', 'refant' : 'DA59', 'datasets' : ['uid___A002_Xa018c4_X3b92.ms', 'uid___A002_Xa018c4_X3de3.ms', 'uid___A002_Xa018c4_X402a.ms'], 'size' : '11x11', 'imsize' : [300], 'band' : '7', 'dir' : 'B7-1', 'freq_act' : 343.47895e9, 'freq_nom' : 345.796e9, } def beamPlots(trial) : # Check directory if os.path.abspath('.').split('/')[-1] != trial['dir'] : raise Exception("You're not in the correct directory!") # Use latex rcParams['text.usetex'] = True ##### Beam definition # beam = 1.22*(clight/frequency)/antennaDiameter # 1.22 * lambda / D is the distance to the first minimum (1/2 of the full beam width) clight = 299792458. antennaDiameter = 12. ##### THIS IS THE THEORETICAL BEAM USED BY THE POINTING OFFSETS ##### ##### Calculate beam (radius to first null, using meters and m/s units) beam = 1.22 * (clight / trial['freq_nom']) / antennaDiameter * 206265 ##### THIS IS THE ACTUAL BEAM ##### beam_act = 1.22 * (clight / trial['freq_act']) / antennaDiameter * 206265 ##### Offset distance # Each step is (1/5) * (1/2) * (1/2) = 1/20 of the full beam width in Az and El # -Or (1/5) * (1/2) * beam [the variable] ##### Delta is calculated with "beam"...but data are subject to beam_act delta = (1/2.) * (1/5.) * beam # Grid of points with appropriate extents if trial['size'] == '11x11' : radius = 5 az,el = np.mgrid[-radius*delta:(radius+1)*delta:delta, -radius*delta:(radius+1)*delta:delta] # Close previous plot windows plt.close('all') # Variables I = [] Q = [] U = [] V = [] pos = [] # Get values from text files I_file = '../{0}/Ipeak.B{1}.txt'.format(trial['dir'],trial['band']) Q_file = '../{0}/Qpeak.B{1}.txt'.format(trial['dir'],trial['band']) U_file = '../{0}/Upeak.B{1}.txt'.format(trial['dir'],trial['band']) V_file = '../{0}/Vpeak.B{1}.txt'.format(trial['dir'],trial['band']) for i in [I_file, Q_file, U_file, V_file] : f = open(i,'r') lines = f.readlines() f.close() # Get list of I values (units = Jy, so divide by 1000) # Except for corner cases where units are in Jy/beam or uJy/beam for line in lines : if i == I_file : if line.split()[6] == 'uJy/beam' : I.append(float(line.split()[3])/1.e6) elif line.split()[6] == 'mJy/beam' : I.append(float(line.split()[3])/1000.) elif line.split()[6] == 'Jy/beam' : I.append(float(line.split()[3])/1.) elif i == Q_file : if line.split()[6] == 'uJy/beam' : Q.append(float(line.split()[3])/1.e6) elif line.split()[6] == 'mJy/beam' : Q.append(float(line.split()[3])/1000.) elif line.split()[6] == 'Jy/beam' : Q.append(float(line.split()[3])/1.) elif i == U_file : if line.split()[6] == 'uJy/beam' : U.append(float(line.split()[3])/1.e6) elif line.split()[6] == 'mJy/beam' : U.append(float(line.split()[3])/1000.) elif line.split()[6] == 'Jy/beam' : U.append(float(line.split()[3])/1.) elif i == V_file : if line.split()[6] == 'uJy/beam' : V.append(float(line.split()[3])/1.e6) elif line.split()[6] == 'mJy/beam' : V.append(float(line.split()[3])/1000.) elif line.split()[6] == 'Jy/beam' : V.append(float(line.split()[3])/1.) # Grid sizes if trial['size'] == '11x11' : size = 11 center = 5 # Make list of positions pos = list(np.arange(size**2)) ##### Reshape grids Igrid = np.array(I).reshape(size,size).T Qgrid = np.array(Q).reshape(size,size).T Ugrid = np.array(U).reshape(size,size).T Vgrid = np.array(V).reshape(size,size).T posgrid = np.array(pos).reshape(size,size).T # .T makes observations move from upper left to lower right, in standard English-reading style # However, the observations were taken from the upper-left, down vertically, in columns moving to the right (i.e., 1 is upper-left, 11 is lower-left, 111 is upper-right, 121 is lower-right). Note that the orientation plot in ALMA SYS#225 report is incorrect. Ipeak = Igrid[center,center] print 'I peak = {0}'.format(Ipeak) Qpeak = Qgrid[center,center] print 'Q peak = {0}'.format(Qpeak) Upeak = Ugrid[center,center] print 'U peak = {0}'.format(Upeak) Vpeak = Vgrid[center,center] print 'V peak = {0}'.format(Vpeak) # FWHM circle FWHM = beam_act r = 0.5*FWHM x = r * np.sin( np.arange(0, 2*np.pi, 0.01) ) y = r * np.cos( np.arange(0, 2*np.pi, 0.01) ) # Extent of plots extent = [-5*delta,5*delta,-5*delta,5*delta] # Color map cmap_div = cm.PuOr_r cmap_seq = cm.YlOrBr # Lighter color for contours to enhance visibility lighter = '#A8A8A8' ##### Band-based choices if trial == B3_3 : contour_levels_PA = np.arange(-4, 6, 1) contour_colors_PA = [lighter,lighter,'k','k','k','k','k','k',lighter,lighter] contour_linestyles_PA = ['dashed','dashed','dashed','dashed','solid','solid','solid','solid','solid','solid'] contour_levels_frac = np.array([-0.004, -0.001, 0.001, 0.004]) contour_colors_frac = ['k','k','k','k'] contour_linestyles_frac = ['dashed','dashed','solid','solid'] contour_levels_V = np.arange(-0.01, 0.015, 0.005) contour_colors_V = [lighter,'k','k','k',lighter] contour_linestyles_V = ['dashed','dashed','solid','solid','solid'] axis_ticks = np.arange(-30,40,10) if trial == B5_1 : contour_levels_PA = np.arange(-1.5, 2.5, 0.5) contour_colors_PA = [lighter,'k','k','k','k','k',lighter,lighter] contour_linestyles_PA = ['dashed','dashed','dashed','solid','solid','solid','solid','solid'] contour_levels_frac = np.array([-0.004, -0.001, 0.001, 0.004]) contour_colors_frac = [lighter,'k','k',lighter] contour_linestyles_frac = ['dashed','dashed','solid','solid'] contour_levels_V = np.arange(-0.015, 0.02, 0.005) contour_colors_V = [lighter,'k','k','k','k','k',lighter] contour_linestyles_V = ['dashed','dashed','dashed','solid','solid','solid','solid'] axis_ticks = np.arange(-15,20,5) elif trial == B6_1 : contour_levels_PA = np.arange(-1.0, 2.0, 0.5) contour_colors_PA = ['k','k','k','k','k',lighter] contour_linestyles_PA = ['dashed','dashed','solid','solid','solid','solid'] contour_levels_frac = np.array([-0.004, -0.001, 0.001, 0.004]) contour_colors_frac = [lighter,'k','k',lighter] contour_linestyles_frac = ['dashed','dashed','solid','solid'] contour_levels_V = np.arange(-0.015, 0.02, 0.005) contour_colors_V = [lighter,'k','k','k','k','k',lighter] contour_linestyles_V = ['dashed','dashed','dashed','solid','solid','solid','solid'] axis_ticks = np.arange(-10,15,5) elif trial == B7_1 : contour_levels_PA = np.arange(-6, 7, 1) contour_colors_PA = [lighter,lighter,lighter,'k','k','k','k','k','k','k','k',lighter,lighter] contour_linestyles_PA = ['dashed','dashed','dashed','dashed','dashed','dashed','solid','solid','solid','solid','solid','solid','solid'] contour_levels_frac = np.array([-0.004, -0.001, 0.001, 0.004]) contour_colors_frac = ['k','k','k','k'] contour_linestyles_frac = ['dashed','dashed','solid','solid'] contour_levels_V = np.arange(-0.02, 0.025, 0.005) contour_colors_V = [lighter,lighter,'k','k','k','k','k',lighter,lighter] contour_linestyles_V = ['dashed','dashed','dashed','dashed','solid','solid','solid','solid','solid'] axis_ticks = np.arange(-5,10,5) ##### PLOT DATA AND ERRORS OUT TO ROUGHLY THE FWHM ##### Q error plt.figure(1) plt.xlabel(r'$\rm Az\; offset\; (^{\prime\prime})$', fontsize=16) plt.xticks(axis_ticks) plt.ylabel(r'$\rm El\; offset\; (^{\prime\prime})$', fontsize=16) plt.yticks(axis_ticks) # To fix B7 extent issue in plots with contours if trial == B7_1 : plt.xlim(-9.12,9.12) plt.ylim(-9.12,9.12) plt.tick_params(direction='out', labelsize=14) plt.plot(x, y, 'k:', lw=1) plt.plot(x/3, y/3, 'k:', lw=1) # Normalize by I and subtract central value values = (Qgrid/Igrid) - Qpeak # Make max and min symmetric around 0 vmax = 0.75*max([ abs(values.min()),abs(values.max()) ]) vmin = -vmax plt.imshow(values, interpolation='none', extent=extent, origin='upper', vmin=vmin, vmax=vmax, cmap=cmap_div) cbar = plt.colorbar() cbar.set_label(r'$\rm Stokes\; {\it Q}\; error$ $(Q_{\rm off}/I_{\rm off} - Q_{\rm on})$', fontsize=16) cbar.ax.tick_params(direction='out', labelsize=14) # Smooth data for contours values = scipy.ndimage.zoom(values,10) plt.subplots_adjust(bottom=0.15) plt.savefig('Q_error_band_{0}.png'.format(trial['dir']), dpi=450) ##### U error plt.figure(2) plt.xlabel(r'$\rm Az\; offset\; (^{\prime\prime})$', fontsize=16) plt.xticks(axis_ticks) plt.ylabel(r'$\rm El\; offset\; (^{\prime\prime})$', fontsize=16) plt.yticks(axis_ticks) plt.tick_params(direction='out', labelsize=14) plt.plot(x, y, 'k:', lw=1) plt.plot(x/3, y/3, 'k:', lw=1) # Normalize by I and subtract central value values = (Ugrid/Igrid) - Upeak # Make max and min symmetric around 0 vmax = 0.75*max([ abs(values.min()),abs(values.max()) ]) vmin = -vmax plt.imshow(values, interpolation='none', extent=extent, origin='upper', vmin=vmin, vmax=vmax, cmap=cmap_div) cbar = plt.colorbar() cbar.set_label(r'$\rm Stokes\; {\it U}\; error$ $(U_{\rm off}/I_{\rm off} - U_{\rm on})$', fontsize=16) cbar.ax.tick_params(direction='out', labelsize=14) # Smooth data for contours values = scipy.ndimage.zoom(values,10) plt.subplots_adjust(bottom=0.15) plt.savefig('U_error_band_{0}.png'.format(trial['dir']), dpi=450) ##### PA error plt.figure(3) plt.xlabel(r'$\rm Az\; offset\; (^{\prime\prime})$', fontsize=16) plt.xticks(axis_ticks) plt.ylabel(r'$\rm El\; offset\; (^{\prime\prime})$', fontsize=16) plt.yticks(axis_ticks) # To fix B7 extent issue in plots with contours if trial == B7_1 : plt.xlim(-9.12,9.12) plt.ylim(-9.12,9.12) plt.tick_params(direction='out', labelsize=14) plt.plot(x, y, 'k:', lw=1) plt.plot(x/3, y/3, 'k:', lw=1) values = np.degrees(0.5*np.arctan2(Ugrid,Qgrid) - 0.5*np.arctan2(Upeak,Qpeak)) # Make max and min symmetric around 0 vmax = 0.75*max([ abs(values.min()),abs(values.max()) ]) vmin = -vmax plt.imshow(values, interpolation='none', extent=extent, origin='upper', vmin=vmin, vmax=vmax, cmap=cmap_div) cbar = plt.colorbar() cbar.set_label(r'$\chi\; {\rm error}\; (\chi_{\rm off} - \chi_{\rm on})$', fontsize=16) cbar.ax.tick_params(direction='out', labelsize=14) # Smooth data for contours values = scipy.ndimage.zoom(values,10) C = plt.contour( values, contour_levels_PA, linewidths=1, extent=extent, origin='upper', colors=contour_colors_PA, linestyles=contour_linestyles_PA) plt.clabel(C, inline=1, inline_spacing=8, fontsize=13, fmt=r'$%1.1f$') plt.subplots_adjust(bottom=0.15) plt.savefig('PA_error_band_{0}.png'.format(trial['dir']), dpi=450) ##### Pol frac error plt.figure(4) plt.xlabel(r'$\rm Az\; offset\; (^{\prime\prime})$', fontsize=16) plt.xticks(axis_ticks) plt.ylabel(r'$\rm El\; offset\; (^{\prime\prime})$', fontsize=16) plt.yticks(axis_ticks) # To fix B7 extent issue in plots with contours if trial == B7_1 : plt.xlim(-9.12,9.12) plt.ylim(-9.12,9.12) plt.tick_params(direction='out', labelsize=14) plt.plot(x, y, 'k:', lw=1) plt.plot(x/3, y/3, 'k:', lw=1) values = np.sqrt(Qgrid**2+Ugrid**2)/Igrid - np.sqrt(Qpeak**2+Upeak**2)/Ipeak # Make max and min symmetric around 0 vmax = 0.75*max([ abs(values.min()),abs(values.max()) ]) vmin = -vmax plt.imshow(values, interpolation='none', extent=extent, origin='upper', vmin=vmin, vmax=vmax, cmap=cmap_div) cbar = plt.colorbar() cbar.set_label(r'$P_{\rm frac}\; {\rm error}\; (P_{\rm frac,off} - P_{\rm frac,on})$', fontsize=16) cbar.ax.tick_params(direction='out', labelsize=14) # Smooth data for contours values = scipy.ndimage.zoom(values,10) C = plt.contour( values, contour_levels_frac, linewidths=1, extent=extent, origin='upper', colors=contour_colors_frac, linestyles=contour_linestyles_frac) plt.clabel(C, inline=1, inline_spacing=8, fontsize=13, fmt=r'$%1.3f$') plt.subplots_adjust(bottom=0.15) plt.savefig('frac_error_band_{0}.png'.format(trial['dir']), dpi=450) ##### V data plt.figure(5) plt.xlabel(r'$\rm Az\; offset\; (^{\prime\prime})$', fontsize=16) plt.xticks(axis_ticks) plt.ylabel(r'$\rm El\; offset\; (^{\prime\prime})$', fontsize=16) plt.yticks(axis_ticks) # To fix B7 extent issue in plots with contours if trial == B7_1 : plt.xlim(-9.12,9.12) plt.ylim(-9.12,9.12) plt.tick_params(direction='out', labelsize=14) plt.plot(x, y, 'k:', lw=1) values = Vgrid # Make max and min symmetric around 0 vmax = max([ abs(values.min()),abs(values.max()) ]) vmin = -vmax plt.imshow(values, interpolation='none', extent=extent, origin='upper', vmin=vmin, vmax=vmax, cmap=cmap_div) cbar = plt.colorbar() cbar.set_label(r'${\rm Stokes}\; V\; {\rm squint}\; (V_{\rm off}/I_{\rm on},\; {\rm Band\; %s})$'%trial['band'], fontsize=16) cbar.ax.tick_params(direction='out', labelsize=14) # Smooth data for contours values = scipy.ndimage.zoom(values,10) C = plt.contour( values, contour_levels_V, linewidths=1, extent=extent, origin='upper', colors=contour_colors_V, linestyles=contour_linestyles_V) plt.clabel(C, inline=1, inline_spacing=8, fontsize=13, fmt=r'$%1.3f$') plt.subplots_adjust(bottom=0.15) plt.savefig('V_band_{0}.png'.format(trial['dir']), dpi=450)