% make var weighted "dirty maps" of all 8 fields

load dvia
ad=calc_ad(5,256);

figure(1)
for i=1:length(dvi)
  % make the "dirty map" with variance weighting
  map(i).dmap=make_dmap(ad,dvi(i).u,dvi(i).v,dvi(i).vis,dvi(i).var);
  
  % estimate noise in map as map std
  % - actually this is very crude as point sources are strong enough
  % to blow up map std - should do a fit to central region of
  % histogram of pixel values
  map(i).std=std(map(i).dmap(:));
  
  % make "significance map" in these units
  map(i).sigim=map(i).dmap/map(i).std;
  
  % plot the maps
  subplot(2,4,i);
  imagesc(ad.t_val_deg,ad.t_val_deg,map(i).sigim);
  axis xy; axis image
  set(gca,'XDir','reverse')
  title(sprintf('field A%d',i));
  colorbar
end

% DASI obs paper astro-ph/0104488 tab 2 A row RA22h Dec -61
rac=[22,23,0,1,2,3,4,5]*15; decc=-61;

% from ftp://ftp.atnf.csiro.au/pub/data/pmn/PLAINTXT get pmns.txt
[p.name,rh,rm,rs,dd,dm,ds,p.s,p.se]=textread('pmns.txt','%s %f:%f:%f %f:%f:%f %f %f  %*[^\n]');

% remember that -ve sign applies to min/sec part of dec
p.ra =15*(rh+rm/60+rs/60^2);
p.dec=(abs(dd)+dm/60+ds/60^2).*sign(dd);
% or can use matlab provided dms2deg function

% plot the pmn data
figure(2); plot(p.ra,p.dec,'.','MarkerSize',1);
set(gca,'XDir','reverse'); axis tight

% to catch the sources in the field which stradles zero RA we double
% up a bit
ind=p.ra>355;
p.name=[p.name;p.name(ind)];
p.ra=[p.ra;p.ra(ind)-360];
p.dec=[p.dec;p.dec(ind)];
p.s=[p.s;p.s(ind)];
p.se=[p.se;p.se(ind)];

% plot the maps again with sources overplotted
figure(3); clf
for i=1:length(map)
  subplot(2,4,i);
  y_tic=ad.t_val_deg+decc;
  x_tic=ad.t_val_deg/cos(decc*pi/180)+rac(i);
  
  imagesc(x_tic,y_tic,map(i).sigim);
  axis xy; axis square; set(gca,'XDir','reverse')
  title(sprintf('field A%d, max=%.1f',i,max(map(i).sigim(:))));
    
  % overplot the PMN sources and their names
  rai=x_tic(1)500;
  hold on
  scatter(p.ra(pind),p.dec(pind),p.s(pind),'k+');
  text(p.ra(pind),p.dec(pind),p.name(pind),'HorizontalAlignment','center');text(p.ra(pind),p.dec(pind),p.name(pind),'HorizontalAlignment','center');
  hold off
end

% using zoom read source names off map
% sigma is just max of map in all cases but one where we can use
% impixelinfo to read off the value (although appears have to replot
% map6 alone to make it work)
%
% source list
% field  sigma  pmn name
% A3     14.3   PMNJ2358-6054
% A4      8.0   PMNJ0109-6049
% A6      8.3   PMNJ0303-6211
% A6      7.4   PMNJ0309-6058
% A8     14.5   PMNJ0506-6109


function map=make_dmap(ad,u,v,vis,var)
% map=make_dmap(ad,u,v,vis,var)
%
% Make simple "dirty map" by gridding visibilities and fft

if(exist('var')~=1)
  var=complex(1,1)*ones(size(vis));
end

% add conj vis
u=[u;-u]; v=[v,-v];
vis=[vis;conj(vis)];
var=[var;var];

% find low/high f grid edges
n=ad.N_pix;
l=ad.u_val(1)-ad.del_u/2;
h=ad.u_val(end)+ad.del_u/2;

% accumulate mean real and imag in each cell
[x_tic,y_tic,r]=hfill2(u,v,n,l,h,n,l,h,real(vis)./real(var));
[x_tic,y_tic,i]=hfill2(u,v,n,l,h,n,l,h,imag(vis)./imag(var));
z=complex(r,i);

% Normally we would normalize here by the sum of the weights (1/var)
% but since we are not trying to preserve overall normalization we
% don't bother.
% Note that there is no need to accumulate nhits and divide that out -
% each vis counts the same regardless of whether it happens to share a
% fourier plane cell or not.

% force cells where we have no obs to zero
z(isnan(z))=0;

% do the fft
% must use fft2 rather than ifft2 to get point sources in correct
% locations
map=real(fftshift(fft2(fftshift(z))));

return

QUaD is a simple "scanning beam" experiment - each detector accepts radiation from a small patch on the sky and the whole telescope moves around to scan the beams across the sky. Let's load in data which has been subjected to low level processing only - it is basically raw and uncalibrated.

clear
load('060607_tod_uncal');

Remember we can see what variables are available by typing "whos" and we can what's in a strcuture by typing the name of the the structure.

>> d
      frame: [1x1 struct]
     lockin: [1x1 struct]
       bias: [1x1 struct]
     fridge: [1x1 struct]
    tracker: [1x1 struct]
       pmac: [1x1 struct]
    weather: [1x1 struct]
        cal: [1x1 struct]
     thermo: [1x1 struct]
          t: [63209x1 double]
         tf: [1264180x1 double]
         td: [63209x1 double]
         az: [1264180x1 double]
         el: [1264180x1 double]
         ra: [1264180x1 double]
        dec: [1264180x1 double]

The data is arranged in sub-structures by data acquisition "board". Some of these are real pieces of hardware like the weather station and some are virtual - information pulled together from a variety of sources. The arrays at base level have been generated for convenience by the low level processing. Note that there is t ("time") and tf ("time fast") - most quantities get read out once per second, but the most important data like the bolometer readouts themselves are read out 20 times per second. Hence the tf vector is 20 times longer than the t vector.

Let's plot the air temp over the run:

plot(d.t,d.weather.air_temp,'.');

The time t is measured in seconds into the UTC day - the run starts around 5pm (60,000 seconds), wraps around at midnight (86,400 seconds) and ends around 10:30am (37,000 seconds). Yes it really was below -60 degrees centigrade when this data was taken!

Now let's plot the telescope azimuth angle - the rotation angle around the vertical axis:

plot(d.tf,d.az,'.','Markersize',0.2);

Note we need to use tf as az is read at fast rate. Use zoom button (on top bar of plot window) to zoom in on the main part of the run - we see the telescope tracking the source and executing triangle wave azimuth scans as a modulation on top the track. Note that at the Pole a source on the sky just goes round and round in azimuth angle remaining at constant elevation.

plot(d.tf,d.el,'.','Markersize',0.2);

Zooming in on this we can see elevation angle slowly increasing with periodic blips up/down.

There is a lot of other stuff going on at the start and end of the run which is observations of calibration sources etc. We want to strip this away and look only at the "on field scanning" portion of the run. The low level processing has gone through and identified the start/end of each scan and provided this info in structure "fs" (for "field scans").

>> fs
      s: [808x1 double]
      e: [808x1 double]
      d: [808x1 double]
     sf: [808x1 double]
     ef: [808x1 double]
      t: [808x1 double]
    std: [808x70 double]

The sf and ef are the "start fast" and "end fast" of each half scan - the out and back motions in az modulated on top of the track. So we can plot the first half scan like so:

i=1; s=fs.sf(i); e=fs.ef(i);
plot(d.tf(s:e),d.az(s:e))

plot(d.tf(mapind),d.az(mapind),'Markersize',0.2)

And see just the bits where we are scanning on source. To do this we use a for loop but it only loops 808 times so it's OK:

function mapind=make_mapind(d,fs)
% mapind=make_mapind(d,fs)
%
% Build index array pointing from scan start/stop pointers

mapind=false(size(d.lockin.adcData,1),1);
for i=1:length(fs.s)
  s=fs.sf(i); e=fs.ef(i);
  mapind(s:e)=true;
end

return

Now we can plot the az angle as above, or better we can plot dec vs ra and see the raster scan pattern which the telescope executed over the target field:

mapind=make_mapind(d,fs);
plot(d.ra(mapind),d.dec(mapind),'.','MarkerSize',0.2); axis tight

Again zoom in to see detailed structure.

plot(d.tf(mapind),d.lockin.adcData(mapind,1),'.','MarkerSize',0.1);

The first column of adcData is the first channel. We see messy structure around -0.4 on the y-axis. This is because the DC offset of the bolometer channels is arbitrary and drifts over time - we scan the telescope around on the sky to make differential measurements of power level at nearby positions on the sky. The DC offset gets reset by the readout electronics between groups of scans which is why we see large jumps on the plot above. We will proceed to make a map for channel 1 but the DC offsets are going to make the map look all stripey:

To keep things clean we make a function map=make_map(x,y,dat):

function map=make_map(x,y,dat)
% map=make_map(x,y,dat)
%
% Make a 2D map of some data

[map.x_tic,map.y_tic,map.nhit]=hfill2(x,y,100,[],[],100);
[map.x_tic,map.y_tic,map.tot]=hfill2(x,y,100,[],[],100,[],[],dat);
map.map=map.tot./map.nhit;

return

NB: you need to recopy hfill2.m from ~pryke/matlab - the old version assumed its input args were double precision and the (very large) lockin.adcData is stored as single to save space.

Now we can do:

map=make_map(d.ra(mapind),d.dec(mapind),d.lockin.adcData(mapind,1));
imagesc(map.x_tic,map.y_tic,map.map); axis xy; set(gca,'XDir','reverse')

The map looks really horrible due to the un-subtracted DC offsets.

• Make a function d=clean_scans(d,fs,percent) which will loop through each half scan and subtract the given percenticle of the intregral distribution of the points in that scan and channel. Use the function prctile to find the percentile value (prctile is a generalization of median but can find percentile values other than 50%). Note that your function MUST NOT do a second loop over the bolometer channels - instead it must make use of the fact that prctile can work over all columns in a single shot like mean, median, max, min etc. Use the 25th percentile to start with. Now plot the maps - the row to row striping is gone although there is still some bogus looking large scale structure across the map. We now see the target source in the middle of the map.
• From the RA/Dec of the blip in the middle of the map find out its identity.
• The small periodic up/down moves in elevation which we noted earlier are calibration "elnods". Before each group of on field scans we inject a ramp into the data stream using the gradient in atmospheric emission. In class 2 we looked at some of this data and regressed it versus air mass to extract relative calibration coefficients for the bolometers channels. load in the file "060607_calval" to get access to these values as extracted by the low level data reduction - they are stored in "en.g" (for elnod.gain). Write a function d=cal_scans(d,fs,en) which will relative cal each of the channels to the same gain. Remember that each scan should be cal'd according to the relative gains of the preceding elnod. The time at which each scan starts is given by fs.t and the time of each elnod by en.t
• The columns of d.lockin.adcData go A,B,A,B... where A and B are the bolometers in each polarization pair. So by taking A+B and A-B for each pair we are generating the total intensity and polarization data streams. Write a function d=sumdiff_pairs(d,fs) which will convert the columns of d.lockin.adcData from A,B,A,B... to sum,diff,sum,diff... Now make the sum and diff maps for the first bolometer pair - you should see the source in the sum map as before, but a much weaker signature (if any) in the diff map (because the source has fractional polarization much less than one).