Raster analysis start_date = '08-Nov-2012:15:11:09' stop_date = '08-Nov-2012:20:21:09'
0001 function [source_name reduced_data_file nchan plot_labels] = Final_Raster_Analysis(save_path,start_date,stop_date,cbass,IntensityType,reduced_data_file) 0002 % Raster analysis 0003 %start_date = '08-Nov-2012:15:11:09' 0004 %stop_date = '08-Nov-2012:20:21:09' 0005 0006 %save_path ='/angela_cbass/angela_beam/new_rasters/fits/' % Adds this on to path of cbass_analysis 0007 %save_path ='/test_fits_June2013/' % Adds this on to path of cbass_analysis 0008 do_fits = 0 % write out complete data file to fits 0009 do_multi_fits = 0 % write out fits file for each complete raster 0010 do_contour = 1 % Do contour plots of individual rasters 0011 do_big_contour = 0 % make a contour plot of all the data together and fit a gaussian 0012 cell_size = 0.05 % grid cell size in degrees 0013 fit_gauss = 1 % Fit gaussian to each beam 0014 do_all_contours=0% Plot contours of I,Q,U and V (no parallactic angle rotation) 0015 do_shifted_contour = 1 0016 plot_results=1 0017 %08-Nov-2012:12:11:09 08-Nov-2012:20:21:52 /home/cbass/gcpCbass/control/sch/source_raster_8degM.sch( 8h , virgoa ) 0018 %% 0019 % Set up path and read in data 0020 0021 [home,installeddir] = where_am_i(); 0022 addpath([home,'/',installeddir,'/Angelas_Raster_Code/']); 0023 % addpath([home,'/',installeddir,'/angela_cbass/angela_beam/Beam_plotting/']); 0024 % addpath([home,'/',installeddir,'/angela_cbass/angela_beam/rasters/']); 0025 % addpath([home,'/',installeddir,'/angela_cbass/angela_beam/Tsys_stuff/']); 0026 % addpath([home,'/',installeddir,'/angela_cbass/angela_beam/reducs/']); 0027 maindir= [home,'/',installeddir,'/',save_path,'/']; 0028 0029 if(strmatch('S',cbass,'exact')) 0030 logfiles = [home,'/',installeddir,'/','log/obs_logSa.html']; % So we can check the schedule name in S logs %%HACKED AT PRESENT 0031 else 0032 logfiles = [home,'/',installeddir,'/','log/obs_log.html']; % So we can check the schedule name in N logs 0033 end 0034 0035 unix(['mkdir ',maindir]) % Make the directory in case it is not there already 0036 % Grab the schedule name - needed to specify which tiype of raster obs as 2 0037 % different schedules have been used. 0038 [start_dates,end_dates,schedule,source_name] =get_days_obs(logfiles,start_date); 0039 0040 disp(['day_reduce: FOUND ',int2str(length(start_dates)),' observations...']); 0041 0042 for i=1:length(start_dates) 0043 0044 sched_bits = regexp(char(schedule(i)), '/', 'split'); 0045 sched_name = char(sched_bits(length(sched_bits))); 0046 % strip out possible open bracket if there is one. 0047 sched_name = strrep(sched_name, '(',''); 0048 source_name = char(source_name); 0049 disp(['day_reduce: Schedule is ',sched_name]) 0050 disp(['Source name is ',source_name]) 0051 end 0052 %% 0053 0054 % Read in and calibrate data: CBASS-S 0055 0056 if(strmatch('S',cbass,'exact')) 0057 0058 disp(['Reducing a CBASS-S observation']) 0059 %d = read_arcSouth(start_date,stop_date); 0060 %d = determineIndicesSouth(d) 0061 %source_name = (char(dc.antenna0.tracker.source(100))); 0062 disp(['Reducing data for source: ',source_name]) 0063 0064 % Relabel the data in case different flag used 0065 0066 if (strfind(sched_name,'source_raster_8x4x0.05South.sch')) 0067 flagname = 'beammap'; 0068 d.index.radio_point_scan.fast = d.index.beammap.fast; 0069 d.index.radio_point_scan.slow = d.index.beammap.slow; 0070 else 0071 if (strfind(sched_name,'source_raster_20x10x0.1South.sch')) 0072 flagname = 'beammap'; 0073 d.index.radio_point_scan.fast = d.index.beammap.fast; 0074 d.index.radio_point_scan.slow = d.index.beammap.slow; 0075 else 0076 flagname = 'radio_point_scan' 0077 end 0078 end 0079 end 0080 %% 0081 % Read in data and calibrate: CBASS-N 0082 0083 if(strmatch('N',cbass,'exact')) 0084 disp(['Reducing a CBASS-N observation']) 0085 0086 %Make a separate entry point if data already reduced and pipeline calibrated 0087 0088 if (isempty(reduced_data_file)) 0089 disp('here') 0090 % Calibrate the data and extract I,Q,U,V 0091 d = read_arc(start_date,stop_date) 0092 % ACT - new reDscript implemented to use the new calibration as of 7/10/2013 0093 %d = reduceData(d,[home,'/',installeddir,'/Angelas_Raster_Code_South/redScript_beam_Nov12.red'], 0,'Raster_south_analysis')% results in [I,Q1,U1,Q2,U2,Q3,U3,I 0094 %d = reduceData(d,[home,'/',installeddir,'/reduc/redScript_raster.red'], 0,'Raster_analysis_newcal')% results in [I1,Q1,U1,Q2,U2,Q3,U3,I2 0095 %d = getIV(d,1) % Formulate I,Q,U,V for the filtered data ('1' option) [I Q1 U1 Q2 U2 Q3 U3 V] 0096 0097 % Ultra new calibration - don't bother with ground subtraction though as 0098 % rasters not just at elevation 37 or 47 0099 % 0100 d = reduceData(d,[home,'/',installeddir,'/Angelas_Raster_Code/redScript_newcal_raster_2014.red'], 0,save_path)% results in [I1,Q1,U1,Q2,U2,Q3,U3,I2] 0101 0102 fits_name=strcat(start_date,'_',source_name,'_',IntensityType) 0103 reduced_data_file = strcat([fits_name,'_pipelined.mat']) 0104 save([maindir,'/',reduced_data_file],'d','-v7.3') 0105 0106 else 0107 %disp(strcat(['Loading pre-reduced file: ',reduced_data_file'])) 0108 load(strcat(([maindir,reduced_data_file]))) 0109 end 0110 0111 % 0112 % Relabel the data in case different flag used 0113 0114 if (strfind(sched_name,'source_raster_8x4')) 0115 flagname = 'beammap'; 0116 d.index.radio_point_scan.fast = d.index.beammap.fast; 0117 d.index.radio_point_scan.slow = d.index.beammap.slow; 0118 0119 else 0120 flagname = 'radio_point_scan' 0121 end 0122 end 0123 0124 %% 0125 % Write out fits data for this whole data set 0126 if(do_fits==1) 0127 hdf = 0 0128 maindir= [home,'/',installeddir,'/',save_path,'/']; 0129 unix(['mkdir ',maindir]) % Just to check that the directory is there 0130 %fits_name='TauA_4Apr'; 0131 dc = cutObs(d, flagname); 0132 fits_name = strcat(start_date,'_',source_name,'_',IntensityType); 0133 %source_name = (char(dc.antenna0.tracker.source(100))); 0134 olddir = writeOutFits_raster_new(cbass,hdf,maindir, fits_name,dc,'all'); 0135 clear dc; 0136 end 0137 %% 0138 % Write out fits for each individual raster 0139 if(do_multi_fits==1) 0140 0141 [s e] = findStartStop(d.index.radio_point_scan.fast); % Find start and end points of each raster 0142 0143 for m=1:length(s) 0144 0145 % cut the data for each scan 0146 ind = zeros(size(d.index.radio_point_scan.fast)); 0147 ind(s(m):e(m)) = 1; 0148 ind = logical(ind); 0149 dc = framecut(d, ind); 0150 %subplot(5,5,m) 0151 %plot(dc.antenna0.receiver.data) 0152 hdf = 0; 0153 maindir= [home,'/',installeddir,'/',save_path,'/']; 0154 unix(['mkdir ',maindir]); % Just to check that the directory is there 0155 fits_name = strcat(start_date,'_',source_name,'_',IntensityType) 0156 number=num2str(m) 0157 olddir = writeOutFits_raster_new(cbass,hdf,maindir, fits_name,dc,number); 0158 end 0159 0160 end 0161 0162 %% 0163 % If C-BASS North then do the following data mash-up (historic reasons) 0164 if(strmatch('N',cbass,'exact')) 0165 %Apply the flagging so we can use flagged data from here on in in matlab 0166 %d = framecut(d,~d.flags.fast(:,1)) 0167 %d = applyFlags(d); % Sets data to nan where flagged 0168 0169 % % Form V explicity so we have it for later 0170 % 0171 d.antenna0.receiver.data(:,9)= (d.antenna0.receiver.data(:,8)-d.antenna0.receiver.data(:,1))/2; %V 0172 0173 nchan=9 % d.data has [I1,Q1,U1,Q2,U2,Q3,U3,I2,V] 0174 plot_labels = {'I1','Q1','U1','Q2','U2','Q3','U3','I2','V'}; 0175 0176 % % Old - now just keep all 8 channels 0177 % % % Re-order I Q U V into data(1,2,3,4) for ease later...should make this a 0178 % % % saner piece of code so don't have to do this! Relic of history... 0179 % % 0180 % % switch(IntensityType) 0181 % % case 'I' 0182 % % disp('Will analyse I,Q,U') 0183 % % d.antenna0.receiver.data(:,2)= d.antenna0.receiver.data(:,6); %Q3 0184 % % d.antenna0.receiver.data(:,3)= d.antenna0.receiver.data(:,7); %U3 0185 % % d.antenna0.receiver.data(:,4)= (d.antenna0.receiver.data(:,8)-d.antenna0.receiver.data(:,1))/2; %V 0186 % % d.antenna0.receiver.data(:,1)= (d.antenna0.receiver.data(:,1)+d.antenna0.receiver.data(:,8))/2; %I 0187 % % case 'I1' 0188 % % disp('Will analyse I1,Q1,U1') 0189 % % d.antenna0.receiver.data(:,2)= d.antenna0.receiver.data(:,2); %Q1 0190 % % d.antenna0.receiver.data(:,3)= d.antenna0.receiver.data(:,3); %U1 0191 % % d.antenna0.receiver.data(:,4)= (d.antenna0.receiver.data(:,8)-d.antenna0.receiver.data(:,1))/2; %V 0192 % % d.antenna0.receiver.data(:,1)= d.antenna0.receiver.data(:,1); %I1 0193 % % case 'I2' 0194 % % disp('Will analyse I2,Q2,U2') 0195 % % d.antenna0.receiver.data(:,2)= d.antenna0.receiver.data(:,4); %Q3 0196 % % d.antenna0.receiver.data(:,3)= d.antenna0.receiver.data(:,5); %U3 0197 % % d.antenna0.receiver.data(:,4)= (d.antenna0.receiver.data(:,8)-d.antenna0.receiver.data(:,1))/2;%V 0198 % % d.antenna0.receiver.data(:,1)= d.antenna0.receiver.data(:,8); %I2 0199 % % end 0200 0201 0202 end 0203 0204 %% 0205 % If C-BASS South then do the following data mash-up (hack to get going 0206 % with exising code) 0207 if(strmatch('S',cbass,'exact')) 0208 0209 0210 [d, nchan, plot_labels,chan_order] = getSDataready(d) ; 0211 0212 % Quick fix so I don;t edit all the code - make chan 4 RR (instead of 0213 % chan6) 0214 0215 0216 0217 % nchan=4; 0218 % LL = [d.antenna0.roach1.LL,fliplr(d.antenna0.roach2.LL)]; % **Check which order the roaches are in ** 0219 % RR = [d.antenna0.roach1.RR,fliplr(d.antenna0.roach2.RR)]; % **Check which order the roaches are in ** 0220 % QQ = [d.antenna0.roach1.Q,fliplr(d.antenna0.roach2.Q)]; 0221 % UU = [d.antenna0.roach1.U,fliplr(d.antenna0.roach2.U)]; 0222 % plot_labels = {'LL','RR','QQ','UU'}; 0223 % % Make a data structure d.antenna0.recevier.data(:,1:4) LL,RR,QQ,UU 0224 % % averaged over frequency 0225 % % saner piece of code so don't have to do this! Relic of history... 0226 % d.antenna0.receiver.data(:,1)= mean(LL'); %LL 0227 % d.antenna0.receiver.data(:,2)= mean(RR'); %RR 0228 % d.antenna0.receiver.data(:,3)= mean(QQ'); %QQ 0229 % d.antenna0.receiver.data(:,4)= mean(UU'); %UU 0230 % d.antenna0.receiver.utc = d.antenna0.roach1.utc; 0231 % 0232 % % Extra kludge to get apparent Az El --> this is wrong - just so we can 0233 % % check code for now. 0234 % d.antenna0.servo.apparent(:,1)=d.antenna0.servo.az; 0235 % d.antenna0.servo.apparent(:,2)=d.antenna0.servo.el; 0236 % 0237 % % Do some flagging - say when LL>1e4 units set all data to nan 0238 % flagarray = (abs(detrend(d.antenna0.receiver.data(:,1)))>1000); %%% HACK until we get some flagging routines - need to check if this level suits the data 0239 % for chan=1:nchan; 0240 % d.antenna0.receiver.data(flagarray,chan)=nan; 0241 % end 0242 end 0243 %% 0244 % For sanity, plot the data 0245 0246 for i=1:nchan 0247 if i<7 0248 figure(1) 0249 subplot(2,3,i) 0250 else 0251 figure(2) 0252 nchan 0253 nplot=i-6 0254 subplot(2,2,nplot) 0255 end 0256 0257 plot(d.antenna0.receiver.utc,d.antenna0.receiver.data(:,i)) 0258 xlabel(plot_labels{i}) 0259 end 0260 0261 if(strmatch('S',cbass,'exact')) 0262 gtitle('Noise-diode calibrated data for LL,Q,U,RR') 0263 else 0264 gtitle('Pipeline calibrated data') 0265 end 0266 %% 0267 % Detrend the data ready for analysis 0268 0269 0270 [s e] = findStartStop(d.index.radio_point_scan.fast); % Find start and end points of each raster 0271 dgood = d; % make a master copy which will later contain the detrended data (see in for loop below) 0272 n=0 0273 0274 0275 for m= 1:length(s) 0276 n = n+1; 0277 ind = zeros(size(d.index.radio_point_scan.fast)); 0278 ind(s(m):e(m)) = 1; 0279 ind = logical(ind); 0280 0281 if(strmatch('S',cbass,'exact')) 0282 dcut = framecutSa(d, ind); 0283 0284 else 0285 dcut = framecut(d, ind); 0286 end 0287 azApp = interp1(dcut.antenna0.tracker.utc, ... 0288 dcut.antenna0.tracker.horiz_topo(:,1), dcut.antenna0.receiver.utc); 0289 azOffSave = (azApp) - dcut.antenna0.servo.apparent(:,1); %%% Wrap180 put in for SOuth data 23/Aug/2013 - CHECK 0290 azOffSave = -azOffSave; 0291 0292 azOffSave = azOffSave.*cos(dcut.antenna0.servo.el*pi/180);% Scale azimuth by cos(elevation)???CHECK 0293 0294 elApp = interp1(dcut.antenna0.tracker.utc, ... 0295 dcut.antenna0.tracker.horiz_topo(:,2), dcut.antenna0.receiver.utc); 0296 elOffSave = elApp - dcut.antenna0.servo.apparent(:,2); 0297 elOffSave = -elOffSave; 0298 0299 %Detrend the data using data from 2 degrees away from source 0300 angle = sqrt(elOffSave.^2 + azOffSave.^2); 0301 0302 0303 for chan = 1:nchan % 0304 % Continue detrending channel by channel 0305 temp_data = dcut.antenna0.receiver.data(:,chan); 0306 temp_time = dcut.antenna0.receiver.utc; 0307 if(strmatch('S',cbass,'exact')) 0308 temp_y = temp_data(angle>7); % Set to 7 for south data __> CHANGE once pointing fixed 0309 temp_x = temp_time(angle>7); 0310 else 0311 temp_y = temp_data(angle>2); % Set to 2 for north data 0312 temp_x = temp_time(angle>2); 0313 end 0314 0315 % In case there are nans in south data from way we've done the flag kludge 0316 0317 if(strmatch('S',cbass,'exact')) 0318 ignore = isnan(temp_y); 0319 p = polyfit(temp_x(~ignore),temp_y(~ignore),1); 0320 else 0321 p = polyfit(temp_x,temp_y,1); 0322 end 0323 0324 all_time = temp_time; 0325 f = polyval(p,all_time); 0326 0327 % Now subtract f from data channel on that scan 0328 dcut.antenna0.receiver.data(:,chan) = dcut.antenna0.receiver.data(:,chan) - f; 0329 %Keep a master copy too of all the data 0330 dgood.antenna0.receiver.data(ind,chan)=dcut.antenna0.receiver.data(:,chan); 0331 0332 0333 clear temp_data temp_time temp_x temp_y p f all_time 0334 end 0335 end 0336 %% 0337 % Plot the detrended raster scan data 0338 close all 0339 figure(1) 0340 points = (dgood.index.radio_point_scan.fast==1);% so we can identify the raster data 0341 0342 for i=1:nchan 0343 if i<7 0344 figure(1) 0345 subplot(2,3,i) 0346 0347 else 0348 figure(2) 0349 nplot=i-6; 0350 subplot(2,2,nplot) 0351 0352 end 0353 0354 0355 0356 plot(dgood.antenna0.receiver.utc(points),dgood.antenna0.receiver.data(points,i)) 0357 xlabel(plot_labels{i}) 0358 gtitle('Pipeline calibrated and detrended data') 0359 end 0360 0361 % Save as a file 0362 if(exist('fits_name')) 0363 map_name = strcat([maindir,'/',fits_name,'_CalData.png']) 0364 %fig_name = strcat([maindir,'/',fits_name,'_Pointing.fig']) 0365 else 0366 fits_name = strcat(start_date,'_',source_name,'_',IntensityType) 0367 map_name = strcat([maindir,'/',fits_name,'_CalData.png']) 0368 %fig_name = strcat([maindir,'/',fits_name,'_Pointing.fig']) 0369 end 0370 %figure(1) 0371 %set(gcf, 'PaperPositionMode', 'auto'); 0372 %set(gcf,'PaperType','A4'); 0373 print('-f1','-dpng',map_name) 0374 if(nchan>6) 0375 map_name = [map_name,'_2']; 0376 print('-f2','-dpng',map_name) 0377 end 0378 clear dcut 0379 0380 % Save the dgood data for re-use later without having to load and 0381 % calibrate again 0382 0383 0384 save_file = strcat([maindir,'/',char(fits_name),'_dgood.mat']) 0385 save(save_file,'dgood','-v7.3') 0386 0387 0388 0389 close all 0390 end 0391 % %% End the function here - stuff below done elsewhere - only need _dgood%% 0392 % % % % % %% 0393 % % % % % Now make contour maps if you wish 0394 % if (do_contour==1) 0395 % n=0 0396 % % Find start and end points of each raster 0397 % [s e] = findStartStop(dgood.index.radio_point_scan.fast); 0398 % 0399 % %Identify the maximum of I,Q,U,V to use as a normalization for each raster 0400 % 0401 % points = (dgood.index.radio_point_scan.fast==1 | dgood.index.beammap.fast==1); 0402 % 0403 % contour_max(1) = max(dgood.antenna0.receiver.data(points,1)); 0404 % contour_max(2) = max(dgood.antenna0.receiver.data(points,2)); 0405 % contour_max(3) = max(dgood.antenna0.receiver.data(points,3)); 0406 % contour_max(4) = max(dgood.antenna0.receiver.data(points,4)); 0407 % 0408 % for m= 1:length(s) 0409 % n = n+1; 0410 % ind = zeros(size(dgood.index.radio_point_scan.fast)); 0411 % ind(s(m):e(m)) = 1; 0412 % ind = logical(ind); 0413 % dcut = framecut(dgood, ind); 0414 % azApp = interp1(dcut.antenna0.tracker.utc, ... 0415 % dcut.antenna0.tracker.horiz_topo(:,1), dcut.antenna0.receiver.utc); 0416 % azOffSave = azApp - dcut.antenna0.servo.apparent(:,1); 0417 % azOffSave = -azOffSave; 0418 % azOffSave = azOffSave.*cos(dcut.antenna0.servo.el*pi/180);% Scale azimuth by cos(elevation)???CHECK 0419 % 0420 % elApp = interp1(dcut.antenna0.tracker.utc, ... 0421 % dcut.antenna0.tracker.horiz_topo(:,2), dcut.antenna0.receiver.utc); 0422 % elOffSave = elApp - dcut.antenna0.servo.apparent(:,2); 0423 % elOffSave = -elOffSave; 0424 % angle = sqrt(elOffSave.^2 + azOffSave.^2); 0425 % 0426 % %Define a grid in degrees 0427 % ti=-10:cell_size:10; 0428 % 0429 % % Make sure we ignore NaNs 0430 % a = (isnan(azOffSave) | isnan(elOffSave)); 0431 % 0432 % for (chan=1) % only plot the I channel data for speed 0433 % 0434 % % Grid the data 0435 % [XI,YI]=meshgrid(ti,ti); 0436 % ZI = griddata(azOffSave(~a),elOffSave(~a),dcut.antenna0.receiver.data(~a,chan),XI,YI); 0437 % 0438 % % Grab the maximum 0439 % [num idx] = max(ZI(:)); 0440 % [nx ny] = ind2sub(size(ZI),idx); 0441 % max_ZI(m,1) = ZI(nx,ny); % Max value 0442 % max_ZI(m,2) = nx; % X-pos of max 0443 % max_ZI(m,3) = ny; % Y-pos of max 0444 % 0445 % %Check if there are multiple maxima 0446 % [num] = max(ZI(:)); 0447 % [nxs nys] = ind2sub(size(ZI),find(ZI==num)); 0448 % max_ZI(m,4) =length(nxs); 0449 % 0450 % % %Fit a gaussian to each one if you want 0451 % % if(chan==1) 0452 % % if(fit_gauss) 0453 % % [sf gof t] = gauss2D_fit_act(XI,YI,ZI); 0454 % % gfit{m} = sf; 0455 % % % If you want to plot the fit 0456 % % %plot(gfit{5},[XI(:),YI(:)],GoodI_c(:)) 0457 % % clear sf gof t 0458 % % end 0459 % % end 0460 % 0461 % % Plot mini contour maps on 3x3 grid on each page 0462 % fig_no = ceil(m/16); 0463 % subplot_no = m-(fig_no-1)*16; 0464 % 0465 % %subplot(nfigs_x,nfigs_y,n); 0466 % %Make contours based on the max_value in the data set for I,Q,U,V 0467 % VI = [1:-0.1:0]*contour_max(chan); 0468 % 0469 % figure(fig_no+2); 0470 % subplot(4,4,subplot_no) 0471 % contour(XI,YI,ZI,VI); 0472 % axis square; 0473 % grid on; 0474 % xlabel('Az Offset,deg'); 0475 % ylabel('El Offset,deg'); 0476 % title(num2str(m)); 0477 % 0478 % end 0479 % end 0480 % 0481 % 0482 % % % % % % Save first page of plots as an example file 0483 % % % % % if(exist('fits_name')) 0484 % % % % % map_name = strcat([maindir,'/',fits_name,'_Contour.png']) 0485 % % % % % %fig_name = strcat([maindir,'/',fits_name,'_Pointing.fig']) 0486 % % % % % else 0487 % % % % % fits_name = strcat(start_date,'_',char(dgood.antenna0.tracker.source(5))) 0488 % % % % % map_name = strcat([maindir,'/',fits_name,'_Contour.png']) 0489 % % % % % %fig_name = strcat([maindir,'/',fits_name,'_Pointing.fig']) 0490 % % % % % end 0491 % % % % % figure(1) 0492 % % % % % set(gcf, 'PaperPositionMode', 'auto'); 0493 % % % % % %set(gcf,'PaperType','A4'); 0494 % % % % % print('-f1','-dpng',map_name) 0495 % % % % % 0496 % % % % % % Save the gfit data for re-use later without having to load and 0497 % % % % % % calibrate again 0498 % % % % % 0499 % % % % % save_file = strcat([maindir,'/',fits_name,'_gfit.mat']) 0500 % % % % % 0501 % % % % % %save(save_file,'gfit') 0502 % % % % % %close all 0503 % % % % % end 0504 % % % % % %% 0505 % % % % % % Plot some information about the pointing - based on gaussian 2d fitting 0506 % % % % % % 0507 % % % % % %if(chan==1) 0508 % % % % % 0509 % % % % % if(fit_gauss) 0510 % % % % % for chan=1:4 0511 % % % % % figure(chan) 0512 % % % % % % Plot ratio of gaussian widths sigmax./sigmay 0513 % % % % % subplot( 2,2,1) 0514 % % % % % for chan=2:length(gfit) 0515 % % % % % i=1 0516 % % % % % % If only done channel 1 then i=1, chan = number of gfits 0517 % % % % % plot(i,((gfit{i,chan}.sigmax)*cell_size)/((gfit{i,chan}.sigmay)*cell_size),'.') 0518 % % % % % % plot(i,((gfit{i,chan}.sigmax)*cell_size)/((gfit{i}.sigmay)*cell_size),'.') 0519 % % % % % % hold all 0520 % % % % % xlabel('Raster number','fontsize',15) 0521 % % % % % ylabel('Ratio','fontsize',15) 0522 % % % % % title('Ratio of x- and y- best fit gaussian widths','fontsize',20) 0523 % % % % % end 0524 % % % % % 0525 % % % % % % Plot position of centre x-direction 0526 % % % % % subplot( 2,2,3) 0527 % % % % % for chan=2:length(gfit) 0528 % % % % % i=1 0529 % % % % plot(i,((gfit{i,chan}.x0)*cell_size),'.') 0530 % % % % hold all 0531 % % % % xlabel('Raster number','fontsize',15) 0532 % % % % ylabel('x-pos, deg','fontsize',15) 0533 % % % % title('x-Position of centre, deg','fontsize',20) 0534 % % % % end 0535 % % % % 0536 % % % % % Plot position of centre y-direction 0537 % % % % subplot( 2,2,4) 0538 % % % % for chan=2:length(gfit) 0539 % % % % i=1 0540 % % % % plot(i,((gfit{i,chan}.y0)*cell_size),'.') 0541 % % % % hold all 0542 % % % % xlabel('Raster number','fontsize',15) 0543 % % % % ylabel('y-pos, deg','fontsize',15) 0544 % % % % title('y-position of centre,deg','fontsize',20) 0545 % % % % end 0546 % % % % 0547 % % % % % Plot amplitude of maximum of gaussian 0548 % % % % subplot( 2,2,2) 0549 % % % % for chan=2:length(gfit) 0550 % % % % i=1 0551 % % % % plot(i,(gfit{i,chan}.a1),'.') 0552 % % % % hold all 0553 % % % % xlabel('Raster number','fontsize',15) 0554 % % % % ylabel('Amplitude','fontsize',15) 0555 % % % % title(['Amplitude of gaussian fit'],'fontsize',20) 0556 % % % % end 0557 % % % % gtitle('Results from 2D Gaussian fits',0.98,1) 0558 % % % % end 0559 % % % % % Save as a file 0560 % % % % if(exist('fits_name')) 0561 % % % % map_name = strcat([maindir,'/',fits_name,'_Pointing.png']) 0562 % % % % fig_name = strcat([maindir,'/',fits_name,'_Pointing.fig']) 0563 % % % % else 0564 % % % % fits_name = strcat(start_date,'_',char(dgood.antenna0.tracker.source(5))) 0565 % % % % map_name = strcat([maindir,'/',fits_name,'_Pointing.png']) 0566 % % % % %fig_name = strcat([maindir,'/',fits_name,'_Pointing.fig']) 0567 % % % % end 0568 % % % % figure(50) 0569 % % % % set(gcf, 'PaperPositionMode', 'auto'); 0570 % % % % %set(gcf,'PaperType','A4'); 0571 % % % % print('-f50','-dpng',map_name) 0572 % % % % 0573 % % % % % % Plot position of centre x- vs y-direction 0574 % % % % % figure 0575 % % % % % for i=2:length(gfit) 0576 % % % % % plot(((gfit{i}.x0)*cell_size),((gfit{i}.y0)*cell_size),'.') 0577 % % % % % hold all 0578 % % % % % title('xy-position of centre,deg','fontsize',20) 0579 % % % % % end 0580 % % % % 0581 % % % % end 0582 % % % % %end 0583 % % % % %end 0584 % % % % 0585 % % % % %% 0586 % % % % % Try doing one large contour 0587 % % % % if (do_big_contour==1) 0588 % % % % 0589 % % % % %Identify the rasters in the previously detrended data 0590 % % % % 0591 % % % % points = (dgood.index.radio_point_scan.fast==1); 0592 % % % % 0593 % % % % contour_max(1) = max(dgood.antenna0.receiver.data(points,1)); 0594 % % % % contour_max(2) = max(dgood.antenna0.receiver.data(points,2)); 0595 % % % % contour_max(3) = max(dgood.antenna0.receiver.data(points,3)); 0596 % % % % contour_max(4) = max(dgood.antenna0.receiver.data(points,4)); 0597 % % % % 0598 % % % % %for m= 1:length(s) 0599 % % % % % n = n+1; 0600 % % % % % ind = zeros(size(dgood.index.radio_point_scan.fast)); 0601 % % % % % ind(s(m):e(m)) = 1; 0602 % % % % % ind = logical(ind); 0603 % % % % dcut = framecut(dgood,points ); 0604 % % % % azApp = interp1(dcut.antenna0.tracker.utc, ... 0605 % % % % dcut.antenna0.tracker.horiz_topo(:,1), dcut.antenna0.receiver.utc); 0606 % % % % azOffSave = azApp - dcut.antenna0.servo.apparent(:,1); 0607 % % % % azOffSave = -azOffSave; 0608 % % % % azOffSave = azOffSave.*cos(dcut.antenna0.servo.el*pi/180);% Scale azimuth by cos(elevation)???CHECK 0609 % % % % 0610 % % % % elApp = interp1(dcut.antenna0.tracker.utc, ... 0611 % % % % dcut.antenna0.tracker.horiz_topo(:,2), dcut.antenna0.receiver.utc); 0612 % % % % elOffSave = elApp - dcut.antenna0.servo.apparent(:,2); 0613 % % % % elOffSave = -elOffSave; 0614 % % % % angle = sqrt(elOffSave.^2 + azOffSave.^2); 0615 % % % % 0616 % % % % %Define a grid in degrees 0617 % % % % ti=-6:cell_size:6; 0618 % % % % 0619 % % % % % Make sure we ignore NaNs 0620 % % % % a = (isnan(azOffSave) | isnan(elOffSave)); 0621 % % % % 0622 % % % % for (chan=1:4) % only plot the I channel data for speed 0623 % % % % 0624 % % % % % Grid the data 0625 % % % % [XI,YI]=meshgrid(ti,ti); 0626 % % % % ZI = griddata(azOffSave(~a),elOffSave(~a),dcut.antenna0.receiver.data(~a,chan),XI,YI); 0627 % % % % if(max(ZI<0)) 0628 % % % % ZI = -1*ZI; 0629 % % % % end 0630 % % % % 0631 % % % % % Grab the maximum 0632 % % % % [num idx] = max(ZI(:)); 0633 % % % % [nx ny] = ind2sub(size(ZI),idx); 0634 % % % % max_bigZI(1,1) = ZI(nx,ny); % Max value 0635 % % % % max_bigZI(1,2) = nx; % X-pos of max 0636 % % % % max_bigZI(1,3) = ny; % Y-pos of max 0637 % % % % 0638 % % % % %Check if there are multiple maxima 0639 % % % % [num] = max(ZI(:)); 0640 % % % % [nxs nys] = ind2sub(size(ZI),find(ZI==num)); 0641 % % % % max_bigZI(1,4) =length(nxs); 0642 % % % % 0643 % % % % %Fit a gaussian to each one if you want 0644 % % % % if(chan==1) 0645 % % % % if(fit_gauss) 0646 % % % % [sf gof t] = gauss2D_fit_act(XI,YI,ZI); 0647 % % % % big_gfit{1} = sf; 0648 % % % % % If you want to plot the fit 0649 % % % % %plot(gfit{5},[XI(:),YI(:)],GoodI_c(:)) 0650 % % % % clear sf gof t 0651 % % % % end 0652 % % % % end 0653 % % % % 0654 % % % % 0655 % % % % %subplot(nfigs_x,nfigs_y,n); 0656 % % % % %Make contours based on the max_value in the data set for I,Q,U,V 0657 % % % % VI = [1:-0.1:0]*contour_max(chan); 0658 % % % % 0659 % % % % figure(60); 0660 % % % % subplot(2,2,chan) 0661 % % % % contour(XI,YI,ZI,VI); 0662 % % % % axis square; 0663 % % % % grid on; 0664 % % % % xlabel('Az Offset,deg'); 0665 % % % % ylabel('El Offset,deg'); 0666 % % % % text(-2.5,2.7,['FWHM X: ', num2str(2*(sqrt(2*log(2))*big_gfit{1}.sigmax))]) 0667 % % % % text(-2.5,2.4,['FWHM Y: ', num2str(2*(sqrt(2*log(2))*big_gfit{1}.sigmay))]) 0668 % % % % title(num2str(m)); 0669 % % % % 0670 % % % % end 0671 % % % % 0672 % % % % 0673 % % % % 0674 % % % % % Save first page of plots as an example file 0675 % % % % if(exist('fits_name')) 0676 % % % % map_name = strcat([maindir,'/',fits_name,'_BigContour.png']) 0677 % % % % %fig_name = strcat([maindir,'/',fits_name,'_Pointing.fig']) 0678 % % % % else 0679 % % % % fits_name = strcat(start_date,'_',char(dgood.antenna0.tracker.source(5))) 0680 % % % % map_name = strcat([maindir,'/',fits_name,'_BigContour.png']) 0681 % % % % %fig_name = strcat([maindir,'/',fits_name,'_Pointing.fig']) 0682 % % % % end 0683 % % % % figure(60) 0684 % % % % set(gcf, 'PaperPositionMode', 'auto'); 0685 % % % % %set(gcf,'PaperType','A4'); 0686 % % % % print('-f60','-dpng',map_name) 0687 % % % % end 0688 % % % % 0689 % % % % %% 0690 % % % % % Having saved dgood, make big contour maps of the whole beam for I,Q,U,V 0691 % % % % % Try doing one large contour 0692 % % % % if (do_all_contours==1) 0693 % % % % 0694 % % % % %Identify the rasters in the previously detrended data 0695 % % % % 0696 % % % % points = (dgood.index.radio_point_scan.fast==1); 0697 % % % % 0698 % % % % contour_max(1) = max(dgood.antenna0.receiver.data(points,1)); 0699 % % % % contour_max(2) = max(dgood.antenna0.receiver.data(points,2)); 0700 % % % % contour_max(3) = max(dgood.antenna0.receiver.data(points,3)); 0701 % % % % contour_max(4) = max(dgood.antenna0.receiver.data(points,4)); 0702 % % % % 0703 % % % % %for m= 1:length(s) 0704 % % % % % n = n+1; 0705 % % % % % ind = zeros(size(dgood.index.radio_point_scan.fast)); 0706 % % % % % ind(s(m):e(m)) = 1; 0707 % % % % % ind = logical(ind); 0708 % % % % dcut = framecut(dgood,points ); 0709 % % % % azApp = interp1(dcut.antenna0.tracker.utc, ... 0710 % % % % dcut.antenna0.tracker.horiz_topo(:,1), dcut.antenna0.receiver.utc); 0711 % % % % azOffSave = azApp - dcut.antenna0.servo.apparent(:,1); 0712 % % % % azOffSave = -azOffSave; 0713 % % % % azOffSave = azOffSave.*cos(dcut.antenna0.servo.el*pi/180);% Scale azimuth by cos(elevation)???CHECK 0714 % % % % 0715 % % % % elApp = interp1(dcut.antenna0.tracker.utc, ... 0716 % % % % dcut.antenna0.tracker.horiz_topo(:,2), dcut.antenna0.receiver.utc); 0717 % % % % elOffSave = elApp - dcut.antenna0.servo.apparent(:,2); 0718 % % % % elOffSave = -elOffSave; 0719 % % % % angle = sqrt(elOffSave.^2 + azOffSave.^2); 0720 % % % % % Make sure we ignore NaNs 0721 % % % % a = (isnan(azOffSave) | isnan(elOffSave)); 0722 % % % % 0723 % % % % %Define a grid in degrees 0724 % % % % ti=-3:cell_size:3; 0725 % % % % 0726 % % % % 0727 % % % % pols=['I','Q','U','V'];% so we can label the plots 0728 % % % % for chan=1:4 % loop through all the channels 0729 % % % % 0730 % % % % 0731 % % % % % Grid the data 0732 % % % % [XI,YI]=meshgrid(ti,ti); 0733 % % % % ZI = griddata(azOffSave(~a),elOffSave(~a),dcut.antenna0.receiver.data(~a,chan),XI,YI); 0734 % % % % if(max(ZI<0)) 0735 % % % % ZI = -1*ZI; 0736 % % % % end 0737 % % % % 0738 % % % % % Grab the maximum 0739 % % % % [num idx] = max(ZI(:)); 0740 % % % % [nx ny] = ind2sub(size(ZI),idx); 0741 % % % % max_bigZI(1,1) = ZI(nx,ny); % Max value 0742 % % % % max_bigZI(1,2) = nx; % X-pos of max 0743 % % % % max_bigZI(1,3) = ny; % Y-pos of max 0744 % % % % 0745 % % % % %Check if there are multiple maxima 0746 % % % % [num] = max(ZI(:)); 0747 % % % % [nxs nys] = ind2sub(size(ZI),find(ZI==num)); 0748 % % % % max_bigZI(1,4) =length(nxs); 0749 % % % % 0750 % % % % %Fit a gaussian to each one if you want 0751 % % % % if(chan==1) 0752 % % % % if(fit_gauss==1) 0753 % % % % [sf gof t] = gauss2D_fit_act(XI,YI,ZI); 0754 % % % % big_gfit{1} = sf; 0755 % % % % % If you want to plot the fit 0756 % % % % %plot(gfit{5},[XI(:),YI(:)],GoodI_c(:)) 0757 % % % % clear sf gof t 0758 % % % % end 0759 % % % % end 0760 % % % % 0761 % % % % 0762 % % % % 0763 % % % % %subplot(nfigs_x,nfigs_y,n); 0764 % % % % %Make contours based on the max_value in the data set for I,Q,U,V 0765 % % % % %VI = [1:-0.1:0]*contour_max(chan); 0766 % % % % VI = [1:-0.1:0];% don't scale to data - this might be a problem if max value not at centre of map. 0767 % % % % 0768 % % % % figure; 0769 % % % % subplot(2,2,chan) 0770 % % % % contour(XI,YI,ZI); 0771 % % % % axis square; 0772 % % % % grid on; 0773 % % % % xlabel('Az Offset,deg'); 0774 % % % % ylabel('El Offset,deg'); 0775 % % % % %text(-2.5,2.7,['FWHM X: ', num2str(2*(sqrt(2*log(2))*big_gfit{1}.sigmax))]) 0776 % % % % %text(-2.5,2.4,['FWHM Y: ', num2str(2*(sqrt(2*log(2))*big_gfit{1}.sigmay))]) 0777 % % % % title(['Contour plot of ',pols(chan)]); 0778 % % % % 0779 % % % % if(plot_results==1) 0780 % % % % % Save first page of plots as an example file 0781 % % % % if(exist('fits_name')) 0782 % % % % save_path ='/test_fits/' 0783 % % % % [home,installeddir] = where_am_i(); 0784 % % % % maindir= [home,'/',installeddir,'/',save_path,'/']; 0785 % % % % map_name = strcat([maindir,'/',fits_name,'_',pols(chan),'_BigContour.png']) 0786 % % % % %fig_name = strcat([maindir,'/',fits_name,'_Pointing.fig']) 0787 % % % % else 0788 % % % % save_path ='/test_fits/' 0789 % % % % [home,installeddir] = where_am_i(); 0790 % % % % maindir= [home,'/',installeddir,'/',save_path,'/']; 0791 % % % % fits_name = strcat(start_date,'_',char(dgood.antenna0.tracker.source(5))) 0792 % % % % map_name = strcat([maindir,'/',fits_name,'_',pols(chan),'_BigContour.png']) 0793 % % % % %fig_name = strcat([maindir,'/',fits_name,'_Pointing.fig']) 0794 % % % % end 0795 % % % % 0796 % % % % set(gcf, 'PaperPositionMode', 'auto'); 0797 % % % % %set(gcf,'PaperType','A4'); 0798 % % % % print(gcf,'-dpng',map_name) 0799 % % % % end 0800 % % % % end 0801 % % % % 0802 % % % % end 0803 % % % % %% 0804 % % % % 0805 % % % % 0806 % % % % %Plot a shifted contour 0807 % % % % if(do_shifted_contour==1) 0808 % % % % [dgood gfit ZIsp XIsp YIsp B] = plotShiftedRaster(dgood,gfit,0.1,plot_results,start_date,stop_date); 0809 % % % % end