south_beam_maps_quick

 21/2/2014 - ACT
 raster and large az and el tracks on the sun

0001 
0002 function [d,TauAmap, azOffSave, elOffSave,xidx,yidx,LL, RR, L1, L2, QQ, UU, mask,sf] = south_beam_maps_quick(start_date,stop_date,comment)
0003 % 21/2/2014 - ACT
0004 % raster and large az and el tracks on the sun
0005 
0006 %d = read_arcSouth('21-Feb-2014:13:23:03 ','     21-Feb-2014:13:57:26 ')
0007  %d = read_arcSouth('16-Jan-2015:17:16:21 ','     16-Jan-2015:17:31:51')
0008 %d = read_arcSouth('23-Feb-2014:13:20:20','23-Feb-2014:13:54:44')
0009 
0010 %d = read_arcSouth('16-Jan-2015:19:21:12','16-Jan-2015:20:26:09')%TauA
0011 %d = read_arcSouth('16-Jan-2015:20:26:09','16-Jan-2015:20:26:10') % M42
0012 %d = read_arcSouth('17-Jan-2015:05:41:43','17-Jan-2015:06:41:50')%Moon
0013 %d = read_arcSouth('17-Jan-2015:17:45:40','      17-Jan-2015:18:19:42')%M42
0014 %d = read_arcSouth('17-Jan-2015:10:17:02','17-Jan-2015:10:32:32')
0015 d = read_arcSouth(start_date,stop_date)
0016 
0017 
0018 %%%%
0019 %
0020 % HACK - 20/1/2015 while no 1pps
0021 %
0022 %antenna0.servo.ntpUSecond[0][0]=0
0023 %antenna0.servo.ntpUSecond[0][1]=200000
0024 %antenna0.servo.ntpUSecond[0][2]=400000
0025 %antenna0.servo.ntpUSecond[0][3]=600000
0026 %antenna0.servo.ntpUSecond[0][4]=800000
0027 %d.antenna0.servo.ntpUSecond=repmat([0 200000 400000 600000 800000]',int64(size(d.antenna0.servo.ntpUSecond)/5),1);
0028 
0029 
0030 % % Horn at 45 deg, 0dBm-30dB attenuation
0031 % close all;
0032 % clear;
0033 % %d = read_arcSouth('23-Feb-2014:15:09','23-Feb-2014:15:30')
0034 % d = read_arcSouth('23-Feb-2014:15:07:51','23-Feb-2014:15:24:02')
0035 %
0036 % %Horn at 0deg, 0dbm-30dB attenuation
0037 % close all;clear;
0038 % d = read_arcSouth('23-Feb-2014:15:27','23-Feb-2014:15:40')
0039 %
0040 % %Horn at 0deg, 0dbm-60dB attenuation
0041 % close all;clear;
0042 % d = read_arcSouth('23-Feb-2014:16:05','23-Feb-2014:16:15')
0043 
0044 
0045 %%
0046 % Put data into Cbass-N style - first grab just the raster data
0047 
0048 %d = framecut(d,d.index.beammap.fast);
0049 
0050 [d LL RR L1 L2 QQ UU mask] = cbassS2N(d,1);
0051 
0052 nchan=4
0053 plot_labels = {'LL','RR','QQ','UU'};
0054 chan_order=[1 6 2 3];
0055 %%
0056 % Detrend the data ready for analysis
0057 
0058    
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         
0088         %Detrend the data using data from 2 degrees away from source
0089         angle = sqrt(elOffSave.^2 + azOffSave.^2);
0090         
0091         
0092         for chan = chan_order %
0093             % Continue detrending channel by channel
0094                 temp_data = dcut.antenna0.receiver.data(:,chan);
0095                 temp_time = dcut.antenna0.receiver.utc;
0096                 %if(strmatch('S',cbass,'exact'))
0097                     temp_y = temp_data(angle>2); % Set to 7 for south data __> CHANGE once pointing fixed
0098                     temp_x = temp_time(angle>2);
0099                 %else
0100                 %    temp_y = temp_data(angle>2); % Set to 2 for north  data
0101                 %    temp_x = temp_time(angle>2);
0102                 %end
0103                 
0104 % In case there are nans in south data from way we've done the flag kludge
0105                                                                  
0106                %if(strmatch('S',cbass,'exact'))
0107                   ignore = isnan(temp_y);
0108                   p = polyfit(temp_x(~ignore),temp_y(~ignore),1);
0109                %else
0110                %   p = polyfit(temp_x,temp_y,1);
0111                %end
0112 
0113                 all_time = temp_time;
0114                 f = polyval(p,all_time);
0115                 
0116                 % Now subtract f from data channel on that scan
0117                 dcut.antenna0.receiver.data(:,chan) = dcut.antenna0.receiver.data(:,chan) - f;
0118                 %Keep a master copy too of all the data
0119                 dgood.antenna0.receiver.data(ind,chan)=dcut.antenna0.receiver.data(:,chan);
0120                 
0121                 
0122                 clear temp_data temp_time temp_x temp_y p f all_time
0123         end
0124     end
0125     %%
0126     % Plot the detrended raster scan data
0127     close all
0128     figure(1)
0129     nchan=4 % d.data has [I1,Q1,U1,Q2,U2,Q3,U3,I2,V]
0130     plot_labels={'LL','RR','QQ','UU'}
0131     
0132     points = (dgood.index.beammap.fast==1);% so we can identify the raster data
0133     
0134     for i=1:nchan
0135     if i<7
0136         figure(1)
0137         subplot(2,3,i)
0138         
0139     else
0140         figure(2)
0141         nplot=i-6;
0142         subplot(2,2,nplot)
0143         
0144     end
0145         
0146          
0147     
0148     plot(dgood.antenna0.receiver.utc(points),dgood.antenna0.receiver.data(points,chan_order(i)))
0149     xlabel(plot_labels{i})
0150     gtitle('Pipeline calibrated and detrended data')
0151     end
0152 
0153 %% Get offsets in Az and El for rasters
0154 d= framecut(dgood,dgood.index.beammap.fast);
0155     azApp = interp1(d.antenna0.tracker.utc, ...
0156         d.antenna0.tracker.horiz_topo(:,1), d.antenna0.roach2.utc);
0157     azOffSave = (azApp - d.antenna0.servo.apparent(:,1));
0158     azOffSave = -wrap180(azOffSave);
0159     azOffSave = azOffSave.*cosd(d.antenna0.servo.el);% Scale azimuth by cos(elevation)???CHECK
0160     
0161     elApp = interp1(d.antenna0.tracker.utc, ...
0162         d.antenna0.tracker.horiz_topo(:,2), d.antenna0.roach2.utc);
0163     elOffSave = elApp - d.antenna0.servo.apparent(:,2);
0164     elOffSave = -elOffSave;
0165     d.angle = sqrt(elOffSave.^2 + azOffSave.^2);%
0166     d.elOffSave = elOffSave;
0167     d.azOffSave= azOffSave;
0168     
0169 
0170 
0171 % %%
0172 %
0173 % % Now analyse the Raster data
0174 %
0175 % pixel_size = 0.1; % degree
0176 %
0177 % azmin = -2;
0178 % azmax = +2;
0179 % elmin = -2;
0180 % elmax = +2;
0181 % nanflag = isnan(d.antenna0.receiver.data(:,1));
0182 % flag = ~nanflag;
0183 % colLR = 1e5;
0184 % colQU = 1e2;
0185 % % Set colour scale mapping
0186 % %clim = [0 3.5;-0.2 0.2;-0.2 0.2;0 0.2]
0187 %
0188 % %for m= 1:length(s)
0189 %  %   n = n+1;
0190 %  %   ind = zeros(size(FLAG));
0191 %  %   ind(s(m):e(m)) = 1;
0192 %  %   ind = logical(ind);
0193 %  %   flagvector3 = (~nanflag & ~FLAG & ind);
0194 %
0195 %     for i=1:4
0196 %         j=chan_order(i);
0197 %         [TauAmap{i}, xidx, yidx, TauAmapIdx{i}] = bin_quickly2d(azmin, azmax, elmin, elmax, pixel_size,...
0198 %             d.azOffSave(flag), d.elOffSave(flag),d.antenna0.receiver.data(flag,j));
0199 %         figure(3)
0200 %         subplot(2,2,i)
0201 %         if i<3
0202 %         clim=[0.5e4 colLR];
0203 %     else
0204 %         clim=[-colQU colQU];
0205 %     end
0206 %         %imagesc(fliplr(xidx),yidx,TauAmap{i},clim(i,:))
0207 %         %clim=[5e3 1.5e4]
0208 %        imagesc(fliplr(xidx),fliplr(yidx),TauAmap{i})
0209 %         %imagesc(fliplr(xidx),fliplr(yidx),TauAmap{i},clim)
0210 %
0211 %         %title([source_name,' map:', num2str(i)])
0212 %         colormap('default')
0213 %         colorbar
0214 %         title(plot_labels{i})
0215 %         axis square
0216 %         figure(4)
0217 %         subplot(2,2,i)
0218 %         clines = [1:-0.01:0]*max(max(TauAmap{i}));
0219 %         contour(fliplr(xidx),yidx,TauAmap{i})
0220 %         %contour(fliplr(xidx),yidx,TauAmap{i},clines)
0221 %         title(plot_labels{i})
0222 %         axis square
0223 %
0224 %     end
0225 %
0226 %
0227 %
0228 %% Plot the raster
0229 
0230 
0231 
0232 disp([ 'Plotting the raster: '] )
0233 azOffRaster=d.azOffSave;
0234 elOffRaster=d.elOffSave;
0235 pixel_size = 0.15; % degree
0236 
0237 azmin=-5.0
0238 azmax=5.0
0239 elmin=-5.0
0240 elmax=5.0
0241 %azmin = min(azOffRaster);
0242 %azmax = max(azOffRaster);
0243 %elmin = min(elOffRaster);
0244 %elmax = max(elOffRaster);
0245 colmin = 0
0246 colmax = 1e5
0247     
0248 for i=1:4
0249     [TauAmap{i}, xidx, yidx, TauAmapIdx{i}] = bin_quickly2d(azmin, azmax, elmin, elmax, pixel_size,...
0250         azOffRaster, elOffRaster,d.antenna0.receiver.data(:,chan_order(i)));
0251     
0252     % Raster Maps
0253     figure(1)
0254     subplot(2,2,i)
0255     %clim=[colmin colmax]
0256     imagesc(fliplr(xidx),fliplr(yidx),TauAmap{i})
0257     colormap('default')
0258     title(plot_labels{i})
0259     axis equal
0260     ylim( [elmin elmax] )
0261     xlim( [azmin azmax] )
0262     colorbar
0263     
0264     %Contour Plots
0265 %     figure(4)
0266 %     subplot(2,2,i)
0267 %     contour(fliplr(xidx),yidx,TauAmap{i})
0268 %     title(plot_labels{i})
0269 %     axis equal
0270 %     ylim( [elmin elmax] )
0271 %     xlim( [azmin azmax] )
0272 end
0273 
0274 % %% Look at cuts through the raster maps
0275 disp( 'Look at cuts through the beam' )
0276 % Get the user to select the centre of the beam pattern
0277 figure(5)
0278 climI = [0 1e5]
0279 colormap('default')
0280 imagesc(TauAmap{1})
0281 title('Please click on the centre of the LL beam pattern')
0282 [usr_TrigAz1 usr_TrigEl1] = ginput(1);
0283 figure(5)
0284 imagesc(TauAmap{2})
0285 title('Please click on the centre of the RR beam pattern')
0286 [usr_TrigAz2 usr_TrigEl2] = ginput(1);
0287 figure(5)
0288 imagesc(TauAmap{3})
0289 title('Please click on the centre of the QQ beam pattern')
0290 [usr_TrigAz3 usr_TrigEl3] = ginput(1);
0291 figure(5)
0292 imagesc(TauAmap{4})
0293 title('Please click on the centre of the UU beam pattern')
0294 [usr_TrigAz4 usr_TrigEl4] = ginput(1);
0295 
0296 usr_TrigAz = mean( [usr_TrigAz1 usr_TrigAz2 usr_TrigAz3 usr_TrigAz4] );
0297 usr_TrigEl = mean( [usr_TrigEl1 usr_TrigEl2 usr_TrigEl3 usr_TrigEl4] );
0298 
0299 usr_TrigAz = round( usr_TrigAz );
0300 usr_TrigEl = round( usr_TrigEl );
0301 
0302 % Plot the two cuts through the beam
0303 figure(5)
0304 subplot(2,4,1)
0305 plot((TauAmap{1}(usr_TrigEl,:)))
0306 ylabel( 'LL' )
0307 xlabel( 'Azimuth' )
0308 subplot(2,4,5)
0309 plot((TauAmap{1}(:,usr_TrigAz)))
0310 ylabel( 'LL' )
0311 xlabel( 'Elevation' )
0312 
0313 
0314 subplot(2,4,2)
0315 plot(TauAmap{2}(usr_TrigEl,:))
0316 ylabel( 'RR' )
0317 xlabel( 'Azimuth' )
0318 subplot(2,4,6)
0319 plot(TauAmap{2}(:,usr_TrigAz))
0320 ylabel( 'RR' )
0321 xlabel( 'Elevation' )
0322 
0323 subplot(2,4,3)
0324 plot(TauAmap{3}(usr_TrigEl,:))
0325 ylabel( 'QQ' )
0326 xlabel( 'Azimuth' )
0327 subplot(2,4,7)
0328 plot(TauAmap{3}(:,usr_TrigAz))
0329 ylabel( 'QQ' )
0330 xlabel( 'Elevation' )
0331 
0332 subplot(2,4,4)
0333 plot(TauAmap{4}(usr_TrigEl,:))
0334 ylabel( 'UU' )
0335 xlabel( 'Azimuth' )
0336 subplot(2,4,8)
0337 plot(TauAmap{4}(:,usr_TrigAz))
0338 ylabel( 'UU' )
0339 xlabel( 'Elevation' )
0340 %%
0341 % One big figure of I to check sidelobes of Moon
0342 figure(30)
0343 imagesc(fliplr(xidx),fliplr(yidx),TauAmap{1},[0 300])
0344 
0345 xlabel('Az offset, deg');
0346 ylabel('El offset, deg');
0347 title(['Raster of ',num2str(start_date)]);
0348 axis equal
0349 % Save this map
0350 map_name=['MirrorMap_',start_date,'_',comment,'.png'];
0351 print('-f30','-dpng',map_name)
0352 
0353 [X,Y]=meshgrid(xidx,yidx);
0354 [sf gof t] = Egauss2D_fit_act(X,Y,TauAmap{1,1});
0355 
0356 fwhm_x = 2*sqrt(2*log(2)).*sf.sigma_x
0357 fwhm_y = 2*sqrt(2*log(2)).*sf.sigma_y
0358 amplitude=sf.a1
0359 x0=sf.x0
0360 y0=sf.y0
0361 theta=sf.theta
0362 
0363 %% Write out the data to a file
0364 filename = 'Dish_Focussing.dat';
0365 fid = fopen(filename, 'a+');
0366 
0367 data2save={start_date, comment, x0,y0,theta,fwhm_x,fwhm_y};
0368 fprintf( fid,'%s %s %d %d %d %d %d\n', data2save{1,:})
0369 
0370 fclose(fid);
0371 end
0372 
0373