sourceScanMap

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function sourceScanMap(d, pixSize)

    function to take a data from a scan and turn it into a map
    d is the data structure
    pixSize is the size of the pixels to grid (0.1 degrees on a side);

  sjcm
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function [azRet elRet, dataMap, dataRms, azGrid, elGrid] = sourceScanMap(d, pixSize, doPlot)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function sourceScanMap(d, pixSize)
0006 %
0007 %    function to take a data from a scan and turn it into a map
0008 %    d is the data structure
0009 %    pixSize is the size of the pixels to grid (0.1 degrees on a side);
0010 %
0011 %  sjcm
0012 %
0013 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0014 
0015 if(nargin<2)
0016   pixSize = 0.5;
0017   doPlot = 1;
0018 end
0019 
0020 
0021 % steps in doing this:
0022 % 2. sample onto regular grid
0023 % 3. make map
0024 if(doPlot)
0025   plot(d.antenna0.receiver.utc, d.antenna0.receiver.data(:,1))
0026 end
0027 
0028 % now we calculate which azimuth and elevation the source should be at at
0029 % every time stamp in d.
0030 time = d.antenna0.receiver.utc;
0031 az_scan = d.antenna0.servo.az;
0032 el_scan = d.antenna0.servo.el;
0033 [az_src, el_src] = calcAzEl(d);
0034 
0035 deltaAz = az_scan - az_src;
0036 deltaEl = el_scan - el_src;
0037 
0038 azGrid = -1:pixSize:1;
0039 elGrid = -1:pixSize:1;
0040 
0041 %azGrid = min(deltaAz):pixSize:max(deltaAz);
0042 %elGrid = min(deltaEl):pixSize:max(deltaEl);
0043 
0044 dataVals = d.antenna0.receiver.data;
0045 dataMap  = nan(length(azGrid), length(elGrid), 6);
0046 dataRms  = nan(length(azGrid), length(elGrid), 6);
0047 % writing an inefficient loop
0048 for m=1:length(azGrid)
0049   display(sprintf('Calculating az row %i of %i', m, length(azGrid)));
0050   azInd = deltaAz>=azGrid(m)-pixSize/2 & deltaAz<azGrid(m)+pixSize/2;
0051 
0052   for n=1:length(elGrid)
0053     elInd = deltaEl>=elGrid(n)-pixSize/2 & deltaEl<elGrid(n)+pixSize/2;
0054     
0055     ind = azInd & elInd;
0056     f = find(ind);
0057     
0058     if(~isempty(f))
0059       vals = dataVals(f,:);
0060       dataMap(m,n,:) = nanmean(vals);
0061       dataRms(m,n,:) = sqrt(nanvar(vals));
0062     end
0063   end
0064 end
0065 
0066 if(doPlot)
0067   figure(1)
0068   setwinsize(gcf, 575, 1000);
0069   subplot(2,1,1)
0070   imagesc(azGrid, elGrid, -dataMap(:,:,1)./dataRms(:,:,1));
0071   colorbar('horizontal');
0072   xlabel('az offset');
0073   ylabel('el offset');
0074   title('Intensity Channel 1');
0075   
0076   subplot(2,1,2)
0077   imagesc(azGrid, elGrid, -dataMap(:,:,6)./dataRms(:,:,6));
0078   colorbar('horizontal');
0079   xlabel('az offset');
0080   ylabel('el offset');
0081   title('Intensity Channel 2');
0082   
0083   figure(2)
0084   setwinsize(gcf, 575, 1000);
0085   subplot(2,1,1)
0086   imagesc(azGrid, elGrid, -dataMap(:,:,2)./dataRms(:,:,2));
0087   colorbar('horizontal');
0088   xlabel('az offset');
0089   ylabel('el offset');
0090   title('Polarization Channel 1');
0091   
0092   subplot(2,1,2)
0093   imagesc(azGrid, elGrid, -dataMap(:,:,3)./dataRms(:,:,3));
0094   colorbar('horizontal');
0095   xlabel('az offset');
0096   ylabel('el offset');
0097   title('Polarization Channel 2');
0098   
0099   
0100   figure(3)
0101   setwinsize(gcf, 575, 1000);
0102   subplot(2,1,1)
0103   imagesc(azGrid, elGrid, -dataMap(:,:,1));
0104   colorbar('horizontal');
0105   xlabel('az offset');
0106   ylabel('el offset');
0107   title('Total Intensity 1');
0108   
0109   subplot(2,1,2)
0110   
0111   imagesc(azGrid, elGrid, -dataMap(:,:,6));
0112   colorbar('horizontal');
0113   xlabel('az offset');
0114   ylabel('el offset');
0115   title('Total Intensity 2');
0116   
0117   figure(4)
0118   setwinsize(gcf, 575, 1000);
0119   subplot(2,1,1)
0120   a = dataMap(:,:,1)./min(min(dataMap(:,:,1)));
0121   imagesc(azGrid, elGrid, log10(a));
0122   colorbar('horizontal');
0123   xlabel('az offset');
0124   ylabel('el offset');
0125   title('Total Intensity 1 (log plot)');
0126   
0127   subplot(2,1,2)
0128   a = dataMap(:,:,6)./min(min(dataMap(:,:,6)));
0129   imagesc(azGrid, elGrid, log10(a));
0130   colorbar('horizontal');
0131   xlabel('az offset');
0132   ylabel('el offset');
0133   title('Total Intensity 2 (log plot)');
0134   
0135 end
0136 
0137 % now we find the max, but only within 5 degrees of the center.
0138 indEl = abs(elGrid)<5;
0139 indAz = abs(azGrid)<5 | (abs(azGrid)>355 & abs(azGrid)<365);
0140 
0141 elCutVals = elGrid(indEl);
0142 azCutVals = azGrid(indAz);
0143 
0144 
0145 a = -dataMap(indAz, indEl, :)./dataRms(indAz, indEl, :);
0146 [x y] = find(a(:,:,1) == max(max(a(:,:,1))));
0147 
0148 display(sprintf('Peak near center at az offset of %3.2f', azCutVals(y)));
0149 display(sprintf('Peak near center at el offset of %3.2f', elCutVals(x)));
0150 
0151 azRet = azCutVals(y);
0152 elRet = elCutVals(x);
0153 
0154 
0155 return;
0156 
0157 
0158 
0159