read_sim_beam_freq

 Data from Christian is in the form
 Line 1: Field data in cuts
 Line 2: minTheta Theta_step n_points phi_angle 3 1 2

0001 function [mean_power_ZP mean_power_P sim_beam sim_data sim_angles phi vals_sim FWHM_sim fit_angles simfile] = read_beam_sim(n_cuts, fit_lim,object)
0002 
0003 % Data from Christian is in the form
0004 % Line 1: Field data in cuts
0005 % Line 2: minTheta Theta_step n_points phi_angle 3 1 2
0006 
0007 % Number of frequencies and there defintion
0008 
0009 ifreqs = 3;
0010 freqs = {'4.5';'5.0';'5.5'};
0011 %ifreqs = 1;
0012 %freqs = {'5.0'};
0013 weights = ones(1,ifreqs);
0014 alpha = 1;
0015 % make a weighting based on spectrum of TauA
0016 
0017 if (strcmpi('CasA',object))
0018 object
0019 alpha = -0.761; %  CasA Hafez et al. < 15GHz J2000 give -0.761
0020 end
0021 
0022 if (strcmpi('TauA',object))
0023 object
0024 alpha = -0.351; % Tau AHafez et al. 1.8-33GHz give -0.351 (Baars give -0.299 1-35GHz)
0025 end
0026 
0027 freq = [4.5,5.0,5.5];
0028 %freq = [5.0];
0029 weights = weights .*freq(1:ifreqs).^alpha;
0030 
0031 % Now read in each frequency file and get the relevant data
0032 for f=1:ifreqs
0033     freq(f);
0034     simfile =([num2str(freq(f)),'GHz_with_cryo_cone_5mm_up.cut']);
0035     %simfile =([num2str(freq(f)),'GHz_with_cryo_cone_5mm_up_0.4deg_tilt.cut']);
0036     
0037     %simfile =(['/home/cbass/angela/sim_beams/',num2str(freq(f)),'GHz_with_cryo_cone_5mm_up_0.4deg_tilt.cut'])
0038     %simfile =(['/Volumes/Data/CBASS/cbass_analysis_Jan2011/cbass_analysis/angela/cbass_MayJun2011/sim_beams/JPL_galindo_51cm_sec_greg_2.96_ill_5GHz_5mm_up_shifted_secondary.cut'])
0039     
0040 
0041 % First set up the Theta Angles samples
0042 
0043 start_theta = dlmread(simfile,'',[1 0 1 0]); % degrees
0044 theta_step  = dlmread(simfile,'',[1 1 1 1]); % degrees
0045 no_steps    = dlmread(simfile,'',[1 2 1 2]); % degrees
0046 sim_angles  =  linspace(start_theta, -1*start_theta, no_steps);
0047 
0048 
0049 % Now read in each phi cut
0050 
0051 for i=1:n_cuts
0052 phi_position  = [ (1 + (i-1)*(no_steps+2)) 3  (1 + (i-1)*(no_steps+2)) 3  ] ;
0053 data_start = 1+((1 + (i-1)*(no_steps+2))) ;
0054 data_stop = data_start + no_steps - 1;
0055 data_position = [data_start 0 data_stop 3];
0056 phi(i) = dlmread(simfile,'', phi_position); % degrees
0057 sim_data{i} = dlmread(simfile,'',data_position ); % 4 cols: co-polar (Re Im)  x-pol(Re Im)
0058 end
0059 
0060 % For the beam patterns we want to combine data from the 0deg and 90deg phi
0061 % cuts
0062 % Take mean of :
0063 %     zeroPhi_coPol = (Re^2 + Im^2)
0064 %     90Phi_coPol   = (Re^2 + Im^2)
0065 %
0066 % To get what CBASS measures:  0.5*(Ex^2 + Ey^2)
0067 
0068 zeroPhi = find(phi == 0);
0069 Phi90 = find(phi == 90);
0070 
0071 cbass_real_zeroPhi{f} = sim_data{zeroPhi}(:,1);
0072 cbass_imag_zeroPhi{f} = sim_data{zeroPhi}(:,2);
0073 cbass_real_Phi90{f}   = sim_data{Phi90}(:,1);
0074 cbass_imag_Phi90{f}   = sim_data{Phi90}(:,2);
0075 
0076 end
0077 
0078 
0079 % Calculate mean over frequency --> first make a matrix, then get mean across frequencies
0080 % at each angle
0081 %figure
0082 
0083 mat_real_ZP = cell2mat(cbass_real_zeroPhi);
0084 mat_imag_ZP = cell2mat(cbass_imag_zeroPhi);
0085 mat_real_P = cell2mat(cbass_real_Phi90);
0086 mat_imag_P = cell2mat(cbass_imag_Phi90);
0087 
0088 % Multiply data by weights
0089 
0090 for i=1:length(mat_real_ZP);
0091     mat_real_ZP(i,:) = mat_real_ZP(i,:).*weights;
0092     mat_imag_ZP(i,:) = mat_imag_ZP(i,:).*weights;
0093     mat_real_P(i,:) = mat_real_P(i,:).*weights;
0094     mat_imag_P(i,:) = mat_imag_P(i,:).*weights;
0095 end
0096 
0097 mean_real_ZP = sum(mat_real_ZP,2)/sum(weights);
0098 mean_imag_ZP = sum(mat_imag_ZP,2)/sum(weights);
0099 mean_real_P = sum(mat_real_P,2)/sum(weights);
0100 mean_imag_P = sum(mat_imag_P,2)/sum(weights);
0101 
0102 mean_power_ZP = (mean_real_ZP.^2 + mean_imag_ZP.^2);
0103 mean_power_P  = (mean_real_P.^2 + mean_imag_P.^2);
0104 
0105 
0106 sim_beam = [sim_angles' (mean_power_ZP + mean_power_P)/2];
0107 
0108 % Finally, out of interest, let's plot all the phi cuts to see how things
0109 % change
0110 
0111 % figure
0112 % for i = 1:n_cuts
0113 %     plot(sim_angles,(sim_data{i}(:,1).^2 + sim_data{i}(:,2).^2))
0114 %     hold all
0115 % end
0116 %
0117 
0118 % Fit gaussian to the simulated data
0119 fit_lim = 0.8;
0120 fit_angles  = abs(sim_beam(:,1))<fit_lim;
0121 
0122 [best_fit_sim ] = fit_gauss(sim_beam(fit_angles,1),sim_beam(fit_angles,2)./max(sim_beam(fit_angles,2)));
0123 vals_sim = coeffvalues(best_fit_sim);
0124 FWHM_sim = 2*vals_sim(3)*sqrt(log(2));