makeCygAMap_06may2010

0001 clear all
0002 
0003 % Read in the data
0004 %d = pipe_read('05-May-2010:15:33:10','05-May-2010:16:53:10');
0005 %d = pipe_read('06-May-2010:15:30:10','06-May-2010:18:13:42');
0006 %d = pipe_read('06-May-2010:15:30:10','06-May-2010:16:50:10');
0007 %d = pipe_read('07-May-2010:15:27:12','07-May-2010:16:47:12');
0008 d = pipe_read('08-May-2010:15:24:10','08-May-2010:16:44:10');
0009 
0010 %d = pipe_read('12-May-2010:15:04:35','12-May-2010:17:48:08');
0011 %d = pipe_read('12-May-2010:15:04:35','12-May-2010:16:24:08');
0012 
0013 
0014 howMuchData = size(d.antenna0.receiver.data(:,1));
0015 
0016 d.antenna0.receiver.switchData(2401:howMuchData,:) = ...
0017     d.antenna0.receiver.switchData(1:howMuchData-2400,:);
0018 d.index.source.fast(2401:howMuchData) = ...
0019     d.index.source.fast(1:howMuchData-2400);
0020 d.index.noise.fast(2401:howMuchData) = ...
0021     d.index.noise.fast(1:howMuchData-2400);
0022 
0023 
0024 
0025 OnSourceBuffer = 1500; % 15 seconds
0026 %OnSourceBuffer = 0;
0027 
0028 xVals = 1:howMuchData;
0029 
0030 OnSource = d.index.source.fast==1 & d.index.noise.fast==0;
0031 
0032 
0033 % use find.
0034 minOn = find(OnSource,1,'first');
0035 maxOn = find(OnSource,1,'last');
0036 
0037 
0038 OnSource(1:minOn+OnSourceBuffer) = 0;
0039 OnSource(maxOn-OnSourceBuffer+1:howMuchData) = 0;
0040 
0041 
0042 % select the azimuth scan data
0043 DataOnSource = framecut(d,logical(OnSource),'receiver');
0044 
0045 
0046 
0047 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0048 % Do the r-factor correction: %
0049 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0050 [r,A] = rFactorCorrection(d,1);
0051 
0052 rCorrData = zeros(size(DataOnSource.antenna0.receiver.data,1),2);
0053 
0054 rCorrData(:,1) = -0.5*(DataOnSource.antenna0.receiver.switchData(:,2) - ... 
0055     r(1) * DataOnSource.antenna0.receiver.switchData(:,1) + ... 
0056     + DataOnSource.antenna0.receiver.switchData(:,3) - ...
0057     r(2) * DataOnSource.antenna0.receiver.switchData(:,4));
0058 
0059 rCorrData(:,2) = -0.5*(DataOnSource.antenna0.receiver.switchData(:,22) - ... 
0060     r(3) * DataOnSource.antenna0.receiver.switchData(:,21) + ... 
0061     + DataOnSource.antenna0.receiver.switchData(:,23) - ...
0062     r(4) * DataOnSource.antenna0.receiver.switchData(:,24));
0063 
0064 
0065 DataOnSource.antenna0.receiver.data(:,1) = rCorrData(:,1);
0066 DataOnSource.antenna0.receiver.data(:,6) = rCorrData(:,2);
0067 
0068 %%%%%%%%%%%%%%%%%%
0069 % Check for RFI: %
0070 %%%%%%%%%%%%%%%%%%
0071 disp('RFI Flagging...');
0072 
0073 lowRFI = kurtclip(DataOnSource.antenna0.receiver.data(:,1),100,5);
0074 
0075 
0076 
0077 %%%%%%%%%%%%%%%%%%%%%%%%%
0078 % All flags combined... %
0079 %%%%%%%%%%%%%%%%%%%%%%%%%
0080 
0081 
0082 flagArray = lowRFI; % * other flag arrays.
0083 
0084 
0085 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0086 % Remove short mid-timescale trends: %
0087 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0088 
0089 channel = 1;
0090 
0091 % First, find the individual azimuth scans:
0092 
0093 azDiff = diff(DataOnSource.antenna0.servo.apparent(:,1));
0094 convAzDiff = conv(azDiff, gauss([1/sqrt(2*pi*100) 33 10], 1:64), 'same');
0095 posScan = diff(convAzDiff > 0.007);
0096 negScan = diff(convAzDiff < -0.007);
0097 
0098 posBegin = find(posScan > 0);
0099 posEnd = find(posScan < 0);
0100 
0101 negBegin = find(negScan > 0);
0102 negEnd = find(negScan < 0);
0103 
0104 
0105 
0106 % Then do the detrending, one scan at a time!
0107 
0108 dataArray = zeros(howMuchData,1);
0109 raArray = zeros(howMuchData,1);
0110 decArray = zeros(howMuchData,1);
0111 
0112 lastPos = 0;
0113 
0114 % First for the positive scans.
0115 for m = 1:size(posBegin)
0116 
0117     % Grab the scan data and detrend it.
0118     scanData = ...
0119         DataOnSource.antenna0.receiver.data( ...
0120         posBegin(m):posEnd(m),channel);
0121 %    deScanData = detrend(scanData);
0122 %    deScanData = scanData - median(scanData);
0123     deScanData = detrend(scanData) + mean(scanData);
0124 %    deScanData = deScanData - median(deScanData);
0125 
0126     clipData = logical(deScanData > (- 4 * std(deScanData) + max(deScanData)));
0127     clipData = logical(deScanData > (- 5 * std(deScanData(clipData)) + max(deScanData)));
0128 %    clipData = logical(deScanData > (- 3 * std(deScanData(clipData)) + max(deScanData)));
0129 %    clipData = logical(deScanData > (- 3 * std(deScanData(clipData)) + max(deScanData)));
0130 %    clipData = logical(deScanData > (- std(deScanData(clipData)) + max(deScanData)));
0131 %    clipData = logical(deScanData > (- std(deScanData(clipData)) + max(deScanData)));
0132 %    clipData = logical(deScanData > (- std(deScanData(clipData)) + max(deScanData)));
0133     
0134 
0135 %    deScanData = clipData * 1;
0136 
0137     deScanData = deScanData - median(deScanData(clipData));
0138     
0139 
0140     % Grab the RA and Dec values.
0141     raData = ...
0142         DataOnSource.antenna0.servo.equa( ...
0143         posBegin(m):posEnd(m),1);
0144     decData = ...
0145         DataOnSource.antenna0.servo.equa( ...
0146         posBegin(m):posEnd(m),2);
0147     
0148     
0149     % Grab the flags and count them
0150     flagData = flagArray(posBegin(m):posEnd(m));
0151     nDe = sum(flagData);
0152     
0153     % Add the important data to our output arrays.
0154     dataArray(lastPos+1:lastPos+nDe) = deScanData(logical(flagData));
0155     raArray(lastPos+1:lastPos+nDe) = raData(logical(flagData));
0156     decArray(lastPos+1:lastPos+nDe) = decData(logical(flagData));
0157     
0158     % Keeping track...
0159     lastPos = lastPos + nDe;
0160     
0161 end
0162 
0163 % Then the negative scans.
0164 for m = 1:size(negBegin)
0165 
0166     % Grab the scan data and detrend it.
0167     scanData = ...
0168         DataOnSource.antenna0.receiver.data( ...
0169         negBegin(m):negEnd(m),channel);
0170     deScanData = detrend(scanData) + mean(scanData);
0171     
0172     clipData = logical(deScanData > (- 4 * std(deScanData) + max(deScanData)));
0173     clipData = logical(deScanData > (- 5 * std(deScanData(clipData)) + max(deScanData)));
0174 %    clipData = logical(deScanData > (- 3 * std(deScanData(clipData)) + max(deScanData)));
0175 %    clipData = logical(deScanData > (- 3 * std(deScanData(clipData)) + max(deScanData)));
0176 
0177 %    deScanData = clipData * 1;
0178 
0179     deScanData = deScanData - median(deScanData(clipData));
0180 
0181     % Grab the RA and Dec values.
0182     raData = ...
0183         DataOnSource.antenna0.servo.equa( ...
0184         negBegin(m):negEnd(m),1);
0185     decData = ...
0186         DataOnSource.antenna0.servo.equa( ...
0187         negBegin(m):negEnd(m),2);
0188 
0189     % Grab the flags and count them
0190     flagData = flagArray(negBegin(m):negEnd(m));
0191     nDe = sum(flagData);
0192     
0193     % Add the important data to our output arrays.
0194     dataArray(lastPos+1:lastPos+nDe) = deScanData(logical(flagData));
0195     raArray(lastPos+1:lastPos+nDe) = raData(logical(flagData));
0196     decArray(lastPos+1:lastPos+nDe) = decData(logical(flagData));
0197     
0198     % Keeping track...
0199     lastPos = lastPos + nDe;
0200     
0201 end
0202 
0203 
0204 dataArray = dataArray(1:lastPos);
0205 raArray = raArray(1:lastPos);
0206 decArray = decArray(1:lastPos);
0207 
0208 
0209 
0210 
0211 
0212 
0213 
0214 %%%%%%%%%%%%%%%%%
0215 % Make the Map! %
0216 %%%%%%%%%%%%%%%%%
0217 disp('Mapping...');
0218 
0219 
0220 %raSize = 20;
0221 %decSize = 40;
0222 raSize = 10; % degrees
0223 decSize = 20; % degrees
0224 raCen = 20.3; % hours
0225 decCen = 40; % degrees
0226 %pixRes = 0.1;
0227 pixRes = 0.005; % degrees
0228 beamSize = 0.75; % degrees
0229 
0230 
0231 [raRange, decRange, fluxArray] = cartplot( ...
0232     raArray, decArray, dataArray, ...
0233     [raCen decCen], [raSize decSize], pixRes, beamSize);
0234 
0235 
0236 save('08may10_cyga.mat', 'raRange', 'decRange', 'fluxArray');
0237 
0238 %cropmin = 181;
0239 %cropmax = 330;
0240 cropmin = 513;
0241 cropmax = 3584;
0242 cropFlux = fluxArray(cropmin:cropmax,cropmin:cropmax);
0243 cropRA = raRange(cropmin:cropmax);
0244 cropDec = decRange(cropmin:cropmax);
0245 
0246 figure
0247 %imagesc(raRange*12/pi,decRange*180/pi,fluxArray)
0248 imagesc(cropRA*12/pi,cropDec*180/pi,cropFlux)
0249 axis xy
0250 
0251 set(gca,'XDir','Rev')
0252 
0253 %imagesc(raRange,flipud(decRange),fluxArray)
0254 xlabel('Right Ascension [h]')
0255 ylabel('Declination [deg]')
0256 
0257 
0258 
0259