FitsRADECMap

 script to plot data taken over a time series on the Northern RX 
 inputs - start_date, stop_date

0001 function FitsRADECMap(fitsfilepath,listoffiles,do_pa_rot)
0002 % script to plot data taken over a time series on the Northern RX
0003 % inputs - start_date, stop_date
0004 
0005 % ACT 10 July - adapted LongTermTimeSeries to read in Fits files instead
0006 
0007 % takes the following inputs:
0008 % listof files - e.g. files.txt file containing all the fits files you want
0009 % to look at
0010 %
0011 %  do_pa_rot - set to 1 if you want to rotate the data
0012 %%
0013 % Setup some paths
0014 [home,installeddir] = where_am_i();
0015 addpath([home,'/',installeddir,'/Angelas_Raster_Code/']);
0016 % addpath([home,'/',installeddir,'/angela_cbass/angela_beam/Beam_plotting/']);
0017 % addpath([home,'/',installeddir,'/angela_cbass/angela_beam/rasters/']);
0018 % addpath([home,'/',installeddir,'/angela_cbass/angela_beam/Tsys_stuff/']);
0019 % addpath([home,'/',installeddir,'/angela_cbass/angela_beam/reducs/']);
0020 maindir= [home,'/',installeddir,'/'];
0021 
0022 % % %TauA - 17-Jul
0023 %  fitsfilepath='/reduc/raster_tests/17-Jul-2013:13:54:27/fits/'
0024 %  calfile = [maindir,'/',fitsfilepath,'17-Jul-2013:13:53:34_taua_I_caldata.mat'];
0025 % %calfile = [maindir,'/',fitsfilepath,'16-Jul-2013:13:48:20_m42_I_caldata.mat'];
0026 %  source_name = 'TauA'
0027 % listoffiles='fits_files_tomap.txt'
0028 %%%%%%
0029 % % %TauA - 6-Feb-2014
0030 %   fitsfilepath='/reduc/n_raster_elliptical_newcal_apr14/06-Feb-2014:01:25:09/fits/'
0031 %   calfile = [maindir,'/',fitsfilepath,'06-Feb-2014:01:24:02_taua_I_caldata.mat'];
0032 %   listoffiles='fits_files_tomap.txt'
0033 %   source_name = 'TauA'
0034 % %%%%%%
0035 
0036 % %%%%%
0037 % %M42 - 13-Dec-2013
0038 % fitsfilepath='/reduc/n_raster_elliptical_newcal_apr14/13-Dec-2013:08:40:31/fits/'
0039 % calfile = [maindir,'/',fitsfilepath,'13-Dec-2013:08:39:09_m42_I_caldata.mat'];
0040 % source_name = 'M42'
0041 % listoffiles='fits_files_tomap.txt'
0042 % %%%%%
0043 
0044 %%%%
0045 % %M42 - 22-Sep-2013
0046 % fitsfilepath='/reduc/RasterAug2014_rfi20_detrend/';
0047 % %calfile = [maindir,'/',fitsfilepath,'22-Sep-2013:10:48:51_m42_I_caldata.mat'];
0048 % source_name = 'M42'
0049 % listoffiles='m42_filestomap.txt'
0050 %TauA - 22-Sep-2013
0051 fitsfilepath='/reduc/RasterAug2014_rfi20_detrend/';
0052 %calfile = [maindir,'/',fitsfilepath,'22-Sep-2013:10:48:51_m42_I_caldata.mat'];
0053 source_name = 'TauA'
0054 listoffiles='files.txt'
0055 %%%%%
0056 %M42
0057 %source_name = 'M42'
0058 %fitsfilepath='/reduc/raster_tests/16-Jul-2013:13:49:14/fits/';
0059 %calfile = [maindir,'/',fitsfilepath,'16-Jul-2013:13:48:20_m42_I_caldata.mat'];
0060 %listoffiles='files.txt'
0061 %
0062 %load(calfile)
0063 %% Read in list of files
0064 %  Format string for each line of text:
0065 %   column1: text (%s)
0066 %   For more information, see the TEXTSCAN documentation.
0067 formatSpec = '%s%[^\n\r]';
0068 delimiter = ' ';
0069 % Open the text file.
0070 fileID = fopen([maindir,'/',fitsfilepath,'/',listoffiles],'r');
0071 
0072 % Read columns of data according to format string.
0073 % This call is based on the structure of the file used to generate this
0074 % code. If an error occurs for a different file, try regenerating the code
0075 % from the Import Tool.
0076 dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter,  'ReturnOnError', false);
0077 
0078 % Close the text file.
0079 fclose(fileID);
0080 
0081 
0082 % Allocate imported array to column variable names
0083 filenames = dataArray{:, 1};
0084 
0085 % Clear temporary variables
0086 clearvars listoffiles delimiter formatSpec fileID dataArray ans;
0087 
0088 %%
0089 
0090 % Initialise arrays
0091 MJD  = [];
0092 RA   = [];
0093 DEC  = [];
0094 EL_D = [];
0095 AZ_D = [];
0096 I1    = [];
0097 Q1    = [];
0098 U1    = [];
0099 I2    = [];
0100 Q2    = [];
0101 U2    = [];
0102 %V= (data{str2num(cell2mat(data_pos(strcmp(data_labels,'V'))))});
0103 %GLAT_D= [];
0104 %GLON_D= [];
0105 DAYFLAG = logical([]);
0106 FLAG = logical([]);
0107 
0108 %Now read in one file at a time
0109 figure
0110 hold all
0111 for i=1:length(filenames)
0112 data_labels =  { 'MJD', 'RA', 'DEC', 'AZ', 'EL', 'I1', 'Q1', 'U1','I2', 'Q2', 'U2','DAYFLAG','FLAG'};
0113 info =fitsinfo([maindir,'/',fitsfilepath,cell2mat(filenames(i))]);
0114 data = fitsread([maindir,'/',fitsfilepath,cell2mat(filenames(i))],'bintable',1);
0115 values = info.BinaryTable(1).Keywords(:,2);
0116 names = info.BinaryTable(1).Keywords(:,1);
0117 pinfo = info.PrimaryData.Keywords(:,1:2) ;
0118 
0119 % Find the position of the data for each of the data_labels above
0120 
0121 for ilab=1:length(data_labels);
0122     ilab;
0123     pos = (names(strcmp(data_labels{ilab},values)));
0124     pos = regexp(pos{1},'\d','match'); %Find numbers in the string
0125     pos = cell2mat(pos);
0126     data_pos{ilab}=pos;
0127 end
0128 
0129    
0130 
0131 
0132 %Now read in the data from each file
0133 
0134 tMJD= (data{str2num(cell2mat(data_pos(strcmp(data_labels,'MJD'))))});
0135 tRA= (data{str2num(cell2mat(data_pos(strcmp(data_labels,'RA'))))});
0136 tDEC = (data{str2num(cell2mat(data_pos(strcmp(data_labels,'DEC'))))});
0137 tEL_D = (data{str2num(cell2mat(data_pos(strcmp(data_labels,'EL'))))}); % in degrees
0138 tAZ_D = (data{str2num(cell2mat(data_pos(strcmp(data_labels,'AZ'))))}); % in degrees
0139 tI1 = ((data{str2num(cell2mat(data_pos(strcmp(data_labels,'I1'))))}));
0140 tQ1 = ((data{str2num(cell2mat(data_pos(strcmp(data_labels,'Q1'))))}));
0141 tU1 = ((data{str2num(cell2mat(data_pos(strcmp(data_labels,'U1'))))}));
0142 tI2 = ((data{str2num(cell2mat(data_pos(strcmp(data_labels,'I2'))))}));
0143 tQ2 = ((data{str2num(cell2mat(data_pos(strcmp(data_labels,'Q2'))))}));
0144 tU2 = ((data{str2num(cell2mat(data_pos(strcmp(data_labels,'U2'))))}));
0145 %V= (data{str2num(cell2mat(data_pos(strcmp(data_labels,'V'))))});
0146 %tGLAT_D= rad2deg(data{str2num(cell2mat(data_pos(strcmp(data_labels,'LAT'))))}); % Galactic Latitude
0147 %tGLON_D= rad2deg(data{str2num(cell2mat(data_pos(strcmp(data_labels,'LON'))))}); % Galactic Longitude
0148 tDAYFLAG = logical(data{str2num(cell2mat(data_pos(strcmp(data_labels,'DAYFLAG'))))}); % Day flag
0149 tFLAG = logical(data{str2num(cell2mat(data_pos(strcmp(data_labels,'FLAG'))))}); % Pipeline flag
0150 
0151 % % Detrend I Q and U
0152 %
0153 % galcut = 30 % cut on galactic latitude in degrees
0154 % flag_vector = (~tFLAG & (abs(tGLAT_D)>galcut));
0155 % order = 1 % linear fit
0156 %
0157 % tI = flagdetrend(tMJD, tI,flag_vector,order);
0158 % tQ = flagdetrend(tMJD, tQ,flag_vector,order);
0159 % tU = flagdetrend(tMJD, tU,flag_vector,order);
0160   
0161 % Concatenate detrended data from all surveys together
0162 MJD  = [MJD;tMJD];
0163 RA   = [RA;tRA];
0164 DEC  = [DEC;tDEC];
0165 EL_D = [EL_D;tEL_D];
0166 AZ_D = [AZ_D;tAZ_D];
0167 I1    = [I1;tI1];
0168 Q1    = [Q1;tQ1];
0169 U1    = [U1;tU1];
0170 I2    = [I2;tI2];
0171 Q2    = [Q2;tQ2];
0172 U2    = [U2;tU2];
0173 %V= (data{str2num(cell2mat(data_pos(strcmp(data_labels,'V'))))});
0174 %GLAT_D= [GLAT_D;tGLAT_D];
0175 %GLON_D= [GLON_D;tGLON_D];
0176 DAYFLAG = logical([DAYFLAG;tDAYFLAG]);
0177 FLAG = logical([FLAG;tFLAG]);
0178 
0179 % Put RA and DEC in degrees
0180 
0181 DEC = rad2deg(DEC);
0182 RA = rad2deg(RA);
0183 % Plot each data set as we go along to look for weird things!
0184 % figure(1)
0185 % plot((data{str2num(cell2mat(data_pos(strcmp(data_labels,'MJD'))))}),((data{str2num(cell2mat(data_pos(strcmp(data_labels,'I1'))))})))
0186 % hold all
0187 end
0188 %%
0189 % make I vector and Q and U vectors so we can look at all data
0190 Iall = [I1,I2,(I1+I2)/2];
0191 Qall = [Q1,Q2,(Q1+Q2)/2];
0192 Uall = [U1,U2,(U1+U2)/2];
0193 Pall = [sqrt(Qall(:,1).^2 + Uall(:,1).^2),sqrt(Qall(:,2).^2 + Uall(:,2).^2),sqrt(Qall(:,3).^2 + Uall(:,3).^2)];
0194 allIQUP = {Iall,Qall,Uall,Pall};
0195 allIQUP_labels={'I','Q','U','P'};
0196 figure(2)
0197 for i=1:4 % i gives I,Q,U,P
0198     subplot(4,1,i)
0199     for j=1:3 % j gives I1,I2,I
0200         plot(allIQUP{1,i}(:,j))
0201         hold all
0202     end
0203     xlabel('MJD')
0204     ylabel('Kelvin')
0205     title(allIQUP_labels{i})
0206 end
0207 gtitle('Timestream data for I,Q,U,P  - no PA rotation')
0208 
0209 %%
0210 %%%%%
0211 %%%NEED TO CHANGE THIS FOR DIFFERENT SIZE CAL MATRIX and use inverse from
0212 %%%what is used in descart
0213 if(exist('calmat'))
0214     calmat % Order is I1,I2,Q1,Q2,U1,U2
0215     
0216     calmat1=calmat(:,[1,3,5]); % calibration matrix for I1,Q1,U1
0217     calmat2=calmat(:,[2,4,6]); % calibration matrix for I2,Q2,U2
0218     
0219     invcalmat1 = inv(calmat1); % Should calibrate raw observed data --> corrected data
0220     invcalmat2 = inv(calmat2);
0221     
0222     data1 = [Iall(:,1),Qall(:,1),Uall(:,1)];
0223     data2 = [Iall(:,2),Qall(:,2),Uall(:,2)];
0224     
0225     % Apply the inverse matrices to dataP (I1,Q1,U1) and dataP(I2,Q2,U2)
0226     % dataP is [I1,Q1,U1,Q2,U2,Q3,U3,I2,V]
0227     for i=1:length(data1(:,1));
0228         data1(i,[1,2,3])=data1(i,[1,2,3])*invcalmat1;  % Mulitply I,Q,U by the calibration matrix
0229         data2(i,[1,2,3])=data2(i,[1,2,3])*invcalmat2;
0230     end
0231     
0232 disp('Calibration now applied')
0233 
0234 Qall = [data1(:,2),data2(:,2),(data1(:,2)+data2(:,2))/2];
0235 Uall = [data1(:,3),data2(:,3),(data1(:,3)+data2(:,3))/2] ;
0236 Pall = [sqrt(Qall(:,1).^2 + Uall(:,1).^2),sqrt(Qall(:,2).^2 + Uall(:,2).^2),sqrt(Qall(:,3).^2 + Uall(:,3).^2)];
0237 
0238 else
0239     disp('No calibration applied')
0240 end
0241     
0242 
0243 
0244 %%%
0245 
0246 %%
0247 %
0248 % Rotate parallactic angle of polarization data
0249 % Grab latitude and longitude of antenna
0250     
0251 % [antenna_no antenna_name] = identifyTelescope(d);
0252 % telescope_constants
0253 %
0254 % LAT_D = antenna_latitude(antenna_no)
0255 % LON_D = antenna_longitude(antenna_no)
0256 %
0257 
0258 LAT_D = cell2mat(pinfo(strcmp(pinfo(:,1),'LAT-OBS'),2)) % latitude degrees
0259 LON_D = cell2mat(pinfo(strcmp(pinfo(:,1),'LONG-OBS'),2)) % longitude degrees
0260 
0261 LAT = deg2rad(LAT_D); % radians
0262 LON = deg2rad(LON_D); % radians
0263 
0264 JD = mjd2jd(MJD);
0265 LST=lst(JD,deg2rad(LON_D),'m'); % Uses East longitude in radians
0266                                 % Gives LST in fraction of days
0267     
0268 sinPA = -((sind(AZ_D).*cosd(LAT_D))./cosd(DEC)); % degrees
0269 sinPA(abs(sinPA)>=1)=1;
0270 PA_D = asind(sinPA); %Parallactic angle in degrees
0271 
0272 
0273 % Rotate Q and U and calculate P
0274 
0275 for i=1:3 % j gives I1,I2,I etc
0276         Qall_pa(:,i) = Qall(:,i).*cosd(2*PA_D) - Uall(:,i).*sind(2*PA_D) ;
0277         Uall_pa(:,i) = Uall(:,i).*cosd(2*PA_D) + Qall(:,i).*sind(2*PA_D) ;
0278         Pall_pa(:,i) = (sqrt(Qall_pa(:,i).^2 + Uall_pa(:,i).^2));
0279 end
0280 
0281 
0282 %%
0283 % Check the data for each raster and flag out dodgy rasters
0284 close all
0285 flagdodgyQ=interactive_flag(Q1);
0286 nanflag = isnan(I1);
0287 flagvector2 = (~nanflag & ~FLAG & ~flagdodgyQ(1:size(Q1,1)));
0288 
0289 % Plot parallactic rotated data  - ask to identify any rasters to flag
0290 if(do_pa_rot==1)
0291     allIQUP_pa = {Iall,Qall_pa,Uall_pa,Pall_pa};
0292 else
0293     disp('No parallactic angle rotation applied')
0294     allIQUP_pa = {Iall,Qall,Uall,Pall};
0295 end
0296 
0297 figure
0298 for i=1:4 % i gives I,Q,U,P
0299     subplot(4,1,i)
0300     for j=1:3 % j gives I1,I2,I
0301         plot(allIQUP_pa{1,i}(flagvector2,j))
0302         %plot(allIQUP_pa{1,i}(flagvector2,j))
0303         hold all
0304     end
0305     xlabel('MJD')
0306     ylabel('Kelvin')
0307     title(allIQUP_labels{i})
0308 end
0309 if(exist('calmat'))
0310 gtitle('Calibrated timestream data for I,Q,U,P  - with PA rotation')
0311 else
0312    gtitle('Timestream data for I,Q,U,P  - with PA rotation')
0313 end
0314 
0315 raster_flags = input('List the rasters you wish to flag e.g. [1 2 3]: ');
0316 if isempty(raster_flags)
0317     raster_flags = [];
0318 end
0319 
0320 % Set that flag (nanflag) to '1' in dodgy rasters and reset flagvector 2
0321 tsig =(abs(diff(EL_D)));
0322 dsig=(tsig<=1);
0323 [s e]=findStartStop(dsig);
0324 
0325 %Now plot out the I1 data for each raster as a check
0326 n=0
0327 for m=raster_flags
0328     for i=1:4
0329          n = n+1;
0330          ind = zeros(size(FLAG));
0331          ind(s(m):e(m)) = 1;
0332          ind = logical(ind);
0333          nanflag(ind)=1;
0334     end
0335     
0336 end
0337 
0338 flagvector2 = (~nanflag & ~FLAG );
0339 figure
0340 for i=1:4 % i gives I,Q,U,P
0341     subplot(4,1,i)
0342     for j=1:3 % j gives I1,I2,I
0343         plot(allIQUP_pa{1,i}(flagvector2,j),'.')
0344         hold all
0345     end
0346     xlabel('MJD')
0347     ylabel('Kelvin')
0348     title(allIQUP_labels{i})
0349 end
0350 
0351 % Plot I1,Q1,U1  and I2,Q2,U2 --> For raster plots memo
0352 figure;
0353 plot_labels={'I1','Q1','U1','P1','I2','Q2','U2','P2'};
0354 plot_lims=[-0.5 3.6; -0.08 0.08; -0.05 0.3; 0 0.25];
0355 k=1;
0356 tdate=(mjd2date_v2(MJD));% get in date form
0357 aa  = datenum(tdate.year,tdate.month,tdate.day,tdate.hour,tdate.minute,tdate.second);
0358 
0359 for j=1:2 % i gives I1,I2
0360     
0361     for i=1:4 % j gives I,Q,U,P
0362         subplot(2,4,k)
0363         %plot(aa(flagvector2),smooth(allIQUP_pa{1,i}(flagvector2,j),10),'.')
0364         plot(smooth(allIQUP_pa{1,i}(flagvector2,j),10))
0365         xlabel('Time hours')
0366         ylabel('Kelvin')
0367         xlim([min(aa) max(aa)]);
0368         ylim([plot_lims(i,1) plot_lims(i,2)])
0369         datetick('x','hh','keeplimits','keepticks');
0370         title(plot_labels(k));
0371         k=k+1;
0372     end
0373 end
0374 %%
0375 %%%% Make a map of just TauA region
0376 %close all
0377 
0378 %source_name = 'TauA';
0379 [RA_tau DEC_tau] = sourcePos(source_name) % deg
0380 RA_off = RA-RA_tau;
0381 DEC_off = DEC-DEC_tau;
0382 RA_off = RA_off.*(cosd(DEC)); % scale RA to declination
0383 pixel_size = 0.1;
0384 % RA = rad2deg(RA);
0385 % DEC= rad2deg(DEC);
0386 % ramin = RA_tau - 3;
0387 % ramax = RA_tau + 3;
0388 % decmin = DEC_tau - 3;
0389 % decmax = DEC_tau + 3;
0390 ramin = - 2;
0391 ramax = + 2;
0392 decmin = - 2;
0393 decmax = + 2;
0394 
0395 figure
0396 map_labels={'I','Q','U','P'};
0397 for i=1:4
0398     [TauAmap{i}, xidx, yidx, TauAmapIdx{i}] = bin_quickly2d(ramin, ramax, decmin, decmax, pixel_size,...
0399         RA_off(flagvector2), DEC_off(flagvector2),allIQUP_pa{1,i}(flagvector2,3));
0400         subplot(2,2,i)
0401     imagesc(fliplr(xidx),yidx,TauAmap{i})
0402     title([source_name,' :',map_labels(i)])
0403     xlabel('RA offset, deg')
0404     ylabel('DEC offset, deg')
0405     colormap('default')
0406     colorbar
0407     axis square
0408 end
0409 
0410 %%
0411 % To plot each raster individually
0412 % Identify start and stop of rasters via distinct change in elevation
0413 tsig =(abs(diff(EL_D)));
0414 dsig=(tsig<=1);
0415 [s e]=findStartStop(dsig);
0416 
0417 % Now make mini raster map for each one RA DEC
0418 
0419 close all
0420 
0421 [RA_tau DEC_tau] = sourcePos(source_name); % deg
0422 
0423 RA_off = RA-RA_tau;
0424 DEC_off = DEC-DEC_tau;
0425 RA_off = RA_off.*(cosd(DEC)); % scale RA to declination
0426 
0427 pixel_size = 0.1;
0428 
0429 ramin = - 5;
0430 ramax = + 5;
0431 decmin = - 5;
0432 decmax = + 5;
0433 
0434 % Set colour scale mapping
0435 clim = [0 3.5;-0.02 0.02;-0.02 0.02;0 0.02]
0436 
0437 for m= 1:length(s)
0438     n = n+1;
0439     ind = zeros(size(FLAG));
0440     ind(s(m):e(m)) = 1;
0441     ind = logical(ind);
0442     flagvector3 = (~nanflag & ~FLAG & ind & ~flagdodgyQ(1:size(Q1,1)));
0443     
0444     for i=1:4
0445         [TauAmap{i}, xidx, yidx, TauAmapIdx{i}] = bin_quickly2d(ramin, ramax, decmin, decmax, pixel_size,...
0446             RA_off(flagvector3), DEC_off(flagvector3),allIQUP_pa{1,i}(flagvector3,3));
0447         figure(m)
0448         subplot(2,2,i)
0449         
0450         %imagesc(fliplr(xidx),yidx,TauAmap{i},clim(i,:))
0451         imagesc(fliplr(xidx),yidx,TauAmap{i})
0452         title([source_name,' map:', num2str(i)])
0453         colormap('default')
0454         colorbar
0455         axis square
0456     end
0457     
0458     
0459     
0460 end
0461 %col = jet(64); tmp = linspace(0,1,64)';
0462 %for n = 1:3, col(:,n) = interp1( 10.^tmp, col(:,n), 1+9*tmp, 'linear'); end
0463 %colormap(col)
0464 
0465 
0466 %%
0467 % To make binned maps in Az and El
0468 % Uses code in /usersVol1/act/CBASS/cbass_analysis/matutils/angles
0469 % Need to calculate the Az and El of the source given DEC and RA
0470 gfit=[];
0471 % Need JD
0472 [RA_tau DEC_tau] = sourcePos(source_name); % de
0473 JD = mjd2jd(MJD);
0474 
0475 % Now need mean sidereal time
0476 LST=lst(JD,deg2rad(LON_D),'m'); % Uses East longitude in radians
0477                                  % Gives LST in fraction of days
0478 %Now need hour angle
0479 %ha = ra2ha(ra,time); % converts ra in degrees to hour angle in HOURS
0480 %time is lst in frac HOURS !!
0481 
0482 LST_hr = LST*24;
0483 
0484 ha_deg = ra2ha(RA_tau,LST_hr)*15; 
0485 ha_rad = deg2rad(ha_deg);
0486 % Finally calculate Az and El of source
0487 [A,E]=hdl2ae(ha_rad,deg2rad(DEC_tau),deg2rad(LAT_D)) ;% all needs to be in radians
0488 
0489 
0490 
0491 Az_source = rad2deg(A);
0492 Az_source = wrap360(Az_source); % So we arap consistetly with Az drive
0493 El_source = rad2deg(E);
0494 % figure
0495 % plot(Az_source);
0496 % hold all
0497 % plot(AZ_D)
0498 
0499 % Now we calculate the Az_off and El_off from the source
0500 
0501 Azoff = AZ_D - Az_source;
0502 Eloff = EL_D - El_source;
0503 
0504 % Scale Azimuth for Elevation
0505 
0506 Azoff = Azoff.*cosd(EL_D);
0507 
0508 % Do the binning
0509 
0510 tsig =(abs(diff(EL_D)));
0511 dsig=(tsig<=1);
0512 [s e]=findStartStop(dsig);
0513 
0514 pixel_size = 0.1; % degree
0515 
0516 azmin = - 5;
0517 azmax = + 5;
0518 elmin = - 5;
0519 elmax = + 5;
0520 
0521 % Set colour scale mapping
0522 clim = [0 3.5;-0.02 0.02;-0.02 0.02;0 0.04]
0523 n=0
0524 map_labels={'I','Q','U','P'};
0525 
0526 for m= 1:length(s)
0527     n = n+1;
0528     ind = zeros(size(FLAG));
0529     ind(s(m):e(m)) = 1;
0530     ind = logical(ind);
0531     flagvector3 = (~nanflag & ~FLAG & ind & ~flagdodgyQ(1:size(Q1,1)));
0532     
0533     for i=1:4
0534         for j=1:3
0535         [TauAmap{i,j}, xidx, yidx, TauAmapIdx{i,j}] = bin_quickly2d(azmin, azmax, elmin, elmax, pixel_size,...
0536             Azoff(flagvector3), Eloff(flagvector3),allIQUP_pa{1,i}(flagvector3,j));
0537         figure(m)
0538         subplot(2,2,i)
0539         
0540         %imagesc(fliplr(xidx),yidx,TauAmap{i},clim(i,:))
0541         imagesc(fliplr(xidx),yidx,TauAmap{i})
0542         title([source_name,' :',map_labels(i)])
0543         xlabel('AZ offset, deg')
0544         ylabel('EL offset, deg')
0545         colormap('default')
0546         colorbar
0547         axis square
0548     %end
0549 
0550 %
0551 %Do Gaussian fits to these beam plots
0552 
0553 %for chan=1:4
0554     [sf gof t] = Egauss2D_fit_act(XI,YI,TauAmap{i,j},[],i,gfit,m);
0555     gfit{m,i,j} = sf;
0556 end
0557     end  
0558 end
0559 %%
0560 % Make a single map in Az El
0561 figure
0562 for i=1:4
0563         [TauAmap{i}, xidx, yidx, TauAmapIdx{i}] = bin_quickly2d(azmin, azmax, elmin, elmax, pixel_size,...
0564             Azoff(~flagvector3), Eloff(~flagvector3),allIQUP_pa{1,i}(~flagvector3,3));
0565         figure(m)
0566         subplot(2,2,i)
0567         
0568         %imagesc(fliplr(xidx),yidx,TauAmap{i},clim(i,:))
0569         imagesc(fliplr(xidx),yidx,TauAmap{i})
0570         title([source_name,' :',map_labels(i)])
0571         xlabel('AZ offset, deg')
0572         ylabel('EL offset, deg')
0573         colormap('default')
0574         colorbar
0575         axis square
0576     end