plotMainsReg

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 function data=plotMainsReg(d, flags, xc, txt, ax, prop, w, chunk, h, ind1)

 d       - data structure
 flag    - flag structure
 xc      - cell element of x and y data
           x{odd}  - x axis data
           x{even} - y axis data Nx1 or Nx6 or Nx24
 txt     - axis text
           txt{1} - title
           txt{2} - x axis
           txt{3} - y axis
           txt[4} - overall title
           txt{5} - legend
 ax         - manual axis settings.  [xmin xmax ymin ymax]
 prop    - initialize properties of plotting
           p.style = {'flag','feature'}
           p.swap  = swap to display plots in time.
           p.featmask= [bitmask for style='feature']


 internal inputs for recursive calling
 w       - waitforbuttonpress result
 chunk   - time chunk to start on
 h       - plot handle
 ind1    - indices
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function data=plotMainsReg(d, flags, xc, txt, ax, prop, w, chunk, h, ind1)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 % function data=plotMainsReg(d, flags, xc, txt, ax, prop, w, chunk, h, ind1)
0005 %
0006 % d       - data structure
0007 % flag    - flag structure
0008 % xc      - cell element of x and y data
0009 %           x{odd}  - x axis data
0010 %           x{even} - y axis data Nx1 or Nx6 or Nx24
0011 % txt     - axis text
0012 %           txt{1} - title
0013 %           txt{2} - x axis
0014 %           txt{3} - y axis
0015 %           txt[4} - overall title
0016 %           txt{5} - legend
0017 % ax         - manual axis settings.  [xmin xmax ymin ymax]
0018 % prop    - initialize properties of plotting
0019 %           p.style = {'flag','feature'}
0020 %           p.swap  = swap to display plots in time.
0021 %           p.featmask= [bitmask for style='feature']
0022 %
0023 %
0024 % internal inputs for recursive calling
0025 % w       - waitforbuttonpress result
0026 % chunk   - time chunk to start on
0027 % h       - plot handle
0028 % ind1    - indices
0029 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0030 
0031 if(exist('w')~=1)
0032   w = 2;
0033 end
0034 
0035 % check for axis setting
0036 if (~exist('ax'))
0037   ax=[];
0038 end
0039 
0040 % set the plot properties.
0041 if (exist('prop'))
0042   p=setProperties(prop);
0043 else
0044   p=setProperties([]);
0045 end
0046 
0047 % num is the number of pages to be displayed.
0048 % matrix is the size of each of the pages to be displayed.
0049 num = 1;
0050 plotsPerPage = 6;
0051 matrix = [3 2];
0052 
0053 plotoptions = {'b', 'g', 'r', 'c', 'm', 'y', 'k', ...
0054     'b.', 'g.', 'r.', 'c.', 'm.', 'y.', 'k.', ...    
0055     'bo', 'go', 'ro', 'co', 'mo', 'yo', 'ko', ...    
0056     'bx', 'gx', 'rx', 'cx', 'mx', 'yx', 'kx'};
0057 
0058 data = d.antenna0.receiver.dataT;
0059 data(d.flags.fast) = nan;
0060 timeVal = d.antenna0.receiver.utc*24*60*60;
0061 timeVal_1 = circshift(timeVal,1);
0062 
0063 pageNum=1;
0064 fmode=logical(0);
0065 
0066 while(1)
0067   if (w~=0 & w~=1)
0068     figure(1)
0069     setwinsize(gcf, 1000, 750); 
0070     clear leg;
0071     for plotNum=1:plotsPerPage
0072       indices = d.correction.mains.indices(plotNum,:);
0073       subplot(matrix(1),matrix(2),plotNum)
0074       % plot data
0075       optionIndex = 1;
0076       for mm=1:size(d.correction.mains.factors,2)
0077     if(mm==1)
0078       % plot uncorrected data
0079       [boo, coo, h(pageNum, plotNum,optionIndex)] = plot_powerspec(d, ...
0080           data(:,plotNum), indices(1), indices(2), ...
0081           0, plotoptions{optionIndex});
0082       hold on;
0083       leg{optionIndex} = 'Orig';
0084       optionIndex = optionIndex+1;
0085     end
0086     clear boo;
0087     clear coo;
0088     
0089     
0090     % remove the harmonic of the mains
0091     best_freq = d.correction.mains.factors(plotNum,mm,1);
0092     best_phase= d.correction.mains.factors(plotNum,mm,2);
0093     best_amp  = d.correction.mains.factors(plotNum,mm,3);
0094     
0095      best_fit_sin = cos((timeVal*2*pi*(best_freq))+best_phase) ...
0096                     - cos((timeVal_1*2*pi*(best_freq))+best_phase);
0097                 
0098      best_fit_sin = best_fit_sin * best_amp;
0099     
0100     %Old version
0101     %best_fit_sin = sin(timeVal*2*pi*best_freq + ...
0102     %    best_phase)*best_amp;
0103     
0104     
0105     data(:,plotNum) = data(:,plotNum)-best_fit_sin;
0106     
0107     % plot the one with it removed
0108     [boo, coo, h(pageNum, plotNum,optionIndex)] = plot_powerspec(d, ...
0109         data(:,plotNum), indices(1), indices(2), ...
0110         0, plotoptions{optionIndex});
0111     hold on;
0112     leg{optionIndex} = sprintf('Harmonic %u', mm);
0113     optionIndex = optionIndex+1;    
0114       end
0115       xlabel(txt{2});
0116       ylabel(txt{3});
0117       thisTitle = sprintf('%s %u', txt{1}{1}, plotNum);
0118       title(thisTitle);
0119     end
0120     
0121     legend(leg,'Location','Best');
0122     gtitle('Mains Hum Removal');
0123     w=waitforbuttonpress;
0124     if (~w)
0125       % mouse click -> zooming
0126 %      flags = plotd(d,flags,xc,txt,ax,p,w,[],h,pageNum);
0127     else
0128       % evaluate key stroke
0129       % and get next plotting index
0130       [pageNum,format,p]=evalkey(d,pageNum,plotsPerPage,'matrix',p,num);
0131       % negative index means exit
0132       if(pageNum<0)
0133     return
0134       end
0135       if (strcmp(format,'swap'))
0136     numbuf=num;
0137     num=matrix;
0138     matrix=numbuf;
0139     pageNum=1;
0140     p.swap=~p.swap;
0141     clear h;
0142     close(1)
0143       end
0144     end
0145     % set w so we enter this statement again
0146     w=2;
0147   else 
0148 %    if (~w)
0149 %      % w=0.
0150 %      % zoom to plot that was clicked on
0151 %      child=get(gca,'Children');
0152 %      ind2 = find(h(ind1,:)==child(1));
0153 %      if (isempty(ind2))
0154 %    if (isempty(chunk))
0155 %      disp(' Unknown plot.  Default to channel 1');
0156 %      ind2=1;
0157 %    else
0158 %      ind2=chunk;
0159 %    end
0160 %      end
0161 %    else
0162 %      % evalate key stroke
0163 %      % and get band index
0164 %      [ind2,format,p] = evalkey(d,chunk,plotsPerPage,'zoom',p,num);
0165 %      % index==-1 means exit
0166 %      % index==-2 means enter flag mode
0167 %      if(ind2==-1)
0168 %    break;
0169 %      elseif (ind2==-2)
0170 %    fmode = logical(1);
0171 %    ind2 = chunk;
0172 %      end
0173 %    end
0174 %
0175   end
0176 end
0177 
0178     
0179 
0180 
0181 %%%%%%%%%%%%%%%%%%%%%%%%%
0182 function p=setProperties(p)
0183 
0184 if (isempty(p))
0185   p.style='flag';
0186   p.swap=0;
0187   p.featflag=0;
0188   return
0189 end
0190   
0191 if (isfield(p,'style'))
0192 
0193 else
0194   p.style='flag';  
0195   p.featflag=0;
0196 end
0197 
0198 if (~isfield(p,'swap'))
0199   p.swap=0;
0200 end
0201 
0202 return;
0203 
0204 
0205 
0206 
0207 %%%%%%%%%%%%%%%%%%%%%%%%%
0208 function [m,n]=swapIndex(swap,i,j)
0209 
0210 if (swap)
0211   m=j; n=i;
0212 else
0213   m=i; n=j;
0214 end
0215 
0216 
0217 %%%%%%%%%%%%%%%%%%%%%%%%%
0218 function h=plotflag(x,allflags,ax,pageNum,plotNum,txt,w,num,matrix,swap, multLinesPerPlot)
0219 % get auto flags and user flags
0220 [flag,autoflag,scriptflag]=cloneflag(allflags);
0221 
0222 % if x has a fourth dimension, plot that first
0223 if(size(x,1)==4)
0224   plot(x{1,pageNum},x{4,pageNum}(:,plotNum),'c.','MarkerSize',10);
0225   hold on;
0226 end
0227 
0228 
0229 % plot data
0230 h=plot(x{1,pageNum},x{2,pageNum}(:,plotNum),'b.','MarkerSize',10);
0231 manaxis(ax,x{2,pageNum}(:,plotNum))
0232 hold on;
0233 
0234 plotdescription(txt,plotNum,pageNum,num,matrix,swap);
0235 
0236 
0237 % check for manual axis settings
0238 manaxis(ax,x{2,pageNum}(:,plotNum));
0239 
0240 % plot order is IMPORTANT!
0241 numpts = size(x{1,pageNum},1);
0242   
0243 % plot auto flags
0244 if (~isempty(autoflag))
0245   f = find(autoflag{1,pageNum}(:,plotNum));
0246   if (~isempty(f))
0247     h=plot(x{1,pageNum}(f),x{2,pageNum}(f,plotNum),'m.','MarkerSize',10);
0248   end
0249 end           
0250 
0251 % plot script flags
0252 if (~isempty(scriptflag))
0253   f = find(scriptflag{1,pageNum}(:,plotNum));
0254   if (~isempty(f))
0255     h=plot(x{1,pageNum}(f),x{2,pageNum}(f,plotNum),'g.','MarkerSize',10);
0256   end
0257 end           
0258 
0259 % plot flags
0260 if (~isempty(flag)) 
0261   % always plotting fast sampled data.
0262   f = find(flag{1,pageNum}(:,plotNum));
0263   if (~isempty(f))
0264     h=plot(x{1,pageNum}(f),x{2,pageNum}(f,plotNum),'r.','MarkerSize',10);
0265   end
0266 end
0267 
0268 hold off
0269     
0270 return;
0271 
0272 
0273 %%%%%%%%%%%%%%%%%%%%%%%%%
0274 function [h, leg] = plotfeat(x,allflags,ax,pageNum,plotNum,txt,w,num,matrix,swap)
0275 
0276 % get the features to plot
0277 % we only have {source, calibrator, abscal, blank, skydip, noise_event, and noise}
0278 
0279 % get auto flags and user flags
0280 [flag,autoflag,scriptflag]=cloneflag(allflags);
0281 
0282 % get the color and legend
0283 [clr, legall] = getcolor('feature');
0284 
0285 % we only need to plot the ones with legall
0286 index = 1;
0287 leg = [];
0288 for m=1:length(legall)
0289   txt1 = sprintf('a = x{1,pageNum}(x{3,pageNum}.%s);', legall{m});
0290   txt2 = sprintf('b = x{2,pageNum}(x{3,pageNum}.%s, plotNum);', legall{m});
0291   eval(txt1); eval(txt2);
0292   
0293   % do not plot if there are no data points to plot
0294   if(sum(a)>0)
0295     txtPlot = sprintf('h = plot(a, b, ''.'', ''markerfacecolor'', clr(m,:), ''markeredgecolor'', clr(m,:), ''MarkerSize'',10);');
0296     eval(txtPlot);
0297     hold on;
0298     leg{index} = legall{m};
0299     index = index+1;
0300   end
0301 end  
0302 
0303 
0304 % plot auto flags
0305 if (~isempty(autoflag))
0306   f = find(autoflag{1,pageNum}(:,plotNum));
0307   if (~isempty(f))
0308     h=plot(x{1,pageNum}(f),x{2,pageNum}(f,plotNum),'m.','MarkerSize',10);
0309     leg{index} = 'auto flag';
0310     index = index+1;
0311   end
0312 end           
0313 
0314 % plot script flags
0315 if (~isempty(scriptflag))
0316   f = find(scriptflag{1,pageNum}(:,plotNum));
0317   if (~isempty(f))
0318     h=plot(x{1,pageNum}(f),x{2,pageNum}(f,plotNum),'g.','MarkerSize',10);
0319     leg{index} = 'script flag';
0320     index = index+1;
0321   end
0322 end           
0323 
0324 % plot flags
0325 if (~isempty(flag)) 
0326   % always plotting fast sampled data.
0327   f = find(flag{1,pageNum}(:,plotNum));
0328   if (~isempty(f))
0329     h=plot(x{1,pageNum}(f),x{2,pageNum}(f,plotNum),'r.','MarkerSize',10);
0330     leg{index} = 'user flag';
0331     index = index+1;
0332   end
0333 end
0334       
0335 
0336 
0337 manaxis(ax,x{2,pageNum}(:,plotNum))
0338 plotdescription(txt,plotNum,pageNum,num,matrix,swap);
0339 
0340 hold off
0341     
0342 return;
0343 
0344 
0345 %%%%%%%%%%%%%%%%%%%%%%%%%
0346 function plotdescription(txt,n,m,num,matrix,swap)
0347 
0348 if (matrix(1)==2 & matrix(2)==3)
0349   % plot of channels
0350   title(sprintf('Channel: %u', n));
0351 elseif (matrix(1)==4 & matrix(2)==6)
0352   % plot of switch states
0353   title(sprintf('Switch State Number: %u Ant: %u',n));
0354 elseif (num==6)
0355   % channels plotted in time
0356   title(sprintf('Channel Number: %u', m));
0357 elseif(num==24)
0358   % phase switch in time
0359   title(sprintf('Switch State Number: %u', m));
0360 elseif(matrix(1) == 1 & matrix(2)==2)
0361   % opacity plots
0362   title(sprintf('Opacity for Channel: %u', m));
0363 end
0364 
0365 if(~isempty(txt))
0366   xlabel(char(txt{2}));
0367   ylabel(char(txt{3}));
0368 end
0369 
0370 return;
0371 
0372 
0373 %%%%%%%%%%%%%%%%%%%%%%%%%
0374 function manaxis(ax,y)
0375 
0376 if (~isempty(ax))
0377   ax4orig = ax(4);
0378   ymax=max(y)*1.1;
0379   if (max(y)>ax(4))
0380     % if max y is larger than mean
0381     ax(4)=ymax;
0382   elseif (8*nanmean(y)<ax(4) & ax(3)>=0)
0383     % if global mean (ax(4)) is much greater than the
0384     % local y data, rescale
0385     % but if ax(3)<0, that means we're looking at phase
0386     % so don't do anything.
0387     ax(4)=ymax;
0388   elseif (isnan(ax(4)) | isnan(ax(3)))
0389     if (isnan(ax(4)))
0390       ax(4)=1;
0391     end
0392     if (isnan(ax(3)))
0393       ax(3)=0;
0394     end
0395   end
0396   % this is in place until mike comes up with a better way to do
0397   % this - all the previous stuff at times resets the value of
0398   % ax(4) so that it is below that of ax(3) - i don't feel like
0399   % trying to figure out mike's logic, so here's a fix;
0400   if(ax(3) >= ax(4))
0401     if(ax4orig>ax(3))
0402           ax(4) = ax4orig;
0403     else
0404       ax(4) = ax(3)+1;
0405     end
0406   end
0407   axis(ax);
0408 end
0409 
0410 return;
0411 
0412 %%%%%%%%%%%%%%%%%%%%%%%%%
0413 function [inc,format,p]=evalkey(d,inc,num,type,p,numPages)
0414 
0415 format='default';
0416 ch=get(gcf,'CurrentCharacter');
0417 switch ch
0418   case 'p'
0419     inc=inc-1;
0420     
0421   case 'n'
0422     inc=inc+1;
0423     
0424   case 'q'
0425     if (strcmp(type,'zoom'))
0426       close(2)
0427     else
0428       close(1)
0429     end
0430     inc=-1;
0431     
0432  otherwise
0433     inc=inc;
0434 end
0435 
0436 if (strcmp(type,'matrix'))
0437   if(inc>numPages)
0438     inc = 1;
0439   elseif(inc==0)
0440     inc=numPages;
0441   end
0442 else
0443   if (inc>num)
0444     inc=1;
0445   elseif (inc==0)
0446     inc=num;
0447   end
0448 end
0449 
0450   
0451 
0452 return;
0453 
0454 %%%%%%%%%%%%%%%%%%%%%%%%%
0455 function printhelp()
0456 
0457 disp(' ')
0458 disp('******************************')
0459 disp('        Plotting Help         ')
0460 disp('******************************')
0461 disp(' ');
0462 disp('IN PLOT WINDOW (MATRIX and ZOOM)')
0463 disp(' ')
0464 disp('h - command menu')
0465 disp('p - Previous plots')
0466 disp('n - Next plots')
0467 disp('q - exit plot')
0468 disp(' ')
0469 disp(' ')
0470 disp('*******************WARNING*********************')
0471 disp('*******************WARNING*********************')
0472 disp('*******************WARNING*********************')
0473 disp('*******************WARNING*********************')
0474 disp('*******************WARNING*********************')
0475 disp('Do not exit plot by manually killing the window.')
0476 disp('This will cause the program to crash!')
0477 disp('***********************************************')
0478 disp(' ')
0479 
0480 
0481 return;
0482 
0483 
0484 %%%%%%%%%%%%%%%%%%%%%%%%%
0485 function [clr,leg]=getcolor(type)
0486 
0487 clr=[ 0.2     0         0;
0488     1         0.66      0;
0489     1         0         1;
0490     0.5       0.2       1;
0491     0.2       0.2       0.8;
0492     1         0.5       0.12;
0493     0.1       0.5       0.9;
0494     0.75      1         1;
0495     0.2       0.5       0.66;
0496     0.9       0.8       0.7;
0497     0.6       0         0.1;
0498     0         0         1
0499     0.5869    0.1909    0.3461
0500     0.0576    0.8439    0.1660
0501     0.3676    0.1739    0.1556
0502     0.6315    0.1708    0.1911
0503     0.7176    0.9943    0.4225
0504     0.6927    0.4398    0.8560
0505     0.0841    0.3400    0.4902
0506     0.4544    0.3142    0.8159
0507     0.4418    0.3651    0.4608
0508     0.3533    0.3932    0.4574
0509     0.1536    0.5915    0.4507
0510     0.6756    0.1197    0.4122
0511     0.6992    0.0381    0.9016
0512     0.7275    0.4586    0.0056
0513     0.4784    0.8699    0.2974
0514     0.5548    0.9342    0.0492
0515     0.1210    0.2644    0.6932
0516     0.4508    0.1603    0.6501
0517     0.7159    0.8729    0.9830
0518     0         0         0];
0519 
0520 
0521 switch type
0522   case 'feature'
0523     leg = {'source', 'calibrator', 'abscal', 'blank', 'skydip', ...
0524     'noise_event', 'noise'};
0525     
0526   case 'bitmask'
0527     leg={'frame','receiver', 'tracking', 'tracker', 'elevation', ...
0528     'continuous', 'fifo', 'clock', '1pps', 'cryo/servo', 'backend time', 'shortfall', 'sun', 'moon'};
0529 end
0530 
0531 return;
0532 
0533