updateAlphaDatabase

 function updateAlphaDatabase(currentTime)
 function updateAlphaDatabase(currentTime,MODE)
 OR
 function updateAlphaDatabase(currentTime,scanLength)
 function updateAlphaDatabase(currentTime,scanLength,MODE)
 OR
 function updateAlphaDatabase(startTime,endTime)
 function updateAlphaDatabase(startTime,endTime,MODE)
 OR
 function updateAlphaDatabase(startTime,endTime,scanlength)
 function updateAlphaDatabase(startTime,endTime,scanlength,MODE)

 This function updates the alphaDatabase.txt file. If only the current time is
 given then it reads in all the data between the final database entry and
 the given argument, currentTime (MJD). It calculates the alpha values in 
 this time, and allows you to graphically flag the bad alpha values. It 
 then give you the option of updating the database file with the new 
 values.

 You can specify the mode, either as 'INTERACTIVE' or 'NONINTERACTIVE'. In
 the interactive mode, the data will be automatically flagged, but you
 will be asked to click through a few screens of alpha values to identify
 and excise any bad ones missed by the automatic detection algorithms. In
 non-interactive mode, the automatically flagged alpha values will be
 written to the database without human oversight.

 If no mode is specified, then we run in 'INTERACTIVE' mode.

 If two MJD times are given then it interprets this as the start time and
 end to scan over. It will replace all the database values within that
 time period, after flagging of course.

 You may optionally specify the length of time to read in, in hours. The
 default is 1 hours.

 ogk, 13 March 2012

 Edited: OGK, 19 November 2012
         Edited to work with monthly database files instead of one large
         database file.

0001 function updateAlphaDatabase(varargin)
0002 % function updateAlphaDatabase(currentTime)
0003 % function updateAlphaDatabase(currentTime,MODE)
0004 % OR
0005 % function updateAlphaDatabase(currentTime,scanLength)
0006 % function updateAlphaDatabase(currentTime,scanLength,MODE)
0007 % OR
0008 % function updateAlphaDatabase(startTime,endTime)
0009 % function updateAlphaDatabase(startTime,endTime,MODE)
0010 % OR
0011 % function updateAlphaDatabase(startTime,endTime,scanlength)
0012 % function updateAlphaDatabase(startTime,endTime,scanlength,MODE)
0013 %
0014 % This function updates the alphaDatabase.txt file. If only the current time is
0015 % given then it reads in all the data between the final database entry and
0016 % the given argument, currentTime (MJD). It calculates the alpha values in
0017 % this time, and allows you to graphically flag the bad alpha values. It
0018 % then give you the option of updating the database file with the new
0019 % values.
0020 %
0021 % You can specify the mode, either as 'INTERACTIVE' or 'NONINTERACTIVE'. In
0022 % the interactive mode, the data will be automatically flagged, but you
0023 % will be asked to click through a few screens of alpha values to identify
0024 % and excise any bad ones missed by the automatic detection algorithms. In
0025 % non-interactive mode, the automatically flagged alpha values will be
0026 % written to the database without human oversight.
0027 %
0028 % If no mode is specified, then we run in 'INTERACTIVE' mode.
0029 %
0030 % If two MJD times are given then it interprets this as the start time and
0031 % end to scan over. It will replace all the database values within that
0032 % time period, after flagging of course.
0033 %
0034 % You may optionally specify the length of time to read in, in hours. The
0035 % default is 1 hours.
0036 %
0037 % ogk, 13 March 2012
0038 %
0039 % Edited: OGK, 19 November 2012
0040 %         Edited to work with monthly database files instead of one large
0041 %         database file.
0042 %
0043 
0044 % What scenario are we in?
0045 switch size(varargin,2)
0046     case 1
0047         % currentTime
0048         disp('updateAlphaDatabase:: Scenario 1: append new entries to the database')
0049         SCENARIO = 1;
0050         MODE = 'INTERACTIVE';
0051         currentTime = varargin{1};
0052         mjdstop = currentTime;
0053         tscans = 1; % hours
0054     case 2
0055         % The arguments are either:
0056         % currentTime,scanLength OR
0057         % startTime,endTime OR
0058         % currentTime,MODE
0059 %        if ischar(varargin{2})  MAS:  Now that we're entering dates as
0060 %        strings rather than MJDs, this doesn't work.
0061         if strcmpi(varargin{2},'interactive') || strcmpi(varargin{2},'noninteractive')
0062             % currentTime,MODE
0063             disp('updateAlphaDatabase:: Scenario 1: append new entries to the database')
0064             SCENARIO = 1;
0065             MODE = varargin{2};
0066             currentTime = varargin{1};
0067             mjdstop = currentTime;
0068             tscans = 1; % hours
0069         elseif ischar(varargin{2}) || varargin{2} > 10000
0070             % startTime,endTime
0071             disp('updateAlphaDatabase:: Scenario 2: replace existing entries in the database')
0072             SCENARIO = 2;
0073             MODE = 'INTERACTIVE';
0074             mjdstart = varargin{1};
0075             mjdstop  = varargin{2};
0076             tscans = 1; % hours
0077         else
0078             % currentTime,scanLength
0079             % varargin{2} is the number of hours to scan in at a time:
0080             disp('updateAlphaDatabase:: Scenario 1: append new entries to the database')
0081             SCENARIO = 1;
0082             MODE = 'INTERACTIVE';
0083             currentTime = varargin{1};
0084             tscans = varargin{2}; % hours
0085             mjdstop  = currentTime;
0086         end
0087     case 3
0088         % Possible options are:
0089         % startTime,endTime,scanlength OR
0090         % startTime,endTime,MODE OR
0091         % currentTime,scanLength,MODE
0092         
0093         if ischar(varargin{3})
0094             % Either:
0095             % startTime,endTime,MODE OR
0096             % currentTime,scanLength,MODE
0097             MODE = varargin{3};
0098             if ischar(varargin{2}) || varargin{2} > 10000
0099                 % startTime,endTime,MODE
0100                 % varargin{2} is the end time
0101                 disp('updateAlphaDatabase:: Scenario 2: replace existing entries in the database')
0102                 SCENARIO = 2;
0103                 mjdstart = varargin{1};
0104                 mjdstop  = varargin{2};
0105                 tscans = 1; % hours
0106             else
0107                 % currentTime,scanLength,MODE
0108                 % varargin{2} is the number of hours to scan in at a time:
0109                 disp('updateAlphaDatabase:: Scenario 1: append new entries to the database')
0110                 SCENARIO = 1;
0111                 currentTime = varargin{1};
0112                 tscans = varargin{2}; % hours
0113                 mjdstop  = currentTime;
0114             end
0115             
0116         else
0117             % startTime,endTime,scanlength
0118             disp('updateAlphaDatabase:: Scenario 2: replace existing entries in the database')
0119             SCENARIO = 2;
0120             MODE = 'INTERACTIVE';
0121             mjdstart = varargin{1};
0122             mjdstop  = varargin{2};
0123             tscans = varargin{3};
0124         end
0125         
0126     case 4
0127         % startTime,endTime,scanlength,MODE
0128         disp('updateAlphaDatabase:: Scenario 2: replace existing entries in the database')
0129         SCENARIO = 2;
0130         MODE = varargin{4};
0131         mjdstart = varargin{1};
0132         mjdstop  = varargin{2};
0133         tscans = varargin{3};
0134         
0135     otherwise
0136         error('updateAlphaDatabase:WrongNumberOfArgs',...
0137         'Incorrect number of arguments!')
0138 end
0139 
0140 disp(['updateAlphaDatabase:: Your selected mode is: ' MODE])
0141 
0142 if ~(strcmp(MODE,'INTERACTIVE') || strcmp(MODE,'NONINTERACTIVE'))
0143     error('updateAlphaDatabase:InvalidMode',...
0144         'Invalid mode specified!')
0145 end
0146 
0147 if strcmp(MODE,'INTERACTIVE')
0148     disp( 'updateAlphaDatabase:: ---------------------------------------------------------')
0149     input( '| Ready to update the CVS repository. Press enter to continue...');
0150     disp( 'updateAlphaDatabase:: | Now updating the alpha database. Enter your CSV password now if asked:')
0151     cmd = 'cvs update -W constants/alpha*.txt';
0152     disp(cmd);
0153     system(cmd);
0154 else
0155     disp('updateAlphaDatabase:: Assuming that you have updated the database already.')
0156 end
0157 
0158 % Read in the database:
0159 disp('updateAlphaDatabase:: Reading in the database...')
0160 if SCENARIO == 1
0161     % append to the database
0162     xold = readAlphaDatabase(mjdstop);
0163     % if we've started a new month, we will get a error message here. So
0164     % read in the previous month's data:
0165     if isempty(xold{1})
0166         if ischar(mjdstop)
0167             stop = mjd2date_v2(tstr2mjd(mjdstop));
0168             xold = readAlphaDatabase(tstr2mjd(mjdstop)-(stop.day+1));
0169         else
0170             stop = mjd2date_v2(mjdstop);
0171             xold = readAlphaDatabase(mjdstop-(stop.day+1));
0172         end
0173     end
0174 elseif SCENARIO == 2
0175     % replace the database files
0176     xold = readAlphaDatabase(mjdstart,mjdstop);
0177 end
0178 
0179 t = xold{1};
0180 
0181 if SCENARIO == 1
0182     mjdstart = t(end);
0183 end
0184 
0185 if ischar(mjdstart)
0186     mjdstart = tstr2mjd(mjdstart);
0187 end
0188 if ischar(mjdstop)
0189     mjdstop = tstr2mjd(mjdstop);
0190 end
0191 if exist('currentTime','var')
0192     if ischar(currentTime)
0193         currentTime = tstr2mjd(currentTime);
0194     end
0195 else
0196     currentTime = 0;
0197 end
0198 
0199 % Is the database up to date?
0200 if SCENARIO == 1 && currentTime <= t(end)
0201     fprintf('Nothing to do: database is up to date!\n');
0202     fprintf('Current time:          %s\n',mjd2string(currentTime));
0203     fprintf('Latest database entry: %s\n',mjd2string(t(end)));
0204 else
0205     % Read in all the data between the latest database time and the present:
0206     disp(['updateAlphaDatabase::   Reading in the data between ' mjd2string(mjdstart) ' and ' mjd2string(mjdstop)])
0207 
0208     % Take the data in tscans-long hour segments:
0209     Nscans = floor((mjdstop-mjdstart)*24/tscans);
0210     
0211     fprintf('***** Number of %d hour segments to scan: %d\n',tscans,Nscans)
0212     % A possibly necessary offset, in case we start or end in a noise diode
0213     % event:
0214     offset = 0;
0215     % the new parameters we calculate will be stored here:
0216     xnew = cell(15,1);
0217     Iautoflag = [];
0218     
0219     % For automatic flagging:
0220     sL = [];
0221     sS = [];
0222     iDS = [];
0223 
0224     if Nscans == 0
0225         Nscans = 1;
0226     end
0227     for k=1:Nscans
0228         stime = mjdstart + (k-1)*tscans/24 + offset;
0229         offset = 0;
0230         if k == Nscans
0231             ftime = mjdstop;
0232         else
0233             ftime = mjdstart + k*tscans/24;
0234         end
0235         fprintf('***** Reading scan %d of %d:\n',k,Nscans)
0236         disp(['updateAlphaDatabase:: *****      Start:  ' mjd2string(stime)])
0237         disp(['updateAlphaDatabase:: *****      Finish: ' mjd2string(ftime)])
0238 
0239         % Read in the data:
0240     try
0241       d0 = pipe_read(mjd2string(stime),mjd2string(ftime));   
0242             % If we end in a ND event, the next scan will start at the
0243             % beginning of the interrupted noise event (with an extra 1
0244             % second buffer before that).
0245             if d0.index.noise_event.fast(end) == 1
0246                 disp('updateAlphaDatabase:: ***** Ended during a noise diode event, adjusting start time of next scan')
0247                 offset = ...
0248                     (find(d0.index.noise_event.fast == 0, 1, 'last' ) - ...
0249                     length(d0.index.noise_event.fast) - 100) / (24 * 3600 * 100);
0250             end
0251 
0252             % Calculate the database parameters for each mode:
0253             d = d0;
0254             d = assembleAlphaStreams(d,'CLASSIC');
0255             [tDD,ADD,GDD,varADD,varGDD,Td,rd] = calculateAlpha(d,'CLASSIC');
0256             
0257             if mjdstart > tstr2mjd('01-oct-2011:00:00:00')
0258                 d = d0;
0259                 d = assembleAlphaStreams(d,'FILTERED');
0260                 [tF,AF,GF,varAF,varGF] = calculateAlpha(d,'FILTERED');
0261             else
0262                 % filtered data reduction is not valid, so make them all
0263                 % zeros.
0264                 tF = zeros(size(tDD));
0265                 AF = zeros(size(ADD,1),3);
0266                 GF = zeros(size(ADD,1),5);
0267                 varAF = zeros(size(ADD,1),3);
0268                 varGF = zeros(size(ADD,1),5);
0269             end
0270             
0271             d = d0;
0272             d = assembleAlphaStreams(d,'POLONLY');
0273             [tPO,APO,GPO,varAPO,varGPO,T0,r0,hrz,equ,offStPos,onEndPos,offEndPos,onStPos] =...
0274                 calculateAlpha(d,'POLONLY');
0275             
0276             d = d0;
0277             % Interpolate the temperature data
0278             if size(d.antenna0.thermal.dlpTemperatureSensors,1) ~= size(d.antenna0.receiver.utc,1)
0279                 d.antenna0.thermal.dlpTemperatureSensors = interp1(d.array.frame.utc,d.antenna0.thermal.dlpTemperatureSensors,d.antenna0.receiver.utc,'linear');
0280             end
0281             
0282             if ~isempty(tDD)
0283                 
0284                 % For flagging the alpha values we first want to apply the
0285                 % alpha corrections and see whether we actually separate sky
0286                 % and load.
0287                 thisAlpha = [tDD ADD GDD Td rd hrz equ];
0288 
0289                 d.correction.alpha.indices = [offStPos onEndPos offEndPos onStPos];
0290                 d.correction.alpha.values = thisAlpha;
0291                 d = assembleAlphaStreams(d,'CLASSIC');
0292                 d = applyAlpha(d,tDD,ADD,GDD,'CLASSIC');
0293                 
0294                 % These parameters determined through experiment:
0295                 maxInterceptDiff = 0.7; % the maximum allowed intercept difference
0296                 minInterceptDiff = 0.3; % the minimum allowed intercept difference
0297                 maxSlope = 0.004; % the maximum allowed slope
0298                 flagParams = [maxInterceptDiff minInterceptDiff maxSlope];
0299                 % run the flagging routine
0300                 [flags,slopeLoad,slopeSky,interceptDiffSky] = flagAlphaOnSlope(d, flagParams);
0301                                 
0302                 tTemp = [];
0303                 for m=1:length(onStPos)
0304                     tTemp = [tTemp; mean(d.antenna0.thermal.dlpTemperatureSensors(onStPos(m):onEndPos(m),1:7),1)];
0305                 end
0306                 
0307                 xnew{1} = cat(1,xnew{1},tDD);
0308                 xnew{2} = cat(1,xnew{2},ADD);
0309                 xnew{3} = cat(1,xnew{3},GDD);
0310                 xnew{4} = cat(1,xnew{4},varADD);
0311                 xnew{5} = cat(1,xnew{5},varGDD);
0312                 xnew{6} = cat(1,xnew{6},GPO(:,[1,5]));
0313                 xnew{7} = cat(1,xnew{7},varGPO(:,[1,5]));
0314                 xnew{8} = cat(1,xnew{8},GF);
0315                 xnew{9} = cat(1,xnew{9},AF);
0316                 xnew{10} = cat(1,xnew{10},varGF);
0317                 xnew{11} = cat(1,xnew{11},varAF);
0318                 xnew{12} = cat(1,xnew{12},rd);
0319                 xnew{13} = cat(1,xnew{13},Td);
0320                 xnew{14} = cat(1,xnew{14},[hrz*pi/180 equ]);
0321                 xnew{15} = cat(1,xnew{15},tTemp);
0322                 
0323                 sL = [sL; slopeLoad];
0324                 sS = [sS; slopeSky];
0325                 iDS = [iDS; interceptDiffSky];
0326                 
0327                 Iautoflag = [Iautoflag; flags(:)];
0328                 
0329             end
0330         catch
0331             disp('updateAlphaDatabase:: ***** BAD DATA; SKIPPING TO NEXT SET')
0332         end
0333     end
0334     if ~isempty(xnew{1})
0335         Iautoflag = any(sL>0.006,2) | any(any(sS>0.006,3),2) ...
0336             | any(iDS<0.2,2) | any(iDS>0.8,2);
0337         
0338 %         tnI = tnew;
0339 %         AnI = Anew;
0340 %         GnI = Gnew;
0341 %         % SJCM -- COMMENTED OUT THIS FIRST LINE BECAUSE IT CAUSES
0342 %         % FILTERALPHA TO BREAK IF THE FIRST T VALUE IS NAN'ED
0343 %         tnI(Iautoflag) = [];
0344 %         AnI(Iautoflag,:) = [];
0345 %         GnI(Iautoflag,:) = [];
0346 %
0347 %         Iflag = zeros(size(Iautoflag));
0348 %         thisFlag = logical(filterAlpha(tnI,AnI,GnI,0.15));
0349 %         Iflag(~Iautoflag) = thisFlag;
0350         
0351         % Combine all the automatic flags
0352         autoflag = 0;
0353         if autoflag
0354             IflaggedNAN = any(isnan(xnew{2}),2) | Iautoflag | Iflag;
0355         else
0356             IflaggedNAN = any(isnan(xnew{2}),2);
0357         end
0358 
0359         if strcmp(MODE,'INTERACTIVE')
0360             % Give the data a manual check:
0361 %             IflaggedA1 = flagloop(tnew,Anew(:,1),IflaggedNAN,'alpha(1)');
0362 %             IflaggedG5 = flagloop(tnew,Gnew(:,5),IflaggedNAN|IflaggedA1,'gain(5)');
0363 
0364             % Create the cells to loop over with the multiple-flag version.
0365             % The alpha_DD parameters:
0366             data_A = {xnew{2}(:,1) xnew{2}(:,4) xnew{2}(:,5)};
0367             varData_A = {xnew{4}(:,1) xnew{4}(:,4) xnew{4}(:,5)};
0368             ident_A = {'alpha DD I1' 'alpha DD P' 'alpha DD I2'};
0369             % The gain_DD parameters:
0370             data_G = {xnew{3}(:,1) xnew{3}(:,4) xnew{3}(:,5)};
0371             varData_G = {xnew{5}(:,1) xnew{5}(:,4) xnew{5}(:,5)};
0372             ident_G = {'1/gain DD I1' '1/gain DD P' '1/gain DD I2'};
0373             
0374             IflaggedA = flagloop(xnew{1},data_A,varData_A,IflaggedNAN,ident_A);
0375             IflaggedG = flagloop(xnew{1},data_G,varData_G,IflaggedNAN|IflaggedA,ident_G);
0376             Iflagged = IflaggedA | IflaggedG | IflaggedNAN;
0377         else
0378             % Proceed with the automatic flags
0379             Iflagged = IflaggedNAN;
0380         end
0381 
0382         if strcmp(MODE,'INTERACTIVE')
0383             proceed = input('Write to disk? Press 1 (then enter) if yes: ');
0384             TYPE = zeros(size(xnew{1}));
0385         else
0386             TYPE = ones(size(xnew{1}));
0387             proceed = 1;
0388         end
0389         
0390         % Update the database:
0391         if proceed == 1
0392             xnew{16} = TYPE;
0393             xnew{17} = Iflagged;
0394             if SCENARIO == 1
0395                 % first remove any duplicate entries before writing to
0396                 % disk.
0397                 [ts,I,J] = unique(xnew{1});
0398                 for k=1:length(xnew)
0399                     xnew{k} = xnew{k}(I,:);
0400                 end
0401                 
0402                 
0403                 disp('updateAlphaDatabase:: ...Updating database with new entries')
0404                 fnames = writeAlphaDatabase(xnew,'APPEND');
0405                 disp('updateAlphaDatabase:: Done!')
0406                 
0407                 % If we've added new files to the database, we need to add
0408                 % them to the CVS repository now.
0409                 old = mjd2date_v2(xold{1}(end));
0410                 new = mjd2date_v2(xnew{1}(end));
0411                 if new.month ~= old.month
0412                     % add the file name for every month in fnames which is
0413                     % after the time in old:
0414                     if new.month > old.month && new.year == old.year
0415                         ms = (old.month+1):new.month;
0416                         ys = ones(size(ms))*new.year;
0417                     elseif new.month < old.month && new.year > old.year
0418                         ms = (old.month+1):(new.month+12);
0419                         ys = ones(size(ms))*old.year;
0420                         ys(ms>12) = new.year;
0421                         ms(ms>12) = ms(ms>12)-12;
0422                     else
0423                         error('updateAlphaDatabase:CVSUpdateError',...
0424                             'Trying to append to database with old data!')
0425                     end
0426                     
0427                     cmd = 'cvs add';
0428                     for k=1:length(ms)
0429                         cmd = [cmd ' ' sprintf('constants/alphaDatabase_%04d_%02d.txt',ys(k),ms(k))]; 
0430                     end
0431                     input('Ready to add newly-created database file[s] to CVS? Press enter to continue.');
0432                     disp('updateAlphaDatabase:: Enter CVS password now if asked:')
0433                     disp(cmd);
0434                     system(cmd);
0435                 end
0436                     
0437             elseif SCENARIO == 2
0438                 disp('updateAlphaDatabase:: ...Replacing database entries with new values')
0439                 % Remove any existing entries from the database:
0440                 tn = xnew{1};
0441                 Iremove = xold{1} >= (tn(1)-5/3600/24) & xold{1} <= (tn(end)+5/3600/24);
0442                 
0443                 % unsorted data:
0444                 xuns = xold;
0445                 for k=1:length(xold)
0446                     xuns{k} = [xold{k}(~Iremove,:); xnew{k}];
0447                 end
0448 
0449                 % Sort the data and remove duplicate entries
0450                 [ts,Isorted] = sort(xuns{1});
0451                 for k=1:length(xuns)
0452                     xuns{k} = xuns{k}(Isorted,:);
0453                 end
0454                 % xuns is now sorted
0455                 [ts,I,J] = unique(ts);
0456                 for k=1:length(xuns)
0457                     xuns{k} = xuns{k}(I,:);
0458                 end
0459                 % xuns is now sorted and contains only unique entries
0460 
0461                 fnames = writeAlphaDatabase(xuns,'REPLACE');
0462                 disp('updateAlphaDatabase:: Done!')
0463             end
0464    
0465             if strcmp(MODE,'INTERACTIVE')
0466                 % Remind the user to run cvs commit:
0467                 uname = getenv('USER');
0468 
0469                 disp( 'updateAlphaDatabase:: ---------------------------------------------------------')
0470                 disp( 'updateAlphaDatabase:: | Now commit the updated alphaDatabase to the repository.')
0471                 input('| Press return when ready...');
0472                 disp('updateAlphaDatabase:: | Enter CVS password now if asked:')
0473                 
0474                 % create the commit command
0475                 cmd = sprintf('cvs commit -m''%s: alpha database updated''',uname);
0476                 for k = 1:length(fnames)
0477                     fn = fnames{k};
0478                     cmd = [cmd ' ' fn(strfind(fn,'constants'):end)];
0479                 end
0480                 disp(cmd)
0481                 system(cmd);
0482             else
0483                 disp('updateAlphaDatabase:: Assuming that your script will update the cvs repository.')
0484             end
0485         end
0486     else
0487         fprintf('No noise diode events in this time period!\n');
0488     end
0489 end
0490 
0491 end
0492 
0493 function Iflagged = flagloop(t,data,varData,Ipreflagged,title)
0494 % Check the quality of the data:
0495 disp(['updateAlphaDatabase:: Flagging ' title])
0496 redo = 1;
0497 while redo
0498     Iflagged = flagAlphaGraphically_multiple(t,data,varData,...
0499         Ipreflagged,title,1);
0500     rsel = input('updateAlphaDatabase::flagloop:: Redo? Press enter to continue to next, any other key to redo:','s');
0501     close(gcf)
0502     if isempty(rsel)
0503         redo = 0;
0504     else
0505         redo = 1;
0506     end
0507 end
0508 end
0509 
0510 function [flags,slopeLoad,slopeSky,interceptDiffSky] = flagAlphaOnSlope(d, flagParams)
0511 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0512 %
0513 %  function flags = flagAlphaOnSlope(d, flagParams)
0514 %
0515 %  function to automatically flag results in the alpha correction
0516 %
0517 %  sjcm
0518 %
0519 % ogk - updated 12 March 2012
0520 %
0521 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0522 
0523 % what we want to flag are the noise diode events where we don't see any
0524 % improvement.
0525 
0526 %%%%%%%%%%%%%%%%%%%%%%%%%
0527 
0528 if size(d.correction.alpha.indices,2) == 4
0529     % Before 8-Sep-2010
0530     si = d.correction.alpha.indices(:,1); % when the first off noise diode event starts
0531     ei = d.correction.alpha.indices(:,2); % when the on noise diode event ends
0532     esi = d.correction.alpha.indices(:,4); % when the on noise diode event starts
0533 else
0534     % After 8-Sep-2010
0535     si = d.correction.alpha.indices(:,1); % when the first off noise diode event starts
0536     ei = d.correction.alpha.indices(:,3); % when the on noise diode event ends
0537     esi = d.correction.alpha.indices(:,6); % when the on noise diode event starts
0538 end
0539 %%%%%%%%%%%%%%%%%%%%%%%%%
0540 numEvents = size(d.correction.alpha.indices,1);
0541 
0542 maxInterceptDiff = flagParams(1);
0543 minInterceptDiff = flagParams(2);
0544 maxSlope = flagParams(3);
0545 
0546 % split up all the noise diode events
0547 for m=1:numEvents
0548   x{1,m} = d.antenna0.receiver.data(si(m):ei(m),[1,2]);
0549   x{2,m} = d.antenna0.receiver.data(si(m):ei(m),[9,10]);
0550 end
0551 
0552 % index 2 in x{1:2,m} should be flat
0553 % index 1 in x{1:2,m} should have a step in them
0554 
0555 % first we fit a line to the flat ones, and make sure there's no slope
0556 % greater than that specified in flagParams(2);
0557 
0558 for m=1:numEvents
0559   thisX = [1:1:size(x{1,m},1)]';
0560   thisY = x{1,m}(:,2);
0561   thisX = thisX/100;
0562 
0563   % Channel 1 slope and intercept
0564   [slopeVal(m,1), interceptVal(m,1)] = linfit(thisX, thisY);
0565   
0566   % Channel 2 slope and intercept
0567   thisY = x{2,m}(:,2);
0568   [slopeVal(m,2), interceptVal(m,2)] = linfit(thisX, thisY);  
0569 end
0570 % find those slopeVals that are greater than the parameters
0571 slopeVal = abs(slopeVal);
0572 slopeLoad = slopeVal;
0573 badSlope1 = slopeVal>maxSlope;
0574 badSlope1 = sum(badSlope1,2)>0; % flag if either of the channels has a bad slope val
0575 
0576 % next we fit a line to both sides of the step ones, and
0577 % A) make sure there's no slope greater than that specified in flagParams(2)
0578 % B) the y-intercept of the horizontal lines are no more than flagParams(1)
0579 % apart
0580 clear slopeVal;
0581 clear interceptVal;
0582 for m=1:numEvents
0583   for nchan=1:2
0584     numSamples = size(x{nchan,m},1);
0585     midSample  = esi(m)-si(m); % the official start of the on noise diode event
0586     slack = 50; % don't trust it
0587     
0588     thisX1 = [1:(midSample-slack)]';
0589     thisY1 = x{nchan,m}(thisX1,1);
0590     thisX1 = thisX1/100;
0591     
0592     thisX2 = [midSample+slack:numSamples]';
0593     thisY2 = x{nchan,m}(thisX2,1);
0594     thisX2 = thisX2/100;
0595 
0596     [slopeA intA] = linfit(thisX1, thisY1);
0597     [slopeB intB] = linfit(thisX2, thisY2);
0598     
0599     thisSlope = [slopeA slopeB];
0600     thisInt   = [intA intB];
0601     
0602     slopeVal(m,nchan,:) = thisSlope;
0603     interceptVal(m,nchan,:) = thisInt;
0604   end
0605 end
0606 
0607 % first we check for criteria A
0608 slopeVal = abs(slopeVal);
0609 slopeSky = slopeVal;
0610 badSlope2 = slopeVal>maxSlope;
0611 badSlope2 = sum(sum(badSlope2,3),2)>0;
0612 
0613 % next we check for criteria B
0614 interceptDiff = interceptVal(:,:,1) - interceptVal(:,:,2);
0615 interceptDiff = abs(interceptDiff);
0616 interceptDiffSky = interceptDiff;
0617 badIntercept  = sum(sum(interceptDiff>maxInterceptDiff | ...
0618     interceptDiff<minInterceptDiff,3),2)>0;
0619 
0620 % now we combine it all to find out which events should be flagged
0621 badEvents = badSlope1 | badSlope2 | badIntercept;
0622 
0623 flags = badEvents;
0624 
0625 end