Home > cbassSouthFunctions > ACTFunctions > beam_fittingSA.m

beam_fittingSA

PURPOSE ^

CHECK BEFORE USING - WAS FIDDLING WITH NORMALISATION ACT 12/1/2011 %

SYNOPSIS ^

function [d FWHM FWHM_bin FWHM_sim angle_bins all_data corr_data source vals d_on_az filedate Offset_bin sim_beam Offset] = beam_fitting(d,phase, direction, simfile, n_cuts,sim_beam, sim_data, sim_angles ,phi, vals_sim ,FWHM_sim, sim_fit_angles,fake,fitoff,flags,bin_width,maxfit_angle)

DESCRIPTION ^

 CHECK BEFORE USING - WAS FIDDLING WITH NORMALISATION ACT 12/1/2011 %

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [d FWHM FWHM_bin FWHM_sim angle_bins all_data corr_data source vals d_on_az filedate Offset_bin sim_beam Offset] = beam_fitting(d,phase, direction, simfile, n_cuts,sim_beam, sim_data, sim_angles ,phi, vals_sim ,FWHM_sim, sim_fit_angles,fake,fitoff,flags,bin_width,maxfit_angle)
0002 
0003 
0004 % CHECK BEFORE USING - WAS FIDDLING WITH NORMALISATION ACT 12/1/2011 %
0005 
0006 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0007 % function[FWHM] = pipeline_el(d,direction)
0008 %
0009 % Function to analyse beam scan data
0010 % Need to specify the data stream, d
0011 % Also need to specify if it is azimuth or elevation scans
0012 % phase      = 1 (old switching calibration) 2(new switching calibration)
0013 % direction = 1 (elevation)  2 (azimuth)
0014 %
0015 % Outputs are plots of the beam for all channels and gaussian fit to FWHM
0016 % of
0017 % these beams (binned up and raw) plus simulated FWHM from Christians beam
0018 % simulations
0019 %
0020 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0021 
0022 % To look at beam azimuth scans
0023 %
0024 % 16-Sep-2010:17:08:23      16-Sep-2010:17:11:58      source_scans.sch( 1h:12m       , taua )
0025 % 16-Sep-2010:17:11:58      16-Sep-2010:17:12:02      source_scans.sch( 2h   ,
0026 % taua )
0027 % 16-Sep-2010:17:12:18      16-Sep-2010:18:38:45      source_scans.sch( 1h:12m
0028 % , taua )
0029 % 16-Sep-2010:18:38:45      16-Sep-2010:21:04:53      source_scans.sch( 2h   , 3c274 )
0030 % 16-Sep-2010:21:04:53      16-Sep-2010:21:04:59      source_raster.sch( 2h   , 3c274 )
0031 % 16-Sep-2010:21:04:59      16-Sep-2010:23:49:13      source_scans.sch( 2h   , 3c274 )
0032 % 16-Sep-2010:23:49:13      17-Sep-2010:02:04:21      source_raster.sch( 2h   , 3c274 )
0033 % 17-Sep-2010:02:04:21      17-Sep-2010:04:18:23      source_raster.sch( 2h   ,
0034 % casa )
0035 % 17-Sep-2010:04:18:23      17-Sep-2010:07:01:27      source_scans.sch( 2h   , casa )
0036 % 17-Sep-2010:07:01:27      17-Sep-2010:09:15:03      source_raster.sch( 2h   , casa )
0037 % 17-Sep-2010:09:15:03      17-Sep-2010:11:55:04      source_scans.sch( 2h   , casa )
0038 % 17-Sep-2010:11:55:04      17-Sep-2010:11:59:26
0039 % 17-Sep-2010:11:59:26      17-Sep-2010:13:17:31      source_raster.sch( 1h   , casa )
0040 % 17-Sep-2010:13:17:31      17-Sep-2010:14:31:47      source_scans.sch( 1h   , casa )
0041 % 17-Sep-2010:14:31:47      17-Sep-2010:16:45:46      source_raster.sch( 2h   , taua )
0042 % 17-Sep-2010:16:45:46      17-Sep-2010:19:00:32      source_scans.sch( 2h   , taua )
0043 % 17-Sep-2010:19:00:32      17-Sep-2010:21:17:50      source_raster.sch( 2h   , taua )
0044 %d = pipe_read('21-Oct-2010:04:02:53','      21-Oct-2010:06:02:58'); % Cyg A
0045 %d = pipe_read('21-Oct-2010:01:03:48',      '21-Oct-2010:03:01:12'); %CasA
0046 %d = pipe_read('21-Oct-2010:06:02:58','      21-Oct-2010:07:51:06'); % TauA
0047 % TauA
0048 
0049 
0050 %clear;
0051 %close all;
0052 %d = pipe_read('21-Oct-2010:04:02:53','      21-Oct-2010:04:42:58');
0053 %d = pipe_read('21-Oct-2010:04:02:53','      21-Oct-2010:06:02:58'); % Cyg A
0054 %d = pipe_read('17-Sep-2010:16:45:46','17-Sep-2010:19:00:32');   %TauA
0055 %d = pipe_read('17-Sep-2010:09:15:03','17-Sep-2010:11:55:04');  %CassA
0056 %d = pipe_read('29-Sep-2010:13:21:32','      29-Sep-2010:16:10:13 ') %Moon
0057 %d = pipe_read('21-Oct-2010:04:02:53','      21-Oct-2010:06:02:58'); % Cyg A
0058 %d = pipe_read('21-Oct-2010:01:03:48',      '21-Oct-2010:03:01:12'); %CasA
0059 %%%%d = pipe_read('21-Oct-2010:06:02:58','      21-Oct-2010:07:51:06'); % TauA
0060 %direction = 2;
0061 %d = pipe_read('29-Sep-2010:13:21:32','29-Sep-2010:16:10:13');              % Moon
0062 %d = pipe_read('15-Jun-2010:00:06:55','15-Jun-2010:02:25:16');
0063 
0064 
0065 %d = pipe_read('04-Dec-2010:03:25:14','04-Dec-2010:04:38:18'); % CygA
0066 %d = pipe_read('04-Dec-2010:06:52:55','04-Dec-2010:09:06:45'); % CasA
0067 %d = pipe_read('04-Dec-2010:10:56:31','04-Dec-2010:11:52:10'); % TauA
0068 %d = pipe_read('04-Dec-2010:14:36:22','04-Dec-2010:17:30:58'); % TauA
0069 %d = pipe_read('04-Dec-2010:10:56:31','04-Dec-2010:17:30:58'); % All TauA
0070 %d = pipe_read('05-Dec-2010:19:22:18','05-Dec-2010:09:05:35'); % VirA
0071 %d = pipe_read('16-Dec-2010:06:54:39','      16-Dec-2010:08:08:05') % CasA
0072 
0073 %%d = pipe_read('04-Dec-2010:06:52:55','04-Dec-2010:09:06:45'); % CasA
0074 %d = pipe_read('16-Dec-2010:06:54:39','      16-Dec-2010:08:08:05') % CasA
0075 
0076 
0077 %% Let's apply a non-linearty correction to the .switchData
0078 if (fake==0)
0079 %Measured_Power = calibrate_linearity(d);
0080 %d.antenna0.receiver.switchData = Measured_Power*1000;
0081 %clear Measured_Power
0082 
0083 
0084 %%  Original calibration - don't use phase switch averaging
0085 
0086     if (phase==1)
0087         %rval = [1.9,1.8,1.8,1.7];
0088         rval = [2.4, 2.3, 2.0, 2.1];
0089 
0090         [thisAlpha r d meanr] = pol_correction(d,rval);
0091     else
0092 %%
0093  
0094 %%%%%%%d = reduceData(d,'redScript_basic.red',0,'beam_plots')
0095 %%%%%%%%d = getIV(d);
0096 
0097 % d = reduceData(d,'redScript_basic.red',0,'beam_plots');
0098        %d = getIV(d);
0099       %d = dr;
0100        
0101       %d= assembleAlphaStreams_PolOnly(d);
0102       %[t,A,G,T,r,horz,equa,offS,onE,offE,onS]=calculateAlpha_PolOnly(d);
0103       %d = applyAlpha_PolOnly(d,t,A,G);%
0104       
0105 %      d = assembleAlphaStreams_PolOnly(d);
0106 %      d = applyAlpha_database_PolOnly(d);
0107 %      d = calculateStokes_PolOnly(d);
0108 %r(1)=1;
0109 %r(2)=1;
0110 %%d = reduceData(d,'redScript_basic.red',0,'beam_plots');
0111 %d = notchFilter(d);
0112 %%d = assembleAlphaStreams_DD(d);
0113 %%[t,A,G,T,r,horz,equa,offS,onE,offE,onS]=calculateAlpha_DD(d);
0114 %%d = applyAlpha_DD(d,t,A,G);
0115 %%IL1 = d.antenna0.receiver.data(:,1)-(r(1)*d.antenna0.receiver.data(:,2));
0116 %%IL2 = d.antenna0.receiver.data(:,9)-(r(2)*d.antenna0.receiver.data(:,10));
0117 %%I = (IL1 + IL2)/2;
0118 %%V = (IL1 - IL2)/2;
0119 %%Q = d.antenna0.receiver.data(:,6);
0120 %%U = d.antenna0.receiver.data(:,7);
0121 %%P = sqrt(Q.^2 + U.^2);
0122 
0123 %%d.antenna0.receiver.data(:,1) = I;
0124 %%d.antenna0.receiver.data(:,8) = V;
0125 %%d.antenna0.receiver.data(:,5) = P;
0126     end
0127 
0128 
0129     % Only use data on the source  - not slewing
0130     if (direction == 1)
0131         Type = 'Elevation';
0132         dcut = framecut(d,d.index.elscan.fast);
0133         max_angle = 3;
0134         min_angle = 0.8;
0135         %maxfit_angle=1.5;
0136     else
0137         Type = 'Azimuth';
0138         dcut = framecut(d,d.index.source.fast);
0139         min_angle = 0.9; 
0140         max_angle = 3;
0141         %maxfit_angle = 1.5;%% a baseline offset will be fitted out for data inthe range 1deg --5deg offset (ie avoiding the main beam data which can mess things up)
0142     end
0143 else
0144     min_angle = 0.7;
0145     max_angle = 3;
0146     maxfit_angle = 0.9;
0147     dcut = d;
0148   if (direction == 2)
0149       Type = 'Azimuth';
0150       if(direction ==1)
0151           Type = 'Elevation';
0152       end
0153   end
0154 end
0155 %%
0156 %Any extra flags to put ont he data - these are for the 4th June CasA data
0157 if(flags==1)
0158 dcut.flags.fast(3.869e4:4.964e4)=1;
0159 dcut.flags.fast(6.685e4:7.484e4)=1;
0160 dcut= framecut(dcut,~dcut.flags.fast(:,1));
0161 plot(dcut.antenna0.receiver.data(:,1))
0162 end
0163 %%
0164 % Also cut out all the flagged data - based on channel 1
0165 
0166 %dcut = framecut(dcut,~dcut.flags.fast(:,1));
0167 
0168 % Do a cut on elevation % might not be strictly necessary - hangover from
0169 % old fitting routine
0170 
0171 %el_limit = 0;
0172 %dcut = framecut(dcut,dcut.antenna0.servo.el>el_limit);
0173 
0174 %clear d;
0175 
0176 time = dcut.antenna0.roach1.utc;
0177 az_scan = dcut.antenna0.servo.apparent(:,1);
0178 el_scan = dcut.antenna0.servo.apparent(:,2);
0179 
0180 %Point sources
0181 
0182 [az_src, el_src] = calcAzElCs(dcut);
0183 angle = spaceangle(az_src, el_src, az_scan,el_scan,'deg');
0184 
0185 %Moon
0186 %[moon_angle az_src el_src]=calcSourceDist2(dcut,'moon');
0187 %angle = moon_angle;
0188 %deltaAz = cosd(el_scan).*(az_scan - az_src);
0189 
0190 deltaEl = el_scan - el_src;
0191 
0192 % Make angle -ve for scans where az_scan < az_src
0193 
0194 if (direction == 1)
0195     disp('Will analyse elevation scans')
0196     angle(el_scan<el_src) = -1.*angle(el_scan<el_src);
0197 else
0198     disp('Will analyse azimuth scans')
0199     angle(az_scan<az_src) = -1.*angle(az_scan<az_src);
0200 end
0201 
0202 %angle = az_scan-az_src;
0203 
0204 %%
0205 % try this for the offsets - taken from calcOffSan.m
0206 
0207   % get the offsets from the control system.
0208   azApp = interp1(dcut.antenna0.tracker.utc, ...
0209   dcut.antenna0.tracker.horiz_topo(:,1), dcut.antenna0.receiver.utc);
0210   azOffSave = wrap180(azApp) - wrap180(dcut.antenna0.servo.apparent(:,1));
0211   azOffSave = -azOffSave;
0212 
0213   elApp = interp1(dcut.antenna0.tracker.utc, ...
0214   dcut.antenna0.tracker.horiz_topo(:,2), dcut.antenna0.receiver.utc);
0215   elOffSave = elApp - dcut.antenna0.servo.apparent(:,2);
0216   elOffSave = -elOffSave;
0217 
0218   if(direction==2) % azimuth
0219   angle = azOffSave .* cos(dcut.antenna0.servo.el*pi/180);
0220   else
0221       angle = elOffSave;
0222   end
0223      
0224   
0225   % Out of interest, try full angle ie sqrt(az^2 + el^2)
0226   %angle = sqrt((azOffSave .* cos(dcut.antenna0.servo.el*pi/180)).^2 + (elOffSave.^2));
0227    %angle = interp1(dcut.antenna0.tracker.utc,dcut.antenna0.tracker.scan_off(:,1),time);
0228 %%
0229 % Set up the bins for binning the data
0230 %bin_width = 0.1;
0231 az_bins = [-max_angle:bin_width:max_angle];
0232 n_bins = length(az_bins(1,:));
0233 deltaEl = angle;
0234 
0235 % Polarization leakage coeeficeints
0236 % Also have an amp_fix to get all channels going positive
0237 
0238 
0239 if (phase ==1)
0240 ampfix = [1 1 -1 1 -1 1];
0241 C = [0 0.0173 -0.00845 0.0185 -0.01127 0];
0242 else
0243 % If after 4-June:
0244 %ampfix = [-1 -1 1 -1 1 -1 1 1];
0245 ampfix = [1 1 1 1 1 1 1 1 1  ];
0246 %else
0247 %ampfix=[1 1 1 1 1 -1 1 1];
0248 C = [0 0 0 0 0 0 0 0 ];
0249 end
0250 
0251 
0252 
0253 
0254 
0255 %% Loop over channels and bin up the beam data
0256 if (phase == 1)
0257     no_chan = [1:6];
0258     x_plots = 2;
0259 else
0260     %no_chan = [1,6,7,8]
0261     no_chan = [1,2,3,4];
0262     x_plots = 2;
0263     b =0;
0264 end;
0265 
0266 for chan = no_chan
0267    chan
0268    angle;
0269    b = b+1
0270    figure(b)
0271     %if (phase==1)
0272     %    subplot(x_plots,3,chan)
0273     %else
0274     %b = b+1;
0275     %subplot(x_plots,4,b)
0276     %subplot(x_plots,2,b)
0277     %end
0278   %%
0279 
0280 %Plot of the raw scan data versus time
0281 %figure(direction+3)
0282 subplot(x_plots,4,1:2)
0283  plot(dcut.antenna0.receiver.utc,dcut.antenna0.receiver.data(:,chan),'b')
0284  titleText=[Type,' Scan data on ',dcut.antenna0.tracker.source{1},10,mjd2string(dcut.antenna0.receiver.utc(1)),' ',mjd2string(dcut.antenna0.receiver.utc(end))];
0285  title(titleText)
0286  xlabel('UTC');
0287  ylabel('Nominal units')  
0288  hold on   
0289 
0290 %%  New bit to remove a linear trend for each scan individually trend in elevation
0291 if (fitoff ==1)
0292 % [dcut angle_new] = temp_fitoff(dcut,angle,az_scan,el_scan,time,maxfit_angle,min_angle,chan,direction);
0293 % else
0294 %     angle_new = angle;
0295 % end
0296 
0297 % %% Alternative detrending
0298 %
0299 inds = zeros(size(dcut.antenna0.receiver.utc));
0300 inds = logical(inds);
0301 if (direction == 1) % elevation scans
0302     ss = abs(gradient(dcut.antenna0.tracker.scan_off(:,2)));
0303     scan_off = interp1(dcut.antenna0.tracker.utc,ss,dcut.antenna0.receiver.utc);
0304     inds(scan_off>0.49)=1;
0305 else
0306     ss = abs(gradient(dcut.antenna0.tracker.scan_off(:,1)));
0307     scan_off = interp1(dcut.antenna0.tracker.utc,ss,dcut.antenna0.receiver.utc);
0308     inds(scan_off>1)=1;
0309 end
0310 
0311 
0312 %plot(scan_off);
0313 
0314 %plot(inds)
0315 
0316 % We can now identify each scan as I(cumsum==i), where i=1:max(indsum)
0317 %plot(dcut.antenna0.receiver.data(indsum==10,1))
0318 
0319 [s e] = findStartStop(inds);
0320 %dgood = dcut.antenna0.receiver.data;
0321 % Grab each scan in turn and remove a slope
0322 
0323 for i=1:length(s)
0324     %for chan=1:1
0325 
0326         figure(21)
0327         temp_y = dcut.antenna0.receiver.data(s(i):e(i),chan);
0328         temp_x = angle(s(i):e(i));
0329         subplot(9,9,i)
0330         plot(temp_x,temp_y,'r')
0331 
0332         p = polyfit(temp_x,temp_y,1);
0333         f = polyval(p,temp_x);
0334 
0335         hold on
0336         plot(temp_x,f,'b')
0337         plot(temp_x,(temp_y-f),'g')
0338 
0339         dcut.antenna0.receiver.data(s(i):e(i)) = dcut.antenna0.receiver.data(s(i):e(i),chan) - f;
0340     %end
0341 end
0342 angle_new = angle;
0343 %Remove unspecified angles at the beginning of the observation
0344 angle_new(1:s(1))=100;
0345 figure(22)
0346 plot(angle_new,dcut.antenna0.receiver.data(:,chan),'b')
0347 
0348 else
0349      angle_new = angle;
0350  end
0351 %
0352 %%
0353 
0354 figure(b)
0355 subplot(x_plots,4,1:2)
0356 plot(dcut.antenna0.receiver.utc,dcut.antenna0.receiver.data(:,chan),'g')
0357 hold off
0358 %plot(angle,dcut.antenna0.receiver.data(:,1))
0359 %figure
0360 d_on_az = (dcut.antenna0.receiver.data(:,chan) *ampfix(chan));
0361 
0362 % Let us throw away data with deltaEl > max(az_bins)
0363 deltaEl = angle_new;
0364 deltaElshort = deltaEl(abs(deltaEl)<max(az_bins(1,:)));
0365 
0366 
0367 %%
0368     d_on_az = d_on_az(abs(deltaEl)<max(az_bins(1,:)));
0369     %figure;plot(d_on_az)
0370     %d_on_az(d_on_az>2)=nan;
0371     %d_on_az= d_on_az(d_on_az<2);
0372     %deltaElshort = deltaElshort(d_on_az<2);
0373     %figure;
0374     %plot(deltaElshort,d_on_az)
0375     az_bins(2,:) = 0;
0376     az_bins(3,:) = 0;
0377     bins_az = ceil(deltaElshort/bin_width) + ceil(n_bins/2);
0378 
0379     for i = 1:length(az_bins)
0380       
0381         az_bins(2,i) = mean(d_on_az(bins_az==i));
0382         number = length(d_on_az(bins_az==i));
0383         az_bins(3,i) = (std(d_on_az(bins_az==i)))/sqrt(number);
0384                 
0385     end
0386 
0387     %clear d_on_az
0388      
0389     if (min(az_bins(2,:)) < 0) 
0390          az_bins(2,:) = az_bins(2,:)+(abs(min(az_bins(2,:))));
0391      else
0392          az_bins(2,:) = az_bins(2,:)-(abs(min(az_bins(2,:))));
0393      end
0394      
0395     all_data(:,chan) = az_bins(2,:);
0396     %errorbar(az_bins(1,:),az_bins(2,:),az_bins(3,:));
0397     %hold on
0398     
0399     %xlabel('Azimuth offset, degrees')
0400     %ylabel('Nominal units')
0401     
0402 
0403 
0404 %% Lets fit with a gaussian
0405 % Check that we don't have any Nan's in az_bins 2 (ie where no data)
0406 % Let's only fit to the inner 1deg - this is only done for the unbinned
0407 % data xx,yy
0408 
0409 if (chan==1)
0410 fit_lim = 0.8;
0411 fit_angles  = abs(deltaElshort)<fit_lim;
0412 fit_angle_bins = abs(az_bins(1,:))< fit_lim;
0413 
0414 temp = isnan(az_bins(2,:));
0415 az_bins(2,temp) = 0;
0416 temp1 = isnan(az_bins(1,:));
0417 az_bins(1,temp1) = 0;
0418 temp3 = isnan(az_bins(1,:));
0419 az_bins(3,temp3) = 0;
0420 
0421 x = az_bins(1,fit_angle_bins)';
0422 y = az_bins(2,fit_angle_bins)';
0423 we = az_bins(3,fit_angle_bins)';
0424 
0425 xx = deltaElshort(fit_angles);
0426 yy = (d_on_az(fit_angles));
0427 disp('here1')
0428 [best_fit_bin ] = fit_gauss(x,y);
0429 disp('here2')
0430 [best_fit ] = fit_gauss(xx,yy);
0431 disp('here3')
0432 vals_bin = coeffvalues(best_fit_bin);
0433 vals = coeffvalues(best_fit);
0434 
0435 end
0436 a1(chan) = vals_bin(1);
0437 b1(chan) = vals_bin(2);
0438 c1(chan) = vals_bin(3);
0439 FWHM_bin(chan) = 2*vals_bin(3)*sqrt(log(2));
0440 FWHM(chan) = 2*vals(3)*sqrt(log(2));
0441 Offset(chan) = vals(2);
0442 Offset_bin(chan)=vals_bin(2)
0443 
0444 
0445 end
0446 
0447 %%
0448 % Now plot all the polarization corrected data beams on figure direction*2
0449 if (phase==1)
0450     meanI = (all_data(:,1) + all_data(:,10))/2;
0451 else
0452     meanI = all_data(:,chan);
0453     b=0;
0454 end
0455 %%%%%%%%%%%%%%%%%%%%
0456 %Put this back in to try to correct for polarization leakage
0457 %     for chan = no_chan
0458 %         corr_data(:,chan) = (all_data(:,chan) - meanI*C(chan));
0459 %     hold on
0460 %     if (phase ==1)
0461 %         subplot(x_plots,3,chan)
0462 %     else
0463 %         b=b+1;
0464 %         subplot(x_plots,5,b)
0465 %     %subplot(x_plots,2,b)
0466 %     end
0467 %%%%%%%%%%%%%%%%%%
0468 % figure(direction*2)
0469 %     for chan = no_chan
0470 %         corr_data(:,chan) = all_data(:,chan);
0471 %     %if (phase ==1)
0472 %     %     subplot(x_plots,3,chan)
0473 %     % else
0474 %     %     b=b+1;
0475 %     %     subplot(x_plots,4,b)
0476 %     % end
0477 %     subplot(1,6,chan)
0478 %     errorbar(az_bins(1,:),corr_data(:,chan),az_bins(3,:))
0479 %     b=1;
0480 %     %plot(az_bins(1,:),corr_data(:,chan),'r');
0481 %     if (b==1)
0482 %        titleText=[Type,' Scan through ',dcut.antenna0.tracker.source{1},10,mjd2string(dcut.antenna0.receiver.utc(1)),' ',mjd2string(dcut.antenna0.receiver.utc(end)),10, 'FWHM_bin = ', num2str(FWHM_bin(chan)),10,'Offset= ', num2str(Offset_bin(chan))];
0483 %     else
0484 %        titleText=[' ',10,' ',10, 'FWHM_bin = ', num2str(FWHM_bin(chan)),10,'Offset= ', num2str(Offset_bin(chan))];
0485 %     end
0486 %
0487 %     title(titleText)
0488 %     grid on
0489 %     hold on
0490 % % Plot the best fit gaussian too
0491 %
0492 %     gauss = a1(chan)*exp(-((x-b1(chan))/c1(chan)).^2);
0493 %     plot(x,gauss,'g')
0494 %     xlim([-3 3])
0495 %     end
0496 % Also put on the intensity plot
0497 
0498 for chan = no_chan
0499     b=b+1
0500      figure(b)
0501         corr_data(:,chan) = all_data(:,chan);
0502     %if (phase ==1)
0503     %     subplot(x_plots,3,chan)
0504     % else
0505     %     b=b+1;
0506     %     subplot(x_plots,4,b)
0507     % end
0508     subplot(x_plots,4,3:4)
0509     errorbar(az_bins(1,:),corr_data(:,chan),az_bins(3,:))
0510     c = 1
0511     %plot(az_bins(1,:),corr_data(:,chan),'r');
0512     if (c==1)
0513        titleText=[Type,' Scan through ',dcut.antenna0.tracker.source{1},10,mjd2string(dcut.antenna0.receiver.utc(1)),' ',mjd2string(dcut.antenna0.receiver.utc(end)),10, 'FWHM_bin = ', num2str(FWHM_bin(chan)),10,'Offset= ', num2str(Offset_bin(chan))];
0514     else
0515        titleText=[' ',10,' ',10, 'FWHM_bin = ', num2str(FWHM_bin(chan)),10,'Offset= ', num2str(Offset_bin(chan))]; 
0516     end
0517     
0518     title(titleText)
0519     grid on
0520     hold on
0521 % Plot the best fit gaussian too
0522     if(chan==1)
0523     gauss = a1(chan)*exp(-((x-b1(chan))/c1(chan)).^2);
0524     plot(x,gauss,'g')
0525     xlim([-max_angle max_angle])
0526     %end
0527 
0528 
0529 FWHM_bin(chan);
0530 FWHM(chan);
0531     end
0532 angle_bins = az_bins(1,:);
0533 source = dcut.antenna0.tracker.source{1};
0534 filedate = mjd2string(dcut.antenna0.receiver.utc(1));
0535 
0536 %%
0537 % Now plot the measured and simulated beam on the same plot
0538 figure(b)
0539 subplot(x_plots,4,5:6)
0540 chan
0541 hold on
0542 x = sim_beam(sim_fit_angles,1);
0543 gauss_sim = vals_sim(1)*exp(-((x-vals_sim(2))/vals_sim(3)).^2);
0544 plot(x,(gauss_sim),'g')
0545 xlim([-max_angle max_angle])
0546 
0547 
0548 FWHM_sim;
0549 
0550 
0551 plot(sim_beam(:,1),(sim_beam(:,2)./max(sim_beam(:,2))),'b')
0552 plot(x,gauss_sim,'gp')
0553 
0554 title(['Plot of simulated and measured C-BASS beam (Linear)',10,Type,' Scan through ',dcut.antenna0.tracker.source{1},10,mjd2string(dcut.antenna0.receiver.utc(1)),' ',mjd2string(dcut.antenna0.receiver.utc(end))] )
0555 
0556 % Plot the shifted CBASS measured beam
0557 corr_data11 = corr_data(:,1);
0558 corr_data1 = corr_data(:,chan);
0559 min_corrdata = 0;
0560 Offset(chan)
0561 Offset_bin(chan)
0562 if (Offset(chan) < 0)
0563     plot(angle_bins+abs(Offset_bin(chan)),((corr_data1-min_corrdata)/max(corr_data11-min_corrdata)),'rp')
0564 else
0565   plot(angle_bins-abs(Offset_bin(chan)),((corr_data1-min_corrdata)/max(corr_data11-min_corrdata)),'rp')  
0566 end
0567 
0568 xlabel('Offset angle, degrees')
0569 ylabel('Normalised units, I channel')
0570 text(0.1,1.3,['\color{red}Measured Beam'])
0571 text(0.1,1.2,['\color{blue}Simulated Beam'])
0572 text(0.1,1.1,['\color{green}Gaussian fit to simulated beam'])
0573 grid on
0574 xlim([-max_angle max_angle])
0575 
0576 %%zoomed in plot of sidelobes
0577 
0578 subplot(x_plots,4,7:8)
0579 
0580 plot(sim_beam(:,1),10.*log10(sim_beam(:,2)./max(sim_beam(:,2))),'b')
0581 hold on
0582 plot(x,10.*log10(gauss_sim),'gp')
0583 grid on
0584 title(['Plot of simulated and measured C-BASS beam (dB)',10,Type,' Scan through ',dcut.antenna0.tracker.source{1},10,mjd2string(dcut.antenna0.receiver.utc(1)),' ',mjd2string(dcut.antenna0.receiver.utc(end))] )
0585 if (Offset(chan) < 0)
0586     plot(angle_bins+abs(Offset_bin(chan)),10.*log10((corr_data1-min_corrdata)/max((corr_data11-min_corrdata))),'rp')
0587 else
0588   plot(angle_bins-abs(Offset_bin(chan)),10.*log10((corr_data1-min_corrdata)/max((corr_data11-min_corrdata))),'rp')  
0589 end
0590 xlabel('Offset angle, degrees')
0591 ylabel('Normalised units, I channel')
0592 text(1,0,['\color{red}Measured Beam'])
0593 text(1,2,['\color{blue}Simulated Beam'])
0594 text(1,4,['\color{green}Gaussian fit to simulated beam'])
0595 xlim([-max_angle max_angle])
0596 end
0597 FWHM_sim 
0598 FWHM_bin 
0599 FWHM_all=FWHM_bin 
0600 FWHM_bin
0601 FWHM_all
0602 
0603 
0604 
0605 
0606 
0607 
0608 
0609
Generated on Sun 14-Jun-2015 17:12:45 by m2html © 2005