30/3/2015 - ACT
To plot cuts of the CBASS beam as measured from the Moon
Example schedules are:
source_scan_az20deg.sch( 4h , moon )
Saves a .mat file containing most of the outputs and the raw and
detrended data so you can quickly re-load it in.
Also saves the final .png images of the beam and simulated beam cut.
% Some scans from March 20150001 0002 function [d,dcut,TauAmap, azOffSave, elOffSave,xidx,yidx,sim_beam_cell_copol,sim_angles] = north_moon_beam(start_date,stop_date,comment) 0003 % 30/3/2015 - ACT 0004 % To plot cuts of the CBASS beam as measured from the Moon 0005 % Example schedules are: 0006 % source_scan_az20deg.sch( 4h , moon ) 0007 % 0008 % Saves a .mat file containing most of the outputs and the raw and 0009 % detrended data so you can quickly re-load it in. 0010 % Also saves the final .png images of the beam and simulated beam cut. 0011 %% Some scans from March 2015 0012 0013 % If youdon't specify start and stop dates in the call to this function, it 0014 % will use whatever dates are in this list: 0015 0016 if(nargin<2) 0017 dates={...%'13-Mar-2015:13:51:59', '13-Mar-2015:15:54:32';... % /home/cbass/gcpCbass/control/sch/source_scan_az20deg.sch( 3h , moon ) 0018 %'13-Mar-2015:16:51:57', '13-Mar-2015:17:00:24';... /home/cbass/gcpCbass/control/sch/source_scan_az20deg.sch( 7h , virgoa ) 0019 %'14-Mar-2015:13:26:18', '14-Mar-2015:16:41:02';... % /home/cbass/gcpCbass/control/sch/source_scan_az20deg.sch( 3h , moon ) 0020 '16-Mar-2015:14:45:44', '16-Mar-2015:18:59:46';... % /home/cbass/gcpCbass/control/sch/source_scan_az20deg.sch( 4h , moon ) 0021 %'16-Mar-2015:19:15:13', '16-Mar-2015:23:52:33';... % /home/cbass/gcpCbass/control/sch/source_scan_az20deg.sch( 4h , casa ) 0022 %'17-Mar-2015:00:14:19', '17-Mar-2015:03:58:30';... % /home/cbass/gcpCbass/control/sch/source_scan_az20deg.sch( 3h , taua ) 0023 %'17-Mar-2015:04:13:46', '17-Mar-2015:06:46:26';... % /home/cbass/gcpCbass/control/sch/source_scan_az20deg.sch( 2h , taua ) 0024 %'17-Mar-2015:15:34:48', '17-Mar-2015:20:19:22';... % /home/cbass/gcpCbass/control/sch/source_scan_az20deg.sch( 4h , moon ) 0025 %'19-Mar-2015:18:20:46', '19-Mar-2015:19:15:18';... % /home/cbass/gcpCbass/control/sch/source_scan_az20deg.sch( 4h , moon ) 0026 %'19-Mar-2015:19:15:56', '19-Mar-2015:22:55:55'}; % /home/cbass/gcpCbass/control/sch/source_scan_az20deg.sch( 3h , moon ) 0027 } 0028 else 0029 dates={start_date,stop_date}; 0030 end 0031 0032 for rasters=1:size(dates,1) 0033 start_date = dates{rasters,1}; 0034 stop_date = dates{rasters,2}; 0035 0036 % read in the data and reduce using a raster script 0037 d = read_arc(start_date,stop_date) 0038 d = reduceData(d,'redScript_rasters_March15.red') 0039 0040 % re-label the indices so that we can use old code: 0041 0042 d.index.beammap.fast=d.index.source.fast; 0043 0044 % plot the data to check we have something 0045 figure 0046 plot(d.antenna0.receiver.data(:,1)) 0047 source_name=char(d.antenna0.tracker.source(10)) 0048 0049 % Some labelling for later plots 0050 nchan=4 0051 plot_labels = {'I1','I2','Q','U'}; 0052 chan_order=[1 8 6 7]; 0053 %% 0054 % Detrend the data ready for analysis - This is the bit we should probably 0055 % change so that it reomves the sky background from our actual maps... 0056 0057 % detrend data using data beyond a certain offset angle - default is 2deg 0058 offset_angle =10 0059 [s e] = findStartStop(d.index.beammap.fast); % Find start and end points of each raster 0060 dgood = d; % make a master copy which will later contain the detrended data (see in for loop below) 0061 n=0 0062 0063 0064 for m= 1:length(s) 0065 n = n+1; 0066 ind = zeros(size(d.index.beammap.fast)); 0067 ind(s(m):e(m)) = 1; 0068 ind = logical(ind); 0069 0070 %if(strmatch('S',cbass,'exact')) 0071 % dcut = framecutSa(d, ind); 0072 0073 %else 0074 dcut = framecut(d, ind); 0075 %end 0076 azApp = interp1(dcut.antenna0.tracker.utc, ... 0077 dcut.antenna0.tracker.horiz_topo(:,1), dcut.antenna0.receiver.utc); 0078 azOffSave = (azApp) - dcut.antenna0.servo.apparent(:,1); %%% Wrap180 put in for SOuth data 23/Aug/2013 - CHECK 0079 azOffSave = -azOffSave; 0080 0081 azOffSave = azOffSave.*cos(dcut.antenna0.servo.el*pi/180);% Scale azimuth by cos(elevation)???CHECK 0082 0083 elApp = interp1(dcut.antenna0.tracker.utc, ... 0084 dcut.antenna0.tracker.horiz_topo(:,2), dcut.antenna0.receiver.utc); 0085 elOffSave = elApp - dcut.antenna0.servo.apparent(:,2); 0086 elOffSave = -elOffSave; 0087 angle = sqrt(elOffSave.^2 + azOffSave.^2); 0088 0089 %Detrend the data using data from X degrees away from source 0090 for chan = 1%chan_order % 0091 % Continue detrending channel by channel 0092 temp_data = dcut.antenna0.receiver.data(:,chan); 0093 temp_time = dcut.antenna0.receiver.utc; 0094 %if(strmatch('S',cbass,'exact')) 0095 temp_y = temp_data(angle>offset_angle); % Set to 7 for south data __> CHANGE once pointing fixed 0096 temp_x = temp_time(angle>offset_angle); 0097 0098 %else 0099 % temp_y = temp_data(angle>2); % Set to 2 for north data 0100 % temp_x = temp_time(angle>2); 0101 %end 0102 0103 % In case there are nans in south data from way we've done the flag kludge 0104 0105 %if(strmatch('S',cbass,'exact')) 0106 ignore = isnan(temp_y); 0107 p = polyfit(temp_x(~ignore),temp_y(~ignore),1); 0108 %else 0109 % p = polyfit(temp_x,temp_y,1); 0110 %end 0111 0112 all_time = temp_time; 0113 f = polyval(p,all_time); 0114 figure 0115 plot(temp_time,temp_data) 0116 hold all 0117 plot(temp_time,f) 0118 % Now subtract f from data channel on that scan 0119 dcut.antenna0.receiver.data(:,chan) = dcut.antenna0.receiver.data(:,chan) - f; 0120 %Keep a master copy too of all the data 0121 dgood.antenna0.receiver.data(ind,chan)=dcut.antenna0.receiver.data(:,chan); 0122 0123 0124 clear temp_data temp_time temp_x temp_y p f all_time 0125 end 0126 end 0127 clear dcut 0128 %% 0129 % Plot the detrended scan data 0130 close all 0131 figure(1) 0132 nchan=4 % d.data has [I1,Q1,U1,Q2,U2,Q3,U3,I2] 0133 plot_labels={'I1','I2','Q','U'} 0134 0135 points = (dgood.index.beammap.fast==1);% so we can identify the raster data 0136 0137 for i=1:nchan 0138 if i<7 0139 figure(1) 0140 subplot(2,3,i) 0141 0142 else 0143 figure(2) 0144 nplot=i-6; 0145 subplot(2,2,nplot) 0146 0147 end 0148 0149 plot(dgood.antenna0.receiver.utc(points),dgood.antenna0.receiver.data(points,chan_order(i))) 0150 xlabel(plot_labels{i}) 0151 gtitle('Pipeline calibrated and detrended data') 0152 end 0153 0154 %% Get offsets in Az and El for rasters 0155 dcut= framecut(dgood,dgood.index.beammap.fast); 0156 azApp = interp1(dcut.antenna0.tracker.utc, ... 0157 dcut.antenna0.tracker.horiz_topo(:,1), dcut.antenna0.receiver.utc); 0158 azOffSave = (azApp - dcut.antenna0.servo.apparent(:,1)); 0159 azOffSave = -wrap180(azOffSave); 0160 azOffSave = azOffSave.*cosd(dcut.antenna0.servo.el);% Scale azimuth by cos(elevation)???CHECK 0161 0162 elApp = interp1(dcut.antenna0.tracker.utc, ... 0163 dcut.antenna0.tracker.horiz_topo(:,2), dcut.antenna0.receiver.utc); 0164 elOffSave = elApp - dcut.antenna0.servo.apparent(:,2); 0165 elOffSave = -elOffSave; 0166 dcut.angle = sqrt(elOffSave.^2 + azOffSave.^2);% 0167 dcut.elOffSave = elOffSave; 0168 dcut.azOffSave= azOffSave; 0169 clear dgood 0170 0171 0172 0173 %% Plot the raster/data --> taken from raster analysis, so uses 2d binning 0174 0175 disp([ 'Plotting the raster: '] ) 0176 azOffRaster=dcut.azOffSave; 0177 elOffRaster=dcut.elOffSave; 0178 pixel_size = 0.1; % degree 0179 0180 azmin=-20.0 0181 azmax=20.0 0182 elmin=-20.0 0183 elmax=20.0 0184 %azmin = min(azOffRaster); 0185 %azmax = max(azOffRaster); 0186 %elmin = min(elOffRaster); 0187 %elmax = max(elOffRaster); 0188 colmin = 0 0189 colmax = 1e5 0190 0191 for i=1:4 0192 [TauAmap{i}, xidx, yidx, TauAmapIdx{i}] = bin_quickly2d(azmin, azmax, elmin, elmax, pixel_size,... 0193 azOffRaster, elOffRaster,dcut.antenna0.receiver.data(:,chan_order(i))); 0194 0195 % Uncomment if you want to see maps - old stuf related to raster maps 0196 % not az scans 0197 % % % Raster Maps 0198 % % figure(1) 0199 % % subplot(2,2,i) 0200 % % %clim=[colmin colmax] 0201 % % imagesc(fliplr(xidx),fliplr(yidx),TauAmap{i}) 0202 % % colormap('default') 0203 % % title(plot_labels{i}) 0204 % % axis equal 0205 % % ylim( [elmin elmax] ) 0206 % % xlim( [azmin azmax] ) 0207 % % colorbar 0208 % % 0209 end 0210 % % 0211 % % % %% Look at cuts through the raster maps 0212 % % disp( 'Look at cuts through the beam' ) 0213 % % % Get the user to select the centre of the beam pattern 0214 % % figure(5) 0215 % % climI = [0 1e5] 0216 % % colormap('default') 0217 % % imagesc(TauAmap{1}) 0218 % % title('Please click on the centre of the LL beam pattern') 0219 % % [usr_TrigAz1 usr_TrigEl1] = ginput(1); 0220 % % figure(5) 0221 % % imagesc(TauAmap{2}) 0222 % % title('Please click on the centre of the RR beam pattern') 0223 % % [usr_TrigAz2 usr_TrigEl2] = ginput(1); 0224 % % figure(5) 0225 % % imagesc(TauAmap{3}) 0226 % % title('Please click on the centre of the QQ beam pattern') 0227 % % [usr_TrigAz3 usr_TrigEl3] = ginput(1); 0228 % % figure(5) 0229 % % imagesc(TauAmap{4}) 0230 % % title('Please click on the centre of the UU beam pattern') 0231 % % [usr_TrigAz4 usr_TrigEl4] = ginput(1); 0232 % % 0233 % % usr_TrigAz = mean( [usr_TrigAz1 usr_TrigAz2 usr_TrigAz3 usr_TrigAz4] ); 0234 % % usr_TrigEl = mean( [usr_TrigEl1 usr_TrigEl2 usr_TrigEl3 usr_TrigEl4] ); 0235 % % 0236 % % usr_TrigAz = round( usr_TrigAz ); 0237 % % usr_TrigEl = round( usr_TrigEl ); 0238 % % 0239 % % % Plot the two cuts through the beam 0240 % % figure(5) 0241 % % subplot(2,4,1) 0242 % % plot((TauAmap{1}(usr_TrigEl,:))) 0243 % % ylabel( 'I1' ) 0244 % % xlabel( 'Azimuth' ) 0245 % % subplot(2,4,5) 0246 % % plot((TauAmap{1}(:,usr_TrigAz))) 0247 % % ylabel( 'I1' ) 0248 % % xlabel( 'Elevation' ) 0249 % % 0250 % % 0251 % % subplot(2,4,2) 0252 % % plot(TauAmap{2}(usr_TrigEl,:)) 0253 % % ylabel( 'I2' ) 0254 % % xlabel( 'Azimuth' ) 0255 % % subplot(2,4,6) 0256 % % plot(TauAmap{2}(:,usr_TrigAz)) 0257 % % ylabel( 'I2' ) 0258 % % xlabel( 'Elevation' ) 0259 % % 0260 % % subplot(2,4,3) 0261 % % plot(TauAmap{3}(usr_TrigEl,:)) 0262 % % ylabel( 'Q' ) 0263 % % xlabel( 'Azimuth' ) 0264 % % subplot(2,4,7) 0265 % % plot(TauAmap{3}(:,usr_TrigAz)) 0266 % % ylabel( 'Q' ) 0267 % % xlabel( 'Elevation' ) 0268 % % 0269 % % subplot(2,4,4) 0270 % % plot(TauAmap{4}(usr_TrigEl,:)) 0271 % % ylabel( 'U' ) 0272 % % xlabel( 'Azimuth' ) 0273 % % subplot(2,4,8) 0274 % % plot(TauAmap{4}(:,usr_TrigAz)) 0275 % % ylabel( 'U' ) 0276 % % xlabel( 'Elevation' ) 0277 0278 %% Plot the profile in azimuth 0279 close all 0280 0281 figure(4) 0282 subplot(2,1,1) 0283 plot(xidx,TauAmap{1}(201,:)) 0284 xlabel('Az Offset,deg') 0285 ylabel('Calibrated units, K') 0286 title('Calibrated scan data') 0287 0288 subplot(2,1,2) 0289 beam_min=min(min(TauAmap{1}(201,:))) 0290 beam=TauAmap{1}(201,:)+beam_min; 0291 xlabel('Az Offset,deg') 0292 ylabel('Amplitude, dBi') 0293 plot(xidx,10.*log10(beam/max(beam))) 0294 0295 %% Now read in a simulated beam 0296 0297 %addpath('./angela_cbass/angela_beam/sim_beams/') 0298 plot_cbass_beam_profiles 0299 0300 %% 0301 % Plot the simulated beam and measured moon beam on the same plot 0302 0303 figure(5) 0304 plot(xidx,10.*log10(beam/max(beam))) 0305 hold all 0306 plot(sim_angles, 10.*log10(sim_beam_cell_copol{f}./max(sim_beam_cell_copol{f}))) 0307 xlim([-20 20]) 0308 grid on 0309 xlabel('Az Offset, deg') 0310 ylabel('Amplitude, dBi') 0311 title('Simulated and measured co-pol beam') 0312 legend('Measured','Simulated') 0313 % % %% 0314 % % % One big figure of I to check sidelobes of Moon 0315 % % figure(30) 0316 % % imagesc(fliplr(xidx),fliplr(yidx),TauAmap{1}) 0317 % % 0318 % % xlabel('Az offset, deg'); 0319 % % ylabel('El offset, deg'); 0320 % % title(['Raster of ',num2str(start_date)]); 0321 % % axis equal 0322 % % % Save this map 0323 comment=source_name; 0324 map_name=(['Az_scan_Raster_',start_date,'_',comment,'.png']); 0325 print('-f5','-dpng',map_name) 0326 % % 0327 % % [X,Y]=meshgrid(xidx,yidx); 0328 % % [sf gof t] = Egauss2D_fit_act(X,Y,TauAmap{1,1}); 0329 % % 0330 % % fwhm_x = 2*sqrt(2*log(2)).*sf.sigma_x 0331 % % fwhm_y = 2*sqrt(2*log(2)).*sf.sigma_y 0332 % % amplitude=sf.a1 0333 % % x0=sf.x0 0334 % % y0=sf.y0 0335 % % theta=sf.theta 0336 % 0337 % 0338 % %% Write out the data to a file 0339 % filename = 'NorthRaster_Mar10.dat'; 0340 % fid = fopen(filename, 'a+'); 0341 % 0342 % data2save={start_date, comment, x0,y0,theta,fwhm_x,fwhm_y}; 0343 % fprintf( fid,'%s %s %d %d %d %d %d\n', data2save{1,:}) 0344 % 0345 % fclose(fid); 0346 % 0347 %% Save all the variables to a .mat file so we can reload at any point 0348 matname=[start_date,'_Az_Scan_',comment,'.mat']; 0349 save(matname,'xidx','yidx','TauAmap','d','dcut','TauAmap', 'azOffSave', 'elOffSave','sim_beam_cell_copol','sim_angles'); 0350 clear xidx,yidx,TauAmap,d,dcut,TauAmap, azOffSave, elOffSave,sim_beam_cell_copol,sim_angles 0351 close all 0352 0353 end 0354 0355