mosaic_map_fast

 mos=mosaic_map_fast(ad,ads,p,beam,uvcut)

 Make mosaic map by linear addition of individual dirty maps
 Called linear mosaicing. e.g. eqn 9&10 in astro-ph/0205388
 or eqn 20-2 in SIIRA.

 Alternate version doing dmap shift by
 manipulating visibility phase.

 ad is big grid
 ads is per point grid
 beam is stack of primary beams one per band
 uvcut is uvr cut range

0001 function [mos, pmap]=mosaic_map_fast(ad,ads,p,beam,uvcut, noiseType)
0002 % mos=mosaic_map_fast(ad,ads,p,beam,uvcut)
0003 %
0004 % Make mosaic map by linear addition of individual dirty maps
0005 % Called linear mosaicing. e.g. eqn 9&10 in astro-ph/0205388
0006 % or eqn 20-2 in SIIRA.
0007 %
0008 % Alternate version doing dmap shift by
0009 % manipulating visibility phase.
0010 %
0011 % ad is big grid
0012 % ads is per point grid
0013 % beam is stack of primary beams one per band
0014 % uvcut is uvr cut range
0015 
0016 if(~exist('uvcut'))
0017   uvcut=[0,Inf];
0018 end
0019 
0020 if(~exist('noiseType'))
0021   noiseType = 1;
0022 end
0023 
0024 % Take the mean beam to unfold with - in principle one could
0025 % make a dmap for each band separately and unfold band by band...
0026 if(ndims(beam)==3)
0027     beam=mean(beam,3);
0028 end
0029 
0030 % Generate individual maps
0031 a=zeros(ad.N_pix); b=a;
0032 a1=zeros(ad.N_pix); b1=a1;
0033 
0034 for i=1:length(p.px)
0035     i
0036     
0037   % Filter vis
0038   uvr=pyth(p.u{i},p.v{i});
0039   ind=uvr>uvcut(1)&uvr<uvcut(2);
0040   
0041   tic
0042   % Shift phase center
0043   % Cross product of phase center position and baseline vectors
0044   x=2*pi*(p.u{i}(ind)*p.px(i)*pi/180+p.v{i}(ind)*p.py(i)*pi/180);
0045   svis=p.vis{i}(ind).*complex(cos(x),sin(x));
0046   
0047   % Make dirty map
0048   [dmap,s]=difmap(ad,p.u{i}(ind),p.v{i}(ind),svis,p.var{i}(ind));
0049   toc
0050   s1 = s;
0051   % Shift beam only to pointing center - dmap already shifted
0052   % find where to center the beam
0053   [xrm, yrm, xrb, yrb] = find_beam_range(ad, ads, p.px(i), p.py(i));
0054   sbeam = zeros(ad.N_pix);
0055   sbeam(yrm, xrm) = beam(yrb, xrb);
0056   %  beam done is this way can be off by ~.1%...which is nothing.
0057   %  if this is a problem for you, you can uncomment the following
0058   %  lines and be bored while you wait.
0059   %x=ads.x+p.px(i); y=ads.y+p.py(i);
0060   %tic
0061   %sbeam2=interp2(x,y,beam,ad.x,ad.y); sbeam(isnan(sbeam))=0;
0062   %toc
0063 
0064   % let's save some values for future reference
0065   pmap.beam.xrm{i} = xrm;
0066   pmap.beam.yrm{i} = yrm;
0067   pmap.beam.xrb{i} = xrb;
0068   pmap.beam.yrb{i} = yrb;
0069   pmap.s1(i) = s1;
0070   pmap.map{i} = single(dmap(yrm, xrm));
0071 
0072   if(noiseType==1)
0073     s = newNoise(dmap);
0074   end
0075   pmap.s(i) = s;
0076  
0077   s = s1;
0078   % Acculmulate arrays to make mosaic map
0079   a=a+sbeam.^2./s;
0080   b=b+sbeam.*dmap./s;
0081   
0082   %a1=a1+sbeam.^2./s1;
0083   %b1=b1+sbeam.*dmap./s1;
0084 end
0085 
0086 % Calculate noise and signal maps
0087 %mos.noise1=single(sqrt(1./a1));
0088 %mos.signal1=single(b1./a1);
0089 mos.noise=single(sqrt(1./a));
0090 mos.signal=single(b./a);
0091 
0092 % calc "significance image"
0093 %mos.sigma=single(mos.signal1./mos.noise1);
0094 mos.sigma=single(mos.signal./mos.noise);
0095 
0096 return
0097 
0098  
0099 function   [xrm, yrm, xrb, yrb] = find_beam_range(ad, ads, px, py)
0100 
0101   if(isfield(ad, 't_val_x'))
0102       x = ad.t_val_x_deg - px;
0103   else
0104       x = ad.t_val_deg - px;
0105   end
0106   y = ad.t_val_deg - py;
0107   x = find(abs(x)==min(abs(x)));
0108   y = find(abs(y)==min(abs(y)));
0109   center = floor(ads.N_pix/2)+1;
0110   xrb = 1:ads.N_pix;
0111   yrb = 1:ads.N_pix;
0112   xrm = xrb - center + x;
0113   yrm = yrb - center + y;
0114   % now we need to crop;  we crop if:  beam gets cut off due to
0115   % edge of map, map gets cut off due to size of beam.
0116   f = find(xrm <=0);
0117   xrm(f) = [];
0118   xrb(f) = [];
0119   f = find(xrm > ad.N_pix);
0120   xrm(f) = [];
0121   xrb(f) = [];
0122   f = find(yrm <=0);
0123   yrm(f) = [];
0124   yrb(f) = [];
0125   f = find(yrm > ad.N_pix);
0126   yrm(f) = [];
0127   yrb(f) = [];  
0128   
0129   
0130 return;
0131 
0132 function noise = newNoise(map)
0133 
0134 %  this is where we try to calculate the noise from the maps
0135 %  themselves
0136 %   i wrote it once and it was lost.  let's try again.
0137 %   i think all we need to do is find the really bright pixels, and
0138 %   remove them, and then calculate our noise. although it won't
0139 %   matter too much because the map is so large anyway.
0140 %  also...i've demonstrated (not here, of course) that once we
0141 %  remove the really bright sources, that the noise goes to
0142 %  something sensible.
0143 
0144 lowSig      = min(map(:));
0145 hiVals      = find(map(:)>abs(lowSig));
0146 map(hiVals) = nan;
0147 noise = nanvar(map(:));
0148 
0149 return;
0150