Description of quick_dirty

0001 % After reading in data then go through each cell to get a raster image.
0002 % Stope when you see a cell saying 'OTHER STUFF'
0003 % Need to addpath=('./Angelas_Raster_Code')  to use nanstd2.m
0004 
0005 %moon
0006 %d = read_arcSouth('10-Jun-2014:18:53:47','10-Jun-2014:20:23:05')
0007 
0008 LL = [d.antenna0.roach1.LL,fliplr(d.antenna0.roach2.LL)]; % **Check which order the roaches are in **
0009 RR = [d.antenna0.roach1.RR,fliplr(d.antenna0.roach2.RR)]; % **Check which order the roaches are in **
0010 L1 = [d.antenna0.roach1.load1,fliplr(d.antenna0.roach2.load1)]; % **Check which order the roaches are in **
0011 L2 = [d.antenna0.roach1.load2,fliplr(d.antenna0.roach2.load2)]; % **Check which order the roaches are in **
0012 
0013 QQ = [d.antenna0.roach1.Q,fliplr(d.antenna0.roach2.Q)];
0014 UU = [d.antenna0.roach1.U,fliplr(d.antenna0.roach2.U)];
0015 
0016 
0017 
0018 %%
0019 % Quick deglitcher - pick out time regions from plot
0020 % Use cursor to identify ranges of data to flag. Pick points in pairs.Press
0021 % return when done. It'll then plot the unflagged data to see if you are
0022 % happy.
0023 
0024 
0025 happy='n'
0026 while(happy=='n')
0027 
0028     figure
0029 plot(mean(LL'))
0030 hold all
0031 plot(mean(RR'))
0032 disp('Now pick time ranges that you want to exclude from the data. Press return when done')
0033 
0034 [x,y]=ginput
0035 
0036 % Make that data NaNs
0037 for i=1:2:length(x)
0038     LL(x(i):x(i+1),:)=NaN;
0039     RR(x(i):x(i+1),:)=NaN;
0040     QQ(x(i):x(i+1),:)=NaN;
0041     UU(x(i):x(i+1),:)=NaN;
0042     L1(x(i):x(i+1),:)=NaN;
0043     L2(x(i):x(i+1),:)=NaN;
0044 end
0045 figure
0046 plot(mean(LL'))
0047 hold all
0048 plot(mean(RR'))
0049 happy = input('Happy? (y/n) ', 's');
0050 end
0051 
0052 % Now check Q and U
0053 happy='n'
0054 while(happy=='n')
0055 
0056 figure
0057 plot(mean(QQ'))
0058 hold all
0059 plot(mean(UU'))
0060 disp('Now pick time ranges that you want to exclude from the data. Press return when done')
0061 
0062 [x,y]=ginput
0063 
0064 % Make that data NaNs
0065 for i=1:2:length(x)
0066     LL(x(i):x(i+1),:)=NaN;
0067     RR(x(i):x(i+1),:)=NaN;
0068     QQ(x(i):x(i+1),:)=NaN;
0069     UU(x(i):x(i+1),:)=NaN;
0070     L1(x(i):x(i+1),:)=NaN;
0071     L2(x(i):x(i+1),:)=NaN;
0072 end
0073 figure
0074 plot(mean(QQ'))
0075 hold all
0076 plot(mean(UU'))
0077 happy = input('Happy? (y/n) ', 's');
0078 end
0079 
0080 % Also flag very bad Q,U values which are less than -6.5e4
0081 dodgyU=(UU<=-6e4);
0082 UU(dodgyU)=NaN;
0083 dodgyQ=(QQ<=-6e4);
0084 QQ(dodgyQ)=NaN;
0085 %%
0086 goodchans = [8:61,68:125];
0087 mask=ones(1,128);
0088 mask(goodchans) = 0;
0089 mask=logical(mask);
0090 mask = ~mask;
0091 %%
0092 % Plot out std and imagesc of remaining good data - ask for more bad
0093 % channels to mask
0094 % Identify the channels to flag from the plots then enter them as e.g. [1 2 3]
0095 % You will be asked if happy -  I have got the happy question the wrong side
0096 % of a plot so you need to say 'n' to check what has been applied. If you
0097 % are then in fact happy, just give a channel that has already been flagged
0098 % and then answer 'y' to the happy question!!  - it is late :-)
0099 
0100 addpath('./Angelas_Raster_Code')
0101 happy='n'
0102 while(happy=='n')
0103     
0104 LL(:,~mask)=NaN;
0105 RR(:,~mask)=NaN;
0106 QQ(:,~mask)=NaN;
0107 UU(:,~mask)=NaN;
0108 close all
0109 figure
0110 subplot(2,2,1)
0111 plot(nanstd2(RR))
0112 title('RR')
0113 subplot(2,2,2)
0114 plot(nanstd2(LL))
0115 title('LL')
0116 subplot(2,2,3)
0117 plot(nanstd2(UU))
0118 title('UU')
0119 subplot(2,2,4)
0120 plot(nanstd2(QQ))
0121 title('QQ')
0122 badchans = input('identify other bad chans: e.g.[75 96 24] ');
0123 mask(badchans)=0;
0124 
0125 happy = input('Happy? (y/n) ', 's');
0126 close all
0127 end
0128 
0129 %%
0130 % Gets the data into North format
0131 
0132 d.antenna0.receiver.data(:,1)=  nanmean((LL(:,mask)-L2(:,mask))'); %LL
0133 d.antenna0.receiver.data(:,6)=  nanmean((RR(:,mask)-L1(:,mask))'); %RR
0134 d.antenna0.receiver.data(:,2)=  nanmean(QQ(:,mask)'); %QQ
0135 d.antenna0.receiver.data(:,3)=  nanmean(UU(:,mask)'); %UU
0136 d.antenna0.receiver.data(:,4) = d.antenna0.receiver.data(:,2);
0137 d.antenna0.receiver.data(:,5) = d.antenna0.receiver.data(:,3);
0138 d.antenna0.receiver.utc=d.antenna0.roach1.utc;
0139 %%
0140 %necessary to prepare for plotting raster data
0141  nchan=4
0142  plot_labels = {'LL','RR','QQ','UU'};
0143  chan_order=[1 6 2 3];
0144 
0145 trigpoint =0
0146 
0147 if (trigpoint==1)
0148     
0149     % Calculate the az and El offset form the source
0150     %az_trig = 17-1
0151     %el_trig = 11.5-4.5
0152     az_trig=16
0153     el_trig=6.5
0154     
0155     
0156     d.azOffSave = (d.antenna0.servo.apparent(:,1)-az_trig).*cosd(d.antenna0.servo.apparent(:,2));
0157     d.elOffSave = d.antenna0.servo.apparent(:,2)-el_trig;
0158     d.angle = sqrt(d.elOffSave.^2 + d.azOffSave.^2);
0159     
0160 else
0161     azApp = interp1(d.antenna0.tracker.utc, ...
0162         d.antenna0.tracker.horiz_topo(:,1), d.antenna0.roach2.utc);
0163     azOffSave = azApp - d.antenna0.servo.apparent(:,1);
0164     azOffSave = -azOffSave;
0165     azOffSave = azOffSave.*cos(d.antenna0.servo.el*pi/180);% Scale azimuth by cos(elevation)???CHECK
0166     
0167     elApp = interp1(d.antenna0.tracker.utc, ...
0168         d.antenna0.tracker.horiz_topo(:,2), d.antenna0.roach2.utc);
0169     elOffSave = elApp - d.antenna0.servo.apparent(:,2);
0170     elOffSave = -elOffSave;
0171     d.angle = sqrt(elOffSave.^2 + azOffSave.^2);%
0172     d.elOffSave = elOffSave;
0173     d.azOffSave= azOffSave;
0174     
0175 end
0176 
0177 draster.antenna0.receiver.data(:,1)=d.antenna0.receiver.data(d.index.beammap.fast,1);
0178 draster.antenna0.receiver.data(:,6)=d.antenna0.receiver.data(d.index.beammap.fast,6);
0179 draster.antenna0.receiver.data(:,2)=d.antenna0.receiver.data(d.index.beammap.fast,2);
0180 draster.antenna0.receiver.data(:,3)=d.antenna0.receiver.data(d.index.beammap.fast,3);
0181 
0182 draster.elOffSave=d.elOffSave(d.index.beammap.fast);
0183 draster.azOffSave=d.azOffSave(d.index.beammap.fast);
0184 %%
0185  
0186 % Now analyse the Raster data and make a map
0187 
0188 pixel_size = 0.1; % degree
0189 
0190 azmin = -10;
0191 azmax = +10;
0192 elmin = -10;
0193 elmax = +10;
0194 nanflag = isnan(draster.antenna0.receiver.data(:,1));
0195 flag = ~nanflag;
0196 colLR = 1e5;
0197 colQU = 1e2;
0198 % Set colour scale mapping
0199 %clim = [0 3.5;-0.2 0.2;-0.2 0.2;0 0.2]
0200 
0201 %for m= 1:length(s)
0202  %   n = n+1;
0203  %   ind = zeros(size(FLAG));
0204  %   ind(s(m):e(m)) = 1;
0205  %   ind = logical(ind);
0206  %   flagvector3 = (~nanflag & ~FLAG & ind);
0207    
0208     for i=1:4
0209         j=chan_order(i);
0210         [TauAmap{i}, xidx, yidx, TauAmapIdx{i}] = bin_quickly2d(azmin, azmax, elmin, elmax, pixel_size,...
0211             draster.azOffSave(flag), draster.elOffSave(flag),draster.antenna0.receiver.data(flag,j));
0212         figure(3)
0213         subplot(2,2,i)
0214         if i<3
0215         clim=[0.5e4 colLR];
0216     else
0217         clim=[-colQU colQU];
0218     end
0219         %imagesc(fliplr(xidx),yidx,TauAmap{i},clim(i,:))
0220         %clim=[5e3 1.5e4]
0221        imagesc(fliplr(xidx),fliplr(yidx),TauAmap{i})
0222         %imagesc(fliplr(xidx),fliplr(yidx),TauAmap{i},clim)
0223         
0224         %title([source_name,' map:', num2str(i)])
0225         colormap('default')
0226         colorbar
0227         title(plot_labels{i})
0228         axis square
0229         figure(4)
0230         subplot(2,2,i)
0231         %clines = [1:-0.01:0]*max(max(TauAmap{i}));
0232         contour(fliplr(xidx),yidx,TauAmap{i})
0233         %contour(fliplr(xidx),yidx,TauAmap{i},clines)
0234         title(plot_labels{i})
0235         axis square
0236         
0237     end
0238 
0239  %% OTHER STUFF -STOP !!
0240 LL(:,mask)=NaN;
0241 RR(:,mask)=NaN;
0242 QQ(:,mask)=NaN;
0243 UU(:,mask)=NaN;
0244 
0245 mask=flagmask(mask, LL, RR, QQ, UU);
0246 display('Flagging channels: ')
0247 find(mask)
0248 %% OTHER STUFF - STOP!
0249 % Check rms on a bit of data
0250 std_start=2.156e4
0251 std_stop=4.479e4
0252 
0253 stdQQ = std(QQ(std_start:std_stop,:));
0254 stdUU = std(UU(std_start:std_stop,:));
0255 stdLL = std(LL(std_start:std_stop,:)-L2(std_start:std_stop,:));
0256 stdRR = std(RR(std_start:std_stop,:)-L1(std_start:std_stop,:));
0257 
0258 stdI1 = std(d.antenna0.receiver.data(std_start:std_stop,1));
0259 stdI2 = std(d.antenna0.receiver.data(std_start:std_stop,6));
0260 
0261 % Get size of load step
0262 load_start1 = 2.96e4  
0263 load_stop1  = 4.462e4
0264 mean1_L1    = mean(L1(load_start1:load_stop1,:))
0265 mean1_L2     = mean(L2(load_start1:load_stop1,:))
0266 
0267 
0268 load_start2 = 7.02e4
0269 load_stop2  = 8.057e4
0270 mean2_L1    = mean(L1(load_start2:load_stop2,:))
0271 mean2_L2     = mean(L2(load_start2:load_stop2,:))
0272 
0273 load_step_L1 = mean2_L1 - mean1_L1;
0274 load_step_L2 = mean2_L2 - mean1_L2;
0275 
0276 
0277 loadK = 1.6920  % 2K step in load equates to about 1.6920K real temp (see spreadsheet)
0278 
0279 loadK_units_L1 = loadK./load_step_L1;
0280 loadK_units_L2 = loadK./load_step_L2;
0281 
0282 rmsQQ = loadK_units_L1.*stdQQ;
0283 rmsUU = loadK_units_L1.*stdUU;
0284 rmsLL = loadK_units_L2.*stdLL;
0285 rmsRR = loadK_units_L1.*stdRR;
0286 plot(rmsQQ)
quick_dirty_south

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
