Sample script to plot co-pol and x-pol beam patterns from .cut files. ACT
0001 function plot_example_beam 0002 % Sample script to plot co-pol and x-pol beam patterns from .cut files. 0003 % ACT 0004 0005 % Data from Christian is in the form 0006 % Line 1: Field data in cuts 0007 % Line 2: minTheta Theta_step n_points phi_angle 3 1 2 0008 0009 % Number of frequencies and there defintion 0010 %n_cuts = 2; % Get this from the file instead 0011 0012 % Specification for the .cut files used in figure 8 of the optics paper 0013 % BUT here weights set to 1 - optics paper uses spectral index of CasA 0014 % to allow direct comparison with our observation presented 0015 ifreqs = 3; 0016 freqs = {'4.5';'5.0';'5.5'}; % Frequency in GHz 0017 weights = ones(1,ifreqs); % Use this to set band 0018 alpha = 0; % spectral index - other code uses this to scale response to sources with frequency 0019 %alpha = -0.351; % Hafez et al. 1.8-33GHz give -0.351 (Baars give -0.299 1-35GHz) 0020 freq = [4.5 5.0 5.5]; 0021 weights = weights .*freq(1:ifreqs).^alpha; 0022 0023 % Specs for different cut files 0024 % % ifreqs = 11; 0025 % % freqs = {'4.5';'4.6';'4.7';'4.8';'4.9';'5.0';'5.1';'5.2';'5.3';'5.4';'5.5'}; 0026 % % weights = ones(1,ifreqs); 0027 % % alpha=0; 0028 % % %make a weighting based on spectrum of TauA 0029 % % %alpha = -0.351; % Hafez et al. 1.8-33GHz give -0.351 (Baars give -0.299 1-35GHz) 0030 % % freq = [4.5,4.6,4.7,4.8,4.9,5.0,5.1,5.2,5.3,5.4,5.5]; 0031 0032 0033 % Path to where you .cut files are 0034 simpath = './sim_beams/'; 0035 0036 % Now read in each frequecny file and get the relevant data 0037 for f=1:ifreqs 0038 freq(f) 0039 simfile =([simpath,num2str(freq(f)),'GHz_with_cryo_cone_5mm_up.cut']) ;% Fig 8 optics paper files 0040 %simfile =([simpath,'sph_1_',num2str(freqs{f}),'.txt']) 0041 %simfile =([path,'full_sim_coneandtunnels_19cuts.txt']); 0042 %simfile =([path,'with cryo cone absorb subrefl no hole 19cuts lin pol 85dB 4.5GHz.cut']) 0043 0044 0045 % Grab the number of cuts in the file 0046 0047 [status, n_cuts] = system(['grep -c "Field" ',simfile]); 0048 clear status; 0049 n_cuts = str2num(n_cuts); 0050 0051 % First set up the Theta Angles samples 0052 0053 start_theta = dlmread(simfile,'',[1 0 1 0]); % degrees 0054 theta_step = dlmread(simfile,'',[1 1 1 1]); % degrees 0055 no_steps = dlmread(simfile,'',[1 2 1 2]); % degrees 0056 sim_angles = linspace(start_theta, -1*start_theta, no_steps); 0057 0058 0059 % Now read in each phi cut 0060 0061 for i=1:n_cuts 0062 phi_position = [ (1 + (i-1)*(no_steps+2)) 3 (1 + (i-1)*(no_steps+2)) 3 ] ; 0063 data_start = 1+((1 + (i-1)*(no_steps+2))) ; 0064 data_stop = data_start + no_steps - 1; 0065 data_position = [data_start 0 data_stop 3]; 0066 phi(i) = dlmread(simfile,'', phi_position); % degrees 0067 sim_data{i} = dlmread(simfile,'',data_position ); % 4 cols: co-polar (Re Im) x-pol(Re Im) 0068 end 0069 0070 % For the x-pol beam patterns we just want data from the 45deg cut 0071 % Take : 0072 % 0073 % 45Phi_xPol = (Re^2 + Im^2) 0074 % 0075 0076 Phi45 = find(phi == 45); 0077 0078 if((Phi45>0)) 0079 cbass_real_Phi45{f} = sim_data{Phi45}(:,3); % second pair of columns is the xpol 0080 cbass_imag_Phi45{f} = sim_data{Phi45}(:,4); 0081 0082 0083 % Calulate powers for real and imaginary, then take average of ZP and P for 0084 % each frequency 0085 0086 power_P{f} = (cbass_real_Phi45{f}.^2 + cbass_imag_Phi45{f}.^2); 0087 0088 sim_beam_cell_xpol{f} = ((power_P{f})); 0089 0090 % Plot the beams 0091 % Far-out beam 0--> 180 deg 0092 figure(1) 0093 %subplot(ifreqs+1,1,f) 0094 plot(sim_angles, 10.*log10(sim_beam_cell_xpol{f}./1.0)) 0095 hold all 0096 xlim([min(sim_angles) max(sim_angles)]) 0097 ylim([-100 0]) 0098 xlabel('Degrees','Fontsize',24,'Fontname','Helvetica') 0099 ylabel('dBi','Fontsize',24,'Fontname','Helvetica') 0100 title(['Cross-polar beam']) 0101 legend([cell2mat(freqs)]) 0102 else 0103 disp('Cut file does not contain cross-polar data') 0104 end 0105 %figure 0106 0107 %% Co-pol 0108 % Example of what to do if you want e.g the beam pattern from the 0 deg 0109 % and 90deg cut and get what CBASS intensity mesures 0110 % % For the beam patterns we want to combine data from the 0deg and 90deg phi 0111 % % cuts 0112 % % Take mean of : 0113 % % zeroPhi_coPol = (Re^2 + Im^2) 0114 % % 90Phi_coPol = (Re^2 + Im^2) 0115 % % 0116 % % To get what CBASS measures: 0.5*(Ex^2 + Ey^2) 0117 % 0118 0119 zeroPhi = find(phi == 0); 0120 Phi90 = find(phi == 90); 0121 0122 cbass_real_zeroPhi{f} = sim_data{zeroPhi}(:,1); 0123 cbass_imag_zeroPhi{f} = sim_data{zeroPhi}(:,2); 0124 cbass_real_Phi90{f} = sim_data{Phi90}(:,1); 0125 cbass_imag_Phi90{f} = sim_data{Phi90}(:,2); 0126 0127 % Calulate powers for real and imaginary, then take average of ZP and P for 0128 % each frequency 0129 0130 power_ZP{f} = (cbass_real_zeroPhi{f}.^2 + cbass_imag_zeroPhi{f}.^2); 0131 power_P{f} = (cbass_real_Phi90{f}.^2 + cbass_imag_Phi90{f}.^2); 0132 0133 sim_beam_cell_copol{f} = ((power_ZP{f} + power_P{f})./2); 0134 0135 % Plot the beams 0136 % Far-out beam 0--> 180 deg 0137 figure(2) 0138 %subplot(1,ifreqs,f) 0139 plot(sim_angles, 10.*log10(sim_beam_cell_copol{f}./1.0))% Plot in dBi 0140 %plot(sim_angles, 10.*log10(sim_beam_cell_copol{f}./max(sim_beam_cell_copol{f})))%Plot in dB 0141 hold all 0142 grid on 0143 xlim([min(sim_angles) max(sim_angles)]) 0144 ylim([-50 50]) 0145 xlabel('Degrees','Fontsize',24,'Fontname','Helvetica') 0146 ylabel('dBi','Fontsize',24,'Fontname','Helvetica') 0147 title(['Co-polar beam']) 0148 legend([cell2mat(freqs)]) 0149 end 0150 0151 %% 0152 % Calculate weighted mean --> do on real and imaginary separatly and then form powers 0153 % need to first convert cells to matrices, then get mean across frequencies 0154 % at each angle and then form powers as before 0155 %(probably a nicer way of doing this using cunnign code!) 0156 0157 % Co-polar 0158 0159 mat_cbass_real_zeroPhi = cell2mat(cbass_real_zeroPhi); 0160 mat_cbass_imag_zeroPhi = cell2mat(cbass_imag_zeroPhi); 0161 mat_cbass_real_Phi90 = cell2mat(cbass_real_Phi90); 0162 mat_cbass_imag_Phi90 = cell2mat(cbass_imag_Phi90); 0163 0164 % Multiply data by weights 0165 0166 for i=1:length(mat_cbass_real_zeroPhi) 0167 mat_cbass_real_zeroPhi(i,:)= mat_cbass_real_zeroPhi(i,:).*weights; 0168 mat_cbass_imag_zeroPhi(i,:)=mat_cbass_imag_zeroPhi(i,:).*weights; 0169 mat_cbass_real_Phi90(i,:)=mat_cbass_real_Phi90(i,:).*weights; 0170 mat_cbass_imag_Phi90(i,:)=mat_cbass_imag_Phi90(i,:).*weights; 0171 end 0172 0173 mean_cbass_real_zeroPhi = sum(mat_cbass_real_zeroPhi,2)/sum(weights); 0174 mean_cbass_imag_zeroPhi = sum(mat_cbass_imag_zeroPhi,2)/sum(weights); 0175 mean_cbass_real_Phi90 = sum(mat_cbass_real_Phi90,2)/sum(weights); 0176 mean_cbass_imag_Phi90 = sum(mat_cbass_imag_Phi90,2)/sum(weights); 0177 0178 %calculate power for each Phi 0179 mean_power_ZP = (mean_cbass_real_zeroPhi.^2 + mean_cbass_imag_zeroPhi.^2); 0180 mean_power_P = (mean_cbass_real_Phi90.^2 + mean_cbass_imag_Phi90.^2); 0181 0182 mean_copol_power = ((mean_power_ZP + mean_power_P)./2); 0183 0184 figure(3) 0185 hold all 0186 plot(sim_angles', 10.*log10( (mean_copol_power)./max(mean_copol_power)),'g')%Plot in dBm (normalised to max of beam) 0187 %plot(sim_beam_all(:,1), 10.*log10(sim_beam_all(:,2)./max(sim_beam_all(:,2))))%Plot in dB 0188 grid on 0189 xlabel('Degrees','Fontsize',24,'Fontname','Helvetica') 0190 ylabel('dBm','Fontsize',24,'Fontname','Helvetica') 0191 title(['Mean Co-polar beam']) 0192 0193 0194 % Cross-polar 0195 0196 if(Phi45>0) 0197 mat_cbass_real_Phi45 = cell2mat(cbass_real_Phi45); 0198 mat_cbass_imag_Phi45 = cell2mat(cbass_imag_Phi45); 0199 0200 % Multiply data by weights 0201 0202 for i=1:length(mat_cbass_real_zeroPhi) 0203 mat_cbass_real_Phi45(i,:)=mat_cbass_real_Phi45(i,:).*weights; 0204 mat_cbass_imag_Phi45(i,:)=mat_cbass_imag_Phi45(i,:).*weights; 0205 end 0206 0207 0208 mean_cbass_real_Phi45 = sum(mat_cbass_real_Phi45,2)/sum(weights); 0209 mean_cbass_imag_Phi45 = sum(mat_cbass_imag_Phi45,2)/sum(weights); 0210 0211 %calculate power for each Phi 0212 0213 mean_xpol_power = (mean_cbass_real_Phi45.^2 + mean_cbass_imag_Phi45.^2); 0214 0215 0216 figure(4) 0217 hold all 0218 plot(sim_angles', 10.*log10( (mean_xpol_power)./max(mean_xpol_power)),'g')%Plot in dB (normalised to max of beam) 0219 grid on 0220 xlabel('Degrees','Fontsize',24,'Fontname','Helvetica') 0221 ylabel('dB','Fontsize',24,'Fontname','Helvetica') 0222 title(['Mean Cross-polar beam']) 0223 end 0224 0225 0226 %% 0227 %% The code below calculates the mean beam directly by averagin gthe power calculated earlier 0228 % i.e doesnt do real and imag separately 0229 % NEEDS CHECKING and cf above... ACT 0230 % - probably should be weighting real and imaginary separately then forming power 0231 % % Now to plot combined beam over all frequencies 0232 % % Calculate mean over frequency 0233 % % Done simply by summing over frequency and dividing by sum of weights 0234 % % 0235 % % Made generic using weights - code a bit messy here ... tidy up some time. 0236 % 0237 % 0238 % % % Check for 3 frequencies - genericusing weights below 0239 % % % sim_beam_sum =(sim_beam_cell_copol{3} + sim_beam_cell_copol{2} + sim_beam_cell_copol{1})/3; 0240 % % % sim_beam_all = [sim_angles' sim_beam_sum./max(sim_beam_sum)]; %normalised to max of beam 0241 % % % 0242 % % % %subplot(ifreqs+1,1,ifreqs+1) 0243 % % % plot(sim_beam_all(:,1), 10.*log10(sim_beam_all(:,2)))% Plot in dBm 0244 % % % %plot(sim_beam_all(:,1), 10.*log10(sim_beam_all(:,2)./max(sim_beam_all(:,2))))%Plot in dB 0245 % % % grid on 0246 % % % xlabel('Degrees','Fontsize',24,'Fontname','Helvetica') 0247 % % % ylabel('dBm','Fontsize',24,'Fontname','Helvetica') 0248 % % % title(['Co-polar beam']) 0249 % % % hold all 0250 % 0251 % % Need to turn cell array into matrix to make things easier 0252 % % co-pol first 0253 % mat_copol_power = cell2mat(sim_beam_cell_copol) % This is a power 0254 % 0255 % % Multiply array by weights and then take men across frequency dimension at 0256 % % each angle 0257 % for i=1:length(mat_copol_power) 0258 % mat_copol_power(i,:) = mat_copol_power(i,:).*weights; 0259 % end 0260 % 0261 % mean_copol_power = sum(mat_copol_power,2)/sum(weights); 0262 % %mean_xpol_power = sum(mat_copol_power,2)/sum(weights); 0263 % 0264 % %subplot(ifreqs+1,1,ifreqs+1) 0265 % figure(3) 0266 % plot(sim_angles', 10.*log10( (mean_copol_power)./max(mean_copol_power)))%Plot in dBm (normalised to max of beam) 0267 % %plot(sim_beam_all(:,1), 10.*log10(sim_beam_all(:,2)./max(sim_beam_all(:,2))))%Plot in dB 0268 % grid on 0269 % xlabel('Degrees','Fontsize',24,'Fontname','Helvetica') 0270 % ylabel('dBm','Fontsize',24,'Fontname','Helvetica') 0271 % title(['Mean Co-polar beam']) 0272 % 0273 % % Now x-pol if any data available 0274 % if(Phi45>0) 0275 % mat_xpol_power = cell2mat(sim_beam_cell_xpol) % This is a power 0276 % 0277 % % Multiply array by weights and then take men across frequency dimension at 0278 % % each angle 0279 % for i=1:length(mat_copol_power) 0280 % mat_xpol_power(i,:) = mat_xpol_power(i,:).*weights; 0281 % end 0282 % 0283 % 0284 % mean_xpol_power = sum(mat_xpol_power,2)/sum(weights); 0285 % 0286 % %subplot(ifreqs+1,1,ifreqs+1) 0287 % figure(4) 0288 % plot(sim_angles', 10.*log10( (mean_xpol_power)./max(mean_xpol_power)))%Plot in dBm (normalised to max of beam) 0289 % %plot(sim_beam_all(:,1), 10.*log10(sim_beam_all(:,2)./max(sim_beam_all(:,2))))%Plot in dB 0290 % grid on 0291 % xlabel('Degrees','Fontsize',24,'Fontname','Helvetica') 0292 % ylabel('dBm','Fontsize',24,'Fontname','Helvetica') 0293 % title(['Mean Cross-polar beam']) 0294 % end 0295 0296