%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function outImage = ...
cartplot(ra,dec,data,mapCentre,mapSize,mapPix,beamSize)
Output are RA and Dec arrays (in radians) and an image matrix. This
represents the data resampled a projected onto a cartesian map.
Inputs are the original ra, dec, and data vectors, a 2-element vector
representing the central coords (in decimal hours and degrees), a
2-element vector representing the RA and Dec size (in degrees), and a
scalar with the pixel size (in decimal degrees). The final parameter
is a scalar with the beam size in degrees.
12-May-2010 (MAS)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0001 function [raArray, decArray, outImage] = ... 0002 cartplot(ra,dec,data,mapCentre,mapSize,mapPix,beamSize) 0003 0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0005 % 0006 % function outImage = ... 0007 % cartplot(ra,dec,data,mapCentre,mapSize,mapPix,beamSize) 0008 % 0009 % Output are RA and Dec arrays (in radians) and an image matrix. This 0010 % represents the data resampled a projected onto a cartesian map. 0011 % Inputs are the original ra, dec, and data vectors, a 2-element vector 0012 % representing the central coords (in decimal hours and degrees), a 0013 % 2-element vector representing the RA and Dec size (in degrees), and a 0014 % scalar with the pixel size (in decimal degrees). The final parameter 0015 % is a scalar with the beam size in degrees. 0016 % 0017 % 12-May-2010 (MAS) 0018 % 0019 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0020 0021 0022 0023 nRA = 2^ceil(log2(mapSize(1)/mapPix/sind(mapCentre(2)))); 0024 nDec = 2^ceil(log2(mapSize(2)/mapPix)); 0025 0026 raArray = (((1:nRA) - nRA / 2) ... 0027 * mapPix / sind(mapCentre(2)) + mapCentre(1) * 15) * pi / 180; 0028 decArray = (((1:nDec) - nDec / 2) ... 0029 * mapPix + mapCentre(2)) * pi / 180; 0030 0031 0032 sumArray = zeros(nRA,nDec); 0033 timeArray = zeros(nRA,nDec); 0034 0035 for m = 1:size(data) 0036 0037 0038 if (ra(m) > min(raArray) && ra(m) < max(raArray) && ... 0039 dec(m) > min(decArray) && dec(m) < max(decArray)) 0040 0041 [raDiff, raIndex] = min(abs(raArray - ra(m))); 0042 [decDiff, decIndex] = min(abs(decArray - dec(m))); 0043 0044 sumArray(raIndex,decIndex) = sumArray(raIndex,decIndex) - data(m); 0045 timeArray(raIndex,decIndex) = timeArray(raIndex,decIndex) + 1; 0046 0047 end 0048 0049 end 0050 0051 0052 GaussSigma = beamSize / 2.355; 0053 0054 nGauss = 2^ceil(log2(5*GaussSigma/mapPix)); 0055 Gauss1D = gauss([1 nGauss/2 GaussSigma/mapPix], 1:nGauss); 0056 GaussBeam = transpose(Gauss1D) * Gauss1D; 0057 0058 0059 sumConvArray = conv2(sumArray,GaussBeam,'same'); 0060 sumTimeArray = conv2(timeArray,GaussBeam,'same'); 0061 0062 fluxArray = sumConvArray ./ sumTimeArray; 0063 0064 % Should then trim it down to requested size... 0065 % XXX 0066 0067 outImage = transpose(fluxArray); 0068 0069 %timeArray = zeros(size(decRange,2),size(raRange,2)); 0070 %sumArray = zeros(size(decRange,2),size(raRange,2)); 0071 % 0072 % 0073 % Could probably speed this up by making big Gaussian image, and shifting 0074 % it to add. Or just do a convolution after the fact... 0075 %for i = 1:size(DataOnSource.antenna0.receiver.data(:,1)) 0076 % 0077 % decPos = DataOnSource.antenna0.servo.equa(i,2) * 180 / pi; 0078 % raPos = DataOnSource.antenna0.servo.equa(i,1) * 180 / pi; 0079 % 0080 % decGauss = gauss([1 decPos beamSize], decRange); 0081 % raGauss = gauss([1 raPos beamSizeRA], raRange); 0082 % 0083 % pixGauss = transpose(decGauss) * raGauss; 0084 % 0085 % timeArray = timeArray + pixGauss; 0086 % sumArray = sumArray - pixGauss * 0087 % DataOnSource.antenna0.receiver.data(i,1); 0088 % 0089 %end 0090 0091 0092 %f = TriScatteredInterp( ... 0093 % DataOnSource.antenna0.servo.equa(:,1), ... 0094 % DataOnSource.antenna0.servo.equa(:,2), ... 0095 % -DataOnSource.antenna0.receiver.data(:,1)); 0096 %[raRangeA, decRangeA] = meshgrid(raRange,decRange); 0097 %fluxArray = f(raRangeA,decRangeA); 0098 0099 0100 0101 return;