read_beam_simHart

 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_simHart(n_cuts, fit_lim)
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 = 1;
0010 freqs = {'5.0'};
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 (object == 'CasA')
0018 % object
0019 % alpha = -0.761; %  CasA Hafez et al. < 15GHz J2000 give -0.761
0020 % end
0021 %
0022 % if(object == 'TauA')
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 =(['/usersVol1/act/CBASS/cbass_analysis/cbass_s_matlab_scripts/HartRAO_Beam/horn_pattern_5.0GHz.cut'])
0035     %simfile =([num2str(freq(f)),'GHz_with_cryo_cone_5mm_up.cut'])
0036     %simfile =(['/home/cbass/angela/sim_beams/',num2str(freq(f)),'GHz_with_cryo_cone_5mm_up_0.4deg_tilt.cut'])
0037     %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'])
0038     
0039 
0040 % First set up the Theta Angles samples
0041 
0042 start_theta = dlmread(simfile,'',[1 0 1 0]); % degrees
0043 theta_step  = dlmread(simfile,'',[1 1 1 1]); % degrees
0044 no_steps    = dlmread(simfile,'',[1 2 1 2]); % degrees
0045 sim_angles  =  linspace(start_theta, -1*start_theta, no_steps);
0046 
0047 
0048 % Now read in each phi cut
0049 
0050 for i=1:n_cuts
0051 phi_position  = [ (1 + (i-1)*(no_steps+2)) 3  (1 + (i-1)*(no_steps+2)) 3  ] ;
0052 data_start = 1+((1 + (i-1)*(no_steps+2))) ;
0053 data_stop = data_start + no_steps - 1;
0054 data_position = [data_start 0 data_stop 3];
0055 phi(i) = dlmread(simfile,'', phi_position); % degrees
0056 sim_data{i} = dlmread(simfile,'',data_position ); % 4 cols: co-polar (Re Im)  x-pol(Re Im)
0057 end
0058 
0059 % For the beam patterns we want to combine data from the 0deg and 90deg phi
0060 % cuts
0061 % Take mean of :
0062 %     zeroPhi_coPol = (Re^2 + Im^2)
0063 %     90Phi_coPol   = (Re^2 + Im^2)
0064 %
0065 % To get what CBASS measures:  0.5*(Ex^2 + Ey^2)
0066 
0067 zeroPhi = find(phi == 0);
0068 Phi90 = find(phi == 90);
0069 
0070 cbass_real_zeroPhi{f} = sim_data{zeroPhi}(:,1);
0071 cbass_imag_zeroPhi{f} = sim_data{zeroPhi}(:,2);
0072 cbass_real_Phi90{f}   = sim_data{Phi90}(:,1);
0073 cbass_imag_Phi90{f}   = sim_data{Phi90}(:,2);
0074 
0075 end
0076 
0077 
0078 % Calculate mean over frequency --> first make a matrix, then get mean across frequencies
0079 % at each angle
0080 %figure
0081 
0082 mat_real_ZP = cell2mat(cbass_real_zeroPhi);
0083 mat_imag_ZP = cell2mat(cbass_imag_zeroPhi);
0084 mat_real_P = cell2mat(cbass_real_Phi90);
0085 mat_imag_P = cell2mat(cbass_imag_Phi90);
0086 
0087 % Multiply data by weights
0088 
0089 for i=1:length(mat_real_ZP)
0090     mat_real_ZP(i,:) = mat_real_ZP(i,:).*weights;
0091     mat_imag_ZP(i,:) = mat_imag_ZP(i,:).*weights;
0092     mat_real_P(i,:) = mat_real_P(i,:).*weights;
0093     mat_imag_P(i,:) = mat_imag_P(i,:).*weights;
0094 end
0095 
0096 mean_real_ZP = sum(mat_real_ZP,2)/sum(weights);
0097 mean_imag_ZP = sum(mat_imag_ZP,2)/sum(weights);
0098 mean_real_P = sum(mat_real_P,2)/sum(weights);
0099 mean_imag_P = sum(mat_imag_P,2)/sum(weights);
0100 
0101 mean_power_ZP = (mean_real_ZP.^2 + mean_imag_ZP.^2);
0102 mean_power_P  = (mean_real_P.^2 + mean_imag_P.^2);
0103 
0104 
0105 sim_beam = [sim_angles' (mean_power_ZP + mean_power_P)/2];
0106 
0107 % Finally, out of interest, let's plot all the phi cuts to see how things
0108 % change
0109 
0110 % figure
0111 % for i = 1:n_cuts
0112 %     plot(sim_angles,(sim_data{i}(:,1).^2 + sim_data{i}(:,2).^2))
0113 %     hold all
0114 % end
0115 %
0116 
0117 % Fit gaussian to the simulated data
0118 fit_lim = 20
0119 fit_angles  = abs(sim_beam(:,1))<fit_lim;
0120 
0121 [best_fit_sim ] = fit_gauss(sim_beam(fit_angles,1),sim_beam(fit_angles,2)./max(sim_beam(fit_angles,2)));
0122 vals_sim = coeffvalues(best_fit_sim);
0123 FWHM_sim = 2*vals_sim(3)*sqrt(log(2));