sourceScanMap_charles

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

  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 sourceScanMap(d, pixSize)
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 end
0018 
0019 
0020 % steps in doing this:
0021 % 0. interpolate az/el data onto same size as rx data
0022 % 1. convert from az/el to RA/DEC for pointings
0023 % 2. sample onto regular grid
0024 % 3. make map
0025 
0026 % first we reshape our data
0027 d = reshapeRx(d);
0028 d = reshapeAz(d);
0029 d = interpAz(d);
0030 
0031 figure;
0032 plot(d.antenna0.receiver.utc, d.antenna0.receiver.dataNew(:,1))
0033 % now we calculate which azimuth and elevation the source should be at at
0034 % every time stamp in d.
0035 time = d.antenna0.receiver.utc;
0036 az_scan = d.antenna0.servo.interpAz;
0037 el_scan = d.antenna0.servo.interpEl;
0038 [az_src, el_src] = calcAzEl(d);
0039 
0040 %ACT - 26 Jan - added in a line to remove the az el pointing offsets we had
0041 %in
0042 az_poff = -22.42;
0043 el_poff = -1.33;
0044 deltaAz = az_scan - az_src - az_poff;
0045 deltaEl = el_scan - el_src - el_poff;
0046 
0047 figure
0048 
0049 
0050 hold
0051 t=24*60*60*(time-time(1));
0052 plot(t,d.antenna0.receiver.dataNew(:,1)+750,'k')
0053 plot(t,az_src)
0054 plot(t,az_scan-az_poff-360,'r')
0055 legend('RX','Src Az','Track Az')
0056 title('Cross Scans in Negative Azimuth Through CasA')
0057 xlabel('Time [s]')
0058 %plot(time-time(1),d.antenna0.receiver.dataNew(:,1),'r')
0059 figure
0060 plot(time-time(1),el_src,'.')
0061 
0062 %Charles Added this 26Jan
0063 figure
0064 hold
0065 plot(deltaAz,d.antenna0.receiver.dataNew(:,1),'.');
0066 title('Azimuth Offset (Negative Scan) and Intensity Channel')
0067 plot(deltaAz,d.antenna0.receiver.dataNew(:,6),'.r');
0068 legend('Intensity Ch1','Intensity Ch2');
0069 xlabel('Az Offset Deg')
0070 
0071 
0072 azGrid = min(deltaAz):pixSize:max(deltaAz);
0073 elGrid = min(deltaEl):pixSize:max(deltaEl);
0074 
0075 dataVals = d.antenna0.receiver.dataNew;
0076 dataMap  = nan(length(azGrid), length(elGrid), 6);
0077 dataRms  = nan(length(azGrid), length(elGrid), 6);
0078 % writing an inefficient loop
0079 for m=1:length(azGrid)
0080   azInd = deltaAz>=azGrid(m)-pixSize/2 & deltaAz<azGrid(m)+pixSize/2;
0081   for n=1:length(elGrid)
0082     elInd = deltaEl>=elGrid(n)-pixSize/2 & deltaEl<elGrid(n)+pixSize/2;
0083     
0084     ind = azInd & elInd;
0085     f = find(ind);
0086     
0087     if(~isempty(f))
0088       vals = dataVals(f,:);
0089       dataMap(m,n,:) = mean(vals);
0090       dataRms(m,n,:) = sqrt(var(vals));
0091     end
0092   end
0093 end
0094 
0095 
0096 % figure(1)
0097 % setwinsize(gcf, 575, 1000);
0098 % subplot(2,1,1)
0099 % imagesc(azGrid, elGrid, -dataMap(:,:,1)./dataRms(:,:,1));
0100 % colorbar('horizontal');
0101 % xlabel('az offset');
0102 % ylabel('el offset');
0103 % title('Intensity Channel 1');
0104 %
0105 % subplot(2,1,2)
0106 % imagesc(azGrid, elGrid, -dataMap(:,:,6)./dataRms(:,:,6));
0107 % colorbar('horizontal');
0108 % xlabel('az offset');
0109 % ylabel('el offset');
0110 % title('Intensity Channel 2');
0111 %
0112 % figure(2)
0113 % subplot(2,1,1)
0114 % imagesc(azGrid, elGrid, -dataMap(:,:,2)./dataRms(:,:,2));
0115 % colorbar('horizontal');
0116 % xlabel('az offset');
0117 % ylabel('el offset');
0118 % title('Polarization Channel 1');
0119 %
0120 % subplot(2,1,2)
0121 % imagesc(azGrid, elGrid, -dataMap(:,:,3)./dataRms(:,:,3));
0122 % colorbar('horizontal');
0123 % xlabel('az offset');
0124 % ylabel('el offset');
0125 % title('Polarization Channel 2');
0126 %
0127 % % now we find the max, but only within 5 degrees of the center.
0128 % indEl = abs(elGrid)<5;
0129 % indAz = abs(azGrid)<5;
0130 %
0131 % elCutVals = elGrid(indEl);
0132 % azCutVals = azGrid(indAz);
0133 %
0134 %
0135 % a = -dataMap(indAz, indEl, :)./dataRms(indAz, indEl, :);
0136 % [x y] = find(a(:,:,1) == max(max(a(:,:,1))));
0137 %
0138 % display(sprintf('Peak near center at az offset of %3.2f', azCutVals(y)));
0139 % display(sprintf('Peak near center at el offset of %3.2f', elCutVals(x)));
0140 %
0141 %
0142 
0143 return;
0144 
0145 
0146 
0147