Here are the clean_scans, cal_scans and sumdiff_pairs functions:

function d=clean_scans(d,fs,percent)
% d=clean_scans(d,fs,percent)
%
% Remove the given percentile of the integral distribution from each
% half scan specified by fs, and for each bolometer channel in
% d.lockin.adcData

for i=1:length(fs.sf)
  s=fs.sf(i); e=fs.ef(i);
  
  %v=double(d.lockin.adcData(s:e,ind));
  v=d.lockin.adcData(s:e,:);
  
  v=v-repmat(prctile(v,percent),[size(v,1),1]);
    
  d.lockin.adcData(s:e,:)=v;
end

return

function d=cal_scans(d,fs,en)
% d=cal_scans(d,fs,en)
%
% relative cal fs scans using en gains so that all channels have same
% gain

% for each field scan
for j=1:length(fs.t)
  
  % find the preceding el nod
  x=fs.t(j)-en.t; x(x<0)=+inf; [dummy,k]=min(x);
  
  % Scale with rel gain
  s=fs.sf(j); e=fs.ef(j);
  cal=en.g(k,1:size(d.lockin.adcData,2));
  cal=repmat(cal,[(e-s+1),1]);
  d.lockin.adcData(s:e,:)=d.lockin.adcData(s:e,:)./cal;
end

function d=sumdiff_pairs(d,fs)
% d=sumdiff_pairs(d,fs)
%
% Take A/B bolometer data to sum/diff

% for each half scan
for i=1:length(fs.sf)
  s=fs.sf(i); e=fs.ef(i);
  v=d.lockin.adcData(s:e,:);
  
  vp=v;
  vp(:,1:2:end)=v(:,1:2:end)+v(:,2:2:end);
  vp(:,2:2:end)=v(:,1:2:end)-v(:,2:2:end);
  
  d.lockin.adcData(s:e,:)=vp;
end

return

Here is a simple script which uses them to make maps:

% get the "time ordered data"
load 060607_tod_uncal

% get the "map index" which points to the on field scanning portion of
% the run
mapind=make_mapind(d,fs);

% clean the scans to remove DC offset from each
d=clean_scans(d,fs,25);

% make maps for A/B channels of first pair
map(1)=make_map(d.ra(mapind),d.dec(mapind),d.lockin.adcData(mapind,1));
map(2)=make_map(d.ra(mapind),d.dec(mapind),d.lockin.adcData(mapind,2));
subplot(1,2,1)
imagesc(map(1).x_tic,map(1).y_tic,map(1).map);
axis xy; axis square; set(gca,'XDir','reverse'); colorbar
subplot(1,2,2)
imagesc(map(2).x_tic,map(2).y_tic,map(2).map);
axis xy; axis square; set(gca,'XDir','reverse'); colorbar

% apply elnod based relative gain calibration to equalize the gain of
% all channels
load 060607_calval
d=cal_scans(d,fs,en);

% Go from A/B/A/B to sum/diff/sum/diff
d=sumdiff_pairs(d,fs);

% remap to see result
map(1)=make_map(d.ra(mapind),d.dec(mapind),d.lockin.adcData(mapind,1));
map(2)=make_map(d.ra(mapind),d.dec(mapind),d.lockin.adcData(mapind,2));
subplot(1,2,1)
imagesc(map(1).x_tic,map(1).y_tic,map(1).map);
axis xy; axis square; set(gca,'XDir','reverse'); colorbar
title('pair sum')
subplot(1,2,2)
imagesc(map(2).x_tic,map(2).y_tic,map(2).map);
axis xy; axis square; set(gca,'XDir','reverse'); colorbar
title('pair diff')

The end product of this looks like this:

In the left panel we see a map of polarized emission from the CenA radio galaxy - the central core and 2 jets extending to upper left and lower right. The pair diff maps shows one compoent of polarization - note that the blob is not coincident with the core - we are seeing high fractional polarization of the upper left jet. Note that we see extensive atmospheric stipping in the left panel, but closer to white noise in the right as atmospheric emission is laregly un-polarized.

% make maps for all pairs
for i=1:size(d.lockin.adcData,2)
  map(i)=make_map(d.ra(mapind),d.dec(mapind),d.lockin.adcData(mapind,i));
end

% plot maps for all pairs
for i=1:length(map)/2
  subplot(6,6,i)
  m=(i-1)*2+1;
  %m=i*2;
  imagesc(map(m).x_tic,map(m).y_tic,map(m).map);
  axis xy; axis square; set(gca,'XDir','reverse');
end

The plot looks like this - we see the source at a different location in each panel - why? The feeds in the focal plane each see a different spot on the sky for a given telescope pointing direction. Therefore as the telescope raster scans around each feed passes over the source at different pointing offset from the nominal pointing direction. The RA/Dec on the plot axes are only correct for the center feed - which by definition has zero offset angle.

To get info on the offset angles of each feed we load another file containing structures p and ind.

load 060607_arrayinfo.mat
plot(p.ra_off_dos,p.dec_off,'+'); grid
axis image

The QUaD focal plane has 7 inner feeds at 150GHz, a ring of 12 100GHz and 12 more 150's around the edge. We make a function fp_diag which will display the the position, beam size, and grid angles of the A and B grids in each feed.

function fp_diag(p)
% fp_diag(p)
%
% Make pretty picture showing focal plane layout

x=p.ra_off_dos;
y=p.dec_off;

% draw circles with beam width
[xc,yc]=circle(x,y,p.beam_sig_a);
plot(xc,yc,'b'); axis image; grid
set(gca,'XDir','reverse');

% draw lines at grid angle
t=(p.grid_phi+90)*pi/180;
r=p.beam_sig_a;
[xo,yo]=pol2cart(t,r);

xl=[x+xo;x-xo];
yl=[y+yo;y-yo];

line(xl(:,1:2:end),yl(:,1:2:end),'Color','b');
line(xl(:,2:2:end),yl(:,2:2:end),'Color','r');

return

Note that the grid angle is shown as projected on the sky and there are grid pairs of two flavours. The angle shown is the metalized grid direction - photons with their E-field aligned with the metalized direction of the bolometer grid will excite voltage in the grid in the same direction and therefore be absorbed - they see the grid as a bunch of parallel dipole antennas. Photons with their E-field perpendicular to the grid don't see it at all and pass straight through. Photon's at intermediate angles have a probability to be absorbed or pass through depending on their angle to the metalized direction. Note that this is exactly the same way polarized sunglasses work - the explaination sometimes given of photons parallel to the metal bars passing through as in through the bars of a cage is exactly wrong.